lcsim/src/org/lcsim/contrib/SteveMagill
diff -N TSCDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TSCDriver.java 29 May 2008 21:33:37 -0000 1.1
@@ -0,0 +1,685 @@
+package org.lcsim.contrib.compile.SteveMagill;
+
+// This driver extrapolates tracks into CAL, associating clusters to the tracks, testing E/p
+// output : modified hit map and track cal clusters
+
+import java.util.*;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.util.swim.*;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.aida.*;
+import org.lcsim.geometry.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.recon.ztracking.cheater.*;
+import org.lcsim.mc.fast.tracking.*;
+import hep.aida.*;
+import org.lcsim.spacegeom.*;
+import org.lcsim.geometry.util.CalorimeterIDDecoder;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
+
+public class TSCDriver extends Driver
+{
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IAnalysisFactory analysisFactory = IAnalysisFactory.create();
+ private AIDA aida = AIDA.defaultInstance();
+ private boolean _initialized = false;
+ // The number of layers in the EM Calorimeter Barrel, Endcap
+ private int _numbemlayers;
+ private int _numecemlayers;
+// The number of layers in the Hadronic Calorimeter Barrel, Endcap
+ private int _numbhadlayers;
+ private int _numechadlayers;
+// The EM calorimeter hits, by layer in barrel, endcap
+ private List[] _emBLayerHits;
+ private List[] _emECLayerHits;
+// The HAD calorimeter hits, by layer in barrel, endcap
+ private List[] _hadBLayerHits;
+ private List[] _hadECLayerHits;
+// The radii of the barrel calorimeter layers
+ private double[] BRadii = new double[100];
+ private double[] ECZs = new double[100];
+ // Z extent of the central barrel calorimeters.
+ private double _embZ; //EM Barrel Z
+ private double[] _emecZ; //EM Endcap Z
+ private double _hadbZ; //HAD barrel z
+ private double[] _hadecZ; // HAD endcap z
+ private double _dminE;
+ private double _dminH;
+ private double _dminTrC;
+ private double _mindist;
+ private int _nloops;
+ private double _mineop;
+ private double _maxeop;
+ private String[] _clusternames;
+ private String _oclname;
+ private boolean trshdb = false;
+ private double zField;
+// private String[] _hitmapnames;
+// CoreReclusterer cr;
+
+ public TSCDriver(double mindist, int nloops, double mineop, double maxeop)
+ {
+ // add arguments if needed
+ _mindist = mindist; // minimum distance between objects to merge (track, clusters) (0.015)
+ _nloops = nloops; // number of iterations to merge clusters (until min E/p is met) (4)
+ _mineop = mineop; // minimum E/p value (0.65)
+ _maxeop = maxeop; // maximum E/p value (1.25 is 1+sigam/E for 60%/sqrt(E))
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event); // executes all added drivers
+
+// Initialize things here
+ if(!_initialized)
+ {
+ // setup specific detector stuff here
+ Detector det = event.getDetector();
+ Hep3Vector ip = new BasicHep3Vector();
+ zField = det.getFieldMap().getField(ip).z();
+// System.out.println("B Field " + zField + " Tesla");
+
+ // Setup layer counts and distances for barrel, endcap
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ _embZ = emb.getZMin();
+ _numbemlayers = emb.getLayering().getLayerCount();
+// System.out.println("EM Barrel Layers " +_numbemlayers);
+ for (int i=0; i<_numbemlayers; ++i)
+ {
+ BRadii[i]=emb.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("EM Barrel Layer Number " +i+ " EM Barrel Radius " +_emBRadii[i]);
+ }
+
+ CylindricalCalorimeter emec = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ _numecemlayers = emec.getLayering().getLayerCount();
+// System.out.println("EM EC Layers " +_numecemlayers);
+ for (int i=0; i<_numecemlayers; ++i)
+ {
+ ECZs[i] = emec.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("EM Endcap Layer Number " +i+ " EM Endcap Z " +_emecZ[i]);
+ }
+
+ CylindricalCalorimeter hadb = ((CylindricalCalorimeter) det.getSubdetectors().get("HADBarrel"));
+ _hadbZ = hadb.getZMin();
+ _numbhadlayers = hadb.getLayering().getLayerCount();
+// System.out.println("HAD Barrel Layers " +_numbhadlayers);
+ for (int i=0; i<_numbhadlayers; ++i)
+ {
+ BRadii[i+_numbemlayers]=hadb.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("HAD Barrel Layer Number " +i+ " HAD Barrel Radius " +_hadBRadii[i]);
+ }
+
+ CylindricalCalorimeter hadec = ((CylindricalCalorimeter) det.getSubdetectors().get("HADEndcap"));
+ _numechadlayers = hadec.getLayering().getLayerCount();
+// System.out.println("HAD Endcap Layers " +_numechadlayers);
+ for (int i=0; i<_numechadlayers; ++i)
+ {
+ ECZs[i+_numecemlayers] = hadec.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("HAD Endcap Layer Number " +i+ " HAD Endcap Z " +_hadecZ[i]);
+ }
+
+ _initialized = true;
+ } // end of initialization section
+
+ // combine all clusters into a single shower cluster list
+ List<BasicCluster> showclus = new ArrayList<BasicCluster>();
+// System.out.println(" Num of Clus Lists " + _clusternames.length);
+ for (int i=0; i<_clusternames.length; i++)
+ {
+ try
+ {
+ List<BasicCluster> clus = event.get(BasicCluster.class,_clusternames[i]);
+// System.out.println(" Num of Clus in list " + clus.size());
+ for (BasicCluster cl : clus)
+ {
+ BasicCluster shclus = new BasicCluster();
+ if (cl.getSize()>0)
+ {
+ if (trshdb) aida.cloud1D("Num Hits in Show Clus TSC").fill(cl.getSize());
+ aida.cloud2D("Num hits Show Clus vs E TSC").fill(cl.getEnergy(),cl.getSize());
+ shclus.addCluster(cl);
+ showclus.add(shclus);
+ }
+ }
+ }
+ catch(java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("requested object not found in event " + _clusternames[i]);
+ }
+ }
+// System.out.println("Number Shower Clusters before Matching " + showclus.size());
+ if (trshdb) aida.cloud1D("Num Shower Clus before TSCMatch").fill(showclus.size());
+
+
+ // create lists of track-cal associated clusters for display
+ List<BasicCluster> trkcalclus = new ArrayList<BasicCluster>();
+ List<BasicCluster> trkshoclus = new ArrayList<BasicCluster>();
+ List<BasicCluster> trkcorclus = new ArrayList<BasicCluster>();
+
+ // Make a map to link track and final cluster association
+ Map<Track, BasicCluster> trkclusmap = new HashMap<Track, BasicCluster>();
+
+ // Get maps linking IL and mip clusters to tracks, tracks linked to SpacePoints
+ Map<Track, BasicCluster> trmipmap = (Map<Track, BasicCluster>) event.get("TrackMipMap");
+ Map<BasicCluster, Integer> clusilmap = (Map<BasicCluster, Integer>) event.get("MipClusILMap");
+ Map<Track, SpacePoint> tre0map = (Map<Track, SpacePoint>) event.get("TrackXE0Map");
+ Map<Track, SpacePoint> tresmmap = (Map<Track, SpacePoint>) event.get("TrackXEShMaxMap");
+ Map<Track, SpacePoint> trh0map = (Map<Track, SpacePoint>) event.get("TrackXH0Map");
+
+ double TotTrP = 0;
+ // use track-mip map to start, get theta and phi at IL and mip cluster center
+ for (Track itr : trmipmap.keySet())
+ {
+ double TrClE = 0;
+ int ndist = 0;
+ BasicCluster trclus = new BasicCluster();
+ BasicCluster trshclus = new BasicCluster();
+ BasicCluster trcorclus = new BasicCluster();
+ int niter = 0;
+ // get track and associated mip cluster
+ BasicCluster trmclus = trmipmap.get(itr);
+ trclus.addCluster(trmclus); // add mip cluster hits to track cluster
+ TrClE = trmclus.getEnergy();
+ int ClIL = clusilmap.get(trmclus);
+// System.out.println("InteractionLayer for track " + ClIL);
+ double[] trmpos = trmclus.getPosition();
+ Hep3Vector mclpos = new BasicHep3Vector(trmpos);
+ double trpt = Math.sqrt(itr.getPX()*itr.getPX()+itr.getPY()*itr.getPY());
+ double[] trp = itr.getMomentum();
+ Hep3Vector trp3 = new BasicHep3Vector(trp);
+ double[] trrp = itr.getReferencePoint();
+ double[] trpar = itr.getTrackParameters();
+ double[] tror = new double[3];
+ tror[0] = -trpar[0]*Math.sin(trpar[1])+trrp[0];
+ tror[1] = trpar[0]*Math.cos(trpar[1])+trrp[1];
+ tror[2] = trpar[3]+trrp[2];
+ Hep3Vector tror3 = new BasicHep3Vector(tror);
+ SpacePoint trorsp = new SpacePoint(tror3);
+ int trq = itr.getCharge();
+// System.out.println("Track Charge " + trq + " Track P " + trp3 + " Origin " + tror3);
+ double TrP = Math.sqrt(itr.getPX()*itr.getPX()+itr.getPY()*itr.getPY()+itr.getPZ()*itr.getPZ());
+ TotTrP += TrP;
+ HelixSwimmer tswim = new HelixSwimmer(zField);
+ // swim track to cluster position, arguments are momentum, origin, and charge
+// tswim.setTrack(trp3, tror3, trq); // deprecated, uses spacepoint for orogin
+ tswim.setTrack(trp3, trorsp, trq);
+ double tobrad = tswim.getDistanceToRadius(BRadii[ClIL]);
+ double toecz = tswim.getDistanceToZ(ECZs[ClIL]);
+ double trtheta = 0;
+ double trphi = 0;
+ double trilth = 0; // this is theta at IL
+ double trilph = 0; // this is phi at IL
+ double clth = 0; // this is theta at mip cluster center
+ double clph = 0; // this is phi at mip cluster center
+// System.out.println("Distance to B IL " + tobrad + " to EC IL " + toecz);
+ if (tobrad<Math.abs(toecz)) // in barrel
+ {
+// System.out.println(" In Barrel " + tobrad);
+ SpacePoint trSP = tswim.getPointAtLength(tobrad);
+
+ trtheta = trSP.theta(); // track theta at IL
+ trphi = trSP.phi(); // track phi at IL
+ trilth = Math.atan(trSP.rxy()/trSP.z());
+ if (trilth<0) trilth+=Math.PI;
+ trilph = Math.atan2(trSP.y(),trSP.x());
+ if (trilph<0) trilph+=2*Math.PI;
+
+ if (ClIL==0)
+ {
+ clth = trilth; // no mip track, so IL,cluster same
+ clph = trilph;
+ } else
+ {
+ // get mip cluster position coordinates
+ double clr = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
+ clth = Math.atan(clr/trmpos[2]);
+ if (clth<0) clth+=Math.PI;
+ clph = Math.atan2(trmpos[1],trmpos[0]);
+ if (clph<0) clph+=2*Math.PI;
+ }
+
+ if (trshdb)
+ {
+ aida.cloud1D("Theta of Track at IL B").fill(trilth); // this is correct value
+ aida.cloud1D("Theta of mip cluster c in B").fill(clth);
+ aida.cloud1D("Phi of Track at IL B").fill(trilph);
+ aida.cloud1D("Phi of mip cluster c in B").fill(clph);
+ aida.cloud2D("Theta vs Phi of Track at IL B").fill(trilph,trilth);
+ aida.cloud2D("Theta vs Phi of mip cluster c in B").fill(clph,clth);
+ }
+ } else // in endcap
+ {
+// System.out.println(" In Endcap " + toecz);
+ SpacePoint trSP = tswim.getPointAtLength(toecz);
+ trtheta = trSP.theta();
+ trphi = trSP.phi();
+ trilth = Math.atan(trSP.rxy()/trSP.z());
+ if (trilth<0) trilth+=Math.PI;
+ trilph = Math.atan2(trSP.y(),trSP.x());
+ if (trilph<0) trilph+=2*Math.PI;
+
+ if (ClIL==0)
+ {
+ clth = trilth;
+ clph = trilph;
+ } else
+ {
+ double clr = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
+ clth = Math.atan(clr/trmpos[2]);
+ if (clth<0) clth+=Math.PI;
+ clph = Math.atan2(trmpos[1],trmpos[0]);
+ if (clph<0) clph+=2*Math.PI;
+ }
+ if (trshdb)
+ {
+ aida.cloud1D("Theta of mip cluster c in EC").fill(clth);
+ aida.cloud1D("Theta of Track at IL EC").fill(trilph);
+ aida.cloud1D("Phi of Track at IL EC").fill(trilph);
+ aida.cloud1D("Phi of mip cluster c in EC").fill(clph);
+ aida.cloud2D("Theta vs Phi of Track at IL EC").fill(trilph,trilth);
+ aida.cloud2D("Theta vs Phi of mip cluster c in EC").fill(clph,clth);
+ }
+ }
+
+ // find clusters very close to track - core clusters
+ int ncore = 0;
+ double[] trcclth = new double[100];
+ double[] trcclph = new double[100];
+ for (Iterator<BasicCluster> shcl = showclus.iterator(); shcl.hasNext();)
+ {
+ BasicCluster ibcl = shcl.next();
+ // don't consider if too large in E/p
+ // System.out.println(" Ratio for this cluster is " + (TrClE+ibcl.getEnergy())/TrP);
+ if (ibcl.getEnergy()/TrP > _maxeop) continue;
+ // System.out.println("Test this cluster");
+ double[] ccpos = ibcl.getPosition();
+ double ccph = Math.atan2(ccpos[1],ccpos[0]);
+ if (ccph<0) ccph+=2*Math.PI;
+ double ccr = Math.sqrt(ccpos[0]*ccpos[0]+ccpos[1]*ccpos[1]);
+ double ccth = Math.atan(ccr/ccpos[2]);
+ if (ccth<0) ccth+=Math.PI;
+
+ // have cluster theta and phi, now calculate distances for test
+ // will use 2 distances, depending on where IL is
+ double dist1 = 999;
+ double dist2 = 999;
+ if (ClIL==0)
+ {
+ // use maps at E0 and H0
+ SpacePoint trSPe0 = tre0map.get(itr);
+ double tre0th = Math.atan(trSPe0.rxy()/trSPe0.z());
+ if (tre0th<0) tre0th+=Math.PI;
+ double tre0ph = Math.atan2(trSPe0.y(),trSPe0.x());
+ if (tre0ph<0) tre0ph+=2*Math.PI;
+ double d1ccth = Math.abs(tre0th-ccth);
+ double d1ccph = Math.abs(tre0ph-ccph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-ccth);
+ double d2ccph = Math.abs(trh0ph-ccph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ } else if (ClIL>0 && ClIL<30)
+ {
+ // use maps at IL and H0
+ double d1ccth = Math.abs(trilth-ccth);
+ double d1ccph = Math.abs(trilph-ccph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-ccth);
+ double d2ccph = Math.abs(trh0ph-ccph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ } else if (ClIL>29)
+ {
+ // use maps at IL and mip cluster center
+ double d1ccth = Math.abs(trilth-ccth);
+ double d1ccph = Math.abs(trilph-ccph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ double d2ccth = Math.abs(clth-ccth);
+ double d2ccph = Math.abs(clph-ccph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ }
+ // take smallest distance to test
+ if (dist1<dist2) dist2 = dist1;
+ if (dist2<_mindist)
+ {
+ TrClE += ibcl.getEnergy();
+ trcclth[ncore] = ccth;
+ trcclph[ncore] = ccph;
+ ncore++;
+ // add cluster to cands, remove from cluster list
+ trcorclus.addCluster(ibcl);
+ trclus.addCluster(ibcl);
+ shcl.remove();
+ }
+ }
+ if (trshdb) aida.cloud1D("Number of Core Clusters added to track").fill(ncore);
+ if (trshdb) aida.cloud1D("E over P after adding Core Clusters").fill(trclus.getEnergy()/TrP);
+
+ // check all clusters for matches using spacepoint of extrapolated track and cluster coords
+ // for each cluster, loop over all hits in the cluster, checking if any hits are close to track
+ // if so, remove entire cluster from list according to E/p test. If E/p jumps from below min
+ // to greater than 1.5, break up cluster with NN and retry.
+ do {
+ niter++;
+ double showclE = 0;
+ for (Iterator<BasicCluster> shcl = showclus.iterator(); shcl.hasNext();)
+ {
+ BasicCluster ishcl = shcl.next();
+
+ // if this cluster has E/p > _maxeop, then no point in considering it for match
+ if (ishcl.getEnergy()/TrP > _maxeop) continue;
+ // first, calculate distance of hits from track, mip, core
+ List<CalorimeterHit> clhits = ishcl.getCalorimeterHits();
+ int nhmtch = 0;
+ int nhits = 0;
+ for (CalorimeterHit iclhit : clhits)
+ {
+ nhits++;
+ // check distance of hit from track - add up hits to match cluster
+ double[] htpos = iclhit.getPosition();
+ double htph = Math.atan2(htpos[1],htpos[0]);
+ if (htph<0) htph+=2*Math.PI;
+ double htr = Math.sqrt(htpos[0]*htpos[0]+htpos[1]*htpos[1]);
+ double htth = Math.atan(htr/htpos[2]);
+ if (htth<0) htth+=Math.PI;
+
+ // now calculate distances same as for cores
+ double dist1 = 999.;
+ double dist2 = 999.;
+ if (ClIL==0)
+ {
+ // use maps at E0 and H0
+ SpacePoint trSPe0 = tre0map.get(itr);
+ double tre0th = Math.atan(trSPe0.rxy()/trSPe0.z());
+ if (tre0th<0) tre0th+=Math.PI;
+ double tre0ph = Math.atan2(trSPe0.y(),trSPe0.x());
+ if (tre0ph<0) tre0ph+=2*Math.PI;
+ double d1ccth = Math.abs(tre0th-htth);
+ double d1ccph = Math.abs(tre0ph-htph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-htth);
+ double d2ccph = Math.abs(trh0ph-htph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ } else if (ClIL>0 && ClIL<30)
+ {
+ // use maps at IL and H0
+ double d1ccth = Math.abs(trilth-htth);
+ double d1ccph = Math.abs(trilph-htph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-htth);
+ double d2ccph = Math.abs(trh0ph-htph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ } else if (ClIL>29)
+ {
+ // use maps at IL and mip cluster center
+ double d1ccth = Math.abs(trilth-htth);
+ double d1ccph = Math.abs(trilph-htph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ double d2ccth = Math.abs(clth-htth);
+ double d2ccph = Math.abs(clph-htph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ }
+ // if a core cluster was added, also check distance to it
+ double distcore = 999.;
+ if (ncore>0)
+ {
+ double dcccth = Math.abs(trcclth[0]-htth);
+ double dcccph = Math.abs(trcclph[0]-htph);
+ if (dcccph>Math.PI) dcccph = 2*Math.PI-dcccph;
+ distcore = Math.sqrt(dcccth*dcccth+dcccph*dcccph);
+ }
+ // take smallest distance to test
+ double distt = 999.;
+ if (dist1<dist2 && dist1<distcore) distt = dist1;
+ if (dist2<dist1 && dist2<distcore) distt = dist2;
+ if (distcore<dist1 && distcore<dist2) distt = distcore;
+ if (distt<niter*_mindist)
+ {
+ nhmtch++;
+ }
+ }
+ double dnhits = nhits;
+ double dnhmtch = nhmtch;;
+ double drat = dnhmtch/dnhits;
+ if (trshdb)
+ {
+ if (nhmtch>0 && niter==1) aida.cloud1D("Ratio of matches to Cl hits").fill(drat);
+ if (nhmtch>0 && nhits>10 && niter==1) aida.cloud1D("Ratio hit match Clus large Cl").fill(drat);
+ }
+
+ // now do the same for clusters, not hits
+ int nclmtch = 0;
+ double[] clpos = ishcl.getPosition();
+ double cph = Math.atan2(clpos[1],clpos[0]);
+ if (cph<0) cph+=2*Math.PI;
+ double cshr = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]);
+ double cth = Math.atan(cshr/clpos[2]);
+ if (cth<0) cth+=Math.PI;
+ if (trshdb) aida.cloud1D("Calculated shower cluster phi").fill(cph);
+ if (trshdb) aida.cloud1D("Calculated shower cluster theta").fill(cth);
+ if (trshdb) aida.cloud2D("Theta vs Phi shower cluster").fill(cph,cth);
+
+ // now calculate distances same as for cores
+ double dist1 = 999;
+ double dist2 = 999;
+ if (ClIL==0)
+ {
+ // use maps at E0 and H0
+ SpacePoint trSPe0 = tre0map.get(itr);
+ double tre0th = Math.atan(trSPe0.rxy()/trSPe0.z());
+ if (tre0th<0) tre0th+=Math.PI;
+ double tre0ph = Math.atan2(trSPe0.y(),trSPe0.x());
+ if (tre0ph<0) tre0ph+=2*Math.PI;
+ double d1ccth = Math.abs(tre0th-cth);
+ double d1ccph = Math.abs(tre0ph-cph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-cth);
+ double d2ccph = Math.abs(trh0ph-cph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ } else if (ClIL>0 && ClIL<30)
+ {
+ // use maps at IL and H0
+ double d1ccth = Math.abs(trilth-cth);
+ double d1ccph = Math.abs(trilph-cph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-cth);
+ double d2ccph = Math.abs(trh0ph-cph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ } else if (ClIL>29)
+ {
+ // use maps at IL and mip cluster center
+ double d1ccth = Math.abs(trilth-cth);
+ double d1ccph = Math.abs(trilph-cph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ double d2ccth = Math.abs(clth-cth);
+ double d2ccph = Math.abs(clph-cph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ }
+ // if a core cluster was added, also check distance to it
+ double distcore = 999.;
+ if (ncore>0)
+ {
+ double dcccth = Math.abs(trcclth[0]-cth);
+ double dcccph = Math.abs(trcclph[0]-cph);
+ if (dcccph>Math.PI) dcccph = 2*Math.PI-dcccph;
+ distcore = Math.sqrt(dcccth*dcccth+dcccph*dcccph);
+ }
+ // take smallest distance to test
+ double distt = 999.;
+ if (dist1<dist2 && dist1<distcore) distt = dist1;
+ if (dist2<dist1 && dist2<distcore) distt = dist2;
+ if (distcore<dist1 && distcore<dist2) distt = distcore;
+ if (distt<niter*_mindist)
+ {
+ nclmtch++;
+ }
+
+ if (nclmtch>0 || drat>0.20) // its a match
+ {
+ trclus.addCluster(ishcl);
+ trshclus.addCluster(ishcl);
+ if (trshdb) aida.cloud2D("Track Shower Match Cluster Theta Phi").fill(cph,cth);
+ showclE += ishcl.getEnergy();
+ shcl.remove();
+// System.out.println("Cluster matched to track");
+ }
+ } // loop over all clusters, adding up energy in this road
+ TrClE += showclE;
+ if (trshdb && TrClE>0) aida.cloud1D("E over P Track Shower matches").fill(TrClE/TrP);
+ if (trshdb && TrClE/TrP>_maxeop)
+ {
+ aida.cloud1D("Number of hits in too large matched shower cluster").fill(trclus.getCalorimeterHits().size());
+ aida.cloud1D("E over p for too large match").fill(TrClE/TrP);
+ }
+ if (niter == _nloops) break;
+ } while (TrClE/TrP<_mineop);
+ if (trshdb && TrClE>0) aida.cloud1D("Number of iterations per Track").fill(niter);
+ // add shower clusters, make map
+ trkcalclus.add(trclus);
+ trkshoclus.add(trshclus);
+ trkcorclus.add(trcorclus);
+ trkclusmap.put(itr, trclus);
+ if (TrClE>0) aida.cloud1D("Final E over p Track Shower matches").fill(TrClE/TrP);
+ if (trshdb && TrClE>0) aida.cloud1D("Final E of Track Cal Matches").fill(TrClE);
+ }
+ if (trkcalclus.size()>0) event.put(_oclname,trkcalclus);
+ if (trkshoclus.size()>0) event.put("TrackShowerClusters",trkshoclus);
+ if (trkcorclus.size()>0) event.put("TrackCoreClusters",trkcorclus);
+ if (showclus.size()>0) event.put("ShowClustersLeft",showclus);
+ event.put("TrackClusMap",trkclusmap);
+
+ // add section for breaking up large clusters when E/p > max value
+ CoreReclusterer cr = new CoreReclusterer();
+ for (Track jtr : trkclusmap.keySet())
+ {
+ // for this track, get E of cluster / p of track
+ double[] tmc = jtr.getMomentum();
+ double tmp = Math.sqrt(tmc[0]*tmc[0]+tmc[1]*tmc[1]+tmc[2]*tmc[2]);
+ BasicCluster tcl = trkclusmap.get(jtr);
+ double clE = tcl.getEnergy();
+ double eprat = clE/tmp;
+ if (trshdb) aida.cloud1D("E over P from Track Cluster Map").fill(eprat);
+ if (eprat>_maxeop)
+ {
+ // re-cluster hits in this cluster
+ if (trshdb) aida.cloud1D("E over P for recluster candidates").fill(eprat);
+ Cluster oclus = (Cluster) trkclusmap.get(jtr);
+ if (trshdb) aida.cloud1D("Hits in Cluster").fill(oclus.getSize());
+ if (trshdb) aida.cloud1D("Number of clusters in TCCluster").fill(oclus.getClusters().size());
+ for (Cluster tccls : oclus.getClusters())
+ {
+ if (trshdb) aida.cloud1D("Ratio of tcclusterE to Track P").fill(tccls.getEnergy()/tmp);
+ if (tccls.getEnergy()/tmp>_maxeop)
+ {
+ //try to recluster this cluster
+ List<Cluster> screclus = cr.reclusterCluster(tccls);
+ if (trshdb) aida.cloud1D("Number of reclusters for sing clus").fill(screclus.size());
+ if (screclus.size()>0)
+ {
+ for (Cluster scrc : screclus)
+ {
+ if (trshdb) aida.cloud1D("Ratio of screclusE and Track P").fill(scrc.getEnergy()/tmp);
+ }
+ }
+ }
+ }
+ List<Cluster> nreclus = cr.reclusterCluster(oclus);
+ if (trshdb) aida.cloud1D("Number of reclusters").fill(nreclus.size());
+ if (nreclus.size()>0)
+ {
+ event.put("reclusters",nreclus);
+ for (Cluster reclus : nreclus)
+ {
+ if (trshdb) aida.cloud1D("Number of Hits in reclusters").fill(reclus.getSize());
+ if (trshdb) aida.cloud1D("Energy of reclusters").fill(reclus.getEnergy());
+ if (trshdb) aida.cloud1D("Ratio of reclusE and Track P").fill(reclus.getEnergy()/tmp);
+ }
+ }
+ }
+ }
+ }
+
+ public void setClusterNames(String[] names)
+ {
+ _clusternames = names;
+ }
+ public void setOutputClusterList(String outclname)
+ {
+ _oclname = outclname;
+ }
+}
+