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
14#include <cassert>
15#include <cmath>
16
17#include <Math/SpecFuncMathCore.h>
18#include <Math/DistFunc.h>
19
20namespace 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 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: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