Merge branch 'devel' of https://cta-gitlab.irap.omp.eu/ctools/ctools into devel

ldivenere
2 parents 6acd8596 03b31792
Showing 29 changed files with 1306 additions and 724 deletions
ChangeLog
NEWS
README.md
caldb/data/cta/prod2/bcf/North_0.5h/irf_file.fits.gz
caldb/data/cta/prod2/bcf/North_50h/irf_file.fits.gz
caldb/data/cta/prod2/bcf/North_5h/irf_file.fits.gz
caldb/data/cta/prod2/bcf/South_0.5h/irf_file.fits.gz
caldb/data/cta/prod2/bcf/South_50h/irf_file.fits.gz
caldb/data/cta/prod2/bcf/South_5h/irf_file.fits.gz
cscripts/csobsinfo.py
cscripts/csphagen.py
cscripts/csresspec.py
cscripts/csroot2caldb.py
doc/source/admin/install_devel.rst
doc/source/admin/release_history/1.5.1.rst
doc/source/admin/release_history/1.6.rst
doc/source/admin/release_history/index.rst
doc/source/users/reference_manual/ctskymap.rst
pyext/ctobservation.i
src/ctskymap/ctskymap.cpp
-2018-01-25
+2018-02-14
+
+        * Version 1.6.0 released
+        ========================
+
+        Allow specifying of input sky map in ctskymap (#2311)
+        Add ras() and decs() methods to csobsinfo (#2313)
+        Accelerate IRF and RING background methods in ctskymap (#2309)
+
+
+2018-02-16
+
+        * Bug fix version 1.5.1 released
+        ================================
+
+        Take energy dispersion in csresspec modelling into account (#2341)
+        Enable using csphagen for HESS, VERITAS and MAGIC observations (#2307)
+
         * Version 1.5.0 released
         ========================
-25 January 2018
+New Features and Important Changes in ctools 1.6.0
+
+14 February 2018
 Introduction
@@ -9,10 +11,7 @@ ctools release with respect to the precedent release.
 Configuration
 -------------
-In case of an error, test output is now written into log file (#2264).
-
-The cscripts::execute() method no longer creates log files, making the
-cscripts behaviour compliant with that of ctools (#2182).
+None
 Python interface
@@ -27,8 +26,7 @@ None
 Base classes
 ------------
-Add ctools::time_reference constant and remove ctools::time_reference() method
-(#2295).
+None
 ctbin - CTA event binning
@@ -43,9 +41,7 @@ None
 ctbutterfly - Generate butterfly diagrams
 -----------------------------------------
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
-
-Add butterfly() method to return the butterfly diagram as CSV object (#1512).
+None
 ctcubemask - CTA cube filtering
@@ -60,7 +56,7 @@ None
 cterror - Likelihood profile error computation
 ----------------------------------------------
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
+None
 ctexpcube - CTA exposure cube generation
@@ -70,7 +66,7 @@ None
 ctlike - Maximum likelihood fitting
 -----------------------------------
-Rename "stats" parameter to "statistic" (#2237).
+None
 ctmapcube - CTA map cube generation
@@ -80,22 +76,12 @@ None
 ctmodel - CTA counts cube modelling
 -----------------------------------
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
-
-The tool now also checks which models actually overlap with the RoI of a given
-observation and then only use these models to compute the model cube. The
-parallelization of the tool has also been optimised so that diffuse models are
-not loaded over and over again (#2189).
+None
 ctobssim - CTA observation simulation
 -------------------------------------
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
-
-Set the OBS_ID keyword to the number that is used for creating the FITS file
-(#1836).
+None
 ctphase - Event phase computation
@@ -110,51 +96,47 @@ None
 ctpsfcube - CTA point spread function cube generation
 -----------------------------------------------------
-The chatter level is now correctly forwarded to the PSF cube fill method (#2262).
+None
 ctselect - CTA event selection
 ------------------------------
-Any .gz suffix is now correctly stripped from the event filenames when creating
-the output filenames (#2287).
-
-ctselect now informs about non-overlapping time interval; before the time interval
-was set to [0,0] which was disturbing and eventually incorrect (#2272).
-
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
 ctskymap - CTA sky mapping
 --------------------------
-Header keywords informing about the dates and times covered by the observations
-used to generate the skymap, as well as the event energy range, are now appended
-as header keywords to the FITS file (#2194).
-
-Add the RING background method, which is a classical analysis methods for IACT
-data analysis, for background subtraction (#2207).
+Add the "inmap", "iterations" and "threshold" parameters. The hidden
+"inmap" parameter allows specifying a sky map file containing a "COUNTS"
+and "ACCEPTANCE" extension which avoids looping over event lists in an
+observation container. The "iterations" parameter allows the iterative
+determination of an exclusion map from the resulting significance map for
+the RING background method. The exclusion map will flag all pixels with
+significances above the "threshold" value so that they will be excluded
+in the subsequent RING background determination. Since this will impact
+the pixel significances, the procedure needs to be repeated a couple of
+times to stabilise (#2311).
+
+Accelerate the IRF and RING background methods. For the RING background
+method the computation is either performed using a FFT or using a direct
+computation. The FFT computation is faster, but assumed Euclidian
+distances for the sky map grid, the direct computation is slower, but
+accurately computes the distances between sky map pixels (#2309).
 cttsmap - Generation of Test Statistic map
 ------------------------------------------
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
+None
 ctulimit - Compute upper limit
 ------------------------------
-The tool now also computes correct upper flux limits for diffuse map cubes
-for which the spectral component is a normalization of the spectral information
-available in the map cube (#2298).
-
-Add "statistic" parameter to allow setting of optimizer statistic (#2237).
+None
 support - Support classes
 -------------------------
-Change DATE_OBS, DATE_END, TIME_OBS and TIME_END FITS keywords into DATA-OBS,
-DATE-END, TIME-OBS and TIME-END.
-
-Add ctobservation.obs() method to set an observation container.
+None
 cscaldb - Inspect calibration database
@@ -189,16 +171,7 @@ None
 cslightcrv - Generate light curve
 ---------------------------------
-The script now also supports the handling of On/Off observations. For that
-purpose parameters were added that are needed for the execution of the
-csphagen script (#2273).
-
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the light
-curve generation (#2237).
-
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
 csmodelinfo - Return model information
@@ -218,8 +191,7 @@ None
 csmodelsois - Puts selected sources in diffuse model cube
 ---------------------------------------------------------
-Add script that replaces selected models in a model container by a diffuse
-map cube (#2180).
+None
 csobs2caldb - Generate calibration database entry for an IACT observation
@@ -234,47 +206,27 @@ None
 csobsinfo - Return observation definition information
 -----------------------------------------------------
-None
+Add ras() and decs() methods (#2313)
 csobsselect - Select observations
 ---------------------------------
-Add tmin and tmax parameters to allow selection of observation by time (#2297).
+None
 csphagen - Generate On/Off observations
 ---------------------------------------
-The CUSTOM background method was added to allow specifying On and Off regions
-from a region file or FITS WCS map. For multiple observations, a custom Off
-region per observation can be specified in the observation definition XML file
-using the OffRegions parameter (#2284).
-
-The script was added to allow the generation of PHA, ARF and RMF datasets for classical On/Off IACT analysis (#2208).
+None
 csphasecrv - Generate phase curve
 ---------------------------------
-The script now also supports the handling of On/Off observations. For that
-purpose parameters were added that are needed for the execution of the
-csphagen script (#2273).
-
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the light
-curve generation (#2237).
+None
 cspull - Generate pull distributions
 ------------------------------------
-The script now handles On/Off observation. The parameters "onsrc" and "onrad"
-were added to allow specifying the On region test source and On region radius
-(#2260).
-
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the pull
-distribution generation (#2237).
-
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
 csresmap - Generate residual map
@@ -284,7 +236,7 @@ None
 csresspec - Generate residual spectrum
 --------------------------------------
-Add script to generate a residual spectrum (#2289)
+None
 csroot2caldb - Generates calibration database from ROOT file
@@ -294,25 +246,12 @@ None
 cssens - Computes sensitivity curve
 -----------------------------------
-The script now derived from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the sensitivity
-computation (#2237).
+None
 csspec - Generates Spectral Energy Distribution
 -----------------------------------------------
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the spectrum
-generation (#2237).
-
-The script now can also handle On/Off spectra for SED generation (#2231).
-
-The script now can also handle non-CTA data. For that purpose a "method"
-parameter has been introduced that allows to selected between slicing of the data
-(method="SLICE"), which is only supported for CTA observations, and replacing of
-the spectral model by a node function (method="NODES"). If method="AUTO" the tool
-will use "SLICE" for observations containing only CTA data and "NODES" for observations
-containing other data (#1927).
+None
 cssrcdetect - Detects sources in sky map
@@ -322,18 +261,7 @@ None
 cstsdist - Computes Test Statistic distribution
 -----------------------------------------------
-The script now derives from the csobservation base class. A "statistic"
-parameter was added to allow setting the optimizer statistic for the Test
-Statistic distribution generation (#2237).
-
-The parameter interface was considerably revised, in particular the following
-parameters were removed: "ra", "dec", "rad", "coordsys", "proj", "emin",
-"emax", "enumbins", "tmin", "tmax", "npix" and "binsz". The "inobs" parameter
-is now always queried. The script now reads in an observation or an observation
-container and derives the TS distribution for a source based on that observation
-container. So far the script supports unbinned, binned and stacked observations.
-If an observation definition XML file is specified the script assumes that all
-observations are event lists (#1567).
+None
 cstsmapmerge - Merge Test Statistic maps
@@ -343,14 +271,12 @@ None
 cstsmapsplit - Prepare split Test Statistic computation
 -------------------------------------------------------
-A "statistic" parameter was added to allow setting the optimizer statistic for
-the TS map generation (#2237).
+None
 csviscube - Generate visibility cube
 ------------------------------------
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+None
 csworkflow - Execute workflow
 ctools information
 ==================
-* Version:             1.6.0.dev1 (1 August 2017)
+* Version:             1.6.0.dev1
 * Author:              Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
 * GammaLib dependency: 1.6.0.dev1
@@ -59,11 +59,11 @@ class csobsinfo(ctools.csobservation):
         """
         Get parameters from parfile and setup the observation
         """
-        # Get parameters   
+        # Get parameters
         if self.obs().size() == 0:
             self._require_inobs('csobsinfo::get_parameters')
             self.obs(self._get_observations(False))
-        
+
         # Initialise object position
         self._obj_dir = gammalib.GSkyDir()
@@ -83,7 +83,7 @@ class csobsinfo(ctools.csobservation):
         # Return
         return
-        
+
     # Public methods
     def run(self):
@@ -96,7 +96,7 @@ class csobsinfo(ctools.csobservation):
         # Get parameters
         self._get_parameters()
-        
+
         # Initialise arrays to store certain values for reuse
         # Todo, think about using a python dictionary
         self._offsets  = []
@@ -107,7 +107,7 @@ class csobsinfo(ctools.csobservation):
         self._ebounds  = gammalib.GEbounds()
         self._gti      = gammalib.GGti()
         obs_names      = []
-        
+
         # Initialise output to be filled
         ontime         = 0.0
         livetime       = 0.0
@@ -115,12 +115,12 @@ class csobsinfo(ctools.csobservation):
         n_eventbins    = 0
         n_obs_binned   = 0
         n_obs_unbinned = 0
-        
+
         # Write header
         if self._logTerse():
             self._log('\n')
             self._log.header1(gammalib.number('Observation', self.obs().size()))
-        
+
         # Loop over observations
         for obs in self.obs():
@@ -128,44 +128,44 @@ class csobsinfo(ctools.csobservation):
             if not obs.classname() == 'GCTAObservation':
                 self._log('Skipping '+obs.instrument()+' observation\n')
                 continue 
-            
+
             # Use observed object as observation name if name is not given
             obs_name = obs.name()
             if obs_name == '':
                 obs_name = obs.object()
-            
+
             # Logging
             if self._logExplicit():
                 self._log.header2(obs_name)
-                
+
             # Retrieve observation name
             obs_names.append(obs_name)
-            
+
             # Retrieve energy boundaries
             obs_bounds = obs.events().ebounds()
-            
+
             # Retrieve time interval
             obs_gti = obs.events().gti()
-            
+
             # Compute mean time and dead time fraction in percent
             deadfrac = (1.0-obs.deadc())*100.0
-            
+
             # Retrieve pointing and store Ra,Dec
-            pnt_dir = obs.pointing().dir()  
+            pnt_dir = obs.pointing().dir()
             self._pnt_ra.append(pnt_dir.ra_deg())
-            self._pnt_dec.append(pnt_dir.dec_deg())       
-            
+            self._pnt_dec.append(pnt_dir.dec_deg())
+
             # If avaliable append energy boundaries
             if obs_bounds.size() > 0 : 
                 self._ebounds.append(obs_bounds.emin(),obs_bounds.emax())
-            
+
             # Append time interval
             self._gti.append(obs_gti.tstart(), obs_gti.tstop())
-            
+
             # Increment global livetime and ontime
             ontime   += obs.ontime()
             livetime += obs.livetime()
-            
+
             # Bookkeeping
             if obs.eventtype() == 'CountsCube':
                 n_eventbins  += obs.events().size()
@@ -190,10 +190,10 @@ class csobsinfo(ctools.csobservation):
             if self._compute_offset:
                 offset = pnt_dir.dist_deg(self._obj_dir)
                 self._offsets.append(offset)
-                    
+
             # Optionally log details
             if self._logExplicit():
-                
+
                 # Log the observation energy range (if available)
                 self._log.parformat('Energy range')
                 if obs_bounds.size() == 0:
@@ -203,7 +203,7 @@ class csobsinfo(ctools.csobservation):
                     self._log(' - ')
                     self._log(str(obs_bounds.emax()))
                 self._log('\n')
-                
+
                 # Log observation time interval
                 self._log.parformat('Time range (MJD)')
                 if obs_gti.size() == 0:
@@ -213,7 +213,7 @@ class csobsinfo(ctools.csobservation):
                     self._log(' - ')
                     self._log(str(obs_gti.tstop().mjd()))
                 self._log('\n')
-                
+
             # Log observation information
             self._log_value(gammalib.EXPLICIT, 'Ontime', '%.3f s' %
                             obs.ontime())
@@ -233,7 +233,7 @@ class csobsinfo(ctools.csobservation):
                             zenith)
             self._log_value(gammalib.EXPLICIT, 'Azimuth angle', '%.2f deg' %
                             azimuth)
-                    
+
         # Write summary header
         self._log_header1(gammalib.NORMAL, 'Summary')
@@ -258,7 +258,7 @@ class csobsinfo(ctools.csobservation):
             mean_azimuth = '%.2f deg' % (sum(self._azimuths) / len(self._azimuths))
         else:
             mean_azimuth = 'Unknown'
-        
+
         # Log mean offset, azimuth and zenith angle
         self._log_header3(gammalib.NORMAL, 'Pointings')
         self._log_value(gammalib.NORMAL, 'Mean offset angle',  mean_offset)
@@ -323,48 +323,60 @@ class csobsinfo(ctools.csobservation):
         # Check if DS9 file is valid
         if ds9file.url() != 'NONE':
-            
+
             # Write header
             self._log_header1(gammalib.TERSE, 'Save pointings in DS9 file')
             # Log filename
             self._log_value(gammalib.NORMAL, 'DS9 filename', ds9file.url())
-            
-            # Open file   
+
+            # Open file
             f = open(ds9file.url(),'w')
-            
+
             # Write coordinate system
             f.write('fk5\n')
-            
+
             # Loop over pointings
             for i in range(len(self._pnt_ra)):
-                
+
                 # Create string
                 line  = 'point('
                 line += str(self._pnt_ra[i])+','+str(self._pnt_dec[i])+')'
                 line += ' # point=cross 20 width=3\n'
-                
+
                 # Write to file
                 f.write(line)
-                
+
             # Close file
             f.close()
-        
+
         # Return
-        return 
-            
+        return
+
     def zeniths(self):
         """
         Return zenith angles
         """
         return self._zeniths
-    
+
     def azimuths(self):
         """
         Return azimuth angles
         """
         return self._azimuths
-    
+
+    def ras(self):
+        """
+        Return pointings right ascension
+        """
+        return self._pnt_ra
+
+    def decs(self):
+        """
+        Return pointings declination
+        """
+        return self._pnt_dec
+
     def offsets(self):
         """
         Return offset angles
@@ -376,13 +388,13 @@ class csobsinfo(ctools.csobservation):
         Return energy boundaries
         """
         return self._ebounds
-    
+
     def gti(self):
         """
         Return good time intervals
         """
         return self._gti
-        
+
 # ======================== #
 # Main routine entry point #
@@ -391,6 +403,6 @@ if __name__ == &#39;__main__&#39;:
     # Create instance of application
     app = csobsinfo(sys.argv)
-    
+
     # Execute application
     app.execute()
@@ -123,7 +123,8 @@ class csphagen(ctools.csobservation):
         # Off region region the parameter 'bkgregfile' for all CTA observations
         # without Off region
         for obs in self.obs():
-            if obs.instrument() == 'CTA':
+            #if obs.instrument() == 'CTA':
+            if obs.classname() == 'GCTAObservation':
                 if obs.off_regions().is_empty():
                     filename = self['bkgregfile'].filename()
                     regions  = gammalib.GSkyRegions(filename)
@@ -263,8 +264,8 @@ class csphagen(ctools.csobservation):
         """
         Set models in observation container
-        The method replaces all "CTA" background models by "CTAOnOff"
-        background models.
+        The method replaces all "CTA", "HESS", "VERITAS", and "MAGIC"
+        background models by "CTAOnOff" background models.
         Parameters
         ----------
@@ -279,11 +280,15 @@ class csphagen(ctools.csobservation):
         # Initialise model container
         models = gammalib.GModels()
-        # Loop over all models and replace CTA background model by CTAOnOff
-        # background model
+        # Loop over all models and replace CTA/HESS/VERITAS/MAGIC background
+        # models by CTAOnOff background model
         for model in self.obs().models():
-            if 'GCTA' in model.classname() and 'CTA' in model.instruments():
-                model.instruments('CTAOnOff')
+            if 'GCTA' in model.classname():
+                if 'CTA'     in model.instruments() or \
+                   'HESS'    in model.instruments() or \
+                   'VERITAS' in model.instruments() or \
+                   'MAGIC'   in model.instruments():
+                    model.instruments('CTAOnOff')
             models.append(model)
         # Append model to observation container
@@ -399,7 +404,8 @@ class csphagen(ctools.csobservation):
         for obs in self.obs():
             # Skip non CTA observations
-            if obs.instrument() != 'CTA':
+            #if obs.instrument() != 'CTA':
+            if obs.classname() != 'GCTAObservation':
                 self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \
                                  (obs.instrument(), obs.id()))
                 continue
@@ -241,14 +241,14 @@ class csresspec(ctools.csobservation):
         # Bin events
         cntcube = ctools.ctbin(obs)
-        cntcube['xref'] = ra
-        cntcube['yref'] = dec
-        cntcube['binsz'] = 0.02
-        cntcube['nxpix'] = npix
-        cntcube['nypix'] = npix
-        cntcube['proj'] = 'TAN'
+        cntcube['xref']     = ra
+        cntcube['yref']     = dec
+        cntcube['binsz']    = 0.02
+        cntcube['nxpix']    = npix
+        cntcube['nypix']    = npix
+        cntcube['proj']     = 'TAN'
         cntcube['coordsys'] = 'CEL'
-        cntcube['ebinalg'] = self['ebinalg'].string()
+        cntcube['ebinalg']  = self['ebinalg'].string()
         if self['ebinalg'].string() == 'FILE':
             cntcube['ebinfile'] = self['ebinfile'].filename()
         else:
@@ -356,7 +356,7 @@ class csresspec(ctools.csobservation):
     def _cube_to_spectrum(self, cube, evlist_info):
         """
         Derive from event cube a count spectrum
-        
+
         If data come from event list use only the ROI and energy range of
         the original data. Apply user defined mask if requested.
@@ -550,29 +550,30 @@ class csresspec(ctools.csobservation):
             # if 3D observation
             if obs[0].classname() == 'GCTAObservation':
-                # Prepare Observations
-
                 # If already binned set the evlist_info dictionary to have
                 # attribute was_list False
                 if obs[0].eventtype() == 'CountsCube':
                     evlist_info = {'was_list': False}
-                # Otherwise bin now
+
+                # ... otherwise bin now
                 else:
-                    # we remember if we binned an event list
-                    # so that we can mask only the ROI for residual calculation
+                    # we remember if we binned an event list so that we can
+                    # mask only the ROI for residual calculation
                     msg = 'Setting up binned observation'
                     self._log_string(gammalib.NORMAL, msg)
                     obs, evlist_info = self._bin_evlist(obs)
-                # Calculate Model and residuals
-                # If model cube is provided load it
+                # Calculate Model and residuals. If model cube is provided load
+                # it
                 if self._use_maps:
                     modcube = gammalib.GCTAEventCube(self['inmodel'].filename())
-                # Otherwise calculate it now
+
+                # ... otherwise calculate it now
                 else:
                     msg = 'Computing model cube'
                     self._log_string(gammalib.NORMAL, msg)
                     modelcube = ctools.ctmodel(obs)
+                    modelcube['edisp'] = self['edisp'].boolean()
                     modelcube.run()
                     modcube = modelcube.cube().copy()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generate IRFs in CALDB format from a ROOT offaxis performance file
 #
-# Copyright (C) 2016-2017 Juergen Knoedlseder
+# Copyright (C) 2016-2018 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
@@ -618,7 +618,7 @@ class csroot2caldb(ctools.cscript):
             ds['HDU_EA'].card('HI_THRES', emin, '[TeV] High energy threshold')
         # Create "EFFAREA" data column
-        self._make_2D(etrue, ds['HDU_EA'], 'EFFAREA', 'm2')
+        self._make_2D(etrue, ds['HDU_EA'], 'EFFAREA', 'm**2')
         # Return
         return
@@ -745,7 +745,6 @@ class csroot2caldb(ctools.cscript):
                           'Generate Gaussian point spread function extension')
         # Get relevant ROOT histograms
-        #r68 = tfile.Get('AngRes_offaxis')
         r68, _ = self._root2psf_get_r68_and_r80(tfile)
         # Extract number of bins in histogram
@@ -826,8 +825,6 @@ class csroot2caldb(ctools.cscript):
                           'Generate King point spread function extension')
         # Get relevant ROOT histograms
-        #r68 = tfile.Get('AngRes_offaxis')
-        #r80 = tfile.Get('AngRes80_offaxis')
         r68, r80 = self._root2psf_get_r68_and_r80(tfile)
         # Extract number of bins in histogram
@@ -1034,7 +1031,6 @@ class csroot2caldb(ctools.cscript):
         # Get relevant ROOT histograms
         matrix = self._root2edisp_get_migra(tfile)
-        #matrix = tfile.Get('EestOverEtrue_offaxis')
         # Write boundary keywords
         self._set_cif_keywords(ds['HDU_EDISP'], irf['EDISP_NAME'],
@@ -1087,7 +1083,7 @@ class csroot2caldb(ctools.cscript):
                                irf['BGD_BOUNDS'], irf['BGD_DESC'], irf)
         # Create "BGD" data column
-        self._make_3D(array, ds['HDU_BGD'], 'BGD', '1/s/MeV/sr')
+        self._make_3D(array, ds['HDU_BGD'], 'BGD', '1/(MeV s sr)')
         # Return
         return
@@ -1327,7 +1323,7 @@ class csroot2caldb(ctools.cscript):
                             # Log power law replacement
                             self._log_value(gammalib.NORMAL,
                                             'Power law replacement',
-                                            '%e cts/s/deg2 (E=%8.4f TeV (%d) '
+                                            '%e cts/(s deg**2) (E=%8.4f TeV (%d) '
                                             'Off=%.2f deg (%d)' %
                                             (value,engs[ieng],ieng,offset,ioff))
@@ -1479,7 +1475,7 @@ class csroot2caldb(ctools.cscript):
                         # Log background infill
                         self._log_value(gammalib.NORMAL,
                                         'Background infill',
-                                        '%e cts/s/deg2 (E=%8.4f TeV (%d) '
+                                        '%e cts/(s deg**2) (E=%8.4f TeV (%d) '
                                         'Off=%.2f deg (%d)' %
                                         (value,engs[ieng],ieng,offset,ioff))
@@ -7,7 +7,7 @@ The current development version of the code can be downloaded as source
 tarballs, Mac OS X binaries, or directly from the `Git <https://git-scm.com/>`_
 repository.
-The current ctools development release is ``ctools-1.5.0.dev1``.
+The current ctools development release is ``ctools-1.6.0.dev1``.
 Source tarballs
@@ -15,8 +15,8 @@ Source tarballs
 Download the source tarballs from the following links
-- `ctools <http://cta.irap.omp.eu/ctools/releases/ctools/ctools-1.5.0.dev1.tar.gz>`_
-- `GammaLib <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.5.0.dev1.tar.gz>`_
+- `ctools <http://cta.irap.omp.eu/ctools/releases/ctools/ctools-1.6.0.dev1.tar.gz>`_
+- `GammaLib <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.6.0.dev1.tar.gz>`_
 and follow the instructions on :ref:`sec_install_source`.
@@ -26,7 +26,7 @@ Binary packages
 Download the installer image from the following link
-- `Mac OS X (10.7+) <http://cta.irap.omp.eu/ctools/releases/ctools/ctools-1.5.0.dev1-macosx10.7.dmg>`_
+- `Mac OS X (10.7+) <http://cta.irap.omp.eu/ctools/releases/ctools/ctools-1.6.0.dev1-macosx10.7.dmg>`_
 and follow the instructions on :ref:`sec_install_binary`.
+.. _1.5.1:
+
+ctools 1.5.1 release
+====================
+
+Overview
+--------
+
+ctools 1.5.1 is a bug fix release for ctools release 1.5.
+
+The following bugs were fixed:
+
+* [`2341 <https://cta-redmine.irap.omp.eu/issues/2341>`_] -
+  Take energy dispersion in :ref:`csresspec` modelling into account
+
+
+The following improvements were made:
+
+* [`2307 <https://cta-redmine.irap.omp.eu/issues/2307>`_] -
+  Enable using :ref:`csphagen` for HESS, VERITAS and MAGIC observations
+.. _1.6:
+
+ctools 1.6 release
+==================
+
+Overview
+--------
+
+ctools 1.6 is a major release that adds significant functionality.
+
+In particular, this release provides:
+
+* ...
+
+
+Bug fixes
+---------
+
+* ...
+
+
+Improvements
+------------
+
+* [`2313 <https://cta-redmine.irap.omp.eu/issues/2313>`_] -
+  Add `ras()` and `decs()` methods to :ref:`csobsinfo`
+* [`2311 <https://cta-redmine.irap.omp.eu/issues/2311>`_] -
+  Allow specifying of input skymap to :ref:`ctskymap` tool
+* [`2309 <https://cta-redmine.irap.omp.eu/issues/2309>`_] -
+  Accelerate `IRF` and `RING` background methods in :ref:`ctskymap`
+
+
+New features
+------------
+
+* ...
+
+
+Clean up
+--------
+
+* ...
@@ -6,6 +6,8 @@ Release History
 .. toctree::
    :maxdepth: 1
+   ctools 1.6 <1.6>
+   ctools 1.5.1 <1.5.1>
    ctools 1.5 <1.5>
    ctools 1.4.3 <1.4.3>
    ctools 1.4.2 <1.4.2>
@@ -41,6 +41,11 @@ General parameters
     Instrument response function (only required for IRF background subtraction
     if no response information is provided by ``inobs``).
+``(inmap = NONE) [file]``
+    Input sky map file containing a ``COUNTS`` and an ``ACCEPTANCE`` extension.
+    Such extensions are produce by a previous run of the :ref:`ctskymap` tool with
+    the ``IRF`` or ``RING`` background method.
+
 ``outmap [file]``
     Output sky map file.
@@ -82,17 +87,37 @@ Background subtraction configuration parameters
     Background subtraction method.
 ``roiradius [real]``
-    Source region radius for estimating on-counts (in degrees).
+    Source region radius for ``RING`` subtraction (in degrees).
 ``inradius [real]``
-    Inner ring radius for RING subtraction (in degrees).
+    Inner background ring radius for ``RING`` subtraction (in degrees).
 ``outradius [real]``
-    Outer ring radius for RING subtraction (in degrees).
+    Outer background ring radius for ``RING`` subtraction (in degrees).
+
+``iterations [integer]``
+    Number of iterations for the automatic computation of exclusion regions.
+    Sky map pixels with detection significance above the value specified by the
+    ``threshold`` parameter will be defined as exclusion regions, and the ``RING``
+    background will be recomputed based on these exclusion regions. Since this
+    will likely change the significance of the pixels, the procedure can be
+    iteratively repeated, and the ``iterations`` parameter specifies how often
+    the procedure should be repeated.
+
+``threshold [real]``
+    Significance threshold above which pixels should be excluded for the
+    background computation in the ``RING`` background method.
 ``(inexclusion = NONE) [file]``
     Exclusion region file as either a FITS map or DS9 region file.
+``(usefft = yes) [boolean]``
+    Specifies whether a FFT should be used for ring computation. If set to
+    ``no`` the ring background will be computed directly. Direct computation is
+    more accurate since the precise angular distances between pixels are taken
+    into account, but FFT is faster, but less accurate since it assumes
+    Euclidean distances between pixels.
+
 Standard parameters
 -------------------
@@ -94,6 +94,7 @@ public:
     // Make methods private in Python by prepending an underscore
     %rename(_first_unbinned_observation) first_unbinned_observation;
     %rename(_next_unbinned_observation)  next_unbinned_observation;
+    %rename(_read_ogip_keywords)         read_ogip_keywords;
     %rename(_write_ogip_keywords)        write_ogip_keywords;
     %rename(_set_obs_statistic)          set_obs_statistic;
     %rename(_save_events_fits)           save_events_fits;
 /***************************************************************************
  *                       ctskymap - Sky mapping tool                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -30,14 +30,14 @@
 #endif
 #include <cstdio>
 #include "ctskymap.hpp"
-#include "GSkyRegions.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_PARAMETERS                          "ctskymap::get_parameter()"
 #define G_MAP_EVENTS                 "ctskymap::map_events(GCTAObservation*)"
 #define G_MAP_BACKGROUND_IRF "ctskymap::map_background_irf(GCTAObservation*)"
-#define G_MAP_BACKGROUND_RING                "ctskymap::map_background_ring("\
-                                                          "GCTAObservation*)"
+#define G_RING_BOUNDING_BOX  "ctskymap::ring_bounding_box(int&, int&, int&, "\
+                                                                "int&, int&)"
+#define G_RING_KERNEL               "ctskymap::ring_kernel(double&, double&)"
 /* __ Debug definitions __________________________________________________ */
@@ -222,37 +222,18 @@ void ctskymap::run(void)
     // Get parameters
     get_parameters();
-    // Write input observation container into logger
-    log_observations(NORMAL, m_obs, "Input observation");
-
-    // Write header into logger
-    log_header1(TERSE, gammalib::number("Map observation", m_obs.size()));
-
-    // Loop over all unbinned CTA observations in the container
-    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
-         obs = next_unbinned_observation()) {
-
-        // Map events into sky map
-        map_events(obs);
+    // Setup maps
+    setup_maps();
-        // Compute background sky map
-        map_background(obs);
-
-        // Dispose events to free memory
-        obs->dispose_events();
-
-    } // endfor: looped over observations
-
-    // If background subtraction is selected then compute significance map
-    if (m_bkgsubtract != "NONE") {
-
-        // Compute the significance map
-        map_significance();
-
-        // Subtract background map from counts map
-        m_skymap -= m_bkgmap;
+    // If no input map exists then fill counts and acceptance maps from
+    // observation container
+    if (!m_has_inmap) {
+        fill_maps();
     }
+    // Compute maps
+    compute_maps();
+
     // Optionally publish sky map
     if (m_publish) {
         publish();
@@ -266,8 +247,8 @@ void ctskymap::run(void)
 /***********************************************************************//**
  * @brief Save sky map
  *
- * Saves the sky map into a FITS file. The FITS file name is specified by the
- * @p outname parameter.
+ * Saves the sky map into a FITS file. The FITS file name is specified by
+ * the @p outname parameter.
  ***************************************************************************/
 void ctskymap::save(void)
 {
@@ -284,39 +265,31 @@ void ctskymap::save(void)
         GFits fits;
         // Write sky map into FITS file
-        GFitsHDU* hdu = m_skymap.write(fits);
-
-        // Write keywords into sky map extension
-        write_ogip_keywords(hdu);
-        write_hdu_keywords(hdu);
+        write_map(fits, m_skymap, "SKYMAP");
-        // If background subtraction is requested then write background map
-        // and significance map to FITS file
+        // If background subtraction is requested then write the background map,
+        // the significance map, the counts map and the acceptance map to the
+        // FITS file
         if (m_bkgsubtract != "NONE") {
-        
-            // Write background map into FITS file
-            hdu = m_bkgmap.write(fits);
-            // Set background map extension name
-            if (hdu != NULL) {
-                hdu->extname("BACKGROUND");
-            }
+            // Write background map into FITS file
+            write_map(fits, m_bkgmap, "BACKGROUND");
-            // Write keywords into background extension
-            write_ogip_keywords(hdu);
-            write_hdu_keywords(hdu);
             // Write significance map into FITS file
-            hdu = m_sigmap.write(fits);
+            write_map(fits, m_sigmap, "SIGNIFICANCE");
-            // Set significance map extension name
-            if (hdu != NULL) {
-                hdu->extname("SIGNIFICANCE");
-            }
+            // Write counts map into FITS file
+            write_map(fits, m_counts, "COUNTS");
+
+            // Write acceptance map into FITS file
+            write_map(fits, m_acceptance, "ACCEPTANCE");
-            // Write keywords into significance extension
-            write_ogip_keywords(hdu);
-            write_hdu_keywords(hdu);
+            // If RING background subtraction is requested then write also
+            // the exclusion map to the FITS file
+            if (m_bkgsubtract == "RING") {
+                write_map(fits, m_exclmap, "EXCLUSION");
+            }
         } // endif: background subtraction was requested
@@ -375,25 +348,36 @@ void ctskymap::publish(const std::string&amp; name)
  ***************************************************************************/
 void ctskymap::init_members(void)
 {
-    // Initialise members
-    m_skymap.clear();
-    m_bkgmap.clear();
-    m_sigmap.clear();
-    m_exclmap.clear();
-    m_alphamap.clear();
-    m_onmap.clear();
+    // Initialise User parameters
+    m_outmap.clear();
+    m_inexclusion.clear();
     m_emin        = 0.0;
     m_emax        = 0.0;
     m_bkgsubtract = "NONE";
     m_roiradius   = 0.0;
     m_inradius    = 0.0;
     m_outradius   = 0.0;
+    m_iterations  =   0;
+    m_threshold   = 5.0;
+    m_usefft      = true;
     m_publish     = false;
     m_chatter     = static_cast<GChatter>(2);
+    // Initialise members
+    m_has_inmap = false;
+    m_skymap.clear();
+    m_bkgmap.clear();
+    m_sigmap.clear();
+    m_counts.clear();
+    m_acceptance.clear();
+    m_exclmap.clear();
+
     // Initialise cache
     m_solidangle.clear();
     m_dirs.clear();
+    m_cos_roiradius = 1.0;
+    m_cos_inradius  = 1.0;
+    m_cos_outradius = 1.0;
     // Return
     return;
@@ -407,25 +391,36 @@ void ctskymap::init_members(void)
  ***************************************************************************/
 void ctskymap::copy_members(const ctskymap& app)
 {
-    // Copy attributes
-    m_skymap      = app.m_skymap;
-    m_bkgmap      = app.m_bkgmap;
-    m_sigmap      = app.m_sigmap;
-    m_exclmap     = app.m_exclmap;
-    m_alphamap    = app.m_alphamap;
-    m_onmap       = app.m_onmap;
+    // Copy User parameters
+    m_outmap      = app.m_outmap;
+    m_inexclusion = app.m_inexclusion;
     m_emin        = app.m_emin;
     m_emax        = app.m_emax;
     m_bkgsubtract = app.m_bkgsubtract;
     m_roiradius   = app.m_roiradius;
     m_inradius    = app.m_inradius;
     m_outradius   = app.m_outradius;
+    m_iterations  = app.m_iterations;
+    m_threshold   = app.m_threshold;
+    m_usefft      = app.m_usefft;
     m_publish     = app.m_publish;
     m_chatter     = app.m_chatter;
+    // Copy members
+    m_has_inmap  = app.m_has_inmap;
+    m_skymap     = app.m_skymap;
+    m_bkgmap     = app.m_bkgmap;
+    m_sigmap     = app.m_sigmap;
+    m_counts     = app.m_counts;
+    m_acceptance = app.m_acceptance;
+    m_exclmap    = app.m_exclmap;
+
     // Copy cache
-    m_solidangle = app.m_solidangle;
-    m_dirs       = app.m_dirs;
+    m_solidangle    = app.m_solidangle;
+    m_dirs          = app.m_dirs;
+    m_cos_roiradius = app.m_cos_roiradius;
+    m_cos_inradius  = app.m_cos_inradius;
+    m_cos_outradius = app.m_cos_outradius;
     // Return
     return;
@@ -453,41 +448,92 @@ void ctskymap::free_members(void)
  ***************************************************************************/
 void ctskymap::get_parameters(void)
 {
-    // Setup observations from "inobs" parameter. Do not request response
-    // information and do not accept counts cubes.
-    setup_observations(m_obs, false, true, false);
+    // So far we have no valid input map
+    m_has_inmap = false;
+
+    // If an input map was specified then extract try loading the COUNTS
+    // and ACCEPTANCE extensions from that file
+    GFilename inmap = (*this)["inmap"].filename();
+    if (inmap != "NONE") {
+        GFits fits(inmap);
+        if (fits.contains("COUNTS") && fits.contains("ACCEPTANCE")) {
+
+            // Get HDUs
+            GFitsHDU &counts     = *fits["COUNTS"];
+            GFitsHDU &acceptance = *fits["ACCEPTANCE"];
+
+            // Read counts and acceptance map
+            m_counts.read(counts);
+            m_acceptance.read(acceptance);
-    // Create sky map based on task parameters
-    m_skymap = create_map(m_obs);
+            // Extract OGIP parameters
+            read_ogip_keywords(&counts);
-    // Get further parameters
-    m_emin        = (*this)["emin"].real();
-    m_emax        = (*this)["emax"].real();
+            // Extract energy range
+            m_emin = (counts.has_card("E_MIN")) ? counts.real("E_MIN") : 0.0;
+            m_emax = (counts.has_card("E_MAX")) ? counts.real("E_MAX") : 0.0;
+
+            // Signal availability of sky map
+            m_has_inmap = true;
+        }
+        else {
+            std::string msg = "Input sky map \""+inmap.url()+"\" does not "
+                              "contain \"COUNTS\" and \"ACCEPTANCE\" "
+                              "extensions.";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
+        }
+        fits.close();
+    }
+
+    // If no valid input map exists then request observations
+    if (!m_has_inmap) {
+
+        // Setup observations from "inobs" parameter. Do not request response
+        // information and do not accept counts cubes.
+        setup_observations(m_obs, false, true, false);
+
+        // Create counts and acceptance maps from the User parameters
+        m_counts     = create_map(m_obs);
+        m_acceptance = m_counts;
+
+        // Get further parameters
+        m_emin = (*this)["emin"].real();
+        m_emax = (*this)["emax"].real();
+
+    } // endif: no valid input map existed
+
+    // Get background method
     m_bkgsubtract = (*this)["bkgsubtract"].string();
     // Get RING background parameters
     if (m_bkgsubtract == "RING") {
         // Get parameters
-        m_roiradius = (*this)["roiradius"].real();
-        m_inradius  = (*this)["inradius"].real();
-        m_outradius = (*this)["outradius"].real();
+        m_roiradius   = (*this)["roiradius"].real();
+        m_inradius    = (*this)["inradius"].real();
+        m_outradius   = (*this)["outradius"].real();
+        m_iterations  = (*this)["iterations"].integer();
+        if (m_iterations > 0) {
+            m_threshold  = (*this)["threshold"].real();
+        }
+        m_inexclusion = (*this)["inexclusion"].filename();
+        m_usefft      = (*this)["usefft"].boolean();
         // Make sure that (roiradius < inradius < outradius)
         if (m_roiradius > m_inradius) {
-            std::string msg("'roiradius' must be smaller than 'inradius'");
+            std::string msg("\"roiradius\" must be smaller than \"inradius\"");
             throw GException::invalid_value(G_GET_PARAMETERS, msg);
         }
         else if (m_inradius > m_outradius) {
-            std::string msg("'inradius' must be smaller than 'outradius'");
+            std::string msg("\"inradius\" must be smaller than \"outradius\"");
             throw GException::invalid_value(G_GET_PARAMETERS, msg);
         }
     } // endif: read RING background parameters
-    // If IRF background subtraction is requested then make sure that the
-    // CTA observations in the observation container have response information
-    if (m_bkgsubtract != "NONE") {
+    // If background subtraction is requested then make sure that the CTA
+    // observations in the observation container have response information
+    if (m_bkgsubtract != "NONE" && !m_has_inmap) {
         set_response(m_obs);
     }
@@ -497,39 +543,209 @@ void ctskymap::get_parameters(void)
     // Read ahead parameters
     if (read_ahead()) {
-        m_outmap    = (*this)["outmap"].filename();
+        m_outmap = (*this)["outmap"].filename();
     }
+    // Write parameters into logger
+    log_parameters(TERSE);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Setup maps
+ *
+ * Setup maps for sky map generation.
+ ***************************************************************************/
+void ctskymap::setup_maps(void)
+{
+    // Clear cache
+    m_solidangle.clear();
+    m_dirs.clear();
+
     // Create background map and significance map if background subtraction
     // is requested
     if (m_bkgsubtract != "NONE") {
-        // Create backgrond and significance maps
-        m_bkgmap = create_map(m_obs);
-        m_sigmap = create_map(m_obs);
-        
         // Setup the exclusions map
-        map_exclusions((*this)["inexclusion"].filename());
+        setup_exclusion_map(m_inexclusion);
         // Cache the pixel solid angles and sky directions
-        m_solidangle.reserve(m_bkgmap.npix());
-        m_dirs.reserve(m_bkgmap.npix());
-        for (int i = 0; i < m_bkgmap.npix(); ++i) {
-            m_solidangle.push_back(m_bkgmap.solidangle(i));
-            m_dirs.push_back(m_bkgmap.inx2dir(i));
+        m_solidangle.reserve(m_counts.npix());
+        m_dirs.reserve(m_counts.npix());
+        for (int i = 0; i < m_counts.npix(); ++i) {
+            m_solidangle.push_back(m_counts.solidangle(i));
+            m_dirs.push_back(m_counts.inx2dir(i));
         }
-        // If doing a ring background subtraction, generate an alpha map and
-        // make sure to cache the counts for each observation
-        if (m_bkgsubtract == "RING") {
-            m_alphamap = create_map(m_obs);
-            m_onmap    = create_map(m_obs);
+    } // endif: background subtraction selected
+
+    // Compute cosine of ring radii for RING background method
+    if (m_bkgsubtract == "RING") {
+        m_cos_roiradius = std::cos(m_roiradius * gammalib::deg2rad);
+        m_cos_inradius  = std::cos(m_inradius  * gammalib::deg2rad);
+        m_cos_outradius = std::cos(m_outradius * gammalib::deg2rad);
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Generates map of pixel exclusions
+ *
+ * @param[in] filename Exclusion file name.
+ *
+ * Generates a sky map of the pixels that are to be excluded from the
+ * background estimation. Pixels with values different from 0 will be
+ * excluded.
+ ***************************************************************************/
+void ctskymap::setup_exclusion_map(const GFilename& filename)
+{
+    // Create exlusion map
+    m_exclmap = m_counts;
+
+    // Set all pixels to 0 (no pixel excluded)
+    m_exclmap = 0.0;
+
+    // Make sure the exclusions filename is valid
+    if (is_valid_filename(filename)) {
+
+        // Fill the exclusions based on the regions supplied
+        if (filename.is_fits()) {
+            setup_exclusion_map_fits(filename);
+        }
+        else {
+            setup_exclusion_map_region(filename);
         }
-    } // endif: background subtraction selected
+    } // endif: filename was valid
-    // Write parameters into logger
-    log_parameters(TERSE);
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Fills exclusions map from FITS image
+ *
+ * @param[in] filename FITS image file name.
+ *
+ * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in
+ * the exclusion sky map FITS file.
+ ***************************************************************************/
+void ctskymap::setup_exclusion_map_fits(const GFilename& filename)
+{
+    // Load the fits image
+    GSkyMap inmap(filename);
+
+    // Loop through the individual pixels in the exclusion map
+    for (int i = 0; i < m_exclmap.npix(); ++i) {
+
+        // Get the pixel direction
+        GSkyDir dir = m_exclmap.pix2dir(i);
+
+        // Check this sky position in the fits map
+        if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
+
+            // Set the pixel to 1
+            m_exclmap(i) = 1.0;
+
+        } // endif: pixel,region overlap check
+
+    } // endfor: looped over exclusion map pixels
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Fills exclusions map from DS9 region file
+ *
+ * @param[in] filename DS9 region file name.
+ *
+ * Sets all exclusion map pixels to 1 that are contained in any of the DS9
+ * regions.
+ ***************************************************************************/
+void ctskymap::setup_exclusion_map_region(const GFilename& filename)
+{
+    // Load the exclusion regions
+    GSkyRegions regions(filename);
+
+    // Loop through the individual pixels in the exclusion map
+    for (int i = 0; i < m_exclmap.npix(); ++i) {
+
+        // Get the pixel position
+        GSkyDir dir = m_exclmap.pix2dir(i);
+
+        // If pixel position overlaps with the regions
+        if (regions.contains(dir)) {
+            m_exclmap(i) = 1.0;
+        }
+
+    } // endfor: looped over exclusion map pixels
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Fill maps from observation container
+ *
+ * Fill counts and optionally acceptance maps from data in the observation
+ * container.
+ ***************************************************************************/
+void ctskymap::fill_maps(void)
+{
+    // Write input observation container into logger
+    log_observations(NORMAL, m_obs, "Input observation");
+
+    // Write header into logger
+    log_header1(TERSE, gammalib::number("Find unbinned observation",
+                m_obs.size()));
+
+    // Find all unbinned CTA observations in m_obs
+    std::vector<GCTAObservation*> obs_list(0);
+    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+            obs = next_unbinned_observation()) {
+
+        // Push observation into list
+        obs_list.push_back(obs);
+
+        // Write message
+        std::string msg = " Including unbinned "+obs->instrument()+
+                          " observation";
+        log_string(NORMAL, msg);
+    }
+
+    // Write header into logger
+    log_header1(TERSE, gammalib::number("Fill map from observation",
+                m_obs.size()));
+
+    // Loop over all unbinned CTA observations in the container
+    #pragma omp parallel for
+    for (int i = 0; i < obs_list.size(); ++i) {
+
+        // Get pointer to observation
+        GCTAObservation* obs = obs_list[i];
+
+        // Fill events into counts map
+        fill_maps_counts(obs);
+
+        // Optionally compute acceptance sky map
+        if (m_bkgsubtract != "NONE") {
+            fill_maps_acceptance(obs);
+        }
+
+        // Dispose events to free memory
+        obs->dispose_events();
+
+    } // endfor: looped over observations
     // Return
     return;
@@ -537,16 +753,17 @@ void ctskymap::get_parameters(void)
 /***********************************************************************//**
- * @brief Map events into a sky map
+ * @brief Fill events into counts map
  *
  * @param[in] obs CTA observation.
  *
  * @exception GException::no_list
  *            No event list found in observation.
  *
- * This method maps the events found in a CTA events list into a sky map.
+ * Fills the events found in a CTA events list into a sky map. The method
+ * adds the events to the m_counts member.
  ***************************************************************************/
-void ctskymap::map_events(GCTAObservation* obs)
+void ctskymap::fill_maps_counts(GCTAObservation* obs)
 {
     // Get non-const pointer on a CTA event list
     GCTAEventList* events = dynamic_cast<GCTAEventList*>
@@ -589,11 +806,11 @@ void ctskymap::map_events(GCTAObservation* obs)
         // Determine sky pixel
         GCTAInstDir* inst  = (GCTAInstDir*)&(event->dir());
         GSkyDir      dir   = inst->dir();
-        GSkyPixel    pixel = m_skymap.dir2pix(dir);
+        GSkyPixel    pixel = m_counts.dir2pix(dir);
         // Skip if pixel is out of range
-        if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) ||
-            pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) {
+        if (pixel.x() < -0.5 || pixel.x() > (m_counts.nx() - 0.5) ||
+            pixel.y() < -0.5 || pixel.y() > (m_counts.ny() - 0.5)) {
             num_outside_map++;
             continue;
         }
@@ -605,49 +822,25 @@ void ctskymap::map_events(GCTAObservation* obs)
         }
         // Fill event in skymap
-        m_skymap(pixel, 0) += 1.0;
+        #pragma omp critical(ctskymap_map_events)
+        m_counts(pixel) += 1.0;
         num_in_map++;
-        
+
     } // endfor: looped over all events
     // Log binning results
-    log_value(NORMAL, "Events in list", obs->events()->size());
-    log_value(NORMAL, "Events in map", num_in_map);
-    log_value(NORMAL, "Events outside RoI", num_outside_roi);
-    log_value(NORMAL, "Events outside map area", num_outside_map);
-    log_value(NORMAL, "Events outside energies", num_outside_erange);
-
-    // Write sky map into header
-    log_header1(EXPLICIT, "Sky map");
-    log_string(EXPLICIT, m_skymap.print(m_chatter));
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Estimates the background in sky map for observation
- *
- * @param[in] obs CTA observation.
- *
- * Estimates the background in the sky map for a given observation and adds
- * this estimate to the background sky map. The background estimation method
- * is specified by the "bkgsubtract" parameter which can take the following
- * values:
- *
- *     NONE - No background estimation
- *     IRF  - Background estimation based on IRF template
- *     RING - Ring background estimation
- ***************************************************************************/
-void ctskymap::map_background(GCTAObservation* obs)
-{
-    // Dispatch to appropriate background estimation method
-    if (m_bkgsubtract == "IRF") {
-        map_background_irf(obs);
-    }
-    else if (m_bkgsubtract == "RING") {
-        map_background_ring(obs);
+    #pragma omp critical(ctskymap_map_events)
+    {
+        log_header3(TERSE, get_obs_header(obs));
+        log_value(NORMAL, "Events in list", obs->events()->size());
+        log_value(NORMAL, "Events in map", num_in_map);
+        log_value(NORMAL, "Events outside RoI", num_outside_roi);
+        log_value(NORMAL, "Events outside map area", num_outside_map);
+        log_value(NORMAL, "Events outside energies", num_outside_erange);
+
+        // Write sky map into header
+        log_header1(EXPLICIT, "Sky map");
+        log_string(EXPLICIT, m_counts.print(m_chatter));
     }
     // Return
@@ -656,7 +849,7 @@ void ctskymap::map_background(GCTAObservation* obs)
 /***********************************************************************//**
- * @brief Estimates the background in sky map based on IRF template
+ * @brief Compute background acceptance sky map based on IRF template
  *
  * @param[in] obs CTA observation.
  *
@@ -664,14 +857,14 @@ void ctskymap::map_background(GCTAObservation* obs)
  *            No response information available for observation.
  *            No background template available in instrument response function.
  *
- * Estimates the background in the sky map using the IRF template for a given
- * observation and adds this estimate to the background sky map. The IRF
- * template is integrated numerically in energy.
+ * Computes a background acceptance sky map using the IRF template for a
+ * given observation and adds the estimate to the acceptance sky map.
  ***************************************************************************/
-void ctskymap::map_background_irf(GCTAObservation* obs)
+void ctskymap::fill_maps_acceptance(GCTAObservation* obs)
 {
     // Get IRF response
-    const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
+    const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
+                                 (obs->response());
     // Throw an exception if observation has no instrument response function
     if (rsp == NULL) {
@@ -694,60 +887,58 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
         throw GException::invalid_value(G_MAP_BACKGROUND_IRF, msg);
     }
-    // Compute natural logarithm of energy range in MeV
-    double lnEmin = std::log(m_emin * 1.0e6);
-    double lnEmax = std::log(m_emax * 1.0e6);
+    // Set minimum and maximum energy as GEnergy instances
+    GEnergy emin(m_emin, "TeV");
+    GEnergy emax(m_emax, "TeV");
-    // Extract region of interest from observation
-    GCTARoi roi = obs->roi();
+    // Extract region of interest center and radius from observation and
+    // compute the cosine of the RoI radius. This allows to mimimize as
+    // much as possible the trigonometric computations. Note that the cosine
+    // of an angle is maximal for an angle of zero. This explains later why
+    // we chose ">=" to test whether an angle is smaller than a given radius.
+    GSkyDir roi_centre     = obs->roi().centre().dir();
+    double  cos_roi_radius = std::cos(obs->roi().radius() * gammalib::deg2rad);
+
+    // Extract exposure
+    double exposure = obs->livetime();
     // Initialise statistics
-    int    calls = 0;
     double total = 0.0;
-    // Loop over all background map pixels
-    for (int i = 0; i < m_bkgmap.npix(); ++i) {
+    // Loop over all acceptance map pixels
+    for (int i = 0; i < m_acceptance.npix(); ++i) {
         // Get sky direction of pixel
-        GSkyDir skydir = m_bkgmap.inx2dir(i);
-
-        // Convert sky direction in instrument direction
-        GCTAInstDir instdir = obs->pointing().instdir(skydir);
+        GSkyDir& skydir = m_dirs[i];
-        // If RoI is valid and instrument direction is not within RoI then
-        // skip pixel
-        if (roi.is_valid() && !roi.contains(instdir)) {
+        // Skip pixel if it is not within the RoI
+        if (roi_centre.cos_dist(skydir) < cos_roi_radius) {
             continue;
         }
-        // Setup integration function
-        ctskymap::irf_kern integrand(bkg, &instdir);
-        GIntegral          integral(&integrand);
-
-        // Set precision (has been carefully adjusted using a test simulation
-        // over the energy range 20 GeV - 120 TeV)
-        integral.eps(1.0e-6);
+        // Convert sky direction to instrument direction
+        GCTAInstDir instdir = obs->pointing().instdir(skydir);
-        // Do Romberg integration
-        double value = integral.romberg(lnEmin, lnEmax);
+        // Compute background value
+        double value = bkg->rate_ebin(instdir, emin, emax);
-        // Update number of background function calls
-        calls += integral.calls();
-        
         // Multiply background rate with livetime and solid angle
-        value *= obs->livetime() * m_solidangle[i];
+        value *= exposure * m_solidangle[i];
-        // Add number of background events to background map
-        m_bkgmap(i,0) += value;
+        // Add number of background events to acceptance map
+        #pragma omp critical(ctskymap_map_acceptance)
+        m_acceptance(i) += value;
         // Update total number of background events
         total += value;
-    } // endfor: looped over background map pixels
+    } // endfor: looped over acceptance map pixels
-    // Log background subtraction results
-    log_value(NORMAL, "Events in background", int(total+0.5));
-    log_value(NORMAL, "Background evaluations", calls);
+    // Log acceptance results
+    #pragma omp critical(ctskymap_map_events)
+    {
+        log_value(NORMAL, "Events in background", int(total+0.5));
+    }
     // Return
     return;
@@ -755,97 +946,99 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
 /***********************************************************************//**
- * @brief Estimates the background in sky map based on the ring background
- *        method
+ * @brief Compute sky map, background map and significance map
  *
- * @param[in] obs CTA observation.
+ * Computes the sky map, background map and significance map. The following
+ * background subtraction methods are supported:
  *
- * @exception GException::invalid_value
- *            No response information available for observation.
- *            No background template available in instrument response function.
+ * NONE
+ * ====
+ *
+ * The sky map is simply the binned counts map, no background and
+ * significance maps are computed.
+ *
+ * IRF
+ * ===
+ * 
+ * The background map is the acceptance map, and this background map is
+ * subtracted from the counts map and assigned to the sky map. The
+ * significance \f$\sigma_i\f$ is computed assuming Poisson statistic in
+ * the Gaussian limit, using
+ *
+ * \f[\sigma_i = \frac{N_i - B_i}{\sqrt{N_i}}\f]
  *
- * Estimates the background in the sky map by summing the events within a
- * ring centered at a given pixel's position. The pixels in the ring are
- * weighted also by the background IRF value.
+ * where
+ * \f$N_i\f$ is the number of observed counts and
+ * \f$B_i\f$ is the estimated number of background counts.
+ *
+ * RING
+ * ====
+ *
+ * For the @c RING method the Li & Ma significance is computed for each
+ * pixel, using the @c roiradius, @c inradius and @c outradius parameters
+ * to specify the radius of the On region and the Off background ring,
+ * respectively. Depending on the @c usefft parameter, either an FFT is used
+ * to compute the number of events and the acceptance of the On and Off
+ * regions, or a direct computation is performed. While the former is faster
+ * it is less accurate since it assumes Euclidean distances in the sky map.
+ * The latter is slower but compute exact distance between pixels of the sky
+ * map. Use FFT if the sky map is close to a cartesian grid, and use the
+ * direct method if the sky map shows important distortions.
  ***************************************************************************/
-void ctskymap::map_background_ring(GCTAObservation* obs)
+void ctskymap::compute_maps(void)
 {
-    // Get IRF response (to scale background counts)
-    const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
+    // Write header into logger
+    log_header1(TERSE, "Compute maps");
-    // Throw an exception if observation has no instrument response function
-    if (rsp == NULL) {
-        std::string msg = "No response information available for "+
-                          get_obs_header(obs)+" to compute IRF background. "
-                          "Please specify response information or use "
-                          "another background subtraction method.";
-        throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg);
+    // If no background subtraction is requested then simply assign the
+    // counts map to the skymap
+    if (m_bkgsubtract == "NONE") {
+        m_skymap = m_counts;
     }
-    // Get IRF background template
-    const GCTABackground* bkg = rsp->background();
-
-    // Throw an exception if observation has no IRF background template
-    if (bkg == NULL) {
-        std::string msg = "No IRF background template found in instrument "
-                          "response function for "+
-                          get_obs_header(obs)+". Please specify an instrument "
-                          "response function containing a background template.";
-        throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg);
+    // else if IRF background is selected then use the acceptance as background
+    else if (m_bkgsubtract == "IRF") {
+        m_skymap = m_counts - m_acceptance;
+        m_bkgmap = m_acceptance;
+        m_sigmap = (m_counts - m_acceptance) / sqrt(m_counts);
     }
-    // Compute natural logarithm of energy range in MeV
-    double lnEmin = std::log(m_emin * 1.0e6);
-    double lnEmax = std::log(m_emax * 1.0e6);
-
-    // Extract region of interest from observation
-    GCTARoi roi = obs->roi();
-
-    // Initialise statistics
-    int    calls    = 0;
-    double exposure = obs->livetime();
-
-    // Loop over all map pixels
-    for (int i = 0; i < m_bkgmap.npix(); ++i) {
-
-        // Get sky direction of pixel
-        GSkyDir& skydir = m_dirs[i];
+    // else if RING background is selected then compute the On
+    else if (m_bkgsubtract == "RING") {
-        // Convert sky direction in instrument direction
-        GCTAInstDir instdir = obs->pointing().instdir(skydir);
-    
-        // If RoI is valid and instrument direction is not within RoI then
-        // skip pixel
-        if (roi.is_valid() && !roi.contains(instdir)) {
-            continue;
+        // Compute initial RING background
+        if (m_usefft) {
+            compute_maps_ring_fft();
+        }
+        else {
+            compute_maps_ring_direct();
         }
-        // Setup integration function
-        ctskymap::irf_kern integrand(bkg, &instdir);
-        GIntegral          integral(&integrand);
-
-        // Set precision
-        integral.eps(1.0e-6);
-
-        // Do Romberg integration to get the sensitivity (note: this assumes
-        // the background IRF is a good approximation for the radial
-        // sensitivity)
-        m_alphamap(i,0) += integral.romberg(lnEmin, lnEmax) * 
-                           m_solidangle[i] * exposure;
-
-        // Store the counts in this bin
-        m_onmap(i,0) += m_skymap(i,0);
+        // Store copy of exclusion regions
+        GSkyMap exclmap = m_exclmap;
-        // Update number of background function calls
-        calls += integral.calls();
+        // Iterative computation of exclusion regions
+        for (int iter = 0; iter < m_iterations; ++iter) {
-    } // endfor: Loop for caching bkg IRF sensitivity
+            // Set exclusion region from significance map
+            m_exclmap = exclmap;
+            for (int i = 0; i < m_exclmap.npix(); ++i) {
+                if (m_sigmap(i) > m_threshold) {
+                    m_exclmap(i) = 1.0;
+                }
+            }
+            
+            // Re-compute RING background
+            if (m_usefft) {
+                compute_maps_ring_fft();
+            }
+            else {
+                compute_maps_ring_direct();
+            }
-    // Log background subtraction results
-    log_value(NORMAL, "Alpha map evaluations", calls);
+        } // endfor: looped over iterations
-    // Zero out the counts map to prepare it for the next observation
-    m_skymap = 0.0;
+    }
     // Return
     return;
@@ -853,24 +1046,85 @@ void ctskymap::map_background_ring(GCTAObservation* obs)
 /***********************************************************************//**
- * @brief Write keywords in FITS HDU
- *
- * @param[in,out] hdu Pointer to FITS HDU.
+ * @brief Compute the maps for RING background using a FFT
  *
- * Writes keywords in FITS HDU.
+ * Computes the Li & Ma significance for each sky map pixel and replaces
+ * the sky and background maps by the On- and Off-count maps. The computation
+ * is done using a FFT.
  ***************************************************************************/
-void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
+void ctskymap::compute_maps_ring_fft(void)
 {
-    // Continue only if pointer is valid
-    if (hdu != NULL) {
+    // Log message about what is being done
+    log_header3(NORMAL, "Computing ring background map (FFT method)");
-        // Set keywords
-        hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
-        hdu->card("E_MIN",  m_emin, "[TeV] Lower energy boundary");
-        hdu->card("E_MAX",  m_emax, "[TeV] Upper energy boundary");
-        hdu->card("EUNIT",  "TeV",  "Units for E_MIN and E_MAX");
+    // Initialise maps
+    m_skymap = m_counts;
+    m_bkgmap = m_counts;
+    m_sigmap = m_counts;
+    m_skymap = 0.0;
+    m_bkgmap = 0.0;
+    m_sigmap = 0.0;
-    } // endif: pointer was valid
+    // Initialise statistics
+    int num_bad_alpha = 0;
+
+    // Set-up sky and acceptance maps with excluded pixels according to the
+    // exclusion map
+    GSkyMap excl_counts     = m_counts;
+    GSkyMap excl_acceptance = m_acceptance;
+    for (int i = 0; i < m_counts.npix(); ++i) {
+        if (m_exclmap(i) != 0.0) {
+            excl_counts(i)     = 0.0;
+            excl_acceptance(i) = 0.0;
+        }
+    }
+
+    // Convolve maps for On-region and Off-region
+    GSkyMap on_counts  = ring_convolve(m_counts,        0.0,        m_roiradius);
+    GSkyMap on_alpha   = ring_convolve(m_acceptance,    0.0,        m_roiradius);
+    GSkyMap off_counts = ring_convolve(excl_counts,     m_inradius, m_outradius);
+    GSkyMap off_alpha  = ring_convolve(excl_acceptance, m_inradius, m_outradius);
+
+    // Compute maps by loop over sky map pixels
+    for (int i = 0; i < m_sigmap.npix(); ++i) {
+
+        // Get On-counts, Off-counts and alpha value for this bin
+        double n_on  = on_counts(i);
+        double n_off = off_counts(i);
+        double alpha;
+        if (off_alpha(i) == 0.0) {
+            alpha = 0.0;
+            n_on  = 0.0;
+            n_off = 0.0;
+        }
+        else {
+            alpha = on_alpha(i) / off_alpha(i);
+        }
+
+        // Store the On and alpha-weighted Off-counts
+        m_bkgmap(i) = alpha * n_off;
+        m_skymap(i) = n_on - m_bkgmap(i);
+
+        // If alpha is zero then increment the bad alpha counter
+        if (alpha == 0.0) {
+            num_bad_alpha++;
+        }
+
+        // ... otherwise compute and store the results
+        else {
+            m_sigmap(i) = sigma_li_ma(n_on, n_off, alpha);
+        }
+
+    } // endfor: looped over all pixels
+
+    // Log the number of bad-alpha bins
+    log_value(NORMAL, "Total number of pixels", m_sigmap.npix());
+    log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha);
+
+    // Now take the square root. Since some bins can be negative, first
+    // take the absolute value of the map, square root that, then multiply
+    // each bin by its sign to preserve the +/- significance.
+    m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
     // Return
     return;
@@ -878,89 +1132,85 @@ void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
 /***********************************************************************//**
- * @brief Computes the bin-by-bin significance if background subtraction
- *        requested
+ * @brief Compute the pixel significance for RING background
  *
- * Method computes the bin-by bin significance if a background subtraction 
- * method is specified. If method is "IRF", Poisson statistics (in the
- * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance
- * is computed for each bin (as well as a map of the computed alpha values).
+ * Computes the Li & Ma significance for each sky map pixel and replaces
+ * the sky and background maps by the On- and Off-count maps. The computation
+ * is done by computing the exact distances between pixels in the sky map.
  ***************************************************************************/
-void ctskymap::map_significance(void)
+void ctskymap::compute_maps_ring_direct(void)
 {
-    // Compute significance from "RING" method (Li & Ma eq. 17)
-    if (m_bkgsubtract == "RING") {
-
-        // Log message about what is being done
-        log_header1(NORMAL, "Computing Ring background map");
-        log_value(NORMAL, "Total pixels to process", m_onmap.npix());
-
-        // Store the number of bins with an inappropriate alpha parameter
-        int num_bad_alpha = 0;
-
-        // Loop through each bin in the on-counts map
-        for (int i = 0; i < m_onmap.npix(); ++i) {
-
-            // Initialise the on/off counts and alpha for this bin
-            double n_on  = 0.0;
-            double n_off = 0.0;
-            double alpha = 0.0;
+    // Log message about what is being done
+    log_header3(NORMAL, "Computing ring background map (direct method)");
+    log_value(NORMAL, "Total pixels to process", m_counts.npix());
+
+    // Initialise maps
+    m_skymap = m_counts;
+    m_bkgmap = m_counts;
+    m_sigmap = m_counts;
+    m_skymap = 0.0;
+    m_bkgmap = 0.0;
+    m_sigmap = 0.0;
-            // Since this can take a long time, keep the user updated on
-            // the progress when another 10% of pixels is processed
-            if (i % (m_onmap.npix()/10) == 0) {
-                log_value(NORMAL, "Pixels remaining", m_onmap.npix()-i);
-            }
+    // Initialise statistics
+    int num_bad_alpha = 0;
+
+    // Set progress logging frequency. By default the user is updated on the
+    // progress when another 10% of pixels is processed, however, if there
+    // are many pixels the user is updated each 100000 pixels.
+    int n_logpix = m_counts.npix() / 10;
+    if (n_logpix > 100000) {
+        n_logpix = 100000;
+    }
-            // Get bin coordinates
-            GSkyDir& skydir = m_dirs[i];
+    // Compute maps by loop over sky map pixels
+    #pragma omp parallel for
+    for (int i = 0; i < m_sigmap.npix(); ++i) {
-            // Compute the alpha and counts for this bin
-            compute_ring_values(m_onmap, m_alphamap, skydir,
-                                n_on, n_off, alpha);
+        // Initialise the on/off counts and alpha for this bin
+        double n_on  = 0.0;
+        double n_off = 0.0;
+        double alpha = 0.0;
-            // If alpha is zero then increment the bad alpha counter
-            if (alpha == 0.0) {
-                num_bad_alpha++;
-            }
+        // Update user about progress (only if there are more than 10000 pixels)
+        if (i % n_logpix == 0 && n_logpix > 10000) {
+            #pragma omp critical(ctskymap_map_significance_ring)
+            log_value(NORMAL, "Pixels remaining", m_sigmap.npix()-i);
+        }
-            // ... otherwise store the results
-            else {
+        // Compute the alpha and counts for this bin
+        compute_ring_values(i, m_counts, m_acceptance, n_on, n_off, alpha);
-                // Store the on-counts and alpha-weighted off-counts
-                m_skymap(i,0) = n_on;
-                m_bkgmap(i,0) = alpha * n_off;
+        // Store the On and alpha-weighted Off-counts
+        m_bkgmap(i) = alpha * n_off;
+        m_skymap(i) = n_on - m_bkgmap(i);
-                // Compute and store significance (Li & Ma eq. 17)
-                if (n_on == 0.0) {
-                    m_sigmap(i,0) = -2.0 * n_off * std::log(1.0+alpha);
-                }
-                else {
-                    m_sigmap(i,0) = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
-                                    (n_on  * std::log((1.0+alpha) *
-                                     n_on / (alpha * (n_on+n_off))) +
-                                     n_off * std::log((1.0+alpha) *
-                                     n_off / (n_on+n_off)));
-                }
+        // If alpha is zero then increment the bad alpha counter
+        if (alpha == 0.0) {
+            num_bad_alpha++;
+            #pragma omp critical(ctskymap_compute_maps_ring_direct)
+            m_sigmap(i) = 0.0;
+        }
-            } // endelse: alpha was non-zero
+        // ... otherwise store the results
+        else {
-        } // endfor: looped over all pixels
+            // Compute and store significance (Li & Ma eq. 17)
+            double sigma = sigma_li_ma(n_on, n_off, alpha);
+            #pragma omp critical(ctskymap_compute_maps_ring_direct)
+            m_sigmap(i) = sigma;
-        // Log the number of bad-alpha bins
-        log_value(NORMAL, "Bins with alpha=0", num_bad_alpha);
+        } // endelse: alpha was non-zero
-        // Now take the square root. Since some bins can be negative, first
-        // take the absolute value of the map, square root that, then multiply
-        // each bin by its sign to preserve the +/- significance.
-        m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
+    } // endfor: looped over all pixels
-    } // endif: ring background selected
+    // Log the number of bad-alpha bins
+    log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha);
-    // Compute significance using Poisson statistics in the Gaussian limit
-    else {
-        m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap);
-    }
+    // Now take the square root. Since some bins can be negative, first
+    // take the absolute value of the map, square root that, then multiply
+    // each bin by its sign to preserve the +/- significance.
+    m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
     // Return
     return;
@@ -970,66 +1220,107 @@ void ctskymap::map_significance(void)
 /***********************************************************************//**
  * @brief Computes Non, Noff and alpha for a counts map and sensitivity map
  *
- * @param[in]  counts           Counts map
- * @param[in]  sensitivity      Sensitivity map
- * @param[in]  position         Position of the ring & roi centers
- * @param[out] non              Returned estimate of ON counts
- * @param[out] noff             Returned estimate of OFF counts
- * @param[out] alpha            Returned estimate of alpha
+ * @param[in]  ipixel      Sky map pixel to consider.
+ * @param[in]  counts      Counts map.
+ * @param[in]  background  Background map.
+ * @param[out] non         Returned estimate of On-counts.
+ * @param[out] noff        Returned estimate of Off-counts.
+ * @param[out] alpha       Returned estimate of alpha.
  * 
- * This method computes the Non counts Noff values at a given position. It 
- * also computes the alpha parameter from the passed sensitivity map.
+ * Computes the On- and Off-counts values at a given position from in counts
+ * map and the alpha parameter from the background map.
+ *
+ * To speed-up the computations the method considers only the pixels within
+ * a bounding box that comprises the outer background ring radius. The method
+ * considers wrapping around of pixel indices in the Right Ascension /
+ * Galactic longitude direction.
+ *
+ * If the alpha of the ring region is zero the values of @p non, @p noff,
+ * and @p alpha will be set to zero.
  ***************************************************************************/
-void ctskymap::compute_ring_values(const GSkyMap& counts, 
-                                   const GSkyMap& sensitivity,
-                                   const GSkyDir& position,
+void ctskymap::compute_ring_values(const int&     ipixel,
+                                   const GSkyMap& counts,
+                                   const GSkyMap& background,
                                    double&        non,
                                    double&        noff,
                                    double&        alpha)
 {
-    // Reset Non and Noff
-    non  = 0.0;
-    noff = 0.0;
+    // Initialise return values (non, noff and alpha)
+    non   = 0.0;
+    noff  = 0.0;
+    alpha = 0.0;
     // Initialise On and Off alpha
     double alpha_on  = 0.0;
     double alpha_off = 0.0;
-    // Define the regions necessary to do the pixel checks
-    GSkyRegionCircle roi_reg(position, m_roiradius);
-    GSkyRegionCircle inner_reg(position, m_inradius);
-    GSkyRegionCircle outer_reg(position, m_outradius);
+    // Get sky direction of pixel
+    GSkyDir& position = m_dirs[ipixel];
-    // Loop over every pixel in the observation to compute Non, Noff
-    for (int j = 0; j < counts.npix(); ++j) {
+    // Get pixel bounding box
+    int ix_start;
+    int ix_stop;
+    int iy_start;
+    int iy_stop;
+    ring_bounding_box(ipixel, ix_start, ix_stop, iy_start, iy_stop);
-        // Get the index and sky direction of this pixel
-        GSkyDir& skydir = m_dirs[j];
+    // Get number of x pixels in sky map
+    int nx = counts.nx();
-        // Check if pixel is inside the background region
-        if ((m_exclmap(j,0) == 0.0) && outer_reg.contains(skydir) &&
-            !inner_reg.contains(skydir)) {
+    // Loop over x pixels of bounding box
+    for (int ix_nowrap = ix_start; ix_nowrap < ix_stop; ++ix_nowrap) {
-            // Update n_off
-            noff += counts(j,0);
+        // Compute x pixel index by wrapping the pixel index into the
+        // interval [0,nx[
+        int ix = ix_nowrap;
+        if (ix < 0) {
+            ix += nx;
+        }
+        else if (ix >= nx) {
+            ix -= nx;
+        }
-            // Update alpha_off
-            alpha_off += sensitivity(j,0);
+        // Initialise pixel index
+        int i = ix + iy_start * nx;
-        }
+        // Loop over y pixels of bounding box
+        for (int iy = iy_start; iy < iy_stop; ++iy, i += nx) {
+
+            // Get the index and sky direction of this pixel
+            GSkyDir& skydir = m_dirs[i];
+
+            // Only consider pixels within the outer radius of the background
+            // region
+            if (position.cos_dist(skydir) >= m_cos_outradius) {
+
+                // Check if pixel is inside the background region
+                if ((m_exclmap(i) == 0.0) &&
+                    (position.cos_dist(skydir) < m_cos_inradius)) {
+
+                    // Update n_off
+                    noff += counts(i);
+
+                    // Update alpha_off
+                    alpha_off += background(i);
+
+                }
+
+                // ... otherwise check if pixel is inside source region
+                else if (position.cos_dist(skydir) >= m_cos_roiradius) {
+
+                    // Update n_on
+                    non += counts(i);
-        // ... otherwise check if pixel is inside source region
-        else if (roi_reg.contains(skydir)) {
+                    // Update alpha_on
+                    alpha_on += background(i);
-            // Update n_on for significance computation
-            non += counts(j,0);
+                } // endif: source and background region check
-            // Update alpha_on
-            alpha_on += sensitivity(j,0);
+            } // endif: pixel is within outer radius of background region
-        } // endif: source and background region check
+        } // endfor: looped over y pixels
-    } // endfor: looped over pixels
+    } // endfor: looped over x pixels
     // Compute alpha. If the off region does not have any sensitivity then
     // set Non = Noff = 0
@@ -1048,34 +1339,100 @@ void ctskymap::compute_ring_values(const GSkyMap&amp; counts,
 /***********************************************************************//**
- * @brief Generates map of pixel exclusions
+ * @brief Computes bounding box for RING background computation
  *
- * @param[in] filename Exclusion file name.
+ * @param[in]  ipixel Sky map pixel to consider.
+ * @param[out] ix1    Index of first pixel in x.
+ * @param[out] ix2    Index after last pixel in x.
+ * @param[out] iy1    Index of first pixel in y.
+ * @param[out] iy2    Index after last pixel in y.
  *
- * Generates a sky map of the pixels that are to be excluded from the
- * background estimation. Pixels with values different from 0 will be
- * excluded.
+ * Computes the bounding box that contained the background ring for a
+ * specific pixel for the RING background method.
+ *
+ * The method determines the local pixel scale for the requested pixel and
+ * draws a bounding box with 1.5 times the outer background ring radius
+ * around the pixel. In the x direction the pixel indices are unconstrained,
+ * and the client has to assure that the pixel value is comprised within
+ * the validity range (this is required to handle the longitude wrap around).
+ * In the y direction the pixel indices are constrained to [0,ny], where
+ * ny is the number of y pixels in the sky map.
  ***************************************************************************/
-void ctskymap::map_exclusions(const GFilename& filename)
+void ctskymap::ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
+                                                    int& iy1, int& iy2) const
 {
-    // Create exlusion map
-    m_exclmap = create_map(m_obs);
+    // Get number of pixels in x and y direction
+    int nx = m_counts.nx();
+    int ny = m_counts.ny();
+
+    // Get x and y index of pixel to consider
+    int ix0 = ipixel % nx;
+    int iy0 = ipixel / nx;
+
+    // Get pointer on WCS projection
+    const GWcs* wcs = dynamic_cast<const GWcs*>(m_counts.projection());
+    if (wcs == NULL) {
+        std::string msg = "Sky map is not a WCS projection. Method is only "
+                          "valid for WCS projections.";
+        throw GException::invalid_value(G_RING_BOUNDING_BOX, msg);
+    }
-    // Set all pixels to 0 (no pixel excluded)
-    m_exclmap = 0.0;
+    // Get X and Y step size
+    double dx_binsz = wcs->cdelt(0);
+    double dy_binsz = wcs->cdelt(1);
-    // Make sure the exclusions filename is valid
-    if (is_valid_filename(filename)) {
+    // Compute pixel increment in x-direction
+    double dx = 0.0;
+    int    ix = (ix0 > 0) ? ix0 - 1 : ix0 + 1;
+    if (ix < nx) {
+        dx = m_dirs[ipixel].dist_deg(m_dirs[ix+iy0*nx]);
+    }
+    if (dx < dx_binsz) {
+        dx = dx_binsz;
+    }
-        // Fill the exclusions based on the regions supplied
-        if (filename.is_fits()) {
-            map_exclusions_fits(filename);
+    // Compute pixel increment in y-direction
+    double dy = 0.0;
+    int    iy = (iy0 > 0) ? iy0 - 1 : iy0 + 1;
+    if (iy < ny) {
+        dx = m_dirs[ipixel].dist_deg(m_dirs[ix0+iy*nx]);
+    }
+    if (dy < dy_binsz) {
+        dy = dy_binsz;
+    }
+
+    // Compute bounding box half size in x and y. The outer radius is
+    // multiplied by 1.5 to have some margin in case of map distortions
+    int nx_bbox  = int(1.5 * m_outradius / dx);
+    int ny_bbox  = int(1.5 * m_outradius / dy);
+
+    // Compute index range for x. We allow that indices are smaller than 0
+    // or equal or larger than nx to handle wrap around, which needs to be
+    // done in the client method.
+    if (2*nx_bbox < nx) {
+        ix1 = ix0 - nx_bbox;
+        ix2 = ix0 + nx_bbox;
+    }
+    else {
+        ix1 = 0;
+        ix2 = nx;
+    }
+
+    // Compute index range for y. The index range is restricted to [0,ny].
+    if (2*ny_bbox < ny) {
+        iy1 = iy0 - ny_bbox;
+        iy2 = iy0 + ny_bbox;
+        if (iy1 < 0) {
+            iy1 = 0;
         }
-        else {
-            map_exclusions_reg(filename);
+        if (iy2 > ny) {
+            iy2 = ny;
         }
-
-    } // endif: filename was valid
+    }
+    else {
+        iy1 = 0;
+        iy2 = ny;
+    }
     // Return
     return;
@@ -1083,64 +1440,166 @@ void ctskymap::map_exclusions(const GFilename&amp; filename)
 /***********************************************************************//**
- * @brief Fills exclusions map from FITS image
+ * @brief Return FFT kernel for background ring
  *
- * @param[in] filename FITS image file name.
+ * @param[in] map  Sky map to be convolved with a background ring.
+ * @param[in] rmin Minimum ring radius (degrees).
+ * @param[in] rmax Maximum ring radius (degrees).
+ * @return FFT kernel.
  *
- * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in
- * the exclusion sky map FITS file.
+ * Computes the FFT kernel for a backgrund ring.
  ***************************************************************************/
-void ctskymap::map_exclusions_fits(const GFilename& filename)
+GSkyMap ctskymap::ring_convolve(const GSkyMap& map, const double& rmin,
+                                                    const double& rmax) const
 {
-    // Load the fits image
-    GSkyMap inmap(filename);
+    // Copy input map
+    GSkyMap conv_map = map;
-    // Loop through the individual pixels in the exclusion map
-    for (int i = 0; i < m_exclmap.npix(); ++i) {
+    // Get FFT of ring kernel
+    GFft fft_kernel(ring_kernel(rmin, rmax));
-        // Get the pixel direction
-        GSkyDir dir = m_exclmap.pix2dir(i);
+    // Extract sky map into Ndarray
+    GNdarray array(conv_map.nx(), conv_map.ny());
+    double *dst = array.data();
+    for (int i = 0; i < conv_map.npix(); ++i) {
+        *dst++ = conv_map(i);
+    }
-        // Check this sky position in the fits map
-        if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) {
+    // FFT of sky map
+    GFft fft_array = GFft(array);
-            // Set the pixel to 1
-            m_exclmap(i) = 1.0;
+    // Multiply FFT of sky map with FFT of kernel
+    GFft fft_smooth = fft_array * fft_kernel;
-        } // endif: pixel,region overlap check
+    // Backward transform sky map
+    GNdarray smooth = fft_smooth.backward();
-    } // endfor: looped over exclusion map pixels
+    // Insert sky map
+    double *src = smooth.data();
+    for (int i = 0; i < conv_map.npix(); ++i) {
+        conv_map(i) = *src++;
+    }
-    // Return
-    return;
+    // Return convolved map
+    return conv_map;
 }
 /***********************************************************************//**
- * @brief Fills exclusions map from DS9 region file
+ * @brief Return FFT kernel for background ring
  *
- * @param[in] filename DS9 region file name.
+ * @param[in] rmin Minimum ring radius (degrees).
+ * @param[in] rmax Maximum ring radius (degrees).
+ * @return FFT kernel.
  *
- * Sets all exclusion map pixels to 1 that are contained in any of the DS9
- * regions.
+ * @exception GException::invalid_value
+ *            Sky map is not a WCS projection.
+ *
+ * Computes the FFT kernel for a background ring. The computation is done
+ * assuming that the sky map represents a cartesian grid, and distances are
+ * computed in that grid using Euclidean distances. The @p rmin and @p rmax
+ * arguments refer to a ring radius in these Euclidean distances.
  ***************************************************************************/
-void ctskymap::map_exclusions_reg(const GFilename& filename)
+GNdarray ctskymap::ring_kernel(const double& rmin, const double& rmax) const
 {
-    // Load the exclusion regions
-    GSkyRegions regions(filename);
+    // Get pointer on WCS projection
+    const GWcs* wcs = dynamic_cast<const GWcs*>(m_counts.projection());
+    if (wcs == NULL) {
+        std::string msg = "Sky map is not a WCS projection. Method is only "
+                          "valid for WCS projections.";
+        throw GException::invalid_value(G_RING_KERNEL, msg);
+    }
-    // Loop through the individual pixels in the exclusion map
-    for (int i = 0; i < m_exclmap.npix(); ++i) {
+    // Get X and Y step size
+    double dx = wcs->cdelt(0);
+    double dy = wcs->cdelt(1);
+
+    // Initialise kernel
+    GNdarray kern(m_counts.nx(), m_counts.ny());
+
+    // Fill kernel
+    for (int ix1 = 0, ix2 = m_counts.nx(); ix1 < m_counts.nx(); ++ix1, --ix2) {
+        double x   = ix1 * dx;
+        double xqs = x * x;
+        for (int iy1 = 0, iy2 = m_counts.ny(); iy1 < m_counts.ny(); ++iy1, --iy2) {
+            double y = iy1 * dy;
+            double r = std::sqrt(xqs + y*y);
+            if ((r >= rmin) && (r <= rmax)) {
+                kern(ix1,iy1) += 1.0;
+                if (ix2 < m_counts.nx()) {
+                    kern(ix2,iy1) += 1.0;
+                }
+                if (iy2 < m_counts.ny()) {
+                    kern(ix1,iy2) += 1.0;
+                }
+                if ((ix2 < m_counts.nx()) && (iy2 < m_counts.ny())) {
+                    kern(ix2,iy2) += 1.0;
+                }
+            }
+        }
+    }
-        // Get the pixel position
-        GSkyDir dir = m_exclmap.pix2dir(i);
+    // Return kernel
+    return kern;
+}
-        // If pixel position overlaps with the regions
-        if (regions.contains(dir)) {
-            m_exclmap(i) = 1.0;
-        }
-    } // endfor: looped over exclusion map pixels
+/***********************************************************************//**
+ * @brief Compute significance following Li & Ma
+ *
+ * @param[in] n_on Number of On-counts.
+ * @param[in] n_off Number of Off-counts.
+ * @param[in] alpha Alpha parameters.
+ * @return Significance.
+ *
+ * Computes the significance following Li & Ma, Equation (17).
+ ***************************************************************************/
+double ctskymap::sigma_li_ma(const double& n_on,
+                             const double& n_off,
+                             const double& alpha) const
+{
+    // Allocate result
+    double sigma;
+
+    // Handle special case of no On-counts
+    if (n_on == 0.0) {
+        sigma = -2.0 * n_off * std::log(1.0+alpha);
+    }
+
+    // ... otherwise handle general case
+    else {
+        sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
+                (n_on  * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) +
+                 n_off * std::log((1.0+alpha) * n_off / (n_on+n_off)));
+    }
+
+    // Return sigma
+    return sigma;
+}
+
+
+/***********************************************************************//**
+ * @brief Write sky map into FITS file
+ *
+ * @param[in,out] fits FITS file.
+ * @param[in] map Sky map.
+ * @param[in] extname Extension name.
+ *
+ * Write one sky map with the extension name and keywords into the FITS file.
+ ***************************************************************************/
+void ctskymap::write_map(GFits& fits, const GSkyMap& map, const std::string& extname) const
+{
+    // Write map into FITS file
+    GFitsHDU* hdu = map.write(fits);
+            
+    // Set map extension name
+    if (hdu != NULL) {
+        hdu->extname(extname);
+    }
+
+    // Write keywords into map extension
+    write_ogip_keywords(hdu);
+    write_hdu_keywords(hdu);
     // Return
     return;
@@ -1148,21 +1607,36 @@ void ctskymap::map_exclusions_reg(const GFilename&amp; filename)
 /***********************************************************************//**
- * @brief Estimates the background in sky map based on IRF template
+ * @brief Write keywords in FITS HDU
  *
- * @param[in] lnE Natural logarithm of energy in MeV.
+ * @param[in,out] hdu Pointer to FITS HDU.
+ *
+ * Writes keywords in FITS HDU.
  ***************************************************************************/
-double ctskymap::irf_kern::eval(const double& lnE)
+void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
 {
-    // Get log10 of energy in TeV
-    double logE = lnE * gammalib::inv_ln10 - 6;
+    // Continue only if pointer is valid
+    if (hdu != NULL) {
+
+        // Set keywords
+        hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
+        hdu->card("E_MIN",  m_emin, "[TeV] Lower energy boundary");
+        hdu->card("E_MAX",  m_emax, "[TeV] Upper energy boundary");
+        hdu->card("EUNIT",  "TeV",  "Units for E_MIN and E_MAX");
-    // Get function value
-    double value = (*m_bgd)(logE, m_dir->detx(), m_dir->dety());
+        // Write RING background keywords
+        if (m_bkgsubtract == "RING") {
+            hdu->card("USEFFT",  m_usefft, "Use FFT for RING background");
+            hdu->card("ROIRAD",  m_roiradius, "[deg] Source region radius");
+            hdu->card("INRAD",   m_inradius,  "[deg] Inner background ring radius");
+            hdu->card("OUTRAD",  m_outradius, "[deg] Outer background ring radius");
+            hdu->card("ITER",    m_iterations, "Exclusion map iterations");
+            hdu->card("THRES",   m_threshold,  "Exclusion map threshold");
+            hdu->card("EXCLMAP", m_inexclusion.url(),  "Exclusion map name");
+        }
-    // Correct for variable substitution
-    value *= std::exp(lnE);
+    } // endif: pointer was valid
-    // Return value
-    return value;
+    // Return
+    return;
 }
 /***************************************************************************
  *                       ctskymap - Sky mapping tool                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -64,59 +64,69 @@ public:
 protected:
     // Protected methods
-    void init_members(void);
-    void copy_members(const ctskymap& app);
-    void free_members(void);
-    void get_parameters(void);
-    void map_events(GCTAObservation* obs);
-    void map_background(GCTAObservation* obs);
-    void map_background_irf(GCTAObservation* obs);
-    void map_background_ring(GCTAObservation* obs);
-    void map_significance(void);
-    void map_exclusions(const GFilename& filename);
-    void map_exclusions_fits(const GFilename& filename);
-    void map_exclusions_reg(const GFilename& filename);
-    void compute_ring_values(const GSkyMap& counts, 
-                             const GSkyMap& sensitivity,
-                             const GSkyDir& position,
-                             double& non, double& noff, double& alpha);
-    void write_hdu_keywords(GFitsHDU* hdu) const;
-
-    // Background integration kernel
-    class irf_kern : public GFunction {
-    public:
-        irf_kern(const GCTABackground* bgd,
-                 const GCTAInstDir*    dir) :
-                 m_bgd(bgd),
-                 m_dir(dir) { }
-        double eval(const double& lnE);
-    protected:
-        const GCTABackground* m_bgd;   //!< Pointer to background
-        const GCTAInstDir*    m_dir;   //!< Pointer to instrument direction
-    };
+    void     init_members(void);
+    void     copy_members(const ctskymap& app);
+    void     free_members(void);
+    void     get_parameters(void);
+    void     setup_maps(void);
+    void     setup_exclusion_map(const GFilename& filename);
+    void     setup_exclusion_map_fits(const GFilename& filename);
+    void     setup_exclusion_map_region(const GFilename& filename);
+    void     fill_maps(void);
+    void     fill_maps_counts(GCTAObservation* obs);
+    void     fill_maps_acceptance(GCTAObservation* obs);
+    void     compute_maps(void);
+    void     compute_maps_ring_fft(void);
+    void     compute_maps_ring_direct(void);
+    void     compute_ring_values(const int&     ipixel,
+                                 const GSkyMap& counts,
+                                 const GSkyMap& background,
+                                 double&        non,
+                                 double&        noff,
+                                 double&        alpha);
+    void     ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
+                                                  int& iy1, int& iy2) const;
+    GSkyMap  ring_convolve(const GSkyMap& map, const double& rmin,
+                                               const double& rmax) const;
+    GNdarray ring_kernel(const double& rmin, const double& rmax) const;
+    double   sigma_li_ma(const double& n_on,
+                         const double& n_off,
+                         const double& alpha) const;
+    void     write_map(GFits&             fits,
+                       const GSkyMap&     map,
+                       const std::string& extname) const;
+    void     write_hdu_keywords(GFitsHDU* hdu) const;
     // User parameters
     GFilename     m_outmap;      //!< Output file name
+    GFilename     m_inexclusion; //!< Exclusion map file name
     double        m_emin;        //!< Minimum energy (TeV)
     double        m_emax;        //!< Maximum energy (TeV)
     std::string   m_bkgsubtract; //!< Background subtraction method
     double        m_roiradius;   //!< Region of interest radius for RING bkg.
-    double        m_inradius;    //!< Inner ring radius for RING bkg.
-    double        m_outradius;   //!< Outer ring radius for RING bkg.
+    double        m_inradius;    //!< Inner ring radius for RING background
+    double        m_outradius;   //!< Outer ring radius for RING background
+    int           m_iterations;  //!< Number of iterations for RING background
+    double        m_threshold;   //!< Threshold for RING background
+    bool          m_usefft;      //!< Use FFT for RING background
     bool          m_publish;     //!< Publish sky map?
     GChatter      m_chatter;     //!< Chattiness
     // Protected members
+    bool          m_has_inmap;   //!< Has valid input map
     GSkyMap       m_skymap;      //!< Sky map
     GSkyMap       m_bkgmap;      //!< Background map
     GSkyMap       m_sigmap;      //!< Significance map
+    GSkyMap       m_counts;      //!< Counts map
+    GSkyMap       m_acceptance;  //!< Acceptance map
     GSkyMap       m_exclmap;     //!< Exclusion map for RING background
-    GSkyMap       m_alphamap;    //!< Alpha map values for RING background
-    GSkyMap       m_onmap;       //!< On counts map for RING background
     // Caching variables to prevent multiple computations of the same thing
-    std::vector<double>  m_solidangle; //!< Cached pixel solid angles
-    std::vector<GSkyDir> m_dirs;       //!< Cached pixel directions
+    std::vector<double>  m_solidangle;    //!< Cached pixel solid angles
+    std::vector<GSkyDir> m_dirs;          //!< Cached pixel directions
+    double               m_cos_roiradius; //!< Cosine of RoI radius
+    double               m_cos_inradius;  //!< Cosine of inner ring radius
+    double               m_cos_outradius; //!< Cosine of outer ring radius
 };
@@ -19,6 +19,7 @@
 inobs,  f, a, events.fits,,, "Input event list or observation definition XML file"
 caldb,  s, a, prod2,,, "Calibration database"
 irf,    s, a, South_0.5h,,, "Instrument response function"
+inmap,  f, h, NONE,,, "Input skymap file"
 outmap, f, a, skymap.fits,,, "Output skymap file"
 #
@@ -43,10 +44,13 @@ yref,     r, a, 22.01,-90,90, &quot;Second coordinate of image center in degrees (DEC
 # Background subtraction parameters
 #==================================
 bkgsubtract, s, a, NONE, NONE|IRF|RING,, "Background subtraction method"
-roiradius,   r, a, 0.2,,, "Source region radius for estimating on-counts (degrees)"
-inradius,    r, a, 0.6,,, "Inner ring radius (degrees)"
-outradius,   r, a, 0.8,,, "Outer ring radius (degrees)"
+roiradius,   r, a, 0.1,,, "Source region radius for estimating on-counts (degrees)"
+inradius,    r, a, 0.6,,, "Inner background ring radius (degrees)"
+outradius,   r, a, 0.8,,, "Outer background ring radius (degrees)"
+iterations,  i, a, 0,0,100, "Number of iterations for exclusion regions computation"
+threshold,   r, a, 5.0,,, "Significance threshold for exclusion regions computation"
 inexclusion, f, h, NONE,,, "Exclusion map or DS9 region file name"
+usefft,      b, h, yes,,, "Use FFT for ring computation?"
 #
 # Standard parameters
@@ -319,6 +319,48 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const
 /***********************************************************************//**
+ * @brief Read OGIP keywords from FITS HDU
+ *
+ * @param[in,out] hdu Pointer to FITS HDU.
+ *
+ * Read OGIP keywords from FITS HDU.
+ ***************************************************************************/
+void ctobservation::read_ogip_keywords(GFitsHDU* hdu) const
+{
+    // Continue only if pointer is valid
+    if (hdu != NULL) {
+
+        // Get observation information
+        m_ogip_telescope = (hdu->has_card("TELESCOP")) ? hdu->string("TELESCOP") : "";
+
+        // Get observation time information
+        std::string date_obs = (hdu->has_card("DATE-OBS")) ? hdu->string("DATE-OBS") : "";
+        std::string time_obs = (hdu->has_card("TIME-OBS")) ? hdu->string("TIME-OBS") : "";
+        std::string date_end = (hdu->has_card("DATE-END")) ? hdu->string("DATE-END") : "";
+        std::string time_end = (hdu->has_card("TIME-END")) ? hdu->string("TIME-END") : "";
+        m_ogip_telapse       = (hdu->has_card("TELAPSE"))  ? hdu->real("TELAPSE") : 0.0;
+        m_ogip_ontime        = (hdu->has_card("ONTIME"))   ? hdu->real("ONTIME") : 0.0;
+        m_ogip_livetime      = (hdu->has_card("LIVETIME")) ? hdu->real("LIVETIME") : 0.0;
+        m_ogip_exposure      = (hdu->has_card("EXPOSURE")) ? hdu->real("EXPOSURE") : 0.0;
+
+        // Get OGIP start and stop time
+        std::string tstart = date_obs + "T" + time_obs;
+        std::string tstop  = date_end + "T" + time_end;
+        if (tstart != "T") {
+            m_ogip_tstart.utc(tstart);
+        }
+        if (tstop != "T") {
+            m_ogip_tstop.utc(tstop);
+        }
+
+    } // endif: pointer was valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Write OGIP keywords in FITS HDU
  *
  * @param[in,out] hdu Pointer to FITS HDU.
 /***************************************************************************
  *             ctobservation - Base class for observation tools            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2016-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2016-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -73,6 +73,7 @@ protected:
     GCTAObservation*       next_unbinned_observation(void);
     const GCTAObservation* first_unbinned_observation(void) const;
     const GCTAObservation* next_unbinned_observation(void) const;
+    void                   read_ogip_keywords(GFitsHDU* hdu) const;
     void                   write_ogip_keywords(GFitsHDU* hdu) const;
     void                   set_obs_statistic(const std::string& statistic);
     void                   save_events_fits(void);
@@ -650,7 +650,11 @@ GEbounds ctool::create_ebounds(void)
  ***************************************************************************/
 GSkyMap ctool::create_map(const GObservations& obs)
 {
-    // Read task parameters
+    // Read coordinate system and projection
+    std::string coordsys = (*this)["coordsys"].string();
+    std::string proj     = (*this)["proj"].string();
+
+    // Read sky map reference
     double xref   = 0.0;
     double yref   = 0.0;
     bool   usepnt = (*this)["usepnt"].boolean();
@@ -658,11 +662,11 @@ GSkyMap ctool::create_map(const GObservations&amp; obs)
         xref = (*this)["xref"].real();
         yref = (*this)["yref"].real();
     }
-    std::string proj     = (*this)["proj"].string();
-    std::string coordsys = (*this)["coordsys"].string();
-    double      binsz    = (*this)["binsz"].real();
-    int         nxpix    = (*this)["nxpix"].integer();
-    int         nypix    = (*this)["nypix"].integer();
+
+    // Read sky map pixel size and number of bins
+    double binsz = (*this)["binsz"].real();
+    int    nxpix = (*this)["nxpix"].integer();
+    int    nypix = (*this)["nypix"].integer();
     // If requested, get pointing from observations
     if (usepnt) {
@@ -112,7 +112,7 @@ class Test(test):
         obsinfo['ra']         = 83.63
         obsinfo['dec']        = 22.01
         obsinfo['logfile']    = 'csobsinfo_py1.log'
-        obsinfo['chatter']    = 2
+        obsinfo['chatter']    = 3
         # Run csobsinfo script and save XML file
         obsinfo.logFileOpen()   # Make sure we get a log file
@@ -122,21 +122,20 @@ class Test(test):
         # Check DS9 file
         self._check_ds9_file('csobsinfo_py1.reg')
-        # Set-up csobsinfo
-        obsinfo = cscripts.csobsinfo()
-        obsinfo['inobs']      = self._inobs
-        obsinfo['outds9file'] = 'csobsinfo_py2.reg'
-        obsinfo['offset']     = True
-        obsinfo['ra']         = 83.63
-        obsinfo['dec']        = 22.01
-        obsinfo['logfile']    = 'csobsinfo_py2.log'
-        obsinfo['chatter']    = 3
-
-        # Execute csobsinfo script
-        obsinfo.execute()
-
-        # Check DS9 file
-        self._check_ds9_file('csobsinfo_py2.reg')
+        # Check result access methods
+        self.test_value(len(obsinfo.zeniths()), 1, 'Test number of zenith angles')
+        self.test_value(obsinfo.zeniths()[0], 0.0, 'Test zenith angle')
+        self.test_value(len(obsinfo.azimuths()), 1, 'Test number of azimuth angles')
+        self.test_value(obsinfo.azimuths()[0], 0.0, 'Test azimuth angle')
+        self.test_value(len(obsinfo.ras()), 1, 'Test number of Right Ascensions')
+        self.test_value(obsinfo.ras()[0], 83.63, 'Test Right Ascension')
+        self.test_value(len(obsinfo.decs()), 1, 'Test number of Declinations')
+        self.test_value(obsinfo.decs()[0], 22.51, 'Test Declination')
+        self.test_value(len(obsinfo.offsets()), 1, 'Test number of offsets')
+        self.test_value(obsinfo.offsets()[0], 0.5, 'Test offset angle')
+        self.test_value(obsinfo.ebounds().size(), 1, 'Test number of energy boundaries')
+        self.test_value(obsinfo.ebounds().emin().GeV(), 100.0, 'Test minimum energy')
+        self.test_value(obsinfo.ebounds().emax().TeV(), 100.0, 'Test maximum energy')
         # Return
         return
@@ -151,13 +150,12 @@ class Test(test):
         # Expect "fk5" in first line
         line = f.readline().strip('\n')
-        self.test_assert(line == 'fk5',
-                         'Check for "fk5" in first line')
+        self.test_value(line, 'fk5',
+                         'Test for "fk5" in first line of DS9 file')
         # Expect "fk5" in first line
         line = f.readline().strip('\n')
-        self.test_assert(line[0:6] == 'point(',
-                         'Check for "point(" in second line')
+        self.test_value(line[0:6], 'point(', 'Test for "point(" in second line of DS9 file')
         # Close file
         f.close()
@@ -248,6 +248,7 @@ class Test(test):
         skymap['roiradius']   = 0.1
         skymap['inradius']    = 0.6
         skymap['outradius']   = 0.8
+        skymap['iterations']  = 0
         skymap['inexclusion'] = 'NONE'
         skymap['caldb']       = self._caldb
         skymap['irf']         = self._irf
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the obsutils module
 #
-# Copyright (C) 2017 Juergen Knoedlseder
+# Copyright (C) 2017-2018 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
@@ -214,7 +214,7 @@ class Test(test):
         # Check simulation results
         self.test_value(res.size(), 1, 'Check number of observations')
         self.test_value(res.models().size(), 2, 'Check number of models')
-        self.test_value(res.nobserved(), 24, 'Check number of observed events')
+        self.test_value(res.nobserved(), 23, 'Check number of observed events')
         self.test_value(res.npred(), 0.0, 'Check number of predicted events')
         self.test_value(res[0].eventtype(), 'CountsCube', 'Check event type')
         self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
@@ -223,7 +223,7 @@ class Test(test):
                         'Check minimum energy')
         self.test_value(res[0].events().ebounds().size(), 5,
                         'Check number of energy bins')
-        self.test_value(res[0].events().number(), 24,
+        self.test_value(res[0].events().number(), 23,
                         'Check number of events in cube')
         # Check energy dispersion flag