.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/data/veritas.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_data_veritas.py: VERITAS with Gammapy ==================== Explore VERITAS point-like DL3 files, including event lists and IRFs and calculate Li & Ma significance, spectra, and fluxes. Introduction ------------ `VERITAS `__ (Very Energetic Radiation Imaging Telescope Array System) is a ground-based gamma-ray instrument operating at the Fred Lawrence Whipple Observatory (FLWO) in southern Arizona, USA. It is an array of four 12m optical reflectors for gamma-ray astronomy in the ~ 100 GeV to > 30 TeV energy range. VERITAS data are private and lower level analysis is done using either the `Eventdisplay `__ or `VEGAS (internal access only) `__ analysis packages to produce DL3 files (using `V2DL3 `__), which can be used in Gammapy to produce high-level analysis products. A small sub-set of archival Crab nebula data has been publicly released to accompany this tutorial, which provides an introduction to VERITAS data analysis using gammapy for VERITAS members and external users alike. This notebook is only intended for use with these publicly released Crab nebula files and the use of other sources or datasets may require modifications to this notebook. .. GENERATED FROM PYTHON SOURCE LINES 35-59 .. code-block:: Python import numpy as np from matplotlib import pyplot as plt import astropy.units as u from gammapy.maps import MapAxis, WcsGeom, RegionGeom from regions import CircleSkyRegion, PointSkyRegion from gammapy.data import DataStore from gammapy.modeling.models import SkyModel, LogParabolaSpectralModel from gammapy.modeling import Fit from gammapy.datasets import Datasets, SpectrumDataset from gammapy.estimators import FluxPointsEstimator, LightCurveEstimator from gammapy.makers import ( ReflectedRegionsBackgroundMaker, SafeMaskMaker, SpectrumDatasetMaker, WobbleRegionsFinder, ) from astropy.coordinates import SkyCoord from gammapy.visualization import plot_spectrum_datasets_off_regions from gammapy.utils.regions import extract_bright_star_regions .. GENERATED FROM PYTHON SOURCE LINES 60-63 Data exploration ---------------- .. GENERATED FROM PYTHON SOURCE LINES 66-73 Load in files ~~~~~~~~~~~~~ First, we select and load VERITAS observations of the Crab Nebula. These files are processed with EventDisplay, but VEGAS analysis should be identical apart from the integration region size, which is handled by ``RAD_MAX``. .. GENERATED FROM PYTHON SOURCE LINES 73-78 .. code-block:: Python data_store = DataStore.from_dir("$GAMMAPY_DATA/veritas/crab-point-like-ED") data_store.info() .. rst-class:: sphx-glr-script-out .. code-block:: none Data store: HDU index table: BASE_DIR: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev/veritas/crab-point-like-ED Rows: 16 OBS_ID: 64080 -- 64083 HDU_TYPE: [np.str_('aeff'), np.str_('edisp'), np.str_('events'), np.str_('gti')] HDU_CLASS: [np.str_('aeff_2d'), np.str_('edisp_2d'), np.str_('events'), np.str_('gti')] Observation table: Observatory name: 'N/A' Number of observations: 4 .. GENERATED FROM PYTHON SOURCE LINES 79-83 We filter our data by only taking observations within :math:`5^\circ` of the Crab Nebula. Further details on how to filter observations can be found in :doc:`../../user-guide/dl3`. .. GENERATED FROM PYTHON SOURCE LINES 83-92 .. code-block:: Python target_position = SkyCoord(83.6333, 22.0145, unit="deg") selected_obs_table = data_store.obs_table.select_sky_circle(target_position, 5 * u.deg) obs_ids = selected_obs_table["OBS_ID"] observations = data_store.get_observations(obs_id=obs_ids, required_irf="point-like") .. GENERATED FROM PYTHON SOURCE LINES 93-96 Peek the first run in the data release ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 98-100 Peek the events and their time/energy/spatial distributions. .. GENERATED FROM PYTHON SOURCE LINES 100-104 .. code-block:: Python observations[0].events.peek() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_veritas_001.png :alt: veritas :srcset: /tutorials/data/images/sphx_glr_veritas_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 105-111 Peek at the IRFs included. You should verify that the IRFs are filled correctly and that there are no values set to zero within your analysis range. We can also peek at the effective area (``aeff``) or energy migration matrices (``edisp``) with the ``peek()`` method. .. GENERATED FROM PYTHON SOURCE LINES 111-115 .. code-block:: Python observations[0].peek(figsize=(25, 5)) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_veritas_002.png :alt: Effective area, Energy dispersion, Events :srcset: /tutorials/data/images/sphx_glr_veritas_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 116-119 Estimate counts and significance -------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 122-141 Set the energy binning ~~~~~~~~~~~~~~~~~~~~~~ The energy axis will determine the bins in which energy is calculated, while the true energy axis defines the binning of the energy dispersion matrix and the effective area. Generally, the true energy axis should be more finely binned than the energy axis and span a larger range of energies, and the energy axis should be binned to match the needs of spectral reconstruction. Note that if the `~gammapy.makers.SafeMaskMaker` (which we will define later) is set to exclude events below a given percentage of the effective area, it will remove the entire bin containing the energy that corresponds to that percentage. Additionally, spectral bins are determined based on the energy axis and cannot be finer or offset from the energy axis bin edges. See :ref:`Safe Data Range ` for more information on how the safe mask maker works. .. GENERATED FROM PYTHON SOURCE LINES 141-147 .. code-block:: Python energy_axis = MapAxis.from_energy_bounds("0.05 TeV", "100 TeV", nbin=50) energy_axis_true = MapAxis.from_energy_bounds( "0.01 TeV", "110 TeV", nbin=200, name="energy_true" ) .. GENERATED FROM PYTHON SOURCE LINES 148-155 Create an exclusion mask ~~~~~~~~~~~~~~~~~~~~~~~~ Here, we create a spatial mask and append exclusion regions for the source region and stars (< 6th magnitude) contained within the ``exclusion_geom``. We define a star exclusion region of 0.3 deg, which should contain bright stars within the VERITAS optical PSF. .. GENERATED FROM PYTHON SOURCE LINES 155-170 .. code-block:: Python exclusion_geom = WcsGeom.create( skydir=(target_position.ra.value, target_position.dec.value), binsz=0.01, width=(4, 4), frame="icrs", proj="CAR", ) source_exclusion_region = CircleSkyRegion(center=target_position, radius=0.35 * u.deg) exclusion_regions = extract_bright_star_regions(exclusion_geom) exclusion_regions.append(source_exclusion_region) exclusion_mask = ~exclusion_geom.region_mask(exclusion_regions) .. GENERATED FROM PYTHON SOURCE LINES 171-185 Define the integration region ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Point-like DL3 files can only be analyzed using the reflected regions background method and for a pre-determined integration region (which is the :math:`\sqrt{\theta^2}` used in IRF simulations). The default values for moderate/medium cuts are determined by the DL3 file’s ``RAD_MAX`` keyword. For VERITAS data (ED and VEGAS), ``RAD_MAX`` is not energy dependent. Note that full-enclosure files are required to use any non-point-like integration region. .. GENERATED FROM PYTHON SOURCE LINES 185-189 .. code-block:: Python on_region = PointSkyRegion(target_position) geom = RegionGeom.create(region=on_region, axes=[energy_axis]) .. GENERATED FROM PYTHON SOURCE LINES 190-205 Define the `~gammapy.makers.SafeMaskMaker` ------------------------------------------ The `~gammapy.makers.SafeMaskMaker` sets the boundaries of our analysis based on the uncertainties contained in the instrument response functions (IRFs). For VERITAS point-like analysis (both ED and VEGAS), the following methods are strongly recommended: * ``offset-max``: Sets the maximum radial offset from the camera center within which we accept events. This is set to the edge of the VERITAS FoV. * ``edisp-bias``: Removes events which are reconstructed with energies that have :math:`>5\%` energy bias. * ``aeff-max``: Removes events which are reconstructed to :math:`<10\%` of the maximum value of the effective area. These are important to remove for spectral analysis, since they have large uncertainties on their reconstructed energies. .. GENERATED FROM PYTHON SOURCE LINES 205-214 .. code-block:: Python safe_mask_maker = SafeMaskMaker( methods=["offset-max", "aeff-max", "edisp-bias"], aeff_percent=10, bias_percent=5, offset_max=1.75 * u.deg, ) .. GENERATED FROM PYTHON SOURCE LINES 215-227 Data reduction -------------- We will now run the data reduction chain to calculate our ON and OFF counts. To get a significance for the whole energy range (to match VERITAS packages), remove the `~gammapy.makers.SafeMaskMaker` from being applied to ``dataset_on_off``. The parameters of the reflected background regions can be changed using the `~gammapy.makers.WobbleRegionsFinder`, which is passed as an argument to the `~gammapy.makers.ReflectedRegionsBackgroundMaker`. .. GENERATED FROM PYTHON SOURCE LINES 227-245 .. code-block:: Python dataset_maker = SpectrumDatasetMaker(selection=["counts", "exposure", "edisp"]) dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true) region_finder = WobbleRegionsFinder(n_off_regions=16) bkg_maker = ReflectedRegionsBackgroundMaker( exclusion_mask=exclusion_mask, region_finder=region_finder ) datasets = Datasets() for obs_id, observation in zip(obs_ids, observations): dataset = dataset_maker.run(dataset_empty.copy(name=str(obs_id)), observation) dataset = safe_mask_maker.run(dataset, observation) dataset_on_off = bkg_maker.run(dataset, observation) datasets.append(dataset_on_off) .. GENERATED FROM PYTHON SOURCE LINES 246-251 The plot below will show your exclusion regions in black and the center of your background regions with coloured stars. You should check to make sure these regions are sensible and that none of your background regions overlap with your exclusion regions. .. GENERATED FROM PYTHON SOURCE LINES 251-258 .. code-block:: Python plt.figure() ax = exclusion_mask.plot() plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_veritas_003.png :alt: veritas :srcset: /tutorials/data/images/sphx_glr_veritas_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.11/site-packages/gammapy/visualization/datasets.py:84: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. handle = Patch(**plot_kwargs) .. GENERATED FROM PYTHON SOURCE LINES 259-262 Significance analysis results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 265-272 Here, we display the results of the significance analysis. ``info_table`` can be modified with ``cumulative = False`` to display a table with rows that correspond to the values for each run separately. However, ``cumulative = True`` is needed to produce the combined values in the next cell. .. GENERATED FROM PYTHON SOURCE LINES 272-282 .. code-block:: Python info_table = datasets.info_table(cumulative=True) print(info_table) print(f"ON: {info_table['counts'][-1]}") print(f"OFF: {info_table['counts_off'][-1]}") print(f"Significance: {info_table['sqrt_ts'][-1]:.2f} sigma") print(f"Alpha: {info_table['alpha'][-1]:.2f}") .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess ... acceptance_off alpha ... ------- ------ ------------------ ... ------------------ ------------------ stacked 106 102.85713958740234 ... 560.0 0.0714285746216774 stacked 255 249.85714721679688 ... 1120.0 0.0714285746216774 stacked 383 375.78570556640625 ... 1679.9998779296875 0.0714285746216774 stacked 517 507.78570556640625 ... 2239.999755859375 0.071428582072258 ON: 517 OFF: 129 Significance: 46.60 sigma Alpha: 0.07 .. GENERATED FROM PYTHON SOURCE LINES 283-288 We can also plot the cumulative excess counts and significance over time. For a steady source, we generally expect excess to increase linearly with time and for significance to increase as :math:`\sqrt{\textrm{time}}`. .. GENERATED FROM PYTHON SOURCE LINES 288-312 .. code-block:: Python fig, (ax_excess, ax_sqrt_ts) = plt.subplots(figsize=(10, 4), ncols=2, nrows=1) ax_excess.plot( info_table["livetime"].to("h"), info_table["excess"], marker="o", ) ax_excess.set_title("Excess") ax_excess.set_xlabel("Livetime [h]") ax_excess.set_ylabel("Excess events") ax_sqrt_ts.plot( info_table["livetime"].to("h"), info_table["sqrt_ts"], marker="o", ) ax_sqrt_ts.set_title("Significance") ax_sqrt_ts.set_xlabel("Livetime [h]") ax_sqrt_ts.set_ylabel("Significance [sigma]") plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_veritas_004.png :alt: Excess, Significance :srcset: /tutorials/data/images/sphx_glr_veritas_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 313-316 Make a spectrum --------------- .. GENERATED FROM PYTHON SOURCE LINES 318-327 Now, we’ll calculate the source spectrum. This uses a forward-folding approach that will assume a given spectrum and fit the counts calculated above to that spectrum in each energy bin specified by the ``energy_axis``. For this reason, it’s important that spectral model be set as closely as possible to the expected spectrum - for the Crab nebula, this is a `~gammapy.modeling.models.LogParabolaSpectralModel`. .. GENERATED FROM PYTHON SOURCE LINES 327-343 .. code-block:: Python spectral_model = LogParabolaSpectralModel( amplitude=3.75e-11 * u.Unit("cm-2 s-1 TeV-1"), alpha=2.45, beta=0.15, reference=1 * u.TeV, ) model = SkyModel(spectral_model=spectral_model, name="crab") datasets.models = [model] fit_joint = Fit() result_joint = fit_joint.run(datasets=datasets) .. GENERATED FROM PYTHON SOURCE LINES 344-346 The best-fit spectral parameters are shown in this table. .. GENERATED FROM PYTHON SOURCE LINES 346-350 .. code-block:: Python print(datasets.models.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value unit ... min max frozen link prior ----- ---- --------- ---------- -------------- ... --- --- ------ ---- ----- crab amplitude 4.0105e-11 TeV-1 s-1 cm-2 ... nan nan False crab reference 1.0000e+00 TeV ... nan nan True crab alpha 2.5289e+00 ... nan nan False crab beta 2.6446e-01 ... nan nan False .. GENERATED FROM PYTHON SOURCE LINES 351-354 We can inspect how well our data fit the model’s predicted counts in each energy bin. .. GENERATED FROM PYTHON SOURCE LINES 354-360 .. code-block:: Python ax_spectrum, ax_residuals = datasets[0].plot_fit() ax_spectrum.set_ylim(0.1, 40) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_veritas_005.png :alt: veritas :srcset: /tutorials/data/images/sphx_glr_veritas_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 361-372 We can now calculate flux points to get a spectrum by fitting the ``result_joint`` model’s amplitude in selected energy bands (defined by ``energy_edges``). We set ``selection_optional = "all"`` in `~gammapy.estimators.FluxPointsEstimator`, which will include a calculation for the upper limits in bins with a significance :math:`< 2\sigma`. In the case of a non-detection or to obtain better upper limits, consider expanding the scan range for the norm parameter in `~gammapy.estimators.FluxPointsEstimator`. See :doc:`/tutorials/details/estimators` for more details on how to do this. .. GENERATED FROM PYTHON SOURCE LINES 372-381 .. code-block:: Python fpe = FluxPointsEstimator( energy_edges=np.logspace(-0.7, 1.5, 12) * u.TeV, source="crab", selection_optional="all", ) flux_points = fpe.run(datasets=datasets) .. GENERATED FROM PYTHON SOURCE LINES 382-384 Now, we can plot our flux points along with the best-fit spectral model. .. GENERATED FROM PYTHON SOURCE LINES 384-392 .. code-block:: Python ax = flux_points.plot() spectral_model.plot(ax=ax, energy_bounds=(0.1, 30) * u.TeV) spectral_model.plot_error(ax=ax, energy_bounds=(0.1, 30) * u.TeV) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_veritas_006.png :alt: veritas :srcset: /tutorials/data/images/sphx_glr_veritas_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 393-396 Make a lightcurve and calculate integral flux --------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 399-407 Integral flux can be calculated by integrating the spectral model we fit earlier. This will be a model-dependent flux estimate, so the choice of spectral model should match the data as closely as possible. ``e_min`` and ``e_max`` should be adjusted depending on the analysis requirements. Note that the actual energy threshold will use the closest bin defined by the ``energy_axis`` binning. .. GENERATED FROM PYTHON SOURCE LINES 407-419 .. code-block:: Python e_min = 0.25 * u.TeV e_max = 30 * u.TeV flux, flux_errp, flux_errn = result_joint.models["crab"].spectral_model.integral_error( e_min, e_max ) print( f"Integral flux > {e_min}: {flux.value:.2} + {flux_errp.value:.2} {flux.unit} - {flux_errn.value:.2} {flux.unit}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Integral flux > 0.25 TeV: 1.7e-10 + 8e-12 1 / (s cm2) - 8.3e-12 1 / (s cm2) .. GENERATED FROM PYTHON SOURCE LINES 420-426 Finally, we’ll create a run-wise binned light curve. See the :doc:`../analysis-time/light_curve_flare` tutorial for instructions on how to set up sub-run binning. Here, we set our energy edges so that the light curve has an energy threshold of 0.25 TeV and will plot upper limits for time bins with significance :math:`<2 \sigma`. .. GENERATED FROM PYTHON SOURCE LINES 426-435 .. code-block:: Python lc_maker = LightCurveEstimator( energy_edges=[0.25, 30] * u.TeV, source="crab", reoptimize=False ) lc_maker.n_sigma_ul = 2 lc_maker.selection_optional = ["ul"] lc = lc_maker.run(datasets) .. GENERATED FROM PYTHON SOURCE LINES 436-440 We can look at our results by printing the light curve as a table (with each line corresponding to a light curve bin) and plotting the light curve. .. GENERATED FROM PYTHON SOURCE LINES 440-448 .. code-block:: Python fig, ax = plt.subplots() lc.sqrt_ts_threshold_ul = 2 lc.plot(ax=ax, axis_name="time", sed_type="flux") plt.tight_layout() table = lc.to_table(format="lightcurve", sed_type="flux") print(table["time_min", "time_max", "flux", "flux_err"]) .. image-sg:: /tutorials/data/images/sphx_glr_veritas_007.png :alt: veritas :srcset: /tutorials/data/images/sphx_glr_veritas_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none time_min time_max ... flux_err ... 1 / (s cm2) ------------------ ------------------ ... ---------------------- 56213.45364922353 56213.467561260564 ... 1.3949943612472306e-11 56213.46796520147 56213.481877238504 ... 1.652355175420978e-11 56213.48233844736 56213.49625048439 ... 1.5247402003173758e-11 56213.496739178285 56213.51065121532 ... 1.432606578265409e-11 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 34.509 seconds) .. _sphx_glr_download_tutorials_data_veritas.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/main?urlpath=lab/tree/notebooks/dev/tutorials/data/veritas.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: veritas.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: veritas.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: veritas.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_