Added On/Off analysis

Luigi Tibaldo
1 parent f28fae01
Showing 3 changed files with 168 additions and 126 deletions
cscripts/csscs.par
cscripts/csscs.py
cscripts/obsutils.py
@@ -18,14 +18,14 @@
 #============================
 inobs,       f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file"
 inmodel,     f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
-srcnames,    s, a, Crab,,, "Semicolon separated list of sources"
+srcnames,    s, a, Crab,,, "Semicolon-separated list of sources"
 expcube,     f, a, NONE,,, "Input exposure cube file"
 psfcube,     f, a, NONE,,, "Input PSF cube file"
 edispcube,   f, a, NONE,,, "Input energy dispersion cube file"
 bkgcube,     f, a, NONE,,, "Input background cube file"
 caldb,       s, a, prod2,,, "Calibration database"
 irf,         s, a, South_0.5h,,, "Instrument response function"
-inexclusion, f, h, NONE,,, "Exclusion regions"
+inexclusion, f, a, NONE,,, "Exclusion regions TODO: add meaningful default"
 edisp,       b, h, no,,, "Apply energy dispersion?"
 outfile,     f, a, map.fits,,, "Output map"
@@ -39,9 +39,9 @@ yref,          r, a, 22.51,-90,90, &quot;Second coordinate of image center in degrees
 nxpix,         i, a, 20,,, "Size of the X axis in pixels"
 nypix,         i, a, 20,,, "Size of the Y axis in pixels"
 binsz,         r, a, 0.1,,, "Pixel size (deg/pixel)"
-roisz,         r, a, 2.,1.42,1000, "Size of the region of interest (numbers of binsz)"
+rad,           r, a, 0.2,,, "Radius of the region for component separation (deg)"
 statistic,     s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
-method,        s, a, 3D, 3D|ONOFF,, "Analysis method"
+method,        s, a, UNBINNED, UNBINNED|ONOFF,, "Analysis method"
 ebinalg,       s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins"
 emin,          r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,          r, a, 100.0,,, "Upper energy limit (TeV)"
@@ -49,6 +49,7 @@ enumbins,      i, a, 20,,, &quot;Number of energy bins&quot;
 bkgmethod,     s, h, REFLECTED,REFLECTED,, "Method for background estimation"
 srcshape,      s, h, CIRCLE,CIRCLE,, "Shape of the source region"
 bkgregmin,     i, h, 2,1,1000, "Minimum number of background regions"
+bkgregskip,    i, h, 1,0,10, "Number of background regions to skip next to On region"
 use_model_bkg, b, h, yes,,, "Use background model?"
 maxoffset,     r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
 etruemin,      r, h, 0.01,,, "Minimum true energy (TeV)"
@@ -50,8 +50,8 @@ class csscs(ctools.csobservation):
         # Initialise data members
         self._nthreads = 0
-        self._roisz = 0.
         self._srcnames = []
+        self._method = None
         self._fits = None
         self._excl_reg_map = None # Exclusion region map for on/off analysis
@@ -72,7 +72,7 @@ class csscs(ctools.csobservation):
         # Set pickled dictionary
         state = {'base'     : ctools.csobservation.__getstate__(self),
                  'nthreads' : self._nthreads,
-                 'roisz'    : self._roisz,
+                 'method'   : self._method,
                  'srcnames' : self._srcnames,
                  'fits'     : self._fits,
                  'excl_reg_map' : self._excl_reg_map}
@@ -92,8 +92,8 @@ class csscs(ctools.csobservation):
         # Set state
         ctools.csobservation.__setstate__(self, state['base'])
         self._nthreads = state['nthreads']
-        self._roisz = state['roisz']
         self._srcnames = state['srcnames']
+        self._method = state['method']
         self._fits = state['fits']
         self._excl_reg_map = state['excl_reg_map']
@@ -102,6 +102,39 @@ class csscs(ctools.csobservation):
     # Private methods
