hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N ParticleHelixProducer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ParticleHelixProducer.java 9 Oct 2012 01:22:28 -0000 1.1
@@ -0,0 +1,183 @@
+package org.lcsim.hps.users.phansson;
+
+import org.lcsim.hps.recon.tracking.*;
+import hep.aida.*;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Matrix;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.*;
+import org.lcsim.constants.Constants;
+import org.lcsim.hps.users.mgraham.*;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IReadout;
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.Transform3D;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.tracker.silicon.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseRawTrackerHit;
+import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.fit.helicaltrack.MultipleScatter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.lcsim.hps.recon.tracking.*;
+import org.lcsim.hps.users.meeg.NoiselessReadoutChip;
+import org.lcsim.hps.users.phansson.WTrack;
+import org.lcsim.hps.users.phansson.WTrackUtils;
+import org.lcsim.recon.tracking.digitization.sisim.*;
+import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType.CoordinateSystem;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.lcio.LCIOUtil;
+
+/**
+ *
+ * @author phansson
+ */
+public class ParticleHelixProducer extends Driver {
+
+ private boolean debug = false;
+ private int _totalTracks = 0;
+ TrackUtils trackUtils = new TrackUtils();
+ TrackerHitUtils trackerhitutils = new TrackerHitUtils();
+ Hep3Matrix detToTrk;
+ SvtTrackExtrapolator extrapolator = new SvtTrackExtrapolator();
+ Hep3Vector _bfield;
+ WTrackUtils wutils = new WTrackUtils();
+ TrackerHitType trackerType = new TrackerHitType(TrackerHitType.CoordinateSystem.GLOBAL, TrackerHitType.MeasurementType.STRIP_1D);
+ CoordinateSystem coordinate_system = trackerType.getCoordinateSystem();
+
+
+ // Name of StripHit1D output collection.
+ private String trackOutputCollectionName = "MCParticle_HelicalTrackFit";
+
+
+ public void setDebug(boolean debug) {
+ this.debug = debug;
+ }
+
+ public void setTrackOutputCollectionName(String trackOutputCollectionName) {
+ this.trackOutputCollectionName = trackOutputCollectionName;
+ }
+
+ /**
+ * Creates a new instance of TrackerHitDriver.
+ */
+ public ParticleHelixProducer() {
+ }
+
+ /**
+ * Do initialization once we get a Detector.
+ */
+ @Override
+ public void detectorChanged(Detector detector) {
+
+ // Call sub-Driver's detectorChanged methods.
+ super.detectorChanged(detector);
+
+
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+ _bfield = new BasicHep3Vector(0,0,detector.getFieldMap().getField(IP).y());
+ detToTrk = trackerhitutils.detToTrackRotationMatrix();
+ this.wutils.setDebug(true);
+
+
+ }
+
+ /**
+ * Perform the digitization.
+ */
+ @Override
+ public void process(EventHeader event) {
+
+
+ //Make new tracks based on the MC particles
+ List<HelicalTrackFit> tracks = new ArrayList<HelicalTrackFit>();
+ if(event.hasCollection(MCParticle.class)) {
+ List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": Number of MC particles = " + mcparticles.size());
+ List<MCParticle> fsParticles = HPSMCParticlePlotsDriver.makeGenFSParticleList(mcparticles);
+ if(debug) System.out.println(this.getClass().getSimpleName()+": Number of FS MC particles = " + fsParticles.size());
+ double bfield = Math.abs(_bfield.z()); //remove sign from B-field, assumed to be along z-direction
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": bfield " + bfield );
+ for(MCParticle part : fsParticles) {
+ if(ParticleTypeClassifier.isElectron(part.getPDGID()) || ParticleTypeClassifier.isPositron(part.getPDGID())) {
+
+
+ Hep3Vector p = part.getMomentum();
+ Hep3Vector org = part.getOrigin();
+ double q = -1*part.getCharge(); //since I flipped the B-field I need to flip the charge
+ //if(debug) System.out.println(this.getClass().getSimpleName() + ": MC particle p=" + p.toString() +" org=" + org.toString() +" q = " + q);
+ p = VecOp.mult(detToTrk, p);
+ org = VecOp.mult(detToTrk, org);
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": MC particle p=" + p.toString() +" org=" + org.toString() +" q = " + q + " (rotated)");
+
+ if(p.magnitude()<0.3) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": this MC particle had too small momentum p=" + p.toString());
+ continue;
+ }
+ if(debug) {
+ double pt = Math.sqrt(p.x()*p.x() + p.y()*p.y());
+ double Rman = q*pt/(Constants.fieldConversion*bfield);
+ double phi = Math.atan2(p.y(),p.x());
+ double xc = org.x() + Rman*Math.sin(phi);
+ double yc = org.y() - Rman*Math.cos(phi);
+ double Rc = Math.sqrt(xc*xc+yc*yc);
+ double dca = q >0 ? (Rman - Rc) : (Rman + Rc);
+ System.out.println(this.getClass().getSimpleName() + ": pt " + pt +" R " + Rman + " phi " + phi + " xc " + xc + " yc " + yc + " Rc " + Rc + " DCA " + dca );
+ }
+ //HelixParamCalculator hpc = new HelixParamCalculator(part,bfield);
+ HelixParamCalculator hpc = new HelixParamCalculator(p, org, (int)q,bfield); //remove sign from b-field
+ double R = hpc.getRadius();
+ double slope = hpc.getSlopeSZPlane();
+ double d0 = hpc.getDCA();
+ double phi0 = hpc.getPhi0();
+ double z0 = hpc.getZ0();
+ double[] pars = new double[5];
+ pars[0] = d0;
+ pars[1] = phi0;
+ pars[2] = 1/R;
+ pars[3] = z0;
+ pars[4] = slope;
+ HelicalTrackFit htf = this.trackUtils.makeHelicalTrackFit(pars);
+ tracks.add(htf);
+ if(debug) {
+ System.out.println(this.getClass().getSimpleName() + ": MC particle created HelicalTrackFit " + htf.toString());
+
+ }
+
+
+ }
+ }
+ }
+ else {
+ System.out.println(this.getClass().getSimpleName() + ": no MC particles in this event");
+ }
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": created " + tracks.size() + " MC particle helix tracks");
+
+ event.put(this.trackOutputCollectionName, tracks, HelicalTrackFit.class, 0);
+ _totalTracks += tracks.size();
+
+
+
+ }
+
+
+
+
+
+ @Override
+ public void endOfData() {
+
+ System.out.println(this.getClass().getSimpleName() + ": produced " + _totalTracks);
+ }
+}