Merge branch '1335-write-csresmap' into integration

Conflicts: scripts/csresmap.py

Merge branch '1335-write-csresmap' into integration
Conflicts: scripts/csresmap.py
Jürgen Knödlseder
2 parents 6f3d5e5a 561d0ed4
Showing 7 changed files with 182 additions and 90 deletions
ChangeLog
NEWS
pyext/ctmodel.i
scripts/csresmap.py
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
test/test_cscripts.sh.in
-2014-10-10  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-10-15  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.8.0 released
         ========================
         Clean-up public interface
-	Adapt to GammaLib 0.9.0 interface
-        Introduce ctool base class
-        ctbin now always merges event lists into a single counts cube
-        ctmodel now always creates a single model cube
+        Adapt to GammaLib 0.9.0 interface
         Add ctexpcube tool (generates exposure cube)
         Add ctpsfcube tool (generates PSF cube)
         Add ctbkgcube tool (generates background cube)
         Add ctcubemask tool (masks event bins in event cube)
         Add cttsmap tool (generates TS map)
+        Add csresmap script (generates residual map)
+        Introduce ctool base class
+        ctbin now always merges event lists into a single counts cube
+        ctmodel now always creates a single model cube
         Support energy dispersion in ctobssim, ctlike and ctmodel	
         Correct ctobssim for thread safety
         Correctly apply time references in ctselect
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
 Author(s) : Juergen Knoedlseder
-10 October 2014
+15 October 2014
 1. Introduction
@@ -192,6 +192,9 @@ The following functions have been added to obsutils:
 * spectrum() generates spectral points by fitting a source model in
   specified energy bins.
+Add the csresmap script that uses ctbin and ctmodel to generate a
+residual counts cube for a given model.
+
 18. caldb
 ---------
@@ -49,6 +49,7 @@ public:
     void                 save(void);
     const GObservations& obs(void) const;
     const GCTAEventCube& cube(void) const;
+    void                 cube(const GCTAEventCube& cube);
 };
 #! /usr/bin/env python
 # ==========================================================================
-# This script generates a residual map.
+# Residual map generation script.
 #
 # Copyright (C) 2014 Michael Mayer
 #
@@ -18,18 +18,17 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-from ctools import *
-from gammalib import *
+import ctools
+import gammalib
 import sys
-import csv
-import math
 import tempfile
+import os
 # ============== #
 # cstsdist class #
 # ============== #
