diff --git a/interface/CMSHggFormula.h b/interface/CMSHggFormula.h
new file mode 100644
index 00000000000..4b0875d08da
--- /dev/null
+++ b/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
diff --git a/interface/CMSHistFunc.h b/interface/CMSHistFunc.h
index f03a31c42d8..3471ca372a9 100644
--- a/interface/CMSHistFunc.h
+++ b/interface/CMSHistFunc.h
@@ -146,6 +146,7 @@ class CMSHistFunc : public RooAbsReal {
 
   static void EnableFastVertical();
   friend class CMSHistV<CMSHistFunc>;
+  friend class CMSHistSum;
 
   /*
 
diff --git a/interface/CMSHistSum.h b/interface/CMSHistSum.h
new file mode 100644
index 00000000000..ae92954f11c
--- /dev/null
+++ b/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
diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h
index 0721e198f25..d99aec236e3 100644
--- a/interface/ProcessNormalization.h
+++ b/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 ;
diff --git a/interface/RooCheapProduct.h b/interface/RooCheapProduct.h
index 623736cc7c7..434bfe7b0f2 100644
--- a/interface/RooCheapProduct.h
+++ b/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
diff --git a/interface/VerticalInterpPdf.h b/interface/VerticalInterpPdf.h
index 63c358966b3..761b2f1194b 100644
--- a/interface/VerticalInterpPdf.h
+++ b/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
diff --git a/python/DatacardParser.py b/python/DatacardParser.py
index 91500b332e4..9748e862d54 100644
--- a/python/DatacardParser.py
+++ b/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)")
diff --git a/python/ShapeTools.py b/python/ShapeTools.py
index 87b323b7786..d0c3c354111 100644
--- a/python/ShapeTools.py
+++ b/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])
diff --git a/src/CMSHggFormula.cc b/src/CMSHggFormula.cc
new file mode 100644
index 00000000000..cb58bbcaa3d
--- /dev/null
+++ b/src/CMSHggFormula.cc
@@ -0,0 +1,259 @@
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHggFormula.h"
+#include "RooConstVar.h"
+
+CMSHggFormulaA1::CMSHggFormulaA1(const char* name, const char* title, RooAbsReal& p0, RooAbsReal& p1,
+                               RooAbsReal& p2, RooAbsReal& p3, RooArgList const& terms,
+                               std::vector<double> const& coeffs)
+    : RooAbsReal(name, title),
+      p0_("p0", "p0", this, p0),
+      p1_("p1", "p1", this, p1),
+      p2_("p2", "p2", this, p2),
+      p3_("p3", "p3", this, p3),
+      terms_("terms", "terms", this),
+      coeffs_(coeffs) {
+  RooFIter iter = terms.fwdIterator();
+  for (RooAbsArg* a = iter.next(); a != 0; a = iter.next()) {
+    RooAbsReal* rar = dynamic_cast<RooAbsReal*>(a);
+    if (!rar) {
+      throw std::invalid_argument(std::string("Component ") + a->GetName() +
+                                  " of CMSHggFormulaA1 is a " + a->ClassName());
+    }
+    terms_.add(*rar);
+    vterms_.push_back(rar);
+  }
+  if (terms_.getSize() != int(coeffs.size())) {
+    throw std::invalid_argument("Terms and coeffs not of equal size in CMSHggFormulaA1");
+  }
+}
+
+CMSHggFormulaA1::CMSHggFormulaA1(const CMSHggFormulaA1& other, const char* name)
+    : RooAbsReal(other, name),
+      p0_("p0", this, other.p0_),
+      p1_("p1", this, other.p1_),
+      p2_("p2", this, other.p2_),
+      p3_("p3", this, other.p3_),
+      terms_("terms", this, other.terms_),
+      coeffs_(other.coeffs_),
+      vterms_(other.vterms_) {}
+
+Double_t CMSHggFormulaA1::evaluate() const {
+  if (vterms_.empty()) {
+    vterms_.resize(terms_.getSize());
+    for (int i = 0; i < terms_.getSize(); ++i) {
+      vterms_[i] = dynamic_cast<RooAbsReal*>(terms_.at(i));
+    }
+  }
+  // (@0+@1)*(1.+@2+@3+@4*@5+@6*@7+@8*@9+@10*@11+@12*@13+@14*@15+@16*@17+@18*@19+@20*@21+@22*@23)
+  double ret = 1. + p2_ + p3_;
+  for (unsigned i = 0; i < coeffs_.size(); ++i) {
+    ret += coeffs_[i] * vterms_[i]->getVal();
+  }
+  return (p0_ + p1_) * ret;
+}
+
+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)
+    : RooAbsReal(name, title),
+      p0_("p0", "p0", this, p0),
+      p1_(p1),
+      p2_("p2", "p2", this, p2),
+      p3_("p3", "p3", this, p3),
+      terms_("terms", "terms", this),
+      coeffs_(coeffs) {
+  RooFIter iter = terms.fwdIterator();
+  for (RooAbsArg* a = iter.next(); a != 0; a = iter.next()) {
+    RooAbsReal* rar = dynamic_cast<RooAbsReal*>(a);
+    if (!rar) {
+      throw std::invalid_argument(std::string("Component ") + a->GetName() +
+                                  " of CMSHggFormulaA2 is a " + a->ClassName());
+    }
+    terms_.add(*rar);
+    vterms_.push_back(rar);
+  }
+  if (terms_.getSize() != int(coeffs.size())) {
+    throw std::invalid_argument("Terms and coeffs not of equal size in CMSHggFormulaA2");
+  }
+}
+
+CMSHggFormulaA2::CMSHggFormulaA2(const CMSHggFormulaA2& other, const char* name)
+    : RooAbsReal(other, name),
+      p0_("p0", this, other.p0_),
+      p1_(other.p1_),
+      p2_("p2", this, other.p2_),
+      p3_("p3", this, other.p3_),
+      terms_("terms", this, other.terms_),
+      coeffs_(other.coeffs_),
+      vterms_(other.vterms_) {}
+
+Double_t CMSHggFormulaA2::evaluate() const {
+  if (vterms_.empty()) {
+    vterms_.resize(terms_.getSize());
+    for (int i = 0; i < terms_.getSize(); ++i) {
+      vterms_[i] = dynamic_cast<RooAbsReal*>(terms_.at(i));
+    }
+  }
+  // (@0+@1)*(1.+@2+@3+@4*@5+@6*@7+@8*@9+@10*@11+@12*@13+@14*@15+@16*@17+@18*@19+@20*@21+@22*@23)
+  double ret = 1. + p2_ + p3_;
+  for (unsigned i = 0; i < coeffs_.size(); ++i) {
+    ret += coeffs_[i] * vterms_[i]->getVal();
+  }
+  return (p0_ + p1_) * ret;
+}
+
+CMSHggFormulaB1::CMSHggFormulaB1(const char* name, const char* title, RooAbsReal& p0,
+                               RooArgList const& terms, std::vector<double> const& coeffs)
+    : RooAbsReal(name, title),
+      p0_("p0", "p0", this, p0),
+      terms_("terms", "terms", this),
+      coeffs_(coeffs) {
+  RooFIter iter = terms.fwdIterator();
+  for (RooAbsArg* a = iter.next(); a != 0; a = iter.next()) {
+    RooAbsReal* rar = dynamic_cast<RooAbsReal*>(a);
+    if (!rar) {
+      throw std::invalid_argument(std::string("Component ") + a->GetName() +
+                                  " of CMSHggFormulaB1 is a " + a->ClassName());
+    }
+    terms_.add(*rar);
+    vterms_.push_back(rar);
+  }
+  if (terms_.getSize() != int(coeffs.size())) {
+    throw std::invalid_argument("Terms and coeffs not of equal size in CMSHggFormulaB1");
+  }
+}
+
+CMSHggFormulaB1::CMSHggFormulaB1(const CMSHggFormulaB1& other, const char* name)
+    : RooAbsReal(other, name),
+      p0_("p0", this, other.p0_),
+      terms_("terms", this, other.terms_),
+      coeffs_(other.coeffs_),
+      vterms_(other.vterms_) {}
+
+Double_t CMSHggFormulaB1::evaluate() const {
+  if (vterms_.empty()) {
+    vterms_.resize(terms_.getSize());
+    for (int i = 0; i < terms_.getSize(); ++i) {
+      vterms_[i] = dynamic_cast<RooAbsReal*>(terms_.at(i));
+    }
+  }
+  // @0*TMath::Max(1.e-2,(1.+@1*@2+@3*@4+@5*@6+@7*@8+@9*@10+@11*@12+@13*@14+@15*@16+@17*@18+@19*@20+@21*@22+@23*@24+@25*@26+@27*@28+@29*@30+@31*@32+@33*@34+@35*@36+@37*@38+@39*@40+@41*@42+@43*@44))
+  double ret = 1.;
+  for (unsigned i = 0; i < coeffs_.size(); ++i) {
+    ret += coeffs_[i] * vterms_[i]->getVal();
+  }
+  return p0_ * TMath::Max(1.e-2, ret);
+}
+
+CMSHggFormulaB2::CMSHggFormulaB2(const char* name, const char* title, double const& p0,
+                               RooArgList const& terms, std::vector<double> const& coeffs)
+    : RooAbsReal(name, title),
+      p0_(p0),
+      terms_("terms", "terms", this),
+      coeffs_(coeffs) {
+  RooFIter iter = terms.fwdIterator();
+  for (RooAbsArg* a = iter.next(); a != 0; a = iter.next()) {
+    RooAbsReal* rar = dynamic_cast<RooAbsReal*>(a);
+    if (!rar) {
+      throw std::invalid_argument(std::string("Component ") + a->GetName() +
+                                  " of CMSHggFormulaB2 is a " + a->ClassName());
+    }
+    terms_.add(*rar);
+    vterms_.push_back(rar);
+  }
+  if (terms_.getSize() != int(coeffs.size())) {
+    throw std::invalid_argument("Terms and coeffs not of equal size in CMSHggFormulaB2");
+  }
+}
+
+CMSHggFormulaB2::CMSHggFormulaB2(const CMSHggFormulaB2& other, const char* name)
+    : RooAbsReal(other, name),
+      p0_(other.p0_),
+      terms_("terms", this, other.terms_),
+      coeffs_(other.coeffs_),
+      vterms_(other.vterms_) {}
+
+Double_t CMSHggFormulaB2::evaluate() const {
+  if (vterms_.empty()) {
+    vterms_.resize(terms_.getSize());
+    for (int i = 0; i < terms_.getSize(); ++i) {
+      vterms_[i] = dynamic_cast<RooAbsReal*>(terms_.at(i));
+    }
+  }
+  // @0*TMath::Max(1.e-2,(1.+@1*@2+@3*@4+@5*@6+@7*@8+@9*@10+@11*@12+@13*@14+@15*@16+@17*@18+@19*@20+@21*@22+@23*@24+@25*@26+@27*@28+@29*@30+@31*@32+@33*@34+@35*@36+@37*@38+@39*@40+@41*@42+@43*@44))
+  double ret = 1.;
+  for (unsigned i = 0; i < coeffs_.size(); ++i) {
+    ret += coeffs_[i] * vterms_[i]->getVal();
+  }
+  return p0_ * TMath::Max(1.e-2, ret);
+}
+
+CMSHggFormulaC1::CMSHggFormulaC1(const char* name, const char* title, RooArgList const& terms,
+                               std::vector<double> const& coeffs)
+    : RooAbsReal(name, title),
+      terms_("terms", "terms", this),
+      coeffs_(coeffs) {
+  RooFIter iter = terms.fwdIterator();
+  for (RooAbsArg* a = iter.next(); a != 0; a = iter.next()) {
+    RooAbsReal* rar = dynamic_cast<RooAbsReal*>(a);
+    if (!rar) {
+      throw std::invalid_argument(std::string("Component ") + a->GetName() +
+                                  " of CMSHggFormulaC1 is a " + a->ClassName());
+    }
+    terms_.add(*rar);
+    vterms_.push_back(rar);
+  }
+  if (terms_.getSize() != int(coeffs.size())) {
+    throw std::invalid_argument("Terms and coeffs not of equal size in CMSHggFormulaC1");
+  }
+}
+
+CMSHggFormulaC1::CMSHggFormulaC1(const CMSHggFormulaC1& other, const char* name)
+    : RooAbsReal(other, name),
+      terms_("terms", this, other.terms_),
+      coeffs_(other.coeffs_),
+      vterms_(other.vterms_) {}
+
+Double_t CMSHggFormulaC1::evaluate() const {
+  if (vterms_.empty()) {
+    vterms_.resize(terms_.getSize());
+    for (int i = 0; i < terms_.getSize(); ++i) {
+      vterms_[i] = dynamic_cast<RooAbsReal*>(terms_.at(i));
+    }
+  }
+  // (1.+@0*@1+@2*@3+@4*@5+@6*@7+@8*@9+@10*@11+@12*@13) @0[C],@1[V],@2[C],@3[V],@4[C],@5[V],@6[C],@7[V],@8[C],@9[V],@10[C],@11[V],@12[C],@13[V]
+  double ret = 1.;
+  for (unsigned i = 0; i < coeffs_.size(); ++i) {
+    ret += coeffs_[i] * vterms_[i]->getVal();
+  }
+  return ret;
+}
+
+CMSHggFormulaD1::CMSHggFormulaD1(const char* name, const char* title, RooAbsReal& p0, RooAbsReal& p1)
+    : RooAbsReal(name, title),
+      p0_("p0", "p0", this, p0),
+      p1_("p1", "p1", this, p1) { }
+
+CMSHggFormulaD1::CMSHggFormulaD1(const CMSHggFormulaD1& other, const char* name)
+    : RooAbsReal(other, name),
+      p0_("p0", this, other.p0_),
+      p1_("p1", this, other.p1_) {}
+
+Double_t CMSHggFormulaD1::evaluate() const {
+  // TMath::Min(@0+@1,1.0) @0[V],@1[V]
+  return TMath::Min(p0_+p1_,1.0);
+}
+
+CMSHggFormulaD2::CMSHggFormulaD2(const char* name, const char* title, RooAbsReal& p0, double const& p1)
+    : RooAbsReal(name, title),
+      p0_("p0", "p0", this, p0),
+      p1_(p1) { }
+
+CMSHggFormulaD2::CMSHggFormulaD2(const CMSHggFormulaD2& other, const char* name)
+    : RooAbsReal(other, name),
+      p0_("p0", this, other.p0_),
+      p1_(other.p1_) {}
+
+Double_t CMSHggFormulaD2::evaluate() const {
+  // TMath::Min(@0+@1,1.0) @0[V],@1[V]
+  return TMath::Min(p0_+p1_,1.0);
+}
diff --git a/src/CMSHistSum.cc b/src/CMSHistSum.cc
new file mode 100644
index 00000000000..7e065aa9b21
--- /dev/null
+++ b/src/CMSHistSum.cc
@@ -0,0 +1,758 @@
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHistSum.h"
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHistFuncWrapper.h"
+#include <stdexcept>
+#include <vector>
+#include <ostream>
+#include <memory>
+#include "Math/ProbFuncMathCore.h"
+#include "Math/QuantFuncMathCore.h"
+#include "RooRealProxy.h"
+#include "RooArgSet.h"
+#include "RooAbsReal.h"
+#include "RooAbsData.h"
+#include "RooConstVar.h"
+#include "RooGaussian.h"
+#include "RooProduct.h"
+#include "vectorized.h"
+
+#define HFVERBOSE 0
+
+bool CMSHistSum::enable_fast_vertical_ = false;
+
+CMSHistSum::CMSHistSum() : initialized_(false), fast_mode_(0) {}
+
+CMSHistSum::CMSHistSum(const char* name,
+                                               const char* title,
+                                               RooRealVar& x,
+                                               RooArgList const& funcs,
+                                               RooArgList const& coeffs)
+    : RooAbsReal(name, title),
+      x_("x", "", this, x),
+      morphpars_("morphpars", "", this),
+      coeffpars_("coeffpars", "", this),
+      binpars_("binpars", "", this),
+      n_procs_(0),
+      n_morphs_(0),
+      sentry_(TString(name) + "_sentry", ""),
+      binsentry_(TString(name) + "_binsentry", ""),
+      initialized_(false),
+      analytic_bb_(false),
+      fast_mode_(0) {
+
+  n_procs_ = funcs.getSize();
+  assert(n_procs_ == coeffs.getSize());
+
+  coeffpars_.add(coeffs);
+
+  // Resize vectors
+  binerrors_.resize(n_procs_);
+  vtype_.resize(n_procs_);
+  vsmooth_par_.resize(n_procs_);
+
+  // First pass through the list of processes - we don't know the number of
+  // vertical morphing parameters yet.
+  std::map<std::string, RooRealVar const*> all_vmorphs;
+  std::vector<std::map<std::string, int>> vmorph_idx_mapping(n_procs_);
+  int n_storage = 0;
+  for (int ip = 0; ip < n_procs_; ++ip) {
+    CMSHistFunc const* func = dynamic_cast<CMSHistFunc const*>(funcs.at(ip));
+
+    // Verify that each func does not use horizontal morphing or rebinning
+    if (func->hmorphs_.getSize() > 0 || func->rebin_) {
+      throw std::runtime_error("CMSHistSum does not support horizontal morphing or on-the-fly rebinning");
+    }
+
+    // Add this CMSHistFunc to vfuncstmp_ - this is not persisted but will be used to create bbb
+    // parameters
+    vfuncstmp_.push_back(func);
+
+    // Need to initialize our cache_ with the correct binning
+    if (ip == 0) {
+      cache_ = func->cache();
+      cache_.Clear();
+    }
+
+    // Determine the size we'll need for storage_
+    n_storage += func->storage_.size();
+
+    // Construct a list of the unique vertical morphing parameters and
+    // make a per-process index of the local list positions of each vmorph
+    for (int iv = 0; iv < func->vmorphs_.getSize(); ++iv) {
+      RooRealVar const* rrv = static_cast<RooRealVar const*>(func->vmorphs_.at(iv));
+      all_vmorphs[rrv->GetName()] = rrv;
+      vmorph_idx_mapping[ip][rrv->GetName()] = iv;
+    }
+
+    // Copy bin errors, vtype and vsmooth_par
+    binerrors_[ip] = func->binerrors_;
+    vtype_[ip] = func->vtype_;
+    vsmooth_par_[ip] = func->vsmooth_par_;
+  }
+
+  // Populate a list of the morphing parameter names,
+  // and add each parameter to the RooListProxy
+  std::vector<std::string> morph_names;
+  for (auto const& it : all_vmorphs) {
+    morph_names.push_back(it.first);
+    morphpars_.add(*it.second);
+  }
+
+  // Now we know the number of morphing parameters
+  n_morphs_ = morph_names.size();
+
+
+  // Now we populate storage_, by stealing from each CMSHistFunc in turn
+  storage_.clear();
+  storage_.resize(n_storage);
+  auto copy_it = storage_.begin();
+  for (int ip = 0; ip < n_procs_; ++ip) {
+    CMSHistFunc const* func = dynamic_cast<CMSHistFunc const*>(funcs.at(ip));
+#if HFVERBOSE > 0
+    std::cout << "Before copy: size = " << storage_.size() << ", distance = " << std::distance(storage_.begin(), copy_it) << "\n";
+#endif
+    copy_it = std::copy(func->storage_.begin(), func->storage_.end(), copy_it);
+#if HFVERBOSE > 0
+    std::cout << "After copy: size = " << storage_.size() << ", distance = " << std::distance(storage_.begin(), copy_it) << "\n";
+#endif
+  }
+
+  // Next populate the storage_ index positions for fast lookup later
+  process_fields_ = std::vector<int>(n_procs_, 0);
+  vmorph_fields_ = std::vector<int>(n_procs_ * n_morphs_, -1);
+  int current_proc = 0;
+  for (int ip = 0; ip < n_procs_; ++ip) {
+    process_fields_[ip] = current_proc;
+    // current_proc += 1;
+    for (int iv = 0; iv < n_morphs_; ++iv) {
+      auto it = vmorph_idx_mapping[ip].find(morph_names[iv]);
+      if (it != vmorph_idx_mapping[ip].end()) {
+        // Position for this vmorph
+        morphField(ip, iv) = current_proc + 1 + it->second * 2;
+        // Convert to sum/diff
+        unsigned idx = process_fields_[ip];
+        unsigned idxLo = morphField(ip, iv) + 0;
+        unsigned idxHi = morphField(ip, iv) + 1;
+        FastTemplate lo = storage_[idxLo];
+        FastTemplate hi = storage_[idxHi];
+        if (vtype_[ip] == CMSHistFunc::VerticalSetting::QuadLinear) {
+          hi.Subtract(storage_[idx]);
+          lo.Subtract(storage_[idx]);
+        } else if (vtype_[ip] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
+          hi.LogRatio(storage_[idx]);
+          lo.LogRatio(storage_[idx]);
+        }
+        FastTemplate::SumDiff(hi, lo, storage_[idxLo], storage_[idxHi]);
+      }
+    }
+    current_proc += (1 +vmorph_idx_mapping[ip].size() * 2);
+  }
+
+#if HFVERBOSE > 0
+  std::cout << "FUNCTIONS\n";
+  for (int ip = 0; ip < n_procs_; ++ip) {
+    std::cout << ip << ": " << funcs.at(ip)->GetName() << std::endl;
+  }
+  std::cout << "MORPHS\n";
+  for (int imorph = 0; imorph < n_morphs_; ++imorph) {
+    std::cout << imorph << ": " << morph_names[imorph] << std::endl;
+  }
+  for (int ip = 0; ip < n_procs_; ++ip) {
+    std::cout << TString::Format("%4i: %6i |", ip, process_fields_[ip]);
+    for (int iv = 0; iv < n_morphs_; ++iv) {
+      std::cout << TString::Format("%6i", morphField(ip, iv));
+    }
+    std::cout << std::endl;
+  }
+#endif
+
+  // exit(0);
+}
+
+CMSHistSum::CMSHistSum(
+    CMSHistSum const& other, const char* name)
+    : RooAbsReal(other, name),
+      x_("x", this, other.x_),
+      morphpars_("morphpars", this, other.morphpars_),
+      coeffpars_("coeffpars", this, other.coeffpars_),
+      binpars_("binpars", this, other.binpars_),
+      n_procs_(other.n_procs_),
+      n_morphs_(other.n_morphs_),
+      storage_(other.storage_),
+      process_fields_(other.process_fields_),
+      vmorph_fields_(other.vmorph_fields_),
+      binerrors_(other.binerrors_),
+      vtype_(other.vtype_),
+      vsmooth_par_(other.vsmooth_par_),
+      bintypes_(other.bintypes_),
+      cache_(other.cache_),
+      sentry_(name ? TString(name) + "_sentry" : TString(other.sentry_.GetName()), ""),
+      binsentry_(name ? TString(name) + "_binsentry" : TString(other.binsentry_.GetName()), ""),
+      initialized_(false),
+      fast_mode_(0) {
+      initialize();
+}
+
+void CMSHistSum::initialize() const {
+  if (initialized_) return;
+  sentry_.SetName(TString(this->GetName()) + "_sentry");
+  binsentry_.SetName(TString(this->GetName()) + "_binsentry");
+
+#if HFVERBOSE > 0
+  std::cout << "Initialising vectors\n";
+#endif
+
+  vmorphpars_.resize(n_morphs_);
+  vcoeffpars_.resize(n_procs_);
+  compcache_.resize(n_procs_);
+
+  for (int ip = 0; ip < n_procs_; ++ip) {
+    vcoeffpars_[ip] = dynamic_cast<RooAbsReal const*>(coeffpars_.at(ip));
+    compcache_[ip] = cache_;
+    compcache_[ip].Clear();
+  }
+  for (int iv = 0; iv < n_morphs_; ++iv) {
+    vmorphpars_[iv] = dynamic_cast<RooAbsReal const*>(morphpars_.at(iv));
+  }
+
+  unsigned nb = cache_.size();
+  vbinpars_.resize(nb);
+  if (bintypes_.size()) {
+    for (unsigned j = 0, r = 0; j < nb; ++j) {
+      vbinpars_[j].resize(bintypes_[j].size(), nullptr);
+      for (unsigned i = 0; i < bintypes_[j].size(); ++i) {
+        if (bintypes_[j][i] >= 1 && bintypes_[j][i] < 4) {
+          vbinpars_[j][i] = dynamic_cast<RooAbsReal *>(binpars_.at(r));
+          ++r;
+        }
+      }
+    }
+  }
+
+  valsum_ = cache_;
+  staging_ = cache_;
+  valsum_.Clear();
+  cache_.Clear();
+  staging_.Clear();
+  err2sum_.resize(nb, 0.);
+  toterr_.resize(nb, 0.);
+  binmods_.resize(n_procs_, std::vector<double>(nb, 0.));
+  scaledbinmods_.resize(n_procs_, std::vector<double>(nb, 0.));
+  coeffvals_.resize(n_procs_, 0.);
+
+  sentry_.addVars(morphpars_);
+  sentry_.addVars(coeffpars_);
+  binsentry_.addVars(binpars_);
+
+  sentry_.setValueDirty();
+  binsentry_.setValueDirty();
+
+  initialized_ = true;
+}
+
+void CMSHistSum::updateMorphs() const {
+  // If we're not in fast mode, need to reset all the compcache_
+  #if HFVERBOSE > 0
+  std::cout << "fast_mode_ = " << fast_mode_ << std::endl;
+  #endif
+  for (unsigned ip = 0; ip < compcache_.size(); ++ip) {
+    if (fast_mode_ == 0) {
+      compcache_[ip].CopyValues(storage_[process_fields_[ip]]);
+    }
+    if (vtype_[ip] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
+      compcache_[ip].Log();
+    }
+  }
+  int n_morphs = vmorphpars_.size();
+
+  if (vertical_prev_vals_.size() == 0) {
+    vertical_prev_vals_.resize(n_morphs);
+  }
+  // Loop through vmorphs
+  for (unsigned iv = 0; iv < vmorphpars_.size(); ++iv) {
+    double x = vmorphpars_[iv]->getVal();
+    if (fast_mode_ == 1 && (x == vertical_prev_vals_[iv])) {
+      #if HFVERBOSE > 0
+      std::cout << "Skipping " << vmorphpars_[iv]->GetName() << ", prev = now = " << x << std::endl;
+      #endif
+      continue;
+    }
+    #if HFVERBOSE > 0
+    if (fast_mode_ == 1) {
+      std::cout << "Updating " << vmorphpars_[iv]->GetName() << ", prev =  " << vertical_prev_vals_[iv] << ", now = " << x << std::endl;
+    }
+    #endif
+
+
+    // If in fast made, check if the value has changed since the last eval,
+    // and if it hasn't, skip it
+
+    // For each vmorph need to know the list of processes we apply to
+    for (unsigned ip = 0; ip < compcache_.size(); ++ip) {
+      int code = vmorph_fields_[ip * n_morphs + iv];
+      if (code == -1) continue;
+      if (fast_mode_ == 1) {
+        double xold = vertical_prev_vals_[iv];
+        compcache_[ip].DiffMeld(storage_[code + 1], storage_[code + 0], 0.5*x, smoothStepFunc(x, ip), 0.5*xold, smoothStepFunc(xold, ip));
+      } else {
+        compcache_[ip].Meld(storage_[code + 1], storage_[code + 0], 0.5*x, smoothStepFunc(x, ip));
+      }
+    }
+    vertical_prev_vals_[iv] = x;
+  }
+
+  if (enable_fast_vertical_) fast_mode_ = 1;
+}
+
+inline double CMSHistSum::smoothStepFunc(double x, int const& ip) const {
+  if (fabs(x) >= vsmooth_par_[ip]) return x > 0 ? +1 : -1;
+  double xnorm = x / vsmooth_par_[ip];
+  double xnorm2 = xnorm * xnorm;
+  return 0.125 * xnorm * (xnorm2 * (3. * xnorm2 - 10.) + 15);
+}
+
+
+void CMSHistSum::updateCache() const {
+  initialize();
+
+#if HFVERBOSE > 0
+  std::cout << "Sentry: " << sentry_.good() << "\n";
+#endif
+  if (!sentry_.good()) {
+    #if HFVERBOSE > 0
+      std::cout << "Calling updateMorphs\n";
+    #endif
+    updateMorphs();
+    for (unsigned i = 0; i < vcoeffpars_.size(); ++i) {
+      // vfuncs_[i]->updateCache();
+      coeffvals_[i] = vcoeffpars_[i]->getVal();
+    }
+    #if HFVERBOSE > 0
+      std::cout << "Updated coeffs\n";
+    #endif
+
+    valsum_.Clear();
+    std::fill(err2sum_.begin(), err2sum_.end(), 0.);
+    for (unsigned i = 0; i < vcoeffpars_.size(); ++i) {
+      staging_ = compcache_[i];
+      if (vtype_[i] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
+        staging_.Exp();
+      }
+      staging_.CropUnderflows();
+      vectorized::mul_add(valsum_.size(), coeffvals_[i], &(staging_[0]), &valsum_[0]);
+      vectorized::mul_add_sqr(valsum_.size(), coeffvals_[i], &(binerrors_[i][0]), &err2sum_[0]);
+    }
+    vectorized::sqrt(valsum_.size(), &err2sum_[0], &toterr_[0]);
+    cache_ = valsum_;
+    #if HFVERBOSE > 0
+      std::cout << "Updated cache\n";
+    #endif
+    sentry_.reset();
+    binsentry_.setValueDirty();
+  }
+
+
+  if (!binsentry_.good()) {
+    #if HFVERBOSE > 0
+      std::cout << "Calling runBarlowBeeston\n";
+    #endif
+    runBarlowBeeston();
+    // bintypes might have size == 0 if we never ran setupBinPars()
+    #if HFVERBOSE > 0
+      std::cout << "Assigning bin shifts\n";
+    #endif
+    for (unsigned j = 0; j < bintypes_.size(); ++j) {
+      cache_[j] = valsum_[j];
+      if (bintypes_[j][0] == 0) {
+        continue;
+      } else if (bintypes_[j][0] == 1) {
+        double x = vbinpars_[j][0]->getVal();
+        cache_[j] += toterr_[j] * x;
+      } else {
+        for (unsigned i = 0; i < bintypes_[j].size(); ++i) {
+          if (bintypes_[j][i] == 2) {
+            // Poisson: this is a multiplier on the process yield
+            scaledbinmods_[i][j] = ((vbinpars_[j][i]->getVal() - 1.) *
+                 compcache_[i][j]);
+            cache_[j] += (scaledbinmods_[i][j] * coeffvals_[i]);
+          } else if (bintypes_[j][i] == 3) {
+            // Gaussian This is the addition of the scaled error
+            scaledbinmods_[i][j] = vbinpars_[j][i]->getVal() * binerrors_[i][j];
+            cache_[j] += (scaledbinmods_[i][j] * coeffvals_[i]);
+          }
+        }
+      }
+    }
+    #if HFVERBOSE > 0
+      std::cout << "Done assigning bin shifts\n";
+    #endif
+    cache_.CropUnderflows();
+    binsentry_.reset();
+  }
+}
+
+void CMSHistSum::runBarlowBeeston() const {
+  if (!bb_.init) return;
+  RooAbsArg::setDirtyInhibit(true);
+
+  const unsigned n = bb_.use.size();
+  for (unsigned j = 0; j < n; ++j) {
+    bb_.dat[j] = data_[bb_.use[j]];
+    bb_.valsum[j] = valsum_[bb_.use[j]] * cache_.GetWidth(bb_.use[j]);
+    bb_.toterr[j] = toterr_[bb_.use[j]] * cache_.GetWidth(bb_.use[j]);
+  }
+  // This pragma statement tells (modern) gcc that loop can be safely
+  // vectorized
+  #pragma GCC ivdep
+  for (unsigned j = 0; j < n; ++j) {
+    bb_.b[j] = bb_.toterr[j] + (bb_.valsum[j] / bb_.toterr[j]) - bb_.gobs[j];
+    bb_.c[j] = bb_.valsum[j] - bb_.dat[j] - (bb_.valsum[j] / bb_.toterr[j]) * bb_.gobs[j];
+    bb_.tmp[j] = -0.5 * (bb_.b[j] + copysign(1.0, bb_.b[j]) * std::sqrt(bb_.b[j] * bb_.b[j] - 4. * bb_.c[j]));
+    bb_.x1[j] = bb_.tmp[j];
+    bb_.x2[j] = bb_.c[j] / bb_.tmp[j];
+    bb_.res[j] = std::max(bb_.x1[j], bb_.x2[j]);
+  }
+  for (unsigned j = 0; j < n; ++j) {
+    if (toterr_[bb_.use[j]] > 0.) bb_.push_res[j]->setVal(bb_.res[j]);
+  }
+  RooAbsArg::setDirtyInhibit(false);
+  for (RooAbsArg *arg : bb_.dirty_prop) {
+    arg->setValueDirty();
+  }
+}
+
+void CMSHistSum::setAnalyticBarlowBeeston(bool flag) const {
+  // Clear it if it's already initialised
+  if (bb_.init && flag) return;
+  if (bb_.init && !flag) {
+    for (unsigned i = 0; i < bb_.push_res.size(); ++i) {
+      bb_.push_res[i]->setConstant(false);
+    }
+    bb_.use.clear();
+    bb_.dat.clear();
+    bb_.valsum.clear();
+    bb_.toterr.clear();
+    bb_.err.clear();
+    bb_.b.clear();
+    bb_.c.clear();
+    bb_.tmp.clear();
+    bb_.x1.clear();
+    bb_.x2.clear();
+    bb_.res.clear();
+    bb_.gobs.clear();
+    bb_.dirty_prop.clear();
+    bb_.push_res.clear();
+    bb_.init = false;
+  }
+  if (flag && data_.size()) {
+    for (unsigned j = 0; j < bintypes_.size(); ++j) {
+      if (bintypes_[j][0] == 1 && !vbinpars_[j][0]->isConstant()) {
+        bb_.use.push_back(j);
+        double gobs_val = 0.;
+        RooFIter iter = vbinpars_[j][0]->valueClientMIterator();
+        RooAbsArg *arg = nullptr;
+        while((arg = iter.next())) {
+          if (arg == this || arg == &binsentry_) {
+            // std::cout << "Skipping " << this << " " << this->GetName() << "\n";
+          } else {
+            // std::cout << "Adding " << arg << " " << arg->GetName() << "\n";
+            bb_.dirty_prop.insert(arg);
+            auto as_gauss = dynamic_cast<RooGaussian*>(arg);
+            if (as_gauss) {
+              auto gobs = dynamic_cast<RooAbsReal*>(as_gauss->findServer(TString(vbinpars_[j][0]->GetName())+"_In"));
+              if (gobs) gobs_val = gobs->getVal();
+            }
+          }
+        }
+        bb_.gobs.push_back(gobs_val);
+        bb_.push_res.push_back((RooRealVar*)vbinpars_[j][0]);
+        bb_.push_res.back()->setConstant(true);
+      }
+    }
+    unsigned n = bb_.use.size();
+    bb_.dat.resize(n);
+    bb_.valsum.resize(n);
+    bb_.toterr.resize(n);
+    bb_.err.resize(n);
+    bb_.b.resize(n);
+    bb_.c.resize(n);
+    bb_.tmp.resize(n);
+    bb_.x1.resize(n);
+    bb_.x2.resize(n);
+    bb_.res.resize(n);
+    bb_.init = true;
+  }
+}
+
+
+RooArgList * CMSHistSum::setupBinPars(double poissonThreshold) {
+  RooArgList * res = new RooArgList();
+  if (bintypes_.size()) {
+    std::cout << "setupBinPars() already called for " << this->GetName() << "\n";
+    return res;
+  }
+
+  // First initialize all the storage
+  initialize();
+  // Now fill the bin contents and errors
+  updateCache(); // the arg (1) forces updateCache to fill the caches for all bins
+
+  bintypes_.resize(valsum_.size(), std::vector<unsigned>(1, 0));
+
+
+  std::cout << std::string(60, '=') << "\n";
+  std::cout << "Analysing bin errors for: " << this->GetName() << "\n";
+  std::cout << "Poisson cut-off: " << poissonThreshold << "\n";
+  std::set<unsigned> skip_idx;
+  std::vector<std::string> skipped_procs;
+  for (unsigned i = 0; i < vfuncstmp_.size(); ++i) {
+    vfuncstmp_[i]->updateCache();
+    if (vfuncstmp_[i]->attributes().count("skipForErrorSum")) {
+      skipped_procs.push_back(vfuncstmp_[i]->getStringAttribute("combine.process"));
+      skip_idx.insert(i);
+    }
+  }
+  if (skipped_procs.size()) {
+    std::cout << "Processes excluded for sums:";
+    for (auto &s: skipped_procs) std::cout << " " << s;
+    std::cout << "\n";
+  }
+  std::cout << std::string(60, '=') << "\n";
+  std::cout << TString::Format("%-10s %-15s %-15s %-30s\n", "Bin", "Contents", "Error", "Notes");
+
+  for (unsigned j = 0; j < valsum_.size(); ++j) {
+    std::cout << TString::Format("%-10i %-15f %-15f %-30s\n", j, valsum_[j], toterr_[j], "total sum");
+    double sub_sum = 0.;
+    double sub_err = 0.;
+    // Check using a possible sub-set of bins
+    for (unsigned i = 0; i < vfuncstmp_.size(); ++i) {
+      if (skip_idx.count(i)) {
+        continue;
+      }
+      sub_sum += vfuncstmp_[i]->cache()[j] * coeffvals_[i];
+      sub_err += std::pow(vfuncstmp_[i]->errors()[j] * coeffvals_[i], 2.);;
+    }
+    sub_err = std::sqrt(sub_err);
+    if (skipped_procs.size()) {
+      std::cout << TString::Format("%-10i %-15f %-15f %-30s\n", j, sub_sum, sub_err, "excluding marked processes");
+    }
+
+    if (sub_err <= 0.) {
+      std::cout << TString::Format("  %-30s\n", "=> Error is zero, ignore");
+      std::cout << std::string(60, '-') << "\n";
+      continue;
+    }
+
+    // Now check if we are below the poisson threshold
+    double n = int(0.5 + ((sub_sum * sub_sum) / (sub_err * sub_err)));
+    double alpha = valsum_[j] / n;
+    std::cout << TString::Format(
+        "%-10i %-15f %-15f %-30s\n", j, n, std::sqrt(n),
+        TString::Format("Unweighted events, alpha=%f", alpha).Data());
+
+    if (n <= poissonThreshold) {
+      std::cout << TString::Format("  %-30s\n", "=> Number of weighted events is below poisson threshold");
+
+      bintypes_[j].resize(vfuncstmp_.size(), 4);
+
+      for (unsigned i = 0; i < vfuncstmp_.size(); ++i) {
+        std::string proc =
+            vfuncstmp_[i]->stringAttributes().count("combine.process")
+                ? vfuncstmp_[i]->getStringAttribute("combine.process")
+                : vfuncstmp_[i]->GetName();
+        double v_p = vfuncstmp_[i]->cache()[j];
+        double e_p = vfuncstmp_[i]->errors()[j];
+        std::cout << TString::Format("    %-20s %-15f %-15f %-30s\n", proc.c_str(), v_p, e_p, "");
+        // relax the condition of v_p >= e_p slightly due to numerical rounding...
+        // Possibilities:
+        //    v_p = any   e_p <= 0       : skip
+        //    v_p < 0     e_p > 0        : for now, skip but technically we should be able to handle this in the future
+        //    v_p >= 0    e_p > v_p      : Create an additive gaussian constraint for this bin
+        //    v_p > 0     0 < e_p <= v_p : do the poisson
+        if (e_p <= 0.) {
+          std::cout << TString::Format("      %-30s\n", "=> Error is zero, ignore");
+          bintypes_[j][i] = 4;
+        } else if (v_p < 0. && e_p > 0.) {
+          std::cout << TString::Format("      %-30s\n", "=> Cannot handle negative content, ignore");
+          bintypes_[j][i] = 4;
+        } else if (v_p > 0. && e_p > 0. && v_p >= (e_p*0.999)) {
+          double n_p_r = int(0.5 + ((v_p * v_p) / (e_p * e_p)));
+          double alpha_p_r = v_p / n_p_r;
+          std::cout << TString::Format(
+              "    %-20s %-15f %-15f %-30s\n", "", n_p_r, std::sqrt(n_p_r),
+              TString::Format("Unweighted events, alpha=%f", alpha_p_r).Data());
+          if (n_p_r <= poissonThreshold) {
+            double sigma = 7.;
+            double rmin = 0.5*ROOT::Math::chisquared_quantile(ROOT::Math::normal_cdf_c(sigma), n_p_r * 2.);
+            double rmax = 0.5*ROOT::Math::chisquared_quantile(1. - ROOT::Math::normal_cdf_c(sigma), n_p_r * 2. + 2.);
+            RooRealVar *var = new RooRealVar(TString::Format("%s_bin%i_%s", this->GetName(), j, proc.c_str()), "", n_p_r, rmin, rmax);
+            RooConstVar *cvar = new RooConstVar(TString::Format("%g", 1. / n_p_r), "", 1. / n_p_r);
+            RooProduct *prod = new RooProduct(TString::Format("%s_prod", var->GetName()), "", RooArgList(*var, *cvar));
+            var->addOwnedComponents(RooArgSet(*prod, *cvar));
+            var->setAttribute("createPoissonConstraint");
+            res->addOwned(*var);
+            binpars_.add(*prod);
+
+            std::cout << TString::Format(
+                "      => Product of %s[%.2f,%.2f,%.2f] and const [%.4f] to be poisson constrained\n",
+                var->GetName(), var->getVal(), var->getMin(), var->getMax(), cvar->getVal());
+            bintypes_[j][i] = 2;
+          } else {
+            RooRealVar *var = new RooRealVar(TString::Format("%s_bin%i_%s", this->GetName(), j, proc.c_str()), "", 0, -7, 7);
+            std::cout << TString::Format(
+                "      => Parameter %s[%.2f,%.2f,%.2f] to be gaussian constrained\n",
+                var->GetName(), var->getVal(), var->getMin(), var->getMax());
+            var->setAttribute("createGaussianConstraint");
+            res->addOwned(*var);
+            binpars_.add(*var);
+            bintypes_[j][i] = 3;
+          }
+        } else if (v_p >= 0 && e_p > v_p) {
+          RooRealVar *var = new RooRealVar(TString::Format("%s_bin%i_%s", this->GetName(), j, proc.c_str()), "", 0, -7, 7);
+          std::cout << TString::Format(
+              "      => Poisson not viable, %s[%.2f,%.2f,%.2f] to be gaussian constrained\n",
+              var->GetName(), var->getVal(), var->getMin(), var->getMax());
+          var->setAttribute("createGaussianConstraint");
+          res->addOwned(*var);
+          binpars_.add(*var);
+          bintypes_[j][i] = 3;
+        } else{
+          std::cout << "      => ERROR: shouldn't be here\n";
+        }
+        std::cout << "  " << std::string(58, '-') << "\n";
+
+      }
+    } else if (toterr_[j] > 0.) {
+      bintypes_[j][0] = 1;
+      RooRealVar *var = new RooRealVar(TString::Format("%s_bin%i", this->GetName(), j), "", 0, -7, 7);
+      std::cout << TString::Format(
+          "  => Total parameter %s[%.2f,%.2f,%.2f] to be gaussian constrained\n",
+          var->GetName(), var->getVal(), var->getMin(), var->getMax());
+      var->setAttribute("createGaussianConstraint");
+      var->setAttribute("forBarlowBeeston");
+      res->addOwned(*var);
+      binpars_.add(*var);
+    }
+    std::cout << std::string(60, '-') << "\n";
+  }
+
+  // binpars_.add(*res);
+  binsentry_.addVars(binpars_);
+  binsentry_.setValueDirty();
+
+  for (unsigned j = 0, r = 0; j < valsum_.size(); ++j) {
+    vbinpars_[j].resize(bintypes_[j].size());
+    for (unsigned i = 0; i < bintypes_[j].size(); ++i) {
+      if (bintypes_[j][i] >= 1 && bintypes_[j][i] < 4) {
+        vbinpars_[j][i] = dynamic_cast<RooAbsReal *>(binpars_.at(r));
+        ++r;
+      }
+    }
+  }
+  return res;
+}
+
+
+std::unique_ptr<RooArgSet> CMSHistSum::getSentryArgs() const {
+  // We can do this without initialising because we're going to hand over
+  // the sentry directly
+  sentry_.SetName(TString(this->GetName()) + "_sentry");
+  binsentry_.SetName(TString(this->GetName()) + "_binsentry");
+  std::unique_ptr<RooArgSet> args(new RooArgSet(sentry_, binsentry_));
+  return args;
+}
+
+Double_t CMSHistSum::evaluate() const {
+  updateCache();
+  return cache().GetAt(x_);
+}
+
+std::map<std::string, Double_t> CMSHistSum::getProcessNorms() const {
+  std::map<std::string, Double_t> vals_;
+  initialize();
+  if (!sentry_.good()) {
+    updateMorphs();
+
+    for (unsigned i = 0; i < vcoeffpars_.size(); ++i) {
+      Double_t coeffval = vcoeffpars_[i]->getVal();    
+
+      staging_ = compcache_[i];
+      if (vtype_[i] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
+        staging_.Exp();
+      }
+      staging_.CropUnderflows();
+
+      double const * valarray = &(staging_[0]);
+      Double_t valsum = 0;
+      for(unsigned j=0; j < valsum_.size(); ++j ){
+        valsum += coeffval * valarray[j];
+      }
+
+      vals_[vcoeffpars_[i]->GetName()] = valsum;
+    }
+  }
+  return vals_;
+}
+
+void CMSHistSum::printMultiline(std::ostream& os, Int_t contents, Bool_t verbose,
+                    TString indent) const {
+  RooAbsReal::printMultiline(os, contents, verbose, indent);
+  updateCache();
+  if (bintypes_.size()) {
+    std::cout << ">> Bin types are:\n";
+    for (unsigned j = 0; j < bintypes_.size(); ++j) {
+      std::cout << ">> " << j << ":";
+      for (unsigned i = 0; i < bintypes_[j].size(); ++i) {
+        std::cout << " " << bintypes_[j][i];
+      }
+      std::cout << "\n";
+    }
+  }
+  std::cout << ">> Current cache:\n";
+  valsum_.Dump();
+  std::cout << ">> Current cache (bin scaled):\n";
+  cache_.Dump();
+  std::cout << ">> Sentry: " << sentry_.good() << "\n";
+  sentry_.Print("v");
+
+}
+
+Int_t CMSHistSum::getAnalyticalIntegral(RooArgSet& allVars,
+                                         RooArgSet& analVars,
+                                         const char* /*rangeName*/) const {
+  if (matchArgs(allVars, analVars, x_)) return 1;
+  return 0;
+}
+
+Double_t CMSHistSum::analyticalIntegral(Int_t code,
+                                         const char* rangeName) const {
+  // TODO: check how RooHistFunc handles ranges that splice bins
+  switch (code) {
+    case 1: {
+      updateCache();
+      return cache().IntegralWidth();
+    }
+  }
+
+  assert(0);
+  return 0;
+}
+
+void CMSHistSum::setData(RooAbsData const& data) const {
+  updateCache();
+  data_.clear();
+  data_.resize(cache_.fullsize(), 0.); // fullsize is important here is we've used activeBins
+  RooArgSet obs(x_.arg());
+  const RooRealVar& x = static_cast<const RooRealVar&>(*obs.first());
+  for (int i = 0, n = data.numEntries(); i < n; ++i) {
+    obs = *data.get(i);
+    int idx = cache_.FindBin(x.getVal());
+    data_[idx] = data.weight();
+  }
+}
+
+void CMSHistSum::EnableFastVertical() {
+  enable_fast_vertical_ = true;
+}
+
+#undef HFVERBOSE
+
diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc
index de6c051a3f9..171eaa581ba 100644
--- a/src/CachingNLL.cc
+++ b/src/CachingNLL.cc
@@ -13,6 +13,7 @@
 #include <HiggsAnalysis/CombinedLimit/interface/CMSHistV.h>
 #include <HiggsAnalysis/CombinedLimit/interface/CMSHistFunc.h>
 #include <HiggsAnalysis/CombinedLimit/interface/CMSHistErrorPropagator.h>
