lcsim/src/org/lcsim/contrib/SiStripSim
diff -N SiStripDetector.java
--- SiStripDetector.java 3 Aug 2006 23:25:04 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,317 +0,0 @@
-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);
- }
-
- }
-
- }
-
-}
