GammaLib

Michael - sounds good to me. Thanks!

One minor question / comment: why call it “extname” in the code and “hdu” in the XML file ... make it consistent?

One minor question / comment: why call it “extname” in the code and “hdu” in the XML file ... make it consistent?

Don’t have a strong opinion “hdu” is shorter and more handy, “extname” makes it more consistent. Jürgen, any opinion?

Mayer Michael wrote:

One minor question / comment: why call it “extname” in the code and “hdu” in the XML file ... make it consistent?

Don’t have a strong opinion “hdu” is shorter and more handy, “extname” makes it more consistent. Jürgen, any opinion?

I agree to have it consistently named extname would be better.

Ok, I have made the change.
We might need the same behaviour for CTAEventList and GTI?

We might need the same behaviour for CTAEventList and GTI?

Yes, please!

For HESS, IMO it would be nice if I have the option to store the 100k HDUs for all HESS data in one FITS file that’s a few GB in size (depending on the cuts) and distribute that instead of 100k separate files for each HDU or 20k where the HDUs for one obs are bundled.

If you can make Gammalib / ctools so flexible the that HDU names can be specified everywhere / processing data from such a large multi-obs bundled FITS file, that would be great!

Mayer Michael wrote:

Ok, I have made the change.
We might need the same behaviour for CTAEventList and GTI?

I support to have this option always when a FITS file is loaded.

We should create a new feature on this and screen the code to see what needs to be modified (I’m current working on integrating your code, Michael). When working on the new feature, please also think about updating the Doxygen description of the interface (I’m about to add the missing extname parameters).

We should create a new feature on this and screen the code to see what needs to be modified

I have created #1598

When working on the new feature, please also think about updating the Doxygen description of the interface (I’m about to add the missing extname parameters).

Oh right, I forgot about this. Thanks.

I’m in the process to merge the change in, yet I was thinking about another possibility.

The cfitsio convention is to specify an extension name in squared brackets following

filename[extname]

How does this sound (sorry Michael, I should have come up with this though before you start coding).

How does this sound (sorry Michael, I should have come up with this though before you start coding).

No worries
Sounds like a very good idea which just didn’t came to my mind. This would simplify things a lot I guess. Would you then already have this as filename in the xml file?
E.g.

<parameter name="PointSpreadFunction" file="psf.fits.gz[PSF]"/>

• Revert my changes and stick with the current interface in CTA response functions i.e. use load(std::string& filename).
• In CTA response components load() functions, create a GFilename instance from the input string
• Set the default_extname of the GFilename instance to the respective component (e.g. “EFFECTIVE AREA” in GCTAAeff2D)
• Get the FITS table in extension GFilename::extname()
• read the FITS table from the corresponding extension name.
• Remove “extname” from XML interface as this can be handled via the file name parameter

Did I miss something?

Here is an example here how I would modify the code e.g. in GCTAAeff2D::load.

void GCTAAeff2D::load(const std::string& filename)
{
    // Initialise filename
    GFilename fname(filename);

    // Set default extension name
    fname.default_extname("EFFECTIVE AREA");

    // Open effective area FITS file
    GFits file(fname.filename());

    // Initialise FITS table
    GFitsTable& table;

    // Check for extno, otherwise use extname to initialise the table
    if (fname.has_extno()) {
        table = *file.table(fname.extno());
    }
    else {
        table = *file.table(fname.extname());
    }

    // Read effective area
    read(table);

    // Close effective area FITS file
    file.close();

    // Store filename
    m_filename = filename;

    // Return
    return;
}

We could also think about putting some more logic into the function GFits::table, e.g. add GFits::table(GFilename&) or even a new constructor for the table: GFitsTable(GFilename&)?

Mayer Michael wrote:

