source: trunk/examples/extended/electromagnetic/TestEm17/src/MuCrossSections.cc@ 1230

Last change on this file since 1230 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 12.5 KB
Line 
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: geant4-09-03-cand-01 $
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
38using namespace std;
39
40MuCrossSections::MuCrossSections()
41{ }
42
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44
45MuCrossSections::~MuCrossSections()
46{ }
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50G4double 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
70G4double 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
94double 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//***
166label10:
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
173double 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
215double 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
261double 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
Note: See TracBrowser for help on using the repository browser.