COMPTEL data analysis using GammaLib and ctools

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 inMeV 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 beneﬁt from the scientiﬁc advances made in the ﬁeld 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 workﬂows 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 beneﬁt from recent advances in gamma-ray astronomy, and gives the community a means to unveil its still hidden treasures.


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 highenergy gamma rays from a plethora of objects and the discovery of new source classes, such as novae, gamma-ray binaries, starforming 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) 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 1 https://heasarc.gsfc.nasa.gov/docs/cgro/archive/ 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 maximumlikelihood 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 COMP-TEL data analysis tasks. These building blocks can be combined A&A 665, A84 (2022) 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 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 (χ, ψ) was obtained from the interaction locations in D1 and D2, the Compton scattering anglē ϕ was computed from the measured energy deposits E 1 and E 2 in D1 and D2, respectively, usinḡ 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 anticoincidence shield, composed of two veto domes for each layer, that allowed for the reduction of instrumental background 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.

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 COMPTELspecific C++ classes through the gammalib Python module. All COMPTEL-specific classes begin with the letters GCOM.

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 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 GCOMEventList, 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 ToF III values (see Appendix B).

Event cubes
COMPTEL data analysis is performed on binned event data using three-dimensional event cubes spanned by the Compton scatter direction (χ, ψ) and the Compton scattering angleφ. 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  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).
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 Gam-maLib (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 gammaray 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 ζ 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 ≥ EHA min (φ) with whereφ min is the lower boundary of the event cube layer that comprises the Compton scattering angleφ. By default, ζ min = 5 • . 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 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). 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. Only events for modules that are signalled as active are considered for the event binning. The detector module status is provided by the GCOMStatus 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 6 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.

Factorisation
The COMPTEL instrument response function, given in units of events cm 2 photons −1 , is factorised using where DRX(α, δ) is the exposure in units of cm 2 s towards a given celestial direction (α, δ), DRG(χ, ψ,φ) 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φ dependence), 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.
IAQ(φ|ϕ geo , E γ ) is the probability that a gamma-ray photon with energy E γ being Compton scattered in D1 interacts in D2, with ϕ geo = arccos (sin ψ sin δ + cos ψ cos δ cos(χ − α)) being the angular separation between (χ, ψ) and (α, δ), L is the livetime 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.

Exposure
The exposure map DRX(α, δ) is computed in GammaLib by the method GCOMDri::compute_drx using where T 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 θ i factor takes into account the reduction of the effective D1 surface when viewed from an off-axis direction θ 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 τ = 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 2 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.

Geometry function
The geometry function DRG(χ, ψ,φ) is the probability that a photon that was scattered in one of the D1 modules towards the  direction (χ, ψ) 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φ 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.

Compton scattering probabilities
The Compton scattering probabilities IAQ(φ|ϕ geo , E γ ) were computed employing the method GCOMIaq::set using where P eff (ϕ geo , E γ ) is an efficiency factor that is detailed in Appendix G.1, P KN (ϕ geo , E γ ) is the contribution of the Klein-Nishina cross-section to a given bin in ϕ geo that is detailed in Appendix G.2, and R D1 (E 1 |Ê 1 ) and R D2 (E 2 |Ê 2 ) are the energy response functions of the D1 and D2 modules, respectively, that are detailed in Appendix G.3. The transformation from (ϕ geo , E γ ) to (Ê 1 ,Ê 2 ) is done usinĝ gives the D2 energy deposit as a function of the D1 energy deposit (E 1 ) and the Compton scattering angle (φ), and is the Jacobian for the variable transformation. Here K(ϕ geo |ϕ geo , θ) is a Gaussian kernel that provides some smearing of the response in ϕ geo to take into account the event location uncertainties in the D1 and D2 modules. The kernel formally depends on the zenith angle θ of the incoming gamma-ray photons. However, the IAQ is assumed independent of θ, and hence the kernel will be computed for an average zenith angle of θ = 25 • (see Appendix G.4 for details). The Compton scattering probabilities are stored as a twodimensional 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, Gam-maLib 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.

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).

