diff --git a/pyext/ctbkgcube.i b/pyext/ctbkgcube.i
index 57a001a..33e1f11 100644
--- a/pyext/ctbkgcube.i
+++ b/pyext/ctbkgcube.i
@@ -48,6 +48,7 @@ public:
     void                 run(void);
     void                 save(void);
     const GCTAEventCube& cube(void) const;
+    const GModels&       models(void) const;
 };
 
 /***********************************************************************//**
diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp
index aad780b..51232b6 100644
--- a/src/ctbkgcube/ctbkgcube.cpp
+++ b/src/ctbkgcube/ctbkgcube.cpp
@@ -42,6 +42,8 @@
 
 /* __ Coding definitions _________________________________________________ */
 
+/* __ Constants __________________________________________________________ */
+
 
 /*==========================================================================
  =                                                                         =
@@ -233,16 +235,27 @@ void ctbkgcube::run(void)
         log << m_obs << std::endl;
     }
 
+    // Log input model
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Input model");
+        log << m_obs.models() << std::endl;
+    }
+
     // Write header
     if (logTerse()) {
         log << std::endl;
         log.header1("Prepare model");
     }
 
-    // Copy models from observation container
+    // Copy models from observation container and reset output model
+    // container
     m_bkgmdl = m_obs.models();
+    m_outmdl.clear();
+    m_cube_model = -1;
 
-    // Remove all model that are not CTA background models from the
+    // Remove all models that are not CTA background models from the
+    // container and put all removed components in the output
     // container
     int num = m_bkgmdl.size();
     for (int i = num-1; i >= 0; --i) {
@@ -281,8 +294,10 @@ void ctbkgcube::run(void)
             log << std::endl;
         }
 
-        // Remove model if requested
+        // If removal is requested, append model to output container and
+        // remove it from the background model container ...
         if (remove) {
+            m_outmdl.append(*(m_bkgmdl[i]));
             m_bkgmdl.remove(i);
         }
 
@@ -345,6 +360,27 @@ void ctbkgcube::run(void)
         log << n_observations << std::endl;
     }
 
+    // Create a background model for output background cube and append
+    // that model to the input model in place of the original
+    // background model
+    GEnergies energies;
+    for (int i = 0; i < m_bkgcube.ebins(); ++i) {
+        energies.append(m_bkgcube.energy(i));
+    }
+    GModelSpatialDiffuseCube spatial(m_bkgcube.map(), energies, 1.0);
+    GModelSpectralPlaw       spectral(1.0, 0.0, GEnergy(1.0, "TeV"));
+    GCTAModelCubeBackground  model(spatial, spectral);
+    model.name("ctbkgcube default background model");
+    m_cube_model = m_outmdl.size(); // Store the slot number for save()
+    m_outmdl.append(model);
+
+    // Log output model
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Output model");
+        log << m_outmdl << std::endl;
+    }
+
     // Return
     return;
 }
@@ -369,8 +405,9 @@ void ctbkgcube::save(void)
         energies.append(m_ebounds.elogmean(i));
     }
 
-    // Get output filename
-    std::string filename = (*this)["outfile"].filename();
+    // Get output filenames
+    std::string outfile  = (*this)["outfile"].filename();
+    std::string outmodel = (*this)["outmodel"].filename();
 
     // Create empty FITS file
     GFits fits;
@@ -382,7 +419,28 @@ void ctbkgcube::save(void)
     energies.write(fits);
     
     // Save FITS file
-    fits.saveto(filename, clobber());
+    fits.saveto(outfile, clobber());
+
+    // Write output models if filename is valid
+    if ((outmodel.length() > 0) && (gammalib::tolower(outmodel) != "none")) {
+
+        // Set filename of map cube.
+        if (m_cube_model >= 0) {
+            GCTAModelCubeBackground* model =
+                dynamic_cast<GCTAModelCubeBackground*>(m_outmdl[m_cube_model]);
+            if (model != NULL) {
+                GModelSpatialDiffuseCube* spatial =
+                    dynamic_cast<GModelSpatialDiffuseCube*>(model->spatial());
+                if (spatial != NULL) {
+                    spatial->filename(outfile);
+                }
+            }
+        }
+
+        // Save output model container
+        m_outmdl.save(outmodel);
+        
+    }
 
     // Return
     return;
@@ -404,8 +462,10 @@ void ctbkgcube::init_members(void)
     m_obs.clear();
     m_bkgcube.clear();
     m_bkgmdl.clear();
+    m_outmdl.clear();
     m_ebounds.clear();
-    m_livetime = 0.0;
+    m_livetime   = 0.0;
+    m_cube_model = -1;
 
     // Return
     return;
@@ -420,11 +480,13 @@ void ctbkgcube::init_members(void)
 void ctbkgcube::copy_members(const ctbkgcube& app)
 {
     // Copy members
-    m_obs      = app.m_obs;
-    m_bkgcube  = app.m_bkgcube;
-    m_bkgmdl   = app.m_bkgmdl;
-    m_ebounds  = app.m_ebounds;
-    m_livetime = app.m_livetime;
+    m_obs        = app.m_obs;
+    m_bkgcube    = app.m_bkgcube;
+    m_bkgmdl     = app.m_bkgmdl;
+    m_outmdl     = app.m_outmdl;
+    m_ebounds    = app.m_ebounds;
+    m_livetime   = app.m_livetime;
+    m_cube_model = app.m_cube_model;
 
     // Return
     return;
@@ -461,15 +523,15 @@ void ctbkgcube::get_parameters(void)
     // If there are no models associated with the observations then load now
     // the model definition from the XML file
     if (m_obs.models().size() == 0) {
-        if ((*this)["bkgmdl"].is_undefined()) {
-            std::string msg = "No background model definition XML file "
-                              "specified. Please set the \"bkgmdl\" parameter "
-                              "to the background model definition XML "
-                              "filename.";
+        if ((*this)["inmodel"].is_undefined()) {
+            std::string msg = "No model definition XML file specified. "
+                              "Please set the \"inmodel\" parameter to the "
+                              "XML file that contains the background model "
+                              "definition.";
             throw GException::invalid_value(G_GET_PARAMETERS, msg);
         }
-        std::string bkgmdl = (*this)["bkgmdl"].filename();
-        GModels     models(bkgmdl);
+        std::string inmodel = (*this)["inmodel"].filename();
+        GModels     models(inmodel);
         m_obs.models(models);
     }
 
@@ -511,9 +573,10 @@ void ctbkgcube::get_parameters(void)
     
     }
 
-    // Read output filename (if needed)
+    // Read output filenames (if needed)
     if (read_ahead()) {
-        std::string filename = (*this)["outfile"].filename();
+        (*this)["outfile"].filename();
+        (*this)["outmodel"].filename();
     }
 
     // Return
diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp
index 04a27f3..1607e8a 100644
--- a/src/ctbkgcube/ctbkgcube.hpp
+++ b/src/ctbkgcube/ctbkgcube.hpp
@@ -60,6 +60,7 @@ public:
     void                 run(void);
     void                 save(void);
     const GCTAEventCube& cube(void) const;
+    const GModels&       models(void) const;
 
 protected:
     // Protected methods
@@ -75,8 +76,10 @@ protected:
     GObservations m_obs;         //!< Observation container
     GCTAEventCube m_bkgcube;     //!< Background cube
     GModels       m_bkgmdl;      //!< CTA background models
+    GModels       m_outmdl;      //!< Output models
     GEbounds      m_ebounds;     //!< Energy boundaries
     double        m_livetime;    //!< Total livetime for normalization
+    int           m_cube_model;  //!< Cube model index
 };
 
 
@@ -92,4 +95,17 @@ const GCTAEventCube& ctbkgcube::cube(void) const
     return (m_bkgcube);
 }
 
+
+/***********************************************************************//**
+ * @brief Return background model container
+ *
+ * @return Background model container.
+ ***************************************************************************/
+inline
+const GModels& ctbkgcube::models(void) const
+{
+    // Return background model container
+    return (m_outmdl);
+}
+
 #endif /* CTBKGCUBE_HPP */
