AGILESim: Monte Carlo Simulation of the AGILE Gamma-Ray Telescope

Because our goal is to study the reconstruction of pair events, we ignore here the presence of the Super-AGILE instrument, assuming a negligible attenuation efficiency toward gamma-rays. It should be included for simulations of the particle-induced background level, however.

The AGILESim three-dimensional mass model is shown in Figure 1 in all its components: the silicon tracker and the MCAL (top row), the ACS panels, the electronic boards, and the spacecraft equivalent mass (bottom row).

Zoom In Zoom Out Reset image size

3.1.1. Silicon Tracker

3.1.2. MCAL

The MCAL scintillation calorimeter consists of 30 Csi(Tl) bars, 15 × 23 × 375 mm³ each, divided into two orthogonal sets of 15 rows, for a total radiation length of 1.5 X₀ (Figure 1, top right panel). The light produced in the bars is converted into charge by two PIN photodiodes, one at each side of the bar. The longer the path crossed by the light, the higher the self-attenuation induced by the scintillation material. After computing the hit position along the bar, the energy readout by a single diode is emulated by applying experimental attenuation coefficients μ, one for each side, to half the energy deposit in the form

$$\begin{eqnarray}&&{E}_{A}=0.5\times E\times {e}^{-\mu t},\end{eqnarray} \tag{ 1 }$$

where E is the total deposited energy and t is the distance of the energy deposit from the diode.

3.1.3. ACS, Electronic Boards, and Spacecraft

From the Monte Carlo simulation we directly obtain the energy deposit in each ST strip, which in reality is the product of a complex readout process depending on the design of the electronic devices. The readout process can be analog (e.g., AGILE), where the true charge is collected from the strips, or digital (e.g., Fermi), where a flag is assigned to the triggered strips. The analog readout increases the spatial resolution, hence the ability to track the electron–positron pair and consequently the angular resolution.

Including the GRID readout logic in the simulation requires applying the charge-sharing process between adjacent strips during the postprocessing of the BoGEMMS simulation output. Charge-sharing is caused by diffusion during the charge collection and capacitive coupling between the strips. The first effect is due to diffusion during the drift of the electron–hole pairs along the field lines that increases the size of the charge cloud. The capacitive coupling is obtained by inserting additional strips between the readout strips. In the ST, the readout pitch, i.e., the distance between two electronic channels connected to the strips, is 242 μm, while the strip pitch, i.e., the physical distance between strips, is 121 μm. The intermediate strips, or floating strips, are placed between contiguous readout strips. If charge is collected by a floating strip, the capacitive coupling induces image charges in the adjacent readout strips, so that the signal can be measured after applying an interpolation algorithm to the charge distribution. The use of floating strips increases the spatial resolution while reducing the number of electronic channels, hence the power consumption. The effect of charge-sharing in the simulation is obtained by assigning a fraction of the energy deposited in one strip to the neighboring ones (see Cattaneo 1990, 2010 for details). We built for this purpose a tracker look-up table containing only readout strips and saving in each one an energy amount given by the relation published in Longo et al. (2002) and obtained from test beam measurements (Barbiellini et al. 2002). Defining k and i the readout and global strip indexes, and their relation k = 1 + (i − 1)/2, the total energy deposited in the kth readout strip is E_k = E_i + 0.38·(E_i−1 + E_i+1) + 0.115·(E_i−2 + E_i+2) + 0.095·(E_i−3 + E_i+3) + 0.045·(E_i−4 + E_i+4) + 0.035·(E_i−5 + E_i+5).

The instrumental noise rms of the GRID strips is of the order of 30 ADC counts (Bulgarelli et al. 2010). Given a conversion factor of 0.174 keV ADC⁻¹ counts, a σ ∼ 5 keV is added to each strip, although it has negligible effects on the trigger rates (Longo et al. 2002). After summing all energy deposits within the same strip, generated by the same primary simulated photons, we apply a trigger threshold of one-fourth of the energy of a minimum ionizing particle (MIP, ∼27 keV) to the ST strips (Barbiellini et al. 2002), as in the real case.

