View Javadoc

1   /*
2    * CDFSiSensorSim.java
3    *
4    * Created on May 10, 2007, 3:03 PM
5    *
6    * To change this template, choose Tools | Template Manager
7    * and open the template in the editor.
8    */
9   
10  package org.lcsim.recon.tracking.digitization.sisim;
11  
12  import org.lcsim.detector.solids.GeomOp3D;
13  import org.lcsim.event.SimTrackerHit;
14  import org.lcsim.detector.tracker.silicon.SiSensor;
15  import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
16  import org.lcsim.detector.tracker.silicon.ChargeCarrier;
17  import org.lcsim.detector.ITransform3D;
18  import org.lcsim.detector.IRotation3D;
19  import org.lcsim.detector.tracker.silicon.ChargeDistribution;
20  import org.lcsim.detector.tracker.silicon.GaussianDistribution2D;
21  
22  import hep.physics.vec.Hep3Vector;
23  import hep.physics.vec.BasicHep3Vector;
24  import hep.physics.matrix.BasicMatrix;
25  import hep.physics.vec.VecOp;
26  
27  import java.util.ArrayList;
28  import java.util.List;
29  import java.util.Map;
30  import java.util.EnumMap;
31  import java.util.Map.Entry;
32  import java.util.SortedMap;
33  
34  import org.lcsim.detector.Transform3D;
35  import org.lcsim.detector.Translation3D;
36  import org.lcsim.detector.solids.Line3D;
37  import org.lcsim.detector.solids.Point3D;
38  
39  /**
40   *
41   * @author tknelson
42   */
43  public class CDFSiSensorSim implements SiSensorSim
44  {
45      
46      // Fields
47      SiSensor _sensor = null;
48  //    ChargeTransferModel _transfer_model = null;
49      Map<ChargeCarrier,Hep3Vector> _drift_direction = null;
50      Map<ChargeCarrier,SiElectrodeDataCollection> _sense_data = null;
51      Map<ChargeCarrier,SiElectrodeDataCollection> _readout_data = null;
52  
53      // Simple simulation of charge trapping, this is a temporary kludge.
54      // Charge collection efficiency with linear drift distance dependence.
55      // Input is fraction lost per 100um drift: 0.2 is typical for 1E15 NEQ.
56      // FIXME: should be calculated from properties of DopedSilicon (radiation dose)
57  
58      double _trapping = 0.0;
59      
60  //    SiElectrodeSim
61      
62      // Static parameters - not intended to be user modifiable
63      private static double _DEPOSITION_GRANULARITY = 0.10; // 10% of pitch or depleted thickness
64      private static final double DISTANCE_ERROR_THRESHOLD = 0.001; //This is the maximum distance outside of the sensor allowed before an error is thrown. 
65      
66      private final boolean debug = false;
67      
68      /**
69       * Creates a new instance of CDFSiSensorSim
70       */
71      public CDFSiSensorSim()
72      {
73          _drift_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
74          _sense_data = new EnumMap<ChargeCarrier,SiElectrodeDataCollection>(ChargeCarrier.class);
75          _readout_data = new EnumMap<ChargeCarrier,SiElectrodeDataCollection>(ChargeCarrier.class);
76          for (ChargeCarrier carrier : ChargeCarrier.values())
77          {
78              _sense_data.put(carrier,new SiElectrodeDataCollection());
79              _readout_data.put(carrier,new SiElectrodeDataCollection());
80          }
81      }
82  
83      // Setters
84      public void setTrapping(double trapping)
85      {
86          _trapping = trapping;
87      }
88  
89      // Implementation of SiSensorSim interface
90      //==================================
91      
92      // Get charge map on electrodes
93      public SiElectrodeDataCollection getReadoutData(ChargeCarrier carrier)
94      {
95          return _readout_data.get(carrier);
96      }
97      
98      // SiSensorSim interface
99      //======================
100     
101     // Simulate charge deposition
102     public void setSensor(SiSensor sensor)
103     {
104         _sensor = sensor;
105     }
106     
107     public Map<ChargeCarrier,SiElectrodeDataCollection> computeElectrodeData()
108     {
109     	if(debug) {
110     		System.out.printf("%s: computeElectrodeData for sensor %s\n", this.getClass().getSimpleName(),this._sensor.getName());
111     		System.out.printf("%s: # Sense strips: %d\n", this.getClass().getSimpleName(),_sensor.getSenseElectrodes(ChargeCarrier.HOLE).getNCells(0));
112     		System.out.printf("%s: # Readout strips: %d\n", this.getClass().getSimpleName(),_sensor.getReadoutElectrodes(ChargeCarrier.HOLE).getNCells(0));
113     	}
114         depositChargeOnSense();
115         transferChargeToReadout();
116         
117         return _readout_data;
118     }
119     
120     // Clear readout data
121     public void clearReadout()
122     {
123         for (ChargeCarrier carrier : ChargeCarrier.values())
124         {
125             _readout_data.get(carrier).clear();
126         }
127     }
128     
129     public void lorentzCorrect(Hep3Vector position, ChargeCarrier carrier)
130     {
131         Line3D drift_line = new Line3D(new Point3D(position), driftDirection(carrier,new BasicHep3Vector(0,0,0))); // use drift direction at origin for now
132         
133         List<Point3D> intersection_points = new ArrayList<Point3D>();
134         
135         for (ChargeCarrier bias_carrier : ChargeCarrier.values())
136         {
137             intersection_points.add(GeomOp3D.intersection(drift_line,_sensor.getBiasSurface(bias_carrier).getPlane()));
138         }
139         
140         Hep3Vector corrected_position = VecOp.mult(0.5,VecOp.add(intersection_points.get(0),intersection_points.get(1)));
141         
142         ITransform3D transform = new Transform3D(new Translation3D(VecOp.sub(corrected_position,position)));
143         
144         transform.transform(position);
145     }
146     
147     // Private
148     private void clearSense()
149     {
150         for (ChargeCarrier carrier : ChargeCarrier.values())
151         {
152             _sense_data.get(carrier).clear();
153         }
154     }
155     
156     private void depositChargeOnSense()
157     {
158         
159         if(debug)
160     		System.out.printf("%s: depositChargeOnSense for sensor %s\n", this.getClass().getSimpleName(),this._sensor.getName());
161 
162         // Set up drift directions // FIXME: put this in a setup method.
163         for (ChargeCarrier carrier : ChargeCarrier.values())
164         {
165             if (_sensor.hasElectrodesOnSide(carrier))
166             {
167                 _drift_direction.put( carrier, driftDirection(carrier,new BasicHep3Vector(0.0,0.0,0.0)) );
168             }
169         }
170         
171         ITransform3D global_to_sensor = _sensor.getGeometry().getGlobalToLocal();
172         List<SimTrackerHit> hits = _sensor.getReadout().getHits(SimTrackerHit.class);
173         
174         for (SimTrackerHit hit : hits)
175         {
176             if(debug) 
177             		System.out.printf("%s: depositChargeOnSense for sim tracker hit at %s \n", this.getClass().getSimpleName(),hit.getPositionVec().toString());
178 
179 //            System.out.println("Hit point: " + "[" + hit.getPoint()[0] + "," + hit.getPoint()[1] + "," + hit.getPoint()[2] + "]");
180 //            System.out.println("Startpoint: " + "[" + hit.getStartPoint()[0] + "," + hit.getStartPoint()[1] + "," + hit.getStartPoint()[2] + "]");
181 //            System.out.println("Endpoint: " + "[" + hit.getEndPoint()[0] + "," + hit.getEndPoint()[1] + "," + hit.getEndPoint()[2] + "]");
182 //
183 //            System.out.println("\n"+"Startpoint is inside: " + _sensor.getGeometry().inside(new BasicHep3Vector(hit.getStartPoint())));
184 //            System.out.println("Endpoint is inside: " + _sensor.getGeometry().inside(new BasicHep3Vector(hit.getEndPoint())));
185             
186 //            if (_sensor.getGeometry().inside(new BasicHep3Vector(hit.getStartPoint())) == Inside.OUTSIDE ||
187 //                _sensor.getGeometry().inside(new BasicHep3Vector(hit.getEndPoint())) == Inside.OUTSIDE )
188 //            {
189 //                throw new RuntimeException("Endpoints of SimTrackerHit are outside the sensor volume!!");
190 //            }
191             
192             
193             TrackSegment track = new TrackSegment(hit);
194             track.transform(global_to_sensor);
195             
196             int nsegments = 0;
197             
198             // Compute number of segments
199             for (ChargeCarrier carrier : ChargeCarrier.values())
200             {
201                 if (_sensor.hasElectrodesOnSide(carrier))
202                 {
203 //                    _drift_direction.put( carrier, this.driftDirection(carrier,new BasicHep3Vector(0.0,0.0,0.0)) );
204                     nsegments = Math.max(nsegments,nSegments(track,carrier, _DEPOSITION_GRANULARITY));
205                 }
206             }
207             
208 //            System.out.println("Number of subsegments: " + nsegments);
209             
210             // Set up segments
211 //            double segment_length = hit.getPathLength()/nsegments;
212 //            double segment_charge = hit.getdEdx()/nsegments/_sensor.getBulk().ENERGY_EHPAIR;
213             
214 //            Hep3Vector segment_step = VecOp.mult(segment_length,VecOp.unit( new BasicHep3Vector(hit.getMomentum()) ));
215 //            Hep3Vector segment_center = VecOp.add( new BasicHep3Vector(hit.getStartPoint()), VecOp.mult(0.5,segment_step) );
216             
217             
218             
219             // Set up segments
220             double segment_length = track.getLength()/nsegments;
221             double segment_charge = track.getEloss()/nsegments/_sensor.getBulk().ENERGY_EHPAIR;
222 
223 //            System.out.println("length of subsegments: " + segment_length);
224 //            System.out.println("subsegment charge: " + segment_charge);
225             
226             Hep3Vector segment_step = VecOp.mult(segment_length,track.getDirection());
227             Hep3Vector segment_center = VecOp.add( track.getP1(),VecOp.mult(0.5,segment_step) );
228             
229 //            System.out.println("Segment step: " + segment_step);
230             
231             // Loop over segments
232             for (int iseg = 0; iseg < nsegments; iseg++)
233             {
234 //                System.out.println("Segment center: " + segment_center);
235                 
236                 // FIXME: Add correct straggling treatment for thin layers
237                 
238                 // loop over sides of detector
239                 for (ChargeCarrier carrier : ChargeCarrier.values())
240                 {
241                     if (_sensor.hasElectrodesOnSide(carrier))
242                     {
243                         SiSensorElectrodes electrodes = _sensor.getSenseElectrodes(carrier);
244 
245                         // Apply collection inefficiency for charge trapping: require between 0 and 1
246                         double collection_efficiency = 1.0 - 10*_trapping*
247                                 driftVector(segment_center,carrier).magnitude();
248                         collection_efficiency = Math.max(0.0,Math.min(1.0,collection_efficiency));
249                         segment_charge *= (collection_efficiency);
250                         
251                         ChargeDistribution charge_distribution = diffusionDistribution(segment_charge,segment_center,carrier);
252                         charge_distribution.transform(electrodes.getParentToLocal());
253                         
254                         SortedMap<Integer,Integer> sense_charge = electrodes.computeElectrodeData(charge_distribution);
255 //                        System.out.println("Sense charge map: " + sense_charge);
256                         
257                         _sense_data.get(carrier).add(new SiElectrodeDataCollection(sense_charge,hit));
258                     }
259                 }
260                 
261                 // step to next segment
262                 segment_center = VecOp.add(segment_center, segment_step);
263             }
264             
265         }
266    
267         if(debug) { 
268         	System.out.printf("%s: Final sense charge map:\n", this.getClass().getSimpleName());
269         	for(Map.Entry<Integer, Integer> entry :  _sense_data.get(ChargeCarrier.HOLE).getChargeMap().entrySet()) {
270         		System.out.printf("%s: cell %d -> %d \n", this.getClass().getSimpleName(),entry.getKey(),entry.getValue());
271         	}
272         }
273         
274     }
275     
276     private void transferChargeToReadout()
277     {
278         for (ChargeCarrier carrier : ChargeCarrier.values())
279         {
280             if (_sensor.hasElectrodesOnSide(carrier))
281             {
282                 if (_sensor.isACCoupled(carrier))
283                 {
284 //                    System.out.println("Is AC coupled");
285                     
286                     SiSensorElectrodes sense_electrodes = _sensor.getSenseElectrodes(carrier);
287                     SiSensorElectrodes readout_electrodes = _sensor.getReadoutElectrodes(carrier);
288                     BasicMatrix transfer_efficiencies = _sensor.getTransferEfficiencies(carrier);
289                     
290                     SiElectrodeDataCollection sense_data = _sense_data.get(carrier);
291                     SiElectrodeDataCollection readout_data = _readout_data.get(carrier);
292                     
293                     for (Integer sense_cell : sense_data.keySet())
294                     {
295                         
296 //                        System.out.println("Processing sense cell : "+sense_cell);
297                         
298                         SiElectrodeData sense_cell_data = sense_data.get(sense_cell);
299                         
300 //                        System.out.println("Sense cell charge : "+sense_cell_data.getCharge());
301                         
302                         int sense_row = sense_electrodes.getRowNumber(sense_cell);
303                         int sense_col = sense_electrodes.getColumnNumber(sense_cell);
304                         int row_steps = transfer_efficiencies.getNRows()-1;
305                         int col_steps = transfer_efficiencies.getNColumns()-1;
306                         
307 //                        System.out.println("sense_row : "+sense_row);
308 //                        System.out.println("sense_col : "+sense_col);
309 //                        System.out.println("row_steps : "+row_steps);
310 //                        System.out.println("col_steps : "+col_steps);
311                         
312 //                        System.out.println("transfer_efficiencies : "+transfer_efficiencies);
313 
314                         for (int irow = sense_row - row_steps; irow <= sense_row + row_steps; irow++)
315                         {
316                             
317 //                            System.out.println("irow : "+irow);
318                             
319                             for (int icol = sense_col - col_steps; icol <= sense_col + col_steps; icol++)
320                             {
321                                 
322 //                                System.out.println("icol : "+icol);
323                                 
324                                 int sense_id = sense_electrodes.getCellID(irow,icol);
325                                 if (sense_id < 0) continue;
326                                 Hep3Vector sense_position = sense_electrodes.getCellPosition(sense_id);
327                                 int readout_cell = readout_electrodes.getCellID(sense_position);
328                                 
329 //                                System.out.println("sense_id : "+sense_id);
330 //                                System.out.println("sense_position: "+sense_position);
331 //                                System.out.println("readout_cell : "+readout_cell);
332                                 
333 //                                System.out.println("position_in_cell : "+readout_electrodes.getPositionInCell(sense_position));
334 //
335 //                                System.out.println("Sense position : "+sense_electrodes.getCellPosition(sense_id));
336 //                                System.out.println("Readout position : "+readout_electrodes.getCellPosition(readout_cell));
337                                 
338                                 if ( readout_electrodes.isValidCell(readout_cell) &&
339                                         readout_electrodes.getPositionInCell(sense_position).x() == 0.0 &&
340                                         readout_electrodes.getPositionInCell(sense_position).y() == 0.0 )
341                                 {
342 //                                    System.out.println("transferring...");
343                                     double transfer_efficiency = transfer_efficiencies.e(Math.abs(irow-sense_row),Math.abs(icol-sense_col));
344 //                                    System.out.println("transfer efficiency: "+transfer_efficiency);
345 //                                    System.out.println("sense charge: "+sense_cell_data.getCharge());
346                                     
347                                     SiElectrodeData readout_datum  = new SiElectrodeData((int)Math.round(transfer_efficiency*sense_cell_data.getCharge()),
348                                             sense_cell_data.getSimulatedHits());
349                                     readout_data.add(readout_cell,readout_datum);
350 //                                    System.out.println("transferring 3...");
351 //                                    System.out.println("Current readout charge map: " + _readout_data.get(carrier).getChargeMap());
352                                 }
353 //                                System.out.println("4...");
354                             }
355 //                            System.out.println("5...");
356                         }
357 //                        System.out.println("6...");
358                     }
359 //                    System.out.println("7...");
360                 }
361                 else
362                 {
363                     _readout_data.put(carrier,_sense_data.get(carrier));
364                 }
365                 
366 //                System.out.println("Final readout charge map: " + _readout_data.get(carrier).getChargeMap());
367             }
368         }
369         
370         clearSense();
371         
372     }
373     
374     private int nSegments(TrackSegment track, ChargeCarrier carrier, double deposition_granularity)
375     {
376         // Decide how to cut track into pieces as a fraction of strip pitch
377         int nsegments = 0;
378         if (!_sensor.hasElectrodesOnSide(carrier)) return nsegments;
379         SiSensorElectrodes electrodes = _sensor.getSenseElectrodes(carrier);
380         
381 //        System.out.println("Track P1: " + track.getP1());
382 //        System.out.println("Track P2: " + track.getP2());
383 //        System.out.println("Drift Destination of P1: " + driftDestination(track.getP1(),carrier));
384 //        System.out.println("Drift Destination of P2: " + driftDestination(track.getP2(),carrier));
385         
386         nsegments = (int)Math.ceil(Math.abs(VecOp.dot(track.getVector(),electrodes.getGeometry().getNormal()))/(_sensor.getThickness()*deposition_granularity));
387         
388 //        nsegments = (int)Math.ceil(track.getVector().z()/(_sensor.getThickness()*deposition_granularity));  // old way
389         
390         Hep3Vector deposition_line = VecOp.sub( driftDestination(track.getP2(),carrier),
391                 driftDestination(track.getP1(),carrier) );
392         
393         int naxes = electrodes.getNAxes();
394         for (int iaxis = 0; iaxis < naxes; iaxis++)
395         {
396             IRotation3D electrode_to_sensor = electrodes.getParentToLocal().getRotation().inverse();
397             Hep3Vector measured_coordinate = electrode_to_sensor.rotated(electrodes.getMeasuredCoordinate(iaxis));
398             
399             double projected_deposition_length = Math.abs(VecOp.dot(deposition_line,measured_coordinate));
400             
401 //            double projected_deposition_length = Math.abs(VecOp.dot(deposition_line,_sensor.getMeasuredCoordinates(carrier)[iaxis]));
402             
403 //            System.out.println("Projected deposition Length: " + projected_deposition_length);
404             
405             int required_segments = (int)Math.ceil(projected_deposition_length/
406                     (deposition_granularity*electrodes.getPitch(iaxis)));
407             nsegments = Math.max(nsegments,required_segments);
408         }
409         
410 //        if (nsegments > 1000)
411 //        {
412 //            System.out.println("Track P1: " + track.getP1());
413 //            System.out.println("Track P2: " + track.getP2());
414 //            System.out.println("Drift Destination of P1: " + driftDestination(track.getP1(),carrier));
415 //            System.out.println("Drift Destination of P2: " + driftDestination(track.getP2(),carrier));
416 //            
417 //            System.out.println("nsegments: " + nsegments);
418 //        }
419         
420         return nsegments;
421     }
422     
423     
424     private Hep3Vector driftDestination(Hep3Vector origin, ChargeCarrier carrier)
425     {
426 //        System.out.println("Beginning driftDestination for origin: "+origin);
427         return VecOp.add(origin,driftVector(origin, carrier));
428     }
429     
430     private Hep3Vector driftVector(Hep3Vector origin, ChargeCarrier carrier)
431     {
432 //        System.out.println("Beginning driftVector");
433         
434         double drift_vector_scale = _sensor.distanceFromSide(origin,carrier)/VecOp.dot(_drift_direction.get(carrier),_sensor.getBiasSurface(carrier).getNormal());
435 //        double drift_vector_scale = _sensor.distanceFromSide(origin,carrier)/_drift_direction.get(carrier).z();
436         
437         return VecOp.mult(drift_vector_scale,_drift_direction.get(carrier));
438     }
439     
440     private Hep3Vector driftDirection(ChargeCarrier carrier, Hep3Vector local_position)
441     {
442 //        System.out.println("Beginning driftDirection");
443         
444 //        System.out.println("Position: "+local_position);
445         
446 //        System.out.println("Carrier: "+carrier);
447 
448         //  Get the magnetic field - already in Si units (Tesla)
449         Hep3Vector b_field = _sensor.getBField(local_position);
450 //        System.out.println("B field: "+b_field);
451 
452         //  Get the electric field and convert from V/mm to Si units (V/m)
453         Hep3Vector e_field = VecOp.mult(1000., _sensor.electricField(local_position));
454 //        System.out.println("E field: "+e_field);
455 
456         //  Get the mobility and convert from cm^2/V/s to Si units of m^2/V/s
457         double mobility = 1.0e-4 * _sensor.getBulk().mobility(carrier);
458 
459         //  Get the charge (to be consistent with mobility, this should be in units of e)
460         double qmu =  carrier.charge() * mobility;
461 //        System.out.println("mobility: "+mobility+" charge: "+qmu/mobility;
462         
463         //  Calculate the velocity parallel to the E field vpar = q * mu * E
464         Hep3Vector vpar = VecOp.mult(qmu, e_field);
465 //        System.out.println("vpar: "+vpar);
466 
467         //  Calculate the velocity perpendicular to the E field vperp = q^2 * mu^2 * E x B
468         Hep3Vector vperp = VecOp.mult(qmu*qmu, VecOp.cross(e_field, b_field));
469 //        System.out.println("vperp: "+vperp);
470 
471         //  Calculate a unit vector in the drift direction
472         Hep3Vector drift_direction = VecOp.unit(VecOp.add(vpar, vperp));
473         
474 //        System.out.println("Drift direction: "+drift_direction);
475         
476         return drift_direction;
477         
478     }
479     
480     private ChargeDistribution diffusionDistribution(double segment_charge, Hep3Vector origin, ChargeCarrier carrier)
481     {
482         
483 //        System.out.println("\n"+"Calculating charge distribution for carrier: "+carrier);
484         
485         // Local variables and a quick check
486         double distance = _sensor.distanceFromSide(origin,carrier);
487         double thickness = _sensor.getThickness();
488         
489 //        System.out.println("Distance from side: "+distance);
490         
491         
492         if (distance < -DISTANCE_ERROR_THRESHOLD || distance > thickness + DISTANCE_ERROR_THRESHOLD){
493             throw new RuntimeException("Distance is outside of sensor by more than "+DISTANCE_ERROR_THRESHOLD+". Distance = "+distance+
494                     ". If this is an isolated event, then perhaps DISTANCE_ERROR_THRESHOLD must be increased in CDFSiSensorSim"); 
495         }
496         else if (distance < 0) distance = 0.;
497         else if (distance > thickness) distance = thickness;
498 //        {
499 //            throw new RuntimeException("Attempting to drift charge from outside of sensor!");
500 //        }
501         double bias_voltage = _sensor.getBiasVoltage();
502         double depletion_voltage = _sensor.getDepletionVoltage();
503         
504         // Common factors
505         double difference_V = bias_voltage - depletion_voltage;
506         double sum_V = bias_voltage + depletion_voltage;
507         double common_factor = 2.0 * distance * depletion_voltage / thickness;
508         
509 //        System.out.println("Distance from side: "+_sensor.distanceFromSide(origin,carrier));
510 //        System.out.println("Origin: "+origin);
511         
512         // Calculate charge spreading without magnetic field
513         double sigmasq = _sensor.getBulk().K_BOLTZMANN * _sensor.getBulk().getTemperature() *
514                 _sensor.getThickness()*_sensor.getThickness() / _sensor.getDepletionVoltage();
515         if (_sensor.getBulk().isNtype() == (carrier==ChargeCarrier.HOLE))
516         {
517             sigmasq *= Math.log( sum_V / (sum_V - common_factor));
518         }
519         else
520         {
521             sigmasq *= Math.log( (difference_V + common_factor) / difference_V );
522         }
523         
524         double sigma = Math.sqrt(sigmasq);
525         
526 //        System.out.println("sigma: "+sigma);
527         
528         // Corrections for magnetic field -- this is an approximation, may have to be done better for high fields
529         
530         // Special case if field is parallel to drift direction
531         Hep3Vector drift_direction = _drift_direction.get(carrier);
532         Hep3Vector bias_surface_normal = _sensor.getBiasSurface(carrier).getNormal();
533         
534         Hep3Vector major_axis;
535         Hep3Vector minor_axis;
536         
537         double major_axis_length;
538         double minor_axis_length;
539         
540         if (VecOp.cross(drift_direction,bias_surface_normal).magnitude() < GeomOp3D.ANGULAR_TOLERANCE)
541         {
542 //            System.out.println("Drift direction perpendicular to bias plane");
543 //            System.out.println("Drift direction: "+drift_direction);
544 //            System.out.println("bias_surface_normal: "+bias_surface_normal);
545 //            System.out.println("First edge of bias plane: "+_sensor.getBiasSurface(carrier).getEdges().get(0).getDirection());
546             
547             
548             major_axis = VecOp.cross(bias_surface_normal,_sensor.getBiasSurface(carrier).getEdges().get(0).getDirection()); // arbitrary vector in plane
549             minor_axis = VecOp.cross(bias_surface_normal,major_axis);
550             
551 //            System.out.println("Major axis at initialization: "+major_axis);
552 //            System.out.println("Minor axis at initialization: "+minor_axis);
553             
554             major_axis_length = minor_axis_length = sigma;
555             
556         }
557         else
558         {
559 //            System.out.println("Drift direction NOT perpendicular to bias plane!!");
560             
561             major_axis = VecOp.unit(VecOp.cross(bias_surface_normal,VecOp.cross(_drift_direction.get(carrier),bias_surface_normal)));
562             minor_axis = VecOp.cross(bias_surface_normal,major_axis);
563             
564 //        double phi_lorentz = VecOp.phi(_drift_direction.get(carrier));  // FIXME: careful here... calculate axis directions to instantiate distribution.
565             // Project-to-plane would definitely be convenient here!!!
566             
567             double cos_theta_lorentz = VecOp.dot(_drift_direction.get(carrier),_sensor.getBiasSurface(carrier).getNormal());
568 //        double cos_theta_lorentz = VecOp.cosTheta(_drift_direction.get(carrier));                                 // FIXME: careful with cosTheta here!
569 //        System.out.println("Cos theta lorentz: "+cos_theta_lorentz);
570             
571             minor_axis_length = sigma*(1.0/cos_theta_lorentz); // drift time correction
572             major_axis_length = minor_axis_length*(1.0/cos_theta_lorentz); // + drift angle correction
573         }
574         
575         
576         major_axis = VecOp.mult(major_axis_length,major_axis);
577         minor_axis = VecOp.mult(minor_axis_length,minor_axis);
578         
579 //        System.out.println("Major axis: "+major_axis);
580 //        System.out.println("Minor axis: "+minor_axis);
581         
582         // FIXME: this has a Z component!!! (is that OK??  I think the whole thing transforms into the electrode coordinates before integrating charge.)
583         Hep3Vector drift_destination = driftDestination(origin,carrier);
584         
585 //        System.out.println("Drift destination: "+drift_destination);
586         
587         ChargeDistribution distribution = new GaussianDistribution2D(segment_charge, drift_destination,major_axis,minor_axis);
588         
589 //        ChargeDistribution distribution = new GaussianDistribution2D(segment_charge, drift_destination,new BasicHep3Vector(major_axis_length,0.0,0.0),new BasicHep3Vector(0.0,minor_axis_length,0.0));
590 //
591 //        ITransform3D phi_rotation = new Transform3D(new RotationPassiveXYZ(0.0,0.0,-phi_lorentz));
592 //        distribution.transform(phi_rotation);
593         
594 //        System.out.println("Done calculating charge distribution for carrier: "+carrier);
595         
596         return distribution;
597         
598     }
599     
600 }
601 
602 
603 
604 
605 
606 
607 
608 
609 
610 
611 
612 
613 
614 
615 
616 
617 
618