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