Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking on MAIN | |||
DataTrackerFakeHitDriver.java | +531 | added 1.1 |
Driver that produces fake clusters (TrackerHits) on sensors from an input helix.
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]); + } + } + } + } +}
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