New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add offset-dependent table PSF class #556

Merged
merged 10 commits into from Jun 28, 2016

Conversation

Projects
None yet
3 participants
@dltiziani
Contributor

dltiziani commented Jun 7, 2016

  • EnergyDependentTablePSF can be created from psf_table format fits file at a given offset.
  • Added a script for quick PSF containment plots for given runs

dltiziani added some commits Jun 3, 2016

Added psf_table fits format
Signed-off-by: Domenico Tiziani <domenico.tiziani@fau.de>

@cdeil cdeil added the feature label Jun 8, 2016

@cdeil cdeil added this to the 0.5 milestone Jun 8, 2016

@cdeil cdeil self-assigned this Jun 8, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 8, 2016

@dltiziani - Thanks for the pull request!

It's great that you're pushing forward with the HESS PSF investigations and this table PSF format specifically.

So far our strategy has been to add one class per IRF format that does I/O and other methods:
http://docs.gammapy.org/en/latest/irf/index.html#classes

I would like to stick with that pattern also for this PSF format you're adding in this PR:
http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/psf/psf_table/index.html

So my suggestion is that you introduce a new class. I don't know a good name at the moment, maybe just call it PSF3D for now (because it has three axes: RAD, THETA, ENERGY) and then we'll do a cleanup of IRF names in the spec and Gammapy IRF class names soon.

A good example to look at is PSFKing for how the .read is implemented and then there's a .to_table_psf to slice out an EnergyDependentTablePSF for a given FOV offset.

It's also important to have tests that execute the functionality once and assert that it's basically working,
e.g. a test like the one here that reads a King PSF and makes a TablePSF and asserts that the content is sensible.
As a test file for this psf_table format, this one should be used:
https://github.com/gammapy/gammapy-extra/blob/master/datasets/hess-crab4-hd-hap-prod2/run023400-023599/run023523/hess_psf_table_023523.fits.gz

@dltiziani - I realise that this is quite a bit of work, especially when you're new to Gammapy development. Please let me know if you want to and have time to do this soon, or if I or some other Gammapy dev should take over this PR and continue with it.

@dltiziani

This comment has been minimized.

Contributor

dltiziani commented Jun 8, 2016

If it's ok for you, I would really like to implement this new class. I think I can learn quite much from that and I also have some time at the moment.

@cdeil cdeil changed the title from Table PSF reading and containment debug plots to Add offset-dependent table PSF class Jun 13, 2016

@joleroi

This comment has been minimized.

Contributor

joleroi commented Jun 23, 2016

@cdeil
Again, a large part of the functionality implemented in this PR is already part of
http://docs.gammapy.org/en/latest/api/gammapy.utils.nddata.NDDataArray.html
(__init__, evaluate, info, energy_logcenters, rad_logcenters, etc.)

I don't understand why you don't encourage people to use the NDDataArray class. At the moment every class in gammapy.irf repeats the same stuff (interpolation, axis handling). This does not only duplicate a lot of code but is also error-prown.

cc @adonath

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 23, 2016

@joleroi - This PR was started before NDDataArray existed. I didn't even have time to have a look.

I'm not sure if it makes sense to ask @dltiziani to rewrite this class using NDDataArray ... it's his first pull request. I'll have a look at this tonight or tomorrow and comment.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 23, 2016

@dltiziani - Thank you for implementing this!

As @joleroi mentioned, he recently introduced a base class NDDataArray base class, which should be able to handle I/O and evaluation for any IRF. Currently it's only used for one IRF, EffectiveAreaTable2D, but the idea is that we use it for the other IRFs as well to avoid duplicated similar boilerplate code.

My suggestion for this PR would be that you add one or two basic tests to this PR first to establish and show that it works:
https://github.com/gammapy/gammapy/blob/master/gammapy/irf/tests/test_psf_king.py#L70

