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