Highlight
Open Access
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/0004-6361/202243826
Published online 15 September 2022

© J. Knödlseder et al. 2022

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

While most of the COMPTEL data accumulated during the nine-year mission are stored at NASA s High Energy Astrophysics Science Archive Research Center (HEASARC)1, no publicly available software had existed to analyse these data. During CGRO operations the COMPTEL data were analysed using the COMPTEL Processing and Analysis Software System (COMPASS) (de Vries & COMPTEL Collaboration 1994), but this system was only available at the COMPTEL collaboration institutes and was decommissioned after the end of the mission. Nevertheless, a few science results were still obtained from COMPTEL data after the end of the mission thanks to efforts at some COMPTEL collaboration institutes to keep Linux ports of the COMPASS data-analysis system alive (Strong & Collmar 2019; Coleman & McConnell 2020).

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

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

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

2 The COMPTEL telescope

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

COMPTEL was composed of two modular detector layers D1 and D2, separated by 158 cm, where an incident gamma-ray photon was first Compton scattered in one of the 7 modules of D1 and then eventually interacted in one of the 14 modules of D2. The Compton scatter direction (χ, ψ) was obtained from the interaction locations in D1 and D2, the Compton scattering angle φ was computed from the measured energy deposits E1 and E2 in D1 and D2, respectively, using (1)

where E = E1 + E2 is the total energy deposit in the detectors, and mec2 = 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 background3 due to charged particles. A further strong discriminator of instrumental background was the time of flight (ToF) measurement, which is the time difference between the interactions in D1 and D2 (cf. Sect. 3.3).

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

3 Software implementation

3.1 GammaLib plug-in

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

3.2 COMPTEL data

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

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

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

3.3 Event cubes

COMPTEL data analysis is performed on binned event data using three-dimensional event cubes spanned by the Compton scatter direction (χ,ψ) 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 signal-to-noise ratio, and Fig. 1 illustrates that the minimum ToF value has a strong impact on the background discrimination, with larger values removing more background at the expense of rejecting also an increasing number of celestial photons.

To correct for the photon rejection by the ToF selection an energy-dependent correction factor is applied to the instrument response. This ToF correction is performed internally by GammaLib (see Appendix C for details). Consequently, flux values returned by GammaLib are always corrected for ToF selection, in contrast to the COMPASS system where the correction had to be applied by the user posterior to the analysis.

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

Furthermore, the Earth atmosphere being a strong gamma-ray source, an important background reduction is obtained by excluding events that may originate from the Earth atmosphere. This is achieved by requiring that event circles have a minimum angular distance ζmin from the Earth horizon. In the process of event binning, this requirement is translated into a constraint on the so-called Earth horizon angle (EHA), which is defined as the angular distance between the scatter direction of an event and the Earth horizon. Specifically, only the events are retained that satisfy EHA ≥ EHAmin 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 file5.

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

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 three-dimensionalFITS images with the file type designation DRE. Projections of an event cube generated using GammaLib for viewing period 2.06 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.

Table 1

Standard event selection parameters that are suitable for most analysis situations and that were used throughout this paper.

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

thumbnail Fig. 2

Comparison of an event cube generated with GammaLib for viewing period 2.0 and the 1–3 MeV energy band (red solid) and the equivalent event cube MPE-DRE-5∅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 cm2 photons−1, is factorised using (3)

where DRX(a, δ) is the exposure in units of cm2 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 gamma-ray 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 so-called deadtime correction factor, specifying the fraction of time during which the telescope is able to register events. For COMPTEL this fraction is rather constant and was determined by van Dijk (1996) to be L/T = 0.965. This value is hardcoded in GammaLib, and automatically applied to all analysis results.

3.4.2 Exposure

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

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

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

3.4.3 Geometry function

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

The computation of the geometry function is implemented in the method GCOMDri::compute_drg and results are stored as three-dimensional FITS images with the file type designation DRG. Figure 4 illustrates the result of the geometry function computation, again in the form of projections for viewing period 2.0. Equivalent projections for a geometry function obtained using the COMPASS system are superimposed. The projections of the datasets are indistinguishable, illustrating that the geometry functions generated by GammaLib are equivalent to those produced by the COMPASS software.

thumbnail Fig. 3

