Skip to content
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

Fix AREASCAL read for PHA count spectrum #1334

Merged
merged 3 commits into from Mar 2, 2018

Conversation

@registerrier
Copy link
Contributor

@registerrier registerrier commented Mar 1, 2018

Add AREASCAL in the list of columns read from OGIP PHA file and pass it to the PHACountSpectrum constructor.

This should solve the issue https://github.com/gammasky/joint-crab/issues/48 @cdeil and @cosimoNigro

@registerrier registerrier requested review from cdeil and joleroi Mar 1, 2018
@cdeil cdeil added the bug label Mar 1, 2018
@cdeil cdeil added this to the 0.8 milestone Mar 1, 2018
@cdeil
Copy link
Member

@cdeil cdeil commented Mar 1, 2018

This is breaking a few existing tests:
https://travis-ci.org/gammapy/gammapy/jobs/347870845#L3766

@registerrier @joleroi Is raising a KeyError on those files appropriate? Or are they valid cases and there's e.g. a default value that could / should be assumed?

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Mar 2, 2018

Right I missed this test. The issue is that the AREASCAL column does not exist for the spectra stored on GAMMAPY_EXTRA.
There should be a test or a try/except on the presence of columns in the read method

@joleroi
Copy link
Contributor

@joleroi joleroi commented Mar 2, 2018

@registerrier @cosimoNigro @cdeil Thanks a lot for cleaning up my mess 😁

Right I missed this test. The issue is that the AREASCAL column does not exist for the spectra stored on GAMMAPY_EXTRA.

I'd rather update gammapy-extra. The AREASCAL column should always be present and maybe filled with 1. Do you agree @registerrier ? I can do this change.

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Mar 2, 2018

There might be cases where the AREASCAL is not there, e.g. if it has been produced by other means or specific scripts made by some people. In the long run we might want to do MWL fits with X-ray spectra which won't have the column in most of the cases.
I would say it is better to have the possibility to read it and possibly print a warning for more safety.

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Mar 2, 2018

BTW @joleroi , in PHACountSpectrum.to_sherpa areascal is hardcoded at 1.

return DataPHA(
            name=name,
            channel=(table['CHANNEL'].data + 1).astype(SherpaFloat),
            counts=table['COUNTS'].data.astype(SherpaFloat),
            quality=table['QUALITY'].data,
            exposure=self.livetime.to('s').value,
            backscal=backscal,
            areascal=1.,
            syserror=None,
            staterror=None,
            grouping=None,
        )

Is it the expected behavior or should it be:

return DataPHA(
            name=name,
            channel=(table['CHANNEL'].data + 1).astype(SherpaFloat),
            counts=table['COUNTS'].data.astype(SherpaFloat),
            quality=table['QUALITY'].data,
            exposure=self.livetime.to('s').value,
            backscal=backscal,
            areascal=self.areascal,
            syserror=None,
            staterror=None,
            grouping=None,
        )

Thanks.

@joleroi
Copy link
Contributor

@joleroi joleroi commented Mar 2, 2018

BTW, in PHACountSpectrum.to_sherpa areascal is hardcoded at 1. Is it the expected behavior or should it be:

To me, this seems like a 🪲

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Mar 2, 2018

OK.
It seems many tests are broken now. In particular, those testing the fit results.
I assume they were reading from the existing ogip files. How should we proceed?
Replace the existing value in the tests I suppose.

In many cases the relative tolerance of 1e-7 is too strong.
How much would be a reasonable value?
e.g. for spectral indices: is an absolute tolerance of 0.01-0.02 not sufficient?

@joleroi
Copy link
Contributor

@joleroi joleroi commented Mar 2, 2018

@registerrier

The files in gammapy-extra don't have an AREASCAL column, so I'm surprised that the fit results change. Note that this is not due to the tolerance, e.g. here

I am looking in to this now

quality=None
areascal=None
backscal=None
if 'QUALITY' in counts_table:
Copy link
Contributor

@joleroi joleroi Mar 2, 2018

This should be

if 'QUALITY' in counts_table.colnames:

This is also the reasons why the tests fail ...
@registerrier I tried to push the changes to your branch but apparently I don't have access. Either you grant acess to me or you correct this yourself, up to you 😄

@registerrier
Copy link
Contributor Author

@registerrier registerrier commented Mar 2, 2018

Ok thanks for the correction @joleroi. Sorry for that. I have committed the change now all tests run.

@joleroi
Copy link
Contributor

@joleroi joleroi commented Mar 2, 2018

@cdeil please approve, I don't like hitting a red ⏺️

@cdeil cdeil removed their request for review Mar 2, 2018
cdeil
cdeil approved these changes Mar 2, 2018
@cdeil
Copy link
Member

@cdeil cdeil commented Mar 2, 2018

@registerrier @joleroi - Thanks!

@cosimoNigro - Please update to Gammapy master and update the analyses in the joint Crab repo.

@lmohrmann - If you use OGIP spec I/O from Gammapy for the HESS checks, you might want to update to latest Gammapy and re-run.

@cdeil cdeil merged commit b31445f into gammapy:master Mar 2, 2018
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@cdeil
Copy link
Member

@cdeil cdeil commented Mar 2, 2018

I re-ran the spectral analyses in https://github.com/gammasky/HESS-DL3-DR1/commit/b8c26ab3c292e7879b344120789c7bcadaeeac40 via

time ./make.py analysis-spec-gammapy --extract --fit

The fit results are identical.

There I do store OGIP files, read them back and then do the fit.
@joleroi - do you understand why the analysis results there aren't affected?

@cdeil cdeil changed the title Add AREASCAL in the list of columns read from OGIP PHA file and pass … Fix AREASCAL read for PHA count spectrum Apr 11, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants