sed() method not working properly with Mac OS

[ranieremenezes]

The GTAnalysis function sed() eventually returns a fits table with an empty TS column and very small values for dnde_err, e2dnde_err, dnde_ul, e2dnde_ul (and also some other columns). I also noticed that, in these cases, the column "dloglike_scan" contains NaNs.

I highlighted the word eventually above, because sometimes the exact same analysis works fine, although this is much rarer to happen. Furthermore, the same analysis works smoothly on Linux. The analysis I am testing is the following:

Mrk 421
time window: START, 14/10/2008 15:43:00
energy window: 100, 300000

There are no errors shown in the terminal after running the function sed().
I am using a macOS High Sierra version 10.13.6, but I got the exact same issue on a macOS Big Sur version 11.7.4 with chip Apple M1.

[omodei]

Can you provide code/configuration files and everything needed to reproduce thus bug? Whoever will look at this issue should be able to run a script, obtain the failure without having to download data/write a script. Thanks!

[ranieremenezes]

Hi @omodei ,

I cannot attach .py, .yaml, and .fits files here. So I will write the script and configuration file below.

Configuration file for Mrk 421:

data:
  evfile : /Users/xjet/Documents/Mrk421/fermipy_test/list.txt
  scfile : /Users/xjet/Documents/Mrk421/fermipy_test/spacecraft.fits

binning:
  roiwidth   : 15
  binsz      : 0.1
  binsperdec : 8

selection :
  emin : 100.0
  emax : 300000.0
  zmax    : 90
  evclass : 128
  evtype  : 3
  ra: 166.1138083333333
  dec: 38.20883277777778
  tmin: 239557417
  tmax: 245691781

gtlike:
  edisp : True
  irfs : 'P8R3_SOURCE_V3'
  edisp_disable : ['isodiff']
  edisp_bins : -2

model:
  src_roiwidth : 25
  galdiff  : '/Users/xjet/Documents/Mrk421/Diffuse/gll_iem_v07.fits'
  isodiff  : '/Users/xjet/Documents/Mrk421/Diffuse/iso_P8R3_SOURCE_V3_v1.txt'
  catalogs : ['4FGL']

Script to generate an SED for Mrk 421:

import numpy as np
from fermipy.gtanalysis import GTAnalysis

gta = GTAnalysis('config.yaml',logging={'verbosity': 3})
gta.setup()
gta.optimize(npred_threshold=100,shape_ts_threshold=50)
gta.write_roi('roi1')

gta.free_sources(distance=5.0,pars=['norm','shape'])
gta.free_source('galdiff')
gta.free_source('isodiff')
fit_results = gta.fit(optimizer='NEWMINUIT')
print('Fit Quality: ',fit_results['fit_quality'])
gta.write_roi('Results',make_plots=False)
resid = gta.residmap("4FGL J1104.4+3812")

c = np.load('Results.npy', allow_pickle=True).flat[0]
E = np.array(c['sources']["4FGL J1104.4+3812"]['model_flux']['energies'])
dnde = np.array(c['sources']["4FGL J1104.4+3812"]['model_flux']['dnde'])
dnde_hi = np.array(c['sources']["4FGL J1104.4+3812"]['model_flux']['dnde_hi'])
dnde_lo = np.array(c['sources']["4FGL J1104.4+3812"]['model_flux']['dnde_lo'])

ebins = np.linspace(2.0,5.477,num=12)
ebins = ebins.tolist()
sed = gta.sed("4FGL J1104.4+3812",loge_bins=ebins)

print("Energy center: ", sed['e_ctr'])
print("E2dNdE: ", sed['e2dnde'])
print("TS: ",sed['ts'])

The output from this script is something like:

Energy center: [1.43e+02 2.98e+02 ...]
E2dNdE: [1.45e-05 1.84e-05 ...]
TS: [nan nan nan ...]

However, we eventually get normal values of TS! I just did the test here and the script worked fine on the first time it was called. From the second try on, I always get "nan" values for the TS column in the sed.fits file.

The data selection information for the download is:
emin : 100.0
emax : 300000.0
zmax : 90
evclass : 128
evtype : 3
ra: 166.1138
dec: 38.2088
tmin: 239557417
tmax: 245691781

[omodei]

Hi @ranieremenezes, unfortunately I have found this issue other times, when running different fit one after the other. In your case, you have gta.fit followed by gta.residmap and by gta.sed. My impression is that internally the fit gets confused. This is obviously not good. A possible solution, is to rerun the setup before the final sed fit. In the following example, I setup again the gpa analysis, I set free the parameters forth fit. I do the fit, and after it, I fix all the parameters (the sed method will free the normalization of the interested source.

gta.setup()
gta.free_sources(distance=5.0,pars=['norm','shape'])
gta.free_source('galdiff')
gta.free_source('isodiff')
fit_results = gta.fit(optimizer='NEWMINUIT')
gta.fit()
gta.free_sources(free=False)
sed = gta.sed("4FGL J1104.4+3812",loge_bins=ebins,make_plots=True)
print("Energy center: ", sed['e_ctr'])
print("E2dNdE: ", sed['e2dnde'])
print("TS: ",sed['ts'])

This seemed to provide stable results. Let me know if it works. Cheers!

[omodei]

I don't see any more activity, so I guess this solution works. If not, feel free to reopen this issue.

[ranieremenezes]

Hi @omodei, thank you very much for finding the solution to this issue! I am quite busy these days, so I am not sure when I will have time to test it. In any case, I am glad to know we can go around this problem.