Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JPDF.hh
Go to the documentation of this file.
1#ifndef __JPHYSICS__JPDF__
2#define __JPHYSICS__JPDF__
3
4#include <vector>
5#include <algorithm>
6#include <cmath>
7#include <limits>
8
9#include "JLang/JCC.hh"
10#include "JLang/JException.hh"
11#include "JPhysics/JConstants.hh"
12#include "JTools/JQuadrature.hh"
13#include "JPhysics/JDispersionInterface.hh"
14#include "JPhysics/JDispersion.hh"
15#include "JPhysics/JAbstractPMT.hh"
16#include "JPhysics/JAbstractMedium.hh"
17#include "JPhysics/JPDFToolkit.hh"
18#include "JPhysics/JPDFTypes.hh"
19#include "JPhysics/JGeane.hh"
20#include "JPhysics/JGeant.hh"
21#include "JPhysics/JGeanz.hh"
22
23
24/**
25 * \author mdejong
26 */
27
28namespace JPHYSICS {}
29namespace JPP { using namespace JPHYSICS; }
30
31namespace JPHYSICS {
32
33 using JLANG::JException;
34
35 /**
36 * Radius of optical module [m].
37 *
38 * This parameter is used to implement shadowing of the PMT by the optical module.
39 */
40 static double MODULE_RADIUS_M = 0.216;
41
42
43 /**
44 * Low level interface for the calculation of the Probability Density Functions (PDFs) of
45 * the arrival time of Cherenkov light from a muon or an EM-shower on a photo-multiplier tube (PMT).
46 * The calculation of the PDF covers direct light light and single-scattered light.
47 * The direct light is attenuated by absorption and any form of light scattering.\n
48 * The attenuation length is defined as:
49 \f[
50 \frac{1}{\lambda_{att}} ~\equiv~ \frac{1}{\lambda_{abs}} + \frac{1}{\lambda_{s}}
51 \f]
52 * where
53 \f$ \lambda_{abs} \f$ and
54 \f$ \lambda_{s} \f$
55 refer to the absorption and scattering length, respectively.\n
56 * The scattered light is attenuated by absorption and light scattering weighed with the scattering angle \f$ \theta_{s} \f$.
57 * The attenuation length is in this case defined as:
58 \f[
59 \frac{1}{\lambda_{att}} ~\equiv~ \frac{1}{\lambda_{abs}} + \frac{1-\left<\cos\theta_s\right>}{\lambda_{s}}
60 \f]
61 *
62 * The PDFs of a muon are evaluated as a function of:
63 * - distance of closest approach of the muon to the PMT, and
64 * - orientation of the PMT (zenith and azimuth angle).
65 *
66 * The PDFs of an EM-shower are evaluated as a function of:
67 * - distance between the EM-shower and the PMT,
68 * - cosine angle EM-shower direction and EM-shower - PMT position, and
69 * - orientation of the PMT (zenith and azimuth angle).
70 *
71 * The intensity of light from an EM-shower is assumed to be proportional to the energy of the shower.\n
72 * For a given energy of the EM-shower, the longitudinal profile of the EM-shower is taken into account.
73 *
74 * The coordinate system is defined such that z-axis corresponds to the muon trajectory (EM-shower).\n
75 * The orientation of the PMT (i.e.\ zenith and azimuth angle) are defined following a rotation around the z-axis,
76 * such that the PMT is located in the x-z plane.
77 *
78 * Schematics of the PDF of direct light from a muon:
79 *
80 \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)
81
82 \put( 1.0, 0.5){\vector(0,1){9}}
83 \put( 1.0,10.0){\makebox(0,0){$z$}}
84 \put( 1.0, 0.0){\makebox(0,0){muon}}
85
86 \put( 1.0, 8.0){\line(1,0){6}}
87 \put( 4.0, 8.5){\makebox(0,0)[c]{$R$}}
88
89 \multiput( 1.0, 2.0)(0.5, 0.5){12}{\qbezier(0.0,0.0)(-0.1,0.35)(0.25,0.25)\qbezier(0.25,0.25)(0.6,0.15)(0.5,0.5)}
90 \put( 4.5, 4.5){\makebox(0,0)[l]{photon}}
91
92 \put( 1.0, 2.0){\circle*{0.2}}
93 \put( 0.5, 2.0){\makebox(0,0)[r]{$(0,0,z,t_{0})$}}
94
95 \put( 1.0, 8.0){\circle*{0.2}}
96 \put( 0.5, 8.0){\makebox(0,0)[r]{$(0,0,0)$}}
97
98 \put( 7.0, 8.0){\circle*{0.2}}
99 \put( 7.0, 9.0){\makebox(0,0)[c]{PMT}}
100 \put( 7.5, 8.0){\makebox(0,0)[l]{$(R,0,0,t_{1})$}}
101
102 \put( 1.1, 2.1){
103 \put(0.0,1.00){\vector(-1,0){0.1}}
104 \qbezier(0.0,1.0)(0.25,1.0)(0.5,0.75)
105 \put(0.5,0.75){\vector(1,-1){0.1}}
106 \put(0.4,1.5){\makebox(0,0){$\theta_{c}$}}
107 }
108
109 \end{picture}
110 \f}
111 *
112 * For the PDF of a muon, the time is defined such that \f$ t ~=~ 0 \f$
113 * refers to the time when the muon passes through the point at minimal approach to the PMT
114 * (i.e.\ position \f$ (0,0,0) \f$ in the above figure).
115 *
116 * The delay for the Cherenkov cone to actually pass through the PMT is defined as
117 \f[
118 \Delta t ~\equiv~ R \tan(\theta_{c}) / c
119 \f]
120 * where \f$ \tan(\theta_{c}) \f$ is given by method getTanThetaC.
121 *
122 * Given a measured hit time, \f$ t_{hit} \f$, the PDF should be probed at \f$ (t_{hit} - t_1) \f$,
123 * where \f$ t_1 \f$ is given by:
124 *
125 \f[
126 ct_1 = ct_0 - z + R \tan(\theta_{c})
127 \f]
128 *
129 * In this, \f$ t_0 \f$ refers to the time the muon passes through point \f$ (0,0,z) \f$.
130 *
131 * Schematics of single scattered light from a muon:
132 *
133 \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)
134
135 \put( 1.0, 0.5){\vector(0,1){9}}
136 \put( 1.0,10.0){\makebox(0,0){$z$}}
137 \put( 1.0, 0.0){\makebox(0,0){muon}}
138
139 \put( 1.0, 8.0){\line(1,0){2}}
140 \put( 2.0, 8.5){\makebox(0,0)[c]{$R$}}
141
142 \multiput( 1.0, 2.0)( 0.5, 0.5){8}{\qbezier(0.0,0.0)(-0.1,0.35)( 0.25,0.25)\qbezier( 0.25,0.25)( 0.6,0.15)( 0.5,0.5)}
143 \multiput( 5.0, 6.0)(-0.5, 0.5){4}{\qbezier(0.0,0.0)(+0.1,0.35)(-0.25,0.25)\qbezier(-0.25,0.25)(-0.6,0.15)(-0.5,0.5)}
144
145 \put( 5.0, 6.0){\circle*{0.2}}
146
147 \put( 3.5, 3.5){\makebox(0,0)[c]{$\bar{u}$}}
148 \put( 5.0, 7.0){\makebox(0,0)[c]{$\bar{v}$}}
149
150 \put( 1.0, 2.0){\circle*{0.2}}
151 \put( 0.5, 2.0){\makebox(0,0)[r]{$(0,0,z,t_0)$}}
152
153 \put( 1.0, 8.0){\circle*{0.2}}
154 \put( 0.5, 8.0){\makebox(0,0)[r]{$(0,0,0)$}}
155
156 \put( 3.0, 8.0){\circle*{0.2}}
157 \put( 3.0, 9.0){\makebox(0,0)[c]{PMT}}
158 \put( 3.5, 8.0){\makebox(0,0)[l]{$(R,0,0,t_1)$}}
159
160 \end{picture}
161 \f}
162 *
163 * In this,
164 *
165 \f{math}{
166 \bar{u} ~\equiv~ u ~ \hat{u}
167 ~\equiv~ u \left(\begin{array}{c} \sin(\theta_{0}) \cos(\phi_{0}) \\ \sin(\theta_{0}) \sin(\phi_{0}) \\ \cos(\theta_{0})\end{array}\right)
168 \textrm{ and }
169 \bar{v} ~\equiv~ v ~ \hat{v}
170 ~\equiv~ v \left(\begin{array}{c} \sin(\theta_{1}) \cos(\phi_{1}) \\ \sin(\theta_{1}) \sin(\phi_{1}) \\ \cos(\theta_{1})\end{array}\right)
171 \f}
172 *
173 * The vectors \f$ \bar{u} \f$ and \f$ \bar{v} \f$ are subject to the following constraint:
174 \f{math}{
175 \begin{array}{ccccccc}
176
177 \begin{array}{c}
178 \mathrm{start}\\
179 \mathrm{point}
180 \end{array}
181
182 &&
183
184 \begin{array}{c}
185 \mathrm{before}\\
186 \mathrm{scattering}
187 \end{array}
188
189 &&
190
191 \begin{array}{c}
192 \mathrm{after}\\
193 \mathrm{scattering}
194 \end{array}
195
196 &&
197
198 \mathrm{PMT}
199
200 \\
201 \\
202
203 \left(\begin{array}{c} 0 \\ 0 \\ z \end{array}\right)
204 & + &
205 u \left(\begin{array}{c} \sin(\theta_{0}) \cos(\phi_{0}) \\ \sin(\theta_{0}) \sin(\phi_{0}) \\ \cos(\theta_{0})\end{array}\right)
206 & + &
207 v \left(\begin{array}{c} \sin(\theta_{1}) \cos(\phi_{1}) \\ \sin(\theta_{1}) \sin(\phi_{1}) \\ \cos(\theta_{1})\end{array}\right)
208 & = &
209 \left(\begin{array}{c} R \\ 0 \\ 0 \end{array}\right)
210 \end{array}
211 * \f}
212 *
213 * Note that an expression for the cosine of the scattering angle, \f$ \cos(\theta_s) ~\equiv~ \hat{u} \cdot \hat{v} \f$,
214 * can directly be obtained by multiplying the left hand side and the right hand side with \f$ \hat{u} \f$:
215 *
216 \f[
217 \cos\theta_s = \frac{R\sin\theta_0\cos\phi_0 - z\cos\theta_0 - u}{v}
218 \f]
219 *
220 *
221 * The arrival time of the single scattered light can be expressed as:
222 *
223 \f[
224 ct_1 = ct_0 + n (u + v)
225 \f]
226 *
227 * where \f$ n \f$ is the index of refraction.
228 *
229 * Schematics of direct light from a EM-shower:
230 *
231 \f{center}\setlength{\unitlength}{0.7cm}\begin{picture}(8,12)
232
233 \put( 1.0, 2.0){\vector(0,1){1}}
234 \put( 1.0, 3.5){\multiput(0,0)(0,0.5){12}{\line(0,1){0.2}}}
235 \put( 1.0,10.0){\makebox(0,0){$z$}}
236 \put( 1.0, 1.0){\makebox(0,0){EM-shower}}
237
238 \put( 1.0, 8.0){\line(1,0){6}}
239 \put( 4.0, 8.5){\makebox(0,0)[c]{$R$}}
240
241 \multiput( 1.0, 2.0)(0.5, 0.5){12}{\qbezier(0.0,0.0)(-0.1,0.35)(0.25,0.25)\qbezier(0.25,0.25)(0.6,0.15)(0.5,0.5)}
242 \put( 4.5, 4.5){\makebox(0,0)[l]{photon}}
243
244 \put( 1.0, 2.0){\circle*{0.2}}
245 \put( 0.5, 2.0){\makebox(0,0)[r]{$(0,0,z,t_0)$}}
246
247 \put( 1.0, 8.0){\circle*{0.2}}
248 \put( 0.5, 8.0){\makebox(0,0)[r]{$(0,0,0)$}}
249
250 \put( 7.0, 8.0){\circle*{0.2}}
251 \put( 7.0, 9.0){\makebox(0,0)[c]{PMT}}
252 \put( 7.5, 8.0){\makebox(0,0)[l]{$(R,0,0,t_1)$}}
253
254 \put( 1.1, 2.1){
255 \put(0.0,1.00){\vector(-1,0){0.1}}
256 \qbezier(0.0,1.0)(0.25,1.0)(0.5,0.75)
257 \put(0.5,0.75){\vector(1,-1){0.1}}
258 \put(0.4,1.5){\makebox(0,0){$\theta_0$}}
259 }
260
261 \end{picture}
262 \f}
263 *
264 * Note that the emission angle \f$ \theta_0 \f$ is in general not equal to the Cherenkov angle \f$ \theta_c \f$.\n
265 * For an EM-shower, the arrival time of the light can be expressed as:
266 *
267 \f[
268 ct_1 = ct_0 + n D
269 \f]
270 *
271 * where \f$ D \f$ is the total distance traveled by the photon from its emission point to the PMT,
272 * including possible scattering.
273 *
274 * For the computation of the various integrals, it is assumed that the index of refraction
275 * decreases monotonously as a function of the wavelength of the light.
276 */
277 class JPDF :
278 public JTOOLS::JGaussLegendre,
279 public virtual JDispersionInterface,
280 public virtual JAbstractPMT,
281 public virtual JAbstractMedium
282 {
283
284 typedef JTOOLS::JElement2D<double, double> element_type;
285
286 public:
287 /**
288 * Constructor.
289 *
290 * \param Wmin minimal wavelength for integration [nm]
291 * \param Wmax maximal wavelength for integration [nm]
292 * \param numberOfPoints number of points for integration
293 * \param epsilon precision of points for integration
294 */
295 JPDF(const double Wmin,
296 const double Wmax,
297 const int numberOfPoints = 20,
298 const double epsilon = 1e-12) :
299 JTOOLS::JGaussLegendre(numberOfPoints,epsilon),
300 wmin(Wmin),
301 wmax(Wmax)
302 {
303 // phi integration points
304
305 const double dx = PI / numberOfPoints;
306
307 for (double x = 0.5*dx; x < PI; x += dx) {
308 phd.push_back(element_type(cos(x),sin(x)));
309 }
310 }
311
312
313 /**
314 * Virtual destructor.
315 */
316 virtual ~JPDF()
317 {}
318
319
320 /**
321 * Number of Cherenkov photons per unit track length
322 *
323 * \return number of photons per unit track length [m^-1]
324 */
325 double getNumberOfPhotons() const
326 {
327 double value = 0.0;
328
329 const double xmin = 1.0 / wmax;
330 const double xmax = 1.0 / wmin;
331
332 for (const_iterator i = begin(); i != end(); ++i) {
333
334 const double x = 0.5 * (xmax + xmin) + i->getX() * 0.5 * (xmax - xmin);
335 const double dx = i->getY() * 0.5 * (xmax - xmin);
336
337 const double w = 1.0 / x;
338 const double dw = dx * w*w;
339
340 const double n = getIndexOfRefractionPhase(w);
341
342 value += cherenkov(w,n) * dw;
343 }
344
345 return value;
346 }
347
348
349 /**
350 * Number of photo-electrons from direct Cherenkov light from muon.
351 *
352 * \param R_m distance between muon and PMT [m]
353 * \param theta zenith angle orientation PMT [rad]
354 * \param phi azimuth angle orientation PMT [rad]
355 * \return P [npe]
356 */
357 double getDirectLightFromMuon(const double R_m,
358 const double theta,
359 const double phi) const
360 {
361 double value = 0;
362
363 const double R = std::max(R_m, getRmin());
364 const double A = getPhotocathodeArea();
365
366 const double px = sin(theta)*cos(phi);
367 const double pz = cos(theta);
368
369 for (const_iterator m = begin(); m != end(); ++m) {
370
371 const double w = 0.5 * (wmax + wmin) + m->getX() * 0.5 * (wmax - wmin);
372 const double dw = m->getY() * 0.5 * (wmax - wmin);
373
374 const double n = getIndexOfRefractionPhase(w);
375
376 const double l_abs = getAbsorptionLength(w);
377 const double ls = getScatteringLength(w);
378
379 const double npe = cherenkov(w,n) * dw * getQE(w);
380
381 const double ct0 = 1.0 / n;
382 const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
383
384 const double d = R / st0; // distance traveled by photon
385 const double ct = st0*px + ct0*pz; // cosine angle of incidence on PMT
386
387 const double U = getAngularAcceptance(ct); // PMT angular acceptance
388 const double V = exp(-d/l_abs) * exp(-d/ls); // absorption & scattering
389 const double W = A / (2.0*PI*R*st0); // solid angle
390
391 value += npe * U * V * W;
392 }
393
394 return value;
395 }
396
397
398 /**
399 * Probability density function for direct light from muon.
400 *
401 * \param R_m distance between muon and PMT [m]
402 * \param theta zenith angle orientation PMT [rad]
403 * \param phi azimuth angle orientation PMT [rad]
404 * \param t_ns time difference relative to direct Cherenkov light [ns]
405 * \return dP/dt [npe/ns]
406 */
407 double getDirectLightFromMuon(const double R_m,
408 const double theta,
409 const double phi,
410 const double t_ns) const
411 {
412 static const int N = 100; // maximal number of iterations
413 static const double eps = 1.0e-6; // precision index of refraction
414
415 const double R = std::max(R_m, getRmin());
416 const double t = R * getTanThetaC() / C + t_ns; // time [ns]
417 const double a = C*t/R; // target value
418 const double A = getPhotocathodeArea();
419
420 const double px = sin(theta)*cos(phi);
421 const double pz = cos(theta);
422
423 // check validity range for index of refraction
424
425 for (double buffer[] = { wmin, wmax, 0.0 }, *p = buffer; *p != 0.0; ++p) {
426
427 const double n = getIndexOfRefractionPhase(*p);
428 const double ng = getIndexOfRefractionGroup(*p);
429
430 const double ct0 = 1.0 / n;
431 const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
432
433 const double b = (ng - ct0) / st0; // running value
434
435 if (*p == wmin && b < a) { return 0.0; }
436 if (*p == wmax && b > a) { return 0.0; }
437 }
438
439
440 double umin = wmin;
441 double umax = wmax;
442
443 for (int i = 0; i != N; ++i) { // binary search for wavelength
444
445 const double w = 0.5 * (umin + umax);
446
447 const double n = getIndexOfRefractionPhase(w);
448 const double ng = getIndexOfRefractionGroup(w);
449
450 const double ct0 = 1.0 / n;
451 const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
452
453 const double b = (ng - ct0) / st0; // running value
454
455 if (fabs(b-a) < a*eps) {
456
457 const double np = getDispersionPhase(w);
458 const double ngp = getDispersionGroup(w);
459
460 const double l_abs = getAbsorptionLength(w);
461 const double ls = getScatteringLength(w);
462
463 const double d = R / st0; // distance traveled by photon
464 const double ct = st0*px + ct0*pz; // cosine angle of incidence on PMT
465
466 const double i3 = ct0*ct0*ct0 / (st0*st0*st0);
467
468 const double U = getAngularAcceptance(ct); // PMT angular acceptance
469 const double V = exp(-d/l_abs - d/ls); // absorption & scattering
470 const double W = A / (2.0*PI*R*st0); // solid angle
471
472 const double Ja = R * (ngp/st0 + np*(n-ng)*i3) / C; // dt/dlambda
473
474 return cherenkov(w,n) * getQE(w) * U * V * W / fabs(Ja);
475 }
476
477 if (b > a)
478 umin = w;
479 else
480 umax = w;
481 }
482
483 return 0.0;
484 }
485
486
487 /**
488 * Probability density function for scattered light from muon.
489 *
490 * \param R_m distance between muon and PMT [m]
491 * \param theta zenith angle orientation PMT [rad]
492 * \param phi azimuth angle orientation PMT [rad]
493 * \param t_ns time difference relative to direct Cherenkov light [ns]
494 * \return dP/dt [npe/ns]
495 */
496 double getScatteredLightFromMuon(const double R_m,
497 const double theta,
498 const double phi,
499 const double t_ns) const
500 {
501 static const double eps = 1.0e-10;
502
503 double value = 0;
504
505 const double R = std::max(R_m, getRmin());
506 const double t = R * getTanThetaC() / C + t_ns; // time [ns]
507 const double A = getPhotocathodeArea();
508
509 const double px = sin(theta)*cos(phi);
510 const double py = sin(theta)*sin(phi);
511 const double pz = cos(theta);
512
513 const double n0 = getIndexOfRefractionGroup(wmax);
514 const double n1 = getIndexOfRefractionGroup(wmin);
515 const double ni = sqrt(R*R + C*t*C*t) / R; // maximal index of refraction
516
517 if (n0 >= ni) {
518 return value;
519 }
520
521 const double nj = std::min(ni,n1);
522
523 double w = wmax;
524
525 for (const_iterator i = begin(); i != end(); ++i) {
526
527 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
528 const double dn = i->getY() * 0.5 * (nj - n0);
529
530 w = getWavelength(ng, w, 1.0e-5);
531
532 const double dw = dn / fabs(getDispersionGroup(w));
533
534 const double n = getIndexOfRefractionPhase(w);
535
536 const double l_abs = getAbsorptionLength(w);
537 const double ls = getScatteringLength(w);
538
539 const double npe = cherenkov(w,n) * dw * getQE(w);
540
541 if (npe <= 0) { continue; }
542
543 const double Jc = 1.0 / ls; // dN/dx
544
545 const double ct0 = 1.0 / n; // photon direction before scattering
546 const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
547
548 JRoot rz(R, ng, t); // square root
549
550 if (!rz.is_valid) { continue; }
551
552 const double zmin = rz.first;
553 const double zmax = rz.second;
554
555 const double zap = 1.0 / l_abs;
556
557 const double xmin = exp(zap*zmax);
558 const double xmax = exp(zap*zmin);
559
560 for (const_iterator j = begin(); j != end(); ++j) {
561
562 const double x = 0.5 * (xmax + xmin) + j->getX() * 0.5 * (xmax - xmin);
563 const double dx = j->getY() * 0.5 * (xmax - xmin);
564
565 const double z = log(x) / zap;
566 const double dz = -dx / (zap*x);
567
568 const double D = sqrt(z*z + R*R);
569 const double cd = -z / D;
570 const double sd = R / D;
571
572 const double d = (C*t - z) / ng; // photon path
573
574 //const double V = exp(-d/l_abs); // absorption
575
576 const double cta = cd*ct0 + sd*st0;
577 const double dca = d - 0.5*(d+D)*(d-D) / (d - D*cta);
578 const double tip = -log(D*D/(dca*dca) + eps) / PI;
579
580 const double ymin = exp(tip*PI);
581 const double ymax = 1.0;
582
583 for (const_iterator k = begin(); k != end(); ++k) {
584
585 const double y = 0.5 * (ymax + ymin) + k->getX() * 0.5 * (ymax - ymin);
586 const double dy = k->getY() * 0.5 * (ymax - ymin);
587
588 const double phi = log(y) / tip;
589 const double dp = -dy / (tip*y);
590
591 const double cp0 = cos(phi);
592 const double sp0 = sin(phi);
593
594 const double u = (R*R + z*z - d*d) / (2*R*st0*cp0 - 2*z*ct0 - 2*d);
595 const double v = d - u;
596
597 if (u <= 0) { continue; }
598 if (v <= 0) { continue; }
599
600 const double vi = 1.0 / v;
601 const double cts = (R*st0*cp0 - z*ct0 - u)*vi; // cosine scattering angle
602
603 const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
604
605 if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
606
607 const double vx = R - u*st0*cp0;
608 const double vy = - u*st0*sp0;
609 const double vz = -z - u*ct0;
610
611 const double ct[] = { // cosine angle of incidence on PMT
612 (vx*px + vy*py + vz*pz) * vi,
613 (vx*px - vy*py + vz*pz) * vi
614 };
615
616 const double U =
617 getAngularAcceptance(ct[0]) +
618 getAngularAcceptance(ct[1]); // PMT angular acceptance
619 const double W = std::min(A*vi*vi, 2.0*PI); // solid angle
620
621 const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
622 const double Jd = ng * (1.0 - cts) / C; // dt/du
623
624 value += (npe * dz * dp / (2*PI)) * U * V * W * Ja * Jc / fabs(Jd);
625 }
626 }
627 }
628
629 return value;
630 }
631
632
633 /**
634 * Probability density function for direct light from EM-showers.
635 *
636 * \param R_m distance between muon and PMT [m]
637 * \param theta zenith angle orientation PMT [rad]
638 * \param phi azimuth angle orientation PMT [rad]
639 * \param t_ns time difference relative to direct Cherenkov light [ns]
640 * \return dP/dt [npe/ns/GeV]
641 */
642 double getDirectLightFromEMshowers(const double R_m,
643 const double theta,
644 const double phi,
645 const double t_ns) const
646 {
647 double value = 0;
648
649 const double R = std::max(R_m, getRmin());
650 const double t = R * getTanThetaC() / C + t_ns; // time [ns]
651 const double A = getPhotocathodeArea();
652 const double D = 2.0*sqrt(A/PI);
653
654 const double px = sin(theta)*cos(phi);
655 //const double py = sin(theta)*sin(phi);
656 const double pz = cos(theta);
657
658 const double n0 = getIndexOfRefractionGroup(wmax);
659 const double n1 = getIndexOfRefractionGroup(wmin);
660 const double ni = sqrt(R*R + C*t*C*t) / R; // maximal index of refraction
661
662 if (n0 >= ni) {
663 return value;
664 }
665
666 const double nj = std::min(n1, ni);
667
668 double w = wmax;
669
670 for (const_iterator i = begin(); i != end(); ++i) {
671
672 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
673 const double dn = i->getY() * 0.5 * (nj - n0);
674
675 w = getWavelength(ng, w, 1.0e-5);
676
677 const double dw = dn / fabs(getDispersionGroup(w));
678
679 const double n = getIndexOfRefractionPhase(w);
680
681 const double l_abs = getAbsorptionLength(w);
682 const double ls = getScatteringLength(w);
683
684 const double npe = cherenkov(w,n) * dw * getQE(w);
685
686 JRoot rz(R, ng, t); // square root
687
688 if (!rz.is_valid) { continue; }
689
690 for (int j = 0; j != 2; ++j) {
691
692 const double z = rz[j];
693
694 if (C*t <= z) continue;
695
696 const double d = sqrt(z*z + R*R); // distance traveled by photon
697
698 const double ct0 = -z / d;
699 const double st0 = R / d;
700
701 const double dz = D / st0; // average track length
702
703 const double ct = st0*px + ct0*pz; // cosine angle of incidence on PMT
704
705 const double U = getAngularAcceptance(ct); // PMT angular acceptance
706 const double V = exp(-d/l_abs - d/ls); // absorption & scattering
707 const double W = A/(d*d); // solid angle
708
709 const double Ja = geant(n,ct0); // d^2N/dcos/dphi
710 const double Jd = (1.0 - ng*ct0) / C; // d t/ dz
711 const double Je = ng*st0*st0*st0 / (R*C); // d^2t/(dz)^2
712
713 value += gWater() * npe * U * V * W * Ja / (fabs(Jd) + 0.5*Je*dz);
714 }
715 }
716
717 return value;
718 }
719
720
721 /**
722 * Probability density function for scattered light from EM-showers.
723 *
724 * \param R_m distance between muon and PMT [m]
725 * \param theta zenith angle orientation PMT [rad]
726 * \param phi azimuth angle orientation PMT [rad]
727 * \param t_ns time difference relative to direct Cherenkov light [ns]
728 * \return dP/dt [npe/ns/GeV]
729 */
730 double getScatteredLightFromEMshowers(const double R_m,
731 const double theta,
732 const double phi,
733 const double t_ns) const
734 {
735 double value = 0;
736
737 const double R = std::max(R_m, getRmin());
738 const double t = R * getTanThetaC() / C + t_ns; // time [ns]
739 const double A = getPhotocathodeArea();
740
741 const double px = sin(theta)*cos(phi);
742 const double py = sin(theta)*sin(phi);
743 const double pz = cos(theta);
744
745 const double n0 = getIndexOfRefractionGroup(wmax);
746 const double n1 = getIndexOfRefractionGroup(wmin);
747 const double ni = sqrt(R*R + C*t*C*t) / R; // maximal index of refraction
748
749 if (n0 >= ni) {
750 return value;
751 }
752
753 const double nj = std::min(ni,n1);
754
755 double w = wmax;
756
757 for (const_iterator i = begin(); i != end(); ++i) {
758
759 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
760 const double dn = i->getY() * 0.5 * (nj - n0);
761
762 w = getWavelength(ng, w, 1.0e-5);
763
764 const double dw = dn / fabs(getDispersionGroup(w));
765
766 const double n = getIndexOfRefractionPhase(w);
767
768 const double l_abs = getAbsorptionLength(w);
769 const double ls = getScatteringLength(w);
770
771 const double npe = cherenkov(w,n) * dw * getQE(w);
772
773 if (npe <= 0) { continue; }
774
775 const double Jc = 1.0 / ls; // dN/dx
776
777 JRoot rz(R, ng, t); // square root
778
779 if (!rz.is_valid) { continue; }
780
781 const double zmin = rz.first;
782 const double zmax = rz.second;
783
784 const double zap = 1.0 / l_abs;
785
786 const double xmin = exp(zap*zmax);
787 const double xmax = exp(zap*zmin);
788
789 for (const_iterator j = begin(); j != end(); ++j) {
790
791 const double x = 0.5 * (xmax + xmin) + j->getX() * 0.5 * (xmax - xmin);
792 const double dx = j->getY() * 0.5 * (xmax - xmin);
793
794 const double z = log(x) / zap;
795 const double dz = -dx / (zap*x);
796
797 const double D = sqrt(z*z + R*R);
798 const double cd = -z / D;
799 const double sd = R / D;
800
801 const double qx = cd*px + 0 - sd*pz;
802 const double qy = 1*py;
803 const double qz = sd*px + 0 + cd*pz;
804
805 const double d = (C*t - z) / ng; // photon path
806
807 //const double V = exp(-d/l_abs); // absorption
808
809 const double ds = 2.0 / (size() + 1);
810
811 for (double sb = 0.5*ds; sb < 1.0 - 0.25*ds; sb += ds) {
812
813 for (int buffer[] = { -1, +1, 0}, *k = buffer; *k != 0; ++k) {
814
815 const double cb = (*k) * sqrt((1.0 + sb)*(1.0 - sb));
816 const double dcb = (*k) * ds*sb/cb;
817
818 const double v = 0.5 * (d + D) * (d - D) / (d - D*cb);
819 const double u = d - v;
820
821 if (u <= 0) { continue; }
822 if (v <= 0) { continue; }
823
824 const double cts = (D*cb - v) / u; // cosine scattering angle
825
826 const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
827
828 if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
829
830 const double W = std::min(A/(v*v), 2.0*PI); // solid angle
831 const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
832 const double Jd = ng * (1.0 - cts) / C; // dt/du
833
834 const double ca = (D - v*cb) / u;
835 const double sa = v*sb / u;
836
837 const double dp = PI / phd.size();
838 const double dom = dcb*dp * v*v / (u*u);
839
840 for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
841
842 const double cp = l->getX();
843 const double sp = l->getY();
844
845 const double ct0 = cd*ca - sd*sa*cp;
846
847 const double vx = sb*cp * qx;
848 const double vy = sb*sp * qy;
849 const double vz = cb * qz;
850
851 const double U =
852 getAngularAcceptance(vx + vy + vz) +
853 getAngularAcceptance(vx - vy + vz); // PMT angular acceptance
854
855 const double Jb = geant(n,ct0); // d^2N/dcos/dphi
856
857 value += dom * gWater() * npe * dz * U * V * W * Ja * Jb * Jc / fabs(Jd);
858 }
859 }
860 }
861 }
862 }
863
864 return value;
865 }
866
867
868 /**
869 * Probability density function for scattered light from muon.
870 *
871 * \param D_m distance between track segment and PMT [m]
872 * \param cd cosine angle muon direction and track segment - PMT position
873 * \param theta zenith angle orientation PMT [rad]
874 * \param phi azimuth angle orientation PMT [rad]
875 * \param t_ns time difference relative to direct Cherenkov light [ns]
876 * \return d^2P/dt/dx [npe/ns/m]
877 */
878 double getScatteredLightFromMuon(const double D_m,
879 const double cd,
880 const double theta,
881 const double phi,
882 const double t_ns) const
883 {
884 static const double eps = 1.0e-10;
885
886 double value = 0;
887
888 const double sd = sqrt((1.0 + cd)*(1.0 - cd));
889 const double D = std::max(D_m, getRmin());
890 const double R = sd * D; // minimal distance of approach [m]
891 const double Z = -cd * D; // photon emission point
892 const double L = D;
893 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
894 const double A = getPhotocathodeArea();
895
896 const double px = sin(theta)*cos(phi);
897 const double py = sin(theta)*sin(phi);
898 const double pz = cos(theta);
899
900 const double n0 = getIndexOfRefractionGroup(wmax);
901 const double n1 = getIndexOfRefractionGroup(wmin);
902
903 const double ni = C*t / L; // maximal index of refraction
904
905 if (n0 >= ni) {
906 return value;
907 }
908
909 const double nj = std::min(ni,n1);
910
911 double w = wmax;
912
913 for (const_iterator i = begin(); i != end(); ++i) {
914
915 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
916 const double dn = i->getY() * 0.5 * (nj - n0);
917
918 w = getWavelength(ng, w, 1.0e-5);
919
920 const double dw = dn / fabs(getDispersionGroup(w));
921
922 const double n = getIndexOfRefractionPhase(w);
923
924 const double l_abs = getAbsorptionLength(w);
925 const double ls = getScatteringLength(w);
926
927 const double npe = cherenkov(w,n) * dw * getQE(w);
928
929 if (npe <= 0) { continue; }
930
931 const double Jc = 1.0 / ls; // dN/dx
932
933 const double ct0 = 1.0 / n; // photon direction before scattering
934 const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
935
936 const double d = C*t / ng; // photon path
937
938 //const double V = exp(-d/l_abs); // absorption
939
940 const double cta = cd*ct0 + sd*st0;
941 const double dca = d - 0.5*(d+L)*(d-L) / (d - L*cta);
942 const double tip = -log(L*L/(dca*dca) + eps) / PI;
943
944 const double ymin = exp(tip*PI);
945 const double ymax = 1.0;
946
947 for (const_iterator j = begin(); j != end(); ++j) {
948
949 const double y = 0.5 * (ymax + ymin) + j->getX() * 0.5 * (ymax - ymin);
950 const double dy = j->getY() * 0.5 * (ymax - ymin);
951
952 const double phi = log(y) / tip;
953 const double dp = -dy / (tip*y);
954
955 const double cp0 = cos(phi);
956 const double sp0 = sin(phi);
957
958 const double u = (R*R + Z*Z - d*d) / (2*R*st0*cp0 - 2*Z*ct0 - 2*d);
959 const double v = d - u;
960
961 if (u <= 0) { continue; }
962 if (v <= 0) { continue; }
963
964 const double vi = 1.0 / v;
965 const double cts = (R*st0*cp0 - Z*ct0 - u)*vi; // cosine scattering angle
966
967 const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
968
969 if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
970
971 const double vx = R - u*st0*cp0;
972 const double vy = - u*st0*sp0;
973 const double vz = -Z - u*ct0;
974
975 const double ct[] = { // cosine angle of incidence on PMT
976 (vx*px + vy*py + vz*pz) * vi,
977 (vx*px - vy*py + vz*pz) * vi
978 };
979
980 const double U =
981 getAngularAcceptance(ct[0]) +
982 getAngularAcceptance(ct[1]); // PMT angular acceptance
983 const double W = std::min(A*vi*vi, 2.0*PI); // solid angle
984
985 const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
986 const double Jd = ng * (1.0 - cts) / C; // dt/du
987
988 value += (npe * dp / (2*PI)) * U * V * W * Ja * Jc / fabs(Jd);
989 }
990 }
991
992 return value;
993 }
994
995
996 /**
997 * Probability density function for direct light from EM-shower.
998 *
999 * \param D_m distance between EM-shower and PMT [m]
1000 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1001 * \param theta zenith angle orientation PMT [rad]
1002 * \param phi azimuth angle orientation PMT [rad]
1003 * \param t_ns time difference relative to direct Cherenkov light [ns]
1004 * \return d^2P/dt/dE [npe/ns/GeV]
1005 */
1006 double getDirectLightFromEMshower(const double D_m,
1007 const double cd,
1008 const double theta,
1009 const double phi,
1010 const double t_ns) const
1011 {
1012 const double ct0 = cd;
1013 const double st0 = sqrt((1.0 + ct0)*(1.0 - ct0));
1014
1015 const double D = std::max(D_m, getRmin());
1016 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1017 const double A = getPhotocathodeArea();
1018
1019 const double px = sin(theta)*cos(phi);
1020 const double pz = cos(theta);
1021
1022 const double n0 = getIndexOfRefractionGroup(wmax);
1023 const double n1 = getIndexOfRefractionGroup(wmin);
1024 const double ng = C*t/D; // index of refraction
1025
1026 if (n0 >= ng) { return 0.0; }
1027 if (n1 <= ng) { return 0.0; }
1028
1029 const double w = getWavelength(ng);
1030 const double n = getIndexOfRefractionPhase(w);
1031
1032 const double l_abs = getAbsorptionLength(w);
1033 const double ls = getScatteringLength(w);
1034
1035 const double npe = cherenkov(w,n) * getQE(w);
1036
1037 const double ct = st0*px + ct0*pz; // cosine angle of incidence on PMT
1038
1039 const double U = getAngularAcceptance(ct); // PMT angular acceptance
1040 const double V = exp(-D/l_abs - D/ls); // absorption & scattering
1041 const double W = A/(D*D); // solid angle
1042
1043 const double ngp = getDispersionGroup(w);
1044
1045 const double Ja = D * ngp / C; // dt/dlambda
1046 const double Jb = geant(n,ct0); // d^2N/dcos/dphi
1047
1048 return npe * geanc() * U * V * W * Jb / fabs(Ja);
1049 }
1050
1051
1052 /**
1053 * Probability density function for scattered light from EM-shower.
1054 *
1055 * \param D_m distance between EM-shower and PMT [m]
1056 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1057 * \param theta zenith angle orientation PMT [rad]
1058 * \param phi azimuth angle orientation PMT [rad]
1059 * \param t_ns time difference relative to direct Cherenkov light [ns]
1060 * \return d^2P/dt/dE [npe/ns/GeV]
1061 */
1062 double getScatteredLightFromEMshower(const double D_m,
1063 const double cd,
1064 const double theta,
1065 const double phi,
1066 const double t_ns) const
1067 {
1068 double value = 0;
1069
1070 const double sd = sqrt((1.0 + cd)*(1.0 - cd));
1071 const double D = std::max(D_m, getRmin());
1072 const double L = D;
1073 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1074
1075 const double A = getPhotocathodeArea();
1076
1077 const double px = sin(theta)*cos(phi);
1078 const double py = sin(theta)*sin(phi);
1079 const double pz = cos(theta);
1080
1081 const double qx = cd*px + 0 - sd*pz;
1082 const double qy = 1*py;
1083 const double qz = sd*px + 0 + cd*pz;
1084
1085 const double n0 = getIndexOfRefractionGroup(wmax);
1086 const double n1 = getIndexOfRefractionGroup(wmin);
1087
1088 const double ni = C*t / L; // maximal index of refraction
1089
1090 if (n0 >= ni) {
1091 return value;
1092 }
1093
1094 const double nj = std::min(ni,n1);
1095
1096 double w = wmax;
1097
1098 for (const_iterator i = begin(); i != end(); ++i) {
1099
1100 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
1101 const double dn = i->getY() * 0.5 * (nj - n0);
1102
1103 w = getWavelength(ng, w, 1.0e-5);
1104
1105 const double dw = dn / fabs(getDispersionGroup(w));
1106
1107 const double n = getIndexOfRefractionPhase(w);
1108
1109 const double l_abs = getAbsorptionLength(w);
1110 const double ls = getScatteringLength(w);
1111
1112 const double npe = cherenkov(w,n) * dw * getQE(w);
1113
1114 if (npe <= 0) { continue; }
1115
1116 const double Jc = 1.0 / ls; // dN/dx
1117
1118 const double d = C*t / ng; // photon path
1119
1120 //const double V = exp(-d/l_abs); // absorption
1121
1122 const double ds = 2.0 / (size() + 1);
1123
1124 for (double sb = 0.5*ds; sb < 1.0 - 0.25*ds; sb += ds) {
1125
1126 for (int buffer[] = { -1, +1, 0}, *k = buffer; *k != 0; ++k) {
1127
1128 const double cb = (*k) * sqrt((1.0 + sb)*(1.0 - sb));
1129 const double dcb = (*k) * ds*sb/cb;
1130
1131 const double v = 0.5 * (d + L) * (d - L) / (d - L*cb);
1132 const double u = d - v;
1133
1134 if (u <= 0) { continue; }
1135 if (v <= 0) { continue; }
1136
1137 const double cts = (L*cb - v) / u; // cosine scattering angle
1138
1139 const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
1140
1141 if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
1142
1143 const double W = std::min(A/(v*v), 2.0*PI); // solid angle
1144 const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
1145 const double Jd = ng * (1.0 - cts) / C; // dt/du
1146
1147 const double ca = (L - v*cb) / u;
1148 const double sa = v*sb / u;
1149
1150 const double dp = PI / phd.size();
1151 const double dom = dcb*dp * v*v / (u*u);
1152 const double dos = sqrt(dom);
1153
1154 for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
1155
1156 const double cp = l->getX();
1157 const double sp = l->getY();
1158
1159 const double ct0 = cd*ca - sd*sa*cp;
1160
1161 const double vx = -sb*cp * qx;
1162 const double vy = -sb*sp * qy;
1163 const double vz = cb * qz;
1164
1165 const double U =
1166 getAngularAcceptance(vx + vy + vz) +
1167 getAngularAcceptance(vx - vy + vz); // PMT angular acceptance
1168
1169 //const double Jb = geant(n,ct0); // d^2N/dcos/dphi
1170
1171 //value += npe * geanc() * dom * U * V * W * Ja * Jb * Jc / fabs(Jd);
1172
1173 const double Jb = geant(n,
1174 ct0 - 0.5*dos,
1175 ct0 + 0.5*dos); // dN/dphi
1176
1177 value += npe * geanc() * dos * U * V * W * Ja * Jb * Jc / fabs(Jd);
1178 }
1179 }
1180 }
1181 }
1182
1183 return value;
1184 }
1185
1186
1187 /**
1188 * Probability density function for direct light from EM-shower.
1189 *
1190 * \param E EM-shower energy [GeV]
1191 * \param D_m distance between EM-shower and PMT [m]
1192 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1193 * \param theta zenith angle orientation PMT [rad]
1194 * \param phi azimuth angle orientation PMT [rad]
1195 * \param t_ns time difference relative to direct Cherenkov light [ns]
1196 * \return dP/dt [npe/ns]
1197 */
1198 double getDirectLightFromEMshower(const double E,
1199 const double D_m,
1200 const double cd,
1201 const double theta,
1202 const double phi,
1203 const double t_ns) const
1204 {
1205 double value = 0;
1206
1207 const double sd = sqrt((1.0 + cd)*(1.0 - cd));
1208 const double D = std::max(D_m, getRmin());
1209 const double R = D * sd; // minimal distance of approach [m]
1210 const double Z = -D * cd;
1211 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1212
1213 const double n0 = getIndexOfRefractionGroup(wmax);
1214 const double n1 = getIndexOfRefractionGroup(wmin);
1215
1216 double zmin = 0.0; // minimal shower length [m]
1217 double zmax = geanz.getLength(E, 1.0); // maximal shower length [m]
1218
1219 const double d = sqrt((Z+zmax)*(Z+zmax) + R*R);
1220
1221 if (C*t > std::max(n1*D, zmax + n1*d)) {
1222 return value;
1223 }
1224
1225 if (C*t < n0*D) {
1226
1227 JRoot rz(R, n0, t + Z/C); // square root
1228
1229 if (!rz.is_valid) {
1230 return value;
1231 }
1232
1233 if (rz.second > Z) { zmin = rz.second - Z; }
1234 if (rz.first > Z) { zmin = rz.first - Z; }
1235 }
1236
1237 if (C*t > n1*D) {
1238
1239 JRoot rz(R, n1, t + Z/C); // square root
1240
1241 if (!rz.is_valid) {
1242 return value;
1243 }
1244
1245 if (rz.second > Z) { zmin = rz.second - Z; }
1246 if (rz.first > Z) { zmin = rz.first - Z; }
1247 }
1248
1249 if (C*t < zmax + n0*d) {
1250
1251 JRoot rz(R, n0, t + Z/C); // square root
1252
1253 if (!rz.is_valid) {
1254 return value;
1255 }
1256
1257 if (rz.first > Z) { zmax = rz.first - Z; }
1258 if (rz.second > Z) { zmax = rz.second - Z; }
1259 }
1260
1261 if (zmin < 0.0) {
1262 zmin = 0.0;
1263 }
1264
1265 if (zmax > zmin) {
1266
1267 const double ymin = geanz.getIntegral(E, zmin);
1268 const double ymax = geanz.getIntegral(E, zmax);
1269 const double dy = (ymax - ymin) / size();
1270
1271 if (dy > 2*std::numeric_limits<double>::epsilon()) {
1272
1273 for (double y = ymin + 0.5*dy; y < ymax; y += dy) {
1274
1275 const double z = Z + geanz.getLength(E, y);
1276 const double d = sqrt(R*R + z*z);
1277 const double t1 = t + (Z - z) / C - d * getIndexOfRefraction() / C;
1278
1279 value += dy * E * getDirectLightFromEMshower(d, -z/d, theta, phi, t1);
1280 }
1281
1282 }
1283 }
1284
1285 return value;
1286 }
1287
1288
1289 /**
1290 * Probability density function for scattered light from EM-shower.
1291 *
1292 * \param E EM-shower energy [GeV]
1293 * \param D_m distance between EM-shower and PMT [m]
1294 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1295 * \param theta zenith angle orientation PMT [rad]
1296 * \param phi azimuth angle orientation PMT [rad]
1297 * \param t_ns time difference relative to direct Cherenkov light [ns]
1298 * \return dP/dt [npe/ns]
1299 */
1300 double getScatteredLightFromEMshower(const double E,
1301 const double D_m,
1302 const double cd,
1303 const double theta,
1304 const double phi,
1305 const double t_ns) const
1306 {
1307 double value = 0;
1308
1309 const double sd = sqrt((1.0 + cd)*(1.0 - cd));
1310 const double D = std::max(D_m, getRmin());
1311 const double R = D * sd; // minimal distance of approach [m]
1312 const double Z = -D * cd;
1313 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1314
1315 const double n0 = getIndexOfRefractionGroup(wmax);
1316
1317 double zmin = 0.0; // minimal shower length [m]
1318 double zmax = geanz.getLength(E, 1.0); // maximal shower length [m]
1319
1320 const double d = sqrt((Z+zmax)*(Z+zmax) + R*R);
1321
1322 if (C*t < n0*D) {
1323
1324 JRoot rz(R, n0, t + Z/C); // square root
1325
1326 if (!rz.is_valid) {
1327 return value;
1328 }
1329
1330 if (rz.second > Z) { zmin = rz.second - Z; }
1331 if (rz.first > Z) { zmin = rz.first - Z; }
1332 }
1333
1334 if (C*t < zmax + n0*d) {
1335
1336 JRoot rz(R, n0, t + Z/C); // square root
1337
1338 if (!rz.is_valid) {
1339 return value;
1340 }
1341
1342 if (rz.first > Z) { zmax = rz.first - Z; }
1343 if (rz.second > Z) { zmax = rz.second - Z; }
1344 }
1345
1346 const double ymin = geanz.getIntegral(E, zmin);
1347 const double ymax = geanz.getIntegral(E, zmax);
1348 const double dy = (ymax - ymin) / size();
1349
1350 if (dy > 2*std::numeric_limits<double>::epsilon()) {
1351
1352 for (double y = ymin + 0.5*dy; y < ymax; y += dy) {
1353
1354 const double z = Z + geanz.getLength(E, y);
1355 const double d = sqrt(R*R + z*z);
1356 const double t1 = t + (Z - z) / C - d * getIndexOfRefraction() / C;
1357
1358 value += dy * E * getScatteredLightFromEMshower(d, -z/d, theta, phi, t1);
1359 }
1360 }
1361
1362 return value;
1363 }
1364
1365
1366 /**
1367 * Probability density function for direct light from delta-rays.
1368 *
1369 * \param R_m distance between muon and PMT [m]
1370 * \param theta zenith angle orientation PMT [rad]
1371 * \param phi azimuth angle orientation PMT [rad]
1372 * \param t_ns time difference relative to direct Cherenkov light [ns]
1373 * \return d^2P/dt/dE [npe/ns x m/GeV]
1374 */
1375 double getDirectLightFromDeltaRays(const double R_m,
1376 const double theta,
1377 const double phi,
1378 const double t_ns) const
1379 {
1380 double value = 0;
1381
1382 const double R = std::max(R_m, getRmin());
1383 const double t = R * getTanThetaC() / C + t_ns; // time [ns]
1384 const double A = getPhotocathodeArea();
1385 const double D = 2.0*sqrt(A/PI);
1386
1387 const double px = sin(theta)*cos(phi);
1388 const double pz = cos(theta);
1389
1390 const double n0 = getIndexOfRefractionGroup(wmax);
1391 const double n1 = getIndexOfRefractionGroup(wmin);
1392 const double ni = sqrt(R*R + C*t*C*t) / R; // maximal index of refraction
1393
1394 if (n0 >= ni) {
1395 return value;
1396 }
1397
1398 const double nj = std::min(n1, ni);
1399
1400 double w = wmax;
1401
1402 for (const_iterator i = begin(); i != end(); ++i) {
1403
1404 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
1405 const double dn = i->getY() * 0.5 * (nj - n0);
1406
1407 w = getWavelength(ng, w, 1.0e-5);
1408
1409 const double dw = dn / fabs(getDispersionGroup(w));
1410
1411 const double n = getIndexOfRefractionPhase(w);
1412
1413 const double l_abs = getAbsorptionLength(w);
1414 const double ls = getScatteringLength(w);
1415
1416 const double npe = cherenkov(w,n) * dw * getQE(w);
1417
1418 JRoot rz(R, ng, t); // square root
1419
1420 if (!rz.is_valid) { continue; }
1421
1422 for (int j = 0; j != 2; ++j) {
1423
1424 const double z = rz[j];
1425
1426 if (C*t <= z) continue;
1427
1428 const double d = sqrt(z*z + R*R); // distance traveled by photon
1429
1430 const double ct0 = -z / d;
1431 const double st0 = R / d;
1432
1433 const double dz = D / st0; // average track length
1434 const double ct = st0*px + ct0*pz; // cosine angle of incidence on PMT
1435
1436 const double U = getAngularAcceptance(ct); // PMT angular acceptance
1437 const double V = exp(-d/l_abs - d/ls); // absorption & scattering
1438 const double W = A/(d*d); // solid angle
1439
1440 const double Ja = getDeltaRayProbability(ct0); // d^2N/dcos/dphi
1441 const double Jd = (1.0 - ng*ct0) / C; // d t/ dz
1442 const double Je = ng*st0*st0*st0 / (R*C); // d^2t/(dz)^2
1443
1444 value += npe * geanc() * U * V * W * Ja / (fabs(Jd) + 0.5*Je*dz);
1445 }
1446 }
1447
1448 return value;
1449 }
1450
1451
1452 /**
1453 * Probability density function for scattered light from delta-rays.
1454 *
1455 * \param R_m distance between muon and PMT [m]
1456 * \param theta zenith angle orientation PMT [rad]
1457 * \param phi azimuth angle orientation PMT [rad]
1458 * \param t_ns time difference relative to direct Cherenkov light [ns]
1459 * \return d^2P/dt/dx [npe/ns x m/GeV]
1460 */
1461 double getScatteredLightFromDeltaRays(const double R_m,
1462 const double theta,
1463 const double phi,
1464 const double t_ns) const
1465 {
1466 double value = 0;
1467
1468 const double R = std::max(R_m, getRmin());
1469 const double t = R * getTanThetaC() / C + t_ns; // time [ns]
1470 const double A = getPhotocathodeArea();
1471
1472 const double px = sin(theta)*cos(phi);
1473 const double py = sin(theta)*sin(phi);
1474 const double pz = cos(theta);
1475
1476 const double n0 = getIndexOfRefractionGroup(wmax);
1477 const double n1 = getIndexOfRefractionGroup(wmin);
1478 const double ni = sqrt(R*R + C*t*C*t) / R; // maximal index of refraction
1479
1480 if (n0 >= ni) {
1481 return value;
1482 }
1483
1484 const double nj = std::min(ni,n1);
1485
1486 double w = wmax;
1487
1488 for (const_iterator i = begin(); i != end(); ++i) {
1489
1490 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
1491 const double dn = i->getY() * 0.5 * (nj - n0);
1492
1493 w = getWavelength(ng, w, 1.0e-5);
1494
1495 const double dw = dn / fabs(getDispersionGroup(w));
1496
1497 const double n = getIndexOfRefractionPhase(w);
1498
1499 const double l_abs = getAbsorptionLength(w);
1500 const double ls = getScatteringLength(w);
1501
1502 const double npe = cherenkov(w,n) * dw * getQE(w);
1503
1504 if (npe <= 0) { continue; }
1505
1506 const double Jc = 1.0 / ls; // dN/dx
1507
1508 JRoot rz(R, ng, t); // square root
1509
1510 if (!rz.is_valid) { continue; }
1511
1512 const double zmin = rz.first;
1513 const double zmax = rz.second;
1514
1515 const double zap = 1.0 / l_abs;
1516
1517 const double xmin = exp(zap*zmax);
1518 const double xmax = exp(zap*zmin);
1519
1520 for (const_iterator j = begin(); j != end(); ++j) {
1521
1522 const double x = 0.5 * (xmax + xmin) + j->getX() * 0.5 * (xmax - xmin);
1523 const double dx = j->getY() * 0.5 * (xmax - xmin);
1524
1525 const double z = log(x) / zap;
1526 const double dz = -dx / (zap*x);
1527
1528 const double D = sqrt(z*z + R*R);
1529 const double cd = -z / D;
1530 const double sd = R / D;
1531
1532 const double qx = cd*px + 0 - sd*pz;
1533 const double qy = 1*py;
1534 const double qz = sd*px + 0 + cd*pz;
1535
1536 const double d = (C*t - z) / ng; // photon path
1537
1538 //const double V = exp(-d/l_abs); // absorption
1539
1540 const double ds = 2.0 / (size() + 1);
1541
1542 for (double sb = 0.5*ds; sb < 1.0 - 0.25*ds; sb += ds) {
1543
1544 for (int buffer[] = { -1, +1, 0}, *k = buffer; *k != 0; ++k) {
1545
1546 const double cb = (*k) * sqrt((1.0 + sb)*(1.0 - sb));
1547 const double dcb = (*k) * ds*sb/cb;
1548
1549 const double v = 0.5 * (d + D) * (d - D) / (d - D*cb);
1550 const double u = d - v;
1551
1552 if (u <= 0) { continue; }
1553 if (v <= 0) { continue; }
1554
1555 const double cts = (D*cb - v) / u; // cosine scattering angle
1556
1557 const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
1558
1559 if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
1560
1561 const double W = std::min(A/(v*v), 2.0*PI); // solid angle
1562 const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
1563 const double Jd = ng * (1.0 - cts) / C; // dt/du
1564
1565 const double ca = (D - v*cb) / u;
1566 const double sa = v*sb / u;
1567
1568 const double dp = PI / phd.size();
1569 const double dom = dcb*dp * v*v / (u*u);
1570
1571 for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
1572
1573 const double cp = l->getX();
1574 const double sp = l->getY();
1575
1576 const double ct0 = cd*ca - sd*sa*cp;
1577
1578 const double vx = sb*cp * qx;
1579 const double vy = sb*sp * qy;
1580 const double vz = cb * qz;
1581
1582 const double U =
1583 getAngularAcceptance(vx + vy + vz) +
1584 getAngularAcceptance(vx - vy + vz); // PMT angular acceptance
1585
1586 const double Jb = getDeltaRayProbability(ct0); // d^2N/dcos/dphi
1587
1588 value += dom * npe * geanc() * dz * U * V * W * Ja * Jb * Jc / fabs(Jd);
1589 }
1590 }
1591 }
1592 }
1593 }
1594
1595 return value;
1596 }
1597
1598
1599 /**
1600 * Probability density function for direct light from isotropic light source.
1601 *
1602 * \param D_m distance between light source and PMT [m]
1603 * \param ct cosine angle PMT
1604 * \param t_ns time difference relative to direct Cherenkov light [ns]
1605 * \return d^2P/dt/dE [npe/ns/GeV]
1606 */
1607 double getDirectLightFromBrightPoint(const double D_m,
1608 const double ct,
1609 const double t_ns) const
1610 {
1611 const double D = std::max(D_m, getRmin());
1612 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1613 const double A = getPhotocathodeArea();
1614
1615 const double n0 = getIndexOfRefractionGroup(wmax);
1616 const double n1 = getIndexOfRefractionGroup(wmin);
1617 const double ng = C*t/D; // index of refraction
1618
1619 if (n0 >= ng) { return 0.0; }
1620 if (n1 <= ng) { return 0.0; }
1621
1622 const double w = getWavelength(ng);
1623 const double n = getIndexOfRefractionPhase(w);
1624
1625 const double l_abs = getAbsorptionLength(w);
1626 const double ls = getScatteringLength(w);
1627
1628 const double npe = cherenkov(w,n) * getQE(w);
1629
1630 const double U = getAngularAcceptance(ct); // PMT angular acceptance
1631 const double V = exp(-D/l_abs - D/ls); // absorption & scattering
1632 const double W = A/(D*D); // solid angle
1633
1634 const double ngp = getDispersionGroup(w);
1635
1636 const double Ja = D * ngp / C; // dt/dlambda
1637 const double Jb = 1.0 / (4.0*PI); // d^2N/dOmega
1638
1639 return npe * geanc() * U * V * W * Jb / fabs(Ja);
1640 }
1641
1642
1643 /**
1644 * Probability density function for scattered light from isotropic light source.
1645 *
1646 * \param D_m distance between light source and PMT [m]
1647 * \param ct cosine angle PMT
1648 * \param t_ns time difference relative to direct Cherenkov light [ns]
1649 * \return d^2P/dt/dE [npe/ns/GeV]
1650 */
1651 double getScatteredLightFromBrightPoint(const double D_m,
1652 const double ct,
1653 const double t_ns) const
1654 {
1655 double value = 0;
1656
1657 const double D = std::max(D_m, getRmin());
1658 const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1659 const double st = sqrt((1.0 + ct)*(1.0 - ct));
1660
1661 const double A = getPhotocathodeArea();
1662
1663 const double n0 = getIndexOfRefractionGroup(wmax);
1664 const double n1 = getIndexOfRefractionGroup(wmin);
1665
1666 const double ni = C*t / D; // maximal index of refraction
1667
1668 if (n0 >= ni) {
1669 return value;
1670 }
1671
1672 const double nj = std::min(ni,n1);
1673
1674 double w = wmax;
1675
1676 for (const_iterator i = begin(); i != end(); ++i) {
1677
1678 const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
1679 const double dn = i->getY() * 0.5 * (nj - n0);
1680
1681 w = getWavelength(ng, w, 1.0e-5);
1682
1683 const double dw = dn / fabs(getDispersionGroup(w));
1684
1685 const double n = getIndexOfRefractionPhase(w);
1686
1687 const double npe = cherenkov(w,n) * dw * getQE(w);
1688
1689 if (npe <= 0) { continue; }
1690
1691 const double l_abs = getAbsorptionLength(w);
1692 const double ls = getScatteringLength(w);
1693
1694 const double Jc = 1.0 / ls; // dN/dx
1695 const double Jb = 1.0 / (4.0*PI); // d^2N/dcos/dphi
1696
1697 const double d = C*t / ng; // photon path
1698
1699 const double dcb = 2.0 / (size() + 1);
1700
1701 for (double cb = -1.0 + 0.5*dcb; cb < +1.0; cb += dcb) {
1702
1703 const double sb = sqrt((1.0 + cb)*(1.0 - cb));
1704
1705 const double v = 0.5 * (d + D) * (d - D) / (d - D*cb);
1706 const double u = d - v;
1707
1708 if (u <= 0) { continue; }
1709 if (v <= 0) { continue; }
1710
1711 const double cts = (D*cb - v) / u; // cosine scattering angle
1712
1713 const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
1714
1715 if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
1716
1717 const double W = std::min(A/(v*v), 2.0*PI); // solid angle
1718 const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
1719 const double Jd = ng * (1.0 - cts) / C; // dt/du
1720
1721 const double dp = PI / phd.size();
1722 const double dom = dcb*dp * v*v / (u*u); // dOmega
1723
1724 for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
1725
1726 const double cp = l->getX();
1727 const double dot = cb*ct + sb*cp*st;
1728
1729 const double U = 2 * getAngularAcceptance(dot); // PMT angular acceptance
1730
1731 value += npe * geanc() * dom * U * V * W * Ja * Jb * Jc / fabs(Jd);
1732 }
1733 }
1734 }
1735
1736 return value;
1737 }
1738
1739
1740 /**
1741 * Probability density function for light from muon.
1742 *
1743 * \param type PDF type
1744 * \param E_GeV muon energy [GeV]
1745 * \param R_m distance between muon and PMT [m]
1746 * \param theta zenith angle orientation PMT [rad]
1747 * \param phi azimuth angle orientation PMT [rad]
1748 * \param t_ns time difference relative to direct Cherenkov light [ns]
1749 * \return d^2P/dt/dx [npe/ns]
1750 */
1751 double getLightFromMuon(const int type,
1752 const double E_GeV,
1753 const double R_m,
1754 const double theta,
1755 const double phi,
1756 const double t_ns) const
1757 {
1758 switch (type) {
1759
1760 case DIRECT_LIGHT_FROM_MUON:
1761 return getDirectLightFromMuon (R_m, theta, phi, t_ns);
1762
1763 case SCATTERED_LIGHT_FROM_MUON:
1764 return getScatteredLightFromMuon (R_m, theta, phi, t_ns);
1765
1766 case DIRECT_LIGHT_FROM_EMSHOWERS:
1767 return getDirectLightFromEMshowers (R_m, theta, phi, t_ns) * E_GeV;
1768
1769 case SCATTERED_LIGHT_FROM_EMSHOWERS:
1770 return getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV;
1771
1772 case DIRECT_LIGHT_FROM_DELTARAYS:
1773 return getDirectLightFromDeltaRays (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);
1774
1775 case SCATTERED_LIGHT_FROM_DELTARAYS:
1776 return getScatteredLightFromDeltaRays(R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);
1777
1778 default:
1779 return 0.0;
1780 }
1781 }
1782
1783
1784 /**
1785 * Probability density function for light from muon.
1786 *
1787 * \param E_GeV muon energy [GeV]
1788 * \param R_m distance between muon and PMT [m]
1789 * \param theta zenith angle orientation PMT [rad]
1790 * \param phi azimuth angle orientation PMT [rad]
1791 * \param t_ns time difference relative to direct Cherenkov light [ns]
1792 * \return d^2P/dt/dx [npe/ns]
1793 */
1794 double getLightFromMuon(const double E_GeV,
1795 const double R_m,
1796 const double theta,
1797 const double phi,
1798 const double t_ns) const
1799 {
1800 return (getDirectLightFromMuon (R_m, theta, phi, t_ns) +
1801 getScatteredLightFromMuon (R_m, theta, phi, t_ns) +
1802 getDirectLightFromEMshowers (R_m, theta, phi, t_ns) * E_GeV +
1803 getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV +
1804 getDirectLightFromDeltaRays (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV) +
1805 getScatteredLightFromDeltaRays(R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV));
1806 }
1807
1808
1809 /**
1810 * Probability density function for light from EM-shower.
1811 *
1812 * \param type PDF type
1813 * \param D_m distance between EM-shower and PMT [m]
1814 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1815 * \param theta zenith angle orientation PMT [rad]
1816 * \param phi azimuth angle orientation PMT [rad]
1817 * \param t_ns time difference relative to direct Cherenkov light [ns]
1818 * \return d^2P/dt/dE [npe/ns/GeV]
1819 */
1820 double getLightFromEMshower(const int type,
1821 const double D_m,
1822 const double cd,
1823 const double theta,
1824 const double phi,
1825 const double t_ns) const
1826 {
1827 switch (type) {
1828
1829 case DIRECT_LIGHT_FROM_EMSHOWER:
1830 return getDirectLightFromEMshower (D_m, cd, theta, phi, t_ns);
1831
1832 case SCATTERED_LIGHT_FROM_EMSHOWER:
1833 return getScatteredLightFromEMshower(D_m, cd, theta, phi, t_ns);
1834
1835 default:
1836 return 0.0;
1837 }
1838 }
1839
1840
1841 /**
1842 * Probability density function for light from EM-shower.
1843 *
1844 * \param D_m distance between EM-shower and PMT [m]
1845 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1846 * \param theta zenith angle orientation PMT [rad]
1847 * \param phi azimuth angle orientation PMT [rad]
1848 * \param t_ns time difference relative to direct Cherenkov light [ns]
1849 * \return d^2P/dt/dE [npe/ns/GeV]
1850 */
1851 double getLightFromEMshower(const double D_m,
1852 const double cd,
1853 const double theta,
1854 const double phi,
1855 const double t_ns) const
1856 {
1857 return (getDirectLightFromEMshower (D_m, cd, theta, phi, t_ns) +
1858 getScatteredLightFromEMshower(D_m, cd, theta, phi, t_ns));
1859 }
1860
1861
1862 /**
1863 * Probability density function for light from EM-shower.
1864 *
1865 * \param type PDF type
1866 * \param E_GeV EM-shower energy [GeV]
1867 * \param D_m distance between EM-shower and PMT [m]
1868 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1869 * \param theta zenith angle orientation PMT [rad]
1870 * \param phi azimuth angle orientation PMT [rad]
1871 * \param t_ns time difference relative to direct Cherenkov light [ns]
1872 * \return d^2P/dt/dE [npe/ns/GeV]
1873 */
1874 double getLightFromEMshower(const int type,
1875 const double E_GeV,
1876 const double D_m,
1877 const double cd,
1878 const double theta,
1879 const double phi,
1880 const double t_ns) const
1881 {
1882 switch (type) {
1883
1884 case DIRECT_LIGHT_FROM_EMSHOWER:
1885 return getDirectLightFromEMshower (E_GeV, D_m, cd, theta, phi, t_ns);
1886
1887 case SCATTERED_LIGHT_FROM_EMSHOWER:
1888 return getScatteredLightFromEMshower(E_GeV, D_m, cd, theta, phi, t_ns);
1889
1890 default:
1891 return 0.0;
1892 }
1893 }
1894
1895
1896 /**
1897 * Probability density function for light from EM-shower.
1898 *
1899 * \param E_GeV EM-shower energy [GeV]
1900 * \param D_m distance between EM-shower and PMT [m]
1901 * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1902 * \param theta zenith angle orientation PMT [rad]
1903 * \param phi azimuth angle orientation PMT [rad]
1904 * \param t_ns time difference relative to direct Cherenkov light [ns]
1905 * \return d^2P/dt/dE [npe/ns/GeV]
1906 */
1907 double getLightFromEMshower(const double E_GeV,
1908 const double D_m,
1909 const double cd,
1910 const double theta,
1911 const double phi,
1912 const double t_ns) const
1913 {
1914 return (getDirectLightFromEMshower (E_GeV, D_m, cd, theta, phi, t_ns) +
1915 getScatteredLightFromEMshower(E_GeV, D_m, cd, theta, phi, t_ns));
1916 }
1917
1918
1919 /**
1920 * Probability density function for direct light from isotropic light source.
1921 *
1922 * \param type PDF type
1923 * \param D_m distance between light source and PMT [m]
1924 * \param ct cosine angle PMT
1925 * \param t_ns time difference relative to direct Cherenkov light [ns]
1926 * \return d^2P/dt/dE [npe/ns/GeV]
1927 */
1928 double getLightFromBrightPoint(const int type,
1929 const double D_m,
1930 const double ct,
1931 const double t_ns) const
1932 {
1933 switch (type) {
1934
1935 case DIRECT_LIGHT_FROM_BRIGHT_POINT:
1936 return getDirectLightFromBrightPoint (D_m, ct, t_ns);
1937
1938 case SCATTERED_LIGHT_FROM_BRIGHT_POINT:
1939 return getScatteredLightFromBrightPoint(D_m, ct, t_ns);
1940
1941 default:
1942 return 0.0;
1943 }
1944 }
1945
1946
1947 /**
1948 * Probability density function for direct light from isotropic light source.
1949 *
1950 * \param D_m distance between light source and PMT [m]
1951 * \param ct cosine angle PMT
1952 * \param t_ns time difference relative to direct Cherenkov light [ns]
1953 * \return d^2P/dt/dE [npe/ns/GeV]
1954 */
1955 double getLightFromBrightPoint(const double D_m,
1956 const double ct,
1957 const double t_ns) const
1958 {
1959 return (getDirectLightFromBrightPoint (D_m, ct, t_ns) +
1960 getScatteredLightFromBrightPoint(D_m, ct, t_ns));
1961 }
1962
1963
1964 /**
1965 * Auxiliary class to find solution(s) to \f$ z \f$ of the square root expression:
1966 \f{eqnarray*}
1967 ct(z=0) & = & z + n \sqrt(z^2 + R^2)
1968 \f}
1969 * where \f$ n = 1/\cos(\theta_{c}) \f$ is the index of refraction.
1970 */
1971 class JRoot {
1972 public:
1973 /**
1974 * Determine solution(s) of quadratic equation.
1975 *
1976 * \param R minimal distance of approach [m]
1977 * \param n index of refraction
1978 * \param t time at z = 0 [ns]
1979 */
1980 JRoot(const double R,
1981 const double n,
1982 const double t) :
1983 is_valid(false),
1984 first (0.0),
1985 second(0.0)
1986 {
1987 const double a = n*n - 1.0;
1988 const double b = 2*C*t;
1989 const double c = R*n * R*n - C*t * C*t;
1990
1991 const double q = b*b - 4*a*c;
1992
1993 if (q >= 0.0) {
1994
1995 first = (-b - sqrt(q)) / (2*a);
1996 second = (-b + sqrt(q)) / (2*a);
1997
1998 is_valid = C*t > second;
1999 }
2000 }
2001
2002
2003 /**
2004 * Get result by index.
2005 *
2006 * \param i index (0 or 1)
2007 * \return i == 0, first; i == 1, second; else throws exception
2008 */
2009 double operator[](const int i) const
2010 {
2011 switch (i) {
2012
2013 case 0:
2014 return first;
2015
2016 case 1:
2017 return second;
2018
2019 default:
2020 THROW(JException, "JRoot::operator[] invalid index");
2021 }
2022 }
2023
2024 bool is_valid; //!< validity of solutions
2025 double first; //!< most upstream solution
2026 double second; //!< most downstream solution
2027 };
2028
2029
2030 protected:
2031 /**
2032 * Determine wavelength for a given index of refraction corresponding to the group velocity.
2033 *
2034 * \param n index of refraction
2035 * \param eps precision index of refraction
2036 * \return wavelength [nm]
2037 */
2038 virtual double getWavelength(const double n,
2039 const double eps = 1.0e-10) const
2040 {
2041 //return getWavelength(n, 0.5*(wmin + wmax), eps);
2042
2043 double vmin = wmin;
2044 double vmax = wmax;
2045
2046 for (int i = 0; i != 1000; ++i) {
2047
2048 const double v = 0.5 * (vmin + vmax);
2049 const double y = getIndexOfRefractionGroup(v);
2050
2051 if (fabs(y - n) < eps) {
2052 return v;
2053 }
2054
2055 if (y < n)
2056 vmax = v;
2057 else
2058 vmin = v;
2059 }
2060
2061 return 0.5 * (vmin + vmax);
2062 }
2063
2064
2065 /**
2066 * Determine wavelength for a given index of refraction corresponding to the group velocity.
2067 * The estimate of the wavelength is made by successive linear extrapolations.
2068 * The procedure starts from the given wavelength and terminates if the index of refraction
2069 * is equal to the target value within the given precision.
2070 *
2071 * \param n index of refraction
2072 * \param w start value wavelength [nm]
2073 * \param eps precision index of refraction
2074 * \return return value wavelength [nm]
2075 */
2076 virtual double getWavelength(const double n,
2077 const double w,
2078 const double eps) const
2079 {
2080 double v = w;
2081
2082 // determine wavelength by linear extrapolation
2083
2084 for ( ; ; ) {
2085
2086 const double y = getIndexOfRefractionGroup(v);
2087
2088 if (fabs(y - n) < eps) {
2089 break;
2090 }
2091
2092 v += (n - y) / getDispersionGroup(v);
2093 }
2094
2095 return v;
2096 }
2097
2098
2099 static double getRmin() { return 1.0e-1; } //!< minimal distance of approach of muon to PMT [m]
2100
2101
2102 /**
2103 * Get the inverse of the attenuation length.
2104 *
2105 * \param l_abs absorption length [m]
2106 * \param ls scattering length [m]
2107 * \param cts cosine scattering angle
2108 * \return inverse attenuation length [m^-1]
2109 */
2110 virtual double getInverseAttenuationLength(const double l_abs, const double ls, const double cts) const
2111 {
2112 static JTOOLS::JGridPolint1Function1D_t f1;
2113
2114 if (f1.empty()) {
2115
2116 const int nx = 100000;
2117 const double xmin = -1.0;
2118 const double xmax = +1.0;
2119
2120 const double dx = (xmax - xmin) / (nx - 1);
2121
2122 for (double x = xmin, W = 0.0; x < xmax; x += dx) {
2123
2124 f1[x] = W;
2125
2126 W += 2*PI * dx * getScatteringProbability(x + 0.5*dx);
2127 }
2128
2129 f1[xmin] = 0.0;
2130 f1[xmax] = 1.0;
2131
2132 f1.compile();
2133 }
2134
2135 return 1.0/l_abs + f1(cts)/ls;
2136 }
2137
2138
2139 protected:
2140 /**
2141 * Integration limits
2142 */
2143 const double wmin; //!< minimal wavelength for integration [nm]
2144 const double wmax; //!< maximal wavelength for integration [nm]
2145
2146 std::vector<element_type> phd; //!< fast evaluation of phi integral
2147 };
2148
2149
2150 /**
2151 * Probability Density Functions of the time response of a PMT
2152 * with an implementation for the JDispersionInterface interface.
2153 */
2154 class JAbstractPDF :
2155 public JPDF,
2156 public JDispersion
2157 {
2158 public:
2159 /**
2160 * Constructor.
2161 *
2162 * \param P_atm ambient pressure [atm]
2163 * \param Wmin minimal wavelength for integration [nm]
2164 * \param Wmax maximal wavelength for integration [nm]
2165 * \param numberOfPoints number of points for integration
2166 * \param epsilon precision of points for integration
2167 */
2168 JAbstractPDF(const double P_atm,
2169 const double Wmin,
2170 const double Wmax,
2171 const int numberOfPoints = 20,
2172 const double epsilon = 1e-12) :
2173 JPDF(Wmin, Wmax, numberOfPoints, epsilon),
2174 JDispersion(P_atm)
2175 {}
2176 };
2177
2178
2179 /**
2180 * Probability Density Functions of the time response of a PMT
2181 * with an implementation of the JAbstractPMT and JAbstractMedium interfaces via C-like methods.
2182 */
2183 class JPDF_C :
2184 public JAbstractPDF
2185 {
2186 public:
2187 /**
2188 * Constructor.
2189 *
2190 * \param Apmt photo-cathode area [m^2]
2191 * \param pF_qe pointer to function for quantum efficiency of PMT
2192 * \param pF_pmt pointer to function for angular acceptance of PMT
2193 * \param pF_l_abs pointer to function for absorption length [m]
2194 * \param pF_ls pointer to function for scattering length [m]
2195 * \param pF_ps pointer to model specific function to describe light scattering in water
2196 * \param P_atm ambient pressure [atm]
2197 * \param Wmin minimal wavelength for integration [nm]
2198 * \param Wmax maximal wavelength for integration [nm]
2199 * \param numberOfPoints number of points for integration
2200 * \param epsilon precision of points for integration
2201 */
2202 JPDF_C(const double Apmt,
2203 //
2204 // pointers to static functions
2205 //
2206 double (*pF_qe) (const double),
2207 double (*pF_pmt) (const double),
2208 double (*pF_l_abs)(const double),
2209 double (*pF_ls) (const double),
2210 double (*pF_ps) (const double),
2211 //
2212 // parameters for base class
2213 //
2214 const double P_atm,
2215 const double Wmin,
2216 const double Wmax,
2217 const int numberOfPoints = 20,
2218 const double epsilon = 1e-12) :
2219 JAbstractPDF(P_atm, Wmin, Wmax, numberOfPoints, epsilon),
2220 A (Apmt),
2221 qe (pF_qe),
2222 l_abs(pF_l_abs),
2223 ls (pF_ls),
2224 pmt (pF_pmt),
2225 ps (pF_ps)
2226 {}
2227
2228
2229 /**
2230 * Photo-cathode area of PMT.
2231 *
2232 * \return photo-cathode area [m^2]
2233 */
2234 virtual double getPhotocathodeArea() const override
2235 {
2236 return A;
2237 }
2238
2239
2240 /**
2241 * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.).
2242 *
2243 * \param lambda wavelenth [nm]
2244 * \return QE
2245 */
2246 virtual double getQE(const double lambda) const override
2247 {
2248 return qe(lambda);
2249 }
2250
2251
2252 /**
2253 * Angular acceptance of PMT.
2254 *
2255 * \param ct cosine angle of incidence
2256 * \return relative efficiency
2257 */
2258 virtual double getAngularAcceptance(const double ct) const override
2259 {
2260 return pmt(ct);
2261 }
2262
2263
2264 /**
2265 * Absorption length.
2266 *
2267 * \param lambda wavelenth [nm]
2268 * \return absorption length [m]
2269 */
2270 virtual double getAbsorptionLength(const double lambda) const override
2271 {
2272 return l_abs(lambda);
2273 }
2274
2275
2276 /**
2277 * Scattering length.
2278 *
2279 * \param lambda wavelenth [nm]
2280 * \return scattering length [m]
2281 */
2282 virtual double getScatteringLength(const double lambda) const override
2283 {
2284 return ls(lambda);
2285 }
2286
2287
2288 /**
2289 * Model specific function to describe light scattering in water
2290 * (integral over full solid angle normalised to one).
2291 *
2292 * \param ct cosine scattering angle
2293 * \return probability
2294 */
2295 virtual double getScatteringProbability(const double ct) const override
2296 {
2297 return ps(ct);
2298 }
2299
2300 protected:
2301 /**
2302 * photo-cathode area [m2]
2303 */
2304 const double A;
2305
2306 /**
2307 * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.)
2308 *
2309 * \param lambda wavelenth [nm]
2310 * \return QE
2311 */
2312 double (*qe)(const double lambda);
2313
2314 /**
2315 * Absorption length
2316 *
2317 * \param lambda wavelenth [nm]
2318 * \return absorption length [m]
2319 */
2320 double (*l_abs)(const double lambda);
2321
2322 /**
2323 * Scattering length
2324 *
2325 * \param lambda wavelenth [nm]
2326 * \return scattering length [m]
2327 */
2328 double (*ls)(const double lambda);
2329
2330 /**
2331 * Angular acceptance of PMT
2332 *
2333 * \param ct cosine angle of incidence
2334 * \return relative efficiency
2335 */
2336 double (*pmt)(const double ct);
2337
2338 /**
2339 * Model specific function to describe light scattering in water
2340 *
2341 * \param ct cosine scattering angle
2342 * \return probability
2343 */
2344 double (*ps)(const double ct);
2345 };
2346}
2347
2348
2349#endif
JCC.hh
Compiler version dependent expressions, macros, etc.
JException.hh
Exceptions.
THROW
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition JException.hh:712
JGeane.hh
Energy loss of muon.
JGeant.hh
Photon emission profile EM-shower.
JGeanz.hh
Longitudinal emission profile EM-shower.
JPDFToolkit.hh
Auxiliary methods for PDF calculations.
JPDFTypes.hh
Numbering scheme for PDF types.
JConstants.hh
Physics constants.
JQuadrature.hh
Auxiliary classes for numerical integration.
numberOfPoints
int numberOfPoints
Definition JResultPDF.cc:22
JLANG::JException
General exception.
Definition JException.hh:24
JPHYSICS::JAbstractMedium::getAbsorptionLength
virtual double getAbsorptionLength(const double lambda) const =0
Absorption length.
JPHYSICS::JAbstractMedium::getScatteringLength
virtual double getScatteringLength(const double lambda) const =0
Scattering length.
JPHYSICS::JAbstractMedium::getScatteringProbability
virtual double getScatteringProbability(const double ct) const =0
Model specific function to describe light scattering in water (integral over full solid angle normali...
JPHYSICS::JAbstractPDF
Probability Density Functions of the time response of a PMT with an implementation for the JDispersio...
Definition JPDF.hh:2157
JPHYSICS::JAbstractPDF::JAbstractPDF
JAbstractPDF(const double P_atm, const double Wmin, const double Wmax, const int numberOfPoints=20, const double epsilon=1e-12)
Constructor.
Definition JPDF.hh:2168
JPHYSICS::JAbstractPMT::getPhotocathodeArea
virtual double getPhotocathodeArea() const =0
Photo-cathode area of PMT.
JPHYSICS::JAbstractPMT::getQE
virtual double getQE(const double lambda) const =0
Quantum efficiency of PMT (incl.
JPHYSICS::JAbstractPMT::getAngularAcceptance
virtual double getAngularAcceptance(const double ct) const =0
Angular acceptence of PMT.
JPHYSICS::JDispersionInterface::getIndexOfRefractionGroup
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
Definition JDispersionInterface.hh:52
JPHYSICS::JDispersionInterface::getDispersionGroup
virtual double getDispersionGroup(const double lambda) const =0
Dispersion of light for group velocity.
JPHYSICS::JDispersionInterface::getDispersionPhase
virtual double getDispersionPhase(const double lambda) const =0
Dispersion of light for phase velocity.
JPHYSICS::JDispersion
Implementation of dispersion for water in deep sea.
Definition JDispersion.hh:28
JPHYSICS::JGeanz::getIntegral
double getIntegral(const double E, const double z) const
Integral of PDF (starting from 0).
Definition JGeanz.hh:95
JPHYSICS::JGeanz::getLength
double getLength(const double E, const double P, const double eps=1.0e-3) const
Get shower length for a given integrated probability.
Definition JGeanz.hh:122
JPHYSICS::JPDF::JRoot
Auxiliary class to find solution(s) to of the square root expression:
Definition JPDF.hh:1971
JPHYSICS::JPDF::JRoot::first
double first
most upstream solution
Definition JPDF.hh:2025
JPHYSICS::JPDF::JRoot::is_valid
bool is_valid
validity of solutions
Definition JPDF.hh:2024
JPHYSICS::JPDF::JRoot::second
double second
most downstream solution
Definition JPDF.hh:2026
JPHYSICS::JPDF::JRoot::operator[]
double operator[](const int i) const
Get result by index.
Definition JPDF.hh:2009
JPHYSICS::JPDF::JRoot::JRoot
JRoot(const double R, const double n, const double t)
Determine solution(s) of quadratic equation.
Definition JPDF.hh:1980
JPHYSICS::JPDF_C
Probability Density Functions of the time response of a PMT with an implementation of the JAbstractPM...
Definition JPDF.hh:2185
JPHYSICS::JPDF_C::JPDF_C
JPDF_C(const double Apmt, double(*pF_qe)(const double), double(*pF_pmt)(const double), double(*pF_l_abs)(const double), double(*pF_ls)(const double), double(*pF_ps)(const double), const double P_atm, const double Wmin, const double Wmax, const int numberOfPoints=20, const double epsilon=1e-12)
Constructor.
Definition JPDF.hh:2202
JPHYSICS::JPDF_C::getScatteringProbability
virtual double getScatteringProbability(const double ct) const override
Model specific function to describe light scattering in water (integral over full solid angle normali...
Definition JPDF.hh:2295
JPHYSICS::JPDF_C::A
const double A
photo-cathode area [m2]
Definition JPDF.hh:2304
JPHYSICS::JPDF_C::getQE
virtual double getQE(const double lambda) const override
Quantum efficiency of PMT (incl.
Definition JPDF.hh:2246
JPHYSICS::JPDF_C::ls
double(* ls)(const double lambda)
Scattering length.
Definition JPDF.hh:2328
JPHYSICS::JPDF_C::pmt
double(* pmt)(const double ct)
Angular acceptance of PMT.
Definition JPDF.hh:2336
JPHYSICS::JPDF_C::qe
double(* qe)(const double lambda)
Quantum efficiency of PMT (incl.
Definition JPDF.hh:2312
JPHYSICS::JPDF_C::getScatteringLength
virtual double getScatteringLength(const double lambda) const override
Scattering length.
Definition JPDF.hh:2282
JPHYSICS::JPDF_C::getPhotocathodeArea
virtual double getPhotocathodeArea() const override
Photo-cathode area of PMT.
Definition JPDF.hh:2234
JPHYSICS::JPDF_C::getAbsorptionLength
virtual double getAbsorptionLength(const double lambda) const override
Absorption length.
Definition JPDF.hh:2270
JPHYSICS::JPDF_C::getAngularAcceptance
virtual double getAngularAcceptance(const double ct) const override
Angular acceptance of PMT.
Definition JPDF.hh:2258
JPHYSICS::JPDF_C::l_abs
double(* l_abs)(const double lambda)
Absorption length.
Definition JPDF.hh:2320
JPHYSICS::JPDF_C::ps
double(* ps)(const double ct)
Model specific function to describe light scattering in water.
Definition JPDF.hh:2344
JPHYSICS::JPDF
Low level interface for the calculation of the Probability Density Functions (PDFs) of the arrival ti...
Definition JPDF.hh:282
JPHYSICS::JPDF::getScatteredLightFromMuon
double getScatteredLightFromMuon(const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for scattered light from muon.
Definition JPDF.hh:878
JPHYSICS::JPDF::getScatteredLightFromEMshower
double getScatteredLightFromEMshower(const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for scattered light from EM-shower.
Definition JPDF.hh:1062
JPHYSICS::JPDF::getLightFromEMshower
double getLightFromEMshower(const double E_GeV, const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for light from EM-shower.
Definition JPDF.hh:1907
JPHYSICS::JPDF::getLightFromBrightPoint
double getLightFromBrightPoint(const double D_m, const double ct, const double t_ns) const
Probability density function for direct light from isotropic light source.
Definition JPDF.hh:1955
JPHYSICS::JPDF::getLightFromEMshower
double getLightFromEMshower(const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for light from EM-shower.
Definition JPDF.hh:1851
JPHYSICS::JPDF::getLightFromMuon
double getLightFromMuon(const int type, const double E_GeV, const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for light from muon.
Definition JPDF.hh:1751
JPHYSICS::JPDF::getLightFromMuon
double getLightFromMuon(const double E_GeV, const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for light from muon.
Definition JPDF.hh:1794
JPHYSICS::JPDF::wmax
const double wmax
maximal wavelength for integration [nm]
Definition JPDF.hh:2144
JPHYSICS::JPDF::getScatteredLightFromBrightPoint
double getScatteredLightFromBrightPoint(const double D_m, const double ct, const double t_ns) const
Probability density function for scattered light from isotropic light source.
Definition JPDF.hh:1651
JPHYSICS::JPDF::getLightFromBrightPoint
double getLightFromBrightPoint(const int type, const double D_m, const double ct, const double t_ns) const
Probability density function for direct light from isotropic light source.
Definition JPDF.hh:1928
JPHYSICS::JPDF::~JPDF
virtual ~JPDF()
Virtual destructor.
Definition JPDF.hh:316
JPHYSICS::JPDF::getWavelength
virtual double getWavelength(const double n, const double w, const double eps) const
Determine wavelength for a given index of refraction corresponding to the group velocity.
Definition JPDF.hh:2076
JPHYSICS::JPDF::getDirectLightFromMuon
double getDirectLightFromMuon(const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for direct light from muon.
Definition JPDF.hh:407
JPHYSICS::JPDF::getDirectLightFromEMshower
double getDirectLightFromEMshower(const double E, const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for direct light from EM-shower.
Definition JPDF.hh:1198
JPHYSICS::JPDF::getLightFromEMshower
double getLightFromEMshower(const int type, const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for light from EM-shower.
Definition JPDF.hh:1820
JPHYSICS::JPDF::getScatteredLightFromDeltaRays
double getScatteredLightFromDeltaRays(const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for scattered light from delta-rays.
Definition JPDF.hh:1461
JPHYSICS::JPDF::element_type
JTOOLS::JElement2D< double, double > element_type
Definition JPDF.hh:284
JPHYSICS::JPDF::getInverseAttenuationLength
virtual double getInverseAttenuationLength(const double l_abs, const double ls, const double cts) const
Get the inverse of the attenuation length.
Definition JPDF.hh:2110
JPHYSICS::JPDF::getRmin
static double getRmin()
minimal distance of approach of muon to PMT [m]
Definition JPDF.hh:2099
JPHYSICS::JPDF::getWavelength
virtual double getWavelength(const double n, const double eps=1.0e-10) const
Determine wavelength for a given index of refraction corresponding to the group velocity.
Definition JPDF.hh:2038
JPHYSICS::JPDF::phd
std::vector< element_type > phd
fast evaluation of phi integral
Definition JPDF.hh:2146
JPHYSICS::JPDF::getScatteredLightFromEMshower
double getScatteredLightFromEMshower(const double E, const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for scattered light from EM-shower.
Definition JPDF.hh:1300
JPHYSICS::JPDF::getDirectLightFromMuon
double getDirectLightFromMuon(const double R_m, const double theta, const double phi) const
Number of photo-electrons from direct Cherenkov light from muon.
Definition JPDF.hh:357
JPHYSICS::JPDF::JPDF
JPDF(const double Wmin, const double Wmax, const int numberOfPoints=20, const double epsilon=1e-12)
Constructor.
Definition JPDF.hh:295
JPHYSICS::JPDF::getNumberOfPhotons
double getNumberOfPhotons() const
Number of Cherenkov photons per unit track length.
Definition JPDF.hh:325
JPHYSICS::JPDF::getDirectLightFromDeltaRays
double getDirectLightFromDeltaRays(const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for direct light from delta-rays.
Definition JPDF.hh:1375
JPHYSICS::JPDF::getDirectLightFromBrightPoint
double getDirectLightFromBrightPoint(const double D_m, const double ct, const double t_ns) const
Probability density function for direct light from isotropic light source.
Definition JPDF.hh:1607
JPHYSICS::JPDF::getScatteredLightFromEMshowers
double getScatteredLightFromEMshowers(const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for scattered light from EM-showers.
Definition JPDF.hh:730
JPHYSICS::JPDF::getLightFromEMshower
double getLightFromEMshower(const int type, const double E_GeV, const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for light from EM-shower.
Definition JPDF.hh:1874
JPHYSICS::JPDF::wmin
const double wmin
Integration limits.
Definition JPDF.hh:2143
JPHYSICS::JPDF::getDirectLightFromEMshower
double getDirectLightFromEMshower(const double D_m, const double cd, const double theta, const double phi, const double t_ns) const
Probability density function for direct light from EM-shower.
Definition JPDF.hh:1006
JPHYSICS::JPDF::getScatteredLightFromMuon
double getScatteredLightFromMuon(const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for scattered light from muon.
Definition JPDF.hh:496
JPHYSICS::JPDF::getDirectLightFromEMshowers
double getDirectLightFromEMshowers(const double R_m, const double theta, const double phi, const double t_ns) const
Probability density function for direct light from EM-showers.
Definition JPDF.hh:642
JTOOLS::JCollection< JElement2D_t >::const_iterator
container_type::const_iterator const_iterator
Definition JCollection.hh:91
JTOOLS::JGaussLegendre
Numerical integrator for .
Definition JQuadrature.hh:113
JPHYSICS
Auxiliary methods for light properties of deep-sea water.
Definition JAbstractMedium.hh:9
JPHYSICS::MODULE_RADIUS_M
static double MODULE_RADIUS_M
Radius of optical module [m].
Definition JPDF.hh:40
JPHYSICS::getDeltaRaysFromMuon
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
JPHYSICS::getDeltaRayProbability
double getDeltaRayProbability(const double x)
Emission profile of photons from delta-rays.
Definition JPDFToolkit.hh:378
JPHYSICS::geanz
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:134
JPHYSICS::cherenkov
double cherenkov(const double lambda, const double n)
Number of Cherenkov photons per unit track length and per unit wavelength.
Definition JPDFToolkit.hh:59
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition JGeane.hh:381
JPHYSICS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:156
JPHYSICS::geant
static const JGeant geant(geanx, 0.0001)
Function object for the number of photons from EM-shower as a function of emission angle.
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition JPDFTypes.hh:38
JPHYSICS::SCATTERED_LIGHT_FROM_BRIGHT_POINT
@ SCATTERED_LIGHT_FROM_BRIGHT_POINT
scattered light from bright point
Definition JPDFTypes.hh:43
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition JPDFTypes.hh:33
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
@ DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition JPDFTypes.hh:29
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWERS
@ SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition JPDFTypes.hh:30
JPHYSICS::DIRECT_LIGHT_FROM_BRIGHT_POINT
@ DIRECT_LIGHT_FROM_BRIGHT_POINT
direct light from bright point
Definition JPDFTypes.hh:42
JPHYSICS::SCATTERED_LIGHT_FROM_MUON
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition JPDFTypes.hh:27
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition JPDFTypes.hh:37
JPHYSICS::DIRECT_LIGHT_FROM_DELTARAYS
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition JPDFTypes.hh:32
JPHYSICS::DIRECT_LIGHT_FROM_MUON
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition JPDFTypes.hh:26
JPHYSICS::geanc
double geanc()
Equivalent muon track length per unit shower energy.
Definition JGeane.hh:28
JPHYSICS::getIndexOfRefractionPhase
double getIndexOfRefractionPhase()
Get average index of refraction of water corresponding to phase velocity.
Definition JPhysics/JConstants.hh:145
JPHYSICS::C
static const double C
Physics constants.
Definition JPhysics/JConstants.hh:25
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JTOOLS
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
Definition JAbstractCollection.hh:11
JTOOLS::JElement2D
2D Element.
Definition JPolint.hh:1131
JTOOLS::JGridPolint1Function1D_t
Type definition of a 1st degree polynomial interpolation based on a JGridCollection with result type ...
Definition JFunction1D_t.hh:294
Generated by doxygen 1.12.0