1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: MuCrossSections.cc,v 1.2 2006/06/29 16:49:03 gunter Exp $ |
---|
27 | // GEANT4 tag $Name: $ |
---|
28 | // |
---|
29 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
30 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
31 | |
---|
32 | #include "MuCrossSections.hh" |
---|
33 | |
---|
34 | #include "G4Material.hh" |
---|
35 | |
---|
36 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
37 | |
---|
38 | using namespace std; |
---|
39 | |
---|
40 | MuCrossSections::MuCrossSections() |
---|
41 | { } |
---|
42 | |
---|
43 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
44 | |
---|
45 | MuCrossSections::~MuCrossSections() |
---|
46 | { } |
---|
47 | |
---|
48 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
49 | |
---|
50 | G4double MuCrossSections::CR_Macroscopic(const G4String& process, G4Material* material, |
---|
51 | G4double tkin, G4double ep) |
---|
52 | |
---|
53 | // return the macroscopic cross section (1/L) in GEANT4 internal units |
---|
54 | { |
---|
55 | const G4ElementVector* theElementVector = material->GetElementVector(); |
---|
56 | const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume(); |
---|
57 | |
---|
58 | G4double SIGMA = 0 ; |
---|
59 | |
---|
60 | for ( size_t i=0 ; i < material->GetNumberOfElements() ; i++ ) |
---|
61 | { |
---|
62 | G4Element* element = (*theElementVector)[i]; |
---|
63 | SIGMA += NbOfAtomsPerVolume[i] * CR_PerAtom(process, element, tkin, ep); |
---|
64 | } |
---|
65 | return SIGMA; |
---|
66 | } |
---|
67 | |
---|
68 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
69 | |
---|
70 | G4double MuCrossSections::CR_PerAtom(const G4String& process, G4Element* element, |
---|
71 | G4double tkin, G4double ep) |
---|
72 | { |
---|
73 | G4double z = element->GetZ(); |
---|
74 | G4double a = element->GetA(); |
---|
75 | |
---|
76 | G4double sigma = 0.; |
---|
77 | if (process == "muBrems") |
---|
78 | sigma = CRB_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; |
---|
79 | |
---|
80 | else if (process == "muIoni") |
---|
81 | sigma = CRK_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; |
---|
82 | |
---|
83 | else if (process == "muNucl") |
---|
84 | sigma = CRN_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; |
---|
85 | |
---|
86 | else if (process == "muPairProd") |
---|
87 | sigma = CRP_Mephi(z,a/(g/mole),tkin/GeV,ep/GeV)*(cm2/(g*GeV))*a/Avogadro; |
---|
88 | |
---|
89 | return sigma; |
---|
90 | } |
---|
91 | |
---|
92 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
93 | |
---|
94 | double MuCrossSections::CRB_Mephi(double z,double a,double tkin,double ep) |
---|
95 | |
---|
96 | //*********************************************************************** |
---|
97 | //*** crb_g4_1.inc in comparison with crb_.inc, following |
---|
98 | //*** changes are introduced (September 24th, 1998): |
---|
99 | //*** 1) function crb_g4 (Z,A,Tkin,EP), Tkin is kinetic energy |
---|
100 | //*** 2) special case of hidrogen (Z<1.5; b,b1,Dn_star) |
---|
101 | //*** Numerical comparison: 5 decimal digits coincide. |
---|
102 | //*** |
---|
103 | //*** Cross section for bremsstrahlung by fast muon |
---|
104 | //*** By R.P.Kokoulin, September 1998 |
---|
105 | //*** Formulae from Kelner,Kokoulin,Petrukhin 1995, Preprint MEPhI |
---|
106 | //*** (7,18,19,20,21,25,26); Dn (18) is modified to incorporate |
---|
107 | //*** Bugaev's inelatic nuclear correction (28) for Z > 1. |
---|
108 | //*********************************************************************** |
---|
109 | { |
---|
110 | // double Z,A,Tkin,EP; |
---|
111 | double crb_g4; |
---|
112 | double e,v,delta,rab0,z_13,dn,b,b1,dn_star,rab1,fn,epmax1,fe,rab2; |
---|
113 | // |
---|
114 | double ame=0.51099907e-3; // GeV |
---|
115 | double amu=0.105658389; // GeV |
---|
116 | double re=2.81794092e-13; // cm |
---|
117 | double avno=6.022137e23; |
---|
118 | double alpha=1./137.036; |
---|
119 | double rmass=amu/ame; // "207" |
---|
120 | double coeff=16./3.*alpha*avno*(re/rmass)*(re/rmass); // cm^2 |
---|
121 | double sqrte=1.64872127; // sqrt(2.71828...) |
---|
122 | double btf=183.; |
---|
123 | double btf1=1429.; |
---|
124 | double bh=202.4; |
---|
125 | double bh1=446.; |
---|
126 | //*** |
---|
127 | if(ep >= tkin) |
---|
128 | { |
---|
129 | crb_g4=0.; |
---|
130 | return crb_g4; |
---|
131 | } |
---|
132 | e=tkin+amu; |
---|
133 | v=ep/e; |
---|
134 | delta=amu*amu*v/(2.*(e-ep)); // qmin |
---|
135 | rab0=delta*sqrte; |
---|
136 | z_13=pow(z,-0.3333333); // |
---|
137 | //*** nuclear size and excitation, screening parameters |
---|
138 | dn=1.54*pow(a,0.27); |
---|
139 | if(z <= 1.5) // special case for hydrogen |
---|
140 | { |
---|
141 | b=bh; |
---|
142 | b1=bh1; |
---|
143 | dn_star=dn; |
---|
144 | } |
---|
145 | else |
---|
146 | { |
---|
147 | b=btf; |
---|
148 | b1=btf1; |
---|
149 | dn_star=pow(dn,(1.-1./z)); // with Bugaev's correction |
---|
150 | } |
---|
151 | //*** nucleus contribution logarithm |
---|
152 | rab1=b*z_13; |
---|
153 | fn=log(rab1/(dn_star*(ame+rab0*rab1))*(amu+delta*(dn_star*sqrte-2.))); |
---|
154 | if(fn < 0.) fn=0.; |
---|
155 | //*** electron contribution logarithm |
---|
156 | epmax1=e/(1.+amu*rmass/(2.*e)); |
---|
157 | if(ep >= epmax1) |
---|
158 | { |
---|
159 | fe=0.; |
---|
160 | goto label10; |
---|
161 | } |
---|
162 | rab2=b1*z_13*z_13; |
---|
163 | fe=log(rab2*amu/((1.+delta*rmass/(ame*sqrte))*(ame+rab0*rab2))); |
---|
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; |
---|
169 | } |
---|
170 | |
---|
171 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
172 | |
---|
173 | double MuCrossSections::CRK_Mephi(double z,double a,double tkin,double ep) |
---|
174 | |
---|
175 | //*********************************************************************** |
---|
176 | //*** Cross section for knock-on electron production by fast muons |
---|
177 | //*** (including bremsstrahlung e-diagrams and rad. correction). |
---|
178 | //*** Units: cm^2/(g*GeV); Tkin, ep - GeV. |
---|
179 | //*** By R.P.Kokoulin, October 1998 |
---|
180 | //*** Formulae from Kelner,Kokoulin,Petrukhin, Phys.Atom.Nuclei, 1997 |
---|
181 | //*** (a bit simplified Kelner's version of Eq.30 - with 2 logarithms). |
---|
182 | //*** |
---|
183 | { |
---|
184 | // double Z,A,Tkin,EP; |
---|
185 | double crk_g4; |
---|
186 | double e,epmax,v,sigma0,a1,a3; |
---|
187 | // |
---|
188 | double ame=0.51099907e-3; // GeV |
---|
189 | double amu=0.105658389; // GeV |
---|
190 | double re=2.81794092e-13; // cm |
---|
191 | double avno=6.022137e23; |
---|
192 | double alpha=1./137.036; |
---|
193 | double pi=3.141592654; |
---|
194 | double bmu=amu*amu/(2.*ame); |
---|
195 | double coeff0=avno*2.*pi*ame*re*re; |
---|
196 | double coeff1=alpha/(2.*pi); |
---|
197 | //*** |
---|
198 | e=tkin+amu; |
---|
199 | epmax=e/(1.+bmu/e); |
---|
200 | if(ep >= epmax) |
---|
201 | { |
---|
202 | crk_g4=0.; |
---|
203 | return crk_g4; |
---|
204 | } |
---|
205 | v=ep/e; |
---|
206 | sigma0=coeff0*(z/a)*(1.-ep/epmax+0.5*v*v)/(ep*ep); |
---|
207 | a1=log(1.+2.*ep/ame); |
---|
208 | a3=log(4.*e*(e-ep)/(amu*amu)); |
---|
209 | crk_g4=sigma0*(1.+coeff1*a1*(a3-a1)); |
---|
210 | return crk_g4; |
---|
211 | } |
---|
212 | |
---|
213 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
214 | |
---|
215 | double MuCrossSections::CRN_Mephi(double z,double a,double tkin,double ep) |
---|
216 | |
---|
217 | //*********************************************************************** |
---|
218 | //*** Differential cross section for photonuclear muon interaction. |
---|
219 | //*** Formulae from Borog & Petrukhin, 1975 |
---|
220 | //*** Real photon cross section: Caldwell e.a., 1979 |
---|
221 | //*** Nuclear shadowing: Brodsky e.a., 1972 |
---|
222 | //*** Units: cm^2 / g GeV. |
---|
223 | //*** CRN_G4_1.inc January 31st, 1998 R.P.Kokoulin |
---|
224 | //*********************************************************************** |
---|
225 | { |
---|
226 | // double Z,A,Tkin,EP; |
---|
227 | double crn_g4; |
---|
228 | double dummy,e,aeff,sigph,v,v1,v2,amu2,up,down; |
---|
229 | //*** |
---|
230 | double amu=0.105658389; // GeV |
---|
231 | double avno=6.022137e23; |
---|
232 | double amp=0.9382723; // GeV |
---|
233 | double pi=3.14159265; |
---|
234 | double alpha=1./137.036; |
---|
235 | //*** |
---|
236 | double epmin_phn=0.20; // GeV |
---|
237 | double alam2=0.400000; // GeV**2 |
---|
238 | double alam =0.632456; // sqrt(alam2) |
---|
239 | double coeffn=alpha/pi*avno*1e-30; // cm^2/microbarn |
---|
240 | //*** |
---|
241 | dummy=z; // Z is a formal parameter at the moment |
---|
242 | e=tkin+amu; |
---|
243 | crn_g4=0.; |
---|
244 | if(ep >= e-0.5*amp) return crn_g4; |
---|
245 | if(ep <= epmin_phn) return crn_g4; |
---|
246 | aeff=0.22*a+0.78*pow(a,0.89); // shadowing |
---|
247 | sigph=49.2+11.1*log(ep)+151.8/sqrt(ep); // microbarn |
---|
248 | v=ep/e; |
---|
249 | v1=1.-v; |
---|
250 | v2=v*v; |
---|
251 | amu2=amu*amu; |
---|
252 | up=e*e*v1/amu2*(1.+amu2*v2/(alam2*v1)); |
---|
253 | down=1.+ep/alam*(1.+alam/(2.*amp)+ep/alam); |
---|
254 | crn_g4=coeffn*aeff/a*sigph/ep*(-v1+(v1+0.5*v2*(1.+2.*amu2/alam2))*log(up/down)); |
---|
255 | if(crn_g4 < 0.) crn_g4=0.; |
---|
256 | return crn_g4; |
---|
257 | } |
---|
258 | |
---|
259 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
260 | |
---|
261 | double MuCrossSections::CRP_Mephi(double z,double a,double tkin,double ep) |
---|
262 | |
---|
263 | //********************************************************************** |
---|
264 | //*** crp_g4_1.inc in comparison with crp_m.inc, following |
---|
265 | //*** changes are introduced (January 16th, 1998): |
---|
266 | //*** 1) Avno/A, cm^2/gram GeV |
---|
267 | //*** 2) zeta_loss(E,Z) from Kelner 1997, Eqs.(53-54) |
---|
268 | //*** 3) function crp_g4 (Z,A,Tkin,EP), Tkin is kinetic energy |
---|
269 | //*** 4) bbb=183 (Thomas-Fermi) |
---|
270 | //*** 5) special case of hidrogen (Z<1.5), bbb,g1,g2 |
---|
271 | //*** 6) expansions in 'xi' are simplified (Jan.17th,1998) |
---|
272 | //*** |
---|
273 | //*** Cross section for electron pair production by fast muon |
---|
274 | //*** By R.P.Kokoulin, December 1997 |
---|
275 | //*** Formulae from Kokoulin & Petrukhin 1971, Hobart, Eqs.(1,8,9,10) |
---|
276 | { |
---|
277 | // double Z,A,Tkin,EP; |
---|
278 | double crp_g4; |
---|
279 | double bbbtf,bbbh,g1tf,g2tf,g1h,g2h,e,z13,e1,alf,a3,bbb; |
---|
280 | double g1,g2,zeta1,zeta2,zeta,z2,screen0,a0,a1,bet,xi0,del; |
---|
281 | double tmn,sum,a4,a5,a6,a7,a9,xi,xii,xi1,screen,yeu,yed,ye1; |
---|
282 | double ale,cre,be,fe,ymu,ymd,ym1,alm_crm,a10,bm,fm; |
---|
283 | // |
---|
284 | double ame=0.51099907e-3; // GeV |
---|
285 | double amu=0.105658389; // GeV |
---|
286 | double re=2.81794092e-13; // cm |
---|
287 | double avno=6.022137e23; |
---|
288 | double pi=3.14159265; |
---|
289 | double alpha=1./137.036; |
---|
290 | double rmass=amu/ame; // "207" |
---|
291 | double coeff=4./(3.*pi)*(alpha*re)*(alpha*re)*avno; // cm^2 |
---|
292 | double sqrte=1.64872127; // sqrt(2.71828...) |
---|
293 | double c3=3.*sqrte*amu/4.; // for limits |
---|
294 | double c7=4.*ame; // -"- |
---|
295 | double c8=6.*amu*amu; // -"- |
---|
296 | |
---|
297 | double xgi[8]={.0199,.1017,.2372,.4083,.5917,.7628,.8983,.9801}; // Gauss, 8 |
---|
298 | double wgi[8]={.0506,.1112,.1569,.1813,.1813,.1569,.1112,.0506}; // Gauss, 8 |
---|
299 | bbbtf=183.; // for the moment... |
---|
300 | bbbh=202.4; // for the moment... |
---|
301 | g1tf=1.95e-5; |
---|
302 | g2tf=5.3e-5; |
---|
303 | g1h=4.4e-5; |
---|
304 | g2h=4.8e-5; |
---|
305 | |
---|
306 | e=tkin+amu; |
---|
307 | z13=pow(z,0.3333333); |
---|
308 | e1=e-ep; |
---|
309 | crp_g4=0.; |
---|
310 | if(e1 <= c3*z13) return crp_g4; // ep > max |
---|
311 | alf=c7/ep; // 4m/ep |
---|
312 | a3=1.-alf; |
---|
313 | if(a3 <= 0.) return crp_g4; // ep < min |
---|
314 | //*** zeta calculation |
---|
315 | if(z <= 1.5) // special case of hidrogen |
---|
316 | { |
---|
317 | bbb=bbbh; |
---|
318 | g1=g1h; |
---|
319 | g2=g2h; |
---|
320 | } |
---|
321 | else |
---|
322 | { |
---|
323 | bbb=bbbtf; |
---|
324 | g1=g1tf; |
---|
325 | g2=g2tf; |
---|
326 | } |
---|
327 | zeta1=0.073*log(e/(amu+g1*z13*z13*e))-0.26; |
---|
328 | if(zeta1 > 0.) |
---|
329 | { |
---|
330 | zeta2=0.058*log(e/(amu+g2*z13*e))-0.14; |
---|
331 | zeta=zeta1/zeta2; |
---|
332 | } |
---|
333 | else |
---|
334 | { |
---|
335 | zeta=0.; |
---|
336 | } |
---|
337 | z2=z*(z+zeta); // |
---|
338 | //*** just to check (for comparison with crp_m) |
---|
339 | // z2=z*(z+1.) |
---|
340 | // bbb=189. |
---|
341 | //*** |
---|
342 | screen0=2.*ame*sqrte*bbb/(z13*ep); // be careful with "ame" |
---|
343 | a0=e*e1; |
---|
344 | a1=ep*ep/a0; // 2*beta |
---|
345 | bet=0.5*a1; // beta |
---|
346 | xi0=0.25*rmass*rmass*a1; // xi0 |
---|
347 | del=c8/a0; // 6mu^2/EE' |
---|
348 | tmn=log((alf+2.*del*a3)/(1.+(1.-del)*sqrt(a3))); // log(1-rmax) |
---|
349 | sum=0.; |
---|
350 | for(int i=0; i<=7; i++) // integration |
---|
351 | { |
---|
352 | a4=exp(tmn*xgi[i]); // 1-r |
---|
353 | a5=a4*(2.-a4); // 1-r2 |
---|
354 | a6=1.-a5; // r2 |
---|
355 | a7=1.+a6; // 1+r2 |
---|
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; |
---|
394 | } |
---|
395 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
396 | |
---|
397 | |
---|