Skip to content

Commit

Permalink
use li ma significance formula from pyirf
Browse files Browse the repository at this point in the history
  • Loading branch information
@morcuended
morcuended committed Jan 31, 2025
1 parent 7f3979e commit a64e463
Showing 1 changed file with 28 additions and 97 deletions.
125 changes: 28 additions & 97 deletions notebooks/plot_theta2_from_dl3.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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}\")"
]
},
Expand All @@ -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)"
]
},
{
Expand All @@ -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)"
]
},
{
Expand All @@ -118,6 +105,7 @@
"metadata": {},
"outputs": [],
"source": [
"# Show the table of observations\n",
"datastore.obs_table"
]
},
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -576,10 +556,8 @@
]
},
{
"cell_type": "code",
"execution_count": null,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"# 5. Plot theta2 distribution for each observation"
]
Expand Down Expand Up @@ -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",
Expand All @@ -703,7 +634,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.11"
"version": "3.11.6"
}
},
"nbformat": 4,
Expand Down

0 comments on commit a64e463

Please sign in to comment.