Skip to content

Commit

Permalink
Adding getProcessNorms() to CMSHistErrorPropagator
Browse files Browse the repository at this point in the history
  • Loading branch information
Nick Wardle committed Jul 17, 2024
1 parent eea2ce3 commit 694e737
Show file tree
Hide file tree
Showing 3 changed files with 60 additions and 17 deletions.
2 changes: 2 additions & 0 deletions interface/CMSHistErrorPropagator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,8 @@ class CMSHistErrorPropagator : public RooAbsReal {
RooArgList wrapperList() const;
RooArgList const& coefList() const { return coeffs_; }
RooArgList const& funcList() const { return funcs_; }

std::map<std::string, Double_t> getProcessNorms() const;

friend class CMSHistV<CMSHistErrorPropagator>;

Expand Down
35 changes: 35 additions & 0 deletions src/CMSHistErrorPropagator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -549,5 +549,40 @@ RooArgList CMSHistErrorPropagator::wrapperList() const {
return result;
}

std::map<std::string, Double_t> CMSHistErrorPropagator::getProcessNorms() const {

std::map<std::string, Double_t> vals_;
RooArgList clist(coefList());
RooArgList plist(funcList());
/*if (plist.getSize() == 1) {
CMSHistErrorPropagator *err = dynamic_cast<CMSHistErrorPropagator*>(plist.at(0));
if (err) {
clist.removeAll();
plist.removeAll();
clist.add(err->coefList());
plist.add(err->wrapperList());
}
}
*/
for (int i = 0, n = clist.getSize(); i < n; ++i) {
RooAbsReal *coeff = (RooAbsReal *) clist.at(i);
std::string coeffName = coeff->GetName();
RooAbsReal* shape = (RooAbsReal*)plist.at(i);
std::unique_ptr<RooArgSet> myobs(shape->getObservables(*x_));
TString normProdName = TString::Format("%s", coeff->GetName());
RooAbsReal * normProd = nullptr;
if (coeff->ownedComponents()) {
normProd = dynamic_cast<RooAbsReal*>(coeff->ownedComponents()->find(normProdName));
}
if (!normProd) {
RooAbsReal* integral = shape->createIntegral(*myobs);
normProd = new RooProduct(normProdName, "", RooArgList(*integral, *coeff));
normProd->addOwnedComponents(*integral);
coeff->addOwnedComponents(*normProd);
}
vals_[normProdName.Data()] = normProd->getVal();
}
return vals_;
}
#undef HFVERBOSE

40 changes: 23 additions & 17 deletions test/printWorkspaceNormalisations.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,14 +49,14 @@
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(
# "",
# "--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(
"",
"--procFilter",
Expand Down Expand Up @@ -169,17 +169,23 @@ def find_chan_proc(name):
else:
chan_procs[chan] = [[proc, norm, 0, 0]]

use_cms_histsum = False

all_props = ws.allFunctions().selectByName("prop_bin*")
if all_props.getSize()>0: use_cms_histsum=True

# Look for cases where chan stored as CMSHistSum, set flag
if options.use_cms_histsum:
if use_cms_histsum:
chan_CMSHistSum_norms = {}
all_props = ws.allFunctions().selectByName("prop_bin*")
#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:
types = [ROOT.CMSHistSum, ROOT.CMSHistErrorPropagator]
if type(prop) in types :
chan_CMSHistSum_norms[chan] = dict(prop.getProcessNorms())


Expand Down Expand Up @@ -243,7 +249,7 @@ def printEndExpand():
continue
print("<hr>")
print("Top-level normalization for process {proc0} -> {proc1_name}<br>".format(proc0=proc[0], proc1_name=proc[1].GetName()))
if options.use_cms_histsum:
if use_cms_histsum:
if chan in chan_CMSHistSum_norms:
default_val = chan_CMSHistSum_norms[chan][proc[1].GetName()]
else:
Expand All @@ -253,7 +259,7 @@ def printEndExpand():
default_norms[chan][proc[1].GetName()] = default_val

if options.printValueOnly:
print("default value = {default_val}<br>".format(default_val=default_val))
print("default value = {default_val:.5f}<br>".format(default_val=default_val))
else:
if proc[2]:
proc_norm_var = ws.function("n_exp_bin%s_proc_%s" % (chan, proc[0]))
Expand Down Expand Up @@ -289,7 +295,7 @@ def printEndExpand():
print(" ... is a constant (RooRealVar)")
printEndExpand()

print(" default value = ", default_val, "<br>")
print(" default value = %.5f "%default_val, "<br>")
print("</tr>")

print(
Expand Down Expand Up @@ -320,7 +326,7 @@ def printEndExpand():
print("---------------------------------------------------------------------------")
print(" Top-level normalization for process %s -> %s" % (proc[0], proc[1].GetName()))
print("---------------------------------------------------------------------------")
if options.use_cms_histsum:
if use_cms_histsum:
if chan in chan_CMSHistSum_norms:
default_val = chan_CMSHistSum_norms[chan][proc[1].GetName()]
else:
Expand All @@ -330,7 +336,7 @@ def printEndExpand():
default_norms[chan][proc[1].GetName()] = default_val

if options.printValueOnly:
print(" default value = ", default_val)
print(" default value = %.5f "%default_val)
# if options.printValueOnly: print " --xcp %s:%s "%(chan,proc[0]),
else:
if proc[2]:
Expand All @@ -354,7 +360,7 @@ def printEndExpand():
proc[1].Print()
print(" ... is a constant (RooRealVar)")
print(" -------------------------------------------------------------------------")
print(" default value = ", default_val)
print(" default value = %.5f "%default_val)

# Save norms to json file
if options.output_json != "":
Expand Down

0 comments on commit 694e737

Please sign in to comment.