Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMultiplicityK40.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <map>
6 #include <cmath>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TRandom3.h"
12 
13 #include "JTools/JRange.hh"
14 #include "JTools/JWeight.hh"
15 #include "JTools/JQuantile.hh"
16 #include "JPhysics/KM3NeT.hh"
17 #include "JPhysics/JConstants.hh"
18 #include "JGeometry3D/JPosition3D.hh"
19 #include "JGeometry3D/JDirection3D.hh"
20 #include "JMath/JMathToolkit.hh"
21 #include "JDetector/JModule.hh"
22 #include "JDetector/JDetectorToolkit.hh"
23 #include "JDetector/JDetectorAddressMapToolkit.hh"
24 #include "JROOT/JManager.hh"
25 #include "JROOT/JRootToolkit.hh"
26 
27 #include "Jeep/JPrint.hh"
28 #include "Jeep/JParser.hh"
29 #include "Jeep/JMessage.hh"
30 
31 namespace {
32 
33  using JTOOLS::PI;
34 
35  typedef JTOOLS::JRange<double> JRange_t;
36 
37  /**
38  * Generation result.
39  */
40  struct JResult {
41  double D; //!< generated distance [m]
42  double V; //!< generation volume [m^3]
43  };
44 
45 
46  /**
47  * Generation interface.
48  */
49  struct JGenerator {
50  virtual const JResult& next() = 0;
51  };
52 
53 
54  /**
55  * Template class for generation of distance between surface of optical module and decay vertex.\n
56  * The template value corresponds to the power for the generation of this distance.
57  */
58  template<int N> struct JEnigma;
59 
60 
61  /**
62  * Generator according D^2.
63  */
64  template<>
65  struct JEnigma<+2> :
66  public JGenerator,
67  private JResult
68  {
69  /*
70  * Constructor.
71  *
72  * \param D_m distance range [m]
73  */
74  JEnigma(const JRange_t& D_m) :
75  D3_min(D_m.getLowerLimit() * D_m.getLowerLimit() * D_m.getLowerLimit()),
76  D3_max(D_m.getUpperLimit() * D_m.getUpperLimit() * D_m.getUpperLimit())
77  {
78  this->V = (4.0/3.0)*PI*(D3_max - D3_min);
79  }
80 
81  virtual const JResult& next()
82  {
83  const double D3 = gRandom->Uniform(D3_min, D3_max);
84 
85  this->D = cbrt(D3);
86 
87  return *this;
88  }
89 
90  const double D3_min;
91  const double D3_max;
92  };
93 
94 
95  /**
96  * Generator according D^0.
97  */
98  template<>
99  struct JEnigma<0> :
100  public JGenerator,
101  private JResult
102  {
103  /*
104  * Constructor.
105  *
106  * \param D_m distance range [m]
107  */
108  JEnigma(const JRange_t& D_m) :
109  D_min(D_m.getLowerLimit()),
110  D_max(D_m.getUpperLimit())
111  {}
112 
113  virtual const JResult& next()
114  {
115  this->D = gRandom->Uniform(D_min, D_max);
116  this->V = 4.0*PI*D*D * (D_max - D_min);
117 
118  return *this;
119  }
120 
121  const double D_min;
122  const double D_max;
123  };
124 
125 
126  /**
127  * Generator according D^-2.
128  */
129  template<>
130  struct JEnigma<-2> :
131  public JGenerator,
132  private JResult
133  {
134  /*
135  * Constructor.
136  *
137  * \param D_m distance range [m]
138  */
139  JEnigma(const JRange_t& D_m) :
140  Dinv_min(1.0/D_m.getUpperLimit()),
141  Dinv_max(1.0/D_m.getLowerLimit())
142  {}
143 
144  virtual const JResult& next()
145  {
146  const double Dinv = gRandom->Uniform(Dinv_min, Dinv_max);
147 
148  this->D = 1.0/Dinv;
149  this->V = 4.0*PI*D*D * D*D * (Dinv_max - Dinv_min);
150 
151  return *this;
152  }
153 
154  const double Dinv_min;
155  const double Dinv_max;
156  };
157 
158 
159  /**
160  * Auxiliary base class for random number generation.
161  */
162  struct JCDF_t :
163  public std::map<double, double>
164  {
165  /**
166  * Get function value.
167  *
168  * \param x random number
169  * \return function value
170  */
171  double operator()(const double x) const
172  {
173  const_iterator i = this->lower_bound(x);
174 
175  if (i != this->begin()) {
176  --i;
177  }
178 
179  return i->second;
180  }
181  };
182 
183 
184  /**
185  * Auxiliary class to generate number of Cherenkov photons from K40 beta decay.\n
186  * Reference: Antares internal note, ANTARES-PHYS-2006-005, J. Brunner.
187  */
188  class JK40BetaDecay :
189  public JCDF_t
190  {
191  public:
192  /**
193  * Default constructor.
194  */
195  JK40BetaDecay()
196  {
197  (*this)[0.000000] = 0.000;
198  (*this)[0.140820] = 1.000;
199  (*this)[0.165515] = 2.000;
200  (*this)[0.185075] = 3.000;
201  (*this)[0.202024] = 4.000;
202  (*this)[0.216980] = 5.000;
203  (*this)[0.231162] = 6.000;
204  (*this)[0.244796] = 7.000;
205  (*this)[0.257734] = 8.000;
206  (*this)[0.270294] = 9.000;
207  (*this)[0.282401] = 10.000;
208  (*this)[0.294095] = 11.000;
209  (*this)[0.305618] = 12.000;
210  (*this)[0.316936] = 13.000;
211  (*this)[0.328064] = 14.000;
212  (*this)[0.338988] = 15.000;
213  (*this)[0.349740] = 16.000;
214  (*this)[0.360311] = 17.000;
215  (*this)[0.370834] = 18.000;
216  (*this)[0.381182] = 19.000;
217  (*this)[0.391350] = 20.000;
218  (*this)[0.401436] = 21.000;
219  (*this)[0.411309] = 22.000;
220  (*this)[0.421253] = 23.000;
221  (*this)[0.430885] = 24.000;
222  (*this)[0.440576] = 25.000;
223  (*this)[0.450050] = 26.000;
224  (*this)[0.459358] = 27.000;
225  (*this)[0.468670] = 28.000;
226  (*this)[0.477844] = 29.000;
227  (*this)[0.487045] = 30.000;
228  (*this)[0.496031] = 31.000;
229  (*this)[0.505032] = 32.000;
230  (*this)[0.513970] = 33.000;
231  (*this)[0.522793] = 34.000;
232  (*this)[0.531627] = 35.000;
233  (*this)[0.540303] = 36.000;
234  (*this)[0.548891] = 37.000;
235  (*this)[0.557355] = 38.000;
236  (*this)[0.565805] = 39.000;
237  (*this)[0.574089] = 40.000;
238  (*this)[0.582449] = 41.000;
239  (*this)[0.590614] = 42.000;
240  (*this)[0.598731] = 43.000;
241  (*this)[0.606772] = 44.000;
242  (*this)[0.614729] = 45.000;
243  (*this)[0.622647] = 46.000;
244  (*this)[0.630523] = 47.000;
245  (*this)[0.638299] = 48.000;
246  (*this)[0.646020] = 49.000;
247  (*this)[0.653625] = 50.000;
248  (*this)[0.661176] = 51.000;
249  (*this)[0.668671] = 52.000;
250  (*this)[0.675974] = 53.000;
251  (*this)[0.683343] = 54.000;
252  (*this)[0.690566] = 55.000;
253  (*this)[0.697753] = 56.000;
254  (*this)[0.704796] = 57.000;
255  (*this)[0.711778] = 58.000;
256  (*this)[0.718749] = 59.000;
257  (*this)[0.725594] = 60.000;
258  (*this)[0.732284] = 61.000;
259  (*this)[0.738976] = 62.000;
260  (*this)[0.745591] = 63.000;
261  (*this)[0.752171] = 64.000;
262  (*this)[0.758597] = 65.000;
263  (*this)[0.764951] = 66.000;
264  (*this)[0.771278] = 67.000;
265  (*this)[0.777509] = 68.000;
266  (*this)[0.783618] = 69.000;
267  (*this)[0.789677] = 70.000;
268  (*this)[0.795685] = 71.000;
269  (*this)[0.801588] = 72.000;
270  (*this)[0.807409] = 73.000;
271  (*this)[0.813077] = 74.000;
272  (*this)[0.818736] = 75.000;
273  (*this)[0.824251] = 76.000;
274  (*this)[0.829650] = 77.000;
275  (*this)[0.835090] = 78.000;
276  (*this)[0.840348] = 79.000;
277  (*this)[0.845611] = 80.000;
278  (*this)[0.850801] = 81.000;
279  (*this)[0.855803] = 82.000;
280  (*this)[0.860776] = 83.000;
281  (*this)[0.865648] = 84.000;
282  (*this)[0.870441] = 85.000;
283  (*this)[0.875074] = 86.000;
284  (*this)[0.879596] = 87.000;
285  (*this)[0.884102] = 88.000;
286  (*this)[0.888459] = 89.000;
287  (*this)[0.892717] = 90.000;
288  (*this)[0.896964] = 91.000;
289  (*this)[0.901105] = 92.000;
290  (*this)[0.905168] = 93.000;
291  (*this)[0.909097] = 94.000;
292  (*this)[0.912907] = 95.000;
293  (*this)[0.916656] = 96.000;
294  (*this)[0.920329] = 97.000;
295  (*this)[0.923924] = 98.000;
296  (*this)[0.927393] = 99.000;
297  (*this)[0.930810] = 100.000;
298  (*this)[0.934111] = 101.000;
299  (*this)[0.937303] = 102.000;
300  (*this)[0.940378] = 103.000;
301  (*this)[0.943379] = 104.000;
302  (*this)[0.946364] = 105.000;
303  (*this)[0.949168] = 106.000;
304  (*this)[0.951929] = 107.000;
305  (*this)[0.954558] = 108.000;
306  (*this)[0.957099] = 109.000;
307  (*this)[0.959577] = 110.000;
308  (*this)[0.961934] = 111.000;
309  (*this)[0.964206] = 112.000;
310  (*this)[0.966378] = 113.000;
311  (*this)[0.968517] = 114.000;
312  (*this)[0.970553] = 115.000;
313  (*this)[0.972528] = 116.000;
314  (*this)[0.974349] = 117.000;
315  (*this)[0.976112] = 118.000;
316  (*this)[0.977841] = 119.000;
317  (*this)[0.979426] = 120.000;
318  (*this)[0.980973] = 121.000;
319  (*this)[0.982446] = 122.000;
320  (*this)[0.983772] = 123.000;
321  (*this)[0.985085] = 124.000;
322  (*this)[0.986341] = 125.000;
323  (*this)[0.987526] = 126.000;
324  (*this)[0.988621] = 127.000;
325  (*this)[0.989622] = 128.000;
326  (*this)[0.990565] = 129.000;
327  (*this)[0.991450] = 130.000;
328  (*this)[0.992296] = 131.000;
329  (*this)[0.993054] = 132.000;
330  (*this)[0.993792] = 133.000;
331  (*this)[0.994464] = 134.000;
332  (*this)[0.995083] = 135.000;
333  (*this)[0.995633] = 136.000;
334  (*this)[0.996158] = 137.000;
335  (*this)[0.996639] = 138.000;
336  (*this)[0.997088] = 139.000;
337  (*this)[0.997470] = 140.000;
338  (*this)[0.997795] = 141.000;
339  (*this)[0.998098] = 142.000;
340  (*this)[0.998372] = 143.000;
341  (*this)[0.998612] = 144.000;
342  (*this)[0.998825] = 145.000;
343  (*this)[0.999005] = 146.000;
344  (*this)[0.999162] = 147.000;
345  (*this)[0.999314] = 148.000;
346  (*this)[0.999441] = 149.000;
347  (*this)[0.999543] = 150.000;
348  (*this)[0.999620] = 151.000;
349  (*this)[0.999695] = 152.000;
350  (*this)[0.999755] = 153.000;
351  (*this)[0.999812] = 154.000;
352  (*this)[0.999847] = 155.000;
353  (*this)[0.999885] = 156.000;
354  (*this)[0.999915] = 157.000;
355  (*this)[0.999939] = 158.000;
356  (*this)[0.999956] = 159.000;
357  (*this)[0.999966] = 160.000;
358  (*this)[0.999975] = 161.000;
359  (*this)[0.999984] = 162.000;
360  (*this)[0.999986] = 163.000;
361  (*this)[0.999992] = 164.000;
362  (*this)[0.999996] = 165.000;
363  (*this)[0.999997] = 166.000;
364  (*this)[0.999998] = 167.000;
365  (*this)[0.999999] = 168.000;
366  (*this)[0.999999] = 169.000;
367  (*this)[0.999999] = 170.000;
368  (*this)[1.000000] = 171.000;
369  /*
370  (*this)[1.000000] = 172.000;
371  (*this)[1.000000] = 173.000;
372  (*this)[1.000000] = 174.000;
373  (*this)[1.000000] = 175.000;
374  (*this)[1.000000] = 176.000;
375  (*this)[1.000000] = 177.000;
376  (*this)[1.000000] = 178.000;
377  (*this)[1.000000] = 179.000;
378  (*this)[1.000000] = 180.000;
379  (*this)[1.000000] = 181.000;
380  (*this)[1.000000] = 182.000;
381  (*this)[1.000000] = 183.000;
382  (*this)[1.000000] = 184.000;
383  (*this)[1.000000] = 185.000;
384  (*this)[1.000000] = 186.000;
385  (*this)[1.000000] = 187.000;
386  (*this)[1.000000] = 188.000;
387  (*this)[1.000000] = 189.000;
388  (*this)[1.000000] = 190.000;
389  (*this)[1.000000] = 191.000;
390  (*this)[1.000000] = 192.000;
391  (*this)[1.000000] = 193.000;
392  (*this)[1.000000] = 194.000;
393  (*this)[1.000000] = 195.000;
394  (*this)[1.000000] = 196.000;
395  (*this)[1.000000] = 197.000;
396  (*this)[1.000000] = 198.000;
397  (*this)[1.000000] = 199.000;
398  (*this)[1.000000] = 200.000;
399  */
400  //compile();
401  }
402 
403  /**
404  * Get branching ratio.
405  *
406  * \return branching ratio
407  */
408  static double getBranchingRatio()
409  {
410  return 0.8928;
411  }
412  };
413 
414 
415  /**
416  * Auxiliary class to generate number of Cherenkov photons from K40 electron capture.\n
417  * Reference: Antares internal note, ANTARES-PHYS-2006-005, J. Brunner.
418  */
419  class JK40ElectronCapture :
420  public JCDF_t
421  {
422  public:
423  /**
424  * Default constructor.
425  */
426  JK40ElectronCapture()
427  {
428  (*this)[0.000000] = 0.000;
429  (*this)[0.000964] = 1.000;
430  (*this)[0.002391] = 2.000;
431  (*this)[0.004139] = 3.000;
432  (*this)[0.005983] = 4.000;
433  (*this)[0.008092] = 5.000;
434  (*this)[0.010555] = 6.000;
435  (*this)[0.013245] = 7.000;
436  (*this)[0.016417] = 8.000;
437  (*this)[0.019601] = 9.000;
438  (*this)[0.023333] = 10.000;
439  (*this)[0.027037] = 11.000;
440  (*this)[0.031226] = 12.000;
441  (*this)[0.035502] = 13.000;
442  (*this)[0.040167] = 14.000;
443  (*this)[0.045124] = 15.000;
444  (*this)[0.050383] = 16.000;
445  (*this)[0.055996] = 17.000;
446  (*this)[0.061656] = 18.000;
447  (*this)[0.067462] = 19.000;
448  (*this)[0.073877] = 20.000;
449  (*this)[0.080352] = 21.000;
450  (*this)[0.086917] = 22.000;
451  (*this)[0.093883] = 23.000;
452  (*this)[0.100961] = 24.000;
453  (*this)[0.108294] = 25.000;
454  (*this)[0.116022] = 26.000;
455  (*this)[0.123821] = 27.000;
456  (*this)[0.131823] = 28.000;
457  (*this)[0.139890] = 29.000;
458  (*this)[0.147991] = 30.000;
459  (*this)[0.156341] = 31.000;
460  (*this)[0.164641] = 32.000;
461  (*this)[0.173283] = 33.000;
462  (*this)[0.181640] = 34.000;
463  (*this)[0.189983] = 35.000;
464  (*this)[0.198629] = 36.000;
465  (*this)[0.207442] = 37.000;
466  (*this)[0.216109] = 38.000;
467  (*this)[0.225215] = 39.000;
468  (*this)[0.234364] = 40.000;
469  (*this)[0.243404] = 41.000;
470  (*this)[0.252236] = 42.000;
471  (*this)[0.261426] = 43.000;
472  (*this)[0.270904] = 44.000;
473  (*this)[0.280399] = 45.000;
474  (*this)[0.290229] = 46.000;
475  (*this)[0.300124] = 47.000;
476  (*this)[0.309784] = 48.000;
477  (*this)[0.319648] = 49.000;
478  (*this)[0.329575] = 50.000;
479  (*this)[0.339769] = 51.000;
480  (*this)[0.350323] = 52.000;
481  (*this)[0.360685] = 53.000;
482  (*this)[0.371175] = 54.000;
483  (*this)[0.381689] = 55.000;
484  (*this)[0.392278] = 56.000;
485  (*this)[0.402836] = 57.000;
486  (*this)[0.413363] = 58.000;
487  (*this)[0.424129] = 59.000;
488  (*this)[0.435091] = 60.000;
489  (*this)[0.445833] = 61.000;
490  (*this)[0.456412] = 62.000;
491  (*this)[0.466519] = 63.000;
492  (*this)[0.477130] = 64.000;
493  (*this)[0.487492] = 65.000;
494  (*this)[0.497307] = 66.000;
495  (*this)[0.507566] = 67.000;
496  (*this)[0.517633] = 68.000;
497  (*this)[0.527133] = 69.000;
498  (*this)[0.536413] = 70.000;
499  (*this)[0.545223] = 71.000;
500  (*this)[0.554018] = 72.000;
501  (*this)[0.562511] = 73.000;
502  (*this)[0.570839] = 74.000;
503  (*this)[0.579232] = 75.000;
504  (*this)[0.587318] = 76.000;
505  (*this)[0.595435] = 77.000;
506  (*this)[0.603502] = 78.000;
507  (*this)[0.611422] = 79.000;
508  (*this)[0.619045] = 80.000;
509  (*this)[0.626384] = 81.000;
510  (*this)[0.633823] = 82.000;
511  (*this)[0.641268] = 83.000;
512  (*this)[0.648831] = 84.000;
513  (*this)[0.656397] = 85.000;
514  (*this)[0.663693] = 86.000;
515  (*this)[0.671029] = 87.000;
516  (*this)[0.678402] = 88.000;
517  (*this)[0.685922] = 89.000;
518  (*this)[0.693255] = 90.000;
519  (*this)[0.700336] = 91.000;
520  (*this)[0.707653] = 92.000;
521  (*this)[0.714999] = 93.000;
522  (*this)[0.721974] = 94.000;
523  (*this)[0.728990] = 95.000;
524  (*this)[0.736015] = 96.000;
525  (*this)[0.742894] = 97.000;
526  (*this)[0.750246] = 98.000;
527  (*this)[0.757448] = 99.000;
528  (*this)[0.764563] = 100.000;
529  (*this)[0.771738] = 101.000;
530  (*this)[0.778704] = 102.000;
531  (*this)[0.785757] = 103.000;
532  (*this)[0.793025] = 104.000;
533  (*this)[0.800100] = 105.000;
534  (*this)[0.807125] = 106.000;
535  (*this)[0.814274] = 107.000;
536  (*this)[0.821156] = 108.000;
537  (*this)[0.828160] = 109.000;
538  (*this)[0.834846] = 110.000;
539  (*this)[0.841731] = 111.000;
540  (*this)[0.848563] = 112.000;
541  (*this)[0.855346] = 113.000;
542  (*this)[0.862256] = 114.000;
543  (*this)[0.868982] = 115.000;
544  (*this)[0.875899] = 116.000;
545  (*this)[0.882461] = 117.000;
546  (*this)[0.888889] = 118.000;
547  (*this)[0.895478] = 119.000;
548  (*this)[0.901776] = 120.000;
549  (*this)[0.908026] = 121.000;
550  (*this)[0.914094] = 122.000;
551  (*this)[0.920233] = 123.000;
552  (*this)[0.926076] = 124.000;
553  (*this)[0.931717] = 125.000;
554  (*this)[0.937147] = 126.000;
555  (*this)[0.942499] = 127.000;
556  (*this)[0.947630] = 128.000;
557  (*this)[0.952460] = 129.000;
558  (*this)[0.956957] = 130.000;
559  (*this)[0.961478] = 131.000;
560  (*this)[0.965667] = 132.000;
561  (*this)[0.969667] = 133.000;
562  (*this)[0.973330] = 134.000;
563  (*this)[0.976881] = 135.000;
564  (*this)[0.980044] = 136.000;
565  (*this)[0.982943] = 137.000;
566  (*this)[0.985614] = 138.000;
567  (*this)[0.987847] = 139.000;
568  (*this)[0.990126] = 140.000;
569  (*this)[0.991874] = 141.000;
570  (*this)[0.993441] = 142.000;
571  (*this)[0.994695] = 143.000;
572  (*this)[0.995898] = 144.000;
573  (*this)[0.996831] = 145.000;
574  (*this)[0.997633] = 146.000;
575  (*this)[0.998305] = 147.000;
576  (*this)[0.998762] = 148.000;
577  (*this)[0.999114] = 149.000;
578  (*this)[0.999362] = 150.000;
579  (*this)[0.999534] = 151.000;
580  (*this)[0.999705] = 152.000;
581  (*this)[0.999801] = 153.000;
582  (*this)[0.999876] = 154.000;
583  (*this)[0.999919] = 155.000;
584  (*this)[0.999953] = 156.000;
585  (*this)[0.999966] = 157.000;
586  (*this)[0.999972] = 158.000;
587  (*this)[0.999978] = 159.000;
588  (*this)[0.999978] = 160.000;
589  (*this)[0.999984] = 161.000;
590  (*this)[0.999988] = 162.000;
591  (*this)[0.999988] = 163.000;
592  (*this)[0.999988] = 164.000;
593  (*this)[0.999988] = 165.000;
594  (*this)[0.999988] = 166.000;
595  (*this)[0.999988] = 167.000;
596  (*this)[0.999988] = 168.000;
597  (*this)[0.999994] = 169.000;
598  (*this)[0.999994] = 170.000;
599  (*this)[0.999994] = 171.000;
600  (*this)[0.999994] = 172.000;
601  (*this)[0.999994] = 173.000;
602  (*this)[0.999994] = 174.000;
603  (*this)[0.999994] = 175.000;
604  (*this)[0.999994] = 176.000;
605  (*this)[0.999997] = 177.000;
606  (*this)[0.999997] = 178.000;
607  (*this)[0.999997] = 179.000;
608  (*this)[0.999997] = 180.000;
609  (*this)[1.000000] = 181.000;
610  /*
611  (*this)[1.000000] = 182.000;
612  (*this)[1.000000] = 183.000;
613  (*this)[1.000000] = 184.000;
614  (*this)[1.000000] = 185.000;
615  (*this)[1.000000] = 186.000;
616  (*this)[1.000000] = 187.000;
617  (*this)[1.000000] = 188.000;
618  (*this)[1.000000] = 189.000;
619  (*this)[1.000000] = 190.000;
620  (*this)[1.000000] = 191.000;
621  (*this)[1.000000] = 192.000;
622  (*this)[1.000000] = 193.000;
623  (*this)[1.000000] = 194.000;
624  (*this)[1.000000] = 195.000;
625  (*this)[1.000000] = 196.000;
626  (*this)[1.000000] = 197.000;
627  (*this)[1.000000] = 198.000;
628  (*this)[1.000000] = 199.000;
629  (*this)[1.000000] = 200.000;
630  */
631  //compile();
632  }
633 
634  /**
635  * Get branching ratio.
636  *
637  * \return branching ratio
638  */
639  static double getBranchingRatio()
640  {
641  return 0.1072;
642  }
643  };
644 
645  /**
646  * Function objects for K40 beta.
647  */
648  static const JK40BetaDecay k40_beta_decay;
649  static const JK40ElectronCapture k40_electron_capture;
650 
651  /**
652  * Parametrised angular acceptence of KM3NeT PMT.\n
653  * Reference: Antares internal note, ANTARES-PHYS-2006-005, J. Brunner.
654  */
655  double a;
656 
657  /**
658  * Parametrised angular acceptence of KM3NeT PMT.
659  *
660  * \param ct cosine of angle of incidence
661  * \return probability
662  */
663  inline double get_angular_acceptance(const double ct)
664  {
665  const double w = 0.25*a*(1.0 + ct)*(1.0 + ct) - ct;
666 
667  if (w > 0.0)
668  return 1.4 * w;
669  else
670  return 0.0;
671  }
672 }
673 
674 
675 /**
676  * \file
677  *
678  * Example program to calculate multiples rate.
679  * \author mdejong
680  */
681 int main(int argc, char* argv[])
682 {
683  using namespace std;
684  using namespace JPP;
685 
686  typedef long long int counter_type;
687 
688  string outputFile;
689  double bequerel;
690  JRange_t D_m;
691  counter_type numberOfEvents;
692  double QE;
693  UInt_t seed;
694  int generator;
695  double focus;
696  int Aefftype;
697  bool exclusive;
698  int debug;
699 
700  try {
701 
702  JParser<> zap("Example program to calculate multiples rate.");
703 
704  zap['o'] = make_field(outputFile) = "k40.root";
705  zap['b'] = make_field(bequerel) = 13.75e3; // [m^-3 s^-1]
706  zap['n'] = make_field(numberOfEvents) = 1e7;
707  zap['D'] = make_field(D_m) = JRange_t(0.216, 10);
708  zap['Q'] = make_field(QE) = 1.0;
709  zap['S'] = make_field(seed) = 0;
710  zap['G'] = make_field(generator) = 0, +2, -2;
711  zap['F'] = make_field(focus) = 1.0;
712  zap['A'] = make_field(Aefftype) = 1, 2, 3;
713  zap['U'] = make_field(exclusive);
714  zap['a'] = make_field(a) = 0.0;
715  zap['d'] = make_field(debug) = 2;
716 
717  zap(argc, argv);
718  }
719  catch(const exception &error) {
720  FATAL(error.what() << endl);
721  }
722 
723  gRandom->SetSeed(seed);
724 
725  using namespace KM3NET;
726  using namespace JPP;
727 
728  const int id = 1;
729  const JModule module = getModule(getDetectorAddressMap(14).get(id), id);
730 
731  DEBUG(module << endl);
732 
733  bequerel /= focus; // correct decay rate for focussing of photons
734 
735  const double R_m = 17.0 * 2.54 * 0.5e-2; // radius 17'' optical module [m]
736  const double A = PI * R_m * R_m; // cross section optical module [m^2]
737 
738  const double wmin = 280.0; // minimal wavelength [nm]
739  const double wmax = 700.0; // maximal wavelength [nm]
740  double ng = 41.0; // average number of photons per decay in given wavelength range
741 
742  const double WAVELENGTH_EXPANSION = (wmax-wmin) / (wmin*wmax) * (300.0*600.0)/(600.0-300.0);
743 
744  JGenerator* enigma = NULL; // distance generation
745 
746  switch (generator) {
747 
748  case +2: enigma = new JEnigma<+2>(D_m); break;
749  case 0: enigma = new JEnigma< 0>(D_m); break;
750  case -2: enigma = new JEnigma<-2>(D_m); break;
751 
752  default: FATAL("No generator type " << generator << endl);
753  }
754 
755  const double vmin = 1.0 / wmax; // wavelength generation
756  const double vmax = 1.0 / wmin; // as lambda^-2
757 
758  double QEmax = 0.0;
759 
760  for (double w = wmin; w <= wmax; w += 1.0) {
761  if (getQE(w) > QEmax) {
762  QEmax = getQE(w);
763  }
764  }
765 
766  NOTICE("Maximal QE " << FIXED(5,3) << QEmax << endl);
767  NOTICE("Wavelength expansion " << FIXED(5,3) << WAVELENGTH_EXPANSION << endl);
768  NOTICE("Number of photons per decay " << FIXED(5,2) << ng << endl);
769 
770  typedef JManager<int, TH1D> JManager_t;
771 
772  JManager_t H1(new TH1D("M[%]", NULL, 100, D_m.getLowerLimit(), D_m.getUpperLimit()));
773 
774  H1->Sumw2();
775 
776  TH1D pmt("pmt", NULL, 1000, -1.0, +1.0);
777 
778  for (Int_t i = 1; i != pmt.GetNbinsX(); ++i) {
779 
780  const double dot = pmt.GetBinCenter(i);
781 
782  double y = 0.0;
783 
784  switch (Aefftype) {
785 
786  case 1:
787  y = getAngularAcceptance(dot);
788  break;
789 
790  case 3:
791  y = get_angular_acceptance(dot);
792  break;
793  }
794 
795  pmt.SetBinContent(i, y);
796  }
797 
798 
799  map<int, JWeight> h1;
800  map<int, JQuantile> P2;
801 
802 
803  vector<double> pi(module.size());
804  vector<int> buffer;
805 
806  for (counter_type event_count = 0; event_count != numberOfEvents; ++event_count) {
807 
808  if (event_count%10000 == 0) {
809  STATUS("event: " << setw(10) << event_count << '\r'); DEBUG(endl);
810  }
811 
812  const JResult& result = enigma->next();
813 
814  const double D = result.D;
815  const double V = result.V;
816 
817  // check first photons on cross section of optical module at its face
818  // it is assumed that the light from the K40 decay is purely isotropic
819 
820  double W = A / (4*PI*(D-R_m)*(D-R_m));
821 
822  if (W > 0.5) {
823  W = 0.5;
824  }
825 
826  //
827  double x = gRandom->Rndm(); // decay mode
828  double y = 0.0; // number of photons
829 
830  if ((x -= k40_beta_decay .getBranchingRatio()) <= 0.0)
831  y = k40_beta_decay (gRandom->Rndm());
832  else if ((x -= k40_electron_capture.getBranchingRatio()) <= 0.0)
833  y = k40_electron_capture(gRandom->Rndm());
834 
835  const int N = gRandom->Poisson(y * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
836  //
837  /*
838  const int N = gRandom->Poisson(ng * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
839  */
840  if (N >= 2) {
841 
842  // generate vertex
843 
844  const double ct = gRandom->Uniform(-1.0, +1.0);
845  const double phi = gRandom->Uniform(-PI, +PI);
846 
847  const double st = sqrt((1.0 - ct) * (1.0 + ct));
848 
849  const JPosition3D vertex(D * st * cos(phi),
850  D * st * sin(phi),
851  D * ct);
852 
853  buffer.clear();
854 
855  for (int i = 0; i != N; ++i) {
856 
857  // generate wavelength of photon
858 
859  const double v = gRandom->Uniform(vmin, vmax);
860  const double w = 1.0 / v;
861 
862  const double l_abs = getAbsorptionLength(w);
863 
864  double P = 0.0;
865 
866  for (size_t pmt = 0; pmt != module.size(); ++pmt) {
867 
868  JPosition3D pos(module[pmt].getPosition());
869 
870  pos -= vertex;
871 
872  const double d = pos.getLength();
873 
874  pos /= d; // direction of photon
875 
876  if (d < D - R_m) {
877  ERROR("Distance " << d << " < " << D << endl);
878  }
879 
880  const double dot = getDot(pos, module[pmt].getDirection());
881 
882  const double U = 0.5 * (1.0 - d / sqrt(d*d + getPhotocathodeArea() / PI));
883 
884  double p = 0.0;
885 
886  switch (Aefftype) {
887 
888  case 1:
889  p = getAngularAcceptance(dot) * getQE(w);
890  break;
891 
892  case 2:
893  p = getPhotocathodeArea2D(dot, w) / getPhotocathodeArea();
894  break;
895 
896  case 3:
897  p = get_angular_acceptance(dot) * getQE(w);
898  break;
899  }
900 
901  P += pi[pmt] = U * p * exp(-d/l_abs);
902  }
903 
904  if (P > W) {
905  ERROR("Probability " << P << " > " << W << endl);
906  }
907 
908  if (W * QEmax * gRandom->Rndm() < P) {
909 
910  int pmt = 0;
911  double y = gRandom->Uniform(P);
912 
913  for (vector<double>::const_iterator i = pi.begin(); i != pi.end() && (y -= *i) > 0.0; ++i, ++pmt) {}
914 
915  buffer.push_back(pmt);
916  }
917  }
918 
919  if (!buffer.empty()) {
920 
921  int M = buffer.size();
922 
923  if (exclusive) {
924 
925  sort(buffer.begin(), buffer.end());
926 
927  M = distance(buffer.begin(), unique(buffer.begin(), buffer.end()));
928  }
929 
930  h1[M].put(V);
931  H1[M]->Fill(D, V);
932 
933  for (int i = 2; i <= M; ++i) {
934  P2[i].put((double) (buffer.size() - M) / (double) M, V);
935  }
936 
937  }
938  }
939  }
940  STATUS(endl);
941 
942  for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
943  i->second->Scale(bequerel / (double) numberOfEvents);
944  }
945 
946  for (size_t M = 2; M != 7; ++M) {
947  cout << "Rate[" << M << "] = "
948  << FIXED(7,3) << bequerel * h1[M].getTotal() / (double) numberOfEvents
949  << " +/- "
950  << FIXED(7,3) << bequerel * h1[M].getError() / (double) numberOfEvents
951  << " Hz" << endl;
952  }
953 
954  for (size_t M = 2; M != 7; ++M) {
955  if (P2[M].getCount() != 0) {
956  cout << "P2[" << M << "] = " << P2[M].getMean() << endl;
957  }
958  }
959 
960  TFile out(outputFile.c_str(), "recreate");
961 
962  out << H1 << pmt;
963 
964  out.Write();
965  out.Close();
966 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:741
D
do echo Generating $dir eval D
Definition: JDrawLED.sh:50
main
int main(int argc, char *argv[])
Definition: Main.cc:15
ANTARES::getPhotocathodeArea
const double getPhotocathodeArea()
Photo-cathode area 10 inch PMT.
Definition: Antares.hh:41
M
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
JMathToolkit.hh
Auxiliary methods for geometrical methods.
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:57
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JASTRONOMY::getDot
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:409
H1
do rm f tmp H1
Definition: plot-equalizer.sh:38
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
JROOT::counter_type
Long64_t counter_type
Type definition for counter.
Definition: JROOT/JCounter.hh:24
JManager.hh
Dynamic ROOT object management.
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
KM3NET
#define KM3NET
Definition: Jpp.hh:32
JMARKOV::JGenerator
Abstract interface for the generation of points in 3D space.
Definition: JMarkovGenerator.hh:35
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:28
JAANET::generator
Description of Monte Carlo event generation applications.
Definition: JHead.hh:403
KM3NeT.hh
Optical properties of KM3NeT deep-sea site.
JDetectorAddressMapToolkit.hh
Detector specific mapping between logical positions and readout channels of PMTs in optical modules...
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
JPrint.hh
I/O formatting auxiliaries.
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition: JAAnetToolkit.hh:193
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
JDETECTOR::getQE
double getQE(const double R, const double mu)
Get QE for given ratio of hit probabilities and expectation value of the number of photo-electrons...
Definition: JPMTParametersToolkit.hh:89
JTOOLS::result
return result
Definition: JPolint.hh:727
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition: JAAnetToolkit.hh:157
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
JConstants.hh
Physics constants.
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
debug
int debug
debug level
Definition: JSirene.cc:63
JDETECTOR::getDetectorAddressMap
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
Definition: JDetectorAddressMapToolkit.hh:363
JGEOMETRY3D::JVector3D::getLength
double getLength() const
Get length.
Definition: JVector3D.hh:246
JMessage.hh
General purpose messaging.
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
ANTARES::getAbsorptionLength
double getAbsorptionLength(const double lambda)
Absoption length.
Definition: Antares.hh:53
d
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:45
JRange.hh
Auxiliary class to define a range between two values.
a
then JCalibrateToT a
Definition: JTuneHV.sh:116
JParser.hh
Utility class to parse command line options.
DEBUG
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Definition: JCalibratePMT.sh:21
JTRIGGER::getCount
int getCount(const T &hit)
Get hit count.
Definition: JTrigger/JHit.hh:277
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:740
KM3NET::getPhotocathodeArea2D
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for photo-cathode area 3 inch PMT.
Definition: KM3NeT.hh:5820
N
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:35
A
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
Definition: JShellParser.csh:15
getModule
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
Definition: JBeaconSimulator.cc:928
P
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:60
JModule.hh
Data structure for a composite optical module.
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84
Generated by   doxygen 1.8.5