This is a fixed-text formatted version of a Jupyter notebook

Phase computation for pulsar using PINT

This notebook has been done for the following version of Gammapy and PINT:

Gammapy version : 1.2

PINT version : 1.0

This notebook shows how to compute and add the phase information into the events files of pulsar observations. This step is needed to perform the pulsar analysis with Gammapy and should be the first step in the high level analysis. For the pulsar analysis we need two ingredients:

  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.

  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.

For the following steps of this tutorial, we need the original EventLists from the DL3 files, and a model in .par format.

The main software that we will use to make the barycentric corrections and the phase-folding to the model is the PINT python library, Luo J., Ransom S. et al., 2021, ASCL. For more information about this package, see PINT documentation.

0. Dependencies and imports

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.

$ conda env create -n gammapy-pint -f gammapy-pint-environment.yml

$ conda activate gammapy-pint

$ pip install pint-pulsar

Alternatively, one can also run the yaml environement file provided in the folder of this notebook:

$ conda env create -n gammapy-pint -f gammapy-pint-environment.yml

[1]:

import gammapy
import pint

print(f"Gammapy version : {gammapy.__version__}")
print(f"PINT version : {pint.__version__}")

Gammapy version : 1.2
PINT version : 1.0.1

[2]:

import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord
import numpy as np
from pathlib import Path
from gammapy.data import DataStore, EventList, Observation
import logging

log = logging.getLogger(__name__)

We also need some imports from PINT:

[3]:

import pint.models as pmodels
from pint import toa

1. Reading DataStore

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

[4]:

# Define the directory containing the DL3 data
DL3_dir = "$GAMMAPY_DATA/magic/rad_max/data"

[5]:

# Read DataStore from a directory
data_store = DataStore.from_dir(DL3_dir)

Let’s run this tutorial for the Crab pulsar :

[6]:

target_pos = SkyCoord(ra=083.633, dec=+22.014, unit="deg", frame="icrs")

[7]:

selection = dict(
    type="sky_circle",
    frame="icrs",
    lon=target_pos.ra,
    lat=target_pos.dec,
    radius="5 deg",
)
selected_obs_table = data_store.obs_table.select_observations(selection)

[8]:

obs_id = selected_obs_table["OBS_ID"]
print(obs_id)

 OBS_ID
-------
5029748
5029747

For the following we will take the run 5029747.

[9]:

observations = data_store.get_observations(
    [5029747], required_irf="point-like"
)

[10]:

print(observations)

Observations
Number of observations: 1
Observation

        obs id            : 5029747
        tstart            : 56569.18
        tstop             : 56569.19
        duration          : 1188.11 s
        pointing (icrs)   : 84.0 deg, 22.2 deg

        deadtime fraction : 0.8%

2. Phase-folding with PINT for one observation

Let’s extract the times from the observation:

[11]:

# Extract times from EventList
observation = observations[0]
times = observation.events.time

[12]:

print(times)

[56569.18112859 56569.18112878 56569.18112978 ... 56569.19487702
 56569.19487703 56569.19487876]

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.

2.1 An ephemeris file from Fermi-LAT data.

In order to compute the phases of a pulsar, one needs an ephemeris file, typically stored as a .par file.

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. This ephemeris file for the Crab pulsar can be found alongside other pulsar ephemeris files at this confluence page.

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).

[13]:

# Path to the ephemeris file
ephemeris_file = "0534+2200_ApJ_708_1254_2010.par"

Note that Fermi-LAT ephemeris files are created primarily by and for 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.

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.

2.2 Computing pulsar phases

Now that we have the model and the times of arrival for the different events, we can compute the timing corrections and the pulsar phases needed for the pulsar analysis. In this case, we use the PINT package described in the introduction.

First we will explore our model. We print some of the relevant quantities

[14]:

model = pmodels.get_model(ephemeris_file)
print(model.components["AstrometryEquatorial"])
print(model.components["SolarSystemShapiro"])
print(model.components["DispersionDM"])
print(model.components["AbsPhase"])
print(model.components["Spindown"])

2024-08-30 13:43:54.870 | INFO     | pint.models.absolute_phase:validate:77 - TZRFRQ was 0.0 or None. Setting to infinite frequency.
/home/runner/miniconda3/envs/gammapy-recipes/lib/python3.11/site-packages/pint/models/model_builder.py:230: UserWarning: Unrecognized parfile line 'EPHVER 5'
  warnings.warn(f"Unrecognized parfile line '{p_line}'", UserWarning)
