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