Jpp
Antares.hh
Go to the documentation of this file.
1 
2 #ifndef __JPHYSICS__ANTARES__
3 #define __JPHYSICS__ANTARES__
4 
5 #include <map>
6 #include <cmath>
7 
8 #include "JTools/JConstants.hh"
9 #include "JTools/JFunction1D_t.hh"
10 
11 
12 /**
13  * \file
14  *
15  * Optical properties of Antares deep-sea site.
16  * \author mdejong
17  */
18 namespace ANTARES {
19 
20  using JTOOLS::PI;
21  using JTOOLS::JPolint1Function1D_t;
22  using JTOOLS::JGridPolint1Function1D_t;
23 
24 
25  /**
26  * Get ambient pressure.
27  *
28  * \return pressure [Atm]
29  */
30  inline double getAmbientPressure()
31  {
32  return 240.0; // ambient pressure [Atm]
33  }
34 
35 
36  /**
37  * Photo-cathode area 10 inch PMT.
38  *
39  * \return photo-cathode area [m^2]
40  */
41  inline const double getPhotocathodeArea()
42  {
43  return 440e-4; // photo-cathode area [m^2]
44  }
45 
46 
47  /**
48  * Absoption length
49  *
50  * \param lambda wavelength of light [nm]
51  * \return absorption length [m]
52  */
53  inline double getAbsorptionLength(const double lambda)
54  {
55  typedef std::map<double, double> map_t;
56 
57  static map_t zmap;
58 
59  if (zmap.empty()) {
60  /**
61  * values taken from:
62  *
63  * afs/in2p3.fr/home/throng/antares/src/km3/v3r6/src/hit-ini_optic.f
64  */
65  //static const double a0 = 0.971;
66  static const double a0 = 1.0;
67 
68  zmap[620.0] = 0.0;
69  zmap[610.0] = a0/0.2890;
70  zmap[600.0] = a0/0.2440;
71  zmap[590.0] = a0/0.1570;
72  zmap[580.0] = a0/0.1080;
73  zmap[570.0] = a0/0.0799;
74  zmap[560.0] = a0/0.0708;
75  zmap[550.0] = a0/0.0638;
76  zmap[540.0] = a0/0.0558;
77  zmap[530.0] = a0/0.0507;
78  zmap[520.0] = a0/0.0477;
79  zmap[510.0] = a0/0.0357;
80  zmap[500.0] = a0/0.0257;
81  zmap[490.0] = a0/0.0196;
82  zmap[480.0] = a0/0.0182;
83  zmap[470.0] = a0/0.0182;
84  zmap[460.0] = a0/0.0191;
85  zmap[450.0] = a0/0.0200;
86  zmap[440.0] = a0/0.0218;
87  zmap[430.0] = a0/0.0237;
88  zmap[420.0] = a0/0.0255;
89  zmap[410.0] = a0/0.0291;
90  zmap[400.0] = a0/0.0325;
91  zmap[390.0] = a0/0.0363;
92  zmap[380.0] = a0/0.0415;
93  zmap[370.0] = a0/0.0473;
94  zmap[360.0] = a0/0.0528;
95  zmap[350.0] = a0/0.0629;
96  zmap[340.0] = a0/0.0710;
97  zmap[330.0] = a0/0.0792;
98  zmap[320.0] = a0/0.0946;
99  zmap[310.0] = a0/0.1090;
100  zmap[300.0] = a0/0.1390;
101  zmap[290.0] = 0.0;
102  }
103 
104  const double x = lambda;
105 
106  if (x > zmap.begin()->first && x < zmap.rbegin()->first) {
107 
108  map_t::const_iterator i = zmap.lower_bound(x);
109  map_t::const_iterator j = i;
110 
111  --j;
112 
113  if (i == zmap.begin()) {
114  ++i;
115  ++j;
116  }
117 
118  const double x1 = i->first;
119  const double x2 = j->first;
120 
121  const double y1 = i->second;
122  const double y2 = j->second;
123 
124  return y1 + (y2 - y1) * (x - x1) / (x2 - x1);
125 
126  } else
127 
128  return 0.0;
129  }
130 
131 
132  /**
133  * Scattering length
134  *
135  * \param lambda wavelength of light [nm]
136  * \return absorption length [m]
137  */
138  inline double getScatteringLength(const double lambda)
139  {
140  typedef std::map<double, double> map_t;
141 
142  static map_t zmap;
143 
144  if (zmap.empty()) {
145 
146  zmap[200.0] = 0.0;
147  zmap[307.0] = 19.7132;
148  zmap[330.0] = 23.8280;
149  zmap[350.0] = 27.6075;
150  zmap[370.0] = 31.5448;
151  zmap[390.0] = 35.6151;
152  zmap[410.0] = 39.7971;
153  zmap[430.0] = 44.0726;
154  zmap[450.0] = 48.4261;
155  zmap[470.0] = 52.8447;
156  zmap[490.0] = 57.3172;
157  zmap[510.0] = 61.8344;
158  zmap[530.0] = 66.3884;
159  zmap[550.0] = 70.9723;
160  zmap[570.0] = 75.5804;
161  zmap[590.0] = 80.2077;
162  zmap[610.0] = 84.8497;
163  zmap[650.0] = 95.0;
164  }
165 
166  const double x = lambda;
167 
168  if (x > zmap.begin()->first && x < zmap.rbegin()->first) {
169 
170  map_t::const_iterator i = zmap.lower_bound(x);
171  map_t::const_iterator j = i;
172 
173  --j;
174 
175  if (i == zmap.begin()) {
176  ++i;
177  ++j;
178  }
179 
180  const double x1 = i->first;
181  const double x2 = j->first;
182 
183  const double y1 = i->second;
184  const double y2 = j->second;
185 
186  return y1 + (y2 - y1) * (x - x1) / (x2 - x1);
187 
188  } else
189 
190  return 0.0;
191  }
192 
193 
194  /**
195  * Auxiliary method to describe light scattering in water (Heneyey-Greenstein)
196  *
197  * \param g angular dependence parameter
198  * \param x cosine scattering angle
199  * \return probability
200  */
201  inline double henyey_greenstein(const double g,
202  const double x)
203 
204  {
205  const double a0 = (1.0 - g*g) / (4*PI);
206  const double y = 1.0 + g*g - 2.0*g*x;
207 
208  return a0 / (y*sqrt(y));
209  }
210 
211 
212  /**
213  * Auxiliary method to describe light scattering in water (Heneyey-Greenstein)
214  *
215  * \param x cosine scattering angle
216  * \return probability
217  */
218  inline double henyey_greenstein(const double x)
219  {
220  static const double g = 0.924;
221 
222  return henyey_greenstein(g, x);
223  }
224 
225 
226  /**
227  * Auxiliary method to describe light scattering in water (Rayleigh)
228  *
229  * \param a angular dependence parameter
230  * \param x cosine scattering angle
231  * \return probability
232  */
233  inline double rayleigh(const double a,
234  const double x)
235 
236  {
237  const double a0 = 1.0 / (1.0 + a/3.0) / (4*PI);
238 
239  return a0 * (1.0 + a*x*x);
240  }
241 
242 
243  /**
244  * Auxiliary method to describe light scattering in water (Rayleigh)
245  *
246  * \param x cosine scattering angle
247  * \return probability
248  */
249  inline double rayleigh(const double x)
250  {
251  return rayleigh(0.853, x);
252  }
253 
254 
255  /**
256  * Model specific function to describe light scattering in water (f4)
257  *
258  * \param x cosine scattering angle
259  * \return probability
260  */
261  inline double f4(const double x)
262  {
263  static const double g1 = 0.77;
264  static const double g2 = 0.00;
265  static const double f = 1.00;
266 
267  return f * henyey_greenstein(g1,x) + (1.0 - f) * henyey_greenstein(g2,x);
268  }
269 
270 
271  /**
272  * Model specific function to describe light scattering in water (p00075)
273  *
274  * \param x cosine scattering angle
275  * \return probability
276  */
277  inline double p00075(const double x)
278  {
279  static const double g = 0.924;
280  static const double f = 0.17;
281 
282  return f * rayleigh(x) + (1.0 - f) * henyey_greenstein(g,x);
283  }
284 
285 
286  /**
287  * Function to describe light scattering in water.
288  *
289  * \param x cosine scattering angle
290  * \return probability
291  */
292  inline double getScatteringProbability(const double x)
293  {
294  return p00075(x);
295  }
296 
297 
298  /**
299  * Angular acceptence of Antares PMT (Genova)
300  *
301  * \param x cosine of angle of incidence
302  * \return probability
303  */
304  inline double genova(const double x)
305  {
306  static const double a0 = 0.3265;
307  static const double a1 = 0.6144;
308  static const double a2 = -0.0343;
309  static const double a3 = -0.0641;
310  static const double a4 = 0.2988;
311  static const double a5 = -0.1422;
312 
313  const double z = -x;
314 
315  double y = 0.0;
316 
317  if (z < -0.65)
318  y = 0.0;
319  else
320  y = a0 + z*(a1 + z*(a2 + z*(a3 + z*(a4 + z*a5))));
321 
322  return y;
323  }
324 
325 
326  /**
327  * Angular acceptence of Antares PMT (Gamelle)
328  *
329  * \param x cosine angle of incidence
330  * \return probability
331  */
332  inline double gamelle(const double x)
333  {
334  static const double a0 = 59.115;
335  static const double a1 = 0.52258;
336  static const double a2 = 0.60944E-02;
337  static const double a3 = -0.16955E-03;
338  static const double a4 = 0.60929E-06;
339 
340  double y = 0.0;
341 
342  if (x > +0.36)
343  y = 0.0;
344  else if (x < -1.00)
345  y = 1.0;
346  else {
347 
348  const double z = acos(-x)*57.29578 + 57.75;
349 
350  y = (a0 + z*(a1 + z*(a2 + z*(a3 + z*a4)))) / 84.0;
351  }
352 
353  return y;
354  }
355 
356 
357  /**
358  * Angular acceptence of Antares PMT
359  *
360  * \param x cosine of angle of incidence
361  * \return probability
362  */
363  inline double getAngularAcceptance(const double x)
364  {
365  return genova(x);
366  }
367 
368 
369  /**
370  * Quantum efficiency of 10-inch Hamamatsu PMT
371  *
372  * \param lambda wavelength of photon [nm]
373  * \param option include absorption in glass and gel (if true, otherwise not)
374  * \return quantum efficiency
375  */
376  inline double getQE(const double lambda, const bool option)
377  {
378  class tuple {
379  public:
380  tuple(const double __QE,
381  const double __l_gel,
382  const double __l_glass) :
383  QE (__QE),
384  l_gel (__l_gel),
385  l_glass(__l_glass)
386  {}
387 
388  double QE; // Quantum efficiiency
389  double l_gel; // gel absorption length [cm]
390  double l_glass; // glass absorption length [cm]
391  };
392 
393  static const tuple ntuple[] = {
394  tuple(0.000e-2, 0.00, 0.00),
395  tuple(1.988e-2, 100.81, 148.37),
396  tuple(2.714e-2, 99.94, 142.87),
397  tuple(3.496e-2, 99.89, 135.64),
398  tuple(4.347e-2, 96.90, 134.58),
399  tuple(5.166e-2, 96.42, 138.27),
400  tuple(6.004e-2, 94.36, 142.40),
401  tuple(6.885e-2, 89.09, 147.16),
402  tuple(8.105e-2, 90.10, 151.80),
403  tuple(10.13e-2, 86.95, 150.88),
404  tuple(13.03e-2, 85.88, 145.68),
405  tuple(15.29e-2, 84.49, 139.70),
406  tuple(16.37e-2, 81.08, 126.55),
407  tuple(17.11e-2, 78.18, 118.86),
408  tuple(17.86e-2, 76.48, 113.90),
409  tuple(18.95e-2, 74.55, 116.08),
410  tuple(20.22e-2, 72.31, 109.23),
411  tuple(21.26e-2, 68.05, 81.63),
412  tuple(22.10e-2, 66.91, 65.66),
413  tuple(22.65e-2, 64.48, 77.30),
414  tuple(23.07e-2, 62.53, 73.02),
415  tuple(23.14e-2, 59.38, 81.25),
416  tuple(23.34e-2, 56.64, 128.04),
417  tuple(22.95e-2, 53.29, 61.84),
418  tuple(22.95e-2, 48.96, 19.23),
419  tuple(22.74e-2, 45.71, 27.21),
420  tuple(23.48e-2, 41.88, 18.09),
421  tuple(22.59e-2, 37.14, 8.41),
422  tuple(20.61e-2, 30.49, 3.92),
423  tuple(17.68e-2, 23.08, 1.82),
424  tuple(13.18e-2, 15.60, 0.84),
425  tuple(7.443e-2, 8.00, 0.39),
426  tuple(2.526e-2, 0.00, 0.17),
427  tuple(0.000e-2, 0.00, 0.00)
428  };
429 
430  static const double cola = 0.9; // collection efficiency
431  static const double x_glass = 1.5; // glass thickness [cm]
432  static const double x_gel = 1.0; // gel thickness [cm]
433 
434  // ntuple
435  static const int N = sizeof(ntuple) / sizeof(ntuple[0]) - 1;
436 
437  static const double xmax = 620.0; // maximal wavelength [nm] (tuple[ 0 ])
438  static const double xmin = 290.0; // minimal wavelength [nm] (tuple[N-1])
439 
440  const double x = lambda;
441 
442  double y = 0.0;
443 
444  if (x > xmin && x < xmax) {
445 
446  const int i = (int) (N * (x - xmax) / (xmin - xmax));
447  const int j = (i == N ? i - 1 : i + 1);
448 
449  const double x1 = xmax + i * (xmin - xmax) / N;
450  const double x2 = xmax + j * (xmin - xmax) / N;
451 
452  const double dx = (x - x1) / (x1 - x2);
453 
454  const double QE = ntuple[i].QE + (ntuple[i].QE - ntuple[j].QE ) * dx;
455  const double l_gel = ntuple[i].l_gel + (ntuple[i].l_gel - ntuple[j].l_gel ) * dx;
456  const double l_glass = ntuple[i].l_glass + (ntuple[i].l_glass - ntuple[j].l_glass) * dx;
457 
458  y = cola * QE;
459 
460  if (option) {
461 
462  if (l_glass > 0.0 && l_gel > 0.0)
463  y *= exp(-x_glass/l_glass) * exp(-x_gel/l_gel);
464  else
465  y = 0.0;
466  }
467  }
468 
469  return y;
470  }
471 
472 
473  /**
474  * Quantum efficiency of 10-inch Hamamatsu PMT
475  *
476  * \param lambda wavelength of photon [nm]
477  * \return quantum efficiency
478  */
479  inline double getQE(const double lambda)
480  {
481  return getQE(lambda, true);
482  }
483 
484 
485  /**
486  * Photo-cathode area 10 inch PMT.
487  *
488  * \param x cosine of angle of incidence
489  * \param lambda wavelength of photon [nm]
490  * \return photo-cathode area [m^2]
491  */
492  inline double getPhotocathodeArea2D(const double x, const double lambda)
493  {
494  return getPhotocathodeArea() * getAngularAcceptance(x);
495  }
496 }
497 
498 #endif
ANTARES::getQE
double getQE(const double lambda, const bool option)
Quantum efficiency of 10-inch Hamamatsu PMT.
Definition: Antares.hh:376
ANTARES::f4
double f4(const double x)
Model specific function to describe light scattering in water (f4)
Definition: Antares.hh:261
ANTARES::gamelle
double gamelle(const double x)
Angular acceptence of Antares PMT (Gamelle)
Definition: Antares.hh:332
g1
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
ANTARES::getScatteringProbability
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:292
ANTARES::getPhotocathodeArea2D
double getPhotocathodeArea2D(const double x, const double lambda)
Photo-cathode area 10 inch PMT.
Definition: Antares.hh:492
JTOOLS::j
int j
Definition: JPolint.hh:634
JFunction1D_t.hh
ANTARES::getAmbientPressure
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:30
JConstants.hh
std::map
Definition: JSTDTypes.hh:16
ANTARES::getAbsorptionLength
double getAbsorptionLength(const double lambda)
Absoption length.
Definition: Antares.hh:53
ANTARES::henyey_greenstein
double henyey_greenstein(const double g, const double x)
Auxiliary method to describe light scattering in water (Heneyey-Greenstein)
Definition: Antares.hh:201
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
JTOOLS::JPolint1Function1D_t
Type definition of a 1st degree polynomial interpolation with result type double.
Definition: JFunction1D_t.hh:185
JTOOLS::JGridPolint1Function1D_t
Type definition of a 1st degree polynomial interpolation based on a JGridCollection with result type ...
Definition: JFunction1D_t.hh:292
ANTARES::genova
double genova(const double x)
Angular acceptence of Antares PMT (Genova)
Definition: Antares.hh:304
ANTARES::rayleigh
double rayleigh(const double a, const double x)
Auxiliary method to describe light scattering in water (Rayleigh)
Definition: Antares.hh:233
ANTARES::p00075
double p00075(const double x)
Model specific function to describe light scattering in water (p00075)
Definition: Antares.hh:277
ANTARES::getScatteringLength
double getScatteringLength(const double lambda)
Scattering length.
Definition: Antares.hh:138
ANTARES::getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of Antares PMT.
Definition: Antares.hh:363
ANTARES
Name space for Antares.
Definition: JDAQClock.hh:293
ANTARES::getPhotocathodeArea
const double getPhotocathodeArea()
Photo-cathode area 10 inch PMT.
Definition: Antares.hh:41
Generated by   doxygen 1.8.16