+    def _get_onoff_parameters(self):
+        """
+        Get On/Off analysis parameters.
+        This is done here rather then in ctools.is_onoff
+        because the exclusion map needs to be handled differently
+        as an automatic parameter but queried only if not already set
+        and we need to verify is not empty.
+        Also many parameters are already queried for all methods
+        TODO: verify if this can be made more uniform with other scripts
+        """
+        # Exclusion map
+        if (self._excl_reg_map is not None) and (self._excl_reg_map.map().npix() > 0):
+            # Exclusion map set and is not empty
+            pass
+        elif self['inexclusion'].is_valid():
+            inexclusion = self['inexclusion'].filename()
+            self._excl_reg_map = gammalib.GSkyRegionMap(inexclusion)
+        else:
+            msg = 'csscs in On/Off mode requires input exclusion region.'
+            raise RuntimeError(msg)
+
+        # Other csphagen parameters
+        self["enumbins"].integer()
+        if self["bkgmethod"].string() == "REFLECTED":
+            self["bkgregmin"].integer()
+            self['bkgregskip'].integer()
+        self["maxoffset"].real()
+        self["etruemin"].real()
+        self["etruemax"].real()
+        self["etruebins"].integer()
+
+        return
+
     def _get_parameters(self):
         """
         Get parameters from parfile
@@ -111,6 +144,9 @@ class csscs(ctools.csobservation):
         if self.obs().is_empty():
             self.obs(self._get_observations())
+        # Set observation statistic
+        self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
+
         # Collect number of unbinned, binned and On/Off observations in
         # observation container
         n_unbinned = 0
@@ -124,24 +160,25 @@ class csscs(ctools.csobservation):
                     n_unbinned += 1
             elif obs.classname() == 'GCTAOnOffObservation':
                 n_onoff += 1
+        n_cta   = n_unbinned + n_binned + n_onoff
+        n_other = self.obs().size() - n_cta
-        # Check that we have at least one unbinned or binned CTA observation
-        if n_unbinned == 0 and n_binned == 0:
-            msg = 'No unbinned or binned CTA observations found. ' \
-                  'Please provide at least one unbinned or binned ' \
-                  'CTA observation.'
+        # Check if there are On/Off or non-IACT observations
+        if n_onoff > 0 or n_other > 0:
+            msg = 'On/Off or non-CTA observations found. '\
+                  'csscs does not support this type of observations.'
+            raise RuntimeError(msg)
+        # Otherwise check if there is a mix of binned and unbinned
+        elif n_unbinned > 0 and n_binned>0:
+            msg = 'Mix of unbinned and binned CTA observations ' \
+                  'found in observation container. csscs does not ' \
+                  'support this mix.'
             raise RuntimeError(msg)
-
-        # Set observation statistic
-        self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
         # Set models if there are none in the container
         if self.obs().models().size() == 0:
             self.obs().models(self['inmodel'].filename())
-        # Query exclusion region
-        self["inexclusion"].query()
-
         # Query the map definition parameters
         self['xref'].real()
         self['yref'].real()
@@ -149,29 +186,24 @@ class csscs(ctools.csobservation):
         self['proj'].string()
         self['nxpix'].integer()
         self['nypix'].integer()
+        self['binsz'].real()
-        # Compute the size of the ROIs for the analysis in deg
-        self._roisz = self['binsz'].real() * self['roisz'].real()
+        # Radius of ROI for component separation
+        self['rad'].real()
-        # Query energy boundaries
+        # Energy boundaries
         self["emin"].real()
         self["emax"].real()
-        # If we have unbinned observations query the analysis method
+        # If we have unbinned observations query analysis method
         if n_unbinned > 0:
-            # If method is ONOFF query some csphagen parameters
-            if self['method'].string() == 'ONOFF':
-                pass
-                # TODO: implement Onoff analysis
-                # self["enumbins"].integer()
-                # self['use_model_bkg'].boolean()
-                # if self["bkgmethod"].string() == "REFLECTED":
-                #     self["srcshape"].string()
-                #     self["bkgregmin"].integer()
-                #     self["maxoffset"].real()
-                #     self["etruemin"].real()
-                #     self["etruemax"].real()
-                #     self["etruebins"].integer()
+            self._method = self['method'].string()
+            # If method is Onoff query csphagen parameters
+            if self._method == 'ONOFF':
+                self._get_onoff_parameters()
+        # Otherwise set method to binned
+        else:
+            self._method = 'BINNED'
         # Query target source names
         srcnames = self['srcnames'].string()
@@ -190,7 +222,27 @@ class csscs(ctools.csobservation):
             msg = 'Not all target sources are present in input model.'
             raise RuntimeError(msg)
-        # Query the hidden parameters, just in case
+        # For On/Off analysis check the number of gamma-ray sources in the model
+        if self._method == 'ONOFF':
+            nsources = 0
+            # Verify if model is of type GModelSky
+            for model in self.obs().models():
+                try:
+                    gammalib.GModelSky(model)
+                    nsources +=1
+                except:
+                    pass
+            # If there are background gamma-ray sources in the model
+            # throw runtime error
+            if nsources > len(self._srcnames):
+                msg = 'Background gamma-ray sources found in the model. ' \
+                      'On/Off analysis does not support this feature.'
+                raise RuntimeError(msg)
+            # Otherwise continue
+            else:
+                pass
+
+        # Query other hidden parameters, just in case
         self['edisp'].boolean()
         self['calc_ulim'].boolean()
         self['calc_ts'].boolean()
@@ -243,6 +295,24 @@ class csscs(ctools.csobservation):
             # Log model name
             self._log_header3(gammalib.EXPLICIT, model.name())
+            # In On/Off mode we cannot support multiple source morphologies
+            # Thus, if we have more than one target we will set them all
+            # to isotropic
+            if self._method == 'ONOFF' and len(self._srcnames) > 1:
+                # Check if it is a source
+                try:
+                    gammalib.GModelSky(model)
+                    # In this case check if the spatial model is already isotropic
+                    if model.spatial() == 'DiffuseIsotropic':
+                        pass
+                    # Otherwise change it to isotropic
+                    else:
+                        msg = ' Setting spatial model to diffuse isotropic'
+                        self._log_string(gammalib.EXPLICIT, msg)
+                        model.spatial(gammalib.GModelSpatialDiffuseConst())
+                except:
+                    pass
+
             # Freeze all parameters except the normalization
             # which is assumed to be the first spectral parameter
             normpar = model.spectral()[0]
@@ -298,26 +368,29 @@ class csscs(ctools.csobservation):
         Returns
         -------
-        new_obs : `gammalib.GCTAObservation`
-            Output observation with masked cube
+        new_obs : `gammalib.GCTAObservations`
+            Output observations with masked cubes
         """
-        # Put observation in container
-        container = gammalib.GObservations()
-        container.append(obs)
-
         # Filter cube according to ROI and user energy range
-        cubemask = ctools.ctcubemask(container)
+        cubemask = ctools.ctcubemask(obs)
         cubemask['ra'] = ra
         cubemask['dec'] = dec
         cubemask['rad'] = rad
         cubemask['regfile'] = 'NONE'
         cubemask['emin'] = self['emin'].real()
         cubemask['emax'] = self['emax'].real()
+
+        # If chatter level is verbose and debugging is requested then
+        # switch also on the debug model in ctcubemask
+        if self._logVerbose() and self._logDebug():
+            cubemask['debug'] = True
+
+        # Mask cube
         cubemask.run()
-        # Get deep copy of filtered observation
-        new_obs = cubemask.obs()[0].copy()
+        # Get deep copy of filtered observations
+        new_obs = cubemask.obs().copy()
         # Return
         return new_obs
@@ -328,23 +401,19 @@ class csscs(ctools.csobservation):
         Parameters
         ----------
-        obs     : `gammalib.GCTAObservation` Input observation of type cube
+        obs     : `gammalib.GCTAObservations` Input observations of type event list
         ra      : `float` R.A. (deg)
         dec     : `float` Dec (deg)
         rad     : `float` radius (deg)
         Returns
         -------
