diff --git a/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h b/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h index 9872f8ff4669704289f7ea7bf03e7404e194acfd..4c6a2f75de485ed7fe10d8c671ad16a211477d9e 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h @@ -38,6 +38,7 @@ namespace HistFactory{ void setInterpCode(RooAbsReal& param, int code); void setAllInterpCodes(int code); + void setGlobalBoundary(double boundary) {_interpBoundary = boundary;} void printAllInterpCodes(); @@ -48,16 +49,17 @@ namespace HistFactory{ protected: RooListProxy _paramList ; - double _nominal; + Double_t _nominal; vector<double> _low; vector<double> _high; vector<int> _interpCode; - + Double_t _interpBoundary; + TIterator* _paramIter ; //! do not persist Double_t evaluate() const; - ClassDef(RooStats::HistFactory::FlexibleInterpVar,1) // flexible interpolation + ClassDef(RooStats::HistFactory::FlexibleInterpVar,2) // flexible interpolation }; } } diff --git a/roofit/histfactory/src/FlexibleInterpVar.cxx b/roofit/histfactory/src/FlexibleInterpVar.cxx index 8812eb210931fd8bab6730780c9208e14a3b9d97..54a9e1852a806fbef68bce58c068d1e15ef73660 100644 --- a/roofit/histfactory/src/FlexibleInterpVar.cxx +++ b/roofit/histfactory/src/FlexibleInterpVar.cxx @@ -44,6 +44,7 @@ FlexibleInterpVar::FlexibleInterpVar() // Default constructor _paramIter = _paramList.createIterator() ; _nominal = 0; + _interpBoundary=1.; } @@ -53,7 +54,7 @@ FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title, double nominal, vector<double> low, vector<double> high) : RooAbsReal(name, title), _paramList("paramList","List of paramficients",this), - _nominal(nominal), _low(low), _high(high) + _nominal(nominal), _low(low), _high(high), _interpBoundary(1.) { _paramIter = _paramList.createIterator() ; @@ -81,7 +82,7 @@ FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title, vector<int> code) : RooAbsReal(name, title), _paramList("paramList","List of paramficients",this), - _nominal(nominal), _low(low), _high(high), _interpCode(code) + _nominal(nominal), _low(low), _high(high), _interpCode(code), _interpBoundary(1.) { _paramIter = _paramList.createIterator() ; @@ -115,7 +116,7 @@ FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title) : FlexibleInterpVar::FlexibleInterpVar(const FlexibleInterpVar& other, const char* name) : RooAbsReal(other, name), _paramList("paramList",this,other._paramList), - _nominal(other._nominal), _low(other._low), _high(other._high), _interpCode(other._interpCode) + _nominal(other._nominal), _low(other._low), _high(other._high), _interpCode(other._interpCode), _interpBoundary(other._interpBoundary) { // Copy constructor @@ -215,6 +216,64 @@ Double_t FlexibleInterpVar::evaluate() const total += a*pow(param->getVal(),2) + b*param->getVal()+c; } + } else if(_interpCode.at(i)==4){ // Aaron Armbruster - exponential extrapolation, polynomial interpolation + double boundary = _interpBoundary; + // piece-wise log + parabolic + if(param->getVal()>=boundary) + { + total *= pow(_high.at(i)/_nominal, +param->getVal()); + } + else if (param->getVal() < boundary && param->getVal() > -boundary && boundary != 0) + { + double x0 = boundary; + double x = param->getVal(); + +// double pow_up = pow(_high.at(i)/_nominal, x0); +// double pow_down = pow(_low.at(i)/_nominal, x0); +// double pow_up_log = pow_up*TMath::Log(_high.at(i)); +// double pow_down_log =-pow_down*TMath::Log(_low.at(i)); +// double pow_up_log2 = pow_up_log*TMath::Log(_high.at(i)); +// double pow_down_log2=-pow_down*TMath::Log(_low.at(i)); + + +//fcns+der are eq at bd +// double a = 1./(4*pow(x0, 1))*(3*A0 - x0*S1); +// double b = 1./(4*pow(x0, 2))*(4*S0 - x0*A1 - 8); +// double c = -1./(4*pow(x0, 3))*( A0 - x0*S1); +// double d = -1./(4*pow(x0, 4))*(2*S0 - x0*A1 - 4); +// total *= 1 + a*x + b*pow(x, 2) + c*pow(x, 3) + d*pow(x, 4); + + +//fcns+der+2nd_der are eq at bd + + double pow_up = pow(_high.at(i)/_nominal, x0); + double pow_down = pow(_low.at(i)/_nominal, x0); + double pow_up_log = pow_up*TMath::Log(_high.at(i)); + double pow_down_log = -pow_down*TMath::Log(_low.at(i)); + double pow_up_log2 = pow_up_log*TMath::Log(_high.at(i)); + double pow_down_log2= pow_down_log*TMath::Log(_low.at(i)); + + double S0 = (pow_up+pow_down)/2; + double A0 = (pow_up-pow_down)/2; + double S1 = (pow_up_log+pow_down_log)/2; + double A1 = (pow_up_log-pow_down_log)/2; + double S2 = (pow_up_log2+pow_down_log2)/2; + double A2 = (pow_up_log2-pow_down_log2)/2; + +//fcns+der+2nd_der are eq at bd + double a = 1./(8*pow(x0, 1))*( 15*A0 - 7*x0*S1 + x0*x0*A2); + double b = 1./(8*pow(x0, 2))*(-24 + 24*S0 - 9*x0*A1 + x0*x0*S2); + double c = 1./(4*pow(x0, 3))*( - 5*A0 + 5*x0*S1 - x0*x0*A2); + double d = 1./(4*pow(x0, 4))*( 12 - 12*S0 + 7*x0*A1 - x0*x0*S2); + double e = 1./(8*pow(x0, 5))*( + 3*A0 - 3*x0*S1 + x0*x0*A2); + double f = 1./(8*pow(x0, 6))*( -8 + 8*S0 - 5*x0*A1 + x0*x0*S2); + + total *= 1 + a*x + b*pow(x, 2) + c*pow(x, 3) + d*pow(x, 4) + e*pow(x, 5) + f*pow(x, 6); + } + else if (param->getVal()<=-boundary) + { + total *= pow(_low.at(i)/_nominal, -param->getVal()); + } } else { coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: " << param->GetName() << " with unknown interpolation code" << endl ; diff --git a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx index 13e99e50e900d9770bdf1fefa64f9d805738b0ef..23da6568e8a9b3b13104df9a75681213b4e73fc7 100644 --- a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx +++ b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx @@ -2098,6 +2098,7 @@ void HistoToWorkspaceFactoryFast::FormatFrameForLikelihood(RooPlot* frame, strin // Type 2 : RooPoisson RooPoisson pois(constrName.c_str(), constrName.c_str(), constrNom, constrMean); + pois.setNoRounding(true); proto->import( pois, RecycleConflictNodes() );