Select relevant observations in cslightcrv

Jürgen Knödlseder
1 parent c3384d71
Showing 2 changed files with 74 additions and 47 deletions
cscripts/cslightcrv.py
cscripts/csobsselect.py
@@ -21,6 +21,7 @@
 import sys
 import gammalib
 import ctools
+import cscripts
 from cscripts import obsutils
 from cscripts import ioutils
 from cscripts import mputils
@@ -454,57 +455,79 @@ class cslightcrv(ctools.csobservation):
                   'values': {}}
         # Log information
-        self._log_header3(gammalib.EXPLICIT, 'Selecting events')
-
-        # Select events
-        select = ctools.ctselect(self.obs())
-        select['emin'] = self['emin'].real()
-        select['emax'] = self['emax'].real()
-        select['tmin'] = tmin
-        select['tmax'] = tmax
-        select['rad']  = 'UNDEFINED'
-        select['ra']   = 'UNDEFINED'
-        select['dec']  = 'UNDEFINED'
+        self._log_header3(gammalib.EXPLICIT, 'Selecting observations')
+
+        # Select observations
+        select = cscripts.csobsselect(self.obs())
+        select['pntselect'] = 'CIRCLE'
+        select['coordsys']  = 'GAL'
+        select['glon']      =   0.0
+        select['glat']      =   0.0
+        select['rad']       = 180.0
+        select['tmin']      = tmin
+        select['tmax']      = tmax
         select.run()
-        # Retrieve observation
+        # Retrieve observations
         obs = select.obs()
-        # Deal with stacked and On/Off Observations
-        if self._stacked or self._onoff:
+        # If there are observations then select now events from them
+        if obs.size() > 0:
-            # If a stacked analysis is requested bin the events
-            # and compute the stacked response functions and setup
-            # an observation container with a single stacked observation.
-            if self._stacked:
-                new_obs = obsutils.get_stacked_obs(self, obs)
+            # Log information
+            self._log_header3(gammalib.EXPLICIT, 'Selecting events')
-            # ... otherwise if On/Off analysis is requested generate
-            # the On/Off observations and response
-            elif self._onoff:
-                new_obs = obsutils.get_onoff_obs(self, obs, nthreads=1)
+            # Select events
+            select = ctools.ctselect(obs)
+            select['emin']  = self['emin'].real()
+            select['emax']  = self['emax'].real()
+            select['tmin']  = tmin
+            select['tmax']  = tmax
+            select['rad']   = 'UNDEFINED'
+            select['ra']    = 'UNDEFINED'
+            select['dec']   = 'UNDEFINED'
+            select.run()
-            # Extract models
-            models = new_obs.models()
+            # Retrieve observations
+            obs = select.obs()
-            # Fix background models if required
-            if self['fix_bkg'].boolean():
-                for model in models:
-                    if model.classname() != 'GModelSky':
-                        for par in model:
-                            par.fix()
+        # Continue only if there are observations
+        if obs.size() > 0:
-            # Put back models
-            new_obs.models(models)
+            # Deal with stacked and On/Off Observations
+            if self._stacked or self._onoff:
-            # Continue with new oberservation container
-            obs = new_obs
+                # If a stacked analysis is requested bin the events
+                # and compute the stacked response functions and setup
+                # an observation container with a single stacked observation.
+                if self._stacked:
+                    new_obs = obsutils.get_stacked_obs(self, obs)
-        # Header
-        self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
+                # ... otherwise if On/Off analysis is requested generate
+                # the On/Off observations and response
+                elif self._onoff:
+                    new_obs = obsutils.get_onoff_obs(self, obs, nthreads=1)
-        # Do maximum likelihood model fitting
-        if obs.size() > 0:
+                # Extract models
+                models = new_obs.models()
+
+                # Fix background models if required
+                if self['fix_bkg'].boolean():
+                    for model in models:
+                        if model.classname() != 'GModelSky':
+                            for par in model:
+                                par.fix()
+
+                # Put back models
+                new_obs.models(models)
+
+                # Continue with new oberservation container
+                obs = new_obs
+
+            # Header
+            self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
+
+            # Do maximum likelihood model fitting
             like = ctools.ctlike(obs)
             like['edisp']    = self['edisp'].boolean()
             like['nthreads'] = 1  # Avoids OpenMP conflict
@@ -533,7 +556,6 @@ class cslightcrv(ctools.csobservation):
                     result['values']['e_'+par] = source[par].error()
                 # Calculate upper limit (-1 if not computed)
-                #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs())
                 ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs)
                 if ul_diff > 0.0:
                     result['ul_diff']  = ul_diff
@@ -137,16 +137,18 @@ class csobsselect(ctools.csobservation):
             # Get observation start and stop time
             obs_gti = obs.gti()
-            tstart  = obs_gti.tstart()
-            tstop   = obs_gti.tstop()
             # Get User start and stop time in the time reference of the CTA
             # observations
-            gti = self._get_gti(obs_gti.reference())
-
-            # If the observation is not within the time selection then skip it
-            if gti.size() > 0 and not gti.contains(tstart) and \
-                                  not gti.contains(tstop):
+            gti    = self._get_gti(obs_gti.reference())
+            tstart = gti.tstart()
+            tstop  = gti.tstop()
+
+            # If there is a valid User start and stop time and if that start
+            # or stop time is not within the observation GTI, then skip the
+            # observation
+            if gti.size() > 0 and not obs_gti.contains(tstart) and \
+                                  not obs_gti.contains(tstop):
                 self._log_selection(obs, 'outside time interval')
             # ... otherwise select spatially
@@ -301,6 +303,9 @@ class csobsselect(ctools.csobservation):
         # Initialise empty observation container for selected observations
         selected_obs = gammalib.GObservations()
+        # Set models
+        selected_obs.models(self.obs().models())
+
         # Write input observation container into logger
         self._log_observations(gammalib.NORMAL, self.obs(), 'Input observation')