In the simulation of AGILE observations, the correct reconstruction of the conversion process of gamma-rays to electron–positron pairs is not enough. Being an indirect detection process, obtaining the primary gamma-ray energy and direction requires a trigger and filtering algorithm to be applied to the Monte Carlo simulation. We make use of the AGILE/GRID standard processing pipeline (Bulgarelli 2019), which takes as input the trigger list from the instruments (ST, MCAL, and AC) and at the end produces a classified event list for the scientific analysis (e.g., production of count maps, spectra, and light curves of the science target).

The first block of the pipeline is the onboard Payload Data Handling Unit (PDHU, Argan et al. 2008, 2019) that applies the hardware and software veto logic to the scientific data and produces the level 2 telemetry data format to be transmitted to the ground station. When Monte Carlo simulations are used as input, the Data Handling Simulator (DHSim) first converts the Monte Carlo data into the instrument trigger data format and then acts as a virtual Data Handling Unit, producing the telemetry files for the on-ground filtering pipeline. In order to use the DHSim software, we converted the AGILESim reference system and data output into the same format as is used in GAMS, including the X and Y views of the ST with the triggered strips in MIP units, the fired AC panels, and the energy deposited in the MCAL bars.

The on-ground FM3.119 filter is then used for the track identification and reconstruction, which first selects the potential tracks of the pair (Pittori & Tavani 2002) and then fits the track parameters using a recursive Kalman filter algorithm (Giuliani et al. 2006). The track minimizing the sum of the χ² of the ST planes is associated with the primary gamma-ray conversion.

We developed an independent GRID event visualizer to directly display the BoGEMMS tracks and compare them to the tracks loaded by the on-ground filter. In Figure 2 (top panel) the AGILESim GRID visualizer displays as example a 100 MeV gamma-ray with an incoming direction parallel to the telescope axis. The same event is displayed in Figure 2 (bottom panels) by the DHSim and the AGILE filter quick-look software (Bulgarelli et al. 2008). Figure 3 shows a simplified schema of the processing pipeline from the onboard or simulated data readout to the on-ground filter.

Zoom In Zoom Out Reset image size

The ST conversion efficiency, i.e., the probability of converting a gamma-ray into an electron–positron pair, depends on the tracker design, in addition to the input photon energy and direction. We evaluate this efficiency in AGILESim by simulating mono-energetic parallel photon beams with energy 50, 100, 200, 400, and 1000 MeV and with an off-axis angle θ = 30°, then computing the ratio between the gamma-rays, producing a pair and those that do not interact (Compton scatterings can be neglected at this regime), as a function of the number of crossed planes. The pair-production attenuation coefficients μ of the material composing the ST trays are tabulated by the National Institute of Standards and Technology (Berger & Hubbell 1987). They are used in the Beer–Lambert formula to compute the expected conversion efficiency for each layer of thickness t, in the form $(1-{e}^{-\mu t/\cos (\theta )})$. The resulting simulated conversion efficiency is shown in blue in Figure 4 for an input energy of 100 MeV and compared with the corresponding analytical values, plotted in red. The error bars are Poisson statistical fluctuations from the number of simulated photons. The simulation reproduces the theoretical efficiency based on the material cross-section quite well. This comparison verifies not only the correct building of the AGILESim ST mass model, but also the Geant4 pair-production physics models implemented in the BoGEMMS simulation framework.

Zoom In Zoom Out Reset image size

Testing the correct simulation of electrons and positrons interacting with the Si strip detectors is another fundamental step of the physics verification. In the AGILE Assembly, Integration, and Verification (AIV) on-ground activity (Bulgarelli et al. 2010), the interaction of atmospheric muons was used to characterize the ST noise and frontend readout. Atmospheric muons behave as MIPs, i.e., they deposit in the Si layers a minimum amount of energy per path unit that is almost independent of the incident energy. The electrons and positrons interacting with the ST strips behave as MIPs as well, and we can use the cluster charge distribution collected in the AIV stage to verify the AGILESim simulation.