2024-08-30 13:43:54.872 | WARNING  | pint.models.model_builder:__call__:234 - UNITS is not specified. Assuming TDB...

AstrometryEquatorial(
    MJDParameter(   POSEPOCH            50739.0000000000000000 (d) frozen=True),
    floatParameter( PX                  0.0               (mas) frozen=True),
    AngleParameter( RAJ                 5:34:31.94000000  (hourangle) frozen=True),
    AngleParameter( DECJ                22:00:52.10000000 (deg) frozen=True),
    floatParameter( PMRA                -11.8             (mas / yr) +/- 1.0 mas / yr frozen=True),
    floatParameter( PMDEC               4.4               (mas / yr) +/- 1.0 mas / yr frozen=True))
SolarSystemShapiro(
    boolParameter(  PLANET_SHAPIRO      N                 frozen=True))
DispersionDM(
    floatParameter( DM                  56.785579397589822356 (pc / cm3) frozen=True),
    floatParameter( DM1                 0.031279168349770640344 (pc / (yr cm3)) frozen=True),
    MJDParameter(   DMEPOCH             55107.8071585532816240 (d) frozen=True))
AbsPhase(
    MJDParameter(   TZRMJD              55638.1552775999516551 (d) frozen=True),
    strParameter(   TZRSITE             coe               frozen=True),
    floatParameter( TZRFRQ              inf               (MHz) frozen=True))
Spindown(
    floatParameter( F0                  29.716913767510206412 (Hz) +/- 7.440760185772888e-08 Hz frozen=True),
    MJDParameter(   PEPOCH              55555.0000000000000000 (d) frozen=True),
    floatParameter( F1                  -3.7105744257791078768e-10 (Hz / s) +/- 7.182646367277094e-16 Hz / s frozen=True),
    floatParameter( F2                  1.1633675551229037574e-20 (Hz / s2) +/- 2.487306981530037e-23 Hz / s2 frozen=True),
    floatParameter( F3                  -8.6724769747827286635e-31 (Hz / s3) +/- 0.0 Hz / s3 frozen=True))

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 for a list of additional parameters. To obtain the complete set of parameters from the ephemeris file, one can simply print the model: print(model)

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:

[15]:

print(
    f"Ephemeris time definition:\n{model.START.value} - {model.FINISH.value}"
)
print(
    f"Observation time definition:\n{observation.tstart} - {observation.tstop}"
)

Ephemeris time definition:
54686.1526259 - 56583.1591704
Observation time definition:
56569.18112772242 - 56569.19487901596

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:

[16]:

def check_time(observation, timing_model):
    """
    Check that the observation time lies within the time definition of the pulsar
    timing model.

    Parameters
    ----------
    observation: `gammapy.data.Observation`
        Observation to check.
    timing_model: `pint.models.TimingModel`
        The timing model that will be used.
    """
    model_time = Time(
        [model.START.value, model.FINISH.value], scale="tt", format="mjd"
    )
    if (model_time[0].value > observation.tstart.tt.mjd) or (
        model_time[1].value < observation.tstop.tt.mjd
    ):
        log.warning(
            f"Warning: Observation time of observation {observation.obs_id} goes out of timing model validity time."
        )

[17]:

check_time(observation, model)

Now we can compute the phases. For that, we define a list of TOA objects that are the main object of PINT.

[18]:

%%time

# Set these to True is your observatory has clock correction files.
# If it is set to True but your observatory does not have clock correction files, it will be ignored.
include_bipm = False
include_gps = False

# Set this to True or False depending on your ephemeris file.
# Here we can see that the 'PLANET_SHAPIRO' entry is 'N' so we set it to True.
planets = False

# Create a TOA object for each time
toas = toa.get_TOAs_array(
    times=times,
    obs="magic",
    errors=1 * u.microsecond,
    ephem="DE421",
    include_gps=include_gps,
    include_bipm=include_bipm,
    planets=planets,
)

