Logo ROOT   6.10/09
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 
23 RooNovosibirsk implements the Novosibirsk function
24 
25 Function taken from H. Ikeda et al. NIM A441 (2000), p. 401 (Belle Collaboration)
26 
27 **/
28 
29 #include "RooFit.h"
30 
31 #include <math.h>
32 #include "TMath.h"
33 
34 #include "RooNovosibirsk.h"
35 #include "RooRealVar.h"
36 
37 using namespace std;
38 
40 
41 ////////////////////////////////////////////////////////////////////////////////
42 
43 RooNovosibirsk::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 
95 Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
96 {
97  if (matchArgs(allVars,analVars,x)) return 1 ;
98  if (matchArgs(allVars,analVars,peak)) return 2 ;
99 
100  //The other two integrals over tali and width are not integrable
101 
102  return 0 ;
103 }
104 
105 ////////////////////////////////////////////////////////////////////////////////
106 
107 Double_t RooNovosibirsk::analyticalIntegral(Int_t code, const char* rangeName) const
108 {
109  assert(code==1 || code==2) ;
110 
111  //The range is defined as [A,B]
112 
113  //Numerical values need for the evaluation of the integral
114  static const Double_t sqrt2 = 1.4142135623730950; // Sqrt(2)
115  static const Double_t sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
116  static const Double_t sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
117  static const Double_t log4 = 1.38629436111989062; //Log(2)
118  static const Double_t rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
119  static const Double_t sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
120 
121  if (code==1) {
122  Double_t A = x.min(rangeName);
123  Double_t B = x.max(rangeName);
124 
125  Double_t result = 0;
126 
127 
128  //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
129  if (TMath::Abs(tail) < 1.e-7) {
130 
131  Double_t xscale = sqrt2*width;
132 
133  result = rootpiby2*width*(TMath::Erf((B-peak)/xscale)-TMath::Erf((A-peak)/xscale));
134 
135  return result;
136 
137  }
138 
139  // If the range is not defined correctly the function becomes complex
140  Double_t log_argument_A = ( (peak - A)*tail + width ) / width ;
141  Double_t log_argument_B = ( (peak - B)*tail + width ) / width ;
142 
143  //lower limit
144  if ( log_argument_A < 1.e-7) {
145  log_argument_A = 1.e-7;
146  }
147 
148  //upper limit
149  if ( log_argument_B < 1.e-7) {
150  log_argument_B = 1.e-7;
151  }
152 
153  Double_t term1 = TMath::ASinH( tail * sqlog4 );
154  Double_t term1_2 = term1 * term1;
155 
156  //Calculate the error function arguments
157  Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
158  Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
159 
160  result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
161 
162  return result;
163 
164  } else if (code==2) {
165  Double_t A = x.min(rangeName);
166  Double_t B = x.max(rangeName);
167 
168  Double_t result = 0;
169 
170 
171  //If tail==0 the function becomes gaussian, thus we return a Gaussian integral
172  if (TMath::Abs(tail) < 1.e-7) {
173 
174  Double_t xscale = sqrt2*width;
175 
176  result = rootpiby2*width*(TMath::Erf((B-x)/xscale)-TMath::Erf((A-x)/xscale));
177 
178  return result;
179 
180  }
181 
182  // If the range is not defined correctly the function becomes complex
183  Double_t log_argument_A = ( (A - x)*tail + width ) / width;
184  Double_t log_argument_B = ( (B - x)*tail + width ) / width;
185 
186  //lower limit
187  if ( log_argument_A < 1.e-7) {
188  log_argument_A = 1.e-7;
189  }
190 
191  //upper limit
192  if ( log_argument_B < 1.e-7) {
193  log_argument_B = 1.e-7;
194  }
195 
196  Double_t term1 = TMath::ASinH( tail * sqlog4 );
197  Double_t term1_2 = term1 * term1;
198 
199  //Calculate the error function arguments
200  Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
201  Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
202 
203  result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
204 
205  return result;
206 
207  }
208 
209  // Emit fatal error
210  coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;
211 
212  // Put dummy return here to avoid compiler warnings
213  return 1.0 ;
214 }
RooRealProxy peak
static double B[]
Double_t Log(Double_t x)
Definition: TMath.h:649
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooNovosibirsk implements the Novosibirsk function.
int Int_t
Definition: RtypesCore.h:41
STL namespace.
static double A[]
RooRealProxy width
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
Double_t x[n]
Definition: legend1.C:17
Double_t ASinH(Double_t)
Definition: TMath.cxx:67
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:187
Double_t evaluate() const
If tail=eta=0 the Belle distribution becomes gaussian.
RooRealProxy tail
#define coutF(a)
Definition: RooMsgService.h:35
RooRealProxy x
Double_t Exp(Double_t x)
Definition: TMath.h:622
#define ClassImp(name)
Definition: Rtypes.h:336
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
double result[121]
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
double log(double)