diff --git a/scripts/cta_root2caldb.py b/scripts/cta_root2caldb.py index b757081..1d63fb5 100755 --- a/scripts/cta_root2caldb.py +++ b/scripts/cta_root2caldb.py @@ -2,11 +2,9 @@ # ========================================================================== # This script generates IRFs in the CALDB format of HEASARC using ROOT 2D # performance files. -# -# Requirements: -# - pyroot -# -# Copyright (C) 2011 Jurgen Knodlseder +# ------------------------------------------------------------------------- +# Copyright (C) 2011-2012 Juergen Knoedlseder +# ------------------------------------------------------------------------- # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -21,6 +19,10 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # +# ------------------------------------------------------------------------- +# Requirements: +# - gammalib +# - pyroot # ========================================================================== from ROOT import TFile, TH1F, TH2F from gammalib import * @@ -37,18 +39,27 @@ import os class caldb(): """ """ - def __init__(self, conf, path=""): + def __init__(self, inst, obsid, path=""): """ - Constructor. + CALDB constructor. The constructor will create one CALDB entry for + an observation. + + Parameters: + inst - Instrument string. + obsid - Observation ID string. + Keywords: + path - CALDB Root path. """ # Store root path self.path = path # Initialise some members self.cif = None + self.irf = None self.ea = None self.psf = None self.edisp = None + self.bgd = None self.hdu_cif = None # Store UTC date of job execution @@ -56,7 +67,8 @@ class caldb(): # Set CALDB information self.cal_tel = "CTA" - self.cal_inst = conf.upper() + self.cal_inst = inst.upper() + self.cal_obsid = obsid self.cal_det = "NONE" self.cal_flt = "NONE" self.cal_class = "BCF" @@ -84,6 +96,11 @@ class caldb(): self.edisp_doc = "???" self.edisp_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] self.edisp_desc = "CTA energy dispersion" + # + self.bgd_name = "BGD" + self.bgd_doc = "???" + self.bgd_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] + self.bgd_desc = "CTA background" # Create directory structure self.make_dirs() @@ -93,6 +110,7 @@ class caldb(): def __del__(self): """ + Destructor. Close all FITS files. """ # Close calibration self.close() @@ -102,32 +120,50 @@ class caldb(): def make_dirs(self): """ + Generate CALDB directory structure for one observation identifier. + The structure is given by + + data///bcf/ + + where is "cta" and is the instrument specified in the + CALDB constructor (the instrument may be used for different array + configurations). + + Parameters: + None + Keywords: + None """ # Set base path self.base_dir = "data" self.base_dir += "/"+self.cal_tel.lower() self.base_dir += "/"+self.cal_inst.lower() - # Set directory names self.bcf_dir = self.base_dir+"/bcf" - self.ea_dir = self.bcf_dir+"/ea" - self.psf_dir = self.bcf_dir+"/psf" - self.edisp_dir = self.bcf_dir+"/edisp" + self.obs_dir = self.bcf_dir+"/"+self.cal_obsid + self.ea_dir = self.obs_dir + self.psf_dir = self.obs_dir + self.edisp_dir = self.obs_dir + self.bgd_dir = self.obs_dir - # Add path (if required) + # Optionally prefix path if len(self.path) > 0: self.base_path = self.path+"/"+self.base_dir self.bcf_path = self.path+"/"+self.bcf_dir + self.obs_path = self.path+"/"+self.obs_dir self.ea_path = self.path+"/"+self.ea_dir self.psf_path = self.path+"/"+self.psf_dir self.edisp_path = self.path+"/"+self.edisp_dir + self.bgd_path = self.path+"/"+self.bgd_dir else: self.base_path = self.base_dir self.bcf_path = self.bcf_dir + self.obs_path = self.obs_dir self.ea_path = self.ea_dir self.psf_path = self.psf_dir self.edisp_path = self.edisp_dir + self.bgd_path = self.bgd_dir # Create directory structure if not os.path.isdir(self.ea_path): @@ -136,18 +172,36 @@ class caldb(): os.makedirs(self.psf_path) if not os.path.isdir(self.edisp_path): os.makedirs(self.edisp_path) + if not os.path.isdir(self.bgd_path): + os.makedirs(self.bgd_path) # Return return - def open(self, name, clobber=True): + def open(self, version, split=False, clobber=True): """ - Open existing or create new calibration. + Open existing or create new calibration. The actual version will put + all calibrations in the same file, although each part of the response + function will have its own logical name. We can thus easily modify + the script to put each calibration information in a separate file. + + Parameters: + version - Filename version + Keywords: + split - Split IRF over several files? + clobber - Overwrite existing files? """ # Set calibrate file names - self.ea_file = "ea_"+name+".fits" - self.psf_file = "psf_"+name+".fits" - self.edisp_file = "edisp_"+name+".fits" + if split: + self.ea_file = "ea_"+version+".fits" + self.psf_file = "psf_"+version+".fits" + self.edisp_file = "edisp_"+version+".fits" + self.bgd_file = "bgd_"+version+".fits" + else: + self.ea_file = "irf_"+version+".fits" + self.psf_file = "irf_"+version+".fits" + self.edisp_file = "irf_"+version+".fits" + self.bgd_file = "irf_"+version+".fits" # Open calibration database index self.cif = GFits(self.base_path+"/caldb.indx", True) @@ -159,17 +213,42 @@ class caldb(): self.create_cif() self.hdu_cif = self.cif.table("CIF") - # Open calibration files - self.open_ea(self.ea_path+"/"+self.ea_file) - self.open_psf(self.psf_path+"/"+self.psf_file) - self.open_edisp(self.edisp_path+"/"+self.edisp_file) + # Set filenames + ea_filename = self.ea_path+"/"+self.ea_file + psf_filename = self.psf_path+"/"+self.psf_file + edisp_filename = self.edisp_path+"/"+self.edisp_file + bgd_filename = self.bgd_path+"/"+self.bgd_file + + # Open files + if split: + self.ea = GFits(ea_filename, True) + self.psf = GFits(psf_filename, True) + self.edisp = GFits(edisp_filename, True) + self.bgd = GFits(bgd_filename, True) + else: + self.irf = GFits(ea_filename, True) + self.ea = self.irf + self.psf = self.irf + self.edisp = self.irf + self.bgd = self.irf + + # Open HDUs + self.open_ea() + self.open_psf() + self.open_edisp() + self.open_bgd() # Return return def close(self): """ - Close calibration. + Close calibration FITS files. + + Parameters: + None + Keywords: + None """ # Add information to CIF. We do this now as before we did not # necessarily have all the information at hand (in particular @@ -181,31 +260,57 @@ class caldb(): self.cif.save(True) self.cif.close() self.cif = None - - # Close effective area - if self.ea != None: - self.ea.save(True) - self.ea.close() - self.ea = None - - # Close point spread function - if self.psf != None: - self.psf.save(True) - self.psf.close() - self.psf = None - - # Close energy dispersion - if self.edisp != None: - self.edisp.save(True) - self.edisp.close() + + # If IRF file exists then close it now. All file pointers will + # be set to None + if self.irf != None: + self.irf.save(True) + self.irf.close() + self.irf = None + self.ea = None + self.psf = None self.edisp = None + self.bgd = None + + # ... otherwise we have split files, so we have to close them + # all + else: + + # Close effective area file + if self.ea != None: + self.ea.save(True) + self.ea.close() + self.ea = None + + # Close point spread function file + if self.psf != None: + self.psf.save(True) + self.psf.close() + self.psf = None + + # Close energy dispersion file + if self.edisp != None: + self.edisp.save(True) + self.edisp.close() + self.edisp = None + + # Close background file + if self.bgd != None: + self.bgd.save(True) + self.bgd.close() + self.bgd = None # Return return def create_cif(self): """ - Create CIF extension in FITS file. + Create Calibration Index File (CIF) extension in CIF FITS file. + + Parameters: + None + Keywords: + None """ # Create binary table table = GFitsBinTable() @@ -248,6 +353,11 @@ class caldb(): def add_cif_info(self): """ Add information to CIF extension. + + Parameters: + None + Keywords: + None """ # Append 3 rows to CIF extension self.hdu_cif.append_rows(3) @@ -310,20 +420,35 @@ class caldb(): self.hdu_cif["CAL_CBD"][row, n] = "NONE" else: self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n] + + # Add background information + row = self.hdu_cif.nrows()-1 + self.hdu_cif["CAL_DIR"][row] = self.bgd_dir + self.hdu_cif["CAL_FILE"][row] = self.bgd_file + self.hdu_cif["CAL_CNAM"][row] = self.bgd_name + self.hdu_cif["CAL_DESC"][row] = self.bgd_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.bgd_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n] # Return return - def open_ea(self, filename): + def open_ea(self): """ Open effective area extension. If no extension was found then create one. We do not yet set any data as we don't know the dimensions of the table yet. - """ - # Open FITS file - self.ea = GFits(filename, True) - # Get extensions. If they do not exist then create them now + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. try: self.hdu_ea = self.ea.table("EFFECTIVE AREA") except: @@ -346,16 +471,18 @@ class caldb(): # Return return - def open_psf(self, filename): + def open_psf(self): """ Open point spread function extension. If no extension was found then create one. We do not yet set any data as we don't know the dimensions of the table yet. - """ - # Open FITS file - self.psf = GFits(filename, True) - # Get extensions. If they do not exist then create them now + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. try: self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") except: @@ -378,16 +505,18 @@ class caldb(): # Return return - def open_edisp(self, filename): + def open_edisp(self): """ Open energy dispersion extension. If no extension was found then create one. We do not yet set any data as we don't know the dimensions of the table yet. - """ - # Open FITS file - self.edisp = GFits(filename, True) - # Get extensions. If they do not exist then create them now + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. try: self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") except: @@ -410,9 +539,50 @@ class caldb(): # Return return + def open_bgd(self): + """ + Open background extension. If no extension was found + then create one. We do not yet set any data as we don't + know the dimensions of the table yet. + + Parameters: + None + Keywords: + None + """ + # Get extension. If it does not exist then create it now. + try: + self.hdu_bgd = self.bgd.table("BACKGROUND") + except: + # Create binary table + table = GFitsBinTable() + + # Set extension name + table.extname("BACKGROUND") + + # Set OGIP keywords + self.set_ogip_keywords(table, self.bgd_doc, \ + ["RESPONSE", self.bgd_name]) + + # Append table + self.bgd.append(table) + + # Get extension + self.hdu_bgd = self.bgd.table("BACKGROUND") + + # Return + return + def set_ogip_keywords(self, hdu, hdudoc, hduclas): """ Set standard OGIP keywords for extension. + + Parameters: + hdu - Header Data Unit. + hdudoc - Documentation reference string + hduclas - Array of HDUCLAS fields + Keywords: + None """ # Set keywords hdu.card("ORIGIN", "IRAP", "Name of organization making this file") @@ -433,12 +603,17 @@ class caldb(): def make_2D(self, array, hdu, name, unit, scale=1.0): """ Make 2D matrix as function of energy and offset angle from a - ROOT 2D histogram. + ROOT 2D histogram. If the HDU has already the energy and + offset angle columns, this method will simply add another data + column. + If name==None, the method will not append any data column. Parameters: - array - ROOT 2D histgram. + array - ROOT 2D histogram. hdu - FITS HDU. - Keyword: + name - Data column name. + unit - Data unit. + Keywords: scale - Scaling factor for histogram values. """ # Extract energy and offset angle vectors @@ -485,7 +660,7 @@ class caldb(): hdu["THETA_HI"][0,ioff] = o_hi # # "NAME" - if not hdu.hascolumn(name): + if name != None and not hdu.hascolumn(name): hdu.append_column(GFitsTableFloatCol(name, 1, neng*noffset)) hdu[name].unit(unit) hdu[name].dim([neng, noffset]) @@ -515,44 +690,217 @@ class caldb(): Parameters: filename - ROOT 2D performance filename. """ - # Open performance file + # Open ROOT performance file file = TFile(filename) - # Allocate ROOT 2D array - array = TH2F() + # Create effective area + self.root2ea(file) + + # Create point spread function + self.root2psf(file) + + # Create energy dispersion + self.root2edisp(file) + + # Create background + self.root2bgd(file) + + # Return + return - # Make 2D effective 80% area + def root2ea(self, file): + """ + Translate ROOT to CALDB effective area extension. The following ROOT + histograms are used: + + EffectiveAreaEtrue_offaxis -> EFFAREA + EffectiveArea_offaxis -> EFFAREA_RECO + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if effective area HDU has been opened if self.hdu_ea != None: - file.GetObject("EffectiveArea80_offaxis", array) - bounds = self.make_2D(array, self.hdu_ea, "EFFAREA", "m2", scale=1.0/0.8) + + # Allocate ROOT 2D arrays + etrue = TH2F() + ereco = TH2F() + file.GetObject("EffectiveAreaEtrue_offaxis", etrue) + file.GetObject("EffectiveArea_offaxis", ereco) + #print etrue.GetXaxis().GetBinLowEdge(1), etrue.GetXaxis().GetBinUpEdge(etrue.GetXaxis().GetNbins()) + #print ereco.GetXaxis().GetBinLowEdge(1), ereco.GetXaxis().GetBinUpEdge(ereco.GetXaxis().GetNbins()) + + # Rebin etrue histogram + etrue.RebinX(10) + neng = etrue.GetXaxis().GetNbins() + noffset = etrue.GetYaxis().GetNbins() + for ioff in range(noffset): + for ieng in range(neng): + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0 + etrue.SetBinContent(ieng+1,ioff+1,value) + + # Write boundaries (use Ereco boundaries) + bounds = self.make_2D(ereco, self.hdu_ea, None, "m2") for b in bounds: self.ea_bounds.append(b) self.set_cif_keywords(self.hdu_ea, self.ea_name, \ self.ea_bounds, self.ea_desc) - # Make 2D 68% point spread function + # EFFAREA + self.make_2D(etrue, self.hdu_ea, "EFFAREA", "m2") + + # EFFAREA_RECO + self.make_2D(ereco, self.hdu_ea, "EFFAREA_RECO", "m2") + + # Return + return + + def root2psf(self, file): + """ + Translate ROOT to CALDB point spread function extension. The following + ROOT histograms are used: + + 1/(2*pi*SIGMA_1) -> SCALE + AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8) + 0.0 -> AMPL_2 + 0.0 -> SIGMA_2 + 0.0 -> AMPL_3 + 0.0 -> SIGMA_3 + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if point spread function HDU has been opened if self.hdu_psf != None: - file.GetObject("AngRes_offaxis", array) - bounds = self.make_2D(array, self.hdu_psf, "R68", "deg") + + # Allocate ROOT 2D array + r68 = TH2F() + + # Read 68% containment histogram + file.GetObject("AngRes_offaxis", r68) + neng = r68.GetXaxis().GetNbins() + noffset = r68.GetYaxis().GetNbins() + + # Compute scale histogram + scale = r68.Clone() + for ioff in range(noffset): + for ieng in range(neng): + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1) + if integral > 0.0: + value = 1.0 / integral + else: + value = 0.0 + scale.SetBinContent(ieng+1,ioff+1,value) + + # Set zero histogram + zero = r68.Clone() + for ioff in range(noffset): + for ieng in range(neng): + zero.SetBinContent(ieng+1,ioff+1,0.0) + + # Set boundaries + bounds = self.make_2D(r68, self.hdu_psf, None, "deg") for b in bounds: self.psf_bounds.append(b) self.set_cif_keywords(self.hdu_psf, self.psf_name, \ self.psf_bounds, self.psf_desc) - # Make 2D 80% point spread function - if self.hdu_psf != None: - file.GetObject("AngRes80_offaxis", array) - self.make_2D(array, self.hdu_psf, "R80", "deg") + # SCALE + self.make_2D(scale, self.hdu_psf, "SCALE", "") + + # SIGMA_1 + self.make_2D(r68, self.hdu_psf, "SIGMA_1", "deg") - # Make 2D energy resolution + # AMPL_2 + self.make_2D(zero, self.hdu_psf, "AMPL_2", "") + + # SIGMA_2 + self.make_2D(zero, self.hdu_psf, "SIGMA_2", "deg") + + # AMPL_3 + self.make_2D(zero, self.hdu_psf, "AMPL_3", "") + + # SIGMA_3 + self.make_2D(zero, self.hdu_psf, "SIGMA_3", "deg") + + # Return + return + + def root2edisp(self, file): + """ + Translate ROOT to CALDB energy dispersion extension. The following ROOT + histograms are used: + + ERes_offaxis -> ERES + Ebias_offaxis -> EBIAS + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if energy dispersion HDU has been opened if self.hdu_edisp != None: - file.GetObject("ERes", array) - bounds = self.make_2D(array, self.hdu_edisp, "ERES68", "") + + # Allocate ROOT 2D array + eres = TH2F() + ebias = TH2F() + file.GetObject("ERes_offaxis", eres) + file.GetObject("Ebias_offaxis", ebias) + + # Set boundaries + bounds = self.make_2D(eres, self.hdu_edisp, None, "deg") for b in bounds: self.edisp_bounds.append(b) self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \ self.edisp_bounds, self.edisp_desc) - + + # ERES + self.make_2D(eres, self.hdu_edisp, "ERES", "Ereco/Etrue") + + # EBIAS + self.make_2D(ebias, self.hdu_edisp, "EBIAS", "Ereco/Etrue") + + # Return + return + + def root2bgd(self, file): + """ + Translate ROOT to CALDB background extension. The following ROOT + histograms are used: + + BGRatePerSqDeg_offaxis -> BGD + BGRatePerSqDeg_offaxis -> BGD_RECO + + Parameters: + file - ROOT file. + Keywords: + None + """ + # Continue only if background HDU has been opened + if self.hdu_bgd != None: + + # Allocate ROOT 2D array + array = TH2F() + file.GetObject("BGRatePerSqDeg_offaxis", array) + + # Set boundaries + bounds = self.make_2D(array, self.hdu_bgd, None, "deg") + for b in bounds: + self.bgd_bounds.append(b) + self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \ + self.bgd_bounds, self.bgd_desc) + + # BGD + self.make_2D(array, self.hdu_bgd, "BGD", "") + + # BGD_RECO + self.make_2D(array, self.hdu_bgd, "BGD_RECO", "") + # Return return @@ -589,22 +937,41 @@ class caldb(): #==========================# if __name__ == '__main__': """ - This script translates a ROOT 2D performance file in a CALDB compliant - CTA response. + This script translates a CTA ROOT 2D performance files in a CALDB compliant + CTA response. The script will create the response for 3 observations to + illustrate the run-wise file structure. + + The CTA ROOT 2D performance files contain the following histograms: + - DiffSens_offaxis (2D) - Differential sensitivity + - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?) + - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment + - AngRes_offaxis (2D) - Angular resolution 68% containment + - AngRes80_offaxis (2D) - Angular resolution 80% containment + - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree + - BGRate_offaxis (2D) - Background rate + - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins) + - MigMatrix_offaxis (3D) + - EestOverEtrue_offaxis (3D) + - ERes_offaxis (2D) - Energy Resolution + - Ebias_offaxis (2D) - Energy bias """ # Set ROOT filename - #path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxis" - #irf = "SubarrayE_offaxis.root" - path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_Nov2011" - irf = "SubarrayE_IFAE_50hours_20111121_offaxis.root" + path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_May2012" + irf = "SubarrayE_IFAE_50hours_20120510_offaxis.root" filename = path+"/"+irf - # Allocate caldb - db = caldb("E") + # Set observation identifiers + obsids = ["000001", "000002", "000003"] + + # Generate CALDB entries, one for each observation (or run) + for obsid in obsids: + + # Allocate caldb + db = caldb("E", obsid) - # Open calibration - db.open("test") + # Open calibration + db.open("test") - # Translate ROOT to CALDB information - db.root2caldb(filename) + # Translate ROOT to CALDB information + db.root2caldb(filename) \ No newline at end of file