LArSoft  v09_90_00
Liquid Argon Software toolkit - https://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 {}

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 26 of file GFEnergyLossBetheBloch.cxx.

References E, f1, 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)
52  gammaSquare = 1. / gammaSquare;
53  else
54  gammaSquare = 1.E10;
55  double gamma = sqrt(gammaSquare);
56 
57  double massRatio = me / mass;
58  double argument = gammaSquare * beta * beta * me * 1.E3 * 2. /
59  ((1.E-6 * meanExcitationEnergy) *
60  sqrt(1 + 2 * sqrt(gammaSquare) * massRatio + massRatio * massRatio));
61  if (argument <= exp(beta * beta))
62  dedx = 0.;
63  else {
64  dedx *= (log(argument) - beta * beta); // Bethe-Bloch [MeV/cm]
65  dedx *= 1.E-3; // in GeV/cm, hence 1.e-3
66  if (dedx <= 0.)
67  throw GFException(
68  "GFEnergyLossBetheBloch::energyLoss(): non-positive dE/dx", __LINE__, __FILE__)
69  .setFatal();
70  }
71 
72  double DE = step * dedx; //always positive
73  double momLoss =
74  sqrt(mom * mom + 2. * sqrt(mom * mom + mass * mass) * DE + DE * DE) - mom; //always positive
75 
76  //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
77  if (fabs(momLoss) < 1.E-11) momLoss = 1.E-11;
78 
79  if (doNoise) {
80  if (!noise)
81  throw GFException(
82  "GFEnergyLossBetheBloch::energyLoss(): no noise matrix specified!", __LINE__, __FILE__)
83  .setFatal();
84 
85  // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
86  double sigma2E = 0.;
87  double zeta = 153.4E3 * charge * charge / (beta * beta) * matZ / matA * matDensity * step; // eV
88  double Emax = 2.E9 * me * beta * beta * gammaSquare /
89  (1. + 2. * gamma * me / mass + (me / mass) * (me / mass)); // eV
90  double kappa = zeta / Emax;
91 
92  if (kappa > 0.01) { // Vavilov-Gaussian regime
93  sigma2E += zeta * Emax * (1. - beta * beta / 2.); // eV^2
94  }
95  else { // Urban/Landau approximation
96  double alpha = 0.996;
97  double sigmaalpha = 15.76;
98  // calculate number of collisions Nc
99  double I = 16. * pow(matZ, 0.9); // eV
100  double f2 = 0.;
101  if (matZ > 2.) f2 = 2. / matZ;
102  double f1 = 1. - f2;
103  double e2 = 10. * matZ * matZ; // eV
104  double e1 = pow((I / pow(e2, f2)), 1. / f1); // eV
105 
106  double mbbgg2 = 2.E9 * mass * beta * beta * gammaSquare; // eV
107  double Sigma1 = dedx * 1.0E9 * f1 / e1 * (log(mbbgg2 / e1) - beta * beta) /
108  (log(mbbgg2 / I) - beta * beta) * 0.6; // 1/cm
109  double Sigma2 = dedx * 1.0E9 * f2 / e2 * (log(mbbgg2 / e2) - beta * beta) /
110  (log(mbbgg2 / I) - beta * beta) * 0.6; // 1/cm
111  double Sigma3 = dedx * 1.0E9 * Emax / (I * (Emax + I) * log((Emax + I) / I)) * 0.4; // 1/cm
112 
113  double Nc = (Sigma1 + Sigma2 + Sigma3) * step;
114 
115  if (Nc > 50.) { // truncated Landau distribution
116  // calculate sigmaalpha (see GEANT3 manual W5013)
117  double RLAMED = -0.422784 - beta * beta - log(zeta / Emax);
118  double RLAMAX = 0.60715 + 1.1934 * RLAMED +
119  (0.67794 + 0.052382 * RLAMED) * exp(0.94753 + 0.74442 * RLAMED);
120  // from lambda max to sigmaalpha=sigma (empirical polynomial)
121  if (RLAMAX <= 1010.) {
122  sigmaalpha = 1.975560 + 9.898841e-02 * RLAMAX - 2.828670e-04 * RLAMAX * RLAMAX +
123  5.345406e-07 * pow(RLAMAX, 3.) - 4.942035e-10 * pow(RLAMAX, 4.) +
124  1.729807e-13 * pow(RLAMAX, 5.);
125  }
126  else {
127  sigmaalpha = 1.871887E+01 + 1.296254E-02 * RLAMAX;
128  }
129  // alpha=54.6 corresponds to a 0.9996 maximum cut
130  if (sigmaalpha > 54.6) sigmaalpha = 54.6;
131  sigma2E += sigmaalpha * sigmaalpha * zeta * zeta; // eV^2
132  }
133  else { // Urban model
134  double Ealpha = I / (1. - (alpha * Emax / (Emax + I))); // eV
135  double meanE32 = I * (Emax + I) / Emax * (Ealpha - I); // eV^2
136  sigma2E += step * (Sigma1 * e1 * e1 + Sigma2 * e2 * e2 + Sigma3 * meanE32); // eV^2
137  }
138  }
139 
140  sigma2E *= 1.E-18; // eV -> GeV
141 
142  // update noise matrix
143  (*noise)[6][6] += (mom * mom + mass * mass) / pow(mom, 6.) * sigma2E;
144  }
145 
146  return momLoss;
147 }
f2
Float_t f2
Definition: comparison_ascii.C:56
noise
unsigned int noise()
Definition: chem4.cc:261
genf::GFAbsEnergyLoss::getParticleMass
double getParticleMass(const int &pdg)
Returns particle mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:33
E
Float_t E
Definition: plot.C:20
f1
Float_t f1
Definition: comparison_ascii.C:56
genf::GFAbsEnergyLoss::getParticleParameters
void getParticleParameters(const int &pdg, double &charge, double &mass)
Gets particle charge and mass (in GeV)
Definition: GFAbsEnergyLoss.cxx:26
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
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:75
double genf::GFAbsEnergyLoss::getParticleMass ( const int &  pdg)
protectedinherited

Returns particle mass (in GeV)

Definition at line 33 of file GFAbsEnergyLoss.cxx.

References part.

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

34 {
35  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg);
36  return part->Mass();
37 }
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 26 of file GFAbsEnergyLoss.cxx.

References part.

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

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