Bug #2861

Discrepancy between simulated counts and model for radial source

Added by Tiziani Domenico over 3 years ago. Updated almost 2 years ago.

Status:ClosedStart date:03/26/2019
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.7.0
Duration:

Description

If I simulate a source with a GModelSpatialRadialDisk spatial model with IRFs from a run of the H.E.S.S. public data release 1, I experience discrepancies in the number of events between simulated counts and modelled counts.

The simulated number of counts is 91389 and the number of model counts is 85285. I simulate only the source and no background.
The residual map (counts-model) does also not look empty (See attached image).

I also attached a script for reproducing the simulation.

residual.png - Residual map (46 KB) Tiziani Domenico, 03/26/2019 10:53 AM

simulation.py Magnifier - Complete script (4.61 KB) Tiziani Domenico, 03/26/2019 10:54 AM

model.png - Model map (29.6 KB) Tiziani Domenico, 03/26/2019 12:21 PM

model2.png - Model map after fix 2860 (26.1 KB) Tiziani Domenico, 03/27/2019 02:01 PM

Residual Model Model2

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen over 3 years ago

Try taking the model cube definition from the counts cube via the incube parameter. This should make sure that the correct bin weighting is applied.

#2 Updated by Tiziani Domenico over 3 years ago

Knödlseder Jürgen wrote:

Try taking the model cube definition from the counts cube via the incube parameter. This should make sure that the correct bin weighting is applied.

I have tried this but it does not change the number of model events.

#3 Updated by Tiziani Domenico over 3 years ago

Something else that I noticed is that the model map shows some strange ring shapes. This seems to originate from the IRF convolution but doesn’t look physically sensible if I am not missing something.

Model map

#4 Updated by Knödlseder Jürgen over 3 years ago

Have a look at #2860, I just increased the integration precision for extended source models since it was apparently not sufficient for the broad H.E.S.S. PSF. Could your issue be related to this? New code is in the devel branch since yesterday evening.

#5 Updated by Tiziani Domenico over 3 years ago

  • File model2.png added

In my test scenario, the fix in #2860 changes the number of counts in the model cube to 89149. (In the unbinned case: Nobs - Npred = 93267.000 - 93322.253 = -55.253)
The model map now looks much smoother:
Model map after fix 2860

#6 Updated by Tiziani Domenico over 3 years ago

  • File deleted (model2.png)

#7 Updated by Tiziani Domenico over 3 years ago

Updated image.

#8 Updated by Knödlseder Jürgen over 3 years ago

  • Tracker changed from Support to Bug
  • Status changed from New to In Progress
  • Target version set to 1.7.0
  • % Done changed from 0 to 10

I looked a bit in the problem and played around with the integration precisions in GCTAResponseIrf::irf_radial. Below a summary of the results, obtained using a modified simulation.py script where the cube energy range was limited to 0.5-90 TeV (to make sure to stay away from the threshold). I tried a variant where the maximum delta of the PSF was restricted to 99% of the PSF containment, but this led to some artefacts in the model cube at high energies. Setting the number of iterations in rho and phi to 6, as done in #2860, reduces the discrepancy between model and counts to about 3%. 1% is only reached when further increasing the integration precision in phi, but this brings the computation time to 10 min instead of 1 min when compared to the initial precision. In other words, the speed penalty would be a factor of 10, while with 6/6 it would be a factor of 4, which is already quite bad.

rho phi Nobs Nmodel ctmodel time delta_max model cube aspect
5 5 84169 78538 (-6.7%) 1m12sec 99% of PSF Artefacts in the model cube at high energies
6 6 84169 82332.2 (-2.2%) 4m09sec 99% of PSF Artefacts in the model cube at high energies
6 6 84132 81336 (-3.3%) 5m03sec Full range Still faint ring-like structures at high energies
7 6 84132 81355.4 (-3.3%) 9m44sec Full range Model looks now very smooth, extremely little radial structure
6 7  84132 83359.9 (-0.9%) 9m37sec Full range Still faint ring-like structures at high energies
7 7 84132 83381.6 (-0.9%) 19m09sec Full range Almost perfectly smooth
5 7 84132 83625.4 (-0.6%) 4m53sec Full range Strong ring-like structures over all energies

I prefer for the time being (and release 1.6) to fix the precision at 6/6, and defer it to the next release to think about a more precise but fast integration scheme. I keep this issue open for that purpose.

#9 Updated by Tiziani Domenico over 3 years ago

I agree, the computation time of ctmodel should be kept short.
Just to keep in track of it: We should also not forget about GCTAResponseIrf::irf_elliptical, where I see similar problems when going to large semi-axes.

#10 Updated by Knödlseder Jürgen over 3 years ago

Tiziani Domenico wrote:

I agree, the computation time of ctmodel should be kept short.
Just to keep in track of it: We should also not forget about GCTAResponseIrf::irf_elliptical, where I see similar problems when going to large semi-axes.

That’s the method where I discovered the issue. I changed the precision in both methods.

#11 Updated by Knödlseder Jürgen almost 2 years ago

  • Status changed from In Progress to Closed
  • Assigned To set to Knödlseder Jürgen
  • % Done changed from 10 to 100

Also available in: Atom PDF