Bug #1150
Spectral model log parabola not compatible between Fermi science tools and gammalib/ctools
| Status: | Closed | Start date: | 02/21/2014 | |
|---|---|---|---|---|
| Priority: | Normal | Due date: | ||
| Assigned To: | - | % Done: | 100% | |
| Category: | - | |||
| Target version: | - | |||
| Duration: |
Description
Performing further cross checks between the fermi science tools and ctools the following problem came up:
If one takes the xmlfile from the science tools and tries to analyse it with gammalib/ctools the following happens:
input example:
<spectrum normPar="norm" type="LogParabola">
<parameter error="0.0403546564" free="1" max="100000" min="1e-05" name="norm" scale="6.664498894e-12" value="0.330797703"/>
<parameter error="0.1932913856" free="1" max="5" min="-0" name="alpha" scale="1" value="2.762437733"/>
<parameter error="0.9494770106" free="1" max="10" min="-10" name="beta" scale="0.1" value="3.304939119"/>
<parameter free="0" max="300" min="0.02" name="Eb" scale="1000" value="1.444281016"/>
</spectrum>
this means if the source is fitted also the Pivotenergy is free, which leads to the huge errors on the fitted parameters ( I think due to problem of correlated parameters) an example shown here:
<spectrum type="LogParabola">
<parameter name="Prefactor" value="0.001" error="0.0510866" scale="1e-09" min="0.001" max="1000" free="1" />
<parameter name="Index" value="-3.71444e-05" error="878492" scale="1" min="-10" max="10" free="1" />
<parameter name="Curvature" value="-0.319421" error="0.108217" scale="1" min="-10" max="10" free="1" />
<parameter name="Scale" value="511.743" error="7.03715e+08" scale="1" min="20" max="10000" free="1" />
</spectrum>
But also fixing the Pivotenergy does not solve the problem at the moment since the definitions of the parameters in the fermi science tools and gammalib/ctools are not completely compatible. This leads to a fit result like this:
<spectrum type="LogParabola">
<parameter name="Prefactor" value="1.13197" error="0.0615543" scale="6.6645e-12" min="1e-05" max="100000" free="1" />
<parameter name="Index" value="-4.18061" error="0.0833878" scale="1" min="-5" max="5" free="1" />
<parameter name="Curvature" value="-10" error="0" scale="0.1" min="-10" max="10" free="1" />
<parameter name="Scale" value="1.44428" scale="1000" min="0.02" max="300" free="0" />
</spectrum>
Recurrence
No recurrence.
History
#1
Updated by Schulz Anneli over 10 years ago
- Description updated (diff)
#2
Updated by Schulz Anneli over 10 years ago
- Description updated (diff)
#3
Updated by Knödlseder Jürgen over 10 years ago
Are the parameters defined identically? I see that for Fermi
norm = 0.330797703 alpha = 2.762437733 beta = 3.304939119
and for GammaLib
Prefactor = 1.13197 Index = -4.18061 Curvature = -10
In particular, the curvature runs against the minimum, so I’m wondering whether the formulae in ScienceTools and GammaLib are at the end the same?
#4
Updated by Mayer Michael over 10 years ago
in the Fermi Science Tools, the parameters are defined differently compared to gammalib. In order to be consistent within gammalib, we decided to to define Index and Curvature with a different sign than alpha and beta (wikipage). Maybe there got something mixed up when reading the signs.
However, since the scale of the Curvature is 0.1 in the example file, the boundaries of [-10;10] might be too narrow.
#5
Updated by Schulz Anneli over 10 years ago
Hi,
so I spent some more time on it and found out that it’s not the log parabola which causes the issue but the super exponential cutoff power law. The number of predicted events goes of by factors 1e6 or 1e7...
I went back to only diffuse and one source and checked the values:
2014-03-04T19:04:55: === GOptimizerLM === 2014-03-04T19:04:55: Optimized function value ..: -386110.673 2014-03-04T19:04:55: Absolute precision ........: 1e-06 2014-03-04T19:04:55: Optimization status .......: converged 2014-03-04T19:04:55: Number of parameters ......: 14 2014-03-04T19:04:55: Number of free parameters .: 5 2014-03-04T19:04:55: Number of iterations ......: 41 2014-03-04T19:04:55: Lambda ....................: 0.0001 2014-03-04T19:04:55: Maximum log likelihood ....: 386110.673 2014-03-04T19:04:55: Observed events (Nobs) ...: 1057409.000 2014-03-04T19:04:55: Predicted events (Npred) ..: 1057409.000 (Nobs - Npred = -7.26331e-05)
If the parameters are fixed to the 2FGL values (e.g. here for
2FGL J2017.3+0603 Number of spectral par's ..: 5 2014-03-04T19:12:10: Prefactor ................: 4.31505e-12 [1e-13,1e-09] ph/cm2/s/MeV (fixed,scale=1e-11,gradient) 2014-03-04T19:12:10: Index1 ...................: -1.04592 [-0,-5] (fixed,scale=-1,gradient) 2014-03-04T19:12:10: Cutoff ...................: 3743.96 [100,1e+06] MeV (fixed,scale=10000,gradient) 2014-03-04T19:12:10: PivotEnergy ..............: 1525.24 [896.688,1525.24] MeV (fixed,scale=1,gradient) 2014-03-04T19:12:10: Index2 ...................: 1 [0,5] (fixed,scale=1,gradient)
I get the following result:
2014-03-04T19:12:10: === GOptimizerLM === 2014-03-04T19:12:10: Optimized function value ..: 1448179.296 2014-03-04T19:12:10: Absolute precision ........: 1e-06 2014-03-04T19:12:10: Optimization status .......: converged 2014-03-04T19:12:10: Number of parameters ......: 22 2014-03-04T19:12:10: Number of free parameters .: 5 2014-03-04T19:12:10: Number of iterations ......: 41 2014-03-04T19:12:10: Lambda ....................: 1000 2014-03-04T19:12:10: Maximum log likelihood ....: -1448179.296 2014-03-04T19:12:10: Observed events (Nobs) ...: 1057409.000 2014-03-04T19:12:10: Predicted events (Npred) ..: 3042652.294 (Nobs - Npred = -1.98524e+06)
if I free the prefactor for the source it is fit to the minimum (1e-13).
So there is something wrong when reading in xmlfiles containing plsuperexpcutoff from fermi.
#6
Updated by Schulz Anneli over 10 years ago
Since the problem does not seem to be related to the log parabola, I put more results on comparisons in the wiki:
https://cta-redmine.irap.omp.eu/projects/gammalib/wiki/Validation_of_Pass_7REP
#7
Updated by Knödlseder Jürgen over 10 years ago
I was doing a check by fitting the Crab using a log-parabola spectrum. The analysis was done for one week of Crab data, using 20 bins in energy between 200 MeV and 20 GeV.
Using the ScienceTools I get
Crab: norm: 2.32745 +/- 1.27195 alpha: 1.66168 +/- 0.182328 beta: 0.194544 +/- 0.0528259 Eb: 299.806 +/- 97.7888 TS value: 1435.31 Flux: 1.03049e-06 +/- 5.70869e-08 photons/cm^2/s Extragal_diffuse: Normalization: 1.49686 +/- 0.354927 Flux: 0.000104134 +/- 2.46826e-05 photons/cm^2/s Galactic_diffuse: Value: 1.10188 +/- 0.0397298 Flux: 0.000328965 +/- 1.18608e-05 photons/cm^2/s WARNING: Fit may be bad in range [3169.79, 3990.52] (MeV) Total number of observed counts: 4832 Total number of model events: 4831.5 -log(Likelihood): 25972.37267 Elapsed CPU time: 45.78
and using the ctlike I get
2014-03-11T17:40:18: Maximum log likelihood ....: -25973.120 2014-03-11T17:40:18: Observed events (Nobs) ...: 4832.000 2014-03-11T17:40:18: Predicted events (Npred) ..: 4832.000 (Nobs - Npred = 4.23437e-05) 2014-03-11T17:40:18: Prefactor ................: 2.31551e-09 +/- 0.00184344 [1e-16,1e-06] ph/cm2/s/MeV (free,scale=1e-09,gradient) 2014-03-11T17:40:18: Index ....................: -1.66673 +/- 181170 [-5,5] (free,scale=1,gradient) 2014-03-11T17:40:18: Curvature ................: -0.189644 +/- 0.0559916 [-10,10] (free,scale=1,gradient) 2014-03-11T17:40:18: PivotEnergy ..............: 299.557 +/- 1.43086e+08 [10,1e+06] MeV (free,scale=1,gradient)
This looks pretty consistent, at least when the values are concerned (note that alpha is -Index and beta is -Curvature). The errors are not corrected due to correlation of parameters, which is not handled correctly so far (known problem).
So is there still a bug related to the log-parabola spectrum or can this be closed?
#8
Updated by Knödlseder Jürgen over 10 years ago
- Status changed from New to In Progress
#9
Updated by Schulz Anneli over 10 years ago
You are right. As written above, the problem is not about the log parabola, but something else. So this issue can be closed, further discussion on the science tools validation will be in the wiki. Sorry that I blamed a wrong part of the code, it just came up there.
#10
Updated by Knödlseder Jürgen over 10 years ago
- Status changed from In Progress to Closed
- % Done changed from 0 to 100