lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/TrackingTest
diff -N LOIEffFake.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LOIEffFake.java 25 Feb 2009 17:51:32 -0000 1.1
@@ -0,0 +1,383 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.Partridge.TrackingTest;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import java.util.Set;
+import org.lcsim.contrib.Partridge.TrackingTest.FindableTrack.Ignore;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author partridge
+ */
+public class LOIEffFake extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogram1D pTeff1;
+ private IHistogram1D pTeff2;
+ private IHistogram1D thetaeff;
+ private IHistogram1D ctheff;
+ private IHistogram1D d0eff1;
+ private IHistogram1D d0eff2;
+ private IHistogram1D z0eff1;
+ private IHistogram1D z0eff2;
+ private IHistogram1D fakes;
+ private IHistogram1D nfakes;
+ int ntrk = 0;
+ int nevt = 0;
+
+ public LOIEffFake() {
+
+ // Define the efficiency histograms
+ IHistogramFactory hf = aida.histogramFactory();
+ pTeff1 = hf.createHistogram1D("Efficiency vs pT", "", 100, 0., 5., "type=efficiency");
+ pTeff2 = hf.createHistogram1D("Efficiency vs pT full", "", 100, 0., 50., "type=efficiency");
+ thetaeff = hf.createHistogram1D("Efficiency vs theta", "", 72, 0., 180., "type=efficiency");
+ ctheff = hf.createHistogram1D("Efficiency vs cos(theta)", "", 200, -1., 1., "type=efficiency");
+ d0eff1 = hf.createHistogram1D("Efficiency vs d0", "", 50, -5., 5., "type=efficiency");
+ 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");
+ fakes = hf.createHistogram1D("Number of mis-matched hits (unnormalized)", "", 10, 0., 10.);
+ nfakes = hf.createHistogram1D("Number of mis-matched hits (normalized)", "", 10, 0., 10.);
+ }
+
+ @Override
+ public void process(EventHeader event) {
+
+ // Increment the event counter
+ nevt++;
+
+ // Get the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
+ double bfield = event.getDetector().getFieldMap().getField(IP).z();
+
+ // Get the list of strategies being used
+ String sfile = "autogen_ttbar_sid02_vs.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(
+ StrategyXMLUtils.getDefaultStrategiesPrefix() + sfile);
+
+ // Find the minimum pT among the strategies
+ double ptCut = 9999.;
+ for (SeedStrategy s : slist) {
+ if (s.getMinPT() < ptCut) ptCut = s.getMinPT();
+ }
+
+ // 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");
+ for (LCRelation relation : mcrelations) {
+ hittomc.add(relation.getFrom(), relation.getTo());
+ }
+
+ // Instantiate the class that determines if a track is "findable"
+ FindableTrack findable = new FindableTrack(event);
+
+ // Create a map between tracks and the associated MCParticle
+ List<Track> tracklist = event.getTracks();
+ RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ // Analyze the tracks in the event
+ for (Track track : tracklist) {
+
+ // Calculate the track pT and cos(theta)
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double pt = Math.sqrt(px * px + py * py);
+ double cth = pz / Math.sqrt(pt * pt + pz * pz);
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+
+ // Analyze the hits on the track
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+
+ // Calculate purity and make appropriate plots
+ int nbad = tkanal.getNBadHits();
+ int nhits = tkanal.getNHits();
+ 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);
+
+ // Make plots for fake, non-fake, and all tracks
+ if (purity < 0.5) {
+ aida.histogram1D("Hits for fake tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for fake tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for fake tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for fake tracks", 50, -10., 10.).fill(d0);
+ aida.histogram1D("z0 for fake tracks", 50, -10., 10.).fill(z0);
+ } else {
+ aida.histogram1D("Hits for non-fake tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for non-fake tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for non-fake tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for non-fake tracks", 50, -10., 10.).fill(d0);
+ aida.histogram1D("z0 for non-fake tracks", 50, -10., 10.).fill(z0);
+ }
+ aida.histogram1D("Hits for all tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for all tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for all tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for all tracks", 50, -10., 10.).fill(d0);
+ aida.histogram1D("z0 for all tracks", 50, -10., 10.).fill(z0);
+
+ // Now analyze MC Particles on this track
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp != null) {
+
+ // Create a map between the tracks found and the assigned MC particle
+ trktomc.add(track, tkanal.getMCParticle());
+
+ // Calculate the MC momentum and polar angle
+ Hep3Vector pmc = mcp.getMomentum();
+ double pxmc = pmc.x();
+ double pymc = pmc.y();
+ double ptmc = Math.sqrt(pxmc * pxmc + pymc * pymc);
+ double pxtk = track.getPX();
+ double pytk = track.getPY();
+ double pttk = Math.sqrt(pxtk * pxtk + pytk * pytk);
+
+ // Calculate the helix parameters for this MC particle and pulls in pT, d0
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+ double d0tk = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double d0mc = helix.getDCA();
+ double d0err = Math.sqrt(track.getErrorMatrix().diagonal(HelicalTrackFit.dcaIndex));
+ double curv = track.getTrackParameter(HelicalTrackFit.curvatureIndex);
+ double curverr = Math.sqrt(track.getErrorMatrix().diagonal(HelicalTrackFit.curvatureIndex));
+ double pterr = pttk * curverr / curv;
+ double d0pull = (d0tk - d0mc) / d0err;
+ double ptpull = (pttk - ptmc) / pterr;
+
+ // Plot the pt and d0 pulls for various purity intervals
+ if (nbad == 0) {
+ aida.histogram2D("pT MC vs pT Reco for 0 Bad Hits",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D("d0 MC vs d0 Reco for 0 Bad Hits",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.histogram1D("pT Pull for 0 Bad Hits", 100, -10., 10.).fill(ptpull);
+ aida.histogram1D("d0 pull for 0 Bad Hits", 100, -10., 10.).fill(d0pull);
+ } else if (purity > 0.5) {
+ aida.histogram2D("pT MC vs pT Reco for 0.5 < purity < 1",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D("d0 MC vs d0 Reco for 0.5 < purity < 1",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.histogram1D("pT Pull for 0.5 < purity < 1", 100, -10., 10.).fill(ptpull);
+ aida.histogram1D("d0 pull for 0.5 < purity < 1", 100, -10., 10.).fill(d0pull);
+ } else if (purity < 0.5) {
+ aida.histogram2D("pT MC vs pT Reco for purity <= 0.5",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D("d0 MC vs d0 Reco for purity <= 0.5",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.histogram1D("pT Pull for purity <= 0.5", 100, -10., 10.).fill(ptpull);
+ aida.histogram1D("d0 pull for purity <= 0.5", 100, -10., 10.).fill(d0pull);
+ }
+ }
+ }
+
+ // 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);
+ }
+
+ // 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();
+ double pt = Math.sqrt(px * px + py * py);
+ double p = Math.sqrt(pt * pt + pz * pz);
+ double cth = pz / p;
+ double theta = 180. * Math.acos(cth) / Math.PI;
+
+ // Find the number of layers hit by this mc particle
+ int nhits = findable.LayersHit(mcp);
+
+ // Calculate the helix parameters for this MC particle
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+ 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());
+// }
+// }
+
+ // 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);
+ }
+
+ // Select charged MC particles
+ if (mcp.getCharge() == 0) continue;
+
+ // Select mcp that fail the final state requirement
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) {
+ aida.histogram1D("Hits for non-final state particles", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for non-final state particles", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // Make plots for the base sample
+ aida.histogram1D("Hits for base MC selection", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for base MC selection", 100, 0., 10.).fill(pt);
+ 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);
+
+ // Make plots for findable tracks
+ if (findable.isFindable(mcp, slist)) {
+ aida.histogram1D("Hits for findable tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for findable tracks", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // Create the running list of conditions to ignore
+ List<Ignore> ignores = new ArrayList<Ignore>();
+
+ // select mc particles that fail on the z0 cut
+ ignores.add(Ignore.NoZ0Cut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for z0 check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for z0 check failures", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // Select mc particles that fail on the d0 cut
+ ignores.add(Ignore.NoDCACut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for d0 check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for d0 check failures", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // select mc particles that fail the confirm check
+ ignores.add(Ignore.NoConfirmCheck);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for confirm check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for confir check failures", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // select mc particles that fail on the seed check
+ ignores.add(Ignore.NoSeedCheck);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for seed check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for seed check failures", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // Select mc particles that fail the number of hit cut
+ ignores.add(Ignore.NoMinHitCut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for nhit check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for nhit check failures", 100, 0., 10.).fill(pt);
+ 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);
+ continue;
+ }
+
+ // Select mc particles that fail on the pT cut
+ ignores.add(Ignore.NoPTCut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for pT check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for pT check failures", 100, 0., 10.).fill(pt);
+ 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);
+ } else {
+ System.out.println("MC Particle is not findable with all ignores set!!");
+ }
+ }
+ return;
+ }
+}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/TrackingTest
diff -u -r1.1.1.1 -r1.2
--- FindableTrack.java 10 Dec 2008 22:03:06 -0000 1.1.1.1
+++ FindableTrack.java 25 Feb 2009 17:51:32 -0000 1.2
@@ -4,20 +4,17 @@
* Created on October 24, 2008, 9:50 PM
*
*/
-
package org.lcsim.contrib.Partridge.TrackingTest;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import java.util.ArrayList;
+import java.util.HashSet;
import java.util.List;
import java.util.Set;
-import org.lcsim.detector.DetectorElementStore;
import org.lcsim.detector.IDetectorElement;
-import org.lcsim.detector.IDetectorElementContainer;
-import org.lcsim.detector.identifier.IIdentifierHelper;
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.RelationalTable;
@@ -25,7 +22,6 @@
import org.lcsim.event.base.BaseRelationalTable;
import org.lcsim.fit.helicaltrack.HelixParamCalculator;
import org.lcsim.fit.helicaltrack.HitIdentifier;
-import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
import org.lcsim.recon.tracking.seedtracker.SeedLayer;
import org.lcsim.recon.tracking.seedtracker.SeedLayer.SeedType;
import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
@@ -36,27 +32,31 @@
* @version 1.0
*/
public class FindableTrack {
- public enum Ignore {NoPTCut, NoDCACut, NoZ0Cut, NoSeedCheck, NoConfirmCheck, NoMinHitCut};
+
+ public enum Ignore {
+
+ NoPTCut, NoDCACut, NoZ0Cut, NoSeedCheck, NoConfirmCheck, NoMinHitCut
+ };
private double _bfield;
private RelationalTable _hittomc;
private HitIdentifier _ID;
-
+
/** Creates a new instance of FindableTrack */
public FindableTrack(EventHeader event) {
-
+
// Get the magnetic field
Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
_bfield = event.getDetector().getFieldMap().getField(IP).z();
-
+
// Instantiate the hit identifier class
_ID = new HitIdentifier();
-
+
// Create a relational table that maps SimTrackerHits to MCParticles
_hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
-
+
// Get the collections of SimTrackerHits
List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
-
+
// Loop over the SimTrackerHits and fill in the relational table
for (List<SimTrackerHit> simlist : simcols) {
for (SimTrackerHit simhit : simlist) {
@@ -64,34 +64,34 @@
}
}
}
-
+
public boolean isFindable(MCParticle mcp, List<SeedStrategy> slist, Ignore ignore) {
List<Ignore> ignores = new ArrayList<Ignore>();
ignores.add(ignore);
return isFindable(mcp, slist, ignores);
}
-
+
public boolean isFindable(MCParticle mcp, List<SeedStrategy> slist) {
return isFindable(mcp, slist, new ArrayList<Ignore>());
}
-
+
public boolean isFindable(MCParticle mcp, List<SeedStrategy> slist, List<Ignore> ignores) {
-
+
// We can't find neutral particles'
if (mcp.getCharge() == 0) return false;
// Find the helix parameters in the L3 convention used by org.lcsim
HelixParamCalculator helix = new HelixParamCalculator(mcp, _bfield);
-
+
// We haven't yet determined the track is findable
boolean findable = false;
-
+
// Loop over strategies and check if the track is findable
for (SeedStrategy strat : slist) {
-
+
// Check the MC Particle's pT
if (!CheckPT(helix, ignores, strat)) continue;
-
+
// Check the MC Particle's DCA
if (!CheckDCA(helix, ignores, strat)) continue;
@@ -111,80 +111,100 @@
findable = true;
break;
}
-
+
return findable;
}
-
+
+ public int LayersHit(MCParticle mcp) {
+
+ // Get the list of hits associated with the MCParticle
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+
+ // Create a set of the identifiers for the hit layers
+ Set<String> idset = new HashSet<String>();
+
+ // Create the set of identifiers
+ for (SimTrackerHit simhit : hitlist) {
+ String identifier = _ID.Identifier(simhit.getDetectorElement());
+ if (!idset.contains(identifier)) idset.add(identifier);
+ }
+
+ return idset.size();
+ }
+
private boolean CheckPT(HelixParamCalculator helix, List<Ignore> ignores, SeedStrategy strat) {
-
+
// First see if we are skipping this check
if (ignores.contains(Ignore.NoPTCut)) return true;
-
+
return helix.getMCTransverseMomentum() >= strat.getMinPT();
}
-
+
private boolean CheckDCA(HelixParamCalculator helix, List<Ignore> ignores, SeedStrategy strat) {
-
+
// First see if we are skipping this check
if (ignores.contains(Ignore.NoDCACut)) return true;
-
+
return Math.abs(helix.getDCA()) <= strat.getMaxDCA();
}
-
+
private boolean CheckZ0(HelixParamCalculator helix, List<Ignore> ignores, SeedStrategy strat) {
-
+
// First see if we are skipping this check
if (ignores.contains(Ignore.NoZ0Cut)) return true;
-
+
return Math.abs(helix.getZ0()) <= strat.getMaxZ0();
}
-
+
private boolean CheckSeed(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
-
+
// First see if we are skipping this check
if (ignores.contains(Ignore.NoSeedCheck)) return true;
-
+
return HitCount(mcp, strat.getLayers(SeedType.Seed)) == 3;
}
-
+
private boolean CheckConfirm(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
-
+
// First see if we are skipping this check
if (ignores.contains(Ignore.NoConfirmCheck)) return true;
-
+
return HitCount(mcp, strat.getLayers(SeedType.Confirm)) >= strat.getMinConfirm();
}
-
+
private boolean CheckMinHits(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
-
+
// First see if we are skipping this check
if (ignores.contains(Ignore.NoMinHitCut)) return true;
-
+
return HitCount(mcp, strat.getLayerList()) >= strat.getMinHits();
}
-
-
+
private int HitCount(MCParticle mcp, List<SeedLayer> lyrlist) {
-
+
// Get the list of hits associated with the MCParticle
Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
-
- // Count the number of hits in layers specified by the strategy
+
+ // Count the number of layers with hits in them
int hitcount = 0;
- for (SimTrackerHit simhit : hitlist) {
-
- // Get the detector element for this hit
- IDetectorElement de = simhit.getDetectorElement();
-
- // See if this hit is on one of the specified layers
- for (SeedLayer lyr : lyrlist) {
+ for (SeedLayer lyr : lyrlist) {
+
+ // Loop over the hits for this MCParticle
+ for (SimTrackerHit simhit : hitlist) {
+
+ // Get the detector element for this hit
+ IDetectorElement de = simhit.getDetectorElement();
+
+ // See if this hit is on the layer we are checking
if (!lyr.getDetName().equals(_ID.getName(de))) continue;
if (lyr.getLayer() != _ID.getLayer(de)) continue;
- if (!lyr.getBarrelEndcapFlag().equals(_ID.getBarrelEndcapFlag(de))) continue;
+ if (!lyr.getBarrelEndcapFlag().equals(_ID.getBarrelEndcapFlag(de)))
+ continue;
hitcount++;
+ break;
}
}
-
+
return hitcount;
}
}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/TrackingTest
diff -u -r1.2 -r1.3
--- ResolutionAnalysis.java 20 Feb 2009 22:50:15 -0000 1.2
+++ ResolutionAnalysis.java 25 Feb 2009 17:51:32 -0000 1.3
@@ -18,6 +18,8 @@
import org.lcsim.event.RelationalTable;
import org.lcsim.event.Track;
import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -33,8 +35,13 @@
public ResolutionAnalysis() {
}
+ @Override
public void process(EventHeader event) {
-
+
+ // Find the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
+ double bfield = event.getDetector().getFieldMap().getField(IP).z();
+
// 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");
@@ -45,15 +52,7 @@
nevt++;
List<Track> tracklist = event.getTracks();
List<MCParticle> mclist = event.getMCParticles();
- for (MCParticle mcp : mclist) {
- Hep3Vector p = mcp.getMomentum();
- double pt = Math.sqrt(p.x()*p.x()+p.y()*p.y());
- if (pt<0.2) continue;
- double angle = 180. * Math.acos(p.z()/p.magnitude())/Math.PI;
- if (angle < 8.5 || angle > 171.5) continue;
- if (tracklist.size() != 1) System.out.println("Event: "+nevt+" pT: "+pt+" theta: "+angle);
- }
-// if (tracklist.size() != 1) System.out.println("Found "+tracklist.size()+" Tracks!");
+
for (Track track : tracklist) {
TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
MCParticle mcp = tkanal.getMCParticle();
@@ -61,17 +60,30 @@
double ptrk = (new BasicHep3Vector(track.getMomentum())).magnitude();
double pmc = mcp.getMomentum().magnitude();
double pdif = ptrk - pmc;
+ double d0trk = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0trk = track.getTrackParameter(HelicalTrackFit.z0Index);
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+ double d0mc = helix.getDCA();
+ double d0dif = d0trk - d0mc;
+ double z0mc = helix.getDCA();
+ double z0dif = z0trk - z0mc;
double theta = Math.acos(mcp.getMomentum().z()/pmc);
if (mcp.getCharge() > 0.) {
- aida.cloud1D("Momentum error for positive tracks").fill(pdif);
+// aida.cloud1D("Momentum error for positive tracks").fill(pdif);
+// aida.cloud1D("d0 error for positive tracks").fill(d0dif);
+// aida.histogram1D("z0 error for positive tracks",100,-.1,.1).fill(z0dif);
} else {
- aida.cloud1D("Momentum error for negative tracks").fill(pdif);
+// aida.cloud1D("Momentum error for negative tracks").fill(pdif);
+// aida.histogram1D("d0 error for negative tracks",100,-.1,.1).fill(d0dif);
+// aida.histogram1D("z0 error for negative tracks",100,-.1,.1).fill(z0dif);
}
- aida.cloud1D("Momentum error for all tracks").fill(pdif);
- aida.cloud1D("Track momentum for all tracks").fill(ptrk);
- aida.cloud1D("MC momentum for all tracks").fill(pmc);
- aida.cloud1D("Momentum resolution for all tracks").fill(pdif/pmc);
- aida.cloud1D("Polar angle").fill(theta);
+ aida.cloud1D("d0 error for all tracks").fill(d0dif);
+// aida.histogram1D("z0 error for all tracks", 100, -.1, .1).fill(z0dif);
+// aida.cloud1D("Momentum error for all tracks").fill(pdif);
+// aida.cloud1D("Track momentum for all tracks").fill(ptrk);
+// aida.cloud1D("MC momentum for all tracks").fill(pmc);
+// aida.cloud1D("Momentum resolution for all tracks").fill(pdif/pmc);
+// aida.cloud1D("Polar angle").fill(theta);
}
}
