Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 /** \class RooGaussModel
18  \ingroup Roofit
19 
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 **/
24 
25 #include "RooFit.h"
26 
27 #include "TMath.h"
28 #include "Riostream.h"
29 #include "Riostream.h"
30 #include "RooGaussModel.h"
31 #include "RooRealConstant.h"
32 #include "RooRandom.h"
33 
34 #include "TError.h"
35 
36 using namespace std;
37 
38 ClassImp(RooGaussModel);
39 
40 ////////////////////////////////////////////////////////////////////////////////
41 
42 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn,
43  RooAbsReal& _mean, RooAbsReal& _sigma) :
44  RooResolutionModel(name,title,xIn),
45  _flatSFInt(kFALSE),
46  _asympInt(kFALSE),
47  mean("mean","Mean",this,_mean),
48  sigma("sigma","Width",this,_sigma),
49  msf("msf","Mean Scale Factor",this,(RooRealVar&)RooRealConstant::value(1)),
50  ssf("ssf","Sigma Scale Factor",this,(RooRealVar&)RooRealConstant::value(1))
51 {
52 }
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 
56 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn,
57  RooAbsReal& _mean, RooAbsReal& _sigma,
58  RooAbsReal& _msSF) :
59  RooResolutionModel(name,title,xIn),
60  _flatSFInt(kFALSE),
61  _asympInt(kFALSE),
62  mean("mean","Mean",this,_mean),
63  sigma("sigma","Width",this,_sigma),
64  msf("msf","Mean Scale Factor",this,_msSF),
65  ssf("ssf","Sigma Scale Factor",this,_msSF)
66 {
67 }
68 
69 ////////////////////////////////////////////////////////////////////////////////
70 
71 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn,
72  RooAbsReal& _mean, RooAbsReal& _sigma,
73  RooAbsReal& _meanSF, RooAbsReal& _sigmaSF) :
74  RooResolutionModel(name,title,xIn),
75  _flatSFInt(kFALSE),
76  _asympInt(kFALSE),
77  mean("mean","Mean",this,_mean),
78  sigma("sigma","Width",this,_sigma),
79  msf("msf","Mean Scale Factor",this,_meanSF),
80  ssf("ssf","Sigma Scale Factor",this,_sigmaSF)
81 {
82 }
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 
86 RooGaussModel::RooGaussModel(const RooGaussModel& other, const char* name) :
87  RooResolutionModel(other,name),
88  _flatSFInt(other._flatSFInt),
89  _asympInt(other._asympInt),
90  mean("mean",this,other.mean),
91  sigma("sigma",this,other.sigma),
92  msf("msf",this,other.msf),
93  ssf("ssf",this,other.ssf)
94 {
95 }
96 
97 ////////////////////////////////////////////////////////////////////////////////
98 /// Destructor
99 
100 RooGaussModel::~RooGaussModel()
101 {
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 
106 Int_t RooGaussModel::basisCode(const char* name) const
107 {
108  if (!TString("exp(-@0/@1)").CompareTo(name)) return expBasisPlus ;
109  if (!TString("exp(@0/@1)").CompareTo(name)) return expBasisMinus ;
110  if (!TString("exp(-abs(@0)/@1)").CompareTo(name)) return expBasisSum ;
111  if (!TString("exp(-@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisPlus ;
112  if (!TString("exp(@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisMinus ;
113  if (!TString("exp(-abs(@0)/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisSum ;
114  if (!TString("exp(-@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisPlus ;
115  if (!TString("exp(@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisMinus ;
116  if (!TString("exp(-abs(@0)/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisSum ;
117  if (!TString("(@0/@1)*exp(-@0/@1)").CompareTo(name)) return linBasisPlus ;
118  if (!TString("(@0/@1)*(@0/@1)*exp(-@0/@1)").CompareTo(name)) return quadBasisPlus ;
119  if (!TString("exp(-@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisPlus;
120  if (!TString("exp(@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisMinus;
121  if (!TString("exp(-abs(@0)/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisSum;
122  if (!TString("exp(-@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisPlus;
123  if (!TString("exp(@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisMinus;
124  if (!TString("exp(-abs(@0)/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisSum;
125  return 0 ;
126 }
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 
130 Double_t RooGaussModel::evaluate() const
131 {
132  // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
133  static Double_t root2(std::sqrt(2.)) ;
134  static Double_t root2pi(std::sqrt(2.*std::atan2(0.,-1.))) ;
135  static Double_t rootpi(std::sqrt(std::atan2(0.,-1.))) ;
136 
137  BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
138  BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
139 
140  Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
141  if (basisType == coshBasis && _basisCode!=noBasis ) {
142  Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
143  if (dGamma==0) basisType = expBasis;
144  }
145 
146  if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
147  Double_t xprime = (x-(mean*msf))/(sigma*ssf) ;
148  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 1st form" << endl ;
149 
150  Double_t result = std::exp(-0.5*xprime*xprime)/(sigma*ssf*root2pi) ;
151  if (_basisCode!=0 && basisSign==Both) result *= 2 ;
152  return result ;
153  }
154 
155  // *** 2nd form: 0, used for sinBasis, linBasis, and quadBasis with tau=0 ***
156  if (tau==0) {
157  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 2nd form" << endl ;
158  return 0. ;
159  }
160 
161  // *** 3nd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
162  Double_t omega = (basisType==sinBasis || basisType==cosBasis) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
163  Double_t dgamma = (basisType==sinhBasis || basisType==coshBasis) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
164  Double_t _x = omega *tau ;
165  Double_t _y = tau*dgamma/2;
166  Double_t xprime = (x-(mean*msf))/tau ;
167  Double_t c = (sigma*ssf)/(root2*tau) ;
168  Double_t u = xprime/(2*c) ;
169 
170  if (basisType==expBasis || (basisType==cosBasis && _x==0.)) {
171  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 3d form tau=" << tau << endl ;
172  Double_t result(0) ;
173  if (basisSign!=Minus) result += evalCerf(0,-u,c).real();
174  if (basisSign!=Plus) result += evalCerf(0, u,c).real();
175  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 1 " << endl; }
176  return result ;
177  }
178 
179  // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
180  if (basisType==sinBasis) {
181  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 4th form omega = " << omega << ", tau = " << tau << endl ;
182  Double_t result(0) ;
183  if (_x==0.) return result ;
184  if (basisSign!=Minus) result += -evalCerf(-_x,-u,c).imag();
185  if (basisSign!=Plus) result += -evalCerf( _x, u,c).imag();
186  if (TMath::IsNaN(result)) cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 3 " << endl;
187  return result ;
188  }
189 
190  // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
191  if (basisType==cosBasis) {
192  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
193  Double_t result(0) ;
194  if (basisSign!=Minus) result += evalCerf(-_x,-u,c).real();
195  if (basisSign!=Plus) result += evalCerf( _x, u,c).real();
196  if (TMath::IsNaN(result)) cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 4 " << endl;
197  return result ;
198  }
199 
200  // ***8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasisSum ***
201  if (basisType==coshBasis || basisType ==sinhBasis) {
202  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 8th form tau = " << tau << endl ;
203  Double_t result(0);
204  int sgn = ( basisType == coshBasis ? +1 : -1 );
205  if (basisSign!=Minus) result += 0.5*( evalCerf(0,-u,c*(1-_y)).real()+sgn*evalCerf(0,-u,c*(1+_y)).real()) ;
206  if (basisSign!=Plus) result += 0.5*(sgn*evalCerf(0, u,c*(1-_y)).real()+ evalCerf(0, u,c*(1+_y)).real()) ;
207  if (TMath::IsNaN(result)) cxcoutE(Tracing) << "RooGaussModel::getVal(" << GetName() << ") got nan during case 8 " << endl;
208  return result ;
209  }
210 
211  // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
212  if (basisType==linBasis) {
213  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 6th form tau = " << tau << endl ;
214  R__ASSERT(basisSign==Plus); // This should only be for positive times
215 
216  Double_t f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
217  Double_t f1 = std::exp(-u*u);
218  return (xprime - 2*c*c)*f0 + (2*c/rootpi)*f1 ;
219  }
220 
221  // *** 7th form: Convolution with (t/tau)^2*exp(-t/tau), used for quadBasis ***
222  if (basisType==quadBasis) {
223  if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 7th form tau = " << tau << endl ;
224  R__ASSERT(basisSign==Plus); // This should only be for positive times
225 
226  Double_t f0 = std::exp(-xprime+c*c) * RooMath::erfc(-u+c);
227  Double_t f1 = std::exp(-u*u);
228  Double_t x2c2 = xprime - 2*c*c;
229  return ( x2c2*x2c2*f0 + (2*c/rootpi)*x2c2*f1 + 2*c*c*f0 );
230  }
231 
232  R__ASSERT(0) ;
233  return 0 ;
234 }
235 
236 ////////////////////////////////////////////////////////////////////////////////
237 
238 Int_t RooGaussModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
239 {
240  switch(_basisCode) {
241 
242  // Analytical integration capability of raw PDF
243  case noBasis:
244 
245  // Optionally advertise flat integral over sigma scale factor
246  if (_flatSFInt) {
247  if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) return 2 ;
248  }
249 
250  if (matchArgs(allVars,analVars,convVar())) return 1 ;
251  break ;
252 
253  // Analytical integration capability of convoluted PDF
254  case expBasisPlus:
255  case expBasisMinus:
256  case expBasisSum:
257  case sinBasisPlus:
258  case sinBasisMinus:
259  case sinBasisSum:
260  case cosBasisPlus:
261  case cosBasisMinus:
262  case cosBasisSum:
263  case linBasisPlus:
264  case quadBasisPlus:
265  case coshBasisMinus:
266  case coshBasisPlus:
267  case coshBasisSum:
268  case sinhBasisMinus:
269  case sinhBasisPlus:
270  case sinhBasisSum:
271 
272  // Optionally advertise flat integral over sigma scale factor
273  if (_flatSFInt) {
274 
275  if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) {
276  return 2 ;
277  }
278  }
279 
280  if (matchArgs(allVars,analVars,convVar())) return 1 ;
281  break ;
282  }
283 
284  return 0 ;
285 }
286 
287 ////////////////////////////////////////////////////////////////////////////////
288 
289 Double_t RooGaussModel::analyticalIntegral(Int_t code, const char* rangeName) const
290 {
291  static const Double_t root2 = std::sqrt(2.) ;
292  //static Double_t rootPiBy2 = std::sqrt(std::atan2(0.0,-1.0)/2.0);
293  static const Double_t rootpi = std::sqrt(std::atan2(0.0,-1.0));
294  Double_t ssfInt(1.0) ;
295 
296  // Code must be 1 or 2
297  R__ASSERT(code==1||code==2) ;
298  if (code==2) ssfInt = (ssf.max(rangeName)-ssf.min(rangeName)) ;
299 
300  BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
301  BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
302 
303  // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
304  Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
305  if (basisType == coshBasis && _basisCode!=noBasis ) {
306  Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
307  if (dGamma==0) basisType = expBasis;
308  }
309  if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
310  Double_t xscale = root2*(sigma*ssf);
311  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 1st form" << endl ;
312 
313  Double_t xpmin = (x.min(rangeName)-(mean*msf))/xscale ;
314  Double_t xpmax = (x.max(rangeName)-(mean*msf))/xscale ;
315 
316  Double_t result ;
317  if (_asympInt) { // modified FMV, 07/24/03
318  result = 1.0 ;
319  } else {
320  result = 0.5*(RooMath::erf(xpmax)-RooMath::erf(xpmin)) ;
321  }
322 
323  if (_basisCode!=0 && basisSign==Both) result *= 2 ;
324  //cout << "Integral 1st form " << " result= " << result*ssfInt << endl;
325  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 1 " << endl; }
326  return result*ssfInt ;
327  }
328 
329 
330  Double_t omega = ((basisType==sinBasis)||(basisType==cosBasis)) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
331  Double_t dgamma =((basisType==sinhBasis)||(basisType==coshBasis)) ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
332 
333  // *** 2nd form: unity, used for sinBasis and linBasis with tau=0 (PDF is zero) ***
334  if (tau==0) {
335  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 2nd form" << endl ;
336  return 0. ;
337  }
338 
339  // *** 3rd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
340  Double_t c = (sigma*ssf)/(root2*tau) ;
341  Double_t xpmin = (x.min(rangeName)-(mean*msf))/tau ;
342  Double_t xpmax = (x.max(rangeName)-(mean*msf))/tau ;
343  Double_t umin = xpmin/(2*c) ;
344  Double_t umax = xpmax/(2*c) ;
345 
346  if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
347  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 3d form tau=" << tau << endl ;
348  Double_t result(0) ;
349  if (basisSign!=Minus) result += evalCerfInt(+1,0,tau,-umin,-umax,c).real();
350  if (basisSign!=Plus) result += evalCerfInt(-1,0,tau, umin, umax,c).real();
351  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 3 " << endl; }
352  return result*ssfInt ;
353  }
354 
355  // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
356  Double_t _x = omega * tau ;
357  Double_t _y = tau*dgamma/2;
358 
359  if (basisType==sinBasis) {
360  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 4th form omega = " << omega << ", tau = " << tau << endl ;
361  Double_t result(0) ;
362  if (_x==0) return result*ssfInt ;
363  if (basisSign!=Minus) result += -1*evalCerfInt(+1,-_x,tau,-umin,-umax,c).imag();
364  if (basisSign!=Plus) result += -1*evalCerfInt(-1, _x,tau, umin, umax,c).imag();
365  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 4 " << endl; }
366  return result*ssfInt ;
367  }
368 
369  // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
370  if (basisType==cosBasis) {
371  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
372  Double_t result(0) ;
373  if (basisSign!=Minus) result += evalCerfInt(+1,-_x,tau,-umin,-umax,c).real();
374  if (basisSign!=Plus) result += evalCerfInt(-1, _x,tau, umin, umax,c).real();
375  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 5 " << endl; }
376  return result*ssfInt ;
377  }
378 
379  // *** 8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasis ***
380  // *** 9th form: Convolution with exp(-|t|/tau)*sinh(dgamma*t/2), used for sinhBasis ***
381  if (basisType==coshBasis || basisType == sinhBasis) {
382  if (verboseEval()>0) {cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 8th form tau=" << tau << endl ; }
383  Double_t result(0) ;
384  int sgn = ( basisType == coshBasis ? +1 : -1 );
385  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());
386  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());
387  if (TMath::IsNaN(result)) { cxcoutE(Tracing) << "RooGaussModel::analyticalIntegral(" << GetName() << ") got nan during case 6 " << endl; }
388  return result*ssfInt ;
389  }
390 
391  // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
392  if (basisType==linBasis) {
393  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 6th form tau=" << tau << endl ;
394 
395  Double_t f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
396  Double_t f1 = std::exp(-umax*umax) - std::exp(-umin*umin);
397 
398  Double_t tmp1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
399  Double_t tmp2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
400 
401  Double_t f2 = tmp1 - tmp2;
402  Double_t f3 = xpmax*tmp1 - xpmin*tmp2;
403 
404  Double_t expc2 = std::exp(c*c);
405 
406  return -tau*( f0 +
407  (2*c/rootpi)*f1 +
408  (1 - 2*c*c)*expc2*f2 +
409  expc2*f3
410  )*ssfInt;
411  }
412 
413  // *** 7th form: Convolution with (t/tau)*(t/tau)*exp(-t/tau), used for quadBasis ***
414  if (basisType==quadBasis) {
415  if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 7th form tau=" << tau << endl ;
416 
417  Double_t f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
418 
419  Double_t tmpA1 = std::exp(-umax*umax);
420  Double_t tmpA2 = std::exp(-umin*umin);
421 
422  Double_t f1 = tmpA1 - tmpA2;
423  Double_t f2 = umax*tmpA1 - umin*tmpA2;
424 
425  Double_t tmpB1 = std::exp(-xpmax)*RooMath::erfc(-umax + c);
426  Double_t tmpB2 = std::exp(-xpmin)*RooMath::erfc(-umin + c);
427 
428  Double_t f3 = tmpB1 - tmpB2;
429  Double_t f4 = xpmax*tmpB1 - xpmin*tmpB2;
430  Double_t f5 = xpmax*xpmax*tmpB1 - xpmin*xpmin*tmpB2;
431 
432  Double_t expc2 = std::exp(c*c);
433 
434  return -tau*( 2*f0 +
435  (4*c/rootpi)*((1-c*c)*f1 + c*f2) +
436  (2*c*c*(2*c*c-1) + 2)*expc2*f3 - (4*c*c-2)*expc2*f4 + expc2*f5
437  )*ssfInt;
438  }
439  R__ASSERT(0) ;
440  return 0 ;
441 }
442 
443 ////////////////////////////////////////////////////////////////////////////////
444 /// use the approximation: erf(z) = exp(-z*z)/(std::sqrt(pi)*z)
445 /// to explicitly cancel the divergent exp(y*y) behaviour of
446 /// CWERF for z = x + i y with large negative y
447 
448 std::complex<Double_t> RooGaussModel::evalCerfApprox(Double_t _x, Double_t u, Double_t c)
449 {
450  static const Double_t rootpi= std::sqrt(std::atan2(0.,-1.));
451  const std::complex<Double_t> z(_x * c, u + c);
452  const std::complex<Double_t> zc(u + c, - _x * c);
453  const std::complex<Double_t> zsq((z.real() + z.imag()) * (z.real() - z.imag()),
454  2. * z.real() * z.imag());
455  const std::complex<Double_t> v(-zsq.real() - u*u, -zsq.imag());
456  const std::complex<Double_t> ev = std::exp(v);
457  const std::complex<Double_t> mez2zcrootpi = -std::exp(zsq)/(zc*rootpi);
458 
459  return 2. * (ev * (mez2zcrootpi + 1.));
460 }
461 
462 ////////////////////////////////////////////////////////////////////////////////
463 
464 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
465 {
466  std::complex<Double_t> diff(2., 0.);
467  if (!_asympInt) {
468  diff = evalCerf(_x,umin,c);
469  diff -= evalCerf(_x,umax,c);
470  diff += RooMath::erf(umin) - RooMath::erf(umax);
471  diff *= sign;
472  }
473  diff *= std::complex<Double_t>(1., _x);
474  diff *= tau / (1.+_x*_x);
475  return diff;
476 }
477 
478 ////////////////////////////////////////////////////////////////////////////////
479 
480 Int_t RooGaussModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
481 {
482  if (matchArgs(directVars,generateVars,x)) return 1 ;
483  return 0 ;
484 }
485 
486 ////////////////////////////////////////////////////////////////////////////////
487 
488 void RooGaussModel::generateEvent(Int_t code)
489 {
490  R__ASSERT(code==1) ;
491  Double_t xmin = x.min();
492  Double_t xmax = x.max();
493  TRandom *generator = RooRandom::randomGenerator();
494  while(true) {
495  Double_t xgen = generator->Gaus(mean*msf,sigma*ssf);
496  if (xgen<xmax && xgen>xmin) {
497  x = xgen ;
498  return ;
499  }
500  }
501 }