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