Luigi Tibaldo / ctools

Commit c057d37145ca3e06447980fa6c063ffebb34344f

Authored by Jürgen Knödlseder
1 parent cfd3f245

Correctly extract subarray response functions and update documentation

Inline Side-by-side
Showing 3 changed files with 66 additions and 432 deletions
cscripts/csadd2caldb.py deleted
View file @cfd3f24
1   -#! /usr/bin/env python
2   -# ==========================================================================
3   -# Add IRFs to CALDB
4   -#
5   -# Copyright (C) 2021 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 glob
24   -import tarfile
25   -import gammalib
26   -import ctools
27   -
28   -
29   -# ================= #
30   -# csadd2caldb class #
31   -# ================= #
32   -class csadd2caldb(ctools.cscript):
33   - """
34   - Add IRFs to CALDB
35   - """
36   -
37   - # Constructor
38   - def __init__(self, *argv):
39   - """
40   - Constructor
41   - """
42   - # Make sure that parfile exists
43   - self._parfile()
44   -
45   - # Initialise application by calling the base class constructor
46   - self._init_cscript(self.__class__.__name__, ctools.__version__, argv)
47   -
48   - # Return
49   - return
50   -
51   -
52   - # Private methods
53   - def _parfile(self):
54   - """
55   - Check if parfile exists. If parfile does not exist then create a
56   - default parfile. This kluge avoids shipping the cscript with a parfile.
57   - """
58   - # Set parfile name
59   - parfile = self.__class__.__name__+'.par'
60   -
61   - # Load parfile or create
62   - try:
63   - pars = gammalib.GApplicationPars(parfile)
64   - except:
65   - # Signal that parfile was not found
66   - print('Parfile '+parfile+' not found. Create default parfile.')
67   -
68   - # Create default parfile
69   - pars = gammalib.GApplicationPars()
70   - pars.append(gammalib.GApplicationPar('indir','f','a','cta-prod5-zenodo-fitsonly-v0','','','Input IRF folder'))
71   - pars.append(gammalib.GApplicationPar('outdir','f','a','$CALDB','','','Output caldb folder'))
72   - pars.append_standard()
73   - pars.save(parfile)
74   -
75   - # Return
76   - return
77   -
78   - def _get_parameters(self):
79   - """
80   - Get parameters from parfile
81   - """
82   - # Query parameters
83   - self['indir'].filename()
84   - self['outdir'].filename()
85   -
86   - # Write input parameters into logger
87   - self._log_parameters(gammalib.TERSE)
88   -
89   - # Return
90   - return
91   -
92   - def _collect_tarfiles(self, indir):
93   - """
94   - Collect tarfiles
95   -
96   - Parameters
97   - ----------
98   - indir : string
99   - Input directory
100   -
101   - Returns
102   - -------
103   - tarfiles : list of str
104   - List of tarfiles
105   - """
106   - # Initialise list of tarfiles
107   - tarfiles = []
108   -
109   - # Search for tarfiles in current directory
110   - files = glob.glob('%s/*.FITS.tar.gz' % indir.url())
111   - tarfiles.extend(files)
112   -
113   - # Search for tarfiles in 'files' subdirectory
114   - files = glob.glob('%s/fits/*.FITS.tar.gz' % indir.url())
115   - tarfiles.extend(files)
116   -
117   - # Return tarfiles
118   - return tarfiles
119   -
120   - def _extract_irfs(self, file, outdir):
121   - """
122   - Collect IRFs from tarfile
123   -
124   - Parameters
125   - ----------
126   - file : str
127   - Tarfile with IRFs
128   - outdir : str
129   - Output directory
130   -
131   - Returns
132   - -------
133   - irfs : list
134   - List of IRF dictionaries
135   - """
136   - # Initialise list of extracted IRFs
137   - irfs = []
138   -
139   - # Open tarfile
140   - tar = tarfile.open(file)
141   -
142   - # Loop over all members in tarfile
143   - for member in tar.getmembers():
144   -
145   - # Extract Instrument
146   - if 'South' in member.name:
147   - site = 'South'
148   - else:
149   - site = 'North'
150   - if '180000s' in member.name:
151   - duration = '50h'
152   - elif '18000s' in member.name:
153   - duration = '5h'
154   - elif '1800s' in member.name:
155   - duration = '0.5h'
156   - if 'SouthAz' in member.name:
157   - orientation = '_S'
158   - azimuth = 180.0
159   - elif 'NorthAz' in member.name:
160   - orientation = '_N'
161   - azimuth = 0.0
162   - else:
163   - orientation = ''
164   - azimuth = 90.0
165   - elements = member.name.split('-')
166   - zenith = [e for e in elements if 'deg' in e][0].strip('deg')
167   - instrument = '%s_z%s%s_%s' % (site, zenith, orientation, duration)
168   -
169   - # Build output directory
170   - path = '%s/%s' % (outdir, instrument)
171   -
172   - # Write header
173   - self._log_header3(gammalib.TERSE, 'Extract IRFs from "%s"' % member.name)
174   - self._log_value(gammalib.NORMAL, 'Instrument', instrument)
175   -
176   - # Extract file
177   - tar.extract(member, path='%s/%s' % (self['outdir'].filename().url(), path))
178   -
179   - # Create IRF record
180   - record = {'path': path, 'instrument': instrument, 'file': member.name,
181   - 'site': site, 'duration': duration, 'zenith': zenith, 'azimuth': azimuth,
182   - 'orientation': orientation}
183   -
184   - # Add record
185   - irfs.append(record)
186   -
187   - # Close tarfile
188   - tar.close()
189   -
190   - # Return
191   - return irfs
192   -
193   - def _add_irf_to_cif(self, prod, hdu, irf):
194   - """
195   - Add IRF to calibration index file HDU
196   -
197   - Parameters
198   - ----------
199   - prod : str
200   - Production
201   - hdu : ~gammalib.GFitsBinTable
202   - CIF HDU table
203   - irf : dict
204   - IRF dictionary
205   - """
206   - # Set list of IRF component names
207   - #names = ['EA', 'PSF', 'EDISP', 'BKG']
208   - components = [{'name': 'EA', 'CNAM': 'EFF_AREA', 'DESC': 'CTA effective area'},
209   - {'name': 'PSF', 'CNAM': 'RPSF', 'DESC': 'CTA point spread function'},
210   - {'name': 'EDISP', 'CNAM': 'EDISP', 'DESC': 'CTA energy dispersion'},
211   - {'name': 'BKG', 'CNAM': 'BKG', 'DESC': 'CTA background'}]
212   -
213   - # Initialise CIF row index
214   - row = hdu.nrows()
215   -
216   - # Append rows for all components to CIF extension
217   - hdu.append_rows(len(components))
218   -
219   - # Add information for all components
220   - for component in components:
221   -
222   - # Set calibration information
223   - cal_name = 'NAME(%s)' % irf['instrument']
224   - cal_version = 'VERSION(%s)' % prod
225   - cal_cut = 'CLASS(BEST)'
226   - cal_analysis = 'ANALYSIS(CTA)'
227   - cal_zenith = 'ZENITH(%.3f)deg' % float(irf['zenith'])
228   - cal_azimuth = 'AZIMUTH(%.3f)deg' % float(irf['azimuth'])
229   - cal_bounds = [cal_name, cal_version, cal_cut, cal_analysis, \
230   - cal_zenith, cal_azimuth]
231   -
232   - # Set generic information
233   - hdu['TELESCOP'][row] = 'CTA'
234   - hdu['INSTRUME'][row] = gammalib.toupper(prod)
235   - hdu['DETNAM'][row] = 'NONE'
236   - hdu['FILTER'][row] = 'NONE'
237   - hdu['CAL_DEV'][row] = 'ONLINE'
238   - hdu['CAL_CLAS'][row] = 'BCF'
239   - hdu['CAL_DTYP'][row] = 'DATA'
240   - hdu['CAL_VSD'][row] = '2014-01-30'
241   - hdu['CAL_VST'][row] = '00:00:00'
242   - hdu['REF_TIME'][row] = 51544.0
243   - hdu['CAL_QUAL'][row] = 0 # 0=good, 1=bad, 2=dubious, ...
244   - hdu['CAL_DATE'][row] = '14/01/30'
245   -
246   - # Set component specific information
247   - hdu['CAL_DIR'][row] = irf['path']
248   - hdu['CAL_FILE'][row] = irf['file']
249   - hdu['CAL_CNAM'][row] = component['CNAM']
250   - hdu['CAL_DESC'][row] = component['DESC']
251   - hdu['CAL_XNO'][row] = 1
252   - for i in range(9):
253   - if i >= len(cal_bounds):
254   - hdu['CAL_CBD'][row,i] = 'NONE'
255   - else:
256   - hdu['CAL_CBD'][row,i] = cal_bounds[i]
257   -
258   - # Increment row index
259   - row += 1
260   -
261   - # Return
262   - return
263   -
264   - def _create_cif_table(self):
265   - """
266   - Create Calibration Database Index File binary table
267   -
268   - Returns
269   - -------
270   - table : ~gammalib.GFitsBinTable
271   - Calibration Database Index File binary table
272   - """
273   - # Create binary table
274   - table = gammalib.GFitsBinTable()
275   -
276   - # Append columns. Reference: CAL/GEN/92-008
277   - table.append(gammalib.GFitsTableStringCol('TELESCOP', 0, 10))
278   - table.append(gammalib.GFitsTableStringCol('INSTRUME', 0, 10))
279   - table.append(gammalib.GFitsTableStringCol('DETNAM', 0, 20))
280   - table.append(gammalib.GFitsTableStringCol('FILTER', 0, 10))
281   - table.append(gammalib.GFitsTableStringCol('CAL_DEV', 0, 20))
282   - table.append(gammalib.GFitsTableStringCol('CAL_DIR', 0, 70))
283   - table.append(gammalib.GFitsTableStringCol('CAL_FILE', 0, 70)) # Extend beyond standard
284   - table.append(gammalib.GFitsTableStringCol('CAL_CLAS', 0, 3))
285   - table.append(gammalib.GFitsTableStringCol('CAL_DTYP', 0, 4))
286   - table.append(gammalib.GFitsTableStringCol('CAL_CNAM', 0, 20))
287   - table.append(gammalib.GFitsTableStringCol('CAL_CBD', 0, 70, 9))
288   - table.append(gammalib.GFitsTableShortCol('CAL_XNO', 0))
289   - table.append(gammalib.GFitsTableStringCol('CAL_VSD', 0, 10))
290   - table.append(gammalib.GFitsTableStringCol('CAL_VST', 0, 8))
291   - table.append(gammalib.GFitsTableDoubleCol('REF_TIME', 0))
292   - table.append(gammalib.GFitsTableShortCol('CAL_QUAL', 0))
293   - table.append(gammalib.GFitsTableStringCol('CAL_DATE', 0, 8))
294   - table.append(gammalib.GFitsTableStringCol('CAL_DESC', 0, 70))
295   -
296   - # Set keywords. Reference: CAL/GEN/92-008
297   - table.extname('CIF')
298   - table.card('CIFVERSN', '1992a', 'Version of CIF format')
299   -
300   - # Return table
301   - return table
302   -
303   - def _add_irfs_to_cif(self, prod, irfs):
304   - """
305   - Add IRFs to calibration index file
306   -
307   - Parameters
308   - ----------
309   - prod : str
310   - Production
311   - irfs : list
312   - List of IRFs
313   - """
314   - # Build calibration database index name
315   - fname = '%s/data/cta/%s/caldb.indx' % (self['outdir'].filename().url(), prod)
316   -
317   - # Open calibration database index file
318   - cif = gammalib.GFits(fname, True)
319   -
320   - # If file has no CIF extension than create it now
321   - if not cif.contains('CIF'):
322   - cif.append(self._create_cif_table())
323   -
324   - # Get calibration database index HDU
325   - cif_hdu = cif.table('CIF')
326   -
327   - # Loop over all IRFs and add information to CIF
328   - for irf in irfs:
329   - self._add_irf_to_cif(prod, cif_hdu, irf)
330   -
331   - # Save and close CIF
332   - cif.save(True)
333   - cif.close()
334   -
335   - # Return
336   - return
337   -
338   - def _add_irfs(self):
339   - """
340   - Add IRFs to calibration database
341   - """
342   - # Collect tarfiles
343   - tarfiles = self._collect_tarfiles(self['indir'].filename())
344   -
345   - # Add IRFs in tarfiles
346   - for file in tarfiles:
347   -
348   - # Extract Production
349   - fname = os.path.basename(file)
350   - elements = fname.split('-')
351   - prod = [e for e in elements if 'prod' in e][0]
352   - version = [e for e in elements if 'v' in e][0]
353   - prod = '%s-%s' % (prod, version)
354   -
355   - # Build output directory
356   - outdir = 'data/cta/%s/bcf' % (prod)
357   -
358   - # Write header and log parameters
359   - self._log_header2(gammalib.TERSE, 'Extract IRFs')
360   - self._log_value(gammalib.NORMAL, 'Tarfile', file)
361   - self._log_value(gammalib.NORMAL, 'Production', prod)
362   - self._log_value(gammalib.NORMAL, 'Target directory', outdir)
363   -
364   - # Extract IRFs from tarfile
365   - irfs = self._extract_irfs(file, outdir)
366   -
367   - # Add IRFs to calibration index database
368   - self._add_irfs_to_cif(prod, irfs)
369   -
370   - # Return
371   - return
372   -
373   -
374   - # Public methods
375   - def run(self):
376   - """
377   - Run the script
378   - """
379   - # Switch screen logging on in debug mode
380   - if self._logDebug():
381   - self._log.cout(True)
382   -
383   - # Get parameters
384   - self._get_parameters()
385   -
386   - # Write header
387   - self._log_header1(gammalib.TERSE, 'Add IRFs to CALDB')
388   -
389   - # Add IRFs to CALDB
390   - self._add_irfs()
391   -
392   - # Return
393   - return
394   -
395   -
396   -# ======================== #
397   -# Main routine entry point #
398   -# ======================== #
399   -if __name__ == '__main__':
400   -
401   - # Create instance of application
402   - app = csadd2caldb(sys.argv)
403   -
404   - # Execute application
405   - app.execute()
doc/source/users/user_manual/csadd2caldb.py
  View file @c057d37
