Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 #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>
37 using namespace std;
38 
39 ClassImp(RooNovosibirsk);
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 
58 RooNovosibirsk::RooNovosibirsk(const RooNovosibirsk& other, const char *name):
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 
70 Double_t RooNovosibirsk::evaluate() const
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 namespace {
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  */
106 template<class Tx, class Twidth, class Tpeak, class Ttail>
107 void 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 
130 RooSpan<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(),
142  BracketAdapter<double> (peak),
143  BracketAdapter<double> (width),
144  BracketAdapter<double> (tail));
145  }
146  else {
147  compute(info.size, output.data(),
148  BracketAdapterWithMask (x,x.getValBatch(begin,info.size)),
149  BracketAdapterWithMask (peak,peak.getValBatch(begin,info.size)),
150  BracketAdapterWithMask (width,width.getValBatch(begin,info.size)),
151  BracketAdapterWithMask (tail,tail.getValBatch(begin,info.size)));
152  }
153  return output;
154 }
155 
156 ////////////////////////////////////////////////////////////////////////////////
157 
158 Int_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 
170 Double_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 }