From e6ffa329024e7de1d3d0167f8b378c3a53fe5246 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Tue, 6 Feb 2024 12:28:14 +0000 Subject: [PATCH] add DL1b_to_DL2.ipynb --- notebooks/DL1b_to_DL2.ipynb | 1368 +++++++++++++++++++++++++++++++++++ 1 file changed, 1368 insertions(+) create mode 100644 notebooks/DL1b_to_DL2.ipynb diff --git a/notebooks/DL1b_to_DL2.ipynb b/notebooks/DL1b_to_DL2.ipynb new file mode 100644 index 0000000000..7c998d9992 --- /dev/null +++ b/notebooks/DL1b_to_DL2.ipynb @@ -0,0 +1,1368 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e15f5327-2353-4cde-9f42-ce4a51913ded", + "metadata": {}, + "source": [ + "# DL1b to DL2" + ] + }, + { + "cell_type": "markdown", + "id": "2a0a9d3f-d898-43c9-9248-e6ec8299b151", + "metadata": {}, + "source": [ + "#### Purpose of this notebook\n", + "#### 1. Check DL1b data processed by OSA\n", + "#### 2. Check Random forest models produced by lstmcpipe\n", + "#### 3. DL1b to DL2 step (observation data)\n", + "#### 4. Explore DL2 data\n", + "#### 5. Source-dependent analysis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8d9ecb4-d912-4586-9936-5176897f20dc", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import glob\n", + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import colors\n", + "from astropy.coordinates import SkyCoord\n", + "import astropy.units as u\n", + "from gammapy.stats import WStatCountsStatistic" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cf43ae8-489a-4dae-a96e-627e939c0d36", + "metadata": {}, + "outputs": [], + "source": [ + "import lstchain\n", + "from lstchain.reco.utils import get_effective_time #, radec_to_camera\\\n", + "from ctapipe.containers import EventType\n", + "from lstchain.io.io import (\n", + " dl1_params_lstcam_key,\n", + " dl2_params_lstcam_key,\n", + " dl2_params_src_dep_lstcam_key,\n", + " get_dataset_keys\n", + ")\n", + "from lstchain.io import (\n", + " standard_config,\n", + " srcdep_config,\n", + " get_srcdep_params, \n", + " get_srcdep_assumed_positions\n", + ")\n", + "from lstchain.reco.utils import (\n", + " compute_theta2,\n", + " extract_source_position,\n", + " get_effective_time\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfdaaa7a-7579-4cd8-b5c9-0f0c817d7619", + "metadata": {}, + "outputs": [], + "source": [ + "# Test data are stored under the following directory on La Palma IT cluster\n", + "# If you are locally running this notebook, please change the data directory path|\n", + "data_dir = '/fefs/aswg/workspace/analysis-school-2024/DL1b_to_DL2/'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "616f24a5-ee27-4d2c-b78e-180b23643cfa", + "metadata": {}, + "outputs": [], + "source": [ + "lstchain_dir = os.path.dirname(lstchain.__file__)" + ] + }, + { + "cell_type": "markdown", + "id": "0ca678f9-b3f5-4795-921e-af1b74ec7de9", + "metadata": {}, + "source": [ + "# 1. Check DL1b data processed by OSA" + ] + }, + { + "cell_type": "markdown", + "id": "13450864-f7ee-49c9-bd00-b2d75ee25313", + "metadata": {}, + "source": [ + "#### Before checking the observation data" + ] + }, + { + "cell_type": "markdown", + "id": "74c9d7b2-e0df-45e8-b334-56e73d2e3d64", + "metadata": {}, + "source": [ + "Please check the run list and observation conditions:\n", + "- LST observation schedule: https://www.lst1.iac.es/schedule/\n", + "- ELOG: https://www.lst1.iac.es/elog/LST+commissioning/\n", + "- OSA source catalog: https://www.lst1.iac.es/datacheck/lstosa/LST_source_catalog.html\n", + "- DL1 data check: https://www.lst1.iac.es/datacheck/dl1/\n", + "- Observation overview: https://www.lst1.iac.es/datacheck/observation_overview/\n" + ] + }, + { + "cell_type": "markdown", + "id": "cf110ff0-76cc-4ee2-8670-bc3e64f25b0a", + "metadata": {}, + "source": [ + "#### Check OSA directory structure" + ] + }, + { + "cell_type": "markdown", + "id": "ab83fe6c-5fae-4904-b704-be8391a119cb", + "metadata": {}, + "source": [ + "OSA (Onsite analysis, https://github.com/cta-observatory/lstosa) process starts at 8UTC to create DL1 data. The OSA process is running under `/fefs/aswg/data/real/running_analysis`, then the OSA products are moved to `/fefs/aswg/data/real/DL1` once all the processes are finished. The status of the processing can be seen here (https://www.lst1.iac.es/datacheck/lstosa/sequencer.xhtml). After that, you can see the processed observation runs in the source catalog page (https://www.lst1.iac.es/datacheck/lstosa/LST_source_catalog.html)." + ] + }, + { + "cell_type": "markdown", + "id": "d8e85f19-1fef-4805-84b9-8b7528191093", + "metadata": {}, + "source": [ + "OSA DL1 directory is `/fefs/aswg/data/real/DL1/`. There are three directories under the corresponding DATE and LSTCHAIN_VERSION directory: \n", + "- `interleaved`: interleaved calibration file for Cat-B calibration\n", + "- `muons`: muon analysis products\n", + "- `tailcut84`: DL1b data (still containing DL1a image data) after the standard tailcuts cleaning of picture threshold 8 p.e. and boundary threshold 4 p.e + pedestal-noise-based threshold." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cda90414-c3f2-4e6d-8cfe-8ee847317cfc", + "metadata": {}, + "outputs": [], + "source": [ + "date='20231216'\n", + "lstchain_version='v0.10'\n", + "\n", + "#date=\"20221031\"\n", + "#lstchain_version='v0.9'\n", + "\n", + "!ls /fefs/aswg/data/real/DL1/{date}/{lstchain_version}/ " + ] + }, + { + "cell_type": "markdown", + "id": "cad75032-2fe6-452c-92ef-d3f951cfb379", + "metadata": {}, + "source": [ + "There are two steps before the DL1b data are created:\n", + "1. Apply calibration and create DL1a data (pixel-wise image data) (`lstchain_data_r0_to_dl1`) \n", + "2. Create DL1b data by applying the image clearning with the pedestal noise information derived from the interleaved pedestal events (`lstchain_dl1ab`)\n", + "\n", + "First, DL1a data (pixel-wise image data) are created for each subrun. During this step, pixel-wise pedestal noise values are computed using the interleaved pedestal events and stored in the DL1a data. The pixel-wise pedestal noise values correspond to the night sky background levels on each pixel. If there is a bright star inside the camera field of view, this value on that pixel can be higher than other pixels. In this case, we need to apply higher image clearning level (2.5 x pedestal noise in the standard config) to this pixel to avoid noise events. So, we can apply the image clearning with pedestal threshold. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1d01bfe-2da5-4789-9a49-6bd08e45e937", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"The standard image cleaning configurations:\")\n", + "print(standard_config['tailcuts_clean_with_pedestal_threshold'])\n", + "print(standard_config['dynamic_cleaning'])" + ] + }, + { + "cell_type": "markdown", + "id": "7355fa4b-1ee8-44fb-ba0b-a8f3e38e19b2", + "metadata": {}, + "source": [ + "#### Check DL1b file" + ] + }, + { + "cell_type": "markdown", + "id": "f1490b01-f818-496d-ae52-b7c9a1f59a0b", + "metadata": {}, + "source": [ + "DL1b data are stored under `tailcut84` directory. This process (DL1a to DL1b) is performed for each subrun file. Each subrun file contains 5-6k events (depending on EVB version) corresponding to ~10 seconds observation and 20-min obserbations have ~100-200 subrun files. The file name of DL1 data is `dl1-1.Run{run_id}.{subrun_id}.h5`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7cd17ac-6e4b-4d71-89c9-21a14443e33a", + "metadata": {}, + "outputs": [], + "source": [ + "run_id = '16179'\n", + "subrun_id = '0000'\n", + "\n", + "!ls /fefs/aswg/data/real/DL1/{date}/{lstchain_version}/tailcut84/dl1_LST-1.Run{run_id}.{subrun_id}.h5 " + ] + }, + { + "cell_type": "markdown", + "id": "0ea527b2-7d69-4815-a317-37d6d0ff46f9", + "metadata": {}, + "source": [ + "In this directory, the merged DL1b data (DL1a image data are removed) are also stored for each observation run." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8439636f-b311-4c1e-bce3-2fb2e8a2169d", + "metadata": {}, + "outputs": [], + "source": [ + "!ls /fefs/aswg/data/real/DL1/{date}/{lstchain_version}/tailcut84/dl1_LST-1.Run{run_id}.h5 " + ] + }, + { + "cell_type": "markdown", + "id": "44ce716c-a9ee-4169-a3c5-4007bb005998", + "metadata": {}, + "source": [ + "The dataset key of DL1b level is `/dl1/event/telescope/parameters/LST_LSTCam` (=`dl1_params_lstcam_key`). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "309d7fca-5e1c-4029-a73f-888426dd44ef", + "metadata": {}, + "outputs": [], + "source": [ + "input_dl1_file = '/fefs/aswg/data/real/DL1/20231216/v0.10/tailcut84/dl1_LST-1.Run16179.0000.h5'\n", + "#input_dl1_file = f'{data_dir}/DL1b/dl1_LST-1.Run02008.0000.h5'\n", + "get_dataset_keys(input_dl1_file)" + ] + }, + { + "cell_type": "markdown", + "id": "2664dad6-9023-4a9f-87ec-fa8256db4ca1", + "metadata": {}, + "source": [ + "#### Check DL1b parameters" + ] + }, + { + "cell_type": "markdown", + "id": "890db0fb-e880-4067-9f6f-3287f2f4b9ab", + "metadata": {}, + "source": [ + "DL1b data contains shower image parameters after the image cleaning. Even if a given event doesn't survive after the image cleaning, each event is stored in DL1b data to save `delta_t` (delta time between events) for the effective time calculation later. Those un-survived events after the clearning have `NaN` value in each shower parameter." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5dbf5f6-86e9-4278-a644-a0e6a266c78f", + "metadata": {}, + "outputs": [], + "source": [ + "# Read DL1b data frame\n", + "dl1b_data = pd.read_hdf(input_dl1_file, key=dl1_params_lstcam_key)\n", + "dl1b_data" + ] + }, + { + "cell_type": "markdown", + "id": "aae6b32b-3a9c-4384-8fea-835aabee10a3", + "metadata": {}, + "source": [ + "You can also check how much fraction of the interleaved pedestal events survives after the image cleaning. If the fraction is too high (like >10%), the NSB level is too high and you need to apply the higher image cleaning threshold." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5cac65ab-dcf4-469c-ad2c-11c0d75be805", + "metadata": {}, + "outputs": [], + "source": [ + "# Survived pedestal event ratio\n", + "pedestal_data = dl1b_data[dl1b_data.event_type == EventType.SKY_PEDESTAL.value]\n", + "n_all_ped_events = len(pedestal_data)\n", + "n_survived_ped_events = np.sum(~np.isnan(pedestal_data.intensity))\n", + "survived_ped_events_ratio = n_survived_ped_events/n_all_ped_events \n", + "\n", + "print(f\"Survived pedestal events: {survived_ped_events_ratio*100:.1f}% ({n_survived_ped_events}/{n_all_ped_events} events)\")" + ] + }, + { + "cell_type": "markdown", + "id": "b9541f90-be8d-43f4-84e5-c35afb4fbc3a", + "metadata": {}, + "source": [ + "DL1b parameters contain the shower image, time information, trigger information, pointing information etc. You can check the definition of each parameter in `DL1ParametersContainer`\n", + "https://cta-observatory.github.io/cta-lstchain/api/lstchain.io.lstcontainers.DL1ParametersContainer.html" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "804ae344-1c59-42ac-9557-56db766e8a70", + "metadata": {}, + "outputs": [], + "source": [ + "# DL1 parameters keys\n", + "dl1_keys = dl1b_data.keys()\n", + "dl1_keys" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "631df64b-0d55-462b-8498-7d1c20035772", + "metadata": {}, + "outputs": [], + "source": [ + "# release the memory\n", + "del dl1b_data" + ] + }, + { + "cell_type": "markdown", + "id": "1cd79ebe-b492-45d8-b286-779c008172ab", + "metadata": {}, + "source": [ + "# 2. Check random forest models produced by lstmcpipe" + ] + }, + { + "cell_type": "markdown", + "id": "c06a4287-addb-46af-914a-e5611562b79a", + "metadata": {}, + "source": [ + "Random forest (RF) models are created by `lstmcpipe` (https://github.com/cta-observatory/lstmcpipe). The `lstmcpipe` products are stored under the following directory. `PROD_ID` is MC production-based directory (e.g. `20240131_allsky_v0.10.5_all_dec_base`), `DEC` is for the corresponding declination line (e.g. `dec_2276`).\n", + "- `/fefs/aswg/data/mc/DL1/AllSky/PROD_ID/TrainingDataset/DEC/{GammaDiffuse, Protons}`: Training DL1 MC\n", + "- `/fefs/aswg/data/mc/DL1/AllSky/PROD_ID/TestingDataset/DEC/`: Test DL1 MC\n", + "- `/fefs/aswg/data/models/AllSky/PROD_ID/DEC`: RF models\n", + "- `/fefs/aswg/data/mc/DL2/AllSky/PROD_ID/TestingDataset/DEC`: Test DL2 MC" + ] + }, + { + "cell_type": "markdown", + "id": "46b9fcf4-ea0b-41be-b7b3-56b4b76d3eda", + "metadata": {}, + "source": [ + "#### Check declination-based random forest models" + ] + }, + { + "cell_type": "markdown", + "id": "d7d13066-904e-4558-ad4a-2806479bcb7a", + "metadata": {}, + "source": [ + "There are several RF models for the corresponding declination line. You need to choose the best RF models of which declination is close to the target source declination." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e65dbf2-ab1e-46f5-91d3-edde53b55ed4", + "metadata": {}, + "outputs": [], + "source": [ + "production_id = \"20240131_allsky_v0.10.5_all_dec_base\"\n", + "declination_dir = \"dec_2276\" # DEC=22.76 deg\n", + "\n", + "!ls /fefs/aswg/data/models/AllSky/{production_id}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a546bd6-1c35-4b0c-82a0-58f6aef55e3b", + "metadata": {}, + "outputs": [], + "source": [ + "def get_alt_az_for_target_dec(dec_target):\n", + " \n", + " latitude = 28.76139 * u.deg # LST-1\n", + " local_hour_angle = np.linspace(-180, 180, 200) * u.deg\n", + " alt = np.arcsin(np.sin(dec_target) * np.sin(latitude) + np.cos(dec_target) * np.cos(latitude) * np.cos(local_hour_angle))\n", + " az = np.arcsin(- np.sin(local_hour_angle) * np.cos(dec_target)/np.cos(alt))\n", + " zd = 90 * u.deg - alt\n", + " \n", + " culmination_zd = latitude - dec_target\n", + " print(f\"Target declination: {dec_target.degree} deg, Culmination is Zd = {np.abs(culmination_zd):.2f}\")\n", + "\n", + " select_zd90 = zd < 90 * u.deg\n", + " zd_cut = zd[select_zd90]\n", + " az_cut = az[select_zd90]\n", + " \n", + " zd_az90 = zd_cut[np.argmin(np.abs(az_cut - 90 * u.deg))]\n", + " \n", + " # azimuth is coputed within [-pi, pi], so it is needed to convert to the range within [pi, 1.5*pi] for south pointing points\n", + " if culmination_zd > 0 * u.deg:\n", + " select_south = zd_cut < zd_az90\n", + " az_cut[select_south] = (az_cut[select_south] + 180 * u.deg)\n", + " az_cut[select_south] = az_cut[select_south][::-1]\n", + " \n", + " return az_cut, zd_cut" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e32f439-b30a-4a51-990d-bf19a3c84a6d", + "metadata": {}, + "outputs": [], + "source": [ + "dec_dir_list = sorted(glob.glob(f\"/fefs/aswg/data/mc/DL1/AllSky/{production_id}/TrainingDataset/dec_*\"))\n", + "\n", + "fig = plt.figure(figsize=(8, 6))\n", + "ax = fig.add_subplot(projection='polar')\n", + "\n", + "for dec_dir in dec_dir_list:\n", + " dec_dir_basename = os.path.basename(dec_dir)\n", + "\n", + " if 'dec_min' in dec_dir_basename:\n", + " dec_mc = float(dec_dir_basename.replace(\"dec_min_\", \"\")) * 1e-2 * -1\n", + " else:\n", + " dec_mc = float(dec_dir_basename.replace(\"dec_\", \"\")) * 1e-2\n", + "\n", + " node_dir_list = sorted(glob.glob(f\"{dec_dir}/GammaDiffuse/node*\"))\n", + "\n", + " mc_zenith_list = []\n", + " mc_azimuth_list = []\n", + "\n", + " for node_dir in node_dir_list:\n", + " node_dir_basename_split = os.path.basename(node_dir).split('_')\n", + " zd_tel = np.deg2rad(float(node_dir_basename_split[3]))\n", + " az_tel = np.deg2rad(float(node_dir_basename_split[5]))\n", + " \n", + " mc_zenith_list.append(zd_tel)\n", + " mc_azimuth_list.append(az_tel)\n", + " \n", + " plt.scatter(mc_azimuth_list, np.rad2deg(mc_zenith_list), label=dec_dir_basename)\n", + "\n", + "az, zd = get_alt_az_for_target_dec(dec_target)\n", + "#plt.plot(az, zd, label=f\"Target: {dec_target.to_value(u.degree)} deg\")\n", + "plt.plot(az, zd, label=f\"Target: {dec_target.degree} deg\")\n", + "\n", + "plt.grid(ls='--')\n", + "plt.legend(loc='upper left', bbox_to_anchor=(1, 1))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "6078aaf0-3c36-4d86-b6de-7fe63e1be88c", + "metadata": {}, + "source": [ + "#### Check RF models " + ] + }, + { + "cell_type": "markdown", + "id": "66fb6e9a-c98c-459a-a78e-a98e98979942", + "metadata": {}, + "source": [ + "Trained model files (scikit-learn) are stored in this directory. There is also a configuration json file for the RF training. You need to use the same configurations (input parameters etc) for DL2 data production.\n", + "- `reg_energy.sav`: energy regression (`reco_energy`)\n", + "- `cls_disp_sign.sav`: arrival direction reconstruction (`disp_sign`)\n", + "- `reg_disp_norm.sav`: arrival direction reconstruction (`disp_norm`)\n", + "- `cls_gh.sav`: particle type classification (`gammaness`)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f133fe2-d6ee-470e-8986-70547740409e", + "metadata": {}, + "outputs": [], + "source": [ + "!ls /fefs/aswg/data/models/AllSky/{production_id}/{declination_dir}" + ] + }, + { + "cell_type": "markdown", + "id": "fd68a6f0", + "metadata": {}, + "source": [ + "# 3. DL1b to DL2 step (observation data)" + ] + }, + { + "cell_type": "markdown", + "id": "a4b35c4c-2621-4967-ab96-7c1849f47137", + "metadata": {}, + "source": [ + "#### Script usage" + ] + }, + { + "cell_type": "markdown", + "id": "57b6a9c1", + "metadata": {}, + "source": [ + "You can check the usage of the script: `lstchain_dl1_to_dl2 -h`" + ] + }, + { + "cell_type": "markdown", + "id": "84c8087b", + "metadata": {}, + "source": [ + "```\n", + "usage: lstchain_dl1_to_dl2 [-h] --input-files INPUT_FILES [INPUT_FILES ...] [--path-models PATH_MODELS] [--output-dir OUTPUT_DIR]\n", + " [--config CONFIG_FILE]\n", + "\n", + "Run the DL1 to DL2 step: Pipeline for the reconstruction of Energy, disp and gamma/hadron separation of events stored in a DL1 file. It\n", + "takes DL1 file(s) and trained Random Forests as input and outputs DL2 data file(s). Run lstchain_dl1_to_dl2 --help to see the options.\n", + "\n", + "options:\n", + " -h, --help show this help message and exit\n", + " --input-files INPUT_FILES [INPUT_FILES ...], -f INPUT_FILES [INPUT_FILES ...]\n", + " Path (or list of paths) to a DL1 HDF5 file\n", + " --path-models PATH_MODELS, -p PATH_MODELS\n", + " Path where to find the trained RF\n", + " --output-dir OUTPUT_DIR, -o OUTPUT_DIR\n", + " Path where to store the reco dl2 events\n", + " --config CONFIG_FILE, -c CONFIG_FILE\n", + " Path to a configuration file. If none is given, a standard configuration is applied\n", + "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "2936ccc1-2e49-426a-bf80-78d36a9510a5", + "metadata": {}, + "source": [ + "#### Launch the script with small test data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9fe03edf-65d2-4a6a-8022-e4abe16b643b", + "metadata": {}, + "outputs": [], + "source": [ + "input_dl1_file = f'{data_dir}/DL1b/dl1_LST-1.Run02008.0000.h5'\n", + "model_path = f'{data_dir}/models/srcindep/'\n", + "output_dir = './DL2_ouput_srcindep'\n", + "os.makedirs(output_dir, exist_ok=True)\n", + "config_file = f'{lstchain_dir}/data/lstchain_standard_config.json'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f0231f63-0b04-41a4-bbb8-5c781fe94fbf", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "!lstchain_dl1_to_dl2 -f {input_dl1_file} -p {model_path} -o {output_dir} -c {config_file}" + ] + }, + { + "cell_type": "markdown", + "id": "a34fa262", + "metadata": {}, + "source": [ + "# 4. Explore DL2 data" + ] + }, + { + "cell_type": "markdown", + "id": "b5ae6f23-93a4-43e8-b2af-c07f03ffa2d0", + "metadata": {}, + "source": [ + "#### Check DL2 parameters" + ] + }, + { + "cell_type": "markdown", + "id": "6d9fe844-afc1-4fa1-94fa-6b71db06831c", + "metadata": {}, + "source": [ + "Here let's use the already-processed DL2 data (run 16179, merge of 50 subruns)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6599be13", + "metadata": {}, + "outputs": [], + "source": [ + "dl2_file_path = f'{data_dir}/DL2/srcindep/dl2_LST-1.Run16179_small.h5'" + ] + }, + { + "cell_type": "markdown", + "id": "6aaf0a17-7306-467f-8684-3c4546ff673d", + "metadata": {}, + "source": [ + "You can see `/dl2/event/telescope/parameters/LST_LSTCam` (=`dl2_params_lstcam_key`) are added in dataset keys." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f26f0bf8-07b8-4f3d-b97d-150531295e4b", + "metadata": {}, + "outputs": [], + "source": [ + "get_dataset_keys(dl2_file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7331b75", + "metadata": {}, + "outputs": [], + "source": [ + "data = pd.read_hdf(dl2_file_path, key=dl2_params_lstcam_key)\n", + "data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "108d1f82-82de-4088-ac9c-2c2eb0c49a0a", + "metadata": {}, + "outputs": [], + "source": [ + "# Added DL2 keys in addition to DL1 keys\n", + "dl2_keys = data.keys()\n", + "dl2_keys.drop(dl1_keys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05b66475-7ded-448d-b654-6afccc653e24", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(5, 4))\n", + "\n", + "m = ax.hist2d(data.log_reco_energy, data.gammaness,\n", + " range=[[-2, 2], [0, 1]], bins=100, cmin=1, norm=colors.LogNorm())\n", + "ax.set_xlabel('log10(reco energy [TeV])')\n", + "ax.set_ylabel('gammaness')\n", + "ax.grid(ls='--')\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2af4400-f6f5-4238-8bd7-ecf92ffae419", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(5, 4))\n", + "\n", + "intensity_cut_edges = [0, 80, 200, 800, 3200]\n", + "for i in range(len(intensity_cut_edges) - 1):\n", + " intensity_cut = intensity_cut_edges[i:i+2]\n", + "\n", + " ax.hist(data.gammaness[(data.intensity > intensity_cut[0]) & (data.intensity <= intensity_cut[1])], \n", + " range=[0, 1], bins=100, histtype='step', label=f'{intensity_cut[0]} < intensity < {intensity_cut[1]}')\n", + "ax.set_xlabel('gammaness')\n", + "ax.legend(fontsize=8)\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d82829e6-199c-4ebc-b446-3aec992aaaf0", + "metadata": {}, + "source": [ + "#### Effective time calculation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3eb9bcdb-d9ab-4ff5-b793-6e2851a1a3d9", + "metadata": {}, + "outputs": [], + "source": [ + "average_rate = 1/np.mean(data.delta_t)\n", + "deadtime = np.min(data.delta_t[data.delta_t>0]) * 1e6" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "829134c7-24c3-4804-8af6-2dc1259fa926", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(5, 4))\n", + "ax.hist(data.delta_t*1e6, range=[0, 1e3], bins=50,\n", + " label=f'Rate={average_rate:.0f} Hz\\nDead time = {deadtime:.2f} us')\n", + "ax.set_yscale('log')\n", + "ax.set_xlabel('delta_t [us]')\n", + "ax.legend()\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef8f0209-601f-418e-9e38-0311b4e290c9", + "metadata": {}, + "outputs": [], + "source": [ + "t_eff, t_elapsed = get_effective_time(data)\n", + "#obstime_real = get_effective_time(data)[0]\n", + "print(f\"Effective time {t_eff.to(u.min):.1f}\")\n", + "print(f\"Elapsed time {t_elapsed.to(u.min):.1f}\")\n", + "print(f\"Deadtime: {100 - t_eff/t_elapsed*100:.1f}%\")" + ] + }, + { + "cell_type": "markdown", + "id": "5080e6c0-d6d5-4022-90da-9af36253ba87", + "metadata": {}, + "source": [ + "#### Data selection" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f13d680", + "metadata": {}, + "outputs": [], + "source": [ + "gammaness_cut = 0.7\n", + "intensity_cut = 50 #p.e." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff8ce148-e967-4422-b56a-5564964b237b", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(5, 4))\n", + "\n", + "ax.hist(np.log10(data.intensity[data.event_type == EventType.SUBARRAY.value]), \n", + " range=[1, 6], bins=100, histtype='step', label='SUBARRAY')\n", + "ax.hist(np.log10(data.intensity[data.event_type == EventType.FLATFIELD.value]), \n", + " range=[1, 6], bins=100, histtype='step', label='FLATFIELD')\n", + "ax.hist(np.log10(data.intensity[data.event_type == EventType.SKY_PEDESTAL.value]), \n", + " range=[1, 6], bins=100, histtype='step', label='SKY_PEDESTAL')\n", + "\n", + "ax.axvline(np.log10(intensity_cut), ls='--', label=f'{intensity_cut} p.e.')\n", + "ax.set_yscale('log')\n", + "ax.set_xlabel('log10(intensity [p.e.])')\n", + "ax.legend()\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "bfdca4e1-11b6-40c1-b773-94b2aeb3200b", + "metadata": {}, + "source": [ + "## 5. Create theta2 plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "599f669d", + "metadata": {}, + "outputs": [], + "source": [ + "condition = (data.gammaness > gammaness_cut) \\\n", + " & (data.intensity > intensity_cut) \\\n", + " & (data.event_type == EventType.SUBARRAY.value)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04b035d5", + "metadata": {}, + "outputs": [], + "source": [ + "selected_data = data[condition] \n", + "print(f'Total number of original events: {len(data)}')\n", + "print(f'Total number of selected events: {len(selected_data)}')\n", + "\n", + "del data" + ] + }, + { + "cell_type": "markdown", + "id": "ee48be9b-d536-4b2d-bd38-10fed082b93d", + "metadata": {}, + "source": [ + "#### Define theta2 range (ON and normalization region)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82de874a-c22e-43d2-a176-8e2dcc74a8b4", + "metadata": {}, + "outputs": [], + "source": [ + "THETA2_GLOBAL_CUT = 0.04\n", + "theta2_range = (0, 1)\n", + "norm_theta2_range = (0.5, 1)" + ] + }, + { + "cell_type": "markdown", + "id": "0bf5160e-3b36-4f80-8831-5a673090e49a", + "metadata": {}, + "source": [ + "#### True source position" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0dc918d2", + "metadata": {}, + "outputs": [], + "source": [ + "#If its a catalogued source, like the Crab, you can use this lstchain function\n", + "true_source_position = extract_source_position(selected_data, 'Crab')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "005b57ec-3d9a-47d1-b318-ae50888f0ccc", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "from lstchain.reco.utils import radec_to_camera\n", + "\n", + "def extract_source_position_from_coord(\n", + " data, coord, equivalent_focal_length=28 * u.m\n", + "):\n", + " obstime = pd.to_datetime(data[\"dragon_time\"], unit=\"s\")\n", + " pointing_alt = u.Quantity(data[\"alt_tel\"], u.rad, copy=False)\n", + " pointing_az = u.Quantity(data[\"az_tel\"], u.rad, copy=False)\n", + " source_pos_camera = radec_to_camera(\n", + " coord,\n", + " obstime,\n", + " pointing_alt,\n", + " pointing_az,\n", + " focal=equivalent_focal_length,\n", + " )\n", + " source_position = [source_pos_camera.x, source_pos_camera.y]\n", + " return source_position\n", + "\n", + "coordinates = SkyCoord(ra = 83.6287 * u.degree, dec = 22.0147 * u.degree, frame='icrs')\n", + "true_source_position = extract_source_position_from_coord(selected_data, coordinates)\n", + "\"\"\"" + ] + }, + { + "cell_type": "markdown", + "id": "9d828b23-c3a4-47e9-b2a9-73a31e4b195f", + "metadata": {}, + "source": [ + "#### Off source position" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "747d9f8d-8617-4ff4-b154-dd8b346f8dee", + "metadata": {}, + "outputs": [], + "source": [ + "off_source_position = [element * -1 for element in true_source_position]" + ] + }, + { + "cell_type": "markdown", + "id": "eb68fbb5-c0f1-4cc2-bc80-a037cb10ff59", + "metadata": {}, + "source": [ + "#### Compute theta2 for ON/OFF position" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ebd414f", + "metadata": {}, + "outputs": [], + "source": [ + "theta2_on = np.array(compute_theta2(selected_data, true_source_position))\n", + "theta2_off = np.array(compute_theta2(selected_data, off_source_position))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41405d22-c831-46f1-934f-6b9ebf6eca8c", + "metadata": {}, + "outputs": [], + "source": [ + "data_on = theta2_on\n", + "data_off = theta2_off\n", + "signal_cut = THETA2_GLOBAL_CUT\n", + "x_range = theta2_range\n", + "norm_x_range = norm_theta2_range \n", + "x_label = \"$\\\\theta^{2} [deg^{2}]$\"\n", + "nbins=round((x_range[1]/signal_cut)*2) # Make the histogram so there are only two bins before the theta2 cut\n", + "#nbins = 100\n", + "nbins=round((theta2_range[1]/THETA2_GLOBAL_CUT)*2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4949dfa", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_on_off_and_excess(data_on, data_off, signal_cut, nbins, x_range, norm_x_range, x_label,t_eff, title=None):\n", + "\n", + " hist_on, bin_edges_on = np.histogram(data_on, density=False, bins=nbins, range=x_range)\n", + " hist_off, bin_edges_off = np.histogram(data_off, density=False, bins=nbins, range=x_range)\n", + "\n", + " bin_width = bin_edges_on[1]-bin_edges_off[0]\n", + " bin_center = bin_edges_on[:-1]+(bin_width/2)\n", + "\n", + " N_on = np.sum(hist_on[bin_edges_on[1:]<=signal_cut])\n", + " N_off = np.sum(hist_off[bin_edges_off[1:]<=signal_cut])\n", + " \n", + " idx_min = (np.abs(bin_edges_on - norm_x_range[0])).argmin()\n", + " idx_max = (np.abs(bin_edges_on - norm_x_range[1])).argmin()\n", + " \n", + " Non_norm = np.sum(hist_on[idx_min:idx_max])\n", + " Noff_norm = np.sum(hist_off[idx_min:idx_max])\n", + " \n", + " alpha = Noff_norm / Non_norm\n", + " \n", + " stat = WStatCountsStatistic(n_on=N_on, n_off=N_off, alpha=alpha)\n", + " significance_lima = stat.sqrt_ts\n", + "\n", + " excess = hist_on - hist_off\n", + " excess_err = np.sqrt(hist_on + hist_off)\n", + "\n", + " props = dict(boxstyle='round', facecolor='wheat', alpha=0.95)\n", + " textstr = r'N$_{{\\rm on}}$ = {:.0f}, N$_{{\\rm off}}$ = {:.0f} '\\\n", + " f'\\n'\\\n", + " r'Time = {:.1f}'\\\n", + " f'\\n'\\\n", + " r'LiMa Significance = {:.1f} $\\sigma$ '.format(N_on,\n", + " N_off,\n", + " t_eff.to(u.h),\n", + " significance_lima)\n", + "\n", + " fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n", + " if title is not None:\n", + " fig.suptitle(title)\n", + " \n", + " #ax[0].hist(bin_center, weights=hist_on, bins=nbins, range=theta2_range, histtype='step', lw=2, color='tab:blue')\n", + " ax[0].errorbar(bin_center, hist_on, yerr=np.sqrt(hist_on), fmt='o', label='ON data', color='tab:blue')\n", + " ax[0].hist(bin_center, weights=hist_off, bins=nbins, range=x_range, histtype='step', lw=2, color='tab:orange', zorder=1)\n", + " ax[0].errorbar(bin_center, hist_off, yerr=np.sqrt(hist_off),fmt='.',label='Background', color='tab:orange', zorder=1)\n", + " ax[0].grid(ls='dashed')\n", + " ax[0].axvline(signal_cut, color='black',ls='--',alpha=0.75)\n", + " ax[0].set_xlabel(x_label)\n", + " ax[0].set_ylabel(\"Counts\")\n", + " ax[0].legend(loc='lower right')\n", + " ax[0].text(0.45, 0.95, textstr, transform=ax[0].transAxes, verticalalignment='top', bbox=props)\n", + " \n", + " ax[1].errorbar(bin_center, excess, yerr=excess_err,fmt='o',color='forestgreen',label='Excess counts')\n", + " ax[1].bar(bin_edges_on[:-1], excess, width = bin_width, align='edge', color='limegreen',alpha=0.5)\n", + " ax[1].axhline(0, color='darkgray')\n", + " ax[1].grid()\n", + " ax[1].axvline(signal_cut, color='black', ls='--', alpha=0.75)\n", + " ax[1].grid(ls='dashed')\n", + " ax[1].set_xlabel(x_label)\n", + " ax[1].set_ylabel(\"Counts\")\n", + " ax[1].legend(title=f'Significance = {significance_lima:.1f} $\\sigma$')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f248b81-0354-4a18-9b56-25049ca99d55", + "metadata": {}, + "outputs": [], + "source": [ + "plot_on_off_and_excess(data_on, data_off, signal_cut, nbins, x_range, norm_x_range, x_label, t_eff)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "149c0913", + "metadata": {}, + "outputs": [], + "source": [ + "log_reco_energy = np.array(selected_data.log_reco_energy)\n", + "emin=0.01 * u.TeV\n", + "emax=10 * u.TeV\n", + "n_bins_energy=2\n", + "log_energy = np.linspace(np.log10(emin.to_value()),\n", + " np.log10(emax.to_value()),\n", + " n_bins_energy + 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d92515f", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "for i in range(n_bins_energy):\n", + " condition_energy_bin= (selected_data.log_reco_energy < log_energy[i+1]) & (selected_data.log_reco_energy >= log_energy[i])\n", + "\n", + " title_txt = f'{10**log_energy[i]:.2f} TeV to {10**log_energy[i+1]:.2f} TeV'\n", + " plot_on_off_and_excess(data_on[condition_energy_bin], data_off[condition_energy_bin], signal_cut, nbins, x_range, norm_x_range, x_label, t_eff, title_txt)\n" + ] + }, + { + "cell_type": "markdown", + "id": "3f2ef57a", + "metadata": {}, + "source": [ + "# 5. source-dependent DL2" + ] + }, + { + "cell_type": "markdown", + "id": "1b7006ea-22bf-45c5-8e30-1c1ab10f8679", + "metadata": {}, + "source": [ + "Let's analyze the data in a different method: source-dependent approach. For the source-dependent analysis, several parameters are computed with an assumption of the expected source position (e.g. `alpha`, `dist`). Thanks to a pri-ori information, the reconstruction performance (e.g. energy bias) is better than the standard approach although there is not so much difference in the sensitivity curve because background condition (S/N>5%) affected the low-energy region.\n", + "However, this approach is helpful for the pulsar analysis.\n", + "Please note that this approach can be used only for the point source." + ] + }, + { + "cell_type": "markdown", + "id": "ae5fe342-4dbc-487c-a0e2-1e8e03f9764e", + "metadata": {}, + "source": [ + "You need to activate `source_dependent` as `True` in the config file. Then you also need to set `source_name` (or `source_ra` and `source_dec`) and `n_off`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e89506f-58f9-4354-8c31-2944d2740c38", + "metadata": {}, + "outputs": [], + "source": [ + "srcdep_related_config_keys = [\n", + " 'source_dependent',\n", + " 'energy_regression_features',\n", + " 'particle_classification_features',\n", + " 'n_off_wobble',\n", + " 'source_name',\n", + " 'source_ra',\n", + " 'source_dec'\n", + "]\n", + "\n", + "for key in srcdep_related_config_keys:\n", + " print(f\"{key} : {srcdep_config[key]}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f801b06f-5795-4f49-9fab-34016bc55266", + "metadata": {}, + "outputs": [], + "source": [ + "model_path = f'{data_dir}/models/srcdep'\n", + "output_dir = \"DL2_ouput_srcdep\"\n", + "os.makedirs(output_dir, exist_ok=True)\n", + "config_file = f'{lstchain_dir}/data/lstchain_src_dep_config.json'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "417b546d-9443-4c9f-992c-de3771350937", + "metadata": {}, + "outputs": [], + "source": [ + "!lstchain_dl1_to_dl2 -f {input_dl1_file} -p {model_path} -o {output_dir} -c {config_file}" + ] + }, + { + "cell_type": "markdown", + "id": "dcb923ba-131a-42a5-9e19-34742eec33a7", + "metadata": {}, + "source": [ + "## Explore source-dependent DL2 data" + ] + }, + { + "cell_type": "markdown", + "id": "eab794cb-31d5-48f8-9bc7-76fc57d173e7", + "metadata": {}, + "source": [ + "#### Check DL2 data" + ] + }, + { + "cell_type": "markdown", + "id": "2c9ac4cd-7d99-4c4d-ad81-d5feed88a164", + "metadata": {}, + "source": [ + "Here let's use the already-processed DL2 data (run 16179, merge of 50 subruns)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b304cf-53a2-489a-849b-d7346236f2db", + "metadata": {}, + "outputs": [], + "source": [ + "dl2_file_path = f'{data_dir}/DL2/srcdep/dl2_LST-1.Run16179_small.h5'" + ] + }, + { + "cell_type": "markdown", + "id": "fc69ad0d-10b0-48f0-99fa-e022aad2949f", + "metadata": {}, + "source": [ + "DL1 dataset is replaced with `/dl2/event/telescope/parameters/LST_LSTCam`. In addition, source-dependent DL2 parameters are stored in `/dl2/event/telescope/parameters_src_dependent/LST_LSTCam`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14d43ad8-7a5c-44f2-9f90-dfaa3c9e126f", + "metadata": {}, + "outputs": [], + "source": [ + "get_dataset_keys(dl2_file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b72b655a-5f5f-4e3c-b2d3-ce425da6f0c6", + "metadata": {}, + "outputs": [], + "source": [ + "data = pd.read_hdf(dl2_file_path, key=dl2_params_lstcam_key)\n", + "data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "066369cf-4ec9-452f-be70-fa7eddfa99de", + "metadata": {}, + "outputs": [], + "source": [ + "data_srcdep = get_srcdep_params(dl2_file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a264464-3d93-4d98-aa92-ff0f95363f18", + "metadata": {}, + "outputs": [], + "source": [ + "get_srcdep_assumed_positions(dl2_file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02b795e4-8e05-4845-af8a-156b4847837c", + "metadata": {}, + "outputs": [], + "source": [ + "data_srcdep['on'].head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa247a85-c54a-440e-9668-256cf166e8b6", + "metadata": {}, + "outputs": [], + "source": [ + "data_srcdep['off_180'].head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cac58059-e16e-46ea-9690-2e4c05c7dc86", + "metadata": {}, + "outputs": [], + "source": [ + "data_srcdep['on'].keys()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee9c6412-56a0-4694-87c7-1ba93d425e34", + "metadata": {}, + "outputs": [], + "source": [ + "t_eff, t_elapsed = get_effective_time(data)\n", + "\n", + "print(f\"Effective time {t_eff.to(u.min):.1f}, Elapsed time {t_elapsed.to(u.min):.1f}\")\n", + "print(f\"Deadtime: {100 - t_eff/t_elapsed*100:.1f}%\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66bb3b3c-9180-43d3-b29d-90cc8a212faa", + "metadata": {}, + "outputs": [], + "source": [ + "gammaness_cut = 0.7\n", + "intensity_cut = 50 #p.e." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca859c86-4018-4ccf-89fb-40f0d5f7b03b", + "metadata": {}, + "outputs": [], + "source": [ + "condition_on = (data_srcdep['on'].gammaness > gammaness_cut) \\\n", + " & (data.intensity > intensity_cut) \\\n", + " & (data.event_type == EventType.SUBARRAY.value)\n", + "\n", + "condition_off = (data_srcdep['off_180'].gammaness > gammaness_cut) \\\n", + " & (data.intensity > intensity_cut) \\\n", + " & (data.event_type == EventType.SUBARRAY.value)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a76ec44-4e8c-4ad9-a851-f413db6c8e59", + "metadata": {}, + "outputs": [], + "source": [ + "selected_data_on = data_srcdep['on'][condition_on]\n", + "selected_data_off = data_srcdep['off_180'][condition_off]\n", + "\n", + "del data\n", + "del data_srcdep" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32219418-4f09-45b4-b43d-54922d25a531", + "metadata": {}, + "outputs": [], + "source": [ + "# alpha cut\n", + "signal_cut = 10 \n", + "x_range = (0, 90)\n", + "norm_x_range = (50, 90)\n", + "x_label = r\"$\\alpha$ [deg]\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9eb5ee9-2096-4fff-8838-ec285cc2bc41", + "metadata": {}, + "outputs": [], + "source": [ + "data_on = selected_data_on.alpha\n", + "data_off = selected_data_off.alpha\n", + "\n", + "plot_on_off_and_excess(data_on, data_off, signal_cut, nbins, x_range, norm_x_range, x_label, t_eff)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "729b249a-d233-48bc-9da4-a34a44a72dac", + "metadata": {}, + "outputs": [], + "source": [ + "log_reco_energy = np.array(selected_data.log_reco_energy)\n", + "emin=0.01 * u.TeV\n", + "emax=10 * u.TeV\n", + "n_bins_energy=2\n", + "log_energy = np.linspace(np.log10(emin.to_value()),\n", + " np.log10(emax.to_value()),\n", + " n_bins_energy + 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5617c658-aceb-4988-8640-7e009b741689", + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(n_bins_energy):\n", + " condition_on_energy_bin= (selected_data_on.log_reco_energy < log_energy[i+1]) & (selected_data_on.log_reco_energy >= log_energy[i])\n", + " condition_off_energy_bin= (selected_data_off.log_reco_energy < log_energy[i+1]) & (selected_data_off.log_reco_energy >= log_energy[i])\n", + "\n", + " title_txt = f'{10**log_energy[i]:.2f} TeV to {10**log_energy[i+1]:.2f} TeV'\n", + " plot_on_off_and_excess(data_on[condition_on_energy_bin], data_off[condition_off_energy_bin], signal_cut, nbins, x_range, norm_x_range, x_label, t_eff, title_txt)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}