GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PetrukhinShestakovModel.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 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>
8  University of Liverpool - 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 #include <Math/Integrator.h>
19 
22 
23 using namespace genie;
24 using namespace genie::mueloss;
25 using namespace genie::constants;
26 
27 //____________________________________________________________________________
28 PetrukhinShestakovModel::PetrukhinShestakovModel() :
29 MuELossI("genie::mueloss::PetrukhinShestakovModel")
30 {
31 
32 }
33 //____________________________________________________________________________
35 MuELossI("genie::mueloss::PetrukhinShestakovModel", config)
36 {
37 
38 }
39 //____________________________________________________________________________
41 {
42 
43 }
44 //____________________________________________________________________________
45 double PetrukhinShestakovModel::dE_dx(double E, MuELMaterial_t material) const
46 {
47 // Calculate the muon -dE/dx due to muon bremsstrahlung (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(material == eMuUndefined) return 0;
52  if(E<=MuELProcess::Threshold(this->Process()) || E>=kMaxMuE) return 0;
53 
54  // material Z,E
55  double Z = MuELMaterial::Z(material);
56  double A = MuELMaterial::A(material);
57 
58  // calculate (the min,max) fraction of energy, v, carried to the photon
59  double Vmin = 0.;
60  double Vmax = 1. - 0.75*kSqrtNapierConst* (kMuonMass/E) * TMath::Power(Z,1/3.);
61 
62  // integrate the Bethe-Heitler muon bremsstrahlung
63  // cross section v*ds/dv over v
64 
65  ROOT::Math::IBaseFunctionOneDim * integrand =
67  ROOT::Math::IntegrationOneDim::Type ig_type =
69 
70  double abstol = 1; // We mostly care about relative tolerance
71  double reltol = 1E-4;
72  int nmaxeval = 100000;
73  ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval);
74 
75  // calculate the b factor (bE = -dE/dx) in GeV^-3
76  A *= units::g;
77  double bbrem = (kNA/A) * ig.Integral(Vmin, Vmax);
78 
79  delete integrand;
80 
81  // calculate the dE/dx due to muon bremsstrahlung in GeV^-2
82  double de_dx = bbrem*E;
83  return de_dx;
84 }
85 //____________________________________________________________________________
87 ROOT::Math::IBaseFunctionOneDim()
88 {
89  fE = E;
90  fZ = Z;
91 }
92 //____________________________________________________________________________
94 {
95 
96 }
97 //____________________________________________________________________________
99 {
100  return 1;
101 }
102 //____________________________________________________________________________
104 {
105 // Calculate the Bethe-Heitler cross section ds/dv for muon bremsstrahlung
106 // Returns v*(ds/dv)
107 
108  double v = xin; // v, the fraction of energy transfered to the photon
109  double v2 = TMath::Power(v,2.);
110 
111  if (! (v >0)) return 0;
112  if ( v >1) return 0;
113  if (! (fE>0)) return 0;
114 
115  // Some constants...
116  double Z2 = TMath::Power(fZ,2.);
117  double Zm13 = TMath::Power(fZ,-1./3.);
118  double Zm23 = TMath::Power(fZ,-2./3.);
119  double a3 = TMath::Power(kAem,3.); // (e/m coupling const)^3
120  double me = kElectronMass;
121  double mmu = kMuonMass;
122  double mmu2 = kMuonMass2;
123  double mmue = mmu/me;
124  double memu = me/mmu;
125  double memu2 = TMath::Power(memu,2);
126 
127  // Calculate the minimum monentum transfer to the nucleus (in GeV)
128  double delta = (mmu2/fE) * 0.5*v/(1.-v);
129 
130  // Calculate the fi(delta) factor for the bremsstrahlung cross section
131  // ds/dv according to the Petrukhin/Shestakov formula (dimensionless)
132  double a = ( (fZ<10) ? 189.*mmue * Zm13 : 189.*mmue * (2./3.)*Zm23 );
133  double b = 1. + (189./me) * Zm13 * kSqrtNapierConst * delta;
134  double fi = TMath::Log(a/b);
135 
136  // Calculate the Bethe-Heitler cross section ds/dv for muon bremsstrahlung
137  double ds_dv = (a3*memu2*kLe2) * (4*Z2) * (fi) * (4/3.-4*v/3.+v2)/v;
138  double vds_dv = v*ds_dv;
139  return vds_dv; // in GeV^-2
140 }
141 //____________________________________________________________________________
142 ROOT::Math::IBaseFunctionOneDim *
144 {
145  return new gsl::PetrukhinShestakovIntegrand(fE, fZ);
146 }
147 //____________________________________________________________________________
148 
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
static constexpr double g
Definition: Units.h:144
static const double kSqrtNapierConst
static double Z(MuELMaterial_t material)
Definition: MuELMaterial.h:236
static const double kElectronMass
static double Threshold(MuELProcess_t p)
Definition: MuELProcess.h:58
static constexpr double b
Definition: Units.h:78
static constexpr double A
Definition: Units.h:74
const double a
const double kMaxMuE
Definition: MuELossI.h:29
double dE_dx(double E, MuELMaterial_t material) const
Implement the MuELossI interface.
static double A(MuELMaterial_t material)
Definition: MuELMaterial.h:304
enum genie::mueloss::EMuELMaterial MuELMaterial_t