ctbin::cube() now returns a reference to a member

The building of a cuve on the fly raised issue with memory going out of scope. This is fixed by making the CTA event cube a member of ctbin, and having the cube() method returning a reference to that member

ctbin::cube() now returns a reference to a member
The building of a cuve on the fly raised issue with memory going out of scope. This is fixed by making the CTA event cube a member of ctbin, and having the cube() method returning a reference to that member
Jürgen Knödlseder
1 parent f3f1a514
Showing 16 changed files with 98 additions and 72 deletions
cscripts/csinfo.py
cscripts/csobsdef.py
cscripts/obsutils.py
examples/show_butterfly.py
examples/show_model.py
examples/show_obs.py
examples/show_pha.py
examples/show_pull_evolution.py
examples/show_pull_histogram.py
examples/show_response.py
examples/show_sensitivity.py
examples/show_spectrum.py
examples/show_ts_distribution.py
pyext/ctbin.i
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
 #!/usr/bin/env python
 # ==========================================================================
-# Print info about Gammalib / ctools to the console.
+# Print info about Gammalib / ctools to the console
 #
 # Copyright (C) 2015-2016 Christoph Deil
 #
@@ -87,8 +87,8 @@ def get_command_output(cmd):
     output : str
         Command output
     """
-    # Try to import getoutput from subprocess, alternatively get it from
-    # commands
+    # Try to import getoutput from the "subprocess" module, alternatively get
+    # it from the "commands" module
     try:
         from subprocess import getoutput
     except ImportError:
@@ -287,8 +287,7 @@ def get_pkg_config_info(info, library):
     # Execute command
     out = get_command_output(cmd)
-    # If the command returns 'was not found' then pkg-config is not
-    # available
+    # If the command returns "not found" then pkg-config is not available
     if 'not found' in out:
         out = 'Not available'
@@ -322,8 +321,8 @@ def csinfo(argv):
     """
     Print information about Gammalib and ctools to the console
     """
-    # If there are no command line arguments then print help and exit with
-    # success
+    # If there are no command line arguments then print help text and exit
+    # with success
     if len(argv) <= 1:
         csinfo_print_help()
         sys.exit(0)
@@ -28,7 +28,7 @@ import ctools
 # ============== #
 class csobsdef(ctools.cscript):
     """
-    Creates an observation definition file from a pointing list.
+    Creates an observation definition file from a pointing list
     The pointing list is a comma-separated value ASCII file with header
     keywords in the first row followed by a list of pointings (one
@@ -92,8 +92,7 @@ class csobsdef(ctools.cscript):
         self._pntdef = gammalib.GCsv()
         self._tmin   = 0.0
-        # Initialise application by calling the appropriate class
-        # constructor.
+        # Initialise application by calling the appropriate class constructor
         self._init_cscript(argv)
         # Return
@@ -128,6 +127,20 @@ class csobsdef(ctools.cscript):
         The method creates an XML element so that that the response XML
         writer will write the database and response name into the
         observation definition file.
+
+        Parameters
+        ----------
+        obs : `~gammalib.GCTAObservation`
+            CTA observation
+        caldb : str
+            Calibration database
+        irf : str
+            Instrument response function
+
+        Returns
+        obs : `~gammalib.GCTAObservation`
+            CTA observation with response attached
+        -------
         """
         # Create XML element
         xml = gammalib.GXmlElement()
@@ -156,7 +169,7 @@ class csobsdef(ctools.cscript):
         Raises
         ------
         RuntimeError
-            Invalid pointing definition file format.
+            Invalid pointing definition file format
         """
         # Switch screen logging on in debug mode
         if self._logDebug():
@@ -170,21 +183,25 @@ class csobsdef(ctools.cscript):
             self._log('\n')
             self._log.header1('Creating observation definition XML file')
-        # Load pointing definition file if it is not already set
+        # Load pointing definition file if it is not already set. Extract
+        # the number of columns and pointings
         if self._pntdef.size() == 0:
             self._pntdef = gammalib.GCsv(self['inpnt'].filename(), ',')
         ncols = self._pntdef.ncols()
         npnt  = self._pntdef.nrows()-1
-        # Throw an exception is there is no header information
+        # Raise an exception if there is no header information
         if self._pntdef.nrows() < 1:
             raise RuntimeError('No header found in pointing definition file.')
         # Clear observation container
         self._obs.clear()
+
+        # Initialise observation identifier counter
         identifier = 1
-        # Extract header from pointing definition file
+        # Extract header columns from pointing definition file and put them
+        # into a list
         header = []
         for col in range(ncols):
             header.append(self._pntdef[0,col])
@@ -192,28 +209,31 @@ class csobsdef(ctools.cscript):
         # Loop over all pointings
         for pnt in range(npnt):
-            # Set row index
+            # Set pointing definition CSV file row index
             row = pnt + 1
-            # Create CTA observation
+            # Create empty CTA observation
             obs = gammalib.GCTAObservation()
-            # Set observation name
+            # Set observation name. If no observation name was given then
+            # use "None".
             if 'name' in header:
                 name = self._pntdef[row, header.index('name')]
             else:
                 name = 'None'
             obs.name(name)
-            # Set identifier
+            # Set observation identifier. If no observation identified was
+            # given the use the internal counter.
             if 'id' in header:
-                id_ = self._pntdef[row, header.index('id')]
+                obsid = self._pntdef[row, header.index('id')]
             else:
-                id_ = '%6.6d' % identifier
+                obsid = '%6.6d' % identifier
                 identifier += 1
-            obs.id(id_)
+            obs.id(obsid)
-            # Set pointing
+            # Set pointing. Either use "ra" and "dec" or "lon" and "lat".
+            # If none of these pairs are given then raise an exception.
             if 'ra' in header and 'dec' in header:
                 ra     = float(self._pntdef[row, header.index('ra')])
                 dec    = float(self._pntdef[row, header.index('dec')])
@@ -229,7 +249,8 @@ class csobsdef(ctools.cscript):
                                    'found in pointing definition file.')
             obs.pointing(gammalib.GCTAPointing(pntdir))
