Enhance ctobssim.

Jürgen Knödlseder
1 parent a03d24de
Showing 5 changed files with 171 additions and 217 deletions
src/ctbin/ctbin.cpp
src/ctobssim/ctobssim.cpp
src/ctobssim/ctobssim.par
test/test_ctbin.py
test/test_ctobssim.py
@@ -70,7 +70,7 @@ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION)
  *
  * param[in] obs Observation container.
  *
- * Constructs a ctbin tool from an observation container.
+ * Constructs ctbin tool from an observation container.
  ***************************************************************************/
 ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION)
 {
@@ -92,7 +92,8 @@ ctbin::ctbin(const GObservations&amp; obs) : ctool(CTBIN_NAME, CTBIN_VERSION)
  * @param[in] argc Number of arguments in command line.
  * @param[in] argv Array of command line arguments.
  *
- * Constructs ctbin tool from command line arguments.
+ * Constructs ctbin tool using command line arguments for user parameter
+ * setting.
  ***************************************************************************/
 ctbin::ctbin(int argc, char *argv[]) : 
        ctool(CTBIN_NAME, CTBIN_VERSION, argc, argv)
@@ -149,7 +150,7 @@ ctbin::~ctbin(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] app Application.
+ * @param[in] app ctbin tool.
  * @return ctbin tool.
  *
  * Assigns ctbin tool.
@@ -240,14 +241,7 @@ void ctbin::run(void)
         // Write header for observation
         if (logTerse()) {
-            std::string header = m_obs[i]->instrument() + " observation";
-            if (!m_obs[i]->name().empty()) {
-                header += " \"" + m_obs[i]->name() + "\"";
-            }
-            if (!m_obs[i]->id().empty()) {
-                header += " (id=" + m_obs[i]->id() +")";
-            }
-            log.header3(header);
+            log.header3(get_obs_header(m_obs[i]));
         }
         // Get CTA observation
@@ -328,8 +322,8 @@ void ctbin::save(void)
             // Log filename
             if (logTerse()) {
-                log << "Save counts cube into file \""+m_outcube+"\".";
-                log << std::endl;
+                log << gammalib::parformat("Counts cube file");
+                log << m_outcube.url() << std::endl;
             }
             // Save cube
@@ -60,6 +60,8 @@ const double g_roi_margin = 0.5;      //!&lt; Simulation radius margin (degrees)
 /***********************************************************************//**
  * @brief Void constructor
+ *
+ * Constructs an empty ctobssim tool.
  ***************************************************************************/
 ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION)
 {
@@ -76,7 +78,7 @@ ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION)
  *
  * @param[in] obs Observation container.
  *
- * Constructs application by initialising from an observation container.
+ * Constructs ctobssim tool from an observation container.
  ***************************************************************************/
 ctobssim::ctobssim(const GObservations& obs) :
           ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION)
@@ -98,7 +100,8 @@ ctobssim::ctobssim(const GObservations&amp; obs) :
  * @param[in] argc Number of arguments in command line.
  * @param[in] argv Array of command line arguments.
  *
- * Constructs application using command line arguments for parameter setting.
+ * Constructs ctobssim tool using command line arguments for user parameter
+ * setting.
  ***************************************************************************/
 ctobssim::ctobssim(int argc, char *argv[]) :
           ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
@@ -115,6 +118,8 @@ ctobssim::ctobssim(int argc, char *argv[]) :
  * @brief Copy constructor
  *
  * @param[in] app Application.
