Belle II Software development
Public Member Functions | Private Attributes | List of all members
InvariantMassMuMuIntegrator Class Reference

The integrator aims to evaluate convolution of PDFgenLevel and resolution function. More...

#include <InvariantMassMuMuIntegrator.h>

Public Member Functions

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.
 
double eval (double t)
 evaluate the PDF integrand for given t - the integration variable
 
double integralTrap (double a, double b)
 Simple integration of the PDF for a to b based on the Trapezoidal rule (for validation)
 
double integralTrapImp (double Eps, double a)
 Integration of the PDF which avoids steps and uses variable transformation (Trapezoidal rule as back-end)
 
double integralKronrod (double a)
 Integration of the PDF which avoids steps and uses variable transformation (Gauss-Konrod rule as back-end)
 

Private Attributes

double m_mean = 4
 mean position of the resolution function, i.e. (Gaus+Exp tails) conv Gaus
 
double m_sigma = 30
 sigma of the resolution function
 
double m_sigmaK = 30
 sigma of the gaus in the convolution
 
double m_bMean = 0
 (bRight + bLeft)/2 where bLeft, bRight are the transition points between gaus and exp (in sigma)
 
double m_bDelta = 2.6
 (bRight - bLeft)/2 where bLeft, bRight are the transition points between gaus and exp (in sigma)
 
double m_tauL = 60
 1/slope of the left exponential tail
 
double m_tauR = 60
 1/slope of the right exponential tail
 
double m_sigmaE = 30
 sigma of the external gaus
 
double m_frac = 0.1
 fraction of events in the external gaus
 
double m_m0 = 10500
 invariant mass of the collisions
 
double m_eps = 0.01
 cut-off term for the power-spectrum caused by the ISR (in GeV)
 
double m_slope = 0.95
 power in the power-like spectrum from the ISR
 
double m_x = 10400
 the resulting PDF is function of this variable the actual rec-level mass
 
double m_C = 16
 the coefficient between part below eps and above eps cut-off
 

Detailed Description

The integrator aims to evaluate convolution of PDFgenLevel and resolution function.

Definition at line 18 of file InvariantMassMuMuIntegrator.h.

Member Function Documentation

◆ eval()

double eval ( double t)

evaluate the PDF integrand for given t - the integration variable

Definition at line 59 of file InvariantMassMuMuIntegrator.cc.

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 }
K
#define K(x)
macro autogenerated by FFTW
Definition DiscreteCosineTransform_31points.cc:37
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

◆ 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 at line 23 of file InvariantMassMuMuIntegrator.cc.

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 }

◆ integralKronrod()

double integralKronrod ( double a)

Integration of the PDF which avoids steps and uses variable transformation (Gauss-Konrod rule as back-end)

Definition at line 256 of file InvariantMassMuMuIntegrator.cc.

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 }
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

◆ integralTrap()

double integralTrap ( double a,
double b )

Simple integration of the PDF for a to b based on the Trapezoidal rule (for validation)

Definition at line 79 of file InvariantMassMuMuIntegrator.cc.

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 }

◆ integralTrapImp()

double integralTrapImp ( double Eps,
double a )

Integration of the PDF which avoids steps and uses variable transformation (Trapezoidal rule as back-end)

Definition at line 95 of file InvariantMassMuMuIntegrator.cc.

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 }

Member Data Documentation

◆ m_bDelta

double m_bDelta = 2.6
private

(bRight - bLeft)/2 where bLeft, bRight are the transition points between gaus and exp (in sigma)

Definition at line 55 of file InvariantMassMuMuIntegrator.h.

◆ m_bMean

double m_bMean = 0
private

(bRight + bLeft)/2 where bLeft, bRight are the transition points between gaus and exp (in sigma)

Definition at line 54 of file InvariantMassMuMuIntegrator.h.

◆ m_C

double m_C = 16
private

the coefficient between part below eps and above eps cut-off

Definition at line 68 of file InvariantMassMuMuIntegrator.h.

◆ m_eps

double m_eps = 0.01
private

cut-off term for the power-spectrum caused by the ISR (in GeV)

Definition at line 63 of file InvariantMassMuMuIntegrator.h.

◆ m_frac

double m_frac = 0.1
private

fraction of events in the external gaus

Definition at line 60 of file InvariantMassMuMuIntegrator.h.

◆ m_m0

double m_m0 = 10500
private

invariant mass of the collisions

Definition at line 62 of file InvariantMassMuMuIntegrator.h.

◆ m_mean

double m_mean = 4
private

mean position of the resolution function, i.e. (Gaus+Exp tails) conv Gaus

Definition at line 51 of file InvariantMassMuMuIntegrator.h.

◆ m_sigma

double m_sigma = 30
private

sigma of the resolution function

Definition at line 52 of file InvariantMassMuMuIntegrator.h.

◆ m_sigmaE

double m_sigmaE = 30
private

sigma of the external gaus

Definition at line 59 of file InvariantMassMuMuIntegrator.h.

◆ m_sigmaK

double m_sigmaK = 30
private

sigma of the gaus in the convolution

Definition at line 53 of file InvariantMassMuMuIntegrator.h.

◆ m_slope

double m_slope = 0.95
private

power in the power-like spectrum from the ISR

Definition at line 64 of file InvariantMassMuMuIntegrator.h.

◆ m_tauL

double m_tauL = 60
private

1/slope of the left exponential tail

Definition at line 56 of file InvariantMassMuMuIntegrator.h.

◆ m_tauR

double m_tauR = 60
private

1/slope of the right exponential tail

Definition at line 57 of file InvariantMassMuMuIntegrator.h.

◆ m_x

double m_x = 10400
private

the resulting PDF is function of this variable the actual rec-level mass

Definition at line 66 of file InvariantMassMuMuIntegrator.h.

The documentation for this class was generated from the following files: