Commit 2ce711e85ea6dfe4e303515349ce9ac8cfd3059d
1 parent
081224fa
initial version
Showing
1 changed file
with
606 additions
and
0 deletions
scripts/cta_root2caldb.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ========================================================================== | ||
| 3 | +# This script generates IRFs in the CALDB format of HEASARC using ROOT 2D | ||
| 4 | +# performance files. | ||
| 5 | +# | ||
| 6 | +# Requirements: | ||
| 7 | +# - pyroot | ||
| 8 | +# | ||
| 9 | +# Copyright (C) 2011 Jurgen Knodlseder | ||
| 10 | +# | ||
| 11 | +# This program is free software: you can redistribute it and/or modify | ||
| 12 | +# it under the terms of the GNU General Public License as published by | ||
| 13 | +# the Free Software Foundation, either version 3 of the License, or | ||
| 14 | +# (at your option) any later version. | ||
| 15 | +# | ||
| 16 | +# This program is distributed in the hope that it will be useful, | ||
| 17 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 18 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 19 | +# GNU General Public License for more details. | ||
| 20 | +# | ||
| 21 | +# You should have received a copy of the GNU General Public License | ||
| 22 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 23 | +# | ||
| 24 | +# ========================================================================== | ||
| 25 | +from ROOT import TFile, TH1F, TH2F | ||
| 26 | +from gammalib import * | ||
| 27 | +from datetime import datetime | ||
| 28 | +import math | ||
| 29 | +import sys | ||
| 30 | +import glob | ||
| 31 | +import os | ||
| 32 | + | ||
| 33 | + | ||
| 34 | +# =========== # | ||
| 35 | +# CalDB class # | ||
| 36 | +# =========== # | ||
| 37 | +class caldb(): | ||
| 38 | + """ | ||
| 39 | + """ | ||
| 40 | + def __init__(self, conf, path=""): | ||
| 41 | + """ | ||
| 42 | + Constructor. | ||
| 43 | + """ | ||
| 44 | + # Store root path | ||
| 45 | + self.path = path | ||
| 46 | + | ||
| 47 | + # Initialise some members | ||
| 48 | + self.cif = None | ||
| 49 | + self.ea = None | ||
| 50 | + self.psf = None | ||
| 51 | + self.edisp = None | ||
| 52 | + self.hdu_cif = None | ||
| 53 | + | ||
| 54 | + # Store UTC date of job execution | ||
| 55 | + self.utc = datetime.utcnow().isoformat()[:19] | ||
| 56 | + | ||
| 57 | + # Set CALDB information | ||
| 58 | + self.cal_tel = "CTA" | ||
| 59 | + self.cal_inst = conf.upper() | ||
| 60 | + self.cal_det = "NONE" | ||
| 61 | + self.cal_flt = "NONE" | ||
| 62 | + self.cal_class = "BCF" | ||
| 63 | + self.cal_type = "DATA" | ||
| 64 | + self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ... | ||
| 65 | + self.cal_date = "11/01/01" | ||
| 66 | + self.val_date = "2011-01-01" | ||
| 67 | + self.val_time = "00:00:00" | ||
| 68 | + self.ref_time = 51544.0 | ||
| 69 | + self.cal_version = "VERSION(DC1)" | ||
| 70 | + self.cal_cut = "CLASS(BEST)" | ||
| 71 | + self.cal_analysis = "ANALYSIS(IFAE)" | ||
| 72 | + # | ||
| 73 | + self.ea_name = "EFF_AREA" | ||
| 74 | + self.ea_doc = "CAL/GEN/92-019" | ||
| 75 | + self.ea_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] | ||
| 76 | + self.ea_desc = "CTA effective area" | ||
| 77 | + # | ||
| 78 | + self.psf_name = "RPSF" | ||
| 79 | + self.psf_doc = "CAL/GEN/92-020" | ||
| 80 | + self.psf_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] | ||
| 81 | + self.psf_desc = "CTA point spread function" | ||
| 82 | + # | ||
| 83 | + self.edisp_name = "EDISP" | ||
| 84 | + self.edisp_doc = "???" | ||
| 85 | + self.edisp_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] | ||
| 86 | + self.edisp_desc = "CTA energy dispersion" | ||
| 87 | + | ||
| 88 | + # Create directory structure | ||
| 89 | + self.make_dirs() | ||
| 90 | + | ||
| 91 | + # Return | ||
| 92 | + return | ||
| 93 | + | ||
| 94 | + def __del__(self): | ||
| 95 | + """ | ||
| 96 | + """ | ||
| 97 | + # Close calibration | ||
| 98 | + self.close() | ||
| 99 | + | ||
| 100 | + # Return | ||
| 101 | + return | ||
| 102 | + | ||
| 103 | + def make_dirs(self): | ||
| 104 | + """ | ||
| 105 | + """ | ||
| 106 | + # Set base path | ||
| 107 | + self.base_dir = "data" | ||
| 108 | + self.base_dir += "/"+self.cal_tel.lower() | ||
| 109 | + self.base_dir += "/"+self.cal_inst.lower() | ||
| 110 | + | ||
| 111 | + | ||
| 112 | + # Set directory names | ||
| 113 | + self.bcf_dir = self.base_dir+"/bcf" | ||
| 114 | + self.ea_dir = self.bcf_dir+"/ea" | ||
| 115 | + self.psf_dir = self.bcf_dir+"/psf" | ||
| 116 | + self.edisp_dir = self.bcf_dir+"/edisp" | ||
| 117 | + | ||
| 118 | + # Add path (if required) | ||
| 119 | + if len(self.path) > 0: | ||
| 120 | + self.bcf_path = self.path+"/"+self.bcf_dir | ||
| 121 | + self.ea_path = self.path+"/"+self.ea_dir | ||
| 122 | + self.psf_path = self.path+"/"+self.psf_dir | ||
| 123 | + self.edisp_path = self.path+"/"+self.edisp_dir | ||
| 124 | + else: | ||
| 125 | + self.bcf_path = self.bcf_dir | ||
| 126 | + self.ea_path = self.ea_dir | ||
| 127 | + self.psf_path = self.psf_dir | ||
| 128 | + self.edisp_path = self.edisp_dir | ||
| 129 | + | ||
| 130 | + # Create directory structure | ||
| 131 | + if not os.path.isdir(self.ea_path): | ||
| 132 | + os.makedirs(self.ea_path) | ||
| 133 | + if not os.path.isdir(self.psf_path): | ||
| 134 | + os.makedirs(self.psf_path) | ||
| 135 | + if not os.path.isdir(self.edisp_path): | ||
| 136 | + os.makedirs(self.edisp_path) | ||
| 137 | + | ||
| 138 | + # Return | ||
| 139 | + return | ||
| 140 | + | ||
| 141 | + def open(self, name, clobber=True): | ||
| 142 | + """ | ||
| 143 | + Open existing or create new calibration. | ||
| 144 | + """ | ||
| 145 | + # Set calibrate file names | ||
| 146 | + self.ea_file = "ea_"+name+".fits" | ||
| 147 | + self.psf_file = "psf_"+name+".fits" | ||
| 148 | + self.edisp_file = "edisp_"+name+".fits" | ||
| 149 | + | ||
| 150 | + # Open calibration database index | ||
| 151 | + self.cif = GFits(self.bcf_path+"/caldb.indx", True) | ||
| 152 | + | ||
| 153 | + # If file has no CIF extension than create it now | ||
| 154 | + try: | ||
| 155 | + self.hdu_cif = self.cif.table("CIF") | ||
| 156 | + except: | ||
| 157 | + self.create_cif() | ||
| 158 | + self.hdu_cif = self.cif.table("CIF") | ||
| 159 | + | ||
| 160 | + # Open calibration files | ||
| 161 | + self.open_ea(self.ea_path+"/"+self.ea_file) | ||
| 162 | + self.open_psf(self.psf_path+"/"+self.psf_file) | ||
| 163 | + self.open_edisp(self.edisp_path+"/"+self.edisp_file) | ||
| 164 | + | ||
| 165 | + # Return | ||
| 166 | + return | ||
| 167 | + | ||
| 168 | + def close(self): | ||
| 169 | + """ | ||
| 170 | + Close calibration. | ||
| 171 | + """ | ||
| 172 | + # Add information to CIF. We do this now as before we did not | ||
| 173 | + # necessarily have all the information at hand (in particular | ||
| 174 | + # about the boundaries) | ||
| 175 | + self.add_cif_info() | ||
| 176 | + | ||
| 177 | + # Close CIF | ||
| 178 | + if self.cif != None: | ||
| 179 | + self.cif.save(True) | ||
| 180 | + self.cif.close() | ||
| 181 | + self.cif = None | ||
| 182 | + | ||
| 183 | + # Close effective area | ||
| 184 | + if self.ea != None: | ||
| 185 | + self.ea.save(True) | ||
| 186 | + self.ea.close() | ||
| 187 | + self.ea = None | ||
| 188 | + | ||
| 189 | + # Close point spread function | ||
| 190 | + if self.psf != None: | ||
| 191 | + self.psf.save(True) | ||
| 192 | + self.psf.close() | ||
| 193 | + self.psf = None | ||
| 194 | + | ||
| 195 | + # Close energy dispersion | ||
| 196 | + if self.edisp != None: | ||
| 197 | + self.edisp.save(True) | ||
| 198 | + self.edisp.close() | ||
| 199 | + self.edisp = None | ||
| 200 | + | ||
| 201 | + # Return | ||
| 202 | + return | ||
| 203 | + | ||
| 204 | + def create_cif(self): | ||
| 205 | + """ | ||
| 206 | + Create CIF extension in FITS file. | ||
| 207 | + """ | ||
| 208 | + # Create binary table | ||
| 209 | + table = GFitsBinTable() | ||
| 210 | + | ||
| 211 | + # Set boundary | ||
| 212 | + | ||
| 213 | + # Attach columns. Reference: CAL/GEN/92-008 | ||
| 214 | + table.append_column(GFitsTableStringCol("TELESCOP", 0, 10)) | ||
| 215 | + table.append_column(GFitsTableStringCol("INSTRUME", 0, 10)) | ||
| 216 | + table.append_column(GFitsTableStringCol("DETNAM", 0, 20)) | ||
| 217 | + table.append_column(GFitsTableStringCol("FILTER", 0, 10)) | ||
| 218 | + table.append_column(GFitsTableStringCol("CAL_DEV", 0, 20)) | ||
| 219 | + table.append_column(GFitsTableStringCol("CAL_DIR", 0, 70)) | ||
| 220 | + table.append_column(GFitsTableStringCol("CAL_FILE", 0, 40)) | ||
| 221 | + table.append_column(GFitsTableStringCol("CAL_CLAS", 0, 3)) | ||
| 222 | + table.append_column(GFitsTableStringCol("CAL_DTYP", 0, 4)) | ||
| 223 | + table.append_column(GFitsTableStringCol("CAL_CNAM", 0, 20)) | ||
| 224 | + table.append_column(GFitsTableStringCol("CAL_CBD", 0, 70, 9)) | ||
| 225 | + table.append_column(GFitsTableShortCol("CAL_XNO", 0)) | ||
| 226 | + table.append_column(GFitsTableStringCol("CAL_VSD", 0, 10)) | ||
| 227 | + table.append_column(GFitsTableStringCol("CAL_VST", 0, 8)) | ||
| 228 | + table.append_column(GFitsTableDoubleCol("REF_TIME", 0)) | ||
| 229 | + table.append_column(GFitsTableShortCol("CAL_QUAL", 0)) | ||
| 230 | + table.append_column(GFitsTableStringCol("CAL_DATE", 0, 8)) | ||
| 231 | + table.append_column(GFitsTableStringCol("CAL_DESC", 0, 70)) | ||
| 232 | + | ||
| 233 | + # Set keywords. Reference: CAL/GEN/92-008 | ||
| 234 | + table.extname("CIF") | ||
| 235 | + table.card("CIFVERSN", "1992a", "Version of CIF format") | ||
| 236 | + | ||
| 237 | + # Attach table to FITS file. Note that at this moment the | ||
| 238 | + # FITS table gets connected to the FITS file. Yet since | ||
| 239 | + # nothing was yet written to the FITS file, we cannot read | ||
| 240 | + # anything from it. | ||
| 241 | + self.cif.append(table) | ||
| 242 | + | ||
| 243 | + # Return | ||
| 244 | + return | ||
| 245 | + | ||
| 246 | + def add_cif_info(self): | ||
| 247 | + """ | ||
| 248 | + Add information to CIF extension. | ||
| 249 | + """ | ||
| 250 | + # Append 3 rows to CIF extension | ||
| 251 | + self.hdu_cif.append_rows(3) | ||
| 252 | + | ||
| 253 | + # Add generic information for these 3 rows | ||
| 254 | + for i in range(3): | ||
| 255 | + | ||
| 256 | + # Set row number | ||
| 257 | + row = i+self.hdu_cif.nrows()-3 | ||
| 258 | + | ||
| 259 | + # Set element | ||
| 260 | + self.hdu_cif["TELESCOP"][row] = self.cal_tel | ||
| 261 | + self.hdu_cif["INSTRUME"][row] = self.cal_inst | ||
| 262 | + self.hdu_cif["DETNAM"][row] = self.cal_det | ||
| 263 | + self.hdu_cif["FILTER"][row] = self.cal_flt | ||
| 264 | + self.hdu_cif["CAL_DEV"][row] = "ONLINE" | ||
| 265 | + self.hdu_cif["CAL_CLAS"][row] = self.cal_class | ||
| 266 | + self.hdu_cif["CAL_DTYP"][row] = self.cal_type | ||
| 267 | + self.hdu_cif["CAL_VSD"][row] = self.val_date | ||
| 268 | + self.hdu_cif["CAL_VST"][row] = self.val_time | ||
| 269 | + self.hdu_cif["REF_TIME"][row] = self.ref_time | ||
| 270 | + self.hdu_cif["CAL_QUAL"][row] = self.cal_qual | ||
| 271 | + self.hdu_cif["CAL_DATE"][row] = self.cal_date | ||
| 272 | + | ||
| 273 | + # Add effective area information | ||
| 274 | + row = self.hdu_cif.nrows()-3 | ||
| 275 | + self.hdu_cif["CAL_DIR"][row] = self.ea_dir | ||
| 276 | + self.hdu_cif["CAL_FILE"][row] = self.ea_file | ||
| 277 | + self.hdu_cif["CAL_CNAM"][row] = self.ea_name | ||
| 278 | + self.hdu_cif["CAL_DESC"][row] = self.ea_desc | ||
| 279 | + self.hdu_cif["CAL_XNO"][row] = 1 | ||
| 280 | + for n in range(9): | ||
| 281 | + if n >= len(self.ea_bounds): | ||
| 282 | + self.hdu_cif["CAL_CBD"][row, n] = "NONE" | ||
| 283 | + else: | ||
| 284 | + self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n] | ||
| 285 | + | ||
| 286 | + # Add point spread function information | ||
| 287 | + row = self.hdu_cif.nrows()-2 | ||
| 288 | + self.hdu_cif["CAL_DIR"][row] = self.psf_dir | ||
| 289 | + self.hdu_cif["CAL_FILE"][row] = self.psf_file | ||
| 290 | + self.hdu_cif["CAL_CNAM"][row] = self.psf_name | ||
| 291 | + self.hdu_cif["CAL_DESC"][row] = self.psf_desc | ||
| 292 | + self.hdu_cif["CAL_XNO"][row] = 1 | ||
| 293 | + for n in range(9): | ||
| 294 | + if n >= len(self.psf_bounds): | ||
| 295 | + self.hdu_cif["CAL_CBD"][row, n] = "NONE" | ||
| 296 | + else: | ||
| 297 | + self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n] | ||
| 298 | + | ||
| 299 | + # Add energy dispersion information | ||
| 300 | + row = self.hdu_cif.nrows()-1 | ||
| 301 | + self.hdu_cif["CAL_DIR"][row] = self.edisp_dir | ||
| 302 | + self.hdu_cif["CAL_FILE"][row] = self.edisp_file | ||
| 303 | + self.hdu_cif["CAL_CNAM"][row] = self.edisp_name | ||
| 304 | + self.hdu_cif["CAL_DESC"][row] = self.edisp_desc | ||
| 305 | + self.hdu_cif["CAL_XNO"][row] = 1 | ||
| 306 | + for n in range(9): | ||
| 307 | + if n >= len(self.edisp_bounds): | ||
| 308 | + self.hdu_cif["CAL_CBD"][row, n] = "NONE" | ||
| 309 | + else: | ||
| 310 | + self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n] | ||
| 311 | + | ||
| 312 | + # Return | ||
| 313 | + return | ||
| 314 | + | ||
| 315 | + def open_ea(self, filename): | ||
| 316 | + """ | ||
| 317 | + Open effective area extension. If no extension was found | ||
| 318 | + then create one. We do not yet set any data as we don't | ||
| 319 | + know the dimensions of the table yet. | ||
| 320 | + """ | ||
| 321 | + # Open FITS file | ||
| 322 | + self.ea = GFits(filename, True) | ||
| 323 | + | ||
| 324 | + # Get extensions. If they do not exist then create them now | ||
| 325 | + try: | ||
| 326 | + self.hdu_ea = self.ea.table("EFFECTIVE AREA") | ||
| 327 | + except: | ||
| 328 | + # Create binary table | ||
| 329 | + table = GFitsBinTable() | ||
| 330 | + | ||
| 331 | + # Set extension name | ||
| 332 | + table.extname("EFFECTIVE AREA") | ||
| 333 | + | ||
| 334 | + # Set OGIP keywords | ||
| 335 | + self.set_ogip_keywords(table, self.ea_doc, \ | ||
| 336 | + ["RESPONSE", self.ea_name]) | ||
| 337 | + | ||
| 338 | + # Append table | ||
| 339 | + self.ea.append(table) | ||
| 340 | + | ||
| 341 | + # Get extension | ||
| 342 | + self.hdu_ea = self.ea.table("EFFECTIVE AREA") | ||
| 343 | + | ||
| 344 | + # Return | ||
| 345 | + return | ||
| 346 | + | ||
| 347 | + def open_psf(self, filename): | ||
| 348 | + """ | ||
| 349 | + Open point spread function extension. If no extension was | ||
| 350 | + found then create one. We do not yet set any data as we don't | ||
| 351 | + know the dimensions of the table yet. | ||
| 352 | + """ | ||
| 353 | + # Open FITS file | ||
| 354 | + self.psf = GFits(filename, True) | ||
| 355 | + | ||
| 356 | + # Get extensions. If they do not exist then create them now | ||
| 357 | + try: | ||
| 358 | + self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") | ||
| 359 | + except: | ||
| 360 | + # Create binary table | ||
| 361 | + table = GFitsBinTable() | ||
| 362 | + | ||
| 363 | + # Set extension name | ||
| 364 | + table.extname("POINT SPREAD FUNCTION") | ||
| 365 | + | ||
| 366 | + # Set OGIP keywords | ||
| 367 | + self.set_ogip_keywords(table, self.psf_doc, \ | ||
| 368 | + ["RESPONSE", self.psf_name]) | ||
| 369 | + | ||
| 370 | + # Append table | ||
| 371 | + self.psf.append(table) | ||
| 372 | + | ||
| 373 | + # Get extension | ||
| 374 | + self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") | ||
| 375 | + | ||
| 376 | + # Return | ||
| 377 | + return | ||
| 378 | + | ||
| 379 | + def open_edisp(self, filename): | ||
| 380 | + """ | ||
| 381 | + Open energy dispersion extension. If no extension was found | ||
| 382 | + then create one. We do not yet set any data as we don't | ||
| 383 | + know the dimensions of the table yet. | ||
| 384 | + """ | ||
| 385 | + # Open FITS file | ||
| 386 | + self.edisp = GFits(filename, True) | ||
| 387 | + | ||
| 388 | + # Get extensions. If they do not exist then create them now | ||
| 389 | + try: | ||
| 390 | + self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") | ||
| 391 | + except: | ||
| 392 | + # Create binary table | ||
| 393 | + table = GFitsBinTable() | ||
| 394 | + | ||
| 395 | + # Set extension name | ||
| 396 | + table.extname("ENERGY DISPERSION") | ||
| 397 | + | ||
| 398 | + # Set OGIP keywords | ||
| 399 | + self.set_ogip_keywords(table, self.edisp_doc, \ | ||
| 400 | + ["RESPONSE", self.edisp_name]) | ||
| 401 | + | ||
| 402 | + # Append table | ||
| 403 | + self.edisp.append(table) | ||
| 404 | + | ||
| 405 | + # Get extension | ||
| 406 | + self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") | ||
| 407 | + | ||
| 408 | + # Return | ||
| 409 | + return | ||
| 410 | + | ||
| 411 | + def set_ogip_keywords(self, hdu, hdudoc, hduclas): | ||
| 412 | + """ | ||
| 413 | + Set standard OGIP keywords for extension. | ||
| 414 | + """ | ||
| 415 | + # Set keywords | ||
| 416 | + hdu.card("ORIGIN", "IRAP", "Name of organization making this file") | ||
| 417 | + hdu.card("DATE", self.utc, "File creation date (YYYY-MM-DDThh:mm:ss UTC)") | ||
| 418 | + hdu.card("TELESCOP", self.cal_tel, "Name of telescope") | ||
| 419 | + hdu.card("INSTRUME", self.cal_inst, "Name of instrument") | ||
| 420 | + hdu.card("DETNAM", self.cal_det, "Name of detector") | ||
| 421 | + hdu.card("HDUCLASS", "OGIP", "") | ||
| 422 | + hdu.card("HDUDOC", hdudoc, "") | ||
| 423 | + for i, item in enumerate(hduclas): | ||
| 424 | + key = "HDUCLAS%d" % (i+1) | ||
| 425 | + hdu.card(key, item, "") | ||
| 426 | + hdu.card("HDUVERS", "1.0.0", "") | ||
| 427 | + | ||
| 428 | + # Return | ||
| 429 | + return | ||
| 430 | + | ||
| 431 | + def make_2D(self, array, hdu, name, unit, scale=1.0): | ||
| 432 | + """ | ||
| 433 | + Make 2D matrix as function of energy and offset angle from a | ||
| 434 | + ROOT 2D histogram. | ||
| 435 | + | ||
| 436 | + Parameters: | ||
| 437 | + array - ROOT 2D histgram. | ||
| 438 | + hdu - FITS HDU. | ||
| 439 | + Keyword: | ||
| 440 | + scale - Scaling factor for histogram values. | ||
| 441 | + """ | ||
| 442 | + # Extract energy and offset angle vectors | ||
| 443 | + energies = array.GetXaxis() | ||
| 444 | + offsets = array.GetYaxis() | ||
| 445 | + neng = energies.GetNbins() | ||
| 446 | + noffset = offsets.GetNbins() | ||
| 447 | + | ||
| 448 | + # Attach columns to HDU. We do this only if they do not yet | ||
| 449 | + # exist. This allows adding columns to the same HDU in successive | ||
| 450 | + # calls of this method | ||
| 451 | + # | ||
| 452 | + # ENERG_LO | ||
| 453 | + if not hdu.has("ENERG_LO"): | ||
| 454 | + hdu.append_column(GFitsTableFloatCol("ENERG_LO", 1, neng)) | ||
| 455 | + hdu["ENERG_LO"].unit("TeV") | ||
| 456 | + for ieng in range(neng): | ||
| 457 | + e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1)) | ||
| 458 | + hdu["ENERG_LO"][0,ieng] = e_lo | ||
| 459 | + # | ||
| 460 | + # ENERG_HI | ||
| 461 | + if not hdu.has("ENERG_HI"): | ||
| 462 | + hdu.append_column(GFitsTableFloatCol("ENERG_HI", 1, neng)) | ||
| 463 | + hdu["ENERG_HI"].unit("TeV") | ||
| 464 | + for ieng in range(neng): | ||
| 465 | + e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1)) | ||
| 466 | + hdu["ENERG_HI"][0,ieng] = e_hi | ||
| 467 | + # | ||
| 468 | + # THETA_LO | ||
| 469 | + if not hdu.has("THETA_LO"): | ||
| 470 | + hdu.append_column(GFitsTableFloatCol("THETA_LO", 1, noffset)) | ||
| 471 | + hdu["THETA_LO"].unit("deg") | ||
| 472 | + for ioff in range(noffset): | ||
| 473 | + o_lo = offsets.GetBinLowEdge(ioff+1) | ||
| 474 | + hdu["THETA_LO"][0,ioff] = o_lo | ||
| 475 | + | ||
| 476 | + # | ||
| 477 | + # THETA_LO | ||
| 478 | + if not hdu.has("THETA_HI"): | ||
| 479 | + hdu.append_column(GFitsTableFloatCol("THETA_HI", 1, noffset)) | ||
| 480 | + hdu["THETA_HI"].unit("deg") | ||
| 481 | + for ioff in range(noffset): | ||
| 482 | + o_hi = offsets.GetBinUpEdge(ioff+1) | ||
| 483 | + hdu["THETA_HI"][0,ioff] = o_hi | ||
| 484 | + # | ||
| 485 | + # "NAME" | ||
| 486 | + if not hdu.has(name): | ||
| 487 | + hdu.append_column(GFitsTableFloatCol(name, 1, neng*noffset)) | ||
| 488 | + hdu[name].unit(unit) | ||
| 489 | + hdu[name].dim([neng, noffset]) | ||
| 490 | + for ioff in range(noffset): | ||
| 491 | + for ieng in range(neng): | ||
| 492 | + index = ieng + ioff * neng | ||
| 493 | + value = array.GetBinContent(ieng+1,ioff+1) | ||
| 494 | + hdu[name][0,index] = value * scale | ||
| 495 | + | ||
| 496 | + # Collect boundary information | ||
| 497 | + bd_eng = "ENERG(%.4f-%.2f)TeV" % \ | ||
| 498 | + (pow(10.0, energies.GetBinLowEdge(1)), \ | ||
| 499 | + pow(10.0, energies.GetBinUpEdge(neng))) | ||
| 500 | + bd_off = "THETA(%.2f-%.2f)deg" % \ | ||
| 501 | + (offsets.GetBinLowEdge(1), \ | ||
| 502 | + offsets.GetBinUpEdge(noffset)) | ||
| 503 | + bd_phi = "PHI(0-360)deg" | ||
| 504 | + bounds = [bd_eng, bd_off, bd_phi] | ||
| 505 | + | ||
| 506 | + # Return boundary information | ||
| 507 | + return bounds | ||
| 508 | + | ||
| 509 | + def root2caldb(self, filename): | ||
| 510 | + """ | ||
| 511 | + Translate ROOT to CALDB information. | ||
| 512 | + | ||
| 513 | + Parameters: | ||
| 514 | + filename - ROOT 2D performance filename. | ||
| 515 | + """ | ||
| 516 | + # Open performance file | ||
| 517 | + file = TFile(filename) | ||
| 518 | + | ||
| 519 | + # Allocate ROOT 2D array | ||
| 520 | + array = TH2F() | ||
| 521 | + | ||
| 522 | + # Make 2D effective 80% area | ||
| 523 | + if self.hdu_ea != None: | ||
| 524 | + file.GetObject("EffectiveArea80_offaxis", array) | ||
| 525 | + bounds = self.make_2D(array, self.hdu_ea, "EFFAREA", "m2", scale=1.0/0.8) | ||
| 526 | + for b in bounds: | ||
| 527 | + self.ea_bounds.append(b) | ||
| 528 | + self.set_cif_keywords(self.hdu_ea, self.ea_name, \ | ||
| 529 | + self.ea_bounds, self.ea_desc) | ||
| 530 | + | ||
| 531 | + # Make 2D 68% point spread function | ||
| 532 | + if self.hdu_psf != None: | ||
| 533 | + file.GetObject("AngRes_offaxis", array) | ||
| 534 | + bounds = self.make_2D(array, self.hdu_psf, "R68", "deg") | ||
| 535 | + for b in bounds: | ||
| 536 | + self.psf_bounds.append(b) | ||
| 537 | + self.set_cif_keywords(self.hdu_psf, self.psf_name, \ | ||
| 538 | + self.psf_bounds, self.psf_desc) | ||
| 539 | + | ||
| 540 | + # Make 2D 80% point spread function | ||
| 541 | + if self.hdu_psf != None: | ||
| 542 | + file.GetObject("AngRes80_offaxis", array) | ||
| 543 | + self.make_2D(array, self.hdu_psf, "R80", "deg") | ||
| 544 | + | ||
| 545 | + # Make 2D energy resolution | ||
| 546 | + if self.hdu_edisp != None: | ||
| 547 | + file.GetObject("ERes", array) | ||
| 548 | + bounds = self.make_2D(array, self.hdu_edisp, "ERES68", "") | ||
| 549 | + for b in bounds: | ||
| 550 | + self.edisp_bounds.append(b) | ||
| 551 | + self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \ | ||
| 552 | + self.edisp_bounds, self.edisp_desc) | ||
| 553 | + | ||
| 554 | + # Return | ||
| 555 | + return | ||
| 556 | + | ||
| 557 | + def set_cif_keywords(self, hdu, name, bounds,desc): | ||
| 558 | + """ | ||
| 559 | + Set standard CIF keywords for extension. | ||
| 560 | + """ | ||
| 561 | + # Set standard CIF keywords | ||
| 562 | + hdu.card("CSYSNAME", "XMA_POL", "") | ||
| 563 | + hdu.card("CCLS0001", self.cal_class, "Calibration class") | ||
| 564 | + hdu.card("CDTP0001", self.cal_type, "Calibration type") | ||
| 565 | + hdu.card("CCNM0001", name, "Calibration name") | ||
| 566 | + | ||
| 567 | + # Set boundary keywords | ||
| 568 | + for n in range(9): | ||
| 569 | + keyname = "CBD%d0001" % (n+1) | ||
| 570 | + if n >= len(bounds): | ||
| 571 | + value = "NONE" | ||
| 572 | + else: | ||
| 573 | + value = bounds[n] | ||
| 574 | + hdu.card(keyname, value, "Calibration boundary") | ||
| 575 | + | ||
| 576 | + # Set validity keywords | ||
| 577 | + hdu.card("CVSD0001", self.val_date, "Calibration validity start date (UTC)") | ||
| 578 | + hdu.card("CVST0001", self.val_time, "Calibration validity start time (UTC)") | ||
| 579 | + hdu.card("CDEC0001", desc, "Calibration description") | ||
| 580 | + | ||
| 581 | + # Return | ||
| 582 | + return | ||
| 583 | + | ||
| 584 | + | ||
| 585 | +#==========================# | ||
| 586 | +# Main routine entry point # | ||
| 587 | +#==========================# | ||
| 588 | +if __name__ == '__main__': | ||
| 589 | + """ | ||
| 590 | + This script translates a ROOT 2D performance file in a CALDB compliant | ||
| 591 | + CTA response. | ||
| 592 | + """ | ||
| 593 | + # Set ROOT filename | ||
| 594 | + path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxis" | ||
| 595 | + irf = "SubarrayE_offaxis.root" | ||
| 596 | + filename = path+"/"+irf | ||
| 597 | + | ||
| 598 | + # Allocate caldb | ||
| 599 | + db = caldb("E") | ||
| 600 | + | ||
| 601 | + # Open calibration | ||
| 602 | + db.open("test") | ||
| 603 | + | ||
| 604 | + # Translate ROOT to CALDB information | ||
| 605 | + db.root2caldb(filename) | ||
| 606 | + | ||
| 0 | \ No newline at end of file | 607 | \ No newline at end of file |