Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
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
10
11
12/**
13 * \file
14 * Auxiliary methods for PDF calculations.
15 * \author mdejong, bjung
16 */
17
18namespace JPHYSICS {}
19namespace JPP { using namespace JPHYSICS; }
20
21namespace 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
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.
static const double ALPHA_ELECTRO_MAGNETIC
Electro-Magnetic coupling constant.
static const double DELTARAY_TMIN
Minimum allowed delta-ray kinetic energy [GeV].
static const double DELTARAY_TMAX
Maximum allowed delta-ray kinetic energy [GeV].
double getDeltaRayTmin()
Get minimum delta-ray kinetic energy.
double getDeltaRaysFromMuonFit(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
static const double DENSITY_SEA_WATER
Fixed environment values.
const double getRayleighScatteringLength(const double n, const double lambda)
Rayleigh scattering length.
double getDeltaRaysFromTauFit(const double E)
Equivalent EM-shower energy due to delta-rays per unit tau track length.
static const double MASS_ELECTRON
electron mass [GeV]
double getDeltaRayProbability(const double x)
Emission profile of photons from delta-rays.
const double getRayleighCrossSection(const double n, const double lambda)
Rayleigh cross section.
double cherenkov(const double lambda, const double n)
Number of Cherenkov photons per unit track length and per unit wavelength.
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
double getDeltaRayTmax(const double E, const double M)
Get maximum delta-ray kinetic energy for given lepton energy and mass.
static const double MASS_TAU
tau mass [GeV]
double getDeltaRays(const double E, const double M, const double Tmin, const double Tmax, const double Z, const double A)
Get equivalent EM-shower energy due to delta-rays per unit track length for an ionising particle with...
static const double AVOGADRO
Avogadro's number.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given energy and mass.
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
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.
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.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure to list files in directory.