DmpIRFs and DmpST: DAMPE Instrument Response Functions and Science Tools for Gamma-Ray Data Analysis

TL;DR

This paper introduces the DAMPE instrument response functions and a dedicated software tool, DmpST, to improve gamma-ray data analysis for dark matter and astrophysics research, utilizing high-statistics simulations and likelihood methods.

Contribution

It provides accurate instrument response functions and a specialized analysis software for DAMPE gamma-ray data, enhancing analysis capabilities.

Findings

Validated the instrument response functions with simulations.

Developed DmpST software for efficient data analysis.

Applied maximum likelihood method for source disentanglement.

Abstract

GeV gamma ray is an important observation target of DArk Matter Particle Explorer (DAMPE) for indirect dark matter searching and high energy astrophysics. We present in this work a set of accurate instrument response functions of DAMPE (DmpIRFs) including the effective area, point-spread function and energy dispersion that are crucial for the gamma-ray data analysis based on the high statistics simulation data. A dedicated software named DmpST is developed to facilitate the scientific analyses of DAMPE gamma-ray data. Considering the limited number of photons and the angular resolution of DAMPE, the maximum likelihood method is adopted in the DmpST to better disentangle different source components. The basic mathematics and the framework regarding this software are also introduced in this paper.

Figures

Equations

$$r(E^{\prime },{\hat{p}}^{\prime },s)=\int \int \int F(E,\hat{p})A_{\mathrm{eff}}(E,\hat{v}(t;\hat{p}),s)$$

$$\times P({\hat{v}}^{\prime }(t;{\hat{p}}^{\prime });E,\hat{v}(t;\hat{p}),s)D(E^{\prime };E,\hat{v}(t;\hat{p}),s)\mathrm{d}E\mathrm{d}\Omega \mathrm{d}t,$$

$$A_{\mathrm{eff}}(E_i,{\theta }_j,{\varphi }_k,s)=\frac{N_{i,j,k,s}}{N_{\mathrm{sim},i,j,k}}{A}_{\mathrm{sim}},$$

$$P(x)=f_{\mathrm{core}}K(x_p;{\sigma }_{\mathrm{core}},{\gamma }_{\mathrm{core}})+(1-f_{\mathrm{core}})K(x_p;{\sigma }_{\mathrm{tail}},{\gamma }_{\mathrm{tail}}),$$

$$K(x_p;\sigma ,\gamma )=\frac{1}{2\pi {\sigma }^2}\left(1-\frac{1}{\gamma }\right){\left[1+\frac{1}{2\gamma }\frac{x_p^2}{{\sigma }^2}\right]}^{-\gamma },$$

$$x_p=\frac{\delta v}{S_{\mathrm{p}}(E,\theta )}.$$

$$x_D=\frac{E^{\prime }-E}{S_{\mathrm{D}}(E,\theta )E},$$

