.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/astrophysics/pulsar_analysis.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_astrophysics_pulsar_analysis.py: Pulsar analysis =============== Produce a phasogram, phased-resolved maps and spectra for pulsar analysis. Introduction ------------ This notebook shows how to do a simple pulsar analysis with Gammapy. We will produce a phasogram, a phase-resolved map and a phase-resolved spectrum of the Vela pulsar. In order to build these products, we will use the `~gammapy.makers.PhaseBackgroundMaker` which takes into account the on and off phase to compute a `~gammapy.datasets.MapDatasetOnOff` and a `~gammapy.datasets.SpectrumDatasetOnOff` in the phase space. This tutorial uses a simulated run of vel observation from the CTA DC1, which already contains a column for the pulsar phases. The phasing in itself is therefore not show here. It requires specific packages like Tempo2 or `PINT `__. A gammapy recipe shows how to compute phases with PINT in the framework of Gammapy. Opening the data ---------------- Let’s first do the imports and load the only observation containing Vela in the CTA 1DC dataset shipped with Gammapy. .. GENERATED FROM PYTHON SOURCE LINES 30-57 .. code-block:: Python import warnings import numpy as np import astropy.units as u from astropy.coordinates import SkyCoord import matplotlib.pyplot as plt # %matplotlib inline from IPython.display import display from gammapy.data import DataStore from gammapy.datasets import Datasets, FluxPointsDataset, MapDataset, SpectrumDataset from gammapy.estimators import ExcessMapEstimator, FluxPointsEstimator from gammapy.makers import ( MapDatasetMaker, PhaseBackgroundMaker, SafeMaskMaker, SpectrumDatasetMaker, ) from gammapy.maps import MapAxis, RegionGeom, WcsGeom from gammapy.modeling import Fit from gammapy.modeling.models import PowerLawSpectralModel, SkyModel from gammapy.stats import WStatCountsStatistic from gammapy.utils.regions import SphericalCircleSkyRegion warnings.filterwarnings("ignore") .. GENERATED FROM PYTHON SOURCE LINES 58-59 Load the data store (which is a subset of CTA-DC1 data): .. GENERATED FROM PYTHON SOURCE LINES 59-64 .. code-block:: Python data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps") .. GENERATED FROM PYTHON SOURCE LINES 65-66 Define observation ID and print events: .. GENERATED FROM PYTHON SOURCE LINES 66-73 .. code-block:: Python id_obs_vela = [111630] obs_list_vela = data_store.get_observations(id_obs_vela) print(obs_list_vela[0].events) .. rst-class:: sphx-glr-script-out .. code-block:: none EventList --------- Instrument : None Telescope : CTA Obs. ID : 111630 Number of events : 101430 Event rate : 56.350 1 / s Time start : 59300.833333333336 Time stop : 59300.854166666664 Min. energy : 3.00e-02 TeV Max. energy : 1.52e+02 TeV Median energy : 1.00e-01 TeV Max. offset : 5.0 deg .. GENERATED FROM PYTHON SOURCE LINES 74-76 Now that we have our observation, let’s select the events in 0.2° radius around the pulsar position. .. GENERATED FROM PYTHON SOURCE LINES 76-87 .. code-block:: Python pos_target = SkyCoord(ra=128.836 * u.deg, dec=-45.176 * u.deg, frame="icrs") on_radius = 0.2 * u.deg on_region = SphericalCircleSkyRegion(pos_target, on_radius) # Apply angular selection events_vela = obs_list_vela[0].events.select_region(on_region) print(events_vela) .. rst-class:: sphx-glr-script-out .. code-block:: none EventList --------- Instrument : None Telescope : CTA Obs. ID : 111630 Number of events : 843 Event rate : 0.468 1 / s Time start : 59300.833333333336 Time stop : 59300.854166666664 Min. energy : 3.00e-02 TeV Max. energy : 4.33e+01 TeV Median energy : 1.07e-01 TeV Max. offset : 1.7 deg .. GENERATED FROM PYTHON SOURCE LINES 88-89 Let’s load the phases of the selected events in a dedicated array. .. GENERATED FROM PYTHON SOURCE LINES 89-97 .. code-block:: Python phases = events_vela.table["PHASE"] # Let's take a look at the first 10 phases display(phases[:10]) .. rst-class:: sphx-glr-script-out .. code-block:: none PHASE ----------- 0.81847286 0.45646095 0.111507416 0.43416595 0.76837444 0.3639946 0.58693695 0.51095676 0.5606985 0.2505703 .. GENERATED FROM PYTHON SOURCE LINES 98-105 Phasogram --------- Once we have the phases, we can make a phasogram. A phasogram is a histogram of phases. It works exactly like any other histogram (you can set the binning, evaluate the errors based on the counts in each bin, etc). .. GENERATED FROM PYTHON SOURCE LINES 105-136 .. code-block:: Python nbins = 30 phase_min, phase_max = (0, 1) values, bin_edges = np.histogram(phases, range=(phase_min, phase_max), bins=nbins) bin_width = (phase_max - phase_min) / nbins bin_center = (bin_edges[:-1] + bin_edges[1:]) / 2 # Poissonian uncertainty on each bin values_err = np.sqrt(values) fig, ax = plt.subplots() ax.bar( x=bin_center, height=values, width=bin_width, color="orangered", alpha=0.7, edgecolor="black", yerr=values_err, ) ax.set_xlim(0, 1) ax.set_xlabel("Phase") ax.set_ylabel("Counts") ax.set_title(f"Phasogram with angular cut of {on_radius}") plt.show() on_phase_range = (0.5, 0.6) off_phase_range = (0.7, 1) .. image-sg:: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_001.png :alt: Phasogram with angular cut of 0.2 deg :srcset: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 137-139 Now let’s add some fancy additions to our phasogram: a patch on the ON- and OFF-phase regions and one for the background level. .. GENERATED FROM PYTHON SOURCE LINES 139-152 .. code-block:: Python # Evaluate background level mask_off = (off_phase_range[0] < phases) & (phases < off_phase_range[1]) count_bkg = mask_off.sum() print(f"Number of Off events: {count_bkg}") # bkg level normalized by the size of the OFF zone (0.3) bkg = count_bkg / nbins / (off_phase_range[1] - off_phase_range[0]) # error on the background estimation bkg_err = np.sqrt(count_bkg) / nbins / (off_phase_range[1] - off_phase_range[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of Off events: 234 .. GENERATED FROM PYTHON SOURCE LINES 153-154 Let's redo the same plot for the basis .. GENERATED FROM PYTHON SOURCE LINES 154-202 .. code-block:: Python fig, ax = plt.subplots(figsize=(10, 7)) ax.bar( x=bin_center, height=values, width=bin_width, color="orangered", alpha=0.7, edgecolor="black", yerr=values_err, ) # Plot background level x_bkg = np.linspace(0, 1, 50) kwargs = {"color": "black", "alpha": 0.7, "ls": "--", "lw": 2} ax.plot(x_bkg, (bkg - bkg_err) * np.ones_like(x_bkg), **kwargs) ax.plot(x_bkg, (bkg + bkg_err) * np.ones_like(x_bkg), **kwargs) ax.fill_between( x_bkg, bkg - bkg_err, bkg + bkg_err, facecolor="grey", alpha=0.5 ) # grey area for the background level # Let's make patches for the on and off phase zones on_patch = ax.axvspan( on_phase_range[0], on_phase_range[1], alpha=0.5, color="royalblue", ec="black" ) off_patch = ax.axvspan( off_phase_range[0], off_phase_range[1], alpha=0.25, color="white", hatch="x", ec="black", ) # Legends "ON" and "OFF" ax.text(0.55, 5, "ON", color="black", fontsize=17, ha="center") ax.text(0.895, 5, "OFF", color="black", fontsize=17, ha="center") ax.set_xlabel("Phase") ax.set_ylabel("Counts") ax.set_xlim(0, 1) ax.set_title(f"Phasogram with angular cut of {on_radius}") plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_002.png :alt: Phasogram with angular cut of 0.2 deg :srcset: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-207 Make a Li&Ma test over the events --------------------------------- Another thing that we want to do is to compute a Li&Ma test between the on-phase and the off-phase. .. GENERATED FROM PYTHON SOURCE LINES 207-230 .. code-block:: Python # Calculate the ratio between the on-phase and the off-phase alpha = (on_phase_range[1] - on_phase_range[0]) / ( off_phase_range[1] - off_phase_range[0] ) # Select events in the on region region_events = obs_list_vela[0].events.select_region(on_region) # Select events in phase space on_events = region_events.select_parameter("PHASE", band=on_phase_range) off_events = region_events.select_parameter("PHASE", band=off_phase_range) # Apply the WStat (Li&Ma statistic) pulse_stat = WStatCountsStatistic( len(on_events.time), len(off_events.time), alpha=alpha ) print(f"Number of excess counts: {pulse_stat.n_sig}") print(f"TS: {pulse_stat.ts}") print(f"Significance: {pulse_stat.sqrt_ts}") .. rst-class:: sphx-glr-script-out .. code-block:: none Number of excess counts: 44.00000000000003 TS: 15.211770556360534 Significance: 3.9002269877996247 .. GENERATED FROM PYTHON SOURCE LINES 231-233 Phase-resolved map ------------------ .. GENERATED FROM PYTHON SOURCE LINES 236-240 Now that we have an overview of the phasogram of the pulsar, we can do a phase-resolved sky map : a map of the ON-phase events minus alpha times the OFF-phase events. Alpha is the ratio between the size of the ON-phase zone (here 0.1) and the OFF-phase zone (0.3). .. GENERATED FROM PYTHON SOURCE LINES 240-253 .. code-block:: Python e_true = MapAxis.from_energy_bounds( 0.003, 10, 6, per_decade=True, unit="TeV", name="energy_true" ) e_reco = MapAxis.from_energy_bounds( 0.01, 10, 4, per_decade=True, unit="TeV", name="energy" ) geom = WcsGeom.create( binsz=0.02 * u.deg, skydir=pos_target, width="4 deg", axes=[e_reco] ) .. GENERATED FROM PYTHON SOURCE LINES 254-255 Let’s create an ON-map and an OFF-map: .. GENERATED FROM PYTHON SOURCE LINES 255-275 .. code-block:: Python map_dataset_empty = MapDataset.create(geom=geom, energy_axis_true=e_true) map_dataset_maker = MapDatasetMaker() phase_bkg_maker = PhaseBackgroundMaker( on_phase=on_phase_range, off_phase=off_phase_range, phase_column_name="PHASE" ) offset_max = 5 * u.deg safe_mask_maker = SafeMaskMaker(methods=["offset-max"], offset_max=offset_max) map_datasets = Datasets() for obs in obs_list_vela: map_dataset = map_dataset_maker.run(map_dataset_empty, obs) map_dataset = safe_mask_maker.run(map_dataset, obs) map_dataset_on_off = phase_bkg_maker.run(map_dataset, obs) map_datasets.append(map_dataset_on_off) .. GENERATED FROM PYTHON SOURCE LINES 276-280 Once the data reduction is done, we can plot the map of the counts-ON (i.e. in the ON-phase) and the map of the background (i.e. the counts-OFF, selected in the OFF-phase, multiplied by alpha). If one wants to plot the counts-OFF instead, ``background`` should be replaced by ``counts_off`` in the following cell. .. GENERATED FROM PYTHON SOURCE LINES 280-301 .. code-block:: Python counts = ( map_datasets[0].counts.smooth(kernel="gauss", width=0.1 * u.deg).sum_over_axes() ) background = ( map_datasets[0].background.smooth(kernel="gauss", width=0.1 * u.deg).sum_over_axes() ) fig, (ax1, ax2) = plt.subplots( figsize=(11, 4), ncols=2, subplot_kw={"projection": counts.geom.wcs} ) counts.plot(ax=ax1, add_cbar=True) ax1.set_title("Counts") background.plot(ax=ax2, add_cbar=True) ax2.set_title("Background") plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_003.png :alt: Counts, Background :srcset: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 302-303 Finally, we can run an `~gammapy.estimators.ExcessMapEstimator` to compute the excess and significance maps. .. GENERATED FROM PYTHON SOURCE LINES 303-325 .. code-block:: Python excess_map_estimator = ExcessMapEstimator( correlation_radius="0.2 deg", energy_edges=[50 * u.GeV, 10 * u.TeV] ) estimator_results = excess_map_estimator.run(dataset=map_datasets[0]) npred_excess = estimator_results.npred_excess sqrt_ts = estimator_results.sqrt_ts fig, (ax1, ax2) = plt.subplots( figsize=(11, 4), ncols=2, subplot_kw={"projection": npred_excess.geom.wcs} ) npred_excess.plot(ax=ax1, add_cbar=True) ax1.set_title("Excess counts") sqrt_ts.plot(ax=ax2, add_cbar=True) ax2.set_title("Significance") plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_004.png :alt: Excess counts, Significance :srcset: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 326-340 Note that here we are lacking statistic because we only use one run of CTAO. Phase-resolved spectrum ----------------------- We can also make a phase-resolved spectrum. In order to do that, we are going to use the `~gammapy.makers.PhaseBackgroundMaker` to create a `~gammapy.datasets.SpectrumDatasetOnOff` with the ON and OFF taken in the phase space. Note that this maker take the ON and OFF in the same spatial region. Here to create the `~gammapy.datasets.SpectrumDatasetOnOff`, we are going to redo the whole data reduction. However, note that one can use the `~gammapy.datasets.MapDatasetOnOff.to_spectrum_dataset()` method (with the ``containment_correction`` parameter set to True) if such a `~gammapy.datasets.MapDatasetOnOff` has been created as shown above. .. GENERATED FROM PYTHON SOURCE LINES 340-366 .. code-block:: Python e_true = MapAxis.from_energy_bounds(0.003, 10, 100, unit="TeV", name="energy_true") e_reco = MapAxis.from_energy_bounds(0.01, 10, 30, unit="TeV", name="energy") geom = RegionGeom.create(region=on_region, axes=[e_reco]) spectrum_dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=e_true) spectrum_dataset_maker = SpectrumDatasetMaker() phase_bkg_maker = PhaseBackgroundMaker( on_phase=on_phase_range, off_phase=off_phase_range, phase_column_name="PHASE" ) offset_max = 5 * u.deg safe_mask_maker = SafeMaskMaker(methods=["offset-max"], offset_max=offset_max) spectrum_datasets = Datasets() for obs in obs_list_vela: spectrum_dataset = spectrum_dataset_maker.run(spectrum_dataset_empty, obs) spectrum_dataset = safe_mask_maker.run(spectrum_dataset, obs) spectrum_dataset_on_off = phase_bkg_maker.run(spectrum_dataset, obs) spectrum_datasets.append(spectrum_dataset_on_off) .. GENERATED FROM PYTHON SOURCE LINES 367-368 Now let’s take a look at the datasets we just created: .. GENERATED FROM PYTHON SOURCE LINES 368-373 .. code-block:: Python spectrum_datasets[0].peek() plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_005.png :alt: Counts, Exposure, Energy Dispersion :srcset: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 374-378 Now we’ll fit a model to the spectrum with the `~gammapy.modeling.Fit` class. First we load a power law model with an initial value for the index and the amplitude and then wo do a likelihood fit. The fit results are printed below. .. GENERATED FROM PYTHON SOURCE LINES 378-397 .. code-block:: Python spectral_model = PowerLawSpectralModel( index=4, amplitude="1.3e-9 cm-2 s-1 TeV-1", reference="0.02 TeV" ) model = SkyModel(spectral_model=spectral_model, name="vela psr") emin_fit, emax_fit = (0.04 * u.TeV, 0.4 * u.TeV) mask_fit = geom.energy_mask(energy_min=emin_fit, energy_max=emax_fit) for dataset in spectrum_datasets: dataset.models = model dataset.mask_fit = mask_fit joint_fit = Fit() joint_result = joint_fit.run(datasets=spectrum_datasets) print(joint_result) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 101 total stat : 7.07 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 398-401 Now you might want to do the stacking here even if in our case there is only one observation which makes it superfluous. We can compute flux points by fitting the norm of the global model in energy bands. .. GENERATED FROM PYTHON SOURCE LINES 401-423 .. code-block:: Python energy_edges = np.logspace(np.log10(0.04), np.log10(0.4), 7) * u.TeV stack_dataset = spectrum_datasets.stack_reduce() stack_dataset.models = model fpe = FluxPointsEstimator( energy_edges=energy_edges, source="vela psr", selection_optional="all" ) flux_points = fpe.run(datasets=[stack_dataset]) flux_points.meta["ts_threshold_ul"] = 1 amplitude_ref = 0.57 * 19.4e-14 * u.Unit("1 / (cm2 s MeV)") spec_model_true = PowerLawSpectralModel( index=4.5, amplitude=amplitude_ref, reference="20 GeV" ) flux_points_dataset = FluxPointsDataset(data=flux_points, models=model) .. GENERATED FROM PYTHON SOURCE LINES 424-425 Now we can plot. .. GENERATED FROM PYTHON SOURCE LINES 425-444 .. code-block:: Python ax_spectrum, ax_residuals = flux_points_dataset.plot_fit() ax_spectrum.set_ylim([1e-14, 3e-11]) ax_residuals.set_ylim([-1.7, 1.7]) spec_model_true.plot( ax=ax_spectrum, energy_bounds=(emin_fit, emax_fit), label="Reference model", c="black", linestyle="dashed", energy_power=2, ) ax_spectrum.legend(loc="best") plt.show() .. image-sg:: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_006.png :alt: pulsar analysis :srcset: /tutorials/astrophysics/images/sphx_glr_pulsar_analysis_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 445-449 This tutorial suffers a bit from the lack of statistics: there were 9 Vela observations in the CTA DC1 while there is only one here. When done on the 9 observations, the spectral analysis is much better agreement between the input model and the gammapy fit. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.765 seconds) .. _sphx_glr_download_tutorials_astrophysics_pulsar_analysis.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v2.0?urlpath=lab/tree/notebooks/2.0/tutorials/astrophysics/pulsar_analysis.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: pulsar_analysis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: pulsar_analysis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: pulsar_analysis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_