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 Gauss PSF to_table_psf method #523

Merged
merged 2 commits into from May 2, 2016

Conversation

Projects
None yet
3 participants
@adonath
Member

adonath commented Apr 26, 2016

This PR adds a .to_table_psf() method to the EnergyDependentMultiGaussPSF class. The containment radii are tested against the analytical solution, achieving a precision of 3% for 0.5, 1, 10, 30 TeV for the chosen test dataset.

Discussing with @cdeil I realized today that the UnivariateSpline only caused problems for the computation of containment radii. I still included a small work around (aka "hack") to stabilize the spline interpolation for this case. So I could use the containment radii for the testing as well.

Here is notebook with some further minimal checks and an example how to compute the PSF in an energy band:

https://github.com/gammapy/gammapy-extra/blob/master/experiments/three_gauss_to_table_psf.ipynb

@leajouvin @cdeil

@cdeil cdeil added the feature label Apr 27, 2016

@cdeil cdeil added this to the 0.5 milestone Apr 27, 2016

@cdeil cdeil self-assigned this Apr 27, 2016

def to_table_psf(self, theta=None, offset_max=None, binsz=0.01, exposure=None):

This comment has been minimized.

@cdeil

cdeil Apr 27, 2016

Member

@adonath - From the parameter description it's not clear what theta, offset_max and binsz are ... at least I had to read all the code to understand how you use them.

Here's my suggestions:

  • Put this as the description of theta: "Offset in the field of view"
  • Remove the offset_max and binsz parameter. Instead put one offset parameter, which is the offset grid at which the PSF is sampled and that is passed to EnergyDependentTablePSF. You could make an example showing people how to call gammapy.utils.axis.sqrt_space before to create a good offset array. I'm not sure if this should have a default ... you could put 1 deg which would work for all IACTs as far as I know.
  • Add an energy parameter with default None and energy grid taken from self, as you do now. But if the user wants to sample with some other energy grid, then they have the option.

There are alternatives like making this a factory method on the EnergyDependentTablePSF class, or passing a EnergyDependentTablePSF object with pre-defined grid into this method, that have advantages when it comes to duplicating axis grid specs (like the log energy grid and sqrt offset grid). But they require more changes, also to the EnergyDependentTablePSF constructor, and aren't obviously better to me. So my suggestion would be to stick with what you have here, just to improve the parameters a bit like desribed above.

energies = ebounds.log_centers
# Default for exposure
exposure = exposure or Quantity(np.ones(len(energies)), 'cm^2 s')

This comment has been minimized.

@cdeil

cdeil Apr 27, 2016

Member

It would probably be better to move this default for exposure to the EnergyDependentTablePSF class.
It's in EnergyDependentTablePSF because it's there for Fermi.
But for IACTs we probably won't always have it / use it.

So since this weighting is something that's done in EnergyDependentTablePSF, not here, I'd prefer to also have the default there, and then here you can just take exposure=None in and pass it along to the EnergyDependentTablePSF initialiser in case someone does want to set it.

table_psf_at_energy = table_psf.table_psf_at_energy(energy)
psf_at_energy = psf.psf_at_energy_and_theta(energy, theta)
# Plot containment radius vs. containment

This comment has been minimized.

@cdeil

cdeil Apr 27, 2016

Member

Remove this comment? (You're not plotting anything in this test, probably this is a copy & paste remnant)

@requires_data('gammapy-extra')
@pytest.mark.parametrize(('energy'), ENERGIES)
def test_to_table_psf(energy):
filename = gammapy_extra.filename('test_datasets/unbundled/irfs/psf.fits')

This comment has been minimized.

@cdeil

cdeil Apr 27, 2016

Member

The files in test_datasets/unbundled are old and we should get rid of them.
Please use datasets/hess-crab4-hd-hap-prod2/run023400-023599/run023523/hess_psf_3gauss_023523.fits.gz here which matches the files we produce now.
Probably this cleanup will only happen when the HESS public data release is through A&R, i.e files are final.
But still ... let's test the latest files we have for now.

for i, energy in enumerate(energies):
psf_gauss = self.psf_at_energy_and_theta(energy, theta)
psf_value[i] = Quantity(psf_gauss(offsets, np.zeros_like(offsets)), 'deg^-2')

This comment has been minimized.

@JouvinLea

JouvinLea Apr 27, 2016

Contributor

@cdeil @adonath @registerrier
For the moment in the method self.psf_at_energy_and_theta() there is no interpolation. You select the energy bin and offset bin that is the closest from the energy and theta that you gave. You select the parameters of the PSF like sigma and norm for the tripple gauss and then you create an EnergyDependentTablePSF with psf_gauss(offsets, np.zeros_like(offsets)). But in theta, you could have quite large variation in the PSF parametrisation so an interpolation would be better no?
For the interpolation, it could be very tricky to do it directly on the parameters of the triple gauss for example and I think even completely false. Maybe we could define 4 EnergyDependentTablePSF for the two energies and two theta closest to the energy and theta you give to self.psf_at_energy_and_theta. And then create an interpolator in (energy,theta) on these tables for each offset of the table? With this method, the interpolation will not depend on the profil you assumed for the PSF (tripple gauss, king ...etc)
I'm not sure to be really clear, this is difficult by writing...

This comment has been minimized.

@cdeil

cdeil Apr 27, 2016

Member

What Gammalib is doing is interpolating the PSF parameters.
We could do the same.
But I'm not sure this is a good, maintainable way to go, so I would suggest to leave what we have now (you could call it no or nearest neighbor interpolation) for analytical IRFs.

Generally I'm not convinced we need analytical IRFs at DL3 at all.
I think IRF3 producers should just use good binnings and always histogram their IRFs before handing them out (if they had analytical IRFs internally), and then it's really simple for science tools.

So in Gammapy for now we do want to support analytical PSFs, because people use them, but my suggestion to always convert to histogram PSF immediately and then to only develop and test that class well.

So @JouvinLea - We need the equivalent to_table_psf for the King PSF in Gammapy so that we can analyse Model++ data. Would you be willing to implement it in the coming days (after the review on this one is done and you have a good template how it should be implemented)?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 27, 2016

Contributor

Ok so let's forget the interpolation for now!

@cdeil
ok for implementing the king profile.
I will have to develop the same kind of class as the ones in gauss.py in morphology ?
I wait for this PR to be finish. If you have some inputs before you leave in holiday of how I should implement this tell me :-)

@adonath adonath force-pushed the adonath:three_gauss_to_table_psf branch from 55d435b to d75fa96 May 2, 2016

@adonath adonath force-pushed the adonath:three_gauss_to_table_psf branch from d75fa96 to 4f42c46 May 2, 2016

@adonath adonath merged commit 416dfcd into gammapy:master May 2, 2016

2 checks passed

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

@cdeil cdeil changed the title from Add .to_table_psf() method to EnergyDependentMultiGaussPSF to Add Gauss PSF to_table_psf method May 19, 2016

@adonath adonath deleted the adonath:three_gauss_to_table_psf branch Nov 20, 2018

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