lcsim/sandbox/Partridge
diff -N AnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AnalysisDriver.java 16 Oct 2008 22:13:57 -0000 1.1
@@ -0,0 +1,170 @@
+/*
+ * AnalysisDriver.java
+ *
+ * Created on February 11, 2008, 11:47 AM
+ *
+ */
+
+package org.lcsim.contrib.seedtracker.example;
+
+import hep.aida.IHistogram1D;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+
+import org.lcsim.constants.Constants;
+import org.lcsim.contrib.seedtracker.SeedTrack;
+import org.lcsim.contrib.seedtracker.SeedCandidate;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.base.BaseMCParticle;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+
+
+/**
+ *
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class AnalysisDriver extends Driver {
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogram1D h1;
+ private IHistogram1D h2;
+
+ /** Creates a new instance of AnalysisDriver */
+ public AnalysisDriver() {
+ h1 = aida.histogramFactory().createHistogram1D("pT Efficiency", "", 100, 0., 50., "type=efficiency");
+ h2 = aida.histogramFactory().createHistogram1D("theta Efficiency", "", 90, 0., 180., "type=efficiency");
+
+ }
+
+ /**
+ * Process the current event
+ * @param event EventHeader for this event
+ */
+ public void process(EventHeader event) {
+ List<Track> tracklist = event.getTracks();
+ Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
+ for (Track track : tracklist) {
+ List<TrackerHit> hitlist = track.getTrackerHits();
+ 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) 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);
+ }
+ }
+
+ List<MCParticle> mclist = event.getMCParticles();
+ for (MCParticle mcp : mclist) {
+ if (mcp.getCharge() == 0) continue;
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) continue;
+ 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;
+ 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);
+ double rstart = Math.sqrt(Math.pow(mcp.getOriginX(),2)+Math.pow(mcp.getOriginY(),2));
+ double rend = Math.sqrt(Math.pow(mcp.getEndPoint().x(),2)+Math.pow(mcp.getEndPoint().y(),2));
+ if (wgt > 0.5) {
+ aida.cloud2D("Start vs end radii for found MCParticles").fill(rstart,rend);
+ } else {
+ aida.cloud2D("Start vs end radii for unfound MCParticles").fill(rstart,rend);
+// System.out.println("Unreconstructed MC Particle - Type:"+mcp.getType().getName());
+// System.out.println("Starting point: "+mcp.getOrigin().toString());
+// System.out.println("Momentum: "+mcp.getMomentum());
+// System.out.println("Endpoint: "+mcp.getEndPoint());
+// System.out.println("# Daughters: "+mcp.getDaughters().size());
+// for (MCParticle mcp2 : mcp.getDaughters()) {
+// System.out.println("Daughter type:"+mcp2.getType().getName());
+// System.out.println("Starting point: "+mcp2.getOrigin().toString());
+// System.out.println("Momentum: "+mcp2.getMomentum());
+// System.out.println("Endpoint: "+mcp2.getEndPoint());
+// System.out.println("# Daughters: "+mcp2.getDaughters().size());
+// }
+// for (MCParticle mcp2 : mcp.getParents()) {
+// System.out.println("Parent type:"+mcp2.getType().getName());
+// System.out.println("Starting point: "+mcp2.getOrigin().toString());
+// System.out.println("Momentum: "+mcp2.getMomentum());
+// System.out.println("Endpoint: "+mcp2.getEndPoint());
+// System.out.println("# Daughters: "+mcp2.getDaughters().size());
+// }
+ }
+ 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);
+ }
+ 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());
+// }
+ }
+ }
+ }
+ return;
+ }
+}
\ No newline at end of file
lcsim/sandbox/Partridge
diff -N MyTrackerDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MyTrackerDriver.java 16 Oct 2008 22:13:57 -0000 1.1
@@ -0,0 +1,72 @@
+/*
+ * MyTrackerDriver.java
+ *
+ * Created on March 29, 2006, 4:58 PM
+ *
+ */
+
+package org.lcsim.contrib.seedtracker.example;
+
+import java.util.List;
+import org.lcsim.contrib.Pelham.Example1.HistogramAnalysisDriver;
+import org.lcsim.contrib.seedtracker.SeedStrategy;
+import org.lcsim.contrib.seedtracker.SeedTracker;
+//import org.lcsim.contrib.seedtracker.diagnostic.SeedTrackerDiagnostics;
+import org.lcsim.contrib.seedtracker.StrategyXMLUtils;
+import org.lcsim.event.EventHeader;
+import org.lcsim.fit.helicaltrack.HelicalTrackHitDriver;
+import org.lcsim.fit.helicaltrack.HelicalTrackHitDriver.HitType;
+import org.lcsim.util.Driver;
+
+/**
+ * Driver for testing the SeedTracker track finding algorithm. This driver causes
+ * the SmearMCHits and SeedTracker drivers to be invoked for each event.
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class MyTrackerDriver extends Driver
+{
+ public MyTrackerDriver()
+ {
+ // Add the hit digitization driver (this example uses the virtual segment digitization code)
+ add(new VSExampleDriver());
+
+ // Add a driver to create HelicalTrackHits from the digitized hits
+ HelicalTrackHitDriver hitdriver = new HelicalTrackHitDriver();
+ hitdriver.addCollection("NewTrackerHits",HitType.VirtualSegmentation);
+ add(hitdriver);
+
+ // Obtain a trained list of strategies from a saved resource.
+ // Because resources are loaded from inside the compiled JAR file,
+ // it is necessary to recompile when modifying a saved resource.
+ // StrategyXMLUtils.getDefaultStrategiesPrefix() will provide the correct
+ // prefix in the resource tree (i.e. "/org/lcsim/contrib/seedtracker/strategybuilder/strategies/")
+// List<SeedStrategy> stratlist = StrategyXMLUtils.getStrategyListFromResource(StrategyXMLUtils.getDefaultStrategiesPrefix()+"autogen_zpole_sid01.xml");
+
+ // For development / modifications, one may choose to load a strategy
+ // list from a local file instead of a resource. The line below shows
+ // an example of how to do this
+
+// List<SeedStrategy> stratlist = StrategyXMLUtils.getStrategyListFromFile(new File("/path/to/file/strategies.xml"));
+
+
+ List<SeedStrategy> stratlist = (new MyStrategy()).getStrategies();
+ // Add a driver to do the track finding
+ SeedTracker st = new SeedTracker(stratlist);
+
+ //Enable diagnostics
+// st.setDiagnostics(new SeedTrackerDiagnostics());
+
+ //Add the Driver
+ add(st);
+
+ // Add an example analysis driver
+// add (new HistogramAnalysisDriver());
+
+ }
+ public void process(EventHeader event)
+ {
+ super.process(event);
+ return;
+ }
+}
