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