Fitting new lattice-constrained NP model, PDFs from correction files through helicities

wremnants/combine_theory_helper.py

@@ -912,14 +912,18 @@ def add_pdf_uncertainty(self, operation=None, scale=-1.0):
             if scale != -1.0
             else theory_tools.pdf_inflation_factor(pdfInfo, self.args.noi)
         )
+        pdf_hist = pdfName
+        pdf_corr_hist = (
+            f"scetlib_dyturbo{pdf.upper().replace('AN3LO', 'an3lo')}VarsCorr"
+            if self.corr_hist_name == "scetlib_dyturboCorr"
+            else self.corr_hist_name.replace(
+                "Corr", "_CT18ZVarsCorr"
+            )  # TODO how do we get around this?
+        )
         if self.from_hels:
-            pdf_hist = f"{pdfName}UncertByHelicity"
-            pdf_corr_hist = f"{pdfName}UncertByHelicity"
-        else:
-            pdf_hist = pdfName
-            pdf_corr_hist = (
-                f"scetlib_dyturbo{pdf.upper().replace('AN3LO', 'an3lo')}VarsCorr"
-            )
+            pdf_hist += "UncertByHelicity"
+            pdf_corr_hist += "UncertByHelicity"
+        symmetrize = "quadratic"
 
         if self.pdf_from_corr:
             theory_unc = self.datagroups.args_from_metadata("theoryCorr")
@@ -998,7 +1002,13 @@ def add_pdf_alphas_variation(self, noi=False):
             asname = pdf_corr_hist.replace("Vars", "_pdfas")
             # alphaS from correction histograms only available for these sets,
             # so fall back to CT18Z for other sets
-            if not ("MSHT20" in asname or "CT18Z" in asname or "MSHT20an3lo" in asname):
+            if not (
+                "MSHT20" in asname
+                or "CT18Z" in asname
+                or "MSHT20an3lo" in asname
+                or "Lattice" in asname
+                or "nnlojet" in asname
+            ):  # TODO should fix the correction file name for Lattice?
                 asname = "scetlib_dyturboCT18Z_pdfasCorr"
                 as_range = theory_tools.pdfMap["ct18z"]["alphasRange"]
             if asname.replace("Corr", "") not in self.datagroups.args_from_metadata(

wremnants/syst_tools.py

@@ -1747,7 +1747,7 @@ def add_pdfUncertByHelicity_hist(
             ],
         )
     safeTensorName = f"{tensorName}_clamped"
-    renorm = theory_tools.pdfMap[pdf].get("renorm", False)
+    renorm = theory_tools.pdfMap.get(pdf, {}).get("renorm", False)
     if renorm:
         central_event_weight = "nominal_weight"
     else:
@@ -2597,11 +2597,27 @@ def add_theory_hists(
                     cols,
                     **info,
                 )
+        if theory_helpers.get("pdf_from_corr") is not None:
+            pdf_from_corr_helpers = theory_helpers.get("pdf_from_corr")
+            for pdf in pdf_from_corr_helpers:
+                logger.debug(
+                    f"Make PDF (from correction file) uncertainty by helicity histograms for {dataset_name} and PDF from correction {pdf}"
+                )
+                add_pdfUncertByHelicity_hist(
+                    results,
+                    df,
+                    pdf_from_corr_helpers[pdf],
+                    pdf,
+                    pdf,
+                    axes,
+                    cols,
+                    **info,
+                )
         if theory_helpers.get("alphaS") is not None:
-            logger.debug(
-                f"Make AlphaS uncertainty by helicity histograms for {dataset_name}"
-            )
             for k, v in theory_helpers["alphaS"].items():
