sergiohcdna / ctadmtool

Analysis scripts for DM cluster searches
Other
1 stars 3 forks source link

Management of `GCTAOnOffObservations` #12

Closed sergiohcdna closed 3 years ago

sergiohcdna commented 3 years ago

There is an exception when trying to use csdmatter with On/Off observations generated with csphagen. The problems seems to be related to the fact that the container of the models is missing(?)

sergiohcdna commented 3 years ago

I got the following error during the likelihood fitting using GModelSpectralTable:

*** ERROR: GModelSpectralPlaw::eval(srcEng=4.51695226161221 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, m_norm=nan, m_index=nan, m_pivot=1e+06, m_last_power=nan)

I attach the log-file obtained with ctlike ctlike.log

sergiohcdna commented 3 years ago

I checked if the same issue is present with other model tables. I used the following script to test a model table describing a PowerLaw with exponential cutoff. Description of the model table:

=== GModelSpectralTable ===
 Table file ................: model_table.fits
 Number of parameters ......: 3
  Normalization ............: 1 +/- 0 [0,1000]  (free,scale=1,gradient)
  Index ....................: -2 +/- 0  (free,scale=-2)
  Cutoff ...................: 1000000 +/- 0  (free,scale=1000000)
 Index values ..............: 100 [-3, -1.02] (linear)
 Cutoff values .............: 50 [100000, 28183829.3126445] (linear)
 Energies ..................: 50 [10 GeV, 300 TeV]
 Spectra array dimension ...: 3
 Number of spectra .........: 5000
 Number of spectral bins ...: 50

I used the following pipeline to generate the onoff observation container and perform the likelihood fitting

ctobssim --> csphagen --> ctlike

I got the exactly same error as above when I put cutoff values near the boundary of the value (1.1e+5, for example):

*** ERROR: GModelSpectralPlaw::eval(srcEng=2.65914794847249 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, m_norm=nan, m_index=nan, m_pivot=1e+06, m_last_power=nan)

*** ERROR: GModelSpectralTable::eval(srcEng=1.5524587917296 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, norm=nan, func=nan)
sergiohcdna commented 3 years ago

Checking the error above:

Using the table model for the powe-law with exponential cutoff. If I set the cutoff parameter to values near to the maximum value allowed, there is no problem creating and fitting an OnOff observation

sergiohcdna commented 3 years ago

I checked if the same issue is present with other model tables. I used the following script to test a model table describing a PowerLaw with exponential cutoff. Description of the model table:

=== GModelSpectralTable ===
 Table file ................: model_table.fits
 Number of parameters ......: 3
  Normalization ............: 1 +/- 0 [0,1000]  (free,scale=1,gradient)
  Index ....................: -2 +/- 0  (free,scale=-2)
  Cutoff ...................: 1000000 +/- 0  (free,scale=1000000)
 Index values ..............: 100 [-3, -1.02] (linear)
 Cutoff values .............: 50 [100000, 28183829.3126445] (linear)
 Energies ..................: 50 [10 GeV, 300 TeV]
 Spectra array dimension ...: 3
 Number of spectra .........: 5000
 Number of spectral bins ...: 50

I used the following pipeline to generate the onoff observation container and perform the likelihood fitting

ctobssim --> csphagen --> ctlike

I got the exactly same error as above when I put cutoff values near the boundary of the value (1.1e+5, for example):

*** ERROR: GModelSpectralPlaw::eval(srcEng=2.65914794847249 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, m_norm=nan, m_index=nan, m_pivot=1e+06, m_last_power=nan)

*** ERROR: GModelSpectralTable::eval(srcEng=1.5524587917296 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, norm=nan, func=nan)

It seems like the error is related to the interpolation calculation or the normalization value retrieved from the table model during the ctlike with the OnOff Observation.

