ctools

The issue with ctmodel was that a different computation of the model is needed for a stacked and an unbinned observation. For a stacked observation, the model value has to be weighted with the correct bin width, while for an unbinned observation this weighting needs not to be done since each observation will be evaluated using its proper IRF. I therefore implemented the following code change in ctmodel::fill_cube:

    // Flag if the CTA observation is a counts cube
    bool obs_is_cube = (obs->eventtype() == "CountsCube");
    ...
            // If observation is a counts cube then set the proper weight for
            // this bin
            if (obs_is_cube) {
                bin.weight(m_cube[ibin]->weight());
            }

Below the residual spectrum, residual profiles in Right Ascension and Declination, and a residual map for 20 energy bins. Overall the residuals look now reasonable.

The background model used in ctbkgmodel was determined by fitting all observations, including the source regions. Consequently, in the ctlike fit the normalisation of the background model is reduced, as can be see in the log file, shown below for the 20 energy bins fit. In this case the background model was reduced by about 7%, the index is compatible with zero. However the shape of the background model, and the relative weighting between the backgrounds of the different observations, is preserved in a stacked analysis, introducing some bias.

2018-11-07T10:30:09: === GCTAModelCubeBackground ===
2018-11-07T10:30:09:  Name ......................: BackgroundModel
2018-11-07T10:30:09:  Instruments ...............: CTA, HESS, MAGIC, VERITAS
2018-11-07T10:30:09:  Instrument scale factors ..: unity
2018-11-07T10:30:09:  Observation identifiers ...: all
2018-11-07T10:30:09:  Model type ................: "PowerLaw" * "Constant" 
2018-11-07T10:30:09:  Number of parameters ......: 4
2018-11-07T10:30:09:  Number of spectral par's ..: 3
2018-11-07T10:30:09:   Prefactor ................: 0.92688008398116 +/- 0.0138008183074852 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)
2018-11-07T10:30:09:   Index ....................: 0.0113778958948462 +/- 0.0144446074027589 [-5,5]  (free,scale=1,gradient)
2018-11-07T10:30:09:   PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
2018-11-07T10:30:09:  Number of temporal par's ..: 1
2018-11-07T10:30:09:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

