Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /// \file electromagnetic/TestEm17/src/MuCross << 26 // $Id: MuCrossSections.cc,v 1.2 2006/06/29 16:49:03 gunter Exp $ 27 /// \brief Implementation of the MuCrossSectio << 27 // GEANT4 tag $Name: geant4-08-03-patch-01 $ 28 // << 28 // 29 // << 30 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 31 33 #include "MuCrossSections.hh" 32 #include "MuCrossSections.hh" 34 33 35 #include "G4DataVector.hh" << 36 #include "G4Exp.hh" << 37 #include "G4Log.hh" << 38 #include "G4Material.hh" 34 #include "G4Material.hh" 39 #include "G4MuonMinus.hh" << 40 #include "G4NistManager.hh" << 41 #include "G4PhysicalConstants.hh" << 42 #include "G4SystemOfUnits.hh" << 43 35 44 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 37 46 using namespace std; 38 using namespace std; 47 << 39 48 MuCrossSections::MuCrossSections() 40 MuCrossSections::MuCrossSections() 49 { << 41 { } 50 fNist = G4NistManager::Instance(); << 42 51 fMuonMass = G4MuonMinus::MuonMinus()->GetPDG << 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 52 fMueRatio = fMuonMass / CLHEP::electron_mass << 44 53 } << 45 MuCrossSections::~MuCrossSections() >> 46 { } 54 47 55 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 49 57 G4double MuCrossSections::CR_Macroscopic(const << 50 G4double MuCrossSections::CR_Macroscopic(const G4String& process, G4Material* material, 58 G4dou 51 G4double tkin, G4double ep) >> 52 59 // return the macroscopic cross section (1/L) 53 // return the macroscopic cross section (1/L) in GEANT4 internal units 60 { 54 { 61 const G4ElementVector* theElementVector = ma 55 const G4ElementVector* theElementVector = material->GetElementVector(); 62 const G4double* NbOfAtomsPerVolume = materia 56 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume(); 63 57 64 G4double SIGMA = 0.; << 58 G4double SIGMA = 0 ; 65 G4int nelm = material->GetNumberOfElements() << 59 66 for (G4int i = 0; i < nelm; ++i) { << 60 for ( size_t i=0 ; i < material->GetNumberOfElements() ; i++ ) 67 const G4Element* element = (*theElementVec << 61 { >> 62 G4Element* element = (*theElementVector)[i]; 68 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAto 63 SIGMA += NbOfAtomsPerVolume[i] * CR_PerAtom(process, element, tkin, ep); 69 } 64 } 70 return SIGMA; 65 return SIGMA; 71 } 66 } 72 67 73 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 74 69 75 G4double MuCrossSections::CR_PerAtom(const G4S << 70 G4double MuCrossSections::CR_PerAtom(const G4String& process, G4Element* element, 76 G4double 71 G4double tkin, G4double ep) 77 { 72 { 78 G4double z = element->GetZ(); << 73 G4double z = element->GetZ(); 79 G4double a = element->GetA(); << 74 G4double a = element->GetA(); 80 G4double sigma = 0.; << 75 81 if (process == "muBrems") << 76 G4double sigma = 0.; 82 sigma = CRB_Mephi(z, a / (g / mole), tkin << 77 if (process == "muBrems") 83 << 78 sigma = CRB_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; 84 else if (process == "muIoni") << 79 85 sigma = CRK_Mephi(z, a / (g / mole), tkin << 80 else if (process == "muIoni") 86 << 81 sigma = CRK_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; 87 // else if (process == "muNucl") << 82 88 else if (process == "muonNuclear") << 83 else if (process == "muNucl") 89 sigma = CRN_Mephi(z, a / (g / mole), tkin << 84 sigma = CRN_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; 90 << 85 91 else if (process == "muPairProd") << 86 else if (process == "muPairProd") 92 sigma = CRP_Mephi(z, a / (g / mole), tkin << 87 sigma = CRP_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; 93 << 88 94 else if (process == "muToMuonPairProd") << 89 return sigma; 95 sigma = CRM_Mephi(z, tkin, ep); << 96 << 97 return sigma; << 98 } 90 } 99 91 100 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 101 93 102 G4double MuCrossSections::CRB_Mephi(G4double z << 94 double MuCrossSections::CRB_Mephi(double z,double a,double tkin,double ep) 103 95 104 //******************************************** 96 //*********************************************************************** 105 //*** crb_g4_1.inc in comparison << 97 //*** crb_g4_1.inc in comparison with crb_.inc, following 106 //*** changes are introduced (September << 98 //*** changes are introduced (September 24th, 1998): 107 //*** 1) function crb_g4 (Z,A,T << 99 //*** 1) function crb_g4 (Z,A,Tkin,EP), Tkin is kinetic energy 108 //*** 2) special case of hidrog << 100 //*** 2) special case of hidrogen (Z<1.5; b,b1,Dn_star) 109 //*** Numerical comparison: 5 d << 101 //*** Numerical comparison: 5 decimal digits coincide. 110 //*** << 102 //*** 111 //*** Cross section for bremsstrahlung << 103 //*** Cross section for bremsstrahlung by fast muon 112 //*** By R.P.Kokoulin, September 1998 << 104 //*** By R.P.Kokoulin, September 1998 113 //*** Formulae from Kelner,Kokoulin,Pet << 105 //*** Formulae from Kelner,Kokoulin,Petrukhin 1995, Preprint MEPhI 114 //*** (7,18,19,20,21,25,26); Dn (18) is << 106 //*** (7,18,19,20,21,25,26); Dn (18) is modified to incorporate 115 //*** Bugaev's inelatic nuclear correct << 107 //*** Bugaev's inelatic nuclear correction (28) for Z > 1. 116 //******************************************** 108 //*********************************************************************** 117 { 109 { 118 // G4double Z,A,Tkin,EP; << 110 // double Z,A,Tkin,EP; 119 G4double crb_g4; << 111 double crb_g4; 120 G4double e, v, delta, rab0, z_13, dn, b, b1, << 112 double e,v,delta,rab0,z_13,dn,b,b1,dn_star,rab1,fn,epmax1,fe,rab2; 121 // << 113 // 122 G4double ame = 0.51099907e-3; // GeV << 114 double ame=0.51099907e-3; // GeV 123 G4double lamu = 0.105658389; // GeV << 115 double amu=0.105658389; // GeV 124 G4double re = 2.81794092e-13; // cm << 116 double re=2.81794092e-13; // cm 125 G4double avno = 6.022137e23; << 117 double avno=6.022137e23; 126 G4double alpha = 1. / 137.036; << 118 double alpha=1./137.036; 127 G4double rmass = lamu / ame; // "207" << 119 double rmass=amu/ame; // "207" 128 G4double coeff = 16. / 3. * alpha * avno * ( << 120 double coeff=16./3.*alpha*avno*(re/rmass)*(re/rmass); // cm^2 129 G4double sqrte = 1.64872127; // sqrt(2.7182 << 121 double sqrte=1.64872127; // sqrt(2.71828...) 130 G4double btf = 183.; << 122 double btf=183.; 131 G4double btf1 = 1429.; << 123 double btf1=1429.; 132 G4double bh = 202.4; << 124 double bh=202.4; 133 G4double bh1 = 446.; << 125 double bh1=446.; 134 //*** << 126 //*** 135 if (ep >= tkin) { << 127 if(ep >= tkin) 136 return 0.; << 128 { 137 } << 129 crb_g4=0.; 138 e = tkin + lamu; << 130 return crb_g4; 139 v = ep / e; << 131 } 140 delta = lamu * lamu * v / (2. * (e - ep)); << 132 e=tkin+amu; 141 rab0 = delta * sqrte; << 133 v=ep/e; 142 G4int Z = G4lrint(z); << 134 delta=amu*amu*v/(2.*(e-ep)); // qmin 143 z_13 = 1. / fNist->GetZ13(z); // << 135 rab0=delta*sqrte; 144 //*** nuclear size and excitation, scr << 136 z_13=pow(z,-0.3333333); // 145 dn = 1.54 * fNist->GetA27(Z); << 137 //*** nuclear size and excitation, screening parameters 146 if (z <= 1.5) // special case for hydrogen << 138 dn=1.54*pow(a,0.27); 147 { << 139 if(z <= 1.5) // special case for hydrogen 148 b = bh; << 140 { 149 b1 = bh1; << 141 b=bh; 150 dn_star = dn; << 142 b1=bh1; 151 } << 143 dn_star=dn; 152 else { << 144 } 153 b = btf; << 145 else 154 b1 = btf1; << 146 { 155 dn_star = pow(dn, (1. - 1. / z)); // with << 147 b=btf; 156 } << 148 b1=btf1; 157 //*** nucleus contribution lo << 149 dn_star=pow(dn,(1.-1./z)); // with Bugaev's correction 158 rab1 = b * z_13; << 150 } 159 fn = G4Log(rab1 / (dn_star * (ame + rab0 * r << 151 //*** nucleus contribution logarithm 160 if (fn < 0.) fn = 0.; << 152 rab1=b*z_13; 161 //*** electron contribution l << 153 fn=log(rab1/(dn_star*(ame+rab0*rab1))*(amu+delta*(dn_star*sqrte-2.))); 162 epmax1 = e / (1. + lamu * rmass / (2. * e)); << 154 if(fn < 0.) fn=0.; 163 if (ep >= epmax1) { << 155 //*** electron contribution logarithm 164 fe = 0.; << 156 epmax1=e/(1.+amu*rmass/(2.*e)); 165 } << 157 if(ep >= epmax1) 166 else { << 158 { 167 rab2 = b1 * z_13 * z_13; << 159 fe=0.; 168 fe = G4Log(rab2 * lamu / ((1. + delta * rm << 160 goto label10; 169 if (fe < 0.) fe = 0.; << 161 } 170 } << 162 rab2=b1*z_13*z_13; 171 crb_g4 = coeff * (1. - v * (1. - 0.75 * v)) << 163 fe=log(rab2*amu/((1.+delta*rmass/(ame*sqrte))*(ame+rab0*rab2))); 172 return crb_g4; << 164 if(fe < 0.) fe=0.; >> 165 //*** >> 166 label10: >> 167 crb_g4=coeff*(1.-v*(1.-0.75*v))*z*(z*fn+fe)/(a*ep); >> 168 return crb_g4; 173 } 169 } 174 170 175 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 176 172 177 G4double MuCrossSections::CRK_Mephi(G4double z << 173 double MuCrossSections::CRK_Mephi(double z,double a,double tkin,double ep) >> 174 178 //******************************************** 175 //*********************************************************************** 179 //*** Cross section for knock-on electr << 176 //*** Cross section for knock-on electron production by fast muons 180 //*** (including bremsstrahlung e-diagr << 177 //*** (including bremsstrahlung e-diagrams and rad. correction). 181 //*** Units: cm^2/(g*GeV); Tkin, ep - G << 178 //*** Units: cm^2/(g*GeV); Tkin, ep - GeV. 182 //*** By R.P.Kokoulin, October 1998 << 179 //*** By R.P.Kokoulin, October 1998 183 //*** Formulae from Kelner,Kokoulin,Pet << 180 //*** Formulae from Kelner,Kokoulin,Petrukhin, Phys.Atom.Nuclei, 1997 184 //*** (a bit simplified Kelner's versio << 181 //*** (a bit simplified Kelner's version of Eq.30 - with 2 logarithms). 185 //*** 182 //*** 186 { 183 { 187 // G4double Z,A,Tkin,EP; << 184 // double Z,A,Tkin,EP; 188 G4double crk_g4; << 185 double crk_g4; 189 G4double v, sigma0, a1, a3; << 186 double e,epmax,v,sigma0,a1,a3; 190 // << 187 // 191 G4double ame = 0.51099907e-3; // GeV << 188 double ame=0.51099907e-3; // GeV 192 G4double lamu = 0.105658389; // GeV << 189 double amu=0.105658389; // GeV 193 G4double re = 2.81794092e-13; // cm << 190 double re=2.81794092e-13; // cm 194 G4double avno = 6.022137e23; << 191 double avno=6.022137e23; 195 G4double alpha = 1. / 137.036; << 192 double alpha=1./137.036; 196 G4double lpi = 3.141592654; << 193 double pi=3.141592654; 197 G4double bmu = lamu * lamu / (2. * ame); << 194 double bmu=amu*amu/(2.*ame); 198 G4double coeff0 = avno * 2. * lpi * ame * re << 195 double coeff0=avno*2.*pi*ame*re*re; 199 G4double coeff1 = alpha / (2. * lpi); << 196 double coeff1=alpha/(2.*pi); 200 //*** << 197 //*** 201 << 198 e=tkin+amu; 202 G4double e = tkin + lamu; << 199 epmax=e/(1.+bmu/e); 203 G4double epmax = e / (1. + bmu / e); << 200 if(ep >= epmax) 204 if (ep >= epmax) { << 201 { 205 return 0.0; << 202 crk_g4=0.; 206 } << 203 return crk_g4; 207 v = ep / e; << 204 } 208 sigma0 = coeff0 * (z / a) * (1. - ep / epmax << 205 v=ep/e; 209 a1 = G4Log(1. + 2. * ep / ame); << 206 sigma0=coeff0*(z/a)*(1.-ep/epmax+0.5*v*v)/(ep*ep); 210 a3 = G4Log(4. * e * (e - ep) / (lamu * lamu) << 207 a1=log(1.+2.*ep/ame); 211 crk_g4 = sigma0 * (1. + coeff1 * a1 * (a3 - << 208 a3=log(4.*e*(e-ep)/(amu*amu)); 212 return crk_g4; << 209 crk_g4=sigma0*(1.+coeff1*a1*(a3-a1)); >> 210 return crk_g4; 213 } 211 } 214 212 215 //....oooOO0OOooo........oooOO0OOooo........oo 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 216 214 217 G4double MuCrossSections::CRN_Mephi(G4double, << 215 double MuCrossSections::CRN_Mephi(double z,double a,double tkin,double ep) 218 216 219 //******************************************** 217 //*********************************************************************** 220 //*** Differential cross section for ph << 218 //*** Differential cross section for photonuclear muon interaction. 221 //*** Formulae from Borog & Petrukhin, << 219 //*** Formulae from Borog & Petrukhin, 1975 222 //*** Real photon cross section: Caldwe << 220 //*** Real photon cross section: Caldwell e.a., 1979 223 //*** Nuclear shadowing: Brodsky e.a., << 221 //*** Nuclear shadowing: Brodsky e.a., 1972 224 //*** Units: cm^2 / g GeV. << 222 //*** Units: cm^2 / g GeV. 225 //*** CRN_G4_1.inc January 31st, << 223 //*** CRN_G4_1.inc January 31st, 1998 R.P.Kokoulin 226 //******************************************** 224 //*********************************************************************** 227 { 225 { 228 // G4double Z,A,Tkin,EP; << 226 // double Z,A,Tkin,EP; 229 G4double crn_g4; << 227 double crn_g4; 230 G4double e, aeff, sigph, v, v1, v2, amu2, up << 228 double dummy,e,aeff,sigph,v,v1,v2,amu2,up,down; 231 //*** << 229 //*** 232 G4double lamu = 0.105658389; // GeV << 230 double amu=0.105658389; // GeV 233 G4double avno = 6.022137e23; << 231 double avno=6.022137e23; 234 G4double amp = 0.9382723; // GeV << 232 double amp=0.9382723; // GeV 235 G4double lpi = 3.14159265; << 233 double pi=3.14159265; 236 G4double alpha = 1. / 137.036; << 234 double alpha=1./137.036; 237 //*** << 235 //*** 238 G4double epmin_phn = 0.20; // GeV << 236 double epmin_phn=0.20; // GeV 239 G4double alam2 = 0.400000; // GeV**2 << 237 double alam2=0.400000; // GeV**2 240 G4double alam = 0.632456; // sqrt(alam2) << 238 double alam =0.632456; // sqrt(alam2) 241 G4double coeffn = alpha / lpi * avno * 1e-30 << 239 double coeffn=alpha/pi*avno*1e-30; // cm^2/microbarn 242 //*** << 240 //*** 243 e = tkin + lamu; << 241 dummy=z; // Z is a formal parameter at the moment 244 crn_g4 = 0.; << 242 e=tkin+amu; 245 if (ep >= e - 0.5 * amp || ep <= epmin_phn) << 243 crn_g4=0.; 246 aeff = 0.22 * a + 0.78 * pow(a, 0.89); // s << 244 if(ep >= e-0.5*amp) return crn_g4; 247 sigph = 49.2 + 11.1 * G4Log(ep) + 151.8 / st << 245 if(ep <= epmin_phn) return crn_g4; 248 v = ep / e; << 246 aeff=0.22*a+0.78*pow(a,0.89); // shadowing 249 v1 = 1. - v; << 247 sigph=49.2+11.1*log(ep)+151.8/sqrt(ep); // microbarn 250 v2 = v * v; << 248 v=ep/e; 251 amu2 = lamu * lamu; << 249 v1=1.-v; 252 up = e * e * v1 / amu2 * (1. + amu2 * v2 / ( << 250 v2=v*v; 253 down = 1. + ep / alam * (1. + alam / (2. * a << 251 amu2=amu*amu; 254 crn_g4 = coeffn * aeff / a * sigph / ep << 252 up=e*e*v1/amu2*(1.+amu2*v2/(alam2*v1)); 255 * (-v1 + (v1 + 0.5 * v2 * (1. + 2. << 253 down=1.+ep/alam*(1.+alam/(2.*amp)+ep/alam); 256 if (crn_g4 < 0.) crn_g4 = 0.; << 254 crn_g4=coeffn*aeff/a*sigph/ep*(-v1+(v1+0.5*v2*(1.+2.*amu2/alam2))*log(up/down)); 257 return crn_g4; << 255 if(crn_g4 < 0.) crn_g4=0.; >> 256 return crn_g4; 258 } 257 } 259 258 260 //....oooOO0OOooo........oooOO0OOooo........oo 259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 261 260 262 G4double MuCrossSections::CRP_Mephi(G4double z << 261 double MuCrossSections::CRP_Mephi(double z,double a,double tkin,double ep) 263 262 264 //******************************************** 263 //********************************************************************** 265 //*** crp_g4_1.inc in comparison << 264 //*** crp_g4_1.inc in comparison with crp_m.inc, following 266 //*** changes are introduced (January 1 << 265 //*** changes are introduced (January 16th, 1998): 267 //*** 1) Avno/A, cm^2/gram GeV << 266 //*** 1) Avno/A, cm^2/gram GeV 268 //*** 2) zeta_loss(E,Z) << 267 //*** 2) zeta_loss(E,Z) from Kelner 1997, Eqs.(53-54) 269 //*** 3) function crp_g4 (Z,A,T << 268 //*** 3) function crp_g4 (Z,A,Tkin,EP), Tkin is kinetic energy 270 //*** 4) bbb=183 (Thomas << 269 //*** 4) bbb=183 (Thomas-Fermi) 271 //*** 5) special case of hidrog << 270 //*** 5) special case of hidrogen (Z<1.5), bbb,g1,g2 272 //*** 6) expansions in 'xi' are << 271 //*** 6) expansions in 'xi' are simplified (Jan.17th,1998) 273 //*** << 272 //*** 274 //*** Cross section for electron pair p << 273 //*** Cross section for electron pair production by fast muon 275 //*** By R.P.Kokoulin, December 1997 << 274 //*** By R.P.Kokoulin, December 1997 276 //*** Formulae from Kokoulin & Petrukhi << 275 //*** 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, << 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 { 276 { 393 /* << 277 // double Z,A,Tkin,EP; 394 ||Cross section for pair production of m << 278 double crp_g4; 395 ||By Siddharth Yajaman 19-07-2022 << 279 double bbbtf,bbbh,g1tf,g2tf,g1h,g2h,e,z13,e1,alf,a3,bbb; 396 ||Based on the formulae from Kelner, Kok << 280 double g1,g2,zeta1,zeta2,zeta,z2,screen0,a0,a1,bet,xi0,del; 397 ||Physics of Atomic Nuclei, Vol. 63, No. << 281 double tmn,sum,a4,a5,a6,a7,a9,xi,xii,xi1,screen,yeu,yed,ye1; 398 ||Equations (15) - (22) << 282 double ale,cre,be,fe,ymu,ymd,ym1,alm_crm,a10,bm,fm; 399 */ << 283 // 400 << 284 double ame=0.51099907e-3; // GeV 401 const G4double xgi[] = {0.0198550717512320, << 285 double amu=0.105658389; // GeV 402 0.4082826787521750, << 286 double re=2.81794092e-13; // cm 403 0.8983332387068135, << 287 double avno=6.022137e23; 404 << 288 double pi=3.14159265; 405 const G4double wgi[] = {0.0506142681451880, << 289 double alpha=1./137.036; 406 0.1813418916891810, << 290 double rmass=amu/ame; // "207" 407 0.1111905172266870, << 291 double coeff=4./(3.*pi)*(alpha*re)*(alpha*re)*avno; // cm^2 408 << 292 double sqrte=1.64872127; // sqrt(2.71828...) 409 static const G4double factorForCross = << 293 double c3=3.*sqrte*amu/4.; // for limits 410 2. / (3 * CLHEP::pi) << 294 double c7=4.*ame; // -"- 411 * pow(CLHEP::fine_structure_const * CLHEP: << 295 double c8=6.*amu*amu; // -"- 412 << 296 413 if (pairEnergy <= 2. * fMuonMass) return 0.0 << 297 double xgi[8]={.0199,.1017,.2372,.4083,.5917,.7628,.8983,.9801}; // Gauss, 8 414 << 298 double wgi[8]={.0506,.1112,.1569,.1813,.1813,.1569,.1112,.0506}; // Gauss, 8 415 G4double totalEnergy = tkin + fMuonMass; << 299 bbbtf=183.; // for the moment... 416 G4double residEnergy = totalEnergy - pairEne << 300 bbbh=202.4; // for the moment... 417 << 301 g1tf=1.95e-5; 418 if (residEnergy <= fMuonMass) return 0.0; << 302 g2tf=5.3e-5; 419 << 303 g1h=4.4e-5; 420 G4double a0 = 1.0 / (totalEnergy * residEner << 304 g2h=4.8e-5; 421 G4double rhomax = 1.0 - 2 * fMuonMass / pair << 305 422 G4double tmnexp = 1. - rhomax; << 306 e=tkin+amu; 423 << 307 z13=pow(z,0.3333333); 424 if (tmnexp >= 1.0) { << 308 e1=e-ep; 425 return 0.0; << 309 crp_g4=0.; 426 } << 310 if(e1 <= c3*z13) return crp_g4; // ep > max 427 << 311 alf=c7/ep; // 4m/ep 428 G4double tmn = G4Log(tmnexp); << 312 a3=1.-alf; 429 << 313 if(a3 <= 0.) return crp_g4; // ep < min 430 G4double z2 = Z * Z; << 314 //*** zeta calculation 431 G4double beta = 0.5 * pairEnergy * pairEnerg << 315 if(z <= 1.5) // special case of hidrogen 432 G4double xi0 = 0.5 * beta; << 316 { 433 << 317 bbb=bbbh; 434 // Gaussian integration in ln(1-ro) ( with 8 << 318 g1=g1h; 435 G4double rho[8]; << 319 g2=g2h; 436 G4double rho2[8]; << 320 } 437 G4double xi[8]; << 321 else 438 G4double xi1[8]; << 322 { 439 G4double xii[8]; << 323 bbb=bbbtf; 440 << 324 g1=g1tf; 441 for (G4int i = 0; i < 8; ++i) { << 325 g2=g2tf; 442 rho[i] = G4Exp(tmn * xgi[i]) - 1.0; // rh << 326 } 443 rho2[i] = rho[i] * rho[i]; << 327 zeta1=0.073*log(e/(amu+g1*z13*z13*e))-0.26; 444 xi[i] = xi0 * (1.0 - rho2[i]); << 328 if(zeta1 > 0.) 445 xi1[i] = 1.0 + xi[i]; << 329 { 446 xii[i] = 1.0 / xi[i]; << 330 zeta2=0.058*log(e/(amu+g2*z13*e))-0.14; 447 } << 331 zeta=zeta1/zeta2; 448 << 332 } 449 G4double ximax = xi0 * (1. - rhomax * rhomax << 333 else 450 << 334 { 451 G4double Y = 10 * sqrt(fMuonMass / totalEner << 335 zeta=0.; 452 G4double U[8]; << 336 } 453 << 337 z2=z*(z+zeta); // 454 for (G4int i = 0; i < 8; ++i) { << 338 //*** just to check (for comparison with crp_m) 455 U[i] = U_func(Z, rho2[i], xi[i], Y, pairEn << 339 // z2=z*(z+1.) 456 } << 340 // bbb=189. 457 << 341 //*** 458 G4double UMax = U_func(Z, rhomax * rhomax, x << 342 screen0=2.*ame*sqrte*bbb/(z13*ep); // be careful with "ame" 459 << 343 a0=e*e1; 460 G4double sum = 0.0; << 344 a1=ep*ep/a0; // 2*beta 461 << 345 bet=0.5*a1; // beta 462 for (G4int i = 0; i < 8; ++i) { << 346 xi0=0.25*rmass*rmass*a1; // xi0 463 G4double X = 1 + U[i] - UMax; << 347 del=c8/a0; // 6mu^2/EE' 464 G4double lnX = G4Log(X); << 348 tmn=log((alf+2.*del*a3)/(1.+(1.-del)*sqrt(a3))); // log(1-rmax) 465 G4double phi = ((2 + rho2[i]) * (1 + beta) << 349 sum=0.; 466 - 3 * rho2[i] + beta * (1 - << 350 for(int i=0; i<=7; i++) // integration 467 + ((1 + rho2[i]) * (1 + 1.5 << 351 { 468 * G4Log(xi1[i]); << 352 a4=exp(tmn*xgi[i]); // 1-r 469 sum += wgi[i] * (1.0 + rho[i]) * phi * lnX << 353 a5=a4*(2.-a4); // 1-r2 470 } << 354 a6=1.-a5; // r2 471 << 355 a7=1.+a6; // 1+r2 472 return -tmn * sum * factorForCross * z2 * re << 356 a9=3.+a6; // 3+r2 >> 357 xi=xi0*a5; >> 358 xii=1./xi; >> 359 xi1=1.+xi; >> 360 screen=screen0*xi1/a5; >> 361 yeu=5.-a6+4.*bet*a7; >> 362 yed=2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6); >> 363 ye1=1.+yeu/yed; >> 364 ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)); >> 365 cre=0.5*log(1.+(1.5/rmass*z13)*(1.5/rmass*z13)*xi1*ye1); >> 366 if(xi <= 1e3) // >> 367 { >> 368 be=((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9; >> 369 } >> 370 else >> 371 { >> 372 be=(3.-a6+a1*a7)/(2.*xi); // -(6.-5.*a6+3.*bet*a6)/(6.*xi*xi); >> 373 } >> 374 fe=max(0.,(ale-cre)*be); >> 375 ymu=4.+a6+3.*bet*a7; >> 376 ymd=a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6; >> 377 ym1=1.+ymu/ymd; >> 378 alm_crm=log(bbb*rmass/(1.5*z13*z13*(1.+screen*ym1))); >> 379 if(xi >= 1e-3) // >> 380 { >> 381 a10=(1.+a1)*a5; // (1+2b)(1-r2) >> 382 bm=(a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10; >> 383 } >> 384 else >> 385 { >> 386 bm=(5.-a6+bet*a9)*(xi/2.); // -(11.-5.*a6+.5*bet*(5.+a6))*(xi*xi/6.) >> 387 } >> 388 fm=max(0.,(alm_crm)*bm); >> 389 //*** >> 390 sum=sum+a4*(fe+fm/(rmass*rmass))*wgi[i]; >> 391 } >> 392 crp_g4=-tmn*sum*(z2/a)*coeff*e1/(e*ep); >> 393 return crp_g4; 473 } 394 } 474 << 475 //....oooOO0OOooo........oooOO0OOooo........oo 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 476 396 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 397 491 //....oooOO0OOooo........oooOO0OOooo........oo << 492 398