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