Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-23 09:20:55

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file electromagnetic/TestEm17/src/MuCrossSections.cc
0027 /// \brief Implementation of the MuCrossSections class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056 
0057 G4double MuCrossSections::CR_Macroscopic(const G4String& process, const G4Material* material,
0058                                          G4double tkin, G4double ep)
0059 // return the macroscopic cross section (1/L) in GEANT4 internal units
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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   // else if (process == "muNucl")
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0101 
0102 G4double MuCrossSections::CRB_Mephi(G4double z, G4double a, G4double tkin, G4double ep)
0103 
0104 //***********************************************************************
0105 //***        crb_g4_1.inc        in comparison with crb_.inc, following
0106 //***        changes are introduced (September 24th, 1998):
0107 //***                1) function crb_g4 (Z,A,Tkin,EP), Tkin is kinetic energy
0108 //***                2) special case of hidrogen (Z<1.5; b,b1,Dn_star)
0109 //***                Numerical comparison: 5 decimal digits coincide.
0110 //***
0111 //***        Cross section for bremsstrahlung by fast muon
0112 //***        By R.P.Kokoulin, September 1998
0113 //***        Formulae from Kelner,Kokoulin,Petrukhin 1995, Preprint MEPhI
0114 //***        (7,18,19,20,21,25,26); Dn (18) is modified to incorporate
0115 //***        Bugaev's inelatic nuclear correction (28) for Z > 1.
0116 //***********************************************************************
0117 {
0118   //    G4double Z,A,Tkin,EP;
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;  // GeV
0123   G4double lamu = 0.105658389;  // GeV
0124   G4double re = 2.81794092e-13;  // cm
0125   G4double avno = 6.022137e23;
0126   G4double alpha = 1. / 137.036;
0127   G4double rmass = lamu / ame;  // "207"
0128   G4double coeff = 16. / 3. * alpha * avno * (re / rmass) * (re / rmass);  // cm^2
0129   G4double sqrte = 1.64872127;  // sqrt(2.71828...)
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));  // qmin
0141   rab0 = delta * sqrte;
0142   G4int Z = G4lrint(z);
0143   z_13 = 1. / fNist->GetZ13(z);  //
0144   //***       nuclear size and excitation, screening parameters
0145   dn = 1.54 * fNist->GetA27(Z);
0146   if (z <= 1.5)  // special case for hydrogen
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));  // with Bugaev's correction
0156   }
0157   //***                nucleus contribution logarithm
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   //***                electron contribution logarithm
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0176 
0177 G4double MuCrossSections::CRK_Mephi(G4double z, G4double a, G4double tkin, G4double ep)
0178 //***********************************************************************
0179 //***        Cross section for knock-on electron production by fast muons
0180 //***        (including bremsstrahlung e-diagrams and rad. correction).
0181 //***        Units: cm^2/(g*GeV); Tkin, ep - GeV.
0182 //***        By R.P.Kokoulin, October 1998
0183 //***        Formulae from Kelner,Kokoulin,Petrukhin, Phys.Atom.Nuclei, 1997
0184 //***        (a bit simplified Kelner's version of Eq.30 - with 2 logarithms).
0185 //***
0186 {
0187   //    G4double Z,A,Tkin,EP;
0188   G4double crk_g4;
0189   G4double v, sigma0, a1, a3;
0190   //
0191   G4double ame = 0.51099907e-3;  // GeV
0192   G4double lamu = 0.105658389;  // GeV
0193   G4double re = 2.81794092e-13;  // cm
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0216 
0217 G4double MuCrossSections::CRN_Mephi(G4double, G4double a, G4double tkin, G4double ep)
0218 
0219 //***********************************************************************
0220 //***        Differential cross section for photonuclear muon interaction.
0221 //***        Formulae from Borog & Petrukhin, 1975
0222 //***        Real photon cross section: Caldwell e.a., 1979
0223 //***        Nuclear shadowing: Brodsky e.a., 1972
0224 //***        Units: cm^2 / g GeV.
0225 //***        CRN_G4_1.inc        January 31st, 1998        R.P.Kokoulin
0226 //***********************************************************************
0227 {
0228   //    G4double Z,A,Tkin,EP;
0229   G4double crn_g4;
0230   G4double e, aeff, sigph, v, v1, v2, amu2, up, down;
0231   //***
0232   G4double lamu = 0.105658389;  // GeV
0233   G4double avno = 6.022137e23;
0234   G4double amp = 0.9382723;  // GeV
0235   G4double lpi = 3.14159265;
0236   G4double alpha = 1. / 137.036;
0237   //***
0238   G4double epmin_phn = 0.20;  // GeV
0239   G4double alam2 = 0.400000;  // GeV**2
0240   G4double alam = 0.632456;  // sqrt(alam2)
0241   G4double coeffn = alpha / lpi * avno * 1e-30;  // cm^2/microbarn
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);  // shadowing
0247   sigph = 49.2 + 11.1 * G4Log(ep) + 151.8 / std::sqrt(ep);  // microbarn
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0261 
0262 G4double MuCrossSections::CRP_Mephi(G4double z, G4double a, G4double tkin, G4double ep)
0263 
0264 //**********************************************************************
0265 //***        crp_g4_1.inc        in comparison with crp_m.inc, following
0266 //***        changes are introduced (January 16th, 1998):
0267 //***                1) Avno/A, cm^2/gram GeV
0268 //***                2) zeta_loss(E,Z)         from Kelner 1997, Eqs.(53-54)
0269 //***                3) function crp_g4 (Z,A,Tkin,EP), Tkin is kinetic energy
0270 //***                4) bbb=183        (Thomas-Fermi)
0271 //***                5) special case of hidrogen (Z<1.5), bbb,g1,g2
0272 //***                6) expansions in 'xi' are simplified        (Jan.17th,1998)
0273 //***
0274 //***        Cross section for electron pair production by fast muon
0275 //***        By R.P.Kokoulin, December 1997
0276 //***        Formulae from Kokoulin & Petrukhin 1971, Hobart, Eqs.(1,8,9,10)
0277 {
0278   //    G4double Z,A,Tkin,EP;
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;  // GeV
0286   G4double lamu = 0.105658389;  // GeV
0287   G4double re = 2.81794092e-13;  // cm
0288   G4double avno = 6.022137e23;
0289   G4double lpi = 3.14159265;
0290   G4double alpha = 1. / 137.036;
0291   G4double rmass = lamu / ame;  // "207"
0292   G4double coeff = 4. / (3. * lpi) * (alpha * re) * (alpha * re) * avno;  // cm^2
0293   G4double sqrte = 1.64872127;  // sqrt(2.71828...)
0294   G4double c3 = 3. * sqrte * lamu / 4.;  // for limits
0295   G4double c7 = 4. * ame;  // -"-
0296   G4double c8 = 6. * lamu * lamu;  // -"-
0297 
0298   G4double xgi[8] = {.0199, .1017, .2372, .4083, .5917, .7628, .8983, .9801};  // Gauss, 8
0299   G4double wgi[8] = {.0506, .1112, .1569, .1813, .1813, .1569, .1112, .0506};  // Gauss, 8
0300   bbbtf = 183.;  // for the moment...
0301   bbbh = 202.4;  // for the moment...
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;  // ep > max
0312   alf = c7 / ep;  // 4m/ep
0313   a3 = 1. - alf;
0314   if (a3 <= 0.) return crp_g4;  // ep < min
0315   //***                zeta calculation
0316   if (z <= 1.5)  // special case of hidrogen
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   //***                just to check (for comparison with crp_m)
0337   //        z2=z*(z+1.)
0338   //        bbb=189.
0339   //***
0340   screen0 = 2. * ame * sqrte * bbb / (z13 * ep);  // be careful with "ame"
0341   a0 = e * e1;
0342   a1 = ep * ep / a0;  // 2*beta
0343   bet = 0.5 * a1;  // beta
0344   xi0 = 0.25 * rmass * rmass * a1;  // xi0
0345   del = c8 / a0;  // 6mu^2/EE'
0346   tmn = G4Log((alf + 2. * del * a3) / (1. + (1. - del) * sqrt(a3)));  // log(1-rmax)
0347   sum = 0.;
0348   for (G4int i = 0; i < 8; ++i) {
0349     a4 = G4Exp(tmn * xgi[i]);  // 1-r
0350     a5 = a4 * (2. - a4);  // 1-r2
0351     a6 = 1. - a5;  // r2
0352     a7 = 1. + a6;  // 1+r2
0353     a9 = 3. + a6;  // 3+r2
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);  // -(6.-5.*a6+3.*bet*a6)/(6.*xi*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;  // (1+2b)(1-r2)
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.);  // -(11.-5.*a6+.5*bet*(5.+a6))*(xi*xi/6.)
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0390 
0391 G4double MuCrossSections::CRM_Mephi(G4double Z, G4double tkin, G4double pairEnergy)
0392 {
0393   /*
0394       ||Cross section for pair production of muons by fast muon
0395       ||By Siddharth Yajaman 19-07-2022
0396       ||Based on the formulae from Kelner, Kokoulin and Petrukhin,
0397       ||Physics of Atomic Nuclei, Vol. 63, No. 9, 2000, pp. 1603-1611
0398       ||Equations (15) - (22)
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   // Gaussian integration in ln(1-ro) ( with 8 points)
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;  // rho = -asymmetry
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......