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