All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
SiEnergyFluct.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * DISCLAIMER *
4 // * *
5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * of contributors to this software. The specific disclaimers,which *
7 // * govern, are listed with their locations in: *
8 // * http://cern.ch/geant4/license *
9 // * *
10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assume any liability for its *
14 // * use. *
15 // * *
16 // * This code implementation is the intellectual property of the *
17 // * GEANT4 collaboration. *
18 // * By copying, distributing or modifying the Program (or any work *
19 // * based on the Program) you indicate your acceptance of this *
20 // * statement, and all its terms. *
21 // ********************************************************************
22 //
23 // -------------------------------------------------------------------
24 //
25 // This product includes software developed by Members of the Geant4
26 // Collaboration ( http://cern.ch/geant4 ).
27 //
28 // GEANT4 Classes utilized:
29 //
30 // File names: G4UniversalFluctuation (by Vladimir Ivanchenko)
31 // G4MollerBhabhaModel (by Vladimir Ivanchenko)
32 // G4MuBetheBlochModel (by Vladimir Ivanchenko)
33 // G4BetheBlochModel (by Vladimir Ivanchenko)
34 //
35 // -------------------------------------------------------------------
36 
37 #include "SiEnergyFluct.h"
38 #include "PhysicalConstants.h"
39 #include <algorithm>
40 #include <cmath>
41 
42 // Include CLHEP classes
43 #include "CLHEP/Random/RandGaussQ.h"
44 #include "CLHEP/Random/RandPoisson.h"
45 #include "CLHEP/Random/RandFlat.h"
46 
47 // Include LCIO classes
48 #include <EVENT/MCParticle.h>
49 
50 // Namespaces
51 using namespace CLHEP ;
52 
53 namespace sistrip {
54 
55 //
56 // Constructor setups various constants for dE/dx calculations & universal fluctuations (for Si only)
57 //
58 SiEnergyFluct::SiEnergyFluct(double cutOnDeltaRays) : _cutOnDeltaRays(cutOnDeltaRays)
59 {
60  // Variables used to reduce recalculation of mean dE/dx
61  _prevMCPart = 0;
62  _prevMeanLoss = 0.;
63 
64  // Cut on secondary electrons
65  _cutOnDeltaRays = cutOnDeltaRays;
66 
67  // Constants related to dE/dx
68  _twoln10 = 2.0*log(10.0);
69 
70  // Constants related to Si material & dE/dx
71  _eexc = 173 * eV; // material->GetIonisation()->GetMeanExcitationEnergy();
72  _eexc2 = _eexc*_eexc;
73  _Zeff = 14; // material->GetElectronDensity()/material->GetTotNbOfAtomsPerVolume();
74  _th = 0.935414 * keV; // 0.25*sqrt(Zeff)
75 
76  // Constants related to Si material & dEdx -> density effect
77  _aden = 0.160077; // material->GetIonisation()->GetAdensity();
78  _cden = 4.43505; // material->GetIonisation()->GetCdensity();
79  _mden = 3; // material->GetIonisation()->GetMdensity();
80  _x0den = 0.2; // material->GetIonisation()->GetX0density();
81  _x1den = 3; // material->GetIonisation()->GetX1density();
82  _KBethe = 0.178328 * MeV/cm; // material->GetElectronDensity()*twopi_mc2_rcl2;
83 
84  // Constants related to dEdx in Si material & muons
85  _xgi[0]=0.0199;_xgi[1]=0.1017;_xgi[2]=0.2372;_xgi[3]=0.4083;_xgi[4]=0.5917;_xgi[5]=0.7628;_xgi[6]=0.8983;_xgi[7]=0.9801;
86  _wgi[0]=0.0506;_wgi[1]=0.1112;_wgi[2]=0.1569;_wgi[3]=0.1813;_wgi[4]=0.1813;_wgi[5]=0.1569;_wgi[6]=0.1112;_wgi[7]=0.0506;
87 
88  _limitKinEnergy = 100.* keV;
89  _logLimitKinEnergy = log(_limitKinEnergy);
90  _alphaPrime = fine_str_const/2./pi;
91 
92  // Constants related to universal fluctuations
93  _minLoss = 10 * eV;
94  _minNumberInteractionsBohr = 10.0;
95  _nmaxCont1 = 4.;
96  _nmaxCont2 = 16.;
97 
98  _facwidth = 1. / keV;
99  _f1Fluct = 0.857143; // material->GetIonisation()->GetF1fluct();
100  _f2Fluct = 0.142857; // material->GetIonisation()->GetF2fluct();
101  _e1Fluct = 0.115437 * keV;// material->GetIonisation()->GetEnergy1fluct();
102  _e2Fluct = 1.96 * keV; // material->GetIonisation()->GetEnergy2fluct();
103  _e1LogFluct = log(_e1Fluct);
104  _e2LogFluct = log(_e2Fluct);
105  _ipotFluct = 173. * eV; // material->GetIonisation()->GetLogMeanExcEnergy();
106  _ipotLogFluct = log(_ipotFluct);
107  _e0 = 10 * eV; // material->GetIonisation()->GetEnergy0fluct();
108 }
109 
110 //
111 // Destructor
112 //
113 SiEnergyFluct::~SiEnergyFluct() {}
114 
115 //
116 // Main method providing energy loss fluctuations, the model used to get the
117 // fluctuations is essentially the same as in Glandz in Geant3 (Cern program
118 // library W5013, phys332). L. Urban et al. NIM A362, p.416 (1995) and Geant4
119 // Physics Reference Manual
120 //
121 double SiEnergyFluct::SampleFluctuations(const MCParticle * part, const double length)
122 {
123  // Calculate particle related quantities
124  double mass = part->getMass() * GeV;
125  double momentum2 = ( part->getMomentum()[0]*part->getMomentum()[0] +
126  part->getMomentum()[1]*part->getMomentum()[1] +
127  part->getMomentum()[2]*part->getMomentum()[2] ) * GeV*GeV;
128  double kineticEnergy = sqrt(momentum2 + mass*mass) - mass;
129  double ratio = e_mass/mass;
130  double tau = kineticEnergy/mass;
131  double gam = tau + 1.0;
132  double gam2 = gam*gam;
133  double bg2 = tau * (tau+2.0);
134  double beta2 = bg2/(gam2);
135  int pdg = part->getPDG();
136  double chargeSquare = part->getCharge()*part->getCharge();
137 
138  double maxT = 2.0*e_mass*tau*(tau + 2.) / (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
139  maxT = std::min(_cutOnDeltaRays,maxT);
140 
141  // Recalculate mean loss - mean dE/dx, if necessary
142  double meanLoss = 0.;
143 
144  if (part != _prevMCPart) {
145 
146  // Electron or positron
147  if (pdg == 11) meanLoss = getElectronDEDX(part);
148 
149  // Muon
150  else if (pdg == 13) meanLoss = getMuonDEDX(part);
151 
152  // Hadron
153  else meanLoss = getHadronDEDX(part);
154  }
155  // not needed - same particle
156  else {
157 
158  meanLoss = _prevMeanLoss;
159  }
160 
161  // Save values for next call
162  _prevMCPart = part;
163  _prevMeanLoss = meanLoss;
164 
165  // Get mean loss in MeV
166  meanLoss *= length;
167 
168  //
169  // Start calculation
170  double loss = 0.;
171  double sigma = 0.;
172 
173  // Trick for very very small loss ( out of model validity)
174  if (meanLoss < _minLoss)
175  {
176  return meanLoss;
177  }
178 
179  //
180  // Gaussian regime for heavy particles only
181  if ((mass > e_mass) && (meanLoss >= _minNumberInteractionsBohr*maxT)) {
182 
183  double massrate = e_mass/mass ;
184  double tmaxkine = 2.*e_mass*beta2*gam2/(1.+massrate*(2.*gam+massrate));
185 
186  if (tmaxkine <= 2.*maxT) {
187 
188  // Sigma
189  sigma = (1.0/beta2 - 0.5) * _KBethe * maxT * length * chargeSquare;
190  sigma = sqrt(sigma);
191 
192  double twomeanLoss = meanLoss + meanLoss;
193 
194  if (twomeanLoss < sigma) {
195 
196  double x;
197  do {
198 
199  loss = twomeanLoss*RandFlat::shoot();
200  x = (loss - meanLoss)/sigma;
201 
202  } while (1.0 - 0.5*x*x < RandFlat::shoot());
203 
204  } else {
205 
206  do {
207 
208  loss = RandGaussQ::shoot(meanLoss,sigma);
209 
210  } while (loss < 0. || loss > twomeanLoss);
211 
212  }
213  return loss;
214  }
215  }
216 
217  //
218  // Glandz regime : initialisation
219 
220  // Trick for very small step or low-density material
221  if (maxT <= _e0) return meanLoss;
222 
223  double a1 = 0.;
224  double a2 = 0.;
225  double a3 = 0.;
226 
227  // Correction to get better width even using stepmax
228  // if(abs(meanLoss- oldloss) < 1.*eV)
229  // samestep += 1;
230  // else
231  // samestep = 1.;
232  // oldloss = meanLoss;
233 
234  double width = 1.+_facwidth*meanLoss;
235  if (width > 4.50) width = 4.50;
236  double e1 = width*_e1Fluct;
237  double e2 = width*_e2Fluct;
238 
239  // Cut and material dependent rate
240  double rate = 1.0;
241  if (maxT > _ipotFluct) {
242 
243  double w2 = log(2.*e_mass*beta2*gam2)-beta2;
244 
245  if (w2 > _ipotLogFluct && w2 > _e2LogFluct) {
246 
247  rate = 0.03+0.23*log(log(maxT/_ipotFluct));
248  double C = meanLoss*(1.-rate)/(w2-_ipotLogFluct);
249 
250  a1 = C*_f1Fluct*(w2 -_e1LogFluct)/e1;
251  a2 = C*_f2Fluct*(w2 -_e2LogFluct)/e2;
252 
253  }
254  }
255 
256  double w1 = maxT/_e0;
257  if(maxT > _e0) a3 = rate*meanLoss*(maxT-_e0)/(_e0*maxT*log(w1));
258 
259  // 'Nearly' Gaussian fluctuation if a1>nmaxCont2&&a2>nmaxCont2&&a3>nmaxCont2
260  double emean = 0.;
261  double sig2e = 0.;
262  double sige = 0.;
263  double p1 = 0.;
264  double p2 = 0.;
265  double p3 = 0.;
266 
267  // Excitation of type 1
268  if (a1 > _nmaxCont2) {
269 
270  emean += a1*e1;
271  sig2e += a1*e1*e1;
272  }
273  else if (a1 > 0.) {
274 
275  p1 = double(RandPoisson::shoot(a1));
276  loss += p1*e1;
277 
278  if(p1 > 0.) loss += (1.-2.*RandFlat::shoot())*e1;
279  }
280 
281  // Excitation of type 2
282  if (a2 > _nmaxCont2) {
283 
284  emean += a2*e2;
285  sig2e += a2*e2*e2;
286  }
287  else if (a2 > 0.) {
288 
289  p2 = double(RandPoisson::shoot(a2));
290  loss += p2*e2;
291 
292  if(p2 > 0.) loss += (1.-2.*RandFlat::shoot())*e2;
293  }
294 
295  // Ionisation
296  double lossc = 0.;
297 
298  if (a3 > 0.) {
299 
300  p3 = a3;
301  double alfa = 1.;
302 
303  if (a3 > _nmaxCont2) {
304 
305  alfa = w1*(_nmaxCont2+a3)/(w1*_nmaxCont2+a3);
306  double alfa1 = alfa*log(alfa)/(alfa-1.);
307  double namean = a3*w1*(alfa-1.)/((w1-1.)*alfa);
308  emean += namean*_e0*alfa1;
309  sig2e += _e0*_e0*namean*(alfa-alfa1*alfa1);
310  p3 = a3-namean;
311  }
312 
313  double w2 = alfa*_e0;
314  double w = (maxT-w2)/maxT;
315 
316  int nb = RandPoisson::shoot(p3);
317 
318  if (nb > 0) for (int k=0; k<nb; k++) lossc += w2/(1.-w*RandFlat::shoot());
319  }
320 
321  if (emean > 0.) {
322 
323  sige = sqrt(sig2e);
324  loss += std::max(0.,RandGaussQ::shoot(emean,sige));
325  }
326 
327  loss += lossc;
328 
329  return loss;
330 }
331 
332 //
333 // Method calculating actual dE/dx for hadrons
334 //
335 double SiEnergyFluct::getHadronDEDX(const MCParticle * part)
336 {
337  // Calculate particle related quantities
338  double mass = part->getMass() * GeV;
339  double momentum2 = ( part->getMomentum()[0]*part->getMomentum()[0] +
340  part->getMomentum()[1]*part->getMomentum()[1] +
341  part->getMomentum()[2]*part->getMomentum()[2] ) * GeV*GeV;
342  double kineticEnergy = sqrt(momentum2 + mass*mass) - mass;
343  double spin = 0.5;
344  double charge2 = part->getCharge()*part->getCharge();
345  double ratio = e_mass/mass;
346  double tau = kineticEnergy/mass;
347  double gam = tau + 1.0;
348  double bg2 = tau * (tau+2.0);
349  double beta2 = bg2/(gam*gam);
350 
351  double maxT = 2.0*e_mass*tau*(tau + 2.) / (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
352  double cutEnergy = std::min(_cutOnDeltaRays,maxT);
353 
354  // Start with dE/dx
355  double dedx = log(2.0*e_mass*bg2*cutEnergy/_eexc2) - (1.0 + cutEnergy/maxT)*beta2;
356 
357  // Spin 0.5 particles
358  if(0.5 == spin) {
359 
360  double del = 0.5*cutEnergy/(kineticEnergy + mass);
361  dedx += del*del;
362  }
363 
364  // Density correction
365  double x = log(bg2)/_twoln10;
366 
367  if ( x >= _x0den ) {
368 
369  dedx -= _twoln10*x - _cden;
370  if ( x < _x1den ) dedx -= _aden*pow((_x1den-x),_mden) ;
371  }
372 
373  // Shell correction --> not used
374  //dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
375 
376  // Now compute the total ionization loss
377  if (dedx < 0.0) dedx = 0.0 ;
378 
379  dedx *= _KBethe*charge2/beta2;
380 
381  // High order correction only for hadrons --> not used
382  //if(!isIon) dedx += corr->HighOrderCorrections(p,material,kineticEnergy);
383 
384  return dedx;
385 }
386 
387 //
388 // Method calculating actual dE/dx for muons
389 //
390 double SiEnergyFluct::getMuonDEDX(const MCParticle * part)
391 {
392  // Calculate particle related quantities
393  double mass = part->getMass() * GeV;
394  double momentum2 = ( part->getMomentum()[0]*part->getMomentum()[0] +
395  part->getMomentum()[1]*part->getMomentum()[1] +
396  part->getMomentum()[2]*part->getMomentum()[2] ) * GeV*GeV;
397  double kineticEnergy = sqrt(momentum2 + mass*mass) - mass;
398  double totEnergy = kineticEnergy + mass;
399  double ratio = e_mass/mass;
400  double tau = kineticEnergy/mass;
401  double gam = tau + 1.0;
402  double bg2 = tau * (tau+2.0);
403  double beta2 = bg2/(gam*gam);
404 
405  double maxT = 2.0*e_mass*tau*(tau + 2.) / (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
406  double cutEnergy = std::min(_cutOnDeltaRays,maxT);
407 
408  // Start with dE/dx
409  double dedx = log(2.0*e_mass*bg2*cutEnergy/_eexc2) - (1.0 + cutEnergy/maxT)*beta2;
410 
411  double del = 0.5*cutEnergy/totEnergy;
412  dedx += del*del;
413 
414  // Density correction
415  double x = log(bg2)/_twoln10;
416 
417  if ( x >= _x0den ) {
418 
419  dedx -= _twoln10*x - _cden ;
420  if ( x < _x1den ) dedx -= _aden*pow((_x1den-x),_mden) ;
421  }
422 
423  // Shell correction --> not used
424  //dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
425 
426 
427  // Now compute the total ionization loss
428  if (dedx < 0.0) dedx = 0.0 ;
429 
430  // Use radiative corrections of R. Kokoulin
431  if (cutEnergy > _limitKinEnergy) {
432 
433  double logtmax = log(cutEnergy);
434  double logstep = logtmax - _logLimitKinEnergy;
435  double dloss = 0.0;
436  double ftot2 = 0.5/(totEnergy*totEnergy);
437 
438  for (int ll=0; ll<8; ll++)
439  {
440  double ep = exp(_logLimitKinEnergy + _xgi[ll]*logstep);
441  double a1 = log(1.0 + 2.0*ep/e_mass);
442  double a3 = log(4.0*totEnergy*(totEnergy - ep)/mass/mass);
443  dloss += _wgi[ll]*(1.0 - beta2*ep/maxT + ep*ep*ftot2)*a1*(a3 - a1);
444  }
445  dedx += dloss*logstep*_alphaPrime;
446  }
447 
448  dedx *= _KBethe/beta2;
449 
450  //High order corrections --> not used
451  //dedx += corr->HighOrderCorrections(p,material,kineticEnergy);
452 
453  return dedx;
454 }
455 
456 //
457 // Method calculating actual dE/dx for electrons & positrons
458 //
459 double SiEnergyFluct::getElectronDEDX(const MCParticle * part)
460 {
461  // Calculate particle related quantities
462  double mass = part->getMass() * GeV;
463  double momentum = ( part->getMomentum()[0]*part->getMomentum()[0] +
464  part->getMomentum()[1]*part->getMomentum()[1] +
465  part->getMomentum()[2]*part->getMomentum()[2] ) * GeV*GeV;
466  double kineticEnergy = sqrt(momentum*momentum + mass*mass) - mass;
467  double charge = part->getCharge();
468  double tau = kineticEnergy/mass;
469  double gam = tau + 1.0;
470  double gamma2 = gam*gam;
471  double bg2 = tau * (tau+2.0);
472  double beta2 = bg2/(gamma2);
473 
474  // Calculate the dE/dx due to the ionization by Seltzer-Berger formula
475  bool lowEnergy = false;
476  double tkin = kineticEnergy;
477 
478  if (kineticEnergy < _th) {
479  tkin = _th;
480  lowEnergy = true;
481  }
482  double lowLimit = 0.2*keV;
483 
484  double maxT = kineticEnergy;
485 
486  if(charge < 0.) maxT *= 0.5;
487 
488  double eexc = _eexc/e_mass;
489  double eexc2 = eexc*eexc;
490 
491  double d = std::min(_cutOnDeltaRays, maxT)/e_mass;
492  double dedx;
493 
494  // Electron
495  if (charge < 0.) {
496 
497  dedx = log(2.0*(tau + 2.0)/eexc2) - 1.0 - beta2 + log((tau-d)*d) + tau/(tau-d)
498  + (0.5*d*d + (2.0*tau + 1.)*log(1. - d/tau))/gamma2;
499  }
500  // Positron
501  else {
502 
503  double d2 = d*d*0.5;
504  double d3 = d2*d/1.5;
505  double d4 = d3*d*3.75;
506  double y = 1.0/(1.0 + gam);
507  dedx = log(2.0*(tau + 2.0)/eexc2) + log(tau*d) - beta2*(tau + 2.0*d - y*(3.0*d2
508  + y*(d - d3 + y*(d2 - tau*d3 + d4))))/tau;
509  }
510 
511  // Density correction
512  double x = log(bg2)/_twoln10;
513 
514  if (x >= _x0den) {
515 
516  dedx -= _twoln10*x - _cden;
517  if (x < _x1den) dedx -= _aden*pow(_x1den-x, _mden);
518  }
519 
520  // Now compute the total ionization loss
521  dedx *= _KBethe/beta2;
522  if (dedx < 0.0) dedx = 0.0;
523 
524  // Lowenergy extrapolation
525 
526  if (lowEnergy) {
527 
528  if (kineticEnergy >= lowLimit) dedx *= sqrt(tkin/kineticEnergy);
529  else dedx *= sqrt(tkin*kineticEnergy)/lowLimit;
530 
531  }
532 
533  return dedx;
534 }
535 
536 } // Namespace
sistrip::SiEnergyFluct::getElectronDEDX
double getElectronDEDX(const MCParticle *part)
Method calculating actual dEdx for electrons &amp; positrons - based on ComputeDEDXPerVolume method G4Mol...
Definition: SiEnergyFluct.cc:459
sistrip::MeV
static const float MeV
Definition: PhysicalConstants.h:43
sistrip::k
static const float k
Definition: PhysicalConstants.h:73
sistrip::C
static const float C
Definition: PhysicalConstants.h:31
sistrip::keV
static const float keV
Definition: PhysicalConstants.h:42
sistrip::e_mass
static const double e_mass
Definition: PhysicalConstants.h:79
sistrip::SiEnergyFluct::_cutOnDeltaRays
double _cutOnDeltaRays
Cut on secondary electrons - must be the same as in Geant4.
Definition: SiEnergyFluct.h:100
sistrip::SiEnergyFluct::SampleFluctuations
double SampleFluctuations(const MCParticle *part, const double length)
Method providing energy loss fluctuations, the model used to get the fluctuations is essentially the ...
Definition: SiEnergyFluct.cc:121
sistrip::SiEnergyFluct::_prevMCPart
const MCParticle * _prevMCPart
Definition: SiEnergyFluct.h:94
sistrip::SiEnergyFluct::getMuonDEDX
double getMuonDEDX(const MCParticle *part)
Method calculating actual dEdx for muons - based on ComputeDEDXPerVolume method from G4MuBetheBlochMo...
Definition: SiEnergyFluct.cc:390
sistrip::GeV
static const float GeV
Definition: PhysicalConstants.h:44
sistrip::eV
static const float eV
Definition: PhysicalConstants.h:41
sistrip::cm
static const float cm
Definition: PhysicalConstants.h:35
sistrip::pi
static const double pi
Definition: PhysicalConstants.h:69
sistrip::SiEnergyFluct::~SiEnergyFluct
~SiEnergyFluct()
Destructor.
Definition: SiEnergyFluct.cc:113
sistrip::SiEnergyFluct::getHadronDEDX
double getHadronDEDX(const MCParticle *part)
Method calculating actual dEdx for hadrons - based on ComputeDEDXPerVolume method from G4BetheBlochMo...
Definition: SiEnergyFluct.cc:335
sistrip::fine_str_const
static const double fine_str_const
Definition: PhysicalConstants.h:80