This page is a tutorial for the (GammaBayes)[https://github.com/lpin0002/GammaBayes] and the original notebook can be found within the docs/tutorials folder within the GitHub repository.

Exposure and GammaLogExposure

Within GammaBayes when we mention exposure we primarily mean the product between the effective area of a given telescope and observation time. However, that is not to say you can’t use the more general definition of taking the astrophysical flux and converting it into an observational flux which includes the use of masking as we will discussed later in this tutorial.

Because this will involve the effective area this tutorial will need to make use of the IRF_Loglikelihood class (which has it’s own tutorial if want to skip to that).However, for the sake of simplicitly all we need to know is that once we have an instance of the IRF_Loglikelihood class, then it has a log_aeff function that outputs the log of the effective area for the given inputs to IRF_Loglikelihood.

For instance, let’s say we have an observation run where the central pointing direction is the Galactic Centre that lasts for 5 hours and we are interested in energy values between 0.1 TeV and 100 TeV with 10 bins per decade in energy, and galactic longitude and latitudes between -4 and 4 degrees with an angular resolution of 0.02 deg then…

import numpy as np
from astropy import units as u
from gammabayes.likelihoods.irfs import IRF_LogLikelihood
from gammabayes import GammaBinning
pointing_dir = np.array([0,0])*u.deg
observation_time = 5*u.hr

bin_geometry = GammaBinning(energy_axis=np.logspace(-1,2, 31)*u.TeV, lon_axis=np.linspace(-4, 4, 401)*u.deg, lat_axis=np.linspace(-4, 4, 401)*u.deg)

# Don't worry what the difference between axes and dependent axes is yet
irf_loglike_instance = IRF_LogLikelihood(axes=bin_geometry.axes, dependent_axes=bin_geometry.axes, pointing_dir=pointing_dir, observation_time=5*u.hr)

log_aeff_function = irf_loglike_instance.log_aeff

/Users/lpin0002/anaconda3/envs/testofwest/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Let’s have a look at a slice the effective area for an energy value of 1 TeV.

spatial_mesh = np.meshgrid(*bin_geometry.spatial_axes, indexing='ij')

log_aeff_function_mesh_values = log_aeff_function(energy=(spatial_mesh[0].value.flatten()*0+1)*u.TeV,
                                                  longitude=spatial_mesh[0].flatten(),
                                                  latitude=spatial_mesh[1].flatten()).reshape(spatial_mesh[0].shape)

from matplotlib import pyplot as plt

fig = plt.figure()
plt.pcolormesh(bin_geometry.lon_axis.value, bin_geometry.lat_axis.value, np.exp(log_aeff_function_mesh_values).T)
plt.colorbar(label=r"Effective Area ["+(irf_loglike_instance.aeff_units).to_string('latex')+r"]")
plt.xlabel(r"Galactic Longitude ["+bin_geometry.lon_axis.unit.to_string('latex')+"]")
plt.ylabel(r"Galactic Latitude ["+bin_geometry.lat_axis.unit.to_string('latex')+"]")
ax = fig.get_axes()[0]
ax.invert_xaxis()

So a simple exposure would be the product of this with the observation time.

from matplotlib import pyplot as plt

fig = plt.figure()
plt.pcolormesh(bin_geometry.lon_axis.value, bin_geometry.lat_axis.value, np.exp(log_aeff_function_mesh_values).T*observation_time.to("s").value, norm='linear')
plt.colorbar(label=r"Exposure ["+(irf_loglike_instance.aeff_units*observation_time.to("s").unit).to_string('latex')+r"]")
plt.xlabel(r"Galactic Longitude ["+bin_geometry.lon_axis.unit.to_string('latex')+"]")
plt.ylabel(r"Galactic Latitude ["+bin_geometry.lat_axis.unit.to_string('latex')+"]")
ax = fig.get_axes()[0]
ax.invert_xaxis()

This quantity is pretty fundamental to GammaBayes both in terms of computation speed and purely not wanting to mess it up because issues with the values here would propagate through every stage of the calculations involving GammaBayes. So we have a dedicated class for this called GammaLogExposure.

from gammabayes import GammaLogExposure

log_exposure = GammaLogExposure(binning_geometry=bin_geometry, irfs=irf_loglike_instance,
                                pointing_dir=pointing_dir, observation_time=observation_time,
                                observation_time_unit=u.s)

• The irfs input is required so that the GammaLogExposure instance can have access to a standardized function for the effective area.

• observation_time_unit is where you can specify a time that you would like your exposures to be converted in regardless of input observation time unit. This is handy as observation runs are typically measured in hours while astrophysical flux models are typically in units of seconds (seconds is the default for observation_time_unit)

We can “peek” at the exposure using the peek method

log_exposure.peek(dpi=70)
plt.show()

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

Multiple exposures

Now let’s say you have 9 observations centred around the Galactic Centre with equal observation times with a total of 525 hours (mimicing the planned GC survey by CTAO).

pointing_dirs = [
    np.array([-1, 1])*u.deg,
    np.array([0, 1])*u.deg,
    np.array([1, 1])*u.deg,

    np.array([-1., 0.])*u.deg,
    np.array([0., 0.])*u.deg,
    np.array([1., 0.])*u.deg,

    np.array([-1, -1])*u.deg,
    np.array([0, -1])*u.deg,
    np.array([1, -1])*u.deg,
    ]
observation_times = [525/9*u.hr for obs in range(len(pointing_dirs))]

plt.figure()
for coord in pointing_dirs:
    plt.scatter(x=coord[0], y=coord[1], marker='x', c='tab:red')
plt.xlim([bin_geometry.lon_axis.value.min(), bin_geometry.lon_axis.value.max()])
plt.ylim([bin_geometry.lat_axis.value.min(), bin_geometry.lat_axis.value.max()])
plt.xlabel("Galactic Longitude [deg]")
plt.ylabel("Galactic Latitude [deg]")
plt.grid(which='both', c='grey', ls='--', alpha=0.4)
plt.show()

Each observation run would have it’s own exposure, which we can construct with the following.

log_exposures = [GammaLogExposure(binning_geometry=bin_geometry,
                                 irfs=irf_loglike_instance,
                                 pointing_dir=pointing,
                                 observation_time=obs_time,
                                 ) for pointing, obs_time in zip(pointing_dirs, observation_times)]

Quite commonly we will then want to combine these exposures into a single combined exposure, this can be done a couple ways.

Firstly, the manual way by accessing the log exposure values directly. We do not recommend this as you can have several issues with the handling of units and the like but it would go something like this.

from scipy.special import logsumexp

combined_log_exposure = logsumexp([log_exposure.log_exposure_map for log_exposure in log_exposures], axis=0)

And then I’m going to put this back into a GammaLogExposure class essentially for the plotting.

combined_log_exposure_class_instance = GammaLogExposure(binning_geometry=bin_geometry,
                                                        log_exposure_map=combined_log_exposure,
                                                        irfs=irf_loglike_instance, pointing_dir=np.array([0,0])*u.deg)

fig, ax = combined_log_exposure_class_instance.peek(pcolormesh_kwargs={"norm":"log"}, dpi=100)
for coord in pointing_dirs:
    ax[0, 1].scatter(x=coord[0], y=coord[1], marker='x', c='tab:red')
    ax[1,1].scatter(x=coord[0], y=coord[1], marker='x', c='tab:red')

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

The second way is to use the “+” behaviour of GammaLogExposure such that when an instance is added to something else, if that something else is also GammaLogExposure then it combined the two in the same manner as above. Or if you have a list of them you can just call sum. e.g.

combined_log_exposure_class_instance_2 = sum(log_exposures)

combined_log_exposure_class_instance_2.pointing_dir = np.array([0,0])*u.deg


fig, ax = combined_log_exposure_class_instance_2.peek(pcolormesh_kwargs={"norm":"log"}, figsize=(12, 6), dpi=100)
for coord in pointing_dirs:
    ax[0, 1].scatter(x=coord[0], y=coord[1], marker='x', c='tab:red')
    ax[1,1].scatter(x=coord[0], y=coord[1], marker='x', c='tab:red')

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

But the key thing we will typically use is a single observation’s exposure and extracting the log_exposure_map attribute of the above class.

The next tutorial will deal with IRFs in more detail.