Allow energy boundaries and energy value FITS file as input to ebinfile parameters (#1978)

Jürgen Knödlseder
1 parent 9617e53e
Showing 7 changed files with 98 additions and 108 deletions
ChangeLog
NEWS
src/ctmapcube/ctmapcube.cpp
src/ctmapcube/ctmapcube.hpp
src/support/ctool.cpp
test/test_cscript.py
test/test_ctmapcube.py
-2017-03-22
+2017-03-24
  
         * Version 1.3.0 released
         ========================
  
+        Support energy boundary FITS file in ctmapcube (#1978)
         Add csebins script (#1976)
-        Take into account RoI in ctskymap and use numerical integration for IRF (#1975)
+        Consider RoI in ctskymap and use numerical integration for IRF (#1975)
         Set local language to english (#1370)
         Enhance logfile output of csroot2caldb script
         Set default deadc value to 0.98 (CTA has <2% deadtime)
 New Features and Important Changes in ctools 1.3.0.dev1
  
-22 March 2017
+24 March 2017
  
  
 Introduction
@@ -87,7 +87,8 @@ None
  
 ctmapcube - CTA map cube generation
 -----------------------------------
-None
+The tool now supports loading the energy binning definition from an energy
+boundary FITS file, behaving now similar to the other cube tools (#1978).
  
  
 ctcubemask - CTA cube filtering
@@ -349,7 +349,7 @@ void ctmapcube::get_parameters(void)
     } // endif: there were no models
  
     // Get energy binning parameters
-    GEnergies energies = create_energies();
+    GEnergies energies(create_ebounds());
  
     // Get spatial binning parameters
     double      xref     = (*this)["xref"].real();
@@ -386,68 +386,6 @@ void ctmapcube::get_parameters(void)
  
  
 /***********************************************************************//**
- * @brief Create map cube energies from user parameters
- *
- * Creates the map cube energies according to the user parameters. If the
- * @p ebinfile parameter is "FILE" the method loads the energy information
- * from a FITS file. The name of the FITS file is given by the @p ebinfile
- * parameter. If no extension name is specified the method will use assume
- * the name to be `ENERGIES`.
- *
- * If the @p ebinfile parameter is "LIN" the method will create a linearly
- * spaced series of @p enumbins energies in the interval @p emin to @p emax.
- * Otherwise, the @p ebinfile parameter is assumed to be "LOG" and the method
- * will create a logarithmically spaced series of @p enumbins energies in the
- * interval @p emin to @p emax.
- ***************************************************************************/
-GEnergies ctmapcube::create_energies(void)
-{
-    // Allocate energies
-    GEnergies energies;
-
-    // Get energy binning algorithm
-    std::string ebinalg = (*this)["ebinalg"].string();
-
-    // If energy binning algorithm is of type "FILE" (case sensitive), then
-    // read energies from FITS file ...
-    if (ebinalg == "FILE") {
-
-        // Get filename
-        GFilename ebinfile = (*this)["ebinfile"].filename();
-
-        // Load energies
-        energies.load(ebinfile);
-
-    } // endif: ebinalg was "FILE"
-
-    // ... otherwise use a linear or a logarithmically-spaced energy binning
-    else {
-
-        // Get task parameters
-        double emin     = (*this)["emin"].real();
-        double emax     = (*this)["emax"].real();
-        int    enumbins = (*this)["enumbins"].integer();
-
-        // Initialise log mode for ebinning
-        bool log = true;
-
-        // Check if algorithm is linear
-        if (ebinalg == "LIN") {
-            log = false;
-        }
-
-        // Setup energies
-        energies = GEnergies(enumbins, GEnergy(emin, "TeV"),
-                                       GEnergy(emax, "TeV"), log);
-
-    } // endelse: ebinalg was not "FILE"
-
-    // Return energies
-    return energies;
-}
-
-
-/***********************************************************************//**
  * @brief Generate map cube
  *
  * Generate map cube by looping over all sky models in the container.
 /***************************************************************************
  *                  ctmapcube - Map cube generation tool                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2016 by Juergen Knoedlseder                              *
+ *  copyright (C) 2016-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -34,7 +34,7 @@
  
 /* __Definitions _________________________________________________________ */
 #define CTMAPCUBE_NAME    "ctmapcube"
-#define CTMAPCUBE_VERSION "1.2.0"
+#define CTMAPCUBE_VERSION "1.3.0"
  
  
 /***********************************************************************//**
@@ -64,17 +64,16 @@ public:
  
 protected:
     // Protected methods
-    void      init_members(void);
-    void      copy_members(const ctmapcube& app);
-    void      free_members(void);
-    void      get_parameters(void);
-    GEnergies create_energies(void);
-    void      create_cube(void);
-    void      add_model(GModelSky* model);
-    void      add_ptsrc_model(GModelSky* model);
-    void      get_bounding_circle(GModelSky* model,
-                                  GSkyDir*   dir,
-                                  double*    radius) const;
+    void init_members(void);
+    void copy_members(const ctmapcube& app);
+    void free_members(void);
+    void get_parameters(void);
+    void create_cube(void);
+    void add_model(GModelSky* model);
+    void add_ptsrc_model(GModelSky* model);
+    void get_bounding_circle(GModelSky* model,
+                             GSkyDir*   dir,
+                             double*    radius) const;
  
     // User parameters
     GFilename m_outcube;                //!< Output map cube filename
@@ -503,15 +503,19 @@ void ctool::setup_models(GObservations&amp;     obs,
  * @brief Create energy boundaries from user parameters
  *
  * @exception GException::invalid_value
- *            No valid energy boundary extension found.
+ *            No valid energy boundary extension found or invalid extension
+ *            name.
  *
  * Get the energy boundaries according to the user parameters. The method
  * supports loading of energy boundary information from the `EBOUNDS` or
- * `ENERGYBINS` extension, or setting energy boundaries using a linear
- * or logarithmical spacing.
+ * `ENERGYBINS` extension, loading of energy value information from the
+ * `ENERGIES` extension, or setting energy boundaries using a linear or
+ * logarithmical spacing.
  *
  * The following parameters are read:
+ *
  *      ebinalg - Energy binning algorithm
+ *      ebinfile - Name of file with energy binning
  *      emin - Minimum energy (if ebinalg != FILE)
  *      emax - Maximum energy (if ebinalg != FILE)
  *      enumbins - Number of energy bins (if ebinalg != FILE)
@@ -529,43 +533,72 @@ GEbounds ctool::create_ebounds(void)
     if (ebinalg == "FILE") {
  
         // Get filename
-        std::string ebinfile = (*this)["ebinfile"].filename();
-
-        // Create filename instance
-        GFilename fname(ebinfile);
+        GFilename ebinfile = (*this)["ebinfile"].filename();
  
-        // Check if extension name was provided
-        if (!fname.has_extname()) {
+        // If no extension name was provided then use default extension names
+        if (!ebinfile.has_extname()) {
  
             // Open energy boundary file using the EBOUNDS or ENERGYBINS
-            // extension. Throw an exception if opening fails.
-            GFits file(ebinfile);
+            // extension or energy values using the ENERGIES extension.
+            // Throw an exception if opening fails.
+            GFits file(ebinfile.url());
             if (file.contains(gammalib::extname_ebounds)) {
                 file.close();
                 ebounds.load(ebinfile);
             }
             else if (file.contains("ENERGYBINS")) {
                 file.close();
-                ebinfile += "[ENERGYBINS]";
+                ebinfile = ebinfile.url() + "[ENERGYBINS]";
                 ebounds.load(ebinfile);
             }
+            else if (file.contains(gammalib::extname_energies)) {
+                file.close();
+                ebounds.set(GEnergies(ebinfile));
+            }
             else {
                 file.close();
                 std::string msg = "No extension with name \""+
-                                  gammalib::extname_ebounds+"\" or"
-                                  " \"ENERGYBINS\" found in FITS file"
-                                  " \""+ebinfile+"\". An \""+
-                                  gammalib::extname_ebounds+"\" or "
-                                  "\"ENERGYBINS\" extension is required if "
-                                  "the parameter \"ebinalg\" is set to "
-                                  "\"FILE\".";
+                                  gammalib::extname_ebounds+"\", "
+                                  "\"ENERGYBINS\" or \""+
+                                  gammalib::extname_energies+"\" found in FITS "
+                                  "file \""+ebinfile+"\". Please specify a "
+                                  "valid energy binning file.";
                 throw GException::invalid_value(G_CREATE_EBOUNDS, msg);
             }
         }
+
+        // ... otherwise load energy boundaries from filename including
+        // extension
         else {
-            // Load ebounds from filename including extension
-            ebounds.load(ebinfile);
-        }
+
+            // Open energy boundary file
+            GFits file(ebinfile.url());
+
+            // If FITS file does not contain requested extension then throw
+            // an exception
+            if (!file.contains(ebinfile.extname())) {
+                std::string msg = "No extension \""+ebinfile.extname()+"\" "
+                                  "found in energy binning file \""+
+                                  ebinfile.url()+"\". Please provide a valid "
+                                  "extension name.";
+                throw GException::invalid_value(G_CREATE_EBOUNDS, msg);
+            }
+
+            // Get table extension
+            GFitsTable& table = *file.table(ebinfile.extname());
+
+            // If table contains one column then load energies, otherwise load
+            // energy boundaries.
+            if (table.ncols() == 1) {
+                file.close();
+                ebounds.set(GEnergies(ebinfile));
+            }
+            else {
+                file.close();
+                ebounds.load(ebinfile);
+            }
+
+        } // endelse: loaded energy boundaries from table extension
  
     } // endif: ebinalg was "FILE"
  
@@ -238,7 +238,13 @@ class Test(test):
                                        gammalib.GEnergy(10.0,'TeV'))
         ebounds.save('test_ebinfile_ebounds.fits[EBOUNDS]', True)
         ebounds.save('test_ebinfile_energybins.fits[ENERGYBINS]', True)
-        ebounds.save('test_ebinfile_energies.fits[ENERGIES]', True)
+        ebounds.save('test_ebinfile_bins.fits[BINS]', True)
+
+        # Create energy bin file
+        energies = gammalib.GEnergies(2, gammalib.GEnergy(1.0,'TeV'),
+                                         gammalib.GEnergy(10.0,'TeV'))
+        energies.save('test_ebinfile_energies.fits[ENERGIES]', True)
+        energies.save('test_ebinfile_engs.fits[ENGS]', True)
  
         # Return
         return
@@ -335,16 +341,28 @@ class Test(test):
         self.test_value(ebounds.emin().TeV(), 1.0, 'Check minimum energy')
         self.test_value(ebounds.emax().TeV(), 10.0, 'Check maximum energy')
  
+        # Test create_ebounds() method for "ENERGIES" extension
+        ebounds = empty.check_create_ebounds('test_ebinfile_energies.fits')
+        self.test_value(ebounds.size(), 1, 'Check number of energy boundaries')
+        self.test_value(ebounds.emin().TeV(), 1.0, 'Check minimum energy')
+        self.test_value(ebounds.emax().TeV(), 10.0, 'Check maximum energy')
+
         # Test create_ebounds() method for invalid extension
         self.test_try('Check create_ebounds() method for invalid extension')
         try:
-            empty.check_create_ebounds('test_ebinfile_energies.fits')
+            empty.check_create_ebounds('test_ebinfile_energies.fits[INVALID]')
             self.test_try_failure('Exception not thrown')
         except ValueError:
             self.test_try_success()
  
-        # Test create_ebounds() method for "ENERGIES" extension
-        ebounds = empty.check_create_ebounds('test_ebinfile_energies.fits[ENERGIES]')
+        # Test create_ebounds() method for "BINS" extension
+        ebounds = empty.check_create_ebounds('test_ebinfile_bins.fits[BINS]')
+        self.test_value(ebounds.size(), 1, 'Check number of energy boundaries')
+        self.test_value(ebounds.emin().TeV(), 1.0, 'Check minimum energy')
+        self.test_value(ebounds.emax().TeV(), 10.0, 'Check maximum energy')
+
+        # Test create_ebounds() method for "ENGS" extension
+        ebounds = empty.check_create_ebounds('test_ebinfile_engs.fits[ENGS]')
         self.test_value(ebounds.size(), 1, 'Check number of energy boundaries')
         self.test_value(ebounds.emin().TeV(), 1.0, 'Check minimum energy')
         self.test_value(ebounds.emax().TeV(), 10.0, 'Check maximum energy')
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the ctmapcube tool.
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -242,7 +242,7 @@ class Test(test):
         return
  
     # Check map cube
-    def _check_cube(self, cube, nmaps=10, npixels=10000):
+    def _check_cube(self, cube, nmaps=11, npixels=10000):
         """
         Check map cube