implemented first version of ctulimit

Michael Mayer
1 parent e5a1a87e
Showing 10 changed files with 613 additions and 70 deletions
doc/source/reference_manual/ctulimit.rst
doc/source/reference_manual/index.rst
pyext/ctulimit.i
src/ctulimit/ctulimit.cpp
src/ctulimit/ctulimit.hpp
src/ctulimit/ctulimit.par
src/ctulimit/main.cpp
test/test_ctools.sh
test/test_ctulimit.py
test/test_python.py
+.. _ctulimit:
+
+ctulimit
+===========
+
+Computes upper limit.
+
+
+Synopsis
+--------
+
+Calculates the upper limit flux of a given source and stores the value in ascii file.
+The input model has to be optimised for ctulimit to converge.
+
+General parameters
+------------------
+
+``inobs [file]``
+    Input event list, counts cube or observation definition XML file.
+ 	 	 
+``inmodel [file]``
+    Model XML file containing the optimised source and background definitions.
+ 	 	 
+``srcname [string]``
+    Name of model component for which upper limit should be computed.
+ 	 	 
+``expcube = NONE [file]``
+    Exposure cube file (only needed for stacked analysis).
+
+``psfcube = NONE [file]``
+    PSF cube file (only needed for stacked analysis).
+
+``caldb [string]``
+    Calibration database.
+ 	 	 
+``irf [string]``
+    Instrumental response function.
+ 	 	 
+``outfile [file]``
+    Output ASCII file name.
+
+``(cl = 0.95) [real]``
+    Confidence Level of upper limit
+    
+``(sigma_min = 0.0) [real]``
+    Minimum boundary to start searching for upper limit value. Number of standard deviations above best fit value
+    
+``(sigma_max = 10.0) [real]``
+    Maximum boundary to start searching for upper limit value. Number of standard deviations above best fit value  
+ 	 	 
+``(emin = 1.0) [real])``
+    Minimum energy of flux upper limit (in TeV).
+ 	 	 
+``(emax = 100) [real]``
+    Maximum energy of flux upper limit (in TeV).
+ 	 	 
+``(tol = 1e-5) [real]``
+    Computation tolerance for minimum
+   
+``(max_iter = 50) [integer]``
+    Maximum number of iterations before throwing an exception
+
+
+Standard parameters
+-------------------
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     chatter = 0: no information will be logged
+     
+     chatter = 1: only errors will be logged
+     
+     chatter = 2: errors and actions will be logged
+     
+     chatter = 3: report about the task execution
+     
+     chatter = 4: detailed report about the task execution
+ 	 	 
+``(clobber = yes) [boolean]``
+    Specifies whether an existing output counts cube should be overwritten.
+ 	 	 
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+ 	 	 
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = ctskymap.log) [string]``
+    Name of log file.
+
+
+Related tools
+-------------
+
+None
@@ -27,6 +27,7 @@ ctools
    ctselect --- Selects event data <ctselect>
    ctskymap --- Generates sky map <ctskymap>
    cttsmap --- Generates Test Statistics map <cttsmap>
+   ctulimit --- Calculates upper limit <ctulimit>
  
  
 cscripts
+/***************************************************************************
+ *                  ctulimit - ctulimit calculation tool                   *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2015 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctulimit
+ * @brief upper limit calculation tool
+ * @author Michael Mayer
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctulimit.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctulimit
+ *
+ * @brief upper limit calculation tool
+ ***************************************************************************/
+class ctulimit : public ctool {
+public:
+    // Constructors and destructors
+    ctulimit(void);
+    explicit ctulimit(const GObservations& obs);
+    ctulimit(int argc, char *argv[]);
+    ctulimit(const ctulimit& app);
+    virtual ~ctulimit(void);
+    
+    // Methods
+    void                 clear(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    const double&     diff_ulimit(void) const;
+    const double&     flux_ulimit(void) const;
+    const double&    eflux_ulimit(void) const;
+};
+
+
+/***********************************************************************//**
+ * @brief upper limit calculation tool Python extensions
+ ***************************************************************************/
+%extend ctulimit {
+    ctulimit copy() {
+        return (*self);
+    }
+}
 /***************************************************************************
- *                    ctulimitt - upper limit calculation tool                    *
+ *                    ctulimit - upper limit calculation tool                    *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
@@ -31,10 +31,12 @@
 #include <cstdio>
 #include "ctulimit.hpp"
 #include "GTools.hpp"
+#include "GOptimizer.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_PARAMETERS                          "ctulimit::get_parameters()"
-
+#define G_UL_BISECTION        "ctulimit::ul_bisection(const double&, const double&)"
+#define G_EVALUATE        "ctulimit::evaluate(const double&)"
 /* __ Debug definitions __________________________________________________ */
  
 /* __ Coding definitions _________________________________________________ */
@@ -225,17 +227,200 @@ void ctulimit::run(void)
         log << m_obs << std::endl;
     }
  
