Logo 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 *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
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 "RooBatchCompute.h"
29
30#include "TMath.h"
31
32#include <cmath>
33
34using namespace std;
35
37
38////////////////////////////////////////////////////////////////////////////////
39
40double RooCBShape::ApproxErf(double arg) const
41{
42 static const double erflim = 5.0;
43 if( arg > erflim )
44 return 1.0;
45 if( arg < -erflim )
46 return -1.0;
47
48 return RooMath::erf(arg);
49}
50
51////////////////////////////////////////////////////////////////////////////////
52
53RooCBShape::RooCBShape(const char *name, const char *title,
54 RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma,
55 RooAbsReal& _alpha, RooAbsReal& _n) :
56 RooAbsPdf(name, title),
57 m("m", "Dependent", this, _m),
58 m0("m0", "M0", this, _m0),
59 sigma("sigma", "Sigma", this, _sigma),
60 alpha("alpha", "Alpha", this, _alpha),
61 n("n", "Order", this, _n)
62{
63}
64
65////////////////////////////////////////////////////////////////////////////////
66
67RooCBShape::RooCBShape(const RooCBShape& other, const char* name) :
68 RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0),
69 sigma("sigma", this, other.sigma), alpha("alpha", this, other.alpha),
70 n("n", this, other.n)
71{
72}
73
74////////////////////////////////////////////////////////////////////////////////
75
76double RooCBShape::evaluate() const {
77 double t = (m-m0)/sigma;
78 if (alpha < 0) t = -t;
79
80 double absAlpha = fabs((double)alpha);
81
82 if (t >= -absAlpha) {
83 return exp(-0.5*t*t);
84 }
85 else {
86 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
87 double b= n/absAlpha - absAlpha;
88
89 return a/TMath::Power(b - t, n);
90 }
91}
92
93////////////////////////////////////////////////////////////////////////////////
94/// Compute multiple values of Crystal ball Shape distribution.
95void RooCBShape::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
96{
98 dispatch->compute(stream, RooBatchCompute::CBShape, output, nEvents,
99 {dataMap.at(m), dataMap.at(m0), dataMap.at(sigma), dataMap.at(alpha), dataMap.at(n)});
100}
101
102////////////////////////////////////////////////////////////////////////////////
103
104Int_t RooCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
105{
106 if( matchArgs(allVars,analVars,m) )
107 return 1 ;
108
109 return 0;
110}
111
112////////////////////////////////////////////////////////////////////////////////
113
114double RooCBShape::analyticalIntegral(Int_t code, const char* rangeName) const
115{
116 static const double sqrtPiOver2 = 1.2533141373;
117 static const double sqrt2 = 1.4142135624;
118
119 R__ASSERT(code==1);
120 double result = 0.0;
121 bool useLog = false;
122
123 if( fabs(n-1.0) < 1.0e-05 )
124 useLog = true;
125
126 double sig = fabs((double)sigma);
127
128 double tmin = (m.min(rangeName)-m0)/sig;
129 double tmax = (m.max(rangeName)-m0)/sig;
130
131 if(alpha < 0) {
132 double tmp = tmin;
133 tmin = -tmax;
134 tmax = -tmp;
135 }
136
137 double absAlpha = fabs((double)alpha);
138
139 if( tmin >= -absAlpha ) {
140 result += sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
141 - ApproxErf(tmin/sqrt2) );
142 }
143 else if( tmax <= -absAlpha ) {
144 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
145 double b = n/absAlpha - absAlpha;
146
147 if(useLog) {
148 result += a*sig*( log(b-tmin) - log(b-tmax) );
149 }
150 else {
151 result += a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
152 - 1.0/(TMath::Power(b-tmax,n-1.0)) );
153 }
154 }
155 else {
156 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
157 double b = n/absAlpha - absAlpha;
158
159 double term1 = 0.0;
160 if(useLog) {
161 term1 = a*sig*( log(b-tmin) - log(n/absAlpha));
162 }
163 else {
164 term1 = a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
165 - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
166 }
167
168 double term2 = sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
169 - ApproxErf(-absAlpha/sqrt2) );
170
171
172 result += term1 + term2;
173 }
174
175 return result != 0. ? result : 1.E-300;
176}
177
178////////////////////////////////////////////////////////////////////////////////
179/// Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
180
182{
183 RooArgSet dummy ;
184
185 if (matchArgs(vars,dummy,m)) {
186 return 1 ;
187 }
188 return 0 ;
189}
190
191////////////////////////////////////////////////////////////////////////////////
192
193double RooCBShape::maxVal(Int_t code) const
194{
195 R__ASSERT(code==1) ;
196
197 // The maximum value for given (m0,alpha,n,sigma)
198 // is 1./ Integral in the variable range
199 return 1.0/analyticalIntegral(1) ;
200}
#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
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
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
PDF implementing the Crystal Ball line shape.
Definition: RooCBShape.h:24
double ApproxErf(double arg) const
Definition: RooCBShape.cxx:40
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooCBShape.cxx:114
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of Crystal ball Shape distribution.
Definition: RooCBShape.cxx:95
RooRealProxy n
Definition: RooCBShape.h:51
RooRealProxy m0
Definition: RooCBShape.h:48
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooCBShape.cxx:104
RooRealProxy m
Definition: RooCBShape.h:47
RooRealProxy sigma
Definition: RooCBShape.h:49
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooCBShape.cxx:76
RooRealProxy alpha
Definition: RooCBShape.h:50
Int_t getMaxVal(const RooArgSet &vars) const override
Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
Definition: RooCBShape.cxx:181
double maxVal(Int_t code) const override
Return maximum value for set of observables identified by code assigned in getMaxVal.
Definition: RooCBShape.cxx:193
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:60
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.
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
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
R__EXTERN RooBatchComputeInterface * dispatchCUDA
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:718
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
static void output()