ctools

First I tried to reproduce Masha’s results and I managed to do so:

2021-12-10T16:11:40: +=========================================+
2021-12-10T16:11:40: | Maximum likelihood optimisation results |
2021-12-10T16:11:40: +=========================================+
2021-12-10T16:11:40: === GOptimizerLM ===
2021-12-10T16:11:40:  Optimized function value ..: 348541.336
2021-12-10T16:11:40:  Absolute precision ........: 0.005
2021-12-10T16:11:40:  Acceptable value decrease .: 2
2021-12-10T16:11:40:  Optimization status .......: converged
2021-12-10T16:11:40:  Number of parameters ......: 231
2021-12-10T16:11:40:  Number of free parameters .: 4
2021-12-10T16:11:40:  Number of iterations ......: 1
2021-12-10T16:11:40:  Lambda ....................: 0.0001
2021-12-10T16:11:40:  Total number of iterations : 8
2021-12-10T16:11:40:  Maximum log likelihood ....: -348541.336
2021-12-10T16:11:40:  Observed events  (Nobs) ...: 52107.000
2021-12-10T16:11:40:  Predicted events (Npred) ..: 51863.589 (Nobs - Npred = 243.411000453176)
...
2021-12-10T16:11:40: === GModelSky ===
2021-12-10T16:11:40:  Name ......................: R000085_00001_J1737-3837
2021-12-10T16:11:40:  Instruments ...............: all
2021-12-10T16:11:40:  Test Statistic ............: 34.7896858851891
2021-12-10T16:11:40:  Observation identifiers ...: all
2021-12-10T16:11:40:  Model type ................: PointSource
2021-12-10T16:11:40:  Model components ..........: "PointSource" * "ExponentialCutoffPowerLaw" * "Constant" 
2021-12-10T16:11:40:  Number of parameters ......: 7
2021-12-10T16:11:40:  Number of spatial par's ...: 2
2021-12-10T16:11:40:   RA .......................: 264.370455982881 deg (fixed,scale=1)
2021-12-10T16:11:40:   DEC ......................: -38.6193821817504 deg (fixed,scale=1)
2021-12-10T16:11:40:  Number of spectral par's ..: 4
2021-12-10T16:11:40:   Prefactor ................: 9.37556142183426e-19 +/- 2.55043578929073e-19 [0,infty[ ph/cm2/s/MeV (free,scale=1e-18,gradient)
2021-12-10T16:11:40:   Index ....................: -1.93913156282528 [-9,0]  (fixed,scale=-1,gradient)
2021-12-10T16:11:40:   CutoffEnergy .............: 2194778.42864015 [10000,500000000] MeV (fixed,scale=1000000,gradient)
2021-12-10T16:11:40:   PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
2021-12-10T16:11:40:  Number of temporal par's ..: 1
2021-12-10T16:11:40:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2021-12-10T16:11:40:  Number of scale par's .....: 0

Then I tried fitting the data of the first light curve bin with different combinations of model components. In all cases, the predicted number of counts correspond to the observed one. Even with only a few components, PWN 2059 got a low flux, comparable to Masha’s run.

Here are smoothed background subtracted sky maps of the two time bins 192 and 201 (files gps_standard_000192.fits and gps_standard_000201.fits). Obviously PWN 2059 appears much brighter in the second time bin compared to the first time bin.

In the event file, PWN 2059 has the Monte Carlo identifier 1242. In the first time slice there are 50 Monte Carlo events from PWN 2059, while in the second time slice there are 86 events, hence a factor of 1.72 more. A large fraction of the flux discrepancy can therefore be explained by the different number of simulated events.

Here a summary of the information on the number of simulated photons and events found in the ctobssim log files for the corresponding event files. What is shown is “Number of simulated events / Number of simulated photons”. These numbers are not fully identically to the number of events determined from the event files due to the 0.1-100 TeV energy selection. The original gps_standard_000192.fits file contained 102933 events, the energy selected file contains 52107 events. The original gps_standard_000201.fits file contained 101877 events, the energy selected file contains 51365 events.

In the 3rd to 6th energy slices, the number of simulated events is systematically larger in the 201 time bin compared to the 192 time bin. Maybe this is just “bad luck”?

The two leftmost plots are so far only for 324 simulations. I will try to make 1000 simulations to verify the match with the expectation, yet already with 324 simulations there seems to be a reasonable match, and no significant difference to the faster simulations where only PWN 2059 is used.

I investigated whether the observed low number of events is explained by pure statistical fluctuations, or whether there is a possible issue with the random number generator.

To do this I run the simulations using ctools version 1.7.0 by changing the seed parameter of ctobssim which sets the seed value of the random number generator, and extract the number of photons and events that are simulated for PWN 2059. Repeating the process 1000 times for different seed values gives then a distribution of number of photons and events that can be compared to the expectation based on statistical fluctuations only.

The outcome of this process are show below, with the left two panels corresponding to the simulation of the full GPS model, reproducing exactly the simulation that was done for the GPS, while the right two panels correspond to simulations of only PWN 2059, speeding up considerably the simulations. The table below gives the moments of the photon and event distributions. The observed number of 53 events for PWN 2059 in run 192 is a deviation of -18.89 events from the mean, corresponding to -2.01 standard deviations.

Obviously, the resulting histograms correspond to the expectations for statistical fluctuations only, hence the low flux of PWN 2059 in the first time bin can be explained by a statistical fluctuation of about -2 sigma.

As indicated at https://cta.cloud.xwiki.com/xwiki/wiki/phys/view/Galactic%20WG/Galactic%20Plane%20Survey/BinariesAnalysis/ the average prefactor of PWN 2059 is about 2.25-18 ph cm-2 s-1 MeV-1 while the fitted result for run 192 is 9.38e-19 ± 2.55e-19, corresponding to an offset of -5.1 sigma.

So the question is why a statistical fluctuation of about -2 sigma translates into a prefactor offset of -5.1 sigma?

Interestingly, 53 counts is only a factor of 0.74 lower than 72 counts, in other words, the mean prefactor expected for run 192 would be 9.38e-19/0.74 = 1.27e-18 which is almost a factor of 2 below the mean value found by Masha for the other bins.

Below is the result for 10000 simulations of PWN 2059 only. The moments of the distributions are given in the Table below. The observed number of 53 events for PWN 2059 in run 192 is a deviation of -18.44 events from the mean, corresponding to -2.18 standard deviations.

Below the results of fitting various models to the runs 192 and 201. For illustration, the fitted SEDs for the exponential cut-off model (left) and of the MC file function model (right) are shown in the two panels below.

Exponential cut-off

File function

Note that while the model fit is driven by the events which have the largest statistical significance (which are potentially events at high energies), the overall event statistics is driven by the low-energy events which have less weight in the likelihood computation. This may explain the apparent discrepancy between the event histograms and the model fit results.

I now investigated the distribution of events with changing random number generator seed for energies above 2.44 TeV, where in run 192 2 events are simulated while in run 201 8 events are simulated. The distribution for 1000 simulations are shown in the plots below, moments of the distribution are given in the following table. While 2 events are on the left tail of the distribution, 8 events are on the right tail of the distribution. Note that for some reason, 8 out of 1000 seed values produce zero photons (seed values 81, 192, 222, 294, 465, 557, 703, and 888).

Repeating the simulations for ctools-2.0.0.dev does not show the zero photon bins. The mean photon number is slightly larger, its standard deviation is slightly reduced. Not clear what the difference between ctools-1.7.0 and ctools-2.0.0.dev regarding simulations is.

I now removed the PWN 2059 events from the simulated events for run 192 and 201 using ctselect with the parameter expr="MC_ID!=1242". This reveals the background underlying the source events. Sky maps using the selected events are shown in the panels below, left is for run 192, right is for run 201. Apparently, there is a lack of background events for run 192 in the region of PWN 2059 which may explain why the fitted flux is so low.

Unfortunately fitting of sources with negative flux does not work (probably due to internal positivity constraints), hence I could not determine the flux level to which the background hole corresponds.

• no anomaly was detected, the data are as expected based on statistics
• the number of simulated events for run 192 is on the low side, but this number is compatible with expectations from statistical fluctuations
• there is a negative background fluctuation underlying PWN 2059 that reduces the number of events attributed to the source

So the low flux point is produced by a coincidence between a small number of events simulated for PWN 2059 and a hole in the background at the location of PWN 2059 for run 192. Looks like the low flux point is really just bad luck!