LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Protected Member Functions | List of all members
genf::GFEnergyLossBetheBloch Class Reference

#include "GFEnergyLossBetheBloch.h"

Inheritance diagram for genf::GFEnergyLossBetheBloch:
genf::GFAbsEnergyLoss

Public Member Functions

double energyLoss (const double &step, const double &mom, const int &pdg, const double &matDensity, const double &matZ, const double &matA, const double &radiationLength, const double &meanExcitationEnergy, const bool &doNoise=false, TMatrixT< Double_t > *noise=NULL, const TMatrixT< Double_t > *jacobian=NULL, const TVector3 *directionBefore=NULL, const TVector3 *directionAfter=NULL)
 Returns energy loss, optional calculation of energy loss straggeling. More...
 
virtual ~GFEnergyLossBetheBloch ()
 

Protected Member Functions

void getParticleParameters (const int &pdg, double &charge, double &mass)
 Gets particle charge and mass (in GeV) More...
 
double getParticleMass (const int &pdg)
 Returns particle mass (in GeV) More...
 

Detailed Description

Definition at line 45 of file GFEnergyLossBetheBloch.h.

Constructor & Destructor Documentation

genf::GFEnergyLossBetheBloch::~GFEnergyLossBetheBloch ( )
virtual

Definition at line 24 of file GFEnergyLossBetheBloch.cxx.

24  {
25 }

Member Function Documentation

double genf::GFEnergyLossBetheBloch::energyLoss ( const double &  step,
const double &  mom,
const int &  pdg,
const double &  matDensity,
const double &  matZ,
const double &  matA,
const double &  radiationLength,
const double &  meanExcitationEnergy,
const bool &  doNoise = false,
TMatrixT< Double_t > *  noise = NULL,
const TMatrixT< Double_t > *  jacobian = NULL,
const TVector3 *  directionBefore = NULL,
const TVector3 *  directionAfter = NULL 
)
virtual

Returns energy loss, optional calculation of energy loss straggeling.

Uses Bethe Bloch formula to calculate energy loss. For the energy loss straggeling, different formulas are used for different regions:

  • Vavilov-Gaussian regime
  • Urban/Landau approximation
  • truncated Landau distribution
  • Urban model

Implements genf::GFAbsEnergyLoss.

Definition at line 27 of file GFEnergyLossBetheBloch.cxx.

References beta, E, larg4::f1, larg4::f2, genf::GFAbsEnergyLoss::getParticleMass(), genf::GFAbsEnergyLoss::getParticleParameters(), and GFException::setFatal().