A cluster of strips is a set of contiguous readout strips triggered by the passage of a particle, and its identification is performed on board by the PDHU. The ST pull is defined as the ratio between the maximum signal in the cluster strips and the instrumental noise rms, having a mean value of 5 keV (see Section 3.2). The number of strips populating the clusters, defined as cluster width, and their pull were measured during the GRID on-ground noise tests for a set of inclination angles of atmospheric muons with respect to the GRID axis. In order to compare such distributions with the AGILESim output, we simulated the on-ground atmospheric muons as beams of parallel particles hitting on the tracker surface at two inclination values (θ = 1° and 30°) and an averaged azimuthal angle of 45° (225° in the BoGEMMS system of reference). According to Shukla & Sankrith (2018), the on-ground spectral distribution I(E) of muons with momentum in the range from 1 GeV/c to 1 TeV/c can be modeled as

$$\begin{eqnarray}&&I(E)={I}_{0}\,N\,{\left({E}_{0}+E\right)}^{-n}\ {\left(1+\displaystyle \frac{E}{\epsilon }\right)}^{-1},\end{eqnarray} \tag{ 2 }$$

where E is the kinetic energy of the muon expressed in GeV, I₀ = 72.5 m⁻² s⁻¹ sr⁻¹, n = 3.06, E₀ = 3.87 GeV, and  = 854 GeV at sea level.¹⁵

The resulting cluster width and pull distributions are shown in Figure 5 (top and bottom panels, respectively) and compared to the real on-ground measurement (Bulgarelli et al. 2010) at similar muon incident angles.

Zoom In Zoom Out Reset image size

For large incidence angles, the probability for a muon of crossing at the same time both reading and floating strips is quite high. For this reason, at θ ∼ 30° most of the clusters are composed by two strips, whereas at inclination near to zero, the probability of having single-strip or double-strip clusters is almost the same. For both inclinations, the simulated cluster width distributions reproduce the real ones very well.

When a muon hits a floating strip, the signal is given by the charge induced on the adjacent readout strips, and its amplitude is lower because only a fraction of the generated charge can contribute to the signal. In contrast, when a muon hits a readout strip, the signal amplitude is given by the full amount of the produced charge. These two mechanisms determine the shape of the pull distribution. Indeed, for a small incidence angle, when floating and readout strips are most likely individually hit, the pull distribution shows two separated peaks, where the one at lower pull values is due to the floating strip contribution. Conversely, when the incidence angle is higher, each muon can affect both the floating and readout strips, and the two peaks merge into a single peak. In both cases, the simulated response reproduces the real data quite well.

The angular resolution of a gamma-ray telescope defines its ability to distinguish between two point sources placed in proximity to each other. In addition to avoiding source confusion in crowded sky regions (e.g., the galactic center), the angular resolution is a key factor affecting the instrument performance because the smaller the imaged shape of the point source, the lower the background contribution within, hence the better the overall sensitivity for detecting faint sources. The point-spread function (PSF) is a probability distribution for photons from a point source that is reconstructed at a certain spherical distance from their true direction. Assuming azimuthal symmetry of such function, as is observed in gamma-ray telescopes below an off-axis angle of 30° (Sabatini et al. 2015), we describe the PSF in terms of surface brightness, i.e., counts per solid angle, as a function of the spherical distance from the true direction. The output of the AGILE filter and reconstruction pipeline for a certain number of simulated gamma-rays at (θ, ϕ) = (30°, 45°) and an input true energy E gives a list of events with reconstructed direction θ_rec and ϕ_rec and reconstructed energy E_rec. A G flag labels all the events that are reconstructed as a likely gamma-ray (Chen et al. 2013). All results here refer to G-flagged events, as those used in the in-flight scientific analysis. We simulate parallel beams of photons hitting the full surface of the tracker, each with a constant energy of 50, 100, 200, 400, and 1000 MeV. For our analysis we select all the photons with a reconstructed energy within defined energy bands (30–70, 70–140, 140–300, 300–700, and 700–1700 MeV) to compare BoGEMMS with the GAMS simulations reported in Sabatini et al. (2015).

The spherical distance d_sph is given by

$$\begin{eqnarray}&&{d}_{\mathrm{sph}}={\cos }^{-1}(\sin ({\theta }_{\mathrm{rec}})\sin (\theta )\cos ({\phi }_{\mathrm{rec}}-\phi )+\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&+\cos ({\theta }_{\mathrm{rec}})\cos (\theta )).\end{eqnarray} \tag{ 4 }$$

