# COMPTEL data analysis using GammaLib and ctools

Jürgen Knödlseder, Werner Collmar, Manon Jarry, Mark Mcconnell

## ► To cite this version:

Jürgen Knödlseder, Werner Collmar, Manon Jarry, Mark Mcconnell. COMPTEL data analysis using GammaLib and ctools. Astronomy & Astrophysics - A&A, 2022, 665, pp.A84. 10.1051/0004-6361/202243826. [hal-03739671](https://hal.science/hal-03739671v1)

HAL Id: hal-03739671

https://hal.science/hal-03739671v1

Submitted on 17 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# COMPTEL data analysis using GammaLib and ctools

J. Knödlseder<fnref n="1" />, W. Collmar<fnref n="2" />, M. Jarry<fnref n="1" />, and M. McConnell<fnref n="3" />,<fnref n="4" />

Received 20 April 2022 / Accepted 3 June 2022

## ABSTRACT

More than 20 yr after the end of NASA's Compton Gamma-Ray Observatory mission, the data collected by its Imaging Compton Telescope (COMPTEL) still provide the most comprehensive and deepest view of our Universe in MeV gamma rays. While most of the COMPTEL data are archived at NASA's High Energy Astrophysics Science Archive Research Center (HEASARC), the absence of any publicly available software for their analysis means the data cannot benefit from the scientific advances made in the field of gamma-ray astronomy at higher energies. To make this unique treasure again accessible for science, we developed open source software that enables a comprehensive and modern analysis of the archived COMPTEL telescope data. Our software is based on a dedicated plug-in to the GammaLib library, a community-developed toolbox for the analysis of astronomical gamma-ray data. We implemented high-level scripts for building science analysis workflows in ctools, a community-developed gamma-ray astronomy science analysis software framework. We describe the implementation of our software and provide the underlying algorithms. Using data from the HEASARC archive, we demonstrate that our software reproduces derived data products that were obtained in the past using the proprietary COMPTEL software. We furthermore demonstrate that our software reproduces COMPTEL science results published in the literature. This brings the COMPTEL telescope data back into life, allowing them to benefit from recent advances in gamma-ray astronomy, and gives the community a means to unveil its still hidden treasures.

**Key words.** methods: data analysis – gamma rays: general – stars: neutron – binaries: general – ISM: atoms

## 1. Introduction

The Imaging Compton Telescope (COMPTEL) was one of the four telescopes aboard NASA's Compton Gamma-Ray Observatory (CGRO) satellite, which was operated in low Earth orbit from April 1991 to May 2000 (Schönfelder et al. 1993). The telescope scrutinised the gamma-ray sky in the 0.75–30 MeV energy range, and its observations still present the most sensitive survey of the MeV sky ever performed. Since the end of the CGRO mission, our view of the gamma-ray Universe has evolved considerably, including observations of very high-energy gamma rays from a plethora of objects and the discovery of new source classes, such as novae, gamma-ray binaries, star-forming regions, globular clusters, and starburst galaxies. Today, the latest catalogue based on data from the Fermi Large Area Telescope lists more than 5000 gamma-ray sources (Abdollahi et al. 2020), and exploring the COMPTEL archive in light of this knowledge has the potential to provide new insights into the MeV sky (Collmar & Zhang 2014).

While most of the COMPTEL data accumulated during the nine-year mission are stored at NASA's High Energy Astrophysics Science Archive Research Center (HEASARC)<fnref n="1" />, no publicly available software had existed to analyse these data. During CGRO operations the COMPTEL data were analysed using the COMPTEL Processing and Analysis Software System (COMPASS) (de Vries & COMPTEL Collaboration 1994), but

this system was only available at the COMPTEL collaboration institutes and was decommissioned after the end of the mission. Nevertheless, a few science results were still obtained from COMPTEL data after the end of the mission thanks to efforts at some COMPTEL collaboration institutes to keep Linux ports of the COMPASS data-analysis system alive (Strong & Collmar 2019; Coleman & McConnell 2020).

In order to preserve the capability to analyse the unique COMPTEL archive and to make COMPTEL data analysis accessible to the astronomical community at large, we have implemented a dedicated COMPTEL plug-in for the GammaLib library, a community-developed toolbox for the analysis of astronomical gamma-ray data (Knödlseder et al. 2011). The plug-in enables a reanalysis of COMPTEL data from the HEASARC archive, including the computation of the COMPTEL instrument response functions that until now had not been publicly available. The plug-in benefits from generic data analysis capabilities provided by GammaLib, including the joint maximum-likelihood fitting of multiple energy bands and the combination of COMPTEL data with data from other gamma-ray telescopes, enabling novel analysis approaches that go beyond the capabilities of the COMPASS system.

We furthermore extended the ctools gamma-ray astronomy science analysis software (Knödlseder et al. 2016) with the addition of several dedicated Python scripts. These scripts provide basic building blocks that each perform well-defined COMPTEL data analysis tasks. These building blocks can be combined

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, CNES, 9 avenue Colonel Roche, 31028 Toulouse, Cedex 4, France
e-mail: jknodlseder@irap.omp.eu</footnote>

<footnote n="2">² Max-Planck-Institut für extraterrestrische Physik, Postfach 1603, 85740 Garching, Germany</footnote>

<footnote n="3">³ University of New Hampshire, Space Science Center, Durham, NH 03824, USA</footnote>

<footnote n="4">⁴ Southwest Research Institute, Dept. of Earth, Oceans, and Space, Durham, NH 03824, USA</footnote>

<footnote n="1">¹ https://heasarc.gsfc.nasa.gov/docs/cgro/archive/</footnote>
to create flexible COMPTEL data analysis workflows for the exploration of the MeV sky.

In this paper we present the GammaLib plug-in and the ctools COMPTEL science analysis scripts, including the algorithms that were implemented, so that any science analysis that makes use of the software has a solid reference. We demonstrate that COMPTEL data analysis products derived using GammaLib and ctools, such as data cubes and response functions, are identical to equivalent products produced with the COMPASS software. We show that the background model that is implemented in GammaLib and ctools provides a reliable description of the COMPTEL instrumental background, and we discuss potential biases and limitations. Finally, we apply the software to a number of science cases and demonstrate that our software reproduces results published in the literature. All the algorithms and results presented in this paper were produced with GammaLib and ctools version 2.0.0. The analysis scripts and data presented in this work are available for download from Zenodo<fnref n="2" />.

## 2. The COMPTEL telescope

Before describing the software implementation, we recall some COMPTEL fundamentals that are important for understanding the remainder of the paper. COMPTEL was an imaging spectrometer that was sensitive to gamma rays in the energy range 0.75–30 MeV with an energy-dependent energy and angular resolution of 5–8% (full width half maximum) and 1.7°–4.4° (full width half maximum), respectively. The telescope had a large field of view of about one steradian and was sensitive to detected gamma-ray sources at a 1–30 MeV flux level of $10^{-9}$ erg cm$^{-2}$ s$^{-1}$ within an observing time of $10^6$ s (Schönfelder et al. 1993).

COMPTEL was composed of two modular detector layers D1 and D2, separated by 158 cm, where an incident gamma-ray photon was first Compton scattered in one of the 7 modules of D1 and then eventually interacted in one of the 14 modules of D2. The Compton scatter direction $(\chi, \psi)$ was obtained from the interaction locations in D1 and D2, the Compton scattering angle $\bar{\varphi}$ was computed from the measured energy deposits $E_1$ and $E_2$ in D1 and D2, respectively, using

$$ \bar{\varphi} = \arccos \left( 1 - m_e c^2 \left( \frac{1}{E_2} - \frac{1}{E} \right) \right), \quad (1) $$

where $E = E_1 + E_2$ is the total energy deposit in the detectors, and $m_e c^2 = 0.511$ MeV is the rest energy of the electron.

The detector layers were surrounded by an active anti-coincidence shield, composed of two veto domes for each layer, that allowed for the reduction of instrumental background<fnref n="3" /> due to charged particles. A further strong discriminator of instrumental background was the time of flight (ToF) measurement, which is the time difference between the interactions in D1 and D2 (cf. Sect. 3.3).

COMPTEL data were often analysed in four standard energy bands, covering 0.75–1, 1–3, 3–10 and 10–30 MeV, and if not otherwise stated, the same energy bands will be adopted in the present paper. Event selection parameters used in this paper are discussed in Sect. 3.3.

## 3. Software implementation

### 3.1. GammaLib plug-in

The COMPTEL support was implemented in GammaLib as an instrument plug-in that provides instrument-specific implementations of abstract virtual C++ base classes defining the data format and the instrument response functions. The plug-in also comprises wrapper functions providing access to all COMPTEL-specific C++ classes through the `gammalib` Python module. All COMPTEL-specific classes begin with the letters `GCOM`.

### 3.2. COMPTEL data

COMPTEL data are available in the HEASARC archive and are grouped by so-called viewing periods with typical durations of two weeks during which the CGRO satellite had a stable pointing. In total, the HEASARC archive contains 255 exploitable viewing periods, which is about 71% of the total number of 359 viewing periods that were executed by COMPTEL<fnref n="4" />. All data are provided in FITS format.

The input data that are relevant to GammaLib for a given viewing period are event files, good time intervals, and orbit aspect data. The orbit aspect data comprise satellite orbit and telescope pointing information that are needed for coordinate transformations. COMPTEL data types are identified by a three-letter code, which is `EVP` for event files, `TIM` for good time intervals, and `OAD` for orbit aspect data. Each viewing period comprises a single `EVP` and `TIM` file, `OAD` files are provided per day. `EVP` data are handled by the GammaLib class `GCOMEEventList`, `TIM` data by `GCOMTim`, and `OAD` data by `GCOMOad` and `GCOMOads`, where `GCOMOad` implements a single `OAD` record (or superpacket), and `GCOMOads` implements a collection of `OAD` records. These data structures are combined in an unbinned COMPTEL observation, implemented by the GammaLib class `GCOMObservation`. To construct an unbinned COMPTEL observation, the input data for a viewing period can be specified either by their file names or via a so-called observation definition XML file.

The HEASARC archive mixes different versions of `EVP` files that have different levels of processing for the ToF values. GammaLib will automatically correct for these differences, assuring that ToF values accessed through GammaLib are always ToFIII values (see Appendix B).

### 3.3. Event cubes

COMPTEL data analysis is performed on binned event data using three-dimensional event cubes spanned by the Compton scatter direction $(\chi, \psi)$ and the Compton scattering angle $\bar{\varphi}$. The binning is performed for events within a given interval of total energy, and each energy interval is treated as an individual COMPTEL observation. Combining multiple energy intervals for spectral analysis can be achieved by collecting the relevant COMPTEL observations in an observation container. The same holds for the combination of multiple viewing periods for a joint maximum likelihood analysis.

The COMPTEL data space is implemented by the `GCOMDri` class, and the `GCOMDri::compute_dre` method bins the events found in an `EVP` file into an event cube. In this process, event

<footnote n="2">² https://doi.org/10.5281/zenodo.6462842</footnote>

<footnote n="3">³ Throughout this paper we denote by instrumental background all events that are not related to celestial gamma-ray emission.</footnote>

<footnote n="4">⁴ Some viewing periods in the HEASARC archive are unreadable, and data from the end of the mission were not archived. Details about issues encountered in the HEASARC COMPTEL data archive are provided in Appendix A.</footnote>
<table><caption>**Table 1.** Standard event selection parameters that are suitable for most analysis situations and that were used throughout this paper.</caption><thead><tr><td>Parameter</td><td>Minimum</td><td>Maximum</td></tr></thead><tbody><tr><td>D1 energy deposit</td><td>70 keV</td><td>20 MeV</td></tr><tr><td>D2 energy deposit</td><td>650 keV</td><td>30 MeV</td></tr><tr><td>ToF value</td><td>115</td><td>130</td></tr><tr><td>PSD value</td><td>0</td><td>110</td></tr><tr><td>Rejection flag</td><td>1</td><td>1000</td></tr><tr><td>Veto flag</td><td>0</td><td>0</td></tr></tbody></table>

Notes. ToF channel values are given in units of 0.25 ns.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Time of flight spectrum of COMPTEL events registered during ground calibration (**left**) and in orbit (**right**). One ToF channel corresponds to 0.25 ns. Downward moving photons interacting first in D1 and then in D2 lead to a peak around channel 120, while upward moving photons interacting first in D2 and then in D1 produce a peak around channel 80. In the calibration data the upward moving photons can be clearly distinguished from the downward moving photons, while in the flight data the downward moving photons are blended with the tail of the dominating upward moving photons. The figure is adapted from Knödlseder (1994).</figcaption></figure>

selection parameters are specified through the GCOMSelection class, comprising selection intervals for D1 energy deposit, D2 energy deposit, ToF channel value, pulse shape discriminator (PSD) channel value, rejection flag and veto flag. Standard event selection parameters that were used throughout this paper are summarised in Table 1.

The ToF value measures the time between the interaction of the photons in the D1 and D2 modules, and provides a strong discriminator against instrumental background. Figure 1 illustrates that the ToF distribution of the events shows two distinct features: a forward peak associated with celestial photons that first interact in D1 followed by an interaction in D2, and a background peak due to events arising from a first interaction in D2 followed by an interaction in D1. While both peaks were clearly distinguishable in calibration data taken on ground (left panel of Fig. 1), they are heavily blended in orbit due to the dominance of upward moving background events (right panel of Fig. 1). Selecting only events from the forward peak considerably improves the signal-to-noise ratio, and Fig. 1 illustrates that the minimum ToF value has a strong impact on the background discrimination, with larger values removing more background at the expense of rejecting also an increasing number of celestial photons.

To correct for the photon rejection by the ToF selection an energy-dependent correction factor is applied to the instrument

response. This ToF correction is performed internally by GammaLib (see Appendix C for details). Consequently, flux values returned by GammaLib are always corrected for ToF selection, in contrast to the COMPASS system where the correction had to be applied by the user posterior to the analysis.

An additional, but less important, event selection parameter is the PSD value, which discriminates between neutron and photon interactions in D1, with gamma rays found around channel 80 and neutron-induced events above channel 100 (Schönfelder et al. 1993). By default, the PSD selection interval is sufficiently large so that no flux correction factor needs to be applied, yet more stringent selection intervals that lead to improved performance (Collmar et al. 1997) eventually require the application of a PSD selection correction factor that is not implemented in GammaLib.

Furthermore, the Earth atmosphere being a strong gamma-ray source, an important background reduction is obtained by excluding events that may originate from the Earth atmosphere. This is achieved by requiring that event circles have a minimum angular distance $\zeta_{min}$ from the Earth horizon. In the process of event binning, this requirement is translated into a constraint on the so-called Earth horizon angle (EHA), which is defined as the angular distance between the scatter direction of an event and the Earth horizon. Specifically, only the events are retained that satisfy $EHA \ge EHA_{min}(\bar{\varphi})$ with

$$EHA_{min}(\bar{\varphi}) = \bar{\varphi}_{min} + \zeta_{min}, \qquad (2)$$

where $\bar{\varphi}_{min}$ is the lower boundary of the event cube layer that comprises the Compton scattering angle $\bar{\varphi}$. By default, $\zeta_{min} = 5^{\circ}$.

Event selection and binning in GammaLib is done per superpacket, which is a bunch of eight spacecraft telemetry packets of 16.384 s duration that define the temporal granularity on which spacecraft orbital data information is available. Orbital data per superpacket are provided by the orbit aspect data (OAD) files, and events will be used only if their times are contained in the validity interval of any of the available superpackets. Furthermore, superpackets are considered valid only if their start and stop times are fully enclosed in one of the good time intervals, specified by the TIM file<fnref n="5" />.

Finally, the GCOMSelection class also supports the specification of phase intervals, so that events can be selected according to the orbital phase of a gamma-ray binary or the rotational phase of a pulsar. Owing to the different timescales involved, the handling of orbital phases differs from that of pulsar phases. For orbital phases, the phase as a function of event time is defined via an instance of the GModelTemporalPhaseCurve class, and the phase selection is done at the level of superpackets, implying that only orbital periods that are significantly longer than the superpacket duration of 16.384 s will produce meaningful results. We illustrate this capability below by a phase-resolved analysis of the gamma-ray binary LS 5039 (cf. Sect. 4.3). For pulsar phases, the phase selection is done at the level of individual events, using pulsar ephemerides and a Solar System barycentre time correction that is applied to the onboard time of the events. Details of the implementation are given in Appendix D, and we illustrate this capability below by a phase-resolved analysis of the Crab pulsar and pulsar wind nebula (cf. Sect. 4.2.2).

<footnote n="5">⁵ We note that COMPTEL times are specified by two values: the truncated Julian days (TJD), which is defined as the number of modified Julian days (MJD) minus 40 000, and the COMPTEL tics, which are 1/8 milliseconds long. As an example, TJD=8393 and zero tics corresponds to 1991-05-17T00:00:00 UT.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Comparison of an event cube generated with GammaLib for viewing period 2.0 and the 1–3 MeV energy band (red solid) and the equivalent event cube MPE-DRE-50607 produced by the COMPASS software (blue dashed). The left panel shows the distribution of $\chi$ values, obtained by summing over all $\psi$ and $\bar{\psi}$, the middle panel shows the distribution of $\psi$ values, obtained by summing over all $\chi$ and $\bar{\psi}$, and the right panel shows the distribution of $\bar{\psi}$ values, obtained by summing over all $\chi$ and $\psi$.</figcaption></figure>

Only events for modules that are signalled as active are considered for the event binning. The detector module status is provided by the GCOMS status class that relies on a database of daily status values that is provided with GammaLib. During the COMPTEL operations a certain number of photomultiplier tubes (PMTs) of the D2 modules failed, degrading the localisation precision of the event interactions in the respective modules. One option is to exclude these modules from the analysis, reducing the number of available events by about 25%. Alternatively, events from zones around the faulty PMTs can be excluded in the analysis, allowing recovery of a large fraction of the events for the analysis. Both analysis options are implemented in GammaLib (see Appendix E) and are explored in the analysis of LS 5039 (cf. Sect. 4.3).

COMPTEL event cubes are stored as three-dimensional FITS images with the file type designation DRE. Projections of an event cube generated using GammaLib for viewing period 2.0⁶ and total energies in the interval 1–3 MeV are shown in Fig. 2 where for comparison the projections of an equivalent DRE obtained using the COMPASS software are shown. The projections of the datasets are nearly indistinguishable, illustrating that the event cubes generated by GammaLib are equivalent to those produced by the COMPASS software.

## 3.4. Instrument response function

### 3.4.1. Factorisation

The COMPTEL instrument response function, given in units of events cm² photons⁻¹, is factorised using

$$R(\chi, \psi, \bar{\psi}|\alpha, \delta, E_{\gamma}) = \frac{L \text{DRX}(\alpha, \delta) \text{DRG}(\chi, \psi, \bar{\psi}) \text{IAQ}(\bar{\psi}|\varphi_{\text{geo}}, E_{\gamma})}{T \sum_{i \in \{S\}} 7\pi r_1^2 \cos \theta_i} \quad (3)$$

where DRX(α, δ) is the exposure in units of cm² s towards a given celestial direction (α, δ), DRG(χ, ψ, $\bar{\psi}$) is a geometry function that specifies the probability that a photon that was scattered in one of the D1 modules into direction (χ, ψ) will encounter one of the active D2 modules (the term also accounts for the Earth horizon event selection, which introduces the $\bar{\psi}$ dependence),

IAQ($\bar{\psi}|\varphi_{\text{geo}}, E_{\gamma}$) is the probability that a gamma-ray photon with energy $E_{\gamma}$ being Compton scattered in D1 interacts in D2, with

$$\varphi_{\text{geo}} = \arccos(\sin \psi \sin \delta + \cos \psi \cos \delta \cos(\chi - \alpha)) \quad (4)$$

being the angular separation between (χ, ψ) and (α, δ), L is the lifetime in units of s, and T is the exposure time in units of s.

The fraction L/T is the so-called deadtime correction factor, specifying the fraction of time during which the telescope is able to register events. For COMPTEL this fraction is rather constant and was determined by van Dijk (1996) to be L/T = 0.965. This value is hardcoded in GammaLib, and automatically applied to all analysis results.

### 3.4.2. Exposure

The exposure map DRX(α, δ) is computed in GammaLib by the method GCOMDri::compute_drx using

$$\text{DRX}(\alpha, \delta) = T_{\text{sp}} \sum_{i \in \{S\}} 7\pi r_1^2 \cos \theta_i \frac{1 - \exp(-\tau / \cos \theta_i)}{1 - \exp(-\tau)}, \quad (5)$$

where $T_{\text{sp}} = 16.384$ s is the duration of a superpacket and $r_1 = 13.8$ cm is the radius of a D1 module. The $\cos \theta_i$ factor takes into account the reduction of the effective D1 surface when viewed from an off-axis direction $\theta_i$, measured as the angle between the celestial direction (α, δ) and the COMPTEL pointing direction for a given superpacket i. The fraction accounts for the increased interaction length of off-axis photons within the D1 modules, where $\tau = 0.2$ is the typical thickness of a D1 module in radiation lengths. The sum is taken over all superpackets {$S$} that satisfy the same selection criteria that were applied to the event selection (cf. Sect. 3.3).

The exposure map is given in units of cm² s and is stored as a two-dimensional FITS image with the file type designation DRX. For illustration, Fig. 3 compares projections of the DRX obtained by GammaLib for viewing period 2.0 to the projections for an equivalent DRX that was obtained by COMPASS. The projections of the datasets are indistinguishable, illustrating that the exposure maps generated by GammaLib are equivalent to those produced by the COMPASS software.

### 3.4.3. Geometry function

The geometry function DRG(χ, ψ, $\bar{\psi}$) is the probability that a photon that was scattered in one of the D1 modules towards the

<footnote n="6">⁶ We have chosen the second viewing period in the HEASARC archive to illustrate the data cube comparisons in our paper since for the first viewing period we encountered percent-level differences that are plausibly attributed to differences in our and the COMPASS processing configuration that we were not able to track down.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of an exposure map obtained with GammaLib for viewing period 2.0 (red solid) to the exposure map MPE-DRX-32302 obtained with COMPASS (blue dashed). The upper panel shows the distribution of Galactic longitude values, obtained by summing over all Galactic latitudes, and the lower panel shows the distribution of Galactic latitude values, obtained by summing over all Galactic longitudes.</figcaption></figure>

direction $(\chi, \psi)$ will encounter one of the active modules of D2. It is computed by determining the geometrical area of the shadow of the D1 modules that is cast on the D2 modules relative to the total area of all seven D1 modules as a function of the scatter direction. Optionally, zones around failed D2 module PMTs can be excluded during this computation. The geometry function also accounts for the Earth horizon event selection, which introduces a $\bar{\phi}$ dependence in the probabilities. As the position of the Earth with respect to the telescope changes continuously, the geometry function is recomputed for each superpacket, and the results are then averaged to provide an effective geometry function that applies to a given event cube. Similar to the computation of the exposure map, the superpackets $\{S\}$ considered are those that satisfy the same selection criteria that were applied to the event selection. Details are provided in Appendix F.

The computation of the geometry function is implemented in the method `GCOMDri::compute_drg` and results are stored as three-dimensional FITS images with the file type designation DRG. Figure 4 illustrates the result of the geometry function computation, again in the form of projections for viewing period 2.0. Equivalent projections for a geometry function obtained using the COMPASS system are superimposed. The projections of

the datasets are indistinguishable, illustrating that the geometry functions generated by GammaLib are equivalent to those produced by the COMPASS software.

### 3.4.4. Compton scattering probabilities

The Compton scattering probabilities $IAQ(\bar{\phi}|\varphi_{\text{geo}}, E_{\gamma})$ were computed employing the method `GCOMIaq::set` using

$$ \text{IAQ}(\bar{\phi}|\varphi_{\text{geo}}, E_{\gamma}) = \int_{\varphi_{\text{geo,min}}}^{\varphi_{\text{geo,max}}} P_{\text{eff}}(\varphi'_{\text{geo}}, E_{\gamma}) \times P_{\text{KN}}(\varphi'_{\text{geo}}, E_{\gamma}) \times \int_{E_{\gamma}^{\text{min}}}^{E_{\gamma}^{\text{max}}} R_{\text{D1}}(E_1|\hat{E}_1) R_{\text{D2}}(E_2|\hat{E}_2) J(\bar{\phi}|E_1) dE_1, \times K(\varphi'_{\text{geo}}|\varphi_{\text{geo}}, \theta) d\varphi'_{\text{geo}} \quad (6) $$

where $P_{\text{eff}}(\varphi'_{\text{geo}}, E_{\gamma})$ is an efficiency factor that is detailed in Appendix G.1, $P_{\text{KN}}(\varphi'_{\text{geo}}, E_{\gamma})$ is the contribution of the Klein-Nishina cross-section to a given bin in $\varphi'_{\text{geo}}$ that is detailed in Appendix G.2, and $R_{\text{D1}}(E_1|\hat{E}_1)$ and $R_{\text{D2}}(E_2|\hat{E}_2)$ are the energy response functions of the D1 and D2 modules, respectively, that are detailed in Appendix G.3. The transformation from $(\varphi'_{\text{geo}}, E_{\gamma})$ to $(\hat{E}_1, \hat{E}_2)$ is done using

$$ \hat{E}_2 = \frac{E_{\gamma}}{(1 - \cos \varphi'_{\text{geo}}) \frac{E_{\gamma}}{m_e c^2} + 1} \quad (7) $$

and $\hat{E}_1 = E_{\gamma} - \hat{E}_2$;

$$ E_2 = \frac{1}{2} \times \left( \sqrt{\frac{4m_e c^2 E_1}{1 - \cos \bar{\phi}}} + E_1^2 - E_1 \right) \quad (8) $$

gives the D2 energy deposit as a function of the D1 energy deposit ($E_1$) and the Compton scattering angle ($\bar{\phi}$), and

$$ J(\bar{\phi}|E_1) = \frac{m_e c^2 \sqrt{E_1} \sin \bar{\phi}}{(1 - \cos \bar{\phi})^2 \sqrt{\frac{4m_e c^2}{1 - \cos \bar{\phi}}} + E_1} \quad (9) $$

is the Jacobian for the variable transformation. Here $K(\varphi'_{\text{geo}}|\varphi_{\text{geo}}, \theta)$ is a Gaussian kernel that provides some smearing of the response in $\varphi_{\text{geo}}$ to take into account the event location uncertainties in the D1 and D2 modules. The kernel formally depends on the zenith angle $\theta$ of the incoming gamma-ray photons. However, the IAQ is assumed independent of $\theta$, and hence the kernel will be computed for an average zenith angle of $\theta = 25^{\circ}$ (see Appendix G.4 for details).

The Compton scattering probabilities are stored as a two-dimensional FITS image with the file type designation IAQ. For illustration, Fig. 5 compares projections of the Compton scattering probabilities obtained by GammaLib for the energy band 1–3 MeV to the projections for an equivalent file that was obtained using COMPASS. Both projections are indistinguishable, and the same is true for the other energy bands that were investigated. As demonstrated, the GammaLib implementation accurately reproduces the response computations that were implemented in COMPASS. We note that the implementation of the COMPTEL response computation within GammaLib is crucial for the preservation of COMPTEL data analysis capabilities, since response functions are lacking in the HEASARC archive. In addition to the internal response computation, GammaLib also includes several response functions that were derived using simulations (Stacy et al. 1996), and that can be used as an alternative to the internally computed response functions.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of a geometry function obtained with GammaLib for viewing period 2.0 (red solid) to the geometry function MPE-DRG-35128 obtained with COMPASS (blue dashed). The left panel shows the distribution of $\chi$ values, obtained by summing overall $\psi$ and $\bar{\phi}$, the middle panel shows the distribution of $\psi$ values, obtained by summing over all $\chi$ and $\bar{\phi}$, and the right panel shows the distribution of $\bar{\phi}$ values, obtained by summing overall $\chi$ and $\psi$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of projections of an IAQ obtained with GammaLib for the energy band 1–3 MeV (red solid) to projections for ROL-IAQ-755 obtained with COMPASS (blue dashed). The left panel shows the distribution of $\phi_{geo}$ values, obtained by summing overall $\bar{\phi}$, the middle panel shows the distribution of $\bar{\phi}$ values, obtained by summing overall $\phi_{geo}$, and the right panel shows angular resolution measure $\bar{\phi} - \phi_{geo}$.</figcaption></figure>

### 3.5. Background modelling

A considerable effort was undertaken by the COMPTEL collaboration to understand and model the instrumental background in the three-dimensional COMPTEL data space (e.g. Bloemen et al. 1994; Knödlseder 1994; van Dijk 1996; Weidenspointner et al. 2001). Satisfactory, although not perfect, results were obtained using the so-called SRCLIX method that was developed by Bloemen et al. (1994) and that has undergone several evolutions (van Dijk 1996). The SRCLIX method iteratively applies the BGDLIX algorithm, of which we implemented two variants in GammaLib. For reference, we also implemented the PHINOR algorithm (described below), and all background modelling methods are invoked in GammaLib via the GCOMObservation::compute_drb method. The iterative SRCLIX method is implemented as a ctools script (cf. Sect. 3.6).

#### 3.5.1. PHINOR

The motivation for the PHINOR algorithm is the observation that the ratio of DRE and DRG multiplied by the solid angle $\Omega(\chi, \psi)$ of the event cube bins is to first order independent of $\chi$ and $\psi$. This leads to a background model given by

$$ \mathrm{DRB}_{\mathrm{phinor}}(\chi, \psi, \bar{\phi}) = \mathrm{DRG}(\chi, \psi, \bar{\phi}) \Omega(\chi, \psi) \times \frac{\sum_{\chi', \psi'} \mathrm{DRE}(\chi', \psi', \bar{\phi}) - \mathrm{DRM}(\chi', \psi', \bar{\phi})}{\sum_{\chi', \psi'} \mathrm{DRG}(\chi', \psi', \bar{\phi}) \Omega(\chi', \psi')} \quad (10) $$

where

$$ \mathrm{DRM}(\chi', \psi', \bar{\phi}) = \int_{\alpha, \delta, E_{\gamma}} I(\alpha, \delta, E_{\gamma}) R(\chi', \psi', \bar{\phi}|\alpha, \delta, E_{\gamma}) d\alpha d\delta dE_{\gamma} \quad (11) $$

is the convolution of a model of celestial sources $I(\alpha, \delta, E_{\gamma})$ with the instrument response function used to subtract any source contribution from the data before normalisation. The evaluation of DRM($\chi', \psi', \bar{\phi}$) is implemented by the method GCOMObservation::drm.

Despite its simplicity, the PHINOR algorithm yields generally empty residual maps above 10 MeV. However, for lower energies significant large-scale residuals are frequently observed (van Dijk 1996).

#### 3.5.2. BGDLIXA

The BGDLIXA algorithm applies a correction to the PHINOR background model by adjusting $\bar{\phi}$ templates TPL($\chi, \psi, \bar{\phi}$) over a limited interval of $\bar{\phi}$ values to a normalisation function $n(\chi, \psi, \bar{\phi})$ that reflects the background event distribution DRE($\chi, \psi, \bar{\phi}$) – DRM($\chi, \psi, \bar{\phi}$) smoothed in $\chi$ and $\psi$ using a running average (see below). The BGDLIXA algorithm evaluates the expression

$$ \mathrm{DRB}'_{\mathrm{bgdlixa}}(\chi, \psi, \bar{\phi}) = \mathrm{TPL}(\chi, \psi, \bar{\phi}) \frac{\sum_{\bar{\phi}' \in \{R_{\bar{\phi}}\}} n(\chi, \psi, \bar{\phi}')}{\sum_{\bar{\phi}' \in \{R_{\bar{\phi}}\}} \mathrm{TPL}(\chi, \psi, \bar{\phi}')}, \quad (12) $$

where the summation is performed over the range $R_{\bar{\phi}} = [\bar{\phi} - \Delta\bar{\phi}_{\mathrm{incl}}, \bar{\phi} - \Delta\bar{\phi}_{\mathrm{excl}}] \cup [\bar{\phi} + \Delta\bar{\phi}_{\mathrm{excl}}, \bar{\phi} + \Delta\bar{\phi}_{\mathrm{incl}}]$ with
$\Delta\bar{\varphi}_{\text{incl}} = \Delta\bar{\varphi}|(N_{\text{incl}} - 1)/2]$, $\Delta\bar{\varphi}_{\text{excl}} = \Delta\bar{\varphi}|(N_{\text{excl}} - 1)/2]$, $N_{\text{incl}}$ and $N_{\text{excl}}$ are either odd integers or zero, and $\lfloor \cdot / \cdot \rfloor$ designates the integer division operator. $\Delta\bar{\varphi}$ is the bin size in $\bar{\varphi}$ of the three-dimensional event cube, which typically is 2°. While in early days $N_{\text{excl}} > 0$ was used, later analyses set $N_{\text{excl}} = 0$, which simplifies the summation range to $R_{\bar{\varphi}} = [\bar{\varphi} - \Delta\bar{\varphi}_{\text{incl}}, \bar{\varphi} + \Delta\bar{\varphi}_{\text{incl}}]$. Past COMPTEL analyses generally used $N_{\text{incl}} = 13$ resulting in $R_{\bar{\varphi}} = [\bar{\varphi} - 12^\circ, \bar{\varphi} + 12^\circ]$. As we will see later, $N_{\text{incl}}$ is the most critical parameter of the BGDLIXA algorithm, and specifies the number of $\bar{\varphi}$ layers of the three-dimensional event cube over which the templates are adjusted. If $N_{\text{incl}}$ is large the fraction in Eq. (12) varies little with $\bar{\varphi}$, leading to a background model that essentially follows the templates. For small values of $N_{\text{incl}}$ the fraction in Eq. (12) accommodates for differences between the templates and the normalisation function, leading to a background model that follows more closely the background event distribution at the expense of including also some of the source events.

The $\bar{\varphi}$ templates are computed using

$$ \text{TPL}(\chi, \psi, \bar{\varphi}) = \text{DRB}_{\text{phinor}}(\chi, \psi, \bar{\varphi}) \times \frac{\sum_{\chi' \in \{R_\chi\}} \sum_{\psi' \in \{R_\psi\}} \sum_{\bar{\varphi}' \in \{\bar{\varphi}'\}} \text{DRE}(\chi', \psi', \bar{\varphi}') - \text{DRM}(\chi', \psi', \bar{\varphi}')}{\sum_{\chi' \in \{R_\chi\}} \sum_{\psi' \in \{R_\psi\}} \sum_{\bar{\varphi}' \in \{\bar{\varphi}'\}} \text{DRB}_{\text{phinor}}(\chi', \psi', \bar{\varphi}')}, \quad (13) $$

which are an adjustment of the PHINOR background model to match the $\bar{\varphi}$-integrated $\chi, \psi$ distribution of the data after subtracting the contributions of celestial sources. To reduce the impact of statistical fluctuations of the data on the model, the match is performed by a running-average summation over a small range in $\chi$ and $\psi$, with $R_\chi = [\chi - \Delta\chi_{\text{runav}}, \chi + \Delta\chi_{\text{runav}}]$ and $R_\psi = [\psi - \Delta\psi_{\text{runav}}, \psi + \Delta\psi_{\text{runav}}]$ where $\Delta\chi_{\text{runav}} = \Delta\chi N_{\text{runav}}$ and $\Delta\psi_{\text{runav}} = \Delta\psi N_{\text{runav}}$. Here, $2N_{\text{runav}} + 1$ specifies the number of $\chi$ and $\psi$ bins over which the running-average summation is performed, with $\Delta\chi$ and $\Delta\psi$ being the bin size in $\chi$ and $\psi$ of the three-dimensional event cube, and $N_{\text{runav}}$ is an integer number. Usually, $\Delta\chi = \Delta\psi = 1^\circ$, and past analyses employed $N_{\text{runav}} = 3$ resulting in a running averaging of $\pm 3^\circ$ around each $\chi, \psi$ pixel (van Dijk 1996).

The adjustment of the $\bar{\varphi}$ templates is done using a normalisation function computed using

$$ n(\chi, \psi, \bar{\varphi}') = \text{DRG}(\chi, \psi, \bar{\varphi}') \Omega(\chi, \psi) \times \frac{\sum_{\chi' \in \{A_\chi\}} \sum_{\psi' \in \{A_\psi\}} \text{DRE}(\chi', \psi', \bar{\varphi}') - \text{DRM}(\chi', \psi', \bar{\varphi}')}{\sum_{\chi' \in \{A_\chi\}} \sum_{\psi' \in \{A_\psi\}} \text{DRG}(\chi', \psi', \bar{\varphi}') \Omega(\chi', \psi')}, \quad (14) $$

which is an analogue of Eq. (10) that is limited to a small range of $\chi, \psi$ around the pixel of interest. Specifically, $A_\chi = [\chi - \Delta\chi_{\text{avgr}}, \chi + \Delta\chi_{\text{avgr}}]$ and $A_\psi = [\psi - \Delta\psi_{\text{avgr}}, \psi + \Delta\psi_{\text{avgr}}]$ where $\Delta\chi_{\text{avgr}} = \Delta\chi |(N_{\text{avgr}} - 1)/2]$ and $\Delta\psi_{\text{avgr}} = \Delta\psi |(N_{\text{avgr}} - 1)/2]$. Here, $N_{\text{avgr}}$ is an odd integer that specifies the number of $\chi$ and $\psi$ bins over which the adjustment is performed, with $\Delta\chi$ and $\Delta\psi$ being the bin size in $\chi$ and $\psi$ of the three-dimensional event cube. Past analyses usually employed $N_{\text{avgr}} = 3$ resulting in an averaging of $A_\chi = [\chi - 1^\circ, \chi + 1^\circ]$ and $A_\psi = [\psi - 1^\circ, \psi + 1^\circ]$ (van Dijk 1996).

We want to point out that the COMPASS analysis system implemented

$$ n(\chi, \psi, \bar{\varphi}') = \text{DRG}(\chi, \psi, \bar{\varphi}') \Omega(\chi, \psi) \times \sum_{\chi' \in \{A_\chi\}} \sum_{\psi' \in \{A_\psi\}} \frac{\text{DRE}(\chi', \psi', \bar{\varphi}')}{\text{DRG}(\chi', \psi', \bar{\varphi}') \Omega(\chi', \psi')} \quad (15) $$

instead of Eq. (14), which can lead to unreasonably large normalisations at the edge of the three-dimensional event cube where DRG values are small. Using Eq. (14) avoids such problems, and in fact leads to a simplified background modelling algorithm that we implemented in GammaLib under the acronym BGDLIXE (see below). Furthermore, we added a global $\bar{\varphi}$ normalisation step,

$$ \text{DRB}_{\text{bgdlixa}}(\chi, \psi, \bar{\varphi}) = \text{DRB}'_{\text{bgdlixa}}(\chi, \psi, \bar{\varphi}) \times \frac{\sum_{\chi', \psi'} \text{DRE}(\chi', \psi', \bar{\varphi}) - \text{DRM}(\chi', \psi', \bar{\varphi})}{\sum_{\chi', \psi'} \text{DRB}'_{\text{bgdlixa}}(\chi', \psi', \bar{\varphi})}, \quad (16) $$

at the end of the computation so that the model is normalised to the number of background events in each $\bar{\varphi}$ layer.

### 3.5.3. BGDLIXE

Recognising that Eqs. (13) and (14) both do the same thing (that is, they normalise the solid-angle-weighted geometry function to the measured event distribution), the BGDLIX background model can be simplified to

$$ \text{DRB}'_{\text{bgdlixa}}(\chi, \psi, \bar{\varphi}) = \text{DRB}_{\text{phinor}}(\chi, \psi, \bar{\varphi}) \times \frac{\sum_{\chi' \in \{A_\chi\}} \sum_{\psi' \in \{A_\psi\}} \sum_{\bar{\varphi}' \in \{\bar{\varphi}'\}} \text{DRE}(\chi', \psi', \bar{\varphi}') - \text{DRM}(\chi', \psi', \bar{\varphi}')}{\sum_{\chi' \in \{A_\chi\}} \sum_{\psi' \in \{A_\psi\}} \sum_{\bar{\varphi}' \in \{\bar{\varphi}'\}} \text{DRB}_{\text{phinor}}(\chi', \psi', \bar{\varphi}')}, \quad (17) $$

which is a local fit of the PHINOR background model to the three-dimensional event distribution after subtracting any celestial sources. We note that $N_{\text{runav}}$ is no longer a parameter of the model, and the number of pixels in $\chi$ and $\psi$ that are used for the local fit is determined by $N_{\text{avgr}}$. We also included a $\bar{\varphi}$ normalisation at the end of the computation equivalent to Eq. (16).

Figure 6 compares the event cube projections to the projections for background models obtained using PHINOR, BGDLIXA and BGDLIXE using the parameters $N_{\text{runav}} = 3$, $N_{\text{avgr}} = 3$, $N_{\text{incl}} = 13$ and $N_{\text{excl}} = 0$ proposed by van Dijk (1996). While the PHINOR model provides a first-order description of the event distribution, it shows significant differences in the details. The BGDLIXA and BGDLIXE models follow the event distribution very closely, and both models are in fact indistinguishable from each other, demonstrating that their results are equivalent for the parameter values chosen.

### 3.6. ctools scripts

To support the science analysis of COMPTEL data, we extended the ctools package with a number of dedicated Python scripts, providing basic building blocks that each perform well-defined science data analysis tasks. This includes the generation of a database based on the data available in the HEASARC archive, the selection of COMPTEL viewing periods from this database, event binning and response computation, data combination, background modelling, model fitting, generation of test statistic (TS) maps, residual inspection, observation simulations, and source detection, as well as generation of pulsar pulse profiles. These scripts complement the already existing generic science analysis tools in ctools that may be used in combination with the new scripts. Table 2 summarises the COMPTEL-specific scripts that we added to ctools.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Comparison of PHINOR (blue dotted), BGDLIXA (dark green dashed), and BGDLIXE (light green dashed-dotted) background model for viewing period 1.0 and the energy band 1–3 MeV to the DRE event distribution (red solid). The left panel shows the distribution of $\chi$ values, obtained by summing overall $\psi$ and $\bar{\psi}$, the middle panel shows the distribution of $\psi$ values, obtained by summing over all $\chi$ and $\bar{\psi}$, and the right panel shows the distribution of $\bar{\psi}$ values, obtained by summing overall $\chi$ and $\psi$.</figcaption></figure>

<table><caption>**Table 2.** COMPTEL-specific ctools scripts.</caption><thead><tr><th>Script</th><th>Task</th></tr></thead><tbody><tr><td>comgendb</td><td>Generate database</td></tr><tr><td>comobsselect</td><td>Select observations</td></tr><tr><td>comobsbin</td><td>Bin COMPTEL observations</td></tr><tr><td>comobsback</td><td>Generate background model</td></tr><tr><td>comobsadd</td><td>Combine observations</td></tr><tr><td>comobsmodel</td><td>Generate XML model for analysis</td></tr><tr><td>comlixfit</td><td>Fit model to data using SRCLIX method</td></tr><tr><td>comlixmap</td><td>Create TS map using SRCLIX method</td></tr><tr><td>comsrcdetect</td><td>Detect sources in TS map</td></tr><tr><td>comobsres</td><td>Generate residuals</td></tr><tr><td>comobssim</td><td>Simulate COMPTEL observations</td></tr><tr><td>compulbin</td><td>Generate pulsar pulse profile</td></tr></tbody></table>

Before starting a COMPTEL data analysis, a database needs to be generated from the data that are available in the HEASARC archive. This analysis step, which only needs to be performed once on a given computer, is performed by the `comgendb` script.

Once the database is generated, a typical COMPTEL analysis starts with executing `comobsselect` to select the relevant viewing periods for a given source or source region and observing time interval from the database. Results are provided as an observation definition file in XML format, which contains all the relevant information for the subsequent analysis steps. Following selection, the data are binned into three-dimensional event cubes and the corresponding DRX and DRG cubes are generated using `comobsbin`. The binning can be done for an arbitrary number of energy bands, enabling a joint spectral-spatial analysis that was not supported by the original COMPASS software. `comobsbin` also computes an initial background model DRB using the PHINOR algorithm (see Sect. 4.1), as well as the IAQ response matrices for the relevant energy bands. All output files are stored as FITS files in a data store, avoiding the recompulation of identical data files in subsequent analyses. Alternative background models using the algorithms described in Sect. 3.5 can be generated using `comobsback`.

Individual viewing periods can be combined into single event cubes for each energy band {$E$} using the `comobsadd` script, leading to a considerable speed-up of the data analysis. The combination of viewing periods is however not required, and alternatively the data can be analysed using a joint maximum-likelihood analysis of the selected viewing periods, a possibility

that was not supported by the original COMPASS software. Combination of the viewing periods $i$ is done using

$$ \text{DRE}(\chi, \psi, \bar{\psi}|\{E\}) = \sum_i \text{DRE}_i(\chi, \psi, \bar{\psi}|\{E\}) \quad (18) $$

$$ \text{DRB}(\chi, \psi, \bar{\psi}|\{E\}) = \sum_i \text{DRB}_i(\chi, \psi, \bar{\psi}|\{E\}) \quad (19) $$

$$ \text{DRX}(\alpha, \delta) = \sum_i \max_{\alpha', \delta'} \text{DRX}_i(\alpha', \delta') \quad (20) $$

$$ \text{DRG}(\chi, \psi, \bar{\psi}) = \frac{1}{T} \sum_i T_i \text{DRG}_i(\chi, \psi, \bar{\psi}), \quad (21) $$

where $T_i$ is the exposure time of viewing period $i$ and $T = \sum_i T_i$ is the total exposure time of the combined data. We note that Eq. (20) leads to an exposure map that is independent of $\alpha, \delta$ and that is given by the maximum of the DRX for each viewing period, an approximation that is justified by the fact that the DRX is rather flat (see Fig. 3) and that the zenith angle variation in the response computation (cf. Eq. (3)) is dominated by the geometry function DRG (see also Knödlseder 1994).

In order to perform a maximum-likelihood analysis of the data, a model definition file in XML format needs to be generated. We automatised this task with the `comobsmodel` script, which in particular generates an adequate model definition for fitting the background model DRB to the data. COMPTEL background model fitting is done using the `DRBPhibarBins` model type, which introduces one scaling factor for each $\bar{\psi}$ layer for all viewing periods and energy bands. In addition, `comobsmodel` supports adding of a point source and diffuse model components to the model definition XML file.

The main script for maximum-likelihood model fitting is `comlixfit`, which implements the iterative SRCLIX algorithm. The algorithm starts from an input model definition XML file and computes a DRM model of the celestial sources using Eq. (11), which is then used to compute an initial background model using one of the PHINOR, BGDLIXA, or BGDLIXE algorithms. The script then uses the `ctlike` tool to fit the source and background models to the binned data. The celestial source model that results from the fit is then used to update the DRM model and to regenerate a background model using the selected algorithm. A new `ctlike` fit is then performed using the updated model, and the procedure is repeated until the maximum log-likelihood change between subsequent iterations becomes negligible, typically less than 0.05. Usually, the SRCLIX algorithm converges after a few iterations.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Test statistic maps of viewing period 1.0 (during which the Crab pulsar and pulsar wind nebula was observed) for the combined analysis of the four standard COMPTEL energy bands as a function of the number of SRCLIX iterations, indicated in the upper-right corner of each panel. The BGDLIXE algorithm with $N_{avgr} = 5$, $N_{incl} = 15$, and $N_{excl} = 0$ was employed. The lower-right panel shows the TS map obtained using the comlixmap script. Colour maps are shown in units of $\sqrt{\text{TS}}$. Negative values correspond to negative fluxes of the test source.</figcaption></figure>

Figure 7 illustrates the algorithm by showing TS maps for subsequent SRCLIX iterations for viewing period 1.0 that were obtained using the cttsmap tool. The TS is defined as

$$TS = 2 \ln L(M_s + M_b) - 2 \ln L(M_b) \qquad (22)$$

(Mattox et al. 1996), where $\ln L(M_s + M_b)$ is the log-likelihood value obtained when fitting the source and the background models together to the data, and $\ln L(M_b)$ is the log-likelihood value obtained when fitting only the background model to the data. Under the hypothesis that the background model $M_b$ provides a satisfactory fit of the data, TS follows a $\chi_n^2$ distribution with n degrees of freedom, where n is the number of free parameters in the source model component. Therefore,

$$p = \int_{\text{TS}}^{+\infty} \chi_n^2(x) \, dx \qquad (23)$$

gives the chance probability (p-value) that the log-likelihood improves by TS/2 when adding the source model $M_s$ due to statistical fluctuations only (Cash 1979).

We note that viewing period 1.0 is an observation of the Crab pulsar and pulsar wind nebula, which is by far the brightest source of gamma rays in the COMPTEL energy range. For the SRCLIX analysis, the data of the four standard energy bands were analysed jointly, and the Crab was modelled as a point source with power-law spectrum. The source location, prefactor and spectral index were free parameters in the fit. After each iteration of the SRCLIX algorithm, a TS map<fnref n="7" /> was generated by replacing the Crab in the fitted model by a test source with fixed power-law spectral index of $\Gamma = 2$. The source visible at

the centre of the TS maps is the Crab, which is already significantly detected after the first SRCLIX iteration. Subsequent iterations slightly increase the source significance, but overall the TS maps change little. The SRCLIX algorithm converged after six iterations.

The TS maps show a halo of marginal significance around the Crab, which can be explained by the fact that the maps were obtained using a background model that assumed the presence of a point source at the location of the Crab. The event cone of a test source placed a few degrees away from the Crab will pick up some of the excess counts left by the background model, which explains the halo in maps.

An alternative way to generate TS maps is provided by the comlixmap script, which applies the SRCLIX algorithm to each test source position, and which is the algorithm that was implemented by the COMPASS task SRCLIX. Here each pixel in the TS map corresponds to a different background model, and when moving away from the source location, no excess counts will remain in the data. This reduces the halo around the sources, as illustrated in the last panel of Fig. 7 that shows the TS map that was obtained using comlixmap for viewing period 1.0. At the same time, negative residuals are amplified, which is explained by the fact that the events of the Crab for test source positions offset from the source will be included as a smoothed event cone in the background model. We note, however, that this is only a feature of the TS maps, as ultimately, an adequate model of celestial sources that describes the COMPTEL event distribution should be derived from the data. Using that model as DRM in the BGDLIX algorithm will exclude any source events from the smoothing algorithm, and hence provides a reliable background model. For this approach the comsrcdetect script can be used, which extracts significant sources from TS maps and adds them to a model definition XML file. A subsequent run

<footnote n="7">⁷ The TS map comprises 50 × 50 pixels of size 1° × 1° situated around the Crab position.</footnote>
of the comlixmap script will then show whether any additional celestial sources remain in the data, building up iteratively an adequate model of celestial sources.

Following model fitting an inspection of the fit residuals is crucial. The comobsres script enables such an inspection by projecting the residual between the event and model cubes onto the sky by summing their content along the event cone using

$$N(\alpha, \delta) = \sum_{\chi, \psi, \bar{\varphi}} \begin{cases} \text{DRE}(\chi, \psi, \bar{\varphi}), & \text{if } \text{ARM}_{\min} \le \bar{\varphi} - \varphi_{\text{geo}} \le \text{ARM}_{\max} \\ 0, & \text{otherwise} \end{cases} \qquad (24)$$

and

$$M(\alpha, \delta) = \sum_{\chi, \psi, \bar{\varphi}} \begin{cases} \text{DRM}(\chi, \psi, \bar{\varphi}), & \text{if } \text{ARM}_{\min} \le \bar{\varphi} - \varphi_{\text{geo}} \le \text{ARM}_{\max}, \\ 0, & \text{otherwise} \end{cases} \qquad (25)$$

where $\varphi_{\text{geo}}$ is the angular distance between a sky position $(\alpha, \delta)$ and the Compton scatter direction $(\chi, \psi)$, and $[\text{ARM}_{\min}, \text{ARM}_{\max}]$ defines a selection window for the so-called angular resolution measure that is typically taken to be $[-3^\circ, 3^\circ]$. We note that DRM designates here the model cube that comprises both the source and background components. By default, comobsres uses

$$R(\alpha, \delta) = \operatorname{sgn}(N(\alpha, \delta) - M(\alpha, \delta)) \times \sqrt{2 \left( N(\alpha, \delta) \ln \frac{N(\alpha, \delta)}{M(\alpha, \delta)} + M(\alpha, \delta) - N(\alpha, \delta) \right)} \qquad (26)$$

to compute the significance of the residuals $R(\alpha, \delta)$ in Gaussian $\sigma$, where the sign term indicates whether the measured number of counts is larger or smaller than the number of counts predicted by the model. Some special cases need to be treated separately. Namely, if $N(\alpha, \delta) = 0$ the residual significance is

$$R(\alpha, \delta) = \operatorname{sgn}(N(\alpha, \delta) - M(\alpha, \delta)) \sqrt{2M(\alpha, \delta)}, \qquad (27)$$

while if $M(\alpha, \delta) = 0$ the significance cannot be computed and we set $R(\alpha, \delta) = 0$.

The comobssim script enables the simulation of DRE event cubes by sampling the events according to a Poisson distribution using the expectation given by a model. Specifically, simulated events for a given celestial source model can be added by comobssim to existing observations, allowing the study of celestial sources with known properties in real data, a possibility that we use extensively in the next section.

Finally, compulbin will generate pulsar phase profiles by applying the algorithms described in Appendix D to individual events. Only events that satisfy the ARM selection according to Eq. (24) will be retained in the phase profiles, with typical values for the ARM window being $[-3.5^\circ, 3.5^\circ]$.

## 4. Science validation

Having verified that GammaLib and ctools reproduce data products that are identical to the ones that were generated with the COMPASS system, we now verify that the use of our software for an analysis of the data provided by HEASARC reproduces COMPTEL science results that were published in the literature. If not stated otherwise, for the analyses that follow we apply the event selection parameters specified in Table 1, a value of $\zeta_{\min} = 5^\circ$, and we exclude D2 modules for which there were faulty photomultipliers.

### 4.1. Background model validation

An important step prior to any science analysis is the validation of the background model. An accurate background model predicts the background event distribution within statistical fluctuations, and allows for a reliable and accurate determination of the contribution of celestial gamma-ray sources to the measured events. As is obvious from Fig. 6, the PHINOR model certainly does not satisfy the first criterion; hence, we no longer consider this model here. On the other hand, the BGDLIXA and BGDLIXE models look promising, and were generally used in the past for science analysis of COMPTEL data. Since the BGDLIXA and BGDLIXE algorithms are equivalent as long as $N_{\text{runav}} \le N_{\text{avgr}}$, we limit our study here to the simpler BGDLIXE algorithm. Specifically, we investigate which values of $N_{\text{avgr}}$ and $N_{\text{incl}}$ provide reliable background models without introducing a significant bias in the reconstruction of celestial gamma-ray source characteristics. In agreement with previous studies of the algorithms (cf. van Dijk 1996) we always set $N_{\text{excl}} = 0$.

#### 4.1.1. Residual maps

We started with modelling the background for all COMPTEL viewing periods and the four standard COMPTEL energy bands using comobsback and the BGDLIXE algorithm under variation of the parameters $N_{\text{avgr}}$ and $N_{\text{incl}}$. For each viewing period and energy band we created residual maps and histograms using comobsres that we inspected visually. We note that due to the dominance of the instrumental background in COMPTEL data, we do not expect to see celestial gamma-ray emission in the residual maps of individual viewing periods, with the exception of emission from the Crab pulsar and pulsar wind nebula, which is the strongest gamma-ray source at MeV energies. We find a general trend of stronger residuals at lower energies, with the largest residuals observed in the 1–3 MeV band. Notably, residuals are stronger for viewing periods for which the Earth horizon selection Eq. (2) introduces a strong $\bar{\varphi}$ dependence in the $\chi, \psi$ distribution of the events. The amplitude of the residuals changes strongly under variation of $N_{\text{incl}}$, with larger residuals for larger values of $N_{\text{incl}}$, while variations of $N_{\text{avgr}}$ impact the residuals only moderately.

For illustration we show residual maps obtained using the BGDLIXE algorithm for different values of $N_{\text{avgr}}$ and $N_{\text{incl}}$ and for three representative viewing periods and the 1–3 MeV energy band in Figs. 8 and 9. Viewing period 1.0 is an observation of the Crab pulsar and pulsar wind nebula, which is clearly visible in the residuals maps. Viewing period 223.0 is an observation of 1E 1740.7–2942, a low-mass X-ray binary that is also known as the ‘Great Annihilator’ and that is situated near the Galactic centre. The residuals in this viewing period are relatively modest. Viewing period 518.5 is an observation of the BL Lacertae object S5 0716+714, which is among the viewing periods with the worst residuals, featuring large zones of significant excess counts and negative depressions. The amplitude of the residuals clearly decreases with decreasing value of $N_{\text{incl}}$, while at the same time smaller value of $N_{\text{incl}}$ also reduce the signal from the Crab. As we will show later, some of this signal can be recovered using the iterative SRCLIX algorithm, yet small values of $N_{\text{incl}}$ tend to lead to an underestimation of source fluxes. Therefore, the choice of $N_{\text{incl}}$ is necessarily a trade-off between the amplitude of background residuals and the suppression of source flux. On the other hand, the amplitude of the residuals changes only moderately with $N_{\text{avgr}}$, with a slight increase of the residual amplitudes for increasing values of $N_{\text{avgr}}$.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Residual maps for three viewing periods and the 1–3 MeV energy band with background modelled using the BGDLIXE algorithm for $N_{\text{avgr}} = 5$ and different values of $N_{\text{incl}}$. The colour scale is limited to the range –5σ (blue) to +5σ (red). Viewing period 1.0 was selected because it contains the Crab, viewing period 223.0 because it has a rather flat residual, and viewing period 518.5 because it is among the viewing periods with the worst residuals.</figcaption></figure>

To systematically quantify the residuals for a given choice of BGDLIXE parameters we computed for all viewing periods and the four standard COMPTEL energy bands the mean and random mean scatter (rms) of the residual maps. The results for the BGDLIXE algorithm as a function of $N_{\text{incl}}$ for $N_{\text{avgr}} = 5$ and as a function of $N_{\text{avgr}}$ for $N_{\text{incl}} = 13$ are shown in Fig. 10. The violins represent the density distribution of the mean and rms of the residual maps for all viewing periods. Viewing periods including the Crab were excluded to avoid any bias due to the presence of a strong source.

The plots confirm that the largest spread in the mean and rms are observed for lower energies, with a particularly large spread for the 1–3 MeV energy band. A large spread indicates that for some viewing periods the background model results in important residuals, while for other viewing periods the algorithm performs rather well, as illustrated in Figs. 8 and 9. Reducing

$N_{\text{incl}}$ considerably reduces the spread in the mean and rms values for all energy bands, yet at some point the rms drops below the expected value of 1, indicating that the background model partially follows the statistical fluctuations of the data. This is also the regime where the source fluxes start to get underestimated.

This overfitting can be slightly reduced by increasing $N_{\text{avgr}}$ so that more events get included in the $\chi, \psi$ averaging, as indicated in the lower row of Fig. 10. On the other hand, increasing $N_{\text{avgr}}$ leads to a slight increase of the mean and rms distributions; hence, the selected value of $N_{\text{avgr}}$ should also not be too large.

### 4.1.2. Flux reconstruction

As the next step we studied the impact of the BGDLIXE parameters on the fitted values of celestial source parameters, such as source position, extent, flux and spectral index. We do this by
<figure><img src="image_10.png" /><figcaption>Fig. 9. Same as Fig. 8 but for $N_{\text{incl}} = 13$ and different values of $N_{\text{avgr}}$.</figcaption></figure>

using `comobssim` to add a simulated source at an offset angle of 20° with respect to the pointing axis to the data of each viewing period for the four standard energy bands. In that way, our study relies essentially on the observed event distribution, and hence is representative for a real analysis situation, while the characteristics of the celestial source are known. As simulated source model we use a spatially extended disk component with radius of 3° combined with a spectral power-law component with an energy flux of $1.068 \times 10^{-8}$ erg cm$^{-2}$ s$^{-1}$ within 0.75–30 MeV and a spectral index of 2.1, which roughly corresponds to the spectral parameters that are observed for the Crab. The simulated data of each viewing period were then fit jointly for the four energy bands using `comlixfit`, determining the maximum likelihood right ascension, declination, disk radius, energy flux, and spectral index of the source. Initial values for the iterative fitting procedure were offset from the simulated values since in a real data analysis the true source parameters are generally not known

in advance. Viewing periods including the Crab were excluded from the analysis to avoid any interference with a known strong source of gamma rays. We furthermore assume that no other source of gamma rays is significantly detected in any of the remaining individual viewing periods, considering these viewing periods as empty fields for the purpose of this analysis.

For each fitted source parameter, $i$, we determine the pull

$$ \text{Pull}(p_i) = \frac{p_i - v_i}{\sigma_i}, \qquad (28) $$

where $p_i$ is the fitted value, $v_i$ is the simulated value, and $\sigma_i$ is the statistical parameter uncertainty as determined by `comlixfit`. In the absence of systematic uncertainties, and under the assumption that the statistical uncertainties are following a Gaussian distribution, Pull($p_i$) follows a Gaussian distribution with a mean of zero and a standard deviation of $\sigma = 1$.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Violin plots of the mean and rms, in units of $\sigma$, of the residual maps with background modelled using the BGDLIXE algorithm for all viewing periods except the ones that include the Crab. The mean and rms are shown as a function of $N_{\text{incl}}$ for $N_{\text{avgr}} = 5$ in the *upper row* and as a function of $N_{\text{avgr}}$ for $N_{\text{incl}} = 13$ in the *lower row*. The groups of four violins correspond to the four standard energy bands; *from left to right*: 0.75–1 MeV (blue), 1–3 MeV (orange), 3–10 MeV (green), and 10–30 MeV (red). Horizontal black bars show the maximum, the median, and the minimum values, respectively.</figcaption></figure>

Figure 11 summarises the results of the analysis, showing the pull distributions for all considered viewing periods as violin plots, with violins for right ascension, declination, disk radius, energy flux, and spectral index. The expectations for purely statistical parameter fluctuations are indicated by semi-transparent violins with grey borders, the horizontal black bars indicate the maximum, median and minimum pull of the distributions. The upper row shows results as a function of $N_{\text{incl}}$ for $N_{\text{avgr}} = 5$, the lower row shows results as a function of $N_{\text{avgr}}$ for $N_{\text{incl}} = 13$.

Figure 11 indicates that large values of $N_{\text{incl}}$ lead to pull distributions that are broader than expected from statistical fluctuations only, in particular for the energy flux, but also for the spectral index and to a lesser extent for the other source parameters. The broadening is due to background residuals in some of the viewing periods for large $N_{\text{incl}}$, as illustrated in Fig. 8, which impact the reconstructed source parameters. Reducing $N_{\text{incl}}$ brings the pull distributions more in line with the expectations, yet for $N_{\text{incl}} \le 13$ the median pull of the fitted energy flux drops below zero, indicating a systematic bias towards too low fluxes. Specifically, for $N_{\text{incl}} = 5$, where background residuals are very small (cf. Fig. 8), the flux reconstruction is strongly biassed, resulting is a significant underestimation of gamma-ray fluxes. The optimum value for $N_{\text{incl}}$ is hence a trade-off between reduction of background residuals and biassing the flux determination.

Flux biassing is also observed with increasing value of $N_{\text{avgr}}$, with a minimum bias that occurs for $N_{\text{avgr}} = 5$. Using hence $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 15$ for the BGDLIXE background modelling leads to results that are basically free from any systematic

bias, although some significant background residuals may persist for this choice of values. As illustrated in Fig. 11, these background residuals translate into an additional uncertainty beyond the statistical fluctuations only. In the present case, the standard deviation of the energy flux pull distribution for $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 15$ is about twice as large as expected from statistical fluctuations only. In other words, when analysing individual viewing periods using BGDLIXE parameters $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 15$, the uncertainties in the energy flux related to the background modelling roughly doubles the uncertainties due to statistical fluctuations only.

In the following we use $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 15$ for the analysis in our paper, and we generally recommend to use these parameters for COMPTEL data analysis with GammaLib and ctools. We note that these values apply for a binning of 1° in $\chi$ and $\psi$ and 2° in $\tilde{\phi}$, and that for a different binning the parameters need to be adjusted accordingly. Specifically, we used an equivalent value of $N_{\text{incl}} = 29$ for analyses in our paper for which a binning of 1° is applied in $\tilde{\phi}$.

## 4.2. Gamma-ray emission from the Crab

We began the science validation of our software with a spectral analysis of the gamma-ray emission from the Crab pulsar and pulsar wind nebula, which is the brightest source of gamma rays at MeV energies. The MeV flux is dominated by emission from the pulsar wind nebula, yet using pulsar ephemerides derived from radio observations the emission from the pulsar is also clearly detectable over the entire COMPTEL energy range. The
<figure><img src="image_12.png" /><figcaption>Fig. 11. Violin plots of the fitted parameter pull distributions for a simulated source at a 20° off-axis angle for all viewing periods except those that include the Crab. The upper row shows results as a function of $N_{\text{incl}}$ for $N_{\text{avgr}} = 5$; the lower row shows results as a function of $N_{\text{avgr}}$ for $N_{\text{incl}} = 13$. For all simulations the source was simulated as a spatially extended disk with a 3° radius and with a power-law spectrum of index $\Gamma = 2.1$. Horizontal black bars show the maximum, the median, and the minimum values, respectively. Semi-transparent violins with grey borders indicate the expected pull distributions from statistical fluctuations only.</figcaption></figure>

emission from the Crab pulsar and pulsar wind nebula was studied extensively by COMPTEL in the past (e.g. Much et al. 1995a,b, 1996; van der Meulen et al. 1998; Kuiper et al. 2001).

### 4.2.1. Total spectrum

We first considered the combined emission from the Crab pulsar and pulsar wind nebula. In their analysis of five years of COMPTEL observations, van der Meulen et al. (1998) analysed the spectrum of the total Crab emission within the energy range 0.78–30 MeV in 30 narrow energy bins. Since this is the only work that quotes total flux values for the Crab (cf. Table 2 of the publication), we used this study as reference.

We analysed the same data that was used by van der Meulen et al. (1998) with GammaLib and ctools, except for viewing period 0 that is not available at HEASARC and viewing period 426.0 for which the EVP file in the HEASARC archive has a corrupted content. We binned the data according to the 30 energy bins defined by van der Meulen et al. (1998) and combined the data for all viewing periods using comobsadd using 80 bins in $\chi$ and $\psi$ that were centred on the position of the Crab pulsar, taken here to be 83.6331° in right ascension and 22.0145° in declination (J2000). Similar to van der Meulen et al. (1998) we used 50 bins in $\bar{\phi}$, bin sizes of 1° in all three data space dimensions, and instrument response functions derived by analytical modelling.

We modelled the Crab using a point source with fixed position, as given above, and using power-law, exponentially cut-off power-law or curved power-law spectral models. We fitted the data jointly for the 30 energy bins using comlixfit

with BGDLIXE background model parameters $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 29$. The best fit was obtained using the exponentially cut-off power law

$$I(E_{\gamma}) = k \left( \frac{E_{\gamma}}{3.5 \text{ MeV}} \right)^{-\Gamma} \exp \left( -\frac{E_{\gamma}}{E_{c}} \right), \quad (29)$$

with $k = (1.81 \pm 0.06) \times 10^{-4} \text{ ph cm}^{-2} \text{ s}^{-1} \text{ MeV}^{-1}$, $\Gamma = 2.00 \pm 0.03$, and $E_c = 39.1 \pm 9.7 \text{ MeV}$. van der Meulen et al. (1998) do not quote the parameters of a fitted spectral model to the total Crab emission data, so we used the GammaLib multi-wavelength interface to adjust the same spectral models to the data of their Table 2, which also favoured the exponentially cut-off power law with best fitting parameters $k = (1.70 \pm 0.05) \times 10^{-4} \text{ ph cm}^{-2} \text{ s}^{-1} \text{ MeV}^{-1}$, $\Gamma = 1.97 \pm 0.02$, and $E_c = 29.8 \pm 5.0 \text{ MeV}$. While our prefactor $k$ is about 6% larger than the one obtained from fitting the spectrum of van der Meulen et al. (1998), the other spectral parameters are equivalent within statistical uncertainties.

We then used ctbutterfly to determine uncertainty bands for the spectral models and cspec to derive flux points for the 30 energy bins. The results of our analysis are compared to those of van der Meulen et al. (1998) in Fig. 12. The uncertainty bands for the data of van der Meulen et al. (1998) were determined using ctbutterfly. Overall the agreement between both analyses is quite good, yet as mentioned earlier, our flux points lie somewhat above the ones determined by van der Meulen et al. (1998). Possibly this discrepancy may be related to correction factors that were applied posterior to COMPASS analyses at the time that were not automatically taken
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Spectral energy distribution of the total emission from the Crab pulsar and pulsar wind nebula as measured by COMPTEL in 30 bins covering the energy band 0.78–30 MeV. Filled red dots correspond to results obtained with ctools, and open blue dots correspond to results obtained by van der Meulen et al. (1998) using COMPASS. The shaded regions correspond to the $1\sigma$ uncertainty bands of the fitted exponentially cut-off power-law models.</figcaption></figure>

into account by the software. These correction factors include an energy-dependent ToF correction factor (cf. Table C.1), an energy-independent deadtime correction factor of 0.965 as well as an energy-independent flux correction factor that was eventually applied to SRCLIX analyses to correct for a flux suppression that arose from the modification of the instrument response function (van Dijk 1996). We recall that GammaLib automatically applies the ToF and deadtime correction factors to the results. Whether or not such correction factors were applied by van der Meulen et al. (1998) is not known, yet they may plausibly explain the 6% discrepancy observed between the analyses.

### 4.2.2. Pulsar and nebula components

We now turn to an analysis that separates the emission from the Crab pulsar and the pulsar wind nebula to validate the implementation of the pulsar phase computations. The most comprehensive analysis of the emission from the Crab pulsar and pulsar wind nebula components using COMPTEL data was performed by Kuiper et al. (2001) using data collected over the nine years mission duration of CGRO. While Kuiper et al. (2001) used data from 33 viewing periods with pointing axis within 30° of the Crab pulsar, we analysed 24 viewing periods that we found with the same selection criteria in the HEASARC database, covering the viewing periods specified in Table 1 of Kuiper et al. (2001) between viewing period 1.0 and viewing period 616.1. As in the analysis above, viewing period 426.0 was excluded from the list since no usable EVP file exists for this observation in the HEASARC database. Similarly to Kuiper et al. (2001) we used ephemerides for the Crab pulsar from the Princeton radio pulsar database, provided in the form of an ASCII file named psrtime.dat<fnref n="8" /> that is part of the reference data of the X-ray Timing Explorer (XTE) module of the HEASoft software (version 6.29)<fnref n="9" />.

We first used `compulbin` with an angular resolution measure of $\pm 3.5^\circ$ to produce pulse profiles for the four standard

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Pulse profiles of the Crab pulsar derived using ctools for the four standard energy bands (red) compared to pulse profiles derived by Kuiper et al. (2001) using the COMPASS software (blue). The latter profiles were scaled to match our results in amplitude, and a phase shift of $-0.03$ was applied to match the profiles in phase.</figcaption></figure>

COMPTEL energy bands, as displayed in Fig. 2 of Kuiper et al. (2001). The results of this analysis are shown in Fig. 13, on which we overlay for comparison the pulse profiles obtained by Kuiper et al. (2001). Since the Kuiper et al. (2001) profiles were obtained for a larger dataset and angular resolution measures that were not specified in their paper, we scaled the profiles so that the minimum and maximum number of events in the profiles matches the numbers that we obtained in our analysis. We also applied a phase shift of $-0.03$ to the pulse profiles of Kuiper et al. (2001) to match them to our profiles. While we do not know the origin of this small discrepancy in the pulse phases, we note that a phase shift of $-0.03$ corresponds to a difference of about 0.5 arcsec in the assumed right ascension of the Crab pulsar. In GammaLib, the pulsar position is taken from the ephemerides file, which in the present case is the radio position in the Princeton database, while Kuiper et al. (2001) do not specify the position that they assumed for the Crab pulsar. We note that the radio position in the Princeton database differs by about 0.5 arcsec from the International Celestial Reference System (ICRS) position provided by the Set of Identifications, Measurements, and Bibliography for Astronomical Data (SIMBAD) service<fnref n="10" />, which could be at the origin of the observed phase shift.

<footnote n="8">⁸ Details are provided at https://heasarc.gsfc.nasa.gov/lheasoft/ftools/fhelp/fasebin.html</footnote>

<footnote n="9">⁹ The HEASoft software can be downloaded from https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/</footnote>

<footnote n="10">¹⁰ See http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=Crab+pulsar&submit=SIMBAD+search</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Spectral energy distributions of the Crab pulsar and pulsar wind nebula components as determined using GammaLib and ctools (red) and by Kuiper et al. (2001) using COMPASS (blue). Results for the Crab pulsar are shown as dots, and results for the Crab pulsar wind nebula are shown as triangles. The figure also shows the $1\sigma$ uncertainty bands of the best fitting spectral models for both components.</figcaption></figure>

As the next step we determined the spectra of the Crab pulsar and pulsar wind nebula to compare them to those given in Table 3 of Kuiper et al. (2001). For this purpose, we binned the events using `comobsbin` for the 12 energy bins specified in that table. The data were binned separately for the Off Pulse and Total Pulse phase intervals as defined in Table 2 of Kuiper et al. (2001), shifted by -0.03 to accommodate for the observed phase shift. Specifically, the Off Pulse interval comprises phases 0.49–0.85 while the Total Pulse interval comprises phases 0.85–1 and 0–0.49. We fitted the data for both intervals jointly using `comlixfit` with the BGDLIXE parameters $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 29$. We used two point-source model components in our fit, one for the Crab nebula that was fitted to the data of both phase intervals, and one for the Crab pulsar that was only fitted to the data of the Total Pulse interval. Consequently, the Crab pulsar component modelled only events that were in excess of the pulsar wind nebula component. Both point-source model components were located at the position of the Crab pulsar, as defined above, and had a spectral model with a free flux value for each of the 12 energy bins. Similar to Kuiper et al. (2001), we used simulated instrument response functions and 50 bins in $\bar{\phi}$ with bin sizes of 1° for our analysis.

The spectra obtained with our analysis are shown in Fig. 14 together with the spectra obtained by Kuiper et al. (2001). The agreement between the results is quite satisfactory and differences are generally well within statistical uncertainties. We note that there is no general flux offset between ours and the COMPASS analysis, as observed above for the total Crab spectrum, and that differences are plausibly explained by the use of a different lists of viewing periods. Kuiper et al. (2001) noticed an enhanced emission in the 10–15 MeV energy interval for the Crab pulsar, and also in our analysis we found an equivalent feature. We note, however, that by shifting the phase interval definition by +0.03 (i.e. using the original phase interval definition of Kuiper et al. 2001) this spectral enhancement is considerably reduced in our analysis, suggesting that the enhancement is probably a statistical fluctuation rather than a physical feature.

We also fitted different spectral models to the data of the Crab pulsar and pulsar wind nebula, including power laws, exponentially cut-off power laws and curved power laws. We

determined the corresponding uncertainty bands using `ctbutterfly` and overlay them on the spectral points in Fig. 14. Using an exponentially cut-off power law for the nebula component instead of a simple power law improved the TS value of the nebula component by 6.5, corresponding to a detection significance of $2.5\sigma$ for the spectral cutoff. For the pulsar component no improvement was achieved when allowing for a cutoff or a curvature in the power law. For the Crab pulsar wind nebula, the best fitting parameters of Eq. (29) were $k = (13.7 \pm 0.5) \times 10^{-5}$ ph cm$^{-2}$ s$^{-1}$ MeV$^{-1}$, $\Gamma = 2.15 \pm 0.03$, and $E_c = 53.3 \pm 15.5$ MeV. For the Crab pulsar the best fitting power-law parameters were $k = (5.1 \pm 0.3) \times 10^{-5}$ ph cm$^{-2}$ s$^{-1}$ MeV$^{-1}$ and $\Gamma = 2.29 \pm 0.06$. This compares to the spectral indices of $2.227 \pm 0.013$ and $2.24 \pm 0.04$ determined by Kuiper et al. (2001) for the nebula and pulsar components using power-law models, respectively. While our index for the pulsar component is compatible with their result, our index for the nebula component is flatter, which can be explained by the spectral cutoff in our model. Using a simple power law for the nebula component, as Kuiper et al. (2001) did, we obtained a steeper index of $2.24 \pm 0.02$ that is compatible with their result.

### 4.3. Phase-resolved analysis of LS 5039

We now turn to an analysis of COMPTEL observations of the gamma-ray binary LS 5039 in order to validate the ability to conduct phase-resolved analyses with GammaLib and ctools. Using an orbit-resolved analysis, Collmar & Zhang (2014) found strong evidence that the MeV flux of GRO J1823-12, the strongest unidentified COMPTEL source in the Galactic plane, is modulated along the binary orbit of about 3.9 days of LS 5039. Specifically, using maximum likelihood significance maps, the authors demonstrated that GRO J1823-12 shows a more significant signal for the inferior conjunction period of LS 5039 as compared to the superior conjunction period. The same trend was also observed in their spectral analysis.

We repeated the analysis of Collmar & Zhang (2014) by choosing all viewing periods with pointing within 35° of $(l, b) = (17.5°, -0.5°)$ from the HEASARC database. In total this resulted in a list of 41 viewing periods, starting with viewing period 5.0 and ending with viewing period 712.0. Up to viewing period 712.0 our list is identical to Table 1 of Collmar & Zhang (2014), yet the authors included 12 more viewing periods in their analysis that are not available in the HEASARC archive. For our analysis we combined the data in a data space of $140 \times 123 \times 25$ bins of $1° \times 1° \times 2°$ in size and centred on $(l, b) = (15.0°, -4.5°)$, which corresponds to the same dimensions that were used by Collmar & Zhang (2014) in some of their analyses. We used the four standard COMPTEL energy bands for our analysis. Similar to Collmar & Zhang (2014) we used the binary ephemeris of Casares et al. (2005), which is an orbital period of 3.90603 days with periastron passage (corresponding to phase 0) at JD 2451943.09, and we define phases $0.45 \le \Phi < 0.9$ as the inferior conjunction interval (INFC) and phases $\Phi \ge 0.9$ and $\Phi < 0.45$ as the superior conjunction interval (SUPC).

As the first step we created TS maps of the region around LS 5039 using `comlixmap` for the INFC and SUPC phase intervals. The data for the four standard energy bands were analysed jointly, using a model composed of a test point source and an additional point source at the location of the quasar PKS 1830-210 that is spatially close to LS 5039. The spectra of both components were modelled using power laws. In addition, the source model included components for Galactic diffuse emission based
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Test statistic maps of LS 5039 for inferior conjunction (**left**) and superior conjunction (**right**) phase intervals derived by jointly fitting the four standard COMPTEL energy bands, covering 0.75–30 MeV. The location of LS 5039 is shown by a white plus symbol, and contours show the $1\sigma$, $2\sigma$, and $3\sigma$ location uncertainties. The quasar PKS 1830–210 was included in the source model and hence is not visible in the maps.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Spectral energy distributions and $1\sigma$ uncertainty bands of fitted power-law models for LS 5039 for the INFC (solid lines and filled area) and the SUPC (dashed lines and hatched area). The left panel compares the ctools and GammaLib results (red) to the results obtained by Collmar & Zhang (2014) for an identical event selection (blue). The right panel compares the ctools and GammaLib results (red) to the results obtained using the COMPASS software by excluding D2 modules with faulty PMTs (blue). Data points from Collmar & Zhang (2014) and derived using COMPASS were displaced by 3% in energy for clarity.</figcaption></figure>

on template maps for bremsstrahlung and inverse Compton emission with free scaling factors for each energy bin$^{11}$. Furthermore, an isotropic component was included to account for the cosmic gamma-ray background emission, with intensity fixed according to $I(E_\gamma) = 1.12 \times 10^{-4} (E_\gamma/5 \text{ MeV})^{-2.2} \text{ ph cm}^{-2} \text{ s}^{-1} \text{ MeV}^{-1} \text{ sr}^{-1}$ as suggested by Weidenspointner (1999). The background was modelled using the BGDLIXE algorithm with parameters $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 15$.

The resulting TS maps are shown in Fig. 15 for the INFC (left) and SUPC (right) phase intervals. The maps can be compared to those shown in Fig. 6 of Collmar & Zhang (2014), which

were determined separately for the 3–10 MeV and 10–30 MeV energy bands. In both analyses, LS 5039 is considerably more significant in the INFC phase interval but only weakly detected in the SUPC phase interval. In the latter interval, the emission maximum seems to be displaced towards the north-east with respect to the position of LS 5039 in ours and the analysis of Collmar & Zhang (2014), yet the emission location is still compatible within the $3\sigma$ uncertainty contour with the position of LS 5039.

As the next step we derived the spectral energy distribution of LS 5039 for the four standard energy bands using `csspec` for both phase intervals to reproduce Fig. 7 of Collmar & Zhang (2014). In addition, we also derived for both phase intervals the uncertainty band of the fitted power-law model using `ctbutterfly`. The results are shown in Fig. 16. For the left panel we used exactly the same event selection that was

<footnote n="11">11 We used the bremsstrahlung emission map with the COMPASS identifier MPE-MIS-13006 and the inverse Compton emission map with the COMPASS identifier ROL-MIS-163 that were also used by Collmar & Zhang (2014).</footnote>
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Orbital light curve of LS 5039 in the 10–30 MeV energy band. Data points for the ctools analysis were displaced by ±0.02 in phase for clarity. Analysis results obtained under analysis conditions identical to those in Collmar & Zhang (2014) ($\zeta_{min} = 0^{\circ}$ and *fpmtflag* = 2) are shown as triangles and dashed error bars, and results obtained using $\zeta_{min} = 5^{\circ}$ and *fpmtflag* = 0 are shown as dots and solid error bars. Vertical lines indicate the definitions of the INFC and SUPC phase intervals.</figcaption></figure>

used by Collmar & Zhang (2014) which is a minimum distance from the Earth horizon of $\zeta_{min} = 0^{\circ}$ and the use of circular exclusion regions to handle D2 modules with failed PMTs (i.e. *fpmtflag* = 2; cf. Appendix E). Only with this event selection we were able to reproduce the 1–3 MeV flux point of Collmar & Zhang (2014), while using our standard setting of $\zeta_{min} = 5^{\circ}$ and *fpmtflag* = 0 that excludes D2 modules with failed PMTs produced a notable variation of the 1–3 MeV flux point between inferior and superior conjunction. To verify that this variation can indeed be attributed to differences in the event selection, we also did an equivalent analysis with COMPASS using $\zeta_{min} = 5^{\circ}$ and *fpmtflag* = 0. This resulted in a 1–3 MeV flux variation between inferior and superior conjunction that was similar to that observed in our analysis, confirming that the spectral differences are due to differences in the event selection. The corresponding results are summarised in the right panel of Fig. 16.

As illustrated by the uncertainty bands of the fitted power-law model, the use of circular exclusion regions for D2 modules with faulty PMTs leads to a flatter (or softer) spectrum in superior conjunction, yet the variation seems still to remain within statistical uncertainties, given the broadness of the uncertainty band. After all, LS 5039 is a very faint source in SUPC, and hence its spectral properties are only poorly constrained in this phase interval.

Finally, we derived the orbital flux variation in the 10–30 MeV energy band for the two different event selections to reproduce Fig. 8 of Collmar & Zhang (2014). For this purpose we split the 10–30 MeV data into five phase intervals of equal length and fitted the data using comlixfit with a source model comprising components for LS 5039, PKS 1830–210, bremsstrahlung emission, inverse Compton emission and cosmic background emission (see above for details). The background was modelled using the BGDLIXE algorithm with parameters $N_{avgr} = 5$ and $N_{incl} = 15$. The results are shown in Fig. 17.

Our analysis reproduces the orbital light curve of Collmar & Zhang (2014) within statistical uncertainties. Differences between their flux points and ours can be explained by differences in the event selection, as Collmar & Zhang (2014) use a larger number of viewing periods compared to our analysis.

<table><caption>**Table 3.** Energy bins for 1.8 MeV gamma-ray line analysis.</caption><thead><tr><td>Number</td><td>Energies (MeV)</td></tr></thead><tbody><tr><td>1</td><td>1.000–1.584</td></tr><tr><td>2</td><td>1.584–1.634</td></tr><tr><td>3</td><td>1.634–1.684</td></tr><tr><td>4</td><td>1.684–1.734</td></tr><tr><td>5</td><td>1.734–1.784</td></tr><tr><td>6</td><td>1.784–1.834</td></tr><tr><td>7</td><td>1.834–1.884</td></tr><tr><td>8</td><td>1.884–1.934</td></tr><tr><td>9</td><td>1.934–1.984</td></tr><tr><td>10</td><td>1.984–2.034</td></tr><tr><td>11</td><td>2.034–3.000</td></tr></tbody></table>

**Notes.** Bins 1 and 11 serve primarily to constrain continuum gamma-ray emission, while bins 2–10 serve to trace the shape of the 1.8 MeV gamma-ray line.

Variations of the same size are also observed in our analysis for the two different event selections, which, however, are well within statistical uncertainties, as expected.

### 4.4. $^{26}$Al line emission from Carina

As the next step we validated the capacities of GammaLib and ctools for gamma-ray line emission analysis together with its abilities to assess the spatial morphology of the emission. As reference, we chose the COMPTEL detection of point-like 1.8 MeV line emission from the Carina region for this validation, as reported by Knödlseder (1994) and Knödlseder et al. (1996a), which are to our knowledge the only published studies where an analysis using a parametric spatial model was performed with COMPTEL.

In their studies the authors analysed COMPTEL data from viewing periods 1 to 301 combined in a data space of $100 \times 100 \times 25$ bins of $1^{\circ} \times 1^{\circ} \times 2^{\circ}$ in size, centred on $(l, b) = (286.5^{\circ}, 0.5^{\circ})$, which corresponds to the peak position of the observed 1.8 MeV line emission feature. The emission feature was found to be compatible with a point-like source with an 1.8 MeV flux of $(3.1 \pm 0.8) \times 10^{-5} \text{ ph cm}^{-2} \text{ s}^{-1}$ that was determined through model fitting. Using fits of models with uniform intensity within a circular region centred on $(l, b) = (286.5^{\circ}, 0.5^{\circ})$, Knödlseder et al. (1996a) determined a $2\sigma$ upper limit of $5.6^{\circ}$ for the diameter of the 1.8 MeV emission region. The analysis was done in a single energy bin covering $1.7–1.9$ MeV, and the instrumental background was modelled using measurements in adjacent energy intervals that were corrected for the energy dependence of the Compton scatter angle $\bar{\phi}$. This method suppresses to first order emission from continuum gamma-ray sources (Knödlseder et al. 1996b).

We analysed the same data and adopted the same data space definition that was used by Knödlseder et al. (1996a), yet we tested an alternative analysis method that jointly handles the $^{26}$Al line signal and any underlying continuum emission. This is more in line with the GammaLib and ctools philosophy of explicitly modelling all emission components and provides a more accurate handling of underlying continuum emission. Specifically, we split the data within the 1–3 MeV energy interval into 11 energy bins, specified in Table 3, and analysed them jointly using a combination of model components describing the 1.8 MeV line signal, any underlying continuum gamma-ray emission and the instrumental background. To model the 1.8 MeV line emission
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Count spectrum of the Carina region determined for 11 bins covering the energy band 1–3 MeV. The top panel shows the measured number of counts per MeV together with the best fitting background model (green dashed), composed of instrumental background (green dotted) and cosmic gamma-ray background. The combined source and background model is shown as a black solid line. The bottom panel shows the background model-subtracted count spectrum together with the best fitting 1.8 MeV line model (red solid) on top of the combined Galactic continuum components (blue dashed), composed of bremsstrahlung (blue dotted) and inverse Compton emission (blue dash-dotted).</figcaption></figure>

spectrum we used a Gaussian spectral component with a fixed mean of 1.809 MeV and a standard deviation of $\sigma = 58.9$ keV that corresponds to COMPTEL's instrumental energy resolution at that energy. Similarly to our analysis of LS 5039 we modelled continuum emission using template maps for Galactic bremsstrahlung and inverse Compton emission and an isotropic component for the cosmic gamma-ray background. The spectra of the three continuum components were model using power laws, where the prefactors and indices were free for the Galactic components, while the prefactor and index were fixed to $I(E_\gamma) = 1.12 \times 10^{-4}(E_\gamma/5 \text{ MeV})^{-2.2} \text{ ph cm}^{-2} \text{ s}^{-1} \text{ MeV}^{-1} \text{ sr}^{-1}$ for the cosmic gamma-ray background component, as suggested by Weidenspointner (1999). The TS map was generated using `comlixmap` and model fitting was done using `comlixfit` with the standard BGDLIX parameters $N_{\text{avgr}} = 5$ and $N_{\text{incl}} = 15$.

We illustrate our analysis procedure in Fig. 18, which shows the counts spectrum determined using Eq. (24) for the position $(l, b) = (286.5^\circ, 0.5^\circ)$ at which Knödlseder (1994) found the 1.8 MeV line emission maximum and an ARM window of $\pm 3^\circ$. We also show the model components that were fitted using `comlixfit` and that we extracted using Eq. (25). We used a point source located at the fixed position $(l, b) = (286.5^\circ, 0.5^\circ)$ as spatial model for the 1.8 MeV line component. Figure 18 illustrates that the data are dominated by instrumental background, followed by cosmic gamma-ray background. The bottom panel illustrates that, once these components are subtracted, a clear line signal becomes apparent that is compatible with the expected signature of the $^{26}$Al decay line. In addition, a continuum signal is detected that is dominated by Galactic bremsstrahlung emission.

In their Fig. 2, Knödlseder et al. (1996a) present a 1.8 MeV line emission maximum likelihood map of the Carina region, and we produced an equivalent TS map with the same spatial binning using `comlixmap`. Specifically, we fitted point-source models for the 1.8 MeV line emission and the 1–3 MeV continuum emission for a grid of source positions, producing hence TS maps for both emission components. In the fitting the fluxes of

both point sources were constrained to non-negative values. The resulting maps are shown in Fig. 19, where the left map can be compared to the maximum likelihood map presented in Fig. 2 of Knödlseder et al. (1996a). Both maps show qualitatively comparable features, yet we obtained a maximum TS value of 23.1 at $(l, b) = (285.5^\circ, 1.5^\circ)$ while Knödlseder et al. (1996a) found a lower maximum TS value of 14.7 at $(l, b) = (286.5^\circ, 0.5^\circ)$. Eventually, these differences may be explained by the background modelling techniques and analysis methods that differ significantly between the studies. We recall that Knödlseder et al. (1996a) analyse the 1.8 MeV line data in a single energy bin covering 1.7–1.9 MeV and used a background model derived from adjacent energy bands, which to first order includes also continuum emission, but which does not properly account for spectral differences between instrumental background and diffuse emission components as well as differences in their $\bar{\phi}$ distributions (Bloemen et al. 1999).

It is actually rather reassuring that both approaches produce qualitatively comparable maps, as already suggested by Bloemen et al. (1999) who implemented a comparable analysis method to ours. The continuum TS map indicates that the continuum emission is located towards the Galactic plane, suggesting that it originates from our Galaxy. We emphasise, however, that the statistical significance of the emission features is modest, and consequently the appearance of the map is notably affected by the statistical fluctuations of the data. Nevertheless, we note that observed TS maxima are not too far from maxima in the Galactic bremsstrahlung emission, which eventually may be the dominant MeV emission component near the Galactic plane (Strong et al. 1996).

As the next step we fitted the data with a source model composed of a 1.8 MeV line component modelled using a point source and a 1–3 MeV continuum component modelled using a combination of bremsstrahlung and inverse Compton spatial maps as well as an isotropic component for the cosmic gamma-ray background. The position of the point-source model as well as the spectral parameters for the Galactic continuum power-law components were free parameters in the fit. Fitting the data using `comlixfit` gave a best fitting position of $(l, b) = (285.4^\circ \pm 0.8^\circ, 1.3^\circ \pm 0.7^\circ)$, a flux of $(3.1 \pm 0.6) \times 10^{-5} \text{ ph cm}^{-2} \text{ s}^{-1}$ and a TS value of 22.3 for the 1.8 MeV line component. Our 1.8 MeV line flux is consistent with the value of $(3.1 \pm 0.8) \times 10^{-5} \text{ ph cm}^{-2} \text{ s}^{-1}$ found by Knödlseder et al. (1996a), yet our best fitting position is offset by about $1.4^\circ$ from the one found in their analysis. Replacing the point source for the 1.8 MeV line component by a radial disk model did not improve the fit and led to a best fitting disk radius of $0.003^\circ$, which is near the minimum value of $0.001^\circ$ that we allowed in the analysis. Using the `ctulimit` tool we derived a $2\sigma$ upper limit of $5.1^\circ$ for the diameter of the disk, a bit smaller than the value of $5.6^\circ$ determined by Knödlseder et al. (1996a). This difference is plausibly explained by the larger detection significance of the 1.8 MeV line emission signal in our analysis, allowing us to put a stronger constraint on the extent of the emission region.

Finally, we tried to reproduce Fig. 5.5 of Knödlseder (1994) that shows the variation of the maximum log-likelihood value and the fitted 1.8 MeV line flux as a function of spatial disk radius. For that purpose, we determined the maximum likelihood solution for a disk model with position fixed at our maximum likelihood solution $(l, b) = (285.4^\circ, 1.3^\circ)$ for a set of disk radius values starting with $0.001^\circ$, and followed by $1^\circ$ up to $10^\circ$ with a step size of $1^\circ$. The results of our analysis are shown in Fig. 20, which also includes the results shown in Fig. 5.5 of Knödlseder (1994) for comparison. The log-likelihood profile between both
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Test statistic maps of 1.8 MeV line emission (**left**) and 1–3 MeV continuum emission (**right**) obtained using comlixmap for the Carina region. The map on the left is equivalent to Fig. 2 of Knödlseder et al. (1996a). Dashed contours in the map on the right reflect the intensity of a combination of Galactic bremsstrahlung and inverse Compton emission maps as fitted in an independent analysis where a model of a 1.8 MeV line point source was fitted together with a combination of bremsstrahlung and inverse Compton spatial maps to the data (see text).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Variation in maximum log-likelihood with increasing radius of the disk model (**left**) and fitted 1.8 MeV line flux (**right**). Data shown as red filled dots and solid lines are those derived in this work, and data shown as blue open dots and dashed lines are taken from Knödlseder (1994).</figcaption></figure>

analyses is very similar, yet our analysis turns out to be a bit more constraining, probably owing to the more significant detection of the 1.8 MeV line emission feature with respect to the analysis of Knödlseder (1994). The flux attributed to the 1.8 MeV line component changes actually very little with disk radius in our analysis, while in Knödlseder (1994) the flux decreases with increasing disk radius. This difference is probably partly due to fact that our analysis detects the 1.8 MeV line signal more significantly, but may also be related to our analysis method that gives freedom to the continuum emission model components to adjust as a function of the 1.8 MeV disk model radius, while in Knödlseder (1994) the continuum emission was implicitly subtracted by the background model, resulting in a much more constrained overall model that leaves little freedom for the model to adjust.

**5. Conclusion**

We have implemented a comprehensive science analysis framework for COMPTEL gamma-ray data in the GammaLib and

ctools software packages, and we have demonstrated that our software reliably reproduces published analysis results that were derived using the COMPASS software. Having public, free, and validated software for COMPTEL science data analysis now opens the HEASARC COMPTEL archive to the community for scientific exploration. In the medium-energy gamma-ray band, from 1–30 MeV, the COMPTEL archive still contains the most sensitive observations ever performed, and a unique dataset for exploring the non-thermal Universe and nuclear transition lines. Since the 1990s when the COMPTEL data were taken, the field of gamma-ray astronomy has made impressive progress thanks to satellites such as the International Gamma-Ray Astrophysics Laboratory (INTEGRAL) and *Fermi* and ground-based observatories such as the High Energy Stereoscopic System (H.E.S.S.), the Major Atmospheric Gamma-ray Imaging Cherenkov Telescope (MAGIC), and the Very Energetic Radiation Imaging Telescope Array System (VERITAS). Many source classes that are known today were not established as gamma-ray emitters during the COMPTEL era, and the COMPTEL data were never comprehensively analysed with the current knowledge in the field. An exception to this is the post-mission discovery of the orbital modulation of MeV gamma-ray emission of LS 5039 in the COMPTEL archival data by Collmar & Zhang (2014), a source that was not an established gamma-ray emitter in the 1990s. Thanks to GammaLib and ctools, such discoveries are now achievable by the community at large.

We want to stress that our work was also motivated by the goal of reducing the carbon footprint of astronomical research. As recently pointed out by Knödlseder et al. (2022), the current deployment rate of new astronomical observatories is not compatible with the imperative of reducing the carbon footprint across all activity sectors of modern societies, and this calls for fundamental changes in astronomical practices in the future. Among the many possible alternatives to the building of ever more and ever bigger new astronomical facilities is the exploitation of archival data from past missions that may have scientific treasures yet to be discovered. Since version 2.0.0, ctools estimates the carbon footprint of its use based on the assessment of Berthoud et al. (2020) for the GRICAD computing centre, and, using this feature, we estimate that the work presented in this publication resulted in the generation of 600 ± 300 kgCO₂e of
greenhouse gases (GHGs) due to the computing related to the
analysis of archival data. This is about 40 times less than the
median per-publication emissions associated with the analysis of
data from an active astronomical observatory (Knödlseder et al.
2022). Taking all sources of emissions related to this work into
account, we estimate the carbon footprint of this research to be
1.9 ± 0.3 tCO₂e (see Appendix H).

This assessment illustrates that the exploitation of archival
data instead of the development of new astronomical observato-
ries has the potential to dramatically reduce the carbon footprint
of astronomical research, which would help realise the reduc-
tions that are needed to limit global warming and reach the goals
of the Paris Climate Agreement. Past missions can be seen as
our scientific carbon legacy since the GHGs that were emitted
during their construction and operations are to a large extent still
present in the Earth's atmosphere. To assure that these GHGs
were at least not emitted in vain, conservation of the archival
data and the development and maintenance of software for their
exploitation should therefore be of the utmost importance.

Acknowledgements. We would like to thank the anonymous referee for the very careful reading of the manuscript and the many constructive suggestions. This research made use of ctools, a community-developed gamma-ray astronomy science analysis software (Knödlseder et al. 2016). ctools is based on GammaLib, a community-developed toolbox for the scientific analysis of astronomical gamma-ray data (Knödlseder et al. 2011). This work has made use of the Python 2D plotting library matplotlib (Hunter 2007). This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research is part of the #LowCarbon-Science initiative that aims in reducing the carbon footprint of scientific research, and has benefitted from discussions held in the context of the GDR Labos 1point5 and the Astro4Earth initiative.

References

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Berthoud, F., Bzeznik, B., Gibelin, N., et al. 2020, Estimation de l'empreinte carbone d'une heure.coeur de calcul, Research report, UGA - Université Grenoble Alpes ; CNRS ; INP Grenoble ; INRIA

Bloemen, H., Diehl, R., Hermsen, W., Knödlseder, J., & Schönfelder, V. 1999, Astrophys. Lett. Commun., 38, 391

Bloemen, H., Hermsen, W., Swanenburg, B. N., et al. 1994, ApJS, 92, 419

Casares, J., Ribó, M., Ribas, I., et al. 2005, MNRAS, 364, 899

Cash, W. 1979, ApJ, 228, 939

Coleman, A., & McConnell, M. 2020, in American Astronomical Society Meeting Abstracts, 235, 206.04

Collmar, W., & Zhang, S. 2014, A&A, 565, A38

Collmar, W., Wessolowski, U., Schönfelder, V., et al. 1997, in American Institute of Physics Conference Series, 410, Proceedings of the Fourth Compton Symposium, eds. C. D. Dermer, M. S. Strickman, & J. D. Kurfess, 1587

de Vries, C. P., & COMPTEL collaboration. 1994, in Astronomical Society of the Pacific Conference Series, 61, Astronomical Data Analysis Software and Systems III, eds. D. R. Crabtree, R. J. Hanisch, & J. Barnes, 399

Diehl, R., Aarts, H., Bennett, K., et al. 1992, in Data Analysis in Astronomy, 201

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

IPBES 2019, Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, eds. E. S. Brondizio, J. Settele, S. Díaz, & H. T. Ngo (IPBES secretariat, Bonn, Germany)

IPCC 2021, Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, eds. V. Masson-Delmotte, P. Zhai, A. Pirani, et al. (Cambridge University Press)

IPCC 2022, Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, eds. H.-O. Pörtner, M. D.C., Roberts, E. Tignor, et al. (Cambridge University Press)

Knödlseder, J. 1994, PhD thesis, Technical University Munich, Germany

Knödlseder, J., Bennett, K., Bloemen, H., et al. 1996a, A&As, 120, 327

Knödlseder, J., von Ballmoos, P., Diehl, R., et al. 1996b, SPIE Conf. Ser., 2806, 386

Knödlseder, J., Mayer, M., Deil, C., et al. 2016, Astrophysics Source Code Library [ascl:1601.0005]

Knödlseder, J., Mayer, M., Deil, C., et al. 2011, Astrophysics Source Code Library [ascl:1110.0007]

Knödlseder, J., Brau-Nogué, S., Coriat, M., et al. 2022, Nat. Astron., 6, 503

Kuiper, L., Hermsen, W., Cusumano, G., et al. 2001, A&A, 378, 918

Landrigan, P., Fuller, R., Acosta, N., et al. 2017, Lancet, 391

Martin, P., Brau-Nogué, S., Coriat, M., et al. 2022, ArXiv e-prints, [arXiv:2204.12362]

Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396

Much, R., Bennett, K., Buccheri, R., et al. 1995a, A&A, 299, 435

Much, R., Bennett, K., Buccheri, R., et al. 1995b, Adv. Space Res., 15, 81

Much, R., Harmon, B. A., Nolan, P., et al. 1996, A&As, 120, 703

Pence, W. D., Chiappetti, L., Page, C. G., Shaw, R. A., & Stobie, E. 2010, A&A, 524, A42

Schönfelder, V., Aarts, H., Bennett, K., et al. 1993, ApJS, 86, 657

Stacy, J. G., Kippen, R. M., Kappadath, S. C., et al. 1996, A&As, 120, 691

Strong, A., & Collmar, W. 2019, Mem. Soc. Astron. Italiana, 90, 297

Strong, A. W., Bennett, K., Bloemen, H., et al. 1996, A&As, 120, 381

van der Meulen, R. D., Bloemen, H., Bennett, K., et al. 1998, A&A, 330, 321

van Dijk, R. C. A. 1996, PhD thesis, University of Amsterdam, The Netherlands

Weidenspointner, G. 1999, PhD thesis, Technical University Munich, Germany

Weidenspointner, G., Varendorff, M., Oberlack, U., et al. 2001, A&A, 368, 347
## Appendix A: HEASARC archive

COMPTEL data in the HEASARC archive are grouped by so-called viewing periods with typical durations of two weeks during which the CGRO satellite had a stable pointing. In total, there exist 359 viewing periods, of which 278 were archived at HEASARC, covering viewing periods 1.0 to 719.0, where the first digit indicates the mission year or phase, with viewing periods smaller than 100 corresponding to mission year one. Some of the viewing periods that are available in the HEASARC archive cannot be exploited due to corrupted or missing data files, or spacecraft operations that prevented data taking. The 23 viewing periods in the range 1.0 to 719.0 that are not exploitable are listed in Table A.1. This leaves 255 exploitable viewing periods in the HEASARC archive, which is about 71% of the total number of 359 viewing periods that were executed during the CGRO mission.

<table><caption>**Table A.1.** Inexploitable viewing periods in the range 1.0 to 719.0 in the COMPTEL HEASARC archive.</caption><thead><tr><th>VP</th><th>Target</th><th>Reason</th></tr></thead><tbody><tr><td>22.0</td><td>Mrk 279</td><td>Missing EVP file</td></tr><tr><td>29.0</td><td>G 224 - 40</td><td>EVP file unreadable</td></tr><tr><td>225.0</td><td>CGRO reboost</td><td>Missing EVP file</td></tr><tr><td>229.3</td><td>Perseids shower</td><td>Missing EVP file</td></tr><tr><td>229.5</td><td>Gal 5 + 5</td><td>Missing EVP file</td></tr><tr><td>303.4</td><td>Not executed</td><td>Missing EVP file</td></tr><tr><td>303.5</td><td>Reboost</td><td>Missing EVP file</td></tr><tr><td>303.7</td><td>Nova Cygnus 92</td><td>Missing EVP file</td></tr><tr><td>308.3</td><td>Reboost</td><td>Missing EVP file</td></tr><tr><td>308.6</td><td>Virgo 283 + 75</td><td>Missing EVP file</td></tr><tr><td>311.3</td><td>Reboost</td><td>Missing EVP file</td></tr><tr><td>311.6</td><td>Virgo 284 + 75</td><td>Missing EVP file</td></tr><tr><td>426.0</td><td>Anticenter</td><td>Invalid event times in EVP</td></tr><tr><td>617.2</td><td>CGRO reboost test</td><td>Missing EVP file</td></tr><tr><td>617.4</td><td>CGRO reboost test</td><td>Missing EVP file</td></tr><tr><td>617.6</td><td>CGRO reboost</td><td>Missing EVP file</td></tr><tr><td>619.0</td><td>Cir X-1</td><td>Missing EVP file</td></tr><tr><td>619.5</td><td>CGRO reboost</td><td>Missing EVP file</td></tr><tr><td>627.0</td><td>PSR 1055 - 52</td><td>Missing EVP file</td></tr><tr><td>632.0</td><td>Not executed</td><td>Missing EVP file</td></tr><tr><td>632.1</td><td>PKS 0235 + 164</td><td>Missing EVP file</td></tr><tr><td>711.0</td><td>2EG J1835 + 5919</td><td>Missing OAD files</td></tr><tr><td>718.0</td><td>Cen X-3</td><td>Missing EVP file</td></tr></tbody></table>

The HEASARC archive was created from the original COMPTEL data using dedicated file conversion software that generated FITS files following recognised standards (Pence et al. 2010). While interpreting these FITS files using GammaLib we found that the units in some of the table columns were not correct. We summarise the inconsistencies that we encountered in Table A.2. We furthermore found that the good time interval dataset for viewing period 8.0 with the identifier MPE-TIM-11481 only contains about half of the good time intervals that were defined in COMPASS, probably as the result of a file truncation error in the input ASCII file that was used for the generation of the FITS file in the HEASARC archive.

The HEASARC archive also comprises binned data products for standard energy bands that can be used for analysis using GammaLib, yet we note that the world coordinate system information attached to these data products is incorrect. While the FITS headers suggest that the data cubes are provided in Mercator projection, the event cubes are presented in a cartesian pro-

<table><caption>**Table A.2.** Incorrect header units in FITS files in the COMPTEL HEASARC archive.</caption><thead><tr><th>File</th><th>Column</th><th>Written unit</th><th>True unit</th></tr></thead><tbody><tr><td>EVP</td><td>TOF</td><td>ns</td><td>channels</td></tr><tr><td>EVP</td><td>PHIBAR</td><td>deg</td><td>rad</td></tr><tr><td>EVP</td><td>GLON_SCAT</td><td>deg</td><td>rad</td></tr><tr><td>EVP</td><td>GLAT_SCAT</td><td>deg</td><td>rad</td></tr><tr><td>EVP</td><td>AZIMUTH_SCAT</td><td>deg</td><td>rad</td></tr><tr><td>EVP</td><td>ZENITH_SCAT</td><td>deg</td><td>rad</td></tr><tr><td>EVP</td><td>EARTH_HORIZON</td><td>deg</td><td>rad</td></tr><tr><td>BVC</td><td>SSB_X</td><td>km</td><td>milli light seconds</td></tr><tr><td>BVC</td><td>SSB_Y</td><td>km</td><td>milli light seconds</td></tr><tr><td>BVC</td><td>SSB_Z</td><td>km</td><td>milli light seconds</td></tr></tbody></table>

jection with a reference latitude value of 0. For all event cubes given in Mercator projection GammaLib assumes that they originate from the HEASARC archive and corrects the world coordinate system information upon reading of the event cubes.

The HEASARC archive comprises further data products such as telescope housekeeping data, gamma-ray burst detector data and sky maps that are not used by GammaLib.

## Appendix B: Time of flight conversion

The HEASARC archive mixes different versions of EVP files that have different levels of processing for the ToF values. The versions can be distinguished by the header keyword DSD_REP in the EVP file, specifying either 2 for ToF<sub>II</sub> or 3 for ToF<sub>III</sub>, where the latter corrects for energy dependent effects, aligning the forward peak at channel number 120 (see van Dijk 1996 and Weidenspointner 1999 for an explanation of the ToF corrections). If a version 2 EVP is encountered, the GOMEEventList class will automatically convert ToF<sub>II</sub> values into ToF<sub>III</sub> using

$$ \text{ToF}_{\text{III}} = \text{ToF}_{\text{II}} + 120 - \left( \sum_{i=0}^{6} a_i E_1^i + \sum_{i=0}^{4} b_i E_2^i - 118.3 \right) \quad (\text{B.1}) $$

with the coefficients given by Tables B.1 and B.1 and the D1 and D2 energy deposits $E_1$ and $E_2$ given in MeV. For the HEASARC archive, ToF values accessed through GammaLib are therefore always ToF<sub>III</sub> values.

<table><caption>**Table B.1.** Coefficients, $a_i$, for the conversion from ToF<sub>II</sub> to ToF<sub>III</sub>.</caption><thead><tr><th>E<sub>1</sub></th><th>&le; 2.25 MeV</th><th>&gt; 2.25 MeV</th></tr></thead><tbody><tr><td>a<sub>0</sub></td><td>111.74858</td><td>116.25374</td></tr><tr><td>a<sub>1</sub></td><td>28.280247</td><td>0.500407092</td></tr><tr><td>a<sub>2</sub></td><td>-45.024305</td><td>0.3818272</td></tr><tr><td>a<sub>3</sub></td><td>35.18321</td><td>-0.080145513</td></tr><tr><td>a<sub>4</sub></td><td>-14.639463</td><td>0.006556979</td></tr><tr><td>a<sub>5</sub></td><td>3.1342536</td><td>-0.00024650067</td></tr><tr><td>a<sub>6</sub></td><td>-0.2711735</td><td>3.507724 &times; 10<sup>-6</sup></td></tr></tbody></table>

## Appendix C: Flux correction due to time of flight selection

To correct for the photon rejection by the ToF selection an energy-dependent correction factor needs to be applied to the instrument response. Since this correction factor depends on the ToF selection interval, it is computed in GammaLib during the
<table><caption>**Table B.1.** Coefficients, $b_i$, for the conversion from $ToF_{II}$ to $ToF_{III}$.</caption><thead><tr><th>E<sub>2</sub></th><th>&le; 1.4 MeV</th><th>1.4 - 5.5 MeV</th><th>&gt; 5.5 MeV</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>181.77024</td><td>120.91608</td><td>119.24278</td></tr><tr><td>b<sub>1</sub></td><td>-252.4107</td><td>-0.1504849</td><td>-0.43134699</td></tr><tr><td>b<sub>2</sub></td><td>371.09898</td><td>-0.45526025</td><td>0.06018308</td></tr><tr><td>b<sub>3</sub></td><td>-232.83985</td><td>0.11710009</td><td>-0.002684779</td></tr><tr><td>b<sub>4</sub></td><td>52.918785</td><td>-0.0082172427</td><td>3.7720986 &times; 10<sup>5</sup></td></tr></tbody></table>

event binning and stored using the header keyword TOFCOR in the DRE FITS files. The correction factor is determined from a linear interpolation of the values given in Table C.1, evaluated at the geometric mean energy $\sqrt{E_{min}E_{max}}$ of the energy interval $[E_{min}, E_{max}]$ for which the event cube is generated. Only corrections for $ToF_{max}$ = 130 are supported.

<table><caption>**Table C.1.** Time of flight correction factors for $ToF_{max}$ = 130.</caption><thead><tr><th rowspan="2">ToF<sub>min</sub></th><th colspan="4">Geometric mean energy (MeV)</th></tr><tr><th>0.8660</th><th>1.7321</th><th>5.4772</th><th>17.3205</th></tr></thead><tbody><tr><td>110</td><td>1.14</td><td>1.07</td><td>1.02</td><td>1.01</td></tr><tr><td>111</td><td>1.17</td><td>1.09</td><td>1.03</td><td>1.01</td></tr><tr><td>112</td><td>1.21</td><td>1.11</td><td>1.05</td><td>1.02</td></tr><tr><td>113</td><td>1.26</td><td>1.15</td><td>1.07</td><td>1.04</td></tr><tr><td>114</td><td>1.32</td><td>1.20</td><td>1.11</td><td>1.06</td></tr><tr><td>115</td><td>1.40</td><td>1.27</td><td>1.17</td><td>1.11</td></tr><tr><td>116</td><td>1.50</td><td>1.36</td><td>1.24</td><td>1.17</td></tr><tr><td>117</td><td>1.63</td><td>1.47</td><td>1.35</td><td>1.28</td></tr><tr><td>118</td><td>1.79</td><td>1.63</td><td>1.51</td><td>1.43</td></tr><tr><td>119</td><td>2.01</td><td>1.85</td><td>1.73</td><td>1.67</td></tr></tbody></table>

## Appendix D: Pulsar timing

GammaLib supports generation of event cubes, geometry functions and exposure maps for phase-resolved pulsar analysis. For this purpose a specific processing is implemented in the methods `GCOMDri::compute_dre`, `GCOMDri::compute_drg`, and `GCOMDri::compute_drx` that is used if the specified `GCOMSelection` instance includes pulsar ephemerides data and the specification of pulsar phase intervals. All three methods will first trim the good time intervals of the observation so that they only cover periods for which the specified pulsar ephemerides are valid. In that way, GammaLib assures that only data will be used for an analysis that cover intervals with valid pulsar ephemeris information.

The remaining pulsar-specific code is implemented in `GCOMDri::compute_dre`. First, the method converts arrival times $t_{CGRO}$ of COMPTEL events at the CGRO satellite, specified in the Coordinated Universal Time (UTC) time system, into arrival times $t_{SSB}$ at the Solar System barycentre, specified in the barycentric dynamical time (TDB) system.

As a side note, before 1992-06-25T01:00:00 UTC the CGRO onboard clock was early by 2.042144 seconds, and this time difference needs to be subtracted from the measured onboard time to get the true arrival time in UTC. In GammaLib, this subtraction is automatically performed when converting onboard times, given in truncated Julian days (TJDs) and tics, into `GTime` objects using the `gammalib::com_time` function, and the conversion is undone when using the inverse functions `gammalib::com_tjd` and `gammalib::com_tics`.

After applying the clock correction, event times are converted using

$$t_{SSB} = t_{CGRO} + \Delta t_{CGRO \to SSB} - \Delta t_{Shapiro} + \Delta t_{UTC \to TT} + \Delta t_{TT \to TBD}, \quad (D.1)$$

where $\Delta t_{CGRO \to SSB}$ corrects for the light travel time from CGRO to the Solar System barycentre, $\Delta t_{Shapiro}$ corrects for the gravitational time delay near the Sun, also known as Shapiro delay, $\Delta t_{UTC \to TT}$ converts from UTC to terrestrial time (TT), and $\Delta t_{TT \to TBD}$ converts from TT to the TDB time system. We note that all time correction terms are themselves time dependent, yet we ignore this time dependence in Eq. (D.1) to simplify the notation. All terms are given in units of seconds.

The computation of the three terms $\Delta t_{CGRO \to SSB} - \Delta t_{Shapiro} + \Delta t_{TT \to TBD}$ is implemented by the `Gephemerides::geo2ssb` method, which regroups all time corrections that depend on planetary ephemerides. GammaLib deals with planetary ephemerides through the `Gephemerides` class, and the software includes the Jet Propulsion Laboratory (JPL) Development Ephemeris 200 including information about the Sun position, the Earth position and its first three time derivatives, and the time difference between TT and TDB in seconds on a daily basis between Julian day (JD) 2436913 (10 December 1959) and 2469807 (31 December 2049). Specifically, the method computes

$$\Delta t_{CGRO \to SSB} = \hat{\mathbf{n}} \cdot (\mathbf{e} + \mathbf{r}) \quad (D.2)$$

and

$$\Delta t_{Shapiro} = -2t_{\odot} \ln \left( 1 + \frac{\hat{\mathbf{n}} \cdot (\mathbf{e} + \mathbf{r} - \mathbf{s})}{|\mathbf{e} + \mathbf{r} - \mathbf{s}|} \right), \quad (D.3)$$

where

$$\hat{\mathbf{n}} = \begin{bmatrix} \cos \delta \cos \alpha \\ \cos \delta \sin \alpha \\ \sin \delta \end{bmatrix} \quad (D.4)$$

is the unit vector in celestial coordinates of a pulsar with right ascension $\alpha$ and declination $\delta$, $\mathbf{e}$ is the vector from the Solar System barycentre to the centre of the Earth, $\mathbf{s}$ is the vector from the Solar System barycentre to the centre of the Sun, $\mathbf{r}$ is the vector from the centre of the Earth to the CGRO spacecraft, and

$$t_{\odot} = 4.92549 \times 10^{-6} \quad (D.5)$$

is half the Schwarzschild radius of the Sun divided by the speed of light. All vectors are given in units of light seconds and are specified in the celestial coordinate system.

The Earth vector $\mathbf{e}$ at a given time $t_{CGRO}$ (specified in the TT time system) is computed using the Taylor expansion

$$\mathbf{e} = \mathbf{e}(i) + \dot{\mathbf{e}}(i)\Delta t + \frac{1}{2}\ddot{\mathbf{e}}(i)\Delta t^2 + \frac{1}{6}\dddot{\mathbf{e}}(i)\Delta t^3, \quad (D.6)$$

where $i$ is the index of the nearest entry in time in the JPL database, $\mathbf{e}(i)$, $\dot{\mathbf{e}}(i)$, $\ddot{\mathbf{e}}(i)$, and $\dddot{\mathbf{e}}(i)$ is the Earth position and its first three time derivatives for this entry, and

$$\Delta t = t_{CGRO} - t(i) \quad (D.7)$$

is the time difference between the event time and the nearest entry in the JPL database in units of days, which by definition is in the range -0.5 and +0.5. Since the Sun moves only a little around the Solar System barycentre, it is sufficient to take for the Sun vector the nearest entry in the JPL database, which is
s = s(i). Finally, the conversion from the TT to the TDB time system at a given time t<sub>CGRO</sub> is computed using

$$
\Delta t_{TT \to TDB} = \Delta t_{TT \to TDB}(i) + (\Delta t_{TT \to TDB}(i+1) - \Delta t_{TT \to TDB}(i))\Delta t, \quad (D.8)
$$

where $\Delta t_{TT \to TDB}(i)$ is the nearest entry in the JPL database.

The last term in the time correction, $\Delta t_{UTC \to TT}$, does not
depend on planetary ephemerides and is given by

$$
\Delta t_{UTC \to TT} = 32.184 + n_{\text{leap}}, \quad (\text{D.9})
$$

where $n_{\text{leap}}$ is the number of leap seconds. The computation of
$\Delta t_{\text{UTC}\to\text{TT}}$ is implemented by the `GTime::utc2tt` method.

The HEASARC archive includes for most of the viewing
periods files that provide for each superpacket the vectors **e** + **r**
as well as the correction terms Δt<sub>UTC→TT</sub> + Δt<sub>TT→TBD</sub>, avoid-
ing the need for planetary ephemerides. These so-called BVC
data can be handled by GammaLib through the classes GCOMBvc
and GCOMBvcs that manage individual data records as well as
entire files. Specifically, the method GCOMBvcs::tdelta com-
putes Δt<sub>CGRO→SSB</sub> − Δt<sub>Shapiro</sub> + Δt<sub>UTC→TT</sub> + Δt<sub>TT→TBD</sub>, and if the
observation that should be binned includes BVC information,
GCOMDri::compute_dre will use this method instead of the
algorithm described above for the barycentric time correction.
The formulae used by GCOMBvcs::tdelta are identical to those
described above, except for the Shapiro time delay for which
the displacement of the Sun from the Solar System barycentre
is neglected:

$$
\Delta t_{\text{Shapiro}} = -2t_{\odot} \ln \left( 1 + \frac{\hat{\mathbf{n}} \cdot (\mathbf{e} + \mathbf{r})}{|\mathbf{e} + \mathbf{r}|} \right). \qquad (\text{D.10})
$$

Following the time correction, the pulsar phase $\Phi$ is
computed using the `GPulsarEphemeris::phase` method that
implements

$$
\Phi = \Phi_0 + \nu\Delta t + \frac{1}{2}\dot{\nu}\Delta t^2 + \frac{1}{6}\ddot{\nu}\Delta t^3, \quad (\text{D.11})
$$

where $\nu$, $\dot{\nu}$ and $\ddot{\nu}$ is the pulsar frequency and its first two time derivatives, $\Delta t = t_{SSB} - t_0$ is the elapsed time since the reference time of the ephemeris, and $\Phi_0$ is the pulsar phase at the reference time. Only the fractional part of the pulsar phase is retained so that its value is within the interval $[0, 1)$.

Pulsar information, including specifically the ephemerides of
a pulsar, is handled by the `GPulsar` class that supports reading
of ephemeris information from various file formats. The format
most relevant to COMPTEL is the `psrtime` format, which is an
ASCII file format containing the radio pulsar database as main-
tained during the CGRO mission by the pulsar group at Princeton
and nowadays by the Jodrell Bank Centre for Astrophysics.<fnref n="12" />
Upon loading of a `psrtime` file, the `GPulsar::load_psrtime`
method computes from the data the pulsar phase at the reference
time using

$$
\Phi_0 = \nu\Delta t + \frac{1}{2}\dot{\nu}\Delta t^2 + \frac{1}{6}\ddot{\nu}\Delta t^3, \quad (\text{D.12})
$$

where

$$
\Delta t = (\text{MJD}_0 - |\text{MJD}_0|) \times 86400 + \Delta t_{\text{CGRO}\to\text{SSB}} - \Delta t_{\text{Shapiro}} + \\
\Delta t_{\text{UTC}\to\text{TT}} + \Delta t_{\text{TT}\to\text{TBD}} \tag{D.13}
$$

and MJD₀ is the modified Julian day of the ephemeris reference
(the GPhemerides::geo2ssb method is used for this compu-
tation, and hence the Shapiro delay in Eq. (D.13) includes the
displacement of the Sun around the Solar System barycentre,
cf. Eq. (D.3)).

Appendix E: Handling of failed photomultiplier tubes

COMPTEL comprised two detector layers, composed of 7 and 14 circular modules for D1 and D2, respectively. Modules of the first layer were composed of the liquid scintillator NE 213A while modules of the second layer were made of NaI(Tl) scintillator crystals. Each D1 module was viewed by eight PMTs while each D2 module was viewed by seven PMTs. The relative amplitudes of the PMT signals for a given module allowed localisation of the interactions within the module with an average $1\sigma$ spatial resolution of 2.3 cm for D1 modules and 1.5 cm for D2 modules (Schönfelder et al. 1993).

During the operations of COMPTEL a certain number of D2
module PMTs failed (cf. Table E.1), degrading the interaction
localisation capabilities within the concerned modules and hence
increasing the uncertainties in the determination of the event
scatter directions (χ, ψ). GammaLib implements several options
for handling data from D2 modules with failed PMTs, controlled
through the `GCOMSelection::fpmtflag` method that takes an
integer value of 0, 1, or 2. This integer value is a user parameter
of the `ctools` scripts `comobsbin` and `compulbin` that, by default,
is assumed to be 0.

For `fpmtflag = 0`, events registered in D2 modules with failed PMTs are ignored, and the corresponding modules are also excluded in the computation of the geometry function (cf. Appendix F). Conversely, for `fpmtflag = 1` the failure of the PMTs is ignored, and events from D2 modules are treated as if the failed PMTs were still operating. Finally, for `fpmtflag = 2`, circular exclusion regions are applied around the zones of the failed PMTs for dates after their failure, as defined in Table E.1. Events localised within these regions are ignored and the regions are removed in the computation of the geometry function (cf. Appendix F). Using circular exclusion regions was the default for most of the COMPTEL analysis published in the literature in the past. The impact of the `fpmtflag` value on the analysis results is illustrated for the case of the spectral energy distribution of LS 5039 in Sect. 4.3.

<table><caption>Table E.1. D2 module exclusion circles applied following the failure of PMTs at the specified dates.</caption>
  <thead>
    <tr>
      <th>TJD</th>
      <th>VP</th>
      <th>Date</th>
      <th>Module</th>
      <th>x<sub>e</sub></th>
      <th>y<sub>e</sub></th>
      <th>r<sub>e</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>8718</td>
      <td>24.0</td>
      <td>6-4-1992</td>
      <td>13</td>
      <td>9.0</td>
      <td>41.2</td>
      <td>9.0</td>
    </tr>
    <tr>
      <td>8737</td>
      <td>26.0</td>
      <td>25-4-1992</td>
      <td>11</td>
      <td>-51.7</td>
      <td>8.7</td>
      <td>9.0</td>
    </tr>
    <tr>
      <td>8756</td>
      <td>28.0</td>
      <td>14-5-1992</td>
      <td>14</td>
      <td>-34.7</td>
      <td>49.0</td>
      <td>9.0</td>
    </tr>
    <tr>
      <td>8981</td>
      <td>204.0</td>
      <td>25-12-1992</td>
      <td>2</td>
      <td colspan="3">entire module</td>
    </tr>
  </tbody>
</table>

Notes. Module numbers are counted from 1. The circle centre locations *x* and *y* and the circle radii *r* are given in units of cm.

**Appendix F: Computation of the geometry function**

The geometry function DRG is computed in GammaLib by the
method `GCOMDri::compute_drg` using

$$
\mathrm{DRG}(\chi, \psi, \bar{\varphi}) = \frac{1}{N} \sum_{i \in \{S\}} \tilde{G}_i(\chi, \psi, \bar{\varphi}) \quad (\mathrm{F.1})
$$

<footnote n="12">¹² https://www.jodrellbank.manchester.ac.uk</footnote>
with

$$ \tilde{G}_i(\chi, \psi, \bar{\phi}) = \begin{cases} G_i(\chi, \psi), & \text{if EHA}(\chi, \psi) \ge \text{EHA}_{\min}(\bar{\phi}) \\ 0, & \text{otherwise} \end{cases}, \quad (F.2) $$

where EHA(χ, ψ) is the distance between scatter direction (χ, ψ) and the Earth horizon, EHA_min(φ) is given by Eq. (2), and N is the number of selected superpackets.

The G_i(χ, ψ) is the geometry factor for a given superpacket and corresponds to the area of the shadow that is cast by all active D1 modules on all active D2 modules for a given scatter angle (χ, ψ) divided by the total area of all D1 modules. The computation of the geometry factor is implemented in GCOMDri::compute_geometry and calculated using

$$ G_i(\chi, \psi) = \frac{1}{7} \sum_{k=0}^{6} \delta_k^{\text{D1}} \sum_{l=0}^{13} \delta_l^{\text{D2}} \tilde{\sigma}_{kl}(\theta, \phi), \quad (F.3) $$

where

$$ \delta_k^{\text{D1}} = \begin{cases} 1, & \text{D1 module } k \text{ active} \\ 0, & \text{otherwise} \end{cases}, \quad (F.4) $$

and

$$ \delta_l^{\text{D2}} = \begin{cases} 1, & \text{D2 module } l \text{ active} \\ 0, & \text{otherwise} \end{cases}, \quad (F.5) $$

are determined using the GCOMStatus class. In case of fpmtflag = 0, failed modules according to Table E.1 are also considered as inactive for superpacket dates after the dates of PMT failure. Furthermore,

$$ \tilde{\sigma}_{kl}(\theta, \phi) = \begin{cases} 0, & \text{if } d_{kl} \ge r_1 + r_2 \\ 1 - f_{kle}(\theta, \phi), & \text{if } d_{kl} \le r_2 - r_1 + 0.1, \\ o_{kl}(\theta, \phi) - f_{kle}(\theta, \phi), & \text{otherwise} \end{cases}, \quad (F.6) $$

with $r_1 = 13.8$ cm being the radius of a D1 module and $r_2 = 14.085$ cm being the radius of a D2 module. We note the margin of 0.1 cm in Eq. (F.6), which assures numerical stability with respect to rounding errors. The $d_{kl}$ is the projected distance between the centres of the D1 and D2 modules, given by

$$ d_{kl} = \sqrt{(x_l - x_k + h \tan \theta \cos \phi)^2 + (y_l - y_k + h \tan \theta \sin \phi)^2}, \quad (F.7) $$

where $x_k$ and $y_k$ are the geometric positions of the D1 modules and $x_l$ and $y_l$ the positions of the D2 modules with respect to the optical axis in cm, $h = 158$ cm being the vertical separation between D1 and D2 modules, $(\theta, \phi)$ being the zenith and azimuth angles of the Compton scatter direction (χ, ψ) with respect to the COMPTEL pointing direction, and

$$ o_{kl}(\theta, \phi) = \frac{r_1^2(\alpha_{kl} - \sin \alpha_{kl} \cos \alpha_{kl}) + r_2^2(\beta_{kl} - \sin \beta_{kl} \cos \beta_{kl})}{\pi r_1^2} \quad (F.8) $$

being the projected overlap of a D1 module and a D2 module with

$$ \cos \alpha_{kl} = \frac{d_{kl}^2 + (r_1^2 - r_2^2)}{2d_{kl}r_1} \quad (F.9) $$

and

$$ \cos \beta_{kl} = \frac{d_{kl}^2 - (r_1^2 - r_2^2)}{2d_{kl}r_2}. \quad (F.10) $$

The term $f_{kle}(\theta, \phi)$ accounts for the exclusion of circular regions around failed PMTs and differs from zero only for fpmtflag = 2. It quantifies the fractional overlap between the projected D1 module k and the part of the exclusion region for failed PMTs that is contained within D2 module l. The exclusion region is circular, and specified by a geometric centre position $x_e$ and $y_e$ and a radius $r_e$ as given in Table E.1. Specifically, for fpmtflag = 2

$$ f_{kle}(\theta, \phi) = \begin{cases} 0, & \text{if } d_{kl} > r_1 + r_2 \text{ or } d_{ke} < r_1 + r_e \\ f_{le}^{\text{contained}}(\theta, \phi), & \text{if } d_{ke} + r_e < r_1 \\ f_{kle}^{\text{partial}}(\theta, \phi), & \text{otherwise} \end{cases}, \quad (F.11) $$

with

$$ d_{ke} = \sqrt{(x_e - x_k + h \tan \theta \cos \phi)^2 + (y_e - y_k + h \tan \theta \sin \phi)^2} \quad (F.12) $$

being the distance between the projected D1 module centres $x_k$ and $y_k$ and the centre $x_e$ and $y_e$ of the exclusion region.

The $f_{le}^{\text{contained}}(\theta, \phi)$ is the overlap for the case that the exclusion circle is fully contained in the projected D1 module circumference. In this case, the relevant overlap to take into account is the overlap between the exclusion circle and the D2 module, given by

$$ f_{le}^{\text{contained}}(\theta, \phi) = \begin{cases} r_2^2/r_1^2, & \text{if } r_1 > d_{le} + r_2 \\ 1, & \text{if } r_2 > d_{le} + r_1 \\ \tilde{f}_{le}^{\text{contained}}(\theta, \phi), & \text{otherwise} \end{cases}, \quad (F.13) $$

with

$$ d_{le} = \sqrt{(x_e - x_l)^2 + (y_e - y_l)^2} \quad (F.14) $$

being the distance between the centres of the exclusion circle and the D2 module, and

$$ \tilde{f}_{le}^{\text{contained}}(\theta, \phi) = \frac{r_1^2(\alpha_{le} - \sin \alpha_{le} \cos \alpha_{le}) + r_2^2(\beta_{le} - \sin \beta_{le} \cos \beta_{le})}{\pi r_1^2} \quad (F.15) $$

being the overlap of the exclusion region and the D2 module with

$$ \cos \alpha_{le} = \frac{d_{le}^2 + (r_1^2 - r_2^2)}{2d_{le}r_1} \quad (F.16) $$

and

$$ \cos \beta_{le} = \frac{d_{le}^2 - (r_1^2 - r_2^2)}{2d_{le}r_2} \quad (F.17) $$

The evaluation of Eq. (F.13) is implemented by the method GCOMDri::compute.

Finally, $f_{kle}^{\text{partial}}(\theta, \phi)$ specifies the fractional overlap between the projected D1 module, the D2 module and the exclusion circle. This quantity is evaluated numerically by testing a grid of 25 × 25 x and y positions around the exclusion circle. The numerical evaluation is implemented by the method GCOMDri::compute_overlap.
# Appendix G: Response computation

The following sections describe some details of the response computations implemented in GammaLib.

## Appendix G.1: Efficiency factors

The efficiency factor $P_{\text{eff}}(\varphi'_{\text{geo}}, E_{\gamma})$ in Eq. (6) is factorised according to

$$ P_{\text{eff}}(\varphi'_{\text{geo}}, E_{\gamma}) = P_{\text{Al}}(E_{\gamma}) \times P_{\text{V1}}(E_{\gamma}) \times P_{\text{D1}}(E_{\gamma}) \times P_{\text{C}}(E_{\gamma}) \times \\ P_{\text{MH}}(E_{\gamma}) \times P_{\text{SV}}(E_{\gamma}) \times \\ P_{\text{PSD}}(\hat{E}_{1}) \times P_{\text{A2}}(\hat{E}_{2}) \times P_{\text{V23}}(\hat{E}_{2}) \times P_{\text{D2}}(\hat{E}_{2}) \times \\ P_{\text{MS}}(\varphi'_{\text{geo}}, E_{\gamma}), \qquad (G.1) $$

with

$$ \hat{E}_2 = \frac{E_\gamma}{(1 - \cos \varphi'_{\text{geo}}) \frac{E_\gamma}{m_e c^2} + 1} \qquad (G.2) $$

being the energy of the photon entering the D2 module for a true scatter angle $\varphi'_{\text{geo}}$ and an incident photon energy of $E_{\gamma}$, and $\hat{E}_1 = E_{\gamma} - \hat{E}_2$ being the true energy deposit in D1 for a single Compton scattering by an angle of $\varphi'_{\text{geo}}$. The GammaLib methods that implement the computation of the efficiency factors are summarised in Table G.1.

<table><caption>**Table G.1.** GammaLib methods that implement the computation of the efficiency factors.</caption><thead><tr><th>Factor</th><th>Method</th></tr></thead><tbody><tr><td>$P_{\text{eff}}(\varphi'_{\text{geo}}, E_{\gamma})$</td><td>GCOMIaq::weight_iaq</td></tr><tr><td>$P_{\text{Al}}(E_{\gamma})$</td><td>GCOMInstChars::trans_D1</td></tr><tr><td>$P_{\text{V1}}(E_{\gamma})$</td><td>GCOMInstChars::trans_V1</td></tr><tr><td>$P_{\text{D1}}(E_{\gamma})$</td><td>GCOMInstChars::prob_D1inter</td></tr><tr><td>$P_{\text{C}}(E_{\gamma})$</td><td>GCOMIaq::weight_iaq</td></tr><tr><td>$P_{\text{MH}}(E_{\gamma})$</td><td>GCOMInstChars::prob_no_multihit</td></tr><tr><td>$P_{\text{SV}}(E_{\gamma})$</td><td>GCOMInstChars::prob_no AMDVTO</td></tr><tr><td>$P_{\text{PSD}}(\hat{E}_1)$</td><td>GCOMInstChars::psdCorrection</td></tr><tr><td>$P_{\text{A2}}(\hat{E}_2)$</td><td>GCOMInstChars::trans_D2</td></tr><tr><td>$P_{\text{V23}}(\hat{E}_2)$</td><td>GCOMInstChars::trans_V23</td></tr><tr><td>$P_{\text{D2}}(\hat{E}_2)$</td><td>GCOMInstChars::prob_D2inter</td></tr><tr><td>$P_{\text{MS}}(\varphi'_{\text{geo}}, E_{\gamma})$</td><td>GCOMInstChars::multi_scatter</td></tr></tbody></table>

The $P_{\text{Al}}(E_{\gamma})$ is the transmission probability for photons for the material above D1, which is composed essentially of aluminium, and is computed using

$$ P_{\text{Al}}(E_{\gamma}) = e^{-\mu_{\text{Al}}(E_{\gamma}) l_{\text{above}}}, \qquad (G.3) $$

where $\mu_{\text{Al}}(E_{\gamma})$ is the energy-dependent interaction coefficient of aluminium in units of cm⁻¹ that is interpolated using a log-log interpolation of the values given in Table G.2 and $l_{\text{above}} = 0.147$ cm is the thickness of the material above D1.

The $P_{\text{V1}}(E_{\gamma})$ is the transmission probability for photons for the first Veto dome and is computed using

$$ P_{\text{V1}}(E_{\gamma}) = e^{-\mu_{\text{Veto}}(E_{\gamma}) l_{\text{V1}}}, \qquad (G.4) $$

where $\mu_{\text{Veto}}(E_{\gamma})$ is the energy-dependent interaction coefficient for the Veto dome in units of cm⁻¹ that is interpolated using a

log-log interpolation of the values given in Table G.2 and $l_{\text{V1}} = 1.721$ cm is the thickness of the first Veto dome.

$P_{\text{D1}}(E_{\gamma})$ is the interaction probability in D1 and is computed using

$$ P_{\text{D1}}(E_{\gamma}) = 1 - e^{-\mu_{\text{D1}}(E_{\gamma}) l_{\text{D1}}}, \qquad (G.5) $$

where $\mu_{\text{D1}}(E_{\gamma})$ is the energy-dependent D1 attenuation coefficient in units of cm⁻¹ that is interpolated using a log-log interpolation of the values given in Table G.2 and $l_{\text{D1}} = 8.5$ cm is the thickness of the D1 modules.

$P_{\text{C}}(E_{\gamma})$ is the fraction of Compton scatter interactions among all possible photon interactions within D1 and is given by

$$ P_{\text{Compton}}(E_{\gamma}) = \begin{cases} 1.067 - 0.0295E_{\gamma} + 3.4 \times 10^{-4}E_{\gamma}^2, & \text{if } E_{\gamma} \ge 2 \text{ MeV} \\ 1, & \text{otherwise} \end{cases}, \qquad (G.6) $$

with the additional constraint of $0 \le P_{\text{Compton}}(E_{\gamma}) \le 1$.

The $P_{\text{MH}}(E_{\gamma})$ is the probability that there is no multi-hit. This probability is computed using a log-log interpolation of the values of Table G.2. $P_{\text{SV}}(E_{\gamma})$ is the probability that the photon was not self-vetoed. This probability is computed using a linear interpolation of the values of Table G.2.

The $P_{\text{PSD}}(\hat{E}_1)$ is the D1 energy-dependent PSD correction and is given by

$$ P_{\text{PSD}}(\hat{E}_1) = 1 - \frac{1}{1727.9 \times \hat{E}_1^{2.53} + 1}, \qquad (G.7) $$

where $\hat{E}_1$ is in units of MeV. This correction applies to a standard PSD selection of 0–110.

$P_{\text{A2}}(\hat{E}_2)$ gives the transmission probability of the aluminium below the D1 modules and is computed using

$$ P_{\text{A2}}(\hat{E}_2) = e^{-\mu_{\text{Al}}(\hat{E}_2) l_{\text{T23}}}, \qquad (G.8) $$

where $\mu_{\text{Al}}(\hat{E}_2)$ is the energy-dependent interaction coefficient for aluminium in units of cm⁻¹ that is interpolated using a log-log interpolation of the values given in Table G.2 and $l_{\text{between}} = 0.574$ cm is the thickness of the aluminium plate.

$P_{\text{V23}}(\hat{E}_2)$ gives the transmission probability of the second and third veto domes that are situated between the D1 and D2 modules, and is computed using

$$ P_{\text{V23}}(\hat{E}_2) = e^{-\mu_{\text{Veto}}(\hat{E}_2) l_{\text{V23}}}, \qquad (G.9) $$

where $\mu_{\text{Veto}}(\hat{E}_2)$ is the energy-dependent interaction coefficient of the Veto domes in units of cm⁻¹ that is interpolated using a log-log interpolation of the values given in Table G.2 and $l_{\text{V23}} = 3.442$ cm is the combined thickness of the second and third veto domes.

$P_{\text{D2}}(E_{\gamma})$ gives the interaction probability in D2 and is computed using

$$ P_{\text{D2}}(E_{\gamma}) = 1 - e^{-\mu_{\text{D2}}(\hat{E}_2) l_{\text{D2}}}, \qquad (G.10) $$

where $\mu_{\text{D2}}(\hat{E}_2)$ is the energy-dependent D2 attenuation coefficient in units of cm⁻¹ that is interpolated using a log-log interpolation of the values given in Table G.2 and $l_{\text{D2}} = 7.525$ cm is the thickness of the D2 modules.

The $P_{\text{MS}}(E_{\gamma}, \hat{E}_2)$ is the probability that a photon that has interacted in D1 leaves the D1 module without any further interaction. This probability is computed using

$$ P_{\text{MS}}(\varphi'_{\text{geo}}, E_{\gamma}) = \frac{\int_{0}^{r_{\text{D1}}} \int_{0}^{l_{\text{D1}}} \mu_{\text{D1}}(E_{\gamma}) e^{-\mu_{\text{D1}}(E_{\gamma}) z} \int_{0}^{\pi} e^{-\mu_{\text{D1}}(\hat{E}_2) l(r,z,\phi)} d\phi dz dr}{\pi r_{\text{D1}}^2 (1 - e^{-\mu_{\text{D1}}(E_{\gamma}) l_{\text{D1}}})}, $$
<table><caption>**Table G.2.** Efficiency parameters as a function of energy.</caption><thead><tr><th>E (MeV)</th><th>μ<sub>Veto</sub>(E)</th><th>μ<sub>Al</sub>(E)</th><th>μ<sub>D1</sub>(E)</th><th>μ<sub>D2</sub>(E)</th><th>P<sub>MS</sub>(E)</th><th>P<sub>SV</sub>(E)</th></tr></thead><tbody><tr><td>0.05</td><td>0.205944</td><td></td><td>0.178262</td><td></td><td></td><td></td></tr><tr><td>0.1</td><td>0.169298</td><td>0.432</td><td>0.146770</td><td>5.9821</td><td></td><td></td></tr><tr><td>0.15</td><td>0.150704</td><td>0.3618</td><td>0.130690</td><td>2.15796</td><td></td><td></td></tr><tr><td>0.2</td><td>0.137444</td><td>0.324</td><td>0.119207</td><td>1.14871</td><td>1.0</td><td></td></tr><tr><td>0.25</td><td>0.127823</td><td></td><td>0.110860</td><td></td><td></td><td></td></tr><tr><td>0.3</td><td>0.119531</td><td>0.2781</td><td>0.103667</td><td>0.57986</td><td>1.0</td><td></td></tr><tr><td>0.35</td><td>0.112835</td><td></td><td>0.097861</td><td></td><td></td><td></td></tr><tr><td>0.4</td><td>0.106932</td><td>0.24894</td><td>0.092743</td><td>0.41471</td><td>1.0</td><td></td></tr><tr><td>0.45</td><td>0.101949</td><td></td><td>0.088426</td><td></td><td></td><td></td></tr><tr><td>0.5</td><td>0.097504</td><td>0.22707</td><td>0.084574</td><td>0.336906</td><td>1.0</td><td></td></tr><tr><td>0.55</td><td>0.093717</td><td></td><td>0.081289</td><td></td><td></td><td></td></tr><tr><td>0.6</td><td>0.090226</td><td>0.21006</td><td>0.078260</td><td>0.293233</td><td>1.0</td><td></td></tr><tr><td>0.65</td><td>0.087195</td><td></td><td>0.075629</td><td></td><td></td><td></td></tr><tr><td>0.7</td><td>0.084350</td><td></td><td>0.073160</td><td></td><td></td><td></td></tr><tr><td>0.75</td><td>0.081692</td><td></td><td>0.070853</td><td></td><td></td><td></td></tr><tr><td>0.8</td><td>0.079220</td><td>0.18468</td><td>0.068708</td><td>0.242954</td><td>1.0</td><td></td></tr><tr><td>0.85</td><td>0.077065</td><td></td><td>0.066840</td><td></td><td></td><td></td></tr><tr><td>0.9</td><td>0.075009</td><td></td><td>0.065057</td><td></td><td></td><td></td></tr><tr><td>0.95</td><td>0.073052</td><td></td><td>0.063360</td><td></td><td></td><td></td></tr><tr><td>1.0</td><td>0.071194</td><td>0.16578</td><td>0.061749</td><td>0.212493</td><td>0.993</td><td></td></tr><tr><td>1.25</td><td>0.064013</td><td></td><td>0.055524</td><td></td><td></td><td></td></tr><tr><td>1.5</td><td>0.058035</td><td></td><td>0.050340</td><td>0.170655</td><td>0.991</td><td></td></tr><tr><td>1.75</td><td>0.053580</td><td></td><td>0.046471</td><td></td><td></td><td></td></tr><tr><td>2.0</td><td>0.049689</td><td>0.1166</td><td>0.043092</td><td>0.150837</td><td>0.985</td><td></td></tr><tr><td>2.5</td><td>0.044222</td><td></td><td>0.038345</td><td></td><td></td><td></td></tr><tr><td>3.0</td><td>0.039767</td><td>0.0953</td><td>0.034476</td><td>0.134322</td><td>0.976</td><td></td></tr><tr><td>3.5</td><td>0.036580</td><td></td><td>0.031708</td><td></td><td></td><td></td></tr><tr><td>4.0</td><td>0.033893</td><td>0.0837</td><td>0.029373</td><td>0.128817</td><td>0.975</td><td></td></tr><tr><td>4.4</td><td></td><td></td><td></td><td></td><td></td><td>0.940</td></tr><tr><td>4.5</td><td>0.031814</td><td></td><td>0.027565</td><td></td><td></td><td></td></tr><tr><td>5.0</td><td>0.030016</td><td>0.07641</td><td>0.026002</td><td>0.127349</td><td>0.973</td><td></td></tr><tr><td>5.5</td><td>0.028562</td><td></td><td>0.024737</td><td></td><td></td><td></td></tr><tr><td>6.0</td><td>0.027268</td><td>0.07155</td><td>0.023611</td><td>0.127716</td><td>0.971</td><td></td></tr><tr><td>6.1</td><td></td><td></td><td></td><td></td><td></td><td>0.912</td></tr><tr><td>6.5</td><td>0.026254</td><td></td><td>0.022727</td><td></td><td></td><td></td></tr><tr><td>7.0</td><td>0.025319</td><td></td><td>0.021912</td><td></td><td></td><td></td></tr><tr><td>7.5</td><td>0.024462</td><td></td><td>0.021166</td><td></td><td></td><td></td></tr><tr><td>8.0</td><td>0.023684</td><td>0.06534</td><td>0.020488</td><td>0.130285</td><td>0.950</td><td></td></tr><tr><td>8.5</td><td>0.023043</td><td></td><td>0.019929</td><td></td><td></td><td></td></tr><tr><td>9.0</td><td>0.022442</td><td></td><td>0.019405</td><td></td><td></td><td></td></tr><tr><td>9.5</td><td>0.021880</td><td></td><td>0.018915</td><td></td><td></td><td></td></tr><tr><td>10.0</td><td>0.021357</td><td>0.0621</td><td>0.018459</td><td>0.134689</td><td>0.929</td><td></td></tr><tr><td>12.0</td><td></td><td></td><td></td><td></td><td></td><td>0.850</td></tr><tr><td>12.5</td><td>0.019618</td><td></td><td>0.016940</td><td></td><td></td><td></td></tr><tr><td>15.0</td><td>0.018280</td><td>0.05859</td><td>0.015770</td><td>0.147534</td><td>0.912</td><td></td></tr><tr><td>17.5</td><td>0.017477</td><td></td><td>0.015066</td><td></td><td></td><td></td></tr><tr><td>20.0</td><td>0.016805</td><td>0.05805</td><td>0.014475</td><td>0.158544</td><td>0.913</td><td></td></tr><tr><td>20.5</td><td></td><td></td><td></td><td></td><td></td><td>0.718</td></tr><tr><td>25.0</td><td>0.015988</td><td></td><td>0.013756</td><td></td><td>0.915</td><td></td></tr><tr><td>30.0</td><td>0.015407</td><td>0.05859</td><td>0.013242</td><td>0.176894</td><td>0.917</td><td></td></tr><tr><td>35.0</td><td>0.015158</td><td></td><td>0.013019</td><td></td><td></td><td></td></tr><tr><td>40.0</td><td>0.014948</td><td>0.06048</td><td>0.012830</td><td>0.190840</td><td></td><td></td></tr><tr><td>45.0</td><td>0.014756</td><td></td><td>0.012659</td><td></td><td></td><td></td></tr><tr><td>50.0</td><td>0.014646</td><td>0.06210</td><td>0.012558</td><td>0.201116</td><td></td><td></td></tr></tbody></table>

**Notes.** Interaction coefficients, $\mu$, are given in units of cm$^{-1}$.

(G.11)

$r_{D1}$ = 13.8 cm and $l_{D1}$ = 8.5 cm is the radius and thickness of a D1 module. $l(r, z, \phi)$ gives the remaining interaction length for a photon that has interacted at radius $r$, vertical depth $z$ and

which integrates all possible second interaction locations within a given D1 module assuming vertically incident photons, where
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Fig. G.1. Energy dependence of the efficiency factor $P_{\text{eff}}(\phi'_{\text{geo}}, E_{\gamma})$ for the three geometrical scatter angles $\phi'_{\text{geo}} = 5^{\circ}, 25^{\circ}$, and $50^{\circ}$ (left) and the energy dependence of its components for $\phi'_{\text{geo}} = 25^{\circ}$ (right).</figcaption></figure>

azimuth angle $\phi$ of the module and is given by

$$l(r, z, \phi) = \begin{cases} \frac{l_{\text{D1}} - z}{\cos \phi'_{\text{geo}}}, & \text{if } l(r, z, \phi) \sin \phi'_{\text{geo}} \le R(r, \phi) \\ \frac{R(r, \phi)}{\sin \phi'_{\text{geo}}}, & \text{otherwise} \end{cases} \qquad (\text{G.12})$$

with

$$R(r, \phi) = \sqrt{r_{\text{D1}}^2 - r^2 \sin^2 \phi - r \cos \phi}. \qquad (\text{G.13})$$

Figure G.1 illustrates the energy dependence of the efficiency factor $P_{\text{eff}}(\phi'_{\text{geo}}, E_{\gamma})$ and its components.

## Appendix G.2: Klein-Nishina contribution

$$P_{\text{KN}}(\phi'_{\text{geo}}, E_{\gamma}) = \frac{\int_{\phi'_{\text{geo,min}}}^{\phi'_{\text{geo,max}}} \sigma_{\text{KN}}(\phi''_{\text{geo}}, E_{\gamma}) d\phi''_{\text{geo}}}{\int_{0}^{\pi} \sigma_{\text{KN}}(\phi''_{\text{geo}}, E_{\gamma}) d\phi''_{\text{geo}}} \qquad (\text{G.14})$$

is the contribution of the Klein-Nishina cross-section $\sigma_{\text{KN}}(\phi''_{\text{geo}}, E_{\gamma})$ to the $\phi'_{\text{geo}}$ bin spanned by $[\phi'_{\text{geo,min}}, \phi'_{\text{geo,max}}]$, where the integrals are computed using

$$\int_{\phi'_{\text{geo,min}}}^{\phi'_{\text{geo,max}}} \sigma_{\text{KN}}(\phi''_{\text{geo}}, E_{\gamma}) d\phi''_{\text{geo}} = \int_{0}^{E_{1}^{\text{up}}} \sigma_{\text{KN}}(E'_{1}, E_{\gamma}) dE'_{1} - \int_{0}^{E_{1}^{\text{low}}} \sigma_{\text{KN}}(E'_{1}, E_{\gamma}) dE'_{1} \qquad (\text{G.15})$$

with

$$\int_{0}^{E_{1}} \sigma_{\text{KN}}(E'_{1}, E_{\gamma}) dE'_{1} = \frac{E_{1}}{E_{\gamma}} \left( \frac{m_{e}c^{2}}{E_{1}} + 1 \right)^{2} + \left( 2 \left( \frac{m_{e}c^{2}}{E_{\gamma}} \right)^{2} + 2 \frac{m_{e}c^{2}}{E_{\gamma}} - 1 \right) \times \log \left( 1 - \frac{E_{1}}{E_{\gamma}} \right) - \frac{1}{2} \left( \frac{E_{1}}{E_{\gamma}} \right)^{2} + \left( \frac{m_{e}c^{2}}{E_{\gamma}} \right)^{2} \frac{1}{1 - \frac{E_{1}}{E_{\gamma}}}. \qquad (\text{G.16})$$

## Appendix G.3: Module spectral response functions

### Appendix G.3.1: D1 spectral response

The spectral response $R_{\text{D1}}(E_1|\hat{E}_1)$ provides the probability for measuring an energy $E_1$ when an incident photon with energy $E_\gamma$ was Compton scattered by an angle $\phi_{\text{geo}}$, producing a true energy deposit of $\hat{E}_1$ in the detector module. The D1 response is implemented by the GCOMD1Response class and the response evaluation is done using the GCOMD1Response::operator() operator. The spectral D1 response is computed using a Gaussian

$$R_{\text{D1}}(E_1|\hat{E}_1) = \frac{T_1(E_1|\hat{E}_1)}{\sigma(\hat{E}_1)\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\mu_p(\hat{E}_1) - E_1}{\sigma(\hat{E}_1)}\right)^2\right), \qquad (\text{G.17})$$

where $\mu_p(\hat{E}_1)$ is the position of the Gaussian peak for true energy deposit of $\hat{E}_1$, with $\mu_p(\hat{E}_1) \approx \hat{E}_1$, and $\sigma(\hat{E}_1)$ is the standard deviation of the Gaussian (Diehl et al. 1992). $T_1(E_1|\hat{E}_1)$ is a threshold function, defined by

$$T_1(E_1|\hat{E}_1) = \begin{cases} 0, & \text{if } E_1 \le E_1^{\text{thres}}(\hat{E}_1) \\ \frac{E_1 - E_1^{\text{thres}}(\hat{E}_1)}{\Delta E_1(\hat{E}_1)}, & \text{if } 0 < E_1 - E_1^{\text{thres}}(\hat{E}_1) \le \Delta E_1(\hat{E}_1) \\ 0, & \text{if } E_1 > E_1^{\text{max}}(\hat{E}_1) \\ 1, & \text{otherwise} \end{cases}, \qquad (\text{G.18})$$

where $E_1^{\text{thres}}(\hat{E}_1) = E_1^{\text{min}}(\hat{E}_1) - \frac{1}{2}\Delta E_1(\hat{E}_1)$. $\mu_p(\hat{E}_1)$ and $\sigma(\hat{E}_1)$ as well as $E_1^{\text{min}}(\hat{E}_1)$, $E_1^{\text{max}}(\hat{E}_1)$ and $\Delta E_1(\hat{E}_1)$ are tabulated as a function of $\hat{E}_1$ in a calibration file of type SDA that is included in GammaLib. Values for a given energy $\hat{E}_1$ are obtained by linear interpolation of the tabulated values.

The spectral response $R_{\text{D1}}(E_1|\hat{E}_1)$ is illustrated in the left panel of Fig. G.2 for a selected number of input energies $\hat{E}_1$. The figure can be compared with Fig. 3c in Diehl et al. (1992), which presents the response functions as implemented in COMPASS at the time of the CGRO launch for comparable energies.
<figure><img src="image_24.png" /><figcaption>Fig. G.2. Detector response functions $R_{D1}(E_1|\hat{E}_1)$ (left) and $R_{D2}(E_2|\hat{E}_2)$ (right) for a selected number of energies that have been chosen to allow for a comparison with Figs. 3c and 4b in Diehl et al. (1992).</figcaption></figure>

<figure><img src="image_25.png" /></figure>

## Appendix G.3.2: D2 spectral response

The spectral response $R_{D2}(E_2|\hat{E}_2)$ provides the probability for measuring an energy $E_2$ when a photon with energy $\hat{E}_2$ hit a D2 detector module. The D2 response is implemented by the GCOMD2Response class and the response evaluation is done using the GCOMD2Response::operator() operator. The spectral D2 response is computed using

$$R_{D2}(E_2|\hat{E}_2) = (P_{pp}(E_2|\hat{E}_2) + P_{se}(E_2|\hat{E}_2) + P_{de}(E_2|\hat{E}_2) + C_{com}(E_2|\hat{E}_2) + C_{bkg}(E_2|\hat{E}_2)) \times T_2(E_2|\hat{E}_2), \quad (G.19)$$

which is a combination of three Gaussian and two continuum components multiplied by a threshold function $T_2(E_2|\hat{E}_2)$ (Diehl et al. 1992). Here, $\hat{E}_2$ is the energy of the gamma-ray photon incident on a D2 module, and $E_2$ is the measured energy deposit.

The first Gaussian,

$$P_{pp}(E_2|\hat{E}_2) = \frac{A_{pp}(\hat{E}_2)}{\sigma(\hat{E}_2)\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\mu_{pp}(\hat{E}_2) - E_2}{\sigma(\hat{E}_2)}\right)^2\right), \quad (G.20)$$

describes the so-called photo-peak which represents photons that were fully absorbed in the D2 module. We note that the position $\mu_{pp}(\hat{E}_2)$ of the photo-peak is not necessarily identical to the incident energy, although below ~ 10 MeV $\mu_{pp}(\hat{E}_2) \approx \hat{E}_2$ is well satisfied.

The two subsequent Gaussians describe the single and double escape peak, respectively, which are generated by photons above 1.022 MeV that undergo pair creation resulting in the production of two 511 keV gamma rays. The escape of one or both of these 511 keV photons from the D2 module leads to characteristic peaks at energies 0.511 MeV and 1.022 MeV below the photo-peak energy. The escape peaks are modelled using

$$P_{se}(E_2|\hat{E}_2) = \frac{A_{se}(\hat{E}_2)}{\sigma(\mu_{se}(\hat{E}_2))\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\mu_{se}(\hat{E}_2) - E_2}{\sigma(\mu_{se}(\hat{E}_2))}\right)^2\right) \quad (G.21)$$

and

$$P_{de}(E_2|\hat{E}_2) = \frac{A_{de}(\hat{E}_2)}{\sigma(\mu_{de}(\hat{E}_2))\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\mu_{de}(\hat{E}_2) - E_2}{\sigma(\mu_{de}(\hat{E}_2))}\right)^2\right), \quad (G.22)$$

with $\mu_{se}(\hat{E}_2) = \mu_{pp}(\hat{E}_2) - m_e c^2$ and $\mu_{de}(\hat{E}_2) = \mu_{pp}(\hat{E}_2) - 2m_e c^2$, where $m_e c^2 \approx 511$ keV is the rest energy of the electron.

The first continuum component,

$$C_{com}(E_2|\hat{E}_2) = A_{com}(\hat{E}_2) \times \int_{E_2 - 3\sigma(E_2)}^{\max(E_2 + 3\sigma(E_2), E_c)} \sigma_{KN}(E|\hat{E}_2) f(E|E_2, \hat{E}_2) K(E|E_2) dE' \quad (G.23)$$

describes the so-called Compton tail which represents photons that underwent a Compton scattering before escaping the D2 module. The component is the product of the Klein-Nishina cross-section for Compton scattering,

$$\sigma_{KN}(E|\hat{E}_2) = \left( \frac{E/\hat{E}_2}{1 - E/\hat{E}_2} \frac{m_e c^2}{\hat{E}_2} \right)^2 - \frac{E}{\hat{E}_2} + \frac{1}{1 - E/\hat{E}_2}, \quad (G.24)$$

and the probability that the Compton scattered photon does not undergo a second Compton scattering before escaping, empirically modelled using

$$f(E|\hat{E}_2) = \begin{cases} \exp(-\mu(E|\hat{E}_2)l(\hat{E}_2)), & \text{if } \hat{E}_2 > 12.14 \text{ MeV} \\ 0, & \text{otherwise} \end{cases}, \quad (G.25)$$

where

$$\mu(E|\hat{E}_2) = 0.72 e^{-1.28(\hat{E}_2 - E)^{0.35}} + 0.01(\hat{E}_2 - E) + 0.014(\hat{E}_2 - E)^{-2.5} \quad (G.26)$$

is the total linear attenuation coefficient in NaI, which is the scintillator material of the D2 modules and

$$l(\hat{E}_2) = 2.9 \ln(\hat{E}_2 - 11.14) \quad (G.27)$$

is an empirical path length in the D2 module, where $\hat{E}_2$ is in units of MeV. The product is convolved with a Gaussian kernel,

$$K(E|E_2) = \frac{1}{\sigma(E_2)\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{E - E_2}{\sigma(E_2)}\right)^2\right), \quad (G.28)$$

to take the energy resolution of the D2 module into account, where $\sigma(E_2)$ is the standard deviation of the detector response at the reconstructed energy $E_2$.
The second continuum component,

$$C_{\text{bkg}}(E_2|\hat{E}_2) = A_{\text{bkg}}(\hat{E}_2) \int_{E_2-3\sigma(E_2)}^{\max(E_2+3\sigma(E_2), \hat{E}_2)} K(E|E_2) dE, \quad (\text{G.29})$$

is a flat background component that is convolved by a Gaussian kernel and that takes into account any higher-order scatterings in the D2 modules.

The threshold function is defined by

$$T_2(E_2|\hat{E}_2) = \begin{cases} \exp\left(-\frac{1}{2}\left(\frac{E_2^{\text{thres}}(\hat{E}_2) - E_2}{\sigma_{\text{thres}}(\hat{E}_2)}\right)^2\right), & \text{if } E_2 < E_2^{\text{thres}}(\hat{E}_2) \\ 0, & \text{if } E_2 > E_2^{\text{max}}(\hat{E}_2) \\ 1, & \text{otherwise} \end{cases}, \quad (\text{G.30})$$

where $E_2^{\text{thres}}(\hat{E}_2) = E_2^{\text{min}}(\hat{E}_2) + \frac{1}{2}\Delta E_2(\hat{E}_2)$ and $\sigma_{\text{thres}}(\hat{E}_2) = \Delta E_2(\hat{E}_2)/\sqrt{2\log 2}$

All five response components have amplitude parameters ($A_{\text{pp}}$, $A_{\text{se}}$, $A_{\text{de}}$, $A_{\text{com}}$ and $A_{\text{bkg}}$) that are tabulated as a function of incident energy $\hat{E}_2$ in a calibration file of type SDB. In addition, this file contains also tabulated values for the photo-peak energy $\mu_{\text{pp}}(\hat{E}_2)$, the standard deviation $\sigma(\hat{E}_2)$ and the minimum and maximum threshold energies $E_2^{\text{min}}(\hat{E}_2)$ and $E_2^{\text{max}}(\hat{E}_2)$ and width $\Delta E_2(\hat{E}_2)$. Values for a given energy $\hat{E}_2$ are obtained by linear interpolation of the tabulated values.

The spectral response $R_{\text{D2}}(E_2|\hat{E}_2)$ is illustrated in the right panel of Fig. G.2 for a selected number of input energies $\hat{E}_2$. The figure can be compared with Fig. 4b in Diehl et al. (1992), which presents the response functions as implemented in COMPASS at the time of the CGRO launch for comparable energies.

## Appendix G.4: Response convolution

The response is convolved with the Gaussian kernel

$$K(\varphi'_{\text{geo}}|\varphi_{\text{geo}}, \theta) = \frac{1}{\sigma(\varphi'_{\text{geo}}, \theta) \sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\varphi_{\text{geo}} - \varphi'_{\text{geo}}}{\sigma(\varphi'_{\text{geo}}, \theta)}\right)^2\right) \quad (\text{G.31})$$

