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