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"
20 #include "JMath/JMathToolkit.hh"
21 #include "JDetector/JModule.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 
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:
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 }
Utility class to parse command line options.
Definition: JParser.hh:1500
data_type w[N+1][M+1]
Definition: JPolint.hh:741
do echo Generating $dir eval D
Definition: JDrawLED.sh:50
int main(int argc, char *argv[])
Definition: Main.cc:15
const double getPhotocathodeArea()
Photo-cathode area 10 inch PMT.
Definition: Antares.hh:41
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
Auxiliary methods for geometrical methods.
Data structure for a composite optical module.
Definition: JModule.hh:57
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:409
do rm f tmp H1
#define STATUS(A)
Definition: JMessage.hh:63
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
Long64_t counter_type
Type definition for counter.
Dynamic ROOT object management.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
#define KM3NET
Definition: Jpp.hh:32
Abstract interface for the generation of points in 3D space.
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
Description of Monte Carlo event generation applications.
Definition: JHead.hh:403
Optical properties of KM3NeT deep-sea site.
Detector specific mapping between logical positions and readout channels of PMTs in optical modules...
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
I/O formatting auxiliaries.
JDirection3D getDirection(const Vec &dir)
Get direction.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
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...
return result
Definition: JPolint.hh:727
JPosition3D getPosition(const Vec &pos)
Get position.
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Physics constants.
static const double PI
Mathematical constants.
int debug
debug level
Definition: JSirene.cc:63
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
double getLength() const
Get length.
Definition: JVector3D.hh:246
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
double getAbsorptionLength(const double lambda)
Absoption length.
Definition: Antares.hh:53
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:45
Auxiliary class to define a range between two values.
then JCalibrateToT a
Definition: JTuneHV.sh:116
Utility class to parse command line options.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
int getCount(const T &hit)
Get hit count.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
data_type v[N+1][M+1]
Definition: JPolint.hh:740
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for photo-cathode area 3 inch PMT.
Definition: KM3NeT.hh:5820
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:35
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:60
Data structure for a composite optical module.
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84