| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file electromagnetic/TestEm17/src/MuCross
27 /// \brief Implementation of the MuCrossSectio
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "MuCrossSections.hh"
34
35 #include "G4DataVector.hh"
36 #include "G4Exp.hh"
37 #include "G4Log.hh"
38 #include "G4Material.hh"
39 #include "G4MuonMinus.hh"
40 #include "G4NistManager.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 using namespace std;
47
48 MuCrossSections::MuCrossSections()
49 {
50 fNist = G4NistManager::Instance();
51 fMuonMass = G4MuonMinus::MuonMinus()->GetPDG
52 fMueRatio = fMuonMass / CLHEP::electron_mass
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 G4double MuCrossSections::CR_Macroscopic(const
58 G4dou
59 // return the macroscopic cross section (1/L)
60 {
61 const G4ElementVector* theElementVector = ma
62 const G4double* NbOfAtomsPerVolume = materia
63
64 G4double SIGMA = 0.;
65 G4int nelm = material->GetNumberOfElements()
66 for (G4int i = 0; i < nelm; ++i) {
67 const G4Element* element = (*theElementVec
68 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAto
69 }
70 return SIGMA;
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oo
74
75 G4double MuCrossSections::CR_PerAtom(const G4S
76 G4double
77 {
78 G4double z = element->GetZ();
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
91 else if (process == "muPairProd")
92 sigma = CRP_Mephi(z, a / (g / mole), tkin
93
94 else if (process == "muToMuonPairProd")
95 sigma = CRM_Mephi(z, tkin, ep);
96
97 return sigma;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
102 G4double MuCrossSections::CRB_Mephi(G4double z
103
104 //********************************************
105 //*** crb_g4_1.inc in comparison
106 //*** changes are introduced (September
107 //*** 1) function crb_g4 (Z,A,T
108 //*** 2) special case of hidrog
109 //*** Numerical comparison: 5 d
110 //***
111 //*** Cross section for bremsstrahlung
112 //*** By R.P.Kokoulin, September 1998
113 //*** Formulae from Kelner,Kokoulin,Pet
114 //*** (7,18,19,20,21,25,26); Dn (18) is
115 //*** Bugaev's inelatic nuclear correct
116 //********************************************
117 {
118 // G4double Z,A,Tkin,EP;
119 G4double crb_g4;
120 G4double e, v, delta, rab0, z_13, dn, b, b1,
121 //
122 G4double ame = 0.51099907e-3; // GeV
123 G4double lamu = 0.105658389; // GeV
124 G4double re = 2.81794092e-13; // cm
125 G4double avno = 6.022137e23;
126 G4double alpha = 1. / 137.036;
127 G4double rmass = lamu / ame; // "207"
128 G4double coeff = 16. / 3. * alpha * avno * (
129 G4double sqrte = 1.64872127; // sqrt(2.7182
130 G4double btf = 183.;
131 G4double btf1 = 1429.;
132 G4double bh = 202.4;
133 G4double bh1 = 446.;
134 //***
135 if (ep >= tkin) {
136 return 0.;
137 }
138 e = tkin + lamu;
139 v = ep / e;
140 delta = lamu * lamu * v / (2. * (e - ep));
141 rab0 = delta * sqrte;
142 G4int Z = G4lrint(z);
143 z_13 = 1. / fNist->GetZ13(z); //
144 //*** nuclear size and excitation, scr
145 dn = 1.54 * fNist->GetA27(Z);
146 if (z <= 1.5) // special case for hydrogen
147 {
148 b = bh;
149 b1 = bh1;
150 dn_star = dn;
151 }
152 else {
153 b = btf;
154 b1 = btf1;
155 dn_star = pow(dn, (1. - 1. / z)); // with
156 }
157 //*** nucleus contribution lo
158 rab1 = b * z_13;
159 fn = G4Log(rab1 / (dn_star * (ame + rab0 * r
160 if (fn < 0.) fn = 0.;
161 //*** electron contribution l
162 epmax1 = e / (1. + lamu * rmass / (2. * e));
163 if (ep >= epmax1) {
164 fe = 0.;
165 }
166 else {
167 rab2 = b1 * z_13 * z_13;
168 fe = G4Log(rab2 * lamu / ((1. + delta * rm
169 if (fe < 0.) fe = 0.;
170 }
171 crb_g4 = coeff * (1. - v * (1. - 0.75 * v))
172 return crb_g4;
173 }
174
175 //....oooOO0OOooo........oooOO0OOooo........oo
176
177 G4double MuCrossSections::CRK_Mephi(G4double z
178 //********************************************
179 //*** Cross section for knock-on electr
180 //*** (including bremsstrahlung e-diagr
181 //*** Units: cm^2/(g*GeV); Tkin, ep - G
182 //*** By R.P.Kokoulin, October 1998
183 //*** Formulae from Kelner,Kokoulin,Pet
184 //*** (a bit simplified Kelner's versio
185 //***
186 {
187 // G4double Z,A,Tkin,EP;
188 G4double crk_g4;
189 G4double v, sigma0, a1, a3;
190 //
191 G4double ame = 0.51099907e-3; // GeV
192 G4double lamu = 0.105658389; // GeV
193 G4double re = 2.81794092e-13; // cm
194 G4double avno = 6.022137e23;
195 G4double alpha = 1. / 137.036;
196 G4double lpi = 3.141592654;
197 G4double bmu = lamu * lamu / (2. * ame);
198 G4double coeff0 = avno * 2. * lpi * ame * re
199 G4double coeff1 = alpha / (2. * lpi);
200 //***
201
202 G4double e = tkin + lamu;
203 G4double epmax = e / (1. + bmu / e);
204 if (ep >= epmax) {
205 return 0.0;
206 }
207 v = ep / e;
208 sigma0 = coeff0 * (z / a) * (1. - ep / epmax
209 a1 = G4Log(1. + 2. * ep / ame);
210 a3 = G4Log(4. * e * (e - ep) / (lamu * lamu)
211 crk_g4 = sigma0 * (1. + coeff1 * a1 * (a3 -
212 return crk_g4;
213 }
214
215 //....oooOO0OOooo........oooOO0OOooo........oo
216
217 G4double MuCrossSections::CRN_Mephi(G4double,
218
219 //********************************************
220 //*** Differential cross section for ph
221 //*** Formulae from Borog & Petrukhin,
222 //*** Real photon cross section: Caldwe
223 //*** Nuclear shadowing: Brodsky e.a.,
224 //*** Units: cm^2 / g GeV.
225 //*** CRN_G4_1.inc January 31st,
226 //********************************************
227 {
228 // G4double Z,A,Tkin,EP;
229 G4double crn_g4;
230 G4double e, aeff, sigph, v, v1, v2, amu2, up
231 //***
232 G4double lamu = 0.105658389; // GeV
233 G4double avno = 6.022137e23;
234 G4double amp = 0.9382723; // GeV
235 G4double lpi = 3.14159265;
236 G4double alpha = 1. / 137.036;
237 //***
238 G4double epmin_phn = 0.20; // GeV
239 G4double alam2 = 0.400000; // GeV**2
240 G4double alam = 0.632456; // sqrt(alam2)
241 G4double coeffn = alpha / lpi * avno * 1e-30
242 //***
243 e = tkin + lamu;
244 crn_g4 = 0.;
245 if (ep >= e - 0.5 * amp || ep <= epmin_phn)
246 aeff = 0.22 * a + 0.78 * pow(a, 0.89); // s
247 sigph = 49.2 + 11.1 * G4Log(ep) + 151.8 / st
248 v = ep / e;
249 v1 = 1. - v;
250 v2 = v * v;
251 amu2 = lamu * lamu;
252 up = e * e * v1 / amu2 * (1. + amu2 * v2 / (
253 down = 1. + ep / alam * (1. + alam / (2. * a
254 crn_g4 = coeffn * aeff / a * sigph / ep
255 * (-v1 + (v1 + 0.5 * v2 * (1. + 2.
256 if (crn_g4 < 0.) crn_g4 = 0.;
257 return crn_g4;
258 }
259
260 //....oooOO0OOooo........oooOO0OOooo........oo
261
262 G4double MuCrossSections::CRP_Mephi(G4double z
263
264 //********************************************
265 //*** crp_g4_1.inc in comparison
266 //*** changes are introduced (January 1
267 //*** 1) Avno/A, cm^2/gram GeV
268 //*** 2) zeta_loss(E,Z)
269 //*** 3) function crp_g4 (Z,A,T
270 //*** 4) bbb=183 (Thomas
271 //*** 5) special case of hidrog
272 //*** 6) expansions in 'xi' are
273 //***
274 //*** Cross section for electron pair p
275 //*** By R.P.Kokoulin, December 1997
276 //*** Formulae from Kokoulin & Petrukhi
277 {
278 // G4double Z,A,Tkin,EP;
279 G4double crp_g4;
280 G4double bbbtf, bbbh, g1tf, g2tf, g1h, g2h,
281 G4double g1, g2, zeta1, zeta2, zeta, z2, scr
282 G4double tmn, sum, a4, a5, a6, a7, a9, xi, x
283 G4double ale, cre, be, fe, ymu, ymd, ym1, al
284 //
285 G4double ame = 0.51099907e-3; // GeV
286 G4double lamu = 0.105658389; // GeV
287 G4double re = 2.81794092e-13; // cm
288 G4double avno = 6.022137e23;
289 G4double lpi = 3.14159265;
290 G4double alpha = 1. / 137.036;
291 G4double rmass = lamu / ame; // "207"
292 G4double coeff = 4. / (3. * lpi) * (alpha *
293 G4double sqrte = 1.64872127; // sqrt(2.7182
294 G4double c3 = 3. * sqrte * lamu / 4.; // fo
295 G4double c7 = 4. * ame; // -"-
296 G4double c8 = 6. * lamu * lamu; // -"-
297
298 G4double xgi[8] = {.0199, .1017, .2372, .408
299 G4double wgi[8] = {.0506, .1112, .1569, .181
300 bbbtf = 183.; // for the moment...
301 bbbh = 202.4; // for the moment...
302 g1tf = 1.95e-5;
303 g2tf = 5.3e-5;
304 g1h = 4.4e-5;
305 g2h = 4.8e-5;
306
307 e = tkin + lamu;
308 z13 = fNist->GetZ13(G4lrint(z));
309 e1 = e - ep;
310 crp_g4 = 0.;
311 if (e1 <= c3 * z13) return crp_g4; // ep >
312 alf = c7 / ep; // 4m/ep
313 a3 = 1. - alf;
314 if (a3 <= 0.) return crp_g4; // ep < min
315 //*** zeta calculation
316 if (z <= 1.5) // special case of hidrogen
317 {
318 bbb = bbbh;
319 g1 = g1h;
320 g2 = g2h;
321 }
322 else {
323 bbb = bbbtf;
324 g1 = g1tf;
325 g2 = g2tf;
326 }
327 zeta1 = 0.073 * G4Log(e / (lamu + g1 * z13 *
328 if (zeta1 > 0.) {
329 zeta2 = 0.058 * log(e / (lamu + g2 * z13 *
330 zeta = zeta1 / zeta2;
331 }
332 else {
333 zeta = 0.;
334 }
335 z2 = z * (z + zeta); //
336 //*** just to check (for comp
337 // z2=z*(z+1.)
338 // bbb=189.
339 //***
340 screen0 = 2. * ame * sqrte * bbb / (z13 * ep
341 a0 = e * e1;
342 a1 = ep * ep / a0; // 2*beta
343 bet = 0.5 * a1; // beta
344 xi0 = 0.25 * rmass * rmass * a1; // xi0
345 del = c8 / a0; // 6mu^2/EE'
346 tmn = G4Log((alf + 2. * del * a3) / (1. + (1
347 sum = 0.;
348 for (G4int i = 0; i < 8; ++i) {
349 a4 = G4Exp(tmn * xgi[i]); // 1-r
350 a5 = a4 * (2. - a4); // 1-r2
351 a6 = 1. - a5; // r2
352 a7 = 1. + a6; // 1+r2
353 a9 = 3. + a6; // 3+r2
354 xi = xi0 * a5;
355 xii = 1. / xi;
356 xi1 = 1. + xi;
357 screen = screen0 * xi1 / a5;
358 yeu = 5. - a6 + 4. * bet * a7;
359 yed = 2. * (1. + 3. * bet) * G4Log(3. + xi
360 ye1 = 1. + yeu / yed;
361 ale = G4Log(bbb / z13 * std::sqrt(xi1 * ye
362 cre = 0.5 * G4Log(1. + (1.5 / rmass * z13)
363 if (xi <= 1e3) {
364 be = ((2. + a6) * (1. + bet) + xi * a9)
365 }
366 else {
367 be = (3. - a6 + a1 * a7) / (2. * xi); /
368 }
369 fe = std::max(0., (ale - cre) * be);
370 ymu = 4. + a6 + 3. * bet * a7;
371 ymd = a7 * (1.5 + a1) * G4Log(3. + xi) + 1
372 ym1 = 1. + ymu / ymd;
373 alm_crm = G4Log(bbb * rmass / (1.5 * z13 *
374 if (xi >= 1e-3) {
375 a10 = (1. + a1) * a5; // (1+2b)(1-r2)
376 bm = (a7 * (1. + 1.5 * bet) - a10 * xii)
377 }
378 else {
379 bm = (5. - a6 + bet * a9) * (xi / 2.);
380 }
381 fm = max(0., (alm_crm)*bm);
382 //***
383 sum = sum + a4 * (fe + fm / (rmass * rmass
384 }
385 crp_g4 = -tmn * sum * (z2 / a) * coeff * e1
386 return crp_g4;
387 }
388
389 //....oooOO0OOooo........oooOO0OOooo........oo
390
391 G4double MuCrossSections::CRM_Mephi(G4double Z
392 {
393 /*
394 ||Cross section for pair production of m
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 }
474
475 //....oooOO0OOooo........oooOO0OOooo........oo
476
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
491 //....oooOO0OOooo........oooOO0OOooo........oo
492