diff --git a/scripts/Makefile.am b/scripts/Makefile.am
index 7360197..810837b 100644
--- a/scripts/Makefile.am
+++ b/scripts/Makefile.am
@@ -4,10 +4,13 @@
 binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py
 
 # Scripts to be distributed in pythondir
-pythonscripts = cspull.py cssens.py cstsdist.py obsutils.py
+pythonscripts = cspull.py cssens.py cstsdist.py \
+                obsutils.py \
+                lsmodel.py
 
 # Scripts for which symbolic links will be created in bindir
-linkscripts = cspull cssens cstsdist
+linkscripts = cspull cssens cstsdist \
+              lsmodel
 
 # Distribute scripts
 EXTRA_DIST         = $(binscripts) $(pythonscripts)
diff --git a/scripts/cta_root2caldb.py b/scripts/cta_root2caldb.py
index 14fddaa..711defa 100755
--- a/scripts/cta_root2caldb.py
+++ b/scripts/cta_root2caldb.py
@@ -650,7 +650,6 @@ class caldb():
             for ioff in range(noffset):
                 o_lo = offsets.GetBinLowEdge(ioff+1)
                 hdu["THETA_LO"][0,ioff] = o_lo
-
         #
         # THETA_LO
         if not hdu.contains("THETA_HI"):
@@ -683,8 +682,144 @@ class caldb():
 
         # Return boundary information
         return bounds
+
+    def make_3D(self, array, hdu, name, unit, scale=1.0):
+        """
+        Make 3D cube as function of DETX, DETY and energy from a
+        ROOT 2D histogram. If the HDU has already the energy and
+        offset angle columns, this method will simply add another data
+        column.
+        If name==None, the method will not append any data column.
+        
+        Parameters:
+         array - ROOT 2D histogram.
+         hdu   - FITS HDU.
+         name  - Data column name.
+         unit  - Data unit.
+        Keywords:
+         scale - Scaling factor for histogram values.
+        """
+        # Set constants
+        deg2sr = 0.01745329*0.01745329
+        
+        # Extract energy and offset angle vectors
+        energies = array.GetXaxis()
+        neng     = energies.GetNbins()
+        offsets  = array.GetYaxis()
+        noffset  = offsets.GetNbins()
+        ewidth   = [] # in MeV
+        for ieng in range(neng):
+            ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \
+                          pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0))
+
+        # Build DETX and DETY axes
+        ndet     = array.GetYaxis().GetNbins()
+        ndets    = 2*ndet
+        det_max  = array.GetYaxis().GetBinUpEdge(ndet)
+        det_bin  = det_max/float(ndet)
+        dets_lo  = []
+        dets_hi  = []
+        dets2    = []
+        for i in range(ndets):
+            det_lo  = -det_max + i*det_bin
+            det_hi  = det_lo + det_bin
+            det_val = 0.5*(det_lo + det_hi)
+            dets_lo.append(det_lo)
+            dets_hi.append(det_hi)
+            dets2.append(det_val*det_val)
+
+        # Attach columns to HDU. We do this only if they do not yet
+        # exist. This allows adding columns to the same HDU in successive
+        # calls of this method
+        #
+        # DETX_LO
+        if not hdu.contains("DETX_LO"):
+            hdu.append(GFitsTableFloatCol("DETX_LO", 1, ndets))
+            hdu["DETX_LO"].unit("deg")
+            for i in range(ndets):
+                hdu["DETX_LO"][0,i] = dets_lo[i]
+        #
+        # DETX_HI
+        if not hdu.contains("DETX_HI"):
+            hdu.append(GFitsTableFloatCol("DETX_HI", 1, ndets))
+            hdu["DETX_HI"].unit("deg")
+            for i in range(ndets):
+                hdu["DETX_HI"][0,i] = dets_hi[i]
+        #
+        # DETY_LO
+        if not hdu.contains("DETY_LO"):
+            hdu.append(GFitsTableFloatCol("DETY_LO", 1, ndets))
+            hdu["DETY_LO"].unit("deg")
+            for i in range(ndets):
+                hdu["DETY_LO"][0,i] = dets_lo[i]
+        #
+        # DETY_HI
+        if not hdu.contains("DETY_HI"):
+            hdu.append(GFitsTableFloatCol("DETY_HI", 1, ndets))
+            hdu["DETY_HI"].unit("deg")
+            for i in range(ndets):
+                hdu["DETY_HI"][0,i] = dets_hi[i]
+        #
+        # ENERG_LO
+        if not hdu.contains("ENERG_LO"):
+            hdu.append(GFitsTableFloatCol("ENERG_LO", 1, neng))
+            hdu["ENERG_LO"].unit("TeV")
+            for ieng in range(neng):
+                e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1))
+                hdu["ENERG_LO"][0,ieng] = e_lo
+        #
+        # ENERG_HI
+        if not hdu.contains("ENERG_HI"):
+            hdu.append(GFitsTableFloatCol("ENERG_HI", 1, neng))
+            hdu["ENERG_HI"].unit("TeV")
+            for ieng in range(neng):
+                e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1))
+                hdu["ENERG_HI"][0,ieng] = e_hi
+        #
+        # "NAME"
+        if name != None and not hdu.contains(name):
+            hdu.append(GFitsTableFloatCol(name, 1, ndets*ndets*neng))
+            hdu[name].unit(unit)
+            hdu[name].dim([ndets,ndets,neng])
+            for ix in range(ndets):
+                for iy in range(ndets):
+                    for ieng in range(neng):
+                        index = ix + (iy + ieng * ndets) * ndets
+                        theta = math.sqrt(dets2[ix] + dets2[iy])
+                        ioff  = offsets.FindBin(theta)
+                        if ioff > 0 and ioff <= noffset:
+                            binsq = array.GetBinContent(ieng+1,ioff)
+                            if binsq >= 0:
+                                value = binsq / deg2sr / (ewidth[ieng])
+                                hdu[name][0,index] = value * scale
+        #
+        # Debugging
+        if False:
+            for ieng in range(neng):
+                sum = 0.0
+                for ioff in range(noffset):
+                    index = ieng + ioff * neng
+                    binsq = array.GetBinContent(ieng+1,ioff+1)
+                    value = binsq / deg2sr / ewidth[ieng]
+                    sum  += value
+                    if ioff == 0:
+                        peak = value
+                        peaksq = binsq
+                print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \
+                      pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq)
+
+        # Collect boundary information
+        bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max)
+        bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max)
+        bd_eng  = "ENERG(%.4f-%.2f)TeV" % \
+                  (pow(10.0, energies.GetBinLowEdge(1)), \
+                   pow(10.0, energies.GetBinUpEdge(neng)))
+        bounds  = [bd_detx, bd_dety, bd_eng]
+
+        # Return boundary information
+        return bounds
     
