Alexander Ziegler / ctools

Browse Code »

Commit 90e04b2dbfd13b664b6c09b60a5c5327ce115ca7

Authored by Jürgen Knödlseder 2014-02-20 11:52:00 +0100

1 parent 8d74cdab

Add script that converts Konrad's root performance files into a caldb format.

Inline Side-by-side

Showing 1 changed file with 668 additions and 0 deletions

scripts/cta_kb_root2caldb.py 0 → 100755

View file @90e04b2

	1	+#! /usr/bin/env python
	2	+# ==========================================================================
	3	+# This script generates ASCII performance files and file functions for
	4	+# background modelling from MC ROOT files. On input it takes the filename
	5	+# of the ROOT performance file while on output it produced 2 ASCII files.
	6	+# The file "performance.dat" contains the performance data in the usual
	7	+# ASCII format, the file "background.txt" contains a file function that can
	8	+# be used within GammaLib for spectral background modelling.
	9	+#
	10	+# Note that the effective area given in the ASCII file is the full effective
	11	+# area (i.e. the area before applying any theta cut). This is required by
	12	+# GammaLib as the maximum likelihood analysis does not imply any theta cut.
	13	+#
	14	+# Copyright (C) 2011-2014 Juergen Knoedlseder
	15	+#
	16	+# This program is free software: you can redistribute it and/or modify
	17	+# it under the terms of the GNU General Public License as published by
	18	+# the Free Software Foundation, either version 3 of the License, or
	19	+# (at your option) any later version.
	20	+#
	21	+# This program is distributed in the hope that it will be useful,
	22	+# but WITHOUT ANY WARRANTY; without even the implied warranty of
	23	+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	24	+# GNU General Public License for more details.
	25	+#
	26	+# You should have received a copy of the GNU General Public License
	27	+# along with this program. If not, see <http://www.gnu.org/licenses/>.
	28	+#
	29	+# ==========================================================================
	30	+from ROOT import TFile, TH1F, TH2F
	31	+from datetime import datetime
	32	+import gammalib
	33	+import math
	34	+import sys
	35	+import glob
	36	+import os
	37	+
	38	+
	39	+# =========== #
	40	+# CalDB class #
	41	+# =========== #
	42	+class caldb():
	43	+ """
	44	+ """
	45	+ def __init__(self, inst, path=""):
	46	+ """
	47	+ CALDB constructor. The constructor will create one CALDB entry for
	48	+ an observation.
	49	+
	50	+ Parameters:
	51	+ inst - Instrument string.
	52	+ Keywords:
	53	+ path - CALDB Root path.
	54	+ """
	55	+ # Store root path
	56	+ self.path = path
	57	+
	58	+ # Initialise some members
	59	+ self.cif = None
	60	+ self.irf = None
	61	+ self.ea = None
	62	+ self.psf = None
	63	+ self.edisp = None
	64	+ self.bgd = None
	65	+ self.hdu_cif = None
	66	+
	67	+ # Store UTC date of job execution
	68	+ self.utc = datetime.utcnow().isoformat()[:19]
	69	+
	70	+ # Set CALDB information
	71	+ self.cal_tel = "CTA"
	72	+ self.cal_inst = inst.upper()
	73	+ self.cal_det = "NONE"
	74	+ self.cal_flt = "NONE"
	75	+ self.cal_class = "BCF"
	76	+ self.cal_type = "DATA"
	77	+ self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ...
	78	+ self.cal_date = "14/01/01"
	79	+ self.val_date = "2014-01-01"
	80	+ self.val_time = "00:00:00"
	81	+ self.ref_time = 51544.0
	82	+ self.cal_name = "NAME()"
	83	+ self.cal_version = "VERSION(Prod1)"
	84	+ self.cal_cut = "CLASS(BEST)"
	85	+ self.cal_analysis = "ANALYSIS(MPIK)"
	86	+ #
	87	+ self.ea_name = "EFF_AREA"
	88	+ self.ea_doc = "???"
	89	+ self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	90	+ self.ea_desc = "CTA effective area"
	91	+ #
	92	+ self.psf_name = "RPSF"
	93	+ self.psf_doc = "???"
	94	+ self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	95	+ self.psf_desc = "CTA point spread function"
	96	+ #
	97	+ self.edisp_name = "EDISP"
	98	+ self.edisp_doc = "???"
	99	+ self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	100	+ self.edisp_desc = "CTA energy dispersion"
	101	+ #
	102	+ self.bgd_name = "BGD"
	103	+ self.bgd_doc = "???"
	104	+ self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	105	+ self.bgd_desc = "CTA background"
	106	+
	107	+ # Create directory structure
	108	+ self.make_dirs()
	109	+
	110	+ # Return
	111	+ return
	112	+
	113	+ def __del__(self):
	114	+ """
	115	+ Destructor. Close all FITS files.
	116	+ """
	117	+ # Close calibration
	118	+ self.close()
	119	+
	120	+ # Return
	121	+ return
	122	+
	123	+ def make_dirs(self):
	124	+ """
	125	+ Generate CALDB directory structure for one observation identifier.
	126	+ The structure is given by
	127	+
	128	+ data/<tel>/<inst>/bcf/
	129	+
	130	+ where <tel> is "cta" and <inst> is the instrument specified in the
	131	+ CALDB constructor (the instrument may be used for different array
	132	+ configurations).
	133	+
	134	+ Parameters:
	135	+ None
	136	+ Keywords:
	137	+ None
	138	+ """
	139	+ # Set base path
	140	+ self.base_dir = "data"
	141	+ self.base_dir += "/"+self.cal_tel.lower()
	142	+ self.base_dir += "/"+self.cal_inst.lower()
	143	+
	144	+ # Set directory names
	145	+ self.bcf_dir = self.base_dir+"/bcf"
	146	+
	147	+ # Optionally prefix path
	148	+ if len(self.path) > 0:
	149	+ self.base_path = self.path+"/"+self.base_dir
	150	+ self.bcf_path = self.path+"/"+self.bcf_dir
	151	+ else:
	152	+ self.base_path = self.base_dir
	153	+ self.bcf_path = self.bcf_dir
	154	+
	155	+ # Create directory structure
	156	+ if not os.path.isdir(self.bcf_path):
	157	+ os.makedirs(self.bcf_path)
	158	+
	159	+ # Return
	160	+ return
	161	+
	162	+ def add(self, rspname, split=False, clobber=True):
	163	+ """
	164	+ Add new calibration. The actual version will put
	165	+ all calibrations in the same file, although each part of the response
	166	+ function will have its own logical name. We can thus easily modify
	167	+ the script to put each calibration information in a separate file.
	168	+
	169	+ Parameters:
	170	+ rspname - Response name
	171	+ Keywords:
	172	+ split - Split IRF over several files?
	173	+ clobber - Overwrite existing files?
	174	+ """
	175	+ # Set calibrate file names
	176	+ if split:
	177	+ self.ea_file = "ea_"+rspname+".dat"
	178	+ self.psf_file = "psf_"+rspname+".dat"
	179	+ self.edisp_file = "edisp_"+rspname+".dat"
	180	+ self.bgd_file = "bgd_"+rspname+".dat"
	181	+ else:
	182	+ self.ea_file = "irf_"+rspname+".dat"
	183	+ self.psf_file = "irf_"+rspname+".dat"
	184	+ self.edisp_file = "irf_"+rspname+".dat"
	185	+ self.bgd_file = "irf_"+rspname+".dat"
	186	+
	187	+ # Open calibration database index
	188	+ if self.cif == None:
	189	+ self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True)
	190	+
	191	+ # If file has no CIF extension than create it now
	192	+ try:
	193	+ self.hdu_cif = self.cif.table("CIF")
	194	+ except:
	195	+ self.create_cif()
	196	+ self.hdu_cif = self.cif.table("CIF")
	197	+
	198	+ # Set response name
	199	+ self.cal_name = "NAME("+rspname+")"
	200	+
	201	+ # Return
	202	+ return
	203	+
	204	+ def close(self):
	205	+ """
	206	+ Close calibration FITS files.
	207	+
	208	+ Parameters:
	209	+ None
	210	+ Keywords:
	211	+ None
	212	+ """
	213	+ # Add information to CIF. We do this now as before we did not
	214	+ # necessarily have all the information at hand (in particular
	215	+ # about the boundaries)
	216	+ #self.add_cif_info()
	217	+
	218	+ # Close CIF
	219	+ if self.cif != None:
	220	+ self.cif.save(True)
	221	+ self.cif.close()
	222	+ self.cif = None
	223	+
	224	+ # Free
	225	+ self.cif = None
	226	+
	227	+ # Return
	228	+ return
	229	+
	230	+ def create_cif(self):
	231	+ """
	232	+ Create Calibration Index File (CIF) extension in CIF FITS file.
	233	+
	234	+ Parameters:
	235	+ None
	236	+ Keywords:
	237	+ None
	238	+ """
	239	+ # Create binary table
	240	+ table = gammalib.GFitsBinTable()
	241	+
	242	+ # Set boundary
	243	+
	244	+ # Attach columns. Reference: CAL/GEN/92-008
	245	+ table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
	246	+ table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
	247	+ table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20))
	248	+ table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10))
	249	+ table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20))
	250	+ table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70))
	251	+ table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40))
	252	+ table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3))
	253	+ table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4))
	254	+ table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20))
	255	+ table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9))
	256	+ table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0))
	257	+ table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10))
	258	+ table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8))
	259	+ table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0))
	260	+ table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
	261	+ table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
	262	+ table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
	263	+
	264	+ # Set keywords. Reference: CAL/GEN/92-008
	265	+ table.extname("CIF")
	266	+ table.card("CIFVERSN", "1992a", "Version of CIF format")
	267	+
	268	+ # Attach table to FITS file. Note that at this moment the
	269	+ # FITS table gets connected to the FITS file. Yet since
	270	+ # nothing was yet written to the FITS file, we cannot read
	271	+ # anything from it.
	272	+ self.cif.append(table)
	273	+
	274	+ # Return
	275	+ return
	276	+
	277	+ def add_cif_info(self):
	278	+ """
	279	+ Add information to CIF extension.
	280	+
	281	+ Parameters:
	282	+ None
	283	+ Keywords:
	284	+ None
	285	+ """
	286	+ # Set boundaries
	287	+ self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	288	+ self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	289	+ self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	290	+ self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
	291	+
	292	+ # Append 4 rows to CIF extension
	293	+ self.hdu_cif.append_rows(4)
	294	+
	295	+ # Add generic information for these 4 rows
	296	+ for i in range(4):
	297	+
	298	+ # Set row number
	299	+ row = i+self.hdu_cif.nrows()-4
	300	+
	301	+ # Set element
	302	+ self.hdu_cif["TELESCOP"][row] = self.cal_tel
	303	+ self.hdu_cif["INSTRUME"][row] = self.cal_inst
	304	+ self.hdu_cif["DETNAM"][row] = self.cal_det
	305	+ self.hdu_cif["FILTER"][row] = self.cal_flt
	306	+ self.hdu_cif["CAL_DEV"][row] = "ONLINE"
	307	+ self.hdu_cif["CAL_CLAS"][row] = self.cal_class
	308	+ self.hdu_cif["CAL_DTYP"][row] = self.cal_type
	309	+ self.hdu_cif["CAL_VSD"][row] = self.val_date
	310	+ self.hdu_cif["CAL_VST"][row] = self.val_time
	311	+ self.hdu_cif["REF_TIME"][row] = self.ref_time
	312	+ self.hdu_cif["CAL_QUAL"][row] = self.cal_qual
	313	+ self.hdu_cif["CAL_DATE"][row] = self.cal_date
	314	+
	315	+ # Add effective area information
	316	+ row = self.hdu_cif.nrows()-4
	317	+ self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
	318	+ self.hdu_cif["CAL_FILE"][row] = self.ea_file
	319	+ self.hdu_cif["CAL_CNAM"][row] = self.ea_name
	320	+ self.hdu_cif["CAL_DESC"][row] = self.ea_desc
	321	+ self.hdu_cif["CAL_XNO"][row] = 1
	322	+ for n in range(9):
	323	+ if n >= len(self.ea_bounds):
	324	+ self.hdu_cif["CAL_CBD"][row, n] = "NONE"
	325	+ else:
	326	+ self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n]
	327	+
	328	+ # Add point spread function information
	329	+ row = self.hdu_cif.nrows()-3
	330	+ self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
	331	+ self.hdu_cif["CAL_FILE"][row] = self.psf_file
	332	+ self.hdu_cif["CAL_CNAM"][row] = self.psf_name
	333	+ self.hdu_cif["CAL_DESC"][row] = self.psf_desc
	334	+ self.hdu_cif["CAL_XNO"][row] = 1
	335	+ for n in range(9):
	336	+ if n >= len(self.psf_bounds):
	337	+ self.hdu_cif["CAL_CBD"][row, n] = "NONE"
	338	+ else:
	339	+ self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n]
	340	+
	341	+ # Add energy dispersion information
	342	+ row = self.hdu_cif.nrows()-2
	343	+ self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
	344	+ self.hdu_cif["CAL_FILE"][row] = self.edisp_file
	345	+ self.hdu_cif["CAL_CNAM"][row] = self.edisp_name
	346	+ self.hdu_cif["CAL_DESC"][row] = self.edisp_desc
	347	+ self.hdu_cif["CAL_XNO"][row] = 1
	348	+ for n in range(9):
	349	+ if n >= len(self.edisp_bounds):
	350	+ self.hdu_cif["CAL_CBD"][row, n] = "NONE"
	351	+ else:
	352	+ self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n]
	353	+
	354	+ # Add background information
	355	+ row = self.hdu_cif.nrows()-1
	356	+ self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
	357	+ self.hdu_cif["CAL_FILE"][row] = self.bgd_file
	358	+ self.hdu_cif["CAL_CNAM"][row] = self.bgd_name
	359	+ self.hdu_cif["CAL_DESC"][row] = self.bgd_desc
	360	+ self.hdu_cif["CAL_XNO"][row] = 1
	361	+ for n in range(9):
	362	+ if n >= len(self.bgd_bounds):
	363	+ self.hdu_cif["CAL_CBD"][row, n] = "NONE"
	364	+ else:
	365	+ self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
	366	+
	367	+ # Return
	368	+ return
	369	+
	370	+ def root2caldb(self, filename):
	371	+ """
	372	+ Generate files.
	373	+
	374	+ Parameters:
	375	+ filename - Performance file name
	376	+ """
	377	+ # Open ROOT performance file
	378	+ rootfile = TFile(filename)
	379	+
	380	+ # Set filenames
	381	+ ea_filename = self.bcf_path+"/"+self.ea_file
	382	+ psf_filename = self.bcf_path+"/"+self.psf_file
	383	+ edisp_filename = self.bcf_path+"/"+self.edisp_file
	384	+ bgd_filename = self.bcf_path+"/"+self.bgd_file
	385	+
	386	+ # Open ASCII performance tables
	387	+ if self.ea_file == self.psf_file and \
	388	+ self.ea_file == self.edisp_file and \
	389	+ self.ea_file == self.bgd_file :
	390	+ split = False
	391	+ fprf = open(ea_filename, "w")
	392	+ files = [fprf]
	393	+ else:
	394	+ split = True
	395	+ fea = open(ea_filename, "w")
	396	+ fpsf = open(psf_filename, "w")
	397	+ fedisp = open(edisp_filename, "w")
	398	+ fbgd = open(bgd_filename, "w")
	399	+ files = [fea, fpsf, fedisp, fbgd]
	400	+
	401	+ # Write header
	402	+ for file in files:
	403	+ file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
	404	+ file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
	405	+ file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
	406	+ file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
	407	+
	408	+ # Get histograms. Konrad's root files do not have an EffectiveArea80
	409	+ # histogram, but his effective area is for 80% containment radius
	410	+ # anyways, so we simply read this histogram into aeff80.
	411	+ sens = TH1F()
	412	+ bgrate = TH1F()
	413	+ bgratesqdeg = TH1F()
	414	+ aeff = TH1F()
	415	+ aeff80 = TH1F()
	416	+ angres = TH1F()
	417	+ angres80 = TH1F()
	418	+ eres = TH1F()
	419	+ rootfile.GetObject("DiffSens", sens)
	420	+ rootfile.GetObject("BGRate", bgrate)
	421	+ rootfile.GetObject("BGRatePerSqDeg", bgratesqdeg)
	422	+ rootfile.GetObject("EffectiveArea", aeff)
	423	+ try:
	424	+ rootfile.GetObject("EffectiveArea80", aeff80)
	425	+ except:
	426	+ rootfile.GetObject("EffectiveArea", aeff80)
	427	+ rootfile.GetObject("AngRes", angres)
	428	+ rootfile.GetObject("AngRes80", angres80)
	429	+ rootfile.GetObject("ERes", eres)
	430	+
	431	+ # Setup vectors
	432	+ energies = bgrate.GetXaxis()
	433	+ nbins_eng = energies.GetNbins()
	434	+ for i in range(nbins_eng):
	435	+
	436	+ # Get logE and background rate per squared degree
	437	+ logE = energies.GetBinCenter(i+1)
	438	+ bgd = bgrate.GetBinContent(i+1,1)
	439	+ bgdsq = bgratesqdeg.GetBinContent(i+1,1)
	440	+ area = aeff.GetBinContent(i+1,1)
	441	+ area80 = aeff80.GetBinContent(i+1,1)
	442	+ r68 = angres.GetBinContent(i+1,1)
	443	+ r80 = angres80.GetBinContent(i+1,1)
	444	+ dE = eres.GetBinContent(i+1,1)
	445	+ diffs = sens.GetBinContent(i+1,1)
	446	+
	447	+ # Skip any NaNs
	448	+ if math.isnan(area80):
	449	+ continue
	450	+
	451	+ # Compute energy (in MeV)
	452	+ energy = math.pow(10.0, logE)*1.0e6
	453	+ emin = math.pow(10.0, logE-0.1)*1.0e6
	454	+ emax = math.pow(10.0, logE+0.1)*1.0e6
	455	+ ewidth = emax-emin
	456	+
	457	+ # Convert into background rate per steradian and MeV
	458	+ bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
	459	+
	460	+ # Compute control background rate (this only works for KB files
	461	+ # as they are for 80% containment radius). But we don't use this
	462	+ # anyways, it's just for control on display ...
	463	+ omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
	464	+ if omega > 0.0:
	465	+ bkg_rate2 = bgd / omega / ewidth
	466	+ else:
	467	+ bkg_rate2 = 0.0
	468	+
	469	+ # Compute full effective area by dividing area80/0.80
	470	+ area = area80/0.80
	471	+
	472	+ # Write results in file
	473	+ line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \
	474	+ (logE, area, r68, r80, dE, bgd, diffs)
	475	+ for file in files:
	476	+ file.write(line+"\n")
	477	+ #fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
	478	+
	479	+ # Write trailer
	480	+ for file in files:
	481	+ file.write("---------------------------------------------\n")
	482	+ file.write("Notes\n")
	483	+ file.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
	484	+ file.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
	485	+ file.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
	486	+ file.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
	487	+ file.write(" 5) Fractional Energy Resolution (rms)\n")
	488	+ file.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
	489	+ file.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
	490	+ file.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
	491	+ file.write(" systematics and statistics and at least 10 photons.\n")
	492	+
	493	+ # Close files
	494	+ for file in files:
	495	+ file.close()
	496	+
	497	+ # Close ROOT performance file
	498	+ rootfile.Close()
	499	+
	500	+ # Add information to CIF. We do this now as before we did not
	501	+ # necessarily have all the information at hand (in particular
	502	+ # about the boundaries)
	503	+ self.add_cif_info()
	504	+
	505	+ # Return
	506	+ return
	507	+
	508	+
	509	+# ============== #
	510	+# Generate files #
	511	+# ============== #
	512	+def make_files(filename, prfname, bkgname):
	513	+ """
	514	+ Generate files.
	515	+
	516	+ Parameters:
	517	+ filename - Performance file name
	518	+ prfname - Performance file name
	519	+ bkgname - File function filename
	520	+ """
	521	+ # Open performance file and file function
	522	+ file = TFile(filename)
	523	+ fprf = open(prfname, "w")
	524	+ fbkg = open(bkgname, "w")
	525	+
	526	+ # Write header
	527	+ fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
	528	+ #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n")
	529	+
	530	+ # Get histograms. Konrad's root files do not have an EffectiveArea80
	531	+ # histogram, but his effective area is for 80% containment radius
	532	+ # anyways, so we simply read this histogram into aeff80.
	533	+ sens = TH1F()
	534	+ bgrate = TH1F()
	535	+ bgratesqdeg = TH1F()
	536	+ aeff = TH1F()
	537	+ aeff80 = TH1F()
	538	+ angres = TH1F()
	539	+ angres80 = TH1F()
	540	+ eres = TH1F()
	541	+ file.GetObject("DiffSens", sens)
	542	+ file.GetObject("BGRate", bgrate)
	543	+ file.GetObject("BGRatePerSqDeg", bgratesqdeg)
	544	+ file.GetObject("EffectiveArea", aeff)
	545	+ try:
	546	+ file.GetObject("EffectiveArea80", aeff80)
	547	+ except:
	548	+ file.GetObject("EffectiveArea", aeff80)
	549	+ file.GetObject("AngRes", angres)
	550	+ file.GetObject("AngRes80", angres80)
	551	+ file.GetObject("ERes", eres)
	552	+
	553	+ # Setup vectors
	554	+ energies = bgrate.GetXaxis()
	555	+ nbins_eng = energies.GetNbins()
	556	+ for i in range(nbins_eng):
	557	+
	558	+ # Get logE and background rate per squared degree
	559	+ logE = energies.GetBinCenter(i+1)
	560	+ bgd = bgrate.GetBinContent(i+1,1)
	561	+ bgdsq = bgratesqdeg.GetBinContent(i+1,1)
	562	+ area = aeff.GetBinContent(i+1,1)
	563	+ area80 = aeff80.GetBinContent(i+1,1)
	564	+ r68 = angres.GetBinContent(i+1,1)
	565	+ r80 = angres80.GetBinContent(i+1,1)
	566	+ dE = eres.GetBinContent(i+1,1)
	567	+ diffs = sens.GetBinContent(i+1,1)
	568	+
	569	+ # Skip any NaNs
	570	+ if math.isnan(area80):
	571	+ continue
	572	+
	573	+ # Compute energy (in MeV)
	574	+ energy = math.pow(10.0, logE)*1.0e6
	575	+ emin = math.pow(10.0, logE-0.1)*1.0e6
	576	+ emax = math.pow(10.0, logE+0.1)*1.0e6
	577	+ ewidth = emax-emin
	578	+
	579	+ # Convert into background rate per steradian and MeV
	580	+ bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
	581	+
	582	+ # Compute control background rate (this only works for KB files
	583	+ # as they are for 80% containment radius). But we don't use this
	584	+ # anyways, it's just for control on display ...
	585	+ omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
	586	+ if omega > 0.0:
	587	+ bkg_rate2 = bgd / omega / ewidth
	588	+ else:
	589	+ bkg_rate2 = 0.0
	590	+
	591	+ # Compute full effective area by dividing area80/0.80
	592	+ area = area80/0.80
	593	+
	594	+ # Write results in file
	595	+ line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \
	596	+ (logE, area, r68, r80, dE, bgd, diffs)
	597	+ fprf.write(line+"\n")
	598	+ fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
	599	+
	600	+ # Show performance file
	601	+ #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \
	602	+ # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2)
	603	+
	604	+ # Write trailer
	605	+ fprf.write("---------------------------------------------\n")
	606	+ fprf.write("Notes\n")
	607	+ fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
	608	+ fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
	609	+ fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
	610	+ fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
	611	+ fprf.write(" 5) Fractional Energy Resolution (rms)\n")
	612	+ fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
	613	+ fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
	614	+ fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
	615	+ fprf.write(" systematics and statistics and at least 10 photons.\n")
	616	+
	617	+ # Close files
	618	+ fprf.close()
	619	+ fbkg.close()
	620	+
	621	+ # Return
	622	+ return
	623	+
	624	+
	625	+#==========================#
	626	+# Main routine entry point #
	627	+#==========================#
	628	+if __name__ == '__main__':
	629	+ """
	630	+ """
	631	+ # Get ROOT filename from command line argument
	632	+ if (len(sys.argv) > 1):
	633	+ filename = sys.argv[1]
	634	+ else:
	635	+ db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_20deg_root"
	636	+ #db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root"
	637	+
	638	+ # Get list of all performance files
	639	+ files = glob.glob(db+"/*.root")
	640	+ files.sort()
	641	+
	642	+ # Allocate caldb
	643	+ irf = caldb("kb")
	644	+
	645	+ # Loop over all files
	646	+ for file in files:
	647	+
	648	+ # Get filename
	649	+ head, tail = os.path.split(file)
	650	+
	651	+ # Build name
	652	+ name = tail.strip(".root")
	653	+ if name.find("kb_") != -1:
	654	+ name = name.replace("_20deg", "")
	655	+ elif name.find("IFAE_") != -1:
	656	+ name = name.replace("Subarray", "")
	657	+ name = name.replace("IFAE_", "")
	658	+ name = "IFAE_"+name
	659	+ i = name.find("hours")
	660	+ if i != -1:
	661	+ name = name[0:i+1]
	662	+ print(tail)
	663	+
	664	+ # Open calibration file
	665	+ irf.add(name)
	666	+
	667	+ # Translate ROOT to CALDB information
	668	+ irf.root2caldb(file)
...	...