Merge branch '1579-add-csobs2caldb' into integration

# Conflicts: # cscripts/Makefile.am # doc/source/reference_manual/index.rst

Merge branch '1579-add-csobs2caldb' into integration
# Conflicts: # cscripts/Makefile.am # doc/source/reference_manual/index.rst
Jürgen Knödlseder
2 parents 44a68b8f 22bfe352
Showing 8 changed files with 499 additions and 6 deletions
.gitignore
cscripts/Makefile.am
cscripts/csobs2caldb.py
doc/source/reference_manual/csobs2caldb.rst
doc/source/reference_manual/index.rst
test/data/obs_unbinned.xml
test/irf/irf_file.fits.gz
test/test_cscripts.sh.in
@@ -63,6 +63,7 @@ src/ctbutterfly/ctbutterfly
 src/ctulimit/ctulimit
 src/cterror/cterror
 test/pfiles/
+test/caldb/
 test/*.fits
 test/*.dat
 test/*.xml
@@ -5,12 +5,14 @@ pythonscripts = __init__.py \
                 cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
                 cshessobs.py csobsdef.py csspec.py cslightcrv.py csobsinfo.py \
                 csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \
+                csobs2caldb.py \
                 obsutils.py
  
 # Scripts for which symbolic links will be created in bindir
 linkscripts = cspull cssens cstsdist cscaldb csresmap \
               cshessobs csobsdef csspec cslightcrv csobsinfo \
               csinfo csmodelinfo csmodelmerge cstsmapmerge
+              csobs2caldb
  
 # Distribute scripts
 EXTRA_DIST         = $(binscripts) $(pythonscripts)
+#!/usr/bin/env python
+# ==========================================================================
+# Generation of an caldb entry from general IACT observation.
+#
+# Copyright (C) 2015 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
+# 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
+import os
+import datetime
+
+
+# ================= #
+# csobs2caldb class #
+# ================= #
+class csobs2caldb(ctools.cscript):
+    """
+    This class implements the creation of a caldb entry for a particular
+    observation, which might be helpful for running simulations for current
+    IACTs. It takes an observation definition file as input and uses the
+    first observation to create a caldb entry which can be used for
+    simulations and sensitivity studies later on.
+    """
+    def __init__(self, *argv):
+        """
+        Constructor.
+        """
+        # Set name
+        self.name    = "csobs2caldb"
+        self.version = "1.1.0"
+
+        # 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"
+
+        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("inobs","f","a","obs.xml","","","Input observation definition file"))
+            pars.append(gammalib.GApplicationPar("rspname","s","a","NONE","","","Response output name (e.g. Zenith50)"))            
+            pars.append(gammalib.GApplicationPar("outirfs","f","h","irf_file.fits","","","Output IRF file name"))            
+            pars.append(gammalib.GApplicationPar("index","i","h","0","","","Index of observation to be used from XML file"))            
+            pars.append_standard()
+            pars.append(gammalib.GApplicationPar("logfile","f","h","csobs2caldb.log","","","Log filename"))
+            pars.save(parfile)
+
+        # Return
+        return
+
+    def get_parameters(self):
+        """
+        Get parameters from parfile and setup the observation.
+        """
+        # Get parameters     
+        self.inobs = self["inobs"].filename()
+        self.obs   = gammalib.GObservations(self.inobs)
+        
+        # Get index of observation to be used (0=default)
+        inx = self["index"].integer()
+        
+        # Get observation
+        self.observation = self.obs[inx]
+        
+        # Set instrument mission to cta
+        # This is to be compliant with standard gammalib caldb
+        self.m_mission = "CTA"
+        
+        # Get instrument name
+        self.m_instrument = self.observation.instrument()
+        
+        # Get config and file name
+        self.m_rspname = self["rspname"].string()
+        self.m_outfile = self["outirfs"].filename()
+        
+        # Make sure we have a cta observation
+        if not self.observation.classname() == "GCTAObservation":
+            raise gammalib.Exception.runtime_error("Input observation not of type \"GCTAObservation\"")
+        
+        # Make sure we have an irf response
+        if not self.observation.response().classname() == "GCTAResponseIrf":
+            raise gammalib.Exception.runtime_error("Response of input observation not of type \"GCTAResponseIrf\"")
+
+        # Get other parameters
+        self.m_clobber = self["clobber"].boolean()
+
+        # Return
+        return
+
+    def make_irf_file(self):
+        
+        # retrieve response component
+        rsp = self.observation.response()
+    
+        # Initialise FITS file
+        fits = gammalib.GFits()
+    
+        # Open FITS files of response components
+        fits_aeff  = gammalib.GFits(rsp.aeff().filename())
+        fits_psf   = gammalib.GFits(rsp.psf().filename())
+        fits_edisp = gammalib.GFits(rsp.edisp().filename())
+        fits_bkg   = gammalib.GFits(rsp.background().filename())
+        
+        # Bundle IRFs into one file
+        fits.append(fits_aeff["EFFECTIVE AREA"])
+        fits.append(fits_psf["POINT SPREAD FUNCTION"])
+        fits.append(fits_edisp["ENERGY DISPERSION"])
+        fits.append(fits_bkg["BACKGROUND"])
+        
+        # Return fits file
+        return fits
+
+    def make_caldb_index(self):
+        
+        # Open calibration database index
+        indx_file = self.base_dir+"/caldb.indx"
+        
+        # Check if indx file exists (create if not existent)
+        if os.path.isfile(indx_file):
+            cif = gammalib.GFits(indx_file)
+        else:
+            cif = gammalib.GFits(indx_file, True)
+        
+        # Retrieve "CIF" table
+        if cif.contains("CIF"):
+            table = cif.table("CIF")
+        else:       
+            # Create binary table if not available
+            bintable = gammalib.GFitsBinTable()
+            bintable.extname("CIF")
+            bintable.card("CIFVERSN", "1992a", "Version of CIF format")     
+            cif.append(bintable)
+            
+            table = cif.table("CIF")
+            
+            # Copied from cta_root2caldb.py
+            # Attach columns. Reference: CAL/GEN/92-008
+            table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
+            table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
+            table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20))
+            table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10))
+            table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20))
+            table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70))
+            table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40))
+            table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3))
+            table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4))
+            table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20))
+            table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9))
+            table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0))
+            table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10))
+            table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8))
+            table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0))
+            table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
+            table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
+            table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
+
+        # Check if output config already exist
+        has_config = False
+        row_index = -1
+        for caldir in table["CAL_DIR"]:
+            row_index += 1
+            if caldir == self.cal_dir:
+                has_config = True
+                break
+         
+        # Create columns if not available   
+        if not has_config:           
+            # Append 4 rows to CIF extension
+            table.append_rows(4)
+            row_index = table.nrows()
+        else:
+            row_index += 4
+
+        # Add generic information for these 4 rows
+        for i in range(4):
+
+            # Set row number
+            row = i + row_index-4
+
+            # Set date
+            now = str(datetime.datetime.now())
+            today, time = now.split()
+            
+            # Set element
+            table["TELESCOP"][row] = self.m_mission
+            table["INSTRUME"][row] = self.m_instrument
+            table["DETNAM"][row]   = "NONE"
+            table["FILTER"][row]   = "NONE"
+            table["CAL_DEV"][row]  = "ONLINE"
+            table["CAL_CLAS"][row] = "BCF"
+            table["CAL_DTYP"][row] = "DATA"
+            table["CAL_VSD"][row]  = today
+            table["CAL_VST"][row]  = time.split(".")[0]
+            table["REF_TIME"][row] = 51544.0
+            table["CAL_QUAL"][row] = 0
+            table["CAL_CBD"][row]  = "NAME("+self.m_rspname+")"
+            table["CAL_DATE"][row] = today.replace("-","/")[2:]
+            table["CAL_DIR"][row]  = self.cal_dir
+            table["CAL_FILE"][row] = self.m_outfile
+
+        # Add effective area information
+        row = row_index-4
+        table["CAL_CNAM"][row]  = "EFF_AREA"
+        table["CAL_DESC"][row]  = self.m_instrument+" effective area"
+
+        # Add point spread function information
+        row = row_index-3
+        table["CAL_CNAM"][row]  = "RPSF"
+        table["CAL_DESC"][row]  = self.m_instrument+" point spread function"
+        
+        # Add energy dispersion information
+        row = row_index-2
+        table["CAL_CNAM"][row]  = "EDISP"
+        table["CAL_DESC"][row]  = self.m_instrument+" energy dispersion"
+
+        # Add background information
+        row = row_index-1
+        table["CAL_CNAM"][row]  = "BGD"
+        table["CAL_DESC"][row]  = self.m_instrument+" background"
+        
+        # Return CIF FITS file
+        return cif
+        
+    def make_dirs(self):
+
+        # Create gammalib Caldb        
+        caldb = gammalib.GCaldb()
+        
+        # Set calibration directory 
+        self.cal_dir  = "data"
+        self.cal_dir += "/"+self.m_mission.lower()
+        self.cal_dir += "/"+self.m_instrument.lower()
+        self.cal_dir += "/bcf/"+self.m_rspname
+        
+        # Set absolute path
+        self.base_dir = caldb.rootdir() +"/data"
+        self.base_dir += "/"+self.m_mission.lower()
+        self.base_dir += "/"+self.m_instrument.lower()
+        
+        # Set directory for irf file
+        self.rsp_dir = caldb.rootdir() + "/" + self.cal_dir
+        
+        # Create directories and log information
+        if not os.path.isdir(self.rsp_dir):
+            if self.logExplicit():
+                self.log(gammalib.parformat("Directory"))
+                self.log(self.rsp_dir)
+                self.log("\n")
+            os.makedirs(self.rsp_dir)
+        else:
+            if self.logExplicit():
+                self.log(gammalib.parformat("Directory (existing)"))
+                self.log(self.rsp_dir)
+                self.log("\n")
+
+        # 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\n")
+            self.log.header1("Creating directory structure")
+        
+        # Create directory structure
+        self.make_dirs()
+
+        # Logging
+        if self.logTerse():
+            self.log("\n")
+            self.log.header1("Updating calibration database")
+         
+        # Create/update calibration database   
+        self.caldb_inx = self.make_caldb_index()
+        
+        # Logging
+        if self.logTerse():
+            self.log("\n")
+            self.log.header1("Creating response file")
+            
+        # Create IRF file
+        self.irf_fits = self.make_irf_file()
+                
+        # Return
+        return
+
+    def save(self):
+
+        # Save caldb index file
+        self.caldb_inx.save(self.m_clobber)
+
+        # Save response file
+        self.irf_fits.saveto(self.rsp_dir + "/" + self.m_outfile, self.m_clobber)
+        
+        # Return
+        return
+        
+    def execute(self):
+        """
+        Execute the script.
+        """
+        # Run the script
+        self.run()
+
+        # Save to disk
+        self.save()
+
+        # Return
+        return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Generates caldb entry.
+    """
+    # Create instance of application
+    app = csobs2caldb(sys.argv)
+    
+    # Open logfile
+    app.logFileOpen()
+
+    # Execute application
+    app.execute()
+.. _csobs2caldb:
+
+csobs2caldb
+===========
+
+Creates a caldb entry from an input observation.
+
+
+Synopsis
+--------
+
+This script creates an entry in the local calibration database from an
+input observation container. The caldb entry can be used for simulations
+with :doc:`ctobssim` later on. On default csobs2caldb uses the first observation
+in the given container. This can be modified using the hidden parameter "index"
+(e.g. index=3 instructs the tool to use the third observation in the container).
+Note that for the moment, the script only works for unbinned observations which 
+use the latest CTA IRF format (i.e. 2D/3D FITS binary tables).  
+
+In particular this script might be useful for simulations of current IACT data.
+Provided a specific observation with IRFs e.g. at certain zenith/azimuth angle, 
+the user can create a caldb entry using this tool and subsequently simulate e.g. 
+50 hours using this configuration. 
+
+General parameters
+------------------
+
+``inobs [file]``
+    Input event list, counts cube or observation definition XML file.
+
+``rspname [string]``
+    Response output name (e.g. Zenith50).
+
+``(outirfs = irf_file.fits) [file]``
+    Output IRF file name.
+
+``(index = 0) [integer]``
+    Index of observation to be used from XML file.
+
+
+Standard parameters
+-------------------
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     chatter = 0: no information will be logged
+     
+     chatter = 1: only errors will be logged
+     
+     chatter = 2: errors and actions will be logged
+     
+     chatter = 3: report about the task execution
+     
+     chatter = 4: detailed report about the task execution
+ 	 	 
+``(clobber = yes) [boolean]``
+    Specifies whether an existing output counts cube should be overwritten.
+ 	 	 
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+ 	 	 
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = csobs2caldb.log) [filename]``
+    Log filename.
+
+
+Related tools or scripts
+------------------------
+
+:doc:`cscaldb`
@@ -49,3 +49,4 @@ cscripts
    csmodelinfo --- Shows model container content <csmodelinfo>
    csmodelmerge --- Merges several model containers into one file <csmodelmerge>
    cstsmapmerge --- Merges slices from ts map computations <cstsmapmerge>
+   csobs2caldb --- Creates a caldb entry from an input observation <csobs2caldb>
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<observation_list title="observation list">
+  <observation name="" id="00001" instrument="CTA">
+    <parameter name="EventList"           file="data/crab_events.fits" />
+    <parameter name="EffectiveArea"       file="irf/irf_file.fits" />
+    <parameter name="PointSpreadFunction" file="irf/irf_file.fits" />
+    <parameter name="EnergyDispersion"    file="irf/irf_file.fits" />
+    <parameter name="Background"          file="irf/irf_file.fits" />
+  </observation>
+</observation_list>
@@ -42,11 +42,13 @@ rm -rf *.fits *.xml *.dat *.reg
 #
 # Creates pfiles directory
 # ========================
-rm -rf pfiles
-mkdir -p pfiles
-cp -r ../src/*/*.par pfiles/
 export PFILES=pfiles