-    def root2caldb(self, filename, rebin=False, psftype="Gauss"):
+    def root2caldb(self, filename, rebin=False, psftype="Gauss", scale=1.0):
         """
         Translate ROOT to CALDB information.
         
@@ -706,7 +841,8 @@ class caldb():
         self.root2edisp(file)
 
         # Create background
-        self.root2bgd(file)
+        #self.root2bgd(file)
+        self.root2bgd3D(file, scale=scale)
         
         # Return
         return
@@ -1057,6 +1193,42 @@ class caldb():
         # Return
         return
 
+    def root2bgd3D(self, file, scale=1.0):
+        """
+        Translate ROOT to CALDB background extension. The following ROOT
+        histograms are used:
+
+        BGRatePerSqDeg_offaxis -> BGD
+        BGRatePerSqDeg_offaxis -> BGD_RECO
+
+        Parameters:
+         file - ROOT file.
+        Keywords:
+         scale - Background rate scaling factor
+        """
+        # Continue only if background HDU has been opened
+        if self.hdu_bgd != None:
+
+            # Allocate ROOT 2D array
+            array = TH2F()
+            file.GetObject("BGRatePerSqDeg_offaxis", array)
+
+            # Set boundaries
+            bounds = self.make_3D(array, self.hdu_bgd, None, "deg")
+            for b in bounds:
+                self.bgd_bounds.append(b)
+            self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \
+                                  self.bgd_bounds, self.bgd_desc)
+
+            # BGD
+            self.make_3D(array, self.hdu_bgd, "BGD", "1/s/MeV/sr", scale=scale)
+
+            # BGD_RECO
+            self.make_3D(array, self.hdu_bgd, "BGD_RECO", "1/s/MeV/sr", scale=scale)
+            
+        # Return
+        return
+
     def set_cif_keywords(self, hdu, name, bounds, desc):
         """
         Set standard CIF keywords for extension.
@@ -1095,13 +1267,14 @@ def set_test():
     # Set database attributes
     path    = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
     rebin   = True
