script skeleton

Luigi Tibaldo
1 parent e9867c1d
Showing 2 changed files with 137 additions and 7 deletions
cscripts/csphagen.par
cscripts/csphagen.py
@@ -16,16 +16,21 @@
 #
 # Input and output parameters
 #============================
-inobs,  f, a, obs.xml,,, "Input event list or observation definition XML file"
+inobs,      f, a, obs.xml,,, "Input event list or observation definition XML file"
 exclusion,  f,  h, NONE,,, "FITS WCS map of exclusion regions"
-outroot, f, a, spectrum.pha,,, "Root of the name for output PHA and XML files"
-outreg, f, h, NONE,,, "Output ds9 reg file showing source/background regions"
+outroot,    f, a, spectrum,,, "Root of the name for output PHA and XML files"
+outreg,     f, h, NONE,,, "Output ds9 reg file showing source/background regions"
 #
 # Script parameters
 #==================
+emin,      r, a, 0.1,,, "Lower energy limit (TeV)"
+emax,      r, a, 100.0,,, "Upper energy limit (TeV)"
+enumbins,  i, a, 120,,, "Number of energy bins"
+ebinalg,   s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
+ebinfile,  f, a, NONE,,, "Name of the file containing the energy bin definition"
 bkgmethod, s, a, REFLECTED,REFLECTED,, "Method for background estimation"
-srcshape, s, a, CIRCLE,CIRCLE,, "Shape of the source region"
+srcshape,  s, a, CIRCLE,CIRCLE,, "Shape of the source region"
 coordsys,  s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
 ra,        r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
 dec,       r, a, 22.01,-90,90, "Declination of source region centre (deg)"
@@ -33,7 +38,7 @@ glon,      r, a, 184.56,0,360, &quot;Galactic longitude of source region centre (deg)
 glat,      r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
 rad,       r, a, 0.2,0,180, "Radius of source region circle (deg)"
 bkgregmin, i, h, 2,1,, "Minimum number of background region for REFLECTED method"
-stacking, b, a, yes,,, "Stack all runs into single pha or produce run-by-run files"
+stack,     b, a, yes,,, "Stack all runs into single pha (if no produce run-by-run files)"
 #
 # Standard parameters
@@ -19,5 +19,130 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-import sys
-import gammalib
 \ No newline at end of file
+import gammalib
+import ctools
+
+
+# =============== #
+# csfindobs class #
+# =============== #
+class csphagen(ctools.cscript):
+    """
+    Generate pha files for classical IACT spectral analysis
+    """
+
+    # Constructor
+    def __init__(self, *argv):
+        """
+        Constructor
+        """
+        # Set name, version
+        self._name = 'csfindobs'
+        self._version = ctools.__version__
+
+        # Initialise observation container from constructor arguments
+        self._obs, argv = self._set_input_obs(argv)
+
+        # Initialise application by calling the appropriate class constructor
+        self._init_cscript(argv)
+
+        # Return
+        return
+
+    # Private methods
+    def _get_parameters(self):
+        """
+        Get parameters from parfile and setup observations
+        """
+        # If there are no observations in container then query the inobs
+        # parameter
+        if self._obs.size() == 0:
+            self['inobs'].filename()
+
+        # Query energy binning parameters
+        self['ebinalg'].string()
+        if self['ebinalg'].string() == 'FILE':
+            self['ebinfile'].filename()
+        else:
+            self['emin'].real()
+            self['emax'].real()
+            self['enumbins'].integer()
+
+        # Initialise background estimation method
+        self._bkgmethod = self["bkgmethod"].string()
+
+        # Initialise source position/region querying relevant parameters
+        self._src_dir = gammalib.GSkyDir()
+        self._src_reg = gammalib.GSkyRegions()
+        coordsys = self['coordsys'].string()
+        if coordsys == "CEL":
+            ra = self['ra'].real()
+            dec = self['dec'].real()
+            self._src_dir.radec_deg(ra, dec)
+        elif coordsys == "GAL":
+            glon = self['glon'].real()
+            glat = self['glat'].real()
+            self._src_dir.lb_deg(glon, glat)
+        self._srcshape = self['srcshape'].string()
+        if self._srcshape == "CIRCLE":
+            self._rad = self['rad'].real()
+            self._src_reg.append(
+                gammalib.GSkyRegionCircle(self._src_dir, self._rad))
+
+        # Stacking
+        self._stack = self['stack'].bool()
+
+        # Query ahead output parameters
+        if (self._read_ahead()):
+            self['outroot'].string()
+
+        # Write input parameters into logger
+        self._log_parameters(gammalib.TERSE)
+
+        # Return
+        return
+
+    def _reflected_regions(self, obs):
+        """
+        Calculate list of reflected regions for a single observation
+        :param obs: observation
+        :return: list of reflected regions
+        """
+        pnt_dir = obs.pointing().dir()
+        offset = pnt_dir.dist_deg(self._src_dir)
+        outregions = []
+
+        return outregions
+
+    # Public methods
+    def run(self):
+        # Switch screen logging on in debug mode
+        if self._logDebug():
+            self._log.cout(True)
+
+        # Write observation into logger
+        self._log_observations(gammalib.NORMAL, self._obs, 'Observation')
+
+        # Set energy binning
+        if self['ebinalg'] == "FILE":
+            pass
+            ## need to load binning from file
+        else:
+            etrue = gammalib.GEbounds(self["enumbins"],
+                                      gammalib.GEnergy(self["emin"], 'TeV'),
+                                      gammalib.GEnergy(self["emax"], 'TeV'))
+            ereco = gammalib.GEbounds(self["enumbins"],
+                                      gammalib.GEnergy(self["emin"], 'TeV'),
+                                      gammalib.GEnergy(self["emax"], 'TeV'))
+
+        # Loop through observations and generate pha, arf, rmf files
+        outobs = gammalib.GObservations()
+        for obs in self._obs:
+            bkg_reg = gammalib.GSkyRegions()
+            if self._bkgmethod == "REFLECTED":
+                regions = self.reflected_regions(obs)
+            for region in regions:
+                bkg_reg.append(region)
+            onoff = gammalib.GCTAOnOffObservation(obs, etrue, ereco,
+                                                  self._src_region, bkg_reg)
+            outobs.append(onoff)