diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index bc9ff33..adecb3a 100755
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -43,78 +43,42 @@
 # Scripts to be distributed in $(pyexecdir)/cscripts
 
 pythonscripts = $(top_builddir)/cscripts/__init__.py \
-
                 $(srcdir)/csbkgmodel.py \
-
                 $(srcdir)/cscaldb.py \
-
                 $(srcdir)/csebins.py \
-
                 $(srcdir)/csinfo.py \
-
                 $(srcdir)/cslightcrv.py \
-
                 $(srcdir)/csmodelinfo.py \
-
                 $(srcdir)/csmodelmerge.py \
-
                 $(srcdir)/csmodelselect.py \
-
                 $(srcdir)/csmodelsois.py \
-
                 $(srcdir)/csobs2caldb.py \
-
                 $(srcdir)/csobsdef.py \
-
                 $(srcdir)/csobsinfo.py \
-
                 $(srcdir)/csobsselect.py \
-
                 $(srcdir)/csphagen.py \
-
                 $(srcdir)/csphasecrv.py \
-
                 $(srcdir)/cspull.py \
-
                 $(srcdir)/csresmap.py \
-
                 $(srcdir)/csresspec.py \
-
                 $(srcdir)/csroot2caldb.py  \
-
                 $(srcdir)/cssens.py \
-
                 $(srcdir)/csspec.py \
-
                 $(srcdir)/cssrcdetect.py \
-
                 $(srcdir)/cstsdist.py \
-
                 $(srcdir)/cstsmapmerge.py \
-
                 $(srcdir)/cstsmapsplit.py \
-
                 $(srcdir)/csviscube.py \
-
                 $(srcdir)/csscs.py \
                 $(srcdir)/csworkflow.py \
-
                 $(srcdir)/csfindobs.py \
-
                 $(srcdir)/csiactcopy.py \
-
                 $(srcdir)/csiactdata.py \
-
                 $(srcdir)/csiactobs.py \
-
                 $(srcdir)/obsutils.py \
-
                 $(srcdir)/modutils.py \
-
                 $(srcdir)/calutils.py \
-
                 $(srcdir)/ioutils.py \
-
                 $(srcdir)/mputils.py 	
 
 
@@ -122,64 +86,35 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
 # Parfiles to be distributed in $(prefix)/syspfiles
 
 parfiles = $(srcdir)/csbkgmodel.par \
-
            $(srcdir)/cscaldb.par \
-
            $(srcdir)/csebins.par \
-
            $(srcdir)/cslightcrv.par \
-
            $(srcdir)/csmodelinfo.par \
-
            $(srcdir)/csmodelmerge.par \
-
            $(srcdir)/csmodelselect.par \
-
            $(srcdir)/csmodelsois.par \
-
            $(srcdir)/csobs2caldb.par \
-
            $(srcdir)/csobsdef.par \
-
            $(srcdir)/csobsinfo.par \
-
            $(srcdir)/csobsselect.par \
-
            $(srcdir)/csphagen.par \
-
            $(srcdir)/csphasecrv.par \
-
            $(srcdir)/cspull.par \
-
            $(srcdir)/csresmap.par \
-
            $(srcdir)/csresspec.par \
-
            $(srcdir)/csroot2caldb.par \
-
            $(srcdir)/cssens.par \
-
            $(srcdir)/csspec.par \
-
            $(srcdir)/cssrcdetect.par \
-
            $(srcdir)/cstsdist.par \
-
            $(srcdir)/cstsmapmerge.par \
-
            $(srcdir)/cstsmapsplit.par \
-
            $(srcdir)/csviscube.par \
-
            $(srcdir)/csscs.par \
            $(srcdir)/csworkflow.par \
-
            $(srcdir)/csfindobs.par \
-
            $(srcdir)/csiactcopy.par \
-
            $(srcdir)/csiactdata.par \
-
            $(srcdir)/csiactobs.par
 
 
@@ -187,66 +122,36 @@ parfiles = $(srcdir)/csbkgmodel.par \
 # Scripts for which symbolic links will be created in bindir
 
 linkscripts = csbkgmodel \
-
               cscaldb \
-
               csebins \
-
               csinfo \
-
               cslightcrv \
-
               csmodelinfo \
-
               csmodelmerge \
-
               csmodelselect \
-
               csmodelsois \
-
               csobs2caldb \
-
               csobsdef \
-
               csobsinfo \
-
               csobsselect \
-
               csphagen \
-
               csphasecrv \
-
               cspull \
-
               csresmap \
-
               csresspec \
-
               csroot2caldb \
-
               cssens \
-
               csspec \
-
               cssrcdetect \
-
               cstsdist \
-
               cstsmapmerge \
-
               cstsmapsplit \
-
               csviscube \
-
               csscs \
               csworkflow \
-
               csfindobs \
-
               csiactcopy \
