Logo ROOT  
Reference Guide
RooCrystalBall.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Authors: *
4 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
5 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
6 * *
7 * Copyright (c) 2000-2019, Regents of the University of California *
8 * and Stanford University. All rights reserved. *
9 * *
10 * Redistribution and use in source and binary forms, *
11 * with or without modification, are permitted according to the terms *
12 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
13 *****************************************************************************/
14
15// Authors of this class:
16// T. Skwarnicki:
17// - modify RooCBShape to Asymmetrical Double-Sided CB
18// Michael Wilkinson
19// - add to RooFit source
20// Jonas Rembser, CERN 02/2021:
21// - merging RooDSCBShape with RooSDSCBShape to RooCrystalBall
22// - implement possibility to have asymmetrical Gaussian core
23// - complete rewrite of evaluation and integral code to reduce code
24// duplication
25
26/** \class RooCrystalBall
27 \ingroup Roofit
28
29PDF implementing the generalized Asymmetrical Double-Sided Crystall Ball line shape.
30\f[
31 f(m;m_0,\sigma,\alpha_L,n_L,\alpha_R,n_R) =
32 \begin{cases}
33 A_L \cdot (B_L - \frac{m - m_0}{\sigma_L})^{-n_L}, & \mbox{for }\frac{m - m_0}{\sigma_L} < -\alpha_L \\
34 \exp \left( - \frac{1}{2} \cdot \left[ \frac{m - m_0}{\sigma_L} \right]^2 \right), & \mbox{for }\frac{m - m_0}{\sigma_L} \leq 0 \\
35 \exp \left( - \frac{1}{2} \cdot \left[ \frac{m - m_0}{\sigma_R} \right]^2 \right), & \mbox{for }\frac{m - m_0}{\sigma_R} \leq \alpha_R \\
36 A_R \cdot (B_R + \frac{m - m_0}{\sigma_R})^{-n_R}, & \mbox{otherwise}, \\
37 \end{cases}
38\f]
39times some normalization factor,
40where
41\f[
42 \begin{align}
43 A_i &= \left(\frac{n_i}{\left| \alpha_i \right|}\right)^{n_i} \cdot \exp\left(- \frac {\left| \alpha_i \right|^2}{2}\right) \\
44 B_i &= \frac{n_i}{\left| \alpha_i \right|} - \left| \alpha_i \right| \\
45 \end{align}
46\f]
47**/
48
49#include "RooCrystalBall.h"
50#include "RooAbsReal.h"
51#include "RooHelpers.h"
52#include "TError.h"
53
54#include <cmath>
55#include <limits>
56#include <memory>
57#include <utility>
58
60
61////////////////////////////////////////////////////////////////////////////////
62/// Creates the fully parametrized crystal ball shape with asymmetric Gaussian core and asymmetric tails.
63///
64/// \param name Name that identifies the PDF in computations.
65/// \param title Title for plotting.
66/// \param x The variable of the PDF.
67/// \param x0 Location parameter of the Gaussian component.
68/// \param sigmaL Width parameter of the left side of the Gaussian component.
69/// \param sigmaR Width parameter of the right side of the Gaussian component.
70/// \param alphaL Location of transition to a power law on the left, in standard deviations away from the mean.
71/// \param nL Exponent of power-law tail on the left.
72/// \param alphaR Location of transition to a power law on the right, in standard deviations away from the mean.
73/// \param nR Exponent of power-law tail on the right.
74RooCrystalBall::RooCrystalBall(const char *name, const char *title, RooAbsReal &x, RooAbsReal &x0, RooAbsReal &sigmaL,
75 RooAbsReal &sigmaR, RooAbsReal &alphaL, RooAbsReal &nL, RooAbsReal &alphaR,
76 RooAbsReal &nR)
77 : RooAbsPdf(name, title), x_("x", "Dependent", this, x), x0_("x0", "X0", this, x0),
78 sigmaL_("sigmaL", "Left Sigma", this, sigmaL), sigmaR_("sigmaR", "Right Sigma", this, sigmaR),
79 alphaL_{"alphaL", "Left Alpha", this, alphaL}, nL_{"nL", "Left Order", this, nL},
80 alphaR_{std::make_unique<RooRealProxy>("alphaR", "Right Alpha", this, alphaR)},
81 nR_{std::make_unique<RooRealProxy>("nR", "Right Order", this, nR)}
82{
83 RooHelpers::checkRangeOfParameters(this, {&sigmaL}, 0.0);
84 RooHelpers::checkRangeOfParameters(this, {&sigmaR}, 0.0);
85 RooHelpers::checkRangeOfParameters(this, {&alphaL}, 0.0);
86 RooHelpers::checkRangeOfParameters(this, {&alphaR}, 0.0);
87 RooHelpers::checkRangeOfParameters(this, {&nL}, 0.0);
88 RooHelpers::checkRangeOfParameters(this, {&nR}, 0.0);
89}
90
91////////////////////////////////////////////////////////////////////////////////
92/// Create a crystal ball shape with symmetric Gaussian core and asymmetric tails (just like `RooDSCBShape`).
93///
94/// \param name Name that identifies the PDF in computations.
95/// \param title Title for plotting.
96/// \param x The variable of the PDF.
97/// \param x0 Location parameter of the Gaussian component.
98/// \param sigmaLR Width parameter of the Gaussian component.
99/// \param alphaL Location of transition to a power law on the left, in standard deviations away from the mean.
100/// \param nL Exponent of power-law tail on the left.
101/// \param alphaR Location of transition to a power law on the right, in standard deviations away from the mean.
102/// \param nR Exponent of power-law tail on the right.
103RooCrystalBall::RooCrystalBall(const char *name, const char *title, RooAbsReal &x, RooAbsReal &x0, RooAbsReal &sigmaLR,
104 RooAbsReal &alphaL, RooAbsReal &nL, RooAbsReal &alphaR, RooAbsReal &nR)
105 : RooAbsPdf(name, title), x_("x", "Dependent", this, x), x0_("x0", "X0", this, x0),
106 sigmaL_("sigmaL", "Left Sigma", this, sigmaLR), sigmaR_("sigmaR", "Right Sigma", this, sigmaLR),
107 alphaL_{"alphaL", "Left Alpha", this, alphaL}, nL_{"nL", "Left Order", this, nL},
108 alphaR_{std::make_unique<RooRealProxy>("alphaR", "Right Alpha", this, alphaR)},
109 nR_{std::make_unique<RooRealProxy>("nR", "Right Order", this, nR)}
110{
111 RooHelpers::checkRangeOfParameters(this, {&sigmaLR}, 0.0);
112 RooHelpers::checkRangeOfParameters(this, {&alphaL}, 0.0);
113 RooHelpers::checkRangeOfParameters(this, {&alphaR}, 0.0);
114 RooHelpers::checkRangeOfParameters(this, {&nL}, 0.0);
115 RooHelpers::checkRangeOfParameters(this, {&nR}, 0.0);
116}
117
118////////////////////////////////////////////////////////////////////////////////
119/// Create a crystal ball shape with symmetric Gaussian core and only a tail on
120/// one side (just like `RooCBShape`) or two symmetric tails (like `RooSDSCBShape`).
121///
122/// \param name Name that identifies the PDF in computations.
123/// \param title Title for plotting.
124/// \param x The variable of the PDF.
125/// \param x0 Location parameter of the Gaussian component.
126/// \param sigmaLR Width parameter of the Gaussian component.
127/// \param alpha Location of transition to a power law, in standard deviations away from the mean.
128/// \param n Exponent of power-law tail.
129/// \param doubleSided Whether the tail is only on one side or on both sides
130RooCrystalBall::RooCrystalBall(const char *name, const char *title, RooAbsReal &x, RooAbsReal &x0, RooAbsReal &sigmaLR,
131 RooAbsReal &alpha, RooAbsReal &n, bool doubleSided)
132 : RooAbsPdf(name, title), x_("x", "Dependent", this, x), x0_("x0", "X0", this, x0),
133 sigmaL_{"sigmaL", "Left Sigma", this, sigmaLR}, sigmaR_{"sigmaR", "Right Sigma", this, sigmaLR},
134 alphaL_{"alphaL", "Left Alpha", this, alpha},
135 nL_{"nL", "Left Order", this, n}
136{
137 if (doubleSided) {
138 alphaR_ = std::make_unique<RooRealProxy>("alphaR", "Right Alpha", this, alpha);
139 nR_ = std::make_unique<RooRealProxy>("nR", "Right Order", this, n);
140 }
141
142 RooHelpers::checkRangeOfParameters(this, {&sigmaLR}, 0.0);
144 if (doubleSided) {
145 RooHelpers::checkRangeOfParameters(this, {&alpha}, 0.0);
146 }
147}
148
149////////////////////////////////////////////////////////////////////////////////
150/// Copy a RooCrystalBall.
152 : RooAbsPdf(other, name), x_("x", this, other.x_), x0_("x0", this, other.x0_),
153 sigmaL_("sigmaL", this, other.sigmaL_),
154 sigmaR_("sigmaR", this, other.sigmaR_), alphaL_{"alphaL", this, other.alphaL_},
155 nL_{"nL", this, other.nL_},
156 alphaR_{other.alphaR_ ? std::make_unique<RooRealProxy>("alphaR", this, *other.alphaR_) : nullptr},
157 nR_{other.nR_ ? std::make_unique<RooRealProxy>("nR", this, *other.nR_) : nullptr}
158{
159}
160
161////////////////////////////////////////////////////////////////////////////////
162
163namespace {
164
165inline double evaluateCrystalBallTail(double t, double alpha, double n)
166{
167 double a = std::pow(n / alpha, n) * std::exp(-0.5 * alpha * alpha);
168 double b = n / alpha - alpha;
169
170 return a / std::pow(b - t, n);
171}
172
173inline double integrateGaussian(double sigmaL, double sigmaR, double tmin, double tmax)
174{
175 constexpr double sqrtPiOver2 = 1.2533141373;
176 constexpr double sqrt2 = 1.4142135624;
177
178 const double sigmaMin = tmin < 0 ? sigmaL : sigmaR;
179 const double sigmaMax = tmax < 0 ? sigmaL : sigmaR;
180
181 return sqrtPiOver2 * (sigmaMax * std::erf(tmax / sqrt2) - sigmaMin * std::erf(tmin / sqrt2));
182}
183
184inline double integrateTailLogVersion(double sigma, double alpha, double n, double tmin, double tmax)
185{
186 double a = std::pow(n / alpha, n) * exp(-0.5 * alpha * alpha);
187 double b = n / alpha - alpha;
188
189 return a * sigma * (log(b - tmin) - log(b - tmax));
190}
191
192inline double integrateTailRegular(double sigma, double alpha, double n, double tmin, double tmax)
193{
194 double a = std::pow(n / alpha, n) * exp(-0.5 * alpha * alpha);
195 double b = n / alpha - alpha;
196
197 return a * sigma / (1.0 - n) * (1.0 / (std::pow(b - tmin, n - 1.0)) - 1.0 / (std::pow(b - tmax, n - 1.0)));
198}
199
200} // namespace
201
202////////////////////////////////////////////////////////////////////////////////
203
205{
206 const double x = x_;
207 const double x0 = x0_;
208 const double sigmaL = std::abs(sigmaL_);
209 const double sigmaR = std::abs(sigmaR_);
210 double alphaL = std::abs(alphaL_);
211 double nL = nL_;
212 double alphaR = alphaR_ ? std::abs(*alphaR_) : std::numeric_limits<double>::infinity();
213 double nR = nR_ ? *nR_ : 0.0;
214
215 // If alphaL is negative, then the tail will be on the right side.
216 // Like this, we follow the convention established by RooCBShape.
217 if(!alphaR_ && alphaL_ < 0.0) {
218 std::swap(alphaL, alphaR);
219 std::swap(nL, nR);
220 }
221
222 const double t = (x - x0) / (x < x0 ? sigmaL : sigmaR);
223
224 if (t < -alphaL) {
225 return evaluateCrystalBallTail(t, alphaL, nL);
226 } else if (t <= alphaR) {
227 return std::exp(-0.5 * t * t);
228 } else {
229 return evaluateCrystalBallTail(-t, alphaR, nR);
230 }
231}
232
233////////////////////////////////////////////////////////////////////////////////
234
235Int_t RooCrystalBall::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
236{
237 return matchArgs(allVars, analVars, x_) ? 1 : 0;
238}
239
240////////////////////////////////////////////////////////////////////////////////
241
242double RooCrystalBall::analyticalIntegral(Int_t code, const char *rangeName) const
243{
244 R__ASSERT(code == 1);
245
246 const double x0 = x0_;
247 const double sigmaL = std::abs(sigmaL_);
248 const double sigmaR = std::abs(sigmaR_);
249 double alphaL = std::abs(alphaL_);
250 double nL = nL_;
251 double alphaR = alphaR_ ? std::abs(*alphaR_) : std::numeric_limits<double>::infinity();
252 double nR = nR_ ? *nR_ : 0.0;
253
254 // If alphaL is negative, then the tail will be on the right side.
255 // Like this, we follow the convention established by RooCBShape.
256 if(!alphaR_ && alphaL_ < 0.0) {
257 std::swap(alphaL, alphaR);
258 std::swap(nL, nR);
259 }
260
261 constexpr double switchToLogThreshold = 1.0e-05;
262
263 const double xmin = x_.min(rangeName);
264 const double xmax = x_.max(rangeName);
265 const double tmin = (xmin - x0) / (xmin < x0 ? sigmaL : sigmaR);
266 const double tmax = (xmax - x0) / (xmax < x0 ? sigmaL : sigmaR);
267
268 double result = 0.0;
269
270 if (tmin < -alphaL) {
271 auto integrateTailL = std::abs(nL - 1.0) < switchToLogThreshold ? integrateTailLogVersion : integrateTailRegular;
272 result += integrateTailL(sigmaL, alphaL, nL, tmin, std::min(tmax, -alphaL));
273 }
274 if (tmax > alphaR) {
275 auto integrateTailR = std::abs(nR - 1.0) < switchToLogThreshold ? integrateTailLogVersion : integrateTailRegular;
276 result += integrateTailR(sigmaR, alphaR, nR, -tmax, std::min(-tmin, -alphaR));
277 }
278 if (tmin < alphaR && tmax > -alphaL) {
279 result += integrateGaussian(sigmaL, sigmaR, std::max(tmin, -alphaL), std::min(tmax, alphaR));
280 }
281
282 return result;
283}
284
285////////////////////////////////////////////////////////////////////////////////
286/// Advertise that we know the maximum of self for given (m0,alpha,n,sigma).
287
289{
290 RooArgSet dummy;
291 return matchArgs(vars, dummy, x_) ? 1 : 0;
292}
293
294////////////////////////////////////////////////////////////////////////////////
295
297{
298 R__ASSERT(code == 1);
299
300 // The maximum value for given (m0,alpha,n,sigma) is 1./ Integral in the variable range
301 // For the crystal ball, the maximum is 1.0 in the current implementation,
302 // but it's maybe better to keep this general in case the implementation changes.
303 return 1.0 / analyticalIntegral(code);
304}
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
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.
Definition: RooArgSet.h:57
PDF implementing the generalized Asymmetrical Double-Sided Crystall Ball line shape.
std::unique_ptr< RooRealProxy > nR_
RooRealProxy sigmaR_
RooRealProxy nL_
RooRealProxy x_
std::unique_ptr< RooRealProxy > alphaR_
RooRealProxy x0_
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooRealProxy alphaL_
double maxVal(Int_t code) const override
Return maximum value for set of observables identified by code assigned in getMaxVal.
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooRealProxy sigmaL_
Int_t getMaxVal(const RooArgSet &vars) const override
Advertise that we know the maximum of self for given (m0,alpha,n,sigma).
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
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.
double erf(double x)
Error function encountered in integrating the normal distribution.
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1739
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
void swap(RDirectoryEntry &e1, RDirectoryEntry &e2) noexcept
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
auto * a
Definition: textangle.C:12