Logo ROOT  
Reference Guide
RooNovosibirsk.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * DB, Dieter Best, UC Irvine, best@slac.stanford.edu *
7 * HT, Hirohisa Tanaka SLAC tanaka@slac.stanford.edu *
8 * *
9 * Updated version with analytical integral *
10 * MP, Marko Petric, J. Stefan Institute, marko.petric@ijs.si *
11 * *
12 * Copyright (c) 2000-2013, Regents of the University of California *
13 * and Stanford University. All rights reserved. *
14 * *
15 * Redistribution and use in source and binary forms, *
16 * with or without modification, are permitted according to the terms *
17 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
18 *****************************************************************************/
19
20/** \class RooNovosibirsk
21 \ingroup Roofit
22
23RooNovosibirsk implements the Novosibirsk function
24
25Function taken from H. Ikeda et al. NIM A441 (2000), p. 401 (Belle Collaboration)
26
27**/
28#include "RooNovosibirsk.h"
29#include "RooRealVar.h"
30#include "RooBatchCompute.h"
31
32#include "TMath.h"
33
34#include <cmath>
35using namespace std;
36
38
39////////////////////////////////////////////////////////////////////////////////
40
41RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
42 RooAbsReal& _x, RooAbsReal& _peak,
43 RooAbsReal& _width, RooAbsReal& _tail) :
44 // The two addresses refer to our first dependent variable and
45 // parameter, respectively, as declared in the rdl file
46 RooAbsPdf(name, title),
47 x("x","x",this,_x),
48 width("width","width",this,_width),
49 peak("peak","peak",this,_peak),
50 tail("tail","tail",this,_tail)
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55
57 RooAbsPdf(other,name),
58 x("x",this,other.x),
59 width("width",this,other.width),
60 peak("peak",this,other.peak),
61 tail("tail",this,other.tail)
62{
63}
64
65////////////////////////////////////////////////////////////////////////////////
66///If tail=eta=0 the Belle distribution becomes gaussian
67
69{
70 if (TMath::Abs(tail) < 1.e-7) {
71 return TMath::Exp( -0.5 * TMath::Power( ( (x - peak) / width ), 2 ));
72 }
73
74 double arg = 1.0 - ( x - peak ) * tail / width;
75
76 if (arg < 1.e-7) {
77 //Argument of logarithm negative. Real continuation -> function equals zero
78 return 0.0;
79 }
80
81 double log = TMath::Log(arg);
82 static const double xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
83
84 double width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
85 double width_zero2 = width_zero * width_zero;
86 double exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
87
88 return TMath::Exp(exponent) ;
89}
90////////////////////////////////////////////////////////////////////////////////
91/// Compute multiple values of Novosibirsk distribution.
92void RooNovosibirsk::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
93{
95 dispatch->compute(stream, RooBatchCompute::Novosibirsk, output, nEvents,
96 {dataMap.at(x), dataMap.at(peak), dataMap.at(width), dataMap.at(tail)});
97}
98
99////////////////////////////////////////////////////////////////////////////////
100
101Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
102{
103 if (matchArgs(allVars,analVars,x)) return 1 ;
104 if (matchArgs(allVars,analVars,peak)) return 2 ;
105
106 //The other two integrals over tali and width are not integrable
107
108 return 0 ;
109}
110
111////////////////////////////////////////////////////////////////////////////////
112
113double RooNovosibirsk::analyticalIntegral(Int_t code, const char* rangeName) const
114{
115 assert(code==1 || code==2) ;
116
117 //The range is defined as [A,B]
118
119 //Numerical values need for the evaluation of the integral
120 static const double sqrt2 = 1.4142135623730950; // Sqrt(2)
121 static const double sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
122 static const double sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
123 static const double log4 = 1.38629436111989062; //Log(2)
124 static const double rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
125 static const double sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
126
127 if (code==1) {
128 double A = x.min(rangeName);
129 double B = x.max(rangeName);
130
131 double result = 0;
132
133
134 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
135 if (TMath::Abs(tail) < 1.e-7) {
136
137 double xscale = sqrt2*width;
138
139 result = rootpiby2*width*(TMath::Erf((B-peak)/xscale)-TMath::Erf((A-peak)/xscale));
140
141 return result;
142
143 }
144
145 // If the range is not defined correctly the function becomes complex
146 double log_argument_A = ( (peak - A)*tail + width ) / width ;
147 double log_argument_B = ( (peak - B)*tail + width ) / width ;
148
149 //lower limit
150 if ( log_argument_A < 1.e-7) {
151 log_argument_A = 1.e-7;
152 }
153
154 //upper limit
155 if ( log_argument_B < 1.e-7) {
156 log_argument_B = 1.e-7;
157 }
158
159 double term1 = TMath::ASinH( tail * sqlog4 );
160 double term1_2 = term1 * term1;
161
162 //Calculate the error function arguments
163 double erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
164 double erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
165
166 result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
167
168 return result;
169
170 } else if (code==2) {
171 double A = x.min(rangeName);
172 double B = x.max(rangeName);
173
174 double result = 0;
175
176
177 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
178 if (TMath::Abs(tail) < 1.e-7) {
179
180 double xscale = sqrt2*width;
181
182 result = rootpiby2*width*(TMath::Erf((B-x)/xscale)-TMath::Erf((A-x)/xscale));
183
184 return result;
185
186 }
187
188 // If the range is not defined correctly the function becomes complex
189 double log_argument_A = ( (A - x)*tail + width ) / width;
190 double log_argument_B = ( (B - x)*tail + width ) / width;
191
192 //lower limit
193 if ( log_argument_A < 1.e-7) {
194 log_argument_A = 1.e-7;
195 }
196
197 //upper limit
198 if ( log_argument_B < 1.e-7) {
199 log_argument_B = 1.e-7;
200 }
201
202 double term1 = TMath::ASinH( tail * sqlog4 );
203 double term1_2 = term1 * term1;
204
205 //Calculate the error function arguments
206 double erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
207 double erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
208
209 result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
210
211 return result;
212
213 }
214
215 // Emit fatal error
216 coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;
217
218 // Put dummy return here to avoid compiler warnings
219 return 1.0 ;
220}
#define coutF(a)
Definition: RooMsgService.h:38
#define ClassImp(name)
Definition: Rtypes.h:375
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
Option_t Option_t width
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
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
RooNovosibirsk implements the Novosibirsk function.
RooRealProxy width
RooRealProxy peak
RooRealProxy tail
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double evaluate() const override
If tail=eta=0 the Belle distribution becomes gaussian.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of Novosibirsk distribution.
RooRealProxy x
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
Double_t x[n]
Definition: legend1.C:17
static double B[]
static double A[]
R__EXTERN RooBatchComputeInterface * dispatchCUDA
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Definition: TMath.h:707
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:190
Double_t ASinH(Double_t)
Returns the area hyperbolic sine of x.
Definition: TMath.cxx:67
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition: TMath.h:754
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
static void output()