dist studies

src/HHbbVV/postprocessing/TopAnalysis.ipynb

@@ -14,7 +14,7 @@
     "from pathlib import Path\n",
     "\n",
     "from pandas.errors import SettingWithCopyWarning\n",
-    "from HHbbVV.hh_vars import data_key\n",
+    "from HHbbVV.hh_vars import data_key, years\n",
     "import postprocessing\n",
     "\n",
     "# ignore these because they don't seem to apply\n",
@@ -51,7 +51,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot_dir = Path(\"../../../plots/ttsfs/24Jul24BLDistortion\")\n",
+    "plot_dir = Path(\"../../../plots/ttsfs/24Jul28Distortion\")\n",
     "plot_dir.mkdir(parents=True, exist_ok=True)"
    ]
   },
@@ -362,55 +362,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Tesing distortion uncertainty"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pickle\n",
-    "\n",
-    "# package_path = Path(__file__).parent.parent.resolve()\n",
-    "package_path = Path(\"../\").resolve()\n",
-    "\n",
-    "# for key in [\n",
-    "#     \"TTToSemiLeptonic\",\n",
-    "#     # \"ST_tW_antitop_5f_NoFullyHadronicDecays\",\n",
-    "#     # \"ST_tW_top_5f_NoFullyHadronicDecays\",\n",
-    "#     # \"ST_s-channel_4f_leptonDecays\",\n",
-    "#     # \"ST_t-channel_antitop_4f_InclusiveDecays\",\n",
-    "#     # \"ST_t-channel_top_4f_InclusiveDecays\",\n",
-    "# ]:\n",
-    "#     sig_lp_hist = utils.get_pickles(f\"{signal_data_dir}/{year}/{key}/pickles\", year, key)[\"lp_hist\"]\n",
-    "\n",
-    "#     # remove negatives\n",
-    "#     sig_lp_hist.values()[sig_lp_hist.values() < 0] = 0\n",
-    "\n",
-    "#     with (package_path / f\"corrections/lp_ratios/signals/{year}_{key}.hist\").open(\"wb\") as f:\n",
-    "#         pickle.dump(sig_lp_hist, f)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sig_lp_hist.values()[sig_lp_hist.values() <= 0] = 1"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sig_mc_ratio = np.nan_to_num((sig_lp_hist.values() / sig_tot) / (mc_nom / mc_tot), nan=1.0)\n",
-    "sig_mc_ratio[sig_mc_ratio == 0] = 1.0\n",
-    "sig_mc_ratio = np.clip(sig_mc_ratio, 0.5, 2.0)"
+    "### Testing distortion uncertainty"
    ]
   },
   {
@@ -420,86 +372,45 @@
    "outputs": [],
    "source": [
     "import uproot\n",
+    "from tqdm import tqdm\n",
     "\n",
-    "# initialize lund plane scale factors lookups\n",
-    "# f = uproot.open(package_path / f\"corrections/lp_ratios/ratio_{year[:4]}.root\")\n",
-    "f = uproot.open(package_path / f\"corrections/lp_ratios/ratio_2016.root\")\n",
+    "for dist_year in tqdm(years[1:]):\n",
+    "    f = uproot.open(package_path / f\"corrections/lp_ratios/ratio_{dist_year}.root\")\n",
     "\n",
-    "# 3D histogram: [subjet_pt, ln(0.8/Delta), ln(kT/GeV)]\n",
-    "mc_nom = f[\"mc_nom\"].to_numpy()\n",
-    "ratio_edges = mc_nom[1:]\n",
-    "mc_nom = mc_nom[0]\n",
+    "    # 3D histogram: [subjet_pt, ln(0.8/Delta), ln(kT/GeV)]\n",
+    "    mc_nom = f[\"mc_nom\"].to_numpy()\n",
+    "    ratio_edges = mc_nom[1:]\n",
+    "    mc_nom = mc_nom[0]\n",
     "\n",
-    "ratio_nom = f[\"ratio_nom\"].to_numpy()[0]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sig_old_lp_hist = utils.get_pickles(\n",
-    "    f\"../../../temp/2016/GluGluToHHTobbVV_node_cHHH1/pickles\", \"2016\", \"GluGluToHHTobbVV_node_cHHH1\"\n",
-    ")[\"lp_hist\"]\n",
-    "sig_lp_hist = utils.get_pickles(f\"../../../outfiles\", \"2016\", \"GluGluToHHTobbVV_node_cHHH1\")[\n",
-    "    \"lp_hist\"\n",
-    "]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "with (package_path / f\"corrections/lp_ratios/signals/2018_GluGluToHHTobbVV_node_cHHH1.hist\").open(\n",
-    "    \"rb\"\n",
-    ") as f:\n",
-    "    sig_lp_hist = pickle.load(f)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# plotting.plot_lund_plane(sig_old_lp_hist.values()[0], show=True)\n",
+    "    mc_tot_pt = np.sum(mc_nom, axis=(1, 2), keepdims=True)\n",
+    "    mc_density = mc_nom / mc_tot_pt\n",
+    "    plotting.plot_lund_plane_six(\n",
+    "        mc_density, ratio_edges, name=f\"{plot_dir}/{dist_year}_MC.pdf\", show=False\n",
+    "    )\n",
     "\n",
-    "# extent = [-1, 8, -5, 7]\n",
-    "# plt.figure(figsize=(12, 13))\n",
-    "# plt.imshow(sig_old_lp_hist.values()[0].T[::-1], extent=extent)\n",
-    "# plt.colorbar()\n",
-    "# plt.show()\n",
+    "    # ratio_nom = f[\"ratio_nom\"].to_numpy()[0]\n",
     "\n",
-    "plt.subplots(\n",
-    "    1,\n",
-    "    1,\n",
-    ")\n",
-    "plt.imshow(sig_old_lp_hist.values()[0])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.subplots(\n",
-    "    1,\n",
-    "    1,\n",
-    ")\n",
-    "plt.imshow(sig_lp_hist.values()[0])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.imshow(mc_nom[0])"
+    "    for sig in [\"GluGluToHHTobbVV_node_cHHH1\", \"VBF_HHTobbVV_CV_1_C2V_2_C3_1\"]:\n",
+    "        with (package_path / f\"corrections/lp_ratios/signals/{dist_year}_{sig}.hist\").open(\n",
+    "            \"rb\"\n",
+    "        ) as f:\n",
+    "            sig_lp_hist = pickle.load(f)\n",
+    "            sig_tot_pt = np.sum(sig_lp_hist.values(), axis=(1, 2), keepdims=True)\n",
+    "            sig_density = sig_lp_hist.values() / sig_tot_pt\n",
+    "\n",
+    "            mc_sig_ratio = np.nan_to_num(mc_density / sig_density, nan=1.0)\n",
+    "            mc_sig_ratio[mc_sig_ratio == 0] = 1.0\n",
+    "            mc_sig_ratio = np.clip(mc_sig_ratio, 0.2, 5.0)\n",
+    "\n",
+    "            plotting.plot_lund_plane_six(\n",
+    "                sig_density, ratio_edges, name=f\"{plot_dir}/{dist_year}_{sig}.pdf\", show=False\n",
+    "            )\n",
+    "            plotting.plot_lund_plane_six(\n",
+    "                mc_sig_ratio,\n",
+    "                ratio_edges,\n",
+    "                name=f\"{plot_dir}/{dist_year}_{sig}_ratio.pdf\",\n",
+    "                show=False,\n",
+    "            )"
    ]
   },
   {

src/HHbbVV/postprocessing/plotting.py

@@ -31,9 +31,9 @@
 formatter.set_powerlimits((-3, 3))
 
 # this is needed for some reason to update the font size for the first plot
-fig, ax = plt.subplots(1, 1, figsize=(12, 12))
-plt.rcParams.update({"font.size": 24})
-plt.close()
+# fig, ax = plt.subplots(1, 1, figsize=(12, 12))
+# plt.rcParams.update({"font.size": 24})
+# plt.close()
 
 
 bg_order = ["Diboson", "HH", "HWW", "Hbb", "ST", "W+Jets", "Z+Jets", "TT", "QCD"]
@@ -1274,13 +1274,27 @@ def cutsLinePlot(
         plt.close()
 
 
-def plot_lund_plane(h: np.ndarray, ax=None, name: str = None, show: bool = False):
+def plot_lund_plane(
+    h: np.ndarray, title: str = "", ax=None, fig=None, name: str = "", show: bool = False
+):
     if ax is None:
-        fig, ax = plt.subplots(figsize=(12, 14))
+        fig, ax = plt.subplots(figsize=(10, 8))
+    else:
+        assert fig is not None, "Must provide fig if providing ax."
 
     extent = [-1, 8, -5, 7]
-    ax.imshow(h.T[::-1], extent=extent)
-    ax.colorbar()
+    im = ax.imshow(h.T, origin="lower", extent=extent, cmap="viridis")
+    ax.set_aspect("auto")
+    fig.colorbar(im, ax=ax)
+    # cbar.set_label('Density')
+
+    ax.set_xlabel(r"ln$(0.8/\Delta)$")
+    ax.set_ylabel(r"ln$(k_T/GeV)$")
+
+    if len(title):
+        ax.set_title(title, fontsize=24)
+
+    plt.tight_layout()
 
     if ax is None:
         if len(name):
@@ -1290,3 +1304,22 @@ def plot_lund_plane(h: np.ndarray, ax=None, name: str = None, show: bool = False
             plt.show()
         else:
             plt.close()
+
+
+def plot_lund_plane_six(hists: np.ndarray, edges: np.ndarray, name: str = "", show: bool = False):
+    fig, axs = plt.subplots(2, 3, figsize=(20, 12), sharex=True, sharey=True)
+    for i, ax in enumerate(axs.flat):
+        plot_lund_plane(
+            hists[i],
+            title=rf"$p_T$: [{edges[0][i]:.0f}, {edges[0][i + 1]:.0f}] GeV",
+            ax=ax,
+            fig=fig,
+        )
+
+    if len(name):
+        plt.savefig(name, bbox_inches="tight")
+
+    if show:
+        plt.show()
+    else:
+        plt.close()

src/HHbbVV/scale_factors/accumulate_lp_hists.py

@@ -25,14 +25,18 @@
 for year in years:
     print(year)
     for key in tqdm(nonres_samples.values()):
-        sig_lp_hist = utils.get_pickles(f"{args.data_path}/{year}/{key}/pickles", year, key)[
-            "lp_hist"
-        ]
+        try:
+            sig_lp_hist = utils.get_pickles(f"{args.data_path}/{year}/{key}/pickles", year, key)[
+                "lp_hist"
+            ]
+        except FileNotFoundError as e:
+            print(e)
+            continue
 
         # remove negatives
         sig_lp_hist.values()[sig_lp_hist.values() < 0] = 0
 
         with (package_path / f"corrections/lp_ratios/signals/{year}_{key}.hist").open("wb") as f:
             pickle.dump(sig_lp_hist, f)
 
-        break
+        # break

src/condor/check_jobs.py

@@ -110,6 +110,22 @@
 
     if not Path(f"{eosdir}/{sample}/pickles").exists():
         print_red(f"No pickles directory for {sample}!")
+        if not args.check_parquet:
+            for i in range(jdl_dict[sample]):
+                if f"{args.year}_{sample}_{i}" in running_jobs:
+                    print(f"Job #{i} for sample {sample} is running.")
+                    continue
+
+                jdl_file = (
+                    f"{user_condor_dir}/{args.processor}/{args.tag}/{args.year}_{sample}_{i}.jdl"
+                )
+                err_file = f"{user_condor_dir}/{args.processor}/{args.tag}/logs/{args.year}_{sample}_{i}.err"
+                print(jdl_file)
+                missing_files.append(jdl_file)
+                err_files.append(err_file)
+                if args.submit_missing:
+                    os.system(f"condor_submit {jdl_file}")
+
         continue
 
     outs_pickles = [