+    GModels models_orig = m_obs.models();
+
+
+    m_bestloglike = m_obs.logL();
+    if (m_bestloglike == 0.0) {
+
+        // Write header
+        if (logTerse()) {
+            log << std::endl;
+            log.header1("Recompute best-fit likelihood");
+        }
+
+        // Reoptimize if likelihood was not given before
+        GOptimizerLM* opt = new GOptimizerLM();
+        m_obs.optimize(*opt);
+        m_obs.errors(*opt);
+        m_bestloglike = m_obs.logL();
+
+    } //endif: likelihood was 0.0
+
+    // Store optimized models and pointer to skymodel of interest
+    m_models = m_obs.models();
+    m_skymodel = dynamic_cast<GModelSky*>(m_models[m_srcname])->clone();
+
+    // Get reference energy for differential upper limit
+    GEnergy eref = GEnergy(m_eref, "TeV");
+
+    // Create energy range for flux limits
+    GEnergy emin = GEnergy(m_emin, "TeV");
+    GEnergy emax = GEnergy(m_emax, "TeV");
+
+    // Get value, scale and error of parameter
+    double value = (*m_skymodel->spectral())[0].value();
+    double error = (*m_skymodel->spectral())[0].error();
+    double scale = (*m_skymodel->spectral())[0].scale();
+
+    // Compute starting boundaries
+    double parmin = value / scale +  m_sigma_min * error / scale;
+    double parmax = value / scale +  m_sigma_max * error / scale;
+
     // Write header
     if (logTerse()) {
         log << std::endl;
         log.header1("Compute upper limit");
+        log << "Searching for upper limit between ";
+        log << parmin * scale;
+        log << " and ";
+        log << parmax * scale;
+        log << ": "<<std::endl;
+    }
+
+    // compute upper limit
+    ulimit_bisection(parmin, parmax, scale);
+
+    m_diff_ulimit = m_skymodel->spectral()->eval(eref, GTime());
+    m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax);
+    m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax);
+
+    // Write results to logfile
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Upper limit computation finished");
+        log << "Upper limits are:" << std::endl;
+        log << " Differential at ";
+        log << m_eref <<" TeV: ";
+        log << m_diff_ulimit << " ph/cm2/s/MeV" << std::endl;
+        log << " Integral [";
+        log << m_emin << " - " << m_emax;
+        log << "] TeV: ";
+        log << m_flux_ulimit << " ph/cm2/s" <<std::endl;
+        log << " Energy flux [";
+        log << m_emin << " - " << m_emax;
+        log << "] TeV: ";
+        log << m_eflux_ulimit << " erg/cm2/s" <<std::endl;
+    }
+
+    // Recover original models
+    m_obs.models(models_orig);
+
+    // Return
+    return;
+}
+
+/***********************************************************************//**
+ * @brief Evaluates the likelihood
+ *
+ * @param[in] value Factorised parameter value of the parameter of interest
+ * @return likelihood value
+ *
+ * This method evaluates the likelihood as a function of the parameter of interest
+ ***************************************************************************/
+double ctulimit::evaluate(const double& value)
+{
+    // Initialise likelihood value
+    double LogL = 0.0;
+
+    // Check if given parameter is within boundaries
+    if (value > (*m_skymodel->spectral())[0].factor_min() && value < (*m_skymodel->spectral())[0].factor_max()) {
+
+        // Remove source model
+        m_models.remove(m_srcname);
+
+        // Change spectral model value
+        (*m_skymodel->spectral())[0].factor_value(value);
+
+        // Re-append model to container
+        m_models.append(*m_skymodel);
+
+        // Assign new models to observations
+        m_obs.models(m_models);
+
+        // Evaluate likelihood for new model container
+        m_obs.eval();
+
+        // Retrieve likelihood
+        LogL = m_obs.logL();
+
+    } // endif: value was inside allowed range
+
+    else {
+        // throw exception if value out of range
+        std::string msg = "Value out of range requested: "
+                "To omit this error, you could increase the allowed "
+                "parameter range in the model xml file.";
+        throw GException::invalid_value(G_EVALUATE, msg);
     }
  
+    // Compute function value
+    double difflike = LogL - m_bestloglike - m_dloglike;
  
-    // Compute upper limit here
+    // Return
+    return difflike;
+}
+
+/***********************************************************************//**
+ * @brief Induces upper limit computation by using a bisection method
+ *
+ * This method calculates the upper limit using a bisection method
+ ***************************************************************************/
+void ctulimit::ulimit_bisection(const double& min, const double& max, const double& scale)
+{
+    // copy values to working values
+    double wrk_min = min;
+    double wrk_max = max;
+
+    // Initialise counter
+    int iter=0;
+
+    // Loop until breaking condition is reached
+    while (true) {
+
+        // Log information
+            if (logExplicit()) {
+                log << " Iteration ";
+                log << iter;
+                log << ": Parameter range reduced to [";
+                log << wrk_min * scale;
+                log <<"; ";
+                log << wrk_max * scale;
+                log << "]" << std::endl;
+            }
+
+        // Throw exception if maximum iterations are reached
+        if( iter > m_max_iter) {
+            throw GException::invalid_value(G_UL_BISECTION, "Maximum iterations reached");;
+        }
+
+        // compute center of boundary
+        double mid = (wrk_min + wrk_max) / 2.0;
+
+        // Calculate function value
+        double eval_mid = evaluate(mid);
+
+        // Check for convergence inside tolerance
+        if (std::abs(eval_mid) < m_tol) {
+            break;
+        }
+
+        // change boundaries for further iteration
+        if (eval_mid > 0.0) {
+            wrk_max = mid;
+        }
+        else if (eval_mid < 0.0) {
+            wrk_min = mid;
+        }
+
+        // increment counter
+        iter++;
+
+    } // endwhile
  
     // Return
     return;
+
 }
  
  
