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