Merge pull request #33 from LPC-HH/fix-pu-weights

Updates to PU and v12

.copier-answers.yml

@@ -9,4 +9,4 @@ org: LPC-HH
 project_name: HH4b
 project_short_description: HH4b
 url: https://github.com/LPC-HH/HH4b
-vcs: true
+vcs: False

.pre-commit-config.yaml

@@ -41,7 +41,7 @@ repos:
     rev: "v3.0.3"
     hooks:
       - id: prettier
-        types_or: [yaml, markdown, html, css, scss, javascript, json]
+        types_or: [markdown, html, css, scss, javascript, json]
         args: [--prose-wrap=always]
 
   - repo: https://github.com/astral-sh/ruff-pre-commit

data/rucio_utils.py

@@ -142,7 +142,7 @@ def get_dataset_files(
                     found = True
 
             if not found:
-                raise Exception(f"No SITE available in the whitelist for file {filedata['name']}")
+                print(f"No SITE available in the whitelist for file {filedata['name']}")
         else:
             possible_sites = list(rses.keys())
             if blacklist_sites:
@@ -151,7 +151,7 @@ def get_dataset_files(
                 )
 
             if len(possible_sites) == 0:
-                raise Exception(f"No SITE available for file {filedata['name']}")
+                print(f"No SITE available for file {filedata['name']}")
 
             # now check for regex
             for site in possible_sites:
@@ -187,7 +187,7 @@ def get_dataset_files(
                     found = True
 
         if not found:
-            raise Exception(f"No SITE available for file {filedata['name']}")
+            print(f"No SITE available for file {filedata['name']}")
         else:
             if output == "all":
                 outfiles.append(outfile)

pyproject.toml

@@ -1,5 +1,5 @@
 [build-system]
-requires = ["hatchling", "hatch-vcs"]
+requires = ["hatchling"]
 build-backend = "hatchling.build"
 
 
@@ -56,8 +56,7 @@ Changelog = "https://github.com/LPC-HH/HH4b/releases"
 
 
 [tool.hatch]
-version.source = "vcs"
-build.hooks.vcs.version-file = "src/HH4b/_version.py"
+version.path = "src/HH4b/__init__.py"
 envs.default.dependencies = [
   "pytest",
   "pytest-cov",
@@ -145,6 +144,8 @@ ignore = [
   "F811",   # redefinition of unused variable
   "PGH001", # eval
   "RET503",
+  "RET504", # unnecessary assignment before return
+  "RET505", # unnecessary else after return statement
 ]
 isort.required-imports = ["from __future__ import annotations"]
 # Uncomment if using a _compat.typing backport

src/HH4b/__init__.py

@@ -7,6 +7,6 @@
 
 from __future__ import annotations
 
-from ._version import version as __version__
+__version__ = "0.1.0"
 
 __all__ = ["__version__"]

src/HH4b/corrections/README.md

@@ -10,5 +10,12 @@ https://twiki.cern.ch/twiki/bin/view/CMS/LumiRecommendationsRun3
 
 ## Pileup reweighting
 
-Following https://gitlab.cern.ch/HiggsDNA-project/HiggsDNA/-/merge_requests/47
-Note that here, 80mb was used as the inelastic proton scattering cross-section.
+No official recommendation:
+https://cms-talk.web.cern.ch/t/pu-reweighing-in-2022mc/28786/9
+
+- Following https://gitlab.cern.ch/HiggsDNA-project/HiggsDNA/-/merge_requests/47
+- Updated with
+  https://gitlab.cern.ch/HiggsDNA-project/HiggsDNA/-/merge_requests/74
+
+Using the "CMS value" of 69.2 mb instead of 80 mb, which better matches the
+data.

src/HH4b/corrections/makePUReWeightJSON.py

@@ -17,7 +17,7 @@
     # ========#
     # https://github.com/cms-sw/cmssw/pull/34460
     # https://github.com/cms-sw/cmssw/blob/master/SimGeneral/MixingModule/python/Run3_2022_LHC_Simulation_10h_2h_cfi.py#L10-L54
-    "2022Prompt_25ns": (
+    "2022_LHC_Simulation_10h_2h": (
         np.linspace(0.0, 100.0, 101),
         [
             7.075550618391933e-8,
@@ -212,11 +212,13 @@ def main():
     parser.add_argument(
         "--up",
         type=str,
+        default=None,
         help="File with the data (true) pileup distribution histogram assuming the nominal+1sigma minimum bias cross-section value",
     )
     parser.add_argument(
         "--down",
         type=str,
+        default=None,
         help="File with the data (true) pileup distribution histogram assuming the nominal-1sigma minimum bias cross-section value",
     )
     parser.add_argument("--rebin", type=int, help="Factor to rebin the data histograms by")
@@ -434,7 +436,12 @@ def main():
         rBinCenters = 0.5 * (ratioBins[:-1] + ratioBins[1:])
         ax.hist(dBinCenters, bins=mcPUBins, weights=mcPUVals, histtype="step", label="MC")
         ax.hist(
-            nBinCenters, bins=nomBins, weights=nomCont, histtype="step", label="Nominal", color="k"
+            nBinCenters,
+            bins=nomBins,
+            weights=nomCont,
+            histtype="step",
+            label="Nominal (data)",
+            color="k",
         )
         rax.hist(rBinCenters, bins=ratioBins, weights=nomRatio, histtype="step", color="k")
         if upCont is not None:
@@ -453,12 +460,16 @@ def main():
             rax.hist(rBinCenters, bins=ratioBins, weights=downRatio, histtype="step", color="b")
         rax.axhline(1.0)
         ax.legend()
-        rax.set_ylim(0.02, 2.0)
+        rax.grid()
+        ax.grid()
+        rax.set_ylim(0.02, 20)
         rax.set_xlim(ratioBins[0], ratioBins[-1])
         if args.mcfiles:
             rax.set_xlabel(args.mcreweightvar)
         elif args.mcprofile:
             ax.set_title(args.mcprofile)
+        ax.set_ylabel("PU profile")
+        rax.set_ylabel("PU weight")
         if args.output.endswith(".json"):
             plt.savefig(args.output.replace(".json", ".png"))
 

src/HH4b/corrections/pu_correction_per_file.py

@@ -0,0 +1,104 @@
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+import hist
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mticker
+import mplhep as hep
+import numpy as np
+import uproot
+from coffea.nanoevents import NanoAODSchema, NanoEventsFactory
+
+hep.style.use(["CMS", "firamath"])
+formatter = mticker.ScalarFormatter(useMathText=True)
+formatter.set_powerlimits((-3, 3))
+plt.rcParams.update({"font.size": 12})
+plt.rcParams["lines.linewidth"] = 2
+plt.rcParams["grid.color"] = "#CCCCCC"
+plt.rcParams["grid.linewidth"] = 0.5
+plt.rcParams["figure.edgecolor"] = "none"
+
+# PU profile is gaussian around 70
+# https://cms-pdmv-prod.web.cern.ch/mcm/public/restapi/requests/get_setup/TSG-Run3Summer22EEDR-00168
+# https://github.com/cms-sw/cmssw/blob/master/SimGeneral/MixingModule/python/mix_POISSON_average_cfi.py
+files = {
+    "Pu70": [
+        # 2022EE, v10, kl1
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv10/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p00_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson70KeepRAW_124X_mcRun3_2022_realistic_postEE_v1-v1/30000/d00363f4-0cac-410d-8fc7-bb6f60ccb6cd.root",
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv10/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p00_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson70KeepRAW_124X_mcRun3_2022_realistic_postEE_v1-v1/30000/ee9ea1b0-ee93-46b1-b2b9-164cf499ef22.root",
+    ],
+    "Pu70_1": [
+        # 2022EE, v10, kl5
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv10/GluGlutoHHto4B_kl-5p00_kt-1p00_c2-0p00_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson70KeepRAW_124X_mcRun3_2022_realistic_postEE_v1-v1/30000/ff50a677-99ec-4148-9c88-648dad33766b.root"
+    ],
+    "Pu60": [
+        # 2022EE, v12, kl1, Pu60
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv12/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p00_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson60KeepRAW_130X_mcRun3_2022_realistic_postEE_v6-v2/2540000/00d98799-ada3-4a26-8558-5052891a8d23.root",
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv12/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p00_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson60KeepRAW_130X_mcRun3_2022_realistic_postEE_v6-v2/30000/70a2f3a2-1cac-4633-bb22-9bb48ea4114c.root",
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv12/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p00_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson60KeepRAW_130X_mcRun3_2022_realistic_postEE_v6-v2/50000/fc603037-ef65-4bbf-9cef-934ecec40bbe.root",
+    ],
+    "Pu60_1": [
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv12/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p10_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson60KeepRAW_130X_mcRun3_2022_realistic_postEE_v6-v2/2540000/a8f839de-fe30-4f04-be63-3284b91f23ce.root",
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv12/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p10_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson60KeepRAW_130X_mcRun3_2022_realistic_postEE_v6-v2/40000/7bd67f96-8b9e-465a-920d-3ee9e18131f7.root",
+        "root://cmsdcadisk.fnal.gov//dcache/uscmsdisk/store/mc/Run3Summer22EENanoAODv12/GluGlutoHHto4B_kl-1p00_kt-1p00_c2-0p10_TuneCP5_13p6TeV_powheg-pythia8/NANOAODSIM/Poisson60KeepRAW_130X_mcRun3_2022_realistic_postEE_v6-v2/50000/b3ba2704-adc1-436e-b111-1dece13a5de7.root",
+    ],
+}
+
+
+def main():
+    npu_axis = hist.axis.Regular(100, 0, 100, name="npu", label=r"nPUInt")
+
+    fig, ax = plt.subplots(1, 1, figsize=(6, 4))
+
+    h_npu = hist.Hist(npu_axis)
+    bins = h_npu.to_numpy()[1]
+
+    for sample in files:
+        h_npu = hist.Hist(npu_axis)
+        for fname in files[sample]:
+            events = NanoEventsFactory.from_root(fname, schemaclass=NanoAODSchema).events()
+            # h_npu.fill(events.Pileup.nTrueInt.to_numpy())
+            # print(events.Pileup.nTrueInt.to_numpy())
+            h_npu.fill(events.Pileup.nPU.to_numpy())
+        pileup_MC = h_npu.to_numpy()[0].astype("float64")
+        # avoid division by zero
+        pileup_MC[pileup_MC == 0.0] = 1
+        # normalize
+        pileup_MC /= pileup_MC.sum()
+        print(np.round(pileup_MC, 3))
+        with Path(f"data/pileup/{sample}.npy").open("wb") as f:
+            np.save(f, pileup_MC)
+
+        # plot
+        hep.histplot(
+            h_npu / pileup_MC.sum(), histtype="step", ax=ax, label=f"MC {sample}", density=True
+        )
+
+    # plot profile for data
+    path_pileup = "data/MyDataPileupHistogram2022FG.root"
+    pileup_profile = uproot.open(path_pileup)["pileup"]
+    pileup_profile = pileup_profile.to_numpy()[0]
+    # normalise
+    pileup_profile /= pileup_profile.sum()
+    print("data ", np.round(pileup_profile, 3))
+    print("bins ", bins)
+    nomBins = bins
+    nBinCenters = 0.5 * (nomBins[:-1] + nomBins[1:])
+    ax.hist(
+        nBinCenters,
+        weights=pileup_profile,
+        bins=nomBins,
+        histtype="step",
+        label="data",
+    )
+
+    ax.legend()
+    ax.set_ylabel("Density")
+    ax.set_xlabel("PU profile")
+    fig.savefig("pileup_comparison_perfile.png")
+
+
+if __name__ == "__main__":
+    sys.exit(main())