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 "RooFit.h"
30#include "RooRealVar.h"
31#include "RooBatchCompute.h"
32
33#include "TMath.h"
34
35#include <cmath>
36using namespace std;
37
39
40////////////////////////////////////////////////////////////////////////////////
41
42RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
43 RooAbsReal& _x, RooAbsReal& _peak,
44 RooAbsReal& _width, RooAbsReal& _tail) :
45 // The two addresses refer to our first dependent variable and
46 // parameter, respectively, as declared in the rdl file
47 RooAbsPdf(name, title),
48 x("x","x",this,_x),
49 width("width","width",this,_width),
50 peak("peak","peak",this,_peak),
51 tail("tail","tail",this,_tail)
52{
53}
54
55////////////////////////////////////////////////////////////////////////////////
56
58 RooAbsPdf(other,name),
59 x("x",this,other.x),
60 width("width",this,other.width),
61 peak("peak",this,other.peak),
62 tail("tail",this,other.tail)
63{
64}
65
66////////////////////////////////////////////////////////////////////////////////
67///If tail=eta=0 the Belle distribution becomes gaussian
68
70{
71 if (TMath::Abs(tail) < 1.e-7) {
72 return TMath::Exp( -0.5 * TMath::Power( ( (x - peak) / width ), 2 ));
73 }
74
75 Double_t arg = 1.0 - ( x - peak ) * tail / width;
76
77 if (arg < 1.e-7) {
78 //Argument of logarithm negative. Real continuation -> function equals zero
79 return 0.0;
80 }
81
82 Double_t log = TMath::Log(arg);
83 static const Double_t xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
84
85 Double_t width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
86 Double_t width_zero2 = width_zero * width_zero;
87 Double_t exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
88
89 return TMath::Exp(exponent) ;
90}
91////////////////////////////////////////////////////////////////////////////////
92/// Compute multiple values of Novosibirsk distribution.
93void RooNovosibirsk::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
94{
96 dispatch->compute(stream, RooBatchCompute::Novosibirsk, output, nEvents, dataMap, {&*x,&*peak,&*width,&*tail,&*_norm});
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_t 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_t sqrt2 = 1.4142135623730950; // Sqrt(2)
121 static const Double_t sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
122 static const Double_t sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
123 static const Double_t log4 = 1.38629436111989062; //Log(2)
124 static const Double_t rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
125 static const Double_t sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
126
127 if (code==1) {
128 Double_t A = x.min(rangeName);
129 Double_t B = x.max(rangeName);
130
131 Double_t 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_t 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_t log_argument_A = ( (peak - A)*tail + width ) / width ;
147 Double_t 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_t term1 = TMath::ASinH( tail * sqlog4 );
160 Double_t term1_2 = term1 * term1;
161
162 //Calculate the error function arguments
163 Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
164 Double_t 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_t A = x.min(rangeName);
172 Double_t B = x.max(rangeName);
173
174 Double_t 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_t 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_t log_argument_A = ( (A - x)*tail + width ) / width;
190 Double_t 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_t term1 = TMath::ASinH( tail * sqlog4 );
203 Double_t term1_2 = term1 * term1;
204
205 //Calculate the error function arguments
206 Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
207 Double_t 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:34
int Int_t
Definition: RtypesCore.h:45
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
include TDocParser_001 C image html pict1_TDocParser_001 png width
Definition: TDocParser.cxx:121
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
RooNovosibirsk implements the Novosibirsk function.
RooRealProxy width
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Double_t evaluate() const
If tail=eta=0 the Belle distribution becomes gaussian.
RooRealProxy peak
RooRealProxy tail
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const
Compute multiple values of Novosibirsk distribution.
RooRealProxy x
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.
Double_t x[n]
Definition: legend1.C:17
static double B[]
static double A[]
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...
Double_t Exp(Double_t x)
Definition: TMath.h:677
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:184
Double_t ASinH(Double_t)
Definition: TMath.cxx:64
Double_t Log(Double_t x)
Definition: TMath.h:710
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:685
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226