ROOT_528-00b_version: roofit/roofit/src/RooGaussModel.cxx Source File

00001 /*****************************************************************************
00002  * Project: RooFit                                                           *
00003  * Package: RooFitModels                                                     *
00004  * @(#)root/roofit:$Id: RooGaussModel.cxx 30333 2009-09-21 15:39:17Z wouter $
00005  * Authors:                                                                  *
00006  *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
00007  *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
00008  *                                                                           *
00009  * Copyright (c) 2000-2005, Regents of the University of California          *
00010  *                          and Stanford University. All rights reserved.    *
00011  *                                                                           *
00012  * Redistribution and use in source and binary forms,                        *
00013  * with or without modification, are permitted according to the terms        *
00014  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
00015  *****************************************************************************/
00016 
00017 //////////////////////////////////////////////////////////////////////////////
00018 //
00019 // BEGIN_HTML
00020 // Class RooGaussModel implements a RooResolutionModel that models a Gaussian
00021 // distribution. Object of class RooGaussModel can be used
00022 // for analytical convolutions with classes inheriting from RooAbsAnaConvPdf
00023 // END_HTML
00024 //
00025 
00026 #include "RooFit.h"
00027 
00028 #include "Riostream.h"
00029 #include "Riostream.h"
00030 #include "RooGaussModel.h"
00031 #include "RooMath.h"
00032 #include "RooRealConstant.h"
00033 #include "RooRandom.h"
00034 
00035 ClassImp(RooGaussModel) 
00036 ;
00037 
00038 
00039 
00040 //_____________________________________________________________________________
00041 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn, 
00042                              RooAbsReal& _mean, RooAbsReal& _sigma) :
00043   RooResolutionModel(name,title,xIn), 
00044   _flatSFInt(kFALSE),
00045   _asympInt(kFALSE),
00046   mean("mean","Mean",this,_mean),
00047   sigma("sigma","Width",this,_sigma),
00048   msf("msf","Mean Scale Factor",this,(RooRealVar&)RooRealConstant::value(1)),
00049   ssf("ssf","Sigma Scale Factor",this,(RooRealVar&)RooRealConstant::value(1))
00050 {  
00051 }
00052 
00053 
00054 
00055 //_____________________________________________________________________________
00056 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn, 
00057                              RooAbsReal& _mean, RooAbsReal& _sigma, 
00058                              RooAbsReal& _msSF) : 
00059   RooResolutionModel(name,title,xIn), 
00060   _flatSFInt(kFALSE),
00061   _asympInt(kFALSE),
00062   mean("mean","Mean",this,_mean),
00063   sigma("sigma","Width",this,_sigma),
00064   msf("msf","Mean Scale Factor",this,_msSF),
00065   ssf("ssf","Sigma Scale Factor",this,_msSF)
00066 {  
00067 }
00068 
00069 
00070 
00071 //_____________________________________________________________________________
00072 RooGaussModel::RooGaussModel(const char *name, const char *title, RooRealVar& xIn, 
00073                              RooAbsReal& _mean, RooAbsReal& _sigma, 
00074                              RooAbsReal& _meanSF, RooAbsReal& _sigmaSF) : 
00075   RooResolutionModel(name,title,xIn), 
00076   _flatSFInt(kFALSE),
00077   _asympInt(kFALSE),
00078   mean("mean","Mean",this,_mean),
00079   sigma("sigma","Width",this,_sigma),
00080   msf("msf","Mean Scale Factor",this,_meanSF),
00081   ssf("ssf","Sigma Scale Factor",this,_sigmaSF)
00082 {  
00083 }
00084 
00085 
00086 
00087 //_____________________________________________________________________________
00088 RooGaussModel::RooGaussModel(const RooGaussModel& other, const char* name) : 
00089   RooResolutionModel(other,name),
00090   _flatSFInt(other._flatSFInt),
00091   _asympInt(other._asympInt),
00092   mean("mean",this,other.mean),
00093   sigma("sigma",this,other.sigma),
00094   msf("msf",this,other.msf),
00095   ssf("ssf",this,other.ssf)
00096 {
00097 }
00098 
00099 
00100 
00101 //_____________________________________________________________________________
00102 RooGaussModel::~RooGaussModel()
00103 {
00104   // Destructor
00105 }
00106 
00107 
00108 
00109 //_____________________________________________________________________________
00110 Int_t RooGaussModel::basisCode(const char* name) const 
00111 {
00112   if (!TString("exp(-@0/@1)").CompareTo(name)) return expBasisPlus ;
00113   if (!TString("exp(@0/@1)").CompareTo(name)) return expBasisMinus ;
00114   if (!TString("exp(-abs(@0)/@1)").CompareTo(name)) return expBasisSum ;
00115   if (!TString("exp(-@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisPlus ;
00116   if (!TString("exp(@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisMinus ;
00117   if (!TString("exp(-abs(@0)/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisSum ;
00118   if (!TString("exp(-@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisPlus ;
00119   if (!TString("exp(@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisMinus ;
00120   if (!TString("exp(-abs(@0)/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisSum ;
00121   if (!TString("(@0/@1)*exp(-@0/@1)").CompareTo(name)) return linBasisPlus ;
00122   if (!TString("(@0/@1)*(@0/@1)*exp(-@0/@1)").CompareTo(name)) return quadBasisPlus ;
00123   if (!TString("exp(-@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisPlus;
00124   if (!TString("exp(@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisMinus;
00125   if (!TString("exp(-abs(@0)/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisSum;
00126   if (!TString("exp(-@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisPlus;
00127   if (!TString("exp(@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisMinus;
00128   if (!TString("exp(-abs(@0)/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisSum;
00129   return 0 ;
00130 } 
00131 
00132 
00133 
00134 //_____________________________________________________________________________
00135 Double_t RooGaussModel::evaluate() const 
00136 {  
00137   //cout << "RooGaussModel::evaluate(" << GetName() << ") basisCode = " << _basisCode << endl ;
00138   
00139   // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
00140   static Double_t root2(sqrt(2.)) ;
00141   static Double_t root2pi(sqrt(2*atan2(0.,-1.))) ;
00142   static Double_t rootpi(sqrt(atan2(0.,-1.))) ;
00143 
00144   BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
00145   BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
00146 
00147   Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
00148   if (basisType == coshBasis && _basisCode!=noBasis ) {
00149      Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
00150      if (dGamma==0) basisType = expBasis;
00151   }
00152 
00153   if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
00154     Double_t xprime = (x-(mean*msf))/(sigma*ssf) ;
00155     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 1st form" << endl ;
00156     
00157     Double_t result = exp(-0.5*xprime*xprime)/(sigma*ssf*root2pi) ;
00158     if (_basisCode!=0 && basisSign==Both) result *= 2 ;
00159     //cout << "1st form " << "x= " << x << " result= " << result << endl;
00160     return result ;
00161   }
00162 
00163   // *** 2nd form: 0, used for sinBasis, linBasis, and quadBasis with tau=0 ***
00164   if (tau==0) {
00165     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 2nd form" << endl ;
00166     return 0. ;
00167   }
00168 
00169   // *** 3nd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
00170   Double_t omega = (basisType==sinBasis || basisType==cosBasis)
00171     ? ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
00172   Double_t xprime = (x-(mean*msf))/tau ;
00173   Double_t c = (sigma*ssf)/(root2*tau) ; 
00174   Double_t u = xprime/(2*c) ;
00175 
00176   if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
00177     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 3d form tau=" << tau << endl ;
00178     Double_t result(0) ;
00179     if (basisSign!=Minus) result += exp(-xprime+c*c) * RooMath::erfc(-u+c) ;
00180     if (basisSign!=Plus)  result += exp(xprime+c*c) * RooMath::erfc(u+c) ;
00181     // equivalent form, added FMV, 07/24/03
00182     //if (basisSign!=Minus) result += evalCerfRe(-u,c) ; 
00183     //if (basisSign!=Plus)  result += evalCerfRe( u,c) ; 
00184     //cout << "3rd form " << "x= " << x << " result= " << result << endl;
00185     return result ;
00186   }
00187   
00188   // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
00189   Double_t wt = omega *tau ;
00190   if (basisType==sinBasis) {
00191     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() << ") 4th form omega = " 
00192                              << omega << ", tau = " << tau << endl ;
00193     Double_t result(0) ;
00194     if (wt==0.) return result ;
00195     if (basisSign!=Minus) result += -1*evalCerfIm(-wt,-u,c) ; 
00196     if (basisSign!=Plus) result += -1*evalCerfIm(wt,u,c) ; 
00197     //cout << "4th form " << "x= " << x << " result= " << result << endl;
00198     return result ;
00199   }
00200 
00201   // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
00202   if (basisType==cosBasis) {
00203     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() 
00204                              << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
00205     Double_t result(0) ;
00206     if (basisSign!=Minus) result += evalCerfRe(-wt,-u,c) ; 
00207     if (basisSign!=Plus) result += evalCerfRe(wt,u,c) ; 
00208     //cout << "5th form " << "x= " << x << " result= " << result << endl;
00209     return result ;  
00210   }
00211 
00212   // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
00213   if (basisType==linBasis) {
00214     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() 
00215                              << ") 6th form tau = " << tau << endl ;
00216 
00217     assert(basisSign==Plus);  // This should only be for positive times
00218 
00219     Double_t f0 = exp(-xprime+c*c) * RooMath::erfc(-u+c);
00220     Double_t f1 = exp(-u*u);
00221 
00222     return (xprime - 2*c*c)*f0 + (2*c/rootpi)*f1 ; 
00223   }
00224 
00225   // *** 7th form: Convolution with (t/tau)^2*exp(-t/tau), used for quadBasis ***
00226   if (basisType==quadBasis) {
00227     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() 
00228                              << ") 7th form tau = " << tau << endl ;
00229 
00230     assert(basisSign==Plus);  // This should only be for positive times
00231 
00232     Double_t f0 = exp(-xprime+c*c) * RooMath::erfc(-u+c);
00233     Double_t f1 = exp(-u*u);
00234 
00235     Double_t x2c2 = xprime - 2*c*c; 
00236 
00237     return ( x2c2*x2c2*f0 + (2*c/rootpi)*x2c2*f1 + 2*c*c*f0 );
00238   }
00239 
00240   // ***8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for         coshBasisSum ***
00241   if (basisType==coshBasis) {
00242     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() 
00243                              << ") 8th form tau = " << tau << endl ;
00244 
00245     Double_t dgamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
00246     Double_t tau1 = 1/(1/tau-dgamma/2) ; 
00247     Double_t tau2 = 1/(1/tau+dgamma/2) ;
00248     Double_t xprime1 = (x-(mean*msf))/tau1 ;
00249     Double_t c1 = (sigma*ssf)/(root2*tau1) ; 
00250     Double_t u1 = xprime1/(2*c1) ;
00251     Double_t xprime2 = (x-(mean*msf))/tau2 ;
00252     Double_t c2 = (sigma*ssf)/(root2*tau2) ; 
00253     Double_t u2 = xprime2/(2*c2) ;
00254     //Double_t c12 = c1*c1;
00255     //Double_t c22 = c2*c2;
00256 
00257     Double_t result(0);   
00258     //if (basisSign!=Minus) result += 0.5*(exp(-xprime1+c12) * RooMath::erfc(-u1+c1)+exp(-xprime2+c22) * RooMath::erfc(-u2+c2)) ;
00259     //if (basisSign!=Plus)  result += 0.5*(exp(xprime1+c12) * RooMath::erfc(u1+c1)+exp(xprime2+c22) * RooMath::erfc(u2+c2)) ;
00260     // equivalent form, added FMV, 07/24/03
00261     if (basisSign!=Minus) result += 0.5*(evalCerfRe(-u1,c1)+evalCerfRe(-u2,c2)) ; 
00262     if (basisSign!=Plus)  result += 0.5*(evalCerfRe( u1,c1)+evalCerfRe( u2,c2)) ; 
00263     //cout << "8th form " << "x= " << x << " result= " << result << endl;
00264     return result ;
00265   }
00266 
00267   // *** 9th form: Convolution with exp(-|t|/tau)*sinh(dgamma*t/2), used for        sinhBasisSum ***
00268   if (basisType==sinhBasis) {
00269     if (verboseEval()>2) cout << "RooGaussModel::evaluate(" << GetName() 
00270                              << ") 9th form tau = " << tau << endl ;
00271 
00272     Double_t dgamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
00273     Double_t tau1 = 1/(1/tau-dgamma/2) ; 
00274     Double_t tau2 = 1/(1/tau+dgamma/2) ;
00275     Double_t xprime1 = (x-(mean*msf))/tau1 ;
00276     Double_t c1 = (sigma*ssf)/(root2*tau1) ; 
00277     Double_t u1 = xprime1/(2*c1) ;
00278     Double_t xprime2 = (x-(mean*msf))/tau2 ;
00279     Double_t c2 = (sigma*ssf)/(root2*tau2) ; 
00280     Double_t u2 = xprime2/(2*c2) ;
00281     //Double_t c12 = c1*c1;
00282     //Double_t c22 = c2*c2;
00283 
00284     Double_t result(0);   
00285     //if (basisSign!=Minus) result += 0.5*(exp(-xprime1+c12) * RooMath::erfc(-u1+c1)-exp(-xprime2+c22) * RooMath::erfc(-u2+c2)) ;
00286     //if (basisSign!=Plus)  result += 0.5*(-exp(xprime1+c12) * RooMath::erfc(u1+c1)+exp(xprime2+c22) * RooMath::erfc(u2+c2)) ;
00287     // equivalent form, added FMV, 07/24/03
00288     if (basisSign!=Minus) result += 0.5*(evalCerfRe(-u1,c1)-evalCerfRe(-u2,c2)) ; 
00289     if (basisSign!=Plus)  result += 0.5*(evalCerfRe( u2,c2)-evalCerfRe( u1,c1)) ; 
00290     //cout << "9th form " << "x= " << x << " result= " << result << endl;
00291     return result ;
00292   }
00293 
00294 
00295   assert(0) ;
00296   return 0 ;
00297 }
00298 
00299 
00300 
00301 //_____________________________________________________________________________
00302 Int_t RooGaussModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const 
00303 {
00304   switch(_basisCode) {    
00305 
00306   // Analytical integration capability of raw PDF
00307   case noBasis:
00308     if (matchArgs(allVars,analVars,convVar())) return 1 ;
00309     break ;
00310 
00311   // Analytical integration capability of convoluted PDF
00312   case expBasisPlus:
00313   case expBasisMinus:
00314   case expBasisSum:
00315   case sinBasisPlus:
00316   case sinBasisMinus:
00317   case sinBasisSum:
00318   case cosBasisPlus:
00319   case cosBasisMinus:
00320   case cosBasisSum:
00321   case linBasisPlus:
00322   case quadBasisPlus:
00323   case coshBasisMinus:
00324   case coshBasisPlus:
00325   case coshBasisSum:
00326   case sinhBasisMinus:
00327   case sinhBasisPlus:
00328   case sinhBasisSum:
00329 
00330     // Optionally advertise flat integral over sigma scale factor
00331     if (_flatSFInt) {
00332       if (matchArgs(allVars,analVars,RooArgSet(convVar(),ssf.arg()))) {
00333         return 2 ;
00334       }
00335     }
00336     
00337     if (matchArgs(allVars,analVars,convVar())) return 1 ;
00338     break ;
00339   }
00340   
00341   return 0 ;
00342 }
00343 
00344 
00345 
00346 //_____________________________________________________________________________
00347 Double_t RooGaussModel::analyticalIntegral(Int_t code, const char* rangeName) const 
00348 {
00349   static Double_t root2 = sqrt(2.) ;
00350   //static Double_t rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
00351   static Double_t rootpi = sqrt(atan2(0.0,-1.0));
00352   Double_t ssfInt(1.0) ;
00353 
00354   // Code must be 1 or 2
00355   assert(code==1||code==2) ;
00356   if (code==2) {
00357     ssfInt = (ssf.max(rangeName)-ssf.min(rangeName)) ;
00358   }
00359 
00360   BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
00361   BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
00362 
00363   // *** 1st form: Straight Gaussian, used for unconvoluted PDF or expBasis with 0 lifetime ***
00364   Double_t tau = (_basisCode!=noBasis)?((RooAbsReal*)basis().getParameter(1))->getVal():0 ;
00365   if (basisType == coshBasis && _basisCode!=noBasis ) {
00366      Double_t dGamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
00367      if (dGamma==0) basisType = expBasis;
00368   }
00369   if (basisType==none || ((basisType==expBasis || basisType==cosBasis) && tau==0.)) {
00370     Double_t xscale = root2*(sigma*ssf);
00371     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 1st form" << endl ;
00372     
00373     Double_t xpmin = (x.min(rangeName)-(mean*msf))/xscale ;
00374     Double_t xpmax = (x.max(rangeName)-(mean*msf))/xscale ;
00375  
00376     Double_t result ;
00377     if (_asympInt) { // modified FMV, 07/24/03
00378       result = 1.0 ;
00379     } else {
00380       if (xpmin<-6 && xpmax>6) {
00381         // If integral is over >6 sigma, approximate with full integral
00382         result = 1.0 ;
00383       } else {
00384         result = 0.5*(RooMath::erf(xpmax)-RooMath::erf(xpmin)) ;
00385       }
00386     }
00387 
00388     if (_basisCode!=0 && basisSign==Both) result *= 2 ;    
00389     //cout << "Integral 1st form " << " result= " << result*ssfInt << endl;
00390     return result*ssfInt ;
00391   }
00392 
00393 
00394   Double_t omega = ((basisType==sinBasis)||(basisType==cosBasis)) ?
00395     ((RooAbsReal*)basis().getParameter(2))->getVal() : 0 ;
00396 
00397   // *** 2nd form: unity, used for sinBasis and linBasis with tau=0 (PDF is zero) ***
00398   if (tau==0) {
00399     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 2nd form" << endl ;
00400     return 0. ;
00401   }
00402 
00403   
00404 
00405   // *** 3rd form: Convolution with exp(-t/tau), used for expBasis and cosBasis(omega=0) ***
00406   Double_t c = (sigma*ssf)/(root2*tau) ; 
00407   Double_t xpmin = (x.min(rangeName)-(mean*msf))/tau ;
00408   Double_t xpmax = (x.max(rangeName)-(mean*msf))/tau ;
00409   Double_t umin = xpmin/(2*c) ;
00410   Double_t umax = xpmax/(2*c) ;
00411 
00412   if (basisType==expBasis || (basisType==cosBasis && omega==0.)) {
00413     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 3d form tau=" << tau << endl ;
00414 
00415     Double_t result(0) ;
00416     if (_asympInt) {   // modified FMV, 07/24/03
00417       if (basisSign!=Minus) result += 2 * tau ;
00418       if (basisSign!=Plus)  result += 2 * tau ;      
00419     } else {
00420       if (basisSign!=Minus) result += -1 * tau * ( RooMath::erf(-umax) - RooMath::erf(-umin) + 
00421                                                    exp(c*c) * ( exp(-xpmax)*RooMath::erfc(-umax+c)
00422                                                                 - exp(-xpmin)*RooMath::erfc(-umin+c) )) ;
00423       if (basisSign!=Plus)  result +=      tau * ( RooMath::erf(umax) - RooMath::erf(umin) + 
00424                                                    exp(c*c) * ( exp(xpmax)*RooMath::erfc(umax+c)
00425                                                                 - exp(xpmin)*RooMath::erfc(umin+c) )) ;     
00426       // equivalent form, added FMV, 07/24/03
00427       //if (basisSign!=Minus) result += evalCerfInt(+1,tau,-umin,-umax,c).re();   
00428       //if (basisSign!=Plus) result += evalCerfInt(-1,tau,umin,umax,c).re();
00429     }
00430     //cout << "Integral 3rd form " << " result= " << result*ssfInt << endl;
00431     return result*ssfInt ;
00432   }
00433 
00434   // *** 4th form: Convolution with exp(-t/tau)*sin(omega*t), used for sinBasis(omega<>0,tau<>0) ***
00435   Double_t wt = omega * tau ;
00436     
00437   if (basisType==sinBasis) {    
00438     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() << ") 4th form omega = " 
00439                              << omega << ", tau = " << tau << endl ;
00440     Double_t result(0) ;
00441     if (wt==0) return result*ssfInt ;
00442     if (basisSign!=Minus) {
00443       RooComplex evalDif(evalCerf(-wt,-umax,c) - evalCerf(-wt,-umin,c)) ;
00444       //result += -tau/(1+wt*wt) * ( -evalDif.im() +   -wt*evalDif.re() -   -wt*(RooMath::erf(-umax) - RooMath::erf(-umin)) ) ; 
00445       // FMV, fixed wrong sign, 07/24/03
00446       result += -tau/(1+wt*wt) * ( -evalDif.im() +   wt*evalDif.re() -   -wt*(RooMath::erf(-umax) - RooMath::erf(-umin)) ) ; 
00447     }
00448     if (basisSign!=Plus) {
00449       RooComplex evalDif(evalCerf(wt,umax,c) - evalCerf(wt,umin,c)) ;
00450       //result +=  tau/(1+wt*wt) * ( -evalDif.im() +    wt*evalDif.re() -    wt*(RooMath::erf(umax) - RooMath::erf(umin)) ) ;
00451       // FMV, fixed wrong sign, 07/24/03
00452       result +=  tau/(1+wt*wt) * ( -evalDif.im() +   -wt*evalDif.re() -    wt*(RooMath::erf(umax) - RooMath::erf(umin)) ) ;
00453     }
00454     // equivalent form, added FMV, 07/24/03
00455     //if (basisSign!=Minus) result += -1*evalCerfInt(+1,-wt,tau,-umin,-umax,c).im();
00456     //if (basisSign!=Plus) result += -1*evalCerfInt(-1,wt,tau,umin,umax,c).im();      
00457 
00458     //cout << "Integral 4th form " << " result= " << result*ssfInt << endl;
00459     return result*ssfInt ;
00460   }
00461 
00462   // *** 5th form: Convolution with exp(-t/tau)*cos(omega*t), used for cosBasis(omega<>0) ***
00463   if (basisType==cosBasis) {
00464     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName() 
00465                              << ") 5th form omega = " << omega << ", tau = " << tau << endl ;
00466     Double_t result(0) ;
00467     if (basisSign!=Minus) {
00468       RooComplex evalDif(evalCerf(-wt,-umax,c) - evalCerf(-wt,-umin,c)) ;
00469       //result += -tau/(1+wt*wt) * ( evalDif.re() + -wt*evalDif.im() + RooMath::erf(-umax) - RooMath::erf(-umin) ) ;
00470       // FMV, fixed wrong sign, 07/24/03
00471       result += -tau/(1+wt*wt) * ( evalDif.re() + wt*evalDif.im() + RooMath::erf(-umax) - RooMath::erf(-umin) ) ;
00472     }
00473     if (basisSign!=Plus) {
00474       RooComplex evalDif(evalCerf(wt,umax,c) - evalCerf(wt,umin,c)) ;
00475       //result +=  tau/(1+wt*wt) * ( evalDif.re() +  wt*evalDif.im() + RooMath::erf(umax) - RooMath::erf(umin) ) ;
00476       // FMV, fixed wrong sign, 07/24/03
00477       result +=  tau/(1+wt*wt) * ( evalDif.re() + -wt*evalDif.im() + RooMath::erf(umax) - RooMath::erf(umin) ) ;
00478     }
00479     // equivalent form, added FMV, 07/24/03
00480     //if (basisSign!=Minus) result += evalCerfInt(+1,-wt,tau,-umin,-umax,c).re();
00481     //if (basisSign!=Plus) result += evalCerfInt(-1,wt,tau,umin,umax,c).re();      
00482 
00483     //cout << "Integral 5th form " << " result= " << result*ssfInt << endl;
00484     return result*ssfInt ;
00485   }
00486 
00487   // *** 6th form: Convolution with (t/tau)*exp(-t/tau), used for linBasis ***
00488   if (basisType==linBasis) {
00489     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName()
00490                              << ") 6th form tau=" << tau << endl ;
00491 
00492     Double_t f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
00493     Double_t f1 = exp(-umax*umax) - exp(-umin*umin);
00494 
00495     Double_t tmp1 = exp(-xpmax)*RooMath::erfc(-umax + c);
00496     Double_t tmp2 = exp(-xpmin)*RooMath::erfc(-umin + c);
00497 
00498     Double_t f2 = tmp1 - tmp2;
00499     Double_t f3 = xpmax*tmp1 - xpmin*tmp2;
00500 
00501     Double_t expc2 = exp(c*c);
00502 
00503     return -tau*(              f0 +
00504                   (2*c/rootpi)*f1 +
00505              (1 - 2*c*c)*expc2*f2 +
00506                          expc2*f3
00507                 )*ssfInt;
00508   }
00509 
00510   // *** 7th form: Convolution with (t/tau)*(t/tau)*exp(-t/tau), used for quadBasis ***
00511   if (basisType==quadBasis) {
00512     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName()
00513                              << ") 7th form tau=" << tau << endl ;
00514 
00515     Double_t f0 = RooMath::erf(-umax) - RooMath::erf(-umin);
00516 
00517     Double_t tmpA1 = exp(-umax*umax);
00518     Double_t tmpA2 = exp(-umin*umin);
00519 
00520     Double_t f1 = tmpA1 - tmpA2;
00521     Double_t f2 = umax*tmpA1 - umin*tmpA2;
00522 
00523     Double_t tmpB1 = exp(-xpmax)*RooMath::erfc(-umax + c);
00524     Double_t tmpB2 = exp(-xpmin)*RooMath::erfc(-umin + c);
00525 
00526     Double_t f3 = tmpB1 - tmpB2;
00527     Double_t f4 = xpmax*tmpB1 - xpmin*tmpB2;
00528     Double_t f5 = xpmax*xpmax*tmpB1 - xpmin*xpmin*tmpB2;
00529 
00530     Double_t expc2 = exp(c*c);
00531 
00532     return -tau*( 2*f0 +
00533                   (4*c/rootpi)*((1-c*c)*f1 + c*f2) +
00534                   (2*c*c*(2*c*c-1) + 2)*expc2*f3 - (4*c*c-2)*expc2*f4 + expc2*f5
00535                 )*ssfInt;
00536   }
00537 
00538   // *** 8th form: Convolution with exp(-|t|/tau)*cosh(dgamma*t/2), used for coshBasis ***
00539   if (basisType==coshBasis) {
00540     if (verboseEval()>0) {cout << "RooGaussModel::analyticalIntegral(" << GetName()                             << ") 8th form tau=" << tau << endl ; }
00541     
00542     Double_t dgamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
00543     Double_t tau1 = 1/(1/tau-dgamma/2) ; 
00544     Double_t tau2 = 1/(1/tau+dgamma/2) ;
00545     Double_t c1 = (sigma*ssf)/(root2*tau1) ; 
00546     Double_t xpmin1 = (x.min(rangeName)-(mean*msf))/tau1 ;
00547     Double_t xpmax1 = (x.max(rangeName)-(mean*msf))/tau1 ;
00548     Double_t umin1 = xpmin1/(2*c1) ;
00549     Double_t umax1 = xpmax1/(2*c1) ;
00550     Double_t c2 = (sigma*ssf)/(root2*tau2) ; 
00551     Double_t xpmin2 = (x.min(rangeName)-(mean*msf))/tau2 ;
00552     Double_t xpmax2 = (x.max(rangeName)-(mean*msf))/tau2 ;
00553     Double_t umin2 = xpmin2/(2*c2) ;
00554     Double_t umax2 = xpmax2/(2*c2) ;
00555     //Double_t c12 = c1*c1;
00556     //Double_t c22 = c2*c2;
00557     //Double_t ec12 = exp(c12);
00558     //Double_t ec22 = exp(c22);
00559     
00560     Double_t result(0) ;
00561     
00562     /*
00563     if (basisSign!=Minus) result += -0.5*(tau1 * ( RooMath::erf(-umax1) - RooMath::erf(-umin1) + 
00564                                                    ec12 * ( exp(-xpmax1)*RooMath::erfc(-umax1+c1)
00565                                                                 - exp(-xpmin1)*RooMath::erfc(-umin1+c1) )) +   
00566                                                    tau2 * ( RooMath::erf(-umax2) - RooMath::erf(-umin2) + 
00567                                                    ec22 * ( exp(-xpmax2)*RooMath::erfc(-umax2+c2)
00568                                                                 - exp(-xpmin2)*RooMath::erfc(-umin2+c2) ))) ;
00569       if (basisSign!=Plus)  result +=  0.5*(tau1 * ( RooMath::erf(umax1) - RooMath::erf(umin1) + 
00570                                                    ec12 * ( exp(xpmax1)*RooMath::erfc(umax1+c1)
00571                                                                 - exp(xpmin1)*RooMath::erfc(umin1+c1) ))+ 
00572                                                    tau2 * ( RooMath::erf(umax2) - RooMath::erf(umin2) + 
00573                                                    ec22 * ( exp(xpmax2)*RooMath::erfc(umax2+c2)
00574                                                                 - exp(xpmin2)*RooMath::erfc(umin2+c2) ))) ;  
00575     */
00576     // equivalent form, added FMV, 07/24/03
00577     if (basisSign!=Minus) result += 0.5*(evalCerfInt(+1,tau1,-umin1,-umax1,c1)+
00578                                          evalCerfInt(+1,tau2,-umin2,-umax2,c2));
00579     if (basisSign!=Plus)  result += 0.5*(evalCerfInt(-1,tau1,umin1,umax1,c1)+
00580                                          evalCerfInt(-1,tau2,umin2,umax2,c2));
00581     
00582     //cout << "Integral 8th form " << " result= " << result*ssfInt << endl;
00583     return result*ssfInt ;
00584   }
00585    
00586   // *** 9th form: Convolution with exp(-|t|/tau)*sinh(dgamma*t/2), used for sinhBasis ***
00587   if (basisType==sinhBasis) {
00588     if (verboseEval()>0) cout << "RooGaussModel::analyticalIntegral(" << GetName()                             << ") 9th form tau=" << tau << endl ; 
00589     
00590     Double_t dgamma = ((RooAbsReal*)basis().getParameter(2))->getVal();
00591     Double_t tau1 = 1/(1/tau-dgamma/2) ; 
00592     Double_t tau2 = 1/(1/tau+dgamma/2) ;
00593     Double_t c1 = (sigma*ssf)/(root2*tau1) ; 
00594     Double_t xpmin1 = (x.min(rangeName)-(mean*msf))/tau1 ;
00595     Double_t xpmax1 = (x.max(rangeName)-(mean*msf))/tau1 ;
00596     Double_t umin1 = xpmin1/(2*c1) ;
00597     Double_t umax1 = xpmax1/(2*c1) ;
00598     Double_t c2 = (sigma*ssf)/(root2*tau2) ; 
00599     Double_t xpmin2 = (x.min(rangeName)-(mean*msf))/tau2 ;
00600     Double_t xpmax2 = (x.max(rangeName)-(mean*msf))/tau2 ;
00601     Double_t umin2 = xpmin2/(2*c2) ;
00602     Double_t umax2 = xpmax2/(2*c2) ;
00603     //Double_t c12 = c1*c1;
00604     //Double_t c22 = c2*c2;
00605     //Double_t ec12 = exp(c12);
00606     //Double_t ec22 = exp(c22);
00607  
00608     Double_t result(0) ;
00609 
00610     /*
00611     if (basisSign!=Minus) result += 0.5*(-tau1 * ( RooMath::erf(-umax1) - RooMath::erf(-umin1) + 
00612                                                    ec12 * ( exp(-xpmax1)*RooMath::erfc(-umax1+c1)
00613                                                             - exp(-xpmin1)*RooMath::erfc(-umin1+c1) )) +   
00614                                                    tau2 * ( RooMath::erf(-umax2) - RooMath::erf(-umin2) + 
00615                                                             ec22 * ( exp(-xpmax2)*RooMath::erfc(-umax2+c2)
00616                                                                      - exp(-xpmin2)*RooMath::erfc(-umin2+c2) ))) ;
00617     if (basisSign!=Plus)  result += 0.5*(-tau1 * ( RooMath::erf(umax1) - RooMath::erf(umin1) + 
00618                                                    ec12 * ( exp(xpmax1)*RooMath::erfc(umax1+c1)
00619                                                             - exp(xpmin1)*RooMath::erfc(umin1+c1) ))+ 
00620                                          tau2 * ( RooMath::erf(umax2) - RooMath::erf(umin2) + 
00621                                                   ec22 * ( exp(xpmax2)*RooMath::erfc(umax2+c2)
00622                                                            - exp(xpmin2)*RooMath::erfc(umin2+c2) ))) ; 
00623     */
00624     // equivalent form, added FMV, 07/24/03
00625     if (basisSign!=Minus) result += 0.5*(evalCerfInt(+1,tau1,-umin1,-umax1,c1)-
00626                                          evalCerfInt(+1,tau2,-umin2,-umax2,c2));
00627     if (basisSign!=Plus)  result += 0.5*(evalCerfInt(-1,tau2,umin2,umax2,c2)-
00628                                          evalCerfInt(-1,tau1,umin1,umax1,c1));
00629 
00630     //cout << "Integral 9th form " << " result= " << result*ssfInt << endl;
00631     return result*ssfInt ;
00632     
00633   }
00634   assert(0) ;
00635   return 0 ;
00636 }
00637 
00638 
00639 
00640 //_____________________________________________________________________________
00641 RooComplex RooGaussModel::evalCerfApprox(Double_t swt, Double_t u, Double_t c) const
00642 {
00643   // use the approximation: erf(z) = exp(-z*z)/(sqrt(pi)*z)
00644   // to explicitly cancel the divergent exp(y*y) behaviour of
00645   // CWERF for z = x + i y with large negative y
00646 
00647   static Double_t rootpi= sqrt(atan2(0.,-1.));
00648   RooComplex z(swt*c,u+c);  
00649   RooComplex zc(u+c,-swt*c);
00650   RooComplex zsq= z*z;
00651   RooComplex v= -zsq - u*u;
00652 
00653   return v.exp()*(-zsq.exp()/(zc*rootpi) + 1)*2 ;
00654 }
00655 
00656 
00657 
00658 // added FMV, 07/24/03
00659 //_____________________________________________________________________________
00660 RooComplex RooGaussModel::evalCerfInt(Double_t sign, Double_t wt, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
00661 {
00662   RooComplex diff;
00663   if (_asympInt) {
00664     diff = RooComplex(2,0) ;
00665   } else {
00666     diff = RooComplex(sign,0.)*(evalCerf(wt,umin,c) - evalCerf(wt,umax,c) + RooMath::erf(umin) - RooMath::erf(umax));
00667   }
00668   return RooComplex(tau/(1.+wt*wt),0)*RooComplex(1,wt)*diff;
00669 }
00670 // added FMV, 08/17/03
00671 
00672 //_____________________________________________________________________________
00673 Double_t RooGaussModel::evalCerfInt(Double_t sign, Double_t tau, Double_t umin, Double_t umax, Double_t c) const
00674 {
00675   Double_t diff;
00676   if (_asympInt) {
00677     diff = 2. ;
00678   } else {
00679     if ((umin<-8 && umax>8)||(umax<-8 && umin>8)) {
00680       // If integral is over >8 sigma, approximate with full integral
00681       diff = 2. ;
00682     } else {
00683       diff = sign*(evalCerfRe(umin,c) - evalCerfRe(umax,c) + RooMath::erf(umin) - RooMath::erf(umax));
00684     }
00685   }
00686   return tau*diff;
00687 }
00688 
00689 
00690 
00691 //_____________________________________________________________________________
00692 Int_t RooGaussModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
00693 {
00694   if (matchArgs(directVars,generateVars,x)) return 1 ;  
00695   return 0 ;
00696 }
00697 
00698 
00699 
00700 //_____________________________________________________________________________
00701 void RooGaussModel::generateEvent(Int_t code)
00702 {
00703   assert(code==1) ;
00704   Double_t xgen ;
00705   while(1) {    
00706     xgen = RooRandom::randomGenerator()->Gaus((mean*msf),(sigma*ssf));
00707     if (xgen<x.max() && xgen>x.min()) {
00708       x = xgen ;
00709       return ;
00710     }
00711   }
00712 }
00713 
00714 
00715 
00716