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