diff --git a/.gitignore b/.gitignore index 51710b2..75cdccf 100644 --- a/.gitignore +++ b/.gitignore @@ -63,6 +63,7 @@ src/ctbutterfly/ctbutterfly src/ctulimit/ctulimit src/cterror/cterror test/pfiles/ +test/caldb/ test/*.fits test/*.dat test/*.xml diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index f2a8d6a..a398e35 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -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) diff --git a/cscripts/csobs2caldb.py b/cscripts/csobs2caldb.py new file mode 100755 index 0000000..13b43cc --- /dev/null +++ b/cscripts/csobs2caldb.py @@ -0,0 +1,387 @@ +#!/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 . +# +# ========================================================================== +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() diff --git a/doc/source/reference_manual/csobs2caldb.rst b/doc/source/reference_manual/csobs2caldb.rst new file mode 100644 index 0000000..d97e572 --- /dev/null +++ b/doc/source/reference_manual/csobs2caldb.rst @@ -0,0 +1,72 @@ +.. _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` diff --git a/doc/source/reference_manual/index.rst b/doc/source/reference_manual/index.rst index 407b27f..86cc9ed 100644 --- a/doc/source/reference_manual/index.rst +++ b/doc/source/reference_manual/index.rst @@ -49,3 +49,4 @@ cscripts csmodelinfo --- Shows model container content csmodelmerge --- Merges several model containers into one file cstsmapmerge --- Merges slices from ts map computations + csobs2caldb --- Creates a caldb entry from an input observation diff --git a/test/data/obs_unbinned.xml b/test/data/obs_unbinned.xml new file mode 100644 index 0000000..4adfdbb --- /dev/null +++ b/test/data/obs_unbinned.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/test/irf/irf_file.fits.gz b/test/irf/irf_file.fits.gz new file mode 100644 index 0000000..e9e0aa4 --- /dev/null +++ b/test/irf/irf_file.fits.gz diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index 48691bf..a38d2ca 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -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 \