Logo ROOT  
Reference Guide
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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
36
37////////////////////////////////////////////////////////////////////////////////
38
39RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
42 // The two addresses refer to our first dependent variable and
43 // parameter, respectively, as declared in the rdl file
44 RooAbsPdf(name, title),
45 x("x","x",this,_x),
46 width("width","width",this,_width),
47 peak("peak","peak",this,_peak),
48 tail("tail","tail",this,_tail)
49{
50}
51
52////////////////////////////////////////////////////////////////////////////////
53
56 x("x",this,other.x),
57 width("width",this,other.width),
58 peak("peak",this,other.peak),
59 tail("tail",this,other.tail)
60{
61}
62
63////////////////////////////////////////////////////////////////////////////////
64///If tail=eta=0 the Belle distribution becomes gaussian
65
67{
68 if (std::abs(tail) < 1.e-7) {
69 return std::exp( -0.5 * std::pow( ( (x - peak) / width ), 2 ));
70 }
71
72 double arg = 1.0 - ( x - peak ) * tail / width;
73
74 if (arg < 1.e-7) {
75 //Argument of logarithm negative. Real continuation -> function equals zero
76 return 0.0;
77 }
78
79 double log = std::log(arg);
80 static const double xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
81
82 double width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
84 double exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
85
86 return std::exp(exponent) ;
87}
88////////////////////////////////////////////////////////////////////////////////
89/// Compute multiple values of Novosibirsk distribution.
91{
93 {ctx.at(x), ctx.at(peak), ctx.at(width), ctx.at(tail)});
94}
95
96////////////////////////////////////////////////////////////////////////////////
97
99{
100 if (matchArgs(allVars,analVars,x)) return 1 ;
101 if (matchArgs(allVars,analVars,peak)) return 2 ;
102
103 //The other two integrals over tali and width are not integrable
104
105 return 0 ;
106}
107
108////////////////////////////////////////////////////////////////////////////////
109
111{
112 assert(code==1 || code==2) ;
113
114 //The range is defined as [A,B]
115
116 //Numerical values need for the evaluation of the integral
117 static const double sqrt2 = 1.4142135623730950; // Sqrt(2)
118 static const double sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
119 static const double sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
120 static const double log4 = 1.38629436111989062; //Log(2)
121 static const double rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
122 static const double sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
123
124 if (code==1) {
125 double A = x.min(rangeName);
126 double B = x.max(rangeName);
127
128 double result = 0;
129
130
131 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
132 if (std::abs(tail) < 1.e-7) {
133
134 double xscale = sqrt2*width;
135
136 result = rootpiby2*width*(std::erf((B-peak)/xscale)-std::erf((A-peak)/xscale));
137
138 return result;
139
140 }
141
142 // If the range is not defined correctly the function becomes complex
143 double log_argument_A = ( (peak - A)*tail + width ) / width ;
144 double log_argument_B = ( (peak - B)*tail + width ) / width ;
145
146 //lower limit
147 if ( log_argument_A < 1.e-7) {
148 log_argument_A = 1.e-7;
149 }
150
151 //upper limit
152 if ( log_argument_B < 1.e-7) {
153 log_argument_B = 1.e-7;
154 }
155
156 double term1 = TMath::ASinH( tail * sqlog4 );
157 double term1_2 = term1 * term1;
158
159 //Calculate the error function arguments
160 double erf_termA = ( term1_2 - log4 * std::log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
161 double erf_termB = ( term1_2 - log4 * std::log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
162
163 result = 0.5 / tail * width * term1 * ( std::erf(erf_termB) - std::erf(erf_termA)) * sqpibylog2;
164
165 return result;
166
167 } else if (code==2) {
168 double A = x.min(rangeName);
169 double B = x.max(rangeName);
170
171 double result = 0;
172
173
174 //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
175 if (std::abs(tail) < 1.e-7) {
176
177 double xscale = sqrt2*width;
178
179 result = rootpiby2*width*(std::erf((B-x)/xscale)-std::erf((A-x)/xscale));
180
181 return result;
182
183 }
184
185 // If the range is not defined correctly the function becomes complex
186 double log_argument_A = ( (A - x)*tail + width ) / width;
187 double log_argument_B = ( (B - x)*tail + width ) / width;
188
189 //lower limit
190 if ( log_argument_A < 1.e-7) {
191 log_argument_A = 1.e-7;
192 }
193
194 //upper limit
195 if ( log_argument_B < 1.e-7) {
196 log_argument_B = 1.e-7;
197 }
198
199 double term1 = TMath::ASinH( tail * sqlog4 );
200 double term1_2 = term1 * term1;
201
202 //Calculate the error function arguments
203 double erf_termA = ( term1_2 - log4 * std::log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
204 double erf_termB = ( term1_2 - log4 * std::log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
205
206 result = 0.5 / tail * width * term1 * ( std::erf(erf_termB) - std::erf(erf_termA)) * sqpibylog2;
207
208 return result;
209
210 }
211
212 // Emit fatal error
213 coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;
214
215 // Put dummy return here to avoid compiler warnings
216 return 1.0 ;
217}
#define coutF(a)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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