Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
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
32namespace {
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 */
682int 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 ULong_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 = 37.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
getAbsorptionLength
double getAbsorptionLength(const double lambda)
Definition JDrawPD0.cc:27
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:72
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
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
GAUSS_LAGUERRE::a
const double a
Definition JQuadrature.cc:42
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
KM3NET2D::getPhotocathodeArea2D
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for effective photo-cathode area of PMT.
Definition KM3NeT2D.hh:5235
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.12.0