nonres bias test

rkansal47 · Apr 9, 2024 · 8587564 · 8587564
1 parent d4028d7
commit 8587564
Show file tree

Hide file tree

Showing 8 changed files with 252 additions and 84 deletions.
diff --git a/src/HHbbVV/VBF_binder/VBFKinematicsStudyRK.ipynb b/src/HHbbVV/VBF_binder/VBFKinematicsStudyRK.ipynb
@@ -339,10 +339,10 @@
    "source": [
     "sel = ak.fill_none(\n",
     "    (\n",
-    "        (txbb[bb_mask] > 0.97)\n",
-    "        * (fatjets.particleNet_H4qvsQCD[~bb_mask] > 0.6)\n",
-    "        * (fatjets.pt[:, 0] > 500)\n",
-    "        * (fatjets.pt[:, 1] > 400)\n",
+    "        # (txbb[bb_mask] > 0.97)\n",
+    "        # * (fatjets.particleNet_H4qvsQCD[~bb_mask] > 0.6)\n",
+    "        (fatjets.pt[:, 0] > 300)\n",
+    "        * (fatjets.pt[:, 1] > 300)\n",
     "        * (np.abs(fatjets[:, 0].delta_phi(fatjets[:, 1])) > 2.6)\n",
     "        * (np.abs(fatjets[:, 0].eta - fatjets[:, 1].eta) < 2.0)\n",
     "    ),\n",

diff --git a/src/HHbbVV/combine/binder/BiasTest.ipynb b/src/HHbbVV/combine/binder/BiasTest.ipynb
diff --git a/src/HHbbVV/combine/binder/PlotScan.ipynb b/src/HHbbVV/combine/binder/PlotScan.ipynb
diff --git a/src/HHbbVV/combine/nonres_scan.sh b/src/HHbbVV/combine/nonres_scan.sh
@@ -15,30 +15,28 @@
 #   )
 # done
 
+templates_tag="24Apr8VBFBDTScanNodEta"
 
-templates_tag="24Apr2ggFk2v0Scan"
-
-# for bdt_cut in 0.99 0.997 0.998 0.999
-for bdt_cut in 0.9995 0.9999
+for bdt_cut in 0.99 0.997 0.998 0.999 0.9997 0.9999
 do
-  for txbb_cut in "HP"
+  for txbb_cut in "MP" "HP"
   do
     cutstr="txbb_${txbb_cut}_bdt_${bdt_cut}_lepton_veto_Hbb"
     echo $cutstr
 
-    python3 -u postprocessing/CreateDatacard.py --templates-dir templates/$templates_tag/$cutstr \
-    --model-name $templates_tag/SM/$cutstr --no-do-jshifts --nTF 0 --only-sm
+    # python3 -u postprocessing/CreateDatacard.py --templates-dir templates/$templates_tag/$cutstr \
+    # --model-name $templates_tag/SM/$cutstr --no-do-jshifts --nTF 0 --only-sm
 
-    (
-      cd cards/$templates_tag/SM/$cutstr || exit
-      /uscms/home/rkansal/nobackup/HHbbVV/src/HHbbVV/combine/run_blinded.sh -wbl
-    )
+    # (
+    #   cd cards/$templates_tag/SM/$cutstr || exit
+    #   /uscms/home/rkansal/nobackup/HHbbVV/src/HHbbVV/combine/run_blinded.sh -wbl
+    # )
 
     python3 -u postprocessing/CreateDatacard.py --templates-dir templates/$templates_tag/$cutstr \
-    --model-name $templates_tag/k2v0/$cutstr --no-do-jshifts --nTF 0 --sig-sample qqHH_CV_1_C2V_0_kl_1_HHbbVV
+    --model-name $templates_tag/$cutstr --no-do-jshifts --nTF 1 --sig-sample qqHH_CV_1_C2V_0_kl_1_HHbbVV
 
     (
-      cd cards/$templates_tag/k2v0/$cutstr || exit
+      cd cards/$templates_tag/$cutstr || exit
       /uscms/home/rkansal/nobackup/HHbbVV/src/HHbbVV/combine/run_blinded.sh -wbl
     )
   done

diff --git a/src/HHbbVV/combine/run_blinded.sh b/src/HHbbVV/combine/run_blinded.sh
@@ -298,7 +298,7 @@ if [ $dfit = 1 ]; then
 
     echo "Fit Shapes"
     PostFitShapesFromWorkspace --dataset "$dataset" -w ${wsm}.root --output FitShapes.root \
-    -m 125 -f fitDiagnosticsBlinded.root:fit_b --postfit --print 2>&1 | tee $outsdir/FitShapes.txt
+    -m 125 -f fitDiagnosticsBlinded.root:fit_b --postfit --sampling --print 2>&1 | tee $outsdir/FitShapes.txt
 fi
 
 
@@ -395,7 +395,7 @@ if [ "$bias" != -1 ]; then
     combine -M FitDiagnostics --trackParameters r --trackErrors r --justFit \
     -m 125 -n "bias${bias}" -d ${wsm_snapshot}.root --rMin "-15" --rMax 15 \
     --snapshotName MultiDimFit --bypassFrequentistFit --toysFrequentist --expectSignal "$bias" \
-    "${unblindedparams},r=$bias" --floatParameters "${freezeparamsblinded}" \
+    ${unblindedparams},r=$bias --floatParameters ${freezeparamsblinded} \
     --robustFit=1 -t "$numtoys" -s "$seed" \
     --X-rtd MINIMIZER_MaxCalls=1000000 --cminDefaultMinimizerTolerance "$mintol" 2>&1 | tee "$outsdir/bias${bias}seed${seed}.txt"
 fi
diff --git a/src/HHbbVV/combine/submit_bias_nonres_loop.sh b/src/HHbbVV/combine/submit_bias_nonres_loop.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+# shellcheck disable=SC2043
+
+#######################################################################
+# Script to submit bias jobs, needs to be run from inside the datacards directory
+#####################################################
+
+seed=$1
+TAG=$2
+
+for bias in 0.0 0.3 1.0 1.5
+do
+  python3 /uscms_data/d1/rkansal/HHbbVV/src/HHbbVV/combine/submit/submit_bias.py --seed "$seed" --num-jobs 100 --toys-per-job 10 --bias "$bias" --submit --tag "$TAG" --mintol 20
+done
+
+# # need to submit extra jobs for these because of high fit failures
+# sample=NMSSM_XToYHTo2W2BTo4Q2B_MX-3000_MY-250
+# cd $sample
+# bias=0.0
+# python3 /uscms_data/d1/rkansal/HHbbVV/src/HHbbVV/combine/submit/submit_bias.py --seed $((seed + 1000)) --num-jobs 100 --toys-per-job 10 --bias $bias --submit --tag $TAG
+
+# bias=0.15
+# python3 /uscms_data/d1/rkansal/HHbbVV/src/HHbbVV/combine/submit/submit_bias.py --seed $((seed + 1000)) --num-jobs 50 --toys-per-job 10 --bias $bias --submit --tag $TAG
+# cd -
diff --git a/src/HHbbVV/postprocessing/PlotFits.ipynb b/src/HHbbVV/postprocessing/PlotFits.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -26,7 +26,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -36,13 +36,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "MAIN_DIR = Path(\"../../../\")\n",
     "\n",
-    "plot_dir = MAIN_DIR / \"plots/PostFit/24Apr3NonresScan\"\n",
+    "# plot_dir = MAIN_DIR / \"plots/PostFit/24Apr4VBFBDTScan/txbb_MP_bdt_0.9999_lepton_veto_Hbb/\"\n",
+    "plot_dir = MAIN_DIR / \"plots/PostFit/24Apr8Uncs1st/\"\n",
     "plot_dir.mkdir(exist_ok=True, parents=True)"
    ]
   },
