diff --git a/recipes/pulsar_phase/env.yml b/recipes/pulsar_phase/env.yml index 4d44158..76a388c 100644 --- a/recipes/pulsar_phase/env.yml +++ b/recipes/pulsar_phase/env.yml @@ -5,10 +5,10 @@ channels: - conda-forge dependencies: - - gammapy=1.0.2 - - python=3.9 + - gammapy=1.2 + - python=3.11 - scipy<1.12 - jupyter - pip - pip: - - pint-pulsar~=0.9.3 + - pint-pulsar~=1.0 diff --git a/recipes/pulsar_phase/gammapy-pint-environment.yml b/recipes/pulsar_phase/gammapy-pint-environment.yml deleted file mode 100644 index d109cba..0000000 --- a/recipes/pulsar_phase/gammapy-pint-environment.yml +++ /dev/null @@ -1,12 +0,0 @@ -#Declare your specific environment -name: pulsar-gammapy - -channels: - - conda-forge - -dependencies: - - gammapy=1.0.1 - - python=3.9 - - pip - - pip: - - pint-pulsar~=0.9.5 diff --git a/recipes/pulsar_phase/gammapy-pint-environment.yml b/recipes/pulsar_phase/gammapy-pint-environment.yml new file mode 120000 index 0000000..066b96f --- /dev/null +++ b/recipes/pulsar_phase/gammapy-pint-environment.yml @@ -0,0 +1 @@ +env.yml \ No newline at end of file diff --git a/recipes/pulsar_phase/pulsar_phase_computation.ipynb b/recipes/pulsar_phase/pulsar_phase_computation.ipynb index 5a19455..1e0d2f2 100644 --- a/recipes/pulsar_phase/pulsar_phase_computation.ipynb +++ b/recipes/pulsar_phase/pulsar_phase_computation.ipynb @@ -15,9 +15,9 @@ "source": [ "This notebook has been done for the following version of Gammapy and PINT:\n", "\n", - "Gammapy version : 1.0.1\n", + "Gammapy version : 1.2\n", "\n", - "PINT version : 0.9.5" + "PINT version : 1.0" ] }, { @@ -30,7 +30,7 @@ "1. The time of arrivals (TOAs). These times should have very high precision due to the common fast periods of pulsars. Usually these times are already stored in the EventList. For the computation of pulsar timing, times must be corrected in order to be referenced in the Solar System barycenter (SSB) because this system can nearly be regarded as an inertial reference frame with respect to the pulsar.\n", "\n", "\n", - "2. The model of rotation of the pulsar, also known as ephemeris, at the epoch of the observations. These ephemerides are stored in an specific format and saved as .par files and contain informations on the periods, derivatives of the periods, coordinates, glitches, etc.\n", + "2. The model of rotation of the pulsar, also known as ephemeris, at the epoch of the observations. These ephemerides are stored in a specific format and saved as .par files which contain the periods, derivatives of the periods, coordinates, glitches, etc.\n", "\n", "__For the following steps of this tutorial, we need the original EventLists from the DL3 files, and a model in .par format.__\n", "\n", @@ -51,10 +51,10 @@ "id": "a9c72d26", "metadata": {}, "source": [ - "In order to run this notebook, one needs to have installed Gammapy as well as PINT (see documentation above) in the same environment. We recommend to first install Gammapy and then install PINT using your prefered package manager.\n", + "To run this notebook, you must have Gammapy and PINT (see documentation above) installed in the same environment. We recommend installing Gammapy first and then installing PINT using your preferred package manager.\n", "\n", "\n", - "`$ conda env create -n gammapy-pint -f gammapy-1.0-environment.yml`\n", + "`$ conda env create -n gammapy-pint -f gammapy-pint-environment.yml`\n", "\n", "`$ conda activate gammapy-pint`\n", "\n", @@ -97,7 +97,10 @@ "from astropy.coordinates import SkyCoord\n", "import numpy as np\n", "from pathlib import Path\n", - "from gammapy.data import DataStore, EventList, Observation" + "from gammapy.data import DataStore, EventList, Observation\n", + "import logging\n", + "\n", + "log = logging.getLogger(__name__)" ] }, { @@ -105,7 +108,7 @@ "id": "81ae3a58", "metadata": {}, "source": [ - "And we also need some imports from PINT:" + "We also need some imports from PINT:" ] }, { @@ -132,7 +135,7 @@ "id": "b486ec45", "metadata": {}, "source": [ - "First we neeed to define the data sample. In this notebook we will use two runs from the MAGIC gammapy data sample available in https://github.com/gammapy/gammapy-data" + "First we need to define the data sample. In this notebook we will use two runs from the MAGIC gammapy data sample available in https://github.com/gammapy/gammapy-data" ] }, { @@ -143,7 +146,7 @@ "outputs": [], "source": [ "# Define the directory containing the DL3 data\n", - "DL3_direc = \"$GAMMAPY_DATA/magic/rad_max/data\"" + "DL3_dir = \"$GAMMAPY_DATA/magic/rad_max/data\"" ] }, { @@ -154,7 +157,7 @@ "outputs": [], "source": [ "# Read DataStore from a directory\n", - "data_store = DataStore.from_dir(DL3_direc)" + "data_store = DataStore.from_dir(DL3_dir)" ] }, { @@ -173,7 +176,7 @@ "outputs": [], "source": [ "target_pos = SkyCoord(\n", - " ra=083.6331144560900, dec=+22.0144871383400, unit=\"deg\", frame=\"icrs\"\n", + " ra=083.633, dec=+22.014, unit=\"deg\", frame=\"icrs\"\n", ")" ] }, @@ -278,7 +281,7 @@ "id": "e93bfbeb", "metadata": {}, "source": [ - "Now we have the TOAs of the events in the system of the telescope. Please note that the actual precision of the times is higher than the diplayed output (and we really need this precision for the pulsar analysis!). In the next step, the timing in the SSB and the phase for each TOA has to be created. " + "Now we have the TOAs for the events in the system of the telescope. Please note: the actual precision of the times is higher than the displayed output (and we really need this precision for the pulsar analysis!). In the next step, the timing in the SSB and the phase for each TOA has to be created." ] }, { @@ -294,11 +297,11 @@ "id": "fb536b7e", "metadata": {}, "source": [ - "In order to compute the phases of a pulsar, one needs an ephemeris file, usually store as a .par file. \n", + "In order to compute the phases of a pulsar, one needs an ephemeris file, typically stored as a .par file.\n", "\n", "In the following, we will use an ephemeris file for the Crab provided by Fermi-LAT, see [Kerr, M.; Ray, P. S.; et al; 2015](https://arxiv.org/abs/1510.05099). This ephemeris file for the Crab pulsar can be found alongside other pulsar ephemeris files at this [confluence page]( https://confluence.slac.stanford.edu/display/GLAMCOG/LAT+Gamma-ray+Pulsar+Timing+Models). \n", "\n", - "However, be aware that most of these ephemeris files are not up-to-date. Therefore they could give bad results on the phase computation. In particular, one should always checked that the MJD of the observations one wants to phased lies between the `START`and `FINISH`entry of the ephemeris file." + "However, it is important to note that many of the ephemeris files are not up-to-date. Therefore, they could give bad results on the phase computation. In particular, you should always check that the MJD of the observations one wants to phase lies between the `START` and `FINISH` entries of the ephemeris file (see next section)." ] }, { @@ -309,7 +312,7 @@ "outputs": [], "source": [ "# Path to the ephemeris file\n", - "ephemeris_file = \"./0534+2200_ApJ_708_1254_2010.par\"" + "ephemeris_file = \"0534+2200_ApJ_708_1254_2010.par\"" ] }, { @@ -317,7 +320,9 @@ "id": "52eb3086", "metadata": {}, "source": [ - "Note that sometimes one needs to change some of the parameters of the ephemeris file that are not used in gamma-ray astronomy by hand. For instance, here we have removed the 'JUMP' line since it does not have any effect in our computation and raise an error in PINT. The ephemeris file provided with this notebook does not have this line. " + "Note that *Fermi*-LAT ephemeris files are created primarily by and for [Tempo2](https://www.pulsarastronomy.net/pulsar/software/tempo2). Most of the time, using such ephemeris file with PINT will not raise any issues. However, in a few cases, PINT does not support features from Tempo2. \n", + "\n", + "In our case, an error occurs when using the ephemeris file with PINT. This is due to the `JUMP` line. To proceed, simply comment out the line (with #) or remove it. Note that this line is not important for the gamma-ray instruments, so it is acceptable to disregard it." ] }, { @@ -348,9 +353,7 @@ "cell_type": "code", "execution_count": null, "id": "d24e1c92", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "model = pmodels.get_model(ephemeris_file)\n", @@ -366,7 +369,69 @@ "id": "0f8cd0d8", "metadata": {}, "source": [ - "There are multiple parameters such as the name of the source, the interval of validity of the model (START to FINISH), the frequencies of rotation and its derivatives (F0,F1,F2). There are other additional parameters that can be checked in the [PINT documentation](https://nanograv-pint.readthedocs.io)" + "There are multiple parameters such as the name of the source, the frequencies of rotation and its derivatives (F0,F1,F2), the dispersion measure, etc. Check the [PINT documentation](https://nanograv-pint.readthedocs.io) for a list of additional parameters. To obtain the complete set of parameters from the ephemeris file, one can simply print the model:\n", + "`print(model)`" + ] + }, + { + "cell_type": "markdown", + "id": "f02f0ca5-09f2-44d4-ad76-518f8964b6a1", + "metadata": {}, + "source": [ + "As mentioned previously, we should ensure the time of the observation lies within the ephemeris time definition. In our example, we only have one run, so we can check that manually:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2584b571-93a9-4fb6-b8e2-60b8810190f6", + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Ephemeris time definition:\\n{model.START.value} - {model.FINISH.value}\")\n", + "print(f\"Observation time definition:\\n{observation.tstart} - {observation.tstop}\")" + ] + }, + { + "cell_type": "markdown", + "id": "86041174-5ec4-43ba-9035-cf567b47a3bb", + "metadata": {}, + "source": [ + "If you have several observations that are sorted by time, you can manually check for the start time of the first observation and the stop time of the last one. Otherwise, you can create a small function like the following one:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8838ee6c-5b53-4d34-8fd7-94d32daa6bb3", + "metadata": {}, + "outputs": [], + "source": [ + "def check_time(observation, timing_model):\n", + " \"\"\"\n", + " Check that the observation time lies within the time definition of the pulsar\n", + " timing model.\n", + "\n", + " Parameters\n", + " ----------\n", + " observation: `gammapy.data.Observation`\n", + " Observation to check.\n", + " timing_model: `pint.models.TimingModel`\n", + " The timing model that will be used.\n", + " \"\"\"\n", + " model_time = Time([model.START.value, model.FINISH.value], scale=\"tt\", format='mjd')\n", + " if (model_time[0].value > observation.tstart.tt.mjd) or (model_time[1].value < observation.tstop.tt.mjd):\n", + " log.warning(f\"Warning: Observation time of observation {observation.obs_id} goes out of timing model validity time.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d37e6cc-2cae-4192-87d9-88d52216e83f", + "metadata": {}, + "outputs": [], + "source": [ + "check_time(observation, model)" ] }, { @@ -381,14 +446,12 @@ "cell_type": "code", "execution_count": null, "id": "a145cdde", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", - "# Put these to True is your observatory has clock correction files.\n", + "# Set these to True is your observatory has clock correction files.\n", "# If it is set to True but your observatory does not have clock correction files, it will be ignored.\n", "include_bipm = False\n", "include_gps = False\n", @@ -414,7 +477,7 @@ "id": "db97ea5b", "metadata": {}, "source": [ - "Once we have the TOAs object and the model, the phases are easily computed using the model.phase() method. Note that the phases are computed in the interval [-0.5,0.5]. Most of the times, we use the phases in the interval [0,1] so we have to shift the negative ones." + "Once we have the TOAs object and the model, the phases are easily computed using the model.phase() method. Note that the phases are computed in the interval [-0.5,0.5]. Most of the time, we use the phases in the interval [0,1] so we have to shift the negative ones." ] }, { @@ -462,13 +525,11 @@ "cell_type": "code", "execution_count": null, "id": "f3269a2c", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "# Show original table\n", - "print(table)" + "table" ] }, { @@ -494,9 +555,7 @@ "cell_type": "code", "execution_count": null, "id": "7f29187e", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "# Show table with phases\n", @@ -516,18 +575,16 @@ "id": "ea9ba223", "metadata": {}, "source": [ - "At this point, we also want to add meta data to the table. It is very useful to keep track of what has been done to the file. For instance, if we have multiple pulsars in the same file, we want to be able to know quickly which column correspond to which pulsar. Moreover, experience shows that one often use different ephemeris file for the same pulsar. Therefore, it is very useful to have several phase columns in the same file and to be able to know which column correspond to which ephemeris file, parameters, etc.\n", + "At this point, we also want to add metadata to the table. It is very useful to keep track of what has been done to the file. For instance, if a file contains multiple pulsars, we want identify quickly which column corresponds to each pulsar. Moreover, experience has shown that it is common to have different ephemeris files for the same pulsar. Therefore, it is useful to have several phase columns in the same file to easily identify which column corresponds to each ephemeris file, parameters, etc.\n", "\n", - "Since there is not yet a \"standard\" format for such metadata, we propose a template for the essential informations that one wants to save in the header of the event file. First, we look at the present meta info on the table." + "Since there is currently no \"standard\" format for such metadata, we propose a template for the essential information that one wants to save in the header of the event file. First, we look at the present meta info on the table." ] }, { "cell_type": "code", "execution_count": null, "id": "af220557", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "table.meta" @@ -651,9 +708,9 @@ "id": "e369d2dc", "metadata": {}, "source": [ - "In the following, we show how to write the files in a directory contained in the original datastore directory. This follows the logic of DL3 data store and facilitate the manipulation of the HDU table.\n", + "In the following, we show how to write the files in a directory contained in the original datastore directory. This follows the logic of DL3 data store and facilitates the manipulation of the HDU table.\n", "\n", - "If one does not want to save the events files and directly perform the pulsar analysis, this step is not required as well as the step of the meta data handling. However, be aware that for large dataset, the computation of phases can take tens of minutes. " + "If you do not want to save the events files bur rather directly perform the pulsar analysis, you can skip both this step and the step of the handling metadata. However, be aware that for large datasets, the computation of the phases can take tens of minutes." ] }, { @@ -697,15 +754,13 @@ "cell_type": "code", "execution_count": null, "id": "f339c797", - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "# Save the observation object in the specified file_path\n", - "print(\"Writing outputfile in \" + str(file_path))\n", - "observation.events.write(\n", - " filename=file_path, gti=observation.gti, overwrite=True\n", + "print(\"Writing output file in \" + str(file_path))\n", + "new_obs.write(\n", + " path=file_path, include_irfs=False, overwrite=True\n", ")" ] }, @@ -737,7 +792,7 @@ "outputs": [], "source": [ "for entry in new_hdu:\n", - " if entry[\"HDU_NAME\"] == \"EVENTS\" and entry[\"OBS_ID\"] == observation.obs_id:\n", + " if (entry[\"HDU_NAME\"] == \"EVENTS\") and (entry[\"OBS_ID\"] == observation.obs_id):\n", " entry[\"FILE_DIR\"] = \"./\" + str(output_directory)\n", " entry[\"FILE_NAME\"] = filename" ] @@ -765,7 +820,7 @@ "id": "e08bf493", "metadata": {}, "source": [ - "Finally, we need to save the new HDU table in the origianl DL3 directory. Here one should be very careful to not name the new HDU file with the same name as the original HDU file of the data store. Otherwise, the original HDU file will be overwrited. " + "Finally, we need to save the new HDU table in the original DL3 directory. One must be very careful with naming the new HDU file, such that it does not have the same name as the original HDU file of the data store. Otherwise, the original HDU file will be overwritten." ] }, { @@ -785,7 +840,7 @@ "id": "e92cf68c", "metadata": {}, "source": [ - "**Note: Here we use only one approach that could be useful, showing the steps to save the new Event files in a random directory and generate a new modified HDU index table. However, the user is free to chose the absolute path of the EventList and DataStore. For instance, another approach could be making a full copy of the DataStore, or changing the location of the pulsar event files to one that could be more convinient for the user.**" + "**Note: Here we demonstrate only one approach that could be useful, showing the steps to save the new Event files in a directory and generate a new modified HDU index table. However, the user is free to choose the absolute path of the EventList and DataStore. Another approach, for instance, could be making a full copy of the DataStore, or changing the location of the pulsar event files to one that is more convenient for the user.**" ] }, { @@ -801,7 +856,7 @@ "id": "e1036448", "metadata": {}, "source": [ - "Once all of this is done, we just have to open the data store using `DataStore.from_dir()`and passing the pulsar HDU table to it :" + "Once all of this is done, we just have to open the data store using DataStore.from_dir() and pass the pulsar HDU table to it :" ] }, { @@ -812,7 +867,7 @@ "outputs": [], "source": [ "pulsar_datastore = DataStore.from_dir(\n", - " DL3_direc, hdu_table_filename=\"hdu-index-pulsar.fits.gz\"\n", + " DL3_dir, hdu_table_filename=\"hdu-index-pulsar.fits.gz\"\n", ")" ] }, @@ -860,7 +915,7 @@ "id": "1fe89d60", "metadata": {}, "source": [ - "Once we have the corret DataStore and the modified EventList with the phase information, we can do the pulsar analysis using different tools for Gammapy to compute the phaseogram, maps, SED, lightcurve, etc... To do so, one can check the following [Gammapy tutorial](https://docs.gammapy.org/1.0/tutorials/analysis-time/pulsar_analysis.html#sphx-glr-tutorials-analysis-time-pulsar-analysis-py)." + "Once we have the correct DataStore and the modified EventList with the phase information, we can perform the pulsar analysis using different tools available in Gammapy. Allowing us to compute the phaseogram, maps, SED, lightcurve and more. To do so, please refer to the following [Gammapy tutorial](https://docs.gammapy.org/1.0/tutorials/analysis-time/pulsar_analysis.html#sphx-glr-tutorials-analysis-time-pulsar-analysis-py)." ] }, { @@ -889,7 +944,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.16" + "version": "3.11.9" } }, "nbformat": 4,