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