CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HitDigitizerFP420.cc
Go to the documentation of this file.
1 // File: HitDigitizerFP420.cc
3 // Date: 08.2008
4 // Description: HitDigitizerFP420 for FP420
5 // Modifications:
8 //#include "SimG4CMS/FP420/interface/FP420G4HitCollection.h"
9 //#include "SimG4CMS/FP420/interface/FP420G4Hit.h"
10 
15 
19 
20 using namespace std;
21 #include<vector>
22 
23 #define CBOLTZ (1.38E-23)
24 #define e_SI (1.6E-19)
25 
26 //HitDigitizerFP420::HitDigitizerFP420(float in,float inp,float inpx,float inpy){
27 //HitDigitizerFP420::HitDigitizerFP420(float in,float inp,float inpx,float inpy,float ild,float ildx,float ildy){
28 HitDigitizerFP420::HitDigitizerFP420(float in,float ild,float ildx,float ildy,float in0,float in2,float in3,int verbosity){
29  moduleThickness =in;
30  double bz420 = in0;
31  double bzD2 = in2;
32  double bzD3 = in3;
33  // double pitch =inp;
34  // double pitchX =inpx;
35  // double pitchY =inpy;
36  double ldrift =ild;
37  double ldriftX =ildx;
38  double ldriftY =ildy;
39  //
40  // Construct default classes
41  //
42 
43  //theCDividerFP420 = new ChargeDividerFP420(pitch);
44  theCDividerFP420 = new ChargeDividerFP420(moduleThickness, bz420, bzD2, bzD3, verbosity);
45 
46  depletionVoltage=20.0; //
47  appliedVoltage=25.0; // a bit bigger than depletionVoltage to have positive value A for logA, A=1-2*Tfract.*Vd/(Vd+Vb)
48 
49  // chargeMobility=480.0;// = 480.0 !holes mobility [cm**2/V/sec] p-side; = 1350.0 !electron mobility - n-side
50  chargeMobility=1350.0;// = 480.0 !holes mobility [cm**2/V/sec] p-side; = 1350.0 !electron mobility - n-side
51  //temperature=297.; // 24 C degree +273 = 297 ---->diffusion const for electrons= (1.38E-23/1.6E-19)*1350.0*297=34.6
52  // diffusion const for holes = 12.3 [cm**2/sec]
53  temperature=263.; // -10 C degree +273 = 263 ---->diffusion const for electrons= (1.38E-23/1.6E-19)*1350.0*263=30.6
54  double diffusionConstant = CBOLTZ/e_SI * chargeMobility * temperature;
55  // noDiffusion=true; // true if no Diffusion
56  noDiffusion=false; // false if Diffusion
57  if (noDiffusion) diffusionConstant *= 1.0e-3;
58 
59  chargeDistributionRMS=6.5e-10;
60 
61  // arbitrary:
62  tanLorentzAnglePerTesla = 0. ; // try =0.106 if B field exist
63 
64  // double timeNormalisation = pow(moduleThickness,2)/(2.*depletionVoltage*chargeMobility);
65  // double timeNormalisation = pow(pitch,2)/(2.*depletionVoltage*chargeMobility);
66 
67  double timeNormalisation = pow(ldrift,2)/(2.*depletionVoltage*chargeMobility);//
68  // double timeNormalisation = pow(ldrift,2)/(2.*depletionVoltage*chargeMobility);//
69  timeNormalisation = timeNormalisation*0.01; // because ldrift in [mm] but mu_e in [cm2/V/sec
70 
71  // double timeNormalisation = pow(pitch/2.,2)/(2.*depletionVoltage*chargeMobility);// i entered pitch as a distance between 2 read-out strips, not real distance between any (p or n) electrods which will be in 2 times less. But in expression for timeNormalisation the real distance of charge collection must be, so use pitch/2. !
72 
73 
74  double clusterWidth=5.;// was = 3
75  gevperelectron= 3.61e-09; // double GevPerElectron = 3.61e-09
76 
77  // GevPerElectron AZ:average deposited energy per e-h pair [keV]??? =0.0036
78  if(verbosity>0) {
79  std::cout << "HitDigitizerFP420: constructor ldrift= " << ldrift << std::endl;
80  std::cout << "ldriftY= " <<ldriftY << "ldriftX= " <<ldriftX << std::endl;
81  std::cout << "depletionVoltage" << depletionVoltage << "appliedVoltage" << appliedVoltage << "chargeMobility" << chargeMobility << "temperature" << temperature << "diffusionConstant" << diffusionConstant << "chargeDistributionRMS" << chargeDistributionRMS << "moduleThickness" << moduleThickness << "timeNormalisation" << timeNormalisation << "gevperelectron" << gevperelectron << std::endl;
82  }
83  //ndif
84 
85  theCDrifterFP420 =
86  new ChargeDrifterFP420(moduleThickness,
87  timeNormalisation,
88  diffusionConstant,
89  temperature,
90  chargeDistributionRMS,
91  depletionVoltage,
92  appliedVoltage,
93  ldriftX,
94  ldriftY, verbosity);
95  // pitchX,
96  // pitchY);
97 
98  theIChargeFP420 = new InduceChargeFP420(clusterWidth,gevperelectron);
99 }
100 
101 
103  delete theCDividerFP420;
104  delete theCDrifterFP420;
105  delete theIChargeFP420;
106 }
107 
108 
109 //HitDigitizerFP420::hit_map_type HitDigitizerFP420::processHit(const PSimHit& hit, G4ThreeVector bfield, int xytype,int numStrips, double pitch){
110 HitDigitizerFP420::hit_map_type HitDigitizerFP420::processHit(const PSimHit& hit, const G4ThreeVector& bfield, int xytype,int numStrips, double pitch, int numStripsW, double pitchW, double moduleThickness, int verbosity){
111 
112  // use chargePosition just for cross-check in "induce" method
113  // hit center in 3D-detector r.f.
114 
115  float middlex = (hit.exitPoint().x() + hit.entryPoint().x() )/2.;
116  float middley = (hit.exitPoint().y() + hit.entryPoint().y() )/2.;
117 
118 
119  float chargePosition= -100.;
120  // Y:
121  if(xytype == 1) {
122  // chargePosition = fabs(-numStrips/2. - ( int(middle.x()/pitch) +1.) );
123  //chargePosition = fabs(int(middle.x()/pitch+0.5*(numStrips+1)) + 1.);
124  // chargePosition = fabs(int(middle.y()/pitch+0.5*(numStrips+1)) + 1.);
125  // local and global reference frames are rotated in 90 degree, so global X and local Y are collinear
126  // chargePosition = int(fabs(middle.x()/pitch + 0.5*numStrips + 1.));// charge in strip coord
127  chargePosition = 0.5*(numStrips) + middlex/pitch ;// charge in strip coord 0 - numStrips-1
128 
129 
130  }
131  // X:
132  else if(xytype == 2) {
133  // chargePosition = fabs(-numStrips/2. - ( int(middle.y()/pitch) +1.) );
134  //chargePosition = fabs(int(middle.y()/pitch+0.5*(numStrips+1)) + 1.);
135  // chargePosition = fabs(int(middle.x()/pitch+0.5*(numStrips+1)) + 1.);
136  // local and global reference frames are rotated in 90 degree, so global X and local Y are collinear
137  // chargePosition = int(fabs(middle.y()/pitch + 0.5*numStrips + 1.));
138  chargePosition = 0.5*(numStrips) + middley/pitch ;// charge in strip coord 0 - numStrips-1
139 
140  // std::cout << " chargePosition SiHitDi... = " << chargePosition << std::endl;
141  }
142  else {
143  std::cout <<"================================================================"<<std::endl;
144  std::cout << "**** HitDigitizerFP420: !!! ERROR: you have not to be here !!! xytype=" << xytype << std::endl;
145  // std::cout << "**** HitDigitizerFP420: !!! ERROR: you have not to be here !!! xytype=" << xytype << std::endl;
146 
147 
148  // break;
149  }
150  // if(chargePosition > numStrips || chargePosition<1) {
151  if(chargePosition > numStrips || chargePosition < 0) {
152  std::cout << "**** HitDigitizerFP420: !!! ERROR: check correspondence of XY detector dimensions in XML and here !!! chargePosition = " << chargePosition << std::endl;
153  // break;
154  }
155 
156  if(verbosity>0) {
157  std::cout << " ====== **** HitDigitizerFP420: !!! processHit !!! : input: xytype=" << xytype << " numStrips= " << numStrips << " pitch= " << pitch << " Calculated chargePosition= " << chargePosition << std::endl;
158  std::cout << "The middle of hit point on input was: middlex = " << middlex << std::endl;
159  std::cout << "The middle of hit point on input was: middley = " << middley << std::endl;
160  // std::cout << "For checks: hit point Entry = " << hit.getEntry() << std::endl;
161  std::cout << " ====== **** HitDigitizerFP420:processHit: start CDrifterFP420 divide" << std::endl;
162  }
163  //
164  // Fully process one SimHit
165  //
166 
167  // CDrifterFP420::ionization_type ion = theCDividerFP420->divide(hit,pitch);
168  CDrifterFP420::ionization_type ion = theCDividerFP420->divide(hit,moduleThickness);
169  //
170  // Compute the drift direction for this det
171  //
172 
173  // G4ThreeVector driftDir = DriftDirection(&det,bfield);
174  G4ThreeVector driftDir = DriftDirection(bfield,xytype,verbosity);
175  if(verbosity>0) {
176  std::cout << " ====== **** HitDigitizerFP420:processHit: driftDir= " << driftDir << std::endl;
177  std::cout << " ====== **** HitDigitizerFP420:processHit: start induce , CDrifterFP420 drift " << std::endl;
178  }
179 
180  // if(driftDir.z() ==0.) {
181  // std::cout << " pxlx: drift in z is zero " << std::endl;
182  // } else
183  //
184 
185  return theIChargeFP420->induce(theCDrifterFP420->drift(ion,driftDir,xytype), numStrips, pitch, numStripsW, pitchW, xytype, verbosity);
186 
187  //
188 }
189 
190 
191 
192 G4ThreeVector HitDigitizerFP420::DriftDirection(const G4ThreeVector& _bfield, int xytype, int verbosity){
193 
194  // LOCAL hit: exchange xytype: 1 <-> 2
195 
196  // Frame detFrame(_detp->surface().position(),_detp->surface().rotation());
197  // G4ThreeVector Bfield=detFrame.toLocal(_bfield);
198  G4ThreeVector Bfield=_bfield;
199  float dir_x, dir_y, dir_z;
200  // this lines with dir_... have to be specified(sign?) in dependence of field direction:
201  /*
202  float dir_x = tanLorentzAnglePerTesla * Bfield.y();
203  float dir_y = -tanLorentzAnglePerTesla * Bfield.x();
204  float dir_z = 1.; // if E field is in z direction
205  */
206  // for electrons:
207  // E field is either in X or Y direction with change vector to opposite
208 
209  // global Y or localX
210  // E field is in Xlocal direction with change vector to opposite
211  if(xytype == 2) {
212  dir_x = 1.; // E field in Xlocal direction
213  dir_y = +tanLorentzAnglePerTesla * Bfield.z();
214  dir_z = -tanLorentzAnglePerTesla * Bfield.y();
215  }
216  // global X
217  // E field is in Ylocal direction with change vector to opposite
218  else if(xytype == 1) {
219  dir_x = +tanLorentzAnglePerTesla * Bfield.z();
220  dir_y = 1.; // E field in Ylocal direction
221  dir_z = -tanLorentzAnglePerTesla * Bfield.x();
222  }
223  else{
224  dir_x = 0.;
225  dir_y = 0.;
226  dir_z = 0.;
227  std::cout << "HitDigitizerFP420: ERROR - wrong xytype=" << xytype << std::endl;
228  }
229 
230 
231  // G4ThreeVector theDriftDirection = LocalVector(dir_x,dir_y,dir_z);
232  G4ThreeVector theDriftDirection(dir_x,dir_y,dir_z);
233  // Local3DPoint EntryPo(aHit->getEntry().x(),aHit->getEntry().y(),aHit->getEntry().z());
234  if(verbosity>0) {
235  std::cout << "HitDigitizerFP420:DriftDirection tanLorentzAnglePerTesla= " << tanLorentzAnglePerTesla << std::endl;
236  std::cout << "HitDigitizerFP420:DriftDirection The drift direction in local coordinate is " << theDriftDirection << std::endl;
237  }
238 
239  return theDriftDirection;
240 
241 }
#define CBOLTZ
std::map< int, float, std::less< int > > hit_map_type
G4ThreeVector DriftDirection(const G4ThreeVector &, int, int)
T y() const
Definition: PV3DBase.h:63
#define e_SI
HitDigitizerFP420(float in, float ild, float ildx, float ildy, float in0, float in2, float in3, int verbosity)
Local3DPoint exitPoint() const
Exit point in the local Det frame.
Definition: PSimHit.h:38
std::vector< EnergySegmentFP420 > ionization_type
Definition: CDrifterFP420.h:15
tuple cout
Definition: gather_cfg.py:121
T x() const
Definition: PV3DBase.h:62
Local3DPoint entryPoint() const
Entry point in the local Det frame.
Definition: PSimHit.h:35
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
hit_map_type processHit(const PSimHit &, const G4ThreeVector &, int, int, double, int, double, double, int)