+ *
+ * Constructs ctobssim tool from another ctobssim instance.
  ***************************************************************************/
 ctobssim::ctobssim(const ctobssim& app) : ctool(app)
 {
@@ -131,6 +136,8 @@ ctobssim::ctobssim(const ctobssim&amp; app) : ctool(app)
 /***********************************************************************//**
  * @brief Destructor
+ *
+ * Destructs ctobssim tool.
  ***************************************************************************/
 ctobssim::~ctobssim(void)
 {
@@ -151,8 +158,10 @@ ctobssim::~ctobssim(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] app Application.
- * @return Application.
+ * @param[in] app ctobssim tool.
+ * @return ctobssim tool.
+ *
+ * Assigns ctobssim tool.
  ***************************************************************************/
 ctobssim& ctobssim::operator=(const ctobssim& app)
 {
@@ -185,7 +194,9 @@ ctobssim&amp; ctobssim::operator=(const ctobssim&amp; app)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Clear application
+ * @brief Clear ctobssim tool
+ *
+ * Clears ctobssim tool.
  ***************************************************************************/
 void ctobssim::clear(void)
 {
@@ -205,10 +216,10 @@ void ctobssim::clear(void)
 /***********************************************************************//**
- * @brief Simulate event data
+ * @brief Run the ctobssim tool.
  *
- * This method runs the simulation. Results are not saved by this method.
- * Invoke "save" to save the results.
+ * Gets the user parameters, loops over all CTA observations in the
+ * observation container, and simulate events for each observation. 
  ***************************************************************************/
 void ctobssim::run(void)
 {
@@ -286,24 +297,14 @@ void ctobssim::run(void)
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Observations");
-        }
-        else {
-            log.header1("Observation");
-        }
+        log.header1(gammalib::number("Observation", m_obs.size()));
         log << m_obs << std::endl;
     }
     // Write header
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Simulate observations");
-        }
-        else {
-            log.header1("Simulate observation");
-        }
+        log.header1(gammalib::number("Simulate observation", m_obs.size()));
     }
     // From here on the code can be parallelized if OpenMP support
@@ -334,14 +335,7 @@ void ctobssim::run(void)
             // Write header for observation
             if (logTerse()) {
-                std::string header = m_obs[i]->instrument() + " observation";
-                if (m_obs[i]->name().length() > 1) {
-                    header += " \"" + m_obs[i]->name() + "\"";
-                }
-                if (m_obs[i]->id().length() > 1) {
-                    header += " (id=" + m_obs[i]->id() +")";
-                }
-                wrklog.header3(header);
+                wrklog.header3(get_obs_header(m_obs[i]));
             }
             // Get pointer on CTA observation
@@ -1455,8 +1449,8 @@ void ctobssim::save_fits(void)
         // Log filename
         if (logTerse()) {
-            log << "Save event list into file \""+m_outevents+"\".";
-            log << std::endl;
+            log << gammalib::parformat("Event list file");
+            log << m_outevents << std::endl;
         }
         // Save observation into FITS file
@@ -1483,7 +1477,7 @@ void ctobssim::save_xml(void)
 {
     // Get output filename and prefix
     m_outevents = (*this)["outevents"].filename();
-    m_prefix  = (*this)["prefix"].string();
+    m_prefix    = (*this)["prefix"].string();
     // Issue warning if output filename has no .xml suffix
     std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4));
