From d16893f961382c67a7c2a8cfbc7431b1227fc453 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 18 Jan 2016 23:23:29 +0100 Subject: [PATCH] A gti_extname() method has been added to GCTAObservation and GCTAEventList to set the name of the Good Time Intervals FITS extension. The data sub-space keyword is now also adjusted to the specified GTI extension name when the data are written to disk. This satisfies the feature request #1598. --- NEWS | 5 +++++ inst/cta/include/GCTAEventList.hpp | 45 ++++++++++++++++++++++++++++++++++++++------- inst/cta/include/GCTAObservation.hpp | 2 +- inst/cta/pyext/GCTAEventList.i | 16 +++++++++------- inst/cta/pyext/GCTAObservation.i | 1 + inst/cta/src/GCTAEventList.cpp | 104 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------- inst/cta/src/GCTAObservation.cpp | 178 +++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------- inst/cta/test/test_CTA.cpp | 44 +++++++++++++++++++++++++++++++++++++++----- 8 files changed, 205 insertions(+), 190 deletions(-) diff --git a/NEWS b/NEWS index 2d3316b..b9ac204 100644 --- a/NEWS +++ b/NEWS @@ -110,6 +110,11 @@ GTIs (and also event lists) in a single FITS file. We do however require that the event and GTI extension reside in the same FITS file. This satisfies the feature request #1598. +A gti_extname() method has been added to GCTAObservation and GCTAEventList to +set the name of the Good Time Intervals FITS extension. The data sub-space +keyword is now also adjusted to the specified GTI extension name when the +data are written to disk. This satisfies the feature request #1598. + 21. COMPTEL interface --------------------- diff --git a/inst/cta/include/GCTAEventList.hpp b/inst/cta/include/GCTAEventList.hpp index 24676a9..87ad7e1 100644 --- a/inst/cta/include/GCTAEventList.hpp +++ b/inst/cta/include/GCTAEventList.hpp @@ -69,19 +69,21 @@ public: virtual void load(const std::string& filename); virtual void save(const std::string& filename, const bool& clobber = false) const; - virtual void read(const GFits& file); - virtual void write(GFits& file) const; + virtual void read(const GFits& fits); + virtual void write(GFits& fits) const; virtual int number(void) const; virtual void roi(const GRoi& roi); virtual const GCTARoi& roi(void) const; std::string print(const GChatter& chatter = NORMAL) const; // Implement other methods - void append(const GCTAEventAtom& event); - void reserve(const int& number); - double irf_cache(const std::string& name, const int& index) const; - void irf_cache(const std::string& name, const int& index, - const double& irf) const; + void append(const GCTAEventAtom& event); + void reserve(const int& number); + double irf_cache(const std::string& name, const int& index) const; + void irf_cache(const std::string& name, const int& index, + const double& irf) const; + void gti_extname(const std::string& extname); + std::string gti_extname(void) const; protected: // Protected methods @@ -172,4 +174,33 @@ void GCTAEventList::reserve(const int& number) return; } + +/***********************************************************************//** + * @brief Return Good Time Intervals FITS extension name + * + * @return Good Time Intervals FITS extension name. + * + * Returns the FITS extension name of the Good Time Intervals. + ***************************************************************************/ +inline +std::string GCTAEventList::gti_extname(void) const +{ + return ((m_gti_extname.empty()) ? "GTI" : m_gti_extname); +} + + +/***********************************************************************//** + * @brief Set Good Time Intervals FITS extension name + * + * @param[in] extname Good Time Intervals FITS extension name. + * + * Sets the FITS extension name for the Good Time Intervals. + ***************************************************************************/ +inline +void GCTAEventList::gti_extname(const std::string& extname) +{ + m_gti_extname = extname; + return; +} + #endif /* GCTAEVENTLIST_HPP */ diff --git a/inst/cta/include/GCTAObservation.hpp b/inst/cta/include/GCTAObservation.hpp index 5bb271d..657bc41 100644 --- a/inst/cta/include/GCTAObservation.hpp +++ b/inst/cta/include/GCTAObservation.hpp @@ -128,6 +128,7 @@ public: void eventfile(const std::string& filename); const std::string& eventfile(void) const; const std::string& eventtype(void) const; + void gti_extname(const std::string& extname); void dispose_events(void); const double& lo_user_thres(void) const; const double& hi_user_thres(void) const; @@ -146,7 +147,6 @@ protected: // Protected members std::string m_instrument; //!< Instrument name std::string m_eventfile; //!< Event filename - //std::string m_gtifile; //!< GTI filename std::string m_eventtype; //!< Event type (for XML file) GCTAResponse* m_response; //!< Pointer to instrument response functions GCTAPointing m_pointing; //!< Pointing direction diff --git a/inst/cta/pyext/GCTAEventList.i b/inst/cta/pyext/GCTAEventList.i index 2900dbb..36a28c2 100644 --- a/inst/cta/pyext/GCTAEventList.i +++ b/inst/cta/pyext/GCTAEventList.i @@ -51,18 +51,20 @@ public: virtual void load(const std::string& filename); virtual void save(const std::string& filename, const bool& clobber = false) const; - virtual void read(const GFits& file); - virtual void write(GFits& file) const; + virtual void read(const GFits& fits); + virtual void write(GFits& fits) const; virtual int number(void) const; virtual void roi(const GRoi& roi); virtual const GCTARoi& roi(void) const; // Implement other methods - void append(const GCTAEventAtom& event); - void reserve(const int& number); - double irf_cache(const std::string& name, const int& index) const; - void irf_cache(const std::string& name, const int& index, - const double& irf) const; + void append(const GCTAEventAtom& event); + void reserve(const int& number); + double irf_cache(const std::string& name, const int& index) const; + void irf_cache(const std::string& name, const int& index, + const double& irf) const; + void gti_extname(const std::string& extname); + std::string gti_extname(void) const; }; diff --git a/inst/cta/pyext/GCTAObservation.i b/inst/cta/pyext/GCTAObservation.i index d654444..a89e01e 100644 --- a/inst/cta/pyext/GCTAObservation.i +++ b/inst/cta/pyext/GCTAObservation.i @@ -99,6 +99,7 @@ public: void eventfile(const std::string& filename); const std::string& eventfile(void) const; const std::string& eventtype(void) const; + void gti_extname(const std::string& extname); void dispose_events(void); const double& lo_user_thres(void) const; const double& hi_user_thres(void) const; diff --git a/inst/cta/src/GCTAEventList.cpp b/inst/cta/src/GCTAEventList.cpp index 12cb144..b3e36f7 100644 --- a/inst/cta/src/GCTAEventList.cpp +++ b/inst/cta/src/GCTAEventList.cpp @@ -293,14 +293,14 @@ void GCTAEventList::load(const std::string& filename) void GCTAEventList::save(const std::string& filename, const bool& clobber) const { - // Create empty FITS file - GFits fits; + // Open or create FITS file + GFits fits(filename, true); // Write event list write(fits); // Save FITS file - fits.saveto(filename, clobber); + fits.save(clobber); // Return return; @@ -335,11 +335,8 @@ void GCTAEventList::read(const GFits& fits) // Clear object clear(); - // Initialise filename from fits file - GFilename fname = GFilename(fits.filename()); - - // Initialise extension name - std::string extname = fname.extname("EVENTS"); + // Initialise events extension name + std::string extname = fits.filename().extname("EVENTS"); // Get event list HDU const GFitsTable& events = *fits.table(extname); @@ -400,31 +397,37 @@ void GCTAEventList::read(const GFits& fits) /***********************************************************************//** * @brief Write CTA events into FITS file. * - * @param[in] file FITS file. + * @param[in] fits FITS file. + * + * Writes the CTA event list into a FITS file. The events will be written by + * default into the extension "EVENTS" unless an extension name is explicitly + * specified in the FITS file name. The method also writes the data sub-space + * keywords in the FITS header of the events table. * - * Write the CTA event list into FITS file. + * In addition, the method will also append a table containing the Good Time + * Intervals of the events to the FITS file. The extension name for the Good + * Time Intervals is either taken from the m_gti_extname member, or if empty, + * is set to "GTI". ***************************************************************************/ -void GCTAEventList::write(GFits& file) const +void GCTAEventList::write(GFits& fits) const { - // Allocate FITS binary table HDU - GFitsBinTable* events = new GFitsBinTable; + // Allocate empty FITS binary table + GFitsBinTable table; - // Write events - write_events(*events); + // Write events into binary table + write_events(table); + + // Set FITS extension name (default: "EVENTS") + table.extname(fits.filename().extname("EVENTS")); // Write data selection keywords - write_ds_keys(*events); + write_ds_keys(table); // Append event table to FITS file - file.append(*events); + fits.append(table); - // Free binary table - delete events; - - // Write GTI extension if we have an extension name - if (!m_gti_extname.empty()) { - gti().write(file, m_gti_extname); - } + // Write GTI extension + gti().write(fits, gti_extname()); // Return return; @@ -895,16 +898,32 @@ void GCTAEventList::write_events(GFitsBinTable& hdu) const /***********************************************************************//** - * @brief Write data selection keywords into FITS HDU + * @brief Write data sub-space keywords into FITS HDU * * @param[in] hdu FITS HDU. * - * This method does nothing if the HDU pointer is NULL. + * Writes the data sub-space keywords for an event list into the FITS HDU. + * The following keywords will be written: + * + * DSTYP1 = "TIME" / Data sub-space type + * DSUNI1 = "s" / Data sub-space unit + * DSVAL1 = "TABLE" / Data sub-space value + * DSREF1 = ":[extname]" / Data sub-space reference + * DSTYP2 = "ENERGY" / Data sub-space type + * DSUNI2 = "TeV" / Data sub-space unit + * DSVAL2 = "[emin]:[emax]" / Data sub-space value + * DSTYP3 = "POS(RA,DEC)" / Data sub-space type + * DSUNI3 = "deg" / Data sub-space unit + * DSVAL3 = "CIRCLE([ra],[dec],[rad])" / Data sub-space value * - * @todo This is a very dumb data selection keyword writing routine that does - * not take into account any existing keywords. We definitely want a - * more secure logic that checks for existing keywords and possible - * conflicts. But for the moment, this code does the job. + * where + * + * [extname] is the GTI extension (default: GTI) + * [emin] is the minimum event energy in TeV + * [emax] is the maximum event energy in TeV + * [ra] is the Right Ascension of the Region of Interest centre in degrees + * [dec] is the Declination of the Region of Interest centre in degrees + * [rad] is the radius of the Region of Interest in degrees ***************************************************************************/ void GCTAEventList::write_ds_keys(GFitsHDU& hdu) const { @@ -921,16 +940,19 @@ void GCTAEventList::write_ds_keys(GFitsHDU& hdu) const std::string dsval2 = gammalib::str(e_min) + ":" + gammalib::str(e_max); + // Set Good Time Intervals extension name + std::string dsref1 = ":"+gti_extname(); + // Add time selection keywords - hdu.card("DSTYP1", "TIME", "Data selection type"); - hdu.card("DSUNI1", "s", "Data selection unit"); - hdu.card("DSVAL1", "TABLE", "Data selection value"); - hdu.card("DSREF1", ":GTI", "Data selection reference"); + hdu.card("DSTYP1", "TIME", "Data sub-space type"); + hdu.card("DSUNI1", "s", "Data sub-space unit"); + hdu.card("DSVAL1", "TABLE", "Data sub-space value"); + hdu.card("DSREF1", dsref1, "Data sub-space reference"); // Add energy range selection - hdu.card("DSTYP2", "ENERGY", "Data selection type"); - hdu.card("DSUNI2", "TeV", "Data selection unit"); - hdu.card("DSVAL2", dsval2, "Data selection value"); + hdu.card("DSTYP2", "ENERGY", "Data sub-space type"); + hdu.card("DSUNI2", "TeV", "Data sub-space unit"); + hdu.card("DSVAL2", dsval2, "Data sub-space value"); // Initialise number of NDSKEYS int ndskeys = 2; @@ -945,15 +967,15 @@ void GCTAEventList::write_ds_keys(GFitsHDU& hdu) const gammalib::str(rad) + ")"; // Write DS keywords - hdu.card("DSTYP3", "POS(RA,DEC)", "Data selection type"); - hdu.card("DSUNI3", "deg", "Data selection unit"); - hdu.card("DSVAL3", dsval3, "Data selection value"); + hdu.card("DSTYP3", "POS(RA,DEC)", "Data sub-space type"); + hdu.card("DSUNI3", "deg", "Data sub-space unit"); + hdu.card("DSVAL3", dsval3, "Data sub-space value"); ndskeys++; } // endif: RoI was valid // Set number of data selection keys - hdu.card("NDSKEYS", ndskeys, "Number of data selections"); + hdu.card("NDSKEYS", ndskeys, "Number of data sub-space keys"); // Return return; diff --git a/inst/cta/src/GCTAObservation.cpp b/inst/cta/src/GCTAObservation.cpp index ee2dc78..566046e 100644 --- a/inst/cta/src/GCTAObservation.cpp +++ b/inst/cta/src/GCTAObservation.cpp @@ -480,37 +480,9 @@ GGti GCTAObservation::gti(void) const // ... otherwise try reading the GTIs from file else { - // If we know in which file the GTIs reside then extract the GTIs - // from that file ... - /* - if (!m_gtifile.empty()) { - - // Initialise filename - GFilename fname(m_gtifile); - - // Open FITS file - GFits fits(fname.filename()); - - // Get extension name - std::string extname = fname.extname("GTI"); - - // If the FITS file contains the specified extension, then read - // Good Time Intervals from that extension - if (fits.contains(extname)) { - const GFitsTable& table = *fits.table(extname); - gti.read(table); - } - - // Close FITS file - fits.close(); - - } // endif: gti file name was not not empty - - // ... otherwise we extract the GTIs on basis on the DSS keywords, - // and if these keywords are not present, we build the GTIs from - // the start and stop time that is provided in the event header - else if (!m_eventfile.empty()) { - */ + // Extract the GTIs on basis on the DSS keywords, and if these + // keywords are not present, we build the GTIs from the start and + // stop time that is provided in the event header if (!m_eventfile.empty()) { // Initialise filename @@ -596,13 +568,6 @@ GGti GCTAObservation::gti(void) const // Check for GTI validity if (gti.size() == 0) { - /* - std::string msg = "Could not find Good Time Interval information " - "in the event file \""+m_eventfile+"\" or in " - "the GTI file "+m_gtifile+"\". Please " - "check that Good Time Interval information is " - "properly set in the event file."; - */ std::string msg = "Could not find Good Time Interval information " "in the event file \""+m_eventfile+"\". Please " "check that Good Time Interval information is " @@ -858,16 +823,6 @@ void GCTAObservation::read(const GXmlElement& xml) } - // Read optional Good Time Intervals file name - /* - else if (par->attribute("name") == "GoodTimeIntervals") { - - // Read GTI filename - m_gtifile = par->attribute("file"); - - } - */ - } // endfor: looped over observation parameters // Analyse parameters @@ -1113,15 +1068,6 @@ void GCTAObservation::write(GXmlElement& xml) const GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, m_eventtype); par->attribute("file", m_eventfile); - // If we had a GTI file name upon reading from a XML file then write - // back the GTI file name information into the XML file - /* - if (!m_gtifile.empty()) { - GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "GoodTimeIntervals"); - par->attribute("file", m_gtifile); - } - */ - } // ... otherwise write the observation definition information @@ -1155,9 +1101,9 @@ void GCTAObservation::write(GXmlElement& xml) const /***********************************************************************//** - * @brief Read data from FITS object + * @brief Read data from FITS file * - * @param[in] fits FITS object. + * @param[in] fits FITS file. ***************************************************************************/ void GCTAObservation::read(const GFits& fits) { @@ -1218,11 +1164,23 @@ void GCTAObservation::read(const GFits& fits) /***********************************************************************//** - * @brief Write CTA observation into FITS object. + * @brief Write observation into FITS file. * - * @param[in] fits FITS object. + * @param[in] fits FITS file. + * @param[in] extname FITS extension name for event list (default: "EVENTS"). * - * This method does nothing if no events are within the CTA observation. + * Writes the observation into a FITS file. + * + * If the observation contains an event list, the event list and Good Time + * Interval will be written to the FITS file. The FITS extension name for the + * event list can be optionally specified by the @p extname argument. If the + * extension name of the GTI should be changed to a non-default value, the + * gti_extname() method should be called prior to calling the write() method. + * + * If the observation contains an event cube, the event cube will be written + * into the primary extension of the FITS file. + * + * This method does nothing if no events are in the CTA observation. ***************************************************************************/ void GCTAObservation::write(GFits& fits, const std::string& extname) const { @@ -1237,8 +1195,9 @@ void GCTAObservation::write(GFits& fits, const std::string& extname) const // the GTI if provided on loading. list->write(fits); - // Get reference to EVENTS HDU - GFitsHDU& hdu = *fits.at("EVENTS"); + // Get reference to EVENTS HDU (see GCTAEventList::write() for + // consistency) + GFitsHDU& hdu = *fits.at(fits.filename().extname("EVENTS")); // Change extension name hdu.extname(extname); @@ -1294,34 +1253,6 @@ void GCTAObservation::load(const std::string& filename) /***********************************************************************//** - * @brief Load GTIs FITS file. - * - * @param[in] filename Name of FITS file from which GTIs are loaded. - * - * Loads CTA GTIs from a separate file. - ***************************************************************************/ -/* -void GCTAObservation::load_gti(const std::string& filename) -{ - - // Check if we have an event list - GCTAEventList* list = dynamic_cast(m_events); - - // Load GTI - if (list != NULL) { - list->load_gti(filename); - } - - // Store GTI file name - m_gtifile = filename; - - // Return - return; -} -*/ - - -/***********************************************************************//** * @brief Load stacked cube from FITS files * * @param[in] cntcube Counts cube file name. @@ -1368,7 +1299,7 @@ void GCTAObservation::load(const std::string& cntcube, * @brief Save CTA observation into FITS file. * * @param[in] filename FITS filename. - * @param[in] clobber Overwrite existing FITS file (default=false). + * @param[in] clobber Overwrite existing FITS file? (default: false) ***************************************************************************/ void GCTAObservation::save(const std::string& filename, const bool& clobber) const { @@ -1393,35 +1324,6 @@ void GCTAObservation::save(const std::string& filename, const bool& clobber) con /***********************************************************************//** - * @brief Save CTA GTIs into FITS file. - * - * @param[in] filename FITS filename. - * @param[in] clobber Overwrite existing FITS file (default=false). - ***************************************************************************/ -/* -void GCTAObservation::save_gti(const std::string& filename, const bool& clobber) const -{ - // Initialise filename - GFilename fname(filename); - - // Initialise extension name - std::string extname = fname.extname("GTI"); - - // Initialise FITS file - GFits fits(fname.filename(), true); - - // Write GTIs to fits files - gti().write(fits, extname); - - // Save FITS file - fits.save(clobber); - - // Return - return; -} -*/ - -/***********************************************************************//** * @brief Set event container * * @param[in] events Event container. @@ -1468,11 +1370,6 @@ const GEvents* GCTAObservation::events(void) const { try { const_cast(this)->load(m_eventfile); - /* - if (!m_gtifile.empty()) { - const_cast(this)->load_gti(m_gtifile); - } - */ } catch (...) { ; @@ -1503,6 +1400,31 @@ const GEvents* GCTAObservation::events(void) const /***********************************************************************//** + * @brief Set Good Time Intervals extension name + * + * @param[in] extname Good Time Intervals extension name. + * + * Sets the FITS extension name of the Good Time Intervals. This method + * should be used prior to saving or writing of the FITS file if the Good + * Time Intervals should be used in an extension that differs from the + * default extension ("GTI") + ***************************************************************************/ +void GCTAObservation::gti_extname(const std::string& extname) +{ + // Get non-const pointer to event list + GCTAEventList* list = dynamic_cast(const_cast(events())); + + // If we have a CTA event list then set the extension name + if (list != NULL) { + list->gti_extname(extname); + } + + // Return + return; +} + + +/***********************************************************************//** * @brief Dispose events * * Drops the events from the observation. Be careful with using this method @@ -1619,7 +1541,6 @@ void GCTAObservation::init_members(void) m_instrument = "CTA"; m_object.clear(); m_eventfile.clear(); - //m_gtifile.clear(); m_eventtype.clear(); m_bgdfile.clear(); m_response = NULL; @@ -1650,7 +1571,6 @@ void GCTAObservation::copy_members(const GCTAObservation& obs) m_instrument = obs.m_instrument; m_object = obs.m_object; m_eventfile = obs.m_eventfile; - //m_gtifile = obs.m_gtifile; m_eventtype = obs.m_eventtype; m_bgdfile = obs.m_bgdfile; m_pointing = obs.m_pointing; diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 8f2fef7..6772e4d 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -1297,10 +1297,6 @@ void TestGCTAModel::test_model_aeff_bgd(void) ***************************************************************************/ void TestGCTAObservation::test_unbinned_obs(void) { - // Set filenames - const std::string file1 = "test_cta_obs_unbinned.xml"; - const std::string file2 = "test_cta_obs_unbinned_gti.xml"; - // Declare observations GObservations obs; GCTAObservation run; @@ -1382,7 +1378,7 @@ void TestGCTAObservation::test_unbinned_obs(void) test_try("Test XML loading and saving"); try { obs = GObservations(cta_unbin_xml); - obs.save(file1); + obs.save("test_cta_obs_unbinned.xml"); GCTAObservation* run = dynamic_cast(obs[0]); test_value(run->roi().centre().dir().ra_deg(), 83.63); test_value(run->roi().centre().dir().dec_deg(), 22.01); @@ -1405,6 +1401,44 @@ void TestGCTAObservation::test_unbinned_obs(void) test_try_failure(e); } + // Test loading of the event file and saving into a different extension + run.load(cta_events); + run.gti_extname("GTI2"); + run.save("test_cta_events1.fits[EVENTS2]", true); + GFits fits("test_cta_events1.fits"); + test_value(fits.size(), 3, "FITS file should contain 3 HDUs"); + test_assert(fits.contains("EVENTS2"), "FITS should contain \"EVENTS2\" HDU"); + test_assert(fits.contains("GTI2"), "FITS should contain \"GTI2\" HDU"); + fits.close(); + run.gti_extname("GTI3"); + run.save("test_cta_events1.fits[EVENTS3]", true); + fits.open("test_cta_events1.fits"); + test_value(fits.size(), 3, "FITS file should contain 3 HDUs"); + test_assert(fits.contains("EVENTS3"), "FITS should contain \"EVENTS3\" HDU"); + test_assert(fits.contains("GTI3"), "FITS should contain \"GTI3\" HDU"); + fits.close(); + + // Test writing of multiple event files into FITS file + run.load(cta_events); + GFits fits2; + run.gti_extname("GTI1"); + run.write(fits2, "EVENTS1"); + run.gti_extname("GTI2"); + run.write(fits2, "EVENTS2"); + run.gti_extname("GTI3"); + run.write(fits2, "EVENTS3"); + fits2.saveto("test_cta_events2.fits", true); + fits2.close(); + fits.open("test_cta_events2.fits"); + test_value(fits.size(), 7, "FITS file should contain 7 HDUs"); + test_assert(fits.contains("EVENTS1"), "FITS should contain \"EVENTS1\" HDU"); + test_assert(fits.contains("GTI1"), "FITS should contain \"GTI1\" HDU"); + test_assert(fits.contains("EVENTS2"), "FITS should contain \"EVENTS2\" HDU"); + test_assert(fits.contains("GTI2"), "FITS should contain \"GTI2\" HDU"); + test_assert(fits.contains("EVENTS3"), "FITS should contain \"EVENTS3\" HDU"); + test_assert(fits.contains("GTI3"), "FITS should contain \"GTI3\" HDU"); + fits.close(); + // Exit test return; -- libgit2 0.21.4