GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BezrukovBugaevModel.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 #include <Math/Integrator.h>
13 
15 //#include "Numerical/IntegratorI.h"
17 
18 using namespace genie;
19 using namespace genie::mueloss;
20 using namespace genie::constants;
21 
22 //____________________________________________________________________________
23 BezrukovBugaevModel::BezrukovBugaevModel() :
24 MuELossI("genie::mueloss::BezrukovBugaevModel")
25 {
26 
27 }
28 //____________________________________________________________________________
30 MuELossI("genie::mueloss::BezrukovBugaevModel", config)
31 {
32 
33 }
34 //____________________________________________________________________________
36 {
37 
38 }
39 //____________________________________________________________________________
40 double BezrukovBugaevModel::dE_dx(double E, MuELMaterial_t material) const
41 {
42 // Calculate the muon -dE/dx due to muon nuclear interaction (in GeV^-2).
43 // To convert the result to more handly units, eg MeV/(gr/cm^2), just write:
44 // dE_dx /= (units::MeV/(units::g/units::cm2));
45 
46  if(material == eMuUndefined) return 0;
47  if(E<=MuELProcess::Threshold(this->Process()) || E>=kMaxMuE) return 0;
48 
49  // material Z,E
50  double Z = MuELMaterial::Z(material);
51  double A = MuELMaterial::A(material);
52 
53  // calculate (the min,max) fraction of energy, v, carried to the photon
54  double Vmin = 0.;
55  double Vmax = 1. - 0.75*kSqrtNapierConst* (kMuonMass/E) * TMath::Power(Z,1/3.);
56 
57  // integrate the Bezrukov-Bugaev differential cross section v*ds/dv for
58  // muon nuclear interaction over v
59 
60  ROOT::Math::IBaseFunctionOneDim * integrand =
62  ROOT::Math::IntegrationOneDim::Type ig_type =
64 
65  double abstol = 1; // We mostly care about relative tolerance
66  double reltol = 1E-4;
67  int nmaxeval = 100000;
68  ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval);
69 
70  // calculate the b factor (bE = -dE/dx) in GeV^-3
71  A *= units::g;
72  double bnucl = (kNA/A) * ig.Integral(Vmin, Vmax);
73 
74  delete integrand;
75 
76  // calculate the dE/dx due to muon nuclear interaction in GeV^-2
77  double de_dx = bnucl*E;
78  return de_dx;
79 }
80 //____________________________________________________________________________
81 //void BezrukovBugaevModel::Configure(const Registry & config)
82 //{
83 // Algorithm::Configure(config);
84 // this->LoadConfig();
85 //}
86 ////____________________________________________________________________________
87 //void BezrukovBugaevModel::Configure(string config)
88 //{
89 // Algorithm::Configure(config);
90 // this->LoadConfig();
91 //}
92 ////____________________________________________________________________________
93 //void BezrukovBugaevModel::LoadConfig(void)
94 //{
95 // //fIntegrator =
96 //// dynamic_cast<const IntegratorI *> (this->SubAlg("Integrator"));
97 //// assert(fIntegrator);
98 //}
99 //____________________________________________________________________________
100 //
101 // BezrukovBugaevIntegrand
102 //
103 //____________________________________________________________________________
105 ROOT::Math::IBaseFunctionOneDim()
106 {
107  fE = E;
108  fA = A;
109 }
110 //____________________________________________________________________________
112 {
113 
114 }
115 //____________________________________________________________________________
116 unsigned int gsl::BezrukovBugaevIntegrand::NDim(void) const
117 {
118  return 1;
119 }
120 //____________________________________________________________________________
121 double gsl::BezrukovBugaevIntegrand::DoEval(double xin) const
122 {
123 // Returns v*(ds/dv)
124 
125  double v = xin; // v, the fraction of energy transfered to the photon
126 
127  if (! (v >0)) return 0;
128  if ( v >1) return 0;
129  if (! (fE>0)) return 0;
130 
131  double a = kAem;
132  double pi = kPi;
133  double mmu2 = kMuonMass2;
134  double v2 = TMath::Power(v,2.);
135  double t = mmu2 *v2/(1-v);
136  double k = 1. - 2./v + 2./v2;
137  double A13 = TMath::Power(fA,1./3.);
138  double M1_2 = 0.54; // m1^2 in photonuclear diff. xsec formula (in GeV^2)
139  double M2_2 = 1.80; // m2^2 in photonuclear diff. xsec formula (in GeV^2)
140  double M1_2_t = M1_2 / t;
141  double M2_2_t = M2_2 / t;
142  double mmu2_t = mmu2 / t;
143  double d = M1_2 / (t + M1_2);
144 
145  // Calculate the cross section (in ub) for photonuclear interaction
146  double Ep = v*fE; // photon energy (GeV)
147  double loge = TMath::Log(0.0213*Ep); // factor 0.0213 has units of GeV^-1
148  double sig = 114.3 + 1.647 * loge*loge; // in ub
149 
150  // Calculate the factor G (dimensionless) appearing in the differential
151  // photonuclear interaction cross section
152  double x = 0.00282*A13*sig; // factor 0.00282 has units of ub^-1
153  double x2 = x*x;
154  double x3 = x2*x;
155  double G = 3*(0.5*x2 - 1. + (1.+x)*TMath::Exp(-x)) /x3;
156 
157  // Calculate the differential cross section ds/dv for muon nuclear
158  // interaction based on the Bezrukov-Bugaev formula.
159  double bbA = 0.5*(a/pi) * fA * sig * v;
160  double bbB = 0.75*G * ( k*TMath::Log(1.+M1_2_t) - k*d - 2.*mmu2_t );
161  double bbC = 0.25 * ( k*TMath::Log(1.+M2_2_t) - 2.*mmu2_t );
162  double bbD = 0.5*mmu2_t * ( 0.75*G*d + 0.25*M2_2_t*TMath::Log(1.+1./M2_2_t) );
163 
164  double ds_dv = bbA*(bbB+bbC+bbD); // in um (microbarns)
165  double vds_dv = v*ds_dv;
166 
167  vds_dv *= units::ub; // ub -> GeV^-2
168  return vds_dv;
169 }
170 //____________________________________________________________________________
171 ROOT::Math::IBaseFunctionOneDim *
173 {
174  return new gsl::BezrukovBugaevIntegrand(fE, fA);
175 }
176 //____________________________________________________________________________
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
static constexpr double g
Definition: Units.h:144
static const double kSqrtNapierConst
static double Z(MuELMaterial_t material)
Definition: MuELMaterial.h:236
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
static constexpr double ub
Definition: Units.h:80
static double Threshold(MuELProcess_t p)
Definition: MuELProcess.h:58
MuELProcess_t Process(void) const
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