Bug #2708
Problems with CR simulation
| Status: | Closed | Start date: | 11/01/2018 | |
|---|---|---|---|---|
| Priority: | Normal | Due date: | ||
| Assigned To: |  Knödlseder Jürgen | % Done: | 100% | |
| Category: | - | |||
| Target version: | 1.6.0 | |||
| Duration: |
Description
We have observed a strange behaviour in the number of counts as a function of energy in the data simulated with ctools.
We simulated CR only with the exposure time T_{obs} = 50 h in 2.5 degrees and observed that the number of simulated photons has a strange sharp kink at 0.1 TeV (see the attached file “Nobs_2p5deg.pdf”). Naively we would expect to have a power law. However at small energies (from 60 to 100 GeV) the number of photons increases and after 100 GeV it drops. Our first guess was that there is some feature in IRFs, however, we checked that IRFs have a smooth behaviour at least for small enough RA (“irf.pdf”). We made a further test and repeated the same simulations inside the 1.5-degree ring (“Nobs_1p5deg.pdf”) where IRF is very smooth. We see that in this case instead of kink we have a plato that also seems to be not physical. We checked that this strange behaviour appears for a point source as well (crab) but the effect is a bit softer.
Another way to see a strange behaviour at the same energy:
Again, we simulate CR only. If we fit this data with the same model (CR only), the residuals of such a fit have a feature for the same bin where we have the kink in N(E) (0.1 TeV) (“residuals.pdf”). In the bin containing 0.1 TeV the residual is especially bad.
We also checked that it exists also in the output of ctmodel (“Nmod.pdf”)
irf.pdf
(23.6 KB)
Nmod.pdf
(35.3 KB)
Nobs_1p5deg.pdf
(35.2 KB)
Nobs_2p5deg.pdf
(35.2 KB)
residuals.pdf
(39.4 KB)
residual-corrected-code.png (21.6 KB)
Recurrence
No recurrence.
History
#1
Updated by Knödlseder Jürgen about 6 years ago
- File residual-corrected-code.png added
- Subject changed from problems with CR simulation to Problems with CR simulation
- Status changed from New to Pull request
- Target version set to 1.6.0
- % Done changed from 0 to 100
The feature does not disappear when changing the ctobssim seed, which means that it is not a statistical fluctuation. Either the model computation is not exact, or the Monte Carlo simulation is the issue.
I investigated first the model computation by changing the integration precision in GObservation::npred_spec. As the table shows below, the Npred computation changes only marginally when increasing the integration precision. Note that the simulated number of events is 338218.
| iter | Npred | Nobs - Npred |
| 8 (default) | 340114.616 | -1896.616 |
| 9 | 340114.889 | -1896.889 |
| 10 | 340114.895 | -1896.895 |
I then investigated the number of energy bins that is used in the Monte-Carlo pre-computation cache in GCTABackground3D::init_mc_cache. This changes the number of simulated events significantly, bringing the number of observed events closer to the predicted number of events with an increasing number of energy bins. It is interesting to notice that going from 500 to 600 energy bins reduced the number of observed events, indicating that the precise number of energy bins influences the result.
| nbins | Nobs |
| 210 | 338218 |
| 400 | 338932 |
| 500 | 339616 |
| 600 | 339224 |
| 1000 | 339645 |
| 2000 | 339627 |
After looking into the code, it turns out that the actual node energies of the background cube where not necessarily used in the Monte-Carlo pre-computation cache, meaning that sharp features in energy could be missed. I therefore changed the code in GCTABackground3D::init_mc_cache so that a number of sub bins is inserted between the energy nodes of the background cube, making sure that the computation is always done for all energy nodes. I selected 10 sub bins, which resulted in the above example in a total of 189 energy bins, i.e. less bins than before. The number of observed events is now 339791, very close to the predicted number of 340114.616.
The figure below shows the residual plot after correcting the code. The outlier noticed earlier has disappeared.
#2
Updated by Knödlseder Jürgen about 6 years ago
- Status changed from Pull request to Closed
Code merged into devel.
#3
Updated by Knödlseder Jürgen about 6 years ago
- Project changed from ctools to GammaLib
- Target version deleted (
1.6.0)
#4
Updated by Knödlseder Jürgen about 6 years ago
- Target version set to 1.6.0