diff --git a/pyext/cterror.i b/pyext/cterror.i
index 07b9e0d..5ac43b2 100644
--- a/pyext/cterror.i
+++ b/pyext/cterror.i
@@ -1,5 +1,5 @@
 /***************************************************************************
- *                   cterror - Parameter error calculation tool            *
+ *                 cterror - Parameter error calculation tool              *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2015 by Florent Forest                                   *
  * ----------------------------------------------------------------------- *
@@ -38,7 +38,7 @@
  ***************************************************************************/
 class cterror : public ctool {
 
-    public:
+public:
     // Constructors and destructors
     cterror(void);
     explicit cterror(const GObservations& obs);
@@ -52,9 +52,9 @@ class cterror : public ctool {
     void                 save(void);
     const GObservations& obs(void) const;
     const GOptimizer*    opt(void) const;
-
 };
 
+
 /***********************************************************************//**
  * @brief Parameter error calculation tool Python extensions
  ***************************************************************************/
diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py
index a9ab5f8..054a68f 100755
--- a/pyext/ctools/obsutils.py
+++ b/pyext/ctools/obsutils.py
@@ -140,9 +140,10 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
     Parameters:
      obs   - Observation container
     Keywords:
-     log   - Create log file(s)
-     debug - Create screen dump
-     edisp - Apply energy dispersion?
+     log     - Create log file(s)
+     debug   - Create screen dump
+     chatter - Chatter level
+     edisp   - Apply energy dispersion?
     """
     # Allocate ctlike application
     like = ctools.ctlike(obs)
@@ -167,19 +168,21 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
     # Return observations
     return like
 
-# ================ #
-# Run cterror      #
-# ================ #
-def cterror(obs, log=False, debug=False):
+
+# ============================================================== #
+# Fit observations and determine errors using likelihood profile #
+# ============================================================== #
+def cterror(obs, log=False, debug=False, chatter=2):
     """
     Perform maximum likelihood fitting of observations in the container.
     
     Parameters:
      obs   - Observation container
     Keywords:
-     log   - Create log file(s)
-     debug - Create screen dump
-     edisp - Apply energy dispersion?
+     log     - Create log file(s)
+     debug   - Create screen dump
+     chatter - Chatter level
+     edisp   - Apply energy dispersion?
     """
     # Allocate cterror application
     error = ctools.cterror(obs)
@@ -192,6 +195,9 @@ def cterror(obs, log=False, debug=False):
     if debug:
         error["debug"].boolean(True)
 
+    # Set chatter level
+    error["chatter"].chatter = chatter
+
     # Run cterror application.
     error.run()
     
diff --git a/scripts/cspull.py b/scripts/cspull.py
index 24c553b..34c29b2 100755
--- a/scripts/cspull.py
+++ b/scripts/cspull.py
@@ -47,6 +47,8 @@ class cspull(ctools.cscript):
         self.obs       = None
         self.model     = None
         self.m_inmodel = None
+        self.m_edisp   = False
+        self.m_profile = False
         
         # Make sure that parfile exists
         file = self.parfile()
@@ -103,6 +105,7 @@ class cspull(ctools.cscript):
             pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function"))
             pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor"))
             pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
+            pars.append(gammalib.GApplicationPar("profile","b","h","no","","","Use likelihood profile method for errors?"))
             pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
             pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
             pars.append(gammalib.GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
@@ -160,6 +163,7 @@ class cspull(ctools.cscript):
         self.m_ntrials = self["ntrials"].integer()   
         self.m_edisp   = self["edisp"].boolean()
         self.m_offset  = self["offset"].real()   
+        self.m_profile = self["profile"].boolean()
 
         # Set some fixed parameters
         self.m_log   = False                    # Logging in client tools
@@ -297,9 +301,12 @@ class cspull(ctools.cscript):
             self.log("\n")
 
         # Fit model
-        #like = obsutils.fit(obs, edisp=self.m_edisp, \
-        #                    log=self.m_log, debug=self.m_debug)
-        like  = obsutils.cterror(obs, log=self.m_log, debug=self.m_debug)
+        if self.m_profile:
+            like  = obsutils.cterror(obs, \
+                                     log=self.m_log, debug=self.m_debug)
+        else:
+            like = obsutils.fit(obs, edisp=self.m_edisp, \
+                                log=self.m_log, debug=self.m_debug)
 
         # Store results
         logL   = like.opt().value()
diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp
index 95aaf02..8c3fe82 100644
--- a/src/cterror/cterror.cpp
+++ b/src/cterror/cterror.cpp
@@ -1,5 +1,5 @@
 /***************************************************************************
- *                   cterror - Parameter error calculation tool            *
+ *                 cterror - Parameter error calculation tool              *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2015 by Florent Forest                                   *
  * ----------------------------------------------------------------------- *
@@ -72,7 +72,7 @@ cterror::cterror(void) : ctool(CTERROR_NAME, CTERROR_VERSION)
  * observations container.
  ***************************************************************************/
 cterror::cterror(const GObservations& obs) :
-          ctool(CTERROR_NAME, CTERROR_VERSION)
+         ctool(CTERROR_NAME, CTERROR_VERSION)
 {
     // Initialise members
     init_members();
@@ -93,7 +93,7 @@ cterror::cterror(const GObservations& obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 cterror::cterror(int argc, char *argv[]) :
-          ctool(CTERROR_NAME, CTERROR_VERSION, argc, argv)
+         ctool(CTERROR_NAME, CTERROR_VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -198,7 +198,6 @@ void cterror::clear(void)
 
 /***********************************************************************//**
  * @brief Computes 
- *
  ***************************************************************************/
 void cterror::run(void)
 {
@@ -228,11 +227,9 @@ void cterror::run(void)
         log << m_obs << std::endl;
     }
 
-    // Save original models
-    GModels models_orig = m_obs.models();
-
     // Allocate optimizer
     m_opt = new GOptimizerLM();
+
     // Optimize and save best log-likelihood
     m_obs.optimize(*m_opt);
     m_obs.errors(*m_opt);
@@ -249,35 +246,44 @@ void cterror::run(void)
         log << m_obs.models() << std::endl;
     }
 
-
     // Continue only if source model exists
     if (m_obs.models().contains(m_srcname)) {
 
-        GModels& models = const_cast<GModels&>(m_obs.models());
-        m_skymodel      = dynamic_cast<GModelSky*>(models[m_srcname]);
-        if (m_skymodel == NULL) {
+        // Save best fitting models
+        GModels models_best = m_obs.models();
+
+        // Get pointer on sky model
+        GModelSky* model  = dynamic_cast<GModelSky*>(models_best[m_srcname]);
+        if (model == NULL) {
             std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
-                    "Please specify the name of a sky model for "
-                    "parameter error computation.";
+                              "Please specify the name of a sky model for "
+                              "parameter error computation.";
             throw GException::invalid_value(G_RUN, msg);
         }
 
-        int npars = m_skymodel->spectral()->size();
+        // Get number of spectral parameters
+        int npars = model->spectral()->size();
 
+        // Loop over spectral parameters of sky model
+        for (int i = 0; i < npars; ++i) {
 
-        for (int i = 0; i < npars; ++i)
-        {
+            // Skip parameter if it is fixed
+            if (model->spectral()->at(i).is_fixed()) {
+                continue;
+            }
 
-            // Re-optimize
-            m_obs.optimize(*m_opt);
+            // Initialise with best fitting models
+            m_obs.models(models_best);
 
-            if(m_skymodel->spectral()->at(i).is_fixed())
-                continue;
+            // Get pointer on sky model parameter
+            GModels&   current_models = const_cast<GModels&>(m_obs.models());
+            GModelSky* current_model  = dynamic_cast<GModelSky*>(current_models[m_srcname]);
+            m_model_par               = &(current_model->spectral()->at(i));
 
-            m_model_par = &(m_skymodel->spectral()->at(i));
+            // Extract current value
+            m_value = m_model_par->factor_value();
 
             // Compute parameter bracketing
-            m_value  = m_model_par->factor_value();
             double parmin = std::max(m_model_par->factor_min(),
                                      m_value - 10*m_model_par->factor_error());
             double parmax = std::min(m_model_par->factor_max(),
@@ -288,7 +294,7 @@ void cterror::run(void)
                 log << std::endl;
                 log.header1("Compute parameter error");
                 log << gammalib::parformat("Model name");
-                log << m_skymodel->name() << std::endl;
+                log << model->name() << std::endl;
                 log << gammalib::parformat("Parameter name");
                 log << m_model_par->name() << std::endl;
                 log << gammalib::parformat("Confidence level");
@@ -304,43 +310,38 @@ void cterror::run(void)
             }
 
             // Compute lower boundary
-            error_bisection(parmin, m_value);
-            double value_lo = m_model_par->factor_value();
+            double value_lo = error_bisection(parmin, m_value);
 
+            // Dump lower parameter value
             if (logTerse()) {
                 log << gammalib::parformat("Final lower parameter value");
                 log << value_lo << std::endl;
             }
 
-            // Reset parameter to original value
-            m_model_par->factor_value(m_value);
-
             // Compute upper boundary
-            error_bisection(m_value, parmax);
-            double value_hi = m_model_par->factor_value();
+            double value_hi = error_bisection(m_value, parmax);
 
+            // Dump upper parameter value
             if (logTerse()) {
                 log << gammalib::parformat("Final upper parameter");
                 log << value_hi << std::endl;
             }
 
-            // Reset parameter to original value
-            m_model_par->factor_value(m_value);
-
             // Compute errors
-            m_error = 0.5*(value_hi - value_lo);
-            double error_max = std::max(value_hi-m_value, m_value-value_lo);
-            double error_min = std::min(value_hi-m_value, m_value-value_lo);
-
+            m_error          = 0.5 * (value_hi - value_lo);
+            double error_neg = value_lo - m_value;
+            double error_pos = value_hi - m_value;
+            //double error_max = std::max(value_hi-m_value, m_value-value_lo);
+            //double error_min = std::min(value_hi-m_value, m_value-value_lo);
 
             // Write results to logfile
             if (logTerse()) {
                 log << std::endl;
                 log.header1("Error results");
-                log << gammalib::parformat("Max error");
-                log << error_max;
-                log << gammalib::parformat("Min error");
-                log << error_min;
+                log << gammalib::parformat("Negative error");
+                log << error_neg;
+                log << gammalib::parformat("Positive error");
+                log << error_pos;
                 log << gammalib::parformat("Mean parameter error");
                 log << m_error;
                 log << " " << m_model_par->unit();
@@ -351,24 +352,24 @@ void cterror::run(void)
             }
 
             // Save error result
-            m_model_par->factor_error(m_error);
+            model->spectral()->at(i).factor_error(m_error);
 
-        }
+        } // endfor: looped over spectral parameters
 
-    }
+        // Restore best fitting models (now with new errors computed)
+        m_obs.models(models_best);
 
+    } // endif: source model exists
 
     // Return
     return;
 }
 
+
 /***********************************************************************//**
- * @brief Save maps
+ * @brief Save source model
  *
- * This method saves the error to an ascii file.
- *
- * @todo No yet implemented as we have no clear use case yet for saving
- * the result.
+ * This method saves the soutce model to an XML file.
  ***************************************************************************/
 void cterror::save(void)
 {
@@ -378,9 +379,8 @@ void cterror::save(void)
         log.header1("Save results");
     }
 
-
     // Get output filename
-    std::string outmodel = (*this)["outfile"].filename();
+    std::string outmodel = (*this)["outmodel"].filename();
 
      // Write results out as XML model
     if (gammalib::toupper(outmodel) != "NONE") {
@@ -391,6 +391,7 @@ void cterror::save(void)
     return;
 }
 
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -405,8 +406,6 @@ void cterror::init_members(void)
     // Initialise user parameters
     m_srcname.clear();
     m_confidence = 0.68;
-    m_sigma_min  = 0.0;
-    m_sigma_max  = 0.0;
     m_tol        = 1.0e-3;
     m_max_iter   = 50;
     m_value      = 0.0;
@@ -415,9 +414,8 @@ void cterror::init_members(void)
     // Initialise protected members
     m_obs.clear();
     m_dlogL        = 0.0;
-    m_skymodel     = NULL;
-    m_model_par    = NULL;
     m_best_logL    = 0.0;
+    m_model_par    = NULL;
 
     // Return
     return;
@@ -434,15 +432,14 @@ void cterror::copy_members(const cterror& app)
     // Copy user parameters
     m_srcname    = app.m_srcname;
     m_confidence = app.m_confidence;
-    m_sigma_min  = app.m_sigma_min;
-    m_sigma_max  = app.m_sigma_max;
     m_tol        = app.m_tol;
     m_max_iter   = app.m_max_iter;
 
     // Copy protected members
-    m_obs          = app.m_obs;
-    m_dlogL        = app.m_dlogL;
-    m_best_logL    = app.m_best_logL;
+    m_obs       = app.m_obs;
+    m_dlogL     = app.m_dlogL;
+    m_best_logL = app.m_best_logL;
+    m_model_par = NULL;
 
     // Return
     return;
@@ -506,19 +503,13 @@ void cterror::get_parameters(void)
 
     // Get confidence level and transform into log-likelihood difference
     m_confidence = (*this)["confidence"].real();
-
     double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
     m_dlogL      = (sigma*sigma) / 2.0;
 
-    // Read starting boundaries for bisection
-    m_sigma_min = (*this)["sigma_min"].real();
-    m_sigma_max = (*this)["sigma_max"].real();
-
-    // Read precision
+    // Read other parameters
     m_tol      = (*this)["tol"].real();
     m_max_iter = (*this)["max_iter"].integer();
 
-
     // Return
     return;
 }
@@ -532,7 +523,7 @@ void cterror::get_parameters(void)
  *
  * This method calculates the error using a bisection method.
  ***************************************************************************/
-void cterror::error_bisection(const double& min, const double& max)
+double cterror::error_bisection(const double& min, const double& max)
 {
     // Copy values to working values
     double wrk_min = min;
@@ -541,6 +532,9 @@ void cterror::error_bisection(const double& min, const double& max)
     // Initialise iteration counter
     int iter = 0;
 
+    // Initialize mid value
+    double mid = (wrk_min + wrk_max) / 2.0;
+
     // Loop until breaking condition is reached
     while (true) {
 
@@ -556,8 +550,7 @@ void cterror::error_bisection(const double& min, const double& max)
 
         // Throw exception if maximum iterations are reached
         if (iter > m_max_iter) {
-            if(wrk_min - m_model_par->factor_min() < m_tol)
-            {
+            if(wrk_min - m_model_par->factor_min() < m_tol) {
                 std::string msg = "The "+m_model_par->name()+" parameter minimum has been"
                                   " reached during error calculation. To obtain accurate "
                                   " errors, consider to set the minimum to a lower value,"
@@ -567,8 +560,7 @@ void cterror::error_bisection(const double& min, const double& max)
                 }
                 break;
             }
-            else if(m_model_par->factor_max() - wrk_max < m_tol)
-            {
+            else if(m_model_par->factor_max() - wrk_max < m_tol) {
                 std::string msg = "The "+m_model_par->name()+" parameter maximum has been"
                                   " reached during error calculation. To obtain accurate "
                                   " errors, consider to set the maxmimum to a higher value"
@@ -578,8 +570,7 @@ void cterror::error_bisection(const double& min, const double& max)
                 }
                 break;
             }
-            else
-            {
+            else {
                 std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+
                                   " has been reached. You may consider to increase"
                                   " the \"max_iter\" parameter, and re-run cterror.";
@@ -588,7 +579,7 @@ void cterror::error_bisection(const double& min, const double& max)
         }
 
         // Compute center of boundary
-        double mid = (wrk_min + wrk_max) / 2.0;
+        mid = (wrk_min + wrk_max) / 2.0;
 
         // Calculate function value
         double eval_mid = evaluate(mid);
@@ -598,9 +589,9 @@ void cterror::error_bisection(const double& min, const double& max)
             break;
         }
 
-        // If we are on the crescent side of the parabola
-        if (mid > m_value)
-        {
+        // If we are on the crescent side of the parabola ...
+        if (mid > m_value) {
+
             // Change boundaries for further iteration
             if (eval_mid > 0.0) {
                 wrk_max = mid;
@@ -609,9 +600,10 @@ void cterror::error_bisection(const double& min, const double& max)
                 wrk_min = mid;
             }
         }
-        // If we are on the decrescent side of the parabola
-        else
-        {
+
+        // ... otherwise we are on the decrescent side of the parabola
+        else {
+
             // Change boundaries for further iteration
             if (eval_mid > 0.0) {
                 wrk_min = mid;
@@ -626,8 +618,8 @@ void cterror::error_bisection(const double& min, const double& max)
 
     } // endwhile
 
-    // Return
-    return;
+    // Return mid value
+    return mid;
 
 }
 
@@ -652,15 +644,12 @@ double cterror::evaluate(const double& value)
         // Change parameter factor
         m_model_par->factor_value(value);
 
-        //Fix parameter
+        // Fix parameter
         m_model_par->fix();
 
         // Re-optimize
         m_obs.optimize(*m_opt);
 
-        // Evaluate likelihood for new model container
-        m_obs.eval();
-
         // Retrieve likelihood
         logL = m_obs.logL();
 
@@ -673,7 +662,7 @@ double cterror::evaluate(const double& 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);
     }
@@ -681,7 +670,6 @@ double cterror::evaluate(const double& value)
     // Compute function value
     double logL_difference = logL - m_best_logL - m_dlogL;
 
-
     // Return
     return logL_difference;
 }
diff --git a/src/cterror/cterror.hpp b/src/cterror/cterror.hpp
index 96b3e00..d4eecf1 100644
--- a/src/cterror/cterror.hpp
+++ b/src/cterror/cterror.hpp
@@ -1,5 +1,5 @@
 /***************************************************************************
- *                   cterror - Parameter error calculation tool            *
+ *                 cterror - Parameter error calculation tool              *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2015 by Florent Forest                                   *
  * ----------------------------------------------------------------------- *
@@ -70,13 +70,11 @@ protected:
     void   free_members(void);
     void   get_parameters(void);
     double evaluate(const double& value);
-    void   error_bisection(const double& min, const double& max);
+    double error_bisection(const double& min, const double& max);
 
     // User parameters
     std::string   m_srcname;      //!< Name of source which is moved around
     double        m_confidence;   //!< Confidence level
-    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_tol;          //!< Tolerance for limit determination
     int           m_max_iter;     //!< Maximum number of iterations
     double        m_value;        //!< Parameter value 
@@ -85,7 +83,6 @@ protected:
     // Protected members
     GObservations m_obs;          //!< Observation container
     double        m_dlogL;        //!< Likelihood difference for upper limit computation
-    GModelSky*    m_skymodel;     //!< Pointer to sky model
     GModelPar*    m_model_par;    //!< Pointer to model parameter
     double        m_best_logL;    //!< Best fit log likelihood of given model
     GOptimizerLM* m_opt;
diff --git a/src/cterror/cterror.par b/src/cterror/cterror.par
index f719705..8c1e776 100644
--- a/src/cterror/cterror.par
+++ b/src/cterror/cterror.par
@@ -19,20 +19,18 @@
 #===================
 inobs,    f, a, "events.fits",,, "Input event list, counts cube or observation definition file"
 inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
-outfile,  f, a, "results.xml",,, "Save results"
+outmodel, f, a, "results.xml",,, "Source model output file"
 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)"
 bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
-caldb,    s, a, "dummy",,, "Calibration database"
-irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
+caldb,    s, a, "prod2",,, "Calibration database"
+irf,      s, a, "South_50h",,, "Instrument response function"
 
 #
 # Upper limit calculation parameters
 #===================================
 confidence, r, h, 0.68,0.0,1.0, "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)"
 tol,        r, h, 1e-3,,, "Computation tolerance"
 max_iter,   i, h, 50,,, "Maximum number of iterations"
 
diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in
index 5194b18..17b0cd1 100755
--- a/test/test_cscripts.sh.in
+++ b/test/test_cscripts.sh.in
@@ -42,6 +42,7 @@ rm -rf *.fits *.xml *.dat
 #
 # Creates pfiles directory
 # ========================
+rm -rf pfiles
 mkdir -p pfiles
 export PFILES=pfiles
 
@@ -316,3 +317,5 @@ else
   exit 1
 fi
 $ECHO " ok"
+
+ls -la pfiles
diff --git a/test/test_ctools.sh b/test/test_ctools.sh
index d3ef41c..dff866d 100755
--- a/test/test_ctools.sh
+++ b/test/test_ctools.sh
@@ -47,6 +47,7 @@ ctbkgcube=../src/ctbkgcube/ctbkgcube
 ctcubemask=../src/ctcubemask/ctcubemask
 ctbutterfly=../src/ctbutterfly/ctbutterfly
 ctulimit=../src/ctulimit/ctulimit
+cterror=../src/cterror/cterror
 
 
 #
@@ -61,7 +62,7 @@ rm -rf ulimit.dat butterfly.txt
 # ========================
 rm -rf pfiles
 mkdir -p pfiles
-cp -r ../src/*/*.par pfiles/
+#cp -r ../src/*/*.par pfiles/
 export PFILES=pfiles
 
 
@@ -595,3 +596,24 @@ else
 fi
 $ECHO " ok"
 
+
+#
+# Test cterror
+# ============
+$ECHO -n "Test cterror: "
+$cterror inobs="data/crab_events.fits.gz" \
+         inmodel="data/crab.xml" \
+         outmodel="results_error.xml" \
+         srcname="Crab" \
+         caldb="irf" \
+         irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "cterror.log" ]
+then
+    $ECHO -n "."
+else
+    $ECHO " cterror.log file is not found"
+    exit 1
+fi
+$ECHO " ok"
+
diff --git a/test/test_python.py b/test/test_python.py
index 7da0400..911178b 100755
--- a/test/test_python.py
+++ b/test/test_python.py
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the ctools package.
 #
-# Copyright (C) 2012-2014 Juergen Knoedlseder
+# Copyright (C) 2012-2015 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
@@ -35,6 +35,7 @@ import test_ctbkgcube
 import test_ctcubemask
 import test_ctbutterfly
 import test_ctulimit
+import test_cterror
 import test_pipelines
 
 
@@ -62,6 +63,7 @@ if __name__ == '__main__':
     suite_ctcubemask  = test_ctcubemask.Test()
     suite_ctbutterfly = test_ctbutterfly.Test()
     suite_ctulimit    = test_ctulimit.Test()
+    suite_cterror     = test_cterror.Test()
     suite_pipelines   = test_pipelines.Test()
 
     # Setup unit tests
@@ -78,6 +80,7 @@ if __name__ == '__main__':
     suite_ctcubemask.set()
     suite_ctbutterfly.set()
     suite_ctulimit.set()
+    suite_cterror.set()
     suite_pipelines.set()
 
     # Append tests to container
@@ -94,6 +97,7 @@ if __name__ == '__main__':
     suites.append(suite_ctcubemask)
     suites.append(suite_ctbutterfly)
     suites.append(suite_ctulimit)
+    suites.append(suite_cterror)
     suites.append(suite_pipelines)
 
     # Set PFILES environment variable