gammalib / gammalib

Commit 708e00f9262712cf2acacf5a7373cc3e9ee41f07

Authored by Jürgen Knödlseder
1 parent a4cacc7c

Add script to test BGDLIXF background model

Inline Side-by-side
Showing 1 changed file with 169 additions and 0 deletions
inst/com/test/dev/test_bgdlixf.py 0 → 100755
  View file @708e00f
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Test BGDLIXF rate fitting
  4 +# --------------------------------------------------------------------------
  5 +#
  6 +# Copyright (C) 2023 Juergen Knoedlseder
  7 +#
  8 +# This program is free software: you can redistribute it and/or modify
  9 +# it under the terms of the GNU General Public License as published by
  10 +# the Free Software Foundation, either version 3 of the License, or
  11 +# (at your option) any later version.
  12 +#
  13 +# This program is distributed in the hope that it will be useful,
  14 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +# GNU General Public License for more details.
  17 +#
  18 +# You should have received a copy of the GNU General Public License
  19 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +#
  21 +# ==========================================================================
  22 +import gammalib
  23 +import glob
  24 +import math
  25 +import matplotlib.pyplot as plt
  26 +
  27 +
  28 +# =================== #
  29 +# Get Veto dome rates #
  30 +# =================== #
  31 +def get_veto_rates(xmlname):
  32 + """
  33 + Get Veto dome rates
  34 +
  35 + Parameters
  36 + ----------
  37 + xmlname : str
  38 + Observation definition XML file name
  39 +
  40 + Returns
  41 + -------
  42 + rates : dict
  43 + Veto rates
  44 + """
  45 + # Initialise rates dictionary
  46 + rates = {'times': [], 'rates': []}
  47 +
  48 + # Load observations
  49 + obs = gammalib.GObservations(xmlname)
  50 +
  51 + # Loop over observations
  52 + for com in obs:
  53 +
  54 + # Get identifier
  55 + rates['id'] = com.id()
  56 +
  57 + # Get 2nd Veto dome rates housekeeping data
  58 + hkd = com.hkds()['SCV2M']
  59 + print(hkd)
  60 +
  61 + # Extract data
  62 + for i in range(hkd.size()):
  63 + time = hkd.time(i).copy()
  64 + rate = hkd.value(i) * 12.8 / 2.048
  65 + if rate > 0.0:
  66 + rates['times'].append(time)
  67 + rates['rates'].append(rate)
  68 +
  69 + # Return rates
  70 + return rates
  71 +
  72 +
  73 +# ====================================== #
  74 +# Check time ordering of Veto dome rates #
  75 +# ====================================== #
  76 +def check_veto_rates(xmlname):
  77 + """
  78 + Check time ordering of Veto dome rates
  79 +
  80 + Parameters
  81 + ----------
  82 + xmlname : str
  83 + Observation definition XML file name
  84 + """
  85 + # Load observations
  86 + obs = gammalib.GObservations(xmlname)
  87 +
  88 + # Loop over observations
  89 + for com in obs:
  90 +
  91 + # Get 2nd Veto dome rates housekeeping data
  92 + hkd = com.hkds()['SCV2M']
  93 +
  94 + # Check times
  95 + nbad = 0
  96 + for i in range(hkd.size()):
  97 + if i == 0:
  98 + last_time = hkd.time(i)
  99 + else:
  100 + if hkd.time(i) <= last_time:
  101 + nbad += 1
  102 + print('%s: %.5f <= %.5f' % (com.id(), hkd.time(i).mjd(), last_time.mjd()))
  103 + if nbad == 0:
  104 + print('%s: ok.' % (com.id()))
  105 +
  106 + # Return
  107 + return
  108 +
  109 +
  110 +# ==================== #
  111 +# Show Veto dome rates #
  112 +# ==================== #
  113 +def show_veto_rates(rates):
  114 + """
  115 + Show Veto dome rates
  116 + """
  117 + # Create figure
  118 + fig = plt.figure(figsize=(14,4))
  119 + fig.subplots_adjust(left=0.07, right=0.99, top=0.93, bottom=0.12)
  120 + plt.title('Veto dome rates for %s' % rates['id'])
  121 +
  122 + # Setup vectors
  123 + x = [time.mjd() for time in rates['times']]
  124 + y = [rate for rate in rates['rates']]
  125 +
  126 + # Plot vectors
  127 + plt.plot(x, y, 'r-')
  128 +
  129 + # Set attributes
  130 + plt.xlabel(r'MJD (days)')
  131 + plt.ylabel(r'Rate (triggers/sec)')
  132 +
  133 + # Show plot
  134 + plt.show()
  135 +
  136 + # Return
  137 + return
  138 +
  139 +
  140 +# ======================== #
  141 +# Main routine entry point #
  142 +# ======================== #
  143 +if __name__ == '__main__':
  144 + # Dump header
  145 + print('')
  146 + print('*****************************')
  147 + print('* Test BGDLIXF rate fitting *')
  148 + print('*****************************')
  149 +
  150 + # Set file names
  151 + #hkdname = '/project-data/comptel/data/phase01/vp0001_0/m20074_hkd.fits'
  152 + #hkdname = '/project-data/comptel/data/phase07/vp0701_0/m26689_hkd.fits'
  153 + #xmlname = 'dbase/xml/vp0001_0.xml'
  154 +
  155 + # Get all XML files
  156 + xmlnames = glob.glob('dbase/xml/vp*.xml')
  157 + xmlnames.sort()
  158 +
  159 + # Loop over XML files
  160 + for xmlname in xmlnames:
  161 +
  162 + # Check Veto rates
  163 + check_veto_rates(xmlname)
  164 +
  165 + # Get Veto dome rates
  166 + #rates = get_veto_rates(xmlname)
  167 +
  168 + # Show Veto dome rates
  169 + #show_veto_rates(rates)
... ...