ctools

1. Provide unbinned observations and do an unbinned analyses in each energy bin (as is working now)
2. Provide an unbinned observation and run binned analysis in each energy bin (as is working now). This might be important if the unbinned analysis takes much too long (even in a small energy bin)
3. Provide a binned observation: I am not sure how to run an analysis in this case as an internal rebinning is not possible (only by merging bins, is that what you mean?). But if we assume a spectral shape in each energy bin, we would probably need more than 1 energy bin per spectrum point, right?

If I understand correctly, your plan is to remove 2) and implement 3) instead, right? I agree this would very much simplify the parameter interface. The only concern I have is that 2) will be important for large data sets (e.g. in Fermi-LAT this is done quite frequently).

Mayer Michael wrote:

You are right, there is an inconsistency in the parameters. I recall when writing the code, I wanted to implement that an exception is thrown if a binned observation is given. Don't know why but apparently I forgot this.
Your change request sounds great, but I am still a bit uncertain, how to achieve an internal binning from a binned observation.
I just want to outline the different needed use cases from my point of view.

1. Provide unbinned observations and do an unbinned analyses in each energy bin (as is working now)
2. Provide an unbinned observation and run binned analysis in each energy bin (as is working now). This might be important if the unbinned analysis takes much too long (even in a small energy bin)
3. Provide a binned observation: I am not sure how to run an analysis in this case as an internal rebinning is not possible (only by merging bins, is that what you mean?). But if we assume a spectral shape in each energy bin, we would probably need more than 1 energy bin per spectrum point, right?

I won’t merge bins. You can blank out bins from a counts cube by setting their values to negative values (see what ctcubemask does). So we can walk over all layers of a cube, blank out all layers other then the selected layer, and do the fit.

One could generalize this to using always 2 layers, 3 layers, etc. dependent on what is provided by the enumbins parameter.

If I understand correctly, your plan is to remove 2) and implement 3) instead, right? I agree this would very much simplify the parameter interface. The only concern I have is that 2) will be important for large data sets (e.g. in Fermi-LAT this is done quite frequently).

Yes, as 3 gives you more control than 2 (in 2 things are done internally, hence you cannot for example blank out specific regions in a counts cube).

Knödlseder Jürgen wrote:

Mayer Michael wrote:

You are right, there is an inconsistency in the parameters. I recall when writing the code, I wanted to implement that an exception is thrown if a binned observation is given. Don't know why but apparently I forgot this.
Your change request sounds great, but I am still a bit uncertain, how to achieve an internal binning from a binned observation.
I just want to outline the different needed use cases from my point of view.

1. Provide unbinned observations and do an unbinned analyses in each energy bin (as is working now)
2. Provide an unbinned observation and run binned analysis in each energy bin (as is working now). This might be important if the unbinned analysis takes much too long (even in a small energy bin)
3. Provide a binned observation: I am not sure how to run an analysis in this case as an internal rebinning is not possible (only by merging bins, is that what you mean?). But if we assume a spectral shape in each energy bin, we would probably need more than 1 energy bin per spectrum point, right?

I won’t merge bins. You can blank out bins from a counts cube by setting their values to negative values (see what ctcubemask does). So we can walk over all layers of a cube, blank out all layers other then the selected layer, and do the fit.

But I admit that this is probably time consuming, hence we could also think about extracting bins.

One could generalize this to using always 2 layers, 3 layers, etc. dependent on what is provided by the enumbins parameter.

If I understand correctly, your plan is to remove 2) and implement 3) instead, right? I agree this would very much simplify the parameter interface. The only concern I have is that 2) will be important for large data sets (e.g. in Fermi-LAT this is done quite frequently).

Yes, as 3 gives you more control than 2 (in 2 things are done internally, hence you cannot for example blank out specific regions in a counts cube).

I won’t merge bins. You can blank out bins from a counts cube by setting their values to negative values (see what ctcubemask does). So we can walk over all layers of a cube, blank out all layers other then the selected layer, and do the fit.

