Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents a0d4125d d5430cfa
Showing 17 changed files with 557 additions and 163 deletions
cscripts/Makefile.am
cscripts/cslightcrv.par
cscripts/cslightcrv.py
cscripts/csphasecrv.par
cscripts/csphasecrv.py
cscripts/obsutils.py
doc/source/users/reference_manual/cslightcrv.rst
doc/source/users/reference_manual/csphasecrv.rst
doc/source/users/tutorials/howto/howto_lightcurve.rst
pyext/ctool.i
src/support/ctool.cpp
src/support/ctool.hpp
test/data/crab_offaxis1.fits.gz
test/data/crab_onoff.xml
test/test_cscript.py
test/test_cslightcrv.py
test/test_csphasecrv.py
@@ -164,6 +164,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/crab_psfcube.fits.gz \
            $(top_srcdir)/test/data/crab_edispcube.fits.gz \
            $(top_srcdir)/test/data/crab_bkgcube.fits.gz \
+           $(top_srcdir)/test/data/crab_onoff.xml \
            $(top_srcdir)/test/data/crab_offaxis1.fits.gz \
            $(top_srcdir)/test/data/crab_offaxis2.fits.gz \
            $(top_srcdir)/test/data/crab_exclusion.fits.gz \
@@ -16,12 +16,13 @@
 #
 # Input and output parameters
 #============================
-inobs,     f, a, events.fits,,, "Input event list or observation definition XML file"
-inmodel,   f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
-srcname,   s, a, Crab,,, "Source name"
-caldb,     s, a, prod2,,, "Calibration database"
-irf,       s, a, South_0.5h,,, "Instrument response function"
-edisp,     b, h, no,,, "Apply energy dispersion?"
+inobs,        f, a, events.fits,,, "Input event list or observation definition XML file"
+inmodel,      f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
+srcname,      s, a, Crab,,, "Source name"
+caldb,        s, a, prod2,,, "Calibration database"
+irf,          s, a, South_0.5h,,, "Instrument response function"
+inexclusion,  f, h, NONE,,, "Sky map of exclusion regions for On/Off analysis"
+edisp,     b, h, no,,, "Apply energy dispersion in Cube analysis?"
 outfile,   f, a, lightcurve.fits,,, "Output light curve file"
 #
@@ -32,16 +33,25 @@ tmin,      t, a, MJD 51544.5,,,&quot;Lightcurve start time (UTC string, JD, MJD or ME
 tmax,      t, a, MJD 51544.6,,,"Lightcurve stop time (UTC string, JD, MJD or MET in seconds)"
 tbins,     i, a, 5,1,10000, "Number of time bins"
 tbinfile,  f, a, NONE,,, "File defining the time binning"
+method,    s, a, 3D, 3D|ONOFF,, "Analysis method"
 emin,      r, a, 0.1,,, "Lower energy limit of events (TeV)"
 emax,      r, a, 100.0,,, "Upper energy limit of events (TeV)"
-enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)"
+enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned for 3D analysis only)"
 coordsys,  s, a, CEL,CEL|GAL,, "Coordinate System"
 proj,      s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
-xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
-yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+xref,      r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
+yref,      r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
 nypix,     i, a, 200,,, "Size of the Y axis in pixels"
 binsz,     r, a, 0.02,,, "Pixel size (deg/pixel)"
+srcshape,  s, h, CIRCLE,CIRCLE,, "Shape of the source region"
+rad,       r, a, 0.2,0,180, "Radius of source region circle (deg)"
+bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
+bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
+maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
+etruemin,  r, h, 0.01,,, "Minimum true energy (TeV)"
+etruemax,  r, h, 300.0,,, "Maximum true energy (TeV)"
+etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
 statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
 calc_ts,   b, h, yes,,, "Compute TS value for each bin?"
 calc_ulim, b, h, yes,,, "Compute upper limit for each bin?"
@@ -74,6 +74,7 @@ class cslightcrv(ctools.csobservation):
         # Initialise some members
         self._srcname = ''
         self._tbins   = gammalib.GGti()
+        self._onoff   = False
         self._stacked = False
         self._fits    = gammalib.GFits()
@@ -103,8 +104,13 @@ class cslightcrv(ctools.csobservation):
         # Get time boundaries
         self._tbins = self._create_tbounds()
+        # Set On/Off analysis flag and query relevant user parameters
+        self._onoff = self._is_onoff()
+
+        # If cube analysis is selected
         # Set stacked analysis flag and query relevant user parameters
-        self._stacked = self._is_stacked()
+        if not self._onoff:
+            self._stacked = self._is_stacked()
         # Query the hidden parameters, just in case
         self['edisp'].boolean()
@@ -357,40 +363,6 @@ class cslightcrv(ctools.csobservation):
         # Return upper limit tuple
         return ul_diff, ul_flux, ul_eflux
