diff --git a/notebooks/plot_theta2_from_dl3.ipynb b/notebooks/plot_theta2_from_dl3.ipynb
index afc74da5b..afbcf0a72 100644
--- a/notebooks/plot_theta2_from_dl3.ipynb
+++ b/notebooks/plot_theta2_from_dl3.ipynb
@@ -6,7 +6,9 @@
    "source": [
     "# This is an example notebook to plot theta² using DL3 files in custom energy bins\n",
     "\n",
-    "It customizes the Gammapy functions make_theta_squared_table\n",
+    "It customizes the Gammapy function [`make_theta_squared_table`](https://docs.gammapy.org/1.3/api/gammapy.makers.utils.make_theta_squared_table.html).\n",
+    "\n",
+    "Theta² plots can also be done from DL2 files, see [explore_DL2.ipynb notebook](https://github.com/cta-observatory/cta-lstchain/blob/main/notebooks/explore_DL2.ipynb).\n",
     "\n",
     "Content:\n",
     "\n",
@@ -44,42 +46,27 @@
    "outputs": [],
    "source": [
     "%matplotlib inline\n",
+    "\n",
+    "import warnings\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import astropy.units as u\n",
     "import matplotlib.pyplot as plt\n",
     "import matplotlib.style as style\n",
-    "import matplotlib.colors as colors\n",
-    "\n",
     "import numpy as np\n",
-    "import regions\n",
-    "from regions import CircleSkyRegion\n",
-    "from scipy.stats import chi2\n",
-    "from pathlib import Path\n",
-    "import os\n",
-    "from glob import glob"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from gammapy.data import DataStore, EventList\n",
-    "from gammapy.maps import Map, MapAxis, WcsNDMap, WcsGeom\n",
+    "from astropy.coordinates import Angle, SkyCoord\n",
+    "from astropy.table import Table\n",
     "from gammapy.data import DataStore\n",
-    "from gammapy.stats import WStatCountsStatistic\n",
+    "from gammapy.maps import MapAxis\n",
     "from gammapy.maps.utils import edges_from_lo_hi\n",
+    "from gammapy.stats import WStatCountsStatistic\n",
+    "from pyirf.statistics import li_ma_significance\n",
+    "from scipy.stats import chi2\n",
     "\n",
-    "import astropy.units as u\n",
-    "from astropy.table import Table, Column, vstack\n",
-    "from astropy.io import fits\n",
-    "from astropy.time import Time\n",
-    "from astropy.coordinates import SkyCoord, Angle\n",
-    "from astropy.io.fits.verify import VerifyWarning\n",
-    "\n",
-    "import warnings\n",
     "warnings.filterwarnings(\"ignore\")\n",
     "\n",
     "from gammapy import __version__ as gammapy_version\n",
+    "\n",
     "print(f\"Using Gammapy {gammapy_version}\")"
    ]
   },
@@ -97,7 +84,7 @@
    "outputs": [],
    "source": [
     "dl3_path = Path('/fefs/aswg/workspace/analysis-school-2024/DL3/Crab')\n",
-    "datastore = DataStore.from_dir(dl3_path)\n"
+    "datastore = DataStore.from_dir(dl3_path)"
    ]
   },
   {
@@ -109,7 +96,7 @@
     "# Assuming there is only 1 unique source in the DL3 files in this directory\n",
     "obj_name = np.unique(datastore.obs_table[\"OBJECT\"])[0]\n",
     "\n",
-    "target_position = SkyCoord.from_name(\"Crab\")"
+    "target_position = SkyCoord.from_name(obj_name)"
    ]
   },
   {
@@ -118,6 +105,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Show the table of observations\n",
     "datastore.obs_table"
    ]
   },
@@ -155,9 +143,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Selecting only some runs based on previous filters\n",
+    "# Select only some runs based on previous filters\n",
     "obs_table = datastore.obs_table[d_zen[0]*d_obj[0]*d_time[0]]\n",
-    "#obs_table = datastore.obs_table[d_obj[0]*d_time[0]]\n",
     "obs_list = datastore.obs_table[d_zen[0]*d_obj[0]*d_time[0]][\"OBS_ID\"]\n",
     "\n",
     "observations = datastore.get_observations(\n",
@@ -441,18 +428,11 @@
     "            N_off = np.sum(table[\"counts_off\"][0][e_reco_index][:bin_cut])\n",
     "            N_ex = np.sum(table[\"n_sig\"][0][e_reco_index][:bin_cut])\n",
     "            \n",
-    "            # Using Li&Ma's formula as stated above\n",
-    "            if N_off >0:\n",
-    "                S_formula = (N_on - alpha * N_off) / np.sqrt(alpha*(N_off+N_on))\n",
-    "                S_lima = np.sqrt(\n",
-    "                    2 * (\n",
-    "                        N_on * np.log((1+alpha)*N_on /(N_on + N_off)/alpha) + \n",
-    "                        N_off * np.log((1+alpha) * N_off / (N_on + N_off))\n",
-    "                    )\n",
-    "                )\n",
-    "            else:\n",
-    "                S_formula = np.nan\n",
-    "                S_lima = np.nan\n",
+    "            # Using Li&Ma's formula:\n",
+    "            S_lima = li_ma_significance(N_on, N_off, alpha)\n",
+    "            \n",
+    "            # This pyIRF function li_ma_significance is analog to Gammapy's WStatCountsStatistic(N_on, N_off, alpha).sqrt_ts\n",
+    "            # but we keep to be consistent with the way explore_DL2.ipynb notebook calculates the significance.\n",
     "                \n",
     "            # Using gammapy method, but normalizing the binned significance to get total significance in the cut region\n",
     "            table[\"Sig_LiMa\"][0][e_reco_index] = S_lima\n",
@@ -502,7 +482,7 @@
     "                    )\n",
     "                    axs[i,j].grid(ls='dashed')\n",
     "                    axs[i,j].axvline(theta2_edges[bin_cut], ls = '--', color='black', alpha=0.75)\n",
-    "                    axs[i,j].legend(title=f'Significance = {S_lima:.2f} $\\sigma$')\n",
+    "                    axs[i,j].legend(title=f'Significance = {S_lima:.1f} $\\sigma$')\n",
     "\n",
     "            else: # Only excess plot\n",
     "                if not custom_style:\n",
@@ -531,7 +511,7 @@
     "                    axs[i,j].axhline(0, color='darkgray')\n",
     "                    axs[i,j].grid(ls='dashed')\n",
     "                    axs[i,j].axvline(theta2_edges[bin_cut], ls = '--', color='black', alpha=0.75)\n",
-    "                    axs[i,j].legend(title=f'Significance = {S_lima:.2f} $\\sigma$')\n",
+    "                    axs[i,j].legend(title=f'Significance = {S_lima:.1f} $\\sigma$')\n",
     "\n",
     "            axs[i,j].set_ylabel(\"Counts\")\n",
     "            axs[i,j].set_title(f'Energy {ereco_edges[e_reco_index]:.2f}:{ereco_edges[e_reco_index+1]:.2f} TeV')\n",
@@ -576,10 +556,8 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
     "# 5. Plot theta2 distribution for each observation"
    ]
@@ -643,56 +621,9 @@
    "source": [
     "plot_theta_squared_table(theta2_total, theta2_cut, custom_style=custom_style, counts_only=False)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "language": "python",
-   "name": "python3"
-  },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
@@ -703,7 +634,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.11"
+   "version": "3.11.6"
   }
  },
  "nbformat": 4,