2024-08-30 13:43:55.618 | DEBUG    | pint.toa:__init__:1377 - No pulse number flags found in the TOAs
2024-08-30 13:43:55.633 | DEBUG    | pint.toa:apply_clock_corrections:2224 - Applying clock corrections (include_bipm = False)
2024-08-30 13:43:56.294 | INFO     | pint.observatory:gps_correction:248 - Applying GPS to UTC clock correction (~few nanoseconds)
2024-08-30 13:43:56.294 | DEBUG    | pint.observatory:_load_gps_clock:122 - Loading global GPS clock file
2024-08-30 13:43:56.297 | INFO     | pint.observatory.global_clock_corrections:get_file:128 - File index.txt to be downloaded due to download policy if_expired: https://raw.githubusercontent.com/ipta/pulsar-clock-corrections/main/index.txt
2024-08-30 13:43:56.451 | DEBUG    | pint.observatory.clock_file:__init__:812 - Global clock file gps2utc.clk saving kwargs={'bogus_last_correction': False, 'valid_beyond_ends': False}
2024-08-30 13:43:56.452 | INFO     | pint.observatory.global_clock_corrections:get_file:128 - File T2runtime/clock/gps2utc.clk to be downloaded due to download policy if_missing: https://raw.githubusercontent.com/ipta/pulsar-clock-corrections/main/T2runtime/clock/gps2utc.clk
2024-08-30 13:43:56.668 | DEBUG    | pint.observatory.clock_file:read_tempo2_clock_file:463 - Loading TEMPO2-format observatory clock correction file gps2utc.clk (/home/runner/.astropy/cache/download/url/d3c81b5766f4bfb84e65504c8a453085/contents) with bogus_last_correction=False
2024-08-30 13:43:56.687 | INFO     | pint.observatory:find_clock_file:994 - Using global clock file for gps2utc.clk with bogus_last_correction=False
2024-08-30 13:43:56.689 | INFO     | pint.observatory.topo_obs:clock_corrections:354 - Observatory magic requires no clock corrections.
2024-08-30 13:43:59.900 | DEBUG    | pint.toa:compute_TDBs:2270 - Computing TDB columns.
2024-08-30 13:43:59.901 | DEBUG    | pint.toa:compute_TDBs:2291 - Using EPHEM = DE421 for TDB calculation.
2024-08-30 13:44:00.968 | DEBUG    | pint.toa:compute_posvels:2371 - Computing PosVels of observatories and Earth, using DE421
2024-08-30 13:44:02.516 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-08-30 13:44:15.633 | DEBUG    | pint.toa:compute_posvels:2424 - SSB obs pos [1.47007462e+11 2.56889811e+10 1.11245045e+10] m
2024-08-30 13:44:15.688 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-08-30 13:44:15.694 | DEBUG    | pint.toa:compute_posvels:2438 - Adding columns ssb_obs_pos ssb_obs_vel obs_sun_pos

CPU times: user 7.15 s, sys: 199 ms, total: 7.35 s
Wall time: 20.8 s

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.

[19]:

# Compute phases
phases = model.phase(toas, abs_phase=True)[1]

# Shift phases to the interval (0,1]
phases = np.where(phases < 0.0, phases + 1.0, phases)

2024-08-30 13:44:15.729 | DEBUG    | pint.models.glitch:glitch_phase:221 - Glitch phase for glitch 1: 0.0 
2024-08-30 13:44:15.735 | DEBUG    | pint.models.absolute_phase:get_TZR_toa:100 - Creating and dealing with the single TZR_toa for absolute phase
2024-08-30 13:44:15.738 | DEBUG    | pint.toa:__init__:1377 - No pulse number flags found in the TOAs
2024-08-30 13:44:15.738 | DEBUG    | pint.toa:apply_clock_corrections:2224 - Applying clock corrections (include_bipm = False)
2024-08-30 13:44:15.760 | DEBUG    | pint.toa:compute_TDBs:2270 - Computing TDB columns.
2024-08-30 13:44:15.761 | DEBUG    | pint.toa:compute_TDBs:2291 - Using EPHEM = DE421 for TDB calculation.
2024-08-30 13:44:15.764 | DEBUG    | pint.toa:compute_posvels:2371 - Computing PosVels of observatories and Earth, using DE421
2024-08-30 13:44:15.834 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-08-30 13:44:15.836 | DEBUG    | pint.toa:compute_posvels:2424 - SSB obs pos [-1.49278181e+08  7.07659442e+06  3.07113250e+06] km
2024-08-30 13:44:15.907 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-08-30 13:44:15.908 | DEBUG    | pint.toa:compute_posvels:2438 - Adding columns ssb_obs_pos ssb_obs_vel obs_sun_pos
2024-08-30 13:44:15.909 | DEBUG    | pint.models.absolute_phase:get_TZR_toa:121 - Done with TZR_toa
2024-08-30 13:44:15.916 | DEBUG    | pint.models.glitch:glitch_phase:221 - Glitch phase for glitch 1: 0.0

