9#include <framework/logging/Logger.h>
10#include <generators/cosmics/SGCosmic.h>
15#include <framework/dataobjects/Helix.h>
16#include <framework/geometry/BFieldManager.h>
25ofstream ofs_cosmic1g(
"cosmic1g.out");
26ofstream ofs_cosmic1c(
"cosmic1c.out");
27ofstream ofs_cosmic2g(
"cosmic2g.out");
28ofstream ofs_cosmic2c(
"cosmic2c.out");
34 m_parameters.
level = 1;
36 m_parameters.
ipdr = 3.;
37 m_parameters.
ipdz = 3.;
38 m_parameters.
ptmin = 0.7;
50 if (p.level != 1 && p.level != 2) {
51 B2ERROR(
"Wrong generator flag");
68 for (
int i = 0; i < num; i++) {
69 if (max < dim[i]) max = dim[i];
77 double dr, pt, phi, dz, tanl;
88 }
else if (-1 == charge) {
94 p.setFirstDaughter(0);
98 if (1 ==
m_params.
level) ofs_cosmic1g << charge <<
" " << dr <<
" " << phi <<
" "
99 << pt <<
" " << dz <<
" " << tanl << endl;
100 else if (2 ==
m_params.
level) ofs_cosmic2g << charge <<
" " << dr <<
" " << phi <<
" "
101 << pt <<
" " << dz <<
" " << tanl << endl;
106 const float EPS = 0.0001;
107 if (fabs(bz) < EPS) bz = EPS;
108 float d0, phi0, omega, z0, tanLambda;
114 Helix CosmicMCHelix(d0, phi0, omega, z0, tanLambda);
119 arcLength = CosmicMCHelix.getArcLength2DAtCylindricalR(cylindricalR);
120 if (isnan(arcLength))
continue;
123 ROOT::Math::XYZVector vector;
124 vector = CosmicMCHelix.getPositionAtArcLength2D(- arcLength);
125 double vx = vector.
x();
126 double vy = vector.y();
127 double vz = vector.z();
128 vector = CosmicMCHelix.getMomentumAtArcLength2D(- arcLength, bz);
129 double px = vector.
x();
130 double py = vector.y();
131 double pz = vector.z();
132 double m = p.getMass();
133 double e =
sqrt(px * px + py * py + pz * pz + m * m);
136 if (1 ==
m_params.
level) ofs_cosmic1c << vx <<
" " << vy <<
" " << vz <<
" "
137 << px <<
" " << py <<
" " << pz << endl;
138 else if (2 ==
m_params.
level) ofs_cosmic2c << vx <<
" " << vy <<
" " << vz <<
" "
139 << px <<
" " << py <<
" " << pz << endl;
142 p.setMomentum(px, py, pz);
144 p.setProductionVertex(vx, vy, vz);
161 double muPluse = 361917;
162 double muMinus = 285667;
163 double fraction = muPluse / (muPluse + muMinus);
169 double parameters[5];
175 tanl = parameters[4];
180 double muPluse = 143391;
181 double muMinus = 143302;
182 double fraction = muPluse / (muPluse + muMinus);
188 double parameters[5];
194 tanl = parameters[4];
203 rand1 = gRandom->Uniform(0., 1.);
204 if ((
double)rand1 > fraction)
return -1;
214 for (
int i = 0; i < 5; i++) {
218 rand1 = gRandom->Uniform(0., 1.);
219 rand2 = gRandom->Uniform(0., 1.);
222 parameters[0] = rand1 * 40. - 20.;
234 parameters[1] = - rand1 * M_PI;
249 parameters[3] = rand1 * 150. - 50.;
261 parameters[4] = (rand1 - 0.5) * 6;
267 }
while (success == 0);
276 const int width = 40;
277 double dr_pos[nch] = {0.0, 0.0, 0.0, 0.0, 0.0,
278 0.0, 0.0, 0.0, 0.0, 0.0,
279 0.0, 0.0, 0.0, 0.0, 0.0,
280 0.0, 0.0, 0.0, 0.0, 0.0,
281 0.0, 0.0, 0.0, 0.0, 0.0,
282 0.0, 0.0, 0.0, 0.0, 0.0,
283 3284.0, 4110.0, 4545.0, 5140.0, 5601.0,
284 6304.0, 6840.0, 7142.0, 7714.0, 8166.0,
285 8795.0, 9309.0, 9600.0, 9759.0, 9892.0,
286 9995.0, 10198.0, 10379.0, 10428.0, 10646.0,
287 10791.0, 11002.0, 11172.0, 11332.0, 11334.0,
288 11381.0, 11226.0, 11241.0, 10958.0, 10840.0,
289 10619.0, 10246.0, 10135.0, 9640.0, 9477.0,
290 9215.0, 8689.0, 8610.0, 8150.0, 7679.0,
291 359.0, 2.0, 0.0, 0.0, 0.0,
292 0.0, 0.0, 0.0, 0.0, 0.0,
293 0.0, 0.0, 0.0, 0.0, 0.0,
294 0.0, 0.0, 0.0, 0.0, 0.0,
295 0.0, 0.0, 0.0, 0.0, 0.0,
296 0.0, 0.0, 0.0, 0.0, 0.0
299 double max =
findMax(dr_pos, nch);
300 for (
int i = 0; i < nch; i++) dr_pos[i] /= max;
301 int index = (int)((
double)nch / (double)width * (dr + 20.));
303 if (dr_pos[index] > rndm)
return 1;
310 const int width = 40;
311 double dr_neg[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
312 0.0, 0.0, 0.0, 0.0, 0.0,
313 0.0, 0.0, 0.0, 0.0, 0.0,
314 0.0, 0.0, 0.0, 0.0, 0.0,
315 0.0, 0.0, 0.0, 0.0, 0.0,
316 0.0, 0.0, 0.0, 0.0, 0.0,
317 2735.0, 3081.0, 3568.0, 3837.0, 4445.0,
318 4938.0, 5338.0, 5701.0, 5985.0, 6404.0,
319 7022.0, 7375.0, 7438.0, 7662.0, 7833.0,
320 7951.0, 8166.0, 8186.0, 8101.0, 8159.0,
321 8541.0, 8687.0, 8855.0, 8908.0, 8998.0,
322 9014.0, 8962.0, 8915.0, 8833.0, 8698.0,
323 8399.0, 8285.0, 7875.0, 7726.0, 7346.0,
324 7058.0, 7074.0, 6887.0, 6546.0, 5910.0,
325 263.0, 0.0, 0.0, 0.0, 0.0,
326 0.0, 0.0, 0.0, 0.0, 0.0,
327 0.0, 0.0, 0.0, 0.0, 0.0,
328 0.0, 0.0, 0.0, 0.0, 0.0,
329 0.0, 0.0, 0.0, 0.0, 0.0,
330 0.0, 0.0, 0.0, 0.0, 0.0
333 double max =
findMax(dr_neg, nch);
334 for (
int i = 0; i < nch; i++) dr_neg[i] /= max;
335 int index = (int)((
double)nch / (double)width * (dr + 20.));
337 if (dr_neg[index] > rndm)
return 1;
344 const double width = M_PI;
345 double phi_pos[nch] = { 95.0, 194.0, 287.0, 433.0, 636.0,
346 839.0, 1190.0, 1486.0, 1795.0, 2090.0,
347 2192.0, 2437.0, 2831.0, 3114.0, 3407.0,
348 3589.0, 3706.0, 3796.0, 4125.0, 4292.0,
349 4650.0, 4898.0, 4877.0, 5003.0, 5037.0,
350 5409.0, 5488.0, 5704.0, 5733.0, 5547.0,
351 5704.0, 5769.0, 5919.0, 5993.0, 5976.0,
352 5965.0, 6011.0, 6036.0, 6069.0, 6450.0,
353 6260.0, 6103.0, 6252.0, 6277.0, 6219.0,
354 6364.0, 6287.0, 6298.0, 6033.0, 5894.0,
355 5883.0, 6036.0, 6032.0, 5894.0, 5748.0,
356 5229.0, 5382.0, 5374.0, 5504.0, 5283.0,
357 5074.0, 4700.0, 4651.0, 4847.0, 4698.0,
358 4446.0, 4410.0, 4104.0, 3870.0, 3790.0,
359 3826.0, 3750.0, 3488.0, 3220.0, 2973.0,
360 2862.0, 2718.0, 2616.0, 2400.0, 2168.0,
361 1927.0, 1860.0, 1708.0, 1667.0, 1334.0,
362 1162.0, 965.0, 787.0, 629.0, 554.0,
363 384.0, 333.0, 243.0, 159.0, 139.0,
364 88.0, 59.0, 50.0, 27.0, 9.0
367 double max =
findMax(phi_pos, nch);
368 for (
int i = 0; i < nch; i++) phi_pos[i] /= max;
369 int index = (int)((
double)nch / width * (phi + M_PI));
371 if (phi_pos[index] > rndm)
return 1;
378 const double width = M_PI;
379 double phi_neg[nch] = { 351.0, 483.0, 520.0, 506.0, 541.0,
380 497.0, 475.0, 503.0, 495.0, 532.0,
381 594.0, 638.0, 787.0, 904.0, 1094.0,
382 1186.0, 1250.0, 1349.0, 1494.0, 1658.0,
383 1807.0, 1929.0, 2031.0, 2077.0, 2292.0,
384 2402.0, 2592.0, 2790.0, 2988.0, 2908.0,
385 2968.0, 3107.0, 3264.0, 3534.0, 3505.0,
386 3582.0, 3608.0, 3649.0, 3772.0, 3827.0,
387 4023.0, 3960.0, 3956.0, 3965.0, 4095.0,
388 4392.0, 4504.0, 4428.0, 4364.0, 4510.0,
389 4633.0, 4767.0, 4968.0, 4866.0, 4889.0,
390 4753.0, 4862.0, 4928.0, 5015.0, 4958.0,
391 4960.0, 4921.0, 4763.0, 4838.0, 4841.0,
392 4779.0, 4627.0, 4679.0, 4478.0, 4528.0,
393 4470.0, 4474.0, 4199.0, 4145.0, 3896.0,
394 3882.0, 3812.0, 3839.0, 3466.0, 3274.0,
395 3067.0, 3043.0, 2914.0, 2776.0, 2677.0,
396 2425.0, 2195.0, 1892.0, 1777.0, 1685.0,
397 1413.0, 1296.0, 1123.0, 1017.0, 795.0,
398 739.0, 665.0, 493.0, 395.0, 251.0
401 double max =
findMax(phi_neg, nch);
402 for (
int i = 0; i < nch; i++) phi_neg[i] /= max;
403 int index = (int)((
double)nch / width * (phi + M_PI));
405 if (phi_neg[index] > rndm)
return 1;
412 const int width = 30;
413 double pt_pos[nch] = { 47.0, 22547.0, 28678.0, 26159.0, 25145.0,
414 23668.0, 21581.0, 19197.0, 16847.0, 14902.0,
415 13020.0, 11836.0, 10397.0, 9414.0, 8356.0,
416 7492.0, 6640.0, 6173.0, 5426.0, 5023.0,
417 4612.0, 4340.0, 3860.0, 3656.0, 3222.0,
418 3042.0, 2846.0, 2655.0, 2422.0, 2288.0,
419 2172.0, 1984.0, 1841.0, 1745.0, 1599.0,
420 1495.0, 1494.0, 1366.0, 1305.0, 1186.0,
421 1153.0, 1063.0, 1051.0, 996.0, 952.0,
422 906.0, 883.0, 783.0, 733.0, 682.0
425 double max =
findMax(pt_pos, nch);
426 for (
int i = 0; i < nch; i++) pt_pos[i] /= max;
427 int index = (int)((
double)nch / width * pt);
429 if (pt_pos[index] > rndm)
return 1;
436 const int width = 30;
437 double pt_neg[nch] = { 538.0, 614.0, 600.0, 626.0, 648.0,
438 729.0, 716.0, 760.0, 828.0, 882.0,
439 914.0, 922.0, 1009.0, 1105.0, 1113.0,
440 1246.0, 1305.0, 1345.0, 1523.0, 1578.0,
441 1649.0, 1856.0, 1968.0, 2126.0, 2319.0,
442 2526.0, 2716.0, 3005.0, 3317.0, 3427.0,
443 3869.0, 4287.0, 4727.0, 5184.0, 5882.0,
444 6299.0, 7191.0, 8074.0, 9017.0, 9968.0,
445 11425.0, 12821.0, 14601.0, 16351.0, 17885.0,
446 19838.0, 20766.0, 23697.0, 22861.0, 116.0
449 double max =
findMax(pt_neg, nch);
450 for (
int i = 0; i < nch; i++) pt_neg[i] /= max;
451 int index = (int)((
double)nch / width * (width - pt));
453 if (pt_neg[index] > rndm)
return 1;
460 const int width = 250;
461 double dz_pos[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
462 0.0, 0.0, 0.0, 0.0, 0.0,
463 0.0, 0.0, 0.0, 0.0, 0.0,
464 0.0, 0.0, 0.0, 0.0, 0.0,
465 1294.0, 1796.0, 2018.0, 2146.0, 2280.0,
466 2422.0, 2601.0, 2913.0, 3076.0, 3564.0,
467 5147.0, 6284.0, 6958.0, 7133.0, 7465.0,
468 7481.0, 7878.0, 7947.0, 8297.0, 8589.0,
469 8366.0, 8540.0, 8717.0, 8825.0, 8884.0,
470 8985.0, 9131.0, 9105.0, 9279.0, 9333.0,
471 9329.0, 9501.0, 9442.0, 9484.0, 9471.0,
472 9322.0, 9201.0, 8344.0, 6855.0, 5879.0,
473 5095.0, 5012.0, 5044.0, 4967.0, 4717.0,
474 4928.0, 4917.0, 4762.0, 4579.0, 4556.0,
475 4461.0, 4366.0, 4324.0, 4185.0, 4185.0,
476 4085.0, 4019.0, 3974.0, 3748.0, 2739.0,
477 0.0, 0.0, 0.0, 0.0, 0.0,
478 0.0, 0.0, 0.0, 0.0, 0.0,
479 0.0, 0.0, 0.0, 0.0, 0.0,
480 0.0, 0.0, 0.0, 0.0, 0.0
483 double max =
findMax(dz_pos, nch);
484 for (
int i = 0; i < nch; i++) dz_pos[i] /= max;
485 int index = (int)((
double)nch / (double)width * (dz + 100.));
487 if (dz_pos[index] > rndm)
return 1;
494 const int width = 250;
495 double dz_neg[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
496 0.0, 0.0, 0.0, 0.0, 0.0,
497 0.0, 0.0, 0.0, 0.0, 0.0,
498 0.0, 0.0, 0.0, 0.0, 0.0,
499 967.0, 1357.0, 1507.0, 1618.0, 1736.0,
500 1926.0, 2090.0, 2182.0, 2390.0, 2748.0,
501 4261.0, 5004.0, 5605.0, 5819.0, 5932.0,
502 5946.0, 6214.0, 6281.0, 6521.0, 6626.0,
503 6743.0, 6830.0, 6912.0, 7217.0, 7089.0,
504 7194.0, 7340.0, 7283.0, 7343.0, 7357.0,
505 7141.0, 7545.0, 7476.0, 7608.0, 7415.0,
506 7369.0, 7365.0, 6827.0, 5364.0, 4520.0,
507 4052.0, 3994.0, 3885.0, 3820.0, 3743.0,
508 3747.0, 3543.0, 3604.0, 3580.0, 3479.0,
509 3467.0, 3395.0, 3409.0, 3365.0, 3231.0,
510 3247.0, 3192.0, 3074.0, 2995.0, 2215.0,
511 0.0, 0.0, 0.0, 0.0, 0.0,
512 0.0, 0.0, 0.0, 0.0, 0.0,
513 0.0, 0.0, 0.0, 0.0, 0.0,
514 0.0, 0.0, 0.0, 0.0, 0.0
517 double max =
findMax(dz_neg, nch);
518 for (
int i = 0; i < nch; i++) dz_neg[i] /= max;
519 int index = (int)((
double)nch / (double)width * (dz + 100.));
521 if (dz_neg[index] > rndm)
return 1;
529 double tanl_pos[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
530 0.0, 0.0, 0.0, 0.0, 0.0,
531 0.0, 0.0, 0.0, 0.0, 0.0,
532 0.0, 0.0, 0.0, 0.0, 0.0,
533 0.0, 0.0, 0.0, 0.0, 1.0,
534 1.0, 1.0, 2.0, 6.0, 15.0,
535 32.0, 30.0, 103.0, 375.0, 906.0,
536 1465.0, 2088.0, 2858.0, 4007.0, 5324.0,
537 7025.0, 8801.0, 11256.0, 13767.0, 16230.0,
538 18402.0, 20313.0, 21764.0, 22369.0, 23219.0,
539 22975.0, 23058.0, 22121.0, 20633.0, 18402.0,
540 16315.0, 14098.0, 11450.0, 8844.0, 6793.0,
541 5207.0, 3937.0, 2850.0, 1990.0, 1464.0,
542 921.0, 384.0, 83.0, 23.0, 18.0,
543 6.0, 3.0, 2.0, 2.0, 4.0,
544 1.0, 0.0, 1.0, 0.0, 0.0,
545 0.0, 0.0, 0.0, 0.0, 0.0,
546 0.0, 0.0, 0.0, 0.0, 0.0,
547 0.0, 0.0, 0.0, 0.0, 0.0,
548 0.0, 0.0, 0.0, 0.0, 0.0
551 double max =
findMax(tanl_pos, nch);
552 for (
int i = 0; i < nch; i++) tanl_pos[i] /= max;
553 int index = (int)((
double)nch / (double)width * (tanl + 3.));
555 if (tanl_pos[index] > rndm)
return 1;
563 double tanl_neg[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
564 0.0, 0.0, 0.0, 0.0, 0.0,
565 0.0, 0.0, 0.0, 0.0, 0.0,
566 0.0, 0.0, 0.0, 0.0, 0.0,
567 0.0, 0.0, 0.0, 0.0, 0.0,
568 1.0, 1.0, 0.0, 7.0, 10.0,
569 11.0, 24.0, 94.0, 335.0, 760.0,
570 1210.0, 1713.0, 2233.0, 3242.0, 4236.0,
571 5530.0, 7004.0, 8768.0, 10923.0, 12809.0,
572 14457.0, 16043.0, 17321.0, 18064.0, 18237.0,
573 18129.0, 17823.0, 17077.0, 16272.0, 14539.0,
574 13015.0, 11031.0, 9000.0, 7219.0, 5411.0,
575 4163.0, 3047.0, 2120.0, 1551.0, 1215.0,
576 672.0, 266.0, 73.0, 23.0, 5.0,
577 4.0, 8.0, 5.0, 3.0, 0.0,
578 0.0, 0.0, 0.0, 0.0, 0.0,
579 0.0, 0.0, 0.0, 0.0, 0.0,
580 0.0, 0.0, 0.0, 0.0, 0.0,
581 0.0, 0.0, 0.0, 0.0, 0.0,
582 0.0, 0.0, 0.0, 0.0, 0.0
585 double max =
findMax(tanl_neg, nch);
586 for (
int i = 0; i < nch; i++) tanl_neg[i] /= max;
587 int index = (int)((
double)nch / (double)width * (tanl + 3.));
589 if (tanl_neg[index] > rndm)
return 1;
599 for (
int i = 0; i < 5; i++) {
603 rand1 = gRandom->Uniform(0., 1.);
604 rand2 = gRandom->Uniform(0., 1.);
606 parameters[0] = rand1 * 6. - 3.;
613 parameters[1] = - rand1 * M_PI;
620 parameters[2] = 50. * rand1;
627 parameters[3] = rand1 * 80. - 30.;
634 parameters[4] = rand1 * 4 - 2;
640 }
while (success == 0);
649 double dr_pos[nch] = {0., 0., 0., 0., 0.,
650 0., 2403.000, 3562.000, 3685.000, 3619.000,
651 3667.000, 3724.000, 3823.000, 3824.000, 3881.000,
652 3764.000, 3867.000, 3918.000, 3964.000, 3916.000,
653 4002.000, 4001.000, 3980.000, 3848.000, 4220.000,
654 4206.000, 3856.000, 3920.000, 3857.000, 3751.000,
655 3753.000, 3898.000, 3898.000, 3822.000, 3872.000,
656 3905.000, 3896.000, 3874.000, 3767.000, 3734.000,
657 3744.000, 3666.000, 3757.000, 2458.000, 0.,
661 double max =
findMax(dr_pos, nch);
662 for (
int i = 0; i < nch; i++) dr_pos[i] /= max;
663 int index = (int)((
double)nch / (double)width * (dr + 3.));
665 if (dr_pos[index] > rndm)
return 1;
673 double dr_neg[nch] = { 0., 0., 0., 0., 0.,
674 0., 2508.000, 3752.000, 3713.000, 3695.000,
675 3758.000, 3727.000, 3847.000, 3887.000, 3896.000,
676 3862.000, 3858.000, 3905.000, 3894.000, 3806.000,
677 3736.000, 3797.000, 3921.000, 3922.000, 4203.000,
678 4167.000, 3927.000, 3952.000, 3983.000, 3982.000,
679 3914.000, 3962.000, 3983.000, 3817.000, 3831.000,
680 3811.000, 3837.000, 3852.000, 3728.000, 3684.000,
681 3629.000, 3663.000, 3588.000, 2394.000, 0.,
685 double max =
findMax(dr_neg, nch);
686 for (
int i = 0; i < nch; i++) dr_neg[i] /= max;
687 int index = (int)((
double)nch / (double)width * (dr + 3.));
689 if (dr_neg[index] > rndm)
return 1;
696 const double width = 8;
697 double phi_pos[nch] = { 0., 0., 0., 0., 0.,
698 1073.000, 1911.000, 2905.000, 4165.000, 5205.000,
699 6062.000, 6955.000, 7382.000, 7262.000, 7142.000,
700 6548.000, 5759.000, 5015.000, 4399.000, 3392.000,
701 2365.000, 1533.000, 801.0000, 245.0000, 25.00000,
702 39.00000, 119.0000, 398.0000, 841.0000, 1487.000,
703 2331.000, 3138.000, 3726.000, 4430.000, 5106.000,
704 5548.000, 5766.000, 5699.000, 5408.000, 5011.000,
705 4439.000, 3732.000, 2883.000, 2034.000, 1023.000,
709 double max =
findMax(phi_pos, nch);
710 for (
int i = 0; i < nch; i++) phi_pos[i] /= max;
711 int index = (int)(nch / width * (phi + 4.));
713 if (phi_pos[index] > rndm)
return 1;
720 const double width = 8;
721 double phi_neg[nch] = { 0., 0., 0., 0., 0.,
722 16.00000, 73.00000, 262.0000, 677.0000, 1252.000,
723 1973.000, 2896.000, 3520.000, 4173.000, 4876.000,
724 5389.000, 5730.000, 5798.000, 5553.000, 5120.000,
725 4676.000, 4036.000, 3185.000, 2320.000, 1688.000,
726 1741.000, 2164.000, 3325.000, 4588.000, 5522.000,
727 6427.000, 7159.000, 7385.000, 7282.000, 6980.000,
728 6157.000, 5496.000, 4854.000, 4066.000, 3042.000,
729 2032.000, 1236.000, 586.0000, 123.0000, 13.00000,
733 double max =
findMax(phi_neg, nch);
734 for (
int i = 0; i < nch; i++) phi_neg[i] /= max;
735 int index = (int)(nch / width * (phi + 4.));
737 if (phi_neg[index] > rndm)
return 1;
744 const int width = 50;
745 double pt_pos[nch] = { 3813.000, 14309.00, 14985.00, 13103.00, 10807.00,
746 9118.000, 7793.000, 6718.000, 5788.000, 5177.000,
747 4458.000, 3889.000, 3480.000, 3191.000, 2798.000,
748 2484.000, 2245.000, 2053.000, 1883.000, 1691.000,
749 1481.000, 1416.000, 1311.000, 1207.000, 1171.000,
750 1017.000, 953.0000, 886.0000, 790.0000, 816.0000,
751 712.0000, 685.0000, 627.0000, 577.0000, 526.0000,
752 502.0000, 442.0000, 453.0000, 415.0000, 351.0000,
753 340.0000, 358.0000, 305.0000, 330.0000, 278.0000,
754 272.0000, 250.0000, 273.0000, 243.0000, 209.0000,
755 203.0000, 176.0000, 212.0000, 182.0000, 153.0000,
756 168.0000, 155.0000, 170.0000, 129.0000, 117.0000,
757 120.0000, 114.0000, 126.0000, 123.0000, 106.0000,
758 90.00000, 97.00000, 103.0000, 82.00000, 87.00000,
759 76.00000, 104.0000, 75.00000, 82.00000, 67.00000,
760 64.00000, 63.00000, 75.00000, 65.00000, 62.00000,
761 54.00000, 86.00000, 44.00000, 44.00000, 53.00000,
762 54.00000, 53.00000, 51.00000, 49.00000, 39.00000,
763 37.00000, 42.00000, 40.00000, 40.00000, 24.00000,
764 38.00000, 39.00000, 30.00000, 28.00000, 32.00000
767 double max =
findMax(pt_pos, nch);
768 for (
int i = 0; i < nch; i++) pt_pos[i] /= max;
769 int index = (int)(nch / width * pt);
771 if (pt_pos[index] > rndm)
return 1;
778 const int width = 50;
779 double pt_neg[nch] = { 3794.000, 14288.00, 15015.00, 13120.00, 10865.00,
780 9167.000, 7826.000, 6733.000, 5760.000, 5253.000,
781 4506.000, 3796.000, 3512.000, 3197.000, 2848.000,
782 2463.000, 2232.000, 2083.000, 1827.000, 1679.000,
783 1556.000, 1414.000, 1299.000, 1163.000, 1116.000,
784 1069.000, 964.0000, 922.0000, 818.0000, 764.0000,
785 723.0000, 658.0000, 614.0000, 541.0000, 552.0000,
786 520.0000, 491.0000, 382.0000, 380.0000, 392.0000,
787 346.0000, 352.0000, 321.0000, 298.0000, 272.0000,
788 277.0000, 246.0000, 238.0000, 237.0000, 241.0000,
789 189.0000, 208.0000, 196.0000, 194.0000, 169.0000,
790 149.0000, 143.0000, 139.0000, 130.0000, 134.0000,
791 125.0000, 116.0000, 110.0000, 104.0000, 122.0000,
792 110.0000, 99.00000, 95.00000, 91.00000, 87.00000,
793 91.00000, 75.00000, 87.00000, 76.00000, 81.00000,
794 81.00000, 67.00000, 64.00000, 55.00000, 59.00000,
795 52.00000, 53.00000, 53.00000, 51.00000, 46.00000,
796 40.00000, 49.00000, 41.00000, 37.00000, 44.00000,
797 44.00000, 42.00000, 47.00000, 30.00000, 30.00000,
798 42.00000, 28.00000, 28.00000, 32.00000, 26.00000
801 double max =
findMax(pt_neg, nch);
802 for (
int i = 0; i < nch; i++) pt_neg[i] /= max;
803 int index = (int)(nch / width * pt);
805 if (pt_neg[index] > rndm)
return 1;
812 const int width = 100;
813 double dz_pos[nch] = { 0., 0., 0., 1.000000, 0.,
820 0., 0., 0., 1.000000, 0.,
821 1.000000, 0., 0., 4.000000, 1.000000,
822 1.000000, 2.000000, 1.000000, 3.000000, 6.000000,
823 20.00000, 11.00000, 30.00000, 45.00000, 64.00000,
824 139.0000, 243.0000, 391.0000, 460.0000, 586.0000,
825 720.0000, 790.0000, 931.0000, 967.0000, 1087.000,
826 1204.000, 1231.000, 1272.000, 1251.000, 1355.000,
827 1453.000, 1447.000, 1370.000, 1364.000, 1460.000,
828 1405.000, 1545.000, 1454.000, 1536.000, 1499.000,
829 1547.000, 1500.000, 1465.000, 1512.000, 1515.000,
830 1547.000, 1548.000, 1554.000, 1579.000, 1599.000,
831 1590.000, 1609.000, 1618.000, 1710.000, 1599.000,
832 1746.000, 1669.000, 1709.000, 1814.000, 1904.000,
833 1896.000, 1634.000, 1598.000, 1680.000, 1727.000,
834 1714.000, 1702.000, 1737.000, 1693.000, 1743.000,
835 1713.000, 1663.000, 1639.000, 1569.000, 1695.000,
836 1697.000, 1745.000, 1708.000, 1682.000, 1686.000,
837 1676.000, 1638.000, 1660.000, 1655.000, 1592.000,
838 1617.000, 1581.000, 1514.000, 1601.000, 1553.000,
839 1451.000, 1528.000, 1456.000, 1489.000, 1465.000,
840 1455.000, 1489.000, 1375.000, 1387.000, 1330.000,
841 1298.000, 1304.000, 1237.000, 1140.000, 1111.000,
842 1038.000, 942.0000, 959.0000, 867.0000, 822.0000,
843 773.0000, 698.0000, 649.0000, 591.0000, 556.0000,
844 513.0000, 474.0000, 424.0000, 408.0000, 414.0000,
845 324.0000, 296.0000, 265.0000, 259.0000, 188.0000,
846 170.0000, 195.0000, 163.0000, 140.0000, 131.0000,
847 122.0000, 113.0000, 98.00000, 92.00000, 68.00000,
848 61.00000, 51.00000, 44.00000, 46.00000, 30.00000,
849 32.00000, 25.00000, 20.00000, 13.00000, 11.00000,
850 3.000000, 2.000000, 1.000000, 1.000000, 0.,
855 double max =
findMax(dz_pos, nch);
856 for (
int i = 0; i < nch; i++) dz_pos[i] /= max;
857 int index = (int)((
double)nch / (double)width * (dz + 50.));
859 if (dz_pos[index] > rndm)
return 1;
867 double dz_neg[nch] = { 0., 1.000000, 0., 0., 0.,
870 0., 0., 0., 0., 1.000000,
871 2.000000, 1.000000, 2.000000, 5.000000, 12.00000,
872 25.00000, 66.00000, 244.0000, 609.0000, 1092.000,
873 1592.000, 1849.000, 2264.000, 2537.000, 2692.000,
874 2889.000, 2711.000, 2901.000, 3015.000, 2976.000,
875 3061.000, 3001.000, 3066.000, 3058.000, 3146.000,
876 3198.000, 3308.000, 3349.000, 3391.000, 3720.000,
877 3552.000, 3291.000, 3404.000, 3459.000, 3447.000,
878 3318.000, 3248.000, 3416.000, 3408.000, 3408.000,
879 3334.000, 3345.000, 3251.000, 3127.000, 2999.000,
880 3038.000, 2944.000, 2970.000, 2921.000, 2603.000,
881 2633.000, 2304.000, 2209.000, 1866.000, 1660.000,
882 1511.000, 1269.000, 1065.000, 893.0000, 765.0000,
883 628.0000, 456.0000, 416.0000, 378.0000, 273.0000,
884 221.0000, 166.0000, 134.0000, 116.0000, 62.00000,
885 46.00000, 35.00000, 11.00000, 5.000000, 0.,
886 0., 0., 0., 1.000000, 0.
889 double max =
findMax(dz_neg, nch);
890 for (
int i = 0; i < nch; i++) dz_neg[i] /= max;
892 int index = (int)(dz + 50.);
893 if (dz_neg[index] > rndm)
return 1;
901 double tanl_pos[nch] = { 0., 0., 0., 0., 0.,
902 1.000000, 0., 2.000000, 0., 7.000000,
903 2.000000, 6.000000, 4.000000, 17.00000, 20.00000,
904 21.00000, 18.00000, 34.00000, 36.00000, 47.00000,
905 64.00000, 85.00000, 126.0000, 116.0000, 163.0000,
906 214.0000, 262.0000, 330.0000, 415.0000, 477.0000,
907 556.0000, 641.0000, 852.0000, 956.0000, 1243.000,
908 1370.000, 1699.000, 1995.000, 2415.000, 3040.000,
909 3352.000, 3889.000, 4390.000, 4823.000, 5399.000,
910 5862.000, 6201.000, 6506.000, 6748.000, 6788.000,
911 6767.000, 6575.000, 6479.000, 6112.000, 5940.000,
912 5459.000, 4986.000, 4504.000, 3879.000, 3460.000,
913 2981.000, 2502.000, 2121.000, 1803.000, 1477.000,
914 1248.000, 1041.000, 927.0000, 708.0000, 605.0000,
915 461.0000, 424.0000, 341.0000, 257.0000, 230.0000,
916 167.0000, 141.0000, 119.0000, 84.00000, 80.00000,
917 46.00000, 44.00000, 39.00000, 30.00000, 17.00000,
918 18.00000, 13.00000, 6.000000, 6.000000, 5.000000,
919 2.000000, 2.000000, 2.000000, 0., 0.,
920 0., 0., 0., 0., 1.000000
923 double max =
findMax(tanl_pos, nch);
924 for (
int i = 0; i < nch; i++) tanl_pos[i] /= max;
925 int index = (int)((
double)nch / (double)width * (tanl + 2.5));
927 if (tanl_pos[index] > rndm)
return 1;
935 double tanl_neg[nch] = { 0., 0., 0., 0., 0.,
936 0., 0., 1.000000, 4.000000, 0.,
937 3.000000, 4.000000, 9.000000, 8.000000, 14.00000,
938 21.00000, 25.00000, 41.00000, 44.00000, 49.00000,
939 73.00000, 84.00000, 107.0000, 139.0000, 166.0000,
940 242.0000, 268.0000, 316.0000, 412.0000, 469.0000,
941 595.0000, 727.0000, 900.0000, 1002.000, 1267.000,
942 1454.000, 1793.000, 2126.000, 2500.000, 2966.000,
943 3492.000, 3829.000, 4505.000, 4980.000, 5518.000,
944 5864.000, 6157.000, 6440.000, 6664.000, 6747.000,
945 6841.000, 6757.000, 6454.000, 6224.000, 5776.000,
946 5458.000, 4855.000, 4399.000, 3887.000, 3323.000,
947 3098.000, 2445.000, 1993.000, 1695.000, 1381.000,
948 1200.000, 993.0000, 844.0000, 640.0000, 590.0000,
949 484.0000, 412.0000, 352.0000, 255.0000, 207.0000,
950 160.0000, 121.0000, 119.0000, 103.0000, 57.00000,
951 53.00000, 37.00000, 27.00000, 26.00000, 25.00000,
952 23.00000, 8.000000, 7.000000, 9.000000, 3.000000,
953 5.000000, 4.000000, 2.000000, 3.000000, 3.000000,
957 double max =
findMax(tanl_neg, nch);
958 for (
int i = 0; i < nch; i++) tanl_neg[i] /= max;
959 int index = (int)((
double)nch / (double)width * (tanl + 2.5));
961 if (tanl_neg[index] > rndm)
return 1;
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Class to represent Particle data in graph.
Class to build, validate and sort a particle decay chain.
@ c_PrimaryParticle
bit 0: Particle is primary particle.
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
SGCosmic()
Default constructor.
int mkPhi_pos_v2(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for positively charged particles.
int mkTanl_neg_v1(const double tanl, const float rndm)
Generates tangent of the polar angle tanl distributions by accept-reject method for negatively charge...
int mkDz_pos_v1(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for positively charged particles.
int mkTanl_pos_v1(const double tanl, const float rndm)
Generates tangent of the polar angle tanl distributions by accept-reject method for positively charge...
int mkPt_neg_v1(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for negatively charged particl...
void genCosmic(const int level, int &charge, double &dr, double &phi, double &Pt, double &dz, double &tanl)
Generates cosmic events according to tabulated distributions in 5-dimensional space: dr,...
int mkPhi_neg_v1(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for negatively charged particles.
int mkDr_neg_v1(const double dr, const float rndm)
Generates vertex distributions in the radial direction dr by accept-reject method for negatively.
void mkdist_v2(const int charge, double *)
Generates distributions in 5-parameter space for different particle charges.
int mkPt_pos_v2(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for positively charged particl...
Parameters m_params
All relevant parameters.
bool generateEvent(MCParticleGraph &graph)
Generates the next event and store the result in the given MCParticle graph.
int mkDz_neg_v1(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for negatively charged particles.
int mkPt_pos_v1(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for positively charged particl...
int mkPhi_neg_v2(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for negatively charged particles.
int mkDr_neg_v2(const double dr, const float rndm)
Generates vertex distributions in the radial direction dr by accept-reject method for negatively char...
int mkDr_pos_v2(const double dr, const float rndm)
Generates vertex distributions in the radial direction dr by accept-reject method for positively char...
int mkTanl_pos_v2(const double tanl, const float rndm)
Generates tangent of the polar angle tanl distributions by accept-reject method for positively charge...
int mkTanl_neg_v2(const double tanl, const float rndm)
Generates tangent of the polar angle tanl distributions by accept-reject method for negatively charge...
int mkDz_neg_v2(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for negatively charged particles.
bool setParameters(const Parameters ¶meters)
Sets the parameters for generating the Particles.
void mkdist_v1(const int charge, double *)
Generates distributions in 5-parameter space for different particle charges.
int muChargeFlag(const double)
Generates the muon charge according to the positively/negatively charged muon ratio.
int mkPt_neg_v2(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for negatively charged particl...
int mkDz_pos_v2(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for positively charged particles.
int mkDr_pos_v1(const double dr, const float rndm)
Generates vertex distributions in the radial direction dr by accept-reject method for positively char...
int mkPhi_pos_v1(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for positively charged particles.
double findMax(const double *dim, const int num)
Finds maximum value in an array.
static const double T
[tesla]
static void getField(const double *pos, double *field)
return the magnetic field at a given position.
double sqrt(double a)
sqrt for double
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.
Struct to keep all necessary parameters for the cosmic generator.
int ipRequirement
Restrict the vertex to IP or not (default)
double ipdr
Vertex restriction in the radial direction.
double cylindricalR
Cylindrical radius of generation.
double ptmin
Minimum value of the transverse momentum.
int level
Generator version: level 1 (default) or 2.
double ipdz
Vertex restriction in z direction.