@@ -249,12 +434,21 @@ void ctulimit::save(void)
     // Write header
     if (logTerse()) {
         log << std::endl;
-        log.header1("Save TS map");
+        log.header1("Save upper limit");
     }
  
     // Get output filename
     m_outfile = (*this)["outfile"].filename();
  
+    // Create CSV table with 3 columns
+    GCsv table(1, 3);
+
+    table.real(0, 0, m_diff_ulimit);
+    table.real(0, 1, m_flux_ulimit);
+    table.real(0, 2, m_eflux_ulimit);
+
+    // Save CSV table
+    table.save(m_outfile, " ", clobber());
  
     // Return
     return;
@@ -273,13 +467,25 @@ void ctulimit::save(void)
 void ctulimit::init_members(void)
 {
     // Initialise members
-    m_infile.clear();
     m_outfile.clear();
-
+    m_srcname.clear();
+    m_bestloglike = 0.0;
+    m_flux_ulimit = 0.0;
+    m_diff_ulimit = 0.0;
+    m_eflux_ulimit = 0.0;
+    m_dloglike = 0.0;
+    m_tol = 1e-6;
+    m_max_iter = 50;
+    m_eref = 0.0;
+    m_emin = 0.0;
+    m_emax = 0.0;
+    m_sigma_min = 0.0;
+    m_sigma_max  = 0.0;
+    m_skymodel = NULL;
  
     // Initialise protected members
     m_obs.clear();
-
+    m_models.clear();
  
     // Return
     return;
@@ -294,15 +500,27 @@ void ctulimit::init_members(void)
 void ctulimit::copy_members(const ctulimit& app)
 {
     // Copy attributes
-    m_infile   = app.m_infile;
     m_outfile  = app.m_outfile;
-    m_modelfile = app.m_modelfile;
     m_srcname = app.m_srcname;
-    m_parname = app.m_parname;
-    m_algorithm = app.m_algorithm;
+    m_diff_ulimit = app.m_diff_ulimit;
+    m_flux_ulimit = app.m_flux_ulimit;
+    m_eflux_ulimit = app.m_eflux_ulimit;
+    m_bestloglike = app.m_bestloglike;
+    m_dloglike = app.m_dloglike;
+    m_tol = app.m_tol;
+    m_max_iter = app.m_max_iter;
+    m_eref = app.m_eref;
+    m_emin = app.m_emin;
+    m_emax = app.m_emax;
+    m_sigma_min = app.m_sigma_min;
+    m_sigma_max = app.m_sigma_max;
  
     // Copy protected members
     m_obs        = app.m_obs;
+    m_models  = app.m_models;
+
+    // Clone protected members
+    m_skymodel = (app.m_skymodel != NULL) ? app.m_skymodel->clone() : NULL;
  
     // Return
     return;
@@ -334,51 +552,25 @@ void ctulimit::free_members(void)
  ***************************************************************************/
 void ctulimit::get_parameters(void)
 {
-    // If there are no observations in container then add a single CTA
-    // observation using the parameters from the parameter file
+    // If there are no observations in container then load them via user
+    // parameters
     if (m_obs.size() == 0) {
  
-        // Allocate CTA observation
-        GCTAObservation obs;
-
-        // Get event file name
-        std::string filename = (*this)["infile"].filename();
-
-        // Try first to open as FITS file
-        try {
-
-            // Load data
-            obs.load(filename);
-
-            // Set response
-            set_obs_response(&obs);
-
-            // Append observation to container
-            m_obs.append(obs);
-
-        }
-
-        // ... otherwise try to open as XML file
-        catch (GException::fits_open_error &e) {
-
-            // Load observations from XML file
-            m_obs.load(filename);
-
-            // Check if all observations have response information. If
-            // not, get the calibration database parameters and set
-            // the response properly
-            set_response(m_obs);
+        // Throw exception if no input observation file is given
+        require_inobs(G_GET_PARAMETERS);
  
-        } // endcatch: file was an XML file
+        // Build observation container
+        m_obs = get_observations();
  
     } // endif: there was no observation in the container
  
+
     // If there is are no models associated with the observations then
     // load now the model definition
     if (m_obs.models().size() == 0) {
  
         // Get models XML filename
-        std::string filename = (*this)["srcmdl"].filename();
+        std::string filename = (*this)["inmodel"].filename();
  
         // Setup models for optimizing.
         m_obs.models(GModels(filename));
@@ -390,10 +582,34 @@ void ctulimit::get_parameters(void)
     if (!m_obs.models().contains(m_srcname)) {
         std::string msg = "Source \""+m_srcname+"\" not found in model "
                           "container. Please add a source with that name "
-                          "or check for a possible typos.";
-    	throw GException::invalid_value(G_GET_PARAMETERS, msg);
+                          "or check for possible typos.";
+        throw GException::invalid_value(G_GET_PARAMETERS, msg);
+    }
+
+    double CL = (*this)["cl"].real();
+    if (CL != 0.95) {
+        std::string msg = "Confidence level different from 95% requested"
+                          "Currently only 95% CL is possible";
+        throw GException::invalid_value(G_GET_PARAMETERS, msg);
+    }
+    else {
+        // Set Likelihood difference for 95% CL.
+        // See Minuit Handbook
+        m_dloglike = 3.84 / 2.0;
     }
  
+    // Read starting boundaries for bisection
+    m_sigma_min = (*this)["sigma_min"].real();
+    m_sigma_max = (*this)["sigma_max"].real();
+
+    // Read energy values
+    m_eref = (*this)["eref"].real();
+    m_emin = (*this)["emin"].real();
+    m_emax = (*this)["emax"].real();
+
+    // Read precision
+    m_tol = (*this)["tol"].real();
+    m_max_iter = (*this)["max_iter"].integer();
  
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
@@ -401,8 +617,6 @@ void ctulimit::get_parameters(void)
         m_outfile = (*this)["outfile"].filename();
     }
  
-    // TODO read other parameters here
-
     // Return
     return;
 }
@@ -69,7 +69,9 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
-    const double&     ulimit(void) const;
+    const double&     diff_ulimit(void) const;
+    const double&     flux_ulimit(void) const;
+    const double&    eflux_ulimit(void) const;
  
 protected:
     // Protected methods
@@ -77,19 +79,32 @@ protected:
     void copy_members(const ctulimit& app);
     void free_members(void);
     void get_parameters(void);
+    void ulimit_bisection(const double& parmin, const double& parmax, const double& scale);
+    double evaluate(const double& value);
  
     // User parameters
-    std::string              m_infile;     //!< Input file
-    std::string              m_modelfile;//!< Name of model optimised model file
     std::string              m_srcname;    //!< Name of source which is moved around
-    std::string              m_parname; //!< Name of parameter upper limit should be computed
-    std::string              m_outfile;    //!< Output counts map or XML file
-    std::string              m_algorithm; //!< algorithm of upper limit computation
+    std::string              m_outfile;    //!< Output ascii file
+
+    GModelSky*            m_skymodel; //!< Pointer to variable spectral model
+
+    int                m_max_iter; //!< maximum number of iterations
+    double                m_dloglike; //!< Likelihood difference for upper limit computation
+    double                m_bestloglike ; //!< Best fit log likelihood of given model
+    double                m_tol; //!< tolerance for limit determination
+    double                m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values)
+    double                m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values)
+    double                m_eref; //!< reference energy for flux limits
+    double                m_emin; //!< minimum energy for flux limits
+    double                m_emax; //!< maximum energy for flux limits
+
+    double                m_diff_ulimit; //!<Differential upper limit value
+    double                m_flux_ulimit; //!< Flux upper limit value
+    double                m_eflux_ulimit; //!< Energy flux upper limits
  
     // Protected members
     GObservations            m_obs;        //!< Observation container
