diff --git a/src/HHbbVV/corrections/ratio_feb20.root b/src/HHbbVV/corrections/ratio_feb20.root deleted file mode 100644 index 145d2e47..00000000 Binary files a/src/HHbbVV/corrections/ratio_feb20.root and /dev/null differ diff --git a/src/HHbbVV/corrections/ratio_june9.root b/src/HHbbVV/corrections/ratio_june9.root new file mode 100644 index 00000000..e0fd9b80 Binary files /dev/null and b/src/HHbbVV/corrections/ratio_june9.root differ diff --git a/src/HHbbVV/postprocessing/corrections.py b/src/HHbbVV/postprocessing/corrections.py index 2059a85d..5981d572 100644 --- a/src/HHbbVV/postprocessing/corrections.py +++ b/src/HHbbVV/postprocessing/corrections.py @@ -316,9 +316,13 @@ def get_lpsf( # fraction of subjets > 350 * 0.21 measured by CASE uncs["sj_pt_unc"] = (np.sum(events[f"{jet}_lp_sf_num_sjpt_gt350"][0]) / tot_matched) * 0.21 + # TODO: check double counting + if VV: num_prongs = events["ak8FatJetHVVNumProngs"][0] + # TODO: re-write as an OR over double matched, unmatched, and near jet boundary + sj_matching_unc = np.sum(events[f"{jet}_lp_sf_double_matched_event"][0]) for nump in range(2, 5): sj_matching_unc += ( diff --git a/src/HHbbVV/processors/corrections.py b/src/HHbbVV/processors/corrections.py index 1fd3cbeb..07bdc927 100644 --- a/src/HHbbVV/processors/corrections.py +++ b/src/HHbbVV/processors/corrections.py @@ -624,7 +624,7 @@ def add_trig_effs(weights: Weights, fatjets: FatJetArray, year: str, num_jets: i # ------------------- Lund plane reweighting ------------------- # -MAX_PT_FPARAMS = 4 # max order (+1) of pt extrapolation functions +MAX_PT_FPARAMS = 3 # max order (+1) of pt extrapolation functions MAX_PT_BIN = 350 # have to use subjet pt extrapolation for subjet pT > this ( @@ -647,7 +647,7 @@ def _get_lund_lookups(seed: int = 42, lnN: bool = True, trunc_gauss: bool = Fals import uproot # initialize lund plane scale factors lookups - f = uproot.open(package_path + "/corrections/lp_ratio_jan20.root") + f = uproot.open(package_path + "/corrections/ratio_june9.root") # 3D histogram: [subjet_pt, ln(0.8/Delta), ln(kT/GeV)] ratio_nom = f["ratio_nom"].to_numpy() @@ -743,7 +743,12 @@ def _np_pad(arr: np.ndarray, target: int = MAX_PT_FPARAMS): ) -def _get_lund_arrays(events: NanoEventsArray, fatjet_idx: Union[int, ak.Array], num_prongs: int): +def _get_lund_arrays( + events: NanoEventsArray, + jec_fatjets: FatJetArray, + fatjet_idx: Union[int, ak.Array], + num_prongs: int, +): """ Gets the ``num_prongs`` subjet pTs and Delta and kT per primary LP splitting of fatjets at ``fatjet_idx`` in each event. @@ -753,6 +758,7 @@ def _get_lund_arrays(events: NanoEventsArray, fatjet_idx: Union[int, ak.Array], Args: events (NanoEventsArray): nano events + jec_fatjets (FatJetArray): post-JEC fatjets, used to update subjet pTs. fatjet_idx (int | ak.Array): fatjet index num_prongs (int): number of prongs / subjets per jet to reweight @@ -850,7 +856,7 @@ def _calc_lund_SFs( high_pt_sel = flat_subjet_pt > max_pt_bin hpt_logD = flat_logD[high_pt_sel] hpt_logkt = flat_logkt[high_pt_sel] - hpt_sjpt = flat_subjet_pt[high_pt_sel] # change this to 1/sjpt for next iteration + hpt_sjpt = 1 / flat_subjet_pt[high_pt_sel] # store polynomial orders for pT extrapolation sj_pt_orders = np.array([np.power(hpt_sjpt, i) for i in range(max_fparams)]).T @@ -885,13 +891,13 @@ def _calc_lund_SFs( ) ) - return np.array( - sf_vals - ).T # output shape: ``[n_jets, len(ratio_lookups) x len(pt_extrap_lookups)]`` + # output shape: ``[n_jets, len(ratio_lookups) x len(pt_extrap_lookups)]`` + return np.array(sf_vals).T def get_lund_SFs( events: NanoEventsArray, + jec_fatjets: FatJetArray, fatjet_idx: Union[int, ak.Array], num_prongs: int, gen_quarks: GenParticleArray, @@ -906,6 +912,7 @@ def get_lund_SFs( Args: events (NanoEventsArray): nano events + jec_fatjets (FatJetArray): post-JEC fatjets, used to update subjet pTs. fatjet_idx (int | ak.Array): fatjet index num_prongs (int): number of prongs / subjets per jet to r seed (int, optional): seed for random smearings. Defaults to 42. @@ -930,8 +937,10 @@ def get_lund_SFs( pt_extrap_lookups_dict, ) = _get_lund_lookups(seed, lnN, trunc_gauss) + jec_fatjet = jec_fatjets[np.arange(len(jec_fatjet)), fatjet_idx] + flat_logD, flat_logkt, flat_subjet_pt, ld_offsets, kt_subjets_vec = _get_lund_arrays( - events, fatjet_idx, num_prongs + events, jec_fatjet, fatjet_idx, num_prongs ) sfs = {} @@ -1011,15 +1020,20 @@ def get_lund_SFs( sj_matched_idx_mask = np.copy(sj_matched_idx) sj_matched_idx_mask[~sj_matched] = -1 + j_q_dr = gen_quarks.delta_r(jec_fatjet) + q_boundary = (j_q_dr > 0.7) * (j_q_dr < 0.9) + # events with quarks at the boundary of the jet + sfs["lp_sf_boundary_quarks"] = np.any(q_boundary, axis=1, keepdims=True) + # events which have more than one quark matched to the same subjet sfs["lp_sf_double_matched_event"] = np.any( - [np.sum(sj_matched_idx_mask == i, axis=1) > 1 for i in range(3)], axis=0 + [np.sum(sj_matched_idx_mask == i, axis=1) > 1 for i in range(num_prongs)], axis=0 ).astype(int)[:, np.newaxis] # number of quarks per event which aren't matched sfs["lp_sf_unmatched_quarks"] = np.sum(~sj_matched, axis=1, keepdims=True) - # old pT extrapolation uncertainty + # OLD pT extrapolation uncertainty sfs["lp_sf_num_sjpt_gt350"] = np.sum(kt_subjets_vec.pt > 350, axis=1, keepdims=True).to_numpy() return sfs diff --git a/src/HHbbVV/scale_factors/corrections.py b/src/HHbbVV/scale_factors/corrections.py index 54666bf0..fa004896 100644 --- a/src/HHbbVV/scale_factors/corrections.py +++ b/src/HHbbVV/scale_factors/corrections.py @@ -385,7 +385,7 @@ def add_top_pt_weight(weights: Weights, events: NanoEventsArray): # find corrections path using this file's path try: - with gzip.open(package_path + "/data/jec_compiled.pkl.gz", "rb") as filehandler: + with open(package_path + "/corrections/jec_compiled.pkl", "rb") as filehandler: jmestuff = pickle.load(filehandler) fatjet_factory = jmestuff["fatjet_factory"] diff --git a/src/HHbbVV/scale_factors/top_reweighting.ipynb b/src/HHbbVV/scale_factors/top_reweighting.ipynb index 88c98ef4..48b2ac44 100644 --- a/src/HHbbVV/scale_factors/top_reweighting.ipynb +++ b/src/HHbbVV/scale_factors/top_reweighting.ipynb @@ -2,14 +2,14 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 61, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Failed loading compiled JECs\n" + "/uscms_data/d3/rkansal/HHbbVV/src/HHbbVV\n" ] } ], @@ -78,27 +78,27 @@ "name": "stderr", "output_type": "stream", "text": [ - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15SubJet_nBHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15SubJet_nBHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15SubJet_nCHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15SubJet_nCHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15_nBHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15_nBHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15_nCHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJetAK15_nCHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_btagDDBvLV2 in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_btagDDBvLV2 in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_btagDDCvBV2 in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_btagDDCvBV2 in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_btagDDCvLV2 in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_btagDDCvLV2 in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_nBHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_nBHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_nCHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch FatJet_nCHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch SubJet_nBHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch SubJet_nBHadrons in , taking first instance\n", " warnings.warn(\n", - "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch SubJet_nCHadrons in , taking first instance\n", + "/uscms_data/d3/rkansal/mambaforge/envs/python39/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:75: UserWarning: Found duplicate branch SubJet_nCHadrons in , taking first instance\n", " warnings.warn(\n" ] } @@ -258,21 +258,9 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "NameError", - "evalue": "name 'metfilters' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m/uscms/home/rkansal/nobackup/HHbbVV/src/HHbbVV/scale_factors/top_reweighting.ipynb Cell 10\u001b[0m in \u001b[0;36m1\n\u001b[1;32m 3\u001b[0m met_selection \u001b[39m=\u001b[39m met\u001b[39m.\u001b[39mpt \u001b[39m>\u001b[39m\u001b[39m=\u001b[39m met_selection[\u001b[39m\"\u001b[39m\u001b[39mpt\u001b[39m\u001b[39m\"\u001b[39m]\n\u001b[1;32m 5\u001b[0m \u001b[39m# metfilters = np.ones(len(events), dtype=\"bool\")\u001b[39;00m\n\u001b[1;32m 6\u001b[0m \u001b[39m# metfilterkey = \"data\" if isData else \"mc\"\u001b[39;00m\n\u001b[1;32m 7\u001b[0m \u001b[39m# for mf in self.metfilters[year][metfilterkey]:\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[39m# if mf in events.Flag.fields:\u001b[39;00m\n\u001b[1;32m 9\u001b[0m \u001b[39m# metfilters = metfilters & events.Flag[mf]\u001b[39;00m\n\u001b[0;32m---> 11\u001b[0m add_selection(\u001b[39m\"\u001b[39m\u001b[39mmet\u001b[39m\u001b[39m\"\u001b[39m, met_selection \u001b[39m*\u001b[39m metfilters, \u001b[39m*\u001b[39mselection_args)\n\u001b[1;32m 13\u001b[0m \u001b[39m# leptonic W selection\u001b[39;00m\n\u001b[1;32m 14\u001b[0m \u001b[39m# add_selection(\"lepW\", (met + muon).pt >= self.lepW_selection[\"pt\"], *selection_args)\u001b[39;00m\n\u001b[1;32m 15\u001b[0m add_selection(\u001b[39m\"\u001b[39m\u001b[39mlepW\u001b[39m\u001b[39m\"\u001b[39m, met\u001b[39m.\u001b[39mpt \u001b[39m+\u001b[39m muon\u001b[39m.\u001b[39mpt \u001b[39m>\u001b[39m\u001b[39m=\u001b[39m \u001b[39m100\u001b[39m, \u001b[39m*\u001b[39mselection_args)\n", - "\u001b[0;31mNameError\u001b[0m: name 'metfilters' is not defined" - ] - } - ], + "outputs": [], "source": [ "# met\n", "met = events.MET\n", @@ -492,12 +480,12 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 79, "metadata": {}, "outputs": [], "source": [ - "fatjets = events.FatJet # CHANGE TO JEC JETS\n", - "# fatjets = corrections.get_jec_jets(events, \"2018\")\n", + "# fatjets = events.FatJet # CHANGE TO JEC JETS\n", + "fatjets = corrections.get_jec_jets(events, \"2018\")\n", "\n", "num_jets = 1\n", "\n", @@ -517,7 +505,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 80, "metadata": {}, "outputs": [], "source": [ @@ -532,7 +520,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 81, "metadata": {}, "outputs": [], "source": [ @@ -555,7 +543,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 82, "metadata": {}, "outputs": [], "source": [ @@ -602,7 +590,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 84, "metadata": {}, "outputs": [], "source": [ @@ -647,7 +635,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 85, "metadata": {}, "outputs": [], "source": [ @@ -656,7 +644,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 86, "metadata": {}, "outputs": [], "source": [ @@ -675,7 +663,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 87, "metadata": {}, "outputs": [], "source": [ @@ -690,7 +678,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 88, "metadata": {}, "outputs": [], "source": [ @@ -711,7 +699,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 89, "metadata": {}, "outputs": [], "source": [ @@ -725,6 +713,42 @@ "merged_pfcands = merged_top_events.PFCands[merged_ak8_pfcands.pFCandsIdx]" ] }, + { + "cell_type": "code", + "execution_count": 91, + "metadata": {}, + "outputs": [], + "source": [ + "jec_fatjets = fatjets[select]\n", + "jec_correction = (\n", + " jec_fatjets.pt[merged_top_jet_match][fatjet_idx[merged_top_jet_match]]\n", + " / presel_events.FatJet.pt[merged_top_jet_match][fatjet_idx[merged_top_jet_match]]\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 110, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 110, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fatjet = ak.pad_none(jec_fatjets[merged_top_jet_match], 2, axis=1)[\n", + " np.arange(np.sum(merged_top_jet_match)), np.array(fatjet_idx[merged_top_jet_match])\n", + "]\n", + "fatjet" + ] + }, { "cell_type": "code", "execution_count": 18, @@ -816,7 +840,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Subjet matching and pT extrapolation uncertainty" + "### Subjet matching uncertainty" ] }, { @@ -843,20 +867,45 @@ "cell_type": "code", "execution_count": 24, "metadata": {}, + "outputs": [], + "source": [ + "sj_matched = []\n", + "sj_matched_idx = []\n", + "\n", + "for i in range(num_prongs):\n", + " sj_q_dr = kt_subjets_vec.delta_r(gen_quarks[:, i])\n", + " sj_matched.append(ak.min(sj_q_dr, axis=1) <= matching_dR)\n", + " sj_matched_idx.append(ak.argmin(sj_q_dr, axis=1))" + ] + }, + { + "cell_type": "code", + "execution_count": 116, + "metadata": {}, + "outputs": [], + "source": [ + "j_q_dr = gen_quarks.delta_r(fatjet)\n", + "q_boundary = (j_q_dr > 0.7) * (j_q_dr < 0.9)" + ] + }, + { + "cell_type": "code", + "execution_count": 119, + "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "21" ] }, - "execution_count": 24, + "execution_count": 119, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "gen_quarks" + "np.sum(np.any(q_boundary, axis=1))" ] }, { @@ -864,21 +913,6 @@ "execution_count": 25, "metadata": {}, "outputs": [], - "source": [ - "sj_matched = []\n", - "sj_matched_idx = []\n", - "\n", - "for i in range(num_prongs):\n", - " sj_q_dr = kt_subjets_vec.delta_r(gen_quarks[:, i])\n", - " sj_matched.append(ak.min(sj_q_dr, axis=1) <= matching_dR)\n", - " sj_matched_idx.append(ak.argmin(sj_q_dr, axis=1))" - ] - }, - { - "cell_type": "code", - "execution_count": 26, - "metadata": {}, - "outputs": [], "source": [ "sj_matched = np.array(sj_matched).T\n", "sj_matched_idx = np.array(sj_matched_idx).T\n", @@ -901,16 +935,16 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 53, "metadata": {}, "outputs": [], "source": [ - "f = uproot.open(\"../corrections/ratio_feb20.root\")\n", + "f = uproot.open(\"../corrections/ratio_june9.root\")\n", "\n", "# 3D histogram: [subjet_pt, ln(0.8/Delta), ln(kT/GeV)]\n", "ratio_nom = f[\"ratio_nom\"].to_numpy()\n", "ratio_nom_errs = f[\"ratio_nom\"].errors()\n", - "ratio_nom_edges = ratio_nom[1:]\n", + "ratio_edges = ratio_nom[1:]\n", "ratio_nom = ratio_nom[0]\n", "\n", "zero_vals = ratio_nom == 0\n", @@ -922,11 +956,11 @@ }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 125, "metadata": {}, "outputs": [], "source": [ - "max_pt_bin = ratio_nom_edges[0][-2]" + "bl_ratio = f[\"h_bl_ratio\"].to_numpy()[0]" ] }, { @@ -936,35 +970,43 @@ "outputs": [], "source": [ "plt.figure(figsize=(10, 6))\n", - "plt.imshow(ratio_nom[5].T[::-1], extent=[-1, 8, -5, 6.9], cmap=\"turbo\")\n", + "plt.imshow(bl_ratio[5].T[::-1], extent=[-1, 8, -5, 6.9], cmap=\"viridis\")\n", "plt.xlabel(r\"ln$(0.8/\\Delta)$\")\n", "plt.ylabel(r\"ln$(k_T/GeV)$\")\n", "_ = plt.colorbar()" ] }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [], + "source": [ + "max_pt_bin = ratio_nom_edges[0][-2]\n", + "n_LP_sf_toys = 100" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# maxp = 0\n", - "# for i in range(ratio_shape[1]):\n", - "# for j in range(ratio_shape[2]):\n", - "# if len(pt_extrap_params[i][j]) > maxp:\n", - "# maxp = len(pt_extrap_params[i][j])\n", - "\n", - "# maxp" + "plt.figure(figsize=(10, 6))\n", + "plt.imshow(ratio_nom[5].T[::-1], extent=[-1, 8, -5, 6.9], cmap=\"turbo\")\n", + "plt.xlabel(r\"ln$(0.8/\\Delta)$\")\n", + "plt.ylabel(r\"ln$(k_T/GeV)$\")\n", + "_ = plt.colorbar()" ] }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 54, "metadata": {}, "outputs": [], "source": [ - "pt_extrap_lookups = {\"params\": [], \"errs\": [], \"sys_up_params\": [], \"sys_down_params\": []}\n", - "max_fparams = 4 # max order (+1) of functions\n", + "pt_extrap_lookups_dict = {\"params\": [], \"errs\": [], \"sys_up_params\": [], \"sys_down_params\": []}\n", + "max_fparams = 3 # max order (+1) of functions\n", "\n", "\n", "def _np_pad(arr: np.ndarray, target: int = max_fparams):\n", @@ -972,23 +1014,23 @@ "\n", "\n", "for i in range(ratio_shape[1]):\n", - " for key in pt_extrap_lookups:\n", - " pt_extrap_lookups[key].append([])\n", + " for key in pt_extrap_lookups_dict:\n", + " pt_extrap_lookups_dict[key].append([])\n", "\n", " for j in range(ratio_shape[2]):\n", " func = f[\"pt_extrap\"][f\"func_{i + 1}_{j + 1}\"]\n", - " pt_extrap_lookups[\"params\"][-1].append(\n", + " pt_extrap_lookups_dict[\"params\"][-1].append(\n", " _np_pad(func._members[\"fFormula\"]._members[\"fClingParameters\"])\n", " )\n", - " pt_extrap_lookups[\"errs\"][-1].append(_np_pad(func._members[\"fParErrors\"]))\n", - " pt_extrap_lookups[\"sys_up_params\"][-1].append(\n", + " pt_extrap_lookups_dict[\"errs\"][-1].append(_np_pad(func._members[\"fParErrors\"]))\n", + " pt_extrap_lookups_dict[\"sys_up_params\"][-1].append(\n", " _np_pad(\n", " f[\"pt_extrap\"][f\"func_sys_tot_up_{i + 1}_{j + 1}\"]\n", " ._members[\"fFormula\"]\n", " ._members[\"fClingParameters\"]\n", " )\n", " )\n", - " pt_extrap_lookups[\"sys_down_params\"][-1].append(\n", + " pt_extrap_lookups_dict[\"sys_down_params\"][-1].append(\n", " _np_pad(\n", " f[\"pt_extrap\"][f\"func_sys_tot_down_{i + 1}_{j + 1}\"]\n", " ._members[\"fFormula\"]\n", @@ -996,13 +1038,38 @@ " )\n", " )\n", "\n", - "for key in pt_extrap_lookups:\n", - " pt_extrap_lookups[key] = dense_lookup(np.array(pt_extrap_lookups[key]), ratio_nom_edges[1:])" + "for key in pt_extrap_lookups_dict:\n", + " pt_extrap_lookups_dict[key] = np.array(pt_extrap_lookups_dict[key])\n", + "\n", + "for arr in pt_extrap_lookups_dict.values():\n", + " print(np.max(np.sum(arr != 0, axis=2)))\n", + "\n", + "# smear parameters according to errors for pt extrap unc.\n", + "rand_noise = np.random.normal(size=[n_LP_sf_toys, *pt_extrap_lookups_dict[\"params\"].shape])\n", + "smeared_pt_params = pt_extrap_lookups_dict[\"params\"] + (pt_extrap_lookups_dict[\"errs\"] * rand_noise)\n", + "\n", + "for key in pt_extrap_lookups_dict:\n", + " pt_extrap_lookups_dict[key] = dense_lookup(pt_extrap_lookups_dict[key], ratio_edges[1:])\n", + "\n", + "pt_extrap_lookups_dict[\"smeared_params\"] = [\n", + " dense_lookup(smeared_pt_params[i], ratio_edges[1:]) for i in range(n_LP_sf_toys)\n", + "]" ] }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 39, + "metadata": {}, + "outputs": [], + "source": [ + "n_sf_toys = 100\n", + "np.random.seed(42)\n", + "rand_noise = np.random.normal(size=[n_sf_toys, *pt_extrap_lookups[\"params\"].shape])" + ] + }, + { + "cell_type": "code", + "execution_count": 42, "metadata": {}, "outputs": [], "source": [ @@ -1019,7 +1086,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 43, "metadata": {}, "outputs": [], "source": [ @@ -1033,36 +1100,53 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "high_pt_sel = flat_subjet_pt > max_pt_bin\n", "hpt_logD = flat_logD[high_pt_sel]\n", "hpt_logkt = flat_logkt[high_pt_sel]\n", - "hpt_sjpt = flat_subjet_pt[high_pt_sel] # change this to 1/sjpt for next iteration" + "hpt_sjpt = 1 / flat_subjet_pt[high_pt_sel]\n", + "\n", + "# store polynomial orders for pT extrapolation\n", + "sj_pt_orders = np.array([np.power(hpt_sjpt, i) for i in range(max_fparams)]).T" ] }, { "cell_type": "code", - "execution_count": 59, + "execution_count": 57, "metadata": {}, "outputs": [], "source": [ - "# store polynomial orders for pT extrapolation\n", - "sj_pt_orders = np.array([np.power(hpt_sjpt, i) for i in range(max_fparams)]).T" + "clip_max, clip_min = 10, 0.1\n", + "pt_lookup = pt_extrap_lookups_dict[\"params\"]\n", + "params = pt_lookup(hpt_logD, hpt_logkt)\n", + "pt_extrap_vals = np.maximum(np.minimum(np.sum(params * sj_pt_orders, axis=1), clip_max), clip_min)" ] }, { "cell_type": "code", - "execution_count": 78, + "execution_count": 58, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAh7klEQVR4nO3df1TUVeL/8Re/GZEB1AAxNPyxaompFYS10UlO2lprR8+m5bZaHm2LLMW1ZDd1xcofmXFyScujaOdklp3U2tRaOWlroRZpmZpJWZIuuKUyggoK9/vHfp1PE6gMDpcfPh/nzDnxnjuXe3kz8GycGfyMMUYAAACW+Df2AgAAwOWF+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVgY29gF+rrq7W4cOHFR4eLj8/v8ZeDgAAqANjjE6cOKG4uDj5+1/4sY0mFx+HDx9WfHx8Yy8DAADUQ1FRka688soLjmly8REeHi7pf4t3Op2NvBoAAFAXLpdL8fHx7t/jF9Lk4uPcP7U4nU7iAwCAZqYuT5ngCacAAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArGpyf9UWQNNw6PgpHSuv9Pm8UWHB6hDp8Pm8AJoP4gNADYeOn1La85t16kyVz+d2BAVo46RUAgS4jBEfAGo4Vl6pU2eqlD28j7pGt/bZvIVHyjThjZ06Vl5JfACXMeIDwHl1jW6tXh0iGnsZAFoYnnAKAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVXsVHVVWVpk6dqoSEBDkcDnXp0kUzZ86UMcY9xhijadOmqX379nI4HEpLS9P+/ft9vnAAANA8eRUfc+bM0cKFC/WPf/xDe/fu1Zw5czR37lwtWLDAPWbu3Ll68cUXtWjRIm3btk1hYWEaOHCgTp8+7fPFAwCA5ifQm8GffPKJhgwZosGDB0uSrrrqKr3++uvavn27pP896pGdna2nnnpKQ4YMkSS9+uqriomJ0Zo1azRixAgfLx8AADQ3Xj3y0b9/f+Xl5embb76RJH3xxRfasmWL7rjjDknSgQMHVFxcrLS0NPdtIiIilJycrPz8/FrnrKiokMvl8rgAAICWy6tHPqZMmSKXy6UePXooICBAVVVVeuaZZzRy5EhJUnFxsSQpJibG43YxMTHu635t1qxZmjFjRn3WDgAAmiGvHvl488039dprr2nFihX6/PPPtXz5cs2bN0/Lly+v9wIyMzNVWlrqvhQVFdV7LgAA0PR59cjH5MmTNWXKFPdzNxITE/XDDz9o1qxZGjVqlGJjYyVJJSUlat++vft2JSUl6tOnT61zhoSEKCQkpJ7LBwAAzY1Xj3ycPHlS/v6eNwkICFB1dbUkKSEhQbGxscrLy3Nf73K5tG3bNqWkpPhguQAAoLnz6pGPu+66S88884w6duyoa665Rjt27ND8+fP14IMPSpL8/Pw0YcIEPf300+rWrZsSEhI0depUxcXF6e67726I9QMAgGbGq/hYsGCBpk6dqkceeURHjhxRXFycHnroIU2bNs095oknnlB5ebnGjRun48eP6+abb9aGDRsUGhrq88UDAIDmx6v4CA8PV3Z2trKzs887xs/PT1lZWcrKyrrUtQEAgBaIv+0CAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACs8jo+Dh06pD/+8Y9q27atHA6HEhMT9dlnn7mvN8Zo2rRpat++vRwOh9LS0rR//36fLhoAADRfXsXHsWPHdNNNNykoKEjr16/Xnj179PzzzysqKso9Zu7cuXrxxRe1aNEibdu2TWFhYRo4cKBOnz7t88UDAIDmJ9CbwXPmzFF8fLxyc3PdxxISEtz/bYxRdna2nnrqKQ0ZMkSS9OqrryomJkZr1qzRiBEjfLRsAADQXHn1yMc777yj66+/Xn/4wx8UHR2tvn37avHixe7rDxw4oOLiYqWlpbmPRUREKDk5Wfn5+bXOWVFRIZfL5XEBAAAtl1fx8d1332nhwoXq1q2b3n//fT388MN67LHHtHz5cklScXGxJCkmJsbjdjExMe7rfm3WrFmKiIhwX+Lj4+uzDwAA0Ex4FR/V1dXq16+fnn32WfXt21fjxo3T2LFjtWjRonovIDMzU6Wlpe5LUVFRvecCAABNn1fx0b59e1199dUex3r27KmDBw9KkmJjYyVJJSUlHmNKSkrc1/1aSEiInE6nxwUAALRcXsXHTTfdpH379nkc++abb9SpUydJ/3vyaWxsrPLy8tzXu1wubdu2TSkpKT5YLgAAaO68erXLxIkT1b9/fz377LO65557tH37dr3yyit65ZVXJEl+fn6aMGGCnn76aXXr1k0JCQmaOnWq4uLidPfddzfE+gEAQDPjVXzccMMNWr16tTIzM5WVlaWEhARlZ2dr5MiR7jFPPPGEysvLNW7cOB0/flw333yzNmzYoNDQUJ8vHgAAND9exYck3XnnnbrzzjvPe72fn5+ysrKUlZV1SQsDAAAtE3/bBQAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALCK+AAAAFYRHwAAwCriAwAAWEV8AAAAq4gPAABgFfEBAACsIj4AAIBVxAcAALAqsLEXAODyU3ikrEHmjQoLVodIR4PMDcB3iA8A1kSFBcsRFKAJb+xskPkdQQHaOCmVAAGaOOIDgDUdIh3aOClVx8orfT534ZEyTXhjp46VVxIfQBNHfACwqkOkgzgALnM84RQAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWXVJ8zJ49W35+fpowYYL72OnTp5Wenq62bduqdevWGjZsmEpKSi51nQAAoIWod3x8+umnevnll9W7d2+P4xMnTtS7776rVatWafPmzTp8+LCGDh16yQsFAAAtQ73io6ysTCNHjtTixYsVFRXlPl5aWqolS5Zo/vz5uu2223TdddcpNzdXn3zyibZu3eqzRQMAgOarXvGRnp6uwYMHKy0tzeN4QUGBzpw543G8R48e6tixo/Lz82udq6KiQi6Xy+MCAABarkBvb7By5Up9/vnn+vTTT2tcV1xcrODgYEVGRnocj4mJUXFxca3zzZo1SzNmzPB2GQAAoJny6pGPoqIiPf7443rttdcUGhrqkwVkZmaqtLTUfSkqKvLJvAAAoGnyKj4KCgp05MgR9evXT4GBgQoMDNTmzZv14osvKjAwUDExMaqsrNTx48c9bldSUqLY2Nha5wwJCZHT6fS4AACAlsurf3YZMGCAdu3a5XHsgQceUI8ePfTkk08qPj5eQUFBysvL07BhwyRJ+/bt08GDB5WSkuK7VQMAgGbLq/gIDw9Xr169PI6FhYWpbdu27uNjxoxRRkaG2rRpI6fTqfHjxyslJUU33nij71YNAACaLa+fcHoxL7zwgvz9/TVs2DBVVFRo4MCBeumll3z9aQAAQDN1yfGxadMmj49DQ0OVk5OjnJycS50aAAC0QPxtFwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFgV2NgLAFB/h46f0rHySp/PW3ikzOdzAsA5xAfQTB06fkppz2/WqTNVDTK/IyhAUWHBDTI3gMsb8QE0U8fKK3XqTJWyh/dR1+jWPp8/KixYHSIdPp8XAIgPoJnrGt1avTpENPYyAKDOeMIpAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKu8io9Zs2bphhtuUHh4uKKjo3X33Xdr3759HmNOnz6t9PR0tW3bVq1bt9awYcNUUlLi00UDAIDmy6v42Lx5s9LT07V161b961//0pkzZ3T77bervLzcPWbixIl69913tWrVKm3evFmHDx/W0KFDfb5wAADQPAV6M3jDhg0eHy9btkzR0dEqKCjQLbfcotLSUi1ZskQrVqzQbbfdJknKzc1Vz549tXXrVt14442+WzkAAGiWLuk5H6WlpZKkNm3aSJIKCgp05swZpaWlucf06NFDHTt2VH5+/qV8KgAA0EJ49cjHL1VXV2vChAm66aab1KtXL0lScXGxgoODFRkZ6TE2JiZGxcXFtc5TUVGhiooK98cul6u+SwIAFR4p8/mcUWHB6hDp8Pm8wOWq3vGRnp6ur776Slu2bLmkBcyaNUszZsy4pDkAICosWI6gAE14Y6fP53YEBWjjpFQCBPCResXHo48+qn/+85/66KOPdOWVV7qPx8bGqrKyUsePH/d49KOkpESxsbG1zpWZmamMjAz3xy6XS/Hx8fVZFoDLWIdIhzZOStWx8kqfzlt4pEwT3tipY+WVxAfgI17FhzFG48eP1+rVq7Vp0yYlJCR4XH/dddcpKChIeXl5GjZsmCRp3759OnjwoFJSUmqdMyQkRCEhIfVcPgD8nw6RDgIBaAa8io/09HStWLFCa9euVXh4uPt5HBEREXI4HIqIiNCYMWOUkZGhNm3ayOl0avz48UpJSeGVLgAAQJKX8bFw4UJJ0q233upxPDc3V6NHj5YkvfDCC/L399ewYcNUUVGhgQMH6qWXXvLJYgEAQPPn9T+7XExoaKhycnKUk5NT70UBAICWi7/tAgAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMCqwMZeAHA5OHT8lI6VV/p0zsIjZT6dDxfWUF/vqLBgdYh0NMjcQFNFfAAN7NDxU0p7frNOnany+dyOoABFhQX7fF78n6iwYDmCAjThjZ0NMr8jKEAbJ6USILisEB9AAztWXqlTZ6qUPbyPuka39unc/F9zw+sQ6dDGSak+f+RK+t+jKRPe2Klj5ZWcR1xWiA/Akq7RrdWrQ0RjLwP10CHSQRwAPsQTTgEAgFXEBwAAsIr4AAAAVhEfAADAKuIDAABYRXwAAACriA8AAGDVZfc+Hw3xNtcSb/YEAEBdNVh85OTk6LnnnlNxcbGuvfZaLViwQElJSQ316eqkod/mmrdIBgDg4hokPt544w1lZGRo0aJFSk5OVnZ2tgYOHKh9+/YpOjq6IT5lnTTU21zzFskAANRdg8TH/PnzNXbsWD3wwAOSpEWLFum9997T0qVLNWXKlIb4lF7hba4BAGg8Po+PyspKFRQUKDMz033M399faWlpys/PrzG+oqJCFRUV7o9LS0slSS6Xy9dLU9kJl6orTqrshEsul5/P5/3yu/+o7ITv143m7bv/ljfI9x2aP352oLFc0TpEVzhDfTrnud/bxpiLDzY+dujQISPJfPLJJx7HJ0+ebJKSkmqMnz59upHEhQsXLly4cGkBl6Kioou2QqO/2iUzM1MZGRnuj6urq3X06FG1bdtWfn6+/b9El8ul+Ph4FRUVyel0+nTupqCl709q+Xtkf81fS98j+2v+GmqPxhidOHFCcXFxFx3r8/ho166dAgICVFJS4nG8pKREsbGxNcaHhIQoJCTE41hkZKSvl+XB6XS22G8qqeXvT2r5e2R/zV9L3yP7a/4aYo8RERF1GufzNxkLDg7Wddddp7y8PPex6upq5eXlKSUlxdefDgAANDMN8s8uGRkZGjVqlK6//nolJSUpOztb5eXl7le/AACAy1eDxMfw4cP13//+V9OmTVNxcbH69OmjDRs2KCYmpiE+XZ2FhIRo+vTpNf6Zp6Vo6fuTWv4e2V/z19L3yP6av6awRz9j6vKaGAAAAN/gD8sBAACriA8AAGAV8QEAAKwiPgAAgFUtKj6eeeYZ9e/fX61atarzG5UZYzRt2jS1b99eDodDaWlp2r9/v8eYo0ePauTIkXI6nYqMjNSYMWNUVlbWADu4OG/X8v3338vPz6/Wy6pVq9zjart+5cqVNrbkoT5f61tvvbXG2v/85z97jDl48KAGDx6sVq1aKTo6WpMnT9bZs2cbciu18nZ/R48e1fjx49W9e3c5HA517NhRjz32mPtvIJ3TmOcvJydHV111lUJDQ5WcnKzt27dfcPyqVavUo0cPhYaGKjExUevWrfO4vi73SZu82d/ixYv129/+VlFRUYqKilJaWlqN8aNHj65xrgYNGtTQ27ggb/a4bNmyGusPDfX8GyHN+RzW9vPEz89PgwcPdo9pSufwo48+0l133aW4uDj5+flpzZo1F73Npk2b1K9fP4WEhKhr165atmxZjTHe3q+95oM/59JkTJs2zcyfP99kZGSYiIiIOt1m9uzZJiIiwqxZs8Z88cUX5ve//71JSEgwp06dco8ZNGiQufbaa83WrVvNv//9b9O1a1dz7733NtAuLszbtZw9e9b85z//8bjMmDHDtG7d2pw4ccI9TpLJzc31GPfLr4Et9flap6ammrFjx3qsvbS01H392bNnTa9evUxaWprZsWOHWbdunWnXrp3JzMxs6O3U4O3+du3aZYYOHWreeecdU1hYaPLy8ky3bt3MsGHDPMY11vlbuXKlCQ4ONkuXLjW7d+82Y8eONZGRkaakpKTW8R9//LEJCAgwc+fONXv27DFPPfWUCQoKMrt27XKPqct90hZv93ffffeZnJwcs2PHDrN3714zevRoExERYX788Uf3mFGjRplBgwZ5nKujR4/a2lIN3u4xNzfXOJ1Oj/UXFxd7jGnO5/Dnn3/22NtXX31lAgICTG5urntMUzqH69atM3/729/M22+/bSSZ1atXX3D8d999Z1q1amUyMjLMnj17zIIFC0xAQIDZsGGDe4y3X7P6aFHxcU5ubm6d4qO6utrExsaa5557zn3s+PHjJiQkxLz++uvGGGP27NljJJlPP/3UPWb9+vXGz8/PHDp0yOdrvxBfraVPnz7mwQcf9DhWl2/ahlbf/aWmpprHH3/8vNevW7fO+Pv7e/yAXLhwoXE6naaiosIna68LX52/N9980wQHB5szZ864jzXW+UtKSjLp6enuj6uqqkxcXJyZNWtWrePvueceM3jwYI9jycnJ5qGHHjLG1O0+aZO3+/u1s2fPmvDwcLN8+XL3sVGjRpkhQ4b4eqn15u0eL/bztaWdwxdeeMGEh4ebsrIy97Gmdg7PqcvPgSeeeMJcc801HseGDx9uBg4c6P74Ur9mddGi/tnFWwcOHFBxcbHS0tLcxyIiIpScnKz8/HxJUn5+viIjI3X99de7x6Slpcnf31/btm2zul5frKWgoEA7d+7UmDFjalyXnp6udu3aKSkpSUuXLq3bn0X2oUvZ32uvvaZ27dqpV69eyszM1MmTJz3mTUxM9HiTu4EDB8rlcmn37t2+38h5+Op7qbS0VE6nU4GBnu8RaPv8VVZWqqCgwOP+4+/vr7S0NPf959fy8/M9xkv/OxfnxtflPmlLffb3aydPntSZM2fUpk0bj+ObNm1SdHS0unfvrocfflg///yzT9deV/XdY1lZmTp16qT4+HgNGTLE437U0s7hkiVLNGLECIWFhXkcbyrn0FsXuw/64mtWF43+V20bU3FxsSTVeOfVmJgY93XFxcWKjo72uD4wMFBt2rRxj7HFF2tZsmSJevbsqf79+3scz8rK0m233aZWrVrpgw8+0COPPKKysjI99thjPlv/xdR3f/fdd586deqkuLg4ffnll3ryySe1b98+vf322+55azvH566zxRfn76efftLMmTM1btw4j+ONcf5++uknVVVV1fq1/frrr2u9zfnOxS/vb+eOnW+MLfXZ3689+eSTiouL8/hBPmjQIA0dOlQJCQn69ttv9de//lV33HGH8vPzFRAQ4NM9XEx99ti9e3ctXbpUvXv3VmlpqebNm6f+/ftr9+7duvLKK1vUOdy+fbu++uorLVmyxON4UzqH3jrffdDlcunUqVM6duzYJX/f10WTj48pU6Zozpw5Fxyzd+9e9ejRw9KKfK+ue7xUp06d0ooVKzR16tQa1/3yWN++fVVeXq7nnnvOJ7+8Gnp/v/xFnJiYqPbt22vAgAH69ttv1aVLl3rPW1e2zp/L5dLgwYN19dVX6+9//7vHdQ15/lA/s2fP1sqVK7Vp0yaPJ2SOGDHC/d+JiYnq3bu3unTpok2bNmnAgAGNsVSvpKSkePyR0P79+6tnz556+eWXNXPmzEZcme8tWbJEiYmJSkpK8jje3M9hU9Dk42PSpEkaPXr0Bcd07ty5XnPHxsZKkkpKStS+fXv38ZKSEvXp08c95siRIx63O3v2rI4ePeq+/aWq6x4vdS1vvfWWTp48qT/96U8XHZucnKyZM2eqoqLikt//39b+zklOTpYkFRYWqkuXLoqNja3xTO2SkhJJ8sk5tLG/EydOaNCgQQoPD9fq1asVFBR0wfG+PH/n065dOwUEBLi/lueUlJScdz+xsbEXHF+X+6Qt9dnfOfPmzdPs2bO1ceNG9e7d+4JjO3furHbt2qmwsND6L65L2eM5QUFB6tu3rwoLCyW1nHNYXl6ulStXKisr66KfpzHPobfOdx90Op1yOBwKCAi45O+JOvHZs0eaEG+fcDpv3jz3sdLS0lqfcPrZZ5+5x7z//vuN+oTT+q4lNTW1xqskzufpp582UVFR9V5rffjqa71lyxYjyXzxxRfGmP97wukvn6n98ssvG6fTaU6fPu27DVxEffdXWlpqbrzxRpOammrKy8vr9Llsnb+kpCTz6KOPuj+uqqoyHTp0uOATTu+8806PYykpKTWecHqh+6RN3u7PGGPmzJljnE6nyc/Pr9PnKCoqMn5+fmbt2rWXvN76qM8ef+ns2bOme/fuZuLEicaYlnEOjfnf75GQkBDz008/XfRzNPY5PEd1fMJpr169PI7de++9NZ5weinfE3Vaq89magJ++OEHs2PHDvdLSXfs2GF27Njh8ZLS7t27m7ffftv98ezZs01kZKRZu3at+fLLL82QIUNqfalt3759zbZt28yWLVtMt27dGvWlthday48//mi6d+9utm3b5nG7/fv3Gz8/P7N+/foac77zzjtm8eLFZteuXWb//v3mpZdeMq1atTLTpk1r8P38mrf7KywsNFlZWeazzz4zBw4cMGvXrjWdO3c2t9xyi/s2515qe/vtt5udO3eaDRs2mCuuuKLRXmrrzf5KS0tNcnKySUxMNIWFhR4v7Tt79qwxpnHP38qVK01ISIhZtmyZ2bNnjxk3bpyJjIx0v7Lo/vvvN1OmTHGP//jjj01gYKCZN2+e2bt3r5k+fXqtL7W92H3SFm/3N3v2bBMcHGzeeustj3N17mfQiRMnzF/+8heTn59vDhw4YDZu3Gj69etnunXrZjWEL2WPM2bMMO+//7759ttvTUFBgRkxYoQJDQ01u3fvdo9pzufwnJtvvtkMHz68xvGmdg5PnDjh/l0nycyfP9/s2LHD/PDDD8YYY6ZMmWLuv/9+9/hzL7WdPHmy2bt3r8nJyan1pbYX+pr5QouKj1GjRhlJNS4ffvihe4z+//shnFNdXW2mTp1qYmJiTEhIiBkwYIDZt2+fx7w///yzuffee03r1q2N0+k0DzzwgEfQ2HSxtRw4cKDGno0xJjMz08THx5uqqqoac65fv9706dPHtG7d2oSFhZlrr73WLFq0qNaxDc3b/R08eNDccsstpk2bNiYkJMR07drVTJ482eN9Powx5vvvvzd33HGHcTgcpl27dmbSpEkeL1W1xdv9ffjhh7V+T0syBw4cMMY0/vlbsGCB6dixowkODjZJSUlm69at7utSU1PNqFGjPMa/+eab5je/+Y0JDg4211xzjXnvvfc8rq/LfdImb/bXqVOnWs/V9OnTjTHGnDx50tx+++3miiuuMEFBQaZTp05m7NixPv2hXh/e7HHChAnusTExMeZ3v/ud+fzzzz3ma87n0Bhjvv76ayPJfPDBBzXmamrn8Hw/I87tadSoUSY1NbXGbfr06WOCg4NN586dPX4nnnOhr5kv+Blj+fWUAADgsnZZv88HAACwj/gAAABWER8AAMAq4gMAAFhFfAAAAKuIDwAAYBXxAQAArCI+AACAVcQHAACwivgAAABWER8AAMAq4gMAAFj1/wAz1cD3mj9elgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "clip_max, clip_min = 10, 0.1\n", - "pt_lookup = pt_extrap_lookups[\"params\"]\n", - "params = pt_lookup(hpt_logD, hpt_logkt)\n", - "pt_extrap_vals = np.maximum(np.minimum(np.sum(params * sj_pt_orders, axis=1), clip_max), clip_min)" + "_ = plt.hist(\n", + " np.log10(pt_extrap_vals),\n", + " np.linspace(-1, 1, 21),\n", + " histtype=\"step\",\n", + ")" ] }, { @@ -1082,8 +1166,9 @@ } ], "source": [ + "# old (non 1/pt) version\n", "_ = plt.hist(\n", - " np.log10(np.maximum(np.sum(params * sj_pt_orders, axis=1), 0.1)),\n", + " np.log10(pt_extrap_vals),\n", " np.linspace(-1, 1, 21),\n", " histtype=\"step\",\n", ")"