5 added files
lcsim/src/org/lcsim/contrib/SiStripSim
diff -N SiStripDetector.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SiStripDetector.java 3 Aug 2006 23:25:04 -0000 1.1
@@ -0,0 +1,317 @@
+package org.lcsim.contrib.SiStripSim;
+/*
+ * SiStripDetector.java
+ *
+ * Created on July 20, 2005, 5:55 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+import org.apache.commons.math.special.Erf;
+import org.apache.commons.math.MathException;
+
+import java.util.*;
+
+/**
+ *
+ * @author tknelson
+ */
+public class SiStripDetector {
+
+ // Enumerated types
+ //=================
+ protected enum Orientation { PINSIDE, POUTSIDE };
+
+ // Fields
+ //=======
+
+ // primary properties
+ //-------------------
+ // strips
+ private Map<ChargeCarrier,SiStrips> _strips = new EnumMap<ChargeCarrier,SiStrips>(ChargeCarrier.class);
+ private Map<ChargeCarrier,Double> _strip_angles = new EnumMap<ChargeCarrier,Double>(ChargeCarrier.class);
+ private Orientation _orientation = Orientation.POUTSIDE;
+
+ // bulk
+ private DopedSilicon _bulk = null;
+ private double _thickness = 0.300;
+
+ // global
+ private Hep3Vector _b_field = new BasicHep3Vector(0.0,0.0,0.0);
+ private double _depletion_voltage = 100.0;
+ private double _bias_voltage = 1.1 * _depletion_voltage;
+
+ // derived properties
+ //-------------------
+ private EnumMap<ChargeCarrier,Hep3Vector> _measured_coordinate = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
+ private EnumMap<ChargeCarrier,Hep3Vector> _strip_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
+ private EnumMap<ChargeCarrier,Hep3Vector> _drift_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
+ private List<TrackSegment> _track_list = new ArrayList<TrackSegment>();
+
+ // Constructors
+ //=============
+
+ // Default, 300 micron, p-in-n silicon; 0-degree, 50 micron single-sided readout with one floating strip
+ public SiStripDetector()
+ {
+ _strips.put(ChargeCarrier.HOLE, new SiStrips());
+ _strip_angles.put(ChargeCarrier.HOLE, 0.0);
+ _bulk = new DopedSilicon();
+ initialize();
+ }
+
+ // Pass p strips
+ public SiStripDetector(SiStrips pstrips, double pstrip_angle)
+ {
+ _strips.put(ChargeCarrier.HOLE,pstrips);
+ _strip_angles.put(ChargeCarrier.HOLE,pstrip_angle);
+ _orientation = Orientation.POUTSIDE;
+ _bulk = new DopedSilicon();
+ initialize();
+ }
+
+ // Pass bulk
+ public SiStripDetector(DopedSilicon bulk)
+ {
+ this();
+ _bulk = bulk;
+ initialize();
+ }
+
+ // Pass p strips and bulk
+ public SiStripDetector(SiStrips pstrips, double pstrip_angle, DopedSilicon bulk)
+ {
+ this(pstrips, pstrip_angle);
+ _bulk = bulk;
+ initialize();
+ }
+
+ // Accessors
+ //==========
+ public void setBField(Hep3Vector b_field)
+ {
+ _b_field = b_field;
+ }
+
+ // Operators
+ //==========
+ private void initialize()
+ {
+ // Store various important directions
+ for (ChargeCarrier carrier : ChargeCarrier.values())
+ {
+ if (hasStripsOnSide(carrier)) {
+ _drift_direction.put(carrier,driftDirection(carrier));
+ double strip_angle = _strip_angles.get(carrier);
+ _measured_coordinate.put(carrier,measuredCoordinate(strip_angle));
+ _strip_direction.put(carrier,stripDirection(strip_angle));
+ }
+ }
+ }
+
+ private double zOfSide(ChargeCarrier carrier)
+ {
+ if ( (carrier == ChargeCarrier.HOLE) == (_orientation == Orientation.POUTSIDE) ) return _thickness;
+ else return 0;
+ }
+
+ private double distanceFromSide(Hep3Vector point, ChargeCarrier carrier)
+ {
+ return point.z() - zOfSide(carrier);
+ }
+
+ private boolean hasStripsOnSide(ChargeCarrier carrier)
+ {
+ if (_strips.get(carrier) == null) return false;
+ else return true;
+ }
+
+ private Hep3Vector driftVector(Hep3Vector origin, ChargeCarrier carrier)
+ {
+ 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)
+ {
+ return VecOp.add(origin,driftVector(origin, carrier));
+ }
+
+ 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
+ double sigmasq = _bulk.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 strip_angle)
+ {
+ return new BasicHep3Vector(Math.cos(strip_angle),Math.sin(strip_angle),0.0);
+ }
+
+ private Hep3Vector stripDirection(double strip_angle)
+ {
+ return measuredCoordinate(strip_angle + (Math.PI/4.0)) ;
+ }
+
+ private Hep3Vector driftDirection(ChargeCarrier carrier)
+ {
+ double carrier_mobility = _bulk.mobility(carrier);
+
+ // 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 void clearStrips()
+ {
+ for (ChargeCarrier carrier : ChargeCarrier.values())
+ {
+ if (hasStripsOnSide(carrier)) _strips.get(carrier).clear();
+ }
+ }
+
+// Delta-ray code now 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 (!hasStripsOnSide(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*_strips.get(carrier).getPitch()));
+ }
+
+ 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) );
+ int base_strip = _strips.get(carrier).baseStrip(drift_destination);
+ double interstrip_position = _strips.get(carrier).interstripPosition(drift_destination);
+
+ // put charge on strips in window 3-sigma strips on each side of base strip
+ double pitch = _strips.get(carrier).getPitch();
+ double diffusion_sigma = diffusionSigma(origin,carrier);
+ int window_size = (int)Math.ceil(3.0*diffusion_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( (_strips.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 );
+ _strips.get(carrier).addCharge(istrip,strip_charge);
+
+ erf_min = erf_max;
+ }
+ }
+
+ public void depositCharge()
+ {
+
+ for (TrackSegment track : _track_list)
+ {
+
+ // Decide how to cut track into pieces - use 5% of pitch
+ int nsegments = 0;
+ for (ChargeCarrier carrier : ChargeCarrier.values())
+ {
+ if (!hasStripsOnSide(carrier)) continue;
+ nsegments = Math.max(nsegments,nSegments(track,carrier, 0.05));
+ }
+
+ // Set up segments
+ double segment_length = track.getLength()/nsegments;
+ double segment_charge = track.getEloss()/nsegments/_bulk.E_PAIR;
+
+ Hep3Vector segment_step = VecOp.mult(segment_length,track.getDirection());
+ Hep3Vector segment_center = VecOp.add( track.getP1(),VecOp.mult(0.5,segment_step) );
+
+ // Loop over segments
+ for (int iseg = 0; iseg < nsegments; iseg++)
+ {
+ // THROW POISSON WITH SEGMENT CHARGE?
+
+ // loop over sides of detector
+ for (ChargeCarrier carrier : ChargeCarrier.values())
+ {
+ if (!hasStripsOnSide(carrier)) continue;
+ depositCharge(segment_charge, segment_center, carrier);
+ }
+
+ // step to next segment
+ segment_center = VecOp.add(segment_center, segment_step);
+ }
+
+ }
+
+ }
+
+}
lcsim/src/org/lcsim/contrib/SiStripSim
diff -N ChargeCarrier.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ChargeCarrier.java 3 Aug 2006 23:25:04 -0000 1.1
@@ -0,0 +1,54 @@
+package org.lcsim.contrib.SiStripSim;
+/*
+ * ChargeCarrier.java
+ *
+ * Created on October 13, 2005, 3:41 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+/**
+ *
+ * @author tknelson
+ */
+
+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);
+
+ private final double _mu_0_factor;
+ private final double _mu_0_exponent;
+ private final double _mu_min_factor;
+ private final double _mu_min_exponent;
+ private final double _N_ref_factor;
+ private final double _N_ref_exponent;
+ 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,
+ double N_ref_factor, double N_ref_exponent, double alpha_factor, double alpha_exponent)
+ {
+ _mu_0_factor = mu_0_factor;
+ _mu_0_exponent = mu_0_exponent;
+ _mu_min_factor = mu_min_factor;
+ _mu_min_exponent = mu_min_exponent;
+ _N_ref_factor = N_ref_factor;
+ _N_ref_exponent = N_ref_exponent;
+ _alpha_factor = alpha_factor;
+ _alpha_exponent = alpha_exponent;
+ }
+
+ // Methods
+ double mu0(double temperature)
+ {return _mu_0_factor * Math.pow( (temperature/300.0), _mu_0_exponent);}
+ double muMin(double temperature)
+ {return _mu_min_factor * Math.pow( (temperature/300.0), _mu_min_exponent);}
+ double nRef(double temperature)
+ {return _N_ref_factor * Math.pow( (temperature/300.0), _N_ref_exponent);}
+ double alpha(double temperature)
+ {return _alpha_factor * Math.pow( (temperature/300.0), _alpha_exponent);}
+
+}
lcsim/src/org/lcsim/contrib/SiStripSim
diff -N TrackSegment.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackSegment.java 3 Aug 2006 23:25:04 -0000 1.1
@@ -0,0 +1,118 @@
+package org.lcsim.contrib.SiStripSim;
+/*
+ * TrackSegment.java
+ *
+ * Created on July 27, 2005, 3:34 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.geometry.IDDecoder;
+
+import java.util.*;
+
+/**
+ *
+ * @author tknelson
+ */
+public class TrackSegment
+{
+
+ // Fields
+ private Hep3Vector _p1;
+ private Hep3Vector _p2;
+ private double _energy_loss;
+
+ /**
+ * Creates a new instance of TrackSegment
+ */
+
+ public TrackSegment(Hep3Vector p1, Hep3Vector p2, double energy_loss)
+ {
+ _p1 = p1;
+ _p2 = p2;
+ _energy_loss = energy_loss;
+ }
+
+// Construct from a SimTrackerHit: currently broken
+//
+// public TrackSegment(SimTrackerHit simulated_hit, IDDecoder decoder)
+// {
+// Hep3Vector midpoint = new BasicHep3Vector(simulated_hit.getPoint());
+// Hep3Vector direction = VecOp.unit(new BasicHep3Vector(simulated_hit.getMomentum()));
+//
+// BasicHep3Vector normal = new BasicHep3Vector();
+// decoder.setID(simulated_hit.getCellID());
+// if (decoder.getBarrelEndcapFlag().isBarrel())
+// {
+// normal.setV(midpoint.x(), midpoint.y(), 0.0);
+// normal = (BasicHep3Vector)VecOp.unit(normal);
+// }
+// else if (decoder.getBarrelEndcapFlag().isEndcap())
+// {
+// normal.setV(0.0,0.0,1.0);
+// }
+//
+// Hep3Vector half_length = VecOp.mult(simulated_hit.getPathLength(),direction);
+//
+// }
+
+ // Accessors
+ public Hep3Vector getP1()
+ {
+ return _p1;
+ }
+
+ public Hep3Vector getP2()
+ {
+ return _p2;
+ }
+
+ public double getEloss()
+ {
+ return _energy_loss;
+ }
+
+ public Hep3Vector getVector()
+ {
+ return VecOp.sub(_p2,_p1);
+ }
+
+ public Hep3Vector getDirection()
+ {
+ return VecOp.unit(getVector());
+ }
+
+ public double getLength()
+ {
+ return getVector().magnitude();
+ }
+
+ public double getDedx()
+ {
+ return _energy_loss/getLength();
+ }
+
+ public void rotateMedium(Hep3Vector axis, double center, double angle)
+ {
+ double old_length = getLength();
+ rotateMedium(_p1,axis,center,angle);
+ rotateMedium(_p2,axis,center,angle);
+ _energy_loss *= getLength()/old_length;
+ }
+
+ private void rotateMedium(Hep3Vector point, Hep3Vector axis, double center, double angle)
+ {
+ assert (axis.z() == 0.0); // assumed that we are rotating silicon z-planes
+ point = VecOp.add(point,VecOp.mult(Math.tan(angle)*(center - point.z()),
+ VecOp.unit(VecOp.cross(axis, new BasicHep3Vector(0,0,1)))));
+ }
+
+}
lcsim/src/org/lcsim/contrib/SiStripSim
diff -N SiStrips.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SiStrips.java 3 Aug 2006 23:25:04 -0000 1.1
@@ -0,0 +1,143 @@
+package org.lcsim.contrib.SiStripSim;
+/*
+ * SiStrips.java
+ *
+ * Created on July 22, 2005, 4:07 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+import java.util.SortedMap;
+import java.util.Arrays;
+
+/**
+ *
+ * @author tknelson
+ */
+public class SiStrips
+{
+
+ private double _pitch = 0.025; // 25 micron sense pitch
+ private int _floating_strips = 1; // 50 micron readout
+ private double _capacitance = 20.0; // 20 pF capacitance
+ private int _nstrips = 100; // simulate 100 strips by default
+// 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;
+
+ // Constructors
+ //=============
+ public SiStrips()
+ {
+ _strip_charge = new double[_nstrips];
+ }
+
+ public SiStrips(int nstrips)
+ {
+ _nstrips = nstrips;
+ _strip_charge = new double[_nstrips];
+ }
+
+ // Settors/Accessors
+ //==================
+ public void setPitch(double pitch)
+ {
+ _pitch = pitch;
+ }
+
+ public void setFloatingStrips(int floating_strips)
+ {
+ _floating_strips = floating_strips;
+ }
+
+ public void setCapacitance(double capacitance)
+ {
+ _capacitance = capacitance;
+ }
+
+// public void setInactiveWidth(double inactive_width)
+// {
+// _inactive_width = inactive_width;
+// }
+
+ public double getNstrips()
+ {
+ return _nstrips;
+ }
+
+ public double getPitch()
+ {
+ return _pitch;
+ }
+
+ public int getFloatingStrips()
+ {
+ return _floating_strips;
+ }
+
+ public double getCapacitance()
+ {
+ return _capacitance;
+ }
+
+// public double getInactiveWidth()
+// {
+// return _inactive_width;
+// }
+
+ // Operators
+ //==========
+ public int baseStrip(double position)
+ {
+ return (int)Math.floor(position/_pitch);
+ }
+
+ public double interstripPosition(double position)
+ {
+ return position - baseStrip(position)*_pitch;
+ }
+
+ public double stripPosition(int strip_number)
+ {
+ return strip_number*_pitch;
+ }
+
+ public double stripPosition(double strip_coordinate)
+ {
+ return strip_coordinate*_pitch;
+ }
+
+ public void clear()
+ {
+ Arrays.fill(_strip_charge,0.0);
+ }
+
+ public void addCharge(int strip, double charge)
+ {
+ _strip_charge[strip] += charge;
+ }
+
+// public void generateReadout()
+// {
+//
+// }
+//
+// public void findClusters()
+// {
+//
+// }
+//
+// public void findTruthCluster()
+// {
+//
+// }
+
+}
lcsim/src/org/lcsim/contrib/SiStripSim
diff -N DopedSilicon.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DopedSilicon.java 3 Aug 2006 23:25:04 -0000 1.1
@@ -0,0 +1,84 @@
+package org.lcsim.contrib.SiStripSim;
+/*
+ * DopedSilicon.java
+ *
+ * Created on July 26, 2005, 3:31 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+import java.util.*;
+
+/**
+ *
+ * @author tknelson
+ */
+public class DopedSilicon
+{
+
+ // Fields
+ //=======
+
+ // Static
+ static public double K_BOLTZMANN = 8.617385E-5; // eV/degK
+ static public double E_PAIR = 3.6E-9; // 3.6E-9 GeV/e-h pair
+
+ // Member
+ private double _temperature = 293.0;
+ private double _doping_concentration = 6.0E+11;
+ private EnumMap<ChargeCarrier,Double> _carrier_concentration = new EnumMap<ChargeCarrier,Double>(ChargeCarrier.class);
+
+ // Constructors
+ //=============
+ public DopedSilicon()
+ {
+ setElectronConcentration(1.0E+14);
+ setHoleConcentration(1.0E+14);
+ }
+
+ // Setters/Accessors
+ //==================
+ public void setTemperature(double temperature){_temperature = temperature;}
+ public void setDopingConcentration(double doping_concentration)
+ {
+ _doping_concentration = doping_concentration;
+ }
+ public void setElectronConcentration(double electron_concentration)
+ {
+ _carrier_concentration.put(ChargeCarrier.ELECTRON,electron_concentration);
+ }
+ public void setHoleConcentration(double hole_concentration)
+ {
+ _carrier_concentration.put(ChargeCarrier.HOLE,hole_concentration);
+ }
+
+ public double getTemperature(){return _temperature;}
+ public double getCarrierConcentration(ChargeCarrier charge_carrier)
+ {
+ return _carrier_concentration.get(charge_carrier);
+ }
+
+ // Methods
+ //========
+
+ // Lorentz angle calculation for silicon sensors
+ public double tanLorentzAngle(double b_field, ChargeCarrier charge_carrier)
+ {
+ return b_field * mobility(charge_carrier) * 1.0E-4;
+ }
+
+ // Mobility calculation with correction for irradiated sensors
+ public double mobility(ChargeCarrier charge_carrier)
+ {
+ return charge_carrier.muMin(_temperature) + charge_carrier.mu0(_temperature) /
+ (1.0 + Math.pow(_carrier_concentration.get(charge_carrier)/charge_carrier.nRef(_temperature),
+ charge_carrier.alpha(_temperature)));
+ }
+
+ // Silicon type
+ public boolean isNtype(){return _doping_concentration > 0.0;}
+
+
+}
CVSspam 0.2.8
