Step1_DL2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59f4400b",
   "metadata": {},
   "source": [
    "# Step 1: DL2 (data and MCs)\n",
    "\n",
    "In this notebook, we will run the magic-cta-pipe (MCP) scripts on a small DL1 data sample. Due to time constraints, it is unfeasible to run the pipeline on a full dataset needed to produce meaningful plots, so we have  provided a complete dataset to get 'nice' plots and a few MCs and *.h5* data files to try to run the pipeline.  \n",
    "\n",
    "\n",
    "### First we import a few basic modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba1c0aab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import logging\n",
    "import os\n",
    "import subprocess\n",
    "import sys\n",
    "log = logging.getLogger()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9172ebf9",
   "metadata": {},
   "source": [
    "### Converting DL1 to DL2 \n",
    "\n",
    "Here we will use the RFs to convert DL1 into DL2 (both MCs and real data). DL2 data contain events whose energy, direction and gammaness are evaluated by the telescope-wise RFs, such that, if the three telescopes see the same event, this event will have three energy/direction/gammaness values, one for each telescope.  \n",
    "\n",
    "Since the events are separated into `combo_types`, the script uses the appropriate RFs (i.e.: right telescope, right combo_type).\n",
    "\n",
    "The reconstructed arrival direction is defined by the MCP script using the MARS-like DISP method, which looks for the minimum angular distance between all the head and tail candidates, as shown in the figure below:\n",
    "\n",
    "![head](./figures/head_tail.png)\n",
    "\n",
    "\n",
    "We can run the script directly in the terminal, or here in the jupyter notebook, with:\n",
    "\n",
    ">$ python lst1_magic_dl1_stereo_to_dl2.py -d Path1 -r Path2 -o Path3\n",
    "\n",
    "with the following options:  \n",
    "\n",
    "-d: input file (DL1 stereo, MCs or real data; test sample gammas are ring-wobble ($0.4^{\\circ}$), test sample protons are diffuse)\n",
    "\n",
    "-r: directory where you stored your RFs\n",
    "\n",
    "-o: output directory, to store DL2 files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebe43b6a",
   "metadata": {},
   "source": [
    "#### Let's start with the MC\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2f26fbfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_dl2_g=('...../DL2/gamma')\n",
    "dir_dl2_p=('...../DL2/proton')\n",
    "dir_rf=('/fefs/aswg/workspace/analysis-school-2024/MCP/dl2/rf/')\n",
    "#dir_dl2_g=('/fefs/aswg/workspace/elisa.visentin/sw_school_2024/DL2/gamma')\n",
    "#dir_dl2_p=('/fefs/aswg/workspace/elisa.visentin/sw_school_2024/DL2/proton')\n",
    "os.makedirs(dir_dl2_g)\n",
    "os.makedirs(dir_dl2_p)\n",
    "\n",
    "f_g=open(f'{dir_dl2_g}/DL2_mc_gamma.log','w')\n",
    "f_p=open(f'{dir_dl2_p}/DL2_mc_proton.log','w')\n",
    "\n",
    "dir_dl1_gamma_test=('/fefs/aswg/workspace/analysis-school-2024/MCP/dl2/input/mc/dl1_stereo_gamma*.h5')   #test gammas\n",
    "dir_dl1_proton_test=('/fefs/aswg/workspace/analysis-school-2024/MCP/dl2/input/mc/dl1_stereo_proton*.h5') #test protons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b024759a",
   "metadata": {},
   "source": [
    "Lines to get files from gamma/proton folders "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2f1c7bbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file_gamma = glob.glob(dir_dl1_gamma_test)\n",
    "input_file_gamma.sort()\n",
    "input_file_proton = glob.glob(dir_dl1_proton_test)\n",
    "input_file_proton.sort()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d02d6f0",
   "metadata": {},
   "source": [
    "Here we use python `subprocess.run()` to run the script and get a log file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b044bc09",
   "metadata": {},
   "outputs": [],
   "source": [
    "for input_file in input_file_gamma: \n",
    "    b=subprocess.run(['lst1_magic_dl1_stereo_to_dl2', f'-d{input_file}', f'-r{dir_rf}',\\\n",
    "        f'-o{dir_dl2_g}'], stdout=f_g, stderr=f_g)     \n",
    "        \n",
    "for input_file in input_file_proton:\n",
    "    c=subprocess.run(['lst1_magic_dl1_stereo_to_dl2', f'-d{input_file}', f'-r{dir_rf}',\\\n",
    "        f'-o{dir_dl2_p}'], stdout=f_p, stderr=f_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3021ed18",
   "metadata": {},
   "source": [
    "To check the log file we do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d541eadb",
   "metadata": {},
   "outputs": [],
   "source": [
    "more $dir_dl2_g/DL2_mc_gamma.log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e936eea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "more $dir_dl2_p/DL2_mc_proton.log"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "070cbc8f",
   "metadata": {},
   "source": [
    "This process just created DL2 files for the MC gammas and protons and put them in the DL2 directory."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b25146b",
   "metadata": {},
   "source": [
    "### Real data\n",
    "\n",
    "Let's repeat the process for the real data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c20272a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "dir_dl2_real=('...../DL2/real')\n",
    "#dir_dl2_real=('/fefs/aswg/workspace/elisa.visentin/sw_school_2024/real')\n",
    "os.makedirs(dir_dl2_real)\n",
    "\n",
    "f=open(f'{dir_dl2_real}/DL2_data.log','w')\n",
    "\n",
    "input_data=('/fefs/aswg/workspace/analysis-school-2024/MCP/dl2/input/real/*.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d45fac1",
   "metadata": {},
   "source": [
    "Lines to get files from data folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "69896208",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file_data = glob.glob(input_data)\n",
    "input_file_data.sort()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d087ce",
   "metadata": {},
   "source": [
    "Here we use python `subprocess.run()` to run the script and also get a log file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "acda4137",
   "metadata": {},
   "outputs": [],
   "source": [
    "for input_file in input_file_data: \n",
    "    d=subprocess.run(['lst1_magic_dl1_stereo_to_dl2', f'-d{input_file}', f'-r{dir_rf}',\\\n",
    "        f'-o{dir_dl2_real}'], stdout=f, stderr=f)     "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "952811a0",
   "metadata": {},
   "source": [
    "To check the log file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9cb9d635",
   "metadata": {},
   "outputs": [],
   "source": [
    "more $dir_dl2_real/DL2_data.log"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48de1b18",
   "metadata": {},
   "source": [
    "This process just created DL2 files for the real data and put them in the DL2 directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e98980e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}