Adapt to ctools coding conventions, clean-up and test code, and add a unit test.

Jürgen Knödlseder
1 parent 6124ca81
Showing 11 changed files with 170 additions and 154 deletions
AUTHORS
ChangeLog
NEWS
configure.ac
pyext/ctcubemask.i
src/Makefile.am
src/ctcubemask/ctcubemask.cpp
src/ctcubemask/ctcubemask.hpp
src/ctcubemask/ctcubemask.par
test/data/exclusion.reg
test/test_ctools.sh
@@ -2,3 +2,4 @@ Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
 Michael Mayer (michael.mayer@desy.de)
+Chia-Chun Lu (Chia-Chun.Lu@mpi-hd.mpg.de)
-2014-07-18  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-07-20  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.8.0 released
         ========================
@@ -6,6 +6,7 @@
 	Adapt to GammaLib 0.9.0 interface
         ctbin now always merges event lists into a single counts cube
         Add ctexpcube tool
+        Add ctcubemask tool
         Support energy dispersion in ctobssim, ctlike and ctmodel	
         Correct ctobssim for thread safety
         Enhance documentation
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
 Author(s) : Juergen Knoedlseder
-18 July 2014
+20 July 2014
 1. Introduction
@@ -88,12 +88,19 @@ file. The cube definition may be taken from a counts map or may be
 specified by user parameters.
-11. Examples
+11. ctcubemask - CTA cube filtering
+-----------------------------------
+The first version of the ctcubemask tool is included in the release.
+ctcubemask allows filtering of bins in a counts cube by setting the
+bin content to -1. Filtering can be done by energy, for a circular
+Region of Interest (ROI) or using ds9 exclusion regions.
+
+12. Examples
 ------------
 None
-12. Scripts
+13. Scripts
 -----------
 Add logL and event statistics to cspull output file.
@@ -102,11 +109,11 @@ Convert all tabs into spaces.
 The obsutils script has now been moved into the ctools module namespace.
-13. caldb
+14. caldb
 ---------
 None
-14. models
+15. models
 ----------
 None