You could also add a script examples/test_psf3.py that's similar to
https://github.com/gammapy/gammapy/blob/master/examples/test_psf_king.py
i.e. contains examples I can quickly run to play with your new class.

Eventually PSF3D should most likely be rewritten as an NDDataArray sub-class.
Whether you want to do this in this PR, or one of us does it in a follow-up PR is up to you, @dltiziani .

ax.set_xlabel('Energy (TeV)')
ax.set_ylabel('Containment radius (deg)')
def plot_psf_vs_rad(self, filename=None, theta=Angle(0, 'deg'), energy=Quantity(1, 'TeV')):

This comment has been minimized.

@cdeil

cdeil Jun 23, 2016

Member

Is it useful to have plot_psf_vs_rad on the PSF3D class?
Can't you use or improve the EnergyDependentTablePSF.plot_psf_vs_theta method to make these plots?

plt.ylabel('PSF (1e-6 sr^-1)')
plt.tight_layout()
if filename != None:

This comment has been minimized.

@cdeil

cdeil Jun 23, 2016

Member

Don't put this option to save the plot here.
We don't do this for any other plotting methods in Gammapy.
Just leave it up to the caller.

class PSF3D(object):
"""Table PSF.

This comment has been minimized.

@cdeil

cdeil Jun 23, 2016

Member

Can you add a link to
http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/psf/psf_table/index.html
to this docstring?

Here's an example how to create a Sphinx link to another Sphinx project:

"""Evaluate formula from here: :ref:`gadf:psf_king`.

In your case the label is this:
https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/blame/master/source/irfs/psf/psf_table/index.rst#L3

Putting this in your docstring should work:

This class implements the format described here: :ref:`gadf:psf_table`.

You can check by running this locally:

python setup.py build_sphinx
# get a coffee
open docs/_build/html/index.html
# Browse to docs for your class
@@ -0,0 +1,18 @@
import numpy as np

This comment has been minimized.

@cdeil

cdeil Jun 27, 2016

Member

For examples/test_psf3d.py:

  • remove numpy import
  • add one line docstring at the top what this is
  • Put print(psf.info()) also before doing the plots.
@@ -753,3 +759,414 @@ def _get_1d_table_psf(self, energy_index, **kwargs):
self._table_psf_cache[energy_index] = table_psf
return self._table_psf_cache[energy_index]
class PSF3D(object):
"""This class implements the format described here: :ref:`gadf:psf_table`.

This comment has been minimized.

@cdeil

cdeil Jun 27, 2016

Member

This intersphinx reference doesn't resolve and causes the docs build to fail:
https://travis-ci.org/gammapy/gammapy/jobs/139788155#L2134

I think this will work:

:ref:`gadf:psf-table`

(You have an underscore, not a minus:)

:ref:`gadf:psf_table`
PSF (3-dim with axes: psf[rad_index, offset_index, energy_index]
energy_thresh_lo : `~astropy.units.Quantity`
Lower energy threshold. Default energy_thresh_lo = 100 GeV
energy_thqresh_hi : `~astropy.units.Quantity`

This comment has been minimized.

@cdeil

cdeil Jun 27, 2016

Member

Typo: energy_thqresh_hi -> energy_thresh_hi

# Set up and compute data
containment = self.containment_radius(self.energy_logcenter(), self.offset, fraction)
print(containment)

This comment has been minimized.

@cdeil

cdeil Jun 27, 2016

Member

Remove print statement.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 27, 2016

@dltiziani - It would be great to have this in master soon, so that we can start using and debugging this PSF format more widely.

The main thing before it can go in would be to fix the travis-ci issue with the Sphinx reference.
And it would be nice to have some more tests that just establish the current behaviour, i.e. to call evaluate and compute some containment radius and to assert on the current values.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 28, 2016

@dltiziani - I'm merging this now. Thank you!

@cdeil cdeil merged commit d0c375c into gammapy:master Jun 28, 2016

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Jun 28, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment