| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file electromagnetic/TestEm17/src/MuCross 26 /// \file electromagnetic/TestEm17/src/MuCrossSections.cc
27 /// \brief Implementation of the MuCrossSectio 27 /// \brief Implementation of the MuCrossSections class
28 // 28 //
29 // << 29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "MuCrossSections.hh" 33 #include "MuCrossSections.hh"
34 34
35 #include "G4DataVector.hh" <<
36 #include "G4Exp.hh" <<
37 #include "G4Log.hh" <<
38 #include "G4Material.hh" 35 #include "G4Material.hh"
39 #include "G4MuonMinus.hh" <<
40 #include "G4NistManager.hh" <<
41 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
43 38
44 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 40
46 using namespace std; 41 using namespace std;
47 << 42
48 MuCrossSections::MuCrossSections() 43 MuCrossSections::MuCrossSections()
49 { << 44 { }
50 fNist = G4NistManager::Instance(); <<
51 fMuonMass = G4MuonMinus::MuonMinus()->GetPDG <<
52 fMueRatio = fMuonMass / CLHEP::electron_mass <<
53 } <<
54 45
55 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 47
57 G4double MuCrossSections::CR_Macroscopic(const << 48 MuCrossSections::~MuCrossSections()
>> 49 { }
>> 50
>> 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 52
>> 53 G4double MuCrossSections::CR_Macroscopic(const G4String& process,
>> 54 const G4Material* material,
58 G4dou 55 G4double tkin, G4double ep)
>> 56
59 // return the macroscopic cross section (1/L) 57 // return the macroscopic cross section (1/L) in GEANT4 internal units
60 { 58 {
61 const G4ElementVector* theElementVector = ma 59 const G4ElementVector* theElementVector = material->GetElementVector();
62 const G4double* NbOfAtomsPerVolume = materia 60 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
63 61
64 G4double SIGMA = 0.; 62 G4double SIGMA = 0.;
65 G4int nelm = material->GetNumberOfElements() << 63 G4int nelm = material->GetNumberOfElements();
66 for (G4int i = 0; i < nelm; ++i) { << 64 for (G4int i=0; i < nelm; ++i)
>> 65 {
67 const G4Element* element = (*theElementVec 66 const G4Element* element = (*theElementVector)[i];
68 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAto 67 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAtom(process, element, tkin, ep);
69 } 68 }
70 return SIGMA; 69 return SIGMA;
71 } 70 }
72 71
73 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 73
75 G4double MuCrossSections::CR_PerAtom(const G4S << 74 G4double MuCrossSections::CR_PerAtom(const G4String& process,
>> 75 const G4Element* element,
76 G4double 76 G4double tkin, G4double ep)
77 { 77 {
78 G4double z = element->GetZ(); << 78 G4double z = element->GetZ();
79 G4double a = element->GetA(); << 79 G4double a = element->GetA();
80 G4double sigma = 0.; <<
81 if (process == "muBrems") <<
82 sigma = CRB_Mephi(z, a / (g / mole), tkin <<
83 <<
84 else if (process == "muIoni") <<
85 sigma = CRK_Mephi(z, a / (g / mole), tkin <<
86 <<
87 // else if (process == "muNucl") <<
88 else if (process == "muonNuclear") <<
89 sigma = CRN_Mephi(z, a / (g / mole), tkin <<
90 80
91 else if (process == "muPairProd") << 81 G4double sigma = 0.;
92 sigma = CRP_Mephi(z, a / (g / mole), tkin << 82 if (process == "muBrems")
93 << 83 sigma = CRB_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro;
94 else if (process == "muToMuonPairProd") << 84
95 sigma = CRM_Mephi(z, tkin, ep); << 85 else if (process == "muIoni")
96 << 86 sigma = CRK_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro;
97 return sigma; << 87
>> 88 //else if (process == "muNucl")
>> 89 else if (process == "muonNuclear")
>> 90 sigma = CRN_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro;
>> 91
>> 92 else if (process == "muPairProd")
>> 93 sigma = CRP_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro;
>> 94
>> 95 return sigma;
98 } 96 }
99 97
100 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 99
102 G4double MuCrossSections::CRB_Mephi(G4double z << 100 double MuCrossSections::CRB_Mephi(double z,double a,double tkin,double ep)
103 101
104 //******************************************** 102 //***********************************************************************
105 //*** crb_g4_1.inc in comparison 103 //*** crb_g4_1.inc in comparison with crb_.inc, following
106 //*** changes are introduced (September 104 //*** changes are introduced (September 24th, 1998):
107 //*** 1) function crb_g4 (Z,A,T 105 //*** 1) function crb_g4 (Z,A,Tkin,EP), Tkin is kinetic energy
108 //*** 2) special case of hidrog 106 //*** 2) special case of hidrogen (Z<1.5; b,b1,Dn_star)
109 //*** Numerical comparison: 5 d 107 //*** Numerical comparison: 5 decimal digits coincide.
110 //*** 108 //***
111 //*** Cross section for bremsstrahlung 109 //*** Cross section for bremsstrahlung by fast muon
112 //*** By R.P.Kokoulin, September 1998 110 //*** By R.P.Kokoulin, September 1998
113 //*** Formulae from Kelner,Kokoulin,Pet 111 //*** Formulae from Kelner,Kokoulin,Petrukhin 1995, Preprint MEPhI
114 //*** (7,18,19,20,21,25,26); Dn (18) is 112 //*** (7,18,19,20,21,25,26); Dn (18) is modified to incorporate
115 //*** Bugaev's inelatic nuclear correct 113 //*** Bugaev's inelatic nuclear correction (28) for Z > 1.
116 //******************************************** 114 //***********************************************************************
117 { 115 {
118 // G4double Z,A,Tkin,EP; << 116 // G4double Z,A,Tkin,EP;
119 G4double crb_g4; << 117 G4double crb_g4;
120 G4double e, v, delta, rab0, z_13, dn, b, b1, << 118 G4double e,v,delta,rab0,z_13,dn,b,b1,dn_star,rab1,fn,epmax1,fe,rab2;
121 // << 119 //
122 G4double ame = 0.51099907e-3; // GeV << 120 G4double ame=0.51099907e-3; // GeV
123 G4double lamu = 0.105658389; // GeV << 121 G4double lamu=0.105658389; // GeV
124 G4double re = 2.81794092e-13; // cm << 122 G4double re=2.81794092e-13; // cm
125 G4double avno = 6.022137e23; << 123 G4double avno=6.022137e23;
126 G4double alpha = 1. / 137.036; << 124 G4double alpha=1./137.036;
127 G4double rmass = lamu / ame; // "207" << 125 G4double rmass=lamu/ame; // "207"
128 G4double coeff = 16. / 3. * alpha * avno * ( << 126 G4double coeff=16./3.*alpha*avno*(re/rmass)*(re/rmass); // cm^2
129 G4double sqrte = 1.64872127; // sqrt(2.7182 << 127 G4double sqrte=1.64872127; // sqrt(2.71828...)
130 G4double btf = 183.; << 128 G4double btf=183.;
131 G4double btf1 = 1429.; << 129 G4double btf1=1429.;
132 G4double bh = 202.4; << 130 G4double bh=202.4;
133 G4double bh1 = 446.; << 131 G4double bh1=446.;
134 //*** << 132 //***
135 if (ep >= tkin) { << 133 if(ep >= tkin)
136 return 0.; << 134 {
137 } << 135 crb_g4=0.;
138 e = tkin + lamu; << 136 return crb_g4;
139 v = ep / e; << 137 }
140 delta = lamu * lamu * v / (2. * (e - ep)); << 138 e=tkin+lamu;
141 rab0 = delta * sqrte; << 139 v=ep/e;
142 G4int Z = G4lrint(z); << 140 delta=lamu*lamu*v/(2.*(e-ep)); // qmin
143 z_13 = 1. / fNist->GetZ13(z); // << 141 rab0=delta*sqrte;
144 //*** nuclear size and excitation, scr << 142 z_13=pow(z,-0.3333333); //
145 dn = 1.54 * fNist->GetA27(Z); << 143 //*** nuclear size and excitation, screening parameters
146 if (z <= 1.5) // special case for hydrogen << 144 dn=1.54*pow(a,0.27);
147 { << 145 if(z <= 1.5) // special case for hydrogen
148 b = bh; << 146 {
149 b1 = bh1; << 147 b=bh;
150 dn_star = dn; << 148 b1=bh1;
151 } << 149 dn_star=dn;
152 else { << 150 }
153 b = btf; << 151 else
154 b1 = btf1; << 152 {
155 dn_star = pow(dn, (1. - 1. / z)); // with << 153 b=btf;
156 } << 154 b1=btf1;
157 //*** nucleus contribution lo << 155 dn_star=pow(dn,(1.-1./z)); // with Bugaev's correction
158 rab1 = b * z_13; << 156 }
159 fn = G4Log(rab1 / (dn_star * (ame + rab0 * r << 157 //*** nucleus contribution logarithm
160 if (fn < 0.) fn = 0.; << 158 rab1=b*z_13;
161 //*** electron contribution l << 159 fn=log(rab1/(dn_star*(ame+rab0*rab1))*(lamu+delta*(dn_star*sqrte-2.)));
162 epmax1 = e / (1. + lamu * rmass / (2. * e)); << 160 if(fn < 0.) fn=0.;
163 if (ep >= epmax1) { << 161 //*** electron contribution logarithm
164 fe = 0.; << 162 epmax1=e/(1.+lamu*rmass/(2.*e));
165 } << 163 if(ep >= epmax1)
166 else { << 164 {
167 rab2 = b1 * z_13 * z_13; << 165 fe=0.;
168 fe = G4Log(rab2 * lamu / ((1. + delta * rm << 166 goto label10;
169 if (fe < 0.) fe = 0.; << 167 }
170 } << 168 rab2=b1*z_13*z_13;
171 crb_g4 = coeff * (1. - v * (1. - 0.75 * v)) << 169 fe=log(rab2*lamu/((1.+delta*rmass/(ame*sqrte))*(ame+rab0*rab2)));
172 return crb_g4; << 170 if(fe < 0.) fe=0.;
>> 171 //***
>> 172 label10:
>> 173 crb_g4=coeff*(1.-v*(1.-0.75*v))*z*(z*fn+fe)/(a*ep);
>> 174 return crb_g4;
173 } 175 }
174 176
175 //....oooOO0OOooo........oooOO0OOooo........oo 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 178
177 G4double MuCrossSections::CRK_Mephi(G4double z << 179 G4double MuCrossSections::CRK_Mephi(G4double z,G4double a,G4double tkin,G4double ep)
>> 180
178 //******************************************** 181 //***********************************************************************
179 //*** Cross section for knock-on electr 182 //*** Cross section for knock-on electron production by fast muons
180 //*** (including bremsstrahlung e-diagr 183 //*** (including bremsstrahlung e-diagrams and rad. correction).
181 //*** Units: cm^2/(g*GeV); Tkin, ep - G 184 //*** Units: cm^2/(g*GeV); Tkin, ep - GeV.
182 //*** By R.P.Kokoulin, October 1998 185 //*** By R.P.Kokoulin, October 1998
183 //*** Formulae from Kelner,Kokoulin,Pet 186 //*** Formulae from Kelner,Kokoulin,Petrukhin, Phys.Atom.Nuclei, 1997
184 //*** (a bit simplified Kelner's versio 187 //*** (a bit simplified Kelner's version of Eq.30 - with 2 logarithms).
185 //*** 188 //***
186 { 189 {
187 // G4double Z,A,Tkin,EP; << 190 // G4double Z,A,Tkin,EP;
188 G4double crk_g4; << 191 G4double crk_g4;
189 G4double v, sigma0, a1, a3; << 192 G4double v,sigma0,a1,a3;
190 // << 193 //
191 G4double ame = 0.51099907e-3; // GeV << 194 G4double ame=0.51099907e-3; // GeV
192 G4double lamu = 0.105658389; // GeV << 195 G4double lamu=0.105658389; // GeV
193 G4double re = 2.81794092e-13; // cm << 196 G4double re=2.81794092e-13; // cm
194 G4double avno = 6.022137e23; << 197 G4double avno=6.022137e23;
195 G4double alpha = 1. / 137.036; << 198 G4double alpha=1./137.036;
196 G4double lpi = 3.141592654; << 199 G4double lpi=3.141592654;
197 G4double bmu = lamu * lamu / (2. * ame); << 200 G4double bmu=lamu*lamu/(2.*ame);
198 G4double coeff0 = avno * 2. * lpi * ame * re << 201 G4double coeff0=avno*2.*lpi*ame*re*re;
199 G4double coeff1 = alpha / (2. * lpi); << 202 G4double coeff1=alpha/(2.*lpi);
200 //*** << 203 //***
201 << 204
202 G4double e = tkin + lamu; << 205 G4double e=tkin+lamu;
203 G4double epmax = e / (1. + bmu / e); << 206 if(e < 0.) {
204 if (ep >= epmax) { << 207 G4cout << "CRK: " << tkin << " " << ep << " " << e << G4endl;
205 return 0.0; << 208 }
206 } << 209 G4double epmax=e/(1.+bmu/e);
207 v = ep / e; << 210 if(ep >= epmax)
208 sigma0 = coeff0 * (z / a) * (1. - ep / epmax << 211 {
209 a1 = G4Log(1. + 2. * ep / ame); << 212 crk_g4=0.;
210 a3 = G4Log(4. * e * (e - ep) / (lamu * lamu) << 213 return crk_g4;
211 crk_g4 = sigma0 * (1. + coeff1 * a1 * (a3 - << 214 }
212 return crk_g4; << 215 v=ep/e;
>> 216 sigma0=coeff0*(z/a)*(1.-ep/epmax+0.5*v*v)/(ep*ep);
>> 217 a1=std::log(1.+2.*ep/ame);
>> 218 a3=std::log(4.*e*(e-ep)/(lamu*lamu));
>> 219 crk_g4=sigma0*(1.+coeff1*a1*(a3-a1));
>> 220 return crk_g4;
213 } 221 }
214 222
215 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216 224
217 G4double MuCrossSections::CRN_Mephi(G4double, << 225 G4double MuCrossSections::CRN_Mephi(G4double /* z */,G4double a,G4double tkin,G4double ep)
218 226
219 //******************************************** 227 //***********************************************************************
220 //*** Differential cross section for ph 228 //*** Differential cross section for photonuclear muon interaction.
221 //*** Formulae from Borog & Petrukhin, 229 //*** Formulae from Borog & Petrukhin, 1975
222 //*** Real photon cross section: Caldwe 230 //*** Real photon cross section: Caldwell e.a., 1979
223 //*** Nuclear shadowing: Brodsky e.a., 231 //*** Nuclear shadowing: Brodsky e.a., 1972
224 //*** Units: cm^2 / g GeV. 232 //*** Units: cm^2 / g GeV.
225 //*** CRN_G4_1.inc January 31st, 233 //*** CRN_G4_1.inc January 31st, 1998 R.P.Kokoulin
226 //******************************************** 234 //***********************************************************************
227 { 235 {
228 // G4double Z,A,Tkin,EP; << 236 // G4double Z,A,Tkin,EP;
229 G4double crn_g4; << 237 G4double crn_g4;
230 G4double e, aeff, sigph, v, v1, v2, amu2, up << 238 G4double e,aeff,sigph,v,v1,v2,amu2,up,down;
231 //*** << 239 //***
232 G4double lamu = 0.105658389; // GeV << 240 G4double lamu=0.105658389; // GeV
233 G4double avno = 6.022137e23; << 241 G4double avno=6.022137e23;
234 G4double amp = 0.9382723; // GeV << 242 G4double amp=0.9382723; // GeV
235 G4double lpi = 3.14159265; << 243 G4double lpi=3.14159265;
236 G4double alpha = 1. / 137.036; << 244 G4double alpha=1./137.036;
237 //*** << 245 //***
238 G4double epmin_phn = 0.20; // GeV << 246 G4double epmin_phn=0.20; // GeV
239 G4double alam2 = 0.400000; // GeV**2 << 247 G4double alam2=0.400000; // GeV**2
240 G4double alam = 0.632456; // sqrt(alam2) << 248 G4double alam =0.632456; // sqrt(alam2)
241 G4double coeffn = alpha / lpi * avno * 1e-30 << 249 G4double coeffn=alpha/lpi*avno*1e-30; // cm^2/microbarn
242 //*** << 250 //***
243 e = tkin + lamu; << 251 e=tkin+lamu;
244 crn_g4 = 0.; << 252 crn_g4=0.;
245 if (ep >= e - 0.5 * amp || ep <= epmin_phn) << 253 if(ep >= e-0.5*amp) return crn_g4;
246 aeff = 0.22 * a + 0.78 * pow(a, 0.89); // s << 254 if(ep <= epmin_phn) return crn_g4;
247 sigph = 49.2 + 11.1 * G4Log(ep) + 151.8 / st << 255 aeff=0.22*a+0.78*pow(a,0.89); // shadowing
248 v = ep / e; << 256 sigph=49.2+11.1*log(ep)+151.8/sqrt(ep); // microbarn
249 v1 = 1. - v; << 257 v=ep/e;
250 v2 = v * v; << 258 v1=1.-v;
251 amu2 = lamu * lamu; << 259 v2=v*v;
252 up = e * e * v1 / amu2 * (1. + amu2 * v2 / ( << 260 amu2=lamu*lamu;
253 down = 1. + ep / alam * (1. + alam / (2. * a << 261 up=e*e*v1/amu2*(1.+amu2*v2/(alam2*v1));
254 crn_g4 = coeffn * aeff / a * sigph / ep << 262 down=1.+ep/alam*(1.+alam/(2.*amp)+ep/alam);
255 * (-v1 + (v1 + 0.5 * v2 * (1. + 2. << 263 crn_g4=coeffn*aeff/a*sigph/ep*(-v1+(v1+0.5*v2*(1.+2.*amu2/alam2))*log(up/down));
256 if (crn_g4 < 0.) crn_g4 = 0.; << 264 if(crn_g4 < 0.) crn_g4=0.;
257 return crn_g4; << 265 return crn_g4;
258 } 266 }
259 267
260 //....oooOO0OOooo........oooOO0OOooo........oo 268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 269
262 G4double MuCrossSections::CRP_Mephi(G4double z << 270 G4double MuCrossSections::CRP_Mephi(G4double z,G4double a,G4double tkin,G4double ep)
263 271
264 //******************************************** 272 //**********************************************************************
265 //*** crp_g4_1.inc in comparison 273 //*** crp_g4_1.inc in comparison with crp_m.inc, following
266 //*** changes are introduced (January 1 274 //*** changes are introduced (January 16th, 1998):
267 //*** 1) Avno/A, cm^2/gram GeV 275 //*** 1) Avno/A, cm^2/gram GeV
268 //*** 2) zeta_loss(E,Z) 276 //*** 2) zeta_loss(E,Z) from Kelner 1997, Eqs.(53-54)
269 //*** 3) function crp_g4 (Z,A,T 277 //*** 3) function crp_g4 (Z,A,Tkin,EP), Tkin is kinetic energy
270 //*** 4) bbb=183 (Thomas 278 //*** 4) bbb=183 (Thomas-Fermi)
271 //*** 5) special case of hidrog 279 //*** 5) special case of hidrogen (Z<1.5), bbb,g1,g2
272 //*** 6) expansions in 'xi' are 280 //*** 6) expansions in 'xi' are simplified (Jan.17th,1998)
273 //*** 281 //***
274 //*** Cross section for electron pair p 282 //*** Cross section for electron pair production by fast muon
275 //*** By R.P.Kokoulin, December 1997 283 //*** By R.P.Kokoulin, December 1997
276 //*** Formulae from Kokoulin & Petrukhi 284 //*** Formulae from Kokoulin & Petrukhin 1971, Hobart, Eqs.(1,8,9,10)
277 { 285 {
278 // G4double Z,A,Tkin,EP; << 286 // G4double Z,A,Tkin,EP;
279 G4double crp_g4; << 287 G4double crp_g4;
280 G4double bbbtf, bbbh, g1tf, g2tf, g1h, g2h, << 288 G4double bbbtf,bbbh,g1tf,g2tf,g1h,g2h,e,z13,e1,alf,a3,bbb;
281 G4double g1, g2, zeta1, zeta2, zeta, z2, scr << 289 G4double g1,g2,zeta1,zeta2,zeta,z2,screen0,a0,a1,bet,xi0,del;
282 G4double tmn, sum, a4, a5, a6, a7, a9, xi, x << 290 G4double tmn,sum,a4,a5,a6,a7,a9,xi,xii,xi1,screen,yeu,yed,ye1;
283 G4double ale, cre, be, fe, ymu, ymd, ym1, al << 291 G4double ale,cre,be,fe,ymu,ymd,ym1,alm_crm,a10,bm,fm;
284 // << 292 //
285 G4double ame = 0.51099907e-3; // GeV << 293 G4double ame=0.51099907e-3; // GeV
286 G4double lamu = 0.105658389; // GeV << 294 G4double lamu=0.105658389; // GeV
287 G4double re = 2.81794092e-13; // cm << 295 G4double re=2.81794092e-13; // cm
288 G4double avno = 6.022137e23; << 296 G4double avno=6.022137e23;
289 G4double lpi = 3.14159265; << 297 G4double lpi=3.14159265;
290 G4double alpha = 1. / 137.036; << 298 G4double alpha=1./137.036;
291 G4double rmass = lamu / ame; // "207" << 299 G4double rmass=lamu/ame; // "207"
292 G4double coeff = 4. / (3. * lpi) * (alpha * << 300 G4double coeff=4./(3.*lpi)*(alpha*re)*(alpha*re)*avno; // cm^2
293 G4double sqrte = 1.64872127; // sqrt(2.7182 << 301 G4double sqrte=1.64872127; // sqrt(2.71828...)
294 G4double c3 = 3. * sqrte * lamu / 4.; // fo << 302 G4double c3=3.*sqrte*lamu/4.; // for limits
295 G4double c7 = 4. * ame; // -"- << 303 G4double c7=4.*ame; // -"-
296 G4double c8 = 6. * lamu * lamu; // -"- << 304 G4double c8=6.*lamu*lamu; // -"-
297 << 305
298 G4double xgi[8] = {.0199, .1017, .2372, .408 << 306 G4double xgi[8]={.0199,.1017,.2372,.4083,.5917,.7628,.8983,.9801}; // Gauss, 8
299 G4double wgi[8] = {.0506, .1112, .1569, .181 << 307 G4double wgi[8]={.0506,.1112,.1569,.1813,.1813,.1569,.1112,.0506}; // Gauss, 8
300 bbbtf = 183.; // for the moment... << 308 bbbtf=183.; // for the moment...
301 bbbh = 202.4; // for the moment... << 309 bbbh=202.4; // for the moment...
302 g1tf = 1.95e-5; << 310 g1tf=1.95e-5;
303 g2tf = 5.3e-5; << 311 g2tf=5.3e-5;
304 g1h = 4.4e-5; << 312 g1h=4.4e-5;
305 g2h = 4.8e-5; << 313 g2h=4.8e-5;
306 << 314
307 e = tkin + lamu; << 315 e=tkin+lamu;
308 z13 = fNist->GetZ13(G4lrint(z)); << 316 z13=pow(z,0.3333333);
309 e1 = e - ep; << 317 e1=e-ep;
310 crp_g4 = 0.; << 318 crp_g4=0.;
311 if (e1 <= c3 * z13) return crp_g4; // ep > << 319 if(e1 <= c3*z13) return crp_g4; // ep > max
312 alf = c7 / ep; // 4m/ep << 320 alf=c7/ep; // 4m/ep
313 a3 = 1. - alf; << 321 a3=1.-alf;
314 if (a3 <= 0.) return crp_g4; // ep < min << 322 if(a3 <= 0.) return crp_g4; // ep < min
315 //*** zeta calculation << 323 //*** zeta calculation
316 if (z <= 1.5) // special case of hidrogen << 324 if(z <= 1.5) // special case of hidrogen
317 { << 325 {
318 bbb = bbbh; << 326 bbb=bbbh;
319 g1 = g1h; << 327 g1=g1h;
320 g2 = g2h; << 328 g2=g2h;
321 } << 329 }
322 else { << 330 else
323 bbb = bbbtf; << 331 {
324 g1 = g1tf; << 332 bbb=bbbtf;
325 g2 = g2tf; << 333 g1=g1tf;
326 } << 334 g2=g2tf;
327 zeta1 = 0.073 * G4Log(e / (lamu + g1 * z13 * << 335 }
328 if (zeta1 > 0.) { << 336 zeta1=0.073*log(e/(lamu+g1*z13*z13*e))-0.26;
329 zeta2 = 0.058 * log(e / (lamu + g2 * z13 * << 337 if(zeta1 > 0.)
330 zeta = zeta1 / zeta2; << 338 {
331 } << 339 zeta2=0.058*log(e/(lamu+g2*z13*e))-0.14;
332 else { << 340 zeta=zeta1/zeta2;
333 zeta = 0.; << 341 }
334 } << 342 else
335 z2 = z * (z + zeta); // << 343 {
336 //*** just to check (for comp << 344 zeta=0.;
337 // z2=z*(z+1.) << 345 }
338 // bbb=189. << 346 z2=z*(z+zeta); //
339 //*** << 347 //*** just to check (for comparison with crp_m)
340 screen0 = 2. * ame * sqrte * bbb / (z13 * ep << 348 // z2=z*(z+1.)
341 a0 = e * e1; << 349 // bbb=189.
342 a1 = ep * ep / a0; // 2*beta << 350 //***
343 bet = 0.5 * a1; // beta << 351 screen0=2.*ame*sqrte*bbb/(z13*ep); // be careful with "ame"
344 xi0 = 0.25 * rmass * rmass * a1; // xi0 << 352 a0=e*e1;
345 del = c8 / a0; // 6mu^2/EE' << 353 a1=ep*ep/a0; // 2*beta
346 tmn = G4Log((alf + 2. * del * a3) / (1. + (1 << 354 bet=0.5*a1; // beta
347 sum = 0.; << 355 xi0=0.25*rmass*rmass*a1; // xi0
348 for (G4int i = 0; i < 8; ++i) { << 356 del=c8/a0; // 6mu^2/EE'
349 a4 = G4Exp(tmn * xgi[i]); // 1-r << 357 tmn=log((alf+2.*del*a3)/(1.+(1.-del)*sqrt(a3))); // log(1-rmax)
350 a5 = a4 * (2. - a4); // 1-r2 << 358 sum=0.;
351 a6 = 1. - a5; // r2 << 359 for(int i=0; i<=7; i++) // integration
352 a7 = 1. + a6; // 1+r2 << 360 {
353 a9 = 3. + a6; // 3+r2 << 361 a4=exp(tmn*xgi[i]); // 1-r
354 xi = xi0 * a5; << 362 a5=a4*(2.-a4); // 1-r2
355 xii = 1. / xi; << 363 a6=1.-a5; // r2
356 xi1 = 1. + xi; << 364 a7=1.+a6; // 1+r2
357 screen = screen0 * xi1 / a5; << 365 a9=3.+a6; // 3+r2
358 yeu = 5. - a6 + 4. * bet * a7; << 366 xi=xi0*a5;
359 yed = 2. * (1. + 3. * bet) * G4Log(3. + xi << 367 xii=1./xi;
360 ye1 = 1. + yeu / yed; << 368 xi1=1.+xi;
361 ale = G4Log(bbb / z13 * std::sqrt(xi1 * ye << 369 screen=screen0*xi1/a5;
362 cre = 0.5 * G4Log(1. + (1.5 / rmass * z13) << 370 yeu=5.-a6+4.*bet*a7;
363 if (xi <= 1e3) { << 371 yed=2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6);
364 be = ((2. + a6) * (1. + bet) + xi * a9) << 372 ye1=1.+yeu/yed;
365 } << 373 ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1));
366 else { << 374 cre=0.5*log(1.+(1.5/rmass*z13)*(1.5/rmass*z13)*xi1*ye1);
367 be = (3. - a6 + a1 * a7) / (2. * xi); / << 375 if(xi <= 1e3) //
368 } << 376 {
369 fe = std::max(0., (ale - cre) * be); << 377 be=((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9;
370 ymu = 4. + a6 + 3. * bet * a7; << 378 }
371 ymd = a7 * (1.5 + a1) * G4Log(3. + xi) + 1 << 379 else
372 ym1 = 1. + ymu / ymd; << 380 {
373 alm_crm = G4Log(bbb * rmass / (1.5 * z13 * << 381 be=(3.-a6+a1*a7)/(2.*xi); // -(6.-5.*a6+3.*bet*a6)/(6.*xi*xi);
374 if (xi >= 1e-3) { << 382 }
375 a10 = (1. + a1) * a5; // (1+2b)(1-r2) << 383 fe=max(0.,(ale-cre)*be);
376 bm = (a7 * (1. + 1.5 * bet) - a10 * xii) << 384 ymu=4.+a6+3.*bet*a7;
377 } << 385 ymd=a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6;
378 else { << 386 ym1=1.+ymu/ymd;
379 bm = (5. - a6 + bet * a9) * (xi / 2.); << 387 alm_crm=log(bbb*rmass/(1.5*z13*z13*(1.+screen*ym1)));
380 } << 388 if(xi >= 1e-3) //
381 fm = max(0., (alm_crm)*bm); << 389 {
382 //*** << 390 a10=(1.+a1)*a5; // (1+2b)(1-r2)
383 sum = sum + a4 * (fe + fm / (rmass * rmass << 391 bm=(a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10;
384 } << 392 }
385 crp_g4 = -tmn * sum * (z2 / a) * coeff * e1 << 393 else
386 return crp_g4; << 394 {
387 } << 395 bm=(5.-a6+bet*a9)*(xi/2.); // -(11.-5.*a6+.5*bet*(5.+a6))*(xi*xi/6.)
388 << 396 }
389 //....oooOO0OOooo........oooOO0OOooo........oo << 397 fm=max(0.,(alm_crm)*bm);
390 << 398 //***
391 G4double MuCrossSections::CRM_Mephi(G4double Z << 399 sum=sum+a4*(fe+fm/(rmass*rmass))*wgi[i];
392 { << 400 }
393 /* << 401 crp_g4=-tmn*sum*(z2/a)*coeff*e1/(e*ep);
394 ||Cross section for pair production of m << 402 return crp_g4;
395 ||By Siddharth Yajaman 19-07-2022 <<
396 ||Based on the formulae from Kelner, Kok <<
397 ||Physics of Atomic Nuclei, Vol. 63, No. <<
398 ||Equations (15) - (22) <<
399 */ <<
400 <<
401 const G4double xgi[] = {0.0198550717512320, <<
402 0.4082826787521750, <<
403 0.8983332387068135, <<
404 <<
405 const G4double wgi[] = {0.0506142681451880, <<
406 0.1813418916891810, <<
407 0.1111905172266870, <<
408 <<
409 static const G4double factorForCross = <<
410 2. / (3 * CLHEP::pi) <<
411 * pow(CLHEP::fine_structure_const * CLHEP: <<
412 <<
413 if (pairEnergy <= 2. * fMuonMass) return 0.0 <<
414 <<
415 G4double totalEnergy = tkin + fMuonMass; <<
416 G4double residEnergy = totalEnergy - pairEne <<
417 <<
418 if (residEnergy <= fMuonMass) return 0.0; <<
419 <<
420 G4double a0 = 1.0 / (totalEnergy * residEner <<
421 G4double rhomax = 1.0 - 2 * fMuonMass / pair <<
422 G4double tmnexp = 1. - rhomax; <<
423 <<
424 if (tmnexp >= 1.0) { <<
425 return 0.0; <<
426 } <<
427 <<
428 G4double tmn = G4Log(tmnexp); <<
429 <<
430 G4double z2 = Z * Z; <<
431 G4double beta = 0.5 * pairEnergy * pairEnerg <<
432 G4double xi0 = 0.5 * beta; <<
433 <<
434 // Gaussian integration in ln(1-ro) ( with 8 <<
435 G4double rho[8]; <<
436 G4double rho2[8]; <<
437 G4double xi[8]; <<
438 G4double xi1[8]; <<
439 G4double xii[8]; <<
440 <<
441 for (G4int i = 0; i < 8; ++i) { <<
442 rho[i] = G4Exp(tmn * xgi[i]) - 1.0; // rh <<
443 rho2[i] = rho[i] * rho[i]; <<
444 xi[i] = xi0 * (1.0 - rho2[i]); <<
445 xi1[i] = 1.0 + xi[i]; <<
446 xii[i] = 1.0 / xi[i]; <<
447 } <<
448 <<
449 G4double ximax = xi0 * (1. - rhomax * rhomax <<
450 <<
451 G4double Y = 10 * sqrt(fMuonMass / totalEner <<
452 G4double U[8]; <<
453 <<
454 for (G4int i = 0; i < 8; ++i) { <<
455 U[i] = U_func(Z, rho2[i], xi[i], Y, pairEn <<
456 } <<
457 <<
458 G4double UMax = U_func(Z, rhomax * rhomax, x <<
459 <<
460 G4double sum = 0.0; <<
461 <<
462 for (G4int i = 0; i < 8; ++i) { <<
463 G4double X = 1 + U[i] - UMax; <<
464 G4double lnX = G4Log(X); <<
465 G4double phi = ((2 + rho2[i]) * (1 + beta) <<
466 - 3 * rho2[i] + beta * (1 - <<
467 + ((1 + rho2[i]) * (1 + 1.5 <<
468 * G4Log(xi1[i]); <<
469 sum += wgi[i] * (1.0 + rho[i]) * phi * lnX <<
470 } <<
471 <<
472 return -tmn * sum * factorForCross * z2 * re <<
473 } 403 }
474 <<
475 //....oooOO0OOooo........oooOO0OOooo........oo 404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
476 405
477 G4double MuCrossSections::U_func(G4double Z, G <<
478 G4double pair <<
479 { <<
480 G4int z = G4lrint(Z); <<
481 G4double A27 = fNist->GetA27(z); <<
482 G4double Z13 = fNist->GetZ13(z); <<
483 static const G4double sqe = 2 * std::sqrt(G4 <<
484 G4double res = (0.65 * B / (A27 * Z13) * fMu <<
485 / (1 <<
486 + (sqe * B * (1 + xi) * (1 <<
487 / (Z13 * CLHEP::electr <<
488 return res; <<
489 } <<
490 406
491 //....oooOO0OOooo........oooOO0OOooo........oo <<
492 407