... ... @@ -67,7 +67,7 @@ class csadd2caldb(ctools.cscript):
67 67  
68 68 # Create default parfile
69 69 pars = gammalib.GApplicationPars()
70   - pars.append(gammalib.GApplicationPar('indir','f','a','cta-prod5-zenodo-fitsonly-v0','','','Input IRF folder'))
  70 + pars.append(gammalib.GApplicationPar('indir','f','a','fits','','','Input IRF folder'))
71 71 pars.append(gammalib.GApplicationPar('outdir','f','a','$CALDB','','','Output caldb folder'))
72 72 pars.append_standard()
73 73 pars.save(parfile)
... ... @@ -110,7 +110,7 @@ class csadd2caldb(ctools.cscript):
110 110 files = glob.glob('%s/*.FITS.tar.gz' % indir.url())
111 111 tarfiles.extend(files)
112 112  
113   - # Search for tarfiles in 'files' subdirectory
  113 + # Search for tarfiles in 'fits' subdirectory
114 114 files = glob.glob('%s/fits/*.FITS.tar.gz' % indir.url())
115 115 tarfiles.extend(files)
116 116  
... ... @@ -142,8 +142,23 @@ class csadd2caldb(ctools.cscript):
142 142 # Loop over all members in tarfile
143 143 for member in tar.getmembers():
144 144  
  145 + # Get subarray
  146 + subarrays = ''
  147 + nsubarrays = 0
  148 + if 'LST' in member.name:
  149 + subarrays += '_LST'
  150 + nsubarrays += 1
  151 + if 'MST' in member.name:
  152 + subarrays += '_MST'
  153 + nsubarrays += 1
  154 + if 'SST' in member.name:
  155 + subarrays += '_SST'
  156 + nsubarrays += 1
  157 + if nsubarrays > 1:
  158 + subarrays = ''
  159 +
