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

Last change on this file since 1281 was 1230, checked in by garnier, 15 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.