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