GammaLib

We have an existing C++ implementation (using GSL via ROOT for heavy numerical lifting) of XSPEC-style forward folding with the energy resolution matrix in HESS.
The plan is to have this open-sourced on Github before the coding sprint, but in any case we can probably re-use some parts for Gammalib.

For testing we can compare to XSPEC.

Jürgen, can you briefly comment on the status in gammalib and possible sub-tasks?

Is there an RMF reader? Which classes would we work on? ...

Deil Christoph wrote:

Jürgen, can you briefly comment on the status in gammalib and possible sub-tasks?

Is there an RMF reader? Which classes would we work on? ...

Nothing has really be done for energy resolution, but the structure is in principle there to plug-in a probability density function taking true and reconstructed energies as parameters. So maybe a starting point would be to implement a very simple energy dispersion function (e.g. a Gaussian), and see if we can get the code working with such a function. The code would be probably awefully slow in a first run, but can then be optimised once we get the functionality working

Work is needed here on the response side, but also on the MC side, so that we can do some testing.

Concerning RMFs, there is the GRmf class available, and I recently added support for variable-length columns to the FITS module so that I can read the files that are currently produced by HESS. My plan is to introduce a generic XSpec observation interface, so that any data in XSpec format can be used in any GammaLib analysis. This should also support energy redistribution, which would be easier as it’s just a multiplication of matrix elements to a data vector.

• implement a class derived from GCTAEdisp (for example a simple Gaussian), to be used for unbinned and binned analysis
• implement full XSpec support, using the GRmf class for Xspec analysis

• GModelSky::integrate_energy: implement integration over true energy
• GCTAResponse: add a method that returns the true energy integration boundaries for a given measured energy; this probably needs a similar method to be implement for GCTAEdisp.
• GCTAResponse::nedisp: implement method (integration of energy dispersion over [emin,emax] of dataspace)
• GCTAEdispPerfTable: implement Gaussian energy dispersion based on rms values in performance table
• GCTAResponse::load_edisp: implement method that loads energy dispersion information (similar to load_psf()), and call this method from GCTAResponse::load (so far, GCTAResponse::load_edisp is defined empty in GCTAResponse.hpp)
• Add MC sampling to GCTAResponse:mc

I hope that’s it ...

Jürgen, thanks for making this list, very helpful!

We also want energy resolution in binned analysis.
Can you please also make a list of things needed for that?

Deil Christoph wrote:

Jürgen, thanks for making this list, very helpful!

We also want energy resolution in binned analysis.
Can you please also make a list of things needed for that?

The changes I listed should in fact do the job for unbinned and binned analysis. But for binned you probably would like to precompute stuff beforehand so that computation is faster. Hence some additional thinking is required here. This is probably linked to issue #803 which asks for implementing a interface that is also optimized for binned analysis (so far the interfaces are only optimized for unbinned).

However, if you’d like to use the RMFs, as defined during the last coding sprint, additional work is needed. But we even have not a fitter yet for Xspec like analysis, hence when we implement the fitter we probably could do the RMFs in one shot.

Just one thing: when implementing energy dispersion we should always have the possibility to switch it off, mainly for computational speed reasons, but also for comparison reasons. Shouldn’t be complicated to implement, by maybe needs a flag that can be toggled by a specific GCTAResponse method.

I think we’d like to use RMFs for energy resolution, at least to start with.
RMFs are already available in some of the HESS exporters and we haven’t talked about another data format yet.

Isn’t the use of RMFs to store energy resolution info independent of the question whether we do an XSPEC-style analysis (i.e. for a source region) or a cube-style likelihood analysis? For now we just want to get correct spectra for small sources in the FOV, so I think it’s OK to apply the correct energy resolution from that position to all positions in the FOV, no?

Deil Christoph wrote:

I think we’d like to use RMFs for energy resolution, at least to start with.
RMFs are already available in some of the HESS exporters and we haven’t talked about another data format yet.

Isn’t the use of RMFs to store energy resolution info independent of the question whether we do an XSPEC-style analysis (i.e. for a source region) or a cube-style likelihood analysis? For now we just want to get correct spectra for small sources in the FOV, so I think it’s OK to apply the correct energy resolution from that position to all positions in the FOV, no?

I think we should implement the full framework, and then we could add a GCTAEdispRmf method which is equivalent to the GCTAAeffArf method. But in fact, this is not really the issue. The issue is that for the standard binned/unbinned analysis, you have to implement the set of methods described above. For an Xspec-style analysis, you won’t need any of these methods, as energy dispersion then boils done to a mere vector multiplication. And this one would implement at the OnOff fitter level I guess.

That’s our working branch for this feature: 1036-implement-energy-resolution

Jürgen, could you please review (and merge if OK) the following branches into the issue_1036 integration branch?
https://github.com/cdeil/gammalib/compare/action_1109
https://github.com/cdeil/ctools/compare/action_1109

Currently the ctools test fail indicating what to work on next:

