Enhance ctools code coverage (#1828)

Jürgen Knödlseder
1 parent 49530689
Showing 14 changed files with 387 additions and 215 deletions
ChangeLog
NEWS
README.md
pyext/Makefile.am
pyext/ctselect.i
src/ctbkgcube/ctbkgcube.cpp
src/ctselect/ctselect.cpp
src/ctselect/ctselect.hpp
src/cttsmap/cttsmap.cpp
src/ctulimit/ctulimit.cpp
test/data/invalid_event_list.fits.gz
test/test_ctselect.py
test/test_cttsmap.py
test/test_ctulimit.py
-2016-08-06
+2016-08-12
         * Version 1.2.0 released
         ========================
+        Remove ctselect::time_reference() method
         Add ctool and cscript pars() methods (#1838)
         Fix ctbin::save() bug for empty observation container (#1835)
         Fix exception when calling ctools clear() methods (#1835)
 New Features and Important Changes in ctools 1.2.0
-6 August 2016
+12 August 2016
 Introduction
@@ -38,6 +38,9 @@ ctselect - CTA event selection
 The clear() method now calls GApplication::clear() to conserve the properties
 of the tool (#1835).
+The time_reference() method has been removed since the method is part of the
+ctool base class, and exposed in Python as _time_reference().
+
 ctbin - CTA event binning
 -------------------------
 ctools information
 ==================
-* Version:             1.2.0.dev1 (6 August 2016)
+* Version:             1.2.0.dev1 (12 August 2016)
 * Author:              Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
 * GammaLib dependency: 1.2.0.dev1
@@ -64,6 +64,7 @@ testscripts = ctools/tests/__init__.py \
 testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/crab_events.fits.gz \
            $(top_srcdir)/test/data/crab_cntmap.fits.gz \
+           $(top_srcdir)/test/data/invalid_event_list.fits.gz \
            $(top_srcdir)/test/data/exclusion.reg
 # Install test scripts and test data. Note that we don't need to prepend
@@ -44,11 +44,10 @@ public:
     virtual ~ctselect(void);
     // Methods
-    void                  clear(void);
-    void                  run(void);
-    void                  save(void);
-    const GObservations&  obs(void) const;
-    const GTimeReference& time_reference(void) const;
+    void                 clear(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
 };
@@ -246,20 +246,12 @@ void ctbkgcube::run(void)
         // Flag removal
         bool remove = true;
-        // Do we have a CTA, HESS, MAGIC or VERITAS specific model?
-        if (m_bkgmdl[i]->is_valid("CTA", "") ||
-            m_bkgmdl[i]->is_valid("HESS", "") ||
-            m_bkgmdl[i]->is_valid("MAGIC", "") ||
-            m_bkgmdl[i]->is_valid("VERITAS", "")) {
-
-            // Do we have a background model?
-            if (dynamic_cast<GModelData*>(m_bkgmdl[i]) != NULL) {
-
-                // ... then keep model
-                remove = false;
-
-            } // endif: had a background model
-        } // endif: had a CTA, HESS, MAGIC or VERITAS model
+        // If we have a data space model with "GCTA" classname then we
+        // have a CTA background model and we want to keep the model
+        if ((dynamic_cast<GModelData*>(m_bkgmdl[i]) != NULL) &&
+            (m_bkgmdl[i]->classname().substr(0,4) == "GCTA")) {
+            remove = false;
+        }
         // Log model removal or keeping
         std::string what  = (remove) ? "Remove model" : "Keep model";
@@ -291,8 +291,8 @@ void ctselect::run(void)
         // If we have an input file then check it
         if (!m_infiles[i].empty()) {
             std::string message = check_infile(m_infiles[i], m_evtname[i]);
-            if (message.length() > 0) {
-                throw GException::app_error(G_RUN, message);
+            if (!message.empty()) {
+                throw GException::invalid_value(G_RUN, message);
             }
         }
@@ -322,8 +322,8 @@ void ctselect::run(void)
         // If we have a temporary file then check it
         if (!filename.empty()) {
             std::string message = check_infile(filename, m_evtname[i]);
-            if (message.length() > 0) {
-                throw GException::app_error(G_RUN, message);
+            if (!message.empty()) {
+                throw GException::invalid_value(G_RUN, message);
             }
         }
@@ -611,8 +611,9 @@ void ctselect::select_events(GCTAObservation*   obs,
     // Write header into logger
     log_header3(NORMAL, "Events selection");
-    // Initialise selection string
+    // Initialise selection and addition strings
     std::string selection;
+    std::string add;
     // Initialise selection flags
     bool remove_all = false;
@@ -679,6 +680,7 @@ void ctselect::select_events(GCTAObservation*   obs,
         sprintf(cmin, "%.8f", tmin);
         sprintf(cmax, "%.8f", tmax);
         selection = "TIME >= "+std::string(cmin)+" && TIME <= "+std::string(cmax);
+        add       = " && ";
         log_value(NORMAL, "Time range",
                   gammalib::str(tmin)+" - "+gammalib::str(tmax)+" s");
@@ -687,11 +689,6 @@ void ctselect::select_events(GCTAObservation*   obs,
     // Make energy selection
     if (m_select_energy) {
-        // If we have aready a selection than add an "&&" operator
-        if (selection.length() > 0) {
-            selection += " && ";
-        }
-
         // Log the requested energy selection
         if (logNormal()) {
             log << gammalib::parformat("Selected energy range");
@@ -725,43 +722,40 @@ void ctselect::select_events(GCTAObservation*   obs,
                 char cmax[80];
                 sprintf(cmin, "%.8f", emin);
                 sprintf(cmax, "%.8f", emax);
-                selection += "ENERGY >= "+std::string(cmin)+" && ENERGY <= "+
-                             std::string(cmax);
+                selection += add + "ENERGY >= " + std::string(cmin)+
+                                   " && ENERGY <= " + std::string(cmax);
+                add        = " && ";
             }
         }
         else if (emin > 0.0) {
             char cmin[80];
             sprintf(cmin, "%.8f", emin);
-            selection += "ENERGY >= "+std::string(cmin);
+            selection += add + "ENERGY >= " + std::string(cmin);
+            add        = " && ";
         }
         else if (emax > 0.0) {
             char cmax[80];
             sprintf(cmax, "%.8f", emax);
-            selection += "ENERGY <= "+std::string(cmax);
+            selection += add + "ENERGY <= " + std::string(cmax);
+            add        = " && ";
         }
         else {
-            log << "None" << std::endl;
             remove_all = true;
         }
     } // endif: made energy selection
-    // Make ROI selection
+    // Make RoI selection
     if (m_select_roi) {
-        // If we have aready a selection than add an "&&" operator
-        if (selection.length() > 0) {
-            selection += " && ";
-        }
-
-        // Log the requested ROI
+        // Log the requested RoI
         log_value(NORMAL, "Requested RoI",
                   "Centre(RA,DEC)=("+gammalib::str(ra)+", "+
                   gammalib::str(dec)+") deg, Radius="+gammalib::str(m_rad)+
                   " deg");
-        // If we have already an ROI then make sure that the selected
-        // ROI overlaps with the existing ROI
+        // If we have already an RoI then make sure that the selected
+        // ROI overlaps with the existing RoI
         double roi_radius = list->roi().radius();
         if (roi_radius > 0.0) {
             GSkyDir roi_centre = list->roi().centre().dir();
@@ -777,7 +771,7 @@ void ctselect::select_events(GCTAObservation*   obs,
             }
         }
-        // Log the selected ROI
+        // Log the selected RoI
         if (logNormal()) {
             log << gammalib::parformat("Selected RoI");
             if (rad <= 0.0) {
@@ -790,29 +784,25 @@ void ctselect::select_events(GCTAObservation*   obs,
             }
         }
-        // Format ROI selection
+        // Format RoI selection
         char cra[80];
         char cdec[80];
         char crad[80];
         sprintf(cra,  "%.6f", ra);
         sprintf(cdec, "%.6f", dec);
         sprintf(crad, "%.6f", rad);
-        selection += "ANGSEP("+std::string(cra)+"," +
-                     std::string(cdec)+",RA,DEC) <= " +
-                     std::string(crad);
+        selection += add + "ANGSEP("+std::string(cra)+"," +
+                           std::string(cdec)+",RA,DEC) <= " +
+                           std::string(crad);
+        add        = " && ";
     } // endif: made ROI selection
     // Make an expression selection
     if (select_expr) {
-        // If we have aready a selection than add an "&&" operator
-        if (selection.length() > 0) {
-            selection += " && ";
-        }
-
         // Append the expression
-        selection += "("+gammalib::strip_whitespace(m_expr)+")";
+        selection += add + "("+gammalib::strip_whitespace(m_expr)+")";
     } // endif: made expression selection
@@ -875,11 +865,10 @@ void ctselect::select_events(GCTAObservation*   obs,
         GCTAInstDir instdir;
         instdir.dir().radec_deg(ra, dec);
         list->roi(GCTARoi(instdir, rad));
-    } // endif: Roi selection was performed
+    }
+    // ... otherwise restore old RoI information if it existed
     else if (old_roi.is_valid()) {
-        // Restore old RoI information in case no selection was performed
-        // and RoI was existing before
         list->roi(old_roi);
     }
@@ -891,11 +880,10 @@ void ctselect::select_events(GCTAObservation*   obs,
         GEbounds ebounds;
         ebounds.append(GEnergy(emin, "TeV"), GEnergy(emax, "TeV"));
         list->ebounds(ebounds);
-    } //endif: energy selection was performed
+    }
+    // ... otherwise restore old Ebounds if they existed
     else if (old_ebounds.size() > 0) {
-        // Restore old Ebounds in case no energy selection was performed
-        // and observation already had valid Ebounds
         list->ebounds(old_ebounds);
     }
@@ -1342,65 +1330,60 @@ void ctselect::save_event_list(const GCTAObservation* obs,
                                const std::string&     gtiname,
                                const std::string&     outfile) const
 {
-    // Save only if observation is valid
-    if (obs != NULL) {
-
-        // Save only if we have an event list
-        if (obs->eventtype() == "EventList") {
-
-            // Set output FITS file event extension names
-            GFilename   outname(outfile);
-            std::string outevt = evtname;
-            std::string outgti = gtiname;
-            if (outname.has_extname()) {
-                std::vector<std::string> extnames =
-                           gammalib::split(outname.extname(), ";");
-                if (extnames.size() > 0) {
-                    std::string extname = gammalib::strip_whitespace(extnames[0]);
-                    if (!extname.empty()) {
-                        outevt = extname;
-                    }
-                }
-                if (extnames.size() > 1) {
-                    std::string extname = gammalib::strip_whitespace(extnames[1]);
-                    if (!extname.empty()) {
-                        outgti = extname;
-                    }
+    // Save only if we have an event list
+    if (obs->eventtype() == "EventList") {
+
+        // Set output FITS file event extension names
+        GFilename   outname(outfile);
+        std::string outevt = evtname;
+        std::string outgti = gtiname;
+        if (outname.has_extname()) {
+            std::vector<std::string> extnames =
+                       gammalib::split(outname.extname(), ";");
+            if (extnames.size() > 0) {
+                std::string extname = gammalib::strip_whitespace(extnames[0]);
+                if (!extname.empty()) {
+                    outevt = extname;
                 }
             }
-
-            // Create output FITS file
-            GFits outfits;
-
-            // Write observation into FITS file
-            obs->write(outfits, outevt, outgti);
-
-            // Copy all extensions other than evtname and gtiname extensions
-            // from the input to the output event list. The evtname and
-            // gtiname extensions are written by the save method, all others
-            // that may eventually be present have to be copied over
-            // explicitly.
-            GFits infits(infile);
-            for (int extno = 1; extno < infits.size(); ++extno) {
-                GFitsHDU* hdu = infits.at(extno);
-                if (hdu->extname() != evtname &&
-                    hdu->extname() != gtiname &&
-                    hdu->extname() != outevt  &&
-                    hdu->extname() != outgti) {
-                    outfits.append(*hdu);
+            if (extnames.size() > 1) {
+                std::string extname = gammalib::strip_whitespace(extnames[1]);
+                if (!extname.empty()) {
+                    outgti = extname;
                 }
             }
+        }
-            // Close input file
-            infits.close();
+        // Create output FITS file
+        GFits outfits;
+
+        // Write observation into FITS file
+        obs->write(outfits, outevt, outgti);
+
+        // Copy all extensions other than evtname and gtiname extensions
+        // from the input to the output event list. The evtname and
+        // gtiname extensions are written by the save method, all others
+        // that may eventually be present have to be copied over
+        // explicitly.
+        GFits infits(infile);
+        for (int extno = 1; extno < infits.size(); ++extno) {
+            GFitsHDU* hdu = infits.at(extno);
+            if (hdu->extname() != evtname &&
+                hdu->extname() != gtiname &&
+                hdu->extname() != outevt  &&
+                hdu->extname() != outgti) {
+                outfits.append(*hdu);
+            }
+        }
-            // Save file to disk and close it (we need both operations)
-            outfits.saveto(outname.url(), clobber());
-            outfits.close();
+        // Close input file
+        infits.close();
-        } // endif: observation was unbinned
+        // Save file to disk and close it (we need both operations)
+        outfits.saveto(outname.url(), clobber());
+        outfits.close();
-    } // endif: observation was valid
+    } // endif: observation was unbinned
     // Return
     return;
@@ -62,7 +62,6 @@ public:
     void                  run(void);
     void                  save(void);
     const GObservations&  obs(void) const;
-    const GTimeReference& time_reference(void) const;
 protected:
     // Protected methods
@@ -128,16 +127,4 @@ const GObservations&amp; ctselect::obs(void) const
     return m_obs;
 }
-
-/***********************************************************************//**
- * @brief Return CTA time reference
- *
- * @return CTA time reference
- ***************************************************************************/
-inline
-const GTimeReference& ctselect::time_reference(void) const
-{
-    return m_cta_ref;
-}
-
 #endif /* CTSELECT_HPP */
@@ -234,7 +234,7 @@ void cttsmap::run(void)
     // Initialise optimizer
     GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(&log)
-        	                            : new GOptimizerLM();
+                                        : new GOptimizerLM();
     // Store initial models
     GModels models_orig = m_obs.models();
@@ -246,15 +246,15 @@ void cttsmap::run(void)
     models.remove(m_srcname);
     // Get likelihood of null hypothesis if not given before
-	if (m_logL0 == 0.0) {
+    if (m_logL0 == 0.0) {
         // Write header
         log_header1(TERSE, "Compute NULL Hypothesis for TS computation");
-		// Compute likelihood without the test source
-		m_obs.models(models);
-		m_obs.optimize(*opt);
-		m_logL0 = -(opt->value());
+        // Compute likelihood without the test source
+        m_obs.models(models);
+        m_obs.optimize(*opt);
+        m_logL0 = -(opt->value());
 	}
@@ -272,11 +272,9 @@ void cttsmap::run(void)
                     " at "+bincentre.print());
     	// Add test source at current bin position
-        if (m_testsource != NULL) {
-            (*m_testsource)["RA"].value(bincentre.ra_deg());
-            (*m_testsource)["DEC"].value(bincentre.dec_deg());
-            models.append(*m_testsource);
-        }
+        (*m_testsource)["RA"].value(bincentre.ra_deg());
+        (*m_testsource)["DEC"].value(bincentre.dec_deg());
+        models.append(*m_testsource);
     	// Assign models to observations
     	m_obs.models(models);
@@ -292,16 +290,11 @@ void cttsmap::run(void)
     	// Log information
     	if (logExplicit()) {
-    		log << " TS value ....: ";
-            log << ts << std::endl;
+            log_value(EXPLICIT, "TS value", ts);
     	}
-        else if (logTerse()) {
-    		log << "TS for bin number ";
-            log << i;
-            log << " at ";
-            log << bincentre.print();
-            log << ": ";
-            log << ts << std::endl;
+        else if (logNormal()) {
+            log_value(NORMAL, "TS value (bin "+gammalib::str(i)+")",
+                      gammalib::str(ts)+" ("+bincentre.print()+")");
         }
     	// Get test source model instance
@@ -312,11 +305,9 @@ void cttsmap::run(void)
     	m_tsmap(i) = ts;
     	// Extract fitted test source parameters
-        if (m_testsource != NULL) {
-            for (int j = 0; j < m_mapnames.size(); ++j) {
-                m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
-            }
-    	}
+        for (int j = 0; j < m_mapnames.size(); ++j) {
+            m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+        }
     	// Set status of bin to true
     	m_statusmap(i) = 1.0;
@@ -479,7 +470,8 @@ void cttsmap::copy_members(const cttsmap&amp; app)
     m_statusmap = app.m_statusmap;
     // Clone protected members
-    m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL;
+    m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone()
+                                              : NULL;
     // Return
     return;
@@ -491,6 +483,12 @@ void cttsmap::copy_members(const cttsmap&amp; app)
  ***************************************************************************/
 void cttsmap::free_members(void)
 {
+    // Free test source
+    if (m_testsource != NULL) {
+        delete m_testsource;
+        m_testsource = NULL;
+    }
+
     // Return
     return;
 }
@@ -519,26 +517,24 @@ void cttsmap::get_parameters(void)
     // Get name of test source
     m_srcname = (*this)["srcname"].string();
-    // Get Model and store it as protected member
+    // Get model and store it as protected member. Note that the clone
+    // method will never return a NULL pointer, hence we do not need
+    // to check the pointer.
     m_testsource = m_obs.models()[m_srcname]->clone();
     // Check if model has RA and DEC parameters which are necessary to
     // continue
-    if (m_testsource != NULL) {
-        if (!m_testsource->has_par("RA") || !m_testsource->has_par("DEC")) {
-            std::string msg = "Source \""+m_srcname+"\" has no \"RA\" and "
-                              "\"DEC\" parameters. Only sources with \"RA\" "
-                              " and \"DEC\" parameters can be used as test "
-                              "sources.";
-            throw GException::invalid_value(G_GET_PARAMETERS, msg);
-        }
+    if (!m_testsource->has_par("RA") || !m_testsource->has_par("DEC")) {
+        std::string msg = "Source \""+m_srcname+"\" has no \"RA\" and "
+                          "\"DEC\" parameters. Only sources with \"RA\" "
+                          "and \"DEC\" parameters can be used as test "
+                          "sources.";
+        throw GException::invalid_value(G_GET_PARAMETERS, msg);
     }
     // Fix the spatial parameters of the test source
-    if (m_testsource != NULL) {
-        (*m_testsource)["RA"].fix();
-        (*m_testsource)["DEC"].fix();
-    }
+    (*m_testsource)["RA"].fix();
+    (*m_testsource)["DEC"].fix();
     // Create sky map based on task parameters
     GSkyMap map = create_map(m_obs);
@@ -583,32 +579,27 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
     m_tsmap = GSkyMap(map);
     // Initialise map information
-	m_statusmap.clear();
-
-	// Create status map
-	m_statusmap = GSkyMap(map);
-
-	// Initialise maps of free parameters
-    if (m_testsource != NULL) {
+    m_statusmap.clear();
-        // Loop over all model parameters
-        for (int i = 0; i < m_testsource->size(); i++) {
-        
-            // Skip fixed parameters
-            if ((*m_testsource)[i].is_fixed()) {
-                continue;
-            }
+    // Create status map
+    m_statusmap = GSkyMap(map);
-            // Add sky map for free parameters. Note that push_back will
-            // create a copy.
-            m_maps.push_back(map);
+    // Loop over all model parameters
+    for (int i = 0; i < m_testsource->size(); i++) {
+    
+        // Skip fixed parameters
+        if ((*m_testsource)[i].is_fixed()) {
+            continue;
+        }
-            // Store parameter name
-            m_mapnames.push_back((*m_testsource)[i].name());
-        
-        } // endfor: looped over all model parameters
+        // Add sky map for free parameters. Note that push_back will
+        // create a copy.
+        m_maps.push_back(map);
-	} // endif: test source was valid
+        // Store parameter name
+        m_mapnames.push_back((*m_testsource)[i].name());
+    
+    } // endfor: looped over all model parameters
     // Return
     return;
@@ -618,10 +618,10 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
         // Change boundaries for further iteration
         if (eval_mid > 0.0) {
-            wrk_max = mid;
+            wrk_max = mid;     // logL too large, new range = [wrk_min, mid]
         }
-        else if (eval_mid < 0.0) {
-            wrk_min = mid;
+        else {
+            wrk_min = mid;     // logL too small, new range = [mid, wrk_max]
         }
         // Increment counter
@@ -672,7 +672,7 @@ double ctulimit::evaluate(const double&amp; value)
     else {
         std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
                           "outside of validity range requested. To omit "
-                          "this error please enlarge the  parameter range "
+                          "this error please enlarge the parameter range "
                           "in the model XML file.";
         throw GException::invalid_value(G_EVALUATE, msg);
     }
@@ -39,6 +39,9 @@ class Test(test):
         # Call base class constructor
         test.__init__(self)
+        # Set test data
+        self._invalid_events = self._datadir + '/invalid_event_list.fits'
+
         # Return
         return
@@ -104,7 +107,7 @@ class Test(test):
         """
         Test ctselect from Python
         """
-        # Allocate ctselect
+        # Allocate empty ctselect tool
         select = ctools.ctselect()
         # Check that empty ctselect tool holds an empty observation
@@ -122,7 +125,6 @@ class Test(test):
         select.clear()
         # Now set ctselect parameters
-        select = ctools.ctselect()
         select['inobs']   = self._events
         select['ra']      = 83.63
         select['dec']     = 22.01
@@ -130,7 +132,7 @@ class Test(test):
         select['tmin']    = 0
         select['tmax']    = 1800
         select['emin']    = 0.1
-        select['emax']    = 100
+        select['emax']    = 100.0
         select['outobs']  = 'ctselect_py1.fits'
         select['logfile'] = 'ctselect_py1.log'
         select['chatter'] = 2
@@ -166,7 +168,6 @@ class Test(test):
         # Check result file
         self._check_result_file('ctselect_py2.fits')
-        print('Step3')
         # Execute again the copy of ctselect tool again, now using the pointing
         # information in the input observation
@@ -186,19 +187,75 @@ class Test(test):
         # Check that cleared ctselect tool holds no observations and events
         self._check_obs(cpy_select.obs(), nobs=0)
+        # Get mixed observation container
+        obs = self._obs_mixed()
+
+        # Now setup ctselect tool from observation container. We make an
+        # event selection that will result in zero events in the file.
+        select = ctools.ctselect(obs)
+        select['ra']      = 83.63
+        select['dec']     = 22.01
+        select['rad']     = 3.0
+        select['tmin']    = 0
+        select['tmax']    = 1800
+        select['emin']    = 120.0
+        select['emax']    = 130.0
+        select['expr']    = 'DETX == 0'
+        select['outobs']  = 'ctselect_py4.fits'
+        select['logfile'] = 'ctselect_py4.log'
+        select['chatter'] = 3
+
+        # Execute tool
+        select.logFileOpen()  # Needed to get a new log file
+        select.execute()
+
+        # Check result file
+        self._check_result_file('ctselect_py4.fits', nevents=0)
+
+        # Setup ctselect tool for an invalid event file
+        select = ctools.ctselect()
+        select['inobs']   = self._invalid_events
+        select['ra']      = 83.63
+        select['dec']     = 22.01
+        select['rad']     = 3
+        select['tmin']    = 0
+        select['tmax']    = 1800
+        select['emin']    = 0.1
+        select['emax']    = 100.0
+        select['outobs']  = 'ctselect_py5.fits'
+        select['logfile'] = 'ctselect_py5.log'
+        select['chatter'] = 3
+
+        # Execute tool
+        select.logFileOpen()  # Needed to get a new log file
+        self.test_try('Test invalid event file')
+        try:
+            select.execute()
+            self.test_try_failure('Exception not thrown')
+        except ValueError:
+            self.test_try_success()
+
+
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename):
+    def _check_result_file(self, filename, nevents=6127):
         """
         Check result file
+
+        Parameters
+        ----------
+        filename : str
+            Event list file name
+        nevents : int, optional
+            Expected number of events
         """
         # Open result file
         events = gammalib.GCTAEventList(filename)
         # Check event list
-        self._check_events(events, nevents=6127)
+        self._check_events(events, nevents=nevents)
         # Return
         return
 #! /usr/bin/env python
 # ==========================================================================
-# This scripts performs unit tests for the cttsmap tool.
+# This scripts performs unit tests for the cttsmap tool
 #
 # Copyright (C) 2014-2016 Juergen Knoedlseder
 #
@@ -18,6 +18,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import os
 import gammalib
 import ctools
 from testing import test
@@ -105,14 +106,31 @@ class Test(test):
         """
         Test cttsmap from Python
         """
-        # Set-up cttsmap
+        # Allocate empty cttsmap tool
         tsmap = ctools.cttsmap()
+
+        # Check that empty cttsmap tool holds an empty observation and an
+        # empty TS map
+        self._check_obs(tsmap.obs(), nobs=0)
+        self._check_ts_map(tsmap.tsmap(), nx=0, ny=0)
+
+        # Check that saving does not nothing
+        tsmap['outmap']  = 'cttsmap_py0.fits'
+        tsmap['logfile'] = 'cttsmap_py0.log'
+        tsmap.logFileOpen()
+        tsmap.save()
+        self.test_assert(not os.path.isfile('cttsmap_py0.fits'),
+             'Check that no TS map has been created')
+
+        # Check that clearing does not lead to an exception or segfault
+        tsmap.clear()
+
+        # Now set cttsmap parameters
         tsmap['inobs']    = self._events
         tsmap['inmodel']  = self._model
         tsmap['srcname']  = 'Crab'
         tsmap['caldb']    = self._caldb
         tsmap['irf']      = self._irf
-        tsmap['outmap']   = 'cttsmap_py1.fits'
         tsmap['nxpix']    = 3
         tsmap['nypix']    = 3
         tsmap['binsz']    = 0.02
@@ -120,6 +138,7 @@ class Test(test):
         tsmap['proj']     = 'CAR'
         tsmap['xref']     = 83.63
         tsmap['yref']     = 22.01
+        tsmap['outmap']   = 'cttsmap_py1.fits'
         tsmap['logfile']  = 'cttsmap_py1.log'
         tsmap['chatter']  = 2
@@ -131,32 +150,116 @@ class Test(test):
         # Check result file
         self._check_result_file('cttsmap_py1.fits')
+        # Copy cttsmap tool
+        cpy_tsmap = tsmap.copy()
+
+        # Check observation container and TS map of copy
+        self._check_obs(cpy_tsmap.obs())
+        self._check_ts_map(cpy_tsmap.tsmap())
+
+        # Execute copy of cttsmap tool again, now with a higher chatter
+        # level than before
+        cpy_tsmap['binmin']  = 0
+        cpy_tsmap['binmax']  = 1
+        cpy_tsmap['outmap']  = 'cttsmap_py2.fits'
+        cpy_tsmap['logfile'] = 'cttsmap_py2.log'
+        cpy_tsmap['chatter'] = 3
+        cpy_tsmap['publish'] = True
+        cpy_tsmap.logFileOpen()  # Needed to get a new log file
+        cpy_tsmap.execute()
+
+        # Check result file
+        self._check_result_file('cttsmap_py2.fits')
+
+        # Now clear copy of cttsmap tool
+        cpy_tsmap.clear()
+
+        # Check that the cleared copy has also cleared the observation
+        # container and TS map
+        self._check_obs(cpy_tsmap.obs(), nobs=0)
+        self._check_ts_map(cpy_tsmap.tsmap(), nx=0, ny=0)
+
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename):
+    def _check_result_file(self, filename, nx=3, ny=3):
         """
         Check result file
+
+        Parameters
+        ----------
+        filename : str
+            TS map file name
+        nx : int, optional
+            Number of X pixels
+        ny : int, optional
+            Number of Y pixels
         """
         # Open result file
         fits = gammalib.GFits(filename)
         # Get HDUs
-        ts        = fits['Primary']
-        prefactor = fits['Prefactor']
-        index     = fits['Index']
+        ts    = fits['Primary']
+        pre   = fits['Prefactor']
+        index = fits['Index']
         # Check dimensions
-        self.test_value(ts.naxis(), 2, 'Check for 2 dimensions')
-        self.test_value(ts.naxes(0), 3, 'Check for 3 pixels in X')
-        self.test_value(ts.naxes(1), 3, 'Check for 3 pixels in Y')
-        self.test_value(prefactor.naxis(), 2, 'Check for 2 dimensions')
-        self.test_value(prefactor.naxes(0), 3, 'Check for 3 pixels in X')
-        self.test_value(prefactor.naxes(1), 3, 'Check for 3 pixels in Y')
-        self.test_value(index.naxis(), 2, 'Check for 2 dimensions')
-        self.test_value(index.naxes(0), 3, 'Check for 3 pixels in X')
-        self.test_value(index.naxes(1), 3, 'Check for 3 pixels in Y')
+        self.test_value(ts.naxis(), 2, 'Check for TS map dimensions')
+        self.test_value(ts.naxes(0), nx, 'Check for TS map pixels in X')
+        self.test_value(ts.naxes(1), ny, 'Check for TS map pixels in Y')
+        self.test_value(pre.naxis(), 2, 'Check for Prefactor map dimensions')
+        self.test_value(pre.naxes(0), nx, 'Check for Prefactor map pixels in X')
+        self.test_value(pre.naxes(1), ny, 'Check for Prefactor map pixels in Y')
+        self.test_value(index.naxis(), 2, 'Check for Index map dimensions')
+        self.test_value(index.naxes(0), nx, 'Check for Index map pixels in X')
+        self.test_value(index.naxes(1), ny, 'Check for Index map pixels in Y')
+
+        # Return
+        return
+
+    # Check observation
+    def _check_obs(self, obs, nobs=1):
+        """
+        Check observation
+
+        Parameters
+        ----------
+        obs : `~gammalib.GObservations`
+            Models
+        nobs : int, optional
+            Expected number of observations
+        """
+        # Check number of observations
+        self.test_value(obs.size(), nobs, 'Check number of observations')
+
+        # Return
+        return
+
+    # Check TS map
+    def _check_ts_map(self, tsmap, nx=3, ny=3):
+        """
+        Check TS map
+
+        Parameters
+        ----------
+        tsmap : `~gammalib.GSkyMap`
+            TS map
+        nx : int, optional
+            Number of X pixels
+        ny : int, optional
+            Number of Y pixels
+        """
+        # Determine number of maps
+        if nx > 0 and ny > 0:
+            nmaps = 1
+        else:
+            nmaps = 0
+        
+        # Check dimensions
+        self.test_value(tsmap.nmaps(), nmaps, 'Check number of maps')
+        self.test_value(tsmap.nx(), nx, 'Check for number of X pixels')
+        self.test_value(tsmap.ny(), ny, 'Check for number of Y pixels')
         # Return
         return
@@ -118,7 +118,6 @@ class Test(test):
         ulimit.clear()
         # Now set ctulimit parameters
-        ulimit = ctools.ctulimit()
         ulimit['inobs']   = self._events
         ulimit['inmodel'] = self._model
         ulimit['srcname'] = 'Crab'
@@ -140,6 +139,13 @@ class Test(test):
         self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
                         'Check upper limit on energy flux')
+        # Check obs() method
+        self.test_value(ulimit.obs().size(), 1,
+                        'Check number of observations in container')
+
+        # Check opt() method
+        self.test_value(ulimit.opt().status(), 0, 'Check optimizer status')
+
         # Copy ctulimit tool
         cpy_ulimit = ulimit.copy()
@@ -166,6 +172,9 @@ class Test(test):
         self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
                         'Check upper limit on energy flux')
+        # Save copy of observation container for later
+        obs = cpy_ulimit.obs().copy()
+
         # Now clear copy of ctulimit tool
         cpy_ulimit.clear()
@@ -177,6 +186,52 @@ class Test(test):
         self.test_value(cpy_ulimit.eflux_ulimit(), 0.0, 1.0e-16,
                         'Check upper limit on energy flux')
+        # Now set ctulimit tool using the observation container from last
+        # time. This should avoid the necessity to recompute the maximum
+        # likelihood
+        ulimit = ctools.ctulimit(obs)
+        #ulimit['inobs']   = self._events
+        #ulimit['inmodel'] = self._model
+        ulimit['srcname'] = 'Crab'
+        #ulimit['caldb']   = self._caldb
+        #ulimit['irf']     = self._irf
+        ulimit['logfile'] = 'ctulimit_py3.log'
+        ulimit['chatter'] = 4
+
+        # Execute ctulimit tool
+        ulimit.logFileOpen()  # Needed to get a new log file
+        ulimit.execute()
+
+        # Check results
+        self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+                        'Check differential upper limit')
+        self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
+                        'ulimit upper limit on photon flux')
+        self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
+                        'Check upper limit on energy flux')
+
+        # Test invalid model name
+        self.test_try('Test invalid model name')
+        try:
+            ulimit['srcname'] = 'Weihnachtsstern'
+            ulimit['logfile'] = 'ctulimit_py4.log'
+            ulimit.logFileOpen()
+            ulimit.execute()
+            self.test_try_failure('Exception not thrown')
+        except ValueError:
+            self.test_try_success()
+
+        # Test specification of background model
+        ulimit['srcname'] = 'Background'
+        ulimit['logfile'] = 'ctulimit_py5.log'
+        ulimit.logFileOpen()
+        self.test_try('Test invalid model name')
+        try:
+            ulimit.execute()
+            self.test_try_failure('Exception not thrown')
+        except ValueError:
+            self.test_try_success()
+
         # Return
         return