+    scale   = 1.0
 
     # Set database content
     db = [{'inst': "e", 'id':   "IFAE20120510_50h",
-           'path': path, 'rebin': rebin, 'psftype': "Gauss",
+           'path': path, 'rebin': rebin, 'psftype': "Gauss", 'scale': scale,
            'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
           {'inst': "e", 'id':   "IFAE20120510_50h_King",
-           'path': path, 'rebin': rebin, 'psftype': "King",
+           'path': path, 'rebin': rebin, 'psftype': "King", 'scale': scale,
            'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}]
 
     # Return database
@@ -1119,16 +1292,17 @@ def set_prod1_ifae():
     path    = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
     rebin   = True
     psftype = "Gauss"
+    scale   = 1.0
 
     # Set database content
     db = [{'inst': "b", 'id':   "IFAE20120510_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "SubarrayB_IFAE_50hours_20120510_offaxis.root"},
           {'inst': "e", 'id':   "IFAE20120510_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
           {'inst': "i", 'id':   "IFAE20120510_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "SubarrayI_IFAE_50hours_20120510_offaxis.root"}]
 
     # Return database
@@ -1146,76 +1320,77 @@ def set_prod2_desy():
     path    = "/project-data/cta/performance/prod2/Performance_DESY_20140128"
     rebin   = False
     psftype = "Gauss"
+    scale   = 1.0/180000.0
 
     # Set database content
     db = [{'inst': "aar", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"},
           {'inst': "aar", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
           {'inst': "aar", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
 
           {'inst': "aar500", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
           {'inst': "aar500", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
           {'inst': "aar500", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
 
           {'inst': "leoncito", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-LeoncitoPP-NS.S.2a.180000s.root"},
           {'inst': "leoncito", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-LeoncitoPP-NS.S.2a.180000s.root"},
           {'inst': "leoncito", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-LeoncitoPP-NS.S.2a.180000s.root"},
 
           {'inst': "sac", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-SAC100-NS.S.2a.180000s.root"},
           {'inst': "sac", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-SAC100-NS.S.2a.180000s.root"},
           {'inst': "sac", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-SAC100-NS.S.2a.180000s.root"},
 
           {'inst': "spm", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-SPM-NS.N.2NN.180000s.root"},
           {'inst': "spm", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-SPM-NS.N.2NN.180000s.root"},
           {'inst': "spm", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-SPM-NS.N.2NN.180000s.root"},
 
           {'inst': "tenerife", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
           {'inst': "tenerife", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
           {'inst': "tenerife", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
 
           {'inst': "us", 'id':   "DESY20140105_50h",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-US-NS.N.2NN.180000s.root"},
           {'inst': "us", 'id':   "DESY20140105_50h_0deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-US-NS.N.2NN.180000s.root"},
           {'inst': "us", 'id':   "DESY20140105_50h_180deg",
-           'path': path, 'rebin': rebin, 'psftype': psftype,
+           'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
            'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-US-NS.N.2NN.180000s.root"}]
 
     # Return database
@@ -1259,6 +1434,7 @@ if __name__ == '__main__':
         id       = entry['id']
         rebin    = entry['rebin']
         psftype  = entry['psftype']
+        scale    = entry['scale']
         filename = path+"/"+entry['file']
 
         # Allocate caldb entry
@@ -1268,4 +1444,4 @@ if __name__ == '__main__':
         irf.open("file")
     
         # Translate ROOT to CALDB information
-        irf.root2caldb(filename, rebin=rebin, psftype=psftype)
+        irf.root2caldb(filename, rebin=rebin, psftype=psftype, scale=scale)
diff --git a/scripts/lsmodel.py b/scripts/lsmodel.py
new file mode 100755
index 0000000..6cfd031
--- /dev/null
+++ b/scripts/lsmodel.py
@@ -0,0 +1,203 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This script generates a model map from a counts map and a source model.
+#
+# 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 sys
+
+
+# ============= #
+# lsmodel class #
+# ============= #
+class lsmodel(gammalib.GApplication):
+    """
+    This class implements the model map generation script. It derives from
+    the GammaLib::GApplication class which provides support for parameter
+    files, command line arguments, and logging. In that way the Python
+    script behaves just as a regular ctool.
+    """
+    def __init__(self, *argv):
+        """
+        Constructor.
+        """
+        # Set name
+        self.name    = "lsmodel"
+        self.version = "0.1.0"
+        
+        # Make sure that parfile exists
+        file = self.parfile()
+
+        # Initialise application
+        if len(argv) == 0:
+            gammalib.GApplication.__init__(self, self.name, self.version)
+        elif len(argv) ==1:
+            gammalib.GApplication.__init__(self, self.name, self.version, *argv)
+        else:
+            raise TypeError("Invalid number of arguments given.")
+
+        # Set logger properties
+        self.log_header()
+        self.log.date(True)
+
+        # Return
+        return
+    
+    def __del__(self):
+        """
+        Destructor.
+        """
+        #  Write separator into logger
+        if self.logTerse():
+            self.log("\n")
+        
+        # Return
+        return
+
+    def parfile(self):
+        """
+        Check if parfile exists. If parfile does not exist then create a
+        default parfile. This kluge avoids shipping the script with a parfile.
+        """
+        # Set parfile name
+        parfile = self.name+".par"
+        
+        try:
+            pars = gammalib.GApplicationPars(parfile)
+        except:
+            # Signal if parfile was not found
+            sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
+            
+            # Create default parfile
+            pars = gammalib.GApplicationPars()
+            pars.append(gammalib.GApplicationPar("cntmap","f","a","cntmap.fits","","","Counts map"))
+            pars.append(gammalib.GApplicationPar("expmap","f","a","binned_expmap.fits","","","Exposure map"))
+            pars.append(gammalib.GApplicationPar("ltcube","f","a","ltcube.fits","","","Livetime cube"))
+            pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab_model.xml","","","Source model"))
+            pars.append(gammalib.GApplicationPar("caldb","s","h","","","","Calibration database"))
+            pars.append(gammalib.GApplicationPar("irf","s","a","P7REP_SOURCE_V15","","","Instrument response function"))
+            pars.append(gammalib.GApplicationPar("outfile","f","a","model_cube.fits","","","Output file name"))
+            pars.append_standard()
+            pars.save(parfile)
+        
+        # Return
+        return
+        
+    def get_parameters(self):
+        """
+        Get parameters from parfile and setup the observation.
+        """
+        # Get parameters
+        self.m_cntmap  = self["cntmap"].filename()
+        self.m_expmap  = self["expmap"].filename()
+        self.m_ltcube  = self["ltcube"].filename()
+        self.m_srcmdl  = self["srcmdl"].filename()
+        self.m_caldb   = self["caldb"].string()
+        self.m_irf     = self["irf"].string()
+        self.m_outfile = self["outfile"].filename()
+
+        # Set some fixed parameters
+        self.m_log   = False # Logging in client tools
+        self.m_debug = False # Debugging in client tools
+        
+        # Return
+        return
+    
+    def execute(self):
+        """
+        Execute the script.
+        """
+        # Run the script
+        self.run()
+        
+        # Return
+        return
+
+    def run(self):
+        """
+        Run the script.
+        """
+        # Switch screen logging on in debug mode
+        if self.logDebug():
+            self.log.cout(True)
+
+        # Get parameters
+        self.get_parameters()
+        
+        #  Write input parameters into logger
+        if self.logTerse():
+            self.log_parameters()
+            self.log("\n")
+        
+        # Write header
+        if self.logTerse():
+            self.log("\n")
+            self.log.header1("Generate model cube")
+
+        # Set LAT observation
+        obs = gammalib.GLATObservation()
+        obs.load_binned(self.m_cntmap, self.m_expmap, self.m_ltcube)
+        obs.response(self.m_irf, self.m_caldb)
+
+        # Set model (use first model for the moment)
+        models = gammalib.GModels(self.m_srcmdl)
+        model  = models[0]
+
+        # Get deep copy of GLATEventCube from observation
+        cube = gammalib.GLATEventCube(obs.events())
+
+        # Fill cube with model values
+        for i in range(cube.size()):
+            event = cube[i]
+            value = model.eval(event, obs) * event.size()
+            event.counts(value)
+            cube[i] = event
+
+        # Log cube
+        if self.logTerse():
+            self.log(str(cube))
+            self.log("\n")
+
+        # Write header
+        if self.logTerse():
+            self.log("\n")
+            self.log.header1("Save model cube")
+
+        # Save model cube
+        cube.save(self.m_outfile, self.clobber())
+
+        # Return
+        return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Generates a model map from a counts map and a source model.
+    """
+    # Create instance of application
+    app = lsmodel(sys.argv)
+    
+    # Open logfile
+    app.logFileOpen()
+    
+    # Execute application
+    app.execute()
+    
\ No newline at end of file