$$D(x_D)=\{\begin{array}{ll}{N}_{L}R(x_D,x_0,{\sigma }_L,{\gamma }_L)& \text{if (xD−x0)<−xˉ}\\ N_{l}R(x_D,x_0,{\sigma }_l,{\gamma }_l)& \text{if (xD−x0)∈[−xˉ,0]}\\ N_{R}R(x_D,x_0,{\sigma }_R,{\gamma }_R)& \text{if (xD−x0)>0}\end{array}$$

$$R(x_D,x_0,\sigma ,\gamma )=N\exp \left(-\frac{1}{2}{\left|\frac{x_D-x_0}{\sigma }\right|}^{\gamma }\right).$$

$$ϵ(E,\hat{p})=\sum _s\int A_{\mathrm{eff}}(E,\hat{v}(t,\hat{p}),s)\mathrm{d}t.$$

$$ϵ(E,\hat{p})$$

$$F(E,\hat{p};\vec{\lambda })=S(E;\vec{\lambda })M(\hat{p}).$$

$$\overline{P}(\delta v;E)=\frac{{\sum }_s\int P(\delta v;E,\theta ,s)A_{\mathrm{eff}}(E,\theta ,\varphi ,s)t_{\mathrm{obs}}(\theta ,\varphi )\mathrm{d}\Omega }{{\sum }_s\int A_{\mathrm{eff}}(E,\theta ,\varphi ,s)t_{\mathrm{obs}}(\theta ,\varphi )\mathrm{d}\Omega },$$

$$\overline{D}(E^{\prime };E)=\frac{{\sum }_s\int D(E^{\prime };E,\theta ,s)A_{\mathrm{eff}}(E,\theta ,\varphi ,s)t_{\mathrm{obs}}(\theta ,\varphi )\mathrm{d}\Omega }{{\sum }_s\int A_{\mathrm{eff}}(E,\theta ,\varphi ,s)t_{\mathrm{obs}}(\theta ,\varphi )\mathrm{d}\Omega }.$$

$$r_{ij}(E_i,{\hat{p}}_i^{\prime };{\vec{\lambda }}_j)=\int \mathrm{d}\Omega F_{ij}(E_i,\hat{p};{\vec{\lambda }}_j)ϵ(E_i,\hat{p})\bar{P}({\hat{p}}_i^{\prime };\hat{p},E_i).$$

$$p(\{N_i\};\{R_i\})=\prod _{i=1}^{N_{\mathrm{bins}}}\frac{R_i^{N_i}}{N_i!}\exp (-R_i).$$

$$\log L(\vec{\lambda })$$

$$\log L(\vec{\lambda })=-\int \mathrm{d}t\int \mathrm{d}E{\int }_{\mathrm{ROI}}\mathrm{d}{\Omega }^{\prime }\sum _{j=1}^{N_{\mathrm{s}}}{r}_j({\vec{\lambda }}_j)+\sum _{i=1}^{N_{\mathrm{events}}}\log \sum _{j=1}^{N_{\mathrm{s}}}{r}_j({\vec{\lambda }}_j).$$

$$T{S}_j=-2(\log L({\hat{\vec{\lambda }}}_{0,j})-\log L(\hat{\vec{\lambda }})),$$

Full text

1 Introduction

DArk Matter Particle Explorer (DAMPE) is a high energy cosmic-ray and gamma-ray observatory (Chang 2014, Chang et al. 2017). It contains four sub-detectors: a Plastic Scintillation Detector (PSD), a Sillicon-Tungsten tracKer-converter (STK), a BGO calorimeter (BGO) and a NeUtron Detector (NUD). The PSD that measures the charge of particles also acts as anti-coincidence detector for gamma-ray observation. The STK measures the trajectories of charged particles, as well as the photons that are converted into e*+e-* pairs. The BGO calorimeter measures the energies of incidence particles and is also able to distinguish the electron and hadron efficiently. The NUD provides an independent measurement and further improvement for the electron/hadron identification. The on-orbit calibration have adopted for DAMPE and is expected to operate stably during the next few years (Ambrosi et al. 2019, Ma et al. 2018, Ding et al. 2019, Jiang et al. 2019).

Based on the photon selection algorithm described in Xu et al. (2018), valuable gamma-ray data have been accumulated. Further scientific analysis of high-level gamma-ray data, however, requires detailed knowledge about the instrument response functions (IRFs) of DAMPE, i.e., the effective area, the point-spread function (PSF) and the energy dispersion function. Using the high-statistics simulation data, we have constructed the IRFs for DAMPE gamma-ray observation in the energy range from 1 GeV to 10 TeV and with the incidence angle from 0*∘* to 60*∘*.

Limited by the relatively low statistics of DAMPE gamma-ray data, the chi square method is not suitable for the data analysis, maximum likelihood method (Mattox et al. 1996) is adopted. Combining the IRFs and the model of gamma-ray sources, we can calculate the expected photon number recorded by the detector. The values, and also the uncertainties, of the parameters in the gamma-ray source model then can be estimated by comparing with the real DAMPE observation using the maximum likelihood method.

The data preparation, the convolution with the IRFs and the parameter inference are realized for DAMPE data analysis using a dedicated software named DmpST, which is also developed to facilitate the scientific analysis. In this paper, we introduce both the DAMPE IRFs and the DmpST software.

This paper is structured as the following. We first introduce the IRFs of DAMPE in Section 2. The observing time and exposure of DAMPE are then described in Section 3. In Section 4, we introduce the maximum likelihood method for DAMPE gamma-ray data analysis, followed by a description of the code structures in Section 5. We summarize this work in Section 6.

2 Instrument Response Functions

Instrument response functions (IRFs) are the parameterized representations of the instrument performance. The DAMPE IRFs can be factorized into three parts (Ackermann et al. 2012). The effective area, $A_{\mathrm{eff}}(E,\hat{v},s)$, is the product of the geometrical cross-section area, the probability of gamma-ray conversion and the efficiency of photon selection for a gamma-ray with energy $E$ and direction $\hat{v}$ in the detector reference frame. The $s$ denotes the trigger type (see below). The Point-spread function (PSF), $P(\hat{v}^{\prime};E,\hat{v},s)$ and the energy dispersion function, $D(E^{\prime};E,\hat{v},s)$ are the probability distributions of the reconstructed direction $\hat{v}^{\prime}$ and the reconstructed energy $E^{\prime}$ for a gamma-ray with energy $E$ and direction $\hat{v}$.

Given the spatial and spectral model of the incidence gamma-ray sources, $F(E,\hat{p})$, where $\hat{p}$ refers to the celestial directions of the gamma-ray sources, we can convolve the model with the IRFs to predict the distribution of observed photons:

[TABLE]

where $\hat{p}^{\prime}$ is the reconstructed celestial directions of the gamma-rays. The integrals are over the time and energy range of interest and the solid angle in the celestial reference frame.

To evalute the DAMPE IRFs, we perform Geant4-based Monte Carlo detector simulation to generate pseudo-photons of DAMPE (MC data hereafter). We simulate gamma-rays with uniform distribution of incidence direction, that can be used to explore the instrument response across the entire field of view (FoV) of DAMPE. The MC data are generated with an $E^{-1}$ counts spectrum uniformly in the logarithm energy, and from a sphere with 6 m2 cross-sectional area centered on the detector to cover the whole energy range and the whole detector of DAMPE. The directions of the gamma rays are sampled uniformly in solid angle with downward-going directions, leading to a semi-isotropic incidence flux of the simulated gamma rays. Here we ignore the back-entering events, because these events would have to traverse a large amount of material and thus presumably lose a lot of energy along their way. Through the same reconstruction and gamma-ray selection algorithm as the on-orbit data, the MC data can describe the response of DAMPE for gamma-ray observation accurately (Xu et al. 2018).

DAMPE uses two sets of trigger directives for physics data: the pre-scaled Low Energy Trigger (LET) and the Higt Energy Trigger (HET). The pre-scale factors of LET are different when the detector is in different geographic latitude (Chang et al. 2017). When the detector is in the low latitude region ($|\phi_{\rm g}|<20^{\circ}$), the Lower Energy trigger is pre-scaled with a factor of 8; and at high latitude region ($|\phi_{\rm g}|>20^{\circ}$) it is 64 pre-scaled. The IRFs is also divided into two sub-sets, LET IRFs and HET IRFs.

2.1 Effective area

Effective area is a numerical function varying with the energy of gamma-ray photon and its incidence direction in the instrument reference frame. We binned the MC data according to the event energy, incidence angle and trigger type. The effective area for each bin centered at $E_{i}$, $\theta_{j}$, $\phi_{k}$ with trigger type $s$ is

[TABLE]

where $N_{\mathrm{sim},i,j,k}$ is the number of photons generated in the simulation in each bin, and $N_{i,j,k,s}$ is the number of photons that passing the selection algorithm with trigger type $s={\rm LET}$ or ${\rm HET}$. The $A_{\mathrm{sim}}$ is the cross-section area of the generated sphere in the simulation.

We divide the MC data into 20 energy bins from 1GeV to 100GeV (40 energy bins from 1 GeV to 10 TeV) and 10 angular bins from 0*∘* to 60*∘* for LET (HET) data. Fig. 1 shows the effective area of DAMPE gamma-ray observation as a function of the energy and incidence direction.

2.2 Point-spread function (PSF)

The reconstructed direction ($\hat{v}^{\prime}$) of the photon may deviate from its true value ($\hat{v}$), and the probability distribution of the deviation $\delta v=|\hat{v}^{\prime}-\hat{v}|$ is parameterized by the PSF. The PSF for the DAMPE is related to the inclination angle $\theta$ and the azimuth angle $\phi$ of the incidence photon in the detector reference frame, and also the photon’s energy and trigger type. Because the $\phi$ dependence of the PSF is much weaker than the $\theta$ dependence, we ignore the $\phi$ dependence in the current version of the PSF.

Based on the MC data, we construct a histogram of the angular deviations of the selected gamma-rays for each energy and incidence angle bin and for each trigger type. We find that the form of the Fermi-LAT PSF (Ackermann et al. 2012) can accommodate DAMPE simulation data well. Accordingly, the PSF histogram is fitted with a double King function,

[TABLE]

where $K(x_{p};\sigma,\gamma)$ is King function defined as

[TABLE]

and $x_{p}$ is the scaled angular deviation

[TABLE]

The $S_{\mathrm{p}}(E,\theta)$ is the angular resolution (defined as 68$\%$ containment of the angular deviation) at energy $E$ and incidence angle $\theta$. The functional form of the King profile originates from XMM Newton (Kirsch et al. 2004, Read et al. 2011) and was later adapted for the Fermi-LAT. Note that the King function is normalized, i.e., $\int_{0}^{\infty}2\pi xK(x;\sigma,\gamma)\mathrm{d}x=1$.

We divide the MC data into 4 energy bins from 1GeV to 100GeV (8 energy bins from 1 GeV to 10 TeV) and 5 angular bins from 0*∘* to 60*∘* for LET (HET) data. Fig. 2 shows the angular resolution of DAMPE for gamma-ray observation at different energy and incidence direction. For each bin, the MC data are fitted with above functions and the best-fit parameters of them are derived and stored in the DmpST. Fig. 3 shows an example of the best fit to the scaled angular deviation with the double King function in the bin of $E\in[3.16,10]\,{\rm GeV}$ and $\theta\in[25.84^{\circ},36.87^{\circ}]$ for HET photons.

2.3 Energy dispersion

Energy dispersion function gives the probability of a photon with true energy ($E$) being allocated an energy ($E^{\prime}$) after the events reconstruction. Similar to the PSF, we ignore the $\phi$ dependence and parameterize the energy dispersion as function of scaled energy deviation

[TABLE]

where the scale $S_{\mathrm{D}}(E,\theta)$ is the energy resolution (defined as the half-width of the 68$\%$ containment range of the energy deviation) at the bin center of energy $E$ and incidence angle $\theta$. We fit the MC data with three piecewise functions of the form

[TABLE]

[TABLE]

We divided the MC data with the same binned method with the PSF. Fig. 4 shows the energy resolution of DAMPE for gamma-ray observation at different energy and incidence direction. And fit the energy dispersion with above function in each bin. Fig. 5 shows an example of energy dispersion fitted with the function in the bin of $E\in[3.16,10]\,{\rm GeV}$ and $\theta\in[25.84^{\circ},36.87^{\circ}]$ for HET photons.

3 Observing time and Exposure

For a particular source in the sky, its direction in the detector reference frame varies with the time. Since the IRFs various appreciably across the DAMPE field of view (FoV), we define the exposure $\epsilon$ for any given energy $E$ and direction in the sky $\hat{p}$ as the integral of the effective area over the time range of interest,

[TABLE]

The exposure can also be expressed as an integral over the solid angle in the detector reference frame,

[TABLE]

here the $t_{\mathrm{obs}}(\hat{v};\hat{p})$ is named observing time and defined as the total time in the range of interest during which DAMPE have observed the direction $\hat{p}$ with detector frame direction $\hat{v}$. The $A_{\mathrm{eff}}^{\mathrm{LET}}$ and $A_{\mathrm{eff}}^{\mathrm{HET}}$ in Eq. (9) are the effective area for LET and HET photons, respectively. As an example, we show the observing time map in the detector reference frame for the Vela pulsar in Fig. 6. With the observing time map and the DAMPE effective area, the exposure then can be calculated according to Eq. (9). Fig. 7 shows the all-sky exposure map of DAMPE at 10 GeV for the first year of operation. Because DAMPE is in a sun-synchronous orbit, we can see that the exposure is not uniform over the sky.

4 Maximum Likelihood Analysis

Analyzing the gamma-ray data from DAMPE requires the maximum likelihood method due to the limited number of photons and the angular resolution. We characterize a source by its photon flux density $F(E,\hat{p},t;\vec{\lambda})$. In order to reduce the computational burden, we assume the source is stationary during the time range in each likelihood analysis. 444For the variable source, the time dependence of the flux can be achieved by repeating the analysis in finer time bins. The model of gamma-ray source then can be modeled by

[TABLE]

Here the $M(\hat{p})$ is a normalized function describing the spatial morphology of the source. For the point source, the spatial distribution can be described with the Dirac delta function, $M(\hat{p})=\delta(\hat{p}-\hat{p}_{0})$, where $\hat{p}_{0}$ is the direction of the point source. The $S(E;\vec{\lambda})$ in Equation (10) is the spectrum of the source with its parameters $\vec{\lambda}$.

To remove the $\theta$ dependence of the PSF and energy dispersion, we calculate the exposure-weighted PSF and energy dispersion for any sources included in the analysis:

[TABLE]

[TABLE]

Considering the excellent energy resolution of DAMPE (i.e., $\sim 5\,\%$ at 1 GeV and $\sim 1\,\%$ at 100 GeV (Chang et al. 2017)), the influence of energy dispersion can be ignored for most gamma-ray science analysis. The only exception is the case of searching for narrow-line feature in the gamma-ray spectrum (Ackermann et al. 2015; Liang et al. 2016; Li et al. 2018), which will be performed with other dedicated code. So in the DmpST we ignore the energy dispersion and regard the measured energy as the true photon energy in current version, and it will be considered in the future if the statistic allows.

With the parameterized source model, the exposure and the exposure-weighted PSF, we can calculate the model predicted photon rate in the bin $i$ (centered on $E_{i},\hat{p}^{\prime}_{i}$) from the source $j$:

[TABLE]

The predicted photon rates are compared to the observation data to determine the model parameters. The information we can get from the DAMPE observation is the energy ($E$), the direction ($\hat{p}^{\prime}$) and the time of arrival ($t$) of each photon. We bin the photons in the region-of-interest (ROI) into a counts cube according to their measured energies and directions. For each bin, the photon number $N$ follows the Poisson distribution with unknown mean $R$: $p(N;R)={R^{N}}/{N!}\cdot\exp(-R)$. Taking into account all the bins with numbers $\{N_{i}\}$, the Poisson distribution becomes

[TABLE]

Because of the broad PSF of DAMPE and the strong Galactic diffuse background, the photons in each bin may originate from multiple sources, the parameters of which should be determined simultaneously utilizing the likelihood fitting. With the model predicted photon rates and the real observed data, and based on the Poisson statistics, we construct the binned likelihood function (in logarithm form) by summing over all $N_{\mathrm{bins}}$ bins and all $N_{\mathrm{s}}$ sources:

[TABLE]

where the $R_{ij}$ is the model expected photon number in the bin $i$ from source $j$ and the integral is calculated in the corresponding bin $i$ as well.

When the bin widths are taken to be infinitesimal such that only 0 or 1 photon in each bin, the summation over $N_{\mathrm{bins}}$ bins becomes to an integral over the whole energy range and the ROI. Then we get the unbinned form of the likelihood function:

[TABLE]

By maximizing the likelihood function of (15) or (16), we can get the best-fit values of all the free parameters in the source models.

5 Implementation

The code is coded with Python, based on NumPy (van der Walt et al. 2011), SciPy 555http://www.scipy.org, AstroPy (Robitaille et al. 2013) and iminuit (James & Roos 1984) packages. The structure of DmpST is shown in Fig. 8.

The input modules are Events, SpaceCraft, DmpIRFs, Spatial Model, Spectrum and Model (shown as blue in Fig. 8). The Events module stores the information of photons that are selected from all the events detected by DAMPE using the photon selection algorithm (Xu et al. 2018). The information of a photon includes the arrive time ($t$), the reconstructed energy ($E$), the reconstructed direction in the celestial coordinates ($\alpha_{2000}$, $\delta_{2000}$, $l$, $b$) and in the detector reference frame ($\theta$, $\phi$), and the trigger type (s). The photons of interest in the analysis can be selected according to their times, energies or directions utilizing the Events module and can be binned into a counts map or a counts cube which is managed by the Sky Map module. The SpaceCraft module stores the position, direction and livetime of DAMPE along with time, and can be used to calculate the observing time of DAMPE for any direction in the sky (see Section 3). The DmpIRFs module is used to manage information of instrument response functions (IRFs), including the effective area matrix, the parameters of PSF and energy dispersion function. With these parameters and the fitting functions described in Section 2, the distributions of PSF and energy dispersion can be reconstructed. The Spatial Model and Spectrum modules provide different kinds of spatial and spectral models of gamma-ray sources, respectively. The Model module includes all the models of sources those will contribute photons to the ROI.

The process modules comprise Sky Map, Exposure, Source, Likelihood Base, Binned Likelihood and Unbinned Likelihood (white parts in Figure 8). The Sky Map module manages the information of counts map or counts cube from the Events module, such as the photon number and celestial coordinates of each bin. The Exposure module calculates the observing time, exposure and the exposure-weighted PSF and energy dispersion based on the information in SpaceCraft and DmpIRFs modules. The Source module combines spatial and spectral models based on the Spatial Model and Spactrum modules for each source in the Model module. The Likelihood Base module convolves the PSF with the spatial model, integrals the spectrum over the energy to calculate the expected photons number for each source based on the Sky Map, Exposure and Source modules. The Binned/Unbinned Likelihood modules construct the likelihood function described in Section 4.

Finally, the Likelihood Analysis module implements the maximum likelihood estimation with the Minuit algorithm and the basic outputs are the best-fit values ($\hat{\vec{\lambda}}$) of source parameters, the source fluxes and the corresponding statistic uncertainties. Also we can obtain the confidence level of each source defined as

[TABLE]

where $\hat{\vec{\lambda}}_{0,j}$ is the best-fit parameters without source $j$ included in the model. The $TS_{j}$ follows $\chi^{2}$ distribution with $h-m$ degrees of freedom (Wilks 1938), where $h$ and $m$ are the number of free parameters in the model with/out source $j$. The DmpSkyObsSimu module simulates photons observed by DAMPE with the DmpIRFs, SpaceCraft and Source modules.

Monte Carlo simulation has been done with the DmpSkyObsSimu module with the Galactic diffuse emission and isotropic emission. With the Likelihood Analysis module, we analyze the simulated data to confirm the distribution of the $TS$. The null hypothesis is there is no point source, only the background including the Galactic diffuse emission and isotropic emission. The alternative hypothesis is the converse: there is a point source with Power-Law spectrum with free normalization parameter. For most point source analysis of DAMPE, the radii of ROI is $\approx 2\times S_{\mathrm{p}}$ and the typical number of photons $N$ in the ROI is about 25. Fig 9 shows that for $TS>0$, the distribution of $TS$ is following $\chi^{2}_{1}/2$, and the one-half of the simulations have TS = 0 (Mattox et al. 1996).

6 Summary

The GeV gamma-ray sky is an important observation target of DAMPE. To facilitate analyzing the DAMPE gamma-ray data, we have developed a dedicated software named DmpST, which implements maximum likelihood analysis to extract the parameters of sources that attribute to the observed gamma-rays. The DAMPE IRFs that are essential to the gamma-ray data analysis, including the effective area, the PSF and the energy dispersion, are also derived based on high-statistics simulation data. Making use of the DmpIRFs and DmpST that are detailed in this paper, scientific analyses of the gamma-ray data could be carried out to obtain the best-fit spectral parameters, fluxes and corresponding statistic uncertainties, and further the spectral energy distribution and light curve of the gamma-ray sources, promoting our understanding the nature of high energy gamma-ray phenomena.

Acknowledgements.

This work is supported in part by National Key Program for Research and Development (No. 2016YFA0400200), the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDB23040000), the 13th Five-year Informatization Plan of Chinese Academy of Sciences (No. XXH13506), the National Natural Science Foundation of China (Nos. U1631111, U1738123, U1738136, U1738210), Youth Innovation Promotion Association of Chinese Academy of Sciences, and the Young Elite Scientists Sponsorship Program. In Europe the activities and the data analysis are supported by the Swiss National Science Foundation (SNSF), Switzerland; the National Institute for Nuclear Physics (INFN), Italy.