+                logger.debug(
+                    f"Make alphaS uncertainty by helicity histogram for {dataset_name} and alphaS from correction {k}"
+                )
                 add_pdfAlphaSByHelicity_hist(results, df, v, axes, cols, name=k, **info)
 
         add_breit_wigner_mass_weights_hist(

wremnants/theory_corrections.py

@@ -512,7 +512,9 @@ def make_corr_by_helicity(
 
 
 def make_theory_helpers(
-    args, procs=["Z", "W"], corrs=["qcdScale", "pdf", "alphaS", "pdf_central"]
+    args,
+    procs=["Z", "W"],
+    corrs=["qcdScale", "pdf", "pdf_from_corr", "alphaS", "pdf_central"],
 ):
 
     theory_helpers_procs = {p: {} for p in procs}
@@ -545,6 +547,14 @@ def make_theory_helpers(
                     pdfs=args.pdfs,
                 )
             )
+        if "pdf_from_corr" in corrs:
+            pdf_from_corrs = [x + "Corr" for x in args.theoryCorr if "Vars" in x]
+            theory_helpers_procs[proc]["pdf_from_corr"] = (
+                make_pdfs_from_corrs_uncertainties_helper_by_helicity(
+                    proc=proc,
+                    pdfs_from_corrs=pdf_from_corrs,
+                )
+            )
         if "alphaS" in corrs:
             as_vars = [x + "Corr" for x in args.theoryCorr if "pdfas" in x]
             theory_helpers_procs[proc]["alphaS"] = (
@@ -725,6 +735,32 @@ def make_pdfs_uncertainties_helper_by_helicity(
     return pdf_helpers
 
 
+def make_pdfs_from_corrs_uncertainties_helper_by_helicity(
+    proc,
+    pdfs_from_corrs,
+    return_tensor=True,
+):
+    pdf_file_template = (
+        common.data_dir
+        + "/TheoryCorrectionsByHelicity/PDFsFromCorrs/w_z_gen_dists_{pdf}_maxFiles_m1_skimmed.hdf5"
+    )
+    pdf_helpers = {}
+    for pdf in pdfs_from_corrs:
+        logger.debug(f"Making PDF uncertainty helper by helicity for theory corr {pdf}")
+        pdf_helper = make_uncertainty_helper_by_helicity(
+            proc=proc,
+            nom=pdf,
+            den="pdf_uncorr",
+            central_weights=False,
+            var_ax_name="vars",
+            filename=pdf_file_template.format(pdf=pdf),
+            return_tensor=return_tensor,
+        )
+        if pdf_helper is not None:
+            pdf_helpers[pdf] = pdf_helper
+    return pdf_helpers
+
+
 def make_alphaS_uncertainties_helper_by_helicity(
     proc,
     as_vars,
@@ -739,11 +775,12 @@ def make_alphaS_uncertainties_helper_by_helicity(
         logger.debug(
             f"Making alphaS uncertainty helper by helicity for theory corr {as_var}"
         )
+        fname = alphas_file_template.format(as_var=as_var)
         as_helper = make_uncertainty_helper_by_helicity(
             proc=proc,
             nom=as_var,
             den="theory_uncorr",
-            filename=alphas_file_template.format(as_var=as_var),
+            filename=fname,
             var_ax_name="vars",
             return_tensor=return_tensor,
         )
@@ -779,12 +816,17 @@ def _collect_hist(hist_name, filename):
         hist_key = f"nominal_gen_{hist_name}"
         hists = []
         for process in proc_map.get(proc, ()):
+            if not os.path.exists(filename):
+                logger.warning(
+                    f"File {filename} does not exist. Not creating histogram of variations by helicities for process {proc} and variation {nom}."
+                )
+                return None
             with h5py.File(filename, "r") as h5file:
                 results = input_tools.load_results_h5py(h5file)
-                outputs = results[process]["output"]
+                outputs = results.get(process, {}).get("output", {})
                 if hist_key not in outputs:
                     logger.warning(
-                        f"Did not find {hist_key} in {filename}. Not creating histogram of PDF variations by helicities for this set."
+                        f"Did not find {hist_key} in {filename}. Not creating histogram of variations by helicities for process {proc} and variation {nom}."
                     )
                     return None
                 hists.append(outputs[hist_key].get())

wremnants/theory_tools.py

@@ -56,7 +56,7 @@
         "alphas": ["LHEPdfWeight[0]", "LHEPdfWeight[101]", "LHEPdfWeight[102]"],
         "alphasRange": "002",
         "inflation_factor_wmass": 3.0,
-        "inflation_factor_alphaS": 3.5,
+        "inflation_factor_alphaS": 3.0,
     },
     "ct18": {
         "name": "pdfCT18",
@@ -71,7 +71,7 @@
         "alphasRange": "002",
         "scale": 1 / 1.645,  # Convert from 90% CL to 68%
         "inflation_factor_wmass": 1.0,
-        "inflation_factor_alphaS": 1.0,
+        "inflation_factor_alphaS": 1.2,
     },
     "nnpdf30": {
         "name": "pdfNNPDF30",
@@ -113,7 +113,7 @@
         ],
         "alphasRange": "001",
         "inflation_factor_wmass": 1.0,
-        "inflation_factor_alphaS": 1.2,
+        "inflation_factor_alphaS": 1.5,
     },
     "msht20": {
         "name": "pdfMSHT20",
@@ -127,7 +127,7 @@
         ],
         "alphasRange": "002",
         "inflation_factor_wmass": 1.5,
-        "inflation_factor_alphaS": 1.7,
+        "inflation_factor_alphaS": 2.0,
     },
     "msht20mcrange": {
         "name": "pdfMSHT20mcrange",
@@ -211,7 +211,7 @@
         ],  # alphas 116-120
         "alphasRange": "002",
         "inflation_factor_wmass": 4.0,
-        "inflation_factor_alphaS": 3.0,
+        "inflation_factor_alphaS": 3.5,
     },
     "herapdf20ext": {
         "name": "pdfHERAPDF20ext",
@@ -285,6 +285,7 @@ def define_scale_tensor(df):
     "scetlib_dyturboMSHT20_pdfas": expand_pdf_entries("msht20", alphas=True),
     "scetlib_dyturboMSHT20Vars": expand_pdf_entries("msht20"),
     "scetlib_dyturboCT18ZVars": expand_pdf_entries("ct18z"),
+    "scetlib_dyturboN3p0LL_LatticeNP_CT18ZVars": expand_pdf_entries("ct18z"),
     "scetlib_dyturboCT18Z_pdfas": expand_pdf_entries("ct18z", alphas=True, renorm=True),
     "scetlib_dyturboN3p1LL_pdfas": expand_pdf_entries(
         "ct18z", alphas=True, renorm=True