-
-    double                m_ulimit; //!<Upper limit value
+    GModels m_models;
  
 };
  
@@ -106,15 +121,36 @@ const GObservations&amp; ctulimit::obs(void) const
 }
  
 /***********************************************************************//**
- * @brief Return observation container
+ * @brief Return differential upper limit
  *
- * @return Reference to observation container
+ * @return differential upper limit
  ***************************************************************************/
 inline
-const double& ctulimit::ulimit(void) const
+const double& ctulimit::diff_ulimit(void) const
 {
-    return m_ulimit;
+    return m_diff_ulimit;
 }
  
+/***********************************************************************//**
+ * @brief Return flux upper limit
+ *
+ * @return flux upper limit
+ ***************************************************************************/
+inline
+const double& ctulimit::flux_ulimit(void) const
+{
+    return m_flux_ulimit;
+}
+
+/***********************************************************************//**
+ * @brief return energy flux upper limit
+ *
+ * @return energy flux upper limits
+ ***************************************************************************/
+inline
+const double& ctulimit::eflux_ulimit(void) const
+{
+    return m_eflux_ulimit;
+}
  
 #endif /* CTULIMIT_HPP */
@@ -9,7 +9,7 @@
 #                                                                           #
 #  Author:    Michael Mayer                                                 #
 #             michael.mayer@physik.hu-berlin.de                             #
-#             DESY                                                          #
+#             HU Berlin                                                     #
 #                                                                           #
 #  Purpose:   Parameter file for the ctulimit application                   #
 #                                                                           #
@@ -17,15 +17,26 @@
 #
 # General parameters
 #===================
-infile,   f, a, "events.fits",,,"Input event list, cube or observation definition file"
-srcmdl,   f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
-srcname,  s, a, "Crab",,,"Test source"
-parname,  s, a, "Prefactor",,,"Parameter name for which upper limit should be computed"
-algorithm,  s, a, "Bayesian",,,"Limit calculation algorithm"
+inobs,    f, a, "events.fits",,,"Input event list, counts cube or observation definition file"
+inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+srcname,  s, a, "Crab",,,"Source of interest"
+expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
+psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
-outfile,  f, a, "crab_limit.txt",,, "Output file"
+outfile,  f, a, "upperlimit.txt",,, "Output ascii file"
  
+#
+# Upper limit calculation parameters
+#===================================
+cl,       r, h, 0.95,,,"Confidence Level"
+sigma_min,r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)"
+sigma_max,r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)"
+eref,     r, h, 1.0,,, "Reference energy for differential limits"
+emin,     r, h, 1.0,,, "Minimum energy for flux limits"
+emax,     r, h, 100.0,,, "Maximum energy for flux limits"
+tol,      r, h, 1e-5,,, "Calculation tolerance"
+max_iter, i, h, 50,,, "Maximum number of iterations"
  
 #
 # Standard parameters
@@ -34,4 +45,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
 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, "cttsmap.log",,, "Log filename"
+logfile, f, h, "ctulimit.log",,, "Log filename"
@@ -62,7 +62,7 @@ int main (int argc, char *argv[])
         // Extract error message
         std::string message = e.what();
         std::string signal  = "*** ERROR encounterted in the execution of"
-                              " cttsmap. Run aborted ...";
+                              " ctulimit. Run aborted ...";
  
         // Write error in logger
         application.log << signal  << std::endl;
@@ -46,12 +46,14 @@ ctpsfcube=../src/ctpsfcube/ctpsfcube
 ctbkgcube=../src/ctbkgcube/ctbkgcube
 ctcubemask=../src/ctcubemask/ctcubemask
 ctbutterfly=../src/ctbutterfly/ctbutterfly
+ctulimit=../src/ctulimit/ctulimit
  
  
 #
 # Remove any existing result files
 # ================================
 rm -rf *.fits *.xml
+rm -rf ulimit.dat butterfly.txt
  
  
 #
@@ -597,3 +599,25 @@ else
     exit 1
 fi
 $ECHO " ok"
+
+
+#
+# Test ctulimit
+# ===============
+$ECHO -n "Test ctulimit: "
+$ctulimit inobs="selected_events.fits" \
+    inmodel="results_unbinned.xml" \
+    srcname="Crab" \
+    outfile="ulimit.dat" \
+    caldb="irf" \
+    irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "ulimit.dat" ]
+then
+    $ECHO -n "."
+else
+    $ECHO " ulimit.dat file is not found"
+    exit 1
+fi
+$ECHO " ok"
+
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cttsmap tool.
+#
+# Copyright (C) 2014 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import ctools
+
+
+# ========================== #
+# Test class for cttsmap tool #
+# ========================== #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for cttsmap tool.
+    """
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Set members
+        self.events_name = "data/crab_events.fits"
+        self.model_name  = "result.xml"
+        self.caldb       = "irf"
+        self.irf         = "cta_dummy_irf"
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name("ctulimit")
+
+        # Append tests
+        self.append(self.test_functional, "Test ctulimit functionality")
+
+        # Return
+        return
+
+    # Test cttsmap functionnality
+    def test_functional(self):
+        """
+        Test cttsmap functionnality.
+        """
+        # Set-up cttsmap
+        ulimit = ctools.ctulimit()
+        ulimit["inobs"].filename(self.events_name)
+        ulimit["inmodel"].filename(self.model_name)
+        ulimit["srcname"].string("Crab")
+        ulimit["caldb"].string(self.caldb)
+        ulimit["irf"].string(self.irf)
+        ulimit["outfile"].filename("ulimit.dat")
+        
+        # Run tool
+        self.test_try("Run ctulimit")
+        try:
+            ulimit.run()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in ctulimit.")
+
+        # Save results
+        self.test_try("Save results")
+        try:
+            ulimit.save()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in saving results.")
+        
+        # Return
+        return
@@ -34,6 +34,7 @@ import test_ctpsfcube
 import test_ctbkgcube
 import test_ctcubemask
 import test_ctbutterfly
+import test_ctulimit
 import test_pipelines
  
  
@@ -60,6 +61,7 @@ if __name__ == &#39;__main__&#39;:
     suite_ctbkgcube   = test_ctbkgcube.Test()
     suite_ctcubemask  = test_ctcubemask.Test()
     suite_ctbutterfly = test_ctbutterfly.Test()
+    suite_ctulimit    = test_ctulimit.Test()
     suite_pipelines   = test_pipelines.Test()
  
     # Setup unit tests
@@ -75,6 +77,7 @@ if __name__ == &#39;__main__&#39;:
     suite_ctbkgcube.set()
     suite_ctcubemask.set()
     suite_ctbutterfly.set()
+    suite_ctulimit.set()
     suite_pipelines.set()
  
     # Append tests to container
@@ -90,6 +93,7 @@ if __name__ == &#39;__main__&#39;:
     suites.append(suite_ctbkgcube)
     suites.append(suite_ctcubemask)
     suites.append(suite_ctbutterfly)
+    suites.append(suite_ctulimit)
     suites.append(suite_pipelines)
  
     # Set PFILES environment variable