// @(#)root/roostats:$Id: cranmer $
// Author: Kyle Cranmer, Akira Shibata
// Author: Giovanni Petrucciani (UCSD) (log-interpolation)
/*************************************************************************
* Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
//_________________________________________________
/*
BEGIN_HTML
END_HTML
*/
//
#include "RooFit.h"
#include "Riostream.h"
#include
#include "TMath.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooMsgService.h"
#include "RooTrace.h"
#include "TMath.h"
#include "RooStats/HistFactory/FlexibleInterpVar.h"
using namespace std;
ClassImp(RooStats::HistFactory::FlexibleInterpVar)
using namespace RooStats;
using namespace HistFactory;
//_____________________________________________________________________________
FlexibleInterpVar::FlexibleInterpVar()
{
// Default constructor
_paramIter = _paramList.createIterator() ;
_nominal = 0;
_interpBoundary=1.;
_logInit = kFALSE ;
TRACE_CREATE
}
//_____________________________________________________________________________
FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
const RooArgList& paramList,
Double_t nominal, vector low, vector high) :
RooAbsReal(name, title),
_paramList("paramList","List of paramficients",this),
_nominal(nominal), _low(low), _high(high), _interpBoundary(1.)
{
_logInit = kFALSE ;
_paramIter = _paramList.createIterator() ;
TIterator* paramIter = paramList.createIterator() ;
RooAbsArg* param ;
while((param = (RooAbsArg*)paramIter->Next())) {
if (!dynamic_cast(param)) {
coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
<< " is not of type RooAbsReal" << endl ;
assert(0) ;
}
_paramList.add(*param) ;
_interpCode.push_back(0); // default code
}
delete paramIter ;
TRACE_CREATE
}
//_____________________________________________________________________________
FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
const RooArgList& paramList,
double nominal, const RooArgList& low, const RooArgList& high) :
RooAbsReal(name, title),
_paramList("paramList","List of paramficients",this),
_nominal(nominal), _interpBoundary(1.)
{
RooFIter lowIter = low.fwdIterator() ;
RooAbsReal* val ;
while ((val = (RooAbsReal*) lowIter.next())) {
_low.push_back(val->getVal()) ;
}
RooFIter highIter = high.fwdIterator() ;
while ((val = (RooAbsReal*) highIter.next())) {
_high.push_back(val->getVal()) ;
}
_logInit = kFALSE ;
_paramIter = _paramList.createIterator() ;
TIterator* paramIter = paramList.createIterator() ;
RooAbsArg* param ;
while((param = (RooAbsArg*)paramIter->Next())) {
if (!dynamic_cast(param)) {
coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
<< " is not of type RooAbsReal" << endl ;
assert(0) ;
}
_paramList.add(*param) ;
_interpCode.push_back(0); // default code
}
delete paramIter ;
TRACE_CREATE
}
//_____________________________________________________________________________
FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
const RooArgList& paramList,
double nominal, vector low, vector high,
vector code) :
RooAbsReal(name, title),
_paramList("paramList","List of paramficients",this),
_nominal(nominal), _low(low), _high(high), _interpCode(code), _interpBoundary(1.)
{
_logInit = kFALSE ;
_paramIter = _paramList.createIterator() ;
TIterator* paramIter = paramList.createIterator() ;
RooAbsArg* param ;
while((param = (RooAbsArg*)paramIter->Next())) {
if (!dynamic_cast(param)) {
coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
<< " is not of type RooAbsReal" << endl ;
assert(0) ;
}
_paramList.add(*param) ;
}
delete paramIter ;
TRACE_CREATE
}
//_____________________________________________________________________________
FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title) :
RooAbsReal(name, title),
_paramList("paramList","List of coefficients",this),
_nominal(0), _interpBoundary(1.)
{
// Constructor of flat polynomial function
_logInit = kFALSE ;
_paramIter = _paramList.createIterator() ;
TRACE_CREATE
}
//_____________________________________________________________________________
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), _interpBoundary(other._interpBoundary)
{
// Copy constructor
_logInit = kFALSE ;
_paramIter = _paramList.createIterator() ;
TRACE_CREATE
}
//_____________________________________________________________________________
FlexibleInterpVar::~FlexibleInterpVar()
{
// Destructor
delete _paramIter ;
TRACE_DESTROY
}
//_____________________________________________________________________________
void FlexibleInterpVar::setInterpCode(RooAbsReal& param, int code){
int index = _paramList.index(¶m);
if(index<0){
coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName()
<< " is not in list" << endl ;
} else {
coutW(InputArguments) << "FlexibleInterpVar::setInterpCode : " << param.GetName()
<< " is now " << code << endl ;
_interpCode.at(index) = code;
}
// GHL: Adding suggestion by Swagato:
_logInit = kFALSE ;
setValueDirty();
}
//_____________________________________________________________________________
void FlexibleInterpVar::setAllInterpCodes(int code){
for(unsigned int i=0; i<_interpCode.size(); ++i){
_interpCode.at(i) = code;
}
// GHL: Adding suggestion by Swagato:
_logInit = kFALSE ;
setValueDirty();
}
//_____________________________________________________________________________
void FlexibleInterpVar::setNominal(Double_t newNominal){
coutW(InputArguments) << "FlexibleInterpVar::setNominal : nominal is now " << newNominal << endl ;
_nominal = newNominal;
_logInit = kFALSE ;
setValueDirty();
}
//_____________________________________________________________________________
void FlexibleInterpVar::setLow(RooAbsReal& param, Double_t newLow){
int index = _paramList.index(¶m);
if(index<0){
coutE(InputArguments) << "FlexibleInterpVar::setLow ERROR: " << param.GetName()
<< " is not in list" << endl ;
} else {
coutW(InputArguments) << "FlexibleInterpVar::setLow : " << param.GetName()
<< " is now " << newLow << endl ;
_low.at(index) = newLow;
}
_logInit = kFALSE ;
setValueDirty();
}
//_____________________________________________________________________________
void FlexibleInterpVar::setHigh(RooAbsReal& param, Double_t newHigh){
int index = _paramList.index(¶m);
if(index<0){
coutE(InputArguments) << "FlexibleInterpVar::setHigh ERROR: " << param.GetName()
<< " is not in list" << endl ;
} else {
coutW(InputArguments) << "FlexibleInterpVar::setHigh : " << param.GetName()
<< " is now " << newHigh << endl ;
_high.at(index) = newHigh;
}
_logInit = kFALSE ;
setValueDirty();
}
//_____________________________________________________________________________
void FlexibleInterpVar::printAllInterpCodes(){
for(unsigned int i=0; i<_interpCode.size(); ++i){
coutI(InputArguments) <<"interp code for " << _paramList.at(i)->GetName() << " = " << _interpCode.at(i) <GetName() << ": low value = " << _low.at(i) << endl;
if( _high.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": high value = " << _high.at(i) << endl;
}
}
//_____________________________________________________________________________
double FlexibleInterpVar::PolyInterpValue(int i, double x) const {
// code for polynomial interpolation used when interpCode=4
double boundary = _interpBoundary;
double x0 = boundary;
// cache the polynomial coefficient values
// which do not dpened on x but on the boundaries values
if (!_logInit) {
_logInit=kTRUE ;
unsigned int n = _low.size();
assert(n == _high.size() );
_polCoeff.resize(n*6) ;
for (unsigned int j = 0; j < n ; j++) {
// location of the 6 coefficient for the j-th variable
double * coeff = &_polCoeff[j * 6];
// GHL: Swagato's suggestions
double pow_up = std::pow(_high[j]/_nominal, x0);
double pow_down = std::pow(_low[j]/_nominal, x0);
double logHi = std::log(_high[j]) ;
double logLo = std::log(_low[j] );
double pow_up_log = _high[j] <= 0.0 ? 0.0 : pow_up * logHi;
double pow_down_log = _low[j] <= 0.0 ? 0.0 : -pow_down * logLo;
double pow_up_log2 = _high[j] <= 0.0 ? 0.0 : pow_up_log * logHi;
double pow_down_log2= _low[j] <= 0.0 ? 0.0 : -pow_down_log* logLo;
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
// cache coefficient of the polynomial
coeff[0] = 1./(8*x0) *( 15*A0 - 7*x0*S1 + x0*x0*A2);
coeff[1] = 1./(8*x0*x0) *(-24 + 24*S0 - 9*x0*A1 + x0*x0*S2);
coeff[2] = 1./(4*pow(x0, 3))*( - 5*A0 + 5*x0*S1 - x0*x0*A2);
coeff[3] = 1./(4*pow(x0, 4))*( 12 - 12*S0 + 7*x0*A1 - x0*x0*S2);
coeff[4] = 1./(8*pow(x0, 5))*( + 3*A0 - 3*x0*S1 + x0*x0*A2);
coeff[5] = 1./(8*pow(x0, 6))*( -8 + 8*S0 - 5*x0*A1 + x0*x0*S2);
}
}
// GHL: Swagato's suggestions
// if( _low[i] == 0 ) _low[i] = 0.0001;
// if( _high[i] == 0 ) _high[i] = 0.0001;
// get pointer to location of coefficients in the vector
const double * coeff = &_polCoeff.front() + 6*i;
double a = coeff[0];
double b = coeff[1];
double c = coeff[2];
double d = coeff[3];
double e = coeff[4];
double f = coeff[5];
// evaluate the 6-th degree polynomial using Horner's method
double value = 1. + x * (a + x * ( b + x * ( c + x * ( d + x * ( e + x * f ) ) ) ) );
return value;
}
//_____________________________________________________________________________
Double_t FlexibleInterpVar::evaluate() const
{
// Calculate and return value of polynomial
Double_t total(_nominal) ;
_paramIter->Reset() ;
RooAbsReal* param ;
//const RooArgSet* nset = _paramList.nset() ;
int i=0;
// TString name = GetName();
// if (name == TString("ZHWW_ll12_vzll_epsilon") )
// // std::cout << "evaluate flexible interp var - init flag is " << _logInit << std::endl;
while((param=(RooAbsReal*)_paramIter->Next())) {
// param->Print("v");
Int_t icode = _interpCode[i] ;
switch(icode) {
case 0: {
// piece-wise linear
if(param->getVal()>0)
total += param->getVal()*(_high[i] - _nominal );
else
total += param->getVal()*(_nominal - _low[i]);
break ;
}
case 1: {
// pice-wise log
if(param->getVal()>=0)
total *= pow(_high[i]/_nominal, +param->getVal());
else
total *= pow(_low[i]/_nominal, -param->getVal());
break ;
}
case 2: {
// parabolic with linear
double a = 0.5*(_high[i]+_low[i])-_nominal;
double b = 0.5*(_high[i]-_low[i]);
double c = 0;
if(param->getVal()>1 ){
total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
} else if(param->getVal()<-1 ) {
total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
} else {
total += a*pow(param->getVal(),2) + b*param->getVal()+c;
}
break ;
}
case 3: {
//parabolic version of log-normal
double a = 0.5*(_high[i]+_low[i])-_nominal;
double b = 0.5*(_high[i]-_low[i]);
double c = 0;
if(param->getVal()>1 ){
total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
} else if(param->getVal()<-1 ) {
total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
} else {
total += a*pow(param->getVal(),2) + b*param->getVal()+c;
}
break ;
}
case 4: {
double boundary = _interpBoundary;
double x = param->getVal();
//std::cout << icode << " param " << param->GetName() << " " << param->getVal() << " boundary " << boundary << std::endl;
if(x >= boundary)
{
total *= std::pow(_high[i]/_nominal, +param->getVal());
}
else if (x <= -boundary)
{
total *= std::pow(_low[i]/_nominal, -param->getVal());
}
else if (x != 0)
{
total *= PolyInterpValue(i, x);
}
break ;
}
default: {
coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: " << param->GetName()
<< " with unknown interpolation code" << endl ;
}
}
++i;
}
if(total<=0) {
total= TMath::Limits::Min();
}
return total;
}
void FlexibleInterpVar::printMultiline(ostream& os, Int_t contents,
Bool_t verbose, TString indent) const
{
RooAbsReal::printMultiline(os,contents,verbose,indent);
os << indent << "--- FlexibleInterpVar ---" << endl;
printFlexibleInterpVars(os);
}
void FlexibleInterpVar::printFlexibleInterpVars(ostream& os) const
{
_paramIter->Reset();
for (int i=0;i<(int)_low.size();i++) {
RooAbsReal* param=(RooAbsReal*)_paramIter->Next();
os << setw(36) << param->GetName()<<": "<