Johan Bregeon / ctools

#! /usr/bin/env python
# ==========================================================================
# Generation of an observation definition file.
#
# Copyright (C) 2015 Juergen Knoedlseder
#
# 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
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# ==========================================================================
import gammalib
import ctools
import sys


# ============== #
# csobsdef class #
# ============== #
class csobsdef(ctools.cscript):
    """
    The csobsdef class generates an observation definition file from a
    pointing list. The pointing list is a colon separated value ASCII
    file with header keywords in the first row followed by a list of
    pointings (one pointing per row). The following header keywords
    are supported (case sensitive, column order irrelevant):
    
    name     - Observation name string
    id       - Unique observation identifier string
    ra       - Right Ascension of pointing (deg)
    dec      - Declination of pointing (deg)
    lon      - Galactic longitude of pointing (deg)
    lat      - Galactic latitude of pointing (deg)
    duration - Duration of pointing (seconds)
    emin     - Lower energy limit (TeV)
    emax     - Upper energy limit (TeV)
    rad      - Radius of region of interest (deg)
    deadc    - Deadtime correction factor [0-1]
    caldb    - Calibration database
    irf      - Response function name
            
    Only the pairs (ra,dec) or (lon,lat) are mandatory header keywords.
    All other keywords are optional and can be specified when calling
    csobsdef as user parameters. The only exception is the "duration"
    keyword that will automatically queried.
    
    Here some usage examples:
    
    csobsdef
      Creates minimal observation definition file.

    csobsdef emin=0.1 emax=100.0
      Creates observation definition file with an energy range 100 GeV - 100 TeV.

    csobsdef rad=5
      Creates observation definition file with a ROI radius of 5 deg.

    csobsdef caldb=dummy irf=cta_dummy_irf
      Creates observation definition file using the "cta_dummy_irf" IRF in the
      "dummy" calibration database.
    """
    def __init__(self, *argv):
        """
        Constructor.
        """
        # Set name and version
        self.name    = "csobsdef"
        self.version = "1.0.0"

        # Initialise class members
        self.obs        = gammalib.GObservations()
        self.inpnt      = ""
        self.tmin       = 0.0
        self.outfile    = ""
        self.chatter    = 2
        self.clobber    = True
        self.debug      = False
        self.read_ahead = False
        
        # Make sure that parfile exists
        file = self.parfile()

        # Initialise application
        if len(argv) == 0:
            ctools.cscript.__init__(self, self.name, self.version)
        elif len(argv) ==1:
            ctools.cscript.__init__(self, self.name, self.version, *argv)
        else:
            raise TypeError("Invalid number of arguments given.")

        # Set logger properties
        self.log_header()
        self.log.date(True)

        # Return
        return
    
    def __del__(self):
        """
        Destructor.
        """
        # Return
        return

    def parfile(self):
        """
        Check if parfile exists. If parfile does not exist then create a
        default parfile. This kluge avoids shipping the cscript with a parfile.
        """
        # Set parfile name
        parfile = self.name+".par"

        # Test parfile existence, and create parfile if it does not exist
        try:
            pars = gammalib.GApplicationPars(parfile)
        except:
            # Signal that parfile was not found
            print("Parfile \""+parfile+"\" not found. Create default parfile.")
            
            # Create default parfile
            pars = gammalib.GApplicationPars()
            pars.append(gammalib.GApplicationPar("inpnt","f","a","NONE","","","Pointing definition file"))
            pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Output observation definition file"))
            pars.append(gammalib.GApplicationPar("duration","r","a","1800.0","","","Pointing duration (seconds)"))
            pars.append(gammalib.GApplicationPar("caldb","s","h","","","","Calibration database"))
            pars.append(gammalib.GApplicationPar("irf","s","h","","","","Instrument response function"))
            pars.append(gammalib.GApplicationPar("emin","r","h","UNDEF","","","Lower energy limit (TeV)"))
            pars.append(gammalib.GApplicationPar("emax","r","h","UNDEF","","","Upper energy limit (TeV)"))
            pars.append(gammalib.GApplicationPar("rad","r","h","UNDEF","","","Radius of ROI (degrees)"))
            pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
            pars.append_standard()
            pars.append(gammalib.GApplicationPar("logfile","f","h","csobsdef.log","","","Log filename"))
            pars.save(parfile)
        
        # Return
        return
        
    def get_parameters(self):
        """
        Get parameters from parfile.
        """
        # Read parameters    
        self.inpnt = self["inpnt"].filename()

        # Read ahead parameters
        if self.read_ahead:
            self.outfile = self["outobs"].filename()

        # Set some fixed parameters
        self.chatter = self["chatter"].integer()
        self.clobber = self["clobber"].boolean()
        self.debug   = self["debug"].boolean()

        # Return
        return
    
    def execute(self):
        """
        Execute the script.
        """
        # Set read ahead flag
        self.read_ahead = True
        
        # Run the script
        self.run()

        # Save observation definition file
        self.obs.save(self.outfile)
        
        # Return
        return

    def run(self):
        """
        Run the script.
        """
        # Switch screen logging on in debug mode
        if self.logDebug():
            self.log.cout(True)

        # Get parameters
        self.get_parameters()
        
        # Write input parameters into logger
        if self.logTerse():
            self.log_parameters()
            self.log("\n")

        # Load pointing definition file
        pntdef = gammalib.GCsv(self.inpnt, ",")
        ncols  = pntdef.ncols()
        npnt   = pntdef.nrows()-1

        # Throw an exception is there is no header information
        if pntdef.nrows() < 1:
            raise RuntimeError("No header found in pointing definition file.")

        # Clear observation container
        self.obs.clear()
        self.id = 1

        # Extract header from pointing definition file
        header = []
        for col in range(ncols):
            header.append(pntdef[0,col])

        # Loop over all pointings
        for pnt in range(npnt):

            # Set row index
            row = pnt + 1

            # Create CTA observation
            obs = gammalib.GCTAObservation()

            # Set observation name
            if "name" in header:
                name = pntdef[row, header.index("name")]
            else:
                name = "None"
            obs.name(name)

            # Set identifier
            if "id" in header:
                id = pntdef[row, header.index("id")]
            else:
                id       = "%6.6d" % self.id
                self.id += 1
            obs.id(id)

            # Set pointing
            if "ra" in header and "dec" in header:
                ra     = float(pntdef[row, header.index("ra")])
                dec    = float(pntdef[row, header.index("dec")])
                pntdir = gammalib.GSkyDir()
                pntdir.radec_deg(ra,dec)
            elif "lon" in header and "lat" in header:
                lon    = float(pntdef[row, header.index("lon")])
                lat    = float(pntdef[row, header.index("lat")])
                pntdir = gammalib.GSkyDir()
                pntdir.lb_deg(lon,lat)
            else:
                raise gammalib.GException.invalid_value(self.name,
                      "No (ra,dec) or (lon,lat) columns found in pointing"
                      " definition file.")
            obs.pointing(gammalib.GCTAPointing(pntdir))

            # Set response function
            if "caldb" in header:
                caldb = pntdef[row, header.index("caldb")]
            else:
                caldb = self["caldb"].string()
            if "irf" in header:
                irf = pntdef[row, header.index("irf")]
            else:
                irf = self["irf"].string()
            if caldb != "" and irf != "":
                obs = self.__set_response(obs, caldb, irf)

            # Set deadtime correction factor
            if "deadc" in header:
                deadc = float(pntdef[row, header.index("deadc")])
            else:
                deadc = self["deadc"].real()
            obs.deadc(deadc)

            # Set Good Time Interval
            if "duration" in header:
                duration = float(pntdef[row, header.index("duration")])
            else:
                duration = self["duration"].real()
            tmin      = self.tmin
            tmax      = self.tmin + duration
            gti       = gammalib.GGti(self.time_reference())
            tstart    = gammalib.GTime(tmin, self.time_reference())
            tstop     = gammalib.GTime(tmax, self.time_reference())
            self.tmin = tmax
            gti.append(tstart, tstop)
            obs.ontime(gti.ontime())
            obs.livetime(gti.ontime()*deadc)

            # Set Energy Boundaries
            has_emin = False
            has_emax = False
            if "emin" in header:
                emin     = float(pntdef[row, header.index("emin")])
                has_emin = True
            else:
                if self["emin"].is_valid():
                    emin     = self["emin"].real()
                    has_emin = True
            if "emax" in header:
                emax     = float(pntdef[row, header.index("emax")])
                has_emax = True
            else:
                if self["emax"].is_valid():
                    emax     = self["emax"].real()
                    has_emax = True
            has_ebounds = has_emin and has_emax
            if has_ebounds:
                ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"),
                                            gammalib.GEnergy(emax, "TeV"))

            # Set ROI
            has_roi = False
            if "rad" in header:
                rad     = float(pntdef[row, header.index("rad")])
                has_roi = True
            else:
                if self["rad"].is_valid():
                    rad     = self["rad"].real()
                    has_roi = True
            if has_roi:
                roi = gammalib.GCTARoi(gammalib.GCTAInstDir(pntdir), rad)

            # Create an empty event list
            list = gammalib.GCTAEventList()
            list.gti(gti)
            
            # Set optional information
            if has_ebounds:
                list.ebounds(ebounds)
            if has_roi:
                list.roi(roi)
            
            # Attach event list to CTA observation
            obs.events(list)

            # Append observation
            self.obs.append(obs)

        # Return
        return

    def __set_response(self, obs, caldb, irf):
        """
        Set response for an observation. We create an XML element for that
        so that the response XML writer will write the database and response
        name into the observation definition file.
        """
        # Create XML element
        xml = gammalib.GXmlElement()

        # Append parameter
        parameter = "parameter name=\"Calibration\" database=\""+\
                    caldb+"\" response=\""+irf+"\""
        xml.append(gammalib.GXmlElement(parameter))

        # Create CTA response
        response = gammalib.GCTAResponseIrf()
        response.read(xml)
        
        # Attach response to observation
        obs.response(response)
        
        # Return observation
        return obs


# ======================== #
# Main routine entry point #
# ======================== #
if __name__ == '__main__':
    """
    Generates observation definition file.
    """
    # Create instance of application
    app = csobsdef(sys.argv)
    
    # Open logfile
    app.logFileOpen()
    
    # Execute application
    app.execute()