Also compute TS in csspec for On/Off analysis using wstat statistic (#2312)

Jürgen Knödlseder
1 parent d9ea7da3
Showing 4 changed files with 30 additions and 16 deletions
ChangeLog
NEWS
cscripts/csspec.py
doc/source/admin/release_history/1.6.rst
-2018-06-14
+2018-06-15
  
         * Version 1.6.0 released
         ========================
@@ -10,8 +10,9 @@
         Keep existing library path in installation script on Mac OS X (#2420)
         csphagen produces Xspec compliant data (#2404)
         Add show_rmf.py example script
-        Allow specifying of input sky map in ctskymap (#2311)
         Add ras() and decs() methods to csobsinfo (#2313)
+        Also compute TS in csspec for On/Off analysis using wstat statistic (#2312)
+        Allow specifying of input sky map in ctskymap (#2311)
         Accelerate IRF and RING background methods in ctskymap (#2309)
         Define MET with respect to the observation time reference or "mjdref"
           parameter (#2054)
 New Features and Important Changes in ctools 1.6.0
  
-14 June 2018
+15 June 2018
  
  
 Introduction
@@ -282,7 +282,9 @@ None
  
 csspec - Generates Spectral Energy Distribution
 -----------------------------------------------
-None
+csspec now properly computes the Test Statistic for an On/Off analysis using
+the "wstat" statistic in the case that the input model definition XML file
+did not contain a background model (#2312).
  
  
 cssrcdetect - Detects sources in sky map
@@ -369,13 +369,15 @@ class csspec(ctools.csobservation):
                 spectrum = gammalib.GModelSpectralNodes(model.spectral(), energies)
                 spectrum.autoscale()
  
-                # Make sure that all nodes are positive
+                # Make sure that all nodes are positive. Autoscale all
+                # parameters so that their nominal value is unity.
                 for i in range(spectrum.nodes()):
                     par     = spectrum[i*2+1]
+                    par.autoscale()
                     value   = par.value()
                     minimum = 1.0e-20 * value
                     if minimum <= 0.0:
-                        minimum = 1.0e-40
+                        minimum = 1.0e-20
                         if minimum < value:
                             value = minimum
                     par.value(value)
@@ -384,6 +386,10 @@ class csspec(ctools.csobservation):
                 # Set spectral component of source model
                 model.spectral(spectrum)
  
+                # Make sure that TS computation is disabled (makes computation
+                # faster)
+                model.tscalc(False)
+
                 # Append model
                 models.append(model)
  
@@ -590,13 +596,14 @@ class csspec(ctools.csobservation):
         # Write header for fitting
         self._log_header3(gammalib.EXPLICIT, 'Performing fit')
  
-        # Perform maximum likelihood fit
+        # Perform maximum likelihood fit. Skip bins that have no events.
         like          = ctools.ctlike(obs)
         like['edisp'] = self['edisp'].boolean()
-        like.run()
+        if like.obs().nobserved() > 0:
+            like.run()
  
         # Continue only if log-likelihood is non-zero
-        if like.obs().logL() != 0.0:
+        if like.obs().logL() != 0.0 and like.obs().nobserved() > 0:
  
             # Get results
             fitted_models = like.obs().models()
@@ -726,7 +733,7 @@ class csspec(ctools.csobservation):
         # Extract fit results
         model    = like.obs().models()[self['srcname'].string()]
         spectrum = model.spectral()
-        ts       = model.ts()
+        logL0    = like.obs().logL()
  
         # Loop over all nodes
         for i in range(spectrum.nodes()):
@@ -757,12 +764,13 @@ class csspec(ctools.csobservation):
                 # Copy observation container
                 obs = like.obs().copy()
  
-                # Set intensity of node to tiny value
+                # Set intensity of node to tiny value by scaling the value
+                # by a factor 1e-8.
                 par = obs.models()[self['srcname'].string()].spectral()[i*2+1]
-                value = 1.0e-20 * par.value()
-                if par.min() > value:
-                    par.min(value)
-                par.value(value)
+                par.autoscale()
+                par.factor_min(1.0e-8)
+                par.factor_value(1.0e-8)
+                par.autoscale()
                 par.fix()
  
                 # Perform maximum likelihood fit
@@ -772,7 +780,8 @@ class csspec(ctools.csobservation):
  
                 # Store Test Statistic
                 model        = tslike.obs().models()[self['srcname'].string()]
-                result['TS'] = ts - model.ts()
+                logL1        = tslike.obs().logL()
+                result['TS'] = 2.0 * (logL1 - logL0)
  
             # Log information
             value = '%e +/- %e' % (result['flux'], result['flux_err'])
@@ -35,6 +35,8 @@ Improvements
 * Add ``show_rmf.py`` example script
 * [`2313 <https://cta-redmine.irap.omp.eu/issues/2313>`_] -
   Add `ras()` and `decs()` methods to :ref:`csobsinfo`
+* [`2312 <https://cta-redmine.irap.omp.eu/issues/2312>`_] -
+  Also compute TS in :ref:`csspec` for On/Off analysis using wstat statistic
 * [`2311 <https://cta-redmine.irap.omp.eu/issues/2311>`_] -
   Allow specifying of input skymap to :ref:`ctskymap` tool
 * [`2309 <https://cta-redmine.irap.omp.eu/issues/2309>`_] -