3. Adding phases and metadata to an EventList and put it in a new Observation.

Once the phases are computed we need to create a new EventList table that includes both the original information of the events and the phase information in extra columns. This is necessary for Gammapy to read the phases and use them as an extra variable of each event.

[20]:

# Extract the table of the EventList
table = observation.events.table

[21]:

# Show original table
table

[21]:
Table length=11189
EVENT_IDTIMERADECENERGY
sdegdegTeV
int64float64float32float32float32
2402333778852.509924984.5945722.030880.18194601
2408333778852.526715384.2146223.449140.08397394
2434333778852.6131505483.52470422.7257920.10596932
2445333778852.669014283.7695722.4510060.19733498
2478333778852.762793983.47851623.4845940.08522219
2481333778852.777854983.7151721.9851151.0020943
2513333778852.864446782.42119622.5676520.14374068
2544333778852.982606483.6413622.0413150.10316629
2559333778853.026941484.06917622.973370.047184493
...............
356223333780039.460049284.1161522.5575050.08110082
356227333780039.4710536683.4153421.673440.2096362
356242333780039.517909583.5516522.7729850.17672835
356282333780039.6299751484.4613321.693570.05068718
356473333780040.338647984.4544121.1598280.1831569
356478333780040.354892683.6833623.4449880.06305169
356485333780040.374132284.3351721.283380.060539745
356486333780040.375515984.8588622.1162220.123453744
356526333780040.5247600784.8692921.2909160.13630114
[22]:

# Add a column for the phases to the table
table["PHASE"] = phases.astype("float64")

Note that you can add multiple columns to a same file, only the name of the column has to be unique, eg table['PHASE_SRC1'], table['PHASE_SRC2'] etc”

[23]:

# Show table with phases
table

[23]:
Table length=11189
EVENT_IDTIMERADECENERGYPHASE
sdegdegTeV
int64float64float32float32float32float64
2402333778852.509924984.5945722.030880.181946010.3934918682770723
2408333778852.526715384.2146223.449140.083973940.8919525879696997
2434333778852.6131505483.52470422.7257920.105969320.45797881185987477
2445333778852.669014283.7695722.4510060.197334980.11641851092409164
2478333778852.762793983.47851623.4845940.085222190.900480594501521
2481333778852.777854983.7151721.9851151.00209430.34760107795590983
2513333778852.864446782.42119622.5676520.143740680.9182740088921875
2544333778852.982606483.6413622.0413150.103166290.42611240294578184
2559333778853.026941484.06917622.973370.0471844930.7422974209419148
..................
356223333780039.460049284.1161522.5575050.081100820.6914919573678725
356227333780039.4710536683.4153421.673440.20963620.018183853280195963
356242333780039.517909583.5516522.7729850.176728350.4092061182772071
356282333780039.6299751484.4613321.693570.050687180.7361270045119108
356473333780040.338647984.4544121.1598280.18315690.7746788606122589
356478333780040.354892683.6833623.4449880.063051690.25693976618589837
356485333780040.374132284.3351721.283380.0605397450.8281108809408599
356486333780040.375515984.8588622.1162220.1234537440.8691880225176731
356526333780040.5247600784.8692921.2909160.136301140.29983892821377583

Now we can see that the ‘PHASE’ column has been added to the table

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.

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.

[24]:

table.meta

[24]:

OrderedDict([('EXTNAME', 'EVENTS'),
             ('HDUCLASS', 'GADF'),
             ('HDUDOC',
              'https://github.com/open-gamma-ray-astro/gamma-astro-data-formats'),
             ('HDUVERS', '0.2'),
             ('HDUCLAS1', 'EVENTS'),
             ('OBS_ID', 5029747),
             ('TSTART', 333778852.435217),
             ('TSTOP', 333780040.546979),
             ('ONTIME', 1188.111761868),
             ('LIVETIME', 1178.06621791733),
             ('DEADC', 0.991545),
             ('EQUINOX', 2000.0),
             ('RADECSYS', 'FK5'),
             ('ORIGIN', 'MAGIC Collaboration'),
             ('TELESCOP', 'MAGIC'),
             ('INSTRUME', 'MAGIC stereo'),
             ('CREATOR', 'MAGIC_DL3'),
             ('VERSION', '0.1.8'),
             ('MJDREFI', 52706),
             ('MJDREFF', 0.0),
             ('TIMEUNIT', 's'),
             ('TIMESYS', 'UTC'),
             ('TIMEREF', 'local'),
             ('OBJECT', 'CrabNebula'),
             ('RA_OBJ', 83.63333),
             ('DEC_OBJ', 22.01444),
             ('OBS_MODE', 'POINTING'),
             ('RA_PNT', 83.98333),
             ('DEC_PNT', 22.24389),
             ('GEOLON', -17.89),
             ('GEOLAT', 28.761944),
             ('ALTITUDE', 2231.28),
             ('ALT_PNT', 69.69974),
             ('AZ_PNT', 103.8848)])

[25]:

def get_log(ephemeris_file, phase_column_name="PHASE"):
    return (
        "COLUMN_PHASE: "
        + str(phase_column_name)
        + "; PINT_VERS: "
        + pint.__version__
        + "; GAMMAPY_VERS: "
        + gammapy.__version__
        + "; EPHEM_FILE: "
        + ephemeris_file
        + "; PSRJ :"
        + str(model.PSR.value)
        + "; START: "
        + str(model.START.value)
        + "; FINISH: "
        + str(model.FINISH.value)
        + "; TZRMJD: "
        + str(model.TZRMJD.value)
        + "; TZRSITE: "
        + str(model.TZRSITE.value)
        + "; TZRFREQ: "
        + str(model.TZRFRQ.value)
        + "; EPHEM: "
        + str(model.EPHEM.value)
        + "; EPHEM_RA: "
        + str(model.RAJ.value)
        + "; EPHEM_DEC: "
        + str(model.DECJ.value)
        + "; PHASE_OFFSET: "
        + "default = 0"
        + "; DATE: "
        + str(Time.now().mjd)
        + ";"
    )

[26]:

phase_log = get_log(ephemeris_file=ephemeris_file, phase_column_name="PHASE")
print(phase_log)

COLUMN_PHASE: PHASE; PINT_VERS: 1.0.1; GAMMAPY_VERS: 1.2; EPHEM_FILE: 0534+2200_ApJ_708_1254_2010.par; PSRJ :J0534+2200; START: 54686.1526259; FINISH: 56583.1591704; TZRMJD: 55638.155277599951656; TZRSITE: coe; TZRFREQ: inf; EPHEM: DE405; EPHEM_RA: 5.575538888888889; EPHEM_DEC: 22.01447222222222; PHASE_OFFSET: default = 0; DATE: 60552.572407022526;

[27]:

# Add the generated string to the meta data of the table
table.meta["PH_LOG"] = phase_log

[28]:

table.meta

[28]:

OrderedDict([('EXTNAME', 'EVENTS'),
             ('HDUCLASS', 'GADF'),
             ('HDUDOC',
              'https://github.com/open-gamma-ray-astro/gamma-astro-data-formats'),
             ('HDUVERS', '0.2'),
             ('HDUCLAS1', 'EVENTS'),
             ('OBS_ID', 5029747),
             ('TSTART', 333778852.435217),
             ('TSTOP', 333780040.546979),
             ('ONTIME', 1188.111761868),
             ('LIVETIME', 1178.06621791733),
             ('DEADC', 0.991545),
             ('EQUINOX', 2000.0),
             ('RADECSYS', 'FK5'),
             ('ORIGIN', 'MAGIC Collaboration'),
             ('TELESCOP', 'MAGIC'),
             ('INSTRUME', 'MAGIC stereo'),
             ('CREATOR', 'MAGIC_DL3'),
             ('VERSION', '0.1.8'),
             ('MJDREFI', 52706),
             ('MJDREFF', 0.0),
             ('TIMEUNIT', 's'),
             ('TIMESYS', 'UTC'),
             ('TIMEREF', 'local'),
             ('OBJECT', 'CrabNebula'),
             ('RA_OBJ', 83.63333),
             ('DEC_OBJ', 22.01444),
             ('OBS_MODE', 'POINTING'),
             ('RA_PNT', 83.98333),
             ('DEC_PNT', 22.24389),
             ('GEOLON', -17.89),
             ('GEOLAT', 28.761944),
             ('ALTITUDE', 2231.28),
             ('ALT_PNT', 69.69974),
             ('AZ_PNT', 103.8848),
             ('PH_LOG',
              'COLUMN_PHASE: PHASE; PINT_VERS: 1.0.1; GAMMAPY_VERS: 1.2; EPHEM_FILE: 0534+2200_ApJ_708_1254_2010.par; PSRJ :J0534+2200; START: 54686.1526259; FINISH: 56583.1591704; TZRMJD: 55638.155277599951656; TZRSITE: coe; TZRFREQ: inf; EPHEM: DE405; EPHEM_RA: 5.575538888888889; EPHEM_DEC: 22.01447222222222; PHASE_OFFSET: default = 0; DATE: 60552.572407022526;')])