+#include <HiggsAnalysis/CombinedLimit/interface/CMSHistSum.h>
 #include <HiggsAnalysis/CombinedLimit/interface/CMSHistFuncWrapper.h>
 #include <HiggsAnalysis/CombinedLimit/interface/VectorizedGaussian.h>
 #include <HiggsAnalysis/CombinedLimit/interface/VectorizedCB.h>
@@ -30,6 +31,7 @@ namespace cacheutils {
     typedef OptimizedCachingPdfT<CMSHistFunc, CMSHistV<CMSHistFunc>> CachingCMSHistFunc;
     typedef OptimizedCachingPdfT<CMSHistFuncWrapper, CMSHistV<CMSHistFuncWrapper>> CachingCMSHistFuncWrapper;
     typedef OptimizedCachingPdfT<CMSHistErrorPropagator, CMSHistV<CMSHistErrorPropagator>> CachingCMSHistErrorPropagator;
+    typedef OptimizedCachingPdfT<CMSHistSum, CMSHistV<CMSHistSum>> CachingCMSHistSum;
     typedef OptimizedCachingPdfT<RooGaussian,VectorizedGaussian> CachingGaussPdf;
     typedef OptimizedCachingPdfT<RooCBShape,VectorizedCBShape> CachingCBPdf;
     typedef OptimizedCachingPdfT<RooExponential,VectorizedExponential> CachingExpoPdf;
@@ -248,33 +250,41 @@ std::pair<std::vector<Double_t> *, bool> cacheutils::ValuesCache::get()
     return std::pair<std::vector<Double_t> *, bool>(&items[found]->values, good);
 }
 
