Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooHypatia2.cxx
Go to the documentation of this file.
1// Author: Stephan Hageboeck, CERN, Oct 2019
2// Based on RooIpatia2 by Diego Martinez Santos, Nikhef, Diego.Martinez.Santos@cern.ch
3/*****************************************************************************
4 * Project: RooFit *
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2019, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18 * \class RooHypatia2
19 *
20 * RooHypatia2 is the two-sided version of the Hypatia distribution described in https://arxiv.org/abs/1312.5000.
21 * \image html RooHypatia2.png
22 *
23 * It has a hyperbolic core of a crystal-ball-like function \f$ G \f$ and two tails:
24 * \f[
25 * \mathrm{Hypatia2}(x, \mu, \sigma, \lambda, \zeta, \beta, a_l, n_l, a_r, n_r) =
26 * \begin{cases}
27 * \frac{ G(\mu - a_l \sigma, \mu, \sigma, \lambda, \zeta, \beta) }
28 * { \left( 1 - \frac{x}{n_l G(\ldots)/G'(\ldots) - a_l\sigma } \right)^{n_l} }
29 * & \text{if } \frac{x-\mu}{\sigma} < -a_l \\
30 * \left( (x-\mu)^2 + A^2_\lambda(\zeta)\sigma^2 \right)^{\frac{1}{2}\lambda-\frac{1}{4}} e^{\beta(x-\mu)} K_{\lambda-\frac{1}{2}}
31 * \left( \zeta \sqrt{1+\left( \frac{x-\mu}{A_\lambda(\zeta)\sigma} \right)^2 } \right) \equiv G(x, \mu, \ldots)
32 * & \text{otherwise} \\
33 * \frac{ G(\mu + a_r \sigma, \mu, \sigma, \lambda, \zeta, \beta) }
34 * { \left( 1 - \frac{x}{-n_r G(\ldots)/G'(\ldots) - a_r\sigma } \right)^{n_r} }
35 * & \text{if } \frac{x-\mu}{\sigma} > a_r \\
36 * \end{cases}
37 * \f]
38 * Here, \f$ K_\lambda \f$ are the modified Bessel functions of the second kind
39 * ("irregular modified cylindrical Bessel functions" from the gsl,
40 * "special Bessel functions of the third kind"),
41 * and \f$ A^2_\lambda(\zeta) \f$ is a ratio of these:
42 * \f[
43 * A_\lambda^{2}(\zeta) = \frac{\zeta K_\lambda(\zeta)}{K_{\lambda+1}(\zeta)}
44 * \f]
45 *
46 * \if false
47 * TODO Enable once analytic integrals work.
48 * ### Analytical Integration
49 * The Hypatia distribution can be integrated analytically if \f$ \beta = \zeta = 0 \f$ and
50 * \f$ \lambda < 0 \f$. An analytic integral will only be used, though, if the parameters are **constant**
51 * at zero, and if \f$ \lambda < 0 \f$. This can be ensured as follows:
52 * ```
53 * RooRealVar beta("beta", "beta", 0.); // NOT beta("beta", "beta", 0., -1., 1.) This would allow it to float.
54 * RooRealVar zeta("zeta", "zeta", 0.);
55 * RooRealVar lambda("lambda", "lambda", -1., -10., -0.00001);
56 * ```
57 * In all other cases, slower / less accurate numeric integration will be used.
58 * Note that including `0.` in the value range of lambda excludes using analytic integrals.
59 * \endif
60 *
61 * ### Concavity
62 * Note that unless the parameters \f$ a_l,\ a_r \f$ are very large, the function has non-hyperbolic tails. This requires
63 * \f$ G \f$ to be strictly concave, *i.e.*, peaked, as otherwise the tails would yield imaginary numbers. Choosing \f$ \lambda,
64 * \beta, \zeta \f$ inappropriately will therefore lead to evaluation errors.
65 *
66 * Further, the original paper establishes that to keep the tails from rising,
67 * \f[
68 * \begin{split}
69 * \beta^2 &< \alpha^2 \\
70 * \Leftrightarrow \beta^2 &< \frac{\zeta^2}{\delta^2} = \frac{\zeta^2}{\sigma^2 A_{\lambda}^2(\zeta)}
71 * \end{split}
72 * \f]
73 * needs to be satisfied, unless the fit range is very restricted, because otherwise, the function rises in the tails.
74 *
75 *
76 * In case of evaluation errors, it is advisable to choose very large values for \f$ a_l,\ a_r \f$, tweak the parameters of the core region to
77 * make it concave, and re-enable the tails. Especially \f$ \beta \f$ needs to be close to zero.
78 *
79 * ## Relation to RooIpatia2
80 * This implementation is largely based on RooIpatia2, https://gitlab.cern.ch/lhcb/Urania/-/blob/master/PhysFit/B2DXFitters/src/RooIpatia2.cxx,
81 * but there are differences:
82 * - At construction time, the Hypatia implementation checks if the range of parameters extends into regions where
83 * the function might be undefined.
84 * - Hypatia supports I/O to ROOT files.
85 * - Hypatia will support faster batched function evaluations.
86 * - Hypatia might support analytical integration for the case \f$ \zeta = \beta = 0, \lambda < 1 \f$.
87 *
88 * Because of these differences, and to avoid name clashes for setups where RooFit is used in an environment that also has
89 * RooIpatia2, class names need to be different.
90 */
91
92#include "RooHypatia2.h"
93#include "BracketAdapters.h"
94#include "RooAbsReal.h"
95#include "RooHelpers.h"
96
97#include "TMath.h"
98#include "Math/SpecFunc.h"
99#include "ROOT/RConfig.hxx"
100
101#include <cmath>
102#include <algorithm>
103
104///////////////////////////////////////////////////////////////////////////////////////////
105/// Construct a new Hypatia2 PDF.
106/// \param[in] name Name of this instance.
107/// \param[in] title Title (for plotting)
108/// \param[in] x The variable of this distribution
109/// \param[in] lambda Shape parameter. Note that \f$ \lambda < 0 \f$ is required if \f$ \zeta = 0 \f$.
110/// \param[in] zeta Shape parameter (\f$ \zeta >= 0 \f$).
111/// \param[in] beta Asymmetry parameter \f$ \beta \f$. Symmetric case is \f$ \beta = 0 \f$,
112/// choose values close to zero.
113/// \param[in] argSigma Width parameter. If \f$ \beta = 0, \ \sigma \f$ is the RMS width.
114/// \param[in] mu Location parameter. Shifts the distribution left/right.
115/// \param[in] a Start of the left tail (\f$ a \geq 0 \f$, to the left of the peak). Note that when setting \f$ a = \sigma = 1 \f$,
116/// the tail region is to the left of \f$ x = \mu - 1 \f$, so a should be positive.
117/// \param[in] n Shape parameter of left tail (\f$ n \ge 0 \f$). With \f$ n = 0 \f$, the function is constant.
118/// \param[in] a2 Start of right tail.
119/// \param[in] n2 Shape parameter of right tail (\f$ n2 \ge 0 \f$). With \f$ n2 = 0 \f$, the function is constant.
120RooHypatia2::RooHypatia2(const char *name, const char *title, RooAbsReal& x, RooAbsReal& lambda,
121 RooAbsReal& zeta, RooAbsReal& beta, RooAbsReal& argSigma, RooAbsReal& mu, RooAbsReal& a,
122 RooAbsReal& n, RooAbsReal& a2, RooAbsReal& n2) :
123 RooAbsPdf(name, title),
124 _x("x", "x", this, x),
125 _lambda("lambda", "Lambda", this, lambda),
126 _zeta("zeta", "zeta", this, zeta),
127 _beta("beta", "Asymmetry parameter beta", this, beta),
128 _sigma("sigma", "Width parameter sigma", this, argSigma),
129 _mu("mu", "Location parameter mu", this, mu),
130 _a("a", "Left tail location a", this, a),
131 _n("n", "Left tail parameter n", this, n),
132 _a2("a2", "Right tail location a2", this, a2),
133 _n2("n2", "Right tail parameter n2", this, n2)
134{
135 RooHelpers::checkRangeOfParameters(this, {&argSigma}, 0.);
136 RooHelpers::checkRangeOfParameters(this, {&zeta, &n, &n2, &a, &a2}, 0., std::numeric_limits<double>::max(), true);
137 if (zeta.getVal() == 0. && zeta.isConstant()) {
138 RooHelpers::checkRangeOfParameters(this, {&lambda}, -std::numeric_limits<double>::max(), 0., false,
139 std::string("Lambda needs to be negative when ") + _zeta.GetName() + " is zero.");
140 }
141}
142
143
144///////////////////////////////////////////////////////////////////////////////////////////
145/// Copy a new Hypatia2 PDF.
146/// \param[in] other Original to copy from.
147/// \param[in] name Optional new name.
148RooHypatia2::RooHypatia2(const RooHypatia2& other, const char* name) :
149 RooAbsPdf(other, name),
150 _x("x", this, other._x),
151 _lambda("lambda", this, other._lambda),
152 _zeta("zeta", this, other._zeta),
153 _beta("beta", this, other._beta),
154 _sigma("sigma", this, other._sigma),
155 _mu("mu", this, other._mu),
156 _a("a", this, other._a),
157 _n("n", this, other._n),
158 _a2("a2", this, other._a2),
159 _n2("n2", this, other._n2)
160{
161}
162
163namespace {
164const double sq2pi_inv = 1./std::sqrt(TMath::TwoPi());
165const double logsq2pi = std::log(std::sqrt(TMath::TwoPi()));
166const double ln2 = std::log(2.);
167
168double low_x_BK(double nu, double x){
169 return TMath::Gamma(nu)*std::pow(2., nu-1.)*std::pow(x, -nu);
170}
171
172double low_x_LnBK(double nu, double x){
173 return std::log(TMath::Gamma(nu)) + ln2*(nu-1.) - std::log(x) * nu;
174}
175
176double besselK(double ni, double x) {
177 const double nu = std::abs(ni);
178 if ((x < 1.e-06 && nu > 0.) ||
179 (x < 1.e-04 && nu > 0. && nu < 55.) ||
180 (x < 0.1 && nu >= 55.) )
181 return low_x_BK(nu, x);
182
183 return ROOT::Math::cyl_bessel_k(nu, x);
184}
185
186double LnBesselK(double ni, double x) {
187 const double nu = std::abs(ni);
188 if ((x < 1.e-06 && nu > 0.) ||
189 (x < 1.e-04 && nu > 0. && nu < 55.) ||
190 (x < 0.1 && nu >= 55.) )
191 return low_x_LnBK(nu, x);
192
193 return std::log(ROOT::Math::cyl_bessel_k(nu, x));
194}
195
196
197double LogEval(double d, double l, double alpha, double beta, double delta) {
198 const double gamma = alpha;//std::sqrt(alpha*alpha-beta*beta);
199 const double dg = delta*gamma;
200 const double thing = delta*delta + d*d;
201 const double logno = l*std::log(gamma/delta) - logsq2pi - LnBesselK(l, dg);
202
203 return std::exp(logno + beta*d
204 + (0.5-l)*(std::log(alpha)-0.5*std::log(thing))
205 + LnBesselK(l-0.5, alpha*std::sqrt(thing)) );// + std::log(std::abs(beta)+0.0001) );
206
207}
208
209
210double diff_eval(double d, double l, double alpha, double beta, double delta){
211 const double gamma = alpha;
212 const double dg = delta*gamma;
213
214 const double thing = delta*delta + d*d;
215 const double sqrthing = std::sqrt(thing);
216 const double alphasq = alpha*sqrthing;
217 const double no = std::pow(gamma/delta, l)/besselK(l, dg)*sq2pi_inv;
218 const double ns1 = 0.5-l;
219
220 return no * std::pow(alpha, ns1) * std::pow(thing, l/2. - 1.25)
221 * ( -d * alphasq * (besselK(l - 1.5, alphasq)
222 + besselK(l + 0.5, alphasq))
223 + (2.*(beta*thing + d*l) - d) * besselK(ns1, alphasq) )
224 * std::exp(beta*d) * 0.5;
225}
226
227/*
228double Gauss2F1(double a, double b, double c, double x){
229 if (std::abs(x) <= 1.) {
230 return ROOT::Math::hyperg(a, b, c, x);
231 } else {
232 return ROOT::Math::hyperg(c-a, b, c, 1-1/(1-x))/std::pow(1-x, b);
233 }
234}
235
236double stIntegral(double d1, double delta, double l){
237 return d1 * Gauss2F1(0.5, 0.5-l, 3./2, -d1*d1/(delta*delta));
238}
239*/
240}
241
243{
244 const double d = _x-_mu;
245 const double cons0 = std::sqrt(_zeta);
246 const double asigma = _a*_sigma;
247 const double a2sigma = _a2*_sigma;
248 const double beta = _beta;
249 double out = 0.;
250
251 if (_zeta > 0.) {
252 // careful if zeta -> 0. You can implement a function for the ratio,
253 // but careful again that |nu + 1 | != |nu| + 1 so you have to deal with the signs
254 const double phi = besselK(_lambda+1., _zeta) / besselK(_lambda, _zeta);
255 const double cons1 = _sigma/std::sqrt(phi);
256 const double alpha = cons0/cons1;
257 const double delta = cons0*cons1;
258
259 if (d < -asigma){
260 const double k1 = LogEval(-asigma, _lambda, alpha, beta, delta);
261 const double k2 = diff_eval(-asigma, _lambda, alpha, beta, delta);
262 const double B = -asigma + _n*k1/k2;
263 const double A = k1 * std::pow(B+asigma, _n);
264
265 out = A * std::pow(B-d, -_n);
266 }
267 else if (d>a2sigma) {
268 const double k1 = LogEval(a2sigma, _lambda, alpha, beta, delta);
269 const double k2 = diff_eval(a2sigma, _lambda, alpha, beta, delta);
270 const double B = -a2sigma - _n2*k1/k2;
271 const double A = k1 * std::pow(B+a2sigma, _n2);
272
273 out = A * std::pow(B+d, -_n2);
274 }
275 else {
276 out = LogEval(d, _lambda, alpha, beta, delta);
277 }
278 }
279 else if (_zeta < 0.) {
280 coutE(Eval) << "The parameter " << _zeta.GetName() << " of the RooHypatia2 " << GetName() << " cannot be < 0." << std::endl;
281 }
282 else if (_lambda < 0.) {
283 const double delta = _sigma;
284
285 if (d < -asigma ) {
286 const double cons1 = std::exp(-beta*asigma);
287 const double phi = 1. + _a * _a;
288 const double k1 = cons1 * std::pow(phi, _lambda-0.5);
289 const double k2 = beta*k1 - cons1*(_lambda-0.5) * std::pow(phi, _lambda-1.5) * 2.*_a/delta;
290 const double B = -asigma + _n*k1/k2;
291 const double A = k1*std::pow(B+asigma, _n);
292
293 out = A*std::pow(B-d, -_n);
294 }
295 else if (d > a2sigma) {
296 const double cons1 = std::exp(beta*a2sigma);
297 const double phi = 1. + _a2*_a2;
298 const double k1 = cons1 * std::pow(phi, _lambda-0.5);
299 const double k2 = beta*k1 + cons1*(_lambda-0.5) * std::pow(phi, _lambda-1.5) * 2.*_a2/delta;
300 const double B = -a2sigma - _n2*k1/k2;
301 const double A = k1*std::pow(B+a2sigma, _n2);
302
303 out = A*std::pow(B+d, -_n2);
304 }
305 else {
306 out = std::exp(beta*d) * std::pow(1. + d*d/(delta*delta), _lambda - 0.5);
307 }
308 }
309 else {
310 coutE(Eval) << "zeta = 0 only supported for lambda < 0. lambda = " << double(_lambda) << std::endl;
311 }
312
313 return out;
314}
315
316namespace {
317//////////////////////////////////////////////////////////////////////////////////////////
318/// The generic compute function that recalculates everything for every loop iteration.
319/// This is only needed in the rare case where a parameter is used as an observable.
320template<typename Tx, typename Tl, typename Tz, typename Tb, typename Ts, typename Tm, typename Ta, typename Tn,
321typename Ta2, typename Tn2>
322void compute(std::span<double> output, Tx x, Tl lambda, Tz zeta, Tb beta, Ts sigma, Tm mu, Ta a, Tn n, Ta2 a2, Tn2 n2) {
323 const auto N = output.size();
324 const bool zetaIsAlwaysLargerZero = !zeta.isBatch() && zeta[0] > 0.;
325 const bool zetaIsAlwaysZero = !zeta.isBatch() && zeta[0] == 0.;
326
327 for (std::size_t i = 0; i < N; ++i) {
328
329 const double d = x[i] - mu[i];
330 const double cons0 = std::sqrt(zeta[i]);
331 const double asigma = a[i]*sigma[i];
332 const double a2sigma = a2[i]*sigma[i];
333// const double beta = beta[i];
334
335 if (zetaIsAlwaysLargerZero || zeta[i] > 0.) {
336 const double phi = besselK(lambda[i]+1., zeta[i]) / besselK(lambda[i], zeta[i]);
337 const double cons1 = sigma[i]/std::sqrt(phi);
338 const double alpha = cons0/cons1;
339 const double delta = cons0*cons1;
340
341 if (d < -asigma){
342 const double k1 = LogEval(-asigma, lambda[i], alpha, beta[i], delta);
343 const double k2 = diff_eval(-asigma, lambda[i], alpha, beta[i], delta);
344 const double B = -asigma + n[i]*k1/k2;
345 const double A = k1 * std::pow(B+asigma, n[i]);
346
347 output[i] = A * std::pow(B-d, -n[i]);
348 }
349 else if (d>a2sigma) {
350 const double k1 = LogEval(a2sigma, lambda[i], alpha, beta[i], delta);
351 const double k2 = diff_eval(a2sigma, lambda[i], alpha, beta[i], delta);
352 const double B = -a2sigma - n2[i]*k1/k2;
353 const double A = k1 * std::pow(B+a2sigma, n2[i]);
354
355 output[i] = A * std::pow(B+d, -n2[i]);
356 }
357 else {
358 output[i] = LogEval(d, lambda[i], alpha, beta[i], delta);
359 }
360 }
361
362 if ((!zetaIsAlwaysLargerZero && zetaIsAlwaysZero) || (zeta[i] == 0. && lambda[i] < 0.)) {
363 const double delta = sigma[i];
364
365 if (d < -asigma ) {
366 const double cons1 = std::exp(-beta[i]*asigma);
367 const double phi = 1. + a[i] * a[i];
368 const double k1 = cons1 * std::pow(phi, lambda[i]-0.5);
369 const double k2 = beta[i]*k1 - cons1*(lambda[i]-0.5) * std::pow(phi, lambda[i]-1.5) * 2.*a[i]/delta;
370 const double B = -asigma + n[i]*k1/k2;
371 const double A = k1*std::pow(B+asigma, n[i]);
372
373 output[i] = A*std::pow(B-d, -n[i]);
374 }
375 else if (d > a2sigma) {
376 const double cons1 = std::exp(beta[i]*a2sigma);
377 const double phi = 1. + a2[i]*a2[i];
378 const double k1 = cons1 * std::pow(phi, lambda[i]-0.5);
379 const double k2 = beta[i]*k1 + cons1*(lambda[i]-0.5) * std::pow(phi, lambda[i]-1.5) * 2.*a2[i]/delta;
380 const double B = -a2sigma - n2[i]*k1/k2;
381 const double A = k1*std::pow(B+a2sigma, n2[i]);
382
383 output[i] = A*std::pow(B+d, -n2[i]);
384 }
385 else {
386 output[i] = std::exp(beta[i]*d) * std::pow(1. + d*d/(delta*delta), lambda[i] - 0.5);
387 }
388 }
389 }
390}
391
392template<bool right>
393std::pair<double, double> computeAB_zetaNZero(double asigma, double lambda, double alpha, double beta, double delta, double n) {
394 const double k1 = LogEval(right ? asigma : -asigma, lambda, alpha, beta, delta);
395 const double k2 = diff_eval(right ? asigma : -asigma, lambda, alpha, beta, delta);
396 const double B = -asigma + (right ? -1 : 1.) * n*k1/k2;
397 const double A = k1 * std::pow(B+asigma, n);
398
399 return {A, B};
400}
401
402template<bool right>
403std::pair<double, double> computeAB_zetaZero(double beta, double asigma, double a, double lambda, double delta, double n) {
404 const double cons1 = std::exp((right ? 1. : -1.) * beta*asigma);
405 const double phi = 1. + a * a;
406 const double k1 = cons1 * std::pow(phi, lambda-0.5);
407 const double k2 = beta*k1 + (right ? 1. : -1) * cons1*(lambda-0.5) * std::pow(phi, lambda-1.5) * 2.*a/delta;
408 const double B = -asigma + (right ? -1. : 1.) * n*k1/k2;
409 const double A = k1*std::pow(B+asigma, n);
410
411 return {A, B};
412}
413
415//////////////////////////////////////////////////////////////////////////////////////////
416/// A specialised compute function where x is an observable, and all parameters are used as
417/// parameters. Since many things can be calculated outside of the loop, it is faster.
418void compute(std::span<double> output, std::span<const double> x,
419 BracketAdapter<double> lambda, BracketAdapter<double> zeta, BracketAdapter<double> beta,
420 BracketAdapter<double> sigma, BracketAdapter<double> mu,
421 BracketAdapter<double> a, BracketAdapter<double> n, BracketAdapter<double> a2, BracketAdapter<double> n2) {
422 const auto N = output.size();
423
424 const double cons0 = std::sqrt(zeta);
425 const double asigma = a*sigma;
426 const double a2sigma = a2*sigma;
427
428 if (zeta > 0.) {
429 const double phi = besselK(lambda+1., zeta) / besselK(lambda, zeta);
430 const double cons1 = sigma/std::sqrt(phi);
431 const double alpha = cons0/cons1;
432 const double delta = cons0*cons1;
433
434 const auto AB_l = computeAB_zetaNZero<false>(asigma, lambda, alpha, beta, delta, n);
435 const auto AB_r = computeAB_zetaNZero<true>(a2sigma, lambda, alpha, beta, delta, n2);
436
437 for (std::size_t i = 0; i < N; ++i) {
438 const double d = x[i] - mu[i];
439
440 if (d < -asigma){
441 output[i] = AB_l.first * std::pow(AB_l.second - d, -n);
442 }
443 else if (d>a2sigma) {
444 output[i] = AB_r.first * std::pow(AB_r.second + d, -n2);
445 }
446 else {
447 output[i] = LogEval(d, lambda, alpha, beta, delta);
448 }
449 }
450 } else if (zeta == 0. && lambda < 0.) {
451 const double delta = sigma;
452
453 const auto AB_l = computeAB_zetaZero<false>(beta, asigma, a, lambda, delta, n);
454 const auto AB_r = computeAB_zetaZero<true>(beta, a2sigma, a2, lambda, delta, n2);
455
456 for (std::size_t i = 0; i < N; ++i) {
457 const double d = x[i] - mu[i];
458
459 if (d < -asigma ) {
460 output[i] = AB_l.first*std::pow(AB_l.second - d, -n);
461 }
462 else if (d > a2sigma) {
463 output[i] = AB_r.first * std::pow(AB_r.second + d, -n2);
464 }
465 else {
466 output[i] = std::exp(beta*d) * std::pow(1. + d*d/(delta*delta), lambda - 0.5);
467 }
468 }
469 }
470}
471
472}
473
475{
476 using namespace RooBatchCompute;
477
478 auto x = ctx.at(_x);
479 auto lambda = ctx.at(_lambda);
480 auto zeta = ctx.at(_zeta);
481 auto beta = ctx.at(_beta);
482 auto sig = ctx.at(_sigma);
483 auto mu = ctx.at(_mu);
484 auto a = ctx.at(_a);
485 auto n = ctx.at(_n);
486 auto a2 = ctx.at(_a2);
487 auto n2 = ctx.at(_n2);
488
489 size_t paramSizeSum=0;
490 for (const auto& i:{lambda, zeta, beta, sig, mu, a, n, a2, n2}) {
491 paramSizeSum += i.size();
492 }
493
494 // Run high performance compute if only x has multiple values
495 if (x.size()>1 && paramSizeSum==9) {
496 compute(ctx.output(), x,
501 } else {
502 compute(ctx.output(), BracketAdapterWithMask(_x, x),
508 }
509}
510
511
512/* Analytical integrals need testing.
513
514Int_t RooHypatia2::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char*) const
515{
516 if (matchArgs(allVars, analVars, _x)
517 && _beta == 0. && _beta.arg().isConstant()
518 && _zeta == 0. && _zeta.arg().isConstant()
519 && _lambda.max() < 0.) return 1;
520 return 0 ;
521}
522
523
524
525double RooHypatia2::analyticalIntegral(Int_t code, const char* rangeName) const
526{
527 if (_beta != 0. || _zeta != 0. || _lambda >= 0) {
528 auto& logstream = coutF(Integration) << "When the PDF " << GetName()
529 << " was constructed, beta,zeta were 0, lambda<0 and all three were constant.\n"
530 << "This allowed for analytic integration, but now, numeric integration would be required.\n"
531 << "These parameters must either be kept constant, or be floating from the beginning. "
532 << " Terminating fit ..."
533 << std::endl;
534 RooArgSet vars;
535 vars.add(_beta.arg());
536 vars.add(_zeta.arg());
537 vars.add(_lambda.arg());
538 vars.printStream(logstream, vars.defaultPrintContents(nullptr), RooPrintable::kVerbose);
539 throw std::runtime_error("RooHypatia2 cannot be integrated analytically since parameters changed.");
540 }
541
542 // The formulae to follow still use beta and zeta to facilitate comparisons with the
543 // evaluate function. Since beta and zeta are zero, all relevant terms will be disabled
544 // by defining these two constexpr:
545 constexpr double beta = 0.;
546 constexpr double cons1 = 1.;
547
548 if (code != 1) {
549 throw std::logic_error("Trying to compute analytic integral with unknown configuration.");
550 }
551
552 const double asigma = _a * _sigma;
553 const double a2sigma = _a2 * _sigma;
554 const double d0 = _x.min(rangeName) - _mu;
555 const double d1 = _x.max(rangeName) - _mu;
556
557
558 double delta;
559 if (_lambda <= -1.0) {
560 delta = _sigma * sqrt(-2. + -2.*_lambda);
561 }
562 else {
563 delta = _sigma;
564 }
565 const double deltaSq = delta*delta;
566
567 if ((d0 > -asigma) && (d1 < a2sigma)) {
568 return stIntegral(d1, delta, _lambda) - stIntegral(d0, delta, _lambda);
569 }
570
571 if (d0 > a2sigma) {
572 const double phi = 1. + a2sigma*a2sigma/deltaSq;
573 const double k1 = cons1*std::pow(phi,_lambda-0.5);
574 const double k2 = beta*k1+ cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2.*a2sigma/deltaSq;
575 const double B = -a2sigma - _n2*k1/k2;
576 const double A = k1*std::pow(B+a2sigma,_n2);
577 return A*(std::pow(B+d1,1-_n2)/(1-_n2) -std::pow(B+d0,1-_n2)/(1-_n2) ) ;
578
579 }
580
581 if (d1 < -asigma) {
582 const double phi = 1. + asigma*asigma/deltaSq;
583 const double k1 = cons1*std::pow(phi,_lambda-0.5);
584 const double k2 = beta*k1- cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2*asigma/deltaSq;
585 const double B = -asigma + _n*k1/k2;
586 const double A = k1*std::pow(B+asigma,_n);
587 const double I0 = A*std::pow(B-d0,1-_n)/(_n-1);
588 const double I1 = A*std::pow(B-d1,1-_n)/(_n-1);
589 return I1 - I0;
590 }
591
592
593 double I0;
594 double I1a = 0;
595 double I1b = 0;
596 if (d0 <-asigma) {
597 const double phi = 1. + asigma*asigma/deltaSq;
598 const double k1 = cons1*std::pow(phi,_lambda-0.5);
599 const double k2 = beta*k1- cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2*asigma/deltaSq;
600 const double B = -asigma + _n*k1/k2;
601 const double A = k1*std::pow(B+asigma,_n);
602 I0 = A*std::pow(B-d0,1-_n)/(_n-1);
603 I1a = A*std::pow(B+asigma,1-_n)/(_n-1) - stIntegral(-asigma, delta, _lambda);
604 }
605 else {
606 I0 = stIntegral(d0, delta, _lambda);
607 }
608
609 if (d1 > a2sigma) {
610 const double phi = 1. + a2sigma*a2sigma/deltaSq;
611 const double k1 = cons1*std::pow(phi,_lambda-0.5);
612 const double k2 = beta*k1+ cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2.*a2sigma/deltaSq;
613 const double B = -a2sigma - _n2*k1/k2;
614 const double A = k1*std::pow(B+a2sigma,_n2);
615 I1b = A*(std::pow(B+d1,1-_n2)/(1-_n2) -std::pow(B+a2sigma,1-_n2)/(1-_n2) ) - stIntegral(d1, delta, _lambda) + stIntegral(a2sigma,delta, _lambda) ;
616 }
617
618 const double I1 = stIntegral(d1, delta, _lambda) + I1a + I1b;
619 return I1 - I0;
620}
621
622*/
#define d(i)
Definition RSha256.hxx:102
#define a(i)
Definition RSha256.hxx:99
#define coutE(a)
#define N
char name[80]
Definition TGX11.cxx:110
bool isConstant() const
Check if the "Constant" attribute is set.
Definition RooAbsArg.h:361
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
Little adapter that gives a bracket operator to types that don't have one.
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
std::span< double > output()
RooHypatia2 is the two-sided version of the Hypatia distribution described in https://arxiv....
Definition RooHypatia2.h:25
RooRealProxy _beta
Definition RooHypatia2.h:46
RooRealProxy _n
Definition RooHypatia2.h:50
RooRealProxy _a
Definition RooHypatia2.h:49
void doEval(RooFit::EvalContext &) const override
Base function for computing multiple values of a RooAbsReal.
RooRealProxy _zeta
Definition RooHypatia2.h:45
RooRealProxy _n2
Definition RooHypatia2.h:52
RooRealProxy _a2
Definition RooHypatia2.h:51
RooRealProxy _mu
Definition RooHypatia2.h:48
RooRealProxy _sigma
Definition RooHypatia2.h:47
RooRealProxy _lambda
Definition RooHypatia2.h:44
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooRealProxy _x
Definition RooHypatia2.h:43
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
double beta(double x, double y)
Calculates the beta function.
double cyl_bessel_k(double nu, double x)
Calculates the modified Bessel functions of the second kind (also called irregular modified (cylindri...
const Double_t sigma
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
double gamma(double x)
Namespace for dispatching RooFit computations to various backends.
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan extraArgs={})
void checkRangeOfParameters(const RooAbsReal *callingClass, std::initializer_list< const RooAbsReal * > pars, double min=-std::numeric_limits< double >::max(), double max=std::numeric_limits< double >::max(), bool limitsInAllowedRange=false, std::string const &extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition TMath.cxx:353
constexpr Double_t TwoPi()
Definition TMath.h:44
TLine l
Definition textangle.C:4
static void output()