Issue 
A&A
Volume 665, September 2022



Article Number  A84  
Number of page(s)  31  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202243826  
Published online  15 September 2022 
COMPTEL data analysis using GammaLib and ctools
^{1}
Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, CNES,
9 avenue Colonel Roche,
31028
Toulouse, Cedex 4, France
email: jknodlseder@irap.omp.eu
^{2}
MaxPlanckInstitut für extraterrestrische Physik,
Postfach 1603,
85740
Garching, Germany
^{3}
University of New Hampshire, Space Science Center,
Durham, NH
03824, USA
^{4}
Southwest Research Institute, Dept. of Earth, Oceans, and Space,
Durham, NH
03824, USA
Received:
20
April
2022
Accepted:
3
June
2022
More than 20 yr after the end of NASA’s Compton GammaRay Observatory mission, the data collected by its Imaging Compton Telescope (COMPTEL) still provide the most comprehensive and deepest view of our Universe in MeV gamma rays. While most of the COMPTEL data are archived at NASA’s High Energy Astrophysics Science Archive Research Center (HEASARC), the absence of any publicly available software for their analysis means the data cannot benefit from the scientific advances made in the field of gammaray 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 plugin to the GammaLib library, a communitydeveloped toolbox for the analysis of astronomical gammaray data. We implemented highlevel scripts for building science analysis workflows in ctools, a communitydeveloped gammaray astronomy science analysis software framework. We describe the implementation of our software and provide the underlying algorithms. Using data from the HEASARC archive, we demonstrate that our software reproduces derived data products that were obtained in the past using the proprietary COMPTEL software. We furthermore demonstrate that our software reproduces COMPTEL science results published in the literature. This brings the COMPTEL telescope data back into life, allowing them to benefit from recent advances in gammaray astronomy, and gives the community a means to unveil its still hidden treasures.
Key words: methods: data analysis / gamma rays: general / stars: neutron / binaries: general / ISM: atoms
© J. Knödlseder et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
The Imaging Compton Telescope (COMPTEL) was one of the four telescopes aboard NASA’s Compton GammaRay 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 gammaray sky in the 0.75–30 MeV energy range, and its observations still present the most sensitive survey of theMeV sky ever performed. Since the end of the CGRO mission, our view of the gammaray 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, gammaray 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 gammaray 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 nineyear 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 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 dataanalysis 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 plugin for the GammaLib library, a communitydeveloped toolbox for the analysis of astronomical gammaray data (Knödlseder et al. 2011). The plugin 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 plugin 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 gammaray telescopes, enabling novel analysis approaches that go beyond the capabilities of the COMPASS system.
We furthermore extended the ctools gammaray 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 welldefined ComPTEL data analysis tasks. These building blocks can be combined to create flexible COMPTEL data analysis workflows for the exploration of the MeV sky.
In this paper we present the GammaLib plugin 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}.
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 energydependent 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 gammaray 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 gammaray 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 angle φ was computed from the measured energy deposits E_{1} and E_{2} in D1 and D2, respectively, using (1)
where E = E_{1} + E_{2} is the total energy deposit in the detectors, and m_{e}c^{2} = 0.511 MeV is the rest energy of the electron.
The detector layers were surrounded by an active 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.
3 Software implementation
3.1 GammaLib plugin
The COMPTEL support was implemented in GammaLib as an instrument plugin that provides instrumentspecific implementations of abstract virtual C++ base classes defining the data format and the instrument response functions. The plugin also comprises wrapper functions providing access to all COMPTELspecific C++ classes through the gammalib Python module. All COMPTELspecific classes begin with the letters GCOM.
3.2 COMPTEL data
COMPTEL data are available in the HEASARC archive and are grouped by socalled 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 threeletter 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 socalled 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).
3.3 Event cubes
COMPTEL data analysis is performed on binned event data using threedimensional 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 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 signaltonoise 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 energydependent correction factor is applied to the instrument response. This ToF correction is performed internally by GammaLib (see Appendix C for details). Consequently, flux values returned by GammaLib are always corrected for ToF selection, in contrast to the COMPASS system where the correction had to be applied by the user posterior to the analysis.
An additional, but less important, event selection parameter is the PSD value, which discriminates between neutron and photon interactions in D1, with gamma rays found around channel 80 and neutroninduced 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 socalled 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 (2)
where 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 gammaray 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 phaseresolved analysis of the gammaray 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 phaseresolved analysis of the Crab pulsar and pulsar wind nebula (cf. Sect. 4.2.2).
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).
COMPTELeventcubes arestoredas threedimensionalFITS 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.
Standard event selection parameters that are suitable for most analysis situations and that were used throughout this paper.
Fig. 1 Time of flight spectrum of COMPTEL events registered during ground calibration (left) and in orbit (right). One ToF channel corresponds to 0.25 ns. Downward moving photons interacting first in D1 and then in D2 lead to a peak around channel 120, while upward moving photons interacting first in D2 and then in D1 produce a peak around channel 80. In the calibration data the upward moving photons can be clearly distinguished from the downward moving photons, while in the flight data the downward moving photons are blended with the tail of the dominating upward moving photons. The figure is adapted from Knödlseder (1994). 
Fig. 2 Comparison of an event cube generated with GammaLib for viewing period 2.0 and the 1–3 MeV energy band (red solid) and the equivalent event cube MPEDRE5∅6∅7 produced by the COMPASS software (blue dashed). The left panel shows the distribution of χ values, obtained by summing over all ψ and , the middle panel shows the distribution of ψ values, obtained by summing over all χ and , and the right panel shows the distribution of values, obtained by summing over all χ and ψ. 
3.4 Instrument response function
3.4.1 Factorisation
The COMPTEL instrument response function, given in units of events cm^{2} photons^{−1}, is factorised using (3)
where DRX(a, δ) is the exposure in units of cm^{2} s towards a given celestial direction (α, δ), DRG(x, ψ, ) 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), IAQ is the probability that a gammaray photon with energy E_{γ} being Compton scattered in D1 interacts in D2, with (4)
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 socalled deadtime correction factor, specifying the fraction of time during which the telescope is able to register events. For COMPTEL this fraction is rather constant and was determined by van Dijk (1996) to be L/T = 0.965. This value is hardcoded in GammaLib, and automatically applied to all analysis results.
3.4.2 Exposure
The exposure map DRX(α, δ) is computed in GammaLib by the method GCOMDri::compute_drx using (5)
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 offaxis direction θ_{i}, measured as the angle between the celestial direction (α, δ) and the COMPTEL pointing direction for a given superpacket;. The fraction accounts for the increased interaction length of offaxis 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 twodimensional FITS image with the file type designation DRX. For illustration, Fig. 3 compares projections of the DRX obtained by GammaLib for viewing period 2.0 to the projections for an equivalent DRX that was obtained by COMPASS. The projections of the datasets are indistinguishable, illustrating that the exposure maps generated by GammaLib are equivalent to those produced by the COMPASS software.
3.4.3 Geometry function
The geometry function DRG(χ, ψ, ) 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 threedimensional 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.
Fig. 3 Comparison of an exposure map obtained with GammaLib for viewing period 2.0 (red solid) to the exposure map MPEDRX32382 obtained with COMPASS (blue dashed). The upper panel shows the distribution of Galactic longitude values, obtained by summing over all Galactic latitudes, and the lower panel shows the distribution of Galactic latitude values, obtained by summing over all Galactic longitudes. 
3.4.4 Compton scattering probabilities
The Compton scattering probabilities IAQ were computed employing the method GCOMIaq::set using (6)
where (P_{eff}, E_{γ}) is an efficiency factor that is detailed in Appendix G.1, P_{KN}(, E_{γ}) is the contribution of the KleinNishina crosssection to a given bin in 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 (,E_{γ}) to (Ê_{1}, Ê_{2}) is done using (7)
and Ê_{1} = E_{γ} – Ê_{2}; (8)
gives the D2 energy deposit as a function of the D1 energy deposit (E_{1}) and the Compton scattering angle , and (9)
is the Jacobian for the variable transformation. Here 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 gammaray 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, GammaLib also includes several response functions that were derived using simulations (Stacy et al. 1996), and that can be used as an alternative to the internally computed response functions.
Fig. 4 Comparison of a geometry function obtained with GammaLib for viewing period 2.0 (red solid) to the geometry function MPEDRG35128 obtained with COMPASS (blue dashed). The left panel shows the distribution of χ values, obtained by summing overall ψ and , the middle panel shows the distribution of ψ values, obtained by summing over all χ and , and the right panel shows the distribution of values, obtained by summing overall χ and ψ. 
Fig. 5 Comparison of projections of an IAQ obtained with GammaLib for the energy band 1–3 MeV (red solid) to projections for ROLIAQ75 5 obtained with COMPASS (blue dashed). The left panel shows the distribution of φ_{geo} values, obtained by summing overall , the middle panel shows the distribution of values, obtained by summing overall φ_{geo}, and the right panel shows angular resolution measure . 
3.5 Background modelling
A considerable effort was undertaken by the COMPTEL collaboration to understand and model the instrumental background in the threedimensional 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 socalled SRCLIX method that was developed by Bloemen et al. (1994) and that has undergone several evolutions (van Dijk 1996). The SRCLIX method iteratively applies the BGDLIX algorithm, of which we implemented two variants in GammaLib. For reference, we also implemented the PHINOR algorithm (described below), and all background modelling methods are invoked in GammaLib via the GCOMObservation::compute_drb method. The iterative SRCLIX method is implemented as a ctools script (cf. Sect. 3.6).
3.5.1 PHINOR
The motivation for the PHINOR algorithm is the observation that the ratio of DRE and DRG multiplied by the solid angle Ω(χ, ψ) of the event cube bins is to first order independent of χ and ψ. This leads to a background model given by (10)
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 10MeV. However, for lower energies significant largescale residuals are frequently observed (van Dijk 1996).
3.5.2 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 (12)
where the summation is performed over the range with , , N_{incl} and N_{excl} are either odd integers or zero, and ⌊/⌋ designates the integer division operator. is the bin size in of the threedimensional 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 . Past COMPTEL analyses generally used N_{incl} = 13 resulting in . As we will see later, N_{incl} is the most critical parameter of the BGDLIXA algorithm, and specifies the number of layers of the threedimensional 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 (13)
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 runningaverage summation over a small range in χ and ψ, with R_{x} = [χ – Δχ_{runav},χ + Δχ_{runav}] and R_{ψ} = [ψ – Δψ_{runav},ψ + Δψ_{runav}] where Δχ_{runav} = Δχ N_{ranav} and Δψ_{runav} = Δψ N_{runav}. Here, 2N_{runav} + 1 specifies the number of χ and ψ bins over which the runningaverage summation is performed, with Δχ and Δψ being the bin size in χ and ψ of the threedimensional 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 (14)
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_{ψ} = [ψ – Δψ_{agvar}, ψ + Δψ_{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 threedimensional event cube. Past analyses usually employed N_{avgr} = 3 resulting in an averaging of A_{χ} = [χ – 1°, χ + 1°] and A_{ψ} = [ψ – 1°, ψ + 1°] (van Dijk 1996).
We want to point out that the COMPASS analysis system implemented (15)
instead of Eq. (14), which can lead to unreasonably large normalisations at the edge of the threedimensional 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, (16)
at the end of the computation so that the model is normalised to the number of background events in each layer.
3.5.3 BGDLIXE
Recognising that Eqs. (13) and (14) both do the same thing (that is, they normalise the solidangleweighted geometry function to the measured event distribution), the BGDLIX background model can be simplified to (17)
which is a local fit of the PHINOR background model to the threedimensional 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 a φ 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_{ranav} = 3, N_{avgr} = 3, N_{incl} = 13 and N_{excl} = 0 proposed by van Dijk (1996). While the PHINOR model provides a firstorder description of the event distribution, it shows significant differences in the details. The BGDLIXA and BGDLIXE models follow the event distribution very closely, and both models are in fact indistinguishable from each other, demonstrating that their results are equivalent for the parameter values chosen.
3.6 ctools scripts
To support the science analysis of COMPTEL data, we extended the ctools package with a number of dedicated Python scripts, providing basic building blocks that each perform welldefined 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 COMPTELspecific 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 threedimensional 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 spectralspatial 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 speedup 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; is done using (18) (19) (20) (21)
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 maximumlikelihood 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 maximumlikelihood 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 loglikelihood 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 (22)
(Mattox et al. 1996), where ln L(M_{s} + M_{b}) is the loglikelihood value obtained when fitting the source and the background models together to the data, and ln L(M_{b}) is the loglikelihood 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 distribution with n degrees of freedom, where π is the number of free parameters in the source model component. Therefore, (23)
gives the chance probability (pvalue) that the loglikelihood 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 powerlaw 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 powerlaw spectral index of Γ = 2. The source visible at the centre of the TS maps is the Crab, which is already significantly detected after the first SRCLIX iteration. Subsequent iterations slightly increase the source significance, but overall the TS maps change little. The SRCLIX algorithm converged after six iterations.
The TS maps show a halo of marginal significance around the Crab, which can be explained by the fact that the maps were obtained using a background model that assumed the presence of a point source at the location of the Crab. The event cone of a test source placed a few degrees away from the Crab will pick up some of the excess counts left by the background model, which explains the halo in maps.
An alternative way to generate TS maps is provided by the comlixmap script, which applies the SRCLIX algorithm to each test source position, and which is the algorithm that was implemented by the COMPASS task SRCLIX. Here each pixel in the TS map corresponds to a different background model, and when moving away from the source location, no excess counts will remain in the data. This reduces the halo around the sources, as illustrated in the last panel of Fig. 7 that shows the TS map that was obtained using comlixmap for viewing period 1.0. At the same time, negative residuals are amplified, which is explained by the fact that the events of the Crab for test source positions offset from the source will be included as a smoothed event cone in the background model. We note, however, that this is only a feature of the TS maps, as ultimately, an adequate model of celestial sources that describes the COMPTEL event distribution should be derived from the data. Using that model as DRM in the BGDLIX algorithm will exclude any source events from the smoothing algorithm, and hence provides a reliable background model. For this approach the comsrcdetect script can be used, which extracts significant sources from TS maps and adds them to a model definition XML file. A subsequent run 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 (24)
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 socalled 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 (26)
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 (27)
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°].
Fig. 6 Comparison of PHINOR (blue dotted), BGDLIXA (dark green dashed), and BGDLIXE (light green dasheddotted) background model for viewing period 1.0 and the energy band 1–3 MeV to the DRE event distribution (red solid). The left panel shows the distribution of χ values, obtained by summing overall ψ and , the middle panel shows the distribution of ψ values, obtained by summing over all χ and , and the right panel shows the distribution of values, obtained by summing overall χ and ψ. 
COMPTELspecific ctools scripts.
Fig. 7 Test statistic maps of viewing period 1.0 (during which the Crab pulsar and pulsar wind nebula was observed) for the combined analysis of the four standard COMPTEL energy bands as a function of the number of SRCLIX iterations, indicated in the upperright corner of each panel. The BGDLIXE algorithm with N_{avgr} = 5, N_{incl} = 15, and N_{excl} = 0 was employed. The lowerright panel shows the TS map obtained using the comlixmap script. Colour maps are shown in units of . Negative values correspond to negative fluxes of the test source. 
4 Science validation
Having verified that GammaLib and ctools reproduce data products that are identical to the ones that were generated with the COMPASS system, we now verify that the use of our software for an analysis of the data provided by HEASARC reproduces COMPTEL science results that were published in the literature. If not stated otherwise, for the analyses that follow we apply the event selection parameters specified in Table 1, a value of ζ_{min} = 5°, and we exclude D2 modules for which there were faulty photomultipliers.
4.1 Background model validation
An important step prior to any science analysis is the validation of the background model. An accurate background model predicts the background event distribution within statistical fluctuations, and allows for a reliable and accurate determination of the contribution of celestial gammaray 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 gammaray source characteristics. In agreement with previous studies of the algorithms (cf. van Dijk 1996) we always set N_{excl} = 0.
4.1.1 Residual maps
We started with modelling the background for all COMPTEL viewing periods and the four standard COMPTEL energy bands using comobsback and the BGDLIXE algorithm under variation of the parameters N_{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 gammaray 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 gammaray 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 lowmass Xray 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 tradeoff 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 overntting 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.
Fig. 8 Residual maps for three viewing periods and the 1–3 MeV energy band with background modelled using the BGDLIXE algorithm for N_{avgr} = 5 and different values of N_{incl}. The colour scale is limited to the range −5σ (blue) to +5σ (red). Viewing period 1.0 was selected because it contains the Crab, viewing period 223.0 because it has a rather flat residual, and viewing period 518.5 because it is among the viewing periods with the worst residuals. 
4.1.2 Flux reconstruction
As the next step we studied the impact of the BGDLIXE parameters on the fitted values of celestial source parameters, such as source position, extent, flux and spectral index. We do this by 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 powerlaw 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 comlixf it, 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 (28)
where p_{i} is the fitted value, v_{t} 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 semitransparent 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 tradeoff 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 ψ.
Fig. 10 Violin plots of the mean and rms, in units of σ, of the residual maps with background modelled using the BGDLIXE algorithm for all viewing periods except the ones that include the Crab. The mean and rms are shown as a function of N_{incl} for N_{avgr} = 5 in the upper row and as a function of N_{avgr} for N_{incl} = 13 in the lower row. The groups of four violins correspond to the four standard energy bands; from left to right: 0.75–1 MeV (blue), 1–3 MeV (orange), 3–10 MeV (green), and 10–30 MeV (red). Horizontal black bars show the maximum, the median, and the minimum values, respectively. 
4.2 Gammaray emission from the Crab
We began the science validation of our software with a spectral analysis of the gammaray emission from the Crab pulsar and pulsar wind nebula, which is the brightest source of gamma rays atMeV 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. 1995a,b, 1996; van der Meulen et al. 1998; Kuiper et al. 2001).
Fig. 11 Violin plots of the fitted parameter pull distributions for a simulated source at a 20° offaxis angle for all viewing periods except those that include the Crab. 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. For all simulations the source was simulated as a spatially extended disk with a 3° radius and with a powerlaw spectrum of index Γ = 2.1. Horizontal black bars show the maximum, the median, and the minimum values, respectively. Semitransparent violins with grey borders indicate the expected pull distributions from statistical fluctuations only. 
4.2.1 Total spectrum
We first considered the combined emission from the Crab pulsar and pulsar wind nebula. In their analysis of five years of COMPTEL observations, van der Meulen et al. (1998) analysed the spectrum of the total Crab emission within the energy range 0.78–30 MeV in 30 narrow energy bins. Since this is the only work that quotes total flux values for the Crab (cf. Table 2 of the publication), we used this study as reference.
We analysed the same data that was used by van der Meulen et al. (1998) with GammaLib and ctools, except for viewing period 0 that is not available at HEASARC and viewing period 426.0 for which the EVP file in the HEASARC archive has a corrupted content. We binned the data according to the 30 energy bins defined by van der Meulen et al. (1998) and combined the data for all viewing periods using comobsadd using 80 bins in χ 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 powerlaw, exponentially cutoff powerlaw or curved powerlaw 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 cutoff power law (29)
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 multiwavelength interface to adjust the same spectral models to the data of their Table 2, which also favoured the exponentially cutoff 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 COMPASS analyses at the time that were not automatically taken into account by the software. These correction factors include an energydependent ToF correction factor (cf. Table C.1), an energyindependent deadtime correction factor of 0.965 as well as an energyindependent 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.
Fig. 12 Spectral energy distribution of the total emission from the Crab pulsar and pulsar wind nebula as measured by COMPTEL in 30 bins covering the energy band 0.78–30 MeV. Filled red dots correspond to results obtained with ctools, and open blue dots correspond to results obtained by van der Meulen et al. (1998) using COMPASS. The shaded regions correspond to the 1σ uncertainty bands of the fitted exponentially cutoff powerlaw models. 
4.2.2 Pulsar and nebula components
We now turn to an analysis that separates the emission from the Crab pulsar and the pulsar wind nebula to validate the implementation of the pulsar phase computations. The most comprehensive analysis of the emission from the Crab pulsar and pulsar wind nebula components using COMPTEL data was performed by Kuiper et al. (2001) using data collected over the nine years mission duration of CGRO. While Kuiper et al. (2001) used data from 33 viewing periods with pointing axis within 30° of the Crab pulsar, we analysed 24 viewing periods that we found with the same selection criteria in the HEASARC database, covering the viewing periods specified in Table 1 of Kuiper et al. (2001) between viewing period 1.0 and viewing period 616.1. As in the analysis above, viewing period 426.0 was excluded from the list since no usable EVP file exists for this observation in the HEASARC database. Similarly to Kuiper et al. (2001) we used ephemerides for the Crab pulsar from the Princeton radio pulsar database, provided in the form of an ASCII file named psrtime. dat^{8} that is part of the reference data of the Xray 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.
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 pointsource 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 pointsource 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 COMPASS analysis, as observed above for the total Crab spectrum, and that differences are plausibly explained by the use of a different lists of viewing periods. Kuiper et al. (2001) noticed an enhanced emission in the 10−15 MeV energy interval for the Crab pulsar, and also in our analysis we found an equivalent feature. We note, however, that by shifting the phase interval definition by +0.03 (i.e. using the original phase interval definition of Kuiper et al. 2001) this spectral enhancement is considerably reduced in our analysis, suggesting that the enhancement is probably a statistical fluctuation rather than a physical feature.
We also fitted different spectral models to the data of the Crab pulsar and pulsar wind nebula, including power laws, exponentially cutoff 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 cutoff 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} phem^{−2} s^{−1} MeV^{−1}, Γ = 2.15 ± 0.03, and E_{c} = 53.3 ± 15.5 MeV. For the Crab pulsar the best fitting powerlaw 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 powerlaw 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.
Fig. 13 Pulse profiles of the Crab pulsar derived using ctools for the four standard energy bands (red) compared to pulse profiles derived by Kuiper et al. (2001) using the COMPASS software (blue). The latter profiles were scaled to match our results in amplitude, and a phase shift of −0.03 was applied to match the profiles in phase. 
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. 
4.3 Phaseresolved analysis of LS 5039
We now turn to an analysis of COMPTEL observations of the gammaray binary LS 5039 in order to validate the ability to conduct phaseresolved analyses with GammaLib and ctools. Using an orbitresolved analysis, Collmar & Zhang (2014) found strong evidence that theMeV flux of GRO J182312, 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 J182312 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 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 gammaray 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 northeast 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 3cr 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 powerlaw model using ctbutterfly. The results are shown in Fig. 16. For the left panel we used exactly the same event selection that was 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 13 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 1830210, 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.
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.
Fig. 15 Test statistic maps of LS 5039 for inferior conjunction (left) and superior conjunction (right) phase intervals derived by jointly fitting the four standard COMPTEL energy bands, covering 0.75–30 MeV. The location of LS 5039 is shown by a white plus symbol, and contours show the 1σ, 2σ% and 3σ location uncertainties. The quasar PKS 1830210 was included in the source model and hence is not visible in the maps. 
Fig. 16 Spectral energy distributions and 1er uncertainty bands of fitted powerlaw models for LS 5039 for the INFC (solid lines and filled area) and the SUPC (dashed lines and hatched area). The left panel compares the ctools and GammaLib results (red) to the results obtained by Collmar & Zhang (2014) for an identical event selection (blue). The right panel compares the ctools and GammaLib results (red) to the results obtained using the COMPASS software by excluding D2 modules with faulty PMTs (blue). Data points from Collmar & Zhang (2014) and derived using COMPASS were displaced by 3% in energy for clarity. 
Fig. 17 Orbital light curve of LS 5039 in the 1030 MeV energy band. Data points for the ctools analysis were displaced by ±0.02 in phase for clarity. Analysis results obtained under analysis conditions identical to those in Collmar & Zhang (2014) (ζ_{miτl} = 0° and fpmtflag = 2) are shown as triangles and dashed error bars, and results obtained using ζ_{min} = 5° and fpmtflag = 0 are shown as dots and solid error bars. Vertical lines indicate the definitions of the INFC and SUPC phase intervals. 
Energy bins for 1.8 MeV gammaray line analysis.
4.4 ^{26}AI line emission from Carina
As the next step we validated the capacities of GammaLib and ctools for gammaray line emission analysis together with its abilities to assess the spatial morphology of the emission. As reference, we chose the COMPTEL detection of pointlike 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 pointlike 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 gammaray 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 gammaray emission and the instrumental background. To model the 1.8 MeV line emission 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 gammaray 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 gammaray 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 gammaray 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 pointsource 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 nonnegative 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 gammaray background. The position of the pointsource model as well as the spectral parameters for the Galactic continuum powerlaw 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 theiranalysis. Replacing the point source forthe 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 Fig. 5.5 of Knödlseder (1994) that shows the variation of the maximum loglikelihood value and the fitted 1.8 MeV line flux as a function of spatial disk radius. For that purpose, we determined the maximum likelihood solution for a disk model with position fixed at our maximum likelihood solution (l, b) = (285.4°, 1.3°) for a set of disk radius values starting with 0.001°, and followed by 1° up to 10° with a step size of 1 ^{°}. The results of our analysis are shown in Fig. 20, which also includes the results shown in Fig. 5.5 of Knödlseder (1994) for comparison. The loglikelihood profile between both 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.
Fig. 18 Count spectrum of the Carina region determined for 11 bins covering the energy band 1–3 MeV. The top panel shows the measured number of counts per MeV together with the best fitting background model (green dashed), composed of instrumental background (green dotted) and cosmic gammaray background. The combined source and background model is shown as a black solid line. The bottom panel shows the background modelsubtracted count spectrum together with the best fitting 1.8 MeV line model (red solid) on top of the combined Galactic continuum components (blue dashed), composed of bremsstrahlung (blue dotted) and inverse Compton emission (blue dashdotted). 
Fig. 19 Test statistic maps of 1.8 MeV line emission (left) and 1–3 MeV continuum emission (right) obtained using comlixmap for the Carina region. The map on the left is equivalent to Fig. 2 of Knödlseder et al. (1996a). Dashed contours in the map on the right reflect the intensity of a combination of Galactic bremsstrahlung and inverse Compton emission maps as fitted in an independent analysis where a model of a 1.8 MeV line point source was fitted together with a combination of bremsstrahlung and inverse Compton spatial maps to the data (see text). 
Fig. 20 Variation in maximum loglikelihood with increasing radius of the disk model (left) and fitted 1.8 MeV line flux (right). Data shown as red filled dots and solid lines are those derived in this work, and data shown as blue open dots and dashed lines are taken from Knödlseder (1994). 
5 Conclusion
We have implemented a comprehensive science analysis framework for COMPTEL gammaray 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 mediumenergy gammaray band, from 1–30 MeV, the COMPTEL archive still contains the most sensitive observations ever performed, and a unique dataset for exploring the nonthermal Universe and nuclear transition lines. Since the 1990s when the COMPTEL data were taken, the field of gammaray astronomy has made impressive progress thanks to satellites such as the International GammaRay Astrophysics Laboratory (INTEGRAL) and Fermi and groundbased observatories such as the High Energy Stereoscopic System (H.E.S.S.), the Major Atmospheric Gammaray 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 gammaray 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 postmission discovery of the orbital modulation of MeV gammaray emission of LS 5039 in the COMPTEL archival data by Collmar & Zhang (2014), a source that was not an established gammaray 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 perpublication 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.
Acknowledgements
We would like to thank the anonymous referee for the very careful reading of the manuscript and the many constructive suggestions. This research made use of ctools, a communitydeveloped gammaray astronomy science analysis software (Knödlseder et al. 2016). ctools is based on GammaLib, a communitydeveloped toolbox for the scientific analysis of astronomical gammaray data (Knödlseder et al. 2011). This work has made use of the Python 2D plotting library matplotlib (Hunter 2007). This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research is part of the #LowCarbonScience initiative that aims in reducing the carbon footprint of scientific research, and has benefitted from discussions held in the context of the GDR Labos 1point5 and the Astro4Earth initiative.
Appendix A HEASARC archive
COMPTEL data in the HEASARC archive are grouped by socalled viewing periods with typical durations of two weeks during which the CGRO satellite had a stable pointing. In total, there exist 359 viewing periods, of which 278 were archived at HEASARC, covering viewing periods 1.0 to 719.0, where the first digit indicates the mission year or phase, with viewing periods smaller than 100 corresponding to mission year one. Some of the viewing periods that are available in the HEASARC archive cannot be exploited due to corrupted or missing data files, or spacecraft operations that prevented data taking. The 23 viewing periods in the range 1.0 to 719.0 that are not exploitable are listed in Table A.1. This leaves 255 exploitable viewing periods in the HEASARC archive, which is about 71% of the total number of 359 viewing periods that were executed during the CGRO mission.
Inexploitable viewing periods in the range 1.0 to 719.0 in the COMPTEL HEASARC archive.
The HEASARC archive was created from the original COMPTEL data using dedicated file conversion software that generated FITS files following recognised standards (Pence et al. 2010). While interpreting these FITS files using GammaLib we found that the units in some of the table columns were not correct. We summarise the inconsistencies that we encountered in Table A.2. We furthermore found that the good time interval dataset for viewing period 8.0 with the identifier MPETIM11481 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 projection 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.
Incorrect header units in FITS files in the COMPTEL HEASARC archive.
The HEASARC archive comprises further data products such as telescope housekeeping data, gammaray 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 and Weidenspointner 1999 for an explanation of the ToF corrections). If a version 2 EVP is encountered, the GCOMEventList class will automatically convert ToFii values into ToF_{III} using (B.1)
with the coefficients given by Tables B.1 and B.1 and the D1 and D2 energy deposits E_{1} and E_{2} given in MeV. For the HEASARC archive, ToF values accessed through GammaLib are therefore always ToF_{III} values.
Coefficients, α_{;}, for the conversion from ToF_{II} to ToF_{III}.
Coefficients, b_{i}, for the conversion from ToF_{II} to ToF_{III}.
Appendix C Flux correction due to time of flight selection
To correct for the photon rejection by the ToF selection an energydependent 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 event binning and stored using the header keyword TOFCOR in the DRE FITS files. The correction factor is determined from a linear interpolation of the values given in Table C.1, evaluated at the geometric mean energy of the energy interval [E_{min}, E_{max}] for which the event cube is generated. Only corrections for ToF_{max} = 130 are supported.
Time of flight correction factors for ToF_{max} = 130.
Appendix D Pulsar timing
GammaLib supports generation of event cubes, geometry functions and exposure maps for phaseresolved 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 pulsarspecific code is implemented in GCOMDri::compute_dre. First, the method converts arrival times fcGRO 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 19920625T01: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 (D.1)
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 (D.2)
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 (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 (D.6)
where i is the index of the nearest entry in time in the JPL database, e(i), is the Earth position and its first three time derivatives for this entry, and (D.7)
is the time difference between the event time and the nearest entry in the JPL database in units of days, which by definition is in the range −0.5 and +0.5. Since the Sun moves only a little around the Solar System barycentre, it is sufficient to take for the Sun vector the nearest entry in the JPL database, which is s = s(i). Finally, the conversion from the TT to the TDB time system at a given time t_{CGRO} is computed using (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 (D.9)
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 socalled 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: (D.10)
Following the time correction, the pulsar phase Φ is computed using the GPulsarEphemeris::phase method that implements (D.11)
where v, 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 (D.12)
and MJD_{0} is the modified Julian day of the ephemeris reference (the GEphemerides::geo2ssb method is used for this computation, and hence the Shapiro delay in Eq. (D.13) includes the displacement of the Sun around the Solar System barycentre, cf. Eq. (D.3)).
Appendix E Handling of failed photomultiplier tubes
COMPTEL comprised two detector layers, composed of 7 and 14 circular modules for D1 and D2, respectively. Modules of the first layer were composed of the liquid scintillator NE 213A while modules of the second layer were made of NaI(Tl) scintillator crystals. Each D1 module was viewed by eight PMTs while each D2 module was viewed by seven PMTs. The relative amplitudes of the PMT signals for a given module allowed localisation of the interactions within the module with an average 1 σ 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.
D2 module exclusion circles applied following the failure of PMTs at the specified dates.
Appendix F Computation of the geometry function
The geometry function DRG is computed in GammaLib by the method GCOMDri::compute_drg using (F.1)
where EHA(χ, ψ) is the distance between scatter direction (χ, ψ) and the Earth horizon, 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 (F.3)
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, (F.6)
with r_{1} = 13.8 cm being the radius of a D1 module and r_{2} = 14.085 cm being the radius of a D2 module. We note the margin of 0.1 cm in Eq. (F.6), which assures numerical stability with respect to rounding errors. The d_{kl} is the projected distance between the centres of the D1 and D2 modules, given by (F.7)
where x_{k} and y_{k} are the geometric positions of the D1 modules and x_{l} and y_{l} the positions of the D2 modules with respect to the optical axis in cm, h = 158 cm being the vertical separation between D1 and D2 modules, (θ, ϕ) being the zenith and azimuth angles of the Compton scatter direction (χ, ψ) with respect to the COMPTEL pointing direction, and (F.8)
being the projected overlap of a D1 module and a D2 module with (F.9)
The term f_{k1e}(θ, ϕ) accounts for the exclusion of circular regions around failed PMTs and differs from zero only for fpmtflag = 2. It quantifies the fractional overlap between the projected D1 module k and the part of the exclusion region for failed PMTs that is contained within D2 module l. The exclusion region is circular, and specified by a geometric centre position x_{e} and y_{e} and a radius r_{e} as given in Table E.1. Specifically, for fpmtflag = 2 (F.11)
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 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)
being the distance between the centres of the exclusion circle and the D2 module, and (F.15)
being the overlap of the exclusion region and the D2 module with (F.16)
The evaluation of Eq. (F.13) is implemented by the method GCOMDri::compute_surface.
Finally, 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 (G.1)
being the energy of the photon entering the D2 module for a true scatter angle φ′_{geo} and an incident photon energy of E_{γ}, and Ê_{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.
GammaLib methods that implement the computation of the efficiency factors.
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 (G.3)
where μ_{Al}(E_{γ}) is the energydependent interaction coefficient of aluminium in units of cm^{−1} that is interpolated using a loglog 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 (G.4)
where μ_{Veto}(E_{γ}) is the energydependent interaction coefficient for the Veto dome in units of cm^{−1} that is interpolated using a loglog 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 (G.5)
where μ_{D1}(E_{γ}) is the energydependent D1 attenuation coefficient in units of cm^{−1} that is interpolated using a loglog 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 (G.6)
with the additional constraint of 0 ≤ P_{Compton}(E_{γ}) ≤ 1.
The P_{MH}(E_{γ}) is the probability that there is no multihit. This probability is computed using a loglog interpolation of the values of Table G.2. P_{SV}(E_{γ}) is the probability that the photon was not selfvetoed. This probability is computed using a linear interpolation of the values of Table G.2.
The P_{PSD}(Ê_{1}) is the D1 energydependent PSD correction and is given by (G.7)
where Ê_{1} is in units of MeV. This correction applies to a standard PSD selection of 0110.
P_{A2}(Ê_{2}) gives the transmission probability of the aluminium below the D1 modules and is computed using (G.8)
where μ_{A1}(Ê_{2}) is the energydependent 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 (G.9)
where μ_{Veto}(Ê_{2}) is the energydependent interaction coefficient of the Veto domes in units of cm^{1} that is interpolated using a loglog 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 (G.10)
where μ_{D2}(Ê_{2}) is the energydependent D2 attenuation coefficient in units of cm^{−1} that is interpolated using a loglog 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 (G.11)
which integrates all possible second interaction locations within a given D1 module assuming vertically incident photons, where r_{D1} = 13.8 cm and l_{D1} = 8.5 cm is the radius and thickness of a D1 module. 1(r, z,ϕ) gives the remaining interaction length for a photon that has interacted at radius r, vertical depth z and azimuth angle ϕ of the module and is given by
Efficiency parameters as a function of energy.
Fig. G.1 Energy dependence of the efficiency factor P_{eff}(φ′_{geo}, E_{γ}) for the three geometrical scatter angles φ′_{geo} = 5°, 25°, and 50° (left) and the energy dependence of its components for φ′_{geo} = 25° (right). 
Figure G.1 illustrates the energy dependence of the efficiency factor P_{eff}(φ′_{geo}, E_{γ}) and its components.
Appendix G.2 KleinNishina contribution
is the contribution of the KleinNishina crosssection σ_{KN}(φ″_{geo}, E_{γ}) to the φ′_{geo} bin spanned by [φ′_{geo,min}, φ′_{geo,max}] where the integrals are computed using (G.15)
Appendix G.3 Module spectral response functions
Appendix G.3.1 D1 spectral response
The spectral response R_{D1}(E_{1}Ê_{1}) provides the probability for measuring an energy E1 when an incident photon with energy E_{γ} was Compton scattered by an angle φ_{geo}, producing a true energy deposit of Ê_{1} in the detector module. The D1 response is implemented by the GCOMD1Response class and the response evaluation is done using the GCOMD1Response::operator() operator. The spectral D1 response is computed using a Gaussian (G.17)
where μ_{p}(Ê_{1}) is the position of the Gaussian peak for true energy deposit of Ê_{1}, with μ_{p}(Ê_{1}) ≈ Ê_{1}, and σ(Ê_{1}) is the standard deviation of the Gaussian (Diehl et al. 1992). T_{1}(E_{1} Ê_{1}) is a threshold function, defined by (G.18)
where . μ_{p}(Ê_{1}) and σ(Ê_{1}) as well as and ∆E_{1}(Ê_{1}) are tabulated as a function of Ê_{1} in a calibration file of type SDA that is included in GammaLib. Values for a given energy Ê_{1} are obtained by linear interpolation of the tabulated values.
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.
Fig. G.2 Detector response functions R_{D1}(E_{1}Ê_{1}) (left) and R_{D2}(E_{2}Ê_{2}) (right) for a selected number of energies that have been chosen to allow for a comparison with Figs. 3c and 4b in Diehl et al. (1992). 
Appendix G.3.2 D2 spectral response
The spectral response R_{D2}(E_{2}Ê_{2}) provides the probability for measuring an energy E_{2} when an photon with energy Ê_{2} hit a D2 detector module. The D2 response is implemented by the GCOMD2Response class and the response evaluation is done using the GCOMD2Response::operator() operator. The spectral D2 response is computed using (G.19)
which is a combination of three Gaussian and two continuum components multiplied by a threshold function T_{2}(E_{2}Ê_{2}) (Diehl et al. 1992). Here, Ê_{2} is the energy of the gammaray photon incident on a D2 module, and E_{2} is the measured energy deposit.
describes the socalled photopeak which represents photons that were fully absorbed in the D2 module We note that the position μ_{PP}(Ê_{2}) of the photopeak is not necessarily identical to the incident energy, although below ~ 10 MeV μ_{pp}(Ê_{2}) ≈ Ê_{2} is well satisfied.
The two subsequent Gaussians describe the single and double escape peak, respectively, which are generated by photons above 1.022 MeV that undergo pair creation resulting in the production of two 511 keV gamma rays. The escape of one or both of these 511 keV photons from the D2 module leads to characteristic peaks at energies 0.511 MeV and 1.022 MeV below the photopeak energy. The escape peaks are modelled using (G.21)
with μ_{Se}(Ê_{2}) = μ_{pp}(Ê_{2}) − m_{e}c^{2} and μ_{de}(Ê_{2}) = μ_{pp}(Ê_{2}) − 2m_{e}c^{2}, where m_{e}c^{2} ≈ 511 keV is the rest energy of the electron.
The first continuum component, (G.23)
describes the socalled Compton tail which represents photons that underwent a Compton scattering before escaping the D2 module. The component is the product of the KleinNishina crosssection for Compton scattering, (G.24)
and the probability that the Compton scattered photon does not undergo a second Compton scattering before escaping, empirically modelled using (G.25)
is the total linear attenuation coefficient in NaI, which is the scintillator material of the D2 modules and (G.27)
is an empirical path length in the D2 module, where Ê_{2} is in units of MeV. The product is convolved with a Gaussian kernel, (G.28)
to take the energy resolution of the D2 module into account, where σ(E_{2}) is the standard deviation of the detector response at the reconstructed energy E_{2}.
The second continuum component, (G.29)
is a flat background component that is convolved by a Gaussian kernel and that takes into account any higherorder scatterings in the D2 modules.
The threshold function is defined by (G.30)
where and
All five response components have amplitude parameters (A_{pp}, A_{se}, A_{de}, A_{com} and A_{bkg}) that are tabulated as a function of incident energy E_{2} in a calibration file of type SDB. In addition, this file contains also tabulated values for the photopeak energy μ_{PP}(Ê_{2}), the standard deviation σ(Ê_{2}) and the minimum and maximum threshold energies and and width ∆E_{2}(Ê_{2}). Values for a given energy Ê_{2} are obtained by linear interpolation of the tabulated values.
The spectral response R_{D2}(E_{2}Ê_{2}) is illustrated in the right panel of Fig. G.2 for a selected number of input energies Ê_{2}. The figure can be compared with Fig. 4b in Diehl et al. (1992), which presents the response functions as implemented in COMPASS at the time of the CGRO launch for comparable energies.
Appendix G.4 Response convolution
The response is convolved with the Gaussian kernel (G.31)
to take into account the location uncertainties of the events in the D1 and D2 modules. The width σ(φ′_{geo}, θ) of the Gaussian kernel is a function of the scatter angle φ′_{geo} and the zenith angle θ of the incoming photon, measured with respect to the COMPTEL pointing axis. The width (G.32)
is the squared mean of the uncertainty (G.33)
in the horizontal direction and the uncertainty (G.34)
in the vertical direction, where σ_{x} = 2.3 cm and σ_{y} = 1.96 cm is the average horizontal location spread in a D1 module with respect to the D2 module in x and y direction, respectively, σ_{z1} = 2.45 cm and σ_{z2} = 2.17 cm is the vertical location spread in a D1 and D2 module, respectively, assuming a homogenous event distribution in the zdirection, and h = 158 cm is the distance between the D1 and D2 modules. The functions (G.35)
are the geometrical relations that transform from the uncertainties within the modules to uncertainties in φ′_{geo} for a given zenith angle θ.
Fig. G.3 Width σ(φ′_{geo}, θ) of the Gaussian smoothing kernel as a function of φ′_{geo} for zenith angles θ of 0° (solid blue), 25° (dashed orange), 35^{°} (dasheddotted green), and 50^{°} (dotted purple). 
Since COMPTEL response functions have no explicit zenith angle dependence, the computations are performed for an average zenith angle of θ = 25°. The uncertainty introduced by this approximation is illustrated in Fig. G.3, and amounts to less than 10% for zenith angle θ < 35°.
Appendix H Carbon footprint of this research
Given the repeated alarms of the scientific community on global warming (IPCC 2021, 2022), biodiversity loss (IPBES 2019) and air, water and soil pollution (Landrigan et al. 2017), the environmental impact of any human activity, including obviously research, can no longer be ignored. Here we assess the GHG emissions related to our research work, including contributions from office use, research activities and research environment. We based our assessment on the carbon footprint estimate of the Institut de Recherche en astrophysique et planétologie (IRAP; Martin et al. 2022) where most of the work for this research was done. To attribute a share of the total footprint to our work, we estimated that roughly one full time equivalent (FTE) was necessary to accomplish this research, corresponding to 220 working days. Consequently, we use the annual emissions per source from Martin et al. (2022) divided by the IRAP staff of 263 persons to assess the contributions from the various sources of GHG emissions.
Carbon footprint estimate by emission source of the research work behind this paper.
The results of our assessment are summarised in Table H.1. We estimate the total carbon footprint of our research work to be 1.9 ± 0.3 tCO_{2}e, composed of 1.2 ± 0.2 tCO_{2}e (65%) for running the office building in which the work was done, 0.4 ± 0.2 tCO_{2}e (21%) for the computing infrastructure and data flow, and 0.3 ± 0.1 tCO_{2}e (14%) for lunch meals and hometooffice commuting.
We based our assessment for computing on the carbon footprint estimation implemented in ctools, yet when we started this research the relevant functionality was not yet available. We note, however, that the carbon footprint estimate for the final analyses for this publication amounted to 60 kgCO_{2}e, of which 20 kgCO_{2}e were attributed to power consumption and 40 kgCO_{2}e to the computing infrastructure. Making the conservative estimate that ten analysis iterations were necessary for this research work we obtain a total carbon footprint due to computing of 600 kgCO_{2}e, of which 200 kgCO_{2}e were due to power consumption and 400 kgCO_{2}e due to the computing infrastructure. To avoid double counting, we disregarded the contribution of power consumption since it is already accounted for in the electricity to power the office building.
For emission from lunch meals and hometooffice commuting we took the specific habits of the coauthors weighted by their contribution to the research into account. We note that our research work did not require any professional travelling. The footprint for a couple of video conferences that we organised for this research is included under the ‘Data flow’ source.
References
 Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 [Google Scholar]
 Berthoud, F., Bzeznik, B., Gibelin, N., et al. 2020, Estimation de l’empreinte carbone d’une heure.coeur de calcul, Research report, UGA  Université Grenoble Alpes; CNRS; INP Grenoble; INRIA [Google Scholar]
 Bloemen, H., Diehl, R., Hermsen, W., Knödlseder, J., & Schönfelder, V. 1999, Astrophys. Lett. Commun., 38, 391 [NASA ADS] [Google Scholar]
 Bloemen, H., Hermsen, W., Swanenburg, B. N., et al. 1994, ApJS, 92, 419 [NASA ADS] [CrossRef] [Google Scholar]
 Casares, J., Ribó, M., Ribas, I., et al. 2005, MNRAS, 364, 899 [Google Scholar]
 Cash, W. 1979, ApJ, 228, 93 [Google Scholar]
 Coleman, A., & McConnell, M. 2020, in American Astronomical Society Meeting Abstracts, 235, 206.04 [NASA ADS] [Google Scholar]
 Collmar, W., & Zhang, S. 2014, A&A, 565, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collmar, W., Wessolowski, U., Schönfelder, V., et al. 1997, in American Institute of Physics Conference Series, 410, Proceedings of the Fourth Compton Symposium, eds. C. D. Dermer, M. S. Strickman, & J. D. Kurfess, 1587 [NASA ADS] [CrossRef] [Google Scholar]
 de Vries, C. P., & COMPTEL collaboration. 1994, in Astronomical Society of the Pacific Conference Series, 61, Astronomical Data Analysis Software and Systems III, eds. D. R. Crabtree, R. J. Hanisch, & J. Barnes, 399 [Google Scholar]
 Diehl, R., Aarts, H., Bennett, K., et al. 1992, in Data Analysis in Astronomy, 201 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 IPBES 2019, Global assessment report on biodiversity and ecosystem services of the Intergovernmental SciencePolicy Platform on Biodiversity and Ecosystem Services, eds. E. S. Brondizio, J. Settele, S. Díaz, & H. T. Ngo (IPBES secretariat, Bonn, Germany) [Google Scholar]
 IPCC 2021, Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, eds. V. MassonDelmotte, P. Zhai, A. Pirani, et al. (Cambridge University Press) [Google Scholar]
 IPCC 2022, Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, eds. H.O. Pörtner, M.D.C. Roberts, E. Tignor, et al. (Cambridge University Press) [Google Scholar]
 Knödlseder, J. 1994, PhD thesis, Technical University Munich, Germany [Google Scholar]
 Knödlseder, J., Bennett, K., Bloemen, H., et al. 1996a, A&As, 120, 327 [Google Scholar]
 Knödlseder, J., von Ballmoos, P., Diehl, R., et al. 1996b, SPIE Conf. Ser., 2806, 386 [Google Scholar]
 Knödlseder, J., Mayer, M., Deil, C., et al. 2016, Astrophysics Source Code Library [ascl:1601.005] [Google Scholar]
 Knödlseder, J., Mayer, M., Deil, C., et al. 2011, Astrophysics Source Code Library [ascl:1110.007] [Google Scholar]
 Knödlseder, J., BrauNogué, S., Coriat, M., et al. 2022, Nat. Astron., 6, 503 [CrossRef] [Google Scholar]
 Kuiper, L., Hermsen, W., Cusumano, G., et al. 2001, A&A, 378, 918 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Landrigan, P., Fuller, R., Acosta, N., et al. 2017, Lancet, 391 [Google Scholar]
 Martin, P., BrauNogué, S., Coriat, M., et al. 2022, ArXiv eprints, [arXiv:2204.12362] [Google Scholar]
 Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 [Google Scholar]
 Much, R., Bennett, K., Buccheri, R., et al. 1995a, A&A, 299, 435 [NASA ADS] [Google Scholar]
 Much, R., Bennett, K., Buccheri, R., et al. 1995b, Adv. Space Res., 15, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Much, R., Harmon, B. A., Nolan, P., et al. 1996, A&As, 120, 703 [NASA ADS] [Google Scholar]
 Pence, W. D., Chiappetti, L., Page, C. G., Shaw, R. A., & Stobie, E. 2010, A&A, 524, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schönfelder, V., Aarts, H., Bennett, K., et al. 1993, ApJS, 86, 657 [NASA ADS] [CrossRef] [Google Scholar]
 Stacy, J. G., Kippen, R. M., Kappadath, S. C., et al. 1996, A&As, 120, 691 [NASA ADS] [Google Scholar]
 Strong, A., & Collmar, W. 2019, Mem. Soc. Astron. Italiana, 90, 297 [Google Scholar]
 Strong, A. W., Bennett, K., Bloemen, H., et al. 1996, A&As, 120, 381 [NASA ADS] [Google Scholar]
 van der Meulen, R. D., Bloemen, H., Bennett, K.,et al. 1998, A&A, 330, 321 [NASA ADS] [Google Scholar]
 van Dijk, R. C. A. 1996, PhD thesis, University of Amsterdam, The Netherlands [Google Scholar]
 Weidenspointner, G. 1999, PhD thesis, Technical University Munich, Germany [Google Scholar]
 Weidenspointner, G., Varendorff, M., Oberlack, U., et al. 2001, A&A, 368, 347 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Some viewing periods in the HEASARC archive are unreadable, and data from the end of the mission were not archived. Details about issues encountered in the HEASARC COMPTEL data archive are provided in Appendix A.
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.
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 percentlevel differences that are plausibly attributed to differences in our and the COMPASS processing configuration that we were not able to track down.
Details are provided at https://heasarc.gsfc.nasa.gov/lheasoft/ftools/fhelp/fasebin.html
The HEASoft software can be downloaded from https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/
We used the bremsstrahlung emission map with the COMPASS identifier MPEMIS13006 and the inverse Compton emission map with the COMPASS identifier ROLMIS163 that were also used by Collmar & Zhang (2014).
All Tables
Standard event selection parameters that are suitable for most analysis situations and that were used throughout this paper.
Inexploitable viewing periods in the range 1.0 to 719.0 in the COMPTEL HEASARC archive.
D2 module exclusion circles applied following the failure of PMTs at the specified dates.
Carbon footprint estimate by emission source of the research work behind this paper.
All Figures
Fig. 1 Time of flight spectrum of COMPTEL events registered during ground calibration (left) and in orbit (right). One ToF channel corresponds to 0.25 ns. Downward moving photons interacting first in D1 and then in D2 lead to a peak around channel 120, while upward moving photons interacting first in D2 and then in D1 produce a peak around channel 80. In the calibration data the upward moving photons can be clearly distinguished from the downward moving photons, while in the flight data the downward moving photons are blended with the tail of the dominating upward moving photons. The figure is adapted from Knödlseder (1994). 

In the text 
Fig. 2 Comparison of an event cube generated with GammaLib for viewing period 2.0 and the 1–3 MeV energy band (red solid) and the equivalent event cube MPEDRE5∅6∅7 produced by the COMPASS software (blue dashed). The left panel shows the distribution of χ values, obtained by summing over all ψ and , the middle panel shows the distribution of ψ values, obtained by summing over all χ and , and the right panel shows the distribution of values, obtained by summing over all χ and ψ. 

In the text 
Fig. 3 Comparison of an exposure map obtained with GammaLib for viewing period 2.0 (red solid) to the exposure map MPEDRX32382 obtained with COMPASS (blue dashed). The upper panel shows the distribution of Galactic longitude values, obtained by summing over all Galactic latitudes, and the lower panel shows the distribution of Galactic latitude values, obtained by summing over all Galactic longitudes. 

In the text 
Fig. 4 Comparison of a geometry function obtained with GammaLib for viewing period 2.0 (red solid) to the geometry function MPEDRG35128 obtained with COMPASS (blue dashed). The left panel shows the distribution of χ values, obtained by summing overall ψ and , the middle panel shows the distribution of ψ values, obtained by summing over all χ and , and the right panel shows the distribution of values, obtained by summing overall χ and ψ. 

In the text 
Fig. 5 Comparison of projections of an IAQ obtained with GammaLib for the energy band 1–3 MeV (red solid) to projections for ROLIAQ75 5 obtained with COMPASS (blue dashed). The left panel shows the distribution of φ_{geo} values, obtained by summing overall , the middle panel shows the distribution of values, obtained by summing overall φ_{geo}, and the right panel shows angular resolution measure . 

In the text 
Fig. 6 Comparison of PHINOR (blue dotted), BGDLIXA (dark green dashed), and BGDLIXE (light green dasheddotted) background model for viewing period 1.0 and the energy band 1–3 MeV to the DRE event distribution (red solid). The left panel shows the distribution of χ values, obtained by summing overall ψ and , the middle panel shows the distribution of ψ values, obtained by summing over all χ and , and the right panel shows the distribution of values, obtained by summing overall χ and ψ. 

In the text 
Fig. 7 Test statistic maps of viewing period 1.0 (during which the Crab pulsar and pulsar wind nebula was observed) for the combined analysis of the four standard COMPTEL energy bands as a function of the number of SRCLIX iterations, indicated in the upperright corner of each panel. The BGDLIXE algorithm with N_{avgr} = 5, N_{incl} = 15, and N_{excl} = 0 was employed. The lowerright panel shows the TS map obtained using the comlixmap script. Colour maps are shown in units of . Negative values correspond to negative fluxes of the test source. 

In the text 
Fig. 8 Residual maps for three viewing periods and the 1–3 MeV energy band with background modelled using the BGDLIXE algorithm for N_{avgr} = 5 and different values of N_{incl}. The colour scale is limited to the range −5σ (blue) to +5σ (red). Viewing period 1.0 was selected because it contains the Crab, viewing period 223.0 because it has a rather flat residual, and viewing period 518.5 because it is among the viewing periods with the worst residuals. 

In the text 
Fig. 9 Same as Fig. 8 but for N_{incl} = 13 and different values of N_{avgr}. 

In the text 
Fig. 10 Violin plots of the mean and rms, in units of σ, of the residual maps with background modelled using the BGDLIXE algorithm for all viewing periods except the ones that include the Crab. The mean and rms are shown as a function of N_{incl} for N_{avgr} = 5 in the upper row and as a function of N_{avgr} for N_{incl} = 13 in the lower row. The groups of four violins correspond to the four standard energy bands; from left to right: 0.75–1 MeV (blue), 1–3 MeV (orange), 3–10 MeV (green), and 10–30 MeV (red). Horizontal black bars show the maximum, the median, and the minimum values, respectively. 

In the text 
Fig. 11 Violin plots of the fitted parameter pull distributions for a simulated source at a 20° offaxis angle for all viewing periods except those that include the Crab. 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. For all simulations the source was simulated as a spatially extended disk with a 3° radius and with a powerlaw spectrum of index Γ = 2.1. Horizontal black bars show the maximum, the median, and the minimum values, respectively. Semitransparent violins with grey borders indicate the expected pull distributions from statistical fluctuations only. 

In the text 
Fig. 12 Spectral energy distribution of the total emission from the Crab pulsar and pulsar wind nebula as measured by COMPTEL in 30 bins covering the energy band 0.78–30 MeV. Filled red dots correspond to results obtained with ctools, and open blue dots correspond to results obtained by van der Meulen et al. (1998) using COMPASS. The shaded regions correspond to the 1σ uncertainty bands of the fitted exponentially cutoff powerlaw models. 

In the text 
Fig. 13 Pulse profiles of the Crab pulsar derived using ctools for the four standard energy bands (red) compared to pulse profiles derived by Kuiper et al. (2001) using the COMPASS software (blue). The latter profiles were scaled to match our results in amplitude, and a phase shift of −0.03 was applied to match the profiles in phase. 

In the text 
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. 

In the text 
Fig. 15 Test statistic maps of LS 5039 for inferior conjunction (left) and superior conjunction (right) phase intervals derived by jointly fitting the four standard COMPTEL energy bands, covering 0.75–30 MeV. The location of LS 5039 is shown by a white plus symbol, and contours show the 1σ, 2σ% and 3σ location uncertainties. The quasar PKS 1830210 was included in the source model and hence is not visible in the maps. 

In the text 
Fig. 16 Spectral energy distributions and 1er uncertainty bands of fitted powerlaw models for LS 5039 for the INFC (solid lines and filled area) and the SUPC (dashed lines and hatched area). The left panel compares the ctools and GammaLib results (red) to the results obtained by Collmar & Zhang (2014) for an identical event selection (blue). The right panel compares the ctools and GammaLib results (red) to the results obtained using the COMPASS software by excluding D2 modules with faulty PMTs (blue). Data points from Collmar & Zhang (2014) and derived using COMPASS were displaced by 3% in energy for clarity. 

In the text 
Fig. 17 Orbital light curve of LS 5039 in the 1030 MeV energy band. Data points for the ctools analysis were displaced by ±0.02 in phase for clarity. Analysis results obtained under analysis conditions identical to those in Collmar & Zhang (2014) (ζ_{miτl} = 0° and fpmtflag = 2) are shown as triangles and dashed error bars, and results obtained using ζ_{min} = 5° and fpmtflag = 0 are shown as dots and solid error bars. Vertical lines indicate the definitions of the INFC and SUPC phase intervals. 

In the text 
Fig. 18 Count spectrum of the Carina region determined for 11 bins covering the energy band 1–3 MeV. The top panel shows the measured number of counts per MeV together with the best fitting background model (green dashed), composed of instrumental background (green dotted) and cosmic gammaray background. The combined source and background model is shown as a black solid line. The bottom panel shows the background modelsubtracted count spectrum together with the best fitting 1.8 MeV line model (red solid) on top of the combined Galactic continuum components (blue dashed), composed of bremsstrahlung (blue dotted) and inverse Compton emission (blue dashdotted). 

In the text 
Fig. 19 Test statistic maps of 1.8 MeV line emission (left) and 1–3 MeV continuum emission (right) obtained using comlixmap for the Carina region. The map on the left is equivalent to Fig. 2 of Knödlseder et al. (1996a). Dashed contours in the map on the right reflect the intensity of a combination of Galactic bremsstrahlung and inverse Compton emission maps as fitted in an independent analysis where a model of a 1.8 MeV line point source was fitted together with a combination of bremsstrahlung and inverse Compton spatial maps to the data (see text). 

In the text 
Fig. 20 Variation in maximum loglikelihood with increasing radius of the disk model (left) and fitted 1.8 MeV line flux (right). Data shown as red filled dots and solid lines are those derived in this work, and data shown as blue open dots and dashed lines are taken from Knödlseder (1994). 

In the text 
Fig. G.1 Energy dependence of the efficiency factor P_{eff}(φ′_{geo}, E_{γ}) for the three geometrical scatter angles φ′_{geo} = 5°, 25°, and 50° (left) and the energy dependence of its components for φ′_{geo} = 25° (right). 

In the text 
Fig. G.2 Detector response functions R_{D1}(E_{1}Ê_{1}) (left) and R_{D2}(E_{2}Ê_{2}) (right) for a selected number of energies that have been chosen to allow for a comparison with Figs. 3c and 4b in Diehl et al. (1992). 

In the text 
Fig. G.3 Width σ(φ′_{geo}, θ) of the Gaussian smoothing kernel as a function of φ′_{geo} for zenith angles θ of 0° (solid blue), 25° (dashed orange), 35^{°} (dasheddotted green), and 50^{°} (dotted purple). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.