From f506071d2a3fa6e57ea4417035b66d1599758cf7 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 25 Jan 2019 23:09:44 +0100 Subject: [PATCH] Select relevant observations in cslightcrv --- cscripts/cslightcrv.py | 102 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------- cscripts/csobsselect.py | 19 ++++++++++++------- 2 files changed, 74 insertions(+), 47 deletions(-) diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index b43f870..14b81ce 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.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 diff --git a/cscripts/csobsselect.py b/cscripts/csobsselect.py index d157cf9..ab797f7 100755 --- a/cscripts/csobsselect.py +++ b/cscripts/csobsselect.py @@ -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') -- libgit2 0.21.4