Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
20 Class RooGExpModel is a RooResolutionModel implementation that models
21 a resolution function that is the convolution of a Gaussian with
22 a one-sided exponential. Object of class RooGExpModel can be used
23 for 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 
39 using namespace std;
40 
41 ClassImp(RooGExpModel);
42 
43 ////////////////////////////////////////////////////////////////////////////////
44 
45 RooGExpModel::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 
59 RooGExpModel::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 
77 RooGExpModel::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 
95 RooGExpModel::RooGExpModel(const RooGExpModel& other, const char* name) :
96  RooResolutionModel(other,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 
111 RooGExpModel::~RooGExpModel()
112 {
113 }
114 
115 ////////////////////////////////////////////////////////////////////////////////
116 
117 Int_t RooGExpModel::basisCode(const char* name) const
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 
139 Double_t RooGExpModel::evaluate() const
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 
153  Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0. ;
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 
279 Double_t RooGExpModel::logErfC(Double_t xx) const
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 
297 std::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 
323 Double_t RooGExpModel::calcSinConv(Double_t sign, Double_t sig, Double_t tau, Double_t rtau, Double_t fsign) const
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 
345 Double_t RooGExpModel::calcDecayConv(Double_t sign, Double_t tau, Double_t sig, Double_t rtau, Double_t fsign) const
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 
428 Double_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 
463 Double_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 
496 Int_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 
538 Double_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 
553  Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
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 
692 std::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 
720 Double_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 
752 std::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 
766 Double_t RooGExpModel::evalCerfInt(Double_t sign, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
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 
787 std::complex<Double_t> RooGExpModel::evalCerfApprox(Double_t swt, Double_t u, Double_t c)
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 
800 Int_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 
808 void RooGExpModel::generateEvent(Int_t code)
809 {
810  R__ASSERT(code==1) ;
811  Double_t xgen ;
812  while(1) {
813  Double_t xgau = RooRandom::randomGenerator()->Gaus(0,(sigma*ssf));
814  Double_t xexp = RooRandom::uniform();
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 }