Merge branch '102x' into 102x_port675

interface/CMSHggFormula.h

@@ -0,0 +1,156 @@
+#ifndef HiggsAnalysis_CombinedLimit_CMSHggFormula_h
+#define HiggsAnalysis_CombinedLimit_CMSHggFormula_h
+#include "RooAbsReal.h"
+#include "RooListProxy.h"
+#include "RooRealProxy.h"
+#include <vector>
+
+class CMSHggFormulaA1 : public RooAbsReal {
+ public:
+  CMSHggFormulaA1() {}
+  CMSHggFormulaA1(const char* name, const char* title, RooAbsReal & p0, RooAbsReal & p1,
+                 RooAbsReal & p2, RooAbsReal & p3, RooArgList const& terms,
+                 std::vector<double> const& coeffs);
+  CMSHggFormulaA1(const CMSHggFormulaA1& other, const char* name = 0);
+  virtual ~CMSHggFormulaA1() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaA1(*this, newname); }
+
+ protected:
+  RooRealProxy p0_;
+  RooRealProxy p1_;
+  RooRealProxy p2_;
+  RooRealProxy p3_;
+  RooListProxy terms_;
+  std::vector<double> coeffs_;
+  mutable std::vector<RooAbsReal*> vterms_; //! not to be serialized
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaA1,1)
+
+};
+
+class CMSHggFormulaA2 : public RooAbsReal {
+ public:
+  CMSHggFormulaA2() {}
+  CMSHggFormulaA2(const char* name, const char* title, RooAbsReal & p0, double const& p1,
+                 RooAbsReal & p2, RooAbsReal & p3, RooArgList const& terms,
+                 std::vector<double> const& coeffs);
+  CMSHggFormulaA2(const CMSHggFormulaA2& other, const char* name = 0);
+  virtual ~CMSHggFormulaA2() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaA2(*this, newname); }
+
+ protected:
+  RooRealProxy p0_;
+  double p1_;
+  RooRealProxy p2_;
+  RooRealProxy p3_;
+  RooListProxy terms_;
+  std::vector<double> coeffs_;
+  mutable std::vector<RooAbsReal*> vterms_; //! not to be serialized
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaA2,1)
+
+};
+
+class CMSHggFormulaB1 : public RooAbsReal {
+ public:
+  CMSHggFormulaB1() {}
+  CMSHggFormulaB1(const char* name, const char* title, RooAbsReal & p0, RooArgList const& terms,
+                 std::vector<double> const& coeffs);
+  CMSHggFormulaB1(const CMSHggFormulaB1& other, const char* name = 0);
+  virtual ~CMSHggFormulaB1() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaB1(*this, newname); }
+
+ protected:
+  RooRealProxy p0_;
+  RooListProxy terms_;
+  std::vector<double> coeffs_;
+  mutable std::vector<RooAbsReal*> vterms_; //! not to be serialized
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaB1,1)
+
+};
+
+class CMSHggFormulaB2 : public RooAbsReal {
+ public:
+  CMSHggFormulaB2() {}
+  CMSHggFormulaB2(const char* name, const char* title, double const& p0, RooArgList const& terms,
+                 std::vector<double> const& coeffs);
+  CMSHggFormulaB2(const CMSHggFormulaB2& other, const char* name = 0);
+  virtual ~CMSHggFormulaB2() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaB2(*this, newname); }
+
+ protected:
+  double p0_;
+  RooListProxy terms_;
+  std::vector<double> coeffs_;
+  mutable std::vector<RooAbsReal*> vterms_; //! not to be serialized
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaB2,1)
+
+};
+
+class CMSHggFormulaC1 : public RooAbsReal {
+ public:
+  CMSHggFormulaC1() {}
+  CMSHggFormulaC1(const char* name, const char* title, RooArgList const& terms,
+                 std::vector<double> const& coeffs);
+  CMSHggFormulaC1(const CMSHggFormulaC1& other, const char* name = 0);
+  virtual ~CMSHggFormulaC1() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaC1(*this, newname); }
+
+ protected:
+  RooListProxy terms_;
+  std::vector<double> coeffs_;
+  mutable std::vector<RooAbsReal*> vterms_; //! not to be serialized
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaC1,1)
+
+};
+
+class CMSHggFormulaD1 : public RooAbsReal {
+ public:
+  CMSHggFormulaD1() {}
+  CMSHggFormulaD1(const char* name, const char* title, RooAbsReal & p0, RooAbsReal & p1);
+  CMSHggFormulaD1(const CMSHggFormulaD1& other, const char* name = 0);
+  virtual ~CMSHggFormulaD1() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaD1(*this, newname); }
+
+ protected:
+  RooRealProxy p0_;
+  RooRealProxy p1_;
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaD1,1)
+
+};
+
+class CMSHggFormulaD2 : public RooAbsReal {
+ public:
+  CMSHggFormulaD2() {}
+  CMSHggFormulaD2(const char* name, const char* title, RooAbsReal & p0, double const& p1);
+  CMSHggFormulaD2(const CMSHggFormulaD2& other, const char* name = 0);
+  virtual ~CMSHggFormulaD2() {}
+  virtual TObject* clone(const char* newname) const { return new CMSHggFormulaD2(*this, newname); }
+
+ protected:
+  RooRealProxy p0_;
+  double p1_;
+  virtual Double_t evaluate() const;
+
+ private:
+  ClassDef(CMSHggFormulaD2,1)
+
+};
+
+#endif

