ctools

*******************************
* ctools science verification *
*******************************
Test background model: ..... ok
Test power law model: ......... ok
Test power law model with energy dispersion: F........ NOK
Test power law model with stacked analysis: FFF.F.F.. NOK
Test power law model with On/Off analysis: E NOK
Test power law 2 model: ......... ok
Test smoothly broken power law model: ............. ok
Test exponentially cut off power law model: ........... ok
Test super exponentially cut off power law model: ............. ok
Test log parabola model: ........... ok
Test Gaussian model: ........... ok
Test file function model: ....... ok
Test nodes model: ............. ok
Test exponential model: ........... ok
Test point source model: ............. ok
Test radial disk model: ............... ok
Test radial Gaussian model: ............... ok
Test radial shell model: ................. ok
Test elliptical disk model: 

The On/Off analysis simply fails, probably due to the invalid instrument direction, as indicated above.

It turned out that the On/Off problem was related to missing information in the observation definition XML file, such as the energy range and the region of interest. After adding the information, the science verification works:

Test power law model with On/Off analysis: ......... ok

I now investigated the power law fit with the energy dispersion enabled. Switching off the smoothing of the energy dispersion resolves the issue. The plots below show the Prefactor pull distribution with smoothing of the energy dispersion (left) and without smoothing of the energy dispersion (right).

I then simulated 5 hours of Crab observation using the Prod2 IRF with energy dispersion and checked the results when fitting without and with energy dispersion enabled. The test were done for smoothing enabled or disabled. Events were simulated for reconstructed energies between 0.1-50 TeV.

The analysis confirms the trend that with a smoothed IRF the Crab fluxes are underestimated. With smoothing enabled, the true energy range for a given reconstructed energy range becomes broader, leading to more photons that are simulated. Specifically, with smoothing enabled 608490 Crab photons were thrown while with smoothing disabled, 441251 Crab photons were thrown. The number of Crab events with smoothing were 37807, the number without smoothing was 37843.

Independently of what is right and what is wrong, the question is why there is a mismatch between simulations and Monte Carlo when smoothing is applied to the energy dispersion.

The figures below show the original (left) and smoothed (right) energy dispersion for the Prod2 South_5h IRF. The original IRF is irregular, but the smoothed IRF certainly poorly captures what happens at the low energy end. Below the 2D figures a comparison between the MC simulation and the model (left panels) and the true energy PDF (right panel). Not obvious why the smoothed IRF does not correctly reproduce the simulations.

I created a script to compare the simulated event distribution to the modeled one (show_model_mc.py).

Below the results for a run without energy dispersion (top panel), with energy dispersion but no smoothing (mid panel), and with energy dispersion and smoothing (bottom panel). The simulation was done for 50h of the Crab nebula, Prod2 South_50h IRF.

Without energy dispersion, everything looks reasonable.

With energy dispersion but no smoothing, there are spectral outliers. Such spectral outliers were also observed for the 1DC data, and for exactly that reason, the smoothing was introduced.

With energy dispersion and smoothing the spectral outliers are gone, but now the low energy points are negative, i.e. there are more modeled than simulated events.

I checked whether the integration precision in GResponse::convolve has an impact on the result by setting iter=7 instead of the default iter=6. The top panel below shows the result for iter=7, the bottom panel for iter=6. No significant difference is observed.

So far the energy dispersion integration is done in logarithmic energy, but since now everything is per linear energy the integration should probably also be done in linear energy. I changed the code in GResponse::convolve to linear energy and run again the computations for iter=6. Below the result (top panel) in comparison to the logarithmic integration (bottom panel). Things look much better now.

After having changed also the energy dispersion integration in GCTAResponseIrf::nroi and cta_nroi_kern::eval to linear (this is needed for an unbinned analysis), the science verification succeeds:

Test power law model with energy dispersion: ......... ok

• ctbin no longer computed the average pointing direction for a stacked cube if usepnt=yes. I added some code so that the average pointing direction is computed as before.
• obsutils.sim() did not pass a counts cube to the response computation, and consequently the background cube was not correctly weighted. I changed this and added the change request #2763 that imposes usage of a counts cube

After these changes the science verification works:

Test power law model with stacked analysis: ......... ok

Below the resulting Crab SED. On the left without energy dispersion, on the right with energy dispersion.

Code was merged into devel.