-    def _bin_observation(self, obs):
-        """
-        Bin an observation if a binned analysis was requested
-
-        Parameters
-        ----------
-        obs : `~gammalib.GObservations`
-            Observation container
-
-        Returns
-        -------
-        obs : `~gammalib.GObservations`
-            Observation container where the first observation is a binned observation
-        """
-        # Get new observation container with binned observation
-        new_obs = obsutils.get_stacked_obs(self, obs)
-        
-        # Extract models
-        models = new_obs.models()
-
-        # Fix background models if required
-        if self['fix_bkg'].boolean():
-            for model in models:
-                if model.classname() != 'GModelSky':
-                    for par in model:
-                        par.fix()
-
-        # Put back models
-        new_obs.models(models)
-
-        # Return new oberservation container
-        return new_obs
-
-
     # Public methods
     def run(self):
         """
@@ -483,31 +455,116 @@ class cslightcrv(ctools.csobservation):
             # Retrieve observation
             obs = select.obs()
-            # If a stacked analysis is requested then bin the events
-            # and compute the stacked response functions and setup
-            # an observation container with a single stacked observation.
-            if self._stacked:
-                obs = self._bin_observation(obs)
+            # Deal with stacked and On/Off Observations
+            if self._stacked or self._onoff:
+
+                # If a stacked analysis is requested bin the events
+                # and compute the stacked response functions and setup
+                # an observation container with a single stacked observation.
+                if self._stacked:
+                    new_obs = obsutils.get_stacked_obs(self, obs)
+
+                # ... otherwise if On/Off analysis is requested generate
+                # the On/Off observations and response
+                elif self._onoff:
+                    new_obs = obsutils.get_onoff_obs(self, obs)
+
+                # Extract models
+                models = new_obs.models()
+
+                # Fix background models if required
+                if self['fix_bkg'].boolean():
+                    for model in models:
+                        if model.classname() != 'GModelSky':
+                            for par in model:
+                                par.fix()
+
+                # Put back models
+                new_obs.models(models)
+
+                # Continue with new oberservation container
+                obs = new_obs
             # Header
             self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
             # Do maximum likelihood model fitting
-            like = ctools.ctlike(obs)
-            like['edisp'] = self['edisp'].boolean()
-            like.run()
+            if obs.size() > 0:
+                like = ctools.ctlike(obs)
+                like['edisp'] = self['edisp'].boolean()
+                like.run()
+
+                # Skip bin if no event was present
+                if like.obs().logL() == 0.0:
-            # Skip bin if no event was present
-            if like.obs().logL() == 0.0:
+                    # Signal skipping of bin
+                    self._log_value(gammalib.TERSE, 'Warning',
+                                    'No event in this time bin, skip bin.')
-                # Signal skipping of bin
+                    # Set all results to 0
+                    for par in pars:
+                        result['values'][par]      = 0.0
+                        result['values']['e_'+par] = 0.0
+
+                    # Append result
+                    results.append(result)
+
+                    # Continue with next time bin
+                    continue
+
+                # Retrieve model fitting results for source of interest
+                source = like.obs().models()[self._srcname]
+
+                # Extract parameter values
+                for par in pars:
+                    result['values'][par]      = source[par].value()
+                    result['values']['e_'+par] = source[par].error()
+
+                # Calculate upper limit (-1 if not computed)
+                #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs())
+                ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs)
+                if ul_diff > 0.0:
+                    result['ul_diff']  = ul_diff
+                    result['ul_flux']  = ul_flux
+                    result['ul_eflux'] = ul_eflux
+
+                # Extract Test Statistic value
+                if self['calc_ts'].boolean():
+                    result['ts'] = source.ts()
+
+                # Append result to list of dictionaries
+                results.append(result)
+
+                # Log results for this time bin
+                self._log.header3('Results')
+                pars = self._get_free_par_names()
+                for par in pars:
+                    value = source[par].value()
+                    error = source[par].error()
+                    unit  = source[par].unit()
+                    self._log_value(gammalib.NORMAL, par,
+                                    str(value)+' +/- '+str(error)+' '+unit)
+                if ul_diff > 0.0:
+                    self._log_value(gammalib.NORMAL, 'Upper flux limit',
+                                    str(result['ul_diff'])+' ph/cm2/s/MeV')
+                    self._log_value(gammalib.NORMAL, 'Upper flux limit',
+                                    str(result['ul_flux'])+' ph/cm2/s')
+                    self._log_value(gammalib.NORMAL, 'Upper flux limit',
+                                    str(result['ul_eflux'])+' erg/cm2/s')
+                if self['calc_ts'].boolean():
+                    self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts'])
+
+            # Otherwise, if observations size is 0, signal bin is skipped and
+            # fill results table with zeros
+            else:
                 self._log_value(gammalib.TERSE, 'Warning',
-                                'No event in this time bin, skip bin.')
+                                'No observations available in this time bin',
+                                'Skip bin.')
                 # Set all results to 0
                 for par in pars:
-                    result['values'][par]      = 0.0
-                    result['values']['e_'+par] = 0.0
+                    result['values'][par] = 0.0
+                    result['values']['e_' + par] = 0.0
                 # Append result
                 results.append(result)
@@ -515,48 +572,6 @@ class cslightcrv(ctools.csobservation):
                 # Continue with next time bin
                 continue
-            # Retrieve model fitting results for source of interest
-            source = like.obs().models()[self._srcname]
-
-            # Extract parameter values
-            for par in pars:
-                result['values'][par]      = source[par].value()
-                result['values']['e_'+par] = source[par].error()
-
-            # Calculate upper limit (-1 if not computed)
-            #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs())
-            ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs)
-            if ul_diff > 0.0:
-                result['ul_diff']  = ul_diff
-                result['ul_flux']  = ul_flux
-                result['ul_eflux'] = ul_eflux
-
-            # Extract Test Statistic value
-            if self['calc_ts'].boolean():
-                result['ts'] = source.ts() 
-
-            # Append result to list of dictionaries
-            results.append(result)
-
-            # Log results for this time bin
-            self._log.header3('Results')
-            pars = self._get_free_par_names()
-            for par in pars:
-                value = source[par].value()
-                error = source[par].error()
-                unit  = source[par].unit()
-                self._log_value(gammalib.NORMAL, par,
-                                str(value)+' +/- '+str(error)+' '+unit)
-            if ul_diff > 0.0:
-                self._log_value(gammalib.NORMAL, 'Upper flux limit',
-                                str(result['ul_diff'])+' ph/cm2/s/MeV')
-                self._log_value(gammalib.NORMAL, 'Upper flux limit',
-                                str(result['ul_flux'])+' ph/cm2/s')
-                self._log_value(gammalib.NORMAL, 'Upper flux limit',
-                                str(result['ul_eflux'])+' erg/cm2/s')
-            if self['calc_ts'].boolean():
-                self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts'])
-
         # Create FITS table from results
         table = self._create_fits_table(results)
@@ -21,7 +21,8 @@ inmodel,   f, a, $CTOOLS/share/models/crab.xml,,, &quot;Input model definition XML fi
 srcname,   s, a, Crab,,, "Source name"
 caldb,     s, a, prod2,,, "Calibration database"
 irf,       s, a, South_0.5h,,, "Instrument response function"
-edisp,     b, h, no,,, "Apply energy dispersion?"
+inexclusion,  f, h, NONE,,, "Sky map of exclusion regions for On/Off analysis"
+edisp,     b, h, no,,, "Apply energy dispersion in Cube analysis?"
 outfile,   f, a, phasecrv.fits,,, "Output fits file name."
 #
@@ -30,16 +31,25 @@ outfile,   f, a, phasecrv.fits,,, &quot;Output fits file name.&quot;
 phbinalg,  s, a, LIN,LIN|FILE,, "Algorithm for defining phase bins"
 phbins,    i, a, 5,1,10000, "Number of phase bins"
 phbinfile, f, a, NONE,,, "File defining the phase binning"
+method,    s, a, 3D, 3D|ONOFF,, "Analysis method"
 emin,      r, a, 0.1,,, "Lower energy limit of events (TeV)"
 emax,      r, a, 100.0,,, "Upper energy limit of events (TeV)"
-enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)"
+enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned for 3D analysis only))"
 coordsys,  s, a, CEL,CEL|GAL,, "Coordinate System"
 proj,      s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
-xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
-yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+xref,      r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
+yref,      r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
 nypix,     i, a, 200,,, "Size of the Y axis in pixels"
 binsz,     r, a, 0.02,,, "Pixel size (deg/pixel)"
+srcshape,  s, h, CIRCLE,CIRCLE,, "Shape of the source region"
+rad,       r, a, 0.2,0,180, "Radius of source region circle (deg)"
+bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
+bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
+maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
+etruemin,  r, h, 0.01,,, "Minimum true energy (TeV)"
+etruemax,  r, h, 300.0,,, "Maximum true energy (TeV)"
+etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
 statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
 #
@@ -63,6 +63,7 @@ class csphasecrv(ctools.csobservation):
         # [[ph1min,ph1max], [ph2min,ph2max],..]
         self._srcname   = ''
         self._phbins    = [[0.0,1.0]]
+        self._onoff     = False
         self._stacked   = False
         self._fits      = gammalib.GFits()
         self._fitmodels = {}
@@ -92,8 +93,13 @@ class csphasecrv(ctools.csobservation):
         # Get phase boundaries
         self._phbins = self._create_tbounds()
+        # Set On/Off analysis flag and query relevant user parameters
+        self._onoff = self._is_onoff()
+
+        # If cube analysis is selected
         # Set stacked analysis flag and query relevant user parameters
-        self._stacked = self._is_stacked()
+        if not self._onoff:
+            self._stacked = self._is_stacked()
         # Query the hidden parameters, just in case
         self['edisp'].boolean()
@@ -236,11 +242,20 @@ class csphasecrv(ctools.csobservation):
             # Store fit results
             phname = str(phmin)+'-'+str(phmax)
-            source = self._fitmodels[phname][self._srcname]
-            for par in pars:
-                result['values'][par]      = source[par].value()
-                result['values']['e_'+par] = source[par].error()
-            
+
+            # If the model contains the source of interest fill results
+            try:
+                source = self._fitmodels[phname][self._srcname]
+                for par in pars:
+                    result['values'][par]      = source[par].value()
+                    result['values']['e_'+par] = source[par].error()
+
+            # ... otherwise fills with zeros
+            except:
+                for par in pars:
+                    result['values'][par]      = 0.
+                    result['values']['e_'+par] = 0.
+
             # Append result to list of dictionaries
             results.append(result)
@@ -348,36 +363,57 @@ class csphasecrv(ctools.csobservation):
             select['dec']  = 'UNDEFINED'
             select['expr'] = 'PHASE>'+str(phbin[0])+' && PHASE<'+str(phbin[1])
             select.run()
-            
+
+            # Set phase string
+            phstr = str(phbin[0]) + '-' + str(phbin[1])
+
             # Add phase to observation id
             for i in range(0,select.obs().size()):
                 oldid = select.obs()[i].id()
-                select.obs()[i].id(oldid+'_'+str(phbin[0])+'-'+str(phbin[1]))
+                select.obs()[i].id(oldid+'_'+phstr)
             obs = select.obs()
-            # If a stacked analysis is requested then bin the events
-            # and compute the stacked response functions and setup
+            # If an On/Off analysis is requested generate the On/Off observations
+            if self._onoff:
+                obs = obsutils.get_onoff_obs(self, select.obs())
+
+            # ... otherwise, if stacked analysis is requested then bin the
+            # events and compute the stacked response functions and setup
             # an observation container with a single stacked observation.
-            if self._stacked:
+            elif self._stacked:
                 obs = obsutils.get_stacked_obs(self, select.obs())
-    
+
             # Header
             self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
-    
-            # Do maximum likelihood model fitting
-            like = ctools.ctlike(obs)
-            like['edisp'] = self['edisp'].boolean()
-            like.run()
-            
-            # Renormalize models to phase selection
-            # TODO move the scaling from the temporal to the spectral component
-            for model in like.obs().models():
-                scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0])
-                model['Normalization'].value(scaled_norm)
-            
-            # Store fit model
-            self._fitmodels[str(phbin[0])+'-'+str(phbin[1])] = \
-                            like.obs().models().copy()
+
+            # The On/Off analysis can produce empty observation containers,
+            # e.g., when on-axis observations are used. To avoid ctlike asking
+            # for a new observation container (or hang, if in interactive mode)
+            # we'll run ctlike only if the size is >0
+            if obs.size() > 0:
+
+                # Do maximum likelihood model fitting
+                like = ctools.ctlike(obs)
+                like['edisp'] = self['edisp'].boolean()
+                like.run()
+
+                # Renormalize models to phase selection
+                # TODO move the scaling from the temporal to the spectral component
+                for model in like.obs().models():
+                    scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0])
+                    model['Normalization'].value(scaled_norm)
+
+                # Store fit model
+                self._fitmodels[phstr] = like.obs().models().copy()
+
+            # ... otherwise we set an empty model container
+            else:
+                self._log_string(gammalib.TERSE,
+                                 'PHASE %f - %f: no observations available'
+                                 ' for fitting' %(phbin[0], phbin[1]))
+
+                # Set empty models container
+                self._fitmodels[phstr] = gammalib.GModels()
         # Create FITS file
         self._create_fits()
@@ -760,3 +760,71 @@ def get_stacked_obs(cls, obs):
     # Return new oberservation container
     return new_obs
+
+
+# ================================ #
+# Get On/Off observation container #
+# ================================ #
+def get_onoff_obs(cls, obs):
+    """
+    Create On/Off observations container from given observations
+
+    Parameters
+    ----------
+    cls : `~ctools.cscript`
+        cscript class
+    obs : `~gammalib.GObservations`
+        Observation container
+
+    Returns
+    -------
+    onoff_obs : `~gammalib.GObservations`
+        Observation container with On/Off observations
+    """
+    # Write header
+    if cls._logExplicit():
+        cls._log.header3('Creating On/Off observations')
+
+    # Create On/Off observations
+    phagen = cscripts.csphagen(obs)
+    phagen['inexclusion'] = cls['inexclusion'].value()
+    phagen['emin']        = cls['emin'].real()
+    phagen['emax']        = cls['emax'].real()
+    phagen['enumbins']    = cls['enumbins'].integer()
+    phagen['ebinalg']     = 'LOG'
+    phagen['srcshape']    = cls['srcshape'].string()
+    phagen['coordsys']    = cls['coordsys'].string()
+    if cls['coordsys'].string() == 'CEL':
+        phagen['ra']  = cls['xref'].real()
+        phagen['dec'] = cls['yref'].real()
+    elif cls['coordsys'].string() == 'GAL':
+        phagen['glon'] = cls['xref'].real()
+        phagen['glat'] = cls['yref'].real()
+    if cls['srcshape'].string() == 'CIRCLE':
+        phagen['rad'] = cls['rad'].real()
+    phagen['bkgmethod'] = cls['bkgmethod'].string()
+    if cls['bkgmethod'].string() == 'REFLECTED':
+        phagen['bkgregmin'] = cls['bkgregmin'].integer()
+    phagen['maxoffset'] = cls['maxoffset'].real()
+    phagen['stack']     = True
+    phagen['etruemin']  = cls['etruemin'].real()
+    phagen['etruemax']  = cls['etruemax'].real()
+    phagen['etruebins'] = cls['etruebins'].integer()
+    phagen['chatter']   = cls['chatter'].integer()
+    phagen['clobber']   = cls['clobber'].boolean()
+    phagen['debug']     = cls['debug'].boolean()
+    phagen.run()
+
+    # Clone resulting observation container
+    onoff_obs = phagen.obs().copy()
+
+    # On/Off observations are created with CSTAT as default statistic
+    # We will change this to the user choice
+    if cls['statistic'].string() == 'DEFAULT':
+        pass
+    else:
+        for observation in onoff_obs:
+            observation.statistic(cls['statistic'].string())
+
+    # Return On/Off oberservation container
+    return onoff_obs
@@ -42,13 +42,19 @@ General parameters
 ``irf [string]``
     Instrumental response function.
+``(inexclusion = NONE) [file]``
+    Optional FITS file containing a WCS map in the first hdu that defines sky
+    regions not to be used for background estimation in ``ONOFF`` analysis (where
+    map value != 0).
+
 ``(edisp = no) [boolean]``
-    Applies energy dispersion to response computation.
+    Applies energy dispersion to response computation (for ``3D`` analysis only,
+    energy dispersion is always taken into account in ``ONOFF`` analysis).
 ``outfile [file]``
     Name of the light curve output file.
-``(tbinalg = GTI) <FILE|LIN|GTI> [string]``
+``tbinalg <FILE|LIN|GTI> [string]``
     Algorithm for defining time bins.
 ``tmin [time]``
@@ -63,6 +69,10 @@ General parameters
 ``tbinfile [file]``
     File defining the time binning.
+``method  <3D|ONOFF> [string]``
+    Selects between ``3D`` analysis (3D spatial/energy likelihood) and ``ONOFF``
+    analysis (1D likelihood with background from Off regions).
+
 ``emin [real]``
     Lower energy limit of events (in TeV).
@@ -70,7 +80,7 @@ General parameters
     Upper energy limit of events (in TeV).
 ``enumbins [integer]``
-    Number of energy bins per light curve bin (0=unbinned).
+    Number of energy bins per light curve bin (0=unbinned for ``3D`` analysis only).
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
@@ -79,19 +89,52 @@ General parameters
     Projection method.
 ``xref [real]``
-    Right Ascension / Galactic longitude of image centre (J2000, in degrees).
+    Right Ascension / Galactic longitude of image centre for ``3D`` analysis or
+    source region centre for ``ONOFF`` analysis (J2000, in degrees).
 ``yref [real]``
-    Declination / Galactic latitude of image centre (J2000, in degrees).
+    Declination / Galactic latitude of image centre for ``3D`` analysis or
+    source region centre for ``ONOFF`` analysis (J2000, in degrees).
 ``nxpix [integer]``
-    Size of the Right Ascension / Galactic longitude axis (in pixels).
+    Size of the Right Ascension / Galactic longitude axis for ``3D`` analysis (in pixels).
 ``nypix [integer]``
-    Size of the Declination / Galactic latitude axis (in pixels).
+    Size of the Declination / Galactic latitude axis for ``3D`` analysis (in pixels).
 ``binsz [real]``
-    Pixel size (in degrees/pixel).
+    Pixel size for ``3D`` analysis (in degrees/pixel).
+
+``(srcshape = CIRCLE) [string]``
+    Shape of the source region for ``ONOFF`` analysis.
+    ``CIRCLE``: circular region around given position.
+
+``rad [real]``
+    Radius of source region circle for ``ONOFF`` analysis (deg)
+
+``(bkgmethod = REFLECTED) [string]``
+    Method for background estimation in ``ONOFF`` analysis.
+    ``REFLECTED:`` background evaluated in regions with the same shape as
+    source region reflected w.r.t. pointing direction for each observation.
+
+``(bkgregmin = 2) [integer]``
+    Minimum number of background regions that are required for an observation in
+    ``ONOFF`` analysis. If this number of background regions is not available the
+    observation is skipped.
+
+``(maxoffset = 4.0) [real]``
+    Maximum offset in degrees of source from camera center to accept the
+    observation for ``ONOFF`` analysis.
+
+``(etruemin = 0.01) [real]``
+    Minimum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV).
+
+``(etruemax = 0.01) [real]``
+    Maximum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV).
+
+``(etruebins = 30) [integer]``
+    Number of bins per decade for true energy bins to evaluate instrumental
+    response in ``ONOFF`` analysis.
 ``(statistic = DEFAULT) <DEFAULT|CSTAT|WSTAT|CHI2> [string]``
     Optimization statistic. ``DEFAULT`` uses the default statistic for all
@@ -10,14 +10,16 @@ Synopsis
 --------
 This script computes spectra by performing a maximum likelihood fit using
-:doc:`ctlike` in a series of phase bins for pulsars. The phase bins can be
+:doc:`ctlike` in a series of phase bins for periodic sources. The input event
+list must have a ``PHASE`` column containing the phase values, that can be
+assigned using :doc:`ctphase`. The phase bins can be
 either specified in an ASCII file or as an interval divided into equally sized
 phase bins. The format of the ASCII file is one row per phase bin, each
 specifying the start of stop value of the phase bin, separated by a whitespace.
 The phase goes from 0.0 to 1.0.
-On output the script writes the fitting results into a FITS file. The script
-also produces one XML file per phase bin that contains the fitting results for
+The script writes the fit results into a FITS file. The script
+also produces one XML file per phase bin that contains the best-fit model for
 that bin.
@@ -39,8 +41,14 @@ General parameters
 ``irf [string]``
     Instrumental response function.
+``(inexclusion = NONE) [file]``
+    Optional FITS file containing a WCS map in the first hdu that defines sky
+    regions not to be used for background estimation in ``ONOFF`` analysis (where
+    map value != 0).
+
 ``(edisp = no) [boolean]``
-    Applies energy dispersion to response computation.
+    Applies energy dispersion to response computation (for ``3D`` analysis only,
+    energy dispersion is always taken into account in ``ONOFF`` analysis).
 ``outfile [file]``
     Name of the XML output file. The phase interval will be automatically
@@ -55,6 +63,10 @@ General parameters
 ``phbinfile [file]``
     File defining the phase binning.
+``method <3D|ONOFF> [string]``
+    Selects between 3D analysis (3D spatial/energy likelihood) and ONOFF
+    analysis (1D likelihood with background from Off regions).
+
 ``emin [real]``
     Lower energy limit of events (in TeV).
@@ -71,19 +83,52 @@ General parameters
     Projection method.
 ``xref [real]``
-    Right Ascension / Galactic longitude of image centre (J2000, in degrees).
+    Right Ascension / Galactic longitude of image centre for ``3D`` analysis or
+    source region centre for ``ONOFF`` analysis (J2000, in degrees).
 ``yref [real]``
-    Declination / Galactic latitude of image centre (J2000, in degrees).
+    Declination / Galactic latitude of image centre for ``3D`` analysis or
+    source region centre for ``ONOFF`` analysis (J2000, in degrees).
 ``nxpix [integer]``
-    Size of the Right Ascension / Galactic longitude axis (in pixels).
+    Size of the Right Ascension / Galactic longitude axis for ``3D`` analysis (in pixels).
 ``nypix [integer]``
-    Size of the Declination / Galactic latitude axis (in pixels).
+    Size of the Declination / Galactic latitude axis for ``3D`` analysis (in pixels).
 ``binsz [real]``
-    Pixel size (in degrees/pixel).
+    Pixel size for ``3D`` analysis (in degrees/pixel).
+
+``(srcshape = CIRCLE) [string]``
+    Shape of the source region for ``ONOFF`` analysis.
+    ``CIRCLE``: circular region around given position.
+
+``rad [real]``
+    Radius of source region circle for ``ONOFF`` analysis (deg)
+
+``(bkgmethod = REFLECTED) [string]``
+    Method for background estimation in ``ONOFF`` analysis.
+    ``REFLECTED:`` background evaluated in regions with the same shape as
+    source region reflected w.r.t. pointing direction for each observation.
+
+``(bkgregmin = 2) [integer]``
+    Minimum number of background regions that are required for an observation in
+    ``ONOFF`` analysis. If this number of background regions is not available the
+    observation is skipped.
+
+``(maxoffset = 4.0) [real]``
+    Maximum offset in degrees of source from camera center to accept the
+    observation for ``ONOFF`` analysis.
+
+``(etruemin = 0.01) [real]``
+    Minimum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV).
+
+``(etruemax = 0.01) [real]``
+    Maximum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV).
+
+``(etruebins = 30) [integer]``
+    Number of bins per decade for true energy bins to evaluate instrumental
+    response in ``ONOFF`` analysis.
 ``(statistic = DEFAULT) <DEFAULT|CSTAT|WSTAT|CHI2> [string]``
     Optimization statistic. ``DEFAULT`` uses the default statistic for all
@@ -129,4 +174,5 @@ Standard parameters
 Related tools or scripts
 ------------------------
+:doc:`ctphase`
 :doc:`ctlike`
@@ -34,12 +34,13 @@ The time interval is linearly divided into 20 time bins.
    Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml
    Source name [Crab] Src001
    Algorithm for defining time bins (FILE|LIN|GTI) [GTI] LIN
-   Lightcurve start time (MJD) [51544.5] 59235.5
-   Lightcurve stop time (MJD) [51544.6] 59236.0
+   Lightcurve start time (UTC string, JD, MJD or MET in seconds) [MJD 51544.5] MJD 59235.5
+   Lightcurve stop time (UTC string, JD, MJD or MET in seconds) [MJD 51544.6] MJD 59236.0
    Number of time bins (1-10000) [5] 20
-   Number of energy bins for binned (0=unbinned) (0-100) [0]
+   Analysis method (3D|ONOFF) [3D]
    Lower energy limit of events (TeV) [0.1]
    Upper energy limit of events (TeV) [100.0]
+   Number of energy bins for binned (0=unbinned for Cube analysis only) (0-100) [0]
    Output light curve file [lightcurve.fits]
 The resulting light curve is shown in the figure below.
@@ -100,6 +100,7 @@ public:
     %rename(_get_skydir)               get_skydir;
     %rename(_set_outfile_name)         set_outfile_name;
     %rename(_is_stacked)               is_stacked;
+    %rename(_is_onoff)                 is_onoff;
     %rename(_set_response)             set_response;
     %rename(_set_edisp)                set_edisp;
     %rename(_restore_edisp)            restore_edisp;
@@ -146,6 +147,7 @@ public:
     // Protected methods that query user parameters
     bool              is_stacked(void);
+    bool              is_onoff(void);
     // Protected support methods
     void              set_response(GObservations& obs);
@@ -1195,6 +1195,54 @@ bool ctool::is_stacked(void)
 /***********************************************************************//**
+ * @brief Query user parameters for On/Off analysis
+ *
+ * @return True if an On/Off analysis should be done.
+ *
+ * Queries the parameter "method", and if this is "ONOFF" signal that an On/
+ * Off analysis is requested. In that case, all necessary input parameters of
+ * csphagen are queried: "inexclusion", "emin", "emax", "enumbins",
+ * "coordsys", "xref", "yref", "srcshape", "rad", "bkgmethod", "bkgregmin",
+ * "maxoffset", "etruemin", "etruemax", and "etruebins".
+ *
+ * Using this method assures that the parameters are always queried in the
+ * same order.
+ ***************************************************************************/
+bool ctool::is_onoff(void)
+{
+    // Set onoff flag
+    bool onoff = ((*this)["method"].string() == "ONOFF");
+    
+    // Query analysis method chosen
+    if (onoff) {
+    
+        // Query csphagen parameters
+        (*this)["inexclusion"].filename();
+        (*this)["emin"].real();
+        (*this)["emax"].real();
+        (*this)["enumbins"].integer();
+        (*this)["coordsys"].string();
+        (*this)["xref"].real();
+        (*this)["yref"].real();
+        if ((*this)["srcshape"].string() == "CIRCLE"){
+            (*this)["rad"].real();
+        }
+        if ((*this)["bkgmethod"].string() == "REFLECTED"){
+            (*this)["bkgregmin"].integer();
+        }
+        (*this)["maxoffset"].real();
+        (*this)["etruemin"].real();
+        (*this)["etruemax"].real();
+        (*this)["etruebins"].integer();
+
+    } // endif: method was ONOFF
+
+    // Return onoff flag
+    return onoff;
+}
+
+
+/***********************************************************************//**
  * @brief Set response for all CTA observations in container
  *
  * @param[in,out] obs Observation container
@@ -104,6 +104,7 @@ protected:
     // Protected methods that query user parameters
     bool              is_stacked(void);
+    bool              is_onoff(void);
     // Protected methods for logging
     void              log_observations(const GChatter&      chatter,
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="PointSource">
+    <spectrum type="PowerLaw">
+       <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+       <parameter name="Index"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+       <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
+      <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="0"/>
+    </spatialModel>
+  </source>
+  <source name="Background" type="RadialAcceptance" instrument="CTAOnOff">
+    <spectrum type="PowerLaw">
+       <parameter name="Prefactor"   scale="1e-6" value="61.8" min="0.0"  max="1000.0" free="1"/>
+       <parameter name="Index"       scale="-1"   value="1.85" min="0.0"  max="+5.0"   free="1"/>
+       <parameter name="PivotEnergy" scale="1e6"  value="1.0"  min="0.01" max="1000.0" free="0"/>
+    </spectrum>
+    <radialModel type="Gaussian">
+       <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
+    </radialModel>
+  </source>
+</source_library>
@@ -455,6 +455,9 @@ class Test(test):
         # Test is_stacked() method
         # TODO: implement
+        # Test is_onoff() method
+        # TODO: implement
+
         # Test set_response() method
         # TODO: implement
@@ -42,6 +42,10 @@ class Test(test):
         # Call base class constructor
         test.__init__(self)
+        # Add off-axis events for classical analysis
+        self._offaxis_events = self._datadir + '/crab_offaxis1.fits'
+        self._model_onoff    = self._datadir + '/crab_onoff.xml'
+
         # Return
         return
@@ -73,7 +77,7 @@ class Test(test):
                          ' inmodel="'+self._model+'" srcname="Crab"'+ \
                          ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                          ' tbinalg="LIN" tmin="MJD 51544.50" tmax="MJD 51544.53"'+ \
-                         ' tbins=3 enumbins=0 emin=0.1 emax=100.0'+ \
+                         ' tbins=3 method="3D" enumbins=0 emin=0.1 emax=100.0'+ \
                          ' outfile="lightcurve_cmd1.fits"'+ \
                          ' logfile="cslightcrv_cmd1.log" chatter=1'
@@ -86,14 +90,14 @@ class Test(test):
              'Check successful execution from command line')
         # Check light curve
-        self._check_light_curve('lightcurve_cmd1.fits', 3)
+        self._check_light_curve('lightcurve_cmd1.fits', 3, 1.6e-16)
         # Setup cslightcrv command
         cmd = cslightcrv+' inobs="events_that_do_not_exist.fits"'+ \
                          ' inmodel="'+self._model+'" srcname="Crab"'+ \
                          ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                          ' tbinalg="LIN" tmin="MJD 51544.50" tmax="MJD 51544.53"'+ \
-                         ' tbins=3 enumbins=0 emin=0.1 emax=100.0'+ \
+                         ' tbins=3 method="3D" enumbins=0 emin=0.1 emax=100.0'+ \
                          ' outfile="lightcurve_cmd1.fits"'+ \
                          ' logfile="cslightcrv_cmd2.log" debug=yes'
@@ -123,6 +127,7 @@ class Test(test):
         lcrv['tmin']     = 'MJD 51544.50'
         lcrv['tmax']     = 'MJD 51544.53'
         lcrv['tbins']    = 3
+        lcrv['method']   = '3D'
         lcrv['enumbins'] = 0
         lcrv['emin']     = 0.1
         lcrv['emax']     = 100.0
@@ -137,7 +142,7 @@ class Test(test):
         lcrv.save()
         # Check light curve
-        self._check_light_curve('cslightcrv_py1.fits', 3)
+        self._check_light_curve('cslightcrv_py1.fits', 3, 1.6e-16)
         # Now use FILE as time bin algorithm. For this we need first to
         # create an ASCII file. We use now 6 time bins. The ASCII file
@@ -159,6 +164,7 @@ class Test(test):
         lcrv['irf']      = self._irf
         lcrv['tbinalg']  = 'FILE'
         lcrv['tbinfile'] = 'cslightcrv_py2.dat'
+        lcrv['method']   = '3D'
         lcrv['enumbins'] = 0
         lcrv['emin']     = 0.1
         lcrv['emax']     = 100.0
@@ -171,7 +177,7 @@ class Test(test):
         lcrv.execute()
         # Check light curve
-        self._check_light_curve('cslightcrv_py2.fits', 2)
+        self._check_light_curve('cslightcrv_py2.fits', 2, 1.6e-16)
         # Now we setup an observation container on input. We attached the
         # model to the observation container so that cslightcrv should
@@ -188,6 +194,7 @@ class Test(test):
         lcrv['caldb']    = self._caldb
         lcrv['irf']      = self._irf
         lcrv['tbinalg']  = 'GTI'
+        lcrv['method']   = '3D'
         lcrv['enumbins'] = 0
         lcrv['emin']     = 0.1
         lcrv['emax']     = 100.0
@@ -199,9 +206,9 @@ class Test(test):
         lcrv.execute()
         # Check light curve
-        self._check_light_curve('cslightcrv_py3.fits', 1)
+        self._check_light_curve('cslightcrv_py3.fits', 1, 1.6e-16)
-        # Finally we set-up a binned cslightcrv
+        # Binned cslightcrv
         lcrv = cscripts.cslightcrv()
         lcrv['inobs']    = self._events
         lcrv['inmodel']  = self._model
@@ -212,6 +219,7 @@ class Test(test):
         lcrv['tmin']     = 'MJD 51544.50'
         lcrv['tmax']     = 'MJD 51544.53'
         lcrv['tbins']    = 2
+        lcrv['method']   = '3D'
         lcrv['emin']     = 0.1
         lcrv['emax']     = 100.0
         lcrv['enumbins'] = 10
@@ -230,13 +238,43 @@ class Test(test):
         lcrv.execute()
         # Check light curve
-        self._check_light_curve('cslightcrv_py4.fits', 2)
+        self._check_light_curve('cslightcrv_py4.fits', 2, 1.6e-16)
+
+        # cslightcrv with classical analysis
+        lcrv = cscripts.cslightcrv()
+        lcrv['inobs']     = self._offaxis_events
+        lcrv['inmodel']   = self._model_onoff
+        lcrv['srcname']   = 'Crab'
+        lcrv['caldb']     = self._caldb
+        lcrv['irf']       = self._irf
+        lcrv['tbinalg']   = 'LIN'
+        lcrv['tmin']      = 'MJD 58849.0008'
+        lcrv['tmax']      = 'MJD 58849.0112'
+        lcrv['tbins']     = 2
+        lcrv['method']    = 'ONOFF'
+        lcrv['emin']      = 0.1
+        lcrv['emax']      = 100.0
+        lcrv['enumbins']  = 10
+        lcrv['coordsys']  = 'CEL'
+        lcrv['xref']      = 83.63
+        lcrv['yref']      = 22.01
+        lcrv['rad']       = 0.2
+        lcrv['statistic'] = 'WSTAT'
+        lcrv['outfile']   = 'cslightcrv_py5.fits'
+        lcrv['logfile']   = 'cslightcrv_py5.log'
+        lcrv['chatter']   = 4
+
+        # Execute cslightcrv script
+        lcrv.execute()
+
+        # Check light curve
+        self._check_light_curve('cslightcrv_py5.fits', 2, 6.e-16)
         # Return
         return
     # Check light curve result file
-    def _check_light_curve(self, filename, bins):
+    def _check_light_curve(self, filename, bins,pref):
         """
         Check light curve file
         """
@@ -266,6 +304,13 @@ class Test(test):
             self.test_assert(table.contains(col),
                  'FITS file contains "'+col+'" column')
+        # Check that table has been filled
+        # Prefactor has right order of magnitude
+        for s in range(table.nrows()):
+            self.test_value(table['Prefactor'][s], pref, 0.2*pref,
+                            'Check prefactor value')
+
+
         # Close FITS file
         fits.close()
@@ -44,6 +44,8 @@ class Test(test):
         # Set test data
         self._phased_events  = self._datadir + '/phased_events.fits'
+        self._offaxis_events = self._datadir + '/crab_offaxis1.fits'
+        self._model_onoff = self._datadir + '/crab_onoff.xml'
         # Return
         return
@@ -75,7 +77,7 @@ class Test(test):
         cmd = csphasecrv+' inobs="'+self._phased_events+'"'+ \
                          ' inmodel="'+self._model+'" srcname="Crab"'+ \
                          ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
-                         ' phbinalg="LIN" phbins=5'+ \
+                         ' phbinalg="LIN" phbins=5 method="3D"'+ \
                          ' enumbins=0 emin=0.1 emax=100.0'+ \
                          ' outfile="csphasecrv_cmd1.fits"'+ \
                          ' logfile="csphasecrv_cmd1.log" chatter=1'
@@ -95,7 +97,7 @@ class Test(test):
         cmd = csphasecrv+' inobs="events_that_do_not_exist.fits"'+ \
                          ' inmodel="'+self._model+'" srcname="Crab"'+ \
                          ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
-                         ' phbinalg="LIN" phbins=5'+ \
+                         ' phbinalg="LIN" phbins=5 method="3D"'+ \
                          ' enumbins=0 emin=0.1 emax=100.0'+ \
                          ' outfile="csphasecrv_cmd2.fits"'+ \
                          ' logfile="csphasecrv_cmd2.log" debug=yes chatter=1'
@@ -124,6 +126,7 @@ class Test(test):
         pcrv['irf']      = self._irf
         pcrv['phbinalg'] = 'LIN'
         pcrv['phbins']   = 5
+        pcrv['method']   = '3D'
         pcrv['enumbins'] = 0
         pcrv['emin']     = 0.1
         pcrv['emax']     = 100.0
@@ -160,6 +163,7 @@ class Test(test):
         pcrv['irf']       = self._irf
         pcrv['phbinalg']  = 'FILE'
         pcrv['phbinfile'] = 'csphasecrv_py2.dat'
+        pcrv['method']    = '3D'
         pcrv['enumbins']  = 0
         pcrv['emin']      = 0.1
         pcrv['emax']      = 100.0
@@ -190,6 +194,7 @@ class Test(test):
         pcrv['irf']      = self._irf
         pcrv['phbinalg'] = 'LIN'
         pcrv['phbins']   = 5
+        pcrv['method']   = '3D'
         pcrv['enumbins'] = 0
         pcrv['emin']     = 0.1
         pcrv['emax']     = 100.0
@@ -204,7 +209,7 @@ class Test(test):
         # Check phase curve
         self._check_phase_curve('csphasecrv_py3.fits', 5)
-        # Finally we set-up a binned csphasecrv
+        # Binned csphasecrv
         pcrv = cscripts.csphasecrv()
         pcrv['inobs']    = self._phased_events
         pcrv['inmodel']  = self._model
@@ -213,6 +218,7 @@ class Test(test):
         pcrv['irf']      = self._irf
         pcrv['phbinalg'] = 'LIN'
         pcrv['phbins']   = 2
+        pcrv['method']   = '3D'
         pcrv['emin']     = 0.1
         pcrv['emax']     = 100.0
         pcrv['enumbins'] = 10
@@ -234,6 +240,35 @@ class Test(test):
         # Check phase curve
         self._check_phase_curve('csphasecrv_py4.fits', 2)
+        # On/Off csphasecrv
+        pcrv = cscripts.csphasecrv()
+        pcrv['inobs']     = self._offaxis_events
+        pcrv['inmodel']   = self._model_onoff
+        pcrv['srcname']   = 'Crab'
+        pcrv['caldb']     = self._caldb
+        pcrv['irf']       = self._irf
+        pcrv['phbinalg']  = 'LIN'
+        pcrv['phbins']    = 2
+        pcrv['method']    = 'ONOFF'
+        pcrv['statistic'] = 'WSTAT'
+        pcrv['emin']      = 0.1
+        pcrv['emax']      = 100.0
+        pcrv['enumbins']  = 10
+        pcrv['coordsys']  = 'CEL'
+        pcrv['xref']      = 83.63
+        pcrv['yref']      = 22.01
+        pcrv['rad']       = 0.2
+        pcrv['outfile']   = 'csphasecrv_py5.fits'
+        pcrv['logfile']   = 'csphasecrv_py5.log'
+        pcrv['chatter']   = 4
+
+        # Execute csphasecrv script
+        pcrv.logFileOpen()   # Make sure we get a log file
+        pcrv.execute()
+
+        # Check phase curve
+        self._check_phase_curve('csphasecrv_py5.fits', 2)
+
         # Return
         return
@@ -267,6 +302,12 @@ class Test(test):
             self.test_assert(table.contains(col),
                  'FITS file contains "'+col+'" column')
+        # Check that table has been filled
+        # Prefactor has right order of magnitude
+        for s in range(table.nrows()):
+            self.test_value(table['Prefactor'][s], 4.e-16, 3.e-16,
+                            'Check prefactor value')
+
         # Close FITS file
         fits.close()