Bug #1150

Spectral model log parabola not compatible between Fermi science tools and gammalib/ctools

Added by Schulz Anneli about 10 years ago. Updated almost 10 years ago.

Status:ClosedStart date:02/21/2014
Priority:NormalDue 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 about 10 years ago

  • Description updated (diff)

#2 Updated by Schulz Anneli about 10 years ago

  • Description updated (diff)

#3 Updated by Knödlseder Jürgen about 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 about 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 about 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 about 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 about 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 about 10 years ago

  • Status changed from New to In Progress

#9 Updated by Schulz Anneli about 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 almost 10 years ago

  • Status changed from In Progress to Closed
  • % Done changed from 0 to 100

Also available in: Atom PDF