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