Add upper limit computation is csspec for NODES method (#2233)

Jürgen Knödlseder
1 parent a92dd304
Showing 4 changed files with 44 additions and 1 deletions
ChangeLog
NEWS
cscripts/csspec.py
doc/source/admin/release_history/1.6.rst
-2018-06-15
+2018-06-20
  
         * Version 1.6.0 released
         ========================
@@ -14,6 +14,7 @@
         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)
+        Add upper limit computation is csspec for NODES method (#2233)
         Define MET with respect to the observation time reference or "mjdref"
           parameter (#2054)
  
@@ -288,6 +288,10 @@ 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).
  
+csspec now computes upper limits also for the NODES method by fixing all node
+intensities and varying each node until the log-likelihood increases by a
+certain amount (#2233).
+
  
 cssrcdetect - Detects sources in sky map
 ----------------------------------------
@@ -758,6 +758,42 @@ class csspec(ctools.csobservation):
             result['flux']     = spectrum[i*2+1].value() * norm
             result['flux_err'] = spectrum[i*2+1].error() * norm
  
+            # Compute upper flux limit
+            ulimit_value = -1.0
+            if self['calc_ulim'].boolean():
+
+                # Logging information
+                self._log_header3(gammalib.EXPLICIT, 'Computing upper limit')
+
+                # Copy observation container
+                obs = like.obs().copy()
+
+                # Fix intensities of all nodes
+                spectral = obs.models()[self['srcname'].string()].spectral()
+                for par in spectral:
+                    par.fix()
+
+                # Create upper limit object  
+                ulimit = ctools.ctulimit(obs)
+                ulimit['srcname'] = self['srcname'].string()
+                ulimit['parname'] = 'Intensity%d' % i
+                ulimit['eref']    = elogmean.TeV()
+                ulimit['tol']     = 1.0e-3
+
+                # Try to run upper limit and catch exceptions
+                try:
+                    ulimit.run()
+                    ulimit_value = ulimit.diff_ulimit()
+                except:
+                    self._log_string(gammalib.EXPLICIT, 'Upper limit '
+                                     'calculation failed.')
+                    ulimit_value = -1.0
+
+                # Compute upper limit
+                if ulimit_value > 0.0:
+                    result['ulimit'] = ulimit_value * elogmean.MeV() * \
+                                       elogmean.MeV() * gammalib.MeV2erg
+
             # Compute TS
             if self['calc_ts'].boolean():
  
@@ -41,6 +41,8 @@ Improvements
   Allow specifying of input skymap to :ref:`ctskymap` tool
 * [`2309 <https://cta-redmine.irap.omp.eu/issues/2309>`_] -
   Accelerate `IRF` and `RING` background methods in :ref:`ctskymap`
+* [`2233 <https://cta-redmine.irap.omp.eu/issues/2233>`_] -
+  Add upper limit computation is :ref:`csspec` for ``NODES`` method
 * [`2045 <https://cta-redmine.irap.omp.eu/issues/2045>`_] -
   Define MET with respect to the observation time reference or the ``mjdref``
   user parameter (#2054)