145 160 # Extract Instrument
146   - if 'South' in member.name:
  161 + if 'South-' in member.name:
147 162 site = 'South'
148 163 else:
149 164 site = 'North'
... ... @@ -162,9 +177,10 @@ class csadd2caldb(ctools.cscript):
162 177 else:
163 178 orientation = ''
164 179 azimuth = 90.0
  180 +
165 181 elements = member.name.split('-')
166 182 zenith = [e for e in elements if 'deg' in e][0].strip('deg')
167   - instrument = '%s_z%s%s_%s' % (site, zenith, orientation, duration)
  183 + instrument = '%s_z%s%s_%s%s' % (site, zenith, orientation, duration, subarrays)
168 184  
169 185 # Build output directory
170 186 path = '%s/%s' % (outdir, instrument)
... ...
doc/source/users/user_manual/irf_cta.rst
  View file @c057d37
... ... @@ -80,6 +80,52 @@ In total, the following six instrument response functions are available:
80 80 ``North_0.5h``, ``North_5h``, ``North_50h``, ``South_0.5h``,
81 81 ``South_5h``, and ``South_50h``.
82 82  
  83 +The latest CTA reponse functions (``prod5-v0.1``) can be downloaded from
  84 +`this link <https://zenodo.org/record/5499840#.YVV1di0itZo>`_.
  85 +
  86 +Addition of CTA response functions to the calibration database can be done
  87 +using the script :download:`csadd2caldb.py <csadd2caldb.py>` that is available
  88 +for download.
  89 +
  90 +To download and install the ``prod5-v0.1`` response functions you should execute
  91 +the following commands in a terminal:
  92 +
  93 +.. code-block:: bash
  94 +
  95 + $ mkdir temp
  96 + $ cd temp
  97 + $ wget https://zenodo.org/record/5499840/files/cta-prod5-zenodo-fitsonly-v0.1.zip
  98 + $ unzip cta-prod5-zenodo-fitsonly-v0.1.zip
  99 + $ ./csadd2caldb.py
  100 + Parfile csadd2caldb.par not found. Create default parfile.
  101 + Input IRF folder [fits]
  102 + Output caldb folder [$CALDB]
  103 +
  104 +This will download the file ``cta-prod5-zenodo-fitsonly-v0.1.zip`` in a
  105 +temporary working directory, unzip the file to get, among others, the folder ``fits``
  106 +with the CTA response in FITS format, and install the response in your calibration
  107 +database.
  108 +
  109 +Use ``sudo`` to run ``csadd2caldb.py`` if the installation of the IRFs requires
  110 +root privileges.
  111 +
  112 +To use the ``prod5-v0.1`` response functions, specify ``prod5-v0.1`` if a ctool
  113 +or cscript queries for the calibration database.
  114 +
  115 +If you want to install the older ``prod3b-v2`` response functions you should use
  116 +the following procedure:
  117 +
  118 +.. code-block:: bash
  119 +
  120 + $ mkdir temp
  121 + $ cd temp
  122 + $ wget https://zenodo.org/record/5163273/files/CTA-Performance-IRFs-prod3b-v2-v1.0.0.zip
  123 + $ unzip CTA-Performance-IRFs-prod3b-v2-v1.0.0.zip
  124 + $ tar xvfz CTA-Prod3-Zenodo-main/fits/CTA-Performance-prod3b-v2-FITS.tar.gz -C $CTOOLS/share
  125 +
  126 +To use the ``prod3b-v2`` response functions, specify ``prod3b-v2`` if a ctool
  127 +or cscript queries for the calibration database.
  128 +