-
               csiactdata \
-
               csiactobs
 
 
@@ -254,78 +159,42 @@ linkscripts = csbkgmodel \
 # Test scripts to be distributed in $(pyexecdir)/cscripts/tests
 
 testscripts = $(srcdir)/tests/__init__.py \
-
               $(top_srcdir)/test/testing.py \
-
               $(top_srcdir)/test/test_python_cscripts.py \
-
               $(top_srcdir)/test/test_cscript.py \
-
               $(top_srcdir)/test/test_cscaldb.py \
-
               $(top_srcdir)/test/test_csbkgmodel.py \
-
               $(top_srcdir)/test/test_csebins.py \
-
               $(top_srcdir)/test/test_csfindobs.py \
-
               $(top_srcdir)/test/test_csinfo.py \
-
               $(top_srcdir)/test/test_cslightcrv.py \
-
               $(top_srcdir)/test/test_csmodelinfo.py \
-
               $(top_srcdir)/test/test_csmodelmerge.py \
-
               $(top_srcdir)/test/test_csmodelselect.py \
-
               $(top_srcdir)/test/test_csmodelsois.py \
-
               $(top_srcdir)/test/test_csobs2caldb.py \
-
               $(top_srcdir)/test/test_csobsdef.py \
-
               $(top_srcdir)/test/test_csobsinfo.py \
-
               $(top_srcdir)/test/test_csobsselect.py \
-
               $(top_srcdir)/test/test_csphagen.py \
-
               $(top_srcdir)/test/test_csphasecrv.py \
-
               $(top_srcdir)/test/test_cspull.py \
-
               $(top_srcdir)/test/test_csresmap.py \
-
               $(top_srcdir)/test/test_csresspec.py \
-
               $(top_srcdir)/test/test_csroot2caldb.py \
-
               $(top_srcdir)/test/test_cssens.py \
-
               $(top_srcdir)/test/test_csspec.py \
-
               $(top_srcdir)/test/test_cssrcdetect.py \
-
               $(top_srcdir)/test/test_cstsdist.py \
-
               $(top_srcdir)/test/test_cstsmapmerge.py \
-
               $(top_srcdir)/test/test_cstsmapsplit.py \
-
               $(top_srcdir)/test/test_csviscube.py \
-
               $(top_srcdir)/test/test_csscs.py \
               $(top_srcdir)/test/test_csworkflow.py \
-
               $(top_srcdir)/test/test_obsutils.py \
-
               $(top_srcdir)/test/test_mputils.py \
-
               $(top_srcdir)/test/test_csiactcopy.py \
-
               $(top_srcdir)/test/test_csiactdata.py \
-
               $(top_srcdir)/test/test_csiactobs.py
 
 
@@ -333,153 +202,79 @@ testscripts = $(srcdir)/tests/__init__.py \
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data
 
 testdata = $(top_srcdir)/test/data/crab.xml \
-
            $(top_srcdir)/test/data/obs_unbinned.xml \
-
            $(top_srcdir)/test/data/obs_unbinned_two.xml \
-
            $(top_srcdir)/test/data/obs_stacked.xml \
-
            $(top_srcdir)/test/data/obs_stacked_two.xml \
-
            $(top_srcdir)/test/data/crab_events.fits \
-
            $(top_srcdir)/test/data/crab_cntmap.fits \
-
            $(top_srcdir)/test/data/crab_expcube.fits.gz \
-
            $(top_srcdir)/test/data/crab_psfcube.fits.gz \
-
            $(top_srcdir)/test/data/crab_edispcube.fits.gz \
-
            $(top_srcdir)/test/data/crab_bkgcube.fits.gz \
-
            $(top_srcdir)/test/data/crab_onoff.xml \
-
            $(top_srcdir)/test/data/crab_offaxis1.fits.gz \
-
            $(top_srcdir)/test/data/crab_offaxis2.fits.gz \
-
            $(top_srcdir)/test/data/crab_exclusion.fits.gz \
-
            $(top_srcdir)/test/data/crab_src_reg.reg \
-
            $(top_srcdir)/test/data/crab_bkg_reg.reg \
-
            $(top_srcdir)/test/data/tsmap_0.fits \
-
            $(top_srcdir)/test/data/tsmap_1.fits \
-
            $(top_srcdir)/test/data/exclusion.reg \
-
            $(top_srcdir)/test/data/pntdef_minimal.dat \
-
            $(top_srcdir)/test/data/pntdef_maximal.dat \
-
            $(top_srcdir)/test/data/models.txt \
-
            $(top_srcdir)/test/data/tsmaps.txt \
-
            $(top_srcdir)/test/data/model_cube_background1.xml \
-
            $(top_srcdir)/test/data/model_cube_background2.xml \
-
            $(top_srcdir)/test/data/crab_bkgcube.xml \
-
            $(top_srcdir)/test/data/irf.root \