-class csresmap(GApplication):
+class csresmap(gammalib.GApplication):
     """
     This class implements the creation of a residual map. It derives
     from the GammaLib::GApplication class which provides support for parameter
@@ -45,14 +44,17 @@ class csresmap(GApplication):
         self.version = "0.1.0"
         # Initialise some members
-        self.obs        = None 
-        self.algorithm = "SUB" 
+        self.obs       = None 
+        self.algorithm = "SUB"
+        self.resmap    = None
         # Initialise some members
-        if isinstance(argv[0],GObservations):
+        if isinstance(argv[0],gammalib.GObservations):
             self.obs = argv[0]
+            argv = argv[1:]
         else:      
-            self.obs      = GObservations()   
+            self.obs      = gammalib.GObservations()
+            self.obs.clear()   
         self.outfile = ""
         # Make sure that parfile exists
@@ -60,9 +62,9 @@ class csresmap(GApplication):
         # Initialise application
         if len(argv) == 0:
-            GApplication.__init__(self, self.name, self.version)
+            gammalib.GApplication.__init__(self, self.name, self.version)
         elif len(argv) ==1:
-            GApplication.__init__(self, self.name, self.version, argv[1:])
+            gammalib.GApplication.__init__(self, self.name, self.version, *argv)
         else:
             raise TypeError("Invalid number of arguments given.")
@@ -93,32 +95,32 @@ class csresmap(GApplication):
         parfile = self.name+".par"
         try:
-            pars = GApplicationPars(parfile)
+            pars = gammalib.GApplicationPars(parfile)
         except:
             # Signal if parfile was not found
             sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
             # Create default parfile
-            pars = GApplicationPars()
-            pars.append(GApplicationPar("infile","f","h","obs.xml","","","Observation definition"))
-            pars.append(GApplicationPar("srcmdl","f","h","crab.xml","","","model file name"))
-            pars.append(GApplicationPar("outfile","f","h","resmap.fits","","","Output file name"))
-            pars.append(GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name"))
-            pars.append(GApplicationPar("modfile","f","h","model.fits","","","Output model file name"))
-            pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
-            pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
-            pars.append(GApplicationPar("emin","r","h","0.1","0.0","","Lower energy limit (TeV)"))
-            pars.append(GApplicationPar("emax","r","h","100.0","0.0","","Upper energy limit (TeV)"))
-            pars.append(GApplicationPar("enumbins","i","h","20","","","Number of energy bins"))
-            pars.append(GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
-            pars.append(GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
-            pars.append(GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|STG|TAN","","Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"))
-            pars.append(GApplicationPar("xref","r","a","83.6331","","","First coordinate of image center in degrees (RA or galactic l)"))
-            pars.append(GApplicationPar("yref","r","a","22.01","","","Second coordinate of image center in degrees (DEC or galactic b)"))
-            pars.append(GApplicationPar("nxpix","i","h","200","","","Size of the X axis in pixels"))
-            pars.append(GApplicationPar("nypix","i","h","200","","","Size of the Y axis in pixels"))
-            pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size (deg/pixel)"))
-            pars.append(GApplicationPar("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residualmap algorithm"))
+            pars = gammalib.GApplicationPars()
+            pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition"))
+            pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name"))
+            pars.append(gammalib.GApplicationPar("outfile","f","a","resmap.fits","","","Output file name"))
+            pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name"))
+            pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name"))
+            pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
+            pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
+            pars.append(gammalib.GApplicationPar("emin","r","h","0.1","0.0","","Lower energy limit (TeV)"))
+            pars.append(gammalib.GApplicationPar("emax","r","h","100.0","0.0","","Upper energy limit (TeV)"))
+            pars.append(gammalib.GApplicationPar("enumbins","i","h","20","","","Number of energy bins"))
+            pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
+            pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
+            pars.append(gammalib.GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|STG|TAN","","Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"))
+            pars.append(gammalib.GApplicationPar("xref","r","a","83.6331","","","First coordinate of image center in degrees (RA or galactic l)"))
+            pars.append(gammalib.GApplicationPar("yref","r","a","22.01","","","Second coordinate of image center in degrees (DEC or galactic b)"))
+            pars.append(gammalib.GApplicationPar("nxpix","i","h","200","","","Size of the X axis in pixels"))
+            pars.append(gammalib.GApplicationPar("nypix","i","h","200","","","Size of the Y axis in pixels"))
+            pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size (deg/pixel)"))
+            pars.append(gammalib.GApplicationPar("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residualmap algorithm"))
             pars.append_standard()
             pars.save(parfile)
@@ -133,9 +135,27 @@ class csresmap(GApplication):
         # Set observation if not done before
         if self.obs.size() == 0:
-            self.obs = GObservations(self["infile"].filename)
-            self.m_irf       = self["irf"].string()
-            self.m_caldb     = self["caldb"].string()
+            
+            obsfile = self["infile"].filename()
+            try: 
+                
+                self.obs = gammalib.GObservations(obsfile)
+            except:
+                self.obs.clear()
+                observation = gammalib.GCTAObservation()
+                observation.load(obsfile)
+                
+                self.m_irf   = self["irf"].string()
+                self.m_caldb = self["caldb"].string()
+                
+                caldb = gammalib.GCaldb()
+                if os.path.isdir(self.m_caldb):
+                    caldb.rootdir(self.m_caldb)    
+                else:
+                    caldb.open("cta",self. m_caldb)
+
+                observation.response(self.m_irf, caldb);
+                self.obs.append(observation)
         # Check for models in the container
         # read models from file if there were none
@@ -185,6 +205,9 @@ class csresmap(GApplication):
         """
         # Run the script
         self.run()
+
+        # Save residual map
+        self.resmap.save(self.m_outfile, self.m_clobber)
         # Return
         return
@@ -218,7 +241,8 @@ class csresmap(GApplication):
             self.log.header1("Generate binned map")
         # Create countsmap
-        bin = ctbin(self.obs)
+        bin = ctools.ctbin(self.obs)
+        bin["debug"].boolean(True)
         bin["nxpix"].integer(self.m_nxpix)
         bin["nypix"].integer(self.m_nypix)
         bin["proj"].string(self.m_proj)