interface/CMSHistFunc.h

@@ -146,6 +146,7 @@ class CMSHistFunc : public RooAbsReal {
 
   static void EnableFastVertical();
   friend class CMSHistV<CMSHistFunc>;
+  friend class CMSHistSum;
 
   /*
 

interface/CMSHistSum.h

@@ -0,0 +1,149 @@
+#ifndef CMSHistSum_h
+#define CMSHistSum_h
+#include <ostream>
+#include <vector>
+#include <memory>
+#include "RooAbsReal.h"
+#include "RooArgSet.h"
+#include "RooListProxy.h"
+#include "RooRealProxy.h"
+#include "Rtypes.h"
+#include "TH1F.h"
+#include "HiggsAnalysis/CombinedLimit/interface/FastTemplate_Old.h"
+#include "HiggsAnalysis/CombinedLimit/interface/SimpleCacheSentry.h"
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHistFunc.h"
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHistV.h"
+
+class CMSHistSum : public RooAbsReal {
+private:
+  struct BarlowBeeston {
+    bool init = false;
+    std::vector<unsigned> use;
+    std::vector<double> dat;
+    std::vector<double> valsum;
+    std::vector<double> toterr;
+    std::vector<double> err;
+    std::vector<double> b;
+    std::vector<double> c;
+    std::vector<double> tmp;
+    std::vector<double> x1;
+    std::vector<double> x2;
+    std::vector<double> res;
+    std::vector<double> gobs;
+    std::set<RooAbsArg*> dirty_prop;
+    std::vector<RooRealVar*> push_res;
+  };
+public:
+
+  CMSHistSum();
+
+  CMSHistSum(const char* name, const char* title, RooRealVar& x,
+                         RooArgList const& funcs, RooArgList const& coeffs);
+
+  CMSHistSum(CMSHistSum const& other, const char* name = 0);
+
+  virtual TObject* clone(const char* newname) const {
+    return new CMSHistSum(*this, newname);
+  }
+
+  virtual ~CMSHistSum() {;}
+
+  Double_t evaluate() const;
+
+  std::map<std::string, Double_t> getProcessNorms() const;
+
+  RooArgList * setupBinPars(double poissonThreshold);
+
+  std::unique_ptr<RooArgSet> getSentryArgs() const;
+
+  void printMultiline(std::ostream& os, Int_t contents, Bool_t verbose,
+                      TString indent) const;
+
+  Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars,
+                              const char* rangeName = 0) const;
+
+  Double_t analyticalIntegral(Int_t code, const char* rangeName = 0) const;
+
+  void setData(RooAbsData const& data) const;
+
+  void setAnalyticBarlowBeeston(bool flag) const;
+
+  inline FastHisto const& cache() const { return cache_; }
+
+  RooArgList const& coefList() const { return coeffpars_; }
+  // RooArgList const& funcList() const { return funcs_; }
+
+  static void EnableFastVertical();
+  friend class CMSHistV<CMSHistSum>;
+
+ protected:
+  RooRealProxy x_;
+
+  RooListProxy morphpars_;
+  RooListProxy coeffpars_;
+  RooListProxy binpars_;
+
+  int n_procs_;
+  int n_morphs_;
+
+  std::vector<FastTemplate> storage_;  // All nominal and vmorph templates
+  std::vector<int> process_fields_; // Indicies for process templates in storage_
+  std::vector<int> vmorph_fields_; // Indicies for vmorph templates in storage_
+
+  std::vector<FastTemplate> binerrors_; // Bin errors for each process
+
+  std::vector<CMSHistFunc::VerticalSetting> vtype_; // Vertical morphing type for each process
+  std::vector<double> vsmooth_par_; // Vertical morphing smooth region for each process
+
+  mutable std::vector<CMSHistFunc const*> vfuncstmp_; //!
+  mutable std::vector<RooAbsReal const*> vcoeffpars_; //!
+  mutable std::vector<RooAbsReal const*> vmorphpars_; //!
+  mutable std::vector<std::vector<RooAbsReal *>> vbinpars_; //!
+  std::vector<std::vector<unsigned>> bintypes_;
+
+  mutable std::vector<double> coeffvals_; //!
+
+  mutable std::vector<FastHisto> compcache_; //!
+  mutable FastHisto staging_; //!
+  mutable FastHisto valsum_; //!
+  mutable FastHisto cache_;
+
+  mutable std::vector<double> err2sum_; //!
+  mutable std::vector<double> toterr_; //!
+  mutable std::vector<std::vector<double>> binmods_; //!
+  mutable std::vector<std::vector<double>> scaledbinmods_; //!
+
+  mutable SimpleCacheSentry sentry_; //!
+  mutable SimpleCacheSentry binsentry_; //!
+
+  mutable std::vector<double> data_; //!
+
+  mutable BarlowBeeston bb_; //!
+
+  mutable bool initialized_; //! not to be serialized
+
+  mutable bool analytic_bb_; //! not to be serialized
+
+  mutable std::vector<double> vertical_prev_vals_; //! not to be serialized
+  mutable int fast_mode_; //! not to be serialized
+  static bool enable_fast_vertical_; //! not to be serialized
+
+  inline int& morphField(int const& ip, int const& iv) {
+    return vmorph_fields_[ip * n_morphs_ + iv];
+  }
+
+  void initialize() const;
+  void updateCache() const;
+  inline double smoothStepFunc(double x, int const& ip) const;
+
+
+  void runBarlowBeeston() const;
+
+  void updateMorphs() const;
+
+
+ private:
+  ClassDef(CMSHistSum,1)
+};
+
+#endif

interface/ProcessNormalization.h

@@ -40,6 +40,9 @@ class ProcessNormalization : public RooAbsReal {
         std::vector<std::pair<double,double> > logAsymmKappa_; // Logarithm of asymmetric kappas (low, high)
         RooListProxy asymmThetaList_;                           // List of nuisances for asymmetric kappas
         RooListProxy otherFactorList_;     // Other multiplicative terms 
+        std::vector<RooAbsReal *> thetaListVec_; //! Don't serialize me
+        std::vector<RooAbsReal *> asymmThetaListVec_; //! Don't serialize me
+        std::vector<RooAbsReal *> otherFactorListVec_; //! Don't serialize me
 
         // get the kappa for the appropriate x
         Double_t logKappaForX(double x, const std::pair<double,double> &logKappas ) const ;

interface/RooCheapProduct.h

@@ -14,9 +14,11 @@ class RooCheapProduct : public RooAbsReal {
         const RooArgList & components() const { return terms_; }
     protected:
         RooListProxy terms_;
-        std::vector<RooAbsReal *> vterms_;
+        std::vector<RooAbsReal *> vterms_; //! not to be serialized
         double offset_;
         virtual Double_t evaluate() const ;
+    private:
+        ClassDef(RooCheapProduct,1)
 };
 
 #endif

interface/VerticalInterpPdf.h

@@ -28,6 +28,8 @@ class VerticalInterpPdf : public RooAbsPdf {
   const RooArgList& funcList() const { return _funcList ; }
   const RooArgList& coefList() const { return _coefList ; }
 
+  void setFloorVals(Double_t const& pdf_val, Double_t const& integral_val);
+
 protected:
 
   class CacheElem : public RooAbsCacheElement {
@@ -47,13 +49,15 @@ class VerticalInterpPdf : public RooAbsPdf {
   TIterator* _funcIter ;     //! Iterator over FUNC list
   TIterator* _coefIter ;    //! Iterator over coefficient list
 
+  Double_t _pdfFloorVal; // PDF floor should be customizable, default is 1e-15
+  Double_t _integralFloorVal; // PDF integral floor should be customizable, default is 1e-10
+
   Double_t interpolate(Double_t coeff, Double_t central, RooAbsReal *fUp, RooAbsReal *fDown) const ; 
 
   bool isConditionalProdPdf(RooAbsReal *pdf) const;
   RooAbsReal* makeConditionalProdPdfIntegral(RooAbsPdf* pdf, RooArgSet const& analVars) const;
 
-private:
-  ClassDef(VerticalInterpPdf,2) // PDF constructed from a sum of (non-pdf) functions
+  ClassDef(VerticalInterpPdf,3) // PDF constructed from a sum of (non-pdf) functions
 };
 
 #endif

python/DatacardParser.py

@@ -19,6 +19,7 @@ def addDatacardParserOptions(parser):
     parser.add_option("--default-morphing",  dest="defMorph", type="string", default="shape", help="Default template morphing algorithm (to be used when the datacard has just 'shape')")
     parser.add_option("--no-b-only","--for-fits",    dest="noBOnly", default=False, action="store_true", help="Do not save the background-only pdf (saves time)")
     parser.add_option("--no-wrappers", dest="noHistFuncWrappers", default=False, action="store_true", help="Do not create and save the CMSHistFuncWrapper objects for autoMCStats-based models (saves time)")
+    parser.add_option("--use-histsum", dest="useCMSHistSum", default=False, action="store_true", help="Use memory-optimized CMSHistSum instead of CMSHistErrorPropagator")
     parser.add_option("--no-optimize-pdfs",    dest="noOptimizePdf", default=False, action="store_true", help="Do not save the RooSimultaneous as RooSimultaneousOpt and Gaussian constraints as SimpleGaussianConstraint")
     parser.add_option("--optimize-simpdf-constraints",    dest="moreOptimizeSimPdf", default="none", type="string", help="Handling of constraints in simultaneous pdf: 'none' = add all constraints on all channels (default); 'lhchcg' = add constraints on only the first channel; 'cms' = add constraints to the RooSimultaneousOpt.")
     #parser.add_option("--use-HistPdf",  dest="useHistPdf", type="string", default="always", help="Use RooHistPdf for TH1s: 'always' (default), 'never', 'when-constant' (i.e. not when doing template morphing)")

python/ShapeTools.py

@@ -142,8 +142,12 @@ def doIndividualModels(self):
                 else:
                     sigcoeffs.append(coeff)
             if self.options.verbose > 1: print "Creating RooAddPdf %s with %s elements" % ("pdf_bin"+b, coeffs.getSize())
-            if channelBinParFlag: 
-                prop = self.addObj(ROOT.CMSHistErrorPropagator, "prop_bin%s" % b, "", pdfs.at(0).getXVar(), pdfs, coeffs)
+            if channelBinParFlag:
+                if self.options.useCMSHistSum:
+                    prop = self.addObj(ROOT.CMSHistSum, "prop_bin%s" % b, "", pdfs.at(0).getXVar(), pdfs, coeffs)
+                    prop.setAttribute('CachingPdf_NoClone', True)
+                else:
+                    prop = self.addObj(ROOT.CMSHistErrorPropagator, "prop_bin%s" % b, "", pdfs.at(0).getXVar(), pdfs, coeffs)
                 prop.setAttribute('CachingPdf_Direct', True)
                 if self.DC.binParFlags[b][0] >= 0.:
                     bbb_args = prop.setupBinPars(self.DC.binParFlags[b][0])