Taking a look to the full output of ctlike, I found that in some (almost in the majority) of the cases, there is a result (different from NaN or infty) for the interpolated value, but the norm and value entries still have a NaN. Then checking the source code for GModelSpectralTable, looks like the norm parameter is causing the error.

sergiohcdna commented 3 years ago

I checked if the same issue is present with other model tables. I used the following script to test a model table describing a PowerLaw with exponential cutoff. Description of the model table:

=== GModelSpectralTable ===
 Table file ................: model_table.fits
 Number of parameters ......: 3
  Normalization ............: 1 +/- 0 [0,1000]  (free,scale=1,gradient)
  Index ....................: -2 +/- 0  (free,scale=-2)
  Cutoff ...................: 1000000 +/- 0  (free,scale=1000000)
 Index values ..............: 100 [-3, -1.02] (linear)
 Cutoff values .............: 50 [100000, 28183829.3126445] (linear)
 Energies ..................: 50 [10 GeV, 300 TeV]
 Spectra array dimension ...: 3
 Number of spectra .........: 5000
 Number of spectral bins ...: 50

I used the following pipeline to generate the onoff observation container and perform the likelihood fitting

ctobssim --> csphagen --> ctlike

I got the exactly same error as above when I put cutoff values near the boundary of the value (1.1e+5, for example):

*** ERROR: GModelSpectralPlaw::eval(srcEng=2.65914794847249 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, m_norm=nan, m_index=nan, m_pivot=1e+06, m_last_power=nan)

*** ERROR: GModelSpectralTable::eval(srcEng=1.5524587917296 TeV, srcTime=0 s (TT)): NaN/Inf encountered (value=nan, norm=nan, func=nan)

It seems like the error is related to the interpolation calculation or the normalization value retrieved from the table model during the ctlike with the OnOff Observation.

Taking a look to the full output of ctlike, I found that in some (almost in the majority) of the cases, there is a result (different from NaN or infty) for the interpolated value, but the norm and value entries still have a NaN. Then checking the source code for GModelSpectralTable, looks like the norm parameter is causing the error.

But, getting the spectral model from the onoff Observation container seems to work well:

onoff_spectral = phagen.obs().models()[srcname].spectral()
print(onoff_spectral.eval(gammalib.GEnergy(9.2,'TeV'), gammalib.GTime()))

with output:

2.1648981409403132e-29

It seems like its also working at the energies reported in the Error during ctlike running

sergiohcdna commented 3 years ago

Inspection of the source code for GCTAOnOffObservation reveals that probably the problem is during the calculation of number of predicted counts for the On region (N_gamma) that is required for the likelihood method.

Checking more info

sergiohcdna commented 3 years ago

The error duringctlike persists whenc changing statistic used to compute the likelihood. The error comes from calculation of preficted for On region, again in the N_gamma method

sergiohcdna commented 3 years ago

Probably the problem is related to numerical precision in the likelihood used. For OnOff Observations, the likelihood has some extra terms related to the off regions and the background counts. One possibly solution is to generate table model with a reduced ranges of masses with the mass of interest in the middle of the range.

sergiohcdna commented 3 years ago

Probably the problem is related to numerical precision in the likelihood used. For OnOff Observations, the likelihood has some extra terms related to the off regions and the background counts. One possibly solution is to generate table model with a reduced ranges of masses with the mass of interest in the middle of the range.

Seems that this solutions doesn't work. One more test is to try to add scaling for the normalization parameter I also reduced the number of true energy bins to check if this hellps

sergiohcdna commented 3 years ago

No, not working, the error still persists.

I will check if creating tables only for specific channels works

sergiohcdna commented 3 years ago

Problem solved

I removed the parameter Channel from the table-model and only create the table-model for one specific channel (say b) and the likelihood fitting of GCTAOnOffObservation works. Probably this is related to the fact that parametel Channel is a discontinuos variable (?)

Work to do:

Add class to create table-models for specific channels and incorporate in the csdmatter analysis properly

sergiohcdna commented 3 years ago

See Merge