Bug #1848

Modelcube from Table-PSF

Added by Tiziani Domenico over 6 years ago. Updated over 6 years ago.

Status:ClosedStart date:08/29/2016
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.2.0
Duration:

Description

When generating a model cube for an observation that uses a GCTAPsfTable as IRF,
there seems to be a problem with the calculation of expected counts.
For the model of a point source, the counts at the position of the
source are zero and reach a maximum at some distance.
An exemplary model map looks like this:

Using a gaussian PSF parametrisation for the same PSF-data leads to a more sensible map:

map_table.png (66.3 KB) Tiziani Domenico, 08/29/2016 10:21 AM

map_gauss.png (55.5 KB) Tiziani Domenico, 08/29/2016 10:31 AM

psf_table_test.fits.gz (259 KB) Tiziani Domenico, 08/30/2016 12:29 PM

psf_table_original.png (84.1 KB) Knödlseder Jürgen, 08/30/2016 01:24 PM

psf_table_modified.png (89.9 KB) Knödlseder Jürgen, 08/30/2016 01:24 PM

Map_table Map_gauss Psf_table_original Psf_table_modified

Recurrence

No recurrence.

History

#1 Updated by Mayer Michael over 6 years ago

Hi Domenico, thanks for sharing this.
Can you reproduce this problem also with the example table PSF in the gammalib test directory:
- $GAMMALIB/inst/cta/test/caldb/psf_table.fits.gz?

This way one could narrow down if it is a problem in the PSF code or a problem with the PSF file you are using.

#2 Updated by Tiziani Domenico over 6 years ago

OK, with the test-PSF-file in the gammalib repository, the problem does not occur.

The differences to the PSF that I am using is that my file uses the axis ordering described here:
http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/psf/psf_table/index.html,
so ENERGY/THETA/RAD instead of RAD/THETA/ENERGY and that the PSF value is stored in units of
1/sr instead of 1/(deg^2).

#3 Updated by Mayer Michael over 6 years ago

Ok great, this might suggest a problem with your file, or gammalib is interpreting the table a bit differently than you expect. From the ring-like shape on the model maps, I could imagine that the problem could arise from a missing solidangle correction factor per bin?

#4 Updated by Tiziani Domenico over 6 years ago

Mayer Michael wrote:

From the ring-like shape on the model maps, I could imagine that the problem could arise from a missing solidangle correction factor per bin?

I would think so too, but I don’t see why it would then work for the “old” file since the conversion from 1/(deg^2) to 1/sr is just a constant factor.
So I guess that there is a problem with the different axis order, although I have not been able to locate an issue in the gammalib code yet.

#5 Updated by Mayer Michael over 6 years ago

I would think so too, but I don’t see why it would then work for the “old” file since the conversion from 1/(deg^2) to 1/sr is just a constant factor.

Good point, maybe there is an issue with the increasing solid angle per “RAD” bin? Jürgen should have an idea how the exact requirement in gammalib is.

So I guess that there is a problem with the different axis order, although I have not been able to locate an issue in the gammalib code yet.

The axes are read in the GCTAPsfTable::read() function ($GAMMALIB/inst/cta/src/GCTAPsfTable.cpp, line 287 in current devel branch). Since gammalib loads the axes by names, the order shouldn’t be an issue.

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

I would have the same suspicion as Michael. Could you post the PSF file you’re using to check the content?

To see if the units are right you may sum over all bins as follows:

sum = sum + PSF(RAD) * sin(RAD) * (RAD_HI - RAD_LO) * 2 * pi
where
RAD = 0.5 * (RAD_LO + RAD_HI)
and the sum should be 1. In this formula sin(RAD) comes from the transformation to polar coordinates. Maybe you have already included the sin(RAD) in your file?

#7 Updated by Tiziani Domenico over 6 years ago

Unfortunately I can not post my PSF file but I changed the file from the gammalib repository accordingly.
You can find the result in the attachment. The problem still occurs with this file.

I checked the sum as described by Knödlseder Jürgen and obtained a value of 1 for each energy and offset bin.

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

I checked the file and confirm that there is a problem with the order of the columns. See the attached plots. The first is with RAD_LO and RAD_HI as first and second columns, the second with RAD_LO and RAD_HI as fifth and sixth columns. Only the later is okay.

So a quick fix on your side is to put the columns in the ENERGY / THETA / RAD order. I change this issue to a bug.

#9 Updated by Tiziani Domenico over 6 years ago

Thank you, Knödlseder Jürgen! I think this is exactly the problem that I encountered.

#10 Updated by Knödlseder Jürgen over 6 years ago

There is one thing coming into my mind: the order of the axis columns relate to the encoding of the 3D array in the RPSF column. In other words, if RAD is the first column, the radial axis must be the most rapidly varying axis of the RPSF array. At least in the example file you sent this is not the case, and that’s the reason why the file was wrong.

How about the file you are using? Have you encoded the array in the order you defined the axis columns?

#11 Updated by Tiziani Domenico over 6 years ago

The file that I am using has exactly the same array encoding and the same order of the axis columns.
You are right, if I change the order to ENERGY / THETA / RAD, everything works fine.

If this relation is a convention for fits files, I should rather adjust the generation of the Table-PSF files.

#12 Updated by Knödlseder Jürgen over 6 years ago

Indeed, this is a convention. The order of the axis columns need to reflect the dimension of the table.

#13 Updated by Knödlseder Jürgen over 6 years ago

  • Status changed from New to Closed
  • Assigned To set to Knödlseder Jürgen
  • Target version set to 1.2.0
  • % Done changed from 0 to 100

I added a check to the GCTAResponseTable class that verifies that the data columns are consistent with the axis columns. I check on the that now an exception is thrown for the test data provided, and that are correct table passes the tests.

Merged into devel.

Also available in: Atom PDF