@@ -883,6 +883,7 @@ AC_CONFIG_FILES([Makefile
                  src/ctskymap/Makefile
                  src/ctselect/Makefile
                  src/ctexpcube/Makefile
+                 src/ctcubemask/Makefile
                  src/ctools-setup
                  test/Makefile
                  test/test_cscripts.sh
 /***************************************************************************
- *                    ctcubemask - CTA mask cube tool                      *
+ *                   ctcubemask.i - CTA mask cube tool                     *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 Chia-Chun Lu                                        *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctcubemask.i
- * @brief CTA data selection tool Python interface definition
+ * @brief CTA mask cube interface definition
  * @author Chia-Chun Lu
  */
 %{
@@ -49,9 +49,9 @@ public:
     void           execute(void);
     void           run(void);
     void           save(void);
-    GObservations& obs(void) { return m_obs; }
+    GObservations& obs(void);
     void           get_parameters(void);
-    void           apply_mask(GCTAObservation* obs, const std::string& filename);
+    void           apply_mask(GCTAObservation* obs);
 };
@@ -12,7 +12,7 @@ SUBDIRS = ctbin \
           ctselect \
           ctskymap \
           ctexpcube \
-	  ctcubemask
+          ctcubemask
 # Scripts to be distributed
 scripts = ctools-setup \
@@ -31,7 +31,7 @@ libctools_la_LIBADD = ctbin/libctbin.la \
                       ctselect/libctselect.la \
                       ctskymap/libctskymap.la \
                       ctexpcube/libctexpcube.la \
-		      ctcubemask/libctcubemask.la
+                      ctcubemask/libctcubemask.la
 # Distribute scripts
 EXTRA_DIST   = $(scripts)
 /***************************************************************************
- *                    ctfilter - CTA cube filter tool                      *
+ *                  ctcubemask.cpp - CTA cube filter tool                  *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -34,9 +34,8 @@
 /* __ Method name definitions ____________________________________________ */
-#define G_RUN                                               "ctcubemask::run()"
-#define G_APPLY_MASK                "ctcubemask::apply_mask(GCTAObservation*,"\
-                                                              "std::string&)"
+#define G_RUN                                             "ctcubemask::run()"
+#define G_APPLY_MASK               "ctcubemask::apply_mask(GCTAObservation*)"
 /* __ Debug definitions __________________________________________________ */
@@ -96,7 +95,7 @@ ctcubemask::ctcubemask(GObservations obs) : GApplication(CTCUBEMASK_NAME, CTCUBE
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctcubemask::ctcubemask(int argc, char *argv[]) : 
-                   GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv)
+            GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -151,7 +150,7 @@ ctcubemask::~ctcubemask(void)
  *
  * @param[in] app Application.
  ***************************************************************************/
-ctcubemask& ctcubemask::operator= (const ctcubemask& app)
+ctcubemask& ctcubemask::operator=(const ctcubemask& app)
 {
     // Execute only if object is not identical
     if (this != &app) {
@@ -287,33 +286,9 @@ void ctcubemask::run(void)
             // Save event file name (for possible saving)
             m_infiles[i] = obs->eventfile();
-            // Get temporary file name
-            std::string filename = std::tmpnam(NULL);
-
-            // Save observation in temporary file
-            obs->save(filename, true);
-
-            // Log saved FITS file
-            if (logExplicit()) {
-                GFits tmpfile(filename);
-                log << std::endl;
-                log.header1("FITS file content of temporary file");
-                log << tmpfile << std::endl;
-                tmpfile.close();
-            }
-
-            // Check temporary file
-            std::string message = check_infile(filename);
-            if (message.length() > 0) {
-                throw GException::app_error(G_RUN, message);
-            }
-
             // Load observation from temporary file, including event selection
-            apply_mask(obs, filename);
+            apply_mask(obs);
-            // Remove temporary file
-            std::remove(filename.c_str());
-            
         } // endif: had a CTA observation
     } // endfor: looped over all observations
@@ -409,7 +384,7 @@ void ctcubemask::get_parameters(void)
         try {
             // Load event list in CTA observation
-            obs.load_binned(m_infile);
+            obs.load(m_infile);
             // Append CTA observation to container
             m_obs.append(obs);
@@ -455,115 +430,119 @@ void ctcubemask::get_parameters(void)
 /***********************************************************************//**
- * @brief Select events
+ * @brief Apply mask to counts cube
  *
  * @param[in] obs CTA observation.
- * @param[in] filename File name.
  *
- * Select events from a FITS file by making use of the selection possibility
- * of the cfitsio library on loading a file. A selection string is created
- * from the specified criteria that is appended to the filename so that
- * cfitsio will automatically filter the event data. This selection string
- * is then applied when opening the FITS file. The opened FITS file is then
- * saved into a temporary file which is the loaded into the actual CTA
- * observation, overwriting the old CTA observation. The ROI, GTI and EBounds
- * of the CTA event list are then set accordingly to the specified selection.
- * Finally, the temporary file created during this process is removed.
- *
- * Good Time Intervals of the observation will be limited to the time
- * interval [m_tmin, m_tmax]. If m_tmin=m_tmax=0, no time selection is
- * performed.
- *
- * @todo Use INDEF instead of 0.0 for pointing as RA/DEC selection
+ * Masks out event bins by setting their content to -1.
  ***************************************************************************/
-void ctcubemask::apply_mask(GCTAObservation* obs, const std::string& filename)
+void ctcubemask::apply_mask(GCTAObservation* obs)
 {
-  std::cout<< "start appply filter" << std::endl;
-    // Open FITS file
-    GFits file(filename);
-    GCTAEventCube* cube = (GCTAEventCube*) obs->events();
-    GEbounds ebounds = cube->ebounds();
-    GSkymap map = cube->map();
-
-    // Apply energy selection
-    int e_idx1 = 0;
-    int e_idx2 = 0;
-    int n_ebin = ebounds.size();
-
-    // Loop over ebounds to find the first energy band
-    for ( int i = 0 ; i < ebounds.size() ; i++){
-      double emin = ebounds.emin(i).TeV();
-      if ( emin > m_emin ){
-	e_idx1 = i;
-	break;
-      }
-    }
+    // Get pointer to CTA event cube. Continue only if we have a cube
+    GCTAEventCube* cube = dynamic_cast<GCTAEventCube*>(const_cast<GEvents*>(obs->events()));
+    if (cube != NULL) {
+
+        // Extract event cube and energy boundaries
+        GSkymap         map     = cube->map();
+        const GEbounds& ebounds = cube->ebounds();
+
+        // Initialise energy selection
+        int n_ebin = ebounds.size();
+        int e_idx1 = 0;
+        int e_idx2 = n_ebin;
+
+        // Loop over ebounds to find the first valid energy band
+        for (int i = 0; i < n_ebin; ++i) {
+            double emin = ebounds.emin(i).TeV() + 1.0e-6; // Rounding tolerance
+            if (emin >= m_emin) {
+                e_idx1 = i;
+                break;
+            }
+        }
-    std::cout << "finish looping over ebounds to find 1st" << std::endl;
-    std::cout << "e_ix1 " << e_idx1 << std::endl;
+        // Loop over ebounds to find the last valid energy band
+        for (int i = n_ebin-1; i >= 0; i--) {
+            double emax = ebounds.emax(i).TeV() - 1.0e-6; // Rounding tolerance
+            if (emax <= m_emax) {
+                e_idx2 = i;
+                break;
+            }
+        }
-    // Loop over ebounds to find the last energy band
-    for ( int i = ebounds.size()-1 ; i < ebounds.size() ; i--){
-      double emin = ebounds.emin(i).TeV();
-      if ( emin < m_emax ){
-	e_idx2 = i-1;
-	break;
-      }
-    }
-     std::cout << "finish looping over ebounds to find 1st" << std::endl;
-     std::cout << "e_ix2 " << e_idx2 << std::endl;
-
-    // Set all pixels outside the desired energy bands negative
-    int npix = map.npix();
-    for ( int i = 0 ; i < e_idx1 ; i++){
-      for ( int pixel = 0 ; pixel < npix ; pixel++){
-		  map(pixel,i) = -1;
-      }
-    }
-    for ( int i = e_idx2 ; i < ebounds.size() ; i++){
-      for ( int pixel = 0 ; pixel < npix ; pixel++){
-		  map(pixel,i) = -1;
-      }
-    }
-    std::cout << "finish energy band filter" << std::endl;
-
-
-   
-    // Set all pixels inside the desired energy bands 
-    // but outside ROI or inside exlusion regions negative
-    GSkyRegions regs = GSkyRegions(m_regfile);
-	GSkyRegionCircle roi = GSkyRegionCircle(m_ra, m_dec, m_rad);
-    for ( int i = e_idx1 ; i <= e_idx2 ; i++){
-      for ( int pixel = 0 ; pixel < npix ; pixel++){
-		  GSkyDir dir = map.inx2dir(pixel);
-		  if ( roi.contains(dir) == false || 
-			   regs.contains(dir) == true){
-			  map(pixel,i) = -1;
-		  }
-      }
-    }
-    std::cout << "finish roi selection" << std::endl;
-   
-    GCTAEventCube newcube( map, ebounds, obs->events()->gti());
-    obs->events(newcube);
+        // Set all pixels outside the desired energy bands to -1.0
+        int npix = map.npix();
+        for (int i = 0; i < e_idx1; ++i) {
+            for (int pixel = 0; pixel < npix; ++pixel) {
+                map(pixel,i) = -1.0;
+            }
+        }
+        for (int i = e_idx2 + 1; i < n_ebin; ++i) {
+            for (int pixel = 0; pixel < npix; ++pixel) {
+                map(pixel,i) = -1.0;
+            }
+        }
-    // Log selected FITS file
-    if (logExplicit()) {
-        log << std::endl;
-        log.header1("FITS file content after filter");
-        log << file << std::endl;
-    }
+        // Log selected energy band
+        if (logTerse()) {
+            log << gammalib::parformat("Selected energy band");
+            log << ebounds.emin(e_idx1).TeV() << " - ";
+            log << ebounds.emax(e_idx2).TeV() << " TeV";
+            log << std::endl;
+        }
-    // Check if we have an IMAGE HDU
-    if (!file.contains("IMAGE")) {
-        std::string message = "No \"IMAGE\" extension found in FITS file. ";
-        throw GException::app_error(G_APPLY_MASK, message);
-    }
+        // Set all pixels inside selected energy bands but outside ROI
+        // to -1.0
+        GSkyRegionCircle roi(m_ra, m_dec, m_rad);
+        for (int i = e_idx1; i <= e_idx2; ++i) {
+            for (int pixel = 0; pixel < npix; ++pixel) {
+                GSkyDir dir = map.inx2dir(pixel);
+                if (!roi.contains(dir)) {
+                    map(pixel,i) = -1.0;
+                }
+            }
+        }
+
+        // Log selected energy band
+        if (logTerse()) {
+            log << gammalib::parformat("Selected ROI");
+            log << "RA=" << m_ra << " deg, ";
+            log << "DEC=" << m_dec << " deg, ";
+            log << "Radius=" << m_rad << " deg";
+            log << std::endl;
+        }
+
+        // Set all pixels inside selected energy bands and inside exclusion
+        // regions to -1.0
+        if (m_regfile != "NONE") {
+            GSkyRegions regions = GSkyRegions(m_regfile);
+            for (int i = e_idx1; i <= e_idx2; ++i) {
+                for (int pixel = 0; pixel < npix; ++pixel) {
+                    GSkyDir dir = map.inx2dir(pixel);
+                    if (regions.contains(dir)) {
+                        map(pixel,i) = -1.0;
+                    }
+                }
+            }
+            if (logTerse()) {
+                log << gammalib::parformat("Exclusion regions");
+                log << std::endl;
+                log.indent(1);
+                log << regions.print();
+                log << std::endl;
+                log.indent(0);
+            }
+        }
+        else {
+            if (logTerse()) {
+                log << gammalib::parformat("Exclusion regions");
+                log << "None" << std::endl;
+            }
+        }
+        // Put back map into the event cube
+        cube->map(map);
-    // Get CTA event list pointer
-    GCTAEventList* list =
-        static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
+    } // endif: observation contained an event cube
     // Return
     return;
@@ -798,7 +777,7 @@ void ctcubemask::save_xml(void)
  * nothing if the observation pointer is not valid.
  ***************************************************************************/
 void ctcubemask::save_counts_map(const GCTAObservation* obs,
-                            const std::string&     outfile) const
+                                 const std::string&     outfile) const
 {
     // Save only if observation is valid
     if (obs != NULL) {
 /***************************************************************************
- *                    ctcubemask - CTA cube filter tool                      *
+ *                  ctcubemask.hpp - CTA cube filter tool                  *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -35,7 +35,7 @@
 /* __Definitions _________________________________________________________ */
 #define CTCUBEMASK_NAME    "ctcubemask"
-#define CTCUBEMASK_VERSION "00-05-00"
+#define CTCUBEMASK_VERSION "00-01-00"
 /***********************************************************************//**
@@ -53,7 +53,7 @@ public:
     virtual ~ctcubemask(void);
     // Operators
-    ctcubemask& operator= (const ctcubemask& app);
+    ctcubemask& operator=(const ctcubemask& app);
     // Methods
     void           clear(void);
@@ -62,7 +62,7 @@ public:
     void           save(void);
     GObservations& obs(void) { return m_obs; }
     void           get_parameters(void);
-    void           apply_mask(GCTAObservation* obs, const std::string& filename);
+    void           apply_mask(GCTAObservation* obs);
 protected:
     // Protected methods
@@ -4,23 +4,23 @@
 #                                                                           #
 # ------------------------------------------------------------------------- #
 #                                                                           #
-#  File:      ctfilter.par                                                  #
-#  Version:   00-04-00                                                      #
+#  File:      ctcubemask.par                                                #
+#  Version:   00-01-00                                                      #
 #                                                                           #
 #  Author:    Chia-Chun Lu                                                  #
 #             c.lu@mpi-hd.mpg.de                                            #
 #             IRAP                                                          #
 #                                                                           #
-#  Purpose:   Parameter file for the ctselect application                   #
+#  Purpose:   Parameter file for the ctcubemask application                 #
 #                                                                           #
 #############################################################################
 #
 # File information
 #=================
-infile,  f, a, "test/data/crab_events.fits.gz",,,"Input event list or observation definition file"
-regfile, f, a, "exlusion.reg",,,"exclusion region file in ds9 format"
-outfile, f, a, "filtered_cube.fits",,,"Output event list or observation definition file"
-prefix,  s, h, "filtered_",,,"Prefix for counts maps in observation definition file"
+infile,  f, a, "cntmap.fits",,,"Input counts cube or observation definition file"
+regfile, f, a, "NONE",,,"Exclusion region file in ds9 format"
+outfile, f, a, "filtered_cube.fits",,,"Output event cube or observation definition file"
+prefix,  s, h, "filtered_",,,"Prefix for counts cube in observation definition file"
 #
 # Selection parameters
 #=====================
+# Region file format: DS9 version 4.1
+global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1
+fk5
+circle(83.63, 21.5, 0.4)
@@ -41,6 +41,7 @@ ctobssim=../src/ctobssim/ctobssim
 ctselect=../src/ctselect/ctselect
 ctskymap=../src/ctskymap/ctskymap
 ctexpcube=../src/ctexpcube/ctexpcube
+ctcubemask=../src/ctcubemask/ctcubemask
 #
@@ -387,3 +388,25 @@ else
   $ECHO " expcube2.fits file is not valid"
   exit 1
 fi
+
+
+#
+# Test ctcubemask
+# ===============
+$ECHO -n "Test ctcubemask: "
+$ctcubemask infile="data/crab_events.fits" \
+            regfile="data/exclusion.reg" \
+            outfile="filtered_cube.fits" \
+            ra=83.63 \
+            dec=22.01 \
+            rad=2.0 \
+            emin=0.1 \
+            emax=100.0
+$ECHO -n "."
+if [ -s "filtered_cube.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " filtered_cube.fits file is not valid"
+  exit 1
+fi