Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
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"
19 #include "JGeometry3D/JPosition3D.hh"
20 #include "JGeometry3D/JDirection3D.hh"
21 #include "JMath/JMathToolkit.hh"
22 #include "JDetector/JModule.hh"
23 #include "JDetector/JDetectorToolkit.hh"
24 #include "JDetector/JDetectorSupportkit.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<JKM3NeT_t>(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 
800  map<int, JWeight> h1;
801  map<int, JQuantile> P2;
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:
894  p = KM3NET2D::getPhotocathodeArea2D(dot, w) / getPhotocathodeArea();
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 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JDetectorSupportkit.hh
Detector support kit.
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
JManager.hh
Dynamic ROOT object management.
JMathToolkit.hh
Auxiliary methods for geometrical methods.
JMessage.hh
General purpose messaging.
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
JModule.hh
Data structure for optical module.
main
int main(int argc, char *argv[])
Definition: JMultiplicityK40.cc:682
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
JConstants.hh
Physics constants.
JPrint.hh
I/O formatting auxiliaries.
JRange.hh
Auxiliary class to define a range between two values.
KM3NeT2D.hh
Properties of KM3NeT PMT and deep-sea water.
KM3NeT.hh
Properties of KM3NeT PMT and deep-sea water.
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
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:75
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
JGEOMETRY3D::JVector3D::getLength
double getLength() const
Get length.
Definition: JVector3D.hh:246
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1698
JROOT::JManager
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition: JManager.hh:47
ANTARES::getPhotocathodeArea
double getPhotocathodeArea()
Get photo-cathode area of PMT.
Definition: Antares.hh:51
ANTARES::getAbsorptionLength
double getAbsorptionLength(const double lambda)
Get absorption length.
Definition: Antares.hh:63
GAUSS_LAGUERRE::a
const double a
Definition: JQuadrature.cc:42
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition: JAAnetToolkit.hh:214
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition: JAAnetToolkit.hh:178
JASTRONOMY::getDot
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:676
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
JLANG::getCount
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition: JVectorize.hh:261
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JROOT::counter_type
Long64_t counter_type
Type definition for counter.
Definition: JROOT/JCounter.hh:24
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:867
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:866
KM3NET2D::getPhotocathodeArea2D
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for effective photo-cathode area of PMT.
Definition: KM3NeT2D.hh:5235
std
Definition: JSTDTypes.hh:14
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
JAANET::generator
Description of Monte Carlo event generation applications.
Definition: JHead.hh:469
Generated by doxygen 1.9.1