83 129 .. note::
84 130 If you use the CTA instrument response functions for a research project
85 131 and specifically for a presentation or publication, we ask to add the
... ... @@ -88,29 +134,6 @@ In total, the following six instrument response functions are available:
88 134 *This research has made use of the CTA instrument response functions
89 135 provided by the CTA Consortium and Observatory.*
90 136  
91   - The latest CTA reponse functions can be downloaded from
92   - `this link <https://zenodo.org/record/5499840/files/cta-prod5-zenodo-fitsonly-v0.1.zip?download=1>`_.
93   -
94   - This will download the file ``cta-prod5-zenodo-fitsonly-v0.1.zip`` in your
95   - current working directory that you should unzip to get a folder named
96   - ``cta-prod5-zenodo-fitsonly-v0``.
97   -
98   - Download the script :download:`csadd2caldb.py <csadd2caldb.py>` and execute
99   - it as follows to add the IRFs to your calibration database:
100   -
101   - .. code-block:: bash
102   -
103   - $ ./csadd2caldb.py debug=yes
104   - Parfile csadd2caldb.par not found. Create default parfile.
105   - Input IRF folder [cta-prod5-zenodo-fitsonly-v0]
106   - Output caldb folder [$CALDB]
107   -
108   - Use ``sudo`` to run ``csadd2caldb.py`` if the installation of the IRFs requires
109   - root privileges.
110   -
111   - To use the prod5-v0.1 functions, specify ``prod5-v0.1`` if a tool
112   - or script queries for the calibration database.
113   -
114 137  
115 138 Specifying CTA response functions
116 139 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
... ...