@@ -231,75 +255,59 @@ class csresmap(GApplication):
         bin["emax"].real(self.m_emax)
         bin["binsz"].real(self.m_binsz)
         bin["clobber"].boolean(self.m_clobber)
-        if not self.m_binfile == "NONE":
-            bin["outfile"].filename(self.m_binfile)
-        else:
-            os,filename = tempfile.mkstemp()
-            bin["outfile"].filename(filename)
-        bin.execute()
+        bin.run()
+
+        # Store counts map as residual map. Note that we need a
+        # special construct here to avoid memory leaks. This seems
+        # to be a SWIG feature as SWIG creates a new object when
+        # calling bin.cube()
+        #residualmap = bin.cube().map()
+        counts_cube = bin.cube()
+        self.resmap = counts_cube.map().copy()
+        self.resmap.stack_maps()
         # Write header
         if self.logTerse():
             self.log("\n")
             self.log.header1("Generate model map")
-        # create model map
-        model = ctmodel(self.obs)
-        model["infile"].filename(bin["outfile"].filename())
+        # Create model map
+        model = ctools.ctmodel(self.obs)
+        model["debug"].boolean(True)
+        model.cube(bin.cube())
         model["clobber"].boolean(self.m_clobber)
-        if not self.m_modfile == "NONE":
-            model["outfile"].filename(self.m_modfile)
-        else:
-            os,filename = tempfile.mkstemp()
-            model["outfile"].filename(filename)
-        model.execute()
+        model.run()
         # Get model map into GSkymap object
-        modelmap = GSkymap(model["outfile"].filename())
+        modelmap = model.cube().map()
         modelmap.stack_maps()
-
-        # Initialise resiudalmap from countmap
-        residualmap = GSkymap(bin["outfile"].filename())
-        residualmap.stack_maps()
-       
-        
-        
-        # continue calculations depending on given algorithm
+        # Continue calculations depending on given algorithm
         if self.m_algorithm == "SUB":
             # Subtract maps 
-            residualmap -= modelmap
+            self.resmap -= modelmap
         elif self.m_algorithm == "SUBDIV":
             # Subtract and divide by model map
-            residualmap -= modelmap
-            residualmap /= modelmap
+            self.resmap -= modelmap
+            self.resmap /= modelmap
         elif self.m_algorithm == "SUBDIVSQRT":
             # subtract and divide by sqrt of model map
-            residualmap -= modelmap
-            residualmap /= modelmap.sqrt()
+            self.resmap -= modelmap
+            self.resmap /= modelmap.sqrt()
         else:
-            # Rais error if algorithm is unkown
+            # Raise error if algorithm is unkown
             raise TypeError("Algorithm \""+self.m_algorithm+"\" not known")
-        # Save files
-        residualmap.save(self.m_outfile,self.m_clobber)
-        
-        # Remove binned and model files in case they shouldnt be saved
-        if self.m_binfile == "NONE":
-            os.remove(bin["outfile"].filename())
-        if self.m_modfile == "NONE":
-            os.remove(model["outfile"].filename())
-        
         # Return
         return
-
+    
 # ======================== #
 # Main routine entry point #
@@ -246,8 +246,10 @@ void ctmodel::run(void)
         log.header1("Generate model cube");
     }
-    // Initialise counts cube
-    init_cube();
+    // Initialise counts cube if cube has not yet been set
+    if (!has_cube()) {
+        init_cube();
+    }
     // Loop over all observations in the container
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -320,6 +322,31 @@ void ctmodel::save(void)
 }
