diff --git a/notebooks/notebook_sampling_ify_u235th_cea_c1.ipynb b/notebooks/notebook_sampling_ify_u235th_cea_c1.ipynb
index cb7913ed..22cfefec 100644
--- a/notebooks/notebook_sampling_ify_u235th_cea_c1.ipynb
+++ b/notebooks/notebook_sampling_ify_u235th_cea_c1.ipynb
@@ -16,6 +16,7 @@
    "outputs": [],
    "source": [
     "import sandy\n",
+    "import matplotlib.pyplot as plt\n",
     "from os.path import join, dirname\n",
     "import pandas as pd\n",
     "import numpy as np\n",
@@ -30,10 +31,18 @@
     "## Extract FYs and covariance data for U235 thermal fission"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "d572351b-4450-4c30-8d81-b3cd1a0cb310",
+   "metadata": {},
+   "source": [
+    "If you don't find these files, download them from the [NEA GITLAB repository](https://git.oecd-nea.org/databank/nds/jeff/jeff-lab/fission-yields-u-235-th)."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "ae456535-ba24-4747-a180-e60926d26fa2",
+   "id": "4e7b4c6d-5c73-4e80-b90a-1ea90137cc87",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -42,16 +51,22 @@
     "file_eval = r\"mixt_cea-jeff33_u_235_th_eval_c1.stn\"\n",
     "file_cov = r\"mixt_cea-jeff33_u_235_th_ind_cov_mat_c1\"\n",
     "tape = sandy.Endf6.from_file(join(path, file_eval))\n",
-    "nfpy = sandy.Fy.from_endf6(tape)\n",
-    "\n",
+    "nfpy = sandy.Fy.from_endf6(tape)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cb8852d0-e541-4580-a4ea-e4139267f513",
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "idx = nfpy.data.query(f\"E==0.0253 & MT==454\").index\n",
     "ifyth = nfpy.data.loc[idx]\n",
-    "IFY_np = np.array([ifyth['FY'].to_numpy()])\n",
     "\n",
-    "cov = pd.read_csv(join(path, file_cov), sep=r\"\\s+\", header=None)\n",
-    "cov = np.divide(cov.data.to_numpy(), IFY_np.T @ IFY_np) # convert to relative terms\n",
-    "cov = pd.DataFrame(cov, index = pd.Index(ifyth.ZAP), columns = pd.Index(ifyth.ZAP))\n",
-    "cov = sandy.CategoryCov(cov)"
+    "acov = pd.read_csv(join(path, file_cov), sep=r\"\\s+\", header=None)                             # absolute covariance matrix\n",
+    "cov = np.divide(acov.values, ifyth.FY.values.reshape(-1, 1) @ ifyth.FY.values.reshape(1, -1)) # convert to relative terms\n",
+    "cov = sandy.CategoryCov(pd.DataFrame(cov, index=ifyth.ZAP.values, columns=ifyth.ZAP.values))"
    ]
   },
   {
@@ -71,7 +86,7 @@
    "source": [
     "seed = random.randrange(2**32 - 1)    # create a seed\n",
     "seed = 1444271359                     # or use always the same\n",
-    "nsmp = 10                             # sample size\n",
+    "nsmp = 1000                           # sample size\n",
     "smp = cov.sampling(nsmp, seed=seed)"
    ]
   },
@@ -121,11 +136,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "path = join(dirname(sandy.__file__), 'appendix', 'fission_yields')\n",
-    "\n",
-    "file_eval = r\"mixt_cea-jeff33_u_235_th_eval_c1.stn\"\n",
-    "tape = sandy.Endf6.from_file(join(path, file_eval))\n",
-    "nfpy = sandy.Fy.from_endf6(tape)\n",
     "idx_ify = nfpy.data.query(f\"E==0.0253 & MT==454\").index\n",
     "idx_cfy = nfpy.data.query(f\"E==0.0253 & MT==459\").index"
    ]