@@ -52,8 +53,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cards_dir = \"24Mar15UpdateData\"\n",
-    "cards_dir = \"24Apr2ggFk2v0Scan/SM/txbb_HP_bdt_0.9995_lepton_veto_Hbb\"\n",
+    "cards_dir = \"24Mar15UpdateDatanTF1\"\n",
+    "# cards_dir = \"24Apr4VBFBDTScan/k2v0/txbb_MP_bdt_0.9999_lepton_veto_Hbb\"\n",
     "asimov = False\n",
     "\n",
     "asimov_label = \"Asimov\" if asimov else \"\"\n",
@@ -69,7 +70,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "templates_dir = Path(\"templates/24Mar15UpdateData\")\n",
+    "templates_dir = Path(f\"templates/{cards_dir[:-4]}\")\n",
+    "# templates_dir = Path(f\"templates/24Apr4VBFBDTScan/txbb_MP_bdt_0.9999_lepton_veto_Hbb\")\n",
     "templates_dict = {}\n",
     "for year in years:\n",
     "    with (templates_dir / f\"{year}_templates.pkl\").open(\"rb\") as f:\n",
@@ -103,6 +105,7 @@
     "hist_label_map = {val: key for key, val in hist_label_map_inverse.items()}\n",
     "\n",
     "sig_keys = [\"HHbbVV\", \"VBFHHbbVV\"]\n",
