lcsim/src/org/lcsim/contrib/Partridge/TrackingTest
diff -u -r1.7 -r1.8
--- AnalysisDriver.java 18 Oct 2008 01:36:01 -0000 1.7
+++ AnalysisDriver.java 10 Nov 2008 19:42:24 -0000 1.8
@@ -8,28 +8,22 @@
package org.lcsim.contrib.Partridge.TrackingTest;
import hep.aida.IHistogram1D;
-import hep.physics.vec.Hep3Vector;
import java.util.List;
-import java.util.Map;
import java.util.HashMap;
+import java.util.Map;
-import org.lcsim.constants.Constants;
-import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackHit;
import org.lcsim.event.EventHeader;
import org.lcsim.event.LCRelation;
import org.lcsim.event.Track;
import org.lcsim.event.MCParticle;
import org.lcsim.event.RelationalTable;
-import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
import org.lcsim.util.aida.AIDA;
import org.lcsim.util.Driver;
-
/**
*
* @author Richard Partridge
@@ -59,82 +53,36 @@
* @param event EventHeader for this event
*/
public void process(EventHeader event) {
- RelationalTable reltab = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
- List<LCRelation> hitrelations = event.get(LCRelation.class, "HelicalTrackHitRelations");
- for (LCRelation relation : hitrelations) {
- reltab.add(relation.getFrom(), relation.getTo());
- }
- List<Track> tracklist = event.getTracks();
- Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
+
+ // Get the list of strategies being used
+ String sfile = "autogen_ttbar_sid02_vs.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(
+ StrategyXMLUtils.getDefaultStrategiesPrefix() + sfile);
+
+
+
+ // 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");
-// System.out.println("MC Relations: "+mcrelations.size());
for (LCRelation relation : mcrelations) {
- HelicalTrackHit hit = (HelicalTrackHit) relation.getFrom();
- MCParticle mcp = (MCParticle) relation.getTo();
- if (!hit.getMCParticles().contains(mcp)) System.out.println("Relation error");
-// if (mcp.getMass() > 0.9) System.out.println(mcp.toString()+hit.getMCParticles().get(0).toString());
+ hittomc.add(relation.getFrom(), relation.getTo());
}
-// List<LCRelation> hitrelations = event.get(LCRelation.class, "HelicalTrackHitRelations");
+ FindableTrack findable = new FindableTrack(event);
+
+ // Create a map between tracks and the associated MCParticle
+ List<Track> tracklist = event.getTracks();
+ Map<Track, TrackAnalysis> tkmap = new HashMap<Track, TrackAnalysis>();
+ RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
for (Track track : tracklist) {
- Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
- if (track instanceof SeedTrack) {
- SeedCandidate seed = ((SeedTrack) track).getSeedCandidate();
- for (HelicalTrackHit hit : seed.getHits()) {
- double x = hit.getCorrectedPosition().x();
- double y = hit.getCorrectedPosition().y();
- HelicalTrackFit helix = seed.getHelix();
- if (hit.getMCParticles().size() == 0) {
- System.out.println("Hit with no MCParticles:"+hit.toString());
- continue;
- }
- MCParticle particle = hit.getMCParticles().get(0);
- Hep3Vector p = particle.getMomentum();
- double phi = Math.atan2(p.y(), p.x());
- double R = Math.sqrt(p.x()*p.x() + p.y()*p.y()) / (5. * Constants.fieldConversion);
- double RS = R;
- if (particle.getCharge() < 0.) RS = -1. * RS;
- Hep3Vector start = particle.getOrigin();
- double xc = start.x() + RS * Math.sin(phi);
- double yc = start.y() - RS * Math.cos(phi);
-// double xc = helix.xc();
-// double yc = helix.yc();
-// double R = Math.abs(helix.R());
- double RHit = Math.sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
- double drphi_ms = 0.;
- if (helix.ScatterMap().containsKey(hit)) drphi_ms = helix.ScatterMap().get(hit).drphi();
- double drphi_res = hit.drphi();
- double drphi = Math.sqrt(drphi_ms*drphi_ms + drphi_res*drphi_res);
- aida.cloud1D("r-phi residual for layer "+hit.getLayerIdentifier()).fill(RHit-R);
- aida.cloud1D("Pull for layer "+hit.getLayerIdentifier()).fill((RHit - R)/drphi);
- aida.cloud1D("Hit resolution for layer "+hit.getLayerIdentifier()).fill(drphi_res);
- aida.cloud2D("Hit MS resolution for layer "+hit.getLayerIdentifier()).fill(drphi_ms, helix.pT(5.));
- List<MCParticle> mclist = hit.getMCParticles();
- for (MCParticle mcp : mclist) {
- if (mcmap.containsKey(mcp)) {
- int entries = mcmap.get(mcp) + 1;
- mcmap.put(mcp, entries);
- } else
- mcmap.put(mcp, 1);
- }
- }
- MCParticle mcmatch = null;
- int nmatch = 0;
- for (Map.Entry<MCParticle, Integer> me : mcmap.entrySet()) {
- if (me.getValue() > nmatch) {
- nmatch = me.getValue();
- mcmatch = me.getKey();
- }
- }
- if (trkmap.containsKey(mcmatch)) System.out.println("more than one track associated with an MCParticle!!");
- else trkmap.put(mcmatch, track);
- }
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+ tkmap.put(track, tkanal);
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp != null) trktomc.add(track, tkanal.getMCParticle());
}
+
List<MCParticle> mclist = event.getMCParticles();
for (MCParticle mcp : mclist) {
- if (mcp.getCharge() == 0) continue;
- if (mcp.getOrigin().magnitude() > 10.) continue;
- if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) continue;
double px = mcp.getPX();
double py = mcp.getPY();
double pz = mcp.getPZ();
@@ -142,27 +90,19 @@
double p = Math.sqrt(pt*pt + pz*pz);
double cth = pz / p;
double theta = 180. * Math.acos(cth) / Math.PI;
- double wgt = 0.;
- if (trkmap.containsKey(mcp)) wgt = 1.;
- if (pt > 1.1) {
- aida.profile1D("Efficiency vs theta", 90, 0., 180.).fill(theta, wgt);
-// h2.fill(theta, wgt);
- aida.histogram1D("MC angle", 90, 0., 180.).fill(theta);
- }
- if (Math.abs(cth) < 0.985) {
- aida.profile1D("Efficiency vs pT", 100, 0., 50.).fill(pt, wgt);
-// h1.fill(pt,wgt);
- aida.histogram1D("MC pT", 100, 0., 50.).fill(pt);
+ int ntrk = trktomc.allTo(mcp).size();
+ if (ntrk > 1) System.out.println("2 tracks associated with a single MC Particle");
+ if (findable.isFindable(mcp, slist, FindableTrack.Ignore.NoPTCut)) {
+ System.out.println("Findable MC Particle found: P = "+mcp.getMomentum().toString());
+ double wgt = 0.;
+ if (ntrk > 0) wgt = 1.;
+ pTeff1.fill(pt, wgt);
+ pTeff2.fill(pt, wgt);
}
- if (!trkmap.containsKey(mcp)) {
- List<HelicalTrackHit> hits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
-// if (hits.size() > 6) {
-// System.out.println("Failed to find track. Found "+hits.size()+" hits @ theta = "
-// +180.*Math.acos(mcp.getMomentum().z() / mcp.getMomentum().magnitude())/Math.PI);
-// for (HelicalTrackHit hit : hits) {
-// System.out.println("Hit in "+hit.getLayerIdentifier()+" at x= "+hit.x()+" y= "+hit.y()+" z= "+hit.z());
-// }
-// }
+ if (findable.isFindable(mcp, slist)) {
+ double wgt = 0.;
+ if (ntrk > 0) wgt = 1.;
+ thetaeff.fill(theta, wgt);
}
}
return;
lcsim/src/org/lcsim/contrib/Partridge/TrackingTest
diff -u -r1.1 -r1.2
--- FindableTrack.java 4 Nov 2008 21:18:54 -0000 1.1
+++ FindableTrack.java 10 Nov 2008 19:42:24 -0000 1.2
@@ -9,15 +9,23 @@
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
+
import java.util.ArrayList;
-import java.util.HashMap;
import java.util.List;
-import java.util.Map;
-import org.lcsim.event.EventHeader;
+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.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
+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;
@@ -30,29 +38,39 @@
public class FindableTrack {
public enum Ignore {NoPTCut, NoDCACut, NoZ0Cut, NoSeedCheck, NoConfirmCheck, NoMinHitCut};
private double _bfield;
- private List<HelicalTrackHit> _hitlist;
- private Map<MCParticle, List<HelicalTrackHit>> _hitmap;
+ private RelationalTable _hittomc;
+ private HitIdentifier _ID;
/** Creates a new instance of FindableTrack */
- public FindableTrack(EventHeader event, String hitcolname) {
+ public FindableTrack(EventHeader event) {
// Get the magnetic field
Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
_bfield = event.getDetector().getFieldMap().getField(IP).z();
- // Get the HelicalTrackHits
- _hitlist = (List<HelicalTrackHit>) event.get(HelicalTrackHit.class, hitcolname);
+ // Instantiate the hit identifier class
+ _ID = new HitIdentifier();
- // Make a map between the MCParticles and the hits on the MCParticle
- _hitmap = new HashMap<MCParticle, List<HelicalTrackHit>>();
- for (HelicalTrackHit hit : _hitlist) {
- for (MCParticle mcp : hit.getMCParticles()) {
- if (!_hitmap.containsKey(mcp)) _hitmap.put(mcp, new ArrayList<HelicalTrackHit>());
- _hitmap.get(mcp).add(hit);
+ // 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) {
+ _hittomc.add(simhit, simhit.getMCParticle());
}
}
}
+ 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>());
}
@@ -61,7 +79,7 @@
// 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);
@@ -73,22 +91,22 @@
// Check the MC Particle's pT
if (!CheckPT(helix, ignores, strat)) continue;
-
+
// Check the MC Particle's DCA
if (!CheckDCA(helix, ignores, strat)) continue;
-
+
// Check the MC Particle's Z0
if (!CheckZ0(helix, ignores, strat)) continue;
-
+
// Check that we have hits on the seed layers
if (!CheckSeed(mcp, ignores, strat)) continue;
-
+
// Check that we have the required confirmation hits
if (!CheckConfirm(mcp, ignores, strat)) continue;
-
+
// Check for the minimum number of hits
if (!CheckMinHits(mcp, ignores, strat)) continue;
-
+
// Passed all the checks - track is findable
findable = true;
break;
@@ -129,14 +147,14 @@
return HitCount(mcp, strat.getLayers(SeedType.Seed)) == 3;
}
- private boolean CheckConfirm(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
+ 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
@@ -149,17 +167,20 @@
private int HitCount(MCParticle mcp, List<SeedLayer> lyrlist) {
// Get the list of hits associated with the MCParticle
- List<HelicalTrackHit> hitlist = _hitmap.get(mcp);
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
// Count the number of hits in layers specified by the strategy
int hitcount = 0;
- for (HelicalTrackHit hit : hitlist) {
+ 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) {
- if (!lyr.getDetName().equals(hit.Detector())) continue;
- if (lyr.getLayer() != hit.Layer()) continue;
- if (!lyr.getBarrelEndcapFlag().equals(hit.BarrelEndcapFlag())) continue;
+ if (!lyr.getDetName().equals(_ID.getName(de))) continue;
+ if (lyr.getLayer() != _ID.getLayer(de)) continue;
+ if (!lyr.getBarrelEndcapFlag().equals(_ID.getBarrelEndcapFlag(de))) continue;
hitcount++;
}
}