Print

Print


Commit in GeomConverter/src/org/lcsim/detector/tracker/silicon on MAIN
DopedSilicon.java+1-11.2 -> 1.3
ErrorEllipse2D.java+24-71.1 -> 1.2
ChargeCarrier.java+7-31.1 -> 1.2
SiSensor.java+269-1351.2 -> 1.3
SiSensorElectrodes.java+1-11.2 -> 1.3
SiStrips.java+177-1161.3 -> 1.4
+479-263
6 modified files
Updates and bugfixes to silicon strip simulation: lots of debug code.

GeomConverter/src/org/lcsim/detector/tracker/silicon
DopedSilicon.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- DopedSilicon.java	19 Apr 2007 22:37:41 -0000	1.2
+++ DopedSilicon.java	7 May 2007 21:22:49 -0000	1.3
@@ -30,7 +30,7 @@
     // Member
     private double _temperature = 293.0; // * kelvin;
     private double _doping_concentration = 6.0E+11; // / cm3; 
-    private EnumMap<ChargeCarrier,Double> _carrier_concentration = new EnumMap<ChargeCarrier,Double>(ChargeCarrier.class);    
+    private EnumMap<ChargeCarrier, Double> _carrier_concentration = new EnumMap<ChargeCarrier,Double>(ChargeCarrier.class);    
     
     // Constructors
     //=============

GeomConverter/src/org/lcsim/detector/tracker/silicon
ErrorEllipse2D.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- ErrorEllipse2D.java	19 Apr 2007 22:37:41 -0000	1.1
+++ ErrorEllipse2D.java	7 May 2007 21:22:49 -0000	1.2
@@ -45,22 +45,39 @@
         return this;
     }
     
-    public double erf1D(double integration_limit, double axis_angle)
+    public double erfc1D(double integration_limit, double axis_angle)
     {
 //        double angle_diff = axis_angle - _major_axis_angle;
 //        double sigma = Math.sqrt( Math.pow(_major_axis*Math.cos(angle_diff),2) + Math.pow(_minor_axis*Math.sin(angle_diff),2) );
         
-        double erf = 0.0;
+        double erfc=0.0;
         
-        try
+        double normalized_integration_limit = integration_limit/sigma1D(axis_angle);
+        if (normalized_integration_limit < -5.0)
         {
-            erf = Erf.erf( integration_limit/sigma1D(axis_angle) );
+            erfc = 1.0;
         }
-        catch (MathException no_convergence)
+        else if (normalized_integration_limit > 5.0)
         {
-            System.out.println("erf fails to converge!!");
+            erfc = 0.0;
         }
-        return erf;        
+        else
+        {
+            try
+                    
+            {
+                erfc = (1.0-Erf.erf( normalized_integration_limit/Math.sqrt(2.0)))/2.0;
+            }
+            catch (MathException no_convergence)
+            {
+                System.out.println("Warning: erf fails to converge!! ");
+                System.out.println("    integration limit: "+integration_limit);
+                System.out.println("    sigma: "+sigma1D(axis_angle));
+                System.out.println("    limit/sigma: "+normalized_integration_limit);
+            }
+        }
+
+        return erfc;        
     }
     
     public double sigma1D(double axis_angle)