-cacheutils::CachingPdf::CachingPdf(RooAbsReal *pdf, const RooArgSet *obs) :
-    obs_(obs),
-    pdfOriginal_(pdf),
-    pdfPieces_(),
-    pdf_(runtimedef::get("CACHINGPDF_NOCLONE") ? pdfOriginal_ : (runtimedef::get("CACHINGPDF_NOCHEAPCLONE") ? utils::fullCloneFunc(pdfOriginal_, pdfPieces_) : utils::fullCloneFunc(pdfOriginal_, *obs_, pdfPieces_))),
-    lastData_(0),
-    cache_(*pdf_,*obs_),
-    includeZeroWeights_(false)
-{
-    if (runtimedef::get("CACHINGPDF_DIRECT") || pdf->getAttribute("CachingPdf_Direct")) {
-        cache_.setDirectMode(true);
-    }
+cacheutils::CachingPdf::CachingPdf(RooAbsReal *pdf, const RooArgSet *obs)
+    : obs_(obs),
+      pdfOriginal_(pdf),
+      pdfPieces_(),
+      pdf_((runtimedef::get("CACHINGPDF_NOCLONE") &&
+            pdfOriginal_->getAttribute("CachingPdf_NoClone"))
+               ? pdfOriginal_
+               : (runtimedef::get("CACHINGPDF_NOCHEAPCLONE")
+                      ? utils::fullCloneFunc(pdfOriginal_, pdfPieces_)
+                      : utils::fullCloneFunc(pdfOriginal_, *obs_, pdfPieces_))),
+      lastData_(0),
+      cache_(*pdf_, *obs_),
+      includeZeroWeights_(false) {
+  if (runtimedef::get("CACHINGPDF_DIRECT") || pdf->getAttribute("CachingPdf_Direct")) {
+    cache_.setDirectMode(true);
+  }
 }
 
-cacheutils::CachingPdf::CachingPdf(const CachingPdf &other) :
-    obs_(other.obs_),
-    pdfOriginal_(other.pdfOriginal_),
-    pdfPieces_(),
-    pdf_(runtimedef::get("CACHINGPDF_NOCLONE") ? pdfOriginal_ : (runtimedef::get("CACHINGPDF_NOCHEAPCLONE") ? utils::fullCloneFunc(pdfOriginal_, pdfPieces_) : utils::fullCloneFunc(pdfOriginal_, *obs_, pdfPieces_))),
-    lastData_(0),
-    cache_(*pdf_,*obs_),
-    includeZeroWeights_(other.includeZeroWeights_)
-{
-    if (runtimedef::get("CACHINGPDF_DIRECT") || other.pdfOriginal_->getAttribute("CachingPdf_Direct")) {
-        cache_.setDirectMode(true);
-    }
-
+cacheutils::CachingPdf::CachingPdf(const CachingPdf &other)
+    : obs_(other.obs_),
+      pdfOriginal_(other.pdfOriginal_),
+      pdfPieces_(),
+      pdf_((runtimedef::get("CACHINGPDF_NOCLONE") &&
+            other.pdfOriginal_->getAttribute("CachingPdf_NoClone"))
+               ? pdfOriginal_
+               : (runtimedef::get("CACHINGPDF_NOCHEAPCLONE")
+                      ? utils::fullCloneFunc(pdfOriginal_, pdfPieces_)
+                      : utils::fullCloneFunc(pdfOriginal_, *obs_, pdfPieces_))),
+      lastData_(0),
+      cache_(*pdf_, *obs_),
+      includeZeroWeights_(other.includeZeroWeights_) {
+  if (runtimedef::get("CACHINGPDF_DIRECT") ||
+      other.pdfOriginal_->getAttribute("CachingPdf_Direct")) {
+    cache_.setDirectMode(true);
+  }
 }
 
 cacheutils::CachingPdf::~CachingPdf() 
@@ -507,6 +517,8 @@ cacheutils::makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) {
         return new CachingCMSHistFuncWrapper(pdf, obs);
     } else if (histfuncNll && typeid(*pdf) == typeid(CMSHistErrorPropagator)) {
         return new CachingCMSHistErrorPropagator(pdf, obs);
+    } else if (histfuncNll && typeid(*pdf) == typeid(CMSHistSum)) {
+        return new CachingCMSHistSum(pdf, obs);
     } else {
         if (verb) {
             std::cout << "I don't have an optimized implementation for " << pdf->ClassName() << " (" << pdf->GetName() << ")" << std::endl;
@@ -806,6 +818,10 @@ void cacheutils::CachingAddNLL::propagateData() {
             // printf("Passing data to %s\n", funci.pdf()->GetName());
             (static_cast<CMSHistErrorPropagator const*>(funci.pdf()))->setData(*data_);
         }
+        if (typeid(*(funci.pdf())) == typeid(CMSHistSum)) {
+            // printf("Passing data to %s\n", funci.pdf()->GetName());
+            (static_cast<CMSHistSum const*>(funci.pdf()))->setData(*data_);
+        }
     }
 }
 
@@ -815,6 +831,10 @@ void cacheutils::CachingAddNLL::setAnalyticBarlowBeeston(bool flag) {
         if (typeid(*(funci.pdf())) == typeid(CMSHistErrorPropagator)) {
             (static_cast<CMSHistErrorPropagator const*>(funci.pdf()))->setAnalyticBarlowBeeston(flag);
         }
+        if (typeid(*(funci.pdf())) == typeid(CMSHistSum)) {
+            (static_cast<CMSHistSum const*>(funci.pdf()))->setAnalyticBarlowBeeston(flag);
+        }
+
     }
 }
 
