Karim Sabri / ctools

Commit 838aa04770568e80f6041512867d8606fd31497f

Authored by Jürgen Knödlseder
1 parent bab52d03

Add comgendb script (#3983)

Inline Side-by-side
Showing 14 changed files with 2307 additions and 10 deletions
.gitignore
  View file @838aa04
... ... @@ -91,6 +91,7 @@ test/comobsback_*
91 91 test/comobsadd_*
92 92 test/comobssim_*
93 93 test/comlixfit_*
  94 +test/comgendb_*
94 95 test/*.fits
95 96 test/*.dat
96 97 test/*.xml
... ...
ChangeLog
  View file @838aa04
1   -2022-02-22
  1 +2022-02-25
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
  6 + Add comgendb script (#3983)
6 7 Rename virtual abstract ctool::run() method to ctool::process() and
7 8 add non-abstract ctool::run() method (#3592) (#3592)
8 9 Add methods to access ctlike attributes
... ...
1 1 New Features and Important Changes in ctools 2.0.0
2 2  
3   -22 February 2022
  3 +25 February 2022
4 4  
5 5  
6 6 Introduction
... ... @@ -402,6 +402,11 @@ added to specify the width, height and position angle of the rectangular sky
402 402 regions (#3248).
403 403  
404 404  
  405 +comgendb - Generate COMPTEL database
  406 +------------------------------------
  407 +Added script (#3983).
  408 +
  409 +
405 410 comlixfit - COMPTEL maximum likelihood fitting using SRCLIX algorithm
406 411 ---------------------------------------------------------------------
407 412 Added script (#3738).
... ...
README.md
  View file @838aa04
1 1 ctools information
2 2 ==================
3   -* Version: 2.0.0.dev (22 February 2022)
  3 +* Version: 2.0.0.dev (25 February 2022)
4 4 * GammaLib dependency: 2.0.0.dev
5 5  
6 6 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/)
... ...
doc/source/admin/release_history/2.0.rst
  View file @838aa04
1 1 .. _2.0:
2 2  
3   -ctools 2.0 release (22 February 2022)
  3 +ctools 2.0 release (25 February 2022)
4 4 =====================================
5 5  
6 6 Overview
... ... @@ -59,6 +59,8 @@ Improvements
59 59 New features
60 60 ------------
61 61  
  62 +* [`3983 <https://cta-redmine.irap.omp.eu/issues/3983>`_] -
  63 + Add :ref:`comgendb` script
62 64 * [`3738 <https://cta-redmine.irap.omp.eu/issues/3738>`_] -
63 65 Add COMPTEL analysis scripts
64 66 * Add :ref:`csadd2caldb` script
... ...
doc/source/users/reference_manual/comgendb.rst 0 → 100755
  View file @838aa04
  1 +.. _comgendb:
  2 +
  3 +comgendb
  4 +========
  5 +
  6 +Generates COMPTEL database.
  7 +
  8 +Synopsis
  9 +--------
  10 +
  11 +This script generates a COMPTEL database from the COMPTEL HEASARC archive. The
  12 +COMPTEL database will then subsequently be used by ``comobsselect``.
  13 +
  14 +On input the script expects a directory as downloaded from the HEASARC archive
  15 +containing subdirectories named ``phase01``, ``phase02`` etc. containing the
  16 +individual viewing periods for each CGRO mission phase.
  17 +
  18 +On output the script will create a COMPTEL database directory that will contains
  19 +the database FITS file ``dbase.fits`` as well as subdirectories names ``tim``
  20 +and ``xml`` containing Good Time Interval FITS files and observation definition
  21 +XML files for each COMPTEL viewing period, respectively.
  22 +
  23 +Optionally a reference database in form of an ASCII file can be specified to
  24 +extract status information for each viewing period.
  25 +
  26 +Furthermore, if ``download=yes`` is specified the script will automatically
  27 +download the COMPTEL HEASARC archive and store the data into the ``archive``
  28 +directory. Data files that do already exists in the ``archive`` directory are
  29 +not downloaded again.
  30 +
  31 +
  32 +General parameters
  33 +------------------
  34 +
  35 +``archive [string]``
  36 + COMPTEL HEASARC archive directory. This directory contains the archive as
  37 + downloaded from https://heasarc.gsfc.nasa.gov/FTP/compton/data/comptel/.
  38 +
  39 +``dbase [string]``
  40 + COMPTEL database directory. This directory will contain the COMPTEL database
  41 + on output.
  42 +
  43 +``(refdata = NONE) [file]``
  44 + COMPTEL reference database to extract status information.
  45 +
  46 +``(quality = 200) [integer]``
  47 + Target quality flag of data.
  48 +
  49 +``(download = no) [boolean]``
  50 + Download COMPTEL HEASARC archive.
  51 +
  52 +
  53 +Standard parameters
  54 +-------------------
  55 +
  56 +``(chatter = 2) [integer]``
  57 + Verbosity of the executable:
  58 + ``chatter = 0``: no information will be logged
  59 +
  60 + ``chatter = 1``: only errors will be logged
  61 +
  62 + ``chatter = 2``: errors and actions will be logged
  63 +
  64 + ``chatter = 3``: report about the task execution
  65 +
  66 + ``chatter = 4``: detailed report about the task execution
  67 +
  68 +``(clobber = yes) [boolean]``
  69 + Specifies whether an existing energy boundaries output file should be overwritten.
  70 +
  71 +``(debug = no) [boolean]``
  72 + Enables debug mode. In debug mode the executable will dump any log file output to the console.
  73 +
  74 +``(mode = ql) [string]``
  75 + Mode of automatic parameters (default is ``ql``, i.e. "query and learn").
  76 +
  77 +``(logfile = comgendb.log) [filename]``
  78 + Log filename.
  79 +
  80 +
  81 +Related tools or scripts
  82 +------------------------
  83 +
  84 +:ref:`comobsselect`
... ...
doc/source/users/reference_manual/index.rst
  View file @838aa04
... ... @@ -98,6 +98,7 @@ COMPTEL Science Analysis scripts
98 98 .. toctree::
99 99 :maxdepth: 1
100 100  
  101 + comgendb --- Generate COMPTEL database <comgendb>
101 102 comlixfit --- Fit model to data using SRCLIX algorithm <comlixfit>
102 103 comlixmap --- Create TS map using SRCLIX algorithm <comlixmap>
103 104 comobsadd --- Combine observations <comobsadd>
... ...
modules/comscripts/Makefile.am
  View file @838aa04
... ... @@ -21,6 +21,7 @@
21 21  
22 22 # Scripts to be distributed in $(pyexecdir)/comscripts
23 23 pythonscripts = $(top_builddir)/modules/comscripts/__init__.py \
  24 + $(srcdir)/comgendb.py \
24 25 $(srcdir)/comlixfit.py \
25 26 $(srcdir)/comlixmap.py \
26 27 $(srcdir)/comobsadd.py \
... ... @@ -34,7 +35,8 @@ pythonscripts = $(top_builddir)/modules/comscripts/__init__.py \
34 35 $(srcdir)/comsrcdetect.py
35 36  
36 37 # Parfiles to be distributed in $(prefix)/syspfiles
37   -parfiles = $(srcdir)/comlixfit.par \
  38 +parfiles = $(srcdir)/comgendb.par \
  39 + $(srcdir)/comlixfit.par \
38 40 $(srcdir)/comlixmap.par \
39 41 $(srcdir)/comobsadd.par \
40 42 $(srcdir)/comobsback.par \
... ... @@ -47,7 +49,8 @@ parfiles = $(srcdir)/comlixfit.par \
47 49 $(srcdir)/comsrcdetect.par
48 50  
49 51 # Scripts for which symbolic links will be created in bindir
50   -linkscripts = comlixfit \
  52 +linkscripts = comgendb \
  53 + comlixfit \
51 54 comlixmap \
52 55 comobsadd \
53 56 comobsback \
... ...
modules/comscripts/__init__.py.in
  View file @838aa04
... ... @@ -17,7 +17,8 @@
17 17 # along with this program. If not, see <http://www.gnu.org/licenses/>.
18 18 #
19 19 # ==========================================================================
20   -__all__ = ['comlixfit',
  20 +__all__ = ['comgendb',
  21 + 'comlixfit',
21 22 'comlixmap',
22 23 'comobsadd',
23 24 'comobsselect',
... ... @@ -30,6 +31,7 @@ __all__ = [&#39;comlixfit&#39;,
30 31 'comsrcdetect']
31 32  
32 33 # Import
  34 +from comscripts.comgendb import comgendb
33 35 from comscripts.comlixfit import comlixfit
34 36 from comscripts.comlixmap import comlixmap
35 37 from comscripts.comobsadd import comobsadd
... ...
modules/comscripts/comgendb.par 0 → 100644
  View file @838aa04
  1 +#############################################################################
  2 +# #
  3 +# COMPTEL Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: comgendb.par #
  8 +# #
  9 +# Author: Juergen Knoedlseder #
  10 +# Jurgen.Knodlseder@irap.omp.eu #
  11 +# IRAP #
  12 +# #
  13 +# Purpose: Parameter file for the comgendb script #
  14 +# #
  15 +#############################################################################
  16 +#
  17 +# Input and output parameters
  18 +#============================
  19 +archive, s, a, $COMDATA,,, "COMPTEL HEASARC archive directory"
  20 +dbase, s, a, $COMDATA/dbase,,, "COMPTEL database directory"
  21 +refdata, f, h, NONE,,, "COMPTEL reference database to extract status information"
  22 +
  23 +#
  24 +# Script parameters
  25 +#==================
  26 +quality, i, h, 200,,, "Target quality flag of data"
  27 +download, b, h, no,,, "Download COMPTEL HEASARC archive?"
  28 +
  29 +#
  30 +# Standard parameters
  31 +#====================
  32 +chatter, i, h, 2,0,4, "Chattiness of output"
  33 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  34 +debug, b, h, no,,, "Activate debugging mode?"
  35 +mode, s, h, ql,,, "Mode of automatic parameters"
  36 +logfile, f, h, comgendb.log,,, "Log filename"
... ...
modules/comscripts/comgendb.py 0 → 100755
  View file @838aa04
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Generate COMPTEL database from HEASARC archive
  4 +#
  5 +# Copyright (C) 2022 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 sys
  22 +import glob
  23 +import os
  24 +from distutils.spawn import find_executable
  25 +import gammalib
  26 +import ctools
  27 +from cscripts import modutils
  28 +
  29 +
  30 +# ============== #
  31 +# comgendb class #
  32 +# ============== #
  33 +class comgendb(ctools.cscript):
  34 + """
  35 + Generate COMPTEL database from HEASARC archive
  36 + """
  37 + # Constructor
  38 + def __init__(self, *argv):
  39 + """
  40 + Constructor
  41 + """
  42 + # Initialise application by calling the base class constructor
  43 + self._init_cscript(self.__class__.__name__, ctools.__version__, argv)
  44 +
  45 + # Initialise members
  46 + self._evp_list = []
  47 + self._oad_list = []
  48 + self._tim_list = []
  49 + self._xml_list = []
  50 +
  51 + # Return
  52 + return
  53 +
  54 + # Private methods
  55 + def _get_parameters(self):
  56 + """
  57 + Get parameters from parfile
  58 + """
  59 + # Query parameters
  60 + self['archive'].string()
  61 + self['refdata'].query()
  62 + self['quality'].integer()
  63 + self['download'].boolean()
  64 +
  65 + # Query ahead output model filename
  66 + if self._read_ahead():
  67 + self['dbase'].string()
  68 +
  69 + # Write input parameters into logger
  70 + self._log_parameters(gammalib.TERSE)
  71 +
  72 + # Return
  73 + return
  74 +
  75 + def _download(self):
  76 + """
  77 + Download COMPTEL HEASARC archive
  78 + """
  79 + # Log header
  80 + self._log_header1(gammalib.TERSE, 'Download HEASARC archive')
  81 +
  82 + # Check whether wget executable exists
  83 + if find_executable('wget') is None:
  84 + msg = ('"wget" command not found on system, unable to download '
  85 + 'COMPTEL HEASARC archive.')
  86 + raise RuntimeError(msg)
  87 +
  88 + # Set archive directory name
  89 + dirname = gammalib.expand_env(self['archive'].string())
  90 +
  91 + # Build wget command string
  92 + cmd = ('wget -r -N -nH --cut-dirs=4 -A "*.fits,*.fits.gz" -e robots=off '
  93 + '-P %s -q https://heasarc.gsfc.nasa.gov/FTP/compton/data/comptel/' %
  94 + (dirname))
  95 +
  96 + # Launch download
  97 + os.system(cmd)
  98 +
  99 + # Return
  100 + return
  101 +
  102 + def _get_phases(self):
  103 + """
  104 + Get available phases
  105 + """
  106 + # Set archive directory name
  107 + dirname = gammalib.expand_env(self['archive'].string())
  108 +
  109 + # Raise an exception if the archive does not exist
  110 + if not os.path.isdir(dirname):
  111 + msg = ('Archive directory "%s" does not exist.' % dirname)
  112 + raise RuntimeError(msg)
  113 +
  114 + # Set search pattern
  115 + pattern = '%s/phase*' % (gammalib.expand_env(self['archive'].string()))
  116 +
  117 + # Get available phases
  118 + phases = glob.glob(pattern)
  119 + phases.sort()
  120 +
  121 + # Return phases
  122 + return phases
  123 +
  124 + def _get_evp_list(self):
  125 + """
  126 + Get EVP list
  127 +
  128 + Returns
  129 + -------
  130 + evp_list : list of dict
  131 + EVP database
  132 + """
  133 + # Initialse EVP list
  134 + evp_list = []
  135 +
  136 + # Get available phases
  137 + phases = self._get_phases()
  138 +
  139 + # Loop over phases
  140 + for phase in phases:
  141 +
  142 + # Get available VPs
  143 + vps = glob.glob('%s/vp*' % phase)
  144 + vps.sort()
  145 +
  146 + # Loop over VPs
  147 + for vp in vps:
  148 +
  149 + # Get EVPs for VP
  150 + evps = self._get_evps_for_vp(vp)
  151 +
  152 + # Append EVPs to list
  153 + evp_list.extend(evps)
  154 +
  155 + # Sort EVP list by start day
  156 + evp_list = sorted(evp_list, key=lambda k: k['tjd_start'])
  157 +
  158 + # Return EVP list
  159 + return evp_list
  160 +
  161 + def _get_evps_for_vp(self, vp):
  162 + """
  163 + Build EVP database entries for Viewing Period
  164 +
  165 + Parameters
  166 + ----------
  167 + vp : str
  168 + VP path
  169 +
  170 + Returns
  171 + -------
  172 + evp_list : list of dict
  173 + EVP database entries
  174 + """
  175 + # Get archive name and expanded path
  176 + archive = self['archive'].string()
  177 + archive_path = gammalib.expand_env(archive)
  178 +
  179 + # Get VP name
  180 + vpname = os.path.basename(vp)
  181 +
  182 + # Initialse EVP list
  183 + evp_list = []
  184 +
  185 + # Get EVP files
  186 + evps = glob.glob('%s/*evp.fits*' % vp)
  187 +
  188 + # Loop over all EVP files
  189 + for evp in evps:
  190 +
  191 + # Set EVP filename
  192 + evpname = archive+'/'+os.path.relpath(evp, archive_path).strip('.gz')
  193 +
  194 + # Try opening file
  195 + try:
  196 + fits = gammalib.GFits(evpname)
  197 + except:
  198 + self._log_string(gammalib.NORMAL, '*** Unable to open EVP file'
  199 + ' "%s". Skip file.' % (evpname))
  200 + continue
  201 +
  202 + # Get information from EVP file
  203 + evp_dict = self._get_evp_info(evpname)
  204 +
  205 + # Skip file if it has a negative quality
  206 + if evp_dict['quality'] < 0:
  207 + self._log_string(gammalib.NORMAL, '*** EVP file "%s" has negative '
  208 + 'quality %d. Skip file.' % (evpname, evp_dict['quality']))
  209 + continue
  210 +
  211 + # Skip file if the duration derived from the events is not consistent
  212 + # with the duration derived from the header
  213 + duration_evt = evp_dict['tjd_stop'] - evp_dict['tjd_start']
  214 + duration_hdr = evp_dict['vieday'] - evp_dict['visday']
  215 + difference = abs(duration_evt-duration_hdr)
  216 + if difference > 3:
  217 + self._log_string(gammalib.NORMAL, '*** EVP file "%s" has inconsistent '
  218 + 'duration between events (%d days) and header '
  219 + '(%d days). Skip file.' % \
  220 + (evpname, duration_evt, duration_hdr))
  221 + continue
  222 + if duration_hdr > duration_evt:
  223 + self._log_string(gammalib.NORMAL, '--- EVP file "%s" header indicates '
  224 + 'duration of %d days, but events cover only %d days.' % \
  225 + (evpname, duration_hdr, duration_evt))
  226 + elif duration_hdr < duration_evt:
  227 + self._log_string(gammalib.NORMAL, '--- EVP file "%s" header indicates '
  228 + 'duration of %d days, but events cover %d days.' % \
  229 + (evpname, duration_hdr, duration_evt))
  230 +
  231 + # Set VP
  232 + evp_dict['vp'] = vpname
  233 +
  234 + # Append dictionary
  235 + evp_list.append(evp_dict)
  236 +
  237 + # Sort EVP list by start day
  238 + evp_list = sorted(evp_list, key=lambda k: k['tjd_start'])
  239 +
  240 + # Return EVP list
  241 + return evp_list
  242 +
  243 + def _get_evp_info(self, evpname):
  244 + """
  245 + Build EVP database entry from EVP file
  246 +
  247 + Parameters
  248 + ----------
  249 + evpname : str
  250 + EVP filename
  251 +
  252 + Returns
  253 + -------
  254 + evp_dict : dict
  255 + EVP database entry
  256 + """
  257 + # Open EVP file
  258 + fits = gammalib.GFits(evpname)
  259 +
  260 + # Get event table
  261 + events = fits.table(1)
  262 + nevents = events.nrows()
  263 +
  264 + # Extract attributes for EVP file
  265 + obs_id = events.real('OBS_ID')
  266 + object = events.string('OBJECT')
  267 + visday = events.integer('VISDAY')
  268 + vistim = events.integer('VISTIM')
  269 + vieday = events.integer('VIEDAY')
  270 + vietim = events.integer('VIETIM')
  271 + glon_scz = events.real('GLON_SCZ')
  272 + glat_scz = events.real('GLAT_SCZ')
  273 + verno = events.integer('DSD_REP')
  274 + quality = events.integer('DSD_QUA')
  275 + if events.has_card('TITLE'):
  276 + title = events.string('TITLE')
  277 + elif events.has_card('TITLE1') and events.has_card('TITLE2'):
  278 + title = events.string('TITLE1') + events.string('TITLE2')
  279 + dsd_id = '%s-%s-%d' % (events.string('DSD_SB'), events.string('DSTID'),
  280 + events.integer('DSD_NO'))
  281 +
  282 + # Get time columns
  283 + tjd = events['TJD']
  284 + tics = events['TICS']
  285 +
  286 + # Get start and end time. Catch any errors since files may be corrupt.
  287 + # In case of errors, use times from file header.
  288 + try:
  289 + tjd_start = tjd[0]
  290 + tics_start = tics[0]
  291 + tjd_stop = tjd[nevents-1]
  292 + tics_stop = tics[nevents-1]
  293 + except:
  294 + self._log_string(gammalib.NORMAL, '*** Corrupt EVP file "%s" encountered. '
  295 + 'Set EVP quality to -200.' % (evpname))
  296 + tjd_start = visday
  297 + tics_start = vistim
  298 + tjd_stop = vieday
  299 + tics_stop = vietim
  300 + quality = -200
  301 +
  302 + # Compute GTime of first and last event
  303 + tstart = gammalib.com_time(tjd_start, tics_start)
  304 + tstop = gammalib.com_time(tjd_stop, tics_stop)
  305 +
  306 + # Check that start time is before stop time
  307 + if tstart >= tstop:
  308 + quality = -200
  309 + gti = gammalib.GGti()
  310 + self._log_string(gammalib.NORMAL, '*** Corrupt EVP file "%s" encountered. '
  311 + 'Start time %s is later than stop time %s. Set EVP '
  312 + 'quality to -200.' % (evpname, tstart.utc(), tstop.utc()))
  313 +
  314 + # Set GTI for EVP
  315 + else:
  316 + gti = gammalib.GGti(tstart, tstop)
  317 +
  318 + # Build dictionary
  319 + evp_dict = {'filename': evpname,
  320 + 'nevents': nevents,
  321 + 'obs_id': obs_id,
  322 + 'dsd_id': dsd_id,
  323 + 'object': object,
  324 + 'glon_scz': glon_scz,
  325 + 'glat_scz': glat_scz,
  326 + 'date_start': tstart.utc(),
  327 + 'date_stop': tstop.utc(),
  328 + 'tjd_start': tjd_start,
  329 + 'tics_start': tics_start,
  330 + 'tjd_stop': tjd_stop,
  331 + 'tics_stop': tics_stop,
  332 + 'tstart': tstart,
  333 + 'tstop': tstop,
  334 + 'visday': visday,
  335 + 'vistim': vistim,
  336 + 'vieday': vieday,
  337 + 'vietim': vietim,
  338 + 'gti': gti,
  339 + 'quality': quality,
  340 + 'verno': verno,
  341 + 'title': title}
  342 +
  343 + # Close FITS file
  344 + fits.close()
  345 +
  346 + # Correct quality flags for some files according to database reports
  347 + if evp_dict['dsd_id'] == 'M-EVP-26342' or \
  348 + evp_dict['dsd_id'] == 'M-EVP-24955' or \
  349 + evp_dict['dsd_id'] == 'M-EVP-24963':
  350 + self._log_string(gammalib.NORMAL, '*** Quality flag for "%s" changed '
  351 + 'from %s to 50 according to database reports.' %
  352 + (evp_dict['filename'], evp_dict['quality']))
  353 + evp_dict['quality'] = 50
  354 +
  355 + # Return information
  356 + return evp_dict
  357 +
  358 + def _get_oad_list(self):
  359 + """
  360 + Get OAD list
  361 +
  362 + Returns
  363 + -------
  364 + oad_list : list of dict
  365 + OAD database
  366 + """
  367 + # Initialse OAD list
  368 + oad_list = []
  369 +
  370 + # Get available phases
  371 + phases = self._get_phases()
  372 +
  373 + # Loop over phases
  374 + for phase in phases:
  375 +
  376 + # Get available VPs
  377 + vps = glob.glob('%s/vp*' % phase)
  378 + vps.sort()
  379 +
  380 + # Loop over VPs
  381 + for vp in vps:
  382 +
  383 + # Get OADs for VP
  384 + oads = self._get_oads_for_vp(vp)
  385 +
  386 + # Append OADs to list
  387 + oad_list.extend(oads)
  388 +
  389 + # Sort OAD list by start day
  390 + oad_list = sorted(oad_list, key=lambda k: k['tjd_start'])
  391 +
  392 + # Return OAD list
  393 + return oad_list
  394 +
  395 + def _get_oads_for_vp(self, vp):
  396 + """
  397 + Build OAD database entries for Viewing Period
  398 +
  399 + Parameters
  400 + ----------
  401 + vp : str
  402 + VP path
  403 +
  404 + Returns
  405 + -------
  406 + oad_list : list of dict
  407 + OAD database entries
  408 + """
  409 + # Get archive name and expanded path
  410 + archive = self['archive'].string()
  411 + archive_path = gammalib.expand_env(archive)
  412 +
  413 + # Get VP name
  414 + vpname = os.path.basename(vp)
  415 +
  416 + # Initialse OAD list
  417 + oad_list = []
  418 +
  419 + # Get OAD files
  420 + oads = glob.glob('%s/*oad.fits*' % vp)
  421 +
  422 + # Loop over all OAD files
  423 + for oad in oads:
  424 +
  425 + # Set OAD filename
  426 + oadname = archive+'/'+os.path.relpath(oad, archive_path).strip('.gz')
  427 +
  428 + # Try opening file
  429 + try:
  430 + fits = gammalib.GFits(oadname)
  431 + except:
  432 + self._log_string(gammalib.NORMAL, '*** Unable to open OAD file "%s". '
  433 + 'Skip file.' % (oadname))
  434 + continue
  435 +
  436 + # Get information from OAD file
  437 + oad_dict = self._get_oad_info(oadname)
  438 +
  439 + # Skip file if it has a negative quality
  440 + if oad_dict['quality'] < 0:
  441 + self._log_string(gammalib.NORMAL, '*** OAD file "%s" has quality %d. '
  442 + 'Skip file.' % (oadname, oad_dict['quality']))
  443 + continue
  444 +
  445 + # Set VP
  446 + oad_dict['vp'] = vpname
  447 +
  448 + # Append dictionary
  449 + oad_list.append(oad_dict)
  450 +
  451 + # Sort OAD list by start day
  452 + oad_list = sorted(oad_list, key=lambda k: k['tjd_start'])
  453 +
  454 + # Return OAD list
  455 + return oad_list
  456 +
  457 + def _get_oad_info(self, oadname):
  458 + """
  459 + Build OAD database entry from OAD file
  460 +
  461 + Parameters
  462 + ----------
  463 + oadname : str
  464 + OAD filename
  465 +
  466 + Returns
  467 + -------
  468 + oad_dict : dict
  469 + OAD database entry
  470 + """
  471 + # Open OAD file
  472 + fits = gammalib.GFits(oadname)
  473 +
  474 + # Get superpacket table
  475 + records = fits.table(1)
  476 + nrecords = records.nrows()
  477 +
  478 + # Extract attributes for OAD file
  479 + obs_id = records.real('OBS_ID')
  480 + object = records.string('OBJECT')
  481 + if records.has_card('GLON_SCZ'):
  482 + glon_scz = records.real('GLON_SCZ')
  483 + else:
  484 + glon_scz = -1000.0
  485 + if records.has_card('GLAT_SCZ'):
  486 + glat_scz = records.real('GLAT_SCZ')
  487 + else:
  488 + glat_scz = -1000.0
  489 + quality = records.integer('DSD_QUA')
  490 + dsd_id = '%s-%s-%d' % (records.string('DSD_SB'), records.string('DSTID'),
  491 + records.integer('DSD_NO'))
  492 +
  493 + # Get time columns
  494 + tjd = records['TJD']
  495 + tics = records['TICS']
  496 +
  497 + # Get start and end time. Catch any errors since files may be corrupt.
  498 + # In case of errors, use times from file header.
  499 + try:
  500 + tjd_start = tjd[0]
  501 + tics_start = tics[0]
  502 + tjd_stop = tjd[nrecords-1]
  503 + tics_stop = tics[nrecords-1]
  504 + except:
  505 + self._log_string(gammalib.NORMAL, '*** Corrupt OAD file "%s" encountered. '
  506 + 'Set OAD quality to -200.' % (evpname))
  507 + tjd_start = records.integer('VISDAY')
  508 + tics_start = records.integer('VISTIM')
  509 + tjd_stop = records.integer('VIEDAY')
  510 + tics_stop = records.integer('VIETIM')
  511 + quality = -200
  512 +
  513 + # Compute GTime of first and last superpacket
  514 + tstart = gammalib.com_time(tjd_start, tics_start)
  515 + tstop = gammalib.com_time(tjd_stop, tics_stop)
  516 +
  517 + # Set GTI for OAD
  518 + gti = gammalib.GGti(tstart, tstop)
  519 +
  520 + # Build dictionary
  521 + oad_dict = {'filename': oadname,
  522 + 'superpackets': nrecords,
  523 + 'obs_id': obs_id,
  524 + 'dsd_id': dsd_id,
  525 + 'object': object,
  526 + 'glon_scz': glon_scz,
  527 + 'glat_scz': glat_scz,
  528 + 'date_start': tstart.utc(),
  529 + 'date_stop': tstop.utc(),
  530 + 'tjd_start': tjd_start,
  531 + 'tics_start': tics_start,
  532 + 'tjd_stop': tjd_stop,
  533 + 'tics_stop': tics_stop,
  534 + 'tstart': tstart,
  535 + 'tstop': tstop,
  536 + 'gti': gti,
  537 + 'quality': quality}
  538 +
  539 + # Close FITS file
  540 + fits.close()
  541 +
  542 + # Return information
  543 + return oad_dict
  544 +
  545 + def _get_tim_list(self):
  546 + """
  547 + Get TIM list
  548 +
  549 + Returns
  550 + -------
  551 + tim_list : list of dict
  552 + TIM database
  553 + """
  554 + # Initialse TIM list
  555 + tim_list = []
  556 +
  557 + # Get available phases
  558 + phases = self._get_phases()
  559 +
  560 + # Loop over phases
  561 + for phase in phases:
  562 +
  563 + # Get available VPs
  564 + vps = glob.glob('%s/vp*' % phase)
  565 + vps.sort()
  566 +
  567 + # Loop over VPs
  568 + for vp in vps:
  569 +
  570 + # Get TIMs for VP
  571 + tims = self._get_tims_for_vp(vp)
  572 +
  573 + # Append TIMs to list
  574 + tim_list.extend(tims)
  575 +
  576 + # Sort TIM list by start day
  577 + tim_list = sorted(tim_list, key=lambda k: k['tjd_start'])
  578 +
  579 + # Return TIM list
  580 + return tim_list
  581 +
  582 + def _get_tims_for_vp(self, vp):
  583 + """
  584 + Build TIM database entries for Viewing Period
  585 +
  586 + Parameters
  587 + ----------
  588 + vp : str
  589 + VP path
  590 +
  591 + Returns
  592 + -------
  593 + tim_list : list of dict
  594 + TIM database entries
  595 + """
  596 + # Get archive name and expanded path
  597 + archive = self['archive'].string()
  598 + archive_path = gammalib.expand_env(archive)
  599 +
  600 + # Get VP name
  601 + vpname = os.path.basename(vp)
  602 +
  603 + # Initialse TIM list
  604 + tim_list = []
  605 +
  606 + # Get TIM files
  607 + tims = glob.glob('%s/*tim.fits*' % vp)
  608 +
  609 + # Loop over all TIM files
  610 + for tim in tims:
  611 +
  612 + # Set TIM filename
  613 + timname = archive+'/'+os.path.relpath(tim, archive_path).strip('.gz')
  614 +
  615 + # Try opening file
  616 + try:
  617 + fits = gammalib.GFits(timname)
  618 + except:
  619 + self._log_string(gammalib.NORMAL, '*** Unable to open TIM file "%s". '
  620 + 'Skip file.' % (timname))
  621 + continue
  622 +
  623 + # Get information from TIM file
  624 + tim_dict = self._get_tim_info(timname)
  625 +
  626 + # Skip file if it has a negative quality
  627 + if tim_dict['quality'] < 0:
  628 + self._log_string(gammalib.NORMAL, '*** TIM file "%s" has quality %d. '
  629 + 'Skip file.' % (timname, tim_dict['quality']))
  630 + continue
  631 +
  632 + # Set VP
  633 + tim_dict['vp'] = vpname
  634 +
  635 + # Append dictionary
  636 + tim_list.append(tim_dict)
  637 +
  638 + # Sort OAD list by start day
  639 + tim_list = sorted(tim_list, key=lambda k: k['tjd_start'])
  640 +
  641 + # Return TIM list
  642 + return tim_list
  643 +
  644 + def _get_tim_info(self, timname):
  645 + """
  646 + Build TIM database entry from TIM file
  647 +
  648 + Parameters
  649 + ----------
  650 + timname : str
  651 + TIM filename
  652 +
  653 + Returns
  654 + -------
  655 + tim_dict : dict
  656 + TIM database entry
  657 + """
  658 + # Open TIM file
  659 + fits = gammalib.GFits(timname)
  660 +
  661 + # Get superpacket table
  662 + records = fits.table(1)
  663 + nrecords = records.nrows()
  664 +
  665 + # Extract attributes for TIM file
  666 + obs_id = records.real('OBS_ID')
  667 + object = records.string('OBJECT')
  668 + visday = records.integer('VISDAY')
  669 + vistim = records.integer('VISTIM')
  670 + vieday = records.integer('VIEDAY')
  671 + vietim = records.integer('VIETIM')
  672 + glon_scz = records.real('GLON_SCZ')
  673 + glat_scz = records.real('GLAT_SCZ')
  674 + quality = records.integer('DSD_QUA')
  675 + dsd_id = '%s-%s-%d' % (records.string('DSD_SB'), records.string('DSTID'),
  676 + records.integer('DSD_NO'))
  677 +
  678 + # Get time columns
  679 + start_tjd = records['START_TJD']
  680 + start_tic = records['START_TIC']
  681 + end_tjd = records['END_TJD']
  682 + end_tic = records['END_TIC']
  683 +
  684 + # Get start and end time. Catch any errors since files may be corrupt.
  685 + # In case of errors, use times from file header.
  686 + try:
  687 + tjd_start = start_tjd[0]
  688 + tics_start = start_tic[0]
  689 + tjd_stop = end_tjd[nrecords-1]
  690 + tics_stop = end_tic[nrecords-1]
  691 + except:
  692 + self._log_string(gammalib.NORMAL, '*** Corrupt TIM file "%s" encountered. '
  693 + 'Set TIM quality to -200.' % (evpname))
  694 + tjd_start = visday
  695 + tics_start = vistim
  696 + tjd_stop = vieday
  697 + tics_stop = vietim
  698 + quality = -200
  699 +
  700 + # Close FITS file
  701 + fits.close()
  702 +
  703 + # Compute GTime of first and last superpacket
  704 + tstart = gammalib.com_time(tjd_start, tics_start)
  705 + tstop = gammalib.com_time(tjd_stop, tics_stop)
  706 +
  707 + # Set GTI for TIM
  708 + gti = gammalib.GGti(tstart, tstop)
  709 +
  710 + # Get telapse and ontime
  711 + tim = gammalib.GCOMTim(timname)
  712 + telapse = tim.gti().telapse()
  713 + ontime = tim.gti().ontime()
  714 +
  715 + # Build dictionary
  716 + tim_dict = {'filename': timname,
  717 + 'intervals': nrecords,
  718 + 'obs_id': obs_id,
  719 + 'dsd_id': dsd_id,
  720 + 'object': object,
  721 + 'glon_scz': glon_scz,
  722 + 'glat_scz': glat_scz,
  723 + 'date_start': tstart.utc(),
  724 + 'date_stop': tstop.utc(),
  725 + 'tjd_start': tjd_start,
  726 + 'tics_start': tics_start,
  727 + 'tjd_stop': tjd_stop,
  728 + 'tics_stop': tics_stop,
  729 + 'visday': visday,
  730 + 'vistim': vistim,
  731 + 'vieday': vieday,
  732 + 'vietim': vietim,
  733 + 'tstart': tstart,
  734 + 'tstop': tstop,
  735 + 'telapse': telapse,
  736 + 'ontime': ontime,
  737 + 'gti': gti,
  738 + 'quality': quality}
  739 +
  740 + # Return information
  741 + return tim_dict
  742 +
  743 + def _get_xml_list(self, evp_list, tim_list, oad_list):
  744 + """
  745 + Build XML database from EVP, TIM and OAD database
  746 +
  747 + Parameters
  748 + ----------
  749 + evp_list : list of dict
  750 + EVP database
  751 + tim_list : list of dict
  752 + TIM database
  753 + oad_list : list of dict
  754 + OAD database
  755 +
  756 + Returns
  757 + -------
  758 + xml_list : list of dict
  759 + XML database
  760 + """
  761 + # Initialse XML list
  762 + xml_list = []
  763 +
  764 + # Get filtered TIM list
  765 + filtered_tim_list = self._filter_tim_list(tim_list)
  766 +
  767 + # Get available phases
  768 + phases = self._get_phases()
  769 +
  770 + # Loop over phases
  771 + for phase in phases:
  772 +
  773 + # Get available VPs
  774 + vps = glob.glob('%s/vp*' % phase)
  775 + vps.sort()
  776 +
  777 + # Loop over VPs
  778 + for vp in vps:
  779 +
  780 + # Get XML for VP
  781 + xml = self._get_xml_for_vp(vp, evp_list, filtered_tim_list, oad_list)
  782 +
  783 + # Append XML to list
  784 + if xml is not None:
  785 + xml_list.append(xml)
  786 +
  787 + # Sort XML list by start day
  788 + xml_list = sorted(xml_list, key=lambda k: k['tjd_start'])
  789 +
  790 + # Return XML list
  791 + return xml_list
  792 +
  793 + def _filter_tim_list(self, tim_list):
  794 + """
  795 + Removes all invalid TIMs from TIM database
  796 +
  797 + For a given VP only the TIMs with the maximum quality will be kept. This
  798 + is needed since some of the TIM dataset cover periods well beyond a given
  799 + VP, and these TIMs are marked with data quality 135. On the other hand,
  800 + there are VPs with valid TIMs of data quality 135, hence we cannot simply
  801 + filter on data quality.
  802 +
  803 + Parameters
  804 + ----------
  805 + tim_list : list of dict
  806 + TIM database
  807 +
  808 + Returns
  809 + -------
  810 + filtered_tim_list : list of dict
  811 + Filtered TIM database
  812 + """
  813 + # Initialise filtered TIM list
  814 + filtered_tim_list = []
  815 +
  816 + # Get available phases
  817 + phases = self._get_phases()
  818 +
  819 + # Loop over phases
  820 + for phase in phases:
  821 +
  822 + # Get available VPs
  823 + vps = glob.glob('%s/vp*' % phase)
  824 + vps.sort()
  825 +
  826 + # Loop over VPs
  827 + for vp in vps:
  828 +
  829 + # Get VP name
  830 + vpname = os.path.basename(vp)
  831 +
  832 + # Get all TIMs for VP
  833 + tims = filter(lambda k: k['vp'] == vpname, tim_list)
  834 +
  835 + # Continue if list is empty
  836 + if len(tims) < 1:
  837 + continue
  838 +
  839 + # Get maximum quality of all TIMs in VP
  840 + quality = max(tim['quality'] for tim in tims)
  841 +
  842 + # Get all TIMs for VP with maximum quality
  843 + best_tims = filter(lambda k: k['quality'] == quality, tims)
  844 +
  845 + # Append these TIMs to the filtered list
  846 + filtered_tim_list.extend(best_tims)
  847 +
  848 + # Return
  849 + return filtered_tim_list
  850 +
  851 + def _get_xml_for_vp(self, vp, evp_list, tim_list, oad_list):
  852 + """
  853 + Build XML database entry for Viewing Period
  854 +
  855 + Parameters
  856 + ----------
  857 + vp : str
  858 + Viewing Period directory
  859 + evp_list : list of dict
  860 + EVP database
  861 + tim_list : list of dict
  862 + TIM database
  863 + oad_list : list of dict
  864 + OAD database
  865 +
  866 + Returns
  867 + -------
  868 + xml_dict : dictionary
  869 + XML database entry
  870 + """
  871 + # Get database directory
  872 + path = self['dbase'].string()
  873 +
  874 + # Set VP name
  875 + vpname = os.path.basename(vp)
  876 +
  877 + # Get EVP for VP
  878 + evp = self._get_evp(vp, evp_list)
  879 + if evp is None:
  880 + if vpname == 'vp0225_0' or \
  881 + vpname == 'vp0303_5' or \
  882 + vpname == 'vp0308_3' or \
  883 + vpname == 'vp0311_3' or \
  884 + vpname == 'vp0617_2' or \
  885 + vpname == 'vp0617_4' or \
  886 + vpname == 'vp0617_6' or \
  887 + vpname == 'vp0619_5':
  888 + self._log_string(gammalib.NORMAL, '--- No EVP file found for "%s due '
  889 + 'to Reboost". Skip VP.' % (vp))
  890 + elif vpname == 'vp0229_3':
  891 + self._log_string(gammalib.NORMAL, '--- No EVP file found for "%s due '
  892 + 'to Perseids". Skip VP.' % (vp))
  893 + elif vpname == 'vp0833_7':
  894 + self._log_string(gammalib.NORMAL, '--- No EVP file found for "%s due '
  895 + 'to Leonid Shower". Skip VP.' % (vp))
  896 + else:
  897 + self._log_string(gammalib.NORMAL, '*** No EVP file found for "%s". '
  898 + 'Skip VP.' % (vp))
  899 + return None
  900 +
  901 + # Get TIM for EVP
  902 + tim = self._get_tim(evp, tim_list)
  903 + if tim is None:
  904 + self._log_string(gammalib.NORMAL, '*** No TIM file found for "%s". '
  905 + 'Skip VP.' % (vp))
  906 + return None
  907 +
  908 + # Get OADs for EVP
  909 + oads = self._get_oads(evp, oad_list)
  910 + if len(oads) == 0:
  911 + self._log_string(gammalib.NORMAL, '*** No OAD files found for "%s". '
  912 + 'Skip VP.' % (vp))
  913 + return None
  914 +
  915 + # Get oadnames and oad qualities
  916 + oadnames = [oad['filename'] for oad in oads]
  917 + oadqualities = [oad['quality'] for oad in oads]
  918 +
  919 + # Build dictionary
  920 + xml_dict = {'vp': vpname,
  921 + 'evpname': evp['filename'],
  922 + 'timname': tim['filename'],
  923 + 'oadnames': oadnames,
  924 + 'xmlname': '%s/xml/%s.xml' % (path, vpname),
  925 + 'nevents': evp['nevents'],
  926 + 'obs_id': evp['obs_id'],
  927 + 'object': evp['object'],
  928 + 'glon_scz': evp['glon_scz'],
  929 + 'glat_scz': evp['glat_scz'],
  930 + 'date_start': tim['date_start'],
  931 + 'date_stop': tim['date_stop'],
  932 + 'tjd_start': tim['tjd_start'],
  933 + 'tics_start': tim['tics_start'],
  934 + 'tjd_stop': tim['tjd_stop'],
  935 + 'tics_stop': tim['tics_stop'],
  936 + 'evp_tstart': evp['tstart'],
  937 + 'evp_tstop': evp['tstop'],
  938 + 'visday': tim['visday'],
  939 + 'vistim': tim['vistim'],
  940 + 'vieday': tim['vieday'],
  941 + 'vietim': tim['vietim'],
  942 + 'telapse': tim['telapse'],
  943 + 'ontime': tim['ontime'],
  944 + 'evp_verno': evp['verno'],
  945 + 'evp_quality': evp['quality'],
  946 + 'tim_quality': tim['quality'],
  947 + 'oad_quality': oadqualities}
  948 +
  949 + # Return XML list
  950 + return xml_dict
  951 +
  952 + def _get_evp(self, vp, evp_list):
  953 + """
  954 + Get EVP database entry for Viewing Period
  955 +
  956 + Parameters
  957 + ----------
  958 + vp : str
  959 + Viewing Period directory
  960 + evp_list : list of dict
  961 + EVP database
  962 +
  963 + Returns
  964 + -------
  965 + best_evp : dict
  966 + EVP database entry
  967 + """
  968 + # Get quality flag
  969 + qf = self['quality'].integer()
  970 +
  971 + # Get VP name
  972 + vpname = os.path.basename(vp)
  973 +
  974 + # Get all EVPs for VP
  975 + evps = filter(lambda k: k['vp'] == vpname, evp_list)
  976 +
  977 + # Initialise best EVP
  978 + best_verno_evp = None
  979 + best_quality_evp = None
  980 + best_verno = 0
  981 + best_quality = 0
  982 +
  983 + # Initialise choosen flag
  984 + choosen = False
  985 +
  986 + # First find EVP with target QF and highest version number
  987 + for evp in evps:
  988 + if evp['quality'] == qf:
  989 + if evp['verno'] > best_verno:
  990 + choosen = True
  991 + best_evp = evp
  992 + best_verno = evp['verno']
  993 +
  994 + # Continue only if no EVP was found
  995 + if not choosen:
  996 +
  997 + # Find EVP with highest version number
  998 + for evp in evps:
  999 + if evp['verno'] > best_verno:
  1000 + best_verno = evp['verno']
  1001 + best_verno_evp = evp
  1002 +
  1003 + # Find EVP with best quality
  1004 + for evp in evps:
  1005 + if evp['quality'] > best_quality:
  1006 + best_quality = evp['quality']
  1007 + best_quality_evp = evp
  1008 +
  1009 + # Initialise best EVP with best quality EVP
  1010 + best_evp = best_quality_evp
  1011 +
  1012 + # If highest version number is not the EVP with the best quality then
  1013 + # select the EVP with the highest version number
  1014 + if best_verno_evp != None and best_quality_evp != None:
  1015 + if best_verno_evp['filename'] != best_quality_evp['filename']:
  1016 + if best_quality_evp['verno'] < best_verno:
  1017 +
  1018 + # Signal issue
  1019 + self._log_string(gammalib.NORMAL, '*** EVP "%s" with best '
  1020 + 'quality has version %d which is not the '
  1021 + 'highest version %d.' %
  1022 + (best_quality_evp['filename'],
  1023 + best_quality_evp['verno'],
  1024 + best_verno))
  1025 +
  1026 + # If there are several EVPs with highest version then select
  1027 + # the one with the best quality
  1028 + evps = filter(lambda k: k['vp'] == vpname and k['verno'] == best_verno, evp_list)
  1029 + if len(evps) > 1:
  1030 + self._log_string(gammalib.NORMAL, '--- Several EVPs for %s '
  1031 + 'have highest version' % (vpname))
  1032 + best = 0
  1033 + for evp in evps:
  1034 + if evp['quality'] > best:
  1035 + best = evp['quality']
  1036 + best_evp = evp
  1037 + for evp in evps:
  1038 + if best_evp['filename'] == evp['filename']:
  1039 + choose = '+'
  1040 + choosen = True
  1041 + else:
  1042 + choose = ' '
  1043 + self._log_string(gammalib.NORMAL, ' %s "%s" (verno=%d, '
  1044 + 'quality=%d)' %
  1045 + (choose,
  1046 + evp['filename'],
  1047 + evp['verno'],
  1048 + evp['quality']))
  1049 +
  1050 + # ... otherwise use the EVP with the best version number
  1051 + else:
  1052 + best_evp = best_verno_evp
  1053 + choosen = True
  1054 + self._log_string(gammalib.NORMAL, '--- Choose "%s" (verno=%d, '
  1055 + 'quality=%d)' %
  1056 + (best_evp['filename'],
  1057 + best_evp['verno'],
  1058 + best_evp['quality']))
  1059 +
  1060 + # Continue only if no file was choosen
  1061 + if not choosen:
  1062 +
  1063 + # If there are several EVPs with the best quality then select the one
  1064 + # with the highest version
  1065 + evps = filter(lambda k: k['vp'] == vpname and k['quality'] == best_quality, evp_list)
  1066 + if len(evps) > 1:
  1067 + self._log_string(gammalib.NORMAL, '--- Several EVPs for %s have '
  1068 + 'best quality' % (vpname))
  1069 + best = 0
  1070 + for evp in evps:
  1071 + if evp['verno'] > best:
  1072 + best = evp['verno']
  1073 + best_evp = evp
  1074 + for evp in evps:
  1075 + if best_evp['filename'] == evp['filename']:
  1076 + choose = '+'
  1077 + else:
  1078 + choose = ' '
  1079 + self._log_string(gammalib.NORMAL, ' %s "%s" (verno=%d, '
  1080 + 'quality=%d)' %
  1081 + (choose,
  1082 + evp['filename'],
  1083 + evp['verno'],
  1084 + evp['quality']))
  1085 +
  1086 + # Return best EVP
  1087 + return best_evp
  1088 +
  1089 + def _get_tim(self, evp, tim_list):
  1090 + """
  1091 + Get TIM database entry for EVP database entry
  1092 +
  1093 + Parameters
  1094 + ----------
  1095 + evp : dict
  1096 + EVP dictionary
  1097 + tim_list : list of dict
  1098 + TIM database
  1099 +
  1100 + Returns
  1101 + -------
  1102 + best_tim : dict
  1103 + TIM database entry for EVP database entry
  1104 + """
  1105 + # Initialise maximum overlap
  1106 + max_overlap = 0.0
  1107 + min_dtabs = 1.0e30
  1108 + best_tim = None
  1109 +
  1110 + # Loop over all TIMs
  1111 + for tim in tim_list:
  1112 +
  1113 + # Compute time difference to EVP start time and stop time
  1114 + dtstart = evp['tstart'] - tim['tstart']
  1115 + dtstop = evp['tstop'] - tim['tstop']
  1116 +
  1117 + # Compute absolute time difference
  1118 + dtabs = abs(dtstart) + abs(dtstop)
  1119 +
  1120 + # Keep smallest time difference
  1121 + if dtabs < min_dtabs:
  1122 + min_dtabs = dtabs
  1123 + best_tim = tim
  1124 +
  1125 + # Return best TIM
  1126 + return best_tim
  1127 +
  1128 + def _get_oads(self, evp, oad_list):
  1129 + """
  1130 + Get OAD database entry for EVP database entry
  1131 +
  1132 + Parameters
  1133 + ----------
  1134 + evp : dict
  1135 + EVP dictionary
  1136 + oad_list : list of dict
  1137 + OAD database
  1138 +
  1139 + Returns
  1140 + -------
  1141 + oads : list of dict
  1142 + OAD database entries for EVP database entry
  1143 + """
  1144 + # Initialise list of OADs
  1145 + oads = []
  1146 +
  1147 + # Loop over all OADs
  1148 + for oad in oad_list:
  1149 +
  1150 + # Compute overlap between EVP and OAD
  1151 + overlap = evp['gti'].overlap(oad['tstart'], oad['tstop'])
  1152 +
  1153 + # If we have overlap then check if the file already exists
  1154 + if overlap > 0.0:
  1155 +
  1156 + # First check if OAD is already in OAD list. If this is the case,
  1157 + # and the current OAD corresponds to the VP of the EVP, then replace
  1158 + # the existing OAD by the current OAD since we prefer having the
  1159 + # OAD in the same folder as the EVP. Otherwise we keep the
  1160 + exists = False
  1161 + for i, oad_exist in enumerate(oads):
  1162 + if oad_exist['dsd_id'] == oad['dsd_id']:
  1163 + exists = True
  1164 + if oad['vp'] == evp['vp']:
  1165 + oads[i] = oad
  1166 + break
  1167 + else:
  1168 + break
  1169 +
  1170 + # If the OAD does not yet exist then append it now
  1171 + if not exists:
  1172 + oads.append(oad)
  1173 +
  1174 + # Sort OAD list by start day
  1175 + oads = sorted(oads, key=lambda k: k['tjd_start'])
  1176 +
  1177 + # Return OADs
  1178 + return oads
  1179 +
  1180 + def _check_xml_list(self, xml_list):
  1181 + """
  1182 + Check XML list against COMPASS database
  1183 +
  1184 + Returns
  1185 + -------
  1186 + xml_list : list of dict
  1187 + XML database with mode and status information
  1188 + """
  1189 + # Load COMPASS database
  1190 + database = self._get_compass_database()
  1191 +
  1192 + # Loop over XML list
  1193 + for xml_dict in xml_list:
  1194 +
  1195 + # Add mode and status to XML list
  1196 + xml_dict['mode'] = '???'
  1197 + xml_dict['status'] = 'unknown'
  1198 +
  1199 + # Continue only if there are entries
  1200 + if len(database) > 0:
  1201 +
  1202 + # Find VP in COMPASS database
  1203 + compass = None
  1204 + for entry in database:
  1205 + if entry['vp'] == xml_dict['vp']:
  1206 + compass = entry
  1207 + break
  1208 +
  1209 + # Signal if VP was not found in COMPASS database
  1210 + if compass == None:
  1211 + xml_dict['status'] = 'not found'
  1212 + self._log_string(gammalib.NORMAL, '--- VP "%s" not found in COMPASS '
  1213 + 'database.' % (xml_dict['vp']))
  1214 + continue
  1215 +
  1216 + # Set mode and status
  1217 + xml_dict['mode'] = compass['mode']
  1218 + xml_dict['status'] = 'found'
  1219 +
  1220 + # Check EVP filename
  1221 + evp = os.path.basename(xml_dict['evpname'])[1:].rstrip('_evp.fits')
  1222 + if evp != compass['evp']:
  1223 + self._log_string(gammalib.NORMAL, '--- VP "%s": EVP "%s" in COMPASS '
  1224 + 'database differs from EVP "%s" in XML database.' %
  1225 + (xml_dict['vp'], compass['evp'], evp))
  1226 +
  1227 + # Check TIM filename
  1228 + tim = os.path.basename(xml_dict['timname'])[1:].rstrip('_tim.fits')
  1229 + if tim != compass['tim']:
  1230 + self._log_string(gammalib.NORMAL, '--- VP "%s": TIM "%s" in COMPASS '
  1231 + 'database differs from TIM "%s" in XML database.' %
  1232 + (xml_dict['vp'], compass['tim'], tim))
  1233 +
  1234 + # Return
  1235 + return xml_list
  1236 +
  1237 + def _get_compass_database(self):
  1238 + """
  1239 + Get COMPASS database
  1240 +
  1241 + Returns
  1242 + -------
  1243 + database : list of dict
  1244 + COMPASS database
  1245 + """
  1246 + # Initialise database
  1247 + database = []
  1248 +
  1249 + # Continue only if database was specified
  1250 + if self['refdata'].is_valid():
  1251 +
  1252 + # Read ASCII file
  1253 + f = open(gammalib.expand_env(self['refdata'].filename().url()), 'r')
  1254 + for line in f:
  1255 +
  1256 + # Skip comment lines
  1257 + if line[0:1] == 'C':
  1258 + continue
  1259 +
  1260 + # Format VP
  1261 + vp = ('vp%06.1f' % (float(line[1:6]))).replace('.', '_')
  1262 +
  1263 + # Extract fields
  1264 + record = {'vp': vp,
  1265 + 'evp': gammalib.strip_whitespace(line[13:19]),
  1266 + 'tim': gammalib.strip_whitespace(line[21:27]),
  1267 + 'oad': gammalib.strip_whitespace(line[29:35]),
  1268 + 'mode': gammalib.strip_whitespace(line[66:70])}
  1269 +
  1270 + # Append record to database
  1271 + database.append(record)
  1272 +
  1273 + # Return database
  1274 + return database
  1275 +
  1276 + def _create_dbase_dirs(self):
  1277 + """
  1278 + Create database directory structure
  1279 + """
  1280 + # Get database directory
  1281 + path = self['dbase'].string()
  1282 +
  1283 + # Create database directory
  1284 + try:
  1285 + os.makedirs(gammalib.expand_env(path))
  1286 + except OSError:
  1287 + self._log_string(gammalib.NORMAL, '--- Directory "%s" exists' % path)
  1288 + try:
  1289 + os.makedirs('%s/xml' % gammalib.expand_env(path))
  1290 + except OSError:
  1291 + self._log_string(gammalib.NORMAL, '--- Directory "%s/xml" exists' % path)
  1292 + try:
  1293 + os.makedirs('%s/tim' % gammalib.expand_env(path))
  1294 + except OSError:
  1295 + self._log_string(gammalib.NORMAL, '--- Directory "%s/tim" exists' % path)
  1296 +
  1297 + # Return
  1298 + return
  1299 +
  1300 + def _build_tim_file(self, xml_dict):
  1301 + """
  1302 + Build TIM file for XML database entry
  1303 +
  1304 + Parameters
  1305 + ----------
  1306 + xml_dict : dict
  1307 + XML database entry
  1308 + """
  1309 + # Get database directory
  1310 + path = self['dbase'].string()
  1311 +
  1312 + # Get TIM from database
  1313 + tim = gammalib.GCOMTim(xml_dict['timname'])
  1314 +
  1315 + # Extract GTIs
  1316 + gti_heasarc = tim.gti()
  1317 +
  1318 + # Get tstart and tstop from first and last event
  1319 + tstart = xml_dict['evp_tstart']
  1320 + tstop = xml_dict['evp_tstop']
  1321 +
  1322 + # Kluge for VP vp0411_1 which overlaps with vp0411_5
  1323 + if xml_dict['vp'] == 'vp0411_1':
  1324 + tstop = gammalib.com_time(9769, 444782835)
  1325 + self._log_string(gammalib.NORMAL, '*** Stop date for "%s" changed from %s '
  1326 + 'to %s.' %
  1327 + (xml_dict['vp'], xml_dict['evp_tstop'].utc(), tstop.utc()))
  1328 +
  1329 + # Initialise GTI
  1330 + gti = gammalib.GGti()
  1331 +
  1332 + # Select GTIs
  1333 + search_start = True
  1334 + for i in range(len(gti_heasarc)):
  1335 +
  1336 + # Case A: search for start
  1337 + if search_start:
  1338 +
  1339 + # If event start time is after GTI stop time then skip this GTI
  1340 + if tstart > gti_heasarc.tstop(i):
  1341 + continue
  1342 +
  1343 + # ... otherwise stop searching for start and set first GTI
  1344 + else:
  1345 + # Stop searching for start
  1346 + search_start = False
  1347 +
  1348 + # Get first GTI
  1349 + if tstart < gti_heasarc.tstart(i):
  1350 + t_start = gti_heasarc.tstart(i)
  1351 + t_stop = gti_heasarc.tstop(i)
  1352 + else:
  1353 + t_start = tstart
  1354 + t_stop = gti_heasarc.tstop(i)
  1355 +
  1356 + # Append GTI
  1357 + gti.append(t_start, t_stop)
  1358 +
  1359 + # Case B: search for stop
  1360 + else:
  1361 +
  1362 + # If event stop time after GTI stop time then and this GTI
  1363 + if tstop > gti_heasarc.tstop(i):
  1364 + t_start = gti_heasarc.tstart(i)
  1365 + t_stop = gti_heasarc.tstop(i)
  1366 + gti.append(t_start, t_stop)
  1367 +
  1368 + # ... otherwise, if the event stop time is before the GTI start
  1369 + # time then stop appending GTIs
  1370 + elif tstop < gti_heasarc.tstart(i):
  1371 + break
  1372 +
  1373 + # ... otherwise append last GTI
  1374 + else:
  1375 + t_start = gti_heasarc.tstart(i)
  1376 + t_stop = tstop
  1377 + gti.append(t_start, t_stop)
  1378 + break
  1379 +
  1380 + # Set new GTIs for TIM
  1381 + tim.gti(gti)
  1382 +
  1383 + # Set new TIM filename
  1384 + xml_dict['timname'] = '%s/tim/%s_tim.fits' % (path, xml_dict['vp'])
  1385 +
  1386 + # Save TIM
  1387 + tim.save(xml_dict['timname'], True)
  1388 +
  1389 + # Return
  1390 + return
  1391 +
  1392 + def _build_xml_file(self, xml_dict):
  1393 + """
  1394 + Build observation definition XML file for XML database entry
  1395 +
  1396 + Parameters
  1397 + ----------
  1398 + xml_dict : dict
  1399 + XML database entry
  1400 + """
  1401 + # Create XML object
  1402 + xml = gammalib.GXml()
  1403 +
  1404 + # Append observation list
  1405 + obslist = xml.append('observation_list title="observation list"')
  1406 +
  1407 + # Append observation
  1408 + obs = obslist.append('observation name="%s" id="%s" instrument="COM"' % \
  1409 + (xml_dict['object'], xml_dict['vp']))
  1410 +
  1411 + # Append EVP file
  1412 + obs.append('parameter name="EVP" file="%s"' % (xml_dict['evpname']))
  1413 +
  1414 + # Append TIM file
  1415 + obs.append('parameter name="TIM" file="%s"' % (xml_dict['timname']))
  1416 +
  1417 + # Append OAD files
  1418 + for oadname in xml_dict['oadnames']:
  1419 + obs.append('parameter name="OAD" file="%s"' % (oadname))
  1420 +
  1421 + # Save XML file
  1422 + xml.save(xml_dict['xmlname'])
  1423 +
  1424 + # Return
  1425 + return
  1426 +
  1427 + def _check_database(self, xml_list):
  1428 + """
  1429 + Check database
  1430 +
  1431 + Parameters
  1432 + ----------
  1433 + xml_list : list of dict
  1434 + XML database
  1435 + """
  1436 + # Determine number of database entries
  1437 + nrows = len(xml_list)
  1438 +
  1439 + # Check that stop dates are after start dates
  1440 + for i in range(nrows):
  1441 + tstart = gammalib.GTime()
  1442 + tstop = gammalib.GTime()
  1443 + tstart.utc(xml_list[i]['date_start'])
  1444 + tstop.utc(xml_list[i]['date_stop'])
  1445 + if tstart > tstop:
  1446 + self._log_string(gammalib.NORMAL, '*** Start date %s of VP %s later '
  1447 + 'than stop date %s.' %
  1448 + (xml_list[i]['date_start'],
  1449 + xml_list[i]['vp'],
  1450 + xml_list[i]['date_stop']))
  1451 +
  1452 + # Check that stop dates of VP is before start date of next VP. If this is
  1453 + # not the case, events may be double counted.
  1454 + for i in range(nrows-1):
  1455 + tstart = gammalib.GTime()
  1456 + tstop = gammalib.GTime()
  1457 + tstop.utc(xml_list[i]['date_stop'])
  1458 + tstart.utc(xml_list[i+1]['date_start'])
  1459 + if tstop >= tstart:
  1460 + self._log_string(gammalib.NORMAL, '*** Stop date %s of VP %s later '
  1461 + 'than start date %s of VP %s.' %
  1462 + (xml_list[i]['date_stop'],
  1463 + xml_list[i]['vp'],
  1464 + xml_list[i+1]['date_start'],
  1465 + xml_list[i+1]['vp']))
  1466 +
  1467 + # Return
  1468 + return
  1469 +
  1470 + def _build_database_file(self, evp_list, oad_list, tim_list, xml_list):
  1471 + """
  1472 + Build database file
  1473 +
  1474 + Parameters
  1475 + ----------
  1476 + evp_list : list of dict
  1477 + EVP database
  1478 + oad_list : list of dict
  1479 + OAD database
  1480 + tim_list : list of dict
  1481 + TIM database
  1482 + xml_list : list of dict
  1483 + XML database
  1484 + """
  1485 + # Setup filename
  1486 + filename = '%s/dbase.fits' % (self['dbase'].string())
  1487 +
  1488 + # Create FITS file
  1489 + fits = gammalib.GFits()
  1490 +
  1491 + # Get tables
  1492 + evp_table = self._build_evp_table(evp_list)
  1493 + oad_table = self._build_oad_table(oad_list)
  1494 + tim_table = self._build_tim_table(tim_list)
  1495 + xml_table = self._build_xml_table(xml_list)
  1496 +
  1497 + # Append tables
  1498 + fits.append(evp_table)
  1499 + fits.append(oad_table)
  1500 + fits.append(tim_table)
  1501 + fits.append(xml_table)
  1502 +
  1503 + # Save database
  1504 + fits.saveto(filename, True)
  1505 +
  1506 + # Close FITS file
  1507 + fits.close()
  1508 +
  1509 + # Return
  1510 + return
  1511 +
  1512 + def _build_evp_table(self, evp_list):
  1513 + """
  1514 + Build EVP FITS table
  1515 +
  1516 + Parameters
  1517 + ----------
  1518 + evp_list : list of dict
  1519 + EVP database
  1520 +
  1521 + Returns
  1522 + -------
  1523 + table : `~gammalib.GFitsBinTable`
  1524 + EVP FITS table
  1525 + """
  1526 + # Determine number of rows in table
  1527 + nrows = len(evp_list)
  1528 +
  1529 + # Create FITS table columns
  1530 + col_vp = gammalib.GFitsTableStringCol('VP', nrows, 9)
  1531 + col_object = gammalib.GFitsTableStringCol('OBJECT', nrows, 20)
  1532 + col_obs_id = gammalib.GFitsTableStringCol('OBS_ID', nrows, 9)
  1533 + col_dsd_id = gammalib.GFitsTableStringCol('DSD_ID', nrows, 15)
  1534 + col_glon_scz = gammalib.GFitsTableFloatCol('GLON_SCZ', nrows)
  1535 + col_glat_scz = gammalib.GFitsTableFloatCol('GLAT_SCZ', nrows)
  1536 + col_date_start = gammalib.GFitsTableStringCol('DATE_START', nrows, 20)
  1537 + col_date_stop = gammalib.GFitsTableStringCol('DATE_STOP', nrows, 20)
  1538 + col_tjd_start = gammalib.GFitsTableShortCol('TJD_START', nrows)
  1539 + col_tics_start = gammalib.GFitsTableLongCol('TICS_START', nrows)
  1540 + col_tjd_stop = gammalib.GFitsTableShortCol('TJD_STOP', nrows)
  1541 + col_tics_stop = gammalib.GFitsTableLongCol('TICS_STOP', nrows)
  1542 + col_visday = gammalib.GFitsTableShortCol('VISDAY', nrows)
  1543 + col_vistim = gammalib.GFitsTableLongCol('VISTIM', nrows)
  1544 + col_vieday = gammalib.GFitsTableShortCol('VIEDAY', nrows)
  1545 + col_vietim = gammalib.GFitsTableLongCol('VIETIM', nrows)
  1546 + col_verno = gammalib.GFitsTableShortCol('VERNO', nrows)
  1547 + col_quality = gammalib.GFitsTableShortCol('QUALITY', nrows)
  1548 + col_evp = gammalib.GFitsTableStringCol('EVP', nrows, 60)
  1549 + col_nevents = gammalib.GFitsTableLongCol('NEVENTS', nrows)
  1550 +
  1551 + # Set column units
  1552 + col_glon_scz.unit('deg')
  1553 + col_glat_scz.unit('deg')
  1554 + col_tjd_start.unit('days')
  1555 + col_tics_start.unit('0.000125 s')
  1556 + col_tjd_stop.unit('days')
  1557 + col_tics_stop.unit('0.000125 s')
  1558 + col_visday.unit('days')
  1559 + col_vistim.unit('0.000125 s')
  1560 + col_vieday.unit('days')
  1561 + col_vietim.unit('0.000125 s')
  1562 +
  1563 + # Fill columns
  1564 + for i, evp_record in enumerate(evp_list):
  1565 + col_vp[i] = evp_record['vp']
  1566 + col_object[i] = evp_record['object']
  1567 + col_obs_id[i] = '%.1f' % evp_record['obs_id']
  1568 + col_dsd_id[i] = evp_record['dsd_id']
  1569 + col_glon_scz[i] = evp_record['glon_scz']
  1570 + col_glat_scz[i] = evp_record['glat_scz']
  1571 + col_date_start[i] = evp_record['date_start']
  1572 + col_date_stop[i] = evp_record['date_stop']
  1573 + col_tjd_start[i] = evp_record['tjd_start']
  1574 + col_tics_start[i] = evp_record['tics_start']
  1575 + col_tjd_stop[i] = evp_record['tjd_stop']
  1576 + col_tics_stop[i] = evp_record['tics_stop']
  1577 + col_visday[i] = evp_record['visday']
  1578 + col_vistim[i] = evp_record['vistim']
  1579 + col_vieday[i] = evp_record['vieday']
  1580 + col_vietim[i] = evp_record['vietim']
  1581 + col_verno[i] = evp_record['verno']
  1582 + col_quality[i] = evp_record['quality']
  1583 + col_evp[i] = evp_record['filename']
  1584 + col_nevents[i] = evp_record['nevents']
  1585 +
  1586 + # Create FITS table
  1587 + table = gammalib.GFitsBinTable()
  1588 +
  1589 + # Append columns to table
  1590 + table.append(col_vp)
  1591 + table.append(col_object)
  1592 + table.append(col_obs_id)
  1593 + table.append(col_dsd_id)
  1594 + table.append(col_glon_scz)
  1595 + table.append(col_glat_scz)
  1596 + table.append(col_date_start)
  1597 + table.append(col_date_stop)
  1598 + table.append(col_tjd_start)
  1599 + table.append(col_tics_start)
  1600 + table.append(col_tjd_stop)
  1601 + table.append(col_tics_stop)
  1602 + table.append(col_visday)
  1603 + table.append(col_vistim)
  1604 + table.append(col_vieday)
  1605 + table.append(col_vietim)
  1606 + table.append(col_verno)
  1607 + table.append(col_quality)
  1608 + table.append(col_evp)
  1609 + table.append(col_nevents)
  1610 +
  1611 + # Set extension name
  1612 + table.extname('EVP')
  1613 +
  1614 + # Return table
  1615 + return table
  1616 +
  1617 + def _build_oad_table(self, oad_list):
  1618 + """
  1619 + Build OAD FITS table
  1620 +
  1621 + Parameters
  1622 + ----------
  1623 + oad_list : list of dict
  1624 + OAD database
  1625 +
  1626 + Returns
  1627 + -------
  1628 + table : `~gammalib.GFitsBinTable`
  1629 + OAD FITS table
  1630 + """
  1631 + # Determine number of rows in table
  1632 + nrows = len(oad_list)
  1633 +
  1634 + # Create FITS table columns
  1635 + col_vp = gammalib.GFitsTableStringCol('VP', nrows, 9)
  1636 + col_object = gammalib.GFitsTableStringCol('OBJECT', nrows, 20)
  1637 + col_obs_id = gammalib.GFitsTableStringCol('OBS_ID', nrows, 9)
  1638 + col_dsd_id = gammalib.GFitsTableStringCol('DSD_ID', nrows, 15)
  1639 + col_glon_scz = gammalib.GFitsTableFloatCol('GLON_SCZ', nrows)
  1640 + col_glat_scz = gammalib.GFitsTableFloatCol('GLAT_SCZ', nrows)
  1641 + col_date_start = gammalib.GFitsTableStringCol('DATE_START', nrows, 20)
  1642 + col_date_stop = gammalib.GFitsTableStringCol('DATE_STOP', nrows, 20)
  1643 + col_tjd_start = gammalib.GFitsTableShortCol('TJD_START', nrows)
  1644 + col_tics_start = gammalib.GFitsTableLongCol('TICS_START', nrows)
  1645 + col_tjd_stop = gammalib.GFitsTableShortCol('TJD_STOP', nrows)
  1646 + col_tics_stop = gammalib.GFitsTableLongCol('TICS_STOP', nrows)
  1647 + col_quality = gammalib.GFitsTableShortCol('QUALITY', nrows)
  1648 + col_oad = gammalib.GFitsTableStringCol('OAD', nrows, 60)
  1649 + col_superpackets = gammalib.GFitsTableLongCol('SUPERPACKETS', nrows)
  1650 +
  1651 + # Set column units
  1652 + col_glon_scz.unit('deg')
  1653 + col_glat_scz.unit('deg')
  1654 + col_tjd_start.unit('days')
  1655 + col_tics_start.unit('0.000125 s')
  1656 + col_tjd_stop.unit('days')
  1657 + col_tics_stop.unit('0.000125 s')
  1658 +
  1659 + # Fill columns
  1660 + for i, oad_record in enumerate(oad_list):
  1661 + col_vp[i] = oad_record['vp']
  1662 + col_object[i] = oad_record['object']
  1663 + col_obs_id[i] = '%.1f' % oad_record['obs_id']
  1664 + col_dsd_id[i] = oad_record['dsd_id']
  1665 + col_glon_scz[i] = oad_record['glon_scz']
  1666 + col_glat_scz[i] = oad_record['glat_scz']
  1667 + col_date_start[i] = oad_record['date_start']
  1668 + col_date_stop[i] = oad_record['date_stop']
  1669 + col_tjd_start[i] = oad_record['tjd_start']
  1670 + col_tics_start[i] = oad_record['tics_start']
  1671 + col_tjd_stop[i] = oad_record['tjd_stop']
  1672 + col_tics_stop[i] = oad_record['tics_stop']
  1673 + col_quality[i] = oad_record['quality']
  1674 + col_oad[i] = oad_record['filename']
  1675 + col_superpackets[i] = oad_record['superpackets']
  1676 +
  1677 + # Create FITS table
  1678 + table = gammalib.GFitsBinTable()
  1679 +
  1680 + # Append columns to table
  1681 + table.append(col_vp)
  1682 + table.append(col_object)
  1683 + table.append(col_obs_id)
  1684 + table.append(col_dsd_id)
  1685 + table.append(col_glon_scz)
  1686 + table.append(col_glat_scz)
  1687 + table.append(col_date_start)
  1688 + table.append(col_date_stop)
  1689 + table.append(col_tjd_start)
  1690 + table.append(col_tics_start)
  1691 + table.append(col_tjd_stop)
  1692 + table.append(col_tics_stop)
  1693 + table.append(col_quality)
  1694 + table.append(col_oad)
  1695 + table.append(col_superpackets)
  1696 +
  1697 + # Set extension name
  1698 + table.extname('OAD')
  1699 +
  1700 + # Return table
  1701 + return table
  1702 +
  1703 + def _build_tim_table(self, tim_list):
  1704 + """
  1705 + Build TIM FITS table
  1706 +
  1707 + Parameters
  1708 + ----------
  1709 + tim_list : list of dict
  1710 + TIM database
  1711 +
  1712 + Returns
  1713 + -------
  1714 + table : `~gammalib.GFitsBinTable`
  1715 + TIM FITS table
  1716 + """
  1717 + # Determine number of rows in table
  1718 + nrows = len(tim_list)
  1719 +
  1720 + # Create FITS table columns
  1721 + col_vp = gammalib.GFitsTableStringCol('VP', nrows, 9)
  1722 + col_object = gammalib.GFitsTableStringCol('OBJECT', nrows, 20)
  1723 + col_obs_id = gammalib.GFitsTableStringCol('OBS_ID', nrows, 9)
  1724 + col_dsd_id = gammalib.GFitsTableStringCol('DSD_ID', nrows, 15)
  1725 + col_glon_scz = gammalib.GFitsTableFloatCol('GLON_SCZ', nrows)
  1726 + col_glat_scz = gammalib.GFitsTableFloatCol('GLAT_SCZ', nrows)
  1727 + col_date_start = gammalib.GFitsTableStringCol('DATE_START', nrows, 20)
  1728 + col_date_stop = gammalib.GFitsTableStringCol('DATE_STOP', nrows, 20)
  1729 + col_tjd_start = gammalib.GFitsTableShortCol('TJD_START', nrows)
  1730 + col_tics_start = gammalib.GFitsTableLongCol('TICS_START', nrows)
  1731 + col_tjd_stop = gammalib.GFitsTableShortCol('TJD_STOP', nrows)
  1732 + col_tics_stop = gammalib.GFitsTableLongCol('TICS_STOP', nrows)
  1733 + col_telapse = gammalib.GFitsTableDoubleCol('TELAPSE', nrows)
  1734 + col_ontime = gammalib.GFitsTableDoubleCol('ONTIME', nrows)
  1735 + col_quality = gammalib.GFitsTableShortCol('QUALITY', nrows)
  1736 + col_tim = gammalib.GFitsTableStringCol('TIM', nrows, 60)
  1737 + col_intervals = gammalib.GFitsTableLongCol('INTERVALS', nrows)
  1738 +
  1739 + # Set column units
  1740 + col_glon_scz.unit('deg')
  1741 + col_glat_scz.unit('deg')
  1742 + col_tjd_start.unit('days')
  1743 + col_tics_start.unit('0.000125 s')
  1744 + col_tjd_stop.unit('days')
  1745 + col_tics_stop.unit('0.000125 s')
  1746 + col_telapse.unit('s')
  1747 + col_ontime.unit('s')
  1748 +
  1749 + # Fill columns
  1750 + for i, tim_record in enumerate(tim_list):
  1751 + col_vp[i] = tim_record['vp']
  1752 + col_object[i] = tim_record['object']
  1753 + col_obs_id[i] = '%.1f' % tim_record['obs_id']
  1754 + col_dsd_id[i] = tim_record['dsd_id']
  1755 + col_glon_scz[i] = tim_record['glon_scz']
  1756 + col_glat_scz[i] = tim_record['glat_scz']
  1757 + col_date_start[i] = tim_record['date_start']
  1758 + col_date_stop[i] = tim_record['date_stop']
  1759 + col_tjd_start[i] = tim_record['tjd_start']
  1760 + col_tics_start[i] = tim_record['tics_start']
  1761 + col_tjd_stop[i] = tim_record['tjd_stop']
  1762 + col_tics_stop[i] = tim_record['tics_stop']
  1763 + col_telapse[i] = tim_record['telapse']
  1764 + col_ontime[i] = tim_record['ontime']
  1765 + col_quality[i] = tim_record['quality']
  1766 + col_tim[i] = tim_record['filename']
  1767 + col_intervals[i] = tim_record['intervals']
  1768 +
  1769 + # Create FITS table
  1770 + table = gammalib.GFitsBinTable()
  1771 +
  1772 + # Append columns to table
  1773 + table.append(col_vp)
  1774 + table.append(col_object)
  1775 + table.append(col_obs_id)
  1776 + table.append(col_dsd_id)
  1777 + table.append(col_glon_scz)
  1778 + table.append(col_glat_scz)
  1779 + table.append(col_date_start)
  1780 + table.append(col_date_stop)
  1781 + table.append(col_tjd_start)
  1782 + table.append(col_tics_start)
  1783 + table.append(col_tjd_stop)
  1784 + table.append(col_tics_stop)
  1785 + table.append(col_telapse)
  1786 + table.append(col_ontime)
  1787 + table.append(col_quality)
  1788 + table.append(col_tim)
  1789 + table.append(col_intervals)
  1790 +
  1791 + # Set extension name
  1792 + table.extname('TIM')
  1793 +
  1794 + # Return table
  1795 + return table
  1796 +
  1797 + def _build_xml_table(self, xml_list):
  1798 + """
  1799 + Build XML FITS table
  1800 +
  1801 + Parameters
  1802 + ----------
  1803 + xml_list : list of dict
  1804 + XML database
  1805 +
  1806 + Returns
  1807 + -------
  1808 + table : `~gammalib.GFitsBinTable`
  1809 + XML FITS table
  1810 + """
  1811 + # Determine number of rows in table
  1812 + nrows = len(xml_list)
  1813 +
  1814 + # Create FITS table columns
  1815 + col_vp = gammalib.GFitsTableStringCol('VP', nrows, 9)
  1816 + col_object = gammalib.GFitsTableStringCol('OBJECT', nrows, 20)
  1817 + col_obs_id = gammalib.GFitsTableStringCol('OBS_ID', nrows, 9)
  1818 + col_glon_scz = gammalib.GFitsTableFloatCol('GLON_SCZ', nrows)
  1819 + col_glat_scz = gammalib.GFitsTableFloatCol('GLAT_SCZ', nrows)
  1820 + col_date_start = gammalib.GFitsTableStringCol('DATE_START', nrows, 20)
  1821 + col_date_stop = gammalib.GFitsTableStringCol('DATE_STOP', nrows, 20)
  1822 + col_tjd_start = gammalib.GFitsTableShortCol('TJD_START', nrows)
  1823 + col_tics_start = gammalib.GFitsTableLongCol('TICS_START', nrows)
  1824 + col_tjd_stop = gammalib.GFitsTableShortCol('TJD_STOP', nrows)
  1825 + col_tics_stop = gammalib.GFitsTableLongCol('TICS_STOP', nrows)
  1826 + col_visday = gammalib.GFitsTableShortCol('VISDAY', nrows)
  1827 + col_vistim = gammalib.GFitsTableLongCol('VISTIM', nrows)
  1828 + col_vieday = gammalib.GFitsTableShortCol('VIEDAY', nrows)
  1829 + col_vietim = gammalib.GFitsTableLongCol('VIETIM', nrows)
  1830 + col_telapse = gammalib.GFitsTableDoubleCol('TELAPSE', nrows)
  1831 + col_ontime = gammalib.GFitsTableDoubleCol('ONTIME', nrows)
  1832 + col_evp_verno = gammalib.GFitsTableShortCol('EVP_VERNO', nrows)
  1833 + col_evp_quality = gammalib.GFitsTableShortCol('EVP_QUALITY', nrows)
  1834 + col_tim_quality = gammalib.GFitsTableShortCol('TIM_QUALITY', nrows)
  1835 + col_oad_quality = gammalib.GFitsTableShortCol('OAD_QUALITY', nrows, 43)
  1836 + col_xml = gammalib.GFitsTableStringCol('XML', nrows, 60)
  1837 + col_evp = gammalib.GFitsTableStringCol('EVP', nrows, 60)
  1838 + col_tim = gammalib.GFitsTableStringCol('TIM', nrows, 60)
  1839 + col_oad = gammalib.GFitsTableStringCol('OAD', nrows, 60, 43)
  1840 + col_nevents = gammalib.GFitsTableLongCol('NEVENTS', nrows)
  1841 + col_model = gammalib.GFitsTableStringCol('MODE', nrows, 8)
  1842 +
  1843 + # Set column units
  1844 + col_glon_scz.unit('deg')
  1845 + col_glat_scz.unit('deg')
  1846 + col_tjd_start.unit('days')
  1847 + col_tics_start.unit('0.000125 s')
  1848 + col_tjd_stop.unit('days')
  1849 + col_tics_stop.unit('0.000125 s')
  1850 + col_visday.unit('days')
  1851 + col_vistim.unit('0.000125 s')
  1852 + col_vieday.unit('days')
  1853 + col_vietim.unit('0.000125 s')
  1854 + col_telapse.unit('s')
  1855 + col_ontime.unit('s')
  1856 +
  1857 + # Fill columns
  1858 + for i, xml_record in enumerate(xml_list):
  1859 + col_vp[i] = xml_record['vp']
  1860 + col_object[i] = xml_record['object']
  1861 + col_obs_id[i] = '%.1f' % xml_record['obs_id']
  1862 + col_glon_scz[i] = xml_record['glon_scz']
  1863 + col_glat_scz[i] = xml_record['glat_scz']
  1864 + col_date_start[i] = xml_record['date_start']
  1865 + col_date_stop[i] = xml_record['date_stop']
  1866 + col_tjd_start[i] = xml_record['tjd_start']
  1867 + col_tics_start[i] = xml_record['tics_start']
  1868 + col_tjd_stop[i] = xml_record['tjd_stop']
  1869 + col_tics_stop[i] = xml_record['tics_stop']
  1870 + col_visday[i] = xml_record['visday']
  1871 + col_vistim[i] = xml_record['vistim']
  1872 + col_vieday[i] = xml_record['vieday']
  1873 + col_vietim[i] = xml_record['vietim']
  1874 + col_telapse[i] = xml_record['telapse']
  1875 + col_ontime[i] = xml_record['ontime']
  1876 + col_evp_verno[i] = xml_record['evp_verno']
  1877 + col_evp_quality[i] = xml_record['evp_quality']
  1878 + col_tim_quality[i] = xml_record['tim_quality']
  1879 + col_xml[i] = xml_record['xmlname']
  1880 + col_evp[i] = xml_record['evpname']
  1881 + col_tim[i] = xml_record['timname']
  1882 + col_nevents[i] = xml_record['nevents']
  1883 + for k, oad_quality in enumerate(xml_record['oad_quality']):
  1884 + col_oad_quality[i,k] = oad_quality
  1885 + for k, oadname in enumerate(xml_record['oadnames']):
  1886 + col_oad[i,k] = oadname
  1887 + mode = 'UNKNOWN'
  1888 + if xml_record['mode'] == 'STD':
  1889 + mode = 'STANDARD'
  1890 + elif xml_record['mode'] == 'S80':
  1891 + mode = 'SOLAR80'
  1892 + col_model[i] = mode
  1893 +
  1894 + # Create FITS table
  1895 + table = gammalib.GFitsBinTable()
  1896 +
  1897 + # Append columns to table
  1898 + table.append(col_vp)
  1899 + table.append(col_object)
  1900 + table.append(col_obs_id)
  1901 + table.append(col_glon_scz)
  1902 + table.append(col_glat_scz)
  1903 + table.append(col_date_start)
  1904 + table.append(col_date_stop)
  1905 + table.append(col_tjd_start)
  1906 + table.append(col_tics_start)
  1907 + table.append(col_tjd_stop)
  1908 + table.append(col_tics_stop)
  1909 + table.append(col_visday)
  1910 + table.append(col_vistim)
  1911 + table.append(col_vieday)
  1912 + table.append(col_vietim)
  1913 + table.append(col_telapse)
  1914 + table.append(col_ontime)
  1915 + table.append(col_nevents)
  1916 + table.append(col_evp_verno)
  1917 + table.append(col_evp_quality)
  1918 + table.append(col_tim_quality)
  1919 + table.append(col_oad_quality)
  1920 + table.append(col_xml)
  1921 + table.append(col_evp)
  1922 + table.append(col_tim)
  1923 + table.append(col_oad)
  1924 +
  1925 + # Set extension name
  1926 + table.extname('XML')
  1927 +
  1928 + # Return table
  1929 + return table
  1930 +
  1931 +
  1932 + # Public methods
  1933 + def process(self):
  1934 + """
  1935 + Process the script
  1936 + """
  1937 + # Get parameters
  1938 + self._get_parameters()
  1939 +
  1940 + # Optionally download COMPTEL HEASARC archive
  1941 + if self['download'].boolean():
  1942 + self._download()
  1943 +
  1944 + # Log header
  1945 + self._log_header1(gammalib.TERSE, 'Get list of EVP files')
  1946 +
  1947 + # Get EVP list
  1948 + self._evp_list = self._get_evp_list()
  1949 +
  1950 + # Log header
  1951 + self._log_header1(gammalib.TERSE, 'Get list of OAD files')
  1952 +
  1953 + # Get OAD list
  1954 + self._oad_list = self._get_oad_list()
  1955 +
  1956 + # Log header
  1957 + self._log_header1(gammalib.TERSE, 'Get list of TIM files')
  1958 +
  1959 + # Get TIM list
  1960 + self._tim_list = self._get_tim_list()
  1961 +
  1962 + # Log header
  1963 + self._log_header1(gammalib.TERSE, 'Generate list of XML files')
  1964 +
  1965 + # Get XML list
  1966 + xml_list = self._get_xml_list(self._evp_list, self._tim_list, self._oad_list)
  1967 +
  1968 + # Log header
  1969 + self._log_header1(gammalib.TERSE, 'Compare with COMPASS database')
  1970 +
  1971 + # Check XML list
  1972 + self._xml_list = self._check_xml_list(xml_list)
  1973 +
  1974 + # Return
  1975 + return
  1976 +
  1977 + def save(self):
  1978 + """
  1979 + Save database file
  1980 + """
  1981 + # Write header
  1982 + self._log_header1(gammalib.TERSE, 'Save database')
  1983 +
  1984 + # Create database directory structure
  1985 + self._create_dbase_dirs()
  1986 +
  1987 + # Generate TIM and XML files
  1988 + for xml in self._xml_list:
  1989 + self._build_tim_file(xml)
  1990 + self._build_xml_file(xml)
  1991 +
  1992 + # Check database
  1993 + self._check_database(self._xml_list)
  1994 +
  1995 + # Build database file
  1996 + self._build_database_file(self._evp_list, self._oad_list,
  1997 + self._tim_list, self._xml_list)
  1998 +
  1999 + # Return
  2000 + return
  2001 +
  2002 +
  2003 +# ======================== #
  2004 +# Main routine entry point #
  2005 +# ======================== #
  2006 +if __name__ == '__main__':
  2007 +
  2008 + # Create instance of application
  2009 + app = comgendb(sys.argv)
  2010 +
  2011 + # Execute application
  2012 + app.execute()
... ...
test/Makefile.am
  View file @838aa04
1 1 #############################################################################
2 2 # Makefile.am for ctools unit tests #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2021 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2022 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -71,6 +71,7 @@ clean-local:
71 71 rm -rf $(top_builddir)/test/comobsback_*
72 72 rm -rf $(top_builddir)/test/comobsbin_*
73 73 rm -rf $(top_builddir)/test/comobssim_*
  74 + rm -rf $(top_builddir)/test/comgendb_*
74 75 rm -rf $(top_builddir)/test/.DS_Store
75 76  
76 77 # Tests to be done on installed version
... ...
test/test_comgendb.py 0 → 100755
  View file @838aa04
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the comgendb script.
  4 +#
  5 +# Copyright (C) 2022 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 comscripts
  24 +from testing import test
  25 +
  26 +
  27 +# ============================== #
  28 +# Test class for comgendb script #
  29 +# ============================== #
  30 +class Test(test):
  31 + """
  32 + Test class for comgendb script
  33 +
  34 + This test class makes unit tests for the comgendb script by using it
  35 + from the command line and from Python.
  36 + """
  37 +
  38 + # Constructor
  39 + def __init__(self):
  40 + """
  41 + Constructor
  42 + """
  43 + # Call base class constructor
  44 + test.__init__(self)
  45 +
  46 + # Set test datasets and parameters
  47 + self._archive = self._datadir + '/comptel/data'
  48 +
  49 + # Return
  50 + return
  51 +
  52 + # Set test functions
  53 + def set(self):
  54 + """
  55 + Set all test functions.
  56 + """
  57 + # Set test name
  58 + self.name('comgendb')
  59 +
  60 + # Append tests
  61 + self.append(self._test_cmd, 'Test comgendb on command line')
  62 + self.append(self._test_python, 'Test comgendb from Python')
  63 +
  64 + # Return
  65 + return
  66 +
  67 + # Test command line use
  68 + def _test_cmd(self):
  69 + """
  70 + Test command line use.
  71 + """
  72 + # Set script name
  73 + comgendb = self._comscript('comgendb')
  74 +
  75 + # Setup command
  76 + cmd = comgendb+' archive="'+self._archive+'" dbase="comgendb_cmd1" '+ \
  77 + 'logfile="comgendb_cmd1.log" chatter=1'
  78 +
  79 + # Check if execution was successful
  80 + self.test_assert(self._execute(cmd) == 0,
  81 + 'Check successful execution from command line')
  82 +
  83 + # Check result
  84 + self._check_result('comgendb_cmd1')
  85 +
  86 + # Setup command
  87 + cmd = comgendb+' archive="folder_that_does_not_exist" dbase="comgendb_cmd2" '+ \
  88 + 'logfile="comgendb_cmd2.log" chatter=1'
  89 +
  90 + # Check if execution was successful
  91 + self.test_assert(self._execute(cmd, success=False) != 0,
  92 + 'Check invalid input file when executed from command line')
  93 +
  94 + # Check --help option
  95 + self._check_help(comgendb)
  96 +
  97 + # Return
  98 + return
  99 +
  100 + # Test from Python
  101 + def _test_python(self):
  102 + """
  103 + Test from Python
  104 + """
  105 + # Set-up script
  106 + db = comscripts.comgendb()
  107 + db['archive'] = self._archive
  108 + db['dbase'] = 'comgendb_py1'
  109 + db['logfile'] = 'comgendb_py1.log'
  110 + db['chatter'] = 2
  111 +
  112 + # Run script and save result
  113 + db.logFileOpen() # Make sure we get a log file
  114 + db.run()
  115 + db.save()
  116 +
  117 + # Check result
  118 + self._check_result('comgendb_py1')
  119 +
  120 + # Return
  121 + return
  122 +
  123 + # Check result
  124 + def _check_result(self, dbase):
  125 + """
  126 + Check result
  127 + """
  128 + # Set database filename
  129 + filename = '%s/dbase.fits' % (dbase)
  130 +
  131 + # Load database file
  132 + fits = gammalib.GFits(filename)
  133 +
  134 + # Check existence of headers
  135 + self.test_assert(fits.contains('EVP'), 'Check for "EVP" extension')
  136 + self.test_assert(fits.contains('OAD'), 'Check for "OAD" extension')
  137 + self.test_assert(fits.contains('TIM'), 'Check for "TIM" extension')
  138 + self.test_assert(fits.contains('XML'), 'Check for "XML" extension')
  139 +
  140 + # Check number of rows in extension
  141 + self.test_value(fits.table('EVP').nrows(), 1, 'Number of rows in "EVP" table')
  142 + self.test_value(fits.table('OAD').nrows(), 1, 'Number of rows in "OAD" table')
  143 + self.test_value(fits.table('TIM').nrows(), 0, 'Number of rows in "TIM" table')
  144 + self.test_value(fits.table('XML').nrows(), 0, 'Number of rows in "XML" table')
  145 +
  146 + # Return
  147 + return
... ...
test/test_python_comscripts.py
  View file @838aa04
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # This scripts performs unit tests for comscripts
4 4 #
5   -# Copyright (C) 2021 Juergen Knoedlseder
  5 +# Copyright (C) 2021-2022 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -23,6 +23,7 @@ import sys
23 23 import gammalib
24 24 import comscripts
25 25 sys.path.append(os.path.dirname(os.path.abspath(__file__)))
  26 +import test_comgendb
26 27 import test_comobsselect
27 28 import test_comobsbin
28 29 import test_comobsback
... ... @@ -99,7 +100,8 @@ def test(installed=False, debug=False):
99 100 os.system('chmod u+w pfiles/*')
100 101  
101 102 # Define list of test suites
102   - tests = [test_comobsselect.Test(),
  103 + tests = [test_comgendb.Test(),
  104 + test_comobsselect.Test(),
103 105 test_comobsbin.Test(),
104 106 test_comobsback.Test(),
105 107 test_comobsmodel.Test(),
... ...