From f27aa021dc52d7d7efb5766ae45618eec8ab0bb3 Mon Sep 17 00:00:00 2001 From: nckw Date: Thu, 14 Dec 2023 08:58:48 +0000 Subject: [PATCH] Remove overlaps with IvyFramework (#878) Co-authored-by: Nick Wardle --- interface/AsymQuad.h | 50 --- interface/RooFuncPdf.h | 44 --- interface/RooNCSplineCore.h | 108 ------ interface/RooNCSplineFactory_1D.h | 65 ---- interface/RooNCSplineFactory_2D.h | 101 ----- interface/RooNCSplineFactory_3D.h | 104 ----- interface/RooNCSpline_1D_fast.h | 68 ---- interface/RooNCSpline_2D_fast.h | 84 ----- interface/RooNCSpline_3D_fast.h | 105 ------ interface/RooPiecewisePolynomial.h | 50 --- interface/RooRealFlooredSumPdf.h | 95 ----- src/AsymQuad.cc | 128 ------- src/RooFuncPdf.cc | 3 - src/RooNCSplineCore.cc | 284 -------------- src/RooNCSplineFactory_1D.cc | 88 ----- src/RooNCSplineFactory_2D.cc | 124 ------ src/RooNCSplineFactory_3D.cc | 146 -------- src/RooNCSpline_1D_fast.cc | 226 ----------- src/RooNCSpline_2D_fast.cc | 394 ------------------- src/RooNCSpline_3D_fast.cc | 550 --------------------------- src/RooPiecewisePolynomial.cc | 215 ----------- src/RooRealFlooredSumPdf.cc | 583 ----------------------------- src/classes.h | 8 - src/classes_def.xml | 9 - 24 files changed, 3632 deletions(-) delete mode 100644 interface/AsymQuad.h delete mode 100644 interface/RooFuncPdf.h delete mode 100644 interface/RooNCSplineCore.h delete mode 100644 interface/RooNCSplineFactory_1D.h delete mode 100644 interface/RooNCSplineFactory_2D.h delete mode 100644 interface/RooNCSplineFactory_3D.h delete mode 100644 interface/RooNCSpline_1D_fast.h delete mode 100644 interface/RooNCSpline_2D_fast.h delete mode 100644 interface/RooNCSpline_3D_fast.h delete mode 100644 interface/RooPiecewisePolynomial.h delete mode 100755 interface/RooRealFlooredSumPdf.h delete mode 100644 src/AsymQuad.cc delete mode 100644 src/RooFuncPdf.cc delete mode 100644 src/RooNCSplineCore.cc delete mode 100644 src/RooNCSplineFactory_1D.cc delete mode 100644 src/RooNCSplineFactory_2D.cc delete mode 100644 src/RooNCSplineFactory_3D.cc delete mode 100644 src/RooNCSpline_1D_fast.cc delete mode 100644 src/RooNCSpline_2D_fast.cc delete mode 100644 src/RooNCSpline_3D_fast.cc delete mode 100644 src/RooPiecewisePolynomial.cc delete mode 100755 src/RooRealFlooredSumPdf.cc diff --git a/interface/AsymQuad.h b/interface/AsymQuad.h deleted file mode 100644 index 56dee3d63cd..00000000000 --- a/interface/AsymQuad.h +++ /dev/null @@ -1,50 +0,0 @@ -#ifndef HiggsAnalysis_CombinedLimit_AsymQuad_h -#define HiggsAnalysis_CombinedLimit_AsymQuad_h - -#include "RooFit.h" -#include "Riostream.h" -#include "TIterator.h" -#include "TList.h" -#include -#include "RooRealVar.h" -#include -#include "RooListProxy.h" -#include "RooMsgService.h" - - -//_________________________________________________ -/* -BEGIN_HTML -

-AsymQuad is helper class for implementing asymmetric additive interpolation. -

-END_HTML -*/ -// -class AsymQuad : public RooAbsReal { - -public: - AsymQuad(); - AsymQuad(const char *name, const char *title, const RooArgList& inFuncList, const RooArgList& inCoefList, Double_t smoothRegion=1., Int_t smoothAlgo=0); - AsymQuad(const AsymQuad& other, const char* name=0); - ~AsymQuad(); - - TObject* clone(const char* newname) const { return new AsymQuad(*this, newname); } - -protected: - Double_t evaluate() const; - - RooListProxy _funcList; // List of component functions - RooListProxy _coefList; // List of coefficients - Double_t smoothRegion_; - Int_t smoothAlgo_; - TIterator* _funcIter; //! Iterator over FUNC list - TIterator* _coefIter; //! Iterator over coefficient list - -private: - Double_t interpolate(Double_t theta_, Double_t valueCenter_, Double_t valueHigh_, Double_t valueLow_) const; - - ClassDef(AsymQuad, 1) // Asymmetric power -}; - -#endif diff --git a/interface/RooFuncPdf.h b/interface/RooFuncPdf.h deleted file mode 100644 index 96aab546912..00000000000 --- a/interface/RooFuncPdf.h +++ /dev/null @@ -1,44 +0,0 @@ -#ifndef RooFUNCPDF -#define RooFUNCPDF - -#include "RooRealProxy.h" -#include "RooAbsPdf.h" - -class RooFuncPdf : public RooAbsPdf{ -protected: - RooRealProxy theFunc; - -public: - RooFuncPdf() : RooAbsPdf(){} - RooFuncPdf( - const char* name, - const char* title - ) : RooAbsPdf(name, title), theFunc("theFunc","theFunc",this){} - RooFuncPdf( - const char* name, - const char* title, - RooAbsReal& inFunc - ) : RooAbsPdf(name, title), theFunc("theFunc", "theFunc", this, inFunc){} - RooFuncPdf(const RooFuncPdf& other, const char* name=0) : RooAbsPdf(other, name), theFunc("theFunc", this, other.theFunc){} - TObject* clone(const char* newname)const{ return new RooFuncPdf(*this, newname); } - inline virtual ~RooFuncPdf(){} - - Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0)const{ return dynamic_cast(theFunc.absArg())->getAnalyticalIntegral(allVars, analVars, rangeName); } - Double_t analyticalIntegral(Int_t code, const char* rangeName=0)const{ return dynamic_cast(theFunc.absArg())->analyticalIntegral(code, rangeName); } - - Int_t getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet *normSet, const char* rangeName=0)const{ return dynamic_cast(theFunc.absArg())->getAnalyticalIntegralWN(allVars, analVars, normSet, rangeName); } - Double_t analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName=0)const{ - RooAbsPdf* pdfcast = dynamic_cast(theFunc.absArg()); - if (pdfcast!=0) return pdfcast->analyticalIntegralWN(code, normSet, rangeName); - else return RooAbsPdf::analyticalIntegralWN(code, normSet, rangeName); - } - -protected: - Double_t evaluate()const{ return theFunc; } - - - ClassDef(RooFuncPdf, 0) - -}; - -#endif diff --git a/interface/RooNCSplineCore.h b/interface/RooNCSplineCore.h deleted file mode 100644 index ead0aa44c55..00000000000 --- a/interface/RooNCSplineCore.h +++ /dev/null @@ -1,108 +0,0 @@ -#ifndef ROONCSPLINECORE -#define ROONCSPLINECORE - -#include -#include "Accumulators.h" -#include "TMatrixD.h" -#include "TVectorD.h" -#include "RooAbsReal.h" -#include "RooRealVar.h" -#include "RooRealProxy.h" -#include "RooConstVar.h" -#include "RooArgList.h" -#include "RooMsgService.h" -#include "RooListProxy.h" - -class RooNCSplineCore : public RooAbsReal{ -public: - typedef Float_t T; - typedef TMatrixT TMatrix_t; - typedef TVectorT TVector_t; - - enum VerbosityLevel{ - kSilent, - kError, - kVerbose - }; - - enum BoundaryCondition{ - bcApproximatedSlope, // D1 is the same as deltaY/deltaX in the first/last segment - bcClamped, // D1=0 at endpoint - bcApproximatedSecondDerivative, // D2 is approximated - bcNaturalSpline, // D2=0 at endpoint - bcQuadratic, // Coefficient D=0 - bcQuadraticWithNullSlope, // Coefficient D=0 and D1=0 at endpoint - NBoundaryConditions - }; - - RooNCSplineCore(); - RooNCSplineCore( - const char* name, - const char* title - ); - RooNCSplineCore( - const char* name, - const char* title, - RooAbsReal& inXVar, - const std::vector& inXList, - Bool_t inUseFloor=true, - T inFloorEval=1e-15, - T inFloorInt=1e-10 - ); - RooNCSplineCore(const RooNCSplineCore& other, const char* name=0); - virtual TObject* clone(const char* newname)const = 0; - inline virtual ~RooNCSplineCore(){} - - virtual void setVerbosity(VerbosityLevel flag); - void setEvalFloor(T val); - void setIntFloor(T val); - void doFloor(Bool_t flag); - - virtual void setRangeValidity(const T valmin, const T valmax, const Int_t whichDirection) = 0; - - virtual Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0)const = 0; - virtual Double_t analyticalIntegral(Int_t code, const char* rangeName=0)const = 0; - -protected: - VerbosityLevel verbosity; - Bool_t useFloor; - T floorEval; - T floorInt; - - T rangeXmin; - T rangeXmax; - - RooRealProxy theXVar; - RooListProxy leafDepsList; - std::vector XList; - - virtual void emptyFcnList() = 0; - - void getLeafDependents(RooRealProxy& proxy, RooArgSet& set); - void addLeafDependents(RooArgSet& set); - - unsigned int npointsX()const{ return XList.size(); } - - virtual Int_t getWhichBin(const T& val, const Int_t whichDirection)const = 0; - virtual void getKappas(std::vector& kappas, const Int_t whichDirection) = 0; - virtual T getTVar(const std::vector& kappas, const T& val, const Int_t& bin, const Int_t whichDirection)const = 0; - - virtual Bool_t testRangeValidity(const T& val, const Int_t whichDirection)const = 0; - virtual void cropValueForRange(T& val, const Int_t whichDirection)const = 0; - - virtual T interpolateFcn(Int_t code, const char* rangeName=0)const = 0; - virtual Double_t evaluate()const = 0; - - virtual void getBArray(const std::vector& kappas, const std::vector& fcnList, std::vector& BArray, BoundaryCondition const& bcBegin, BoundaryCondition const& bcEnd)const; - virtual void getAArray(const std::vector& kappas, std::vector>& AArray, BoundaryCondition const& bcBegin, BoundaryCondition const& bcEnd)const; - virtual std::vector> getCoefficientsAlongDirection(const std::vector& kappas, const TMatrix_t& Ainv, const std::vector& fcnList, BoundaryCondition const& bcBegin, BoundaryCondition const& bcEnd, const Int_t pickBin)const; - virtual std::vector getCoefficients(const TVector_t& S, const std::vector& kappas, const std::vector& fcnList, const Int_t& bin)const; - - virtual T evalSplineSegment(const std::vector& coefs, const T& kappa, const T& tup, const T& tdn, Bool_t doIntegrate=false)const; - - - ClassDef(RooNCSplineCore, 3) - -}; - -#endif diff --git a/interface/RooNCSplineFactory_1D.h b/interface/RooNCSplineFactory_1D.h deleted file mode 100644 index bd1e2c61175..00000000000 --- a/interface/RooNCSplineFactory_1D.h +++ /dev/null @@ -1,65 +0,0 @@ -#ifndef ROONCSPLINEFACTORY_1D -#define ROONCSPLINEFACTORY_1D - -#include -#include -#include -#include "TGraph.h" -#include "TTree.h" -#include "RooNCSpline_1D_fast.h" -#include "RooFuncPdf.h" - -class RooNCSplineFactory_1D{ -protected: - TString appendName; - - RooNCSplineCore::BoundaryCondition bcBeginX; - RooNCSplineCore::BoundaryCondition bcEndX; - - RooAbsReal* splineVar; - RooNCSpline_1D_fast* fcn; - RooFuncPdf* PDF; - - const std::vector> getPoints(const std::vector& XList, const std::vector& FcnList); - - void destroyPDF(); - void initPDF(const std::vector>& pList); - -public: - RooNCSplineFactory_1D( - RooAbsReal& splineVar_, TString appendName_="", - RooNCSplineCore::BoundaryCondition const bcBeginX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndX_=RooNCSplineCore::bcNaturalSpline - ); - ~RooNCSplineFactory_1D(); - - RooNCSpline_1D_fast* getFunc(){ return fcn; } - RooFuncPdf* getPDF(){ return PDF; } - - void setEndConditions( - RooNCSplineCore::BoundaryCondition const bcBegin, - RooNCSplineCore::BoundaryCondition const bcEnd, - const unsigned int direction=0 - ); - - void setPoints(TTree* tree); - void setPoints(TGraph* tg); - void setPoints(const std::vector>& pList){ initPDF(pList); } - template void setPoints(const std::vector& XList, const std::vector& FcnList){ - std::vector transXList; - std::vector transFcnList; - for (unsigned int ip=0; ip> pList = getPoints(transXList, transFcnList); - initPDF(pList); - } - -}; - -template void RooNCSplineFactory_1D::setPoints(const std::vector& XList, const std::vector& FcnList); -template void RooNCSplineFactory_1D::setPoints(const std::vector& XList, const std::vector& FcnList); - -#endif - - - diff --git a/interface/RooNCSplineFactory_2D.h b/interface/RooNCSplineFactory_2D.h deleted file mode 100644 index 902c874c3df..00000000000 --- a/interface/RooNCSplineFactory_2D.h +++ /dev/null @@ -1,101 +0,0 @@ -#ifndef ROONCSPLINEFACTORY_2D -#define ROONCSPLINEFACTORY_2D - -#include -#include -#include -#include "TTree.h" -#include "RooNCSpline_2D_fast.h" -#include "RooFuncPdf.h" - - -namespace NumUtils{ - template struct triplet{ - T value[3]; - triplet(T i1, T i2, T i3){ - value[0]=i1; - value[1]=i2; - value[2]=i3; - - } - triplet(T i1){ - value[0]=i1; - value[1]=i1; - value[2]=i1; - - } - triplet(){} - T& operator[](std::size_t ipos){ return value[ipos]; } // Return by reference - const T& operator[](std::size_t ipos)const{ return value[ipos]; } // Return by const reference - }; - - typedef triplet intTriplet_t; - typedef triplet floatTriplet_t; - typedef triplet doubleTriplet_t; -} - -typedef NumUtils::triplet splineTriplet_t; - - -class RooNCSplineFactory_2D{ -protected: - TString appendName; - - RooNCSplineCore::BoundaryCondition bcBeginX; - RooNCSplineCore::BoundaryCondition bcEndX; - RooNCSplineCore::BoundaryCondition bcBeginY; - RooNCSplineCore::BoundaryCondition bcEndY; - - RooAbsReal* XVar; - RooAbsReal* YVar; - RooNCSpline_2D_fast* fcn; - RooFuncPdf* PDF; - - const std::vector getPoints(const std::vector& XList, const std::vector& YList, const std::vector& FcnList); - - void destroyPDF(); - void initPDF(const std::vector& pList); - - void addUnique(std::vector& list, RooNCSplineCore::T val); - -public: - RooNCSplineFactory_2D( - RooAbsReal& XVar_, RooAbsReal& YVar_, TString appendName_="", - RooNCSplineCore::BoundaryCondition const bcBeginX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcBeginY_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndY_=RooNCSplineCore::bcNaturalSpline - ); - ~RooNCSplineFactory_2D(); - - RooNCSpline_2D_fast* getFunc(){ return fcn; } - RooFuncPdf* getPDF(){ return PDF; } - - void setEndConditions( - RooNCSplineCore::BoundaryCondition const bcBegin, - RooNCSplineCore::BoundaryCondition const bcEnd, - const unsigned int direction - ); - - void setPoints(TTree* tree); - void setPoints(const std::vector& pList){ initPDF(pList); } - template void setPoints(const std::vector& XList, const std::vector& YList, const std::vector& FcnList){ - std::vector transXList; - std::vector transYList; - std::vector transFcnList; - for (unsigned int ip=0; ip pList = getPoints(transXList, transYList, transFcnList); - setPoints(pList); - } - -}; - -template void RooNCSplineFactory_2D::setPoints(const std::vector& XList, const std::vector& YList, const std::vector& FcnList); -template void RooNCSplineFactory_2D::setPoints(const std::vector& XList, const std::vector& YList, const std::vector& FcnList); - -#endif - - - diff --git a/interface/RooNCSplineFactory_3D.h b/interface/RooNCSplineFactory_3D.h deleted file mode 100644 index 1f0b68c0e4b..00000000000 --- a/interface/RooNCSplineFactory_3D.h +++ /dev/null @@ -1,104 +0,0 @@ -#ifndef ROONCSPLINEFACTORY_3D -#define ROONCSPLINEFACTORY_3D - -#include -#include -#include -#include "TTree.h" -#include "RooNCSpline_3D_fast.h" -#include "RooFuncPdf.h" - - -namespace NumUtils{ - template struct quadruplet{ - T value[4]; - quadruplet(T i1, T i2, T i3, T i4){ - value[0]=i1; - value[1]=i2; - value[2]=i3; - value[3]=i4; - } - quadruplet(T i1){ for (unsigned int idim=0; idim<4; idim++) value[idim] = i1; } - quadruplet(){} - T& operator[](std::size_t ipos){ return value[ipos]; } // Return by reference - const T& operator[](std::size_t ipos)const{ return value[ipos]; } // Return by const reference - }; - - typedef quadruplet intQuad_t; - typedef quadruplet floatQuad_t; - typedef quadruplet doubleQuad_t; -} - -typedef NumUtils::quadruplet splineQuadruplet_t; - - -class RooNCSplineFactory_3D{ -protected: - TString appendName; - - RooNCSplineCore::BoundaryCondition bcBeginX; - RooNCSplineCore::BoundaryCondition bcEndX; - RooNCSplineCore::BoundaryCondition bcBeginY; - RooNCSplineCore::BoundaryCondition bcEndY; - RooNCSplineCore::BoundaryCondition bcBeginZ; - RooNCSplineCore::BoundaryCondition bcEndZ; - - RooAbsReal* XVar; - RooAbsReal* YVar; - RooAbsReal* ZVar; - RooNCSpline_3D_fast* fcn; - RooFuncPdf* PDF; - - const std::vector getPoints(const std::vector& XList, const std::vector& YList, const std::vector& ZList, const std::vector& FcnList); - - void destroyPDF(); - void initPDF(const std::vector& pList); - - void addUnique(std::vector& list, RooNCSplineCore::T val); - -public: - RooNCSplineFactory_3D( - RooAbsReal& XVar_, RooAbsReal& YVar_, RooAbsReal& ZVar_, TString appendName_="", - RooNCSplineCore::BoundaryCondition const bcBeginX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcBeginY_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndY_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcBeginZ_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndZ_=RooNCSplineCore::bcNaturalSpline - ); - ~RooNCSplineFactory_3D(); - - RooNCSpline_3D_fast* getFunc(){ return fcn; } - RooFuncPdf* getPDF(){ return PDF; } - - void setEndConditions( - RooNCSplineCore::BoundaryCondition const bcBegin, - RooNCSplineCore::BoundaryCondition const bcEnd, - const unsigned int direction - ); - - void setPoints(TTree* tree); - void setPoints(const std::vector& pList){ initPDF(pList); } - template void setPoints(const std::vector& XList, const std::vector& YList, const std::vector& ZList, const std::vector& FcnList){ - std::vector transXList; - std::vector transYList; - std::vector transZList; - std::vector transFcnList; - for (unsigned int ip=0; ip pList = getPoints(transXList, transYList, transZList, transFcnList); - setPoints(pList); - } - -}; - -template void RooNCSplineFactory_3D::setPoints(const std::vector& XList, const std::vector& YList, const std::vector& ZList, const std::vector& FcnList); -template void RooNCSplineFactory_3D::setPoints(const std::vector& XList, const std::vector& YList, const std::vector& ZList, const std::vector& FcnList); - - -#endif - - - diff --git a/interface/RooNCSpline_1D_fast.h b/interface/RooNCSpline_1D_fast.h deleted file mode 100644 index 46cbbb1dd6b..00000000000 --- a/interface/RooNCSpline_1D_fast.h +++ /dev/null @@ -1,68 +0,0 @@ -#ifndef ROONCSPLINE_1D_FAST -#define ROONCSPLINE_1D_FAST - -#include -#include "RooAbsPdf.h" -#include "RooRealProxy.h" -#include "RooRealVar.h" -#include "RooAbsReal.h" -#include "RooNCSplineCore.h" - - -class RooNCSpline_1D_fast : public RooNCSplineCore{ -protected: - BoundaryCondition const bcBeginX; - BoundaryCondition const bcEndX; - - std::vector FcnList; // List of function values - - std::vector kappaX; - std::vector> coefficients; - -public: - RooNCSpline_1D_fast(); - RooNCSpline_1D_fast( - const char* name, - const char* title - ); - RooNCSpline_1D_fast( - const char* name, - const char* title, - RooAbsReal& inXVar, - const std::vector& inXList, - const std::vector& inFcnList, - RooNCSplineCore::BoundaryCondition const bcBeginX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndX_=RooNCSplineCore::bcNaturalSpline, - Bool_t inUseFloor=true, - T inFloorEval=0, - T inFloorInt=0 - ); - RooNCSpline_1D_fast(const RooNCSpline_1D_fast& other, const char* name=0); - virtual TObject* clone(const char* newname)const { return new RooNCSpline_1D_fast(*this, newname); } - inline virtual ~RooNCSpline_1D_fast(){} - - void setRangeValidity(const T valmin, const T valmax, const Int_t whichDirection=0); - - virtual Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0)const; - virtual Double_t analyticalIntegral(Int_t code, const char* rangeName=0)const; - -protected: - virtual void emptyFcnList(){ std::vector tmp; FcnList.swap(tmp); } - - virtual Int_t getWhichBin(const T& val, const Int_t whichDirection)const; - virtual T getTVar(const std::vector& kappas, const T& val, const Int_t& bin, const Int_t whichDirection)const; - virtual void getKappas(std::vector& kappas, const Int_t whichDirection); - - Bool_t testRangeValidity(const T& val, const Int_t whichDirection=0)const; - void cropValueForRange(T& val, const Int_t whichDirection=0)const; - - virtual T interpolateFcn(Int_t code, const char* rangeName=0)const; - - virtual Double_t evaluate()const; - - - ClassDef(RooNCSpline_1D_fast, 2) - -}; - -#endif diff --git a/interface/RooNCSpline_2D_fast.h b/interface/RooNCSpline_2D_fast.h deleted file mode 100644 index 8f8ae2b4a16..00000000000 --- a/interface/RooNCSpline_2D_fast.h +++ /dev/null @@ -1,84 +0,0 @@ -#ifndef ROONCSPLINE_2D_FAST -#define ROONCSPLINE_2D_FAST - -#include -#include "RooAbsPdf.h" -#include "RooRealProxy.h" -#include "RooRealVar.h" -#include "RooAbsReal.h" -#include "RooNCSplineCore.h" - -class RooNCSpline_2D_fast : public RooNCSplineCore{ -protected: - T rangeYmin; - T rangeYmax; - - BoundaryCondition const bcBeginX; - BoundaryCondition const bcEndX; - BoundaryCondition const bcBeginY; - BoundaryCondition const bcEndY; - - RooRealProxy theYVar; - std::vector YList; - - std::vector> FcnList; - - std::vector kappaX; - std::vector kappaY; - std::vector>>> coefficients; // [ix][A_x,B_x,C_x,D_x][iy][A_x_y,B_x_y,C_x_y,D_x_y] - -public: - RooNCSpline_2D_fast(); - RooNCSpline_2D_fast( - const char* name, - const char* title - ); - RooNCSpline_2D_fast( - const char* name, - const char* title, - RooAbsReal& inXVar, - RooAbsReal& inYVar, - const std::vector& inXList, - const std::vector& inYList, - const std::vector>& inFcnList, - RooNCSplineCore::BoundaryCondition const bcBeginX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcBeginY_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndY_=RooNCSplineCore::bcNaturalSpline, - Bool_t inUseFloor=true, - T inFloorEval=0, - T inFloorInt=0 - ); - RooNCSpline_2D_fast(const RooNCSpline_2D_fast& other, const char* name=0); - virtual TObject* clone(const char* newname)const { return new RooNCSpline_2D_fast(*this, newname); } - inline virtual ~RooNCSpline_2D_fast(){} - - void setRangeValidity(const T valmin, const T valmax, const Int_t whichDirection); - - virtual Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const; - virtual Double_t analyticalIntegral(Int_t code, const char* rangeName=0)const; - -protected: - virtual void emptyFcnList(){ std::vector> tmp; FcnList.swap(tmp); } - - unsigned int npointsY()const{ return YList.size(); } - - virtual Int_t getWhichBin(const T& val, const Int_t whichDirection)const; - virtual T getTVar(const std::vector& kappas, const T& val, const Int_t& bin, const Int_t whichDirection)const; - virtual void getKappas(std::vector& kappas, const Int_t whichDirection); - - Bool_t testRangeValidity(const T& val, const Int_t whichDirection)const; - void cropValueForRange(T& val, const Int_t whichDirection)const; - - virtual std::vector> getCoefficientsPerY(const std::vector& kappaX, const TMatrix_t& xAinv, const Int_t& ybin, RooNCSplineCore::BoundaryCondition const& bcBegin, RooNCSplineCore::BoundaryCondition const& bcEnd, const Int_t xbin)const; // xbin can be -1, which means push all of them - - virtual T interpolateFcn(Int_t code, const char* rangeName=0)const; - - virtual Double_t evaluate()const; - - - ClassDef(RooNCSpline_2D_fast, 2) - -}; - -#endif diff --git a/interface/RooNCSpline_3D_fast.h b/interface/RooNCSpline_3D_fast.h deleted file mode 100644 index d32b8dcbf64..00000000000 --- a/interface/RooNCSpline_3D_fast.h +++ /dev/null @@ -1,105 +0,0 @@ -#ifndef ROONCSPLINE_3D_FAST -#define ROONCSPLINE_3D_FAST - -#include -#include "RooAbsPdf.h" -#include "RooRealProxy.h" -#include "RooRealVar.h" -#include "RooAbsReal.h" -#include "RooNCSplineCore.h" - -class RooNCSpline_3D_fast : public RooNCSplineCore{ -protected: - T rangeYmin; - T rangeYmax; - T rangeZmin; - T rangeZmax; - - BoundaryCondition const bcBeginX; - BoundaryCondition const bcEndX; - BoundaryCondition const bcBeginY; - BoundaryCondition const bcEndY; - BoundaryCondition const bcBeginZ; - BoundaryCondition const bcEndZ; - - RooRealProxy theYVar; - RooRealProxy theZVar; - std::vector YList; - std::vector ZList; - - std::vector>> FcnList; - - std::vector kappaX; - std::vector kappaY; - std::vector kappaZ; - std::vector> - >> - >> coefficients; // [ix][A_x,B_x,C_x,D_x][iy][A_x_y,B_x_y,C_x_y,D_x_y][iz][A_x_y_z,B_x_y_z,C_x_y_z,D_x_y_z] - -public: - RooNCSpline_3D_fast(); - RooNCSpline_3D_fast( - const char* name, - const char* title - ); - RooNCSpline_3D_fast( - const char* name, - const char* title, - RooAbsReal& inXVar, - RooAbsReal& inYVar, - RooAbsReal& inZVar, - const std::vector& inXList, - const std::vector& inYList, - const std::vector& inZList, - const std::vector>>& inFcnList, - RooNCSplineCore::BoundaryCondition const bcBeginX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndX_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcBeginY_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndY_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcBeginZ_=RooNCSplineCore::bcNaturalSpline, - RooNCSplineCore::BoundaryCondition const bcEndZ_=RooNCSplineCore::bcNaturalSpline, - Bool_t inUseFloor=true, - T inFloorEval=0, - T inFloorInt=0 - ); - RooNCSpline_3D_fast(const RooNCSpline_3D_fast& other, const char* name=0); - virtual TObject* clone(const char* newname)const { return new RooNCSpline_3D_fast(*this, newname); } - inline virtual ~RooNCSpline_3D_fast(){} - - void setRangeValidity(const T valmin, const T valmax, const Int_t whichDirection); - - virtual Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const; - virtual Double_t analyticalIntegral(Int_t code, const char* rangeName=0)const; - -protected: - virtual void emptyFcnList(){ std::vector>> tmp; FcnList.swap(tmp); } - - unsigned int npointsY()const{ return YList.size(); } - unsigned int npointsZ()const{ return ZList.size(); } - - virtual Int_t getWhichBin(const T& val, const Int_t whichDirection)const; - virtual T getTVar(const std::vector& kappas, const T& val, const Int_t& bin, const Int_t whichDirection)const; - virtual void getKappas(std::vector& kappas, const Int_t whichDirection); - - Bool_t testRangeValidity(const T& val, const Int_t whichDirection)const; - void cropValueForRange(T& val, const Int_t whichDirection)const; - - virtual std::vector> getCoefficientsPerYPerZ( - const std::vector& kappaX, const TMatrix_t& xAinv, - const Int_t& ybin, const Int_t& zbin, - RooNCSplineCore::BoundaryCondition const& bcBegin, RooNCSplineCore::BoundaryCondition const& bcEnd, - const Int_t xbin - )const; // xbin can be -1, which means push all of them - - virtual T interpolateFcn(Int_t code, const char* rangeName=0)const; - - virtual Double_t evaluate()const; - - - ClassDef(RooNCSpline_3D_fast, 2) - -}; - -#endif diff --git a/interface/RooPiecewisePolynomial.h b/interface/RooPiecewisePolynomial.h deleted file mode 100644 index a548af454cb..00000000000 --- a/interface/RooPiecewisePolynomial.h +++ /dev/null @@ -1,50 +0,0 @@ -#ifndef ROOPIECEWISEPOLYNOMIAL_H -#define ROOPIECEWISEPOLYNOMIAL_H -#include -#include "RooAbsReal.h" -#include "RooArgList.h" -#include "RooRealProxy.h" -#include "RooListProxy.h" - -// Piecewise polynomial that looks like -// -- . -- ... -- . -- -class RooPiecewisePolynomial : public RooAbsReal{ -protected: - RooRealProxy xvar; - - // First [0,...,nnodes-1] parameters are nodes - // There should be 2*ndof_endfcn+(nfcn-2)*ndof_middlefcn more parameters for the free dofs in the functions - // where ndof_endfcn=polyndof-1 and ndof_middlefcn=polyndof-2 - // For nfcn=4, polyndof=4, the order of parameters is - // node 0, node 1, node 2 (nnodes=nfcn-1), and then - // a0 + a1*x + a2*x^2 + not_a_par*x^3 (= fcn 0) - // b0 + b1*x + not_a_par*x^2 + not_a_par*x^3 (= fcn 1) - // c0 + c1*x + not_a_par*x^2 + not_a_par*x^3 (= fcn 2) - // d0 + d1*x + not_a_par*x^2 + d3*x^3 (= fcn 3) - // Continuity between the fcns is guaranteed, but notice that smoothness is not. - RooListProxy parList; - - const int nfcn; // Number of piecewise functions - const int polyndof; // Ndof of the polynomial (e.g. 4: cubic) - const int nnodes; // How many nodes are in between (= nfcn-1) - const int ndof_endfcn; // Number of degrees of freedom in fcns in the middle of the nodes (= polyndof-1) - const int ndof_middlefcn; // Number of degrees of freedom in fcns outside the nodes (= polyndof-2) - - double eval(double x, std::vector const& par)const; - -public: - RooPiecewisePolynomial(const int nfcn_=1, const int polyndof_=1); - RooPiecewisePolynomial(const char* name, const char* title, const int nfcn_=1, const int polyndof_=1); - RooPiecewisePolynomial(const char* name, const char* title, RooAbsReal& xvar_, RooArgList const& parList_, const int nfcn_, const int polyndof_); - RooPiecewisePolynomial(RooPiecewisePolynomial const& other, const char* name=0); - virtual TObject* clone(const char* newname)const{ return new RooPiecewisePolynomial(*this, newname); } - inline virtual ~RooPiecewisePolynomial(){} - - double evaluate()const; - double evaluate(double* x, double* p)const; // For calling in a TF1 object - - ClassDef(RooPiecewisePolynomial, 1) - -}; - -#endif diff --git a/interface/RooRealFlooredSumPdf.h b/interface/RooRealFlooredSumPdf.h deleted file mode 100755 index 660fa784c01..00000000000 --- a/interface/RooRealFlooredSumPdf.h +++ /dev/null @@ -1,95 +0,0 @@ -/***************************************************************************** - * Project: RooFit * - * * - * This code was autogenerated by RooClassFactory * - *****************************************************************************/ - -#ifndef ROOREALFLOOREDSUMPDF -#define ROOREALFLOOREDSUMPDF - - -#include -#include "RooAbsReal.h" -#include "RooAbsPdf.h" -#include "RooRealVar.h" -#include "RooRealProxy.h" -#include "RooAbsCategory.h" -#include "RooCategoryProxy.h" -#include "RooListProxy.h" -#include "RooAICRegistry.h" -#include "RooObjCacheManager.h" -#include "TH3F.h" -#include "TH1.h" -#include "RooDataHist.h" -#include "RooHistFunc.h" - - -class RooRealFlooredSumPdf : public RooAbsPdf { -public: - - RooRealFlooredSumPdf(); - RooRealFlooredSumPdf(const char *name, const char *title); - RooRealFlooredSumPdf(const char *name, const char *title, const RooArgList& funcList, const RooArgList& coefList, Bool_t extended = kFALSE); - RooRealFlooredSumPdf(const RooRealFlooredSumPdf& other, const char* name = 0); - virtual TObject* clone(const char* newname) const { return new RooRealFlooredSumPdf(*this, newname); } - virtual ~RooRealFlooredSumPdf(); - - Double_t evaluate() const; - virtual Bool_t checkObservables(const RooArgSet* nset) const; - - virtual Bool_t forceAnalyticalInt(const RooAbsArg&) const { return kTRUE; } - Int_t getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& numVars, const RooArgSet* normSet, const char* rangeName = 0) const; - Double_t analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName = 0) const; - - const RooArgList& funcList() const { return _funcList; } - const RooArgList& coefList() const { return _coefList; } - - void setFloor(Double_t val); - - virtual ExtendMode extendMode() const; - - virtual Double_t expectedEvents(const RooArgSet* nset) const; - virtual Double_t expectedEvents(const RooArgSet& nset) const { - // Return expected number of events for extended likelihood calculation - // which is the sum of all coefficients - return expectedEvents(&nset); - } - - void printMetaArgs(std::ostream& os) const; - - - virtual std::list* binBoundaries(RooAbsRealLValue& /*obs*/, Double_t /*xlo*/, Double_t /*xhi*/) const; - virtual std::list* plotSamplingHint(RooAbsRealLValue& /*obs*/, Double_t /*xlo*/, Double_t /*xhi*/) const; - Bool_t isBinnedDistribution(const RooArgSet& obs) const; - - -protected: - - class CacheElem : public RooAbsCacheElement { - public: - CacheElem() {}; - virtual ~CacheElem() {}; - virtual RooArgList containedArgs(Action) { RooArgList ret(_funcIntList); ret.add(_funcNormList); return ret; } - RooArgList _funcIntList; - RooArgList _funcNormList; - }; - mutable RooObjCacheManager _normIntMgr; // The integration cache manager - - - Bool_t _haveLastCoef; - - RooListProxy _funcList; // List of component FUNCs - RooListProxy _coefList; // List of coefficients - TIterator* _funcIter; //! Iterator over FUNC list - TIterator* _coefIter; //! Iterator over coefficient list - Bool_t _extended; // Allow use as extended p.d.f. - Bool_t _doFloor; - Double_t _floorVal; - -private: - - ClassDef(RooRealFlooredSumPdf, 2) // PDF constructed from a sum of (non-pdf) functions -}; - - -#endif diff --git a/src/AsymQuad.cc b/src/AsymQuad.cc deleted file mode 100644 index e9886929c44..00000000000 --- a/src/AsymQuad.cc +++ /dev/null @@ -1,128 +0,0 @@ -#include "../interface/AsymQuad.h" - -#include -#include -#include - -AsymQuad::AsymQuad() : -RooAbsReal(), -_funcList("funcList", "List of functions", this), -_coefList("coefList", "List of coefficients", this), -smoothRegion_(0), -smoothAlgo_(0) -{ - _funcIter = _funcList.createIterator(); - _coefIter = _coefList.createIterator(); -} - -AsymQuad::AsymQuad(const char *name, const char *title, const RooArgList& inFuncList, const RooArgList& inCoefList, Double_t smoothRegion, Int_t smoothAlgo) : -RooAbsReal(name, title), -_funcList("funcList", "List of functions", this), -_coefList("coefList", "List of coefficients", this), -smoothRegion_(smoothRegion), -smoothAlgo_(smoothAlgo) -{ - if (inFuncList.getSize()!=2*inCoefList.getSize()+1) { - coutE(InputArguments) << "AsymQuad::AsymQuad(" << GetName() - << ") number of functions and coefficients inconsistent, must have Nfunc=1+2*Ncoef" << std::endl; - assert(0); - } - - TIterator* funcIter = inFuncList.createIterator(); - RooAbsArg* func; - while ((func = (RooAbsArg*)funcIter->Next())) { - if (!dynamic_cast(func)) { - coutE(InputArguments) << "ERROR: AsymQuad::AsymQuad(" << GetName() << ") function " << func->GetName() << " is not of type RooAbsReal" << std::endl; - assert(0); - } - _funcList.add(*func); - } - delete funcIter; - - TIterator* coefIter = inCoefList.createIterator(); - RooAbsArg* coef; - while ((coef = (RooAbsArg*)coefIter->Next())) { - if (!dynamic_cast(coef)) { - coutE(InputArguments) << "ERROR: AsymQuad::AsymQuad(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal" << std::endl; - assert(0); - } - _coefList.add(*coef); - } - delete coefIter; - - _funcIter = _funcList.createIterator(); - _coefIter = _coefList.createIterator(); -} - -AsymQuad::AsymQuad(const AsymQuad& other, const char* name): -RooAbsReal(other, name), -_funcList("!funcList", this, other._funcList), -_coefList("!coefList", this, other._coefList), -smoothRegion_(other.smoothRegion_), -smoothAlgo_(other.smoothAlgo_) -{ - _funcIter = _funcList.createIterator(); - _coefIter = _coefList.createIterator(); -} - -AsymQuad::~AsymQuad() { - delete _funcIter; - delete _coefIter; -} - -Double_t AsymQuad::evaluate() const { - Double_t result(0); - - _funcIter->Reset(); - _coefIter->Reset(); - RooAbsReal* coef; - RooAbsReal* func = (RooAbsReal*)_funcIter->Next(); - - Double_t central = func->getVal(); - result = central; - - while ((coef=(RooAbsReal*)_coefIter->Next())) { - Double_t coefVal = coef->getVal(); - RooAbsReal* funcUp = (RooAbsReal*)_funcIter->Next(); - RooAbsReal* funcDn = (RooAbsReal*)_funcIter->Next(); - result += interpolate(coefVal, central, funcUp->getVal(), funcDn->getVal()); - } - - return result; -} - -Double_t AsymQuad::interpolate(Double_t theta_, Double_t valueCenter_, Double_t valueHigh_, Double_t valueLow_) const { - if (smoothAlgo_<0) return 0; - else{ - if (fabs(theta_)>=smoothRegion_) return theta_ * (theta_ > 0 ? valueHigh_ - valueCenter_ : valueCenter_ - valueLow_); - if (smoothAlgo_ == 0) { - // Quadratic interpolation null at zero and continuous at boundaries but not smooth at boundaries - Double_t c_up = +theta_ * (smoothRegion_ + theta_) / (2 * smoothRegion_); - Double_t c_dn = -theta_ * (smoothRegion_ - theta_) / (2 * smoothRegion_); - Double_t c_cen = -theta_ * theta_ / smoothRegion_; - return c_up * valueHigh_ + c_dn * valueLow_ + c_cen * valueCenter_; - } - else if (smoothAlgo_ == 1){ - // Quadratic interpolation that is everywhere differentiable but not null at zero - Double_t c_up = (smoothRegion_ + theta_) * (smoothRegion_ + theta_) / (4 * smoothRegion_); - Double_t c_dn = (smoothRegion_ - theta_) * (smoothRegion_ - theta_) / (4 * smoothRegion_); - Double_t c_cen = -c_up - c_dn; - return c_up * valueHigh_ + c_dn * valueLow_ + c_cen * valueCenter_; - } - else/* if (smoothAlgo_ == 2)*/{ - // Quadratic interpolation that is everywhere differentiable and null at zero - Double_t cnorm = theta_/smoothRegion_; - Double_t cnorm2 = pow(cnorm, 2); - Double_t hi = valueHigh_ - valueCenter_; - Double_t lo = valueLow_ - valueCenter_; - Double_t sum = hi+lo; - Double_t diff = hi-lo; - Double_t a = theta_/2.; // cnorm*smoothRegion_ - Double_t b = 0.125 * cnorm * (cnorm2 * (3.*cnorm2 - 10.) + 15.); - Double_t result = a*(diff + b*sum); - return result; - } - } -} - -ClassImp(AsymQuad) diff --git a/src/RooFuncPdf.cc b/src/RooFuncPdf.cc deleted file mode 100644 index 1bfac2e790a..00000000000 --- a/src/RooFuncPdf.cc +++ /dev/null @@ -1,3 +0,0 @@ -#include "../interface/RooFuncPdf.h" - -ClassImp(RooFuncPdf) diff --git a/src/RooNCSplineCore.cc b/src/RooNCSplineCore.cc deleted file mode 100644 index 68df9d9e025..00000000000 --- a/src/RooNCSplineCore.cc +++ /dev/null @@ -1,284 +0,0 @@ -#include "../interface/RooNCSplineCore.h" -#include -#include "TMath.h" -#include "TIterator.h" -#include "Riostream.h" - -using namespace TMath; -using namespace RooFit; -using namespace std; - - -ClassImp(RooNCSplineCore) - -RooNCSplineCore::RooNCSplineCore() : -RooAbsReal(), -verbosity(RooNCSplineCore::kSilent), -useFloor(true), floorEval(0), floorInt(0), -rangeXmin(1), rangeXmax(-1), -theXVar("theXVar", "theXVar", this), -leafDepsList("leafDepsList", "leafDepsList", this) -{} - -RooNCSplineCore::RooNCSplineCore( - const char* name, - const char* title - ) : - RooAbsReal(name, title), - verbosity(RooNCSplineCore::kSilent), - useFloor(true), floorEval(0), floorInt(0), - rangeXmin(1), rangeXmax(-1), - theXVar("theXVar", "theXVar", this), - leafDepsList("leafDepsList", "leafDepsList", this) -{} - -RooNCSplineCore::RooNCSplineCore( - const char* name, - const char* title, - RooAbsReal& inXVar, - const std::vector& inXList, - Bool_t inUseFloor, - T inFloorEval, - T inFloorInt - ) : - RooAbsReal(name, title), - verbosity(RooNCSplineCore::kSilent), - useFloor(inUseFloor), floorEval(inFloorEval), floorInt(inFloorInt), - rangeXmin(1), rangeXmax(-1), - theXVar("theXVar", "theXVar", this, inXVar), - leafDepsList("leafDepsList", "leafDepsList", this), - XList(inXList) -{} - -RooNCSplineCore::RooNCSplineCore( - const RooNCSplineCore& other, - const char* name - ) : - RooAbsReal(other, name), - verbosity(other.verbosity), - useFloor(other.useFloor), floorEval(other.floorEval), floorInt(other.floorInt), - rangeXmin(other.rangeXmin), rangeXmax(other.rangeXmax), - theXVar("theXVar", this, other.theXVar), - leafDepsList("leafDepsList", this, other.leafDepsList), - XList(other.XList) -{} - -void RooNCSplineCore::setVerbosity(VerbosityLevel flag){ verbosity = flag; } -void RooNCSplineCore::setEvalFloor(RooNCSplineCore::T val){ floorEval = val; } -void RooNCSplineCore::setIntFloor(RooNCSplineCore::T val){ floorInt = val; } -void RooNCSplineCore::doFloor(Bool_t flag){ useFloor = flag; } - -void RooNCSplineCore::getBArray(const std::vector& kappas, const vector& fcnList, std::vector& BArray, BoundaryCondition const& bcBegin, BoundaryCondition const& bcEnd)const{ - BArray.clear(); - int npoints=kappas.size(); - if (npoints!=(int)fcnList.size()){ - coutE(InputArguments) << "RooNCSplineCore::getBArray: Dim(kappas)=" << npoints << " != Dim(fcnList)=" << fcnList.size() << endl; - assert(0); - } - if ( - bcBegin==bcQuadraticWithNullSlope || bcEnd==bcQuadraticWithNullSlope - || - bcBegin==bcQuadratic || bcEnd==bcQuadratic - ){ - int nthr=1+(bcBegin==bcQuadraticWithNullSlope || bcBegin==bcQuadratic ? 1 : 0) + (bcEnd==bcQuadraticWithNullSlope || bcEnd==bcQuadratic ? 1 : 0); - if (npoints<=nthr){ - cerr << "Npoints " << npoints << " <= " << nthr << endl; - assert(0); - } - } - if (npoints>1){ - RooNCSplineCore::T bcval=0, ecval=0; - // First point constraint - switch (bcBegin){ - case bcApproximatedSecondDerivative: - bcval=(npoints<3 ? 0 : ((fcnList.at(2)-fcnList.at(1))*kappas.at(1)-(fcnList.at(1)-fcnList.at(0))*kappas.at(0))/(0.5/kappas.at(0)+0.5/kappas.at(1))); - case bcNaturalSpline: - BArray.push_back(3.*(fcnList.at(1)-fcnList.at(0)) - bcval/2./pow(kappas.at(0), 2)); - break; - case bcApproximatedSlope: - bcval=(fcnList.at(1)-fcnList.at(0))*kappas.at(0); - case bcClamped: - BArray.push_back(bcval/kappas.at(0)); - break; - case bcQuadratic: - bcval=2.*(fcnList.at(1)-fcnList.at(0)); - BArray.push_back(bcval); - break; - case bcQuadraticWithNullSlope: - bcval=2.*(fcnList.at(1)-fcnList.at(0)); - BArray.push_back(0); - BArray.push_back(bcval*kappas.at(0)/kappas.at(1)); - break; - default: - cerr << "RooNCSplineCore::getBArray: bcBegin " << bcBegin << " is not implemented!" << endl; - assert(0); - } - // Intermediate point constraint is always D0, D1 and D2 continuous - for (int j=1; j& kappas, vector>& AArray, BoundaryCondition const& bcBegin, BoundaryCondition const& bcEnd)const{ - AArray.clear(); - Int_t npoints = kappas.size(); - for (int i=0; i Ai(npoints, 0.); - if (npoints==1) Ai[0]=1; - else if (i==0){ - switch(bcBegin){ - case bcNaturalSpline: - case bcApproximatedSecondDerivative: - Ai[0]=2; Ai[1]=kappas.at(1)/kappas.at(0); - break; - case bcClamped: - case bcApproximatedSlope: - case bcQuadraticWithNullSlope: - Ai[0]=1; - break; - case bcQuadratic: - Ai[0]=1; Ai[1]=kappas.at(1)/kappas.at(0); - break; - default: - cerr << "RooNCSplineCore::getAArray: bcBegin " << bcBegin << " is not implemented!" << endl; - assert(0); - } - } - else if (i==npoints-1){ - switch(bcEnd){ - case bcNaturalSpline: - case bcApproximatedSecondDerivative: - Ai[npoints-2]=1; Ai[npoints-1]=2.*kappas.at(npoints-1)/kappas.at(npoints-2); - break; - case bcClamped: - case bcApproximatedSlope: - case bcQuadraticWithNullSlope: - Ai[npoints-1]=1; - break; - case bcQuadratic: - Ai[npoints-2]=1; Ai[npoints-1]=kappas.at(npoints-1)/kappas.at(npoints-2); - break; - default: - cerr << "RooNCSplineCore::getAArray: bcEnd " << bcEnd << " is not implemented!" << endl; - assert(0); - } - } - else{ - if ((i==1 && bcBegin==bcQuadraticWithNullSlope) || (i==npoints-2 && bcEnd==bcQuadraticWithNullSlope)) Ai[i]=1; - else{ - RooNCSplineCore::T kappa_j = kappas.at(i); - RooNCSplineCore::T kappa_jmo = kappas.at(i-1); - RooNCSplineCore::T kappa_jpo = kappas.at(i+1); - - Ai[i-1]=1; - Ai[i]=2.*kappa_j/kappa_jmo*(1.+kappa_j/kappa_jmo); - Ai[i+1]=kappa_j*kappa_jpo/pow(kappa_jmo, 2); - } - } - AArray.push_back(Ai); - } -} - -vector RooNCSplineCore::getCoefficients(const TVector_t& S, const vector& kappas, const vector& fcnList, const Int_t& bin)const{ - DefaultAccumulator A, B, C, D; - vector res; - - const int fcnsize = fcnList.size(); - if (fcnsize>bin){ - A=fcnList.at(bin); - B=S[bin]; - if (fcnsize>(bin+1)){ - DefaultAccumulator dFcn = fcnList.at(bin+1); - dFcn -= A; - - C += RooNCSplineCore::T(3.)*dFcn; - C -= 2.*B; - C -= S[bin+1]*kappas.at(bin+1)/kappas.at(bin); - - D += RooNCSplineCore::T(-2.)*dFcn; - D += B; - D += S[bin+1]*kappas.at(bin+1)/kappas.at(bin); - } - } - - res.push_back(A); - res.push_back(B); - res.push_back(C); - res.push_back(D); - return res; -} -vector> RooNCSplineCore::getCoefficientsAlongDirection(const std::vector& kappas, const TMatrix_t& Ainv, const vector& fcnList, BoundaryCondition const& bcBegin, BoundaryCondition const& bcEnd, const Int_t pickBin)const{ - vector BArray; - getBArray(kappas, fcnList, BArray, bcBegin, bcEnd); - - Int_t npoints = BArray.size(); - TVector_t Btrans(npoints); - for (int i=0; i> coefs; - for (Int_t bin=0; bin<(npoints>1 ? npoints-1 : 1); bin++){ - if (pickBin>=0 && bin!=pickBin) continue; - vector coef = getCoefficients(Strans, kappas, fcnList, bin); - coefs.push_back(coef); - } - return coefs; -} - -RooNCSplineCore::T RooNCSplineCore::evalSplineSegment(const std::vector& coefs, const RooNCSplineCore::T& kappa, const RooNCSplineCore::T& tup, const RooNCSplineCore::T& tdn, Bool_t doIntegrate)const{ - DefaultAccumulator res; - for (unsigned int ic=0; icleafNodeServerList(&deps, 0, true); - set.add(deps); -} -void RooNCSplineCore::addLeafDependents(RooArgSet& set){ - TIterator* iter = set.createIterator(); - RooAbsArg* absarg; - while ((absarg = (RooAbsArg*)iter->Next())){ if (dynamic_cast(absarg)) leafDepsList.add(*absarg); } - delete iter; -} diff --git a/src/RooNCSplineFactory_1D.cc b/src/RooNCSplineFactory_1D.cc deleted file mode 100644 index 2555d2e029c..00000000000 --- a/src/RooNCSplineFactory_1D.cc +++ /dev/null @@ -1,88 +0,0 @@ -#include "../interface/RooNCSplineFactory_1D.h" -#include - -using namespace std; - - -RooNCSplineFactory_1D::RooNCSplineFactory_1D( - RooAbsReal& splineVar_, TString appendName_, - RooNCSplineCore::BoundaryCondition const bcBeginX_, - RooNCSplineCore::BoundaryCondition const bcEndX_ -) : - appendName(appendName_), - bcBeginX(bcBeginX_), bcEndX(bcEndX_), - splineVar(&splineVar_), - fcn(0), - PDF(0) -{} -RooNCSplineFactory_1D::~RooNCSplineFactory_1D(){ - destroyPDF(); -} -void RooNCSplineFactory_1D::setPoints(TTree* tree){ - vector> pList; - RooNCSplineCore::T x, fcn; - tree->SetBranchAddress("X", &x); - tree->SetBranchAddress("Fcn", &fcn); - int n = tree->GetEntries(); - for (int ip=0; ipGetEntry(ip); pList.push_back(pair(x, fcn)); } - setPoints(pList); -} -void RooNCSplineFactory_1D::setPoints(TGraph* tg){ - vector> pList; - double* xx = tg->GetX(); - double* yy = tg->GetY(); - int n = tg->GetN(); - for (int ip=0; ip(xx[ip], yy[ip])); - setPoints(pList); -} -const std::vector> RooNCSplineFactory_1D::getPoints(const std::vector& XList, const std::vector& FcnList){ - const unsigned int nX = XList.size(); - const unsigned int n = FcnList.size(); - if (nX!=n){ - cerr << "RooNCSplineFactory_1D::getPoints: nX=" << nX << " != nFcn=" << n << endl; - assert(0); - } - std::vector> pList; pList.reserve(n); - for (unsigned int ip=0; ip(XList.at(ip), FcnList.at(ip))); - return pList; -} - -void RooNCSplineFactory_1D::destroyPDF(){ delete PDF; PDF=0; delete fcn; fcn=0; } -void RooNCSplineFactory_1D::initPDF(const std::vector>& pList){ - destroyPDF(); - - const unsigned int n = pList.size(); - std::vector XList; - std::vector FcnList; - for (unsigned int ip=0; ip& list, RooNCSplineCore::T val){ - for (unsigned int ip=0; ip RooNCSplineFactory_2D::getPoints( - const std::vector& XList, - const std::vector& YList, - const std::vector& FcnList - ){ - const unsigned int nX = XList.size(); - const unsigned int nY = YList.size(); - const unsigned int n = FcnList.size(); - if (nX*nY!=n){ - cerr << "RooNCSplineFactory_2D::getPoints: nX=" << nX << " x nY=" << nY << " != nFcn=" << n << endl; - assert(0); - } - - std::vector pList; pList.reserve(n); - for (unsigned int ix=0; ix& pList){ - destroyPDF(); - - const unsigned int n = pList.size(); - vector XList; - vector YList; - vector> FcnList; - for (unsigned int ip=0; ip dum; dum.reserve(XList.size()); - for (unsigned int ix=0; ix pList; - RooNCSplineCore::T x, y, fcn; - tree->SetBranchAddress("X", &x); - tree->SetBranchAddress("Y", &y); - tree->SetBranchAddress("Fcn", &fcn); - int n = tree->GetEntries(); - for (int ip=0; ipGetEntry(ip); pList.push_back(splineTriplet_t(x, y, fcn)); } - setPoints(pList); -} - -void RooNCSplineFactory_2D::setEndConditions( - RooNCSplineCore::BoundaryCondition const bcBegin, - RooNCSplineCore::BoundaryCondition const bcEnd, - const unsigned int direction -){ - switch (direction){ - case 0: - bcBeginX=bcBegin; - bcEndX=bcEnd; - break; - case 1: - bcBeginY=bcBegin; - bcEndY=bcEnd; - break; - default: - // Do nothing - break; - } -} diff --git a/src/RooNCSplineFactory_3D.cc b/src/RooNCSplineFactory_3D.cc deleted file mode 100644 index 1b51a056bf5..00000000000 --- a/src/RooNCSplineFactory_3D.cc +++ /dev/null @@ -1,146 +0,0 @@ -#include "../interface/RooNCSplineFactory_3D.h" - -using namespace std; - - -RooNCSplineFactory_3D::RooNCSplineFactory_3D( - RooAbsReal& XVar_, RooAbsReal& YVar_, RooAbsReal& ZVar_, TString appendName_, - RooNCSplineCore::BoundaryCondition const bcBeginX_, - RooNCSplineCore::BoundaryCondition const bcEndX_, - RooNCSplineCore::BoundaryCondition const bcBeginY_, - RooNCSplineCore::BoundaryCondition const bcEndY_, - RooNCSplineCore::BoundaryCondition const bcBeginZ_, - RooNCSplineCore::BoundaryCondition const bcEndZ_ -) : - appendName(appendName_), - bcBeginX(bcBeginX_), bcEndX(bcEndX_), - bcBeginY(bcBeginY_), bcEndY(bcEndY_), - bcBeginZ(bcBeginZ_), bcEndZ(bcEndZ_), - XVar(&XVar_), YVar(&YVar_), ZVar(&ZVar_), - fcn(0), - PDF(0) -{} -RooNCSplineFactory_3D::~RooNCSplineFactory_3D(){ - destroyPDF(); -} - -void RooNCSplineFactory_3D::addUnique(std::vector& list, RooNCSplineCore::T val){ - for (unsigned int ip=0; ip RooNCSplineFactory_3D::getPoints( - const std::vector& XList, - const std::vector& YList, - const std::vector& ZList, - const std::vector& FcnList - ){ - const unsigned int nX = XList.size(); - const unsigned int nY = YList.size(); - const unsigned int nZ = ZList.size(); - const unsigned int n = FcnList.size(); - if (nX*nY*nZ!=n){ - cerr << "RooNCSplineFactory_3D::getPoints: nX=" << nX << " x nY=" << nY << " x nZ=" << nZ << " != nFcn=" << n << endl; - assert(0); - } - - std::vector pList; pList.reserve(n); - for (unsigned int ix=0; ix& pList){ - destroyPDF(); - - const unsigned int n = pList.size(); - vector XList; - vector YList; - vector ZList; - vector>> FcnList; - for (unsigned int ip=0; ip> dumz; - dumz.reserve(YList.size()); - for (unsigned int iy=0; iy dumy; - dumy.reserve(XList.size()); - for (unsigned int ix=0; ix pList; - RooNCSplineCore::T x, y, z, fcn; - tree->SetBranchAddress("X", &x); - tree->SetBranchAddress("Y", &y); - tree->SetBranchAddress("Z", &z); - tree->SetBranchAddress("Fcn", &fcn); - int n = tree->GetEntries(); - for (int ip=0; ipGetEntry(ip); pList.push_back(splineQuadruplet_t(x, y, z, fcn)); } - setPoints(pList); -} - -void RooNCSplineFactory_3D::setEndConditions( - RooNCSplineCore::BoundaryCondition const bcBegin, - RooNCSplineCore::BoundaryCondition const bcEnd, - const unsigned int direction -){ - switch (direction){ - case 0: - bcBeginX=bcBegin; - bcEndX=bcEnd; - break; - case 1: - bcBeginY=bcBegin; - bcEndY=bcEnd; - break; - case 2: - bcBeginZ=bcBegin; - bcEndZ=bcEnd; - break; - default: - // Do nothing - break; - } -} diff --git a/src/RooNCSpline_1D_fast.cc b/src/RooNCSpline_1D_fast.cc deleted file mode 100644 index fd0de69a443..00000000000 --- a/src/RooNCSpline_1D_fast.cc +++ /dev/null @@ -1,226 +0,0 @@ -#include "../interface/RooNCSpline_1D_fast.h" -#include -#include "TMath.h" -#include "Riostream.h" -#include "RooAbsReal.h" - -using namespace TMath; -using namespace RooFit; -using namespace std; - - -ClassImp(RooNCSpline_1D_fast) - -RooNCSpline_1D_fast::RooNCSpline_1D_fast() : - RooNCSplineCore(), - bcBeginX(RooNCSplineCore::bcNaturalSpline), bcEndX(RooNCSplineCore::bcNaturalSpline) -{} - -RooNCSpline_1D_fast::RooNCSpline_1D_fast( - const char* name, - const char* title - ) : - RooNCSplineCore(name, title), - bcBeginX(RooNCSplineCore::bcNaturalSpline), bcEndX(RooNCSplineCore::bcNaturalSpline) -{} - -RooNCSpline_1D_fast::RooNCSpline_1D_fast( - const char* name, - const char* title, - RooAbsReal& inXVar, - const std::vector& inXList, - const std::vector& inFcnList, - RooNCSplineCore::BoundaryCondition const bcBeginX_, - RooNCSplineCore::BoundaryCondition const bcEndX_, - Bool_t inUseFloor, - T inFloorEval, - T inFloorInt - ) : - RooNCSplineCore(name, title, inXVar, inXList, inUseFloor, inFloorEval, inFloorInt), - bcBeginX(bcBeginX_), bcEndX(bcEndX_), - FcnList(inFcnList) -{ - if (npointsX()>1){ - int npoints; - - vector> xA; getKappas(kappaX, 0); getAArray(kappaX, xA, bcBeginX, bcEndX); - npoints=kappaX.size(); - TMatrix_t xAtrans(npoints, npoints); - for (int i=0; i res; - - if (verbosity==RooNCSplineCore::kVerbose) cout << "RooNCSpline_1D_fast(" << GetName() << ")::interpolateFcn begin with code: " << code << endl; - - // Get bins - Int_t xbin=-1, xbinmin=-1, xbinmax=-1; - RooNCSplineCore::T tx=0, txmin=0, txmax=0; - if (code==0 || code%2!=0){ // Case to just compute the value at x - if (!testRangeValidity(theXVar)) return 0; - xbin = getWhichBin(theXVar, 0); - tx = getTVar(kappaX, theXVar, xbin, 0); - } - else{ // Case to integrate along x - RooNCSplineCore::T coordmin = theXVar.min(rangeName); cropValueForRange(coordmin); - RooNCSplineCore::T coordmax = theXVar.max(rangeName); cropValueForRange(coordmax); - xbinmin = getWhichBin(coordmin, 0); - txmin = getTVar(kappaX, coordmin, xbinmin, 0); - xbinmax = getWhichBin(coordmax, 0); - txmax = getTVar(kappaX, coordmax, xbinmax, 0); - } - - int nxbins = (int)coefficients.size(); - for (int ix=0; ix=0 && ix!=xbin) - || - (xbinmin>=0 && xbinmax>=xbinmin && !(xbinmin<=ix && ix<=xbinmax)) - ) continue; - - RooNCSplineCore::T txlow=0, txhigh=1; - if (code>0 && code%2==0){ - if (ix==xbinmin) txlow=txmin; - if (ix==xbinmax) txhigh=txmax; - } - else txhigh=tx; - - // Get the x coefficients at bin ix and evaluate value of spline at x - res += evalSplineSegment(coefficients.at(ix), kappaX.at(ix), txhigh, txlow, (code>0 && code%2==0)); - } - - return res; -} - -void RooNCSpline_1D_fast::getKappas(vector& kappas, const Int_t /*whichDirection*/){ - kappas.clear(); - RooNCSplineCore::T kappa=1; - - Int_t npoints; - vector const* coord; - npoints=npointsX(); - coord=&XList; - - for (Int_t j=0; jat(j); - RooNCSplineCore::T val_jpo = coord->at(j+1); - RooNCSplineCore::T val_diff = (val_jpo-val_j); - if (fabs(val_diff)>RooNCSplineCore::T(0)) kappa = 1./val_diff; - else kappa = 0; - kappas.push_back(kappa); - } - kappas.push_back(kappa); // Push the same kappa_(N-1)=kappa_(N-2) at the end point -} -Int_t RooNCSpline_1D_fast::getWhichBin(const RooNCSplineCore::T& val, const Int_t /*whichDirection*/)const{ - Int_t bin=-1; - RooNCSplineCore::T valj, valjpo; - Int_t npoints; - vector const* coord; - coord=&XList; - npoints=npointsX(); - - if (npoints<=1) bin=0; - else{ - valjpo = coord->at(0); - for (Int_t j=0; jat(j); - valjpo = coord->at(j+1); - if (val=valj){ bin=j; break; } - } - if (bin==-1 && val>=valjpo) bin=npoints-2; - else if (bin==-1) bin=0; - } - - return bin; -} -RooNCSplineCore::T RooNCSpline_1D_fast::getTVar(const vector& kappas, const RooNCSplineCore::T& val, const Int_t& bin, const Int_t /*whichDirection*/)const{ - const RooNCSplineCore::T& K=kappas.at(bin); - return (val-XList.at(bin))*K; -} - -Double_t RooNCSpline_1D_fast::evaluate() const{ - Double_t value = interpolateFcn(0); - if (useFloor && value=RooNCSplineCore::kError) coutE(Eval) << "RooNCSpline_1D_fast ERROR::RooNCSpline_1D_fast(" << GetName() << ") evaluation returned " << value << " at x = " << theXVar << endl; - value = floorEval; - } - if (verbosity==RooNCSplineCore::kVerbose){ - cout << "RooNCSpline_1D_fast(" << GetName() << ")::evaluate = " << value << " at x = " << theXVar << endl; - RooArgSet Xdeps; theXVar.absArg()->leafNodeServerList(&Xdeps, 0, true); - TIterator* iter = Xdeps.createIterator(); - RooAbsArg* var; - while ((var = (RooAbsArg*)iter->Next())){ - cout << var->GetName() << " value = " << dynamic_cast(var)->getVal() << endl; - } - delete iter; - cout << endl; - } - return value; -} -Int_t RooNCSpline_1D_fast::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const{ - if (_forceNumInt) return 0; - Int_t code=0; - if (dynamic_cast(theXVar.absArg())!=0){ - if (matchArgs(allVars, analVars, theXVar)) code=2; - } - return code; -} -Double_t RooNCSpline_1D_fast::analyticalIntegral(Int_t code, const char* rangeName) const{ - Double_t value = interpolateFcn(code, rangeName); - if (useFloor && value=RooNCSplineCore::kError) coutE(Integration) << "RooNCSpline_1D_fast ERROR::RooNCSpline_1D_fast(" << GetName() << ") integration returned " << value << " for code = " << code << endl; - value = floorInt; - } - if (verbosity==RooNCSplineCore::kVerbose){ cout << "RooNCSpline_1D_fast(" << GetName() << ")::analyticalIntegral = " << value << " for code = " << code << endl; } - return value; -} - -Bool_t RooNCSpline_1D_fast::testRangeValidity(const T& val, const Int_t /*whichDirection*/) const{ - const T* range[2]; - range[0] = &rangeXmin; - range[1] = &rangeXmax; - return (*(range[0])>*(range[1]) || (val>=*(range[0]) && val<=*(range[1]))); -} -void RooNCSpline_1D_fast::setRangeValidity(const T valmin, const T valmax, const Int_t /*whichDirection*/){ - T* range[2]; - range[0] = &rangeXmin; - range[1] = &rangeXmax; - *(range[0])=valmin; - *(range[1])=valmax; -} -void RooNCSpline_1D_fast::cropValueForRange(T& val, const Int_t /*whichDirection*/)const{ - if (testRangeValidity(val)) return; - const T* range[2]; - range[0] = &rangeXmin; - range[1] = &rangeXmax; - if (val<*(range[0])) val = *(range[0]); - if (val>*(range[1])) val = *(range[1]); -} diff --git a/src/RooNCSpline_2D_fast.cc b/src/RooNCSpline_2D_fast.cc deleted file mode 100644 index 1067b075e88..00000000000 --- a/src/RooNCSpline_2D_fast.cc +++ /dev/null @@ -1,394 +0,0 @@ -#include "../interface/RooNCSpline_2D_fast.h" -#include -#include "TMath.h" -#include "Riostream.h" -#include "RooAbsReal.h" - -using namespace TMath; -using namespace RooFit; -using namespace std; - - -ClassImp(RooNCSpline_2D_fast) - -RooNCSpline_2D_fast::RooNCSpline_2D_fast() : - RooNCSplineCore(), - rangeYmin(1), rangeYmax(-1), - bcBeginX(RooNCSplineCore::bcNaturalSpline), bcEndX(RooNCSplineCore::bcNaturalSpline), - bcBeginY(RooNCSplineCore::bcNaturalSpline), bcEndY(RooNCSplineCore::bcNaturalSpline), - theYVar("theYVar", "theYVar", this) -{} - -RooNCSpline_2D_fast::RooNCSpline_2D_fast( - const char* name, - const char* title - ) : - RooNCSplineCore(name, title), - rangeYmin(1), rangeYmax(-1), - bcBeginX(RooNCSplineCore::bcNaturalSpline), bcEndX(RooNCSplineCore::bcNaturalSpline), - bcBeginY(RooNCSplineCore::bcNaturalSpline), bcEndY(RooNCSplineCore::bcNaturalSpline), - theYVar("theYVar", "theYVar", this) -{} - -RooNCSpline_2D_fast::RooNCSpline_2D_fast( - const char* name, - const char* title, - RooAbsReal& inXVar, - RooAbsReal& inYVar, - const std::vector& inXList, - const std::vector& inYList, - const std::vector>& inFcnList, - RooNCSplineCore::BoundaryCondition const bcBeginX_, - RooNCSplineCore::BoundaryCondition const bcEndX_, - RooNCSplineCore::BoundaryCondition const bcBeginY_, - RooNCSplineCore::BoundaryCondition const bcEndY_, - Bool_t inUseFloor, - T inFloorEval, - T inFloorInt - ) : - RooNCSplineCore(name, title, inXVar, inXList, inUseFloor, inFloorEval, inFloorInt), - rangeYmin(1), rangeYmax(-1), - bcBeginX(bcBeginX_), bcEndX(bcEndX_), - bcBeginY(bcBeginY_), bcEndY(bcEndY_), - theYVar("theYVar", "theYVar", this, inYVar), - YList(inYList), - FcnList(inFcnList) -{ - if (npointsX()>1 && npointsY()>1){ - // Prepare A and kappa arrays for x and y coordinates - int npoints; - Double_t det; - - vector> xA; getKappas(kappaX, 0); getAArray(kappaX, xA, bcBeginX, bcEndX); - npoints=kappaX.size(); - TMatrix_t xAtrans(npoints, npoints); - for (int i=0; i> yA; getKappas(kappaY, 1); getAArray(kappaY, yA, bcBeginY, bcEndY); - npoints=kappaY.size(); - TMatrix_t yAtrans(npoints, npoints); - for (int i=0; i>> coefsAlongY; // [Ax(y),Bx(y),Cx(y),Dx(y)][xbin][ybin] - int npoldim=0; - int nxbins=0; - for (unsigned int j=0; j> xcoefsAtYj = getCoefficientsPerY(kappaX, xAinv, j, bcBeginX, bcEndX, -1); // [ix][Ax,Bx,Cx,Dx] at each y_j - if (j==0){ - nxbins=xcoefsAtYj.size(); - npoldim=xcoefsAtYj.at(0).size(); - for (int ipow=0; ipow> dum_xycoefarray; - for (int ix=0; ix dum_ycoefarray; - dum_xycoefarray.push_back(dum_ycoefarray); - } - coefsAlongY.push_back(dum_xycoefarray); - } - } - if (nxbins!=(int)xcoefsAtYj.size() || npoldim!=(int)xcoefsAtYj.at(0).size()){ - coutE(InputArguments) << "RooNCSpline_2D_fast::interpolateFcn: nxbins!=(int)xcoefsAtYj.size() || npoldim!=(int)xcoefsAtYj.at(0).size()!" << endl; - assert(0); - } - for (int ix=0; ix>> xCoefs; - for (int ic=0; ic> yCoefs = getCoefficientsAlongDirection(kappaY, yAinv, coefsAlongY.at(ic).at(ix), bcBeginY, bcEndY, -1); // [iy][A,B,C,D] - xCoefs.push_back(yCoefs); - } - coefficients.push_back(xCoefs); - } - } - else assert(0); - - RooArgSet leafset; - getLeafDependents(theXVar, leafset); - getLeafDependents(theYVar, leafset); - addLeafDependents(leafset); - - emptyFcnList(); -} - -RooNCSpline_2D_fast::RooNCSpline_2D_fast( - const RooNCSpline_2D_fast& other, - const char* name - ) : - RooNCSplineCore(other, name), - rangeYmin(other.rangeYmin), rangeYmax(other.rangeYmax), - bcBeginX(other.bcBeginX), bcEndX(other.bcEndX), - bcBeginY(other.bcBeginY), bcEndY(other.bcEndY), - theYVar("theYVar", this, other.theYVar), - YList(other.YList), - FcnList(other.FcnList), - kappaX(other.kappaX), - kappaY(other.kappaY), - coefficients(other.coefficients) -{} - - -RooNCSplineCore::T RooNCSpline_2D_fast::interpolateFcn(Int_t code, const char* rangeName)const{ - DefaultAccumulator res; - - if (verbosity==RooNCSplineCore::kVerbose){ cout << "RooNCSpline_2D_fast(" << GetName() << ")::interpolateFcn begin with code: " << code << endl; } - - // Get bins - Int_t xbin=-1, xbinmin=-1, xbinmax=-1, ybin=-1, ybinmin=-1, ybinmax=-1; - RooNCSplineCore::T tx=0, txmin=0, txmax=0, ty=0, tymin=0, tymax=0; - if (code==0 || code%2!=0){ // Case to just compute the value at x - if (!testRangeValidity(theXVar, 0)) return 0; - xbin = getWhichBin(theXVar, 0); - tx = getTVar(kappaX, theXVar, xbin, 0); - } - else{ // Case to integrate along x - RooNCSplineCore::T coordmin = theXVar.min(rangeName); cropValueForRange(coordmin, 0); - RooNCSplineCore::T coordmax = theXVar.max(rangeName); cropValueForRange(coordmax, 0); - xbinmin = getWhichBin(coordmin, 0); - txmin = getTVar(kappaX, coordmin, xbinmin, 0); - xbinmax = getWhichBin(coordmax, 0); - txmax = getTVar(kappaX, coordmax, xbinmax, 0); - } - if (code==0 || code%3!=0){ // Case to just compute the value at y - if (!testRangeValidity(theYVar, 1)) return 0; - ybin = getWhichBin(theYVar, 1); - ty = getTVar(kappaY, theYVar, ybin, 1); - } - else{ // Case to integrate along y - RooNCSplineCore::T coordmin = theYVar.min(rangeName); cropValueForRange(coordmin, 1); - RooNCSplineCore::T coordmax = theYVar.max(rangeName); cropValueForRange(coordmax, 1); - ybinmin = getWhichBin(coordmin, 1); - tymin = getTVar(kappaY, coordmin, ybinmin, 1); - ybinmax = getWhichBin(coordmax, 1); - tymax = getTVar(kappaY, coordmax, ybinmax, 1); - } - - for (int ix=0; ix<(int)coefficients.size(); ix++){ - if ( - (xbin>=0 && ix!=xbin) - || - (xbinmin>=0 && xbinmax>=xbinmin && !(xbinmin<=ix && ix<=xbinmax)) - ) continue; - - RooNCSplineCore::T txlow=0, txhigh=1; - if (code>0 && code%2==0){ - if (ix==xbinmin) txlow=txmin; - if (ix==xbinmax) txhigh=txmax; - } - else txhigh=tx; - - if (verbosity==RooNCSplineCore::kVerbose){ - if (code==0 || code%2!=0) cout << "Evaluating tx=" << txhigh << " in bin " << ix << endl; - else cout << "Evaluating tx[" << txlow << ", " << txhigh << "] in bin " << ix << endl; - } - - // Get the x coefficients interpolated across y - vector xCoefs; - for (int ic=0; ic<(int)coefficients.at(ix).size(); ic++){ - const vector>& yCoefs = coefficients.at(ix).at(ic); - - if (verbosity==RooNCSplineCore::kVerbose) cout << "\tCoefficient " << ic << ":\n"; - - DefaultAccumulator theCoef; - for (int iy=0; iy<(int)yCoefs.size(); iy++){ - if ( - (ybin>=0 && iy!=ybin) - || - (ybinmin>=0 && ybinmax>=ybinmin && !(ybinmin<=iy && iy<=ybinmax)) - ) continue; - - RooNCSplineCore::T tylow=0, tyhigh=1; - if (code>0 && code%3==0){ - if (iy==ybinmin) tylow=tymin; - if (iy==ybinmax) tyhigh=tymax; - } - else tyhigh=ty; - - if (verbosity==RooNCSplineCore::kVerbose){ - if (code==0 || code%3!=0) cout << "\tEvaluating ty=" << tyhigh << " in bin " << iy << endl; - else cout << "\tEvaluating ty[" << tylow << ", " << tyhigh << "] in bin " << iy << endl; - } - - theCoef += evalSplineSegment(yCoefs.at(iy), kappaY.at(iy), tyhigh, tylow, (code>0 && code%3==0)); - } - - //if (code==0) cout << "\tCoefficient is " << theCoef << endl; - - xCoefs.push_back(theCoef); - } - - // Evaluate value of spline at x with coefficients evaluated at y - res += evalSplineSegment(xCoefs, kappaX.at(ix), txhigh, txlow, (code>0 && code%2==0)); - } - - return res; -} - -void RooNCSpline_2D_fast::getKappas(vector& kappas, const Int_t whichDirection){ - kappas.clear(); - RooNCSplineCore::T kappa=1; - - Int_t npoints; - vector const* coord; - if (whichDirection==0){ - npoints=npointsX(); - coord=&XList; - } - else{ - npoints=npointsY(); - coord=&YList; - } - - for (Int_t j=0; jat(j); - RooNCSplineCore::T val_jpo = coord->at(j+1); - RooNCSplineCore::T val_diff = (val_jpo-val_j); - if (fabs(val_diff)>RooNCSplineCore::T(0)) kappa = 1./val_diff; - else kappa = 0; - kappas.push_back(kappa); - } - kappas.push_back(kappa); // Push the same kappa_(N-1)=kappa_(N-2) at the end point -} -Int_t RooNCSpline_2D_fast::getWhichBin(const RooNCSplineCore::T& val, const Int_t whichDirection)const{ - Int_t bin=-1; - RooNCSplineCore::T valj, valjpo; - Int_t npoints; - vector const* coord; - if (whichDirection==0){ - coord=&XList; - npoints=npointsX(); - } - else{ - coord=&YList; - npoints=npointsY(); - } - - if (npoints<=1) bin=0; - else{ - valjpo = coord->at(0); - for (Int_t j=0; jat(j); - valjpo = coord->at(j+1); - if (val=valj){ bin=j; break; } - } - if (bin==-1 && val>=valjpo) bin=npoints-2; - else if (bin==-1) bin=0; - } - - return bin; -} -RooNCSplineCore::T RooNCSpline_2D_fast::getTVar(const vector& kappas, const RooNCSplineCore::T& val, const Int_t& bin, const Int_t whichDirection)const{ - const RooNCSplineCore::T& K=kappas.at(bin); - vector const* coord; - if (whichDirection==0) coord=&XList; - else coord=&YList; - return (val-coord->at(bin))*K; -} - -vector> RooNCSpline_2D_fast::getCoefficientsPerY(const std::vector& kappaX, const TMatrix_t& xAinv, const Int_t& ybin, RooNCSplineCore::BoundaryCondition const& bcBegin, RooNCSplineCore::BoundaryCondition const& bcEnd, const Int_t xbin)const{ - vector fcnList; - for (unsigned int bin=0; bin> coefs = getCoefficientsAlongDirection(kappaX, xAinv, fcnList, bcBegin, bcEnd, xbin); - return coefs; -} - -Double_t RooNCSpline_2D_fast::evaluate() const{ - Double_t value = interpolateFcn(0); - if (useFloor && value=RooNCSplineCore::kError) coutE(Eval) << "RooNCSpline_2D_fast ERROR::RooNCSpline_2D_fast(" << GetName() << ") evaluation returned " << value << " at (x, y) = (" << theXVar << ", " << theYVar << ")" << endl; - value = floorEval; - } - if (verbosity==RooNCSplineCore::kVerbose){ cout << "RooNCSpline_2D_fast(" << GetName() << ")::evaluate = " << value << " at (x, y) = (" << theXVar << ", " << theYVar << ")" << endl; } - return value; -} -Int_t RooNCSpline_2D_fast::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const{ - if (_forceNumInt) return 0; - - Int_t code=1; - - RooArgSet Xdeps, Ydeps; - RooRealVar* rrv_x = dynamic_cast(theXVar.absArg()); - RooRealVar* rrv_y = dynamic_cast(theYVar.absArg()); - if (rrv_x==0) theXVar.absArg()->leafNodeServerList(&Xdeps, 0, true); - if (rrv_y==0) theYVar.absArg()->leafNodeServerList(&Ydeps, 0, true); - - if (rrv_x!=0){ - if (Ydeps.find(*rrv_x)==0 || rrv_y!=0){ - if (matchArgs(allVars, analVars, theXVar)) code*=2; - } - } - if (rrv_y!=0){ - if (Xdeps.find(*rrv_y)==0 || rrv_x!=0){ - if (matchArgs(allVars, analVars, theYVar)) code*=3; - } - } - - if (code==1) code=0; - return code; -} -Double_t RooNCSpline_2D_fast::analyticalIntegral(Int_t code, const char* rangeName) const{ - Double_t value = interpolateFcn(code, rangeName); - if (useFloor && value=RooNCSplineCore::kError) coutE(Integration) << "RooNCSpline_2D_fast ERROR::RooNCSpline_2D_fast(" << GetName() << ") integration returned " << value << " for code = " << code << endl; - value = floorInt; - } - if (verbosity==RooNCSplineCore::kVerbose){ cout << "RooNCSpline_2D_fast(" << GetName() << ")::analyticalIntegral = " << value << " for code = " << code << endl; } - return value; -} - -Bool_t RooNCSpline_2D_fast::testRangeValidity(const T& val, const Int_t whichDirection) const{ - const T* range[2]; - if (whichDirection==0){ - range[0] = &rangeXmin; - range[1] = &rangeXmax; - } - else{ - range[0] = &rangeYmin; - range[1] = &rangeYmax; - } - return (*(range[0])>*(range[1]) || (val>=*(range[0]) && val<=*(range[1]))); -} -void RooNCSpline_2D_fast::setRangeValidity(const T valmin, const T valmax, const Int_t whichDirection){ - T* range[2]; - if (whichDirection==0){ - range[0] = &rangeXmin; - range[1] = &rangeXmax; - } - else{ - range[0] = &rangeYmin; - range[1] = &rangeYmax; - } - *(range[0])=valmin; - *(range[1])=valmax; -} -void RooNCSpline_2D_fast::cropValueForRange(T& val, const Int_t whichDirection)const{ - if (testRangeValidity(val, whichDirection)) return; - const T* range[2]; - if (whichDirection==0){ - range[0] = &rangeXmin; - range[1] = &rangeXmax; - } - else{ - range[0] = &rangeYmin; - range[1] = &rangeYmax; - } - if (val<*(range[0])) val = *(range[0]); - if (val>*(range[1])) val = *(range[1]); -} diff --git a/src/RooNCSpline_3D_fast.cc b/src/RooNCSpline_3D_fast.cc deleted file mode 100644 index 19ef55de6c3..00000000000 --- a/src/RooNCSpline_3D_fast.cc +++ /dev/null @@ -1,550 +0,0 @@ -#include "../interface/RooNCSpline_3D_fast.h" -#include -#include "TMath.h" -#include "Riostream.h" -#include "RooAbsReal.h" - -using namespace TMath; -using namespace RooFit; -using namespace std; - - -ClassImp(RooNCSpline_3D_fast) - -RooNCSpline_3D_fast::RooNCSpline_3D_fast() : - RooNCSplineCore(), - rangeYmin(1), rangeYmax(-1), - rangeZmin(1), rangeZmax(-1), - bcBeginX(RooNCSplineCore::bcNaturalSpline), bcEndX(RooNCSplineCore::bcNaturalSpline), - bcBeginY(RooNCSplineCore::bcNaturalSpline), bcEndY(RooNCSplineCore::bcNaturalSpline), - bcBeginZ(RooNCSplineCore::bcNaturalSpline), bcEndZ(RooNCSplineCore::bcNaturalSpline), - theYVar("theYVar", "theYVar", this), - theZVar("theZVar", "theZVar", this) -{} - -RooNCSpline_3D_fast::RooNCSpline_3D_fast( - const char* name, - const char* title - ) : - RooNCSplineCore(name, title), - rangeYmin(1), rangeYmax(-1), - rangeZmin(1), rangeZmax(-1), - bcBeginX(RooNCSplineCore::bcNaturalSpline), bcEndX(RooNCSplineCore::bcNaturalSpline), - bcBeginY(RooNCSplineCore::bcNaturalSpline), bcEndY(RooNCSplineCore::bcNaturalSpline), - bcBeginZ(RooNCSplineCore::bcNaturalSpline), bcEndZ(RooNCSplineCore::bcNaturalSpline), - theYVar("theYVar", "theYVar", this), - theZVar("theZVar", "theZVar", this) -{} - -RooNCSpline_3D_fast::RooNCSpline_3D_fast( - const char* name, - const char* title, - RooAbsReal& inXVar, - RooAbsReal& inYVar, - RooAbsReal& inZVar, - const std::vector& inXList, - const std::vector& inYList, - const std::vector& inZList, - const std::vector>>& inFcnList, - RooNCSplineCore::BoundaryCondition const bcBeginX_, - RooNCSplineCore::BoundaryCondition const bcEndX_, - RooNCSplineCore::BoundaryCondition const bcBeginY_, - RooNCSplineCore::BoundaryCondition const bcEndY_, - RooNCSplineCore::BoundaryCondition const bcBeginZ_, - RooNCSplineCore::BoundaryCondition const bcEndZ_, - Bool_t inUseFloor, - T inFloorEval, - T inFloorInt - ) : - RooNCSplineCore(name, title, inXVar, inXList, inUseFloor, inFloorEval, inFloorInt), - rangeYmin(1), rangeYmax(-1), - rangeZmin(1), rangeZmax(-1), - bcBeginX(bcBeginX_), bcEndX(bcEndX_), - bcBeginY(bcBeginY_), bcEndY(bcEndY_), - bcBeginZ(bcBeginZ_), bcEndZ(bcEndZ_), - theYVar("theYVar", "theYVar", this, inYVar), - theZVar("theZVar", "theZVar", this, inZVar), - YList(inYList), - ZList(inZList), - FcnList(inFcnList) -{ - if (npointsX()>1 && npointsY()>1 && npointsZ()>1){ - // Prepare A and kappa arrays for x, y and z coordinates - int npoints; - Double_t det; - - vector> xA; getKappas(kappaX, 0); getAArray(kappaX, xA, bcBeginX, bcEndX); - npoints=kappaX.size(); - TMatrix_t xAtrans(npoints, npoints); - for (int i=0; i> yA; getKappas(kappaY, 1); getAArray(kappaY, yA, bcBeginY, bcEndY); - npoints=kappaY.size(); - TMatrix_t yAtrans(npoints, npoints); - for (int i=0; i> zA; getKappas(kappaZ, 2); getAArray(kappaZ, zA, bcBeginZ, bcEndZ); - npoints=kappaZ.size(); - TMatrix_t zAtrans(npoints, npoints); - for (int i=0; i> - >> - > coefsAlongZ; // [ix][A_x,B_x,C_x,D_x][iy][A_x_y,B_x_y,C_x_y,D_x_y][z_k] at each z_k - for (unsigned int k=0; k> - >> coefficients_perZ; // [ix][A_x,B_x,C_x,D_x][iy][A_x_y,B_x_y,C_x_y,D_x_y] at each z_k - - vector>> coefsAlongY; // [xbin][Ax(y),Bx(y),Cx(y),Dx(y)][ybin] in each z - for (unsigned int j=0; j> xcoefsAtYjZk = getCoefficientsPerYPerZ(kappaX, xAinv, j, k, bcBeginX, bcEndX, -1); // [ix][Ax,Bx,Cx,Dx] at each y_j z_k - //cout << "\tCoefficients in y line " << j << " are found" << endl; - if (j==0){ - if (k==0){ - nxbins=xcoefsAtYjZk.size(); - nxpoldim=xcoefsAtYjZk.at(0).size(); - } - for (int ix=0; ix> dum_xycoefarray; - for (int icx=0; icx dum_ycoefarray; - dum_xycoefarray.push_back(dum_ycoefarray); - } - coefsAlongY.push_back(dum_xycoefarray); - } - } - if (nxbins!=(int)xcoefsAtYjZk.size() || nxpoldim!=(int)xcoefsAtYjZk.at(0).size()){ - coutE(InputArguments) << "RooNCSpline_3D_fast::interpolateFcn: nxbins!=(int)xcoefsAtYjZk.size() || nxpoldim!=(int)xcoefsAtYjZk.at(0).size()!" << endl; - assert(0); - } - for (int ix=0; ix>> xCoefs; - for (int icx=0; icx> yCoefs = getCoefficientsAlongDirection(kappaY, yAinv, coefsAlongY.at(ix).at(icx), bcBeginY, bcEndY, -1); // [iy][A,B,C,D] - xCoefs.push_back(yCoefs); - } - coefficients_perZ.push_back(xCoefs); - } - - if (k==0){ - nybins = coefficients_perZ.at(0).at(0).size(); - nypoldim = coefficients_perZ.at(0).at(0).at(0).size(); - for (int ix=0; ix>>> xCoefs; - for (int icx=0; icx>> xCoefsAlongY; - for (int iy=0; iy> yCoefs; - for (int icy=0; icy yCoefAtZj; - yCoefs.push_back(yCoefAtZj); - } - xCoefsAlongY.push_back(yCoefs); - } - xCoefs.push_back(xCoefsAlongY); - } - coefsAlongZ.push_back(xCoefs); - } - } - for (int ix=0; ix>>>> xCoefs; - for (int icx=0; icx>>> xCoefsAlongY; - for (int iy=0; iy>> yCoefs; - for (int icy=0; icy> yCoefsAlongZ = getCoefficientsAlongDirection(kappaZ, zAinv, coefsAlongZ.at(ix).at(icx).at(iy).at(icy), bcBeginZ, bcEndZ, -1); // [iz][A,B,C,D] - yCoefs.push_back(yCoefsAlongZ); - } - xCoefsAlongY.push_back(yCoefs); - } - xCoefs.push_back(xCoefsAlongY); - } - coefficients.push_back(xCoefs); - } - - } - else assert(0); - - RooArgSet leafset; - getLeafDependents(theXVar, leafset); - getLeafDependents(theYVar, leafset); - getLeafDependents(theZVar, leafset); - addLeafDependents(leafset); - - emptyFcnList(); -} - -RooNCSpline_3D_fast::RooNCSpline_3D_fast( - const RooNCSpline_3D_fast& other, - const char* name - ) : - RooNCSplineCore(other, name), - rangeYmin(other.rangeYmin), rangeYmax(other.rangeYmax), - rangeZmin(other.rangeZmin), rangeZmax(other.rangeZmax), - bcBeginX(other.bcBeginX), bcEndX(other.bcEndX), - bcBeginY(other.bcBeginY), bcEndY(other.bcEndY), - bcBeginZ(other.bcBeginZ), bcEndZ(other.bcEndZ), - theYVar("theYVar", this, other.theYVar), - theZVar("theZVar", this, other.theZVar), - YList(other.YList), - ZList(other.ZList), - FcnList(other.FcnList), - kappaX(other.kappaX), - kappaY(other.kappaY), - kappaZ(other.kappaZ), - coefficients(other.coefficients) -{} - -RooNCSplineCore::T RooNCSpline_3D_fast::interpolateFcn(Int_t code, const char* rangeName)const{ - DefaultAccumulator res; - - if (verbosity==RooNCSplineCore::kVerbose) cout << "RooNCSpline_3D_fast(" << GetName() << ")::interpolateFcn begin with code: " << code << endl; - - // Get bins - vector varprime; varprime.push_back(2); varprime.push_back(3); varprime.push_back(5); - const Int_t ndims=(const Int_t)varprime.size(); - vector varbin; - vector varbinmin; - vector varbinmax; - vector tvar; - vector tvarmin; - vector tvarmax; - vector*> varkappa; varkappa.push_back(&kappaX); varkappa.push_back(&kappaY); varkappa.push_back(&kappaZ); - vector varcoord; varcoord.push_back(&theXVar); varcoord.push_back(&theYVar); varcoord.push_back(&theZVar); - for (Int_t idim=0; idimmin(rangeName); cropValueForRange(coordmin, idim); - RooNCSplineCore::T coordmax = varcoord.at(idim)->max(rangeName); cropValueForRange(coordmax, idim); - const std::vector& kappadim = *(varkappa.at(idim)); - binmin = getWhichBin(coordmin, idim); - tmin = getTVar(kappadim, coordmin, binmin, idim); - binmax = getWhichBin(coordmax, idim); - tmax = getTVar(kappadim, coordmax, binmax, idim); - } - varbin.push_back(binval); - varbinmin.push_back(binmin); - varbinmax.push_back(binmax); - tvar.push_back(tval); - tvarmin.push_back(tmin); - tvarmax.push_back(tmax); - } - - for (int ix=0; ix<(int)coefficients.size(); ix++){ - if ( - (varbin[0]>=0 && ix!=varbin[0]) - || - (varbinmin[0]>=0 && varbinmax[0]>=varbinmin[0] && !(varbinmin[0]<=ix && ix<=varbinmax[0])) - ) continue; - RooNCSplineCore::T txlow=0, txhigh=1; - if (code>0 && code%varprime[0]==0){ - if (ix==varbinmin[0]) txlow=tvarmin[0]; - if (ix==varbinmax[0]) txhigh=tvarmax[0]; - } - else txhigh=tvar[0]; - // Get the x coefficients interpolated across y - vector xCoefs; - for (int icx=0; icx<(int)coefficients.at(ix).size(); icx++){ - DefaultAccumulator theXCoef; - for (int iy=0; iy<(int)coefficients.at(ix).at(icx).size(); iy++){ - if ( - (varbin[1]>=0 && iy!=varbin[1]) - || - (varbinmin[1]>=0 && varbinmax[1]>=varbinmin[1] && !(varbinmin[1]<=iy && iy<=varbinmax[1])) - ) continue; - RooNCSplineCore::T tylow=0, tyhigh=1; - if (code>0 && code%varprime[1]==0){ - if (iy==varbinmin[1]) tylow=tvarmin[1]; - if (iy==varbinmax[1]) tyhigh=tvarmax[1]; - } - else tyhigh=tvar[1]; - // Get the y coefficients interpolated across z - vector yCoefs; - for (int icy=0; icy<(int)coefficients.at(ix).at(icx).at(iy).size(); icy++){ - DefaultAccumulator theYCoef; - for (int iz=0; iz<(int)coefficients.at(ix).at(icx).at(iy).at(icy).size(); iz++){ - if ( - (varbin[2]>=0 && iz!=varbin[2]) - || - (varbinmin[2]>=0 && varbinmax[2]>=varbinmin[2] && !(varbinmin[2]<=iz && iz<=varbinmax[2])) - ) continue; - - RooNCSplineCore::T tzlow=0, tzhigh=1; - if (code>0 && code%varprime[2]==0){ - if (iz==varbinmin[2]) tzlow=tvarmin[2]; - if (iz==varbinmax[2]) tzhigh=tvarmax[2]; - } - else tzhigh=tvar[2]; - - theYCoef += evalSplineSegment(coefficients.at(ix).at(icx).at(iy).at(icy).at(iz), varkappa[2]->at(iz), tzhigh, tzlow, (code>0 && code%varprime[2]==0)); - } - yCoefs.push_back(theYCoef); - } - theXCoef += evalSplineSegment(yCoefs, varkappa[1]->at(iy), tyhigh, tylow, (code>0 && code%varprime[1]==0)); - } - xCoefs.push_back(theXCoef); - } - // Evaluate value of spline at x with coefficients evaluated at y - res += evalSplineSegment(xCoefs, varkappa[0]->at(ix), txhigh, txlow, (code>0 && code%varprime[0]==0)); - } - - return res; -} - -void RooNCSpline_3D_fast::getKappas(vector& kappas, const Int_t whichDirection){ - kappas.clear(); - RooNCSplineCore::T kappa=1; - - Int_t npoints; - vector const* coord; - if (whichDirection==0){ - npoints=npointsX(); - coord=&XList; - } - else if (whichDirection==1){ - npoints=npointsY(); - coord=&YList; - } - else{ - npoints=npointsZ(); - coord=&ZList; - } - - for (Int_t j=0; jat(j); - RooNCSplineCore::T val_jpo = coord->at(j+1); - RooNCSplineCore::T val_diff = (val_jpo-val_j); - if (fabs(val_diff)>RooNCSplineCore::T(0)) kappa = 1./val_diff; - else kappa = 0; - kappas.push_back(kappa); - } - kappas.push_back(kappa); // Push the same kappa_(N-1)=kappa_(N-2) at the end point -} -Int_t RooNCSpline_3D_fast::getWhichBin(const RooNCSplineCore::T& val, const Int_t whichDirection)const{ - Int_t bin=-1; - RooNCSplineCore::T valj, valjpo; - Int_t npoints; - vector const* coord; - if (whichDirection==0){ - npoints=npointsX(); - coord=&XList; - } - else if (whichDirection==1){ - npoints=npointsY(); - coord=&YList; - } - else{ - npoints=npointsZ(); - coord=&ZList; - } - - if (npoints<=1) bin=0; - else{ - valjpo = coord->at(0); - for (Int_t j=0; jat(j); - valjpo = coord->at(j+1); - if (val=valj){ bin=j; break; } - } - if (bin==-1 && val>=valjpo) bin=npoints-2; - else if (bin==-1) bin=0; - } - - return bin; -} -RooNCSplineCore::T RooNCSpline_3D_fast::getTVar(const vector& kappas, const RooNCSplineCore::T& val, const Int_t& bin, const Int_t whichDirection)const{ - const RooNCSplineCore::T& K=kappas.at(bin); - vector const* coord; - if (whichDirection==0) coord=&XList; - else if (whichDirection==1) coord=&YList; - else coord=&ZList; - return (val-coord->at(bin))*K; -} - -vector> RooNCSpline_3D_fast::getCoefficientsPerYPerZ( - const std::vector& kappaX, const TMatrix_t& xAinv, - const Int_t& ybin, const Int_t& zbin, - RooNCSplineCore::BoundaryCondition const& bcBegin, RooNCSplineCore::BoundaryCondition const& bcEnd, - const Int_t xbin - )const{ - vector fcnList; - for (unsigned int bin=0; bin> coefs = getCoefficientsAlongDirection(kappaX, xAinv, fcnList, bcBegin, bcEnd, xbin); - return coefs; -} - -Double_t RooNCSpline_3D_fast::evaluate() const{ - Double_t value = interpolateFcn(0); - if (useFloor && value=RooNCSplineCore::kError) coutE(Eval) << "RooNCSpline_3D_fast ERROR::RooNCSpline_3D_fast(" << GetName() << ") evaluation returned " << value << " at (x, y, z) = (" << theXVar << ", " << theYVar << ", " << theZVar << ")" << endl; - value = floorEval; - } - if (verbosity==RooNCSplineCore::kVerbose){ cout << "RooNCSpline_3D_fast(" << GetName() << ")::evaluate = " << value << " at (x, y, z) = (" << theXVar << ", " << theYVar << ", " << theZVar << ")" << endl; } - return value; -} -Int_t RooNCSpline_3D_fast::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const{ - if (_forceNumInt) return 0; - - Int_t code=1; - - RooArgSet Xdeps, Ydeps, Zdeps; - RooRealVar* rrv_x = dynamic_cast(theXVar.absArg()); - RooRealVar* rrv_y = dynamic_cast(theYVar.absArg()); - RooRealVar* rrv_z = dynamic_cast(theZVar.absArg()); - if (rrv_x==0) theXVar.absArg()->leafNodeServerList(&Xdeps, 0, true); - if (rrv_y==0) theYVar.absArg()->leafNodeServerList(&Ydeps, 0, true); - if (rrv_z==0) theZVar.absArg()->leafNodeServerList(&Zdeps, 0, true); - - if (rrv_x!=0){ - if ( - (Ydeps.find(*rrv_x)==0 || rrv_y!=0) - && - (Zdeps.find(*rrv_x)==0 || rrv_z!=0) - ){ - if (matchArgs(allVars, analVars, theXVar)) code*=2; - } - } - if (rrv_y!=0){ - if ( - (Xdeps.find(*rrv_y)==0 || rrv_x!=0) - && - (Zdeps.find(*rrv_y)==0 || rrv_z!=0) - ){ - if (matchArgs(allVars, analVars, theYVar)) code*=3; - } - } - if (rrv_z!=0){ - if ( - (Xdeps.find(*rrv_z)==0 || rrv_x!=0) - && - (Ydeps.find(*rrv_z)==0 || rrv_y!=0) - ){ - if (matchArgs(allVars, analVars, theZVar)) code*=5; - } - } - - if (code==1) code=0; - return code; -} -Double_t RooNCSpline_3D_fast::analyticalIntegral(Int_t code, const char* rangeName) const{ - Double_t value = interpolateFcn(code, rangeName); - if (useFloor && value=RooNCSplineCore::kError) coutE(Integration) << "RooNCSpline_3D_fast ERROR::RooNCSpline_3D_fast(" << GetName() << ") integration returned " << value << " for code = " << code << endl; - value = floorInt; - } - if (verbosity==RooNCSplineCore::kVerbose){ cout << "RooNCSpline_3D_fast(" << GetName() << ")::analyticalIntegral = " << value << " for code = " << code << endl; } - return value; -} - -Bool_t RooNCSpline_3D_fast::testRangeValidity(const T& val, const Int_t whichDirection) const{ - const T* range[2]; - if (whichDirection==0){ - range[0] = &rangeXmin; - range[1] = &rangeXmax; - } - else if (whichDirection==1){ - range[0] = &rangeYmin; - range[1] = &rangeYmax; - } - else{ - range[0] = &rangeZmin; - range[1] = &rangeZmax; - } - return (*(range[0])>*(range[1]) || (val>=*(range[0]) && val<=*(range[1]))); -} -void RooNCSpline_3D_fast::setRangeValidity(const T valmin, const T valmax, const Int_t whichDirection){ - T* range[2]; - if (whichDirection==0){ - range[0] = &rangeXmin; - range[1] = &rangeXmax; - } - else if (whichDirection==1){ - range[0] = &rangeYmin; - range[1] = &rangeYmax; - } - else{ - range[0] = &rangeZmin; - range[1] = &rangeZmax; - } - *(range[0])=valmin; - *(range[1])=valmax; -} -void RooNCSpline_3D_fast::cropValueForRange(T& val, const Int_t whichDirection)const{ - if (testRangeValidity(val, whichDirection)) return; - const T* range[2]; - if (whichDirection==0){ - range[0] = &rangeXmin; - range[1] = &rangeXmax; - } - else if (whichDirection==1){ - range[0] = &rangeYmin; - range[1] = &rangeYmax; - } - else{ - range[0] = &rangeZmin; - range[1] = &rangeZmax; - } - if (val<*(range[0])) val = *(range[0]); - if (val>*(range[1])) val = *(range[1]); -} diff --git a/src/RooPiecewisePolynomial.cc b/src/RooPiecewisePolynomial.cc deleted file mode 100644 index d3b9557dd05..00000000000 --- a/src/RooPiecewisePolynomial.cc +++ /dev/null @@ -1,215 +0,0 @@ -#include "../interface/RooPiecewisePolynomial.h" -#include -#include -#include -#include "TVectorD.h" -#include "TMatrixD.h" -#include "TIterator.h" - - -using namespace std; - - -RooPiecewisePolynomial::RooPiecewisePolynomial(const int nfcn_, const int polyndof_) : - RooAbsReal(), - xvar("xvar", "xvar", this), - parList("parList", "parList", this), - nfcn(nfcn_), polyndof(polyndof_), - nnodes(nfcn-1), // How many nodes are in between - ndof_endfcn(polyndof-1), // 1 degree for the function value - ndof_middlefcn(polyndof-2) // +1 for slope of the other node -{ - assert((nfcn>2 && polyndof>=2) || (nfcn==2 && polyndof>=1) || (nfcn==1 && polyndof>=0)); -} -RooPiecewisePolynomial::RooPiecewisePolynomial(const char* name, const char* title, const int nfcn_, const int polyndof_) : - RooAbsReal(name, title), - xvar("xvar", "xvar", this), - parList("parList", "parList", this), - nfcn(nfcn_), polyndof(polyndof_), - nnodes(nfcn-1), // How many nodes are in between - ndof_endfcn(polyndof-1), // 1 degree for the function value - ndof_middlefcn(polyndof-2) // +1 for slope of the other node -{ - assert((nfcn>2 && polyndof>=2) || (nfcn==2 && polyndof>=1) || (nfcn==1 && polyndof>=0)); -} -RooPiecewisePolynomial::RooPiecewisePolynomial(const char* name, const char* title, RooAbsReal& xvar_, RooArgList const& parList_, const int nfcn_, const int polyndof_) : - RooAbsReal(name, title), - xvar("xvar", "xvar", this, xvar_), - parList("parList", "parList", this), - nfcn(nfcn_), polyndof(polyndof_), - nnodes(nfcn-1), // How many nodes are in between - ndof_endfcn(polyndof-1), // 1 degree for the function value - ndof_middlefcn(polyndof-2) // +1 for slope of the other node -{ - assert((nfcn>2 && polyndof>=2) || (nfcn==2 && polyndof>=1) || (nfcn==1 && polyndof>=0)); - - TIterator* coefIter = parList_.createIterator(); - RooAbsArg* func; - while ((func = (RooAbsArg*) coefIter->Next())) { - if (!dynamic_cast(func)) { - cerr << "RooPiecewisePolynomial::RooPiecewisePolynomial(" << GetName() << ") funcficient " << func->GetName() << " is not of type RooAbsReal" << endl; - assert(0); - } - parList.add(*func); - } - delete coefIter; -} -RooPiecewisePolynomial::RooPiecewisePolynomial(RooPiecewisePolynomial const& other, const char* name) : - RooAbsReal(other, name), - xvar("xvar", this, other.xvar), - parList("parList", this, other.parList), - nfcn(other.nfcn), polyndof(other.polyndof), - nnodes(other.nnodes), - ndof_endfcn(other.ndof_endfcn), - ndof_middlefcn(other.ndof_middlefcn) -{} - -double RooPiecewisePolynomial::eval(double x, std::vector const& par)const{ - // If we say the form of the polynomial is [0] + [1]*x + [2]*x2 + [3]*x3 + [4]*x4..., - // use the highest two orders for matching at the nodes and free the rest. - const double d_epsilon = 0; - if ( - (nfcn>=2 && (int) par.size()!=2*ndof_endfcn+(nfcn-2)*ndof_middlefcn+nnodes) - || - (nfcn==1 && (int) par.size()!=polyndof) - || - nfcn<=0 - ) assert(0); - - if (nfcn==1){ - double res = 0; - for (int ip=0; ip node(nnodes, 0); // First [0,...,nnodes-1] parameters are nodes - vector> pars_full(nfcn, vector(polyndof, 0)); // The full coefficients array - - for (int ip=0; ipnode[ip2]) return d_epsilon; - } - } - int pos_ctr = nnodes; - for (int index=0; index> xton(nnodes, vector(polyndof, 0)); // Array of node^power - vector> nxtom(nnodes, vector(polyndof, 0)); // Array of power*node^(power-1) - for (int inode=0; inode ysbar_nodes(2*nnodes, 0); - vector> coeff_ysbar(2*nnodes, vector(2*nnodes, 0)); - int cstart=-1; - for (int inode=0; inode=0){ - coeff_ysbar[inode][cstart] = -sign_i*xton[inode][polyndof-2]; - coeff_ysbar[nnodes + inode][cstart] = -sign_i*nxtom[inode][polyndof-2]; - } - coeff_ysbar[inode][cstart+1] = -sign_i*xton[inode][polyndof-1]; - coeff_ysbar[nnodes + inode][cstart+1] = -sign_i*nxtom[inode][polyndof-1]; - coeff_ysbar[inode][cstart+2] = -sign_j*xton[inode][polyndof-2]; - coeff_ysbar[nnodes + inode][cstart+2] = -sign_j*nxtom[inode][polyndof-2]; - if ((cstart+3)<2*nnodes){ - coeff_ysbar[inode][cstart+3] = -sign_j*xton[inode][polyndof-1]; - coeff_ysbar[nnodes + inode][cstart+3] = -sign_j*nxtom[inode][polyndof-1]; - } - cstart+=2; - } - - TVectorD polyvec(2*nnodes, ysbar_nodes.data()); - TMatrixD polycoeff(2*nnodes, 2*nnodes); - for (int i=0; i<2*nnodes; i++){ for (int j=0; j<2*nnodes; j++) polycoeff(i, j)=coeff_ysbar[i][j]; } - double testdet=0; - TMatrixD polycoeff_inv = polycoeff.Invert(&testdet); - if (testdet!=0){ - TVectorD unknowncoeffs = polycoeff_inv*polyvec; - pos_ctr=0; - for (int index=0; index=node[index] && x=node[nnodes-1]) index_chosen = nfcn-1; - - double res = 0; - for (int ip=0; ip par; par.reserve(parList.getSize()); - for (int ip=0; ip(parList.at(ip)))->getVal()); - return eval(xvar, par); -} -double RooPiecewisePolynomial::evaluate(double* x, double* p)const{ // For calling in a TF1 object - // No size check is done, so be careful! - unsigned int psize=2*ndof_endfcn+(nfcn-2)*ndof_middlefcn+nnodes; - //cout << "RooPiecewisePolynomial::evaluate: N parameters = " << psize << endl; - std::vector par; par.reserve(psize); - for (unsigned int ip=0; ip -#include -#include - -using namespace std; -using namespace RooFit; -using namespace TMath; - - -ClassImp(RooRealFlooredSumPdf) - - - -//_____________________________________________________________________________ -RooRealFlooredSumPdf::RooRealFlooredSumPdf() -{ - // Default constructor - // coverity[UNINIT_CTOR] - _funcIter = _funcList.createIterator() ; - _coefIter = _coefList.createIterator() ; - _extended = kFALSE ; - _doFloor = kTRUE ; - _floorVal = 1e-100 ; -} - - - -//_____________________________________________________________________________ -RooRealFlooredSumPdf::RooRealFlooredSumPdf(const char *name, const char *title) : - RooAbsPdf(name,title), - _normIntMgr(this,10), - _haveLastCoef(kFALSE), - _funcList("!funcList","List of functions",this), - _coefList("!coefList","List of coefficients",this), - _extended(kFALSE), - _doFloor(kTRUE), - _floorVal(1e-100) -{ - // Constructor with name and title - _funcIter = _funcList.createIterator() ; - _coefIter = _coefList.createIterator() ; -} - - - -//_____________________________________________________________________________ -RooRealFlooredSumPdf::RooRealFlooredSumPdf(const char *name, const char *title, const RooArgList& inFuncList, const RooArgList& inCoefList, Bool_t extended) : -RooAbsPdf(name, title), -_normIntMgr(this, 10), -_haveLastCoef(kFALSE), -_funcList("!funcList", "List of functions", this), -_coefList("!coefList", "List of coefficients", this), -_extended(extended), -_doFloor(kTRUE), -_floorVal(1e-100) -{ - // Constructor p.d.f implementing sum_i [ coef_i * func_i ], if N_coef==N_func - // or sum_i [ coef_i * func_i ] + (1 - sum_i [ coef_i ] )* func_N if Ncoef==N_func-1 - // - // All coefficients and functions are allowed to be negative - // but the sum is not, which is enforced at runtime. - - if (!(inFuncList.getSize() == inCoefList.getSize() + 1 || inFuncList.getSize() == inCoefList.getSize())) { - coutE(InputArguments) << "RooRealFlooredSumPdf::RooRealFlooredSumPdf(" << GetName() - << ") number of pdfs and coefficients inconsistent, must have Nfunc=Ncoef or Nfunc=Ncoef+1" << endl; - assert(0); - } - - _funcIter = _funcList.createIterator(); - _coefIter = _coefList.createIterator(); - - // Constructor with N functions and N or N-1 coefs - TIterator* funcIter = inFuncList.createIterator(); - TIterator* coefIter = inCoefList.createIterator(); - RooAbsArg* func; - RooAbsArg* coef; - - while ((coef = (RooAbsArg*)coefIter->Next())) { - func = (RooAbsArg*)funcIter->Next(); - - if (!dynamic_cast(coef)) { - coutW(InputArguments) << "RooRealFlooredSumPdf::RooRealFlooredSumPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal, ignored" << endl; - continue; - } - if (!dynamic_cast(func)) { - coutW(InputArguments) << "RooRealFlooredSumPdf::RooRealFlooredSumPdf(" << GetName() << ") func " << func->GetName() << " is not of type RooAbsReal, ignored" << endl; - continue; - } - _funcList.add(*func); - _coefList.add(*coef); - } - - func = (RooAbsReal*)funcIter->Next(); - if (func) { - if (!dynamic_cast(func)) { - coutE(InputArguments) << "RooRealFlooredSumPdf::RooRealFlooredSumPdf(" << GetName() << ") last func " << func->GetName() << " is not of type RooAbsReal, fatal error" << endl; - assert(0); - } - _funcList.add(*func); - } - else _haveLastCoef = kTRUE; - - - delete funcIter; - delete coefIter; -} - - - - -//_____________________________________________________________________________ -RooRealFlooredSumPdf::RooRealFlooredSumPdf(const RooRealFlooredSumPdf& other, const char* name) : -RooAbsPdf(other, name), -_normIntMgr(other._normIntMgr, this), -_haveLastCoef(other._haveLastCoef), -_funcList("!funcList", this, other._funcList), -_coefList("!coefList", this, other._coefList), -_extended(other._extended), -_doFloor(other._doFloor), -_floorVal(other._floorVal) -{ - // Copy constructor - - _funcIter = _funcList.createIterator(); - _coefIter = _coefList.createIterator(); -} - - -//_____________________________________________________________________________ -void RooRealFlooredSumPdf::setFloor(Double_t val) -{ - _floorVal = val; -} - - - -//_____________________________________________________________________________ -RooRealFlooredSumPdf::~RooRealFlooredSumPdf() -{ - // Destructor - delete _funcIter; - delete _coefIter; -} - - - - - -//_____________________________________________________________________________ -RooAbsPdf::ExtendMode RooRealFlooredSumPdf::extendMode() const -{ - return (_extended && (_funcList.getSize() == _coefList.getSize())) ? CanBeExtended : CanNotBeExtended; -} - - - - -//_____________________________________________________________________________ -Double_t RooRealFlooredSumPdf::evaluate() const -{ - // Calculate the current value - - Double_t value(0); - - // Do running sum of coef/func pairs, calculate lastCoef. - RooFIter funcIter = _funcList.fwdIterator(); - RooFIter coefIter = _coefList.fwdIterator(); - RooAbsReal* coef; - RooAbsReal* func; - - // N funcs, N-1 coefficients - Double_t lastCoef(1); - while ((coef = (RooAbsReal*)coefIter.next())) { - func = (RooAbsReal*)funcIter.next(); - Double_t coefVal = coef->getVal(); - if (coefVal) { - cxcoutD(Eval) << "RooRealFlooredSumPdf::eval(" << GetName() << ") coefVal = " << coefVal << " funcVal = " << func->ClassName() << "::" << func->GetName() << " = " << func->getVal() << endl; - value += func->getVal()*coefVal; - lastCoef -= coef->getVal(); - } - } - - if (!_haveLastCoef) { - // Add last func with correct coefficient - func = (RooAbsReal*)funcIter.next(); - value += func->getVal()*lastCoef; - - cxcoutD(Eval) << "RooRealFlooredSumPdf::eval(" << GetName() << ") lastCoef = " << lastCoef << " funcVal = " << func->getVal() << endl; - - // Warn about coefficient degeneration - if (lastCoef<0 || lastCoef>1) { - coutW(Eval) << "RooRealFlooredSumPdf::evaluate(" << GetName() - << " WARNING: sum of FUNC coefficients not in range [0-1], value=" - << 1 - lastCoef << endl; - } - } - - // Introduce floor - if (value <= 0. && _doFloor) value = _floorVal; - - return value; -} - - - - -//_____________________________________________________________________________ -Bool_t RooRealFlooredSumPdf::checkObservables(const RooArgSet* nset) const -{ - // Check if FUNC is valid for given normalization set. - // Coeffient and FUNC must be non-overlapping, but func-coefficient - // pairs may overlap each other - // - // In the present implementation, coefficients may not be observables or derive - // from observables - - Bool_t ret(kFALSE); - - _funcIter->Reset(); - _coefIter->Reset(); - RooAbsReal* coef; - RooAbsReal* func; - while ((coef = (RooAbsReal*)_coefIter->Next())) { - func = (RooAbsReal*)_funcIter->Next(); - if (func->observableOverlaps(nset, *coef)) { - coutE(InputArguments) << "RooRealFlooredSumPdf::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->GetName() - << " and FUNC " << func->GetName() << " have one or more observables in common" << endl; - ret = kTRUE; - } - if (coef->dependsOn(*nset)) { - coutE(InputArguments) << "RooRealPdf::checkObservables(" << GetName() << "): ERROR coefficient " << coef->GetName() - << " depends on one or more of the following observables"; nset->Print("1"); - ret = kTRUE; - } - } - - return ret; -} - - - - -//_____________________________________________________________________________ -Int_t RooRealFlooredSumPdf::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet* normSet2, const char* rangeName) const -{ - //cout << "RooRealFlooredSumPdf::getAnalyticalIntegralWN:"<") << endl; - // Advertise that all integrals can be handled internally. - - // Handle trivial no-integration scenario - if (allVars.getSize() == 0) return 0; - if (_forceNumInt) return 0; - - // Select subset of allVars that are actual dependents - analVars.add(allVars); - RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0; - - - // Check if this configuration was created before - Int_t sterileIdx(-1); - CacheElem* cache = (CacheElem*)_normIntMgr.getObj(normSet, &analVars, &sterileIdx, RooNameReg::ptr(rangeName)); - if (cache) { - //cout << "RooRealFlooredSumPdf("<") << " -> " << _normIntMgr.lastIndex()+1 << " (cached)" << endl; - return _normIntMgr.lastIndex() + 1; - } - - // Create new cache element - cache = new CacheElem; - - // Make list of function projection and normalization integrals - _funcIter->Reset(); - RooAbsReal *func; - while ((func = (RooAbsReal*)_funcIter->Next())) { - RooAbsReal* funcInt = func->createIntegral(analVars, rangeName); - cache->_funcIntList.addOwned(*funcInt); - if (normSet && normSet->getSize() > 0) { - RooAbsReal* funcNorm = func->createIntegral(*normSet); - cache->_funcNormList.addOwned(*funcNorm); - } - } - - // Store cache element - Int_t code = _normIntMgr.setObj(normSet, &analVars, (RooAbsCacheElement*)cache, RooNameReg::ptr(rangeName)); - - if (normSet) delete normSet; - - //cout << "RooRealFlooredSumPdf("<") << " -> " << code+1 << endl; - return code + 1; -} - - - - -//_____________________________________________________________________________ -Double_t RooRealFlooredSumPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet2, const char* rangeName) const -{ - //cout << "RooRealFlooredSumPdf::analyticalIntegralWN:"<") << endl; - // Implement analytical integrations by deferring integration of component - // functions to integrators of components - - // Handle trivial passthrough scenario - if (code == 0) return getVal(normSet2); - - - // WVE needs adaptation for rangeName feature - CacheElem* cache = (CacheElem*)_normIntMgr.getObjByIndex(code - 1); - if (cache == 0) { // revive the (sterilized) cache - //cout << "RooRealFlooredSumPdf("<") << ": reviving cache "<< endl; - std::unique_ptr vars(getParameters(RooArgSet())); - RooArgSet dummy; -#if ROOT_VERSION_CODE < ROOT_VERSION(6,26,0) - std::unique_ptr iset(_normIntMgr.nameSet2ByIndex(code - 1)->select(*vars)); - std::unique_ptr nset(_normIntMgr.nameSet1ByIndex(code - 1)->select(*vars)); - Int_t code2 = getAnalyticalIntegralWN(*iset, dummy, nset.get(), rangeName); -#else - // In ROOT 6.26, the RooNameSet was removed and the "selectFromSet*" - // functions were introduced to replace its functionality - RooArgSet iset{_normIntMgr.selectFromSet2(*vars, code - 1)}; - RooArgSet nset{_normIntMgr.selectFromSet1(*vars, code - 1)}; - Int_t code2 = getAnalyticalIntegralWN(iset, dummy, &nset, rangeName); -#endif - assert(code == code2); // must have revived the right (sterilized) slot... - cache = (CacheElem*)_normIntMgr.getObjByIndex(code - 1); - assert(cache != 0); - } - - RooFIter funcIntIter = cache->_funcIntList.fwdIterator(); - RooFIter coefIter = _coefList.fwdIterator(); - RooFIter funcIter = _funcList.fwdIterator(); - RooAbsReal *coef(0), *funcInt(0), *func(0); - Double_t value(0); - - // N funcs, N-1 coefficients - Double_t lastCoef(1); - while ((coef = (RooAbsReal*)coefIter.next())) { - funcInt = (RooAbsReal*)funcIntIter.next(); - func = (RooAbsReal*)funcIter.next(); - Double_t coefVal = coef->getVal(normSet2); - if (coefVal) { - assert(func); - assert(funcInt); - value += funcInt->getVal()*coefVal; - lastCoef -= coef->getVal(normSet2); - } - } - - if (!_haveLastCoef) { - // Add last func with correct coefficient - funcInt = (RooAbsReal*)funcIntIter.next(); - assert(funcInt); - value += funcInt->getVal()*lastCoef; - - // Warn about coefficient degeneration - if (lastCoef<0 || lastCoef>1) { - coutW(Eval) << "RooRealFlooredSumPdf::integral(" << GetName() - << ") WARNING: Sum of FUNC coefficients not in range [0-1], value=" - << 1 - lastCoef << endl; - } - } - - Double_t normVal(1); - if (normSet2 && normSet2->getSize() > 0) { - normVal = 0; - - // N funcs, N-1 coefficients - RooAbsReal* funcNorm; - RooFIter funcNormIter = cache->_funcNormList.fwdIterator(); - RooFIter coefIter2 = _coefList.fwdIterator(); - while ((coef = (RooAbsReal*)coefIter2.next())) { - funcNorm = (RooAbsReal*)funcNormIter.next(); - Double_t coefVal = coef->getVal(normSet2); - if (coefVal) { - assert(funcNorm); - normVal += funcNorm->getVal()*coefVal; - } - } - - // Add last func with correct coefficient - if (!_haveLastCoef) { - funcNorm = (RooAbsReal*)funcNormIter.next(); - assert(funcNorm); - normVal += funcNorm->getVal()*lastCoef; - } - } - - Double_t result = 0; - if(normVal>0) result = value / normVal; - if (result<=0. && _doFloor){ - coutW(Eval) << "RooRealFlooredSumPdf::integral(" << GetName() - << ") WARNING: Integral " << result << " below threshold." << endl; - result = _floorVal; - } - return result; -} - - -//_____________________________________________________________________________ -Double_t RooRealFlooredSumPdf::expectedEvents(const RooArgSet* nset) const -{ - Double_t n = getNorm(nset); - if (n < 0) { - logEvalError("Expected number of events is negative"); - } - return n; -} - - -//_____________________________________________________________________________ -std::list* RooRealFlooredSumPdf::binBoundaries(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const -{ - - list* sumBinB = 0; - Bool_t needClean(kFALSE); - - RooFIter iter = _funcList.fwdIterator(); - RooAbsReal* func; - // Loop over components pdf - while ((func = (RooAbsReal*)iter.next())) { - - list* funcBinB = func->binBoundaries(obs, xlo, xhi); - - // Process hint - if (funcBinB) { - if (!sumBinB) { - // If this is the first hint, then just save it - sumBinB = funcBinB; - } - else { - - list* newSumBinB = new list(sumBinB->size() + funcBinB->size()); - - // Merge hints into temporary array - merge(funcBinB->begin(), funcBinB->end(), sumBinB->begin(), sumBinB->end(), newSumBinB->begin()); - - // Copy merged array without duplicates to new sumBinBArrau - delete sumBinB; - delete funcBinB; - sumBinB = newSumBinB; - needClean = kTRUE; - } - } - } - - // Remove consecutive duplicates - if (needClean) { - list::iterator new_end = unique(sumBinB->begin(), sumBinB->end()); - sumBinB->erase(new_end, sumBinB->end()); - } - - return sumBinB; -} - - - -//_____________________________________________________________________________B -Bool_t RooRealFlooredSumPdf::isBinnedDistribution(const RooArgSet& obs) const -{ - // If all components that depend on obs are binned that so is the product - - RooFIter iter = _funcList.fwdIterator(); - RooAbsReal* func; - while ((func = (RooAbsReal*)iter.next())) { - if (func->dependsOn(obs) && !func->isBinnedDistribution(obs)) { - return kFALSE; - } - } - - return kTRUE; -} - - - - - -//_____________________________________________________________________________ -std::list* RooRealFlooredSumPdf::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const -{ - list* sumHint = 0; - Bool_t needClean(kFALSE); - - RooFIter iter = _funcList.fwdIterator(); - RooAbsReal* func; - // Loop over components pdf - while ((func = (RooAbsReal*)iter.next())) { - - list* funcHint = func->plotSamplingHint(obs, xlo, xhi); - - // Process hint - if (funcHint) { - if (!sumHint) { - - // If this is the first hint, then just save it - sumHint = funcHint; - - } - else { - - list* newSumHint = new list(sumHint->size() + funcHint->size()); - - // Merge hints into temporary array - merge(funcHint->begin(), funcHint->end(), sumHint->begin(), sumHint->end(), newSumHint->begin()); - - // Copy merged array without duplicates to new sumHintArrau - delete sumHint; - sumHint = newSumHint; - needClean = kTRUE; - } - } - } - - // Remove consecutive duplicates - if (needClean) { - list::iterator new_end = unique(sumHint->begin(), sumHint->end()); - sumHint->erase(new_end, sumHint->end()); - } - - return sumHint; -} - - - - -//_____________________________________________________________________________ -void RooRealFlooredSumPdf::printMetaArgs(ostream& os) const -{ - // Customized printing of arguments of a RooRealFlooredSumPdf to more intuitively reflect the contents of the - // product operator construction - - _funcIter->Reset(); - _coefIter->Reset(); - - Bool_t first(kTRUE); - - RooAbsArg* coef, *func; - if (_coefList.getSize() != 0) { - while ((coef = (RooAbsArg*)_coefIter->Next())) { - if (!first) { - os << " + "; - } - else { - first = kFALSE; - } - func = (RooAbsArg*)_funcIter->Next(); - os << coef->GetName() << " * " << func->GetName(); - } - func = (RooAbsArg*)_funcIter->Next(); - if (func) { - os << " + [%] * " << func->GetName(); - } - } - else { - - while ((func = (RooAbsArg*)_funcIter->Next())) { - if (!first) { - os << " + "; - } - else { - first = kFALSE; - } - os << func->GetName(); - } - } - - os << " "; -} diff --git a/src/classes.h b/src/classes.h index 867ffeee1c9..ae1a054339a 100644 --- a/src/classes.h +++ b/src/classes.h @@ -3,7 +3,6 @@ #include "HiggsAnalysis/CombinedLimit/interface/VerticalInterpPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/VerticalInterpHistPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/AsymPow.h" -#include "HiggsAnalysis/CombinedLimit/interface/AsymQuad.h" #include "HiggsAnalysis/CombinedLimit/interface/CombDataSetFactory.h" #include "HiggsAnalysis/CombinedLimit/interface/TH1Keys.h" #include "HiggsAnalysis/CombinedLimit/interface/RooSimultaneousOpt.h" @@ -19,7 +18,6 @@ #include "HiggsAnalysis/CombinedLimit/interface/HZGRooPdfs.h" #include "HiggsAnalysis/CombinedLimit/interface/SequentialMinimizer.h" #include "HiggsAnalysis/CombinedLimit/interface/ProcessNormalization.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooRealFlooredSumPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/RooSpline1D.h" #include "HiggsAnalysis/CombinedLimit/interface/RooSplineND.h" #include "HiggsAnalysis/CombinedLimit/interface/RooScaleLOSM.h" @@ -56,13 +54,7 @@ #include "HiggsAnalysis/CombinedLimit/interface/RooTaylorExpansion.h" #include "HiggsAnalysis/CombinedLimit/interface/SimpleTaylorExpansion1D.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooPiecewisePolynomial.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooNCSplineCore.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooNCSpline_1D_fast.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooNCSpline_2D_fast.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooNCSpline_3D_fast.h" -#include "HiggsAnalysis/CombinedLimit/interface/RooFuncPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/RooCheapProduct.h" #include "HiggsAnalysis/CombinedLimit/interface/CMSHggFormula.h" #include "HiggsAnalysis/CombinedLimit/interface/SimpleProdPdf.h" diff --git a/src/classes_def.xml b/src/classes_def.xml index 3bf3eb4f689..6656870b2e7 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -36,7 +36,6 @@ - @@ -88,13 +87,6 @@ - - - - - - - @@ -154,7 +146,6 @@ -