Belle II Software  release-08-01-10
InvariantMassMuMuIntegrator.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 
10 
11 #include <tracking/calibration/InvariantMassMuMuIntegrator.h>
12 #include <tracking/calibration/InvariantMassMuMuStandAlone.h>
13 
14 #include <cassert>
15 #include <cmath>
16 
17 #include <Math/SpecFuncMathCore.h>
18 #include <Math/DistFunc.h>
19 
20 namespace Belle2::InvariantMassMuMuCalib {
21 
22  void InvariantMassMuMuIntegrator::init(double Mean,
23  double Sigma,
24  double SigmaK,
25  double BMean,
26  double BDelta,
27  double Tau,
28  double SigmaE,
29  double Frac,
30  double M0,
31  double Eps,
32  double CC,
33  double Slope,
34  double X)
35  {
36  m_mean = Mean;
37  m_sigma = Sigma;
38  m_sigmaK = SigmaK;
39  m_bMean = BMean;
40  m_bDelta = BDelta;
41  m_tauL = Tau;
42  m_tauR = Tau;
43  m_sigmaE = SigmaE;
44  m_frac = Frac;
45 
46  m_m0 = M0;
47  m_eps = Eps;
48  m_C = CC;
49  m_slope = Slope;
50 
51  m_x = X;
52  }
53 
54 
55 
56 
58  double InvariantMassMuMuIntegrator::eval(double t)
59  {
60  double CoreC = gausExpConv(m_mean, m_sigma, m_bMean, m_bDelta, m_tauL, m_tauR, m_sigmaK, m_x + t - m_m0);
61  double CoreE = 1. / (sqrt(2 * M_PI) * m_sigmaE) * exp(-1. / 2 * pow((m_x + t - m_m0 - m_mean) / m_sigmaE, 2));
62  double Core = (1 - m_frac) * CoreC + m_frac * CoreE;
63 
64  assert(t >= 0);
65 
66  double K = 0;
67  if (t >= m_eps)
68  K = pow(t, -m_slope);
69  else if (t >= 0)
70  K = (1 + (m_C - 1) * (m_eps - t) / m_eps) * pow(m_eps, -m_slope);
71  else
72  K = 0;
73 
74  return Core * K;
75  }
76 
77  // Simple integration of the PDF based on the Trapezoidal rule
78  double InvariantMassMuMuIntegrator::integralTrap(double a, double b)
79  {
80  const int N = 14500000;
81  double step = (b - a) / N;
82  double s = 0;
83  for (int i = 0; i <= N; ++i) {
84  double K = (i == 0 || i == N) ? 0.5 : 1;
85  double t = a + i * step;
86  s += eval(t) * step * K;
87  }
88 
89 
90  return s;
91  }
92 
93  // Integration of the PDF which avoids steps and uses variable transformation (Trapezoidal rule as back-end)
94  double InvariantMassMuMuIntegrator::integralTrapImp(double Eps, double a)
95  {
96  double tMin = m_slope * m_tauL;
97 
98  //only one function type
99  if (m_x - m_m0 >= 0) {
100 
101  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
102  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, Eps / m_tauL);
103 
104 
105  const int N = 15000;
106  double step = (r2 - r1) / N;
107  double s = 0;
108  for (int i = 0; i <= N; ++i) {
109  double r = r1 + step * i;
110  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
111  double K = (i == 0 || i == N) ? 0.5 : 1;
112 
113  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
114 
115  s += eval(t) / est * step * K;
116  }
117 
118  s *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
119 
120  return s;
121 
122  }
123 
124  else if (m_x - m_m0 >= -tMin) {
125 
126  //integrate from m_m0 - m_x to a
127  double s1 = 0;
128 
129  {
130  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
131  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, (m_m0 - m_x) / m_tauL);
132 
133 
134  const int N = 150000;
135  double step = (r2 - r1) / N;
136  for (int i = 0; i <= N; ++i) {
137  double r = r1 + step * i;
138  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
139  double K = (i == 0 || i == N) ? 0.5 : 1;
140 
141  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
142 
143  s1 += eval(t) / est * step * K;
144  }
145 
146  s1 *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
147  }
148 
149 
150  // integrate from eps to (m_m0 - m_x)
151  double s2 = 0;
152 
153  {
154  double r1 = pow(Eps, -m_slope + 1) / (1 - m_slope);
155  double r2 = pow(m_m0 - m_x, -m_slope + 1) / (1 - m_slope);
156 
157 
158  const int N = 150000;
159  double step = (r2 - r1) / N;
160  for (int i = 0; i <= N; ++i) {
161  double r = r1 + step * i;
162  double t = pow(r * (1 - m_slope), 1. / (1 - m_slope));
163  double K = (i == 0 || i == N) ? 0.5 : 1;
164 
165  double est = pow(t, -m_slope);
166 
167  s2 += eval(t) / est * step * K;
168  }
169 
170  }
171 
172  return (s1 + s2);
173 
174  } else {
175 
176  //integrate from m_m0 - m_x to a
177  double s1 = 0;
178 
179  {
180  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
181  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, (m_m0 - m_x) / m_tauL);
182 
183 
184  const int N = 150000;
185  double step = (r2 - r1) / N;
186  for (int i = 0; i <= N; ++i) {
187  double r = r1 + step * i;
188  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
189  double K = (i == 0 || i == N) ? 0.5 : 1;
190 
191  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
192 
193  s1 += eval(t) / est * step * K;
194  }
195 
196  s1 *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
197  }
198 
199 
200  // integrate from eps to tMin
201  double s2 = 0;
202 
203  {
204  double r1 = pow(Eps, -m_slope + 1) / (1 - m_slope);
205  double r2 = pow(tMin, -m_slope + 1) / (1 - m_slope);
206 
207 
208  const int N = 150000;
209  double step = (r2 - r1) / N;
210  for (int i = 0; i <= N; ++i) {
211  double r = r1 + step * i;
212  double t = pow(r * (1 - m_slope), 1. / (1 - m_slope));
213  double K = (i == 0 || i == N) ? 0.5 : 1;
214 
215  double est = pow(t, -m_slope);
216 
217  s2 += eval(t) / est * step * K;
218  }
219 
220  }
221 
222  //integrate from tMin to m_m0 - m_x
223  double s3 = 0;
224 
225  {
226  double r1 = exp(tMin / m_tauR) * m_tauR;
227  double r2 = exp((m_m0 - m_x) / m_tauR) * m_tauR;
228 
229 
230  const int N = 150000;
231  double step = (r2 - r1) / N;
232  for (int i = 0; i <= N; ++i) {
233  double r = r1 + step * i;
234  double t = log(r / m_tauR) * m_tauR;
235  double K = (i == 0 || i == N) ? 0.5 : 1;
236 
237  double est = exp(t / m_tauR);
238 
239  s3 += eval(t) / est * step * K;
240  }
241 
242 
243  }
244 
245  return (s1 + s2 + s3);
246 
247  }
248 
249  return 0;
250 
251  }
252 
253 
254  // Integration of the PDF which avoids steps and uses variable transformation (Gauss-Konrod rule as back-end)
255  double InvariantMassMuMuIntegrator::integralKronrod(double a)
256  {
257 
258  double s0 = integrate([&](double t) {
259  return eval(t);
260  }, 0, m_eps);
261 
262 
263  double tMin = m_slope * m_tauL;
264 
265  //only one function type
266  if (m_x - m_m0 >= 0) {
267 
268  assert(m_eps <= a);
269 
270  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
271  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, m_eps / m_tauL);
272 
273  double s = integrate([&](double r) {
274 
275  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
276  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
277  return eval(t) / est;
278 
279 
280  }, r1, r2);
281 
282  s *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
283 
284  return (s0 + s);
285 
286  }
287 
288  //only one function type
289  else if (m_x - m_m0 >= -m_eps) {
290 
291  assert(0 <= m_m0 - m_x);
292  assert(m_m0 - m_x <= m_eps);
293  assert(m_eps <= a);
294 
295 
296  double s01 = integrate([&](double t) {
297  return eval(t);
298  }, 0, m_m0 - m_x);
299 
300 
301  double s02 = integrate([&](double t) {
302  return eval(t);
303  }, m_m0 - m_x, m_eps);
304 
305 
306  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
307  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, m_eps / m_tauL);
308 
309  double s = integrate([&](double r) {
310 
311  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
312  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
313  return eval(t) / est;
314 
315 
316  }, r1, r2);
317 
318  s *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
319 
320  return (s01 + s02 + s);
321 
322  }
323 
324 
325  //two function types
326  else if (m_x - m_m0 >= -tMin) {
327 
328  assert(m_eps <= m_m0 - m_x);
329  assert(m_m0 - m_x <= a);
330 
331  //integrate from m_m0 - m_x to a
332  double s1 = 0;
333 
334  {
335  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
336  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, (m_m0 - m_x) / m_tauL);
337 
338 
339  s1 = integrate([&](double r) {
340 
341  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
342  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
343  return eval(t) / est;
344 
345 
346  }, r1, r2);
347 
348  s1 *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
349 
350  }
351 
352 
353  // integrate from eps to (m_m0 - m_x)
354  double s2 = 0;
355 
356  {
357  double r1 = pow(m_eps, -m_slope + 1) / (1 - m_slope);
358  double r2 = pow(m_m0 - m_x, -m_slope + 1) / (1 - m_slope);
359 
360 
361 
362  s2 = integrate([&](double r) {
363  double t = pow(r * (1 - m_slope), 1. / (1 - m_slope));
364  double est = pow(t, -m_slope);
365  return eval(t) / est;
366  }, r1, r2);
367  }
368 
369  return (s0 + s1 + s2);
370 
371  } else {
372 
373  assert(m_eps <= tMin);
374  assert(tMin <= m_m0 - m_x);
375  assert(m_m0 - m_x <= a);
376 
377  //integrate from m_m0 - m_x to a
378  double s1 = 0;
379 
380  {
381  double r1 = ROOT::Math::inc_gamma_c(1 - m_slope, a / m_tauL);
382  double r2 = ROOT::Math::inc_gamma_c(1 - m_slope, (m_m0 - m_x) / m_tauL);
383 
384 
385  s1 = integrate([&](double r) {
386 
387  double t = m_tauL * ROOT::Math::gamma_quantile_c(r, 1 - m_slope, 1);
388  double est = pow(t, -m_slope) * exp(- (m_x - m_m0 + t) / m_tauL);
389  return eval(t) / est;
390 
391 
392  }, r1, r2);
393 
394  s1 *= exp((m_m0 - m_x) / m_tauL) * pow(m_tauL, -m_slope + 1) * ROOT::Math::tgamma(1 - m_slope);
395 
396  }
397 
398 
399  // integrate from eps to tMin
400  double s2 = 0;
401 
402  {
403  double r1 = pow(m_eps, -m_slope + 1) / (1 - m_slope);
404  double r2 = pow(tMin, -m_slope + 1) / (1 - m_slope);
405 
406  s2 = integrate([&](double r) {
407  double t = pow(r * (1 - m_slope), 1. / (1 - m_slope));
408  double est = pow(t, -m_slope);
409  return eval(t) / est;
410  }, r1, r2);
411 
412  }
413 
414  //integrate from tMin to m_m0 - m_x
415  double s3 = 0;
416 
417  {
418  double r1 = exp(tMin / m_tauR) * m_tauR;
419  double r2 = exp((m_m0 - m_x) / m_tauR) * m_tauR;
420 
421 
422  s3 = integrate([&](double r) {
423  double t = log(r / m_tauR) * m_tauR;
424  double est = exp(t / m_tauR);
425  return eval(t) / est;
426  }, r1, r2);
427 
428  }
429 
430  return (s0 + s1 + s2 + s3);
431 
432 
433  }
434 
435  return 0;
436 
437  }
438 
439 }
K
#define K(x)
macro autogenerated by FFTW
Definition: DiscreteCosineTransform_31points.cc:37
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::integralKronrod
double integralKronrod(double a)
Integration of the PDF which avoids steps and uses variable transformation (Gauss-Konrod rule as back...
Definition: InvariantMassMuMuIntegrator.cc:255
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_eps
double m_eps
cut-off term for the power-spectrum caused by the ISR (in GeV)
Definition: InvariantMassMuMuIntegrator.h:63
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_C
double m_C
the coeficient between part bellow eps and above eps cut-off
Definition: InvariantMassMuMuIntegrator.h:68
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_x
double m_x
the resulting PDF is function of this variable the actual rec-level mass
Definition: InvariantMassMuMuIntegrator.h:66
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::init
void init(double Mean, double Sigma, double SigmaK, double BMean, double BDelta, double Tau, double SigmaE, double Frac, double M0, double Eps, double CC, double Slope, double X)
Init the parameters of the PDF integrator.
Definition: InvariantMassMuMuIntegrator.cc:22
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_bMean
double m_bMean
(bRight + bLeft)/2 where bLeft, bRight are the transition points between gaus and exp (in sigma)
Definition: InvariantMassMuMuIntegrator.h:54
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_slope
double m_slope
power in the power-like spectrum from the ISR
Definition: InvariantMassMuMuIntegrator.h:64
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::eval
double eval(double t)
evaluate the PDF integrand for given t - the integration variable
Definition: InvariantMassMuMuIntegrator.cc:58
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::integralTrap
double integralTrap(double a, double b)
Simple integration of the PDF for a to b based on the Trapezoidal rule (for validation)
Definition: InvariantMassMuMuIntegrator.cc:78
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_mean
double m_mean
mean position of the resolution function, i.e. (Gaus+Exp tails) conv Gaus
Definition: InvariantMassMuMuIntegrator.h:51
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::m_bDelta
double m_bDelta
(bRight - bLeft)/2 where bLeft, bRight are the transition points between gaus and exp (in sigma)
Definition: InvariantMassMuMuIntegrator.h:55
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::integralTrapImp
double integralTrapImp(double Eps, double a)
Integration of the PDF which avoids steps and uses variable transformation (Trapezoidal rule as back-...
Definition: InvariantMassMuMuIntegrator.cc:94
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28