Thanks for adding GFilename (#1599) - does it need to be GFileName though?.
Just to clarify what is necessary now:

• Revert my changes and stick with the current interface in CTA response functions i.e. use load(std::string& filename).
• In CTA response components load() functions, create a GFilename instance from the input string
• Set the default_extname of the GFilename instance to the respective component (e.g. “EFFECTIVE AREA” in GCTAAeff2D)
• Get the FITS table in extension GFilename::extname()
• read the FITS table from the corresponding extension name.
• Remove “extname” from XML interface as this can be handled via the file name parameter

Did I miss something?

I’m still fiddling around with the GFilename interface. I basically try to use the interface in GEbounds. Still wait a bit (I just removed default_extname).

Mayer Michael wrote:

Thanks for adding GFilename (#1599) - does it need to be GFileName though?.

I have not a very strong opinion about the name.

Sometimes I write filename as a single word, some times as two words (see https://en.wikipedia.org/wiki/Filename).

We should decide before using it if we keep it as is or we want to change. Any votes?

Any votes?

Didn’t know that it can be both one word or two. I also don’t have a strong opinion. I guess I would vote for keeping GFilename. I recall there was a similar discussion about GSkymap vs GSkyMap (which is why I brought up the point in the first place).
You make the call

Here’s how you can use now GFilename (I changed GEbounds as a test case):

void GEbounds::load(const std::string& filename)
{
    // Create file name
    GFilename fname(filename);

    // Allocate FITS file
    GFits file;

    // Open FITS file
    file.open(fname.filename());

    // Get energy boundary table
    const GFitsTable& table = *file.table(fname.extname("EBOUNDS"));

    // Read energy boundaries from table
    read(table);

    // Close FITS file
    file.close();

    // Return
    return;
}

For saving you would do

void GEbounds::save(const std::string& filename,
                    const bool&        clobber,
                    const std::string& unit) const
{
    // Create file name
    GFilename fname(filename);

    // Allocate FITS file
    GFits file;

    // Write energy boundaries to FITS file
    write(file, fname.extname("EBOUNDS"), unit);

    // Save to file
    file.saveto(fname.filename(), clobber);

    // Return
    return;
}

Thanks a lot for the input. I will checkout a fresh branch and make the changes to the CTA response methods.
One question though:
In the example I provided above, we would also take into account the possibility to specify an HDU number instead of a name:

// Check for extno, otherwise use extname to initialise the table
    if (fname.has_extno()) {
        table = *file.table(fname.extno());
    }
    else {
        table = *file.table(fname.extname());
    }

    // Get PSF table
    const GFitsTable& table = *file.table(1);

You are right that we may have the two cases, and I agree that ultimately we should support both. I’m wondering what we should do for release 1.0 and what for beyond. I would limit release 1.0 to what is strictly necessary (extno seems a bit more low level, and maybe no one needs it now), and do other things later with more time for thinking about the best way to do it.

With the actual GFilename class you would need to add a lot of if’s in many methods. This can be avoided by introducing a GFitsExtension class that handles extname and extno (and possibly also extver) in a unified way. GFilename would then fill a GFitsExtension object as needed, and you pass this to the FITS classes instead of having methods for extension number and extension name. But we should do this beyond release 1.0 (otherwise we never meet the date).

I recognized that we should also add default values to the extno and extver methods. If you want you can do that (in the way how it’s done for extname). Make sure you commit this separately from any other change (makes it easier to deal with stuff for release 1.0 and other stuff).

P.S. Will be offline now until tonight.

I added the default arguments to the functions extno() and extver on branch 1599-add-defaults-to-GFilename

Mayer Michael wrote:

I added the default arguments to the functions extno() and extver on branch 1599-add-defaults-to-GFilename

Merged into devel.

I have realised that I forgot to store the filename in the load functions of the response components.
Therefore,I have updated the branch 1597-use-GFilename-in-CTAresponse-components to the latest version and included the changes.
Could you merge that in, too? Sorry for the inconvenience.

• GCTAResponseIrf::load_aeff()
• GCTAResponseIrf::load_psf()
• GCTAResponseIrf::load_edisp()
• GCTAResponseIrf::load_background()

I have updated them on the same branch (including the file inst/cta/test/data/test_1dc.xml, now specifying extension names).

Mayer Michael wrote:

I have realised that I forgot to store the filename in the load functions of the response components.
Therefore,I have updated the branch 1597-use-GFilename-in-CTAresponse-components to the latest version and included the changes.
Could you merge that in, too? Sorry for the inconvenience.

Merged into devel.

Mayer Michael wrote:

I also realised that the following functions need adaption, too:

• GCTAResponseIrf::load_aeff()
• GCTAResponseIrf::load_psf()
• GCTAResponseIrf::load_edisp()
• GCTAResponseIrf::load_background()

I have updated them on the same branch (including the file inst/cta/test/data/test_1dc.xml, now specifying extension names).

I went over your proposed changes which required to add extension names to the inst/cta/test/data/test_1dc.xml file.

I made a slightly different change so that this is not needed. Here an example of how I modified the GCTAResponseIrf::load_aeff method. The main difference is that the default extension name is created on basis of the extensions that exist in the FITS file. Thus if a 1DC file is read without providing the extension, the default extension name will be set to SPECRESP in that case. Please let me know if this looks also fine for you, or if you have trouble with this new implementation.

    // If file is a FITS file ...
    if (gammalib::is_fits(filename)) {

        // Open FITS file
        GFits file(fname.filename());

        // Get the extension name. If an extension name has been specified
        // then use this name, otherwise use either the "EFFECTIVE AREA" 
        // or the "SPECRESP" extension.
        std::string extname = "";
        if (fname.has_extname()) {
            extname = fname.extname();
        }
        else {
            if (file.contains("EFFECTIVE AREA")) {
                extname = "EFFECTIVE AREA";
            }
            else if (file.contains("SPECRESP")) {
                extname = "SPECRESP";
            }
        }

        // Continue only if extension name is not empty
        if (!extname.empty()) {

            // Get FITS table
            const GFitsTable& table = *file.table(extname);

            // Read save energy thresholds if available
            if (table.has_card("LO_THRES")) {
                m_lo_save_thres = table.real("LO_THRES");
            }
            if (table.has_card("HI_THRES")) {
                m_hi_save_thres = table.real("HI_THRES");
            }

            // Check for specific table column
            if (table.contains("SPECRESP")) {

                // Close file
                file.close();

                // Allocate Aeff from file
                m_aeff = new GCTAAeffArf(filename);

            } // endif: load as GCTAAeffArf

            else {

                // Close file
                file.close();

                // Allocate Aeff from file
                m_aeff = new GCTAAeff2D(filename);

            } // endelse: load as GCTAAeff2D

        } // endif: extension name is not empty

        // Signal that no extension was found
        else {
            std::string msg = "FITS file \""+filename+"\" does not " 
                              "contain a valid effective area table. " 
                              "Please specify a valid effective area " 
                              "response file.";
            throw GException::file_error(G_LOAD_AEFF, msg);
        }

    } // endif: file was FITS file

Excellent, the code looks fine with me. Much better now as we do not have to specify extension names for the 1DC format. Thanks!