Commit in lcsim/src/org/lcsim/recon/tracking/digitization/sisim on MAIN
StripHitMaker.java+159-1911.5 -> 1.6
Put in settors for position errors based on size of cluster

lcsim/src/org/lcsim/recon/tracking/digitization/sisim
StripHitMaker.java 1.5 -> 1.6
diff -u -r1.5 -r1.6
--- StripHitMaker.java	8 May 2009 18:56:36 -0000	1.5
+++ StripHitMaker.java	30 Jun 2009 16:33:29 -0000	1.6
@@ -6,7 +6,6 @@
  * To change this template, choose Tools | Template Manager
  * and open the template in the editor.
  */
-
 package org.lcsim.recon.tracking.digitization.sisim;
 
 import hep.physics.matrix.SymmetricMatrix;
@@ -38,40 +37,37 @@
  *
  * @author tknelson
  */
-public class StripHitMaker implements Clusterer
-{
-    
+public class StripHitMaker implements Clusterer {
+
     private static String _NAME = "StripClusterer";
 
     // Clustering algorithm
     ClusteringAlgorithm _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<RawTrackerHit,Integer> _strip_map = new HashMap<RawTrackerHit,Integer>();
+    Map<RawTrackerHit, Integer> _strip_map = new HashMap<RawTrackerHit, Integer>();
+    double _oneClusterErr = 1 / Math.sqrt(12);
+    double _twoClusterErr = 1 / 5;
+    double _threeClusterErr = 1 / 3;
+    double _fourClusterErr = 1 / 2;
+    double _fiveClusterErr = 1;
 
-        /**
+    /**
      * Create a new instance of StripHitMaker that uses the default clustering algorithm
      *
      * @param simulation charge collection simulation
      * @param readout_chip readout chip
      */
-    public StripHitMaker(SiSensorSim simulation, ReadoutChip readout_chip)
-    {
+    public StripHitMaker(SiSensorSim simulation, ReadoutChip readout_chip) {
         this(simulation, readout_chip, new NearestNeighbor());
     }
 
@@ -92,96 +88,104 @@
         _clustering = clustering_algorithm;
     }
 
-    public void setCentralStripAveragingThreshold(int max_noaverage_nstrips)
-    {
+    public void setCentralStripAveragingThreshold(int max_noaverage_nstrips) {
         _max_noaverage_nstrips = max_noaverage_nstrips;
     }
-    
-    public void setMaxClusterSize(int max_cluster_nstrips)
-    {
+
+    public void setMaxClusterSize(int max_cluster_nstrips) {
         _max_cluster_nstrips = max_cluster_nstrips;
     }
-    
-    public String getName()
-    {
+
+    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());
+    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())
-            {
+        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());
-    	
+    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();
+        _sid_helper = (SiTrackerIdentifierHelper) sensor.getIdentifierHelper();
         _strip_map.clear();
-        
+
         // Get hits for this sensor
         IReadout ro = sensor.getReadout();
         List<RawTrackerHit> raw_hits = ro.getHits(RawTrackerHit.class);
-        
-        Map<SiSensorElectrodes,List<RawTrackerHit>> electrode_hits = new HashMap<SiSensorElectrodes,List<RawTrackerHit>>();
-        
-        for (RawTrackerHit raw_hit : raw_hits)
-        {
-            
+
+        Map<SiSensorElectrodes, List<RawTrackerHit>> electrode_hits = new HashMap<SiSensorElectrodes, List<RawTrackerHit>>();
+
+        for (RawTrackerHit raw_hit : raw_hits) {
+
             // get id and create strip map, get electrodes.
             IIdentifier id = raw_hit.getIdentifier();
-            _strip_map.put(raw_hit,_sid_helper.getElectrodeValue(id));
-            
+            _strip_map.put(raw_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)raw_hit.getDetectorElement()).getReadoutElectrodes(carrier);
-            if ( !(electrodes instanceof SiStrips) ) continue;
-            
-            if (electrode_hits.get(electrodes) == null)
-            {
-                electrode_hits.put(electrodes,new ArrayList<RawTrackerHit>());
+            SiSensorElectrodes electrodes = ((SiSensor) raw_hit.getDetectorElement()).getReadoutElectrodes(carrier);
+            if (!(electrodes instanceof SiStrips)) {
+                continue;
+            }
+
+            if (electrode_hits.get(electrodes) == null) {
+                electrode_hits.put(electrodes, new ArrayList<RawTrackerHit>());
             }
-            
+
             electrode_hits.get(electrodes).add(raw_hit);
         }
-        
-        for (Entry entry : electrode_hits.entrySet())
-        {
-            hits.addAll(makeHits(sensor, (SiStrips)entry.getKey(),(List<RawTrackerHit>)entry.getValue()));
+
+        for (Entry entry : electrode_hits.entrySet()) {
+            hits.addAll(makeHits(sensor, (SiStrips) entry.getKey(), (List<RawTrackerHit>) entry.getValue()));
         }
-        
+
         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 methods
     //=========================
-    
+
     // Make hits for an electrode
-    public List<SiTrackerHit> makeHits(SiSensor sensor, SiSensorElectrodes electrodes, List<RawTrackerHit> raw_hits)
-    {
+    public List<SiTrackerHit> makeHits(SiSensor sensor, SiSensorElectrodes electrodes, List<RawTrackerHit> raw_hits) {
 
         //  Call the clustering algorithm to make clusters
         List<List<RawTrackerHit>> cluster_list = _clustering.findClusters(electrodes, _readout_chip, raw_hits);
@@ -193,87 +197,75 @@
         for (List<RawTrackerHit> 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);
+            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);
+                ((SiSensor) electrodes.getDetectorElement()).getReadout().addHit(hit);
                 hits.add(hit);
             }
         }
 
         return hits;
     }
-              
+
     //Make the hit
-    private SiTrackerHitStrip1D makeTrackerHit(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes)
-    {
-        Hep3Vector position = getPosition(cluster,electrodes);
-        SymmetricMatrix covariance = getCovariance(cluster,electrodes);
+    private SiTrackerHitStrip1D makeTrackerHit(List<RawTrackerHit> 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);
-        
+        TrackerHitType type = new TrackerHitType(TrackerHitType.CoordinateSystem.GLOBAL, TrackerHitType.MeasurementType.STRIP_1D);
+
         SiTrackerHitStrip1D hit = new SiTrackerHitStrip1D(position, covariance, energy, time, cluster, type);
         return hit;
     }
-    
-    
-    private List<SimTrackerHit> getSimulatedHits(List<RawTrackerHit> cluster)
-    {
+
+    private List<SimTrackerHit> getSimulatedHits(List<RawTrackerHit> cluster) {
         Set<SimTrackerHit> simulated_hits = new HashSet<SimTrackerHit>();
-        for (RawTrackerHit hit : cluster)
-        {
+        for (RawTrackerHit hit : cluster) {
             simulated_hits.addAll(hit.getSimTrackerHits());
         }
         return new ArrayList<SimTrackerHit>(simulated_hits);
     }
-    
-    
-    private Hep3Vector getPosition(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes)
-    {
+
+    private Hep3Vector getPosition(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes) {
         List<Double> signals = new ArrayList<Double>();
         List<Hep3Vector> positions = new ArrayList<Hep3Vector>();
-        
-        for (RawTrackerHit hit : cluster)
-        {
+
+        for (RawTrackerHit hit : cluster) {
             signals.add(_readout_chip.decodeCharge(hit));
-            positions.add( ((SiStrips)electrodes).getStripCenter(_strip_map.get(hit)) );
+            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;
-            
+        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++)
-            {
+            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 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++)
-        {
+        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.add(position, VecOp.mult(signal, positions.get(istrip)));
         }
-        position = VecOp.mult(1/total_charge,position);
-        
+        position = VecOp.mult(1 / total_charge, position);
+
 //        double total_charge = 0;
 //        Hep3Vector position = new BasicHep3Vector(0,0,0);
 //        
@@ -286,57 +278,52 @@
 //            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());
-        
+        _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 ((SiSensor) electrodes.getDetectorElement()).getGeometry().getLocalToGlobal().transformed(position);
 //        return electrodes.getLocalToGlobal().transformed(position);
     }
-    
-    
-    private double getTime(List<RawTrackerHit> cluster)
-    {
+
+    private double getTime(List<RawTrackerHit> cluster) {
         int time_sum = 0;
         int signal_sum = 0;
-        
-        for (RawTrackerHit hit : cluster)
-        {
-            
+
+        for (RawTrackerHit hit : cluster) {
+
             int strip_number = _strip_map.get(hit);
             double signal = _readout_chip.decodeCharge(hit);
             double time = _readout_chip.decodeTime(hit);
-            
+
             time_sum += time_sum * signal;
             signal_sum += signal;
-            
+
         }
-        return (double)time_sum/(double)signal_sum;
+        return (double) time_sum / (double) signal_sum;
     }
-    
-    
-    private SymmetricMatrix getCovariance(List<RawTrackerHit> cluster, SiSensorElectrodes electrodes)
-    {
+
+    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);
-        
+        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();
@@ -351,70 +338,51 @@
 //
 //        return new SymmetricMatrix((Matrix)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 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 / Math.sqrt(12);
-        }
-        else if (cluster.size() == 2)
-        {
-            measured_resolution = sense_pitch / 5;
-        }
-        else if (cluster.size() == 3)
-        {
-            measured_resolution = sense_pitch / 3;
-        }
-        else if (cluster.size() == 4)
-        {
-            measured_resolution = sense_pitch / 2;
-        }
-        else
-        {
-            measured_resolution = sense_pitch;
+
+        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)
-    {
+
+    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(_strip_map.get(hit)) );
-        }
-        return hit_length/Math.sqrt(12);
-    }
-    
-    
-    private double getEnergy(List<RawTrackerHit> cluster)
-    {
+        for (RawTrackerHit 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<RawTrackerHit> cluster) {
         double total_charge = 0.0;
-        for (RawTrackerHit hit : cluster)
-        {
-            
+        for (RawTrackerHit hit : cluster) {
+
             int strip_number = _strip_map.get(hit);
             double signal = _readout_chip.decodeCharge(hit);
-            
+
             total_charge += signal;
         }
         return total_charge * DopedSilicon.ENERGY_EHPAIR;
     }
-    
-    
 }
CVSspam 0.2.8