Now I understand. I wasn’t aware ctcubemask can also make an energy selection (I thought it works only spatial.

But I admit that this is probably time consuming, hence we could also think about extracting bins.

One would need to check how long ctcubemask runs on a counts cube (but I guess it is not longer than ctselect on an observation definition file with a lot of event lists)

One could generalize this to using always 2 layers, 3 layers, etc. dependent on what is provided by the enumbins parameter.

Ok, so if the user provides a counts cube with 20 energy bins and wants to create 25 spectral points, we would run into troubles of course. On the other hand, 20 energy bins would only allow for 20 (or 10 or 5) spectral points without losing information, right? We would have to think how to handle cases like 20 energy bins in the cube and the user wants to use 3 layers (do we create 6 spectral points and omit the 2 last energy bins in the cube?).

Nevertheless, I am a big fan of making the spectral points binning dependent on the binning of the input counts cube.

Mayer Michael wrote:

One could generalize this to using always 2 layers, 3 layers, etc. dependent on what is provided by the enumbins parameter.

Ok, so if the user provides a counts cube with 20 energy bins and wants to create 25 spectral points, we would run into troubles of course. On the other hand, 20 energy bins would only allow for 20 (or 10 or 5) spectral points without losing information, right? We would have to think how to handle cases like 20 energy bins in the cube and the user wants to use 3 layers (do we create 6 spectral points and omit the 2 last energy bins in the cube?).

Right, the translation can only be done approximate. The finest binning is driven by the binning of the counts cube, the number of layers used will be computed using

int n_layers = n_bins_in_counts_cube / enumbins;

You are right that we have to define the behavior for the case where the division has a remainder. I probably would just cut off the thing at the end (i.e. use for example only 1 bin even if the division above would suggest 3).

Ok that sounds reasonable. Accordingly, if the user really wants a finer binning, the data products (countcube, psf, bkg, exp) need to be precomputed with a finer energy binning, too. This makes sense to me.

I have implemented the change you proposed on branch 1566-csspec-update-for-binned-mode.

A slight modification to ctcubemask was necessary to allow for undefined parameters, i.e. making only an energy selection and no RoI selection. This is analogous to ctselect.

I also created a test case for the binned mode. Note that I had to change the header of test/test_cscripts.sh.in: it now only deletes .fits files that were created from the cscripts. The reason is that binned csspec needs input cubes that are not available unless the ctools tests were run beforehand. We might think of permanently adding such cubes to $CTOOLS/test/data?

Mayer Michael wrote:

I have implemented the change you proposed on branch 1566-csspec-update-for-binned-mode.

A slight modification to ctcubemask was necessary to allow for undefined parameters, i.e. making only an energy selection and no RoI selection. This is analogous to ctselect.

I also created a test case for the binned mode. Note that I had to change the header of test/test_cscripts.sh.in: it now only deletes .fits files that were created from the cscripts. The reason is that binned csspec needs input cubes that are not available unless the ctools tests were run beforehand. We might think of permanently adding such cubes to $CTOOLS/test/data?

Indeed, I’d prefer to add test data as otherwise one test depends on another and when something breaks we don’t know who is the actual culprit.

Indeed, I’d prefer to add test data as otherwise one test depends on another and when something breaks we don’t know who is the actual culprit.

We would need a counts cube, exposure cube, psf cube and background cube and a corresponding model XML file. Ideally the cubes would have a rather coarse binning to speed up the unit tests.

Mayer Michael wrote:

Indeed, I’d prefer to add test data as otherwise one test depends on another and when something breaks we don’t know who is the actual culprit.

We would need a counts cube, exposure cube, psf cube and background cube and a corresponding model XML file. Ideally the cubes would have a rather coarse binning to speed up the unit tests.

There are these cubes in gammalib/inst/cta/test/data but I have not checked if the cubes are adapted. Also the exposure cube is a bit large (I probably should try to reduce the size). We should try to make cubes as small as possible as we only want to test the functionality.

The cubes which come out of the unit tests are reasonable in size. After running test_ctools.sh there are files called cntcube2.fits, expcube2.fits, psfcube2.fits, bkgcube2.fits, bkgcube2.xml. The cubes are in the 60KB region and seem to be suitable.

Mayer Michael wrote:

The cubes which come out of the unit tests are reasonable in size. After running test_ctools.sh there are files called cntcube2.fits, expcube2.fits, psfcube2.fits, bkgcube2.fits, bkgcube2.xml. The cubes are in the 60KB region and seem to be suitable.

Indeed. I propose to rename and gzip the files necessary and add them to the test/data folder.

I’ve added the files to the branch 1566-csspec-update-for-binned-mode

Great.
For the moment, we should also note in the documentation that the Npred column will not be filled for binned analyses. The reason is that GObservation::npred (which calls GCTAResponsCube::nroi) is not implemented for binned CTA observations.

Mayer Michael wrote:

Great.
For the moment, we should also note in the documentation that the Npred column will not be filled for binned analyses. The reason is that GObservation::npread (which calls GCTAResponsCube::nroi) is not implemented for binned CTA observations.

Indeed, I have not really changed the text but removed the unused parameters. Do not hesitate to also update the text. User documentation is as important as code!

I added some sentences to the doc.
I have pushed it to the above branch.

Merged into devel.