 |
 |
| Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill on MAIN | |||
| TrShowerDriver.java | +435 | added 1.1 | |
diff -N TrShowerDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrShowerDriver.java 14 Feb 2012 19:26:24 -0000 1.1 @@ -0,0 +1,435 @@
+ package org.lcsim.contrib.SteveMagill;
+
+// Modified from TSCForDriver, this algorithm matches clusters to track extrapolations,
+// iterating in distance from the extrapolated track, using E/p to limit the size
+// of the cluster
+// 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 TrShowerDriver extends Driver
+{
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IAnalysisFactory analysisFactory = IAnalysisFactory.create();
+ private AIDA aida = AIDA.defaultInstance();
+ private double _mindist;
+ private int _nloops;
+ private double _mineop;
+ private double _maxeop;
+ private double _hres;
+ private double _mintrp;
+ private String _inclusname;
+ private String _oclname;
+ private String _tmname;
+ private String _tcname;
+ private String _tilname;
+ private String _tispname;
+ private String _tiptname;
+ private String _tfcname;
+ private boolean trshdb = true;
+
+ public TrShowerDriver(double mindist, int nloops, double hadres, double mintrp)
+ {
+ // add arguments if needed
+ _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
+ }
+
+ protected void process(EventHeader event)
+ {
+ // get list of cal clusters
+ List<BasicCluster> inclus = event.get(BasicCluster.class,_inclusname);
+ // make copy to preserve original list
+ List<BasicCluster> copyclus = new ArrayList<BasicCluster>();
+ for (BasicCluster icl : inclus)
+ {
+ BasicCluster shclus = new BasicCluster();
+ if (icl.getSize()>0)
+ {
+ shclus.addCluster(icl);
+ copyclus.add(shclus);
+ }
+ }
+ // create lists of track-cal associated clusters for display
+ List<BasicCluster> trkcalclus = new ArrayList<BasicCluster>();
+ List<BasicCluster> trkshoclus = new ArrayList<BasicCluster>();
+ List<BasicCluster> unmcalclus = new ArrayList<BasicCluster>();
+
+ // Make a map to link track and final cluster association
+ Map<Track, BasicCluster> trkclusmap = new HashMap<Track, BasicCluster>();
+ // and a list of the tracks in the map
+ List<Track> trkcaltrack = new ArrayList<Track>();
+
+ // Get maps linking IL and mip clusters to tracks, tracks linked to SpacePoints
+ Map<Track, BasicCluster> trmipmap = (Map<Track, BasicCluster>) event.get(_tmname);
+ Map<Track, BasicCluster> trcoremap = (Map<Track, BasicCluster>) event.get(_tcname);
+// Map<BasicCluster, Integer> clusilmap = (Map<BasicCluster, Integer>) event.get("MipClusILMap");
+ Map<Track, SpacePoint> trilmap = (Map<Track, SpacePoint>) event.get(_tilname);
+ Map<Track, BasicCluster> trilspmap = (Map<Track, BasicCluster>) event.get(_tispname);
+ Map<Track, BasicCluster> triptrmap = (Map<Track, BasicCluster>) event.get(_tiptname);
+// Map<Track, SpacePoint> trh0map = (Map<Track, SpacePoint>) event.get("TrackXH0Map");
+
+ double TotTrP = 0; // gets sum of all track momentum in event
+ int ntracks = 0; // count tracks in event
+ int nmip = 0;
+ int ncore = 0;
+ int nilsp = 0;
+ int niptr = 0;
+ // loop over all tracks from the mip map
+ for (Track itr : trmipmap.keySet())
+ {
+ boolean abort = false;
+ double TrClE = 0;
+ BasicCluster trclus = new BasicCluster(); // includes mips, cores, and showers
+ BasicCluster trshclus = new BasicCluster(); // just showers
+ int niter = 0;
+ ntracks++;
+
+ // 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
+ // used 1.25 for mineop and 2.6 for maxeop originally
+ if (TrP<_mintrp)
+ {
+ _mineop = 1.0-(1.00*_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.00*_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++;
+ 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]);
+ TrClE += trmclus.getEnergy();
+
+ // 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;
+
+ // add core clusters and get coordinates of core
+ double trcclth = 0.;
+ double trcclph = 0.;
+ if (trcoremap.get(itr) != null)
+ {
+ BasicCluster coreclus = trcoremap.get(itr);
+// if ((TrClE+coreclus.getEnergy())/TrP > _maxeop) abort = true;
+// if (abort) continue;
+ trclus.addCluster(coreclus); // add no matter what test later
+ if (coreclus.getSize()>0)
+ {
+ ncore++;
+ double[] trcpos = coreclus.getPosition();
+ double trcxym = Math.sqrt(trcpos[0]*trcpos[0]+trcpos[1]*trcpos[1]);
+ trcclth = Math.atan(trcxym/trcpos[2]);
+ if (trcclth<0) trcclth+=Math.PI;
+ trcclph = Math.atan2(trcpos[1],trcpos[0]);
+ if (trcclph<0) trcclph+=2*Math.PI;
+ double trcxyzm = Math.sqrt(trcpos[0]*trcpos[0]+trcpos[1]*trcpos[1]+trcpos[2]*trcpos[2]);
+ TrClE += coreclus.getEnergy();
+ }
+ }
+// if (abort) continue; // go to next track if this is true
+
+ // add ILSP clusters - don't need coord info for these
+ if (trilspmap.get(itr) != null)
+ {
+ BasicCluster ilspclus = trilspmap.get(itr);
+// if ((TrClE+ilspclus.getEnergy())/TrP > _maxeop) abort = true;
+// if (abort) continue;
+ trclus.addCluster(ilspclus);
+ if (ilspclus.getSize()>0)
+ {
+ nilsp++;
+ TrClE += ilspclus.getEnergy();
+ }
+ }
+// if (abort) continue;
+
+ // add IPTrack clusters - don't need coord info for these
+ if (triptrmap.get(itr) != null)
+ {
+ BasicCluster iptrclus = triptrmap.get(itr);
+// if ((TrClE+iptrclus.getEnergy())/TrP > _maxeop) abort = true;
+// if (abort) continue;
+ trclus.addCluster(iptrclus);
+ if (iptrclus.getSize()>0)
+ {
+ niptr++;
+ TrClE += iptrclus.getEnergy();
+ }
+ }
+// if (abort) continue;
+
+ // 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 Core ILSP").fill(TrClE/TrP);
+ }
+
+ // check all remaining 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.
+ do {
+ niter++;
+ double showclE = 0; // sum of all cluster match energies
+ for (Iterator<BasicCluster> shcl = copyclus.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 makes 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 from ILSP
+ 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-htth);
+ double dcccph = Math.abs(trcclph-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;
+
+ // now do the same for clusters, not hits
+ int nclmtch = 0;
+ 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-cth);
+ double dcccph = Math.abs(trcclph-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.30) // 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) 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, these numbers have been tuned for single pions
+ // first look - use 0.86 - 1.16
+ if (TrClE/TrP>0.86 && TrClE/TrP<1.16)
+ {
+ trkcalclus.add(trclus);
+ if (trshclus.getSize()>0) trkshoclus.add(trshclus);
+ trkclusmap.put(itr, trclus);
+ trkcaltrack.add(itr);
+ } else
+ {
+ if (trclus.getSize()>0) unmcalclus.add(trclus);
+ }
+ if (TrClE>0 && trshdb) aida.cloud1D("Final E over p Track Shower matches").fill(TrClE/TrP);
+ if (TrClE>0 && trshdb && nmip==0) aida.cloud1D("Final EoP nomiphits").fill(TrClE/TrP);
+ if (TrClE>0 && trshdb) aida.cloud2D("Final E over p vs p TrShower").fill(TrP,TrClE/TrP);
+ if (TrClE>0 && trshdb) aida.cloud2D("Final E vs p of Track Cal Matches").fill(TrP,TrClE);
+ if (TrClE>TrP && trshdb) aida.cloud1D("Number sigmas E gt p").fill((TrClE-TrP)/(_hres*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 (trkcalclus.size()>0) event.put("TrCalTracks",trkcaltrack);
+ if (trkshoclus.size()>0) event.put("TrackShowerClusters",trkshoclus);
+// if (trkcalclus.size()>0) event.put(_tfcname,trkclusmap);
+ event.put(_tfcname,trkclusmap); // write out even if no entries
+ event.put("UnAssocCalClusters",unmcalclus);
+// List<Double> qqpT = event.get(Double.class, "DJpT");
+ aida.cloud1D("Number of Track Cluster Matches per Event").fill(trkcalclus.size());
+ double ntm = (double) trkcalclus.size();
+ double ntt = (double) ntracks;
+ double tmfrac = 0.;
+ if (ntt>0.) tmfrac = ntm/ntt;
+// double jpT = 0.;
+// for (Double qp : qqpT)
+// {
+// jpT = qp.doubleValue();
+// }
+// if (jpT>0.) aida.cloud2D("Num TrCAL Matches vs pT").fill(jpT, trkcalclus.size());
+// if (jpT>0. && tmfrac>0.) aida.cloud2D("FracTrCAL Matches vs ET").fill(jpT,tmfrac);
+
+ if (copyclus.size()>0)
+ {
+ event.put("ShowerClustersLeft",copyclus);
+ aida.cloud1D("Number Clusters after TSC match").fill(copyclus.size());
+ int nchit = 0;
+ int nphit = 0;
+ int nnhit = 0;
+ for (BasicCluster shclus : copyclus)
+ {
+ 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++;
+ }
+ }
+ }
+ if (trshdb)
+ {
+ 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);
+ }
+ }
+ }
+
+ public void setInputClusterList(String inname)
+ {
+ _inclusname = inname;
+ }
+ public void setOutputClusterList(String outclname)
+ {
+ _oclname = outclname;
+ }
+ public void setTrackMipClusMap(String tmname)
+ {
+ _tmname = tmname;
+ }
+ public void setTrackILPosMap(String tilname)
+ {
+ _tilname = tilname;
+ }
+ public void setTrackCoreClusMap(String tcname)
+ {
+ _tcname = tcname;
+ }
+ public void setTrackILSPClusMap(String tispname)
+ {
+ _tispname = tispname;
+ }
+ public void setTrackIPTrClusMap(String tiptname)
+ {
+ _tiptname = tiptname;
+ }
+ public void setTrackCalClusMap(String tfcname)
+ {
+ _tfcname = tfcname;
+ }
+}
+
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1