Once this is done, we can put back the table in a new EventList object and in a new Observation object.

[29]:

# Create new event list and add it to observation object
new_event_list = EventList(table)
new_obs = observation.copy(in_memory=True, events=new_event_list)

'THETA' axis is stored as a scalar -- converting to 1D array.
'THETA' axis is stored as a scalar -- converting to 1D array.
'THETA' axis is stored as a scalar -- converting to 1D array.

[30]:

new_obs.events.table

[30]:
Table length=11189
EVENT_IDTIMERADECENERGYPHASE
sdegdegTeV
int64float64float32float32float32float64
2402333778852.509924984.5945722.030880.181946010.3934918682770723
2408333778852.526715384.2146223.449140.083973940.8919525879696997
2434333778852.6131505483.52470422.7257920.105969320.45797881185987477
2445333778852.669014283.7695722.4510060.197334980.11641851092409164
2478333778852.762793983.47851623.4845940.085222190.900480594501521
2481333778852.777854983.7151721.9851151.00209430.34760107795590983
2513333778852.864446782.42119622.5676520.143740680.9182740088921875
2544333778852.982606483.6413622.0413150.103166290.42611240294578184
2559333778853.026941484.06917622.973370.0471844930.7422974209419148
..................
356223333780039.460049284.1161522.5575050.081100820.6914919573678725
356227333780039.4710536683.4153421.673440.20963620.018183853280195963
356242333780039.517909583.5516522.7729850.176728350.4092061182772071
356282333780039.6299751484.4613321.693570.050687180.7361270045119108
356473333780040.338647984.4544121.1598280.18315690.7746788606122589
356478333780040.354892683.6833623.4449880.063051690.25693976618589837
356485333780040.374132284.3351721.283380.0605397450.8281108809408599
356486333780040.375515984.8588622.1162220.1234537440.8691880225176731
356526333780040.5247600784.8692921.2909160.136301140.29983892821377583

4. Save new Event List and writing a modify HDU index table

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.

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.

[31]:

data_store.hdu_table.base_dir

[31]:

PosixPath('/home/runner/work/gammapy-recipes/gammapy-recipes/gammapy-datasets/1.2/magic/rad_max/data')

[32]:

# Define output directory and filename
datastore_dir = str(data_store.hdu_table.base_dir) + "/"
output_directory = "pulsar_events_file/"
output_path = datastore_dir + output_directory
filename = f"dl3_pulsar_{observation.obs_id:04d}.fits.gz"
file_path = output_path + filename

Path(output_path).mkdir(parents=True, exist_ok=True)

[33]:

output_path

[33]:

'/home/runner/work/gammapy-recipes/gammapy-recipes/gammapy-datasets/1.2/magic/rad_max/data/pulsar_events_file/'

[34]:

# Save the observation object in the specified file_path
print("Writing output file in " + str(file_path))
new_obs.write(path=file_path, include_irfs=False, overwrite=True)

Writing output file in /home/runner/work/gammapy-recipes/gammapy-recipes/gammapy-datasets/1.2/magic/rad_max/data/pulsar_events_file/dl3_pulsar_5029747.fits.gz

Once the file has been written, we want to write a modified version of the HDU table. This is mandatory if we want to open the phased events file together with its associated IRFs.

[35]:

# Print the current data store HDU table.
new_hdu = data_store.hdu_table.copy()
new_hdu

