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 58 of file InvariantMassMuMuIntegrator.cc.

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

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 }

◆ 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 255 of file InvariantMassMuMuIntegrator.cc.

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 }
Belle2::InvariantMassMuMuCalib::InvariantMassMuMuIntegrator::eval
double eval(double t)
evaluate the PDF integrand for given t - the integration variable
Definition: InvariantMassMuMuIntegrator.cc:58

◆ 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 78 of file InvariantMassMuMuIntegrator.cc.

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 }

◆ 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 94 of file InvariantMassMuMuIntegrator.cc.

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 }

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: