lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N GTrMipCompDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ GTrMipCompDriver.java 12 Aug 2009 18:30:39 -0000 1.1
@@ -0,0 +1,158 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.util.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.Track;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.spacegeom.*;
+import org.lcsim.event.base.*;
+import org.lcsim.event.MCParticle;
+
+/**
+ * Defines perfect particle flow objects - perfect charged particles, perfect photons and perfect neutral hadrons
+ */
+
+public class GTrMipCompDriver extends Driver
+{
+ public GTrMipCompDriver() {
+ }
+
+ private String _tmclusname;
+ private AIDA aida = AIDA.defaultInstance();
+ private boolean perfPFAD = true;
+ private boolean muons = false;
+
+ public void process(EventHeader event)
+ {
+
+ // compare mip endpoint and MC track endpoint
+// List<ReconstructedParticle> perfrps = event.get(ReconstructedParticle.class,"PerfectRecoParticles");
+ List<Track> evtracks = event.get(Track.class, "PerfectTracks");
+ Map<Track, SpacePoint> trMEPmap = (Map<Track, SpacePoint>) event.get("TrackILPosMap");
+ Map<MCParticle, Track> mcptr = (Map<MCParticle, Track>) event.get("MCPTrMap");
+ Map<Track, SpacePoint> treSP = (Map<Track, SpacePoint>) event.get("TrackEndPointSP");
+
+ // check endpoints of tracks to compare to IL finding
+ for (Track itrack : evtracks)
+ {
+ double dedx = itrack.getdEdx();
+ double[] p = itrack.getMomentum();
+ double ptot = Math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+ double pt = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
+// aida.cloud2D("dEdx vs P").fill(ptot,dedx);
+// BaseTrackMC ptrac = (BaseTrackMC) itrack;
+// Hep3Vector endp = ptrac.getMCParticle().getEndPoint();
+// System.out.println(" Track EndPt " + endp);
+ SpacePoint endpSP = treSP.get(itrack);
+ Hep3Vector endp = new SpacePoint(endpSP);
+ double endxyz = Math.sqrt(endp.x()*endp.x()+endp.y()*endp.y()+endp.z()*endp.z());
+ double endxy = Math.sqrt(endp.x()*endp.x()+endp.y()*endp.y());
+ double cosepth = endxy/endxyz;
+
+ SpacePoint MEPsp = trMEPmap.get(itrack);
+ double mepxyz = MEPsp.rxyz();
+ double mepxy = MEPsp.rxy();
+// System.out.println("Tr EP" + endxyz + "IL SP" + mepxyz);
+ double trildif = Math.abs(endxyz-mepxyz);
+ double meptrep = mepxyz-endxyz;
+ double trmepdif = Math.abs(endxy-mepxy);
+ if (evtracks.size()==1)
+ {
+ if (endxy>1270. || Math.abs(endp.z())>1680.) aida.cloud1D("Track MEPxyz Difference in mm").fill(trildif);
+ if (endxy>1270. || Math.abs(endp.z())>1680.) aida.cloud1D("Track MEPxy Diff in mm").fill(trmepdif);
+ if (endxy>1270. || Math.abs(endp.z())>1680.) aida.cloud1D("MEP TrEP Diff in mm").fill(meptrep);
+ if (endxy>1270. || Math.abs(endp.z())>1680.) aida.cloud2D("MEP TrEP vs TrP").fill(ptot,meptrep);
+ if (endxy>1270. || Math.abs(endp.z())>1680.) aida.cloud2D("MEP TrEP vs Trpt").fill(pt,meptrep);
+// System.out.println(" track endpt " + endxyz + " track pt " + pt);
+ if (endxy>1270. || Math.abs(endp.z())>1680.) aida.cloud2D("Mip endpoint vs track pt").fill(pt,mepxyz);
+ aida.cloud1D("Track endpoint z").fill(endp.z());
+ if (Math.abs(endp.z())<1500.) aida.cloud1D("Track endpoint rxy in barrel").fill(endxy);
+ if (Math.abs(endp.z())<1500. && endxy>1250.) aida.cloud1D("Modified Tr endprxy").fill((endxy-1269.25)/3.5);
+ if (Math.abs(endp.z())<1500. && endxy>1250.) aida.cloud1D("Modified mip endprxy").fill((mepxy-1269.25)/3.5);
+ if (Math.abs(endp.z())<1500.) aida.cloud1D("Track endpoint in barrel").fill(endxyz);
+ }
+ }
+
+ int NMipCl = 0;
+ int nhitstm = 0; // number of hits in all mip clusters
+ int nchhitstm = 0; // number of charged particle hits in all mip clusters
+ int nnhhitstm = 0; // number of neutral hadron hits in mip clusters
+ int nphohitstm = 0; // number of photon hits in mip clusters
+ try
+ {
+ List<BasicCluster> tmclusters = event.get(BasicCluster.class,_tmclusname);
+// Map<Long, CalorimeterHit> scidcmap = (Map<Long, CalorimeterHit>) event.get("ScEMIDMap");
+
+ for (BasicCluster tmclus : tmclusters)
+ {
+ NMipCl++;
+ int numhits = tmclus.getSize();
+ nhitstm += tmclus.getSize(); // sum hits from all mips in event
+ double tmclE = tmclus.getEnergy();
+ aida.cloud1D("Tr Mip Clus E").fill(tmclE);
+ aida.cloud1D("Tr Mip Clus Size").fill(numhits);
+ aida.cloud2D("Tr Mip NumHits vs E").fill(tmclE,numhits);
+ List<CalorimeterHit> tmcalhits = tmclus.getCalorimeterHits();
+ int numchits = 0;
+ for (CalorimeterHit tmhit : tmcalhits)
+ {
+ SimCalorimeterHit hit = (SimCalorimeterHit) tmhit;
+ double hitch = Math.abs(hit.getMCParticle(0).getCharge());
+ double hitm = hit.getMCParticle(0).getMass();
+ if (hitch>0) numchits++;
+ if (hitch>0) nchhitstm++;
+ if (hitch==0 && hitm==0) nphohitstm++;
+ if (hitch==0 && hitm>0) nnhhitstm++;
+// double che = tmhit.getCorrectedEnergy();
+// long chid = tmhit.getCellID();
+// boolean scid = scidcmap.containsKey(chid);
+// if (scid)
+// {
+// double sce = scidcmap.get(chid).getCorrectedEnergy();
+// aida.cloud1D("C over S for EM Mip Clus Hits").fill(che/sce);
+// }
+ }
+ if (numhits>0)
+ {
+ double dnumhits = numhits;
+ double dnumchits = numchits;
+ double hrat = dnumchits/dnumhits;
+ aida.cloud2D("NClus hits vs NMC hits per TrMipClus").fill(numchits,numhits);
+ aida.cloud1D("MC Particle Purity per TrMipClus").fill(hrat);
+ }
+ }
+ double dnhitstm = nhitstm;
+ double dnchhitstm = nchhitstm;
+ double mipTot = dnchhitstm/dnhitstm;
+ double dnphohitstm = nphohitstm;
+ double dnnhhitstm = nnhhitstm;
+ if (nhitstm>0)
+ {
+ aida.cloud1D("MC Particle Purity TrMips per event").fill(mipTot);
+ aida.cloud1D("Pho conf TrMips per event").fill(dnphohitstm/dnhitstm);
+ aida.cloud1D("NeuH conf TrMips per event").fill(dnnhhitstm/dnhitstm);
+ }
+
+ }
+ catch(java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("requested object not found in event " + _tmclusname);
+ }
+
+ if (NMipCl == 1)
+ {
+ aida.cloud1D("Num Hits Mip Clus single cluster").fill(nhitstm);
+ }
+ }
+
+ public void setTrMipClusName(String tmclusname)
+ {
+ _tmclusname = tmclusname;
+ }
+
+}