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