{ "cells": [ { "cell_type": "markdown", "id": "7d57cbdf", "metadata": {}, "source": [ "### Implement packages and Load Data" ] }, { "cell_type": "code", "execution_count": 1, "id": "ce0fa7b6", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/billyli/miniforge_x86_new/lib/python3.9/site-packages/coffea/util.py:154: FutureWarning: In coffea version v0.8.0 (target date: 31 Dec 2022), this will be an error.\n", "(Set coffea.deprecations_as_errors = True to get a stack trace now.)\n", "ImportError: coffea.hist is deprecated\n", " warnings.warn(message, FutureWarning)\n" ] } ], "source": [ "import itertools\n", "import logging\n", "from pathlib import Path\n", "import numba as nb\n", "\n", "import awkward as ak\n", "import click\n", "import h5py as h5\n", "import numpy as np\n", "import vector\n", "\n", "from coffea.hist.plot import clopper_pearson_interval\n", "import matplotlib.pyplot as plt\n", "\n", "# from src.data.cms.convert_to_h5 import MIN_JETS, N_JETS, N_FJETS\n", "\n", "vector.register_awkward()\n", "\n", "logging.basicConfig(level=logging.INFO)" ] }, { "cell_type": "code", "execution_count": 2, "id": "72060c9f", "metadata": {}, "outputs": [], "source": [ "# read test target file\n", "test_file = \"//Users/billyli/UCSD/hhh/reports/bv2/hhh_test.h5\"\n", "test_h5 = h5.File(test_file)\n", "\n", "# read baseline prediction\n", "baseline_file = \"//Users/billyli/UCSD/hhh/reports/bv2/pred_baseline.h5\"\n", "b_h5 = h5.File(baseline_file)\n", "\n", "# read spanet prediction\n", "spanet_file = \"//Users/billyli/UCSD/hhh/reports/bv2/dp_on/pred_v53.h5\"\n", "s_h5 = h5.File(spanet_file)" ] }, { "cell_type": "markdown", "id": "c93aa3a9", "metadata": {}, "source": [ "### Abstract\n", "Overlapped-jets removal is motivated by the fact that some resolved Higgs (rH) can have the corresponding boosted Higgs (bH). As those rH and bH are two different representations of the same Higgs, we should only choose one of them to avoid double counting when evaluating purity/efficiency. As a natural solution to prioritize bH, overlapped-jets removal is to remove the jets overlapped with the bH by a delta_R condiction. One will calculate dR between each jet and the bH candidate. The jets that have a dR less than dR_min will be removed. The removal is carried out after reconstructing bH as we prioritize the bH reconstruction.\n", "\n", "It is possible that some rH jets having no bH correspondance but still removed if we set the delta_R_min too large. On the other hand, some bHs' corresponding rHs' jets could be left remaining if the delta_R_min is too small.\n", "\n", "The goal of this notebook is not to pick the optimal deltaR_min, but to have a grasp of how well the default deltaR_min = 0.8 does." ] }, { "cell_type": "markdown", "id": "c6098ee3", "metadata": {}, "source": [ "### Methods" ] }, { "cell_type": "markdown", "id": "6fbad25c", "metadata": {}, "source": [ "### Gen-level Overlap Removal Study\n", "Below is to provide some info about how often a resolved Higgs and a boosted Higgs overlap\n", "\n", "* N1: number of resolved H without overlapping removal\n", "* N2: number of resolved H with overlapping removal\n", "* N3: number of Higgs that are bothe resolved and boosted\n", "\n", "N1-N2-N3 = Number of Higgs that were removed - Number of Higgs that are both resolved and boosted\n", "\n", "Also we can generate a confusion matrix of bi-cat (a higgs that is both bH and rH) and removal. " ] }, { "cell_type": "markdown", "id": "48831e09", "metadata": {}, "source": [ "### General Outline\n", "\n", "1. Load functions to reco bH\n", "2. Load functions to remove overlap\n", "3. Find the necessary information" ] }, { "cell_type": "markdown", "id": "e9642e14", "metadata": {}, "source": [ "### 1." ] }, { "cell_type": "code", "execution_count": 3, "id": "f70c30f1", "metadata": {}, "outputs": [], "source": [ "def sel_pred_bH_by_dp(dps, aps, bb_ps, dp_cut, ap_cut=1/13):\n", " # parse predicted bb assignment by DP\n", " dp_filter = dps>dp_cut\n", " ap_filter = aps>ap_cut\n", " ak8_filter = bb_ps>9\n", " filter = dp_filter&ak8_filter\n", " \n", " bb_ps_passed = bb_ps.mask[filter]\n", " bb_ps_passed = ak.drop_none(bb_ps_passed)\n", " \n", " return bb_ps_passed" ] }, { "cell_type": "code", "execution_count": 4, "id": "b49cc7ba", "metadata": {}, "outputs": [], "source": [ "def sel_target_bH_by_mask(bb_ts, bh_pts, bh_masks):\n", " bb_ts_selected = bb_ts.mask[bh_masks]\n", " bb_ts_selected = ak.drop_none(bb_ts_selected)\n", " \n", " bh_selected_pts = bh_pts.mask[bh_masks]\n", " bh_selected_pts = ak.drop_none(bh_selected_pts)\n", " \n", " return bb_ts_selected, bh_selected_pts" ] }, { "cell_type": "code", "execution_count": 5, "id": "5bb0057d", "metadata": {}, "outputs": [], "source": [ "# A pred look up table is in shape\n", "# [event,\n", "# pred_H, \n", "# [correct, pred_H_pt]]\n", "def gen_pred_bH_LUT(bb_ps_passed, bb_ts_selected, fj_pts):\n", " LUT = []\n", " # for each event\n", " for bb_t_event, bb_p_event, fj_pt_event in zip(bb_ts_selected, bb_ps_passed, fj_pts):\n", " # for each predicted bb assignment, check if any target H have a same bb assignment\n", " LUT_event = []\n", " for i, bb_p in enumerate(bb_p_event):\n", " correct = 0\n", " predH_pt = fj_pt_event[bb_p-10]\n", " for bb_t in bb_t_event:\n", " if bb_p == bb_t+10:\n", " correct = 1\n", " LUT_event.append([correct, predH_pt])\n", " LUT.append(LUT_event)\n", " return LUT" ] }, { "cell_type": "code", "execution_count": 6, "id": "fb653736", "metadata": {}, "outputs": [], "source": [ "# A target look up table is in shape\n", "# [event,\n", "# target_H, \n", "# target_bb_assign,\n", "# [retrieved, targetH_pt]]\n", "def gen_target_bH_LUT(bb_ps_passed, bb_ts_selected, targetH_pts):\n", " LUT = []\n", " # for each event\n", " for bb_t_event, bb_p_event, targetH_pts_event in zip(bb_ts_selected, bb_ps_passed, targetH_pts):\n", " # for each target fatjet, check if the predictions have a p fatject same with the t fatjet\n", " LUT_event = []\n", " for i, bb_t in enumerate(bb_t_event):\n", " retrieved = 0\n", " targetH_pt = targetH_pts_event[i]\n", " for bb_p in bb_p_event:\n", " if bb_p == bb_t+10:\n", " retrieved = 1\n", " LUT_event.append([retrieved, targetH_pt])\n", " LUT.append(LUT_event)\n", " return LUT" ] }, { "cell_type": "code", "execution_count": 7, "id": "d436a2aa", "metadata": {}, "outputs": [], "source": [ "# generate pred/target LUT\n", "# each entry corresponds to [recoH correct or not, reco H pt]\n", "# or \n", "# [targetH retrieved or not, target H pt]\n", "def parse_boosted_w_target(testfile, predfile, dp_cut=0.8):\n", " # Collect H pt, mask, target and predicted jet and fjets for 3 Hs in each event\n", " # h pt\n", " bh1_pt = np.array(testfile['TARGETS']['bh1']['pt'])\n", " bh2_pt = np.array(testfile['TARGETS']['bh2']['pt'])\n", " bh3_pt = np.array(testfile['TARGETS']['bh3']['pt'])\n", "\n", " # mask\n", " bh1_mask = np.array(testfile['TARGETS']['bh1']['mask'])\n", " bh2_mask = np.array(testfile['TARGETS']['bh2']['mask'])\n", " bh3_mask = np.array(testfile['TARGETS']['bh3']['mask'])\n", "\n", " # target assignment\n", " bb_bh1_t = np.array(testfile[\"TARGETS\"][\"bh1\"]['bb'])\n", " bb_bh2_t = np.array(testfile[\"TARGETS\"][\"bh2\"]['bb'])\n", " bb_bh3_t = np.array(testfile[\"TARGETS\"][\"bh3\"]['bb'])\n", "\n", " try:\n", " # pred assignment\n", " bb_bh1_p = np.array(predfile[\"TARGETS\"][\"bh1\"]['bb'])\n", " bb_bh2_p = np.array(predfile[\"TARGETS\"][\"bh2\"]['bb'])\n", " bb_bh3_p = np.array(predfile[\"TARGETS\"][\"bh3\"]['bb'])\n", " \n", " # boosted Higgs detection probability\n", " dp_bh1 = np.array(predfile[\"TARGETS\"][\"bh1\"]['detection_probability'])\n", " dp_bh2 = np.array(predfile[\"TARGETS\"][\"bh2\"]['detection_probability'])\n", " dp_bh3 = np.array(predfile[\"TARGETS\"][\"bh3\"]['detection_probability'])\n", "\n", " # fatjet assignment probability\n", " ap_bh1 = np.array(predfile[\"TARGETS\"][\"bh1\"]['assignment_probability'])\n", " ap_bh2 = np.array(predfile[\"TARGETS\"][\"bh2\"]['assignment_probability'])\n", " ap_bh3 = np.array(predfile[\"TARGETS\"][\"bh3\"]['assignment_probability'])\n", " except:\n", " # pred assignment\n", " bb_bh1_p = np.array(predfile[\"TARGETS\"][\"bh1\"]['bb'])+10\n", " bb_bh2_p = np.array(predfile[\"TARGETS\"][\"bh2\"]['bb'])+10\n", " bb_bh3_p = np.array(predfile[\"TARGETS\"][\"bh3\"]['bb'])+10\n", " \n", " # boosted Higgs detection probability\n", " dp_bh1 = np.array(predfile[\"TARGETS\"][\"bh1\"]['mask']).astype('float')\n", " dp_bh2 = np.array(predfile[\"TARGETS\"][\"bh2\"]['mask']).astype('float')\n", " dp_bh3 = np.array(predfile[\"TARGETS\"][\"bh3\"]['mask']).astype('float')\n", "\n", " # fatjet assignment probability\n", " ap_bh1 = np.array(predfile[\"TARGETS\"][\"bh1\"]['mask']).astype('float')\n", " ap_bh2 = np.array(predfile[\"TARGETS\"][\"bh2\"]['mask']).astype('float')\n", " ap_bh3 = np.array(predfile[\"TARGETS\"][\"bh3\"]['mask']).astype('float')\n", " \n", " # collect fatjet pt\n", " fj_pt = np.array(testfile['INPUTS']['BoostedJets']['fj_pt'])\n", " \n", " # convert some arrays to ak array\n", " dps = np.concatenate((dp_bh1.reshape(-1, 1), dp_bh2.reshape(-1, 1), dp_bh3.reshape(-1, 1)), axis=1)\n", " dps = ak.Array(dps)\n", " aps = np.concatenate((ap_bh1.reshape(-1, 1), ap_bh2.reshape(-1, 1), ap_bh3.reshape(-1, 1)), axis=1)\n", " aps = ak.Array(aps)\n", " bb_ps = np.concatenate((bb_bh1_p.reshape(-1, 1), bb_bh2_p.reshape(-1, 1), bb_bh3_p.reshape(-1, 1)), axis=1)\n", " bb_ps = ak.Array(bb_ps)\n", " bb_ts = np.concatenate((bb_bh1_t.reshape(-1, 1), bb_bh2_t.reshape(-1, 1), bb_bh3_t.reshape(-1, 1)), axis=1)\n", " bb_ts = ak.Array(bb_ts)\n", " fj_pt = ak.Array(fj_pt)\n", " bh_masks = np.concatenate((bh1_mask.reshape(-1, 1), bh2_mask.reshape(-1, 1), bh3_mask.reshape(-1, 1)), axis=1)\n", " bh_masks = ak.Array(bh_masks)\n", " bh_pts = np.concatenate((bh1_pt.reshape(-1, 1), bh2_pt.reshape(-1, 1), bh3_pt.reshape(-1, 1)), axis=1)\n", " bh_pts = ak.Array(bh_pts)\n", " \n", " # select predictions and targets\n", " bb_ts_selected, targetH_selected_pts = sel_target_bH_by_mask(bb_ts, bh_pts, bh_masks)\n", " bb_ps_selected = sel_pred_bH_by_dp(dps, aps, bb_ps, dp_cut)\n", " \n", " # generate correct/retrieved LUT for pred/target respectively\n", " LUT_pred = gen_pred_bH_LUT(bb_ps_selected, bb_ts_selected, fj_pt)\n", " LUT_target = gen_target_bH_LUT(bb_ps_selected, bb_ts_selected, targetH_selected_pts)\n", " \n", " # reconstruct bH to remove overlapped ak4 jets\n", " fj_eta = np.array(testfile['INPUTS']['BoostedJets']['fj_eta'])\n", " fj_phi = np.array(testfile['INPUTS']['BoostedJets']['fj_phi'])\n", " fj_mass = np.array(testfile['INPUTS']['BoostedJets']['fj_mass'])\n", " \n", " fjs = ak.zip(\n", " {\n", " \"pt\": fj_pt,\n", " \"eta\": fj_eta,\n", " \"phi\": fj_phi,\n", " \"mass\": fj_mass,\n", " },\n", " with_name=\"Momentum4D\"\n", " )\n", " fj_reco = fjs[bb_ps_selected-10]\n", " \n", " return LUT_pred, LUT_target, fj_reco" ] }, { "cell_type": "markdown", "id": "df17aea2", "metadata": {}, "source": [ "### 2." ] }, { "cell_type": "code", "execution_count": 8, "id": "9e92b322", "metadata": {}, "outputs": [], "source": [ "def get_unoverlapped_jet_index(fjs, js, dR_min=0.8):\n", " overlapped = ak.sum(js[:, np.newaxis].deltaR(fjs)0\n", " jet_index_passed = ak.local_index(js).mask[~overlapped]\n", " jet_index_passed = ak.drop_none(jet_index_passed)\n", " return jet_index_passed" ] }, { "cell_type": "markdown", "id": "3fdf6e32", "metadata": {}, "source": [ "### 3." ] }, { "cell_type": "markdown", "id": "b85f6680", "metadata": {}, "source": [ "#### Getting N1" ] }, { "cell_type": "code", "execution_count": 9, "id": "debd2c80", "metadata": {}, "outputs": [], "source": [ "def get_unmasked_arr(arr, mask):\n", " return ak.drop_none(arr.mask[mask])" ] }, { "cell_type": "code", "execution_count": 10, "id": "e042b235", "metadata": {}, "outputs": [], "source": [ "testfile = test_h5\n", "# N1\n", "h1_pt = ak.Array(testfile['TARGETS']['h1']['pt'])\n", "h2_pt = ak.Array(testfile['TARGETS']['h2']['pt'])\n", "h3_pt = ak.Array(testfile['TARGETS']['h3']['pt'])\n", "\n", "# mask\n", "h1_mask = ak.Array(testfile['TARGETS']['h1']['mask'])\n", "h2_mask = ak.Array(testfile['TARGETS']['h2']['mask'])\n", "h3_mask = ak.Array(testfile['TARGETS']['h3']['mask'])\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "cd6cdf50", "metadata": {}, "outputs": [], "source": [ "N1 = ak.count(get_unmasked_arr(h1_pt, h1_mask))+ak.count(get_unmasked_arr(h2_pt, h2_mask))+ak.count(get_unmasked_arr(h3_pt, h3_mask))" ] }, { "cell_type": "markdown", "id": "95a0dce8", "metadata": {}, "source": [ "#### Getting N2" ] }, { "cell_type": "code", "execution_count": 12, "id": "b53b9b62", "metadata": {}, "outputs": [], "source": [ "# assume perfect reconstruction\n", "_, _, fjs_reco = parse_boosted_w_target(test_h5, test_h5)" ] }, { "cell_type": "code", "execution_count": 13, "id": "d99c0a16", "metadata": {}, "outputs": [], "source": [ "# target assignments\n", "b1_h1_t = np.array(testfile[\"TARGETS\"][\"h1\"]['b1']).astype('int')\n", "b1_h2_t = np.array(testfile[\"TARGETS\"][\"h2\"]['b1']).astype('int')\n", "b1_h3_t = np.array(testfile[\"TARGETS\"][\"h3\"]['b1']).astype('int')\n", "\n", "b2_h1_t = np.array(testfile[\"TARGETS\"][\"h1\"]['b2']).astype('int')\n", "b2_h2_t = np.array(testfile[\"TARGETS\"][\"h2\"]['b2']).astype('int')\n", "b2_h3_t = np.array(testfile[\"TARGETS\"][\"h3\"]['b2']).astype('int')\n", "\n", "\n", "b1_ts = np.concatenate((b1_h1_t.reshape(-1, 1), b1_h2_t.reshape(-1, 1), b1_h3_t.reshape(-1, 1)), axis=1)\n", "b1_ts = ak.Array(b1_ts)\n", "\n", "b2_ts = np.concatenate((b2_h1_t.reshape(-1, 1), b2_h2_t.reshape(-1, 1), b2_h3_t.reshape(-1, 1)), axis=1)\n", "b2_ts = ak.Array(b2_ts)\n", "\n", "# reconstruct jet 4-momentum objects\n", "j_pt = np.array(testfile['INPUTS']['Jets']['pt'])\n", "j_eta = np.array(testfile['INPUTS']['Jets']['eta'])\n", "j_phi = np.array(testfile['INPUTS']['Jets']['phi'])\n", "j_mass = np.array(testfile['INPUTS']['Jets']['mass'])\n", "js = ak.zip(\n", " {\n", " \"pt\": j_pt,\n", " \"eta\": j_eta,\n", " \"phi\": j_phi,\n", " \"mass\": j_mass,\n", " },\n", " with_name=\"Momentum4D\"\n", ")\n", " \n", "goodJetIdx = get_unoverlapped_jet_index(fjs_reco, js, dR_min=0.8)\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "662317c0", "metadata": {}, "outputs": [], "source": [ "N2 = 0\n", "for tb1_e, tb2_e, goodJetIdx_e in zip(b1_ts, b2_ts, goodJetIdx):\n", " for b1Hx, b2Hx in zip(tb1_e, tb2_e):\n", " if (b1Hx in goodJetIdx_e) & (b2Hx in goodJetIdx_e):\n", " N2 += 1" ] }, { "cell_type": "markdown", "id": "9ba3e441", "metadata": {}, "source": [ "#### Getting N3" ] }, { "cell_type": "code", "execution_count": 15, "id": "5bf115bc", "metadata": {}, "outputs": [], "source": [ "# mask\n", "# mask\n", "h1_mask = np.array(testfile['TARGETS']['h1']['mask'])\n", "h2_mask = np.array(testfile['TARGETS']['h2']['mask'])\n", "h3_mask = np.array(testfile['TARGETS']['h3']['mask'])\n", "h_masks = np.concatenate((h1_mask.reshape(-1, 1), h2_mask.reshape(-1, 1), h3_mask.reshape(-1, 1)), axis=1)\n", "\n", "bh1_mask = np.array(testfile['TARGETS']['bh1']['mask'])\n", "bh2_mask = np.array(testfile['TARGETS']['bh2']['mask'])\n", "bh3_mask = np.array(testfile['TARGETS']['bh3']['mask'])\n", "\n", "bh_masks = np.concatenate((bh1_mask.reshape(-1, 1), bh2_mask.reshape(-1, 1), bh3_mask.reshape(-1, 1)), axis=1)" ] }, { "cell_type": "code", "execution_count": 16, "id": "96b3648c", "metadata": {}, "outputs": [], "source": [ "N3=np.sum(h_masks&bh_masks)" ] }, { "cell_type": "markdown", "id": "d0142eed", "metadata": {}, "source": [ "#### Summary" ] }, { "cell_type": "code", "execution_count": 17, "id": "715c5977", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.017112822393350446" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# this percents of resolved Higgs should not be removed\n", "(N1-N2-N3)/N1" ] }, { "cell_type": "code", "execution_count": 18, "id": "2bf46477", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7944" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# this number of resolved Higgs removed for overlapping with a boosted Higgs\n", "N1-N2" ] }, { "cell_type": "code", "execution_count": 19, "id": "460fd60a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6544" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# this number of resolved Higgs can also be reconstructed as a boosted Higgs\n", "N3" ] }, { "cell_type": "code", "execution_count": 20, "id": "d855337e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.07999022124434665" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "N3/N1" ] }, { "cell_type": "markdown", "id": "6e88329b", "metadata": {}, "source": [ "#### Confusion Matrix" ] }, { "cell_type": "code", "execution_count": 21, "id": "62037f53", "metadata": {}, "outputs": [], "source": [ "N_overlap_bicat = 0\n", "N_overlap_resolve = 0\n", "N_unoverlap_bicat = 0\n", "N_unoverlap_resolve = 0\n", "\n", "bAndr = h_masks&bh_masks\n", "for tb1_e, tb2_e, goodJetIdx_e, h_masks_e, bh_masks_e, in zip(b1_ts, b2_ts, goodJetIdx, h_masks, bh_masks):\n", " for b1Hx, b2Hx, resolve, boost in zip(tb1_e, tb2_e, h_masks_e, bh_masks_e):\n", " \n", " if (b1Hx in goodJetIdx_e) & (b2Hx in goodJetIdx_e):\n", " overlap = False\n", " else:\n", " overlap=True\n", " \n", " if resolve & boost:\n", " bicat = True\n", " else:\n", " bicat = False\n", " \n", " if overlap & bicat:\n", " N_overlap_bicat += 1\n", " elif ~overlap & bicat:\n", " N_unoverlap_bicat += 1\n", " elif overlap & resolve:\n", " N_overlap_resolve += 1\n", " elif ~overlap & resolve:\n", " N_unoverlap_resolve += 1" ] }, { "cell_type": "code", "execution_count": 22, "id": "cb310983", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Resolved+Boosted | Resolved only\n", "Removed: 6544 1400\n", "Not Removed: 0 73866\n" ] } ], "source": [ "print(\" \", \"Resolved+Boosted |\", \"Resolved only\")\n", "print(\"Removed: \", N_overlap_bicat, \" \",N_overlap_resolve)\n", "print(\"Not Removed: \", N_unoverlap_bicat, \" \", N_unoverlap_resolve)" ] }, { "cell_type": "code", "execution_count": null, "id": "1a2454bc", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.10" } }, "nbformat": 4, "nbformat_minor": 5 }