-            # Set response function
+            # Set response function. If no "caldb" or "irf" information is
+            # provided then use the user parameter values.
             if 'caldb' in header:
                 caldb = self._pntdef[row, header.index('caldb')]
             else:
@@ -241,14 +262,16 @@ class csobsdef(ctools.cscript):
             if caldb != '' and irf != '':
                 obs = self._set_response(obs, caldb, irf)
-            # Set deadtime correction factor
+            # Set deadtime correction factor. If no information is provided
+            # then use the user parameter value "deadc".
             if 'deadc' in header:
                 deadc = float(self._pntdef[row, header.index('deadc')])
             else:
                 deadc = self['deadc'].real()
             obs.deadc(deadc)
-            # Set Good Time Interval
+            # Set Good Time Interval. If no information is provided then use
+            # the user parameter value "duration".
             if 'duration' in header:
                 duration = float(self._pntdef[row, header.index('duration')])
             else:
@@ -263,7 +286,9 @@ class csobsdef(ctools.cscript):
             obs.ontime(gti.ontime())
             obs.livetime(gti.ontime()*deadc)
-            # Set Energy Boundaries
+            # Set Energy Boundaries. If no "emin" or "emax" information is
+            # provided then use the user parameter values in case they are
+            # valid.
             has_emin = False
             has_emax = False
             if 'emin' in header:
@@ -285,7 +310,8 @@ class csobsdef(ctools.cscript):
                 ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, 'TeV'),
                                             gammalib.GEnergy(emax, 'TeV'))
-            # Set ROI
+            # Set ROI. If no ROI radius is provided then use the user
+            # parameters "rad".
             has_roi = False
             if 'rad' in header:
                 rad     = float(self._pntdef[row, header.index('rad')])
@@ -298,17 +324,17 @@ class csobsdef(ctools.cscript):
                 roi = gammalib.GCTARoi(gammalib.GCTAInstDir(pntdir), rad)
             # Create an empty event list
-            list_ = gammalib.GCTAEventList()
-            list_.gti(gti)
+            event_list = gammalib.GCTAEventList()
+            event_list.gti(gti)
-            # Set optional information
+            # If available, set the energy boundaries and the ROI
             if has_ebounds:
-                list_.ebounds(ebounds)
+                event_list.ebounds(ebounds)
             if has_roi:
-                list_.roi(roi)
+                event_list.roi(roi)
             # Attach event list to CTA observation
-            obs.events(list_)
+            obs.events(event_list)
             # Write observation into logger
             if self._logExplicit():
@@ -136,15 +136,11 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
         # cubes and append them to the observation
         if len(sim.obs()) > 1:
-            # First get the cube. We need this because ctbin builds a cube on
-            # the fly. We should avoid this and make the cube part of ctbin
-            cube = bin.cube()
-
             # Get xref, yref, coordsys and proj from counts cube
-            _xref     = cube.counts().projection().crval(0)
-            _yref     = cube.counts().projection().crval(1)
-            _coordsys = cube.counts().projection().coordsys()
-            _proj     = cube.counts().projection().code()
+            _xref     = bin.cube().counts().projection().crval(0)
+            _yref     = bin.cube().counts().projection().crval(1)
+            _coordsys = bin.cube().counts().projection().coordsys()
+            _proj     = bin.cube().counts().projection().code()
             # Kluge: map "EQU" to "CEL" coordinate system (not sure why the
             # GSkyProjection::coordsys returns "EQU" instead of "CEL"