*****************
* Test cscripts *
*****************
Test cspull: 
Traceback (most recent call last):
  File "/Users/deil/code/ctools/scripts/cspull.py", line 360, in <module>
    app.execute()
  File "/Users/deil/code/ctools/scripts/cspull.py", line 173, in execute
    self.run()
  File "/Users/deil/code/ctools/scripts/cspull.py", line 210, in run
    result = self.trial(seed)
  File "/Users/deil/code/ctools/scripts/cspull.py", line 286, in trial
    like = obsutils.fit(obs, log=self.m_log, debug=self.m_debug)
  File "/Users/deil/code/ctools/scripts/obsutils.py", line 147, in fit
    like.run()
  File "/Users/deil/code/ctools/pyext/ctools.py", line 319, in run
    return _ctools.ctlike_run(self)
RuntimeError: *** ERROR in GModelSky::integrate_energy(GEvent&, GTime&, GObservation&, bool): Feature not implemented. In case that you need this feature for your application please submit a feature request on https://sourceforge.net/projects/gammalib/, join this error message and provide a detailed description of your needs.
.
 pull.dat file is not valid

Both Ellis and I won’t have much time to work on GammaLib / ctools in the next month, so we’d like to wrap this up for now.

Jürgen, can you merge the existing code into devel (I don’t know if this is already done) and make separate issues for future work so that we can come back to this ~ in April?

I think some things that still need to be done:
- Add option to ctlike to apply energy resolution in the fit (I think I did that already and if not it should be easy to do).
- More unit tests? (where can se see the coverage of the code we wrote? I doubt it’s well-tested at the moment.)
- Some numerical precision / computation speed tests?
- Make high-level tests (pull distributions) simulating / fitting spectra with energy dispersion.
- Write high-level documentation an energy dispersion (including maybe a plot or two) in the user guide.

Deil Christoph wrote:

Both Ellis and I won’t have much time to work on GammaLib / ctools in the next month, so we’d like to wrap this up for now.

Jürgen, can you merge the existing code into devel (I don’t know if this is already done) and make separate issues for future work so that we can come back to this ~ in April?

Most of the code is aleady in devel, I just started to look in the new code using the RMFs. Some work is needed here (see #1120).

I think some things that still need to be done:
- Add option to ctlike to apply energy resolution in the fit (I think I did that already and if not it should be easy to do).

It’s there aleady (parameter edisp). I made some first tests and everything seems to work. But we need science validation now!

- More unit tests? (where can se see the coverage of the code we wrote? I doubt it’s well-tested at the moment.)

• GCTAEdisp: 61.3%
• GCTAEdispPerfTable: 62.8%

Is not so bad, but can of course be improved.

- Some numerical precision / computation speed tests?

Indeed, needed.

- Make high-level tests (pull distributions) simulating / fitting spectra with energy dispersion.

Also. I can start doing this at some point.

- Write high-level documentation an energy dispersion (including maybe a plot or two) in the user guide.

Agree.

Knödlseder Jürgen wrote:

Deil Christoph wrote:

- Add option to ctlike to apply energy resolution in the fit (I think I did that already and if not it should be easy to do).

It’s there aleady (parameter edisp). I made some first tests and everything seems to work. But we need science validation now!

I tried a simple example and energy dispersion doesn’t seem to work with ctlike in devel at the moment!?
https://github.com/cdeil/gammalib-tutorials/tree/master/energy_dispersion

Maybe I’m making a simple mistake?
Can you share the test case where it works for you?

Jürgen, do you have a working example of a ctlike fit with energy resolution?

Thanks for making the check.

No, the Gaussian model has no bias as the energy dispersion works well for performance tables. I suspect that the problem is that there is a lot of noise in the actual energy dispersion matrix and that the integration method is not accurate enough. I’m working on that.

Also note that the Gaussian model is much more sensitive to energy dispersion problems than the power law model is, so it may just be that you cannot detect the effects.

Mayer Michael wrote:

I understand. I therefore also tried to use a Gaussian function as spectral model. The used model and the pull distributions are attached. The pull distribution don’t look too bad. I assume I do not use enough statistics to make these effects visible.
The noise in the computation might also depend on the parameter values, right? What values were you using?

I use a line at 5 TeV with a sigma of 1 TeV (see above).

I was using a line at 2 TeV with 0.4 TeV width. So mine would be more pronounced. Could that influence the pull distributions?

Maybe. What response function are you using (I still suspect that the noise in the response function is influencing the results)?

I am using 'prod2’ and 'South_50h’. Of course I can try to do the same exercise using the HESS IRFs in order to see if there is a difference.

Mayer Michael wrote:

I am using 'prod2’ and 'South_50h’. Of course I can try to do the same exercise using the HESS IRFs in order to see if there is a difference.

I confirm that using the 'prod2’ and 'South_50h’ I get much better results. I will post pull distributions once they have a sufficient number of samples.