diff --git a/src/ctbkgcube/ctbkgcube.par b/src/ctbkgcube/ctbkgcube.par
index a875c14..918e9ed 100644
--- a/src/ctbkgcube/ctbkgcube.par
+++ b/src/ctbkgcube/ctbkgcube.par
@@ -18,9 +18,10 @@
 # General parameters
 #===================
 infile,   f, a, "NONE",,, "Event list or observation definition file"
-bkgmdl,   f, a, "NONE",,, "Background model XML file"
+inmodel,  f, a, "NONE",,, "Input (background) model XML file"
 cntmap,   f, a, "NONE",,, "Counts cube for background cube definition"
 outfile,  f, a, "bkgcube.fits",,, "Output background cube file"
+outmodel, f, a, "NONE",,,"Output (background) model XML file"
 #
 # Energy binning parameters
 #==========================
diff --git a/test/test_ctbkgcube.py b/test/test_ctbkgcube.py
index 6c9b27d..f17a87d 100755
--- a/test/test_ctbkgcube.py
+++ b/test/test_ctbkgcube.py
@@ -66,9 +66,10 @@ class Test(gammalib.GPythonTestSuite):
         # Set-up ctbkgcube
         bkgcube = ctools.ctbkgcube()
         bkgcube["infile"].filename(self.events_name)
