Belle II Software  release-08-01-10
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 
20 using namespace std;
21 using namespace Belle2;
22 
23 //#define DEBUG
24 #ifdef DEBUG
25 ofstream ofs_cosmic1g("cosmic1g.out");
26 ofstream ofs_cosmic1c("cosmic1c.out");
27 ofstream ofs_cosmic2g("cosmic2g.out");
28 ofstream ofs_cosmic2c("cosmic2c.out");
29 #endif
30 
31 SGCosmic::SGCosmic()
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 
44 bool SGCosmic::setParameters(const Parameters& p)
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
65 double 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 
74 bool SGCosmic::generateEvent(MCParticleGraph& graph)
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 
83  MCParticleGraph::GraphParticle& p = graph.addParticle();
84  p.setStatus(MCParticle::c_PrimaryParticle);
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);
145  p.addStatus(MCParticle::c_StableInGenerator);
146 
147  return true;
148  }
149  return false;
150 }
151 
152 void 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 
200 int 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 --------
209 void 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 
273 int 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 
307 int 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 
341 int 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 
375 int 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 
409 int 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 
433 int 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 
457 int 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 
491 int 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 
525 int 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 
559 int 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 --------
594 void 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 
645 int 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 
669 int 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 
693 int 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 
717 int 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 
741 int 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 
775 int 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 
809 int 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 
863 int 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 
897 int 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 
931 int 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 }
Belle2::CDC::Helix
Helix parameter class.
Definition: Helix.h:48
Belle2::CDC::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Helix.cc:197
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:75
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:37
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:294
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
Belle2::SGCosmic::Parameters
Struct to keep all necessary parameters for the cosmic generator.
Definition: SGCosmic.h:30
Belle2::SGCosmic::Parameters::ipRequirement
int ipRequirement
Restrict the vertex to IP or not (default)
Definition: SGCosmic.h:38
Belle2::SGCosmic::Parameters::ipdr
double ipdr
Vertex restriction in the radial direction.
Definition: SGCosmic.h:42
Belle2::SGCosmic::Parameters::cylindricalR
double cylindricalR
Cylindrical radius of generation.
Definition: SGCosmic.h:54
Belle2::SGCosmic::Parameters::ptmin
double ptmin
Minimum value of the transverse momentum.
Definition: SGCosmic.h:50
Belle2::SGCosmic::Parameters::level
int level
Generator version: level 1 (default) or 2.
Definition: SGCosmic.h:34
Belle2::SGCosmic::Parameters::ipdz
double ipdz
Vertex restriction in z direction.
Definition: SGCosmic.h:46