Iss135: create SiliconResolutionModel for use in strip making

tracking/src/main/java/org/hps/recon/tracking/DataTrackerHitDriver.java

@@ -49,6 +49,8 @@ public class DataTrackerHitDriver extends Driver {
     private double threeClusterErr = clusterErrorMultiplier / 3.0;
     private double fourClusterErr = clusterErrorMultiplier / 2.0;
     private double fiveClusterErr = clusterErrorMultiplier / 1.0;
+    // weight the hits in a cluster by charge? (if not, all hits have equal weight)
+    private boolean useWeights = true;
     // Various data lists required by digitization.
     private List<String> processPaths = new ArrayList<String>();
     private List<IDetectorElement> processDEs = new ArrayList<IDetectorElement>();
@@ -127,6 +129,10 @@ public void setFiveClusterErr(double fiveClusterErr) {
         this.fiveClusterErr = fiveClusterErr;
     }
 
+    public void setUseWeights(boolean useWeights){
+        this.useWeights = useWeights;
+    }
+
     /**
      * Creates a new instance of TrackerHitDriver.
      */
@@ -180,11 +186,18 @@ public void detectorChanged(Detector detector) {
         stripClusterer.setCentralStripAveragingThreshold(clusterCentralStripAveragingThreshold);
 
         // Set the cluster errors.
-        stripClusterer.SetOneClusterErr(oneClusterErr);
-        stripClusterer.SetTwoClusterErr(twoClusterErr);
-        stripClusterer.SetThreeClusterErr(threeClusterErr);
-        stripClusterer.SetFourClusterErr(fourClusterErr);
-        stripClusterer.SetFiveClusterErr(fiveClusterErr);
+
+        DefaultSiliconResolutionModel model = new DefaultSiliconResolutionModel();
+
+        model.setOneClusterErr(oneClusterErr);
+        model.setTwoClusterErr(twoClusterErr);
+        model.setThreeClusterErr(threeClusterErr);
+        model.setFourClusterErr(fourClusterErr);
+        model.setFiveClusterErr(fiveClusterErr);
+        model.setUseWeights(useWeights);
+
+        stripClusterer.setResolutionModel(model);
+
 
         // Set the detector to process.
         processPaths.add(subdetectorName);

tracking/src/main/java/org/hps/recon/tracking/DefaultSiliconResolutionModel.java

@@ -0,0 +1,108 @@
+package org.hps.recon.tracking;
+
+import java.util.List;
+import java.util.Map;
+
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.detector.tracker.silicon.SiStrips;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+public class DefaultSiliconResolutionModel implements SiliconResolutionModel{
+
+    @Override
+    public double getMeasuredResolution(List<FittedRawTrackerHit> cluster, SiSensorElectrodes electrodes) 
+
+    {
+        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;
+    }
+
+    public double getUnmeasuredResolution(List<FittedRawTrackerHit> cluster, SiSensorElectrodes electrodes, Map<FittedRawTrackerHit, Integer> strip_map) {
+        // Get length of longest strip in hit
+        double hit_length = 0;
+        for (FittedRawTrackerHit hit : cluster) {
+            hit_length = Math.max(hit_length, ((SiStrips) electrodes).getStripLength(strip_map.get(hit)));
+        }
+        return hit_length / Math.sqrt(12);
+    }
+
+
+    //perhaps the best values for these are .19, .12 and .2?    
+    private double _oneClusterErr = 1 / Math.sqrt(12);
+    private double _twoClusterErr = 1 / 5.;
+    private double _threeClusterErr = 1 / 3.;
+    private double _fourClusterErr = 1 / 2.;
+    private double _fiveClusterErr = 1;
+
+    private boolean _useWeights = true;
+
+    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;
+    }
+
+    public void setUseWeights(boolean useWeights){
+        _useWeights = useWeights;
+    }
+
+    @Override
+    public Hep3Vector weightedAveragePosition(List<Double> signals, List<Hep3Vector> positions) {
+        double total_weight = 0;
+        Hep3Vector position = new BasicHep3Vector(0, 0, 0);
+
+        for (int istrip = 0; istrip < signals.size(); istrip++) {
+            double signal = signals.get(istrip);
+
+            double weight = _useWeights ? signal : 1;
+            total_weight += weight;
+            position = VecOp.add(position, VecOp.mult(weight, positions.get(istrip)));
+            /*if (_debug) {
+                System.out.println(this.getClass().getSimpleName() + "strip " + istrip + ": signal " + signal + " position " + positions.get(istrip) + " -> total_position " + position.toString() + " ( total charge " + total_charge + ")");
+            }*/
+
+        }
+
+        return VecOp.mult(1 / total_weight, position);
+    }
+
+}

tracking/src/main/java/org/hps/recon/tracking/SiliconResolutionModel.java

@@ -0,0 +1,14 @@
+package org.hps.recon.tracking;
+
+import java.util.List;
+import java.util.Map;
+
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+
+import hep.physics.vec.Hep3Vector;
+
+public interface SiliconResolutionModel {
+    public double getMeasuredResolution(List<FittedRawTrackerHit> cluster, SiSensorElectrodes electrodes);
+    public double getUnmeasuredResolution(List<FittedRawTrackerHit> cluster, SiSensorElectrodes electrodes, Map<FittedRawTrackerHit, Integer> strip_map);
+    public Hep3Vector weightedAveragePosition(List<Double> signals, List<Hep3Vector> positions);
+}

tracking/src/main/java/org/hps/recon/tracking/StripMaker.java

@@ -1,9 +1,7 @@
 package org.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;
@@ -46,14 +44,15 @@ public class StripMaker {
     SiTrackerIdentifierHelper _sid_helper;
     // Temporary map connecting hits to strip numbers for sake of speed (reset once per sensor)
     Map<FittedRawTrackerHit, Integer> _strip_map = new HashMap<FittedRawTrackerHit, Integer>();
-    double _oneClusterErr = 1 / Math.sqrt(12);
-    double _twoClusterErr = 1 / 5;
-    double _threeClusterErr = 1 / 3;
-    double _fourClusterErr = 1 / 2;
-    double _fiveClusterErr = 1;
+
 
     boolean _debug = false;
+    private SiliconResolutionModel _res_model = new DefaultSiliconResolutionModel();
 
+    public void setResolutionModel(SiliconResolutionModel model){
+        _res_model = model;
+    }
+
     public StripMaker(ClusteringAlgorithm algo) {
         _clustering = algo;
     }
@@ -154,25 +153,7 @@ public List<SiTrackerHit> makeHits(SiSensor sensor, SiSensorElectrodes electrode
         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;
-    }
+
 
     public void setCentralStripAveragingThreshold(int max_noaverage_nstrips) {
         _max_noaverage_nstrips = max_noaverage_nstrips;
@@ -242,20 +223,9 @@ private Hep3Vector getPosition(List<FittedRawTrackerHit> cluster, SiSensorElectr
             System.out.println(this.getClass().getSimpleName() + " Calculate charge weighted mean for " + signals.size() + " signals");
         }
 
-        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)));
-            if (_debug) {
-                System.out.println(this.getClass().getSimpleName() + "strip " + istrip + ": signal " + signal + " position " + positions.get(istrip) + " -> total_position " + position.toString() + " ( total charge " + total_charge + ")");
-            }
-
-        }
-        position = VecOp.mult(1 / total_charge, position);
+
+        Hep3Vector position = _res_model.weightedAveragePosition(signals, positions);
+
         if (_debug) {
             System.out.println(this.getClass().getSimpleName() + " charge weighted position " + position.toString() + " (before trans)");
         }
@@ -303,8 +273,8 @@ private double getTime(List<FittedRawTrackerHit> cluster) {
 
     private SymmetricMatrix getCovariance(List<FittedRawTrackerHit> 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(0, 0, Math.pow(_res_model.getMeasuredResolution(cluster, electrodes), 2));
+        covariance.setElement(1, 1, Math.pow(_res_model.getUnmeasuredResolution(cluster, electrodes, _strip_map), 2));
         covariance.setElement(2, 2, 0.0);
 
         SymmetricMatrix covariance_global = electrodes.getLocalToGlobal().transformed(covariance);
@@ -329,57 +299,7 @@ private SymmetricMatrix getCovariance(List<FittedRawTrackerHit> cluster, SiSenso
         // return new SymmetricMatrix((Matrix)covariance_global);
     }
 
-    private double getMeasuredResolution(List<FittedRawTrackerHit> 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<FittedRawTrackerHit> cluster, SiSensorElectrodes electrodes) {
-        // Get length of longest strip in hit
-        double hit_length = 0;
-        for (FittedRawTrackerHit 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<FittedRawTrackerHit> cluster) {
         double total_charge = 0.0;