diff --git a/.gitignore b/.gitignore index f9266d6..711bd04 100644 --- a/.gitignore +++ b/.gitignore @@ -72,6 +72,7 @@ src/ctulimit/ctulimit src/cterror/cterror test/pfiles/ test/caldb/ +test/csroot2caldb_* test/iactdata_* test/*.fits test/*.dat diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 5d3aebb..34e95b1 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -5,7 +5,7 @@ pythonscripts = __init__.py \ cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \ csobsdef.py csspec.py cslightcrv.py csobsinfo.py \ csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \ - csobs2caldb.py \ + csroot2caldb.py csobs2caldb.py \ csfindobs.py csiactcopy.py csiactdata.py csiactobs.py \ obsutils.py @@ -13,14 +13,14 @@ pythonscripts = __init__.py \ parfiles = cspull.par cssens.par cstsdist.par cscaldb.par csresmap.par \ csobsdef.par csspec.par cslightcrv.par csobsinfo.par \ csmodelinfo.par csmodelmerge.par cstsmapmerge.par \ - csobs2caldb.par \ + csroot2caldb.par csobs2caldb.par \ csfindobs.par csiactcopy.par csiactdata.par csiactobs.par # Scripts for which symbolic links will be created in bindir linkscripts = cspull cssens cstsdist cscaldb csresmap \ csobsdef csspec cslightcrv csobsinfo \ csinfo csmodelinfo csmodelmerge cstsmapmerge \ - csobs2caldb \ + csroot2caldb csobs2caldb \ csfindobs csiactcopy csiactdata csiactobs # Test scripts to be distributed in $(pyexecdir)/cscripts/tests @@ -34,6 +34,7 @@ testscripts = tests/__init__.py \ ../test/test_csobsdef.py \ ../test/test_csobsinfo.py \ ../test/test_cspull.py \ + ../test/test_csroot2caldb.py \ ../test/test_csfindobs.py \ ../test/test_csiactcopy.py \ ../test/test_csiactdata.py \ @@ -53,7 +54,8 @@ testdata = ../test/data/crab.xml \ ../test/data/models.txt \ ../test/data/model_cube_background.xml \ ../test/data/model_cube_background2.xml \ - ../test/data/crab_bkgcube.xml + ../test/data/crab_bkgcube.xml \ + ../test/data/irf.root # Test data to be distributed in $(pyexecdir)/cscripts/tests/iactdata testiactdata = ../test/iactdata/hdu-index.fits.gz \ @@ -61,7 +63,7 @@ testiactdata = ../test/iactdata/hdu-index.fits.gz \ ../test/iactdata/obs-index.fits.gz # Test IRFs to be distributed in $(pyexecdir)/cscripts/tests/irf -testirf = ../test/irf/irf_file.fits.gz +testirf = ../test/irf/irf_file.fits.gz # Distribute scripts and parfiles EXTRA_DIST = $(binscripts) $(pythonscripts) diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in index 2a1064d..e4260c5 100644 --- a/cscripts/__init__.py.in +++ b/cscripts/__init__.py.in @@ -34,6 +34,7 @@ __all__ = [ "csobsinfo", "cspull", "csresmap", + "csroot2caldb", "cssens", "csspec", "cstsdist", @@ -58,6 +59,7 @@ from cscripts.csobsdef import csobsdef from cscripts.csobsinfo import csobsinfo from cscripts.cspull import cspull from cscripts.csresmap import csresmap +from cscripts.csroot2caldb import csroot2caldb from cscripts.cssens import cssens from cscripts.csspec import csspec from cscripts.cstsdist import cstsdist diff --git a/cscripts/csroot2caldb.par b/cscripts/csroot2caldb.par new file mode 100644 index 0000000..1ed70a6 --- /dev/null +++ b/cscripts/csroot2caldb.par @@ -0,0 +1,38 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csroot2caldb.par # +# Version: 1.0.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csroot2caldb script # +# # +############################################################################# +# +# Input parameters +#================= +infile, f, a, irf.root,,, "Input ROOT file" +outdir, s, h, ,,, "Output CALDB directory" +inst, s, a, prod3,,, "Instrument name" +id, s, a, South_50h,,, "Instrument Response Function identifier" +version, s, a, prod3,,, "Instrument Response Function version" +analysis, s, a, DESY,,, "Instrument Response Function analysis" +psftype, s, h, Gauss,Gauss|King,, "Point spread function type" +rebin, b, h, no,,, "Rebin" +eascale, r, h, 1.0,,, "Scaling factor for effective areas" +bgdscale, r, h, 1.0,,, "Scaling factor for background rates" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, csroot2caldb.log,,, "Log filename" diff --git a/cscripts/csroot2caldb.py b/cscripts/csroot2caldb.py new file mode 100755 index 0000000..9447a98 --- /dev/null +++ b/cscripts/csroot2caldb.py @@ -0,0 +1,1277 @@ +#! /usr/bin/env python +# ========================================================================== +# Generate IRFs in CALDB format from a ROOT offaxis performance file +# +# Copyright (C) 2016 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 . +# +# ========================================================================== +import os +import sys +import math +import gammalib +import ctools +from datetime import datetime + +# Optional ROOT import +try: + from ROOT import TFile, TH2F, TH3F + _has_root = True +except: + print('ROOT module not present, so csroot2caldb script will not work.') + _has_root = False + + +# ================== # +# csroot2caldb class # +# ================== # +class csroot2caldb(ctools.cscript): + """ + Generates IRFs in CALDB from a ROOT offaxis performance file. + """ + + # Constructor + def __init__(self, *argv): + """ + CALDB constructor. + + """ + # Set name + self._name = "csroot2caldb" + self._version = "1.0.0" + + # Initialise some members + self._cif = None + self._irf = None + self._ea = None + self._psf = None + self._edisp = None + self._bgd = None + self._hdu_cif = None + + # Initialise application by calling the appropriate class constructor + self._init_cscript(argv) + + # Return + return + + + # Private methods + def _get_parameters(self): + """ + Get parameters from parfile. + """ + # Query parameters + self["infile"].filename() + self["outdir"].string() + self["inst"].string() + self["id"].string() + self["analysis"].string() + self["version"].string() + self["psftype"].string() + self["rebin"].boolean() + self["eascale"].real() + self["bgdscale"].real() + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") + + # Return + return + + def _initialise(self): + + # Set CALDB information + self._cal_tel = "CTA" + self._cal_inst = self["inst"].string().upper() + self._cal_obsid = self["id"].string() + self._cal_det = "NONE" + self._cal_flt = "NONE" + self._cal_class = "BCF" + self._cal_type = "DATA" + self._cal_qual = 0 # 0=good, 1=bad, 2=dubious, ... + self._cal_date = "14/01/30" + self._val_date = "2014-01-30" + self._val_time = "00:00:00" + self._ref_time = 51544.0 + self._cal_name = "NAME("+self["id"].string()+")" + self._cal_version = "VERSION("+self["version"].string()+")" + self._cal_cut = "CLASS(BEST)" + self._cal_analysis = "ANALYSIS("+self["analysis"].string()+")" + # + self._ea_name = "EFF_AREA" + self._ea_doc = "CAL/GEN/92-019" + self._ea_bounds = [self._cal_name, self._cal_version, + self._cal_cut, self._cal_analysis] + self._ea_desc = "CTA effective area" + # + self._psf_name = "RPSF" + self._psf_doc = "CAL/GEN/92-020" + self._psf_bounds = [self._cal_name, self._cal_version, + self._cal_cut, self._cal_analysis] + self._psf_desc = "CTA point spread function" + # + self._edisp_name = "EDISP" + self._edisp_doc = "???" + self._edisp_bounds = [self._cal_name, self._cal_version, + self._cal_cut, self._cal_analysis] + self._edisp_desc = "CTA energy dispersion" + # + self._bgd_name = "BGD" + self._bgd_doc = "???" + self._bgd_bounds = [self._cal_name, self._cal_version, + self._cal_cut, self._cal_analysis] + self._bgd_desc = "CTA background" + + # Create directory structure + self._make_dirs() + + # Return + return + + def _make_dirs(self): + """ + Generate CALDB directory structure for one observation identifier. + The structure is given by + + data///bcf/ + + where is "cta" and is the instrument specified in the + CALDB constructor (the instrument may be used for different array + configurations). + """ + # Set base directory + base_dir = "data" + base_dir += "/"+self._cal_tel.lower() + base_dir += "/"+self._cal_inst.lower() + + # Set directory names + obs_dir = base_dir+"/bcf/"+self._cal_obsid + self._ea_dir = obs_dir + self._psf_dir = obs_dir + self._edisp_dir = obs_dir + self._bgd_dir = obs_dir + + # Optionally prefix path + outdir = self["outdir"].string() + if len(outdir) > 0: + self._base_path = outdir+"/"+base_dir + self._ea_path = outdir+"/"+self._ea_dir + self._psf_path = outdir+"/"+self._psf_dir + self._edisp_path = outdir+"/"+self._edisp_dir + self._bgd_path = outdir+"/"+self._bgd_dir + else: + self._base_path = base_dir + self._ea_path = self._ea_dir + self._psf_path = self._psf_dir + self._edisp_path = self._edisp_dir + self._bgd_path = self._bgd_dir + + # Create directory structure + if not os.path.isdir(self._ea_path): + os.makedirs(self._ea_path) + if not os.path.isdir(self._psf_path): + os.makedirs(self._psf_path) + if not os.path.isdir(self._edisp_path): + os.makedirs(self._edisp_path) + if not os.path.isdir(self._bgd_path): + os.makedirs(self._bgd_path) + + # Return + return + + def _open(self, version, split=False, clobber=True): + """ + Open existing or create new calibration. The actual version will put + all calibrations in the same file, although each part of the response + function will have its own logical name. We can thus easily modify + the script to put each calibration information in a separate file. + + Parameters: + version - Filename version + Keywords: + split - Split IRF over several files? + clobber - Overwrite existing files? + """ + # Set calibrate file names + if split: + self._ea_file = "ea_"+version+".fits" + self._psf_file = "psf_"+version+".fits" + self._edisp_file = "edisp_"+version+".fits" + self._bgd_file = "bgd_"+version+".fits" + else: + self._ea_file = "irf_"+version+".fits" + self._psf_file = "irf_"+version+".fits" + self._edisp_file = "irf_"+version+".fits" + self._bgd_file = "irf_"+version+".fits" + + # Open calibration database index + self._cif = gammalib.GFits(self._base_path+"/caldb.indx", True) + + # If file has no CIF extension than create it now + try: + self._hdu_cif = self._cif.table("CIF") + except: + self._create_cif() + self._hdu_cif = self._cif.table("CIF") + + # Set filenames + ea_filename = self._ea_path+"/"+self._ea_file + psf_filename = self._psf_path+"/"+self._psf_file + edisp_filename = self._edisp_path+"/"+self._edisp_file + bgd_filename = self._bgd_path+"/"+self._bgd_file + + # Open files + if split: + self._ea = gammalib.GFits(ea_filename, True) + self._psf = gammalib.GFits(psf_filename, True) + self._edisp = gammalib.GFits(edisp_filename, True) + self._bgd = gammalib.GFits(bgd_filename, True) + else: + self._irf = gammalib.GFits(ea_filename, True) + self._ea = self._irf + self._psf = self._irf + self._edisp = self._irf + self._bgd = self._irf + + # Open HDUs + self._hdu_ea = self._open_hdu(self._ea, "EFFECTIVE AREA", + self._ea_name, self._ea_doc) + self._hdu_psf = self._open_hdu(self._psf, "POINT SPREAD FUNCTION", + self._psf_name, self._psf_doc) + self._hdu_edisp = self._open_hdu(self._edisp, "ENERGY DISPERSION", + self._edisp_name, self._edisp_doc) + self._hdu_bgd = self._open_hdu(self._bgd, "BACKGROUND", + self._bgd_name, self._bgd_doc) + + # Return + return + + def _close(self): + """ + Close calibration FITS files. + """ + # Add information to CIF. We do this now as before we did not + # necessarily have all the information at hand (in particular + # about the boundaries) + self._add_cif_info() + + # Close CIF + self._close_file(self._cif) + + # If IRF file exists then close it now. All file pointers will + # be set to None + if self._irf != None: + self._close_file(self._irf) + self._ea = None + self._psf = None + self._edisp = None + self._bgd = None + + # ... otherwise we have split files, so we have to close them all + else: + self._close_file(self._ea) + self._close_file(self._psf) + self._close_file(self._edisp) + self._close_file(self._bgd) + + # Return + return + + def _close_file(self, file): + """ + Close one file. + """ + # If file is open then close it + if file != None: + file.save(True) + file.close() + file = None + + # Return + return + + def _open_hdu(self, fits, extname, name, doc): + """ + Open HDU + + If HDU does not exist then create one. + + Parameters: + fits: FITS file + extname: Extension name + name: Name string + doc: Document string + + Returns: + FITS table + """ + # Create table if it does not yet exist + if not fits.contains(extname): + + # Create binary table + table = gammalib.GFitsBinTable() + + # Set extension name + table.extname(extname) + + # Set OGIP keywords + self._set_ogip_keywords(table, doc, ["RESPONSE", name]) + + # Append table to FITS file + fits.append(table) + + # Return FITS table + return fits.table(extname) + + def _set_ogip_keywords(self, hdu, hdudoc, hduclas): + """ + Set standard OGIP keywords for extension. + + Parameters: + hdu: Header Data Unit. + hdudoc: Documentation reference string + hduclas: Array of HDUCLAS fields + """ + # Set UTC date of file creation + utc = datetime.utcnow().isoformat()[:19] + + # Set keywords + hdu.card("ORIGIN", "IRAP", "Name of organization making this file") + hdu.card("DATE", utc, "File creation date (YYYY-MM-DDThh:mm:ss UTC)") + hdu.card("TELESCOP", self._cal_tel, "Name of telescope") + hdu.card("INSTRUME", self._cal_inst, "Name of instrument") + hdu.card("DETNAM", self._cal_det, "Name of detector") + hdu.card("HDUCLASS", "OGIP", "HDU class") + hdu.card("HDUDOC", hdudoc, "HDU documentation") + for i, item in enumerate(hduclas): + key = "HDUCLAS%d" % (i+1) + hdu.card(key, item, "HDU class") + hdu.card("HDUVERS", "1.0.0", "HDU version") + + # Return + return + + def _create_cif(self): + """ + Create Calibration Index File (CIF) extension in CIF FITS file. + + Parameters: + None + Keywords: + None + """ + # Create binary table + table = gammalib.GFitsBinTable() + + # Set boundary + + # 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)) + + # Set keywords. Reference: CAL/GEN/92-008 + table.extname("CIF") + table.card("CIFVERSN", "1992a", "Version of CIF format") + + # Attach table to FITS file. Note that at this moment the + # FITS table gets connected to the FITS file. Yet since + # nothing was yet written to the FITS file, we cannot read + # anything from it. + self._cif.append(table) + + # Return + return + + def _add_cif_info(self): + """ + Add information to CIF extension. + """ + # Append 4 rows to CIF extension + self._hdu_cif.append_rows(4) + + # Add generic information for these 4 rows + for i in range(4): + + # Set row number + row = i + self._hdu_cif.nrows() - 4 + + # Set element + self._hdu_cif["TELESCOP"][row] = self._cal_tel + self._hdu_cif["INSTRUME"][row] = self._cal_inst + self._hdu_cif["DETNAM"][row] = self._cal_det + self._hdu_cif["FILTER"][row] = self._cal_flt + self._hdu_cif["CAL_DEV"][row] = "ONLINE" + self._hdu_cif["CAL_CLAS"][row] = self._cal_class + self._hdu_cif["CAL_DTYP"][row] = self._cal_type + self._hdu_cif["CAL_VSD"][row] = self._val_date + self._hdu_cif["CAL_VST"][row] = self._val_time + self._hdu_cif["REF_TIME"][row] = self._ref_time + self._hdu_cif["CAL_QUAL"][row] = self._cal_qual + self._hdu_cif["CAL_DATE"][row] = self._cal_date + + # Add effective area information + row = self._hdu_cif.nrows() - 4 + self._hdu_cif["CAL_DIR"][row] = self._ea_dir + self._hdu_cif["CAL_FILE"][row] = self._ea_file + self._hdu_cif["CAL_CNAM"][row] = self._ea_name + self._hdu_cif["CAL_DESC"][row] = self._ea_desc + self._hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self._ea_bounds): + self._hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self._hdu_cif["CAL_CBD"][row, n] = self._ea_bounds[n] + + # Add point spread function information + row = self._hdu_cif.nrows() - 3 + self._hdu_cif["CAL_DIR"][row] = self._psf_dir + self._hdu_cif["CAL_FILE"][row] = self._psf_file + self._hdu_cif["CAL_CNAM"][row] = self._psf_name + self._hdu_cif["CAL_DESC"][row] = self._psf_desc + self._hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self._psf_bounds): + self._hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self._hdu_cif["CAL_CBD"][row, n] = self._psf_bounds[n] + + # Add energy dispersion information + row = self._hdu_cif.nrows() - 2 + self._hdu_cif["CAL_DIR"][row] = self._edisp_dir + self._hdu_cif["CAL_FILE"][row] = self._edisp_file + self._hdu_cif["CAL_CNAM"][row] = self._edisp_name + self._hdu_cif["CAL_DESC"][row] = self._edisp_desc + self._hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self._edisp_bounds): + self._hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self._hdu_cif["CAL_CBD"][row, n] = self._edisp_bounds[n] + + # Add background information + row = self._hdu_cif.nrows() - 1 + self._hdu_cif["CAL_DIR"][row] = self._bgd_dir + self._hdu_cif["CAL_FILE"][row] = self._bgd_file + self._hdu_cif["CAL_CNAM"][row] = self._bgd_name + self._hdu_cif["CAL_DESC"][row] = self._bgd_desc + self._hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self._bgd_bounds): + self._hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self._hdu_cif["CAL_CBD"][row, n] = self._bgd_bounds[n] + + # Return + return + + def _set_cif_keywords(self, hdu, name, bounds, desc): + """ + Set standard CIF keywords for extension. + """ + # Set standard CIF keywords + hdu.card("CSYSNAME", "XMA_POL", "") + hdu.card("CCLS0001", self._cal_class, "Calibration class") + hdu.card("CDTP0001", self._cal_type, "Calibration type") + hdu.card("CCNM0001", name, "Calibration name") + + # Set boundary keywords + for n in range(9): + keyname = "CBD%d0001" % (n+1) + if n >= len(bounds): + value = "NONE" + else: + value = bounds[n] + hdu.card(keyname, value, "Calibration boundary") + + # Set validity keywords + hdu.card("CVSD0001", self._val_date, "Calibration validity start date (UTC)") + hdu.card("CVST0001", self._val_time, "Calibration validity start time (UTC)") + hdu.card("CDEC0001", desc, "Calibration description") + + # Return + return + + + # ROOT specific private members + def _root2caldb(self): + """ + Translate ROOT to CALDB information. + + Parameters: + filename - ROOT 2D performance filename. + """ + # Open ROOT performance file + file = TFile(self["infile"].filename().url()) + + # Create effective area + self._root2ea(file, rebin=self["rebin"].boolean(), + eascale=self["eascale"].real()) + + # Create point spread function + self._root2psf(file, self["psftype"].string()) + + # Create energy dispersion + self._root2edisp(file) + + # Create background + #self._root2bgd(file) + self._root2bgd3D(file, bgdscale=self["bgdscale"].real()) + + # Return + return + + def _root2ea(self, file, rebin=False, eascale=1.0): + """ + Translate ROOT to CALDB effective area extension. The following ROOT + histograms are used: + + EffectiveAreaEtrue_offaxis -> EFFAREA + EffectiveArea_offaxis -> EFFAREA_RECO + + Parameters: + file: ROOT file. + + Keywords: + rebin: Rebin Etrue histogram (useful for Prod1 IFAE runs) + eascale: Effective area scaling factor + """ + # Continue only if effective area HDU has been opened + if self._hdu_ea != None: + + # Allocate ROOT 2D arrays + etrue = TH2F() + ereco = TH2F() + file.GetObject("EffectiveAreaEtrue_offaxis", etrue) + file.GetObject("EffectiveArea_offaxis", ereco) + + # Rebin etrue histogram + if rebin: + etrue.RebinX(10) + neng = etrue.GetXaxis().GetNbins() + noffset = etrue.GetYaxis().GetNbins() + for ioff in range(noffset): + for ieng in range(neng): + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0 + etrue.SetBinContent(ieng+1,ioff+1,value) + + # Scale histograms (if needed) + if eascale != 1.0: + neng = etrue.GetXaxis().GetNbins() + noffset = etrue.GetYaxis().GetNbins() + for ioff in range(noffset): + for ieng in range(neng): + value = etrue.GetBinContent(ieng+1,ioff+1) * eascale + etrue.SetBinContent(ieng+1,ioff+1,value) + neng = ereco.GetXaxis().GetNbins() + noffset = ereco.GetYaxis().GetNbins() + for ioff in range(noffset): + for ieng in range(neng): + value = ereco.GetBinContent(ieng+1,ioff+1) * eascale + ereco.SetBinContent(ieng+1,ioff+1,value) + + # Write boundaries (use Ereco boundaries) + bounds = self._make_2D(ereco, self._hdu_ea, None, "m2") + for b in bounds: + self._ea_bounds.append(b) + self._set_cif_keywords(self._hdu_ea, self._ea_name, \ + self._ea_bounds, self._ea_desc) + + # EFFAREA + self._make_2D(etrue, self._hdu_ea, "EFFAREA", "m2") + + # EFFAREA_RECO + self._make_2D(ereco, self._hdu_ea, "EFFAREA_RECO", "m2") + + # Return + return + + def _root2psf(self, file, psftype): + """ + Translate ROOT to CALDB point spread function extension. + Parameters: + file - ROOT file. + psftype - PSF type (Gauss, King) + Keywords: + None + """ + # Continue only if point spread function HDU has been opened + if self._hdu_psf != None: + + # King profile PSF + if psftype == "King": + self._root2psf_king(file) + else: + self._root2psf_gauss(file) + + # Return + return + + def _root2psf_gauss(self, file): + """ + Translate ROOT to CALDB point spread function extension. The following + ROOT histograms are used: + + 1/(2*pi*SIGMA_1) -> SCALE + AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8) + 0.0 -> AMPL_2 + 0.0 -> SIGMA_2 + 0.0 -> AMPL_3 + 0.0 -> SIGMA_3 + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if point spread function HDU has been opened + if self._hdu_psf != None: + + # Allocate ROOT 2D array + r68 = TH2F() + + # Read 68% containment histogram + file.GetObject("AngRes_offaxis", r68) + neng = r68.GetXaxis().GetNbins() + noffset = r68.GetYaxis().GetNbins() + + # Converts 68% -> 1 sigma + r68_to_sigma = 0.6624305 + for ioff in range(noffset): + for ieng in range(neng): + sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma + r68.SetBinContent(ieng+1,ioff+1,sigma) + + # Compute scale histogram + scale = r68.Clone() + for ioff in range(noffset): + for ieng in range(neng): + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1) + if integral > 0.0: + value = 1.0 / integral + else: + value = 0.0 + scale.SetBinContent(ieng+1,ioff+1,value) + + # Set zero histogram + zero = r68.Clone() + for ioff in range(noffset): + for ieng in range(neng): + zero.SetBinContent(ieng+1,ioff+1,0.0) + + # Set boundaries + bounds = self._make_2D(r68, self._hdu_psf, None, "deg") + for b in bounds: + self._psf_bounds.append(b) + self._psf_bounds.append("PSF(GAUSS)") + self._set_cif_keywords(self._hdu_psf, self._psf_name, \ + self._psf_bounds, self._psf_desc) + + # SCALE + #self._make_2D(scale, self._hdu_psf, "SCALE", "") + + # SIGMA_1 + self._make_2D(r68, self._hdu_psf, "SIGMA_1", "deg") + + # AMPL_2 + self._make_2D(zero, self._hdu_psf, "AMPL_2", "") + + # SIGMA_2 + self._make_2D(zero, self._hdu_psf, "SIGMA_2", "deg") + + # AMPL_3 + self._make_2D(zero, self._hdu_psf, "AMPL_3", "") + + # SIGMA_3 + self._make_2D(zero, self._hdu_psf, "SIGMA_3", "deg") + + # Return + return + + def _root2psf_king(self, file): + """ + Translate ROOT to CALDB point spread function extension. The following + ROOT histograms are used: + + AngRes_offaxis + AngRes80_offaxis + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if point spread function HDU has been opened + if self._hdu_psf != None: + + # Allocate ROOT 2D arrays + r68 = TH2F() + r80 = TH2F() + + # Read 68% and 80% containment histograms + file.GetObject("AngRes_offaxis", r68) + file.GetObject("AngRes80_offaxis", r80) + neng = r68.GetXaxis().GetNbins() + noffset = r68.GetYaxis().GetNbins() + + # Initialise parameter maps by cloning the r68 2D map + gamma2D = r68.Clone() + sigma2D = r68.Clone() + + # Compute gamma and sigma values + for ioff in range(noffset): + + # Initialise last results + last_gamma = 0.0 + last_sigma = 0.0 + + for ieng in range(neng): + + # Extract radii + r_68 = r68.GetBinContent(ieng+1,ioff+1) + r_80 = r80.GetBinContent(ieng+1,ioff+1) + + # Initialise results + gamma = 0.0 + sigma = 0.0 + + # Continue only if both radii are positive + if r_68 > 0 and r_80 > 0: + + # Derive constants for equation to solve + a = 1.0 - 0.68 + b = 1.0 - 0.80 + c = r_68*r_68/(r_80*r_80) + + # Solve equation (a^x-1)/(b^x-1)=c for x using secant + # method. Stop when we are better than 1e-6. + x1 = -0.5 + x2 = -1.0 + f1 = (math.pow(a,x1) - 1.0)/(math.pow(b,x1) - 1.0) - c + f2 = (math.pow(a,x2) - 1.0)/(math.pow(b,x2) - 1.0) - c + iter = 0 + while True: + x = x1 - f1 * (x1-x2)/(f1-f2) + f = (math.pow(a,x) - 1.0)/(math.pow(b,x) - 1.0) - c + iter += 1 + if abs(f) < 1.0e-6: + break + else: + f2 = f1 + x2 = x1 + f1 = f + x1 = x + + # Compute gamma. + if x < 0.0: + gamma = 1.0 - 1.0/x + else: + gamma = 1.0 + + # Compute sigma + denominator = 2.0 * gamma * (math.pow(a, x) - 1.0) + if denominator > 0.0: + sigma = r_68 * math.sqrt(1.0/denominator) + else: + denominator = 2.0 * gamma * (math.pow(b, x) - 1.0) + if denominator > 0.0: + sigma = r_80 * math.sqrt(1.0/denominator) + else: + gamma = 0.0 + sigma = 0.0 + + # Handle special case that no results were found. + # This takes care of pixels that are ill defined + # in the MC file. + if gamma == 0.0 and sigma == 0.0: + gamma = last_gamma + sigma = last_sigma + + # Store surrent result as last result + last_gamma = gamma + last_sigma = sigma + + # Show results on console + #print(ioff,ieng,gamma,sigma,x,f,iter) + + # Set bin contents + gamma2D.SetBinContent(ieng+1,ioff+1,gamma) + sigma2D.SetBinContent(ieng+1,ioff+1,sigma) + + # Set boundaries + bounds = self._make_2D(r68, self._hdu_psf, None, "deg") + for b in bounds: + self._psf_bounds.append(b) + self._psf_bounds.append("PSF(KING)") + self._set_cif_keywords(self._hdu_psf, self._psf_name, \ + self._psf_bounds, self._psf_desc) + + # GAMMA + self._make_2D(gamma2D, self._hdu_psf, "GAMMA", "") + + # SIGMA + self._make_2D(sigma2D, self._hdu_psf, "SIGMA", "deg") + + # Return + return + + def _root2edisp(self, file): + """ + Translate ROOT to CALDB energy dispersion extension. The following ROOT + histograms are used: + + EestOverEtrue_offaxis -> MATRIX + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if energy dispersion HDU has been opened + if self._hdu_edisp != None: + + # Allocate ROOT 3D array + matrix = TH3F() + file.GetObject("EestOverEtrue_offaxis", matrix) + + # Set boundaries + bounds = self._make_3D_migra(matrix, self._hdu_edisp, None, "") + for b in bounds: + self._edisp_bounds.append(b) + self._set_cif_keywords(self._hdu_edisp, self._edisp_name, \ + self._edisp_bounds, self._edisp_desc) + + # MATRIX + self._make_3D_migra(matrix, self._hdu_edisp, "MATRIX", "") + + # Return + return + + def _root2bgd(self, file): + """ + Translate ROOT to CALDB background extension. The following ROOT + histograms are used: + + BGRatePerSqDeg_offaxis -> BGD + BGRatePerSqDeg_offaxis -> BGD_RECO + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if background HDU has been opened + if self._hdu_bgd != None: + + # Allocate ROOT 2D array + array = TH2F() + file.GetObject("BGRatePerSqDeg_offaxis", array) + + # Set boundaries + bounds = self._make_2D(array, self._hdu_bgd, None, "deg") + for b in bounds: + self._bgd_bounds.append(b) + self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \ + self._bgd_bounds, self._bgd_desc) + + # BGD + self._make_2D(array, self._hdu_bgd, "BGD", "") + + # BGD_RECO + self._make_2D(array, self._hdu_bgd, "BGD_RECO", "") + + # Return + return + + def _root2bgd3D(self, file, bgdscale=1.0): + """ + Translate ROOT to CALDB background extension. The following ROOT + histograms are used: + + BGRatePerSqDeg_offaxis -> BGD + + Parameters: + file - ROOT file. + Keywords: + scale - Background rate scaling factor + """ + # Continue only if background HDU has been opened + if self._hdu_bgd != None: + + # Allocate ROOT 2D array + array = TH2F() + file.GetObject("BGRatePerSqDeg_offaxis", array) + + # Set boundaries + bounds = self._make_3D(array, self._hdu_bgd, None, "deg") + for b in bounds: + self._bgd_bounds.append(b) + self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \ + self._bgd_bounds, self._bgd_desc) + + # BGD (reconstructed energy) + self._make_3D(array, self._hdu_bgd, "BGD", "1/s/MeV/sr", scale=bgdscale) + + # Return + return + + def _append_column(self, hdu, name, unit, axis, log=False): + """ + Append column to HDU. + + Parameters: + hdu: HDU + name: Column name + unit: Column unit + axis: ROOT histogram axis + + Keywords: + log: Logarithmic axis + """ + # Continue only if columns does not yet exist + if not hdu.contains(name): + + # Get number of axis bins + nbins = axis.GetNbins() + + # Append column and set column unit + hdu.append(gammalib.GFitsTableFloatCol(name, 1, nbins)) + hdu[name].unit(unit) + + # Do we have a LO axis? If not we have a HI axis + if name[-2:] == "LO": + low = True + else: + low = False + + # Fill column + for i in range(nbins): + if low: + value = axis.GetBinLowEdge(i+1) + else: + value = axis.GetBinUpEdge(i+1) + if log: + value = pow(10.0, value) + hdu[name][0,i] = value + + # Return + return + + def _make_2D(self, array, hdu, name, unit, scale=1.0): + """ + Make 2D matrix as function of energy and offset angle from a + ROOT 2D histogram. If the HDU has already the energy and + offset angle columns, this method will simply add another data + column. + If name==None, the method will not append any data column. + + Parameters: + array - ROOT 2D histogram. + hdu - FITS HDU. + name - Data column name. + unit - Data unit. + Keywords: + scale - Scaling factor for histogram values. + """ + # Extract energy and offset angle vectors + energies = array.GetXaxis() + offsets = array.GetYaxis() + neng = energies.GetNbins() + noffset = offsets.GetNbins() + + # Append axis columns to HDU + self._append_column(hdu, "ENERG_LO", "TeV", energies, log=True) + self._append_column(hdu, "ENERG_HI", "TeV", energies, log=True) + self._append_column(hdu, "THETA_LO", "deg", offsets) + self._append_column(hdu, "THETA_HI", "deg", offsets) + + # Append array column + if name != None and not hdu.contains(name): + hdu.append(gammalib.GFitsTableFloatCol(name, 1, neng*noffset)) + hdu[name].unit(unit) + hdu[name].dim([neng, noffset]) + for ioff in range(noffset): + for ieng in range(neng): + index = ieng + ioff * neng + value = array.GetBinContent(ieng+1,ioff+1) + hdu[name][0,index] = value * scale + + # Collect boundary information + bd_eng = "ENERG(%.4f-%.2f)TeV" % \ + (pow(10.0, energies.GetBinLowEdge(1)), \ + pow(10.0, energies.GetBinUpEdge(neng))) + bd_off = "THETA(%.2f-%.2f)deg" % \ + (offsets.GetBinLowEdge(1), \ + offsets.GetBinUpEdge(noffset)) + bd_phi = "PHI(0-360)deg" + bounds = [bd_eng, bd_off, bd_phi] + + # Return boundary information + return bounds + + def _make_3D(self, array, hdu, name, unit, scale=1.0): + """ + Make 3D cube as function of DETX, DETY and energy from a + ROOT 2D histogram. If the HDU has already the energy and + offset angle columns, this method will simply add another data + column. + If name==None, the method will not append any data column. + + Parameters: + array - ROOT 2D histogram. + hdu - FITS HDU. + name - Data column name. + unit - Data unit. + Keywords: + scale - Scaling factor for histogram values. + """ + # Set constants + deg2sr = 0.01745329*0.01745329 + + # Extract energy and offset angle vectors + energies = array.GetXaxis() + neng = energies.GetNbins() + offsets = array.GetYaxis() + noffset = offsets.GetNbins() + ewidth = [] # in MeV + for ieng in range(neng): + ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \ + pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0)) + + # Build DETX and DETY axes + ndet = array.GetYaxis().GetNbins() + ndets = 2*ndet + det_max = array.GetYaxis().GetBinUpEdge(ndet) + det_bin = det_max/float(ndet) + dets_lo = [] + dets_hi = [] + dets2 = [] + for i in range(ndets): + det_lo = -det_max + i*det_bin + det_hi = det_lo + det_bin + det_val = 0.5*(det_lo + det_hi) + dets_lo.append(det_lo) + dets_hi.append(det_hi) + dets2.append(det_val*det_val) + + # Append DETX_LO column + if not hdu.contains("DETX_LO"): + hdu.append(gammalib.GFitsTableFloatCol("DETX_LO", 1, ndets)) + hdu["DETX_LO"].unit("deg") + for i in range(ndets): + hdu["DETX_LO"][0,i] = dets_lo[i] + + # Append DETX_HI column + if not hdu.contains("DETX_HI"): + hdu.append(gammalib.GFitsTableFloatCol("DETX_HI", 1, ndets)) + hdu["DETX_HI"].unit("deg") + for i in range(ndets): + hdu["DETX_HI"][0,i] = dets_hi[i] + + # Append DETY_LO column + if not hdu.contains("DETY_LO"): + hdu.append(gammalib.GFitsTableFloatCol("DETY_LO", 1, ndets)) + hdu["DETY_LO"].unit("deg") + for i in range(ndets): + hdu["DETY_LO"][0,i] = dets_lo[i] + + # Append DETY_HI column + if not hdu.contains("DETY_HI"): + hdu.append(gammalib.GFitsTableFloatCol("DETY_HI", 1, ndets)) + hdu["DETY_HI"].unit("deg") + for i in range(ndets): + hdu["DETY_HI"][0,i] = dets_hi[i] + + # Append ENERG_LO and ENERG_HI columns + self._append_column(hdu, "ENERG_LO", "TeV", energies, log=True) + self._append_column(hdu, "ENERG_HI", "TeV", energies, log=True) + + # Append array column + if name != None and not hdu.contains(name): + hdu.append(gammalib.GFitsTableFloatCol(name, 1, ndets*ndets*neng)) + hdu[name].unit(unit) + hdu[name].dim([ndets,ndets,neng]) + for ix in range(ndets): + for iy in range(ndets): + for ieng in range(neng): + index = ix + (iy + ieng * ndets) * ndets + theta = math.sqrt(dets2[ix] + dets2[iy]) + ioff = offsets.FindBin(theta) + if ioff > 0 and ioff <= noffset: + binsq = array.GetBinContent(ieng+1,ioff) + if binsq >= 0: + value = binsq / deg2sr / (ewidth[ieng]) + hdu[name][0,index] = value * scale + + # Debugging + if False: + for ieng in range(neng): + sum = 0.0 + for ioff in range(noffset): + index = ieng + ioff * neng + binsq = array.GetBinContent(ieng+1,ioff+1) + value = binsq / deg2sr / ewidth[ieng] + sum += value + if ioff == 0: + peak = value + peaksq = binsq + print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \ + pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq) + + # Collect boundary information + bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max) + bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max) + bd_eng = "ENERG(%.4f-%.2f)TeV" % \ + (pow(10.0, energies.GetBinLowEdge(1)), + pow(10.0, energies.GetBinUpEdge(neng))) + bounds = [bd_detx, bd_dety, bd_eng] + + # Return boundary information + return bounds + + def _make_3D_migra(self, array, hdu, name, unit, scale=1.0): + """ + Make 3D cube as function of ETRUE, MIGRA and THETA form ROOT 3D + histogram. If the HDU has already the columns, this method will + simply add another data column. + If name==None, the method will not append any data column. + + Parameters: + array - ROOT 3D histogram. + hdu - FITS HDU. + name - Data column name. + unit - Data unit. + Keywords: + scale - Scaling factor for histogram values. + """ + # Extract Etrue, Eobs/Etrue and offset angle vectors + etrue = array.GetXaxis() + netrue = etrue.GetNbins() + migra = array.GetYaxis() + nmigra = migra.GetNbins() + offsets = array.GetZaxis() + noffset = offsets.GetNbins() + ewidth = [] # in MeV + for ieng in range(netrue): + ewidth.append(pow(10.0, etrue.GetBinUpEdge(ieng+1) +6.0) - \ + pow(10.0, etrue.GetBinLowEdge(ieng+1)+6.0)) + + # Append axis columns to HDU + self._append_column(hdu, "ETRUE_LO", "TeV", etrue, log=True) + self._append_column(hdu, "ETRUE_HI", "TeV", etrue, log=True) + self._append_column(hdu, "MIGRA_LO", "", migra) + self._append_column(hdu, "MIGRA_HI", "", migra) + self._append_column(hdu, "THETA_LO", "deg", offsets) + self._append_column(hdu, "THETA_HI", "deg", offsets) + + # Append migration matrix to HDU + if name != None and not hdu.contains(name): + hdu.append(gammalib.GFitsTableFloatCol(name, 1, netrue*nmigra*noffset)) + hdu[name].unit(unit) + hdu[name].dim([netrue, nmigra, noffset]) + for ioff in range(noffset): + for imigra in range(nmigra): + for ieng in range(netrue): + index = ieng + (imigra + ioff * nmigra) * netrue + value = array.GetBinContent(ieng+1,imigra+1,ioff+1) + hdu[name][0,index] = value * scale + + # Collect boundary information + bd_eng = "ETRUE(%.4f-%.2f)TeV" % \ + (pow(10.0, etrue.GetBinLowEdge(1)), + pow(10.0, etrue.GetBinUpEdge(netrue))) + bd_migra = "MIGRA(%.3f-%.3f)" % \ + (migra.GetBinLowEdge(1), + migra.GetBinUpEdge(nmigra)) + bd_off = "THETA(%.2f-%.2f)deg" % \ + (offsets.GetBinLowEdge(1), + offsets.GetBinUpEdge(noffset)) + bd_phi = "PHI(0-360)deg" + bounds = [bd_eng, bd_migra, bd_off, bd_phi] + + # Return boundary information + return bounds + + + # Public methods + 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() + + # Initialise IRF + self._initialise() + + # Open calibration file + self._open("file") + + # Translate ROOT to CALDB information + if _has_root: + self._root2caldb() + + # Close calibration file + self._close() + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Open logfile + self.logFileOpen() + + # Run the script + self.run() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Create instance of application + app = csroot2caldb(sys.argv) + + # Execute application + app.execute() diff --git a/cscripts/cta_root2caldb.py b/cscripts/cta_root2caldb.py deleted file mode 100755 index 4b88cf4..0000000 --- a/cscripts/cta_root2caldb.py +++ /dev/null @@ -1,1601 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script generates IRFs in the CALDB format of HEASARC using ROOT -# offaxis performance files. -# ------------------------------------------------------------------------- -# Copyright (C) 2011-2016 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 . -# -# ------------------------------------------------------------------------- -# Requirements: -# - gammalib -# - pyroot -# ========================================================================== -#from ROOT import TFile, TH1F, TH2F, TH3F -from ROOT import TFile, TH2F, TH3F -import gammalib -from datetime import datetime -import math -#import sys -#import glob -import os - - -# =========== # -# CalDB class # -# =========== # -class caldb(): - """ - """ - def __init__(self, inst, obsid, path=""): - """ - CALDB constructor. The constructor will create one CALDB entry for - an observation. - - Parameters: - inst - Instrument string. - obsid - Observation ID string. - Keywords: - path - CALDB Root path. - """ - # Store root path - self.path = path - - # Initialise some members - self.cif = None - self.irf = None - self.ea = None - self.psf = None - self.edisp = None - self.bgd = None - self.hdu_cif = None - - # Store UTC date of job execution - self.utc = datetime.utcnow().isoformat()[:19] - - # Set CALDB information - self.cal_tel = "CTA" - self.cal_inst = inst.upper() - self.cal_obsid = obsid - self.cal_det = "NONE" - self.cal_flt = "NONE" - self.cal_class = "BCF" - self.cal_type = "DATA" - self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ... - self.cal_date = "14/01/30" - self.val_date = "2014-01-30" - self.val_time = "00:00:00" - self.ref_time = 51544.0 - self.cal_name = "NAME("+obsid+")" - self.cal_version = "VERSION(Prod2)" - self.cal_cut = "CLASS(BEST)" - self.cal_analysis = "ANALYSIS(DESY)" - # - self.ea_name = "EFF_AREA" - self.ea_doc = "CAL/GEN/92-019" - self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.ea_desc = "CTA effective area" - # - self.psf_name = "RPSF" - self.psf_doc = "CAL/GEN/92-020" - self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.psf_desc = "CTA point spread function" - # - self.edisp_name = "EDISP" - self.edisp_doc = "???" - self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.edisp_desc = "CTA energy dispersion" - # - self.bgd_name = "BGD" - self.bgd_doc = "???" - self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.bgd_desc = "CTA background" - - # Create directory structure - self.make_dirs() - - # Return - return - - def __del__(self): - """ - Destructor. Close all FITS files. - """ - # Close calibration - self.close() - - # Return - return - - def make_dirs(self): - """ - Generate CALDB directory structure for one observation identifier. - The structure is given by - - data///bcf/ - - where is "cta" and is the instrument specified in the - CALDB constructor (the instrument may be used for different array - configurations). - - Parameters: - None - Keywords: - None - """ - # Set base path - self.base_dir = "data" - self.base_dir += "/"+self.cal_tel.lower() - self.base_dir += "/"+self.cal_inst.lower() - - # Set directory names - self.bcf_dir = self.base_dir+"/bcf" - self.obs_dir = self.bcf_dir+"/"+self.cal_obsid - self.ea_dir = self.obs_dir - self.psf_dir = self.obs_dir - self.edisp_dir = self.obs_dir - self.bgd_dir = self.obs_dir - - # Optionally prefix path - if len(self.path) > 0: - self.base_path = self.path+"/"+self.base_dir - self.bcf_path = self.path+"/"+self.bcf_dir - self.obs_path = self.path+"/"+self.obs_dir - self.ea_path = self.path+"/"+self.ea_dir - self.psf_path = self.path+"/"+self.psf_dir - self.edisp_path = self.path+"/"+self.edisp_dir - self.bgd_path = self.path+"/"+self.bgd_dir - else: - self.base_path = self.base_dir - self.bcf_path = self.bcf_dir - self.obs_path = self.obs_dir - self.ea_path = self.ea_dir - self.psf_path = self.psf_dir - self.edisp_path = self.edisp_dir - self.bgd_path = self.bgd_dir - - # Create directory structure - if not os.path.isdir(self.ea_path): - os.makedirs(self.ea_path) - if not os.path.isdir(self.psf_path): - os.makedirs(self.psf_path) - if not os.path.isdir(self.edisp_path): - os.makedirs(self.edisp_path) - if not os.path.isdir(self.bgd_path): - os.makedirs(self.bgd_path) - - # Return - return - - def open(self, version, split=False, clobber=True): - """ - Open existing or create new calibration. The actual version will put - all calibrations in the same file, although each part of the response - function will have its own logical name. We can thus easily modify - the script to put each calibration information in a separate file. - - Parameters: - version - Filename version - Keywords: - split - Split IRF over several files? - clobber - Overwrite existing files? - """ - # Set calibrate file names - if split: - self.ea_file = "ea_"+version+".fits" - self.psf_file = "psf_"+version+".fits" - self.edisp_file = "edisp_"+version+".fits" - self.bgd_file = "bgd_"+version+".fits" - else: - self.ea_file = "irf_"+version+".fits" - self.psf_file = "irf_"+version+".fits" - self.edisp_file = "irf_"+version+".fits" - self.bgd_file = "irf_"+version+".fits" - - # Open calibration database index - self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True) - - # If file has no CIF extension than create it now - try: - self.hdu_cif = self.cif.table("CIF") - except: - self.create_cif() - self.hdu_cif = self.cif.table("CIF") - - # Set filenames - ea_filename = self.ea_path+"/"+self.ea_file - psf_filename = self.psf_path+"/"+self.psf_file - edisp_filename = self.edisp_path+"/"+self.edisp_file - bgd_filename = self.bgd_path+"/"+self.bgd_file - - # Open files - if split: - self.ea = gammalib.GFits(ea_filename, True) - self.psf = gammalib.GFits(psf_filename, True) - self.edisp = gammalib.GFits(edisp_filename, True) - self.bgd = gammalib.GFits(bgd_filename, True) - else: - self.irf = gammalib.GFits(ea_filename, True) - self.ea = self.irf - self.psf = self.irf - self.edisp = self.irf - self.bgd = self.irf - - # Open HDUs - self.open_ea() - self.open_psf() - self.open_edisp() - self.open_bgd() - - # Return - return - - def close(self): - """ - Close calibration FITS files. - - Parameters: - None - Keywords: - None - """ - # Add information to CIF. We do this now as before we did not - # necessarily have all the information at hand (in particular - # about the boundaries) - self.add_cif_info() - - # Close CIF - if self.cif != None: - self.cif.save(True) - self.cif.close() - self.cif = None - - # If IRF file exists then close it now. All file pointers will - # be set to None - if self.irf != None: - self.irf.save(True) - self.irf.close() - self.irf = None - self.ea = None - self.psf = None - self.edisp = None - self.bgd = None - - # ... otherwise we have split files, so we have to close them - # all - else: - - # Close effective area file - if self.ea != None: - self.ea.save(True) - self.ea.close() - self.ea = None - - # Close point spread function file - if self.psf != None: - self.psf.save(True) - self.psf.close() - self.psf = None - - # Close energy dispersion file - if self.edisp != None: - self.edisp.save(True) - self.edisp.close() - self.edisp = None - - # Close background file - if self.bgd != None: - self.bgd.save(True) - self.bgd.close() - self.bgd = None - - # Return - return - - def create_cif(self): - """ - Create Calibration Index File (CIF) extension in CIF FITS file. - - Parameters: - None - Keywords: - None - """ - # Create binary table - table = gammalib.GFitsBinTable() - - # Set boundary - - # 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)) - - # Set keywords. Reference: CAL/GEN/92-008 - table.extname("CIF") - table.card("CIFVERSN", "1992a", "Version of CIF format") - - # Attach table to FITS file. Note that at this moment the - # FITS table gets connected to the FITS file. Yet since - # nothing was yet written to the FITS file, we cannot read - # anything from it. - self.cif.append(table) - - # Return - return - - def add_cif_info(self): - """ - Add information to CIF extension. - - Parameters: - None - Keywords: - None - """ - # Append 4 rows to CIF extension - self.hdu_cif.append_rows(4) - - # Add generic information for these 4 rows - for i in range(4): - - # Set row number - row = i+self.hdu_cif.nrows()-4 - - # Set element - self.hdu_cif["TELESCOP"][row] = self.cal_tel - self.hdu_cif["INSTRUME"][row] = self.cal_inst - self.hdu_cif["DETNAM"][row] = self.cal_det - self.hdu_cif["FILTER"][row] = self.cal_flt - self.hdu_cif["CAL_DEV"][row] = "ONLINE" - self.hdu_cif["CAL_CLAS"][row] = self.cal_class - self.hdu_cif["CAL_DTYP"][row] = self.cal_type - self.hdu_cif["CAL_VSD"][row] = self.val_date - self.hdu_cif["CAL_VST"][row] = self.val_time - self.hdu_cif["REF_TIME"][row] = self.ref_time - self.hdu_cif["CAL_QUAL"][row] = self.cal_qual - self.hdu_cif["CAL_DATE"][row] = self.cal_date - - # Add effective area information - row = self.hdu_cif.nrows()-4 - self.hdu_cif["CAL_DIR"][row] = self.ea_dir - self.hdu_cif["CAL_FILE"][row] = self.ea_file - self.hdu_cif["CAL_CNAM"][row] = self.ea_name - self.hdu_cif["CAL_DESC"][row] = self.ea_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.ea_bounds): - self.hdu_cif["CAL_CBD"][row, n] = "NONE" - else: - self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n] - - # Add point spread function information - row = self.hdu_cif.nrows()-3 - self.hdu_cif["CAL_DIR"][row] = self.psf_dir - self.hdu_cif["CAL_FILE"][row] = self.psf_file - self.hdu_cif["CAL_CNAM"][row] = self.psf_name - self.hdu_cif["CAL_DESC"][row] = self.psf_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.psf_bounds): - self.hdu_cif["CAL_CBD"][row, n] = "NONE" - else: - self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n] - - # Add energy dispersion information - row = self.hdu_cif.nrows()-2 - self.hdu_cif["CAL_DIR"][row] = self.edisp_dir - self.hdu_cif["CAL_FILE"][row] = self.edisp_file - self.hdu_cif["CAL_CNAM"][row] = self.edisp_name - self.hdu_cif["CAL_DESC"][row] = self.edisp_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.edisp_bounds): - self.hdu_cif["CAL_CBD"][row, n] = "NONE" - else: - self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n] - - # Add background information - row = self.hdu_cif.nrows()-1 - self.hdu_cif["CAL_DIR"][row] = self.bgd_dir - self.hdu_cif["CAL_FILE"][row] = self.bgd_file - self.hdu_cif["CAL_CNAM"][row] = self.bgd_name - self.hdu_cif["CAL_DESC"][row] = self.bgd_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.bgd_bounds): - self.hdu_cif["CAL_CBD"][row, n] = "NONE" - else: - self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n] - - # Return - return - - def open_ea(self): - """ - Open effective area extension. If no extension was found - then create one. We do not yet set any data as we don't - know the dimensions of the table yet. - - Parameters: - None - Keywords: - None - """ - # Get extension. If it does not exist then create it now. - try: - self.hdu_ea = self.ea.table("EFFECTIVE AREA") - except: - # Create binary table - table = gammalib.GFitsBinTable() - - # Set extension name - table.extname("EFFECTIVE AREA") - - # Set OGIP keywords - self.set_ogip_keywords(table, self.ea_doc, \ - ["RESPONSE", self.ea_name]) - - # Append table - self.ea.append(table) - - # Get extension - self.hdu_ea = self.ea.table("EFFECTIVE AREA") - - # Return - return - - def open_psf(self): - """ - Open point spread function extension. If no extension was - found then create one. We do not yet set any data as we don't - know the dimensions of the table yet. - - Parameters: - None - Keywords: - None - """ - # Get extension. If it does not exist then create it now. - try: - self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") - except: - # Create binary table - table = gammalib.GFitsBinTable() - - # Set extension name - table.extname("POINT SPREAD FUNCTION") - - # Set OGIP keywords - self.set_ogip_keywords(table, self.psf_doc, \ - ["RESPONSE", self.psf_name]) - - # Append table - self.psf.append(table) - - # Get extension - self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") - - # Return - return - - def open_edisp(self): - """ - Open energy dispersion extension. If no extension was found - then create one. We do not yet set any data as we don't - know the dimensions of the table yet. - - Parameters: - None - Keywords: - None - """ - # Get extension. If it does not exist then create it now. - try: - self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") - except: - # Create binary table - table = gammalib.GFitsBinTable() - - # Set extension name - table.extname("ENERGY DISPERSION") - - # Set OGIP keywords - self.set_ogip_keywords(table, self.edisp_doc, \ - ["RESPONSE", self.edisp_name]) - - # Append table - self.edisp.append(table) - - # Get extension - self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") - - # Return - return - - def open_bgd(self): - """ - Open background extension. If no extension was found - then create one. We do not yet set any data as we don't - know the dimensions of the table yet. - - Parameters: - None - Keywords: - None - """ - # Get extension. If it does not exist then create it now. - try: - self.hdu_bgd = self.bgd.table("BACKGROUND") - except: - # Create binary table - table = gammalib.GFitsBinTable() - - # Set extension name - table.extname("BACKGROUND") - - # Set OGIP keywords - self.set_ogip_keywords(table, self.bgd_doc, \ - ["RESPONSE", self.bgd_name]) - - # Append table - self.bgd.append(table) - - # Get extension - self.hdu_bgd = self.bgd.table("BACKGROUND") - - # Return - return - - def set_ogip_keywords(self, hdu, hdudoc, hduclas): - """ - Set standard OGIP keywords for extension. - - Parameters: - hdu - Header Data Unit. - hdudoc - Documentation reference string - hduclas - Array of HDUCLAS fields - Keywords: - None - """ - # Set keywords - hdu.card("ORIGIN", "IRAP", "Name of organization making this file") - hdu.card("DATE", self.utc, "File creation date (YYYY-MM-DDThh:mm:ss UTC)") - hdu.card("TELESCOP", self.cal_tel, "Name of telescope") - hdu.card("INSTRUME", self.cal_inst, "Name of instrument") - hdu.card("DETNAM", self.cal_det, "Name of detector") - hdu.card("HDUCLASS", "OGIP", "") - hdu.card("HDUDOC", hdudoc, "") - for i, item in enumerate(hduclas): - key = "HDUCLAS%d" % (i+1) - hdu.card(key, item, "") - hdu.card("HDUVERS", "1.0.0", "") - - # Return - return - - def make_2D(self, array, hdu, name, unit, scale=1.0): - """ - Make 2D matrix as function of energy and offset angle from a - ROOT 2D histogram. If the HDU has already the energy and - offset angle columns, this method will simply add another data - column. - If name==None, the method will not append any data column. - - Parameters: - array - ROOT 2D histogram. - hdu - FITS HDU. - name - Data column name. - unit - Data unit. - Keywords: - scale - Scaling factor for histogram values. - """ - # Extract energy and offset angle vectors - energies = array.GetXaxis() - offsets = array.GetYaxis() - neng = energies.GetNbins() - noffset = offsets.GetNbins() - - # Attach columns to HDU. We do this only if they do not yet - # exist. This allows adding columns to the same HDU in successive - # calls of this method - # - # ENERG_LO - if not hdu.contains("ENERG_LO"): - hdu.append(gammalib.GFitsTableFloatCol("ENERG_LO", 1, neng)) - hdu["ENERG_LO"].unit("TeV") - for ieng in range(neng): - e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1)) - hdu["ENERG_LO"][0,ieng] = e_lo - # - # ENERG_HI - if not hdu.contains("ENERG_HI"): - hdu.append(gammalib.GFitsTableFloatCol("ENERG_HI", 1, neng)) - hdu["ENERG_HI"].unit("TeV") - for ieng in range(neng): - e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1)) - hdu["ENERG_HI"][0,ieng] = e_hi - # - # THETA_LO - if not hdu.contains("THETA_LO"): - hdu.append(gammalib.GFitsTableFloatCol("THETA_LO", 1, noffset)) - hdu["THETA_LO"].unit("deg") - for ioff in range(noffset): - o_lo = offsets.GetBinLowEdge(ioff+1) - hdu["THETA_LO"][0,ioff] = o_lo - # - # THETA_LO - if not hdu.contains("THETA_HI"): - hdu.append(gammalib.GFitsTableFloatCol("THETA_HI", 1, noffset)) - hdu["THETA_HI"].unit("deg") - for ioff in range(noffset): - o_hi = offsets.GetBinUpEdge(ioff+1) - hdu["THETA_HI"][0,ioff] = o_hi - # - # "NAME" - if name != None and not hdu.contains(name): - hdu.append(gammalib.GFitsTableFloatCol(name, 1, neng*noffset)) - hdu[name].unit(unit) - hdu[name].dim([neng, noffset]) - for ioff in range(noffset): - for ieng in range(neng): - index = ieng + ioff * neng - value = array.GetBinContent(ieng+1,ioff+1) - hdu[name][0,index] = value * scale - - # Collect boundary information - bd_eng = "ENERG(%.4f-%.2f)TeV" % \ - (pow(10.0, energies.GetBinLowEdge(1)), \ - pow(10.0, energies.GetBinUpEdge(neng))) - bd_off = "THETA(%.2f-%.2f)deg" % \ - (offsets.GetBinLowEdge(1), \ - offsets.GetBinUpEdge(noffset)) - bd_phi = "PHI(0-360)deg" - bounds = [bd_eng, bd_off, bd_phi] - - # Return boundary information - return bounds - - def make_3D(self, array, hdu, name, unit, scale=1.0): - """ - Make 3D cube as function of DETX, DETY and energy from a - ROOT 2D histogram. If the HDU has already the energy and - offset angle columns, this method will simply add another data - column. - If name==None, the method will not append any data column. - - Parameters: - array - ROOT 2D histogram. - hdu - FITS HDU. - name - Data column name. - unit - Data unit. - Keywords: - scale - Scaling factor for histogram values. - """ - # Set constants - deg2sr = 0.01745329*0.01745329 - - # Extract energy and offset angle vectors - energies = array.GetXaxis() - neng = energies.GetNbins() - offsets = array.GetYaxis() - noffset = offsets.GetNbins() - ewidth = [] # in MeV - for ieng in range(neng): - ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \ - pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0)) - - # Build DETX and DETY axes - ndet = array.GetYaxis().GetNbins() - ndets = 2*ndet - det_max = array.GetYaxis().GetBinUpEdge(ndet) - det_bin = det_max/float(ndet) - dets_lo = [] - dets_hi = [] - dets2 = [] - for i in range(ndets): - det_lo = -det_max + i*det_bin - det_hi = det_lo + det_bin - det_val = 0.5*(det_lo + det_hi) - dets_lo.append(det_lo) - dets_hi.append(det_hi) - dets2.append(det_val*det_val) - - # Attach columns to HDU. We do this only if they do not yet - # exist. This allows adding columns to the same HDU in successive - # calls of this method - # - # DETX_LO - if not hdu.contains("DETX_LO"): - hdu.append(gammalib.GFitsTableFloatCol("DETX_LO", 1, ndets)) - hdu["DETX_LO"].unit("deg") - for i in range(ndets): - hdu["DETX_LO"][0,i] = dets_lo[i] - # - # DETX_HI - if not hdu.contains("DETX_HI"): - hdu.append(gammalib.GFitsTableFloatCol("DETX_HI", 1, ndets)) - hdu["DETX_HI"].unit("deg") - for i in range(ndets): - hdu["DETX_HI"][0,i] = dets_hi[i] - # - # DETY_LO - if not hdu.contains("DETY_LO"): - hdu.append(gammalib.GFitsTableFloatCol("DETY_LO", 1, ndets)) - hdu["DETY_LO"].unit("deg") - for i in range(ndets): - hdu["DETY_LO"][0,i] = dets_lo[i] - # - # DETY_HI - if not hdu.contains("DETY_HI"): - hdu.append(gammalib.GFitsTableFloatCol("DETY_HI", 1, ndets)) - hdu["DETY_HI"].unit("deg") - for i in range(ndets): - hdu["DETY_HI"][0,i] = dets_hi[i] - # - # ENERG_LO - if not hdu.contains("ENERG_LO"): - hdu.append(gammalib.GFitsTableFloatCol("ENERG_LO", 1, neng)) - hdu["ENERG_LO"].unit("TeV") - for ieng in range(neng): - e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1)) - hdu["ENERG_LO"][0,ieng] = e_lo - # - # ENERG_HI - if not hdu.contains("ENERG_HI"): - hdu.append(gammalib.GFitsTableFloatCol("ENERG_HI", 1, neng)) - hdu["ENERG_HI"].unit("TeV") - for ieng in range(neng): - e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1)) - hdu["ENERG_HI"][0,ieng] = e_hi - # - # "NAME" - if name != None and not hdu.contains(name): - hdu.append(gammalib.GFitsTableFloatCol(name, 1, ndets*ndets*neng)) - hdu[name].unit(unit) - hdu[name].dim([ndets,ndets,neng]) - for ix in range(ndets): - for iy in range(ndets): - for ieng in range(neng): - index = ix + (iy + ieng * ndets) * ndets - theta = math.sqrt(dets2[ix] + dets2[iy]) - ioff = offsets.FindBin(theta) - if ioff > 0 and ioff <= noffset: - binsq = array.GetBinContent(ieng+1,ioff) - if binsq >= 0: - value = binsq / deg2sr / (ewidth[ieng]) - hdu[name][0,index] = value * scale - # - # Debugging - if False: - for ieng in range(neng): - sum = 0.0 - for ioff in range(noffset): - index = ieng + ioff * neng - binsq = array.GetBinContent(ieng+1,ioff+1) - value = binsq / deg2sr / ewidth[ieng] - sum += value - if ioff == 0: - peak = value - peaksq = binsq - print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \ - pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq) - - # Collect boundary information - bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max) - bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max) - bd_eng = "ENERG(%.4f-%.2f)TeV" % \ - (pow(10.0, energies.GetBinLowEdge(1)), \ - pow(10.0, energies.GetBinUpEdge(neng))) - bounds = [bd_detx, bd_dety, bd_eng] - - # Return boundary information - return bounds - - def make_3D_migra(self, array, hdu, name, unit, scale=1.0): - """ - Make 3D cube as function of ETRUE, MIGRA and THETA form ROOT 3D - histogram. If the HDU has already the columns, this method will - simply add another data column. - If name==None, the method will not append any data column. - - Parameters: - array - ROOT 3D histogram. - hdu - FITS HDU. - name - Data column name. - unit - Data unit. - Keywords: - scale - Scaling factor for histogram values. - """ - # Extract Etrue, Eobs/Etrue and offset angle vectors - etrue = array.GetXaxis() - netrue = etrue.GetNbins() - migra = array.GetYaxis() - nmigra = migra.GetNbins() - offsets = array.GetZaxis() - noffset = offsets.GetNbins() - ewidth = [] # in MeV - for ieng in range(netrue): - ewidth.append(pow(10.0, etrue.GetBinUpEdge(ieng+1) +6.0) - \ - pow(10.0, etrue.GetBinLowEdge(ieng+1)+6.0)) - - # Attach columns to HDU. We do this only if they do not yet - # exist. This allows adding columns to the same HDU in successive - # calls of this method - # - # ETRUE_LO - if not hdu.contains("ETRUE_LO"): - hdu.append(gammalib.GFitsTableFloatCol("ETRUE_LO", 1, netrue)) - hdu["ETRUE_LO"].unit("TeV") - for ieng in range(netrue): - e_lo = pow(10.0, etrue.GetBinLowEdge(ieng+1)) - hdu["ETRUE_LO"][0,ieng] = e_lo - # - # ETRUE_HI - if not hdu.contains("ETRUE_HI"): - hdu.append(gammalib.GFitsTableFloatCol("ETRUE_HI", 1, netrue)) - hdu["ETRUE_HI"].unit("TeV") - for ieng in range(netrue): - e_hi = pow(10.0, etrue.GetBinUpEdge(ieng+1)) - hdu["ETRUE_HI"][0,ieng] = e_hi - # - # MIGRA_LO - if not hdu.contains("MIGRA_LO"): - hdu.append(gammalib.GFitsTableFloatCol("MIGRA_LO", 1, nmigra)) - hdu["MIGRA_LO"].unit("") - for imigra in range(nmigra): - migra_lo = migra.GetBinLowEdge(imigra+1) - hdu["MIGRA_LO"][0,imigra] = migra_lo - # - # MIGRA_HI - if not hdu.contains("MIGRA_HI"): - hdu.append(gammalib.GFitsTableFloatCol("MIGRA_HI", 1, nmigra)) - hdu["MIGRA_HI"].unit("") - for imigra in range(nmigra): - migra_hi = migra.GetBinUpEdge(imigra+1) - hdu["MIGRA_HI"][0,imigra] = migra_hi - # - # THETA_LO - if not hdu.contains("THETA_LO"): - hdu.append(gammalib.GFitsTableFloatCol("THETA_LO", 1, noffset)) - hdu["THETA_LO"].unit("deg") - for ioff in range(noffset): - o_lo = offsets.GetBinLowEdge(ioff+1) - hdu["THETA_LO"][0,ioff] = o_lo - # - # THETA_HI - if not hdu.contains("THETA_HI"): - hdu.append(gammalib.GFitsTableFloatCol("THETA_HI", 1, noffset)) - hdu["THETA_HI"].unit("deg") - for ioff in range(noffset): - o_hi = offsets.GetBinUpEdge(ioff+1) - hdu["THETA_HI"][0,ioff] = o_hi - # - # "NAME" - if name != None and not hdu.contains(name): - hdu.append(gammalib.GFitsTableFloatCol(name, 1, netrue*nmigra*noffset)) - hdu[name].unit(unit) - hdu[name].dim([netrue, nmigra, noffset]) - for ioff in range(noffset): - for imigra in range(nmigra): - for ieng in range(netrue): - index = ieng + (imigra + ioff * nmigra) * netrue - value = array.GetBinContent(ieng+1,imigra+1,ioff+1) - hdu[name][0,index] = value * scale - - # Collect boundary information - bd_eng = "ETRUE(%.4f-%.2f)TeV" % \ - (pow(10.0, etrue.GetBinLowEdge(1)), \ - pow(10.0, etrue.GetBinUpEdge(netrue))) - bd_migra = "MIGRA(%.4f-%.2f)TeV" % \ - (pow(10.0, migra.GetBinLowEdge(1)), \ - pow(10.0, migra.GetBinUpEdge(nmigra))) - bd_off = "THETA(%.2f-%.2f)deg" % \ - (offsets.GetBinLowEdge(1), \ - offsets.GetBinUpEdge(noffset)) - bd_phi = "PHI(0-360)deg" - bounds = [bd_eng, bd_migra, bd_off, bd_phi] - - # Return boundary information - return bounds - - def root2caldb(self, filename, rebin=False, psftype="Gauss", eascale=1.0, bgdscale=1.0): - """ - Translate ROOT to CALDB information. - - Parameters: - filename - ROOT 2D performance filename. - Keywords: - rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs) - psftype - PSF type - eascale - Effective area scaling factor - bgdscale - Background scaling factor - """ - # Open ROOT performance file - file = TFile(filename) - - # Create effective area - self.root2ea(file, rebin=rebin, eascale=eascale) - - # Create point spread function - self.root2psf(file, psftype) - - # Create energy dispersion - self.root2edisp_migra(file) - - # Create background - #self.root2bgd(file) - self.root2bgd3D(file, bgdscale=bgdscale) - - # Return - return - - def root2ea(self, file, rebin=False, eascale=1.0): - """ - Translate ROOT to CALDB effective area extension. The following ROOT - histograms are used: - - EffectiveAreaEtrue_offaxis -> EFFAREA - EffectiveArea_offaxis -> EFFAREA_RECO - - Parameters: - file - ROOT file. - Keywords: - rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs) - eascale - Effective area scaling factor - """ - # Continue only if effective area HDU has been opened - if self.hdu_ea != None: - - # Allocate ROOT 2D arrays - etrue = TH2F() - ereco = TH2F() - file.GetObject("EffectiveAreaEtrue_offaxis", etrue) - file.GetObject("EffectiveArea_offaxis", ereco) - - # Rebin etrue histogram - if rebin: - etrue.RebinX(10) - neng = etrue.GetXaxis().GetNbins() - noffset = etrue.GetYaxis().GetNbins() - for ioff in range(noffset): - for ieng in range(neng): - value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0 - etrue.SetBinContent(ieng+1,ioff+1,value) - - # Scale histograms (if needed) - if eascale != 1.0: - neng = etrue.GetXaxis().GetNbins() - noffset = etrue.GetYaxis().GetNbins() - for ioff in range(noffset): - for ieng in range(neng): - value = etrue.GetBinContent(ieng+1,ioff+1) * eascale - etrue.SetBinContent(ieng+1,ioff+1,value) - neng = ereco.GetXaxis().GetNbins() - noffset = ereco.GetYaxis().GetNbins() - for ioff in range(noffset): - for ieng in range(neng): - value = ereco.GetBinContent(ieng+1,ioff+1) * eascale - ereco.SetBinContent(ieng+1,ioff+1,value) - - # Write boundaries (use Ereco boundaries) - bounds = self.make_2D(ereco, self.hdu_ea, None, "m2") - for b in bounds: - self.ea_bounds.append(b) - self.set_cif_keywords(self.hdu_ea, self.ea_name, \ - self.ea_bounds, self.ea_desc) - - # EFFAREA - self.make_2D(etrue, self.hdu_ea, "EFFAREA", "m2") - - # EFFAREA_RECO - self.make_2D(ereco, self.hdu_ea, "EFFAREA_RECO", "m2") - - # Return - return - - def root2psf(self, file, psftype): - """ - Translate ROOT to CALDB point spread function extension. - Parameters: - file - ROOT file. - psftype - PSF type (Gauss, King) - Keywords: - None - """ - # Continue only if point spread function HDU has been opened - if self.hdu_psf != None: - - # King profile PSF - if psftype == "King": - self.root2psf_king(file) - else: - self.root2psf_gauss(file) - - # Return - return - - def root2psf_gauss(self, file): - """ - Translate ROOT to CALDB point spread function extension. The following - ROOT histograms are used: - - 1/(2*pi*SIGMA_1) -> SCALE - AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8) - 0.0 -> AMPL_2 - 0.0 -> SIGMA_2 - 0.0 -> AMPL_3 - 0.0 -> SIGMA_3 - - Parameters: - file - ROOT file. - Keywords: - None - """ - # Continue only if point spread function HDU has been opened - if self.hdu_psf != None: - - # Allocate ROOT 2D array - r68 = TH2F() - - # Read 68% containment histogram - file.GetObject("AngRes_offaxis", r68) - neng = r68.GetXaxis().GetNbins() - noffset = r68.GetYaxis().GetNbins() - - # Converts 68% -> 1 sigma - r68_to_sigma = 0.6624305 - for ioff in range(noffset): - for ieng in range(neng): - sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma - r68.SetBinContent(ieng+1,ioff+1,sigma) - - # Compute scale histogram - scale = r68.Clone() - for ioff in range(noffset): - for ieng in range(neng): - integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1) - if integral > 0.0: - value = 1.0 / integral - else: - value = 0.0 - scale.SetBinContent(ieng+1,ioff+1,value) - - # Set zero histogram - zero = r68.Clone() - for ioff in range(noffset): - for ieng in range(neng): - zero.SetBinContent(ieng+1,ioff+1,0.0) - - # Set boundaries - bounds = self.make_2D(r68, self.hdu_psf, None, "deg") - for b in bounds: - self.psf_bounds.append(b) - self.psf_bounds.append("PSF(GAUSS)") - self.set_cif_keywords(self.hdu_psf, self.psf_name, \ - self.psf_bounds, self.psf_desc) - - # SCALE - self.make_2D(scale, self.hdu_psf, "SCALE", "") - - # SIGMA_1 - self.make_2D(r68, self.hdu_psf, "SIGMA_1", "deg") - - # AMPL_2 - self.make_2D(zero, self.hdu_psf, "AMPL_2", "") - - # SIGMA_2 - self.make_2D(zero, self.hdu_psf, "SIGMA_2", "deg") - - # AMPL_3 - self.make_2D(zero, self.hdu_psf, "AMPL_3", "") - - # SIGMA_3 - self.make_2D(zero, self.hdu_psf, "SIGMA_3", "deg") - - # Return - return - - def root2psf_king(self, file): - """ - Translate ROOT to CALDB point spread function extension. The following - ROOT histograms are used: - - AngRes_offaxis - AngRes80_offaxis - - Parameters: - file - ROOT file. - Keywords: - None - """ - # Continue only if point spread function HDU has been opened - if self.hdu_psf != None: - - # Allocate ROOT 2D arrays - r68 = TH2F() - r80 = TH2F() - - # Read 68% and 80% containment histograms - file.GetObject("AngRes_offaxis", r68) - file.GetObject("AngRes80_offaxis", r80) - neng = r68.GetXaxis().GetNbins() - noffset = r68.GetYaxis().GetNbins() - - # Initialise parameter maps by cloning the r68 2D map - gamma2D = r68.Clone() - sigma2D = r68.Clone() - - # Compute gamma and sigma values - for ioff in range(noffset): - - # Initialise last results - last_gamma = 0.0 - last_sigma = 0.0 - - for ieng in range(neng): - - # Extract radii - r_68 = r68.GetBinContent(ieng+1,ioff+1) - r_80 = r80.GetBinContent(ieng+1,ioff+1) - - # Initialise results - gamma = 0.0 - sigma = 0.0 - - # Continue only if both radii are positive - if r_68 > 0 and r_80 > 0: - - # Derive constants for equation to solve - a = 1.0 - 0.68 - b = 1.0 - 0.80 - c = r_68*r_68/(r_80*r_80) - - # Solve equation (a^x-1)/(b^x-1)=c for x using secant - # method. Stop when we are better than 1e-6. - x1 = -0.5 - x2 = -1.0 - f1 = (math.pow(a,x1) - 1.0)/(math.pow(b,x1) - 1.0) - c - f2 = (math.pow(a,x2) - 1.0)/(math.pow(b,x2) - 1.0) - c - iter = 0 - while True: - x = x1 - f1 * (x1-x2)/(f1-f2) - f = (math.pow(a,x) - 1.0)/(math.pow(b,x) - 1.0) - c - iter += 1 - if abs(f) < 1.0e-6: - break - else: - f2 = f1 - x2 = x1 - f1 = f - x1 = x - - # Compute gamma. - if x < 0.0: - gamma = 1.0 - 1.0/x - else: - gamma = 1.0 - - # Compute sigma - denominator = 2.0 * gamma * (math.pow(a, x) - 1.0) - if denominator > 0.0: - sigma = r_68 * math.sqrt(1.0/denominator) - else: - denominator = 2.0 * gamma * (math.pow(b, x) - 1.0) - if denominator > 0.0: - sigma = r_80 * math.sqrt(1.0/denominator) - else: - gamma = 0.0 - sigma = 0.0 - - # Handle special case that no results were found. - # This takes care of pixels that are ill defined - # in the MC file. - if gamma == 0.0 and sigma == 0.0: - gamma = last_gamma - sigma = last_sigma - - # Store surrent result as last result - last_gamma = gamma - last_sigma = sigma - - # Show results on console - #print(ioff,ieng,gamma,sigma,x,f,iter) - - # Set bin contents - gamma2D.SetBinContent(ieng+1,ioff+1,gamma) - sigma2D.SetBinContent(ieng+1,ioff+1,sigma) - - # Set boundaries - bounds = self.make_2D(r68, self.hdu_psf, None, "deg") - for b in bounds: - self.psf_bounds.append(b) - self.psf_bounds.append("PSF(KING)") - self.set_cif_keywords(self.hdu_psf, self.psf_name, \ - self.psf_bounds, self.psf_desc) - - # GAMMA - self.make_2D(gamma2D, self.hdu_psf, "GAMMA", "") - - # SIGMA - self.make_2D(sigma2D, self.hdu_psf, "SIGMA", "deg") - - # Return - return - - def root2edisp_2D(self, file): - """ - Translate ROOT to CALDB energy dispersion extension. The following ROOT - histograms are used: - - ERes_offaxis -> ERES - Ebias_offaxis -> EBIAS - - Parameters: - file - ROOT file. - Keywords: - None - """ - # Continue only if energy dispersion HDU has been opened - if self.hdu_edisp != None: - - # Allocate ROOT 2D array - eres = TH2F() - ebias = TH2F() - file.GetObject("ERes_offaxis", eres) - file.GetObject("Ebias_offaxis", ebias) - - # Set boundaries - bounds = self.make_2D(eres, self.hdu_edisp, None, "deg") - for b in bounds: - self.edisp_bounds.append(b) - self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \ - self.edisp_bounds, self.edisp_desc) - - # ERES - self.make_2D(eres, self.hdu_edisp, "ERES", "Ereco/Etrue") - - # EBIAS - self.make_2D(ebias, self.hdu_edisp, "EBIAS", "Ereco/Etrue") - - # Return - return - - def root2edisp_migra(self, file): - """ - Translate ROOT to CALDB energy dispersion extension. The following ROOT - histograms are used: - - EestOverEtrue_offaxis -> MATRIX - - Parameters: - file - ROOT file. - Keywords: - None - """ - # Continue only if energy dispersion HDU has been opened - if self.hdu_edisp != None: - - # Allocate ROOT 3D array - matrix = TH3F() - file.GetObject("EestOverEtrue_offaxis", matrix) - - # Set boundaries - bounds = self.make_3D_migra(matrix, self.hdu_edisp, None, "") - for b in bounds: - self.edisp_bounds.append(b) - self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \ - self.edisp_bounds, self.edisp_desc) - - # MATRIX - self.make_3D_migra(matrix, self.hdu_edisp, "MATRIX", "") - - # Return - return - - def root2bgd(self, file): - """ - Translate ROOT to CALDB background extension. The following ROOT - histograms are used: - - BGRatePerSqDeg_offaxis -> BGD - BGRatePerSqDeg_offaxis -> BGD_RECO - - Parameters: - file - ROOT file. - Keywords: - None - """ - # Continue only if background HDU has been opened - if self.hdu_bgd != None: - - # Allocate ROOT 2D array - array = TH2F() - file.GetObject("BGRatePerSqDeg_offaxis", array) - - # Set boundaries - bounds = self.make_2D(array, self.hdu_bgd, None, "deg") - for b in bounds: - self.bgd_bounds.append(b) - self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \ - self.bgd_bounds, self.bgd_desc) - - # BGD - self.make_2D(array, self.hdu_bgd, "BGD", "") - - # BGD_RECO - self.make_2D(array, self.hdu_bgd, "BGD_RECO", "") - - # Return - return - - def root2bgd3D(self, file, bgdscale=1.0): - """ - Translate ROOT to CALDB background extension. The following ROOT - histograms are used: - - BGRatePerSqDeg_offaxis -> BGD - - Parameters: - file - ROOT file. - Keywords: - scale - Background rate scaling factor - """ - # Continue only if background HDU has been opened - if self.hdu_bgd != None: - - # Allocate ROOT 2D array - array = TH2F() - file.GetObject("BGRatePerSqDeg_offaxis", array) - - # Set boundaries - bounds = self.make_3D(array, self.hdu_bgd, None, "deg") - for b in bounds: - self.bgd_bounds.append(b) - self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \ - self.bgd_bounds, self.bgd_desc) - - # BGD (reconstructed energy) - self.make_3D(array, self.hdu_bgd, "BGD", "MeV^-1 s^-1 sr^-1", scale=bgdscale) - - # Return - return - - def set_cif_keywords(self, hdu, name, bounds, desc): - """ - Set standard CIF keywords for extension. - """ - # Set standard CIF keywords - hdu.card("CSYSNAME", "XMA_POL", "") - hdu.card("CCLS0001", self.cal_class, "Calibration class") - hdu.card("CDTP0001", self.cal_type, "Calibration type") - hdu.card("CCNM0001", name, "Calibration name") - - # Set boundary keywords - for n in range(9): - keyname = "CBD%d0001" % (n+1) - if n >= len(bounds): - value = "NONE" - else: - value = bounds[n] - hdu.card(keyname, value, "Calibration boundary") - - # Set validity keywords - hdu.card("CVSD0001", self.val_date, "Calibration validity start date (UTC)") - hdu.card("CVST0001", self.val_time, "Calibration validity start time (UTC)") - hdu.card("CDEC0001", desc, "Calibration description") - - # Return - return - - -# ================= # -# Set test database # -# ================= # -def set_test(): - """ - Set Prod1 IFAE database. - """ - # Set database attributes - path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012" - rebin = True - - # Set database content - db = [{'inst': "e", 'id': "IFAE20120510_50h", - 'path': path, 'rebin': rebin, 'psftype': "Gauss", - 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}, - {'inst': "e", 'id': "IFAE20120510_50h_King", - 'path': path, 'rebin': rebin, 'psftype': "King", - 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}] - - # Return database - return db - - -# ======================= # -# Set Prod1 IFAE database # -# ======================= # -def set_prod1_ifae(): - """ - Set Prod1 IFAE database. - """ - # Set database attributes - path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012" - rebin = True - psftype = "Gauss" - - # Set database content - db = [{'inst': "b", 'id': "IFAE20120510_50h", - 'path': path, 'rebin': rebin, 'psftype': psftype, - 'file': "SubarrayB_IFAE_50hours_20120510_offaxis.root"}, - {'inst': "e", 'id': "IFAE20120510_50h", - 'path': path, 'rebin': rebin, 'psftype': psftype, - 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}, - {'inst': "i", 'id': "IFAE20120510_50h", - 'path': path, 'rebin': rebin, 'psftype': psftype, - 'file': "SubarrayI_IFAE_50hours_20120510_offaxis.root"}] - - # Return database - return db - - -# ======================= # -# Set Prod2 DESY database # -# ======================= # -def set_prod2_desy(): - """ - Set Prod2 DESY database. - """ - # Set database attributes - path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar" - path_ten = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Tenerife_50h" - rebin = False - psftype = "Gauss" - eascale = 1.0 / 0.68 # Prod 2 DESY effective areas are for 68% containment - - # Set database content - db = [{'inst': "aar", 'id': "DESY20140105_50h", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_0deg", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_180deg", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_noLST", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a-noLST.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_0deg_noLST", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a-noLST.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_180deg_noLST", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a-noLST.180000s.root"}, - - {'inst': "aar500", 'id': "DESY20140105_50h", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"}, - {'inst': "aar500", 'id': "DESY20140105_50h_0deg", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"}, - {'inst': "aar500", 'id': "DESY20140105_50h_180deg", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_noLST", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a-noLST.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_0deg_noLST", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a-noLST.180000s.root"}, - {'inst': "aar", 'id': "DESY20140105_50h_180deg_noLST", - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a-noLST.180000s.root"}, - - {'inst': "tenerife", 'id': "DESY20140105_50h", - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}, - {'inst': "tenerife", 'id': "DESY20140105_50h_0deg", - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}, - {'inst': "tenerife", 'id': "DESY20140105_50h_180deg", - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}, - {'inst': "tenerife", 'id': "DESY20140105_50h_noLST", - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN-MST.180000s.root"}, - {'inst': "tenerife", 'id': "DESY20140105_50h_0deg_noLST", - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN-MST.180000s.root"}, - {'inst': "tenerife", 'id': "DESY20140105_50h_180deg_noLST", - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale, - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN-MST.180000s.root"} - ] - - # Return database - return db - - -# ======================== # -# Main routine entry point # -# ======================== # -if __name__ == '__main__': - """ - This script translates a CTA ROOT 2D performance files into a CALDB - compliant CTA response. The script will create the response for 3 - observations to illustrate the run-wise file structure. - - The CTA ROOT 2D performance files contain the following histograms: - - DiffSens_offaxis (2D) - Differential sensitivity - - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?) - - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment - - AngRes_offaxis (2D) - Angular resolution 68% containment - - AngRes80_offaxis (2D) - Angular resolution 80% containment - - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree [counts/s/deg2] - - BGRate_offaxis (2D) - Background rate - - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins) - - MigMatrix_offaxis (3D) - Migration matrix - - EestOverEtrue_offaxis (3D) - Migration matrix - - ERes_offaxis (2D) - Energy Resolution - - Ebias_offaxis (2D) - Energy bias - """ - # Get database - #entries = set_test() - #entries = set_prod1_ifae() - entries = set_prod2_desy() - - # Loop over entries - for entry in entries: - - # Get attributes - inst = entry['inst'] - path = entry['path'] - id = entry['id'] - rebin = entry['rebin'] - psftype = entry['psftype'] - if entry.has_key('eascale'): - eascale = entry['eascale'] - else: - eascale = 1.0 - if entry.has_key('bgdscale'): - bgdscale = entry['bgdscale'] - else: - bgdscale = 1.0 - filename = path+"/"+entry['file'] - - # Allocate caldb entry - irf = caldb(inst, id) - - # Open calibration file - irf.open("file") - - # Translate ROOT to CALDB information - irf.root2caldb(filename, rebin=rebin, psftype=psftype, eascale=eascale, bgdscale=bgdscale) diff --git a/examples/generate_prod3_irfs.py b/examples/generate_prod3_irfs.py new file mode 100755 index 0000000..fbb263d --- /dev/null +++ b/examples/generate_prod3_irfs.py @@ -0,0 +1,88 @@ +#! /usr/bin/env python +# ========================================================================== +# Generates Prod3 IRFs from ROOT performance files. +# +# Copyright (C) 2016 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 . +# +# ========================================================================== +import gammalib +import cscripts + + +# ======================= # +# Set Prod3 DESY database # +# ======================= # +def set_prod3_desy(): + """ + Set Prod3 DESY database. + """ + # Set database attributes + path = "/project-data/cta/performance/prod3" + psftype = "Gauss" + + # Set database content + db = [{'inst': "prod3", 'id': "South_50h", 'psftype': psftype, 'path': path, + 'file': "DESY.d20160304.V5.ID0NIM2.prod3-paranalp05-NN.S.3HB1-ND-3.180000s.root"} + ] + + # Return database + return db + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Get database + entries = set_prod3_desy() + + # Loop over entries + for entry in entries: + + # Set optional attributes + if entry.has_key('rebin'): + rebin = entry['rebin'] + else: + rebin = False + if entry.has_key('eascale'): + eascale = entry['eascale'] + else: + eascale = 1.0 + if entry.has_key('bgdscale'): + bgdscale = entry['bgdscale'] + else: + bgdscale = 1.0 + + # Set filename + filename = entry['path']+"/"+entry['file'] + + # Set-up csroot2caldb + caldb = cscripts.csroot2caldb() + caldb["infile"] = filename + caldb["outdir"] = "prod3" + caldb["inst"] = entry['inst'] + caldb["id"] = entry['id'] + caldb["version"] = entry['inst'] + caldb["analysis"] = "DESY" + caldb["psftype"] = entry['psftype'] + caldb["rebin"] = rebin + caldb["eascale"] = eascale + caldb["bgdscale"] = bgdscale + caldb["logfile"] = "generate_prod3_irfs.log" + + # Add CALDB entry + caldb.execute() diff --git a/test/data/irf.root b/test/data/irf.root new file mode 100644 index 0000000..9364218 --- /dev/null +++ b/test/data/irf.root diff --git a/test/test_csroot2caldb.py b/test/test_csroot2caldb.py new file mode 100644 index 0000000..004aae0 --- /dev/null +++ b/test/test_csroot2caldb.py @@ -0,0 +1,118 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the csroot2caldb script. +# +# Copyright (C) 2016 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 . +# +# ========================================================================== +import os +import gammalib +import cscripts + + +# ================================== # +# Test class for csroot2caldb script # +# ================================== # +class Test(gammalib.GPythonTestSuite): + """ + Test class for csroot2caldb script. + + This test class makes unit tests for the cslightcrv script by using it + from the command line and from Python. + """ + + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + gammalib.GPythonTestSuite.__init__(self) + + # Clean calibration database + os.system("rm -rf csroot2caldb_* >/dev/null 2>&1") + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name("csroot2caldb") + + # Append tests + self.append(self._test_cmd, "Test csroot2caldb on command line") + self.append(self._test_python, "Test csroot2caldb from Python") + + # Return + return + + # Test csroot2caldb on command line + def _test_cmd(self): + """ + Test csroot2caldb on the command line. + """ + # Kluge to set the command (installed version has no README file) + if os.path.isfile("README"): + csroot2caldb = "../cscripts/csroot2caldb.py" + else: + csroot2caldb = "csroot2caldb" + + # Setup csroot2caldb command + cmd = csroot2caldb+' infile="data/irf.root"'+ \ + ' outdir="csroot2caldb_cmd1"'+ \ + ' inst=prod3 id=South_50h version=prod3'+ \ + ' analysis=DESY'+ \ + ' logfile="csroot2caldb_cmd1.log" chatter=1' + + # Execute csroot2caldb, make sure we catch any exception + try: + rc = os.system(cmd+" >/dev/null 2>&1") + except: + pass + + # Check if execution was successful + self.test_assert(rc == 0, + "Successful csroot2caldb execution on command line") + + # Return + return + + # Test csroot2caldb from Python + def _test_python(self): + """ + Test csroot2caldb from Python. + """ + # Set-up csroot2caldb + caldb = cscripts.csroot2caldb() + caldb["infile"] = "data/irf.root" + caldb["outdir"] = "csroot2caldb_py1" + caldb["inst"] = "prod3" + caldb["id"] = "South_50h" + caldb["version"] = "prod3" + caldb["analysis"] = "DESY" + caldb["logfile"] = "csroot2caldb_py1.log" + caldb["chatter"] = 2 + + # Run script + caldb.logFileOpen() # Make sure we get a log file + caldb.run() + + # Return + return diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py index 415ca2e..6ac74d8 100755 --- a/test/test_python_cscripts.py +++ b/test/test_python_cscripts.py @@ -31,6 +31,7 @@ import test_csobs2caldb import test_csobsdef import test_csobsinfo import test_cspull +import test_csroot2caldb import test_csfindobs import test_csiactcopy import test_csiactdata @@ -92,6 +93,7 @@ def test(installed=False, debug=False): suite_csobsdef = test_csobsdef.Test() suite_csobsinfo = test_csobsinfo.Test() suite_cspull = test_cspull.Test() + suite_csroot2caldb = test_csroot2caldb.Test() # Setup unit tests suite_cscaldb.set() @@ -102,6 +104,7 @@ def test(installed=False, debug=False): suite_csobsdef.set() suite_csobsinfo.set() suite_cspull.set() + suite_csroot2caldb.set() # Append tests to container suites.append(suite_cscaldb) @@ -112,6 +115,7 @@ def test(installed=False, debug=False): suites.append(suite_csobsdef) suites.append(suite_csobsinfo) suites.append(suite_cspull) + suites.append(suite_csroot2caldb) # Append tests for Python 2.6+ (the IACT cscripts depend on the json # module which is only available since Python 2.6+