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