diff --git a/src/HHbbVV/hh_vars.py b/src/HHbbVV/hh_vars.py index ae1cb108..dc1414eb 100644 --- a/src/HHbbVV/hh_vars.py +++ b/src/HHbbVV/hh_vars.py @@ -259,7 +259,10 @@ ("lp_sf_pt_extrap_vars", 100), ("lp_sf_sys_down", 1), ("lp_sf_sys_up", 1), + ("lp_sf_np_down", 1), + ("lp_sf_np_up", 1), ("lp_sf_double_matched_event", 1), - ("lp_sf_boundary_quarks", 1), + ("lp_sf_inside_boundary_quarks", 1), + ("lp_sf_outside_boundary_quarks", 1), ("lp_sf_unmatched_quarks", 1), ] diff --git a/src/HHbbVV/processors/bbVVSkimmer.py b/src/HHbbVV/processors/bbVVSkimmer.py index f4916c83..7df17501 100644 --- a/src/HHbbVV/processors/bbVVSkimmer.py +++ b/src/HHbbVV/processors/bbVVSkimmer.py @@ -120,6 +120,7 @@ class bbVVSkimmer(SkimmerABC): jecs = hh_vars.jecs # only the branches necessary for templates and post processing + # IMPORTANT: Add Lund plane branches in hh_vars.py min_branches = [ # noqa: RUF012 "ak8FatJetPhi", "ak8FatJetEta", @@ -673,15 +674,16 @@ def process(self, events: ak.Array): # Lund plane SFs ################ + lp_hist = None + if isSignal and self._systematics and self._lp_sfs: - # (var, # columns) logging.info("Starting LP SFs and saving: " + str(hh_vars.lp_sf_vars)) if len(skimmed_events["weight"]): genbb = genbb[sel_all] genq = genq[sel_all] - sf_dicts = [] + sf_dicts, lp_hists = [], [] lp_num_jets = num_jets if self._save_all else 1 for i in range(lp_num_jets): @@ -713,7 +715,7 @@ def process(self, events: ak.Array): for key, (selector, gen_quarks, num_prongs) in selectors.items(): if np.sum(selector) > 0: - selected_sfs[key] = get_lund_SFs( + selected_sfs[key], lp_hist = get_lund_SFs( year, events[sel_all][selector], fatjets[sel_all][selector], @@ -724,10 +726,13 @@ def process(self, events: ak.Array): ), # giving HVV jet index if only doing LP SFs for HVV jet num_prongs, gen_quarks[selector], + weights_dict["weight"][sel_all][selector], trunc_gauss=False, lnN=True, ) + lp_hists.append(lp_hist) + sf_dict = {} # collect all the scale factors, fill in 1s for unmatched jets @@ -743,6 +748,7 @@ def process(self, events: ak.Array): sf_dicts.append(sf_dict) sf_dicts = concatenate_dicts(sf_dicts) + lp_hist = sum(lp_hists) else: logging.info("No signal events selected") @@ -790,7 +796,7 @@ def process(self, events: ak.Array): fname = events.behavior["__events_factory__"]._partition_key.replace("/", "_") + ".parquet" self.dump_table(pddf, fname) - return {year: {dataset: {"totals": totals_dict, "cutflow": cutflow}}} + return {year: {dataset: {"totals": totals_dict, "cutflow": cutflow, "lp_hist": lp_hist}}} def postprocess(self, accumulator): return accumulator diff --git a/src/HHbbVV/processors/corrections.py b/src/HHbbVV/processors/corrections.py index 02d4d1c4..9f196c52 100644 --- a/src/HHbbVV/processors/corrections.py +++ b/src/HHbbVV/processors/corrections.py @@ -16,6 +16,7 @@ import awkward as ak import correctionlib +import hist import numpy as np from coffea import util as cutil from coffea.analysis_tools import Weights @@ -672,7 +673,8 @@ def add_trig_effs(weights: Weights, fatjets: FatJetArray, year: str, num_jets: i ratio_sys_down, pt_extrap_lookups_dict, bratio, -) = (None, None, None, None, None, None, None) + ratio_edges, +) = (None, None, None, None, None, None, None, None) def _get_lund_lookups(year: str, seed: int = 42, lnN: bool = True, trunc_gauss: bool = False): @@ -785,16 +787,17 @@ def _np_pad(arr: np.ndarray, target: int = MAX_PT_FPARAMS): ratio_sys_down, pt_extrap_lookups_dict, bratio, + ratio_edges, ) def _get_flat_lp_vars(lds, kt_subjets_pt): if len(lds) != 1: # flatten and save offsets to unflatten afterwards - if type(lds.layout) == ak._ext.ListOffsetArray64: + if type(lds.layout) is ak._ext.ListOffsetArray64: ld_offsets = lds.kt.layout.offsets flat_subjet_pt = ak.flatten(kt_subjets_pt) - elif type(lds.layout) == ak._ext.ListArray64: + elif type(lds.layout) is ak._ext.ListArray64: ld_offsets = lds.layout.toListOffsetArray64(False).offsets flat_subjet_pt = kt_subjets_pt else: @@ -820,9 +823,6 @@ def _get_lund_arrays( Gets the ``num_prongs`` subjet pTs and Delta and kT per primary LP splitting of fatjets at ``fatjet_idx`` in each event. - Features are flattened (for now), and offsets are saved in ``ld_offsets`` to recover the event - structure. - Args: events (NanoEventsArray): nano events jec_fatjets (FatJetArray): post-JEC fatjets, used to update subjet pTs. @@ -830,7 +830,7 @@ def _get_lund_arrays( num_prongs (int): number of prongs / subjets per jet to reweight Returns: - flat_logD, flat_logkt, flat_subjet_pt, ld_offsets, kt_subjets_vec + lds, kt_subjets_vec, kt_subjets_pt: lund declusterings, subjet 4-vectors, JEC-corrected subjet pTs """ # jet definitions for LP SFs @@ -880,6 +880,27 @@ def _get_lund_arrays( return lds, kt_subjets_vec, kt_subjets_pt +def _get_flat_lund_arrays(events, jec_fatjet, fatjet_idx, num_prongs): + """Wrapper for the _get_lund_arrays and _get_flat_lp_vars functions + + returns: lds - lund declusterings, + kt_subjets_vec - subjet 4-vectors, + kt_subjets_pt - JEC-corrected subjet pTs, + ld_offsets - offsets for jagged structure, + flat_logD - flattened log(0.8/Delta), + flat_logkt - flattened log(kT/GeV), + flat_subjet_pt - flattened JEC-corrected subjet pTs + """ + lds, kt_subjets_vec, kt_subjets_pt = _get_lund_arrays( + events, jec_fatjet, fatjet_idx, num_prongs + ) + + lds_flat = ak.flatten(lds, axis=1) + ld_offsets, flat_logD, flat_logkt, flat_subjet_pt = _get_flat_lp_vars(lds_flat, kt_subjets_pt) + + return lds, kt_subjets_vec, kt_subjets_pt, ld_offsets, flat_logD, flat_logkt, flat_subjet_pt + + def _calc_lund_SFs( flat_logD: np.ndarray, flat_logkt: np.ndarray, @@ -969,6 +990,7 @@ def get_lund_SFs( fatjet_idx: int | ak.Array, num_prongs: int, gen_quarks: GenParticleArray, + weights: np.ndarray, seed: int = 42, trunc_gauss: bool = False, lnN: bool = True, @@ -984,6 +1006,8 @@ def get_lund_SFs( 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 + gen_quarks (GenParticleArray): gen quarks + weights (np.ndarray): event weights, for filling the LP histogram seed (int, optional): seed for random smearings. Defaults to 42. trunc_gauss (bool, optional): use truncated gaussians for smearing. Defaults to False. lnN (bool, optional): use log normals for smearings. Defaults to True. @@ -995,7 +1019,7 @@ def get_lund_SFs( """ # global variable to not have to load + smear LP ratios each time - global ratio_smeared_lookups, ratio_lnN_smeared_lookups, ratio_sys_up, ratio_sys_down, pt_extrap_lookups_dict, bratio, lp_year # noqa: PLW0603 + global ratio_smeared_lookups, ratio_lnN_smeared_lookups, ratio_sys_up, ratio_sys_down, pt_extrap_lookups_dict, bratio, ratio_edges, lp_year # noqa: PLW0603 if ( (lnN and ratio_lnN_smeared_lookups is None) @@ -1009,6 +1033,7 @@ def get_lund_SFs( ratio_sys_down, pt_extrap_lookups_dict, bratio, + ratio_edges, ) = _get_lund_lookups(year, seed, lnN, trunc_gauss) lp_year = year @@ -1016,16 +1041,38 @@ def get_lund_SFs( jec_fatjet = jec_fatjets[np.arange(len(jec_fatjets)), fatjet_idx] - lds, kt_subjets_vec, kt_subjets_pt = _get_lund_arrays( - events, jec_fatjet, fatjet_idx, num_prongs + # get lund plane declusterings, subjets, and flattened LP vars + lds, kt_subjets_vec, kt_subjets_pt, ld_offsets, flat_logD, flat_logkt, flat_subjet_pt = ( + _get_flat_lund_arrays(events, jec_fatjet, fatjet_idx, num_prongs) ) - lds_flat = ak.flatten(lds, axis=1) - ld_offsets, flat_logD, flat_logkt, flat_subjet_pt = _get_flat_lp_vars(lds_flat, kt_subjets_pt) + ################################################################################################ + # ---- Fill LP histogram for signal for distortion uncertainty ---- # + ################################################################################################ - sfs = {} + lp_hist = hist.Hist( + hist.axis.Variable(ratio_edges[0], name="subjet_pt", label="Subjet pT [GeV]"), + hist.axis.Variable(ratio_edges[1], name="logD", label="ln(0.8/Delta)"), + hist.axis.Variable(ratio_edges[2], name="logkt", label="ln(kT/GeV)"), + ) + + # repeat weights for each LP splitting + flat_weights = np.repeat( + np.repeat(weights, num_prongs), ak.count(ak.flatten(lds.kt, axis=1), axis=1) + ) + lp_hist.fill( + subjet_pt=flat_subjet_pt, + logD=flat_logD, + logkt=flat_logkt, + weight=flat_weights, + ) + + ################################################################################################ # ---- get scale factors per jet + smearings for stat unc. + syst. variations + pt extrap unc. ---- # + ################################################################################################ + + sfs = {} if trunc_gauss: sfs["lp_sf"] = _calc_lund_SFs( @@ -1079,7 +1126,29 @@ def get_lund_SFs( pt_extrap_lookups_dict["smeared_params"], ) + ################################################################################################ + # ---- get scale factors after re-clusteing with +/- one prong, for subjet matching uncs. ---- # + ################################################################################################ + + for shift, nps in [("down", num_prongs - 1), ("up", num_prongs + 1)]: + # get lund plane declusterings, subjets, and flattened LP vars + _, _, _, np_ld_offsets, np_flat_logD, np_flat_logkt, np_flat_subjet_pt = ( + _get_flat_lund_arrays(events, jec_fatjet, fatjet_idx, nps) + ) + + sfs[f"lp_sf_np_{shift}"] = _calc_lund_SFs( + np_flat_logD, + np_flat_logkt, + np_flat_subjet_pt, + np_ld_offsets, + nps, + ratio_lnN_smeared_lookups, + [pt_extrap_lookups_dict["params"]], + ) + + ################################################################################################ # ---- b-quark related uncertainties ---- # + ################################################################################################ if gen_bs is not None: assert ak.all( @@ -1125,7 +1194,9 @@ def get_lund_SFs( # weird edge case where b-subjet has no splittings sfs["lp_sfs_bl_ratio"] = 1.0 + ################################################################################################ # ---- subjet matching uncertainties ---- # + ################################################################################################ matching_dR = 0.2 sj_matched = [] @@ -1147,9 +1218,12 @@ def get_lund_SFs( 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.array(np.any(q_boundary, axis=1, keepdims=True)) + # events with quarks at the inside boundary of the jet + q_boundary = (j_q_dr > 0.7) * (j_q_dr <= 0.8) + sfs["lp_sf_inside_boundary_quarks"] = np.array(np.any(q_boundary, axis=1, keepdims=True)) + # events with quarks at the outside boundary of the jet + q_boundary = (j_q_dr > 0.8) * (j_q_dr <= 0.9) + sfs["lp_sf_outside_boundary_quarks"] = np.array(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( @@ -1162,4 +1236,4 @@ def get_lund_SFs( # 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 + return sfs, lp_hist diff --git a/src/HHbbVV/scale_factors/top_reweighting.ipynb b/src/HHbbVV/scale_factors/top_reweighting.ipynb index fcff82dd..94fd1a04 100644 --- a/src/HHbbVV/scale_factors/top_reweighting.ipynb +++ b/src/HHbbVV/scale_factors/top_reweighting.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -36,14 +36,12 @@ "import os\n", "\n", "from HHbbVV.processors import corrections\n", - "import correctionlib\n", - "\n", - "import utils" + "import correctionlib" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -53,48 +51,19 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "plot_dir = \"../../../plots/ScaleFactors/24Feb26\"\n", + "plot_dir = \"../../../plots/ScaleFactors/24Jul22\"\n", "_ = os.system(f\"mkdir -p {plot_dir}\")" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJetAK15SubJet_nBHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJetAK15SubJet_nCHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJetAK15_nBHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJetAK15_nCHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJet_btagDDBvLV2 in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJet_btagDDCvBV2 in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJet_btagDDCvLV2 in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJet_nBHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch FatJet_nCHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch SubJet_nBHadrons in , taking first instance\n", - " warnings.warn(\n", - "/uscms/home/rkansal/.local/lib/python3.9/site-packages/coffea/nanoevents/mapping/uproot.py:86: UserWarning: Found duplicate branch SubJet_nCHadrons in , taking first instance\n", - " warnings.warn(\n" - ] - } - ], + "outputs": [], "source": [ "events = nanoevents.NanoEventsFactory.from_root(\n", " # \"/eos/uscms/store/user/lpcpfnano/cmantill/v2_3/2017/HH_gen/GluGluToHHTobbVV_node_cHHH1_TuneCP5_13TeV-powheg-pythia8/GluGluToHHTobbVV_node_cHHH1/221017_221918/0000/nano_mc2017_100.root\",\n", @@ -118,7 +87,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -163,7 +132,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -181,7 +150,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -227,7 +196,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -259,7 +228,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -283,7 +252,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -337,7 +306,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -348,7 +317,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -365,7 +334,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -529,7 +498,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -554,7 +523,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -569,7 +538,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -592,7 +561,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -635,7 +604,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -680,7 +649,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -689,27 +658,7 @@ }, { "cell_type": "code", - "execution_count": 42, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 42, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "leading_fatjets" - ] - }, - { - "cell_type": "code", - "execution_count": 20, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -728,7 +677,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -743,7 +692,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -764,103 +713,57 @@ }, { "cell_type": "code", - "execution_count": 39, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "45000" - ] - }, - "execution_count": 39, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "len(events)" - ] - }, - { - "cell_type": "code", - "execution_count": 38, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "60" - ] - }, - "execution_count": 38, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "np.sum(top_match_dict[\"top_matched\"])" + "print(len(events))\n", + "print(np.sum(top_match_dict[\"top_matched\"]))\n", + "print(np.sum(top_match_dict[\"w_matched\"]))\n", + "print(np.sum(top_match_dict[\"unmatched\"]))" ] }, { "cell_type": "code", - "execution_count": 40, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "180" - ] - }, - "execution_count": 40, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "np.sum(top_match_dict[\"w_matched\"])" + "merged_top_events = presel_events[merged_top_jet_match]\n", + "had_top_jets = merged_top_events.FatJet[:, 0]\n", + "\n", + "merged_ak8_pfcands = merged_top_events.FatJetPFCands\n", + "merged_ak8_pfcands = merged_ak8_pfcands[\n", + " merged_ak8_pfcands.jetIdx == fatjet_idx[merged_top_jet_match]\n", + "]\n", + "merged_pfcands = merged_top_events.PFCands[merged_ak8_pfcands.pFCandsIdx]" ] }, { "cell_type": "code", - "execution_count": 41, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "855" - ] - }, - "execution_count": 41, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "np.sum(top_match_dict[\"unmatched\"])" + "merged_pfcands_vector_ptetaphi = ak.Array(\n", + " [\n", + " [{kin_key: cand[kin_key] for kin_key in skim_vars} for cand in event_cands]\n", + " for event_cands in merged_pfcands\n", + " ],\n", + " with_name=\"PtEtaPhiMLorentzVector\",\n", + ")" ] }, { - "cell_type": "code", - "execution_count": 23, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "merged_top_events = presel_events[merged_top_jet_match]\n", - "had_top_jets = merged_top_events.FatJet[:, 0]\n", - "\n", - "merged_ak8_pfcands = merged_top_events.FatJetPFCands\n", - "merged_ak8_pfcands = merged_ak8_pfcands[\n", - " merged_ak8_pfcands.jetIdx == fatjet_idx[merged_top_jet_match]\n", - "]\n", - "merged_pfcands = merged_top_events.PFCands[merged_ak8_pfcands.pFCandsIdx]" + "JEC Correction" ] }, { "cell_type": "code", - "execution_count": 24, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -887,22 +790,7 @@ }, { "cell_type": "code", - "execution_count": 25, - "metadata": {}, - "outputs": [], - "source": [ - "merged_pfcands_vector_ptetaphi = ak.Array(\n", - " [\n", - " [{kin_key: cand[kin_key] for kin_key in skim_vars} for cand in event_cands]\n", - " for event_cands in merged_pfcands\n", - " ],\n", - " with_name=\"PtEtaPhiMLorentzVector\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 26, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -916,36 +804,9 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "#--------------------------------------------------------------------------\n", - "# FastJet release 3.4.0\n", - "# M. Cacciari, G.P. Salam and G. Soyez \n", - "# A software package for jet finding and analysis at colliders \n", - "# http://fastjet.fr \n", - "#\t \n", - "# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package\n", - "# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. \n", - "# \n", - "# FastJet is provided without warranty under the GNU GPL v2 or higher. \n", - "# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code,\n", - "# CGAL and 3rd party plugin jet algorithms. See COPYING file for details.\n", - "#--------------------------------------------------------------------------\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "UserWarning: dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.\n" - ] - } - ], + "outputs": [], "source": [ "# cluster first with kT\n", "kt_clustering = fastjet.ClusterSequence(merged_pfcands_vector_ptetaphi, ktdef)\n", @@ -964,7 +825,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -983,7 +844,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -998,20 +859,9 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 30, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "closest_sjidx = np.argmin(subjet_bs_dr, axis=1).to_numpy()\n", "kt_subjets_pt[np.arange(len(kt_subjets_pt)), closest_sjidx]" @@ -1019,7 +869,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1035,40 +885,18 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 32, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "kt_subjets_pt" ] }, { "cell_type": "code", - "execution_count": 33, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 33, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "kt_subjets_pt[np.arange(len(kt_subjets_pt)), closest_sjidx]" ] @@ -1165,44 +993,6 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "f = uproot.open(\"../corrections/lp_ratios/ratio_2017.root\")\n", - "f.keys()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "f[\"ratio_sys_tot_down\"]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "f[\"pt_extrap\"][\"func_16_12\"]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fold = uproot.open(\"../corrections/lp_ratios/old/ratio_june9.root\")\n", - "fold.keys()" - ] - }, - { - "cell_type": "code", - "execution_count": 34, - "metadata": {}, - "outputs": [], "source": [ "f = uproot.open(\"../corrections/lp_ratios/ratio_2017.root\")\n", "# f = uproot.open(\"../corrections/ratio_june9.root\")\n", @@ -1218,211 +1008,9 @@ "ratio_nom_errs[zero_vals] = 0\n", "\n", "ratio_shape = f[\"ratio_nom\"].values().shape\n", - "bl_ratio = f[\"h_bl_ratio\"].to_numpy()[0]" - ] - }, - { - "cell_type": "code", - "execution_count": 36, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[[0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " ...,\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ]],\n", - "\n", - " [[0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " ...,\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ]],\n", - "\n", - " [[0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " ...,\n", - " [2. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ]],\n", - "\n", - " [[0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " ...,\n", - " [0. , 1.840388 , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ]],\n", - "\n", - " [[0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " ...,\n", - " [1.4767264 , 1.3192567 , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [1.0331354 , 1.4990366 , 1.9544313 , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ]],\n", - "\n", - " [[0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " ...,\n", - " [0. , 0.53056574, 2. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0.42190078, 0. , ..., 0. ,\n", - " 0. , 0. ],\n", - " [0. , 0. , 0. , ..., 0. ,\n", - " 0. , 0. ]]], dtype=float32)" - ] - }, - "execution_count": 36, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "bl_ratio" - ] - }, - { - "cell_type": "code", - "execution_count": 35, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[[1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " ...,\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ]],\n", - "\n", - " [[1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " ...,\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ]],\n", - "\n", - " [[1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " ...,\n", - " [0.81918097, 1.6875561 , 0.16781662, ..., 1. ,\n", - " 1. , 1. ],\n", - " [1.6282549 , 3.4878435 , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ]],\n", - "\n", - " [[1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " ...,\n", - " [0.20374994, 1. , 1.4421521 , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 0.6945209 , 2.0556872 , ..., 1. ,\n", - " 1. , 1. ],\n", - " [6. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ]],\n", - "\n", - " [[1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " ...,\n", - " [0.2418148 , 1.459049 , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1.5244532 , 0.805354 , 1.6780305 , ..., 1. ,\n", - " 1. , 1. ],\n", - " [4. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ]],\n", - "\n", - " [[1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " ...,\n", - " [1. , 1. , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [0.11765861, 2.4962053 , 1. , ..., 1. ,\n", - " 1. , 1. ],\n", - " [1. , 2. , 1. , ..., 1. ,\n", - " 1. , 1. ]]], dtype=float32)" - ] - }, - "execution_count": 35, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ratio_nom" + "bl_ratio = f[\"h_bl_ratio\"].to_numpy()[0]\n", + "\n", + "n_LP_sf_toys = 100" ] }, { @@ -1438,16 +1026,6 @@ "_ = plt.colorbar()" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# max_pt_bin = ratio_nom_edges[0][-2]\n", - "n_LP_sf_toys = 100" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1455,7 +1033,11 @@ "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(\n", + " ratio_nom[5].T[::-1],\n", + " extent=[-ratio_edges[1][0], ratio_edges[1][-1], ratio_edges[2][0], ratio_edges[2][-1]],\n", + " cmap=\"turbo\",\n", + ")\n", "plt.xlabel(r\"ln$(0.8/\\Delta)$\")\n", "plt.ylabel(r\"ln$(k_T/GeV)$\")\n", "_ = plt.colorbar()" @@ -1463,24 +1045,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 59, "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" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2\n", + "2\n", + "2\n", + "2\n" + ] + } + ], "source": [ "pt_extrap_lookups_dict = {\"params\": [], \"errs\": [], \"sys_up_params\": [], \"sys_down_params\": []}\n", "max_fparams = 4 # max order (+1) of functions\n", @@ -1545,9 +1123,52 @@ "\n", "# produces array of shape ``[n_sf_toys, subjet_pt bins, ln(0.8/Delta) bins, ln(kT/GeV) bins]``\n", "ratio_nom_smeared = ratio_nom + (ratio_nom_errs * rand_noise)\n", - "ratio_smeared_lookups = [\n", - " dense_lookup(ratio_nom_smeared[i], ratio_nom_edges) for i in range(n_sf_toys)\n", - "]" + "ratio_smeared_lookups = [dense_lookup(ratio_nom_smeared[i], ratio_edges) for i in range(n_sf_toys)]" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": {}, + "outputs": [], + "source": [ + "def _get_flat_lp_vars(lds, kt_subjets_pt):\n", + " if len(lds) != 1:\n", + " # flatten and save offsets to unflatten afterwards\n", + " if type(lds.layout) == ak._ext.ListOffsetArray64:\n", + " ld_offsets = lds.kt.layout.offsets\n", + " flat_subjet_pt = ak.flatten(kt_subjets_pt)\n", + " elif type(lds.layout) == ak._ext.ListArray64:\n", + " ld_offsets = lds.layout.toListOffsetArray64(False).offsets\n", + " flat_subjet_pt = kt_subjets_pt\n", + " else:\n", + " # edge case of single subjet...\n", + " ld_offsets = [0]\n", + " flat_subjet_pt = kt_subjets_pt\n", + "\n", + " # repeat subjet pt for each lund declustering\n", + " flat_subjet_pt = np.repeat(flat_subjet_pt, ak.count(lds.kt, axis=1)).to_numpy()\n", + " flat_logD = np.log(0.8 / ak.flatten(lds).Delta).to_numpy()\n", + " flat_logkt = np.log(ak.flatten(lds).kt).to_numpy()\n", + "\n", + " return ld_offsets, flat_logD, flat_logkt, flat_subjet_pt" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [], + "source": [ + "lds_flat = ak.flatten(lds, axis=1)\n", + "ld_offsets, flat_logD, flat_logkt, flat_subjet_pt = _get_flat_lp_vars(lds_flat, kt_subjets_pt)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Fill LP histogram for signal" ] }, { @@ -1556,20 +1177,563 @@ "metadata": {}, "outputs": [], "source": [ - "# save offsets to recover awkward structure later\n", - "ld_offsets = lds.kt.layout.offsets\n", - "flat_logD = np.log(0.8 / ak.flatten(lds).Delta).to_numpy()\n", - "flat_logkt = np.log(ak.flatten(lds).kt).to_numpy()\n", - "# repeat subjet pt for each lund declustering\n", - "flat_subjet_pt = np.repeat(ak.flatten(kt_subjets_pt), ak.count(lds.kt, axis=1)).to_numpy()" + "import hist\n", + "\n", + "lp_hist = hist.Hist(\n", + " hist.axis.Variable(ratio_edges[0], name=\"subjet_pt\", label=\"Subjet pT [GeV]\"),\n", + " hist.axis.Variable(ratio_edges[1], name=\"logD\", label=\"ln(0.8/Delta)\"),\n", + " hist.axis.Variable(ratio_edges[2], name=\"logkt\", label=\"ln(kT/GeV)\"),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "metadata": {}, + "outputs": [], + "source": [ + "weight = merged_top_events.genWeight.to_numpy()\n", + "flat_weight = np.repeat(np.repeat(weight, num_prongs), ak.count(lds_flat.kt, axis=1))" + ] + }, + { + "cell_type": "code", + "execution_count": 93, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Hist(\n", + " Variable([10, 65, 110, 175, 240, 300, 99999], name='subjet_pt', label='Subjet pT [GeV]'),\n", + " Variable([-1, -0.55, -0.1, 0.35, 0.8, 1.25, 1.7, 2.15, 2.6, 3.05, 3.5, 3.95, 4.4, 4.85, 5.3, 5.75, 6.2, 6.65, 7.1, 7.55, 8], name='logD', label='ln(0.8/Delta)'),\n", + " Variable([-5, -4.40461, -3.80922, -3.21384, -2.61845, -2.02306, -1.42767, -0.832286, -0.236898, 0.35849, 0.953878, 1.54927, 2.14465, 2.74004, 3.33543, 3.93082, 4.5262, 5.12159, 5.71698, 6.31237, 6.90776], name='logkt', label='ln(kT/GeV)'),\n", + " storage=Double()) # Sum: 435622.09045410156 (435925.44845581055 with flow)" + ] + }, + "execution_count": 93, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lp_hist.fill(\n", + " subjet_pt=flat_subjet_pt,\n", + " logD=flat_logD,\n", + " logkt=flat_logkt,\n", + " weight=flat_weight,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 97, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + "
\n", + "\n", + "\n", + "-1\n", + "\n", + "\n", + "8\n", + "\n", + "\n", + "-5\n", + "\n", + "\n", + "6.91\n", + "\n", + "\n", + "ln(0.8/Delta)\n", + "\n", + "\n", + "ln(kT/GeV)\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
\n", + "
\n", + "Variable([-1, -0.55, -0.1, 0.35, 0.8, 1.25, 1.7, 2.15, 2.6, 3.05, 3.5, 3.95, 4.4, 4.85, 5.3, 5.75, 6.2, 6.65, 7.1, 7.55, 8], name='logD', label='ln(0.8/Delta)')
\n", + "Variable([-5, -4.40461, -3.80922, -3.21384, -2.61845, -2.02306, -1.42767, -0.832286, -0.236898, 0.35849, 0.953878, 1.54927, 2.14465, 2.74004, 3.33543, 3.93082, 4.5262, 5.12159, 5.71698, 6.31237, 6.90776], name='logkt', label='ln(kT/GeV)')
\n", + "
\n", + "Double() Σ=39436.54022216797\n", + "\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + "Hist(\n", + " Variable([-1, -0.55, -0.1, 0.35, 0.8, 1.25, 1.7, 2.15, 2.6, 3.05, 3.5, 3.95, 4.4, 4.85, 5.3, 5.75, 6.2, 6.65, 7.1, 7.55, 8], name='logD', label='ln(0.8/Delta)'),\n", + " Variable([-5, -4.40461, -3.80922, -3.21384, -2.61845, -2.02306, -1.42767, -0.832286, -0.236898, 0.35849, 0.953878, 1.54927, 2.14465, 2.74004, 3.33543, 3.93082, 4.5262, 5.12159, 5.71698, 6.31237, 6.90776], name='logkt', label='ln(kT/GeV)'),\n", + " storage=Double()) # Sum: 39436.54022216797" + ] + }, + "execution_count": 97, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lp_hist[0, ...]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Testing summing over histograms\n", + "from copy import deepcopy\n", + "\n", + "lp_hist2 = deepcopy(lp_hist)\n", + "lp_hist2.values()[:] = 1\n", + "lp_hist2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Hist(\n", + " Variable([10, 65, 110, 175, 240, 300, 99999], name='subjet_pt', label='Subjet pT [GeV]'),\n", + " Variable([-1, -0.55, -0.1, 0.35, 0.8, 1.25, 1.7, 2.15, 2.6, 3.05, 3.5, 3.95, 4.4, 4.85, 5.3, 5.75, 6.2, 6.65, 7.1, 7.55, 8], name='logD', label='ln(0.8/Delta)'),\n", + " Variable([-5, -4.40461, -3.80922, -3.21384, -2.61845, -2.02306, -1.42767, -0.832286, -0.236898, 0.35849, 0.953878, 1.54927, 2.14465, 2.74004, 3.33543, 3.93082, 4.5262, 5.12159, 5.71698, 6.31237, 6.90776], name='logkt', label='ln(kT/GeV)'),\n", + " storage=Double()) # Sum: 438022.09045410156 (438628.80645751953 with flow)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sum([lp_hist, lp_hist2])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "pT Extrapolation" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": {}, "outputs": [], "source": [ + "max_pt_bin = ratio_edges[0][-2]\n", "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", @@ -1581,11 +1745,11 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 62, "metadata": {}, "outputs": [], "source": [ - "clip_max, clip_min = 10, 0.1\n", + "clip_max, clip_min = 5.0, 1 / 5.0\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)" @@ -1803,7 +1967,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.11" + "version": "3.9.15" }, "orig_nbformat": 4, "vscode": {