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
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
77 Double_t t = (m-m0)/sigma;
78 if (alpha < 0) t = -t;
79
80 Double_t absAlpha = fabs((Double_t)alpha);
81
82 if (t >= -absAlpha) {
83 return exp(-0.5*t*t);
84 }
85 else {
86 Double_t a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
87 Double_t 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, RooBatchCompute::DataMap& dataMap) const
96{
98 dispatch->compute(stream, RooBatchCompute::CBShape, output, nEvents, dataMap, {&*m,&*m0,&*sigma,&*alpha,&*n,&*_norm});
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_t 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_t)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_t)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
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 b(i)
Definition: RSha256.hxx:100
#define ClassImp(name)
Definition: Rtypes.h:364
#define R__ASSERT(e)
Definition: TError.h:118
char name[80]
Definition: TGX11.cxx:110
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
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.
Definition: RooArgSet.h:35
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const DataMap &, const VarVector &, const ArgVector &={})=0
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:192
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooCBShape.cxx:113
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:76
Double_t ApproxErf(Double_t arg) const
Definition: RooCBShape.cxx:40
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:103
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:180
RooRealProxy sigma
Definition: RooCBShape.h:49
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const
Compute multiple values of Crystal ball Shape distribution.
Definition: RooCBShape.cxx:95
RooRealProxy alpha
Definition: RooCBShape.h:50
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:575
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.
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
std::map< DataKey, RooSpan< const double > > DataMap
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)
Definition: TMath.h:685
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
static void output(int code)
Definition: gifencode.c:226