Logo ROOT  
Reference Guide
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/P2VV/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#include "RunContext.h"
97
98#include "TMath.h"
99#include "Math/SpecFunc.h"
100#include "ROOT/RConfig.hxx"
101
102#include <cmath>
103#include <algorithm>
104
105///////////////////////////////////////////////////////////////////////////////////////////
106/// Construct a new Hypatia2 PDF.
107/// \param[in] name Name of this instance.
108/// \param[in] title Title (for plotting)
109/// \param[in] x The variable of this distribution
110/// \param[in] lambda Shape parameter. Note that \f$ \lambda < 0 \f$ is required if \f$ \zeta = 0 \f$.
111/// \param[in] zeta Shape parameter (\f$ \zeta >= 0 \f$).
112/// \param[in] beta Asymmetry parameter \f$ \beta \f$. Symmetric case is \f$ \beta = 0 \f$,
113/// choose values close to zero.
114/// \param[in] argSigma Width parameter. If \f$ \beta = 0, \ \sigma \f$ is the RMS width.
115/// \param[in] mu Location parameter. Shifts the distribution left/right.
116/// \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$,
117/// the tail region is to the left of \f$ x = \mu - 1 \f$, so a should be positive.
118/// \param[in] n Shape parameter of left tail (\f$ n \ge 0 \f$). With \f$ n = 0 \f$, the function is constant.
119/// \param[in] a2 Start of right tail.
120/// \param[in] n2 Shape parameter of right tail (\f$ n2 \ge 0 \f$). With \f$ n2 = 0 \f$, the function is constant.
121RooHypatia2::RooHypatia2(const char *name, const char *title, RooAbsReal& x, RooAbsReal& lambda,
122 RooAbsReal& zeta, RooAbsReal& beta, RooAbsReal& argSigma, RooAbsReal& mu, RooAbsReal& a,
123 RooAbsReal& n, RooAbsReal& a2, RooAbsReal& n2) :
124 RooAbsPdf(name, title),
125 _x("x", "x", this, x),
126 _lambda("lambda", "Lambda", this, lambda),
127 _zeta("zeta", "zeta", this, zeta),
128 _beta("beta", "Asymmetry parameter beta", this, beta),
129 _sigma("sigma", "Width parameter sigma", this, argSigma),
130 _mu("mu", "Location parameter mu", this, mu),
131 _a("a", "Left tail location a", this, a),
132 _n("n", "Left tail parameter n", this, n),
133 _a2("a2", "Right tail location a2", this, a2),
134 _n2("n2", "Right tail parameter n2", this, n2)
135{
136 RooHelpers::checkRangeOfParameters(this, {&argSigma}, 0.);
137 RooHelpers::checkRangeOfParameters(this, {&zeta, &n, &n2, &a, &a2}, 0., std::numeric_limits<double>::max(), true);
138 if (zeta.getVal() == 0. && zeta.isConstant()) {
139 RooHelpers::checkRangeOfParameters(this, {&lambda}, -std::numeric_limits<double>::max(), 0., false,
140 std::string("Lambda needs to be negative when ") + _zeta.GetName() + " is zero.");
141 }
142
143#ifndef R__HAS_MATHMORE
144 throw std::logic_error("RooHypatia2 needs ROOT with mathmore enabled to access gsl functions.");
145#endif
146}
147
148
149///////////////////////////////////////////////////////////////////////////////////////////
150/// Copy a new Hypatia2 PDF.
151/// \param[in] other Original to copy from.
152/// \param[in] name Optional new name.
153RooHypatia2::RooHypatia2(const RooHypatia2& other, const char* name) :
154 RooAbsPdf(other, name),
155 _x("x", this, other._x),
156 _lambda("lambda", this, other._lambda),
157 _zeta("zeta", this, other._zeta),
158 _beta("beta", this, other._beta),
159 _sigma("sigma", this, other._sigma),
160 _mu("mu", this, other._mu),
161 _a("a", this, other._a),
162 _n("n", this, other._n),
163 _a2("a2", this, other._a2),
164 _n2("n2", this, other._n2)
165{
166#ifndef R__HAS_MATHMORE
167 throw std::logic_error("RooHypatia2 needs ROOT with mathmore enabled to access gsl functions.");
168#endif
169}
170
171namespace {
172const double sq2pi_inv = 1./std::sqrt(TMath::TwoPi());
173const double logsq2pi = std::log(std::sqrt(TMath::TwoPi()));
174const double ln2 = std::log(2.);
175
176double low_x_BK(double nu, double x){
177 return TMath::Gamma(nu)*std::pow(2., nu-1.)*std::pow(x, -nu);
178}
179
180double low_x_LnBK(double nu, double x){
181 return std::log(TMath::Gamma(nu)) + ln2*(nu-1.) - std::log(x) * nu;
182}
183
184double besselK(double ni, double x) {
185 const double nu = std::fabs(ni);
186 if ((x < 1.e-06 && nu > 0.) ||
187 (x < 1.e-04 && nu > 0. && nu < 55.) ||
188 (x < 0.1 && nu >= 55.) )
189 return low_x_BK(nu, x);
190
191#ifdef R__HAS_MATHMORE
192 return ROOT::Math::cyl_bessel_k(nu, x);
193#else
194 return std::numeric_limits<double>::signaling_NaN();
195#endif
196
197}
198
199double LnBesselK(double ni, double x) {
200 const double nu = std::fabs(ni);
201 if ((x < 1.e-06 && nu > 0.) ||
202 (x < 1.e-04 && nu > 0. && nu < 55.) ||
203 (x < 0.1 && nu >= 55.) )
204 return low_x_LnBK(nu, x);
205
206#ifdef R__HAS_MATHMORE
208#else
209 return std::numeric_limits<double>::signaling_NaN();
210#endif
211}
212
213
214double LogEval(double d, double l, double alpha, double beta, double delta) {
215 const double gamma = alpha;//std::sqrt(alpha*alpha-beta*beta);
216 const double dg = delta*gamma;
217 const double thing = delta*delta + d*d;
218 const double logno = l*std::log(gamma/delta) - logsq2pi - LnBesselK(l, dg);
219
220 return std::exp(logno + beta*d
221 + (0.5-l)*(std::log(alpha)-0.5*std::log(thing))
222 + LnBesselK(l-0.5, alpha*std::sqrt(thing)) );// + std::log(std::fabs(beta)+0.0001) );
223
224}
225
226
227double diff_eval(double d, double l, double alpha, double beta, double delta){
228 const double gamma = alpha;
229 const double dg = delta*gamma;
230
231 const double thing = delta*delta + d*d;
232 const double sqrthing = std::sqrt(thing);
233 const double alphasq = alpha*sqrthing;
234 const double no = std::pow(gamma/delta, l)/besselK(l, dg)*sq2pi_inv;
235 const double ns1 = 0.5-l;
236
237 return no * std::pow(alpha, ns1) * std::pow(thing, l/2. - 1.25)
238 * ( -d * alphasq * (besselK(l - 1.5, alphasq)
239 + besselK(l + 0.5, alphasq))
240 + (2.*(beta*thing + d*l) - d) * besselK(ns1, alphasq) )
241 * std::exp(beta*d) * 0.5;
242}
243
244/*
245double Gauss2F1(double a, double b, double c, double x){
246 if (fabs(x) <= 1.) {
247 return ROOT::Math::hyperg(a, b, c, x);
248 } else {
249 return ROOT::Math::hyperg(c-a, b, c, 1-1/(1-x))/std::pow(1-x, b);
250 }
251}
252
253double stIntegral(double d1, double delta, double l){
254 return d1 * Gauss2F1(0.5, 0.5-l, 3./2, -d1*d1/(delta*delta));
255}
256*/
257}
258
260{
261 const double d = _x-_mu;
262 const double cons0 = std::sqrt(_zeta);
263 const double asigma = _a*_sigma;
264 const double a2sigma = _a2*_sigma;
265 const double beta = _beta;
266 double out = 0.;
267
268 if (_zeta > 0.) {
269 // careful if zeta -> 0. You can implement a function for the ratio,
270 // but careful again that |nu + 1 | != |nu| + 1 so you have to deal with the signs
271 const double phi = besselK(_lambda+1., _zeta) / besselK(_lambda, _zeta);
272 const double cons1 = _sigma/std::sqrt(phi);
273 const double alpha = cons0/cons1;
274 const double delta = cons0*cons1;
275
276 if (d < -asigma){
277 const double k1 = LogEval(-asigma, _lambda, alpha, beta, delta);
278 const double k2 = diff_eval(-asigma, _lambda, alpha, beta, delta);
279 const double B = -asigma + _n*k1/k2;
280 const double A = k1 * std::pow(B+asigma, _n);
281
282 out = A * std::pow(B-d, -_n);
283 }
284 else if (d>a2sigma) {
285 const double k1 = LogEval(a2sigma, _lambda, alpha, beta, delta);
286 const double k2 = diff_eval(a2sigma, _lambda, alpha, beta, delta);
287 const double B = -a2sigma - _n2*k1/k2;
288 const double A = k1 * std::pow(B+a2sigma, _n2);
289
290 out = A * std::pow(B+d, -_n2);
291 }
292 else {
293 out = LogEval(d, _lambda, alpha, beta, delta);
294 }
295 }
296 else if (_zeta < 0.) {
297 coutE(Eval) << "The parameter " << _zeta.GetName() << " of the RooHypatia2 " << GetName() << " cannot be < 0." << std::endl;
298 }
299 else if (_lambda < 0.) {
300 const double delta = _sigma;
301
302 if (d < -asigma ) {
303 const double cons1 = std::exp(-beta*asigma);
304 const double phi = 1. + _a * _a;
305 const double k1 = cons1 * std::pow(phi, _lambda-0.5);
306 const double k2 = beta*k1 - cons1*(_lambda-0.5) * std::pow(phi, _lambda-1.5) * 2.*_a/delta;
307 const double B = -asigma + _n*k1/k2;
308 const double A = k1*std::pow(B+asigma, _n);
309
310 out = A*std::pow(B-d, -_n);
311 }
312 else if (d > a2sigma) {
313 const double cons1 = std::exp(beta*a2sigma);
314 const double phi = 1. + _a2*_a2;
315 const double k1 = cons1 * std::pow(phi, _lambda-0.5);
316 const double k2 = beta*k1 + cons1*(_lambda-0.5) * std::pow(phi, _lambda-1.5) * 2.*_a2/delta;
317 const double B = -a2sigma - _n2*k1/k2;
318 const double A = k1*std::pow(B+a2sigma, _n2);
319
320 out = A*std::pow(B+d, -_n2);
321 }
322 else {
323 out = std::exp(beta*d) * std::pow(1. + d*d/(delta*delta), _lambda - 0.5);
324 }
325 }
326 else {
327 coutE(Eval) << "zeta = 0 only supported for lambda < 0. lambda = " << double(_lambda) << std::endl;
328 }
329
330 return out;
331}
332
333namespace {
334//////////////////////////////////////////////////////////////////////////////////////////
335/// The generic compute function that recalculates everything for every loop iteration.
336/// This is only needed in the rare case where a parameter is used as an observable.
337template<typename Tx, typename Tl, typename Tz, typename Tb, typename Ts, typename Tm, typename Ta, typename Tn,
338typename Ta2, typename Tn2>
339void compute(RooSpan<double> output, Tx x, Tl lambda, Tz zeta, Tb beta, Ts sigma, Tm mu, Ta a, Tn n, Ta2 a2, Tn2 n2) {
340 const auto N = output.size();
341 const bool zetaIsAlwaysLargerZero = !zeta.isBatch() && zeta[0] > 0.;
342 const bool zetaIsAlwaysZero = !zeta.isBatch() && zeta[0] == 0.;
343
344 for (std::size_t i = 0; i < N; ++i) {
345
346 const double d = x[i] - mu[i];
347 const double cons0 = std::sqrt(zeta[i]);
348 const double asigma = a[i]*sigma[i];
349 const double a2sigma = a2[i]*sigma[i];
350// const double beta = beta[i];
351
352 if (zetaIsAlwaysLargerZero || zeta[i] > 0.) {
353 const double phi = besselK(lambda[i]+1., zeta[i]) / besselK(lambda[i], zeta[i]);
354 const double cons1 = sigma[i]/std::sqrt(phi);
355 const double alpha = cons0/cons1;
356 const double delta = cons0*cons1;
357
358 if (d < -asigma){
359 const double k1 = LogEval(-asigma, lambda[i], alpha, beta[i], delta);
360 const double k2 = diff_eval(-asigma, lambda[i], alpha, beta[i], delta);
361 const double B = -asigma + n[i]*k1/k2;
362 const double A = k1 * std::pow(B+asigma, n[i]);
363
364 output[i] = A * std::pow(B-d, -n[i]);
365 }
366 else if (d>a2sigma) {
367 const double k1 = LogEval(a2sigma, lambda[i], alpha, beta[i], delta);
368 const double k2 = diff_eval(a2sigma, lambda[i], alpha, beta[i], delta);
369 const double B = -a2sigma - n2[i]*k1/k2;
370 const double A = k1 * std::pow(B+a2sigma, n2[i]);
371
372 output[i] = A * std::pow(B+d, -n2[i]);
373 }
374 else {
375 output[i] = LogEval(d, lambda[i], alpha, beta[i], delta);
376 }
377 }
378
379 if ((!zetaIsAlwaysLargerZero && zetaIsAlwaysZero) || (zeta[i] == 0. && lambda[i] < 0.)) {
380 const double delta = sigma[i];
381
382 if (d < -asigma ) {
383 const double cons1 = std::exp(-beta[i]*asigma);
384 const double phi = 1. + a[i] * a[i];
385 const double k1 = cons1 * std::pow(phi, lambda[i]-0.5);
386 const double k2 = beta[i]*k1 - cons1*(lambda[i]-0.5) * std::pow(phi, lambda[i]-1.5) * 2.*a[i]/delta;
387 const double B = -asigma + n[i]*k1/k2;
388 const double A = k1*std::pow(B+asigma, n[i]);
389
390 output[i] = A*std::pow(B-d, -n[i]);
391 }
392 else if (d > a2sigma) {
393 const double cons1 = std::exp(beta[i]*a2sigma);
394 const double phi = 1. + a2[i]*a2[i];
395 const double k1 = cons1 * std::pow(phi, lambda[i]-0.5);
396 const double k2 = beta[i]*k1 + cons1*(lambda[i]-0.5) * std::pow(phi, lambda[i]-1.5) * 2.*a2[i]/delta;
397 const double B = -a2sigma - n2[i]*k1/k2;
398 const double A = k1*std::pow(B+a2sigma, n2[i]);
399
400 output[i] = A*std::pow(B+d, -n2[i]);
401 }
402 else {
403 output[i] = std::exp(beta[i]*d) * std::pow(1. + d*d/(delta*delta), lambda[i] - 0.5);
404 }
405 }
406 }
407}
408
409template<bool right>
410std::pair<double, double> computeAB_zetaNZero(double asigma, double lambda, double alpha, double beta, double delta, double n) {
411 const double k1 = LogEval(right ? asigma : -asigma, lambda, alpha, beta, delta);
412 const double k2 = diff_eval(right ? asigma : -asigma, lambda, alpha, beta, delta);
413 const double B = -asigma + (right ? -1 : 1.) * n*k1/k2;
414 const double A = k1 * std::pow(B+asigma, n);
415
416 return {A, B};
417}
418
419template<bool right>
420std::pair<double, double> computeAB_zetaZero(double beta, double asigma, double a, double lambda, double delta, double n) {
421 const double cons1 = std::exp((right ? 1. : -1.) * beta*asigma);
422 const double phi = 1. + a * a;
423 const double k1 = cons1 * std::pow(phi, lambda-0.5);
424 const double k2 = beta*k1 + (right ? 1. : -1) * cons1*(lambda-0.5) * std::pow(phi, lambda-1.5) * 2.*a/delta;
425 const double B = -asigma + (right ? -1. : 1.) * n*k1/k2;
426 const double A = k1*std::pow(B+asigma, n);
427
428 return {A, B};
429}
430
432//////////////////////////////////////////////////////////////////////////////////////////
433/// A specialised compute function where x is an observable, and all parameters are used as
434/// parameters. Since many things can be calculated outside of the loop, it is faster.
436 BracketAdapter<double> lambda, BracketAdapter<double> zeta, BracketAdapter<double> beta,
437 BracketAdapter<double> sigma, BracketAdapter<double> mu,
438 BracketAdapter<double> a, BracketAdapter<double> n, BracketAdapter<double> a2, BracketAdapter<double> n2) {
439 const auto N = output.size();
440
441 const double cons0 = std::sqrt(zeta);
442 const double asigma = a*sigma;
443 const double a2sigma = a2*sigma;
444
445 if (zeta > 0.) {
446 const double phi = besselK(lambda+1., zeta) / besselK(lambda, zeta);
447 const double cons1 = sigma/std::sqrt(phi);
448 const double alpha = cons0/cons1;
449 const double delta = cons0*cons1;
450
451 const auto AB_l = computeAB_zetaNZero<false>(asigma, lambda, alpha, beta, delta, n);
452 const auto AB_r = computeAB_zetaNZero<true>(a2sigma, lambda, alpha, beta, delta, n2);
453
454 for (std::size_t i = 0; i < N; ++i) {
455 const double d = x[i] - mu[i];
456
457 if (d < -asigma){
458 output[i] = AB_l.first * std::pow(AB_l.second - d, -n);
459 }
460 else if (d>a2sigma) {
461 output[i] = AB_r.first * std::pow(AB_r.second + d, -n2);
462 }
463 else {
464 output[i] = LogEval(d, lambda, alpha, beta, delta);
465 }
466 }
467 } else if (zeta == 0. && lambda < 0.) {
468 const double delta = sigma;
469
470 const auto AB_l = computeAB_zetaZero<false>(beta, asigma, a, lambda, delta, n);
471 const auto AB_r = computeAB_zetaZero<true>(beta, a2sigma, a2, lambda, delta, n2);
472
473 for (std::size_t i = 0; i < N; ++i) {
474 const double d = x[i] - mu[i];
475
476 if (d < -asigma ) {
477 output[i] = AB_l.first*std::pow(AB_l.second - d, -n);
478 }
479 else if (d > a2sigma) {
480 output[i] = AB_r.first * std::pow(AB_r.second + d, -n2);
481 }
482 else {
483 output[i] = std::exp(beta*d) * std::pow(1. + d*d/(delta*delta), lambda - 0.5);
484 }
485 }
486 }
487}
488
489}
490
492 using namespace RooBatchCompute;
493
494 auto x = _x->getValues(evalData, normSet);
495 auto lambda = _lambda->getValues(evalData, normSet);
496 auto zeta = _zeta->getValues(evalData, normSet);
497 auto beta = _beta->getValues(evalData, normSet);
498 auto sig = _sigma->getValues(evalData, normSet);
499 auto mu = _mu->getValues(evalData, normSet);
500 auto a = _a->getValues(evalData, normSet);
501 auto n = _n->getValues(evalData, normSet);
502 auto a2 = _a2->getValues(evalData, normSet);
503 auto n2 = _n2->getValues(evalData, normSet);
504
505 size_t paramSizeSum=0, batchSize = x.size();
506 for (const auto& i:{lambda, zeta, beta, sig, mu, a, n, a2, n2}) {
507 paramSizeSum += i.size();
508 batchSize = std::max(batchSize, i.size());
509 }
510 RooSpan<double> output = evalData.makeBatch(this, batchSize);
511
512 // Run high performance compute if only x has multiple values
513 if (x.size()>1 && paramSizeSum==9) {
514 compute(output, x,
519 } else {
526 }
527 return output;
528}
529
530
531/* Analytical integrals need testing.
532
533Int_t RooHypatia2::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char*) const
534{
535 if (matchArgs(allVars, analVars, _x)
536 && _beta == 0. && _beta.arg().isConstant()
537 && _zeta == 0. && _zeta.arg().isConstant()
538 && _lambda.max() < 0.) return 1;
539 return 0 ;
540}
541
542
543
544double RooHypatia2::analyticalIntegral(Int_t code, const char* rangeName) const
545{
546 if (_beta != 0. || _zeta != 0. || _lambda >= 0) {
547 auto& logstream = coutF(Integration) << "When the PDF " << GetName()
548 << " was constructed, beta,zeta were 0, lambda<0 and all three were constant.\n"
549 << "This allowed for analytic integration, but now, numeric integration would be required.\n"
550 << "These parameters must either be kept constant, or be floating from the beginning. "
551 << " Terminating fit ..."
552 << std::endl;
553 RooArgSet vars;
554 vars.add(_beta.arg());
555 vars.add(_zeta.arg());
556 vars.add(_lambda.arg());
557 vars.printStream(logstream, vars.defaultPrintContents(nullptr), RooPrintable::kVerbose);
558 throw std::runtime_error("RooHypatia2 cannot be integrated analytically since parameters changed.");
559 }
560
561 // The formulae to follow still use beta and zeta to facilitate comparisons with the
562 // evaluate function. Since beta and zeta are zero, all relevant terms will be disabled
563 // by defining these two constexpr:
564 constexpr double beta = 0.;
565 constexpr double cons1 = 1.;
566
567 if (code != 1) {
568 throw std::logic_error("Trying to compute analytic integral with unknown configuration.");
569 }
570
571 const double asigma = _a * _sigma;
572 const double a2sigma = _a2 * _sigma;
573 const double d0 = _x.min(rangeName) - _mu;
574 const double d1 = _x.max(rangeName) - _mu;
575
576
577 double delta;
578 if (_lambda <= -1.0) {
579 delta = _sigma * sqrt(-2. + -2.*_lambda);
580 }
581 else {
582 delta = _sigma;
583 }
584 const double deltaSq = delta*delta;
585
586 if ((d0 > -asigma) && (d1 < a2sigma)) {
587 return stIntegral(d1, delta, _lambda) - stIntegral(d0, delta, _lambda);
588 }
589
590 if (d0 > a2sigma) {
591 const double phi = 1. + a2sigma*a2sigma/deltaSq;
592 const double k1 = cons1*std::pow(phi,_lambda-0.5);
593 const double k2 = beta*k1+ cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2.*a2sigma/deltaSq;
594 const double B = -a2sigma - _n2*k1/k2;
595 const double A = k1*std::pow(B+a2sigma,_n2);
596 return A*(std::pow(B+d1,1-_n2)/(1-_n2) -std::pow(B+d0,1-_n2)/(1-_n2) ) ;
597
598 }
599
600 if (d1 < -asigma) {
601 const double phi = 1. + asigma*asigma/deltaSq;
602 const double k1 = cons1*std::pow(phi,_lambda-0.5);
603 const double k2 = beta*k1- cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2*asigma/deltaSq;
604 const double B = -asigma + _n*k1/k2;
605 const double A = k1*std::pow(B+asigma,_n);
606 const double I0 = A*std::pow(B-d0,1-_n)/(_n-1);
607 const double I1 = A*std::pow(B-d1,1-_n)/(_n-1);
608 return I1 - I0;
609 }
610
611
612 double I0;
613 double I1a = 0;
614 double I1b = 0;
615 if (d0 <-asigma) {
616 const double phi = 1. + asigma*asigma/deltaSq;
617 const double k1 = cons1*std::pow(phi,_lambda-0.5);
618 const double k2 = beta*k1- cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2*asigma/deltaSq;
619 const double B = -asigma + _n*k1/k2;
620 const double A = k1*std::pow(B+asigma,_n);
621 I0 = A*std::pow(B-d0,1-_n)/(_n-1);
622 I1a = A*std::pow(B+asigma,1-_n)/(_n-1) - stIntegral(-asigma, delta, _lambda);
623 }
624 else {
625 I0 = stIntegral(d0, delta, _lambda);
626 }
627
628 if (d1 > a2sigma) {
629 const double phi = 1. + a2sigma*a2sigma/deltaSq;
630 const double k1 = cons1*std::pow(phi,_lambda-0.5);
631 const double k2 = beta*k1+ cons1*(_lambda-0.5)*std::pow(phi,_lambda-1.5)*2.*a2sigma/deltaSq;
632 const double B = -a2sigma - _n2*k1/k2;
633 const double A = k1*std::pow(B+a2sigma,_n2);
634 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) ;
635 }
636
637 const double I1 = stIntegral(d1, delta, _lambda) + I1a + I1b;
638 return I1 - I0;
639}
640
641*/
#define d(i)
Definition: RSha256.hxx:102
#define coutE(a)
Definition: RooMsgService.h:37
#define N
char name[80]
Definition: TGX11.cxx:110
bool isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:380
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:94
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
Definition: RooAbsReal.cxx:311
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
Little adapter that gives a bracket operator to types that don't have one.
RooHypatia2 is the two-sided version of the Hypatia distribution described in https://arxiv....
Definition: RooHypatia2.h:25
RooRealProxy _beta
Definition: RooHypatia2.h:47
RooRealProxy _n
Definition: RooHypatia2.h:51
RooRealProxy _a
Definition: RooHypatia2.h:50
RooRealProxy _zeta
Definition: RooHypatia2.h:46
RooRealProxy _n2
Definition: RooHypatia2.h:53
RooRealProxy _a2
Definition: RooHypatia2.h:52
RooRealProxy _mu
Definition: RooHypatia2.h:49
RooRealProxy _sigma
Definition: RooHypatia2.h:48
RooRealProxy _lambda
Definition: RooHypatia2.h:45
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooRealProxy _x
Definition: RooHypatia2.h:44
RooSpan< double > evaluateSpan(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet) const override
Evaluate this object for a batch/span of data points.
A simple container to hold a batch of data values.
Definition: RooSpan.h:34
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...
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1753
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1748
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
const Double_t sigma
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
static double B[]
static double A[]
double gamma(double x)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Namespace for dispatching RooFit computations to various backends.
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.
Definition: RooHelpers.cxx:119
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
This struct enables passing computation data around between elements of a computation graph.
Definition: RunContext.h:32
RooSpan< double > makeBatch(const RooAbsArg *owner, std::size_t size)
Create a writable batch.
Definition: RunContext.cxx:89
TLine l
Definition: textangle.C:4
TArc a
Definition: textangle.C:12
static void output()