Logo ROOT   6.08/07
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 /**
21 \file RooNovosibirsk.cxx
22 \class RooNovosibirsk
23 \ingroup Roofit
24 
25 
26 RooNovosibirsk implements the Novosibirsk function
27 
28 Function taken from H. Ikeda et al. NIM A441 (2000), p. 401 (Belle Collaboration)
29 
30 **/
31 
32 
33 #include "RooFit.h"
34 
35 #include <math.h>
36 #include "TMath.h"
37 
38 #include "RooNovosibirsk.h"
39 #include "RooRealVar.h"
40 
41 using namespace std;
42 
44 
45 
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 
49 RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
50  RooAbsReal& _x, RooAbsReal& _peak,
51  RooAbsReal& _width, RooAbsReal& _tail) :
52  // The two addresses refer to our first dependent variable and
53  // parameter, respectively, as declared in the rdl file
54  RooAbsPdf(name, title),
55  x("x","x",this,_x),
56  width("width","width",this,_width),
57  peak("peak","peak",this,_peak),
58  tail("tail","tail",this,_tail)
59 {
60 }
61 
62 
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 
67  RooAbsPdf(other,name),
68  x("x",this,other.x),
69  width("width",this,other.width),
70  peak("peak",this,other.peak),
71  tail("tail",this,other.tail)
72 {
73 }
74 
75 
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 ///If tail=eta=0 the Belle distribution becomes gaussian
79 
81 {
82  if (TMath::Abs(tail) < 1.e-7) {
83  return TMath::Exp( -0.5 * TMath::Power( ( (x - peak) / width ), 2 ));
84  }
85 
86  Double_t arg = 1.0 - ( x - peak ) * tail / width;
87 
88  if (arg < 1.e-7) {
89  //Argument of logaritem negative. Real continuation -> function equals zero
90  return 0.0;
91  }
92 
93  Double_t log = TMath::Log(arg);
94  static const Double_t xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
95 
96  Double_t width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
97  Double_t width_zero2 = width_zero * width_zero;
98  Double_t exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
99 
100  return TMath::Exp(exponent) ;
101 }
102 
103 
104 
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 
108 Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
109 {
110  if (matchArgs(allVars,analVars,x)) return 1 ;
111  if (matchArgs(allVars,analVars,peak)) return 2 ;
112 
113  //The other two integrals over tali and width are not integrable
114 
115  return 0 ;
116 }
117 
118 
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 
122 Double_t RooNovosibirsk::analyticalIntegral(Int_t code, const char* rangeName) const
123 {
124  assert(code==1 || code==2) ;
125 
126  //The range is defined as [A,B]
127 
128  //Numerical values need for the evaluation of the integral
129  static const Double_t sqrt2 = 1.4142135623730950; // Sqrt(2)
130  static const Double_t sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
131  static const Double_t sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
132  static const Double_t log4 = 1.38629436111989062; //Log(2)
133  static const Double_t rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
134  static const Double_t sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
135 
136  if (code==1) {
137  Double_t A = x.min(rangeName);
138  Double_t B = x.max(rangeName);
139 
140  Double_t result = 0;
141 
142 
143  //If tail==0 the function becomes gaussian, thus we return a gassian integral
144  if (TMath::Abs(tail) < 1.e-7) {
145 
146  Double_t xscale = sqrt2*width;
147 
148  result = rootpiby2*width*(TMath::Erf((B-peak)/xscale)-TMath::Erf((A-peak)/xscale));
149 
150  return result;
151 
152  }
153 
154  // If the range is not defined correctly the function becomes complex
155  Double_t log_argument_A = ( (peak - A)*tail + width ) / width ;
156  Double_t log_argument_B = ( (peak - B)*tail + width ) / width ;
157 
158  //lower limit
159  if ( log_argument_A < 1.e-7) {
160  log_argument_A = 1.e-7;
161  }
162 
163  //upper limit
164  if ( log_argument_B < 1.e-7) {
165  log_argument_B = 1.e-7;
166  }
167 
168  Double_t term1 = TMath::ASinH( tail * sqlog4 );
169  Double_t term1_2 = term1 * term1;
170 
171  //Calculate the error function arguments
172  Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
173  Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
174 
175  result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
176 
177  return result;
178 
179  } else if (code==2) {
180  Double_t A = x.min(rangeName);
181  Double_t B = x.max(rangeName);
182 
183  Double_t result = 0;
184 
185 
186  //If tail==0 the function becomes gaussian, thus we return a gassian integral
187  if (TMath::Abs(tail) < 1.e-7) {
188 
189  Double_t xscale = sqrt2*width;
190 
191  result = rootpiby2*width*(TMath::Erf((B-x)/xscale)-TMath::Erf((A-x)/xscale));
192 
193  return result;
194 
195  }
196 
197  // If the range is not defined correctly the function becomes complex
198  Double_t log_argument_A = ( (A - x)*tail + width ) / width;
199  Double_t log_argument_B = ( (B - x)*tail + width ) / width;
200 
201  //lower limit
202  if ( log_argument_A < 1.e-7) {
203  log_argument_A = 1.e-7;
204  }
205 
206  //upper limit
207  if ( log_argument_B < 1.e-7) {
208  log_argument_B = 1.e-7;
209  }
210 
211  Double_t term1 = TMath::ASinH( tail * sqlog4 );
212  Double_t term1_2 = term1 * term1;
213 
214  //Calculate the error function arguments
215  Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
216  Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
217 
218  result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
219 
220  return result;
221 
222  }
223 
224  // Emit fatal error
225  coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;
226 
227  // Put dummy return here to avoid compiler warnings
228  return 1.0 ;
229 }
RooRealProxy peak
static double B[]
Double_t Log(Double_t x)
Definition: TMath.h:526
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:110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
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:36
RooRealProxy x
Double_t Exp(Double_t x)
Definition: TMath.h:495
#define ClassImp(name)
Definition: Rtypes.h:279
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. ...
char name[80]
Definition: TGX11.cxx:109
double log(double)