-        new_obs : `gammalib.GCTAObservation`
-            Output observation with masked event list
+        new_obs : `gammalib.GCTAObservations`
+            Output observations with masked event lists
         """
-        # Put observation in container
-        container = gammalib.GObservations()
-        container.append(obs)
-
         # Filter event list according to ROI and user energy range
-        select = ctools.ctselect(container)
+        select = ctools.ctselect(obs)
         select['ra'] = ra
         select['dec'] = dec
         select['rad'] = rad
@@ -352,40 +421,21 @@ class csscs(ctools.csobservation):
         select['emax'] = self['emax'].real()
         select['tmin'] = 'INDEF'
         select['tmax'] = 'INDEF'
+
+        # If chatter level is verbose and debugging is requested then
+        # switch also on the debug model in ctcubemask
+        if self._logVerbose() and self._logDebug():
+            select['debug'] = True
+
+        # Select event list
         select.run()
         # Get deep copy of filtered observation
-        new_obs = select.obs()[0].copy()
+        new_obs = select.obs().copy()
         # Return
         return new_obs
-    def _mask_onoff(self,obs,ra,dec,rad):
-        """
-        Create On/Off observation
-        with On region matching the required mask
-
-        Parameters
-        ----------
-        obs     : `gammalib.GCTAObservation` Input observation of type cube
-        ra      : `float` R.A. (deg)
-        dec     : `float` Dec (deg)
-        rad     : `float` radius (deg)
-
-        Returns
-        -------
-        new_obs : `gammalib.GCTAOnoffObservation`
-            Output On/Off observation
-        """
-
-        # Put observation in container
-        container = gammalib.GObservations()
-        container.append(obs)
-
-        onoff_obs = obsutils.get_onoff_obs(self, container, ra=ra, dec=dec, rad=rad)
-
-        return onoff_obs[0]
-
     def _mask_observations(self,ra,dec,rad):
         """
         Create observations restricted to circular ROI
@@ -402,29 +452,16 @@ class csscs(ctools.csobservation):
             Observations in circular ROI
         """
-        # Create output observation container
-        new_obs = gammalib.GObservations()
-
-        # Loop over input observations
-        for obs in self.obs():
-            if obs.classname() == 'GCTAObservation':
-                if obs.eventtype() == 'CountsCube':
-                    masked_obs = self._mask_cube(obs,ra,dec,rad)
-                else:
-                    if self['method'].string() == '3D':
-                        masked_obs = self._mask_evlist(obs,ra,dec,rad)
-                    elif self['method'].string() == 'ONOFF':
-                        pass
-                        pass
-                        # TODO: implement Onoff analysis
-                        # masked_obs = self._mask_onoff(obs, ra, dec, rad)
-                new_obs.append(masked_obs)
-            # Skip On-Off and non-CTA observations
-            else:
-                pass
-
-        # Append models to output observations
-        new_obs.models(self.obs().models())
+        # Determine type of masking according
+        # to observation type and analysis method
+        if self._method == 'UNBINNED':
+            new_obs = self._mask_evlist(self.obs(), ra, dec, rad)
+        elif self._method == 'BINNED':
+            new_obs = self._mask_cube(self.obs(), ra, dec, rad)
+        elif self._method == 'ONOFF':
+            new_obs = obsutils.get_onoff_obs(self,self.obs(),nthreads=1,
+                                             ra = ra, dec = dec,
+                                             srcname = self._srcnames[0])
         # Return
         return new_obs
@@ -465,7 +502,7 @@ class csscs(ctools.csobservation):
         # Mask observations
         self._log_header3(gammalib.EXPLICIT, 'Masking observations')
-        masked_obs = self._mask_observations(ra,dec,self._roisz)
+        masked_obs = self._mask_observations(ra,dec,self['rad'].real())
         # Set up likelihood analysis
         self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin')
@@ -488,7 +525,7 @@ class csscs(ctools.csobservation):
         # ROI
         centre = gammalib.GSkyDir()
         centre.radec_deg(ra,dec)
-        roi = gammalib.GSkyRegionCircle(centre,self._roisz)
+        roi = gammalib.GSkyRegionCircle(centre,self['rad'].real())
         # Energy boundaries
         emin = gammalib.GEnergy(self['emin'].real(),'TeV')
         emax = gammalib.GEnergy(self['emax'].real(), 'TeV')
@@ -508,7 +545,7 @@ class csscs(ctools.csobservation):
                 # Multiply by flux from spatial component in ROI
                 flux *= source.spatial().flux(roi)
                 # Divide by solid angle of ROI
-                flux /= gammalib.twopi * (1 - math.cos(math.radians(self._roisz)))
+                flux /= gammalib.twopi * (1 - math.cos(math.radians(self['rad'].real())))
                 result[name]['flux'] = flux
                 # Get flux error
@@ -561,7 +598,7 @@ class csscs(ctools.csobservation):
         # Prepare some header cards
         # Flux units
         fcard = gammalib.GFitsHeaderCard('BUNIT', 'ph/cm2/s/sr' , 'Photon flux')
-        txt = 'Correlation radius %f deg' %(self._roisz)
+        txt = 'Correlation radius %f deg' %(self['rad'].real())
         roiszcard = gammalib.GFitsHeaderCard('COMMENT', '' , txt)
         # Loop over target sources
@@ -646,14 +683,6 @@ class csscs(ctools.csobservation):
             for i in range(npix):
                 results.append(self._pixel_analysis(i))
-        # # Loop over pixels and extract results
-        # results = []
-        # for inx in range(self._create_map().npix()):
-        #
-        #     # Analyse ROI
-        #     result = self._pixel_analysis(inx)
-        #     results.append(result)
-
         # Fill output Fits
         self._fits = self._fill_fits(results)
@@ -745,7 +745,7 @@ def get_stacked_obs(cls, obs, nthreads=0):
 # ================================ #
 # Get On/Off observation container #
 # ================================ #
-def get_onoff_obs(cls, obs, nthreads=0):
+def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = ''):
     """
     Create On/Off observations container from given observations
