From 666c073f6c9cf181f1f4c321580f1f684b5581fc Mon Sep 17 00:00:00 2001 From: Farouk Date: Thu, 27 Nov 2025 10:18:43 +0100 Subject: [PATCH] updating Figure 4 post-CWR --- .../cms/cmssw-validation-particles.ipynb | 2745 ++++++++++------- 1 file changed, 1582 insertions(+), 1163 deletions(-) diff --git a/notebooks/cms/cmssw-validation-particles.ipynb b/notebooks/cms/cmssw-validation-particles.ipynb index 2bfe40669..448c8a4c3 100644 --- a/notebooks/cms/cmssw-validation-particles.ipynb +++ b/notebooks/cms/cmssw-validation-particles.ipynb @@ -30,6 +30,7 @@ "import boost_histogram as bh\n", "import numba\n", "import mplhep\n", + "import mplhep as hep\n", "import sklearn\n", "import sklearn.metrics\n", "import matplotlib\n", @@ -52,6 +53,26 @@ "from plot_utils import SAMPLE_LABEL_CMS, pid_to_text, EXPERIMENT_LABELS" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "55aad96d", + "metadata": {}, + "outputs": [], + "source": [ + "pid_to_text[22] = \"photon\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "77597084", + "metadata": {}, + "outputs": [], + "source": [ + "EVALUATION_DATASET_NAMES[\"cms_pf_qcd_nopu\"] = \"QCD multijet, no PU\"" + ] + }, { "cell_type": "code", "execution_count": null, @@ -59,12 +80,16 @@ "metadata": {}, "outputs": [], "source": [ - "def sample_label(ax, sample, additional_text=\"\", x=0.03, y=0.97, fontsize=20):\n", + "def sample_label(ax, sample, additional_text=\"\", x=0.03, y=0.97, fontsize=20, color=\"black\"):\n", " text = EVALUATION_DATASET_NAMES[sample]\n", - " plt.text(x, y, text + additional_text, ha=\"left\", va=\"top\", transform=ax.transAxes, fontsize=fontsize)\n", + " plt.text(x, y, text + additional_text, ha=\"left\", va=\"top\", transform=ax.transAxes, fontsize=fontsize, color=color)\n", + " \n", + " \n", + "def cms_label(ax): \n", + " hep.cms.lumitext(\"Run 3 (13.6 TeV)\", ax=ax, fontsize=30)\n", + " hep.cms.text(\"Simulation\", ax=ax, fontsize=25)\n", "\n", - "def cms_label(ax):\n", - " return experiment_label(ax, experiment=\"CMS\", tag1=\"Simulation Preliminary\", tag2=\"Run 3 (13.6 TeV)\", x1=0.13)" + "# return experiment_label(ax, experiment=\"CMS\", tag1=\"Simulation\", tag2=\"Run 3 (13.6 TeV)\", x1=0.13)" ] }, { @@ -255,24 +280,6 @@ " return RelIso" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "54348d3e", - "metadata": {}, - "outputs": [], - "source": [ - "dataset = \"cms\"\n", - "\n", - "save_as = {\n", - " \"cms_pf_qcd_nopu\": \"QCD_noPU_13p6\",\n", - " \"cms_pf_ttbar_nopu\": \"TTbar_noPU_13p6\",\n", - "}\n", - "\n", - "sample = \"cms_pf_ttbar_nopu\"\n", - "# sample = \"cms_pf_qcd_nopu\"" - ] - }, { "cell_type": "code", "execution_count": null, @@ -299,19 +306,21 @@ " # Compute masks\n", " status = gen[\"status\"]\n", " pdgId = gen[\"pdgId\"]\n", - " mask_final = (status == 1)\n", - " mask_nu = (abs(pdgId) == vELE_PDGID) | \\\n", + " \n", + " msk_status = (status == 1)\n", + " msk_nu = (abs(pdgId) == vELE_PDGID) | \\\n", " (abs(pdgId) == vMU_PDGID) | \\\n", " (abs(pdgId) == vTAU_PDGID)\n", - " mask = mask_final & ~mask_nu\n", + " \n", + " msk = msk_status & ~msk_nu\n", "\n", " # Slice and zip\n", " sliced = {\n", - " \"GenCands_pt\": gen.pt[mask],\n", - " \"GenCands_eta\": gen.eta[mask],\n", - " \"GenCands_phi\": gen.phi[mask],\n", - " \"GenCands_energy\": gen.energy[mask],\n", - " \"GenCands_pdgId\": gen.pdgId[mask],\n", + " \"GenCands_pt\": gen.pt[msk],\n", + " \"GenCands_eta\": gen.eta[msk],\n", + " \"GenCands_phi\": gen.phi[msk],\n", + " \"GenCands_energy\": gen.energy[msk],\n", + " \"GenCands_pdgId\": gen.pdgId[msk],\n", " \"PFCands_pt\": pf.pt,\n", " \"PFCands_eta\": pf.eta,\n", " \"PFCands_phi\": pf.phi,\n", @@ -322,29 +331,29 @@ " remap_pid_gen(sliced)\n", "\n", "# get rel iso\n", - "# gen_vecs = ak.zip(\n", - "# {\n", - "# \"pt\": sliced[\"GenCands_pt\"],\n", - "# \"eta\": sliced[\"GenCands_eta\"],\n", - "# \"phi\": sliced[\"GenCands_phi\"],\n", - "# \"energy\": sliced[\"GenCands_energy\"],\n", - "# \"pid\": sliced[\"GenCands_pid\"],\n", - "# },\n", - "# with_name=\"Momentum4D\",\n", - "# )\n", + " gen_vecs = ak.zip(\n", + " {\n", + " \"pt\": sliced[\"GenCands_pt\"],\n", + " \"eta\": sliced[\"GenCands_eta\"],\n", + " \"phi\": sliced[\"GenCands_phi\"],\n", + " \"energy\": sliced[\"GenCands_energy\"],\n", + " \"pid\": sliced[\"GenCands_pid\"],\n", + " },\n", + " with_name=\"Momentum4D\",\n", + " )\n", "# sliced[\"GenCands_reliso03\"] = compute_relative_isolation(gen_vecs, dr=0.3)\n", "# sliced[\"GenCands_reliso05\"] = compute_relative_isolation(gen_vecs, dr=0.5)\n", "\n", - "# pf_vecs = ak.zip(\n", - "# {\n", - "# \"pt\": sliced[\"PFCands_pt\"],\n", - "# \"eta\": sliced[\"PFCands_eta\"], \n", - "# \"phi\": sliced[\"PFCands_phi\"], \n", - "# \"energy\": sliced[\"PFCands_energy\"], \n", - "# \"pid\": sliced[\"PFCands_pid\"],\n", - "# },\n", - "# with_name=\"Momentum4D\",\n", - "# ) \n", + " pf_vecs = ak.zip(\n", + " {\n", + " \"pt\": sliced[\"PFCands_pt\"],\n", + " \"eta\": sliced[\"PFCands_eta\"], \n", + " \"phi\": sliced[\"PFCands_phi\"], \n", + " \"energy\": sliced[\"PFCands_energy\"], \n", + " \"pid\": sliced[\"PFCands_pid\"],\n", + " },\n", + " with_name=\"Momentum4D\",\n", + " ) \n", "# sliced[\"PFCands_reliso03\"] = compute_relative_isolation(pf_vecs, dr=0.3)\n", "# sliced[\"PFCands_reliso05\"] = compute_relative_isolation(pf_vecs, dr=0.5)\n", " \n", @@ -369,6 +378,24 @@ " return ak.concatenate(results)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "54348d3e", + "metadata": {}, + "outputs": [], + "source": [ + "dataset = \"cms\"\n", + "\n", + "save_as = {\n", + " \"cms_pf_qcd_nopu\": \"QCD_noPU_13p6\",\n", + " \"cms_pf_ttbar_nopu\": \"TTbar_noPU_13p6\",\n", + "}\n", + "\n", + "sample = \"cms_pf_ttbar_nopu\"\n", + "# sample = \"cms_pf_qcd_nopu\"" + ] + }, { "cell_type": "code", "execution_count": null, @@ -376,7 +403,7 @@ "metadata": {}, "outputs": [], "source": [ - "! ls ../../data/CMSSW_15_0_5_mlpf_v2.5.0_p01_f8ae2f_test/cuda_False/" + "! ls ../../data/CMSSW_15_0_5_mlpf_v2.6.0pre1_p05_ccd0c7/cuda_False/" ] }, { @@ -386,8 +413,8 @@ "metadata": {}, "outputs": [], "source": [ - "path_mlpf = f\"../../data/CMSSW_15_0_5_mlpf_v2.5.0_p01_f8ae2f_test/cuda_False/{save_as[sample]}_mlpfpu/step3_MINI_*.pkl\"\n", - "path_pf = f\"../../data/CMSSW_15_0_5_mlpf_v2.5.0_p01_f8ae2f_test/cuda_False/{save_as[sample]}_pf/step3_MINI_*.pkl\"\n", + "path_mlpf = f\"../../data/CMSSW_15_0_5_mlpf_v2.6.0pre1_p05_ccd0c7/cuda_False/{save_as[sample]}_mlpfpu/step3_MINI_*.pkl\"\n", + "path_pf = f\"../../data/CMSSW_15_0_5_mlpf_v2.6.0pre1_p05_ccd0c7/cuda_False/{save_as[sample]}_pf/step3_MINI_*.pkl\"\n", "\n", "files_to_skip = []\n", "# mlpf\n", @@ -410,35 +437,35 @@ "metadata": {}, "outputs": [], "source": [ - "def remap_pid_mlpf(data):\n", - " \"\"\"\n", - " MLPF reconstructs...\n", - " \"\"\"\n", - " # define eta cut and pid values\n", - " eta_cut = 2.5\n", - " old_pid_to_new = {\n", - " 22: 1,\n", - " 130: 2,\n", - " 211: 2, \n", - " }\n", - "\n", - " # original arrays\n", - " eta = data[\"PFCands_eta\"]\n", - " pid = data[\"PFCands_pid\"]\n", - "\n", - " # start with the original pid array\n", - " new_pid = pid\n", - "\n", - " # apply remapping\n", - " for old_pid, new_pid_val in old_pid_to_new.items():\n", - " mask = (abs(eta) > eta_cut) & (pid == old_pid)\n", - " new_pid = ak.where(mask, new_pid_val, new_pid)\n", - "\n", - " # now new_pid is your remapped array\n", - " # you can either use it directly or attach it as a new field\n", - " data[\"PFCands_pid\"] = new_pid\n", + "# def remap_pid_mlpf(data):\n", + "# \"\"\"\n", + "# MLPF reconstructs...\n", + "# \"\"\"\n", + "# # define eta cut and pid values\n", + "# eta_cut = 2.4\n", + "# old_pid_to_new = {\n", + "# 22: 1,\n", + "# 130: 2,\n", + "# 211: 2, \n", + "# }\n", + "\n", + "# # original arrays\n", + "# eta = data[\"PFCands_eta\"]\n", + "# pid = data[\"PFCands_pid\"]\n", + "\n", + "# # start with the original pid array\n", + "# new_pid = pid\n", + "\n", + "# # apply remapping\n", + "# for old_pid, new_pid_val in old_pid_to_new.items():\n", + "# mask = (abs(eta) > eta_cut) & (pid == old_pid)\n", + "# new_pid = ak.where(mask, new_pid_val, new_pid)\n", + "\n", + "# # now new_pid is your remapped array\n", + "# # you can either use it directly or attach it as a new field\n", + "# data[\"PFCands_pid\"] = new_pid\n", " \n", - "remap_pid_mlpf(data_mlpf)" + "# remap_pid_mlpf(data_mlpf)" ] }, { @@ -457,9 +484,9 @@ "outputs": [], "source": [ "color_code = {\n", - " \"Gen\": \"tab:blue\",\n", - " \"PF\": \"tab:orange\",\n", - " \"MLPF\": \"tab:red\",\n", + " \"Gen\": \"#648DF4\",\n", + " \"PF\": \"#F3A041\",\n", + " \"MLPF\": \"#D23B3D\",\n", "}\n", "\n", "bins_pt = {\n", @@ -481,18 +508,18 @@ "\n", "bins_eta = {\n", " \"cms_pf_qcd_nopu\": {\n", - " 211: np.linspace(-2.5,2.5,41),\n", - " 130: np.linspace(-2.5,2.5,41),\n", - " 22: np.linspace(-2.5,2.5,41),\n", - " 11: np.linspace(-2.5,2.5,41),\n", - " 13: np.linspace(-2.5,2.5,41),\n", + " 211: np.linspace(-2.4,2.4,41),\n", + " 130: np.linspace(-2.4,2.4,41),\n", + " 22: np.linspace(-2.4,2.4,41),\n", + " 11: np.linspace(-2.4,2.4,41),\n", + " 13: np.linspace(-2.4,2.4,41),\n", " },\n", " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(-2.5,2.5,41),\n", - " 130: np.linspace(-2.5,2.5,41),\n", - " 22: np.linspace(-2.5,2.5,41),\n", - " 11: np.linspace(-2.5,2.5,41),\n", - " 13: np.linspace(-2.5,2.5,41),\n", + " 211: np.linspace(-2.4,2.4,41),\n", + " 130: np.linspace(-2.4,2.4,41),\n", + " 22: np.linspace(-2.4,2.4,41),\n", + " 11: np.linspace(-2.4,2.4,41),\n", + " 13: np.linspace(-2.4,2.4,41),\n", " },\n", "}\n", "\n", @@ -558,15 +585,15 @@ }, { "cell_type": "markdown", - "id": "0f9f69bd", + "id": "aa0832c2", "metadata": {}, "source": [ - "# Kinematic distributions" + "# Stacked hist" ] }, { "cell_type": "markdown", - "id": "00cef323", + "id": "06a4d170", "metadata": {}, "source": [ "## pT" @@ -575,631 +602,318 @@ { "cell_type": "code", "execution_count": null, - "id": "cb9b3b38", + "id": "f3104db6", "metadata": {}, "outputs": [], "source": [ - "ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " 2: 10**6,\n", - " 1: 10**6, \n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " 2: 10**6,\n", - " 1: 10**6, \n", - " }, \n", - "}\n", + "from matplotlib.lines import Line2D\n", "\n", - "rax_ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: (0,2),\n", - " 130: (0,2.2),\n", - " 22: (0,4),\n", - " 11: (0,8),\n", - " 13: (0,4),\n", - " 2: (0,4),\n", - " 1: (0,4), \n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: (0.5,1.5),\n", - " 130: (0,2),\n", - " 22: (0,2.5),\n", - " 11: (0,3),\n", - " 13: (0,2.4),\n", - " 2: (0,2.4),\n", - " 1: (0,2.4), \n", - " }, \n", + "ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": 10**12,\n", + " \"cms_pf_ttbar_nopu\": 5*10**8,\n", "}\n", "\n", - "bins_pt = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: np.linspace(1,500,41),\n", - " 130: np.linspace(1,200,41),\n", - " 22: np.linspace(1,200,41),\n", - " 11: np.linspace(1,50,21),\n", - " 13: np.linspace(1,50,21),\n", - " 2: np.linspace(1,50,21),\n", - " 1: np.linspace(1,50,21), \n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(1,100,31),\n", - " 130: np.linspace(1,40,31),\n", - " 22: np.linspace(1,40,31),\n", - " 11: np.linspace(1,100,31),\n", - " 13: np.linspace(1,100,7),\n", - " 2: np.linspace(1,100,7),\n", - " 1: np.linspace(1,100,7), \n", - " }, \n", - "}\n", + "bins = np.linspace(1,100,30)\n", + "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", "\n", - "eta_cut = 2.5\n", + "fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(12, 12), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + ")\n", "\n", + "gen_entries = []\n", + "pf_entries = []\n", + "mlpf_entries = []\n", "for pid in [\n", - "# 211,\n", + " 13,\n", + " 11,\n", " 130,\n", - "# 22,\n", - "# 11,\n", - "# 13,\n", - "# 2,\n", - "# 1,\n", + " 22, \n", + " 211,\n", "]:\n", - " fig, (ax, rax) = plt.subplots(\n", - " 2, 1, figsize=(10, 10), sharex=True,\n", - " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", - " )\n", "\n", - " # Extract flattened arrays\n", - " sel = (abs(data_pf[\"GenCands_pid\"]) == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + " # Gen\n", + " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1) & (abs(data_pf[\"GenCands_eta\"]) < 2.4)\n", " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_pt\"][sel]))\n", + " \n", + " total_yield = len(X_gen)\n", + " gen_entries.append((total_yield, X_gen, pid_to_text[pid], col_pid[pid]))\n", "\n", - " sel = (abs(data_pf[\"PFCands_pid\"]) == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) \n", + " # PF\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) < 2.4)\n", " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_pt\"][sel]))\n", "\n", - " sel = (abs(data_mlpf[\"PFCands_pid\"]) == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + " total_yield = len(X_pf)\n", + " pf_entries.append((total_yield, X_pf, pid_to_text[pid], col_pid[pid]))\n", + " \n", + " # MLPF\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.4)\n", " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_pt\"][sel]))\n", + " \n", + " total_yield = len(X_mlpf)\n", + " mlpf_entries.append((total_yield, X_mlpf, pid_to_text[pid], col_pid[pid]))\n", "\n", - " bins = bins_pt[sample][pid]\n", - " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", - " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", "\n", - " # Histograms\n", - " h_gen, _ = np.histogram(X_gen, bins=bins)\n", - " h_pf, _ = np.histogram(X_pf, bins=bins)\n", - " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", + "#### Gen legend\n", + "gen_data_list = [entry[1] for entry in gen_entries]\n", + "gen_labels = [f\"Gen {entry[2]}\" for entry in gen_entries]\n", + "gen_colors = [entry[3] for entry in gen_entries]\n", "\n", - " # Main plot\n", - " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "ax.hist(\n", + " gen_data_list,\n", + " bins=bins,\n", + " stacked=True,\n", + " histtype='stepfilled',\n", + " color=gen_colors,\n", + " label=gen_labels,\n", + " alpha=0.3,\n", + " edgecolor='k',\n", + " linewidth=0.8\n", + ")\n", "\n", - " ax.set_yscale(\"log\")\n", - " ax.set_ylabel(\"Counts\")\n", - " ax.set_ylim(1, ylim_dict[sample][pid])\n", - " ax.set_xlim(1)\n", - " ax.legend()\n", - " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", - " \n", - " plt.text(0.03, 0.82, pid_to_text[pid] + f\", $p_T>1$GeV\", transform=ax.transAxes, fontsize=25)\n", + "#### PF legend\n", + "pf_data_list = [entry[1] for entry in pf_entries]\n", + "pf_labels = [entry[2] for entry in pf_entries]\n", + "pf_colors = [entry[3] for entry in pf_entries]\n", "\n", - " # Ratio + error bars\n", - " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", - " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + "for i in range(len(pf_labels)):\n", + " if i == 0:\n", + " cum_sum = pf_data_list[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,pf_data_list[i]])\n", "\n", - " # Error bars: Poisson error propagation\n", - " err_pf = np.sqrt(h_pf) / h_gen\n", - " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", - " err_pf[h_gen == 0] = np.nan\n", - " err_mlpf[h_gen == 0] = np.nan\n", + " h_pf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + " ax.errorbar(\n", + " bin_centers, h_pf, yerr=np.sqrt(h_pf),\n", + "# linestyle=linestyle[\"PF\"],\n", + " fmt=marker_style[\"PF\"], markersize=7, markeredgewidth=2,\n", + " label=f\"PF {pf_labels[i]}\",\n", + " color=pf_colors[i],\n", + " )\n", "\n", - " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", - " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "#### MLPF legend\n", + "mlpf_data_list = [entry[1] for entry in mlpf_entries]\n", + "mlpf_labels = [entry[2] for entry in mlpf_entries]\n", + "mlpf_colors = [entry[3] for entry in mlpf_entries]\n", "\n", - " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", - " rax.set_ylim(rax_ylim_dict[sample][pid])\n", - " rax.set_ylabel(\"Ratio\")\n", - " \n", - " if sample==\"cms_pf_ttbar_nopu\":\n", - " if ( (pid==211) ):\n", - " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", - " else:\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + "for i in range(len(mlpf_labels)):\n", + " if i == 0:\n", + " cum_sum = mlpf_data_list[i]\n", " else:\n", - "# if ( (pid==22) ):\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates \n", - " \n", - " \n", - " rax.set_xlabel(\"$p_T$ (GeV)\")\n", - " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + " cum_sum = np.concatenate([cum_sum,mlpf_data_list[i]])\n", "\n", - " plt.savefig(f\"./plots/{save_as[sample]}/particle_pt_dist/{sample}_pt_{pid}.pdf\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fd8664ff", - "metadata": {}, - "outputs": [], - "source": [ - "ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " }, \n", - "}\n", - "\n", - "rax_ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: (0,2),\n", - " 130: (0,2.2),\n", - " 22: (0,4),\n", - " 11: (0,8),\n", - " 13: (0,4),\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: (0.5,1.5),\n", - " 130: (0,2),\n", - " 22: (0,2.5),\n", - " 11: (0,3),\n", - " 13: (0,2.4),\n", - " }, \n", - "}\n", - "\n", - "bins_pt = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: np.linspace(1,500,41),\n", - " 130: np.linspace(1,200,41),\n", - " 22: np.linspace(1,200,41),\n", - " 11: np.linspace(1,50,21),\n", - " 13: np.linspace(1,50,21),\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(1,100,31),\n", - " 130: np.linspace(1,40,31),\n", - " 22: np.linspace(1,40,31),\n", - " 11: np.linspace(1,100,31),\n", - " 13: np.linspace(1,100,100),\n", - " }, \n", - "}\n", - "\n", - "eta_cut = 2.5\n", - "\n", - "for pid in [\n", - "# 211,\n", - "# 130,\n", - "# 22,\n", - "# 11,\n", - " 13,\n", - "]:\n", - " fig, (ax, rax) = plt.subplots(\n", - " 2, 1, figsize=(10, 10), sharex=True,\n", - " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", - " )\n", - "\n", - " # Extract flattened arrays\n", - " sel = (abs(data_pf[\"GenCands_pid\"]) == pid) & (abs(data_pf[\"GenCands_eta\"]) < eta_cut) & (abs(data_pf[\"GenCands_pt\"]) > 10)\n", - " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_pt\"][sel]))\n", - "\n", - " sel = (abs(data_pf[\"PFCands_pid\"]) == pid) & (abs(data_pf[\"PFCands_eta\"]) < eta_cut) & (abs(data_pf[\"PFCands_pt\"]) > 10) \n", - " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_pt\"][sel]))\n", - "\n", - " sel = (abs(data_mlpf[\"PFCands_pid\"]) == pid) & (abs(data_mlpf[\"PFCands_eta\"]) < eta_cut) & (abs(data_mlpf[\"PFCands_pt\"]) > 10)\n", - " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_pt\"][sel]))\n", - "\n", - " bins = bins_pt[sample][pid]\n", - " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", - " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", - "\n", - " # Histograms\n", - " h_gen, _ = np.histogram(X_gen, bins=bins)\n", - " h_pf, _ = np.histogram(X_pf, bins=bins)\n", - " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", - "\n", - " # Main plot\n", - " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", - "\n", - " ax.set_yscale(\"log\")\n", - " ax.set_ylabel(\"Counts\")\n", - " ax.set_ylim(1, ylim_dict[sample][pid])\n", - " ax.set_xlim(1)\n", - " ax.legend()\n", - " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", - " \n", - " plt.text(0.03, 0.82, pid_to_text[pid] + f\", $p_T>1$GeV, $|\\eta| < {eta_cut}$\", transform=ax.transAxes, fontsize=25)\n", - "\n", - " # Ratio + error bars\n", - " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", - " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", - "\n", - " # Error bars: Poisson error propagation\n", - " err_pf = np.sqrt(h_pf) / h_gen\n", - " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", - " err_pf[h_gen == 0] = np.nan\n", - " err_mlpf[h_gen == 0] = np.nan\n", - "\n", - " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", - " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", - "\n", - " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", - " rax.set_ylim(rax_ylim_dict[sample][pid])\n", - " rax.set_ylabel(\"Ratio\")\n", + " h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", " \n", - " if sample==\"cms_pf_ttbar_nopu\":\n", - " if ( (pid==211) ):\n", - " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", - " else:\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", - " else:\n", - "# if ( (pid==22) ):\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates \n", - " \n", - " \n", - " rax.set_xlabel(\"$p_T$ (GeV)\")\n", - " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", - "\n", - " plt.savefig(f\"./plots/{save_as[sample]}/particle_pt_dist/{sample}_pt_{pid}.pdf\")" - ] - }, - { - "cell_type": "markdown", - "id": "23923816", - "metadata": {}, - "source": [ - "## eta" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "777d1c43", - "metadata": {}, - "outputs": [], - "source": [ - "bins_eta = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: np.linspace(-2.5,2.5,41),\n", - " 130: np.linspace(-2.5,2.5,41),\n", - " 22: np.linspace(-2.5,2.5,41),\n", - " 11: np.linspace(-2.5,2.5,41),\n", - " 13: np.linspace(-2.5,2.5,41),\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(-5,5,41),\n", - " 130: np.linspace(-5,5,41),\n", - " 22: np.linspace(-5,5,41),\n", - " 11: np.linspace(-5,5,41),\n", - " 13: np.linspace(-5,5,41),\n", - " 2: np.linspace(-5,5,31),\n", - " 1: np.linspace(-5,5,31), \n", - " },\n", - "}\n", - "\n", - "ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " 2: 10**6,\n", - " 1: 10**6, \n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " 2: 10**6,\n", - " 1: 10**6, \n", - " }, \n", - "}\n", + " ax.errorbar(\n", + " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", + "# linestyle=linestyle[\"PF\"],\n", + " fmt=marker_style[\"MLPF\"], markersize=6, markeredgewidth=1,\n", + " label=f\"MLPF {mlpf_labels[i]}\",\n", + " color=mlpf_colors[i],\n", + " )\n", "\n", - "rax_ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: (0.8,1.6),\n", - " 130: (0,3),\n", - " 22: (0,2),\n", - " 11: (0,2.5),\n", - " 13: (0,5),\n", - " 2: (0,5),\n", - " 1: (0,5), \n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: (0.4,1.6),\n", - " 130: (0,2),\n", - " 22: (0,1.7),\n", - " 11: (0,2),\n", - " 13: (0,2),\n", - " 2: (0,2),\n", - " 1: (0,2), \n", - " }, \n", - "}\n", "\n", - "for pid in [\n", - " 211,\n", - " 130,\n", - " 22,\n", - " 11,\n", - " 13,\n", - " 2,\n", - " 1,\n", - "]:\n", - " fig, (ax, rax) = plt.subplots(\n", - " 2, 1, figsize=(10, 10), sharex=True,\n", - " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + "############## add legend\n", + "labels_Gen, labels_PF, labels_MLPF = [], [], []\n", + "handles_Gen, handles_PF, handles_MLPF = [], [], []\n", + "for i, lab in enumerate(labels):\n", + " if lab.startswith(\"Gen\"):\n", + " labels_Gen += [lab.replace(\"Gen \", \"\")]\n", + " handles_Gen += [handles[i]]\n", + " if lab.startswith(\"PF\"):\n", + " labels_PF += [lab.replace(\"PF \", \"\")]\n", + " handles_PF += [handles[i]]\n", + " if lab.startswith(\"MLPF\"):\n", + " labels_MLPF += [lab.replace(\"MLPF \", \"\")]\n", + " handles_MLPF += [handles[i]]\n", + "\n", + "# Build marker-only legend handles for PF and MLPF\n", + "pf_marker_handles = [\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"PF\"],\n", + " linestyle=\"None\",\n", + " color=h[0].get_color(), # use same color as the plotted PF series\n", + " markersize=7,\n", + " markeredgewidth=2,\n", " )\n", + " for h in handles_PF\n", + "]\n", + "\n", + "mlpf_marker_handles = [\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"MLPF\"],\n", + " linestyle=\"None\",\n", + " color=h[0].get_color(), # same color as MLPF series\n", + " markersize=6,\n", + " markeredgewidth=1,\n", + " )\n", + " for h in handles_MLPF\n", + "] \n", "\n", - " # Extract flattened arrays\n", - " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", - " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_eta\"][sel]))\n", - "\n", - " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", - " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_eta\"][sel]))\n", + "if \"qcd\" in sample:\n", + " bbox_to_anchor_xy = (0.64, 0.92)\n", + "else:\n", + " bbox_to_anchor_xy = (0.64, 0.98)\n", + " \n", + " \n", + "leg1 = ax.legend(\n", + " handles=handles_Gen,\n", + " labels=labels_Gen,\n", + " title=\"Gen\",\n", + " ncol=1, # Two-column layout\n", + " handletextpad=0.5, # Space between marker and text\n", + " title_fontsize=26,\n", + " fontsize=22,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0], bbox_to_anchor_xy[1]),\n", + ")\n", + "ax.add_artist(leg1)\n", "\n", - " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", - " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_eta\"][sel]))\n", + "handles_PF.reverse()\n", + "labels_PF.reverse()\n", + "leg2 = ax.legend(\n", + "# handles=handles_PF,\n", + " handles=pf_marker_handles,\n", + " labels=labels_PF,\n", + " title=\"PF\",\n", + " ncol=1, # Two-column layout\n", + " handletextpad=0.5, # Space between marker and text\n", + " title_fontsize=26,\n", + " fontsize=22,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*1, bbox_to_anchor_xy[1]),\n", + ")\n", + "ax.add_artist(leg2)\n", "\n", - " bins = bins_eta[sample][pid]\n", - " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", - " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", + "handles_MLPF.reverse()\n", + "labels_MLPF.reverse()\n", + "leg3 = ax.legend(\n", + "# handles=handles_MLPF,\n", + " handles=mlpf_marker_handles,\n", + " labels=labels_MLPF,\n", + " title=\"MLPF\",\n", + " ncol=1, # Two-column layout\n", + " handletextpad=0.5, # Space between marker and text\n", + " title_fontsize=26, \n", + " fontsize=22, # Font size of labels\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*2, bbox_to_anchor_xy[1]), \n", + ")\n", + "ax.add_artist(leg3)\n", "\n", - " # Histograms\n", - " h_gen, _ = np.histogram(X_gen, bins=bins)\n", - " h_pf, _ = np.histogram(X_pf, bins=bins)\n", - " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", + "ax.set_yscale(\"log\")\n", + "ax.set_ylabel(\"Total number of particles / bin\", fontsize=35)\n", + "ax.set_ylim(2*10**2, ylim_dict[sample])\n", + "ax.set_xlim(1,100)\n", "\n", - " # Main plot\n", - " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "cms_label(ax)\n", + "sample_label(ax, sample, fontsize=35)\n", "\n", - " ax.set_yscale(\"log\")\n", - " ax.set_ylabel(\"Counts\")\n", - " ax.set_ylim(1, ylim_dict[sample][pid])\n", - " ax.set_xlim(-3,3)\n", - " ax.legend()\n", - " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", + "########################### Ratio + error bars\n", + "for i in range(len(gen_labels)):\n", + " if i == 0:\n", + " cum_sum = gen_data_list[i]\n", " else:\n", - " sample_label(ax, sample, fontsize=30)\n", - " \n", - " plt.text(0.03, 0.82, pid_to_text[pid] + \", $p_T>1$GeV\", transform=ax.transAxes, fontsize=25)\n", - "\n", - " # Ratio + error bars\n", - " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", - " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", - "\n", - " # Error bars: Poisson error propagation\n", - " err_pf = np.sqrt(h_pf) / h_gen\n", - " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", - " err_pf[h_gen == 0] = np.nan\n", - " err_mlpf[h_gen == 0] = np.nan\n", - "\n", - " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", - " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + " cum_sum = np.concatenate([cum_sum,gen_data_list[i]])\n", + "h_gen, _ = np.histogram(cum_sum, bins=bins)\n", "\n", - " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", - " rax.set_ylim(rax_ylim_dict[sample][pid])\n", - " rax.set_ylabel(\"Ratio\")\n", + "for i in range(len(pf_labels)):\n", + " if i == 0:\n", + " cum_sum = pf_data_list[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,pf_data_list[i]])\n", + "h_pf, _ = np.histogram(cum_sum, bins=bins)\n", " \n", - " if sample==\"cms_pf_qcd_nopu\":\n", - " if (pid==211):\n", - " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", - " else:\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", - "\n", - " if sample==\"cms_pf_ttbar_nopu\":\n", - " if ( (pid==211) ):\n", - " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates \n", - " else:\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", - " \n", - " rax.set_xlabel(r\"$\\eta$\")\n", - " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", - "\n", - " plt.xlim(-5, 5)\n", - " \n", - " plt.savefig(f\"./plots/{save_as[sample]}/particle_eta_dist/{sample}_eta_{pid}.pdf\")" - ] - }, - { - "cell_type": "markdown", - "id": "812e9548", - "metadata": {}, - "source": [ - "## phi" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c59fc6d9", - "metadata": {}, - "outputs": [], - "source": [ - "ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: 10**10,\n", - " 130: 10**8,\n", - " 22: 10**10,\n", - " 11: 10**7,\n", - " 13: 10**6,\n", - " }, \n", - "}\n", - "\n", - "rax_ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: (0.5,2),\n", - " 130: (0.5,1.5),\n", - " 22: (0,2),\n", - " 11: (0,2),\n", - " 13: (0,4),\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: (0.5,1.6),\n", - " 130: (0.3,1.8),\n", - " 22: (0,2),\n", - " 11: (0,2),\n", - " 13: (0.5,1.6),\n", - " }, \n", - "}\n", - "\n", - "for pid in [\n", - " 211,\n", - " 130,\n", - " 22,\n", - " 11,\n", - " 13,\n", - "]:\n", - " fig, (ax, rax) = plt.subplots(\n", - " 2, 1, figsize=(10, 10), sharex=True,\n", - " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", - " )\n", - "\n", - " # Extract flattened arrays\n", - " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_eta\"]) < 2.5) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", - " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_phi\"][sel]))\n", - "\n", - " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_eta\"]) < 2.5) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", - " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_phi\"][sel]))\n", - "\n", - " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.5) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", - " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_phi\"][sel]))\n", - "\n", - " bins = bins_phi[sample][pid]\n", - " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", - " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", - "\n", - " # Histograms\n", - " h_gen, _ = np.histogram(X_gen, bins=bins)\n", - " h_pf, _ = np.histogram(X_pf, bins=bins)\n", - " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", - "\n", - " # Main plot\n", - " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", - " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", - "\n", - " ax.set_yscale(\"log\")\n", - " ax.set_ylabel(\"Counts\")\n", - " ax.set_ylim(1, ylim_dict[sample][pid])\n", - " ax.set_xlim(-3,3)\n", - " ax.legend()\n", - " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", + "for i in range(len(mlpf_labels)):\n", + " if i == 0:\n", + " cum_sum = mlpf_data_list[i]\n", " else:\n", - " sample_label(ax, sample, fontsize=30)\n", - " \n", - " plt.text(0.03, 0.82, pid_to_text[pid] + \", $p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=25)\n", - "\n", - " # Ratio + error bars\n", - " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", - " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + " cum_sum = np.concatenate([cum_sum,mlpf_data_list[i]])\n", + "h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + "ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + "ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", "\n", - " # Error bars: Poisson error propagation\n", - " err_pf = np.sqrt(h_pf) / h_gen\n", - " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", - " err_pf[h_gen == 0] = np.nan\n", - " err_mlpf[h_gen == 0] = np.nan\n", + "# Error bars: Poisson error propagation\n", + "err_pf = np.sqrt(h_pf) / h_gen\n", + "err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + "err_pf[h_gen == 0] = np.nan\n", + "err_mlpf[h_gen == 0] = np.nan\n", "\n", - " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", - " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=7, markeredgewidth=2, marker=marker_style[\"PF\"], label=\"PF\", color=color_code[\"PF\"])\n", + "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=6, markeredgewidth=1, marker=marker_style[\"MLPF\"], label=\"MLPF\", color=color_code[\"MLPF\"])\n", "\n", - " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", - " rax.set_ylim(rax_ylim_dict[sample][pid])\n", - " rax.set_ylabel(\"Ratio\")\n", + "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + "rax.set_ylabel(\"Reco / Gen\")\n", "\n", - " if sample==\"cms_pf_qcd_nopu\":\n", - " if ( (pid==211) | (pid==130) ):\n", - " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates \n", - " else:\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + "rax.set_ylim(0, 2)\n", + "# rax.legend(fontsize=25, ncol=2, loc=\"upper center\")\n", + "\n", + "# --- Build custom legend handles (marker only, no error bars) ---\n", + "ratio_legend_handles = [\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"PF\"],\n", + " linestyle=\"None\",\n", + " markersize=7,\n", + " markeredgewidth=2,\n", + " markeredgecolor=color_code[\"PF\"],\n", + " markerfacecolor=color_code[\"PF\"],\n", + " label=\"PF\",\n", + " ),\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"MLPF\"],\n", + " linestyle=\"None\",\n", + " markersize=6,\n", + " markeredgewidth=1,\n", + " markeredgecolor=color_code[\"MLPF\"], \n", + " markerfacecolor=color_code[\"MLPF\"],\n", + " label=\"MLPF\",\n", + " ),\n", + "]\n", + "\n", + "rax.legend(\n", + " handles=ratio_legend_handles,\n", + " fontsize=25,\n", + " ncol=2,\n", + " loc=\"upper center\"\n", + ")\n", "\n", - " if sample==\"cms_pf_ttbar_nopu\":\n", - " if ( (pid==211) | (pid==130) | (pid==13) ):\n", - " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", - " else:\n", - " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates \n", + "rax.set_xlabel(\"Particle $p_T$ (GeV)\", fontsize=35)\n", + "rax.yaxis.set_label_coords(-0.03, 0.9) # move left/right and up/down\n", "\n", - " rax.set_xlabel(r\"$\\phi$\")\n", - " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + "plt.text(0.03, 0.82, r\"$p_T>1~$GeV, $|\\eta| < 2.4$\", transform=ax.transAxes, fontsize=26)\n", "\n", - " plt.savefig(f\"./plots/{save_as[sample]}/particle_phi_dist/{sample}_phi_{pid}.pdf\")" - ] - }, - { - "cell_type": "markdown", - "id": "aa0832c2", - "metadata": {}, - "source": [ - "# Stacked hist" + "plt.tight_layout()\n", + " \n", + "plt.savefig(f\"./plots/{save_as[sample]}/particle_pt_dist/{sample}_pt_stacked.pdf\")" ] }, { "cell_type": "markdown", - "id": "06a4d170", + "id": "13aad0e2", "metadata": {}, "source": [ - "## pT" + "## eta" ] }, { "cell_type": "code", "execution_count": null, - "id": "f3104db6", - "metadata": {}, + "id": "5e5f398c", + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ "ylim_dict = {\n", " \"cms_pf_qcd_nopu\": 10**12,\n", - " \"cms_pf_ttbar_nopu\": 10**9, \n", - "}\n", - "bbox_to_anchor_xy = {\n", - " \"cms_pf_qcd_nopu\": (0.55, 0.92),\n", - " \"cms_pf_ttbar_nopu\": (0.57, 0.98),\n", + " \"cms_pf_ttbar_nopu\": 5*10**8,\n", + "# \"cms_pf_ttbar_nopu\": 5*10**7, \n", "}\n", "\n", - "bins = np.linspace(1,100,30)\n", + "bins = np.linspace(-2.4,2.4,30)\n", "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", "\n", "fig, (ax, rax) = plt.subplots(\n", @@ -1219,22 +933,22 @@ "]:\n", "\n", " # Gen\n", - " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1) & (abs(data_pf[\"GenCands_eta\"]) < 2.5)\n", - " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_pt\"][sel]))\n", + " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1) & (abs(data_pf[\"GenCands_eta\"]) < 2.4)\n", + " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_eta\"][sel]))\n", " \n", " total_yield = len(X_gen)\n", " gen_entries.append((total_yield, X_gen, pid_to_text[pid], col_pid[pid]))\n", "\n", " # PF\n", - " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) < 2.5)\n", - " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_pt\"][sel]))\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) < 2.4)\n", + " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_eta\"][sel]))\n", "\n", " total_yield = len(X_pf)\n", " pf_entries.append((total_yield, X_pf, pid_to_text[pid], col_pid[pid]))\n", " \n", " # MLPF\n", - " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.5)\n", - " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_pt\"][sel]))\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.4)\n", + " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_eta\"][sel]))\n", " \n", " total_yield = len(X_mlpf)\n", " mlpf_entries.append((total_yield, X_mlpf, pid_to_text[pid], col_pid[pid]))\n", @@ -1273,7 +987,7 @@ " ax.errorbar(\n", " bin_centers, h_pf, yerr=np.sqrt(h_pf),\n", "# linestyle=linestyle[\"PF\"],\n", - " fmt=marker_style[\"PF\"], markersize=4, markeredgewidth=1,\n", + " fmt=marker_style[\"PF\"], markersize=7, markeredgewidth=2,\n", " label=f\"PF {pf_labels[i]}\",\n", " color=pf_colors[i],\n", " )\n", @@ -1294,15 +1008,13 @@ " ax.errorbar(\n", " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", "# linestyle=linestyle[\"PF\"],\n", - " fmt=marker_style[\"MLPF\"], markersize=4, markeredgewidth=1,\n", + " fmt=marker_style[\"MLPF\"], markersize=6, markeredgewidth=1,\n", " label=f\"MLPF {mlpf_labels[i]}\",\n", " color=mlpf_colors[i],\n", " )\n", "\n", "\n", "############## add legend\n", - "handles, labels = ax.get_legend_handles_labels()\n", - "\n", "labels_Gen, labels_PF, labels_MLPF = [], [], []\n", "handles_Gen, handles_PF, handles_MLPF = [], [], []\n", "for i, lab in enumerate(labels):\n", @@ -1314,17 +1026,48 @@ " handles_PF += [handles[i]]\n", " if lab.startswith(\"MLPF\"):\n", " labels_MLPF += [lab.replace(\"MLPF \", \"\")]\n", - " handles_MLPF += [handles[i]] \n", + " handles_MLPF += [handles[i]]\n", + "\n", + "# Build marker-only legend handles for PF and MLPF\n", + "pf_marker_handles = [\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"PF\"],\n", + " linestyle=\"None\",\n", + " color=h[0].get_color(), # use same color as the plotted PF series\n", + " markersize=7,\n", + " markeredgewidth=2,\n", + " )\n", + " for h in handles_PF\n", + "]\n", + "\n", + "mlpf_marker_handles = [\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"MLPF\"],\n", + " linestyle=\"None\",\n", + " color=h[0].get_color(), # same color as MLPF series\n", + " markersize=6,\n", + " markeredgewidth=1,\n", + " )\n", + " for h in handles_MLPF\n", + "] \n", "\n", + "if \"qcd\" in sample:\n", + " bbox_to_anchor_xy = (0.64, 0.92)\n", + "else:\n", + " bbox_to_anchor_xy = (0.64, 0.98)\n", + " \n", + " \n", "leg1 = ax.legend(\n", " handles=handles_Gen,\n", " labels=labels_Gen,\n", " title=\"Gen\",\n", " ncol=1, # Two-column layout\n", " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25,\n", - " fontsize=20,\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0], bbox_to_anchor_xy[sample][1]),\n", + " title_fontsize=26,\n", + " fontsize=22,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0], bbox_to_anchor_xy[1]),\n", ")\n", "ax.add_artist(leg1)\n", "\n", @@ -1336,9 +1079,9 @@ " title=\"PF\",\n", " ncol=1, # Two-column layout\n", " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25,\n", - " fontsize=20,\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0]+0.22*1, bbox_to_anchor_xy[sample][1]),\n", + " title_fontsize=26,\n", + " fontsize=22,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*1, bbox_to_anchor_xy[1]),\n", ")\n", "ax.add_artist(leg2)\n", "\n", @@ -1350,22 +1093,20 @@ " title=\"MLPF\",\n", " ncol=1, # Two-column layout\n", " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25, \n", - " fontsize=20, # Font size of labels\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0]+0.22*2, bbox_to_anchor_xy[sample][1]), \n", + " title_fontsize=26, \n", + " fontsize=22, # Font size of labels\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*2, bbox_to_anchor_xy[1]), \n", ")\n", "ax.add_artist(leg3)\n", "\n", "ax.set_yscale(\"log\")\n", - "ax.set_ylabel(\"Total number of particles / bin\")\n", - "ax.set_ylim(1, ylim_dict[sample])\n", - "ax.set_xlim(1,100)\n", + "ax.set_ylabel(\"Total number of particles / bin\", fontsize=35)\n", + "# ax.set_ylim(1, ylim_dict[sample]) # TODO\n", + "ax.set_ylim(900, ylim_dict[sample])\n", + "ax.set_xlim(-2.4,2.4)\n", "\n", "cms_label(ax)\n", - "if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=30)\n", - "else:\n", - " sample_label(ax, sample, fontsize=30)\n", + "sample_label(ax, sample, fontsize=35)\n", "\n", "########################### Ratio + error bars\n", "for i in range(len(gen_labels)):\n", @@ -1398,33 +1139,54 @@ "err_pf[h_gen == 0] = np.nan\n", "err_mlpf[h_gen == 0] = np.nan\n", "\n", - "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"PF\"], label=\"PF\", color=color_code[\"PF\"])\n", - "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"MLPF\"], label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=7, markeredgewidth=2, marker=marker_style[\"PF\"], label=\"PF\", color=color_code[\"PF\"])\n", + "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=6, markeredgewidth=1, marker=marker_style[\"MLPF\"], label=\"MLPF\", color=color_code[\"MLPF\"])\n", "\n", "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", "rax.set_ylabel(\"Reco / Gen\")\n", "\n", - "if sample == \"cms_pf_qcd_nopu\":\n", - " rax.set_ylim(0, 2.2)\n", - " rax.legend(fontsize=23, ncol=2, loc=\"upper center\")\n", - "else:\n", - " rax.set_ylim(0,2)\n", - " rax.legend(fontsize=23, ncol=2, loc=\"upper center\")\n", + "rax.set_ylim(0, 2)\n", + "# rax.legend(fontsize=25, ncol=2, loc=\"upper center\")\n", + "\n", + "# --- Build custom legend handles (marker only, no error bars) ---\n", + "ratio_legend_handles = [\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"PF\"],\n", + " linestyle=\"None\",\n", + " markersize=7,\n", + " markeredgewidth=2,\n", + " markeredgecolor=color_code[\"PF\"],\n", + " markerfacecolor=color_code[\"PF\"],\n", + " label=\"PF\",\n", + " ),\n", + " Line2D(\n", + " [0], [0],\n", + " marker=marker_style[\"MLPF\"],\n", + " linestyle=\"None\",\n", + " markersize=6,\n", + " markeredgewidth=1,\n", + " markeredgecolor=color_code[\"MLPF\"], \n", + " markerfacecolor=color_code[\"MLPF\"],\n", + " label=\"MLPF\",\n", + " ),\n", + "]\n", + "\n", + "rax.legend(\n", + " handles=ratio_legend_handles,\n", + " fontsize=25,\n", + " ncol=2,\n", + " loc=\"upper center\"\n", + ")\n", "\n", - "rax.set_xlabel(\"Particle $p_T$ (GeV)\")\n", + "rax.set_xlabel(\"Particle $\\eta$\", fontsize=35)\n", "rax.yaxis.set_label_coords(-0.03, 0.9) # move left/right and up/down\n", "\n", - "plt.text(0.03, 0.82, r\"$p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=22)\n", - " \n", - "plt.savefig(f\"./plots/{save_as[sample]}/particle_pt_dist/{sample}_pt_stacked.pdf\")" - ] - }, - { - "cell_type": "markdown", - "id": "13aad0e2", - "metadata": {}, - "source": [ - "## eta" + "plt.text(0.03, 0.82, r\"$p_T>1~$GeV, $|\\eta| < 2.4$\", transform=ax.transAxes, fontsize=26)\n", + "\n", + "plt.tight_layout()\n", + "\n", + "plt.savefig(f\"./plots/{save_as[sample]}/particle_eta_dist/{sample}_eta_stacked.pdf\")" ] }, { @@ -1444,11 +1206,7 @@ "\n", "ylim_dict = {\n", " \"cms_pf_qcd_nopu\": 10**16,\n", - " \"cms_pf_ttbar_nopu\": 10**14, \n", - "}\n", - "bbox_to_anchor_xy = {\n", - " \"cms_pf_qcd_nopu\": (0.55, 0.92),\n", - " \"cms_pf_ttbar_nopu\": (0.30, 0.55),\n", + " \"cms_pf_ttbar_nopu\": 10**14,\n", "}\n", "\n", "# bins = np.linspace(-5,5,50)\n", @@ -1503,14 +1261,14 @@ " \n", " ########### central\n", " # PF\n", - " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) < 2.5)\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) < 2.4)\n", " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_eta\"][sel]))\n", "\n", " total_yield = len(X_pf)\n", " pf_entries_central.append((total_yield, X_pf, pid_to_text[pid], col_pid[pid]))\n", "\n", " # MLPF\n", - " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.5)\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.4)\n", " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_eta\"][sel]))\n", "\n", " total_yield = len(X_mlpf)\n", @@ -1518,14 +1276,14 @@ " \n", " ########### forward\n", " # PF\n", - " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) > 2.5)\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) & (abs(data_pf[\"PFCands_eta\"]) > 2.4)\n", " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_eta\"][sel]))\n", "\n", " total_yield = len(X_pf)\n", " pf_entries_forward.append((total_yield, X_pf, pid_to_text[pid], col_pid[pid]))\n", "\n", " # MLPF\n", - " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) > 2.5)\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1) & (abs(data_mlpf[\"PFCands_eta\"]) > 2.4)\n", " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_eta\"][sel]))\n", "\n", " total_yield = len(X_mlpf)\n", @@ -1596,72 +1354,312 @@ "pf_labels_central = [entry[2] for entry in pf_entries_central]\n", "pf_colors_central = [entry[3] for entry in pf_entries_central]\n", "\n", - "for i in range(len(pf_labels_central)):\n", + "for i in range(len(pf_labels_central)):\n", + " \n", + " if i == 0:\n", + " cum_sum = pf_data_list_central[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,pf_data_list_central[i]])\n", + "\n", + " h_pf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + " ax.errorbar(\n", + " bin_centers, h_pf, yerr=np.sqrt(h_pf),\n", + " fmt=marker_style[\"PF\"], markersize=4, markeredgewidth=1, \n", + " label=f\"PF {pf_labels_central[i]}\",\n", + " color=pf_colors_central[i],\n", + " )\n", + " \n", + "#### MLPF legend\n", + "mlpf_data_list_forward = [entry[1] for entry in mlpf_entries_forward]\n", + "mlpf_labels_forward = [entry[2] for entry in mlpf_entries_forward]\n", + "mlpf_colors_forward = [entry[3] for entry in mlpf_entries_forward]\n", + "\n", + "# mlpf_data_list_forward.reverse()\n", + "# mlpf_labels_forward.reverse()\n", + "# mlpf_colors_forward.reverse()\n", + "\n", + "tmp = [mlpf_data_list_forward[i] for i in forward_order]\n", + "mlpf_data_list_forward = tmp\n", + "\n", + "tmp = [mlpf_labels_forward[i] for i in forward_order]\n", + "mlpf_labels_forward = tmp\n", + "\n", + "tmp = [mlpf_colors_forward[i] for i in forward_order]\n", + "mlpf_colors_forward = tmp\n", + "\n", + "for i in range(len(mlpf_labels_forward)):\n", + " if i == 0:\n", + " cum_sum = mlpf_data_list_forward[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,mlpf_data_list_forward[i]])\n", + "\n", + " h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + " ax.errorbar(\n", + " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", + " fmt=marker_style[\"MLPF\"], markersize=4, markeredgewidth=1,\n", + "# label=f\"MLPF {mlpf_labels_forward[i]}\",\n", + " color=mlpf_colors_forward[i],\n", + " )\n", + "\n", + "mlpf_data_list_central = [entry[1] for entry in mlpf_entries_central]\n", + "mlpf_labels_central = [entry[2] for entry in mlpf_entries_central]\n", + "mlpf_colors_central = [entry[3] for entry in mlpf_entries_central]\n", + "\n", + "for i in range(len(mlpf_labels_central)):\n", + " if i == 0:\n", + " cum_sum = mlpf_data_list_central[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,mlpf_data_list_central[i]])\n", + "\n", + " h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + " ax.errorbar(\n", + " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", + " fmt=marker_style[\"MLPF\"], markersize=4, markeredgewidth=1,\n", + " label=f\"MLPF {mlpf_labels_central[i]}\",\n", + " color=mlpf_colors_central[i],\n", + " )\n", + "\n", + "############## add legend\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "\n", + "labels_Gen, labels_PF, labels_MLPF = [], [], []\n", + "handles_Gen, handles_PF, handles_MLPF = [], [], []\n", + "for i, lab in enumerate(labels):\n", + " if lab.startswith(\"Gen\"):\n", + " labels_Gen += [lab.replace(\"Gen \", \"\")]\n", + " handles_Gen += [handles[i]]\n", + " if lab.startswith(\"PF\"):\n", + " labels_PF += [lab.replace(\"PF \", \"\")]\n", + " handles_PF += [handles[i]]\n", + " if lab.startswith(\"MLPF\"):\n", + " labels_MLPF += [lab.replace(\"MLPF \", \"\")]\n", + " handles_MLPF += [handles[i]] \n", + "\n", + "if \"qcd\" in sample:\n", + " bbox_to_anchor_xy = (0.64, 0.92)\n", + "else:\n", + " bbox_to_anchor_xy = (0.64, 0.98)\n", + " \n", + "leg1 = ax.legend(\n", + " handles=handles_Gen,\n", + " labels=labels_Gen,\n", + " title=\"Gen\",\n", + " ncol=1, # Two-column layout\n", + " handletextpad=0.5, # Space between marker and text\n", + " title_fontsize=25,\n", + " fontsize=20,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0], bbox_to_anchor_xy[1]),\n", + ")\n", + "ax.add_artist(leg1)\n", + "\n", + "handles_PF.reverse()\n", + "labels_PF.reverse()\n", + "leg2 = ax.legend(\n", + " handles=handles_PF,\n", + " labels=labels_PF,\n", + " title=\"PF\",\n", + " ncol=1, # Two-column layout\n", + " handletextpad=0.5, # Space between marker and text\n", + " title_fontsize=25,\n", + " fontsize=20,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*1, bbox_to_anchor_xy[1]),\n", + ")\n", + "ax.add_artist(leg2)\n", + "\n", + "handles_MLPF.reverse()\n", + "labels_MLPF.reverse()\n", + "leg3 = ax.legend(\n", + " handles=handles_MLPF,\n", + " labels=labels_MLPF,\n", + " title=\"MLPF\",\n", + " ncol=1, # Two-column layout\n", + " handletextpad=0.5, # Space between marker and text\n", + " title_fontsize=25, \n", + " fontsize=20, # Font size of labels\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*2, bbox_to_anchor_xy[1]), \n", + ")\n", + "ax.add_artist(leg3)\n", + "\n", + "ax.set_yscale(\"log\")\n", + "ax.set_ylabel(\"Total number of particles / bin\", fontsize=35)\n", + "ax.set_ylim(1, ylim_dict[sample])\n", + "ax.set_xlim(-5,5)\n", + "\n", + "cms_label(ax)\n", + "sample_label(ax, sample, fontsize=35)\n", + "\n", + "########################### Ratio + error bars\n", + "for i in range(len(gen_labels)):\n", + " if i == 0:\n", + " cum_sum = gen_data_list[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,gen_data_list[i]])\n", + "h_gen, _ = np.histogram(cum_sum, bins=bins)\n", + "\n", + "for i in range(len(pf_labels)):\n", + " if i == 0:\n", + " cum_sum = pf_data_list[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,pf_data_list[i]])\n", + "h_pf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + "for i in range(len(mlpf_labels)):\n", + " if i == 0:\n", + " cum_sum = mlpf_data_list[i]\n", + " else:\n", + " cum_sum = np.concatenate([cum_sum,mlpf_data_list[i]])\n", + "h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", + "ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + "ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + "\n", + "# Error bars: Poisson error propagation\n", + "err_pf = np.sqrt(h_pf) / h_gen\n", + "err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + "err_pf[h_gen == 0] = np.nan\n", + "err_mlpf[h_gen == 0] = np.nan\n", + "\n", + "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"PF\"], label=\"PF\", color=color_code[\"PF\"])\n", + "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"MLPF\"], label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + "rax.set_ylabel(\"Reco / Gen\")\n", + "\n", + "rax.set_ylim(0, 2)\n", + "rax.legend(fontsize=25, ncol=2, loc=\"upper center\")\n", + "\n", + "rax.set_xlabel(r\"Particle $\\eta$\", fontsize=35)\n", + "rax.yaxis.set_label_coords(-0.03, 0.9) # move left/right and up/down\n", + "\n", + "plt.text(0.03, 0.82, r\"$p_T>1~$GeV\", transform=ax.transAxes, fontsize=26)\n", + "\n", + "plt.tight_layout()\n", + "\n", + "plt.savefig(f\"./plots/{save_as[sample]}/particle_eta_dist/{sample}_eta_stacked.pdf\")" + ] + }, + { + "cell_type": "markdown", + "id": "a93cf564", + "metadata": {}, + "source": [ + "## phi" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7dbd370", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": 10**16,\n", + " \"cms_pf_ttbar_nopu\": 10**13, \n", + "}\n", + "\n", + "bins = np.linspace(-3,3,50)\n", + "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + "\n", + "fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(12, 12), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + ")\n", + "\n", + "gen_entries = []\n", + "pf_entries = []\n", + "mlpf_entries = []\n", + "for pid in [\n", + " 13,\n", + " 11,\n", + " 130,\n", + " 22, \n", + " 211,\n", + "]:\n", + "\n", + " # Gen\n", + " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_eta\"]) < 2.4) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_phi\"][sel]))\n", " \n", - " if i == 0:\n", - " cum_sum = pf_data_list_central[i]\n", - " else:\n", - " cum_sum = np.concatenate([cum_sum,pf_data_list_central[i]])\n", + " total_yield = len(X_gen)\n", + " gen_entries.append((total_yield, X_gen, pid_to_text[pid], col_pid[pid]))\n", "\n", - " h_pf, _ = np.histogram(cum_sum, bins=bins)\n", + " # PF\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_eta\"]) < 2.4) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", + " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_phi\"][sel]))\n", + "\n", + " total_yield = len(X_pf)\n", + " pf_entries.append((total_yield, X_pf, pid_to_text[pid], col_pid[pid]))\n", " \n", - " ax.errorbar(\n", - " bin_centers, h_pf, yerr=np.sqrt(h_pf),\n", - " fmt=marker_style[\"PF\"], markersize=4, markeredgewidth=1, \n", - " label=f\"PF {pf_labels_central[i]}\",\n", - " color=pf_colors_central[i],\n", - " )\n", + " # MLPF\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.4) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_phi\"][sel]))\n", " \n", - "#### MLPF legend\n", - "mlpf_data_list_forward = [entry[1] for entry in mlpf_entries_forward]\n", - "mlpf_labels_forward = [entry[2] for entry in mlpf_entries_forward]\n", - "mlpf_colors_forward = [entry[3] for entry in mlpf_entries_forward]\n", - "\n", - "# mlpf_data_list_forward.reverse()\n", - "# mlpf_labels_forward.reverse()\n", - "# mlpf_colors_forward.reverse()\n", - "\n", - "tmp = [mlpf_data_list_forward[i] for i in forward_order]\n", - "mlpf_data_list_forward = tmp\n", + " total_yield = len(X_mlpf)\n", + " mlpf_entries.append((total_yield, X_mlpf, pid_to_text[pid], col_pid[pid]))\n", + " \n", + "#### Gen legend\n", + "gen_data_list = [entry[1] for entry in gen_entries]\n", + "gen_labels = [f\"Gen {entry[2]}\" for entry in gen_entries]\n", + "gen_colors = [entry[3] for entry in gen_entries]\n", "\n", - "tmp = [mlpf_labels_forward[i] for i in forward_order]\n", - "mlpf_labels_forward = tmp\n", + "ax.hist(\n", + " gen_data_list,\n", + " bins=bins,\n", + " stacked=True,\n", + " histtype='stepfilled',\n", + " color=gen_colors,\n", + " label=gen_labels,\n", + " alpha=0.3,\n", + " edgecolor='k',\n", + " linewidth=0.8\n", + ")\n", "\n", - "tmp = [mlpf_colors_forward[i] for i in forward_order]\n", - "mlpf_colors_forward = tmp\n", + "#### PF legend\n", + "pf_data_list = [entry[1] for entry in pf_entries]\n", + "pf_labels = [entry[2] for entry in pf_entries]\n", + "pf_colors = [entry[3] for entry in pf_entries]\n", "\n", - "for i in range(len(mlpf_labels_forward)):\n", + "for i in range(len(pf_labels)):\n", " if i == 0:\n", - " cum_sum = mlpf_data_list_forward[i]\n", + " cum_sum = pf_data_list[i]\n", " else:\n", - " cum_sum = np.concatenate([cum_sum,mlpf_data_list_forward[i]])\n", + " cum_sum = np.concatenate([cum_sum,pf_data_list[i]])\n", "\n", - " h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", - " \n", + " h_pf, _ = np.histogram(cum_sum, bins=bins)\n", + " \n", " ax.errorbar(\n", - " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", - " fmt=marker_style[\"MLPF\"], markersize=4, markeredgewidth=1,\n", - "# label=f\"MLPF {mlpf_labels_forward[i]}\",\n", - " color=mlpf_colors_forward[i],\n", + " bin_centers, h_pf, yerr=np.sqrt(h_pf),\n", + "# linestyle=linestyle[\"PF\"],\n", + " fmt=marker_style[\"PF\"], markersize=4, markeredgewidth=1,\n", + " label=f\"PF {pf_labels[i]}\",\n", + " color=pf_colors[i],\n", " )\n", "\n", - "mlpf_data_list_central = [entry[1] for entry in mlpf_entries_central]\n", - "mlpf_labels_central = [entry[2] for entry in mlpf_entries_central]\n", - "mlpf_colors_central = [entry[3] for entry in mlpf_entries_central]\n", + "#### MLPF legend\n", + "mlpf_data_list = [entry[1] for entry in mlpf_entries]\n", + "mlpf_labels = [entry[2] for entry in mlpf_entries]\n", + "mlpf_colors = [entry[3] for entry in mlpf_entries]\n", "\n", - "for i in range(len(mlpf_labels_central)):\n", + "for i in range(len(mlpf_labels)):\n", " if i == 0:\n", - " cum_sum = mlpf_data_list_central[i]\n", + " cum_sum = mlpf_data_list[i]\n", " else:\n", - " cum_sum = np.concatenate([cum_sum,mlpf_data_list_central[i]])\n", + " cum_sum = np.concatenate([cum_sum,mlpf_data_list[i]])\n", "\n", " h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", - " \n", + " \n", " ax.errorbar(\n", " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", + "# linestyle=linestyle[\"PF\"],\n", " fmt=marker_style[\"MLPF\"], markersize=4, markeredgewidth=1,\n", - " label=f\"MLPF {mlpf_labels_central[i]}\",\n", - " color=mlpf_colors_central[i],\n", + " label=f\"MLPF {mlpf_labels[i]}\",\n", + " color=mlpf_colors[i],\n", " )\n", "\n", "############## add legend\n", @@ -1680,15 +1678,17 @@ " labels_MLPF += [lab.replace(\"MLPF \", \"\")]\n", " handles_MLPF += [handles[i]] \n", "\n", + "bbox_to_anchor_xy = (0.64, 0.98)\n", + " \n", "leg1 = ax.legend(\n", " handles=handles_Gen,\n", " labels=labels_Gen,\n", " title=\"Gen\",\n", " ncol=1, # Two-column layout\n", " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25,\n", - " fontsize=20,\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0], bbox_to_anchor_xy[sample][1]),\n", + " title_fontsize=26,\n", + " fontsize=22,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0], bbox_to_anchor_xy[1]),\n", ")\n", "ax.add_artist(leg1)\n", "\n", @@ -1700,9 +1700,9 @@ " title=\"PF\",\n", " ncol=1, # Two-column layout\n", " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25,\n", - " fontsize=20,\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0]+0.22*1, bbox_to_anchor_xy[sample][1]-0.12),\n", + " title_fontsize=26,\n", + " fontsize=22,\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*1, bbox_to_anchor_xy[1]),\n", ")\n", "ax.add_artist(leg2)\n", "\n", @@ -1714,16 +1714,17 @@ " title=\"MLPF\",\n", " ncol=1, # Two-column layout\n", " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25, \n", - " fontsize=20, # Font size of labels\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0]+0.22*2, bbox_to_anchor_xy[sample][1]-0.12), \n", + " title_fontsize=26, \n", + " fontsize=22, # Font size of labels\n", + " bbox_to_anchor=(bbox_to_anchor_xy[0]+0.18*2, bbox_to_anchor_xy[1]), \n", ")\n", "ax.add_artist(leg3)\n", "\n", "ax.set_yscale(\"log\")\n", - "ax.set_ylabel(\"Total number of particles / bin\")\n", + "ax.set_ylabel(\"Total number of particles / bin\", fontsize=35)\n", + "\n", "ax.set_ylim(1, ylim_dict[sample])\n", - "ax.set_xlim(-5,5)\n", + "ax.set_xlim(-3,3)\n", "\n", "cms_label(ax)\n", "if \"qcd\" in sample:\n", @@ -1766,262 +1767,486 @@ "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"MLPF\"], label=\"MLPF\", color=color_code[\"MLPF\"])\n", "\n", "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + "rax.set_ylim(0, 2)\n", "rax.set_ylabel(\"Reco / Gen\")\n", "\n", - "if sample == \"cms_pf_qcd_nopu\":\n", - " rax.set_ylim(0, 2.2)\n", - " rax.legend(fontsize=23, ncol=2, loc=\"upper center\")\n", - "else:\n", - " rax.set_ylim(0,2)\n", - " rax.legend(fontsize=23, ncol=2, loc=\"upper center\")\n", + "rax.set_ylim(0, 2)\n", + "rax.legend(fontsize=25, ncol=2, loc=\"upper center\")\n", + "\n", + "rax.set_xlabel(r\"Particle $\\phi$\", fontsize=35)\n", + "\n", + "plt.text(0.03, 0.82, r\"$p_T>1~$GeV, $|\\eta| < 2.4$\", transform=ax.transAxes, fontsize=26)\n", + "\n", + "plt.savefig(f\"./plots/{save_as[sample]}/particle_phi_dist/{sample}_phi_stacked.pdf\")" + ] + }, + { + "cell_type": "markdown", + "id": "f28230b7", + "metadata": {}, + "source": [ + "# Kinematic distributions" + ] + }, + { + "cell_type": "markdown", + "id": "e2aee193", + "metadata": {}, + "source": [ + "## pT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27043976", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: 10**10,\n", + " 130: 10**8,\n", + " 22: 10**10,\n", + " 11: 10**7,\n", + " 13: 10**6,\n", + " 2: 10**6,\n", + " 1: 10**6, \n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: 10**10,\n", + " 130: 10**8,\n", + " 22: 10**10,\n", + " 11: 10**7,\n", + " 13: 10**6,\n", + " 2: 10**6,\n", + " 1: 10**6, \n", + " }, \n", + "}\n", + "\n", + "rax_ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: (0,2),\n", + " 130: (0,2.2),\n", + " 22: (0,4),\n", + " 11: (0,8),\n", + " 13: (0,4),\n", + " 2: (0,4),\n", + " 1: (0,4), \n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: (0.5,1.5),\n", + " 130: (0,2),\n", + " 22: (0,2.5),\n", + " 11: (0,3),\n", + " 13: (0,2.4),\n", + " 2: (0,2.4),\n", + " 1: (0,2.4), \n", + " }, \n", + "}\n", + "\n", + "bins_pt = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: np.linspace(1,500,41),\n", + " 130: np.linspace(1,200,41),\n", + " 22: np.linspace(1,200,41),\n", + " 11: np.linspace(1,50,21),\n", + " 13: np.linspace(1,50,21),\n", + " 2: np.linspace(1,50,21),\n", + " 1: np.linspace(1,50,21), \n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: np.linspace(1,100,31),\n", + " 130: np.linspace(1,40,31),\n", + " 22: np.linspace(1,40,31),\n", + " 11: np.linspace(1,100,31),\n", + " 13: np.linspace(1,100,31), \n", + " }, \n", + "}\n", + "\n", + "eta_cut = 2.4\n", + "\n", + "for pid in [\n", + " 211,\n", + " 130,\n", + " 22,\n", + " 11,\n", + " 13,\n", + "]:\n", + " fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(10, 10), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + " )\n", + "\n", + " # Extract flattened arrays\n", + " sel = (abs(data_pf[\"GenCands_pid\"]) == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_pt\"][sel]))\n", + "\n", + " sel = (abs(data_pf[\"PFCands_pid\"]) == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1) \n", + " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_pt\"][sel]))\n", + "\n", + " sel = (abs(data_mlpf[\"PFCands_pid\"]) == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_pt\"][sel]))\n", + "\n", + " bins = bins_pt[sample][pid]\n", + " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", + "\n", + " # Histograms\n", + " h_gen, _ = np.histogram(X_gen, bins=bins)\n", + " h_pf, _ = np.histogram(X_pf, bins=bins)\n", + " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", + "\n", + " # Main plot\n", + " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", + " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", + " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + " ax.set_yscale(\"log\")\n", + " ax.set_ylabel(\"Counts\")\n", + " ax.set_ylim(1, ylim_dict[sample][pid])\n", + " ax.set_xlim(1)\n", + " ax.legend()\n", + " cms_label(ax)\n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=24)\n", + " else:\n", + " sample_label(ax, sample, fontsize=30)\n", + " \n", + " plt.text(0.03, 0.82, pid_to_text[pid] + f\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", + "\n", + " # Ratio + error bars\n", + " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + "\n", + " # Error bars: Poisson error propagation\n", + " err_pf = np.sqrt(h_pf) / h_gen\n", + " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + " err_pf[h_gen == 0] = np.nan\n", + " err_mlpf[h_gen == 0] = np.nan\n", "\n", - "rax.set_xlabel(r\"Particle $\\eta$\")\n", - "rax.yaxis.set_label_coords(-0.03, 0.9) # move left/right and up/down\n", + " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", + " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", "\n", - "plt.text(0.03, 0.82, r\"$p_T>1$GeV\", transform=ax.transAxes, fontsize=22)\n", + " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + " rax.set_ylim(rax_ylim_dict[sample][pid])\n", + " rax.set_ylabel(\"Ratio\")\n", " \n", - "plt.savefig(f\"./plots/{save_as[sample]}/particle_eta_dist/{sample}_eta_stacked.pdf\")" + " if sample==\"cms_pf_ttbar_nopu\":\n", + " if ( (pid==211) ):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", + " else:\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + " else:\n", + "# if ( (pid==22) ):\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates \n", + " \n", + " \n", + " rax.set_xlabel(\"$p_T$ (GeV)\")\n", + " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + "\n", + " plt.savefig(f\"./plots/{save_as[sample]}/particle_pt_dist/{sample}_pt_{pid}.pdf\")" ] }, { "cell_type": "markdown", - "id": "a93cf564", + "id": "776e048e", "metadata": {}, "source": [ - "## phi" + "## eta" ] }, { "cell_type": "code", "execution_count": null, - "id": "c7dbd370", - "metadata": {}, + "id": "9e5b80d8", + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ - "ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": 10**16,\n", - " \"cms_pf_ttbar_nopu\": 10**13, \n", - "}\n", - "bbox_to_anchor_xy = {\n", - "# \"cms_pf_qcd_nopu\": (0.25, 0.4),\n", - " \"cms_pf_qcd_nopu\": (0.55, 0.92), \n", - " \"cms_pf_ttbar_nopu\": (0.30, 0.45),\n", + "bins_eta = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: np.linspace(-2.4,2.4,41),\n", + " 130: np.linspace(-2.4,2.4,41),\n", + " 22: np.linspace(-2.4,2.4,41),\n", + " 11: np.linspace(-2.4,2.4,41),\n", + " 13: np.linspace(-2.4,2.4,41),\n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: np.linspace(-5,5,41),\n", + " 130: np.linspace(-5,5,41),\n", + " 22: np.linspace(-5,5,41),\n", + " 11: np.linspace(-5,5,41),\n", + " 13: np.linspace(-5,5,41),\n", + " 2: np.linspace(-5,5,31),\n", + " 1: np.linspace(-5,5,31), \n", + " },\n", "}\n", "\n", - "bins = np.linspace(-3,3,50)\n", - "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + "ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: 10**10,\n", + " 130: 10**8,\n", + " 22: 10**10,\n", + " 11: 10**7,\n", + " 13: 10**6,\n", + " 2: 10**6,\n", + " 1: 10**6, \n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: 10**10,\n", + " 130: 10**8,\n", + " 22: 10**10,\n", + " 11: 10**7,\n", + " 13: 10**6,\n", + " 2: 10**6,\n", + " 1: 10**6, \n", + " }, \n", + "}\n", "\n", - "fig, (ax, rax) = plt.subplots(\n", - " 2, 1, figsize=(10, 10), sharex=True,\n", - " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", - ")\n", + "rax_ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: (0.8,1.6),\n", + " 130: (0,3),\n", + " 22: (0,2),\n", + " 11: (0,2.5),\n", + " 13: (0,5),\n", + " 2: (0,5),\n", + " 1: (0,5), \n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: (0.4,1.6),\n", + " 130: (0,2),\n", + " 22: (0,1.7),\n", + " 11: (0,2),\n", + " 13: (0,2),\n", + " 2: (0,2),\n", + " 1: (0,2), \n", + " }, \n", + "}\n", "\n", - "gen_entries = []\n", - "pf_entries = []\n", - "mlpf_entries = []\n", "for pid in [\n", - " 13,\n", - " 11,\n", - " 130,\n", - " 22, \n", " 211,\n", + " 130,\n", + " 22,\n", + " 11,\n", + " 13,\n", "]:\n", + " fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(10, 10), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + " )\n", "\n", - " # Gen\n", - " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_eta\"]) < 2.5) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", - " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_phi\"][sel]))\n", - " \n", - " total_yield = len(X_gen)\n", - " gen_entries.append((total_yield, X_gen, pid_to_text[pid], col_pid[pid]))\n", + " # Extract flattened arrays\n", + " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_eta\"][sel]))\n", "\n", - " # PF\n", - " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_eta\"]) < 2.5) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", - " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_phi\"][sel]))\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", + " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_eta\"][sel]))\n", "\n", - " total_yield = len(X_pf)\n", - " pf_entries.append((total_yield, X_pf, pid_to_text[pid], col_pid[pid]))\n", - " \n", - " # MLPF\n", - " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_eta\"]) < 2.5) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", - " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_phi\"][sel]))\n", - " \n", - " total_yield = len(X_mlpf)\n", - " mlpf_entries.append((total_yield, X_mlpf, pid_to_text[pid], col_pid[pid]))\n", - " \n", - "#### Gen legend\n", - "gen_data_list = [entry[1] for entry in gen_entries]\n", - "gen_labels = [f\"Gen {entry[2]}\" for entry in gen_entries]\n", - "gen_colors = [entry[3] for entry in gen_entries]\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_eta\"][sel]))\n", "\n", - "ax.hist(\n", - " gen_data_list,\n", - " bins=bins,\n", - " stacked=True,\n", - " histtype='stepfilled',\n", - " color=gen_colors,\n", - " label=gen_labels,\n", - " alpha=0.3,\n", - " edgecolor='k',\n", - " linewidth=0.8\n", - ")\n", + " bins = bins_eta[sample][pid]\n", + " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", "\n", - "#### PF legend\n", - "pf_data_list = [entry[1] for entry in pf_entries]\n", - "pf_labels = [entry[2] for entry in pf_entries]\n", - "pf_colors = [entry[3] for entry in pf_entries]\n", + " # Histograms\n", + " h_gen, _ = np.histogram(X_gen, bins=bins)\n", + " h_pf, _ = np.histogram(X_pf, bins=bins)\n", + " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", "\n", - "for i in range(len(pf_labels)):\n", - " if i == 0:\n", - " cum_sum = pf_data_list[i]\n", + " # Main plot\n", + " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", + " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", + " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + " ax.set_yscale(\"log\")\n", + " ax.set_ylabel(\"Counts\")\n", + " ax.set_ylim(1, ylim_dict[sample][pid])\n", + " ax.set_xlim(-3,3)\n", + " ax.legend()\n", + " cms_label(ax)\n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=24)\n", " else:\n", - " cum_sum = np.concatenate([cum_sum,pf_data_list[i]])\n", + " sample_label(ax, sample, fontsize=30)\n", + " \n", + " plt.text(0.03, 0.82, pid_to_text[pid] + \", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", "\n", - " h_pf, _ = np.histogram(cum_sum, bins=bins)\n", + " # Ratio + error bars\n", + " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + "\n", + " # Error bars: Poisson error propagation\n", + " err_pf = np.sqrt(h_pf) / h_gen\n", + " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + " err_pf[h_gen == 0] = np.nan\n", + " err_mlpf[h_gen == 0] = np.nan\n", + "\n", + " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", + " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + " rax.set_ylim(rax_ylim_dict[sample][pid])\n", + " rax.set_ylabel(\"Ratio\")\n", " \n", - " ax.errorbar(\n", - " bin_centers, h_pf, yerr=np.sqrt(h_pf),\n", - "# linestyle=linestyle[\"PF\"],\n", - " fmt=marker_style[\"PF\"], markersize=4, markeredgewidth=1,\n", - " label=f\"PF {pf_labels[i]}\",\n", - " color=pf_colors[i],\n", - " )\n", + " if sample==\"cms_pf_qcd_nopu\":\n", + " if (pid==211):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", + " else:\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", "\n", - "#### MLPF legend\n", - "mlpf_data_list = [entry[1] for entry in mlpf_entries]\n", - "mlpf_labels = [entry[2] for entry in mlpf_entries]\n", - "mlpf_colors = [entry[3] for entry in mlpf_entries]\n", + " if sample==\"cms_pf_ttbar_nopu\":\n", + " if ( (pid==211) ):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates \n", + " else:\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + " \n", + " rax.set_xlabel(r\"$\\eta$\")\n", + " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + "\n", + " plt.xlim(-5, 5)\n", + " \n", + " plt.savefig(f\"./plots/{save_as[sample]}/particle_eta_dist/{sample}_eta_{pid}.pdf\")" + ] + }, + { + "cell_type": "markdown", + "id": "3e6adb3c", + "metadata": {}, + "source": [ + "## phi" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8e71f97", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: 10**10,\n", + " 130: 10**8,\n", + " 22: 10**10,\n", + " 11: 10**7,\n", + " 13: 10**6,\n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: 10**10,\n", + " 130: 10**8,\n", + " 22: 10**10,\n", + " 11: 10**7,\n", + " 13: 10**6,\n", + " }, \n", + "}\n", "\n", - "for i in range(len(mlpf_labels)):\n", - " if i == 0:\n", - " cum_sum = mlpf_data_list[i]\n", - " else:\n", - " cum_sum = np.concatenate([cum_sum,mlpf_data_list[i]])\n", + "rax_ylim_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: (0.5,2),\n", + " 130: (0.5,1.5),\n", + " 22: (0,2),\n", + " 11: (0,2),\n", + " 13: (0,4),\n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: (0.5,1.6),\n", + " 130: (0.3,1.8),\n", + " 22: (0,2),\n", + " 11: (0,2),\n", + " 13: (0.5,1.6),\n", + " }, \n", + "}\n", "\n", - " h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", - " \n", - " ax.errorbar(\n", - " bin_centers, h_mlpf, yerr=np.sqrt(h_mlpf),\n", - "# linestyle=linestyle[\"PF\"],\n", - " fmt=marker_style[\"MLPF\"], markersize=4, markeredgewidth=1,\n", - " label=f\"MLPF {mlpf_labels[i]}\",\n", - " color=mlpf_colors[i],\n", + "for pid in [\n", + " 211,\n", + " 130,\n", + " 22,\n", + " 11,\n", + " 13,\n", + "]:\n", + " fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(10, 10), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", " )\n", "\n", - "############## add legend\n", - "handles, labels = ax.get_legend_handles_labels()\n", - "\n", - "labels_Gen, labels_PF, labels_MLPF = [], [], []\n", - "handles_Gen, handles_PF, handles_MLPF = [], [], []\n", - "for i, lab in enumerate(labels):\n", - " if lab.startswith(\"Gen\"):\n", - " labels_Gen += [lab.replace(\"Gen \", \"\")]\n", - " handles_Gen += [handles[i]]\n", - " if lab.startswith(\"PF\"):\n", - " labels_PF += [lab.replace(\"PF \", \"\")]\n", - " handles_PF += [handles[i]]\n", - " if lab.startswith(\"MLPF\"):\n", - " labels_MLPF += [lab.replace(\"MLPF \", \"\")]\n", - " handles_MLPF += [handles[i]] \n", + " # Extract flattened arrays\n", + " sel = (data_pf[\"GenCands_pid\"] == pid) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + " X_gen = ak.to_numpy(ak.flatten(data_pf[\"GenCands_phi\"][sel]))\n", "\n", - "leg1 = ax.legend(\n", - " handles=handles_Gen,\n", - " labels=labels_Gen,\n", - " title=\"Gen\",\n", - " ncol=1, # Two-column layout\n", - " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25,\n", - " fontsize=20,\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0], bbox_to_anchor_xy[sample][1]),\n", - ")\n", - "ax.add_artist(leg1)\n", + " sel = (data_pf[\"PFCands_pid\"] == pid) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", + " X_pf = ak.to_numpy(ak.flatten(data_pf[\"PFCands_phi\"][sel]))\n", "\n", - "handles_PF.reverse()\n", - "labels_PF.reverse()\n", - "leg2 = ax.legend(\n", - " handles=handles_PF,\n", - " labels=labels_PF,\n", - " title=\"PF\",\n", - " ncol=1, # Two-column layout\n", - " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25,\n", - " fontsize=20,\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0]+0.22*1, bbox_to_anchor_xy[sample][1]),\n", - ")\n", - "ax.add_artist(leg2)\n", + " sel = (data_mlpf[\"PFCands_pid\"] == pid) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + " X_mlpf = ak.to_numpy(ak.flatten(data_mlpf[\"PFCands_phi\"][sel]))\n", "\n", - "handles_MLPF.reverse()\n", - "labels_MLPF.reverse()\n", - "leg3 = ax.legend(\n", - " handles=handles_MLPF,\n", - " labels=labels_MLPF,\n", - " title=\"MLPF\",\n", - " ncol=1, # Two-column layout\n", - " handletextpad=0.5, # Space between marker and text\n", - " title_fontsize=25, \n", - " fontsize=20, # Font size of labels\n", - " bbox_to_anchor=(bbox_to_anchor_xy[sample][0]+0.22*2, bbox_to_anchor_xy[sample][1]), \n", - ")\n", - "ax.add_artist(leg3)\n", + " bins = bins_phi[sample][pid]\n", + " bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + " bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", "\n", - "ax.set_yscale(\"log\")\n", - "ax.set_ylabel(\"Total number of particles / bin\")\n", - "ax.set_ylim(1, ylim_dict[sample])\n", - "ax.set_xlim(-3,3)\n", + " # Histograms\n", + " h_gen, _ = np.histogram(X_gen, bins=bins)\n", + " h_pf, _ = np.histogram(X_pf, bins=bins)\n", + " h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", "\n", - "cms_label(ax)\n", - "if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=30)\n", - "else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " # Main plot\n", + " ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", + " ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", + " ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", "\n", - "########################### Ratio + error bars\n", - "for i in range(len(gen_labels)):\n", - " if i == 0:\n", - " cum_sum = gen_data_list[i]\n", + " ax.set_yscale(\"log\")\n", + " ax.set_ylabel(\"Counts\")\n", + " ax.set_ylim(1, ylim_dict[sample][pid])\n", + " ax.set_xlim(-3,3)\n", + " ax.legend()\n", + " cms_label(ax)\n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=24)\n", " else:\n", - " cum_sum = np.concatenate([cum_sum,gen_data_list[i]])\n", - "h_gen, _ = np.histogram(cum_sum, bins=bins)\n", + " sample_label(ax, sample, fontsize=30)\n", + " \n", + " plt.text(0.03, 0.82, pid_to_text[pid] + \", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", "\n", - "for i in range(len(pf_labels)):\n", - " if i == 0:\n", - " cum_sum = pf_data_list[i]\n", - " else:\n", - " cum_sum = np.concatenate([cum_sum,pf_data_list[i]])\n", - "h_pf, _ = np.histogram(cum_sum, bins=bins)\n", - " \n", - "for i in range(len(mlpf_labels)):\n", - " if i == 0:\n", - " cum_sum = mlpf_data_list[i]\n", - " else:\n", - " cum_sum = np.concatenate([cum_sum,mlpf_data_list[i]])\n", - "h_mlpf, _ = np.histogram(cum_sum, bins=bins)\n", - " \n", - "ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", - "ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + " # Ratio + error bars\n", + " ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + " ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", "\n", - "# Error bars: Poisson error propagation\n", - "err_pf = np.sqrt(h_pf) / h_gen\n", - "err_mlpf = np.sqrt(h_mlpf) / h_gen\n", - "err_pf[h_gen == 0] = np.nan\n", - "err_mlpf[h_gen == 0] = np.nan\n", + " # Error bars: Poisson error propagation\n", + " err_pf = np.sqrt(h_pf) / h_gen\n", + " err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + " err_pf[h_gen == 0] = np.nan\n", + " err_mlpf[h_gen == 0] = np.nan\n", "\n", - "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"PF\"], label=\"PF\", color=color_code[\"PF\"])\n", - "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, marker=marker_style[\"MLPF\"], label=\"MLPF\", color=color_code[\"MLPF\"])\n", + " rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", + " rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", "\n", - "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", - "rax.set_ylim(0, 2)\n", - "rax.set_ylabel(\"Reco / Gen\")\n", + " rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + " rax.set_ylim(rax_ylim_dict[sample][pid])\n", + " rax.set_ylabel(\"Ratio\")\n", "\n", - "if sample == \"cms_pf_qcd_nopu\":\n", - " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", - "else:\n", - " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + " if sample==\"cms_pf_qcd_nopu\":\n", + " if ( (pid==211) | (pid==130) ):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates \n", + " else:\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", "\n", - "rax.set_xlabel(r\"Particle $\\phi$\")\n", + " if sample==\"cms_pf_ttbar_nopu\":\n", + " if ( (pid==211) | (pid==130) | (pid==13) ):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", + " else:\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates \n", "\n", - "plt.text(0.03, 0.82, r\"$p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=18)\n", + " rax.set_xlabel(r\"$\\phi$\")\n", + " rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", "\n", - "plt.savefig(f\"./plots/{save_as[sample]}/particle_phi_dist/{sample}_phi_stacked.pdf\")" + " plt.savefig(f\"./plots/{save_as[sample]}/particle_phi_dist/{sample}_phi_{pid}.pdf\")" ] }, { @@ -2258,9 +2483,6 @@ "# )\n", "# msk_gen = (data_mlpf[\"GenCands_pt\"] > 48) & (data_mlpf[\"GenCands_pt\"] < 54) & (abs(data_mlpf[\"GenCands_pid\"]) == 13 )\n", "\n", - "\n", - "\n", - "\n", " genpid = 211\n", " msk_gen = (abs(data[\"GenCands_pid\"]) == genpid ) \n", " plt.scatter(\n", @@ -2337,7 +2559,7 @@ "metadata": {}, "outputs": [], "source": [ - "sel_eta_inclusive = True\n", + "sel_eta_inclusive = False\n", "\n", "sel_isolation_gen = False # change that\n", "sel_isolation_reco = False\n", @@ -2349,7 +2571,7 @@ "if sel_eta_inclusive:\n", " eta_cut = float('inf') \n", "else:\n", - " eta_cut = 2.5\n", + " eta_cut = 2.4\n", "\n", "# define baseline kinematic cuts\n", "msk_pf_gen_pt = data_pf[\"GenCands_pt\"] > pt_cut\n", @@ -2513,21 +2735,23 @@ " else:\n", " plt.ylabel(\"Efficiency\", fontsize=35)\n", " plt.xlabel(\"Particle $p_T^{gen}$ (GeV)\", fontsize=35)\n", - " plt.legend(loc=\"upper right\", fontsize=30)\n", + " plt.legend(loc=\"center right\", fontsize=30)\n", " cms_label(ax)\n", " if \"qcd\" in sample:\n", " sample_label(ax, sample, fontsize=24)\n", " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=35)\n", "\n", - " more_text = \", $p_T>1$GeV\"\n", + " more_text = \", $p_T>1~$GeV\"\n", " if not sel_eta_inclusive:\n", - " more_text += r\", $|\\eta|<2.5$\"\n", + " more_text += r\", $|\\eta|<2.4$\"\n", " if sel_isolation_gen:\n", " more_text += r\", $\\mathrm{RelIso03}^{gen}<0.2$\"\n", "\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + more_text, transform=ax.transAxes, fontsize=25)\n", + " plt.text(0.03, 0.85, pid_to_text[pid] + more_text, transform=ax.transAxes, fontsize=30)\n", " \n", + " plt.tight_layout()\n", + " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_pt/{sample}_efficiency_pt_{pid}_matchingPID.pdf\")\n", " else:\n", @@ -2561,7 +2785,7 @@ "if sel_eta_inclusive:\n", " eta_cut = float('inf')\n", "else:\n", - " eta_cut = 2.5\n", + " eta_cut = 2.4\n", "\n", "# define baseline kinematic cuts\n", "msk_pf_gen_pt = data_pf[\"GenCands_pt\"] > pt_cut\n", @@ -2623,7 +2847,6 @@ " \"cms_pf_ttbar_nopu\": {\n", " 211: np.linspace(1,100,31),\n", " 130: np.linspace(1,60,21),\n", - "# 130: bins,\n", " 22: np.linspace(1,40,21),\n", " 11: np.linspace(1,200,21),\n", " 13: np.linspace(1,200,21),\n", @@ -2709,21 +2932,23 @@ " else:\n", " plt.ylabel(\"Fake rate\", fontsize=35) \n", " plt.xlabel(\"Particle $p_T^{reco}$ (GeV)\", fontsize=35)\n", - " plt.legend(loc=\"upper right\", fontsize=30)\n", + " plt.legend(loc=\"center right\", fontsize=30)\n", " cms_label(ax)\n", " if \"qcd\" in sample:\n", " sample_label(ax, sample, fontsize=24)\n", " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=35)\n", " \n", - " more_text = \", $p_T>1$GeV\"\n", + " more_text = \", $p_T>1~$GeV\"\n", " if not sel_eta_inclusive:\n", - " more_text += r\", $|\\eta|<2.5$\"\n", + " more_text += r\", $|\\eta|<2.4$\"\n", " if sel_isolation_reco:\n", " more_text += r\", $\\mathrm{RelIso03}^{reco}<0.2$\"\n", "\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + more_text, transform=ax.transAxes, fontsize=25)\n", + " plt.text(0.03, 0.85, pid_to_text[pid] + more_text, transform=ax.transAxes, fontsize=30)\n", "\n", + " plt.tight_layout()\n", + " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_pt/{sample}_fakerate_pt_{pid}_matchingPID.pdf\")\n", " else:\n", @@ -2745,9 +2970,9 @@ "metadata": {}, "outputs": [], "source": [ - "sel_eta_inclusive = True\n", + "sel_eta_inclusive = False\n", "\n", - "# if you want ti impose any isolation requirement before the matching\n", + "# if you want to impose any isolation requirement before the matching\n", "sel_isolation_gen = False\n", "sel_isolation_reco = False\n", "\n", @@ -2758,7 +2983,7 @@ "if sel_eta_inclusive:\n", " eta_cut = float('inf') \n", "else:\n", - " eta_cut = 2.5\n", + " eta_cut = 2.4\n", "\n", "# define baseline kinematic cuts\n", "msk_pf_gen_pt = data_pf[\"GenCands_pt\"] > pt_cut\n", @@ -2821,33 +3046,34 @@ "cell_type": "code", "execution_count": null, "id": "582ea637-62cb-4dd0-a0b5-6004f888f485", - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ "bins_pt = {\n", " \"cms_pf_qcd_nopu\": {\n", " 211: np.linspace(1,500,41),\n", - " 130: np.linspace(1,200,41),\n", + " 130: np.linspace(1,60,21),\n", " 22: np.linspace(1,200,41),\n", " 11: np.linspace(1,50,21),\n", " 13: np.linspace(1,50,21),\n", " },\n", " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(1,100,11),\n", - " 130: np.linspace(1,20,11),\n", - " 22: np.linspace(1,40,21),\n", - " 11: np.linspace(1,200,21),\n", - " 13: np.linspace(1,200,21),\n", + " 211: np.linspace(1,100,31),\n", + " 130: np.linspace(1,100,31),\n", + " 22: np.linspace(1,60,21),\n", + " 11: np.linspace(1,100,21),\n", + " 13: np.linspace(1,100,21),\n", " }, \n", "}\n", "\n", - "\n", "match_pid = False\n", "\n", "for pid in [\n", - " 211,\n", + "# 211,\n", " 130,\n", - " 22,\n", + "# 22,\n", "# 11,\n", "# 13,\n", "]:\n", @@ -2939,16 +3165,15 @@ " plt.xlabel(\"Particle $p_T^{gen}$ (GeV)\", fontsize=35)\n", " plt.legend(loc=\"best\", fontsize=30)\n", " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=30)\n", "\n", - " if eta_cut == 2.5:\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta|>2.5$\", transform=ax.transAxes, fontsize=25)\n", + " if sel_eta_inclusive:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", " else:\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV\", transform=ax.transAxes, fontsize=25)\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta|<2.4$\", transform=ax.transAxes, fontsize=25)\n", " \n", + "# plt.text(0.03, 0.77, r\"matching condition ($\\Delta R<0.05$)\", transform=ax.transAxes, fontsize=20)\n", + " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_pt/{sample}_efficiency_pt_{pid}_matchingPID.pdf\")\n", " else:\n", @@ -2967,19 +3192,19 @@ " if match_pid:\n", " plt.ylabel(\"Misidentification\", fontsize=35)\n", " else:\n", - " plt.ylabel(\"Fake rate\", fontsize=35) \n", + "# plt.ylabel(\"Fake rate\", fontsize=35)\n", + " plt.ylabel(\"Misidentification rate\", fontsize=35) \n", " plt.xlabel(\"Particle $p_T^{reco}$ (GeV)\", fontsize=35)\n", " plt.legend(loc=\"best\", fontsize=30)\n", " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=30)\n", " \n", " if sel_eta_inclusive:\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV\", transform=ax.transAxes, fontsize=25)\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", " else:\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta|>2.5$\", transform=ax.transAxes, fontsize=25)\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta|<2.4$\", transform=ax.transAxes, fontsize=25)\n", + " \n", + "# plt.text(0.03, 0.77, r\"matching condition ($\\Delta R<0.05$)\", transform=ax.transAxes, fontsize=20)\n", " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_pt/{sample}_fakerate_pt_{pid}_matchingPID.pdf\")\n", @@ -3006,10 +3231,10 @@ "\n", "for pid in [\n", " 211,\n", - "# 130,\n", - "# 22,\n", - "# 11,\n", - "# 13,\n", + " 130,\n", + " 22,\n", + " 11,\n", + " 13,\n", "]:\n", "\n", " #pick genparticles, and genparticles matched to reco\n", @@ -3097,12 +3322,12 @@ " plt.xlabel(\"Particle $\\eta^{gen}$\", fontsize=35)\n", " plt.legend(loc=\"upper right\", fontsize=30)\n", " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=30)\n", " \n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=25)\n", + " if sel_eta_inclusive:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", + " else:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta|>2.4$\", transform=ax.transAxes, fontsize=25)\n", " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_eta/{sample}_efficiency_eta_{pid}_matchingPID.pdf\")\n", @@ -3129,12 +3354,12 @@ " plt.xlabel(\"Particle $\\eta^{reco}$\", fontsize=35)\n", " plt.legend(loc=\"upper right\", fontsize=30)\n", " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=30)\n", " \n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=25)\n", + " if sel_eta_inclusive:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", + " else:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta|>2.4$\", transform=ax.transAxes, fontsize=25)\n", " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_eta/{sample}_fakerate_eta_{pid}_matchingPID.pdf\")\n", @@ -3157,7 +3382,7 @@ "metadata": {}, "outputs": [], "source": [ - "match_pid = True\n", + "match_pid = False\n", "\n", "for pid in [\n", " 211,\n", @@ -3250,14 +3475,15 @@ " else:\n", " plt.ylabel(\"Efficiency\", fontsize=35)\n", " plt.xlabel(\"Particle $\\phi^{gen}$\")\n", - " plt.legend(loc=\"best\")\n", - " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=25)\n", + " plt.legend(loc=\"best\")\n", + " cms_label(ax)\n", + " sample_label(ax, sample, fontsize=30)\n", " \n", + " if sel_eta_inclusive:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", + " else:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta|>2.4$\", transform=ax.transAxes, fontsize=25)\n", + " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_phi/{sample}_efficiency_phi_{pid}_matchingPID.pdf\")\n", " else:\n", @@ -3280,14 +3506,13 @@ " plt.xlabel(\"Particle $\\phi^{reco}$\")\n", " plt.legend(loc=\"best\")\n", " cms_label(ax)\n", - " \n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", - " else:\n", - " sample_label(ax, sample, fontsize=30)\n", + " sample_label(ax, sample, fontsize=30)\n", "\n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=25)\n", - " \n", + " if sel_eta_inclusive:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", + " else:\n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta|>2.4$\", transform=ax.transAxes, fontsize=25)\n", + " \n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_phi/{sample}_fakerate_phi_{pid}_matchingPID.pdf\")\n", " else:\n", @@ -3363,18 +3588,18 @@ "\n", "bins_eta_2dhist = {\n", " \"cms_pf_qcd_nopu\": {\n", - " 211: np.linspace(-5,5,41),\n", - " 130: np.linspace(-5,5,41),\n", - " 22: np.linspace(-5,5,41),\n", - " 11: np.linspace(-5,5,41),\n", - " 13: np.linspace(-5,5,21),\n", + " 211: np.linspace(-2.4,2.4,41),\n", + " 130: np.linspace(-2.4,2.4,41),\n", + " 22: np.linspace(-2.4,2.4,41),\n", + " 11: np.linspace(-2.4,2.4,41),\n", + " 13: np.linspace(-2.4,2.4,41),\n", " },\n", " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(-5,5,41),\n", - " 130: np.linspace(-5,5,41),\n", - " 22: np.linspace(-5,5,41),\n", - " 11: np.linspace(-5,5,41),\n", - " 13: np.linspace(-5,5,41),\n", + " 211: np.linspace(-2.4,2.4,41),\n", + " 130: np.linspace(-2.4,2.4,41),\n", + " 22: np.linspace(-2.4,2.4,41),\n", + " 11: np.linspace(-2.4,2.4,41),\n", + " 13: np.linspace(-2.4,2.4,41),\n", " },\n", "}" ] @@ -3389,9 +3614,9 @@ "match_pid = False\n", "\n", "for pid in [\n", - "# 211,\n", - " 130,\n", - "# 22,\n", + " 211,\n", + "# 130,\n", + " 22,\n", "]:\n", " \n", " #################### PF eff\n", @@ -3425,20 +3650,25 @@ " X, Y = np.meshgrid(pt_edges, eta_edges, indexing='ij') # shape (41, 31)\n", "\n", " # Now everything matches → C.shape = (40,30), X.shape = (41,31)\n", - " fig, ax = plt.subplots(figsize=(12, 10))\n", + " fig, ax = plt.subplots(figsize=(11, 9))\n", " pcm = ax.pcolormesh(X, Y, eff_vals, cmap=\"viridis\", vmin=0, vmax=1, shading=\"auto\", edgecolor='face')\n", " cbar = fig.colorbar(pcm, ax=ax)\n", " if match_pid:\n", - " cbar.set_label(\"PF \" + pid_to_text[pid] + \" efficiency (w/ matched PID)\", fontsize=30) \n", + " cbar.set_label(\"PF \" + pid_to_text[pid] + \" efficiency (w/ matched PID)\", fontsize=35) \n", " else:\n", - " cbar.set_label(\"PF \" + pid_to_text[pid] + \" efficiency\", fontsize=30)\n", + " cbar.set_label(\"PF \" + pid_to_text[pid] + \" efficiency\", fontsize=35)\n", " \n", - " ax.set_xlabel(r\"Particle $p_T^{gen}$ (GeV)\", fontsize=35)\n", - " ax.set_ylabel(r\"Particle $\\eta^{gen}$\", fontsize=35)\n", + " ax.set_xlabel(r\"Particle $p_T^{gen}$ (GeV)\", fontsize=40)\n", + " ax.set_ylabel(r\"Particle $\\eta^{gen}$\", fontsize=40)\n", "\n", " cms_label(ax)\n", " ax.grid(False)\n", " \n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=35)\n", + " else:\n", + " sample_label(ax, sample, fontsize=35)\n", + " \n", " plt.tight_layout()\n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_2d/PF_{sample}_eff_pteta_{pid}_matchingPID.pdf\")\n", @@ -3475,20 +3705,25 @@ " X, Y = np.meshgrid(pt_edges, eta_edges, indexing='ij') # shape (41, 31)\n", "\n", " # Now everything matches → C.shape = (40,30), X.shape = (41,31)\n", - " fig, ax = plt.subplots(figsize=(12, 10))\n", + " fig, ax = plt.subplots(figsize=(11, 9))\n", " pcm = ax.pcolormesh(X, Y, fake_vals, cmap=\"viridis\", vmin=0, vmax=1, shading=\"auto\", edgecolor='face')\n", " cbar = fig.colorbar(pcm, ax=ax)\n", "\n", " if match_pid:\n", - " cbar.set_label(\"PF \" + pid_to_text[pid] + \" misidentification\", fontsize=30) \n", + " cbar.set_label(\"PF \" + pid_to_text[pid] + \" misidentification\", fontsize=35) \n", " else:\n", - " cbar.set_label(\"PF \" + pid_to_text[pid] + \" fake rate\", fontsize=30) \n", + " cbar.set_label(\"PF \" + pid_to_text[pid] + \" fake rate\", fontsize=35) \n", "\n", - " ax.set_xlabel(r\"Particle $p_T^{reco}$ (GeV)\", fontsize=35)\n", - " ax.set_ylabel(r\"Particle $\\eta^{reco}$\", fontsize=35)\n", + " ax.set_xlabel(r\"Particle $p_T^{reco}$ (GeV)\", fontsize=40)\n", + " ax.set_ylabel(r\"Particle $\\eta^{reco}$\", fontsize=40)\n", " \n", " cms_label(ax)\n", " ax.grid(False)\n", + " \n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=35, color=\"white\")\n", + " else:\n", + " sample_label(ax, sample, fontsize=35, color=\"white\")\n", " \n", " plt.tight_layout()\n", " if match_pid:\n", @@ -3526,19 +3761,24 @@ " X, Y = np.meshgrid(pt_edges, eta_edges, indexing='ij') # shape (41, 31)\n", "\n", " # Now everything matches → C.shape = (40,30), X.shape = (41,31)\n", - " fig, ax = plt.subplots(figsize=(12, 10))\n", + " fig, ax = plt.subplots(figsize=(11, 9))\n", " pcm = ax.pcolormesh(X, Y, eff_vals, cmap=\"viridis\", vmin=0, vmax=1, shading=\"auto\", edgecolor='face')\n", " cbar = fig.colorbar(pcm, ax=ax)\n", " if match_pid:\n", - " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" efficiency (w/ matched PID)\", fontsize=30) \n", + " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" efficiency (w/ matched PID)\", fontsize=35) \n", " else:\n", - " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" efficiency\", fontsize=30)\n", + " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" efficiency\", fontsize=35)\n", " \n", - " ax.set_xlabel(r\"Particle $p_T^{gen}$ (GeV)\", fontsize=35)\n", - " ax.set_ylabel(r\"Particle $\\eta^{gen}$\", fontsize=35)\n", + " ax.set_xlabel(r\"Particle $p_T^{gen}$ (GeV)\", fontsize=40)\n", + " ax.set_ylabel(r\"Particle $\\eta^{gen}$\", fontsize=40)\n", " \n", " cms_label(ax)\n", " ax.grid(False)\n", + " \n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=35)\n", + " else:\n", + " sample_label(ax, sample, fontsize=35)\n", " \n", " plt.tight_layout()\n", " if match_pid:\n", @@ -3576,22 +3816,28 @@ " X, Y = np.meshgrid(pt_edges, eta_edges, indexing='ij') # shape (41, 31)\n", "\n", " # Now everything matches → C.shape = (40,30), X.shape = (41,31)\n", - " fig, ax = plt.subplots(figsize=(12, 10))\n", + " fig, ax = plt.subplots(figsize=(11, 9))\n", " pcm = ax.pcolormesh(X, Y, fake_vals, cmap=\"viridis\", vmin=0, vmax=1, shading=\"auto\", edgecolor='face')\n", " cbar = fig.colorbar(pcm, ax=ax)\n", "\n", " if match_pid:\n", - " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" misidentification\", fontsize=30) \n", + " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" misidentification\", fontsize=35) \n", " else:\n", - " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" fake rate\", fontsize=30) \n", + " cbar.set_label(\"MLPF \" + pid_to_text[pid] + \" fake rate\", fontsize=35) \n", "\n", - " ax.set_xlabel(r\"Particle $p_T^{reco}$ (GeV)\", fontsize=35)\n", - " ax.set_ylabel(r\"Particle $\\eta^{reco}$\", fontsize=35)\n", + " ax.set_xlabel(r\"Particle $p_T^{reco}$ (GeV)\", fontsize=40)\n", + " ax.set_ylabel(r\"Particle $\\eta^{reco}$\", fontsize=40)\n", " \n", " cms_label(ax)\n", " ax.grid(False)\n", - " \n", + " \n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=35, color=\"white\")\n", + " else:\n", + " sample_label(ax, sample, fontsize=35, color=\"white\")\n", + "\n", " plt.tight_layout()\n", + "\n", " if match_pid:\n", " plt.savefig(f\"./plots/{save_as[sample]}/particle_eff_fakerate_2d/MLPF_{sample}_fakerate_pteta_{pid}_matchingPID.pdf\")\n", " else:\n", @@ -3606,6 +3852,82 @@ "## pT resolution" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "a9646a9a", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "########## checking the isolation distributons first\n", + "\n", + "b_dict = {\n", + " \"cms_pf_qcd_nopu\": {\n", + " 211: np.linspace(0,10,200),\n", + " 130: np.linspace(0,10,200),\n", + " 22: np.linspace(0,20,200),\n", + " 11: np.linspace(0,5,100),\n", + " 13: np.linspace(0,2,100),\n", + " },\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: np.linspace(0,2,100),\n", + " 130: np.linspace(0,2,100),\n", + " 22: np.linspace(0,2,100),\n", + " 11: np.linspace(0,2,100),\n", + " 13: np.linspace(0,2,100),\n", + " }, \n", + "}\n", + "\n", + "ylim_dict = {\n", + " \"cms_pf_ttbar_nopu\": {\n", + " 211: 10**8,\n", + " 130: 5*10**7,\n", + " 22: 5*10**7,\n", + " 11: 10**5,\n", + " 13: 10**5,\n", + " }, \n", + "}\n", + "\n", + "for pid in [\n", + " 211,\n", + " 130,\n", + " 22,\n", + " 11,\n", + " 13,\n", + "]:\n", + " \n", + " fig = plt.figure()\n", + " ax = plt.axes()\n", + " \n", + " msk_gen = (data_pf[\"GenCands_pid\"]==pid) & (data_pf[\"GenCands_pt\"]>1)# & (abs(data_pf[\"GenCands_eta\"])<2.4)\n", + " plt.hist(ak.flatten(data_pf[\"GenCands_reliso03\"][msk_gen]), bins=b_dict[sample][pid], histtype=\"step\",lw=2, label=\"Gen\", color=color_code[\"Gen\"]); \n", + " \n", + " msk_pf = (data_pf[\"PFCands_pid\"]==pid) & (data_pf[\"PFCands_pt\"]>1)# & (abs(data_pf[\"PFCands_eta\"])<2.4)\n", + " plt.hist(ak.flatten(data_pf[\"PFCands_reliso03\"][msk_pf]), bins=b_dict[sample][pid], histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"]);\n", + " \n", + " msk_mlpf = (data_mlpf[\"PFCands_pid\"]==pid) & (data_mlpf[\"PFCands_pt\"]>1)# & (abs(data_mlpf[\"PFCands_eta\"])<2.4)\n", + " plt.hist(ak.flatten(data_mlpf[\"PFCands_reliso03\"][msk_mlpf]), bins=b_dict[sample][pid], histtype=\"step\",lw=2, label=\"MLPF\", color=color_code[\"MLPF\"]);\n", + "\n", + " ax.axvline(0.2, color=\"grey\", linestyle=\"--\", label=\"RelIso03=0.2\") \n", + " \n", + " plt.yscale(\"log\")\n", + " plt.ylim(1, ylim_dict[sample][pid])\n", + " plt.xlabel(\"RelIso03\")\n", + " plt.legend()\n", + " cms_label(ax)\n", + " if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=24)\n", + " else:\n", + " sample_label(ax, sample, fontsize=30)\n", + " \n", + " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV\", transform=ax.transAxes, fontsize=25)\n", + "# plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1~$GeV, $|\\eta| < 2.4$\", transform=ax.transAxes, fontsize=25)\n", + " \n", + " plt.savefig(f\"./plots/{save_as[sample]}/lepisolation/{sample}_{pid}.pdf\")" + ] + }, { "cell_type": "code", "execution_count": null, @@ -3615,7 +3937,7 @@ "source": [ "dr_circle_size = 0.15 ###### specify dR threshold\n", "pid = 211\n", - "pt_cut = 0.5\n", + "pt_cut = 1\n", "\n", "iev = 4\n", "\n", @@ -3670,7 +3992,7 @@ "if sel_eta_inclusive:\n", " eta_cut = float('inf') \n", "else:\n", - " eta_cut = 2.5\n", + " eta_cut = 2.4\n", "\n", "# define baseline kinematic cuts\n", "msk_pf_gen_pt = data_pf[\"GenCands_pt\"] > pt_cut\n", @@ -3742,11 +4064,11 @@ "if sel_isolation is True:\n", " ylim_dict = {\n", " \"cms_pf_ttbar_nopu\": {\n", - " 211: 10**6,\n", - " 130: 10**4,\n", - " 22: 10**4,\n", - " 11: 10**3,\n", - " 13: 10**3,\n", + " 211: 10**7,\n", + " 130: 10**5,\n", + " 22: 10**5,\n", + " 11: 10**5,\n", + " 13: 10**5,\n", " }, \n", " }\n", "else:\n", @@ -3765,8 +4087,8 @@ " 211,\n", " 130,\n", " 22,\n", - "# 11,\n", - "# 13,\n", + " 11,\n", + " 13,\n", "]:\n", " \n", " fig = plt.figure()\n", @@ -3782,106 +4104,203 @@ " else:\n", " sample_label(ax, sample, fontsize=30)\n", " \n", - " more_text = \", $p_T>1$GeV\"\n", - " save_ = \"\"\n", + " more_text = \", $p_T>1~$GeV\"\n", + " save_ = \"_\"\n", " \n", " if sel_eta_inclusive is not True:\n", - " more_text += \", $|\\eta| < 2.5$\"\n", - " save_ += \"_central\"\n", + " more_text += \", $|\\eta| < 2.4$\"\n", + " save_ += \"central\"\n", " else:\n", - " save_ += \"_inclusive\"\n", + " save_ += \"inclusive\"\n", " \n", " if sel_isolation is True:\n", " more_text += \", $\\mathrm{RelIso03}^{(gen)}<0.2$\"\n", - " plt.legend(loc=\"lower left\")\n", + " plt.legend(loc=\"center right\")\n", " save_ += \"_isolated\"\n", " else:\n", " plt.legend()\n", " \n", " plt.text(0.03, 0.87, pid_to_text[pid] + more_text, transform=ax.transAxes, fontsize=25)\n", "\n", - " plt.savefig(f\"./plots/{save_as[sample]}/particle_ptres/{sample}_ptres_{pid}_{save_}.pdf\")" + " plt.savefig(f\"./plots/{save_as[sample]}/particle_ptres/{sample}_ptres_{pid}{save_}.pdf\")" + ] + }, + { + "cell_type": "markdown", + "id": "011fcadf", + "metadata": {}, + "source": [ + "# Forward region" ] }, { "cell_type": "code", "execution_count": null, - "id": "a9646a9a", + "id": "a1f99306", "metadata": {}, "outputs": [], "source": [ - "########## checking the isolation distributons\n", + "bins_pt = np.linspace(1,40,31)\n", "\n", - "b_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: np.linspace(0,10,200),\n", - " 130: np.linspace(0,10,200),\n", - " 22: np.linspace(0,20,200),\n", - " 11: np.linspace(0,5,100),\n", - " 13: np.linspace(0,2,100),\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: np.linspace(0,2,100),\n", - " 130: np.linspace(0,2,100),\n", - " 22: np.linspace(0,2,100),\n", - " 11: np.linspace(0,2,100),\n", - " 13: np.linspace(0,2,100),\n", - " }, \n", - "}\n", + "ylim_dict = 10**8\n", + "rax_ylim_dict = (0,2)\n", "\n", - "ylim_dict = {\n", - " \"cms_pf_qcd_nopu\": {\n", - " 211: 10**8,\n", - " 130: 10**6,\n", - " 22: 10**6,\n", - " 11: 10**6,\n", - " 13: 10**6,\n", - " },\n", - " \"cms_pf_ttbar_nopu\": {\n", - " 211: 10**6,\n", - " 130: 5*10**5,\n", - " 22: 5*10**5,\n", - " 11: 10**3,\n", - " 13: 10**3,\n", - " }, \n", - "}\n", + "fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(10, 10), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + ")\n", "\n", - "for pid in [\n", - " 211,\n", - " 130,\n", - " 22,\n", - "# 11,\n", - "# 13,\n", - "]:\n", - " \n", - " fig = plt.figure()\n", - " ax = plt.axes()\n", - " \n", - " msk_gen = (data_pf[\"GenCands_pid\"]==pid) & (data_pf[\"GenCands_pt\"]>1)# & (abs(data_pf[\"GenCands_eta\"])<2.5)\n", - " plt.hist(ak.flatten(data_pf[\"GenCands_reliso03\"][msk_gen]), bins=b_dict[sample][pid], histtype=\"step\",lw=2, label=\"Gen\", color=color_code[\"Gen\"]); \n", - " \n", - " msk_pf = (data_pf[\"PFCands_pid\"]==pid) & (data_pf[\"PFCands_pt\"]>1)# & (abs(data_pf[\"PFCands_eta\"])<2.5)\n", - " plt.hist(ak.flatten(data_pf[\"PFCands_reliso03\"][msk_pf]), bins=b_dict[sample][pid], histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"]);\n", - " \n", - " msk_mlpf = (data_mlpf[\"PFCands_pid\"]==pid) & (data_mlpf[\"PFCands_pt\"]>1)# & (abs(data_mlpf[\"PFCands_eta\"])<2.5)\n", - " plt.hist(ak.flatten(data_mlpf[\"PFCands_reliso03\"][msk_mlpf]), bins=b_dict[sample][pid], histtype=\"step\",lw=2, label=\"MLPF\", color=color_code[\"MLPF\"]);\n", + "# Extract flattened arrays\n", + "sel = (abs(data_pf[\"GenCands_eta\"]) > 3) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + "X_gen = ak.to_numpy(ak.flatten(abs(data_pf[\"GenCands_pt\"][sel])))\n", "\n", - " ax.axvline(0.2, color=\"grey\", linestyle=\"--\", label=\"RelIso03=0.2\") \n", - " \n", - " plt.yscale(\"log\")\n", - " plt.ylim(1, ylim_dict[sample][pid])\n", - " plt.xlabel(\"RelIso03\")\n", - " plt.legend()\n", - " cms_label(ax)\n", - " if \"qcd\" in sample:\n", - " sample_label(ax, sample, fontsize=24)\n", + "sel = (abs(data_pf[\"PFCands_eta\"]) > 3) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", + "X_pf = ak.to_numpy(ak.flatten(abs(data_pf[\"PFCands_pt\"][sel])))\n", + "\n", + "sel = (abs(data_mlpf[\"PFCands_eta\"]) > 3)& (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + "X_mlpf = ak.to_numpy(ak.flatten(abs(data_mlpf[\"PFCands_pt\"][sel])))\n", + "\n", + "bins = bins_pt\n", + "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + "bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", + "\n", + "# Histograms\n", + "h_gen, _ = np.histogram(X_gen, bins=bins)\n", + "h_pf, _ = np.histogram(X_pf, bins=bins)\n", + "h_mlpf, _ = np.histogram(X_mlpf, bins=bins)\n", + "\n", + "# Main plot\n", + "ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", + "ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", + "ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + "ax.set_yscale(\"log\")\n", + "ax.set_ylabel(\"Counts\")\n", + "ax.set_ylim(1, ylim_dict)\n", + "ax.set_xlim(0)\n", + "ax.legend(loc=\"lower left\")\n", + "cms_label(ax)\n", + "\n", + "sample_label(ax, sample, fontsize=30)\n", + "\n", + "plt.text(0.03, 0.82, \"all particles\" + \", $p_T>1~$GeV, $|\\eta| > 3$\", transform=ax.transAxes, fontsize=25)\n", + "\n", + "# Ratio + error bars\n", + "ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + "ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + "\n", + "# Error bars: Poisson error propagation\n", + "err_pf = np.sqrt(h_pf) / h_gen\n", + "err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + "err_pf[h_gen == 0] = np.nan\n", + "err_mlpf[h_gen == 0] = np.nan\n", + "\n", + "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", + "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + "rax.set_ylim(rax_ylim_dict)\n", + "rax.set_ylabel(\"Ratio\")\n", + "\n", + "rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + "\n", + "rax.set_xlabel(\"$p_T$ (GeV)\")\n", + "rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + "\n", + "plt.savefig(f\"./plots/{save_as[sample]}/forward_region/{sample}_pt_all.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81bcb894", + "metadata": {}, + "outputs": [], + "source": [ + "bins_eta = np.linspace(3,5,12)\n", + "\n", + "ylim_dict = 10**8\n", + "rax_ylim_dict = (0,2)\n", + "\n", + "fig, (ax, rax) = plt.subplots(\n", + " 2, 1, figsize=(10, 10), sharex=True,\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.05}\n", + ")\n", + "\n", + "# Extract flattened arrays\n", + "sel = (abs(data_pf[\"GenCands_eta\"]) > 3) & (abs(data_pf[\"GenCands_pt\"]) > 1)\n", + "X_gen = ak.to_numpy(ak.flatten(abs(data_pf[\"GenCands_eta\"][sel])))\n", + "w_gen = ak.to_numpy(ak.flatten(abs(data_pf[\"GenCands_pt\"][sel])))\n", + "\n", + "sel = (abs(data_pf[\"PFCands_eta\"]) > 3) & (abs(data_pf[\"PFCands_pt\"]) > 1)\n", + "X_pf = ak.to_numpy(ak.flatten(abs(data_pf[\"PFCands_eta\"][sel])))\n", + "w_pf = ak.to_numpy(ak.flatten(abs(data_pf[\"PFCands_pt\"][sel]))) \n", + "\n", + "sel = (abs(data_mlpf[\"PFCands_eta\"]) > 3) & (abs(data_mlpf[\"PFCands_pt\"]) > 1)\n", + "X_mlpf = ak.to_numpy(ak.flatten(abs(data_mlpf[\"PFCands_eta\"][sel])))\n", + "w_mlpf = ak.to_numpy(ak.flatten(abs(data_mlpf[\"PFCands_pt\"][sel])))\n", + "\n", + "bins = bins_eta\n", + "bin_centers = 0.5 * (bins[:-1] + bins[1:])\n", + "bin_widths = 0.5 * (bins[1:] - bins[:-1])\n", + "\n", + "# Histograms\n", + "h_gen, _ = np.histogram(X_gen, bins=bins, weights=w_gen)\n", + "h_pf, _ = np.histogram(X_pf, bins=bins, weights=w_pf)\n", + "h_mlpf, _ = np.histogram(X_mlpf, bins=bins, weights=w_mlpf)\n", + "\n", + "# Main plot\n", + "ax.hist(bins[:-1], bins=bins, weights=h_gen, histtype=\"step\", lw=2, label=\"Gen\", linestyle=\"--\", color=color_code[\"Gen\"])\n", + "ax.hist(bins[:-1], bins=bins, weights=h_pf, histtype=\"step\", lw=2, label=\"PF\", linestyle=\"--\", color=color_code[\"PF\"])\n", + "ax.hist(bins[:-1], bins=bins, weights=h_mlpf, histtype=\"step\", lw=2, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + "ax.set_yscale(\"log\")\n", + "ax.set_ylabel(\"Counts (weighted by pT)\")\n", + "ax.set_ylim(1, ylim_dict)\n", + "ax.legend()\n", + "cms_label(ax)\n", + "if \"qcd\" in sample:\n", + " sample_label(ax, sample, fontsize=24)\n", + "else:\n", + " sample_label(ax, sample, fontsize=30)\n", + "\n", + "plt.text(0.03, 0.82, \"all particles\" + \", $p_T>1~$GeV, $|\\eta| > 3$\", transform=ax.transAxes, fontsize=25) \n", + "\n", + "# Ratio + error bars\n", + "ratio_pf = np.divide(h_pf, h_gen, out=np.full_like(h_pf, np.nan, dtype=float), where=h_gen > 0)\n", + "ratio_mlpf = np.divide(h_mlpf, h_gen, out=np.full_like(h_mlpf, np.nan, dtype=float), where=h_gen > 0)\n", + "\n", + "# Error bars: Poisson error propagation\n", + "err_pf = np.sqrt(h_pf) / h_gen\n", + "err_mlpf = np.sqrt(h_mlpf) / h_gen\n", + "err_pf[h_gen == 0] = np.nan\n", + "err_mlpf[h_gen == 0] = np.nan\n", + "\n", + "rax.errorbar(bin_centers, ratio_pf, yerr=err_pf, fmt='o', markersize=4, markeredgewidth=1, label=\"PF\", color=color_code[\"PF\"])\n", + "rax.errorbar(bin_centers, ratio_mlpf, yerr=err_mlpf, fmt='o', markersize=4, markeredgewidth=1, label=\"MLPF\", color=color_code[\"MLPF\"])\n", + "\n", + "rax.axhline(1.0, color=\"gray\", linestyle=\"--\")\n", + "rax.set_ylim(rax_ylim_dict)\n", + "rax.set_ylabel(\"Ratio\")\n", + "\n", + "if sample==\"cms_pf_qcd_nopu\":\n", + " if (pid==211):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates\n", " else:\n", - " sample_label(ax, sample, fontsize=30)\n", - " \n", - " plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV\", transform=ax.transAxes, fontsize=25)\n", - "# plt.text(0.03, 0.87, pid_to_text[pid] + r\", $p_T>1$GeV, $|\\eta| < 2.5$\", transform=ax.transAxes, fontsize=25)\n", - " \n", - " plt.savefig(f\"./plots/{save_as[sample]}/lepisolation/{sample}_{pid}.pdf\")" + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + "\n", + "if sample==\"cms_pf_ttbar_nopu\":\n", + " if ( (pid==211) ):\n", + " rax.yaxis.set_label_coords(-0.065, 0.7) # x, y in axis coordinates \n", + " else:\n", + " rax.yaxis.set_label_coords(-0.04, 0.7) # x, y in axis coordinates\n", + "\n", + "rax.set_xlabel(r\"Particle $|\\eta|$\")\n", + "rax.legend(fontsize=20, ncol=2, loc=\"upper center\")\n", + "\n", + "# plt.xlim(-2.5, 2.5)\n", + "\n", + "plt.savefig(f\"./plots/{save_as[sample]}/forward_region/{sample}_eta_all.pdf\")" ] }, { @@ -4011,7 +4430,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.10" + "version": "3.9.0" } }, "nbformat": 4,