Handle On/Off observations in cspull (#2260)

Jürgen Knödlseder
1 parent b0341d4a
Showing 6 changed files with 184 additions and 124 deletions
ChangeLog
NEWS
cscripts/cspull.par
cscripts/cspull.py
cscripts/obsutils.py
doc/source/admin/release_history/1.5.rst
-2017-10-22
+2017-10-26
         * Version 1.5.0 released
         ========================
+        Handle On/Off observations by cspull (#2260)
         Add "statistic" parameter to all likelihood tools and scripts (#2237)
         Rename ctlike "stat" parameter to "statistic" (#2237)
         Enable handling of On/Off observations in csspec (#2231)
 New Features and Important Changes in ctools 1.5.0
-21 October 2017
+26 October 2017
 Introduction
@@ -227,6 +227,10 @@ curve generation (#2237).
 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 derived from the csobservation base class. A "statistic"
 parameter was added to allow setting the optimizer statistic for the pull
 distribution generation (#2237).
@@ -18,6 +18,8 @@
 #============================
 inobs,   f, h, NONE,,, "Input event list, counts cube, or observation definition XML file"
 inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
+onsrc,   s, a, ,,, "Name of On source (only for On/Off analysis)"
+onrad,   r, a, 0.2,,, "Radius of On region (deg)"
 caldb,   s, a, prod2,,, "Calibration database"
 irf,     s, a, South_0.5h,,, "Instrument response function"
 edisp,   b, h, no,,, "Apply energy dispersion?"
@@ -40,18 +40,6 @@ class cspull(ctools.csobservation):
         # Initialise application by calling the appropriate class constructor
         self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
-        # Initialise some members
-        self._edisp       = False
-        self._coordsys    = 'CEL'
-        self._proj        = 'TAN'
-        self._log_clients = False  # Static parameter
-        self._ntrials     = 0
-        self._npix        = 0
-        self._binsz       = 0.0
-        self._enumbins    = 0
-        self._seed        = 1
-        self._chatter     = 2
-
         # Return
         return
@@ -62,7 +50,7 @@ class cspull(ctools.csobservation):
         Get parameters from parfile
         """
         # If there are no observations in container then get some ...
-        if self.obs().size() == 0:
+        if self.obs().is_empty():
             self.obs(self._get_observations())
         # ... otherwise add response information and energy boundaries
@@ -73,25 +61,26 @@ class cspull(ctools.csobservation):
         # Set observation statistic
         self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
-        # Get number of energy bins
-        self._enumbins = self['enumbins'].integer()
+        # Query parameters for On/Off observation
+        if len(self['onsrc'].string()) > 0:
+            self['onrad'].real()
-        # Read parameters for binned if requested
-        if self._enumbins != 0:
-            self._npix     = self['npix'].integer()
-            self._binsz    = self['binsz'].real()
-            self._coordsys = self['coordsys'].string()
-            self._proj     = self['proj'].string()
+        # Query parameters for binned if requested
+        elif 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:
+        if self.obs().models().is_empty():
             self.obs().models(self['inmodel'].filename())
-        # Read other parameters    
-        self._ntrials = self['ntrials'].integer()
-        self._edisp   = self['edisp'].boolean()
-        self._seed    = self['seed'].integer()
-        self._chatter = self['chatter'].integer()
+        # Query other parameters
+        self['ntrials'].integer()
+        self['edisp'].boolean()
+        self['seed'].integer()
+        self['chatter'].integer()
         # Query some parameters
         self['outfile'].filename()
@@ -118,14 +107,19 @@ class cspull(ctools.csobservation):
             Summary string
         """
         # Extract information from observation
-        emin   = obs.events().ebounds().emin().TeV()
-        emax   = obs.events().ebounds().emax().TeV()
-        events = obs.events().number()
-        binned = (obs.events().classname() == 'GCTAEventCube')
-        if binned:
-            mode = 'binned'
+        if obs.classname() == 'GCTAOnOffObservation':
+            emin = obs.on_spec().ebounds().emin().TeV()
+            emax = obs.on_spec().ebounds().emin().TeV()
+            mode = 'On/Off'
         else:
-            mode = 'unbinned'
+            emin = obs.events().ebounds().emin().TeV()
+            emax = obs.events().ebounds().emax().TeV()
+            binned = (obs.events().classname() == 'GCTAEventCube')
+            if binned:
+                mode = 'binned'
+            else:
+                mode = 'unbinned'
+        events = obs.nobserved()
         # Compose summary string
         if events > 0:
@@ -151,38 +145,55 @@ class cspull(ctools.csobservation):
             Dictionary of results
         """
         # Write header
-        if self._logNormal():
-            self._log.header2('Trial '+str(seed-self._seed+1))
-
-        # If we have a binned obeservation then specify the lower and
-        # upper energy limit in TeV
-        if self._enumbins > 0:
-            emin = self['emin'].real()
-            emax = self['emax'].real()
+        self._log_header2(gammalib.NORMAL, 'Trial %d' %
+                          (seed-self['seed'].integer()+1))
+
+        # Get number of energy bins and On source name and initialise
+        # some parameters
+        nbins    = self['enumbins'].integer()
+        onsrc    = self['onsrc'].string()
+        emin     = None
+        emax     = None
+        binsz    = 0.0
+        npix     = 0
+        proj     = 'TAN'
+        coordsys = 'CEL'
+
+        # If we have a On source name then set On region radius
+        if len(onsrc) > 0:
+            onrad = self['onrad'].real()
+            emin  = self['emin'].real()
+            emax  = self['emax'].real()
         else:
-            emin = None
-            emax = None
+
+            # Reset On region source name and radius
+            onrad = 0.0
+            onsrc = None
+
+            # If we have a binned obeservation then specify the lower and
+            # upper energy limit in TeV
+            if nbins > 0:
+                emin     = self['emin'].real()
+                emax     = self['emax'].real()
+                binsz    = self['binsz'].real()
+                npix     = self['npix'].integer()
+                proj     = self['proj'].string()
+                coordsys = self['coordsys'].string()
         # Simulate events
         obs = obsutils.sim(self.obs(),
-                           emin=emin,
-                           emax=emax,
-                           nbins=self._enumbins,
-                           addbounds=True,
-                           seed=seed,
-                           binsz=self._binsz,
-                           npix=self._npix,
-                           proj=self._proj,
-                           coord=self._coordsys,
-                           edisp=self._edisp,
-                           log=self._log_clients,
-                           debug=self._logDebug(),
-                           chatter=self._chatter)
+                           emin=emin, emax=emax, nbins=nbins,
+                           onsrc=onsrc, onrad=onrad,
+                           addbounds=True, seed=seed,
+                           binsz=binsz, npix=npix, proj=proj, coord=coordsys,
+                           edisp=self['edisp'].boolean(),
+                           log=False, debug=self._logDebug(),
+                           chatter=self['chatter'].integer())
         # Determine number of events in simulation
         nevents = 0.0
         for run in obs:
-            nevents += run.events().number()
+            nevents += run.nobserved()
         # Write simulation results
         self._log_header3(gammalib.NORMAL, 'Simulation')
@@ -200,15 +211,15 @@ class cspull(ctools.csobservation):
             for model in models:
                 like = ctools.cterror(obs)
                 like['srcname'] = model.name()
-                like['edisp']   = self._edisp
+                like['edisp']   = self['edisp'].boolean()
                 like['debug']   = self._logDebug()
-                like['chatter'] = self._chatter
+                like['chatter'] = self['chatter'].integer()
                 like.run()
         else:
             like = ctools.ctlike(obs)
-            like['edisp']   = self._edisp
+            like['edisp']   = self['edisp'].boolean()
             like['debug']   = self._logDebug()
-            like['chatter'] = self._chatter
+            like['chatter'] = self['chatter'].integer()
             like.run()
         # Store results
@@ -217,8 +228,7 @@ class cspull(ctools.csobservation):
         models = like.obs().models()
         # Write result header
-        if self._logNormal():
-            self._log.header3('Pulls')
+        self._log_header3(gammalib.NORMAL, 'Pulls')
         # Gather results in form of a list of result columns and a
         # dictionary containing the results. The result contains the
@@ -294,10 +304,10 @@ class cspull(ctools.csobservation):
         self._log_header1(gammalib.TERSE, 'Generate pull distribution')
         # Loop over trials
-        for seed in range(self._ntrials):
+        for seed in range(self['ntrials'].integer()):
             # Make a trial and add initial seed
-            result = self._trial(seed + self._seed)
+            result = self._trial(seed + self['seed'].integer())
             # Write out trial result
             ioutils.write_csv_row(self['outfile'].filename().url(), seed,
@@ -20,13 +20,14 @@
 import math
 import gammalib
 import ctools
+import cscripts
 # ===================== #
 # Simulate observations #
 # ===================== #
 def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
-        emin=None, emax=None, nbins=0, addbounds=False,
+        emin=None, emax=None, nbins=0, onsrc=None, onrad=0.2, addbounds=False,
         binsz=0.05, npix=200, proj='TAN', coord='GAL'):
     """
     Simulate events for all observations in the container
@@ -59,6 +60,10 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
         Maximum energy of counts cube for binned (TeV)
     nbins : int, optional
         Number of energy bins (0=unbinned)
+    onsrc : str, optional
+        Name of source for On region (None if no On/Off obs. is used)
+    onrad : float, optional
+        Radius for On region (deg)
     addbounds : bool, optional
         Add boundaries at observation energies
     binsz : float, optional
@@ -105,61 +110,97 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
                 emin = min(run.events().ebounds().emin().TeV(), emin)
                 emax = max(run.events().ebounds().emax().TeV(), emax)
-        # Allocate ctbin application and set parameters
-        binning = ctools.ctbin(obssim.obs())
-        binning['ebinalg']  = 'LOG'
-        binning['emin']     = emin
-        binning['emax']     = emax
-        binning['enumbins'] = nbins
-        binning['usepnt']   = True # Use pointing for map centre
-        binning['nxpix']    = npix
-        binning['nypix']    = npix
-        binning['binsz']    = binsz
-        binning['coordsys'] = coord
-        binning['proj']     = proj
-        binning['chatter']  = chatter
-        binning['debug']    = debug
-
-        # Optionally open the log file
-        if log:
-            binning.logFileOpen()
-
-        # Run ctbin application. This will loop over all observations in
-        # the container and bin the events in counts maps
-        binning.run()
-
-        # If we have multiple input observations then create stacked response
-        # cubes and append them to the observation
-        if len(obssim.obs()) > 1:
-
-            # Get stacked response (use pointing for map centre)
-            response = get_stacked_response(obssim.obs(), None, None,
-                                            binsz=binsz, nxpix=npix, nypix=npix,
-                                            emin=emin, emax=emax, enumbins=nbins,
-                                            edisp=edisp,
-                                            coordsys=coord, proj=proj,
-                                            addbounds=addbounds,
-                                            log=log, debug=debug,
-                                            chatter=chatter)
-
-            # Set stacked response
-            if edisp:
-                binning.obs()[0].response(response['expcube'],
-                                          response['psfcube'],
-                                          response['edispcube'],
-                                          response['bkgcube'])
-            else:
-                binning.obs()[0].response(response['expcube'],
-                                          response['psfcube'],
-                                          response['bkgcube'])
-
-            # Set new models
-            binning.obs().models(response['models'])
-
-        # Make a deep copy of the observation that will be returned
-        # (the ctbin object will go out of scope one the function is
-        # left)
-        obs = binning.obs().copy()
+        # If a On source is specified then create On/Off observations
+        if onsrc != None:
+
+            # Extract source position from model
+            model = obssim.obs().models()[onsrc]
+            ra    = model['RA'].value()
+            dec   = model['DEC'].value()
+
+            # Allocate csphagen application and set parameters
+            phagen = cscripts.csphagen(obssim.obs())
+            phagen['ebinalg']     = 'LOG'
+            phagen['emin']        = emin
+            phagen['emax']        = emax
+            phagen['enumbins']    = nbins
+            phagen['coordsys']    = 'CEL'
+            phagen['ra']          = ra
+            phagen['dec']         = dec
+            phagen['rad']         = onrad
+            phagen['stack']       = False
+            phagen['inexclusion'] = 'NONE'
+
+            # Optionally open the log file
+            if log:
+                phagen.logFileOpen()
+
+            # Run csphagen application
+            phagen.run()
+
+            # Make a deep copy of the observation that will be returned
+            # (the csphagen object will go out of scope one the function is
+            # left)
+            obs = phagen.obs().copy()
+
+        # ... otherwise use binned observations
+        else:
+
+            # Allocate ctbin application and set parameters
+            binning = ctools.ctbin(obssim.obs())
+            binning['ebinalg']  = 'LOG'
+            binning['emin']     = emin
+            binning['emax']     = emax
+            binning['enumbins'] = nbins
+            binning['usepnt']   = True # Use pointing for map centre
+            binning['nxpix']    = npix
+            binning['nypix']    = npix
+            binning['binsz']    = binsz
+            binning['coordsys'] = coord
+            binning['proj']     = proj
+            binning['chatter']  = chatter
+            binning['debug']    = debug
+
+            # Optionally open the log file
+            if log:
+                binning.logFileOpen()
+
+            # Run ctbin application. This will loop over all observations in
+            # the container and bin the events in counts maps
+            binning.run()
+
+            # If we have multiple input observations then create stacked response
+            # cubes and append them to the observation
+            if len(obssim.obs()) > 1:
+
+                # Get stacked response (use pointing for map centre)
+                response = get_stacked_response(obssim.obs(), None, None,
+                                                binsz=binsz, nxpix=npix, nypix=npix,
+                                                emin=emin, emax=emax, enumbins=nbins,
+                                                edisp=edisp,
+                                                coordsys=coord, proj=proj,
+                                                addbounds=addbounds,
+                                                log=log, debug=debug,
+                                                chatter=chatter)
+
+                # Set stacked response
+                if edisp:
+                    binning.obs()[0].response(response['expcube'],
+                                             response['psfcube'],
+                                             response['edispcube'],
+                                             response['bkgcube'])
+                else:
+                    binning.obs()[0].response(response['expcube'],
+                                              response['psfcube'],
+                                              response['bkgcube'])
+
+                # Set new models
+                binning.obs().models(response['models'])
+
+            # Make a deep copy of the observation that will be returned
+            # (the ctbin object will go out of scope one the function is
+            # left)
+            obs = binning.obs().copy()
     else:
@@ -23,6 +23,8 @@ Bug fixes
 Improvements
 ------------
+* [`2260 <https://cta-redmine.irap.omp.eu/issues/2260>`_] -
+  Handle On/Off observations by :ref:`cspull`
 * [`2237 <https://cta-redmine.irap.omp.eu/issues/2237>`_] -
   Add ``statistic`` parameters to :ref:`ctbutterfly`, :ref:`cterror`,
   :ref:`cttsmap`, :ref:`ctulimit`, :ref:`cslightcrv, :ref:`csphasecrv`,