Public 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 ()
 
- Public Member Functions inherited from genf::GFAbsEnergyLoss
virtual ~GFAbsEnergyLoss ()
 

Additional Inherited Members

- Protected Member Functions inherited from genf::GFAbsEnergyLoss
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 46 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.

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 }
pionana::beta
double beta(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:59
reco_momentum_tuples.f2
f2
Definition: reco_momentum_tuples.py:79
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
genf::GFAbsEnergyLoss::getParticleMass
double getParticleMass(const int &pdg)
Returns particle mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:35
alpha
double alpha
Definition: doAna.cpp:15
pionana::gamma
double gamma(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:55
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:48
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:78
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
reco_momentum_tuples.f1
f1
Definition: reco_momentum_tuples.py:78
The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11