BetheBlochModel.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab - December 10, 2003
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14 */
15 //____________________________________________________________________________
16 
17 #include <TMath.h>
18 
22 
23 using namespace genie;
24 using namespace genie::mueloss;
25 using namespace genie::constants;
26 
27 //____________________________________________________________________________
29 MuELossI("genie::mueloss::BetheBlochModel")
30 {
31 
32 }
33 //____________________________________________________________________________
35 MuELossI("genie::mueloss::BetheBlochModel", config)
36 {
37 
38 }
39 //____________________________________________________________________________
41 {
42 
43 }
44 //____________________________________________________________________________
45 double BetheBlochModel::dE_dx(double E, MuELMaterial_t mt) const
46 {
47 // Calculates ionization dE/dx for muons via Bethe-Bloch formula (in GeV^-2)
48 // To convert the result to more handly units, eg MeV/(gr/cm^2), just write:
49 // dE_dx /= (units::MeV/(units::g/units::cm2));
50 
51  if(mt == eMuUndefined) return 0;
52  if(E<=MuELProcess::Threshold(this->Process()) || E>=kMaxMuE) return 0;
53 
54  double Z = MuELMaterial::Z(mt);
55  double A = MuELMaterial::A(mt);
56  double Z_A = Z/A; // in mol/gr
57  double a2 = kAem2; // (em coupling const)^2
58  double Na = kNA; // Avogadro's number
59  double lamda2 = kLe2/units::cm2; // (e compton wavelength)^2 in cm^2
60  double me = kElectronMass; // in GeV
61  double me2 = kElectronMass2;
62  double mmu = kMuonMass; // in GeV
63  double mmu2 = kMuonMass2;
64  double E2 = TMath::Power(E,2);
65  double beta = TMath::Sqrt(E2-mmu2)/E;
66  double beta2 = TMath::Power(beta,2);
67  double gamma = E/mmu;
68  double gamma2 = TMath::Power(gamma,2);
70  double I2 = TMath::Power(I,2); // in GeV^2
71 
72  // Calculate the maximum energy transfer to the electron (in GeV)
73 
74  double p2 = E2-mmu2;
75  double Emaxt = 2*me*p2 / (me2 + mmu2 + 2*me*E);
76  double Emaxt2 = TMath::Power(Emaxt,2);
77 
78  // Calculate the density correction factor delta
79 
85  double X = TMath::Log10(beta*gamma);
86 
87  double delta = 0;
88  if(X0<X && X<X1) delta = 4.6052*X + a*TMath::Power(X1-X,m) + C;
89  if(X>X1) delta = 4.6052*X + C;
90 
91  LOG("MuELoss", pDEBUG) << "density correction factor = " << delta;
92  LOG("MuELoss", pDEBUG) << "max energy transfer (GeV) = " << Emaxt;
93  LOG("MuELoss", pDEBUG) << "ionization potential (GeV)= " << I;
94  LOG("MuELoss", pDEBUG) << "E = " << E << ", p2 = " << p2;
95  LOG("MuELoss", pDEBUG) << "beta = " << beta << ", gamma = " << gamma;
96 
97  // Calculate the -dE/dx
98  double de_dx = a2 * (2*kPi*Na*lamda2) * Z_A * (me/beta2) *
99  (TMath::Log( 2*me*beta2*gamma2*Emaxt/I2 ) -
100  2*beta2 + 0.25*(Emaxt2/E2) - delta);
101 
102  de_dx *= (units::GeV/(units::g/units::cm2));
103  return de_dx; // in GeV^-2
104 }
105 //____________________________________________________________________________
106 
107 
const double kPi
Basic constants.
TH1F * a2
Definition: f2_nu.C:545
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static const double kLe2
Definition: Constants.h:52
double delta
Definition: runWimpSim.h:98
const int Na
static double Z(MuELMaterial_t material)
Definition: MuELMaterial.h:236
static double DensityCorrection_m(MuELMaterial_t material)
The MuELoss utility package that computes muon energy losses in the energy range from 1 GeV to 10 TeV...
static double DensityCorrection_C(MuELMaterial_t material)
static const double eV
Definition: Units.h:128
Definition: config.py:1
static const double g
Definition: Units.h:145
MuELProcess_t Process(void) const
Double_t beta
static const double cm2
Definition: Units.h:77
static const double kElectronMass
Definition: Constants.h:71
static double Threshold(MuELProcess_t p)
Definition: MuELProcess.h:58
Float_t Z
Definition: plot.C:38
const double C
static double DensityCorrection_a(MuELMaterial_t material)
Definition: NueSkimmer.h:24
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static const double kElectronMass2
Definition: Constants.h:84
static double IonizationPotential(MuELMaterial_t material)
Float_t E
Definition: plot.C:20
Double_t X1
Definition: plot.C:264
static const double kAem2
Definition: Constants.h:58
const double a
static const double kMuonMass2
Definition: Constants.h:85
static const double kNA
Definition: Constants.h:50
const double kMaxMuE
Definition: MuELossI.h:29
static double DensityCorrection_X0(MuELMaterial_t material)
static const double A
Definition: Units.h:82
static const double kMuonMass
Definition: Constants.h:72
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double dE_dx(double E, MuELMaterial_t material) const
implement the MuELossI interface
static const double GeV
Definition: Units.h:29
Float_t X
Definition: plot.C:38
static double DensityCorrection_X1(MuELMaterial_t material)
static double A(MuELMaterial_t material)
Definition: MuELMaterial.h:304
#define pDEBUG
Definition: Messenger.h:64
enum genie::mueloss::EMuELMaterial MuELMaterial_t