ROOT   Reference Guide
RooCBShape.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
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-2005, Regents of the University of California *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17/** \class RooCBShape
18 \ingroup Roofit
19
20PDF implementing the Crystal Ball line shape.
21**/
22
23#include "RooCBShape.h"
24
25#include "RooAbsReal.h"
26#include "RooRealVar.h"
27#include "RooMath.h"
28#include "BatchHelpers.h"
30
31#include "TMath.h"
32
33#include <cmath>
34
35using namespace std;
36
38
39////////////////////////////////////////////////////////////////////////////////
40
42{
43 static const double erflim = 5.0;
44 if( arg > erflim )
45 return 1.0;
46 if( arg < -erflim )
47 return -1.0;
48
49 return RooMath::erf(arg);
50}
51
52////////////////////////////////////////////////////////////////////////////////
53
54RooCBShape::RooCBShape(const char *name, const char *title,
55 RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma,
56 RooAbsReal& _alpha, RooAbsReal& _n) :
57 RooAbsPdf(name, title),
58 m("m", "Dependent", this, _m),
59 m0("m0", "M0", this, _m0),
60 sigma("sigma", "Sigma", this, _sigma),
61 alpha("alpha", "Alpha", this, _alpha),
62 n("n", "Order", this, _n)
63{
64}
65
66////////////////////////////////////////////////////////////////////////////////
67
68RooCBShape::RooCBShape(const RooCBShape& other, const char* name) :
69 RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0),
70 sigma("sigma", this, other.sigma), alpha("alpha", this, other.alpha),
71 n("n", this, other.n)
72{
73}
74
75////////////////////////////////////////////////////////////////////////////////
76
78 Double_t t = (m-m0)/sigma;
79 if (alpha < 0) t = -t;
80
81 Double_t absAlpha = fabs((Double_t)alpha);
82
83 if (t >= -absAlpha) {
84 return exp(-0.5*t*t);
85 }
86 else {
87 Double_t a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
88 Double_t b= n/absAlpha - absAlpha;
89
90 return a/TMath::Power(b - t, n);
91 }
92}
93
94////////////////////////////////////////////////////////////////////////////////
95
96namespace {
97//Author: Emmanouil Michalainas, CERN 21 August 2019
98
99template<class Tm, class Tm0, class Tsigma, class Talpha, class Tn>
100void compute( size_t batchSize,
101 double * __restrict output,
102 Tm M, Tm0 M0, Tsigma S, Talpha A, Tn N)
103{
104 for (size_t i=0; i<batchSize; i++) {
105 const double t = (M[i]-M0[i]) / S[i];
106 if ((A[i]>0 && t>=-A[i]) || (A[i]<0 && -t>=A[i])) {
107 output[i] = -0.5*t*t;
108 } else {
109 output[i] = N[i] / (N[i] -A[i]*A[i] -A[i]*t);
110 output[i] = _rf_fast_log(output[i]);
111 output[i] *= N[i];
112 output[i] -= 0.5*A[i]*A[i];
113 }
114 }
115
116 for (size_t i=0; i<batchSize; i++) {
117 output[i] = _rf_fast_exp(output[i]);
118 }
119}
120};
121
122RooSpan<double> RooCBShape::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
123 using namespace BatchHelpers;
124
125 EvaluateInfo info = getInfo( {&m, &m0, &sigma, &alpha, &n}, begin, batchSize );
126 if (info.nBatches == 0) {
127 return {};
128 }
129 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
130 auto mData = m.getValBatch(begin, info.size);
131
132 if (info.nBatches==1 && !mData.empty()) {
133 compute(info.size, output.data(), mData.data(),
138 }
139 else {
140 compute(info.size, output.data(),
146 }
147 return output;
148}
149
150////////////////////////////////////////////////////////////////////////////////
151
152Int_t RooCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
153{
154 if( matchArgs(allVars,analVars,m) )
155 return 1 ;
156
157 return 0;
158}
159
160////////////////////////////////////////////////////////////////////////////////
161
162Double_t RooCBShape::analyticalIntegral(Int_t code, const char* rangeName) const
163{
164 static const double sqrtPiOver2 = 1.2533141373;
165 static const double sqrt2 = 1.4142135624;
166
167 R__ASSERT(code==1);
168 double result = 0.0;
169 bool useLog = false;
170
171 if( fabs(n-1.0) < 1.0e-05 )
172 useLog = true;
173
174 double sig = fabs((Double_t)sigma);
175
176 double tmin = (m.min(rangeName)-m0)/sig;
177 double tmax = (m.max(rangeName)-m0)/sig;
178
179 if(alpha < 0) {
180 double tmp = tmin;
181 tmin = -tmax;
182 tmax = -tmp;
183 }
184
185 double absAlpha = fabs((Double_t)alpha);
186
187 if( tmin >= -absAlpha ) {
188 result += sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
189 - ApproxErf(tmin/sqrt2) );
190 }
191 else if( tmax <= -absAlpha ) {
192 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
193 double b = n/absAlpha - absAlpha;
194
195 if(useLog) {
196 result += a*sig*( log(b-tmin) - log(b-tmax) );
197 }
198 else {
199 result += a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
200 - 1.0/(TMath::Power(b-tmax,n-1.0)) );
201 }
202 }
203 else {
204 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
205 double b = n/absAlpha - absAlpha;
206
207 double term1 = 0.0;
208 if(useLog) {
209 term1 = a*sig*( log(b-tmin) - log(n/absAlpha));
210 }
211 else {
212 term1 = a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
213 - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
214 }
215
216 double term2 = sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
217 - ApproxErf(-absAlpha/sqrt2) );
218
219
220 result += term1 + term2;
221 }
222
223 return result != 0. ? result : 1.E-300;
224}
225
226////////////////////////////////////////////////////////////////////////////////
227/// Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
228
230{
232
233 if (matchArgs(vars,dummy,m)) {
234 return 1 ;
235 }
236 return 0 ;
237}
238
239////////////////////////////////////////////////////////////////////////////////
240
242{
243 R__ASSERT(code==1) ;
244
245 // The maximum value for given (m0,alpha,n,sigma)
246 return 1.0 ;
247}
#define b(i)
Definition: RSha256.hxx:100
static RooMathCoreReg dummy
double _rf_fast_exp(double x)
double _rf_fast_log(double x)
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
#define R__ASSERT(e)
Definition: TError.h:96
#define N
char name[80]
Definition: TGX11.cxx:109
double exp(double)
double log(double)
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
Little adapter that gives a bracket operator to types that don't have one.
Definition: BatchHelpers.h:58
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
PDF implementing the Crystal Ball line shape.
Definition: RooCBShape.h:24
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
Definition: RooCBShape.cxx:241
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooCBShape.cxx:162
RooRealProxy n
Definition: RooCBShape.h:51
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooCBShape.cxx:77
Double_t ApproxErf(Double_t arg) const
Definition: RooCBShape.cxx:41
RooRealProxy m0
Definition: RooCBShape.h:48
RooRealProxy m
Definition: RooCBShape.h:47
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooCBShape.cxx:152
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
Definition: RooCBShape.cxx:229
RooRealProxy sigma
Definition: RooCBShape.h:49
RooRealProxy alpha
Definition: RooCBShape.h:50
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
Definition: RooCBShape.cxx:122
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:580
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
const Double_t sigma
const Int_t n
Definition: legend1.C:16
EvaluateInfo getInfo(std::vector< const RooRealProxy * > parameters, size_t begin, size_t batchSize)
static double A[]
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
RooArgSet S(const RooAbsArg &v1)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
static void output(int code)
Definition: gifencode.c:226