Michael Mayer / ctools

Commit 0a8e572995a9eb22174ddf61abaf3ed741e44a3e

Authored by Jürgen Knödlseder
1 parent 7c3dedd7

Replace cta_root2caldb.py script by csroot2caldb cscript. 

Convert the cta_root2caldb.py script into a cscript and add unit test. The script only works if PyROOT is installed.

Instead of the cta_root2caldb.py script a generate_prod3_irf.py script was added to the example folder.
Inline Side-by-side
Showing 10 changed files with 1535 additions and 1606 deletions
.gitignore
  View file @0a8e572
... ... @@ -72,6 +72,7 @@ src/ctulimit/ctulimit
72 72 src/cterror/cterror
73 73 test/pfiles/
74 74 test/caldb/
  75 +test/csroot2caldb_*
75 76 test/iactdata_*
76 77 test/*.fits
77 78 test/*.dat
... ...
cscripts/Makefile.am
  View file @0a8e572
... ... @@ -5,7 +5,7 @@ pythonscripts = __init__.py \
5 5 cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
6 6 csobsdef.py csspec.py cslightcrv.py csobsinfo.py \
7 7 csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \
8   - csobs2caldb.py \
  8 + csroot2caldb.py csobs2caldb.py \
9 9 csfindobs.py csiactcopy.py csiactdata.py csiactobs.py \
10 10 obsutils.py
11 11  
... ... @@ -13,14 +13,14 @@ pythonscripts = __init__.py \
13 13 parfiles = cspull.par cssens.par cstsdist.par cscaldb.par csresmap.par \
14 14 csobsdef.par csspec.par cslightcrv.par csobsinfo.par \
15 15 csmodelinfo.par csmodelmerge.par cstsmapmerge.par \
16   - csobs2caldb.par \
  16 + csroot2caldb.par csobs2caldb.par \
17 17 csfindobs.par csiactcopy.par csiactdata.par csiactobs.par
18 18  
19 19 # Scripts for which symbolic links will be created in bindir
20 20 linkscripts = cspull cssens cstsdist cscaldb csresmap \
21 21 csobsdef csspec cslightcrv csobsinfo \
22 22 csinfo csmodelinfo csmodelmerge cstsmapmerge \
23   - csobs2caldb \
  23 + csroot2caldb csobs2caldb \
24 24 csfindobs csiactcopy csiactdata csiactobs
25 25  
26 26 # Test scripts to be distributed in $(pyexecdir)/cscripts/tests
... ... @@ -34,6 +34,7 @@ testscripts = tests/__init__.py \
34 34 ../test/test_csobsdef.py \
35 35 ../test/test_csobsinfo.py \
36 36 ../test/test_cspull.py \
  37 + ../test/test_csroot2caldb.py \
37 38 ../test/test_csfindobs.py \
38 39 ../test/test_csiactcopy.py \
39 40 ../test/test_csiactdata.py \
... ... @@ -53,7 +54,8 @@ testdata = ../test/data/crab.xml \
53 54 ../test/data/models.txt \
54 55 ../test/data/model_cube_background.xml \
55 56 ../test/data/model_cube_background2.xml \
56   - ../test/data/crab_bkgcube.xml
  57 + ../test/data/crab_bkgcube.xml \
  58 + ../test/data/irf.root
57 59  
58 60 # Test data to be distributed in $(pyexecdir)/cscripts/tests/iactdata
59 61 testiactdata = ../test/iactdata/hdu-index.fits.gz \
... ... @@ -61,7 +63,7 @@ testiactdata = ../test/iactdata/hdu-index.fits.gz \
61 63 ../test/iactdata/obs-index.fits.gz
62 64  
63 65 # Test IRFs to be distributed in $(pyexecdir)/cscripts/tests/irf
64   -testirf = ../test/irf/irf_file.fits.gz
  66 +testirf = ../test/irf/irf_file.fits.gz
65 67  
66 68 # Distribute scripts and parfiles
67 69 EXTRA_DIST = $(binscripts) $(pythonscripts)
... ...
cscripts/__init__.py.in
  View file @0a8e572
... ... @@ -34,6 +34,7 @@ __all__ = [
34 34 "csobsinfo",
35 35 "cspull",
36 36 "csresmap",
  37 + "csroot2caldb",
37 38 "cssens",
38 39 "csspec",
39 40 "cstsdist",
... ... @@ -58,6 +59,7 @@ from cscripts.csobsdef import csobsdef
58 59 from cscripts.csobsinfo import csobsinfo
59 60 from cscripts.cspull import cspull
60 61 from cscripts.csresmap import csresmap
  62 +from cscripts.csroot2caldb import csroot2caldb
61 63 from cscripts.cssens import cssens
62 64 from cscripts.csspec import csspec
63 65 from cscripts.cstsdist import cstsdist
... ...
cscripts/csroot2caldb.par 0 → 100644
  View file @0a8e572
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csroot2caldb.par #
  8 +# Version: 1.0.0 #
  9 +# #
  10 +# Author: Juergen Knoedlseder #
  11 +# Jurgen.Knodlseder@irap.omp.eu #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the csroot2caldb script #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# Input parameters
  19 +#=================
  20 +infile, f, a, irf.root,,, "Input ROOT file"
  21 +outdir, s, h, ,,, "Output CALDB directory"
  22 +inst, s, a, prod3,,, "Instrument name"
  23 +id, s, a, South_50h,,, "Instrument Response Function identifier"
  24 +version, s, a, prod3,,, "Instrument Response Function version"
  25 +analysis, s, a, DESY,,, "Instrument Response Function analysis"
  26 +psftype, s, h, Gauss,Gauss|King,, "Point spread function type"
  27 +rebin, b, h, no,,, "Rebin"
  28 +eascale, r, h, 1.0,,, "Scaling factor for effective areas"
  29 +bgdscale, r, h, 1.0,,, "Scaling factor for background rates"
  30 +
  31 +#
  32 +# Standard parameters
  33 +#====================
  34 +chatter, i, h, 2,0,4, "Chattiness of output"
  35 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  36 +debug, b, h, no,,, "Debugging mode activated"
  37 +mode, s, h, ql,,, "Mode of automatic parameters"
  38 +logfile, f, h, csroot2caldb.log,,, "Log filename"
... ...
cscripts/csroot2caldb.py 0 → 100755
  View file @0a8e572
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Generate IRFs in CALDB format from a ROOT offaxis performance file
  4 +#
  5 +# Copyright (C) 2016 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import os
  22 +import sys
  23 +import math
  24 +import gammalib
  25 +import ctools
  26 +from datetime import datetime
  27 +
  28 +# Optional ROOT import
  29 +try:
  30 + from ROOT import TFile, TH2F, TH3F
  31 + _has_root = True
  32 +except:
  33 + print('ROOT module not present, so csroot2caldb script will not work.')
  34 + _has_root = False
  35 +
  36 +
  37 +# ================== #
  38 +# csroot2caldb class #
  39 +# ================== #
  40 +class csroot2caldb(ctools.cscript):
  41 + """
  42 + Generates IRFs in CALDB from a ROOT offaxis performance file.
  43 + """
  44 +
  45 + # Constructor
  46 + def __init__(self, *argv):
  47 + """
  48 + CALDB constructor.
  49 +
  50 + """
  51 + # Set name
  52 + self._name = "csroot2caldb"
  53 + self._version = "1.0.0"
  54 +
  55 + # Initialise some members
  56 + self._cif = None
  57 + self._irf = None
  58 + self._ea = None
  59 + self._psf = None
  60 + self._edisp = None
  61 + self._bgd = None
  62 + self._hdu_cif = None
  63 +
  64 + # Initialise application by calling the appropriate class constructor
  65 + self._init_cscript(argv)
  66 +
  67 + # Return
  68 + return
  69 +
  70 +
  71 + # Private methods
  72 + def _get_parameters(self):
  73 + """
  74 + Get parameters from parfile.
  75 + """
  76 + # Query parameters
  77 + self["infile"].filename()
  78 + self["outdir"].string()
  79 + self["inst"].string()
  80 + self["id"].string()
  81 + self["analysis"].string()
  82 + self["version"].string()
  83 + self["psftype"].string()
  84 + self["rebin"].boolean()
  85 + self["eascale"].real()
  86 + self["bgdscale"].real()
  87 +
  88 + # Write input parameters into logger
  89 + if self._logTerse():
  90 + self._log_parameters()
  91 + self._log("\n")
  92 +
  93 + # Return
  94 + return
  95 +
  96 + def _initialise(self):
  97 +
  98 + # Set CALDB information
  99 + self._cal_tel = "CTA"
  100 + self._cal_inst = self["inst"].string().upper()
  101 + self._cal_obsid = self["id"].string()
  102 + self._cal_det = "NONE"
  103 + self._cal_flt = "NONE"
  104 + self._cal_class = "BCF"
  105 + self._cal_type = "DATA"
  106 + self._cal_qual = 0 # 0=good, 1=bad, 2=dubious, ...
  107 + self._cal_date = "14/01/30"
  108 + self._val_date = "2014-01-30"
  109 + self._val_time = "00:00:00"
  110 + self._ref_time = 51544.0
  111 + self._cal_name = "NAME("+self["id"].string()+")"
  112 + self._cal_version = "VERSION("+self["version"].string()+")"
  113 + self._cal_cut = "CLASS(BEST)"
  114 + self._cal_analysis = "ANALYSIS("+self["analysis"].string()+")"
  115 + #
  116 + self._ea_name = "EFF_AREA"
  117 + self._ea_doc = "CAL/GEN/92-019"
  118 + self._ea_bounds = [self._cal_name, self._cal_version,
  119 + self._cal_cut, self._cal_analysis]
  120 + self._ea_desc = "CTA effective area"
  121 + #
  122 + self._psf_name = "RPSF"
  123 + self._psf_doc = "CAL/GEN/92-020"
  124 + self._psf_bounds = [self._cal_name, self._cal_version,
  125 + self._cal_cut, self._cal_analysis]
  126 + self._psf_desc = "CTA point spread function"
  127 + #
  128 + self._edisp_name = "EDISP"
  129 + self._edisp_doc = "???"
  130 + self._edisp_bounds = [self._cal_name, self._cal_version,
  131 + self._cal_cut, self._cal_analysis]
  132 + self._edisp_desc = "CTA energy dispersion"
  133 + #
  134 + self._bgd_name = "BGD"
  135 + self._bgd_doc = "???"
  136 + self._bgd_bounds = [self._cal_name, self._cal_version,
  137 + self._cal_cut, self._cal_analysis]
  138 + self._bgd_desc = "CTA background"
  139 +
  140 + # Create directory structure
  141 + self._make_dirs()
  142 +
  143 + # Return
  144 + return
  145 +
  146 + def _make_dirs(self):
  147 + """
  148 + Generate CALDB directory structure for one observation identifier.
  149 + The structure is given by
  150 +
  151 + data/<tel>/<inst>/bcf/<obsid>
  152 +
  153 + where <tel> is "cta" and <inst> is the instrument specified in the
  154 + CALDB constructor (the instrument may be used for different array
  155 + configurations).
  156 + """
  157 + # Set base directory
  158 + base_dir = "data"
  159 + base_dir += "/"+self._cal_tel.lower()
  160 + base_dir += "/"+self._cal_inst.lower()
  161 +
  162 + # Set directory names
  163 + obs_dir = base_dir+"/bcf/"+self._cal_obsid
  164 + self._ea_dir = obs_dir
  165 + self._psf_dir = obs_dir
  166 + self._edisp_dir = obs_dir
  167 + self._bgd_dir = obs_dir
  168 +
  169 + # Optionally prefix path
  170 + outdir = self["outdir"].string()
  171 + if len(outdir) > 0:
  172 + self._base_path = outdir+"/"+base_dir
  173 + self._ea_path = outdir+"/"+self._ea_dir
  174 + self._psf_path = outdir+"/"+self._psf_dir
  175 + self._edisp_path = outdir+"/"+self._edisp_dir
  176 + self._bgd_path = outdir+"/"+self._bgd_dir
  177 + else:
  178 + self._base_path = base_dir
  179 + self._ea_path = self._ea_dir
  180 + self._psf_path = self._psf_dir
  181 + self._edisp_path = self._edisp_dir
  182 + self._bgd_path = self._bgd_dir
  183 +
  184 + # Create directory structure
  185 + if not os.path.isdir(self._ea_path):
  186 + os.makedirs(self._ea_path)
  187 + if not os.path.isdir(self._psf_path):
  188 + os.makedirs(self._psf_path)
  189 + if not os.path.isdir(self._edisp_path):
  190 + os.makedirs(self._edisp_path)
  191 + if not os.path.isdir(self._bgd_path):
  192 + os.makedirs(self._bgd_path)
  193 +
  194 + # Return
  195 + return
  196 +
  197 + def _open(self, version, split=False, clobber=True):
  198 + """
  199 + Open existing or create new calibration. The actual version will put
  200 + all calibrations in the same file, although each part of the response
  201 + function will have its own logical name. We can thus easily modify
  202 + the script to put each calibration information in a separate file.
  203 +
  204 + Parameters:
  205 + version - Filename version
  206 + Keywords:
  207 + split - Split IRF over several files?
  208 + clobber - Overwrite existing files?
  209 + """
  210 + # Set calibrate file names
  211 + if split:
  212 + self._ea_file = "ea_"+version+".fits"
  213 + self._psf_file = "psf_"+version+".fits"
  214 + self._edisp_file = "edisp_"+version+".fits"
  215 + self._bgd_file = "bgd_"+version+".fits"
  216 + else:
  217 + self._ea_file = "irf_"+version+".fits"
  218 + self._psf_file = "irf_"+version+".fits"
  219 + self._edisp_file = "irf_"+version+".fits"
  220 + self._bgd_file = "irf_"+version+".fits"
  221 +
  222 + # Open calibration database index
  223 + self._cif = gammalib.GFits(self._base_path+"/caldb.indx", True)
  224 +
  225 + # If file has no CIF extension than create it now
  226 + try:
  227 + self._hdu_cif = self._cif.table("CIF")
  228 + except:
  229 + self._create_cif()
  230 + self._hdu_cif = self._cif.table("CIF")
  231 +
  232 + # Set filenames
  233 + ea_filename = self._ea_path+"/"+self._ea_file
  234 + psf_filename = self._psf_path+"/"+self._psf_file
  235 + edisp_filename = self._edisp_path+"/"+self._edisp_file
  236 + bgd_filename = self._bgd_path+"/"+self._bgd_file
  237 +
  238 + # Open files
  239 + if split:
  240 + self._ea = gammalib.GFits(ea_filename, True)
  241 + self._psf = gammalib.GFits(psf_filename, True)
  242 + self._edisp = gammalib.GFits(edisp_filename, True)
  243 + self._bgd = gammalib.GFits(bgd_filename, True)
  244 + else:
  245 + self._irf = gammalib.GFits(ea_filename, True)
  246 + self._ea = self._irf
  247 + self._psf = self._irf
  248 + self._edisp = self._irf
  249 + self._bgd = self._irf
  250 +
  251 + # Open HDUs
  252 + self._hdu_ea = self._open_hdu(self._ea, "EFFECTIVE AREA",
  253 + self._ea_name, self._ea_doc)
  254 + self._hdu_psf = self._open_hdu(self._psf, "POINT SPREAD FUNCTION",
  255 + self._psf_name, self._psf_doc)
  256 + self._hdu_edisp = self._open_hdu(self._edisp, "ENERGY DISPERSION",
  257 + self._edisp_name, self._edisp_doc)
  258 + self._hdu_bgd = self._open_hdu(self._bgd, "BACKGROUND",
  259 + self._bgd_name, self._bgd_doc)
  260 +
  261 + # Return
  262 + return
  263 +
  264 + def _close(self):
  265 + """
  266 + Close calibration FITS files.
  267 + """
  268 + # Add information to CIF. We do this now as before we did not
  269 + # necessarily have all the information at hand (in particular
  270 + # about the boundaries)
  271 + self._add_cif_info()
  272 +
  273 + # Close CIF
  274 + self._close_file(self._cif)
  275 +
  276 + # If IRF file exists then close it now. All file pointers will
  277 + # be set to None
  278 + if self._irf != None:
  279 + self._close_file(self._irf)
  280 + self._ea = None
  281 + self._psf = None
  282 + self._edisp = None
  283 + self._bgd = None
  284 +
  285 + # ... otherwise we have split files, so we have to close them all
  286 + else:
  287 + self._close_file(self._ea)
  288 + self._close_file(self._psf)
  289 + self._close_file(self._edisp)
  290 + self._close_file(self._bgd)
  291 +
  292 + # Return
  293 + return
  294 +
  295 + def _close_file(self, file):
  296 + """
  297 + Close one file.
  298 + """
  299 + # If file is open then close it
  300 + if file != None:
  301 + file.save(True)
  302 + file.close()
  303 + file = None
  304 +
  305 + # Return
  306 + return
  307 +
  308 + def _open_hdu(self, fits, extname, name, doc):
  309 + """
  310 + Open HDU
  311 +
  312 + If HDU does not exist then create one.
  313 +
  314 + Parameters:
  315 + fits: FITS file
  316 + extname: Extension name
  317 + name: Name string
  318 + doc: Document string
  319 +
  320 + Returns:
  321 + FITS table
  322 + """
  323 + # Create table if it does not yet exist
  324 + if not fits.contains(extname):
  325 +
  326 + # Create binary table
  327 + table = gammalib.GFitsBinTable()
  328 +
  329 + # Set extension name
  330 + table.extname(extname)
  331 +
  332 + # Set OGIP keywords
  333 + self._set_ogip_keywords(table, doc, ["RESPONSE", name])
  334 +
  335 + # Append table to FITS file
  336 + fits.append(table)
  337 +
  338 + # Return FITS table
  339 + return fits.table(extname)
  340 +
  341 + def _set_ogip_keywords(self, hdu, hdudoc, hduclas):
  342 + """
  343 + Set standard OGIP keywords for extension.
  344 +
  345 + Parameters:
  346 + hdu: Header Data Unit.
  347 + hdudoc: Documentation reference string
  348 + hduclas: Array of HDUCLAS fields
  349 + """
  350 + # Set UTC date of file creation
  351 + utc = datetime.utcnow().isoformat()[:19]
  352 +
  353 + # Set keywords
  354 + hdu.card("ORIGIN", "IRAP", "Name of organization making this file")
  355 + hdu.card("DATE", utc, "File creation date (YYYY-MM-DDThh:mm:ss UTC)")
  356 + hdu.card("TELESCOP", self._cal_tel, "Name of telescope")
  357 + hdu.card("INSTRUME", self._cal_inst, "Name of instrument")
  358 + hdu.card("DETNAM", self._cal_det, "Name of detector")
  359 + hdu.card("HDUCLASS", "OGIP", "HDU class")
  360 + hdu.card("HDUDOC", hdudoc, "HDU documentation")
  361 + for i, item in enumerate(hduclas):
  362 + key = "HDUCLAS%d" % (i+1)
  363 + hdu.card(key, item, "HDU class")
  364 + hdu.card("HDUVERS", "1.0.0", "HDU version")
  365 +
  366 + # Return
  367 + return
  368 +
  369 + def _create_cif(self):
  370 + """
  371 + Create Calibration Index File (CIF) extension in CIF FITS file.
  372 +
  373 + Parameters:
  374 + None
  375 + Keywords:
  376 + None
  377 + """
  378 + # Create binary table
  379 + table = gammalib.GFitsBinTable()
  380 +
  381 + # Set boundary
  382 +
  383 + # Attach columns. Reference: CAL/GEN/92-008
  384 + table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
  385 + table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
  386 + table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20))
  387 + table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10))
  388 + table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20))
  389 + table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70))
  390 + table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40))
  391 + table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3))
  392 + table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4))
  393 + table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20))
  394 + table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9))
  395 + table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0))
  396 + table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10))
  397 + table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8))
  398 + table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0))
  399 + table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
  400 + table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
  401 + table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
  402 +
  403 + # Set keywords. Reference: CAL/GEN/92-008
  404 + table.extname("CIF")
  405 + table.card("CIFVERSN", "1992a", "Version of CIF format")
  406 +
  407 + # Attach table to FITS file. Note that at this moment the
  408 + # FITS table gets connected to the FITS file. Yet since
  409 + # nothing was yet written to the FITS file, we cannot read
  410 + # anything from it.
  411 + self._cif.append(table)
  412 +
  413 + # Return
  414 + return
  415 +
  416 + def _add_cif_info(self):
  417 + """
  418 + Add information to CIF extension.
  419 + """
  420 + # Append 4 rows to CIF extension
  421 + self._hdu_cif.append_rows(4)
  422 +
  423 + # Add generic information for these 4 rows
  424 + for i in range(4):
  425 +
  426 + # Set row number
  427 + row = i + self._hdu_cif.nrows() - 4
  428 +
  429 + # Set element
  430 + self._hdu_cif["TELESCOP"][row] = self._cal_tel
  431 + self._hdu_cif["INSTRUME"][row] = self._cal_inst
  432 + self._hdu_cif["DETNAM"][row] = self._cal_det
  433 + self._hdu_cif["FILTER"][row] = self._cal_flt
  434 + self._hdu_cif["CAL_DEV"][row] = "ONLINE"
  435 + self._hdu_cif["CAL_CLAS"][row] = self._cal_class
  436 + self._hdu_cif["CAL_DTYP"][row] = self._cal_type
  437 + self._hdu_cif["CAL_VSD"][row] = self._val_date
  438 + self._hdu_cif["CAL_VST"][row] = self._val_time
  439 + self._hdu_cif["REF_TIME"][row] = self._ref_time
  440 + self._hdu_cif["CAL_QUAL"][row] = self._cal_qual
  441 + self._hdu_cif["CAL_DATE"][row] = self._cal_date
  442 +
  443 + # Add effective area information
  444 + row = self._hdu_cif.nrows() - 4
  445 + self._hdu_cif["CAL_DIR"][row] = self._ea_dir
  446 + self._hdu_cif["CAL_FILE"][row] = self._ea_file
  447 + self._hdu_cif["CAL_CNAM"][row] = self._ea_name
  448 + self._hdu_cif["CAL_DESC"][row] = self._ea_desc
  449 + self._hdu_cif["CAL_XNO"][row] = 1
  450 + for n in range(9):
  451 + if n >= len(self._ea_bounds):
  452 + self._hdu_cif["CAL_CBD"][row, n] = "NONE"
  453 + else:
  454 + self._hdu_cif["CAL_CBD"][row, n] = self._ea_bounds[n]
  455 +
  456 + # Add point spread function information
  457 + row = self._hdu_cif.nrows() - 3
  458 + self._hdu_cif["CAL_DIR"][row] = self._psf_dir
  459 + self._hdu_cif["CAL_FILE"][row] = self._psf_file
  460 + self._hdu_cif["CAL_CNAM"][row] = self._psf_name
  461 + self._hdu_cif["CAL_DESC"][row] = self._psf_desc
  462 + self._hdu_cif["CAL_XNO"][row] = 1
  463 + for n in range(9):
  464 + if n >= len(self._psf_bounds):
  465 + self._hdu_cif["CAL_CBD"][row, n] = "NONE"
  466 + else:
  467 + self._hdu_cif["CAL_CBD"][row, n] = self._psf_bounds[n]
  468 +
  469 + # Add energy dispersion information
  470 + row = self._hdu_cif.nrows() - 2
  471 + self._hdu_cif["CAL_DIR"][row] = self._edisp_dir
  472 + self._hdu_cif["CAL_FILE"][row] = self._edisp_file
  473 + self._hdu_cif["CAL_CNAM"][row] = self._edisp_name
  474 + self._hdu_cif["CAL_DESC"][row] = self._edisp_desc
  475 + self._hdu_cif["CAL_XNO"][row] = 1
  476 + for n in range(9):
  477 + if n >= len(self._edisp_bounds):
  478 + self._hdu_cif["CAL_CBD"][row, n] = "NONE"
  479 + else:
  480 + self._hdu_cif["CAL_CBD"][row, n] = self._edisp_bounds[n]
  481 +
  482 + # Add background information
  483 + row = self._hdu_cif.nrows() - 1
  484 + self._hdu_cif["CAL_DIR"][row] = self._bgd_dir
  485 + self._hdu_cif["CAL_FILE"][row] = self._bgd_file
  486 + self._hdu_cif["CAL_CNAM"][row] = self._bgd_name
  487 + self._hdu_cif["CAL_DESC"][row] = self._bgd_desc
  488 + self._hdu_cif["CAL_XNO"][row] = 1
  489 + for n in range(9):
  490 + if n >= len(self._bgd_bounds):
  491 + self._hdu_cif["CAL_CBD"][row, n] = "NONE"
  492 + else:
  493 + self._hdu_cif["CAL_CBD"][row, n] = self._bgd_bounds[n]
  494 +
  495 + # Return
  496 + return
  497 +
  498 + def _set_cif_keywords(self, hdu, name, bounds, desc):
  499 + """
  500 + Set standard CIF keywords for extension.
  501 + """
  502 + # Set standard CIF keywords
  503 + hdu.card("CSYSNAME", "XMA_POL", "")
  504 + hdu.card("CCLS0001", self._cal_class, "Calibration class")
  505 + hdu.card("CDTP0001", self._cal_type, "Calibration type")
  506 + hdu.card("CCNM0001", name, "Calibration name")
  507 +
  508 + # Set boundary keywords
  509 + for n in range(9):
  510 + keyname = "CBD%d0001" % (n+1)
  511 + if n >= len(bounds):
  512 + value = "NONE"
  513 + else:
  514 + value = bounds[n]
  515 + hdu.card(keyname, value, "Calibration boundary")
  516 +
  517 + # Set validity keywords
  518 + hdu.card("CVSD0001", self._val_date, "Calibration validity start date (UTC)")
  519 + hdu.card("CVST0001", self._val_time, "Calibration validity start time (UTC)")
  520 + hdu.card("CDEC0001", desc, "Calibration description")
  521 +
  522 + # Return
  523 + return
  524 +
  525 +
  526 + # ROOT specific private members
  527 + def _root2caldb(self):
  528 + """
  529 + Translate ROOT to CALDB information.
  530 +
  531 + Parameters:
  532 + filename - ROOT 2D performance filename.
  533 + """
  534 + # Open ROOT performance file
  535 + file = TFile(self["infile"].filename().url())
  536 +
  537 + # Create effective area
  538 + self._root2ea(file, rebin=self["rebin"].boolean(),
  539 + eascale=self["eascale"].real())
  540 +
  541 + # Create point spread function
  542 + self._root2psf(file, self["psftype"].string())
  543 +
  544 + # Create energy dispersion
  545 + self._root2edisp(file)
  546 +
  547 + # Create background
  548 + #self._root2bgd(file)
  549 + self._root2bgd3D(file, bgdscale=self["bgdscale"].real())
  550 +
  551 + # Return
  552 + return
  553 +
  554 + def _root2ea(self, file, rebin=False, eascale=1.0):
  555 + """
  556 + Translate ROOT to CALDB effective area extension. The following ROOT
  557 + histograms are used:
  558 +
  559 + EffectiveAreaEtrue_offaxis -> EFFAREA
  560 + EffectiveArea_offaxis -> EFFAREA_RECO
  561 +
  562 + Parameters:
  563 + file: ROOT file.
  564 +
  565 + Keywords:
  566 + rebin: Rebin Etrue histogram (useful for Prod1 IFAE runs)
  567 + eascale: Effective area scaling factor
  568 + """
  569 + # Continue only if effective area HDU has been opened
  570 + if self._hdu_ea != None:
  571 +
  572 + # Allocate ROOT 2D arrays
  573 + etrue = TH2F()
  574 + ereco = TH2F()
  575 + file.GetObject("EffectiveAreaEtrue_offaxis", etrue)
  576 + file.GetObject("EffectiveArea_offaxis", ereco)
  577 +
  578 + # Rebin etrue histogram
  579 + if rebin:
  580 + etrue.RebinX(10)
  581 + neng = etrue.GetXaxis().GetNbins()
  582 + noffset = etrue.GetYaxis().GetNbins()
  583 + for ioff in range(noffset):
  584 + for ieng in range(neng):
  585 + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
  586 + etrue.SetBinContent(ieng+1,ioff+1,value)
  587 +
  588 + # Scale histograms (if needed)
  589 + if eascale != 1.0:
  590 + neng = etrue.GetXaxis().GetNbins()
  591 + noffset = etrue.GetYaxis().GetNbins()
  592 + for ioff in range(noffset):
  593 + for ieng in range(neng):
  594 + value = etrue.GetBinContent(ieng+1,ioff+1) * eascale
  595 + etrue.SetBinContent(ieng+1,ioff+1,value)
  596 + neng = ereco.GetXaxis().GetNbins()
  597 + noffset = ereco.GetYaxis().GetNbins()
  598 + for ioff in range(noffset):
  599 + for ieng in range(neng):
  600 + value = ereco.GetBinContent(ieng+1,ioff+1) * eascale
  601 + ereco.SetBinContent(ieng+1,ioff+1,value)
  602 +
  603 + # Write boundaries (use Ereco boundaries)
  604 + bounds = self._make_2D(ereco, self._hdu_ea, None, "m2")
  605 + for b in bounds:
  606 + self._ea_bounds.append(b)
  607 + self._set_cif_keywords(self._hdu_ea, self._ea_name, \
  608 + self._ea_bounds, self._ea_desc)
  609 +
  610 + # EFFAREA
  611 + self._make_2D(etrue, self._hdu_ea, "EFFAREA", "m2")
  612 +
  613 + # EFFAREA_RECO
  614 + self._make_2D(ereco, self._hdu_ea, "EFFAREA_RECO", "m2")
  615 +
  616 + # Return
  617 + return
  618 +
  619 + def _root2psf(self, file, psftype):
  620 + """
  621 + Translate ROOT to CALDB point spread function extension.
  622 + Parameters:
  623 + file - ROOT file.
  624 + psftype - PSF type (Gauss, King)
  625 + Keywords:
  626 + None
  627 + """
  628 + # Continue only if point spread function HDU has been opened
  629 + if self._hdu_psf != None:
  630 +
  631 + # King profile PSF
  632 + if psftype == "King":
  633 + self._root2psf_king(file)
  634 + else:
  635 + self._root2psf_gauss(file)
  636 +
  637 + # Return
  638 + return
  639 +
  640 + def _root2psf_gauss(self, file):
  641 + """
  642 + Translate ROOT to CALDB point spread function extension. The following
  643 + ROOT histograms are used:
  644 +
  645 + 1/(2*pi*SIGMA_1) -> SCALE
  646 + AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8)
  647 + 0.0 -> AMPL_2
  648 + 0.0 -> SIGMA_2
  649 + 0.0 -> AMPL_3
  650 + 0.0 -> SIGMA_3
  651 +
  652 + Parameters:
  653 + file - ROOT file.
  654 + Keywords:
  655 + None
  656 + """
  657 + # Continue only if point spread function HDU has been opened
  658 + if self._hdu_psf != None:
  659 +
  660 + # Allocate ROOT 2D array
  661 + r68 = TH2F()
  662 +
  663 + # Read 68% containment histogram
  664 + file.GetObject("AngRes_offaxis", r68)
  665 + neng = r68.GetXaxis().GetNbins()
  666 + noffset = r68.GetYaxis().GetNbins()
  667 +
  668 + # Converts 68% -> 1 sigma
  669 + r68_to_sigma = 0.6624305
  670 + for ioff in range(noffset):
  671 + for ieng in range(neng):
  672 + sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma
  673 + r68.SetBinContent(ieng+1,ioff+1,sigma)
  674 +
  675 + # Compute scale histogram
  676 + scale = r68.Clone()
  677 + for ioff in range(noffset):
  678 + for ieng in range(neng):
  679 + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1)
  680 + if integral > 0.0:
  681 + value = 1.0 / integral
  682 + else:
  683 + value = 0.0
  684 + scale.SetBinContent(ieng+1,ioff+1,value)
  685 +
  686 + # Set zero histogram
  687 + zero = r68.Clone()
  688 + for ioff in range(noffset):
  689 + for ieng in range(neng):
  690 + zero.SetBinContent(ieng+1,ioff+1,0.0)
  691 +
  692 + # Set boundaries
  693 + bounds = self._make_2D(r68, self._hdu_psf, None, "deg")
  694 + for b in bounds:
  695 + self._psf_bounds.append(b)
  696 + self._psf_bounds.append("PSF(GAUSS)")
  697 + self._set_cif_keywords(self._hdu_psf, self._psf_name, \
  698 + self._psf_bounds, self._psf_desc)
  699 +
  700 + # SCALE
  701 + #self._make_2D(scale, self._hdu_psf, "SCALE", "")
  702 +
  703 + # SIGMA_1
  704 + self._make_2D(r68, self._hdu_psf, "SIGMA_1", "deg")
  705 +
  706 + # AMPL_2
  707 + self._make_2D(zero, self._hdu_psf, "AMPL_2", "")
  708 +
  709 + # SIGMA_2
  710 + self._make_2D(zero, self._hdu_psf, "SIGMA_2", "deg")
  711 +
  712 + # AMPL_3
  713 + self._make_2D(zero, self._hdu_psf, "AMPL_3", "")
  714 +
  715 + # SIGMA_3
  716 + self._make_2D(zero, self._hdu_psf, "SIGMA_3", "deg")
  717 +
  718 + # Return
  719 + return
  720 +
  721 + def _root2psf_king(self, file):
  722 + """
  723 + Translate ROOT to CALDB point spread function extension. The following
  724 + ROOT histograms are used:
  725 +
  726 + AngRes_offaxis
  727 + AngRes80_offaxis
  728 +
  729 + Parameters:
  730 + file - ROOT file.
  731 + Keywords:
  732 + None
  733 + """
  734 + # Continue only if point spread function HDU has been opened
  735 + if self._hdu_psf != None:
  736 +
  737 + # Allocate ROOT 2D arrays
  738 + r68 = TH2F()
  739 + r80 = TH2F()
  740 +
  741 + # Read 68% and 80% containment histograms
  742 + file.GetObject("AngRes_offaxis", r68)
  743 + file.GetObject("AngRes80_offaxis", r80)
  744 + neng = r68.GetXaxis().GetNbins()
  745 + noffset = r68.GetYaxis().GetNbins()
  746 +
  747 + # Initialise parameter maps by cloning the r68 2D map
  748 + gamma2D = r68.Clone()
  749 + sigma2D = r68.Clone()
  750 +
  751 + # Compute gamma and sigma values
  752 + for ioff in range(noffset):
  753 +
  754 + # Initialise last results
  755 + last_gamma = 0.0
  756 + last_sigma = 0.0
  757 +
  758 + for ieng in range(neng):
  759 +
  760 + # Extract radii
  761 + r_68 = r68.GetBinContent(ieng+1,ioff+1)
  762 + r_80 = r80.GetBinContent(ieng+1,ioff+1)
  763 +
  764 + # Initialise results
  765 + gamma = 0.0
  766 + sigma = 0.0
  767 +
  768 + # Continue only if both radii are positive
  769 + if r_68 > 0 and r_80 > 0:
  770 +
  771 + # Derive constants for equation to solve
  772 + a = 1.0 - 0.68
  773 + b = 1.0 - 0.80
  774 + c = r_68*r_68/(r_80*r_80)
  775 +
  776 + # Solve equation (a^x-1)/(b^x-1)=c for x using secant
  777 + # method. Stop when we are better than 1e-6.
  778 + x1 = -0.5
  779 + x2 = -1.0
  780 + f1 = (math.pow(a,x1) - 1.0)/(math.pow(b,x1) - 1.0) - c
  781 + f2 = (math.pow(a,x2) - 1.0)/(math.pow(b,x2) - 1.0) - c
  782 + iter = 0
  783 + while True:
  784 + x = x1 - f1 * (x1-x2)/(f1-f2)
  785 + f = (math.pow(a,x) - 1.0)/(math.pow(b,x) - 1.0) - c
  786 + iter += 1
  787 + if abs(f) < 1.0e-6:
  788 + break
  789 + else:
  790 + f2 = f1
  791 + x2 = x1
  792 + f1 = f
  793 + x1 = x
  794 +
  795 + # Compute gamma.
  796 + if x < 0.0:
  797 + gamma = 1.0 - 1.0/x
  798 + else:
  799 + gamma = 1.0
  800 +
  801 + # Compute sigma
  802 + denominator = 2.0 * gamma * (math.pow(a, x) - 1.0)
  803 + if denominator > 0.0:
  804 + sigma = r_68 * math.sqrt(1.0/denominator)
  805 + else:
  806 + denominator = 2.0 * gamma * (math.pow(b, x) - 1.0)
  807 + if denominator > 0.0:
  808 + sigma = r_80 * math.sqrt(1.0/denominator)
  809 + else:
  810 + gamma = 0.0
  811 + sigma = 0.0
  812 +
  813 + # Handle special case that no results were found.
  814 + # This takes care of pixels that are ill defined
  815 + # in the MC file.
  816 + if gamma == 0.0 and sigma == 0.0:
  817 + gamma = last_gamma
  818 + sigma = last_sigma
  819 +
  820 + # Store surrent result as last result
  821 + last_gamma = gamma
  822 + last_sigma = sigma
  823 +
  824 + # Show results on console
  825 + #print(ioff,ieng,gamma,sigma,x,f,iter)
  826 +
  827 + # Set bin contents
  828 + gamma2D.SetBinContent(ieng+1,ioff+1,gamma)
  829 + sigma2D.SetBinContent(ieng+1,ioff+1,sigma)
  830 +
  831 + # Set boundaries
  832 + bounds = self._make_2D(r68, self._hdu_psf, None, "deg")
  833 + for b in bounds:
  834 + self._psf_bounds.append(b)
  835 + self._psf_bounds.append("PSF(KING)")
  836 + self._set_cif_keywords(self._hdu_psf, self._psf_name, \
  837 + self._psf_bounds, self._psf_desc)
  838 +
  839 + # GAMMA
  840 + self._make_2D(gamma2D, self._hdu_psf, "GAMMA", "")
  841 +
  842 + # SIGMA
  843 + self._make_2D(sigma2D, self._hdu_psf, "SIGMA", "deg")
  844 +
  845 + # Return
  846 + return
  847 +
  848 + def _root2edisp(self, file):
  849 + """
  850 + Translate ROOT to CALDB energy dispersion extension. The following ROOT
  851 + histograms are used:
  852 +
  853 + EestOverEtrue_offaxis -> MATRIX
  854 +
  855 + Parameters:
  856 + file - ROOT file.
  857 + Keywords:
  858 + None
  859 + """
  860 + # Continue only if energy dispersion HDU has been opened
  861 + if self._hdu_edisp != None:
  862 +
  863 + # Allocate ROOT 3D array
  864 + matrix = TH3F()
  865 + file.GetObject("EestOverEtrue_offaxis", matrix)
  866 +
  867 + # Set boundaries
  868 + bounds = self._make_3D_migra(matrix, self._hdu_edisp, None, "")
  869 + for b in bounds:
  870 + self._edisp_bounds.append(b)
  871 + self._set_cif_keywords(self._hdu_edisp, self._edisp_name, \
  872 + self._edisp_bounds, self._edisp_desc)
  873 +
  874 + # MATRIX
  875 + self._make_3D_migra(matrix, self._hdu_edisp, "MATRIX", "")
  876 +
  877 + # Return
  878 + return
  879 +
  880 + def _root2bgd(self, file):
  881 + """
  882 + Translate ROOT to CALDB background extension. The following ROOT
  883 + histograms are used:
  884 +
  885 + BGRatePerSqDeg_offaxis -> BGD
  886 + BGRatePerSqDeg_offaxis -> BGD_RECO
  887 +
  888 + Parameters:
  889 + file - ROOT file.
  890 + Keywords:
  891 + None
  892 + """
  893 + # Continue only if background HDU has been opened
  894 + if self._hdu_bgd != None:
  895 +
  896 + # Allocate ROOT 2D array
  897 + array = TH2F()
  898 + file.GetObject("BGRatePerSqDeg_offaxis", array)
  899 +
  900 + # Set boundaries
  901 + bounds = self._make_2D(array, self._hdu_bgd, None, "deg")
  902 + for b in bounds:
  903 + self._bgd_bounds.append(b)
  904 + self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \
  905 + self._bgd_bounds, self._bgd_desc)
  906 +
  907 + # BGD
  908 + self._make_2D(array, self._hdu_bgd, "BGD", "")
  909 +
  910 + # BGD_RECO
  911 + self._make_2D(array, self._hdu_bgd, "BGD_RECO", "")
  912 +
  913 + # Return
  914 + return
  915 +
  916 + def _root2bgd3D(self, file, bgdscale=1.0):
  917 + """
  918 + Translate ROOT to CALDB background extension. The following ROOT
  919 + histograms are used:
  920 +
  921 + BGRatePerSqDeg_offaxis -> BGD
  922 +
  923 + Parameters:
  924 + file - ROOT file.
  925 + Keywords:
  926 + scale - Background rate scaling factor
  927 + """
  928 + # Continue only if background HDU has been opened
  929 + if self._hdu_bgd != None:
  930 +
  931 + # Allocate ROOT 2D array
  932 + array = TH2F()
  933 + file.GetObject("BGRatePerSqDeg_offaxis", array)
  934 +
  935 + # Set boundaries
  936 + bounds = self._make_3D(array, self._hdu_bgd, None, "deg")
  937 + for b in bounds:
  938 + self._bgd_bounds.append(b)
  939 + self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \
  940 + self._bgd_bounds, self._bgd_desc)
  941 +
  942 + # BGD (reconstructed energy)
  943 + self._make_3D(array, self._hdu_bgd, "BGD", "1/s/MeV/sr", scale=bgdscale)
  944 +
  945 + # Return
  946 + return
  947 +
  948 + def _append_column(self, hdu, name, unit, axis, log=False):
  949 + """
  950 + Append column to HDU.
  951 +
  952 + Parameters:
  953 + hdu: HDU
  954 + name: Column name
  955 + unit: Column unit
  956 + axis: ROOT histogram axis
  957 +
  958 + Keywords:
  959 + log: Logarithmic axis
  960 + """
  961 + # Continue only if columns does not yet exist
  962 + if not hdu.contains(name):
  963 +
  964 + # Get number of axis bins
  965 + nbins = axis.GetNbins()
  966 +
  967 + # Append column and set column unit
  968 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, nbins))
  969 + hdu[name].unit(unit)
  970 +
  971 + # Do we have a LO axis? If not we have a HI axis
  972 + if name[-2:] == "LO":
  973 + low = True
  974 + else:
  975 + low = False
  976 +
  977 + # Fill column
  978 + for i in range(nbins):
  979 + if low:
  980 + value = axis.GetBinLowEdge(i+1)
  981 + else:
  982 + value = axis.GetBinUpEdge(i+1)
  983 + if log:
  984 + value = pow(10.0, value)
  985 + hdu[name][0,i] = value
  986 +
  987 + # Return
  988 + return
  989 +
  990 + def _make_2D(self, array, hdu, name, unit, scale=1.0):
  991 + """
  992 + Make 2D matrix as function of energy and offset angle from a
  993 + ROOT 2D histogram. If the HDU has already the energy and
  994 + offset angle columns, this method will simply add another data
  995 + column.
  996 + If name==None, the method will not append any data column.
  997 +
  998 + Parameters:
  999 + array - ROOT 2D histogram.
  1000 + hdu - FITS HDU.
  1001 + name - Data column name.
  1002 + unit - Data unit.
  1003 + Keywords:
  1004 + scale - Scaling factor for histogram values.
  1005 + """
  1006 + # Extract energy and offset angle vectors
  1007 + energies = array.GetXaxis()
  1008 + offsets = array.GetYaxis()
  1009 + neng = energies.GetNbins()
  1010 + noffset = offsets.GetNbins()
  1011 +
  1012 + # Append axis columns to HDU
  1013 + self._append_column(hdu, "ENERG_LO", "TeV", energies, log=True)
  1014 + self._append_column(hdu, "ENERG_HI", "TeV", energies, log=True)
  1015 + self._append_column(hdu, "THETA_LO", "deg", offsets)
  1016 + self._append_column(hdu, "THETA_HI", "deg", offsets)
  1017 +
  1018 + # Append array column
  1019 + if name != None and not hdu.contains(name):
  1020 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, neng*noffset))
  1021 + hdu[name].unit(unit)
  1022 + hdu[name].dim([neng, noffset])
  1023 + for ioff in range(noffset):
  1024 + for ieng in range(neng):
  1025 + index = ieng + ioff * neng
  1026 + value = array.GetBinContent(ieng+1,ioff+1)
  1027 + hdu[name][0,index] = value * scale
  1028 +
  1029 + # Collect boundary information
  1030 + bd_eng = "ENERG(%.4f-%.2f)TeV" % \
  1031 + (pow(10.0, energies.GetBinLowEdge(1)), \
  1032 + pow(10.0, energies.GetBinUpEdge(neng)))
  1033 + bd_off = "THETA(%.2f-%.2f)deg" % \
  1034 + (offsets.GetBinLowEdge(1), \
  1035 + offsets.GetBinUpEdge(noffset))
  1036 + bd_phi = "PHI(0-360)deg"
  1037 + bounds = [bd_eng, bd_off, bd_phi]
  1038 +
  1039 + # Return boundary information
  1040 + return bounds
  1041 +
  1042 + def _make_3D(self, array, hdu, name, unit, scale=1.0):
  1043 + """
  1044 + Make 3D cube as function of DETX, DETY and energy from a
  1045 + ROOT 2D histogram. If the HDU has already the energy and
  1046 + offset angle columns, this method will simply add another data
  1047 + column.
  1048 + If name==None, the method will not append any data column.
  1049 +
  1050 + Parameters:
  1051 + array - ROOT 2D histogram.
  1052 + hdu - FITS HDU.
  1053 + name - Data column name.
  1054 + unit - Data unit.
  1055 + Keywords:
  1056 + scale - Scaling factor for histogram values.
  1057 + """
  1058 + # Set constants
  1059 + deg2sr = 0.01745329*0.01745329
  1060 +
  1061 + # Extract energy and offset angle vectors
  1062 + energies = array.GetXaxis()
  1063 + neng = energies.GetNbins()
  1064 + offsets = array.GetYaxis()
  1065 + noffset = offsets.GetNbins()
  1066 + ewidth = [] # in MeV
  1067 + for ieng in range(neng):
  1068 + ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \
  1069 + pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0))
  1070 +
  1071 + # Build DETX and DETY axes
  1072 + ndet = array.GetYaxis().GetNbins()
  1073 + ndets = 2*ndet
  1074 + det_max = array.GetYaxis().GetBinUpEdge(ndet)
  1075 + det_bin = det_max/float(ndet)
  1076 + dets_lo = []
  1077 + dets_hi = []
  1078 + dets2 = []
  1079 + for i in range(ndets):
  1080 + det_lo = -det_max + i*det_bin
  1081 + det_hi = det_lo + det_bin
  1082 + det_val = 0.5*(det_lo + det_hi)
  1083 + dets_lo.append(det_lo)
  1084 + dets_hi.append(det_hi)
  1085 + dets2.append(det_val*det_val)
  1086 +
  1087 + # Append DETX_LO column
  1088 + if not hdu.contains("DETX_LO"):
  1089 + hdu.append(gammalib.GFitsTableFloatCol("DETX_LO", 1, ndets))
  1090 + hdu["DETX_LO"].unit("deg")
  1091 + for i in range(ndets):
  1092 + hdu["DETX_LO"][0,i] = dets_lo[i]
  1093 +
  1094 + # Append DETX_HI column
  1095 + if not hdu.contains("DETX_HI"):
  1096 + hdu.append(gammalib.GFitsTableFloatCol("DETX_HI", 1, ndets))
  1097 + hdu["DETX_HI"].unit("deg")
  1098 + for i in range(ndets):
  1099 + hdu["DETX_HI"][0,i] = dets_hi[i]
  1100 +
  1101 + # Append DETY_LO column
  1102 + if not hdu.contains("DETY_LO"):
  1103 + hdu.append(gammalib.GFitsTableFloatCol("DETY_LO", 1, ndets))
  1104 + hdu["DETY_LO"].unit("deg")
  1105 + for i in range(ndets):
  1106 + hdu["DETY_LO"][0,i] = dets_lo[i]
  1107 +
  1108 + # Append DETY_HI column
  1109 + if not hdu.contains("DETY_HI"):
  1110 + hdu.append(gammalib.GFitsTableFloatCol("DETY_HI", 1, ndets))
  1111 + hdu["DETY_HI"].unit("deg")
  1112 + for i in range(ndets):
  1113 + hdu["DETY_HI"][0,i] = dets_hi[i]
  1114 +
  1115 + # Append ENERG_LO and ENERG_HI columns
  1116 + self._append_column(hdu, "ENERG_LO", "TeV", energies, log=True)
  1117 + self._append_column(hdu, "ENERG_HI", "TeV", energies, log=True)
  1118 +
  1119 + # Append array column
  1120 + if name != None and not hdu.contains(name):
  1121 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, ndets*ndets*neng))
  1122 + hdu[name].unit(unit)
  1123 + hdu[name].dim([ndets,ndets,neng])
  1124 + for ix in range(ndets):
  1125 + for iy in range(ndets):
  1126 + for ieng in range(neng):
  1127 + index = ix + (iy + ieng * ndets) * ndets
  1128 + theta = math.sqrt(dets2[ix] + dets2[iy])
  1129 + ioff = offsets.FindBin(theta)
  1130 + if ioff > 0 and ioff <= noffset:
  1131 + binsq = array.GetBinContent(ieng+1,ioff)
  1132 + if binsq >= 0:
  1133 + value = binsq / deg2sr / (ewidth[ieng])
  1134 + hdu[name][0,index] = value * scale
  1135 +
  1136 + # Debugging
  1137 + if False:
  1138 + for ieng in range(neng):
  1139 + sum = 0.0
  1140 + for ioff in range(noffset):
  1141 + index = ieng + ioff * neng
  1142 + binsq = array.GetBinContent(ieng+1,ioff+1)
  1143 + value = binsq / deg2sr / ewidth[ieng]
  1144 + sum += value
  1145 + if ioff == 0:
  1146 + peak = value
  1147 + peaksq = binsq
  1148 + print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \
  1149 + pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq)
  1150 +
  1151 + # Collect boundary information
  1152 + bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max)
  1153 + bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max)
  1154 + bd_eng = "ENERG(%.4f-%.2f)TeV" % \
  1155 + (pow(10.0, energies.GetBinLowEdge(1)),
  1156 + pow(10.0, energies.GetBinUpEdge(neng)))
  1157 + bounds = [bd_detx, bd_dety, bd_eng]
  1158 +
  1159 + # Return boundary information
  1160 + return bounds
  1161 +
  1162 + def _make_3D_migra(self, array, hdu, name, unit, scale=1.0):
  1163 + """
  1164 + Make 3D cube as function of ETRUE, MIGRA and THETA form ROOT 3D
  1165 + histogram. If the HDU has already the columns, this method will
  1166 + simply add another data column.
  1167 + If name==None, the method will not append any data column.
  1168 +
  1169 + Parameters:
  1170 + array - ROOT 3D histogram.
  1171 + hdu - FITS HDU.
  1172 + name - Data column name.
  1173 + unit - Data unit.
  1174 + Keywords:
  1175 + scale - Scaling factor for histogram values.
  1176 + """
  1177 + # Extract Etrue, Eobs/Etrue and offset angle vectors
  1178 + etrue = array.GetXaxis()
  1179 + netrue = etrue.GetNbins()
  1180 + migra = array.GetYaxis()
  1181 + nmigra = migra.GetNbins()
  1182 + offsets = array.GetZaxis()
  1183 + noffset = offsets.GetNbins()
  1184 + ewidth = [] # in MeV
  1185 + for ieng in range(netrue):
  1186 + ewidth.append(pow(10.0, etrue.GetBinUpEdge(ieng+1) +6.0) - \
  1187 + pow(10.0, etrue.GetBinLowEdge(ieng+1)+6.0))
  1188 +
  1189 + # Append axis columns to HDU
  1190 + self._append_column(hdu, "ETRUE_LO", "TeV", etrue, log=True)
  1191 + self._append_column(hdu, "ETRUE_HI", "TeV", etrue, log=True)
  1192 + self._append_column(hdu, "MIGRA_LO", "", migra)
  1193 + self._append_column(hdu, "MIGRA_HI", "", migra)
  1194 + self._append_column(hdu, "THETA_LO", "deg", offsets)
  1195 + self._append_column(hdu, "THETA_HI", "deg", offsets)
  1196 +
  1197 + # Append migration matrix to HDU
  1198 + if name != None and not hdu.contains(name):
  1199 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, netrue*nmigra*noffset))
  1200 + hdu[name].unit(unit)
  1201 + hdu[name].dim([netrue, nmigra, noffset])
  1202 + for ioff in range(noffset):
  1203 + for imigra in range(nmigra):
  1204 + for ieng in range(netrue):
  1205 + index = ieng + (imigra + ioff * nmigra) * netrue
  1206 + value = array.GetBinContent(ieng+1,imigra+1,ioff+1)
  1207 + hdu[name][0,index] = value * scale
  1208 +
  1209 + # Collect boundary information
  1210 + bd_eng = "ETRUE(%.4f-%.2f)TeV" % \
  1211 + (pow(10.0, etrue.GetBinLowEdge(1)),
  1212 + pow(10.0, etrue.GetBinUpEdge(netrue)))
  1213 + bd_migra = "MIGRA(%.3f-%.3f)" % \
  1214 + (migra.GetBinLowEdge(1),
  1215 + migra.GetBinUpEdge(nmigra))
  1216 + bd_off = "THETA(%.2f-%.2f)deg" % \
  1217 + (offsets.GetBinLowEdge(1),
  1218 + offsets.GetBinUpEdge(noffset))
  1219 + bd_phi = "PHI(0-360)deg"
  1220 + bounds = [bd_eng, bd_migra, bd_off, bd_phi]
  1221 +
  1222 + # Return boundary information
  1223 + return bounds
  1224 +
  1225 +
  1226 + # Public methods
  1227 + def run(self):
  1228 + """
  1229 + Run the script.
  1230 + """
  1231 + # Switch screen logging on in debug mode
  1232 + if self._logDebug():
  1233 + self._log.cout(True)
  1234 +
  1235 + # Get parameters
  1236 + self._get_parameters()
  1237 +
  1238 + # Initialise IRF
  1239 + self._initialise()
  1240 +
  1241 + # Open calibration file
  1242 + self._open("file")
  1243 +
  1244 + # Translate ROOT to CALDB information
  1245 + if _has_root:
  1246 + self._root2caldb()
  1247 +
  1248 + # Close calibration file
  1249 + self._close()
  1250 +
  1251 + # Return
  1252 + return
  1253 +
  1254 + def execute(self):
  1255 + """
  1256 + Execute the script.
  1257 + """
  1258 + # Open logfile
  1259 + self.logFileOpen()
  1260 +
  1261 + # Run the script
  1262 + self.run()
  1263 +
  1264 + # Return
  1265 + return
  1266 +
  1267 +
  1268 +# ======================== #
  1269 +# Main routine entry point #
  1270 +# ======================== #
  1271 +if __name__ == '__main__':
  1272 +
  1273 + # Create instance of application
  1274 + app = csroot2caldb(sys.argv)
  1275 +
  1276 + # Execute application
  1277 + app.execute()
... ...
cscripts/cta_root2caldb.py deleted
View file @7c3dedd
1   -#! /usr/bin/env python
2   -# ==========================================================================
3   -# This script generates IRFs in the CALDB format of HEASARC using ROOT
4   -# offaxis performance files.
5   -# -------------------------------------------------------------------------
6   -# Copyright (C) 2011-2016 Juergen Knoedlseder
7   -# -------------------------------------------------------------------------
8   -#
9   -# This program is free software: you can redistribute it and/or modify
10   -# it under the terms of the GNU General Public License as published by
11   -# the Free Software Foundation, either version 3 of the License, or
12   -# (at your option) any later version.
13   -#
14   -# This program is distributed in the hope that it will be useful,
15   -# but WITHOUT ANY WARRANTY; without even the implied warranty of
16   -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17   -# GNU General Public License for more details.
18   -#
19   -# You should have received a copy of the GNU General Public License
20   -# along with this program. If not, see <http://www.gnu.org/licenses/>.
21   -#
22   -# -------------------------------------------------------------------------
23   -# Requirements:
24   -# - gammalib
25   -# - pyroot
26   -# ==========================================================================
27   -#from ROOT import TFile, TH1F, TH2F, TH3F
28   -from ROOT import TFile, TH2F, TH3F
29   -import gammalib
30   -from datetime import datetime
31   -import math
32   -#import sys
33   -#import glob
34   -import os
35   -
36   -
37   -# =========== #
38   -# CalDB class #
39   -# =========== #
40   -class caldb():
41   - """
42   - """
43   - def __init__(self, inst, obsid, path=""):
44   - """
45   - CALDB constructor. The constructor will create one CALDB entry for
46   - an observation.
47   -
48   - Parameters:
49   - inst - Instrument string.
50   - obsid - Observation ID string.
51   - Keywords:
52   - path - CALDB Root path.
53   - """
54   - # Store root path
55   - self.path = path
56   -
57   - # Initialise some members
58   - self.cif = None
59   - self.irf = None
60   - self.ea = None
61   - self.psf = None
62   - self.edisp = None
63   - self.bgd = None
64   - self.hdu_cif = None
65   -
66   - # Store UTC date of job execution
67   - self.utc = datetime.utcnow().isoformat()[:19]
68   -
69   - # Set CALDB information
70   - self.cal_tel = "CTA"
71   - self.cal_inst = inst.upper()
72   - self.cal_obsid = obsid
73   - self.cal_det = "NONE"
74   - self.cal_flt = "NONE"
75   - self.cal_class = "BCF"
76   - self.cal_type = "DATA"
77   - self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ...
78   - self.cal_date = "14/01/30"
79   - self.val_date = "2014-01-30"
80   - self.val_time = "00:00:00"
81   - self.ref_time = 51544.0
82   - self.cal_name = "NAME("+obsid+")"
83   - self.cal_version = "VERSION(Prod2)"
84   - self.cal_cut = "CLASS(BEST)"
85   - self.cal_analysis = "ANALYSIS(DESY)"
86   - #
87   - self.ea_name = "EFF_AREA"
88   - self.ea_doc = "CAL/GEN/92-019"
89   - self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
90   - self.ea_desc = "CTA effective area"
91   - #
92   - self.psf_name = "RPSF"
93   - self.psf_doc = "CAL/GEN/92-020"
94   - self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
95   - self.psf_desc = "CTA point spread function"
96   - #
97   - self.edisp_name = "EDISP"
98   - self.edisp_doc = "???"
99   - self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
100   - self.edisp_desc = "CTA energy dispersion"
101   - #
102   - self.bgd_name = "BGD"
103   - self.bgd_doc = "???"
104   - self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
105   - self.bgd_desc = "CTA background"
106   -
107   - # Create directory structure
108   - self.make_dirs()
109   -
110   - # Return
111   - return
112   -
113   - def __del__(self):
114   - """
115   - Destructor. Close all FITS files.
116   - """
117   - # Close calibration
118   - self.close()
119   -
120   - # Return
121   - return
122   -
123   - def make_dirs(self):
124   - """
125   - Generate CALDB directory structure for one observation identifier.
126   - The structure is given by
127   -
128   - data/<tel>/<inst>/bcf/<obsid>
129   -
130   - where <tel> is "cta" and <inst> is the instrument specified in the
131   - CALDB constructor (the instrument may be used for different array
132   - configurations).
133   -
134   - Parameters:
135   - None
136   - Keywords:
137   - None
138   - """
139   - # Set base path
140   - self.base_dir = "data"
141   - self.base_dir += "/"+self.cal_tel.lower()
142   - self.base_dir += "/"+self.cal_inst.lower()
143   -
144   - # Set directory names
145   - self.bcf_dir = self.base_dir+"/bcf"
146   - self.obs_dir = self.bcf_dir+"/"+self.cal_obsid
147   - self.ea_dir = self.obs_dir
148   - self.psf_dir = self.obs_dir
149   - self.edisp_dir = self.obs_dir
150   - self.bgd_dir = self.obs_dir
151   -
152   - # Optionally prefix path
153   - if len(self.path) > 0:
154   - self.base_path = self.path+"/"+self.base_dir
155   - self.bcf_path = self.path+"/"+self.bcf_dir
156   - self.obs_path = self.path+"/"+self.obs_dir
157   - self.ea_path = self.path+"/"+self.ea_dir
158   - self.psf_path = self.path+"/"+self.psf_dir
159   - self.edisp_path = self.path+"/"+self.edisp_dir
160   - self.bgd_path = self.path+"/"+self.bgd_dir
161   - else:
162   - self.base_path = self.base_dir
163   - self.bcf_path = self.bcf_dir
164   - self.obs_path = self.obs_dir
165   - self.ea_path = self.ea_dir
166   - self.psf_path = self.psf_dir
167   - self.edisp_path = self.edisp_dir
168   - self.bgd_path = self.bgd_dir
169   -
170   - # Create directory structure
171   - if not os.path.isdir(self.ea_path):
172   - os.makedirs(self.ea_path)
173   - if not os.path.isdir(self.psf_path):
174   - os.makedirs(self.psf_path)
175   - if not os.path.isdir(self.edisp_path):
176   - os.makedirs(self.edisp_path)
177   - if not os.path.isdir(self.bgd_path):
178   - os.makedirs(self.bgd_path)
179   -
180   - # Return
181   - return
182   -
183   - def open(self, version, split=False, clobber=True):
184   - """
185   - Open existing or create new calibration. The actual version will put
186   - all calibrations in the same file, although each part of the response
187   - function will have its own logical name. We can thus easily modify
188   - the script to put each calibration information in a separate file.
189   -
190   - Parameters:
191   - version - Filename version
192   - Keywords:
193   - split - Split IRF over several files?
194   - clobber - Overwrite existing files?
195   - """
196   - # Set calibrate file names
197   - if split:
198   - self.ea_file = "ea_"+version+".fits"
199   - self.psf_file = "psf_"+version+".fits"
200   - self.edisp_file = "edisp_"+version+".fits"
201   - self.bgd_file = "bgd_"+version+".fits"
202   - else:
203   - self.ea_file = "irf_"+version+".fits"
204   - self.psf_file = "irf_"+version+".fits"
205   - self.edisp_file = "irf_"+version+".fits"
206   - self.bgd_file = "irf_"+version+".fits"
207   -
208   - # Open calibration database index
209   - self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True)
210   -
211   - # If file has no CIF extension than create it now
212   - try:
213   - self.hdu_cif = self.cif.table("CIF")
214   - except:
215   - self.create_cif()
216   - self.hdu_cif = self.cif.table("CIF")
217   -
218   - # Set filenames
219   - ea_filename = self.ea_path+"/"+self.ea_file
220   - psf_filename = self.psf_path+"/"+self.psf_file
221   - edisp_filename = self.edisp_path+"/"+self.edisp_file
222   - bgd_filename = self.bgd_path+"/"+self.bgd_file
223   -
224   - # Open files
225   - if split:
226   - self.ea = gammalib.GFits(ea_filename, True)
227   - self.psf = gammalib.GFits(psf_filename, True)
228   - self.edisp = gammalib.GFits(edisp_filename, True)
229   - self.bgd = gammalib.GFits(bgd_filename, True)
230   - else:
231   - self.irf = gammalib.GFits(ea_filename, True)
232   - self.ea = self.irf
233   - self.psf = self.irf
234   - self.edisp = self.irf
235   - self.bgd = self.irf
236   -
237   - # Open HDUs
238   - self.open_ea()
239   - self.open_psf()
240   - self.open_edisp()
241   - self.open_bgd()
242   -
243   - # Return
244   - return
245   -
246   - def close(self):
247   - """
248   - Close calibration FITS files.
249   -
250   - Parameters:
251   - None
252   - Keywords:
253   - None
254   - """
255   - # Add information to CIF. We do this now as before we did not
256   - # necessarily have all the information at hand (in particular
257   - # about the boundaries)
258   - self.add_cif_info()
259   -
260   - # Close CIF
261   - if self.cif != None:
262   - self.cif.save(True)
263   - self.cif.close()
264   - self.cif = None
265   -
266   - # If IRF file exists then close it now. All file pointers will
267   - # be set to None
268   - if self.irf != None:
269   - self.irf.save(True)
270   - self.irf.close()
271   - self.irf = None
272   - self.ea = None
273   - self.psf = None
274   - self.edisp = None
275   - self.bgd = None
276   -
277   - # ... otherwise we have split files, so we have to close them
278   - # all
279   - else:
280   -
281   - # Close effective area file
282   - if self.ea != None:
283   - self.ea.save(True)
284   - self.ea.close()
285   - self.ea = None
286   -
287   - # Close point spread function file
288   - if self.psf != None:
289   - self.psf.save(True)
290   - self.psf.close()
291   - self.psf = None
292   -
293   - # Close energy dispersion file
294   - if self.edisp != None:
295   - self.edisp.save(True)
296   - self.edisp.close()
297   - self.edisp = None
298   -
299   - # Close background file
300   - if self.bgd != None:
301   - self.bgd.save(True)
302   - self.bgd.close()
303   - self.bgd = None
304   -
305   - # Return
306   - return
307   -
308   - def create_cif(self):
309   - """
310   - Create Calibration Index File (CIF) extension in CIF FITS file.
311   -
312   - Parameters:
313   - None
314   - Keywords:
315   - None
316   - """
317   - # Create binary table
318   - table = gammalib.GFitsBinTable()
319   -
320   - # Set boundary
321   -
322   - # Attach columns. Reference: CAL/GEN/92-008
323   - table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
324   - table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
325   - table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20))
326   - table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10))
327   - table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20))
328   - table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70))
329   - table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40))
330   - table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3))
331   - table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4))
332   - table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20))
333   - table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9))
334   - table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0))
335   - table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10))
336   - table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8))
337   - table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0))
338   - table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
339   - table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
340   - table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
341   -
342   - # Set keywords. Reference: CAL/GEN/92-008
343   - table.extname("CIF")
344   - table.card("CIFVERSN", "1992a", "Version of CIF format")
345   -
346   - # Attach table to FITS file. Note that at this moment the
347   - # FITS table gets connected to the FITS file. Yet since
348   - # nothing was yet written to the FITS file, we cannot read
349   - # anything from it.
350   - self.cif.append(table)
351   -
352   - # Return
353   - return
354   -
355   - def add_cif_info(self):
356   - """
357   - Add information to CIF extension.
358   -
359   - Parameters:
360   - None
361   - Keywords:
362   - None
363   - """
364   - # Append 4 rows to CIF extension
365   - self.hdu_cif.append_rows(4)
366   -
367   - # Add generic information for these 4 rows
368   - for i in range(4):
369   -
370   - # Set row number
371   - row = i+self.hdu_cif.nrows()-4
372   -
373   - # Set element
374   - self.hdu_cif["TELESCOP"][row] = self.cal_tel
375   - self.hdu_cif["INSTRUME"][row] = self.cal_inst
376   - self.hdu_cif["DETNAM"][row] = self.cal_det
377   - self.hdu_cif["FILTER"][row] = self.cal_flt
378   - self.hdu_cif["CAL_DEV"][row] = "ONLINE"
379   - self.hdu_cif["CAL_CLAS"][row] = self.cal_class
380   - self.hdu_cif["CAL_DTYP"][row] = self.cal_type
381   - self.hdu_cif["CAL_VSD"][row] = self.val_date
382   - self.hdu_cif["CAL_VST"][row] = self.val_time
383   - self.hdu_cif["REF_TIME"][row] = self.ref_time
384   - self.hdu_cif["CAL_QUAL"][row] = self.cal_qual
385   - self.hdu_cif["CAL_DATE"][row] = self.cal_date
386   -
387   - # Add effective area information
388   - row = self.hdu_cif.nrows()-4
389   - self.hdu_cif["CAL_DIR"][row] = self.ea_dir
390   - self.hdu_cif["CAL_FILE"][row] = self.ea_file
391   - self.hdu_cif["CAL_CNAM"][row] = self.ea_name
392   - self.hdu_cif["CAL_DESC"][row] = self.ea_desc
393   - self.hdu_cif["CAL_XNO"][row] = 1
394   - for n in range(9):
395   - if n >= len(self.ea_bounds):
396   - self.hdu_cif["CAL_CBD"][row, n] = "NONE"
397   - else:
398   - self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n]
399   -
400   - # Add point spread function information
401   - row = self.hdu_cif.nrows()-3
402   - self.hdu_cif["CAL_DIR"][row] = self.psf_dir
403   - self.hdu_cif["CAL_FILE"][row] = self.psf_file
404   - self.hdu_cif["CAL_CNAM"][row] = self.psf_name
405   - self.hdu_cif["CAL_DESC"][row] = self.psf_desc
406   - self.hdu_cif["CAL_XNO"][row] = 1
407   - for n in range(9):
408   - if n >= len(self.psf_bounds):
409   - self.hdu_cif["CAL_CBD"][row, n] = "NONE"
410   - else:
411   - self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n]
412   -
413   - # Add energy dispersion information
414   - row = self.hdu_cif.nrows()-2
415   - self.hdu_cif["CAL_DIR"][row] = self.edisp_dir
416   - self.hdu_cif["CAL_FILE"][row] = self.edisp_file
417   - self.hdu_cif["CAL_CNAM"][row] = self.edisp_name
418   - self.hdu_cif["CAL_DESC"][row] = self.edisp_desc
419   - self.hdu_cif["CAL_XNO"][row] = 1
420   - for n in range(9):
421   - if n >= len(self.edisp_bounds):
422   - self.hdu_cif["CAL_CBD"][row, n] = "NONE"
423   - else:
424   - self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n]
425   -
426   - # Add background information
427   - row = self.hdu_cif.nrows()-1
428   - self.hdu_cif["CAL_DIR"][row] = self.bgd_dir
429   - self.hdu_cif["CAL_FILE"][row] = self.bgd_file
430   - self.hdu_cif["CAL_CNAM"][row] = self.bgd_name
431   - self.hdu_cif["CAL_DESC"][row] = self.bgd_desc
432   - self.hdu_cif["CAL_XNO"][row] = 1
433   - for n in range(9):
434   - if n >= len(self.bgd_bounds):
435   - self.hdu_cif["CAL_CBD"][row, n] = "NONE"
436   - else:
437   - self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
438   -
439   - # Return
440   - return
441   -
442   - def open_ea(self):
443   - """
444   - Open effective area extension. If no extension was found
445   - then create one. We do not yet set any data as we don't
446   - know the dimensions of the table yet.
447   -
448   - Parameters:
449   - None
450   - Keywords:
451   - None
452   - """
453   - # Get extension. If it does not exist then create it now.
454   - try:
455   - self.hdu_ea = self.ea.table("EFFECTIVE AREA")
456   - except:
457   - # Create binary table
458   - table = gammalib.GFitsBinTable()
459   -
460   - # Set extension name
461   - table.extname("EFFECTIVE AREA")
462   -
463   - # Set OGIP keywords
464   - self.set_ogip_keywords(table, self.ea_doc, \
465   - ["RESPONSE", self.ea_name])
466   -
467   - # Append table
468   - self.ea.append(table)
469   -
470   - # Get extension
471   - self.hdu_ea = self.ea.table("EFFECTIVE AREA")
472   -
473   - # Return
474   - return
475   -
476   - def open_psf(self):
477   - """
478   - Open point spread function extension. If no extension was
479   - found then create one. We do not yet set any data as we don't
480   - know the dimensions of the table yet.
481   -
482   - Parameters:
483   - None
484   - Keywords:
485   - None
486   - """
487   - # Get extension. If it does not exist then create it now.
488   - try:
489   - self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION")
490   - except:
491   - # Create binary table
492   - table = gammalib.GFitsBinTable()
493   -
494   - # Set extension name
495   - table.extname("POINT SPREAD FUNCTION")
496   -
497   - # Set OGIP keywords
498   - self.set_ogip_keywords(table, self.psf_doc, \
499   - ["RESPONSE", self.psf_name])
500   -
501   - # Append table
502   - self.psf.append(table)
503   -
504   - # Get extension
505   - self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION")
506   -
507   - # Return
508   - return
509   -
510   - def open_edisp(self):
511   - """
512   - Open energy dispersion extension. If no extension was found
513   - then create one. We do not yet set any data as we don't
514   - know the dimensions of the table yet.
515   -
516   - Parameters:
517   - None
518   - Keywords:
519   - None
520   - """
521   - # Get extension. If it does not exist then create it now.
522   - try:
523   - self.hdu_edisp = self.edisp.table("ENERGY DISPERSION")
524   - except:
525   - # Create binary table
526   - table = gammalib.GFitsBinTable()
527   -
528   - # Set extension name
529   - table.extname("ENERGY DISPERSION")
530   -
531   - # Set OGIP keywords
532   - self.set_ogip_keywords(table, self.edisp_doc, \
533   - ["RESPONSE", self.edisp_name])
534   -
535   - # Append table
536   - self.edisp.append(table)
537   -
538   - # Get extension
539   - self.hdu_edisp = self.edisp.table("ENERGY DISPERSION")
540   -
541   - # Return
542   - return
543   -
544   - def open_bgd(self):
545   - """
546   - Open background extension. If no extension was found
547   - then create one. We do not yet set any data as we don't
548   - know the dimensions of the table yet.
549   -
550   - Parameters:
551   - None
552   - Keywords:
553   - None
554   - """
555   - # Get extension. If it does not exist then create it now.
556   - try:
557   - self.hdu_bgd = self.bgd.table("BACKGROUND")
558   - except:
559   - # Create binary table
560   - table = gammalib.GFitsBinTable()
561   -
562   - # Set extension name
563   - table.extname("BACKGROUND")
564   -
565   - # Set OGIP keywords
566   - self.set_ogip_keywords(table, self.bgd_doc, \
567   - ["RESPONSE", self.bgd_name])
568   -
569   - # Append table
570   - self.bgd.append(table)
571   -
572   - # Get extension
573   - self.hdu_bgd = self.bgd.table("BACKGROUND")
574   -
575   - # Return
576   - return
577   -
578   - def set_ogip_keywords(self, hdu, hdudoc, hduclas):
579   - """
580   - Set standard OGIP keywords for extension.
581   -
582   - Parameters:
583   - hdu - Header Data Unit.
584   - hdudoc - Documentation reference string
585   - hduclas - Array of HDUCLAS fields
586   - Keywords:
587   - None
588   - """
589   - # Set keywords
590   - hdu.card("ORIGIN", "IRAP", "Name of organization making this file")
591   - hdu.card("DATE", self.utc, "File creation date (YYYY-MM-DDThh:mm:ss UTC)")
592   - hdu.card("TELESCOP", self.cal_tel, "Name of telescope")
593   - hdu.card("INSTRUME", self.cal_inst, "Name of instrument")
594   - hdu.card("DETNAM", self.cal_det, "Name of detector")
595   - hdu.card("HDUCLASS", "OGIP", "")
596   - hdu.card("HDUDOC", hdudoc, "")
597   - for i, item in enumerate(hduclas):
598   - key = "HDUCLAS%d" % (i+1)
599   - hdu.card(key, item, "")
600   - hdu.card("HDUVERS", "1.0.0", "")
601   -
602   - # Return
603   - return
604   -
605   - def make_2D(self, array, hdu, name, unit, scale=1.0):
606   - """
607   - Make 2D matrix as function of energy and offset angle from a
608   - ROOT 2D histogram. If the HDU has already the energy and
609   - offset angle columns, this method will simply add another data
610   - column.
611   - If name==None, the method will not append any data column.
612   -
613   - Parameters:
614   - array - ROOT 2D histogram.
615   - hdu - FITS HDU.
616   - name - Data column name.
617   - unit - Data unit.
618   - Keywords:
619   - scale - Scaling factor for histogram values.
620   - """
621   - # Extract energy and offset angle vectors
622   - energies = array.GetXaxis()
623   - offsets = array.GetYaxis()
624   - neng = energies.GetNbins()
625   - noffset = offsets.GetNbins()
626   -
627   - # Attach columns to HDU. We do this only if they do not yet
628   - # exist. This allows adding columns to the same HDU in successive
629   - # calls of this method
630   - #
631   - # ENERG_LO
632   - if not hdu.contains("ENERG_LO"):
633   - hdu.append(gammalib.GFitsTableFloatCol("ENERG_LO", 1, neng))
634   - hdu["ENERG_LO"].unit("TeV")
635   - for ieng in range(neng):
636   - e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1))
637   - hdu["ENERG_LO"][0,ieng] = e_lo
638   - #
639   - # ENERG_HI
640   - if not hdu.contains("ENERG_HI"):
641   - hdu.append(gammalib.GFitsTableFloatCol("ENERG_HI", 1, neng))
642   - hdu["ENERG_HI"].unit("TeV")
643   - for ieng in range(neng):
644   - e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1))
645   - hdu["ENERG_HI"][0,ieng] = e_hi
646   - #
647   - # THETA_LO
648   - if not hdu.contains("THETA_LO"):
649   - hdu.append(gammalib.GFitsTableFloatCol("THETA_LO", 1, noffset))
650   - hdu["THETA_LO"].unit("deg")
651   - for ioff in range(noffset):
652   - o_lo = offsets.GetBinLowEdge(ioff+1)
653   - hdu["THETA_LO"][0,ioff] = o_lo
654   - #
655   - # THETA_LO
656   - if not hdu.contains("THETA_HI"):
657   - hdu.append(gammalib.GFitsTableFloatCol("THETA_HI", 1, noffset))
658   - hdu["THETA_HI"].unit("deg")
659   - for ioff in range(noffset):
660   - o_hi = offsets.GetBinUpEdge(ioff+1)
661   - hdu["THETA_HI"][0,ioff] = o_hi
662   - #
663   - # "NAME"
664   - if name != None and not hdu.contains(name):
665   - hdu.append(gammalib.GFitsTableFloatCol(name, 1, neng*noffset))
666   - hdu[name].unit(unit)
667   - hdu[name].dim([neng, noffset])
668   - for ioff in range(noffset):
669   - for ieng in range(neng):
670   - index = ieng + ioff * neng
671   - value = array.GetBinContent(ieng+1,ioff+1)
672   - hdu[name][0,index] = value * scale
673   -
674   - # Collect boundary information
675   - bd_eng = "ENERG(%.4f-%.2f)TeV" % \
676   - (pow(10.0, energies.GetBinLowEdge(1)), \
677   - pow(10.0, energies.GetBinUpEdge(neng)))
678   - bd_off = "THETA(%.2f-%.2f)deg" % \
679   - (offsets.GetBinLowEdge(1), \
680   - offsets.GetBinUpEdge(noffset))
681   - bd_phi = "PHI(0-360)deg"
682   - bounds = [bd_eng, bd_off, bd_phi]
683   -
684   - # Return boundary information
685   - return bounds
686   -
687   - def make_3D(self, array, hdu, name, unit, scale=1.0):
688   - """
689   - Make 3D cube as function of DETX, DETY and energy from a
690   - ROOT 2D histogram. If the HDU has already the energy and
691   - offset angle columns, this method will simply add another data
692   - column.
693   - If name==None, the method will not append any data column.
694   -
695   - Parameters:
696   - array - ROOT 2D histogram.
697   - hdu - FITS HDU.
698   - name - Data column name.
699   - unit - Data unit.
700   - Keywords:
701   - scale - Scaling factor for histogram values.
702   - """
703   - # Set constants
704   - deg2sr = 0.01745329*0.01745329
705   -
706   - # Extract energy and offset angle vectors
707   - energies = array.GetXaxis()
708   - neng = energies.GetNbins()
709   - offsets = array.GetYaxis()
710   - noffset = offsets.GetNbins()
711   - ewidth = [] # in MeV
712   - for ieng in range(neng):
713   - ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \
714   - pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0))
715   -
716   - # Build DETX and DETY axes
717   - ndet = array.GetYaxis().GetNbins()
718   - ndets = 2*ndet
719   - det_max = array.GetYaxis().GetBinUpEdge(ndet)
720   - det_bin = det_max/float(ndet)
721   - dets_lo = []
722   - dets_hi = []
723   - dets2 = []
724   - for i in range(ndets):
725   - det_lo = -det_max + i*det_bin
726   - det_hi = det_lo + det_bin
727   - det_val = 0.5*(det_lo + det_hi)
728   - dets_lo.append(det_lo)
729   - dets_hi.append(det_hi)
730   - dets2.append(det_val*det_val)
731   -
732   - # Attach columns to HDU. We do this only if they do not yet
733   - # exist. This allows adding columns to the same HDU in successive
734   - # calls of this method
735   - #
736   - # DETX_LO
737   - if not hdu.contains("DETX_LO"):
738   - hdu.append(gammalib.GFitsTableFloatCol("DETX_LO", 1, ndets))
739   - hdu["DETX_LO"].unit("deg")
740   - for i in range(ndets):
741   - hdu["DETX_LO"][0,i] = dets_lo[i]
742   - #
743   - # DETX_HI
744   - if not hdu.contains("DETX_HI"):
745   - hdu.append(gammalib.GFitsTableFloatCol("DETX_HI", 1, ndets))
746   - hdu["DETX_HI"].unit("deg")
747   - for i in range(ndets):
748   - hdu["DETX_HI"][0,i] = dets_hi[i]
749   - #
750   - # DETY_LO
751   - if not hdu.contains("DETY_LO"):
752   - hdu.append(gammalib.GFitsTableFloatCol("DETY_LO", 1, ndets))
753   - hdu["DETY_LO"].unit("deg")
754   - for i in range(ndets):
755   - hdu["DETY_LO"][0,i] = dets_lo[i]
756   - #
757   - # DETY_HI
758   - if not hdu.contains("DETY_HI"):
759   - hdu.append(gammalib.GFitsTableFloatCol("DETY_HI", 1, ndets))
760   - hdu["DETY_HI"].unit("deg")
761   - for i in range(ndets):
762   - hdu["DETY_HI"][0,i] = dets_hi[i]
763   - #
764   - # ENERG_LO
765   - if not hdu.contains("ENERG_LO"):
766   - hdu.append(gammalib.GFitsTableFloatCol("ENERG_LO", 1, neng))
767   - hdu["ENERG_LO"].unit("TeV")
768   - for ieng in range(neng):
769   - e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1))
770   - hdu["ENERG_LO"][0,ieng] = e_lo
771   - #
772   - # ENERG_HI
773   - if not hdu.contains("ENERG_HI"):
774   - hdu.append(gammalib.GFitsTableFloatCol("ENERG_HI", 1, neng))
775   - hdu["ENERG_HI"].unit("TeV")
776   - for ieng in range(neng):
777   - e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1))
778   - hdu["ENERG_HI"][0,ieng] = e_hi
779   - #
780   - # "NAME"
781   - if name != None and not hdu.contains(name):
782   - hdu.append(gammalib.GFitsTableFloatCol(name, 1, ndets*ndets*neng))
783   - hdu[name].unit(unit)
784   - hdu[name].dim([ndets,ndets,neng])
785   - for ix in range(ndets):
786   - for iy in range(ndets):
787   - for ieng in range(neng):
788   - index = ix + (iy + ieng * ndets) * ndets
789   - theta = math.sqrt(dets2[ix] + dets2[iy])
790   - ioff = offsets.FindBin(theta)
791   - if ioff > 0 and ioff <= noffset:
792   - binsq = array.GetBinContent(ieng+1,ioff)
793   - if binsq >= 0:
794   - value = binsq / deg2sr / (ewidth[ieng])
795   - hdu[name][0,index] = value * scale
796   - #
797   - # Debugging
798   - if False:
799   - for ieng in range(neng):
800   - sum = 0.0
801   - for ioff in range(noffset):
802   - index = ieng + ioff * neng
803   - binsq = array.GetBinContent(ieng+1,ioff+1)
804   - value = binsq / deg2sr / ewidth[ieng]
805   - sum += value
806   - if ioff == 0:
807   - peak = value
808   - peaksq = binsq
809   - print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \
810   - pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq)
811   -
812   - # Collect boundary information
813   - bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max)
814   - bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max)
815   - bd_eng = "ENERG(%.4f-%.2f)TeV" % \
816   - (pow(10.0, energies.GetBinLowEdge(1)), \
817   - pow(10.0, energies.GetBinUpEdge(neng)))
818   - bounds = [bd_detx, bd_dety, bd_eng]
819   -
820   - # Return boundary information
821   - return bounds
822   -
823   - def make_3D_migra(self, array, hdu, name, unit, scale=1.0):
824   - """
825   - Make 3D cube as function of ETRUE, MIGRA and THETA form ROOT 3D
826   - histogram. If the HDU has already the columns, this method will
827   - simply add another data column.
828   - If name==None, the method will not append any data column.
829   -
830   - Parameters:
831   - array - ROOT 3D histogram.
832   - hdu - FITS HDU.
833   - name - Data column name.
834   - unit - Data unit.
835   - Keywords:
836   - scale - Scaling factor for histogram values.
837   - """
838   - # Extract Etrue, Eobs/Etrue and offset angle vectors
839   - etrue = array.GetXaxis()
840   - netrue = etrue.GetNbins()
841   - migra = array.GetYaxis()
842   - nmigra = migra.GetNbins()
843   - offsets = array.GetZaxis()
844   - noffset = offsets.GetNbins()
845   - ewidth = [] # in MeV
846   - for ieng in range(netrue):
847   - ewidth.append(pow(10.0, etrue.GetBinUpEdge(ieng+1) +6.0) - \
848   - pow(10.0, etrue.GetBinLowEdge(ieng+1)+6.0))
849   -
850   - # Attach columns to HDU. We do this only if they do not yet
851   - # exist. This allows adding columns to the same HDU in successive
852   - # calls of this method
853   - #
854   - # ETRUE_LO
855   - if not hdu.contains("ETRUE_LO"):
856   - hdu.append(gammalib.GFitsTableFloatCol("ETRUE_LO", 1, netrue))
857   - hdu["ETRUE_LO"].unit("TeV")
858   - for ieng in range(netrue):
859   - e_lo = pow(10.0, etrue.GetBinLowEdge(ieng+1))
860   - hdu["ETRUE_LO"][0,ieng] = e_lo
861   - #
862   - # ETRUE_HI
863   - if not hdu.contains("ETRUE_HI"):
864   - hdu.append(gammalib.GFitsTableFloatCol("ETRUE_HI", 1, netrue))
865   - hdu["ETRUE_HI"].unit("TeV")
866   - for ieng in range(netrue):
867   - e_hi = pow(10.0, etrue.GetBinUpEdge(ieng+1))
868   - hdu["ETRUE_HI"][0,ieng] = e_hi
869   - #
870   - # MIGRA_LO
871   - if not hdu.contains("MIGRA_LO"):
872   - hdu.append(gammalib.GFitsTableFloatCol("MIGRA_LO", 1, nmigra))
873   - hdu["MIGRA_LO"].unit("")
874   - for imigra in range(nmigra):
875   - migra_lo = migra.GetBinLowEdge(imigra+1)
876   - hdu["MIGRA_LO"][0,imigra] = migra_lo
877   - #
878   - # MIGRA_HI
879   - if not hdu.contains("MIGRA_HI"):
880   - hdu.append(gammalib.GFitsTableFloatCol("MIGRA_HI", 1, nmigra))
881   - hdu["MIGRA_HI"].unit("")
882   - for imigra in range(nmigra):
883   - migra_hi = migra.GetBinUpEdge(imigra+1)
884   - hdu["MIGRA_HI"][0,imigra] = migra_hi
885   - #
886   - # THETA_LO
887   - if not hdu.contains("THETA_LO"):
888   - hdu.append(gammalib.GFitsTableFloatCol("THETA_LO", 1, noffset))
889   - hdu["THETA_LO"].unit("deg")
890   - for ioff in range(noffset):
891   - o_lo = offsets.GetBinLowEdge(ioff+1)
892   - hdu["THETA_LO"][0,ioff] = o_lo
893   - #
894   - # THETA_HI
895   - if not hdu.contains("THETA_HI"):
896   - hdu.append(gammalib.GFitsTableFloatCol("THETA_HI", 1, noffset))
897   - hdu["THETA_HI"].unit("deg")
898   - for ioff in range(noffset):
899   - o_hi = offsets.GetBinUpEdge(ioff+1)
900   - hdu["THETA_HI"][0,ioff] = o_hi
901   - #
902   - # "NAME"
903   - if name != None and not hdu.contains(name):
904   - hdu.append(gammalib.GFitsTableFloatCol(name, 1, netrue*nmigra*noffset))
905   - hdu[name].unit(unit)
906   - hdu[name].dim([netrue, nmigra, noffset])
907   - for ioff in range(noffset):
908   - for imigra in range(nmigra):
909   - for ieng in range(netrue):
910   - index = ieng + (imigra + ioff * nmigra) * netrue
911   - value = array.GetBinContent(ieng+1,imigra+1,ioff+1)
912   - hdu[name][0,index] = value * scale
913   -
914   - # Collect boundary information
915   - bd_eng = "ETRUE(%.4f-%.2f)TeV" % \
916   - (pow(10.0, etrue.GetBinLowEdge(1)), \
917   - pow(10.0, etrue.GetBinUpEdge(netrue)))
918   - bd_migra = "MIGRA(%.4f-%.2f)TeV" % \
919   - (pow(10.0, migra.GetBinLowEdge(1)), \
920   - pow(10.0, migra.GetBinUpEdge(nmigra)))
921   - bd_off = "THETA(%.2f-%.2f)deg" % \
922   - (offsets.GetBinLowEdge(1), \
923   - offsets.GetBinUpEdge(noffset))
924   - bd_phi = "PHI(0-360)deg"
925   - bounds = [bd_eng, bd_migra, bd_off, bd_phi]
926   -
927   - # Return boundary information
928   - return bounds
929   -
930   - def root2caldb(self, filename, rebin=False, psftype="Gauss", eascale=1.0, bgdscale=1.0):
931   - """
932   - Translate ROOT to CALDB information.
933   -
934   - Parameters:
935   - filename - ROOT 2D performance filename.
936   - Keywords:
937   - rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs)
938   - psftype - PSF type
939   - eascale - Effective area scaling factor
940   - bgdscale - Background scaling factor
941   - """
942   - # Open ROOT performance file
943   - file = TFile(filename)
944   -
945   - # Create effective area
946   - self.root2ea(file, rebin=rebin, eascale=eascale)
947   -
948   - # Create point spread function
949   - self.root2psf(file, psftype)
950   -
951   - # Create energy dispersion
952   - self.root2edisp_migra(file)
953   -
954   - # Create background
955   - #self.root2bgd(file)
956   - self.root2bgd3D(file, bgdscale=bgdscale)
957   -
958   - # Return
959   - return
960   -
961   - def root2ea(self, file, rebin=False, eascale=1.0):
962   - """
963   - Translate ROOT to CALDB effective area extension. The following ROOT
964   - histograms are used:
965   -
966   - EffectiveAreaEtrue_offaxis -> EFFAREA
967   - EffectiveArea_offaxis -> EFFAREA_RECO
968   -
969   - Parameters:
970   - file - ROOT file.
971   - Keywords:
972   - rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs)
973   - eascale - Effective area scaling factor
974   - """
975   - # Continue only if effective area HDU has been opened
976   - if self.hdu_ea != None:
977   -
978   - # Allocate ROOT 2D arrays
979   - etrue = TH2F()
980   - ereco = TH2F()
981   - file.GetObject("EffectiveAreaEtrue_offaxis", etrue)
982   - file.GetObject("EffectiveArea_offaxis", ereco)
983   -
984   - # Rebin etrue histogram
985   - if rebin:
986   - etrue.RebinX(10)
987   - neng = etrue.GetXaxis().GetNbins()
988   - noffset = etrue.GetYaxis().GetNbins()
989   - for ioff in range(noffset):
990   - for ieng in range(neng):
991   - value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
992   - etrue.SetBinContent(ieng+1,ioff+1,value)
993   -
994   - # Scale histograms (if needed)
995   - if eascale != 1.0:
996   - neng = etrue.GetXaxis().GetNbins()
997   - noffset = etrue.GetYaxis().GetNbins()
998   - for ioff in range(noffset):
999   - for ieng in range(neng):
1000   - value = etrue.GetBinContent(ieng+1,ioff+1) * eascale
1001   - etrue.SetBinContent(ieng+1,ioff+1,value)
1002   - neng = ereco.GetXaxis().GetNbins()
1003   - noffset = ereco.GetYaxis().GetNbins()
1004   - for ioff in range(noffset):
1005   - for ieng in range(neng):
1006   - value = ereco.GetBinContent(ieng+1,ioff+1) * eascale
1007   - ereco.SetBinContent(ieng+1,ioff+1,value)
1008   -
1009   - # Write boundaries (use Ereco boundaries)
1010   - bounds = self.make_2D(ereco, self.hdu_ea, None, "m2")
1011   - for b in bounds:
1012   - self.ea_bounds.append(b)
1013   - self.set_cif_keywords(self.hdu_ea, self.ea_name, \
1014   - self.ea_bounds, self.ea_desc)
1015   -
1016   - # EFFAREA
1017   - self.make_2D(etrue, self.hdu_ea, "EFFAREA", "m2")
1018   -
1019   - # EFFAREA_RECO
1020   - self.make_2D(ereco, self.hdu_ea, "EFFAREA_RECO", "m2")
1021   -
1022   - # Return
1023   - return
1024   -
1025   - def root2psf(self, file, psftype):
1026   - """
1027   - Translate ROOT to CALDB point spread function extension.
1028   - Parameters:
1029   - file - ROOT file.
1030   - psftype - PSF type (Gauss, King)
1031   - Keywords:
1032   - None
1033   - """
1034   - # Continue only if point spread function HDU has been opened
1035   - if self.hdu_psf != None:
1036   -
1037   - # King profile PSF
1038   - if psftype == "King":
1039   - self.root2psf_king(file)
1040   - else:
1041   - self.root2psf_gauss(file)
1042   -
1043   - # Return
1044   - return
1045   -
1046   - def root2psf_gauss(self, file):
1047   - """
1048   - Translate ROOT to CALDB point spread function extension. The following
1049   - ROOT histograms are used:
1050   -
1051   - 1/(2*pi*SIGMA_1) -> SCALE
1052   - AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8)
1053   - 0.0 -> AMPL_2
1054   - 0.0 -> SIGMA_2
1055   - 0.0 -> AMPL_3
1056   - 0.0 -> SIGMA_3
1057   -
1058   - Parameters:
1059   - file - ROOT file.
1060   - Keywords:
1061   - None
1062   - """
1063   - # Continue only if point spread function HDU has been opened
1064   - if self.hdu_psf != None:
1065   -
1066   - # Allocate ROOT 2D array
1067   - r68 = TH2F()
1068   -
1069   - # Read 68% containment histogram
1070   - file.GetObject("AngRes_offaxis", r68)
1071   - neng = r68.GetXaxis().GetNbins()
1072   - noffset = r68.GetYaxis().GetNbins()
1073   -
1074   - # Converts 68% -> 1 sigma
1075   - r68_to_sigma = 0.6624305
1076   - for ioff in range(noffset):
1077   - for ieng in range(neng):
1078   - sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma
1079   - r68.SetBinContent(ieng+1,ioff+1,sigma)
1080   -
1081   - # Compute scale histogram
1082   - scale = r68.Clone()
1083   - for ioff in range(noffset):
1084   - for ieng in range(neng):
1085   - integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1)
1086   - if integral > 0.0:
1087   - value = 1.0 / integral
1088   - else:
1089   - value = 0.0
1090   - scale.SetBinContent(ieng+1,ioff+1,value)
1091   -
1092   - # Set zero histogram
1093   - zero = r68.Clone()
1094   - for ioff in range(noffset):
1095   - for ieng in range(neng):
1096   - zero.SetBinContent(ieng+1,ioff+1,0.0)
1097   -
1098   - # Set boundaries
1099   - bounds = self.make_2D(r68, self.hdu_psf, None, "deg")
1100   - for b in bounds:
1101   - self.psf_bounds.append(b)
1102   - self.psf_bounds.append("PSF(GAUSS)")
1103   - self.set_cif_keywords(self.hdu_psf, self.psf_name, \
1104   - self.psf_bounds, self.psf_desc)
1105   -
1106   - # SCALE
1107   - self.make_2D(scale, self.hdu_psf, "SCALE", "")
1108   -
1109   - # SIGMA_1
1110   - self.make_2D(r68, self.hdu_psf, "SIGMA_1", "deg")
1111   -
1112   - # AMPL_2
1113   - self.make_2D(zero, self.hdu_psf, "AMPL_2", "")
1114   -
1115   - # SIGMA_2
1116   - self.make_2D(zero, self.hdu_psf, "SIGMA_2", "deg")
1117   -
1118   - # AMPL_3
1119   - self.make_2D(zero, self.hdu_psf, "AMPL_3", "")
1120   -
1121   - # SIGMA_3
1122   - self.make_2D(zero, self.hdu_psf, "SIGMA_3", "deg")
1123   -
1124   - # Return
1125   - return
1126   -
1127   - def root2psf_king(self, file):
1128   - """
1129   - Translate ROOT to CALDB point spread function extension. The following
1130   - ROOT histograms are used:
1131   -
1132   - AngRes_offaxis
1133   - AngRes80_offaxis
1134   -
1135   - Parameters:
1136   - file - ROOT file.
1137   - Keywords:
1138   - None
1139   - """
1140   - # Continue only if point spread function HDU has been opened
1141   - if self.hdu_psf != None:
1142   -
1143   - # Allocate ROOT 2D arrays
1144   - r68 = TH2F()
1145   - r80 = TH2F()
1146   -
1147   - # Read 68% and 80% containment histograms
1148   - file.GetObject("AngRes_offaxis", r68)
1149   - file.GetObject("AngRes80_offaxis", r80)
1150   - neng = r68.GetXaxis().GetNbins()
1151   - noffset = r68.GetYaxis().GetNbins()
1152   -
1153   - # Initialise parameter maps by cloning the r68 2D map
1154   - gamma2D = r68.Clone()
1155   - sigma2D = r68.Clone()
1156   -
1157   - # Compute gamma and sigma values
1158   - for ioff in range(noffset):
1159   -
1160   - # Initialise last results
1161   - last_gamma = 0.0
1162   - last_sigma = 0.0
1163   -
1164   - for ieng in range(neng):
1165   -
1166   - # Extract radii
1167   - r_68 = r68.GetBinContent(ieng+1,ioff+1)
1168   - r_80 = r80.GetBinContent(ieng+1,ioff+1)
1169   -
1170   - # Initialise results
1171   - gamma = 0.0
1172   - sigma = 0.0
1173   -
1174   - # Continue only if both radii are positive
1175   - if r_68 > 0 and r_80 > 0:
1176   -
1177   - # Derive constants for equation to solve
1178   - a = 1.0 - 0.68
1179   - b = 1.0 - 0.80
1180   - c = r_68*r_68/(r_80*r_80)
1181   -
1182   - # Solve equation (a^x-1)/(b^x-1)=c for x using secant
1183   - # method. Stop when we are better than 1e-6.
1184   - x1 = -0.5
1185   - x2 = -1.0
1186   - f1 = (math.pow(a,x1) - 1.0)/(math.pow(b,x1) - 1.0) - c
1187   - f2 = (math.pow(a,x2) - 1.0)/(math.pow(b,x2) - 1.0) - c
1188   - iter = 0
1189   - while True:
1190   - x = x1 - f1 * (x1-x2)/(f1-f2)
1191   - f = (math.pow(a,x) - 1.0)/(math.pow(b,x) - 1.0) - c
1192   - iter += 1
1193   - if abs(f) < 1.0e-6:
1194   - break
1195   - else:
1196   - f2 = f1
1197   - x2 = x1
1198   - f1 = f
1199   - x1 = x
1200   -
1201   - # Compute gamma.
1202   - if x < 0.0:
1203   - gamma = 1.0 - 1.0/x
1204   - else:
1205   - gamma = 1.0
1206   -
1207   - # Compute sigma
1208   - denominator = 2.0 * gamma * (math.pow(a, x) - 1.0)
1209   - if denominator > 0.0:
1210   - sigma = r_68 * math.sqrt(1.0/denominator)
1211   - else:
1212   - denominator = 2.0 * gamma * (math.pow(b, x) - 1.0)
1213   - if denominator > 0.0:
1214   - sigma = r_80 * math.sqrt(1.0/denominator)
1215   - else:
1216   - gamma = 0.0
1217   - sigma = 0.0
1218   -
1219   - # Handle special case that no results were found.
1220   - # This takes care of pixels that are ill defined
1221   - # in the MC file.
1222   - if gamma == 0.0 and sigma == 0.0:
1223   - gamma = last_gamma
1224   - sigma = last_sigma
1225   -
1226   - # Store surrent result as last result
1227   - last_gamma = gamma
1228   - last_sigma = sigma
1229   -
1230   - # Show results on console
1231   - #print(ioff,ieng,gamma,sigma,x,f,iter)
1232   -
1233   - # Set bin contents
1234   - gamma2D.SetBinContent(ieng+1,ioff+1,gamma)
1235   - sigma2D.SetBinContent(ieng+1,ioff+1,sigma)
1236   -
1237   - # Set boundaries
1238   - bounds = self.make_2D(r68, self.hdu_psf, None, "deg")
1239   - for b in bounds:
1240   - self.psf_bounds.append(b)
1241   - self.psf_bounds.append("PSF(KING)")
1242   - self.set_cif_keywords(self.hdu_psf, self.psf_name, \
1243   - self.psf_bounds, self.psf_desc)
1244   -
1245   - # GAMMA
1246   - self.make_2D(gamma2D, self.hdu_psf, "GAMMA", "")
1247   -
1248   - # SIGMA
1249   - self.make_2D(sigma2D, self.hdu_psf, "SIGMA", "deg")
1250   -
1251   - # Return
1252   - return
1253   -
1254   - def root2edisp_2D(self, file):
1255   - """
1256   - Translate ROOT to CALDB energy dispersion extension. The following ROOT
1257   - histograms are used:
1258   -
1259   - ERes_offaxis -> ERES
1260   - Ebias_offaxis -> EBIAS
1261   -
1262   - Parameters:
1263   - file - ROOT file.
1264   - Keywords:
1265   - None
1266   - """
1267   - # Continue only if energy dispersion HDU has been opened
1268   - if self.hdu_edisp != None:
1269   -
1270   - # Allocate ROOT 2D array
1271   - eres = TH2F()
1272   - ebias = TH2F()
1273   - file.GetObject("ERes_offaxis", eres)
1274   - file.GetObject("Ebias_offaxis", ebias)
1275   -
1276   - # Set boundaries
1277   - bounds = self.make_2D(eres, self.hdu_edisp, None, "deg")
1278   - for b in bounds:
1279   - self.edisp_bounds.append(b)
1280   - self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \
1281   - self.edisp_bounds, self.edisp_desc)
1282   -
1283   - # ERES
1284   - self.make_2D(eres, self.hdu_edisp, "ERES", "Ereco/Etrue")
1285   -
1286   - # EBIAS
1287   - self.make_2D(ebias, self.hdu_edisp, "EBIAS", "Ereco/Etrue")
1288   -
1289   - # Return
1290   - return
1291   -
1292   - def root2edisp_migra(self, file):
1293   - """
1294   - Translate ROOT to CALDB energy dispersion extension. The following ROOT
1295   - histograms are used:
1296   -
1297   - EestOverEtrue_offaxis -> MATRIX
1298   -
1299   - Parameters:
1300   - file - ROOT file.
1301   - Keywords:
1302   - None
1303   - """
1304   - # Continue only if energy dispersion HDU has been opened
1305   - if self.hdu_edisp != None:
1306   -
1307   - # Allocate ROOT 3D array
1308   - matrix = TH3F()
1309   - file.GetObject("EestOverEtrue_offaxis", matrix)
1310   -
1311   - # Set boundaries
1312   - bounds = self.make_3D_migra(matrix, self.hdu_edisp, None, "")
1313   - for b in bounds:
1314   - self.edisp_bounds.append(b)
1315   - self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \
1316   - self.edisp_bounds, self.edisp_desc)
1317   -
1318   - # MATRIX
1319   - self.make_3D_migra(matrix, self.hdu_edisp, "MATRIX", "")
1320   -
1321   - # Return
1322   - return
1323   -
1324   - def root2bgd(self, file):
1325   - """
1326   - Translate ROOT to CALDB background extension. The following ROOT
1327   - histograms are used:
1328   -
1329   - BGRatePerSqDeg_offaxis -> BGD
1330   - BGRatePerSqDeg_offaxis -> BGD_RECO
1331   -
1332   - Parameters:
1333   - file - ROOT file.
1334   - Keywords:
1335   - None
1336   - """
1337   - # Continue only if background HDU has been opened
1338   - if self.hdu_bgd != None:
1339   -
1340   - # Allocate ROOT 2D array
1341   - array = TH2F()
1342   - file.GetObject("BGRatePerSqDeg_offaxis", array)
1343   -
1344   - # Set boundaries
1345   - bounds = self.make_2D(array, self.hdu_bgd, None, "deg")
1346   - for b in bounds:
1347   - self.bgd_bounds.append(b)
1348   - self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \
1349   - self.bgd_bounds, self.bgd_desc)
1350   -
1351   - # BGD
1352   - self.make_2D(array, self.hdu_bgd, "BGD", "")
1353   -
1354   - # BGD_RECO
1355   - self.make_2D(array, self.hdu_bgd, "BGD_RECO", "")
1356   -
1357   - # Return
1358   - return
1359   -
1360   - def root2bgd3D(self, file, bgdscale=1.0):
1361   - """
1362   - Translate ROOT to CALDB background extension. The following ROOT
1363   - histograms are used:
1364   -
1365   - BGRatePerSqDeg_offaxis -> BGD
1366   -
1367   - Parameters:
1368   - file - ROOT file.
1369   - Keywords:
1370   - scale - Background rate scaling factor
1371   - """
1372   - # Continue only if background HDU has been opened
1373   - if self.hdu_bgd != None:
1374   -
1375   - # Allocate ROOT 2D array
1376   - array = TH2F()
1377   - file.GetObject("BGRatePerSqDeg_offaxis", array)
1378   -
1379   - # Set boundaries
1380   - bounds = self.make_3D(array, self.hdu_bgd, None, "deg")
1381   - for b in bounds:
1382   - self.bgd_bounds.append(b)
1383   - self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \
1384   - self.bgd_bounds, self.bgd_desc)
1385   -
1386   - # BGD (reconstructed energy)
1387   - self.make_3D(array, self.hdu_bgd, "BGD", "MeV^-1 s^-1 sr^-1", scale=bgdscale)
1388   -
1389   - # Return
1390   - return
1391   -
1392   - def set_cif_keywords(self, hdu, name, bounds, desc):
1393   - """
1394   - Set standard CIF keywords for extension.
1395   - """
1396   - # Set standard CIF keywords
1397   - hdu.card("CSYSNAME", "XMA_POL", "")
1398   - hdu.card("CCLS0001", self.cal_class, "Calibration class")
1399   - hdu.card("CDTP0001", self.cal_type, "Calibration type")
1400   - hdu.card("CCNM0001", name, "Calibration name")
1401   -
1402   - # Set boundary keywords
1403   - for n in range(9):
1404   - keyname = "CBD%d0001" % (n+1)
1405   - if n >= len(bounds):
1406   - value = "NONE"
1407   - else:
1408   - value = bounds[n]
1409   - hdu.card(keyname, value, "Calibration boundary")
1410   -
1411   - # Set validity keywords
1412   - hdu.card("CVSD0001", self.val_date, "Calibration validity start date (UTC)")
1413   - hdu.card("CVST0001", self.val_time, "Calibration validity start time (UTC)")
1414   - hdu.card("CDEC0001", desc, "Calibration description")
1415   -
1416   - # Return
1417   - return
1418   -
1419   -
1420   -# ================= #
1421   -# Set test database #
1422   -# ================= #
1423   -def set_test():
1424   - """
1425   - Set Prod1 IFAE database.
1426   - """
1427   - # Set database attributes
1428   - path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
1429   - rebin = True
1430   -
1431   - # Set database content
1432   - db = [{'inst': "e", 'id': "IFAE20120510_50h",
1433   - 'path': path, 'rebin': rebin, 'psftype': "Gauss",
1434   - 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
1435   - {'inst': "e", 'id': "IFAE20120510_50h_King",
1436   - 'path': path, 'rebin': rebin, 'psftype': "King",
1437   - 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}]
1438   -
1439   - # Return database
1440   - return db
1441   -
1442   -
1443   -# ======================= #
1444   -# Set Prod1 IFAE database #
1445   -# ======================= #
1446   -def set_prod1_ifae():
1447   - """
1448   - Set Prod1 IFAE database.
1449   - """
1450   - # Set database attributes
1451   - path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
1452   - rebin = True
1453   - psftype = "Gauss"
1454   -
1455   - # Set database content
1456   - db = [{'inst': "b", 'id': "IFAE20120510_50h",
1457   - 'path': path, 'rebin': rebin, 'psftype': psftype,
1458   - 'file': "SubarrayB_IFAE_50hours_20120510_offaxis.root"},
1459   - {'inst': "e", 'id': "IFAE20120510_50h",
1460   - 'path': path, 'rebin': rebin, 'psftype': psftype,
1461   - 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
1462   - {'inst': "i", 'id': "IFAE20120510_50h",
1463   - 'path': path, 'rebin': rebin, 'psftype': psftype,
1464   - 'file': "SubarrayI_IFAE_50hours_20120510_offaxis.root"}]
1465   -
1466   - # Return database
1467   - return db
1468   -
1469   -
1470   -# ======================= #
1471   -# Set Prod2 DESY database #
1472   -# ======================= #
1473   -def set_prod2_desy():
1474   - """
1475   - Set Prod2 DESY database.
1476   - """
1477   - # Set database attributes
1478   - path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar"
1479   - path_ten = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Tenerife_50h"
1480   - rebin = False
1481   - psftype = "Gauss"
1482   - eascale = 1.0 / 0.68 # Prod 2 DESY effective areas are for 68% containment
1483   -
1484   - # Set database content
1485   - db = [{'inst': "aar", 'id': "DESY20140105_50h",
1486   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1487   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"},
1488   - {'inst': "aar", 'id': "DESY20140105_50h_0deg",
1489   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1490   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1491   - {'inst': "aar", 'id': "DESY20140105_50h_180deg",
1492   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1493   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1494   - {'inst': "aar", 'id': "DESY20140105_50h_noLST",
1495   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1496   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a-noLST.180000s.root"},
1497   - {'inst': "aar", 'id': "DESY20140105_50h_0deg_noLST",
1498   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1499   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a-noLST.180000s.root"},
1500   - {'inst': "aar", 'id': "DESY20140105_50h_180deg_noLST",
1501   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1502   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a-noLST.180000s.root"},
1503   -
1504   - {'inst': "aar500", 'id': "DESY20140105_50h",
1505   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1506   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1507   - {'inst': "aar500", 'id': "DESY20140105_50h_0deg",
1508   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1509   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1510   - {'inst': "aar500", 'id': "DESY20140105_50h_180deg",
1511   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1512   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1513   - {'inst': "aar", 'id': "DESY20140105_50h_noLST",
1514   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1515   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a-noLST.180000s.root"},
1516   - {'inst': "aar", 'id': "DESY20140105_50h_0deg_noLST",
1517   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1518   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a-noLST.180000s.root"},
1519   - {'inst': "aar", 'id': "DESY20140105_50h_180deg_noLST",
1520   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1521   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a-noLST.180000s.root"},
1522   -
1523   - {'inst': "tenerife", 'id': "DESY20140105_50h",
1524   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1525   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1526   - {'inst': "tenerife", 'id': "DESY20140105_50h_0deg",
1527   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1528   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1529   - {'inst': "tenerife", 'id': "DESY20140105_50h_180deg",
1530   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1531   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1532   - {'inst': "tenerife", 'id': "DESY20140105_50h_noLST",
1533   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1534   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN-MST.180000s.root"},
1535   - {'inst': "tenerife", 'id': "DESY20140105_50h_0deg_noLST",
1536   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1537   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN-MST.180000s.root"},
1538   - {'inst': "tenerife", 'id': "DESY20140105_50h_180deg_noLST",
1539   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1540   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN-MST.180000s.root"}
1541   - ]
1542   -
1543   - # Return database
1544   - return db
1545   -
1546   -
1547   -# ======================== #
1548   -# Main routine entry point #
1549   -# ======================== #
1550   -if __name__ == '__main__':
1551   - """
1552   - This script translates a CTA ROOT 2D performance files into a CALDB
1553   - compliant CTA response. The script will create the response for 3
1554   - observations to illustrate the run-wise file structure.
1555   -
1556   - The CTA ROOT 2D performance files contain the following histograms:
1557   - - DiffSens_offaxis (2D) - Differential sensitivity
1558   - - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?)
1559   - - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment
1560   - - AngRes_offaxis (2D) - Angular resolution 68% containment
1561   - - AngRes80_offaxis (2D) - Angular resolution 80% containment
1562   - - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree [counts/s/deg2]
1563   - - BGRate_offaxis (2D) - Background rate
1564   - - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins)
1565   - - MigMatrix_offaxis (3D) - Migration matrix
1566   - - EestOverEtrue_offaxis (3D) - Migration matrix
1567   - - ERes_offaxis (2D) - Energy Resolution
1568   - - Ebias_offaxis (2D) - Energy bias
1569   - """
1570   - # Get database
1571   - #entries = set_test()
1572   - #entries = set_prod1_ifae()
1573   - entries = set_prod2_desy()
1574   -
1575   - # Loop over entries
1576   - for entry in entries:
1577   -
1578   - # Get attributes
1579   - inst = entry['inst']
1580   - path = entry['path']
1581   - id = entry['id']
1582   - rebin = entry['rebin']
1583   - psftype = entry['psftype']
1584   - if entry.has_key('eascale'):
1585   - eascale = entry['eascale']
1586   - else:
1587   - eascale = 1.0
1588   - if entry.has_key('bgdscale'):
1589   - bgdscale = entry['bgdscale']
1590   - else:
1591   - bgdscale = 1.0
1592   - filename = path+"/"+entry['file']
1593   -
1594   - # Allocate caldb entry
1595   - irf = caldb(inst, id)
1596   -
1597   - # Open calibration file
1598   - irf.open("file")
1599   -
1600   - # Translate ROOT to CALDB information
1601   - irf.root2caldb(filename, rebin=rebin, psftype=psftype, eascale=eascale, bgdscale=bgdscale)
examples/generate_prod3_irfs.py 0 → 100755
  View file @0a8e572
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Generates Prod3 IRFs from ROOT performance files.
  4 +#
  5 +# Copyright (C) 2016 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import cscripts
  23 +
  24 +
  25 +# ======================= #
  26 +# Set Prod3 DESY database #
  27 +# ======================= #
  28 +def set_prod3_desy():
  29 + """
  30 + Set Prod3 DESY database.
  31 + """
  32 + # Set database attributes
  33 + path = "/project-data/cta/performance/prod3"
  34 + psftype = "Gauss"
  35 +
  36 + # Set database content
  37 + db = [{'inst': "prod3", 'id': "South_50h", 'psftype': psftype, 'path': path,
  38 + 'file': "DESY.d20160304.V5.ID0NIM2.prod3-paranalp05-NN.S.3HB1-ND-3.180000s.root"}
  39 + ]
  40 +
  41 + # Return database
  42 + return db
  43 +
  44 +
  45 +# ======================== #
  46 +# Main routine entry point #
  47 +# ======================== #
  48 +if __name__ == '__main__':
  49 +
  50 + # Get database
  51 + entries = set_prod3_desy()
  52 +
  53 + # Loop over entries
  54 + for entry in entries:
  55 +
  56 + # Set optional attributes
  57 + if entry.has_key('rebin'):
  58 + rebin = entry['rebin']
  59 + else:
  60 + rebin = False
  61 + if entry.has_key('eascale'):
  62 + eascale = entry['eascale']
  63 + else:
  64 + eascale = 1.0
  65 + if entry.has_key('bgdscale'):
  66 + bgdscale = entry['bgdscale']
  67 + else:
  68 + bgdscale = 1.0
  69 +
  70 + # Set filename
  71 + filename = entry['path']+"/"+entry['file']
  72 +
  73 + # Set-up csroot2caldb
  74 + caldb = cscripts.csroot2caldb()
  75 + caldb["infile"] = filename
  76 + caldb["outdir"] = "prod3"
  77 + caldb["inst"] = entry['inst']
  78 + caldb["id"] = entry['id']
  79 + caldb["version"] = entry['inst']
  80 + caldb["analysis"] = "DESY"
  81 + caldb["psftype"] = entry['psftype']
  82 + caldb["rebin"] = rebin
  83 + caldb["eascale"] = eascale
  84 + caldb["bgdscale"] = bgdscale
  85 + caldb["logfile"] = "generate_prod3_irfs.log"
  86 +
  87 + # Add CALDB entry
  88 + caldb.execute()
... ...
test/data/irf.root 0 → 100644
View file @0a8e572
No preview for this file type
test/test_csroot2caldb.py 0 → 100644
  View file @0a8e572
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the csroot2caldb script.
  4 +#
  5 +# Copyright (C) 2016 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import os
  22 +import gammalib
  23 +import cscripts
  24 +
  25 +
  26 +# ================================== #
  27 +# Test class for csroot2caldb script #
  28 +# ================================== #
  29 +class Test(gammalib.GPythonTestSuite):
  30 + """
  31 + Test class for csroot2caldb script.
  32 +
  33 + This test class makes unit tests for the cslightcrv script by using it
  34 + from the command line and from Python.
  35 + """
  36 +
  37 + # Constructor
  38 + def __init__(self):
  39 + """
  40 + Constructor.
  41 + """
  42 + # Call base class constructor
  43 + gammalib.GPythonTestSuite.__init__(self)
  44 +
  45 + # Clean calibration database
  46 + os.system("rm -rf csroot2caldb_* >/dev/null 2>&1")
  47 +
  48 + # Return
  49 + return
  50 +
  51 + # Set test functions
  52 + def set(self):
  53 + """
  54 + Set all test functions.
  55 + """
  56 + # Set test name
  57 + self.name("csroot2caldb")
  58 +
  59 + # Append tests
  60 + self.append(self._test_cmd, "Test csroot2caldb on command line")
  61 + self.append(self._test_python, "Test csroot2caldb from Python")
  62 +
  63 + # Return
  64 + return
  65 +
  66 + # Test csroot2caldb on command line
  67 + def _test_cmd(self):
  68 + """
  69 + Test csroot2caldb on the command line.
  70 + """
  71 + # Kluge to set the command (installed version has no README file)
  72 + if os.path.isfile("README"):
  73 + csroot2caldb = "../cscripts/csroot2caldb.py"
  74 + else:
  75 + csroot2caldb = "csroot2caldb"
  76 +
  77 + # Setup csroot2caldb command
  78 + cmd = csroot2caldb+' infile="data/irf.root"'+ \
  79 + ' outdir="csroot2caldb_cmd1"'+ \
  80 + ' inst=prod3 id=South_50h version=prod3'+ \
  81 + ' analysis=DESY'+ \
  82 + ' logfile="csroot2caldb_cmd1.log" chatter=1'
  83 +
  84 + # Execute csroot2caldb, make sure we catch any exception
  85 + try:
  86 + rc = os.system(cmd+" >/dev/null 2>&1")
  87 + except:
  88 + pass
  89 +
  90 + # Check if execution was successful
  91 + self.test_assert(rc == 0,
  92 + "Successful csroot2caldb execution on command line")
  93 +
  94 + # Return
  95 + return
  96 +
  97 + # Test csroot2caldb from Python
  98 + def _test_python(self):
  99 + """
  100 + Test csroot2caldb from Python.
  101 + """
  102 + # Set-up csroot2caldb
  103 + caldb = cscripts.csroot2caldb()
  104 + caldb["infile"] = "data/irf.root"
  105 + caldb["outdir"] = "csroot2caldb_py1"
  106 + caldb["inst"] = "prod3"
  107 + caldb["id"] = "South_50h"
  108 + caldb["version"] = "prod3"
  109 + caldb["analysis"] = "DESY"
  110 + caldb["logfile"] = "csroot2caldb_py1.log"
  111 + caldb["chatter"] = 2
  112 +
  113 + # Run script
  114 + caldb.logFileOpen() # Make sure we get a log file
  115 + caldb.run()
  116 +
  117 + # Return
  118 + return
... ...
test/test_python_cscripts.py
  View file @0a8e572
... ... @@ -31,6 +31,7 @@ import test_csobs2caldb
31 31 import test_csobsdef
32 32 import test_csobsinfo
33 33 import test_cspull
  34 +import test_csroot2caldb
34 35 import test_csfindobs
35 36 import test_csiactcopy
36 37 import test_csiactdata
... ... @@ -92,6 +93,7 @@ def test(installed=False, debug=False):
92 93 suite_csobsdef = test_csobsdef.Test()
93 94 suite_csobsinfo = test_csobsinfo.Test()
94 95 suite_cspull = test_cspull.Test()
  96 + suite_csroot2caldb = test_csroot2caldb.Test()
95 97  
96 98 # Setup unit tests
97 99 suite_cscaldb.set()
... ... @@ -102,6 +104,7 @@ def test(installed=False, debug=False):
102 104 suite_csobsdef.set()
103 105 suite_csobsinfo.set()
104 106 suite_cspull.set()
  107 + suite_csroot2caldb.set()
105 108  
106 109 # Append tests to container
107 110 suites.append(suite_cscaldb)
... ... @@ -112,6 +115,7 @@ def test(installed=False, debug=False):
112 115 suites.append(suite_csobsdef)
113 116 suites.append(suite_csobsinfo)
114 117 suites.append(suite_cspull)
  118 + suites.append(suite_csroot2caldb)
115 119  
116 120 # Append tests for Python 2.6+ (the IACT cscripts depend on the json
117 121 # module which is only available since Python 2.6+
... ...