Skip to content

GitLab

Commit ba9ff106 authored by Marcel Rieger's avatar Marcel Rieger
Browse files

Fix scaling of eft cross sections and limits.

parent af3a78e2
Hide whitespace changes
Inline Side-by-side
dhi/eft_tools.py
View file @ ba9ff106
......@@ -14,7 +14,7 @@ from dhi.util import make_list
class EFTCrossSectionProvider(object):
"""
Helper class to calculate HH cross sections in EFT, as usualy in units of pb.
Helper class to calculate HH cross sections in EFT in units of pb.
Coefficients and formulae are taken from
https://github.com/fabio-mon/HHStatAnalysis/blob/c8fc33d2ae3f7e04cfc83e773e2880657ffdce3b/AnalyticalModels/python/NonResonantModelNLO.py
with credits to F. Monti and P. Mandrik.
......@@ -23,21 +23,22 @@ class EFTCrossSectionProvider(object):
def __init__(self):
super(EFTCrossSectionProvider, self).__init__()
# from https://github.com/pmandrik/VSEVA/blob/f7224649297f900a4ae25cf721d65cae8bd7b408/HHWWgg/reweight/reweight_HH.C#L117
self.coeffs_ggf_nlo_13tev = [
# ggf nlo coefficients in pb, converted from fb values
# https://github.com/pmandrik/VSEVA/blob/f7224649297f900a4ae25cf721d65cae8bd7b408/HHWWgg/reweight/reweight_HH.C#L117
self.coeffs_ggf_nlo_13tev = [0.001 * c for c in [
62.5088, 345.604, 9.63451, 4.34841, 39.0143, -268.644, -44.2924, 96.5595, 53.515,
-155.793, -23.678, 54.5601, 12.2273, -26.8654, -19.3723, -0.0904439, 0.321092, 0.452381,
-0.0190758, -0.607163, 1.27408, 0.364487, -0.499263,
]
]]
self.ggf_xsec_sm_nnlo = 0.03105
self.ggf_xsec_sm_nnlo = 0.03105 # pb
def get_ggf_xsec_nlo(self, kl=1., kt=1., c2=0., cg=0., c2g=0., coeffs=None):
if coeffs is None:
coeffs = self.coeffs_ggf_nlo_13tev
return 0.001 * (
coeffs[0] * kt**4 + \
return (
coeffs[0] * kt**4 +
coeffs[1] * c2**2 +
coeffs[2] * kt**2 * kl**2 +
coeffs[3] * cg**2 * kl**2 +
......@@ -63,13 +64,10 @@ class EFTCrossSectionProvider(object):
)
def get_ggf_xsec_nnlo(self, kl=1., kt=1., c2=0., cg=0., c2g=0., coeffs=None):
if coeffs is None:
coeffs = self.coeffs_ggf_nlo_13tev
xsec_bsm_nlo = self.get_ggf_xsec_nlo(kl=kl, kt=kt, c2=c2, cg=cg, c2g=c2g, coeffs=coeffs)
xsec_sm_nlo = self.get_ggf_xsec_nlo(kl=1., kt=1., c2=0., cg=0., c2g=0., coeffs=coeffs)
return self.ggf_xsec_sm_nnlo * xsec_bsm_nlo / xsec_sm_nlo
return xsec_bsm_nlo * self.ggf_xsec_sm_nnlo / xsec_sm_nlo
#: EFTCrossSectionProvider singleton.
......
dhi/plots/limits.py
View file @ ba9ff106
......@@ -709,6 +709,7 @@ def plot_benchmark_limits(
y_max=None,
y_log=False,
xsec_unit="fb",
hh_process=None,
campaign=None,
bar_width=0.66,
):
......@@ -739,9 +740,12 @@ def plot_benchmark_limits(
*y_min* and *y_max* define the y axis range and default to the range of the given values. When
*y_log* is *True*, the y-axis is plotted with a logarithmic scale. *xsec_unit* denotes the unit
of the passed values and defaults to fb. *campaign* should refer to the name of a campaign label
defined in dhi.config.campaign_labels. The *bar_width* should be a value between 0 and 1 and
controls the fraction of the limit bar width relative to the bin width.
of the passed values and defaults to fb. *hh_process* can be the name of a HH subprocess
configured in *dhi.config.br_hh_names* and is inserted to the process name in the title of the
y-axis and indicates that the plotted cross section data was (e.g.) scaled by a branching ratio.
*campaign* should refer to the name of a campaign label defined in dhi.config.campaign_labels.
The *bar_width* should be a value between 0 and 1 and controls the fraction of the limit bar
width relative to the bin width.
Example: https://cms-hh.web.cern.ch/tools/inference/tasks/eft.html#benchmark-limits TODO
"""
......@@ -787,7 +791,7 @@ def plot_benchmark_limits(
# dummy histogram to control axes
x_title = "Shape benchmark"
y_title = "Upper 95% CLs limit on #sigma({}) / {}".format(
create_hh_process_label(poi), to_root_latex(xsec_unit))
create_hh_process_label(poi, hh_process), to_root_latex(xsec_unit))
h_dummy = ROOT.TH1F("dummy", ";{};{}".format(x_title, y_title), n, -0.5, n - 0.5)
r.setup_hist(h_dummy, pad=pad, props={"LineWidth": 0, "Minimum": y_min, "Maximum": y_max})
r.setup_x_axis(h_dummy.GetXaxis(), pad=pad, props={"Ndivisions": n})
......
dhi/tasks/eft.py
View file @ ba9ff106
......@@ -22,6 +22,7 @@ from dhi.eft_tools import (
get_eft_ggf_xsec_nnlo, sort_eft_benchmark_names, sort_eft_scan_names,
extract_eft_scan_parameter,
)
from dhi.config import br_hh
class EFTBase(MultiDatacardTask):
......@@ -304,6 +305,17 @@ class MergeEFTUpperLimits(EFTScanBase):
class PlotEFTBenchmarkLimits(EFTBenchmarkBase, PlotTask):
xsec = luigi.ChoiceParameter(
default="fb",
choices=["pb", "fb"],
description="convert limits to cross sections in this unit; choices: pb,fb; default: fb",
)
br = luigi.ChoiceParameter(
default=law.NO_STR,
choices=[law.NO_STR, ""] + list(br_hh.keys()),
description="name of a branching ratio defined in dhi.config.br_hh to scale the cross "
"section when xsec is used; choices: '',{}; no default".format(",".join(br_hh.keys())),
)
y_log = luigi.BoolParameter(
default=False,
description="apply log scaling to the y-axis; default: False",
......@@ -318,11 +330,15 @@ class PlotEFTBenchmarkLimits(EFTBenchmarkBase, PlotTask):
return MergeEFTBenchmarkLimits.req(self)
def output(self):
parts = self.get_output_postfix(join=False)
# additional postfix
parts = []
parts.append(self.xsec)
if self.br and self.br != law.NO_STR:
parts.append(self.br)
if self.y_log:
parts.append("log")
name = self.create_plot_name(["eftlimits"] + parts)
name = self.create_plot_name(["eftlimits", self.get_output_postfix(), parts])
return self.local_target(name)
@law.decorator.log
......@@ -337,15 +353,18 @@ class PlotEFTBenchmarkLimits(EFTBenchmarkBase, PlotTask):
# load limit values
limit_values = self.input().load(formatter="numpy")["data"]
# prepare conversion scale
scale = br_hh.get(self.br, 1.) * {"fb": 1., "pb": 0.001}[self.xsec]
# fill data entries as expected by the plot function
data = []
for name, record in zip(self.benchmark_datacards.keys(), limit_values):
entry = {
"name": name,
"expected": record.tolist()[:5],
"expected": [v * scale for v in record.tolist()[:5]],
}
if self.unblinded:
entry["observed"] = float(record[5])
entry["observed"] = float(record[5]) * scale
data.append(entry)
# some printing
......@@ -360,13 +379,25 @@ class PlotEFTBenchmarkLimits(EFTBenchmarkBase, PlotTask):
y_min=self.get_axis_limit("y_min"),
y_max=self.get_axis_limit("y_max"),
y_log=self.y_log,
xsec_unit="fb",
xsec_unit=self.xsec,
hh_process=self.br if self.br in br_hh else None,
campaign=self.campaign if self.campaign != law.NO_STR else None,
)
class PlotEFTUpperLimits(EFTScanBase, PlotTask):
xsec = luigi.ChoiceParameter(
default=law.NO_STR,
choices=[law.NO_STR, "", "pb", "fb"],
description="convert limits to cross sections in this unit; choices: '',pb,fb; no default",
)
br = luigi.ChoiceParameter(
default=law.NO_STR,
choices=[law.NO_STR, ""] + list(br_hh.keys()),
description="name of a branching ratio defined in dhi.config.br_hh to scale the cross "
"section when xsec is used; choices: '',{}; no default".format(",".join(br_hh.keys())),
)
y_log = luigi.BoolParameter(
default=False,
description="apply log scaling to the y-axis; default: False",
......@@ -379,11 +410,16 @@ class PlotEFTUpperLimits(EFTScanBase, PlotTask):
return MergeEFTUpperLimits.req(self)
def output(self):
parts = self.get_output_postfix(join=False)
# additional postfix
parts = []
if self.xsec in ["pb", "fb"]:
parts.append(self.xsec)
if self.br not in (law.NO_STR, ""):
parts.append(self.br)
if self.y_log:
parts.append("log")
name = self.create_plot_name(["eftlimits"] + parts)
name = self.create_plot_name(["eftlimits", self.get_output_postfix(), parts])
return self.local_target(name)
@law.decorator.log
......@@ -391,20 +427,50 @@ class PlotEFTUpperLimits(EFTScanBase, PlotTask):
@view_output_plots
@law.decorator.safe_output
def run(self):
import numpy as np
# prepare the output
output = self.output()
output.parent.touch()
# load limit values
# load limit values, given in fb
limit_values = self.input().load(formatter="numpy")["data"]
# get theory values
# get nnlo cross sections in pb
xsecs = OrderedDict(
(v, get_eft_ggf_xsec_nnlo(**{self.scan_parameter: v}))
for v in limit_values[self.scan_parameter]
)
# prepare the br value
br_value = br_hh.get(self.br, 1.)
# get scaling factors for xsec and limit values
# apply scaling by xsec and br
n = len(xsecs)
if self.xsec == "fb":
xsec_unit = "fb"
limit_scales = n * [br_value]
xsec_scales = n * [1000. * br_value]
elif self.xsec == "pb":
xsec_unit = "pb"
limit_scales = n * [0.001 * br_value]
xsec_scales = n * [br_value]
else:
xsec_unit = None
limit_scales = (np.array(xsecs.values()) * 1000.)**-1
xsec_scales = [xsec**-1 for xsec in xsecs.values()]
# apply scales
for name in limit_values.dtype.names:
if name in ["limit", "limit_p1", "limit_m1", "limit_p2", "limit_m2", "observed"]:
limit_values[name] *= limit_scales
xsecs = OrderedDict((v, xsec * s) for (v, xsec), s in zip(xsecs.items(), xsec_scales))
# preate theory values in the correct structure
thy_values = {
self.scan_parameter: limit_values[self.scan_parameter],
"xsec": [
get_eft_ggf_xsec_nnlo(**{self.scan_parameter: v}) * 0.001
for v in limit_values[self.scan_parameter]
],
"xsec": list(xsecs.values()),
}
# prepare observed values
......@@ -432,7 +498,8 @@ class PlotEFTUpperLimits(EFTScanBase, PlotTask):
y_min=self.get_axis_limit("y_min"),
y_max=self.get_axis_limit("y_max"),
y_log=self.y_log,
xsec_unit="fb",
xsec_unit=xsec_unit,
hh_process=self.br if xsec_unit and self.br in br_hh else None,
model_parameters=model_parameters,
campaign=self.campaign if self.campaign != law.NO_STR else None,
)
dhi/tasks/limits.py
View file @ ba9ff106
......@@ -182,7 +182,7 @@ class PlotUpperLimits(UpperLimitsBase, POIPlotTask):
parts = []
if self.xsec in ["pb", "fb"]:
parts.append(self.xsec)
if self.br not in (law.NO_STR, ""):
if self.br and self.br != law.NO_STR:
parts.append(self.br)
if self.y_log:
parts.append("log")
......@@ -268,7 +268,7 @@ class PlotUpperLimits(UpperLimitsBase, POIPlotTask):
y_max=self.get_axis_limit("y_max"),
y_log=self.y_log,
xsec_unit=xsec_unit,
hh_process=None if self.br in (None, law.NO_STR) else self.br,
hh_process=self.br if xsec_unit and self.br in br_hh else None,
model_parameters=self.get_shown_parameters(),
campaign=self.campaign if self.campaign != law.NO_STR else None,
)
......@@ -318,7 +318,7 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
parts = []
if self.xsec in ["pb", "fb"]:
parts.append(self.xsec)
if self.br not in (law.NO_STR, ""):
if self.br and self.br != law.NO_STR:
parts.append(self.br)
if self.y_log:
parts.append("log")
......@@ -405,7 +405,7 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
y_max=self.get_axis_limit("y_max"),
y_log=self.y_log,
xsec_unit=xsec_unit,
hh_process=None if self.br in (None, law.NO_STR) else self.br,
hh_process=self.br if xsec_unit and self.br in br_hh else None,
model_parameters=self.get_shown_parameters(),
campaign=self.campaign if self.campaign != law.NO_STR else None,
)
......@@ -429,7 +429,7 @@ class PlotMultipleUpperLimitsByModel(PlotUpperLimits, MultiHHModelTask):
parts = []
if self.xsec in ["pb", "fb"]:
parts.append(self.xsec)
if self.br not in (law.NO_STR, ""):
if self.br and self.br != law.NO_STR:
parts.append(self.br)
if self.y_log:
parts.append("log")
......@@ -523,7 +523,7 @@ class PlotMultipleUpperLimitsByModel(PlotUpperLimits, MultiHHModelTask):
y_max=self.get_axis_limit("y_max"),
y_log=self.y_log,
xsec_unit=xsec_unit,
hh_process=None if self.br in (None, law.NO_STR) else self.br,
hh_process=self.br if xsec_unit and self.br in br_hh else None,
model_parameters=self.get_shown_parameters(),
campaign=self.campaign if self.campaign != law.NO_STR else None,
)
......@@ -591,7 +591,7 @@ class PlotUpperLimitsAtPoint(POIPlotTask, MultiDatacardTask):
parts = []
if self.xsec in ["pb", "fb"]:
parts.append(self.xsec)
if self.br not in (law.NO_STR, ""):
if self.br and self.br != law.NO_STR:
parts.append(self.br)
if self.x_log:
parts.append("log")
......@@ -680,7 +680,7 @@ class PlotUpperLimitsAtPoint(POIPlotTask, MultiDatacardTask):
x_max=self.get_axis_limit("x_max"),
x_log=self.x_log,
xsec_unit=xsec_unit,
hh_process=None if self.br in (None, law.NO_STR) else self.br,
hh_process=self.br if xsec_unit and self.br in br_hh else None,
left_margin=None if self.left_margin == law.NO_INT else self.left_margin,
model_parameters=self.get_shown_parameters(),
h_lines=self.h_lines,
......
dhi/tasks/test.py
View file @ ba9ff106
......@@ -263,11 +263,13 @@ class TestPlots(six.with_metaclass(TestRegister, AnalysisTask)):
if self.check_enabled("eft_benchmark_limits"):
reqs["eft_benchmark_limits"] = PlotEFTBenchmarkLimits.req(self,
multi_datacards=(eft_cards_bm,),
xsec="fb",
)
if self.check_enabled("eft_upper_limits"):
reqs["eft_upper_limits"] = PlotEFTUpperLimits.req(self,
multi_datacards=(eft_cards_c2,),
xsec="fb",
)
if self.check_enabled("morphing_scales"):
......
law @ c43ef2ed
Compare 8e3ecd42...c43ef2ed
Subproject commit 8e3ecd4281418d846179045de6fa1e5badb6ecbb
Subproject commit c43ef2ed3e03e57580f00ff1f6edecada8d91be8
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment