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