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 "BatchHelpers.h"
32#include "RooVDTHeaders.h"
33
34#include "TMath.h"
35
36#include <cmath>
37using namespace std;
38
40
41////////////////////////////////////////////////////////////////////////////////
42
43RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
44 RooAbsReal& _x, RooAbsReal& _peak,
45 RooAbsReal& _width, RooAbsReal& _tail) :
46 // The two addresses refer to our first dependent variable and
47 // parameter, respectively, as declared in the rdl file
48 RooAbsPdf(name, title),
49 x("x","x",this,_x),
50 width("width","width",this,_width),
51 peak("peak","peak",this,_peak),
52 tail("tail","tail",this,_tail)
53{
54}
55
56////////////////////////////////////////////////////////////////////////////////
57
59 RooAbsPdf(other,name),
60 x("x",this,other.x),
61 width("width",this,other.width),
62 peak("peak",this,other.peak),
63 tail("tail",this,other.tail)
64{
65}
66
67////////////////////////////////////////////////////////////////////////////////
68///If tail=eta=0 the Belle distribution becomes gaussian
69
71{
72 if (TMath::Abs(tail) < 1.e-7) {
73 return TMath::Exp( -0.5 * TMath::Power( ( (x - peak) / width ), 2 ));
74 }
75
76 Double_t arg = 1.0 - ( x - peak ) * tail / width;
77
78 if (arg < 1.e-7) {
79 //Argument of logarithm negative. Real continuation -> function equals zero
80 return 0.0;
81 }
82
83 Double_t log = TMath::Log(arg);
84 static const Double_t xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
85
86 Double_t width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
87 Double_t width_zero2 = width_zero * width_zero;
88 Double_t exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
89
90 return TMath::Exp(exponent) ;
91}
92
93////////////////////////////////////////////////////////////////////////////////
94
95namespace {
96//Author: Emmanouil Michalainas, CERN 10 September 2019
97
98/* TMath::ASinH(x) needs to be replaced with ln( x + sqrt(x^2+1))
99 * argasinh -> the argument of TMath::ASinH()
100 * argln -> the argument of the logarithm that replaces AsinH
101 * asinh -> the value that the function evaluates to
102 *
103 * ln is the logarithm that was solely present in the initial
104 * formula, that is before the asinh replacement
105 */
106template<class Tx, class Twidth, class Tpeak, class Ttail>
107void compute( size_t batchSize, double * __restrict output,
108 Tx X, Tpeak P, Twidth W, Ttail T)
109{
110 constexpr double xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
111 for (size_t i=0; i<batchSize; i++) {
112 double argasinh = 0.5*xi*T[i];
113 double argln = argasinh + 1/_rf_fast_isqrt(argasinh*argasinh +1);
114 double asinh = _rf_fast_log(argln);
115
116 double argln2 = 1 -(X[i]-P[i])*T[i]/W[i];
117 double ln = _rf_fast_log(argln2);
118 output[i] = ln/asinh;
119 output[i] *= -0.125*xi*xi*output[i];
120 output[i] -= 2.0/xi/xi*asinh*asinh;
121 }
122
123 //faster if you exponentiate in a seperate loop (dark magic!)
124 for (size_t i=0; i<batchSize; i++) {
125 output[i] = _rf_fast_exp(output[i]);
126 }
127}
128};
129
130RooSpan<double> RooNovosibirsk::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
131 using namespace BatchHelpers;
132
133 EvaluateInfo info = getInfo( {&x, &peak, &width, &tail}, begin, batchSize );
134 if (info.nBatches == 0) {
135 return {};
136 }
137 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
138 auto xData = x.getValBatch(begin, info.size);
139
140 if (info.nBatches==1 && !xData.empty()) {
141 compute(info.size, output.data(), xData.data(),
145 }
146 else {
147 compute(info.size, output.data(),
152 }
153 return output;
154}
155
156////////////////////////////////////////////////////////////////////////////////
157
158Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
159{
160 if (matchArgs(allVars,analVars,x)) return 1 ;
161 if (matchArgs(allVars,analVars,peak)) return 2 ;
162
163 //The other two integrals over tali and width are not integrable
164
165 return 0 ;
166}
167
168////////////////////////////////////////////////////////////////////////////////
169
170Double_t RooNovosibirsk::analyticalIntegral(Int_t code, const char* rangeName) const
171{
172 assert(code==1 || code==2) ;
173
174 //The range is defined as [A,B]
175
176 //Numerical values need for the evaluation of the integral
177 static const Double_t sqrt2 = 1.4142135623730950; // Sqrt(2)
178 static const Double_t sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
179 static const Double_t sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
180 static const Double_t log4 = 1.38629436111989062; //Log(2)
181 static const Double_t rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
182 static const Double_t sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
183
184 if (code==1) {
185 Double_t A = x.min(rangeName);
186 Double_t B = x.max(rangeName);
187
188 Double_t result = 0;
189
190
191 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
192 if (TMath::Abs(tail) < 1.e-7) {
193
194 Double_t xscale = sqrt2*width;
195
196 result = rootpiby2*width*(TMath::Erf((B-peak)/xscale)-TMath::Erf((A-peak)/xscale));
197
198 return result;
199
200 }
201
202 // If the range is not defined correctly the function becomes complex
203 Double_t log_argument_A = ( (peak - A)*tail + width ) / width ;
204 Double_t log_argument_B = ( (peak - B)*tail + width ) / width ;
205
206 //lower limit
207 if ( log_argument_A < 1.e-7) {
208 log_argument_A = 1.e-7;
209 }
210
211 //upper limit
212 if ( log_argument_B < 1.e-7) {
213 log_argument_B = 1.e-7;
214 }
215
216 Double_t term1 = TMath::ASinH( tail * sqlog4 );
217 Double_t term1_2 = term1 * term1;
218
219 //Calculate the error function arguments
220 Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
221 Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
222
223 result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
224
225 return result;
226
227 } else if (code==2) {
228 Double_t A = x.min(rangeName);
229 Double_t B = x.max(rangeName);
230
231 Double_t result = 0;
232
233
234 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
235 if (TMath::Abs(tail) < 1.e-7) {
236
237 Double_t xscale = sqrt2*width;
238
239 result = rootpiby2*width*(TMath::Erf((B-x)/xscale)-TMath::Erf((A-x)/xscale));
240
241 return result;
242
243 }
244
245 // If the range is not defined correctly the function becomes complex
246 Double_t log_argument_A = ( (A - x)*tail + width ) / width;
247 Double_t log_argument_B = ( (B - x)*tail + width ) / width;
248
249 //lower limit
250 if ( log_argument_A < 1.e-7) {
251 log_argument_A = 1.e-7;
252 }
253
254 //upper limit
255 if ( log_argument_B < 1.e-7) {
256 log_argument_B = 1.e-7;
257 }
258
259 Double_t term1 = TMath::ASinH( tail * sqlog4 );
260 Double_t term1_2 = term1 * term1;
261
262 //Calculate the error function arguments
263 Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
264 Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
265
266 result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
267
268 return result;
269
270 }
271
272 // Emit fatal error
273 coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;
274
275 // Put dummy return here to avoid compiler warnings
276 return 1.0 ;
277}
#define coutF(a)
Definition: RooMsgService.h:35
double _rf_fast_exp(double x)
VDT headers for RooFit.
Definition: RooVDTHeaders.h:47
double _rf_fast_log(double x)
Definition: RooVDTHeaders.h:51
double _rf_fast_isqrt(double x)
Definition: RooVDTHeaders.h:55
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
include TDocParser_001 C image html pict1_TDocParser_001 png width
Definition: TDocParser.cxx:121
char name[80]
Definition: TGX11.cxx:109
double log(double)
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
Little adapter that gives a bracket operator to types that don't have one.
Definition: BatchHelpers.h:58
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooNovosibirsk implements the Novosibirsk function.
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
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
RooRealProxy x
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t 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
EvaluateInfo getInfo(std::vector< const RooRealProxy * > parameters, size_t begin, size_t batchSize)
static double B[]
static double P[]
static double A[]
double T(double x)
Definition: ChebyshevPol.h:34
Double_t Exp(Double_t x)
Definition: TMath.h:717
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:750
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226