Revise cstsdist parameters (#1567)

The cstsdist parameter interface is now much simpler, requiring on input either an unbinned, binned or stacked observation. If an observation definition XML file is provided it is assumed that the file only contains event lists.

Revise cstsdist parameters (#1567)
The cstsdist parameter interface is now much simpler, requiring on input either an unbinned, binned or stacked observation. If an observation definition XML file is provided it is assumed that the file only contains event lists.
Jürgen Knödlseder
1 parent b030e64d
Showing 7 changed files with 161 additions and 116 deletions
ChangeLog
NEWS
cscripts/cstsdist.par
cscripts/cstsdist.py
doc/source/users/reference_manual/cstsdist.rst
src/Makefile.am
test/test_cstsdist.py
@@ -26,6 +26,7 @@
         Add "NODES" method to csspec to support non-CTA analysis (#1927)
         Use t-type parameters for all "tmin" and "tmax" parameters (#1864)
         Set OBS_ID in ctobssim (#1836)
+        Revise cstsdist parameters (#1567)
         Add ctbutterfly::butterfly() method (#1512)
@@ -319,8 +319,14 @@ The script now derives from the csobservation base class. A &quot;statistic&quot;
 parameter was added to allow setting the optimizer statistic for the Test
 Statistic distribution generation (#2237).
-Use t-type parameters for "tmin" and "tmax" to all specifying times in various
-formats (#1864).
+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).
 cstsmapmerge - Merge Test Statistic maps
@@ -16,7 +16,7 @@
 #
 # Input and output parameters
 #============================
-inobs,     f, h, NONE,,, "Input event list, counts cube, or observation definition XML file"
+inobs,     f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file"
 inmodel,   f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
 srcname,   s, a, Crab,,, "Source name"
 expcube,   f, a, NONE,,, "Input exposure cube file"
@@ -33,18 +33,6 @@ outfile,   f, a, ts.dat,,, &quot;Output Test Statistic distribution file&quot;
 # Script parameters
 #==================
 ntrials,   i, a, 10,,, "Number of trials"
-ra,        r, a, 83.63,0,360, "RA of pointing (degrees)"
-dec,       r, a, 22.01,-90,90, "Dec of pointing (degrees)"
-rad,       r, h, 5.0,0,180, "Radius of ROI (degrees)"
-coordsys,  s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
-proj,      s, a, TAN,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
-emin,      r, a, 0.1,,, "Lower energy limit (TeV)"
-emax,      r, a, 100.0,,, "Upper energy limit (TeV)"
-enumbins,  i, a, 0,,, "Number of energy bins (0=unbinned)"
-tmin,      t, a, 2020-01-01T00:00:00,,, "Start time (UTC string, JD, MJD or MET in seconds)"
-tmax,      t, a, 2020-01-01T00:30:00,,, "Stop time (UTC string, JD, MJD or MET in seconds)"
-npix,      i, a, 200,,, "Number of pixels for binned"
-binsz,     r, a, 0.05,,, "Pixel size for binned (degrees/pixel)"
 statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
 #
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates the TS distribution for a particular model
 #
-# Copyright (C) 2011-2017 Juergen Knoedlseder
+# Copyright (C) 2011-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
@@ -37,7 +37,7 @@ class cstsdist(ctools.csobservation):
     # Constructor
     def __init__(self, *argv):
         """
-        Constructor.
+        Constructor
         """
         # Initialise application by calling the appropriate class constructor
         self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
@@ -45,6 +45,7 @@ class cstsdist(ctools.csobservation):
         # Initialise some members
         self._srcname     = ''
         self._log_clients = False
+        self._model       = None
         # Return
         return
@@ -65,13 +66,6 @@ class cstsdist(ctools.csobservation):
         # Get source name
         self._srcname = self['srcname'].string()
-        # Query parameters for binned if requested
-        if self['enumbins'].integer() != 0:
-            self['npix'].integer()
-            self['binsz'].real()
-            self['coordsys'].string()
-            self['proj'].string()
-
         # Set models if we have none
         if self.obs().models().size() == 0:
             self.obs().models(self['inmodel'].filename())
@@ -79,15 +73,73 @@ class cstsdist(ctools.csobservation):
         # Query parameters
         self['edisp'].boolean()
         self['ntrials'].integer()
-        self['outfile'].filename()
+        
         self['debug'].boolean()
+        # Read ahead output parameters
+        #if self._read_ahead():
+        #    self['outfile'].filename()
+        self['outfile'].filename()
+
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
         # Return
         return
+    def _sim(self, seed):
+        """
+        Return a simulated observation container
+
+        Parameters
+        ----------
+        seed : int
+            Random number generator seed
+
+        Returns
+        -------
+        sim : `~gammalib.GObservations`
+            Simulated observation container
+        """
+        # If observation is a counts cube then simulate events from the counts
+        # cube model ...
+        if self.obs().size() == 1 and self.obs()[0].eventtype() == 'CountsCube':
+
+            # If no counts cube model exists then compute it now
+            if self._model == None:
+                model            = ctools.ctmodel(self.obs())
+                model['debug']   = self['debug'].boolean()
+                model['chatter'] = self['chatter'].integer()
+                model.run()
+                self._model = model.cube().copy() # Save copy for persistence
+
+            # Allocate random number generator
+            ran = gammalib.GRan()
+
+            # Get copy of model map
+            counts = self._model.counts().copy()
+
+            # Randomize counts
+            for i in range(counts.npix()):
+                counts[i] = ran.poisson(counts[i])
+
+            # Copy observations
+            sim = self.obs().copy()
+
+            # Set counts map
+            sim[0].events().counts(counts)
+
+        # ... otherwise simuate events from the observation container (works
+        # only for event lists
+        else:
+            sim = obsutils.sim(self.obs(),
+                               seed  = seed,
+                               log   = self._log_clients,
+                               debug = self['debug'].boolean())
+
+        # Return simulated observation
+        return sim
+
     def _trial(self, seed):
         """
         Create the TS for a single trial
@@ -105,29 +157,8 @@ class cstsdist(ctools.csobservation):
         # Write header
         self._log_header2(gammalib.EXPLICIT, 'Trial %d' % (seed+1))
-        # Set default binned parameters
-        coordsys = 'CEL'
-        proj     = 'TAN'
-        npix     = 0
-        binsz    = 0.0
-
-        # If binned analysis is requested the read the binned parameters
-        if self['enumbins'].integer() != 0:
-            npix     = self['npix'].integer()
-            binsz    = self['binsz'].real()
-            coordsys = self['coordsys'].string()
-            proj     = self['proj'].string()
-
         # Simulate events
-        sim = obsutils.sim(self.obs(),
-                           nbins = self['enumbins'].integer(),
-                           seed  = seed,
-                           proj  = proj,
-                           coord = coordsys,
-                           binsz = binsz,
-                           npix  = npix,
-                           log   = self._log_clients,
-                           debug = self['debug'].boolean())
+        sim = self._sim(seed)
         # Determine number of events in simulation
         nevents = 0.0
@@ -160,7 +191,7 @@ class cstsdist(ctools.csobservation):
             for model in models:
                 self._log_value(gammalib.EXPLICIT, 'Model', model.name())
                 for par in model:
-                    self._log_string(gammalib.EXPLICIT, str(par)+'\n')
+                    self._log_string(gammalib.EXPLICIT, str(par))
         elif self._logNormal():
             name  = 'Trial %d' % seed
             value = 'TS=%.3f  Prefactor=%e +/- %e' % \
@@ -15,9 +15,6 @@ sets. The Test Statistics is defined as twice the log-likelihood difference
 that is obtained when fitting simulated data with and without a given
 source model component.
-This script either performs an unbinned or a binned simulation, stacked 
-analysis is not yet supported.
-
 cstsdist will create an ASCII file in comma-separated value (CSV) format,
 containing one row per TS computation. The first row is a header row providing
 the column names. The following rows give the TS value, the log-likelihood 
@@ -26,13 +23,18 @@ fitted events, as well as the values and errors for all fitted parameters.
 From the output file, TS distribution plots can be generated using for
 example the ``show_ts_distribution.py`` script in the examples folder. The
-script require matplotlib for plotting.
+script requires matplotlib for plotting.
+
+.. warning::
+   This script does not work for On/Off observations. If an observation
+   definition XML file is specified the script assumes that all observations
+   are event lists.
 General parameters
 ------------------
-``(inobs = NONE) [file]``
+``inobs [file]``
     Input event list, counts cube or observation definition XML file.
 ``inmodel [file]``
@@ -72,42 +74,6 @@ General parameters
 ``ntrials [integer]``
     Number of Monte Carlo samples.
-``ra [real]``
-    Right Ascension of CTA pointing (J2000, in degrees).
-
-``dec [real]``
-    Declination of CTA pointing (J2000, in degrees).
-
-``(rad = 5.0) [real]``
-    ROI radius (in degrees).
-
-``coordsys <CEL|GAL> [string]``
-    Coordinate system (CEL - celestial, GAL - galactic).
-
-``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
-    Projection method.
-
-``emin [real]``
-    Lower energy limit (in TeV).
-
-``emax [real]``
-    Upper energy limit (in TeV).
-
-``enumbins [integer]``
-    Number of energy bins (0=unbinned).
-
-``tmin [time]``
-    Start time (UTC string, JD, MJD or MET in seconds).
-
-``tmax [time]``
-    Stop time (UTC string, JD, MJD or MET in seconds).
-
-``npix [integer]``
-    Number of pixels for binned analysis.
-
-``binsz [real]``
-    Pixel size for binned analysis.
-
 ``(statistic = DEFAULT) <DEFAULT|CSTAT|WSTAT|CHI2> [string]``
     Optimization statistic. ``DEFAULT`` uses the default statistic for all
     observations, which is ``CSTAT`` or the statistic specified in the
@@ -100,7 +100,7 @@ help-install:
         echo "$(INSTALL_DATA) $$file $(DESTDIR)$(datadir)/help/$$helpname" ; \
         $(INSTALL_DATA) "$$file" $(DESTDIR)$(datadir)/help/"$$helpname" ; \
         sed -i.bak -e 's/``//g' -e 's/`//g' -e 's/:ref://g' -e 's/:doc://g' \
-            -e 's/.. note::/Note/g' -e '1d' \
+            -e 's/.. note::/Note/g' -e 's/.. warning::/Warning/g' -e '1d' \
             $(DESTDIR)$(datadir)/help/"$$helpname" ; \
         rm -rf $(DESTDIR)$(datadir)/help/"$$helpname".bak ; \
     done
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the cstsdist script.
 #
-# 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
@@ -42,6 +42,10 @@ class Test(test):
         # Call base class constructor
         test.__init__(self)
+        # Set members
+        self._stacked_inobs = self._datadir + '/obs_stacked.xml'
+        self._stacked_model = self._datadir + '/crab_bkgcube.xml'
+
         # Return
         return
@@ -55,7 +59,12 @@ class Test(test):
         # Append tests
         self.append(self._test_cmd, 'Test cstsdist on command line')
-        self.append(self._test_python, 'Test cstsdist from Python')
+        self.append(self._test_python_unbinned,
+                    'Test cstsdist from Python in unbinned mode')
+        self.append(self._test_python_binned,
+                    'Test cstsdist from Python in binned mode')
+        self.append(self._test_python_stacked,
+                    'Test cstsdist from Python in stacked mode')
         # Return
         return
@@ -69,12 +78,10 @@ class Test(test):
         cstsdist = self._script('cstsdist')
         # Setup cstsdist command
-        cmd = cstsdist+' inmodel="'+self._model+'"'+ \
-                       ' srcname="Crab"'+ \
+        cmd = cstsdist+' inobs="'+self._events+'"'+ \
+                       ' inmodel="'+self._model+'" srcname="Crab"'+ \
                        ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
-                       ' ntrials=1 ra=83.63 dec=22.01 emin=0.1 emax=100.0'+ \
-                       ' enumbins=0 tmin=0.0 tmax=1800.0 rad=5.0 npix=200 npix=200'+ \
-                       ' binsz=0.05'+ \
+                       ' ntrials=2 '+ \
                        ' outfile="cstsdist_cmd1.dat"'+ \
                        ' logfile="cstsdist_cmd1.log" chatter=2'
@@ -90,14 +97,12 @@ class Test(test):
         self._check_result_file('cstsdist_cmd1.dat')
         # Setup cstsdist command
-        cmd = cstsdist+' inmodel="model_that_does_not_exist.xml"'+ \
-                       ' srcname="Crab"'+ \
+        cmd = cstsdist+' inobs="event_file_that_does_not_exist.xml"'+ \
+                       ' inmodel="'+self._model+'" srcname="Crab"'+ \
                        ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
-                       ' ntrials=1 ra=83.63 dec=22.01 emin=0.1 emax=100.0'+ \
-                       ' enumbins=0 tmin=0.0 tmax=1800.0 rad=5.0 npix=200 npix=200'+ \
-                       ' binsz=0.05'+ \
+                       ' ntrials=1 '+ \
                        ' outfile="cstsdist_cmd2.dat"'+ \
-                       ' logfile="cstsdist_cmd2.log" debug=yes chatter=2'
+                       ' logfile="cstsdist_cmd2.log" chatter=2'
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
@@ -110,25 +115,18 @@ class Test(test):
         return
     # Test cstsdist from Python
-    def _test_python(self):
+    def _test_python_unbinned(self):
         """
         Test cstsdist from Python
         """
-        # Set-up cstsdist
+        # Set-up cstsdist for event list
         tsdist = cscripts.cstsdist()
+        tsdist['inobs']    = self._events
         tsdist['inmodel']  = self._model
         tsdist['srcname']  = 'Crab'
         tsdist['caldb']    = self._caldb
         tsdist['irf']      = self._irf
-        tsdist['ntrials']  = 1
-        tsdist['ra']       = 83.63
-        tsdist['dec']      = 22.01
-        tsdist['emin']     = 0.1
-        tsdist['emax']     = 100.0
-        tsdist['enumbins'] = 0
-        tsdist['tmin']     = 0.0
-        tsdist['tmax']     = 1800.0
-        tsdist['rad']      = 5.0
+        tsdist['ntrials']  = 2
         tsdist['outfile']  = 'cstsdist_py1.dat'
         tsdist['logfile']  = 'cstsdist_py1.log'
         tsdist['chatter']  = 3
@@ -144,8 +142,61 @@ class Test(test):
         # Return
         return
+    def _test_python_binned(self):
+        """
+        Test cstsdist from Python
+        """
+        # Set-up cstsdist for counts cube
+        tsdist = cscripts.cstsdist()
+        tsdist['inobs']    = self._cntcube
+        tsdist['inmodel']  = self._model
+        tsdist['srcname']  = 'Crab'
+        tsdist['expcube']  = 'NONE'
+        tsdist['psfcube']  = 'NONE'
+        tsdist['bkgcube']  = 'NONE'
+        tsdist['caldb']    = self._caldb
+        tsdist['irf']      = self._irf
+        tsdist['ntrials']  = 2
+        tsdist['outfile']  = 'cstsdist_py2.dat'
+        tsdist['logfile']  = 'cstsdist_py2.log'
+        tsdist['chatter']  = 4
+
+        # Run cstsdist script
+        tsdist.logFileOpen()   # Make sure we get a log file
+        tsdist.execute()
+
+        # Check pull distribution file
+        self._check_result_file('cstsdist_py2.dat')
+
+        # Return
+        return
+
+    def _test_python_stacked(self):
+        """
+        Test cstsdist from Python
+        """
+        # Set-up cstsdist for stacked cube
+        tsdist = cscripts.cstsdist()
+        tsdist['inobs']    = self._stacked_inobs
+        tsdist['inmodel']  = self._stacked_model
+        tsdist['srcname']  = 'Crab'
+        tsdist['ntrials']  = 2
+        tsdist['outfile']  = 'cstsdist_py3.dat'
+        tsdist['logfile']  = 'cstsdist_py3.log'
+        tsdist['chatter']  = 4
+
+        # Run cstsdist script
+        tsdist.logFileOpen()   # Make sure we get a log file
+        tsdist.execute()
+
+        # Check pull distribution file
+        self._check_result_file('cstsdist_py3.dat')
+
+        # Return
+        return
+
     # Check result file
-    def _check_result_file(self, filename, ncols=9):
+    def _check_result_file(self, filename, nrows=3, ncols=9):
         """
         Check result file
@@ -153,6 +204,8 @@ class Test(test):
         ----------
         filename : str
             Name of result file
+        nrows : int, optional
+            Required number of rows
         ncols : int, optional
             Required number of columns
         """
@@ -160,7 +213,7 @@ class Test(test):
         results = gammalib.GCsv(filename, ',')
         # Check dimensions
-        self.test_value(results.nrows(), 2,     'Check rows in TS file')
+        self.test_value(results.nrows(), nrows, 'Check rows in TS file')
         self.test_value(results.ncols(), ncols, 'Check columns in TS file')
         # Return