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

00001 /*****************************************************************************
00002  * Project: RooFit                                                           *
00003  * Package: RooFitCore                                                       *
00004  * @(#)root/roofitcore:$Id: RooTruthModel.cxx 24285 2008-06-16 15:05:15Z 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 // RooTruthModel is an implementation of RooResolution
00021 // model that provides a delta-function resolution model
00022 // <p>
00023 // The truth model supports <i>all</i> basis functions because it evaluates each basis function as  
00024 // as a RooFormulaVar.  The 6 basis functions used in B mixing and decay and 2 basis
00025 // functions used in D mixing have been hand coded for increased execution speed.
00026 // END_HTML
00027 //
00028 
00029 #include "RooFit.h"
00030 
00031 #include "Riostream.h"
00032 #include "Riostream.h"
00033 #include "RooTruthModel.h"
00034 
00035 ClassImp(RooTruthModel) 
00036 ;
00037 
00038 
00039 
00040 //_____________________________________________________________________________
00041 RooTruthModel::RooTruthModel(const char *name, const char *title, RooRealVar& xIn) : 
00042   RooResolutionModel(name,title,xIn)
00043 {  
00044   // Constructor of a truth resolution model, i.e. a delta function in observable 'xIn'
00045 
00046 }
00047 
00048 
00049 
00050 //_____________________________________________________________________________
00051 RooTruthModel::RooTruthModel(const RooTruthModel& other, const char* name) : 
00052   RooResolutionModel(other,name)
00053 {
00054   // Copy constructor
00055 }
00056 
00057 
00058 
00059 //_____________________________________________________________________________
00060 RooTruthModel::~RooTruthModel()
00061 {
00062   // Destructor
00063 }
00064 
00065 
00066 
00067 //_____________________________________________________________________________
00068 Int_t RooTruthModel::basisCode(const char* name) const 
00069 {
00070   // Return basis code for given basis definition string. Return special
00071   // codes for 'known' bases for which compiled definition exists. Return
00072   // generic bases code if implementation relies on TFormula interpretation
00073   // of basis name
00074 
00075   // Check for optimized basis functions
00076   if (!TString("exp(-@0/@1)").CompareTo(name)) return expBasisPlus ;
00077   if (!TString("exp(@0/@1)").CompareTo(name)) return expBasisMinus ;
00078   if (!TString("exp(-abs(@0)/@1)").CompareTo(name)) return expBasisSum ;
00079   if (!TString("exp(-@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisPlus ;
00080   if (!TString("exp(@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisMinus ;
00081   if (!TString("exp(-abs(@0)/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisSum ;
00082   if (!TString("exp(-@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisPlus ;
00083   if (!TString("exp(@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisMinus ;
00084   if (!TString("exp(-abs(@0)/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisSum ;
00085   if (!TString("(@0/@1)*exp(-@0/@1)").CompareTo(name)) return linBasisPlus ;
00086   if (!TString("(@0/@1)*(@0/@1)*exp(-@0/@1)").CompareTo(name)) return quadBasisPlus ;
00087   if (!TString("exp(-@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisPlus;
00088   if (!TString("exp(@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisMinus;
00089   if (!TString("exp(-abs(@0)/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisSum;
00090   if (!TString("exp(-@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisPlus;
00091   if (!TString("exp(@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisMinus;
00092   if (!TString("exp(-abs(@0)/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisSum;
00093 
00094   // Truth model is delta function, i.e. convolution integral
00095   // is basis function, therefore we can handle any basis function
00096   return genericBasis ;
00097 }
00098 
00099 
00100 
00101 
00102 //_____________________________________________________________________________
00103 void RooTruthModel::changeBasis(RooFormulaVar* inBasis) 
00104 {
00105   // Changes associated bases function to 'inBasis'
00106 
00107   // Process change basis function. Since we actually
00108   // evaluate the basis function object, we need to
00109   // adjust our client-server links to the basis function here
00110 
00111   // Remove client-server link to old basis
00112   if (_basis) {
00113     removeServer(*_basis) ;
00114   }
00115 
00116   // Change basis pointer and update client-server link
00117   _basis = inBasis ;
00118   if (_basis) {
00119     addServer(*_basis,kTRUE,kFALSE) ;
00120   }
00121 
00122   _basisCode = inBasis?basisCode(inBasis->GetTitle()):0 ;
00123 }
00124 
00125 
00126 
00127 
00128 
00129 //_____________________________________________________________________________
00130 Double_t RooTruthModel::evaluate() const 
00131 {
00132   // Evaluate the truth model: a delta function when used as PDF,
00133   // the basis function itself, when convoluted with a basis function.
00134 
00135   // No basis: delta function
00136   if (_basisCode == noBasis) {
00137     if (x==0) return 1 ;
00138     return 0 ;
00139   }
00140 
00141   // Generic basis: evaluate basis function object
00142   if (_basisCode == genericBasis) {
00143     return basis().getVal() ;
00144   }
00145 
00146   // Precompiled basis functions
00147   BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
00148   BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
00149 
00150   // Enforce sign compatibility
00151   if ((basisSign==Minus && x>0) || 
00152       (basisSign==Plus  && x<0)) return 0 ;
00153 
00154 
00155   Double_t tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ;
00156   // Return desired basis function
00157   switch(basisType) {    
00158   case expBasis: {
00159     return exp(-fabs((Double_t)x)/tau) ;
00160   }
00161   case sinBasis: {
00162     Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ; 
00163     return exp(-fabs((Double_t)x)/tau)*sin(x*dm) ;
00164   }
00165   case cosBasis: {
00166     Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ; 
00167     return exp(-fabs((Double_t)x)/tau)*cos(x*dm) ;
00168   }
00169   case linBasis: {
00170     Double_t tscaled = fabs((Double_t)x)/tau;
00171     return exp(-tscaled)*tscaled ;
00172   }
00173   case quadBasis: {
00174     Double_t tscaled = fabs((Double_t)x)/tau;
00175     return exp(-tscaled)*tscaled*tscaled;
00176   }  
00177   case sinhBasis: {
00178     Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ; 
00179     return exp(-fabs((Double_t)x)/tau)*sinh(x*dg/2) ;
00180   }
00181   case coshBasis: {
00182     Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ; 
00183     return exp(-fabs((Double_t)x)/tau)*cosh(x*dg/2) ;
00184   }
00185   default:
00186     assert(0) ;
00187   }
00188 
00189   return 0 ;
00190 }
00191 
00192 
00193 
00194 //_____________________________________________________________________________
00195 Int_t RooTruthModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const 
00196 {
00197   // Advertise analytical integrals for compiled basis functions and when used
00198   // as p.d.f without basis function.
00199 
00200   switch(_basisCode) {
00201 
00202   // Analytical integration capability of raw PDF
00203   case noBasis:
00204     if (matchArgs(allVars,analVars,convVar())) return 1 ;
00205     break ;
00206 
00207   // Analytical integration capability of convoluted PDF
00208   case expBasisPlus: 
00209   case expBasisMinus:
00210   case expBasisSum:
00211   case sinBasisPlus:
00212   case sinBasisMinus:
00213   case sinBasisSum:
00214   case cosBasisPlus:
00215   case cosBasisMinus:
00216   case cosBasisSum:
00217   case linBasisPlus:
00218   case quadBasisPlus:
00219   case sinhBasisPlus:
00220   case sinhBasisMinus:
00221   case sinhBasisSum:
00222   case coshBasisPlus:
00223   case coshBasisMinus:
00224   case coshBasisSum:
00225     if (matchArgs(allVars,analVars,convVar())) return 1 ;
00226     break ;
00227   }
00228 
00229   return 0 ;
00230 }
00231 
00232 
00233 
00234 //_____________________________________________________________________________
00235 Double_t RooTruthModel::analyticalIntegral(Int_t code, const char* rangeName) const 
00236 {
00237   // Implement analytical integrals when used as p.d.f and for compiled
00238   // basis functions.
00239 
00240 
00241   // Code must be 1
00242   assert(code==1) ;
00243 
00244   // Unconvoluted PDF
00245   if (_basisCode==noBasis) return 1 ;
00246 
00247   // Precompiled basis functions
00248   BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
00249   BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
00250   //cout << " calling RooTruthModel::analyticalIntegral with basisType " << basisType << endl; 
00251 
00252   Double_t tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ;
00253   switch (basisType) {
00254   case expBasis:
00255     {
00256       Double_t result(0) ;
00257       if (tau==0) return 1 ;
00258       if (basisSign != Minus) result += tau*(1-exp(-x.max(rangeName)/tau)) ;
00259       if (basisSign != Plus) result += tau*(1-exp(x.min(rangeName)/tau)) ;
00260       return result ;
00261     }
00262   case sinBasis:
00263     {
00264       Double_t result(0) ;
00265       if (tau==0) return 0 ;
00266       Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
00267       //if (basisSign != Minus) result += exp(-x.max(rangeName)/tau)*(-1/tau*sin(dm*x.max(rangeName)) - dm*cos(dm*x.max(rangeName))) + 1/tau;
00268       //if (basisSign != Plus)  result -= exp( x.min(rangeName)/tau)*(-1/tau*sin(dm*(-x.min(rangeName))) - dm*cos(dm*(-x.min(rangeName)))) + 1/tau ;
00269       if (basisSign != Minus) result += exp(-x.max(rangeName)/tau)*(-1/tau*sin(dm*x.max(rangeName)) - dm*cos(dm*x.max(rangeName))) + dm;  // fixed FMV 08/29/03
00270       if (basisSign != Plus)  result -= exp( x.min(rangeName)/tau)*(-1/tau*sin(dm*(-x.min(rangeName))) - dm*cos(dm*(-x.min(rangeName)))) + dm ;  // fixed FMV 08/29/03
00271       return result / (1/(tau*tau) + dm*dm) ;
00272     }
00273   case cosBasis:
00274     {
00275       Double_t result(0) ;
00276       if (tau==0) return 1 ;
00277       Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
00278       if (basisSign != Minus) result += exp(-x.max(rangeName)/tau)*(-1/tau*cos(dm*x.max(rangeName)) + dm*sin(dm*x.max(rangeName))) + 1/tau ;
00279       //if (basisSign != Plus)  result += exp( x.min(rangeName)/tau)*(-1/tau*cos(dm*(-x.min(rangeName))) - dm*sin(dm*(-x.min(rangeName)))) + 1/tau ;
00280       if (basisSign != Plus)  result += exp( x.min(rangeName)/tau)*(-1/tau*cos(dm*(-x.min(rangeName))) + dm*sin(dm*(-x.min(rangeName)))) + 1/tau ; // fixed FMV 08/29/03
00281       return result / (1/(tau*tau) + dm*dm) ;
00282     }
00283   case linBasis:
00284     {
00285       if (tau==0) return 0 ;
00286       Double_t t_max = x.max(rangeName)/tau ;
00287       return tau*( 1 - (1 + t_max)*exp(-t_max) ) ;
00288     }
00289   case quadBasis:
00290     {
00291       if (tau==0) return 0 ;
00292       Double_t t_max = x.max(rangeName)/tau ;
00293       return tau*( 2 - (2 + (2 + t_max)*t_max)*exp(-t_max) ) ;
00294     }
00295   case sinhBasis:
00296     {
00297       Double_t result(0) ;
00298       if (tau==0) return 0 ;
00299       Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
00300       Double_t taup = 2*tau/(2-tau*dg);
00301       Double_t taum = 2*tau/(2+tau*dg);
00302       if (basisSign != Minus) result += 0.5*( taup*(1-exp(-x.max(rangeName)/taup)) - taum*(1-exp(-x.max(rangeName)/taum)) ) ;
00303       if (basisSign != Plus)  result -= 0.5*( taup*(1-exp( x.min(rangeName)/taup)) - taum*(1-exp( x.min(rangeName)/taum)) ) ;
00304       return result ;
00305     }
00306   case coshBasis:
00307     {
00308       Double_t result(0) ;
00309       if (tau==0) return 1 ;
00310       Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
00311       Double_t taup = 2*tau/(2-tau*dg);
00312       Double_t taum = 2*tau/(2+tau*dg);
00313       if (basisSign != Minus) result += 0.5*( taup*(1-exp(-x.max(rangeName)/taup)) + taum*(1-exp(-x.max(rangeName)/taum)) ) ;
00314       if (basisSign != Plus)  result += 0.5*( taup*(1-exp( x.min(rangeName)/taup)) + taum*(1-exp( x.min(rangeName)/taum)) ) ;
00315       return result ;
00316     }
00317   default:
00318     assert(0) ;
00319   }
00320 
00321   assert(0) ;
00322   return 0 ;
00323 }
00324 
00325 
00326 
00327 //_____________________________________________________________________________
00328 Int_t RooTruthModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
00329 {
00330   // Advertise internal generator for observable x
00331   if (matchArgs(directVars,generateVars,x)) return 1 ;  
00332   return 0 ;
00333 }
00334 
00335 
00336 
00337 //_____________________________________________________________________________
00338 void RooTruthModel::generateEvent(Int_t code)
00339 {
00340   // Implement internal generator for observable x,
00341   // x=0 for all events following definition
00342   // of delta function
00343 
00344   assert(code==1) ;
00345   Double_t zero(0.) ;
00346   x = zero ;
00347   return;
00348 }