lcsim/src/org/lcsim/contrib/SteveMagill
diff -u -r1.2 -r1.3
--- TrMipCompDriver.java 29 May 2008 22:36:47 -0000 1.2
+++ TrMipCompDriver.java 20 Aug 2008 20:20:08 -0000 1.3
@@ -15,7 +15,8 @@
import org.lcsim.event.Track;
import org.lcsim.util.hitmap.HitMap;
import org.lcsim.recon.cluster.util.BasicCluster;
-
+import org.lcsim.spacegeom.*;
+import org.lcsim.event.base.*;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.base.MCReconstructedParticle;
import org.lcsim.event.ParticleID;
@@ -39,14 +40,59 @@
private String _tmclusname;
private AIDA aida = AIDA.defaultInstance();
private boolean perfPFAD = true;
+ private boolean muons = false;
public void process(EventHeader event)
{
- List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
+ // 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("TrackMipEPMap");
+
+ // 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);
+ 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 && perfrps.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 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);
@@ -66,22 +112,32 @@
{
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++;
}
if (numhits>0)
{
double dnumhits = numhits;
double dnumchits = numchits;
double hrat = dnumchits/dnumhits;
- aida.cloud2D("NClus hits vs NMC hits Tr Mips").fill(numchits,numhits);
- aida.cloud1D("MC Particle Purity per Tr Mip Clus").fill(hrat);
+ 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;
- if (nhitstm>0) aida.cloud1D("MC Particle Purity for all Tr Mip Clus").fill(mipTot);
+ 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)
@@ -94,7 +150,11 @@
aida.cloud1D("Num Hits Mip Clus single cluster").fill(nhitstm);
}
- // Check the cheated clusters also - only good for muons
+ // Check the cheated clusters also - only good for muons
+ if (muons)
+ {
+ List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
+
double chclusE = 0;
int NCheatCl = 0;
double PerfChHClE = 0;
@@ -123,6 +183,7 @@
aida.cloud1D("Num Hits in Cheated Cluster single cluster Mips").fill(nhitsingcl);
if (NMipCl == 1) aida.cloud1D("Diff NHits TM Cheat Clusters").fill(nhitstm-nhitsingcl);
}
+ }
}
public void setTrMipClusName(String tmclusname)