diff --git a/src/Combine.cc b/src/Combine.cc
index 163bd820ede..101fb7562ba 100644
--- a/src/Combine.cc
+++ b/src/Combine.cc
@@ -63,6 +63,7 @@
 #include "HiggsAnalysis/CombinedLimit/interface/ProfilingTools.h"
 #include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistFunc.h"
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHistSum.h"
 
 #include "HiggsAnalysis/CombinedLimit/interface/Logger.h"
 
@@ -897,6 +898,7 @@ void Combine::run(TString hlfFile, const std::string &dataset, double &limit, do
 
   if (runtimedef::get("FAST_VERTICAL_MORPH")) {
     CMSHistFunc::EnableFastVertical();
+    CMSHistSum::EnableFastVertical();
   }
 
   // Warn the user that they might be using funky values of POIs 
diff --git a/src/ProcessNormalization.cc b/src/ProcessNormalization.cc
index 1f770ea6075..0216cb9f9ce 100644
--- a/src/ProcessNormalization.cc
+++ b/src/ProcessNormalization.cc
@@ -58,26 +58,51 @@ void ProcessNormalization::addOtherFactor(RooAbsReal &factor) {
 
 Double_t ProcessNormalization::evaluate() const {
     double logVal = 0.0;
+    if (thetaListVec_.empty()) {
+        RooFIter iterTheta = thetaList_.fwdIterator();
+        std::vector<RooAbsReal *> & thetaListVec = const_cast<std::vector<RooAbsReal *>&>(thetaListVec_);
+        thetaListVec.reserve(thetaList_.getSize());
+        for (RooAbsArg *a = iterTheta.next(); a != 0; a = iterTheta.next()) {
+            thetaListVec.push_back(dynamic_cast<RooAbsReal *>(a));
+        }
+    }
+    if (asymmThetaListVec_.empty()) {
+        RooFIter iterTheta = asymmThetaList_.fwdIterator();
+        std::vector<RooAbsReal *> & asymmThetaListVec = const_cast<std::vector<RooAbsReal *>&>(asymmThetaListVec_);
+        asymmThetaListVec.reserve(asymmThetaList_.getSize());
+        for (RooAbsArg *a = iterTheta.next(); a != 0; a = iterTheta.next()) {
+            asymmThetaListVec.push_back(dynamic_cast<RooAbsReal *>(a));
+        }
+    }
+    if (otherFactorListVec_.empty()) {
+        RooFIter iterOther = otherFactorList_.fwdIterator();
+        std::vector<RooAbsReal *> & otherFactorListVec = const_cast<std::vector<RooAbsReal *>&>(otherFactorListVec_);
+        otherFactorListVec.reserve(otherFactorList_.getSize());
+        for (RooAbsArg *a = iterOther.next(); a != 0; a = iterOther.next()) {
+            otherFactorListVec.push_back(dynamic_cast<RooAbsReal *>(a));
+        }
+    }
     if (!logKappa_.empty()) {
-        RooLinkedListIter iterTheta = thetaList_.iterator();
-        std::vector<double>::const_iterator logKappa = logKappa_.begin();
-        for (RooAbsReal *theta = (RooAbsReal*) iterTheta.Next(); theta != 0; theta = (RooAbsReal*) iterTheta.Next(), ++logKappa) {
-            logVal += theta->getVal() * (*logKappa);
+        assert(logKappa_.size()==thetaListVec_.size());
+        for(unsigned int i=0; i < thetaListVec_.size() ; i++){
+            const RooAbsReal *theta = thetaListVec_.at(i);
+            const double logKappa = logKappa_.at(i);
+            logVal += theta->getVal() * (logKappa);
         }
     }
     if (!logAsymmKappa_.empty()) {
-        RooLinkedListIter iterTheta = asymmThetaList_.iterator();
-        std::vector<std::pair<double,double> >::const_iterator logKappas = logAsymmKappa_.begin();
-        for (RooAbsReal *theta = (RooAbsReal*) iterTheta.Next(); theta != 0; theta = (RooAbsReal*) iterTheta.Next(), ++logKappas) {
+        assert(logAsymmKappa_.size()==asymmThetaListVec_.size());
+        for( unsigned int i=0; i < asymmThetaListVec_.size(); i++){
+            const RooAbsReal *theta = asymmThetaListVec_.at(i);
+            const std::pair<double,double> logKappas = logAsymmKappa_.at(i);
             double x = theta->getVal();
-            logVal +=  x * logKappaForX(x, *logKappas);
+            logVal +=  x * logKappaForX(x, logKappas);
         }
     }
     double norm = nominalValue_;
     if (logVal) norm *= std::exp(logVal);
     if (otherFactorList_.getSize()) {
-        RooLinkedListIter iterOther = otherFactorList_.iterator();
-        for (RooAbsReal *fact = (RooAbsReal*) iterOther.Next(); fact != 0; fact = (RooAbsReal*) iterOther.Next()) {
+        for (const RooAbsReal *fact :otherFactorListVec_){
             norm *= fact->getVal();
         }
     }
diff --git a/src/VerticalInterpPdf.cc b/src/VerticalInterpPdf.cc
index 66aff0a30d5..6c1b9c50354 100644
--- a/src/VerticalInterpPdf.cc
+++ b/src/VerticalInterpPdf.cc
@@ -27,6 +27,8 @@ VerticalInterpPdf::VerticalInterpPdf()
   _funcIter  = _funcList.createIterator() ;
   _coefIter  = _coefList.createIterator() ;
   _quadraticRegion = 0;
+  _pdfFloorVal = 1e-15;
+  _integralFloorVal = 1e-10;
 }
 
 
@@ -37,7 +39,9 @@ VerticalInterpPdf::VerticalInterpPdf(const char *name, const char *title, const
   _funcList("!funcList","List of functions",this),
   _coefList("!coefList","List of coefficients",this),
   _quadraticRegion(quadraticRegion),
-  _quadraticAlgo(quadraticAlgo)
+  _quadraticAlgo(quadraticAlgo),
+  _pdfFloorVal(1e-15),
+  _integralFloorVal(1e-10)
 { 
 
   if (inFuncList.getSize()!=2*inCoefList.getSize()+1) {
@@ -90,7 +94,9 @@ VerticalInterpPdf::VerticalInterpPdf(const VerticalInterpPdf& other, const char*
   _funcList("!funcList",this,other._funcList),
   _coefList("!coefList",this,other._coefList),
   _quadraticRegion(other._quadraticRegion),
-  _quadraticAlgo(other._quadraticAlgo)
+  _quadraticAlgo(other._quadraticAlgo),
+  _pdfFloorVal(other._pdfFloorVal),
+  _integralFloorVal(other._integralFloorVal)
 {
   // Copy constructor
 
@@ -141,7 +147,7 @@ Double_t VerticalInterpPdf::evaluate() const
       }
   }
    
-  return ( value > 0 ? value : 1E-15 ); // Last IEEE double precision
+  return ( value > 0. ? value : _pdfFloorVal);
 }
 
 
@@ -340,8 +346,8 @@ Double_t VerticalInterpPdf::analyticalIntegralWN(Int_t code, const RooArgSet* no
   }
 
   Double_t result = 0;
-  if(normVal>0) result = value / normVal;
-  return (result > 0 ? result : 1E-10);
+  if(normVal>0.) result = value / normVal;
+  return (result > 0. ? result : _integralFloorVal);
 }
 
 Double_t VerticalInterpPdf::interpolate(Double_t coeff, Double_t central, RooAbsReal *fUp, RooAbsReal *fDn) const  
@@ -420,3 +426,7 @@ Double_t VerticalInterpPdf::interpolate(Double_t coeff, Double_t central, RooAbs
         }
     }
 }
+
+
+//_____________________________________________________________________________
+void VerticalInterpPdf::setFloorVals(Double_t const& pdf_val, Double_t const& integral_val){ _pdfFloorVal = pdf_val; _integralFloorVal = integral_val; }
diff --git a/src/classes.h b/src/classes.h
index eb3a0dea4e3..c4aca7d6fdc 100644
--- a/src/classes.h
+++ b/src/classes.h
@@ -52,6 +52,7 @@
 #include "HiggsAnalysis/CombinedLimit/interface/RooDoubleCBFast.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistFunc.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistErrorPropagator.h"
+#include "HiggsAnalysis/CombinedLimit/interface/CMSHistSum.h"
 #include "HiggsAnalysis/CombinedLimit/interface/CMSHistFuncWrapper.h"
 #include "HiggsAnalysis/CombinedLimit/interface/RooTaylorExpansion.h"
 #include "HiggsAnalysis/CombinedLimit/interface/SimpleTaylorExpansion1D.h"
diff --git a/src/classes_def.xml b/src/classes_def.xml
index 447a443a6e4..98dc9faddd4 100644
--- a/src/classes_def.xml
+++ b/src/classes_def.xml
@@ -202,6 +202,7 @@
         <function name="function th1fmorph" />
   <class name="CMSHistFunc" />
   <class name="CMSHistErrorPropagator" />
+  <class name="CMSHistSum" />
   <class name="CMSHistFuncWrapper" />
   <class name="SimpleTaylorExpansion1D" />
   <class name="RooTaylorExpansion" />
diff --git a/test/diffNuisances.py b/test/diffNuisances.py
index ecadfe3611f..903293a025c 100644
--- a/test/diffNuisances.py
+++ b/test/diffNuisances.py
@@ -54,9 +54,9 @@
 
 if options.sortBy not in ['correlation','impact']: exit("choose one of [ %s ] for --sortBy"%(",".join()['correlation','impact']))
 
-if args.regex != ".*":
+if options.regex != ".*":
     print("Including only nuisance parameters following this regex query:")
-    print(args.regex)
+    print(options.regex)
 
 setUpString = "diffNuisances run on %s, at %s with the following options ... "%(args[0],datetime.datetime.utcnow())+str(options)
 
@@ -97,19 +97,25 @@ def getGraph(hist,shift):
    return gr
 """
 
+# Compile regex object, since we will be checking regex patterns several times.
+regex_NP_obj = re.compile(options.regex)
+
+np_count = 0
+for i in range(fpf_s.getSize()):
+    nuis_s = fpf_s.at(i)
+    name = nuis_s.GetName();
+    if bool(regex_NP_obj.match(name)): np_count += 1
+
 # Also make histograms for pull distributions:
-hist_fit_b  = ROOT.TH1F("fit_b"   ,"B-only fit Nuisances;;%s "%title,prefit.getSize(),0,prefit.getSize())
-hist_fit_s  = ROOT.TH1F("fit_s"   ,"S+B fit Nuisances   ;;%s "%title,prefit.getSize(),0,prefit.getSize())
-hist_prefit = ROOT.TH1F("prefit_nuisancs","Prefit Nuisances    ;;%s "%title,prefit.getSize(),0,prefit.getSize())
+hist_fit_b  = ROOT.TH1F("fit_b"   ,"B-only fit Nuisances;;%s "%title,np_count,0,np_count)
+hist_fit_s  = ROOT.TH1F("fit_s"   ,"S+B fit Nuisances   ;;%s "%title,np_count,0,np_count)
+hist_prefit = ROOT.TH1F("prefit_nuisancs","Prefit Nuisances    ;;%s "%title,np_count,0,np_count)
 # Store also the *asymmetric* uncertainties
 gr_fit_b    = ROOT.TGraphAsymmErrors(); gr_fit_b.SetTitle("fit_b_g")
 gr_fit_s    = ROOT.TGraphAsymmErrors(); gr_fit_s.SetTitle("fit_b_s")
 
 error_poi = fpf_s.find(options.poi).getError()
 
-# Compile regex object, since we will be checking regex patterns several times.
-regex_NP_obj = re.compile(args.regex)
-
 # loop over all fitted parameters
 for i in range(fpf_s.getSize()):
 
diff --git a/test/printWorkspaceNormalisations.py b/test/printWorkspaceNormalisations.py
index 9be554949bd..bde7d85b063 100644
--- a/test/printWorkspaceNormalisations.py
+++ b/test/printWorkspaceNormalisations.py
@@ -2,6 +2,7 @@
 from sys import argv, stdout, stderr, exit
 import datetime
 from optparse import OptionParser
+from collections import OrderedDict as od
 
 hasHelp = False
 for X in ("-h", "-?", "--help"):
@@ -19,6 +20,8 @@
 parser.add_option("", "--printValueOnly",  dest="printValueOnly", default=False, action='store_true', help="Just print the default value of the normalisation.")
 parser.add_option("", "--min_threshold",  dest="min_threshold", default=-1.0, type='float', help="Only print values if yield is greater than this threshold.")
 parser.add_option("", "--max_threshold",  dest="max_threshold", default=-1.0, type='float', help="Only print values if yield is less than this threshold.")
+parser.add_option("", "--use-cms-histsum",  dest="use_cms_histsum", default=False, action="store_true", help="Set to true if workspace built with --use-cms-histsum option.")
+parser.add_option("", "--output-json",  dest="output_json", default='', help="Output norms in json.")
 
 
 (options, args) = parser.parse_args()
@@ -78,8 +81,23 @@ def find_chan_proc(name):
   chan_procs[chan].append([proc,norm,0,0])
  else: chan_procs[chan]= [[proc,norm,0,0]]
 
-# Now print out information 
+# Look for cases where chan stored as CMSHistSum, set flag
+if options.use_cms_histsum:
+  chan_CMSHistSum_norms = {}
+  all_props = ws.allFunctions().selectByName("prop_bin*")
+  for chan in chan_procs.keys():
+    prop_it = all_props.createIterator()
+    for i in range(all_props.getSize()):
+      prop = prop_it.Next()
+      prop_name = prop.GetName()
+      if chan == prop_name.split("_bin")[-1]:
+        if type(prop) == ROOT.CMSHistSum: chan_CMSHistSum_norms[chan] = dict(prop.getProcessNorms())
+
+# Now print out information
+default_norms = od()
+
 for chan in chan_procs.keys(): 
+ default_norms[chan] = od()
  print "---------------------------------------------------------------------------"
  print "---------------------------------------------------------------------------"
  print "Channel - %s "%chan 
@@ -92,7 +110,16 @@ def find_chan_proc(name):
   print "---------------------------------------------------------------------------"
   print "  Top-level normalisation for process %s -> %s"%(proc[0],proc[1].GetName())
   print "  -------------------------------------------------------------------------"
-  if options.printValueOnly: print "  default value = ",proc[1].getVal()
+  if options.use_cms_histsum:
+    if chan in chan_CMSHistSum_norms:
+      default_val = chan_CMSHistSum_norms[chan][proc[1].GetName()]
+    else:
+      default_val = proc[1].getVal()
+  else:
+    default_val = proc[1].getVal()
+  default_norms[chan][proc[1].GetName()] = default_val
+
+  if options.printValueOnly: print "  default value = ",default_val
   #if options.printValueOnly: print " --xcp %s:%s "%(chan,proc[0]),
   else: 
     if proc[2]: 
@@ -115,4 +142,25 @@ def find_chan_proc(name):
         proc[1].Print()
       	print " ... is a constant (RooRealVar)"
     print "  -------------------------------------------------------------------------"
-    print "  default value = ",proc[1].getVal()
+    print "  default value = ", default_val
+
+# Save norms to json file
+if options.output_json != '':
+  if ".json" not in options.output_json: options.output_json += ".json"
+  with open(options.output_json,"w") as jf:
+    jf.write("{\n")
+    chans = default_norms.keys()
+    for chan in chans:
+      jf.write("    \"%s\":{\n"%chan)
+      procs = default_norms[chan].keys()
+      for proc in default_norms[chan]:
+        if proc == procs[-1]:
+          jf.write("        \"%s\":%.8f\n"%(proc,default_norms[chan][proc]))
+        else:
+          jf.write("        \"%s\":%.8f,\n"%(proc,default_norms[chan][proc]))
+      if chan == chans[-1]:
+        jf.write("    }\n")
+      else:
+        jf.write("    },\n")
+    jf.write("}")
+    
diff --git a/test/systematicsAnalyzer.py b/test/systematicsAnalyzer.py
index 0e8fec4e9ca..07f8f87f8b5 100755
--- a/test/systematicsAnalyzer.py
+++ b/test/systematicsAnalyzer.py
@@ -43,6 +43,7 @@
 options.optimizeMHDependency = False 
 options.optimizeTemplateBins=False
 options.forceNonSimPdf = False
+options.removeMultiPdf = False
 options.physModel = "HiggsAnalysis.CombinedLimit.PhysicsModel:defaultModel"
 # import ROOT with a fix to get batch mode (http://root.cern.ch/phpBB3/viewtopic.php?t=3198)
 import sys