ROOT   Reference Guide
RooGExpModel.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 *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17/** \class RooGExpModel
18 \ingroup Roofit
19
20The RooGExpModel is a RooResolutionModel implementation that models
21a resolution function that is the convolution of a Gaussian with
22a one-sided exponential. Such objects can be used
23for analytical convolutions with classes inheriting from RooAbsAnaConvPdf.
24\f[
25 \mathrm{GExp} = \exp \left( -\frac{1}{2} \left(\frac{x-\mu}{\sigma} \right)^2 \right)^2
26 \otimes \exp\left( -\frac{x}{\tau} \right)
27\f]
28
29**/
30
31#include "RooGExpModel.h"
32
33#include "RooFit.h"
34#include "RooMath.h"
35#include "RooRealConstant.h"
36#include "RooRandom.h"
37#include "TMath.h"
38
39
40using namespace std;
41
43
44////////////////////////////////////////////////////////////////////////////////
45/// Create a Gauss (x) Exp model with mean, sigma and tau parameters and scale factors for each parameter.
46///
47/// \note If scale factors for the parameters are not needed, RooConst(1.) can be passed.
48///
49/// \param[in] name Name of this instance.
50/// \param[in] title Title (e.g. for plotting)
51/// \param[in] x The convolution observable.
52/// \param[in] mean The mean of the Gaussian.
53/// \param[in] sigma Width of the Gaussian.
54/// \param[in] rlife Lifetime constant \f$\tau \f$.
55/// \param[in] meanSF Scale factor for mean.
56/// \param[in] sigmaSF Scale factor for sigma.
57/// \param[in] rlifeSF Scale factor for rlife.
58/// \param[in] nlo Include next-to-leading order for higher accuracy of convolution.
59/// \param[in] type Switch between normal and flipped model.
60RooGExpModel::RooGExpModel(const char *name, const char *title, RooAbsRealLValue& xIn,
61 RooAbsReal& meanIn, RooAbsReal& sigmaIn, RooAbsReal& rlifeIn,
62 RooAbsReal& meanSF, RooAbsReal& sigmaSF, RooAbsReal& rlifeSF,
63 Bool_t nlo, Type type) :
64 RooResolutionModel(name, title, xIn),
65 _mean("mean", "Mean of Gaussian component", this, meanIn),
66 sigma("sigma", "Width", this, sigmaIn),
67 rlife("rlife", "Life time", this, rlifeIn),
68 _meanSF("meanSF", "Scale factor for mean", this, meanSF),
69 ssf("ssf", "Sigma Scale Factor", this, sigmaSF),
70 rsf("rsf", "RLife Scale Factor", this, rlifeSF),
71 _flip(type==Flipped),
72 _nlo(nlo),
73 _flatSFInt(false),
74 _asympInt(false)
75{
76}
77
78
79////////////////////////////////////////////////////////////////////////////////
80/// Create a Gauss (x) Exp model with sigma and tau parameters.
81///
82/// \param[in] name Name of this instance.
83/// \param[in] title Title (e.g. for plotting)
84/// \param[in] x The convolution observable.
85/// \param[in] sigma Width of the Gaussian.
86/// \param[in] rlife Lifetime constant \f$\tau \f$.
87/// \param[in] nlo Include next-to-leading order for higher accuracy of convolution.
88/// \param[in] type Switch between normal and flipped model.
89RooGExpModel::RooGExpModel(const char *name, const char *title, RooAbsRealLValue& xIn,
90 RooAbsReal& _sigma, RooAbsReal& _rlife,
91 Bool_t nlo, Type type) :
92 RooResolutionModel(name,title,xIn),
93 _mean("mean", "Mean of Gaussian component", this, RooRealConstant::value(0.)),
94 sigma("sigma","Width",this,_sigma),
95 rlife("rlife","Life time",this,_rlife),
96 _meanSF("meanSF", "Scale factor for mean", this, RooRealConstant::value(1)),
97 ssf("ssf","Sigma Scale Factor",this,(RooRealVar&)RooRealConstant::value(1)),
98 rsf("rsf","RLife Scale Factor",this,(RooRealVar&)RooRealConstant::value(1)),
99 _flip(type==Flipped),_nlo(nlo), _flatSFInt(kFALSE), _asympInt(kFALSE)
100{
101}
102
103////////////////////////////////////////////////////////////////////////////////
104/// Create a Gauss (x) Exp model with sigma and tau parameters.
105///
106/// \param[in] name Name of this instance.
107/// \param[in] title Title (e.g. for plotting)
108/// \param[in] x The convolution observable.
109/// \param[in] sigma Width of the Gaussian.
110/// \param[in] rlife Lifetime constant \f$\tau \f$.
111/// \param[in] srSF Scale factor for both sigma and tau.
112/// \param[in] nlo Include next-to-leading order for higher accuracy of convolution.
113/// \param[in] type Switch between normal and flipped model.
114RooGExpModel::RooGExpModel(const char *name, const char *title, RooAbsRealLValue& xIn,
115 RooAbsReal& _sigma, RooAbsReal& _rlife,
117 Bool_t nlo, Type type) :
118 RooResolutionModel(name,title,xIn),
119 _mean("mean", "Mean of Gaussian component", this, RooRealConstant::value(0.)),
120 sigma("sigma","Width",this,_sigma),
121 rlife("rlife","Life time",this,_rlife),
122 _meanSF("meanSF", "Scale factor for mean", this, RooRealConstant::value(1)),
125 _flip(type==Flipped),
126 _nlo(nlo),
127 _flatSFInt(kFALSE),
128 _asympInt(kFALSE)
129{
130}
131
132////////////////////////////////////////////////////////////////////////////////
133/// Create a Gauss (x) Exp model with sigma and tau parameters and scale factors.
134///
135/// \param[in] name Name of this instance.
136/// \param[in] title Title (e.g. for plotting)
137/// \param[in] x The convolution observable.
138/// \param[in] sigma Width of the Gaussian.
139/// \param[in] rlife Lifetime constant \f$\tau \f$.
140/// \param[in] sigmaSF Scale factor for sigma.
141/// \param[in] rlifeSF Scale factor for rlife.
142/// \param[in] nlo Include next-to-leading order for higher accuracy of convolution.
143/// \param[in] type Switch between normal and flipped model.
144RooGExpModel::RooGExpModel(const char *name, const char *title, RooAbsRealLValue& xIn,
145 RooAbsReal& _sigma, RooAbsReal& _rlife,
146 RooAbsReal& _sigmaSF, RooAbsReal& _rlifeSF,
147 Bool_t nlo, Type type) :
148 RooResolutionModel(name,title,xIn),
149 _mean("mean", "Mean of Gaussian component", this, RooRealConstant::value(0.)),
150 sigma("sigma","Width",this,_sigma),
151 rlife("rlife","Life time",this,_rlife),
152 _meanSF("meanSF", "Scale factor for mean", this, RooRealConstant::value(1)),
153 ssf("ssf","Sigma Scale Factor",this,_sigmaSF),
154 rsf("rsf","RLife Scale Factor",this,_rlifeSF),
155 _flip(type==Flipped),
156 _nlo(nlo),
157 _flatSFInt(kFALSE),
158 _asympInt(kFALSE)
159{
160}
161
162////////////////////////////////////////////////////////////////////////////////
163
164RooGExpModel::RooGExpModel(const RooGExpModel& other, const char* name) :
166 _mean("mean", this, other._mean),
167 sigma("sigma",this,other.sigma),
168 rlife("rlife",this,other.rlife),
169 _meanSF("meanSf", this, other._meanSF),
170 ssf("ssf",this,other.ssf),
171 rsf("rsf",this,other.rsf),
172 _flip(other._flip),
173 _nlo(other._nlo),
174 _flatSFInt(other._flatSFInt),
175 _asympInt(other._asympInt)
176{
177}
178
179////////////////////////////////////////////////////////////////////////////////
180/// Destructor
181
183{
184}
185
186////////////////////////////////////////////////////////////////////////////////
187
189{
190 if (!TString("exp(-@0/@1)").CompareTo(name)) return expBasisPlus ;
191 if (!TString("exp(@0/@1)").CompareTo(name)) return expBasisMinus ;
192 if (!TString("exp(-abs(@0)/@1)").CompareTo(name)) return expBasisSum ;
193 if (!TString("exp(-@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisPlus ;
194 if (!TString("exp(@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisMinus ;
195 if (!TString("exp(-abs(@0)/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisSum ;
196 if (!TString("exp(-@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisPlus ;
197 if (!TString("exp(@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisMinus ;
198 if (!TString("exp(-abs(@0)/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisSum ;
199 if (!TString("exp(-@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisPlus;
200 if (!TString("exp(@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisMinus;
201 if (!TString("exp(-abs(@0)/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisSum;
202 if (!TString("exp(-@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisPlus;
203 if (!TString("exp(@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisMinus;
204 if (!TString("exp(-abs(@0)/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisSum;
205 return 0 ;
206}
207
208
209namespace {
210////////////////////////////////////////////////////////////////////////////////
211/// Approximation of the log of the complex error function
212Double_t logErfC(Double_t xx)
213{
214 Double_t t,z,ans;
215 z=fabs(xx);
216 t=1.0/(1.0+0.5*z);
217
218 if(xx >= 0.0)
219 ans=log(t)+(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+
220 t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))));
221 else
222 ans=log(2.0-t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+
223 t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277))))))))));
224
225 return ans;
226}
227
228////////////////////////////////////////////////////////////////////////////////
229/// use the approximation: erf(z) = exp(-z*z)/(sqrt(pi)*z)
230/// to explicitly cancel the divergent exp(y*y) behaviour of
231/// CWERF for z = x + i y with large negative y
232
233std::complex<Double_t> evalCerfApprox(Double_t swt, Double_t u, Double_t c)
234{
235 static Double_t rootpi= sqrt(atan2(0.,-1.));
236 std::complex<Double_t> z(swt*c,u+c);
237 std::complex<Double_t> zc(u+c,-swt*c);
238 std::complex<Double_t> zsq= z*z;
239 std::complex<Double_t> v= -zsq - u*u;
240
241 return std::exp(v)*(-std::exp(zsq)/(zc*rootpi) + 1.)*2.;
242}
243
244
245// Calculate exp(-u^2) cwerf(swt*c + i(u+c)), taking care of numerical instabilities
246std::complex<Double_t> evalCerf(Double_t swt, Double_t u, Double_t c)
247{
248 std::complex<Double_t> z(swt*c,u+c);
249 return (z.imag()>-4.0) ? RooMath::faddeeva_fast(z)*std::exp(-u*u) : evalCerfApprox(swt,u,c) ;
250}
251
252
253// Calculate Re(exp(-u^2) cwerf(i(u+c)))
255inline Double_t evalCerfRe(Double_t u, Double_t c) {
256 Double_t expArg = u*2*c+c*c ;
257 if (expArg<300) {
258 return exp(expArg) * RooMath::erfc(u+c);
259 } else {
260 return exp(expArg+logErfC(u+c));
261 }
262}
263
264}
265
266
267
268////////////////////////////////////////////////////////////////////////////////
269
271{
272 static Double_t root2(sqrt(2.)) ;
273
274 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
275 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
276
277 Double_t fsign = _flip?-1:1 ;
278
279 Double_t sig = sigma*ssf ;
280 Double_t rtau = rlife*rsf ;
281
284 if (basisType == coshBasis && _basisCode!=noBasis ) {
285 Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
286 if (dGamma==0) basisType = expBasis;
287 }
288
289 // *** 1st form: Straight GExp, used for unconvoluted PDF or expBasis with 0 lifetime ***
290 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
291 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName() << ") 1st form" << endl ;
292
293 Double_t expArg = sig*sig/(2*rtau*rtau) + fsign*(x - _mean*_meanSF)/rtau ;
294
295 Double_t result ;
296 if (expArg<300) {
297 result = 1/(2*rtau) * exp(expArg) * RooMath::erfc(sig/(root2*rtau) + fsign*(x - _mean*_meanSF)/(root2*sig));
298 } else {
299 // If exponent argument is very large, bring canceling RooMath::erfc() term inside exponent
300 // to avoid floating point over/underflows of intermediate calculations
301 result = 1/(2*rtau) * exp(expArg + logErfC(sig/(root2*rtau) + fsign*(x - _mean*_meanSF)/(root2*sig))) ;
302 }
303
304// Double_t result = 1/(2*rtau)
305// * exp(sig*sig/(2*rtau*rtau) + fsign*x/rtau)
306// * RooMath::erfc(sig/(root2*rtau) + fsign*x/(root2*sig));
307
308 // equivalent form, added FMV, 07/24/03
309 //Double_t xprime = x/rtau ;
310 //Double_t c = sig/(root2*rtau) ;
311 //Double_t u = xprime/(2*c) ;
312 //Double_t result = 0.5*evalCerf(fsign*u,c).real() ; // sign=-1 !
313
314 if (_basisCode!=0 && basisSign==Both) result *= 2 ;
315 //cout << "1st form " << "x= " << x << " result= " << result << endl;
316 return result ;
317 }
318
319 // *** 2nd form: 0, used for sinBasis and cosBasis with tau=0 ***
320 if (tau==0) {
321 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName() << ") 2nd form" << endl ;
322 return 0. ;
323 }
324
325 Double_t omega = (basisType!=expBasis)?((RooAbsReal*)basis().getParameter(2))->getVal():0. ;
326
327 // *** 3nd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
328 if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
329 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName() << ") 3d form tau=" << tau << endl ;
330 Double_t result(0) ;
331 if (basisSign!=Minus) result += calcDecayConv(+1,tau,sig,rtau,fsign) ; // modified FMV,08/13/03
332 if (basisSign!=Plus) result += calcDecayConv(-1,tau,sig,rtau,fsign) ; // modified FMV,08/13/03
333 //cout << "3rd form " << "x= " << x << " result= " << result << endl;
334 return result ;
335 }
336
337 // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
338 Double_t wt = omega *tau ;
339 if (basisType==sinBasis) {
340 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName() << ") 4th form omega = "
341 << omega << ", tau = " << tau << endl ;
342 Double_t result(0) ;
343 if (wt==0.) return result ;
344 if (basisSign!=Minus) result += -1*calcSinConv(+1,sig,tau,omega,rtau,fsign).imag() ;
345 if (basisSign!=Plus) result += -1*calcSinConv(-1,sig,tau,omega,rtau,fsign).imag() ;
346 //cout << "4th form " << "x= " << x << " result= " << result << endl;
347 return result ;
348 }
349
350 // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
351 if (basisType==cosBasis) {
352 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName()
353 << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
354 Double_t result(0) ;
355 if (basisSign!=Minus) result += calcSinConv(+1,sig,tau,omega,rtau,fsign).real() ;
356 if (basisSign!=Plus) result += calcSinConv(-1,sig,tau,omega,rtau,fsign).real() ;
357 //cout << "5th form " << "x= " << x << " result= " << result << endl;
358 return result ;
359 }
360
361
362 // *** 6th form: Convolution with exp(-t/tau)*sinh(dgamma*t/2), used for sinhBasis ***
363 if (basisType==sinhBasis) {
364 Double_t dgamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
365
366 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName()
367 << ") 6th form = " << dgamma << ", tau = " << tau << endl;
368 Double_t result(0);
369 //if (basisSign!=Minus) result += calcSinhConv(+1,+1,-1,tau,dgamma,sig,rtau,fsign);
370 //if (basisSign!=Plus) result += calcSinhConv(-1,-1,+1,tau,dgamma,sig,rtau,fsign);
371 // better form, since it also accounts for the numerical divergence region, added FMV, 07/24/03
372 Double_t tau1 = 1/(1/tau-dgamma/2) ;
373 Double_t tau2 = 1/(1/tau+dgamma/2) ;
374 if (basisSign!=Minus) result += 0.5*(calcDecayConv(+1,tau1,sig,rtau,fsign)-calcDecayConv(+1,tau2,sig,rtau,fsign));
375 // modified FMV,08/13/03
376 if (basisSign!=Plus) result += 0.5*(calcDecayConv(-1,tau2,sig,rtau,fsign)-calcDecayConv(-1,tau1,sig,rtau,fsign));
377 // modified FMV,08/13/03
378 //cout << "6th form " << "x= " << x << " result= " << result << endl;
379 return result;
380 }
381
382 // *** 7th form: Convolution with exp(-t/tau)*cosh(dgamma*t/2), used for coshBasis ***
383 if (basisType==coshBasis) {
384 Double_t dgamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
385
386 if (verboseEval()>2) cout << "RooGExpModel::evaluate(" << GetName()
387 << ") 7th form = " << dgamma << ", tau = " << tau << endl;
388 Double_t result(0);
389 //if (basisSign!=Minus) result += calcCoshConv(+1,tau,dgamma,sig,rtau,fsign);
390 //if (basisSign!=Plus) result += calcCoshConv(-1,tau,dgamma,sig,rtau,fsign);
391 // better form, since it also accounts for the numerical divergence region, added FMV, 07/24/03
392 Double_t tau1 = 1/(1/tau-dgamma/2) ;
393 Double_t tau2 = 1/(1/tau+dgamma/2) ;
394 if (basisSign!=Minus) result += 0.5*(calcDecayConv(+1,tau1,sig,rtau,fsign)+calcDecayConv(+1,tau2,sig,rtau,fsign));
395 // modified FMV,08/13/03
396 if (basisSign!=Plus) result += 0.5*(calcDecayConv(-1,tau1,sig,rtau,fsign)+calcDecayConv(-1,tau2,sig,rtau,fsign));
397 // modified FMV,08/13/03
398 //cout << "7th form " << "x= " << x << " result= " << result << endl;
399 return result;
400 }
401 R__ASSERT(0) ;
402 return 0 ;
403 }
404
405
406////////////////////////////////////////////////////////////////////////////////
407
408std::complex<Double_t> RooGExpModel::calcSinConv(Double_t sign, Double_t sig, Double_t tau, Double_t omega, Double_t rtau, Double_t fsign) const
409{
410 static Double_t root2(sqrt(2.)) ;
411
412 Double_t s1= -sign*(x - _mean*_meanSF)/tau;
413 //Double_t s1= x/tau;
414 Double_t c1= sig/(root2*tau);
415 Double_t u1= s1/(2*c1);
416 Double_t s2= (x - _mean*_meanSF)/rtau;
417 Double_t c2= sig/(root2*rtau);
418 Double_t u2= fsign*s2/(2*c2) ;
419 //Double_t u2= s2/(2*c2) ;
420
421 std::complex<Double_t> eins(1,0);
422 std::complex<Double_t> k(1/tau,sign*omega);
423 //return (evalCerf(-sign*omega*tau,u1,c1)+evalCerf(0,u2,c2)*fsign*sign) / (eins + k*fsign*sign*rtau) ;
424
425 return (evalCerf(-sign*omega*tau,u1,c1)+std::complex<Double_t>(evalCerfRe(u2,c2),0)*fsign*sign) / (eins + k*fsign*sign*rtau) ;
426 // equivalent form, added FMV, 07/24/03
427 //return (evalCerf(-sign*omega*tau,-sign*u1,c1)+evalCerf(0,fsign*u2,c2)*fsign*sign) / (eins + k*fsign*sign*rtau) ;
428}
429
431
432////////////////////////////////////////////////////////////////////////////////
433
435{
436 static Double_t root2(sqrt(2.)) ;
437
438 Double_t s1= -sign*(x - _mean*_meanSF)/tau;
439 //Double_t s1= x/tau;
440 Double_t c1= sig/(root2*tau);
441 Double_t u1= s1/(2*c1);
442 Double_t s2= (x - _mean*_meanSF)/rtau;
443 Double_t c2= sig/(root2*rtau);
444 Double_t u2= fsign*s2/(2*c2) ;
445 //Double_t u2= s2/(2*c2) ;
446
447 Double_t eins(1);
448 Double_t k(1/tau);
449 return (evalCerfRe(u1,c1)+evalCerfRe(u2,c2)*fsign*sign) / (eins + k*fsign*sign*rtau) ;
450 // equivalent form, added FMV, 07/24/03
451 //return (evalCerf(-sign*u1,c1).real()+evalCerf(fsign*u2,c2).real()*fsign*sign) / (eins + k*fsign*sign*rtau) ;
452}
453
454////////////////////////////////////////////////////////////////////////////////
455
457// modified FMV,08/13/03
458{
459 static Double_t root2(sqrt(2.)) ;
460 static Double_t root2pi(sqrt(2*atan2(0.,-1.))) ;
461 static Double_t rootpi(sqrt(atan2(0.,-1.)));
462
463 // Process flip status
464 Double_t xp(x - _mean*_meanSF) ;
465 //if (_flip) {
466 // xp *= -1 ;
467 // sign *= -1 ;
468 //}
469 xp *= fsign ; // modified FMV,08/13/03
470 sign *= fsign ; // modified FMV,08/13/03
471
472 Double_t cFly;
473 if ((sign<0)&&(fabs(tau-rtau)<tau/260)) {
474
475 Double_t MeanTau=0.5*(tau+rtau);
476 if (fabs(xp/MeanTau)>300) {
477 return 0 ;
478 }
479
480 cFly=1./(MeanTau*MeanTau*root2pi) *
481 exp(-(-xp/MeanTau-sig*sig/(2*MeanTau*MeanTau)))
482 *(sig*exp(-1/(2*sig*sig)*TMath::Power((sig*sig/MeanTau+xp),2))
483 -(sig*sig/MeanTau+xp)*(rootpi/root2)*RooMath::erfc(sig/(root2*MeanTau)+xp/(root2*sig)));
484
485 if(_nlo) {
486 Double_t epsilon=0.5*(tau-rtau);
487 Double_t a=sig/(root2*MeanTau)+xp/(root2*sig);
488 cFly += 1./(MeanTau*MeanTau)
489 *exp(-(-xp/MeanTau-sig*sig/(2*MeanTau*MeanTau)))
490 *0.5/MeanTau*epsilon*epsilon*
491 (exp(-a*a)*(sig/MeanTau*root2/rootpi
492 -(4*a*sig*sig)/(2*rootpi*MeanTau*MeanTau)
493 +(-4/rootpi+8*a*a/rootpi)/6
494 *TMath::Power(sig/(root2*MeanTau),3)
495 +2/rootpi*(sig*sig/(MeanTau*MeanTau)+xp/MeanTau)*
496 (sig/(root2*MeanTau)-a*(sig*sig)/(2*MeanTau*MeanTau))
497 +2/rootpi*((3*sig*sig)/(2*MeanTau*MeanTau)+xp/MeanTau+
498 0.5*TMath::Power(sig*sig/(MeanTau*MeanTau)+xp/MeanTau,2))*sig/(root2*MeanTau))
499 -(2*sig*sig/(MeanTau*MeanTau)+xp/MeanTau+(sig*sig/(MeanTau*MeanTau)+xp/MeanTau)*
500 (3*sig*sig/(2*MeanTau*MeanTau)+xp/MeanTau)
501 +TMath::Power(sig*sig/(MeanTau*MeanTau)+xp/MeanTau,3)/6)*RooMath::erfc(a));
502 }
503
504 } else {
505
506 Double_t expArg1 = sig*sig/(2*tau*tau)-sign*xp/tau ;
507 Double_t expArg2 = sig*sig/(2*rtau*rtau)+xp/rtau ;
508
509 Double_t term1, term2 ;
510 if (expArg1<300) {
511 term1 = exp(expArg1) *RooMath::erfc(sig/(root2*tau)-sign*xp/(root2*sig)) ;
512 } else {
513 term1 = exp(expArg1+logErfC(sig/(root2*tau)-sign*xp/(root2*sig))) ; ;
514 }
515 if (expArg2<300) {
516 term2 = exp(expArg2) *RooMath::erfc(sig/(root2*rtau)+xp/(root2*sig)) ;
517 } else {
518 term2 = exp(expArg2+logErfC(sig/(root2*rtau)+xp/(root2*sig))) ; ;
519 }
520
521 cFly=(term1+sign*term2)/(2*(tau+sign*rtau));
522
523 // WVE prevent numeric underflows
524 if (cFly<1e-100) {
525 cFly = 0 ;
526 }
527
528 // equivalent form, added FMV, 07/24/03
529 //cFly = calcSinConv(sign, sig, tau, rtau, fsign)/(2*tau) ;
530 }
531
532 return cFly*2*tau ;
533}
534
535/* commented FMV, 07/24/03
536
537////////////////////////////////////////////////////////////////////////////////
538
539Double_t RooGExpModel::calcCoshConv(Double_t sign, Double_t tau, Double_t dgamma, Double_t sig, Double_t rtau, Double_t fsign) const
540{
541
542
543 static Double_t root2(sqrt(2.)) ;
544 static Double_t root2pi(sqrt(2*atan2(0.,-1.))) ;
545 static Double_t rootpi(sqrt(atan2(0.,-1.)));
546 Double_t tau1 = 1/(1/tau-dgamma/2);
547 Double_t tau2 = 1/(1/tau+dgamma/2);
548 Double_t cFly;
549 Double_t xp(x);
550
551 //if (_flip) {
552 // xp *= -1 ;
553 // sign *= -1 ;
554 //}
555 xp *= fsign ; // modified FMV,08/13/03
556 sign *= fsign ; // modified FMV,08/13/03
557
558 cFly=tau1*(exp(sig*sig/(2*tau1*tau1)-sign*xp/tau1)
559 *RooMath::erfc(sig/(root2*tau1)-sign*xp/(root2*sig))
560 +sign*exp(sig*sig/(2*rtau*rtau)+xp/rtau)
561 *RooMath::erfc(sig/(root2*rtau)+xp/(root2*sig)))/(2*(tau1+sign*rtau))
562 +tau2*(exp(sig*sig/(2*tau2*tau2)-sign*xp/tau2)
563 *RooMath::erfc(sig/(root2*tau2)-sign*xp/(root2*sig))
564 +sign*exp(sig*sig/(2*rtau*rtau)+xp/rtau)
565 *RooMath::erfc(sig/(root2*rtau)+xp/(root2*sig)))/(2*(tau2+sign*rtau));;
566 return cFly;
567}
568*/
569
570/* commented FMV, 07/24/03
571
572////////////////////////////////////////////////////////////////////////////////
573
574Double_t RooGExpModel::calcSinhConv(Double_t sign, Double_t sign1, Double_t sign2, Double_t tau, Double_t dgamma, Double_t sig, Double_t rtau, Double_t fsign) const
575{
576 static Double_t root2(sqrt(2.)) ;
577 static Double_t root2pi(sqrt(2*atan2(0.,-1.))) ;
578 static Double_t rootpi(sqrt(atan2(0.,-1.)));
579 Double_t tau1 = 1/(1/tau-dgamma/2);
580 Double_t tau2 = 1/(1/tau+dgamma/2);
581 Double_t cFly;
582 Double_t xp(x);
583
584 //if (_flip) {
585 // xp *= -1 ;
586 // sign1 *= -1 ;
587 // sign2 *= -1 ;
588 //}
589 xp *= fsign ; // modified FMV,08/13/03
590 sign1 *= fsign ; // modified FMV,08/13/03
591 sign2 *= fsign ; // modified FMV,08/13/03
592
593 cFly=sign1*tau1*(exp(sig*sig/(2*tau1*tau1)-sign*xp/tau1)
594 *RooMath::erfc(sig/(root2*tau1)-sign*xp/(root2*sig))
595 +sign*exp(sig*sig/(2*rtau*rtau)+xp/rtau)
596 *RooMath::erfc(sig/(root2*rtau)+xp/(root2*sig)))/(2*(tau1+sign*rtau))
597 +sign2*tau2*(exp(sig*sig/(2*tau2*tau2)-sign*xp/tau2)
598 *RooMath::erfc(sig/(root2*tau2)-sign*xp/(root2*sig))
599 +sign*exp(sig*sig/(2*rtau*rtau)+xp/rtau)
600 *RooMath::erfc(sig/(root2*rtau)+xp/(root2*sig)))/(2*(tau2+sign*rtau));;
601 return cFly;
602}
603*/
604
605////////////////////////////////////////////////////////////////////////////////
606
607Int_t RooGExpModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
608{
609 switch(_basisCode) {
610
611 // Analytical integration capability of raw PDF
612 case noBasis:
613 if (matchArgs(allVars,analVars,convVar())) return 1 ;
614 break ;
615
616 // Analytical integration capability of convoluted PDF
617 case expBasisPlus:
618 case expBasisMinus:
619 case expBasisSum:
620 case sinBasisPlus:
621 case sinBasisMinus:
622 case sinBasisSum:
623 case cosBasisPlus:
624 case cosBasisMinus:
625 case cosBasisSum:
626 case sinhBasisPlus:
627 case sinhBasisMinus:
628 case sinhBasisSum:
629 case coshBasisPlus:
630 case coshBasisMinus:
631 case coshBasisSum:
632
633 // Optionally advertise flat integral over sigma scale factor
634 if (_flatSFInt) {
635 if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) {
636 return 2 ;
637 }
638 }
639
640 if (matchArgs(allVars,analVars,convVar())) return 1 ;
641 break ;
642 }
643
644 return 0 ;
645}
646
647////////////////////////////////////////////////////////////////////////////////
648
649Double_t RooGExpModel::analyticalIntegral(Int_t code, const char* rangeName) const
650{
651 static Double_t root2 = sqrt(2.) ;
652// static Double_t rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
653 Double_t ssfInt(1.0) ;
654
655 // Code must be 1 or 2
656 R__ASSERT(code==1||code==2) ;
657 if (code==2) {
658 ssfInt = (ssf.max(rangeName)-ssf.min(rangeName)) ;
659 }
660
661 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
662 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
663
665
667 if (basisType == coshBasis && _basisCode!=noBasis ) {
668 Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
669 if (dGamma==0) basisType = expBasis;
670 }
671 Double_t fsign = _flip?-1:1 ;
672 Double_t sig = sigma*ssf ;
673 Double_t rtau = rlife*rsf ;
674
675 // *** 1st form????
676 if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
677 if (verboseEval()>0) cout << "RooGExpModel::analyticalIntegral(" << GetName() << ") 1st form" << endl ;
678
679 //Double_t result = 1.0 ; // WVE inferred from limit(tau->0) of cosBasisNorm
680 // finite+asymtotic normalization, added FMV, 07/24/03
681 Double_t xpmin = (x.min(rangeName) - _mean*_meanSF)/rtau ;
682 Double_t xpmax = (x.max(rangeName) - _mean*_meanSF)/rtau ;
683 Double_t c = sig/(root2*rtau) ;
684 Double_t umin = xpmin/(2*c) ;
685 Double_t umax = xpmax/(2*c) ;
686 Double_t result ;
687 if (_asympInt) {
688 result = 1.0 ;
689 } else {
690 result = 0.5*evalCerfInt(-fsign,rtau,fsign*umin,fsign*umax,c)/rtau ; //WVEFIX add 1/rtau
691 }
692
693 if (_basisCode!=0 && basisSign==Both) result *= 2 ;
694 //cout << "Integral 1st form " << " result= " << result*ssfInt << endl;
695 return result*ssfInt ;
696 }
697
698 Double_t omega = (basisType!=expBasis) ?((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
699
700 // *** 2nd form: unity, used for sinBasis and cosBasis with tau=0 (PDF is zero) ***
701 //if (tau==0&&omega!=0) {
702 if (tau==0) { // modified, FMV 07/24/03
703 if (verboseEval()>0) cout << "RooGExpModel::analyticalIntegral(" << GetName() << ") 2nd form" << endl ;
704 return 0. ;
705 }
706
707 // *** 3rd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
708 if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
709 //Double_t result = 2*tau ;
710 //if (basisSign==Both) result *= 2 ;
711 // finite+asymtotic normalization, added FMV, 07/24/03
712 Double_t result(0.);
713 if (basisSign!=Minus) result += calcSinConvNorm(+1,tau,sig,rtau,fsign,rangeName);
714 if (basisSign!=Plus) result += calcSinConvNorm(-1,tau,sig,rtau,fsign,rangeName);
715 //cout << "Integral 3rd form " << " result= " << result*ssfInt << endl;
716 return result*ssfInt ;
717 }
718
719 // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
720 Double_t wt = omega * tau ;
721 if (basisType==sinBasis) {
722 if (verboseEval()>0) cout << "RooGExpModel::analyticalIntegral(" << GetName() << ") 4th form omega = "
723 << omega << ", tau = " << tau << endl ;
724 //cout << "sin integral" << endl;
725 Double_t result(0) ;
726 if (wt==0) return result ;
727 //if (basisSign!=Minus) result += calcSinConvNorm(+1,tau,omega).imag() ;
728 //if (basisSign!=Plus) result += calcSinConvNorm(-1,tau,omega).imag() ;
729 // finite+asymtotic normalization, added FMV, 07/24/03
730 if (basisSign!=Minus) result += -1*calcSinConvNorm(+1,tau,omega,sig,rtau,fsign,rangeName).imag();
731 if (basisSign!=Plus) result += -1*calcSinConvNorm(-1,tau,omega,sig,rtau,fsign,rangeName).imag();
732 //cout << "Integral 4th form " << " result= " << result*ssfInt << endl;
733 return result*ssfInt ;
734 }
735
736 // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
737 if (basisType==cosBasis) {
738 if (verboseEval()>0) cout << "RooGExpModel::analyticalIntegral(" << GetName()
739 << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
740 //cout << "cos integral" << endl;
741 Double_t result(0) ;
742 //if (basisSign!=Minus) result += calcSinConvNorm(+1,tau,omega).real() ;
743 //if (basisSign!=Plus) result += calcSinConvNorm(-1,tau,omega).real() ;
744 // finite+asymtotic normalization, added FMV, 07/24/03
745 if (basisSign!=Minus) result += calcSinConvNorm(+1,tau,omega,sig,rtau,fsign,rangeName).real();
746 if (basisSign!=Plus) result += calcSinConvNorm(-1,tau,omega,sig,rtau,fsign,rangeName).real();
747 //cout << "Integral 5th form " << " result= " << result*ssfInt << endl;
748 return result*ssfInt ;
749 }
750
751 Double_t dgamma = ((basisType==coshBasis)||(basisType==sinhBasis))?((RooAbsReal*)basis().getParameter(2))->getVal():0 ;
752
753 // *** 6th form: Convolution with exp(-t/tau)*sinh(dgamma*t/2), used for sinhBasis ***
754 if (basisType==sinhBasis) {
755 if (verboseEval()>0) cout << "RooGExpModel::analyticalIntegral(" << GetName()
756 << ") 6th form dgamma = " << dgamma << ", tau = " << tau << endl ;
757 Double_t tau1 = 1/(1/tau-dgamma/2);
758 Double_t tau2 = 1/(1/tau+dgamma/2);
759 //cout << "sinh integral" << endl;
760 Double_t result(0) ;
761 //if (basisSign!=Minus) result += tau1-tau2 ;
762 //if (basisSign!=Plus) result += tau2-tau1 ;
763 // finite+asymtotic normalization, added FMV, 07/24/03
764 if (basisSign!=Minus) result += 0.5*(calcSinConvNorm(+1,tau1,sig,rtau,fsign,rangeName)-
765 calcSinConvNorm(+1,tau2,sig,rtau,fsign,rangeName));
766 if (basisSign!=Plus) result += 0.5*(calcSinConvNorm(-1,tau2,sig,rtau,fsign,rangeName)-
767 calcSinConvNorm(-1,tau1,sig,rtau,fsign,rangeName));
768 //cout << "Integral 6th form " << " result= " << result*ssfInt << endl;
769 return result;
770 }
771
772 // ** 7th form: Convolution with exp(-t/tau)*cosh(dgamma*t/2), used for coshBasis ***
773 if (basisType==coshBasis) {
774 if (verboseEval()>0) cout << "RooGExpModel::analyticalIntegral(" << GetName()
775 << ") 6th form dgamma = " << dgamma << ", tau = " << tau << endl ;
776 //cout << "cosh integral" << endl;
777 Double_t tau1 = 1/(1/tau-dgamma/2);
778 Double_t tau2 = 1/(1/tau+dgamma/2);
779 //Double_t result = (tau1+tau2) ;
780 //if (basisSign==Both) result *= 2 ;
781 // finite+asymtotic normalization, added FMV, 07/24/03
782 Double_t result(0);
783 if (basisSign!=Minus) result += 0.5*(calcSinConvNorm(+1,tau1,sig,rtau,fsign,rangeName)+
784 calcSinConvNorm(+1,tau2,sig,rtau,fsign,rangeName));
785 if (basisSign!=Plus) result += 0.5*(calcSinConvNorm(-1,tau1,sig,rtau,fsign,rangeName)+
786 calcSinConvNorm(-1,tau2,sig,rtau,fsign,rangeName));
787 //cout << "Integral 7th form " << " result= " << result*ssfInt << endl;
788 return result;
789
790 }
791
792 R__ASSERT(0) ;
793 return 1 ;
794}
795
796// modified FMV, 07/24/03. Finite+asymtotic normalization
797
798////////////////////////////////////////////////////////////////////////////////
799/// old code (asymptotic normalization only)
800/// std::complex<Double_t> z(1/tau,sign*omega);
801/// return z*2/(omega*omega+1/(tau*tau));
802
803std::complex<Double_t> RooGExpModel::calcSinConvNorm(Double_t sign, Double_t tau, Double_t omega,
804 Double_t sig, Double_t rtau, Double_t fsign, const char* rangeName) const
805{
806 static Double_t root2(sqrt(2.)) ;
807
808 Double_t smin1= (x.min(rangeName) - _mean*_meanSF)/tau;
809 Double_t smax1= (x.max(rangeName) - _mean*_meanSF)/tau;
810 Double_t c1= sig/(root2*tau);
811 Double_t umin1= smin1/(2*c1);
812 Double_t umax1= smax1/(2*c1);
813 Double_t smin2= (x.min(rangeName) - _mean*_meanSF)/rtau;
814 Double_t smax2= (x.max(rangeName) - _mean*_meanSF)/rtau;
815 Double_t c2= sig/(root2*rtau);
816 Double_t umin2= smin2/(2*c2) ;
817 Double_t umax2= smax2/(2*c2) ;
818
819 std::complex<Double_t> eins(1,0);
820 std::complex<Double_t> k(1/tau,sign*omega);
821 std::complex<Double_t> term1 = evalCerfInt(sign,-sign*omega*tau, tau, -sign*umin1, -sign*umax1, c1);
822 //std::complex<Double_t> term2 = evalCerfInt(-fsign,0., rtau, fsign*umin2, fsign*umax2, c2)*std::complex<Double_t>(fsign*sign,0);
823 std::complex<Double_t> term2 = std::complex<Double_t>(evalCerfInt(-fsign, rtau, fsign*umin2, fsign*umax2, c2)*fsign*sign,0);
824 return (term1+term2)/(eins + k*fsign*sign*rtau) ;
825}
826
828
829////////////////////////////////////////////////////////////////////////////////
830
831Double_t RooGExpModel::calcSinConvNorm(Double_t sign, Double_t tau, Double_t sig, Double_t rtau, Double_t fsign, const char* rangeName) const
832{
833 static Double_t root2(sqrt(2.)) ;
834
835 Double_t smin1= (x.min(rangeName) - _mean*_meanSF)/tau;
836 Double_t smax1= (x.max(rangeName) - _mean*_meanSF)/tau;
837 Double_t c1= sig/(root2*tau);
838 Double_t umin1= smin1/(2*c1);
839 Double_t umax1= smax1/(2*c1);
840 Double_t smin2= (x.min(rangeName) - _mean*_meanSF)/rtau;
841 Double_t smax2= (x.max(rangeName) - _mean*_meanSF)/rtau;
842 Double_t c2= sig/(root2*rtau);
843 Double_t umin2= smin2/(2*c2) ;
844 Double_t umax2= smax2/(2*c2) ;
845
846 Double_t eins(1);
847 Double_t k(1/tau);
848 Double_t term1 = evalCerfInt(sign, tau, -sign*umin1, -sign*umax1, c1);
849 Double_t term2 = evalCerfInt(-fsign, rtau, fsign*umin2, fsign*umax2, c2)*fsign*sign;
850
851 // WVE Handle 0/0 numeric divergence
852 if (fabs(tau-rtau)<1e-10 && fabs(term1+term2)<1e-10) {
853 cout << "epsilon method" << endl ;
854 static Double_t epsilon = 1e-4 ;
855 return calcSinConvNorm(sign,tau+epsilon,sig,rtau-epsilon,fsign,rangeName) ;
856 }
857 return (term1+term2)/(eins + k*fsign*sign*rtau) ;
858}
859
861////////////////////////////////////////////////////////////////////////////////
862
863std::complex<Double_t> RooGExpModel::evalCerfInt(Double_t sign, Double_t wt, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
864{
865 std::complex<Double_t> diff;
866 if (_asympInt) {
867 diff = std::complex<Double_t>(2,0) ;
868 } else {
869 diff = std::complex<Double_t>(sign,0.)*(evalCerf(wt,umin,c) - evalCerf(wt,umax,c) + RooMath::erf(umin) - RooMath::erf(umax));
870 }
871 return std::complex<Double_t>(tau/(1.+wt*wt),0)*std::complex<Double_t>(1,wt)*diff;
872}
873// added FMV, 08/17/03. Modified FMV, 08/30/03
874
875////////////////////////////////////////////////////////////////////////////////
876
878{
879 Double_t diff;
880 if (_asympInt) {
881 diff = 2. ;
882 } else {
883 if ((umin<-8 && umax>8)||(umax<-8 && umin>8)) {
884 // If integral is over >8 sigma, approximate with full integral
885 diff = 2. ;
886 } else {
887 diff = sign*(evalCerfRe(umin,c) - evalCerfRe(umax,c) + RooMath::erf(umin) - RooMath::erf(umax));
888 }
889 }
890 return tau*diff;
891}
892
893////////////////////////////////////////////////////////////////////////////////
894
895Int_t RooGExpModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
896{
897 if (matchArgs(directVars,generateVars,x)) return 1 ;
898 return 0 ;
899}
900
901////////////////////////////////////////////////////////////////////////////////
902
904{
905 R__ASSERT(code==1) ;
906 Double_t xgen ;
907 while (true) {
910 if (!_flip)
911 xgen = xgau + (rlife*rsf)*log(xexp); // modified, FMV 08/13/03
912 else
913 xgen = xgau - (rlife*rsf)*log(xexp);
914
915 if (xgen < (x.max() - _mean*_meanSF) && xgen > (x.min() - _mean*_meanSF)) {
916 x = xgen + _mean*_meanSF;
917 return ;
918 }
919 }
920}
#define c(i)
Definition: RSha256.hxx:101
#define s1(x)
Definition: RSha256.hxx:91
#define e(i)
Definition: RSha256.hxx:103
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
double atan2(double, double)
double sqrt(double)
double exp(double)
double log(double)
friend class RooArgSet
Definition: RooAbsArg.h:572
static int verboseEval()
Return global level of verbosity for p.d.f. evaluations.
Definition: RooAbsPdf.cxx:3303
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:60
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:90
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooAbsArg * getParameter(const char *name) const
Definition: RooFormulaVar.h:40
The RooGExpModel is a RooResolutionModel implementation that models a resolution function that is the...
Definition: RooGExpModel.h:25
RooRealProxy ssf
Definition: RooGExpModel.h:108
virtual Int_t basisCode(const char *name) const
std::complex< Double_t > calcSinConvNorm(Double_t sign, Double_t tau, Double_t omega, Double_t sig, Double_t rtau, Double_t fsign, const char *rangeName) const
old code (asymptotic normalization only) std::complex<Double_t> z(1/tau,sign*omega); return z*2/(omeg...
Bool_t _flatSFInt
Definition: RooGExpModel.h:113
virtual ~RooGExpModel()
Destructor.
RooRealProxy rlife
Definition: RooGExpModel.h:106
RooRealProxy sigma
Definition: RooGExpModel.h:105
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
std::complex< Double_t > evalCerfInt(Double_t sign, Double_t wt, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
RooRealProxy rsf
Definition: RooGExpModel.h:109
void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
RooRealProxy _meanSF
Definition: RooGExpModel.h:107
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.
Bool_t _asympInt
Definition: RooGExpModel.h:114
Double_t calcDecayConv(Double_t sign, Double_t tau, Double_t sig, Double_t rtau, Double_t fsign) const
RooRealProxy _mean
Definition: RooGExpModel.h:104
std::complex< Double_t > calcSinConv(Double_t sign, Double_t sig, Double_t tau, Double_t omega, Double_t rtau, Double_t fsign) const
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
Definition: RooMath.cxx:556
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:580
static std::complex< double > faddeeva_fast(std::complex< double > z)
evaluate Faddeeva function for complex argument (fast version)
Definition: RooMath.cxx:549
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
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:35
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
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:263
Basic string class.
Definition: TString.h:131
const Double_t sigma
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
return c2
Definition: legend2.C:14
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617