| 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-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