240418 dev

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()

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")

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())}")

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:

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] 

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
 

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)

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")

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
+    return result, meta