51 virtual const JResult& next() = 0;
59 template<
int N>
struct JEnigma;
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())
79 this->V = (4.0/3.0)*PI*(D3_max - D3_min);
82 virtual const JResult& next()
84 const double D3 = gRandom->Uniform(D3_min, D3_max);
109 JEnigma(
const JRange_t& D_m) :
110 D_min(D_m.getLowerLimit()),
111 D_max(D_m.getUpperLimit())
114 virtual const JResult& next()
116 this->D = gRandom->Uniform(D_min, D_max);
117 this->V = 4.0*PI*D*D * (D_max - D_min);
140 JEnigma(
const JRange_t& D_m) :
141 Dinv_min(1.0/D_m.getUpperLimit()),
142 Dinv_max(1.0/D_m.getLowerLimit())
145 virtual const JResult& next()
147 const double Dinv = gRandom->Uniform(Dinv_min, Dinv_max);
150 this->V = 4.0*PI*D*D * D*D * (Dinv_max - Dinv_min);
155 const double Dinv_min;
156 const double Dinv_max;
172 double operator()(
const double x)
const
174 const_iterator i = this->lower_bound(x);
176 if (i != this->begin()) {
189 class JK40BetaDecay :
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;
409 static double getBranchingRatio()
420 class JK40ElectronCapture :
427 JK40ElectronCapture()
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;
640 static double getBranchingRatio()
649 static const JK40BetaDecay k40_beta_decay;
650 static const JK40ElectronCapture k40_electron_capture;
664 inline double get_angular_acceptance(
const double ct)
666 const double w = 0.25*
a*(1.0 + ct)*(1.0 + ct) - ct;
682int main(
int argc,
char* argv[])
703 JParser<> zap(
"Example program to calculate multiples rate.");
708 zap[
'D'] =
make_field(D_m) = JRange_t(0.216, 10);
720 catch(
const exception &error) {
721 FATAL(error.what() << endl);
724 gRandom->SetSeed(seed);
726 using namespace NAMESPACE;
730 const JModule module = getModule<JKM3NeT_t>(
id);
732 DEBUG(module << endl);
736 const double R_m = 17.0 * 2.54 * 0.5e-2;
737 const double A = PI * R_m * R_m;
739 const double wmin = 280.0;
740 const double wmax = 700.0;
743 const double WAVELENGTH_EXPANSION = (wmax-wmin) / (wmin*wmax) * (300.0*600.0)/(600.0-300.0);
745 JGenerator* enigma = NULL;
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;
756 const double vmin = 1.0 / wmax;
757 const double vmax = 1.0 / wmin;
761 for (
double w = wmin; w <= wmax; w += 1.0) {
762 if (getQE(w) > QEmax) {
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);
771 typedef JManager<int, TH1D> JManager_t;
773 JManager_t H1(
new TH1D(
"M[%]", NULL, 100, D_m.getLowerLimit(), D_m.getUpperLimit()));
777 TH1D pmt(
"pmt", NULL, 1000, -1.0, +1.0);
779 for (Int_t i = 1; i != pmt.GetNbinsX(); ++i) {
781 const double dot = pmt.GetBinCenter(i);
792 y = get_angular_acceptance(dot);
796 pmt.SetBinContent(i, y);
807 for (
counter_type event_count = 0; event_count != numberOfEvents; ++event_count) {
809 if (event_count%10000 == 0) {
810 STATUS(
"event: " << setw(10) << event_count <<
'\r');
DEBUG(endl);
813 const JResult& result = enigma->next();
815 const double D = result.D;
816 const double V = result.V;
821 double W = A / (4*PI*(D-R_m)*(D-R_m));
828 double x = gRandom->Rndm();
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());
836 const int N = gRandom->Poisson(y * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
845 const double ct = gRandom->Uniform(-1.0, +1.0);
846 const double phi = gRandom->Uniform(-PI, +PI);
848 const double st = sqrt((1.0 - ct) * (1.0 + ct));
856 for (
int i = 0; i != N; ++i) {
860 const double v = gRandom->Uniform(vmin, vmax);
861 const double w = 1.0 / v;
867 for (
size_t pmt = 0; pmt !=
module.size(); ++pmt) {
878 ERROR(
"Distance " << d <<
" < " << D << endl);
881 const double dot = getDot(pos, module[pmt].getDirection());
883 const double U = 0.5 * (1.0 - d / sqrt(d*d + getPhotocathodeArea() / PI));
898 p = get_angular_acceptance(dot) * getQE(w);
902 P += pi[pmt] = U * p * exp(-d/l_abs);
906 ERROR(
"Probability " << P <<
" > " << W << endl);
909 if (W * QEmax * gRandom->Rndm() < P) {
912 double y = gRandom->Uniform(P);
914 for (vector<double>::const_iterator i = pi.begin(); i != pi.end() && (y -= *i) > 0.0; ++i, ++pmt) {}
916 buffer.push_back(pmt);
920 if (!buffer.empty()) {
922 int M = buffer.size();
926 sort(buffer.begin(), buffer.end());
928 M =
distance(buffer.begin(), unique(buffer.begin(), buffer.end()));
934 for (
int i = 2; i <= M; ++i) {
935 P2[i].put((
double) (buffer.size() - M) / (
double) M, V);
943 for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
944 i->second->Scale(bequerel / (
double) numberOfEvents);
947 for (
size_t M = 2; M != 7; ++M) {
948 cout <<
"Rate[" << M <<
"] = "
949 <<
FIXED(8,3) << bequerel * h1[M].getTotal() / (double) numberOfEvents
951 <<
FIXED(7,3) << bequerel * h1[M].getError() / (double) numberOfEvents
955 for (
size_t M = 2; M != 7; ++M) {
956 if (P2[M].getCount() != 0) {
957 cout <<
"P2[" << M <<
"] = " << P2[M].getMean() << endl;
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
double getAbsorptionLength(const double lambda)
Dynamic ROOT object management.
General purpose messaging.
#define DEBUG(A)
Message macros.
Data structure for optical module.
int main(int argc, char *argv[])
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
I/O formatting auxiliaries.
Auxiliary class to define a range between two values.
Properties of KM3NeT PMT and deep-sea water.
Properties of KM3NeT PMT and deep-sea water.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Data structure for a composite optical module.
Data structure for position in three dimensions.
double getLength() const
Get length.
Utility class to parse command line options.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Long64_t counter_type
Type definition for counter.
static const JPhotocathodeArea2D getPhotocathodeArea2D
Function object for effective photo-cathode area of PMT.
Auxiliary data structure for floating point format specification.
Description of Monte Carlo event generation applications.