diff --git a/interface/CMSHistSum.h b/interface/CMSHistSum.h index 4fe3e375ac9..791f737431a 100644 --- a/interface/CMSHistSum.h +++ b/interface/CMSHistSum.h @@ -73,9 +73,14 @@ class CMSHistSum : public RooAbsReal { RooArgList const& coefList() const { return coeffpars_; } // RooArgList const& funcList() const { return funcs_; } + RooAbsReal const& getXVar() const { return x_.arg(); } + static void EnableFastVertical(); friend class CMSHistV; + // TODO: allow any class that implements hasChanged() and batchGetBinValues() + void injectExternalMorph(int idx, CMSInterferenceFunc& morph); + protected: RooRealProxy x_; @@ -128,6 +133,9 @@ class CMSHistSum : public RooAbsReal { mutable int fast_mode_; //! not to be serialized static bool enable_fast_vertical_; //! not to be serialized + RooListProxy external_morphs_; + std::vector external_morph_indices_; + inline int& morphField(int const& ip, int const& iv) { return vmorph_fields_[ip * n_morphs_ + iv]; } @@ -143,7 +151,7 @@ class CMSHistSum : public RooAbsReal { private: - ClassDef(CMSHistSum,1) + ClassDef(CMSHistSum,2) }; #endif diff --git a/python/InterferenceModels.py b/python/InterferenceModels.py index 9726da3edfd..88064370bde 100644 --- a/python/InterferenceModels.py +++ b/python/InterferenceModels.py @@ -139,13 +139,15 @@ def done(self): # if there are no systematics, it ends up with a different name? hfname = "shapeSig_{process}_{channel}_rebinPdf".format(channel=channel, process=process) histfunc = self.modelBuilder.out.function(hfname) + if not histfunc: + # assume this is a CMSHistSum workspace + hfname = f"prop_bin{channel}" + histfunc = self.modelBuilder.out.function(hfname) # TODO: support FastVerticalInterpHistPdf2 - if not isinstance(histfunc, ROOT.CMSHistFunc): - self.modelBuilder.out.Print("v") + if not isinstance(histfunc, (ROOT.CMSHistFunc, ROOT.CMSHistSum)): raise RuntimeError( - "Could not locate the CMSHistFunc for {string}.\nNote that CMSInterferenceFunc currently only supports workspaces that use CMSHistFunc".format( - string=string - ) + "Could not locate the CMSHistFunc or CMSHistSum for {string}.\n".format(string=string) + + "Note that CMSInterferenceFunc currently only supports workspaces that use these classes" ) funcname = hfname + "_externalMorph" @@ -164,7 +166,16 @@ def done(self): self.modelBuilder.out.safe_import(ROOT.CMSInterferenceFunc(funcname, "", histfunc.getXVar(), params, edges, scaling_array)) func = self.modelBuilder.out.function(funcname) - histfunc.injectExternalMorph(func) + if isinstance(histfunc, ROOT.CMSHistFunc): + histfunc.injectExternalMorph(func) + elif isinstance(histfunc, ROOT.CMSHistSum): + coefname = "n_exp_final_bin{channel}_proc_{process}".format(channel=channel, process=process) + for idx, coef in enumerate(histfunc.coefList()): + if coef.GetName() == coefname: + if self.verbose: + print("Injecting external morph for " + coefname) + histfunc.injectExternalMorph(idx, func) + break interferenceModel = InterferenceModel() diff --git a/src/CMSHistFunc.cc b/src/CMSHistFunc.cc index a05a5c1854f..2ef5504f082 100644 --- a/src/CMSHistFunc.cc +++ b/src/CMSHistFunc.cc @@ -323,7 +323,7 @@ void CMSHistFunc::updateCache() const { if (mcache_.size() == 0) mcache_.resize(storage_.size()); } - bool external_morph_updated = (external_morph_.getSize() && dynamic_cast(external_morph_.at(0))->hasChanged()); + bool external_morph_updated = (external_morph_.getSize() && static_cast(external_morph_.at(0))->hasChanged()); if (step1 || external_morph_updated) { fast_vertical_ = false; } @@ -567,7 +567,7 @@ void CMSHistFunc::updateCache() const { std::cout << "Template before external morph update:" << mcache_[idx].step2.Integral() << "\n"; mcache_[idx].step2.Dump(); #endif - auto& extdata = dynamic_cast(external_morph_.at(0))->batchGetBinValues(); + auto& extdata = static_cast(external_morph_.at(0))->batchGetBinValues(); for(size_t i=0; igetParameters({*x_}); + sentry_.addVars(*deps); + delete deps; + } + sentry_.setValueDirty(); binsentry_.setValueDirty(); @@ -250,6 +260,16 @@ void CMSHistSum::initialize() const { } void CMSHistSum::updateMorphs() const { + // set up pointers ahead of time for quick loop + std::vector process_morphs(compcache_.size(), nullptr); + // if any external morphs are dirty, disable fast_mode_ + for(size_t i=0; i < external_morph_indices_.size(); ++i) { + auto* morph = static_cast(external_morphs_.at(i)); + process_morphs[external_morph_indices_[i]] = morph; + if (morph->hasChanged()) { + fast_mode_ = 0; + } + } // If we're not in fast mode, need to reset all the compcache_ #if HFVERBOSE > 0 std::cout << "fast_mode_ = " << fast_mode_ << std::endl; @@ -257,6 +277,12 @@ void CMSHistSum::updateMorphs() const { for (unsigned ip = 0; ip < compcache_.size(); ++ip) { if (fast_mode_ == 0) { compcache_[ip].CopyValues(storage_[process_fields_[ip]]); + if ( process_morphs[ip] != nullptr ) { + auto& extdata = process_morphs[ip]->batchGetBinValues(); + for(size_t ibin=0; ibin= coeffpars_.getSize() ) { + throw std::runtime_error("Process index larger than number of processes in CMSHistSum"); + } + if ( morph.batchGetBinValues().size() != cache_.size() ) { + throw std::runtime_error("Mismatched binning between external morph and CMSHistSum"); + // equal edges are user responsibility for now + } + + for (auto other_idx : external_morph_indices_) { + if ( idx == other_idx ) { + external_morphs_.replace(external_morphs_[idx], morph); + return; + } + } + external_morph_indices_.push_back(idx); + external_morphs_.add(morph); +} + #undef HFVERBOSE diff --git a/test/test_interference.py b/test/test_interference.py index 9712293cfd2..1cb146909c5 100755 --- a/test/test_interference.py +++ b/test/test_interference.py @@ -147,6 +147,10 @@ def to_TH1(name, array): ret = subprocess.call(t2wcmd) assert ret == 0 +histsum_args = ["--for-fits", "--no-wrappers", "--use-histsum", "-o", "card_histsum.root"] +ret = subprocess.call(t2wcmd + histsum_args) +assert ret == 0 + fws = ROOT.TFile.Open("card.root") w = fws.Get("w") @@ -168,5 +172,27 @@ def setvars(x, kl, kv, k2v): setvars(2, 1.1, 0.9, 1.3) assert abs(func.getVal() - 4.372110974178483) < 1e-14, func.getVal() -ret = subprocess.call("combine -M MultiDimFit card.root -t 100".split(" ")) +# toy generation is different between the histsum and histfunc models, somehow +ntoys = 10 +ret = subprocess.call(f"combine -M GenerateOnly card.root -t {ntoys} --saveToys".split(" ")) assert ret == 0 + +ret = subprocess.call(f"combine -M MultiDimFit card.root -t {ntoys} --toysFile higgsCombineTest.GenerateOnly.mH120.123456.root -n HistFunc".split(" ")) +assert ret == 0 + +ret = subprocess.call(f"combine -M MultiDimFit card_histsum.root -t {ntoys} --toysFile higgsCombineTest.GenerateOnly.mH120.123456.root -n HistSum".split(" ")) +assert ret == 0 + +f_histfunc = ROOT.TFile.Open("higgsCombineHistFunc.MultiDimFit.mH120.123456.root") +f_histsum = ROOT.TFile.Open("higgsCombineHistSum.MultiDimFit.mH120.123456.root") + +ndiff = {"kl": 0, "kv": 0, "k2v": 0} +for row1, row2 in zip(f_histfunc.Get("limit"), f_histsum.Get("limit")): + if abs(row1.kl - row2.kl) > 1e-4: + ndiff["kl"] += 1 + if abs(row1.kv - row2.kv) > 1e-4: + ndiff["kv"] += 1 + if abs(row1.k2v - row2.k2v) > 1e-4: + ndiff["k2v"] += 1 + +print(f"Out of {ntoys} toys, {ndiff} are not matching (tolerance: 1e-4) between CMSHistFunc and CMSHistSum")