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

PINT version : 0.9.5

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 an specific format and saved as .par files and contain informations on 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

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.

$ conda env create -n gammapy-pint -f gammapy-1.0-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.0.2
PINT version : 0.9.8

[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

And we also need some imports from PINT:

[3]:

import pint.models as pmodels
from pint import toa

1. Reading DataStore

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

[4]:

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

[5]:

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

Let’s run this tutorial for the Crab pulsar :

[6]:

target_pos = SkyCoord(
    ra=083.6331144560900, dec=+22.0144871383400, 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 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.

2.1 An ephemeris file from Fermi-LAT data.

In order to compute the phases of a pulsar, one needs an ephemeris file, usually store 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, 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 STARTand FINISHentry of the ephemeris file.

[13]:

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

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.

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-03-27 11:26:26.901 | INFO     | pint.models.absolute_phase:validate:74 - TZRFRQ was 0.0 or None. Setting to infinite frequency.
/usr/share/miniconda3/envs/gammapy-recipes/lib/python3.9/site-packages/pint/models/model_builder.py:198: UserWarning: Unrecognized parfile line 'EPHVER 5'
  warnings.warn(f"Unrecognized parfile line '{p_line}'", UserWarning)
2024-03-27 11:26:26.903 | WARNING  | pint.models.model_builder:__call__:202 - 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 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

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

[15]:

%%time

# Put 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-03-27 11:26:27.247 | DEBUG    | pint.toa:__init__:1358 - No pulse number flags found in the TOAs
2024-03-27 11:26:27.263 | DEBUG    | pint.toa:apply_clock_corrections:2200 - Applying clock corrections (include_gps = False, include_bipm = False)
2024-03-27 11:26:27.923 | INFO     | pint.observatory.topo_obs:clock_corrections:365 - Observatory magic requires no clock corrections.
2024-03-27 11:26:31.192 | DEBUG    | pint.toa:compute_TDBs:2251 - Computing TDB columns.
2024-03-27 11:26:31.193 | DEBUG    | pint.toa:compute_TDBs:2272 - Using EPHEM = DE421 for TDB calculation.
2024-03-27 11:26:32.378 | DEBUG    | pint.toa:compute_posvels:2350 - Computing PosVels of observatories and Earth, using DE421
2024-03-27 11:26:33.420 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-03-27 11:26:44.359 | DEBUG    | pint.toa:compute_posvels:2403 - SSB obs pos [1.47007462e+11 2.56889811e+10 1.11245045e+10] m
2024-03-27 11:26:44.587 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-03-27 11:26:44.593 | DEBUG    | pint.toa:compute_posvels:2417 - Adding columns ssb_obs_pos ssb_obs_vel obs_sun_pos

CPU times: user 7.02 s, sys: 188 ms, total: 7.21 s
Wall time: 17.7 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 times, we use the phases in the interval [0,1] so we have to shift the negative ones.

[16]:

# 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-03-27 11:26:44.636 | DEBUG    | pint.models.glitch:glitch_phase:213 - Glitch phase for glitch 1: 0.0 
2024-03-27 11:26:44.640 | DEBUG    | pint.models.absolute_phase:get_TZR_toa:98 - Creating and dealing with the single TZR_toa for absolute phase
2024-03-27 11:26:44.646 | DEBUG    | pint.toa:__init__:1358 - No pulse number flags found in the TOAs
2024-03-27 11:26:44.646 | DEBUG    | pint.toa:apply_clock_corrections:2200 - Applying clock corrections (include_gps = False, include_bipm = False)
2024-03-27 11:26:44.648 | DEBUG    | pint.toa:compute_TDBs:2251 - Computing TDB columns.
2024-03-27 11:26:44.648 | DEBUG    | pint.toa:compute_TDBs:2272 - Using EPHEM = DE421 for TDB calculation.
2024-03-27 11:26:44.651 | DEBUG    | pint.toa:compute_posvels:2350 - Computing PosVels of observatories and Earth, using DE421
2024-03-27 11:26:44.929 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-03-27 11:26:44.931 | DEBUG    | pint.toa:compute_posvels:2403 - SSB obs pos [-1.49278181e+08  7.07659442e+06  3.07113250e+06] km
2024-03-27 11:26:45.150 | INFO     | pint.solar_system_ephemerides:_load_kernel_link:55 - Set solar system ephemeris to de421 from download
2024-03-27 11:26:45.176 | DEBUG    | pint.toa:compute_posvels:2417 - Adding columns ssb_obs_pos ssb_obs_vel obs_sun_pos
2024-03-27 11:26:45.177 | DEBUG    | pint.models.absolute_phase:get_TZR_toa:121 - Done with TZR_toa
2024-03-27 11:26:45.188 | DEBUG    | pint.models.glitch:glitch_phase:213 - 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.

[17]:

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

[18]:

# Show original table
print(table)

EVENT_ID        TIME            RA       DEC       ENERGY
                 s             deg       deg        TeV
-------- ------------------ --------- --------- -----------
    2402  333778852.5099249  84.59457  22.03088  0.18194601
    2408  333778852.5267153  84.21462  23.44914  0.08397394
    2434 333778852.61315054 83.524704 22.725792  0.10596932
    2445  333778852.6690142  83.76957 22.451006  0.19733498
    2478  333778852.7627939 83.478516 23.484594  0.08522219
    2481  333778852.7778549  83.71517 21.985115   1.0020943
    2513  333778852.8644467 82.421196 22.567652  0.14374068
    2544  333778852.9826064  83.64136 22.041315  0.10316629
    2559  333778853.0269414 84.069176  22.97337 0.047184493
    2561  333778853.0339344  84.84237 22.175398 0.118843034
     ...                ...       ...       ...         ...
  356222  333780039.4520397  84.74482 20.894981 0.043312162
  356223  333780039.4600492  84.11615 22.557505  0.08110082
  356227 333780039.47105366  83.41534  21.67344   0.2096362
  356242  333780039.5179095  83.55165 22.772985  0.17672835
  356282 333780039.62997514  84.46133  21.69357  0.05068718
  356473  333780040.3386479  84.45441 21.159828   0.1831569
  356478  333780040.3548926  83.68336 23.444988  0.06305169
  356485  333780040.3741322  84.33517  21.28338 0.060539745
  356486  333780040.3755159  84.85886 22.116222 0.123453744
  356526 333780040.52476007  84.86929 21.290916  0.13630114
Length = 11189 rows

[19]:

# 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”

[20]:

# Show table with phases
table

[20]:
Table length=11189
EVENT_IDTIMERADECENERGYPHASE
sdegdegTeV
int64float64float32float32float32float64
2402333778852.509924984.5945722.030880.181946010.39349195931383935
2408333778852.526715384.2146223.449140.083973940.8919526790064667
2434333778852.6131505483.52470422.7257920.105969320.4579789031294724
2445333778852.669014283.7695722.4510060.197334980.11641860219368935
2478333778852.762793983.47851623.4845940.085222190.9004806857711186
2481333778852.777854983.7151721.9851151.00209430.3476011689926769
2513333778852.864446782.42119622.5676520.143740680.9182740999289545
2544333778852.982606483.6413622.0413150.103166290.42611249421537956
2559333778853.026941484.06917622.973370.0471844930.7422975122115125
..................
356222333780039.452039784.7448220.8949810.0433121620.4537121484157113
356223333780039.460049284.1161522.5575050.081100820.6914920486374702
356227333780039.4710536683.4153421.673440.20963620.018183944316962984
356242333780039.517909583.5516522.7729850.176728350.40920620046641115
356282333780039.6299751484.4613321.693570.050687180.7361270955486778
356473333780040.338647984.4544121.1598280.18315690.7746789518818565
356478333780040.354892683.6833623.4449880.063051690.2569398572226654
356485333780040.374132284.3351721.283380.0605397450.8281109722104576
356486333780040.375515984.8588622.1162220.1234537440.8691881135544401
356526333780040.5247600784.8692921.2909160.136301140.29983901925054285

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

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.

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.

[21]:

table.meta

[21]:

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

[22]:

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)
        + ";"
    )