@@ -149,7 +159,7 @@
    "outputs": [],
    "source": [
     "smp_min = 0   # write ENDF-6 file only in the sample range [smp_min, smp_max]\n",
-    "smp_max = 9\n",
+    "smp_max = 99\n",
     "file_template = \"u235_fy_cea_c1_{}.jeff33\"\n",
     "for ismp in range(smp_min, smp_max+1):\n",
     "    file = file_template.format(ismp)\n",
@@ -157,10 +167,62 @@
     "    f.data.loc[idx_ify, \"DFY\"] = f.data.loc[idx_ify, \"FY\"]            # just for me, i copy the original IFYs where uncertainties should be, so i can compare them to the perturbed ones (anyways I don't use uncertainties)\n",
     "    f.data.loc[idx_cfy, \"DFY\"] = f.data.loc[idx_cfy, \"FY\"]            # same but for CFYs\n",
     "    f.data.loc[idx_ify, \"FY\"] *= smp.data[ismp].values                # IMPORTANT, this does not update the CFYs, which in random ENDF-6 file are inconsistent with the perturbed IFYs\n",
-    "    f = f.apply_qmatrix(922350, 0.0253, rdd, keep_fy_index=True)      # Run this if you want to update the CFYs (slower), or else comment it out\n",
+    "    # f = f.apply_qmatrix(922350, 0.0253, rdd, keep_fy_index=True)      # Run this if you want to update the CFYs (slower), or else comment it out\n",
     "    print(f\"writing file '{file}'...\")\n",
     "    f.to_endf6(tape).to_file(file)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "337a51c3-8c9c-46d8-9ab4-46b7d6d1958d",
+   "metadata": {},
+   "source": [
+    "## Compare sample estimates of mass yields with deterministic uncertainty propagation (sandwich formula) "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b6e97e07-dedb-48f1-b59f-4d0ae97e31c7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mfy = {}\n",
+    "for ismp in range(smp_min, smp_max+1):\n",
+    "    mfy[ismp] = sandy.Fy.from_endf6(sandy.Endf6.from_file(f\"u235_fy_cea_c1_{ismp}.jeff33\")).get_mass_yield(zam=922350, e=0.0253)\n",
+    "mfy = pd.DataFrame(mfy).rename_axis(\"SMP\", axis=1)\n",
+    "\n",
+    "S = sandy.Fy(nfpy.data.query(\"MT==454 and E==0.0253\")).get_mass_yield_sensitivity()\n",
+    "C = pd.DataFrame(acov.values, index=ifyth.ZAP.values, columns=ifyth.ZAP.values)\n",
+    "cov_mfy = sandy.CategoryCov(S @ C @ S.T)\n",
+    "\n",
+    "mu = nfpy.get_mass_yield(zam=922350, e=0.0253)\n",
+    "sigma = cov_mfy.get_std()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "261f5b3a-645b-4966-a06b-2d97862e1bf9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axs = plt.subplots(1, 2, figsize=(10, 5), sharex=True)\n",
+    "\n",
+    "ax = axs[0]\n",
+    "ax.errorbar(x=mfy.index, y=mfy.T.mean(), yerr=mfy.T.std(), marker=\"s\", ms=4, linestyle=\"none\", capsize=2, ecolor=\"blue\", color=\"dodgerblue\", label=\"sample estimate\")\n",
+    "ax.errorbar(x=mu.index, y=mu.values, yerr=sigma.values, marker=\"s\", ms=4, linestyle=\"none\", capsize=2, ecolor=\"k\", color=\"k\", alpha=.4, label=\"original data\")\n",
+    "ax.set(ylabel=\"mass yields\", xlabel=\"A\")\n",
+    "\n",
+    "ax = axs[1]\n",
+    "diff_mean = (mfy.T.mean() / mu  - 1) * 100\n",
+    "diff_std = (mfy.T.std() / sigma - 1) * 100\n",
+    "ax.errorbar(x=diff_mean.index, y=diff_mean.values, marker=\"s\", ms=4, linestyle=\"none\", color=\"dodgerblue\", label=\"$\\\\frac{m- \\\\mu}{\\\\mu} \\\\times 100$\")\n",
+    "ax.errorbar(x=diff_std.index, y=diff_std.values, marker=\"s\", ms=4, linestyle=\"none\", color=\"tomato\", label=\"$\\\\frac{s- \\\\sigma}{\\\\sigma} \\\\times 100$\")\n",
+    "ax.set(ylim=[-100, 100], ylabel=\"relative difference / $\\\\%$\", xlabel=\"A\")\n",
+    "ax.legend()\n",
+    "fig.tight_layout()"
+   ]
   }
  ],
  "metadata": {