Piping for ctskymap implemented.

Andreas Specovius · Jürgen Knödlseder
1 parent 0be0a1a4
Showing 3 changed files with 84 additions and 18 deletions
pyext/ctskymap.i
src/ctskymap/ctskymap.cpp
src/ctskymap/ctskymap.hpp
@@ -48,6 +48,9 @@ public:
     void           run(void);
     void           save(void);
     void           publish(const std::string& name = "");
+
+    void           exclmap(const GSkyMap& exclusionmap);
+    const GSkyMap& exclmap(void) const;
     const GSkyMap& skymap(void) const;
 };
  
@@ -337,6 +337,21 @@ void ctskymap::publish(const std::string&amp; name)
 }
  
  
+/***********************************************************************//**
+ * @brief Set the exclusion map
+ *
+ * @param[in] exclmap Exclusion map object.
+ ***************************************************************************/
+void ctskymap::exclmap(const GSkyMap& exclusionmap)
+{
+    // Assign exclusion map
+    m_exclmap = exclusionmap;
+
+    // Return
+    return;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -522,6 +537,7 @@ void ctskymap::get_parameters(void)
         if (m_iterations > 0) {
             m_threshold  = (*this)["threshold"].real();
         }
+
         (*this)["inexclusion"].query();
         m_usefft      = (*this)["usefft"].boolean();
  
@@ -576,7 +592,12 @@ void ctskymap::setup_maps(void)
     if (m_bkgsubtract != "NONE") {
  
         // Setup the exclusions map
-        setup_exclusion_map();
+        if (m_exclmap.is_empty()) {
+            setup_exclusion_map();
+        }
+        else {
+            adjust_exclusion_map();
+        }
  
         // Cache the pixel solid angles and sky directions
         m_solidangle.reserve(m_counts.npix());
@@ -648,24 +669,11 @@ void ctskymap::setup_exclusion_map(void)
  ***************************************************************************/
 void ctskymap::setup_exclusion_map_fits(const GFilename& filename)
 {
-    // Load the fits image
-    GSkyMap inmap(filename);
-
-    // Loop through the individual pixels in the exclusion map
-    for (int i = 0; i < m_exclmap.npix(); ++i) {
-
-        // Get the pixel direction
-        GSkyDir dir = m_exclmap.pix2dir(i);
+    // Load the fits image to the exclusion map
+    m_exclmap = GSkyMap(filename);
  
-        // Check this sky position in the fits map
-        if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
-
-            // Set the pixel to 1
-            m_exclmap(i) = 1.0;
-
-        } // endif: pixel,region overlap check
-
-    } // endfor: looped over exclusion map pixels
+    // Verify exclusion map
+    adjust_exclusion_map();
  
     // Return
     return;
@@ -704,6 +712,45 @@ void ctskymap::setup_exclusion_map_region(const GFilename&amp; filename)
  
  
 /***********************************************************************//**
+ * @brief Verifys that exclusion map points in fov of counts
+ *
+ * Adjust exclusion map. Reinitialise exclusion map via the counts object 
+ * and transfer the masked areas provided by the current exclusion map to the
+ * updated fov.
+ ***************************************************************************/
+void ctskymap::adjust_exclusion_map(void)
+{
+    // Initialise working map
+    GSkyMap dest_map = m_counts;
+
+    // Set all pixels to 0 (no pixel excluded)
+    dest_map = 0.0;
+
+    // Loop through the individual pixels in the exclusion map
+    for (int i = 0; i < dest_map.npix(); ++i) {
+
+        // Get the pixel direction
+        GSkyDir dir = dest_map.pix2dir(i);
+
+        // Check this sky position in the exclusion map
+        if (m_exclmap.contains(dir) && (m_exclmap(m_exclmap.dir2inx(dir)) != 0.0)) {
+
+            // Set the pixel to 1
+            dest_map(i) = 1.0;
+
+        } // endif: pixel,region overlap check
+
+    } // endfor: looped over exclusion map pixels
+
+    // Set exclusion map
+    m_exclmap = dest_map;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Fill maps from observation container
  *
  * Fill counts and optionally acceptance maps from data in the observation
@@ -60,6 +60,9 @@ public:
     void           run(void);
     void           save(void);
     void           publish(const std::string& name = "");
+
+    void           exclmap(const GSkyMap& exclusionmap);
+    const GSkyMap& exclmap(void) const;
     const GSkyMap& skymap(void) const;
  
 protected:
@@ -72,6 +75,7 @@ protected:
     void     setup_exclusion_map(void);
     void     setup_exclusion_map_fits(const GFilename& filename);
     void     setup_exclusion_map_region(const GFilename& filename);
+    void     adjust_exclusion_map(void);
     void     fill_maps(void);
     void     fill_maps_counts(GCTAObservation* obs);
     void     fill_maps_acceptance(GCTAObservation* obs);
@@ -131,6 +135,18 @@ protected:
  
  
 /***********************************************************************//**
+ * @brief Return exclusion map
+ *
+ * @return Reference to exclusion map
+ ***************************************************************************/
+inline
+const GSkyMap& ctskymap::exclmap(void) const
+{
+    return m_exclmap;
+}
+
+
+/***********************************************************************//**
  * @brief Return observation container
  *
  * @return Reference to observation container