hps-java/src/main/java/org/lcsim/hps/recon/tracking
diff -N DataTrackerFakeHitDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DataTrackerFakeHitDriver.java 9 Oct 2012 01:21:41 -0000 1.1
@@ -0,0 +1,531 @@
+package org.lcsim.hps.recon.tracking;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IPlotter;
+import hep.aida.IProfile1D;
+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.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.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseRawTrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.hps.users.phansson.WTrack;
+import org.lcsim.hps.users.phansson.WTrackUtils;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType;
+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 DataTrackerFakeHitDriver extends Driver {
+
+ private boolean debug = false;
+ 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();
+
+ private String trackCollectionName = "MCParticle_HelicalTrackFit";
+
+
+ // Subdetector name.
+ private String subdetectorName = "Tracker";
+ // Name of StripHit1D output collection.
+ private String stripHitOutputCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+ // Various data lists required by digitization.
+ private List<String> processPaths = new ArrayList<String>();
+ private List<IDetectorElement> processDEs = new ArrayList<IDetectorElement>();
+ private Set<SiSensor> processSensors = new HashSet<SiSensor>();
+
+ //Visualization
+ private boolean hideFrame = false;
+ private AIDA aida = AIDA.defaultInstance();
+ private HashMap<SiSensor,IProfile1D> _delta_histos = new HashMap<SiSensor,IProfile1D>();
+ private HashMap<SiSensor,IHistogram1D> _delta_itercount = new HashMap<SiSensor,IHistogram1D>();
+ IProfile1D _prf_final_deltas;
+ IProfile1D _prf_all_deltas;
+ IAnalysisFactory af = aida.analysisFactory();
+ IPlotter plotter_iter;
+ IPlotter plotter_itercount;
+ IPlotter plotter_iter_final;
+ IPlotter plotter_iter_all;
+
+
+ int[][] counts = new int[2][10];
+
+ public void setDebug(boolean debug) {
+ this.debug = debug;
+ }
+
+// public void setReadoutCollectionName(String readoutCollectionName) {
+// this.readoutCollectionName = readoutCollectionName;
+// }
+ public void setSubdetectorName(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void setStripHitOutputCollectionName(String stripHitOutputCollectionName) {
+ this.stripHitOutputCollectionName = stripHitOutputCollectionName;
+ }
+
+ public void setTrackCollectionName(String trackCollectionName) {
+ this.trackCollectionName = trackCollectionName;
+ }
+
+
+ /**
+ * Creates a new instance of TrackerHitDriver.
+ */
+ public DataTrackerFakeHitDriver() {
+ }
+
+ /**
+ * 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(this.debug);
+
+
+ // Process detectors specified by path, otherwise process entire
+ // detector
+ IDetectorElement deDetector = detector.getDetectorElement();
+
+ for (String path : processPaths) {
+ processDEs.add(deDetector.findDetectorElement(path));
+ }
+
+ if (processDEs.isEmpty()) {
+ processDEs.add(deDetector);
+ }
+
+ for (IDetectorElement detectorElement : processDEs) {
+ processSensors.addAll(detectorElement.findDescendants(SiSensor.class));
+
+
+ }
+
+
+
+ // Set the detector to process.
+ processPaths.add(subdetectorName);
+
+ this.makePlots();
+
+ }
+
+ /**
+ * 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(HelicalTrackFit.class, trackCollectionName)) {
+ tracks = event.get(HelicalTrackFit.class, trackCollectionName);
+ }
+
+
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": found " + tracks.size() + " tracks (" + this.trackCollectionName + ")");
+
+ if(tracks.isEmpty()) return;
+
+ // Make new lists for output.
+ List<SiTrackerHit> stripHits1D = new ArrayList<SiTrackerHit>();
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": Add hits for " + tracks.size() + " tracks (" + this.trackCollectionName + ")");
+
+
+ for(HelicalTrackFit helix : tracks) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": trying to add hits for this track");
+
+ if(debug) {
+ System.out.println(this.getClass().getSimpleName() + helix.toString());
+ System.out.println(this.getClass().getSimpleName() + ": htf x0 " + helix.x0() + " y0 " + helix.y0());
+ }
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": create a WTrack object");
+
+ WTrack wtrack = new WTrack(helix,Math.abs(_bfield.z())); //remove sign from B-field (assumed to go along z-direction)
+
+ if(debug) {
+ System.out.println(this.getClass().getSimpleName() + ": " + wtrack.toString());
+
+ }
+
+
+
+ extrapolator.setTrack(helix.parameters());
+
+ // Make hits if the helix passed through the sensor
+ for (SiSensor sensor : processSensors) {
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": add hits to sensor " + sensor.getName());
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": sensor " + sensor.getName() + " position at " + sensor.getGeometry().getPosition().toString());
+
+ double zposition = sensor.getGeometry().getPosition().z();
+
+ Hep3Vector trackPosAtSensor = extrapolator.extrapolateTrack(zposition);
+
+ if(debug) {
+ System.out.println(this.getClass().getSimpleName() + ": track position at sensor ( z= " + zposition + ") = " + trackPosAtSensor.toString());
+ }
+
+
+
+ boolean isHit = trackUtils.sensorContainsTrack(trackPosAtSensor, sensor);
+
+
+ if(isHit) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": make a tracker hit and add to this sensor");
+
+ stripHits1D.add(this.makeTrackerHit(sensor, wtrack));
+ } else {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": this helix didn't pass within the sensor so no hit was added");
+ }
+ }
+ }
+
+
+ if (debug) {
+ System.out.println("TrackerHit collection " + this.stripHitOutputCollectionName + " has " + stripHits1D.size() + " hits.");
+ }
+
+ // Put output hits into collection.
+ int flag = LCIOUtil.bitSet(0, 31, true); // Turn on 64-bit cell ID.
+// event.put(this.rawTrackerHitOutputCollectionName, rawHits, RawTrackerHit.class, flag, toString());
+ event.put(this.stripHitOutputCollectionName, stripHits1D, SiTrackerHitStrip1D.class, 0, toString());
+ if (debug) {
+ for (int mod = 0; mod < 2; mod++) {
+ for (int layer = 0; layer < 10; layer++) {
+ counts[mod][layer] += SvtUtils.getInstance().getSensor(mod, layer).getReadout().getHits(SiTrackerHit.class).size();
+ }
+ }
+ }
+ }
+
+
+ private SiTrackerHitStrip1D makeTrackerHit(SiSensor sensor, WTrack wtrack) {
+ //private SiTrackerHitStrip1D makeTrackerHit(List<HPSFittedRawTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ //create fake raw tracker hit
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": makeTrackerHit");
+
+ List<RawTrackerHit> rth_cluster = this.makeRawTrackerFakeHit(sensor);
+ if(rth_cluster.size()!=1) {
+ System.out.println(this.getClass().getSimpleName() + ": the fake raw tracker hit cluster is different than one!? " + rth_cluster.size() );
+ System.exit(1);
+ }
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": created a fake raw tracker hit ");
+ RawTrackerHit raw_hit = rth_cluster.get(0);
+ IIdentifier id = raw_hit.getIdentifier();
+
+ //Get the electrode objects
+ SiTrackerIdentifierHelper _sid_helper = (SiTrackerIdentifierHelper) sensor.getIdentifierHelper();
+ ChargeCarrier carrier = ChargeCarrier.getCarrier(_sid_helper.getSideValue(id));
+ SiSensorElectrodes electrodes = ((SiSensor) raw_hit.getDetectorElement()).getReadoutElectrodes(carrier);
+
+
+
+
+ ITransform3D local_to_global;
+ if (coordinate_system == TrackerHitType.CoordinateSystem.GLOBAL) {
+ local_to_global = new Transform3D();
+ }
+ else if (coordinate_system == TrackerHitType.CoordinateSystem.SENSOR) {
+ local_to_global = sensor.getGeometry().getLocalToGlobal();
+ }
+ else {
+ throw new RuntimeException(this.getClass().getSimpleName() + " problem with coord system " + coordinate_system.toString());
+ }
+
+
+ ITransform3D electrodes_to_global = electrodes.getLocalToGlobal();
+ ITransform3D global_to_hit = local_to_global.inverse();
+ ITransform3D electrodes_to_hit = Transform3D.multiply(global_to_hit,electrodes_to_global);
+
+ Hep3Vector u= electrodes_to_hit.rotated(electrodes.getMeasuredCoordinate(0));
+ Hep3Vector v= electrodes_to_hit.rotated(electrodes.getUnmeasuredCoordinate(0));
+ Hep3Vector w = VecOp.cross(u, v);
+ Hep3Vector _orgloc = new BasicHep3Vector(0., 0., 0.);
+ electrodes_to_global.transformed(_orgloc);
+ if(debug) {
+ System.out.println(this.getClass().getSimpleName() + ": electrodes u " + u.toString());
+ System.out.println(this.getClass().getSimpleName() + ": electrodes v " + v.toString());
+ System.out.println(this.getClass().getSimpleName() + ": electrodes w " + w.toString() + "( " + w.magnitude() + ")");
+ }
+
+ electrodes_to_global.getTranslation().translate(_orgloc);
+ if(debug) System.out.print(this.getClass().getSimpleName() + ": orgloc " + _orgloc.toString() + " -> ");
+ _orgloc = VecOp.mult(detToTrk, _orgloc);
+ if(debug) System.out.println(_orgloc.toString());
+
+
+
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": Try to find the interception point");
+
+ //B-field must go along Z direction
+ Hep3Vector h = new BasicHep3Vector(_bfield.x(),_bfield.y(),Math.abs(_bfield.z()));
+ h = VecOp.unit(h);
+ //Rotate into tracking frame
+ Hep3Vector eta = VecOp.mult(detToTrk, w);
+
+
+ Hep3Vector position = wutils.getHelixAndPlaneIntercept(wtrack, _orgloc, eta, h);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": found interception point at position " + position.toString());
+
+ position = VecOp.mult(VecOp.inverse(detToTrk), position);
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": rotate the hit position to the global frame -> " + position.toString());
+
+
+
+ //Fill some debug stuff
+ List<Double> deltas = wutils.getDeltas();
+ this._delta_itercount.get(sensor).fill(deltas.size());
+ IProfile1D dhisto = this._delta_histos.get(sensor);
+ for(int i=0;i<deltas.size();++i) {
+ this._prf_all_deltas.fill(i, deltas.get(i));
+ dhisto.fill(i,deltas.get(i));
+ }
+ _prf_final_deltas.fill(deltas.size(), deltas.get(deltas.size()-1));
+
+
+ SymmetricMatrix covariance = this.getCovariance(rth_cluster, electrodes);
+ double time = getTime(rth_cluster);
+ double energy = getEnergy(rth_cluster);
+
+
+ SiTrackerHitStrip1D hit = new SiTrackerHitStrip1D(position, covariance, energy, time, rth_cluster, trackerType);
+ return hit;
+
+
+ }
+
+
+
+
+
+
+ //private SymmetricMatrix getCovariance() {
+ private SymmetricMatrix getCovariance(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ SymmetricMatrix covariance = new SymmetricMatrix(3);
+ covariance.setElement(0, 0, Math.pow(getMeasuredResolution(cluster, electrodes), 2));
+ covariance.setElement(1, 1, Math.pow(getUnmeasuredResolution(cluster, electrodes), 2));
+ covariance.setElement(2, 2, 0.0);
+
+ SymmetricMatrix covariance_global = electrodes.getLocalToGlobal().transformed(covariance);
+
+// System.out.println("Global covariance matrix: \n"+covariance_global);
+
+ return covariance_global;
+
+ }
+
+
+
+ private double getMeasuredResolution(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes) // should replace this by a ResolutionModel class that gives expected resolution. This could be a big job.
+ {
+ double measured_resolution;
+
+ double sense_pitch = ((SiSensor) electrodes.getDetectorElement()).getSenseElectrodes(electrodes.getChargeCarrier()).getPitch(0);
+
+// double readout_pitch = electrodes.getPitch(0);
+// double noise = _readout_chip.getChannel(strip_number).computeNoise(electrodes.getCapacitance(strip_number));
+// double signal_expected = (0.000280/DopedSilicon.ENERGY_EHPAIR) *
+// ((SiSensor)electrodes.getDetectorElement()).getThickness(); // ~280 KeV/mm for thick Si sensors
+ double _oneClusterErr = 1 / Math.sqrt(12);
+ double _twoClusterErr = 1 / 5;
+ double _threeClusterErr = 1 / 3;
+ double _fourClusterErr = 1 / 2;
+ double _fiveClusterErr = 1;
+
+ if (cluster.size() == 1) {
+ measured_resolution = sense_pitch * _oneClusterErr;
+ } else if (cluster.size() == 2) {
+ measured_resolution = sense_pitch * _twoClusterErr;
+ } else if (cluster.size() == 3) {
+ measured_resolution = sense_pitch * _threeClusterErr;
+ } else if (cluster.size() == 4) {
+ measured_resolution = sense_pitch * _fourClusterErr;
+ } else {
+ measured_resolution = sense_pitch * _fiveClusterErr;
+ }
+
+ return measured_resolution;
+ }
+
+ private double getUnmeasuredResolution(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ // Get length of longest strip in hit
+ double hit_length = 0;
+ for (RawTrackerHit hit : cluster) {
+ hit_length = Math.max(hit_length, ((SiStrips) electrodes).getStripLength(1));
+ }
+ return hit_length / Math.sqrt(12);
+ }
+
+ private double getTime(List<RawTrackerHit> cluster) {
+ return 0;
+ }
+
+ private double getEnergy(List<RawTrackerHit> cluster) {
+ double total_charge = 20000;
+//
+// double total_charge = 0.0;
+// for (RawTrackerHit hit : cluster) {
+// double signal = hit.getAmp();
+// total_charge += signal;
+// }
+ return total_charge * DopedSilicon.ENERGY_EHPAIR;
+ }
+
+
+
+
+
+ public List<RawTrackerHit> makeRawTrackerFakeHit(SiSensor sensor)
+ {
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": makeRawTrackerFakeHit for sensor " + sensor.getName());
+ List<RawTrackerHit> raw_hits = new ArrayList<RawTrackerHit>();
+
+ // Get SimTrackerHits
+ IReadout ro = sensor.getReadout();
+
+
+ // Loop over electrodes and digitize with readout chip
+ for (ChargeCarrier carrier : ChargeCarrier.values())
+ {
+ if (sensor.hasElectrodesOnSide(carrier))
+ {
+
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": creating a dummy hit for sensor " + sensor.getName());
+ //SortedMap<Integer,List<Integer>> digitized_hits = _readout_chip.readout(electrode_data.get(carrier),sensor.getReadoutElectrodes(carrier));
+ //if(debug) System.out.println(this.getClass().getSimpleName() + ": creating a dummy hit for sensor " + sensor.getName());
+
+ int channel = 1;
+ int time = 0;
+ long cell_id = sensor.makeStripId(channel,carrier.charge()).getValue();
+ //List<Integer> readout_data = new ArrayList<Integer>();
+ short[] adc_values = new short[6];
+ for(Integer i=0;i<6;++i) {
+ Integer adc = 50;
+ adc_values[i] = adc.shortValue(); //ADC counts
+ }
+ IDetectorElement detector_element = sensor;
+ RawTrackerHit raw_hit = new BaseRawTrackerHit(time,cell_id,adc_values,new ArrayList<SimTrackerHit>(),detector_element);
+ ro.addHit(raw_hit);
+ raw_hits.add(raw_hit);
+
+ }
+
+ }
+
+ return raw_hits;
+ }
+
+
+
+
+ private void makePlots() {
+
+ for (SiSensor sensor : processSensors) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": Making plots for " + sensor.getName());
+ IProfile1D h = aida.profile1D("deltas " + sensor.getName(),7,0,7);
+ this._delta_histos.put(sensor, h);
+ IHistogram1D h1 = aida.histogram1D("Number of iterations " + sensor.getName(),7,0,7);
+ this._delta_itercount.put(sensor, h1);
+ }
+
+ _prf_all_deltas = aida.profile1D("alldeltas",10,0,10);//,50,-20,20);
+ _prf_final_deltas = aida.profile1D("finaldeltas",10,0,10);//,50,-20,20);
+
+ plotter_iter_final = af.createPlotterFactory().create();
+ plotter_iter_final.createRegions(1,1);
+ plotter_iter_final.region(0).plot(_prf_final_deltas);
+ plotter_iter_final.region(0).style().xAxisStyle().setLabel("Final iteration");
+ plotter_iter_final.region(0).style().yAxisStyle().setLabel("<Distance to xp>");
+
+ plotter_iter_all = af.createPlotterFactory().create();
+ plotter_iter_all.createRegions(1,1);
+ plotter_iter_all.region(0).plot(_prf_all_deltas);
+ plotter_iter_all.region(0).style().xAxisStyle().setLabel("Iteration");
+ plotter_iter_all.region(0).style().yAxisStyle().setLabel("<Distance to xp>");
+
+
+ plotter_iter = af.createPlotterFactory().create();
+ plotter_iter.createRegions(5,4);
+ plotter_iter.style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter_iter.style().statisticsBoxStyle().setVisible(false);
+ plotter_itercount = af.createPlotterFactory().create();
+ plotter_itercount.createRegions(5,4);
+ int i=0;
+ for(SiSensor sensor: this.processSensors) {
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": " + i + ": adding plot to plotter for " + sensor.getName());
+ plotter_iter.region(i).plot(this._delta_histos.get(sensor));
+ plotter_iter.style().setParameter("hist2DStyle", "colorMap");
+ plotter_iter.region(i).style().xAxisStyle().setLabel("Iteration");
+ plotter_iter.region(i).style().yAxisStyle().setLabel("<Distance to xp>");
+
+ plotter_itercount.region(i).plot(this._delta_itercount.get(sensor));
+ plotter_itercount.region(i).style().xAxisStyle().setLabel("# iterations");
+ ++i;
+ }
+ if(!hideFrame) {
+ plotter_iter.show();
+ plotter_itercount.show();
+ plotter_iter_final.show();
+ plotter_iter_all.show();
+ }
+ }
+
+
+ @Override
+ public void endOfData() {
+ if (debug) {
+ for (int mod = 0; mod < 2; mod++) {
+ for (int layer = 0; layer < 10; layer++) {
+ System.out.format("mod %d, layer %d, count %d\n", mod, layer, counts[mod][layer]);
+ }
+ }
+ }
+ }
+}