From fb6981b726c5f9212847a66bfc991ba399d6f81f Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 29 Mar 2019 15:03:10 +0100 Subject: [PATCH] Remove code duplication in csresmap and csresspec --- cscripts/csresmap.py | 27 +++------------------------ cscripts/csresspec.py | 33 ++++++++------------------------- cscripts/obsutils.py | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 62 insertions(+), 49 deletions(-) diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index 3f61273..f244392 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -2,7 +2,7 @@ # ========================================================================== # Generates a residual map. # -# Copyright (C) 2014-2018 Michael Mayer +# Copyright (C) 2014-2019 Michael Mayer # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -165,29 +165,8 @@ class csresmap(ctools.csobservation): # Write header self._log_header1(gammalib.TERSE, 'Generate binned map (ctbin)') - # Create countsmap - binning = ctools.ctbin(self.obs()) - binning['xref'] = self['xref'].real() - binning['yref'] = self['yref'].real() - binning['proj'] = self['proj'].string() - binning['coordsys'] = self['coordsys'].string() - binning['ebinalg'] = self['ebinalg'].string() - binning['nxpix'] = self['nxpix'].integer() - binning['nypix'] = self['nypix'].integer() - binning['binsz'] = self['binsz'].real() - if self['ebinalg'].string() == 'FILE': - binning['ebinfile'] = self['ebinfile'].filename().file() - else: - binning['enumbins'] = self['enumbins'].integer() - binning['emin'] = self['emin'].real() - binning['emax'] = self['emax'].real() - binning['chatter'] = self['chatter'].integer() - binning['clobber'] = self['clobber'].boolean() - binning['debug'] = self['debug'].boolean() - binning.run() - - # Retrieve counts cube - cta_counts_cube = binning.cube() + # Create counts cube + cta_counts_cube = obsutils.create_counts_cube(self, self.obs()) # Assign GCTAEventCube to skymap countmap = cta_counts_cube.counts() diff --git a/cscripts/csresspec.py b/cscripts/csresspec.py index ac34a9d..2277a36 100755 --- a/cscripts/csresspec.py +++ b/cscripts/csresspec.py @@ -18,10 +18,10 @@ # along with this program. If not, see . # # ========================================================================== +import sys import gammalib import ctools from cscripts import obsutils -import sys # =============== # @@ -672,35 +672,18 @@ class csresspec(ctools.csobservation): # Log processing header self._log_header1(gammalib.TERSE, 'Processing Observations') - #If observations are unbinned and we stack + # If observations are unbinned and we stack if self._stack and self.obs()[0].classname() == 'GCTAObservation': + # Log counts cube computation msg = 'Computing count cube from multiple unbinned observations' self._log_string(gammalib.NORMAL, msg) - # Build count cube - binning = ctools.ctbin(self.obs()) - binning['xref'] = self['xref'].real() - binning['yref'] = self['yref'].real() - binning['proj'] = self['proj'].string() - binning['coordsys'] = self['coordsys'].string() - binning['ebinalg'] = self['ebinalg'].string() - binning['nxpix'] = self['nxpix'].integer() - binning['nypix'] = self['nypix'].integer() - binning['binsz'] = self['binsz'].real() - if self['ebinalg'].string() == 'FILE': - binning['ebinfile'] = self['ebinfile'].filename().file() - else: - binning['enumbins'] = self['enumbins'].integer() - binning['emin'] = self['emin'].real() - binning['emax'] = self['emax'].real() - binning['chatter'] = self['chatter'].integer() - binning['clobber'] = self['clobber'].boolean() - binning['debug'] = self['debug'].boolean() - binning.run() - - #compute residuals using cube - table = self._residuals_3D(self.obs(),'',binning.cube()) + # Create counts cube + cta_counts_cube = obsutils.create_counts_cube(self, self.obs()) + + # Compute residuals using cube + table = self._residuals_3D(self.obs(), '', cta_counts_cube) # Append results table to output file self._fits.append(table) diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index 1d110cc..2a1c86b 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -937,3 +937,54 @@ def residuals(cls, counts, model): # Return return residuals + + +# ================= # +# Create counts map # +# ================= # +def create_counts_cube(cls, obs): + """ + Create counts cube from observations + + Parameters + ---------- + cls : `~ctools.cscript` + cscript class + obs : `~gammalib.GObservations` + Observation container + + Returns + ------- + cntcube : `~gammalib.GCTAEventCube` + Counts cube + """ + # Initialise ctbin + ctbin = ctools.ctbin(obs) + + # Set parameters + ctbin['xref'] = cls['xref'].real() + ctbin['yref'] = cls['yref'].real() + ctbin['proj'] = cls['proj'].string() + ctbin['coordsys'] = cls['coordsys'].string() + ctbin['ebinalg'] = cls['ebinalg'].string() + ctbin['nxpix'] = cls['nxpix'].integer() + ctbin['nypix'] = cls['nypix'].integer() + ctbin['binsz'] = cls['binsz'].real() + if cls['ebinalg'].string() == 'FILE': + ctbin['ebinfile'] = cls['ebinfile'].filename().file() + else: + ctbin['enumbins'] = cls['enumbins'].integer() + ctbin['emin'] = cls['emin'].real() + ctbin['emax'] = cls['emax'].real() + ctbin['chatter'] = cls['chatter'].integer() + ctbin['clobber'] = cls['clobber'].boolean() + ctbin['debug'] = cls['debug'].boolean() + + # Run ctbin + ctbin.run() + + # Retrieve counts cube + cntcube = ctbin.cube().copy() + + # Return counts cube + return cntcube -- libgit2 0.21.4