112x semianalytical gradient

interface/CMSHistSum.h

@@ -14,8 +14,15 @@
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistFunc.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistV.h"
 
+class DerivativeRateParamCMSHistSum;
+class DerivativeLogNormalCMSHistSum;
+class SimNLLDerivativesHelper;
+
 class CMSHistSum : public RooAbsReal {
 private:
+  friend DerivativeRateParamCMSHistSum;
+  friend DerivativeLogNormalCMSHistSum;
+  friend SimNLLDerivativesHelper;
   struct BarlowBeeston {
     bool init = false;
     std::vector<unsigned> use;

interface/CachingNLL.h

@@ -21,6 +21,12 @@
 #include <boost/ptr_container/ptr_vector.hpp>
 
 class RooMultiPdf;
+class SimNLLDerivativesHelper;
+class DerivativeAbstract; // ideally we would like only this one. TODO
+class DerivativeLogNormal;
+class DerivativeRateParam;
+class DerivativeLogNormalCMSHistSum;
+class DerivativeRateParamCMSHistSum;
 
 // Part zero: ArgSet checker
 namespace cacheutils {
@@ -112,6 +118,12 @@ class OptimizedCachingPdfT : public CachingPdf {
 CachingPdfBase * makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ;
 
 class CachingAddNLL : public RooAbsReal {
+    friend SimNLLDerivativesHelper; // probably not needed w/ data
+    friend DerivativeAbstract;
+    friend DerivativeLogNormal;
+    friend DerivativeRateParam;
+    friend DerivativeLogNormalCMSHistSum;
+    friend DerivativeRateParamCMSHistSum;
     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 +137,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 +199,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_;

interface/CascadeMinimizer.h

@@ -12,6 +12,9 @@ class RooRealVar;
 #include <boost/program_options.hpp>
 #include <boost/algorithm/string.hpp>
 
+#include "HiggsAnalysis/CombinedLimit/interface/RooMinimizerSemiAnalytic.h"
+#include "HiggsAnalysis/CombinedLimit/interface/SimNLLDerivativesHelper.h"
+#include <stdexcept>
 
 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_; }
@@ -43,7 +46,14 @@ class CascadeMinimizer {
 	std::string algo() {return defaultMinimizerAlgo_;};
     private:
         RooAbsReal & nll_;
+
         std::unique_ptr<RooMinimizer> minimizer_;
+        //
+        bool isSemiAnalyticMinimizer{false};
+        std::unique_ptr<RooMinimizerSemiAnalytic> minimizerSemiAnalytic_;
+        std::map<std::string,RooAbsReal*> derivatives_; // not used
+        std::unique_ptr<SimNLLDerivativesHelper> helper_;
+        //
         Mode         mode_;
         static int          strategy_;
         RooRealVar * poi_; 

interface/ProcessNormalization.h

@@ -14,6 +14,13 @@ ProcessNormalization is helper class for implementing process normalizations
 END_HTML
 */
 //
+class SimNLLDerivativesHelper;
+class DerivativeAbstract;
+class DerivativeLogNormal;
+class DerivativeRateParam;
+class DerivativeLogNormalCMSHistSum;
+class DerivativeRateParamCMSHistSum;
+
 class ProcessNormalization : public RooAbsReal {
    public:
       ProcessNormalization() : nominalValue_(1) {}
@@ -33,6 +40,12 @@ class ProcessNormalization : public RooAbsReal {
         Double_t evaluate() const;
 
     private:
+        friend SimNLLDerivativesHelper;
+        friend DerivativeAbstract;
+        friend DerivativeLogNormal;
+        friend DerivativeRateParam;
+        friend DerivativeLogNormalCMSHistSum;
+        friend DerivativeRateParamCMSHistSum;
         // ---- PERSISTENT ----
         double nominalValue_;                         
         std::vector<double> logKappa_; // Logarithm of symmetric kappas

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, [email protected]       *
+ *   DK, David Kirkby,    UC Irvine,         [email protected]                 *
+ *                                                                           *
+ * 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 <vector>
+#include <string>
+#include <map>
+
+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 &paramList);
+
+  // 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<double>& globalCC, const TMatrixDSym& corrs, const TMatrixDSym& covs) ;
+  void fillLegacyCorrMatrix() const ;
+  void fillPrefitCorrMatrix();
+  void setStatusHistory(std::vector<std::pair<std::string,int> >& 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<std::pair<std::string,int> > _statusHistory ; // History of status codes
+
+  ClassDef(RooFitResult,5) // Container class for fit result
+};
+
+#endif