@@ -1513,15 +1507,16 @@ void ctobssim::save_xml(void)
                 if (obs->events()->size() != 0) {
                     // Set event output file name
-                    std::string outfile = m_prefix + gammalib::str(i) + ".fits";
+                    std::string outfile = m_prefix + gammalib::str(i) +
+                                          ".fits";
                     // Store output file name in observation
                     obs->eventfile(outfile);
                     // Log filename
                     if (logTerse()) {
-                        log << "Save event list into file \""+outfile+"\".";
-                        log << std::endl;
+                        log << gammalib::parformat("Event list file");
+                        log << outfile << std::endl;
                     }
                     // Save observation into FITS file
@@ -7,7 +7,7 @@
 #  File:      ctobssim.par                                                  #
 #  Version:   1.1.0                                                         #
 #                                                                           #
-#  Author:    Jurgen Knodlseder                                             #
+#  Author:    Juergen Knoedlseder                                           #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
 #             IRAP                                                          #
 #                                                                           #
@@ -17,13 +17,13 @@
 #
 # General parameters
 #===================
-inobs,     f, h, "NONE",,, "Input event list or observation definition XML file"
-inmodel,   f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
-caldb,     s, a, "prod2",,, "Calibration database"
-irf,       s, a, "South_0.5h",,, "Instrument response function"
+inobs,     f, h, NONE,,, "Input event list or observation definition XML file"
+inmodel,   f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file"
+caldb,     s, a, prod2,,, "Calibration database"
+irf,       s, a, South_0.5h,,, "Instrument response function"
 edisp,     b, h, no,,, "Apply energy dispersion?"
-outevents, f, a, "events.fits",,, "Output event data file or observation definition XML file"
-prefix,    s, h, "sim_events_",,, "Prefix for event lists in observation definition XML file"
+outevents, f, a, events.fits,,, "Output event data file or observation definition XML file"
+prefix,    s, h, sim_events_,,, "Prefix for event lists in observation definition XML file"
 #
 # Simulation parameters
@@ -38,7 +38,7 @@ emin,    r, a, 0.1,,, &quot;Lower energy limit (TeV)&quot;
 emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
 deadc,   r, h, 0.95,0,1, "Average deadtime correction factor"
 maxrate, r, h, 1.0e6,,, "Maximum photon rate"
-eslices, i, h, 10,1,100,"Number of energy slices"
+eslices, i, h, 10,1,100, "Number of energy slices"
 #
 # Standard parameters
@@ -46,5 +46,5 @@ eslices, i, h, 10,1,100,&quot;Number of energy slices&quot;
 chatter, i, h, 2,0,4, "Chattiness of output"
 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
 debug,   b, h, no,,, "Debugging mode activated"
-mode,    s, h, "ql",,, "Mode of automatic parameters"
-logfile, f, h, "ctobssim.log",,, "Log filename"
+mode,    s, h, ql,,, "Mode of automatic parameters"
+logfile, f, h, ctobssim.log,,, "Log filename"
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ========================= #
-# Test class for ctbin tool #
-# ========================= #
+# =========================== #
+# Test class for ctbin module #
+# =========================== #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctbin module.
@@ -39,7 +39,7 @@ class Test(gammalib.GPythonTestSuite):
         gammalib.GPythonTestSuite.__init__(self)
         # Set members
-        self.events_name = "data/crab_events.fits"
+        self._events_name = "data/crab_events.fits"
         # Return
         return
@@ -65,7 +65,7 @@ class Test(gammalib.GPythonTestSuite):
         """
         # Set-up ctbin
         bin = ctools.ctbin()
-        bin["inobs"]    = self.events_name
+        bin["inobs"]    = self._events_name
         bin["outcube"]  = "cntmap.fits"
         bin["ebinalg"]  = "LOG"
         bin["emin"]     = 0.1
@@ -111,20 +111,14 @@ class Test(gammalib.GPythonTestSuite):
 #Please specify a valid parameter name.
         # Save counts cube
-        self.test_try("Save counts cube")
-        try:
-            bin.save()
-            self.test_try_success()
-        except:
-            msg = "Exception occured in saving counts cube."
-            self.test_try_failure(msg)
+        bin.save()
         # Load counts cube and check content.
         evt = gammalib.GCTAEventCube("cntmap.fits")
         self._test_cube(evt, 5542)
         # Prepare observation container for stacked analysis
-        cta = gammalib.GCTAObservation(self.events_name)
+        cta = gammalib.GCTAObservation(self._events_name)
         obs = gammalib.GObservations()
         cta.id("0001")
         obs.append(cta)
 #! /usr/bin/env python
 # ==========================================================================
-# This scripts performs unit tests for the ctobssim tool.
+# This scripts performs unit tests for the ctobssim module.
 #
 # Copyright (C) 2014-2016 Juergen Knoedlseder
 #
@@ -20,15 +20,17 @@
 # ==========================================================================
 import gammalib
 import ctools
+from cscripts import obsutils
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# ============================== #
+# Test class for ctobssim module #
+# ============================== #
 class Test(gammalib.GPythonTestSuite):
     """
-    Test class for ctobssim tool.
+    Test class for ctobssim module.
     """
+
     # Constructor
     def __init__(self):
         """
@@ -38,9 +40,9 @@ class Test(gammalib.GPythonTestSuite):
         gammalib.GPythonTestSuite.__init__(self)
         # Set members
-        self.model_name = "data/crab.xml"
-        self.caldb      = "irf"
-        self.irf        = "cta_dummy_irf"
+        self._model_name = "data/crab.xml"
+        self._caldb      = "irf"
+        self._irf        = "cta_dummy_irf"
         # Return
         return
@@ -54,23 +56,23 @@ class Test(gammalib.GPythonTestSuite):
         self.name("ctobssim")
         # Append tests
-        self.append(self.test_functional, "Test ctobssim functionality")
-        self.append(self.test_container, "Test ctobssim on observation container")
+        self.append(self._test_ctobssim, "Test ctobssim")
+        #self.append(self._test_container, "Test ctobssim on observation container")
         # Return
         return
-    # Test ctobssim functionnality
-    def test_functional(self):
+    # Test ctobssim
+    def _test_ctobssim(self):
         """
-        Test ctobssim functionnality.
+        Test ctobssim.
         """
         # Set-up ctobssim
         sim = ctools.ctobssim()
-        sim["inmodel"]   = self.model_name
+        sim["inmodel"]   = self._model_name
         sim["outevents"] = "events.fits"
-        sim["caldb"]     = self.caldb
-        sim["irf"]       = self.irf
+        sim["caldb"]     = self._caldb
+        sim["irf"]       = self._irf
         sim["ra"]        = 83.63
         sim["dec"]       = 22.01
         sim["rad"]       = 10.0
@@ -80,154 +82,123 @@ class Test(gammalib.GPythonTestSuite):
         sim["emax"]      = 100.0
         # Run tool
-        self.test_try("Run ctobssim")
-        try:
-            sim.run()
-            self.test_try_success()
-        #except Exception as e:
-        #    msg = "Exception occured in ctobssim: %s." % (e,)
-        except:
-            msg = "Exception occured in ctobssim."
-            self.test_try_failure(msg)
+        sim.run()
-        # Retrieve observation and check content
-        obs = gammalib.GCTAObservation(sim.obs()[0])
-        pnt = obs.pointing()
-        evt = obs.events()
-        self.test_assert(obs.instrument() == "CTA", "Observation not a CTA observation")
-        self.test_value(sim.obs().size(), 1, "There is not a single observation")
-        self.test_value(obs.ontime(), 1800.0, 1.0e-6, "Ontime is not 1800 sec")
-        self.test_value(obs.livetime(), 1710.0, 1.0e-6, "Livetime is not 1710 sec")
-        self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "ROI Right Ascension is not 83.63 deg")
-        self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "ROI Declination is not 22.01 deg")
-        self.test_value(evt.size(), 4119, "Number of events is not 4119")
+        # Check content of observation
+        self._test_observation(sim)
+        self._test_list(sim.obs()[0].events(), 4119)
         # Save events
