lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/TrackingTest
diff -u -r1.1.1.1 -r1.2
--- AnalysisDriver.java 10 Dec 2008 22:03:06 -0000 1.1.1.1
+++ AnalysisDriver.java 30 Jan 2009 15:46:16 -0000 1.2
@@ -4,11 +4,11 @@
* Created on February 11, 2008, 11:47 AM
*
*/
-
package org.lcsim.contrib.Partridge.TrackingTest;
import hep.aida.IHistogram1D;
+import hep.physics.vec.Hep3Vector;
import java.util.List;
import java.util.HashMap;
import java.util.Map;
@@ -19,6 +19,8 @@
import org.lcsim.event.MCParticle;
import org.lcsim.event.RelationalTable;
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.aida.AIDA;
@@ -30,6 +32,7 @@
* @version 1.0
*/
public class AnalysisDriver extends Driver {
+
private AIDA aida = AIDA.defaultInstance();
private IHistogram1D pTeff1;
private IHistogram1D pTeff2;
@@ -37,63 +40,79 @@
private IHistogram1D ctheff;
private IHistogram1D d0eff;
private IHistogram1D z0eff;
-
+
/** Creates a new instance of AnalysisDriver */
public AnalysisDriver() {
- pTeff1 = aida.histogramFactory().createHistogram1D("Efficiency vs pT", "", 100, 0., 5., "type=efficiency");
+ pTeff1 = aida.histogramFactory().createHistogram1D("Efficiency vs pT", "", 40, 0., 5., "type=efficiency");
pTeff2 = aida.histogramFactory().createHistogram1D("Efficiency vs pT full", "", 100, 0., 50., "type=efficiency");
- thetaeff = aida.histogramFactory().createHistogram1D("Efficiency vs theta", "", 360, 0., 180., "type=efficiency");
+ thetaeff = aida.histogramFactory().createHistogram1D("Efficiency vs theta", "", 72, 0., 180., "type=efficiency");
ctheff = aida.histogramFactory().createHistogram1D("Efficiency vs cos(theta)", "", 200, -1., 1., "type=efficiency");
d0eff = aida.histogramFactory().createHistogram1D("Efficiency vs d0", "", 360, -20., 20., "type=efficiency");
z0eff = aida.histogramFactory().createHistogram1D("Efficiency vs z0", "", 360, -20., 20., "type=efficiency");
}
-
+
/**
* Process the current event
* @param event EventHeader for this event
*/
public void process(EventHeader event) {
-
+
// Get the list of strategies being used
- String sfile = "autogen_ttbar_sid02_vs.xml";
+ String sfile = "AtlasU01.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");
for (LCRelation relation : mcrelations) {
hittomc.add(relation.getFrom(), relation.getTo());
}
-
+
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) {
TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+ aida.cloud1D("Purity for " + tkanal.getNHits() + " hits").fill(tkanal.getPurity());
tkmap.put(track, tkanal);
MCParticle mcp = tkanal.getMCParticle();
- if (mcp != null) trktomc.add(track, tkanal.getMCParticle());
+ if (mcp != null) {
+ trktomc.add(track, tkanal.getMCParticle());
+ 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);
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, 2.0);
+ double d0tk = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double d0mc = helix.getDCA();
+ aida.cloud2D("Pt trk vs MC with "+tkanal.getNBadHits()+" bad hits").fill(ptmc,pttk);
+ aida.cloud2D("d0 trk vs MC with "+tkanal.getNBadHits()+" bad hits").fill(d0mc,d0tk);
+ }
}
-
+
List<MCParticle> mclist = event.getMCParticles();
for (MCParticle mcp : mclist) {
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 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;
int ntrk = trktomc.allTo(mcp).size();
- if (ntrk > 1) System.out.println("2 tracks associated with a single MC Particle");
+ 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());
+// System.out.println("Findable MC Particle found: P = "+mcp.getMomentum().toString());
double wgt = 0.;
if (ntrk > 0) wgt = 1.;
pTeff1.fill(pt, wgt);