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)
140 std::string(
"Lambda needs to be negative when ") +
_zeta.
GetName() +
" is zero.");
143#ifndef R__HAS_MATHMORE
144 throw std::logic_error(
"RooHypatia2 needs ROOT with mathmore enabled to access gsl functions.");
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)
166#ifndef R__HAS_MATHMORE
167 throw std::logic_error(
"RooHypatia2 needs ROOT with mathmore enabled to access gsl functions.");
173const double logsq2pi = std::log(std::sqrt(
TMath::TwoPi()));
174const double ln2 = std::log(2.);
176double low_x_BK(
double nu,
double x){
177 return TMath::Gamma(nu)*std::pow(2., nu-1.)*std::pow(
x, -nu);
180double low_x_LnBK(
double nu,
double x){
181 return std::log(
TMath::Gamma(nu)) + ln2*(nu-1.) - std::log(
x) * nu;
184double besselK(
double ni,
double x) {
185 const double nu = std::abs(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);
191#ifdef R__HAS_MATHMORE
194 return std::numeric_limits<double>::signaling_NaN();
199double LnBesselK(
double ni,
double x) {
200 const double nu = std::abs(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);
206#ifdef R__HAS_MATHMORE
209 return std::numeric_limits<double>::signaling_NaN();
214double LogEval(
double d,
double l,
double alpha,
double beta,
double delta) {
215 const double gamma = alpha;
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);
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)) );
227double diff_eval(
double d,
double l,
double alpha,
double beta,
double delta){
228 const double gamma = alpha;
229 const double dg = delta*
gamma;
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;
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;
262 const double cons0 = std::sqrt(
_zeta);
265 const double beta =
_beta;
272 const double cons1 =
_sigma/std::sqrt(phi);
273 const double alpha = cons0/cons1;
274 const double delta = cons0*cons1;
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);
282 out = A * std::pow(B-
d, -
_n);
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);
290 out = A * std::pow(B+
d, -
_n2);
293 out = LogEval(
d,
_lambda, alpha, beta, delta);
296 else if (
_zeta < 0.) {
297 coutE(Eval) <<
"The parameter " <<
_zeta.
GetName() <<
" of the RooHypatia2 " <<
GetName() <<
" cannot be < 0." << std::endl;
300 const double delta =
_sigma;
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);
310 out = A*std::pow(B-
d, -
_n);
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);
320 out = A*std::pow(B+
d, -
_n2);
323 out = std::exp(beta*
d) * std::pow(1. +
d*
d/(delta*delta),
_lambda - 0.5);
327 coutE(Eval) <<
"zeta = 0 only supported for lambda < 0. lambda = " <<
double(
_lambda) << std::endl;
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) {
341 const bool zetaIsAlwaysLargerZero = !zeta.isBatch() && zeta[0] > 0.;
342 const bool zetaIsAlwaysZero = !zeta.isBatch() && zeta[0] == 0.;
344 for (std::size_t i = 0; i <
N; ++i) {
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];
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;
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]);
364 output[i] = A * std::pow(B-
d, -
n[i]);
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]);
372 output[i] = A * std::pow(B+
d, -n2[i]);
375 output[i] = LogEval(
d, lambda[i], alpha, beta[i], delta);
379 if ((!zetaIsAlwaysLargerZero && zetaIsAlwaysZero) || (zeta[i] == 0. && lambda[i] < 0.)) {
380 const double delta =
sigma[i];
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]);
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]);
400 output[i] = A*std::pow(B+
d, -n2[i]);
403 output[i] = std::exp(beta[i]*
d) * std::pow(1. +
d*
d/(delta*delta), lambda[i] - 0.5);
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);
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);
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) {
441 const double cons0 = std::sqrt(zeta);
442 const double asigma =
a*
sigma;
443 const double a2sigma = a2*
sigma;
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;
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);
454 for (std::size_t i = 0; i <
N; ++i) {
455 const double d =
x[i] - mu[i];
458 output[i] = AB_l.first * std::pow(AB_l.second -
d, -
n);
460 else if (
d>a2sigma) {
461 output[i] = AB_r.first * std::pow(AB_r.second +
d, -n2);
464 output[i] = LogEval(
d, lambda, alpha, beta, delta);
467 }
else if (zeta == 0. && lambda < 0.) {
468 const double delta =
sigma;
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);
473 for (std::size_t i = 0; i <
N; ++i) {
474 const double d =
x[i] - mu[i];
477 output[i] = AB_l.first*std::pow(AB_l.second -
d, -
n);
479 else if (
d > a2sigma) {
480 output[i] = AB_r.first * std::pow(AB_r.second +
d, -n2);
483 output[i] = std::exp(beta*
d) * std::pow(1. +
d*
d/(delta*delta), lambda - 0.5);
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());
513 if (
x.size()>1 && paramSizeSum==9) {
bool isConstant() const
Check if the "Constant" attribute is set.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
RooArgSet is a container object that can hold multiple RooAbsArg objects.
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....
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
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.
const char * GetName() const override
Returns name of object.
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...
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.
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
constexpr Double_t TwoPi()
This struct enables passing computation data around between elements of a computation graph.
RooSpan< double > makeBatch(const RooAbsArg *owner, std::size_t size)
Create a writable batch.