-        bkgcube["bkgmdl"].filename(self.bkg_model)
+        bkgcube["inmodel"].filename(self.bkg_model)
         bkgcube["cntmap"].filename("NONE")
-        bkgcube["outfile"].filename("psfcube.fits")
+        bkgcube["outfile"].filename("bkgcube.fits")
+        bkgcube["outmodel"].filename("bkgcube.xml")
         bkgcube["ebinalg"].string("LOG")
         bkgcube["emin"].real(0.1)
         bkgcube["emax"].real(100.0)
diff --git a/test/test_ctools.sh b/test/test_ctools.sh
index 146a395..e7c4e99 100755
--- a/test/test_ctools.sh
+++ b/test/test_ctools.sh
@@ -85,7 +85,7 @@ if [ -s "events.fits" ]
 then
   $ECHO -n "."
 else
-  echo " events.fits file is not valid"
+  echo " events.fits file is not found"
   exit 1
 fi
 #
@@ -106,7 +106,7 @@ if [ -s "events2.fits" ]
 then
   $ECHO -n "."
 else
-  echo " events2.fits file is not valid"
+  echo " events2.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -132,7 +132,7 @@ if [ -s "skymap.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " skymap.fits file is not valid"
+  $ECHO " skymap.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -162,7 +162,7 @@ if [ -s "cntmap.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " cntmap.fits file is not valid"
+  $ECHO " cntmap.fits file is not found"
   exit 1
 fi
 #
