lcsim/src/org/lcsim/contrib/SteveMagill
diff -N TSCForDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TSCForDriver.java 20 Aug 2008 20:20:44 -0000 1.1
@@ -0,0 +1,456 @@
+package org.lcsim.contrib.SteveMagill;
+
+// This driver associates calorimeter clusters to extrapolated tracks, testing both
+// cluster-track distance and E/p, iterating in cone around track until E/p min is reached
+// This is a slightly modified version of TSCDriver, relying on the correct determination of
+// mip cluster by TrackMipClusterDriver, only looking forward from the IL for showers, thereby
+// reducing the volume of the detector that is searched.
+// input : calorimeter clusters
+// output : clusters associated to tracks
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.Driver;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.aida.*;
+import hep.aida.*;
+import org.lcsim.spacegeom.*;
+
+public class TSCForDriver extends Driver
+{
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IAnalysisFactory analysisFactory = IAnalysisFactory.create();
+ private AIDA aida = AIDA.defaultInstance();
+ private double _dminE;
+ private double _dminH;
+ private double _dminTrC;
+ private double _mindist;
+ private int _nloops;
+ private double _mineop;
+ private double _maxeop;
+ private double _mincor;
+ private double _hres;
+ private double _mintrp;
+ private String[] _clusternames;
+ private String _oclname;
+ private boolean trshdb = false;
+
+ public TSCForDriver(double mincor, double mindist, int nloops, double hadres, double mintrp)
+ {
+ // add arguments if needed
+ _mincor = mincor; // minimum distance to core (default 0.01)
+ _mindist = mindist; // minimum distance between objects to merge (track, clusters) (0.0185)
+ _nloops = nloops; // number of iterations to merge clusters (until min E/p is met) (5)
+ _hres = hadres; // stochastic term alpha for hadron showers (0.7 rpc, 0.5 scin)
+ _mintrp = mintrp; // track momentum which determines minimum E/p value
+ // note min E/p = 1.0-_hres/sqrt(mintrp) this value is used for tracks with p<_mintrp
+// _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
+
+ // 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());
+ if (trshdb) 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());
+ 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> trilmap = (Map<Track, SpacePoint>) event.get("TrackILPosMap");
+// 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; // gets sum of all track momentum in event
+ int nmip = 0;
+ // loop over all tracks from the mip map
+ for (Track itr : trmipmap.keySet())
+ {
+ double TrClE = 0;
+ BasicCluster trclus = new BasicCluster(); // includes mips, cores, and showers
+ BasicCluster trshclus = new BasicCluster(); // cores and showers
+ BasicCluster trcorclus = new BasicCluster(); // only cores
+ int niter = 0;
+
+ // get some track properties
+ double TrP = Math.sqrt(itr.getPX()*itr.getPX()+itr.getPY()*itr.getPY()+itr.getPZ()*itr.getPZ());
+ TotTrP += TrP;
+ double Trpt = Math.sqrt(itr.getPX()*itr.getPX()+itr.getPY()*itr.getPY());
+
+ // define mineop and maxeop based on cal sigma and min track momentum
+ if (TrP<_mintrp)
+ {
+ _mineop = 1.0-(1.25*_hres/Math.sqrt(_mintrp));
+ _maxeop = 1.0+(2.6*_hres/Math.sqrt(_mintrp));
+ aida.cloud1D("sigma for eop").fill(_hres*Math.sqrt(_mintrp));
+ } else
+ {
+ _mineop = 1.0-(1.25*_hres/Math.sqrt(TrP));
+ _maxeop = 1.0+(2.6*_hres/Math.sqrt(TrP));
+ aida.cloud1D("sigma for eop").fill(_hres*Math.sqrt(TrP));
+ }
+ aida.cloud1D("Min and Max EoP").fill(_mineop);
+ aida.cloud1D("Min and Max EoP").fill(_maxeop);
+
+ // get associated mip cluster coordinates, add mip cluster to track cluster, from mip map
+ BasicCluster trmclus = trmipmap.get(itr);
+ trclus.addCluster(trmclus); // add mip cluster hits to track cluster
+ if (trmclus.getSize()>0) nmip++;
+ // this should be correct energy based on dEdx for mips - set in TrackMipClusterDriver
+ TrClE = trmclus.getEnergy();
+ double[] trmpos = trmclus.getPosition();
+ double trxym = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
+ double mclth = Math.atan(trxym/trmpos[2]);
+ if (mclth<0) mclth+=Math.PI;
+ double mclph = Math.atan2(trmpos[1],trmpos[0]);
+ if (mclph<0) mclph+=2*Math.PI;
+ double trxyzm = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]+trmpos[2]*trmpos[2]);
+
+ // get associated IL coordinates from IL map
+ SpacePoint trIL = trilmap.get(itr);
+ double trxyzil = trIL.rxyz(); // rxyz of track spacepoint at IL
+ double trtheta = trIL.theta(); // track theta at IL
+ double trphi = trIL.phi(); // track phi at IL
+ double trilth = Math.atan(trIL.rxy()/trIL.z());
+ if (trilth<0) trilth+=Math.PI;
+ double trilph = Math.atan2(trIL.y(),trIL.x());
+ if (trilph<0) trilph+=2*Math.PI;
+
+ // some histograms
+ if (trshdb)
+ {
+ aida.cloud1D("Theta of mip cluster c").fill(mclth);
+ aida.cloud1D("Theta of Track at IL").fill(trilth);
+ aida.cloud1D("Phi of Track at IL").fill(trilph);
+ aida.cloud1D("Phi of mip cluster c").fill(mclph);
+ aida.cloud2D("Theta vs Phi of Track at IL").fill(trilph,trilth);
+ aida.cloud2D("Theta vs Phi of mip cluster c").fill(mclph,mclth);
+ aida.cloud1D("E over P after adding Mip cluster").fill(TrClE/TrP);
+ }
+
+ // first, find clusters very close to track - core clusters using same distance criteria as for mips
+ // this is essentially the mip finder on clusters with no density rquirement
+ int ncore = 0;
+ double[] trcclth = new double[100];
+ double[] trcclph = new double[100];
+ double[] trcclrxyz = new double[100];
+ for (Iterator<BasicCluster> shcl = showclus.iterator(); shcl.hasNext();)
+ {
+ BasicCluster ibcl = shcl.next();
+ // first, check cluster rxyz and compare to trIL rxyz, only considering if
+ // cluster rxyz>trIL rxyz - attaches clusters beyond IL
+ double[] ccpos = ibcl.getPosition();
+ double clrxyz = Math.sqrt(ccpos[0]*ccpos[0]+ccpos[1]*ccpos[1]+ccpos[2]*ccpos[2]);
+ if (TrP>5. && trxyzil>clrxyz) continue; // don't consider this cluster
+ // don't consider if too large in E/p
+ // System.out.println(" Ratio for this cluster is " + (TrClE+ibcl.getEnergy())/TrP);
+// if ((TrClE+ibcl.getEnergy())/TrP > _maxeop) continue;
+ 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 distance for test
+ // uses distance of cluster to IL spacepoint
+ double dist = 999;
+ double dccth = Math.abs(trilth-ccth);
+ double dccph = Math.abs(trilph-ccph);
+ if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
+ dist = Math.sqrt(dccth*dccth+dccph*dccph);
+ // test distance
+ if (dist<_mincor)
+ {
+ TrClE += ibcl.getEnergy();
+ trcclth[ncore] = ccth;
+ trcclph[ncore] = ccph;
+ trcclrxyz[ncore] = clrxyz;
+ ncore++;
+ // add cluster to cands, remove from cluster list
+ trcorclus.addCluster(ibcl);
+ trshclus.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(TrClE/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; // sum of all cluster match energies
+ for (Iterator<BasicCluster> shcl = showclus.iterator(); shcl.hasNext();)
+ {
+ BasicCluster ishcl = shcl.next();
+ // first, check cluster rxyz and compare to trIL rxyz, only considering if
+ // cluster rxyz>trIL rxyz
+ double[] clpos = ishcl.getPosition();
+ double clrxyz = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]+clpos[2]*clpos[2]);
+ if (TrP>5. && trxyzil>clrxyz) continue; // if cluster closer than mip endpoint, don't consider
+ // if this cluster has E/p > _maxeop, then no point in considering it for match
+ if ((TrClE+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 dist = 999.;
+ double dccth = Math.abs(trilth-htth);
+ double dccph = Math.abs(trilph-htph);
+ if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
+ dist = Math.sqrt(dccth*dccth+dccph*dccph);
+
+ // if a core cluster was added, also check distances 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 (dist<distcore) distt = dist;
+ if (distcore<dist) distt = distcore;
+ if (trshdb) aida.cloud1D("Distance of shower hit from track").fill(distt);
+ if (distt<niter*_mindist)
+ {
+ nhmtch++;
+ }
+ }
+ double dnhits = nhits;
+ double dnhmtch = nhmtch;
+ double drat = 0.;
+ if (nhits>0 && nhmtch>0) 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 dist = 999;
+ double dccth = Math.abs(trilth-cth);
+ double dccph = Math.abs(trilph-cph);
+ if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
+ dist = Math.sqrt(dccth*dccth+dccph*dccph);
+
+ // 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 (dist<distcore) distt = dist;
+ if (distcore<dist) distt = distcore;
+ if (trshdb) aida.cloud1D("Distance of shower cluster from track").fill(distt);
+ if (distt<niter*_mindist)
+ {
+ nclmtch++;
+ }
+
+ if (nclmtch>0 || drat>0.15) // its a match
+ {
+ trclus.addCluster(ishcl);
+ trshclus.addCluster(ishcl);
+ if (trshdb) aida.cloud2D("Track Shower Match Cluster Theta Phi").fill(cph,cth);
+ if (trshdb) aida.cloud1D("Distance of matched shower cluster from track").fill(distt);
+ if (trshdb && ncore==0) aida.cloud1D("Dist MatShCl from Tr noCore").fill(distt);
+ if (trshdb) aida.cloud2D("Dist of matched ShCl from Tr vs P").fill(TrP,distt);
+ showclE += ishcl.getEnergy();
+ TrClE += ishcl.getEnergy();
+ shcl.remove();
+// System.out.println("Cluster matched to track");
+ }
+ } // loop over all clusters, adding up energy in this road
+ 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 && niter<_nloops);
+ if (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 (TrClE>0 && trshdb && ncore==0) aida.cloud1D("Final EoP nocores").fill(TrClE/TrP);
+ if (TrClE>0 && trshdb && nmip==0) aida.cloud1D("Final EoP nomips").fill(TrClE/TrP);
+ if (TrClE>0) aida.cloud2D("Final E over p vs p TSCFor").fill(TrP,TrClE/TrP);
+ if (TrClE>0) aida.cloud2D("Final E vs p of Track Cal Matches").fill(TrP,TrClE);
+ if (TrClE>TrP) aida.cloud1D("Num sigmas E gt p").fill((TrClE-TrP)/(0.7*Math.sqrt(TrP)));
+ if (TrClE==0 && trshdb) aida.cloud1D("Track P for no match").fill(TrP);
+ }
+ if (trkcalclus.size()>0) event.put(_oclname,trkcalclus);
+ if (trkshoclus.size()>0) event.put("TrackShowerClusters",trkshoclus);
+ if (trkcorclus.size()>0) event.put("TrackCoreClusters",trkcorclus);
+ List<Track> evtracks = event.get(Track.class, "PerfectTracks");
+// if (evtracks.size()==1)
+// {
+ if (showclus.size()>0)
+ {
+ event.put("ShowClustersLeft",showclus);
+ aida.cloud1D("Num of shower clusters after TSC match").fill(showclus.size());
+ int nchit = 0;
+ int nphit = 0;
+ int nnhit = 0;
+ int nhlmin = 0;
+ for (BasicCluster shclus : showclus)
+ {
+ List<CalorimeterHit> calhits = shclus.getCalorimeterHits();
+ for (CalorimeterHit calhit : calhits)
+ {
+ SimCalorimeterHit mchit = (SimCalorimeterHit) calhit;
+ double mch = Math.abs(mchit.getMCParticle(0).getCharge());
+ double hitm = mchit.getMCParticle(0).getMass();
+ if (shclus.getSize()>5)
+ {
+ if (mch>0) nchit++;
+ if (mch==0 && hitm==0) nphit++;
+ if (mch==0 && hitm>0) nnhit++;
+ }
+ }
+ }
+ aida.cloud1D("Num ChHits in leftover showclus 5 hits").fill(nchit);
+ aida.cloud1D("Num PhoHits in leftover showclus 5 hits").fill(nphit);
+ aida.cloud1D("Num NHHits in leftover showclus 5 hits").fill(nnhit);
+ }
+// }
+ event.put("TrackClusMap",trkclusmap);
+
+ // add section for breaking up large clusters when E/p > max value
+ // not used now
+ 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;
+ }
+}
+
