78 :
RooAbsPdf(
name, title), x_(
"x",
"Dependent", this,
x), x0_(
"x0",
"X0", this, x0),
79 sigmaL_(
"sigmaL",
"Left Sigma", this, sigmaL), sigmaR_(
"sigmaR",
"Right Sigma", this, sigmaR),
80 alphaL_{
"alphaL",
"Left Alpha", this, alphaL}, nL_{
"nL",
"Left Order", this, nL},
81 alphaR_{std::make_unique<
RooRealProxy>(
"alphaR",
"Right Alpha", this, alphaR)},
82 nR_{std::make_unique<
RooRealProxy>(
"nR",
"Right Order", this, nR)}
106 :
RooAbsPdf(
name, title), x_(
"x",
"Dependent", this,
x), x0_(
"x0",
"X0", this, x0),
107 sigmaL_(
"sigmaL",
"Left Sigma", this, sigmaLR), sigmaR_(
"sigmaR",
"Right Sigma", this, sigmaLR),
108 alphaL_{
"alphaL",
"Left Alpha", this, alphaL}, nL_{
"nL",
"Left Order", this, nL},
109 alphaR_{std::make_unique<
RooRealProxy>(
"alphaR",
"Right Alpha", this, alphaR)},
110 nR_{std::make_unique<
RooRealProxy>(
"nR",
"Right Order", this, nR)}
133 :
RooAbsPdf(
name, title), x_(
"x",
"Dependent", this,
x), x0_(
"x0",
"X0", this, x0),
134 sigmaL_{
"sigmaL",
"Left Sigma", this, sigmaLR}, sigmaR_{
"sigmaR",
"Right Sigma", this, sigmaLR},
135 alphaL_{
"alphaL",
"Left Alpha", this, alpha},
136 nL_{
"nL",
"Left Order", this,
n}
139 alphaR_ = std::make_unique<RooRealProxy>(
"alphaR",
"Right Alpha",
this, alpha);
140 nR_ = std::make_unique<RooRealProxy>(
"nR",
"Right Order",
this,
n);
153 :
RooAbsPdf(other,
name), x_(
"x", this, other.x_), x0_(
"x0", this, other.x0_),
154 sigmaL_(
"sigmaL", this, other.sigmaL_),
155 sigmaR_(
"sigmaR", this, other.sigmaR_), alphaL_{
"alphaL", this, other.alphaL_},
156 nL_{
"nL", this, other.nL_},
157 alphaR_{other.alphaR_ ? std::make_unique<
RooRealProxy>(
"alphaR", this, *other.alphaR_) : nullptr},
158 nR_{other.nR_ ? std::make_unique<
RooRealProxy>(
"nR", this, *other.nR_) : nullptr}
166inline double evaluateCrystalBallTail(
double t,
double alpha,
double n)
168 double a = std::pow(
n / alpha,
n) * std::exp(-0.5 * alpha * alpha);
169 double b =
n / alpha - alpha;
171 return a / std::pow(
b - t,
n);
174inline double integrateGaussian(
double sigmaL,
double sigmaR,
double tmin,
double tmax)
176 constexpr double sqrtPiOver2 = 1.2533141373;
177 constexpr double sqrt2 = 1.4142135624;
179 const double sigmaMin = tmin < 0 ? sigmaL : sigmaR;
180 const double sigmaMax = tmax < 0 ? sigmaL : sigmaR;
182 return sqrtPiOver2 * (sigmaMax * std::erf(tmax / sqrt2) - sigmaMin * std::erf(tmin / sqrt2));
185inline double integrateTailLogVersion(
double sigma,
double alpha,
double n,
double tmin,
double tmax)
187 double a = std::pow(
n / alpha,
n) *
exp(-0.5 * alpha * alpha);
188 double b =
n / alpha - alpha;
193inline double integrateTailRegular(
double sigma,
double alpha,
double n,
double tmin,
double tmax)
195 double a = std::pow(
n / alpha,
n) *
exp(-0.5 * alpha * alpha);
196 double b =
n / alpha - alpha;
198 return a *
sigma / (1.0 -
n) * (1.0 / (std::pow(
b - tmin,
n - 1.0)) - 1.0 / (std::pow(
b - tmax,
n - 1.0)));
208 const double x0 =
x0_;
209 const double sigmaL = std::abs(
sigmaL_);
210 const double sigmaR = std::abs(
sigmaR_);
211 double alphaL = std::abs(
alphaL_);
213 double alphaR =
alphaR_ ? std::abs(*
alphaR_) : std::numeric_limits<double>::infinity();
214 double nR =
nR_ ? *
nR_ : 0.0;
219 std::swap(alphaL, alphaR);
223 const double t = (
x - x0) / (
x < x0 ? sigmaL : sigmaR);
226 return evaluateCrystalBallTail(t, alphaL, nL);
227 }
else if (t <= alphaR) {
228 return std::exp(-0.5 * t * t);
230 return evaluateCrystalBallTail(-t, alphaR, nR);
247 const double x0 =
x0_;
248 const double sigmaL = std::abs(
sigmaL_);
249 const double sigmaR = std::abs(
sigmaR_);
250 double alphaL = std::abs(
alphaL_);
252 double alphaR =
alphaR_ ? std::abs(*
alphaR_) : std::numeric_limits<double>::infinity();
253 double nR =
nR_ ? *
nR_ : 0.0;
258 std::swap(alphaL, alphaR);
262 constexpr double switchToLogThreshold = 1.0e-05;
266 const double tmin = (
xmin - x0) / (
xmin < x0 ? sigmaL : sigmaR);
267 const double tmax = (
xmax - x0) / (
xmax < x0 ? sigmaL : sigmaR);
271 if (tmin < -alphaL) {
272 auto integrateTailL = std::abs(nL - 1.0) < switchToLogThreshold ? integrateTailLogVersion : integrateTailRegular;
273 result += integrateTailL(sigmaL, alphaL, nL, tmin, std::min(tmax, -alphaL));
276 auto integrateTailR = std::abs(nR - 1.0) < switchToLogThreshold ? integrateTailLogVersion : integrateTailRegular;
277 result += integrateTailR(sigmaR, alphaR, nR, -tmax, std::min(-tmin, -alphaR));
279 if (tmin < alphaR && tmax > -alphaL) {
280 result += integrateGaussian(sigmaL, sigmaR, std::max(tmin, -alphaL), std::min(tmax, alphaR));
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Bool_t 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.
PDF implementing the generalized Asymmetrical Double-Sided Crystall Ball line shape.
std::unique_ptr< RooRealProxy > nR_
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
std::unique_ptr< RooRealProxy > alphaR_
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise that we know the maximum of self for given (m0,alpha,n,sigma).
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
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 extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.