lcsim/src/org/lcsim/contrib/SteveMagill
diff -N TrackShowerClusterDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackShowerClusterDriver.java 23 Apr 2007 20:08:53 -0000 1.1
@@ -0,0 +1,434 @@
+package org.lcsim.contrib.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.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;
+
+public class TrackShowerClusterDriver 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[] _emBRadii;
+ private double[] _hadBRadii;
+ private double[] BRadii = new double[100];
+ private double[] ECZs = new double[100];
+// modified Barrel helix swimmer for org.lcsim test
+// private Helix[] _emBSwimmers; //EM
+// private Helix[] _hadBSwimmers; //HAD
+// 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
+// The central magnetic field strength
+ private double[] _fieldStrength;
+ private double _dminE;
+ private double _dminH;
+ private double _dminTrC;
+ private String[] _clusternames;
+// private String[] _hitmapnames;
+
+ public TrackShowerClusterDriver()
+ {
+ // add arguments if needed
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event); // executes all added drivers
+
+ boolean mcfltr = true;
+// mcfltr = (Boolean)event.get("MCFilt");
+// System.out.println("MCFilt is " +mcfltr);
+ if (mcfltr)
+ {
+
+// Initialize things here
+ if(!_initialized)
+ {
+ // setup specific detector stuff here
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ _fieldStrength = det.getFieldMap().getField(zero);
+// Hep3Vector h3zero = new BasicHep3Vector(zero);
+// Hep3Vector field = det.getFieldMap().getField(h3zero);
+// System.out.println("B Field " +_fieldStrength[2]+ " 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>();
+ for (int i=0; i<_clusternames.length; i++)
+ {
+ try
+ {
+ List<BasicCluster> clus = event.get(BasicCluster.class,_clusternames[i]);
+ for (BasicCluster cl : clus)
+ {
+ BasicCluster shclus = new BasicCluster();
+ if (cl.getSize()>0)
+ {
+ 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());
+
+ // create a list of track-cal associated clusters for display
+ List<BasicCluster> trkcalclus = new ArrayList<BasicCluster>();
+ // Loop over tracks in event, extrapolate to layers and match hits if mips
+// List<PerfectTrack> evtracks = event.get(PerfectTrack.class, "PerfectTracks");
+ // Get maps linking IL and mip clusters to tracks
+ Map<PerfectTrack, BasicCluster> trmipmap = (Map<PerfectTrack, BasicCluster>) event.get("TrackMipMap");
+ Map<BasicCluster, Integer> clusilmap = (Map<BasicCluster, Integer>) event.get("MipClusILMap");
+ double TotTrP = 0;
+ for (PerfectTrack itr : trmipmap.keySet())
+ {
+ double TrClE = 0;
+ BasicCluster trclus = new BasicCluster();
+ int niter = 0;
+ // get track and associated mip cluster
+// if (trmipmap.isEmpty()) break;
+// PerfectTrack itr = trmipmap.keySet().iterator().next();
+ BasicCluster trmclus = trmipmap.get(itr);
+ trclus.addCluster(trmclus);
+ TrClE += trmclus.getEnergy();
+ Integer 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);
+ Hep3Vector tror3 = itr.getOrigin();
+ 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(_fieldStrength[2]);
+ // swim track to cluster position, arguments are momentum, origin, and charge
+ tswim.setTrack(trp3, tror3, trq);
+ double tobrad = tswim.getDistanceToRadius(BRadii[ClIL]);
+ double toecz = tswim.getDistanceToZ(ECZs[ClIL]);
+ double trtheta = 0;
+ double trphi = 0;
+ double ctrth = 0;
+ double ctrph = 0;
+ if (tobrad<toecz) // in barrel
+ {
+ SpacePoint trSP = tswim.getPointAtDistance(tobrad);
+
+ trtheta += trSP.theta();
+ trphi += trSP.phi();
+ ctrth += Math.atan(trSP.rxy()/trSP.z());
+ if (ctrth<0) ctrth+=Math.PI;
+ ctrph += Math.atan2(trSP.y(),trSP.x());
+ if (ctrph<0) ctrph+=2*Math.PI;
+// aida.cloud1D("Theta of Track at IL B").fill(trtheta);
+ aida.cloud1D("Theta of Track c at IL B").fill(ctrth); // this is correct value
+ double clr = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
+ double clth = Math.atan(clr/trmpos[2]);
+ if (clth<0) clth+=Math.PI;
+ double clph = Math.atan2(trmpos[1],trmpos[0]);
+ if (clph<0) clph+=2*Math.PI;
+ aida.cloud1D("Theta of mip cluster c at clus B").fill(Math.atan(clr/trmpos[2]));
+ aida.cloud1D("Phi of Track at IL B").fill(trphi);
+ aida.cloud1D("Phi of Track c at IL B").fill(ctrph);
+ aida.cloud1D("Phi of mip cluster c at clus B").fill(Math.atan2(trmpos[1],trmpos[0]));
+ aida.cloud2D("Theta vs Phi of Track at IL B").fill(ctrph,ctrth);
+ aida.cloud2D("Theta vs Phi of mip cluster at IL B").fill(clph,clth);
+ } else // in endcap
+ {
+ SpacePoint trSP = tswim.getPointAtDistance(toecz);
+ trtheta += trSP.theta();
+ trphi += trSP.phi();
+ ctrth += Math.atan(trSP.rxy()/trSP.z());
+ if (ctrth<0) ctrth+=Math.PI;
+ ctrph += Math.atan2(trSP.y(),trSP.x());
+ if (ctrph<0) ctrph+=2*Math.PI;
+ aida.cloud1D("Theta of Track c at IL EC").fill(ctrth);
+ double clr = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
+ double clth = Math.atan(clr/trmpos[2]);
+ if (clth<0) clth+=Math.PI;
+ double clph = Math.atan2(trmpos[1],trmpos[0]);
+ if (clph<0) clph+=2*Math.PI;
+ aida.cloud1D("Theta of mip cluster c at clus EC").fill(Math.atan(clr/trmpos[2]));
+ aida.cloud1D("Phi of Track at IL EC").fill(trphi);
+ aida.cloud1D("Phi of Track c at IL EC").fill(ctrph);
+ aida.cloud1D("Phi of mip cluster c at clus EC").fill(Math.atan2(trmpos[1],trmpos[0]));
+ aida.cloud2D("Theta vs Phi of Track at IL EC").fill(ctrph,ctrth);
+ aida.cloud2D("Theta vs Phi of mip cluster at IL EC").fill(clph,clth);
+ }
+ aida.cloud2D("TrackMip Theta vs Phi Extrap").fill(trphi,trtheta);
+ // preselection of clusters near this track - first make a list of clusters that
+ // are candidates to be associated with this track
+ List<BasicCluster> candshclus = new ArrayList<BasicCluster>();
+ for (BasicCluster ibcl : showclus)
+ {
+ 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;
+ double dtrccth = Math.abs(ctrth-ccth);
+ double dtrccph = Math.abs(ctrph-ccph);
+ if (dtrccph>Math.PI) dtrccph = 2*Math.PI-dtrccph;
+ double distrcc = Math.sqrt(dtrccth*dtrccth+dtrccph*dtrccph);
+ if (distrcc<1.0)
+ {
+ for (int i=0; i<ibcl.getCalorimeterHits().size(); i++)
+ {
+ aida.cloud1D("Distance Track Cluster in Shower Match hit wt").fill(distrcc);
+ }
+ aida.cloud1D("Distance Track Cluster in Shower Match cl wt").fill(distrcc);
+ // add cluster to cands
+ candshclus.add(ibcl);
+ }
+ }
+ aida.cloud1D("Number of Cand Clusters added to track").fill(candshclus.size());
+ event.put("CandidateShowerClusters",candshclus);
+ // 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();
+ List<CalorimeterHit> clhits = ishcl.getCalorimeterHits();
+ // double cltheta = Math.abs(ishcl.getITheta());
+ // double clphi = Math.abs(ishcl.getIPhi());
+ int nhmtch = 0;
+ for (CalorimeterHit iclhit : clhits)
+ {
+ // 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;
+ double dtrhth = Math.abs(ctrth-htth);
+ double dtrhph = Math.abs(ctrph-htph);
+ if (dtrhph>Math.PI) dtrhph = 2*Math.PI-dtrhph;
+ double distrh = Math.sqrt(dtrhth*dtrhth+dtrhph*dtrhph);
+ if (distrh<niter*0.0075) nhmtch++;
+ }
+ 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;
+ aida.cloud1D("Calculated shower cluster phi").fill(cph);
+ aida.cloud1D("Calculated shower cluster theta").fill(cth);
+ aida.cloud2D("Theta vs Phi shower cluster").fill(cph,cth);
+ double dthtrcl = Math.abs(cth-ctrth);
+ double dphtrcl = Math.abs(cph-ctrph);
+ if (dphtrcl>Math.PI) dphtrcl = 2*Math.PI-dphtrcl;
+ aida.cloud1D("Track Shower Cluster dtheta").fill(dthtrcl);
+ aida.cloud1D("Track Shower Cluster dphi").fill(dphtrcl);
+ double dist = Math.sqrt(dthtrcl*dthtrcl+dphtrcl*dphtrcl);
+ aida.cloud1D("Track Shower dist").fill(dist);
+ aida.cloud1D("Number of hits matched to Track").fill(nhmtch);
+ // if (dist < niter*0.015)
+ if (nhmtch>5)
+ {
+ trclus.addCluster(ishcl);
+ aida.cloud2D("Track Shower Match Cluster Theta Phi").fill(cph,cth);
+ aida.cloud2D("Track Shower Match Track Theta Phi").fill(ctrph,ctrth);
+ aida.cloud1D("Number of matched hits for matched CLusters").fill(nhmtch);
+ // try changing theta, phi to average of trackcluster coord if matched
+ showclE += ishcl.getEnergy();
+ ctrph = (TrP*ctrph+showclE*cph)/(TrP+showclE);
+ ctrth = (TrP*ctrth+showclE*cth)/(TrP+showclE);
+ shcl.remove();
+// System.out.println("Cluster matched to track");
+ }
+ }
+ TrClE += showclE;
+ if (TrClE>0) aida.cloud1D("E over P Track Shower matches").fill(TrClE/TrP);
+ if (niter == 6) break;
+ } while (TrClE/TrP<0.55);
+ if (TrClE>0) aida.cloud1D("Number of iterations for Track").fill(niter);
+ if (TrClE>0) aida.cloud1D("Final E over p Track Shower matches").fill(TrClE/TrP);
+ trkcalclus.add(trclus);
+ }
+ if (trkcalclus.size()>0) event.put("TrackCALClusters",trkcalclus);
+ // simple check for neutrals - based on min hits in cluster
+ List<BasicCluster> neuclus = new ArrayList<BasicCluster>();
+ for (BasicCluster shc : showclus)
+ {
+ aida.cloud1D("Number of hits in NeuHClusters").fill(shc.getCalorimeterHits().size());
+ if (shc.getCalorimeterHits().size()>20) neuclus.add(shc);
+ }
+ if (neuclus.size()>0) event.put("NeuHClusters",neuclus);
+// System.out.println("Number Shower Clusters after Matching " + showclus.size());
+
+ // make some plots
+ List<ReconstructedParticle> prepars = event.get(ReconstructedParticle.class,"PerfectRecoParticles");
+ double perfpho = 0;
+ double perfneu = 0;
+ double perfchpE = 0;
+ double perfchpP = 0;
+ for (ReconstructedParticle prepar : prepars)
+ {
+ double perfpCh = prepar.getCharge();
+// ParticleID perfpID = perfrp.getParticleIDUsed();
+ int perfpID = Math.abs(prepar.getParticleIDUsed().getPDG());
+// System.out.println("Particle ID = " + perfpID);
+ double perfpE = prepar.getEnergy();
+ Hep3Vector perfpP = prepar.getMomentum();
+ double pP = perfpP.magnitude();
+ if (perfpCh == 0)
+ {
+ if (perfpID == 22)
+ {
+ perfpho += perfpE;
+ }else
+ {
+ perfneu += perfpE;
+ }
+ } else
+ {
+ perfchpE += perfpE;
+ perfchpP += pP;
+ }
+ }
+ List<MCReconstructedParticle> mcpars = event.get(MCReconstructedParticle.class,"CheatReconstructedParticles");
+ double trphoE = 0;
+ for (MCReconstructedParticle mcpar : mcpars)
+ {
+ int mcID = Math.abs(mcpar.getParticleIDUsed().getPDG());
+ double mcE = mcpar.getEnergy();
+ if (mcID == 22) trphoE += mcE;
+ }
+ double TotTrClE = 0;
+ for (BasicCluster trccl : trkcalclus)
+ {
+ TotTrClE += trccl.getEnergy();
+ }
+ aida.cloud1D("PFA Charged Hadron ESum").fill(TotTrClE);
+ aida.cloud1D("PFA Total Track PSum").fill(TotTrP);
+ double TotRemainE = 0;
+ for (BasicCluster nc : neuclus)
+ {
+ TotRemainE += nc.getEnergy()*1.3;
+ }
+ aida.cloud1D("Neutral Cluster ESum").fill(TotRemainE);
+ List<BasicCluster> phobcl = event.get(BasicCluster.class,"PhotonBClusters");
+ List<BasicCluster> phoeccl = event.get(BasicCluster.class,"PhotonECClusters");
+ double TotPhoE = 0;
+ for (BasicCluster phob : phobcl)
+ {
+ TotPhoE += phob.getEnergy()*1.015*1.013;
+ }
+ for (BasicCluster phoec : phoeccl)
+ {
+ TotPhoE += phoec.getEnergy()*1.015*1.013;
+ }
+ aida.cloud1D("PFA Photon ESum").fill(TotPhoE);
+ aida.cloud1D("PFA Track P + PhotonE").fill(TotTrP+TotPhoE);
+ double PFAESum = TotTrP+TotPhoE+TotRemainE;
+ aida.cloud1D("PFA TrP + PhoE + RemE").fill(TotTrP+TotPhoE+TotRemainE);
+ aida.cloud1D("Difference Perf Pho PFA Pho").fill(perfpho-TotPhoE);
+ aida.cloud1D("Difference True Pho Perf Pho").fill((trphoE-perfpho)/Math.sqrt(trphoE));
+ aida.cloud2D("TrPho PerfPho vs PerfPho PFAPho").fill(trphoE-perfpho,perfpho-TotPhoE);
+ aida.cloud1D("Difference Perf Neu PFA Neu").fill(perfneu-TotRemainE);
+ aida.cloud1D("Difference Perf PFA PFA").fill(perfpho+perfneu-TotPhoE-TotRemainE);
+ aida.cloud1D("Perf Pho Neu PFA Pho Neu").fill(perfpho+perfneu-TotPhoE-TotRemainE);
+
+ } // end of filter loop
+ }
+
+ public void setClusterNames(String[] names)
+ {
+ _clusternames = names;
+ }
+}
+
