Jpp
JRadiation.hh
Go to the documentation of this file.
1 #ifndef __JPHYSICS__JRADIATION__
2 #define __JPHYSICS__JRADIATION__
3 
4 #include <cmath>
5 
6 #include <TRandom3.h>
7 
8 #include "JLang/JCC.hh"
9 #include "JTools/JConstants.hh"
10 
11 /**
12  * \file
13  * Muon radiative cross sections.
14  * \author P. Kooijman
15  */
16 
17 namespace JPHYSICS {}
18 namespace JPP { using namespace JPHYSICS; }
19 
20 namespace JPHYSICS {
21 
22  namespace {
23 
24  using JTOOLS::ALPHA_ELECTRO_MAGNETIC;
25  using JTOOLS::AVOGADRO;
26  using JTOOLS::MASS_ELECTRON;
27  using JTOOLS::MASS_MUON;
28  using JTOOLS::MASS_PROTON;
29  using JTOOLS::PI;
30 
31  static const double LAMBDA = 0.65; // [GeV]
32  }
33 
34 
35  /**
36  * Auxiliary class for the calculation of the muon radiative cross sections.
37  * \author P. Kooijman
38  */
39  class JRadiation {
40  public:
41 
42  /**
43  * Virtual desctructor.
44  */
45  virtual ~JRadiation()
46  {}
47 
48 
49  /**
50  * Constructor.
51  *
52  * \param z number of protons
53  * \param a number of nucleons
54  * \param integrsteps number of integration steps
55  * \param eminBrems energy threshold Bremsstrahlung [GeV]
56  * \param eminEErad energy threshold pair production [GeV]
57  * \param eminGNrad energy threshold photo-production [GeV]
58  */
59  JRadiation(const double z,
60  const double a,
61  const int integrsteps,
62  const double eminBrems,
63  const double eminEErad,
64  const double eminGNrad) :
65  Z(z),
66  A(a),
67  Dn (1.54*pow(A,0.27)),
68  DnP (pow(Dn,(1-1/Z))),
69  steps(integrsteps),
70  EminBrems(eminBrems),
71  EminEErad(eminEErad),
72  EminGNrad(eminGNrad)
73  {}
74 
75 
76  /**
77  * Pair production cross section.
78  *
79  * \param E muon energy [GeV]
80  * \param eps shower energy [GeV]
81  * \return cross section [m^2/g]
82  */
83  double SigmaEErad(const double E,
84  const double eps) const
85  {
86  //routine to calculate dsigma/de in cm^2/GeV for radiating a photon of energy eps
87  if(eps<4.*MASS_ELECTRON)return 0.;
88  if(eps>E-0.75*exp(1.0)*pow(Z,1./3.)) return 0.;
89  double zeta;
90  if(E<35.*MASS_MUON)
91  {
92  zeta = 0.;
93  }
94  else
95  {
96  zeta = (0.073*log((E/MASS_MUON)/(1.+1.95e-5*pow(Z,2./3.)*E/MASS_MUON))-0.26);
97  if(zeta<0.)
98  {
99  zeta=0.;
100  }
101  else
102  {
103  zeta = zeta/(0.058*log((E/MASS_MUON)/(1.+5.3e-5*pow(Z,1./3.)*E/MASS_MUON))-0.14);
104  }
105  }
106  double integ = IntegralofG(E,eps);
107  //cout<< eps<<" integ "<< integ<<endl;
108  return integ*(4/(3.*M_PI))*(Z*(Z+zeta)/A)*AVOGADRO*pow((ALPHA_ELECTRO_MAGNETIC*r0()),2.)*(1-eps/E)/eps;
109  }
110 
111 
112  /**
113  * Pair production cross section.
114  *
115  * \param E muon energy [GeV]
116  * \return cross section [m^2/g]
117  */
118  virtual double TotalCrossSectionEErad(const double E) const
119  {
120  //cout<<" entering TotcrsecEErad"<<endl;
121  // radiation goes approx as 1/E so sample in log E
122  // calculate the total cross section as S(dsigma/de)*de
123  // start at E and go backwards for precision
124  double sum(0.);
125  double dle = log(E/EminEErad)/steps;
126  //cout << "dle "<< dle<<endl;
127  for(int i=0;i<steps+1;i++)
128  {
129  double factor = 1.0;
130  if(i==0 || i==steps)factor=0.5;
131  //cout<<"i"<<i<<endl;
132  double eps = EminEErad*exp(i*dle);sum += factor*eps*SigmaEErad( E, eps);
133  //cout<<"eps "<<eps<<" sum "<<sum<<endl;
134  }
135  return sum*dle;
136  }
137 
138 
139  /**
140  * Bremsstrahlung cross section.
141  *
142  * \param E muon energy [GeV]
143  * \return cross section [m^2/g]
144  */
145  double TotalCrossSectionBrems(const double E) const
146  {
147  // double delta = (MASS_MUON*MASS_MUON*eps)/(2*E*(E-eps));
148  double delta = (MASS_MUON*MASS_MUON)/(2*E);
149  //double v = eps/E;
150  double Phin = log((B()*pow(Z,-1./3.)*(MASS_MUON+delta*(DnP*sqrt(exp(1.0))-2.)))/(DnP*(MASS_ELECTRON+delta*sqrt(exp(1.0))*B()*pow(Z,-1./3.))));
151  if(Phin<0.)Phin=0.;
152  double Phie = log((BP()*pow(Z,-2./3.)*MASS_MUON)/((1+delta*MASS_MUON/(MASS_ELECTRON*MASS_ELECTRON*sqrt(exp(1.0))))*(MASS_ELECTRON+delta*sqrt(exp(1.0))*BP()*pow(Z,-2./3.))));
153  if(Phie<0.)Phie=0.;
154  //if(eps>E/(1+MASS_MUON*MASS_MUON/(2.*MASS_ELECTRON*E)))Phie=0.;
155  double sig = (16./3.)*ALPHA_ELECTRO_MAGNETIC*AVOGADRO*pow((MASS_ELECTRON/MASS_MUON*r0()),2.0)*Z*(Z*Phin+Phie)/(A);
156  double epsint = log((E-MASS_MUON)/EminBrems)-(E-MASS_MUON-EminBrems)/E+0.375*(pow(E-MASS_MUON,2.)-pow(EminBrems,2.0))/pow(E,2.0);
157  return epsint*sig;//"cross section" in m^2/g multiplied by density and inverted gives mean free path
158  }
159 
160 
161  /**
162  * Photo-nuclear cross section.
163  *
164  * \param E muon energy [GeV]
165  * \param eps shower energy [GeV]
166  * \return cross section [m^2/g]
167  */
168  virtual double SigmaGNrad(const double E, const double eps) const
169  {
170  //minimum radiated energy is 0.2 GeV, not very critical formulae become invalid for very much lower
171  //result is given in m^2/g GeV
172  if (eps<0.2) return 0.;
173  if (eps>E-MASS_PROTON) return 0.;
174  double Aeff;
175  if (A==1)
176  {Aeff = 1.;}
177  else
178  {Aeff = (0.22*A + 0.78 * pow(A,0.89));}
179  double sigmaGammaPofeps = (49.2 + 11.1 * log(eps) + 151.8/sqrt(eps))*pow(10,-34.);
180  double Psiofeps = ALPHA_ELECTRO_MAGNETIC/PI*Aeff*AVOGADRO/A*sigmaGammaPofeps/eps;
181  double Denom = 1+eps/LAMBDA*(1+LAMBDA/(2*MASS_PROTON)+eps/LAMBDA);
182  double epsoverE = eps/E;
183  double Numerator = E * E * (1 - epsoverE) / (MASS_MUON * MASS_MUON) * (1 + MASS_MUON * MASS_MUON * epsoverE * epsoverE / (LAMBDA * LAMBDA * (1 - epsoverE)));
184  double Factor = 1 - epsoverE + epsoverE * epsoverE / 2 * (1 + 2 * MASS_MUON * MASS_MUON / (LAMBDA * LAMBDA));
185  double PhiofEofeps = epsoverE - 1 + Factor * log (Numerator / Denom);
186  //cout<<"PhiofEofeps "<<PhiofEofeps<<" Psiofeps "<<Psiofeps<<endl;
187  return Psiofeps*PhiofEofeps;
188  }
189 
190 
191  /**
192  * Photo-nuclear cross section.
193  *
194  * \param E muon energy [GeV]
195  * \return cross section [m^2/g]
196  */
197  virtual double TotalCrossSectionGNrad(const double E) const
198  {
199  double epsmin = 0.2;
200  double epsmax = E-MASS_PROTON/2;
201  double sum(0.);
202  double dle = log(epsmax/epsmin)/steps;
203  for(int i=0;i<steps+1;i++)
204  {
205  float factor = 1.0;
206  if(i==0 || i==steps)factor=0.5;
207  //cout<<"i"<<i<<endl;
208  double eps = epsmin*exp(i*dle);sum += factor*eps*SigmaGNrad( E, eps);
209  //cout<<"eps "<<eps<<" sum "<<sum<<endl;
210  }
211 
212  return sum*dle;
213  }
214 
215 
216  /**
217  * Bremsstrahlung shower energy.
218  *
219  * \param E muon energy [GeV]
220  * \return shower energy [GeV]
221  */
222  double EfromBrems(const double E) const
223  {
224  //double EminBrems(0.01);
225  //generate according to 1/k from minimum energy to E
226  double Er = 0.0;
227  for (int i = 1000; i != 0; --i) {
228  Er = EminBrems*exp(gRandom->Rndm()*log(E/EminBrems));
229  //check on the extra factor (1-v+3/4v^2)
230  if(gRandom->Rndm()<(1.-Er/E+0.75*pow(Er/E,2))) break;
231  }
232  return Er;
233  }
234 
235 
236  /**
237  * Pair production shower energy.
238  *
239  * \param E muon energy [GeV]
240  * \return shower energy [GeV]
241  */
242  double EfromEErad(const double E) const
243  {
244  //generate according to 1/k from minimum energy to E
245 
246  const double eps =0.2;
247  const double IntGmax = IntegralofG(E,eps)*2.0;
248 
249  double Er = 0.0;
250  for (int i = 1000; i != 0; --i)
251  {
252  Er = EminEErad*exp(gRandom->Rndm()*log(E/EminEErad));
253  //check on the extra factor, (1-v) and IntofG
254  double factor = (1.-Er/E)*IntegralofG(E,Er)/IntGmax;
255  if(gRandom->Rndm()<factor) break;
256  }
257  return Er;
258  }
259 
260  /**
261  * Photo-nuclear shower energy.
262  *
263  * \param E muon energy [GeV]
264  * \return shower energy [GeV]
265  */
266  double EfromGNrad(const double E) const
267  {
268  //generate according to 1/k from minimum energy to E
269 
270  double cmax = sigmaGammaPparam(EminGNrad);
271  if (cmax < sigmaGammaPparam(E))
272  cmax=sigmaGammaPparam(E);
273 
274  double Pmax = PhiofEofepsparam(E,EminGNrad);
275 
276  double Er = 0.0;
277  for (int i = 1000; i != 0; --i) {
278  Er = EminGNrad*exp(gRandom->Rndm()*log(E/EminGNrad));
279  const double factor = PhiofEofepsparam(E,Er)*sigmaGammaPparam(Er)/(cmax*Pmax);
280  if (gRandom->Rndm() < factor) return Er;
281  }
282 
283  return Er;
284  }
285 
286 
287  protected:
288  double GofZEvrho(const double E,
289  const double v,
290  const double r) const
291  {
292  const double ksi = MASS_MUON*MASS_MUON*v*v*(1-r*r)/(4*MASS_ELECTRON*MASS_ELECTRON*(1-v));//
293  const double b = v*v/(2*(1-v));//
294  const double Be = ((2+r*r)*(1+b)+ksi*(3+r*r))*log(1+1/ksi)+(1-r*r-b)/(1+ksi)-(3+r*r);
295  const double Bm = ((1+r*r)*(1+3*b/2)-(1+2*b)*(1-r*r)/ksi)*log(1+ksi)+ksi*(1-r*r-b)/(1+ksi)+(1+2*b)*(1-r*r);
296  const double Ye = (5-r*r+4*b*(1+r*r))/(2*(1+3*b)*log(3+1/ksi)-r*r-2*b*(2-r*r));
297  const double Ym = (4+r*r+3*b*(1+r*r))/((1+r*r)*(1.5+2*b)*log(3+ksi)+1-1.5*r*r);
298  const double Le = log((Astar()*pow(Z,-1./3.)*sqrt((1+ksi)*(1+Ye)))/
299  (1.+(2.*MASS_ELECTRON*sqrt(exp(1.))*Astar()*pow(Z,-1./3.)*(1+ksi)*(1+Ye))/(E*v*(1-r*r))));
300  const double Lm = log(((MASS_MUON/MASS_ELECTRON)*Astar()*pow(Z,-1./3.)*sqrt((1.+1./ksi)*(1.+Ym)))/
301  (1.+(2.*MASS_ELECTRON*sqrt(exp(1.))*Astar()*pow(Z,-1./3.)*(1+ksi)*(1+Ym))/(E*v*(1-r*r))));
302  double Phie = Be*Le; if(Phie<0.)Phie=0.;
303  double Phim = Bm*Lm; if(Phim<0.)Phim=0.;
304  return Phie+(MASS_ELECTRON*MASS_ELECTRON/(MASS_MUON*MASS_MUON))*Phim;
305  }
306 
307  virtual double IntegralofG(const double E,
308  const double eps) const
309  {
310  const double EP = E-eps;
311  const double v = eps/E;
312  const double tmin = log((4*MASS_ELECTRON/eps+12*pow(MASS_MUON,2)/(E*EP)*(1-4*MASS_ELECTRON/eps))
313  /(1+(1-6*pow(MASS_MUON,2)/(E*EP))*sqrt(1-4*MASS_ELECTRON/eps)));
314 
315  if (tmin < 0.0) {
316  const double dt = -tmin/steps;
317  double sum(0.);
318  for (int i = 0;i<steps+1;i++)
319  {
320  double fac = 1.0;if(i==0 || i==steps)fac=0.5;
321  double t = tmin+i*dt;
322  double r = 1.-exp(t);
323  sum+=fac*GofZEvrho(E,v,r)*exp(t);
324  }
325  return sum*dt;
326  } else
327  return 0.0;
328  }
329 
330  static double sigmaGammaPparam(const double eps)
331  {
332  return (49.2 + 11.1 * log(eps) + 151.8/sqrt(eps));
333  }
334 
335  static double PhiofEofepsparam(const double E,const double eps)
336  {
337  const double Denom = 1+eps/LAMBDA*(1+LAMBDA/(2*MASS_PROTON)+eps/LAMBDA);
338  const double epsoverE = eps/E;
339  const double Numerator = E * E * (1 - epsoverE) / (MASS_MUON * MASS_MUON) * (1 + MASS_MUON * MASS_MUON * epsoverE * epsoverE /
340  (LAMBDA * LAMBDA * (1 - epsoverE)));
341  const double Factor = 1 - epsoverE + epsoverE * epsoverE / 2 * (1 + 2 * MASS_MUON * MASS_MUON / (LAMBDA * LAMBDA));
342  return (epsoverE - 1 + Factor * log (Numerator / Denom));
343  }
344 
345 
346  static double r0 () { return 2.817940e-15; }
347  static double Astar() { return 183.0; }
348  static double B () { return 183.0; }
349  static double BP () { return 1429.0; }
350 
351  const double Z;
352  const double A;
353  const double Dn;
354  const double DnP;
355  const int steps;
356  const double EminBrems;
357  const double EminEErad;
358  const double EminGNrad;
359  };
360 }
361 
362 #endif
363 
JPHYSICS::JRadiation::r0
static double r0()
Definition: JRadiation.hh:346
JPHYSICS::JRadiation::steps
const int steps
Definition: JRadiation.hh:355
JPHYSICS::JRadiation::EminGNrad
const double EminGNrad
Definition: JRadiation.hh:358
JPHYSICS::JRadiation::PhiofEofepsparam
static double PhiofEofepsparam(const double E, const double eps)
Definition: JRadiation.hh:335
JTOOLS::MASS_MUON
static const double MASS_MUON
muon mass [GeV]
Definition: JConstants.hh:59
JPHYSICS::JRadiation::IntegralofG
virtual double IntegralofG(const double E, const double eps) const
Definition: JRadiation.hh:307
JTOOLS::ALPHA_ELECTRO_MAGNETIC
static const double ALPHA_ELECTRO_MAGNETIC
Electro-Magnetic coupling constant.
Definition: JConstants.hh:28
JTOOLS::MASS_PROTON
static const double MASS_PROTON
proton mass [GeV]
Definition: JConstants.hh:67
JPHYSICS::JRadiation::~JRadiation
virtual ~JRadiation()
Virtual desctructor.
Definition: JRadiation.hh:45
JPHYSICS::JRadiation::A
const double A
Definition: JRadiation.hh:352
JPHYSICS::JRadiation::BP
static double BP()
Definition: JRadiation.hh:349
JPHYSICS
Auxiliary classes and methods for calculation of PDF and muon energy loss.
Definition: JAbstractMedium.hh:9
JPHYSICS::JRadiation::Z
const double Z
Definition: JRadiation.hh:351
JTOOLS::AVOGADRO
static const double AVOGADRO
Avogadro's number [gr^-1].
Definition: JConstants.hh:24
JPHYSICS::JRadiation::JRadiation
JRadiation(const double z, const double a, const int integrsteps, const double eminBrems, const double eminEErad, const double eminGNrad)
Constructor.
Definition: JRadiation.hh:59
JPHYSICS::JRadiation::EfromBrems
double EfromBrems(const double E) const
Bremsstrahlung shower energy.
Definition: JRadiation.hh:222
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JConstants.hh
JPHYSICS::JRadiation::TotalCrossSectionBrems
double TotalCrossSectionBrems(const double E) const
Bremsstrahlung cross section.
Definition: JRadiation.hh:145
JPHYSICS::JRadiation::EminBrems
const double EminBrems
Definition: JRadiation.hh:356
JPHYSICS::JRadiation::Dn
const double Dn
Definition: JRadiation.hh:353
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:707
JTOOLS::MASS_ELECTRON
static const double MASS_ELECTRON
electron mass [GeV]
Definition: JConstants.hh:58
JPHYSICS::JRadiation::GofZEvrho
double GofZEvrho(const double E, const double v, const double r) const
Definition: JRadiation.hh:288
JPHYSICS::JRadiation::EfromEErad
double EfromEErad(const double E) const
Pair production shower energy.
Definition: JRadiation.hh:242
JPHYSICS::JRadiation::TotalCrossSectionEErad
virtual double TotalCrossSectionEErad(const double E) const
Pair production cross section.
Definition: JRadiation.hh:118
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
JPHYSICS::JRadiation::B
static double B()
Definition: JRadiation.hh:348
JPHYSICS::JRadiation::EfromGNrad
double EfromGNrad(const double E) const
Photo-nuclear shower energy.
Definition: JRadiation.hh:266
JTOOLS::r
data_type r[M+1]
Definition: JPolint.hh:709
JPHYSICS::JRadiation::sigmaGammaPparam
static double sigmaGammaPparam(const double eps)
Definition: JRadiation.hh:330
JPHYSICS::JRadiation::SigmaGNrad
virtual double SigmaGNrad(const double E, const double eps) const
Photo-nuclear cross section.
Definition: JRadiation.hh:168
JPHYSICS::JRadiation::DnP
const double DnP
Definition: JRadiation.hh:354
JPHYSICS::JRadiation::Astar
static double Astar()
Definition: JRadiation.hh:347
JPHYSICS::JRadiation::SigmaEErad
double SigmaEErad(const double E, const double eps) const
Pair production cross section.
Definition: JRadiation.hh:83
JPHYSICS::JRadiation
Auxiliary class for the calculation of the muon radiative cross sections.
Definition: JRadiation.hh:39
JPHYSICS::JRadiation::EminEErad
const double EminEErad
Definition: JRadiation.hh:357
JCC.hh
JPHYSICS::JRadiation::TotalCrossSectionGNrad
virtual double TotalCrossSectionGNrad(const double E) const
Photo-nuclear cross section.
Definition: JRadiation.hh:197
Generated by   doxygen 1.8.16