@@ -25,7 +25,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -25,7 +25,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -25,7 +25,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -26,7 +26,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -25,7 +25,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
 try:
@@ -26,7 +26,7 @@ try:
     import matplotlib.mlab as mlab
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
 try:
@@ -26,7 +26,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -26,7 +26,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -25,7 +25,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
@@ -26,7 +26,7 @@ try:
     import matplotlib.pyplot as plt
     plt.figure()
     plt.close()
-except ImportError, RuntimeError:
+except (ImportError, RuntimeError):
     print('This script needs the "matplotlib" module')
     sys.exit()
 try:
@@ -49,7 +49,7 @@ public:
     void                 save(void);
     void                 publish(const std::string& name = "");
     const GObservations& obs(void) const;
-    GCTAEventCube        cube(void) const;
+    const GCTAEventCube& cube(void) const;
 };
@@ -281,6 +281,10 @@ void ctbin::run(void)
     } // endfor: looped over observations
+    // Build event cube (needs to come before obs_cube() since this method
+    // relies on correct setting of m_cube)
+    m_cube = GCTAEventCube(m_counts, m_weights, m_ebounds, m_gti);
+
     // Set a single cube in the observation container
     obs_cube();
@@ -400,6 +404,7 @@ void ctbin::init_members(void)
     m_weights.clear();
     m_ebounds.clear();
     m_gti.clear();
+    m_cube.clear();
     m_ontime   = 0.0;
     m_livetime = 0.0;
@@ -427,6 +432,7 @@ void ctbin::copy_members(const ctbin&amp; app)
     m_weights   = app.m_weights;
     m_ebounds   = app.m_ebounds;
     m_gti       = app.m_gti;
+    m_cube      = app.m_cube;
     m_ontime    = app.m_ontime;
     m_livetime  = app.m_livetime;
@@ -81,7 +81,7 @@ public:
     void                 save(void);
     void                 publish(const std::string& name = "");
     const GObservations& obs(void) const;
-    GCTAEventCube        cube(void) const;
+    const GCTAEventCube& cube(void) const;
 protected:
     // Protected methods
@@ -95,25 +95,28 @@ protected:
     std::vector<bool>   cube_layer_usage(const GEbounds& ebounds) const;
     // User parameters
-    GFilename     m_outcube;    //!< Output counts map file name
-    bool          m_usepnt;     //!< Use pointing instead of xref/yref parameters
-    bool          m_publish;    //!< Publish counts cube?
-    GChatter      m_chatter;    //!< Chattiness
+    GFilename     m_outcube;  //!< Output counts map file name
+    bool          m_usepnt;   //!< Use pointing instead of xref/yref parameters
+    bool          m_publish;  //!< Publish counts cube?
+    GChatter      m_chatter;  //!< Chattiness
     // Protected members
-    GObservations m_obs;        //!< Observation container
-    GSkyMap       m_counts;     //!< Event cube counts
-    GSkyMap       m_weights;    //!< Event cube weights
-    GEbounds      m_ebounds;    //!< Energy boundaries
-    GGti          m_gti;        //!< Good time intervals
-    double        m_ontime;     //!< Total ontime
-    double        m_livetime;   //!< Total livetime
+    GObservations m_obs;      //!< Observation container
+    GSkyMap       m_counts;   //!< Event cube counts
+    GSkyMap       m_weights;  //!< Event cube weights
+    GEbounds      m_ebounds;  //!< Energy boundaries
+    GGti          m_gti;      //!< Good time intervals
+    GCTAEventCube m_cube;     //!< Events cube (for cube() method)
+    double        m_ontime;   //!< Total ontime
+    double        m_livetime; //!< Total livetime
 };
 /***********************************************************************//**
  * @brief Return observation container
  *
  * @return Reference to observation container
+ *
+ * Returns a reference to an observation container.
  ***************************************************************************/
 inline
 const GObservations& ctbin::obs(void) const
@@ -125,18 +128,14 @@ const GObservations&amp; ctbin::obs(void) const
 /***********************************************************************//**
  * @brief Return event cube
  *
- * @return Event cube
+ * @return Reference to event cube
  *
- * Returns a event cube built from the ctbin member information.
+ * Returns a reference to the event cube.
  ***************************************************************************/
 inline
-GCTAEventCube ctbin::cube(void) const
+const GCTAEventCube& ctbin::cube(void) const
 {
-    // Build event cube
-    GCTAEventCube cube(m_counts, m_weights, m_ebounds, m_gti);
-    
-    // Return cube
-    return cube;
+    return m_cube;
 }
 #endif /* CTBIN_HPP */