The photon true direction (θ, ϕ) is assumed to be the center of a new system of reference, and each reconstructed photon is characterized by its radial distance from the center. The spherical distance from the true direction is then equally binned in annular regions around the center, subtending an increasing solid angle from the center. The number of counts in each region is divided by its subtended solid angle to compute the surface brightness, and then the distribution maximum is normalized to 1. The uncertainty for each bin is given by the Poisson statistics. The simulated PSF radial profiles at 50, 100, 200, 400, and 1000 MeV are shown in Figure 6. In gamma-ray astronomy, the angular resolution is described as the radius containing a certain percentage of photons in the PSF, usually the 68%, or the full width at half-maximum (FWHM) of its radial profile. The containment radius (CR68 from now on) takes into account the whole distribution, including the length of the tail. The FWHM instead measures the central shape of the profile. For these reasons, the PSF CR68 and FWHM can have different values and a different dependence on the source energy and off-axis angle.

Zoom In Zoom Out Reset image size

Both parameters can be computed analytically, by interpolation, or by fitting the PSF with an appropriate function. We obtain the analytical 68% containment radius by just counting the reconstructed photons: in this case, the 1σ error is given by the statistics of the simulation in terms of Poisson uncertainty $\sqrt{N}$. The normalized radial profile can be linearly interpolated by a function and the spherical distance at which the function is 0.5 defines the half-width at half-maximum. Twice this value gives the FWHM. For each bin of the radial profile, the 1σ lower and the upper limits are also interpolated in the same way to extract the corresponding FWHM upper and lower limits. We also report the CR68 obtained from the interpolation function, with errors computed from the 1σ lower and upper limits as the FWHM.

The PSF radial profile of a point source is often best described not by a Gaussian profile, but by a superimposition of an exponential core and an r² tail, with r indicating the spherical or radial distance. This function, presented in King (1971) and referred to as King function in Read et al. (2011) for the evaluation of XMM-Newton angular resolution, takes the form

$$\begin{eqnarray}&&K(r,\sigma ,\gamma )=\displaystyle \frac{1}{2\pi {\sigma }^{2}}\left(1-\displaystyle \frac{1}{\gamma }\right){\left(1+\displaystyle \frac{{r}^{2}}{2{\sigma }^{2}\gamma }\right)}^{-\gamma },\end{eqnarray} \tag{ 5 }$$

with σ and γ characterizing the size of the angular distribution and the weight of the tail, respectively. In the small-angle approximation $\sin (r)\simeq r$, the King profile integrated over the radial distance is normalized to 1, as a probability distribution. Because we normalize the PSF to its maximum, a multiplicative normalization parameter is added to the King function. The AGILE/GRID calibration (Chen et al. 2013) used a single King function to model the PSF radial profile from Monte Carlo simulations and produce the mission instrument response function. We decided to follow the same procedures as the mission calibration team, and optimization studies of the AGILE scientific performance will be the subject of future works on this subject.

The simulated PSF radial profile is fit by the King profile, and the resulting σ and γ parameters are used for the evaluation of the CR68 radius as derived by Cattaneo et al. (2018),

$$\begin{eqnarray}&&\mathrm{CR}68=\sigma \sqrt{2\gamma \left({\left(1-0.68\right)}^{\tfrac{1}{1-\gamma }}-1\right)}.\end{eqnarray} \tag{ 6 }$$

The 1σ error on the CR68 is obtained from the 68% confidence contour levels of the σ and γ best-fit values. The best-fit King models are superimposed on the radial profiles in Figure 6, and the residuals, in standard deviations, are plotted in the bottom panels. We are able to describe all mono-energetic profiles with a single King profile with residuals below a 3σ deviation from the model.

The resulting AGILESim PSF obtained for the two methods, CR68 and FWHM, and three different procedures (analytical, interpolation, and fitting) are listed in Table 1, which also reports for comparison the angular resolution obtained by the GAMS simulator in the same energy channels (Sabatini et al. 2015). The CR68 and FWHM values from GAMS are obtained from analytical procedures, and the reported errors are derived from the binning applied to the radial profile. The last panel in Figure 6 summarizes the tabulated angular resolutions obtained with AGILESim and GAMS simulations and their evolution with the true input energy.

