hps-java/src/main/java/org/lcsim/hps/recon/tracking
diff -N HPSStripMaker.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HPSStripMaker.java 24 Apr 2012 14:07:08 -0000 1.1
@@ -0,0 +1,328 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+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.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IReadout;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.tracker.silicon.*;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.base.BaseTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.*;
+
+/**
+ *
+ * @author mgraham
+ */
+public class HPSStripMaker {
+
+ private static String _NAME = "HPSStripClusterer";
+ // Clustering algorithm
+ HPSClusteringAlgorithm _clustering;
+ // Number of strips beyond which charge is averaged on center strips
+ int _max_noaverage_nstrips = 4;
+ // Absolute maximum cluster size
+ int _max_cluster_nstrips = 10;
+ // Readout chip needed to decode hit information
+ ReadoutChip _readout_chip;
+ // Sensor simulation needed to correct for Lorentz drift
+ SiSensorSim _simulation;
+ // Identifier helper (reset once per sensor)
+ SiTrackerIdentifierHelper _sid_helper;
+ // Temporary map connecting hits to strip numbers for sake of speed (reset once per sensor)
+ Map<HPSTrackerHit, Integer> _strip_map = new HashMap<HPSTrackerHit, Integer>();
+ double _oneClusterErr = 1 / Math.sqrt(12);
+ double _twoClusterErr = 1 / 5;
+ double _threeClusterErr = 1 / 3;
+ double _fourClusterErr = 1 / 2;
+ double _fiveClusterErr = 1;
+
+ public String getName() {
+ return _NAME;
+ }
+
+
+ // Make hits for all sensors within a DetectorElement
+ public List<SiTrackerHit> makeHits(IDetectorElement detector) {
+ System.out.println("makeHits(IDetectorElement): " + detector.getName());
+ List<SiTrackerHit> hits = new ArrayList<SiTrackerHit>();
+ List<SiSensor> sensors = detector.findDescendants(SiSensor.class);
+
+ // Loop over all sensors
+ for (SiSensor sensor : sensors) {
+ if (sensor.hasStrips()) {
+ hits.addAll(makeHits(sensor));
+ }
+ }
+
+ // Return hit list
+ return hits;
+ }
+
+ // Make hits for a sensor
+ public List<SiTrackerHit> makeHits(SiSensor sensor) {
+
+ //System.out.println("makeHits: " + sensor.getName());
+
+ List<SiTrackerHit> hits = new ArrayList<SiTrackerHit>();
+
+ // Get SiTrackerIdentifierHelper for this sensor and refresh the strip map used to increase speed
+ _sid_helper = (SiTrackerIdentifierHelper) sensor.getIdentifierHelper();
+ _strip_map.clear();
+
+ // Get hits for this sensor
+ IReadout ro = sensor.getReadout();
+ List<HPSTrackerHit> hps_hits = ro.getHits(HPSTrackerHit.class);
+
+ Map<SiSensorElectrodes, List<HPSTrackerHit>> electrode_hits = new HashMap<SiSensorElectrodes, List<HPSTrackerHit>>();
+
+ for (HPSTrackerHit hps_hit : hps_hits) {
+
+ // get id and create strip map, get electrodes.
+ IIdentifier id = hps_hit.getIdentifier();
+ _strip_map.put(hps_hit, _sid_helper.getElectrodeValue(id));
+
+ // Get electrodes and check that they are strips
+ //System.out.println("proc raw hit from: " + DetectorElementStore.getInstance().find(raw_hit.getIdentifier()).get(0).getName());
+ ChargeCarrier carrier = ChargeCarrier.getCarrier(_sid_helper.getSideValue(id));
+ SiSensorElectrodes electrodes = ((SiSensor) hps_hit.getDetectorElement()).getReadoutElectrodes(carrier);
+ if (!(electrodes instanceof SiStrips)) {
+ continue;
+ }
+
+ if (electrode_hits.get(electrodes) == null) {
+ electrode_hits.put(electrodes, new ArrayList<HPSTrackerHit>());
+ }
+
+ electrode_hits.get(electrodes).add(hps_hit);
+ }
+
+ for (Map.Entry entry : electrode_hits.entrySet()) {
+ hits.addAll(makeHits(sensor, (SiStrips) entry.getKey(), (List<HPSTrackerHit>) entry.getValue()));
+ }
+
+ return hits;
+ }
+
+ public List<SiTrackerHit> makeHits(SiSensor sensor, SiSensorElectrodes electrodes,List<HPSTrackerHit> hps_hits) {
+
+
+
+ // Call the clustering algorithm to make clusters
+ List<List<HPSTrackerHit>> cluster_list = _clustering.findClusters(hps_hits);
+
+ // Create an empty list for the pixel hits to be formed from clusters
+ List<SiTrackerHit> hits = new ArrayList<SiTrackerHit>();
+
+ // Make a pixel hit from this cluster
+ for (List<HPSTrackerHit> cluster : cluster_list) {
+
+ // Make a TrackerHit from the cluster if it meets max cluster size requirement
+ if (cluster.size() <= _max_cluster_nstrips) {
+ SiTrackerHitStrip1D hit = makeTrackerHit(cluster, electrodes);
+ // Add to readout and to list of hits
+// ((SiSensor) electrodes.getDetectorElement()).getReadout().addHit(hit);
+ hits.add(hit);
+ }
+ }
+
+ return hits;
+ }
+
+ public void SetOneClusterErr(double err) {
+ _oneClusterErr = err;
+ }
+
+ public void SetTwoClusterErr(double err) {
+ _twoClusterErr = err;
+ }
+
+ public void SetThreeClusterErr(double err) {
+ _threeClusterErr = err;
+ }
+
+ public void SetFourClusterErr(double err) {
+ _fourClusterErr = err;
+ }
+
+ public void SetFiveClusterErr(double err) {
+ _fiveClusterErr = err;
+ }
+
+ private SiTrackerHitStrip1D makeTrackerHit(List<HPSTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ Hep3Vector position = getPosition(cluster, electrodes);
+ SymmetricMatrix covariance = getCovariance(cluster, electrodes);
+ double time = getTime(cluster);
+ double energy = getEnergy(cluster);
+ TrackerHitType type = new TrackerHitType(TrackerHitType.CoordinateSystem.GLOBAL, TrackerHitType.MeasurementType.STRIP_1D);
+ List<RawTrackerHit> rth_cluster = new ArrayList<RawTrackerHit>();
+ for (HPSTrackerHit bth : cluster) {
+ rth_cluster.add((RawTrackerHit) bth);
+ }
+ SiTrackerHitStrip1D hit = new SiTrackerHitStrip1D(position, covariance, energy, time, rth_cluster, type);
+ return hit;
+ }
+
+ private Hep3Vector getPosition(List<HPSTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ List<Double> signals = new ArrayList<Double>();
+ List<Hep3Vector> positions = new ArrayList<Hep3Vector>();
+
+ for (HPSTrackerHit hit : cluster) {
+ signals.add(hit.getAmp());
+ positions.add(((SiStrips) electrodes).getStripCenter(_strip_map.get(hit)));
+ }
+
+ // Average charge on central strips of longer clusters
+ if (signals.size() > _max_noaverage_nstrips) {
+ int nstrips_center = signals.size() - 2;
+
+ // collect sum of charges on center strips
+ double center_charge_sum = 0.0;
+ for (int istrip = 1; istrip < signals.size() - 1; istrip++) {
+ center_charge_sum += signals.get(istrip);
+ }
+
+ // distribute evenly on center strips
+ double center_charge_strip = center_charge_sum / nstrips_center;
+ for (int istrip = 1; istrip < signals.size() - 1; istrip++) {
+ signals.set(istrip, center_charge_strip);
+ }
+ }
+
+ double total_charge = 0;
+ Hep3Vector position = new BasicHep3Vector(0, 0, 0);
+
+ for (int istrip = 0; istrip < signals.size(); istrip++) {
+ double signal = signals.get(istrip);
+
+ total_charge += signal;
+ position = VecOp.add(position, VecOp.mult(signal, positions.get(istrip)));
+ }
+ position = VecOp.mult(1 / total_charge, position);
+
+// double total_charge = 0;
+// Hep3Vector position = new BasicHep3Vector(0,0,0);
+//
+// for (RawTrackerHit hit : cluster)
+// {
+// int strip_number = _strip_map.get(hit);
+// double signal = _readout_chip.decodeCharge(hit);
+//
+// total_charge += signal;
+// position = VecOp.add(position,VecOp.mult(signal,((SiStrips)electrodes).getStripCenter(strip_number)));
+// }
+// position = VecOp.mult(1/total_charge,position);
+
+ // Put position in sensor coordinates
+ electrodes.getParentToLocal().inverse().transform(position);
+
+// System.out.println("Position \n"+position);
+
+ // Swim position back through lorentz drift direction to midpoint between bias surfaces
+ _simulation.setSensor((SiSensor) electrodes.getDetectorElement());
+ _simulation.lorentzCorrect(position, electrodes.getChargeCarrier());
+
+// System.out.println("Lorentz corrected position \n"+position);
+
+ // return position in global coordinates
+ return ((SiSensor) electrodes.getDetectorElement()).getGeometry().getLocalToGlobal().transformed(position);
+// return electrodes.getLocalToGlobal().transformed(position);
+ }
+
+ private double getTime(List<HPSTrackerHit> cluster) {
+ int time_sum = 0;
+ int signal_sum = 0;
+
+ for (HPSTrackerHit hit : cluster) {
+
+ double signal = hit.getAmp();
+ double time = hit.getTime();
+
+ time_sum += time * signal;
+ signal_sum += signal;
+
+ }
+ return (double) time_sum / (double) signal_sum;
+ }
+
+ private SymmetricMatrix getCovariance(List<HPSTrackerHit> 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;
+
+// BasicHep3Matrix rotation_matrix = (BasicHep3Matrix)electrodes.getLocalToGlobal().getRotation().getRotationMatrix();
+// BasicHep3Matrix rotation_matrix_transposed = new BasicHep3Matrix(rotation_matrix);
+// rotation_matrix_transposed.transpose();
+//
+//// System.out.println("Rotation matrix: \n"+rotation_matrix);
+//// System.out.println("Rotation matrix transposed: \n"+rotation_matrix_transposed);
+//// System.out.println("Local covariance matrix: \n"+covariance);
+//
+// BasicHep3Matrix covariance_global = (BasicHep3Matrix)VecOp.mult(rotation_matrix,VecOp.mult(covariance,rotation_matrix_transposed));
+//
+//// System.out.println("Global covariance matrix: \n"+covariance_global);
+//
+// return new SymmetricMatrix((Matrix)covariance_global);
+ }
+
+ private double getMeasuredResolution(List<HPSTrackerHit> 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
+
+ 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<HPSTrackerHit> cluster, SiSensorElectrodes electrodes) {
+ // Get length of longest strip in hit
+ double hit_length = 0;
+ for (HPSTrackerHit hit : cluster) {
+ hit_length = Math.max(hit_length, ((SiStrips) electrodes).getStripLength(_strip_map.get(hit)));
+ }
+ return hit_length / Math.sqrt(12);
+ }
+
+ private double getEnergy(List<HPSTrackerHit> cluster) {
+ double total_charge = 0.0;
+ for (HPSTrackerHit hit : cluster) {
+ double signal = hit.getAmp();
+ total_charge += signal;
+ }
+ return total_charge * DopedSilicon.ENERGY_EHPAIR;
+ }
+}