GammaLib

Computing the model cube with a sampling of 0.001 deg per pixel does not change the situation, hence it’s not a question of spatial sampling, as shown by the figure below.

The problem comes from the discret integration steps in GCTAResponseIrf::irf_diffuse, with the number of Romberg iterations in rho and phi fixed to 5, which corresponds to 16 integration steps in radius and azimuth. Below the resulting model map in the 100-140 GeV bin for rho=phi=5 (left, default), rho=6, phi=5 (centre) and rho=5, phi=6. Increasing the number of steps in the radial direction smoothes the result radially, but still leads to a star-like structure. Increasing the number of steps in the azimuthal direction doubles the number of azimuthal steps and hence reduces the star-like appearance of the model.

The problem arises in fact for sky maps that have a spatial resolution that is smaller than the CTA PSF. Due to the energy dependence this happens actually for many sky maps.

The left image below illustrated the problem. It shows a high-resolution model map for iter_rho=iter_phi=5 generate at 20 GeV, i.e. at an energy where the PSF is worst. At each point that was sampled, the image shows a little replica of the input sky map (a whirly structure around a spot). Increasing the number of iteration steps produces the image on the right. There is still some structure due to the discrete integration steps, but overall the model is now rather smooth.

To understand which integration precision is required, I did an event simulation and fitted the sky map using an unbinned analysis to the simulated events, varying the number of integration steps. An energy range of 20 GeV - 100 TeV was used. The results are summarized in the table below. For iter_rho=iter_phi=8 on the fit does no longer improve.

I then implemented an algorithm where the actual integration precision depends on the spatial resolution of the map, computing iter_rho and iter_phi as follows:

iter_rho = log_2(rho_max/resolution) + 1
iter_phi = log_2(2*pi*rho/resolution) + 1

I then validated the algorithm using the simulation of Michele. The table below shows the results. For reference, the fit with the old scheme is shown in the first row. The next three rows show the results of the new scheme with a maximum number of iterations set to 7, 8 and 9, respectively. Also here, the TS value and fit parameters converge for a value of 8, hence we retain this value for the code. Note that with the new scheme, the TS value for the radio map are now larger than the values for the point source, radial disk and Gaussian disk models.

Merged code into devel.