@@ -185,7 +185,7 @@ if [ -s "cntmap2.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " cntmap2.fits file is not valid"
+  $ECHO " cntmap2.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -208,7 +208,7 @@ if [ -s "modmap1.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " modmap1.fits file is not valid"
+  $ECHO " modmap1.fits file is not found"
   exit 1
 fi
 #
@@ -238,7 +238,7 @@ if [ -s "modmap3.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " modmap3.fits file is not valid"
+  $ECHO " modmap3.fits file is not found"
   exit 1
 fi
 #
@@ -254,7 +254,7 @@ if [ -s "modmap2.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " modmap2.fits file is not valid"
+  $ECHO " modmap2.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -278,7 +278,7 @@ if [ -s "selected_events.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " selected_events.fits file is not valid"
+  $ECHO " selected_events.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -300,7 +300,7 @@ if [ -s "results_binned.xml" ]
 then
   $ECHO -n "."
 else
-  $ECHO " results_binned.xml file is not valid"
+  $ECHO " results_binned.xml file is not found"
   exit 1
 fi
 #
@@ -315,7 +315,7 @@ if [ -s "results_unbinned.xml" ]
 then
   $ECHO -n "."
 else
-  $ECHO " results_unbinned.xml file is not valid"
+  $ECHO " results_unbinned.xml file is not found"
   exit 1
 fi
 #
@@ -330,7 +330,7 @@ if [ -s "results_unbinned_background_cube.xml" ]
 then
   $ECHO -n "."
 else
-  $ECHO " results_unbinned_background_cube.xml file is not valid"
+  $ECHO " results_unbinned_background_cube.xml file is not found"
   exit 1
 fi
 #
@@ -345,7 +345,7 @@ if [ -s "results_binned_background_cube.xml" ]
 then
   $ECHO -n "."
 else
-  $ECHO " results_binned_background_cube.xml file is not valid"
+  $ECHO " results_binned_background_cube.xml file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -373,7 +373,7 @@ if [ -s "tsmap.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " tsmap.fits file is not valid"
+  $ECHO " tsmap.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -405,7 +405,7 @@ if [ -s "expcube1.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " expcube1.fits file is not valid"
+  $ECHO " expcube1.fits file is not found"
   exit 1
 fi
 #
@@ -420,7 +420,7 @@ if [ -s "expcube2.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " expcube2.fits file is not valid"
+  $ECHO " expcube2.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -454,7 +454,7 @@ if [ -s "psfcube1.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " psfcube1.fits file is not valid"
+  $ECHO " psfcube1.fits file is not found"
   exit 1
 fi
 #
@@ -471,7 +471,7 @@ if [ -s "psfcube2.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " psfcube2.fits file is not valid"
+  $ECHO " psfcube2.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -484,9 +484,10 @@ $ECHO -n "Test ctbkgcube: "
 #
 # Run 1
 $ctbkgcube infile="data/crab_events.fits" \
-           bkgmdl="data/crab.xml" \
+           inmodel="data/crab.xml" \
            cntmap="NONE" \
            outfile="bkgcube1.fits" \
+           outmodel="bkgcube1.xml" \
            emin=0.1 \
            emax=100.0 \
            enumbins=20 \
@@ -502,21 +503,36 @@ if [ -s "bkgcube1.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " bkgcube1.fits file is not valid"
+  $ECHO " bkgcube1.fits file is not found"
+  exit 1
+fi
+if [ -s "bkgcube1.xml" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " bkgcube1.xml file is not found"
   exit 1
 fi
 #
 # Run 2
 $ctbkgcube infile="data/crab_events.fits" \
-           bkgmdl="data/crab.xml" \
+           inmodel="data/crab.xml" \
            cntmap="data/crab_cntmap.fits" \
-           outfile="bkgcube2.fits"
+           outfile="bkgcube2.fits" \
+           outmodel="bkgcube2.xml"
 $ECHO -n "."
 if [ -s "bkgcube2.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " bkgcube2.fits file is not valid"
+  $ECHO " bkgcube2.fits file is not found"
+  exit 1
+fi
+if [ -s "bkgcube2.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " bkgcube2.xml file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -539,7 +555,7 @@ if [ -s "filtered_cube.fits" ]
 then
   $ECHO -n "."
 else
-  $ECHO " filtered_cube.fits file is not valid"
+  $ECHO " filtered_cube.fits file is not found"
   exit 1
 fi
 $ECHO " ok"
@@ -562,7 +578,7 @@ if [ -s "butterfly.txt" ]
 then
     $ECHO -n "."
 else
-    $ECHO " butterfly.txt file is not valid"
+    $ECHO " butterfly.txt file is not found"
     exit 1
 fi
 $ECHO " ok"