GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BetheBlochModel.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
16 
17 using namespace genie;
18 using namespace genie::mueloss;
19 using namespace genie::constants;
20 
21 //____________________________________________________________________________
22 BetheBlochModel::BetheBlochModel() :
23 MuELossI("genie::mueloss::BetheBlochModel")
24 {
25 
26 }
27 //____________________________________________________________________________
29 MuELossI("genie::mueloss::BetheBlochModel", config)
30 {
31 
32 }
33 //____________________________________________________________________________
35 {
36 
37 }
38 //____________________________________________________________________________
39 double BetheBlochModel::dE_dx(double E, MuELMaterial_t mt) const
40 {
41 // Calculates ionization dE/dx for muons via Bethe-Bloch formula (in GeV^-2)
42 // To convert the result to more handly units, eg MeV/(gr/cm^2), just write:
43 // dE_dx /= (units::MeV/(units::g/units::cm2));
44 
45  if(mt == eMuUndefined) return 0;
46  if(E<=MuELProcess::Threshold(this->Process()) || E>=kMaxMuE) return 0;
47 
48  double Z = MuELMaterial::Z(mt);
49  double A = MuELMaterial::A(mt);
50  double Z_A = Z/A; // in mol/gr
51  double a2 = kAem2; // (em coupling const)^2
52  double Na = kNA; // Avogadro's number
53  double lamda2 = kLe2/units::cm2; // (e compton wavelength)^2 in cm^2
54  double me = kElectronMass; // in GeV
55  double me2 = kElectronMass2;
56  double mmu = kMuonMass; // in GeV
57  double mmu2 = kMuonMass2;
58  double E2 = TMath::Power(E,2);
59  double beta = TMath::Sqrt(E2-mmu2)/E;
60  double beta2 = TMath::Power(beta,2);
61  double gamma = E/mmu;
62  double gamma2 = TMath::Power(gamma,2);
64  double I2 = TMath::Power(I,2); // in GeV^2
65 
66  // Calculate the maximum energy transfer to the electron (in GeV)
67 
68  double p2 = E2-mmu2;
69  double Emaxt = 2*me*p2 / (me2 + mmu2 + 2*me*E);
70  double Emaxt2 = TMath::Power(Emaxt,2);
71 
72  // Calculate the density correction factor delta
73 
79  double X = TMath::Log10(beta*gamma);
80 
81  double delta = 0;
82  if(X0<X && X<X1) delta = 4.6052*X + a*TMath::Power(X1-X,m) + C;
83  if(X>X1) delta = 4.6052*X + C;
84 
85  LOG("MuELoss", pDEBUG) << "density correction factor = " << delta;
86  LOG("MuELoss", pDEBUG) << "max energy transfer (GeV) = " << Emaxt;
87  LOG("MuELoss", pDEBUG) << "ionization potential (GeV)= " << I;
88  LOG("MuELoss", pDEBUG) << "E = " << E << ", p2 = " << p2;
89  LOG("MuELoss", pDEBUG) << "beta = " << beta << ", gamma = " << gamma;
90 
91  // Calculate the -dE/dx
92  double de_dx = a2 * (2*kPi*Na*lamda2) * Z_A * (me/beta2) *
93  (TMath::Log( 2*me*beta2*gamma2*Emaxt/I2 ) -
94  2*beta2 + 0.25*(Emaxt2/E2) - delta);
95 
96  de_dx *= (units::GeV/(units::g/units::cm2));
97  return de_dx; // in GeV^-2
98 }
99 //____________________________________________________________________________
static constexpr double g
Definition: Units.h:144
static double Z(MuELMaterial_t material)
Definition: MuELMaterial.h:236
static double DensityCorrection_m(MuELMaterial_t material)
static double DensityCorrection_C(MuELMaterial_t material)
MuELProcess_t Process(void) const
static const double kElectronMass
static double Threshold(MuELProcess_t p)
Definition: MuELProcess.h:58
static double DensityCorrection_a(MuELMaterial_t material)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static const double kElectronMass2
static constexpr double eV
Definition: Units.h:127
static double IonizationPotential(MuELMaterial_t material)
static constexpr double A
Definition: Units.h:74
static constexpr double cm2
Definition: Units.h:69
static constexpr double GeV
Definition: Units.h:28
const double a
const double kMaxMuE
Definition: MuELossI.h:29
static double DensityCorrection_X0(MuELMaterial_t material)
double dE_dx(double E, MuELMaterial_t material) const
implement the MuELossI interface
static double DensityCorrection_X1(MuELMaterial_t material)
static constexpr double m
Definition: Units.h:71
static double A(MuELMaterial_t material)
Definition: MuELMaterial.h:304
#define pDEBUG
Definition: Messenger.h:63
enum genie::mueloss::EMuELMaterial MuELMaterial_t