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