39  {
40 
41  static const double me = getParticleMass(11); // electron mass (GeV) 0.000519...
42 
43  double charge, mass;
44  getParticleParameters(pdg, charge, mass);
45 
46  const double beta = mom/sqrt(mass*mass+mom*mom);
47  double dedx = 0.307075*matZ/matA*matDensity/(beta*beta)*charge*charge;
48 
49  //for numerical stability
50  double gammaSquare = 1.-beta*beta;
51  if(gammaSquare>1.E-10) gammaSquare = 1./gammaSquare;
52  else gammaSquare = 1.E10;
53  double gamma = sqrt(gammaSquare);
54 
55  double massRatio = me/mass;
56  double argument = gammaSquare*beta*beta*me*1.E3*2./((1.E-6*meanExcitationEnergy) * sqrt(1+2*sqrt(gammaSquare)*massRatio + massRatio*massRatio));
57  if (argument <= exp(beta*beta))
58  dedx = 0.;
59  else{
60  dedx *= (log(argument)-beta*beta); // Bethe-Bloch [MeV/cm]
61  dedx *= 1.E-3; // in GeV/cm, hence 1.e-3
62  if (dedx <= 0.)
63  throw GFException("GFEnergyLossBetheBloch::energyLoss(): non-positive dE/dx", __LINE__, __FILE__).setFatal();
64  }
65 
66  double DE = step * dedx; //always positive
67  double momLoss = sqrt(mom*mom+2.*sqrt(mom*mom+mass*mass)*DE+DE*DE) - mom; //always positive
68 
69  //in vacuum it can numerically happen that momLoss becomes a small negative number. A cut-off at 0.01 eV for momentum loss seems reasonable
70  if(fabs(momLoss)<1.E-11)momLoss=1.E-11;
71 
72 
73  if (doNoise) {
74  if (!noise)
75  throw GFException("GFEnergyLossBetheBloch::energyLoss(): no noise matrix specified!", __LINE__, __FILE__).setFatal();
76 
77  // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
78  double sigma2E = 0.;
79  double zeta = 153.4E3 * charge*charge/(beta*beta) * matZ/matA * matDensity * step; // eV
80  double Emax = 2.E9*me*beta*beta*gammaSquare / (1. + 2.*gamma*me/mass + (me/mass)*(me/mass) ); // eV
81  double kappa = zeta/Emax;
82 
83  if (kappa > 0.01) { // Vavilov-Gaussian regime
84  sigma2E += zeta*Emax*(1.-beta*beta/2.); // eV^2
85  }
86  else { // Urban/Landau approximation
87  double alpha = 0.996;
88  double sigmaalpha = 15.76;
89  // calculate number of collisions Nc
90  double I = 16. * pow(matZ, 0.9); // eV
91  double f2 = 0.;
92  if (matZ > 2.) f2 = 2./matZ;
93  double f1 = 1. - f2;
94  double e2 = 10.*matZ*matZ; // eV
95  double e1 = pow( (I/pow(e2,f2)), 1./f1); // eV
96 
97  double mbbgg2 = 2.E9*mass*beta*beta*gammaSquare; // eV
98  double Sigma1 = dedx*1.0E9 * f1/e1 * (log(mbbgg2 / e1) - beta*beta) / (log(mbbgg2 / I) - beta*beta) * 0.6; // 1/cm
99  double Sigma2 = dedx*1.0E9 * f2/e2 * (log(mbbgg2 / e2) - beta*beta) / (log(mbbgg2 / I) - beta*beta) * 0.6; // 1/cm
100  double Sigma3 = dedx*1.0E9 * Emax / ( I*(Emax+I)*log((Emax+I)/I) ) * 0.4; // 1/cm
101 
102  double Nc = (Sigma1 + Sigma2 + Sigma3)*step;
103 
104  if (Nc > 50.) { // truncated Landau distribution
105  // calculate sigmaalpha (see GEANT3 manual W5013)
106  double RLAMED = -0.422784 - beta*beta - log(zeta/Emax);
107  double RLAMAX = 0.60715 + 1.1934*RLAMED +(0.67794 + 0.052382*RLAMED)*exp(0.94753+0.74442*RLAMED);
108  // from lambda max to sigmaalpha=sigma (empirical polynomial)
109  if(RLAMAX <= 1010.) {
110  sigmaalpha = 1.975560
111  +9.898841e-02 *RLAMAX
112  -2.828670e-04 *RLAMAX*RLAMAX
113  +5.345406e-07 *pow(RLAMAX,3.)
114  -4.942035e-10 *pow(RLAMAX,4.)
115  +1.729807e-13 *pow(RLAMAX,5.);
116  }
117  else { sigmaalpha = 1.871887E+01 + 1.296254E-02 *RLAMAX; }
118  // alpha=54.6 corresponds to a 0.9996 maximum cut
119  if(sigmaalpha > 54.6) sigmaalpha=54.6;
120  sigma2E += sigmaalpha*sigmaalpha * zeta*zeta; // eV^2
121  }
122  else { // Urban model
123  double Ealpha = I / (1.-(alpha*Emax/(Emax+I))); // eV
124  double meanE32 = I*(Emax+I)/Emax * (Ealpha-I); // eV^2
125  sigma2E += step * (Sigma1*e1*e1 + Sigma2*e2*e2 + Sigma3*meanE32); // eV^2
126  }
127  }
128 
129  sigma2E*=1.E-18; // eV -> GeV
130 
131  // update noise matrix
132  (*noise)[6][6] += (mom*mom+mass*mass)/pow(mom,6.)*sigma2E;
133  }
134 
135  return momLoss;
136 }
E
Float_t E
Definition: plot.C:23
larg4::f2
SynchrotronAndGN f2
Definition: CustomPhysicsBuiltIns.cxx:11
beta
Double_t beta
Definition: TG43_relative_dose.C:15
noise
unsigned int noise()
Definition: chem4.cc:265
genf::GFAbsEnergyLoss::getParticleMass
double getParticleMass(const int &pdg)
Returns particle mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:35
larg4::f1
EmPhysicsFactory f1
Definition: CustomPhysicsBuiltIns.cxx:10
genf::GFAbsEnergyLoss::getParticleParameters
void getParticleParameters(const int &pdg, double &charge, double &mass)
Gets particle charge and mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:27
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
GFException::setFatal
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:80
double genf::GFAbsEnergyLoss::getParticleMass ( const int &  pdg)
protectedinherited

Returns particle mass (in GeV)

Definition at line 35 of file GFAbsEnergyLoss.cxx.

References part.

Referenced by genf::GFEnergyLossBrems::energyLoss(), and energyLoss().

35  {
36  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(pdg);
37  return part->Mass();
38 }
part
TString part[npart]
Definition: Style.C:32
void genf::GFAbsEnergyLoss::getParticleParameters ( const int &  pdg,
double &  charge,
double &  mass 
)
protectedinherited

Gets particle charge and mass (in GeV)

Definition at line 27 of file GFAbsEnergyLoss.cxx.

References part.

Referenced by genf::GFEnergyLossBrems::energyLoss(), energyLoss(), and genf::GFEnergyLossCoulomb::energyLoss().

29  {
30  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(pdg);
31  charge = part->Charge()/(3.);
32  mass = part->Mass();
33 }
part
TString part[npart]
Definition: Style.C:32
The documentation for this class was generated from the following files: