Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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>
35
37
38////////////////////////////////////////////////////////////////////////////////
39
40RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
41 RooAbsReal& _x, RooAbsReal& _peak,
42 RooAbsReal& _width, RooAbsReal& _tail) :
43 // The two addresses refer to our first dependent variable and
44 // parameter, respectively, as declared in the rdl file
45 RooAbsPdf(name, title),
46 x("x","x",this,_x),
47 width("width","width",this,_width),
48 peak("peak","peak",this,_peak),
49 tail("tail","tail",this,_tail)
50{
51}
52
53////////////////////////////////////////////////////////////////////////////////
54
56 RooAbsPdf(other,name),
57 x("x",this,other.x),
58 width("width",this,other.width),
59 peak("peak",this,other.peak),
60 tail("tail",this,other.tail)
61{
62}
63
64////////////////////////////////////////////////////////////////////////////////
65///If tail=eta=0 the Belle distribution becomes gaussian
66
68{
69 if (std::abs(tail) < 1.e-7) {
70 return std::exp( -0.5 * std::pow( ( (x - peak) / width ), 2 ));
71 }
72
73 double arg = 1.0 - ( x - peak ) * tail / width;
74
75 if (arg < 1.e-7) {
76 //Argument of logarithm negative. Real continuation -> function equals zero
77 return 0.0;
78 }
79
80 double log = std::log(arg);
81 static const double xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
82
83 double width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
84 double width_zero2 = width_zero * width_zero;
85 double exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
86
87 return std::exp(exponent) ;
88}
89////////////////////////////////////////////////////////////////////////////////
90/// Compute multiple values of Novosibirsk distribution.
92{
94 {ctx.at(x), ctx.at(peak), ctx.at(width), ctx.at(tail)});
95}
96
97////////////////////////////////////////////////////////////////////////////////
98
99Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
100{
101 if (matchArgs(allVars,analVars,x)) return 1 ;
102 if (matchArgs(allVars,analVars,peak)) return 2 ;
103
104 //The other two integrals over tali and width are not integrable
105
106 return 0 ;
107}
108
109////////////////////////////////////////////////////////////////////////////////
110
111double RooNovosibirsk::analyticalIntegral(Int_t code, const char* rangeName) const
112{
113 assert(code==1 || code==2) ;
114
115 //The range is defined as [A,B]
116
117 //Numerical values need for the evaluation of the integral
118 static const double sqrt2 = 1.4142135623730950; // Sqrt(2)
119 static const double sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
120 static const double sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
121 static const double log4 = 1.38629436111989062; //Log(2)
122 static const double rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
123 static const double sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
124
125 if (code==1) {
126 double A = x.min(rangeName);
127 double B = x.max(rangeName);
128
129 double result = 0;
130
131
132 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
133 if (std::abs(tail) < 1.e-7) {
134
135 double xscale = sqrt2*width;
136
137 result = rootpiby2*width*(std::erf((B-peak)/xscale)-std::erf((A-peak)/xscale));
138
139 return result;
140
141 }
142
143 // If the range is not defined correctly the function becomes complex
144 double log_argument_A = ( (peak - A)*tail + width ) / width ;
145 double log_argument_B = ( (peak - B)*tail + width ) / width ;
146
147 //lower limit
148 if ( log_argument_A < 1.e-7) {
149 log_argument_A = 1.e-7;
150 }
151
152 //upper limit
153 if ( log_argument_B < 1.e-7) {
154 log_argument_B = 1.e-7;
155 }
156
157 double term1 = TMath::ASinH( tail * sqlog4 );
158 double term1_2 = term1 * term1;
159
160 //Calculate the error function arguments
161 double erf_termA = ( term1_2 - log4 * std::log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
162 double erf_termB = ( term1_2 - log4 * std::log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
163
164 result = 0.5 / tail * width * term1 * ( std::erf(erf_termB) - std::erf(erf_termA)) * sqpibylog2;
165
166 return result;
167
168 } else if (code==2) {
169 double A = x.min(rangeName);
170 double B = x.max(rangeName);
171
172 double result = 0;
173
174
175 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
176 if (std::abs(tail) < 1.e-7) {
177
178 double xscale = sqrt2*width;
179
180 result = rootpiby2*width*(std::erf((B-x)/xscale)-std::erf((A-x)/xscale));
181
182 return result;
183
184 }
185
186 // If the range is not defined correctly the function becomes complex
187 double log_argument_A = ( (A - x)*tail + width ) / width;
188 double log_argument_B = ( (B - x)*tail + width ) / width;
189
190 //lower limit
191 if ( log_argument_A < 1.e-7) {
192 log_argument_A = 1.e-7;
193 }
194
195 //upper limit
196 if ( log_argument_B < 1.e-7) {
197 log_argument_B = 1.e-7;
198 }
199
200 double term1 = TMath::ASinH( tail * sqlog4 );
201 double term1_2 = term1 * term1;
202
203 //Calculate the error function arguments
204 double erf_termA = ( term1_2 - log4 * std::log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
205 double erf_termB = ( term1_2 - log4 * std::log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
206
207 result = 0.5 / tail * width * term1 * ( std::erf(erf_termB) - std::erf(erf_termA)) * sqpibylog2;
208
209 return result;
210
211 }
212
213 // Emit fatal error
214 coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;
215
216 // Put dummy return here to avoid compiler warnings
217 return 1.0 ;
218}
#define coutF(a)
#define ClassImp(name)
Definition Rtypes.h:382
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
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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:24
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
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.
RooRealProxy x
void doEval(RooFit::EvalContext &) const override
Compute multiple values of Novosibirsk distribution.
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.
Double_t x[n]
Definition legend1.C:17
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan extraArgs={})
Double_t ASinH(Double_t)
Returns the area hyperbolic sine of x.
Definition TMath.cxx:67