126 _x(
"x",
"x",
this,
x),
127 _lambda(
"lambda",
"Lambda",
this, lambda),
129 _beta(
"beta",
"Asymmetry parameter beta",
this, beta),
130 _sigma(
"sigma",
"Width parameter sigma",
this,
argSigma),
131 _mu(
"mu",
"Location parameter mu",
this, mu),
132 _a(
"a",
"Left tail location a",
this,
a),
133 _n(
"n",
"Left tail parameter n",
this,
n),
134 _a2(
"a2",
"Right tail location a2",
this,
a2),
135 _n2(
"n2",
"Right tail parameter n2",
this,
n2)
139 if (
zeta.getVal() == 0. &&
zeta.isConstant()) {
141 std::string(
"Lambda needs to be negative when ") +
_zeta.
GetName() +
" is zero.");
168const double ln2 = std::log(2.);
179 const double nu = std::abs(
ni);
189 const double nu = std::abs(
ni);
199double LogEval(
double d,
double l,
double alpha,
double beta,
double delta) {
200 const double gamma = alpha;
202 const double thing = delta*delta +
d*
d;
205 return std::exp(
logno + beta*
d
206 + (0.5-
l)*(std::log(alpha)-0.5*std::log(
thing))
212double diff_eval(
double d,
double l,
double alpha,
double beta,
double delta){
213 const double gamma = alpha;
216 const double thing = delta*delta +
d*
d;
220 const double ns1 = 0.5-
l;
222 return no * std::pow(alpha,
ns1) * std::pow(
thing,
l/2. - 1.25)
226 * std::exp(beta*
d) * 0.5;
230 if (std::abs(
x) <= 1.) {
248 const double beta =
_beta;
265 out = A * std::pow(B-
d, -
_n);
273 out = A * std::pow(B+
d, -
_n2);
279 else if (
_zeta < 0.) {
280 coutE(Eval) <<
"The parameter " <<
_zeta.
GetName() <<
" of the RooHypatia2 " <<
GetName() <<
" cannot be < 0." << std::endl;
283 const double delta =
_sigma;
287 const double phi = 1. +
_a *
_a;
293 out = A*std::pow(B-
d, -
_n);
297 const double phi = 1. +
_a2*
_a2;
303 out = A*std::pow(B+
d, -
_n2);
306 out = std::exp(beta*
d) * std::pow(1. +
d*
d/(delta*delta),
_lambda - 0.5);
310 coutE(Eval) <<
"zeta = 0 only supported for lambda < 0. lambda = " <<
double(
_lambda) << std::endl;
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) {
327 for (std::size_t i = 0; i <
N; ++i) {
329 const double d =
x[i] - mu[i];
345 const double A =
k1 * std::pow(B+
asigma,
n[i]);
347 output[i] = A * std::pow(B-
d, -
n[i]);
363 const double delta =
sigma[i];
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;
371 const double A =
k1*std::pow(B+
asigma,
n[i]);
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;
386 output[i] = std::exp(beta[i]*
d) * std::pow(1. +
d*
d/(delta*delta), lambda[i] - 0.5);
393std::pair<double, double>
computeAB_zetaNZero(
double asigma,
double lambda,
double alpha,
double beta,
double delta,
double n) {
396 const double B = -
asigma + (right ? -1 : 1.) *
n*
k1/
k2;
397 const double A =
k1 * std::pow(B+
asigma,
n);
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);
437 for (std::size_t i = 0; i <
N; ++i) {
438 const double d =
x[i] - mu[i];
450 }
else if (
zeta == 0. && lambda < 0.) {
451 const double delta =
sigma;
456 for (std::size_t i = 0; i <
N; ++i) {
457 const double d =
x[i] - mu[i];
466 output[i] = std::exp(beta*
d) * std::pow(1. +
d*
d/(delta*delta), lambda - 0.5);
483 auto mu = ctx.
at(
_mu);
490 for (
const auto& i:{lambda,
zeta, beta, sig, mu,
a,
n,
a2,
n2}) {
528 <<
" was constructed, beta,zeta were 0, lambda<0 and all three were constant.\n"
529 <<
"This allowed for analytic integration, but now, numeric integration would be required.\n"
530 <<
"These parameters must either be kept constant, or be floating from the beginning. "
531 <<
" Terminating fit ..."
538 throw std::runtime_error(
"RooHypatia2 cannot be integrated analytically since parameters changed.");
544 constexpr double beta = 0.;
545 constexpr double cons1 = 1.;
548 throw std::logic_error(
"Trying to compute analytic integral with unknown configuration.");
564 const double deltaSq = delta*delta;
586 const double I0 = A*std::pow(B-
d0,1-
_n)/(
_n-1);
587 const double I1 = A*std::pow(B-
d1,1-
_n)/(
_n-1);
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Int_t defaultPrintContents(Option_t *opt) const override
Define default RooPrinable print options for given Print() flag string For inline printing only show ...
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Abstract interface for all probability density functions.
Abstract base class for objects that represent a real value and implements functionality common to al...
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
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.
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....
void doEval(RooFit::EvalContext &) const override
Base function for computing multiple values of a RooAbsReal.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
const T & arg() const
Return reference to object held in proxy.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
const char * GetName() const override
Returns name of object.
double beta(double x, double y)
Calculates the beta function.
double hyperg(double a, double b, double c, double x)
Calculates Gauss' hypergeometric 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 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.
constexpr Double_t TwoPi()