Bug #1801

Energy dispersion does not work for stacked analysis

Added by Mayer Michael over 8 years ago. Updated over 8 years ago.

Status:ClosedStart date:06/22/2016
Priority:HighDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.1.0
Duration:

Description

I have checked the usage of energy dispersion and I always get the same warnings when edisp is turned on:

Parameter "Prefactor" has zero curvature. Fix parameter.
Parameter "Index" has zero curvature. Fix parameter.

As a result, the fit parameter don’t change in the fit.

I added a script and a file such that this can be reproduced (file and script have to be in the same folder).

script.py Magnifier (2.8 KB) Mayer Michael, 06/22/2016 05:54 PM

pnt.dat Magnifier (55 Bytes) Mayer Michael, 06/22/2016 05:54 PM

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen over 8 years ago

  • Assigned To set to Knödlseder Jürgen
  • Target version set to 1.1.0

I will look into that.

#2 Updated by Knödlseder Jürgen over 8 years ago

It turned our that there were several bugs in GCTACubeEdisp that prevented a correct computation of the energy dispersion for stacked analysis. Probably, these bugs have always existed. They comprised:
  • An incorrect computation of migra in GCTACubeEdisp::set and GCTACubeEdisp::fill (minus sign should be a plus sign)
  • Usage of an incorrect edisp function in GCTACubeEdisp::set and GCTACubeEdisp::fill
  • Omission of the set_eng_axis() method in GCTACubeEdisp::read (a node array was not setup correctly)
  • Wrong looping condition in GCTACubeEdisp::compute_ebounds() (prevented computation of maximum energy)
  • Division instead of multiplication by migra in GCTACubeEdisp::compute_ebounds()

All these bugs have been fixed and the script runs successfully, giving the following results:

