Jpp  19.1.0-rc.1
the software that should make you happy
JPDFToolkit.hh
Go to the documentation of this file.
1 #ifndef __JPHYSICS__JPDFTOOLKIT__
2 #define __JPHYSICS__JPDFTOOLKIT__
3 
4 #include <cmath>
5 
6 #include "JTools/JRange.hh"
7 
8 #include "JPhysics/JConstants.hh"
10 
11 
12 /**
13  * \file
14  * Auxiliary methods for PDF calculations.
15  * \author mdejong, bjung
16  */
17 
18 namespace JPHYSICS {}
19 namespace JPP { using namespace JPHYSICS; }
20 
21 namespace JPHYSICS {
22 
23  using JTOOLS::JRange;
24 
25 
26  static const double DELTARAY_TMIN = 0.000915499; //!< Minimum allowed delta-ray kinetic energy [GeV]
27  static const double DELTARAY_TMAX = 1.0e10; //!< Maximum allowed delta-ray kinetic energy [GeV]
28 
29 
30  /**
31  * Get minimal wavelength for PDF evaluations.
32  *
33  * \return wavelength of light [nm]
34  */
35  inline double getMinimalWavelength()
36  {
37  return 300.0;
38  }
39 
40 
41  /**
42  * Get maximal wavelength for PDF evaluations.
43  *
44  * \return wavelength of light [nm]
45  */
46  inline double getMaximalWavelength()
47  {
48  return 700.0;
49  }
50 
51 
52  /**
53  * Number of Cherenkov photons per unit track length and per unit wavelength.
54  *
55  * \param lambda wavelength of light [nm]
56  * \param n index of refraction
57  * \return number of photons per unit track length and per unit wavelength [m^-1 nm^-1]
58  */
59  inline double cherenkov(const double lambda,
60  const double n)
61  {
62  const double x = n*lambda;
63 
64  return 1.0e9 * 2 * PI * ALPHA_ELECTRO_MAGNETIC * (n*n - 1.0) / (x*x);
65  }
66 
67 
68  /**
69  * Get minimum delta-ray kinetic energy.
70  *
71  * \return minimum delta-ray kinetic energy [GeV]
72  */
73  inline double getDeltaRayTmin()
74  {
75  const double Emin = MASS_ELECTRON / getSinThetaC(); // delta-ray Cherenkov threshold [GeV]
76 
77  return getKineticEnergy(Emin, MASS_ELECTRON);
78  }
79 
80 
81  /**
82  * Get maximum delta-ray kinetic energy for given lepton energy and mass.\n
83  * This formula is taken from reference
84  * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
85  * \n\n
86  * N.B.: This function should not be used for electrons.\n
87  * For electrons use `0.5 * getKineticEnergy(E, MASS_ELECTRON)` instead.
88  *
89  * \param E particle energy [GeV]
90  * \param M particle mass [GeV]
91  * \return maximum delta-ray kinetic energy [GeV]
92  */
93  inline double getDeltaRayTmax(const double E,
94  const double M)
95  {
96  const double ratio = MASS_ELECTRON / M;
97 
98  const double gamma = E / M;
99  const double beta = sqrt((1.0 + 1.0/gamma) * (1.0 - 1.0/gamma));
100 
101  return ( (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) /
102  (1.0 + 2.0*gamma*ratio + ratio*ratio) );
103  }
104 
105 
106  /**
107  * Get equivalent EM-shower energy due to delta-rays per unit track length\n
108  * for an ionising particle with given energy and given mass.\n\n
109  *
110  * N.B: For this function a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$ is assumed, where\n:
111  * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
112  * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
113  * - \f$E\f$ to the ionising particle's energy and
114  * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
115  *
116  * \param E particle energy [GeV]
117  * \param M particle mass [GeV]
118  * \param Tmin minimum delta-ray kinetic energy [GeV]
119  * \param Tmax maximum delta-ray kinetic energy [GeV]
120  * \param Z atomic number [unit]
121  * \param A atomic mass [g/mol]
122  * \return equivalent energy loss [GeV/m]
123  */
124  inline double getDeltaRays(const double E,
125  const double M,
126  const double Tmin,
127  const double Tmax,
128  const double Z,
129  const double A)
130  {
131  if (Tmax > Tmin) {
132 
133  const double K = 0.307075; // [MeV mol^-1 cm^2]
134 
135  const double gamma = E / M;
136  const double beta = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));
137 
138  const double a = 0.25/(E*E);
139  const double b = beta*beta/Tmax;
140  const double c = 1.0;
141 
142  const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]
143 
144  const double sT = Tmax + Tmin;
145  const double dT = Tmax - Tmin;
146  const double rT = Tmax / Tmin;
147 
148  const double weight = W * (a*sT*dT - b*dT + c*log(rT));
149 
150  return weight * DENSITY_SEA_WATER * 1.0e-1; // [GeV m^-1]
151 
152  } else {
153 
154  return 0.0;
155  }
156  }
157 
158 
159  /**
160  * Get equivalent EM-shower energy due to delta-rays per unit track length\n
161  * for an ionising particle with given energy and given mass and for a given form factor.\n
162  * The template parameter corresponds to a class which contains an `operator()(const double)`\n
163  * to compute the form factor corresponding to a given delta-ray kinetic energy.
164  *
165  * \param E particle energy [GeV]
166  * \param M particle mass [GeV]
167  * \param Tmin minimum delta-ray kinetic energy [GeV]
168  * \param Tmax maximum delta-ray kinetic energy [GeV]
169  * \param Z atomic number [unit]
170  * \param A atomic mass [g/mol]
171  * \param F form factor functor
172  * \param N number of points for numeric integration
173  * \return equivalent energy loss [GeV/m]
174  */
175  template<class JFormFactor_t>
176  inline double getDeltaRays(const double E,
177  const double M,
178  const double Tmin,
179  const double Tmax,
180  const double Z,
181  const double A,
182  const JFormFactor_t& F,
183  const int N = 1000000)
184  {
185  using namespace JPP;
186 
187  if (Tmax > Tmin) {
188 
189  const double K = 0.307075; // [MeV mol^-1 cm^2]
190 
191  const double gamma = E / M;
192  const double beta = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));
193 
194  const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]
195 
196  const double xmin = log(Tmin);
197  const double xmax = log(Tmax);
198  const double dx = (xmax - xmin) / ((double) N);
199 
200  double weight = 0.0;
201 
202  for (double x = xmin; x <= xmax; x += dx) {
203 
204  const double T = exp(x);
205  const double y = W * F(T) * dx; // [MeV g^-1 cm^-2]
206 
207  weight += y;
208  }
209 
210  return weight * DENSITY_SEA_WATER * 1.0e-1; // [GeV m^-1]
211 
212  } else {
213 
214  return 0.0;
215  }
216  }
217 
218 
219  /**
220  * Equivalent EM-shower energy due to delta-rays per unit electron track length.\n\n
221  *
222  * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
223  * and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
224  * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
225  * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
226  * - \f$E\f$ to the ionising particle's energy and
227  * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
228  *
229  * \param E electron energy [GeV]
230  * \param T_GeV allowed kinetic energy range of delta-rays [GeV]
231  * \return equivalent energy loss [GeV/m]
232  */
233  inline double getDeltaRaysFromElectron(const double E,
235  {
236  static const double Z = 10.0; // atomic number [unit]
237  static const double A = 18.0; // atomic mass [g/mol]
238 
239  const double Tmin = T_GeV.constrain(getDeltaRayTmin());
240  const double Tmax = T_GeV.constrain(0.5 * getKineticEnergy(E, MASS_ELECTRON));
241 
242  return getDeltaRays(E, MASS_ELECTRON, Tmin, Tmax, Z, A);
243  }
244 
245 
246  /**
247  * Equivalent EM-shower energy due to delta-rays per unit muon track length.\n\n
248  *
249  * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
250  * and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
251  * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
252  * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
253  * - \f$E\f$ to the ionising particle's energy and
254  * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
255  *
256  * \param E muon energy [GeV]
257  * \param T_GeV allowed kinetic energy range of delta-rays [GeV]
258  * \return equivalent energy loss [GeV/m]
259  */
260  inline double getDeltaRaysFromMuon(const double E,
262  {
263  static const double Z = 10.0; // atomic number [unit]
264  static const double A = 18.0; // atomic mass [g/mol]
265 
266  const double Tmin = T_GeV.constrain(getDeltaRayTmin());
267  const double Tmax = T_GeV.constrain(getDeltaRayTmax(E, MASS_MUON));
268 
269  return getDeltaRays(E, MASS_MUON, Tmin, Tmax, Z, A);
270  }
271 
272 
273  /**
274  * Equivalent EM-shower energy due to delta-rays per unit muon track length.\n\n
275  *
276  * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
277  * and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
278  * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
279  * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
280  * - \f$E\f$ to the ionising particle's energy and
281  * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
282  *
283  * \param E muon energy [GeV]
284  * \return equivalent energy loss [GeV/m]
285  */
286  inline double getDeltaRaysFromMuonFit(const double E)
287  {
288  static const double a = 3.195e-01;
289  static const double b = 3.373e-01;
290  static const double c = -2.731e-02;
291  static const double d = 1.610e-03;
292  static const double Emin = 0.13078; // [GeV]
293 
294  if (E > Emin) {
295 
296  const double x = log10(E); //
297  const double y = a + x*(b + x*(c + x*(d))); // [MeV g^-1 cm^2]
298 
299  if (y > 0.0) {
300  return y * DENSITY_SEA_WATER * 1.0e-1; // [GeV/m]
301  }
302  }
303 
304  return 0.0;
305  }
306 
307 
308  /**
309  * Equivalent EM-shower energy due to delta-rays per unit tau track length.\n\n
310  *
311  * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
312  * and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
313  * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
314  * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
315  * - \f$E\f$ to the ionising particle's energy and
316  * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
317  *
318  * \param E tau energy [GeV]
319  * \param T_GeV allowed kinetic energy range of delta-rays [GeV]
320  * \return equivalent energy loss [GeV/m]
321  */
322  inline double getDeltaRaysFromTau(const double E,
324  {
325  static const double Z = 10.0; // atomic number [unit]
326  static const double A = 18.0; // atomic mass [g/mol]
327 
328  const double Tmin = T_GeV.constrain(getDeltaRayTmin());
329  const double Tmax = T_GeV.constrain(getDeltaRayTmax(E, MASS_TAU));
330 
331  return getDeltaRays(E, MASS_TAU, Tmin, Tmax, Z, A);
332  }
333 
334 
335  /**
336  * Equivalent EM-shower energy due to delta-rays per unit tau track length.\n\n
337  *
338  * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
339  * and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
340  * - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
341  * - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
342  * - \f$E\f$ to the ionising particle's energy and
343  * - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
344  *
345  * \param E tau energy [GeV]
346  * \return equivalent energy loss [GeV/m]
347  */
348  inline double getDeltaRaysFromTauFit(const double E)
349  {
350  static const double a = -2.178e-01;
351  static const double b = 4.995e-01;
352  static const double c = -3.906e-02;
353  static const double d = 1.615e-03;
354  static const double Emin = 2.19500; // [GeV]
355 
356  if (E > Emin) {
357 
358  const double x = log10(E); //
359  const double y = a + x*(b + x*(c + x*(d))); // [MeV g^-1 cm^2]
360 
361  if (y > 0) {
362  return y * DENSITY_SEA_WATER * 1.0e-1; // [GeV/m]
363  }
364  }
365 
366  return 0.0;
367  }
368 
369 
370  /**
371  * Emission profile of photons from delta-rays.
372  *
373  * Profile is taken from reference ANTARES-SOFT-2002-015, J.\ Brunner (fig.\ 3).
374  *
375  * \param x cosine emission angle
376  * \return probability
377  */
378  inline double getDeltaRayProbability(const double x)
379  {
380  //return 1.0 / (4.0 * PI);
381  return 0.188 * exp(-1.25 * pow(fabs(x - 0.90), 1.30));
382  }
383 
384 
385  /**
386  * Rayleigh cross section.
387  *
388  * \param n index of refraction
389  * \param lambda wavelength of light [nm]
390  * \return cross section [cm^2]
391  */
392  inline const double getRayleighCrossSection(const double n,
393  const double lambda)
394  {
395  static const double d = 0.36; // size of H2O molecule [nm]
396  static const double U = PI*PI*PI*PI*PI*2.0/3.0;
397  static const double V = d*d*d*d*d*d;
398 
399  const double W = (n*n - 1.0) / (n*n + 2.0);
400  const double sigma = 1.0e-14 * U*V*W*W / (lambda*lambda*lambda*lambda); // [cm^2]
401 
402  return sigma;
403  }
404 
405 
406  /**
407  * Rayleigh scattering length.
408  *
409  * \param n index of refraction
410  * \param lambda wavelength of light [nm]
411  * \return scattering length [m]
412  */
413  inline const double getRayleighScatteringLength(const double n,
414  const double lambda)
415  {
416  static const double amu = 18.01528; // H20 mass in Atomic mass units
417 
418  const double sigma = getRayleighCrossSection(n, lambda);
419  const double ls = 1.0e-2 / (DENSITY_SEA_WATER * AVOGADRO * sigma / amu); // [m]
420 
421  return ls;
422  }
423 }
424 
425 #endif
Auxiliary methods for physics calculations.
Physics constants.
Auxiliary class to define a range between two values.
Range of values.
Definition: JRange.hh:42
static const uint32_t K[64]
Definition: crypt.cc:77
const double sigma[]
Definition: JQuadrature.cc:74
const double a
Definition: JQuadrature.cc:42
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
double getKineticEnergy(const Trk &trk)
Get track kinetic energy.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
static const double PI
Mathematical constants.
Auxiliary methods for light properties of deep-sea water.
double getDeltaRaysFromMuon(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:260
static const double ALPHA_ELECTRO_MAGNETIC
Electro-Magnetic coupling constant.
static const double DELTARAY_TMIN
Minimum allowed delta-ray kinetic energy [GeV].
Definition: JPDFToolkit.hh:26
static const double DELTARAY_TMAX
Maximum allowed delta-ray kinetic energy [GeV].
Definition: JPDFToolkit.hh:27
double getDeltaRayTmin()
Get minimum delta-ray kinetic energy.
Definition: JPDFToolkit.hh:73
double getDeltaRays(const double E, const double M, const double Tmin, const double Tmax, const double Z, const double A, const JFormFactor_t &F, const int N=1000000)
Get equivalent EM-shower energy due to delta-rays per unit track length for an ionising particle with...
Definition: JPDFToolkit.hh:176
double getDeltaRaysFromMuonFit(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:286
static const double DENSITY_SEA_WATER
Fixed environment values.
static const double MASS_MUON
muon mass [GeV]
const double getRayleighScatteringLength(const double n, const double lambda)
Rayleigh scattering length.
Definition: JPDFToolkit.hh:413
double getDeltaRaysFromTauFit(const double E)
Equivalent EM-shower energy due to delta-rays per unit tau track length.
Definition: JPDFToolkit.hh:348
static const double MASS_ELECTRON
electron mass [GeV]
double getDeltaRayProbability(const double x)
Emission profile of photons from delta-rays.
Definition: JPDFToolkit.hh:378
const double getRayleighCrossSection(const double n, const double lambda)
Rayleigh cross section.
Definition: JPDFToolkit.hh:392
double cherenkov(const double lambda, const double n)
Number of Cherenkov photons per unit track length and per unit wavelength.
Definition: JPDFToolkit.hh:59
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:35
double getDeltaRayTmax(const double E, const double M)
Get maximum delta-ray kinetic energy for given lepton energy and mass.
Definition: JPDFToolkit.hh:93
static const double MASS_TAU
tau mass [GeV]
static const double AVOGADRO
Avogadro's number.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:46
double getDeltaRaysFromTau(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit tau track length.
Definition: JPDFToolkit.hh:322
double getDeltaRaysFromElectron(const double E, const JRange< double > T_GeV=JRange< double >(DELTARAY_TMIN, DELTARAY_TMAX))
Equivalent EM-shower energy due to delta-rays per unit electron track length.
Definition: JPDFToolkit.hh:233
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
const int n
Definition: JPolint.hh:786
Auxiliary data structure to list files in directory.
Definition: JFilesystem.hh:20