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