csphagen now creates appropriate model for likelihood fitting (#2711)

Jürgen Knödlseder
1 parent 679001d7
Showing 7 changed files with 173 additions and 68 deletions
ChangeLog
NEWS
cscripts/csphagen.par
cscripts/csphagen.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/csphagen.rst
test/test_csphagen.py
-2018-11-22
+2018-11-23
  
         * Version 1.6.0 released
         ========================
@@ -8,6 +8,7 @@
         Apply event bin weighting in ctmodel for stacked computation (#2721)
         Enable ctbin to produce counts cubes for a joint binned analysis (#2714)
         Add Jupyter notebooks for H.E.S.S. DR1 (#2712)
+        csphagen now creates appropriate model for likelihood fitting (#2711)
         csphagen now supports any kind of background model (#2710)
         Add code generator (#2709)
         Add csbkgmodel script (#2694)
 New Features and Important Changes in ctools 1.6.0
  
-22 November 2018
+23 November 2018
  
  
 Introduction
@@ -328,6 +328,11 @@ Use time reference of input observations if time selection is given in MET
  
 csphagen - Generate On/Off observations
 ---------------------------------------
+A "outmodel" parameter was added and the script now produces an appropriate
+background model for maximum likelihood model fitting. The script now also
+sets the appropriate fit statistic depending on the "use_model_bkg" parameter
+(#2711).
+
 The script now supports any kind of background model, not only IRF background
 models (#2710).
  
@@ -22,7 +22,8 @@ srcname,     s, a, Crab,,, &quot;Source name&quot;
 caldb,       s, a, prod2,,, "Calibration database"
 irf,         s, a, South_0.5h,,, "Instrument response function"
 inexclusion, f, h, NONE,,, "Sky map of exclusion regions"
-outobs,      s, a, onoff_obs.xml,,, "Output observation definition XML file"
+outobs,      f, a, onoff_obs.xml,,, "Output observation definition XML file"
+outmodel,    f, a, onoff_model.xml,,, "Output model definition XML file"
 prefix,      s, h, onoff,,, "Prefix for output files"
  
 #
@@ -49,7 +49,6 @@ class csphagen(ctools.csobservation):
         self._src_reg       = gammalib.GSkyRegions()
         self._models        = gammalib.GModels()
         self._srcname       = ''
-        self._use_model_bkg = True
         self._bkg_regs      = []
         self._excl_reg      = None
         self._has_exclusion = False
@@ -78,7 +77,6 @@ class csphagen(ctools.csobservation):
                  'src_reg'       : self._src_reg,
                  'models'        : self._models,
                  'srcname'       : self._srcname,
-                 'use_model_bkg' : self._use_model_bkg,
                  'bkg_regs'      : self._bkg_regs,
                  'excl_reg'      : self._excl_reg,
                  'has_exclusion' : self._has_exclusion,
@@ -105,7 +103,6 @@ class csphagen(ctools.csobservation):
         self._src_reg       = state['src_reg']
         self._models        = state['models']
         self._srcname       = state['srcname']
-        self._use_model_bkg = state['use_model_bkg']
         self._bkg_regs      = state['bkg_regs']
         self._excl_reg      = state['excl_reg']
         self._has_exclusion = state['has_exclusion']
@@ -261,12 +258,14 @@ class csphagen(ctools.csobservation):
         else:
             self._has_exclusion = False
  
-        # Query remaining parametersStacking
+        # Query remaining parameters
+        self['use_model_bkg'].boolean()
         self['stack'].boolean()
  
         # Query ahead output parameters
         if (self._read_ahead()):
-            self['outobs'].string()
+            self['outobs'].filename()
+            self['outmodel'].filename()
             self['prefix'].string()
  
         # If there are no observations in container then get them from the
@@ -292,7 +291,7 @@ class csphagen(ctools.csobservation):
             model.name('Dummy')
             self._models.append(model)
             self._srcname       = 'Dummy'
-            self._use_model_bkg = False
+            self['use_model_bkg'].boolean(False)
  
         # Return
         return
@@ -361,41 +360,109 @@ class csphagen(ctools.csobservation):
         # Return reflected regions
         return regions
  
-    def _set_models(self, obs):
+    def _set_models(self, results):
         """
-        Set models in observation container
+        Set models for On/Off fitting
  
-        The method replaces all "CTA", "HESS", "VERITAS", and "MAGIC"
-        background models by "CTAOnOff" background models.
+        The method does the following
+        - append "OnOff" to the instrument name of all background models
+        - fix all spatial and temporal parameters
  
         Parameters
         ----------
-        obs : `~gammalib.GObservations()`
-            Observation container
+        results : list of dict
+            Result dictionaries
  
         Returns
         -------
-        obs : `~gammalib.GObservations()`
-            Observation container
+        models : `~gammalib.GModels()`
+            Model container
         """
+        # Write header
+        self._log_header1(gammalib.NORMAL, 'Set models')
+
         # Initialise model container
         models = gammalib.GModels()
  
-        # Loop over all models and replace CTA/HESS/VERITAS/MAGIC background
-        # models by CTAOnOff background model
-        for model in self.obs().models():
+        # Initialies stacked model flag
+        has_stacked_model = False
+
+        # Loop over all models in observation and append "OnOff" to instrument
+        # names
+        for model in self._models:
+
+            # If model is a background model then append "OnOff"
             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')
+
+                # Skip model is background model should not be used
+                if not self['use_model_bkg'].boolean():
+                    self._log_string(gammalib.NORMAL, ' Skip "%s" model "%s" (%s)' % \
+                         (model.instruments(), model.name(), model.ids()))
+                    continue
+
+                # Check if model corresponds to one of the relevant
+                # observations
+                is_onoff = False
+                for result in results:
+                    if model.is_valid(result['instrument'], result['id']):
+                        is_onoff = True
+                        break
+
+                # If stacked analysis is requested then just use for model
+                # and remove instrument ID
+                if self['stack'].boolean():
+
+                    # If there is already a model for stacked analysis then
+                    # skip this one
+                    if has_stacked_model:
+                        msg = ' Skip "%s" model "%s" (%s). There is already ' \
+                              'a model for stacked analysis.' % \
+                              (model.instruments(), model.name(), model.ids())
+                        self._log_string(gammalib.NORMAL, msg)
+                        continue
+
+                    # ... otherwise use model for stacked analysis
+                    else:
+                        has_stacked_model = True
+                        model.ids('')
+
+                # Set instrument name
+                model.instruments(model.instruments()+'OnOff')
+
+            # Log model usage
+            self._log_string(gammalib.NORMAL, ' Use "%s" model "%s" (%s)' % \
+                 (model.instruments(), model.name(), model.ids()))
+
+            # Append model to container
             models.append(model)
  
-        # Append model to observation container
-        obs.models(models)
+        # Return model container
+        return models
  
-        # Return observation container
+    def _set_statistic(self, obs):
+        """
+        Set statistic for observation
+
+        If the "use_model_bkg" is true then set statistic to "cstat",
+        otherwise set it to "wstat"
+
+        Parameters
+        ----------
+        obs : `~gammalib.GObservation()`
+            Observation
+
+        Returns
+        -------
+        obs : `~gammalib.GObservation()`
+            Observation
+        """
+        # Set statistic according to background model usage
+        if self['use_model_bkg'].boolean():
+            obs.statistic('cstat')
+        else:
+            obs.statistic('wstat')
+
+        # Return observation
         return obs
  
     def _etrue_ebounds(self):
@@ -494,7 +561,7 @@ class csphagen(ctools.csobservation):
         # Otherwise calculate On/Off spectra
         else:
             # Log current observation
-            self._log_string(gammalib.NORMAL, 'Process observation "%s"' % \
+            self._log_string(gammalib.NORMAL, ' Process observation "%s"' % \
                              (obs.id()))
  
             # Set background regions for this observation
@@ -503,7 +570,7 @@ class csphagen(ctools.csobservation):
             # Get IRF background template usage flag
             use_model_bkg = self['use_model_bkg'].boolean()
             if not use_model_bkg:
-                self._log_string(gammalib.NORMAL, 'Background model not used, '
+                self._log_string(gammalib.NORMAL, ' Background model not used, '
                                  'assume constant backround rate.')
  
             # If there are background regions then create On/Off observation
@@ -520,13 +587,14 @@ class csphagen(ctools.csobservation):
                 onoff.id(obs.id())
  
             else:
-                self._log_string(gammalib.NORMAL, 'Observation %s not included '
+                self._log_string(gammalib.NORMAL, ' Observation %s not included '
                                  'in spectra generation' % (obs.id()))
  
         # Construct dictionary with results
-        result = {'onoff'   : onoff,
-                  'bkg_reg' : bkg_reg,
-                  'id'      : obs.id()}
+        result = {'onoff'     : onoff,
+                  'bkg_reg'   : bkg_reg,
+                  'instrument': obs.instrument(),
+                  'id'        : obs.id()}
  
         # Return results
         return result
@@ -549,12 +617,15 @@ class csphagen(ctools.csobservation):
         """
         # Continue only if result is valid
         if result['onoff'] != None:
-        
+
             # If the results contain an On/Off observation
             if result['onoff'].classname() == 'GCTAOnOffObservation':
  
+                # Set statistic according to background model usage
+                obs = self._set_statistic(result['onoff'])
+
                 # Append observation to observation container
-                outobs.append(result['onoff'])
+                outobs.append(obs)
  
                 # Append background regions
                 self._bkg_regs.append({'regions': result['bkg_reg'],
@@ -598,6 +669,7 @@ class csphagen(ctools.csobservation):
         # Initialise run variables
         outobs         = gammalib.GObservations()
         self._bkg_regs = []
+        results        = []
  
         # If there is more than one observation and we use multiprocessing
         if self._nthreads > 1 and self.obs().size() > 1:
@@ -609,7 +681,9 @@ class csphagen(ctools.csobservation):
  
             # Construct results
             for i in range(self.obs().size()):
-                outobs = self._unpack_result(outobs,poolresults[i][0])
+                result = poolresults[i][0]
+                outobs = self._unpack_result(outobs, result)
+                results.append(result)
                 self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
  
         # Otherwise, loop through observations and generate pha, arf, rmf files
@@ -618,9 +692,10 @@ class csphagen(ctools.csobservation):
                 # Process individual observation
                 result = self._process_observation(i)
                 outobs = self._unpack_result(outobs, result)
+                results.append(result)
  
         # Stack observations
-        if outobs.size() > 1 and self['stack'].boolean() == True:
+        if outobs.size() > 1 and self['stack'].boolean():
  
             # Write header
             self._log_header1(gammalib.NORMAL, 'Stacking %d observations' %
@@ -629,12 +704,18 @@ class csphagen(ctools.csobservation):
             # Stack observations
             stacked_obs = gammalib.GCTAOnOffObservation(outobs)
  
+            # Set statistic according to background model usage
+            stacked_obs = self._set_statistic(stacked_obs)
+
             # Put stacked observations in output container
             outobs = gammalib.GObservations()
             outobs.append(stacked_obs)
  
+        # Create models that allow On/Off fitting
+        models = self._set_models(results)
+
         # Set models in output container
-        outobs = self._set_models(outobs)
+        outobs.models(models)
  
         # Set observation container
         self.obs(outobs)
@@ -650,9 +731,10 @@ class csphagen(ctools.csobservation):
         self._log_header1(gammalib.TERSE, 'Save data')
  
         # Get XML output filename, prefix and clobber
-        outobs  = self['outobs'].string()
-        prefix  = self['prefix'].string()
-        clobber = self['clobber'].boolean()
+        outobs   = self['outobs'].filename()
+        outmodel = self['outmodel'].filename()
+        prefix   = self['prefix'].string()
+        clobber  = self['clobber'].boolean()
  
         # Loop over all observation in container
         for obs in self.obs():
@@ -694,8 +776,12 @@ class csphagen(ctools.csobservation):
         # Save observation definition XML file
         self.obs().save(outobs)
  
-        # Log file name
-        self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs)
+        # Save model definition XML file
+        self.obs().models().save(outmodel)
+
+        # Log file names
+        self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs.url())
+        self._log_value(gammalib.NORMAL, 'Model definition XML file', outmodel.url())
  
         # Save ds9 On region file
         regname = prefix + '_on.reg'
@@ -50,6 +50,8 @@ Improvements
   Rename ``outcube`` to ``outobs`` parameter.
 * [`2712 <https://cta-redmine.irap.omp.eu/issues/2712>`_] -
   Add Jupyter notebooks for H.E.S.S. DR1
+* [`2711 <https://cta-redmine.irap.omp.eu/issues/2711>`_] -
+  :ref:`csphagen` now creates appropriate model for likelihood fitting
 * [`2685 <https://cta-redmine.irap.omp.eu/issues/2685>`_] -
   Add ``ctbkgcube::cntcube()`` method
   Account for varying energy thresholds in :ref:`ctbin` weight computation
@@ -66,10 +66,6 @@ script will write out:
 General parameters
 ------------------
  
-
-``(nthreads = 0) [integer]``
-    Number of parallel processes (0=use all available CPUs).
-
 ``inobs [file]``
     Input event list or observation definition XML file.
  
@@ -91,9 +87,12 @@ General parameters
     Optional FITS file containing a WCS map in the first extension that defines
     sky regions not to be used for background estimation (where map value != 0).
  
-``outobs [string]``
+``outobs [file]``
     Output observation definition XML file.
  
+``outmodel [file]``
+    Output model XML file.
+
 ``(prefix = onoff) [string]``
     Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files.
  
@@ -203,6 +202,9 @@ General parameters
 Standard parameters
 -------------------
  
+``(nthreads = 0) [integer]``
+    Number of parallel processes (0=use all available CPUs).
+
 ``(chatter = 2) [integer]``
     Verbosity of the executable:
      ``chatter = 0``: no information will be logged
@@ -233,4 +235,3 @@ Related tools or scripts
 ------------------------
  
 None
-
@@ -92,7 +92,8 @@ class Test(test):
                          ' rad=0.2 stack="no" inexclusion="' + \
                          self._exclusion + \
                          '" bkgmethod="REFLECTED" '+ \
-                         'outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
+                         'outobs="csphagen_cmd1_obs.xml" prefix="csphagen_cmd1" ' + \
+                         'outmodel="csphagen_cmd1_model.xml" ' + \
                          'logfile="csphagen_cmd1.log" chatter=2'
  
         # Check if execution was successful
@@ -101,7 +102,7 @@ class Test(test):
  
         # Check output files
         self._check_output('csphagen_cmd1', self._nbins, self._nreg_with_excl)
-        self._check_outobs('csphagen_cmd1', 1)
+        self._check_outobs('csphagen_cmd1_obs.xml', 1)
  
         # Setup csphagen command
         cmd = csphagen + ' inobs="events_that_do_not_exist.fits" ' + \
@@ -114,7 +115,8 @@ class Test(test):
                          ' rad=0.2 stack="no" inexclusion="' + \
                          self._exclusion + \
                          '" bkgmethod="REFLECTED" '+ \
-                         'outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
+                         'outobs="csphagen_cmd2_obs.xml" prefix="csphagen_cmd2" ' + \
+                         'outmodel="csphagen_cmd2_model.xml" ' + \
                          'logfile="csphagen_cmd2.log" debug=yes chatter=1'
  
         # Check if execution failed
@@ -151,7 +153,8 @@ class Test(test):
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
-        phagen['outobs']      = 'csphagen_py1.xml'
+        phagen['outobs']      = 'csphagen_py1_obs.xml'
+        phagen['outmodel']    = 'csphagen_py1_model.xml'
         phagen['prefix']      = 'csphagen_py1'
         phagen['logfile']     = 'csphagen_py1.log'
         phagen['chatter']     = 1
@@ -161,7 +164,7 @@ class Test(test):
  
         # Check output
         self._check_output('csphagen_py1', self._nbins, self._nreg_with_excl)
-        self._check_outobs('csphagen_py1', 1)
+        self._check_outobs('csphagen_py1_obs.xml', 1)
  
         # Now test without exclusion region
         phagen = cscripts.csphagen()
@@ -182,7 +185,8 @@ class Test(test):
         phagen['etruemin']  = 0.05
         phagen['etruemax']  = 150.0
         phagen['etruebins'] = 5
-        phagen['outobs']    = 'csphagen_py2.xml'
+        phagen['outobs']    = 'csphagen_py2_obs.xml'
+        phagen['outmodel']  = 'csphagen_py2_model.xml'
         phagen['prefix']    = 'csphagen_py2'
         phagen['logfile']   = 'csphagen_py2.log'
         phagen['chatter']   = 2
@@ -192,7 +196,7 @@ class Test(test):
  
         # Check output
         self._check_output('csphagen_py2', self._nbins, self._nreg_wo_excl)
-        self._check_outobs('csphagen_py2', 1)
+        self._check_outobs('csphagen_py2_obs.xml', 1)
  
         # Test with multiple input observations, no stacking
         # Create observation container
@@ -220,7 +224,8 @@ class Test(test):
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
-        phagen['outobs']      = 'csphagen_py3.xml'
+        phagen['outobs']      = 'csphagen_py3_obs.xml'
+        phagen['outmodel']    = 'csphagen_py3_model.xml'
         phagen['prefix']      = 'csphagen_py3'
         phagen['logfile']     = 'csphagen_py3.log'
         phagen['chatter']     = 3
@@ -232,7 +237,7 @@ class Test(test):
         for s in range(2):
             self._check_output('csphagen_py3_' + str(s + 1), self._nbins,
                                self._nreg_mul[s])
-        self._check_outobs('csphagen_py3', 2)
+        self._check_outobs('csphagen_py3_obs.xml', 2)
  
         # Setup csphagen for test with multiple input observations and stacking
         phagen = cscripts.csphagen(obs)
@@ -251,7 +256,8 @@ class Test(test):
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
-        phagen['outobs']      = 'csphagen_py4.xml'
+        phagen['outobs']      = 'csphagen_py4_obs.xml'
+        phagen['outmodel']    = 'csphagen_py4_model.xml'
         phagen['prefix']      = 'csphagen_py4'
         phagen['logfile']     = 'csphagen_py4.log'
         phagen['chatter']     = 4
@@ -263,7 +269,7 @@ class Test(test):
         for s in range(2):
             self._check_output('csphagen_py4_stacked', self._nbins,
                                0, check_regions=False)
-        self._check_outobs('csphagen_py4', 1)
+        self._check_outobs('csphagen_py4_obs.xml', 1)
  
         # Setup csphagen for test with custom On and Off regions provided
         phagen = cscripts.csphagen()
@@ -282,7 +288,8 @@ class Test(test):
         phagen['etruemax']   = 150.0
         phagen['etruebins']  = 5
         phagen['stack']      = False
-        phagen['outobs']     = 'csphagen_py5.xml'
+        phagen['outobs']     = 'csphagen_py5_obs.xml'
+        phagen['outmodel']   = 'csphagen_py5_model.xml'
         phagen['prefix']     = 'csphagen_py5'
         phagen['logfile']    = 'csphagen_py5.log'
         phagen['chatter']    = 2
@@ -292,7 +299,7 @@ class Test(test):
  
         # Check output
         self._check_output('csphagen_py5', self._nbins, self._nreg_bkg_reg)
-        self._check_outobs('csphagen_py5', 1)
+        self._check_outobs('csphagen_py5_obs.xml', 1)
  
         # Append off regions to observation container
         for run in obs:
@@ -313,7 +320,8 @@ class Test(test):
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
-        phagen['outobs']      = 'csphagen_py6.xml'
+        phagen['outobs']      = 'csphagen_py6_obs.xml'
+        phagen['outmodel']    = 'csphagen_py6_model.xml'
         phagen['prefix']      = 'csphagen_py6'
         phagen['logfile']     = 'csphagen_py6.log'
         phagen['chatter']     = 4
@@ -325,7 +333,7 @@ class Test(test):
         for s in range(2):
             self._check_output('csphagen_py6_stacked', self._nbins,
                                0, check_regions=False)
-        self._check_outobs('csphagen_py6', 1)
+        self._check_outobs('csphagen_py6_obs.xml', 1)
  
         # Return
         return
@@ -358,7 +366,8 @@ class Test(test):
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
-        phagen['outobs']      = 'csphagen_py1_pickle.xml'
+        phagen['outobs']      = 'csphagen_py1_obs_pickle.xml'
+        phagen['outmodel']    = 'csphagen_py1_model_pickle.xml'
         phagen['prefix']      = 'csphagen_py1_pickle'
         phagen['logfile']     = 'csphagen_py1_pickle.log'
         phagen['chatter']     = 1
@@ -373,7 +382,7 @@ class Test(test):
  
         # Check result files
         self._check_output('csphagen_py1_pickle', self._nbins, self._nreg_with_excl)
-        self._check_outobs('csphagen_py1_pickle', 1)
+        self._check_outobs('csphagen_py1_obs_pickle.xml', 1)
  
         # Return
         return
@@ -544,7 +553,7 @@ class Test(test):
         Check the output XML file containing ON/OFF observations
         """
         # Load observation container
-        obs = gammalib.GObservations(filenameroot + '.xml')
+        obs = gammalib.GObservations(filenameroot)
  
         # Check container size
         self.test_value(obs.size(), nout, 'Check for ' + str(nout) +