+    "sig_keys = [\"qqHH_CV_1_C2V_0_kl_1_HHbbVV\"]\n",
     "samples = bg_keys + sig_keys + [data_key]"
    ]
   },
@@ -133,6 +136,7 @@
    "outputs": [],
    "source": [
     "hists = {}\n",
+    "bgerrs = {}\n",
     "\n",
     "for shape in shapes:\n",
     "    print(shape)\n",
@@ -144,6 +148,7 @@
     "        )\n",
     "        for region in selection_regions\n",
     "    }\n",
+    "    bgerrs[shape] = {}\n",
     "\n",
     "    for region in selection_regions:\n",
     "        h = hists[shape][region]\n",
@@ -165,7 +170,9 @@
     "        data_key_index = np.where(np.array(list(h.axes[0])) == data_key)[0][0]\n",
     "        h.view(flow=False)[data_key_index, :] = np.nan_to_num(\n",
     "            templates[hist_label_map_inverse[data_key]].values()\n",
-    "        )"
+    "        )\n",
+    "\n",
+    "        bgerrs[shape][region] = templates[\"TotalBkg\"].errors()"
    ]
   },
   {
@@ -174,9 +181,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pass_ylim = 22\n",
+    "pass_ylim = 17\n",
     "fail_ylim = 600000\n",
     "title_label = \" Asimov Dataset\" if asimov else \"\"\n",
+    "sig_scale_dict = {\"HHbbVV\": 100, \"VBFHHbbVV\": 2000, \"qqHH_CV_1_C2V_0_kl_1_HHbbVV\": 1}\n",
+    "sig_scale_dict = {key: val for key, val in sig_scale_dict.items() if key in sig_keys}\n",
+    "\n",
     "for shape, shape_label in shapes.items():\n",
     "    for region, region_label in selection_regions.items():\n",
     "        pass_region = region.startswith(\"pass\")\n",
@@ -185,7 +195,8 @@
     "                \"hists\": hists[shape][region],\n",
     "                \"sig_keys\": sig_keys,\n",
     "                \"bg_keys\": bg_keys,\n",
-    "                \"sig_scale_dict\": {\"HHbbVV\": 100, \"VBFHHbbVV\": 2000} if pass_region else None,\n",
+    "                \"bg_err\": bgerrs[shape][region],\n",
+    "                \"sig_scale_dict\": sig_scale_dict if pass_region else None,\n",
     "                \"show\": True,\n",
     "                \"year\": \"all\",\n",
     "                \"ylim\": pass_ylim if pass_region else fail_ylim,\n",
@@ -204,7 +215,51 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "hists[\"postfit\"][\"pass\"][\"QCD\", ...].values()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "edges = hists[\"postfit\"][\"pass\"].axes[1].edges\n",
+    "mps = (edges[1:] + edges[:-1]) / 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "errratio = bgerrs[\"postfit\"][\"pass\"] / sum(\n",
+    "    [hists[\"postfit\"][\"pass\"][sample, :] for sample in bg_keys]\n",
+    ")\n",
+    "plt.scatter(mps, errratio)\n",
+    "plt.ylabel(\"Error / Total\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sum(hists[\"postfit\"][\"pass\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sum([hists[\"postfit\"][\"pass\"][sample, :] for sample in bg_keys])"
+   ]
   }
  ],
  "metadata": {

diff --git a/src/HHbbVV/postprocessing/plotting.py b/src/HHbbVV/postprocessing/plotting.py
@@ -383,6 +383,13 @@ def ratioHistPlot(
             )
 
     if bg_err is not None:
+        if divide_bin_width:
+            raise NotImplementedError("Background error for divide bin width not checked yet")
+
+        if len(np.array(bg_err).shape) == 1:
+            bg_tot = sum([pre_divide_hists[sample, :] for sample in bg_keys])
+            bg_err = [bg_tot - bg_err, bg_tot + bg_err]
+
         ax.fill_between(
             np.repeat(hists.axes[1].edges, 2)[1:-1],
             np.repeat(bg_err[0], 2),