Table 1.  Comparison of the PSF Values for AGILE/GRID Obtained by Simulating a Mono-energetic Source with AGILEsim and GAMS (Sabatini et al. 2015)

PSF—Mono-energetic Source (θ = 30°, ϕ = 45°)
	AGILESim (G-class)					GAMS (G-class)
Energy	Analytical	Interpolation		King Fit		Analytical
(MeV)	CR68	CR68	FWHM	CR68	χ2/dof	CR68	FWHM
	(°)	(°)	(°)	(°)		(°)	(°)
50	8.0 ± 0.1	${7.7}_{-0.1}^{+0.1}$	${8.8}_{-0.2}^{+0.2}$	${8.0}_{-0.2}^{+0.2}$	38/27	8.5 ± 0.5	9.0 ± 1.0
(30–70)
100	4.3 ± 0.03	${4.1}_{-0.1}^{+0.1}$	${4.2}_{-0.1}^{+0.1}$	${4.8}_{-0.2}^{+0.2}$	94/23	5.0 ± 0.25	4.5 ± 0.5
(70–140)
200	2.6 ± 0.02	${2.5}_{-0.03}^{+0.04}$	${2.2}_{-0.03}^{+0.03}$	${2.6}_{-0.1}^{+0.1}$	16/13	2.7 ± 0.1	2.4 ± 0.2
(140–300)
400	1.2 ± 0.01	${1.1}_{-0.03}^{+0.02}$	${1.1}_{-0.02}^{+0.02}$	${1.1}_{-0.01}^{+0.01}$	27/27	1.4 ± 0.1	1.4 ± 0.2
(300–700)
1000	0.58 ± 0.004	${0.51}_{-0.004}^{+0.004}$	${0.65}_{-0.01}^{+0.01}$	${0.56}_{-0.01}^{+0.01}$	19/18	0.7 ± 0.1	0.8 ± 0.2
(700–1700)

Download table as:  ASCIITypeset image

AGILESim features a slightly better angular resolution than the GAMS values, but they are always consistent within a 2σ uncertainty level. We established in general a good agreement between the two simulators, especially if the FWHM values are compared.

We achieve the scientific validation of the AGILESim framework by comparing the simulated PSF profile with the angular resolution obtained from AGILE/GRID observations of the Vela pulsar.

4.4.1. Vela Data Analysis

The Vela pulsar (l = 26359, b = −284) is the brightest nonflaring gamma-ray source in the sky and one of the first targets observed by AGILE (Pellizzoni et al. 2009) and Fermi (Abdo et al. 2010). With no other sources near it and with a steady spectral distribution, the Vela pulsar is an optimal point-like target for tracking gamma-ray telescopes such as AGILE and Fermi and represents a good candidate for the characterization of the GRID angular resolution.

From the second AGILE catalog of gamma-ray sources (Bulgarelli et al. 2019), the Vela pulsar (2AGL J0835−4514) phase-averaged gamma-ray emission is best modeled by a power law with a super-exponential cutoff of the form

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{E}^{-\alpha }\,\exp \left[-{\left(\displaystyle \frac{E}{{E}_{c}}\right)}^{\beta }\right],\end{eqnarray} \tag{ 7 }$$

with α = 1.71, β = 1.3, and E_c = 3.91 GeV. The validation is performed in the core of GRID operative energy range (100–1000 MeV), divided into two wide channels (100–300 and 300–1000 MeV).

We extracted all observations of Vela with θ in the 20°–40° interval to be consistent with simulations obtained at an inclination θ = 30°. The observation period covers almost two years, from 2007 December 1 (the end of the science verification phase) to 2009 October 15, the end of the AGILE operations in pointing mode. The same data set was used for the Vela analysis of the AGILE second catalog. To further reduce gamma-ray albedo contamination from the Earth, the data selection excludes photons coming within 80° from the reconstructed satellite-Earth vector.

