fermi-lat / Fermitools-conda

Conda recipe files for the Fermi Sciencetools software analysis package: Fermitools
https://fermi.gsfc.nasa.gov/ssc/data/analysis/
BSD 3-Clause "New" or "Revised" License
34 stars 16 forks source link

Floating point exception (core dumped) in gtsrcmap #138

Closed mireianievas closed 1 year ago

mireianievas commented 1 year ago

I'm running some Crab Nebula analysis with the following settings:

Fermitools_conda version: 2.2.12 (freshly installed 2 days ago). Photon list: weekly event files Time: MJD59169.0 to MJD59176.0, (1 week of data) Energy range: 100 MeV to 1 TeV (10 energy bins / decade) Event selection: evclass = 128 and evtype for selection = 3. zmax = 90. Summed likelihood with FRONT and BACK

I'm using the Fermitools frontend 'enrico' to run the different stages of the analysis and handle the summed likehood fit.

At gtsrcmaps, running it like:

gtsrcmaps scfile=/net/diva/scratch1/mnievas/enrico//Data/download/lat_spacecraft_merged.fits sctable="SC_DATA" expcube=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONT_ltCube.fits cmap=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONT_CCUBE.fits srcmdl=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LogParabola_model_FRONT.xml bexpmap=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONT_BinnedMap.fits wmap=none outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONT_LogParabola_srcMap.fits.tmpout irfs="P8R3_SOURCE_V3" evtype=1 convol=yes resample=yes rfactor=2 minbinsz=0.1 ptsrc=yes psfcorr=yes emapbnds=no edisp_bins=-2 copyall=no chatter=4 clobber=yes debug=no gui=no

I get the following error

This is gtsrcmaps version HEAD
ResponseFunctions::load: IRF used: P8R3_SOURCE_V3
  event_types:  0
Creating source named IsoDiffModel
Creating source named GalDiffModel
Creating source named Crab
Floating point exception (core dumped)

I never saw this error before.

The event file (before any binning, gti application, etc) contains some events, but there are parts of the image with 0 counts. I guess that is ok, or maybe ... Did for that specific week part of the image have by any chance zero exposure?

image

I am trying again without dividing FRONT and BACK (maybe if the exposure is a bit at the limit it is sorted out this way?). If the problem persists, where can I upload the files to be investigated?

nmirabal commented 1 year ago

Unfortunately, the FSSC does not maintain Enrico so we cannot trace the whole analysis chain. Can you run the following analysis thread to see if the error comes up? https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html

mireianievas commented 1 year ago

Hi @nmirabal

enrico just executes the standard fermi tools, at least until the likelihood part. The only thing that changes is how the energy/time binning is handled for SEDs and LCs, but that happens at a later point.

I tried to fetch directly from the photon server the photon and spacecrafts, to make sure it is not something to do with some buggy/corrupt weekly file from my side. Same error happens at the gtsrcmaps stage.

Please find attached the directory with the photon file, the spacecraft file and the xml file

https://drive.google.com/file/d/1rhRfNhMgdRho1Fi_SC-ne4H5UXkoSs06/view?usp=share_link

These are the commands that are executed:

gtselect (coarse cuts)

gtselect infile=/net/diva/scratch1/mnievas/Crab_asgardpy/event.list outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_EvtCoarse.fits.tmpout ra=83.6333333334 dec=22.0133333333 rad=15.0 tmin=627177601.0 tmax=627782401.0 emin=0.0 emax=10000000.0 zmin=0.0 zmax=90.0 evclass=128 evtype="INDEF" convtype=-1 phasemin=0.0 phasemax=1.0 evtable="EVENTS" chatter=2 clobber=yes debug=no gui=no mode="ql"

gtselect (finer cuts)

gtselect infile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_EvtCoarse.fits outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_Evt.fits.tmpout ra=0.0 dec=0.0 rad=180.0 tmin="INDEF" tmax="INDEF" emin=50.11872336272722 emax=1995262.3149688789 zmin=0.0 zmax=90.0 evclass=128 evtype=3 convtype=-1 phasemin=0.0 phasemax=1.0 evtable="EVENTS" chatter=2 clobber=yes debug=no gui=no mode="ql"

gtmktime

gtmktime scfile=/net/diva/scratch1/mnievas/Crab_asgardpy/L230112061201F357373F89_SC00.fits sctable="SC_DATA" filter="(DATA_QUAL>0)&&(LAT_CONFIG==1)" roicut=no evfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_Evt.fits evtable="EVENTS" outfile="/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_MkTime.fits.tmpout" apply_filter=yes overwrite=no header_obstimes=yes tstart=627177601.0 tstop=627782401.0 gtifile="default" chatter=2 clobber=yes debug=no gui=no mode="ql"

gtbin

