diff --git a/interface/CachingNLL.h b/interface/CachingNLL.h index 35578f4fa31..43e89459253 100644 --- a/interface/CachingNLL.h +++ b/interface/CachingNLL.h @@ -21,6 +21,9 @@ #include class RooMultiPdf; +class SimNLLDerivativesHelper; +class DerivativeLogNormal; +class DerivativeRateParam; // Part zero: ArgSet checker namespace cacheutils { @@ -112,6 +115,9 @@ class OptimizedCachingPdfT : public CachingPdf { CachingPdfBase * makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ; class CachingAddNLL : public RooAbsReal { + friend SimNLLDerivativesHelper; // probably not needed w/ data + friend DerivativeLogNormal; + friend DerivativeRateParam; public: CachingAddNLL(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data, bool includeZeroWeights = false) ; CachingAddNLL(const CachingAddNLL &other, const char *name = 0) ; @@ -125,6 +131,7 @@ class CachingAddNLL : public RooAbsReal { virtual RooArgSet* getParameters(const RooArgSet* depList, Bool_t stripDisconnected = kTRUE) const ; double sumWeights() const { return sumWeights_; } const RooAbsPdf *pdf() const { return pdf_; } + const RooAbsData *data() const {return data_;} void setZeroPoint() ; void clearZeroPoint() ; void clearConstantZeroPoint() ; @@ -186,6 +193,8 @@ class CachingSimNLL : public RooAbsReal { // trap this call, since we don't care about propagating it to the sub-components virtual void constOptimizeTestStatistic(ConstOpCode opcode, Bool_t doAlsoTrackingOpt=kTRUE) { } private: + friend SimNLLDerivativesHelper; + void setup_(); RooSimultaneous *pdfOriginal_; const RooAbsData *dataOriginal_; diff --git a/interface/CascadeMinimizer.h b/interface/CascadeMinimizer.h index 89536f5108c..3d3c61f3cce 100644 --- a/interface/CascadeMinimizer.h +++ b/interface/CascadeMinimizer.h @@ -12,6 +12,9 @@ class RooRealVar; #include #include +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" +#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" +#include class CascadeMinimizer { public: @@ -27,10 +30,10 @@ class CascadeMinimizer { bool improve(int verbose=0, bool cascade=true, bool forceResetMinimizer=false); // declare nuisance parameters for pre-fit void setNuisanceParameters(const RooArgSet *nuis) { nuisances_ = nuis; } - RooMinimizer & minimizer() { return *minimizer_; } + RooMinimizer & minimizer() { if(isSemiAnalyticMinimizer) throw std::runtime_error("unimplemented"); return *minimizer_; } RooFitResult *save() { return minimizer().save(); } void setStrategy(int strategy) { strategy_ = strategy; } - void setErrorLevel(float errorLevel) { minimizer_->setErrorLevel(errorLevel); } + void setErrorLevel(float errorLevel) { (not isSemiAnalyticMinimizer)? minimizer_->setErrorLevel(errorLevel):minimizerSemiAnalytic_->setErrorLevel(errorLevel); } static void initOptions() ; static void applyOptions(const boost::program_options::variables_map &vm) ; static const boost::program_options::options_description & options() { return options_; } @@ -44,6 +47,12 @@ class CascadeMinimizer { private: RooAbsReal & nll_; std::auto_ptr minimizer_; + // + bool isSemiAnalyticMinimizer{false}; + std::auto_ptr minimizerSemiAnalytic_; + std::map derivatives_; // not used + std::unique_ptr helper_; + // Mode mode_; static int strategy_; RooRealVar * poi_; diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 0721e198f25..ca4cd7ab48f 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -14,6 +14,10 @@ ProcessNormalization is helper class for implementing process normalizations END_HTML */ // +class SimNLLDerivativesHelper; +class DerivativeLogNormal; +class DerivativeRateParam; + class ProcessNormalization : public RooAbsReal { public: ProcessNormalization() : nominalValue_(1) {} @@ -33,6 +37,9 @@ class ProcessNormalization : public RooAbsReal { Double_t evaluate() const; private: + friend SimNLLDerivativesHelper; + friend DerivativeLogNormal; + friend DerivativeRateParam; // ---- PERSISTENT ---- double nominalValue_; std::vector logKappa_; // Logarithm of symmetric kappas diff --git a/interface/RooFitResult.h b/interface/RooFitResult.h new file mode 100644 index 00000000000..47224da1ec0 --- /dev/null +++ b/interface/RooFitResult.h @@ -0,0 +1,207 @@ +/***************************************************************************** + * Project: RooFit * + * Package: RooFitCore * + * File: $Id: RooFitResult.h,v 1.28 2007/05/11 09:11:30 verkerke Exp $ + * Authors: * + * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu * + * DK, David Kirkby, UC Irvine, dkirkby@uci.edu * + * * + * Copyright (c) 2000-2005, Regents of the University of California * + * and Stanford University. All rights reserved. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ +#ifndef ROO_FIT_RESULT +#define ROO_FIT_RESULT + +#include "RooAbsArg.h" +#include "RooPrintable.h" +#include "RooDirItem.h" +#include "RooArgList.h" + +#include "RVersion.h" +#include "TMatrixFfwd.h" +#include "TMatrixDSym.h" +#include "TRootIOCtor.h" + +#include +#include +#include + +class RooArgSet ; +class RooAbsPdf ; +class RooPlot; +class TObject ; +class TH2 ; +typedef RooArgSet* pRooArgSet ; + +class RooFitResult : public TNamed, public RooPrintable, public RooDirItem { +public: + + // Constructors, assignment etc. + RooFitResult(const char* name=0, const char* title=0) ; + RooFitResult(const RooFitResult& other) ; + virtual TObject* Clone(const char* newname = 0) const { + RooFitResult* r = new RooFitResult(*this) ; + if (newname && *newname) r->SetName(newname) ; + return r ; + } + virtual TObject* clone() const { return new RooFitResult(*this); } + virtual ~RooFitResult() ; + + static RooFitResult* lastMinuitFit(const RooArgList& varList=RooArgList()) ; + + static RooFitResult *prefitResult(const RooArgList ¶mList); + + // Printing interface (human readable) + virtual void printValue(std::ostream& os) const ; + virtual void printName(std::ostream& os) const ; + virtual void printTitle(std::ostream& os) const ; + virtual void printClassName(std::ostream& os) const ; + virtual void printArgs(std::ostream& os) const ; + void printMultiline(std::ostream& os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const ; + + inline virtual void Print(Option_t *options= 0) const { + // Printing interface + printStream(defaultPrintStream(),defaultPrintContents(options),defaultPrintStyle(options)); + } + + virtual Int_t defaultPrintContents(Option_t* opt) const ; + virtual StyleOption defaultPrintStyle(Option_t* opt) const ; + + RooAbsPdf* createHessePdf(const RooArgSet& params) const ; + + // Accessors + inline Int_t status() const { + // Return MINUIT status code + return _status ; + } + + inline UInt_t numStatusHistory() const { return _statusHistory.size() ; } + Int_t statusCodeHistory(UInt_t icycle) const ; + const char* statusLabelHistory(UInt_t icycle) const ; + + inline Int_t covQual() const { + // Return MINUIT quality code of covariance matrix + return _covQual ; + } + inline Int_t numInvalidNLL() const { + // Return number of NLL evaluations with problems + return _numBadNLL ; + } + inline Double_t edm() const { + // Return estimated distance to minimum + return _edm ; + } + inline Double_t minNll() const { + // Return minimized -log(L) value + return _minNLL ; + } + inline const RooArgList& constPars() const { + // Return list of constant parameters + return *_constPars ; + } + inline const RooArgList& floatParsInit() const { + // Return list of floating parameters before fit + return *_initPars ; + } + inline const RooArgList& floatParsFinal() const { + // Return list of floarting parameters after fit + return *_finalPars ; + } + + TH2* correlationHist(const char* name = "correlation_matrix") const ; + + Double_t correlation(const RooAbsArg& par1, const RooAbsArg& par2) const { + // Return correlation between par1 and par2 + return correlation(par1.GetName(),par2.GetName()) ; + } + const RooArgList* correlation(const RooAbsArg& par) const { + // Return pointer to list of correlations of all parameters with par + return correlation(par.GetName()) ; + } + + Double_t correlation(const char* parname1, const char* parname2) const ; + const RooArgList* correlation(const char* parname) const ; + + + const TMatrixDSym& covarianceMatrix() const ; + const TMatrixDSym& correlationMatrix() const ; + TMatrixDSym reducedCovarianceMatrix(const RooArgList& params) const ; + TMatrixDSym conditionalCovarianceMatrix(const RooArgList& params) const ; + + + // Global correlation accessors + Double_t globalCorr(const RooAbsArg& par) { return globalCorr(par.GetName()) ; } + Double_t globalCorr(const char* parname) ; + const RooArgList* globalCorr() ; + + + // Add objects to a 2D plot + inline RooPlot *plotOn(RooPlot *frame, const RooAbsArg &par1, const RooAbsArg &par2, + const char *options= "ME") const { + // Plot error ellipse in par1 and par2 on frame + return plotOn(frame,par1.GetName(),par2.GetName(),options); + } + RooPlot *plotOn(RooPlot *plot, const char *parName1, const char *parName2, + const char *options= "ME") const; + + // Generate random perturbations of the final parameters using the covariance matrix + const RooArgList& randomizePars() const; + + Bool_t isIdentical(const RooFitResult& other, Double_t tol=5e-5, Double_t tolCorr=1e-4, Bool_t verbose=kTRUE) const ; + + void SetName(const char *name) ; + void SetNameTitle(const char *name, const char* title) ; + +protected: + + friend class RooMinuit ; + friend class RooMinimizer ; + friend class RooMinimizerSemiAnalytic; // :( why with friendship? + + void setCovarianceMatrix(TMatrixDSym& V) ; + void setConstParList(const RooArgList& list) ; + void setInitParList(const RooArgList& list) ; + void setFinalParList(const RooArgList& list) ; + inline void setMinNLL(Double_t val) { _minNLL = val ; } + inline void setEDM(Double_t val) { _edm = val ; } + inline void setStatus(Int_t val) { _status = val ; } + inline void setCovQual(Int_t val) { _covQual = val ; } + inline void setNumInvalidNLL(Int_t val) { _numBadNLL=val ; } + void fillCorrMatrix() ; + void fillCorrMatrix(const std::vector& globalCC, const TMatrixDSym& corrs, const TMatrixDSym& covs) ; + void fillLegacyCorrMatrix() const ; + void fillPrefitCorrMatrix(); + void setStatusHistory(std::vector >& hist) { _statusHistory = hist ; } + + Double_t correlation(Int_t row, Int_t col) const; + Double_t covariance(Int_t row, Int_t col) const; + + Int_t _status ; // MINUIT status code + Int_t _covQual ; // MINUIT quality code of covariance matrix + Int_t _numBadNLL ; // Number calls with bad (zero,negative) likelihood + Double_t _minNLL ; // NLL at minimum + Double_t _edm ; // Estimated distance to minimum + RooArgList* _constPars ; // List of constant parameters + RooArgList* _initPars ; // List of floating parameters with initial values + RooArgList* _finalPars ; // List of floating parameters with final values + + mutable RooArgList* _globalCorr ; //! List of global correlation coefficients + mutable TList _corrMatrix ; //! Correlation matrix (list of RooArgLists) + + mutable RooArgList *_randomPars; //! List of floating parameters with most recent random perturbation applied + mutable TMatrixF* _Lt; //! triangular matrix used for generate random perturbations + + TMatrixDSym* _CM ; // Correlation matrix + TMatrixDSym* _VM ; // Covariance matrix + TVectorD* _GC ; // Global correlation coefficients + + std::vector > _statusHistory ; // History of status codes + + ClassDef(RooFitResult,5) // Container class for fit result +}; + +#endif diff --git a/interface/RooMinimizerFcnSemiAnalytic.h b/interface/RooMinimizerFcnSemiAnalytic.h new file mode 100644 index 00000000000..c73fdb9e000 --- /dev/null +++ b/interface/RooMinimizerFcnSemiAnalytic.h @@ -0,0 +1,119 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on the equivalent roofit code/file. * + * Wed Oct 16 14:04:34 CEST 2019 + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +#ifndef __ROOFIT_NOROOMINIMIZER + +#ifndef ROO_MINIMIZER_FCN_SEMIANALYTIC +#define ROO_MINIMIZER_FCN_SEMIANALYTIC + +#include "Math/IFunction.h" +#include "Math/MultiNumGradFunction.h" +#include "Fit/ParameterSettings.h" +#include "Fit/FitResult.h" + +#include "TMatrixDSym.h" + +#include "RooAbsReal.h" +#include "RooArgList.h" + +#include +#include +#include + +class RooMinimizerSemiAnalytic; + +class RooMinimizerFcnSemiAnalytic : + //public ROOT::Math::MultiNumGradFunction { --> need to be constructed on a function. With the nominal Fcn? + public ROOT::Math::IMultiGradFunction { + //public ROOT::Math::IBaseFunctionMultiDim { + + public: + + RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic *context, std::map* knownDerivatives, bool verbose = false); + RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other); + virtual ~RooMinimizerFcnSemiAnalytic(); + + ROOT::Math::IMultiGradFunction* Clone() const override; + unsigned int NDim() const override{ return _nDim; } + + RooArgList* GetFloatParamList() { return _floatParamList; } + RooArgList* GetConstParamList() { return _constParamList; } + RooArgList* GetInitFloatParamList() { return _initFloatParamList; } + RooArgList* GetInitConstParamList() { return _initConstParamList; } + + void SetEvalErrorWall(Bool_t flag) { _doEvalErrorWall = flag ; } + void SetPrintEvalErrors(Int_t numEvalErrors) { _printEvalErrors = numEvalErrors ; } + Bool_t SetLogFile(const char* inLogfile); + std::ofstream* GetLogFile() { return _logfile; } + void SetVerbose(Bool_t flag=kTRUE) { _verbose = flag ; } + + Double_t& GetMaxFCN() { return _maxFCN; } + Int_t GetNumInvalidNLL() { return _numBadNLL; } + + Bool_t Synchronize(std::vector& parameters, + Bool_t optConst, Bool_t verbose); + void BackProp(const ROOT::Fit::FitResult &results); + void ApplyCovarianceMatrix(TMatrixDSym& V); + + Int_t evalCounter() const { return _evalCounter ; } + void zeroEvalCount() { _evalCounter = 0 ; } + + + private: + + Double_t GetPdfParamVal(Int_t index); + Double_t GetPdfParamErr(Int_t index); + void SetPdfParamErr(Int_t index, Double_t value); + void ClearPdfParamAsymErr(Int_t index); + void SetPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal); + + inline Bool_t SetPdfParamVal(const Int_t &index, const Double_t &value) const; + + + double DoEval(const double * x) const override; + double DoDerivative(const double * x,unsigned int icoord) const override; + // semi-analytic. The following function is called by DoDerivative if analytical derivative is not present. + virtual double DoNumericalDerivative(const double * x,int icoord) const; + void updateFloatVec() ; + +private: + + mutable Int_t _evalCounter ; + + RooAbsReal *_funct; + // for name -> Derivative + std::map *_knownDerivatives; // does not own them. Otherwise needs to design ad hoc copy constructors. + + RooMinimizerSemiAnalytic *_context; + + mutable double _maxFCN; + mutable int _numBadNLL; + mutable int _printEvalErrors; + Bool_t _doEvalErrorWall; + + int _nDim; + std::ofstream *_logfile; + bool _verbose; + + RooArgList* _floatParamList; + std::vector _floatParamVec ; + std::vector _derivParamVec ; // vector of derivative functions - algined with the one above + // + RooArgList* _constParamList; + RooArgList* _initFloatParamList; + RooArgList* _initConstParamList; + + int _useNumDerivatives{1};// 0 = ROOT (not-impl), 1 re-implementation of gsl 5 point, 2 re-implementation of gsl 3 point + +}; + +#endif +#endif diff --git a/interface/RooMinimizerSemiAnalytic.h b/interface/RooMinimizerSemiAnalytic.h new file mode 100644 index 00000000000..d27df39c5dd --- /dev/null +++ b/interface/RooMinimizerSemiAnalytic.h @@ -0,0 +1,141 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on the equivalent roofit code/file. * + * Wed Oct 16 14:04:34 CEST 2019 + * * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + + +/* ~Copy Of RooMINIMIZER, + * the difference is that is templated in RooMinimizerFcn. i + * An instance of RooMinimizerCopy is present as typedef for compatibility. + */ + +#ifndef __ROOFIT_NOROOMINIMIZER + +#ifndef ROO_MINIMIZER_SEMIANALYTIC +#define ROO_MINIMIZER_SEMIANALYTIC + +#include "TObject.h" +#include "TStopwatch.h" +#include +#include "TMatrixDSymfwd.h" + + +#include "Fit/Fitter.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerFcnSemiAnalytic.h" +#include "RooMinimizerFcn.h" + +class RooAbsReal ; +class RooFitResult ; +class RooArgList ; +class RooRealVar ; +class RooArgSet ; +class TH2F ; +class RooPlot ; + +class RooMinimizerSemiAnalytic : public TObject { +public: + + RooMinimizerSemiAnalytic(RooAbsReal& function,std::map* knownDerivatives) ; + virtual ~RooMinimizerSemiAnalytic() ; + + enum Strategy { Speed=0, Balance=1, Robustness=2 } ; + enum PrintLevel { None=-1, Reduced=0, Normal=1, ExtraForProblem=2, Maximum=3 } ; + void setStrategy(Int_t strat) ; + void setErrorLevel(Double_t level) ; + void setEps(Double_t eps) ; + void optimizeConst(Int_t flag) ; + void setEvalErrorWall(Bool_t flag) { fitterFcn()->SetEvalErrorWall(flag); } + void setOffsetting(Bool_t flag) ; + void setMaxIterations(Int_t n) ; + void setMaxFunctionCalls(Int_t n) ; + + RooFitResult* fit(const char* options) ; + + Int_t migrad() ; + Int_t hesse() ; + Int_t minos() ; + Int_t minos(const RooArgSet& minosParamList) ; + Int_t seek() ; + Int_t simplex() ; + Int_t improve() ; + + Int_t minimize(const char* type, const char* alg=0) ; + + RooFitResult* save(const char* name=0, const char* title=0) ; + RooPlot* contour(RooRealVar& var1, RooRealVar& var2, + Double_t n1=1, Double_t n2=2, Double_t n3=0, + Double_t n4=0, Double_t n5=0, Double_t n6=0, unsigned int npoints = 50) ; + + Int_t setPrintLevel(Int_t newLevel) ; + void setPrintEvalErrors(Int_t numEvalErrors) { fitterFcn()->SetPrintEvalErrors(numEvalErrors); } + void setVerbose(Bool_t flag=kTRUE) { _verbose = flag ; fitterFcn()->SetVerbose(flag); } + void setProfile(Bool_t flag=kTRUE) { _profile = flag ; } + Bool_t setLogFile(const char* logf=0) { return fitterFcn()->SetLogFile(logf); } + + void setMinimizerType(const char* type) ; + + static void cleanup() ; + static RooFitResult* lastMinuitFit(const RooArgList& varList=RooArgList()) ; + + void saveStatus(const char* label, Int_t status) { _statusHistory.push_back(std::pair(label,status)) ; } + + Int_t evalCounter() const { return fitterFcn()->evalCounter() ; } + void zeroEvalCount() { fitterFcn()->zeroEvalCount() ; } + + ROOT::Fit::Fitter* fitter() ; + const ROOT::Fit::Fitter* fitter() const ; + +protected: + + friend class RooAbsPdf ; + void applyCovarianceMatrix(TMatrixDSym& V) ; + + void profileStart() ; + void profileStop() ; + + inline Int_t getNPar() const { return fitterFcn()->NDim() ; } + inline std::ofstream* logfile() { return fitterFcn()->GetLogFile(); } + inline Double_t& maxFCN() { return fitterFcn()->GetMaxFCN() ; } + + const RooMinimizerFcnSemiAnalytic* fitterFcn() const { return ( fitter()->GetFCN() ? (dynamic_cast ( fitter()->GetFCN()) ): dynamic_cast(_fcn) ) ; } + RooMinimizerFcnSemiAnalytic* fitterFcn() { return ( fitter()->GetFCN() ? (dynamic_cast (fitter()->GetFCN())) : dynamic_cast(_fcn) ) ; } + +private: + + Int_t _printLevel ; + Int_t _status ; + Bool_t _optConst ; + Bool_t _profile ; + RooAbsReal* _func ; + + Bool_t _verbose ; + TStopwatch _timer ; + TStopwatch _cumulTimer ; + Bool_t _profileStart ; + + TMatrixDSym* _extV ; + + RooMinimizerFcnSemiAnalytic *_fcn; // Here it is the difference + std::string _minimizerType; + + static ROOT::Fit::Fitter *_theFitter ; + + std::vector > _statusHistory ; + + RooMinimizerSemiAnalytic(const RooMinimizerSemiAnalytic&) ; + +public: + ClassDef(RooMinimizerSemiAnalytic,0) // RooMinimizerSemiAnalytic this should be a title? +} ; + + +#endif + +#endif diff --git a/interface/SimNLLDerivativesHelper.h b/interface/SimNLLDerivativesHelper.h new file mode 100644 index 00000000000..ea6dd2dffb5 --- /dev/null +++ b/interface/SimNLLDerivativesHelper.h @@ -0,0 +1,129 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, CERN (CH), andrea.carlo.marini@cern.ch * + * Code based on roofit code. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +#ifndef SimNLLDerivHelper_H +#define SimNLLDerivHelper_H +#include +#include +#include +#include +#include +#include "RooAddition.h" +#include "RooDataSet.h" +#include "RooRealVar.h" +#include "HiggsAnalysis/CombinedLimit/interface/CachingNLL.h" +#include "HiggsAnalysis/CombinedLimit/interface/ProcessNormalization.h" + +/* This class takes in input the CachingSimNLL and returns the derivatives of what it knows + * + */ + +// TODO: Make a common DerivativeAbstract class. Marke DerivativeLogNormal Inherits from this class +class DerivativeAbstract: public RooAbsReal +{ + protected: // Full access to derived classes + const RooDataSet * data_{nullptr}; // not owned. + // RooRealProxy? + cacheutils::CachingAddNLL * pdf_{nullptr}; // CachingAddNLL + RooRealProxy pdfproxy_; + public: + + DerivativeAbstract(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data): + RooAbsReal(name, title), + data_(data), + pdf_(pdf), + pdfproxy_( ( std::string("proxy_")+name).c_str() ,"",pdf) {}; + ~DerivativeAbstract(){}; + virtual Bool_t isDerived() const { return kTRUE; } + const RooDataSet *data() const {return data_;} + cacheutils::CachingAddNLL *pdf() { return pdf_; } + bool verbose{true}; + + virtual Double_t evaluate() const =0 ; // cacheutils::ReminderSum + +}; + +// TODO: Inherit from DerivativeAbstract +class DerivativeLogNormal: public RooAbsReal +{ + const RooDataSet * data_{nullptr}; // not owned. + // RooRealProxy? + cacheutils::CachingAddNLL * pdf_{nullptr}; // CachingAddNLL + RooRealProxy pdfproxy_; + std::string thetaname_{""}; + + RooRealProxy theta_;// useful for numerical derivatives -> terms + + public: + DerivativeLogNormal(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& thetaname,int&found); + ~DerivativeLogNormal(){}; + virtual Bool_t isDerived() const { return kTRUE; } + virtual DerivativeLogNormal *clone(const char *name = 0) const ; + const RooDataSet *data() const {return data_;} + cacheutils::CachingAddNLL *pdf() { return pdf_; } + + virtual Double_t evaluate() const override ; // cacheutils::ReminderSum + + // name-> index association per process + std::vector kappa_pos_; + + bool verbose{true}; + + //bool maskConstraints{false}; + + // for debug purposes. Compute numerically the derivative corresponding to evaluate + //Double_t numericalDerivative() const; +}; + +// Implementation of the derivative for rate params and rate POIs +class DerivativeRateParam : public DerivativeAbstract +{ + std::string ratename_{""}; + RooRealProxy rate_;// useful for numerical derivatives -> terms + public: + DerivativeRateParam(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data,const std::string& thetaname,int&found); + ~DerivativeRateParam(){}; + virtual DerivativeRateParam *clone(const char *name = 0) const ; + + virtual Double_t evaluate() const override ; // cacheutils::ReminderSum + // name-> index association per process in otherFactorList + std::vector rate_pos_; + +}; + +//This class takes the nll and constructs the derivatives terms +class SimNLLDerivativesHelper +{ + private: + //RooRealProxy? + cacheutils::CachingSimNLL* nll_{nullptr}; //not owned, keep pointer + RooRealVar unmask_;//("mask_derivative_constraint_","",1.); + + std::set getServersVars(const RooAbsArg *node); + public: + SimNLLDerivativesHelper( cacheutils::CachingSimNLL * nll ): unmask_("mask_derivative_constraint_","",1.) {nll_=nll;} + ~SimNLLDerivativesHelper(){ + for (auto x : derivatives_) x.second->Delete(); + derivatives_.clear(); + }; + void init(); + bool verbose{true}; // debug + std::map derivatives_; + + void setMaskConstraint(int val=1){ + if (val ==0 or val ==1){ + unmask_.setVal( double(1-val) ); + } + else throw std::invalid_argument("[SimNLLDerivativesHelper]::setMask called with value != 0,1"); + } +}; + +#endif + diff --git a/interface/SimpleGaussianConstraint.h b/interface/SimpleGaussianConstraint.h index f1f2fd4c2dc..a8bf9985b25 100644 --- a/interface/SimpleGaussianConstraint.h +++ b/interface/SimpleGaussianConstraint.h @@ -5,7 +5,7 @@ class SimpleGaussianConstraint : public RooGaussian { public: - SimpleGaussianConstraint() {} ; + SimpleGaussianConstraint() {}; SimpleGaussianConstraint(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean, RooAbsReal& _sigma): RooGaussian(name,title,_x,_mean,_sigma) { init(); } @@ -17,6 +17,8 @@ class SimpleGaussianConstraint : public RooGaussian { inline virtual ~SimpleGaussianConstraint() { } const RooAbsReal & getX() const { return x.arg(); } + const RooAbsReal & getMean() const { return mean.arg(); } + const RooAbsReal & getSigma() const { return sigma.arg(); } double getLogValFast() const { if (_valueDirty) { @@ -24,7 +26,8 @@ class SimpleGaussianConstraint : public RooGaussian { //Double_t sig = sigma ; //return -0.5*arg*arg/(sig*sig); _value = scale_*arg*arg; - _valueDirty = false; + //_valueDirty = false; + _valueDirty = true; // something wrong with cache? } return _value; } diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index 23952be84cd..0440a9017b4 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -611,7 +611,7 @@ cacheutils::CachingAddNLL::evaluate() const boost::ptr_vector::iterator itp = pdfs_.begin();//, edp = pdfs_.end(); std::vector::const_iterator itw; //bgw = weights_.begin();//, edw = weights_.end(); std::vector::iterator its, bgs = partialSum_.begin(), eds = partialSum_.end(); - double sumCoeff = 0; + double sumCoeff = 0.; bool allBasicIntegralsOk = (basicIntegrals_ == 1); //std::cout << "Performing evaluation of " << GetName() << std::endl; for ( ; itc != edc; ++itp, ++itc ) { @@ -985,7 +985,6 @@ cacheutils::CachingSimNLL::setup_() factorizedPdf_.release(); simpdf = dynamic_cast(pdfclone); } - std::auto_ptr catClone((RooAbsCategoryLValue*) simpdf->indexCat().Clone()); pdfs_.resize(catClone->numBins(NULL), 0); @@ -1078,12 +1077,12 @@ cacheutils::CachingSimNLL::evaluate() const double logpdfval = (*it)->getLogValFast(); //std::cout << "pdf " << (*it)->GetName() << " = " << logpdfval << std::endl; ret2 += (logpdfval + *itz); + //std::cout << "[CachingSimNLL]: Simple Gaussian Constraint pdf " << (*it)->GetName() << " = " << logpdfval << " pars="<<(*it)->getMean().getVal()<<" "<<(*it)->getX().getVal()<<"/"<<(*it)->getSigma().getVal() <<" ZP "<<*itz<< std::endl; } /// ============= FAST POISSON CONSTRAINTS ========= itz = constrainZeroPointsFastPoisson_.begin(); for (std::vector::const_iterator it = constrainPdfsFastPoisson_.begin(), ed = constrainPdfsFastPoisson_.end(); it != ed; ++it, ++itz) { double logpdfval = (*it)->getLogValFast(); - //std::cout << "pdf " << (*it)->GetName() << " = " << logpdfval << std::endl; ret2 += (logpdfval + *itz); } } diff --git a/src/CascadeMinimizer.cc b/src/CascadeMinimizer.cc index fafd4c30b65..460dc7bd3dd 100644 --- a/src/CascadeMinimizer.cc +++ b/src/CascadeMinimizer.cc @@ -15,6 +15,14 @@ #include +#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" + +namespace local{ + template + class getter{ + }; +}; + boost::program_options::options_description CascadeMinimizer::options_("Cascade Minimizer options"); std::vector CascadeMinimizer::fallbacks_; bool CascadeMinimizer::preScan_; @@ -59,11 +67,31 @@ CascadeMinimizer::CascadeMinimizer(RooAbsReal &nll, Mode mode, RooRealVar *poi) void CascadeMinimizer::remakeMinimizer() { cacheutils::CachingSimNLL *simnll = dynamic_cast(&nll_); if (simnll) simnll->setHideRooCategories(true); + minimizer_.reset(); // avoid two copies in memory - minimizer_.reset(new RooMinimizer(nll_)); + minimizerSemiAnalytic_.reset(); + + if (runtimedef::get("MINIMIZER_SemiAnalytic")){ + isSemiAnalyticMinimizer=true; + if (simnll==nullptr) { + std::cout<<"[ERROR]::[CascadeMinimizer] I need simnll to have SemiAnalytic derivatives done by SimNLLDerivativesHelper"<init(); + //helper_->setMaskConstraint(1); // Why the numerical derivatives correspond (wrongly) to this? TODO + minimizerSemiAnalytic_.reset(new RooMinimizerSemiAnalytic(nll_, (runtimedef::get("MINIMIZER_SemiAnalytic_NOANALYTIC"))? &derivatives_: &helper_->derivatives_)); + //minimizerSemiAnalytic_->setVerbose(kTRUE); // DEBUG + } + else{ + isSemiAnalyticMinimizer=false; + minimizer_.reset(new RooMinimizer(nll_)); + } + if (simnll) simnll->setHideRooCategories(false); } + bool CascadeMinimizer::freezeDiscParams(const bool freeze) { static bool freezeDisassParams = runtimedef::get(std::string("MINIMIZER_freezeDisassociatedParams")); @@ -100,29 +128,29 @@ bool CascadeMinimizer::improve(int verbose, bool cascade, bool forceResetMinimiz simnllbb->setAnalyticBarlowBeeston(true); forceResetMinimizer = true; } - if (forceResetMinimizer || !minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-1); + if (forceResetMinimizer || !(minimizer_.get() or minimizerSemiAnalytic_.get())) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); strategy_ = ROOT::Math::MinimizerOptions::DefaultStrategy(); // re-configure - minimizer_->setStrategy(strategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(strategy_):minimizerSemiAnalytic_->setStrategy(strategy_); std::string nominalType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string nominalAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - minimizer_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); if (approxPreFitTolerance_ > 0) { double tol = std::max(approxPreFitTolerance_, 10. * nominalTol); do { if (verbose > 1) std::cout << "Running pre-fit with " << nominalType << "," << nominalAlgo << " and tolerance " << tol << std::endl; Significance::MinimizerSentry minimizerConfig(nominalType+","+nominalAlgo, tol); - minimizer_->setEps(tol); - minimizer_->setStrategy(approxPreFitStrategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(approxPreFitStrategy_):minimizerSemiAnalytic_->setStrategy(approxPreFitStrategy_); improveOnce(verbose-1, true); if (runtimedef::get("DBG_QUICKEXIT")) { exit(0); } - minimizer_->setEps(nominalTol); - minimizer_->setStrategy(strategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(approxPreFitStrategy_):minimizerSemiAnalytic_->setStrategy(strategy_); } while (autoBounds_ && !autoBoundsOk(verbose-1)); } bool outcome; @@ -145,8 +173,8 @@ bool CascadeMinimizer::improve(int verbose, bool cascade, bool forceResetMinimiz std::cerr << "Will fallback to minimization using " << it->algo << ", strategy " << myStrategy << " and tolerance " << it->tolerance << std::endl; Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Will fallback to minimization using %s, strategy %d and tolerance %g",__LINE__,(it->algo).c_str(),myStrategy,it->tolerance)),Logger::kLogLevelDebug,__func__); } - minimizer_->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()); - minimizer_->setStrategy(myStrategy); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()):minimizerSemiAnalytic_->setEps(ROOT::Math::MinimizerOptions::DefaultTolerance()); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(myStrategy):minimizerSemiAnalytic_->setStrategy(myStrategy); outcome = improveOnce(verbose-2); if (outcome) break; } @@ -171,37 +199,42 @@ bool CascadeMinimizer::improveOnce(int verbose, bool noHesse) bool outcome = false; double tol = ROOT::Math::MinimizerOptions::DefaultTolerance(); static int maxcalls = runtimedef::get("MINIMIZER_MaxCalls"); - if (!minimizer_.get()) remakeMinimizer(); + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); // freeze non active parameters if MINIMIZER_freezeDisassociatedParams enabled freezeDiscParams(true); if (maxcalls) { - minimizer_->setMaxFunctionCalls(maxcalls); - minimizer_->setMaxIterations(maxcalls); + (not isSemiAnalyticMinimizer)?minimizer_->setMaxFunctionCalls(maxcalls):minimizerSemiAnalytic_->setMaxFunctionCalls(maxcalls); + (not isSemiAnalyticMinimizer)?minimizer_->setMaxIterations(maxcalls):minimizerSemiAnalytic_->setMaxIterations(maxcalls); } if (oldFallback_){ - if (optConst) minimizer_->optimizeConst(std::max(0,optConst)); - if (rooFitOffset) minimizer_->setOffsetting(std::max(0,rooFitOffset)); + if (isSemiAnalyticMinimizer)throw std::logic_error("unimplemented");// check robust minimizer stuff with Minimizer + + if (optConst){ + minimizer_->optimizeConst(std::max(0,optConst)); + } + if (rooFitOffset) { minimizer_->setOffsetting(std::max(0,rooFitOffset));} + outcome = nllutils::robustMinimize(nll_, *minimizer_, verbose, setZeroPoint_); } else { if (verbose+2>0) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Minimisation configured with Type=%s, Algo=%s, strategy=%d, tolerance=%g",__LINE__,myType.c_str(),myAlgo.c_str(),myStrategy,tol)),Logger::kLogLevelInfo,__func__); cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast(&nll_) : 0; if (simnll) simnll->setZeroPoint(); - if ((!simnll) && optConst) minimizer_->optimizeConst(std::max(0,optConst)); - if ((!simnll) && rooFitOffset) minimizer_->setOffsetting(std::max(0,rooFitOffset)); + if ((!simnll) && optConst) (not isSemiAnalyticMinimizer)?minimizer_->optimizeConst(std::max(0,optConst)):minimizerSemiAnalytic_->optimizeConst(std::max(0,optConst)); + if ((!simnll) && rooFitOffset) (not isSemiAnalyticMinimizer)?minimizer_->setOffsetting(std::max(0,rooFitOffset)):minimizerSemiAnalytic_->setOffsetting(std::max(0,rooFitOffset)); if (firstHesse_ && !noHesse) { - minimizer_->setPrintLevel(std::max(0,verbose-3)); - minimizer_->hesse(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(std::max(0,verbose-3)):minimizerSemiAnalytic_->setPrintLevel(std::max(0,verbose-3)); + (not isSemiAnalyticMinimizer)?minimizer_->hesse():minimizerSemiAnalytic_->hesse(); if (simnll) simnll->updateZeroPoint(); - minimizer_->setPrintLevel(verbose-1); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); } - int status = minimizer_->minimize(myType.c_str(), myAlgo.c_str()); + int status = (not isSemiAnalyticMinimizer) ? minimizer_->minimize(myType.c_str(), myAlgo.c_str()):minimizerSemiAnalytic_->minimize(myType.c_str(), myAlgo.c_str()); if (lastHesse_ && !noHesse) { if (simnll) simnll->updateZeroPoint(); - minimizer_->setPrintLevel(std::max(0,verbose-3)); - status = minimizer_->hesse(); - minimizer_->setPrintLevel(verbose-1); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(std::max(0,verbose-3)):minimizerSemiAnalytic_->setPrintLevel(std::max(0,verbose-3)); + status = (not isSemiAnalyticMinimizer)?minimizer_->hesse():minimizerSemiAnalytic_->hesse(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); if (verbose+2>0 ) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Hesse finished with status=%d",__LINE__,status)),Logger::kLogLevelDebug,__func__); } if (simnll) simnll->clearZeroPoint(); @@ -239,8 +272,8 @@ bool CascadeMinimizer::minos(const RooArgSet & params , int verbose ) { utils::setAllConstant(toFreeze, false); remakeMinimizer(); } - if (!minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-1); // for debugging + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); // for debugging std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); @@ -255,8 +288,8 @@ bool CascadeMinimizer::minos(const RooArgSet & params , int verbose ) { // freeze parameters not active under current indexes if MINIMIZER_freezeDisassociatedParams enabled freezeDiscParams(true); // need to re-run Migrad before running minos - minimizer_->minimize(myType.c_str(), "Migrad"); - int iret = minimizer_->minos(params); + (not isSemiAnalyticMinimizer)?minimizer_->minimize(myType.c_str(), "Migrad"):minimizerSemiAnalytic_->minimize(myType.c_str(), "Migrad"); + int iret = (not isSemiAnalyticMinimizer)?minimizer_->minos(params):minimizerSemiAnalytic_->minos(params); if (verbose>0 ) Logger::instance().log(std::string(Form("CascadeMinimizer.cc: %d -- Minos finished with status=%d",__LINE__,iret)),Logger::kLogLevelDebug,__func__); freezeDiscParams(false); @@ -281,12 +314,12 @@ bool CascadeMinimizer::hesse(int verbose ) { // Have to reset and minimize again first to get all parameters in remakeMinimizer(); float nominalTol(ROOT::Math::MinimizerOptions::DefaultTolerance()); - minimizer_->setEps(nominalTol); - minimizer_->setStrategy(strategy_); + (not isSemiAnalyticMinimizer)?minimizer_->setEps(nominalTol):minimizerSemiAnalytic_->setEps(nominalTol); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(strategy_):minimizerSemiAnalytic_->setStrategy(strategy_); improveOnce(verbose - 1); } - if (!minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-1); // for debugging + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-1):minimizerSemiAnalytic_->setPrintLevel(verbose-1); // for debugging std::string myType(ROOT::Math::MinimizerOptions::DefaultMinimizerType()); std::string myAlgo(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo()); @@ -298,7 +331,7 @@ bool CascadeMinimizer::hesse(int verbose ) { } freezeDiscParams(true); - int iret = minimizer_->hesse(); + int iret = (not isSemiAnalyticMinimizer)?minimizer_->hesse():minimizerSemiAnalytic_->hesse(); freezeDiscParams(false); if (setZeroPoint_) { @@ -386,9 +419,9 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) bool doMultipleMini = (CascadeMinimizerGlobalConfigs::O().pdfCategories.getSize()>0); if (runtimedef::get(std::string("MINIMIZER_skipDiscreteIterations"))) doMultipleMini=false; // if ( doMultipleMini ) preFit_ = 1; - if (!minimizer_.get()) remakeMinimizer(); - minimizer_->setPrintLevel(verbose-2); - minimizer_->setStrategy(strategy_); + if (!minimizer_.get() and !minimizerSemiAnalytic_.get()) remakeMinimizer(); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-2):minimizerSemiAnalytic_->setPrintLevel(verbose-2); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(strategy_):minimizerSemiAnalytic_->setStrategy(strategy_); RooArgSet nuisances = CascadeMinimizerGlobalConfigs::O().nuisanceParameters; @@ -399,13 +432,13 @@ bool CascadeMinimizer::minimize(int verbose, bool cascade) freezeDiscParams(true); remakeMinimizer(); - minimizer_->setPrintLevel(verbose-2); - minimizer_->setStrategy(preFit_-1); + (not isSemiAnalyticMinimizer)?minimizer_->setPrintLevel(verbose-2):minimizerSemiAnalytic_->setPrintLevel(verbose-2); + (not isSemiAnalyticMinimizer)?minimizer_->setStrategy(preFit_-1):minimizerSemiAnalytic_->setStrategy(preFit_-1); cacheutils::CachingSimNLL *simnll = setZeroPoint_ ? dynamic_cast(&nll_) : 0; if (simnll) simnll->setZeroPoint(); - if (optConst) minimizer_->optimizeConst(std::max(0,optConst)); - if (rooFitOffset) minimizer_->setOffsetting(std::max(0,rooFitOffset)); - minimizer_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()); + if (optConst) (not isSemiAnalyticMinimizer)?minimizer_->optimizeConst(std::max(0,optConst)):minimizerSemiAnalytic_->optimizeConst(std::max(0,optConst)); + if (rooFitOffset) (not isSemiAnalyticMinimizer)?minimizer_->setOffsetting(std::max(0,rooFitOffset)):minimizerSemiAnalytic_->setOffsetting(std::max(0,rooFitOffset)); + (not isSemiAnalyticMinimizer)?minimizer_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()):minimizerSemiAnalytic_->minimize(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(), ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str()); if (simnll) simnll->clearZeroPoint(); utils::setAllConstant(frozen,false); freezeDiscParams(false); @@ -556,7 +589,8 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters, if (hideConstants && simnll) { simnll->setHideConstants(true); if (maskConstraints) simnll->setMaskConstraints(true); - minimizer_.reset(); // will be recreated when needed by whoever needs it + minimizer_.reset(); + minimizerSemiAnalytic_.reset(); // will be recreated when needed by whoever needs it } std::vector > myCombos; @@ -712,6 +746,7 @@ bool CascadeMinimizer::multipleMinimize(const RooArgSet &reallyCleanParameters, simnll->setHideConstants(false); if (maskConstraints) simnll->setMaskConstraints(false); minimizer_.reset(); // will be recreated when needed by whoever needs it + minimizerSemiAnalytic_.reset(); } diff --git a/src/RooMinimizerFcnSemiAnalytic.cxx b/src/RooMinimizerFcnSemiAnalytic.cxx new file mode 100644 index 00000000000..eb6c735ddff --- /dev/null +++ b/src/RooMinimizerFcnSemiAnalytic.cxx @@ -0,0 +1,777 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on roofit code. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +#ifndef __ROOFIT_NOROOMINIMIZER + +////////////////////////////////////////////////////////////////////////////// +// +// RooMinimizerFcnSemiAnalytic is an interface class to the ROOT::Math function +// for minization. +// + +#include + +#include "RooFit.h" +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerFcnSemiAnalytic.h" + +#include "Riostream.h" + +#include "TIterator.h" +#include "TClass.h" + +#include "RooAbsArg.h" +#include "RooAbsPdf.h" +#include "RooArgSet.h" +#include "RooRealVar.h" +#include "RooAbsRealLValue.h" +#include "RooMsgService.h" + +#include "RooMinimizer.h" +#include + +#include "RooAddition.h" // DEBUG & Print +#include "RooFormulaVar.h" //DEBUG & Print +//#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" // DEBUG to have evaluate +#define DEBUG_DIRTY + + +using namespace std; + +RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(RooAbsReal *funct, RooMinimizerSemiAnalytic* context, std::map * knownDerivatives, bool verbose) : + _funct(funct), + _context(context), + // Reset the *largest* negative log-likelihood value we have seen so far + _maxFCN(-1e30), _numBadNLL(0), + _printEvalErrors(10), _doEvalErrorWall(kTRUE), + _nDim(0), _logfile(0), + _verbose(verbose) +{ + + _knownDerivatives = knownDerivatives; + _evalCounter = 0 ; + + // Examine parameter list + RooArgSet* paramSet = _funct->getParameters(RooArgSet()); + RooArgList paramList(*paramSet); + delete paramSet; + + _floatParamList = (RooArgList*) paramList.selectByAttrib("Constant",kFALSE); + if (_floatParamList->getSize()>1) { + _floatParamList->sort(); + } + _floatParamList->setName("floatParamList"); + + _constParamList = (RooArgList*) paramList.selectByAttrib("Constant",kTRUE); + if (_constParamList->getSize()>1) { + _constParamList->sort(); + } + _constParamList->setName("constParamList"); + + // Remove all non-RooRealVar parameters from list (MINUIT cannot handle them) + TIterator* pIter = _floatParamList->createIterator(); + RooAbsArg* arg; + while ((arg=(RooAbsArg*)pIter->Next())) { + if (!arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic: removing parameter " + << arg->GetName() + << " from list because it is not of type RooRealVar" << endl; + _floatParamList->remove(*arg); + } + } + delete pIter; + + _nDim = _floatParamList->getSize(); + + _verbose=true; // print the variables for which we have derivatives? + updateFloatVec() ; + _verbose=verbose; + + // Save snapshot of initial lists + _initFloatParamList = (RooArgList*) _floatParamList->snapshot(kFALSE) ; + _initConstParamList = (RooArgList*) _constParamList->snapshot(kFALSE) ; + +} + + + +RooMinimizerFcnSemiAnalytic::RooMinimizerFcnSemiAnalytic(const RooMinimizerFcnSemiAnalytic& other) : ROOT::Math::IMultiGradFunction(other), + _evalCounter(other._evalCounter), + _funct(other._funct), + _context(other._context), + _maxFCN(other._maxFCN), + _numBadNLL(other._numBadNLL), + _printEvalErrors(other._printEvalErrors), + _doEvalErrorWall(other._doEvalErrorWall), + _nDim(other._nDim), + _logfile(other._logfile), + _verbose(other._verbose), + _floatParamVec(other._floatParamVec), + _derivParamVec(other._derivParamVec) +{ + _knownDerivatives = other._knownDerivatives; + _floatParamList = new RooArgList(*other._floatParamList) ; + _constParamList = new RooArgList(*other._constParamList) ; + _initFloatParamList = (RooArgList*) other._initFloatParamList->snapshot(kFALSE) ; + _initConstParamList = (RooArgList*) other._initConstParamList->snapshot(kFALSE) ; +} + + +RooMinimizerFcnSemiAnalytic::~RooMinimizerFcnSemiAnalytic() +{ + delete _floatParamList; + delete _initFloatParamList; + delete _constParamList; + delete _initConstParamList; +} + + +ROOT::Math::IMultiGradFunction* RooMinimizerFcnSemiAnalytic::Clone() const +{ + return new RooMinimizerFcnSemiAnalytic(*this) ; +} + + +Bool_t RooMinimizerFcnSemiAnalytic::Synchronize(std::vector& parameters, + Bool_t optConst, Bool_t verbose) +{ + + // Internal function to synchronize TMinimizer with current + // information in RooAbsReal function parameters + + Bool_t constValChange(kFALSE) ; + Bool_t constStatChange(kFALSE) ; + + Int_t index(0) ; + + // Handle eventual migrations from constParamList -> floatParamList + for(index= 0; index < _constParamList->getSize() ; index++) { + + RooRealVar *par= dynamic_cast(_constParamList->at(index)) ; + if (!par) continue ; + + RooRealVar *oldpar= dynamic_cast(_initConstParamList->at(index)) ; + if (!oldpar) continue ; + + // Test if constness changed + if (!par->isConstant()) { + + // Remove from constList, add to floatList + _constParamList->remove(*par) ; + _floatParamList->add(*par) ; + _initFloatParamList->addClone(*oldpar) ; + _initConstParamList->remove(*oldpar) ; + constStatChange=kTRUE ; + _nDim++ ; + + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: parameter " + << par->GetName() << " is now floating." << endl ; + } + } + + // Test if value changed + if (par->getVal()!= oldpar->getVal()) { + constValChange=kTRUE ; + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of constant parameter " + << par->GetName() + << " changed from " << oldpar->getVal() << " to " + << par->getVal() << endl ; + } + } + + } + + // Update reference list + *_initConstParamList = *_constParamList ; + + // Synchronize MINUIT with function state + // Handle floatParamList + for(index= 0; index < _floatParamList->getSize(); index++) { + RooRealVar *par= dynamic_cast(_floatParamList->at(index)) ; + + if (!par) continue ; + + Double_t pstep(0) ; + Double_t pmin(0) ; + Double_t pmax(0) ; + + if(!par->isConstant()) { + + // Verify that floating parameter is indeed of type RooRealVar + if (!par->IsA()->InheritsFrom(RooRealVar::Class())) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic::fit: Error, non-constant parameter " + << par->GetName() + << " is not of type RooRealVar, skipping" << endl ; + _floatParamList->remove(*par); + index--; + _nDim--; + continue ; + } + + // Set the limits, if not infinite + if (par->hasMin() ) + pmin = par->getMin(); + if (par->hasMax() ) + pmax = par->getMax(); + + // Calculate step size + pstep = par->getError(); + if(pstep <= 0) { + // Floating parameter without error estitimate + if (par->hasMin() && par->hasMax()) { + pstep= 0.1*(pmax-pmin); + + // Trim default choice of error if within 2 sigma of limit + if (pmax - par->getVal() < 2*pstep) { + pstep = (pmax - par->getVal())/2 ; + } else if (par->getVal() - pmin < 2*pstep) { + pstep = (par->getVal() - pmin )/2 ; + } + + // If trimming results in zero error, restore default + if (pstep==0) { + pstep= 0.1*(pmax-pmin); + } + + } else { + pstep=1 ; + } + if(verbose) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic::synchronize: WARNING: no initial error estimate available for " + << par->GetName() << ": using " << pstep << endl; + } + } + } else { + pmin = par->getVal() ; + pmax = par->getVal() ; + } + + // new parameter + if (index>=Int_t(parameters.size())) { + + if (par->hasMin() && par->hasMax()) { + parameters.push_back(ROOT::Fit::ParameterSettings(par->GetName(), + par->getVal(), + pstep, + pmin,pmax)); + } + else { + parameters.push_back(ROOT::Fit::ParameterSettings(par->GetName(), + par->getVal(), + pstep)); + if (par->hasMin() ) + parameters.back().SetLowerLimit(pmin); + else if (par->hasMax() ) + parameters.back().SetUpperLimit(pmax); + } + + continue; + + } + + Bool_t oldFixed = parameters[index].IsFixed(); + Double_t oldVar = parameters[index].Value(); + Double_t oldVerr = parameters[index].StepSize(); + Double_t oldVlo = parameters[index].LowerLimit(); + Double_t oldVhi = parameters[index].UpperLimit(); + + if (par->isConstant() && !oldFixed) { + + // Parameter changes floating -> constant : update only value if necessary + if (oldVar!=par->getVal()) { + parameters[index].SetValue(par->getVal()); + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of parameter " + << par->GetName() << " changed from " << oldVar + << " to " << par->getVal() << endl ; + } + } + parameters[index].Fix(); + constStatChange=kTRUE ; + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: parameter " + << par->GetName() << " is now fixed." << endl ; + } + + } else if (par->isConstant() && oldFixed) { + + // Parameter changes constant -> constant : update only value if necessary + if (oldVar!=par->getVal()) { + parameters[index].SetValue(par->getVal()); + constValChange=kTRUE ; + + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of fixed parameter " + << par->GetName() << " changed from " << oldVar + << " to " << par->getVal() << endl ; + } + } + + } else { + // Parameter changes constant -> floating + if (!par->isConstant() && oldFixed) { + parameters[index].Release(); + constStatChange=kTRUE ; + + if (verbose) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: parameter " + << par->GetName() << " is now floating." << endl ; + } + } + + // Parameter changes constant -> floating : update all if necessary + if (oldVar!=par->getVal() || oldVlo!=pmin || oldVhi != pmax || oldVerr!=pstep) { + parameters[index].SetValue(par->getVal()); + parameters[index].SetStepSize(pstep); + if (par->hasMin() && par->hasMax() ) + parameters[index].SetLimits(pmin,pmax); + else if (par->hasMin() ) + parameters[index].SetLowerLimit(pmin); + else if (par->hasMax() ) + parameters[index].SetUpperLimit(pmax); + } + + // Inform user about changes in verbose mode + if (verbose) { + // if ierr<0, par was moved from the const list and a message was already printed + + if (oldVar!=par->getVal()) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: value of parameter " + << par->GetName() << " changed from " << oldVar << " to " + << par->getVal() << endl ; + } + if (oldVlo!=pmin || oldVhi!=pmax) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: limits of parameter " + << par->GetName() << " changed from [" << oldVlo << "," << oldVhi + << "] to [" << pmin << "," << pmax << "]" << endl ; + } + + // If oldVerr=0, then parameter was previously fixed + if (oldVerr!=pstep && oldVerr!=0) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: error/step size of parameter " + << par->GetName() << " changed from " << oldVerr << " to " << pstep << endl ; + } + } + } + } + + if (optConst) { + if (constStatChange) { + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: set of constant parameters changed, rerunning const optimizer" << endl ; + _funct->constOptimizeTestStatistic(RooAbsArg::ConfigChange) ; + } else if (constValChange) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::synchronize: constant parameter values changed, rerunning const optimizer" << endl ; + _funct->constOptimizeTestStatistic(RooAbsArg::ValueChange) ; + } + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + + } + + updateFloatVec() ; + + return 0 ; + +} + +Double_t RooMinimizerFcnSemiAnalytic::GetPdfParamVal(Int_t index) +{ + // Access PDF parameter value by ordinal index (needed by MINUIT) + + return ((RooRealVar*)_floatParamList->at(index))->getVal() ; +} + +Double_t RooMinimizerFcnSemiAnalytic::GetPdfParamErr(Int_t index) +{ + // Access PDF parameter error by ordinal index (needed by MINUIT) + return ((RooRealVar*)_floatParamList->at(index))->getError() ; +} + + +void RooMinimizerFcnSemiAnalytic::SetPdfParamErr(Int_t index, Double_t value) +{ + // Modify PDF parameter error by ordinal index (needed by MINUIT) + + ((RooRealVar*)_floatParamList->at(index))->setError(value) ; +} + + + +void RooMinimizerFcnSemiAnalytic::ClearPdfParamAsymErr(Int_t index) +{ + // Modify PDF parameter error by ordinal index (needed by MINUIT) + + ((RooRealVar*)_floatParamList->at(index))->removeAsymError() ; +} + + +void RooMinimizerFcnSemiAnalytic::SetPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal) +{ + // Modify PDF parameter error by ordinal index (needed by MINUIT) + + ((RooRealVar*)_floatParamList->at(index))->setAsymError(loVal,hiVal) ; +} + + +void RooMinimizerFcnSemiAnalytic::BackProp(const ROOT::Fit::FitResult &results) +{ + // Transfer MINUIT fit results back into RooFit objects + + for (Int_t index= 0; index < _nDim; index++) { + Double_t value = results.Value(index); + SetPdfParamVal(index, value); + + // Set the parabolic error + Double_t err = results.Error(index); + SetPdfParamErr(index, err); + + Double_t eminus = results.LowerError(index); + Double_t eplus = results.UpperError(index); + + if(eplus > 0 || eminus < 0) { + // Store the asymmetric error, if it is available + SetPdfParamErr(index, eminus,eplus); + } else { + // Clear the asymmetric error + ClearPdfParamAsymErr(index) ; + } + } + +} + +Bool_t RooMinimizerFcnSemiAnalytic::SetLogFile(const char* inLogfile) +{ + // Change the file name for logging of a RooMinimizer of all MINUIT steppings + // through the parameter space. If inLogfile is null, the current log file + // is closed and logging is stopped. + + if (_logfile) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::setLogFile: closing previous log file" << endl ; + _logfile->close() ; + delete _logfile ; + _logfile = 0 ; + } + _logfile = new ofstream(inLogfile) ; + if (!_logfile->good()) { + /*oocoutI(_context,Minimization)*/std::cout << "RooMinimizerFcnSemiAnalytic::setLogFile: cannot open file " << inLogfile << endl ; + _logfile->close() ; + delete _logfile ; + _logfile= 0; + } + + return kFALSE ; + +} + + +void RooMinimizerFcnSemiAnalytic::ApplyCovarianceMatrix(TMatrixDSym& V) +{ + // Apply results of given external covariance matrix. i.e. propagate its errors + // to all RRV parameter representations and give this matrix instead of the + // HESSE matrix at the next save() call + + for (Int_t i=0 ; i<_nDim ; i++) { + // Skip fixed parameters + if (_floatParamList->at(i)->isConstant()) { + continue ; + } + SetPdfParamErr(i, sqrt(V(i,i))) ; + } + +} + + +Bool_t RooMinimizerFcnSemiAnalytic::SetPdfParamVal(const Int_t &index, const Double_t &value) const +{ + //RooRealVar* par = (RooRealVar*)_floatParamList->at(index); + RooRealVar* par = (RooRealVar*)_floatParamVec[index] ; + + if (par->getVal()!=value) { + if (_verbose) cout << par->GetName() << "=" << value << ", " ; + + par->setVal(value); + return kTRUE; + } + + return kFALSE; +} + + + +//////////////////////////////////////////////////////////////////////////////// + +void RooMinimizerFcnSemiAnalytic::updateFloatVec() +{ + // derivative vector and float vectors needs to be aligned + _floatParamVec.clear() ; + _derivParamVec.clear(); + + RooFIter iter = _floatParamList->fwdIterator() ; + RooAbsArg* arg ; + _floatParamVec = std::vector(_floatParamList->getSize()) ; + _derivParamVec = std::vector(_floatParamList->getSize()) ; + Int_t i(0) ; + while((arg=iter.next())) { + auto derivative=_knownDerivatives->find(arg->GetName()); + if (derivative == _knownDerivatives->end()){ + if(_verbose) std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "derivative for param "<GetName()<<" not found. I will use numerical derivatives."<GetName()<<" found with name "<second->GetName()<<"."<second; + } + _floatParamVec[i++] = arg ; + } +} + + + +double RooMinimizerFcnSemiAnalytic::DoEval(const double *x) const +{ + //std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "DoEval"<getVal(); + RooAbsReal::setHideOffset(kTRUE) ; + + //if (RooAbsReal::numEvalErrors()>0 || fvalue > 1e30) { REMOVEME + if (RooAbsPdf::evalError() || RooAbsReal::numEvalErrors()>0 || fvalue>1e30) { + + if (_printEvalErrors>=0) { + + if (_doEvalErrorWall) { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic: Minimized function has error status." << endl + << "Returning maximum FCN so far (" << _maxFCN + << ") to force MIGRAD to back out of this region. Error log follows" << endl ; + } else { + std::cout + /*oocoutW(_context,Minimization)*/ << "RooMinimizerFcnSemiAnalytic: Minimized function has error status but is ignored" << endl ; + } + + Bool_t first(kTRUE) ; + /*ooccoutW(_context,Minimization)*/std::cout << "Parameter values: " ; + //for (const auto par : *_floatParamList) { //new loop} + for (int ip =0 ;ip < _floatParamList->getSize();++ip){ + const auto par=_floatParamList->at(ip); + auto var = static_cast(par); + if (first) { first = kFALSE ; } else /*ooccoutW(_context,Minimization)*/std::cout << ", " ; + /*ooccoutW(_context,Minimization)*/std::cout << var->GetName() << "=" << var->getVal() ; + } + /*ooccoutW(_context,Minimization)*/std::cout << endl ; + + RooAbsReal::printEvalErrors(/*ooccoutW(_context,Minimization)*/std::cout,_printEvalErrors) ; + /*ooccoutW(_context,Minimization)*/std::cout << endl ; + } + + if (_doEvalErrorWall) { + fvalue = _maxFCN+1; + } + + RooAbsReal::clearEvalErrorLog() ; + _numBadNLL++ ; + } else { + _maxFCN = std::max(fvalue, _maxFCN); + } + + // Optional logging + if (_logfile) + (*_logfile) << setprecision(15) << fvalue << setprecision(4) << endl; + if (_verbose) { + cout << "\nprevFCN" << (_funct->isOffsetting()?"-offset":"") << " = " << setprecision(10) + << fvalue << setprecision(4) << " " ; + cout.flush() ; + } + + _evalCounter++ ; + + return fvalue; +} + +double RooMinimizerFcnSemiAnalytic::DoDerivative(const double * x, unsigned int icoord) const{ + //std::cout<<"[DEBUG]:["<<__PRETTY_FUNCTION__<<"]:"<< "DoDerivative:"<getVal()<getVal(); + RooAbsReal::setHideOffset(kTRUE) ; + return ret; + +} + +double RooMinimizerFcnSemiAnalytic::DoNumericalDerivative(const double * x,int icoord) const +{ + // ROOT + if (_useNumDerivatives==0) { + //return ROOT::Math::MultiNumGradFunction::DoDerivative(x,icoord); + //copied from ROOT::Math::MultiNumGradFunction::DoDerivative + + // calculate derivative using mathcore derivator class + // step size can be changes using SetDerivPrecision() + // this also use gsl. + + //static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + //double x0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + //double step = (x0>0) ? std::max( fgEps* x0, 8.0*kPrecision*(x0 + kPrecision) ) : kPrecision; + //return ROOT::Math::Derivator::Eval(*fFunc, x, icoord, step); + throw std::logic_error("Unimplemented"); + } + if (_useNumDerivatives==1){ + /* + * + ROOT use this has step above, maybe copy it dynamically? TODO + double MultiNumGradFunction::fgEps = 0.001; + + double MultiNumGradFunction::DoDerivative (const double * x, unsigned int icoord ) const + // calculate derivative using mathcore derivator class + // step size can be changes using SetDerivPrecision() + + static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + double x0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + double step = (x0>0) ? std::max( fgEps* x0, 8.0*kPrecision*(x0 + kPrecision) ) : kPrecision; + */ + + // from gsl_deriv_central -> central_deriv + /* Compute the derivative using the 5-point rule (x-h, x-h/2, x, + x+h/2, x+h). Note that the central point is not used. + + Compute the error using the difference between the 5-point and + the 3-point rule (x-h,x,x+h). Again the central point is not + used. */ + RooAbsReal::setHideOffset(kFALSE) ; + + // smallest number representable in doubles + static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + static double fgEps = 0.001; + + double ax0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + double h = (ax0>0) ? std::max( fgEps* ax0, 8.0*kPrecision*(ax0 + kPrecision) ) : kPrecision; + //static const double h=std::numeric_limits::epsilon()*8.; + const double hhalf=h/2.; + + //double f0 = _funct->getVal(); + //double x0 = x[icoord]; + + SetPdfParamVal(icoord,x[icoord]-h); + double fm1 = _funct->getVal(); + SetPdfParamVal(icoord,x[icoord]+h); + double fp1 = _funct->getVal(); + + SetPdfParamVal(icoord,x[icoord]-hhalf); + double fmh = _funct->getVal(); + SetPdfParamVal(icoord,x[icoord]+hhalf); + double fph = _funct->getVal(); + + RooAbsReal::setHideOffset(kTRUE) ; + + double r3 = 0.5 * (fp1 - fm1); + double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3; + + // don't do counts I'm not using. + //double e3 = (fabs (fp1) + fabs (fm1)) * GSL_DBL_EPSILON; + //double e5 = 2.0 * (fabs (fph) + fabs (fmh)) * GSL_DBL_EPSILON + e3; + + /* The next term is due to finite precision in x+h = O (eps * x) */ + + //double dy = std::max (fabs (r3 / h), fabs (r5 / h)) *(fabs (x) / h) * GSL_DBL_EPSILON; + + /* The truncation error in the r5 approximation itself is O(h^4). + However, for safety, we estimate the error from r5-r3, which is + O(h^2). By scaling h we will minimise this estimated error, not + the actual truncation error in r5. */ + + //*result = r5 / h; + //*abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */ + //*abserr_round = fabs (e5 / h) + dy; /* Rounding error (cancellations) */ + // reset + SetPdfParamVal(icoord,x[icoord]); + return r5 /h ; + } + + if (_useNumDerivatives==2){ + /* Compute the derivative using the 5-point rule (x-h, x-h/2, x, + x+h/2, x+h). Note that the central point is not used. + + Compute the error using the difference between the 5-point and + the 3-point rule (x-h,x,x+h). Again the central point is not + used. */ + + // smallest number representable in doubles + static double kPrecision = std::sqrt ( std::numeric_limits::epsilon() ); + static double fgEps = 0.001; + + double ax0 = std::abs(x[icoord]); + //double step = (x0 > 0) ? kPrecision * x0 : kPrecision; + // this seems to work better than above + double h = (ax0>0) ? std::max( fgEps* ax0, 8.0*kPrecision*(ax0 + kPrecision) ) : kPrecision; + //static const double h=std::numeric_limits::epsilon()*8.; + + SetPdfParamVal(icoord,x[icoord]-h); + double fm1 = _funct->getVal(); + SetPdfParamVal(icoord,x[icoord]+h); + double fp1 = _funct->getVal(); + + double r3 = 0.5 * (fp1 - fm1); + + // reset + SetPdfParamVal(icoord,x[icoord]); + return r3 /h ; + } + throw std::runtime_error(Form("Numerical derivatives method un-implemented: %d",_useNumDerivatives)); +} + + +#endif + diff --git a/src/RooMinimizerSemiAnalytic.cc b/src/RooMinimizerSemiAnalytic.cc new file mode 100644 index 00000000000..13e08fe3ba7 --- /dev/null +++ b/src/RooMinimizerSemiAnalytic.cc @@ -0,0 +1,942 @@ +/***************************************************************************** + * Authors: * + * Andrea Carlo Marini, MIT (MA), andrea.carlo.marini@cern.ch * + * Code based on the equivalent roofit code/file. * + * * + * Redistribution and use in source and binary forms, * + * with or without modification, are permitted according to the terms * + * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * + *****************************************************************************/ + +/** +\file RooMinimizerSemiAnalytic.cxx +\class RooMinimizerSemiAnalytic +\ingroup Roofitcore + +RooMinimizerSemiAnalytic is a copy of RooMinimizer. +RooMinimizerSemiAnalytic is templated in RooMinimizerFcn. +An instance is present as typedef definined on 'RooMinimizerFcn', that should act as an exact copy of RooMinimizer. +**/ + +#ifndef __ROOFIT_NOROOMINIMIZER + +#ifdef ROO_FIT_RESULT + #error I need to import RooFitResult +#else + #include "HiggsAnalysis/CombinedLimit/interface/RooFitResult.h" +#endif + +#include "RooFit.h" +#include "Riostream.h" + +#include "TClass.h" + +#include + +#include "TH1.h" +#include "TH2.h" +#include "TMarker.h" +#include "TGraph.h" +#include "Fit/FitConfig.h" +#include "TStopwatch.h" +#include "TDirectory.h" +#include "TMatrixDSym.h" + +#include "RooArgSet.h" +#include "RooArgList.h" +#include "RooAbsReal.h" +#include "RooAbsRealLValue.h" +#include "RooRealVar.h" +#include "RooAbsPdf.h" +#include "RooSentinel.h" +#include "RooMsgService.h" +#include "RooPlot.h" + + +#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h" +#include "RooFitResult.h" + +#include "Math/Minimizer.h" + +#if (__GNUC__==3&&__GNUC_MINOR__==2&&__GNUC_PATCHLEVEL__==3) +char* operator+( streampos&, char* ); +#endif + +using namespace std; + +//ClassImp(RooMinimizerSemiAnalytic); + + +ROOT::Fit::Fitter *RooMinimizerSemiAnalytic::_theFitter = 0 ; +//ROOT::Fit::Fitter *RooMinimizerSemiAnalytic::_theFitter = 0 ; + + + +//////////////////////////////////////////////////////////////////////////////// +/// Cleanup method called by atexit handler installed by RooSentinel +/// to delete all global heap objects when the program is terminated + +void RooMinimizerSemiAnalytic::cleanup() +{ + if (_theFitter) { + delete _theFitter ; + _theFitter =0 ; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Construct MINUIT interface to given function. Function can be anything, +/// but is typically a -log(likelihood) implemented by RooNLLVar or a chi^2 +/// (implemented by RooChi2Var). Other frequent use cases are a RooAddition +/// of a RooNLLVar plus a penalty or constraint term. This class propagates +/// all RooFit information (floating parameters, their values and errors) +/// to MINUIT before each MINUIT call and propagates all MINUIT information +/// back to the RooFit object at the end of each call (updated parameter +/// values, their (asymmetric errors) etc. The default MINUIT error level +/// for HESSE and MINOS error analysis is taken from the defaultErrorLevel() +/// value of the input function. + +RooMinimizerSemiAnalytic::RooMinimizerSemiAnalytic(RooAbsReal& function, std::map* knownDerivatives) +{ + RooSentinel::activate() ; + + // Store function reference + _extV = 0 ; + _func = &function ; + _optConst = kFALSE ; + _verbose = kFALSE ; + _profile = kFALSE ; + _profileStart = kFALSE ; + _printLevel = 1 ; + _minimizerType = "Minuit"; // default minimizer + + if (_theFitter) delete _theFitter ; + _theFitter = new ROOT::Fit::Fitter; + _fcn = new RooMinimizerFcnSemiAnalytic(_func,this,knownDerivatives,_verbose); + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + setEps(1.0); // default tolerance + // default max number of calls + _theFitter->Config().MinimizerOptions().SetMaxIterations(500*_fcn->NDim()); + _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(500*_fcn->NDim()); + + // Shut up for now + setPrintLevel(-1) ; + + // Use +0.5 for 1-sigma errors + setErrorLevel(_func->defaultErrorLevel()) ; + + // Declare our parameters to MINUIT + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + + // Now set default verbosity + if (RooMsgService::instance().silentMode()) { + setPrintLevel(-1) ; + } else { + setPrintLevel(1) ; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Destructor + +RooMinimizerSemiAnalytic::~RooMinimizerSemiAnalytic() +{ + if (_extV) { + delete _extV ; + } + + if (_fcn) { + delete _fcn; + } + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change MINUIT strategy to istrat. Accepted codes +/// are 0,1,2 and represent MINUIT strategies for dealing +/// most efficiently with fast FCNs (0), expensive FCNs (2) +/// and 'intermediate' FCNs (1) + +void RooMinimizerSemiAnalytic::setStrategy(Int_t istrat) +{ + _theFitter->Config().MinimizerOptions().SetStrategy(istrat); + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change maximum number of MINUIT iterations +/// (RooMinimizerSemiAnalytic default 500 * #parameters) + +void RooMinimizerSemiAnalytic::setMaxIterations(Int_t n) +{ + _theFitter->Config().MinimizerOptions().SetMaxIterations(n); +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change maximum number of likelihood function calss from MINUIT +/// (RooMinimizerSemiAnalytic default 500 * #parameters) + +void RooMinimizerSemiAnalytic::setMaxFunctionCalls(Int_t n) +{ + _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(n); +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Set the level for MINUIT error analysis to the given +/// value. This function overrides the default value +/// that is taken in the RooMinimizerSemiAnalytic constructor from +/// the defaultErrorLevel() method of the input function + +void RooMinimizerSemiAnalytic::setErrorLevel(Double_t level) +{ + _theFitter->Config().MinimizerOptions().SetErrorDef(level); + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change MINUIT epsilon + +void RooMinimizerSemiAnalytic::setEps(Double_t eps) +{ + _theFitter->Config().MinimizerOptions().SetTolerance(eps); + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Enable internal likelihood offsetting for enhanced numeric precision + +void RooMinimizerSemiAnalytic::setOffsetting(Bool_t flag) +{ + _func->enableOffsetting(flag) ; +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Choose the minimzer algorithm. + +void RooMinimizerSemiAnalytic::setMinimizerType(const char* type) +{ + _minimizerType = type; +} + + + + +//////////////////////////////////////////////////////////////////////////////// +/// Return underlying ROOT fitter object + +ROOT::Fit::Fitter* RooMinimizerSemiAnalytic::fitter() +{ + return _theFitter ; +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Return underlying ROOT fitter object + +const ROOT::Fit::Fitter* RooMinimizerSemiAnalytic::fitter() const +{ + return _theFitter ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Parse traditional RooAbsPdf::fitTo driver options +/// +/// m - Run Migrad only +/// h - Run Hesse to estimate errors +/// v - Verbose mode +/// l - Log parameters after each Minuit steps to file +/// t - Activate profile timer +/// r - Save fit result +/// 0 - Run Migrad with strategy 0 + +RooFitResult* RooMinimizerSemiAnalytic::fit(const char* options) +{ + TString opts(options) ; + opts.ToLower() ; + + // Initial configuration + if (opts.Contains("v")) setVerbose(1) ; + if (opts.Contains("t")) setProfile(1) ; + if (opts.Contains("l")) setLogFile(Form("%s.log",_func->GetName())) ; + if (opts.Contains("c")) optimizeConst(1) ; + + // Fitting steps + if (opts.Contains("0")) setStrategy(0) ; + migrad() ; + if (opts.Contains("0")) setStrategy(1) ; + if (opts.Contains("h")||!opts.Contains("m")) hesse() ; + if (!opts.Contains("m")) minos() ; + + return (opts.Contains("r")) ? save() : 0 ; +} + + + + +//////////////////////////////////////////////////////////////////////////////// + +Int_t RooMinimizerSemiAnalytic::minimize(const char* type, const char* alg) +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + + _theFitter->Config().SetMinimizer(type,alg); + + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MINIMIZE",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute MIGRAD. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::migrad() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migrad"); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MIGRAD",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute HESSE. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::hesse() +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::hesse: Error, run Migrad before Hesse!" + << endl ; + _status = -1; + } + else { + + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + bool ret = _theFitter->CalculateHessErrors(); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("HESSE",_status) ; + + } + + return _status ; + +} + +//////////////////////////////////////////////////////////////////////////////// +/// Execute MINOS. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::minos() +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::minos: Error, run Migrad before Minos!" + << endl ; + _status = -1; + } + else { + + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + bool ret = _theFitter->CalculateMinosErrors(); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MINOS",_status) ; + + } + + return _status ; + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute MINOS for given list of parameters. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::minos(const RooArgSet& minosParamList) +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::minos: Error, run Migrad before Minos!" + << endl ; + _status = -1; + } + else if (minosParamList.getSize()>0) { + + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + // get list of parameters for Minos + TIterator* aIter = minosParamList.createIterator() ; + RooAbsArg* arg ; + std::vector paramInd; + while((arg=(RooAbsArg*)aIter->Next())) { + RooAbsArg* par = _fcn->GetFloatParamList()->find(arg->GetName()); + if (par && !par->isConstant()) { + Int_t index = _fcn->GetFloatParamList()->index(par); + paramInd.push_back(index); + } + } + delete aIter ; + + if (paramInd.size()) { + // set the parameter indeces + _theFitter->Config().SetMinosErrors(paramInd); + + _theFitter->Config().SetMinimizer(_minimizerType.c_str()); + bool ret = _theFitter->CalculateMinosErrors(); + _status = ((ret) ? _theFitter->Result().Status() : -1); + // to avoid that following minimization computes automatically the Minos errors + _theFitter->Config().SetMinosErrors(kFALSE); + + } + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("MINOS",_status) ; + + } + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute SEEK. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::seek() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"seek"); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("SEEK",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute SIMPLEX. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::simplex() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"simplex"); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("SEEK",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Execute IMPROVE. Changes in parameter values +/// and calculated errors are automatically +/// propagated back the RooRealVars representing +/// the floating parameters in the MINUIT operation + +Int_t RooMinimizerSemiAnalytic::improve() +{ + _fcn->Synchronize(_theFitter->Config().ParamsSettings(), + _optConst,_verbose) ; + profileStart() ; + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + RooAbsReal::clearEvalErrorLog() ; + + _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migradimproved"); + bool ret = _theFitter->FitFCN(*(ROOT::Math::IMultiGradFunction*)_fcn); + _status = ((ret) ? _theFitter->Result().Status() : -1); + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + profileStop() ; + _fcn->BackProp(_theFitter->Result()); + + saveStatus("IMPROVE",_status) ; + + return _status ; +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Change the MINUIT internal printing level + +Int_t RooMinimizerSemiAnalytic::setPrintLevel(Int_t newLevel) +{ + Int_t ret = _printLevel ; + _theFitter->Config().MinimizerOptions().SetPrintLevel(newLevel+1); + _printLevel = newLevel+1 ; + return ret ; +} + +//////////////////////////////////////////////////////////////////////////////// +/// If flag is true, perform constant term optimization on +/// function being minimized. + +void RooMinimizerSemiAnalytic::optimizeConst(Int_t flag) +{ + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ; + + if (_optConst && !flag){ + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: deactivating const optimization" << endl ; + _func->constOptimizeTestStatistic(RooAbsArg::DeActivate) ; + _optConst = flag ; + } else if (!_optConst && flag) { + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: activating const optimization" << endl ; + _func->constOptimizeTestStatistic(RooAbsArg::Activate,flag>1) ; + _optConst = flag ; + } else if (_optConst && flag) { + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: const optimization already active" << endl ; + } else { + if (_printLevel>-1) coutI(Minimization) << "RooMinimizerSemiAnalytic::optimizeConst: const optimization wasn't active" << endl ; + } + + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ; + +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Save and return a RooFitResult snaphot of current minimizer status. +/// This snapshot contains the values of all constant parameters, +/// the value of all floating parameters at RooMinimizerSemiAnalytic construction and +/// after the last MINUIT operation, the MINUIT status, variance quality, +/// EDM setting, number of calls with evaluation problems, the minimized +/// function value and the full correlation matrix + +RooFitResult* RooMinimizerSemiAnalytic::save(const char* userName, const char* userTitle) +{ + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::save: Error, run minimization before!" + << endl ; + return 0; + } + + TString name,title ; + name = userName ? userName : Form("%s", _func->GetName()) ; + title = userTitle ? userTitle : Form("%s", _func->GetTitle()) ; + RooFitResult* fitRes = new RooFitResult(name,title) ; + + // Move eventual fixed parameters in floatList to constList + Int_t i ; + RooArgList saveConstList(*(_fcn->GetConstParamList())) ; + RooArgList saveFloatInitList(*(_fcn->GetInitFloatParamList())) ; + RooArgList saveFloatFinalList(*(_fcn->GetFloatParamList())) ; + for (i=0 ; i<_fcn->GetFloatParamList()->getSize() ; i++) { + RooAbsArg* par = _fcn->GetFloatParamList()->at(i) ; + if (par->isConstant()) { + saveFloatInitList.remove(*saveFloatInitList.find(par->GetName()),kTRUE) ; + saveFloatFinalList.remove(*par) ; + saveConstList.add(*par) ; + } + } + saveConstList.sort() ; + + fitRes->setConstParList(saveConstList) ; + fitRes->setInitParList(saveFloatInitList) ; + + fitRes->setStatus(_status) ; + fitRes->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ; + fitRes->setMinNLL(_theFitter->Result().MinFcnValue()) ; + fitRes->setNumInvalidNLL(_fcn->GetNumInvalidNLL()) ; + fitRes->setEDM(_theFitter->Result().Edm()) ; + fitRes->setFinalParList(saveFloatFinalList) ; + if (!_extV) { + std::vector globalCC; + TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ; + TMatrixDSym covs(_theFitter->Result().Parameters().size()) ; + for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) { + globalCC.push_back(_theFitter->Result().GlobalCC(ic)); + for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) { + corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii); + covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii); + } + } + fitRes->fillCorrMatrix(globalCC,corrs,covs) ; + } else { + fitRes->setCovarianceMatrix(*_extV) ; + } + + fitRes->setStatusHistory(_statusHistory) ; + + return fitRes ; + +} + +//////////////////////////////////////////////////////////////////////////////// +/// Create and draw a TH2 with the error contours in the parameters `var1` and `var2`. +/// \param[in] var1 The first parameter (x axis). +/// \param[in] var2 The second parameter (y axis). +/// \param[in] n1 First contour. +/// \param[in] n2 Optional contour. 0 means don't draw. +/// \param[in] n3 Optional contour. 0 means don't draw. +/// \param[in] n4 Optional contour. 0 means don't draw. +/// \param[in] n5 Optional contour. 0 means don't draw. +/// \param[in] n6 Optional contour. 0 means don't draw. +/// \param[in] npoints Number of points for evaluating the contour. +/// +/// Up to six contours can be drawn using the arguments `n1` to `n6` to request the desired +/// coverage in units of \f$ \sigma = n^2 \cdot \mathrm{ErrorDef} \f$. +/// See ROOT::Math::Minimizer::ErrorDef(). + +RooPlot* RooMinimizerSemiAnalytic::contour(RooRealVar& var1, RooRealVar& var2, + Double_t n1, Double_t n2, Double_t n3, + Double_t n4, Double_t n5, Double_t n6, unsigned int npoints) +{ + + + RooArgList* params = _fcn->GetFloatParamList() ; + RooArgList* paramSave = (RooArgList*) params->snapshot() ; + + // Verify that both variables are floating parameters of PDF + Int_t index1= _fcn->GetFloatParamList()->index(&var1); + if(index1 < 0) { + coutE(Minimization) << "RooMinimizerSemiAnalytic::contour(" << GetName() + << ") ERROR: " << var1.GetName() + << " is not a floating parameter of " + << _func->GetName() << endl ; + return 0; + } + + Int_t index2= _fcn->GetFloatParamList()->index(&var2); + if(index2 < 0) { + coutE(Minimization) << "RooMinimizerSemiAnalytic::contour(" << GetName() + << ") ERROR: " << var2.GetName() + << " is not a floating parameter of PDF " + << _func->GetName() << endl ; + return 0; + } + + // create and draw a frame + RooPlot* frame = new RooPlot(var1,var2) ; + + // draw a point at the current parameter values + TMarker *point= new TMarker(var1.getVal(), var2.getVal(), 8); + frame->addObject(point) ; + + // check first if a inimizer is available. If not means + // the minimization is not done , so do it + if (_theFitter->GetMinimizer()==0) { + coutW(Minimization) << "RooMinimizerSemiAnalytic::contour: Error, run Migrad before contours!" + << endl ; + return frame; + } + + + // remember our original value of ERRDEF + Double_t errdef= _theFitter->GetMinimizer()->ErrorDef(); + + Double_t n[6] ; + n[0] = n1 ; n[1] = n2 ; n[2] = n3 ; n[3] = n4 ; n[4] = n5 ; n[5] = n6 ; + + for (Int_t ic = 0 ; ic<6 ; ic++) { + if(n[ic] > 0) { + + // set the value corresponding to an n1-sigma contour + _theFitter->GetMinimizer()->SetErrorDef(n[ic]*n[ic]*errdef); + + // calculate and draw the contour + Double_t *xcoor = new Double_t[npoints+1]; + Double_t *ycoor = new Double_t[npoints+1]; + bool ret = _theFitter->GetMinimizer()->Contour(index1,index2,npoints,xcoor,ycoor); + + if (!ret) { + coutE(Minimization) << "RooMinimizerSemiAnalytic::contour(" + << GetName() + << ") ERROR: MINUIT did not return a contour graph for n=" + << n[ic] << endl ; + } else { + xcoor[npoints] = xcoor[0]; + ycoor[npoints] = ycoor[0]; + TGraph* graph = new TGraph(npoints+1,xcoor,ycoor); + + graph->SetName(Form("contour_%s_n%f",_func->GetName(),n[ic])) ; + graph->SetLineStyle(ic+1) ; + graph->SetLineWidth(2) ; + graph->SetLineColor(kBlue) ; + frame->addObject(graph,"L") ; + } + + delete [] xcoor; + delete [] ycoor; + } + } + + + // restore the original ERRDEF + _theFitter->Config().MinimizerOptions().SetErrorDef(errdef); + + // restore parameter values + *params = *paramSave ; + delete paramSave ; + + return frame ; + +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Start profiling timer + +void RooMinimizerSemiAnalytic::profileStart() +{ + if (_profile) { + _timer.Start() ; + _cumulTimer.Start(_profileStart?kFALSE:kTRUE) ; + _profileStart = kTRUE ; + } +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Stop profiling timer and report results of last session + +void RooMinimizerSemiAnalytic::profileStop() +{ + if (_profile) { + _timer.Stop() ; + _cumulTimer.Stop() ; + coutI(Minimization) << "Command timer: " ; _timer.Print() ; + coutI(Minimization) << "Session timer: " ; _cumulTimer.Print() ; + } +} + + + +//////////////////////////////////////////////////////////////////////////////// +/// Apply results of given external covariance matrix. i.e. propagate its errors +/// to all RRV parameter representations and give this matrix instead of the +/// HESSE matrix at the next save() call + +void RooMinimizerSemiAnalytic::applyCovarianceMatrix(TMatrixDSym& V) +{ + _extV = (TMatrixDSym*) V.Clone() ; + _fcn->ApplyCovarianceMatrix(*_extV); + +} + + +RooFitResult* RooMinimizerSemiAnalytic::lastMinuitFit(const RooArgList& varList) +{ + // Import the results of the last fit performed, interpreting + // the fit parameters as the given varList of parameters. + + if (_theFitter==0 || _theFitter->GetMinimizer()==0) { + oocoutE((TObject*)0,InputArguments) << "RooMinimizerSemiAnalytic::save: Error, run minimization before!" + << endl ; + return 0; + } + + // Verify length of supplied varList + if (varList.getSize()>0 && varList.getSize()!=Int_t(_theFitter->Result().NTotalParameters())) { + oocoutE((TObject*)0,InputArguments) + << "RooMinimizerSemiAnalytic::lastMinuitFit: ERROR: supplied variable list must be either empty " << endl + << " or match the number of variables of the last fit (" + << _theFitter->Result().NTotalParameters() << ")" << endl ; + return 0 ; + } + + + // Verify that all members of varList are of type RooRealVar + TIterator* iter = varList.createIterator() ; + RooAbsArg* arg ; + while((arg=(RooAbsArg*)iter->Next())) { + if (!dynamic_cast(arg)) { + oocoutE((TObject*)0,InputArguments) << "RooMinimizerSemiAnalytic::lastMinuitFit: ERROR: variable '" + << arg->GetName() << "' is not of type RooRealVar" << endl ; + return 0 ; + } + } + delete iter ; + + RooFitResult* res = new RooFitResult("lastMinuitFit","Last MINUIT fit") ; + + // Extract names of fit parameters + // and construct corresponding RooRealVars + RooArgList constPars("constPars") ; + RooArgList floatPars("floatPars") ; + + UInt_t i ; + for (i = 0; i < _theFitter->Result().NTotalParameters(); ++i) { + + TString varName(_theFitter->Result().GetParameterName(i)); + Bool_t isConst(_theFitter->Result().IsParameterFixed(i)) ; + + Double_t xlo = _theFitter->Config().ParSettings(i).LowerLimit(); + Double_t xhi = _theFitter->Config().ParSettings(i).UpperLimit(); + Double_t xerr = _theFitter->Result().Error(i); + Double_t xval = _theFitter->Result().Value(i); + + RooRealVar* var ; + if (varList.getSize()==0) { + + if ((xlosetConstant(isConst) ; + } else { + + var = (RooRealVar*) varList.at(i)->Clone() ; + var->setConstant(isConst) ; + var->setVal(xval) ; + if (xlosetRange(xlo,xhi) ; + } + + if (varName.CompareTo(var->GetName())) { + oocoutI((TObject*)0,Eval) << "RooMinimizerSemiAnalytic::lastMinuitFit: fit parameter '" << varName + << "' stored in variable '" << var->GetName() << "'" << endl ; + } + + } + + if (isConst) { + constPars.addOwned(*var) ; + } else { + var->setError(xerr) ; + floatPars.addOwned(*var) ; + } + } + + res->setConstParList(constPars) ; + res->setInitParList(floatPars) ; + res->setFinalParList(floatPars) ; + res->setMinNLL(_theFitter->Result().MinFcnValue()) ; + res->setEDM(_theFitter->Result().Edm()) ; + res->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ; + res->setStatus(_theFitter->Result().Status()) ; + std::vector globalCC; + TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ; + TMatrixDSym covs(_theFitter->Result().Parameters().size()) ; + for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) { + globalCC.push_back(_theFitter->Result().GlobalCC(ic)); + for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) { + corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii); + covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii); + } + } + res->fillCorrMatrix(globalCC,corrs,covs) ; + + return res; + +} + +#endif diff --git a/src/SimNLLDerivativesHelper.cc b/src/SimNLLDerivativesHelper.cc new file mode 100644 index 00000000000..23877773033 --- /dev/null +++ b/src/SimNLLDerivativesHelper.cc @@ -0,0 +1,679 @@ +#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h" +#include "RooDataHist.h" + +//#define DERIVATIVE_RATEPARAM_DEBUG 1 +//#define DERIVATIVE_LOGNORMAL_DEBUG 1 + +DerivativeLogNormal::DerivativeLogNormal(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& thetaname,int&found) : + RooAbsReal(name, title), + data_(data), + pdf_(pdf), + pdfproxy_( ( std::string("proxy_")+name).c_str() ,"",pdf), + thetaname_(thetaname) +{ + //setDirtyInhibit(true); // TODO: understand cache and proxies + setOperMode(RooAbsArg::ADirty,false); // Recursive, bool + + found=0; // keep track if at least one depends on the given kappa + for (unsigned int i=0;ipdfs_.size();++i) // process loop + { + if(verbose) std::cout<<"[DerivativeLogNormal]: loop "<< i <<"/"<pdfs_.size()<<" considering pdf "<GetName() < (pdf_->coeffs_[i]); + if (verbose) std::cout<<"[DerivativeLogNormal]: coeff is of class "<coeffs_[i]->ClassName()<(pdf_->coeffs_[i]); + if (p == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) p = dynamic_cast( &pp->components()[ip]); + } + } + // + if (p != nullptr) { // unable to add the sum for this process + + for (unsigned int j=0 ; jlogKappa_.size();++j) + { + if ( p->thetaList_.at(j)->GetName() == thetaname_ ) {val=int(j);} + } + if (val==-1) { + if(verbose) std::cout<<"[DerivativeLogNormal]: Unable to find kappa corresponding to "<GetName() <GetName() << "because not ProcessNormalization" <GetName() < sum of expectations in the bin + * lambdat_b -> sum of ( expecations times logK) + */ + + double sum=0.0; + // caching: TODO +#ifdef DERIVATIVE_LOGNORMAL_DEBUG + if(verbose) { + std::cout<<"[DerivativeLogNormal]: data"<Print("V"); + std::cout<<"[DerivativeLogNormal]: --"< (data_); + //if (dh==nullptr) std::cout<<"[DerivativeLogNormal]:"<<" NO ROODATAHIST"< reminder(pdf_->pdfs_.size(),0.0); // the sums must do 1, the ib==ndata does that with data=0 + int ndata=data_->numEntries(); + + for (int ib=0;ib<=ndata;++ib) + //for (unsigned int ib=0;ibweights_.size();++ib) + { + auto x= (ibget(ib): data_->get(0) ; + double db = (ibweight() : 0.; + //double bw = (dh != nullptr) ? dh->binVolume(): 1.0; + double bw = 1; //(pdf_->binWidths_.size() > 1) ? pdf_->binWidths_[ib]: (pdf_->binWidths_.size()==1) ?pdf_->binWidths_[0]: 1.0; + // const RooArgSet *obs = data_->get(); + // RooRealVar *xvar = dynamic_cast(obs->first()); + if (x->getSize()==1) { + RooRealVar *xvar=dynamic_cast(x->first()); + const RooAbsBinning &bins = xvar->getBinning(); + bw=bins.binWidth(0); // only costant supported. Need to figure out for zero bins + //if (verbose){ + // std::cout <<"BinWidth: "; + // for(int ii =0 ;iinumBins();++ii){ + // double bc2 = bins.binCenter(ii); + // std::cout<< bins.binWidth(ii); + // } + // std::cout <weights_[ib]; + //double bw = pdf_->binWidths_[ib]; + + //double lambda=pdf_->getVal()*bw; // sum coeff*pdf + double lambdat=0.0; + double lambda=0.;// debug: recompute lambda with loop. not sure that the term before include everything I want + + //bool recursive=true; + //double running_prod=1.0; + for (unsigned int i=0;ipdfs_.size();++i) // loop over processes + { + // find kappa + ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); + + // fix RooProduct + RooProduct *pp = dynamic_cast(pdf_->coeffs_[i]); + if (c == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); + } + } + //double integral = (pdf_->isRooRealSum_ && pdf_->basicIntegrals_ < 2) ? pdf_->integrals_[i]->getVal(x): 1.0; + //--- + double logK = (kappa_pos_[i]>=0)? c -> logKappa_ [kappa_pos_[i]] : 0.0; + if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" pdf is a RooRealSumPdf"<pdfs_.at(i).pdf()->setValueDirty(); + + // TH1 + if (ibpdfs_.at(i).pdf()->getVal(x)*bw; + double expectedEvents= (ibcoeffs_[i]->getVal(x) * pdf_->pdfs_.at(i).pdf()->getVal(x)*bw : + pdf_->coeffs_[i]->getVal(x) * (1.-reminder[i]); + + // sum coeff or sum coeff*integral + //double expectedEvents= dynamic_cast(pdf_->pdfs_.at(i).pdf())->expectedEvents(x); + //bool expEventsNoNorm=false; + //double expectedEvents = (pdf_->isRooRealSum_ && !expEventsNoNorm ? dynamic_cast(pdf_->pdfs_.at(i).pdf())->getNorm(data_->get()) : pdf_->coeffs_[i]->getVal()); + lambdat+= expectedEvents * logK; + lambda += expectedEvents; + + if (pdf_->basicIntegrals_ and verbose)std::cout<<"[DerivativeLogNormal]:DEBUG "<<" PDF has basic integrals:"<basicIntegrals_< 1*/) std::cout<<"[DerivativeLogNormal]: *** ibin="<< ib<< " data="<(pdf_->pdfs_.at(i).pdf())->getValV(x) << "pdfU="<(pdf_->pdfs_.at(i).pdf())->getValV(nullptr) <<"bw="<getVal()*bw <<")" <<" lambdat="<Delete(); + derivatives_.clear(); + //--- + + std::set logNormal; + std::set rateParams; // notice that logNormal cannot be in rateParams and vice-versa. Mixing could in principle be possible, but needs to be handled in logNormal to understand constraint part. + + // Find LogNormal candidates + for (auto f : nll_->constrainPdfsFast_) + { + if (verbose) std::cout<<"inserting candidate: "<< f->getX().GetName()<getX().GetName() ); + } + + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][init] List of constraint pdfs parameters:"<params_); + for (int ir=0;ir::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { + cacheutils::CachingAddNLL * pdf = *it; + if (pdf==nullptr) continue ; // ?!? needed? + const RooAbsData* data= pdf->data(); + bool isWeighted = data->isWeighted() ;//and (dynamic_cast(data) != nullptr); // binned vs unbinned?!? TBC + + if (verbose) { std::cout << "[SimNLLDerivativesHelper][init] > considering pdf "<GetName()<< ((isWeighted)?" isWeighted":" notWeighted") < toRemove; + + for(unsigned proc =0 ; proc < pdf ->pdfs_.size();++proc) + { + RooAbsReal * coeff = pdf -> coeffs_[proc] ; // what to do for RooProduct? + const RooAbsReal* pdfi= pdf->pdfs_[proc].pdf(); + + ProcessNormalization * pn = dynamic_cast (coeff); + RooProduct *pp = dynamic_cast(coeff); + RooRealVar *rv = dynamic_cast(coeff); + //if (rv !=nullptr) std::cout<<"[SimNLLDerivativesHelper]::[FIXME]"<<"derivative for "GetName()<<" is known, but discarded in the loop :( "<components().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) pn = dynamic_cast( &pp->components()[ip]); + else { + if (verbose) { std::cout <<"components:";} + std::set servers= getServersVars(&pp->components()[ip]); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|(LN)"< likely shape + if (verbose) { std::cout <<"removing pdf variables (shape):";} + std::set servers= getServersVars(pdfi); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|(LN)"<ClassName()): "") + << ((pn==nullptr and not isWeighted) ? " and":"" ) + << ((not isWeighted)?" dataset is not weighted (unbinned?)":"") + <Print("V"); + std::cout<<"[SimNLLDerivativesHelper][init]"<< " -- PDF -- "<Print("V"); + std::cout<<"-----------------"< servers= getServersVars(coeff); + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|(LN)"<*/ servers= getServersVars(pdfi); // do I still need this block? + for(auto v : servers) { + if (logNormal.find(v) != logNormal.end()) { + logNormal.erase(v); + if (verbose) std::cout<<"|(LN)"<> Removing all asymmThetaList "<asymmThetaList_.getSize() ;++idx){ + RooAbsArg*v=pn->asymmThetaList_.at(idx); + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName()<<" from the list of known derivatives because appears in an asymmThetaList " + <GetName()) != logNormal.end()) logNormal.erase(v->GetName()); + if (rateParams.find(v->GetName()) != rateParams.end()) {rateParams.erase(v->GetName());} // probably not necessary, since this are constrained + } + + // remove all otherFactorList + std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Removing all otherfactorsList "<otherFactorList_.getSize() ;++idx){ + RooAbsArg*v=pn->otherFactorList_.at(idx); + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName()<<" from the list of known derivatives (LN) because appears in an otherFactors " + <GetName()) != logNormal.end()) logNormal.erase(v->GetName()); + if (dynamic_cast(v)==nullptr and rateParams.find(v->GetName()) != rateParams.end()) { + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName()<< " : "<ClassName()<<" from the list of known derivatives (R) because appears in an otherFactors " + <GetName()); + } + } + + // remove all symmThetaList + for(int idx=0; idx< pn->thetaList_.getSize() ;++idx){ + RooAbsArg*v=pn->thetaList_.at(idx); + if (verbose) { + std::cout<<"[SimNLLDerivativesHelper][INFO]:" + <<"Removing "<GetName() + // << ","<GetX().GetName() + // << ","<GetMean().GetName() + // << ","<GetSigma().GetName() + <<" from the list of known derivatives because appears in an thetaList " + <GetName()) != rateParams.end()) {rateParams.erase(v->GetName());} // probably not necessary, since this are constrained + //if (rateParams.find(v->GetX().GetName()) != rateParams.end()) {rateParams.erase(v->GetX().GetName());} // probably not necessary, since this are constrained + //if (rateParams.find(v->GetMean().GetName()) != rateParams.end()) {rateParams.erase(v->GetMean().GetName());} // probably not necessary, since this are constrained + //if (rateParams.find(v->GetSigma().GetName()) != rateParams.end()) {rateParams.erase(v->GetSigma().GetName());} // probably not necessary, since this are constrained + } + std::cout<<"[SimNLLDerivativesHelper][init]"<< ">> Continue loop "< Continue loop "<::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { + // construct the derivative as sum + cacheutils::CachingAddNLL * pdf = *it; + if (pdf==nullptr) continue ; // ?!? needed? + const RooAbsData* data= pdf->data(); + //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string thetaname) + int found; + DerivativeLogNormal *der= new DerivativeLogNormal((std::string("dln_")+name+Form("_%d",idx)).c_str(),"",pdf,dynamic_cast(data),name,found); + if(found) list.add(*der); // or addOwned or add? + else der->Delete(); + } + + if (list.getSize() == 0){ // old fix for shapes. Leave it, probably useful for debug + std::cout<<"[SimNLLDerivativesHelper][init][ERROR]"<< "Something went wrong: Unable to match "<< name <<" with kappas"<constrainPdfsFast_ + std::cout<<"[SimNLLDerivativesHelper][init]"<< "Adding constraint term for "<< name <constrainPdfsFast_.size();++i) + { + if(verbose) std::cout<<"[SimNLLDerivativesHelper][init]"<< " Considering "<< nll_->constrainPdfsFast_[i]->getX().GetName() <constrainPdfsFast_[i]->getX().GetName() != name) continue; + + if(verbose) std::cout<<"[SimNLLDerivativesHelper][init]"<< "Loop. found constraint "<< nll_->constrainPdfsFast_[i]->getX().GetName() <<" == "<constrainPdfsFast_[i]->getX(),nll_->constrainPdfsFast_[i]->getMean(), nll_->constrainPdfsFast_[i]->getSigma(),unmask_));// (x-mean)/sigma^2 + list.add(*constr) ; // or addOwned or add? + } + if (verbose) {std::cout <<"[SimNLLDerivativesHelper][init] --- LIST OF addition ---"<::const_iterator it = nll_->pdfs_.begin(), ed = nll_->pdfs_.end(); it != ed; ++it, ++idx) { + // construct the derivative as sum + cacheutils::CachingAddNLL * pdf = *it; + if (pdf==nullptr) continue ; // ?!? needed? + const RooAbsData* data= pdf->data(); + //(const char *name, const char *title, RooAbsPdf *pdf, RooAbsData *data,std::string thetaname) + int found; + DerivativeRateParam *der= new DerivativeRateParam((std::string("dln_")+name+Form("_%d",idx)).c_str(),"",pdf,dynamic_cast(data),name,found); + if(found) list.add(*der); // or addOwned or add? + else der->Delete(); + } + + if (list.getSize() == 0){ // old fix for shapes. Leave it, probably useful for debug + std::cout<<"[SimNLLDerivativesHelper][init][ERROR]"<< "Something went wrong: Unable to match "<< name <<" with kappas"<constrainPdfsFast_ + if (verbose) {std::cout <<"[SimNLLDerivativesHelper][init] --- LIST OF addition ---"<first <<" : "<second->getVal(); // DEBUG + } +} + +//def getServers(node): +// servers = [] +// iter = node.serverIterator() +// while True: +// server = iter.Next() +// if server == None: +// break +// servers.append(server) +// return servers + +std::set SimNLLDerivativesHelper::getServersVars(const RooAbsArg *node){ + std::set R; + auto iter = node->serverIterator(); + while (true) + { + auto server = iter->Next(); + if (server==nullptr)break; + + if (dynamic_cast (server) != nullptr) + { + R.insert(std::string(server->GetName())); + } + else{ + std::set r1 = getServersVars((RooAbsArg*)server); + R.merge(r1); + } + } + return R; +} + + + +/// --- +DerivativeRateParam::DerivativeRateParam(const char *name, const char *title, cacheutils::CachingAddNLL *pdf, const RooDataSet *data, const std::string& ratename,int&found) : + DerivativeAbstract(name,title,pdf,data), + ratename_(ratename) +{ + //setDirtyInhibit(true); // TODO: understand cache and proxies + setOperMode(RooAbsArg::ADirty,false); // Recursive, bool + + found=0; // keep track if at least one depends on the given kappa + for (unsigned int i=0;ipdfs_.size();++i) // process loop + { + if(verbose) std::cout<<"[DerivativeRateParam]: loop "<< i <<"/"<pdfs_.size()<<" considering pdf "<GetName() < (pdf_->coeffs_[i]); + if (verbose) std::cout<<"[DerivativeRateParam]: coeff is of class "<coeffs_[i]->ClassName()<(pdf_->coeffs_[i]); + if (p == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) p = dynamic_cast( &pp->components()[ip]); + } + } + // + RooRealVar *rv = dynamic_cast (pdf_->coeffs_[i]); + if (rv !=nullptr){ + val=-2; // the whole coeff is a RooRealVar + found=1; + if(verbose) std::cout<<"[DerivativeRateParam]: All Coeff corresponding to "<GetName() <<" is a RooRealVar"<otherFactorList_.getSize();++j) + { + if ( p->otherFactorList_.at(j)->GetName() == ratename_ ) {val=int(j);} + } + if (val==-1) { + if(verbose) std::cout<<"[DerivativeRateParam]: Unable to find rate corresponding to "<GetName() <GetName() << "because not ProcessNormalization" <GetName() < sum of expectations in the bin + * lambdat_b -> sum of ( expectations with delta(rate param)) + */ + + double sum=0.0; + // caching: TODO +#ifdef DERIVATIVE_RATEPARAM_DEBUG + if(verbose) { + std::cout<<"[DerivativeRateParam]: data"<Print("V"); + std::cout<<"[DerivativeRateParam]: --"< (data_); + //if (dh==nullptr) std::cout<<"[DerivativeRateParam]:"<<" NO ROODATAHIST"< reminder(pdf_->pdfs_.size(),0.0); // the sums must do 1, the ib==ndata does that with data=0 + int ndata=data_->numEntries(); + + for (int ib=0;ib<=ndata;++ib) + //for (unsigned int ib=0;ibweights_.size();++ib) + { + auto x= (ibget(ib): data_->get(0) ; + double db = (ibweight() : 0.; + double bw = 1; + if (x->getSize()==1) { + RooRealVar *xvar=dynamic_cast(x->first()); + const RooAbsBinning &bins = xvar->getBinning(); + bw=bins.binWidth(0); // only costant supported. Need to figure out for zero bins + } + + //double lambda=pdf_->getVal()*bw; // sum coeff*pdf + double lambdat=0.0; + double lambda=0.;// debug: recompute lambda with loop. not sure that the term before include everything I want + + for (unsigned int i=0;ipdfs_.size();++i) // loop over processes + { + // find rate param. TODO make it more efficient using rate_pos + ProcessNormalization *c = dynamic_cast(pdf_->coeffs_[i]); + RooRealVar *rv = dynamic_cast(pdf_->coeffs_[i]); + + // fix RooProduct + RooProduct *pp = dynamic_cast(pdf_->coeffs_[i]); + if (c == nullptr and pp != nullptr){ + for ( int ip =0 ;ipcomponents().getSize(); ++ip) + { + if ( dynamic_cast( &pp->components()[ip]) != nullptr) c = dynamic_cast( &pp->components()[ip]); + } + } + RooRealVar *rate = dynamic_cast( (rate_pos_[i]==-2)? rv: (rate_pos_[i]>=0) ? &c->otherFactorList_[rate_pos_[i]] : nullptr ); + + //--- + int diracDelta = (rate_pos_[i]>=0 or rate_pos_[i]==-2)?1 : 0; + if (pdf_->isRooRealSum_ and verbose)std::cout<<"[DerivativeRateParam]:DEBUG "<<" pdf is a RooRealSumPdf"<basicIntegrals_<getVal()*bw <<")" <<" lambdat="< my_RooBernsteinFast_1; diff --git a/src/classes_def.xml b/src/classes_def.xml index 5bb2496a6ab..4927b4646c1 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -219,4 +219,6 @@ + +