diff --git a/.gitignore b/.gitignore index 92a490e5..2bf8914d 100644 --- a/.gitignore +++ b/.gitignore @@ -34,9 +34,9 @@ src/* !src/HHbbVV/postprocessing/** src/HHbbVV/postprocessing/data/* src/HHbbVV/postprocessing/*.pkl +src/HHbbVV/postprocessing/*.txt src/HHbbVV/postprocessing/cards src/HHbbVV/postprocessing/templates_old -src/HHbbVV/postprocessing/*.txt src/HHbbVV/postprocessing/condor_templates src/HHbbVV/postprocessing/outs diff --git a/README.md b/README.md index 3be1f8f5..eb6eda59 100644 --- a/README.md +++ b/README.md @@ -289,7 +289,7 @@ nohup bash_scripts/res_bias_templates.sh $TAG &> bias.txt & ### Create Datacard -Need `root==6.22.6`, and `square_coef` branch of https://github.com/rkansal47/rhalphalib installed (`pip install -e . --user` after checking out the branch). `CMSSW_11_2_0` recommended. +Need `root==6.22.6`, and `square_coef` branch of https://github.com/rkansal47/rhalphalib installed (`pip install -e . --user` after checking out the branch). ```bash python3 postprocessing/CreateDatacard.py --templates-dir templates/$TAG --model-name $TAG (--resonant) diff --git a/src/HHbbVV/postprocessing/CreateDatacard.py b/src/HHbbVV/postprocessing/CreateDatacard.py index d9c2fcf9..9e80a574 100644 --- a/src/HHbbVV/postprocessing/CreateDatacard.py +++ b/src/HHbbVV/postprocessing/CreateDatacard.py @@ -5,13 +5,12 @@ Based on https://github.com/LPC-HH/combine-hh/blob/master/create_datacard.py -Author: Raghav Kansal +Authors: Raghav Kansal, Andres Nava """ import os, sys from typing import Dict, List, Tuple, Union -from dataclasses import dataclass, field import numpy as np import pickle, json @@ -20,10 +19,15 @@ import hist from hist import Hist + import rhalphalib as rl -rl.util.install_roofit_helpers() -rl.ParametericSample.PreferRooParametricHist = False +try: + rl.util.install_roofit_helpers() + rl.ParametericSample.PreferRooParametricHist = False +except: + print("rootfit install failed (not an issue for VBF)") + # logging.basicConfig(level=logging.DEBUG) logging.basicConfig(level=logging.INFO) adjust_posdef_yields = False @@ -31,73 +35,24 @@ # from utils import add_bool_arg from hh_vars import LUMI, res_sig_keys, qcd_key, data_key, years, jecs, jmsr, res_mps import utils +import datacardHelpers as helpers +from datacardHelpers import ( + Syst, + ShapeVar, + add_bool_arg, + rem_neg, + sum_templates, + combine_templates, + get_effect_updown, + join_with_padding, + rebin_channels, + get_channels, + mxmy, +) import argparse -@dataclass -class Syst: - """For storing info about systematics""" - - name: str = None - prior: str = None # e.g. "lnN", "shape", etc. - - # float if same value in all regions/samples, dictionary of values per region/sample if not - # if both, region should be the higher level of the dictionary - value: Union[float, Dict[str, float]] = None - value_down: Union[float, Dict[str, float]] = None # if None assumes symmetric effect - # if the value is different for different regions or samples - diff_regions: bool = False - diff_samples: bool = False - - samples: List[str] = None # samples affected by it - # in case of uncorrelated unc., which years to split into - uncorr_years: List[str] = field(default_factory=lambda: years) - pass_only: bool = False # is it applied only in the pass regions - - def __post_init__(self): - if isinstance(self.value, dict) and not (self.diff_regions or self.diff_samples): - raise RuntimeError( - f"Value for systematic is a dictionary but neither ``diff_regions`` nor ``diff_samples`` is set." - ) - - -@dataclass -class ShapeVar: - """For storing and calculating info about variables used in fit""" - - name: str = None - bins: np.ndarray = None # bin edges - order: int = None # TF order - - def __post_init__(self): - # use bin centers for polynomial fit - self.pts = self.bins[:-1] + 0.5 * np.diff(self.bins) - # scale to be between [0, 1] - self.scaled = (self.pts - self.bins[0]) / (self.bins[-1] - self.bins[0]) - - -def add_bool_arg(parser, name, help, default=False, no_name=None): - """Add a boolean command line argument for argparse""" - varname = "_".join(name.split("-")) # change hyphens to underscores - group = parser.add_mutually_exclusive_group(required=False) - group.add_argument("--" + name, dest=varname, action="store_true", help=help) - if no_name is None: - no_name = "no-" + name - no_help = "don't " + help - else: - no_help = help - group.add_argument("--" + no_name, dest=varname, action="store_false", help=no_help) - parser.set_defaults(**{varname: default}) - - -def mxmy(sample): - mY = int(sample.split("-")[-1]) - mX = int(sample.split("NMSSM_XToYHTo2W2BTo4Q2B_MX-")[1].split("_")[0]) - - return (mX, mY) - - parser = argparse.ArgumentParser() parser.add_argument( "--templates-dir", @@ -105,6 +60,9 @@ def mxmy(sample): type=str, help="input pickle file of dict of hist.Hist templates", ) + +add_bool_arg(parser, "vbf", "VBF category datacards", default=False) + add_bool_arg(parser, "sig-separate", "separate templates for signals and bgs", default=False) add_bool_arg(parser, "do-jshifts", "Do JEC/JMC corrections.", default=True) @@ -119,7 +77,12 @@ def mxmy(sample): parser.add_argument("--cards-dir", default="cards", type=str, help="output card directory") parser.add_argument("--mcstats-threshold", default=100, type=float, help="mcstats threshold n_eff") -parser.add_argument("--epsilon", default=1e-3, type=float, help="epsilon to avoid numerical errs") +parser.add_argument( + "--epsilon", + default=1e-2, + type=float, + help="epsilon to avoid numerical errs - also used to decide whether to add mc stats error", +) parser.add_argument( "--scale-templates", default=None, type=float, help="scale all templates for bias tests" ) @@ -155,6 +118,7 @@ def mxmy(sample): CMS_PARAMS_LABEL = "CMS_bbWW_hadronic" if not args.resonant else "CMS_XHYbbWW_boosted" +qcd_data_key = "qcd_datadriven" if args.nTF is None: args.nTF = [1, 2] if args.resonant else [0] @@ -177,12 +141,11 @@ def mxmy(sample): "ggHH_kl_0_kt_1_HHbbVV", ] nonres_sig_keys_vbf = [ - "qqHH_CV_1_C2V_1_kl_1_HHbbVV", + "VBFHHbbVV", "qqHH_CV_1_C2V_0_kl_1_HHbbVV", "qqHH_CV_1p5_C2V_1_kl_1_HHbbVV", "qqHH_CV_1_C2V_1_kl_2_HHbbVV", "qqHH_CV_1_C2V_2_kl_1_HHbbVV", - "qqHH_CV_1_C2V_2_kl_1_HHbbVV", "qqHH_CV_1_C2V_1_kl_0_HHbbVV", "qqHH_CV_0p5_C2V_1_kl_1_HHbbVV", ] @@ -204,11 +167,18 @@ def mxmy(sample): print("--sig-sample needs to be specified for resonant cards") sys.exit() else: + # change names to match HH combination convention for key in nonres_sig_keys: - mc_samples[key] = key.replace("HHbbVV", "hbbhww") + # check in case single sig sample is specified + if args.sig_sample is None or key == args.sig_sample: + if key == "HHbbVV": + mc_samples["HHbbVV"] = "ggHH_kl_1_kt_1_hbbhww" + elif key == "VBFHHbbVV": + mc_samples["VBFHHbbVV"] = "qqHH_CV_1_C2V_1_kl_1_HHbbww" + else: + mc_samples[key] = key.replace("HHbbVV", "hbbhww") - mc_samples["HHbbVV"] = "ggHH_kl_1_kt_1_hbbhww" - sig_keys = nonres_sig_keys + sig_keys.append(key) all_mc = list(mc_samples.keys()) @@ -352,186 +322,6 @@ def mxmy(sample): ) -def main(args): - # (pass, fail) x (unblinded, blinded) - regions: List[str] = [ - f"{pf}{blind_str}" for pf in [f"pass", "fail"] for blind_str in ["", "Blinded"] - ] - - # templates per region per year, templates per region summed across years - templates_dict, templates_summed = get_templates( - args.templates_dir, years, args.sig_separate, args.scale_templates, args.combine_lasttwo - ) - - # TODO: check if / how to include signal trig eff uncs. (rn only using bg uncs.) - process_systematics(args.templates_dir, args.sig_separate) - - # random template from which to extract shape vars - sample_templates: Hist = templates_summed[regions[0]] - # [mH(bb)] for nonresonant, [mY, mX] for resonant - shape_vars = [ - ShapeVar(name=axis.name, bins=axis.edges, order=args.nTF[i]) - for i, axis in enumerate(sample_templates.axes[1:]) - ] - - # build actual fit model now - model = rl.Model("HHModel" if not args.resonant else "XHYModel") - - # Fill templates per sample, incl. systematics - # TODO: blinding for resonant - fill_args = [ - model, - regions, - templates_dict, - templates_summed, - mc_samples, - nuisance_params, - nuisance_params_dict, - corr_year_shape_systs, - uncorr_year_shape_systs, - shape_systs_dict, - args.bblite, - ] - - fit_args = [model, shape_vars, templates_summed, args.scale_templates, args.min_qcd_val] - - if args.resonant: - # fill 1 channel per mX bin - for i in range(len(shape_vars[1].pts)): - logging.info(f"\n\nFilling templates for mXbin {i}") - fill_regions(*fill_args, mX_bin=i) - - res_alphabet_fit(*fit_args) - else: - fill_regions(*fill_args) - nonres_alphabet_fit(*fit_args) - - ############################################## - # Save model - ############################################## - - logging.info("rendering combine model") - - os.system(f"mkdir -p {args.cards_dir}") - - out_dir = ( - os.path.join(str(args.cards_dir), args.model_name) - if args.model_name is not None - else args.cards_dir - ) - model.renderCombine(out_dir) - - with open(f"{out_dir}/model.pkl", "wb") as fout: - pickle.dump(model, fout, 2) # use python 2 compatible protocol - - -def _rem_neg(template_dict: Dict): - for sample, template in template_dict.items(): - template.values()[template.values() < 0] = 0 - - return template_dict - - -def _match_samples(template: Hist, match_template: Hist) -> Hist: - """Temporary solution for case where 2018 templates don't have L1 prefiring in their axis - - Args: - template (Hist): template to which samples may need to be added - match_template (Hist): template from which to extract extra samples - """ - samples = list(match_template.axes[0]) - - # if template already has all the samples, don't do anything - if list(template.axes[0]) == samples: - return template - - # otherwise remake template with samples from ``match_template`` - h = hist.Hist(*match_template.axes, storage="weight") - - for sample in template.axes[0]: - sample_index = np.where(np.array(list(h.axes[0])) == sample)[0][0] - h.view()[sample_index] = template[sample, ...].view() - - return h - - -def sum_templates(template_dict: Dict): - """Sum templates across years""" - - ttemplate = list(template_dict.values())[0] # sample templates to extract values from - combined = {} - - for region in ttemplate: - thists = [] - - for year in years: - # temporary solution for case where 2018 templates don't have L1 prefiring in their axis - thists.append( - _match_samples(template_dict[year][region], template_dict["2016"][region]) - ) - - combined[region] = sum(thists) - - return combined - - -def combine_templates( - bg_templates: Dict[str, Hist], sig_templates: List[Dict[str, Hist]] -) -> Dict[str, Hist]: - """ - Combines BG and signal templates into a single Hist (per region). - - Args: - bg_templates (Dict[str, Hist]): dictionary of region -> Hist - sig_templates (List[Dict[str, Hist]]): list of dictionaries of region -> Hist for each - signal samples - """ - ctemplates = {} - - for region, bg_template in bg_templates.items(): - # combined sig + bg samples - csamples = list(bg_template.axes[0]) + [ - s for sig_template in sig_templates for s in list(sig_template[region].axes[0]) - ] - - # new hist with all samples - ctemplate = Hist( - hist.axis.StrCategory(csamples, name="Sample"), - *bg_template.axes[1:], - storage="weight", - ) - - # add background hists - for sample in bg_template.axes[0]: - sample_key_index = np.where(np.array(list(ctemplate.axes[0])) == sample)[0][0] - ctemplate.view(flow=True)[sample_key_index, ...] = bg_template[sample, ...].view( - flow=True - ) - - # add signal hists - for st in sig_templates: - sig_template = st[region] - for sample in sig_template.axes[0]: - sample_key_index = np.where(np.array(list(ctemplate.axes[0])) == sample)[0][0] - ctemplate.view(flow=True)[sample_key_index, ...] = sig_template[sample, ...].view( - flow=True - ) - - ctemplates[region] = ctemplate - - return ctemplates - - -def _combine_last_two_bins(templates_dict): - for year in years: - for region in templates_dict[year]: - templates_dict[year][region] = utils.rebin_hist( - templates_dict[year][region], - "bbFatJetParticleNetMass", - list(range(50, 240, 10)) + [250], - ) - - def get_templates( templates_dir: str, years: List[str], @@ -546,18 +336,18 @@ def get_templates( # signal and background templates in same hist, just need to load and sum across years for year in years: with open(f"{templates_dir}/{year}_templates.pkl", "rb") as f: - templates_dict[year] = _rem_neg(pickle.load(f)) + templates_dict[year] = rem_neg(pickle.load(f)) else: # signal and background in different hists - need to combine them into one hist for year in years: with open(f"{templates_dir}/backgrounds/{year}_templates.pkl", "rb") as f: - bg_templates = _rem_neg(pickle.load(f)) + bg_templates = rem_neg(pickle.load(f)) sig_templates = [] for sig_key in sig_keys: with open(f"{templates_dir}/{hist_names[sig_key]}/{year}_templates.pkl", "rb") as f: - sig_templates.append(_rem_neg(pickle.load(f))) + sig_templates.append(rem_neg(pickle.load(f))) templates_dict[year] = combine_templates(bg_templates, sig_templates) @@ -579,9 +369,9 @@ def get_templates( templates_dict[year][key].variances()[j, ...] = variances * (scale**2) if combine_lasttwo: - _combine_last_two_bins(templates_dict) + helpers.combine_last_two_bins(templates_dict, years) - templates_summed: Dict[str, Hist] = sum_templates(templates_dict) # sum across years + templates_summed: Dict[str, Hist] = sum_templates(templates_dict, years) # sum across years return templates_dict, templates_summed @@ -615,8 +405,6 @@ def process_systematics_combined(systematics: Dict): nuisance_params[f"{CMS_PARAMS_LABEL}_triggerEffSF_uncorrelated"].value = tdict - print("Nuisance Parameters\n", nuisance_params) - def process_systematics_separate(bg_systematics: Dict, sig_systs: Dict[str, Dict]): """Get total uncertainties from per-year systs separated into bg and sig systs""" @@ -672,6 +460,48 @@ def process_systematics(templates_dir: str, sig_separate: bool): process_systematics_separate(bg_systematics, sig_systs) # LP SF and trig effs. +# TODO: separate function for VBF? +def get_year_updown( + templates_dict, sample, region, region_noblinded, blind_str, year, skey, mX_bin=None, vbf=False +): + """ + Return templates with only the given year's shapes shifted up and down by the ``skey`` systematic. + Returns as [up templates, down templates] + """ + updown = [] + + for shift in ["up", "down"]: + sshift = f"{skey}_{shift}" + # get nominal templates for each year + templates = {y: templates_dict[y][region][sample, ...] for y in years} + + # replace template for this year with the shifted template + if skey in jecs or skey in jmsr: + # JEC/JMCs saved as different "region" in dict + reg_name = ( + f"{region_noblinded}_{sshift}{blind_str}" + if mX_bin is None + else f"{region}_{sshift}" + ) + + templates[year] = templates_dict[year][reg_name][sample, ...] + else: + # weight uncertainties saved as different "sample" in dict + templates[year] = templates_dict[year][region][f"{sample}_{sshift}", ...] + + if mX_bin is not None: + for year, template in templates.items(): + templates[year] = template[:, mX_bin] + + # sum templates with year's template replaced with shifted + if vbf: + updown.append(sum([t.value for t in templates.values()])) + else: + updown.append(sum(list(templates.values())).values()) + + return updown + + def fill_regions( model: rl.Model, regions: List[str], @@ -932,7 +762,7 @@ def nonres_alphabet_fit( # was integer, and numpy complained about subtracting float from it initial_qcd = failCh.getObservation().astype(float) for sample in failCh: - if args.resonant and sample.sampletype == rl.Sample.SIGNAL: + if sample.sampletype == rl.Sample.SIGNAL: continue logging.debug("subtracting %s from qcd" % sample._name) initial_qcd -= sample.getExpectation(nominal=True) @@ -1069,108 +899,317 @@ def res_alphabet_fit( passCh.addSample(pass_qcd) -def _shape_checks(values_up, values_down, values_nominal, effect_up, effect_down, logger): - norm_up = np.sum(values_up) - norm_down = np.sum(values_down) - norm_nominal = np.sum(values_nominal) - prob_up = values_up / norm_up - prob_down = values_down / norm_down - prob_nominal = values_nominal / norm_nominal - shapeEffect_up = np.sum( - np.abs(prob_up - prob_nominal) / (np.abs(prob_up) + np.abs(prob_nominal)) - ) - shapeEffect_down = np.sum( - np.abs(prob_down - prob_nominal) / (np.abs(prob_down) + np.abs(prob_nominal)) - ) +def createDatacardAlphabet(args, templates_dict, templates_summed, shape_vars): + # (pass, fail) x (unblinded, blinded) + regions: List[str] = [ + f"{pf}{blind_str}" for pf in [f"pass", "fail"] for blind_str in ["", "Blinded"] + ] - valid = True - if np.allclose(effect_up, 1.0) and np.allclose(effect_down, 1.0): - valid = False - logger.warning("No shape effect") - elif np.allclose(effect_up, effect_down): - valid = False - logger.warning("Up is the same as Down, but different from nominal") - elif np.allclose(effect_up, 1.0) or np.allclose(effect_down, 1.0): - valid = False - logger.warning("Up or Down is the same as nominal (one-sided)") - elif shapeEffect_up < 0.001 and shapeEffect_down < 0.001: - valid = False - logger.warning("No genuine shape effect (just norm)") - elif (norm_up > norm_nominal and norm_down > norm_nominal) or ( - norm_up < norm_nominal and norm_down < norm_nominal - ): - valid = False - logger.warning("Up and Down vary norm in the same direction") - # print("Norm nominal", norm_nominal) - # print("Norm up", norm_up) - # print("Norm down", norm_down) - # print("values nominal", values_nominal) - # print("values up", values_up) - # print("values down", values_down) + # build actual fit model now + model = rl.Model("HHModel" if not args.resonant else "XHYModel") + + # Fill templates per sample, incl. systematics + fill_args = [ + model, + regions, + templates_dict, + templates_summed, + mc_samples, + nuisance_params, + nuisance_params_dict, + corr_year_shape_systs, + uncorr_year_shape_systs, + shape_systs_dict, + args.bblite, + ] + + fit_args = [model, shape_vars, templates_summed, args.scale_templates, args.min_qcd_val] + + if args.resonant: + # fill 1 channel per mX bin + for i in range(len(shape_vars[1].pts)): + logging.info(f"\n\nFilling templates for mXbin {i}") + fill_regions(*fill_args, mX_bin=i) - # if valid: - # logger.info("Shapes are valid") + res_alphabet_fit(*fit_args) + else: + fill_regions(*fill_args) + nonres_alphabet_fit(*fit_args) + ############################################## + # Save model + ############################################## -def get_effect_updown(values_nominal, values_up, values_down, mask, logger): - effect_up = np.ones_like(values_nominal) - effect_down = np.ones_like(values_nominal) + logging.info("rendering combine model") - mask_up = mask & (values_up >= 0) - mask_down = mask & (values_down >= 0) + out_dir = ( + os.path.join(str(args.cards_dir), args.model_name) + if args.model_name is not None + else args.cards_dir + ) + model.renderCombine(out_dir) - effect_up[mask_up] = values_up[mask_up] / values_nominal[mask_up] - effect_down[mask_down] = values_down[mask_down] / values_nominal[mask_down] + with open(f"{out_dir}/model.pkl", "wb") as fout: + pickle.dump(model, fout, 2) # use python 2 compatible protocol - zero_up = values_up == 0 - zero_down = values_down == 0 - effect_up[mask_up & zero_up] = values_nominal[mask_up & zero_up] * args.epsilon - effect_down[mask_down & zero_down] = values_nominal[mask_down & zero_down] * args.epsilon +def fill_yields(channels, channels_summed): + """Fill top-level datacard info for VBF""" - _shape_checks(values_up, values_down, values_nominal, effect_up, effect_down, logger) + # dict storing all the substitutions for the datacard template + datacard_dict = { + "num_bins": len(channels), + "num_bgs": len(bg_keys) + 1, # + 1 for qcd + "bins": join_with_padding(channels, padding=16), + "observations": join_with_padding( + [channels_summed[channel][data_key].value for channel in channels], padding=16 + ), + "qcdlabel": qcd_data_key, + } + + for i, l in enumerate(["A", "B", "C", "D"]): + # channel labels + datacard_dict[f"bin{l}"] = f"{channels[i]:<14}" + + # fill in data - MC yields for ABCD method + if i > 0: + data_obs = channels_summed[channels[i]][data_key].value + mc_bg_yields = sum([channels_summed[channels[i]][key].value for key in bg_keys]) + datacard_dict[f"dataqcd{l}"] = data_obs - mc_bg_yields + + # collecting MC samples and yields per chanel + channel_bins_dict = { + "bins_x_processes": [], + "processes_per_bin": [], + "processes_index": [], + "processes_rates": [], + } + + rates_dict = {} + + for channel in channels: + rates_dict[channel] = {} + for i, key in enumerate(sig_keys + bg_keys + [qcd_data_key]): + channel_bins_dict["bins_x_processes"].append(channel) + channel_bins_dict["processes_index"].append(i + 1 - len(sig_keys)) + if key == qcd_data_key: + channel_bins_dict["processes_per_bin"].append(qcd_data_key) + channel_bins_dict["processes_rates"].append(1) + else: + channel_bins_dict["processes_per_bin"].append(mc_samples[key]) + channel_bins_dict["processes_rates"].append(channels_summed[channel][key].value) + rates_dict[channel][key] = channels_summed[channel][key] + + for key, arr in channel_bins_dict.items(): + datacard_dict[key] = join_with_padding(arr) + + return datacard_dict, rates_dict + + +def get_systematics_abcd(channels, channels_dict, channels_summed, rates_dict): + channel_systs_dict = {} + + for region in channels: + logging.info("starting region: %s" % region) - logging.debug("nominal : {nominal}".format(nominal=values_nominal)) - logging.debug("effect_up : {effect_up}".format(effect_up=effect_up)) - logging.debug("effect_down: {effect_down}".format(effect_down=effect_down)) + channel_systs_dict[region] = {} - return effect_up, effect_down + pass_region = region.startswith("pass") + region_nosidebands = region.split("_sidebands")[0] + sideband_str = "_sidebands" if region.endswith("_sidebands") else "" + # TODO: bblite + channel_systs_dict[region]["mcstats"] = {} -def get_year_updown( - templates_dict, sample, region, region_noblinded, blind_str, year, skey, mX_bin=None -): - """ - Return templates with only the given year's shapes shifted up and down by the ``skey`` systematic. - Returns as [up templates, down templates] - """ - updown = [] + for sample_name, card_name in mc_samples.items(): + systs_dict = {} + channel_systs_dict[region][sample_name] = systs_dict - for shift in ["up", "down"]: - sshift = f"{skey}_{shift}" - # get nominal templates for each year - templates = {y: templates_dict[y][region][sample, ...] for y in years} - # replace template for this year with the shifted tempalte - if skey in jecs or skey in jmsr: - # JEC/JMCs saved as different "region" in dict - reg_name = ( - f"{region_noblinded}_{sshift}{blind_str}" - if mX_bin is None - else f"{region}_{sshift}" + # skip signal nuisances in non-signal region + if sample_name in sig_keys and region != "pass": + continue + + # MC stats + mcstats_err = ( + np.sqrt(rates_dict[region][sample_name].variance) + / rates_dict[region][sample_name].value ) - templates[year] = templates_dict[year][reg_name][sample, :, ...] + if mcstats_err > args.epsilon: + channel_systs_dict[region]["mcstats"][sample_name] = 1.0 + mcstats_err + + # rate systematics + for skey, syst in nuisance_params.items(): + if sample_name not in syst.samples or (not pass_region and syst.pass_only): + continue + + val, val_down = syst.value, syst.value_down + if syst.diff_regions: + val = val[region_nosidebands] + val_down = val_down[region_nosidebands] if val_down is not None else val_down + if syst.diff_samples: + val = val[sample_name] + val_down = val_down[sample_name] if val_down is not None else val_down + + systs_dict[skey] = (val, val_down) + + # correlated shape systematics + for skey, syst in corr_year_shape_systs.items(): + if sample_name not in syst.samples or (not pass_region and syst.pass_only): + continue + + if skey in jecs or skey in jmsr: + # JEC/JMCs saved as different "region" in dict + val_up = channels_summed[f"{region_nosidebands}_{skey}_up{sideband_str}"][ + sample_name + ].value + val_down = channels_summed[f"{region_nosidebands}_{skey}_down{sideband_str}"][ + sample_name + ].value + else: + # weight uncertainties saved as different "sample" in dict + val_up = channels_summed[region][f"{sample_name}_{skey}_up"].value + val_down = channels_summed[region][f"{sample_name}_{skey}_down"].value + + srate = rates_dict[region][sample_name].value + systs_dict[skey] = (val_up / srate, val_down / srate) + + # uncorrelated shape systematics + for skey, syst in uncorr_year_shape_systs.items(): + if sample_name not in syst.samples or (not pass_region and syst.pass_only): + continue + + # TODO: figure out why pileup is going crazy + if skey == "pileup": + continue + + for year in years: + if year not in syst.uncorr_years: + continue + + val_up, val_down = get_year_updown( + channels_dict, + sample_name, + region, + region_nosidebands, + sideband_str, + year, + skey, + vbf=True, + ) + + srate = rates_dict[region][sample_name].value + systs_dict[skey] = (val_up / srate, val_down / srate) + + syst_strs = [] + all_processes = sig_keys + bg_keys + [qcd_data_key] + num_ps = len(all_processes) + + # mc stats + for j, channel in enumerate(channels): + cmcstats = channel_systs_dict[channel]["mcstats"] + + # mc stats error too low to add + if cmcstats is None: + continue + + # bblite criteria sastisfied for this channel - single nuisance for all mc samples + elif isinstance(cmcstats, float): + vals = [] + skey = f"mcstats_{channel}" + sstr = f"{skey:<54}lnN " + vals += ["-"] * (j * num_ps) + # add same nuisance for all processes except qcd + vals += [cmcstats] * (num_ps - 1) + ["-"] + vals += ["-"] * ((len(channels) - j - 1) * num_ps) + sstr += join_with_padding(vals) + syst_strs.append(sstr) + + # bblite not satisifed - separate nuisance for all mc samples else: - # weight uncertainties saved as different "sample" in dict - templates[year] = templates_dict[year][region][f"{sample}_{sshift}", ...] + for i, key in enumerate(sig_keys + bg_keys): + if key in cmcstats: + vals = [] + skey = f"mcstats_{channel}_{key}" + sstr = f"{skey:<54}lnN " + # add single nuisance for this sample + vals += ["-"] * ((j * num_ps) + i) + vals += [cmcstats[key]] + vals += ["-"] * ((len(channels) - j - 1) * num_ps + (num_ps - i - 1)) + sstr += join_with_padding(vals) + syst_strs.append(sstr) + + # all other nuisances + for skey in ( + list(nuisance_params.keys()) + + list(corr_year_shape_systs.keys()) + + list(uncorr_year_shape_systs.keys()) + ): + sstr = f"{skey:<54}lnN " + vals = [] + for channel in channels: + for i, key in enumerate(all_processes): + if key == qcd_data_key or skey not in channel_systs_dict[channel][key]: + vals.append("-") + else: + val, val_down = channel_systs_dict[channel][key][skey] + val_str = val if val_down is None else f"{val_down}/{val}" + vals.append(val_str) - if mX_bin is not None: - for year, template in templates.items(): - templates[year] = template[:, mX_bin] + sstr += join_with_padding(vals) + syst_strs.append(sstr) - # sum templates with year's template replaced with shifted - updown.append(sum(list(templates.values())).values()) + syst_str = "\n".join(syst_strs) - return updown + return syst_str + + +def createDatacardABCD(args, templates_dict, templates_summed, shape_vars): + # A, B, C, D (in order) + channels = ["pass", "pass_sidebands", "fail", "fail_sidebands"] # analogous to regions + channels_dict, channels_summed = get_channels(templates_dict, templates_summed, shape_vars[0]) + + datacard_dict, rates_dict = fill_yields(channels, channels_summed) + datacard_dict["systematics"] = get_systematics_abcd( + channels, channels_dict, channels_summed, rates_dict + ) + + out_dir = ( + os.path.join(str(args.cards_dir), args.model_name) + if args.model_name is not None + else args.cards_dir + ) + + with open(f"{out_dir}/datacard.txt", "w") as f: + f.write(helpers.abcd_datacard_template.substitute(datacard_dict)) + + return + + +def main(args): + # templates per region per year, templates per region summed across years + templates_dict, templates_summed = get_templates( + args.templates_dir, years, args.sig_separate, args.scale_templates, args.combine_lasttwo + ) + + # TODO: check if / how to include signal trig eff uncs. (rn only using bg uncs.) + process_systematics(args.templates_dir, args.sig_separate) + + # random template from which to extract shape vars + sample_templates: Hist = templates_summed[list(templates_summed.keys())[0]] + + # [mH(bb)] for nonresonant, [mY, mX] for resonant + shape_vars = [ + ShapeVar(name=axis.name, bins=axis.edges, order=args.nTF[i]) + for i, axis in enumerate(sample_templates.axes[1:]) + ] + + os.makedirs(args.cards_dir, exist_ok=True) + + dc_args = [args, templates_dict, templates_summed, shape_vars] + if args.vbf: + createDatacardABCD(*dc_args) + else: + createDatacardAlphabet(*dc_args) main(args) diff --git a/src/HHbbVV/postprocessing/PostProcess.ipynb b/src/HHbbVV/postprocessing/PostProcess.ipynb index eeca032e..7fd30fbf 100644 --- a/src/HHbbVV/postprocessing/PostProcess.ipynb +++ b/src/HHbbVV/postprocessing/PostProcess.ipynb @@ -77,7 +77,7 @@ "# bdt_data_dir = \"/eos/uscms/store/user/cmantill/bbVV/skimmer/Jun10/bdt_data/\"\n", "year = \"2018\"\n", "\n", - "date = \"23Aug21\"\n", + "date = \"23Oct12\"\n", "plot_dir = f\"../../../plots/PostProcessing/{date}/\"\n", "templates_dir = f\"templates/{date}\"\n", "_ = os.system(f\"mkdir -p {plot_dir}\")\n", @@ -177,6 +177,40 @@ "postprocessing.load_bdt_preds(events_dict, year, bdt_preds_dir, jec_jmsr_shifts=True)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Aside: checking if regressed mass is actually saving useful jets" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bbmsd = utils.get_feat(events_dict[\"HHbbVV\"], \"bbFatJetMsd\", bb_masks[\"HHbbVV\"])\n", + "vvmsd = utils.get_feat(events_dict[\"HHbbVV\"], \"VVFatJetMsd\", bb_masks[\"HHbbVV\"])\n", + "bbpnetm = utils.get_feat(events_dict[\"HHbbVV\"], \"bbFatJetParticleNetMass\", bb_masks[\"HHbbVV\"])\n", + "vvpnetm = utils.get_feat(events_dict[\"HHbbVV\"], \"VVFatJetParticleNetMass\", bb_masks[\"HHbbVV\"])\n", + "bbpnet = utils.get_feat(events_dict[\"HHbbVV\"], \"bbFatJetParticleNetMD_Xbb\", bb_masks[\"HHbbVV\"])\n", + "vvpart = utils.get_feat(events_dict[\"HHbbVV\"], \"VVFatJetParTMD_THWW4q\", bb_masks[\"HHbbVV\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bbreg = (bbmsd < 50) * (bbpnetm > 50)\n", + "plt.title(\"Jets with Msd < 50 GeV and PNet Mass > 50 GeV\")\n", + "plt.hist(bbpnet[bbreg], np.linspace(0, 1, 50), histtype=\"step\", label=\"bb\")\n", + "plt.xlabel(\"Txbb\")\n", + "plt.show()" + ] + }, { "attachments": {}, "cell_type": "markdown", diff --git a/src/HHbbVV/postprocessing/TrainBDT.py b/src/HHbbVV/postprocessing/TrainBDT.py index a8e8deb4..41f9819b 100644 --- a/src/HHbbVV/postprocessing/TrainBDT.py +++ b/src/HHbbVV/postprocessing/TrainBDT.py @@ -195,6 +195,19 @@ def main(args): data_dict = load_data(args.data_path, args.year, args.all_years) + for year, data in data_dict.items(): + for key in training_keys: + print( + ( + f"{year} {key} Yield: " + f'{np.sum(data[data["Dataset"] == key][weight_key])}, ' + "Number of Events: ", + f'{len(data[data["Dataset"] == key])}, ', + ) + ) + + return + if not args.inference_only: training_data_dict = { year: data[ diff --git a/src/HHbbVV/postprocessing/datacardHelpers.py b/src/HHbbVV/postprocessing/datacardHelpers.py new file mode 100644 index 00000000..f82bcd4b --- /dev/null +++ b/src/HHbbVV/postprocessing/datacardHelpers.py @@ -0,0 +1,324 @@ +from typing import Dict, List, Tuple, Union +from dataclasses import dataclass, field +import logging +from string import Template + +import numpy as np +import hist +from hist import Hist +import utils + +from hh_vars import years as all_years + + +################################################# +# Common +################################################# + + +@dataclass +class Syst: + """For storing info about systematics""" + + name: str = None + prior: str = None # e.g. "lnN", "shape", etc. + + # float if same value in all regions/samples, dictionary of values per region/sample if not + # if both, region should be the higher level of the dictionary + value: Union[float, Dict[str, float]] = None + value_down: Union[float, Dict[str, float]] = None # if None assumes symmetric effect + + # if the value is different for different regions or samples + diff_regions: bool = False + diff_samples: bool = False + + samples: List[str] = None # samples affected by it + # in case of uncorrelated unc., which years to split into + uncorr_years: List[str] = field(default_factory=lambda: all_years) + pass_only: bool = False # is it applied only in the pass regions + + def __post_init__(self): + if isinstance(self.value, dict) and not (self.diff_regions or self.diff_samples): + raise RuntimeError( + f"Value for systematic is a dictionary but neither ``diff_regions`` nor ``diff_samples`` is set." + ) + + +@dataclass +class ShapeVar: + """For storing and calculating info about variables used in fit""" + + name: str = None + bins: np.ndarray = None # bin edges + order: int = None # TF order + + def __post_init__(self): + # use bin centers for polynomial fit + self.pts = self.bins[:-1] + 0.5 * np.diff(self.bins) + # scale to be between [0, 1] + self.scaled = (self.pts - self.bins[0]) / (self.bins[-1] - self.bins[0]) + + +def add_bool_arg(parser, name, help, default=False, no_name=None): + """Add a boolean command line argument for argparse""" + varname = "_".join(name.split("-")) # change hyphens to underscores + group = parser.add_mutually_exclusive_group(required=False) + group.add_argument("--" + name, dest=varname, action="store_true", help=help) + if no_name is None: + no_name = "no-" + name + no_help = "don't " + help + else: + no_help = help + group.add_argument("--" + no_name, dest=varname, action="store_false", help=no_help) + parser.set_defaults(**{varname: default}) + + +################################################# +# Template manipulation +################################################# + + +def combine_last_two_bins(templates_dict: Dict, years: List[str]): + for year in years: + for region in templates_dict[year]: + templates_dict[year][region] = utils.rebin_hist( + templates_dict[year][region], + "bbFatJetParticleNetMass", + list(range(50, 240, 10)) + [250], + ) + + +def rem_neg(template_dict: Dict): + for sample, template in template_dict.items(): + template.values()[template.values() < 0] = 0 + + return template_dict + + +def _match_samples(template: Hist, match_template: Hist) -> Hist: + """Temporary solution for case where 2018 templates don't have L1 prefiring in their axis + + Args: + template (Hist): template to which samples may need to be added + match_template (Hist): template from which to extract extra samples + """ + samples = list(match_template.axes[0]) + + # if template already has all the samples, don't do anything + if list(template.axes[0]) == samples: + return template + + # otherwise remake template with samples from ``match_template`` + h = hist.Hist(*match_template.axes, storage="weight") + + for sample in template.axes[0]: + sample_index = np.where(np.array(list(h.axes[0])) == sample)[0][0] + h.view()[sample_index] = template[sample, ...].view() + + return h + + +def sum_templates(template_dict: Dict, years: List[str]): + """Sum templates across years""" + + ttemplate = list(template_dict.values())[0] # sample templates to extract values from + combined = {} + + for region in ttemplate: + thists = [] + + for year in years: + # temporary solution for case where 2018 templates don't have L1 prefiring in their axis + thists.append( + _match_samples(template_dict[year][region], template_dict["2017"][region]) + ) + + combined[region] = sum(thists) + + return combined + + +def combine_templates( + bg_templates: Dict[str, Hist], sig_templates: List[Dict[str, Hist]] +) -> Dict[str, Hist]: + """ + Combines BG and signal templates into a single Hist (per region). + + Args: + bg_templates (Dict[str, Hist]): dictionary of region -> Hist + sig_templates (List[Dict[str, Hist]]): list of dictionaries of region -> Hist for each + signal samples + """ + ctemplates = {} + + for region, bg_template in bg_templates.items(): + # combined sig + bg samples + csamples = list(bg_template.axes[0]) + [ + s for sig_template in sig_templates for s in list(sig_template[region].axes[0]) + ] + + # new hist with all samples + ctemplate = Hist( + hist.axis.StrCategory(csamples, name="Sample"), + *bg_template.axes[1:], + storage="weight", + ) + + # add background hists + for sample in bg_template.axes[0]: + sample_key_index = np.where(np.array(list(ctemplate.axes[0])) == sample)[0][0] + ctemplate.view(flow=True)[sample_key_index, ...] = bg_template[sample, ...].view( + flow=True + ) + + # add signal hists + for st in sig_templates: + sig_template = st[region] + for sample in sig_template.axes[0]: + sample_key_index = np.where(np.array(list(ctemplate.axes[0])) == sample)[0][0] + ctemplate.view(flow=True)[sample_key_index, ...] = sig_template[sample, ...].view( + flow=True + ) + + ctemplates[region] = ctemplate + + return ctemplates + + +################################################# +# Shape Fills +################################################# + + +def _shape_checks(values_up, values_down, values_nominal, effect_up, effect_down, logger): + norm_up = np.sum(values_up) + norm_down = np.sum(values_down) + norm_nominal = np.sum(values_nominal) + prob_up = values_up / norm_up + prob_down = values_down / norm_down + prob_nominal = values_nominal / norm_nominal + shapeEffect_up = np.sum( + np.abs(prob_up - prob_nominal) / (np.abs(prob_up) + np.abs(prob_nominal)) + ) + shapeEffect_down = np.sum( + np.abs(prob_down - prob_nominal) / (np.abs(prob_down) + np.abs(prob_nominal)) + ) + + valid = True + if np.allclose(effect_up, 1.0) and np.allclose(effect_down, 1.0): + valid = False + logger.warning("No shape effect") + elif np.allclose(effect_up, effect_down): + valid = False + logger.warning("Up is the same as Down, but different from nominal") + elif np.allclose(effect_up, 1.0) or np.allclose(effect_down, 1.0): + valid = False + logger.warning("Up or Down is the same as nominal (one-sided)") + elif shapeEffect_up < 0.001 and shapeEffect_down < 0.001: + valid = False + logger.warning("No genuine shape effect (just norm)") + elif (norm_up > norm_nominal and norm_down > norm_nominal) or ( + norm_up < norm_nominal and norm_down < norm_nominal + ): + valid = False + logger.warning("Up and Down vary norm in the same direction") + + +def get_effect_updown(values_nominal, values_up, values_down, mask, logger): + effect_up = np.ones_like(values_nominal) + effect_down = np.ones_like(values_nominal) + + mask_up = mask & (values_up >= 0) + mask_down = mask & (values_down >= 0) + + effect_up[mask_up] = values_up[mask_up] / values_nominal[mask_up] + effect_down[mask_down] = values_down[mask_down] / values_nominal[mask_down] + + zero_up = values_up == 0 + zero_down = values_down == 0 + + effect_up[mask_up & zero_up] = values_nominal[mask_up & zero_up] * args.epsilon + effect_down[mask_down & zero_down] = values_nominal[mask_down & zero_down] * args.epsilon + + _shape_checks(values_up, values_down, values_nominal, effect_up, effect_down, logger) + + logging.debug("nominal : {nominal}".format(nominal=values_nominal)) + logging.debug("effect_up : {effect_up}".format(effect_up=effect_up)) + logging.debug("effect_down: {effect_down}".format(effect_down=effect_down)) + + return effect_up, effect_down + + +################################################# +# VBF +################################################# + + +abcd_datacard_template = Template( + """ +imax $num_bins +jmax $num_bgs +kmax * +--------------- +bin $bins +observation $observations +------------------------------ +bin $bins_x_processes +process $processes_per_bin +process $processes_index +rate $processes_rates +------------------------------ +$systematics +single_A rateParam $binA $qcdlabel (@0*@2/@1) single_B,single_C,single_D +single_B rateParam $binB $qcdlabel $dataqcdB +single_C rateParam $binC $qcdlabel $dataqcdC +single_D rateParam $binD $qcdlabel $dataqcdD +""" +) + + +def join_with_padding(items, padding: int = 40): + """Joins items into a string, padded by ``padding`` spaces""" + ret = f"{{:<{padding}}}" * len(items) + return ret.format(*items) + + +def rebin_channels(templates: Dict, axis_name: str, mass_window: List[float]): + """Convert templates into single-bin hists per ABCD channel""" + rebinning_edges = [50] + mass_window + [250] # [0] [1e4] + + channels = {} + for region, template in templates.items(): + # rebin into [left sideband, mass window, right window] + rebinned = utils.rebin_hist(template, axis_name, rebinning_edges) + channels[region] = rebinned[..., 1] + # sum sidebands + channels[f"{region}_sidebands"] = rebinned[..., 0] + rebinned[..., 2] + + return channels + + +def get_channels(templates_dict: Dict, templates_summed: Dict, shape_var: ShapeVar): + """Convert templates into single-bin hists per ABCD channel""" + + mass_window = [100, 150] + axis_name = shape_var.name + + channels_summed = rebin_channels(templates_summed, axis_name, mass_window) + channels_dict = {} + for key, templates in templates_dict.items(): + channels_dict[key] = rebin_channels(templates, axis_name, mass_window) + + return channels_dict, channels_summed + + +################################################# +# Resonant +################################################# + + +def mxmy(sample): + mY = int(sample.split("-")[-1]) + mX = int(sample.split("NMSSM_XToYHTo2W2BTo4Q2B_MX-")[1].split("_")[0]) + + return (mX, mY) diff --git a/src/HHbbVV/postprocessing/hh_vars.py b/src/HHbbVV/postprocessing/hh_vars.py index b5e545e9..bb823e5c 100644 --- a/src/HHbbVV/postprocessing/hh_vars.py +++ b/src/HHbbVV/postprocessing/hh_vars.py @@ -230,6 +230,8 @@ "VVFatJetPtOverDijetPt", "VVFatJetPtOverbbFatJetPt", "BDTScore", + "VBFJetPt", + "vbf_Mass_jj", ] diff --git a/src/HHbbVV/postprocessing/templates/23Oct3VBF/2017_templates.pkl b/src/HHbbVV/postprocessing/templates/23Oct3VBF/2017_templates.pkl new file mode 100644 index 00000000..938af030 Binary files /dev/null and b/src/HHbbVV/postprocessing/templates/23Oct3VBF/2017_templates.pkl differ diff --git a/src/HHbbVV/postprocessing/templates/23Oct3VBF/lpsfs.csv b/src/HHbbVV/postprocessing/templates/23Oct3VBF/lpsfs.csv new file mode 100644 index 00000000..19d688c3 --- /dev/null +++ b/src/HHbbVV/postprocessing/templates/23Oct3VBF/lpsfs.csv @@ -0,0 +1,12 @@ +,SF,syst_unc,stat_unc,sj_pt_unc,sj_matching_unc +HHbbVV,0.93 ± 0.21,0.08735995513326773,0.02178232885498517,0.0061149325206791465,0.21056243132864816 +ggHH_kl_2p45_kt_1_HHbbVV,0.94 ± 0.22,0.09825588332373836,0.017728326791750713,0.007576634540588289,0.2109042985659727 +ggHH_kl_5_kt_1_HHbbVV,0.90 ± 0.31,0.26670664011549117,0.029146162489653592,0.009392097264437689,0.21495771560903995 +ggHH_kl_0_kt_1_HHbbVV,0.97 ± 0.22,0.0855042203073505,0.019319390242510937,0.0057113506172409585,0.21249833138903884 +VBFHHbbVV,0.97 ± 0.29,0.21606146805271081,0.022866755372481395,0.011791744840525328,0.19970070579826676 +qqHH_CV_1_C2V_0_kl_1_HHbbVV,0.95 ± 0.18,0.0440400535017882,0.016305841465033288,0.05912205378108328,0.1762753439506925 +qqHH_CV_1p5_C2V_1_kl_1_HHbbVV,0.97 ± 0.19,0.043170552991571264,0.01843468279444452,0.05548925971951003,0.17659871392784582 +qqHH_CV_1_C2V_1_kl_2_HHbbVV,0.93 ± 0.41,0.3741327061402456,0.0365511621336087,0.015097514340344167,0.2173996175908222 +qqHH_CV_1_C2V_2_kl_1_HHbbVV,0.96 ± 0.18,0.04751934156410037,0.016217262585498455,0.06330571028260944,0.17501290974191908 +qqHH_CV_1_C2V_1_kl_0_HHbbVV,0.97 ± 0.22,0.12492900211693345,0.017097719425056068,0.00947242206235012,0.1944825092307107 +qqHH_CV_0p5_C2V_1_kl_1_HHbbVV,0.95 ± 0.18,0.043193276036633484,0.01618427553643857,0.06082940744062617,0.1750516645939186 diff --git a/src/HHbbVV/postprocessing/templates/23Oct3VBF/systematics.json b/src/HHbbVV/postprocessing/templates/23Oct3VBF/systematics.json new file mode 100644 index 00000000..f23fe994 --- /dev/null +++ b/src/HHbbVV/postprocessing/templates/23Oct3VBF/systematics.json @@ -0,0 +1,122 @@ +{ + "2017": { + "pass": { + "trig_total": 87.6234207605359, + "trig_total_err": 1.232538971320039 + }, + "fail": { + "trig_total": 346.5607033736205, + "trig_total_err": 0.7244257707820206 + } + }, + "HHbbVV": { + "lp_sf": 0.9255566155790668, + "lp_sf_unc": 0.2290854894968195, + "lp_sf_uncs": { + "syst_unc": 0.08735995513326773, + "stat_unc": 0.02178232885498517, + "sj_pt_unc": 0.0061149325206791465, + "sj_matching_unc": 0.21056243132864816 + } + }, + "ggHH_kl_2p45_kt_1_HHbbVV": { + "lp_sf": 0.9391472823461275, + "lp_sf_unc": 0.23346635886810285, + "lp_sf_uncs": { + "syst_unc": 0.09825588332373836, + "stat_unc": 0.017728326791750713, + "sj_pt_unc": 0.007576634540588289, + "sj_matching_unc": 0.2109042985659727 + } + }, + "ggHH_kl_5_kt_1_HHbbVV": { + "lp_sf": 0.8985676068767847, + "lp_sf_unc": 0.34391417775434757, + "lp_sf_uncs": { + "syst_unc": 0.26670664011549117, + "stat_unc": 0.029146162489653592, + "sj_pt_unc": 0.009392097264437689, + "sj_matching_unc": 0.21495771560903995 + } + }, + "ggHH_kl_0_kt_1_HHbbVV": { + "lp_sf": 0.9670162916725116, + "lp_sf_unc": 0.22993992889167636, + "lp_sf_uncs": { + "syst_unc": 0.0855042203073505, + "stat_unc": 0.019319390242510937, + "sj_pt_unc": 0.0057113506172409585, + "sj_matching_unc": 0.21249833138903884 + } + }, + "VBFHHbbVV": { + "lp_sf": 0.9679521376909226, + "lp_sf_unc": 0.2953385576267809, + "lp_sf_uncs": { + "syst_unc": 0.21606146805271081, + "stat_unc": 0.022866755372481395, + "sj_pt_unc": 0.011791744840525328, + "sj_matching_unc": 0.19970070579826676 + } + }, + "qqHH_CV_1_C2V_0_kl_1_HHbbVV": { + "lp_sf": 0.9481126483385646, + "lp_sf_unc": 0.19176501481383776, + "lp_sf_uncs": { + "syst_unc": 0.0440400535017882, + "stat_unc": 0.016305841465033288, + "sj_pt_unc": 0.05912205378108328, + "sj_matching_unc": 0.1762753439506925 + } + }, + "qqHH_CV_1p5_C2V_1_kl_1_HHbbVV": { + "lp_sf": 0.971884691277887, + "lp_sf_unc": 0.19097041100787895, + "lp_sf_uncs": { + "syst_unc": 0.043170552991571264, + "stat_unc": 0.01843468279444452, + "sj_pt_unc": 0.05548925971951003, + "sj_matching_unc": 0.17659871392784582 + } + }, + "qqHH_CV_1_C2V_1_kl_2_HHbbVV": { + "lp_sf": 0.93214579437633, + "lp_sf_unc": 0.43451328854827137, + "lp_sf_uncs": { + "syst_unc": 0.3741327061402456, + "stat_unc": 0.0365511621336087, + "sj_pt_unc": 0.015097514340344167, + "sj_matching_unc": 0.2173996175908222 + } + }, + "qqHH_CV_1_C2V_2_kl_1_HHbbVV": { + "lp_sf": 0.9569683205759278, + "lp_sf_unc": 0.19276467248739185, + "lp_sf_uncs": { + "syst_unc": 0.04751934156410037, + "stat_unc": 0.016217262585498455, + "sj_pt_unc": 0.06330571028260944, + "sj_matching_unc": 0.17501290974191908 + } + }, + "qqHH_CV_1_C2V_1_kl_0_HHbbVV": { + "lp_sf": 0.9679770635328802, + "lp_sf_unc": 0.23197577622646612, + "lp_sf_uncs": { + "syst_unc": 0.12492900211693345, + "stat_unc": 0.017097719425056068, + "sj_pt_unc": 0.00947242206235012, + "sj_matching_unc": 0.1944825092307107 + } + }, + "qqHH_CV_0p5_C2V_1_kl_1_HHbbVV": { + "lp_sf": 0.9529558040722325, + "lp_sf_unc": 0.1909735373189062, + "lp_sf_uncs": { + "syst_unc": 0.043193276036633484, + "stat_unc": 0.01618427553643857, + "sj_pt_unc": 0.06082940744062617, + "sj_matching_unc": 0.1750516645939186 + } + } +} \ No newline at end of file