2018-11-07T11:19:47: === GOptimizerLM ===
2018-11-07T11:19:47:  Optimized function value ..: 44648.462
2018-11-07T11:19:47:  Absolute precision ........: 0.005
2018-11-07T11:19:47:  Acceptable value decrease .: 2
2018-11-07T11:19:47:  Optimization status .......: converged
2018-11-07T11:19:47:  Number of parameters ......: 10
2018-11-07T11:19:47:  Number of free parameters .: 6
2018-11-07T11:19:47:  Number of iterations ......: 3
2018-11-07T11:19:47:  Lambda ....................: 1e-06
2018-11-07T11:19:47:  Maximum log likelihood ....: -44648.462
2018-11-07T11:19:47:  Observed events  (Nobs) ...: 8461.000
2018-11-07T11:19:47:  Predicted events (Npred) ..: 8455.000 (Nobs - Npred = 5.99964775683293)
2018-11-07T11:19:47: === GModels ===
2018-11-07T11:19:47:  Number of models ..........: 2
2018-11-07T11:19:47:  Number of parameters ......: 10
2018-11-07T11:19:47: === GModelSky ===
2018-11-07T11:19:47:  Name ......................: Crab
2018-11-07T11:19:47:  Instruments ...............: all
2018-11-07T11:19:47:  Test Statistic ............: 1845.6796865592
2018-11-07T11:19:47:  Instrument scale factors ..: unity
2018-11-07T11:19:47:  Observation identifiers ...: all
2018-11-07T11:19:47:  Model type ................: PointSource
2018-11-07T11:19:47:  Model components ..........: "PointSource" * "PowerLaw" * "Constant" 
2018-11-07T11:19:47:  Number of parameters ......: 6
2018-11-07T11:19:47:  Number of spatial par's ...: 2
2018-11-07T11:19:47:   RA .......................: 83.6209669517135 +/- 0.00252919459519199 deg (free,scale=1)
2018-11-07T11:19:47:   DEC ......................: 22.0248642601462 +/- 0.00233953943144988 deg (free,scale=1)
2018-11-07T11:19:47:  Number of spectral par's ..: 3
2018-11-07T11:19:47:   Prefactor ................: 5.23287376743654e-17 +/- 3.38933067177766e-18 [1e-25,infty[ ph/cm2/s/MeV (free,scale=1e-17,gradient)
2018-11-07T11:19:47:   Index ....................: -2.75931417445264 +/- 0.0757563212467188 [10,-10]  (free,scale=-2,gradient)
2018-11-07T11:19:47:   PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
2018-11-07T11:19:47:  Number of temporal par's ..: 1
2018-11-07T11:19:47:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2018-11-07T11:19:47: === GCTAModelCubeBackground ===
2018-11-07T11:19:47:  Name ......................: BackgroundModel
2018-11-07T11:19:47:  Instruments ...............: CTA, HESS, MAGIC, VERITAS
2018-11-07T11:19:47:  Instrument scale factors ..: unity
2018-11-07T11:19:47:  Observation identifiers ...: all
2018-11-07T11:19:47:  Model type ................: "PowerLaw" * "Constant" 
2018-11-07T11:19:47:  Number of parameters ......: 4
2018-11-07T11:19:47:  Number of spectral par's ..: 3
2018-11-07T11:19:47:   Prefactor ................: 0.999715699285799 +/- 0.0148947550844672 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)
2018-11-07T11:19:47:   Index ....................: 0.00566705715270537 +/- 0.014450942963835 [-5,5]  (free,scale=1,gradient)
2018-11-07T11:19:47:   PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
2018-11-07T11:19:47:  Number of temporal par's ..: 1
2018-11-07T11:19:47:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

I also checked the results obtained with a joint binned analysis. For this purpose I created for each observation a counts cube using the usepnt=yes option. Each cube had 250 x 250 spatial bins with a bin size of 0.02 deg. 20 energy bins were used. Below the results for unbinned, stacked and joint binned analysis.

I inspected the GCTAEdispCube::operator() interface and revised it according to the GCTAEdisp::operator() interface. Now the GCTAEdispCube::operator() takes the reconstructed and true energy in form of GEnergy objects, followed by the true sky direction.

I also corrected the GCTAEdispCube::set() method which performed some incorrect exposure weighting. This has been fixed. I added unit tests to make sure that the energy dispersion cube when integrated over reconstructed energy for a given set of true energies is unity. The tests are performed for the GCTAEdispCube::set() and GCTAEdispCube::fill() methods, hence we can be sure that the methods actually work.

I now also redid the analysis with the new stacked response cubes without cuts. For some reason the results differs from the former result Not clear what happened.

I therefore progressively put back the true energy cuts to check what’s going on. Without energy dispersion, cutting in true energy should be more realistic since it avoids spill-over in energy.

The impact of the PSF cut is negligible, which is understandable since the exact weight of the PSF between the various observations should not be relevant. The important cut is the Aeff cut. Cutting on both gives similar results to the Aeff cut alone, which is expected since the PSF cut is negligible.

Still, we don’t get exactly the same results are before, which may be due a different energy binning used for the PSF.

• spill-over leads to a bias in the prefactor of -4.8% and in the index of -1.7% without energy dispersion
• spill-over leads to a bias in the prefactor of -6.2% and in the index of -1.3% with energy dispersion (reference = unbinned)

I compared the model counts computed using an unbinned observation to the model counts computed using a stacked observation to study the effect of the spill over. The figures below show the residual of the stacked model minus the joint binned model. Left panels show the spectral residuals, central panels the spatial residual profiles in Right Ascension and Declination, and the right panels a residual map integrated over all energies. The top panels show the residuals without clipping, the bottom panel show the residuals with clipping. Obviously, clipping leads to a strong residual at low energies, near the threshold of the observations. Not clipping the stacked response is therefore strongly preferred.

The non-clipped residuals show an energy dependent trend, where only at high energies, the computation is identical. This indicates that the stacked PSF cube may not fully comprise the H.E.S.S. PSF. I therefore increased the maximum PSF angle from 0.2 degrees to 0.7 degrees which covers the angles provided in the H.E.S.S. IRF. I also set the number of angular bins to 300, which gives an ovsersampling by a factor of 2 of the PSF. The image below shows, obtained without clipping, that this solves the residual trend.

Here the plots without and with energy dispersion for no clipping:

And now with clipping enabled:

From these plots it is obvious that the stacked response should not be clipped. I therefore changed the code.

For the record, attached the script to generate the spill over plot: show_spill_over.py

Below the spectra obtained using an unbinned, joint binned and stacked analysis. Note that the energy bins for the binned and stacked analysis are not aligned with the bin boundaries of the counts cube, which explains the gaps in the spectrum. The joined binned and stacked analysis results are very similar.

Method	No energy dispersion	Energy dispersion
Unbinned
Binned
Stacked

I consider that the stacked analysis is now checked.