Commit in hps-java/sandbox on MAIN | |||
TrackerReconDriver.java | +226 | added 1.1 | |
WTrackUtils.java | +258 | added 1.1 | |
+484 |
sandboxing
diff -N TrackerReconDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackerReconDriver.java 27 Nov 2012 02:33:14 -0000 1.1 @@ -0,0 +1,226 @@
+package org.lcsim.hps.recon.tracking; + +import hep.physics.vec.BasicHep3Vector; + +import java.util.ArrayList; +import java.util.List; + +import org.lcsim.event.EventHeader; +import org.lcsim.event.Track; +import org.lcsim.fit.helicaltrack.HelicalTrackHit; +import org.lcsim.geometry.Detector; +import org.lcsim.hps.event.HPSTransformations; +import org.lcsim.recon.tracking.digitization.sisim.config.ReadoutCleanupDriver; +import org.lcsim.recon.tracking.seedtracker.SeedStrategy; +import org.lcsim.recon.tracking.seedtracker.SeedTracker; +import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils; +import org.lcsim.recon.tracking.seedtracker.diagnostic.SeedTrackerDiagnostics; +import org.lcsim.util.Driver; + +/** + * This class runs the Track Reconstruction for the HPS Test Proposal detector. + * The tracker digitization must be run in front of it. It is intended to work + * with the {@link TrackerDigiDriver} digitization Driver. + * + * @author jeremym + * @version $Id: TrackerReconDriver.java,v 1.18 2012/05/03 20:50:29 mgraham Exp + * $ + */ +public final class TrackerReconDriver extends Driver { + + // Debug flag. + private boolean debug = false; + // Tracks found across all events. + int ntracks = 0; + // Number of events processed. + int nevents = 0; + // Cache detector object. + Detector detector = null; + // Default B-field value. + private double bfield = 0.5; + // TrackerHit readout name for readout cleanup. + private String trackerReadoutName = "TrackerHits"; + // Tracking strategies resource path. + private String strategyResource = "HPS-Test-4pt1.xml"; + // Output track collection. + private String trackCollectionName = "MatchedTracks"; + // HelicalTrackHit input collection. + private String stInputCollectionName = "RotatedHelicalTrackHits"; + HPSHelicalTrackHitDriver hthdriver = new HPSHelicalTrackHitDriver(); + private boolean runHelicalTrackHitDriver = true; + + public TrackerReconDriver() { + } + + public void setDebug(boolean debug) { + this.debug = debug; + hthdriver.setDebug(debug); + } + + public void setRunHelicalTrackHitDriver(boolean runHelicalTrackHitDriver) { + this.runHelicalTrackHitDriver = runHelicalTrackHitDriver; + } + + public void setSubdetectorName(String subdetectorName) { + hthdriver.setSubdetectorName(subdetectorName); + } + + public void setStripHitsCollectionName(String stripHitsCollectionName) { + hthdriver.HitRelationName(stripHitsCollectionName); + } + + public void setHelicalTrackHitRelationsCollectionName(String helicalTrackHitRelationsCollectionName) { + hthdriver.HitRelationName(helicalTrackHitRelationsCollectionName); + } + + public void setHelicalTrackMCRelationsCollectionName(String helicalTrackMCRelationsCollectionName) { + hthdriver.MCRelationName(helicalTrackMCRelationsCollectionName); + } + + public void setOutputHitCollectionName(String outputHitCollectionName) { + hthdriver.setOutputCollectionName(outputHitCollectionName); + } + + public void setStripMaxSeparation(double stripMaxSeparation) { + hthdriver.setMaxSeperation(stripMaxSeparation); + } + + public void setStripTolerance(double stripTolerance) { + hthdriver.setTolerance(stripTolerance); + } + + public void setLayerGeometryType(String geomType) { + hthdriver.setLayerGeometryType(geomType); + } + + /** + * Set the tracking strategy resource. + * + * @param strategyResource The absolute path to the strategy resource in the + * hps-java jar. + */ + public void setStrategyResource(String strategyResource) { + this.strategyResource = strategyResource; + } + + public void setInputHitCollectionName(String inputHitCollectionName) { + this.stInputCollectionName = inputHitCollectionName; + } + + public void setTrackCollectionName(String trackCollectionName) { + this.trackCollectionName = trackCollectionName; + } + + /** + * Set the SimTrackerHit collection to be used for tracking. + * + * @param simTrackerHitCollectionName The name of the SimTrackerHit + * collection in the event. + */ + public void setSimTrackerHitCollectionName(String simTrackerHitCollectionName) { + hthdriver.setCollection(simTrackerHitCollectionName); + } + + /** + * This is used to setup the Drivers after XML config. + */ + @Override + public void detectorChanged(Detector detector) { + // Cache Detector object. + this.detector = detector; + + // Get B-field Y with no sign. Seed Tracker doesn't like signed B-field components. + // FIXME Is this always right? + this.bfield = Math.abs((detector.getFieldMap().getField(new BasicHep3Vector(0, 0, 0)).y())); + if (debug) { + System.out.println("Set B-field to " + this.bfield); + } + + initialize(); + + super.detectorChanged(detector); + } + + /** + * Setup all the child Drivers necessary for track reconstruction. + */ + private void initialize() { + // + // 1) Driver to create HelicalTrackHits expected by Seedtracker. + // + // TODO Make this step its own separate Driver??? (Matt) + + if (runHelicalTrackHitDriver) { + add(hthdriver); + } + + // + // 2) Driver to run Seed Tracker. + // + + if (!strategyResource.startsWith("/")) { + strategyResource = "/org/lcsim/hps/recon/tracking/strategies/" + strategyResource; + } + List<SeedStrategy> sFinallist = StrategyXMLUtils.getStrategyListFromInputStream(this.getClass().getResourceAsStream(strategyResource)); + SeedTracker stFinal = new SeedTracker(sFinallist); + HPSTransformations hpstrans = new HPSTransformations(); + stFinal.setMaterialManagerTransform(hpstrans.getTransform()); + stFinal.setInputCollectionName(stInputCollectionName); + stFinal.setTrkCollectionName(trackCollectionName); + stFinal.setBField(bfield); + stFinal.setDiagnostics(new SeedTrackerDiagnostics()); +// stFinal.setSectorParams(false); //this doesn't actually seem to do anything + stFinal.setSectorParams(1, 10000); + add(stFinal); + + // + // 3) Cleanup the readouts for next event. + // + List<String> readoutCleanup = new ArrayList<String>(); + readoutCleanup.add(this.trackerReadoutName); + add(new ReadoutCleanupDriver(readoutCleanup)); + } + + /** + * This method is used to run the reconstruction and print debug + * information. + */ + @Override + public void process(EventHeader event) { + // This call runs the track reconstruction using the sub-Drivers. + super.process(event); + + // Debug printouts. + if (debug) { + if(event.hasCollection(HelicalTrackHit.class,stInputCollectionName)) { + System.out.println(this.getClass().getSimpleName() + ": The HelicalTrackHit collection " + stInputCollectionName + " has " + event.get(HelicalTrackHit.class, stInputCollectionName).size() + " hits."); + } else { + System.out.println(this.getClass().getSimpleName() + ": No HelicalTrackHit collection for this event"); + } + // Check for Tracks. + List<Track> tracks = event.get(Track.class, trackCollectionName); + System.out.println(this.getClass().getSimpleName() + ": The Track collection " + trackCollectionName + " has " + tracks.size() + " tracks."); + + // Print out track info. + for (Track track : tracks) { + System.out.println(this.getClass().getSimpleName() + ": " + track.toString()); + System.out.println(this.getClass().getSimpleName() + ": number of layers = " + track.getTrackerHits().size()); + System.out.println(this.getClass().getSimpleName() + ": chi2 = " + track.getChi2()); + } + } + + // Increment number of events. + ++nevents; + + // Add to tracks found. + ntracks += event.get(Track.class, trackCollectionName).size(); + } + + @Override + public void endOfData() { + if (debug) { + System.out.println("-------------------------------------------"); + System.out.println(this.getName() + " found " + ntracks + " tracks in " + nevents + " events which is " + ((double) ntracks / (double) nevents) + " tracks per event."); + } + } +}
\ No newline at end of file
diff -N WTrackUtils.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ WTrackUtils.java 27 Nov 2012 02:33:14 -0000 1.1 @@ -0,0 +1,258 @@
+/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package org.lcsim.hps.users.phansson; + +import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; +import java.util.ArrayList; +import java.util.List; + +/** + * + * @author phansson + */ +public class WTrackUtils { + private int _debug = 0; + List<Double> _delta_point = new ArrayList<Double>(); + + public WTrackUtils() { + } + + public WTrackUtils(boolean debug) { + _debug = debug?1:0; + } + + public WTrackUtils(int debug) { + _debug = debug; + } + + + public void setDebug(boolean debug) { + _debug = debug?1:0; + } + + public boolean getDebug() { + return _debug>0?true:false; + } + + public List<Double> getDeltas() { + return _delta_point; + } + + private static Hep3Vector getMomentumOnHelix(WTrack track, double s) { + double a = track.getBfieldConstant(); + Hep3Vector p0 = track.getP0(); + double rho = a / p0.magnitude(); + double px = p0.x()*Math.cos(rho*s) - p0.y()*Math.sin(rho*s); + double py = p0.y()*Math.cos(rho*s) + p0.x()*Math.sin(rho*s); + double pz = p0.z(); + return (new BasicHep3Vector(px,py,pz)); + } + + private static Hep3Vector getPointOnHelix(WTrack track, double s) { + double a = track.getBfieldConstant(); + Hep3Vector p0 = track.getP0(); + Hep3Vector x0 = track.getX0(); + double rho = a / p0.magnitude(); + double x = x0.x() + p0.x()/a*Math.sin(rho*s) - p0.y()/a*(1-Math.cos(rho*s)); + double y = x0.y() + p0.y()/a*Math.sin(rho*s) + p0.x()/a*(1-Math.cos(rho*s)); + double z = x0.z() + p0.z()/p0.magnitude()*s; + return (new BasicHep3Vector(x,y,z)); + } + + private static double getPathLengthToPlaneApprox(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) { + /* + * Find the approximate path length to the point xp + * in arbitrary oriented, constant magnetic field with unit vector h + */ + boolean debug = false; + double a = track.getBfieldConstant(); + Hep3Vector p0 = track.getP0(); + Hep3Vector x0 = track.getX0(); + double p = p0.magnitude(); + double rho = a / p; + double A = VecOp.dot(eta,VecOp.cross(p0, h))/p*0.5*rho; + double B = VecOp.dot(p0,eta)/p; + double C = VecOp.dot(VecOp.sub(x0,xp),eta); + double t = B*B-4*A*C; + if(t<0) { + System.out.println(" getPathLengthToPlaneApprox ERROR t is negative (" + t + ")!" ); + System.out.println(" p " + p + " rho " + rho + " A " + A + " B " + B + " C " + C); + System.out.println(" track params\n" + track.toString()); + System.out.println(" xp " + xp.toString()); + System.out.println(" eta " + eta.toString()); + System.out.println(" h " + h.toString()); + System.exit(1); + } + double root1 = (-B + Math.sqrt(t)) /(2*A); + double root2 = (-B - Math.sqrt(t)) /(2*A); + + // choose the smallest positive solution + // if both negative choose the smallest negative ??? + //if(root1==0 || root2==0) root=0; + double root = Math.abs(root1) <= Math.abs(root2) ? root1 : root2; +// else if(Math.signum(root1)>0 && Math.signum(root2)<0) root = root1; +// else if(Math.signum(root2)>0 && Math.signum(root1)<0) root = root2; +// else if(Math.signum(root1)>0 && Math.signum(root2)>0) root = root1 > root2 ? root2 : root1; +// else if(Math.signum(root1)<0 && Math.signum(root2)<0) root = root1 < root2 ? root2 : root1; +// else { +// System.out.println(" I should never get here! (root1=" + root1 + " root2=" + root2+")"); +// System.exit(1); +// } + if((debug)) { + System.out.println(" getPathLengthToPlaneApprox "); + System.out.println(" " + track.toString()); + System.out.println(" xp " + xp.toString()); + System.out.println(" eta " + eta.toString()); + System.out.println(" h " + h.toString()); + System.out.println(" p " + p + " rho " + rho + " t " + t + " A " + A + " B " + B + " C " + C); + System.out.println(" root1 " + root1 + " root2 " + root2 + " -> root " + root); + } + return root; + + } + + + private static Hep3Vector getPointOnHelix(WTrack track, double s, Hep3Vector h) { + /* + * Get point on helix at path lenght s + * in arbitrary oriented, constant magnetic field with unit vector h + */ + + double a = track.getBfieldConstant(); + Hep3Vector p0 = track.getP0(); + double p = p0.magnitude(); + Hep3Vector x0 = track.getX0(); + double rho = a / p0.magnitude(); + double srho = s*rho; + Hep3Vector a1 = VecOp.mult(1/a*Math.sin(srho), p0); + Hep3Vector a2 = VecOp.mult(1/a*(1-Math.cos(srho)),VecOp.cross(p0,h)); + Hep3Vector a3 = VecOp.mult(VecOp.dot(p0, h)/p,h); + Hep3Vector a4 = VecOp.mult(s-Math.sin(srho)/rho, a3); + Hep3Vector x = VecOp.add(x0,a1); + x = VecOp.sub(x,a2); + x = VecOp.add(x,a4); + return x; + } + + private static Hep3Vector getMomentumOnHelix(WTrack track, double s, Hep3Vector h) { + /* + * Get point on helix at path lenght s + * in arbitrary oriented, constant magnetic field with unit vector h + */ + + double a = track.getBfieldConstant(); + Hep3Vector p0 = track.getP0(); + double rho = a / p0.magnitude(); + double srho = s*rho; + Hep3Vector a1 = VecOp.mult(Math.cos(srho), p0); + Hep3Vector a2 = VecOp.cross(p0, VecOp.mult(Math.sin(srho),h)); + Hep3Vector a3 = VecOp.mult(VecOp.dot(p0,h),VecOp.mult(1-Math.cos(srho),h)); + Hep3Vector p = VecOp.sub(a1,a2); + p = VecOp.add(p,a3); + return p; + } + + + private static double[] getHelixParametersAtPathLength(WTrack track, double s, Hep3Vector h) { + /* + * Calculate the exact position of the new helix parameters at path length s + * in an arbitrarily oriented, constant magnetic field + * + * point xp is the point + * h is a unit vector in the direction of the magnetic field + */ + boolean debug = false; + // Find track parameters at that path length + Hep3Vector p = getMomentumOnHelix(track, s, h); + Hep3Vector x = getPointOnHelix(track, s, h); + + Hep3Vector p_tmp = getMomentumOnHelix(track, s); + Hep3Vector x_tmp = getPointOnHelix(track, s); + + if((debug)) { + System.out.println(" point on helix at s"); + System.out.println(" p " + p.toString() + " p_tmp " + p_tmp.toString()); + System.out.println(" x " + x.toString() + " x_tmp " + x_tmp.toString()); + } + + + //Create the new parameter array + double [] pars = new double[7]; + pars[0] = p.x(); + pars[1] = p.y(); + pars[2] = p.z(); + pars[3] = track.getParameters()[3]; //E is unchanged + pars[4] = x.x(); + pars[5] = x.y(); + pars[6] = x.z(); + return pars; + } + + + + public static Hep3Vector getHelixAndPlaneIntercept(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) { + + /* + * Find the interception point between the helix and plane + * track: the helix in W representation + * xp: point on the plane + * eta: unit vector of the plane + * h: unit vector of magnetic field + */ + boolean debug = false; + + int max_it = 10; + double epsilon = 1e-4; + int iteration = 1; + double delta_s = 9999; + double s_prev = 999999; + List<WTrack> tracks = new ArrayList<WTrack>(); + WTrack trk = new WTrack(track); + while(iteration<=max_it && Math.abs(s_prev)>epsilon) { + if((debug)) { + System.out.println("iteration " + iteration + " --- "); + System.out.println(" s_prev " + s_prev + " delta_s " + delta_s); + System.out.println(trk.toString()); + } + // Start by approximating the path length to the point + double s = getPathLengthToPlaneApprox(trk, xp, eta, h); + + if((debug)) { + System.out.println(": Approx. path s " + s); + } + + // Find the track parameters at this point + double[] params = getHelixParametersAtPathLength(trk, s, h); + // update the track parameters + trk.setTrackParameters(params); + + if((debug)) { + System.out.println(" updated params " + trk.toString()); + } + + tracks.add(trk); + iteration++; + delta_s = Math.abs(s-s_prev); + s_prev = s; + + //Save distance between point and estimate + Hep3Vector dpoint = VecOp.sub(xp, trk.getX0()); + + + } + + if((debug)) { + System.out.println(" Final s " + s_prev + " delta_s " + delta_s); + System.out.println(" Final track params " + trk.toString()); + } + + Hep3Vector point = trk.getX0(); + return point; + + } + +}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1