Add parameter and methods for accessing generated GFits object in ctskymap

This enables accessing the generated skymaps without calling 'save()' first

Add parameter and methods for accessing generated GFits object in ctskymap
This enables accessing the generated skymaps without calling 'save()' first
Josh Cardenzana
1 parent 6642ae65
Showing 4 changed files with 62 additions and 32 deletions
pyext/ctskymap.i
src/ctskymap/ctskymap.cpp
src/ctskymap/ctskymap.hpp
test/test_ctskymap.py
@@ -49,6 +49,7 @@ public:
     void           save(void);
     void           publish(const std::string& name = "");
     const GSkyMap& skymap(void) const;
+    const GFits&   fits(void) const;
     void           exclusion_map(const GSkyRegionMap& exclusion_map);
     GSkyRegionMap  exclusion_map(void) const;
 };
@@ -240,6 +240,9 @@ void ctskymap::run(void)
         publish();
     }
  
+    // Create the output fits object
+    construct_fits();
+
     // Return
     return;
 }
@@ -262,51 +265,60 @@ void ctskymap::save(void)
         // Get sky map filename
         m_outmap = (*this)["outmap"].filename();
  
-        // Create empty FITS file
-        GFits fits;
+        // Save FITS file to disk
+        m_fits.saveto(m_outmap, clobber());
  
-        // Write sky map into FITS file
-        write_map(fits, m_skymap, "SKYMAP");
+    } // endif: filename and map were not empty
  
-        // If background subtraction is requested then write the background map,
-        // the significance map, the counts map and the acceptance map to the
-        // FITS file
-        if (m_bkgsubtract != "NONE") {
+    // Write into logger what has been done
+    std::string fname = (m_outmap.is_empty()) ? "NONE" : m_outmap.url();
+    if (m_skymap.is_empty()) {
+        fname.append(" (map is empty, no file created)");
+    }
+    log_value(NORMAL, "Sky map file", fname);
  
-            // Write background map into FITS file
-            write_map(fits, m_bkgmap, "BACKGROUND");
+    // Return
+    return;
+}
  
  
-            // Write significance map into FITS file
-            write_map(fits, m_sigmap, "SIGNIFICANCE");
+/***********************************************************************//**
+ * @brief Construct GFits object consisting of all maps
+ *
+ * Assembles all of the skymaps into a GFits object. This object is saved
+ * when ctskymap::save() is called. the object can also be accessed by 
+ * calling ctskymap::fits().
+ ***************************************************************************/
+void ctskymap::construct_fits(void)
+{
+    // Write sky map into FITS object
+    write_map(m_fits, m_skymap, "SKYMAP");
  
-            // Write counts map into FITS file
-            write_map(fits, m_counts, "COUNTS");
+    // If background subtraction is requested then write the background map,
+    // the significance map, the counts map and the acceptance map to the
+    // FITS file
+    if (m_bkgsubtract != "NONE") {
  
-            // Write acceptance map into FITS file
-            write_map(fits, m_acceptance, "ACCEPTANCE");
+        // Write background map into FITS file
+        write_map(m_fits, m_bkgmap, "BACKGROUND");
  
-            // If RING background subtraction is requested then write also
-            // the exclusion map to the FITS file
-            if (m_bkgsubtract == "RING") {
-                write_map(fits, m_exclmap, "EXCLUSION");
-            }
+        // Write significance map into FITS file
+        write_map(m_fits, m_sigmap, "SIGNIFICANCE");
  
-        } // endif: background subtraction was requested
+        // Write counts map into FITS file
+        write_map(m_fits, m_counts, "COUNTS");
  
-        // Save FITS file to disk
-        fits.saveto(m_outmap, clobber());
+        // Write acceptance map into FITS file
+        write_map(m_fits, m_acceptance, "ACCEPTANCE");
  
-    } // endif: filename and map were not empty
+        // If RING background subtraction is requested then write also
+        // the exclusion map to the FITS file
+        if (m_bkgsubtract == "RING") {
+            write_map(m_fits, m_exclmap, "EXCLUSION");
+        }
  
-    // Write into logger what has been done
-    std::string fname = (m_outmap.is_empty()) ? "NONE" : m_outmap.url();
-    if (m_skymap.is_empty()) {
-        fname.append(" (map is empty, no file created)");
-    }
-    log_value(NORMAL, "Sky map file", fname);
+    } // endif: background subtraction was requested
  
-    // Return
     return;
 }
  
@@ -372,6 +384,7 @@ void ctskymap::init_members(void)
     m_counts.clear();
     m_acceptance.clear();
     m_exclmap.clear();
+    m_fits.clear();
  
     // Initialise cache
     m_solidangle.clear();
@@ -415,6 +428,7 @@ void ctskymap::copy_members(const ctskymap&amp; app)
     m_counts     = app.m_counts;
     m_acceptance = app.m_acceptance;
     m_exclmap    = app.m_exclmap;
+    m_fits       = app.m_fits;
  
     // Copy cache
     m_solidangle    = app.m_solidangle;
@@ -61,6 +61,7 @@ public:
     void           save(void);
     void           publish(const std::string& name = "");
     const GSkyMap& skymap(void) const;
+    const GFits&   fits(void) const;
     void           exclusion_map(const GSkyRegionMap& exclusion_map);
     GSkyRegionMap  exclusion_map(void) const;
  
@@ -87,6 +88,7 @@ protected:
                                  double&        non,
                                  double&        noff,
                                  double&        alpha);
+    void     construct_fits(void);
     void     ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
                                                   int& iy1, int& iy2) const;
     GSkyMap  ring_convolve(const GSkyMap& map, const double& rmin,
@@ -123,6 +125,7 @@ protected:
     GSkyMap       m_counts;      //!< Counts map
     GSkyMap       m_acceptance;  //!< Acceptance map
     GSkyMap       m_exclmap;     //!< Exclusion map for RING background
+    GFits         m_fits;        //!< Output GFits object
  
     // Caching variables to prevent multiple computations of the same thing
     std::vector<double>  m_solidangle;    //!< Cached pixel solid angles
@@ -145,6 +148,13 @@ const GSkyMap&amp; ctskymap::skymap(void) const
 }
  
  
+inline
+const GFits& ctskymap::fits(void) const
+{
+    return m_fits;
+}
+
+
 /***********************************************************************//**
  * @brief Set exclusion region map
  *
@@ -290,6 +290,11 @@ class Test(test):
         self._check_result_file('ctskymap_ring_py1.fits[BACKGROUND]', 10, 10)
         self._check_result_file('ctskymap_ring_py1.fits[SIGNIFICANCE]', 10, 10)
  
+        # Check generated fits object
+        self._check_map(gammalib.GSkyMap(skymap.fits().image('SKYMAP')), nx=10, ny=10)
+        self._check_map(gammalib.GSkyMap(skymap.fits().image('BACKGROUND')), nx=10, ny=10)
+        self._check_map(gammalib.GSkyMap(skymap.fits().image('SIGNIFICANCE')), nx=10, ny=10)
+
         # Test iterative RING background subtraction
         skymap = ctools.ctskymap()
         skymap['inobs']       = self._events