-
+export CALDB=$PWD/caldb
+rm -rf $PFILES
+mkdir -p $PFILES
+cp -r ../src/*/*.par pfiles/
+rm -rf $CALDB
+mkdir -p $CALDB
  
 #
 # Test csobsdef
@@ -72,6 +74,24 @@ else
 fi
 $ECHO " ok"
  
+#
+# Test csobs2caldb
+# ================
+$ECHO -n "Test csobs2caldb: "
+chmod +x @abs_top_srcdir@/cscripts/csobs2caldb.py
+@abs_top_srcdir@/cscripts/csobs2caldb.py \
+  inobs="data/obs_unbinned.xml" \
+  rspname="test_irfs" \
+  chatter=4 #> csobs2caldb.log
+$ECHO -n "."
+if [ -s "$CALDB/data/cta/cta/bcf/test_irfs/irf_file.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " $CALDB/data/cta/cta/bcf/test_irfs/irf_file.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
  
 #
 # Test csobsinfo
@@ -279,8 +299,8 @@ chmod +x @abs_top_srcdir@/cscripts/cspull.py
   inmodel="data/crab.xml" \
   outfile="pull.dat" \
   ntrials="3" \
-  caldb="irf" \
-  irf="cta_dummy_irf" \
+  caldb="cta" \
+  irf="test_irfs" \
   ra=83.6331 \
   dec=22.0145 \
   emin=0.1 \