ROOT  6.06/09
Reference Guide
RooGaussModel.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitModels *
4  * @(#)root/roofit:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 //////////////////////////////////////////////////////////////////////////////
18 //
19 // BEGIN_HTML
20 // Class RooGaussModel implements a RooResolutionModel that models a Gaussian
21 // distribution. Object of class RooGaussModel can be used
22 // for analytical convolutions with classes inheriting from RooAbsAnaConvPdf
23 // END_HTML
24 //
25 
26 #include "RooFit.h"
27 
28 #include "TMath.h"
29 #include "Riostream.h"
30 #include "Riostream.h"
31 #include "RooGaussModel.h"
32 #include "RooRealConstant.h"
33 #include "RooRandom.h"
34 
35 #include "TError.h"
36 
37 using namespace std;
38 
39 ClassImp(RooGaussModel)
40 ;
41 
42 
43 
44 ////////////////////////////////////////////////////////////////////////////////
45 
46 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn,
47  RooAbsReal& _mean, RooAbsReal& _sigma) :
48  RooResolutionModel(name,title,xIn),
49  _flatSFInt(kFALSE),
50  _asympInt(kFALSE),
51  mean("mean","Mean",this,_mean),
52  sigma("sigma","Width",this,_sigma),
53  msf("msf","Mean Scale Factor",this,(RooRealVar&)RooRealConstant::value(1)),
54  ssf("ssf","Sigma Scale Factor",this,(RooRealVar&)RooRealConstant::value(1))
55 {
56 }
57 
58 
59 
60 ////////////////////////////////////////////////////////////////////////////////
61 
62 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn,
63  RooAbsReal& _mean, RooAbsReal& _sigma,
64  RooAbsReal& _msSF) :
65  RooResolutionModel(name,title,xIn),
66  _flatSFInt(kFALSE),
67  _asympInt(kFALSE),
68  mean("mean","Mean",this,_mean),
69  sigma("sigma","Width",this,_sigma),
70  msf("msf","Mean Scale Factor",this,_msSF),
71  ssf("ssf","Sigma Scale Factor",this,_msSF)
72 {
73 }
74 
75 
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 
79 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn,
80  RooAbsReal& _mean, RooAbsReal& _sigma,
81  RooAbsReal& _meanSF, RooAbsReal& _sigmaSF) :
82  RooResolutionModel(name,title,xIn),
83  _flatSFInt(kFALSE),
84  _asympInt(kFALSE),
85  mean("mean","Mean",this,_mean),
86  sigma("sigma","Width",this,_sigma),
87  msf("msf","Mean Scale Factor",this,_meanSF),
88  ssf("ssf","Sigma Scale Factor",this,_sigmaSF)
89 {
90 }
91 
92 
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 
96 RooGaussModel::RooGaussModel(const RooGaussModel& other, const char* name) :
97  RooResolutionModel(other,name),
98  _flatSFInt(other._flatSFInt),
99  _asympInt(other._asympInt),
100  mean("mean",this,other.mean),
101  sigma("sigma",this,other.sigma),
102  msf("msf",this,other.msf),
103  ssf("ssf",this,other.ssf)
104 {
105 }
106 
107 
108 
109 ////////////////////////////////////////////////////////////////////////////////
110 /// Destructor
111 
112 RooGaussModel::~RooGaussModel()
113 {
114 }
115 
116 
117 
118 ////////////////////////////////////////////////////////////////////////////////
119 
120 Int_t RooGaussModel::basisCode(const char* name) const
121 {
122  if (!TString("exp(-@0/@1)").CompareTo(name)) return expBasisPlus ;
123  if (!TString("exp(@0/@1)").CompareTo(name)) return expBasisMinus ;
124  if (!TString("exp(-abs(@0)/@1)").CompareTo(name)) return expBasisSum ;
125  if (!TString("exp(-@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisPlus ;
126  if (!TString("exp(@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisMinus ;
127  if (!TString("exp(-abs(@0)/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisSum ;
128  if (!TString("exp(-@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisPlus ;
129  if (!TString("exp(@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisMinus ;
130  if (!TString("exp(-abs(@0)/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisSum ;
131  if (!TString("(@0/@1)*exp(-@0/@1)").CompareTo(name)) return linBasisPlus ;
132  if (!TString("(@0/@1)*(@0/@1)*exp(-@0/@1)").CompareTo(name)) return quadBasisPlus ;
133  if (!TString("exp(-@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisPlus;
134  if (!TString("exp(@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisMinus;
135  if (!TString("exp(-abs(@0)/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisSum;
136  if (!TString("exp(-@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisPlus;
137  if (!TString("exp(@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisMinus;
138  if (!TString("exp(-abs(@0)/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisSum;
139  return 0 ;
140 }
141 
142 
143 
144 ////////////////////////////////////////////////////////////////////////////////
145 ///cout << "RooGaussModel::evaluate(" << GetName() << ") basisCode = " << _basisCode << endl ;
146 
147 Double_t RooGaussModel::evaluate() const
148 {
149  // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
150  static Double_t root2(std::sqrt(2.)) ;
151  static Double_t root2pi(std::sqrt(2.*std::atan2(0.,-1.))) ;
152  static Double_t rootpi(std::sqrt(std::atan2(0.,-1.))) ;
153 
154  BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
155  BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
156 
157  Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
158  if (basisType == coshBasis && _basisCode!=noBasis ) {
159  Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
160  if (dGamma==0) basisType = expBasis;
161  }
162 
163  if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
164  Double_t xprime = (x-(mean*msf))/(sigma*ssf) ;
165  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 1st form" << endl ;
166 
167  Double_t result = std::exp(-0.5*xprime*xprime)/(sigma*ssf*root2pi) ;
168  if (_basisCode!=0 && basisSign==Both) result *= 2 ;
169  return result ;
170  }
171 
172  // *** 2nd form: 0, used for sinBasis, linBasis, and quadBasis with tau=0 ***
173  if (tau==0) {
174  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 2nd form" << endl ;
175  return 0. ;
176  }
177 
178  // *** 3nd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
179  Double_t omega = (basisType==sinBasis || basisType==cosBasis) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
180  Double_t dgamma = (basisType==sinhBasis || basisType==coshBasis) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
181  Double_t _x = omega *tau ;
182  Double_t _y = tau*dgamma/2;
183  Double_t xprime = (x-(mean*msf))/tau ;
184  Double_t c = (sigma*ssf)/(root2*tau) ;
185  Double_t u = xprime/(2*c) ;
186 
187  if (basisType==expBasis || (basisType==cosBasis && _x==0.)) {
188  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 3d form tau=" << tau << endl ;
189  Double_t result(0) ;
190  if (basisSign!=Minus) result += evalCerf(0,-u,c).real();
191  if (basisSign!=Plus) result += evalCerf(0, u,c).real();
192  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 1 " << endl; }
193  return result ;
194  }
195 
196  // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
197  if (basisType==sinBasis) {
198  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 4th form omega = " << omega << ", tau = " << tau << endl ;
199  Double_t result(0) ;
200  if (_x==0.) return result ;
201  if (basisSign!=Minus) result += -evalCerf(-_x,-u,c).imag();
202  if (basisSign!=Plus) result += -evalCerf( _x, u,c).imag();
203  if (TMath::IsNaN(result)) cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 3 " << endl;
204  return result ;
205  }
206 
207  // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
208  if (basisType==cosBasis) {
209  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
210  Double_t result(0) ;
211  if (basisSign!=Minus) result += evalCerf(-_x,-u,c).real();
212  if (basisSign!=Plus) result += evalCerf( _x, u,c).real();
213  if (TMath::IsNaN(result)) cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 4 " << endl;
214  return result ;
215  }
216 
217  // ***8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasisSum ***
218  if (basisType==coshBasis || basisType ==sinhBasis) {
219  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 8th form tau = " << tau << endl ;
220  Double_t result(0);
221  int sgn = ( basisType == coshBasis ? +1 : -1 );
222  if (basisSign!=Minus) result += 0.5*( evalCerf(0,-u,c*(1-_y)).real()+sgn*evalCerf(0,-u,c*(1+_y)).real()) ;
223  if (basisSign!=Plus) result += 0.5*(sgn*evalCerf(0, u,c*(1-_y)).real()+ evalCerf(0, u,c*(1+_y)).real()) ;
224  if (TMath::IsNaN(result)) cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 8 " << endl;
225  return result ;
226  }
227 
228  // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
229  if (basisType==linBasis) {
230  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 6th form tau = " << tau << endl ;
231  R__ASSERT(basisSign==Plus); // This should only be for positive times
232 
233  Double_t f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
234  Double_t f1 = std::exp(-u*u);
235  return (xprime - 2*c*c)*f0 + (2*c/rootpi)*f1 ;
236  }
237 
238  // *** 7th form: Convolution with (t/tau)^2*exp(-t/tau), used for quadBasis ***
239  if (basisType==quadBasis) {
240  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 7th form tau = " << tau << endl ;
241  R__ASSERT(basisSign==Plus); // This should only be for positive times
242 
243  Double_t f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
244  Double_t f1 = std::exp(-u*u);
245  Double_t x2c2 = xprime - 2*c*c;
246  return ( x2c2*x2c2*f0 + (2*c/rootpi)*x2c2*f1 + 2*c*c*f0 );
247  }
248 
249  R__ASSERT(0) ;
250  return 0 ;
251 }
252 
253 
254 
255 ////////////////////////////////////////////////////////////////////////////////
256 
257 Int_t RooGaussModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
258 {
259  switch(_basisCode) {
260 
261  // Analytical integration capability of raw PDF
262  case noBasis:
263 
264  // Optionally advertise flat integral over sigma scale factor
265  if (_flatSFInt) {
266  if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) return 2 ;
267  }
268 
269  if (matchArgs(allVars,analVars,convVar())) return 1 ;
270  break ;
271 
272  // Analytical integration capability of convoluted PDF
273  case expBasisPlus:
274  case expBasisMinus:
275  case expBasisSum:
276  case sinBasisPlus:
277  case sinBasisMinus:
278  case sinBasisSum:
279  case cosBasisPlus:
280  case cosBasisMinus:
281  case cosBasisSum:
282  case linBasisPlus:
283  case quadBasisPlus:
284  case coshBasisMinus:
285  case coshBasisPlus:
286  case coshBasisSum:
287  case sinhBasisMinus:
288  case sinhBasisPlus:
289  case sinhBasisSum:
290 
291  // Optionally advertise flat integral over sigma scale factor
292  if (_flatSFInt) {
293 
294  if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) {
295  return 2 ;
296  }
297  }
298 
299  if (matchArgs(allVars,analVars,convVar())) return 1 ;
300  break ;
301  }
302 
303  return 0 ;
304 }
305 
306 
307 
308 ////////////////////////////////////////////////////////////////////////////////
309 
310 Double_t RooGaussModel::analyticalIntegral(Int_t code, const char* rangeName) const
311 {
312  static const Double_t root2 = std::sqrt(2.) ;
313  //static Double_t rootPiBy2 = std::sqrt(std::atan2(0.0,-1.0)/2.0);
314  static const Double_t rootpi = std::sqrt(std::atan2(0.0,-1.0));
315  Double_t ssfInt(1.0) ;
316 
317  // Code must be 1 or 2
318  R__ASSERT(code==1||code==2) ;
319  if (code==2) ssfInt = (ssf.max(rangeName)-ssf.min(rangeName)) ;
320 
321  BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
322  BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
323 
324  // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
325  Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
326  if (basisType == coshBasis && _basisCode!=noBasis ) {
327  Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
328  if (dGamma==0) basisType = expBasis;
329  }
330  if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
331  Double_t xscale = root2*(sigma*ssf);
332  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 1st form" << endl ;
333 
334  Double_t xpmin = (x.min(rangeName)-(mean*msf))/xscale ;
335  Double_t xpmax = (x.max(rangeName)-(mean*msf))/xscale ;
336 
337  Double_t result ;
338  if (_asympInt) { // modified FMV, 07/24/03
339  result = 1.0 ;
340  } else {
341  result = 0.5*(RooMath::erf(xpmax)-RooMath::erf(xpmin)) ;
342  }
343 
344  if (_basisCode!=0 && basisSign==Both) result *= 2 ;
345  //cout << "Integral 1st form " << " result= " << result*ssfInt << endl;
346  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 1 " << endl; }
347  return result*ssfInt ;
348  }
349 
350 
351  Double_t omega = ((basisType==sinBasis)||(basisType==cosBasis)) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
352  Double_t dgamma =((basisType==sinhBasis)||(basisType==coshBasis)) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
353 
354  // *** 2nd form: unity, used for sinBasis and linBasis with tau=0 (PDF is zero) ***
355  if (tau==0) {
356  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 2nd form" << endl ;
357  return 0. ;
358  }
359 
360  // *** 3rd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
361  Double_t c = (sigma*ssf)/(root2*tau) ;
362  Double_t xpmin = (x.min(rangeName)-(mean*msf))/tau ;
363  Double_t xpmax = (x.max(rangeName)-(mean*msf))/tau ;
364  Double_t umin = xpmin/(2*c) ;
365  Double_t umax = xpmax/(2*c) ;
366 
367  if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
368  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 3d form tau=" << tau << endl ;
369  Double_t result(0) ;
370  if (basisSign!=Minus) result += evalCerfInt(+1,0,tau,-umin,-umax,c).real();
371  if (basisSign!=Plus) result += evalCerfInt(-1,0,tau, umin, umax,c).real();
372  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 3 " << endl; }
373  return result*ssfInt ;
374  }
375 
376  // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
377  Double_t _x = omega * tau ;
378  Double_t _y = tau*dgamma/2;
379 
380  if (basisType==sinBasis) {
381  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 4th form omega = " << omega << ", tau = " << tau << endl ;
382  Double_t result(0) ;
383  if (_x==0) return result*ssfInt ;
384  if (basisSign!=Minus) result += -1*evalCerfInt(+1,-_x,tau,-umin,-umax,c).imag();
385  if (basisSign!=Plus) result += -1*evalCerfInt(-1, _x,tau, umin, umax,c).imag();
386  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 4 " << endl; }
387  return result*ssfInt ;
388  }
389 
390  // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
391  if (basisType==cosBasis) {
392  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
393  Double_t result(0) ;
394  if (basisSign!=Minus) result += evalCerfInt(+1,-_x,tau,-umin,-umax,c).real();
395  if (basisSign!=Plus) result += evalCerfInt(-1, _x,tau, umin, umax,c).real();
396  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 5 " << endl; }
397  return result*ssfInt ;
398  }
399 
400  // *** 8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasis ***
401  // *** 9th form: Convolution with exp(-|t|/tau)*sinh(dgamma*t/2), used for sinhBasis ***
402  if (basisType==coshBasis || basisType == sinhBasis) {
403  if (verboseEval()>0) {cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 8th form tau=" << tau << endl ; }
404  Double_t result(0) ;
405  int sgn = ( basisType == coshBasis ? +1 : -1 );
406  if (basisSign!=Minus) result += 0.5*( evalCerfInt(+1,0,tau/(1-_y),-umin,-umax,c*(1-_y)).real()+ sgn*evalCerfInt(+1,0,tau/(1+_y),-umin,-umax,c*(1+_y)).real());
407  if (basisSign!=Plus) result += 0.5*(sgn*evalCerfInt(-1,0,tau/(1-_y), umin, umax,c*(1-_y)).real()+ evalCerfInt(-1,0,tau/(1+_y), umin, umax,c*(1+_y)).real());
408  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 6 " << endl; }
409  return result*ssfInt ;
410  }
411 
412  // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
413  if (basisType==linBasis) {
414  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 6th form tau=" << tau << endl ;
415 
416  Double_t f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
417  Double_t f1 = std::exp(-umax*umax) - std::exp(-umin*umin);
418 
419  Double_t tmp1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
420  Double_t tmp2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
421 
422  Double_t f2 = tmp1 - tmp2;
423  Double_t f3 = xpmax*tmp1 - xpmin*tmp2;
424 
425  Double_t expc2 = std::exp(c*c);
426 
427  return -tau*( f0 +
428  (2*c/rootpi)*f1 +
429  (1 - 2*c*c)*expc2*f2 +
430  expc2*f3
431  )*ssfInt;
432  }
433 
434  // *** 7th form: Convolution with (t/tau)*(t/tau)*exp(-t/tau), used for quadBasis ***
435  if (basisType==quadBasis) {
436  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 7th form tau=" << tau << endl ;
437 
438  Double_t f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
439 
440  Double_t tmpA1 = std::exp(-umax*umax);
441  Double_t tmpA2 = std::exp(-umin*umin);
442 
443  Double_t f1 = tmpA1 - tmpA2;
444  Double_t f2 = umax*tmpA1 - umin*tmpA2;
445 
446  Double_t tmpB1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
447  Double_t tmpB2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
448 
449  Double_t f3 = tmpB1 - tmpB2;
450  Double_t f4 = xpmax*tmpB1 - xpmin*tmpB2;
451  Double_t f5 = xpmax*xpmax*tmpB1 - xpmin*xpmin*tmpB2;
452 
453  Double_t expc2 = std::exp(c*c);
454 
455  return -tau*( 2*f0 +
456  (4*c/rootpi)*((1-c*c)*f1 + c*f2) +
457  (2*c*c*(2*c*c-1) + 2)*expc2*f3 - (4*c*c-2)*expc2*f4 + expc2*f5
458  )*ssfInt;
459  }
460  R__ASSERT(0) ;
461  return 0 ;
462 }
463 
464 ////////////////////////////////////////////////////////////////////////////////
465 /// use the approximation: erf(z) = exp(-z*z)/(std::sqrt(pi)*z)
466 /// to explicitly cancel the divergent exp(y*y) behaviour of
467 /// CWERF for z = x + i y with large negative y
468 
469 std::complex<Double_t> RooGaussModel::evalCerfApprox(Double_t _x, Double_t u, Double_t c)
470 {
471  static const Double_t rootpi= std::sqrt(std::atan2(0.,-1.));
472  const std::complex<Double_t> z(_x * c, u + c);
473  const std::complex<Double_t> zc(u + c, - _x * c);
474  const std::complex<Double_t> zsq((z.real() + z.imag()) * (z.real() - z.imag()),
475  2. * z.real() * z.imag());
476  const std::complex<Double_t> v(-zsq.real() - u*u, -zsq.imag());
477  const std::complex<Double_t> ev = std::exp(v);
478  const std::complex<Double_t> mez2zcrootpi = -std::exp(zsq)/(zc*rootpi);
479 
480  return 2. * (ev * (mez2zcrootpi + 1.));
481 }
482 
483 ////////////////////////////////////////////////////////////////////////////////
484 
485 std::complex<Double_t> RooGaussModel::evalCerfInt(Double_t sign, Double_t _x, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
486 {
487  std::complex<Double_t> diff(2., 0.);
488  if (!_asympInt) {
489  diff = evalCerf(_x,umin,c);
490  diff -= evalCerf(_x,umax,c);
491  diff += RooMath::erf(umin) - RooMath::erf(umax);
492  diff *= sign;
493  }
494  diff *= std::complex<Double_t>(1., _x);
495  diff *= tau / (1.+_x*_x);
496  return diff;
497 }
498 
499 ////////////////////////////////////////////////////////////////////////////////
500 
501 Int_t RooGaussModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
502 {
503  if (matchArgs(directVars,generateVars,x)) return 1 ;
504  return 0 ;
505 }
506 
507 ////////////////////////////////////////////////////////////////////////////////
508 
509 void RooGaussModel::generateEvent(Int_t code)
510 {
511  R__ASSERT(code==1) ;
512  Double_t xmin = x.min();
513  Double_t xmax = x.max();
514  TRandom *generator = RooRandom::randomGenerator();
515  while(true) {
516  Double_t xgen = generator->Gaus(mean*msf,sigma*ssf);
517  if (xgen<xmax && xgen>xmin) {
518  x = xgen ;
519  return ;
520  }
521  }
522 }
RooGaussModel::analyticalIntegral
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Definition: RooGaussModel.cxx:310
xmin
float xmin
Definition: THbookFile.cxx:93
RooGaussModel::getAnalyticalIntegral
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
Definition: RooGaussModel.cxx:257
RooGaussModel::mean
RooRealProxy mean
Definition: RooGaussModel.h:82
cxcoutE
#define cxcoutE(a)
Definition: RooMsgService.h:96
TRandom::Gaus
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
R__ASSERT
#define R__ASSERT(e)
Definition: TError.h:98
TString
Basic string class.
Definition: TString.h:137
RooGaussModel::evalCerf
static std::complex< Double_t > evalCerf(Double_t swt, Double_t u, Double_t c)
Definition: RooGaussModel.h:69
Int_t
int Int_t
Definition: RtypesCore.h:41
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
kFALSE
const Bool_t kFALSE
Definition: Rtypes.h:92
std
STL namespace.
RooAbsPdf::verboseEval
static int verboseEval()
Return global level of verbosity for p.d.f. evaluations.
Definition: RooAbsPdf.cxx:2950
RooGaussModel::evalCerfApprox
static std::complex< Double_t > evalCerfApprox(Double_t swt, Double_t u, Double_t c)
use the approximation: erf(z) = exp(-z*z)/(std::sqrt(pi)*z) to explicitly cancel the divergent exp(y*...
Definition: RooGaussModel.cxx:469
RooGaussModel::_asympInt
Bool_t _asympInt
Definition: RooGaussModel.h:80
RooRealProxy::arg
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
sqrt
double sqrt(double)
TRandom
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:29
RooRandom::randomGenerator
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:53
RooAbsArg::RooArgSet
friend class RooArgSet
Definition: RooAbsArg.h:469
RooAbsReal::getVal
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
RooGaussModel::~RooGaussModel
virtual ~RooGaussModel()
Destructor.
Definition: RooGaussModel.cxx:112
RooGaussModel::ssf
RooRealProxy ssf
Definition: RooGaussModel.h:85
RooGaussModel::_flatSFInt
Bool_t _flatSFInt
Definition: RooGaussModel.h:78
RooGaussModel::basisCode
virtual Int_t basisCode(const char *name) const
Definition: RooGaussModel.cxx:120
RooGaussModel::generateEvent
void generateEvent(Int_t code)
Interface for generation of anan event using the algorithm corresponding to the specified code...
Definition: RooGaussModel.cxx:509
return
return
Definition: TBase64.cxx:62
v
SVector< double, 2 > v
Definition: Dict.h:5
RooResolutionModel::convVar
RooRealVar & convVar() const
Return the convolution variable of the resolution model.
Definition: RooResolutionModel.cxx:264
CompareTo
return fString CompareTo(((TObjString *) obj) ->fString)
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
xmax
float xmax
Definition: THbookFile.cxx:93
RooMath::erfc
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
Definition: RooMath.cxx:560
RooGaussModel::evaluate
virtual Double_t evaluate() const
cout << "RooGaussModel::evaluate(" << GetName() << ") basisCode = " << _basisCode << endl ; ...
Definition: RooGaussModel.cxx:147
Double_t
double Double_t
Definition: RtypesCore.h:55
RooAbsReal
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
RooFormulaVar::getParameter
RooAbsArg * getParameter(const char *name) const
Definition: RooFormulaVar.h:39
RooResolutionModel::basis
const RooFormulaVar & basis() const
Definition: RooResolutionModel.h:51
atan2
double atan2(double, double)
name
#define name(a, b)
Definition: linkTestLib0.cpp:5
sample_config.c
tuple c
Definition: sample_config.py:3
RooRealProxy::min
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
f2
double f2(const double *x)
Definition: testIntegration.cxx:16
TMath::IsNaN
Int_t IsNaN(Double_t x)
Definition: TMath.h:617
RooGaussModel::getGenerator
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooGaussModel.cxx:501
f1
TF1 * f1
Definition: legend1.C:11
RooGaussModel::msf
RooRealProxy msf
Definition: RooGaussModel.h:84
RooMath::erf
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:584
RooGaussModel::evalCerfInt
std::complex< Double_t > evalCerfInt(Double_t sign, Double_t wt, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
Definition: RooGaussModel.cxx:485
result
double result[121]
Definition: testSampleQuantiles.cxx:17
RootCsg::sign
Int_t sign(Double_t x)
Definition: CsgOps.cxx:89
RooAbsReal::matchArgs
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Definition: RooAbsReal.cxx:3295
exp
double exp(double)
value
float value
Definition: math.cpp:443
RooRealProxy::max
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
RooGaussModel::sigma
RooRealProxy sigma
Definition: RooGaussModel.h:83
ClassImp
ClassImp(RooGaussModel)