Belle II Software  release-05-01-25
SGCosmic.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Based on the Belle single track generator by KUNIYA Toshio *
7  * Contributors: Sergey Yashchenko *
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <framework/logging/Logger.h>
13 #include <generators/cosmics/SGCosmic.h>
14 
15 #include <cmath>
16 #include <fstream>
17 
18 #include <framework/dataobjects/Helix.h>
19 #include <framework/geometry/BFieldManager.h>
20 #include <TRandom.h>
21 #include <TVector3.h>
22 
23 using namespace std;
24 using namespace Belle2;
25 
26 //#define DEBUG
27 #ifdef DEBUG
28 ofstream ofs_cosmic1g("cosmic1g.out");
29 ofstream ofs_cosmic1c("cosmic1c.out");
30 ofstream ofs_cosmic2g("cosmic2g.out");
31 ofstream ofs_cosmic2c("cosmic2c.out");
32 #endif
33 
34 SGCosmic::SGCosmic()
35 {
36  Parameters m_parameters;
37  m_parameters.level = 1;
38  m_parameters.ipRequirement = 0;
39  m_parameters.ipdr = 3.; // Only relevant for ipRequirement = 1
40  m_parameters.ipdz = 3.; // Only relevant for ipRequirement = 1
41  m_parameters.ptmin = 0.7;
42  m_parameters.cylindricalR = 125.0;
43  setParameters(m_parameters);
44 }
45 
46 
47 bool SGCosmic::setParameters(const Parameters& p)
48 {
49  // Sanity checks
50  bool ok(true);
51 
52  // Check that we have correct generator level
53  if (p.level != 1 && p.level != 2) {
54  B2ERROR("Wrong generator flag");
55  ok = false;
56  }
57 
58  // If everything is ok, set the new parameters, else return false
59  if (ok) {
60  m_params = p;
61  return true;
62  }
63  return false;
64 }
65 
66 
67 // Utility from the Belle cosmics version
68 double SGCosmic::findMax(const double* dim, const int num)
69 {
70  double max = 0.;
71  for (int i = 0; i < num; i++) {
72  if (max < dim[i]) max = dim[i];
73  }
74  return max;
75 }
76 
77 bool SGCosmic::generateEvent(MCParticleGraph& graph)
78 {
79  int charge = 0; // Initialization of charge
80  double dr, pt, phi, dz, tanl;
81 
82  while (true) {
83 
84  genCosmic(m_params.level, charge, dr, phi, pt, dz, tanl);
85 
86  MCParticleGraph::GraphParticle& p = graph.addParticle();
87  p.setStatus(MCParticle::c_PrimaryParticle);
88 
89  if (1 == charge) {
90  p.setPDG(-13);
91  } else if (-1 == charge) {
92  p.setPDG(13);
93  } else {
94  continue;
95  }
96  p.setMassFromPDG();
97  p.setFirstDaughter(0);
98  p.setLastDaughter(0);
99 
100 #ifdef DEBUG
101  if (1 == m_params.level) ofs_cosmic1g << charge << " " << dr << " " << phi << " "
102  << pt << " " << dz << " " << tanl << endl;
103  else if (2 == m_params.level) ofs_cosmic2g << charge << " " << dr << " " << phi << " "
104  << pt << " " << dz << " " << tanl << endl;
105 #endif
106 
107  // Simulate helix parameter at perigee
108  float bz = BFieldManager::getField(0, 0, 0).Z() / Unit::T; // Magnetic field
109  const float EPS = 0.0001; // Avoid using zero magnetic field
110  if (fabs(bz) < EPS) bz = EPS; // Set the value of the magnetic field to a small number
111  float d0, phi0, omega, z0, tanLambda; // The definition is the same as in the Helix class
112  d0 = dr;
113  phi0 = phi;
114  omega = (double)charge / (pt * Helix::getAlpha(bz));
115  z0 = dz;
116  tanLambda = tanl;
117  Helix CosmicMCHelix(d0, phi0, omega, z0, tanLambda);
118 
119  const float cylindricalR = m_params.cylindricalR; // radius (cm) of generation
120  float arcLength;
121  // Get arc length at ToP radius
122  arcLength = CosmicMCHelix.getArcLength2DAtCylindricalR(cylindricalR);
123  if (isnan(arcLength)) continue;
124 
125  // Calculate coordinates and momentum at ToP radius
126  TVector3 vector;
127  vector = CosmicMCHelix.getPositionAtArcLength2D(arcLength);
128  double vx = vector[0];
129  double vy = vector[1];
130  double vz = vector[2];
131  vector = CosmicMCHelix.getMomentumAtArcLength2D(arcLength, bz);
132  double px = - vector[0];
133  double py = - vector[1];
134  double pz = - vector[2];
135  double m = p.getMass();
136  double e = sqrt(px * px + py * py + pz * pz + m * m);
137 
138 #ifdef DEBUG
139  if (1 == m_params.level) ofs_cosmic1c << vx << " " << vy << " " << vz << " "
140  << px << " " << py << " " << pz << endl;
141  else if (2 == m_params.level) ofs_cosmic2c << vx << " " << vy << " " << vz << " "
142  << px << " " << py << " " << pz << endl;
143 #endif
144 
145  p.setMomentum(px, py, pz);
146  p.setEnergy(e);
147  p.setProductionVertex(vx, vy, vz);
148  p.addStatus(MCParticle::c_StableInGenerator);
149 
150  return true;
151  }
152  return false;
153 }
154 
155 void SGCosmic::genCosmic(const int level, int& charge,
156  double& dr,
157  double& phi,
158  double& pt,
159  double& dz,
160  double& tanl)
161 {
162  // generator version 1
163  if (level == 1) {
164  double muPluse = 361917;
165  double muMinus = 285667;
166  double fraction = muPluse / (muPluse + muMinus);
167 
168  // mu+ or mu- --> fraction = f(mu+)
169  charge = muChargeFlag(fraction);
170 
171  // make distribution
172  double parameters[5]; // dr,phi,pt(not Kappa),dz,tanl
173  mkdist_v1(charge, parameters);
174  dr = parameters[0];
175  phi = parameters[1];
176  pt = parameters[2]; // not Kappa
177  dz = parameters[3];
178  tanl = parameters[4];
179  }
180 
181  // generator version 2
182  if (level == 2) {
183  double muPluse = 143391;
184  double muMinus = 143302;
185  double fraction = muPluse / (muPluse + muMinus);
186 
187  // mu+ or mu- --> fraction = f(mu+)
188  charge = muChargeFlag(fraction);
189 
190  // make distribution
191  double parameters[5]; // dr,phi,pt(not Kappa),dz,tanl
192  mkdist_v2(charge, parameters);
193  dr = parameters[0];
194  phi = parameters[1];
195  pt = parameters[2]; // not Kappa
196  dz = parameters[3];
197  tanl = parameters[4];
198  }
199 
200  return;
201 }
202 
203 int SGCosmic::muChargeFlag(const double fraction)
204 {
205  double rand1;
206  rand1 = gRandom->Uniform(0., 1.);
207  if ((double)rand1 > fraction) return -1;
208  else return 1;
209 }
210 
211 // -------- version 1 --------
212 void SGCosmic::mkdist_v1(const int charge, double* parameters)
213 {
214  double rand1;
215  double rand2;
216 
217  for (int i = 0; i < 5; i++) {
218  int success;
219  success = 0;
220  do {
221  rand1 = gRandom->Uniform(0., 1.);
222  rand2 = gRandom->Uniform(0., 1.);
223  if (i == 0) { // dr
224  if (m_params.ipRequirement == 0) {
225  parameters[0] = rand1 * 40. - 20.; // (-20, 20)
226  } else {
227  // IP requirement
228  double widthdr = m_params.ipdr * 2;
229  parameters[0] = rand1 * widthdr - m_params.ipdr; // (-ipdr, ipdr)
230  }
231  if (charge == 1)
232  success = mkDr_pos_v1(parameters[0], rand2);
233  else
234  success = mkDr_neg_v1(parameters[0], rand2);
235  }
236  if (i == 1) { // phi
237  parameters[1] = rand1 * M_PI; // (0, pi)
238  if (charge == 1)
239  success = mkPhi_pos_v1(parameters[1], rand2);
240  else
241  success = mkPhi_neg_v1(parameters[1], rand2);
242  }
243  if (i == 2) { // pt
244  parameters[2] = (30.0 - m_params.ptmin) * rand1 + m_params.ptmin; // (ptmin, 30)
245  if (charge == 1)
246  success = mkPt_pos_v1(parameters[2], rand2);
247  else
248  success = mkPt_neg_v1(parameters[2], rand2);
249  }
250  if (i == 3) { // dz
251  if (m_params.ipRequirement == 0) {
252  parameters[3] = rand1 * 150. - 50.; //(-50, 100)
253  } else {
254  // IP requirement
255  double widthdz = m_params.ipdz * 2;
256  parameters[3] = rand1 * widthdz - m_params.ipdz; //(-ipdz, ipdz)
257  }
258  if (charge == 1)
259  success = mkDz_pos_v1(parameters[3], rand2);
260  else
261  success = mkDz_neg_v1(parameters[3], rand2);
262  }
263  if (i == 4) { // tanl
264  parameters[4] = (rand1 - 0.5) * 6; // (-3, 3)
265  if (charge == 1)
266  success = mkTanl_pos_v1(parameters[4], rand2);
267  else
268  success = mkTanl_neg_v1(parameters[4], rand2);
269  }
270  } while (success == 0);
271  }
272  return;
273 }
274 
275 
276 int SGCosmic::mkDr_pos_v1(const double dr, const float rndm)
277 {
278  const int nch = 100;
279  const int width = 40;
280  double dr_pos[nch] = {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  0.0, 0.0, 0.0, 0.0, 0.0,
284  0.0, 0.0, 0.0, 0.0, 0.0,
285  0.0, 0.0, 0.0, 0.0, 0.0,
286  3284.0, 4110.0, 4545.0, 5140.0, 5601.0,
287  6304.0, 6840.0, 7142.0, 7714.0, 8166.0,
288  8795.0, 9309.0, 9600.0, 9759.0, 9892.0,
289  9995.0, 10198.0, 10379.0, 10428.0, 10646.0,
290  10791.0, 11002.0, 11172.0, 11332.0, 11334.0,
291  11381.0, 11226.0, 11241.0, 10958.0, 10840.0,
292  10619.0, 10246.0, 10135.0, 9640.0, 9477.0,
293  9215.0, 8689.0, 8610.0, 8150.0, 7679.0,
294  359.0, 2.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  0.0, 0.0, 0.0, 0.0, 0.0,
298  0.0, 0.0, 0.0, 0.0, 0.0,
299  0.0, 0.0, 0.0, 0.0, 0.0
300  };
301  // normalize by peak value --> max = 1
302  double max = findMax(dr_pos, nch);
303  for (int i = 0; i < nch; i++) dr_pos[i] /= max;
304  int index = (int)((double)nch / (double)width * (dr + 20.));
305 
306  if (dr_pos[index] > rndm) return 1;
307  return 0;
308 }
309 
310 int SGCosmic::mkDr_neg_v1(const double dr, const float rndm)
311 {
312  const int nch = 100;
313  const int width = 40;
314  double dr_neg[nch] = { 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  0.0, 0.0, 0.0, 0.0, 0.0,
318  0.0, 0.0, 0.0, 0.0, 0.0,
319  0.0, 0.0, 0.0, 0.0, 0.0,
320  2735.0, 3081.0, 3568.0, 3837.0, 4445.0,
321  4938.0, 5338.0, 5701.0, 5985.0, 6404.0,
322  7022.0, 7375.0, 7438.0, 7662.0, 7833.0,
323  7951.0, 8166.0, 8186.0, 8101.0, 8159.0,
324  8541.0, 8687.0, 8855.0, 8908.0, 8998.0,
325  9014.0, 8962.0, 8915.0, 8833.0, 8698.0,
326  8399.0, 8285.0, 7875.0, 7726.0, 7346.0,
327  7058.0, 7074.0, 6887.0, 6546.0, 5910.0,
328  263.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  0.0, 0.0, 0.0, 0.0, 0.0,
332  0.0, 0.0, 0.0, 0.0, 0.0,
333  0.0, 0.0, 0.0, 0.0, 0.0
334  };
335  // normalize by peak value --> max = 1
336  double max = findMax(dr_neg, nch);
337  for (int i = 0; i < nch; i++) dr_neg[i] /= max;
338  int index = (int)((double)nch / (double)width * (dr + 20.));
339 
340  if (dr_neg[index] > rndm) return 1;
341  return 0;
342 }
343 
344 int SGCosmic::mkPhi_pos_v1(const double phi, const float rndm)
345 {
346  const int nch = 100;
347  const double width = M_PI;
348  double phi_pos[nch] = { 95.0, 194.0, 287.0, 433.0, 636.0,
349  839.0, 1190.0, 1486.0, 1795.0, 2090.0,
350  2192.0, 2437.0, 2831.0, 3114.0, 3407.0,
351  3589.0, 3706.0, 3796.0, 4125.0, 4292.0,
352  4650.0, 4898.0, 4877.0, 5003.0, 5037.0,
353  5409.0, 5488.0, 5704.0, 5733.0, 5547.0,
354  5704.0, 5769.0, 5919.0, 5993.0, 5976.0,
355  5965.0, 6011.0, 6036.0, 6069.0, 6450.0,
356  6260.0, 6103.0, 6252.0, 6277.0, 6219.0,
357  6364.0, 6287.0, 6298.0, 6033.0, 5894.0,
358  5883.0, 6036.0, 6032.0, 5894.0, 5748.0,
359  5229.0, 5382.0, 5374.0, 5504.0, 5283.0,
360  5074.0, 4700.0, 4651.0, 4847.0, 4698.0,
361  4446.0, 4410.0, 4104.0, 3870.0, 3790.0,
362  3826.0, 3750.0, 3488.0, 3220.0, 2973.0,
363  2862.0, 2718.0, 2616.0, 2400.0, 2168.0,
364  1927.0, 1860.0, 1708.0, 1667.0, 1334.0,
365  1162.0, 965.0, 787.0, 629.0, 554.0,
366  384.0, 333.0, 243.0, 159.0, 139.0,
367  88.0, 59.0, 50.0, 27.0, 9.0
368  };
369  // normalize by peak value --> max = 1
370  double max = findMax(phi_pos, nch);
371  for (int i = 0; i < nch; i++) phi_pos[i] /= max;
372  int index = (int)((double)nch / width * phi);
373 
374  if (phi_pos[index] > rndm) return 1;
375  return 0;
376 }
377 
378 int SGCosmic::mkPhi_neg_v1(const double phi, const float rndm)
379 {
380  const int nch = 100;
381  const double width = M_PI;
382  double phi_neg[nch] = { 351.0, 483.0, 520.0, 506.0, 541.0,
383  497.0, 475.0, 503.0, 495.0, 532.0,
384  594.0, 638.0, 787.0, 904.0, 1094.0,
385  1186.0, 1250.0, 1349.0, 1494.0, 1658.0,
386  1807.0, 1929.0, 2031.0, 2077.0, 2292.0,
387  2402.0, 2592.0, 2790.0, 2988.0, 2908.0,
388  2968.0, 3107.0, 3264.0, 3534.0, 3505.0,
389  3582.0, 3608.0, 3649.0, 3772.0, 3827.0,
390  4023.0, 3960.0, 3956.0, 3965.0, 4095.0,
391  4392.0, 4504.0, 4428.0, 4364.0, 4510.0,
392  4633.0, 4767.0, 4968.0, 4866.0, 4889.0,
393  4753.0, 4862.0, 4928.0, 5015.0, 4958.0,
394  4960.0, 4921.0, 4763.0, 4838.0, 4841.0,
395  4779.0, 4627.0, 4679.0, 4478.0, 4528.0,
396  4470.0, 4474.0, 4199.0, 4145.0, 3896.0,
397  3882.0, 3812.0, 3839.0, 3466.0, 3274.0,
398  3067.0, 3043.0, 2914.0, 2776.0, 2677.0,
399  2425.0, 2195.0, 1892.0, 1777.0, 1685.0,
400  1413.0, 1296.0, 1123.0, 1017.0, 795.0,
401  739.0, 665.0, 493.0, 395.0, 251.0
402  };
403  // normalize by peak value --> max = 1
404  double max = findMax(phi_neg, nch);
405  for (int i = 0; i < nch; i++) phi_neg[i] /= max;
406  int index = (int)((double)nch / width * phi);
407 
408  if (phi_neg[index] > rndm) return 1;
409  return 0;
410 }
411 
412 int SGCosmic::mkPt_pos_v1(const double pt, const float rndm)
413 {
414  const int nch = 50;
415  const int width = 30;
416  double pt_pos[nch] = { 47.0, 22547.0, 28678.0, 26159.0, 25145.0,
417  23668.0, 21581.0, 19197.0, 16847.0, 14902.0,
418  13020.0, 11836.0, 10397.0, 9414.0, 8356.0,
419  7492.0, 6640.0, 6173.0, 5426.0, 5023.0,
420  4612.0, 4340.0, 3860.0, 3656.0, 3222.0,
421  3042.0, 2846.0, 2655.0, 2422.0, 2288.0,
422  2172.0, 1984.0, 1841.0, 1745.0, 1599.0,
423  1495.0, 1494.0, 1366.0, 1305.0, 1186.0,
424  1153.0, 1063.0, 1051.0, 996.0, 952.0,
425  906.0, 883.0, 783.0, 733.0, 682.0
426  };
427  // normalize by peak value --> max = 1
428  double max = findMax(pt_pos, nch);
429  for (int i = 0; i < nch; i++) pt_pos[i] /= max;
430  int index = (int)((double)nch / width * pt);
431 
432  if (pt_pos[index] > rndm) return 1;
433  return 0;
434 }
435 
436 int SGCosmic::mkPt_neg_v1(const double pt, const float rndm)
437 {
438  const int nch = 50;
439  const int width = 30;
440  double pt_neg[nch] = { 538.0, 614.0, 600.0, 626.0, 648.0,
441  729.0, 716.0, 760.0, 828.0, 882.0,
442  914.0, 922.0, 1009.0, 1105.0, 1113.0,
443  1246.0, 1305.0, 1345.0, 1523.0, 1578.0,
444  1649.0, 1856.0, 1968.0, 2126.0, 2319.0,
445  2526.0, 2716.0, 3005.0, 3317.0, 3427.0,
446  3869.0, 4287.0, 4727.0, 5184.0, 5882.0,
447  6299.0, 7191.0, 8074.0, 9017.0, 9968.0,
448  11425.0, 12821.0, 14601.0, 16351.0, 17885.0,
449  19838.0, 20766.0, 23697.0, 22861.0, 116.0
450  };
451  // normalize by peak value --> max = 1
452  double max = findMax(pt_neg, nch);
453  for (int i = 0; i < nch; i++) pt_neg[i] /= max;
454  int index = (int)((double)nch / width * (width - pt));
455 
456  if (pt_neg[index] > rndm) return 1;
457  return 0;
458 }
459 
460 int SGCosmic::mkDz_pos_v1(const double dz, const float rndm)
461 {
462  const int nch = 100;
463  const int width = 250;
464  double dz_pos[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
465  0.0, 0.0, 0.0, 0.0, 0.0,
466  0.0, 0.0, 0.0, 0.0, 0.0,
467  0.0, 0.0, 0.0, 0.0, 0.0,
468  1294.0, 1796.0, 2018.0, 2146.0, 2280.0,
469  2422.0, 2601.0, 2913.0, 3076.0, 3564.0,
470  5147.0, 6284.0, 6958.0, 7133.0, 7465.0,
471  7481.0, 7878.0, 7947.0, 8297.0, 8589.0,
472  8366.0, 8540.0, 8717.0, 8825.0, 8884.0,
473  8985.0, 9131.0, 9105.0, 9279.0, 9333.0,
474  9329.0, 9501.0, 9442.0, 9484.0, 9471.0,
475  9322.0, 9201.0, 8344.0, 6855.0, 5879.0,
476  5095.0, 5012.0, 5044.0, 4967.0, 4717.0,
477  4928.0, 4917.0, 4762.0, 4579.0, 4556.0,
478  4461.0, 4366.0, 4324.0, 4185.0, 4185.0,
479  4085.0, 4019.0, 3974.0, 3748.0, 2739.0,
480  0.0, 0.0, 0.0, 0.0, 0.0,
481  0.0, 0.0, 0.0, 0.0, 0.0,
482  0.0, 0.0, 0.0, 0.0, 0.0,
483  0.0, 0.0, 0.0, 0.0, 0.0
484  };
485  // normalize by peak value --> max = 1
486  double max = findMax(dz_pos, nch);
487  for (int i = 0; i < nch; i++) dz_pos[i] /= max;
488  int index = (int)((double)nch / (double)width * (dz + 100.));
489 
490  if (dz_pos[index] > rndm) return 1;
491  return 0;
492 }
493 
494 int SGCosmic::mkDz_neg_v1(const double dz, const float rndm)
495 {
496  const int nch = 100;
497  const int width = 250;
498  double dz_neg[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
499  0.0, 0.0, 0.0, 0.0, 0.0,
500  0.0, 0.0, 0.0, 0.0, 0.0,
501  0.0, 0.0, 0.0, 0.0, 0.0,
502  967.0, 1357.0, 1507.0, 1618.0, 1736.0,
503  1926.0, 2090.0, 2182.0, 2390.0, 2748.0,
504  4261.0, 5004.0, 5605.0, 5819.0, 5932.0,
505  5946.0, 6214.0, 6281.0, 6521.0, 6626.0,
506  6743.0, 6830.0, 6912.0, 7217.0, 7089.0,
507  7194.0, 7340.0, 7283.0, 7343.0, 7357.0,
508  7141.0, 7545.0, 7476.0, 7608.0, 7415.0,
509  7369.0, 7365.0, 6827.0, 5364.0, 4520.0,
510  4052.0, 3994.0, 3885.0, 3820.0, 3743.0,
511  3747.0, 3543.0, 3604.0, 3580.0, 3479.0,
512  3467.0, 3395.0, 3409.0, 3365.0, 3231.0,
513  3247.0, 3192.0, 3074.0, 2995.0, 2215.0,
514  0.0, 0.0, 0.0, 0.0, 0.0,
515  0.0, 0.0, 0.0, 0.0, 0.0,
516  0.0, 0.0, 0.0, 0.0, 0.0,
517  0.0, 0.0, 0.0, 0.0, 0.0
518  };
519  // normalize by peak value --> max = 1
520  double max = findMax(dz_neg, nch);
521  for (int i = 0; i < nch; i++) dz_neg[i] /= max;
522  int index = (int)((double)nch / (double)width * (dz + 100.));
523 
524  if (dz_neg[index] > rndm) return 1;
525  return 0;
526 }
527 
528 int SGCosmic::mkTanl_pos_v1(const double tanl, const float rndm)
529 {
530  const int nch = 100;
531  const int width = 6;
532  double tanl_pos[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
533  0.0, 0.0, 0.0, 0.0, 0.0,
534  0.0, 0.0, 0.0, 0.0, 0.0,
535  0.0, 0.0, 0.0, 0.0, 0.0,
536  0.0, 0.0, 0.0, 0.0, 1.0,
537  1.0, 1.0, 2.0, 6.0, 15.0,
538  32.0, 30.0, 103.0, 375.0, 906.0,
539  1465.0, 2088.0, 2858.0, 4007.0, 5324.0,
540  7025.0, 8801.0, 11256.0, 13767.0, 16230.0,
541  18402.0, 20313.0, 21764.0, 22369.0, 23219.0,
542  22975.0, 23058.0, 22121.0, 20633.0, 18402.0,
543  16315.0, 14098.0, 11450.0, 8844.0, 6793.0,
544  5207.0, 3937.0, 2850.0, 1990.0, 1464.0,
545  921.0, 384.0, 83.0, 23.0, 18.0,
546  6.0, 3.0, 2.0, 2.0, 4.0,
547  1.0, 0.0, 1.0, 0.0, 0.0,
548  0.0, 0.0, 0.0, 0.0, 0.0,
549  0.0, 0.0, 0.0, 0.0, 0.0,
550  0.0, 0.0, 0.0, 0.0, 0.0,
551  0.0, 0.0, 0.0, 0.0, 0.0
552  };
553  // normalize by peak value --> max = 1
554  double max = findMax(tanl_pos, nch);
555  for (int i = 0; i < nch; i++) tanl_pos[i] /= max;
556  int index = (int)((double)nch / (double)width * (tanl + 3.));
557 
558  if (tanl_pos[index] > rndm) return 1;
559  return 0;
560 }
561 
562 int SGCosmic::mkTanl_neg_v1(const double tanl, const float rndm)
563 {
564  const int nch = 100;
565  const int width = 6;
566  double tanl_neg[nch] = { 0.0, 0.0, 0.0, 0.0, 0.0,
567  0.0, 0.0, 0.0, 0.0, 0.0,
568  0.0, 0.0, 0.0, 0.0, 0.0,
569  0.0, 0.0, 0.0, 0.0, 0.0,
570  0.0, 0.0, 0.0, 0.0, 0.0,
571  1.0, 1.0, 0.0, 7.0, 10.0,
572  11.0, 24.0, 94.0, 335.0, 760.0,
573  1210.0, 1713.0, 2233.0, 3242.0, 4236.0,
574  5530.0, 7004.0, 8768.0, 10923.0, 12809.0,
575  14457.0, 16043.0, 17321.0, 18064.0, 18237.0,
576  18129.0, 17823.0, 17077.0, 16272.0, 14539.0,
577  13015.0, 11031.0, 9000.0, 7219.0, 5411.0,
578  4163.0, 3047.0, 2120.0, 1551.0, 1215.0,
579  672.0, 266.0, 73.0, 23.0, 5.0,
580  4.0, 8.0, 5.0, 3.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  0.0, 0.0, 0.0, 0.0, 0.0,
584  0.0, 0.0, 0.0, 0.0, 0.0,
585  0.0, 0.0, 0.0, 0.0, 0.0
586  };
587  // normalize by peak value --> max = 1
588  double max = findMax(tanl_neg, nch);
589  for (int i = 0; i < nch; i++) tanl_neg[i] /= max;
590  int index = (int)((double)nch / (double)width * (tanl + 3.));
591 
592  if (tanl_neg[index] > rndm) return 1;
593  return 0;
594 }
595 
596 // -------- version 2 --------
597 void SGCosmic::mkdist_v2(const int charge, double* parameters)
598 {
599  double rand1;
600  double rand2;
601 
602  for (int i = 0; i < 5; i++) {
603  int success;
604  success = 0;
605  do {
606  rand1 = gRandom->Uniform(0., 1.);
607  rand2 = gRandom->Uniform(0., 1.);
608  if (i == 0) { // dr
609  parameters[0] = rand1 * 6. - 3.; // (-3, 3)
610  if (charge == 1)
611  success = mkDr_pos_v2(parameters[0], rand2);
612  else
613  success = mkDr_neg_v2(parameters[0], rand2);
614  }
615  if (i == 1) { // phi
616  parameters[1] = rand1 * M_PI; // (0, pi)
617  if (charge == 1)
618  success = mkPhi_pos_v2(parameters[1], rand2);
619  else
620  success = mkPhi_neg_v2(parameters[1], rand2);
621  }
622  if (i == 2) { // pt
623  parameters[2] = 50. * rand1; // (0., 50)
624  if (charge == 1)
625  success = mkPt_pos_v2(parameters[2], rand2);
626  else
627  success = mkPt_neg_v2(parameters[2], rand2);
628  }
629  if (i == 3) { // dz
630  parameters[3] = rand1 * 80. - 30.; // (-30, 50)
631  if (charge == 1)
632  success = mkDz_pos_v2(parameters[3], rand2);
633  else
634  success = mkDz_neg_v2(parameters[3], rand2);
635  }
636  if (i == 4) { // tanl
637  parameters[4] = rand1 * 4 - 2; // (-2, 2)
638  if (charge == 1)
639  success = mkTanl_pos_v2(parameters[4], rand2);
640  else
641  success = mkTanl_neg_v2(parameters[4], rand2);
642  }
643  } while (success == 0);
644  }
645  return;
646 }
647 
648 int SGCosmic::mkDr_pos_v2(const double dr, const float rndm)
649 {
650  const int nch = 50;
651  const int width = 6;
652  double dr_pos[nch] = {0., 0., 0., 0., 0.,
653  0., 2403.000, 3562.000, 3685.000, 3619.000,
654  3667.000, 3724.000, 3823.000, 3824.000, 3881.000,
655  3764.000, 3867.000, 3918.000, 3964.000, 3916.000,
656  4002.000, 4001.000, 3980.000, 3848.000, 4220.000,
657  4206.000, 3856.000, 3920.000, 3857.000, 3751.000,
658  3753.000, 3898.000, 3898.000, 3822.000, 3872.000,
659  3905.000, 3896.000, 3874.000, 3767.000, 3734.000,
660  3744.000, 3666.000, 3757.000, 2458.000, 0.,
661  0., 0., 0., 0., 0
662  };
663  // normalize by peak value --> max = 1
664  double max = findMax(dr_pos, nch);
665  for (int i = 0; i < nch; i++) dr_pos[i] /= max;
666  int index = (int)((double)nch / (double)width * (dr + 3.));
667 
668  if (dr_pos[index] > rndm) return 1;
669  return 0;
670 }
671 
672 int SGCosmic::mkDr_neg_v2(const double dr, const float rndm)
673 {
674  const int nch = 50;
675  const int width = 6;
676  double dr_neg[nch] = { 0., 0., 0., 0., 0.,
677  0., 2508.000, 3752.000, 3713.000, 3695.000,
678  3758.000, 3727.000, 3847.000, 3887.000, 3896.000,
679  3862.000, 3858.000, 3905.000, 3894.000, 3806.000,
680  3736.000, 3797.000, 3921.000, 3922.000, 4203.000,
681  4167.000, 3927.000, 3952.000, 3983.000, 3982.000,
682  3914.000, 3962.000, 3983.000, 3817.000, 3831.000,
683  3811.000, 3837.000, 3852.000, 3728.000, 3684.000,
684  3629.000, 3663.000, 3588.000, 2394.000, 0.,
685  0., 0., 0., 0., 0.
686  };
687  // normalize by peak value --> max = 1
688  double max = findMax(dr_neg, nch);
689  for (int i = 0; i < nch; i++) dr_neg[i] /= max;
690  int index = (int)((double)nch / (double)width * (dr + 3.));
691 
692  if (dr_neg[index] > rndm) return 1;
693  return 0;
694 }
695 
696 int SGCosmic::mkPhi_pos_v2(const double phi, const float rndm)
697 {
698  const int nch = 50;
699  const double width = 8;
700  double phi_pos[nch] = { 0., 0., 0., 0., 0.,
701  1073.000, 1911.000, 2905.000, 4165.000, 5205.000,
702  6062.000, 6955.000, 7382.000, 7262.000, 7142.000,
703  6548.000, 5759.000, 5015.000, 4399.000, 3392.000,
704  2365.000, 1533.000, 801.0000, 245.0000, 25.00000,
705  39.00000, 119.0000, 398.0000, 841.0000, 1487.000,
706  2331.000, 3138.000, 3726.000, 4430.000, 5106.000,
707  5548.000, 5766.000, 5699.000, 5408.000, 5011.000,
708  4439.000, 3732.000, 2883.000, 2034.000, 1023.000,
709  0., 0., 0., 0., 0.
710  };
711  // normalize by peak value --> max = 1
712  double max = findMax(phi_pos, nch);
713  for (int i = 0; i < nch; i++) phi_pos[i] /= max;
714  int index = (int)(nch / width * (phi + 4.));
715 
716  if (phi_pos[index] > rndm) return 1;
717  return 0;
718 }
719 
720 int SGCosmic::mkPhi_neg_v2(const double phi, const float rndm)
721 {
722  const int nch = 50;
723  const double width = 8;
724  double phi_neg[nch] = { 0., 0., 0., 0., 0.,
725  16.00000, 73.00000, 262.0000, 677.0000, 1252.000,
726  1973.000, 2896.000, 3520.000, 4173.000, 4876.000,
727  5389.000, 5730.000, 5798.000, 5553.000, 5120.000,
728  4676.000, 4036.000, 3185.000, 2320.000, 1688.000,
729  1741.000, 2164.000, 3325.000, 4588.000, 5522.000,
730  6427.000, 7159.000, 7385.000, 7282.000, 6980.000,
731  6157.000, 5496.000, 4854.000, 4066.000, 3042.000,
732  2032.000, 1236.000, 586.0000, 123.0000, 13.00000,
733  0., 0., 0., 0., 0.
734  };
735  // normalize by peak value --> max = 1
736  double max = findMax(phi_neg, nch);
737  for (int i = 0; i < nch; i++) phi_neg[i] /= max;
738  int index = (int)(nch / width * (phi + 4.));
739 
740  if (phi_neg[index] > rndm) return 1;
741  return 0;
742 }
743 
744 int SGCosmic::mkPt_pos_v2(const double pt, const float rndm)
745 {
746  const int nch = 100;
747  const int width = 50;
748  double pt_pos[nch] = { 3813.000, 14309.00, 14985.00, 13103.00, 10807.00,
749  9118.000, 7793.000, 6718.000, 5788.000, 5177.000,
750  4458.000, 3889.000, 3480.000, 3191.000, 2798.000,
751  2484.000, 2245.000, 2053.000, 1883.000, 1691.000,
752  1481.000, 1416.000, 1311.000, 1207.000, 1171.000,
753  1017.000, 953.0000, 886.0000, 790.0000, 816.0000,
754  712.0000, 685.0000, 627.0000, 577.0000, 526.0000,
755  502.0000, 442.0000, 453.0000, 415.0000, 351.0000,
756  340.0000, 358.0000, 305.0000, 330.0000, 278.0000,
757  272.0000, 250.0000, 273.0000, 243.0000, 209.0000,
758  203.0000, 176.0000, 212.0000, 182.0000, 153.0000,
759  168.0000, 155.0000, 170.0000, 129.0000, 117.0000,
760  120.0000, 114.0000, 126.0000, 123.0000, 106.0000,
761  90.00000, 97.00000, 103.0000, 82.00000, 87.00000,
762  76.00000, 104.0000, 75.00000, 82.00000, 67.00000,
763  64.00000, 63.00000, 75.00000, 65.00000, 62.00000,
764  54.00000, 86.00000, 44.00000, 44.00000, 53.00000,
765  54.00000, 53.00000, 51.00000, 49.00000, 39.00000,
766  37.00000, 42.00000, 40.00000, 40.00000, 24.00000,
767  38.00000, 39.00000, 30.00000, 28.00000, 32.00000
768  };
769  // normalize by peak value --> max = 1
770  double max = findMax(pt_pos, nch);
771  for (int i = 0; i < nch; i++) pt_pos[i] /= max;
772  int index = (int)(nch / width * pt);
773 
774  if (pt_pos[index] > rndm) return 1;
775  return 0;
776 }
777 
778 int SGCosmic::mkPt_neg_v2(const double pt, const float rndm)
779 {
780  const int nch = 100;
781  const int width = 50;
782  double pt_neg[nch] = { 3794.000, 14288.00, 15015.00, 13120.00, 10865.00,
783  9167.000, 7826.000, 6733.000, 5760.000, 5253.000,
784  4506.000, 3796.000, 3512.000, 3197.000, 2848.000,
785  2463.000, 2232.000, 2083.000, 1827.000, 1679.000,
786  1556.000, 1414.000, 1299.000, 1163.000, 1116.000,
787  1069.000, 964.0000, 922.0000, 818.0000, 764.0000,
788  723.0000, 658.0000, 614.0000, 541.0000, 552.0000,
789  520.0000, 491.0000, 382.0000, 380.0000, 392.0000,
790  346.0000, 352.0000, 321.0000, 298.0000, 272.0000,
791  277.0000, 246.0000, 238.0000, 237.0000, 241.0000,
792  189.0000, 208.0000, 196.0000, 194.0000, 169.0000,
793  149.0000, 143.0000, 139.0000, 130.0000, 134.0000,
794  125.0000, 116.0000, 110.0000, 104.0000, 122.0000,
795  110.0000, 99.00000, 95.00000, 91.00000, 87.00000,
796  91.00000, 75.00000, 87.00000, 76.00000, 81.00000,
797  81.00000, 67.00000, 64.00000, 55.00000, 59.00000,
798  52.00000, 53.00000, 53.00000, 51.00000, 46.00000,
799  40.00000, 49.00000, 41.00000, 37.00000, 44.00000,
800  44.00000, 42.00000, 47.00000, 30.00000, 30.00000,
801  42.00000, 28.00000, 28.00000, 32.00000, 26.00000
802  };
803  // normalize by peak value --> max = 1
804  double max = findMax(pt_neg, nch);
805  for (int i = 0; i < nch; i++) pt_neg[i] /= max;
806  int index = (int)(nch / width * pt);
807 
808  if (pt_neg[index] > rndm) return 1;
809  return 0;
810 }
811 
812 int SGCosmic::mkDz_pos_v2(const double dz, const float rndm)
813 {
814  const int nch = 200;
815  const int width = 100;
816  double dz_pos[nch] = { 0., 0., 0., 1.000000, 0.,
817  0., 0., 0., 0., 0.,
818  0., 0., 0., 0., 0.,
819  0., 0., 0., 0., 0.,
820  0., 0., 0., 0., 0.,
821  0., 0., 0., 0., 0.,
822  0., 0., 0., 0., 0.,
823  0., 0., 0., 1.000000, 0.,
824  1.000000, 0., 0., 4.000000, 1.000000,
825  1.000000, 2.000000, 1.000000, 3.000000, 6.000000,
826  20.00000, 11.00000, 30.00000, 45.00000, 64.00000,
827  139.0000, 243.0000, 391.0000, 460.0000, 586.0000,
828  720.0000, 790.0000, 931.0000, 967.0000, 1087.000,
829  1204.000, 1231.000, 1272.000, 1251.000, 1355.000,
830  1453.000, 1447.000, 1370.000, 1364.000, 1460.000,
831  1405.000, 1545.000, 1454.000, 1536.000, 1499.000,
832  1547.000, 1500.000, 1465.000, 1512.000, 1515.000,
833  1547.000, 1548.000, 1554.000, 1579.000, 1599.000,
834  1590.000, 1609.000, 1618.000, 1710.000, 1599.000,
835  1746.000, 1669.000, 1709.000, 1814.000, 1904.000,
836  1896.000, 1634.000, 1598.000, 1680.000, 1727.000,
837  1714.000, 1702.000, 1737.000, 1693.000, 1743.000,
838  1713.000, 1663.000, 1639.000, 1569.000, 1695.000,
839  1697.000, 1745.000, 1708.000, 1682.000, 1686.000,
840  1676.000, 1638.000, 1660.000, 1655.000, 1592.000,
841  1617.000, 1581.000, 1514.000, 1601.000, 1553.000,
842  1451.000, 1528.000, 1456.000, 1489.000, 1465.000,
843  1455.000, 1489.000, 1375.000, 1387.000, 1330.000,
844  1298.000, 1304.000, 1237.000, 1140.000, 1111.000,
845  1038.000, 942.0000, 959.0000, 867.0000, 822.0000,
846  773.0000, 698.0000, 649.0000, 591.0000, 556.0000,
847  513.0000, 474.0000, 424.0000, 408.0000, 414.0000,
848  324.0000, 296.0000, 265.0000, 259.0000, 188.0000,
849  170.0000, 195.0000, 163.0000, 140.0000, 131.0000,
850  122.0000, 113.0000, 98.00000, 92.00000, 68.00000,
851  61.00000, 51.00000, 44.00000, 46.00000, 30.00000,
852  32.00000, 25.00000, 20.00000, 13.00000, 11.00000,
853  3.000000, 2.000000, 1.000000, 1.000000, 0.,
854  0., 0., 0., 0., 0.,
855  0., 0., 0., 0., 0.
856  };
857  // normalize by peak value --> max = 1
858  double max = findMax(dz_pos, nch);
859  for (int i = 0; i < nch; i++) dz_pos[i] /= max;
860  int index = (int)((double)nch / (double)width * (dz + 50.));
861 
862  if (dz_pos[index] > rndm) return 1;
863  return 0;
864 }
865 
866 int SGCosmic::mkDz_neg_v2(const double dz, const float rndm)
867 {
868  const int nch = 100;
869  //const int width = 100;
870  double dz_neg[nch] = { 0., 1.000000, 0., 0., 0.,
871  0., 0., 0., 0., 0.,
872  0., 0., 0., 0., 0.,
873  0., 0., 0., 0., 1.000000,
874  2.000000, 1.000000, 2.000000, 5.000000, 12.00000,
875  25.00000, 66.00000, 244.0000, 609.0000, 1092.000,
876  1592.000, 1849.000, 2264.000, 2537.000, 2692.000,
877  2889.000, 2711.000, 2901.000, 3015.000, 2976.000,
878  3061.000, 3001.000, 3066.000, 3058.000, 3146.000,
879  3198.000, 3308.000, 3349.000, 3391.000, 3720.000,
880  3552.000, 3291.000, 3404.000, 3459.000, 3447.000,
881  3318.000, 3248.000, 3416.000, 3408.000, 3408.000,
882  3334.000, 3345.000, 3251.000, 3127.000, 2999.000,
883  3038.000, 2944.000, 2970.000, 2921.000, 2603.000,
884  2633.000, 2304.000, 2209.000, 1866.000, 1660.000,
885  1511.000, 1269.000, 1065.000, 893.0000, 765.0000,
886  628.0000, 456.0000, 416.0000, 378.0000, 273.0000,
887  221.0000, 166.0000, 134.0000, 116.0000, 62.00000,
888  46.00000, 35.00000, 11.00000, 5.000000, 0.,
889  0., 0., 0., 1.000000, 0.
890  };
891  // normalize by peak value --> max = 1
892  double max = findMax(dz_neg, nch);
893  for (int i = 0; i < nch; i++) dz_neg[i] /= max;
894  // nch is equal to width here, not needed to compute the ratio
895  int index = (int)(dz + 50.);
896  if (dz_neg[index] > rndm) return 1;
897  return 0;
898 }
899 
900 int SGCosmic::mkTanl_pos_v2(const double tanl, const float rndm)
901 {
902  const int nch = 100;
903  const int width = 5;
904  double tanl_pos[nch] = { 0., 0., 0., 0., 0.,
905  1.000000, 0., 2.000000, 0., 7.000000,
906  2.000000, 6.000000, 4.000000, 17.00000, 20.00000,
907  21.00000, 18.00000, 34.00000, 36.00000, 47.00000,
908  64.00000, 85.00000, 126.0000, 116.0000, 163.0000,
909  214.0000, 262.0000, 330.0000, 415.0000, 477.0000,
910  556.0000, 641.0000, 852.0000, 956.0000, 1243.000,
911  1370.000, 1699.000, 1995.000, 2415.000, 3040.000,
912  3352.000, 3889.000, 4390.000, 4823.000, 5399.000,
913  5862.000, 6201.000, 6506.000, 6748.000, 6788.000,
914  6767.000, 6575.000, 6479.000, 6112.000, 5940.000,
915  5459.000, 4986.000, 4504.000, 3879.000, 3460.000,
916  2981.000, 2502.000, 2121.000, 1803.000, 1477.000,
917  1248.000, 1041.000, 927.0000, 708.0000, 605.0000,
918  461.0000, 424.0000, 341.0000, 257.0000, 230.0000,
919  167.0000, 141.0000, 119.0000, 84.00000, 80.00000,
920  46.00000, 44.00000, 39.00000, 30.00000, 17.00000,
921  18.00000, 13.00000, 6.000000, 6.000000, 5.000000,
922  2.000000, 2.000000, 2.000000, 0., 0.,
923  0., 0., 0., 0., 1.000000
924  };
925  // normalize by peak value --> max = 1
926  double max = findMax(tanl_pos, nch);
927  for (int i = 0; i < nch; i++) tanl_pos[i] /= max;
928  int index = (int)((double)nch / (double)width * (tanl + 2.5));
929 
930  if (tanl_pos[index] > rndm) return 1;
931  return 0;
932 }
933 
934 int SGCosmic::mkTanl_neg_v2(const double tanl, const float rndm)
935 {
936  const int nch = 100;
937  const int width = 5;
938  double tanl_neg[nch] = { 0., 0., 0., 0., 0.,
939  0., 0., 1.000000, 4.000000, 0.,
940  3.000000, 4.000000, 9.000000, 8.000000, 14.00000,
941  21.00000, 25.00000, 41.00000, 44.00000, 49.00000,
942  73.00000, 84.00000, 107.0000, 139.0000, 166.0000,
943  242.0000, 268.0000, 316.0000, 412.0000, 469.0000,
944  595.0000, 727.0000, 900.0000, 1002.000, 1267.000,
945  1454.000, 1793.000, 2126.000, 2500.000, 2966.000,
946  3492.000, 3829.000, 4505.000, 4980.000, 5518.000,
947  5864.000, 6157.000, 6440.000, 6664.000, 6747.000,
948  6841.000, 6757.000, 6454.000, 6224.000, 5776.000,
949  5458.000, 4855.000, 4399.000, 3887.000, 3323.000,
950  3098.000, 2445.000, 1993.000, 1695.000, 1381.000,
951  1200.000, 993.0000, 844.0000, 640.0000, 590.0000,
952  484.0000, 412.0000, 352.0000, 255.0000, 207.0000,
953  160.0000, 121.0000, 119.0000, 103.0000, 57.00000,
954  53.00000, 37.00000, 27.00000, 26.00000, 25.00000,
955  23.00000, 8.000000, 7.000000, 9.000000, 3.000000,
956  5.000000, 4.000000, 2.000000, 3.000000, 3.000000,
957  0., 0., 0., 0., 0.
958  };
959  // normalize by peak value --> max = 1
960  double max = findMax(tanl_neg, nch);
961  for (int i = 0; i < nch; i++) tanl_neg[i] /= max;
962  int index = (int)((double)nch / (double)width * (tanl + 2.5));
963 
964  if (tanl_neg[index] > rndm) return 1;
965  return 0;
966 }
Belle2::SGCosmic::Parameters::ipdz
double ipdz
Vertex restriction in z direction.
Definition: SGCosmic.h:58
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::CDC::Helix
Helix parameter class.
Definition: Helix.h:51
Belle2::SGCosmic::Parameters::ptmin
double ptmin
Minimum value of the transverse momentum.
Definition: SGCosmic.h:62
Belle2::SGCosmic::Parameters::level
int level
Generator version: level 1 (default) or 2.
Definition: SGCosmic.h:46
Belle2::SGCosmic::Parameters
Struct to keep all necessary parameters for the cosmic generator.
Definition: SGCosmic.h:42
Belle2::SGCosmic::Parameters::ipdr
double ipdr
Vertex restriction in the radial direction.
Definition: SGCosmic.h:54
Belle2::SGCosmic::Parameters::cylindricalR
double cylindricalR
Cylindrical radius of generation.
Definition: SGCosmic.h:66
Belle2::SGCosmic::Parameters::ipRequirement
int ipRequirement
Restrict the vertex to IP or not (default)
Definition: SGCosmic.h:50
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86