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 "RooRealVar.h"
26#include "RooMath.h"
27#include "RooBatchCompute.h"
28
29#include "TMath.h"
30
31#include <cmath>
32
33using namespace std;
34
36
37////////////////////////////////////////////////////////////////////////////////
38
39double RooCBShape::ApproxErf(double arg) const
40{
41 static const double erflim = 5.0;
42 if( arg > erflim )
43 return 1.0;
44 if( arg < -erflim )
45 return -1.0;
46
47 return RooMath::erf(arg);
48}
49
50////////////////////////////////////////////////////////////////////////////////
51
52RooCBShape::RooCBShape(const char *name, const char *title,
53 RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma,
54 RooAbsReal& _alpha, RooAbsReal& _n) :
55 RooAbsPdf(name, title),
56 m("m", "Dependent", this, _m),
57 m0("m0", "M0", this, _m0),
58 sigma("sigma", "Sigma", this, _sigma),
59 alpha("alpha", "Alpha", this, _alpha),
60 n("n", "Order", this, _n)
61{
62}
63
64////////////////////////////////////////////////////////////////////////////////
65
66RooCBShape::RooCBShape(const RooCBShape& other, const char* name) :
67 RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0),
68 sigma("sigma", this, other.sigma), alpha("alpha", this, other.alpha),
69 n("n", this, other.n)
70{
71}
72
73////////////////////////////////////////////////////////////////////////////////
74
75double RooCBShape::evaluate() const {
76 double t = (m-m0)/sigma;
77 if (alpha < 0) t = -t;
78
79 double absAlpha = fabs((double)alpha);
80
81 if (t >= -absAlpha) {
82 return exp(-0.5*t*t);
83 }
84 else {
85 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
86 double b= n/absAlpha - absAlpha;
87
88 return a/TMath::Power(b - t, n);
89 }
90}
91
92////////////////////////////////////////////////////////////////////////////////
93/// Compute multiple values of Crystal ball Shape distribution.
94void RooCBShape::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
95{
97 dispatch->compute(stream, RooBatchCompute::CBShape, output, nEvents,
98 {dataMap.at(m), dataMap.at(m0), dataMap.at(sigma), dataMap.at(alpha), dataMap.at(n)});
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
103Int_t RooCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
104{
105 if( matchArgs(allVars,analVars,m) )
106 return 1 ;
107
108 return 0;
109}
110
111////////////////////////////////////////////////////////////////////////////////
112
113double RooCBShape::analyticalIntegral(Int_t code, const char* rangeName) const
114{
115 static const double sqrtPiOver2 = 1.2533141373;
116 static const double sqrt2 = 1.4142135624;
117
118 R__ASSERT(code==1);
119 double result = 0.0;
120 bool useLog = false;
121
122 if( fabs(n-1.0) < 1.0e-05 )
123 useLog = true;
124
125 double sig = fabs((double)sigma);
126
127 double tmin = (m.min(rangeName)-m0)/sig;
128 double tmax = (m.max(rangeName)-m0)/sig;
129
130 if(alpha < 0) {
131 double tmp = tmin;
132 tmin = -tmax;
133 tmax = -tmp;
134 }
135
136 double absAlpha = fabs((double)alpha);
137
138 if( tmin >= -absAlpha ) {
139 result += sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
140 - ApproxErf(tmin/sqrt2) );
141 }
142 else if( tmax <= -absAlpha ) {
143 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
144 double b = n/absAlpha - absAlpha;
145
146 if(useLog) {
147 result += a*sig*( log(b-tmin) - log(b-tmax) );
148 }
149 else {
150 result += a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
151 - 1.0/(TMath::Power(b-tmax,n-1.0)) );
152 }
153 }
154 else {
155 double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
156 double b = n/absAlpha - absAlpha;
157
158 double term1 = 0.0;
159 if(useLog) {
160 term1 = a*sig*( log(b-tmin) - log(n/absAlpha));
161 }
162 else {
163 term1 = a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
164 - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
165 }
166
167 double term2 = sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
168 - ApproxErf(-absAlpha/sqrt2) );
169
170
171 result += term1 + term2;
172 }
173
174 return result != 0. ? result : 1.E-300;
175}
176
177////////////////////////////////////////////////////////////////////////////////
178/// Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
179
181{
182 RooArgSet dummy ;
183
184 if (matchArgs(vars,dummy,m)) {
185 return 1 ;
186 }
187 return 0 ;
188}
189
190////////////////////////////////////////////////////////////////////////////////
191
192double RooCBShape::maxVal(Int_t code) const
193{
194 R__ASSERT(code==1) ;
195
196 // The maximum value for given (m0,alpha,n,sigma)
197 // is 1./ Integral in the variable range
198 return 1.0/analyticalIntegral(1) ;
199}
#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:62
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:56
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:39
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooCBShape.cxx:103
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:94
RooRealProxy n
Definition: RooCBShape.h:51
RooRealProxy m0
Definition: RooCBShape.h:48
RooRealProxy m
Definition: RooCBShape.h:47
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooCBShape.cxx:113
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:75
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:180
double maxVal(Int_t code) const override
Return maximum value for set of observables identified by code assigned in getMaxVal.
Definition: RooCBShape.cxx:192
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 max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1787
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1783
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:719
TMarker m
Definition: textangle.C:8
TArc a
Definition: textangle.C:12
static void output()