Add "use_irf_bkg" parameter to cslightcrv

Jürgen Knödlseder
1 parent 6540305d
Showing 7 changed files with 82 additions and 55 deletions
ChangeLog
NEWS
cscripts/cslightcrv.par
cscripts/csphagen.py
cscripts/obsutils.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/cslightcrv.rst
-2018-09-11
+2018-09-19
  
         * Version 1.6.0 released
         ========================
  
+        Add "use_irf_bkg" parameter to cslightcrv
         Add "like_accuracy" parameter to likelihood tools and scripts (#2662)
         Do not skip bins with zero counts in csspec (#2671)
         Use "ENERG" instead of "ETRUE" in energy dispersion axis names in csroot2caldb (#2670)
 New Features and Important Changes in ctools 1.6.0
  
-11 September 2018
+19 September 2018
  
  
 Introduction
@@ -208,6 +208,8 @@ None
  
 cslightcrv - Generate light curve
 ---------------------------------
+Add hidden "use_irf_bkg" parameter for forwarding to csphagen.
+
 Add hidden "nthreads" parameter and implement parallelisation of light
 curve computation (#2421).
  
@@ -28,36 +28,37 @@ outfile,     f, a, lightcurve.fits,,, &quot;Output light curve file&quot;
 #
 # Script parameters
 #==================
-tbinalg,   s, a, GTI,FILE|LIN|GTI,, "Algorithm for defining time bins"
-tmin,      t, a, MJD 51544.5,,,"Lightcurve start time (UTC string, JD, MJD or MET in seconds)"
-tmax,      t, a, MJD 51544.6,,,"Lightcurve stop time (UTC string, JD, MJD or MET in seconds)"
-mjdref,    r, h, 51544.5,,, "Reference MJD for MET"
-tbins,     i, a, 5,1,10000, "Number of time bins"
-tbinfile,  f, a, NONE,,, "File defining the time binning"
-method,    s, a, 3D, 3D|ONOFF,, "Analysis method"
-emin,      r, a, 0.1,,, "Lower energy limit of events (TeV)"
-emax,      r, a, 100.0,,, "Upper energy limit of events (TeV)"
-enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned for 3D analysis only)"
-coordsys,  s, a, CEL,CEL|GAL,, "Coordinate System"
-proj,      s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
-xref,      r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
-yref,      r, a, 22.51,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
-nxpix,     i, a, 200,,, "Size of the X axis in pixels"
-nypix,     i, a, 200,,, "Size of the Y axis in pixels"
-binsz,     r, a, 0.02,,, "Pixel size (deg/pixel)"
-srcshape,  s, h, CIRCLE,CIRCLE,, "Shape of the source region"
-rad,       r, a, 0.2,0,180, "Radius of source region circle (deg)"
-bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
-bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
-maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
-etruemin,  r, h, 0.01,,, "Minimum true energy (TeV)"
-etruemax,  r, h, 300.0,,, "Maximum true energy (TeV)"
-etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
-statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
-calc_ts,   b, h, yes,,, "Compute TS value for each bin?"
-calc_ulim, b, h, yes,,, "Compute upper limit for each bin?"
-fix_srcs,  b, h, yes,,, "Fix other sky model parameters?"
-fix_bkg,   b, h, no,,, "Fix background model parameters?"
+tbinalg,     s, a, GTI,FILE|LIN|GTI,, "Algorithm for defining time bins"
+tmin,        t, a, MJD 51544.5,,,"Lightcurve start time (UTC string, JD, MJD or MET in seconds)"
+tmax,        t, a, MJD 51544.6,,,"Lightcurve stop time (UTC string, JD, MJD or MET in seconds)"
+mjdref,      r, h, 51544.5,,, "Reference MJD for MET"
+tbins,       i, a, 5,1,10000, "Number of time bins"
+tbinfile,    f, a, NONE,,, "File defining the time binning"
+method,      s, a, 3D, 3D|ONOFF,, "Analysis method"
+emin,        r, a, 0.1,,, "Lower energy limit of events (TeV)"
+emax,        r, a, 100.0,,, "Upper energy limit of events (TeV)"
+enumbins,    i, a, 0,0,100, "Number of energy bins for binned (0=unbinned for 3D analysis only)"
+coordsys,    s, a, CEL,CEL|GAL,, "Coordinate System"
+proj,        s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
+xref,        r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
+yref,        r, a, 22.51,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
+nxpix,       i, a, 200,,, "Size of the X axis in pixels"
+nypix,       i, a, 200,,, "Size of the Y axis in pixels"
+binsz,       r, a, 0.02,,, "Pixel size (deg/pixel)"
+srcshape,    s, h, CIRCLE,CIRCLE,, "Shape of the source region"
+rad,         r, a, 0.2,0,180, "Radius of source region circle (deg)"
+bkgmethod,   s, h, REFLECTED,REFLECTED,, "Method for background estimation"
+bkgregmin,   i, h, 2,1,1000, "Minimum number of background regions"
+use_irf_bkg, b, h, yes,,, "Use IRF background template?"
+maxoffset,   r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
+etruemin,    r, h, 0.01,,, "Minimum true energy (TeV)"
+etruemax,    r, h, 300.0,,, "Maximum true energy (TeV)"
+etruebins,   i, h, 30,,, "Number of bins per decade for true energy bins"
+statistic,   s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
+calc_ts,     b, h, yes,,, "Compute TS value for each bin?"
+calc_ulim,   b, h, yes,,, "Compute upper limit for each bin?"
+fix_srcs,    b, h, yes,,, "Fix other sky model parameters?"
+fix_bkg,     b, h, no,,, "Fix background model parameters?"
  
 #
 # Standard parameters
@@ -463,14 +463,14 @@ class csphagen(ctools.csobservation):
             On/Off spectra, background regions, observation id
         """
         # Retrieve observation from container
-        obs = self.obs()[i]
+        onoff   = None
+        bkg_reg = None
+        obs     = self.obs()[i]
  
         # Skip non CTA observations
         if obs.classname() != 'GCTAObservation':
             self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \
                              (obs.instrument(), obs.id()))
-            onoff = 0
-            bkg_reg = 0
  
         # Otherwise calculate On/Off spectra
         else:
@@ -527,15 +527,18 @@ class csphagen(ctools.csobservation):
         outobs : `~gammalib.GObservations`
             Observation container with result appended
         """
-        # If the results contain an On/Off observation
-        if result['onoff'].classname() == 'GCTAOnOffObservation':
-
-            # Append observation to observation container
-            outobs.append(result['onoff'])
-
-            # Append background regions
-            self._bkg_regs.append({'regions': result['bkg_reg'],
-                                   'id': result['id']})
+        # Continue only if result is valid
+        if result['onoff'] != None:
+        
+            # If the results contain an On/Off observation
+            if result['onoff'].classname() == 'GCTAOnOffObservation':
+
+                # Append observation to observation container
+                outobs.append(result['onoff'])
+
+                # Append background regions
+                self._bkg_regs.append({'regions': result['bkg_reg'],
+                                      'id': result['id']})
  
         # Return observation container
         return outobs
@@ -799,12 +799,13 @@ def get_onoff_obs(cls, obs, nthreads=0):
     if cls._logExplicit():
         cls._log.header3('Creating On/Off observations')
  
-    # Initialise inmodel, srcname and inexclusion
+    # Initialise inmodel, srcname, inexclusion and use_irf_bkg
     inmodel     = 'NONE'
     srcname     = ''
     inexclusion = 'NONE'
+    use_irf_bkg = True
  
-    # Set inmodel, srcname and inexclusion if possible
+    # Set inmodel, srcname, inexclusion and use_irf_bkg if possible
     if 'inmodel' in cls:
         if cls['inmodel'].is_valid():
             inmodel = cls['inmodel'].value()
@@ -813,6 +814,8 @@ def get_onoff_obs(cls, obs, nthreads=0):
     if 'inexclusion' in cls:
         if cls['inexclusion'].is_valid():
             inexclusion = cls['inexclusion'].value()
+    if cls.has_par('use_irf_bkg'):
+        use_irf_bkg = cls['use_irf_bkg'].boolean()
  
     # Create On/Off observations
     phagen = cscripts.csphagen(obs)
@@ -836,15 +839,16 @@ def get_onoff_obs(cls, obs, nthreads=0):
     phagen['bkgmethod'] = cls['bkgmethod'].string()
     if cls['bkgmethod'].string() == 'REFLECTED':
         phagen['bkgregmin'] = cls['bkgregmin'].integer()
-    phagen['maxoffset'] = cls['maxoffset'].real()
-    phagen['stack']     = True
-    phagen['etruemin']  = cls['etruemin'].real()
-    phagen['etruemax']  = cls['etruemax'].real()
-    phagen['etruebins'] = cls['etruebins'].integer()
-    phagen['chatter']   = cls['chatter'].integer()
-    phagen['clobber']   = cls['clobber'].boolean()
-    phagen['debug']     = cls['debug'].boolean()
-    phagen['nthreads']  = nthreads
+    phagen['use_irf_bkg'] = use_irf_bkg
+    phagen['maxoffset']   = cls['maxoffset'].real()
+    phagen['stack']       = True
+    phagen['etruemin']    = cls['etruemin'].real()
+    phagen['etruemax']    = cls['etruemax'].real()
+    phagen['etruebins']   = cls['etruebins'].integer()
+    phagen['chatter']     = cls['chatter'].integer()
+    phagen['clobber']     = cls['clobber'].boolean()
+    phagen['debug']       = cls['debug'].boolean()
+    phagen['nthreads']    = nthreads
     phagen.run()
  
     # Clone resulting observation container
@@ -32,6 +32,7 @@ Bug fixes
 Improvements
 ------------
  
+* Add ``use_irf_bkg`` parameter to :ref:`cslightcrv`
 * [`2671 <https://cta-redmine.irap.omp.eu/issues/2671>`_] -
   Do not skip bins with zero counts in :ref:`csspec`
 * [`2670 <https://cta-redmine.irap.omp.eu/issues/2670>`_] -
@@ -144,6 +144,21 @@ General parameters
     ``ONOFF`` analysis. If this number of background regions is not available the
     observation is skipped.
  
+``(use_irf_bkg = yes) [boolean]``
+    Specifies whether the IRF background template should be used for the
+    computation of the ``alpha`` parameter and the predicted background rate in
+    the Off region that is stored in the ``BACKRESP`` column of the Off spectrum
+    when using the ``ONOFF`` method.
+
+    If the parameter is set to ``no`` the IRF background template is not used and
+    the background rate is assumed identical within the On and Off regions. This
+    is the classical IACT analysis method that is used when using reflected Off
+    regions. In that case the ``alpha`` parameter becomes independent of energy
+    and only reflects the ratio between the solid angles of the On and Off
+    regions. The ``BACKRESP`` column in the Off spectrum will be filled with zeros.
+    The data need to be fitted with the ``wstat`` statistic, fitting with ``cstat``
+    will not work.
+
 ``(maxoffset = 4.0) [real]``
     Maximum offset in degrees of source from camera center to accept the
     observation for On/Off analysis.