diff --git a/scripts/combine/pullsAndImpacts.py b/scripts/combine/pullsAndImpacts.py index 7b34912ef..34811e7eb 100644 --- a/scripts/combine/pullsAndImpacts.py +++ b/scripts/combine/pullsAndImpacts.py @@ -310,6 +310,7 @@ def parseArgs(): parser.add_argument("--filters", nargs="*", type=str, help="Filter regexes to select nuisances by name") parser.add_argument("--grouping", type=str, default=None, help="Select nuisances by a predefined grouping", choices=groupings.keys()) parser.add_argument("-t","--translate", type=str, default=None, help="Specify .json file to translate labels") + parser.add_argument("--noImpacts", action='store_true', help="Don't show impacts") parsers = parser.add_subparsers(dest='output_mode') interactive = parsers.add_parser("interactive", help="Launch and interactive dash session") interactive.add_argument("-i", "--interface", default="localhost", help="The network interface to bind to.") @@ -319,7 +320,6 @@ def parseArgs(): output.add_argument("--otherExtensions", default=[], type=str, nargs="*", help="Additional output file types to write") output.add_argument("-n", "--num", type=int, help="Number of nuisances to plot") output.add_argument("--noPulls", action='store_true', help="Don't show pulls (not defined for groups)") - output.add_argument("--noImpacts", action='store_true', help="Don't show impacts") output.add_argument("--eoscp", action='store_true', help="Use of xrdcp for eos output rather than the mount") return parser.parse_args() diff --git a/scripts/combine/setupCombine.py b/scripts/combine/setupCombine.py index 1fa836405..35a781ecb 100644 --- a/scripts/combine/setupCombine.py +++ b/scripts/combine/setupCombine.py @@ -87,8 +87,7 @@ def make_parser(parser=None): parser.add_argument("--resumUnc", default="tnp", type=str, choices=["scale", "tnp", "tnp_minnlo", "minnlo", "none"], help="Include SCETlib uncertainties") parser.add_argument("--noTransitionUnc", action="store_true", help="Do not include matching transition parameter variations.") parser.add_argument("--npUnc", default="Delta_Lambda", type=str, choices=combine_theory_helper.TheoryHelper.valid_np_models, help="Nonperturbative uncertainty model") - parser.add_argument("--tnpMagnitude", default=1, type=float, help="Variation size for the TNP") - parser.add_argument("--scaleTNP", default=5, type=float, help="Scale the TNP uncertainties by this factor") + parser.add_argument("--scaleTNP", default=1, type=float, help="Scale the TNP uncertainties by this factor") parser.add_argument("--scalePdf", default=1, type=float, help="Scale the PDF hessian uncertainties by this factor") parser.add_argument("--pdfUncFromCorr", action='store_true', help="Take PDF uncertainty from correction hist (Requires having run that correction)") parser.add_argument("--massVariation", type=float, default=100, help="Variation of boson mass") @@ -158,14 +157,22 @@ def setup(args, inputFile, fitvar, xnorm=False): if not xnorm and (args.axlim or args.rebin or args.absval): datagroups.set_rebin_action(fitvar, args.axlim, args.rebin, args.absval, args.rebinBeforeSelection, rename=False) - wmass = datagroups.mode in ["wmass", "lowpu_w"] + wmass = datagroups.mode in ["wmass", "lowpu_w"] wlike = datagroups.mode == "wlike" lowPU = "lowpu" in datagroups.mode # Detect lowpu dilepton dilepton = "dilepton" in datagroups.mode or any(x in ["ptll", "mll"] for x in fitvar) + genfit = datagroups.mode == "vgen" + + if genfit: + hasw = any("W" in x for x in args.filterProcGroups) + hasz = any("Z" in x for x in args.filterProcGroups) + if hasw and hasz: + raise ValueError("Only W or Z processes are permitted in the gen fit") + wmass = hasw simultaneousABCD = wmass and args.simultaneousABCD and not xnorm - constrainMass = args.forceConstrainMass or args.fitXsec or (dilepton and not "mll" in fitvar) + constrainMass = args.forceConstrainMass or args.fitXsec or (dilepton and not "mll" in fitvar) or genfit logger.debug(f"constrainMass = {constrainMass}") if wmass: @@ -303,8 +310,8 @@ def setup(args, inputFile, fitvar, xnorm=False): label = 'W' if wmass else 'Z' cardTool.setCustomSystGroupMapping({ "theoryTNP" : f".*resum.*|.*TNP.*|mass.*{label}.*", - "resumTheory" : f".*resum.*|.*TNP.*|mass.*{label}.*", - "allTheory" : f"pdf.*|.*QCD.*|.*resum.*|.*TNP.*|mass.*{label}.*", + "resumTheory" : f".*scetlib.*|.*resum.*|.*TNP.*|mass.*{label}.*", + "allTheory" : f".*scetlib.*|pdf.*|.*QCD.*|.*resum.*|.*TNP.*|mass.*{label}.*", "ptTheory" : f".*QCD.*|.*resum.*|.*TNP.*|mass.*{label}.*", }) cardTool.setCustomSystForCard(args.excludeNuisances, args.keepNuisances) @@ -639,9 +646,8 @@ def assertSample(name, startsWith=["W", "Z"], excludeMatch=[]): transitionUnc = not args.noTransitionUnc, propagate_to_fakes=to_fakes, np_model=args.npUnc, - tnp_magnitude=args.tnpMagnitude, tnp_scale = args.scaleTNP, - mirror_tnp=True, + mirror_tnp=False, pdf_from_corr=args.pdfUncFromCorr, scale_pdf_unc=args.scalePdf, minnlo_unc=args.minnloScaleUnc, @@ -657,7 +663,9 @@ def assertSample(name, startsWith=["W", "Z"], excludeMatch=[]): theory_helper.add_all_theory_unc(theorySystSamples, skipFromSignal=args.noPDFandQCDtheorySystOnSignal) - if xnorm or datagroups.mode == "vgen": + if xnorm or genfit: + if genfit: + cardTool.addLnNSystematic("lumi", processes=["signal_samples"], size=1.012, group="luminosity") return cardTool # Below: experimental uncertainties @@ -1007,11 +1015,13 @@ def main(args, xnorm=False): writer.set_fitresult(args.fitresult, mc_stat=not args.noMCStat) if len(outnames) == 1: - outfile, outfolder = outnames[0] + outfolder, outfile = outnames[0] else: - outfile, outfolder = f"{args.outfolder}/Combination{'_statOnly' if args.doStatOnly else ''}{'_'+args.postfix if args.postfix else ''}/", "Combination" + dir_append = '_'.join(['', *filter(lambda x: x, ['statOnly' if args.doStatOnly else '', args.postfix])]) + outfolder = f"{args.outfolder}/Combination_{''.join([o[1] for o in outnames])}{dir_append}/" + outfile = "Combination" logger.info(f"Writing HDF5 output to {outfile}") - writer.write(args, outfile, outfolder) + writer.write(args, outfolder, outfile) else: if len(args.inputFile) > 1: raise IOError(f"Multiple input files only supported within --hdf5 mode") diff --git a/scripts/corrections/make_ptv_corr.py b/scripts/corrections/make_ptv_corr.py index 0187339a8..d15ba058a 100644 --- a/scripts/corrections/make_ptv_corr.py +++ b/scripts/corrections/make_ptv_corr.py @@ -67,6 +67,6 @@ } meta_dict = input_tools.get_metadata(args.inputFile) -outfile = args.outpath+"/dataPtll" +outfile = args.outpath+"/dataRecoPtll" output_tools.write_theory_corr_hist(outfile, args.proc.upper(), output_dict, args, meta_dict) logger.info(f"Average correction is {np.mean(corrh.values())}") diff --git a/scripts/histmakers/mz_dilepton.py b/scripts/histmakers/mz_dilepton.py index 2f7c6d469..ba1d24f6f 100644 --- a/scripts/histmakers/mz_dilepton.py +++ b/scripts/histmakers/mz_dilepton.py @@ -22,8 +22,6 @@ parser.add_argument("--useDileptonTriggerSelection", action='store_true', help="Use dilepton trigger selection (default uses the Wlike one, with one triggering muon and odd/even event selection to define its charge, staying agnostic to the other)") parser.add_argument("--noAuxiliaryHistograms", action="store_true", help="Remove auxiliary histograms to save memory (removed by default with --unfolding or --theoryAgnostic)") -parser = common.set_parser_default(parser, "pt", [34,26.,60.]) -parser = common.set_parser_default(parser, "eta", [48,-2.4,2.4]) parser = common.set_parser_default(parser, "aggregateGroups", ["Diboson", "Top", "Wtaunu", "Wmunu"]) parser = common.set_parser_default(parser, "ewTheoryCorr", ["virtual_ew", "pythiaew_ISR", "horaceqedew_FSR", "horacelophotosmecoffew_FSR",]) parser = common.set_parser_default(parser, "excludeProcs", ["QCD"]) @@ -102,7 +100,7 @@ nominal_axes = [all_axes[a] for a in nominal_cols] if isUnfolding: - unfolding_axes, unfolding_cols, unfolding_selections = differential.get_dilepton_axes(args.genAxes, common.get_gen_axes(isPoiAsNoi, dilepton_ptV_binning), add_out_of_acceptance_axis=isPoiAsNoi) + unfolding_axes, unfolding_cols, unfolding_selections = differential.get_dilepton_axes(args.genAxes, common.get_gen_axes(isPoiAsNoi or inclusive, dilepton_ptV_binning), add_out_of_acceptance_axis=isPoiAsNoi) if not isPoiAsNoi: datasets = unfolding_tools.add_out_of_acceptance(datasets, group = "Zmumu") @@ -173,16 +171,16 @@ def build_graph(df, dataset): cols = nominal_cols if isUnfolding and dataset.name == "ZmumuPostVFP": - df = unfolding_tools.define_gen_level(df, args.genLevel, dataset.name, mode="dilepton") + fidmode = "mz_window" if inclusive else "mz_dilepton" + df = unfolding_tools.define_gen_level(df, args.genLevel, dataset.name, mode=fidmode) + fidargs = unfolding_tools.get_fiducial_args(fidmode, pt_min=args.pt[1], pt_max=args.pt[2]) if hasattr(dataset, "out_of_acceptance"): logger.debug("Reject events in fiducial phase space") - df = unfolding_tools.select_fiducial_space(df, mode="dilepton", pt_min=args.pt[1], pt_max=args.pt[2], - mass_min=mass_min, mass_max=mass_max, selections=unfolding_selections, accept=False) + df = unfolding_tools.select_fiducial_space(df, mode=fidmode, selections=[], accept=False, **fidargs) else: logger.debug("Select events in fiducial phase space") - df = unfolding_tools.select_fiducial_space(df, mode="dilepton", pt_min=args.pt[1], pt_max=args.pt[2], - mass_min=mass_min, mass_max=mass_max, selections=unfolding_selections, select=not isPoiAsNoi, accept=True) + df = unfolding_tools.select_fiducial_space(df, mode="mz_dilepton", selections=[], select=not isPoiAsNoi, accept=True, **fidargs) unfolding_tools.add_xnorm_histograms(results, df, args, dataset.name, corr_helpers, qcdScaleByHelicity_helper, unfolding_axes, unfolding_cols) if not isPoiAsNoi: axes = [*nominal_axes, *unfolding_axes] @@ -305,7 +303,8 @@ def build_graph(df, dataset): logger.debug(f"Creating special histogram '{noiAsPoiHistName}' for unfolding to treat POIs as NOIs") results.append(df.HistoBoost(noiAsPoiHistName, [*nominal_axes, *unfolding_axes], [*nominal_cols, *unfolding_cols, "nominal_weight"])) - for obs in ["ptll", "mll", "yll", "cosThetaStarll", "phiStarll", "etaPlus", "etaMinus", "ptPlus", "ptMinus"]: + #for obs in ["ptll", "mll", "yll", "cosThetaStarll", "phiStarll", "etaPlus", "etaMinus", "ptPlus", "ptMinus"]: + for obs in ["ptll", "mll", "yll", "etaPlus", "etaMinus", "ptPlus", "ptMinus"]: if dataset.is_data: results.append(df.HistoBoost(f"nominal_{obs}", [all_axes[obs]], [obs])) else: diff --git a/scripts/histmakers/mz_wlike_with_mu_eta_pt.py b/scripts/histmakers/mz_wlike_with_mu_eta_pt.py index af26b356e..d4f7c9c45 100644 --- a/scripts/histmakers/mz_wlike_with_mu_eta_pt.py +++ b/scripts/histmakers/mz_wlike_with_mu_eta_pt.py @@ -146,16 +146,16 @@ def build_graph(df, dataset): cols = nominal_cols if isUnfolding and isZ: - df = unfolding_tools.define_gen_level(df, args.genLevel, dataset.name, mode="wlike") + fidmode = "mz_wlike_inclusive" if args.inclusive else "mz_wlike" + df = unfolding_tools.define_gen_level(df, args.genLevel, dataset.name, mode=fidmode) + fidargs = unfolding_tools.get_fiducial_args(fidmode, pt_min=args.pt[1], pt_max=pt_min[2]) if hasattr(dataset, "out_of_acceptance"): logger.debug("Reject events in fiducial phase space") - df = unfolding_tools.select_fiducial_space(df, - mode="wlike", pt_min=args.pt[1], pt_max=args.pt[2], mass_min=mass_min, mass_max=mass_max, mtw_min=mtw_min, accept=False) + df = unfolding_tools.select_fiducial_space(df, mode=fidmode, accept=False, **fidargs) else: logger.debug("Select events in fiducial phase space") - df = unfolding_tools.select_fiducial_space(df, - mode="wlike", pt_min=args.pt[1], pt_max=args.pt[2], mass_min=mass_min, mass_max=mass_max, mtw_min=mtw_min, accept=True) + df = unfolding_tools.select_fiducial_space(df, mode=fidmode, accept=True, **fidargs) unfolding_tools.add_xnorm_histograms(results, df, args, dataset.name, corr_helpers, qcdScaleByHelicity_helper, unfolding_axes, unfolding_cols) axes = [*nominal_axes, *unfolding_axes] diff --git a/scripts/histmakers/w_z_gen_dists.py b/scripts/histmakers/w_z_gen_dists.py index c038a141e..82b9fcea2 100644 --- a/scripts/histmakers/w_z_gen_dists.py +++ b/scripts/histmakers/w_z_gen_dists.py @@ -5,7 +5,7 @@ import narf import wremnants -from wremnants import theory_tools,syst_tools,theory_corrections +from wremnants import theory_tools,syst_tools,theory_corrections,unfolding_tools from wremnants.datasets.dataset_tools import getDatasets import hist import math @@ -21,13 +21,16 @@ parser.add_argument("--photonHists", action='store_true', help="Also store photon kinematics") parser.add_argument("--skipEWHists", action='store_true', help="Also store histograms for EW reweighting. Use with --filter horace") parser.add_argument("--signedY", action='store_true', help="use signed Y") -parser.add_argument("--applySelection", action='store_true', help="Apply selection on leptons") +parser.add_argument("--fiducial", choices=["inclusive", "masswindow", "dilepton", "singlelep"], help="Apply selection on leptons") parser.add_argument("--auxiliaryHistograms", action="store_true", help="Safe auxiliary histograms (mainly for ew analysis)") parser.add_argument("--ptqVgen", action='store_true', help="To store qt by Q variable instead of ptVgen, GEN only ", default=None) +parser.add_argument("--helicity", action='store_true', help="Make qcdScaleByHelicity hist") parser = common.set_parser_default(parser, "filterProcs", common.vprocs) parser = common.set_parser_default(parser, "theoryCorr", []) parser = common.set_parser_default(parser, "ewTheoryCorr", []) +parser = common.set_parser_default(parser, "pt", [34,26.,60.]) +parser = common.set_parser_default(parser, "eta", [48,-2.4,2.4]) args = parser.parse_args() @@ -66,7 +69,8 @@ if not args.useTheoryAgnosticBinning: axis_ptVgen = hist.axis.Variable( - list(range(0,151))+[160., 190.0, 220.0, 250.0, 300.0, 400.0, 500.0, 800.0, 13000.0], + #list(range(0,151))+[160., 190.0, 220.0, 250.0, 300.0, 400.0, 500.0, 800.0, 13000.0], + [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 17, 20, 23, 27, 32, 40, 54, 100], #common.ptV_binning, name = "ptVgen", underflow=False, ) @@ -90,8 +94,8 @@ 0, 1, name="chargeVgen", underflow=False, overflow=False ) -axis_l_eta_gen = hist.axis.Regular(48, -2.4, 2.4, name = "eta") -axis_l_pt_gen = hist.axis.Regular(29, 26., 55., name = "pt") +axis_abseta_gen = hist.axis.Regular(24, 0, 2.4, name = "abseta") +axis_l_pt_gen = hist.axis.Regular(34, 26., 60., name = "pt") # dilepton gen axes used in unfolding gen_axes = common.get_gen_axes(flow=False) @@ -127,20 +131,21 @@ def build_graph(df, dataset): if isZ: nominal_axes = [axis_massZgen, axis_rapidity, axis_ptqVgen if args.ptqVgen else axis_ptVgen, axis_chargeZgen] - lep_axes = [axis_l_eta_gen, axis_l_pt_gen, axis_chargeZgen] + lep_axes = [axis_abseta_gen, axis_l_pt_gen, axis_chargeZgen] else: nominal_axes = [axis_massWgen, axis_rapidity, axis_ptqVgen if args.ptqVgen else axis_ptVgen, axis_chargeWgen] - lep_axes = [axis_l_eta_gen, axis_l_pt_gen, axis_chargeWgen] + lep_axes = [axis_abseta_gen, axis_l_pt_gen, axis_chargeWgen] nominal_cols = ["massVgen", col_rapidity, "ptqVgen" if args.ptqVgen else "ptVgen", "chargeVgen"] - lep_cols = ["etaPrefsrLep", "ptPrefsrLep", "chargeVgen"] + lep_cols = ["absEtaGen", "ptGen", "chargeVgen"] + + if args.fiducial is not None: + fidmode = f"{'mz' if isZ else 'mw'}_{args.fiducial if args.fiducial else 'inclusive'}" + df = unfolding_tools.define_gen_level(df, "preFSR", dataset.name, mode=fidmode) + + fidargs = unfolding_tools.get_fiducial_args(fidmode) + df = unfolding_tools.select_fiducial_space(df, mode=fidmode, **fidargs) if args.singleLeptonHists and (isW or isZ): - if isW: - df = df.Define('ptPrefsrLep', 'genl.pt()') - df = df.Define('etaPrefsrLep', 'genl.eta()') - else: - df = df.Define('ptPrefsrLep', 'genlanti.pt()') - df = df.Define('etaPrefsrLep', 'genlanti.eta()') results.append(df.HistoBoost("nominal_genlep", lep_axes, [*lep_cols, "nominal_weight"], storage=hist.storage.Weight())) if not args.skipEWHists and (isW or isZ): @@ -292,7 +297,10 @@ def build_graph(df, dataset): if 'horace' not in dataset.name and 'winhac' not in dataset.name and \ "LHEScaleWeight" in df.GetColumnNames() and "LHEPdfWeight" in df.GetColumnNames(): - df = syst_tools.add_theory_hists(results, df, args, dataset.name, corr_helpers, None, nominal_axes, nominal_cols, base_name="nominal_gen") + + qcdScaleByHelicity_helper = wremnants.theory_corrections.make_qcd_uncertainty_helper_by_helicity(is_w_like = dataset.name[0] != "W") if args.helicity else None + + df = syst_tools.add_theory_hists(results, df, args, dataset.name, corr_helpers, qcdScaleByHelicity_helper, nominal_axes, nominal_cols, base_name="nominal_gen") return results, weightsum diff --git a/scripts/utilities/fitresult_pois_to_hist.py b/scripts/utilities/fitresult_pois_to_hist.py index f3cdab3c7..c2b0ac609 100644 --- a/scripts/utilities/fitresult_pois_to_hist.py +++ b/scripts/utilities/fitresult_pois_to_hist.py @@ -1,5 +1,6 @@ import h5py import os +import narf from utilities import common, logging from utilities.io_tools.conversion_tools import fitresult_pois_to_hist @@ -38,21 +39,22 @@ logger.info(f"Creating output folder {args.outfolder}") os.makedirs(args.outfolder) +res_dict = { + "results" : result, + "meta_info" : narf.ioutils.make_meta_info_dict(args=args, wd=common.base_dir), +} +if meta is not None: + res_dict["combine_meta"] = meta +if meta_exp is not None: + res_dict["combine_meta_exp"] = meta_exp + if args.h5py: from narf import ioutils with h5py.File(outfile, "w") as f: logger.debug(f"Pickle and dump results") - ioutils.pickle_dump_h5py("results", result, f) - if meta is not None: - ioutils.pickle_dump_h5py("meta", meta, f) - if meta_exp is not None: - ioutils.pickle_dump_h5py("meta_exp", meta_exp, f) + for k,v in res_dict.items(): + ioutils.pickle_dump_h5py(k, v, f) else: import pickle with open(outfile, "wb") as f: - meta_info = {} - if meta is not None: - meta_info["meta"] = meta - if meta_exp is not None: - meta_info["meta_exp"] = meta_exp - pickle.dump({"results": result, **meta_info}, f) + pickle.dump(res_dict, f) diff --git a/utilities/common.py b/utilities/common.py index 42c1d3f94..432888162 100644 --- a/utilities/common.py +++ b/utilities/common.py @@ -173,6 +173,7 @@ def set_subparsers(subparser, name): help="Generator level definition for unfolding") subparser.add_argument("--genBins", type=int, nargs="+", default=[16, 0], help="Number of generator level bins") + subparser.add_argument("--inclusive", action='store_true', help="No fiducial selection (mass window only)") elif "theoryAgnostic" in name: # specific for theory agnostic subparser.add_argument("--genAxes", type=str, nargs="+", default=["ptVgenSig", "absYVgenSig", "helicitySig"], choices=["qGen", "ptVgenSig", "absYVgenSig", "helicitySig"], help="Generator level variable") diff --git a/utilities/io_tools/conversion_tools.py b/utilities/io_tools/conversion_tools.py index 5baeb477f..95667c63d 100644 --- a/utilities/io_tools/conversion_tools.py +++ b/utilities/io_tools/conversion_tools.py @@ -116,6 +116,9 @@ def fitresult_pois_to_hist(infile, result=None, poi_types = ["mu", "pmaskedexp", axes_names = [a.name for a in axes] data = combinetf_input.select_pois(df, axes_names, base_processes=proc, flow=True) + logger.debug(f"-----> Poi {poi_type}") + logger.debug(f"value = {data["value"].values}") + logger.debug(f"Scale = {scale}") values = np.reshape(data["value"].values/channel_scale, shape) variances = np.reshape( (data["err_total"].values/channel_scale)**2, shape) @@ -153,4 +156,4 @@ def fitresult_pois_to_hist(infile, result=None, poi_types = ["mu", "pmaskedexp", logger.warning(f"Histogram {hist_name_syst} already in result, it will be overridden") result[poi_key][channel][proc][hist_name_syst] = h_syst - return result, meta \ No newline at end of file + return result, meta diff --git a/wremnants-data b/wremnants-data index 392211c86..b3691fdd4 160000 --- a/wremnants-data +++ b/wremnants-data @@ -1 +1 @@ -Subproject commit 392211c86f5f2f0f572e1d4b398b230891d05144 +Subproject commit b3691fdd43aa622eff6e8b26ebcfde6aaba7c0ee diff --git a/wremnants/CardTool.py b/wremnants/CardTool.py index da73143e0..7dc922828 100644 --- a/wremnants/CardTool.py +++ b/wremnants/CardTool.py @@ -546,6 +546,7 @@ def systHists(self, hvar, syst, hnom): if "systNameReplace" in systInfo and systInfo["systNameReplace"]: for rep in systInfo["systNameReplace"]: name = name.replace(*rep) + logger.debug(f"Replacement {rep} yields new name {name}") if name and "systNamePrepend" in systInfo and systInfo["systNamePrepend"]: name = systInfo["systNamePrepend"]+name # Obviously there is a nicer way to do this... diff --git a/wremnants/combine_helpers.py b/wremnants/combine_helpers.py index 595e68eae..4bb7cd579 100644 --- a/wremnants/combine_helpers.py +++ b/wremnants/combine_helpers.py @@ -48,7 +48,7 @@ def add_electroweak_uncertainty(card_tool, ewUncs, flavor="mu", samples="single_ for ewUnc in ewUncs: if ewUnc == "default": if z_samples: - if mode in ["wmass", "wlike", "lowpu_w", "vgen"]: + if mode in ["wmass", "wlike", "lowpu_w"]: ewUnc = "virtual_ew_wlike" else: ewUnc = "virtual_ew" diff --git a/wremnants/combine_theory_helper.py b/wremnants/combine_theory_helper.py index 7d3c40e45..220e58233 100644 --- a/wremnants/combine_theory_helper.py +++ b/wremnants/combine_theory_helper.py @@ -26,7 +26,6 @@ def __init__(self, card_tool, hasNonsigSamples=False): self.np_model = "Delta_Lambda" self.pdf_from_corr = False self.scale_pdf_unc = 1. - self.tnp_magnitude = 1. self.mirror_tnp = True self.minnlo_unc = 'byHelicityPt' self.skipFromSignal = False @@ -42,7 +41,6 @@ def sample_label(self, sample_group): def configure(self, resumUnc, np_model, transitionUnc = True, propagate_to_fakes=True, - tnp_magnitude=1, tnp_scale=1., mirror_tnp=True, pdf_from_corr=False, @@ -56,13 +54,11 @@ def configure(self, resumUnc, np_model, self.set_minnlo_unc(minnlo_unc) self.transitionUnc = transitionUnc - self.tnp_magnitude = tnp_magnitude self.tnp_scale = tnp_scale self.mirror_tnp = mirror_tnp self.pdf_from_corr = pdf_from_corr self.pdf_operation = pdf_operation self.scale_pdf_unc = scale_pdf_unc - self.minnlo_unc = minnlo_unc self.samples = [] self.skipFromSignal = False @@ -70,7 +66,7 @@ def add_all_theory_unc(self, samples, skipFromSignal=False): self.samples = samples self.skipFromSignal = skipFromSignal self.add_nonpert_unc(model=self.np_model) - self.add_resum_unc(magnitude=self.tnp_magnitude, mirror=self.mirror_tnp, scale=self.tnp_scale) + self.add_resum_unc(scale=self.tnp_scale) self.add_pdf_uncertainty(from_corr=self.pdf_from_corr, operation=self.pdf_operation, scale=self.scale_pdf_unc) try: self.add_quark_mass_vars() @@ -95,22 +91,32 @@ def set_resum_unc_type(self, resumUnc): self.corr_hist = self.card_tool.getHistsForProcAndSyst(signal_samples[0], self.corr_hist_name) if resumUnc.startswith("tnp"): - self.tnp_nuisances = self.card_tool.match_str_axis_entries(self.corr_hist.axes[self.syst_ax], - ["^gamma_.*[+|-]\d+", "^b_.*[+|-]\d+", "^s[+|-]\d+", "^h_.*\d+"]) + self.tnp_nuisance_names = ['b_qg', + 'b_qqDS', + 'b_qqS', + 'b_qqV', + 'b_qqbarV', + 'gamma_cusp', + 'gamma_mu_q', + 'gamma_nu', + 'h_qqV', + 's' + ] + + self.tnp_nuisances = self.card_tool.match_str_axis_entries(self.corr_hist.axes[self.syst_ax], self.tnp_nuisance_names) + if not self.tnp_nuisances: raise ValueError(f"Did not find TNP uncertainties in hist {self.corr_hist_name}") - self.tnp_names = set([x.split("+")[0].split("-")[0] for x in self.tnp_nuisances]) - self.resumUnc = resumUnc - def add_resum_unc(self, magnitude=1, mirror=False, scale=1): + def add_resum_unc(self, scale=1): if not self.resumUnc: logger.warning("No resummation uncertainty will be applied!") return if self.resumUnc.startswith("tnp"): - self.add_resum_tnp_unc(magnitude, mirror, scale) + self.add_resum_tnp_unc(scale) fo_scale = self.resumUnc == "tnp" self.add_transition_fo_scale_uncertainties(transition = self.transitionUnc, scale = fo_scale) @@ -188,6 +194,13 @@ def add_minnlo_scale_uncertainty(self, sample_group, extra_name="", use_hel_hist hscale = self.card_tool.getHistsForProcAndSyst(signal_samples[0], scale_hist) # A bit janky, but refer to the original ptVgen ax since the alt hasn't been added yet orig_binning = hscale.axes[pt_ax.replace("Alt", "")].edges + if not (np.isclose(binning[0], orig_binning[0]) and np.isclose(binning[-1], orig_binning[-1])): + if len(binning) == 2: + binning = np.array((orig_binning[0], orig_binning[-1])) + else: + binning = binning[(binning >= orig_binning[0]-1e-5) & (binning <= orig_binning[-1]+1e-6)] + logger.warning(f"Adjusting requested binning to {binning}") + if not hh.compatibleBins(orig_binning, binning): logger.warning(f"Requested binning {binning} is not compatible with hist binning {orig_binning}. Will not rebin!") binning = orig_binning @@ -196,11 +209,13 @@ def add_minnlo_scale_uncertainty(self, sample_group, extra_name="", use_hel_hist pt_idx = np.argmax(binning >= pt_min) skip_entries.extend([{pt_ax : complex(0, x)} for x in binning[:pt_idx]]) - func = syst_tools.hist_to_variations + func = syst_tools.gen_hist_to_variations if pt_ax == "ptVgenAlt" else syst_tools.hist_to_variations preop_map = {proc : func for proc in expanded_samples} preop_args["gen_axes"] = [pt_ax] preop_args["rebin_axes"] = [pt_ax] preop_args["rebin_edges"] = [binning] + if pt_ax == "ptVgenAlt": + preop_args["gen_obs"] = ["ptVgen"] # Skip MiNNLO unc. if self.resumUnc and not (pt_binned or helicity): @@ -262,13 +277,17 @@ def add_scetlib_dyturbo_scale_uncertainty(self, extra_name="", transition = True def preop_func(h, *args, **kwargs): hsel = h[{"vars" : ["pdf0"] + sel_vars}] - return syst_tools.hist_to_variations(hsel, *args, **kwargs) + func = syst_tools.gen_hist_to_variations if pt_ax == "ptVgenAlt" else syst_tools.hist_to_variations + return func(hsel, *args, **kwargs) preop_args = {} preop_args["gen_axes"] = [pt_ax] preop_args["rebin_axes"] = [pt_ax] preop_args["rebin_edges"] = [binning] + if pt_ax == "ptVgenAlt": + preop_args["gen_obs"] = ["ptVgen"] + self.card_tool.addSystematic(name=self.scale_hist_name, processes=[sample_group], group="resumTransitionFOScale", @@ -288,36 +307,21 @@ def preop_func(h, *args, **kwargs): def set_propagate_to_fakes(self, to_fakes): self.propagate_to_fakes = to_fakes - def add_resum_tnp_unc(self, magnitude, mirror, scale=1): + def add_resum_tnp_unc(self, magnitude, scale=1): syst_ax = self.corr_hist.axes[self.syst_ax] - np_mag = lambda x: x.split("+")[-1].split("-")[-1] - tnp_magnitudes = set([np_mag(x) for x in self.tnp_nuisances]) tnp_has_mirror = ("+" in self.tnp_nuisances[0] and self.tnp_nuisances[0].replace("+", "-") in syst_ax) or \ - ("-" in self.tnp_nuisances[0] and self.tnp_nuisances[0].replace("-", "+") in syst_ax) + ("-" in self.tnp_nuisances[0] and self.tnp_nuisances[0].replace("-", "") in syst_ax) + if not tnp_has_mirror and not self.mirror_tnp: + logger.warning("Up or down variation missing in TNP histogram. Will use mirroring") + self.mirror_tnp = True - if not any(np.isclose(float(x), magnitude) for x in tnp_magnitudes): - raise ValueError(f"TNP magnitude variation {magnitude} is not present in the histogram {self.corr_hist_name}. Options are {tnp_magnitudes}") - - qqV_mag = magnitude/2. - if qqV_mag == 2.5: - qqV_mag = 2.0 - - selected_tnp_nuisances = [x for x in self.tnp_nuisances if np.isclose(magnitude, float(np_mag(x))) or ("h_qqV" in x and np.isclose(qqV_mag, float(np_mag(x))))] central_var = syst_ax[0] - if mirror: - name_replace = [] - if tnp_has_mirror: - selected_tnp_nuisances = [x for x in selected_tnp_nuisances if "+" in x] - else: - if not tnp_has_mirror: - raise ValueError(f"Uncertainty hist {self.corr_hist_name} does not have double sided TNP nuisances. " \ - "mirror=False is therefore not defined") - name_replace = [(f"+{m}", "Up") for m in tnp_magnitudes]+[(f"-{m}", "Up") for m in tnp_magnitudes] + tnp_magnitudes = ["2.5", "0.5", "1."] + name_replace = [(f"-{x}", "Down") for x in tnp_magnitudes]+[(x, "Up") for x in tnp_magnitudes] - logger.debug(f"TNP nuisances in correction hist: {self.tnp_nuisances}") - logger.debug(f"Selected TNP nuisances: {selected_tnp_nuisances}") + logger.debug(f"Selected TNP nuisances: {self.tnp_nuisance_names}") self.card_tool.addSystematic(name=self.corr_hist_name, processes=['wtau_samples', 'single_v_nonsig_samples'] if self.skipFromSignal else ['single_v_samples'], @@ -326,8 +330,8 @@ def add_resum_tnp_unc(self, magnitude, mirror, scale=1): systAxes=["vars"], passToFakes=self.propagate_to_fakes, systNameReplace=name_replace, - preOp=lambda h: h[{self.syst_ax : [central_var, *selected_tnp_nuisances]}], - mirror=mirror, + preOp=lambda h: h[{self.syst_ax : [central_var, *self.tnp_nuisances]}], + mirror=self.mirror_tnp, scale=scale, skipEntries=[{self.syst_ax : central_var},], rename=f"resumTNP", diff --git a/wremnants/plot_tools.py b/wremnants/plot_tools.py index 01e4c4b51..e160574d3 100644 --- a/wremnants/plot_tools.py +++ b/wremnants/plot_tools.py @@ -353,14 +353,13 @@ def makePlotWithRatioToRef( hists, labels, colors, linestyles=[], xlabel="", ylabel="Events/bin", rlabel="x/nominal", rrange=[0.9, 1.1], ylim=None, xlim=None, nlegcols=2, binwnorm=None, alpha=1., - baseline=True, data=False, autorrange=None, grid = False, extra_text=None, extra_text_loc=(0.8, 0.7), + baseline=True, dataIdx=None, autorrange=None, grid = False, extra_text=None, extra_text_loc=(0.8, 0.7), yerr=False, legtext_size=20, plot_title=None, x_ticks_ndp = None, bin_density = 300, yscale=None, - logy=False, logx=False, fill_between=False, title_padding = 0, cms_label = None + logy=False, logx=False, fill_between=0, title_padding = 0, cms_label = None, ): if len(hists) != len(labels) or len(hists) != len(colors): raise ValueError(f"Number of hists ({len(hists)}), colors ({len(colors)}), and labels ({len(labels)}) must agree!") - # nominal is always at first, data is always at last, if included - ratio_hists = [hh.divideHists(h, hists[0], cutoff=1e-6, flow=False, by_ax_name=False) for h in hists[not baseline:]] + ratio_hists = [hh.divideHists(h, hists[0], cutoff=1e-6, flow=False, rel_unc=True, by_ax_name=False) for h in hists[not baseline:]] fig, ax1, ax2 = figureWithRatio( hists[0], xlabel, ylabel, ylim, rlabel, rrange, xlim=xlim, grid_on_ratio_plot = grid, plot_title = plot_title, title_padding=title_padding, @@ -368,14 +367,15 @@ def makePlotWithRatioToRef( ) linestyles = linestyles+['solid']*(len(hists)-len(linestyles)) + + exclude_data = lambda x,idx=dataIdx: [j for i,j in enumerate(x) if i != dataIdx] - count = len(hists)-data hep.histplot( - hists[:count], + exclude_data(hists), histtype="step", - color=colors[:count], - label=labels[:count], - linestyle=linestyles[:count], + color=exclude_data(colors), + label=exclude_data(labels), + linestyle=exclude_data(linestyles), stack=False, ax=ax1, yerr=yerr, @@ -385,34 +385,31 @@ def makePlotWithRatioToRef( ) if len(hists) > 1: - ratio_hists = [hh.divideHists(h, hists[0], cutoff=0.00001, flow=False, by_ax_name=False) for h in hists[not baseline:]] - if fill_between: - for up,down,color in zip(hists[1::2], hists[2::2], colors[1::2]): - upr = hh.divideHists(up, hists[0], 1e-6, flow=False, by_ax_name=False) - downr = hh.divideHists(down, hists[0], 1e-6,flow=False, by_ax_name=False) + ratio_hists = [hh.divideHists(h, hists[0], cutoff=0.00001, flow=False, rel_unc=True, by_ax_name=False) for h in hists] + if fill_between != 0: + for upr,downr,color in zip(ratio_hists[-fill_between::2], ratio_hists[-fill_between+1::2], colors[-fill_between::2]): ax2.fill_between(upr.axes[0].edges, np.append(upr.values(), upr.values()[-1]), np.append(downr.values(), downr.values()[-1]), step='post', color=color, alpha=0.5) - count = len(ratio_hists) - data if not fill_between else 1 hep.histplot( - ratio_hists[(not baseline):count], + exclude_data(ratio_hists)[not baseline:], histtype="step", - color=colors[(not baseline):count], - linestyle=linestyles[(not baseline):count], + color=exclude_data(colors)[not baseline:], + linestyle=exclude_data(linestyles)[not baseline:], yerr=yerr, stack=False, ax=ax2, alpha=alpha, flow='none', ) - if data: + if dataIdx is not None: hep.histplot( - hists[-1], + hists[dataIdx], histtype="errorbar", - color=colors[-1], - label=labels[-1], + color=colors[dataIdx], + label=labels[dataIdx], stack=False, ax=ax1, binwnorm=binwnorm, @@ -420,9 +417,9 @@ def makePlotWithRatioToRef( flow='none', ) hep.histplot( - hh.divideHists(data, hists[0], cutoff=1.e-8, flow=False, by_ax_name=False), + hh.divideHists(hists[dataIdx], hists[0], cutoff=1.e-8, flow=False, by_ax_name=False, rel_unc=True), histtype="errorbar", - color=colors[-1], + color=colors[dataIdx], xerr=False, yerr=True, stack=False, diff --git a/wremnants/syst_tools.py b/wremnants/syst_tools.py index e7eb12a84..9840872a4 100644 --- a/wremnants/syst_tools.py +++ b/wremnants/syst_tools.py @@ -66,6 +66,10 @@ def projAx(hname): return hname.split("-") resum_tnps = ['pdf0', 'gamma_cusp+1', 'gamma_mu_q+1', 'gamma_nu+1', 'h_qqV-0.5', 's+1', 'b_qqV+1', 'b_qqbarV+1', 'b_qqS+1', 'b_qqDS+1', 'b_qg+1'] + resum_tnps2p1_up = ['pdf0', 'gamma_cusp1.', 'gamma_mu_q1.', 'gamma_nu1.', 's1.', 'b_qqV0.5', 'b_qqV0.5', 'b_qqbarV0.5', 'b_qqS0.5', 'b_qqDS0.5', 'b_qg0.5'] + resum_tnps2p1_down = ['pdf0', 'gamma_cusp-1.', 'gamma_mu_q-1.', 'gamma_nu-1.', 's-1.', 'b_qqV-2.5', 'b_qqV-2.5', 'b_qqbarV-2.5', 'b_qqS-2.5', 'b_qqDS-2.5', 'b_qg-2.5'] + resum_tnps3p0_up = ['pdf0', 'gamma_cusp1.', 'gamma_mu_q1.', 'gamma_nu1.', 's1.', 'b_qqV0.5', 'b_qqV0.5', 'b_qqbarV0.5', 'b_qqS0.5', 'b_qqDS0.5', 'b_qg0.5'] + resum_tnps3p0_down = ['pdf0', 'gamma_cusp-1.', 'gamma_mu_q-1.', 'gamma_nu-1.', 's-1.', 'b_qqV-0.5', 'b_qqV-0.5', 'b_qqbarV-0.5', 'b_qqS-0.5', 'b_qqDS-0.5', 'b_qg-0.5'] transforms.update({ "resumFOScaleUp" : { @@ -86,11 +90,17 @@ def projAx(hname): ["transition_points0.2_0.65_1.1", "transition_points0.4_0.55_0.7", "transition_points0.2_0.45_0.7", "transition_points0.4_0.75_1.1", ], no_flow=["ptVgen"], do_min=True)}, - "resumTNPUp" : { - "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps}], "vars")[0] + "resumTNP2p1Up" : { + "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps3p0_up}], "vars")[0] }, - "resumTNPDown" : { - "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps}], "vars")[1] + "resumTNP3p0Down" : { + "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps3p0_down}], "vars")[1] + }, + "resumTNP3p0Up" : { + "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps2p1_up}], "vars")[0] + }, + "resumTNP2p1Down" : { + "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps2p1_down}], "vars")[1] }, "resumTNPx5Up" : { "action" : lambda h: h if "vars" not in h.axes.name else hh.rssHists(h[{"vars" : resum_tnps}], "vars", scale=5)[0] @@ -303,6 +313,12 @@ def make_fakerate_variation(href, fakerate_axes, fakerate_axes_syst, variation_f # 5) add back to the nominal histogram and broadcase the nominal histogram return hh.addHists(href, hsyst) +def gen_hist_to_variations(hist_in, gen_obs, gen_axes=["ptVgen", "chargeVgen", "helicity"], sum_axes=[], rebin_axes=[], rebin_edges=[]): + for obs in gen_obs: + hist_in = hh.expand_hist_by_duplicate_axis(hist_in, obs, obs+"Alt", swap_axes=True) + + return hist_to_variations(hist_in, gen_axes, sum_axes, rebin_axes, rebin_edges) + def hist_to_variations(hist_in, gen_axes = [], sum_axes = [], rebin_axes=[], rebin_edges=[]): if hist_in.name is None: diff --git a/wremnants/theory_corrections.py b/wremnants/theory_corrections.py index e03dbe7c0..9e3ccdb05 100644 --- a/wremnants/theory_corrections.py +++ b/wremnants/theory_corrections.py @@ -115,10 +115,7 @@ def compute_envelope(h, name, entries, axis_name="vars", slice_axis = None, slic def postprocess_corr_hist(corrh): # extend variations with some envelopes and special kinematic slices - if "vars" not in corrh.axes.name: - return corrh - - if type(corrh.axes["vars"]) != hist.axis.StrCategory: + if "vars" not in corrh.axes.name or type(corrh.axes["vars"]) != hist.axis.StrCategory: return corrh additional_var_hists = {} @@ -131,7 +128,7 @@ def postprocess_corr_hist(corrh): resum_scale_vars_exclusive = [var for var in corrh.axes["vars"] if any(resum_scale in var for resum_scale in resum_scales)] resum_scale_vars = ["pdf0"] + resum_scale_vars_exclusive - if len(resum_scale_vars) == 1: + if len(renorm_fact_scale_vars) == 1: return corrh transition_vars_exclusive = ["transition_points0.2_0.35_1.0", "transition_points0.2_0.75_1.0"] @@ -158,7 +155,7 @@ def postprocess_corr_hist(corrh): additional_var_hists.update(compute_envelope(corrh, "renorm_fact_resum_scale_envelope", renorm_fact_resum_scale_vars)) if all(var in corrh.axes["vars"] for var in renorm_fact_resum_transition_scale_vars): additional_var_hists.update(compute_envelope(corrh, "renorm_fact_resum_transition_scale_envelope", renorm_fact_resum_transition_scale_vars)) - if all(var in corrh.axes["vars"] for var in resum_scale_vars): + if len(resum_scale_vars) > 1 and all(var in corrh.axes["vars"] for var in resum_scale_vars): additional_var_hists.update(compute_envelope(corrh, "resum_scale_envelope", resum_scale_vars)) @@ -184,7 +181,7 @@ def postprocess_corr_hist(corrh): def get_corr_name(generator): # Hack for now label = generator.replace("1D", "") - if "dataPtll" in generator: + if "dataPtll" in generator or "dataRecoPtll" in generator: return "MC_data_ratio" return f"{label}_minnlo_ratio" if "Helicity" not in generator else f"{label.replace('Helicity', '')}_minnlo_coeffs" diff --git a/wremnants/theory_tools.py b/wremnants/theory_tools.py index ebd63a626..06683c613 100644 --- a/wremnants/theory_tools.py +++ b/wremnants/theory_tools.py @@ -320,7 +320,7 @@ def define_postfsr_vars(df, mode=None): if mode is not None: # defition of more complex postfsr object # use fiducial gen met, see: https://twiki.cern.ch/twiki/bin/viewauth/CMS/ParticleLevelProducer - if mode in ["wlike", "dilepton"]: + if "mz" in mode: # find the leading charged lepton and antilepton idx df = df.Define("postfsrLep", "postfsrLeptons && (GenPart_pdgId==11 || GenPart_pdgId==13)") df = df.Define("postfsrAntiLep", "postfsrLeptons && (GenPart_pdgId==-11 || GenPart_pdgId==-13)") @@ -354,7 +354,7 @@ def define_postfsr_vars(df, mode=None): df = df.Define("postfsrLep_absEta", "static_cast(std::fabs(postfsrLep_eta))") - if mode in ["wmass", "wlike"]: + if "singlelep" in mode: if mode == "wlike": # for wlike selection df = df.Define("postfsrMET_wlike", "wrem::get_met_wlike(postfsrOtherLep_pt, postfsrOtherLep_phi, MET_fiducialGenPt, MET_fiducialGenPhi)") @@ -368,7 +368,7 @@ def define_postfsr_vars(df, mode=None): df = df.Define("postfsrDeltaPhiMuonMet", "std::fabs(wrem::deltaPhi(postfsrLep_phi, postfsrMET_phi))") # definition of boson kinematics - if mode in ["dilepton", "wlike"]: + if "mz" in mode: # four vectors df = df.Define("postfsrLep_mom4", "ROOT::Math::PtEtaPhiMVector(postfsrLep_pt, postfsrLep_eta, postfsrLep_phi, postfsrLep_mass)") df = df.Define("postfsrAntiLep_mom4", "ROOT::Math::PtEtaPhiMVector(postfsrOtherLep_pt, postfsrOtherLep_eta, postfsrOtherLep_phi, postfsrOtherLep_mass)") @@ -379,7 +379,7 @@ def define_postfsr_vars(df, mode=None): df = df.Define('postfsrabsYV', 'std::fabs(postfsrYV)') df = df.Define('postfsrPTV', 'postfsrGenV_mom4.pt()') - if mode == "wmass": + if "mw" in mode: df = df.Define("postfsrPTV", "wrem::pt_2(postfsrLep_pt, postfsrLep_phi, postfsrMET_pt, postfsrMET_phi)") return df diff --git a/wremnants/unfolding_tools.py b/wremnants/unfolding_tools.py index e4c08d674..2c99cddba 100644 --- a/wremnants/unfolding_tools.py +++ b/wremnants/unfolding_tools.py @@ -32,9 +32,8 @@ def define_gen_level(df, gen_level, dataset_name, mode="wmass"): if gen_level not in gen_levels: raise ValueError(f"Unknown gen level '{gen_level}'! Supported gen level definitions are '{gen_levels}'.") - modes = ["wmass", "wlike", "dilepton"] - if mode not in modes: - raise ValueError(f"Unknown mode '{mode}'! Supported modes are '{modes}'.") + logger.info(f"Using {gen_level} leptons") + mz = "mz" in mode if gen_level == "preFSR": df = theory_tools.define_prefsr_vars(df, mode=mode) @@ -47,11 +46,11 @@ def define_gen_level(df, gen_level, dataset_name, mode="wmass"): if mode in ["wmass", "wlike"]: df = df.Alias("mTVGen", "mTVgen") - if mode == "wmass": + if "mw" in mode: df = df.Define("ptGen", "chargeVgen < 0 ? genl.pt() : genlanti.pt()") df = df.Define("absEtaGen", "chargeVgen < 0 ? std::fabs(genl.eta()) : std::fabs(genlanti.eta())") - if mode in ["wlike", "dilepton"]: + if mz: df = df.Define("ptGen", "event % 2 == 0 ? genl.pt() : genlanti.pt()") df = df.Define("absEtaGen", "event % 2 == 0 ? std::fabs(genl.eta()) : std::fabs(genlanti.eta())") df = df.Define("ptOtherGen", "event % 2 == 0 ? genlanti.pt() : genl.pt()") @@ -66,7 +65,7 @@ def define_gen_level(df, gen_level, dataset_name, mode="wmass"): if mode in ["wmass", "wlike"]: df = df.Alias("mTVGen", "postfsrMT") - if mode in ["wlike", "dilepton"]: + if mz: df = df.Alias("ptOtherGen", "postfsrOtherLep_pt") df = df.Alias("absEtaOtherGen", f"postfsrOtherLep_absEta") @@ -75,31 +74,44 @@ def define_gen_level(df, gen_level, dataset_name, mode="wmass"): df = df.Alias("ptVGen", "postfsrPTV") - if mode == "wlike": + if "wlike" in mode: df = df.Define("qGen", "event % 2 == 0 ? -1 : 1") return df -def select_fiducial_space(df, select=True, accept=True, mode="wmass", pt_min=None, pt_max=None, mass_min=60, mass_max=120, mtw_min=0, selections=[]): +def get_fiducial_args(mode, pt_min=28, pt_max=60, abseta_max=2.4): + fidargs = {} + if "inclusive" in mode or mode == "mz_masswindow": + fidargs = {"abseta_max" : 100.} + if mode == "mz_inclusive": + fidargs.update({"mass_min" : 60, "mass_max" : 120}) + return fidargs + + fidargs.update({"mtw_min" : 40 if "mw" in mode or "singlelep" in mode else 0, + "pt_min" : pt_min, "pt_max" : pt_max, "abseta_max" : abseta_max}) + + return fidargs + +def select_fiducial_space(df, select=True, accept=True, mode="mw", pt_min=0, pt_max=1300, abseta_max=2.4, mass_min=60, mass_max=120, mtw_min=0, selections=[]): # Define a fiducial phase space and if select=True, either select events inside/outside # accept = True: select events in fiducial phase space # accept = False: reject events in fiducial pahse space - if mode == "wmass": - selection = "(absEtaGen < 2.4)" - elif mode == "wlike": + if "mw" in mode: + selection = f"(absEtaGen < {abseta_max})" + elif "singlelep" in mode: selection = f""" - (absEtaGen < 2.4) && (absEtaOtherGen < 2.4) + (absEtaGen < {abseta_max}) && (absEtaOtherGen < {abseta_max}) && (ptOtherGen > {pt_min}) && (ptOtherGen < {pt_max}) && (massVGen > {mass_min}) && (massVGen < {mass_max}) """ - elif mode == "dilepton": + elif "mz" in mode: selection = f""" - (absEtaGen < 2.4) && (absEtaOtherGen < 2.4) - && (ptGen > {pt_min}) && (ptOtherGen > {pt_min}) - && (ptGen < {pt_max}) && (ptOtherGen < {pt_max}) + (absEtaGen < {abseta_max}) && (absEtaOtherGen < {abseta_max}) + && (ptGen > {pt_min}) && (ptOtherGen > {pt_min}) + && (ptGen < {pt_max}) && (ptOtherGen < {pt_max}) && (massVGen > {mass_min}) && (massVGen < {mass_max}) - """ + """ else: raise NotImplementedError(f"No fiducial phase space definiton found for mode '{mode}'!") @@ -110,6 +122,8 @@ def select_fiducial_space(df, select=True, accept=True, mode="wmass", pt_min=Non logger.debug(f"Add selection {sel} for fiducial phase space") selection += f" && ({sel})" + logger.info(f"Applying fiducial selection {selection}") + df = df.Define("acceptance", selection) if select and accept: