Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update pulsar recipe #22

Merged
merged 9 commits into from
Aug 30, 2024
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions recipes/pulsar_phase/env.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
6 changes: 3 additions & 3 deletions recipes/pulsar_phase/gammapy-pint-environment.yml
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@ channels:
- conda-forge

dependencies:
- gammapy=1.0.1
- python=3.9
- gammapy=1.2
- python=3.11
- pip
- pip:
- pint-pulsar~=0.9.5
- pint-pulsar~=1.0
128 changes: 92 additions & 36 deletions recipes/pulsar_phase/pulsar_phase_computation.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand Down Expand Up @@ -54,7 +54,7 @@
"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",
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
"\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",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this gammapy-pint-environment.yml ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't really know. I think that is was the original env file that was presented to the user. But then we also need this env.yml for the CI. In the end they are the same so I could make gammapy-pint-environment.yml a simlink to env.

"\n",
"`$ conda activate gammapy-pint`\n",
"\n",
Expand Down Expand Up @@ -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__)"
]
},
{
Expand Down Expand Up @@ -143,7 +146,8 @@
"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\"\n",
"DL3_dir = \"/Users/mregeard/Workspace/dev/code/gammapy/gammapy-data/magic/rad_max/data/\""
]
},
{
Expand All @@ -154,7 +158,7 @@
"outputs": [],
"source": [
"# Read DataStore from a directory\n",
"data_store = DataStore.from_dir(DL3_direc)"
"data_store = DataStore.from_dir(DL3_dir)"
]
},
{
Expand All @@ -173,7 +177,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",
")"
]
},
Expand Down Expand Up @@ -298,7 +302,7 @@
"\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, 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 (see next section)."
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
]
},
{
Expand All @@ -309,15 +313,17 @@
"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\""
]
},
{
"cell_type": "markdown",
"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 file 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 feature from Tempo2. \n",
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
"\n",
"In our case, an error occur when using the ephemeris file with PINT. This is due to the `JUMP` line. In order to proceed, one as to either comment (with #) or remove this line. Note that this line is not important for gamma-ray instrument, so it is totally fine to ignore it."
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
]
},
{
Expand Down Expand Up @@ -348,9 +354,7 @@
"cell_type": "code",
"execution_count": null,
"id": "d24e1c92",
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"model = pmodels.get_model(ephemeris_file)\n",
Expand All @@ -366,7 +370,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. There are other additional parameters that can be checked in the [PINT documentation](https://nanograv-pint.readthedocs.io). To get the complete set of parameter of the ephemeris file, one can simply print the model:\n",
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
"`print(model)`"
]
},
{
"cell_type": "markdown",
"id": "f02f0ca5-09f2-44d4-ad76-518f8964b6a1",
"metadata": {},
"source": [
"As said above, we should make sure that 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:"
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
]
},
{
"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": [
"In the case you have several observations, if you know that they 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 build a small function such as the following one:"
MRegeard marked this conversation as resolved.
Show resolved Hide resolved
]
},
{
"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)"
]
},
{
Expand All @@ -381,9 +447,7 @@
"cell_type": "code",
"execution_count": null,
"id": "a145cdde",
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
Expand Down Expand Up @@ -462,13 +526,11 @@
"cell_type": "code",
"execution_count": null,
"id": "f3269a2c",
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"# Show original table\n",
"print(table)"
"table"
]
},
{
Expand All @@ -494,9 +556,7 @@
"cell_type": "code",
"execution_count": null,
"id": "7f29187e",
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"# Show table with phases\n",
Expand Down Expand Up @@ -525,9 +585,7 @@
"cell_type": "code",
"execution_count": null,
"id": "af220557",
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"table.meta"
Expand Down Expand Up @@ -697,15 +755,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",
")"
]
},
Expand Down Expand Up @@ -737,7 +793,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"
]
Expand Down Expand Up @@ -812,7 +868,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",
")"
]
},
Expand Down Expand Up @@ -889,7 +945,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.16"
"version": "3.11.9"
}
},
"nbformat": 4,
Expand Down
Loading