Comparison of an exposure map obtained with GammaLib for viewing period 2.0 (red solid) to the exposure map MPE-DRX-32382 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 (Peff, Eγ) is an efficiency factor that is detailed in Appendix G.1, PKN(, Eγ) is the contribution of the Klein-Nishina cross-section to a given bin in that is detailed in Appendix G.2, and RD1(E1| Ê1) and RD2(E2Ê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 (E1) 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 gamma-ray photons. However, the IAQ is assumed independent of θ, and hence the kernel will be computed for an average zenith angle of θ = 25° (see Appendix G.4 for details).

The Compton scattering probabilities are stored as a two-dimensional FITS image with the file type designation IAQ. For illustration, Fig. 5 compares projections of the Compton scattering probabilities obtained by GammaLib for the energy band 1–3 MeV to the projections for an equivalent file that was obtained using COMPASS. Both projections are indistinguishable, and the same is true for the other energy bands that were investigated. As demonstrated, the GammaLib implementation accurately reproduces the response computations that were implemented in COMPASS. We note that the implementation of the COMPTEL response computation within GammaLib is crucial for the preservation of COMPTEL data analysis capabilities, since response functions are lacking in the HEASARC archive. In addition to the internal response computation, GammaLib also includes several response functions that were derived using simulations (Stacy et al. 1996), and that can be used as an alternative to the internally computed response functions.

thumbnail Fig. 4

Comparison of a geometry function obtained with GammaLib for viewing period 2.0 (red solid) to the geometry function MPE-DRG-35128 obtained with COMPASS (blue dashed). The left panel shows the distribution of χ 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 ψ.

thumbnail Fig. 5

Comparison of projections of an IAQ obtained with GammaLib for the energy band 1–3 MeV (red solid) to projections for ROL-IAQ-75 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 three-dimensional COMPTEL data space (e.g. Bloemen et al. 1994; Knödlseder 1994; van Dijk 1996; Weidenspointner et al. 2001). Satisfactory, although not perfect, results were obtained using the so-called SRCLIX method that was developed by Bloemen et al. (1994) and that has undergone several evolutions (van Dijk 1996). The SRCLIX method iteratively applies the BGDLIX algorithm, of which we implemented two variants in GammaLib. For reference, we also implemented the PHINOR algorithm (described below), and all background modelling methods are invoked in GammaLib via the GCOMObservation::compute_drb method. The iterative SRCLIX method is implemented as a ctools script (cf. Sect. 3.6).

3.5.1 PHINOR

The motivation for the PHINOR algorithm is the observation that the ratio of DRE and DRG multiplied by the solid angle Ω(χ, ψ) of the event cube bins is to first order independent of χ and ψ. This leads to a background model given by (10)

where (11)

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 large-scale residuals are frequently observed (van Dijk 1996).

3.5.2 BGDLIXA

The BGDLIXA algorithm applies a correction to the PHINOR background model by adjusting 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 , , Nincl and Nexcl are either odd integers or zero, and ⌊/⌋ designates the integer division operator. is the bin size in of the three-dimensional event cube, which typically is 2°. While in early days Nexcl > 0 was used, later analyses set Nexcl = 0, which simplifies the summation range to . Past COMPTEL analyses generally used Nincl = 13 resulting in . As we will see later, Nincl is the most critical parameter of the BGDLIXA algorithm, and specifies the number of layers of the three-dimensional event cube over which the templates are adjusted. If Nincl is large the fraction in Eq. (12) varies little with , leading to a background model that essentially follows the templates. For small values of Nincl 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 running-average summation over a small range in χ and ψ, with Rx = [χ – Δχrunav,χ + Δχrunav] and Rψ = [ψ – Δψrunav,ψ + Δψrunav] where Δχrunav = Δχ Nranav and Δψrunav = Δψ Nrunav. Here, 2Nrunav + 1 specifies the number of χ and ψ bins over which the running-average summation is performed, with Δχ and Δψ being the bin size in χ and ψ of the three-dimensional event cube, and Nrunav is an integer number. Usually, Δχ = Δψ = 1°, and past analyses employed Nrunav = 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 = Δχ ⌊(Navgr – 1)/2⌋ and Δψavgr = Δψ ⌊(Navgr – 1)/2⌋.

Here, Navgr is an odd integer that specifies the number of χ and ψ bins over which the adjustment is performed, with Δχ and Δψ being the bin size inχ and ψ of the three-dimensional event cube. Past analyses usually employed Navgr = 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 three-dimensional event cube where DRG values are small. Using Eq. (14) avoids such problems, and in fact leads to a simplified background modelling algorithm that we implemented in GammaLib under the acronym BGDLIXE (see below). Furthermore, we added a global φ normalisation step, (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 solid-angle-weighted 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 three-dimensional event distribution after subtracting any celestial sources. We note that Nrunav 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 Navgr. 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 Nranav = 3, Navgr = 3, Nincl = 13 and Nexcl = 0 proposed by van Dijk (1996). While the PHINOR model provides a first-order description of the event distribution, it shows significant differences in the details. The BGDLIXA and BGDLIXE models follow the event distribution very closely, and both models are in fact indistinguishable from each other, demonstrating that their results are equivalent for the parameter values chosen.

3.6 ctools scripts

To support the science analysis of COMPTEL data, we extended the ctools package with a number of dedicated Python scripts, providing basic building blocks that each perform well-defined science data analysis tasks. This includes the generation of a database based on the data available in the HEASARC archive, the selection of COMPTEL viewing periods from this database, event binning and response computation, data combination, background modelling, model fitting, generation of test statistic (TS) maps, residual inspection, observation simulations, and source detection, as well as generation of pulsar pulse profiles. These scripts complement the already existing generic science analysis tools in ctools that may be used in combination with the new scripts. Table 2 summarises the COMPTEL-specific scripts that we added to ctools.

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

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

Individual viewing periods can be combined into single event cubes for each energy band {E} using the comobsadd script, leading to a considerable speed-up of the data analysis. The combination of viewing periods is however not required, and alternatively the data can be analysed using a joint maximum-likelihood 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 Ti is the exposure time of viewing period i and T = Σi; Ti is the total exposure time of the combined data. We note that Eq. (20) leads to an exposure map that is independent of α, δ and that is given by the maximum of the DRX for each viewing period, an approximation that is justified by the fact that the DRX is rather flat (see Fig. 3) and that the zenith angle variation in the response computation (cf. Eq. (3)) is dominated by the geometry function DRG (see also Knödlseder 1994).

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

The main script for maximum-likelihood model fitting is comlixfit, which implements the iterative SRCLIX algorithm. The algorithm starts from an input model definition XML file and computes a DRM model of the celestial sources using Eq. (11), which is then used to compute an initial background model using one of the PHINOR, BGDLIXA, or BGDLIXE algorithms. The script then uses the ctlike tool to fit the source and background models to the binned data. The celestial source model that results from the fit is then used to update the DRM model and to regenerate a background model using the selected algorithm. A new ctlike fit is then performed using the updated model, and the procedure is repeated until the maximum log-likelihood change between subsequent iterations becomes negligible, typically less than 0.05. Usually, the SRCLIX algorithm converges after a few iterations.

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

(Mattox et al. 1996), where ln L(Ms + Mb) is the log-likelihood value obtained when fitting the source and the background models together to the data, and ln L(Mb) is the log-likelihood value obtained when fitting only the background model to the data. Under the hypothesis that the background model Mb 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 (p-value) that the log-likelihood improves by TS/2 when adding the source model Ms due to statistical fluctuations only (Cash 1979).

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

and (25)

where φgeo is the angular distance between a sky position (α, δ) and the Compton scatter direction (χ, ψ), and [ARMmin, ARMmax] defines a selection window for the so-called angular resolution measure that is typically taken to be [−3°, 3°]. We note that DRM designates here the model cube that comprises both the source and background components. By default, comobsres uses (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°].

thumbnail Fig. 6

Comparison of PHINOR (blue dotted), BGDLIXA (dark green dashed), and BGDLIXE (light green dashed-dotted) background model for viewing period 1.0 and the energy band 1–3 MeV to the DRE event distribution (red solid). The left panel shows the distribution of χ 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 ψ.

Table 2

COMPTEL-specific ctools scripts.

thumbnail Fig. 7

Test statistic maps of viewing period 1.0 (during which the Crab pulsar and pulsar wind nebula was observed) for the combined analysis of the four standard COMPTEL energy bands as a function of the number of SRCLIX iterations, indicated in the upper-right corner of each panel. The BGDLIXE algorithm with Navgr = 5, Nincl = 15, and Nexcl = 0 was employed. The lower-right 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 gamma-ray sources to the measured events. As is obvious from Fig. 6, the PHINOR model certainly does not satisfy the first criterion; hence, we no longer consider this model here. On the other hand, the BGDLIXA and BGDLIXE models look promising, and were generally used in the past for science analysis of COMPTEL data. Since the BGDLIXA and BGDLIXE algorithms are equivalent as long as NrunavNavgr, we limit our study here to the simpler BGDLIXE algorithm. Specifically, we investigate which values of Navgr and Nincl provide reliable background models without introducing a significant bias in the reconstruction of celestial gamma-ray source characteristics. In agreement with previous studies of the algorithms (cf. van Dijk 1996) we always set Nexcl = 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 Navgr and Nincl. For each viewing period and energy band we created residual maps and histograms using comobsres that we inspected visually. We note that due to the dominance of the instrumental background in COMPTEL data, we do not expect to see celestial gamma-ray emission in the residual maps of individual viewing periods, with the exception of emission from the Crab pulsar and pulsar wind nebula, which is the strongest gamma-ray source at MeV energies. We find a general trend of stronger residuals at lower energies, with the largest residuals observed in the 1–3 MeV band. Notably, residuals are stronger for viewing periods for which the Earth horizon selection Eq. (2) introduces a strong dependence in the χ ψ distribution of the events. The amplitude of the residuals changes strongly under variation of Nincl, with larger residuals for larger values of Nincl, while variations of Navgr impact the residuals only moderately.

For illustration we show residual maps obtained using the BGDLIXE algorithm for different values of Navgr and Nincl and for three representative viewing periods and the 1–3 MeV energy band in Figs. 8 and 9. Viewing period 1.0 is an observation of the Crab pulsar and pulsar wind nebula, which is clearly visible in the residuals maps. Viewing period 223.0 is an observation of 1E 1740.7–2942, a low-mass X-ray binary that is also known as the ‘Great Annihilator’ and that is situated near the Galactic centre. The residuals in this viewing period are relatively modest. Viewing period 518.5 is an observation of the BL Lacertae object S5 0716+714, which is among the viewing periods with the worst residuals, featuring large zones of significant excess counts and negative depressions. The amplitude of the residuals clearly decreases with decreasing value of Nincl, while at the same time smaller value of Nincl 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 Nincl tend to lead to an underestimation of source fluxes. Therefore, the choice of Nincl is necessarily a trade-off between the amplitude of background residuals and the suppression of source flux. On the other hand, the amplitude of the residuals changes only moderately with Navgr, with a slight increase of the residual amplitudes for increasing values of Navgr.

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 Nincl for Navgr = 5 and as a function of Navgr for Nincl = 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 Nincl 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 Navgr so that more events get included in the χ, ψ averaging, as indicated in the lower row of Fig. 10. On the other hand, increasing Navgr leads to a slight increase of the mean and rms distributions; hence, the selected value of Navgr should also not be too large.

thumbnail Fig. 8

Residual maps for three viewing periods and the 1–3 MeV energy band with background modelled using the BGDLIXE algorithm for Navgr = 5 and different values of Nincl. The colour scale is limited to the range − (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 power-law component with an energy flux of 1.068 × 10−8 erg cm−2 s−1 within 0.75–30 MeV and a spectral index of 2.1, which roughly corresponds to the spectral parameters that are observed for the Crab. The simulated data of each viewing period were then fit jointly for the four energy bands using 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 pi is the fitted value, vt 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(pi) follows a Gaussian distribution with a mean of zero and a standard deviation of σ = 1.

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

Figure 11 indicates that large values of Nincl 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 Nincl, as illustrated in Fig. 8, which impact the reconstructed source parameters. Reducing Nincl brings the pull distributions more in line with the expectations, yet for Nincl ≤ 13 the median pull of the fitted energy flux drops below zero, indicating a systematic bias towards too low fluxes. Specifically, for Nincl = 5, where background residuals are very small (cf. Fig. 8), the flux reconstruction is strongly biassed, resulting is a significant underestimation of gamma-ray fluxes. The optimum value for Nincl is hence a trade-off between reduction of background residuals and biassing the flux determination.

Flux biassing is also observed with increasing value of Navgr, with a minimum bias that occurs for Navgr = 5. Using hence Navgr = 5 and Nincl = 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 Navgr = 5 and Nincl = 15 is about twice as large as expected from statistical fluctuations only. In other words, when analysing individual viewing periods using BGDLIXE parameters Navgr = 5 and Nincl = 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 Navgr = 5 and Nincl = 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 Nincl = 29 for analyses in our paper for which a binning of 1° is applied in ψ.

thumbnail Fig. 9

Same as Fig. 8 but for Nincl = 13 and different values of Navgr.

thumbnail 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 Nincl for Navgr = 5 in the upper row and as a function of Navgr for Nincl = 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 Gamma-ray emission from the Crab

We began the science validation of our software with a spectral analysis of the gamma-ray emission from the Crab pulsar and pulsar wind nebula, which is the brightest source of gamma rays 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).

thumbnail Fig. 11

Violin plots of the fitted parameter pull distributions for a simulated source at a 20° off-axis angle for all viewing periods except those that include the Crab. The upper row shows results as a function of Nincl for Navgr = 5; the lower row shows results as a function of Navgr for Nincl = 13. For all simulations the source was simulated as a spatially extended disk with a 3° radius and with a power-law spectrum of index Γ = 2.1. Horizontal black bars show the maximum, the median, and the minimum values, respectively. Semi-transparent violins with grey borders indicate the expected pull distributions from statistical fluctuations only.

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 power-law, exponentially cut-off power-law or curved power-law spectral models. We fitted the data jointly for the 30 energy bins using comlixfit with BGDLIXE background model parameters Navgr = 5 and Nincl = 29. The best fit was obtained using the exponentially cut-off power law (29)

with k = (1.81 ± 0.06) × 10−4 ph cm−2 s−1 MeV−1, Γ = 2.00 ± 0.03, and Ec = 39.1 ±9.7 MeV. van der Meulen et al. (1998) do not quote the parameters of a fitted spectral model to the total Crab emission data, so we used the GammaLib multi-wavelength interface to adjust the same spectral models to the data of their Table 2, which also favoured the exponentially cut-off power law with best fitting parameters k = (1.70 ± 0.05) × 10−4 ph cm−2 s−1 MeV−1, Γ = 1.97 ± 0.02, and Ec = 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 energy-dependent ToF correction factor (cf. Table C.1), an energy-independent deadtime correction factor of 0.965 as well as an energy-independent flux correction factor that was eventually applied to SRCLIX analyses to correct for a flux suppression that arose from the modification of the instrument response function (van Dijk 1996). We recall that GammaLib automatically applies the ToF and deadtime correction factors to the results. Whether or not such correction factors were applied by van der Meulen et al. (1998) is not known, yet they may plausibly explain the 6% discrepancy observed between the analyses.

thumbnail 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 cut-off power-law 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. dat8 that is part of the reference data of the X-ray Timing Explorer (XTE) module of the HEASoft software (version 6.29)9.

We first used compulbin with an angular resolution measure of ±3.5° to produce pulse profiles for the four standard COMPTEL energy bands, as displayed in Fig. 2 of Kuiper et al. (2001). The results of this analysis are shown in Fig. 13, on which we overlay for comparison the pulse profiles obtained by Kuiper et al. (2001). Since the Kuiper et al. (2001) profiles were obtained for a larger dataset and angular resolution measures that were not specified in their paper, we scaled the profiles so that the minimum and maximum number of events in the profiles matches the numbers that we obtained in our analysis. We also applied a phase shift of −0.03 to the pulse profiles of Kuiper et al. (2001) to match them to our profiles. While we do not know the origin of this small discrepancy in the pulse phases, we note that a phase shift of −0.03 corresponds to a difference of about 0.5 arcsec in the assumed right ascension of the Crab pulsar. In GammaLib, the pulsar position is taken from the ephemerides file, which in the present case is the radio position in the Princeton database, while Kuiper et al. (2001) do not specify the position that they assumed for the Crab pulsar. We note that the radio position in the Princeton database differs by about 0.5 arcsec from the International Celestial Reference System (ICRS) position provided by the Set of Identifications, Measurements, and Bibliography for Astronomical Data (SIMBAD) service10, 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 Navgr = 5 and Nincl = 29. We used two point-source model components in our fit, one for the Crab nebula that was fitted to the data of both phase intervals, and one for the Crab pulsar that was only fitted to the data of the Total Pulse interval. Consequently, the Crab pulsar component modelled only events that were in excess of the pulsar wind nebula component. Both point-source model components were located at the position of the Crab pulsar, as defined above, and had a spectral model with a free flux value for each of the 12 energy bins. Similar to Kuiper et al. (2001), we used simulated instrument response functions and 50 bins in ψ with bin sizes of 1° for our analysis.

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

We also fitted different spectral models to the data of the Crab pulsar and pulsar wind nebula, including power laws, exponentially cut-off power laws and curved power laws. We determined the corresponding uncertainty bands using ctbutterfly and overlay them on the spectral points in Fig. 14. Using an exponentially cut-off power law for the nebula component instead of a simple power law improved the TS value of the nebula component by 6.5, corresponding to a detection significance of 2.5σ 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 Ec = 53.3 ± 15.5 MeV. For the Crab pulsar the best fitting power-law parameters were k = (5.1 ± 0.3) × 10−5 ph cm−2 s−1 MeV−1 and Γ = 2.29 ± 0.06. This compares to the spectral indices of 2.227 ± 0.013 and 2.24 ± 0.04 determined by Kuiper et al. (2001) for the nebula and pulsar components using power-law models, respectively. While our index for the pulsar component is compatible with their result, our index for the nebula component is flatter, which can be explained by the spectral cutoff in our model. Using a simple power law for the nebula component, as Kuiper et al. (2001) did, we obtained a steeper index of 2.24 ± 0.02 that is compatible with their result.

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

thumbnail 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 Phase-resolved analysis of LS 5039

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

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

As the first step we created TS maps of the region around LS 5039 using comlixmap for the INFC and SUPC phase intervals. The data for the four standard energy bands were analysed jointly, using a model composed of a test point source and an additional point source at the location of the quasar PKS 1830–210 that is spatially close to LS 5039. The spectra of both components were modelled using power laws. In addition, the source model included components for Galactic diffuse emission based on template maps for bremsstrahlung and inverse Compton emission with free scaling factors for each energy bin11. Furthermore, an isotropic component was included to account for the cosmic gamma-ray background emission, with intensity fixed according to I(Eγ) = 1.12 × 10−4(Eγ/5 MeV)−2.2 ph cm−2 s−1 MeV−1 sr−1 as suggested by Weidenspointner (1999). The background was modelled using the BGDLIXE algorithm with parameters Navgr = 5 and Nincl = 15.

The resulting TS maps are shown in Fig. 15 for the INFC (left) and SUPC (right) phase intervals. The maps can be compared to those shown in Fig. 6 of Collmar & Zhang (2014), which were determined separately for the 3–10 MeV and 10–30 MeV energy bands. In both analyses, LS 5039 is considerably more significant in the INFC phase interval but only weakly detected in the SUPC phase interval. In the latter interval, the emission maximum seems to be displaced towards the north-east with respect to the position of LS 5039 in ours and the analysis of Collmar & Zhang (2014), yet the emission location is still compatible within the 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 power-law model using ctbutterfly. The results are shown in Fig. 16. For the left panel we used exactly the same event selection that was used by Collmar & Zhang (2014) which is a minimum distance from the Earth horizon of ζmin = 0° and the use of circular exclusion regions to handle D2 modules with failed PMTs (i.e. fpmtflag = 2; cf. Appendix E). Only with this event selection we were able to reproduce the 1-3 MeV flux point of Collmar & Zhang (2014), while using our standard setting of ζmin = 5° and fpmtflag = 0 that excludes D2 modules with failed PMTs produced a notable variation of the 1–3 MeV flux point between inferior and superior conjunction. To verify that this variation can indeed be attributed to differences in the event selection, we also did an equivalent analysis with COMPASS using ζmin = 5° and fpmtflag = 0. This resulted in a 1–3 MeV flux variation between inferior and superior conjunction that was similar to that observed in our analysis, confirming that the spectral differences are due to differences in the event selection. The corresponding results are summarised in the right panel of Fig. 16.

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

Finally, we derived the orbital flux variation in the 10–30 MeV energy band for the two different event selections to reproduce Fig. 8 of Collmar & Zhang (2014). For this purpose we split the 10–30 MeV data into five phase intervals of equal length and fitted the data using comlixfit with a source model comprising components for LS 5039, PKS 1830-210, bremsstrahlung emission, inverse Compton emission and cosmic background emission (see above for details). The background was modelled using the BGDLIXE algorithm with parameters Navgr = 5 and Nincl = 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.

thumbnail 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 1830-210 was included in the source model and hence is not visible in the maps.

thumbnail Fig. 16

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

thumbnail Fig. 17

Orbital light curve of LS 5039 in the 10-30 MeV energy band. Data points for the ctools analysis were displaced by ±0.02 in phase for clarity. Analysis results obtained under analysis conditions identical to those in Collmar & Zhang (2014) (ζ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.

Table 3

Energy bins for 1.8 MeV gamma-ray line analysis.

4.4 26AI line emission from Carina

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

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

We analysed the same data and adopted the same data space definition that was used by Knödlseder et al. (1996a), yet we tested an alternative analysis method that jointly handles the 26Al line signal and any underlying continuum emission. This is more in line with the GammaLib and ctools philosophy of explicitly modelling all emission components and provides a more accurate handling of underlying continuum emission. Specifically, we split the data within the 1–3 MeV energy interval into 11 energy bins, specified in Table 3, and analysed them jointly using a combination of model components describing the 1.8 MeV line signal, any underlying continuum gamma-ray emission and the instrumental background. To model the 1.8 MeV line emission spectrum we used a Gaussian spectral component with a fixed mean of 1.809 MeV and a standard deviation of σ = 58.9 keV that corresponds to COMPTEL’s instrumental energy resolution at that energy. Similarly to our analysis of LS 5039 we modelled continuum emission using template maps for Galactic bremsstrahlung and inverse Compton emission and an isotropic component for the cosmic gamma-ray background. The spectra of the three continuum components were model using power laws, where the prefactors and indices were free for the Galactic components, while the prefactor and index were fixed to I(Eγ) = 1.12 × 10−4(Eγ/5 MeV)−2.2ph cm−2 s−1 MeV−1 sr−1 for the cosmic gamma-ray background component, as suggested by Weidenspointner (1999). The TS map was generated using comlixmap and model fitting was done using comlixfit with the standard BGDLIX parameters Navgr = 5 and Nincl = 15.

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

In their Fig. 2, Knödlseder et al. (1996a) present a 1.8 MeV line emission maximum likelihood map of the Carina region, and we produced an equivalent TS map with the same spatial binning using comlixmap. Specifically, we fitted point-source models for the 1.8 MeV line emission and the 1–3 MeV continuum emission for a grid of source positions, producing hence TS maps for both emission components. In the fitting the fluxes of both point sources were constrained to non-negative values. The resulting maps are shown in Fig. 19, where the left map can be compared to the maximum likelihood map presented in Fig. 2 of Knödlseder et al. (1996a). Both maps show qualitatively comparable features, yet we obtained a maximum TS value of 23.1 at (l, b) = (285.5°, 1.5°) while Knödlseder et al. (1996a) found a lower maximum TS value of 14.7 at (l, b) = (286.5°, 0.5°). Eventually, these differences may be explained by the background modelling techniques and analysis methods that differ significantly between the studies. We recall that Knödlseder et al. (1996a) analyse the 1.8 MeV line data in a single energy bin covering 1.7–1.9 MeV and used a background model derived from adjacent energy bands, which to first order includes also continuum emission, but which does not properly account for spectral differences between instrumental background and diffuse emission components as well as differences in their distributions (Bloemen et al. 1999).

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

As the next step we fitted the data with a source model composed of a 1.8 MeV line component modelled using a point source and a 1–3 MeV continuum component modelled using a combination of bremsstrahlung and inverse Compton spatial maps as well as an isotropic component for the cosmic gamma-ray background. The position of the point-source model as well as the spectral parameters for the Galactic continuum power-law components were free parameters in the fit. Fitting the data using comlixfit gave a best fitting position of (l, b) = (285.4° ± 0.8°, 1.3° ± 0.7°), a flux of (3.1 ± 0.6) × 10−5 ph cm−2 s−1 and a TS value of 22.3 for the 1.8 MeV line component. Our 1.8 MeV line flux is consistent with the value of (3.1 ±0.8) × 10−5 ph cm−2 s−1 found by Knödlseder et al. (1996a), yet our best fitting position is offset by about 1.4° from the one found in 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 log-likelihood value and the fitted 1.8 MeV line flux as a function of spatial disk radius. For that purpose, we determined the maximum likelihood solution for a disk model with position fixed at our maximum likelihood solution (l, b) = (285.4°, 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 log-likelihood 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.

thumbnail Fig. 18

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

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

thumbnail Fig. 20

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

5 Conclusion

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

We want to stress that our work was also motivated by the goal of reducing the carbon footprint of astronomical research. As recently pointed out by Knödlseder et al. (2022), the current deployment rate of new astronomical observatories is not compatible with the imperative of reducing the carbon footprint across all activity sectors of modern societies, and this calls for fundamental changes in astronomical practices in the future. Among the many possible alternatives to the building of ever more and ever bigger new astronomical facilities is the exploitation of archival data from past missions that may have scientific treasures yet to be discovered. Since version 2.0.0, ctools estimates the carbon footprint of its use based on the assessment of Berthoud et al. (2020) for the GRICAD computing centre, and, using this feature, we estimate that the work presented in this publication resulted in the generation of 600 ± 300 kgCO2e of greenhouse gases (GHGs) due to the computing related to the analysis of archival data. This is about 40 times less than the median per-publication emissions associated with the analysis of data from an active astronomical observatory (Knödlseder et al. 2022). Taking all sources of emissions related to this work into account, we estimate the carbon footprint of this research to be 1.9 ± 0.3 tCO2e (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 community-developed gamma-ray astronomy science analysis software (Knödlseder et al. 2016). ctools is based on GammaLib, a community-developed toolbox for the scientific analysis of astronomical gamma-ray data (Knödlseder et al. 2011). This work has made use of the Python 2D plotting library matplotlib (Hunter 2007). This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research is part of the #LowCarbon-Science initiative that aims in reducing the carbon footprint of scientific research, and has benefitted from discussions held in the context of the GDR Labos 1point5 and the Astro4Earth initiative.

Appendix A HEASARC archive

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

Table A.1

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 Gam-maLib we found that the units in some of the table columns were not correct. We summarise the inconsistencies that we encountered in Table A.2. We furthermore found that the good time interval dataset for viewing period 8.0 with the identifier MPE-TIM-11481 only contains about half of the good time intervals that were defined in COMPASS, probably as the result of a file truncation error in the input ASCII file that was used for the generation of the FITS file in the HEASARC archive.

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

Table A.2

Incorrect header units in FITS files in the COMPTEL HEASARC archive.

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

Appendix B Time of flight conversion

The HEASARC archive mixes different versions of EVP files that have different levels of processing for the ToF values. The versions can be distinguished by the header keyword DSD_REP in the EVP file, specifying either 2 for ToFII or 3 for ToFIII, 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 ToFIII using (B.1)

with the coefficients given by Tables B.1 and B.1 and the D1 and D2 energy deposits E1 and E2 given in MeV. For the HEASARC archive, ToF values accessed through GammaLib are therefore always ToFIII values.

Table B.1

Coefficients, α;, for the conversion from ToFII to ToFIII.

Table B.1

Coefficients, bi, for the conversion from ToFII to ToFIII.

Appendix C Flux correction due to time of flight selection

To correct for the photon rejection by the ToF selection an energy-dependent correction factor needs to be applied to the instrument response. Since this correction factor depends on the ToF selection interval, it is computed in GammaLib during the 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 [Emin, Emax] for which the event cube is generated. Only corrections for ToFmax = 130 are supported.

Table C.1

Time of flight correction factors for ToFmax = 130.

Appendix D Pulsar timing

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

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

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

After applying the clock correction, event times are converted using (D.1)

where ∆tCGRO→SSB corrects for the light travel time from CGRO to the Solar System barycentre, ∆tShapiro corrects for the gravitational time delay near the Sun, also known as Shapiro delay, ∆tUTC→TT converts from UTC to terrestrial time (TT), and ∆tTT→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 ∆tCGRO→SSB − ∆tShapiro + ∆tTT→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)

and (D.3)

where (D.4)

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 tCGRO (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 tCGRO is computed using (D.8)

where ∆tTT→TDB(i) is the nearest entry in the JPL database.

The last term in the time correction, ∆tUTC→TT, does not depend on planetary ephemerides and is given by (D.9)

where nleap is the number of leap seconds. The computation of ∆tUTC→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 ∆tUTC→TT + ∆tTT→TBD, avoiding the need for planetary ephemerides. These so-called BVC data can be handled by GammaLib through the classes GCOMBvc and GCOMBvcs that manage individual data records as well as entire files. Specifically, the method GCOMBvcs::tdelta computes ∆tCGRO→SSB − ∆tShapiro + ∆tUTC→TT + ∆ttt→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 = tSSBt0 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)

where (D.13)

and MJD0 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.

Table E.1

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)

with (F.2)

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 Gi(χ, ψ) 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)

where (F.4)

and (F.5)

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

with r1 = 13.8 cm being the radius of a D1 module and r2 = 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 dkl is the projected distance between the centres of the D1 and D2 modules, given by (F.7)

where xk and yk are the geometric positions of the D1 modules and xl and yl 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)

and (F.10)

The term fk1e(θ, ϕ) 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 xe and ye and a radius re as given in Table E.1. Specifically, for fpmtflag = 2 (F.11)

with (F.12)

being the distance between the projected D1 module centres xk and yk and the centre xe and ye 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)

with (F.14)

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)

and (F.17)

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 Peff(φgeo, Eγ) in Eq. (6) is factorised according to (G.1)

with (G.2)

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.

Table G.1

GammaLib methods that implement the computation of the efficiency factors.

The PA1(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 energy-dependent interaction coefficient of aluminium in units of cm−1 that is interpolated using a log-log interpolation of the values given in Table G.2 and labove = 0.147 cm is the thickness of the material above D1.

The PV1(Eγ) is the transmission probability for photons for the first Veto dome and is computed using (G.4)

where μVeto(Eγ) is the energy-dependent interaction coefficient for the Veto dome in units of cm−1 that is interpolated using a log-log interpolation of the values given in Table G.2 and lV1 = 1.721 cm is the thickness of the first Veto dome.

PD1(Eγ) is the interaction probability in D1 and is computed using (G.5)

where μD1(Eγ) is the energy-dependent D1 attenuation coefficient in units of cm−1 that is interpolated using a log-log interpolation of the values given in Table G.2 and lD1 = 8.5 cm is the thickness of the D1 modules.

PC(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 ≤ PCompton(Eγ) ≤ 1.

The PMH(Eγ) is the probability that there is no multi-hit. This probability is computed using a log-log interpolation of the values of Table G.2. PSV(Eγ) is the probability that the photon was not self-vetoed. This probability is computed using a linear interpolation of the values of Table G.2.

The PPSD(Ê1) is the D1 energy-dependent PSD correction and is given by (G.7)

where Ê1 is in units of MeV. This correction applies to a standard PSD selection of 0-110.

PA2(Ê2) gives the transmission probability of the aluminium below the D1 modules and is computed using (G.8)

where μA1(Ê2) is the energy-dependent interaction coefficient for aluminium in units of cm-1 that is interpolated using a loglog interpolation of the values given in Table G.2 and lbetween = 0.574 cm is the thickness of the aluminium plate.

PV23(Ê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 energy-dependent interaction coefficient of the Veto domes in units of cm-1 that is interpolated using a log-log interpolation of the values given in Table G.2 and lV23 = 3.442 cm is the combined thickness of the second and third veto domes.

PD2(Eγ) gives the interaction probability in D2 and is computed using (G.10)

where μD2(Ê2) is the energy-dependent D2 attenuation coefficient in units of cm−1 that is interpolated using a log-log interpolation of the values given in Table G.2 and lD2 = 7.525 cm is the thickness of the D2 modules.

The PMS(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 rD1 = 13.8 cm and lD1 = 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

Table G.2

Efficiency parameters as a function of energy.

thumbnail Fig. G.1

Energy dependence of the efficiency factor Peff(φ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).

(G.12)

with (G.13)

Figure G.1 illustrates the energy dependence of the efficiency factor Peff(φgeo, Eγ) and its components.

Appendix G.2 Klein-Nishina contribution

(G.14)

is the contribution of the Klein-Nishina cross-section σKN(φgeo, Eγ) to the φgeo bin spanned by [φgeo,min, φgeo,max] where the integrals are computed using (G.15)

with (G.16)

Appendix G.3 Module spectral response functions

Appendix G.3.1 D1 spectral response

The spectral response RD1(E1|Ê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). T1(E1 |Ê1) is a threshold function, defined by (G.18)

where . μp(Ê1) and σ(Ê1) as well as and ∆E1(Ê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 RD1(E1|Ê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.

thumbnail Fig. G.2

Detector response functions RD1(E1|Ê1) (left) and RD2(E2|Ê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 RD2(E2|Ê2) provides the probability for measuring an energy E2 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 T2(E2|Ê2) (Diehl et al. 1992). Here, Ê2 is the energy of the gamma-ray photon incident on a D2 module, and E2 is the measured energy deposit.

The first Gaussian, (G.20)

describes the so-called photo-peak which represents photons that were fully absorbed in the D2 module We note that the position μPP(Ê2) of the photo-peak 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 photo-peak energy. The escape peaks are modelled using (G.21)

and (G.22)

with μSe(Ê2) = μpp(Ê2) − mec2 and μde(Ê2) = μpp(Ê2) − 2mec2, where mec2 ≈ 511 keV is the rest energy of the electron.

The first continuum component, (G.23)

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

and the probability that the Compton scattered photon does not undergo a second Compton scattering before escaping, empirically modelled using (G.25)

where (G.26)

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 σ(E2) is the standard deviation of the detector response at the reconstructed energy E2.

The second continuum component, (G.29)

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

The threshold function is defined by (G.30)

where and

All five response components have amplitude parameters (App, Ase, Ade, Acom and Abkg) that are tabulated as a function of incident energy E2 in a calibration file of type SDB. In addition, this file contains also tabulated values for the photo-peak energy μPP(Ê2), the standard deviation σ(Ê2) and the minimum and maximum threshold energies and and width ∆E2(Ê2). Values for a given energy Ê2 are obtained by linear interpolation of the tabulated values.

The spectral response RD2(E2|Ê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 z-direction, and h = 158 cm is the distance between the D1 and D2 modules. The functions (G.35)

and (G.36)

are the geometrical relations that transform from the uncertainties within the modules to uncertainties in φgeo for a given zenith angle θ.

thumbnail 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° (dashed-dotted 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.

Table H.1

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 tCO2e, composed of 1.2 ± 0.2 tCO2e (65%) for running the office building in which the work was done, 0.4 ± 0.2 tCO2e (21%) for the computing infrastructure and data flow, and 0.3 ± 0.1 tCO2e (14%) for lunch meals and home-to-office commuting.

We based our assessment for computing on the carbon footprint estimation implemented in ctools, yet when we started this research the relevant functionality was not yet available. We note, however, that the carbon footprint estimate for the final analyses for this publication amounted to 60 kgCO2e, of which 20 kgCO2e were attributed to power consumption and 40 kgCO2e 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 kgCO2e, of which 200 kgCO2e were due to power consumption and 400 kgCO2e due to the computing infrastructure. To avoid double counting, we disregarded the contribution of power consumption since it is already accounted for in the electricity to power the office building.

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

References

  1. Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 [Google Scholar]
  2. 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]
  3. Bloemen, H., Diehl, R., Hermsen, W., Knödlseder, J., & Schönfelder, V. 1999, Astrophys. Lett. Commun., 38, 391 [NASA ADS] [Google Scholar]
  4. Bloemen, H., Hermsen, W., Swanenburg, B. N., et al. 1994, ApJS, 92, 419 [NASA ADS] [CrossRef] [Google Scholar]
  5. Casares, J., Ribó, M., Ribas, I., et al. 2005, MNRAS, 364, 899 [Google Scholar]
  6. Cash, W. 1979, ApJ, 228, 93 [Google Scholar]
  7. Coleman, A., & McConnell, M. 2020, in American Astronomical Society Meeting Abstracts, 235, 206.04 [NASA ADS] [Google Scholar]
  8. Collmar, W., & Zhang, S. 2014, A&A, 565, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. 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]
  10. 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]
  11. Diehl, R., Aarts, H., Bennett, K., et al. 1992, in Data Analysis in Astronomy, 201 [Google Scholar]
  12. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  13. IPBES 2019, Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, eds. E. S. Brondizio, J. Settele, S. Díaz, & H. T. Ngo (IPBES secretariat, Bonn, Germany) [Google Scholar]
  14. IPCC 2021, Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, eds. V. Masson-Delmotte, P. Zhai, A. Pirani, et al. (Cambridge University Press) [Google Scholar]
  15. 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]
  16. Knödlseder, J. 1994, PhD thesis, Technical University Munich, Germany [Google Scholar]
  17. Knödlseder, J., Bennett, K., Bloemen, H., et al. 1996a, A&As, 120, 327 [Google Scholar]
  18. Knödlseder, J., von Ballmoos, P., Diehl, R., et al. 1996b, SPIE Conf. Ser., 2806, 386 [Google Scholar]
  19. Knödlseder, J., Mayer, M., Deil, C., et al. 2016, Astrophysics Source Code Library [ascl:1601.005] [Google Scholar]
  20. Knödlseder, J., Mayer, M., Deil, C., et al. 2011, Astrophysics Source Code Library [ascl:1110.007] [Google Scholar]
  21. Knödlseder, J., Brau-Nogué, S., Coriat, M., et al. 2022, Nat. Astron., 6, 503 [CrossRef] [Google Scholar]
  22. Kuiper, L., Hermsen, W., Cusumano, G., et al. 2001, A&A, 378, 918 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Landrigan, P., Fuller, R., Acosta, N., et al. 2017, Lancet, 391 [Google Scholar]
  24. Martin, P., Brau-Nogué, S., Coriat, M., et al. 2022, ArXiv e-prints, [arXiv:2204.12362] [Google Scholar]
  25. Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 [Google Scholar]
  26. Much, R., Bennett, K., Buccheri, R., et al. 1995a, A&A, 299, 435 [NASA ADS] [Google Scholar]
  27. Much, R., Bennett, K., Buccheri, R., et al. 1995b, Adv. Space Res., 15, 81 [NASA ADS] [CrossRef] [Google Scholar]
  28. Much, R., Harmon, B. A., Nolan, P., et al. 1996, A&As, 120, 703 [NASA ADS] [Google Scholar]
  29. 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]
  30. Schönfelder, V., Aarts, H., Bennett, K., et al. 1993, ApJS, 86, 657 [NASA ADS] [CrossRef] [Google Scholar]
  31. Stacy, J. G., Kippen, R. M., Kappadath, S. C., et al. 1996, A&As, 120, 691 [NASA ADS] [Google Scholar]
  32. Strong, A., & Collmar, W. 2019, Mem. Soc. Astron. Italiana, 90, 297 [Google Scholar]
  33. Strong, A. W., Bennett, K., Bloemen, H., et al. 1996, A&As, 120, 381 [NASA ADS] [Google Scholar]
  34. van der Meulen, R. D., Bloemen, H., Bennett, K.,et al. 1998, A&A, 330, 321 [NASA ADS] [Google Scholar]
  35. van Dijk, R. C. A. 1996, PhD thesis, University of Amsterdam, The Netherlands [Google Scholar]
  36. Weidenspointner, G. 1999, PhD thesis, Technical University Munich, Germany [Google Scholar]
  37. Weidenspointner, G., Varendorff, M., Oberlack, U., et al. 2001, A&A, 368, 347 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

3

Throughout this paper we denote by instrumental background all events that are not related to celestial gamma-ray emission.

4

Some viewing periods in the HEASARC archive are unreadable, and data from the end of the mission were not archived. Details about issues encountered in the HEASARC COMPTEL data archive are provided in Appendix A.

5

We note that COMPTEL times are specified by two values: the truncated Julian days (TJD), which is defined as the number of modified Julian days (MJD) minus 40 000, and the COMPTEL tics, which are 1/8 milliseconds long. As an example, TJD = 8393 and zero tics corresponds to 1991–05–17T00:00:00 UT.

6

We have chosen the second viewing period in the HEASARC archive to illustrate the data cube comparisons in our paper since for the first viewing period we encountered percent-level differences that are plausibly attributed to differences in our and the COMPASS processing configuration that we were not able to track down.

7

The TS map comprises 50 × 50 pixels of size 1° × 1° situated around the Crab position.

9

The HEASoft software can be downloaded from https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/

11

We used the bremsstrahlung emission map with the COMPASS identifier MPE-MIS-13006 and the inverse Compton emission map with the COMPASS identifier ROL-MIS-163 that were also used by Collmar & Zhang (2014).

All Tables

Table 1

Standard event selection parameters that are suitable for most analysis situations and that were used throughout this paper.

Table 2

COMPTEL-specific ctools scripts.

Table 3

Energy bins for 1.8 MeV gamma-ray line analysis.

Table A.1

Inexploitable viewing periods in the range 1.0 to 719.0 in the COMPTEL HEASARC archive.

Table A.2

Incorrect header units in FITS files in the COMPTEL HEASARC archive.

Table B.1

Coefficients, α;, for the conversion from ToFII to ToFIII.

Table B.1

Coefficients, bi, for the conversion from ToFII to ToFIII.

Table C.1

Time of flight correction factors for ToFmax = 130.

Table E.1

D2 module exclusion circles applied following the failure of PMTs at the specified dates.

Table G.1

GammaLib methods that implement the computation of the efficiency factors.

Table G.2

Efficiency parameters as a function of energy.

Table H.1

Carbon footprint estimate by emission source of the research work behind this paper.

All Figures

thumbnail 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
thumbnail Fig. 2

Comparison of an event cube generated with GammaLib for viewing period 2.0 and the 1–3 MeV energy band (red solid) and the equivalent event cube MPE-DRE-5∅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
thumbnail Fig. 3

Comparison of an exposure map obtained with GammaLib for viewing period 2.0 (red solid) to the exposure map MPE-DRX-32382 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
thumbnail Fig. 4

Comparison of a geometry function obtained with GammaLib for viewing period 2.0 (red solid) to the geometry function MPE-DRG-35128 obtained with COMPASS (blue dashed). The left panel shows the distribution of χ 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
thumbnail Fig. 5

Comparison of projections of an IAQ obtained with GammaLib for the energy band 1–3 MeV (red solid) to projections for ROL-IAQ-75 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
thumbnail Fig. 6

Comparison of PHINOR (blue dotted), BGDLIXA (dark green dashed), and BGDLIXE (light green dashed-dotted) background model for viewing period 1.0 and the energy band 1–3 MeV to the DRE event distribution (red solid). The left panel shows the distribution of χ 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
thumbnail Fig. 7

Test statistic maps of viewing period 1.0 (during which the Crab pulsar and pulsar wind nebula was observed) for the combined analysis of the four standard COMPTEL energy bands as a function of the number of SRCLIX iterations, indicated in the upper-right corner of each panel. The BGDLIXE algorithm with Navgr = 5, Nincl = 15, and Nexcl = 0 was employed. The lower-right 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
thumbnail Fig. 8

Residual maps for three viewing periods and the 1–3 MeV energy band with background modelled using the BGDLIXE algorithm for Navgr = 5 and different values of Nincl. The colour scale is limited to the range − (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
thumbnail Fig. 9

Same as Fig. 8 but for Nincl = 13 and different values of Navgr.

In the text
thumbnail 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 Nincl for Navgr = 5 in the upper row and as a function of Navgr for Nincl = 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
thumbnail Fig. 11

Violin plots of the fitted parameter pull distributions for a simulated source at a 20° off-axis angle for all viewing periods except those that include the Crab. The upper row shows results as a function of Nincl for Navgr = 5; the lower row shows results as a function of Navgr for Nincl = 13. For all simulations the source was simulated as a spatially extended disk with a 3° radius and with a power-law spectrum of index Γ = 2.1. Horizontal black bars show the maximum, the median, and the minimum values, respectively. Semi-transparent violins with grey borders indicate the expected pull distributions from statistical fluctuations only.

In the text
thumbnail 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 cut-off power-law models.

In the text
thumbnail 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
thumbnail 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
thumbnail 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 1830-210 was included in the source model and hence is not visible in the maps.

In the text
thumbnail Fig. 16

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

In the text
thumbnail Fig. 17

Orbital light curve of LS 5039 in the 10-30 MeV energy band. Data points for the ctools analysis were displaced by ±0.02 in phase for clarity. Analysis results obtained under analysis conditions identical to those in Collmar & Zhang (2014) (ζ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
thumbnail Fig. 18

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

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

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

In the text
thumbnail Fig. G.1

Energy dependence of the efficiency factor Peff(φ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
thumbnail Fig. G.2

Detector response functions RD1(E1|Ê1) (left) and RD2(E2|Ê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
thumbnail 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° (dashed-dotted green), and 50° (dotted purple).

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.