@@ -757,6 +757,12 @@ def get_onoff_obs(cls, obs, nthreads=0):
         Observation container
     nthreads : str, optional
         Number of parallel processes for On/Off spectra computation (0=all available CPUs)
+    ra : float, optional
+        R.A. of source region centre
+    dec : float, optional
+        Dec. of source region centre
+    srcname : str, optional
+        Source name
     Returns
     -------
@@ -767,9 +773,8 @@ def get_onoff_obs(cls, obs, nthreads=0):
     if cls._logExplicit():
         cls._log.header3('Creating On/Off observations')
-    # Initialise inmodel, srcname, inexclusion and use_model_bkg
+    # Initialise inmodel, inexclusion and use_model_bkg
     inmodel       = 'NONE'
-    srcname       = ''
     inexclusion   = 'NONE'
     use_model_bkg = True
@@ -777,7 +782,7 @@ def get_onoff_obs(cls, obs, nthreads=0):
     if cls.has_par('inmodel') and obs.models().size() == 0:
         if cls['inmodel'].is_valid():
             inmodel = cls['inmodel'].value()
-    if cls.has_par('srcname'):
+    if srcname == '' and cls.has_par('srcname'):
         srcname = cls['srcname'].value()
     if cls.has_par('inexclusion'):
         if cls['inexclusion'].is_valid():
@@ -795,13 +800,20 @@ def get_onoff_obs(cls, obs, nthreads=0):
     phagen['enumbins']    = cls['enumbins'].integer()
     phagen['ebinalg']     = 'LOG'
     phagen['srcshape']    = cls['srcshape'].string()
-    phagen['coordsys']    = cls['coordsys'].string()
-    if cls['coordsys'].string() == 'CEL':
-        phagen['ra']  = cls['xref'].real()
-        phagen['dec'] = cls['yref'].real()
-    elif cls['coordsys'].string() == 'GAL':
-        phagen['glon'] = cls['xref'].real()
-        phagen['glat'] = cls['yref'].real()
+    # User has specified custom centre for the source region
+    if ra != None and dec != None:
+        phagen['coordsys'] = 'CEL'
+        phagen['ra'] = ra
+        phagen['dec'] = dec
+    # Otherwise use default in class
+    else:
+        phagen['coordsys']    = cls['coordsys'].string()
+        if cls['coordsys'].string() == 'CEL':
+            phagen['ra']  = cls['xref'].real()
+            phagen['dec'] = cls['yref'].real()
+        elif cls['coordsys'].string() == 'GAL':
+            phagen['glon'] = cls['xref'].real()
+            phagen['glat'] = cls['yref'].real()
     if cls['srcshape'].string() == 'CIRCLE':
         phagen['rad'] = cls['rad'].real()
     phagen['bkgmethod'] = cls['bkgmethod'].string()