Counts from the Vela region are extracted in aperture photometry by selecting from the G-flagged event list all photons with reconstructed energy included in the intervals 100–300 and 300–1000 MeV, and encircled in rings centered at the source, each with bin size of 01. For each bin, a Poisson statistical error is assumed. The surface brightness for each bin is obtained by dividing the counts in the ring by the subtended solid angle. Because Vela observations are not normalized in exposure, we do not take into account the effective area and the energy dispersion of AGILE/GRID to compare the observation with the output of the AGILESim processing pipeline.

The in-flight count maps in the 100–300 and 300–1000 MeV energy bands of the Vela pulsar are displayed in Figure 7. A bin size of 01 and a Gaussian smoothing are applied to both maps. The background surface brightness, in counts sr⁻¹, is extracted from an annular region at 8°–10° radial distance outside the AGILE/GRID PSF, and subtracted from the on-source surface brightness. The resulting background-subtracted radial profiles, with a binning of 04 (100–300 MeV) and 02 (300–1000 MeV), are shown in Figure 8. The error bars are statistical.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A single King function well describes both profiles (${\chi }^{2}/\mathrm{dof}\lt 1.4$), although, because of the large statistical errors in the bins and the reduced number of bins, we obtain large errors in the best-fit parameters. For this reason, CR68 is only evaluated by analytical and interpolation methods.

4.4.2. Comparison with Simulations

Because the Vela pulsar exponential cutoff is well above the energy upper limit of our analysis, the source is simulated in AGILESim as a simple power law with spectral index α = 1.71 for a θ = 30° inclination. We simulate the Vela pulsar with a wider energy range, from 30 MeV to 3 GeV, and select photon with a reconstructed energy in the 100–300 and 300–1000 MeV channels, to take into account a potential contamination from photons with a true energy outside the energy range of analysis. The simulated count maps with superimposed the contours of the in-flight PSF for the two energy ranges are shown in the left panels of Figure 9. The same binning and Gaussian smoothing as for the in-flight maps are applied. The radial profiles of simulated and in-flight Vela pulsar observations are compared in the right panels of Figure 9. The CR68 and FWHM are obtained by interpolation of the PSF profile, with the former also computed analytically from the list of background-free photons provided by the simulation. All the values are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 2.  Comparison of the AGILE/GRID PSF Values Obtained with the Simulation (First Three Columns) and from the In-flight Observation of the Vela Pulsar (Last Two Columns) in the 100–300 MeV and 300–1000 MeV Energy Ranges

PSF—Vela Pulsar ($\theta =30^\circ$, $\phi =45^\circ$)
	AGILESim			AGILE/GRID
Energy	Analytical	Interpolation		Interpolation
(MeV)	CR68	CR68	FWHM	CR68	FWHM
	(°)	(°)	(°)	(°)	(°)
100–300	3.1 ± 0.04	${3.0}_{-0.1}^{+0.1}$	${2.1}_{-0.1}^{+0.1}$	${3.7}_{-0.7}^{+0.4}$	${2.0}_{-0.3}^{+0.2}$
300–1000	2.7 ± 0.2	${2.6}_{-0.9}^{+1.0}$	${0.8}_{-0.1}^{+0.1}$	${2.7}_{-1.0}^{+0.9}$	${0.8}_{-0.1}^{+0.1}$

Download table as:  ASCIITypeset image

We find an in-flight angular resolution (FWHM) for Vela-like sources of ${2\buildrel{\circ}\over{.} 0}_{-0\buildrel{\circ}\over{.} 3}^{+0\buildrel{\circ}\over{.} 2}$ (100–300 MeV) and ${0\buildrel{\circ}\over{.} 8}_{-0\buildrel{\circ}\over{.} 1}^{+0\buildrel{\circ}\over{.} 1}$ (300–1000 MeV). Our values are consistent with those reported in Sabatini et al. (2015) for Crab-like sources, (2.5 ± 0.5)° (100–400 MeV) and (1.2 ± 0.5)° (400–1000 MeV), where a higher spectral index α = 2.1 decreases the angular resolution.

Because CR68 takes into account the full distribution, the presence of large tails because of the wide energy channels causes large uncertainties in CR68 obtained from interpolation, especially above 300 MeV, where the statistic is lower. The FWHM instead represents an unbiased measure of the angular resolution. We obtain a good agreement between the AGILESim and the observed angular resolution of AGILE/GRID in the shape of the radial profile and its measurables.