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