From 69faa938d3fb5bc3d3acc1c4ea0f8841d7dc89bd Mon Sep 17 00:00:00 2001 From: Lucas Corcodilos Date: Thu, 16 May 2019 10:23:01 -0500 Subject: [PATCH 1/5] Changes to allow for fitting of 2D parametric hist --- interface/RooParametricHist2D.h | 54 +++++++ interface/chebyshevBasis.h | 54 +++++++ python/ShapeTools.py | 22 ++- src/RooParametricHist2D.cxx | 250 ++++++++++++++++++++++++++++++++ src/chebyshevBasis.cxx | 180 +++++++++++++++++++++++ src/classes.h | 2 + src/classes_def.xml | 2 + 7 files changed, 561 insertions(+), 3 deletions(-) create mode 100644 interface/RooParametricHist2D.h create mode 100644 interface/chebyshevBasis.h create mode 100644 src/RooParametricHist2D.cxx create mode 100644 src/chebyshevBasis.cxx diff --git a/interface/RooParametricHist2D.h b/interface/RooParametricHist2D.h new file mode 100644 index 00000000000..e658b40ed22 --- /dev/null +++ b/interface/RooParametricHist2D.h @@ -0,0 +1,54 @@ +#ifndef ROOPARAMETRICHIST2D +#define ROOPARAMETRICHIST2D + +#include "RooRealProxy.h" +#include "RooListProxy.h" +#include "RooAbsPdf.h" +#include "RooAddition.h" +#include "RooAbsReal.h" +//#include "RooAbsData.h" +#include "TH1.h" +#include "TH2F.h" +#include "TH2D.h" + +class RooParametricHist2D : public RooAbsPdf { +public: + + RooParametricHist2D() {} ; + RooParametricHist2D (const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _y, RooArgList& _pars, const TH2 &_shape); + + RooParametricHist2D(const RooParametricHist2D& other, const char* name=0); + virtual TObject* clone(const char* newname) const { return new RooParametricHist2D(*this,newname); } + inline virtual ~RooParametricHist2D (){}; + Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char* rangeName=0) const ; + Double_t analyticalIntegral(Int_t, const char* rangeName=0) const ; + + //RooAddition & getYieldVar(){return sum;}; + +protected: + + RooRealProxy x; + RooRealProxy y; + //RooAddition sum; + RooListProxy pars; + mutable int N_bins; + mutable int N_bins_x; + mutable int N_bins_y; + mutable std::vector bins_x; + mutable std::vector bins_y; + mutable std::vector widths_x; + mutable std::vector widths_y; + + void initializeBins(const TH2&) const; + //void initializeNorm(); + + Double_t evaluate() const ; + double getFullSum() const ; + + mutable double cval; + void update_cval(double r){cval=r;}; +private: + ClassDef(RooParametricHist2D, 1) +}; + +#endif diff --git a/interface/chebyshevBasis.h b/interface/chebyshevBasis.h new file mode 100644 index 00000000000..c085ef437dd --- /dev/null +++ b/interface/chebyshevBasis.h @@ -0,0 +1,54 @@ +#ifndef CHEBYSHEVBASIS +#define CHEBYSHEVBASIS + +#include "TH2F.h" +#include "TCanvas.h" +#include "TFile.h" +#include "RooFit.h" +#include "RooFormulaVar.h" +#include "RooConstVar.h" +#include "RooArgList.h" +#include "RooAddition.h" +#include "RooAbsArg.h" +#include "Riostream.h" +#include "RooRealVar.h" +#include + +class chebyshevBasis : public TObject { +public: + chebyshevBasis() {}; + chebyshevBasis( const char *name, + const char *title, + const Int_t _orderX, + const Int_t _orderY, + const TH2 &_binning ); + inline virtual ~chebyshevBasis (){}; + // float getChebyshevVal(float xval, float yval, float thisOrderX, float thisOrderY) const; + RooAddition getBinVal(const float xCenter, const float yCenter) const; + void drawBasis(std::string file_name); + +protected: + // variables + TH2F* binning; + Int_t orderX; + Int_t orderY; + RooArgList coefList; + Double_t xmin; + Double_t xmax; + Double_t ymin; + Double_t ymax; + float slope_x; + float slope_y; + // std::map polyHists; + + // functions + std::pair mapToChebyshev(float ix,float iy) const; + float Eval2DChebyshev(float x, float y, int thisOrderX, int thisOrderY) const; + // TH2F Make2DChebyshev(const int& thisOrderX, const int& thisOrderY) const; + + +private: + ClassDef(chebyshevBasis, 1) +}; + +#endif \ No newline at end of file diff --git a/python/ShapeTools.py b/python/ShapeTools.py index c350af6857b..c1460121554 100644 --- a/python/ShapeTools.py +++ b/python/ShapeTools.py @@ -376,8 +376,15 @@ def prepareAllShapes(self): if self.options.verbose > 1: stderr.write("Observables: %s\n" % str(shapeObs.keys())) if len(shapeObs.keys()) != 1: self.out.binVars = ROOT.RooArgSet() - for obs in shapeObs.values(): - self.out.binVars.add(obs, True) + # Custom code starts - LC 5/14/19 + for obs_key in shapeObs.keys(): + # RooArgSet does not have method to add another RooArgSet (RooCollection does!) so we have to manually loop over the RooArgSet contents (in obs) + if ',' in obs_key: # if multiple vars and looking at a RooArgSet value + for k in obs_key.split(','): + self.out.binVars.add(shapeObs[obs_key].find(k), True) + else: + self.out.binVars.add(shapeObs[obs_key], True) + # Custom code ends else: self.out.binVars = shapeObs.values()[0] self.out._import(self.out.binVars) @@ -638,7 +645,16 @@ def getPdf(self,channel,process,_cache={}): pdfs.Add(self.shape2Pdf(shapeDown,channel,process)) histpdf = nominalPdf if nominalPdf.InheritsFrom("RooDataHist") else nominalPdf.dataHist() xvar = histpdf.get().first() - rhp = ROOT.FastVerticalInterpHistPdf2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, pdfs, coeffs, qrange, qalgo) + # Custom code start - LC 5/14/19 + histpdfSubset = histpdf.get().Clone() + histpdfSubset.remove(xvar) + if histpdfSubset.first() != False: + yvar = histpdfSubset.first() # y is a little more complex + rhp = ROOT.FastVerticalInterpHistPdf2D2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, yvar, False, pdfs, coeffs, qrange, qalgo) + else: + rhp = ROOT.FastVerticalInterpHistPdf2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, pdfs, coeffs, qrange, qalgo) + # Custom code ends (plus the next line was commented out) + #rhp = ROOT.FastVerticalInterpHistPdf2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, pdfs, coeffs, qrange, qalgo) _cache[(channel,process)] = rhp return rhp else: diff --git a/src/RooParametricHist2D.cxx b/src/RooParametricHist2D.cxx new file mode 100644 index 00000000000..9fb10aaa0d8 --- /dev/null +++ b/src/RooParametricHist2D.cxx @@ -0,0 +1,250 @@ +/***************************************************************************** + *****************************************************************************/ + + +#include "Riostream.h" + +#include "RooAbsData.h" +#include "RooAbsPdf.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooParametricHist2D.h" + +#include +#include "TMath.h" +#include "RooFormulaVar.h" +#include "RooAbsReal.h" +#include "RooFit.h" + +#include "TFile.h" +#include + +//using namespace RooFit ; + +ClassImp(RooParametricHist2D) + +RooParametricHist2D::RooParametricHist2D( const char *name, + const char *title, + RooAbsReal& _x, + RooAbsReal& _y, + RooArgList& _pars, // stores all bins - Ex. a 3x3 space [[1,2,3],[4,5,6],[7,8,9]] would be stored [1,2,3,4,5,6,7,8,9] + const TH2 &_shape ) :// need this to sort bins into x and y + RooAbsPdf(name,title), + x("observable_x","observable_x",this,_x), + y("observable_y","observable_y",this,_y), + pars("pars","pars",this) + //SM_shape("SM_shape","SM_shape",this,_SM_shape), +{ + // std::cout << "Constructing..." << std::endl; + // std::cout << "_pars size... " << _pars.getSize() << std::endl; + // std::cout << "_pars type... " << typeid(_pars.first()).name() << std::endl; + + TIterator *varIter=_pars.createIterator(); + RooAbsReal *fVar; + while ( (fVar = (RooAbsReal*)varIter->Next()) ){ + // std::cout << "fVar type... " << typeid(fVar).name() << std::endl; + // std::cout << "*fVar type... " << typeid(*fVar).name() << std::endl; + + pars.add(*fVar); + } + + // std::cout << "Constructed..." << std::endl; + // std::cout << "pars size... " << pars.getSize() << std::endl; + // std::cout << "pars type... " << typeid(pars.first()).name() << std::endl; + + if ( pars.getSize() != _shape.GetNbinsY()*_shape.GetNbinsX() ){ + std::cout << " Warning, number of parameters not equal to number of bins in shape histogram! " << std::endl; + } + + initializeBins(_shape); + // std::cout << "Bins initialized" << std::endl; + +// initializeNorm(); + cval = -1; +} + +//_____________________________________________________________________________ +RooParametricHist2D::RooParametricHist2D( const RooParametricHist2D& other, + const char* name) : + RooAbsPdf(other, name), + x("observable_x",this,other.x), + y("observable_y",this,other.y), + pars("_pars",this,RooListProxy()) +{ + //std::cout << "Cloning RooParametricHist2D with name " <Next()) ){ + // std::cout << fVar->GetName() << ": " << typeid(&fVar).name() << std::endl; + pars.add(*fVar); + } + // std::cout << "Check 2" << std::endl; + + for(int i=0; i<=N_bins_x; ++i) { + bins_x.push_back(other.bins_x[i]); + if (i bins; + //std::vector widths; + N_bins_x = shape.GetNbinsX(); + N_bins_y = shape.GetNbinsY(); + + + for(int i=1; i<=N_bins_x+1; ++i) { + bins_x.push_back(shape.GetXaxis()->GetBinLowEdge(i)); + if (i<=N_bins_x) widths_x.push_back(shape.GetXaxis()->GetBinWidth(i)); + } + + for(int j=1; j<=N_bins_y+1; ++j) { + bins_y.push_back(shape.GetYaxis()->GetBinLowEdge(j)); + if (j<=N_bins_y) widths_y.push_back(shape.GetYaxis()->GetBinWidth(j)); + } + // std::cout << "Initialize bins_x = "; + // for (auto i = bins_x.begin(); i != bins_x.end(); ++i){ + // std::cout << *i << ' '; + // } + // std::cout << std::endl; + + + // std::cout << "Initialize bins_y = "; + // for (auto i = bins_y.begin(); i != bins_y.end(); ++i){ + // std::cout << *i << ' '; + // } + // std::cout << std::endl; + +} + +double RooParametricHist2D::getFullSum() const { + //std::cout << "Getting full sum" << std::endl; + double sum=0; + + TIterator *varIter=pars.createIterator(); + RooAbsReal *fVar; + while ( (fVar = (RooAbsReal*)varIter->Next()) ){ + sum+=fVar->getVal(); + } + + return sum; +} + +Int_t RooParametricHist2D::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet & analVars, const char*) const { + if (matchArgs(allVars,analVars,x,y)){ + return 1; + } + return 0; +} + +Double_t RooParametricHist2D::analyticalIntegral(Int_t code, const char* rangeName) const +{ + assert(code==1) ; + + // Case without range is trivial: p.d.f is by construction normalized + if (!rangeName) { + //return 1;//getFullSum() ; + return getFullSum(); + } else { + std::cout << "Analytical integral for range " << rangeName << " in RooParametricHist2D not yet implemented" << std::endl; + } + + // // Case with ranges, calculate integral explicitly + // double xmin = x.min(rangeName) ; + // double xmax = x.max(rangeName) ; + // double ymin = y.min(rangeName) ; + // double ymax = y.max(rangeName) ; + + // double sum=0 ; + // int i ; + // for (i=1 ; i<=N_bins_x ; i++) { + // double binVal = (static_cast(pars.at(i-1))->getVal())/widths[i-1]; + // if (bins[i-1]>=xmin && bins[i]<=xmax) { + // // Bin fully in the integration domain + // sum += (bins[i]-bins[i-1])*binVal ; + // } else if (bins[i-1]xmax) { + // // Domain is fully contained in this bin + // sum += (xmax-xmin)*binVal ; + // // Exit here, this is the last bin to be processed by construction + // return sum/getFullSum() ; + // } else if (bins[i-1]xmin) { + // // Lower domain boundary is in bin + // sum += (bins[i]-xmin)*binVal ; + // } else if (bins[i-1]>=xmin && bins[i]>xmax && bins[i-1] + +ClassImp(chebyshevBasis) + +chebyshevBasis::chebyshevBasis( const char *name, + const char *title, + const Int_t _orderX, + const Int_t _orderY, + const TH2 &_binning ) : + orderX(_orderX), orderY(_orderY) + +{ + std::string thisName; + binning = static_cast(_binning.Clone()); + // Build list of RooRealVar Coefficients + for (Int_t ix=0; ix <= orderX; ix++){ + for (Int_t iy=0; iy <= orderY; iy++){ + thisName = "ChebCoeff_x"+std::to_string(ix)+"y"+std::to_string(iy); + RooRealVar* thisRRV = new RooRealVar(thisName.c_str(),thisName.c_str(),0.01,-10.0,10.0); + coefList.add(*thisRRV); + } + } + + + // Make the var map variables + xmin = _binning.GetXaxis()->GetXmin(); + xmax = _binning.GetXaxis()->GetXmax(); + ymin = _binning.GetYaxis()->GetXmin(); + ymax = _binning.GetYaxis()->GetXmax(); + slope_x = 2/(xmax-xmin); + slope_y = 2/(ymax-ymin); + +} + +RooAddition chebyshevBasis::getBinVal(const float xCenter, const float yCenter) const { + // Map "real" axis values to [-1,1] range + std::pair mappedvals = mapToChebyshev(xCenter,yCenter); + + Int_t xbin = binning->GetXaxis()->FindBin(xCenter); + Int_t ybin = binning->GetYaxis()->FindBin(yCenter); + + std::string basisName; + std::string coeffName; + std::string formulaName; + + RooArgList toSum("toSum"); + for (Int_t ix=0; ix<=orderX; ix++){ + for (Int_t iy=0; iy<=orderY; iy++){ + // std::cout << "DEBUG: " << ix << ", " << iy << std::endl; + // Bin value + basisName = "ChebBasisVal_x"+std::to_string(ix)+"y"+std::to_string(iy)+"_bin"+std::to_string(xbin)+"-"+std::to_string(ybin); + RooConstVar* basisConst = new RooConstVar(basisName.c_str(),basisName.c_str(),Eval2DChebyshev(mappedvals.first,mappedvals.second,ix,iy)); + // Coefficient + coeffName = "ChebCoeff_x"+std::to_string(ix)+"y"+std::to_string(iy); + RooAbsArg* coeff = coefList.find(coeffName.c_str()); + + RooArgList* formulaInput = new RooArgList(*coeff,*basisConst); + + formulaName = "ChebFormula_x"+std::to_string(ix)+"y"+std::to_string(iy)+"_bin"+std::to_string(xbin)+"-"+std::to_string(ybin); + RooFormulaVar* thisCheb = new RooFormulaVar(formulaName.c_str(),formulaName.c_str(), + "@0*@1+abs(@0)", + *formulaInput); // a_ij * T_i(x)*T_j(y) + abs(a_ij) + + // std::cout << "DEBUG: " << formulaName << std::endl; + + toSum.add(*thisCheb); + } + } + + std::string thisName = "ChebSum_"+std::to_string(xbin)+"-"+std::to_string(ybin); + RooAddition binVal(thisName.c_str(),thisName.c_str(),toSum); + + return binVal; +} + +float chebyshevBasis::Eval2DChebyshev(float x, float y, int thisOrderX, int thisOrderY) const{ + // Initialize polynomials + float Tx_n2 = 0; // T_n-2 + float Tx_n1 = 0; // T_n-1 + float Tx_n = 0; // T_n + float Ty_n2 = 0; // T_n-2 + float Ty_n1 = 0; // T_n-1 + float Ty_n = 0; // T_n + float Tx = 0; + float Ty = 0; + + // Construct polynomial in x + for(Int_t i=0; i<=thisOrderX; ++i){ + // T_0 == 1 + if (i == 0) { + Tx_n = 1; + } + // T_1 = x + else if (i == 1) { + Tx_n1 = 1; + Tx_n = x; + } + // T_n = 2*x*T_n-1 - T_n-2 + else { + Tx_n2 = Tx_n1; + Tx_n1 = Tx_n; + Tx_n = 2*x*Tx_n1-Tx_n2; + } + } + + // Construct polynomial in y + for (Int_t j=0; j<=thisOrderY; ++j){ + // T_0 == 1 + if (j == 0) { + Ty_n = 1; + } + // T_1 = x + else if (j == 1) { + Ty_n1 = 1; + Ty_n = y; + } + // T_n = 2*x*T_n-1 - T_n-2 + else { + Ty_n2 = Ty_n1; + Ty_n1 = Ty_n; + Ty_n = 2*y*Ty_n1-Ty_n2; + } + } + // Only want to shift if order is > 0 + Tx = Tx_n; + Ty = Ty_n; + + return Tx*Ty; +} + + +std::pair chebyshevBasis::mapToChebyshev(float ix,float iy) const { + + float xp = slope_x*(ix-xmax)+1; // Map x + float yp = slope_y*(iy-ymax)+1; // Map y + + return std::make_pair(xp,yp); +} + +void chebyshevBasis::drawBasis(std::string file_name) { + + TFile outfile(file_name.c_str(),"RECREATE"); + + TCanvas basisCan("basisCan","basisCan",800,700); + std::string canName; + outfile.cd(); + + std::string polyKey; + std::string cheb_name; + + for (Int_t oX=0; oX<=orderX; ++oX){ + for (Int_t oY=0; oY<=orderY; ++oY) { + polyKey = "x"+std::to_string(oX)+"y"+std::to_string(oY); + + cheb_name = "cheb_Tx"+std::to_string(oX)+"_Ty"+std::to_string(oY); + const char *cheb_name_cstr = cheb_name.c_str(); + TH2F* cheb2D = new TH2F(cheb_name_cstr, cheb_name_cstr, 100, -1.0, 1.0, 100, -1.0, 1.0); + + // Loop over TH2F bins and evaluate + for (Int_t xbin=1; xbin<=cheb2D->GetNbinsX(); xbin++){ + for (Int_t ybin=1; ybin<=cheb2D->GetNbinsY(); ybin++){ + float xcenter = cheb2D->GetXaxis()->GetBinCenter(xbin); + float ycenter = cheb2D->GetYaxis()->GetBinCenter(ybin); + float val = Eval2DChebyshev(xcenter,ycenter,oX,oY); + + cheb2D->SetBinContent(xbin,ybin,val); + } + } + + cheb2D->Draw("surf"); + cheb2D->Write(); + + canName = "basis_plots/x"+std::to_string(oX)+"y"+std::to_string(oY)+".png"; + basisCan.Print(canName.c_str(),"png"); + } + } + + outfile.Close(); +} \ No newline at end of file diff --git a/src/classes.h b/src/classes.h index e7b434770a5..9275e79751a 100644 --- a/src/classes.h +++ b/src/classes.h @@ -46,6 +46,8 @@ //#include "HiggsAnalysis/CombinedLimit/interface/RooMomentMorphND.h" #include "HiggsAnalysis/CombinedLimit/interface/RooMorphingPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/RooParametricHist.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooParametricHist2D.h" +#include "HiggsAnalysis/CombinedLimit/interface/chebyshevBasis.h" #include "HiggsAnalysis/CombinedLimit/interface/RooParametricShapeBinPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/GaussExp.h" #include "HiggsAnalysis/CombinedLimit/interface/RooDoubleCBFast.h" diff --git a/src/classes_def.xml b/src/classes_def.xml index 490ed437aac..326a2c369f3 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -193,6 +193,8 @@ + + From b539fb84e4b4613c8d497b09c9c8bc0360eb3408 Mon Sep 17 00:00:00 2001 From: lcorcodilos Date: Fri, 17 May 2019 12:57:49 -0400 Subject: [PATCH 2/5] Remove chebyshevBasis --- interface/chebyshevBasis.h | 54 ----------- src/chebyshevBasis.cxx | 180 ------------------------------------- src/classes.h | 1 - src/classes_def.xml | 1 - 4 files changed, 236 deletions(-) delete mode 100644 interface/chebyshevBasis.h delete mode 100644 src/chebyshevBasis.cxx diff --git a/interface/chebyshevBasis.h b/interface/chebyshevBasis.h deleted file mode 100644 index c085ef437dd..00000000000 --- a/interface/chebyshevBasis.h +++ /dev/null @@ -1,54 +0,0 @@ -#ifndef CHEBYSHEVBASIS -#define CHEBYSHEVBASIS - -#include "TH2F.h" -#include "TCanvas.h" -#include "TFile.h" -#include "RooFit.h" -#include "RooFormulaVar.h" -#include "RooConstVar.h" -#include "RooArgList.h" -#include "RooAddition.h" -#include "RooAbsArg.h" -#include "Riostream.h" -#include "RooRealVar.h" -#include - -class chebyshevBasis : public TObject { -public: - chebyshevBasis() {}; - chebyshevBasis( const char *name, - const char *title, - const Int_t _orderX, - const Int_t _orderY, - const TH2 &_binning ); - inline virtual ~chebyshevBasis (){}; - // float getChebyshevVal(float xval, float yval, float thisOrderX, float thisOrderY) const; - RooAddition getBinVal(const float xCenter, const float yCenter) const; - void drawBasis(std::string file_name); - -protected: - // variables - TH2F* binning; - Int_t orderX; - Int_t orderY; - RooArgList coefList; - Double_t xmin; - Double_t xmax; - Double_t ymin; - Double_t ymax; - float slope_x; - float slope_y; - // std::map polyHists; - - // functions - std::pair mapToChebyshev(float ix,float iy) const; - float Eval2DChebyshev(float x, float y, int thisOrderX, int thisOrderY) const; - // TH2F Make2DChebyshev(const int& thisOrderX, const int& thisOrderY) const; - - -private: - ClassDef(chebyshevBasis, 1) -}; - -#endif \ No newline at end of file diff --git a/src/chebyshevBasis.cxx b/src/chebyshevBasis.cxx deleted file mode 100644 index c9d7cbdbda4..00000000000 --- a/src/chebyshevBasis.cxx +++ /dev/null @@ -1,180 +0,0 @@ -#include "HiggsAnalysis/CombinedLimit/interface/chebyshevBasis.h" -#include - -ClassImp(chebyshevBasis) - -chebyshevBasis::chebyshevBasis( const char *name, - const char *title, - const Int_t _orderX, - const Int_t _orderY, - const TH2 &_binning ) : - orderX(_orderX), orderY(_orderY) - -{ - std::string thisName; - binning = static_cast(_binning.Clone()); - // Build list of RooRealVar Coefficients - for (Int_t ix=0; ix <= orderX; ix++){ - for (Int_t iy=0; iy <= orderY; iy++){ - thisName = "ChebCoeff_x"+std::to_string(ix)+"y"+std::to_string(iy); - RooRealVar* thisRRV = new RooRealVar(thisName.c_str(),thisName.c_str(),0.01,-10.0,10.0); - coefList.add(*thisRRV); - } - } - - - // Make the var map variables - xmin = _binning.GetXaxis()->GetXmin(); - xmax = _binning.GetXaxis()->GetXmax(); - ymin = _binning.GetYaxis()->GetXmin(); - ymax = _binning.GetYaxis()->GetXmax(); - slope_x = 2/(xmax-xmin); - slope_y = 2/(ymax-ymin); - -} - -RooAddition chebyshevBasis::getBinVal(const float xCenter, const float yCenter) const { - // Map "real" axis values to [-1,1] range - std::pair mappedvals = mapToChebyshev(xCenter,yCenter); - - Int_t xbin = binning->GetXaxis()->FindBin(xCenter); - Int_t ybin = binning->GetYaxis()->FindBin(yCenter); - - std::string basisName; - std::string coeffName; - std::string formulaName; - - RooArgList toSum("toSum"); - for (Int_t ix=0; ix<=orderX; ix++){ - for (Int_t iy=0; iy<=orderY; iy++){ - // std::cout << "DEBUG: " << ix << ", " << iy << std::endl; - // Bin value - basisName = "ChebBasisVal_x"+std::to_string(ix)+"y"+std::to_string(iy)+"_bin"+std::to_string(xbin)+"-"+std::to_string(ybin); - RooConstVar* basisConst = new RooConstVar(basisName.c_str(),basisName.c_str(),Eval2DChebyshev(mappedvals.first,mappedvals.second,ix,iy)); - // Coefficient - coeffName = "ChebCoeff_x"+std::to_string(ix)+"y"+std::to_string(iy); - RooAbsArg* coeff = coefList.find(coeffName.c_str()); - - RooArgList* formulaInput = new RooArgList(*coeff,*basisConst); - - formulaName = "ChebFormula_x"+std::to_string(ix)+"y"+std::to_string(iy)+"_bin"+std::to_string(xbin)+"-"+std::to_string(ybin); - RooFormulaVar* thisCheb = new RooFormulaVar(formulaName.c_str(),formulaName.c_str(), - "@0*@1+abs(@0)", - *formulaInput); // a_ij * T_i(x)*T_j(y) + abs(a_ij) - - // std::cout << "DEBUG: " << formulaName << std::endl; - - toSum.add(*thisCheb); - } - } - - std::string thisName = "ChebSum_"+std::to_string(xbin)+"-"+std::to_string(ybin); - RooAddition binVal(thisName.c_str(),thisName.c_str(),toSum); - - return binVal; -} - -float chebyshevBasis::Eval2DChebyshev(float x, float y, int thisOrderX, int thisOrderY) const{ - // Initialize polynomials - float Tx_n2 = 0; // T_n-2 - float Tx_n1 = 0; // T_n-1 - float Tx_n = 0; // T_n - float Ty_n2 = 0; // T_n-2 - float Ty_n1 = 0; // T_n-1 - float Ty_n = 0; // T_n - float Tx = 0; - float Ty = 0; - - // Construct polynomial in x - for(Int_t i=0; i<=thisOrderX; ++i){ - // T_0 == 1 - if (i == 0) { - Tx_n = 1; - } - // T_1 = x - else if (i == 1) { - Tx_n1 = 1; - Tx_n = x; - } - // T_n = 2*x*T_n-1 - T_n-2 - else { - Tx_n2 = Tx_n1; - Tx_n1 = Tx_n; - Tx_n = 2*x*Tx_n1-Tx_n2; - } - } - - // Construct polynomial in y - for (Int_t j=0; j<=thisOrderY; ++j){ - // T_0 == 1 - if (j == 0) { - Ty_n = 1; - } - // T_1 = x - else if (j == 1) { - Ty_n1 = 1; - Ty_n = y; - } - // T_n = 2*x*T_n-1 - T_n-2 - else { - Ty_n2 = Ty_n1; - Ty_n1 = Ty_n; - Ty_n = 2*y*Ty_n1-Ty_n2; - } - } - // Only want to shift if order is > 0 - Tx = Tx_n; - Ty = Ty_n; - - return Tx*Ty; -} - - -std::pair chebyshevBasis::mapToChebyshev(float ix,float iy) const { - - float xp = slope_x*(ix-xmax)+1; // Map x - float yp = slope_y*(iy-ymax)+1; // Map y - - return std::make_pair(xp,yp); -} - -void chebyshevBasis::drawBasis(std::string file_name) { - - TFile outfile(file_name.c_str(),"RECREATE"); - - TCanvas basisCan("basisCan","basisCan",800,700); - std::string canName; - outfile.cd(); - - std::string polyKey; - std::string cheb_name; - - for (Int_t oX=0; oX<=orderX; ++oX){ - for (Int_t oY=0; oY<=orderY; ++oY) { - polyKey = "x"+std::to_string(oX)+"y"+std::to_string(oY); - - cheb_name = "cheb_Tx"+std::to_string(oX)+"_Ty"+std::to_string(oY); - const char *cheb_name_cstr = cheb_name.c_str(); - TH2F* cheb2D = new TH2F(cheb_name_cstr, cheb_name_cstr, 100, -1.0, 1.0, 100, -1.0, 1.0); - - // Loop over TH2F bins and evaluate - for (Int_t xbin=1; xbin<=cheb2D->GetNbinsX(); xbin++){ - for (Int_t ybin=1; ybin<=cheb2D->GetNbinsY(); ybin++){ - float xcenter = cheb2D->GetXaxis()->GetBinCenter(xbin); - float ycenter = cheb2D->GetYaxis()->GetBinCenter(ybin); - float val = Eval2DChebyshev(xcenter,ycenter,oX,oY); - - cheb2D->SetBinContent(xbin,ybin,val); - } - } - - cheb2D->Draw("surf"); - cheb2D->Write(); - - canName = "basis_plots/x"+std::to_string(oX)+"y"+std::to_string(oY)+".png"; - basisCan.Print(canName.c_str(),"png"); - } - } - - outfile.Close(); -} \ No newline at end of file diff --git a/src/classes.h b/src/classes.h index 9275e79751a..be30252e0b7 100644 --- a/src/classes.h +++ b/src/classes.h @@ -47,7 +47,6 @@ #include "HiggsAnalysis/CombinedLimit/interface/RooMorphingPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/RooParametricHist.h" #include "HiggsAnalysis/CombinedLimit/interface/RooParametricHist2D.h" -#include "HiggsAnalysis/CombinedLimit/interface/chebyshevBasis.h" #include "HiggsAnalysis/CombinedLimit/interface/RooParametricShapeBinPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/GaussExp.h" #include "HiggsAnalysis/CombinedLimit/interface/RooDoubleCBFast.h" diff --git a/src/classes_def.xml b/src/classes_def.xml index 326a2c369f3..cac938caa80 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -194,7 +194,6 @@ - From 267e2db0a6b40b4ca8d1c4b9864ff502db59dc19 Mon Sep 17 00:00:00 2001 From: Nick Smith Date: Wed, 6 Apr 2022 09:56:44 -0500 Subject: [PATCH 3/5] Use log(N) search --- src/RooParametricHist.cxx | 17 +++++----- src/RooParametricHist2D.cxx | 66 ++++++++++++++----------------------- 2 files changed, 34 insertions(+), 49 deletions(-) diff --git a/src/RooParametricHist.cxx b/src/RooParametricHist.cxx index 80de96fbc27..6df0319d3a2 100644 --- a/src/RooParametricHist.cxx +++ b/src/RooParametricHist.cxx @@ -210,15 +210,16 @@ double RooParametricHist::evaluateMorphFunction(int j) const } double RooParametricHist::evaluatePartial() const { - int bin_i; - if (x < bins[0]) return 0; // should set to 0 instead? - else if (x >= bins[N_bins]) return 0; - - else { - for(bin_i=0; bin_i=bins[bin_i] && x < bins[bin_i+1] ) break; - } + auto it = std::upper_bound(std::begin(bins), std::end(bins), x); + if ( it == std::begin(bins) ) { + // underflow + return 0; + } + else if ( it == std::end(bins) ) { + // overflow + return 0; } + size_t bin_i = std::distance(std::begin(bins), it) - 1; RooAbsReal *retVar = (RooAbsReal*)pars.at(bin_i); double ret = retVar->getVal(); diff --git a/src/RooParametricHist2D.cxx b/src/RooParametricHist2D.cxx index 9fb10aaa0d8..df3c85d62d1 100644 --- a/src/RooParametricHist2D.cxx +++ b/src/RooParametricHist2D.cxx @@ -201,50 +201,34 @@ Double_t RooParametricHist2D::analyticalIntegral(Int_t code, const char* rangeNa Double_t RooParametricHist2D::evaluate() const { + auto itx = std::upper_bound(std::begin(bins_x), std::end(bins_x), x); + if ( itx == std::begin(bins_x) ) { + // underflow + return 0; + } + else if ( itx == std::end(bins_x) ) { + // overflow + return 0; + } + auto ity = std::upper_bound(std::begin(bins_y), std::end(bins_y), y); + if ( ity == std::begin(bins_y) ) { + // underflow + return 0; + } + else if ( ity == std::end(bins_y) ) { + // overflow + return 0; + } + size_t bin_ix = std::distance(std::begin(bins_x), itx) - 1; + size_t bin_iy = std::distance(std::begin(bins_y), ity) - 1; - int bin_ix; - int bin_iy; - - // std::cout << "Evaluating..." << std::endl; - // std::cout << "bins_x size = " << bins_x.size() << std::endl; - // std::cout << "bins_x type = " << typeid(bins_x).name() << std::endl; - // for (auto i = bins_x.begin(); i != bins_x.end(); ++i){ - // std::cout << *i << ' '; - // } - // std::cout << std::endl; - - // std::cout << "bins_y size = " << bins_y.size() << std::endl; - // std::cout << "bins_y type = " << typeid(bins_y).name() << std::endl; - // for (auto i = bins_y.begin(); i != bins_y.end(); ++i){ - // std::cout << *i << ' '; - // } - // std::cout << std::endl; - - if (x < bins_x[0]) { - return 0; // should set to 0 instead? - } else if (y < bins_y[0]) { - return 0; - } else if (x >= bins_x[N_bins_x]) { - return 0; - } else if (y >= bins_y[N_bins_y]) { - return 0; - } else { - for(bin_ix=0; bin_ix=bins_x[bin_ix] && x < bins_x[bin_ix+1] ) break; - } - for(bin_iy=0; bin_iy=bins_y[bin_iy] && y < bins_y[bin_iy+1] ) break; - } - } - - int globalbin = N_bins_x * bin_iy + bin_ix; + int globalbin = N_bins_x * bin_iy + bin_ix; - RooAbsReal *retVar = (RooAbsReal*)pars.at(globalbin); + RooAbsReal *retVar = (RooAbsReal*)pars.at(globalbin); - double ret = retVar->getVal() / (widths_x[bin_ix]*widths_y[bin_iy]); + double ret = retVar->getVal() / (widths_x[bin_ix]*widths_y[bin_iy]); - cval=ret; - // std::cout << "Evaluated" << std::endl; - return ret; + cval=ret; + return ret; } From 8aafddae43bc279fd41f5eb692e33d51e308fd21 Mon Sep 17 00:00:00 2001 From: Nick Smith Date: Wed, 6 Apr 2022 11:04:39 -0500 Subject: [PATCH 4/5] Remove extraneous comments --- python/ShapeTools.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/python/ShapeTools.py b/python/ShapeTools.py index 4b6bd3c659e..b33a9edf1c3 100644 --- a/python/ShapeTools.py +++ b/python/ShapeTools.py @@ -432,7 +432,6 @@ def prepareAllShapes(self): if self.options.verbose > 1: stderr.write("Observables: %s\n" % str(shapeObs.keys())) if len(shapeObs.keys()) != 1: self.out.binVars = ROOT.RooArgSet() - # Custom code starts - LC 5/14/19 for obs_key in shapeObs.keys(): # RooArgSet does not have method to add another RooArgSet (RooCollection does!) so we have to manually loop over the RooArgSet contents (in obs) if ',' in obs_key: # if multiple vars and looking at a RooArgSet value @@ -440,7 +439,6 @@ def prepareAllShapes(self): self.out.binVars.add(shapeObs[obs_key].find(k), True) else: self.out.binVars.add(shapeObs[obs_key], True) - # Custom code ends else: self.out.binVars = shapeObs.values()[0] self.out._import(self.out.binVars) @@ -706,7 +704,6 @@ def getPdf(self,channel,process,_cache={}): pdfs.Add(self.shape2Pdf(shapeDown,channel,process)) histpdf = nominalPdf if nominalPdf.InheritsFrom("RooDataHist") else nominalPdf.dataHist() xvar = histpdf.get().first() - # Custom code start - LC 5/14/19 histpdfSubset = histpdf.get().Clone() histpdfSubset.remove(xvar) if histpdfSubset.first() != False: @@ -714,8 +711,6 @@ def getPdf(self,channel,process,_cache={}): rhp = ROOT.FastVerticalInterpHistPdf2D2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, yvar, False, pdfs, coeffs, qrange, qalgo) else: rhp = ROOT.FastVerticalInterpHistPdf2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, pdfs, coeffs, qrange, qalgo) - # Custom code ends (plus the next line was commented out) - #rhp = ROOT.FastVerticalInterpHistPdf2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, pdfs, coeffs, qrange, qalgo) _cache[(channel,process)] = rhp return rhp elif nominalPdf.InheritsFrom("RooParametricHist") : From 2648c699f4e0f20054a75f1ed117513b47eaebcf Mon Sep 17 00:00:00 2001 From: Nick Smith Date: Wed, 6 Apr 2022 15:14:45 -0500 Subject: [PATCH 5/5] Correct handling of 1D RooParametricHist workspaces --- python/ShapeTools.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/python/ShapeTools.py b/python/ShapeTools.py index b33a9edf1c3..e789f005dfc 100644 --- a/python/ShapeTools.py +++ b/python/ShapeTools.py @@ -703,11 +703,12 @@ def getPdf(self,channel,process,_cache={}): pdfs.Add(self.shape2Pdf(shapeUp,channel,process)) pdfs.Add(self.shape2Pdf(shapeDown,channel,process)) histpdf = nominalPdf if nominalPdf.InheritsFrom("RooDataHist") else nominalPdf.dataHist() - xvar = histpdf.get().first() - histpdfSubset = histpdf.get().Clone() - histpdfSubset.remove(xvar) - if histpdfSubset.first() != False: - yvar = histpdfSubset.first() # y is a little more complex + varIter = histpdf.get().createIterator() + xvar = varIter.Next() + yvar = varIter.Next() + if varIter.Next(): + raise ValueError("No support for 3+ dimensional histpdfs") + elif yvar: rhp = ROOT.FastVerticalInterpHistPdf2D2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, yvar, False, pdfs, coeffs, qrange, qalgo) else: rhp = ROOT.FastVerticalInterpHistPdf2("shape%s_%s_%s_morph" % (postFix,channel,process), "", xvar, pdfs, coeffs, qrange, qalgo)