[35]:
HDUIndexTable length=10
OBS_IDHDU_TYPEHDU_CLASSFILE_DIRFILE_NAMEHDU_NAME
int64bytes30bytes30bytes100bytes50bytes30
5029748eventsevents./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsEVENTS
5029748gtigti./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsGTI
5029748rad_maxrad_max_2d./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsRAD_MAX
5029748aeffaeff_2d./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsEFFECTIVE AREA
5029748edispedisp_2d./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsENERGY DISPERSION
5029747eventsevents./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsEVENTS
5029747gtigti./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsGTI
5029747rad_maxrad_max_2d./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsRAD_MAX
5029747aeffaeff_2d./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsEFFECTIVE AREA
5029747edispedisp_2d./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsENERGY DISPERSION
[36]:

for entry in new_hdu:
    if (entry["HDU_NAME"] == "EVENTS") and (
        entry["OBS_ID"] == observation.obs_id
    ):
        entry["FILE_DIR"] = "./" + str(output_directory)
        entry["FILE_NAME"] = filename

[37]:

new_hdu

[37]:
HDUIndexTable length=10
OBS_IDHDU_TYPEHDU_CLASSFILE_DIRFILE_NAMEHDU_NAME
int64bytes30bytes30bytes100bytes50bytes30
5029748eventsevents./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsEVENTS
5029748gtigti./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsGTI
5029748rad_maxrad_max_2d./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsRAD_MAX
5029748aeffaeff_2d./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsEFFECTIVE AREA
5029748edispedisp_2d./20131004_05029748_DL3_CrabNebula-W0.40+215.fitsENERGY DISPERSION
5029747eventsevents./pulsar_events_file/dl3_pulsar_5029747.fits.gzEVENTS
5029747gtigti./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsGTI
5029747rad_maxrad_max_2d./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsRAD_MAX
5029747aeffaeff_2d./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsEFFECTIVE AREA
5029747edispedisp_2d./20131004_05029747_DL3_CrabNebula-W0.40+035.fitsENERGY DISPERSION

We see that the FILE_DIRand FILE_NAMEentry have been modified for our phased events file.

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.

[38]:

new_hdu.write(
    datastore_dir + "hdu-index-pulsar.fits.gz", format="fits", overwrite=True
)

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.

5. Opening the new DataStore

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 :

[39]:

pulsar_datastore = DataStore.from_dir(
    DL3_dir, hdu_table_filename="hdu-index-pulsar.fits.gz"
)

[40]:

observations = pulsar_datastore.get_observations(
    [5029747], required_irf="point-like"
)
observations[0].available_hdus

[40]:

['events', 'gti', 'aeff', 'edisp', 'rad_max']

[41]:

observations[0].events.table

[41]:
Table length=11189
EVENT_IDTIMERADECENERGYPHASE
sdegdegTeV
int64float64float32float32float32float64
2402333778852.509924984.5945722.030880.181946010.3934918682770723
2408333778852.526715384.2146223.449140.083973940.8919525879696997
2434333778852.6131505483.52470422.7257920.105969320.45797881185987477
2445333778852.669014283.7695722.4510060.197334980.11641851092409164
2478333778852.762793983.47851623.4845940.085222190.900480594501521
2481333778852.777854983.7151721.9851151.00209430.34760107795590983
2513333778852.864446782.42119622.5676520.143740680.9182740088921875
2544333778852.982606483.6413622.0413150.103166290.42611240294578184
2559333778853.026941484.06917622.973370.0471844930.7422974209419148
..................
356223333780039.460049284.1161522.5575050.081100820.6914919573678725
356227333780039.4710536683.4153421.673440.20963620.018183853280195963
356242333780039.517909583.5516522.7729850.176728350.4092061182772071
356282333780039.6299751484.4613321.693570.050687180.7361270045119108
356473333780040.338647984.4544121.1598280.18315690.7746788606122589
356478333780040.354892683.6833623.4449880.063051690.25693976618589837
356485333780040.374132284.3351721.283380.0605397450.8281108809408599
356486333780040.375515984.8588622.1162220.1234537440.8691880225176731
356526333780040.5247600784.8692921.2909160.136301140.29983892821377583

We can see that we recover both the IRFs and the events file with the phase column.

6. Pulsar analysis tools with gammapy

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.

Recipe made by Alvaros Mas, Maxime Regeard, Jan Lukas Schubert.