gtbin evfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_MkTime.fits scfile=/net/diva/scratch1/mnievas/Crab_asgardpy/L230112061201F357373F89_SC00.fits outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_CountMap.fits.tmpout algorithm="CMAP" ebinalg="LOG" emin=50.11872336272722 emax=1995262.3149688789 enumbins=12 denergy=0.0 ebinfile=NONE tbinalg="LIN" tstart=239557418.0 tstop=677656293.699315 dtime=0.0 tbinfile=NONE snratio=0.0 lcemin=0.0 lcemax=0.0 nxpix=212 nypix=212 binsz=0.1 coordsys="CEL" xref=83.6333333334 yref=22.0133333333 axisrot=0.0 rafield="RA" decfield="DEC" proj="AIT" hpx_ordering_scheme="RING" hpx_order=3 hpx_ebin=yes hpx_region="" evtable="EVENTS" sctable="SC_DATA" efield="ENERGY" tfield="TIME" chatter=2 clobber=yes debug=no gui=no mode="ql"

gtltcube

gtltcube evfile="/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_MkTime.fits" evtable="EVENTS" scfile=/net/diva/scratch1/mnievas/Crab_asgardpy/L230112061201F357373F89_SC00.fits sctable="SC_DATA" outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_ltCube.fits.tmpout dcostheta=0.025 binsz=1.0 phibins=0 tmin=0.0 tmax=0.0 file_version="1" zmin=0.0 zmax=90.0 chatter=2 clobber=yes debug=no gui=no mode="ql"
Working on file /net/diva/scratch1/mnievas/Crab_asgardpy/L230112061201F357373F89_SC00.fits

gtbin

gtbin evfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_MkTime.fits scfile=/net/diva/scratch1/mnievas/Crab_asgardpy/L230112061201F357373F89_SC00.fits outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_CCUBE.fits.tmpout algorithm="CCUBE" ebinalg="LOG" emin=50.11872336272722 emax=1995262.3149688789 enumbins=46 denergy=0.0 ebinfile=NONE tbinalg="LIN" tstart=627177601.0 tstop=627782401.0 dtime=0.0 tbinfile=NONE snratio=0.0 lcemin=0.0 lcemax=0.0 nxpix=212 nypix=212 binsz=0.1 coordsys="CEL" xref=83.6333333334 yref=22.0133333333 axisrot=0.0 rafield="RA" decfield="DEC" proj="AIT" hpx_ordering_scheme="RING" hpx_order=3 hpx_ebin=yes hpx_region="" evtable="EVENTS" sctable="SC_DATA" efield="ENERGY" tfield="TIME" chatter=2 clobber=yes debug=no gui=no mode="ql"

gtexpcube2

gtexpcube2 infile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_ltCube.fits cmap=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_CCUBE.fits outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_BinnedMap.fits.tmpout irfs="P8R3_SOURCE_V3" evtype=3 edisp_bins=-1 nxpix="INDEF" nypix="INDEF" binsz="INDEF" coordsys="CEL" xref="INDEF" yref="INDEF" axisrot=0.0 proj="AIT" ebinalg="LOG" emin=50.11872336272722 emax=1995262.3149688789 enumbins=46 ebinfile="NONE" hpx_ordering_scheme="RING" hpx_order=6 bincalc="EDGE" ignorephi=no thmax=180.0 thmin=0.0 table="EXPOSURE" chatter=2 clobber=yes debug=no mode="ql"

gtsrcmaps

gtsrcmaps scfile=/net/diva/scratch1/mnievas/Crab_asgardpy/L230112061201F357373F89_SC00.fits sctable="SC_DATA" expcube=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_ltCube.fits cmap=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_CCUBE.fits srcmdl=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LogParabola_model.xml bexpmap=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_BinnedMap.fits wmap=none outfile=/net/diva/scratch1/mnievas/Crab_asgardpy/Crab_LAT_Analysis_FRONTBACK_LogParabola_srcMap.fits.tmpout irfs="P8R3_SOURCE_V3" evtype=3 convol=yes resample=yes rfactor=2 minbinsz=0.1 ptsrc=yes psfcorr=yes emapbnds=no edisp_bins=-2 copyall=no chatter=2 clobber=yes debug=no gui=no mode="ql"
mireianievas commented 1 year ago

It seems to work just fine with other fields and other periods, so it seems to be related with something in that Crab field or in that period of time (MJD59169.0 to MJD59176.0)

nmirabal commented 1 year ago

So I ran your photon files following https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html and everything work properly for the Crab in that time interval. It is true that enrico runs the Fermitools commands but the value of the input in each command is critical as the Fermitools expect energies, size selection and time interval to always match. That could be one source of the problem. The other is that I noticed that in your isotropic and Galactic diffuse paths in your XML file had an extra slash / between enrico//Data and that could be causing the whole problem. I would start changing that. If that doesn't work, then check the enrico inputs step by step so that they match the example in https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html

mireianievas commented 1 year ago

Found the issue. There was a 'nan' in my xml file.

I wonder if the Fermitools could have caught this simple issue (invalid values in the xml model parameters) and perhaps raise a more meaningful error, as Floating point exception (core dumped) was not very helpful and increasing the verbosity and switching debug to True did not help much either.

mireianievas commented 1 year ago

And found the underlying issue. enrico's xml maker takes the flux from the PL_Flux_Density field (from the 4FGL) to set the normalization of the spectral component. Well, for Crab that returns a NaN for some reason ... Only the LP_Flux_Density has a floating value.

nmirabal commented 1 year ago

Glad that you found the error. The FSSC recommendation is to use make4FGLxml.py for XML generation. Enrico is not maintained by the FSSC.