-
            $(top_srcdir)/test/data/irf_file.fits.gz \
-
            $(top_srcdir)/test/data/workflow.xml \
-
            $(top_srcdir)/test/data/butterfly.txt \
-
            $(top_srcdir)/test/data/pha.fits \
-
            $(top_srcdir)/test/data/pha_channels.fits \
-
            $(top_srcdir)/test/data/pull.dat \
-
            $(top_srcdir)/test/data/sensitivity.dat \
-
            $(top_srcdir)/test/data/sensitivity_integral.dat \
-
            $(top_srcdir)/test/data/skymap.fits \
-
            $(top_srcdir)/test/data/spectrum.fits \
-
            $(top_srcdir)/test/data/residual.fits \
-
            $(top_srcdir)/test/data/residual_onoff.fits \
-
            $(top_srcdir)/test/data/lightcurve.fits \
-
            $(top_srcdir)/test/data/ts.dat \
-
            $(top_srcdir)/test/data/models_spatial.xml \
-
            $(top_srcdir)/test/data/models_spectral.xml \
-
            $(top_srcdir)/test/data/model_filefunction.txt \
-
            $(top_srcdir)/test/data/model_map.fits.gz \
-
            $(top_srcdir)/test/data/model_mapcube.fits.gz \
-
            $(top_srcdir)/test/data/iact_runlist.dat \
-
            $(top_srcdir)/test/data/iact_runlist_1.dat \
-
            $(top_srcdir)/test/data/iact_runlist_2.dat \
-
            $(top_srcdir)/test/data/model_temporal_phasecurve.xml \
-
            $(top_srcdir)/test/data/model_temporal_phasecurve.fits \
-
            $(top_srcdir)/test/data/phased_events.fits.gz \
-
            $(top_srcdir)/test/data/onoff_arf.fits.gz \
-
            $(top_srcdir)/test/data/onoff_arf1.fits.gz \
-
            $(top_srcdir)/test/data/onoff_arf2.fits.gz \
-
            $(top_srcdir)/test/data/onoff_model.xml \
-
            $(top_srcdir)/test/data/onoff_obs.xml \
-
            $(top_srcdir)/test/data/onoff_obs1.xml \
-
            $(top_srcdir)/test/data/onoff_obs2.xml \
-
            $(top_srcdir)/test/data/onoff_obs_stacked.xml \
-
            $(top_srcdir)/test/data/onoff_pha_off.fits.gz \
-
            $(top_srcdir)/test/data/onoff_pha_off1.fits.gz \
-
            $(top_srcdir)/test/data/onoff_pha_off2.fits.gz \
-
            $(top_srcdir)/test/data/onoff_pha_on.fits.gz \
-
            $(top_srcdir)/test/data/onoff_pha_on1.fits.gz \
-
            $(top_srcdir)/test/data/onoff_pha_on2.fits.gz \
-
            $(top_srcdir)/test/data/onoff_rmf.fits.gz \
-
            $(top_srcdir)/test/data/onoff_rmf1.fits.gz \
-
            $(top_srcdir)/test/data/onoff_rmf2.fits.gz \
-
            $(top_srcdir)/test/data/model_crab_radialacceptance.xml \
-
            $(top_srcdir)/test/data/sigdist_map.fits.gz \
-
            $(top_srcdir)/test/data/variability.fits.gz \
-
            $(top_srcdir)/test/data/viscube.fits.gz \
-
            $(top_srcdir)/test/data/hess_bkg_lookup.fits
 
 
@@ -487,9 +282,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata
 
 testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
-
                $(top_srcdir)/test/data/iactdata/master.json \
-
                $(top_srcdir)/test/data/iactdata/obs-index.fits.gz
 
 
@@ -497,23 +290,14 @@ testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata/events
 
 testiactevents = $(top_srcdir)/test/data/iactdata/events/irf_file.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_0.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_1.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_2.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_3.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_4.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_5.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_6.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_7.fits.gz \
-
                  $(top_srcdir)/test/data/iactdata/events/events_8.fits.gz
 
 
@@ -557,9 +341,7 @@ testiactevents_DATA = $(testiactevents)
 # Files to be distributed in the tar file
 
 EXTRA_DIST = $(binscripts) $(pythonscripts) $(parfiles) \
-
              $(testscripts) $(testdata) $(testiactdata) \
-
              $(testiactevents)
 
 
@@ -569,23 +351,15 @@ EXTRA_DIST = $(binscripts) $(pythonscripts) $(parfiles) \
 install-data-hook:
 
 	@binlist='$(binscripts)';\
-
 	for p in $$binlist; do \
-
 	  chmod +x "$(DESTDIR)$(bindir)/$$p"; \
-
 	done
 
 	@linklist='$(linkscripts)';\
-
 	for p in $$linklist; do \