[23]:

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

COLUMN_PHASE: PHASE; PINT_VERS: 0.9.8; GAMMAPY_VERS: 1.0.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: 60396.4769125122;

[24]:

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

[25]:

table.meta

[25]:

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: 0.9.8; GAMMAPY_VERS: 1.0.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: 60396.4769125122;')])

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

[26]:

# 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.
'THETA' axis is stored as a scalar -- converting to 1D array.

[27]:

new_obs.events.table

[27]:
Table length=11189
EVENT_IDTIMERADECENERGYPHASE
sdegdegTeV
int64float64float32float32float32float64
2402333778852.509924984.5945722.030880.181946010.39349195931383935
2408333778852.526715384.2146223.449140.083973940.8919526790064667
2434333778852.6131505483.52470422.7257920.105969320.4579789031294724
2445333778852.669014283.7695722.4510060.197334980.11641860219368935
2478333778852.762793983.47851623.4845940.085222190.9004806857711186
2481333778852.777854983.7151721.9851151.00209430.3476011689926769
2513333778852.864446782.42119622.5676520.143740680.9182740999289545
2544333778852.982606483.6413622.0413150.103166290.42611249421537956
2559333778853.026941484.06917622.973370.0471844930.7422975122115125
..................
356222333780039.452039784.7448220.8949810.0433121620.4537121484157113
356223333780039.460049284.1161522.5575050.081100820.6914920486374702
356227333780039.4710536683.4153421.673440.20963620.018183944316962984
356242333780039.517909583.5516522.7729850.176728350.40920620046641115
356282333780039.6299751484.4613321.693570.050687180.7361270955486778
356473333780040.338647984.4544121.1598280.18315690.7746789518818565
356478333780040.354892683.6833623.4449880.063051690.2569398572226654
356485333780040.374132284.3351721.283380.0605397450.8281109722104576
356486333780040.375515984.8588622.1162220.1234537440.8691881135544401
356526333780040.5247600784.8692921.2909160.136301140.29983901925054285

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 facilitate the manipulation of the HDU table.

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.

[28]:

data_store.hdu_table.base_dir

[28]:

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

[29]:

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

[30]:

output_path

[30]:

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

[31]:

# Save the observation object in the specified file_path
print("Writing outputfile in " + str(file_path))
observation.events.write(
    filename=file_path, gti=observation.gti, overwrite=True
)

Writing outputfile in /home/runner/work/gammapy-recipes/gammapy-recipes/gammapy-datasets/1.0.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.

[32]:

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

[32]:
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
[33]:

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

[34]:

new_hdu

[34]:
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 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.

[35]:

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

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.

5. Opening the new DataStore

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 :

[36]:

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

[37]:

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

[37]:

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

[38]:

observations[0].events.table

[38]:
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
...............
356222333780039.452039784.7448220.8949810.043312162
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

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

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