to take into account the location uncertainties of the events in the D1 and D2 modules. The width $\sigma(\varphi'_{\text{geo}}, \theta)$ of the Gaussian kernel is a function of the scatter angle $\varphi'_{\text{geo}}$ and the zenith angle $\theta$ of the incoming photon, measured with respect to the COMPTEL pointing axis. The width

$$\sigma(\varphi'_{\text{geo}}\theta) = \sqrt{\sigma_{xy}^2(\varphi'_{\text{geo}}, \theta) + \sigma_z^2(\varphi'_{\text{geo}}, \theta)} \quad (\text{G.32})$$

is the squared mean of the uncertainty

$$\sigma_{xy}^2(\varphi'_{\text{geo}}, \theta) = \frac{(\sigma_x^2 + \sigma_y^2)}{h^2} f_{xy}(\varphi'_{\text{geo}}, \theta) \quad (\text{G.33})$$

in the horizontal direction and the uncertainty

$$\sigma_z^2(\varphi'_{\text{geo}}, \theta) = \frac{(\sigma_{z1}^2 + \sigma_{z2}^2)}{h^2} f_z(\varphi'_{\text{geo}}, \theta) \quad (\text{G.34})$$

in the vertical direction, where $\sigma_x = 2.3$ cm and $\sigma_y = 1.96$ cm is the average horizontal location spread in a D1 module with

<figure><img src="image_26.png" /><figcaption>Fig. G.3. Width $\sigma(\varphi'_{\text{geo}}, \theta)$ of the Gaussian smoothing kernel as a function of $\varphi'_{\text{geo}}$ for zenith angles $\theta$ of 0° (solid blue), 25° (dashed orange), 35° (dashed-dotted green), and 50° (dotted purple).</figcaption></figure>

respect to the D2 module in x and y direction, respectively, $\sigma_{z1} = 2.45$ cm and $\sigma_{z2} = 2.17$ cm is the vertical location spread in a D1 and D2 module, respectively, assuming a homogenous event distribution in the z-direction, and $h = 158$ cm is the distance between the D1 and D2 modules. The functions

$$f_{xy}(\varphi'_{\text{geo}}, \theta) = \cos^2 \theta \cos^2 \varphi'_{\text{geo}} (4 \sin^2 \theta \sin^2 \varphi'_{\text{geo}} + \cos^2 \varphi'_{\text{geo}} - \sin^2 \theta) + \frac{1}{2} \sin^2 \theta \times \left( \cos^2 \theta \cos^2 \varphi'_{\text{geo}} + \sin^2 \theta \sin^2 \varphi'_{\text{geo}} \left(1 - \frac{3}{4} \cos^2 \varphi'_{\text{geo}}\right) \right) \quad (\text{G.35})$$

and

$$f_z(\varphi'_{\text{geo}}, \theta) = \cos^4 \theta \sin^2 \varphi'_{\text{geo}} \cos^2 \varphi'_{\text{geo}} + \cos^2 \theta \sin^2 \theta \left(\frac{1}{2} - 3 \sin^2 \varphi'_{\text{geo}} \cos^2 \varphi'_{\text{geo}}\right) + \frac{3}{8} \sin^4 \theta \sin^2 \varphi'_{\text{geo}} \cos^2 \varphi'_{\text{geo}} \quad (\text{G.36})$$

are the geometrical relations that transform from the uncertainties within the modules to uncertainties in $\varphi'_{\text{geo}}$ for a given zenith angle $\theta$.

Since COMPTEL response functions have no explicit zenith angle dependence, the computations are performed for an average zenith angle of $\theta = 25^\circ$. The uncertainty introduced by this approximation is illustrated in Fig. G.3, and amounts to less than 10% for zenith angle $\theta < 35^\circ$.

## Appendix H: Carbon footprint of this research

Given the repeated alarms of the scientific community on global warming (IPCC 2021, 2022), biodiversity loss (IPBES 2019) and air, water and soil pollution (Landrigan et al. 2017), the environmental impact of any human activity, including obviously research, can no longer be ignored. Here we assess the GHG emissions related to our research work, including contributions from office use, research activities and research environment. We based our assessment on the carbon footprint estimate of the Institut de Recherche en astrophysique et planétologie (IRAP;
<table><caption>**Table H.1.** Carbon footprint estimate by emission source of the research work behind this paper.</caption><thead><tr><th>Source</th><th>Emissions kgCO<sub>2</sub>e</th><th>Comments</th></tr></thead><tbody><tr><td>Electricity to power office building</td><td>525 ± 53</td><td></td></tr><tr><td>Heating of office building</td><td>422 ± 118</td><td></td></tr><tr><td>Water usage in office building</td><td>8 ± 2</td><td></td></tr><tr><td>Air conditioning in office building</td><td>80 ± 23</td><td></td></tr><tr><td>Waste management</td><td>209 ± 106</td><td></td></tr><tr><td>Computing infrastructure</td><td>400 ± 200</td><td>Evaluated using csfootprint</td></tr><tr><td>Data flow</td><td>8 ± 8</td><td>Includes videoconferencing</td></tr><tr><td>Lunch meals</td><td>171 ± 103</td><td>10% classical, 10% flexitarian, 80% vegetarian</td></tr><tr><td>Home-to-office commuting</td><td>89 ± 36</td><td>90% cycling, 5% public transport, 5% car</td></tr><tr><td>Total</td><td>1911 ± 284</td><td></td></tr></tbody></table>

**Notes.** The footprint is based on a total work time estimate of one FTE or 220 working days.

Martin et al. 2022) where most of the work for this research was done. To attribute a share of the total footprint to our work, we estimated that roughly one full time equivalent (FTE) was necessary to accomplish this research, corresponding to 220 working days. Consequently, we use the annual emissions per source from Martin et al. (2022) divided by the IRAP staff of 263 persons to assess the contributions from the various sources of GHG emissions.

The results of our assessment are summarised in Table H.1. We estimate the total carbon footprint of our research work to be 1.9 ± 0.3 tCO₂e, composed of 1.2 ± 0.2 tCO₂e (65%) for running the office building in which the work was done, 0.4 ± 0.2 tCO₂e (21%) for the computing infrastructure and data flow, and 0.3 ± 0.1 tCO₂e (14%) for lunch meals and home-to-office commuting.

We based our assessment for computing on the carbon footprint estimation implemented in ctools, yet when we started this research the relevant functionality was not yet available. We note, however, that the carbon footprint estimate for the final analyses for this publication amounted to 60 kgCO₂e, of which 20 kgCO₂e were attributed to power consumption and 40 kgCO₂e to the computing infrastructure. Making the conservative estimate that ten analysis iterations were necessary for this research work we obtain a total carbon footprint due to computing of 600 kgCO₂e, of which 200 kgCO₂e were due to power consumption and 400 kgCO₂e due to the computing infrastructure. To avoid double counting, we disregarded the contribution of power consumption since it is already accounted for in the electricity to power the office building.

For emission from lunch meals and home-to-office commuting we took the specific habits of the co-authors weighted by their contribution to the research into account. We note that our research work did not require any professional travelling. The footprint for a couple of video conferences that we organised for this research is included under the 'Data flow' source.