File indexing completed on 2026-06-13 07:54:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 #include "MuCrossSections.hh"
0030
0031 #include "G4DataVector.hh"
0032 #include "G4Exp.hh"
0033 #include "G4Log.hh"
0034 #include "G4Material.hh"
0035 #include "G4MuonMinus.hh"
0036 #include "G4NistManager.hh"
0037 #include "G4PhysicalConstants.hh"
0038 #include "G4SystemOfUnits.hh"
0039
0040
0041
0042 using namespace std;
0043
0044 MuCrossSections::MuCrossSections()
0045 {
0046 fNist = G4NistManager::Instance();
0047 fMuonMass = G4MuonMinus::MuonMinus()->GetPDGMass();
0048 fMueRatio = fMuonMass / CLHEP::electron_mass_c2;
0049 }
0050
0051
0052
0053 G4double MuCrossSections::CR_Macroscopic(const G4String& process, const G4Material* material,
0054 G4double tkin, G4double ep)
0055
0056 {
0057 const G4ElementVector* theElementVector = material->GetElementVector();
0058 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
0059
0060 G4double SIGMA = 0.;
0061 G4int nelm = material->GetNumberOfElements();
0062 for (G4int i = 0; i < nelm; ++i) {
0063 const G4Element* element = (*theElementVector)[i];
0064 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAtom(process, element, tkin, ep);
0065 }
0066 return SIGMA;
0067 }
0068
0069
0070
0071 G4double MuCrossSections::CR_PerAtom(const G4String& process, const G4Element* element,
0072 G4double tkin, G4double ep)
0073 {
0074 G4double z = element->GetZ();
0075 G4double a = element->GetA();
0076 G4double sigma = 0.;
0077 if (process == "muBrems")
0078 sigma = CRB_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro;
0079
0080 else if (process == "muIoni")
0081 sigma = CRK_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro;
0082
0083
0084 else if (process == "muonNuclear")
0085 sigma = CRN_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro;
0086
0087 else if (process == "muPairProd")
0088 sigma = CRP_Mephi(z, a / (g / mole), tkin / GeV, ep / GeV) * (cm2 / (g * GeV)) * a / Avogadro;
0089
0090 else if (process == "muToMuonPairProd")
0091 sigma = CRM_Mephi(z, tkin, ep);
0092
0093 return sigma;
0094 }
0095
0096
0097
0098 G4double MuCrossSections::CRB_Mephi(G4double z, G4double a, G4double tkin, G4double ep)
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113 {
0114
0115 G4double crb_g4;
0116 G4double e, v, delta, rab0, z_13, dn, b, b1, dn_star, rab1, fn, epmax1, fe, rab2;
0117
0118 G4double ame = 0.51099907e-3;
0119 G4double lamu = 0.105658389;
0120 G4double re = 2.81794092e-13;
0121 G4double avno = 6.022137e23;
0122 G4double alpha = 1. / 137.036;
0123 G4double rmass = lamu / ame;
0124 G4double coeff = 16. / 3. * alpha * avno * (re / rmass) * (re / rmass);
0125 G4double sqrte = 1.64872127;
0126 G4double btf = 183.;
0127 G4double btf1 = 1429.;
0128 G4double bh = 202.4;
0129 G4double bh1 = 446.;
0130
0131 if (ep >= tkin) {
0132 return 0.;
0133 }
0134 e = tkin + lamu;
0135 v = ep / e;
0136 delta = lamu * lamu * v / (2. * (e - ep));
0137 rab0 = delta * sqrte;
0138 G4int Z = G4lrint(z);
0139 z_13 = 1. / fNist->GetZ13(z);
0140
0141 dn = 1.54 * fNist->GetA27(Z);
0142 if (z <= 1.5)
0143 {
0144 b = bh;
0145 b1 = bh1;
0146 dn_star = dn;
0147 }
0148 else {
0149 b = btf;
0150 b1 = btf1;
0151 dn_star = pow(dn, (1. - 1. / z));
0152 }
0153
0154 rab1 = b * z_13;
0155 fn = G4Log(rab1 / (dn_star * (ame + rab0 * rab1)) * (lamu + delta * (dn_star * sqrte - 2.)));
0156 if (fn < 0.) fn = 0.;
0157
0158 epmax1 = e / (1. + lamu * rmass / (2. * e));
0159 if (ep >= epmax1) {
0160 fe = 0.;
0161 }
0162 else {
0163 rab2 = b1 * z_13 * z_13;
0164 fe = G4Log(rab2 * lamu / ((1. + delta * rmass / (ame * sqrte)) * (ame + rab0 * rab2)));
0165 if (fe < 0.) fe = 0.;
0166 }
0167 crb_g4 = coeff * (1. - v * (1. - 0.75 * v)) * z * (z * fn + fe) / (a * ep);
0168 return crb_g4;
0169 }
0170
0171
0172
0173 G4double MuCrossSections::CRK_Mephi(G4double z, G4double a, G4double tkin, G4double ep)
0174
0175
0176
0177
0178
0179
0180
0181
0182 {
0183
0184 G4double crk_g4;
0185 G4double v, sigma0, a1, a3;
0186
0187 G4double ame = 0.51099907e-3;
0188 G4double lamu = 0.105658389;
0189 G4double re = 2.81794092e-13;
0190 G4double avno = 6.022137e23;
0191 G4double alpha = 1. / 137.036;
0192 G4double lpi = 3.141592654;
0193 G4double bmu = lamu * lamu / (2. * ame);
0194 G4double coeff0 = avno * 2. * lpi * ame * re * re;
0195 G4double coeff1 = alpha / (2. * lpi);
0196
0197
0198 G4double e = tkin + lamu;
0199 G4double epmax = e / (1. + bmu / e);
0200 if (ep >= epmax) {
0201 return 0.0;
0202 }
0203 v = ep / e;
0204 sigma0 = coeff0 * (z / a) * (1. - ep / epmax + 0.5 * v * v) / (ep * ep);
0205 a1 = G4Log(1. + 2. * ep / ame);
0206 a3 = G4Log(4. * e * (e - ep) / (lamu * lamu));
0207 crk_g4 = sigma0 * (1. + coeff1 * a1 * (a3 - a1));
0208 return crk_g4;
0209 }
0210
0211
0212
0213 G4double MuCrossSections::CRN_Mephi(G4double, G4double a, G4double tkin, G4double ep)
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223 {
0224
0225 G4double crn_g4;
0226 G4double e, aeff, sigph, v, v1, v2, amu2, up, down;
0227
0228 G4double lamu = 0.105658389;
0229 G4double avno = 6.022137e23;
0230 G4double amp = 0.9382723;
0231 G4double lpi = 3.14159265;
0232 G4double alpha = 1. / 137.036;
0233
0234 G4double epmin_phn = 0.20;
0235 G4double alam2 = 0.400000;
0236 G4double alam = 0.632456;
0237 G4double coeffn = alpha / lpi * avno * 1e-30;
0238
0239 e = tkin + lamu;
0240 crn_g4 = 0.;
0241 if (ep >= e - 0.5 * amp || ep <= epmin_phn) return crn_g4;
0242 aeff = 0.22 * a + 0.78 * pow(a, 0.89);
0243 sigph = 49.2 + 11.1 * G4Log(ep) + 151.8 / std::sqrt(ep);
0244 v = ep / e;
0245 v1 = 1. - v;
0246 v2 = v * v;
0247 amu2 = lamu * lamu;
0248 up = e * e * v1 / amu2 * (1. + amu2 * v2 / (alam2 * v1));
0249 down = 1. + ep / alam * (1. + alam / (2. * amp) + ep / alam);
0250 crn_g4 = coeffn * aeff / a * sigph / ep
0251 * (-v1 + (v1 + 0.5 * v2 * (1. + 2. * amu2 / alam2)) * log(up / down));
0252 if (crn_g4 < 0.) crn_g4 = 0.;
0253 return crn_g4;
0254 }
0255
0256
0257
0258 G4double MuCrossSections::CRP_Mephi(G4double z, G4double a, G4double tkin, G4double ep)
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273 {
0274
0275 G4double crp_g4;
0276 G4double bbbtf, bbbh, g1tf, g2tf, g1h, g2h, e, z13, e1, alf, a3, bbb;
0277 G4double g1, g2, zeta1, zeta2, zeta, z2, screen0, a0, a1, bet, xi0, del;
0278 G4double tmn, sum, a4, a5, a6, a7, a9, xi, xii, xi1, screen, yeu, yed, ye1;
0279 G4double ale, cre, be, fe, ymu, ymd, ym1, alm_crm, a10, bm, fm;
0280
0281 G4double ame = 0.51099907e-3;
0282 G4double lamu = 0.105658389;
0283 G4double re = 2.81794092e-13;
0284 G4double avno = 6.022137e23;
0285 G4double lpi = 3.14159265;
0286 G4double alpha = 1. / 137.036;
0287 G4double rmass = lamu / ame;
0288 G4double coeff = 4. / (3. * lpi) * (alpha * re) * (alpha * re) * avno;
0289 G4double sqrte = 1.64872127;
0290 G4double c3 = 3. * sqrte * lamu / 4.;
0291 G4double c7 = 4. * ame;
0292 G4double c8 = 6. * lamu * lamu;
0293
0294 G4double xgi[8] = {.0199, .1017, .2372, .4083, .5917, .7628, .8983, .9801};
0295 G4double wgi[8] = {.0506, .1112, .1569, .1813, .1813, .1569, .1112, .0506};
0296 bbbtf = 183.;
0297 bbbh = 202.4;
0298 g1tf = 1.95e-5;
0299 g2tf = 5.3e-5;
0300 g1h = 4.4e-5;
0301 g2h = 4.8e-5;
0302
0303 e = tkin + lamu;
0304 z13 = fNist->GetZ13(G4lrint(z));
0305 e1 = e - ep;
0306 crp_g4 = 0.;
0307 if (e1 <= c3 * z13) return crp_g4;
0308 alf = c7 / ep;
0309 a3 = 1. - alf;
0310 if (a3 <= 0.) return crp_g4;
0311
0312 if (z <= 1.5)
0313 {
0314 bbb = bbbh;
0315 g1 = g1h;
0316 g2 = g2h;
0317 }
0318 else {
0319 bbb = bbbtf;
0320 g1 = g1tf;
0321 g2 = g2tf;
0322 }
0323 zeta1 = 0.073 * G4Log(e / (lamu + g1 * z13 * z13 * e)) - 0.26;
0324 if (zeta1 > 0.) {
0325 zeta2 = 0.058 * log(e / (lamu + g2 * z13 * e)) - 0.14;
0326 zeta = zeta1 / zeta2;
0327 }
0328 else {
0329 zeta = 0.;
0330 }
0331 z2 = z * (z + zeta);
0332
0333
0334
0335
0336 screen0 = 2. * ame * sqrte * bbb / (z13 * ep);
0337 a0 = e * e1;
0338 a1 = ep * ep / a0;
0339 bet = 0.5 * a1;
0340 xi0 = 0.25 * rmass * rmass * a1;
0341 del = c8 / a0;
0342 tmn = G4Log((alf + 2. * del * a3) / (1. + (1. - del) * sqrt(a3)));
0343 sum = 0.;
0344 for (G4int i = 0; i < 8; ++i) {
0345 a4 = G4Exp(tmn * xgi[i]);
0346 a5 = a4 * (2. - a4);
0347 a6 = 1. - a5;
0348 a7 = 1. + a6;
0349 a9 = 3. + a6;
0350 xi = xi0 * a5;
0351 xii = 1. / xi;
0352 xi1 = 1. + xi;
0353 screen = screen0 * xi1 / a5;
0354 yeu = 5. - a6 + 4. * bet * a7;
0355 yed = 2. * (1. + 3. * bet) * G4Log(3. + xii) - a6 - a1 * (2. - a6);
0356 ye1 = 1. + yeu / yed;
0357 ale = G4Log(bbb / z13 * std::sqrt(xi1 * ye1) / (1. + screen * ye1));
0358 cre = 0.5 * G4Log(1. + (1.5 / rmass * z13) * (1.5 / rmass * z13) * xi1 * ye1);
0359 if (xi <= 1e3) {
0360 be = ((2. + a6) * (1. + bet) + xi * a9) * G4Log(1. + xii) + (a5 - bet) / xi1 - a9;
0361 }
0362 else {
0363 be = (3. - a6 + a1 * a7) / (2. * xi);
0364 }
0365 fe = std::max(0., (ale - cre) * be);
0366 ymu = 4. + a6 + 3. * bet * a7;
0367 ymd = a7 * (1.5 + a1) * G4Log(3. + xi) + 1. - 1.5 * a6;
0368 ym1 = 1. + ymu / ymd;
0369 alm_crm = G4Log(bbb * rmass / (1.5 * z13 * z13 * (1. + screen * ym1)));
0370 if (xi >= 1e-3) {
0371 a10 = (1. + a1) * a5;
0372 bm = (a7 * (1. + 1.5 * bet) - a10 * xii) * G4Log(xi1) + xi * (a5 - bet) / xi1 + a10;
0373 }
0374 else {
0375 bm = (5. - a6 + bet * a9) * (xi / 2.);
0376 }
0377 fm = max(0., (alm_crm)*bm);
0378
0379 sum = sum + a4 * (fe + fm / (rmass * rmass)) * wgi[i];
0380 }
0381 crp_g4 = -tmn * sum * (z2 / a) * coeff * e1 / (e * ep);
0382 return crp_g4;
0383 }
0384
0385
0386
0387 G4double MuCrossSections::CRM_Mephi(G4double Z, G4double tkin, G4double pairEnergy)
0388 {
0389
0390
0391
0392
0393
0394
0395
0396
0397 const G4double xgi[] = {0.0198550717512320, 0.1016667612931865, 0.2372337950418355,
0398 0.4082826787521750, 0.5917173212478250, 0.7627662049581645,
0399 0.8983332387068135, 0.9801449282487680};
0400
0401 const G4double wgi[] = {0.0506142681451880, 0.1111905172266870, 0.1568533229389435,
0402 0.1813418916891810, 0.1813418916891810, 0.1568533229389435,
0403 0.1111905172266870, 0.0506142681451880};
0404
0405 static const G4double factorForCross =
0406 2. / (3 * CLHEP::pi)
0407 * pow(CLHEP::fine_structure_const * CLHEP::classic_electr_radius / fMueRatio, 2);
0408
0409 if (pairEnergy <= 2. * fMuonMass) return 0.0;
0410
0411 G4double totalEnergy = tkin + fMuonMass;
0412 G4double residEnergy = totalEnergy - pairEnergy;
0413
0414 if (residEnergy <= fMuonMass) return 0.0;
0415
0416 G4double a0 = 1.0 / (totalEnergy * residEnergy);
0417 G4double rhomax = 1.0 - 2 * fMuonMass / pairEnergy;
0418 G4double tmnexp = 1. - rhomax;
0419
0420 if (tmnexp >= 1.0) {
0421 return 0.0;
0422 }
0423
0424 G4double tmn = G4Log(tmnexp);
0425
0426 G4double z2 = Z * Z;
0427 G4double beta = 0.5 * pairEnergy * pairEnergy * a0;
0428 G4double xi0 = 0.5 * beta;
0429
0430
0431 G4double rho[8];
0432 G4double rho2[8];
0433 G4double xi[8];
0434 G4double xi1[8];
0435 G4double xii[8];
0436
0437 for (G4int i = 0; i < 8; ++i) {
0438 rho[i] = G4Exp(tmn * xgi[i]) - 1.0;
0439 rho2[i] = rho[i] * rho[i];
0440 xi[i] = xi0 * (1.0 - rho2[i]);
0441 xi1[i] = 1.0 + xi[i];
0442 xii[i] = 1.0 / xi[i];
0443 }
0444
0445 G4double ximax = xi0 * (1. - rhomax * rhomax);
0446
0447 G4double Y = 10 * sqrt(fMuonMass / totalEnergy);
0448 G4double U[8];
0449
0450 for (G4int i = 0; i < 8; ++i) {
0451 U[i] = U_func(Z, rho2[i], xi[i], Y, pairEnergy);
0452 }
0453
0454 G4double UMax = U_func(Z, rhomax * rhomax, ximax, Y, pairEnergy);
0455
0456 G4double sum = 0.0;
0457
0458 for (G4int i = 0; i < 8; ++i) {
0459 G4double X = 1 + U[i] - UMax;
0460 G4double lnX = G4Log(X);
0461 G4double phi = ((2 + rho2[i]) * (1 + beta) + xi[i] * (3 + rho2[i])) * G4Log(1 + xii[i]) - 1
0462 - 3 * rho2[i] + beta * (1 - 2 * rho2[i])
0463 + ((1 + rho2[i]) * (1 + 1.5 * beta) - xii[i] * (1 + 2 * beta) * (1 - rho2[i]))
0464 * G4Log(xi1[i]);
0465 sum += wgi[i] * (1.0 + rho[i]) * phi * lnX;
0466 }
0467
0468 return -tmn * sum * factorForCross * z2 * residEnergy / (totalEnergy * pairEnergy);
0469 }
0470
0471
0472
0473 G4double MuCrossSections::U_func(G4double Z, G4double rho2, G4double xi, G4double Y,
0474 G4double pairEnergy, const G4double B)
0475 {
0476 G4int z = G4lrint(Z);
0477 G4double A27 = fNist->GetA27(z);
0478 G4double Z13 = fNist->GetZ13(z);
0479 static const G4double sqe = 2 * std::sqrt(G4Exp(1.0)) * fMuonMass * fMuonMass;
0480 G4double res = (0.65 * B / (A27 * Z13) * fMueRatio)
0481 / (1
0482 + (sqe * B * (1 + xi) * (1 + Y))
0483 / (Z13 * CLHEP::electron_mass_c2 * pairEnergy * (1 - rho2)));
0484 return res;
0485 }
0486
0487