PHINOR
The motivation for the PHINOR algorithm is the observation that the ratio of DRE and DRG multiplied by the solid angle Ω(χ, ψ) of the event cube bins is to first order independent of χ and ψ. This leads to a background model given by where is the convolution of a model of celestial sources I(α, δ, E γ ) with the instrument response function used to subtract any source contribution from the data before normalisation. The evaluation of DRM(χ , ψ ,φ) 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).

BGDLIXA
The BGDLIXA algorithm applies a correction to the PHINOR background model by adjustingφ templates TPL(χ, ψ,φ) over a limited interval ofφ values to a normalisation function n(χ, ψ,φ) that reflects the background event distribution DRE(χ, ψ,φ) − DRM(χ, ψ,φ) smoothed in χ and ψ using a running average (see below). The BGDLIXA algorithm evaluates the expression where the summation is performed over the range A84, page 6 of 31 ∆φ incl = ∆φ (N incl − 1)/2 , ∆φ excl = ∆φ (N excl − 1)/2 , N incl and N excl are either odd integers or zero, and ·/· designates the integer division operator. ∆φ is the bin size inφ of the three-dimensional event cube, which typically is 2 • . While in early days N excl > 0 was used, later analyses set N excl = 0, which simplifies the summation range to Rφ = [φ − ∆φ incl ,φ + ∆φ incl ]. Past COMPTEL analyses generally used N incl = 13 resulting in Rφ = [φ − 12 • ,φ + 12 • ]. As we will see later, N incl is the most critical parameter of the BGDLIXA algorithm, and specifies the number ofφ layers of the three-dimensional event cube over which the templates are adjusted. If N incl is large the fraction in Eq. (12) varies little withφ, leading to a background model that essentially follows the templates. For small values of N 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φ templates are computed using which are an adjustment of the PHINOR background model to match theφ-integrated χ, ψ 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 χ and ψ, with R χ = [χ − ∆χ runav , χ + ∆χ runav ] and R ψ = [ψ − ∆ψ runav , ψ + ∆ψ runav ] where ∆χ runav = ∆χ N runav and ∆ψ runav = ∆ψ N runav . Here, 2N runav + 1 specifies the number of χ and ψ bins over which the running-average summation is performed, with ∆χ and ∆ψ being the bin size in χ and ψ of the three-dimensional event cube, and N runav is an integer number. Usually, ∆χ = ∆ψ = 1 • , and past analyses employed N runav = 3 resulting in a running averaging of ±3 • around each χ, ψ pixel (van Dijk 1996). The adjustment of theφ templates is done using a normalisation function computed using which is an analogue of Eq. (10) that is limited to a small range of χ, ψ around the pixel of interest. Specifically, A χ = [χ − ∆χ avgr , χ + ∆χ avgr ] and A ψ = [ψ − ∆ψ avgr , ψ + ∆ψ avgr ] where ∆χ avgr = ∆χ (N avgr − 1)/2 and ∆ψ avgr = ∆ψ (N avgr − 1)/2 . Here, N avgr is an odd integer that specifies the number of χ and ψ bins over which the adjustment is performed, with ∆χ and ∆ψ being the bin size in χ and ψ of the three-dimensional event cube. Past analyses usually employed N avgr = 3 resulting in an averag- 1996).
We want to point out that the COMPASS analysis system implemented 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φ normalisation step, at the end of the computation so that the model is normalised to the number of background events in eachφ layer.

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 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 runav is no longer a parameter of the model, and the number of pixels in χ and ψ that are used for the local fit is determined by N avgr . We also included ā ϕ 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 runav = 3, N avgr = 3, N incl = 13 and N 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.

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.  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 recomputation 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 maximumlikelihood 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 where T i is the exposure time of viewing period i and T = 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 α, δ 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φ 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 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 (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 χ 2 n distribution with n degrees of freedom, where n is the number of free parameters in the source model component. Therefore, 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 7 was generated by replacing the Crab in the fitted model by a test source with fixed power-law spectral index of Γ = 2. The source visible at 7 The TS map comprises 50 × 50 pixels of size 1 • × 1 • situated around the Crab position. 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 A84, page 9 of 31 A&A 665, A84 (2022) 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 and where ϕ geo is the angular distance between a sky position (α, δ) and the Compton scatter direction (χ, ψ), and [ARM min , ARM max ] defines a selection window for the so-called angular resolution measure that is typically taken to be [−3 • , 3 • ]. We note that DRM designates here the model cube that comprises both the source and background components. By default, comobsres uses to compute the significance of the residuals R(α, δ) in Gaussian σ, 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(α, δ) = 0 the residual significance is while if M(α, δ) = 0 the significance cannot be computed and we set R(α, δ) = 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 • , 3.5 • ].

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 ζ min = 5 • , and we exclude D2 modules for which there were faulty photomultipliers.

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 runav ≤ N avgr , we limit our study here to the simpler BGDLIXE algorithm. Specifically, we investigate which values of N avgr and N 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 excl = 0.

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 avgr and N 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φ dependence in the χ, ψ distribution of the events. The amplitude of the residuals changes strongly under variation of N incl , with larger residuals for larger values of N incl , while variations of N avgr impact the residuals only moderately.
For illustration we show residual maps obtained using the BGDLIXE algorithm for different values of N avgr and N 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 incl , while at the same time smaller value of N 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 incl tend to lead to an underestimation of source fluxes. Therefore, the choice of N 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 avgr , with a slight increase of the residual amplitudes for increasing values of N avgr . 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 incl for N avgr = 5 and as a function of N avgr for N 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 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 avgr so that more events get included in the χ, ψ averaging, as indicated in the lower row of Fig. 10. On the other hand, increasing N avgr leads to a slight increase of the mean and rms distributions; hence, the selected value of N avgr should also not be too large.

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 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 × 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 where p i is the fitted value, v i is the simulated value, and σ 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 σ = 1.  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 incl for N avgr = 5, the lower row shows results as a function of N avgr for N incl = 13. Figure 11 indicates that large values of N 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 incl , as illustrated in Fig. 8, which impact the reconstructed source parameters. Reducing N incl brings the pull distributions more in line with the expectations, yet for N incl ≤ 13 the median pull of the fitted energy flux drops below zero, indicating a systematic bias towards too low fluxes. Specifically, for N incl = 5, where background residuals are very small (cf. Fig. 8), the flux reconstruction is strongly biassed, resulting is a significant underestimation of gammaray fluxes. The optimum value for N 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 avgr , with a minimum bias that occurs for N avgr = 5. Using hence N avgr = 5 and N 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 avgr = 5 and N 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 avgr = 5 and N 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 avgr = 5 and N 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 χ and ψ and 2 • inφ, and that for a different binning the parameters need to be adjusted accordingly. Specifically, we used an equivalent value of N incl = 29 for analyses in our paper for which a binning of 1 • is applied inφ.

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 emission from the Crab pulsar and pulsar wind nebula was studied extensively by COMPTEL in the past (e.g. Much et al. 1995aMuch et al. ,b, 1996van der Meulen et al. 1998;Kuiper et al. 2001).

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 χ and ψ 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φ, 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 avgr = 5 and N incl = 29. The best fit was obtained using the exponentially cut-off power law with k = (1.81 ± 0.06) × 10 −4 ph cm −2 s −1 MeV −1 , Γ = 2.00 ± 0.03, and E c = 39.1 ± 9.7 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 ± 0.05) × 10 −4 ph cm −2 s −1 MeV −1 , Γ = 1.97 ± 0.02, and E c = 29.8 ± 5.0 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 csspec 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 COM-PASS analyses at the time that were not automatically taken 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.

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 8 that is part of the reference data of the X-ray Timing Explorer (XTE) module of the HEASoft software (version 6.29) 9 . We first used compulbin with an angular resolution measure of ±3.5 • to produce pulse profiles for the four standard  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 10 , which could be at the origin of the observed phase shift. 10 1 Energy (MeV)  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σ uncertainty bands of the best fitting spectral models for both components.
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 avgr = 5 and N 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φ 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 COM-PASS 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σ 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 ± 0.5) × 10 −5 ph cm −2 s −1 MeV −1 , Γ = 2.15 ± 0.03, and E c = 53.3 ± 15.5 MeV. For the Crab pulsar the best fitting power-law parameters were k = (5.1 ± 0.3) × 10 −5 ph cm −2 s −1 MeV −1 and Γ = 2.29 ± 0.06. This compares to the spectral indices of 2.227 ± 0.013 and 2.24 ± 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 ± 0.02 that is compatible with their result.

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 × 123 × 25 bins of 1 • × 1 • × 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 ≤ Φ < 0.9 as the inferior conjunction interval (INFC) and phases Φ ≥ 0.9 and Φ < 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   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. 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 γ ) = 1.12 × 10 −4 (E γ /5 MeV) −2.2 ph cm −2 s −1 MeV −1 sr −1 as suggested by Weidenspointner (1999). The background was modelled using the BGDLIXE algorithm with parameters N avgr = 5 and N 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σ 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 used by Collmar & Zhang (2014) which is a minimum distance from the Earth horizon of ζ min = 0 • 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 ζ min = 5 • 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 ζ min = 5 • 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 powerlaw 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. Notes. Bins 1 and 11 serve primarily to constrain continuum gammaray 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.

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 × 100 × 25 bins of 1 • × 1 • × 2 • in size, centred on (l, b) = (286.5 • , 0.5 • ), 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 ± 0.8) × 10 −5 ph cm −2 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 • , 0.5 • ), Knödlseder et al. (1996a) determined a 2σ upper limit of 5.6 • 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φ. 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 A84, page 18 of 31  spectrum we used a Gaussian spectral component with a fixed mean of 1.809 MeV and a standard deviation of σ = 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 γ ) = 1.12 × 10 −4 (E γ /5 MeV) −2.2 ph cm −2 s −1 MeV −1 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 avgr = 5 and N 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 • , 0.5 • ) at which Knödlseder (1994) found the 1.8 MeV line emission maximum and an ARM window of ±3 • . 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 • , 0.5 • ) 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 • , 1.5 • ) while Knödlseder et al. (1996a) found a lower maximum TS value of 14.7 at (l, b) = (286.5 • , 0.5 • ). 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φ 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 • ± 0.8 • , 1.3 • ± 0.7 • ), a flux of (3.1 ± 0.6) × 10 −5 ph cm −2 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 ± 0.8) × 10 −5 ph cm −2 s −1 found by Knödlseder et al. (1996a), yet our best fitting position is offset by about 1.4 • 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 • , which is near the minimum value of 0.001 • that we allowed in the analysis. Using the ctulimit tool we derived a 2σ upper limit of 5.1 • for the diameter of the disk, a bit smaller than the value of 5.6 • 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 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.

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 2 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 2 e (see Appendix H).
This assessment illustrates that the exploitation of archival data instead of the development of new astronomical observatories has the potential to dramatically reduce the carbon footprint of astronomical research, which would help realise the reductions 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. 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 Gam-maLib 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- 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 II or 3 for ToF III , where the latter corrects for energy dependent effects, aligning the forward peak at channel number 120 (see van Dijk 1996 andWeidenspointner 1999 for an explanation of the ToF corrections). If a version 2 EVP is encountered, the GCOMEventList class will automatically convert ToF II values into ToF III using

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 A84, page 22 of 31 J. Knödlseder et al.: COMPTEL data analysis using GammaLib and ctools

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 where ∆t CGRO→SSB corrects for the light travel time from CGRO to the Solar System barycentre, ∆t Shapiro corrects for the gravitational time delay near the Sun, also known as Shapiro delay, ∆t UTC→TT converts from UTC to terrestrial time (TT), and ∆t TT→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 ∆t CGRO→SSB − ∆t Shapiro + ∆t TT→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 is the unit vector in celestial coordinates of a pulsar with right ascension α and declination δ, e is the vector from the Solar System barycentre to the centre of the Earth, s is the vector from the Solar System barycentre to the centre of the Sun, r is the vector from the centre of the Earth to the CGRO spacecraft, and t = 4.92549 × 10 −6 (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 e at a given time t CGRO (specified in the TT time system) is computed using the Taylor expansion e = e(i) +ė(i)∆t + 1 2ë (i)∆t 2 + 1 6 ... e (i)∆t 3 , (D.6) where i is the index of the nearest entry in time in the JPL database, e(i),ė(i),ë(i), and ... e (i) is the Earth position and its first three time derivatives for this entry, and 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 A84, page 23 of 31 A&A 665, A84 (2022) s = s(i). Finally, the conversion from the TT to the TDB time system at a given time t CGRO is computed using ∆t TT→TDB = ∆t TT→TDB (i) + (∆t TT→TDB (i + 1) − ∆t TT→TDB (i))∆t, (D.8) where ∆t TT→TDB (i) is the nearest entry in the JPL database. The last term in the time correction, ∆t UTC→TT , does not depend on planetary ephemerides and is given by where n leap is the number of leap seconds. The computation of ∆t UTC→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 UTC→TT + ∆t TT→TBD , avoiding 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 computes ∆t CGRO→SSB − ∆t Shapiro + ∆t UTC→TT + ∆t TT→TBD , 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: Following the time correction, the pulsar phase Φ is computed using the GPulsarEphemeris::phase method that implements where ν,ν andν is the pulsar frequency and its first two time derivatives, ∆t = t SSB − t 0 is the elapsed time since the reference time of the ephemeris, and Φ 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 maintained during the CGRO mission by the pulsar group at Princeton and nowadays by the Jodrell Bank Centre for Astrophysics. 12 Upon loading of a psrtime file, the GPulsar::load_psrtime method computes from the data the pulsar phase at the reference time using

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σ 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.

Appendix F: Computation of the geometry function
The geometry function DRG is computed in GammaLib by the method GCOMDri::compute_drg using 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 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, 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 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, (θ, φ) being the zenith and azimuth angles of the Compton scatter direction (χ, ψ) with respect to the COMPTEL pointing direction, and o kl (θ, φ) = r 2 1 (α kl − sin α kl cos α kl ) + r 2 2 (β kl − sin β kl cos β kl ) πr 2 1 (F.8) being the projected overlap of a D1 module and a D2 module with cos α kl = d 2 kl + (r 2 1 − r 2 2 ) 2d kl r 1 (F.9) and cos β kl = d 2 kl − (r 2 1 − r 2 2 ) 2d kl r 2 . (F.10) The term f kle (θ, φ) 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 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 contained le (θ, φ) 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.13) with d le = (x e − x l ) 2 + (y e − y l ) 2 (F.14) being the distance between the centres of the exclusion circle and the D2 module, and f contained le (θ, φ) = r 2 1 (α le − sin α le cos α le ) + r 2 2 (β le − sin β le cos β le ) πr 2 1 (F.15) being the overlap of the exclusion region and the D2 module with cos α le = d 2 le + (r 2 1 − r 2 2 ) 2d le r 1 (F.16) and cos β le = d 2 le − (r 2 1 − r 2 2 ) 2d le r 2 . (F.17) The evaluation of Eq. (F.13) is implemented by the method GCOMDri::compute_surface. Finally, f partial kle (θ, φ) 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 eff (ϕ geo , E γ ) in Eq. (6) is factorised according to being the energy of the photon entering the D2 module for a true scatter angle ϕ geo and an incident photon energy of E γ , and E 1 = E γ −Ê 2 being the true energy deposit in D1 for a single Compton scattering by an angle of ϕ geo . The GammaLib methods that implement the computation of the efficiency factors are summarised in Table G.1. Factor Method P eff (ϕ geo , E γ ) GCOMIaq::weight_iaq P A1 (E γ ) GCOMInstChars::trans_D1 P V1 (E γ ) GCOMInstChars::trans_V1 P D1 (E γ ) GCOMInstChars::prob_D1inter P C (E γ ) GCOMIaq::weight_iaq P MH (E γ ) GCOMInstChars::prob_no_multihit P SV (E γ ) GCOMInstChars::prob_no_selfveto P PSD (Ê 1 ) GCOMInstChars::psd_correction P A2 (Ê 2 ) GCOMInstChars::trans_D2 P V23 (Ê 2 ) GCOMInstChars::trans_V23 P D2 (Ê 2 ) GCOMInstChars::prob_D2inter P MS (ϕ geo , E γ ) GCOMInstChars::multi_scatter The P A1 (E γ ) is the transmission probability for photons for the material above D1, which is composed essentially of aluminium, and is computed using where µ Al (E γ ) is the energy-dependent interaction coefficient of aluminium in units of cm −1 that is interpolated using a log-log interpolation of the values given in Table G.2 and l above = 0.147 cm is the thickness of the material above D1. The P V1 (E γ ) is the transmission probability for photons for the first Veto dome and is computed using where µ Veto (E γ ) is the energy-dependent interaction coefficient for the Veto dome in units of cm −1 that is interpolated using a log-log interpolation of the values given in Table G.2 and l V1 = 1.721 cm is the thickness of the first Veto dome. P D1 (E γ ) is the interaction probability in D1 and is computed using where µ D1 (E γ ) is the energy-dependent D1 attenuation coefficient in units of cm −1 that is interpolated using a log-log interpolation of the values given in Table G.2 and l D1 = 8.5 cm is the thickness of the D1 modules. P C (E γ ) is the fraction of Compton scatter interactions among all possible photon interactions within D1 and is given by with the additional constraint of 0 ≤ P Compton (E γ ) ≤ 1. The P MH (E γ ) 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 SV (E γ ) 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 PSD (Ê 1 ) is the D1 energy-dependent PSD correction and is given by P PSD (Ê 1 ) = 1 − 1 1727.9 ×Ê 1 2.53 + 1 , (G.7) whereÊ 1 is in units of MeV. This correction applies to a standard PSD selection of 0-110. P A2 (Ê 2 ) gives the transmission probability of the aluminium below the D1 modules and is computed using P A2 (Ê 2 ) = e −µ Al (Ê 2 ) l T23 , (G.8) where µ Al (Ê 2 ) is the energy-dependent interaction coefficient for aluminium in units of cm −1 that is interpolated using a loglog interpolation of the values given in Table G.2 and l between = 0.574 cm is the thickness of the aluminium plate. P V23 (Ê 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 V23 (Ê 2 ) = e −µ Veto (Ê 2 ) l V23 , (G.9) where µ Veto (Ê 2 ) is the energy-dependent interaction coefficient of the Veto domes in units of cm −1 that is interpolated using a log-log interpolation of the values given in Table G.2 and l V23 = 3.442 cm is the combined thickness of the second and third veto domes. P D2 (E γ ) gives the interaction probability in D2 and is computed using P D2 (E γ ) = 1 − e −µ D2 (Ê 2 ) l D2 , (G.10) where µ D2 (Ê 2 ) is the energy-dependent D2 attenuation coefficient in units of cm −1 that is interpolated using a log-log interpolation of the values given in Table G.2 and l D2 = 7.525 cm is the thickness of the D2 modules. The P MS (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 MS (ϕ geo , E γ ) =   azimuth angle φ of the module and is given by with R(r, φ) = r 2 D1 − r 2 sin 2 φ − r cos φ. (G.13) Figure G.1 illustrates the energy dependence of the efficiency factor P eff (ϕ geo , E γ ) and its components.
The spectral response R D1 (E 1 |Ê 1 ) is illustrated in the left panel of Fig. G.2 for a selected number of input energiesÊ 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.