+/***********************************************************************//**
+ * @brief Set model cube
+ *
+ * @param[in] cube Model cube.
+ *
+ * Set model cube and set all cube bins to zero.
+ ***************************************************************************/
+void ctmodel::cube(const GCTAEventCube& cube)
+{
+    // Set cube
+    m_cube = cube;
+
+    // Set all cube bins to zero
+    for (int i = 0; i < m_cube.size(); ++i) {
+        m_cube[i]->counts(0.0);
+    }
+
+    // Signal that cube has been set
+    m_has_cube = true;
+
+    // Return
+    return;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -362,6 +389,7 @@ void ctmodel::init_members(void)
     m_cube.clear();
     m_ebounds.clear();
     m_gti.clear();
+    m_has_cube = false;
     // Return
     return;
@@ -406,6 +434,7 @@ void ctmodel::copy_members(const ctmodel&amp; app)
     m_cube       = app.m_cube;
     m_ebounds    = app.m_ebounds;
     m_gti        = app.m_gti;
+    m_has_cube   = app.m_has_cube;
     // Return
     return;
@@ -431,7 +460,9 @@ void ctmodel::free_members(void)
 void ctmodel::get_parameters(void)
 {
     // Read input and optionally output cube filenames
-    m_infile = (*this)["infile"].filename();
+    if (!has_cube()) {
+        m_infile = (*this)["infile"].filename();
+    }
     if (read_ahead()) {
         m_outfile = (*this)["outfile"].filename();
     }
@@ -482,7 +513,8 @@ void ctmodel::get_parameters(void)
     // If no cube file has been specified then read all parameters that
     // are necessary to create the cube from scratch (see method init_cube)
-    if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) {
+    if ((!has_cube()) &&
+        ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == ""))) {
         m_ebinalg  = (*this)["ebinalg"].string();
         if (m_ebinalg == "FILE") {
             m_ebinfile = (*this)["ebinfile"].filename();
@@ -670,6 +702,9 @@ void ctmodel::init_cube(void)
     }
+    // Signal that cube has been set
+    m_has_cube = true;
+
     // Return
     return;
 }
@@ -738,7 +773,7 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
         if (logTerse()) {
             log << gammalib::parformat("Model events in cube");
             log << sum << std::endl;
-            log << gammalib::parformat("Cube bins outside energy range");
+            log << gammalib::parformat("Bins outside energy range");
             log << num_outside_ebds << std::endl;
         }
@@ -70,6 +70,7 @@ public:
     void                 save(void);
     const GObservations& obs(void) const;
     const GCTAEventCube& cube(void) const;
+    void                 cube(const GCTAEventCube& cube);
 protected:
     // Protected methods
@@ -80,6 +81,7 @@ protected:
     void setup_obs(void);
     void init_cube(void);
     void fill_cube(const GCTAObservation* obs);
+    bool has_cube(void) const;
     // User parameters
     std::string m_infile;      //!< Input counts cube
@@ -112,13 +114,14 @@ protected:
     GCTAEventCube m_cube;       //!< Model cube
     GEbounds      m_ebounds;    //!< Energy boundaries
     GGti          m_gti;        //!< Model cube GTIs
+    bool          m_has_cube;   //!< Signal if cube has been set or loaded
 };
 /***********************************************************************//**
  * @brief Return observation container
  *
- * @return Reference to observation container
+ * @return Reference to observation container.
  ***************************************************************************/
 inline
 const GObservations& ctmodel::obs(void) const
@@ -130,7 +133,7 @@ const GObservations&amp; ctmodel::obs(void) const
 /***********************************************************************//**
  * @brief Return model cube
  *
- * @return Reference to model cube
+ * @return Reference to model cube.
  ***************************************************************************/
 inline
 const GCTAEventCube& ctmodel::cube(void) const
@@ -138,4 +141,16 @@ const GCTAEventCube&amp; ctmodel::cube(void) const
     return m_cube;
 }
+
+/***********************************************************************//**
+ * @brief Signal if cube has been set or loaded
+ *
+ * @return True if cube has been set or loaded.
+ ***************************************************************************/
+inline
+bool ctmodel::has_cube(void) const
+{
+    return m_has_cube;
+}
+
 #endif /* CTMODEL_HPP */
@@ -131,3 +131,32 @@ else
   exit 1
 fi
 $ECHO " ok"
+
+#
+# Test csresmap
+# =============
+$ECHO @ECHO_N@ "Test csresmap: "
+chmod +x @abs_top_srcdir@/scripts/csresmap.py
+@abs_top_srcdir@/scripts/csresmap.py \
+  infile="data/crab_events.fits" \
+  outfile="resmap.fits" \
+  srcmdl="data/crab.xml" \
+  caldb="irf" \
+  irf="cta_dummy_irf" \
+  nxpix=50 \
+  nypix=50 \
+  binsz=0.02 \
+  coordsys="CEL" \
+  xref=83.63 \
+  yref=22.01 \
+  proj="CAR" \
+  binsz=0.05 > csresmap.log
+$ECHO @ECHO_N@ "."
+if [ -s "resmap.fits" ]
+then
+  $ECHO @ECHO_N@ "."
+else
+  $ECHO " resmap.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"