From ea2e100b194ff2279f31e3f993db19a003896ca0 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Tue, 1 Nov 2011 21:55:02 +0000
Subject: [PATCH] merge fix 41700 from Kyle in Histfactory
git-svn-id: http://root.cern.ch/svn/root/trunk@41701 27541ba8-7e3a-0410-8455-c3a389f83636
---
.../src/HistoToWorkspaceFactoryFast.cxx | 6 ++--
.../src/PiecewiseInterpolation.cxx | 30 ++++++++++++++++++-
2 files changed, 33 insertions(+), 3 deletions(-)
diff --git a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
index 6e028ae5486..03642d50e67 100644
--- a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
+++ b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
@@ -233,6 +233,7 @@ namespace HistFactory{
else if (classname.find("TH2")==0) { histndim=2; }
else if (classname.find("TH3")==0) { histndim=3; }
R__ASSERT( histndim==fObsNameVec.size() );
+ // cout <<"In LinInterpWithConstriants and histndim = " << histndim <<endl;
// create roorealvar observables
RooArgList observables;
@@ -269,11 +270,11 @@ namespace HistFactory{
temp = (RooRealVar*) proto->factory(("alpha_"+sourceName.at(j)+range).c_str());
// now add a constraint term for these parameters
- string command=("Gaussian::alpha_"+sourceName.at(j)+"Constraint(alpha_"+sourceName.at(j)+",nom_"+sourceName.at(j)+"[0.,-10,10],1.)");
+ string command=("Gaussian::alpha_"+sourceName.at(j)+"Constraint(alpha_"+sourceName.at(j)+",nom_alpha_"+sourceName.at(j)+"[0.,-10,10],1.)");
cout << command << endl;
constraintTermNames.push_back( proto->factory( command.c_str() )->GetName() );
proto->var(("nom_alpha_"+sourceName.at(j)).c_str())->setConstant();
- const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*proto->var(("nom_"+sourceName.at(j)).c_str()));
+ const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*proto->var(("nom_alpha_"+sourceName.at(j)).c_str()));
}
params.add(* temp );
}
@@ -299,6 +300,7 @@ namespace HistFactory{
PiecewiseInterpolation interp(prefix.c_str(),"",*nominalFunc,lowSet,highSet,params);
interp.setPositiveDefinite();
interp.setAllInterpCodes(0); // MB : change default to 1? = piece-wise log interpolation from pice-wise linear (=0)
+ // KC: interpo codes 1 etc. don't have proper analytic integral.
// cout << "check: " << interp.getVal() << endl;
proto->import(interp); // individual params have already been imported in first loop of this function
diff --git a/roofit/histfactory/src/PiecewiseInterpolation.cxx b/roofit/histfactory/src/PiecewiseInterpolation.cxx
index b0de946c127..56f5c72251d 100644
--- a/roofit/histfactory/src/PiecewiseInterpolation.cxx
+++ b/roofit/histfactory/src/PiecewiseInterpolation.cxx
@@ -240,6 +240,18 @@ Int_t PiecewiseInterpolation::getAnalyticalIntegralWN(RooArgSet& allVars, RooArg
if (_forceNumInt) return 0 ;
+ // KC: check if interCode=0 for all
+ RooFIter paramIterExtra(_paramSet.fwdIterator()) ;
+ int i=0;
+ while( paramIterExtra.next() ) {
+ if(_interpCode.at(i)!=0){
+ // can't factorize integral
+ cout <<"can't factorize integral"<<endl;
+ return 0;
+ }
+ ++i;
+ }
+
// Select subset of allVars that are actual dependents
analVars.add(allVars) ;
// RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0 ;
@@ -271,7 +283,8 @@ Int_t PiecewiseInterpolation::getAnalyticalIntegralWN(RooArgSet& allVars, RooArg
RooFIter highIter(_highSet.fwdIterator()) ;
RooFIter paramIter(_paramSet.fwdIterator()) ;
- int i=0;
+ // int i=0;
+ i=0;
while(paramIter.next() ) {
func = (RooAbsReal*)lowIter.next() ;
funcInt = func->createIntegral(analVars) ;
@@ -320,6 +333,19 @@ Double_t PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSe
i = 0;
RooFIter paramIter(_paramSet.fwdIterator()) ;
+ // KC: old interp code with new iterator
+ while( (param=(RooAbsReal*)paramIter.next()) ) {
+ low = (RooAbsReal*)lowIntIter.next() ;
+ high = (RooAbsReal*)highIntIter.next() ;
+
+ if(param->getVal()>0) {
+ value += param->getVal()*(high->getVal() - nominal );
+ } else {
+ value += param->getVal()*(nominal - low->getVal());
+ }
+ ++i;
+ }
+
/* // MB : old bit of interpolation code
while( (param=(RooAbsReal*)_paramIter->Next()) ) {
low = (RooAbsReal*)lowIntIter->Next() ;
@@ -334,6 +360,7 @@ Double_t PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSe
}
*/
+ /* KC: the code below is wrong. Can't pull out a constant change to a non-linear shape deformation.
while( (param=(RooAbsReal*)paramIter.next()) ) {
low = (RooAbsReal*)lowIntIter.next() ;
high = (RooAbsReal*)highIntIter.next() ;
@@ -381,6 +408,7 @@ Double_t PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSe
}
++i;
}
+ */
return value;
}
--
GitLab