Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 /// \file electromagnetic/TestEm17/src/MuCrossSections.cc 27 /// \brief Implementation of the MuCrossSections class 28 // 29 // 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 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........oooOO0OOooo........oooOO0OOooo...... 45 46 using namespace std; 47 48 MuCrossSections::MuCrossSections() 49 { 50 fNist = G4NistManager::Instance(); 51 fMuonMass = G4MuonMinus::MuonMinus()->GetPDGMass(); 52 fMueRatio = fMuonMass / CLHEP::electron_mass_c2; 53 } 54 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 57 G4double MuCrossSections::CR_Macroscopic(const G4String& process, const G4Material* material, 58 G4double tkin, G4double ep) 59 // return the macroscopic cross section (1/L) in GEANT4 internal units 60 { 61 const G4ElementVector* theElementVector = material->GetElementVector(); 62 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume(); 63 64 G4double SIGMA = 0.; 65 G4int nelm = material->GetNumberOfElements(); 66 for (G4int i = 0; i < nelm; ++i) { 67 const G4Element* element = (*theElementVector)[i]; 68 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAtom(process, element, tkin, ep); 69 } 70 return SIGMA; 71 } 72 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 74 75 G4double MuCrossSections::CR_PerAtom(const G4String& process, const G4Element* element, 76 G4double tkin, G4double ep) 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 / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro; 83 84 else if (process == "muIoni") 85 sigma = CRK_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro; 86 87 // else if (process == "muNucl") 88 else if (process == "muonNuclear") 89 sigma = CRN_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro; 90 91 else if (process == "muPairProd") 92 sigma = CRP_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro; 93 94 else if (process == "muToMuonPairProd") 95 sigma = CRM_Mephi(z, tkin, ep); 96 97 return sigma; 98 } 99 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 101 102 G4double MuCrossSections::CRB_Mephi(G4double z, G4double a, G4double tkin, G4double ep) 103 104 //*********************************************************************** 105 //*** crb_g4_1.inc in comparison with crb_.inc, following 106 //*** changes are introduced (September 24th, 1998): 107 //*** 1) function crb_g4 (Z,A,Tkin,EP), Tkin is kinetic energy 108 //*** 2) special case of hidrogen (Z<1.5; b,b1,Dn_star) 109 //*** Numerical comparison: 5 decimal digits coincide. 110 //*** 111 //*** Cross section for bremsstrahlung by fast muon 112 //*** By R.P.Kokoulin, September 1998 113 //*** Formulae from Kelner,Kokoulin,Petrukhin 1995, Preprint MEPhI 114 //*** (7,18,19,20,21,25,26); Dn (18) is modified to incorporate 115 //*** Bugaev's inelatic nuclear correction (28) for Z > 1. 116 //*********************************************************************** 117 { 118 // G4double Z,A,Tkin,EP; 119 G4double crb_g4; 120 G4double e, v, delta, rab0, z_13, dn, b, b1, dn_star, rab1, fn, epmax1, fe, rab2; 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 * (re / rmass) * (re / rmass); // cm^2 129 G4double sqrte = 1.64872127; // sqrt(2.71828...) 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)); // qmin 141 rab0 = delta * sqrte; 142 G4int Z = G4lrint(z); 143 z_13 = 1. / fNist->GetZ13(z); // 144 //*** nuclear size and excitation, screening parameters 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 Bugaev's correction 156 } 157 //*** nucleus contribution logarithm 158 rab1 = b * z_13; 159 fn = G4Log(rab1 / (dn_star * (ame + rab0 * rab1)) * (lamu + delta * (dn_star * sqrte - 2.))); 160 if (fn < 0.) fn = 0.; 161 //*** electron contribution logarithm 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 * rmass / (ame * sqrte)) * (ame + rab0 * rab2))); 169 if (fe < 0.) fe = 0.; 170 } 171 crb_g4 = coeff * (1. - v * (1. - 0.75 * v)) * z * (z * fn + fe) / (a * ep); 172 return crb_g4; 173 } 174 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 176 177 G4double MuCrossSections::CRK_Mephi(G4double z, G4double a, G4double tkin, G4double ep) 178 //*********************************************************************** 179 //*** Cross section for knock-on electron production by fast muons 180 //*** (including bremsstrahlung e-diagrams and rad. correction). 181 //*** Units: cm^2/(g*GeV); Tkin, ep - GeV. 182 //*** By R.P.Kokoulin, October 1998 183 //*** Formulae from Kelner,Kokoulin,Petrukhin, Phys.Atom.Nuclei, 1997 184 //*** (a bit simplified Kelner's version of Eq.30 - with 2 logarithms). 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 * 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 + 0.5 * v * v) / (ep * ep); 209 a1 = G4Log(1. + 2. * ep / ame); 210 a3 = G4Log(4. * e * (e - ep) / (lamu * lamu)); 211 crk_g4 = sigma0 * (1. + coeff1 * a1 * (a3 - a1)); 212 return crk_g4; 213 } 214 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 216 217 G4double MuCrossSections::CRN_Mephi(G4double, G4double a, G4double tkin, G4double ep) 218 219 //*********************************************************************** 220 //*** Differential cross section for photonuclear muon interaction. 221 //*** Formulae from Borog & Petrukhin, 1975 222 //*** Real photon cross section: Caldwell e.a., 1979 223 //*** Nuclear shadowing: Brodsky e.a., 1972 224 //*** Units: cm^2 / g GeV. 225 //*** CRN_G4_1.inc January 31st, 1998 R.P.Kokoulin 226 //*********************************************************************** 227 { 228 // G4double Z,A,Tkin,EP; 229 G4double crn_g4; 230 G4double e, aeff, sigph, v, v1, v2, amu2, up, down; 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; // cm^2/microbarn 242 //*** 243 e = tkin + lamu; 244 crn_g4 = 0.; 245 if (ep >= e - 0.5 * amp || ep <= epmin_phn) return crn_g4; 246 aeff = 0.22 * a + 0.78 * pow(a, 0.89); // shadowing 247 sigph = 49.2 + 11.1 * G4Log(ep) + 151.8 / std::sqrt(ep); // microbarn 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 / (alam2 * v1)); 253 down = 1. + ep / alam * (1. + alam / (2. * amp) + ep / alam); 254 crn_g4 = coeffn * aeff / a * sigph / ep 255 * (-v1 + (v1 + 0.5 * v2 * (1. + 2. * amu2 / alam2)) * log(up / down)); 256 if (crn_g4 < 0.) crn_g4 = 0.; 257 return crn_g4; 258 } 259 260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 261 262 G4double MuCrossSections::CRP_Mephi(G4double z, G4double a, G4double tkin, G4double ep) 263 264 //********************************************************************** 265 //*** crp_g4_1.inc in comparison with crp_m.inc, following 266 //*** changes are introduced (January 16th, 1998): 267 //*** 1) Avno/A, cm^2/gram GeV 268 //*** 2) zeta_loss(E,Z) from Kelner 1997, Eqs.(53-54) 269 //*** 3) function crp_g4 (Z,A,Tkin,EP), Tkin is kinetic energy 270 //*** 4) bbb=183 (Thomas-Fermi) 271 //*** 5) special case of hidrogen (Z<1.5), bbb,g1,g2 272 //*** 6) expansions in 'xi' are simplified (Jan.17th,1998) 273 //*** 274 //*** Cross section for electron pair production by fast muon 275 //*** By R.P.Kokoulin, December 1997 276 //*** Formulae from Kokoulin & Petrukhin 1971, Hobart, Eqs.(1,8,9,10) 277 { 278 // G4double Z,A,Tkin,EP; 279 G4double crp_g4; 280 G4double bbbtf, bbbh, g1tf, g2tf, g1h, g2h, e, z13, e1, alf, a3, bbb; 281 G4double g1, g2, zeta1, zeta2, zeta, z2, screen0, a0, a1, bet, xi0, del; 282 G4double tmn, sum, a4, a5, a6, a7, a9, xi, xii, xi1, screen, yeu, yed, ye1; 283 G4double ale, cre, be, fe, ymu, ymd, ym1, alm_crm, a10, bm, fm; 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 * re) * (alpha * re) * avno; // cm^2 293 G4double sqrte = 1.64872127; // sqrt(2.71828...) 294 G4double c3 = 3. * sqrte * lamu / 4.; // for limits 295 G4double c7 = 4. * ame; // -"- 296 G4double c8 = 6. * lamu * lamu; // -"- 297 298 G4double xgi[8] = {.0199, .1017, .2372, .4083, .5917, .7628, .8983, .9801}; // Gauss, 8 299 G4double wgi[8] = {.0506, .1112, .1569, .1813, .1813, .1569, .1112, .0506}; // Gauss, 8 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 > max 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 * z13 * e)) - 0.26; 328 if (zeta1 > 0.) { 329 zeta2 = 0.058 * log(e / (lamu + g2 * z13 * e)) - 0.14; 330 zeta = zeta1 / zeta2; 331 } 332 else { 333 zeta = 0.; 334 } 335 z2 = z * (z + zeta); // 336 //*** just to check (for comparison with crp_m) 337 // z2=z*(z+1.) 338 // bbb=189. 339 //*** 340 screen0 = 2. * ame * sqrte * bbb / (z13 * ep); // be careful with "ame" 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. - del) * sqrt(a3))); // log(1-rmax) 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. + xii) - a6 - a1 * (2. - a6); 360 ye1 = 1. + yeu / yed; 361 ale = G4Log(bbb / z13 * std::sqrt(xi1 * ye1) / (1. + screen * ye1)); 362 cre = 0.5 * G4Log(1. + (1.5 / rmass * z13) * (1.5 / rmass * z13) * xi1 * ye1); 363 if (xi <= 1e3) { 364 be = ((2. + a6) * (1. + bet) + xi * a9) * G4Log(1. + xii) + (a5 - bet) / xi1 - a9; 365 } 366 else { 367 be = (3. - a6 + a1 * a7) / (2. * xi); // -(6.-5.*a6+3.*bet*a6)/(6.*xi*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. - 1.5 * a6; 372 ym1 = 1. + ymu / ymd; 373 alm_crm = G4Log(bbb * rmass / (1.5 * z13 * z13 * (1. + screen * ym1))); 374 if (xi >= 1e-3) { 375 a10 = (1. + a1) * a5; // (1+2b)(1-r2) 376 bm = (a7 * (1. + 1.5 * bet) - a10 * xii) * G4Log(xi1) + xi * (a5 - bet) / xi1 + a10; 377 } 378 else { 379 bm = (5. - a6 + bet * a9) * (xi / 2.); // -(11.-5.*a6+.5*bet*(5.+a6))*(xi*xi/6.) 380 } 381 fm = max(0., (alm_crm)*bm); 382 //*** 383 sum = sum + a4 * (fe + fm / (rmass * rmass)) * wgi[i]; 384 } 385 crp_g4 = -tmn * sum * (z2 / a) * coeff * e1 / (e * ep); 386 return crp_g4; 387 } 388 389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 390 391 G4double MuCrossSections::CRM_Mephi(G4double Z, G4double tkin, G4double pairEnergy) 392 { 393 /* 394 ||Cross section for pair production of muons by fast muon 395 ||By Siddharth Yajaman 19-07-2022 396 ||Based on the formulae from Kelner, Kokoulin and Petrukhin, 397 ||Physics of Atomic Nuclei, Vol. 63, No. 9, 2000, pp. 1603-1611 398 ||Equations (15) - (22) 399 */ 400 401 const G4double xgi[] = {0.0198550717512320, 0.1016667612931865, 0.2372337950418355, 402 0.4082826787521750, 0.5917173212478250, 0.7627662049581645, 403 0.8983332387068135, 0.9801449282487680}; 404 405 const G4double wgi[] = {0.0506142681451880, 0.1111905172266870, 0.1568533229389435, 406 0.1813418916891810, 0.1813418916891810, 0.1568533229389435, 407 0.1111905172266870, 0.0506142681451880}; 408 409 static const G4double factorForCross = 410 2. / (3 * CLHEP::pi) 411 * pow(CLHEP::fine_structure_const * CLHEP::classic_electr_radius / fMueRatio, 2); 412 413 if (pairEnergy <= 2. * fMuonMass) return 0.0; 414 415 G4double totalEnergy = tkin + fMuonMass; 416 G4double residEnergy = totalEnergy - pairEnergy; 417 418 if (residEnergy <= fMuonMass) return 0.0; 419 420 G4double a0 = 1.0 / (totalEnergy * residEnergy); 421 G4double rhomax = 1.0 - 2 * fMuonMass / pairEnergy; 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 * pairEnergy * a0; 432 G4double xi0 = 0.5 * beta; 433 434 // Gaussian integration in ln(1-ro) ( with 8 points) 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; // rho = -asymmetry 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 / totalEnergy); 452 G4double U[8]; 453 454 for (G4int i = 0; i < 8; ++i) { 455 U[i] = U_func(Z, rho2[i], xi[i], Y, pairEnergy); 456 } 457 458 G4double UMax = U_func(Z, rhomax * rhomax, ximax, Y, pairEnergy); 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) + xi[i] * (3 + rho2[i])) * G4Log(1 + xii[i]) - 1 466 - 3 * rho2[i] + beta * (1 - 2 * rho2[i]) 467 + ((1 + rho2[i]) * (1 + 1.5 * beta) - xii[i] * (1 + 2 * beta) * (1 - rho2[i])) 468 * G4Log(xi1[i]); 469 sum += wgi[i] * (1.0 + rho[i]) * phi * lnX; 470 } 471 472 return -tmn * sum * factorForCross * z2 * residEnergy / (totalEnergy * pairEnergy); 473 } 474 475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 476 477 G4double MuCrossSections::U_func(G4double Z, G4double rho2, G4double xi, G4double Y, 478 G4double pairEnergy, const G4double B) 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(G4Exp(1.0)) * fMuonMass * fMuonMass; 484 G4double res = (0.65 * B / (A27 * Z13) * fMueRatio) 485 / (1 486 + (sqe * B * (1 + xi) * (1 + Y)) 487 / (Z13 * CLHEP::electron_mass_c2 * pairEnergy * (1 - rho2))); 488 return res; 489 } 490 491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 492