00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 
00023 
00024 
00025 
00026 
00027 
00028 
00029 
00030 #include "RooFit.h"
00031 
00032 #include "Riostream.h"
00033 #include "TArrayD.h"
00034 #include <math.h>
00035 
00036 #include "RooStepFunction.h"
00037 #include "RooAbsReal.h"
00038 #include "RooRealVar.h"
00039 #include "RooArgList.h"
00040 #include "RooMsgService.h"
00041 #include "RooMath.h"
00042 
00043 ClassImp(RooStepFunction)
00044   ;
00045 
00046 
00047 
00048 RooStepFunction::RooStepFunction()
00049 {
00050   
00051   _coefIter = _coefList.createIterator() ;
00052 }
00053                                  
00054 
00055 
00056 
00057 RooStepFunction::RooStepFunction(const char* name, const char* title, 
00058                                  RooAbsReal& x, const RooArgList& coefList, const RooArgList& boundaryList, Bool_t interpolate) :
00059   RooAbsReal(name, title),
00060   _x("x", "Dependent", this, x),
00061   _coefList("coefList","List of coefficients",this),
00062   _boundaryList("boundaryList","List of boundaries",this),
00063   _interpolate(interpolate)
00064 {
00065   
00066 
00067   _coefIter = _coefList.createIterator() ;
00068   TIterator* coefIter = coefList.createIterator() ;
00069   RooAbsArg* coef ;
00070   while((coef = (RooAbsArg*)coefIter->Next())) {
00071     if (!dynamic_cast<RooAbsReal*>(coef)) {
00072       cout << "RooStepFunction::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() 
00073            << " is not of type RooAbsReal" << endl ;
00074       assert(0) ;
00075     }
00076     _coefList.add(*coef) ;
00077   }
00078   delete coefIter ;
00079 
00080   TIterator* boundaryIter = boundaryList.createIterator() ;
00081   RooAbsArg* boundary ;
00082   while((boundary = (RooAbsArg*)boundaryIter->Next())) {
00083     if (!dynamic_cast<RooAbsReal*>(boundary)) {
00084       cout << "RooStepFunction::ctor(" << GetName() << ") ERROR: boundary " << boundary->GetName() 
00085            << " is not of type RooAbsReal" << endl ;
00086       assert(0) ;
00087     }
00088     _boundaryList.add(*boundary) ;
00089   }
00090 
00091   if (_boundaryList.getSize()!=_coefList.getSize()+1) {
00092     coutE(InputArguments) << "RooStepFunction::ctor(" << GetName() << ") ERROR: Number of boundaries must be number of coefficients plus 1" << endl ;
00093     throw string("RooStepFunction::ctor() ERROR: Number of boundaries must be number of coefficients plus 1") ;
00094   }
00095 
00096 }
00097 
00098 
00099 
00100 
00101 RooStepFunction::RooStepFunction(const RooStepFunction& other, const char* name) :
00102   RooAbsReal(other, name), 
00103   _x("x", this, other._x), 
00104   _coefList("coefList",this,other._coefList),
00105   _boundaryList("boundaryList",this,other._boundaryList),
00106   _interpolate(other._interpolate)
00107 {
00108   
00109   _coefIter = _coefList.createIterator();
00110   _boundIter = _boundaryList.createIterator();
00111 }
00112 
00113 
00114 
00115 
00116 RooStepFunction::~RooStepFunction()
00117 {
00118   
00119   delete _coefIter ;
00120   delete _boundIter ;
00121 }
00122 
00123 
00124 
00125 
00126 Double_t RooStepFunction::evaluate() const 
00127 {
00128   
00129   vector<double> b(_boundaryList.getSize()) ;
00130   vector<double> c(_coefList.getSize()+3) ;
00131   Int_t nb(0) ;
00132   _boundIter->Reset() ;
00133   RooAbsReal* boundary ;
00134   while ((boundary=(RooAbsReal*)_boundIter->Next())) {
00135     b[nb++] = boundary->getVal() ;
00136   }
00137 
00138   
00139   if ((_x<b[0]) || (_x>b[nb-1])) return 0 ;
00140 
00141   if (!_interpolate) {
00142 
00143     
00144     for (Int_t i=0;i<nb-1;i++){
00145       if (_x>b[i]&&_x<=b[i+1]) {
00146         return ((RooAbsReal*)_coefList.at(i))->getVal() ;
00147       }
00148     } 
00149     return 0 ;
00150 
00151   } else {
00152 
00153     
00154 
00155     
00156     c[0] = b[0] ; c[nb] = b[nb-1] ;
00157     for (Int_t i=0 ; i<nb-1 ; i++) {
00158       c[i+1] = (b[i]+b[i+1])/2 ;
00159     }
00160 
00161     
00162     Int_t nc(0) ;
00163     _coefIter->Reset() ;
00164     RooAbsReal* coef ;
00165     vector<double> y(_coefList.getSize()+3) ;
00166     y[nc++] = 0 ;
00167     while ((coef=(RooAbsReal*)_coefIter->Next())) {
00168       y[nc++] = coef->getVal() ;
00169     }
00170     y[nc++] = 0 ;
00171 
00172     for (Int_t i=0;i<nc-1;i++){
00173       if (_x>c[i]&&_x<=c[i+1]) {
00174         Double_t xx[2] ; xx[0]=c[i] ; xx[1]=c[i+1] ;
00175         Double_t yy[2] ; yy[0]=y[i] ; yy[1]=y[i+1] ;
00176         return RooMath::interpolate(xx,yy,2,_x) ;
00177       }
00178     } 
00179     return 0;   
00180   }  
00181 }
00182