Belle II Software development
SGCosmic.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <framework/logging/Logger.h>
10#include <generators/cosmics/SGCosmic.h>
11
12#include <cmath>
13#include <fstream>
14
15#include <framework/dataobjects/Helix.h>
16#include <framework/geometry/BFieldManager.h>
17#include <TRandom.h>
18#include <TVector3.h>
19
20using namespace std;
21using namespace Belle2;
22
23//#define DEBUG
24#ifdef DEBUG
25ofstream ofs_cosmic1g("cosmic1g.out");
26ofstream ofs_cosmic1c("cosmic1c.out");
27ofstream ofs_cosmic2g("cosmic2g.out");
28ofstream ofs_cosmic2c("cosmic2c.out");
29#endif
30
32{
33 Parameters m_parameters;
34 m_parameters.level = 1;
35 m_parameters.ipRequirement = 0;
36 m_parameters.ipdr = 3.; // Only relevant for ipRequirement = 1
37 m_parameters.ipdz = 3.; // Only relevant for ipRequirement = 1
38 m_parameters.ptmin = 0.7;
39 m_parameters.cylindricalR = 125.0;
40 setParameters(m_parameters);
41}
42
43
45{
46 // Sanity checks
47 bool ok(true);
48
49 // Check that we have correct generator level
50 if (p.level != 1 && p.level != 2) {
51 B2ERROR("Wrong generator flag");
52 ok = false;
53 }
54
55 // If everything is ok, set the new parameters, else return false
56 if (ok) {
57 m_params = p;
58 return true;
59 }
60 return false;
61}
62
63
64// Utility from the Belle cosmics version
65double SGCosmic::findMax(const double* dim, const int num)
66{
67 double max = 0.;
68 for (int i = 0; i < num; i++) {
69 if (max < dim[i]) max = dim[i];
70 }
71 return max;
72}
73
75{
76 int charge = 0; // Initialization of charge
77 double dr, pt, phi, dz, tanl;
78
79 while (true) {
80
81 genCosmic(m_params.level, charge, dr, phi, pt, dz, tanl);
82
85
86 if (1 == charge) {
87 p.setPDG(-13);
88 } else if (-1 == charge) {
89 p.setPDG(13);
90 } else {
91 continue;
92 }
93 p.setMassFromPDG();
94 p.setFirstDaughter(0);
95 p.setLastDaughter(0);
96
97#ifdef DEBUG
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;
102#endif
103
104 // Simulate helix parameter at perigee
105 float bz = BFieldManager::getField(0, 0, 0).Z() / Unit::T; // Magnetic field
106 const float EPS = 0.0001; // Avoid using zero magnetic field
107 if (fabs(bz) < EPS) bz = EPS; // Set the value of the magnetic field to a small number
108 float d0, phi0, omega, z0, tanLambda; // The definition is the same as in the Helix class
109 d0 = dr;
110 phi0 = phi;
111 omega = (double)charge / (pt * Helix::getAlpha(bz));
112 z0 = dz;
113 tanLambda = tanl;
114 Helix CosmicMCHelix(d0, phi0, omega, z0, tanLambda);
115
116 const float cylindricalR = m_params.cylindricalR; // radius (cm) of generation
117 float arcLength;
118 // Get arc length at ToP radius
119 arcLength = CosmicMCHelix.getArcLength2DAtCylindricalR(cylindricalR);
120 if (isnan(arcLength)) continue;
121
122 // Calculate coordinates and momentum at ToP radius
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);
134
135#ifdef DEBUG
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;
140#endif
141
142 p.setMomentum(px, py, pz);
143 p.setEnergy(e);
144 p.setProductionVertex(vx, vy, vz);
146
147 return true;
148 }
149 return false;
150}
151
152void SGCosmic::genCosmic(const int level, int& charge,
153 double& dr,
154 double& phi,
155 double& pt,
156 double& dz,
157 double& tanl)
158{
159 // generator version 1
160 if (level == 1) {
161 double muPluse = 361917;
162 double muMinus = 285667;
163 double fraction = muPluse / (muPluse + muMinus);
164
165 // mu+ or mu- --> fraction = f(mu+)
166 charge = muChargeFlag(fraction);
167
168 // make distribution
169 double parameters[5]; // dr,phi,pt(not Kappa),dz,tanl
170 mkdist_v1(charge, parameters);
171 dr = parameters[0];
172 phi = parameters[1];
173 pt = parameters[2]; // not Kappa
174 dz = parameters[3];
175 tanl = parameters[4];
176 }
177
178 // generator version 2
179 if (level == 2) {
180 double muPluse = 143391;
181 double muMinus = 143302;
182 double fraction = muPluse / (muPluse + muMinus);
183
184 // mu+ or mu- --> fraction = f(mu+)
185 charge = muChargeFlag(fraction);
186
187 // make distribution
188 double parameters[5]; // dr,phi,pt(not Kappa),dz,tanl
189 mkdist_v2(charge, parameters);
190 dr = parameters[0];
191 phi = parameters[1];
192 pt = parameters[2]; // not Kappa
193 dz = parameters[3];
194 tanl = parameters[4];
195 }
196
197 return;
198}
199
200int SGCosmic::muChargeFlag(const double fraction)
201{
202 double rand1;
203 rand1 = gRandom->Uniform(0., 1.);
204 if ((double)rand1 > fraction) return -1;
205 else return 1;
206}
207
208// -------- version 1 --------
209void SGCosmic::mkdist_v1(const int charge, double* parameters)
210{
211 double rand1;
212 double rand2;
213
214 for (int i = 0; i < 5; i++) {
215 int success;
216 success = 0;
217 do {
218 rand1 = gRandom->Uniform(0., 1.);
219 rand2 = gRandom->Uniform(0., 1.);
220 if (i == 0) { // dr
221 if (m_params.ipRequirement == 0) {
222 parameters[0] = rand1 * 40. - 20.; // (-20, 20)
223 } else {
224 // IP requirement
225 double widthdr = m_params.ipdr * 2;
226 parameters[0] = rand1 * widthdr - m_params.ipdr; // (-ipdr, ipdr)
227 }
228 if (charge == 1)
229 success = mkDr_pos_v1(parameters[0], rand2);
230 else
231 success = mkDr_neg_v1(parameters[0], rand2);
232 }
233 if (i == 1) { // phi
234 parameters[1] = - rand1 * M_PI; // (-pi, 0)
235 if (charge == 1)
236 success = mkPhi_pos_v1(parameters[1], rand2);
237 else
238 success = mkPhi_neg_v1(parameters[1], rand2);
239 }
240 if (i == 2) { // pt
241 parameters[2] = (30.0 - m_params.ptmin) * rand1 + m_params.ptmin; // (ptmin, 30)
242 if (charge == 1)
243 success = mkPt_pos_v1(parameters[2], rand2);
244 else
245 success = mkPt_neg_v1(parameters[2], rand2);
246 }
247 if (i == 3) { // dz
248 if (m_params.ipRequirement == 0) {
249 parameters[3] = rand1 * 150. - 50.; //(-50, 100)
250 } else {
251 // IP requirement
252 double widthdz = m_params.ipdz * 2;
253 parameters[3] = rand1 * widthdz - m_params.ipdz; //(-ipdz, ipdz)
254 }
255 if (charge == 1)
256 success = mkDz_pos_v1(parameters[3], rand2);
257 else
258 success = mkDz_neg_v1(parameters[3], rand2);
259 }
260 if (i == 4) { // tanl
261 parameters[4] = (rand1 - 0.5) * 6; // (-3, 3)
262 if (charge == 1)
263 success = mkTanl_pos_v1(parameters[4], rand2);
264 else
265 success = mkTanl_neg_v1(parameters[4], rand2);
266 }
267 } while (success == 0);
268 }
269 return;
270}
271
272
273int SGCosmic::mkDr_pos_v1(const double dr, const float rndm)
274{
275 const int nch = 100;
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
297 };
298 // normalize by peak value --> max = 1
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.));
302
303 if (dr_pos[index] > rndm) return 1;
304 return 0;
305}
306
307int SGCosmic::mkDr_neg_v1(const double dr, const float rndm)
308{
309 const int nch = 100;
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
331 };
332 // normalize by peak value --> max = 1
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.));
336
337 if (dr_neg[index] > rndm) return 1;
338 return 0;
339}
340
341int SGCosmic::mkPhi_pos_v1(const double phi, const float rndm)
342{
343 const int nch = 100;
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
365 };
366 // normalize by peak value --> max = 1
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));
370
371 if (phi_pos[index] > rndm) return 1;
372 return 0;
373}
374
375int SGCosmic::mkPhi_neg_v1(const double phi, const float rndm)
376{
377 const int nch = 100;
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
399 };
400 // normalize by peak value --> max = 1
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));
404
405 if (phi_neg[index] > rndm) return 1;
406 return 0;
407}
408
409int SGCosmic::mkPt_pos_v1(const double pt, const float rndm)
410{
411 const int nch = 50;
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
423 };
424 // normalize by peak value --> max = 1
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);
428
429 if (pt_pos[index] > rndm) return 1;
430 return 0;
431}
432
433int SGCosmic::mkPt_neg_v1(const double pt, const float rndm)
434{
435 const int nch = 50;
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
447 };
448 // normalize by peak value --> max = 1
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));
452
453 if (pt_neg[index] > rndm) return 1;
454 return 0;
455}
456
457int SGCosmic::mkDz_pos_v1(const double dz, const float rndm)
458{
459 const int nch = 100;
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
481 };
482 // normalize by peak value --> max = 1
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.));
486
487 if (dz_pos[index] > rndm) return 1;
488 return 0;
489}
490
491int SGCosmic::mkDz_neg_v1(const double dz, const float rndm)
492{
493 const int nch = 100;
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
515 };
516 // normalize by peak value --> max = 1
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.));
520
521 if (dz_neg[index] > rndm) return 1;
522 return 0;
523}
524
525int SGCosmic::mkTanl_pos_v1(const double tanl, const float rndm)
526{
527 const int nch = 100;
528 const int width = 6;
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
549 };
550 // normalize by peak value --> max = 1
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.));
554
555 if (tanl_pos[index] > rndm) return 1;
556 return 0;
557}
558
559int SGCosmic::mkTanl_neg_v1(const double tanl, const float rndm)
560{
561 const int nch = 100;
562 const int width = 6;
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
583 };
584 // normalize by peak value --> max = 1
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.));
588
589 if (tanl_neg[index] > rndm) return 1;
590 return 0;
591}
592
593// -------- version 2 --------
594void SGCosmic::mkdist_v2(const int charge, double* parameters)
595{
596 double rand1;
597 double rand2;
598
599 for (int i = 0; i < 5; i++) {
600 int success;
601 success = 0;
602 do {
603 rand1 = gRandom->Uniform(0., 1.);
604 rand2 = gRandom->Uniform(0., 1.);
605 if (i == 0) { // dr
606 parameters[0] = rand1 * 6. - 3.; // (-3, 3)
607 if (charge == 1)
608 success = mkDr_pos_v2(parameters[0], rand2);
609 else
610 success = mkDr_neg_v2(parameters[0], rand2);
611 }
612 if (i == 1) { // phi
613 parameters[1] = - rand1 * M_PI; // (-pi, 0)
614 if (charge == 1)
615 success = mkPhi_pos_v2(parameters[1], rand2);
616 else
617 success = mkPhi_neg_v2(parameters[1], rand2);
618 }
619 if (i == 2) { // pt
620 parameters[2] = 50. * rand1; // (0., 50)
621 if (charge == 1)
622 success = mkPt_pos_v2(parameters[2], rand2);
623 else
624 success = mkPt_neg_v2(parameters[2], rand2);
625 }
626 if (i == 3) { // dz
627 parameters[3] = rand1 * 80. - 30.; // (-30, 50)
628 if (charge == 1)
629 success = mkDz_pos_v2(parameters[3], rand2);
630 else
631 success = mkDz_neg_v2(parameters[3], rand2);
632 }
633 if (i == 4) { // tanl
634 parameters[4] = rand1 * 4 - 2; // (-2, 2)
635 if (charge == 1)
636 success = mkTanl_pos_v2(parameters[4], rand2);
637 else
638 success = mkTanl_neg_v2(parameters[4], rand2);
639 }
640 } while (success == 0);
641 }
642 return;
643}
644
645int SGCosmic::mkDr_pos_v2(const double dr, const float rndm)
646{
647 const int nch = 50;
648 const int width = 6;
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.,
658 0., 0., 0., 0., 0
659 };
660 // normalize by peak value --> max = 1
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.));
664
665 if (dr_pos[index] > rndm) return 1;
666 return 0;
667}
668
669int SGCosmic::mkDr_neg_v2(const double dr, const float rndm)
670{
671 const int nch = 50;
672 const int width = 6;
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.,
682 0., 0., 0., 0., 0.
683 };
684 // normalize by peak value --> max = 1
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.));
688
689 if (dr_neg[index] > rndm) return 1;
690 return 0;
691}
692
693int SGCosmic::mkPhi_pos_v2(const double phi, const float rndm)
694{
695 const int nch = 50;
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,
706 0., 0., 0., 0., 0.
707 };
708 // normalize by peak value --> max = 1
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.));
712
713 if (phi_pos[index] > rndm) return 1;
714 return 0;
715}
716
717int SGCosmic::mkPhi_neg_v2(const double phi, const float rndm)
718{
719 const int nch = 50;
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,
730 0., 0., 0., 0., 0.
731 };
732 // normalize by peak value --> max = 1
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.));
736
737 if (phi_neg[index] > rndm) return 1;
738 return 0;
739}
740
741int SGCosmic::mkPt_pos_v2(const double pt, const float rndm)
742{
743 const int nch = 100;
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
765 };
766 // normalize by peak value --> max = 1
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);
770
771 if (pt_pos[index] > rndm) return 1;
772 return 0;
773}
774
775int SGCosmic::mkPt_neg_v2(const double pt, const float rndm)
776{
777 const int nch = 100;
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
799 };
800 // normalize by peak value --> max = 1
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);
804
805 if (pt_neg[index] > rndm) return 1;
806 return 0;
807}
808
809int SGCosmic::mkDz_pos_v2(const double dz, const float rndm)
810{
811 const int nch = 200;
812 const int width = 100;
813 double dz_pos[nch] = { 0., 0., 0., 1.000000, 0.,
814 0., 0., 0., 0., 0.,
815 0., 0., 0., 0., 0.,
816 0., 0., 0., 0., 0.,
817 0., 0., 0., 0., 0.,
818 0., 0., 0., 0., 0.,
819 0., 0., 0., 0., 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.,
851 0., 0., 0., 0., 0.,
852 0., 0., 0., 0., 0.
853 };
854 // normalize by peak value --> max = 1
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.));
858
859 if (dz_pos[index] > rndm) return 1;
860 return 0;
861}
862
863int SGCosmic::mkDz_neg_v2(const double dz, const float rndm)
864{
865 const int nch = 100;
866 //const int width = 100;
867 double dz_neg[nch] = { 0., 1.000000, 0., 0., 0.,
868 0., 0., 0., 0., 0.,
869 0., 0., 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.
887 };
888 // normalize by peak value --> max = 1
889 double max = findMax(dz_neg, nch);
890 for (int i = 0; i < nch; i++) dz_neg[i] /= max;
891 // nch is equal to width here, not needed to compute the ratio
892 int index = (int)(dz + 50.);
893 if (dz_neg[index] > rndm) return 1;
894 return 0;
895}
896
897int SGCosmic::mkTanl_pos_v2(const double tanl, const float rndm)
898{
899 const int nch = 100;
900 const int width = 5;
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
921 };
922 // normalize by peak value --> max = 1
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));
926
927 if (tanl_pos[index] > rndm) return 1;
928 return 0;
929}
930
931int SGCosmic::mkTanl_neg_v2(const double tanl, const float rndm)
932{
933 const int nch = 100;
934 const int width = 5;
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,
954 0., 0., 0., 0., 0.
955 };
956 // normalize by peak value --> max = 1
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));
960
961 if (tanl_neg[index] > rndm) return 1;
962 return 0;
963}
Helix parameter class.
Definition: Helix.h:48
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Helix.cc:197
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Definition: Helix.cc:109
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.
Definition: MCParticle.h:47
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition: MCParticle.h:49
SGCosmic()
Default constructor.
Definition: SGCosmic.cc:31
int mkPhi_pos_v2(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for positively charged particles.
Definition: SGCosmic.cc:693
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...
Definition: SGCosmic.cc:559
int mkDz_pos_v1(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for positively charged particles.
Definition: SGCosmic.cc:457
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...
Definition: SGCosmic.cc:525
int mkPt_neg_v1(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for negatively charged particl...
Definition: SGCosmic.cc:433
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,...
Definition: SGCosmic.cc:152
int mkPhi_neg_v1(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for negatively charged particles.
Definition: SGCosmic.cc:375
int mkDr_neg_v1(const double dr, const float rndm)
Generates vertex distributions in the radial direction dr by accept-reject method for negatively.
Definition: SGCosmic.cc:307
void mkdist_v2(const int charge, double *)
Generates distributions in 5-parameter space for different particle charges.
Definition: SGCosmic.cc:594
int mkPt_pos_v2(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for positively charged particl...
Definition: SGCosmic.cc:741
Parameters m_params
All relevant parameters.
Definition: SGCosmic.h:83
bool generateEvent(MCParticleGraph &graph)
Generates the next event and store the result in the given MCParticle graph.
Definition: SGCosmic.cc:74
int mkDz_neg_v1(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for negatively charged particles.
Definition: SGCosmic.cc:491
int mkPt_pos_v1(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for positively charged particl...
Definition: SGCosmic.cc:409
int mkPhi_neg_v2(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for negatively charged particles.
Definition: SGCosmic.cc:717
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...
Definition: SGCosmic.cc:669
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...
Definition: SGCosmic.cc:645
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...
Definition: SGCosmic.cc:897
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...
Definition: SGCosmic.cc:931
int mkDz_neg_v2(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for negatively charged particles.
Definition: SGCosmic.cc:863
bool setParameters(const Parameters &parameters)
Sets the parameters for generating the Particles.
Definition: SGCosmic.cc:44
void mkdist_v1(const int charge, double *)
Generates distributions in 5-parameter space for different particle charges.
Definition: SGCosmic.cc:209
int muChargeFlag(const double)
Generates the muon charge according to the positively/negatively charged muon ratio.
Definition: SGCosmic.cc:200
int mkPt_neg_v2(const double pt, const float rndm)
Generates transverse momentum pt distributions by accept-reject method for negatively charged particl...
Definition: SGCosmic.cc:775
int mkDz_pos_v2(const double dz, const float rndm)
Generates z vertex dz distributions by accept-reject method for positively charged particles.
Definition: SGCosmic.cc:809
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...
Definition: SGCosmic.cc:273
int mkPhi_pos_v1(const double phi, const float rndm)
Generates azimuthal angle phi distributions by accept-reject method for positively charged particles.
Definition: SGCosmic.cc:341
double findMax(const double *dim, const int num)
Finds maximum value in an array.
Definition: SGCosmic.cc:65
static const double T
[tesla]
Definition: Unit.h:120
static void getField(const double *pos, double *field)
return the magnetic field at a given position.
Definition: BFieldManager.h:91
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.
STL namespace.
Struct to keep all necessary parameters for the cosmic generator.
Definition: SGCosmic.h:30
int ipRequirement
Restrict the vertex to IP or not (default)
Definition: SGCosmic.h:38
double ipdr
Vertex restriction in the radial direction.
Definition: SGCosmic.h:42
double cylindricalR
Cylindrical radius of generation.
Definition: SGCosmic.h:54
double ptmin
Minimum value of the transverse momentum.
Definition: SGCosmic.h:50
int level
Generator version: level 1 (default) or 2.
Definition: SGCosmic.h:34
double ipdz
Vertex restriction in z direction.
Definition: SGCosmic.h:46