.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-time/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_analysis-time_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-59 .. code-block:: Python # Remove warnings 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 from gammapy.utils.check import check_tutorials_setup warnings.filterwarnings("ignore") .. GENERATED FROM PYTHON SOURCE LINES 60-62 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 62-66 .. code-block:: Python check_tutorials_setup() .. rst-class:: sphx-glr-script-out .. code-block:: none System: python_executable : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/bin/python python_version : 3.9.20 machine : x86_64 system : Linux Gammapy package: version : 1.3 path : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.13.1 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.2 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.2 pandas : not installed healpy : 1.17.3 iminuit : 2.30.1 sherpa : 4.16.1 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.3 ray : 2.39.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3 .. GENERATED FROM PYTHON SOURCE LINES 67-68 Load the data store (which is a subset of CTA-DC1 data): .. GENERATED FROM PYTHON SOURCE LINES 68-73 .. code-block:: Python data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps") .. GENERATED FROM PYTHON SOURCE LINES 74-75 Define observation ID and print events: .. GENERATED FROM PYTHON SOURCE LINES 75-82 .. 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 83-85 Now that we have our observation, let’s select the events in 0.2° radius around the pulsar position. .. GENERATED FROM PYTHON SOURCE LINES 85-96 .. 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 97-98 Let’s load the phases of the selected events in a dedicated array. .. GENERATED FROM PYTHON SOURCE LINES 98-106 .. 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 107-114 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 114-145 .. 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/analysis-time/images/sphx_glr_pulsar_analysis_001.png :alt: Phasogram with angular cut of 0.2 deg :srcset: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 146-148 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 148-161 .. 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 162-163 Let's redo the same plot for the basis .. GENERATED FROM PYTHON SOURCE LINES 163-211 .. 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/analysis-time/images/sphx_glr_pulsar_analysis_002.png :alt: Phasogram with angular cut of 0.2 deg :srcset: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 212-216 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 216-239 .. 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 240-242 Phase-resolved map ------------------ .. GENERATED FROM PYTHON SOURCE LINES 245-249 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 249-262 .. 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 263-264 Let’s create an ON-map and an OFF-map: .. GENERATED FROM PYTHON SOURCE LINES 264-284 .. 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 285-288 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 288-309 .. 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/analysis-time/images/sphx_glr_pulsar_analysis_003.png :alt: Counts, Background :srcset: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 310-311 Finally, we can run an `~gammapy.estimators.ExcessMapEstimator` to compute the excess and significance maps. .. GENERATED FROM PYTHON SOURCE LINES 311-333 .. 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/analysis-time/images/sphx_glr_pulsar_analysis_004.png :alt: Excess counts, Significance :srcset: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 334-348 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 348-374 .. 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 375-376 Now let’s take a look at the datasets we just created: .. GENERATED FROM PYTHON SOURCE LINES 376-381 .. code-block:: Python spectrum_datasets[0].peek() plt.show() .. image-sg:: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_005.png :alt: Counts, Exposure, Energy Dispersion :srcset: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 382-386 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 386-405 .. 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 406-409 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 409-431 .. 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 432-433 Now we can plot. .. GENERATED FROM PYTHON SOURCE LINES 433-452 .. 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/analysis-time/images/sphx_glr_pulsar_analysis_006.png :alt: pulsar analysis :srcset: /tutorials/analysis-time/images/sphx_glr_pulsar_analysis_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 453-457 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. .. _sphx_glr_download_tutorials_analysis-time_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/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/analysis-time/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 `_