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"
12#include "JTools/JQuadrature.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
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 :
279 public virtual JDispersionInterface,
280 public virtual JAbstractPMT,
281 public virtual JAbstractMedium
282 {
283
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) :
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
1761 return getDirectLightFromMuon (R_m, theta, phi, t_ns);
1762
1764 return getScatteredLightFromMuon (R_m, theta, phi, t_ns);
1765
1767 return getDirectLightFromEMshowers (R_m, theta, phi, t_ns) * E_GeV;
1768
1770 return getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV;
1771
1773 return getDirectLightFromDeltaRays (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);
1774
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
1830 return getDirectLightFromEMshower (D_m, cd, theta, phi, t_ns);
1831
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
1885 return getDirectLightFromEMshower (E_GeV, D_m, cd, theta, phi, t_ns);
1886
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
1936 return getDirectLightFromBrightPoint (D_m, ct, t_ns);
1937
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 {
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 */
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
Compiler version dependent expressions, macros, etc.
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Energy loss of muon.
Photon emission profile EM-shower.
Longitudinal emission profile EM-shower.
Auxiliary methods for PDF calculations.
Numbering scheme for PDF types.
Physics constants.
Auxiliary classes for numerical integration.
int numberOfPoints
Definition JResultPDF.cc:22
General exception.
Definition JException.hh:24
virtual double getAbsorptionLength(const double lambda) const =0
Absorption length.
virtual double getScatteringLength(const double lambda) const =0
Scattering length.
virtual double getScatteringProbability(const double ct) const =0
Model specific function to describe light scattering in water (integral over full solid angle normali...
Probability Density Functions of the time response of a PMT with an implementation for the JDispersio...
Definition JPDF.hh:2157
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
virtual double getPhotocathodeArea() const =0
Photo-cathode area of PMT.
virtual double getQE(const double lambda) const =0
Quantum efficiency of PMT (incl.
virtual double getAngularAcceptance(const double ct) const =0
Angular acceptence of PMT.
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
virtual double getDispersionGroup(const double lambda) const =0
Dispersion of light for group velocity.
virtual double getDispersionPhase(const double lambda) const =0
Dispersion of light for phase velocity.
Implementation of dispersion for water in deep sea.
double getIntegral(const double E, const double z) const
Integral of PDF (starting from 0).
Definition JGeanz.hh:95
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
Auxiliary class to find solution(s) to of the square root expression:
Definition JPDF.hh:1971
double first
most upstream solution
Definition JPDF.hh:2025
bool is_valid
validity of solutions
Definition JPDF.hh:2024
double second
most downstream solution
Definition JPDF.hh:2026
double operator[](const int i) const
Get result by index.
Definition JPDF.hh:2009
JRoot(const double R, const double n, const double t)
Determine solution(s) of quadratic equation.
Definition JPDF.hh:1980
Probability Density Functions of the time response of a PMT with an implementation of the JAbstractPM...
Definition JPDF.hh:2185
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
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
const double A
photo-cathode area [m2]
Definition JPDF.hh:2304
virtual double getQE(const double lambda) const override
Quantum efficiency of PMT (incl.
Definition JPDF.hh:2246
double(* ls)(const double lambda)
Scattering length.
Definition JPDF.hh:2328
double(* pmt)(const double ct)
Angular acceptance of PMT.
Definition JPDF.hh:2336
double(* qe)(const double lambda)
Quantum efficiency of PMT (incl.
Definition JPDF.hh:2312
virtual double getScatteringLength(const double lambda) const override
Scattering length.
Definition JPDF.hh:2282
virtual double getPhotocathodeArea() const override
Photo-cathode area of PMT.
Definition JPDF.hh:2234
virtual double getAbsorptionLength(const double lambda) const override
Absorption length.
Definition JPDF.hh:2270
virtual double getAngularAcceptance(const double ct) const override
Angular acceptance of PMT.
Definition JPDF.hh:2258
double(* l_abs)(const double lambda)
Absorption length.
Definition JPDF.hh:2320
double(* ps)(const double ct)
Model specific function to describe light scattering in water.
Definition JPDF.hh:2344
Low level interface for the calculation of the Probability Density Functions (PDFs) of the arrival ti...
Definition JPDF.hh:282
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
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
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
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
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
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
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
const double wmax
maximal wavelength for integration [nm]
Definition JPDF.hh:2144
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
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
virtual ~JPDF()
Virtual destructor.
Definition JPDF.hh:316
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
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
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
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
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
JTOOLS::JElement2D< double, double > element_type
Definition JPDF.hh:284
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
static double getRmin()
minimal distance of approach of muon to PMT [m]
Definition JPDF.hh:2099
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
std::vector< element_type > phd
fast evaluation of phi integral
Definition JPDF.hh:2146
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
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
JPDF(const double Wmin, const double Wmax, const int numberOfPoints=20, const double epsilon=1e-12)
Constructor.
Definition JPDF.hh:295
double getNumberOfPhotons() const
Number of Cherenkov photons per unit track length.
Definition JPDF.hh:325
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
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
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
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
const double wmin
Integration limits.
Definition JPDF.hh:2143
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
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
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
container_type::const_iterator const_iterator
Numerical integrator for .
Auxiliary methods for light properties of deep-sea water.
static double MODULE_RADIUS_M
Radius of optical module [m].
Definition JPDF.hh:40
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.
double getDeltaRayProbability(const double x)
Emission profile of photons from delta-rays.
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
double cherenkov(const double lambda, const double n)
Number of Cherenkov photons per unit track length and per unit wavelength.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition JGeane.hh:381
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
static const JGeant geant(geanx, 0.0001)
Function object for the number of photons from EM-shower as a function of emission angle.
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition JPDFTypes.hh:38
@ SCATTERED_LIGHT_FROM_BRIGHT_POINT
scattered light from bright point
Definition JPDFTypes.hh:43
@ SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition JPDFTypes.hh:33
@ DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition JPDFTypes.hh:29
@ SCATTERED_LIGHT_FROM_EMSHOWERS
scattered light from EM showers
Definition JPDFTypes.hh:30
@ DIRECT_LIGHT_FROM_BRIGHT_POINT
direct light from bright point
Definition JPDFTypes.hh:42
@ SCATTERED_LIGHT_FROM_MUON
scattered light from muon
Definition JPDFTypes.hh:27
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition JPDFTypes.hh:37
@ DIRECT_LIGHT_FROM_DELTARAYS
direct light from delta-rays
Definition JPDFTypes.hh:32
@ DIRECT_LIGHT_FROM_MUON
direct light from muon
Definition JPDFTypes.hh:26
double geanc()
Equivalent muon track length per unit shower energy.
Definition JGeane.hh:28
double getIndexOfRefractionPhase()
Get average index of refraction of water corresponding to phase velocity.
static const double C
Physics constants.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
2D Element.
Definition JPolint.hh:1131
Type definition of a 1st degree polynomial interpolation based on a JGridCollection with result type ...