-        self.test_try("Save events")
-        try:
-            sim.save()
-            self.test_try_success()
-        #except Exception as e:
-        #    msg = "Exception occured in saving events: %s." % (e,)
-        except:
-            msg = "Exception occured in saving events."
-            self.test_try_failure(msg)
+        sim.save()
+
+        # Load counts cube and check content.
+        evt = gammalib.GCTAEventList("events.fits")
+        self._test_list(evt, 4119)
-    # Test ctobssim on observation container
-    def test_container(self):
-        """
-        Test ctobssim on observation container.
-        """
         # Set-up observation container
-        obs = self.set_obs(4)
+        pnts = [{'ra': 83.63, 'dec': 21.01},
+                {'ra': 84.63, 'dec': 22.01},
+                {'ra': 83.63, 'dec': 23.01},
+                {'ra': 82.63, 'dec': 22.01}]
+        obs = obsutils.set_obs_list(pnts, caldb=self._caldb, irf=self._irf)
-        # Set-up ctobssim
+        # Set-up ctobssim from observation container
         sim = ctools.ctobssim(obs)
         sim["outevents"] = "sim_events.xml"
-        sim["inmodel"]   = self.model_name
+        sim["inmodel"]   = self._model_name
+
         # Run tool
-        self.test_try("Run ctobssim")
-        try:
-            sim.run()
-            self.test_try_success()
-        #except Exception as e:
-        #    msg = "Exception occured in ctobssim: %s." % (e,)
-        except:
-            msg = "Exception occured in ctobssim."
-            self.test_try_failure(msg)
+        sim.run()
         # Retrieve observation and check content
-        self.test_value(sim.obs().size(), 4, "There are not 4 observations")
+        self._test_observation(sim, nobs=4, pnts=pnts)
+        self._test_list(sim.obs()[0].events(), 4063)
+        self._test_list(sim.obs()[1].events(), 4059)
+        self._test_list(sim.obs()[2].events(), 4010)
+        self._test_list(sim.obs()[3].events(), 4138)
         # Save events
-        self.test_try("Save events")
+        sim.save()
+
+        # Load events
+        obs = gammalib.GObservations("sim_events.xml")
+
+        # Retrieve observation and check content
+        self._test_list(obs[0].events(), 4063)
+        self._test_list(obs[1].events(), 4059)
+        self._test_list(obs[2].events(), 4010)
+        self._test_list(obs[3].events(), 4138)
+
+        # Set-up ctobssim with invalid event file
+        sim = ctools.ctobssim()
+        sim["inmodel"]   = "model_file_that_does_not_exist.xml"
+        sim["outevents"] = "events.fits"
+        sim["caldb"]     = self._caldb
+        sim["irf"]       = self._irf
+        sim["ra"]        = 83.63
+        sim["dec"]       = 22.01
+        sim["rad"]       = 10.0
+        sim["tmin"]      = 0.0
+        sim["tmax"]      = 1800.0
+        sim["emin"]      = 0.1
+        sim["emax"]      = 100.0
+
+        # Run ctbin tool
+        self.test_try("Run ctobssim with invalid model file")
         try:
-            sim.save()
-            self.test_try_success()
-        #except Exception as e:
-        #    msg = "Exception occured in saving events: %s." % (e,)
+            sim.run()
+            self.test_try_failure()
         except:
-            msg = "Exception occured in saving events."
-            self.test_try_failure(msg)
+            self.test_try_success()
+
+        # Return
+        return
-    # Setup observation container
-    def set_obs(self, number):
+    # Check observation
+    def _test_observation(self, ctobssim, nobs=1,
+                          pnts=[{'ra': 83.63, 'dec': 22.01}]):
         """
-        Setup observation container.
+        Test content of an observation.
+        
+        Args:
+            ctobssim: ctobssim instance
+
+        Kwargs:
+            nobs:     Number of observations
+            pnts:     List of pointing dictionaries
         """
-        # Initialise empty observation container
-        obs = gammalib.GObservations()
-
-        # Initialise first time and identifier
-        tstart = 0.0
-        offset = -float(number)/2.0
-
-        # Loop over number of observations
-        for i in range(number):
-            id      = "%6.6d" % i
-            obs_cta = self.set_one_obs(id, offset)
-            offset += 1.0
-            tstart += 1830.0
-            obs.append(obs_cta)
-
-        # Return observation container
-        return obs
-
-    # Setup one observation
-    def set_one_obs(self, id, offset, \
-                    tstart=0.0, duration=1800.0, deadc=0.95, \
-                    emin=0.1, emax=100.0, rad=5.0, \
-                    irf="cta_dummy_irf", caldb="dummy"):
+        # Test observation container
+        self.test_value(ctobssim.obs().size(), nobs,
+             "There is one observation")
+        for i in range(ctobssim.obs().size()):
+            obs = gammalib.GCTAObservation(ctobssim.obs()[i])
+            pnt = obs.pointing()
+            self.test_assert(obs.instrument() == "CTA",
+                 "Observation is CTA observation")
+            self.test_value(obs.ontime(), 1800.0, 1.0e-6,
+                 "Ontime is 1800 sec")
+            self.test_value(obs.livetime(), 1710.0, 1.0e-6,
+                 "Livetime is 1710 sec")
+            self.test_value(pnt.dir().ra_deg(), pnts[i]['ra'], 1.0e-6,
+                 "Pointing Right Ascension is "+str(pnts[i]['ra'])+" deg")
+            self.test_value(pnt.dir().dec_deg(), pnts[i]['dec'], 1.0e-6,
+                 "Pointing Declination is "+str(pnts[i]['dec'])+" deg")
+
+        # Return
+        return
+
+    # Check event list
+    def _test_list(self, list, nevents):
         """
-        Setup one observation for test purposes.
+        Test content of event list."
+        
+        Args:
+            list:    Event list
+            nevents: Expected number of events.
         """
-        # Allocate CTA observation
-        obs_cta = gammalib.GCTAObservation()
-
-        # Set calibration database
-        db = gammalib.GCaldb()
-        if (gammalib.dir_exists(self.caldb)):
-            db.rootdir(self.caldb)
-        else:
-            db.open("cta", self.caldb)
-
-        # Set pointing direction
-        pntdir = gammalib.GSkyDir()
-        pntdir.radec_deg(83.63, 22.01+offset)
-        pnt = gammalib.GCTAPointing()
-        pnt.dir(pntdir)
-        obs_cta.pointing(pnt)
-
-        # Set ROI
-        roi     = gammalib.GCTARoi()
-        instdir = gammalib.GCTAInstDir()
-        instdir.dir(pntdir)
-        roi.centre(instdir)
-        roi.radius(rad)
-
-        # Set GTI
-        gti = gammalib.GGti()
-        gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
-
-        # Set energy boundaries
-        ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
-                                    gammalib.GEnergy(emax, "TeV"))
-
-        # Allocate event list
-        events = gammalib.GCTAEventList()
-        events.roi(roi)
-        events.gti(gti)
-        events.ebounds(ebounds)
-        obs_cta.events(events)
-
-        # Set instrument response
-        obs_cta.response(self.irf, db)
-
-        # Set ontime, livetime, and deadtime correction factor
-        obs_cta.ontime(duration)
-        obs_cta.livetime(duration*deadc)
-        obs_cta.deadc(deadc)
-        obs_cta.id(id)
-
-        # Return CTA observation
-        return obs_cta
+        # Test event list
+        self.test_value(list.size(), nevents, str(nevents)+" elements")
+        self.test_value(list.number(), nevents, str(nevents)+" events")
+
+        # Return
+        return