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