Skip to content

Commit

Permalink
Implement external morph for CMSHistSum as well
Browse files Browse the repository at this point in the history
  • Loading branch information
nsmith- committed Oct 31, 2023
1 parent e80071c commit 3e25193
Show file tree
Hide file tree
Showing 5 changed files with 102 additions and 12 deletions.
10 changes: 9 additions & 1 deletion interface/CMSHistSum.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<CMSHistSum>;

// TODO: allow any class that implements hasChanged() and batchGetBinValues()
void injectExternalMorph(int idx, CMSInterferenceFunc& morph);

protected:
RooRealProxy x_;

Expand Down Expand Up @@ -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<int> external_morph_indices_;

inline int& morphField(int const& ip, int const& iv) {
return vmorph_fields_[ip * n_morphs_ + iv];
}
Expand All @@ -143,7 +151,7 @@ class CMSHistSum : public RooAbsReal {


private:
ClassDef(CMSHistSum,1)
ClassDef(CMSHistSum,2)
};

#endif
23 changes: 17 additions & 6 deletions python/InterferenceModels.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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()
4 changes: 2 additions & 2 deletions src/CMSHistFunc.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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<CMSInterferenceFunc*>(external_morph_.at(0))->hasChanged());
bool external_morph_updated = (external_morph_.getSize() && static_cast<CMSInterferenceFunc*>(external_morph_.at(0))->hasChanged());
if (step1 || external_morph_updated) {
fast_vertical_ = false;
}
Expand Down Expand Up @@ -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<CMSInterferenceFunc*>(external_morph_.at(0))->batchGetBinValues();
auto& extdata = static_cast<CMSInterferenceFunc*>(external_morph_.at(0))->batchGetBinValues();
for(size_t i=0; i<extdata.size(); ++i) {
mcache_[idx].step2[i] *= extdata[i];
}
Expand Down
49 changes: 47 additions & 2 deletions src/CMSHistSum.cc
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ CMSHistSum::CMSHistSum(const char* name,
binsentry_(TString(name) + "_binsentry", ""),
initialized_(false),
analytic_bb_(false),
fast_mode_(0) {
fast_mode_(0),
external_morphs_("external_morphs", "", this) {

n_procs_ = funcs.getSize();
assert(n_procs_ == coeffs.getSize());
Expand Down Expand Up @@ -188,7 +189,10 @@ CMSHistSum::CMSHistSum(
sentry_(name ? TString(name) + "_sentry" : TString(other.GetName())+"_sentry", ""),
binsentry_(name ? TString(name) + "_binsentry" : TString(other.GetName())+"_binsentry", ""),
initialized_(false),
fast_mode_(0) {
fast_mode_(0),
external_morphs_("external_morphs", this, other.external_morphs_),
external_morph_indices_(other.external_morph_indices_)
{
initialize();
}

Expand Down Expand Up @@ -243,20 +247,42 @@ void CMSHistSum::initialize() const {
sentry_.addVars(coeffpars_);
binsentry_.addVars(binpars_);

for (const auto* morph : external_morphs_) {
RooArgSet* deps = morph->getParameters({*x_});
sentry_.addVars(*deps);
delete deps;
}

sentry_.setValueDirty();
binsentry_.setValueDirty();

initialized_ = true;
}

void CMSHistSum::updateMorphs() const {
// set up pointers ahead of time for quick loop
std::vector<CMSInterferenceFunc*> 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<CMSInterferenceFunc*>(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;
#endif
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<extdata.size(); ++ibin) {
compcache_[ip][ibin] *= extdata[ibin];
}
}
}
if (vtype_[ip] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
compcache_[ip].Log();
Expand Down Expand Up @@ -754,5 +780,24 @@ void CMSHistSum::EnableFastVertical() {
enable_fast_vertical_ = true;
}

void CMSHistSum::injectExternalMorph(int idx, CMSInterferenceFunc& morph) {
if ( idx >= 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

28 changes: 27 additions & 1 deletion test/test_interference.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand All @@ -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")

0 comments on commit 3e25193

Please sign in to comment.