hps-java/sandbox
diff -N HPSHelicalTrackHitDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HPSHelicalTrackHitDriver.java 22 Nov 2012 02:52:04 -0000 1.1
@@ -0,0 +1,507 @@
+/*
+ * HPSHelicalTrackHitDriverOLD.java
+ *
+ * Created on January 4, 2010
+ *
+ */
+package org.lcsim.hps.recon.tracking;
+
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.*;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.tracker.silicon.SiTrackerModule;
+import org.lcsim.digisim.MyLCRelation;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.fit.helicaltrack.*;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
+import org.lcsim.hps.event.HPSTransformations;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType.CoordinateSystem;
+import org.lcsim.recon.tracking.vsegment.hit.DigiTrackerHit;
+import org.lcsim.recon.tracking.vsegment.hit.TrackerCluster;
+import org.lcsim.util.Driver;
+
+/**
+ Create the appropriate HelicalTrackHits for the specified TrackerHit
+ collections. The resulting HelicalTrackHits encapsulate the information needed
+ to perform a helical track hit for either a segmented strip detector, a pixel
+ detector, or cross hits from a stereo detector.
+
+ This driver was developed for the Dark Forces silicon tracker, where digitized
+ hits are in either axial or stereo layers. The pairs of layers to be used in
+ forming stereo hits can be specified using the setStereoPairs method.
+
+ The list of hit collections to be converted must be specified before the process
+ method is executed.
+
+ @author Richard Partridge
+ @version 1.0
+ */
+public class HPSHelicalTrackHitDriverOLD extends Driver {
+// private StereoHitMaker _crosser = new StereoHitMaker(2., 10.);
+
+ private boolean _debug = false;
+ private StereoHitMaker _crosser = new StereoHitMaker(0.01, 10.01);
+ private HitIdentifier _ID = new HitIdentifier();
+ private String _outname = "HelicalTrackHits";
+ private String _hitrelname = "HelicalTrackHitRelations";
+ private String _mcrelname = "HelicalTrackMCRelations";
+ private List<String> _colnames = new ArrayList<String>();
+ private Map<String, String> _stereomap = new HashMap<String, String>();
+ private Hep3Vector _orgloc = new BasicHep3Vector(0., 0., 0.);
+ private boolean _doTranformToTracking = true;
+ private boolean addAllStripsToEventForDebug = false;
+ private HPSTransformations _detToTrk;
+ String subdetectorName = "Tracker";
+
+ /**
+ Creates a new instance of HelicalTrackHitDriver
+ */
+ public HPSHelicalTrackHitDriverOLD() {
+// setDetToTrkMatrix();
+ _detToTrk = new HPSTransformations();
+ _colnames.add("StripClusterer_SiTrackerHitStrip1D");
+ }
+
+ public void setSubdetectorName(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ @Override
+ public void detectorChanged(Detector detector) {
+ // Setup default stereo pairings, which should work for even number of
+ // modules.
+ List<SiTrackerModule> modules = detector.getSubdetector(subdetectorName).getDetectorElement().findDescendants(SiTrackerModule.class);
+ if (modules.isEmpty()) {
+ throw new RuntimeException("No SiTrackerModules found in detector.");
+ }
+ if (_debug)
+ System.out.println("Number of layers = " + detector.getSubdetector(subdetectorName).getLayering().getLayers().getNumberOfLayers());
+ int nlayers = detector.getSubdetector(subdetectorName).getLayering().getLayers().getNumberOfLayers();
+ if (nlayers % 2 != 0) {
+ throw new RuntimeException("Don't know how to do stereo pairing for odd number of modules.");
+ }
+ List<int[]> pairs = new ArrayList<int[]>();
+ for (int i = 1; i <= (nlayers) - 1; i += 2) {
+ int[] pair = {i, i + 1};
+ if (_debug)
+ System.out.println("Adding stereo pair: " + pair[0] + ", " + pair[1]);
+ pairs.add(pair);
+ }
+ for (int[] pair : pairs) {
+ setStereoPair(subdetectorName, pair[0], pair[1]);
+ }
+ }
+
+ /**
+ Create the HelicalTrackHits for the specified hit collections.
+
+ @param event EventHeader of the event to be processed
+ */
+ @Override
+ public void process(EventHeader event) {
+ super.process(event);
+
+ // Initialize the list of HelicalTrackHits
+ List<HelicalTrackHit> helhits = new ArrayList<HelicalTrackHit>();
+
+ // Create a List of LCRelations to relate HelicalTrackHits to the original hits
+ List<LCRelation> hitrelations = new ArrayList<LCRelation>();
+
+ List<HelicalTrackCross> stereoCrosses = new ArrayList<HelicalTrackCross>();
+
+ List<HelicalTrackStrip> allstrips = new ArrayList<HelicalTrackStrip>();
+
+ for (String _colname : _colnames) {
+
+ if (!event.hasCollection(SiTrackerHit.class, _colname)) {
+ continue;
+ }
+
+ // Get the list of SiTrackerHits for this collection
+ List<SiTrackerHit> hitlist = (List<SiTrackerHit>) event.get(_colname);
+
+ // Create collections for strip hits by layer and hit cross references
+ Map<String, List<HelicalTrackStrip>> striplistmap = new HashMap<String, List<HelicalTrackStrip>>();
+ Map<HelicalTrackStrip, SiTrackerHitStrip1D> stripmap = new HashMap<HelicalTrackStrip, SiTrackerHitStrip1D>();
+ if (_debug) {
+ System.out.println("HPSHelicalTrackHitDriverOLD: # of hits found = " + hitlist.size());
+ }
+
+ // Loop over the SiTrackerHits in this collection
+ for (SiTrackerHit hit : hitlist) {
+
+ // Check if we have a 1D strip hit
+ if (hit instanceof SiTrackerHitStrip1D) {
+
+ // Cast the hit as a 1D strip hit and find the identifier for the detector/layer combo
+ SiTrackerHitStrip1D h = (SiTrackerHitStrip1D) hit;
+ IDetectorElement de = h.getSensor();
+ String id = makeID(_ID.getName(de), _ID.getLayer(de));
+
+ // Check if hit is on a stereo pair layer
+ // Why wouldn't it be on a pair?
+ if (_stereomap.containsKey(id) || _stereomap.containsValue(id)) {
+
+ // Create a HelicalTrackStrip for this hit
+ HelicalTrackStrip strip = makeDigiStrip(h);
+
+ // Get the list of strips for this layer - create a new list if one doesn't already exist
+ List<HelicalTrackStrip> lyrhits = striplistmap.get(id);
+ if (lyrhits == null) {
+ lyrhits = new ArrayList<HelicalTrackStrip>();
+ striplistmap.put(id, lyrhits);
+ }
+
+ // Add the strip to the list of strips on this sensor
+ lyrhits.add(strip);
+
+ // Map a reference back to the hit needed to create the stereo hit LC relations
+ stripmap.put(strip, h);
+
+ if( addAllStripsToEventForDebug ) allstrips.add(strip);
+
+ } else {
+
+ // If not a stereo hit, create an axial strip hit and the corresponding LC Relation
+ HelicalTrackHit dah = makeDigiAxialHit(h);
+ helhits.add(dah);
+ hitrelations.add(new MyLCRelation(dah, hit));
+ }
+
+ } else {
+
+ // If not a 1D strip hit, make a pixel hit
+ HelicalTrackHit hit3d = makeDigi3DHit(hit);
+ helhits.add(hit3d);
+ hitrelations.add(new MyLCRelation(hit3d, hit));
+ }
+ }
+
+ if (_debug) {
+ for (Map.Entry<String, List<HelicalTrackStrip>> entry : striplistmap.entrySet()) {
+ String ident = entry.getKey();
+ List<HelicalTrackStrip> hts = entry.getValue();
+ System.out.println(this.getClass().getSimpleName() + ": Total number of HelicalTrackStrips: " + ident + ": " + hts.size());
+ for(HelicalTrackStrip strip : hts) {
+ System.out.println(this.getClass().getSimpleName() + ": strip at origin " + strip.origin().toString() + " umeas " + strip.umeas());
+ }
+ }
+ }
+
+ // Now create the stereo hits
+ // Create a list of stereo hits
+ List<HelicalTrackCross> stereohits = new ArrayList<HelicalTrackCross>();
+
+ if (_debug) System.out.println(this.getClass().getSimpleName() + ": Create stereo hits from " + striplistmap.size() + " strips (map size)");
+
+ // Loop over the stereo layer pairs
+ for (String id1 : _stereomap.keySet()) {
+
+ // Get the second layer
+ String id2 = _stereomap.get(id1);
+
+ if (_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": Form stereo hits from " + id1 + " ( " + striplistmap.get(id1) + ") and " + id2 + " (" + striplistmap.get(id2) );
+ }
+
+ // Form the stereo hits and add them to our hit list
+ stereohits.addAll(_crosser.MakeHits(striplistmap.get(id1), striplistmap.get(id2)));
+
+ if (_debug) System.out.println(this.getClass().getSimpleName() + ": stereo hits so far " + stereohits.size());
+
+
+ }
+
+ if (_debug) {
+ System.out.println(this.getClass().getSimpleName() + ":Number of stereo hits = " + stereohits.size());
+ }
+
+ // Add the stereo hits to our list of HelicalTrackHits
+ helhits.addAll(stereohits);
+ stereoCrosses.addAll(stereohits);
+ //add LCRelation for strip hits
+ for (HelicalTrackCross cross : stereohits) {
+ for (HelicalTrackStrip strip : cross.getStrips()) {
+ hitrelations.add(new MyLCRelation(cross, stripmap.get(strip)));
+ }
+ if (_debug)
+ System.out.println(this.getClass().getSimpleName() + ": Found HTC in Layer : " + cross.Layer());
+ }
+ }
+
+ //mg...6/24/11 add new set of helical track hits that are rotated to from jlab to tracking frame
+ if (_doTranformToTracking) {
+ addRotatedHitsToEvent(event, stereoCrosses);
+
+ }
+
+
+ if( addAllStripsToEventForDebug ) addStripsToEvent(event,allstrips);
+
+
+ // Create the LCRelations between HelicalTrackHits and MC particles
+ List<LCRelation> mcrelations = new ArrayList<LCRelation>();
+ for (HelicalTrackHit hit : helhits) {
+ for (MCParticle mcp : hit.getMCParticles()) {
+ mcrelations.add(new MyLCRelation(hit, mcp));
+ }
+ }
+
+ // Put the HelicalTrackHits back into the event
+ event.put(_outname, helhits, HelicalTrackHit.class, 0);
+ if (_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": The HelicalTrackHit collection " + _outname + " has " + helhits.size() + " hits.");
+ }
+ event.put(_hitrelname, hitrelations, LCRelation.class, 0);
+ event.put(_mcrelname, mcrelations, LCRelation.class, 0);
+ }
+
+ public void addCollection(String colname) {
+ _colnames.add(colname);
+ }
+
+ public void setCollection(String colname) {
+ _colnames.clear();
+ _colnames.add(colname);
+ }
+
+ public void setCollectionNames(String[] collectionNames) {
+ _colnames = Arrays.asList(collectionNames);
+ }
+
+ /**
+ Name of the HelicalTrackHit collection to be put back in the event.
+
+ @param outname Name to use for the HelicalTrackHit collection
+ */
+ public void setOutputCollectionName(String outname) {
+ _outname = outname;
+ }
+
+ public void setHitRelationName(String hitrelname) {
+ _hitrelname = hitrelname;
+ }
+
+ public void setMCRelationName(String mcrelname) {
+ _mcrelname = mcrelname;
+ }
+
+ public void setStripMaxSeparation(double maxsep) {
+ _crosser.setMaxSeparation(maxsep);
+ }
+
+ public void setStripTolerance(double tolerance) {
+ _crosser.setTolerance(tolerance);
+ }
+
+ public void setStereoPair(String detname, int lyr1, int lyr2) {
+ _stereomap.put(makeID(detname, lyr1), makeID(detname, lyr2));
+ }
+
+ public void setTransformToTracking(boolean xfrm) {
+ _doTranformToTracking = xfrm;
+ }
+
+ private List<MCParticle> getMCParticles(TrackerCluster cluster) {
+ List<MCParticle> mcplist = new ArrayList<MCParticle>();
+ for (DigiTrackerHit dhit : cluster.getDigiHits()) {
+ // Get the elemental hits - not sure what the dif is...
+ for (DigiTrackerHit dhit2 : dhit.getElementalHits()) {
+ // Get the MCParticle and add it to the hit
+ MCParticle mcp = dhit2.getMCParticle();
+ if (mcp != null) {
+ mcplist.add(mcp);
+ }
+ }
+ }
+ return mcplist;
+ }
+
+ private HelicalTrackHit makeDigi3DHit(SiTrackerHit h) {
+
+ IDetectorElement de = h.getSensor();
+ int lyr = _ID.getLayer(de);
+ BarrelEndcapFlag be = _ID.getBarrelEndcapFlag(de);
+
+ HelicalTrackHit hit = new HelicalTrack3DHit(h.getPositionAsVector(),
+ h.getCovarianceAsMatrix(), h.getdEdx(), h.getTime(),
+ h.getRawHits(), _ID.getName(de), lyr, be);
+
+ for (MCParticle p : h.getMCParticles()) {
+ hit.addMCParticle(p);
+ }
+
+ return hit;
+
+ }
+
+ private HelicalTrackHit makeDigiAxialHit(SiTrackerHitStrip1D h) {
+
+ double z1 = h.getHitSegment().getEndPoint().z();
+ double z2 = h.getHitSegment().getStartPoint().z();
+ double zmin = Math.min(z1, z2);
+ double zmax = Math.max(z1, z2);
+ IDetectorElement de = h.getSensor();
+
+ HelicalTrackHit hit = new HelicalTrack2DHit(h.getPositionAsVector(),
+ h.getCovarianceAsMatrix(), h.getdEdx(), h.getTime(),
+ h.getRawHits(), _ID.getName(de), _ID.getLayer(de),
+ _ID.getBarrelEndcapFlag(de), zmin, zmax);
+
+ for (MCParticle p : h.getMCParticles()) {
+ hit.addMCParticle(p);
+ }
+ return hit;
+ }
+
+ private HelicalTrackStrip makeDigiStrip(SiTrackerHitStrip1D h) {
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": makeDigiStrip--");
+ System.out.println(this.getClass().getSimpleName() + ": SiTrackerHitStrip1D at " + h.getPositionAsVector().toString());
+ }
+
+ SiTrackerHitStrip1D local = h.getTransformedHit(CoordinateSystem.SENSOR);
+ SiTrackerHitStrip1D global = h.getTransformedHit(CoordinateSystem.GLOBAL);
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": local at " + local.getPositionAsVector().toString());
+ System.out.println(this.getClass().getSimpleName() + ": global at " + global.getPositionAsVector().toString());
+ }
+
+ ITransform3D trans = local.getLocalToGlobal();
+ Hep3Vector org = trans.transformed(_orgloc);
+ Hep3Vector u = global.getMeasuredCoordinate();
+ Hep3Vector v = global.getUnmeasuredCoordinate();
+ Hep3Vector w = VecOp.cross(u, v);
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": local to global using transform:");
+ System.out.println(this.getClass().getSimpleName() + ": rotation:");
+ System.out.println(trans.getRotation().getRotationMatrix().toString());
+ System.out.println(this.getClass().getSimpleName() + ": translation:");
+ System.out.println(trans.getTranslation().getTranslationVector().toString());
+ }
+
+ double umeas = local.getPosition()[0];
+ double vmin = VecOp.dot(local.getUnmeasuredCoordinate(), local.getHitSegment().getStartPoint());
+ double vmax = VecOp.dot(local.getUnmeasuredCoordinate(), local.getHitSegment().getEndPoint());
+ double du = Math.sqrt(local.getCovarianceAsMatrix().diagonal(0));
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": resulting org " + org.toString());
+ System.out.println(this.getClass().getSimpleName() + ": resulting u " + u.toString());
+ System.out.println(this.getClass().getSimpleName() + ": resulting v " + v.toString());
+ System.out.println(this.getClass().getSimpleName() + ": u.cross(v) " + w.toString());
+ System.out.println(this.getClass().getSimpleName() + ": resulting umeas " + umeas);
+ System.out.println(this.getClass().getSimpleName() + ": vmin " + vmin + " vmax " + vmax);
+
+ }
+
+ IDetectorElement de = h.getSensor();
+ String det = _ID.getName(de);
+ int lyr = _ID.getLayer(de);
+ BarrelEndcapFlag be = _ID.getBarrelEndcapFlag(de);
+
+ double dEdx = h.getdEdx();
+ double time = h.getTime();
+ List<RawTrackerHit> rawhits = h.getRawHits();
+ HelicalTrackStrip strip = new HelicalTrackStrip(org, u, v, umeas, du,
+ vmin, vmax, dEdx, time, rawhits, det, lyr, be);
+
+ try {
+ if (h.getMCParticles() != null) {
+ for (MCParticle p : h.getMCParticles()) {
+ strip.addMCParticle(p);
+ }
+ }
+ } catch (RuntimeException e) {
+ // Okay when MC info not present.
+ }
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": produce final strip at origin " + strip.origin().toString());
+ System.out.println(this.getClass().getSimpleName() + ": produce final strip with vmin " + strip.vmin() + " vmax " + strip.vmax());
+ }
+
+ return strip;
+ }
+
+ private String makeID(String detname, int lyr) {
+ return detname + lyr;
+ }
+// public HelicalTrack3DHit(Hep3Vector pos, SymmetricMatrix cov, double dEdx, double time,
+// List rawhits, String detname, int layer, BarrelEndcapFlag beflag) {
+
+ private void addRotatedHitsToEvent(EventHeader event, List<HelicalTrackCross> stereohits) {
+
+ List<HelicalTrackHit> rotatedhits = new ArrayList<HelicalTrackHit>();
+ List<LCRelation> hthrelations = new ArrayList<LCRelation>();
+ List<LCRelation> mcrelations = new ArrayList<LCRelation>();
+ for (HelicalTrackCross cross : stereohits) {
+ List<HelicalTrackStrip> rotatedstriphits = new ArrayList<HelicalTrackStrip>();
+ for (HelicalTrackStrip strip : cross.getStrips()) {
+
+ Hep3Vector origin = strip.origin();
+ Hep3Vector u = strip.u();
+ Hep3Vector v = strip.v();
+ double umeas = strip.umeas();
+ double du = strip.du();
+ double vmin = strip.vmin();
+ double vmax = strip.vmax();
+ double dedx = strip.dEdx();
+ double time = strip.time();
+ List<RawTrackerHit> rthList = strip.rawhits();
+ String detname = strip.detector();
+ int layer = strip.layer();
+ BarrelEndcapFlag bec = strip.BarrelEndcapFlag();
+ Hep3Vector neworigin = getPosTrkSystem(origin);
+ Hep3Vector newu = getPosTrkSystem(u);
+ Hep3Vector newv = getPosTrkSystem(v);
+ HelicalTrackStrip newstrip = new HelicalTrackStrip(neworigin, newu, newv, umeas, du, vmin, vmax, dedx, time, rthList, detname, layer, bec);
+ for (MCParticle p : strip.MCParticles()) {
+ newstrip.addMCParticle(p);
+ }
+ rotatedstriphits.add(newstrip);
+ }
+ HelicalTrackCross newhit = new HelicalTrackCross(rotatedstriphits.get(0), rotatedstriphits.get(1));
+ rotatedhits.add(newhit);
+ hthrelations.add(new MyLCRelation(cross, newhit));
+ for (MCParticle mcp : newhit.getMCParticles()) {
+ mcrelations.add(new MyLCRelation(newhit, mcp));
+ }
+ }
+ event.put("RotatedHelicalTrackHits", rotatedhits, HelicalTrackHit.class, 0);
+ event.put("RotatedHTHRelation", hthrelations, LCRelation.class, 0);
+ event.put("RotatedMCRelations", mcrelations, LCRelation.class, 0);
+
+ // Create the LCRelations between HelicalTrackHits and MC particles
+
+ }
+
+ private void addStripsToEvent(EventHeader event, List<HelicalTrackStrip> strips) {
+ event.put("HelicalTrackStrips", strips, HelicalTrackStrip.class,0);
+ }
+
+ private Hep3Vector getPosTrkSystem(Hep3Vector pos) {
+ return _detToTrk.transformVectorToTracking(pos);
+ }
+
+ private SymmetricMatrix getCovTrkSystem(SymmetricMatrix cov) {
+ return _detToTrk.transformCovarianceToTracking(cov);
+ }
+
+ public void setDebug(boolean debug) {
+ this._debug = debug;
+ }
+}