2016-06-22T21:55:05: +=================================+
2016-06-22T21:55:05: | Maximum likelihood optimisation |
2016-06-22T21:55:05: +=================================+
2016-06-22T21:55:13:  >Iteration   0: -logL=3845.721, Lambda=1.0e-03
2016-06-22T21:55:21:  >Iteration   1: -logL=3844.454, Lambda=1.0e-03, delta=1.267, max(|grad|)=6.198178 [Index:3]
2016-06-22T21:55:30:  >Iteration   2: -logL=3844.451, Lambda=1.0e-04, delta=0.004, max(|grad|)=0.008839 [Index:3]
2016-06-22T21:55:30:  
2016-06-22T21:55:30:  +==================+
2016-06-22T21:55:30:  | Curvature matrix |
2016-06-22T21:55:30:  +==================+
2016-06-22T21:55:30:  === GMatrixSparse ===
2016-06-22T21:55:30:   Number of rows ............: 10
2016-06-22T21:55:30:   Number of columns .........: 10
2016-06-22T21:55:30:   Number of nonzero elements : 16
2016-06-22T21:55:30:   Number of allocated cells .: 528
2016-06-22T21:55:30:   Memory block size .........: 512
2016-06-22T21:55:30:   Sparse matrix fill ........: 0.16
2016-06-22T21:55:30:   Pending element ...........: 0
2016-06-22T21:55:30:   Fill stack size ...........: 0 (none)
2016-06-22T21:55:30:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2016-06-22T21:55:30:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2016-06-22T21:55:30:   0, 0, 71.4554, -729.518, 0, 0, 2.4137, 0.0744688, 0, 0
2016-06-22T21:55:30:   0, 0, -729.518, 9114.55, 0, 0, -17.6776, -4.02596, 0, 0
2016-06-22T21:55:30:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2016-06-22T21:55:30:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2016-06-22T21:55:30:   0, 0, 2.4137, -17.6776, 0, 0, 862.199, 26.5778, 0, 0
2016-06-22T21:55:30:   0, 0, 0.0744688, -4.02596, 0, 0, 26.5778, 218.022, 0, 0
2016-06-22T21:55:30:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2016-06-22T21:55:30:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2016-06-22T21:55:38: 
2016-06-22T21:55:38: +=========================================+
2016-06-22T21:55:38: | Maximum likelihood optimisation results |
2016-06-22T21:55:38: +=========================================+
2016-06-22T21:55:38: === GOptimizerLM ===
2016-06-22T21:55:38:  Optimized function value ..: 3844.451
2016-06-22T21:55:38:  Absolute precision ........: 0.005
2016-06-22T21:55:38:  Acceptable value decrease .: 2
2016-06-22T21:55:38:  Optimization status .......: converged
2016-06-22T21:55:38:  Number of parameters ......: 10
2016-06-22T21:55:38:  Number of free parameters .: 4
2016-06-22T21:55:38:  Number of iterations ......: 2
2016-06-22T21:55:38:  Lambda ....................: 1e-05
2016-06-22T21:55:38:  Maximum log likelihood ....: -3844.451
2016-06-22T21:55:38:  Observed events  (Nobs) ...: 3417.000
2016-06-22T21:55:38:  Predicted events (Npred) ..: 3416.993 (Nobs - Npred = 0.00721418)
2016-06-22T21:55:38: === GModels ===
2016-06-22T21:55:38:  Number of models ..........: 2
2016-06-22T21:55:38:  Number of parameters ......: 10
2016-06-22T21:55:38: === GModelSky ===
2016-06-22T21:55:38:  Name ......................: Crab
2016-06-22T21:55:38:  Instruments ...............: all
2016-06-22T21:55:38:  Instrument scale factors ..: unity
2016-06-22T21:55:38:  Observation identifiers ...: all
2016-06-22T21:55:38:  Model type ................: PointSource
2016-06-22T21:55:38:  Model components ..........: "SkyDirFunction" * "PowerLaw" * "Constant" 
2016-06-22T21:55:38:  Number of parameters ......: 6
2016-06-22T21:55:38:  Number of spatial par's ...: 2
2016-06-22T21:55:38:   RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
2016-06-22T21:55:38:   DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
2016-06-22T21:55:38:  Number of spectral par's ..: 3
2016-06-22T21:55:38:   Prefactor ................: 5.98652e-16 +/- 2.76669e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2016-06-22T21:55:38:   Index ....................: -2.5052 +/- 0.0244963 [-0,-5]  (free,scale=-1,gradient)
2016-06-22T21:55:38:   PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2016-06-22T21:55:38:  Number of temporal par's ..: 1
2016-06-22T21:55:38:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2016-06-22T21:55:38: === GCTAModelCubeBackground ===
2016-06-22T21:55:38:  Name ......................: BackgroundModel
2016-06-22T21:55:38:  Instruments ...............: CTA, HESS, MAGIC, VERITAS
2016-06-22T21:55:38:  Instrument scale factors ..: unity
2016-06-22T21:55:38:  Observation identifiers ...: all
2016-06-22T21:55:38:  Model type ................: "PowerLaw" * "Constant" 
2016-06-22T21:55:38:  Number of parameters ......: 4
2016-06-22T21:55:38:  Number of spectral par's ..: 3
2016-06-22T21:55:38:   Prefactor ................: 0.979871 +/- 0.0341226 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)
2016-06-22T21:55:38:   Index ....................: -0.0655904 +/- 0.0678538 [-5,5]  (free,scale=1,gradient)
2016-06-22T21:55:38:   PivotEnergy ..............: 1e+06 MeV (fixed,scale=1e+06,gradient)
2016-06-22T21:55:38:  Number of temporal par's ..: 1
2016-06-22T21:55:38:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

#3 Updated by Knödlseder Jürgen over 8 years ago

  • Status changed from New to Closed
  • % Done changed from 0 to 100

Merged into devel.

#4 Updated by Mayer Michael over 8 years ago

Thanks for catching this. I guess these bugs are on me :)

Also available in: Atom PDF