KokoulinPetrukhinModel.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
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/IntegratorMultiDim.h>
19 
22 
23 using namespace genie;
24 using namespace genie::mueloss;
25 using namespace genie::constants;
26 
27 //____________________________________________________________________________
29 MuELossI("genie::mueloss::KokoulinPetrukhinModel")
30 {
31 
32 }
33 //____________________________________________________________________________
35 MuELossI("genie::mueloss::KokoulinPetrukhinModel", config)
36 {
37 
38 }
39 //____________________________________________________________________________
41 {
42 
43 }
44 //____________________________________________________________________________
46 {
47 // Calculate the muon -dE/dx due to e+e- pair production (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 = 4.*kElectronMass/E;
60  double Vmax = 1. - 0.75*kSqrtNapierConst* (kMuonMass/E) * TMath::Power(Z,1/3.);
61 
62  // claculate the limits of the asymmetry parameter of the e+e- pair
63  // p = (E(+) - E(-)) / (E(+) + E(-))
64  double Pmin = 0.;
65  double Pmax = 1.;
66 
67  // integrate the Kokulin-Petrukhin differential cross section v*ds/dv for
68  // muon e+e- pair production over v and p
69 
70  ROOT::Math::IBaseFunctionMultiDim * integrand =
74 
75  double abstol = 1; // We mostly care about relative tolerance.
76  double reltol = 1E-4;
77  int nmaxeval = 500000;
78 
79  ROOT::Math::IntegratorMultiDim ig(*integrand, ig_type, abstol, reltol, nmaxeval);
80 
81  double min[2] = { Vmin, Pmin };
82  double max[2] = { Vmax, Pmax };
83 
84  // calculate the b factor (bE = -dE/dx) in GeV^-3
85  A *= units::g;
86  double bpair = (2*kNA/A) * ig.Integral(min, max);
87 
88  delete integrand;
89 
90  // calculate the dE/dx due to muon bremsstrahlung in GeV^-2
91  double de_dx = bpair*E;
92  return de_dx;
93 }
94 //____________________________________________________________________________
95 //void KokoulinPetrukhinModel::Configure(const Registry & config)
96 //{
97 // Algorithm::Configure(config);
98 // this->LoadConfig();
99 //}
100 ////____________________________________________________________________________
101 //void KokoulinPetrukhinModel::Configure(string config)
102 //{
103 // Algorithm::Configure(config);
104 // this->LoadConfig();
105 //}
106 ////____________________________________________________________________________
107 //void KokoulinPetrukhinModel::LoadConfig(void)
108 //{
109 //// fIntegrator = dynamic_cast<const IntegratorI *> (this->SubAlg("Integrator"));
110 //// assert(fIntegrator);
111 //}
112 //____________________________________________________________________________
113 //
114 // KokoulinPetrukhinIntegrand
115 //
116 //____________________________________________________________________________
118 ROOT::Math::IBaseFunctionMultiDim()
119 {
120  fE = E;
121  fZ = Z;
122 }
123 //____________________________________________________________________________
125 {
126 
127 }
128 //____________________________________________________________________________
130 {
131  return 2;
132 }
133 //____________________________________________________________________________
134 double gsl::KokoulinPetrukhinIntegrand::DoEval (const double * xin) const
135 {
136 // Returns v*(d^2s/dvdp)
137 
138  double v = xin[0]; // v, the fraction of energy transfered to the photon
139  double p = xin[1]; //
140 
141  if (! (v >0)) return 0;
142  if ( v >1) return 0;
143  if (! (fE>0)) return 0;
144 
145  double pmax_v = (1. - 6.*kMuonMass2 / (fE*fE*(1.-v)) ) *
146  TMath::Sqrt(1.-4.*kElectronMass/(fE*v));
147  if(p>pmax_v) return 0;
148 
149  double v2 = TMath::Power(v,2.);
150  double p2 = TMath::Power(p,2.);
151  double R = 189.;
152  double a4 = TMath::Power(kAem,4.);
153  double pi = kPi;
154  double ZLe2 = TMath::Power(fZ*kLe,2.);
155  double Zm13 = TMath::Power(fZ,-1./3.);
156  double Zm23 = TMath::Power(fZ,-2./3.);
157  double Z13 = TMath::Power(fZ,1./3.);
158  double me = kElectronMass;
159  double me2 = kElectronMass2;
160  double mmu = kMuonMass;
161  double mmu2 = kMuonMass2;
162  double memu2 = me2/mmu2;
163  double memu = me/mmu;
164  double mume = mmu/me;
165 
166  double b = 0.5*v2/(1.-v);
167  double xi = (1.-p2) * TMath::Power(0.5*v*mume, 2.) / (1.-v);
168 
169  // Approximate the FIm factor (dimensionless) for the Kokoulin Petrukhin
170  // pair production cross section
171 
172  double FImA = ( (1.+p2)*(1.+3.*b/2.) - (1.+2.*b)*(1.-p2)/xi ) * TMath::Log(1.+xi);
173  double FImB = xi*(1.-p2-b) / (1.+xi);
174  double FImC = (1.+2.*b) * (1.-p2);
175  double Ym = (4. + p2 + 3.*b*(1.+p2)) /
176  ((1.+p2)*(1.5+2.*b)*TMath::Log(3.+xi) + 1. - 1.5*p2);
177  double LmA = (2./3.) * mume * R * Zm23;
178  double LmB = 1. + (2.*me*R * kSqrtNapierConst * Zm13 * (1+xi) * (1+Ym)) / (fE*v*(1-p2) );
179  double Lm = TMath::Log(LmA/LmB);
180  double FIm = (FImA+FImB+FImC)*Lm;
181 
182  // Approximate the FIe factor (dimensionless) for the Kokoulin Petrukhin
183  // pair production cross section
184 
185  double FIeA = ( (2.+p2) * (1.+b) + xi*(3.+p2) ) * log(1.+1./xi);
186  double FIeB = (1.-p2-b) / (1.+xi);
187  double FIeC = 3. + p2;
188  double Ye = (5. - p2 + 4*b*(1+p2)) /
189  (2.*(1.+3.*b)*TMath::Log(3.+1./xi) - p2 - 2.*b*(2.-p2));
190  double x_Y = (1+xi)*(1+Ye);
191  double LeA = R*Zm13*TMath::Sqrt(x_Y);
192  double LeB = 1. + (2.*me*R*kSqrtNapierConst*Zm13*x_Y) / (fE*v*(1-p2));
193  double LeC = 1.5 * memu * Z13;
194  double LeD = 1 + TMath::Power(LeC,2.)*x_Y;
195  double Le = TMath::Log(LeA/LeB) - 0.5*TMath::Log(LeD);
196  double FIe = (FIeA+FIeB-FIeC)*Le;
197 
198  double d2s_dvdp = a4 * (2./(3.*pi)) * ZLe2 * ((1.-v)/v) * (FIe + FIm*memu2);
199 
200  double vd2s_dvdp = v*d2s_dvdp;
201  return vd2s_dvdp;
202 }
203 //____________________________________________________________________________
204 ROOT::Math::IBaseFunctionMultiDim *
206 {
208 }
209 //____________________________________________________________________________
210 
const double kPi
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
static const double kSqrtNapierConst
Definition: Constants.h:45
static const double kLe
Definition: Constants.h:51
static double Z(MuELMaterial_t material)
Definition: MuELMaterial.h:236
const char * p
Definition: xmltok.h:285
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 const double kElectronMass
Definition: Constants.h:71
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
Definition: NueSkimmer.h:24
static const double kElectronMass2
Definition: Constants.h:84
double dE_dx(double E, MuELMaterial_t material) const
Implement the MuELossI interface.
Float_t E
Definition: plot.C:20
static const double kMuonMass2
Definition: Constants.h:85
static const double kNA
Definition: Constants.h:50
const double kMaxMuE
Definition: MuELossI.h:29
#define R(x)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
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 integrand(double *x, double *par)
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
TH1F * a4
Definition: f2_nu.C:666
const hit & b
Definition: hits.cxx:21
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
static double A(MuELMaterial_t material)
Definition: MuELMaterial.h:304
enum genie::mueloss::EMuELMaterial MuELMaterial_t