Jpp  15.0.0
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 = 1.0/(4*PI); // 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 dp = PI / phd.size();
1566  const double dom = dcb*dp * v*v / (u*u);
1567 
1568  for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
1569 
1570  const double cp = l->getX();
1571  const double sp = l->getY();
1572 
1573  const double vx = sb*cp * qx;
1574  const double vy = sb*sp * qy;
1575  const double vz = cb * qz;
1576 
1577  const double U =
1578  getAngularAcceptance(vx + vy + vz) +
1579  getAngularAcceptance(vx - vy + vz); // PMT angular acceptance
1580 
1581  const double Jb = 1.0/(4*PI); // d^2N/dcos/dphi
1582 
1583  value += dom * npe * geanc() * dz * U * V * W * Ja * Jb * Jc / fabs(Jd);
1584  }
1585  }
1586  }
1587  }
1588  }
1589 
1590  return value;
1591  }
1592 
1593 
1594  /**
1595  * Probability density function for direct light from isotropic light source.
1596  *
1597  * \param D_m distance between light source and PMT [m]
1598  * \param ct cosine angle PMT
1599  * \param t_ns time difference relative to direct Cherenkov light [ns]
1600  * \return d^2P/dt/dE [npe/ns/GeV]
1601  */
1602  double getDirectLightFromBrightPoint(const double D_m,
1603  const double ct,
1604  const double t_ns) const
1605  {
1606  const double D = std::max(D_m, getRmin());
1607  const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1608  const double A = getPhotocathodeArea();
1609 
1610  const double n0 = getIndexOfRefractionGroup(wmax);
1611  const double n1 = getIndexOfRefractionGroup(wmin);
1612  const double ng = C*t/D; // index of refraction
1613 
1614  if (n0 >= ng) { return 0.0; }
1615  if (n1 <= ng) { return 0.0; }
1616 
1617  const double w = getWavelength(ng);
1618  const double n = getIndexOfRefractionPhase(w);
1619 
1620  const double l_abs = getAbsorptionLength(w);
1621  const double ls = getScatteringLength(w);
1622 
1623  const double npe = cherenkov(w,n) * getQE(w);
1624 
1625  const double U = getAngularAcceptance(ct); // PMT angular acceptance
1626  const double V = exp(-D/l_abs - D/ls); // absorption & scattering
1627  const double W = A/(D*D); // solid angle
1628 
1629  const double ngp = getDispersionGroup(w);
1630 
1631  const double Ja = D * ngp / C; // dt/dlambda
1632  const double Jb = 1.0 / (4.0*PI); // d^2N/dOmega
1633 
1634  return npe * geanc() * U * V * W * Jb / fabs(Ja);
1635  }
1636 
1637 
1638  /**
1639  * Probability density function for scattered light from isotropic light source.
1640  *
1641  * \param D_m distance between light source and PMT [m]
1642  * \param ct cosine angle PMT
1643  * \param t_ns time difference relative to direct Cherenkov light [ns]
1644  * \return d^2P/dt/dE [npe/ns/GeV]
1645  */
1646  double getScatteredLightFromBrightPoint(const double D_m,
1647  const double ct,
1648  const double t_ns) const
1649  {
1650  double value = 0;
1651 
1652  const double D = std::max(D_m, getRmin());
1653  const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
1654  const double st = sqrt((1.0 + ct)*(1.0 - ct));
1655 
1656  const double A = getPhotocathodeArea();
1657 
1658  const double n0 = getIndexOfRefractionGroup(wmax);
1659  const double n1 = getIndexOfRefractionGroup(wmin);
1660 
1661  const double ni = C*t / D; // maximal index of refraction
1662 
1663  if (n0 >= ni) {
1664  return value;
1665  }
1666 
1667  const double nj = std::min(ni,n1);
1668 
1669  double w = wmax;
1670 
1671  for (const_iterator i = begin(); i != end(); ++i) {
1672 
1673  const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
1674  const double dn = i->getY() * 0.5 * (nj - n0);
1675 
1676  w = getWavelength(ng, w, 1.0e-5);
1677 
1678  const double dw = dn / fabs(getDispersionGroup(w));
1679 
1680  const double n = getIndexOfRefractionPhase(w);
1681 
1682  const double npe = cherenkov(w,n) * dw * getQE(w);
1683 
1684  if (npe <= 0) { continue; }
1685 
1686  const double l_abs = getAbsorptionLength(w);
1687  const double ls = getScatteringLength(w);
1688 
1689  const double Jc = 1.0 / ls; // dN/dx
1690  const double Jb = 1.0 / (4.0*PI); // d^2N/dcos/dphi
1691 
1692  const double d = C*t / ng; // photon path
1693 
1694  const double dcb = 2.0 / (size() + 1);
1695 
1696  for (double cb = -1.0 + 0.5*dcb; cb < +1.0; cb += dcb) {
1697 
1698  const double sb = sqrt((1.0 + cb)*(1.0 - cb));
1699 
1700  const double v = 0.5 * (d + D) * (d - D) / (d - D*cb);
1701  const double u = d - v;
1702 
1703  if (u <= 0) { continue; }
1704  if (v <= 0) { continue; }
1705 
1706  const double cts = (D*cb - v) / u; // cosine scattering angle
1707 
1708  const double V = exp(-d*getInverseAttenuationLength(l_abs, ls, cts));
1709 
1710  if (cts < 0.0 && v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) { continue; }
1711 
1712  const double W = std::min(A/(v*v), 2.0*PI); // solid angle
1713  const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
1714  const double Jd = ng * (1.0 - cts) / C; // dt/du
1715 
1716  const double dp = PI / phd.size();
1717  const double dom = dcb*dp * v*v / (u*u); // dOmega
1718 
1719  for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
1720 
1721  const double cp = l->getX();
1722  const double dot = cb*ct + sb*cp*st;
1723 
1724  const double U = 2 * getAngularAcceptance(dot); // PMT angular acceptance
1725 
1726  value += npe * geanc() * dom * U * V * W * Ja * Jb * Jc / fabs(Jd);
1727  }
1728  }
1729  }
1730 
1731  return value;
1732  }
1733 
1734 
1735  /**
1736  * Probability density function for light from muon.
1737  *
1738  * \param type PDF type
1739  * \param E_GeV muon energy [GeV]
1740  * \param R_m distance between muon and PMT [m]
1741  * \param theta zenith angle orientation PMT [rad]
1742  * \param phi azimuth angle orientation PMT [rad]
1743  * \param t_ns time difference relative to direct Cherenkov light [ns]
1744  * \return d^2P/dt/dx [npe/ns]
1745  */
1746  double getLightFromMuon(const int type,
1747  const double E_GeV,
1748  const double R_m,
1749  const double theta,
1750  const double phi,
1751  const double t_ns) const
1752  {
1753  switch (type) {
1754 
1755  case DIRECT_LIGHT_FROM_MUON:
1756  return getDirectLightFromMuon (R_m, theta, phi, t_ns);
1757 
1758  case SCATTERED_LIGHT_FROM_MUON:
1759  return getScatteredLightFromMuon (R_m, theta, phi, t_ns);
1760 
1761  case DIRECT_LIGHT_FROM_EMSHOWERS:
1762  return getDirectLightFromEMshowers (R_m, theta, phi, t_ns) * E_GeV;
1763 
1764  case SCATTERED_LIGHT_FROM_EMSHOWERS:
1765  return getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV;
1766 
1767  case DIRECT_LIGHT_FROM_DELTARAYS:
1768  return getDirectLightFromDeltaRays (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);
1769 
1770  case SCATTERED_LIGHT_FROM_DELTARAYS:
1771  return getScatteredLightFromDeltaRays(R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV);
1772 
1773  default:
1774  return 0.0;
1775  }
1776  }
1777 
1778 
1779  /**
1780  * Probability density function for light from muon.
1781  *
1782  * \param E_GeV muon energy [GeV]
1783  * \param R_m distance between muon and PMT [m]
1784  * \param theta zenith angle orientation PMT [rad]
1785  * \param phi azimuth angle orientation PMT [rad]
1786  * \param t_ns time difference relative to direct Cherenkov light [ns]
1787  * \return d^2P/dt/dx [npe/ns]
1788  */
1789  double getLightFromMuon(const double E_GeV,
1790  const double R_m,
1791  const double theta,
1792  const double phi,
1793  const double t_ns) const
1794  {
1795  return (getDirectLightFromMuon (R_m, theta, phi, t_ns) +
1796  getScatteredLightFromMuon (R_m, theta, phi, t_ns) +
1797  getDirectLightFromEMshowers (R_m, theta, phi, t_ns) * E_GeV +
1798  getScatteredLightFromEMshowers(R_m, theta, phi, t_ns) * E_GeV +
1799  getDirectLightFromDeltaRays (R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV) +
1800  getScatteredLightFromDeltaRays(R_m, theta, phi, t_ns) * getDeltaRaysFromMuon(E_GeV));
1801  }
1802 
1803 
1804  /**
1805  * Probability density function for light from EM-shower.
1806  *
1807  * \param type PDF type
1808  * \param D_m distance between EM-shower and PMT [m]
1809  * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1810  * \param theta zenith angle orientation PMT [rad]
1811  * \param phi azimuth angle orientation PMT [rad]
1812  * \param t_ns time difference relative to direct Cherenkov light [ns]
1813  * \return d^2P/dt/dE [npe/ns/GeV]
1814  */
1815  double getLightFromEMshower(const int type,
1816  const double D_m,
1817  const double cd,
1818  const double theta,
1819  const double phi,
1820  const double t_ns) const
1821  {
1822  switch (type) {
1823 
1824  case DIRECT_LIGHT_FROM_EMSHOWER:
1825  return getDirectLightFromEMshower (D_m, cd, theta, phi, t_ns);
1826 
1827  case SCATTERED_LIGHT_FROM_EMSHOWER:
1828  return getScatteredLightFromEMshower(D_m, cd, theta, phi, t_ns);
1829 
1830  default:
1831  return 0.0;
1832  }
1833  }
1834 
1835 
1836  /**
1837  * Probability density function for light from EM-shower.
1838  *
1839  * \param D_m distance between EM-shower and PMT [m]
1840  * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1841  * \param theta zenith angle orientation PMT [rad]
1842  * \param phi azimuth angle orientation PMT [rad]
1843  * \param t_ns time difference relative to direct Cherenkov light [ns]
1844  * \return d^2P/dt/dE [npe/ns/GeV]
1845  */
1846  double getLightFromEMshower(const double D_m,
1847  const double cd,
1848  const double theta,
1849  const double phi,
1850  const double t_ns) const
1851  {
1852  return (getDirectLightFromEMshower (D_m, cd, theta, phi, t_ns) +
1853  getScatteredLightFromEMshower(D_m, cd, theta, phi, t_ns));
1854  }
1855 
1856 
1857  /**
1858  * Probability density function for light from EM-shower.
1859  *
1860  * \param type PDF type
1861  * \param E_GeV EM-shower energy [GeV]
1862  * \param D_m distance between EM-shower and PMT [m]
1863  * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1864  * \param theta zenith angle orientation PMT [rad]
1865  * \param phi azimuth angle orientation PMT [rad]
1866  * \param t_ns time difference relative to direct Cherenkov light [ns]
1867  * \return d^2P/dt/dE [npe/ns/GeV]
1868  */
1869  double getLightFromEMshower(const int type,
1870  const double E_GeV,
1871  const double D_m,
1872  const double cd,
1873  const double theta,
1874  const double phi,
1875  const double t_ns) const
1876  {
1877  switch (type) {
1878 
1879  case DIRECT_LIGHT_FROM_EMSHOWER:
1880  return getDirectLightFromEMshower (E_GeV, D_m, cd, theta, phi, t_ns);
1881 
1882  case SCATTERED_LIGHT_FROM_EMSHOWER:
1883  return getScatteredLightFromEMshower(E_GeV, D_m, cd, theta, phi, t_ns);
1884 
1885  default:
1886  return 0.0;
1887  }
1888  }
1889 
1890 
1891  /**
1892  * Probability density function for light from EM-shower.
1893  *
1894  * \param E_GeV EM-shower energy [GeV]
1895  * \param D_m distance between EM-shower and PMT [m]
1896  * \param cd cosine angle EM-shower direction and EM-shower - PMT position
1897  * \param theta zenith angle orientation PMT [rad]
1898  * \param phi azimuth angle orientation PMT [rad]
1899  * \param t_ns time difference relative to direct Cherenkov light [ns]
1900  * \return d^2P/dt/dE [npe/ns/GeV]
1901  */
1902  double getLightFromEMshower(const double E_GeV,
1903  const double D_m,
1904  const double cd,
1905  const double theta,
1906  const double phi,
1907  const double t_ns) const
1908  {
1909  return (getDirectLightFromEMshower (E_GeV, D_m, cd, theta, phi, t_ns) +
1910  getScatteredLightFromEMshower(E_GeV, D_m, cd, theta, phi, t_ns));
1911  }
1912 
1913 
1914  /**
1915  * Probability density function for direct light from isotropic light source.
1916  *
1917  * \param type PDF type
1918  * \param D_m distance between light source and PMT [m]
1919  * \param ct cosine angle PMT
1920  * \param t_ns time difference relative to direct Cherenkov light [ns]
1921  * \return d^2P/dt/dE [npe/ns/GeV]
1922  */
1923  double getLightFromBrightPoint(const int type,
1924  const double D_m,
1925  const double ct,
1926  const double t_ns) const
1927  {
1928  switch (type) {
1929 
1930  case DIRECT_LIGHT_FROM_BRIGHT_POINT:
1931  return getDirectLightFromBrightPoint (D_m, ct, t_ns);
1932 
1933  case SCATTERED_LIGHT_FROM_BRIGHT_POINT:
1934  return getScatteredLightFromBrightPoint(D_m, ct, t_ns);
1935 
1936  default:
1937  return 0.0;
1938  }
1939  }
1940 
1941 
1942  /**
1943  * Probability density function for direct light from isotropic light source.
1944  *
1945  * \param D_m distance between light source and PMT [m]
1946  * \param ct cosine angle PMT
1947  * \param t_ns time difference relative to direct Cherenkov light [ns]
1948  * \return d^2P/dt/dE [npe/ns/GeV]
1949  */
1950  double getLightFromBrightPoint(const double D_m,
1951  const double ct,
1952  const double t_ns) const
1953  {
1954  return (getDirectLightFromBrightPoint (D_m, ct, t_ns) +
1955  getScatteredLightFromBrightPoint(D_m, ct, t_ns));
1956  }
1957 
1958 
1959  /**
1960  * Auxiliary class to find solution(s) to \f$z\f$ of the square root expression:
1961  \f{eqnarray*}
1962  ct(z=0) & = & z + n \sqrt(z^2 + R^2)
1963  \f}
1964  * where \f$n = 1/\cos(\theta_{c})\f$ is the index of refraction.
1965  */
1966  class JRoot {
1967  public:
1968  /**
1969  * Determine solution(s) of quadratic equation.
1970  *
1971  * \param R minimal distance of approach [m]
1972  * \param n index of refraction
1973  * \param t time at z = 0 [ns]
1974  */
1975  JRoot(const double R,
1976  const double n,
1977  const double t) :
1978  is_valid(false),
1979  first (0.0),
1980  second(0.0)
1981  {
1982  const double a = n*n - 1.0;
1983  const double b = 2*C*t;
1984  const double c = R*n * R*n - C*t * C*t;
1985 
1986  const double q = b*b - 4*a*c;
1987 
1988  if (q >= 0.0) {
1989 
1990  first = (-b - sqrt(q)) / (2*a);
1991  second = (-b + sqrt(q)) / (2*a);
1992 
1993  is_valid = C*t > second;
1994  }
1995  }
1996 
1997 
1998  /**
1999  * Get result by index.
2000  *
2001  * \param i index (0 or 1)
2002  * \return i == 0, first; i == 1, second; else throws exception
2003  */
2004  double operator[](const int i) const
2005  {
2006  switch (i) {
2007 
2008  case 0:
2009  return first;
2010 
2011  case 1:
2012  return second;
2013 
2014  default:
2015  throw JException("JRoot::operator[] invalid index");
2016  }
2017  }
2018 
2019  bool is_valid; //!< validity of solutions
2020  double first; //!< most upstream solution
2021  double second; //!< most downstream solution
2022  };
2023 
2024 
2025  protected:
2026  /**
2027  * Determine wavelength for a given index of refraction corresponding to the group velocity.
2028  *
2029  * \param n index of refraction
2030  * \param eps precision index of refraction
2031  * \return wavelength [nm]
2032  */
2033  virtual double getWavelength(const double n,
2034  const double eps = 1.0e-10) const
2035  {
2036  //return getWavelength(n, 0.5*(wmin + wmax), eps);
2037 
2038  double vmin = wmin;
2039  double vmax = wmax;
2040 
2041  for (int i = 0; i != 1000; ++i) {
2042 
2043  const double v = 0.5 * (vmin + vmax);
2044  const double y = getIndexOfRefractionGroup(v);
2045 
2046  if (fabs(y - n) < eps) {
2047  return v;
2048  }
2049 
2050  if (y < n)
2051  vmax = v;
2052  else
2053  vmin = v;
2054  }
2055 
2056  return 0.5 * (vmin + vmax);
2057  }
2058 
2059 
2060  /**
2061  * Determine wavelength for a given index of refraction corresponding to the group velocity.
2062  * The estimate of the wavelength is made by successive linear extrapolations.
2063  * The procedure starts from the given wavelength and terminates if the index of refraction
2064  * is equal to the target value within the given precision.
2065  *
2066  * \param n index of refraction
2067  * \param w start value wavelength [nm]
2068  * \param eps precision index of refraction
2069  * \return return value wavelength [nm]
2070  */
2071  virtual double getWavelength(const double n,
2072  const double w,
2073  const double eps) const
2074  {
2075  double v = w;
2076 
2077  // determine wavelength by linear extrapolation
2078 
2079  for ( ; ; ) {
2080 
2081  const double y = getIndexOfRefractionGroup(v);
2082 
2083  if (fabs(y - n) < eps) {
2084  break;
2085  }
2086 
2087  v += (n - y) / getDispersionGroup(v);
2088  }
2089 
2090  return v;
2091  }
2092 
2093 
2094  static double getRmin() { return 1.0e-1; } //!< minimal distance of approach of muon to PMT [m]
2095 
2096 
2097  /**
2098  * Get the inverse of the attenuation length.
2099  *
2100  * \param l_abs absorption length [m]
2101  * \param ls scattering length [m]
2102  * \param cts cosine scattering angle
2103  * \return inverse attenuation length [m^-1]
2104  */
2105  virtual double getInverseAttenuationLength(const double l_abs, const double ls, const double cts) const
2106  {
2107  static JTOOLS::JGridPolint1Function1D_t f1;
2108 
2109  if (f1.empty()) {
2110 
2111  const int nx = 100000;
2112  const double xmin = -1.0;
2113  const double xmax = +1.0;
2114 
2115  const double dx = (xmax - xmin) / (nx - 1);
2116 
2117  for (double x = xmin, W = 0.0; x < xmax; x += dx) {
2118 
2119  f1[x] = W;
2120 
2121  W += 2*PI * dx * getScatteringProbability(x + 0.5*dx);
2122  }
2123 
2124  f1[xmin] = 0.0;
2125  f1[xmax] = 1.0;
2126 
2127  f1.compile();
2128  }
2129 
2130  return 1.0/l_abs + f1(cts)/ls;
2131  }
2132 
2133 
2134  protected:
2135  /**
2136  * Integration limits
2137  */
2138  const double wmin; //!< minimal wavelength for integration [nm]
2139  const double wmax; //!< maximal wavelength for integration [nm]
2140 
2141  std::vector<element_type> phd; //!< fast evaluation of phi integral
2142  };
2143 
2144 
2145  /**
2146  * Probability Density Functions of the time response of a PMT
2147  * with an implementation for the JDispersionInterface interface.
2148  */
2149  class JAbstractPDF :
2150  public JPDF,
2151  public JDispersion
2152  {
2153  public:
2154  /**
2155  * Constructor.
2156  *
2157  * \param P_atm ambient pressure [atm]
2158  * \param Wmin minimal wavelength for integration [nm]
2159  * \param Wmax maximal wavelength for integration [nm]
2160  * \param numberOfPoints number of points for integration
2161  * \param epsilon precision of points for integration
2162  */
2163  JAbstractPDF(const double P_atm,
2164  const double Wmin,
2165  const double Wmax,
2166  const int numberOfPoints = 20,
2167  const double epsilon = 1e-12) :
2168  JPDF(Wmin, Wmax, numberOfPoints, epsilon),
2169  JDispersion(P_atm)
2170  {}
2171  };
2172 
2173 
2174  /**
2175  * Probability Density Functions of the time response of a PMT
2176  * with an implementation of the JAbstractPMT and JAbstractMedium interfaces via C-like methods.
2177  */
2178  class JPDF_C :
2179  public JAbstractPDF
2180  {
2181  public:
2182  /**
2183  * Constructor.
2184  *
2185  * \param Apmt photo-cathode area [m^2]
2186  * \param pF_qe pointer to function for quantum efficiency of PMT
2187  * \param pF_pmt pointer to function for angular acceptance of PMT
2188  * \param pF_l_abs pointer to function for absorption length [m]
2189  * \param pF_ls pointer to function for scattering length [m]
2190  * \param pF_ps pointer to model specific function to describe light scattering in water
2191  * \param P_atm ambient pressure [atm]
2192  * \param Wmin minimal wavelength for integration [nm]
2193  * \param Wmax maximal wavelength for integration [nm]
2194  * \param numberOfPoints number of points for integration
2195  * \param epsilon precision of points for integration
2196  */
2197  JPDF_C(const double Apmt,
2198  //
2199  // pointers to static functions
2200  //
2201  double (*pF_qe) (const double),
2202  double (*pF_pmt) (const double),
2203  double (*pF_l_abs)(const double),
2204  double (*pF_ls) (const double),
2205  double (*pF_ps) (const double),
2206  //
2207  // parameters for base class
2208  //
2209  const double P_atm,
2210  const double Wmin,
2211  const double Wmax,
2212  const int numberOfPoints = 20,
2213  const double epsilon = 1e-12) :
2214  JAbstractPDF(P_atm, Wmin, Wmax, numberOfPoints, epsilon),
2215  A (Apmt),
2216  qe (pF_qe),
2217  l_abs(pF_l_abs),
2218  ls (pF_ls),
2219  pmt (pF_pmt),
2220  ps (pF_ps)
2221  {}
2222 
2223 
2224  /**
2225  * Photo-cathode area of PMT.
2226  *
2227  * \return photo-cathode area [m^2]
2228  */
2229  virtual double getPhotocathodeArea() const override
2230  {
2231  return A;
2232  }
2233 
2234 
2235  /**
2236  * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.).
2237  *
2238  * \param lambda wavelenth [nm]
2239  * \return QE
2240  */
2241  virtual double getQE(const double lambda) const override
2242  {
2243  return qe(lambda);
2244  }
2245 
2246 
2247  /**
2248  * Angular acceptance of PMT.
2249  *
2250  * \param ct cosine angle of incidence
2251  * \return relative efficiency
2252  */
2253  virtual double getAngularAcceptance(const double ct) const override
2254  {
2255  return pmt(ct);
2256  }
2257 
2258 
2259  /**
2260  * Absorption length.
2261  *
2262  * \param lambda wavelenth [nm]
2263  * \return absorption length [m]
2264  */
2265  virtual double getAbsorptionLength(const double lambda) const override
2266  {
2267  return l_abs(lambda);
2268  }
2269 
2270 
2271  /**
2272  * Scattering length.
2273  *
2274  * \param lambda wavelenth [nm]
2275  * \return scattering length [m]
2276  */
2277  virtual double getScatteringLength(const double lambda) const override
2278  {
2279  return ls(lambda);
2280  }
2281 
2282 
2283  /**
2284  * Model specific function to describe light scattering in water
2285  * (integral over full solid angle normalised to one).
2286  *
2287  * \param ct cosine scattering angle
2288  * \return probability
2289  */
2290  virtual double getScatteringProbability(const double ct) const override
2291  {
2292  return ps(ct);
2293  }
2294 
2295  protected:
2296  /**
2297  * photo-cathode area [m2]
2298  */
2299  const double A;
2300 
2301  /**
2302  * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.)
2303  *
2304  * \param lambda wavelenth [nm]
2305  * \return QE
2306  */
2307  double (*qe)(const double lambda);
2308 
2309  /**
2310  * Absorption length
2311  *
2312  * \param lambda wavelenth [nm]
2313  * \return absorption length [m]
2314  */
2315  double (*l_abs)(const double lambda);
2316 
2317  /**
2318  * Scattering length
2319  *
2320  * \param lambda wavelenth [nm]
2321  * \return scattering length [m]
2322  */
2323  double (*ls)(const double lambda);
2324 
2325  /**
2326  * Angular acceptance of PMT
2327  *
2328  * \param ct cosine angle of incidence
2329  * \return relative efficiency
2330  */
2331  double (*pmt)(const double ct);
2332 
2333  /**
2334  * Model specific function to describe light scattering in water
2335  *
2336  * \param ct cosine scattering angle
2337  * \return probability
2338  */
2339  double (*ps)(const double ct);
2340  };
2341 }
2342 
2343 
2344 #endif
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:1789
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:741
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 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: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:2323
JPHYSICS::JPDF::JRoot::second
double second
most downstream solution
Definition: JPDF.hh:2021
JPHYSICS::JPDF::getNumberOfPhotons
double getNumberOfPhotons() const
Number of Cherenkov photons per unit track length.
Definition: JPDF.hh:325
JPHYSICS::geanc
double geanc()
Equivalent muon track length per unit shower energy.
Definition: JGeane.hh:26
JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:127
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:2149
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:357
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:1746
JPHYSICS::JPDF::JRoot::JRoot
JRoot(const double R, const double n, const double t)
Determine solution(s) of quadratic equation.
Definition: JPDF.hh:1975
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWERS
direct light from EM showers
Definition: JPDFTypes.hh:32
JPHYSICS::JPDF_C::getAbsorptionLength
virtual double getAbsorptionLength(const double lambda) const override
Absorption length.
Definition: JPDF.hh:2265
N
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
Definition: JShowerPostfit.sh:95
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:1869
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:323
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:2197
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:45
JPHYSICS::JPDF_C::getPhotocathodeArea
virtual double getPhotocathodeArea() const override
Photo-cathode area of PMT.
Definition: JPDF.hh:2229
JPHYSICS::JPDF_C::getQE
virtual double getQE(const double lambda) const override
Quantum efficiency of PMT (incl.
Definition: JPDF.hh:2241
E
then usage E
Definition: JMuonPostfit.sh:35
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...
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"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
Definition: hydrophone-t0:fit.sh:31
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: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:2163
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:2004
JPHYSICS::JPDF::JRoot
Auxiliary class to find solution(s) to of the square root expression: where is the index of refrac...
Definition: JPDF.hh:1966
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:1950
JPHYSICS::C
static const double C
Physics constants.
Definition: JPhysics/JConstants.hh:25
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:1902
JTOOLS::n
const int n
Definition: JPolint.hh:660
JPHYSICS::JPDF::phd
std::vector< element_type > phd
fast evaluation of phi integral
Definition: JPDF.hh:2141
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: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:2331
JPHYSICS::JPDF_C::l_abs
double(* l_abs)(const double lambda)
Absorption length.
Definition: JPDF.hh:2315
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.
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:2178
JPHYSICS::SCATTERED_LIGHT_FROM_DELTARAYS
scattered light from delta-rays
Definition: JPDFTypes.hh:36
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: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:2033
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: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:1375
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:43
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:2019
JPHYSICS::JPDF_C::ps
double(* ps)(const double ct)
Model specific function to describe light scattering in water.
Definition: JPDF.hh:2339
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:2094
JPHYSICS::JPDF::wmax
const double wmax
maximal wavelength for integration [nm]
Definition: JPDF.hh:2139
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
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
d
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
a
then JCalibrateToT a
Definition: JTuneHV.sh:116
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:2071
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:2105
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:1646
numberOfPoints
int numberOfPoints
Definition: JResultPDF.cc:22
JPHYSICS::JPDF::wmin
const double wmin
Integration limits.
Definition: JPDF.hh:2138
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:1602
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::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:149
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:666
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:2290
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:1815
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:1846
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:740
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:739
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:1923
JPHYSICS::JDispersionInterface
Light dispersion inteface.
Definition: JDispersionInterface.hh:19
JPHYSICS::JPDF_C::getScatteringLength
virtual double getScatteringLength(const double lambda) const override
Scattering length.
Definition: JPDF.hh:2277
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
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:2020
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:2307
JPHYSICS::JPDF_C::A
const double A
photo-cathode area [m2]
Definition: JPDF.hh:2299
JPHYSICS::JPDF_C::getAngularAcceptance
virtual double getAngularAcceptance(const double ct) const override
Angular acceptance of PMT.
Definition: JPDF.hh:2253
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