From 419af8b07af16e379d94b59b48cbbc684ce71fb3 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 26 Jun 2016 22:57:41 +0200 Subject: [PATCH] Move some obsolete cscripts in the test/dev folder --- cscripts/cta_kb_bkgrates.py | 112 ---------------------------------------------------------------------------------------------------------------- cscripts/cta_kb_root2caldb.py | 546 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ cscripts/cta_make_nodes.py | 115 ------------------------------------------------------------------------------------------------------------------- cscripts/cta_root2irf.py | 291 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- test/dev/cta_kb_bkgrates.py | 112 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/dev/cta_kb_root2caldb.py | 546 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/dev/cta_make_nodes.py | 115 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/dev/cta_root2irf.py | 291 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 1064 insertions(+), 1064 deletions(-) delete mode 100755 cscripts/cta_kb_bkgrates.py delete mode 100755 cscripts/cta_kb_root2caldb.py delete mode 100755 cscripts/cta_make_nodes.py delete mode 100755 cscripts/cta_root2irf.py create mode 100755 test/dev/cta_kb_bkgrates.py create mode 100755 test/dev/cta_kb_root2caldb.py create mode 100755 test/dev/cta_make_nodes.py create mode 100755 test/dev/cta_root2irf.py diff --git a/cscripts/cta_kb_bkgrates.py b/cscripts/cta_kb_bkgrates.py deleted file mode 100755 index 559e73d..0000000 --- a/cscripts/cta_kb_bkgrates.py +++ /dev/null @@ -1,112 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script converts the background rates in Konrad's performance files -# into a column separated file function table. The table is automatically -# named "bkg_XXX.txt" where "XXX" is for example "kb_E_50h_v3". -# -# Copyright (C) 2011-2016 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 -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -import os -import glob -import math -import sys - - -# ================================= # -# Generate file function from model # -# ================================= # -def make_file_function(irfname, filename): - """ - Generate file function from model. - - Parameters: - irfname - IRF filename - filename - File function filename - """ - # Open performance file and file function - irf = open(irfname, "r") - f = open(filename, "w") - - # Read data - for row in irf: - - # Skip header and trailer - if row.find("log(E)") != -1: - continue - if row.find("--------") != -1: - break - - # Extract data - fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",") - logE = float(fields[0]) - #aeff = float(fields[1]) - #r68 = float(fields[2]) - r80 = float(fields[3]) - #eres = float(fields[4]) - bgd = float(fields[5]) - #sens = float(fields[6]) - - # Compute energy (in MeV) - energy = math.pow(10.0, logE)*1.0e6 - emin = math.pow(10.0, logE-0.1)*1.0e6 - emax = math.pow(10.0, logE+0.1)*1.0e6 - ewidth = emax-emin - - # Compute solid angle of r80 region in sr - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) - - # Compute background rate per steradian and MeV - if omega > 0.0: - bkg_rate = bgd/omega/ewidth - else: - bkg_rate = 0.0 - - # Write - f.write(str(energy)+" "+str(bkg_rate)+"\n") - - # Close files - irf.close() - f.close() - - # Return - return - - -# ================ # -# Main entry point # -# ================ # -if __name__ == '__main__': - - # Set parameters - caldb = "/usr/local/gamma/share/caldb/data/cta/kb/" - - # Get IRF from command line argument - if (len(sys.argv) > 1): - irfs = [caldb+"/"+sys.argv[1]] - else: - irfs = glob.glob(caldb+"/*.dat") - irfs.sort() - - # Loop over all IRFs - for irfname in irfs: - - # Build background filename - fname = "bkg_"+os.path.basename(irf).strip(".dat")+".txt" - print(fname) - - # Make background file - make_file_function(irf, fname) diff --git a/cscripts/cta_kb_root2caldb.py b/cscripts/cta_kb_root2caldb.py deleted file mode 100755 index fca22eb..0000000 --- a/cscripts/cta_kb_root2caldb.py +++ /dev/null @@ -1,546 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script generates CALDB compliant performance tables from MC ROOT -# files. -# -# Note that the effective area given in the ASCII file is the full effective -# area (i.e. the area before applying any theta cut). This is required by -# GammaLib as the maximum likelihood analysis does not imply any theta cut. -# -# Copyright (C) 2014-2016 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 -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -#from ROOT import TFile, TH1F, TH2F -from ROOT import TFile, TH1F -from datetime import datetime -import gammalib -import math -import sys -import glob -import os - - -# =========== # -# CalDB class # -# =========== # -class caldb(): - """ - """ - def __init__(self, inst, path=""): - """ - CALDB constructor. The constructor will create one CALDB entry for - an observation. - - Parameters: - inst - Instrument 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 - self.utc = datetime.utcnow().isoformat()[:19] - - # Set CALDB information - self.cal_tel = "CTA" - self.cal_inst = inst.upper() - self.cal_det = "NONE" - self.cal_flt = "NONE" - self.cal_class = "BCF" - self.cal_type = "DATA" - self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ... - self.cal_date = "14/01/01" - self.val_date = "2014-01-01" - self.val_time = "00:00:00" - self.ref_time = 51544.0 - self.cal_name = "NAME()" - self.cal_version = "VERSION(Prod1)" - self.cal_cut = "CLASS(BEST)" - self.cal_analysis = "ANALYSIS(MPIK)" - # - self.ea_name = "EFF_AREA" - self.ea_doc = "???" - self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.ea_desc = "CTA effective area" - # - self.psf_name = "RPSF" - self.psf_doc = "???" - self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.psf_desc = "CTA point spread function" - # - self.edisp_name = "EDISP" - self.edisp_doc = "???" - self.edisp_bounds = [self.cal_name, 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_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.bgd_desc = "CTA background" - - # Create directory structure - self.make_dirs() - - # Return - return - - def __del__(self): - """ - Destructor. Close all FITS files. - """ - # Close calibration - self.close() - - # Return - return - - 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" - - # Optionally prefix path - if len(self.path) > 0: - self.base_path = self.path+"/"+self.base_dir - self.bcf_path = self.path+"/"+self.bcf_dir - else: - self.base_path = self.base_dir - self.bcf_path = self.bcf_dir - - # Create directory structure - if not os.path.isdir(self.bcf_path): - os.makedirs(self.bcf_path) - - # Return - return - - def add(self, rspname, split=False, clobber=True): - """ - Add 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: - rspname - Response name - Keywords: - split - Split IRF over several files? - clobber - Overwrite existing files? - """ - # Set calibrate file names - if split: - self.ea_file = "ea_"+rspname+".dat" - self.psf_file = "psf_"+rspname+".dat" - self.edisp_file = "edisp_"+rspname+".dat" - self.bgd_file = "bgd_"+rspname+".dat" - else: - self.ea_file = "irf_"+rspname+".dat" - self.psf_file = "irf_"+rspname+".dat" - self.edisp_file = "irf_"+rspname+".dat" - self.bgd_file = "irf_"+rspname+".dat" - - # Open calibration database index - if self.cif == None: - self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True) - - # If file has no CIF extension than create it now - try: - self.hdu_cif = self.cif.table("CIF") - except: - self.create_cif() - self.hdu_cif = self.cif.table("CIF") - - # Set response name - self.cal_name = "NAME("+rspname+")" - - # Return - return - - def close(self): - """ - 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 - # about the boundaries) - #self.add_cif_info() - - # Close CIF - if self.cif != None: - self.cif.save(True) - self.cif.close() - self.cif = None - - # Free - self.cif = None - - # Return - return - - def create_cif(self): - """ - Create Calibration Index File (CIF) extension in CIF FITS file. - - Parameters: - None - Keywords: - None - """ - # Create binary table - table = gammalib.GFitsBinTable() - - # Set boundary - - # Attach columns. Reference: CAL/GEN/92-008 - table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10)) - table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10)) - table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20)) - table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10)) - table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20)) - table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70)) - table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40)) - table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3)) - table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4)) - table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20)) - table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9)) - table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0)) - table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10)) - table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8)) - table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0)) - table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0)) - table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8)) - table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70)) - - # Set keywords. Reference: CAL/GEN/92-008 - table.extname("CIF") - table.card("CIFVERSN", "1992a", "Version of CIF format") - - # Attach table to FITS file. Note that at this moment the - # FITS table gets connected to the FITS file. Yet since - # nothing was yet written to the FITS file, we cannot read - # anything from it. - self.cif.append(table) - - # Return - return - - def add_cif_info(self): - """ - Add information to CIF extension. - - Parameters: - None - Keywords: - None - """ - # Set boundaries - self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] - - # Append 4 rows to CIF extension - self.hdu_cif.append_rows(4) - - # Add generic information for these 4 rows - for i in range(4): - - # Set row number - row = i+self.hdu_cif.nrows()-4 - - # Set element - self.hdu_cif["TELESCOP"][row] = self.cal_tel - self.hdu_cif["INSTRUME"][row] = self.cal_inst - self.hdu_cif["DETNAM"][row] = self.cal_det - self.hdu_cif["FILTER"][row] = self.cal_flt - self.hdu_cif["CAL_DEV"][row] = "ONLINE" - self.hdu_cif["CAL_CLAS"][row] = self.cal_class - self.hdu_cif["CAL_DTYP"][row] = self.cal_type - self.hdu_cif["CAL_VSD"][row] = self.val_date - self.hdu_cif["CAL_VST"][row] = self.val_time - self.hdu_cif["REF_TIME"][row] = self.ref_time - self.hdu_cif["CAL_QUAL"][row] = self.cal_qual - self.hdu_cif["CAL_DATE"][row] = self.cal_date - - # Add effective area information - row = self.hdu_cif.nrows()-4 - self.hdu_cif["CAL_DIR"][row] = self.bcf_dir - self.hdu_cif["CAL_FILE"][row] = self.ea_file - self.hdu_cif["CAL_CNAM"][row] = self.ea_name - self.hdu_cif["CAL_DESC"][row] = self.ea_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.ea_bounds): - self.hdu_cif["CAL_CBD"][row, n] = "NONE" - else: - self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n] - - # Add point spread function information - row = self.hdu_cif.nrows()-3 - self.hdu_cif["CAL_DIR"][row] = self.bcf_dir - self.hdu_cif["CAL_FILE"][row] = self.psf_file - self.hdu_cif["CAL_CNAM"][row] = self.psf_name - self.hdu_cif["CAL_DESC"][row] = self.psf_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.psf_bounds): - self.hdu_cif["CAL_CBD"][row, n] = "NONE" - else: - self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n] - - # Add energy dispersion information - row = self.hdu_cif.nrows()-2 - self.hdu_cif["CAL_DIR"][row] = self.bcf_dir - self.hdu_cif["CAL_FILE"][row] = self.edisp_file - self.hdu_cif["CAL_CNAM"][row] = self.edisp_name - self.hdu_cif["CAL_DESC"][row] = self.edisp_desc - self.hdu_cif["CAL_XNO"][row] = 1 - for n in range(9): - if n >= len(self.edisp_bounds): - 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.bcf_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 root2caldb(self, filename): - """ - Generate files. - - Parameters: - filename - Performance file name - """ - # Open ROOT performance file - rootfile = TFile(filename) - - # Set filenames - ea_filename = self.bcf_path+"/"+self.ea_file - psf_filename = self.bcf_path+"/"+self.psf_file - edisp_filename = self.bcf_path+"/"+self.edisp_file - bgd_filename = self.bcf_path+"/"+self.bgd_file - - # Open ASCII performance tables - if self.ea_file == self.psf_file and \ - self.ea_file == self.edisp_file and \ - self.ea_file == self.bgd_file : - split = False - fprf = open(ea_filename, "w") - files = [fprf] - else: - split = True - fea = open(ea_filename, "w") - fpsf = open(psf_filename, "w") - fedisp = open(edisp_filename, "w") - fbgd = open(bgd_filename, "w") - files = [fea, fpsf, fedisp, fbgd] - - # Write header - for file in files: - file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") - - # Get histograms. Konrad's root files do not have an EffectiveArea80 - # histogram, but his effective area is for 80% containment radius - # anyways, so we simply read this histogram into aeff80. - sens = TH1F() - bgrate = TH1F() - bgratesqdeg = TH1F() - aeff = TH1F() - aeff80 = TH1F() - angres = TH1F() - angres80 = TH1F() - eres = TH1F() - rootfile.GetObject("DiffSens", sens) - rootfile.GetObject("BGRate", bgrate) - rootfile.GetObject("BGRatePerSqDeg", bgratesqdeg) - rootfile.GetObject("EffectiveArea", aeff) - try: - rootfile.GetObject("EffectiveArea80", aeff80) - except: - rootfile.GetObject("EffectiveArea", aeff80) - rootfile.GetObject("AngRes", angres) - rootfile.GetObject("AngRes80", angres80) - rootfile.GetObject("ERes", eres) - - # Setup vectors - energies = bgrate.GetXaxis() - nbins_eng = energies.GetNbins() - for i in range(nbins_eng): - - # Get logE and background rate per squared degree - logE = energies.GetBinCenter(i+1) - bgd = bgrate.GetBinContent(i+1,1) - bgdsq = bgratesqdeg.GetBinContent(i+1,1) - area = aeff.GetBinContent(i+1,1) - area80 = aeff80.GetBinContent(i+1,1) - r68 = angres.GetBinContent(i+1,1) - r80 = angres80.GetBinContent(i+1,1) - dE = eres.GetBinContent(i+1,1) - diffs = sens.GetBinContent(i+1,1) - - # Skip any NaNs - if math.isnan(area80): - continue - - # Compute energy (in MeV) - energy = math.pow(10.0, logE)*1.0e6 - emin = math.pow(10.0, logE-0.1)*1.0e6 - emax = math.pow(10.0, logE+0.1)*1.0e6 - ewidth = emax-emin - - # Convert into background rate per steradian and MeV - bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth - - # Compute control background rate (this only works for KB files - # as they are for 80% containment radius). But we don't use this - # anyways, it's just for control on display ... - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) - if omega > 0.0: - bkg_rate2 = bgd / omega / ewidth - else: - bkg_rate2 = 0.0 - - # Compute full effective area by dividing area80/0.80 - area = area80/0.80 - - # Write results in file - line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ - (logE, area, r68, r80, dE, bgd, diffs) - for file in files: - file.write(line+"\n") - #fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") - - # Write trailer - for file in files: - file.write("---------------------------------------------\n") - file.write("Notes\n") - file.write(" 1) log(E) = log10(E/TeV) - bin centre\n") - file.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") - file.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") - file.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") - file.write(" 5) Fractional Energy Resolution (rms)\n") - file.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") - file.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") - file.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") - file.write(" systematics and statistics and at least 10 photons.\n") - - # Close files - for file in files: - file.close() - - # Close ROOT performance file - rootfile.Close() - - # Add information to CIF. We do this now as before we did not - # necessarily have all the information at hand (in particular - # about the boundaries) - self.add_cif_info() - - # Return - return - - -#==========================# -# Main routine entry point # -#==========================# -if __name__ == '__main__': - """ - """ - # Get ROOT filename from command line argument - #if (len(sys.argv) > 1): - # filename = sys.argv[1] - #else: - db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_20deg_root" - #db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" - - # Get list of all performance files - files = glob.glob(db+"/*.root") - files.sort() - - # Allocate caldb - irf = caldb("kb") - - # Loop over all files - for f in files: - - # Get filename - head, tail = os.path.split(f) - - # Build name - name = tail.strip(".root") - if name.find("kb_") != -1: - name = name.replace("_20deg", "") - elif name.find("IFAE_") != -1: - name = name.replace("Subarray", "") - name = name.replace("IFAE_", "") - name = "IFAE_"+name - index = name.find("hours") - if index != -1: - name = name[0:index+1] - print(tail) - - # Open calibration file - irf.add(name) - - # Translate ROOT to CALDB information - irf.root2caldb(f) diff --git a/cscripts/cta_make_nodes.py b/cscripts/cta_make_nodes.py deleted file mode 100755 index 9d8a3dd..0000000 --- a/cscripts/cta_make_nodes.py +++ /dev/null @@ -1,115 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This Python script creates the node section of the NodeFunction using -# logarithmically spaced energy bins. The intensity scale is set to the -# HESS Crab intensity (assuming a power law). -# -# Copyright (C) 2012-2016 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 -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -import sys -import math - - -# ============== # -# Write one node # -# ============== # -def write_node(f, energy, scale): - """ - Writes one node to XML file. - """ - # Convert to strings - str_energy = str(energy) - str_scale = str(scale) - - # Write start tag - f.write(' \n') - - # Write energy - f.write(' \n') - - # Write intensity - f.write(' \n') - - # Write end tag - f.write(' \n') - - # Return - return - - -# ============ # -# Create nodes # -# ============ # -def create_nodes(emin, emax, enumbins): - """ - Create nodes (energies in TeV). - """ - # Open file - f = open("nodes.xml", "w") - - # Set node boundaries - elogmin = math.log10(float(emin)) - elogmax = math.log10(float(emax)) - elogbin = (elogmax - elogmin)/(float(enumbins)-1.0) - - # Fill arrays - for i in range(int(enumbins)): - - # Compute energy - energy = math.pow(10.0, i*elogbin+elogmin) - - # Compute scale (HESS Crab spectrum) - scale = 3.45e-17 * math.pow(energy, -2.63) - - # Write node - write_node(f, energy, scale) - - # Debug - #sys.stdout.write(energy) - #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n") - - # Close file - f.close() - - # Return - return - - -# ======================= # -# Main script entry point # -# ======================= # -if __name__ == '__main__': - - # Check command line - usage = "Usage: cta_make_nodes emin emax enumbins" - if len(sys.argv) < 3: - sys.stdout.write(usage+"\n") - sys.exit() - - # Extract parameters - emin = sys.argv[1] - emax = sys.argv[2] - enumbins = sys.argv[3] - - # Create nodes - create_nodes(emin, emax, enumbins) diff --git a/cscripts/cta_root2irf.py b/cscripts/cta_root2irf.py deleted file mode 100755 index e0f6654..0000000 --- a/cscripts/cta_root2irf.py +++ /dev/null @@ -1,291 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script generates ASCII performance files and file functions for -# background modelling from MC ROOT files. On input it takes the filename -# of the ROOT performance file while on output it produced 2 ASCII files. -# The file "performance.dat" contains the performance data in the usual -# ASCII format, the file "background.txt" contains a file function that can -# be used within GammaLib for spectral background modelling. -# -# Note that the effective area given in the ASCII file is the full effective -# area (i.e. the area before applying any theta cut). This is required by -# GammaLib as the maximum likelihood analysis does not imply any theta cut. -# -# Copyright (C) 2011-2014 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 -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -from ROOT import TFile, TH1F, TH2F -import math -import sys -import glob -import os - - -# ================================== # -# Compute Eq. (17) of Li & Ma (1983) # -# ================================== # -def sigma_lima(Non, Noff, alpha=0.2): - """ - Compute Eq. (17) of Li & Ma (1983). - - Parameters: - Non - Number of on counts - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure - """ - # Compute sensitivity - alpha1 = alpha + 1.0 - sum = Non + Noff - arg1 = Non / sum - arg2 = Noff / sum - term1 = Non * math.log((alpha1/alpha)*arg1) - term2 = Noff * math.log(alpha1*arg2) - sigma = math.sqrt(2.0 * (term1 + term2)) - - # Return sensitivity - return sigma - - -# ======================================== # -# Solve Eq. (17) of Li & Ma (1983) for Non # -# ======================================== # -def Non_lima(sigma, Noff, alpha=0.2): - """ - Solve Eq. (17) of Li & Ma (1983) for Non - - Parameters: - sigma - Required significance - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure - """ - # Set initial guess for Non - S = 25.0 - Non = S + alpha*Noff - - # Compute initial sigma estimate - s = sigma_lima(Non, Noff, alpha=alpha) - delta = s - sigma - - # Iterate - while abs(delta) > 1.0e-6: - S *= sigma/s - Non = S + alpha*Noff - s = sigma_lima(Non, Noff, alpha=alpha) - delta = s - sigma - - # Return - return Non - - -# =============================================== # -# Compute effective area assuming a Crab spectrum # -# =============================================== # -def get_aeff(bgrate, emin, emax, \ - duration=180000.0, alpha=0.2, sigma=5.0, norm=2.79e-7, index=-2.57): - """ - Computes effective area assuming a Crab spectrum. This does not - include any hadron cut efficiencies, and thus is not very useful in a first - place ... - - Parameters: - bgrate - Background rate (TeV) - emin - Minimum energy (TeV) - emax - Maximum energy (TeV) - Keywords: - duration - Observing time in seconds (default: 50h) - alpha - Non/Noff ratio (default: 0.2) - sigma - Significance (default: 5.0) - norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7) - index - Crab spectral index (default: -2.57) - """ - # Compute number of source events - B = bgrate*duration - Noff = B/alpha - Non = Non_lima(sigma, Noff) - S = (Non - B) - - # Compute Crab flux in specified energy band - inx = index+1.0 - cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx)) - - # Compute effective area - area = S / duration / cu - - # Return effective area - return area - - -# ============== # -# Generate files # -# ============== # -def make_files(filename, prfname, bkgname): - """ - Generate files. - - Parameters: - filename - Performance file name - prfname - Performance file name - bkgname - File function filename - """ - # Open performance file and file function - file = TFile(filename) - fprf = open(prfname, "w") - fbkg = open(bkgname, "w") - - # Write header - fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") - #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") - - # Get histograms. Konrad's root files do not have an EffectiveArea80 - # histogram, but his effective area is for 80% containment radius - # anyways, so we simply read this histogram into aeff80. - sens = TH1F() - bgrate = TH1F() - bgratesqdeg = TH1F() - aeff = TH1F() - aeff80 = TH1F() - angres = TH1F() - angres80 = TH1F() - eres = TH1F() - file.GetObject("DiffSens", sens) - file.GetObject("BGRate", bgrate) - file.GetObject("BGRatePerSqDeg", bgratesqdeg) - file.GetObject("EffectiveArea", aeff) - try: - file.GetObject("EffectiveArea80", aeff80) - except: - file.GetObject("EffectiveArea", aeff80) - file.GetObject("AngRes", angres) - file.GetObject("AngRes80", angres80) - file.GetObject("ERes", eres) - - # Setup vectors - energies = bgrate.GetXaxis() - nbins_eng = energies.GetNbins() - for i in range(nbins_eng): - - # Get logE and background rate per squared degree - logE = energies.GetBinCenter(i+1) - bgd = bgrate.GetBinContent(i+1,1) - bgdsq = bgratesqdeg.GetBinContent(i+1,1) - area = aeff.GetBinContent(i+1,1) - area80 = aeff80.GetBinContent(i+1,1) - r68 = angres.GetBinContent(i+1,1) - r80 = angres80.GetBinContent(i+1,1) - dE = eres.GetBinContent(i+1,1) - diffs = sens.GetBinContent(i+1,1) - - # Skip any NaNs - if math.isnan(area80): - continue - - # Compute energy (in MeV) - energy = math.pow(10.0, logE)*1.0e6 - emin = math.pow(10.0, logE-0.1)*1.0e6 - emax = math.pow(10.0, logE+0.1)*1.0e6 - ewidth = emax-emin - - # Convert into background rate per steradian and MeV - bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth - - # Compute control background rate (this only works for KB files - # as they are for 80% containment radius). But we don't use this - # anyways, it's just for control on display ... - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) - if omega > 0.0: - bkg_rate2 = bgd / omega / ewidth - else: - bkg_rate2 = 0.0 - - # Compute full effective area by dividing area80/0.80 - area = area80/0.80 - - # Write results in file - line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ - (logE, area, r68, r80, dE, bgd, diffs) - fprf.write(line+"\n") - fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") - - # Show performance file - #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \ - # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2) - - # Write trailer - fprf.write("---------------------------------------------\n") - fprf.write("Notes\n") - fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n") - fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") - fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") - fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") - fprf.write(" 5) Fractional Energy Resolution (rms)\n") - fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") - fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") - fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") - fprf.write(" systematics and statistics and at least 10 photons.\n") - - # Close files - fprf.close() - fbkg.close() - - # Return - return - - -#==========================# -# Main routine entry point # -#==========================# -if __name__ == '__main__': - """ - """ - # Get ROOT filename from command line argument - if (len(sys.argv) > 1): - filename = sys.argv[1] - else: - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010" - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root" - caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" - - # Get list of all performance files - files = glob.glob(caldb+"/*.root") - files.sort() - - # Loop over all files - for file in files: - - # Get filename - head, tail = os.path.split(file) - - # Build name - name = tail.strip(".root") - if name.find("kb_") != -1: - name = name.replace("_20deg", "") - elif name.find("IFAE_") != -1: - name = name.replace("Subarray", "") - name = name.replace("IFAE_", "") - name = "IFAE_"+name - i = name.find("hours") - if i != -1: - name = name[0:i+1] - - # Build output filenames - prfname = name + ".dat" - bkgname = "bkg_"+name+".txt" - #print prfname, bkgname - - # Create files - make_files(file, prfname, bkgname) diff --git a/test/dev/cta_kb_bkgrates.py b/test/dev/cta_kb_bkgrates.py new file mode 100755 index 0000000..559e73d --- /dev/null +++ b/test/dev/cta_kb_bkgrates.py @@ -0,0 +1,112 @@ +#! /usr/bin/env python +# ========================================================================== +# This script converts the background rates in Konrad's performance files +# into a column separated file function table. The table is automatically +# named "bkg_XXX.txt" where "XXX" is for example "kb_E_50h_v3". +# +# Copyright (C) 2011-2016 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 +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import os +import glob +import math +import sys + + +# ================================= # +# Generate file function from model # +# ================================= # +def make_file_function(irfname, filename): + """ + Generate file function from model. + + Parameters: + irfname - IRF filename + filename - File function filename + """ + # Open performance file and file function + irf = open(irfname, "r") + f = open(filename, "w") + + # Read data + for row in irf: + + # Skip header and trailer + if row.find("log(E)") != -1: + continue + if row.find("--------") != -1: + break + + # Extract data + fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",") + logE = float(fields[0]) + #aeff = float(fields[1]) + #r68 = float(fields[2]) + r80 = float(fields[3]) + #eres = float(fields[4]) + bgd = float(fields[5]) + #sens = float(fields[6]) + + # Compute energy (in MeV) + energy = math.pow(10.0, logE)*1.0e6 + emin = math.pow(10.0, logE-0.1)*1.0e6 + emax = math.pow(10.0, logE+0.1)*1.0e6 + ewidth = emax-emin + + # Compute solid angle of r80 region in sr + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) + + # Compute background rate per steradian and MeV + if omega > 0.0: + bkg_rate = bgd/omega/ewidth + else: + bkg_rate = 0.0 + + # Write + f.write(str(energy)+" "+str(bkg_rate)+"\n") + + # Close files + irf.close() + f.close() + + # Return + return + + +# ================ # +# Main entry point # +# ================ # +if __name__ == '__main__': + + # Set parameters + caldb = "/usr/local/gamma/share/caldb/data/cta/kb/" + + # Get IRF from command line argument + if (len(sys.argv) > 1): + irfs = [caldb+"/"+sys.argv[1]] + else: + irfs = glob.glob(caldb+"/*.dat") + irfs.sort() + + # Loop over all IRFs + for irfname in irfs: + + # Build background filename + fname = "bkg_"+os.path.basename(irf).strip(".dat")+".txt" + print(fname) + + # Make background file + make_file_function(irf, fname) diff --git a/test/dev/cta_kb_root2caldb.py b/test/dev/cta_kb_root2caldb.py new file mode 100755 index 0000000..fca22eb --- /dev/null +++ b/test/dev/cta_kb_root2caldb.py @@ -0,0 +1,546 @@ +#! /usr/bin/env python +# ========================================================================== +# This script generates CALDB compliant performance tables from MC ROOT +# files. +# +# Note that the effective area given in the ASCII file is the full effective +# area (i.e. the area before applying any theta cut). This is required by +# GammaLib as the maximum likelihood analysis does not imply any theta cut. +# +# Copyright (C) 2014-2016 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 +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +#from ROOT import TFile, TH1F, TH2F +from ROOT import TFile, TH1F +from datetime import datetime +import gammalib +import math +import sys +import glob +import os + + +# =========== # +# CalDB class # +# =========== # +class caldb(): + """ + """ + def __init__(self, inst, path=""): + """ + CALDB constructor. The constructor will create one CALDB entry for + an observation. + + Parameters: + inst - Instrument 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 + self.utc = datetime.utcnow().isoformat()[:19] + + # Set CALDB information + self.cal_tel = "CTA" + self.cal_inst = inst.upper() + self.cal_det = "NONE" + self.cal_flt = "NONE" + self.cal_class = "BCF" + self.cal_type = "DATA" + self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ... + self.cal_date = "14/01/01" + self.val_date = "2014-01-01" + self.val_time = "00:00:00" + self.ref_time = 51544.0 + self.cal_name = "NAME()" + self.cal_version = "VERSION(Prod1)" + self.cal_cut = "CLASS(BEST)" + self.cal_analysis = "ANALYSIS(MPIK)" + # + self.ea_name = "EFF_AREA" + self.ea_doc = "???" + self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.ea_desc = "CTA effective area" + # + self.psf_name = "RPSF" + self.psf_doc = "???" + self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.psf_desc = "CTA point spread function" + # + self.edisp_name = "EDISP" + self.edisp_doc = "???" + self.edisp_bounds = [self.cal_name, 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_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.bgd_desc = "CTA background" + + # Create directory structure + self.make_dirs() + + # Return + return + + def __del__(self): + """ + Destructor. Close all FITS files. + """ + # Close calibration + self.close() + + # Return + return + + 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" + + # Optionally prefix path + if len(self.path) > 0: + self.base_path = self.path+"/"+self.base_dir + self.bcf_path = self.path+"/"+self.bcf_dir + else: + self.base_path = self.base_dir + self.bcf_path = self.bcf_dir + + # Create directory structure + if not os.path.isdir(self.bcf_path): + os.makedirs(self.bcf_path) + + # Return + return + + def add(self, rspname, split=False, clobber=True): + """ + Add 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: + rspname - Response name + Keywords: + split - Split IRF over several files? + clobber - Overwrite existing files? + """ + # Set calibrate file names + if split: + self.ea_file = "ea_"+rspname+".dat" + self.psf_file = "psf_"+rspname+".dat" + self.edisp_file = "edisp_"+rspname+".dat" + self.bgd_file = "bgd_"+rspname+".dat" + else: + self.ea_file = "irf_"+rspname+".dat" + self.psf_file = "irf_"+rspname+".dat" + self.edisp_file = "irf_"+rspname+".dat" + self.bgd_file = "irf_"+rspname+".dat" + + # Open calibration database index + if self.cif == None: + self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True) + + # If file has no CIF extension than create it now + try: + self.hdu_cif = self.cif.table("CIF") + except: + self.create_cif() + self.hdu_cif = self.cif.table("CIF") + + # Set response name + self.cal_name = "NAME("+rspname+")" + + # Return + return + + def close(self): + """ + 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 + # about the boundaries) + #self.add_cif_info() + + # Close CIF + if self.cif != None: + self.cif.save(True) + self.cif.close() + self.cif = None + + # Free + self.cif = None + + # Return + return + + def create_cif(self): + """ + Create Calibration Index File (CIF) extension in CIF FITS file. + + Parameters: + None + Keywords: + None + """ + # Create binary table + table = gammalib.GFitsBinTable() + + # Set boundary + + # Attach columns. Reference: CAL/GEN/92-008 + table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10)) + table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10)) + table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20)) + table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10)) + table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20)) + table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70)) + table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40)) + table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3)) + table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4)) + table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20)) + table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9)) + table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0)) + table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10)) + table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8)) + table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0)) + table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0)) + table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8)) + table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70)) + + # Set keywords. Reference: CAL/GEN/92-008 + table.extname("CIF") + table.card("CIFVERSN", "1992a", "Version of CIF format") + + # Attach table to FITS file. Note that at this moment the + # FITS table gets connected to the FITS file. Yet since + # nothing was yet written to the FITS file, we cannot read + # anything from it. + self.cif.append(table) + + # Return + return + + def add_cif_info(self): + """ + Add information to CIF extension. + + Parameters: + None + Keywords: + None + """ + # Set boundaries + self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + + # Append 4 rows to CIF extension + self.hdu_cif.append_rows(4) + + # Add generic information for these 4 rows + for i in range(4): + + # Set row number + row = i+self.hdu_cif.nrows()-4 + + # Set element + self.hdu_cif["TELESCOP"][row] = self.cal_tel + self.hdu_cif["INSTRUME"][row] = self.cal_inst + self.hdu_cif["DETNAM"][row] = self.cal_det + self.hdu_cif["FILTER"][row] = self.cal_flt + self.hdu_cif["CAL_DEV"][row] = "ONLINE" + self.hdu_cif["CAL_CLAS"][row] = self.cal_class + self.hdu_cif["CAL_DTYP"][row] = self.cal_type + self.hdu_cif["CAL_VSD"][row] = self.val_date + self.hdu_cif["CAL_VST"][row] = self.val_time + self.hdu_cif["REF_TIME"][row] = self.ref_time + self.hdu_cif["CAL_QUAL"][row] = self.cal_qual + self.hdu_cif["CAL_DATE"][row] = self.cal_date + + # Add effective area information + row = self.hdu_cif.nrows()-4 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.ea_file + self.hdu_cif["CAL_CNAM"][row] = self.ea_name + self.hdu_cif["CAL_DESC"][row] = self.ea_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.ea_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n] + + # Add point spread function information + row = self.hdu_cif.nrows()-3 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.psf_file + self.hdu_cif["CAL_CNAM"][row] = self.psf_name + self.hdu_cif["CAL_DESC"][row] = self.psf_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.psf_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n] + + # Add energy dispersion information + row = self.hdu_cif.nrows()-2 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.edisp_file + self.hdu_cif["CAL_CNAM"][row] = self.edisp_name + self.hdu_cif["CAL_DESC"][row] = self.edisp_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.edisp_bounds): + 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.bcf_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 root2caldb(self, filename): + """ + Generate files. + + Parameters: + filename - Performance file name + """ + # Open ROOT performance file + rootfile = TFile(filename) + + # Set filenames + ea_filename = self.bcf_path+"/"+self.ea_file + psf_filename = self.bcf_path+"/"+self.psf_file + edisp_filename = self.bcf_path+"/"+self.edisp_file + bgd_filename = self.bcf_path+"/"+self.bgd_file + + # Open ASCII performance tables + if self.ea_file == self.psf_file and \ + self.ea_file == self.edisp_file and \ + self.ea_file == self.bgd_file : + split = False + fprf = open(ea_filename, "w") + files = [fprf] + else: + split = True + fea = open(ea_filename, "w") + fpsf = open(psf_filename, "w") + fedisp = open(edisp_filename, "w") + fbgd = open(bgd_filename, "w") + files = [fea, fpsf, fedisp, fbgd] + + # Write header + for file in files: + file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") + + # Get histograms. Konrad's root files do not have an EffectiveArea80 + # histogram, but his effective area is for 80% containment radius + # anyways, so we simply read this histogram into aeff80. + sens = TH1F() + bgrate = TH1F() + bgratesqdeg = TH1F() + aeff = TH1F() + aeff80 = TH1F() + angres = TH1F() + angres80 = TH1F() + eres = TH1F() + rootfile.GetObject("DiffSens", sens) + rootfile.GetObject("BGRate", bgrate) + rootfile.GetObject("BGRatePerSqDeg", bgratesqdeg) + rootfile.GetObject("EffectiveArea", aeff) + try: + rootfile.GetObject("EffectiveArea80", aeff80) + except: + rootfile.GetObject("EffectiveArea", aeff80) + rootfile.GetObject("AngRes", angres) + rootfile.GetObject("AngRes80", angres80) + rootfile.GetObject("ERes", eres) + + # Setup vectors + energies = bgrate.GetXaxis() + nbins_eng = energies.GetNbins() + for i in range(nbins_eng): + + # Get logE and background rate per squared degree + logE = energies.GetBinCenter(i+1) + bgd = bgrate.GetBinContent(i+1,1) + bgdsq = bgratesqdeg.GetBinContent(i+1,1) + area = aeff.GetBinContent(i+1,1) + area80 = aeff80.GetBinContent(i+1,1) + r68 = angres.GetBinContent(i+1,1) + r80 = angres80.GetBinContent(i+1,1) + dE = eres.GetBinContent(i+1,1) + diffs = sens.GetBinContent(i+1,1) + + # Skip any NaNs + if math.isnan(area80): + continue + + # Compute energy (in MeV) + energy = math.pow(10.0, logE)*1.0e6 + emin = math.pow(10.0, logE-0.1)*1.0e6 + emax = math.pow(10.0, logE+0.1)*1.0e6 + ewidth = emax-emin + + # Convert into background rate per steradian and MeV + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth + + # Compute control background rate (this only works for KB files + # as they are for 80% containment radius). But we don't use this + # anyways, it's just for control on display ... + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) + if omega > 0.0: + bkg_rate2 = bgd / omega / ewidth + else: + bkg_rate2 = 0.0 + + # Compute full effective area by dividing area80/0.80 + area = area80/0.80 + + # Write results in file + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ + (logE, area, r68, r80, dE, bgd, diffs) + for file in files: + file.write(line+"\n") + #fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") + + # Write trailer + for file in files: + file.write("---------------------------------------------\n") + file.write("Notes\n") + file.write(" 1) log(E) = log10(E/TeV) - bin centre\n") + file.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") + file.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") + file.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") + file.write(" 5) Fractional Energy Resolution (rms)\n") + file.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") + file.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") + file.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") + file.write(" systematics and statistics and at least 10 photons.\n") + + # Close files + for file in files: + file.close() + + # Close ROOT performance file + rootfile.Close() + + # Add information to CIF. We do this now as before we did not + # necessarily have all the information at hand (in particular + # about the boundaries) + self.add_cif_info() + + # Return + return + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + """ + """ + # Get ROOT filename from command line argument + #if (len(sys.argv) > 1): + # filename = sys.argv[1] + #else: + db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_20deg_root" + #db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" + + # Get list of all performance files + files = glob.glob(db+"/*.root") + files.sort() + + # Allocate caldb + irf = caldb("kb") + + # Loop over all files + for f in files: + + # Get filename + head, tail = os.path.split(f) + + # Build name + name = tail.strip(".root") + if name.find("kb_") != -1: + name = name.replace("_20deg", "") + elif name.find("IFAE_") != -1: + name = name.replace("Subarray", "") + name = name.replace("IFAE_", "") + name = "IFAE_"+name + index = name.find("hours") + if index != -1: + name = name[0:index+1] + print(tail) + + # Open calibration file + irf.add(name) + + # Translate ROOT to CALDB information + irf.root2caldb(f) diff --git a/test/dev/cta_make_nodes.py b/test/dev/cta_make_nodes.py new file mode 100755 index 0000000..9d8a3dd --- /dev/null +++ b/test/dev/cta_make_nodes.py @@ -0,0 +1,115 @@ +#! /usr/bin/env python +# ========================================================================== +# This Python script creates the node section of the NodeFunction using +# logarithmically spaced energy bins. The intensity scale is set to the +# HESS Crab intensity (assuming a power law). +# +# Copyright (C) 2012-2016 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 +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import sys +import math + + +# ============== # +# Write one node # +# ============== # +def write_node(f, energy, scale): + """ + Writes one node to XML file. + """ + # Convert to strings + str_energy = str(energy) + str_scale = str(scale) + + # Write start tag + f.write(' \n') + + # Write energy + f.write(' \n') + + # Write intensity + f.write(' \n') + + # Write end tag + f.write(' \n') + + # Return + return + + +# ============ # +# Create nodes # +# ============ # +def create_nodes(emin, emax, enumbins): + """ + Create nodes (energies in TeV). + """ + # Open file + f = open("nodes.xml", "w") + + # Set node boundaries + elogmin = math.log10(float(emin)) + elogmax = math.log10(float(emax)) + elogbin = (elogmax - elogmin)/(float(enumbins)-1.0) + + # Fill arrays + for i in range(int(enumbins)): + + # Compute energy + energy = math.pow(10.0, i*elogbin+elogmin) + + # Compute scale (HESS Crab spectrum) + scale = 3.45e-17 * math.pow(energy, -2.63) + + # Write node + write_node(f, energy, scale) + + # Debug + #sys.stdout.write(energy) + #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n") + + # Close file + f.close() + + # Return + return + + +# ======================= # +# Main script entry point # +# ======================= # +if __name__ == '__main__': + + # Check command line + usage = "Usage: cta_make_nodes emin emax enumbins" + if len(sys.argv) < 3: + sys.stdout.write(usage+"\n") + sys.exit() + + # Extract parameters + emin = sys.argv[1] + emax = sys.argv[2] + enumbins = sys.argv[3] + + # Create nodes + create_nodes(emin, emax, enumbins) diff --git a/test/dev/cta_root2irf.py b/test/dev/cta_root2irf.py new file mode 100755 index 0000000..e0f6654 --- /dev/null +++ b/test/dev/cta_root2irf.py @@ -0,0 +1,291 @@ +#! /usr/bin/env python +# ========================================================================== +# This script generates ASCII performance files and file functions for +# background modelling from MC ROOT files. On input it takes the filename +# of the ROOT performance file while on output it produced 2 ASCII files. +# The file "performance.dat" contains the performance data in the usual +# ASCII format, the file "background.txt" contains a file function that can +# be used within GammaLib for spectral background modelling. +# +# Note that the effective area given in the ASCII file is the full effective +# area (i.e. the area before applying any theta cut). This is required by +# GammaLib as the maximum likelihood analysis does not imply any theta cut. +# +# Copyright (C) 2011-2014 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 +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +from ROOT import TFile, TH1F, TH2F +import math +import sys +import glob +import os + + +# ================================== # +# Compute Eq. (17) of Li & Ma (1983) # +# ================================== # +def sigma_lima(Non, Noff, alpha=0.2): + """ + Compute Eq. (17) of Li & Ma (1983). + + Parameters: + Non - Number of on counts + Noff - Number of off counts + Keywords: + alpha - Ratio of on-to-off exposure + """ + # Compute sensitivity + alpha1 = alpha + 1.0 + sum = Non + Noff + arg1 = Non / sum + arg2 = Noff / sum + term1 = Non * math.log((alpha1/alpha)*arg1) + term2 = Noff * math.log(alpha1*arg2) + sigma = math.sqrt(2.0 * (term1 + term2)) + + # Return sensitivity + return sigma + + +# ======================================== # +# Solve Eq. (17) of Li & Ma (1983) for Non # +# ======================================== # +def Non_lima(sigma, Noff, alpha=0.2): + """ + Solve Eq. (17) of Li & Ma (1983) for Non + + Parameters: + sigma - Required significance + Noff - Number of off counts + Keywords: + alpha - Ratio of on-to-off exposure + """ + # Set initial guess for Non + S = 25.0 + Non = S + alpha*Noff + + # Compute initial sigma estimate + s = sigma_lima(Non, Noff, alpha=alpha) + delta = s - sigma + + # Iterate + while abs(delta) > 1.0e-6: + S *= sigma/s + Non = S + alpha*Noff + s = sigma_lima(Non, Noff, alpha=alpha) + delta = s - sigma + + # Return + return Non + + +# =============================================== # +# Compute effective area assuming a Crab spectrum # +# =============================================== # +def get_aeff(bgrate, emin, emax, \ + duration=180000.0, alpha=0.2, sigma=5.0, norm=2.79e-7, index=-2.57): + """ + Computes effective area assuming a Crab spectrum. This does not + include any hadron cut efficiencies, and thus is not very useful in a first + place ... + + Parameters: + bgrate - Background rate (TeV) + emin - Minimum energy (TeV) + emax - Maximum energy (TeV) + Keywords: + duration - Observing time in seconds (default: 50h) + alpha - Non/Noff ratio (default: 0.2) + sigma - Significance (default: 5.0) + norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7) + index - Crab spectral index (default: -2.57) + """ + # Compute number of source events + B = bgrate*duration + Noff = B/alpha + Non = Non_lima(sigma, Noff) + S = (Non - B) + + # Compute Crab flux in specified energy band + inx = index+1.0 + cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx)) + + # Compute effective area + area = S / duration / cu + + # Return effective area + return area + + +# ============== # +# Generate files # +# ============== # +def make_files(filename, prfname, bkgname): + """ + Generate files. + + Parameters: + filename - Performance file name + prfname - Performance file name + bkgname - File function filename + """ + # Open performance file and file function + file = TFile(filename) + fprf = open(prfname, "w") + fbkg = open(bkgname, "w") + + # Write header + fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") + #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") + + # Get histograms. Konrad's root files do not have an EffectiveArea80 + # histogram, but his effective area is for 80% containment radius + # anyways, so we simply read this histogram into aeff80. + sens = TH1F() + bgrate = TH1F() + bgratesqdeg = TH1F() + aeff = TH1F() + aeff80 = TH1F() + angres = TH1F() + angres80 = TH1F() + eres = TH1F() + file.GetObject("DiffSens", sens) + file.GetObject("BGRate", bgrate) + file.GetObject("BGRatePerSqDeg", bgratesqdeg) + file.GetObject("EffectiveArea", aeff) + try: + file.GetObject("EffectiveArea80", aeff80) + except: + file.GetObject("EffectiveArea", aeff80) + file.GetObject("AngRes", angres) + file.GetObject("AngRes80", angres80) + file.GetObject("ERes", eres) + + # Setup vectors + energies = bgrate.GetXaxis() + nbins_eng = energies.GetNbins() + for i in range(nbins_eng): + + # Get logE and background rate per squared degree + logE = energies.GetBinCenter(i+1) + bgd = bgrate.GetBinContent(i+1,1) + bgdsq = bgratesqdeg.GetBinContent(i+1,1) + area = aeff.GetBinContent(i+1,1) + area80 = aeff80.GetBinContent(i+1,1) + r68 = angres.GetBinContent(i+1,1) + r80 = angres80.GetBinContent(i+1,1) + dE = eres.GetBinContent(i+1,1) + diffs = sens.GetBinContent(i+1,1) + + # Skip any NaNs + if math.isnan(area80): + continue + + # Compute energy (in MeV) + energy = math.pow(10.0, logE)*1.0e6 + emin = math.pow(10.0, logE-0.1)*1.0e6 + emax = math.pow(10.0, logE+0.1)*1.0e6 + ewidth = emax-emin + + # Convert into background rate per steradian and MeV + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth + + # Compute control background rate (this only works for KB files + # as they are for 80% containment radius). But we don't use this + # anyways, it's just for control on display ... + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) + if omega > 0.0: + bkg_rate2 = bgd / omega / ewidth + else: + bkg_rate2 = 0.0 + + # Compute full effective area by dividing area80/0.80 + area = area80/0.80 + + # Write results in file + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ + (logE, area, r68, r80, dE, bgd, diffs) + fprf.write(line+"\n") + fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") + + # Show performance file + #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \ + # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2) + + # Write trailer + fprf.write("---------------------------------------------\n") + fprf.write("Notes\n") + fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n") + fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") + fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") + fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") + fprf.write(" 5) Fractional Energy Resolution (rms)\n") + fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") + fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") + fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") + fprf.write(" systematics and statistics and at least 10 photons.\n") + + # Close files + fprf.close() + fbkg.close() + + # Return + return + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + """ + """ + # Get ROOT filename from command line argument + if (len(sys.argv) > 1): + filename = sys.argv[1] + else: + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010" + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root" + caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" + + # Get list of all performance files + files = glob.glob(caldb+"/*.root") + files.sort() + + # Loop over all files + for file in files: + + # Get filename + head, tail = os.path.split(file) + + # Build name + name = tail.strip(".root") + if name.find("kb_") != -1: + name = name.replace("_20deg", "") + elif name.find("IFAE_") != -1: + name = name.replace("Subarray", "") + name = name.replace("IFAE_", "") + name = "IFAE_"+name + i = name.find("hours") + if i != -1: + name = name[0:i+1] + + # Build output filenames + prfname = name + ".dat" + bkgname = "bkg_"+name+".txt" + #print prfname, bkgname + + # Create files + make_files(file, prfname, bkgname) -- libgit2 0.21.4