working out pt extrapolation

rkansal47 · Jun 10, 2023 · d453c33 · d453c33
1 parent 53af5ea
commit d453c33
Showing 1 changed file with 62 additions and 9 deletions.
diff --git a/src/HHbbVV/scale_factors/top_reweighting.ipynb b/src/HHbbVV/scale_factors/top_reweighting.ipynb
@@ -1033,26 +1033,78 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": null,
+   "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"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "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"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 78,
+   "metadata": {},
+   "outputs": [],
+   "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)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 74,
    "metadata": {},
    "outputs": [
     {
      "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApCElEQVR4nO3df1BU973/8RciLCgsiAYQRatiNQaJia2ENFWrVCVer7kyc9V0Gsx19CaX2CppftCbH8U0g03zw7ZDTCZjNHcm1sSMP5reqK1OwZsWjCFSNT8csaYaI3hjhBXEBfHz/aNf92YDKgtnP7D4fMycGc/nfPaz788edvfl7jlnw4wxRgAAAJb06e4CAADA9YXwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMCqvt1dwNddunRJn3/+uWJjYxUWFtbd5QAAgA4wxujcuXNKSUlRnz5X/2yjx4WPzz//XKmpqd1dBgAA6IQTJ05o6NChV+3T48JHbGyspH8U73a7u7kaAADQER6PR6mpqb738avpceHj8lctbreb8AEAQIjpyCETHHAKAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACs6nG/agsAnXWyrklnG5sdH3dA/0gNiY92fFzgekX4ANArnKxrUvZzZWpqaXV87OiIcO16cAoBBHAI4QNAr3C2sVlNLa1aPX+C0hJjHBu3+nSDlr9RpbONzYQPwCGEDwC9SlpijNKHxHV3GQCuggNOAQCAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYFWXwseqVasUFham5cuX+9ouXLig/Px8DRw4UDExMcrNzVVtbW1X6wQAAL1Ep8PHvn379PLLLysjI8OvfcWKFXr77be1adMmlZWV6fPPP9e8efO6XCgAAOgdOhU+Ghoa9IMf/ECvvPKKBgwY4Guvr6/X2rVr9fzzz2vatGmaOHGi1q1bp7/85S+qqKhwrGgAABC6OhU+8vPzNXv2bGVnZ/u1V1ZWqqWlxa997NixGjZsmMrLy9sdy+v1yuPx+C0AAKD36hvoDTZu3KgPPvhA+/bta7OtpqZGkZGRio+P92tPSkpSTU1Nu+MVFxerqKgo0DIAAECICuiTjxMnTujHP/6xXn/9dUVFRTlSQGFhoerr633LiRMnHBkXAAD0TAGFj8rKSp0+fVq33nqr+vbtq759+6qsrEy//vWv1bdvXyUlJam5uVl1dXV+t6utrVVycnK7Y7pcLrndbr8FAAD0XgF97TJ9+nQdPHjQr+3ee+/V2LFj9cgjjyg1NVURERHavXu3cnNzJUmHDx/W8ePHlZWV5VzVAAAgZAUUPmJjY5Wenu7X1r9/fw0cONDXvnjxYhUUFCghIUFut1vLli1TVlaWbrvtNueqBgAAISvgA06v5YUXXlCfPn2Um5srr9ermTNn6sUXX3T6bgAAQIjqcvgoLS31W4+KilJJSYlKSkq6OjQAAOiF+G0XAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYJXjl1cHgKs5Wdeks43Njo9bfbrB8TEBBAfhA4A1J+ualP1cmZpaWoMyfnREuAb0jwzK2ACcQ/gAYM3ZxmY1tbRq9fwJSkuMcXz8Af0jNSQ+2vFxATiL8AHAurTEGKUPievuMgB0Ew44BQAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVAYWPNWvWKCMjQ263W263W1lZWdq+fbtv+9SpUxUWFua33HfffY4XDQAAQlffQDoPHTpUq1at0ujRo2WM0Wuvvaa5c+dq//79uummmyRJS5Ys0cqVK3236devn7MVAwCAkBZQ+JgzZ47f+tNPP601a9aooqLCFz769eun5ORk5yoEAAC9SqeP+WhtbdXGjRvV2NiorKwsX/vrr7+uQYMGKT09XYWFhTp//vxVx/F6vfJ4PH4LAADovQL65EOSDh48qKysLF24cEExMTHasmWLxo0bJ0m6++67NXz4cKWkpOjAgQN65JFHdPjwYW3evPmK4xUXF6uoqKjzMwAAACEl4PAxZswYVVVVqb6+Xm+99Zby8vJUVlamcePGaenSpb5+48eP1+DBgzV9+nQdPXpUo0aNane8wsJCFRQU+NY9Ho9SU1M7MRUAABAKAg4fkZGRSktLkyRNnDhR+/bt069+9Su9/PLLbfpmZmZKkqqrq68YPlwul1wuV6BlAACAENXl63xcunRJXq+33W1VVVWSpMGDB3f1bgAAQC8R0CcfhYWFysnJ0bBhw3Tu3Dlt2LBBpaWl2rlzp44ePaoNGzbozjvv1MCBA3XgwAGtWLFCkydPVkZGRrDqBwAAISag8HH69Gndc889OnXqlOLi4pSRkaGdO3fq+9//vk6cOKFdu3Zp9erVamxsVGpqqnJzc/XYY48Fq3YAABCCAgofa9euveK21NRUlZWVdbkgAADQu/HbLgAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsCih8rFmzRhkZGXK73XK73crKytL27dt92y9cuKD8/HwNHDhQMTExys3NVW1treNFAwCA0BVQ+Bg6dKhWrVqlyspKvf/++5o2bZrmzp2rDz/8UJK0YsUKvf3229q0aZPKysr0+eefa968eUEpHAAAhKa+gXSeM2eO3/rTTz+tNWvWqKKiQkOHDtXatWu1YcMGTZs2TZK0bt063XjjjaqoqNBtt93mXNUAACBkdfqYj9bWVm3cuFGNjY3KyspSZWWlWlpalJ2d7eszduxYDRs2TOXl5Vccx+v1yuPx+C0AAKD3Cjh8HDx4UDExMXK5XLrvvvu0ZcsWjRs3TjU1NYqMjFR8fLxf/6SkJNXU1FxxvOLiYsXFxfmW1NTUgCcBAABCR8DhY8yYMaqqqtLevXt1//33Ky8vTx999FGnCygsLFR9fb1vOXHiRKfHAgAAPV9Ax3xIUmRkpNLS0iRJEydO1L59+/SrX/1K8+fPV3Nzs+rq6vw+/aitrVVycvIVx3O5XHK5XIFXDgAAQlKXr/Nx6dIleb1eTZw4UREREdq9e7dv2+HDh3X8+HFlZWV19W4AAEAvEdAnH4WFhcrJydGwYcN07tw5bdiwQaWlpdq5c6fi4uK0ePFiFRQUKCEhQW63W8uWLVNWVhZnugAAAJ+Awsfp06d1zz336NSpU4qLi1NGRoZ27typ73//+5KkF154QX369FFubq68Xq9mzpypF198MSiFAwCA0BRQ+Fi7du1Vt0dFRamkpEQlJSVdKgoAAPRe/LYLAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsCCh/FxcX69re/rdjYWCUmJuquu+7S4cOH/fpMnTpVYWFhfst9993naNEAACB0BRQ+ysrKlJ+fr4qKCv3xj39US0uLZsyYocbGRr9+S5Ys0alTp3zLM88842jRAAAgdPUNpPOOHTv81tevX6/ExERVVlZq8uTJvvZ+/fopOTnZmQoBAECv0qVjPurr6yVJCQkJfu2vv/66Bg0apPT0dBUWFur8+fNXHMPr9crj8fgtAACg9wrok4+vunTpkpYvX67vfOc7Sk9P97XffffdGj58uFJSUnTgwAE98sgjOnz4sDZv3tzuOMXFxSoqKupsGQAAIMR0Onzk5+fr0KFDevfdd/3aly5d6vv3+PHjNXjwYE2fPl1Hjx7VqFGj2oxTWFiogoIC37rH41FqampnywIAAD1cp8LHAw88oN///vfas2ePhg4detW+mZmZkqTq6up2w4fL5ZLL5epMGQAAIAQFFD6MMVq2bJm2bNmi0tJSjRgx4pq3qaqqkiQNHjy4UwUCAIDeJaDwkZ+frw0bNmjbtm2KjY1VTU2NJCkuLk7R0dE6evSoNmzYoDvvvFMDBw7UgQMHtGLFCk2ePFkZGRlBmQAAAAgtAYWPNWvWSPrHhcS+at26dVq0aJEiIyO1a9curV69Wo2NjUpNTVVubq4ee+wxxwoGAAChLeCvXa4mNTVVZWVlXSoIAAD0bvy2CwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCqU79qC6BnOFnXpLONzUEZe0D/SA2Jjw7K2ACub4QPIESdrGtS9nNlamppDcr40RHh2vXgFAIIAMcRPoAQdbaxWU0trVo9f4LSEmMcHbv6dIOWv1Gls43NhA8AjiN8ACEuLTFG6UPiursMAOgwDjgFAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYFFD6Ki4v17W9/W7GxsUpMTNRdd92lw4cP+/W5cOGC8vPzNXDgQMXExCg3N1e1tbWOFg0AAEJXQOGjrKxM+fn5qqio0B//+Ee1tLRoxowZamxs9PVZsWKF3n77bW3atEllZWX6/PPPNW/ePMcLBwAAoalvIJ137Njht75+/XolJiaqsrJSkydPVn19vdauXasNGzZo2rRpkqR169bpxhtvVEVFhW677TbnKgcAACGpS8d81NfXS5ISEhIkSZWVlWppaVF2dravz9ixYzVs2DCVl5d35a4AAEAvEdAnH1916dIlLV++XN/5zneUnp4uSaqpqVFkZKTi4+P9+iYlJammpqbdcbxer7xer2/d4/F0tiQACJrq0w1BGXdA/0gNiY8OythAT9Xp8JGfn69Dhw7p3Xff7VIBxcXFKioq6tIYABAsA/pHKjoiXMvfqArK+NER4dr14BQCCK4rnQofDzzwgH7/+99rz549Gjp0qK89OTlZzc3Nqqur8/v0o7a2VsnJye2OVVhYqIKCAt+6x+NRampqZ8oCAMcNiY/Wrgen6Gxjs+NjV59u0PI3qnS2sZnwgetKQOHDGKNly5Zpy5YtKi0t1YgRI/y2T5w4UREREdq9e7dyc3MlSYcPH9bx48eVlZXV7pgul0sul6uT5QNA8A2JjyYcAA4KKHzk5+drw4YN2rZtm2JjY33HccTFxSk6OlpxcXFavHixCgoKlJCQILfbrWXLlikrK4szXQAAgKQAw8eaNWskSVOnTvVrX7dunRYtWiRJeuGFF9SnTx/l5ubK6/Vq5syZevHFFx0pFgAAhL6Av3a5lqioKJWUlKikpKTTRQEAgN6L33YBAABWdfpUWwAdd7KuyfGzJYJ13QkACDbCBxBkJ+ualP1cmZpaWh0fOzoiXAP6Rzo+LgAEE+EDCLKzjc1qamnV6vkTlJYY4+jYXB0TQCgifACWpCXGKH1IXHeXAQDdjgNOAQCAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFV9u7sAAD1X9emGHj0egNBE+ADQxoD+kYqOCNfyN6ocHzs6IlwD+kc6Pi6A0EH4ANDGkPho7Xpwis42Njs+9oD+kRoSH+34uABCB+EDQLuGxEcTEgAEBQecAgAAqwIOH3v27NGcOXOUkpKisLAwbd261W/7okWLFBYW5rfMmjXLqXoBAECICzh8NDY26uabb1ZJSckV+8yaNUunTp3yLb/97W+7VCQAAOg9Aj7mIycnRzk5OVft43K5lJyc3OmiAABA7xWUYz5KS0uVmJioMWPG6P7779eZM2eCcTcAACAEOX62y6xZszRv3jyNGDFCR48e1U9/+lPl5OSovLxc4eHhbfp7vV55vV7fusfjcbokAADQgzgePhYsWOD79/jx45WRkaFRo0aptLRU06dPb9O/uLhYRUVFTpcBAAB6qKCfajty5EgNGjRI1dXV7W4vLCxUfX29bzlx4kSwSwIAAN0o6BcZ++yzz3TmzBkNHjy43e0ul0sulyvYZQAAgB4i4PDR0NDg9ynGsWPHVFVVpYSEBCUkJKioqEi5ublKTk7W0aNH9fDDDystLU0zZ850tHAAABCaAg4f77//vr73ve/51gsKCiRJeXl5WrNmjQ4cOKDXXntNdXV1SklJ0YwZM/TUU0/x6QYAAJDUifAxdepUGWOuuH3nzp1dKggAAPRu/LYLAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKzq290F2HayrklnG5sdH3dA/0gNiY92fFwAAHqb6yp8nKxrUvZzZWpqaXV87OiIcO16cAoBBACAa7iuwsfZxmY1tbRq9fwJSkuMcWzc6tMNWv5Glc42NhM+AAC4husqfFyWlhij9CFx3V0GAADXJQ44BQAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYFHD727NmjOXPmKCUlRWFhYdq6davfdmOMnnjiCQ0ePFjR0dHKzs7WkSNHnKoXAACEuIDDR2Njo26++WaVlJS0u/2ZZ57Rr3/9a7300kvau3ev+vfvr5kzZ+rChQtdLhYAAIS+gH/bJScnRzk5Oe1uM8Zo9erVeuyxxzR37lxJ0n/9138pKSlJW7du1YIFC7pWLQAACHmOHvNx7Ngx1dTUKDs729cWFxenzMxMlZeXt3sbr9crj8fjtwAAgN7L0fBRU1MjSUpKSvJrT0pK8m37uuLiYsXFxfmW1NRUJ0sCAAA9TLef7VJYWKj6+nrfcuLEie4uCQAABJGj4SM5OVmSVFtb69deW1vr2/Z1LpdLbrfbbwEAAL2Xo+FjxIgRSk5O1u7du31tHo9He/fuVVZWlpN3BQAAQlTAZ7s0NDSourrat37s2DFVVVUpISFBw4YN0/Lly/Xzn/9co0eP1ogRI/T4448rJSVFd911l5N1AwCAEBVw+Hj//ff1ve99z7deUFAgScrLy9P69ev18MMPq7GxUUuXLlVdXZ3uuOMO7dixQ1FRUc5VDQAAQlbA4WPq1Kkyxlxxe1hYmFauXKmVK1d2qTAAANA7dfvZLgAA4PpC+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFb17e4CgJ7iZF2TzjY2Oz5u9ekGx8cEgFBG+AD0j+CR/VyZmlpagzJ+dES4BvSPDMrYABBqCB+ApLONzWpqadXq+ROUlhjj+PgD+kdqSHy04+MCQCgifABfkZYYo/Qhcd1dBgD0ahxwCgAArHI8fPzsZz9TWFiY3zJ27Fin7wYAAISooHztctNNN2nXrl3/dyd9+XYHAAD8Q1BSQd++fZWcnByMoQEAQIgLSvg4cuSIUlJSFBUVpaysLBUXF2vYsGHt9vV6vfJ6vb51j8cTjJJgWbCumSFx5ggAhDrHw0dmZqbWr1+vMWPG6NSpUyoqKtJ3v/tdHTp0SLGxsW36FxcXq6ioyOky0I1sXDNj14NTCCAAEKIcDx85OTm+f2dkZCgzM1PDhw/Xm2++qcWLF7fpX1hYqIKCAt+6x+NRamqq02XBomBeM6P6dIOWv1Gls43NhA8ACFFBPxI0Pj5e3/zmN1VdXd3udpfLJZfLFewy0A24ZgYAoD1Bv85HQ0ODjh49qsGDBwf7rgAAQAhwPHz85Cc/UVlZmT799FP95S9/0b/8y78oPDxcCxcudPquAABACHL8a5fPPvtMCxcu1JkzZ3TDDTfojjvuUEVFhW644Qan7woAAIQgx8PHxo0bnR4SAAD0Ivy2CwAAsIrrngNAN6s+3RCUcbkgH3oqwgcAdJMB/SMVHRGu5W9UBWV8LsiHnorwAQDdZEh8tHY9OCUoP0XABfnQkxE+AKAbDYmPJhzgusMBpwAAwCrCBwAAsIrwAQAArCJ8AAAAqzjg1EGcq2+P0491sPYd0N1C7W+b17vrA+HDAZyrb08wH+voiHAN6B/p+LhAdwj261Kw8Hp3fSB8OIBz9e0J5mPN/7jQmwTzuRIsvN5dPwgfDuFcfXt4rIGO4bmCnooDTgEAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVZztEiJC6UJBoVQrgOvDybqmoJ12zGn6gSN89HChfKEgLtgFoCc4Wdek7OfK1NTSGpTxuTBa4AgfPVwoXihI4n8CAHqOs43Nampp1er5E5SWGOPo2FwYrXMIHyGACwUBQNelJcYofUhcd5cBccApAACwjPABAACsInwAAACrCB8AAMAqDjgFAPQoTl8riGsP9TyEDwBAjxDM6xpx7aGehfABAOgRgnldI6491LMQPgAAPQbXNbo+cMApAACwKmjho6SkRN/4xjcUFRWlzMxMvffee8G6KwAAEEKCEj7eeOMNFRQU6Mknn9QHH3ygm2++WTNnztTp06eDcXcAACCEBCV8PP/881qyZInuvfdejRs3Ti+99JL69eunV199NRh3BwAAQojjB5w2NzersrJShYWFvrY+ffooOztb5eXlbfp7vV55vV7fen19vSTJ4/E4XZoaznl0yXteDec88njCHB8fAHB9ufy+cuBvp9Rwzvn3rWC5IcalG9xRjo55+X3bGHPNvo6Hjy+++EKtra1KSkrya09KStInn3zSpn9xcbGKioratKempjpdmk/W6qANDQC4Dv1gdXdX0HOcO3dOcXFX//Xgbj/VtrCwUAUFBb71S5cu6csvv9TAgQMVFubspxMej0epqak6ceKE3G63o2P3BL19flLvnyPzC329fY7ML/QFa47GGJ07d04pKSnX7Ot4+Bg0aJDCw8NVW1vr115bW6vk5OQ2/V0ul1wul19bfHy802X5cbvdvfaPSur985N6/xyZX+jr7XNkfqEvGHO81icelzl+wGlkZKQmTpyo3bt3+9ouXbqk3bt3Kysry+m7AwAAISYoX7sUFBQoLy9P3/rWtzRp0iStXr1ajY2Nuvfee4NxdwAAIIQEJXzMnz9f//u//6snnnhCNTU1mjBhgnbs2NHmIFTbXC6XnnzyyTZf8/QWvX1+Uu+fI/MLfb19jswv9PWEOYaZjpwTAwAA4BB+2wUAAFhF+AAAAFYRPgAAgFWEDwAAYFWvCh9PP/20br/9dvXr16/DFyozxuiJJ57Q4MGDFR0drezsbB05csSvz5dffqkf/OAHcrvdio+P1+LFi9XQ0BCEGVxboLV8+umnCgsLa3fZtGmTr1972zdu3GhjSn4681hPnTq1Te333XefX5/jx49r9uzZ6tevnxITE/XQQw/p4sWLwZxKuwKd35dffqlly5ZpzJgxio6O1rBhw/SjH/3I9xtIl3Xn/ispKdE3vvENRUVFKTMzU++9995V+2/atEljx45VVFSUxo8fr3feecdve0eekzYFMr9XXnlF3/3udzVgwAANGDBA2dnZbfovWrSozb6aNWtWsKdxVYHMcf369W3qj4ry/42QUN6H7b2ehIWFafbs2b4+PWkf7tmzR3PmzFFKSorCwsK0devWa96mtLRUt956q1wul9LS0rR+/fo2fQJ9XgfM9CJPPPGEef75501BQYGJi4vr0G1WrVpl4uLizNatW81f//pX88///M9mxIgRpqmpyddn1qxZ5uabbzYVFRXmf/7nf0xaWppZuHBhkGZxdYHWcvHiRXPq1Cm/paioyMTExJhz5875+kky69at8+v31cfAls481lOmTDFLlizxq72+vt63/eLFiyY9Pd1kZ2eb/fv3m3feeccMGjTIFBYWBns6bQQ6v4MHD5p58+aZ3/3ud6a6utrs3r3bjB492uTm5vr16679t3HjRhMZGWleffVV8+GHH5olS5aY+Ph4U1tb227/P//5zyY8PNw888wz5qOPPjKPPfaYiYiIMAcPHvT16chz0pZA53f33XebkpISs3//fvPxxx+bRYsWmbi4OPPZZ5/5+uTl5ZlZs2b57asvv/zS1pTaCHSO69atM26326/+mpoavz6hvA/PnDnjN7dDhw6Z8PBws27dOl+fnrQP33nnHfOf//mfZvPmzUaS2bJly1X7/+1vfzP9+vUzBQUF5qOPPjK/+c1vTHh4uNmxY4evT6CPWWf0qvBx2bp16zoUPi5dumSSk5PNL3/5S19bXV2dcblc5re//a0xxpiPPvrISDL79u3z9dm+fbsJCwszJ0+edLz2q3GqlgkTJph/+7d/82vryB9tsHV2flOmTDE//vGPr7j9nXfeMX369PF7gVyzZo1xu93G6/U6UntHOLX/3nzzTRMZGWlaWlp8bd21/yZNmmTy8/N9662trSYlJcUUFxe32/9f//VfzezZs/3aMjMzzb//+78bYzr2nLQp0Pl93cWLF01sbKx57bXXfG15eXlm7ty5TpfaaYHO8Vqvr71tH77wwgsmNjbWNDQ0+Np62j68rCOvAw8//LC56aab/Nrmz59vZs6c6Vvv6mPWEb3qa5dAHTt2TDU1NcrOzva1xcXFKTMzU+Xl5ZKk8vJyxcfH61vf+pavT3Z2tvr06aO9e/dardeJWiorK1VVVaXFixe32Zafn69BgwZp0qRJevXVVzv0s8hO6sr8Xn/9dQ0aNEjp6ekqLCzU+fPn/cYdP36830XuZs6cKY/How8//ND5iVyBU39L9fX1crvd6tvX/xqBtvdfc3OzKisr/Z4/ffr0UXZ2tu/583Xl5eV+/aV/7IvL/TvynLSlM/P7uvPnz6ulpUUJCQl+7aWlpUpMTNSYMWN0//3368yZM47W3lGdnWNDQ4OGDx+u1NRUzZ071+951Nv24dq1a7VgwQL179/fr72n7MNAXes56MRj1hHd/qu23ammpkaS2lx5NSkpybetpqZGiYmJftv79u2rhIQEXx9bnKhl7dq1uvHGG3X77bf7ta9cuVLTpk1Tv3799Ic//EH/8R//oYaGBv3oRz9yrP5r6ez87r77bg0fPlwpKSk6cOCAHnnkER0+fFibN2/2jdvePr68zRYn9t8XX3yhp556SkuXLvVr747998UXX6i1tbXdx/aTTz5p9zZX2hdffb5dbrtSH1s6M7+ve+SRR5SSkuL3Qj5r1izNmzdPI0aM0NGjR/XTn/5UOTk5Ki8vV3h4uKNzuJbOzHHMmDF69dVXlZGRofr6ej377LO6/fbb9eGHH2ro0KG9ah++9957OnTokNauXevX3pP2YaCu9Bz0eDxqamrS2bNnu/x33xE9Pnw8+uij+sUvfnHVPh9//LHGjh1rqSLndXSOXdXU1KQNGzbo8ccfb7Ptq2233HKLGhsb9ctf/tKRN69gz++rb8Tjx4/X4MGDNX36dB09elSjRo3q9LgdZWv/eTwezZ49W+PGjdPPfvYzv23B3H/onFWrVmnjxo0qLS31OyBzwYIFvn+PHz9eGRkZGjVqlEpLSzV9+vTuKDUgWVlZfj8Sevvtt+vGG2/Uyy+/rKeeeqobK3Pe2rVrNX78eE2aNMmvPdT3YU/Q48PHgw8+qEWLFl21z8iRIzs1dnJysiSptrZWgwcP9rXX1tZqwoQJvj6nT5/2u93Fixf15Zdf+m7fVR2dY1dreeutt3T+/Hndc8891+ybmZmpp556Sl6vt8vX/7c1v8syMzMlSdXV1Ro1apSSk5PbHKldW1srSY7sQxvzO3funGbNmqXY2Fht2bJFERERV+3v5P67kkGDBik8PNz3WF5WW1t7xfkkJydftX9HnpO2dGZ+lz377LNatWqVdu3apYyMjKv2HTlypAYNGqTq6mrrb1xdmeNlERERuuWWW1RdXS2p9+zDxsZGbdy4UStXrrzm/XTnPgzUlZ6Dbrdb0dHRCg8P7/LfRIc4dvRIDxLoAafPPvusr62+vr7dA07ff/99X5+dO3d26wGnna1lypQpbc6SuJKf//znZsCAAZ2utTOceqzfffddI8n89a9/Ncb83wGnXz1S++WXXzZut9tcuHDBuQlcQ2fnV19fb2677TYzZcoU09jY2KH7srX/Jk2aZB544AHfemtrqxkyZMhVDzj9p3/6J7+2rKysNgecXu05aVOg8zPGmF/84hfG7Xab8vLyDt3HiRMnTFhYmNm2bVuX6+2Mzszxqy5evGjGjBljVqxYYYzpHfvQmH+8j7hcLvPFF19c8z66ex9epg4ecJqenu7XtnDhwjYHnHblb6JDtTo2Ug/w97//3ezfv993Kun+/fvN/v37/U4pHTNmjNm8ebNvfdWqVSY+Pt5s27bNHDhwwMydO7fdU21vueUWs3fvXvPuu++a0aNHd+uptler5bPPPjNjxowxe/fu9bvdkSNHTFhYmNm+fXubMX/3u9+ZV155xRw8eNAcOXLEvPjii6Zfv37miSeeCPp8vi7Q+VVXV5uVK1ea999/3xw7dsxs27bNjBw50kyePNl3m8un2s6YMcNUVVWZHTt2mBtuuKHbTrUNZH719fUmMzPTjB8/3lRXV/ud2nfx4kVjTPfuv40bNxqXy2XWr19vPvroI7N06VITHx/vO7Pohz/8oXn00Ud9/f/85z+bvn37mmeffdZ8/PHH5sknn2z3VNtrPSdtCXR+q1atMpGRkeatt97y21eXX4POnTtnfvKTn5jy8nJz7Ngxs2vXLnPrrbea0aNHWw3CXZljUVGR2blzpzl69KiprKw0CxYsMFFRUebDDz/09QnlfXjZHXfcYebPn9+mvaftw3Pnzvne6ySZ559/3uzfv9/8/e9/N8YY8+ijj5of/vCHvv6XT7V96KGHzMcff2xKSkraPdX2ao+ZE3pV+MjLyzOS2ix/+tOffH30/6+HcNmlS5fM448/bpKSkozL5TLTp083hw8f9hv3zJkzZuHChSYmJsa43W5z7733+gUam65Vy7Fjx9rM2RhjCgsLTWpqqmltbW0z5vbt282ECRNMTEyM6d+/v7n55pvNSy+91G7fYAt0fsePHzeTJ082CQkJxuVymbS0NPPQQw/5XefDGGM+/fRTk5OTY6Kjo82gQYPMgw8+6Heqqi2Bzu9Pf/pTu3/TksyxY8eMMd2//37zm9+YYcOGmcjISDNp0iRTUVHh2zZlyhSTl5fn1//NN9803/zmN01kZKS56aabzH//93/7be/Ic9KmQOY3fPjwdvfVk08+aYwx5vz582bGjBnmhhtuMBEREWb48OFmyZIljr6od0Ygc1y+fLmvb1JSkrnzzjvNBx984DdeKO9DY4z55JNPjCTzhz/8oc1YPW0fXuk14vKc8vLyzJQpU9rcZsKECSYyMtKMHDnS7z3xsqs9Zk4IM8by+ZQAAOC6dl1f5wMAANhH+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGDV/wMCp7MX1W0yZQAAAABJRU5ErkJggg==",
       "text/plain": [
-       "array([ 78.93258875,  78.93258875,  78.93258875, ..., 151.11306395,\n",
-       "       151.11306395, 151.11306395])"
+       "<Figure size 640x480 with 1 Axes>"
       ]
      },
-     "execution_count": 36,
      "metadata": {},
-     "output_type": "execute_result"
+     "output_type": "display_data"
     }
    ],
    "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]"
+    "_ = plt.hist(\n",
+    "    np.log10(np.maximum(np.sum(params * sj_pt_orders, axis=1), 0.1)),\n",
+    "    np.linspace(-1, 1, 21),\n",
+    "    histtype=\"step\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 75,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ratio_nom_vals = ratio_smeared_lookups[0](flat_subjet_pt, flat_logD, flat_logkt)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 80,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ratio_nom_vals[high_pt_sel] = pt_extrap_vals"
    ]
   },
   {
@@ -1062,6 +1114,7 @@
    "outputs": [],
    "source": [
     "sf_vals = []\n",
+    "\n",
     "# could be parallelised but not sure if memory / time trade-off is worth it\n",
     "for i, ratio_nom_lookup in enumerate(ratio_smeared_lookups):\n",
     "    ratio_nom_vals = ratio_nom_lookup(flat_subjet_pt, flat_logD, flat_logkt)\n",