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