GeomConverter/src/org/lcsim/detector/tracker/silicon
ChargeCarrier.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- ChargeCarrier.java	10 Apr 2007 19:27:44 -0000	1.1
+++ ChargeCarrier.java	7 May 2007 21:22:49 -0000	1.2
@@ -16,9 +16,10 @@
 
 public enum ChargeCarrier
 {
-    ELECTRON(1268.0,-2.33,92.0,-0.57,1.3E+17,2.4,0.91,-0.146),
-    HOLE(406.9,-2.23,54.3,-0.57,2.35E+17,2.4,0.88,-0.146);
+    ELECTRON(-1,1268.0,-2.33,92.0,-0.57,1.3E+17,2.4,0.91,-0.146),
+    HOLE(1,406.9,-2.23,54.3,-0.57,2.35E+17,2.4,0.88,-0.146);
     
+    private final int _charge;
     private final double _mu_0_factor;
     private final double _mu_0_exponent;
     private final double _mu_min_factor;
@@ -28,9 +29,10 @@
     private final double _alpha_factor;
     private final double _alpha_exponent;
     
-    ChargeCarrier(double mu_0_factor, double mu_0_exponent, double mu_min_factor, double mu_min_exponent,
+    ChargeCarrier(int charge, double mu_0_factor, double mu_0_exponent, double mu_min_factor, double mu_min_exponent,
             double N_ref_factor, double N_ref_exponent, double alpha_factor, double alpha_exponent)
     {
+        _charge = charge;
         _mu_0_factor = mu_0_factor;
         _mu_0_exponent = mu_0_exponent;
         _mu_min_factor = mu_min_factor;
@@ -42,6 +44,8 @@
     }
     
     // Methods
+    int charge()
+    {return _charge;}
     double mu0(double temperature)
     {return _mu_0_factor * Math.pow( (temperature/300.0), _mu_0_exponent);}
     double muMin(double temperature)

GeomConverter/src/org/lcsim/detector/tracker/silicon
SiSensor.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- SiSensor.java	20 Apr 2007 06:31:27 -0000	1.2
+++ SiSensor.java	7 May 2007 21:22:49 -0000	1.3
@@ -14,7 +14,7 @@
 import org.lcsim.detector.IPhysicalVolumePath;
 import org.lcsim.detector.solids.Box;
 import org.lcsim.detector.Rotation3D;
-import org.lcsim.geometry.Detector;
+import org.lcsim.detector.converter.compact.DeDetector;
 //import static org.lcsim.units.clhep.SystemOfUnits.*;
 //import static org.lcsim.units.clhep.PhysicalConstants.*;        
 
@@ -30,19 +30,22 @@
  */
 public class SiSensor extends DetectorElement {
     
+//    _pside_direction?
+ //   center at zero!!!
+    
     // Enumerated types
     //=================
-    public enum Orientation { PINSIDE, POUTSIDE };
+    public enum Orientation { PSIDE_NEGATIVE_Z, PSIDE_POSITIVE_Z };
     
     // Fields
     //=======
     // Static defaults - actual values user modifiable
-    private static Orientation _ORIENTATION_DEFAULT = Orientation.POUTSIDE;
+    private static Orientation _ORIENTATION_DEFAULT = Orientation.PSIDE_POSITIVE_Z;
     private static double _DEPLETION_VOLTAGE_DEFAULT = 100;// * volt;
     private static double _BIAS_VOLTAGE_DEFAULT = 110;// * volt;
     
     // Static parameters - not intended to be user modifiable
-    private static double _DEPOSITION_GRANULARITY = 0.05; // 5% of pitch
+    private static double _DEPOSITION_GRANULARITY = 0.10; // 10% of pitch or depleted thickness
     
     // primary properties
     //-------------------
@@ -51,8 +54,8 @@
     private int _sensorid;
     
     // electrodes, electrode angles and orientation of sensor
-    private Map<ChargeCarrier,SiSensorElectrodes> _electrodes = new EnumMap<ChargeCarrier,SiSensorElectrodes>(ChargeCarrier.class);
-    private Map<ChargeCarrier,Double> _electrode_angles = new EnumMap<ChargeCarrier,Double>(ChargeCarrier.class);
+    private Map<ChargeCarrier, SiSensorElectrodes> _electrodes = new EnumMap<ChargeCarrier,SiSensorElectrodes>(ChargeCarrier.class);
+    private Map<ChargeCarrier, Double> _electrode_angles = new EnumMap<ChargeCarrier,Double>(ChargeCarrier.class);
     private Orientation _orientation;
 
     // bulk - propoerties of the bulk
@@ -65,14 +68,11 @@
     
     // derived properties
     //-------------------
-    // measured coordinate in local coordinates
-//    private EnumMap<ChargeCarrier,Hep3Vector> _measured_coordinate = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
-    private EnumMap<ChargeCarrier,Hep3Vector[]> _measured_coordinates = new EnumMap<ChargeCarrier,Hep3Vector[]>(ChargeCarrier.class);
-    // strip direction in local coordinates
-//    private EnumMap<ChargeCarrier,Hep3Vector> _strip_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
+    // measured coordinates in local coordinates
+    private EnumMap<ChargeCarrier, Hep3Vector[]> _measured_coordinates = new EnumMap<ChargeCarrier,Hep3Vector[]>(ChargeCarrier.class);
     // direction of Lorentz drift in local coordinates
-    private EnumMap<ChargeCarrier,Hep3Vector> _drift_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
-    // list of charge depostions to be distributed onto the electrodes
+    private EnumMap<ChargeCarrier, Hep3Vector> _drift_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
+    // list of charge depositions to be distributed onto the electrodes
     private List<TrackSegment> _track_list = new ArrayList<TrackSegment>();
     
     // Constructors
@@ -89,8 +89,9 @@
         super(name,parent,support);
         setSensorID(sensorid);
         setBulk(new DopedSilicon());
-        setDepletionVoltage(_DEPLETION_VOLTAGE_DEFAULT);
-        setBiasVoltage(_BIAS_VOLTAGE_DEFAULT);
+        this.setOrientation(SiSensor._ORIENTATION_DEFAULT);
+        setDepletionVoltage(SiSensor._DEPLETION_VOLTAGE_DEFAULT);
+        setBiasVoltage(SiSensor._BIAS_VOLTAGE_DEFAULT);
     }
     
     // Accessors
@@ -181,39 +182,66 @@
     
     public Hep3Vector getBField(Hep3Vector local_position)
     {
+        
+//        System.out.println("Beginning getBField");                
+        
         IDetectorElement ancestor = this.getParent();
-        while (!(ancestor instanceof Detector) && !(ancestor==null))
+        while (!(ancestor instanceof DeDetector) && !(ancestor==null))
         {
             ancestor = ancestor.getParent();
         }
-        
+        if (ancestor == null) System.out.println("WARNING: SiSensor.getBField CANNOT FIND DETECTOR!!");                
+
         Hep3Vector global_position = getGeometry().getLocalToGlobal().transformed(local_position);
-        Hep3Vector field_global = ((Detector)ancestor).getFieldMap().getField(global_position);
+        Hep3Vector field_global = ((DeDetector)ancestor).getBField(global_position);
+
+//        return new BasicHep3Vector(0.0,-5.0,0.0);
         
         // FIXME - This is silly!!!!!
         return VecOp.mult(getGeometry().getGlobalToLocal().getRotation().getRotationMatrix(),field_global);
+        
     }
         
+    
+    
+    
+    
+    
+    
     // Operators
     //==========
-    private void initialize()
+    public void initialize()
     {
+//        System.out.println("Beginning initialize");                
+
+        // Cache thickness of bulk
+        _thickness = 2.0*((Box)getGeometry().getLogicalVolume().getSolid()).getZHalfLength();
+        
         // Store various important directions
         for (ChargeCarrier carrier : ChargeCarrier.values())
         {
             if (hasElectrodesOnSide(carrier)) {
 
+//                System.out.println("init: Carrier: "+carrier);                        
+                
                 // cache drift direction for electrodes on each side
                 _drift_direction.put(carrier,driftDirection(carrier, new BasicHep3Vector(0.0,0.0,0.0)));
-
+//                System.out.println("init: Drift direction: "+_drift_direction);                        
+                
                 // cache coordinates measured by the pattern of electrodes on each side
-                double electrode_angle = _electrode_angles.get(carrier);
-                int naxes = _electrodes.get(carrier).getNAxes();                
+                double electrode_angle = getElectrodeAngle(carrier); 
+//                System.out.println("Electrode angle: "+electrode_angle);                                        
+               
+                int naxes = _electrodes.get(carrier).getNAxes();
+//                System.out.println("# axes "+naxes);                                        
+
                 Hep3Vector[] measured_coordinates = new Hep3Vector[naxes];
                 
-                for (int iaxis = 1; iaxis < naxes ; iaxis++)
+                for (int iaxis = 0; iaxis < naxes ; iaxis++)
                 {
+//                    System.out.println("Axis number: "+iaxis);                                                            
                     measured_coordinates[iaxis] = measuredCoordinate(electrode_angle + iaxis*Math.PI/naxes);
+//                    System.out.println("Measured coordinate: "+measured_coordinates[iaxis]);                                        
                 }
                 _measured_coordinates.put(carrier,measured_coordinates);
                 
@@ -221,22 +249,27 @@
 //                _strip_direction.put(carrier,stripDirection(electrode_angle));
             }
         }
-        // Cache thickness of bulk
-        _thickness = 2.0*((Box)getGeometry().getLogicalVolume().getSolid()).getZHalfLength();
+
     }
     
     private double zOfSide(ChargeCarrier carrier)
-    {
-        if ( (carrier == ChargeCarrier.HOLE) == (_orientation == Orientation.POUTSIDE) ) return _thickness;
-        else return 0;
+    {   
+//        System.out.println("_thickness: "+_thickness);
+        if ( (carrier == ChargeCarrier.HOLE) == (_orientation == Orientation.PSIDE_POSITIVE_Z) ) return _thickness/2.0;
+        else return -_thickness/2.0;
     }
     
     private double distanceFromSide(Hep3Vector point, ChargeCarrier carrier)
     {
-        return point.z() - zOfSide(carrier);
+//        System.out.println("Beginning distanceFromSide");        
+        double distance =  Math.abs(point.z() - zOfSide(carrier));
+//        System.out.println("point.z(): "+point.z());
+//        System.out.println("side z: "+zOfSide(carrier));        
+//        System.out.println("Distance from side is: "+distance);
+        return distance;
     }
     
-    private boolean hasElectrodesOnSide(ChargeCarrier carrier)
+    public boolean hasElectrodesOnSide(ChargeCarrier carrier)
     {
         if(_electrodes.get(carrier) == null) return false;
         else return true;
@@ -244,12 +277,14 @@
     
     private Hep3Vector driftVector(Hep3Vector origin, ChargeCarrier carrier)
     {
+//        System.out.println("Beginning driftVector");        
         double drift_vector_scale = distanceFromSide(origin,carrier)/_drift_direction.get(carrier).z();
         return VecOp.mult(drift_vector_scale,_drift_direction.get(carrier));
     }
     
     private Hep3Vector driftDestination(Hep3Vector origin, ChargeCarrier carrier)
     {
+//        System.out.println("Beginning driftDestination");        
         return VecOp.add(origin,driftVector(origin, carrier));
     }
     
@@ -261,6 +296,9 @@
         double sum_V = _bias_voltage + _depletion_voltage;
         double common_factor = 2.0*distanceFromSide(point,carrier)*_depletion_voltage/_thickness;
 
+//        System.out.println("sum_V: "+sum_V);  
+//        System.out.println("common_factor: "+common_factor);  
+        
         // Calculate charge spreading without magnetic field
         double sigmasq = _bulk.K_BOLTZMANN * _bulk.getTemperature() * _thickness*_thickness / _depletion_voltage;
         if (_bulk.isNtype() == (carrier==ChargeCarrier.HOLE))
@@ -274,128 +312,109 @@
         
         double sigma = Math.sqrt(sigmasq);
         
+//        System.out.println("Sigma: "+sigma);  
+        
         // Corrections for magnetic field -- this is an approximation, may have to be done better for high fields
         double cos_theta_lorentz = VecOp.cosTheta(_drift_direction.get(carrier));
         double phi_lorentz = VecOp.phi(_drift_direction.get(carrier));
-        double phi_electrode = _electrode_angles.get(carrier);
+        double phi_electrode = getElectrodeAngle(carrier);
         
         double minor_axis = sigma*(1.0/cos_theta_lorentz); // drift time correction
         double major_axis = minor_axis*(1.0/cos_theta_lorentz); // + drift angle correction
         double phi_ellipse = -phi_electrode; // orientation of ellipse, relative to electrode coordinates
         
+//        System.out.println("Minor axis: "+minor_axis);  
+//        System.out.println("Major axis: "+major_axis);  
+                
         // Create error ellipse
         return new ErrorEllipse2D(major_axis, minor_axis, phi_ellipse);
         
     }
-    
-//    private double diffusionSigma(Hep3Vector point, ChargeCarrier carrier)
-//    {
-//        
-//        // Common factors
-//        double difference_V = _bias_voltage - _depletion_voltage;
-//        double sum_V = _bias_voltage + _depletion_voltage;
-//        double common_factor = 2.0*distanceFromSide(point,carrier)*_depletion_voltage/_thickness;
-//
-//        // Calculate charge spreading without magnetic field
-//        double sigmasq = k_Boltzmann * _bulk.getTemperature() * _thickness*_thickness / _depletion_voltage;
-//        if (_bulk.isNtype() == (carrier==ChargeCarrier.HOLE))
-//        {
-//            sigmasq *= Math.log( sum_V / (sum_V - common_factor));
-//        }
-//        else
-//        {
-//            sigmasq *= Math.log( (difference_V + common_factor) / difference_V );
-//        }
-//        
-//        double sigma = Math.sqrt(sigmasq);
-//        
-//        // Corrections for magnetic field -- this is an approximation, may have to be done better for high fields
-//        double cos_theta_lorentz_sq = Math.pow(VecOp.cosTheta(_drift_direction.get(carrier)),2);
-//        double phi_lorentz = VecOp.phi(_drift_direction.get(carrier));
-//        double phi_measured = VecOp.phi(_measured_coordinate.get(carrier));
-//        double cos_phi_diff_sq = Math.pow(Math.cos(phi_measured - phi_lorentz),2);
-//        
-//        sigma *= (1.0/cos_theta_lorentz_sq) *
-//                Math.sqrt(cos_theta_lorentz_sq + cos_phi_diff_sq - cos_theta_lorentz_sq*cos_phi_diff_sq);
-//        
-//        return sigma;
-//        
-//    }
 
+    
     private Hep3Vector measuredCoordinate(double electrode_angle)
     {
         return new BasicHep3Vector(Math.cos(electrode_angle),Math.sin(electrode_angle),0.0);
     }
     
-//    private Hep3Vector stripDirection(double strip_angle)
-//    {
-//        return measuredCoordinate(strip_angle + (Math.PI/2.0)) ;
-//    }
-    
+
     private Hep3Vector driftDirection(ChargeCarrier carrier, Hep3Vector local_position)
     {
-        double carrier_mobility = _bulk.mobility(carrier);        
-        Hep3Vector b_field = getBField(local_position);
+//        System.out.println("Beginning driftDirection");  
+        
+//        System.out.println("Position: "+local_position);
+        
+        double carrier_mobility = _bulk.mobility(carrier);
+//        System.out.println("Carrier: "+carrier);
+        
+        Hep3Vector b_field = this.getBField(local_position);
+//        System.out.println("B field: "+b_field);
+        
+        Hep3Vector e_field = this.electricField(local_position);
+//        System.out.println("E field: "+e_field);
+        
+        double tan_lorentz = _bulk.tanLorentzAngle(b_field.magnitude(), carrier);
+        
+//        System.out.println("Tan lorentz: "+tan_lorentz);                
+        
+        Hep3Vector drift_direction = VecOp.unit(VecOp.add(
+                e_field,VecOp.mult(tan_lorentz, VecOp.cross(e_field,VecOp.unit(b_field)))
+                ));
+        
+//        System.out.println("Drift direction: "+drift_direction);                
         
-        // Use magnetic field in plane of silicon to get total Lorentz angle
-        Hep3Vector b_planar = new BasicHep3Vector(b_field.x(), b_field.y(),0.0);
-        double bplanar_mag = b_planar.magnitude();
-        double tan_lorentz = _bulk.tanLorentzAngle(bplanar_mag, carrier);
-
-        // Drift direction in plane of silicon is direction of magnetic field in plane of silicon - PI/4
-        Hep3Vector drift_direction = 
-                VecOp.unit(new BasicHep3Vector(b_field.y(),-b_field.x(),bplanar_mag/tan_lorentz));
-
         return drift_direction;
         
     }
     
+    public Hep3Vector electricField(Hep3Vector position)
+    {
+        Hep3Vector electric_field_direction = (this._orientation == Orientation.PSIDE_POSITIVE_Z) ?
+            new BasicHep3Vector(0.0,0.0,1.0) : new BasicHep3Vector(0.0,0.0,-1.0);
+        
+//        System.out.println("E field direction: "+electric_field_direction);
+        
+        double electric_field_magnitude = (_bias_voltage-_depletion_voltage)/_thickness + 
+                (2.0*_depletion_voltage)/_thickness * (1.0 - this.distanceFromSide(position,ChargeCarrier.ELECTRON));
+        
+//        System.out.println("E field magnitude: "+electric_field_magnitude);
+        
+        return VecOp.mult(electric_field_magnitude,electric_field_direction);
+        
+    }
+    
     public void clearElectrodes()
     {
         for (ChargeCarrier carrier : ChargeCarrier.values())
         {
             if (hasElectrodesOnSide(carrier)) _electrodes.get(carrier).clear();
-        }
+        }        
     }
-    
-// Delta-ray code was deprecated in favor of letting GEANT do the work    
-//
-//    public void generateDeltaRays()
-//    {
-//        for (TrackSegment track : _track_list)
-//        {
-//            // Uncommitted standalone code exists to do this, which...
-//            // retrieves track from _track_list and calculates delta ray production
-//            // modifies original track in _track_list
-//            // adds delta rays to _track_list
-//        }
-//        return;
-//    }
-
-//    private int nSegments(TrackSegment track, ChargeCarrier carrier, double deposition_granularity)
-//    {
-//        // Decide how to cut track into pieces as a fraction of strip pitch
-//        if (!hasElectrodesOnSide(carrier)) return 0;
-//        Hep3Vector deposition_line = VecOp.sub( driftDestination(track.getP2(),carrier),
-//                driftDestination(track.getP1(),carrier) );
-//
-//        double projected_deposition_length = VecOp.dot(deposition_line,_measured_coordinate.get(carrier));
-//        
-//        return (int)Math.ceil(projected_deposition_length/(deposition_granularity*_electrodes.get(carrier).getPitch()));
-//    }
+   
 
     private int nSegments(TrackSegment track, ChargeCarrier carrier, double deposition_granularity)
     {
         // Decide how to cut track into pieces as a fraction of strip pitch        
         int nsegments = 0;
         if (!hasElectrodesOnSide(carrier)) return nsegments;
+        
+//        System.out.println("Track P1: " + track.getP1());        
+//        System.out.println("Track P2: " + track.getP2());        
+//        System.out.println("Drift Destination of P1: " + driftDestination(track.getP1(),carrier));
+//        System.out.println("Drift Destination of P2: " + driftDestination(track.getP2(),carrier));        
+        
+        nsegments = (int)Math.ceil(track.getVector().z()/(_thickness*deposition_granularity));
+        
         Hep3Vector deposition_line = VecOp.sub( driftDestination(track.getP2(),carrier),
                 driftDestination(track.getP1(),carrier) );
        
         int naxes = _electrodes.get(carrier).getNAxes();                        
         for (int iaxis = 0; iaxis < naxes; iaxis++)
         {
-            double projected_deposition_length = VecOp.dot(deposition_line,_measured_coordinates.get(carrier)[iaxis]);
+            double projected_deposition_length = Math.abs(VecOp.dot(deposition_line,_measured_coordinates.get(carrier)[iaxis]));
+
+//            System.out.println("Projected deposition Length: " + projected_deposition_length);
+            
             int required_segments = (int)Math.ceil(projected_deposition_length/(deposition_granularity*_electrodes.get(carrier).getPitch(iaxis)));
             nsegments = Math.max(nsegments,required_segments);
         }   
@@ -403,62 +422,56 @@
     }
     
     
-//    private void depositCharge(double charge, Hep3Vector origin, ChargeCarrier carrier)
-//    {
-//        if (!hasStripsOnSide(carrier)) return;
-//        
-//        // find center of charge deposition
-//        double drift_destination = VecOp.dot( driftDestination(origin,carrier),
-//                _measured_coordinate.get(carrier) );
-//        double diffusion_sigma = diffusionSigma(origin,carrier);
-//        
-//        _electrodes.get(carrier).depositCharge(new BasicHep3Vector(drift_destination,0.0,0.0),charge,diffusion_sigma);
-//      
-//    }
-    
     public void depositCharge()
     {
-        
+                
         for (TrackSegment track : _track_list)
-        {
+        {  
+            
+//            System.out.println("New TrackSegment... ");
             
             // Decide how to cut track into pieces - use 5% of pitch
             int nsegments = 0;
+//            System.out.println("Number of charge carriers: " + ChargeCarrier.values().length);                        
             for (ChargeCarrier carrier : ChargeCarrier.values())
             {
+//                System.out.println("Charge carrier: " + carrier);
+//                System.out.println("Has strips on side: "+hasElectrodesOnSide(carrier));
                 if (!hasElectrodesOnSide(carrier)) continue;
+                
                 nsegments = Math.max(nsegments,nSegments(track,carrier, _DEPOSITION_GRANULARITY));
             }
             
+//            System.out.println("Number of subsegments: " + nsegments);
+            
             // Set up segments
             double segment_length = track.getLength()/nsegments;
             double segment_charge = track.getEloss()/nsegments/_bulk.ENERGY_EHPAIR;
             
+//            System.out.println("length of subsegments: " + segment_length);
+//            System.out.println("total charge: " + segment_charge);
+            
             Hep3Vector segment_step = VecOp.mult(segment_length,track.getDirection());
             Hep3Vector segment_center = VecOp.add( track.getP1(),VecOp.mult(0.5,segment_step) );
             
+//            System.out.println("Segment step: " + segment_step);
+//            System.out.println("Segment center: " + segment_center);
+            
             // Loop over segments
             for (int iseg = 0; iseg < nsegments; iseg++)
             {
-                // THROW POISSON WITH SEGMENT CHARGE?
+                // FIXME: Add correct straggling treatment for thin layers
                 
                 // loop over sides of detector
                 for (ChargeCarrier carrier : ChargeCarrier.values())
                 {
                     if (hasElectrodesOnSide(carrier))
                      {
-//                      depositCharge(segment_charge, segment_center, carrier);       
- 
-                        // find center of charge deposition
-//                        double drift_destination = VecOp.dot( driftDestination(segment_center,carrier),_measured_coordinate.get(carrier) );
-//                        double diffusion_sigma = diffusionSigma(segment_center,carrier);
-//                        _electrodes.get(carrier).depositCharge(new BasicHep3Vector(drift_destination,0.0,0.0),segment_charge,diffusion_sigma);
-                        
                         Rotation3D sensor_to_electrodes = new Rotation3D();
-                        sensor_to_electrodes.setPassiveXYZ(0.0,0.0,-_electrode_angles.get(carrier));
+                        sensor_to_electrodes.setPassiveXYZ(0.0,0.0,-getElectrodeAngle(carrier));
                         
                         Hep3Vector electrode_drift_destination = VecOp.mult(sensor_to_electrodes.getRotationMatrix(),driftDestination(segment_center,carrier));
-                        ErrorEllipse2D electrode_charge_distribution = diffusionEllipse(segment_center,carrier).rotate(-_electrode_angles.get(carrier));
+                        ErrorEllipse2D electrode_charge_distribution = diffusionEllipse(segment_center,carrier).rotate(-getElectrodeAngle(carrier));
                        _electrodes.get(carrier).depositCharge(segment_charge,electrode_drift_destination,electrode_charge_distribution);
                     }
                 }
@@ -468,7 +481,128 @@
             }
               
         }
-
+        
+        _track_list.clear();
+        
+    }
+    
+    public String toString()
+    {
+        String newline = System.getProperty("line.separator");
+        String output = "SiSensor Object: "+newline+
+                "   Property 1";
+        return output;
     }
+   
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+//    private void depositCharge(double charge, Hep3Vector origin, ChargeCarrier carrier)
+//    {
+//        if (!hasStripsOnSide(carrier)) return;
+//        
+//        // find center of charge deposition
+//        double drift_destination = VecOp.dot( driftDestination(origin,carrier),
+//                _measured_coordinate.get(carrier) );
+//        double diffusion_sigma = diffusionSigma(origin,carrier);
+//        
+//        _electrodes.get(carrier).depositCharge(new BasicHep3Vector(drift_destination,0.0,0.0),charge,diffusion_sigma);
+//      
+//    }    
+    
+//    private double diffusionSigma(Hep3Vector point, ChargeCarrier carrier)
+//    {
+//        
+//        // Common factors
+//        double difference_V = _bias_voltage - _depletion_voltage;
+//        double sum_V = _bias_voltage + _depletion_voltage;
+//        double common_factor = 2.0*distanceFromSide(point,carrier)*_depletion_voltage/_thickness;
+//
+//        // Calculate charge spreading without magnetic field
+//        double sigmasq = k_Boltzmann * _bulk.getTemperature() * _thickness*_thickness / _depletion_voltage;
+//        if (_bulk.isNtype() == (carrier==ChargeCarrier.HOLE))
+//        {
+//            sigmasq *= Math.log( sum_V / (sum_V - common_factor));
+//        }
+//        else
+//        {
+//            sigmasq *= Math.log( (difference_V + common_factor) / difference_V );
+//        }
+//        
+//        double sigma = Math.sqrt(sigmasq);
+//        
+//        // Corrections for magnetic field -- this is an approximation, may have to be done better for high fields
+//        double cos_theta_lorentz_sq = Math.pow(VecOp.cosTheta(_drift_direction.get(carrier)),2);
+//        double phi_lorentz = VecOp.phi(_drift_direction.get(carrier));
+//        double phi_measured = VecOp.phi(_measured_coordinate.get(carrier));
+//        double cos_phi_diff_sq = Math.pow(Math.cos(phi_measured - phi_lorentz),2);
+//        
+//        sigma *= (1.0/cos_theta_lorentz_sq) *
+//                Math.sqrt(cos_theta_lorentz_sq + cos_phi_diff_sq - cos_theta_lorentz_sq*cos_phi_diff_sq);
+//        
+//        return sigma;
+//        
+//    }
+
+    
+// Delta-ray code was deprecated in favor of letting GEANT do the work    
+//
+//    public void generateDeltaRays()
+//    {
+//        for (TrackSegment track : _track_list)
+//        {
+//            // Uncommitted standalone code exists to do this, which...
+//            // retrieves track from _track_list and calculates delta ray production
+//            // modifies original track in _track_list
+//            // adds delta rays to _track_list
+//        }
+//        return;
+//    }
+
+//    private int nSegments(TrackSegment track, ChargeCarrier carrier, double deposition_granularity)
+//    {
+//        // Decide how to cut track into pieces as a fraction of strip pitch
+//        if (!hasElectrodesOnSide(carrier)) return 0;
+//        Hep3Vector deposition_line = VecOp.sub( driftDestination(track.getP2(),carrier),
+//                driftDestination(track.getP1(),carrier) );
+//
+//        double projected_deposition_length = VecOp.dot(deposition_line,_measured_coordinate.get(carrier));
+//        
+//        return (int)Math.ceil(projected_deposition_length/(deposition_granularity*_electrodes.get(carrier).getPitch()));
+//    }
+    
+    
     
 }

GeomConverter/src/org/lcsim/detector/tracker/silicon
SiSensorElectrodes.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- SiSensorElectrodes.java	30 Apr 2007 21:15:33 -0000	1.2
+++ SiSensorElectrodes.java	7 May 2007 21:22:49 -0000	1.3
@@ -43,7 +43,7 @@
     public void depositCharge(double charge, Hep3Vector position, ErrorEllipse2D distribution);
     
     // Get map of charge deposition on electrodes
-    public Map<Integer,Double> getChargeMap();
+    public Map<Integer,Integer> getChargeMap();
     
     // Clear charge from electrodes
     public void clear();

GeomConverter/src/org/lcsim/detector/tracker/silicon
SiStrips.java 1.3 -> 1.4
diff -u -r1.3 -r1.4
--- SiStrips.java	30 Apr 2007 21:15:33 -0000	1.3
+++ SiStrips.java	7 May 2007 21:22:49 -0000	1.4
@@ -29,19 +29,10 @@
     private int _nstrips; // number of strips
     private double _pitch; // sense pitch
     private int _floating_strips; // number of floating strips between readout strips
-    private SortedMap<Integer,Double> _strip_charge = new TreeMap<Integer,Double>();
+    private SortedMap<Integer,Integer> _strip_charge = new TreeMap<Integer,Integer>();
 
-//    private double _capacitance = 20.0; // 20 pF capacitance
-//    private double _inactive_width = 0.0; // no inactive region
-//    private double _strip_charge[] = null;
-
-    
-//    private double _crosstalk_fraction = 0.0;
-//    private double _integration_fraction = 1.0;
-//    private SortedMap<Integer,Integer> _strip_adc;
-    
-//    private SortedMap<Integer,Double> _strip_charge;
-//    private ArrayList<SiStripCluster> _clusters;
+    private double _capacitance_total = 16.0; // 15 pF
+    private double _capacitance_interstrip = 13.0; // 11 pF
     
     // Constructors
     //=============
@@ -53,13 +44,7 @@
         setFloatingStrips(floating_strips);
     }
     
-//    public SiStrips(int nstrips)
-//    {
-//        _nstrips = nstrips;
-//        _strip_charge = new double[_nstrips];
-//    }
-    
-    // Settors/Accessors
+    // Setters
     //==================
     public void setNStrips(int nstrips)
     {
@@ -81,31 +66,84 @@
 //        _capacitance = capacitance;
 //    }
 
-//    public void setInactiveWidth(double inactive_width)
-//    {
-//        _inactive_width = inactive_width;
-//    }
+ 
+    // Getters
+    //===================
+    private int getFloatingStrips()
+    {
+        return _floating_strips;
+    }
+    
+    private boolean isFloatingStrip(int sense_strip)
+    {
+        System.out.println("    Sense strip: "+sense_strip);
+        System.out.println("    Remainder: "+Math.IEEEremainder(sense_strip,_floating_strips+1));
+        System.out.println("Math.IEEEremainder(sense_strip,_floating_strips+1) != 0: " + (Math.IEEEremainder(sense_strip,_floating_strips+1) != 0));
+        return ( Math.IEEEremainder(sense_strip,_floating_strips+1) != 0 );
+    }
     
-    public double getPitch()
+    private int senseIDToReadoutID(int sense_strip_number)
+    {
+        return sense_strip_number/(_floating_strips+1);
+    }
+   
+    private double getSensePitch()
     {
         return _pitch;
     }
     
-    public int getFloatingStrips()
+    private double getReadoutPitch()
     {
-        return _floating_strips;
+        return _pitch*(1+getFloatingStrips());
+    }   
+    
+    private int getNSenseStrips()
+    {
+        return _nstrips;
+    }
+    
+    private int getNReadoutStrips()
+    {
+       return (_nstrips-1)/(_floating_strips+1) + 1;
+    }
+
+    private int getSenseStripID(Hep3Vector position)
+    {
+        return (int)Math.floor(position.x()/getSensePitch());
     }
     
+    private int getReadoutStripID(Hep3Vector position)
+    {
+        return (int)Math.floor(position.x()/getReadoutPitch());
+    }
+    
+    private Hep3Vector getPositionInSenseCell(Hep3Vector position)
+    {
+        return VecOp.add(position,new BasicHep3Vector(-getSenseStripID(position)*getSensePitch(), 0.0, 0.0));
+    }
+    
+    private Hep3Vector getPositionInReadoutCell(Hep3Vector position)
+    {
+        return VecOp.add(position,new BasicHep3Vector(-getReadoutStripID(position)*getReadoutPitch(), 0.0, 0.0));
+    }
+    
+    private Hep3Vector getSenseStripPosition(int sense_strip_number)
+    {
+        return new BasicHep3Vector(sense_strip_number*getSensePitch(),0.0,0.0);
+    }
+    
+    private Hep3Vector getReadoutStripPosition(int readout_strip_number)
+    {
+        return new BasicHep3Vector(readout_strip_number*getReadoutPitch(),0.0,0.0);
+    }
+    
+    
+    
 //    public double getCapacitance()
 //    {
 //        return _capacitance;
 //    }
     
-//    public double getInactiveWidth()
-//    {
-//        return _inactive_width;
-//    }
-    
     // Operators
     //==========
     public int getNAxes()
@@ -115,138 +153,161 @@
     
     public int getNCells()
     {
-        return _nstrips;
+        return getNReadoutStrips();
     }
     
     public int getNCells(int direction)
     {
         if (direction == 0)
         {
-            return _nstrips;
+            return getNReadoutStrips();
         }
         else return 1;
     }
     
     public double getPitch(int direction)
     {
-        return _pitch;
+        return getReadoutPitch();
     }
-    
+     
+    // FIXME: cell ID assignment is not correct in all cases
     public int getCellID(Hep3Vector position)
     {
-        return (int)Math.floor(position.x()/_pitch);
+        return getReadoutStripID(position);
     }
     
     public Hep3Vector getPositionInCell(Hep3Vector position)
     {
-        return VecOp.add(position,new BasicHep3Vector(-getCellID(position)*_pitch, 0.0, 0.0));
+        return getPositionInReadoutCell(position);
     }
    
     public Hep3Vector getCellPosition(int strip_number)
     {
-        return new BasicHep3Vector(strip_number*_pitch,0.0,0.0);
+        return getReadoutStripPosition(strip_number);
     }
     
-    public SortedMap<Integer,Double> getChargeMap()
+    public SortedMap<Integer,Integer> getChargeMap()
     {
         return _strip_charge;
     }
-    
-//    public double stripPosition(double strip_coordinate)
-//    {
-//        return strip_coordinate*_pitch;
-//    }    
-    
-//    public void depositCharge(Hep3Vector position, double charge, double sigma)
-//    {
-//        int base_strip = getCellID(position);
-//        Hep3Vector interstrip_position = getPositionInCell(position);
-//        
-//        // put charge on strips in window 3-sigma strips on each side of base strip
-//        double pitch = getPitch(1);
-//
-//        int window_size = (int)Math.ceil(3.0*sigma/pitch);
-//        
-//        double erf_min = 1.0;
-//        double erf_max = 1.0;
-//        for (int istrip = base_strip-window_size; istrip <= base_strip+window_size; istrip++)
-//        {
-//            
-//            try
-//            {
-//                erf_max = Erf.erf( (_electrodes.get(carrier).stripPosition(istrip+0.5) - drift_destination)/diffusion_sigma );
-//            }
-//            catch (MathException no_convergence)
-//            {
-//                System.out.println("erf fails to converge!!");
-//            }
-//            
-//            int strip_charge = (int)Math.round( (erf_min-erf_max) * charge );
-//            addCharge(istrip,strip_charge);
-//            
-//            erf_min = erf_max;
-//        }                  
-//        return;
-//    }
-    
-    
+        
     public void depositCharge(double charge, Hep3Vector position, ErrorEllipse2D distribution)
     {
-        int base_strip = getCellID(position);
-        Hep3Vector interstrip_position = getPositionInCell(position);
-       
+//        System.out.println("Beginning depositCharge");
+//        System.out.println("    Charge: "+charge);        
+//        System.out.println("    Position: "+position);       
+        
+        int base_strip = getSenseStripID(position);
+        
+//        System.out.println("    Base Strip: "+base_strip);        
+        
+        Hep3Vector interstrip_position = getPositionInSenseCell(position);
+
+//        System.out.println("    Interstrip position: "+interstrip_position);        
+        
         // put charge on strips in window 3-sigma strips on each side of base strip
-        double pitch = getPitch(1);
+        double pitch = getSensePitch();
         double axis_angle = 0.0;
-
+        
+//        System.out.println("    Pitch: "+pitch);                
+        
         int window_size = (int)Math.ceil(3.0*distribution.sigma1D(axis_angle)/pitch);
         
-        double erf_lower = 1.0;
-        double erf_upper = 1.0;
+//        System.out.println("    Sigma: "+distribution.sigma1D(axis_angle));
+//        System.out.println("    Window Size: "+window_size);        
+                
+        double erfc_lower = 1.0;
+        double erfc_upper = 1.0;
     
         for (int istrip = base_strip-window_size; istrip <= base_strip+window_size; istrip++)
         {
-            double cell_edge_upper = getCellPosition(istrip).x() + pitch/2.0;
-            erf_upper = distribution.erf1D(cell_edge_upper,axis_angle);
+//            System.out.println("    istrip: "+istrip);                    
+            
+            double cell_edge_upper = getSenseStripPosition(istrip).x() + pitch/2.0;
+//            System.out.println("    Cell upper edge: "+cell_edge_upper);                                
+            
+            double erfc_limit = cell_edge_upper-position.x();
+            
+//            System.out.println("    erfc_limit: "+erfc_limit);                                            
+            
+            erfc_upper = distribution.erfc1D(erfc_limit,axis_angle);
 
-            int strip_charge = (int)Math.round( (erf_lower-erf_upper) * charge); 
-            addCharge(istrip,strip_charge);
+//            System.out.println("    erfc_lower: "+erfc_lower);
+//            System.out.println("    erfc_upper: "+erfc_upper);
+            
+            if (erfc_lower<erfc_upper) System.out.println("SQUEAL LIKE A PIG!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
             
-            erf_lower = erf_upper;
+            int strip_charge = (int)Math.round( (erfc_lower-erfc_upper) * charge);
+                        
+            if (strip_charge != 0)
+            {
+                addChargeToSense(istrip,strip_charge);
+            }
+            
+            erfc_lower = erfc_upper;
             
         }
         return;
     }
-    
-    
-    public void addCharge(int strip, double charge)
-    { 
-        if (_strip_charge.containsKey(strip)) {
-            charge += _strip_charge.get(strip);            
-        }   
-        _strip_charge.put(strip,charge);
-    }    
-    
+        
     public void clear()
     {
        _strip_charge.clear();
     }
-
-
     
-//    public void generateReadout()
-//    {
-//        
-//    }
-//    
-//    public void findClusters()
-//    {
-//        
-//    }
-//    
-//    public void findTruthCluster()
-//    {
-//        
-//    }
     
+    // Actions
+    //==================
+    private void addChargeToSense(int sense_strip, int charge)
+    {
+        
+        System.out.println(" Adding charge to sense strip: "+sense_strip);
+        
+        if (!isFloatingStrip(sense_strip))
+        {    
+            System.out.println(" Found central readout strip: "+sense_strip);
+            addChargeToReadout(senseIDToReadoutID(sense_strip),charge);
+        }
+        else
+        {
+            
+            System.out.println(" NOT central readout strip: "+sense_strip);
+            
+            int ileft = sense_strip-1;
+            int charge_left = charge/2;
+            while (isFloatingStrip(ileft))
+            {
+                System.out.println(" NOT left readout strip: "+ileft);
+                charge_left *= (_capacitance_interstrip/_capacitance_total);
+                ileft--;
+            }
+            System.out.println(" Found left readout strip: "+ileft);
+            addChargeToReadout(senseIDToReadoutID(ileft),charge_left);
+            
+            int iright = sense_strip+1;
+            int charge_right = charge/2;
+            while (isFloatingStrip(iright))
+            {
+                System.out.println(" NOT right readout strip: "+iright);
+                charge_right *= (_capacitance_interstrip/_capacitance_total);
+                iright++;
+            }
+            System.out.println(" Found right readout strip: "+iright);
+            addChargeToReadout(senseIDToReadoutID(iright),charge_right);
+        }
+    }
+    
+        
+    private void addChargeToReadout(int readout_strip, int charge)
+    {
+        if (charge == 0) return;
+        
+        if (_strip_charge.containsKey(readout_strip)) {
+            charge += _strip_charge.get(readout_strip);            
+        }   
+        _strip_charge.put(readout_strip,charge);
+        return;
+    }    
+    
+
 }
CVSspam 0.2.8