GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
KokoulinPetrukhinModel.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/IntegratorMultiDim.h>
13 
16 
17 using namespace genie;
18 using namespace genie::mueloss;
19 using namespace genie::constants;
20 
21 //____________________________________________________________________________
22 KokoulinPetrukhinModel::KokoulinPetrukhinModel() :
23 MuELossI("genie::mueloss::KokoulinPetrukhinModel")
24 {
25 
26 }
27 //____________________________________________________________________________
29 MuELossI("genie::mueloss::KokoulinPetrukhinModel", config)
30 {
31 
32 }
33 //____________________________________________________________________________
35 {
36 
37 }
38 //____________________________________________________________________________
39 double KokoulinPetrukhinModel::dE_dx(double E, MuELMaterial_t material) const
40 {
41 // Calculate the muon -dE/dx due to e+e- pair production (in GeV^-2).
42 // To convert the result to more handly units, eg MeV/(gr/cm^2), just write:
43 // dE_dx /= (units::MeV/(units::g/units::cm2));
44 
45  if(material == eMuUndefined) return 0;
46  if(E<=MuELProcess::Threshold(this->Process()) || E>=kMaxMuE) return 0;
47 
48  // material Z,E
49  double Z = MuELMaterial::Z(material);
50  double A = MuELMaterial::A(material);
51 
52  // calculate (the min,max) fraction of energy, v, carried to the photon
53  double Vmin = 4.*kElectronMass/E;
54  double Vmax = 1. - 0.75*kSqrtNapierConst* (kMuonMass/E) * TMath::Power(Z,1/3.);
55 
56  // claculate the limits of the asymmetry parameter of the e+e- pair
57  // p = (E(+) - E(-)) / (E(+) + E(-))
58  double Pmin = 0.;
59  double Pmax = 1.;
60 
61  // integrate the Kokulin-Petrukhin differential cross section v*ds/dv for
62  // muon e+e- pair production over v and p
63 
64  ROOT::Math::IBaseFunctionMultiDim * integrand =
66  ROOT::Math::IntegrationMultiDim::Type ig_type =
68 
69  double abstol = 1; // We mostly care about relative tolerance.
70  double reltol = 1E-4;
71  int nmaxeval = 500000;
72 
73  ROOT::Math::IntegratorMultiDim ig(*integrand, ig_type, abstol, reltol, nmaxeval);
74 
75  double min[2] = { Vmin, Pmin };
76  double max[2] = { Vmax, Pmax };
77 
78  // calculate the b factor (bE = -dE/dx) in GeV^-3
79  A *= units::g;
80  double bpair = (2*kNA/A) * ig.Integral(min, max);
81 
82  delete integrand;
83 
84  // calculate the dE/dx due to muon bremsstrahlung in GeV^-2
85  double de_dx = bpair*E;
86  return de_dx;
87 }
88 //____________________________________________________________________________
89 //
90 // KokoulinPetrukhinIntegrand
91 //
92 //____________________________________________________________________________
94 ROOT::Math::IBaseFunctionMultiDim()
95 {
96  fE = E;
97  fZ = Z;
98 }
99 //____________________________________________________________________________
101 {
102 
103 }
104 //____________________________________________________________________________
106 {
107  return 2;
108 }
109 //____________________________________________________________________________
110 double gsl::KokoulinPetrukhinIntegrand::DoEval (const double * xin) const
111 {
112 // Returns v*(d^2s/dvdp)
113 
114  double v = xin[0]; // v, the fraction of energy transfered to the photon
115  double p = xin[1]; //
116 
117  if (! (v >0)) return 0;
118  if ( v >1) return 0;
119  if (! (fE>0)) return 0;
120 
121  double pmax_v = (1. - 6.*kMuonMass2 / (fE*fE*(1.-v)) ) *
122  TMath::Sqrt(1.-4.*kElectronMass/(fE*v));
123  if(p>pmax_v) return 0;
124 
125  double v2 = TMath::Power(v,2.);
126  double p2 = TMath::Power(p,2.);
127  double R = 189.;
128  double a4 = TMath::Power(kAem,4.);
129  double pi = kPi;
130  double ZLe2 = TMath::Power(fZ*kLe,2.);
131  double Zm13 = TMath::Power(fZ,-1./3.);
132  double Zm23 = TMath::Power(fZ,-2./3.);
133  double Z13 = TMath::Power(fZ,1./3.);
134  double me = kElectronMass;
135  double me2 = kElectronMass2;
136  double mmu = kMuonMass;
137  double mmu2 = kMuonMass2;
138  double memu2 = me2/mmu2;
139  double memu = me/mmu;
140  double mume = mmu/me;
141 
142  double b = 0.5*v2/(1.-v);
143  double xi = (1.-p2) * TMath::Power(0.5*v*mume, 2.) / (1.-v);
144 
145  // Approximate the FIm factor (dimensionless) for the Kokoulin Petrukhin
146  // pair production cross section
147 
148  double FImA = ( (1.+p2)*(1.+3.*b/2.) - (1.+2.*b)*(1.-p2)/xi ) * TMath::Log(1.+xi);
149  double FImB = xi*(1.-p2-b) / (1.+xi);
150  double FImC = (1.+2.*b) * (1.-p2);
151  double Ym = (4. + p2 + 3.*b*(1.+p2)) /
152  ((1.+p2)*(1.5+2.*b)*TMath::Log(3.+xi) + 1. - 1.5*p2);
153  double LmA = (2./3.) * mume * R * Zm23;
154  double LmB = 1. + (2.*me*R * kSqrtNapierConst * Zm13 * (1+xi) * (1+Ym)) / (fE*v*(1-p2) );
155  double Lm = TMath::Log(LmA/LmB);
156  double FIm = (FImA+FImB+FImC)*Lm;
157 
158  // Approximate the FIe factor (dimensionless) for the Kokoulin Petrukhin
159  // pair production cross section
160 
161  double FIeA = ( (2.+p2) * (1.+b) + xi*(3.+p2) ) * log(1.+1./xi);
162  double FIeB = (1.-p2-b) / (1.+xi);
163  double FIeC = 3. + p2;
164  double Ye = (5. - p2 + 4*b*(1+p2)) /
165  (2.*(1.+3.*b)*TMath::Log(3.+1./xi) - p2 - 2.*b*(2.-p2));
166  double x_Y = (1+xi)*(1+Ye);
167  double LeA = R*Zm13*TMath::Sqrt(x_Y);
168  double LeB = 1. + (2.*me*R*kSqrtNapierConst*Zm13*x_Y) / (fE*v*(1-p2));
169  double LeC = 1.5 * memu * Z13;
170  double LeD = 1 + TMath::Power(LeC,2.)*x_Y;
171  double Le = TMath::Log(LeA/LeB) - 0.5*TMath::Log(LeD);
172  double FIe = (FIeA+FIeB-FIeC)*Le;
173 
174  double d2s_dvdp = a4 * (2./(3.*pi)) * ZLe2 * ((1.-v)/v) * (FIe + FIm*memu2);
175 
176  double vd2s_dvdp = v*d2s_dvdp;
177  return vd2s_dvdp;
178 }
179 //____________________________________________________________________________
180 ROOT::Math::IBaseFunctionMultiDim *
182 {
183  return new gsl::KokoulinPetrukhinIntegrand(fE,fZ);
184 }
185 //____________________________________________________________________________
static constexpr double g
Definition: Units.h:144
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
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 const double kElectronMass2
double dE_dx(double E, MuELMaterial_t material) const
Implement the MuELossI interface.
static constexpr double A
Definition: Units.h:74
const double kMaxMuE
Definition: MuELossI.h:29
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
static double A(MuELMaterial_t material)
Definition: MuELMaterial.h:304
enum genie::mueloss::EMuELMaterial MuELMaterial_t