-
 		rm -rf "$(DESTDIR)$(bindir)/$$p"; \
-
 		cp "$(DESTDIR)$(pyexecdir)/cscripts/$$p.py" "$(DESTDIR)$(bindir)/$$p"; \
-
 		chmod +x "$(DESTDIR)$(bindir)/$$p"; \
-
 	done
 
 
@@ -595,11 +369,8 @@ install-data-hook:
 uninstall-hook:
 
 	@linklist='$(linkscripts)';\
-
 	for p in $$linklist; do \
-
 		rm -rf "$(DESTDIR)$(bindir)/$$p"; \
-
 	done
 
 
@@ -609,21 +380,13 @@ uninstall-hook:
 clean-local:
 
 	@if ! [ '@top_srcdir@' -ef '@top_builddir@' ]; then \
-
 		rm -rf $(top_builddir)/cscripts/cs*.py; \
-
 		rm -rf $(top_builddir)/cscripts/obsutils.py; \
-
 		rm -rf $(top_builddir)/cscripts/modutils.py; \
-
 		rm -rf $(top_builddir)/cscripts/calutils.py; \
-
 		rm -rf $(top_builddir)/cscripts/ioutils.py; \
-
                 rm -rf $(top_builddir)/cscripts/mputils.py; \
-
 		rm -rf $(top_builddir)/cscripts/*.par; \
-
 	fi
 
 	rm -rf $(top_builddir)/cscripts/*.pyc
diff --git a/cscripts/csscs.par b/cscripts/csscs.par
index 93cafc0..df65ef9 100644
--- a/cscripts/csscs.par
+++ b/cscripts/csscs.par
@@ -14,13 +14,56 @@
 #                                                                           #
 #############################################################################
 #
-# General parameters
-#===================
-# Your parameters go here
+# Input and output parameters
+#============================
+inobs,       f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file"
+inmodel,     f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
+srcnames,    s, a, Crab,,, "Semicolon separated list of sources"
+expcube,     f, a, NONE,,, "Input exposure cube file"
+psfcube,     f, a, NONE,,, "Input PSF cube file"
+edispcube,   f, a, NONE,,, "Input energy dispersion cube file"
+bkgcube,     f, a, NONE,,, "Input background cube file"
+caldb,       s, a, prod2,,, "Calibration database"
+irf,         s, a, South_0.5h,,, "Instrument response function"
+inexclusion, f, h, NONE,,, "Exclusion regions"
+edisp,       b, h, no,,, "Apply energy dispersion?"
+outfile,     f, a, map.fits,,, "Output map"
+
+#
+# Script parameters
+#==================
+coordsys,      s, a, CEL,CEL|GAL,, "Coordinate System"
+proj,          s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
+xref,          r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
+yref,          r, a, 22.51,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+nxpix,         i, a, 20,,, "Size of the X axis in pixels"
+nypix,         i, a, 20,,, "Size of the Y axis in pixels"
+binsz,         r, a, 0.1,,, "Pixel size (deg/pixel)"
+roisz,         r, a, 2.,1.42,1000, "Size of the region of interest (numbers of binsz)"
+statistic,     s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
+method,        s, a, 3D, 3D|ONOFF,, "Analysis method"
+ebinalg,       s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins"
+emin,          r, a, 0.1,,, "Lower energy limit (TeV)"
+emax,          r, a, 100.0,,, "Upper energy limit (TeV)"
+enumbins,      i, a, 20,,, "Number of energy bins"
+bkgmethod,     s, h, REFLECTED,REFLECTED,, "Method for background estimation"
+srcshape,      s, h, CIRCLE,CIRCLE,, "Shape of the source region"
+bkgregmin,     i, h, 2,1,1000, "Minimum number of background regions"
+use_model_bkg, b, h, yes,,, "Use background model?"
+maxoffset,     r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
+etruemin,      r, h, 0.01,,, "Minimum true energy (TeV)"
+etruemax,      r, h, 300.0,,, "Maximum true energy (TeV)"
+etruebins,     i, h, 30,,, "Number of bins per decade for true energy bins"
+calc_ts,       b, h, yes,,, "Compute TS value in each bin"
+calc_ulim,     b, h, yes,,, "Compute upper limit in each bin"
+fix_srcs,      b, h, yes,,, "Fix other skymodel parameters"
+fix_bkg,       b, h, no,,, "Fix background parameters"
 
 #
 # Standard parameters
 #====================
+nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
+publish,  b, h, no,,, "Publish output map on VO Hub?"
 chatter,  i, h, 2,0,4, "Chattiness of output"
 clobber,  b, h, yes,,, "Overwrite existing output files with new output files?"
 debug,    b, h, no,,, "Activate debugging mode?"
diff --git a/cscripts/csscs.py b/cscripts/csscs.py
index 412981e..9aeace7 100755
--- a/cscripts/csscs.py
+++ b/cscripts/csscs.py
@@ -21,6 +21,9 @@
 import sys
 import gammalib
 import ctools
+from cscripts import mputils
+from cscripts import obsutils
+import math
 
 
 # ================= #
@@ -45,6 +48,13 @@ class csscs(ctools.csobservation):
         # Initialise application by calling the base class constructor
         self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
 
+        # Initialise data members
+        self._nthreads = 0
+        self._roisz = 0.
+        self._srcnames = []
+        self._fits = None
+        self._excl_reg_map = None # Exclusion region map for on/off analysis
+
         # Return
         return
 
@@ -59,7 +69,13 @@ class csscs(ctools.csobservation):
             Pickled instance
         """
         # Set pickled dictionary
-        state = {'base' : ctools.csobservation.__getstate__(self)}
+        # Set pickled dictionary
+        state = {'base'     : ctools.csobservation.__getstate__(self),
+                 'nthreads' : self._nthreads,
+                 'roisz'    : self._roisz,
+                 'srcnames' : self._srcnames,
+                 'fits'     : self._fits,
+                 'excl_reg_map' : self._excl_reg_map}
 
         # Return pickled dictionary
         return state
@@ -75,6 +91,11 @@ class csscs(ctools.csobservation):
         """
         # Set state
         ctools.csobservation.__setstate__(self, state['base'])
+        self._nthreads = state['nthreads']
+        self._roisz = state['roisz']
+        self._srcnames = state['srcnames']
+        self._fits = state['fits']
+        self._excl_reg_map = state['excl_reg_map']
 
         # Return
         return
@@ -85,14 +106,497 @@ class csscs(ctools.csobservation):
         """
         Get parameters from parfile
         """
-        # TODO: Add code to query all relevant parameters
+
+        # Set observation if not done before
+        if self.obs().is_empty():
+            self.obs(self._get_observations())
+
+        # Collect number of unbinned, binned and On/Off observations in
+        # observation container
+        n_unbinned = 0
+        n_binned   = 0
+        n_onoff    = 0
+        for obs in self.obs():
+            if obs.classname() == 'GCTAObservation':
+                if obs.eventtype() == 'CountsCube':
+                    n_binned += 1
+                else:
+                    n_unbinned += 1
+            elif obs.classname() == 'GCTAOnOffObservation':
+                n_onoff += 1
+
+        # Check that we have at least one unbinned or binned CTA observation
+        if n_unbinned == 0 and n_binned == 0:
+            msg = 'No unbinned or binned CTA observations found. ' \
+                  'Please provide at least one unbinned or binned ' \
+                  'CTA observation.'
+            raise RuntimeError(msg)
+
+        # Set observation statistic
+        self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
+
+        # Set models if there are none in the container
+        if self.obs().models().size() == 0:
+            self.obs().models(self['inmodel'].filename())
+
+        # Query exclusion region
+        self["inexclusion"].query()
+
+        # Query the map definition parameters
+        self['xref'].real()
+        self['yref'].real()
+        self['coordsys'].string()
+        self['proj'].string()
+        self['nxpix'].integer()
+        self['nypix'].integer()
+
+        # Compute the size of the ROIs for the analysis in deg
+        self._roisz = self['binsz'].real() * self['roisz'].real()
+
+        # Query energy boundaries
+        self["emin"].real()
+        self["emax"].real()
+
+        # If we have unbinned observations query the analysis method
+        if n_unbinned > 0:
+            # If method is ONOFF query some csphagen parameters
+            if self['method'].string() == 'ONOFF':
+                pass
+                # TODO: implement Onoff analysis
+                # self["enumbins"].integer()
+                # self['use_model_bkg'].boolean()
+                # if self["bkgmethod"].string() == "REFLECTED":
+                #     self["srcshape"].string()
+                #     self["bkgregmin"].integer()
+                #     self["maxoffset"].real()
+                #     self["etruemin"].real()
+                #     self["etruemax"].real()
+                #     self["etruebins"].integer()
+
+        # Query target source names
+        srcnames = self['srcnames'].string()
+
+        # Fashion target sources into Python list
+        self._srcnames = srcnames.split(';')
+        # Strip leading and trailing spaces from names
+        for s, name in enumerate(self._srcnames):
+            self._srcnames[s] = name.strip()
+
+        # Verify that the target sources are in the model
+        src_in_model = True
+        for name in self._srcnames:
+            src_in_model = src_in_model and self.obs().models().contains(name)
+        if src_in_model == False:
+            msg = 'Not all target sources are present in input model.'
+            raise RuntimeError(msg)
+
+        # Query the hidden parameters, just in case
+        self['edisp'].boolean()
+        self['calc_ulim'].boolean()
+        self['calc_ts'].boolean()
+        self['fix_bkg'].boolean()
+        self['fix_srcs'].boolean()
+
+        # Read ahead output parameters
+        if self._read_ahead():
+            self['outfile'].query()
 
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
 
+        # Set number of processes for multiprocessing
+        self._nthreads = mputils.nthreads(self)
+
+        # Return
+        return
+
+    def _create_map(self):
+        """
+        Create a gammalib.GSkyMap object to store the output
+
+        Returns
+        -------
+        skymap : `~gammalib.GSkyMap`
+        """
+
+        skymap = gammalib.GSkyMap(self['proj'].string(), self['coordsys'].string(),
+                                  self['xref'].real(), self['yref'].real(),
+                                  -self['binsz'].real(), self['binsz'].real(),
+                                  self['nxpix'].integer(), self['nypix'].integer(),
+                                  1)
+
+        return skymap
+
+    def _adjust_model_pars(self):
+        """
+        Adjust model parameters
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Adjust model parameters')
+
+        # Adjust model parameters dependent on input user parameters
+        for model in self.obs().models():
+
+            # Initialise TS flag for all models to false
+            model.tscalc(False)
+
+            # Log model name
+            self._log_header3(gammalib.EXPLICIT, model.name())
+
+            # Freeze all parameters except the normalization
+            # which is assumed to be the first spectral parameter
+            normpar = model.spectral()[0]
+            for par in model:
+                if par.name() == normpar.name():
+                    pass
+                else:
+                    if par.is_free():
+                        self._log_string(gammalib.EXPLICIT, ' Fixing "' + par.name() + '"')
+                        par.fix()
+
+            # Deal with the target sources
+            if model.name() in self._srcnames:
+
+                # Free the normalisation parameter which is assumed to be
+                # the first spectral parameter
+                if normpar.is_fixed():
+                    self._log_string(gammalib.EXPLICIT, ' Freeing "'+normpar.name()+'"')
+                    normpar.free()
+
+                # Optionally compute Test Statistic value
+                if self['calc_ts'].boolean():
+                    model.tscalc(True)
+
+            # Deal with background models
+            elif self['fix_bkg'].boolean() and not model.classname() == 'GModelSky':
+
+                if normpar.is_free():
+                    self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"')
+                    normpar.fix()
+
+
+            # Deal with background sources
+            elif self['fix_srcs'].boolean() and model.classname() == 'GModelSky':
+
+                if normpar.is_free():
+                    self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"')
+                    normpar.fix()
+
         # Return
         return
 
+    def _mask_cube(self,obs,ra,dec,rad):
+        """
+        Mask cube observation
+
+        Parameters
+        ----------
+        obs     : `gammalib.GCTAObservation` Input observation of type cube
+        ra      : `float` R.A. (deg)
+        dec     : `float` Dec (deg)
+        rad     : `float` radius (deg)
+
+        Returns
+        -------
+        new_obs : `gammalib.GCTAObservation`
+            Output observation with masked cube
+        """
+
+        # Put observation in container
+        container = gammalib.GObservations()
+        container.append(obs)
+
+        # Filter cube according to ROI and user energy range
+        cubemask = ctools.ctcubemask(container)
+        cubemask['ra'] = ra
+        cubemask['dec'] = dec
+        cubemask['rad'] = rad
+        cubemask['regfile'] = 'NONE'
+        cubemask['emin'] = self['emin'].real()
+        cubemask['emax'] = self['emax'].real()
+        cubemask.run()
+
+        # Get deep copy of filtered observation
+        new_obs = cubemask.obs()[0].copy()
+
+        # Return
+        return new_obs
+
+    def _mask_evlist(self,obs,ra,dec,rad):
+        """
+        Mask event list
+
+        Parameters
+        ----------
+        obs     : `gammalib.GCTAObservation` Input observation of type cube
+        ra      : `float` R.A. (deg)
+        dec     : `float` Dec (deg)
+        rad     : `float` radius (deg)
+
+        Returns
+        -------
+        new_obs : `gammalib.GCTAObservation`
+            Output observation with masked event list
+        """
+
+        # Put observation in container
+        container = gammalib.GObservations()
+        container.append(obs)
+
+        # Filter event list according to ROI and user energy range
+        select = ctools.ctselect(container)
+        select['ra'] = ra
+        select['dec'] = dec
+        select['rad'] = rad
+        select['emin'] = self['emin'].real()
+        select['emax'] = self['emax'].real()
+        select['tmin'] = 'INDEF'
+        select['tmax'] = 'INDEF'
+        select.run()
+
+        # Get deep copy of filtered observation
+        new_obs = select.obs()[0].copy()
+
+        # Return
+        return new_obs
+
+    def _mask_onoff(self,obs,ra,dec,rad):
+        """
+        Create On/Off observation
+        with On region matching the required mask
+
+        Parameters
+        ----------
+        obs     : `gammalib.GCTAObservation` Input observation of type cube
+        ra      : `float` R.A. (deg)
+        dec     : `float` Dec (deg)
+        rad     : `float` radius (deg)
+
+        Returns
+        -------
+        new_obs : `gammalib.GCTAOnoffObservation`
+            Output On/Off observation
+        """
+
+        # Put observation in container
+        container = gammalib.GObservations()
+        container.append(obs)
+
+        onoff_obs = obsutils.get_onoff_obs(self, container, ra=ra, dec=dec, rad=rad)
+
+        return onoff_obs[0]
+
+    def _mask_observations(self,ra,dec,rad):
+        """
+        Create observations restricted to circular ROI
+
+        Parameters
+        ----------
+        ra      : `float` R.A. (deg)
+        dec     : `float` Dec (deg)
+        rad     : `float` radius (deg)
+
+        Returns
+        -------
+        new_obs : `~gammalib.GObservations`
+            Observations in circular ROI
+        """
+
+        # Create output observation container
+        new_obs = gammalib.GObservations()
+
+        # Loop over input observations
+        for obs in self.obs():
+            if obs.classname() == 'GCTAObservation':
+                if obs.eventtype() == 'CountsCube':
+                    masked_obs = self._mask_cube(obs,ra,dec,rad)
+                else:
+                    if self['method'].string() == '3D':
+                        masked_obs = self._mask_evlist(obs,ra,dec,rad)
+                    elif self['method'].string() == 'ONOFF':
+                        pass
+                        # TODO: implement Onoff analysis
+                        # masked_obs = self._mask_onoff(obs, ra, dec, rad)
+                new_obs.append(masked_obs)
+            # Skip On-Off and non-CTA observations
+            else:
+                pass
+
+        # Append models to output observations
+        new_obs.models(self.obs().models())
+
+        # Return
+        return new_obs
+
+    def _pixel_analysis(self,ra,dec):
+        """
+        Performs analysis over the region of interest
+        corresponding to a single pixel of output map
+
+        Parameters
+        ----------
+        inx     : `int` pixel index
+        ra      : `float` R.A. (deg)
+        dec     : `float` Dec (deg)
+
+        Returns
+        -------
+        result  : `dict` Results
+        """
+
+        # Write header for spatial bin
+        msg = 'Spatial bin centred on (R.A.,Dec) = (%f,%f) deg' % (ra,dec)
+        self._log_header2(gammalib.EXPLICIT, msg)
+
+        # Initialize results
+        result = {}
+        for name in self._srcnames:
+            result[name] = {'flux'    : 0.0,
+                            'flux_err': 0.0,
+                            'TS'      : 0.0,
+                            'ulimit'  : 0.0}
+
+        # Mask observations
+        self._log_header3(gammalib.EXPLICIT, 'Masking observations')
+        masked_obs = self._mask_observations(ra,dec,self._roisz)
+
+        # Set up likelihood analysis
+        self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin')
+        like = ctools.ctlike(masked_obs)
+        like['edisp']    = self['edisp'].boolean()
+        like['nthreads'] = 1  # Avoids OpenMP conflict
+
+        # If chatter level is verbose and debugging is requested then
+        # switch also on the debug model in ctlike
+        if self._logVerbose() and self._logDebug():
+            like['debug'] = True
+
+        # Perform maximum likelihood fit
+        like.run()
+
+        # Write model results for explicit chatter level
+        self._log_string(gammalib.EXPLICIT, str(like.obs().models()))
+
+        # Prepare objects for flux extraction
+        # ROI
+        centre = gammalib.GSkyDir()
+        centre.radec_deg(ra,dec)
+        roi = gammalib.GSkyRegionCircle(centre,self._roisz)
+        # Energy boundaries
+        emin = gammalib.GEnergy(self['emin'].real(),'TeV')
+        emax = gammalib.GEnergy(self['emax'].real(), 'TeV')
+
+        # Continue only if log-likelihood is non-zero
+        if like.obs().logL() != 0.0:
+
+            # Loop over target sources
+            for name in self._srcnames:
+
+                # Get source
+                source = like.obs().models()[name]
+
+                # Get flux
+                # Integrate spectral model between emin and emax
+                flux = source.spectral().flux(emin,emax)
+                # Multiply by flux from spatial component in ROI
+                flux *= source.spatial().flux(roi)
+                # Divide by solid angle of ROI
+                flux /= 1 - math.cos(math.radians(self._roisz))
+                result[name]['flux'] = flux
+
+                # Get flux error
+                # Normalization parameter
+                normpar = source.spectral()[0]
+                # Only normalization parameter free
+                # Relative error on flux is same as on normalization
+                # Avoid zero division error
+                if normpar.value() > 0.:
+                    flux_error = flux * normpar.error()/normpar.value()
+                else:
+                    flux_error = 0.
+                result[name]['flux_error'] = flux_error
+
+                # If requested get TS
+                if self['calc_ts'].boolean():
+                    result[name]['TS'] = source.ts()
+
+                # Upper limit
+                # TODO
+
+        else:
+            value = 'Likelihood is zero. Bin is skipped.'
+            self._log_value(gammalib.TERSE, '(R.A.,Dec) = (%f,%f) deg' % (ra,dec), value)
+
+        # Return result
+        return result
+
+    def _fill_fits(self,results):
+        """
+        Fill FITS object to store the results
+
+        Parameters
+        ----------
+        result  : `dict` Results
+
+        Returns
+        -------
+        fits    : `~gammalib.GFits` Fits object
+        """
+
+        # Create new Fits object
+        fits = gammalib.GFits()
+
+        # If more than one source of interest
+        # append empty primary HDU
+        if len(self._srcnames) > 1:
+            fits.append(gammalib.GFitsImageDouble())
+
+        # Prepare some header cards
+        # Flux units
+        fcard = gammalib.GFitsHeaderCard('BUNIT', 'ph/cm2/s/sr' , 'Photon flux')
+        txt = 'Correlation radius %f deg' %(self._roisz)
+        roiszcard = gammalib.GFitsHeaderCard('COMMENT', '' , txt)
+
+        # Loop over target sources
+        for s, name in enumerate(self._srcnames):
+
+            # Create minimal set of skymaps to store fluxes and flux errors
+            fmap = self._create_map()
+            errmap = self._create_map()
+
+            #Create map for TS
+            if self['calc_ts'].boolean():
+                tsmap = self._create_map()
+
+            # Loop over pixels and fill maps
+            for inx in range(fmap.npix()):
+
+                # Fill flux and flux error
+                fmap[inx,0] = results[inx][name]['flux']
+                errmap[inx, 0] = results[inx][name]['flux_error']
+
+                # If requested fill TS
+                if self['calc_ts'].boolean():
+                    tsmap[inx, 0] = results[inx][name]['TS']
+
+            # Write maps to Fits
+            fhdu = fmap.write(fits, name + ' FLUX')
+            errhdu = errmap.write(fits, name + ' FLUX ERROR')
+
+            # Add units and comments to headers
+            fhdu.header().append(fcard)
+            errhdu.header().append(fcard)
+            fhdu.header().append(roiszcard)
+            errhdu.header().append(roiszcard)
+
+            # If requested write TS map to fits
+            if self['calc_ts'].boolean():
+                tshdu = tsmap.write(fits, name + ' TS')
+                tshdu.header().append(roiszcard)
+
+        # Return
+        return fits
+
 
     # Public methods
     def run(self):
@@ -106,23 +610,77 @@ class csscs(ctools.csobservation):
         # Get parameters
         self._get_parameters()
 
-        # TODO: Your code goes here
+        # Write input observation container into logger
+        self._log_observations(gammalib.NORMAL, self.obs(), 'Input observation')
+
+        # Create template for output map, to use to extract pixel values
+        outmap = self._create_map()
+
+        # Adjust model parameters dependent on input user parameters
+        self._adjust_model_pars()
+
+        # Loop over pixels and extract results
+        results = {}
+        for inx in range(outmap.npix()):
+
+            # Determine pixel centre
+            pixdir = outmap.inx2dir(inx)
+            ra = pixdir.ra_deg()
+            dec = pixdir.dec_deg()
+
+            # Analyse ROI
+            result = self._pixel_analysis(ra,dec)
+            results[inx] = result
+
+        # Fill output Fits
+        self._fits = self._fill_fits(results)
 
         # Return
         return
 
     def save(self):
         """
-        Save something
+        Save maps to Fits flie
         """
         # Write header
-        self._log_header1(gammalib.TERSE, 'Save something')
+        self._log_header1(gammalib.TERSE, 'Save maps')
+
+        # Continue only if FITS file is valid
+        if self._fits != None:
+
+            # Get outmap parameter
+            outfile = self['outfile'].filename()
 
-        # TODO: Your code goes here
+            # Log file name
+            self._log_value(gammalib.NORMAL, 'Output file', outfile.url())
+
+            # Save spectrum
+            self._fits.saveto(outfile, self['clobber'].boolean())
 
         # Return
         return
 
+    def exclusion_map(self, object=None):
+        """
+        Return and optionally set the exclusion regions map
+
+        Parameters
+        ----------
+        object : `~gammalib.GSkyRegion` or `~gammalib.GSkyMap` or `~gammalib.GFilename`
+            Exclusion regions object
+
+        Returns
+        -------
+        region : `~gammalib.GSkyRegionMap`
+            Exclusion regions map
+        """
+        # If a regions object is provided then set the regions ...
+        if object is not None:
+            self._excl_reg_map = gammalib.GSkyRegionMap(object)
+
+        # Return
+        return self._excl_reg_map
+
 
 # ======================== #
 # Main routine entry point #