lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/TrackingTest
diff -u -r1.1 -r1.2
--- LOIEffFake.java 25 Feb 2009 17:51:32 -0000 1.1
+++ LOIEffFake.java 4 Mar 2009 20:36:33 -0000 1.2
@@ -6,9 +6,11 @@
import hep.aida.IHistogram1D;
import hep.aida.IHistogramFactory;
+import hep.physics.jet.EventShape;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
import java.util.ArrayList;
import java.util.List;
@@ -19,6 +21,7 @@
import org.lcsim.event.MCParticle;
import org.lcsim.event.RelationalTable;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.BaseRelationalTable;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelixParamCalculator;
@@ -34,6 +37,7 @@
public class LOIEffFake extends Driver {
private AIDA aida = AIDA.defaultInstance();
+ private EventShape es = new EventShape();
private IHistogram1D pTeff1;
private IHistogram1D pTeff2;
private IHistogram1D thetaeff;
@@ -42,9 +46,10 @@
private IHistogram1D d0eff2;
private IHistogram1D z0eff1;
private IHistogram1D z0eff2;
+ private IHistogram1D thrusteff;
private IHistogram1D fakes;
private IHistogram1D nfakes;
- int ntrk = 0;
+ int ntrktot = 0;
int nevt = 0;
public LOIEffFake() {
@@ -59,6 +64,7 @@
d0eff2 = hf.createHistogram1D("Efficiency vs d0 full", "", 24, -12., 12., "type=efficiency");
z0eff1 = hf.createHistogram1D("Efficiency vs z0", "", 50, -5., 5., "type=efficiency");
z0eff2 = hf.createHistogram1D("Efficiency vs z0 full", "", 24, -12., 12., "type=efficiency");
+ thrusteff = hf.createHistogram1D("Efficiency vs alpha", "", 100, 0., 90., "type=efficiency");
fakes = hf.createHistogram1D("Number of mis-matched hits (unnormalized)", "", 10, 0., 10.);
nfakes = hf.createHistogram1D("Number of mis-matched hits (normalized)", "", 10, 0., 10.);
}
@@ -84,6 +90,9 @@
if (s.getMinPT() < ptCut) ptCut = s.getMinPT();
}
+ // Get the list of MCParticles
+ List<MCParticle> mclist = event.getMCParticles();
+
// Create a relational table that maps TrackerHits to MCParticles
RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
@@ -101,12 +110,16 @@
// Analyze the tracks in the event
for (Track track : tracklist) {
+ // Increment the total number of tracks
+ ntrktot++;
+
// Calculate the track pT and cos(theta)
- double px = track.getPX();
- double py = track.getPY();
- double pz = track.getPZ();
+ Hep3Vector ptrk = new BasicHep3Vector(track.getMomentum());
+ double px = ptrk.x();
+ double py = ptrk.y();
+ double pz = ptrk.z();
double pt = Math.sqrt(px * px + py * py);
- double cth = pz / Math.sqrt(pt * pt + pz * pz);
+ double cth = pz / ptrk.magnitude();
double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
@@ -119,10 +132,7 @@
double purity = tkanal.getPurity();
aida.histogram1D("Mis-matched hits for all tracks", 10, 0., 10.).fill(nbad);
aida.histogram1D("Mis-matched hits " + nhits + " hit tracks", 10, 0., 10.).fill(nbad);
-
- // Generate a normalized histogram after 1000 events
- ntrk++;
- if (nevt <= 1000) fakes.fill(nbad);
+ fakes.fill(nbad);
// Make plots for fake, non-fake, and all tracks
if (purity < 0.5) {
@@ -197,29 +207,29 @@
}
}
- // Make the normalized fake plot after the specified number of events
- if (nevt == 1000) {
- double wgt = 1. / ntrk;
- for (int i = 0; i < 10; i++) {
- System.out.println(" Entries: " + fakes.binEntries(i) + " for mismatches: " + i);
- for (int j = 0; j < fakes.binHeight(i); j++) nfakes.fill(i, wgt);
- }
- System.out.println("Normalization: " + nfakes.sumAllBinHeights() + " after ntrk = " + ntrk);
- }
+ // Calculate the thrust axis and find it's polar angle
+ Hep3Vector taxis = Thrust(mclist);
+ double cth_thrust = taxis.z() / taxis.magnitude();
+ double thrang = 180. * Math.acos(cth_thrust) / Math.PI;
+ aida.histogram1D("Polar angle of the thrust axis", 90, 0., 180.).fill(thrang);
+ aida.histogram1D("cos(theta) of the thrust axis", 100, -1., 1.).fill(cth_thrust);
// Now loop over all MC Particles
- List<MCParticle> mclist = event.getMCParticles();
for (MCParticle mcp : mclist) {
// Calculate the pT and polar angle of the MC particle
- double px = mcp.getPX();
- double py = mcp.getPY();
- double pz = mcp.getPZ();
+ Hep3Vector pmc = mcp.getMomentum();
+ double px = pmc.x();
+ double py = pmc.y();
+ double pz = pmc.z();
double pt = Math.sqrt(px * px + py * py);
- double p = Math.sqrt(pt * pt + pz * pz);
- double cth = pz / p;
+ double cth = pz / pmc.magnitude();
double theta = 180. * Math.acos(cth) / Math.PI;
+ // Calculate the angle between the MC Particle and the thrust axis
+ double calpha = VecOp.dot(pmc, taxis) / (pmc.magnitude() * taxis.magnitude());
+ double alpha = 180. * Math.acos(Math.abs(calpha)) / Math.PI;
+
// Find the number of layers hit by this mc particle
int nhits = findable.LayersHit(mcp);
@@ -228,56 +238,40 @@
double d0 = helix.getDCA();
double z0 = helix.getZ0();
- // Check cases where we have multiple tracks associated with this MC particle
- Set<Track> trklist = trktomc.allTo(mcp);
- int ntrk = trklist.size();
-// if (ntrk > 1) {
- // Count tracks where the assigned MC particle has more than 1 hit
-// int nmulthits = 0;
-// for (Track trk : trklist) {
-// TrackAnalysis tkanal = new TrackAnalysis(trk, hittomc);
-// if (tkanal.getNBadHits() < tkanal.getNHits() - 1)
-// nmulthits++;
-// }
- // Flag any anomalous cases that we find
-// if (nmulthits > 1) {
-// System.out.println("2 tracks associated with a single MC Particle");
-// for (Track trk : trklist) System.out.println(trk.toString());
-// }
-// }
+ // Calculate the weight for the efficiency plots
+ double wgt = 0.;
+ if (trktomc.allTo(mcp).size() > 0) wgt = 1.;
// Make pT efficiency plot
if (findable.isFindable(mcp, slist, Ignore.NoPTCut)) {
- double wgt = 0.;
- if (ntrk > 0) wgt = 1.;
pTeff1.fill(pt, wgt);
pTeff2.fill(pt, wgt);
}
// Make angular efficiency plot
if (findable.isFindable(mcp, slist)) {
- double wgt = 0.;
- if (ntrk > 0) wgt = 1.;
thetaeff.fill(theta, wgt);
ctheff.fill(cth, wgt);
}
// Make d0 efficiency plot
if (findable.isFindable(mcp, slist, Ignore.NoDCACut)) {
- double wgt = 0.;
- if (ntrk > 0) wgt = 1.;
d0eff1.fill(d0, wgt);
d0eff2.fill(d0, wgt);
}
// Make z0 efficiency plot
if (findable.isFindable(mcp, slist, Ignore.NoZ0Cut)) {
- double wgt = 0.;
- if (ntrk > 0) wgt = 1.;
z0eff1.fill(z0, wgt);
z0eff2.fill(z0, wgt);
}
+ // Make the thrust axis efficiency plot
+ if (findable.isFindable(mcp, slist)) {
+ thrusteff.fill(alpha, wgt);
+ }
+
+ //
// Select charged MC particles
if (mcp.getCharge() == 0) continue;
@@ -288,6 +282,7 @@
aida.histogram1D("cos(theta) for non-final state particles", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for non-final state particles", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for non-final state particles", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for non-final state particles", 90, 0., 90.).fill(alpha);
continue;
}
@@ -297,6 +292,7 @@
aida.histogram1D("cos(theta) for base MC selection", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for base MC selection", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for base MC selection", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for base MC selection", 90, 0., 90.).fill(alpha);
// Make plots for findable tracks
if (findable.isFindable(mcp, slist)) {
@@ -305,6 +301,17 @@
aida.histogram1D("cos(theta) for findable tracks", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for findable tracks", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for findable tracks", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for findable tracks", 90, 0., 90.).fill(alpha);
+
+ // Make plots for findable tracks that are not found
+ if (trktomc.allTo(mcp).size() == 0) {
+ aida.histogram1D("Hits for unfound tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for unfound tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for unfound tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for unfound tracks", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for unfound tracks", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for unfound tracks", 90, 0., 90.).fill(alpha);
+ }
continue;
}
@@ -319,6 +326,7 @@
aida.histogram1D("cos(theta) for z0 check failures", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for z0 check failures", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for z0 check failures", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for z0 check failures", 90, 0., 90.).fill(alpha);
continue;
}
@@ -330,6 +338,7 @@
aida.histogram1D("cos(theta) for d0 check failures", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for d0 check failures", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for d0 check failures", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for d0 check failures", 90, 0., 90.).fill(alpha);
continue;
}
@@ -341,6 +350,7 @@
aida.histogram1D("cos(theta) for confirm check failures", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for seed confirm failures", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for seed confirm failures", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for seed confirm failures", 90, 0., 90.).fill(alpha);
continue;
}
@@ -352,6 +362,7 @@
aida.histogram1D("cos(theta) for seed check failures", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for seed check failures", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for seed check failures", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for seed check failures", 90, 0., 90.).fill(alpha);
continue;
}
@@ -363,6 +374,7 @@
aida.histogram1D("cos(theta) for nhit check failures", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for nhit check failures", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for nhit check failures", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for nhit check failures", 90, 0., 90.).fill(alpha);
continue;
}
@@ -374,10 +386,39 @@
aida.histogram1D("cos(theta) for pT check failures", 100, -1., 1.).fill(cth);
aida.histogram1D("d0 for pT check failures", 100, -100., 100.).fill(d0);
aida.histogram1D("z0 for pT check failures", 100, -100., 100.).fill(z0);
+ aida.histogram1D("alpha for pT check failures", 90, 0., 90.).fill(alpha);
} else {
- System.out.println("MC Particle is not findable with all ignores set!!");
+ System.out.println("**** MC Particle is not findable with all ignores set!!");
}
}
return;
}
+
+ @Override
+ protected void endOfData() {
+
+ // Make the normalized fake plot
+ double wgt = 100. / ntrktot;
+ for (int i = 0; i < 10; i++) {
+ System.out.println(" Entries: " + fakes.binEntries(i) + " for mismatches: " + i);
+ for (int j = 0; j < fakes.binHeight(i); j++) {
+ nfakes.fill(i, wgt);
+ }
+ }
+ System.out.println("Normalization: " + nfakes.sumAllBinHeights() + " after ntrk = " + ntrktot);
+ }
+
+ private Hep3Vector Thrust(List<MCParticle> mclist) {
+
+ // Make a list of final state particle momenta
+ List<Hep3Vector> plist = new ArrayList<Hep3Vector>();
+ for (MCParticle mcp : mclist) {
+ if (mcp.getGeneratorStatus() == MCParticle.FINAL_STATE)
+ plist.add(mcp.getMomentum());
+ }
+
+ // Caclulate the thrust axis
+ es.setEvent(plist);
+ return es.thrustAxis();
+ }
}
\ No newline at end of file