Andreas Specovius / ctools

Commit d5049a1a6fae633f25f4cfd19fdf6903e6831c0b

Authored by Jürgen Knödlseder
1 parent 5ce18429

Update to conform with agreed response format.

Inline Side-by-side
Showing 1 changed file with 456 additions and 89 deletions
scripts/cta_root2caldb.py
  View file @d5049a1
@@ -2,11 +2,9 @@ @@ -2,11 +2,9 @@
2 # ========================================================================== 2 # ==========================================================================
3 # This script generates IRFs in the CALDB format of HEASARC using ROOT 2D 3 # This script generates IRFs in the CALDB format of HEASARC using ROOT 2D
4 # performance files. 4 # performance files.
5 -#  
6 -# Requirements:  
7 -# - pyroot  
8 -#  
9 -# Copyright (C) 2011 Jurgen Knodlseder 5 +# -------------------------------------------------------------------------
  6 +# Copyright (C) 2011-2012 Juergen Knoedlseder
  7 +# -------------------------------------------------------------------------
10 # 8 #
11 # This program is free software: you can redistribute it and/or modify 9 # This program is free software: you can redistribute it and/or modify
12 # it under the terms of the GNU General Public License as published by 10 # it under the terms of the GNU General Public License as published by
@@ -21,6 +19,10 @@ @@ -21,6 +19,10 @@
21 # You should have received a copy of the GNU General Public License 19 # You should have received a copy of the GNU General Public License
22 # along with this program. If not, see <http://www.gnu.org/licenses/>. 20 # along with this program. If not, see <http://www.gnu.org/licenses/>.
23 # 21 #
  22 +# -------------------------------------------------------------------------
  23 +# Requirements:
  24 +# - gammalib
  25 +# - pyroot
24 # ========================================================================== 26 # ==========================================================================
25 from ROOT import TFile, TH1F, TH2F 27 from ROOT import TFile, TH1F, TH2F
26 from gammalib import * 28 from gammalib import *
@@ -37,18 +39,27 @@ import os @@ -37,18 +39,27 @@ import os
37 class caldb(): 39 class caldb():
38 """ 40 """
39 """ 41 """
40 - def __init__(self, conf, path=""): 42 + def __init__(self, inst, obsid, path=""):
41 """ 43 """
42 - Constructor. 44 + CALDB constructor. The constructor will create one CALDB entry for
  45 + an observation.
  46 +
  47 + Parameters:
  48 + inst - Instrument string.
  49 + obsid - Observation ID string.
  50 + Keywords:
  51 + path - CALDB Root path.
43 """ 52 """
44 # Store root path 53 # Store root path
45 self.path = path 54 self.path = path
46 55
47 # Initialise some members 56 # Initialise some members
48 self.cif = None 57 self.cif = None
  58 + self.irf = None
49 self.ea = None 59 self.ea = None
50 self.psf = None 60 self.psf = None
51 self.edisp = None 61 self.edisp = None
  62 + self.bgd = None
52 self.hdu_cif = None 63 self.hdu_cif = None
53 64
54 # Store UTC date of job execution 65 # Store UTC date of job execution
@@ -56,7 +67,8 @@ class caldb(): @@ -56,7 +67,8 @@ class caldb():
56 67
57 # Set CALDB information 68 # Set CALDB information
58 self.cal_tel = "CTA" 69 self.cal_tel = "CTA"
59 - self.cal_inst = conf.upper() 70 + self.cal_inst = inst.upper()
  71 + self.cal_obsid = obsid
60 self.cal_det = "NONE" 72 self.cal_det = "NONE"
61 self.cal_flt = "NONE" 73 self.cal_flt = "NONE"
62 self.cal_class = "BCF" 74 self.cal_class = "BCF"
@@ -84,6 +96,11 @@ class caldb(): @@ -84,6 +96,11 @@ class caldb():
84 self.edisp_doc = "???" 96 self.edisp_doc = "???"
85 self.edisp_bounds = [self.cal_version, self.cal_cut, self.cal_analysis] 97 self.edisp_bounds = [self.cal_version, self.cal_cut, self.cal_analysis]
86 self.edisp_desc = "CTA energy dispersion" 98 self.edisp_desc = "CTA energy dispersion"
  99 + #
  100 + self.bgd_name = "BGD"
  101 + self.bgd_doc = "???"
  102 + self.bgd_bounds = [self.cal_version, self.cal_cut, self.cal_analysis]
  103 + self.bgd_desc = "CTA background"
87 104
88 # Create directory structure 105 # Create directory structure
89 self.make_dirs() 106 self.make_dirs()
@@ -93,6 +110,7 @@ class caldb(): @@ -93,6 +110,7 @@ class caldb():
93 110
94 def __del__(self): 111 def __del__(self):
95 """ 112 """
  113 + Destructor. Close all FITS files.
96 """ 114 """
97 # Close calibration 115 # Close calibration
98 self.close() 116 self.close()
@@ -102,32 +120,50 @@ class caldb(): @@ -102,32 +120,50 @@ class caldb():
102 120
103 def make_dirs(self): 121 def make_dirs(self):
104 """ 122 """
  123 + Generate CALDB directory structure for one observation identifier.
  124 + The structure is given by
  125 +
  126 + data/<tel>/<inst>/bcf/<obsid>
  127 +
  128 + where <tel> is "cta" and <inst> is the instrument specified in the
  129 + CALDB constructor (the instrument may be used for different array
  130 + configurations).
  131 +
  132 + Parameters:
  133 + None
  134 + Keywords:
  135 + None
105 """ 136 """
106 # Set base path 137 # Set base path
107 self.base_dir = "data" 138 self.base_dir = "data"
108 self.base_dir += "/"+self.cal_tel.lower() 139 self.base_dir += "/"+self.cal_tel.lower()
109 self.base_dir += "/"+self.cal_inst.lower() 140 self.base_dir += "/"+self.cal_inst.lower()
110 141
111 -  
112 # Set directory names 142 # Set directory names
113 self.bcf_dir = self.base_dir+"/bcf" 143 self.bcf_dir = self.base_dir+"/bcf"
114 - self.ea_dir = self.bcf_dir+"/ea"  
115 - self.psf_dir = self.bcf_dir+"/psf"  
116 - self.edisp_dir = self.bcf_dir+"/edisp" 144 + self.obs_dir = self.bcf_dir+"/"+self.cal_obsid
  145 + self.ea_dir = self.obs_dir
  146 + self.psf_dir = self.obs_dir
  147 + self.edisp_dir = self.obs_dir
  148 + self.bgd_dir = self.obs_dir
117 149
118 - # Add path (if required) 150 + # Optionally prefix path
119 if len(self.path) > 0: 151 if len(self.path) > 0:
120 self.base_path = self.path+"/"+self.base_dir 152 self.base_path = self.path+"/"+self.base_dir
121 self.bcf_path = self.path+"/"+self.bcf_dir 153 self.bcf_path = self.path+"/"+self.bcf_dir
  154 + self.obs_path = self.path+"/"+self.obs_dir
122 self.ea_path = self.path+"/"+self.ea_dir 155 self.ea_path = self.path+"/"+self.ea_dir
123 self.psf_path = self.path+"/"+self.psf_dir 156 self.psf_path = self.path+"/"+self.psf_dir
124 self.edisp_path = self.path+"/"+self.edisp_dir 157 self.edisp_path = self.path+"/"+self.edisp_dir
  158 + self.bgd_path = self.path+"/"+self.bgd_dir
125 else: 159 else:
126 self.base_path = self.base_dir 160 self.base_path = self.base_dir
127 self.bcf_path = self.bcf_dir 161 self.bcf_path = self.bcf_dir
  162 + self.obs_path = self.obs_dir
128 self.ea_path = self.ea_dir 163 self.ea_path = self.ea_dir
129 self.psf_path = self.psf_dir 164 self.psf_path = self.psf_dir
130 self.edisp_path = self.edisp_dir 165 self.edisp_path = self.edisp_dir
  166 + self.bgd_path = self.bgd_dir
131 167
132 # Create directory structure 168 # Create directory structure
133 if not os.path.isdir(self.ea_path): 169 if not os.path.isdir(self.ea_path):
@@ -136,18 +172,36 @@ class caldb(): @@ -136,18 +172,36 @@ class caldb():
136 os.makedirs(self.psf_path) 172 os.makedirs(self.psf_path)
137 if not os.path.isdir(self.edisp_path): 173 if not os.path.isdir(self.edisp_path):
138 os.makedirs(self.edisp_path) 174 os.makedirs(self.edisp_path)
  175 + if not os.path.isdir(self.bgd_path):
  176 + os.makedirs(self.bgd_path)
139 177
140 # Return 178 # Return
141 return 179 return
142 180
143 - def open(self, name, clobber=True): 181 + def open(self, version, split=False, clobber=True):
144 """ 182 """
145 - Open existing or create new calibration. 183 + Open existing or create new calibration. The actual version will put
  184 + all calibrations in the same file, although each part of the response
  185 + function will have its own logical name. We can thus easily modify
  186 + the script to put each calibration information in a separate file.
  187 +
  188 + Parameters:
  189 + version - Filename version
  190 + Keywords:
  191 + split - Split IRF over several files?
  192 + clobber - Overwrite existing files?
146 """ 193 """
147 # Set calibrate file names 194 # Set calibrate file names
148 - self.ea_file = "ea_"+name+".fits"  
149 - self.psf_file = "psf_"+name+".fits"  
150 - self.edisp_file = "edisp_"+name+".fits" 195 + if split:
  196 + self.ea_file = "ea_"+version+".fits"
  197 + self.psf_file = "psf_"+version+".fits"
  198 + self.edisp_file = "edisp_"+version+".fits"
  199 + self.bgd_file = "bgd_"+version+".fits"
  200 + else:
  201 + self.ea_file = "irf_"+version+".fits"
  202 + self.psf_file = "irf_"+version+".fits"
  203 + self.edisp_file = "irf_"+version+".fits"
  204 + self.bgd_file = "irf_"+version+".fits"
151 205
152 # Open calibration database index 206 # Open calibration database index
153 self.cif = GFits(self.base_path+"/caldb.indx", True) 207 self.cif = GFits(self.base_path+"/caldb.indx", True)
@@ -159,17 +213,42 @@ class caldb(): @@ -159,17 +213,42 @@ class caldb():
159 self.create_cif() 213 self.create_cif()
160 self.hdu_cif = self.cif.table("CIF") 214 self.hdu_cif = self.cif.table("CIF")
161 215
162 - # Open calibration files  
163 - self.open_ea(self.ea_path+"/"+self.ea_file)  
164 - self.open_psf(self.psf_path+"/"+self.psf_file)  
165 - self.open_edisp(self.edisp_path+"/"+self.edisp_file) 216 + # Set filenames
  217 + ea_filename = self.ea_path+"/"+self.ea_file
  218 + psf_filename = self.psf_path+"/"+self.psf_file
  219 + edisp_filename = self.edisp_path+"/"+self.edisp_file
  220 + bgd_filename = self.bgd_path+"/"+self.bgd_file
  221 +
  222 + # Open files
  223 + if split:
  224 + self.ea = GFits(ea_filename, True)
  225 + self.psf = GFits(psf_filename, True)
  226 + self.edisp = GFits(edisp_filename, True)
  227 + self.bgd = GFits(bgd_filename, True)
  228 + else:
  229 + self.irf = GFits(ea_filename, True)
  230 + self.ea = self.irf
  231 + self.psf = self.irf
  232 + self.edisp = self.irf
  233 + self.bgd = self.irf
  234 +
  235 + # Open HDUs
  236 + self.open_ea()
  237 + self.open_psf()
  238 + self.open_edisp()
  239 + self.open_bgd()
166 240
167 # Return 241 # Return
168 return 242 return
169 243
170 def close(self): 244 def close(self):
171 """ 245 """
172 - Close calibration. 246 + Close calibration FITS files.
  247 +
  248 + Parameters:
  249 + None
  250 + Keywords:
  251 + None
173 """ 252 """
174 # Add information to CIF. We do this now as before we did not 253 # Add information to CIF. We do this now as before we did not
175 # necessarily have all the information at hand (in particular 254 # necessarily have all the information at hand (in particular
@@ -181,31 +260,57 @@ class caldb(): @@ -181,31 +260,57 @@ class caldb():
181 self.cif.save(True) 260 self.cif.save(True)
182 self.cif.close() 261 self.cif.close()
183 self.cif = None 262 self.cif = None
184 -  
185 - # Close effective area  
186 - if self.ea != None:  
187 - self.ea.save(True)  
188 - self.ea.close()  
189 - self.ea = None  
190 -  
191 - # Close point spread function  
192 - if self.psf != None:  
193 - self.psf.save(True)  
194 - self.psf.close()  
195 - self.psf = None  
196 -  
197 - # Close energy dispersion  
198 - if self.edisp != None:  
199 - self.edisp.save(True)  
200 - self.edisp.close() 263 +
  264 + # If IRF file exists then close it now. All file pointers will
  265 + # be set to None
  266 + if self.irf != None:
  267 + self.irf.save(True)
  268 + self.irf.close()
  269 + self.irf = None
  270 + self.ea = None
  271 + self.psf = None
201 self.edisp = None 272 self.edisp = None
  273 + self.bgd = None
  274 +
  275 + # ... otherwise we have split files, so we have to close them
  276 + # all
  277 + else:
  278 +
  279 + # Close effective area file
  280 + if self.ea != None:
  281 + self.ea.save(True)
  282 + self.ea.close()
  283 + self.ea = None
  284 +
  285 + # Close point spread function file
  286 + if self.psf != None:
  287 + self.psf.save(True)
  288 + self.psf.close()
  289 + self.psf = None
  290 +
  291 + # Close energy dispersion file
  292 + if self.edisp != None:
  293 + self.edisp.save(True)
  294 + self.edisp.close()
  295 + self.edisp = None
  296 +
  297 + # Close background file
  298 + if self.bgd != None:
  299 + self.bgd.save(True)
  300 + self.bgd.close()
  301 + self.bgd = None
202 302
203 # Return 303 # Return
204 return 304 return
205 305
206 def create_cif(self): 306 def create_cif(self):
207 """ 307 """
208 - Create CIF extension in FITS file. 308 + Create Calibration Index File (CIF) extension in CIF FITS file.
  309 +
  310 + Parameters:
  311 + None
  312 + Keywords:
  313 + None
209 """ 314 """
210 # Create binary table 315 # Create binary table
211 table = GFitsBinTable() 316 table = GFitsBinTable()
@@ -248,6 +353,11 @@ class caldb(): @@ -248,6 +353,11 @@ class caldb():
248 def add_cif_info(self): 353 def add_cif_info(self):
249 """ 354 """
250 Add information to CIF extension. 355 Add information to CIF extension.
  356 +
  357 + Parameters:
  358 + None
  359 + Keywords:
  360 + None
251 """ 361 """
252 # Append 3 rows to CIF extension 362 # Append 3 rows to CIF extension
253 self.hdu_cif.append_rows(3) 363 self.hdu_cif.append_rows(3)
@@ -310,20 +420,35 @@ class caldb(): @@ -310,20 +420,35 @@ class caldb():
310 self.hdu_cif["CAL_CBD"][row, n] = "NONE" 420 self.hdu_cif["CAL_CBD"][row, n] = "NONE"
311 else: 421 else:
312 self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n] 422 self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n]
  423 +
  424 + # Add background information
  425 + row = self.hdu_cif.nrows()-1
  426 + self.hdu_cif["CAL_DIR"][row] = self.bgd_dir
  427 + self.hdu_cif["CAL_FILE"][row] = self.bgd_file
  428 + self.hdu_cif["CAL_CNAM"][row] = self.bgd_name
  429 + self.hdu_cif["CAL_DESC"][row] = self.bgd_desc
  430 + self.hdu_cif["CAL_XNO"][row] = 1
  431 + for n in range(9):
  432 + if n >= len(self.bgd_bounds):
  433 + self.hdu_cif["CAL_CBD"][row, n] = "NONE"
  434 + else:
  435 + self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
313 436
314 # Return 437 # Return
315 return 438 return
316 439
317 - def open_ea(self, filename): 440 + def open_ea(self):
318 """ 441 """
319 Open effective area extension. If no extension was found 442 Open effective area extension. If no extension was found
320 then create one. We do not yet set any data as we don't 443 then create one. We do not yet set any data as we don't
321 know the dimensions of the table yet. 444 know the dimensions of the table yet.
322 - """  
323 - # Open FITS file  
324 - self.ea = GFits(filename, True)  
325 445
326 - # Get extensions. If they do not exist then create them now 446 + Parameters:
  447 + None
  448 + Keywords:
  449 + None
  450 + """
  451 + # Get extension. If it does not exist then create it now.
327 try: 452 try:
328 self.hdu_ea = self.ea.table("EFFECTIVE AREA") 453 self.hdu_ea = self.ea.table("EFFECTIVE AREA")
329 except: 454 except:
@@ -346,16 +471,18 @@ class caldb(): @@ -346,16 +471,18 @@ class caldb():
346 # Return 471 # Return
347 return 472 return
348 473
349 - def open_psf(self, filename): 474 + def open_psf(self):
350 """ 475 """
351 Open point spread function extension. If no extension was 476 Open point spread function extension. If no extension was
352 found then create one. We do not yet set any data as we don't 477 found then create one. We do not yet set any data as we don't
353 know the dimensions of the table yet. 478 know the dimensions of the table yet.
354 - """  
355 - # Open FITS file  
356 - self.psf = GFits(filename, True)  
357 479
358 - # Get extensions. If they do not exist then create them now 480 + Parameters:
  481 + None
  482 + Keywords:
  483 + None
  484 + """
  485 + # Get extension. If it does not exist then create it now.
359 try: 486 try:
360 self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION") 487 self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION")
361 except: 488 except:
@@ -378,16 +505,18 @@ class caldb(): @@ -378,16 +505,18 @@ class caldb():
378 # Return 505 # Return
379 return 506 return
380 507
381 - def open_edisp(self, filename): 508 + def open_edisp(self):
382 """ 509 """
383 Open energy dispersion extension. If no extension was found 510 Open energy dispersion extension. If no extension was found
384 then create one. We do not yet set any data as we don't 511 then create one. We do not yet set any data as we don't
385 know the dimensions of the table yet. 512 know the dimensions of the table yet.
386 - """  
387 - # Open FITS file  
388 - self.edisp = GFits(filename, True)  
389 513
390 - # Get extensions. If they do not exist then create them now 514 + Parameters:
  515 + None
  516 + Keywords:
  517 + None
  518 + """
  519 + # Get extension. If it does not exist then create it now.
391 try: 520 try:
392 self.hdu_edisp = self.edisp.table("ENERGY DISPERSION") 521 self.hdu_edisp = self.edisp.table("ENERGY DISPERSION")
393 except: 522 except:
@@ -410,9 +539,50 @@ class caldb(): @@ -410,9 +539,50 @@ class caldb():
410 # Return 539 # Return
411 return 540 return
412 541
  542 + def open_bgd(self):
  543 + """
  544 + Open background extension. If no extension was found
  545 + then create one. We do not yet set any data as we don't
  546 + know the dimensions of the table yet.
  547 +
  548 + Parameters:
  549 + None
  550 + Keywords:
  551 + None
  552 + """
  553 + # Get extension. If it does not exist then create it now.
  554 + try:
  555 + self.hdu_bgd = self.bgd.table("BACKGROUND")
  556 + except:
  557 + # Create binary table
  558 + table = GFitsBinTable()
  559 +
  560 + # Set extension name
  561 + table.extname("BACKGROUND")
  562 +
  563 + # Set OGIP keywords
  564 + self.set_ogip_keywords(table, self.bgd_doc, \
  565 + ["RESPONSE", self.bgd_name])
  566 +
  567 + # Append table
  568 + self.bgd.append(table)
  569 +
  570 + # Get extension
  571 + self.hdu_bgd = self.bgd.table("BACKGROUND")
  572 +
  573 + # Return
  574 + return
  575 +
413 def set_ogip_keywords(self, hdu, hdudoc, hduclas): 576 def set_ogip_keywords(self, hdu, hdudoc, hduclas):
414 """ 577 """
415 Set standard OGIP keywords for extension. 578 Set standard OGIP keywords for extension.
  579 +
  580 + Parameters:
  581 + hdu - Header Data Unit.
  582 + hdudoc - Documentation reference string
  583 + hduclas - Array of HDUCLAS fields
  584 + Keywords:
  585 + None
416 """ 586 """
417 # Set keywords 587 # Set keywords
418 hdu.card("ORIGIN", "IRAP", "Name of organization making this file") 588 hdu.card("ORIGIN", "IRAP", "Name of organization making this file")
@@ -433,12 +603,17 @@ class caldb(): @@ -433,12 +603,17 @@ class caldb():
433 def make_2D(self, array, hdu, name, unit, scale=1.0): 603 def make_2D(self, array, hdu, name, unit, scale=1.0):
434 """ 604 """
435 Make 2D matrix as function of energy and offset angle from a 605 Make 2D matrix as function of energy and offset angle from a
436 - ROOT 2D histogram. 606 + ROOT 2D histogram. If the HDU has already the energy and
  607 + offset angle columns, this method will simply add another data
  608 + column.
  609 + If name==None, the method will not append any data column.
437 610
438 Parameters: 611 Parameters:
439 - array - ROOT 2D histgram. 612 + array - ROOT 2D histogram.
440 hdu - FITS HDU. 613 hdu - FITS HDU.
441 - Keyword: 614 + name - Data column name.
  615 + unit - Data unit.
  616 + Keywords:
442 scale - Scaling factor for histogram values. 617 scale - Scaling factor for histogram values.
443 """ 618 """
444 # Extract energy and offset angle vectors 619 # Extract energy and offset angle vectors
@@ -485,7 +660,7 @@ class caldb(): @@ -485,7 +660,7 @@ class caldb():
485 hdu["THETA_HI"][0,ioff] = o_hi 660 hdu["THETA_HI"][0,ioff] = o_hi
486 # 661 #
487 # "NAME" 662 # "NAME"
488 - if not hdu.hascolumn(name): 663 + if name != None and not hdu.hascolumn(name):
489 hdu.append_column(GFitsTableFloatCol(name, 1, neng*noffset)) 664 hdu.append_column(GFitsTableFloatCol(name, 1, neng*noffset))
490 hdu[name].unit(unit) 665 hdu[name].unit(unit)
491 hdu[name].dim([neng, noffset]) 666 hdu[name].dim([neng, noffset])
@@ -515,44 +690,217 @@ class caldb(): @@ -515,44 +690,217 @@ class caldb():
515 Parameters: 690 Parameters:
516 filename - ROOT 2D performance filename. 691 filename - ROOT 2D performance filename.
517 """ 692 """
518 - # Open performance file 693 + # Open ROOT performance file
519 file = TFile(filename) 694 file = TFile(filename)
520 695
521 - # Allocate ROOT 2D array  
522 - array = TH2F() 696 + # Create effective area
  697 + self.root2ea(file)
  698 +
  699 + # Create point spread function
  700 + self.root2psf(file)
  701 +
  702 + # Create energy dispersion
  703 + self.root2edisp(file)
  704 +
  705 + # Create background
  706 + self.root2bgd(file)
  707 +
  708 + # Return
  709 + return
523 710
524 - # Make 2D effective 80% area 711 + def root2ea(self, file):
  712 + """
  713 + Translate ROOT to CALDB effective area extension. The following ROOT
  714 + histograms are used:
  715 +
  716 + EffectiveAreaEtrue_offaxis -> EFFAREA
  717 + EffectiveArea_offaxis -> EFFAREA_RECO
  718 +
  719 + Parameters:
  720 + file - ROOT file.
  721 + Keywords:
  722 + None
  723 + """
  724 + # Continue only if effective area HDU has been opened
525 if self.hdu_ea != None: 725 if self.hdu_ea != None:
526 - file.GetObject("EffectiveArea80_offaxis", array)  
527 - bounds = self.make_2D(array, self.hdu_ea, "EFFAREA", "m2", scale=1.0/0.8) 726 +
  727 + # Allocate ROOT 2D arrays
  728 + etrue = TH2F()
  729 + ereco = TH2F()
  730 + file.GetObject("EffectiveAreaEtrue_offaxis", etrue)
  731 + file.GetObject("EffectiveArea_offaxis", ereco)
  732 + #print etrue.GetXaxis().GetBinLowEdge(1), etrue.GetXaxis().GetBinUpEdge(etrue.GetXaxis().GetNbins())
  733 + #print ereco.GetXaxis().GetBinLowEdge(1), ereco.GetXaxis().GetBinUpEdge(ereco.GetXaxis().GetNbins())
  734 +
  735 + # Rebin etrue histogram
  736 + etrue.RebinX(10)
  737 + neng = etrue.GetXaxis().GetNbins()
  738 + noffset = etrue.GetYaxis().GetNbins()
  739 + for ioff in range(noffset):
  740 + for ieng in range(neng):
  741 + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
  742 + etrue.SetBinContent(ieng+1,ioff+1,value)
  743 +
  744 + # Write boundaries (use Ereco boundaries)
  745 + bounds = self.make_2D(ereco, self.hdu_ea, None, "m2")
528 for b in bounds: 746 for b in bounds:
529 self.ea_bounds.append(b) 747 self.ea_bounds.append(b)
530 self.set_cif_keywords(self.hdu_ea, self.ea_name, \ 748 self.set_cif_keywords(self.hdu_ea, self.ea_name, \
531 self.ea_bounds, self.ea_desc) 749 self.ea_bounds, self.ea_desc)
532 750
533 - # Make 2D 68% point spread function 751 + # EFFAREA
  752 + self.make_2D(etrue, self.hdu_ea, "EFFAREA", "m2")
  753 +
  754 + # EFFAREA_RECO
  755 + self.make_2D(ereco, self.hdu_ea, "EFFAREA_RECO", "m2")
  756 +
  757 + # Return
  758 + return
  759 +
  760 + def root2psf(self, file):
  761 + """
  762 + Translate ROOT to CALDB point spread function extension. The following
  763 + ROOT histograms are used:
  764 +
  765 + 1/(2*pi*SIGMA_1) -> SCALE
  766 + AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8)
  767 + 0.0 -> AMPL_2
  768 + 0.0 -> SIGMA_2
  769 + 0.0 -> AMPL_3
  770 + 0.0 -> SIGMA_3
  771 +
  772 + Parameters:
  773 + file - ROOT file.
  774 + Keywords:
  775 + None
  776 + """
  777 + # Continue only if point spread function HDU has been opened
534 if self.hdu_psf != None: 778 if self.hdu_psf != None:
535 - file.GetObject("AngRes_offaxis", array)  
536 - bounds = self.make_2D(array, self.hdu_psf, "R68", "deg") 779 +
  780 + # Allocate ROOT 2D array
  781 + r68 = TH2F()
  782 +
  783 + # Read 68% containment histogram
  784 + file.GetObject("AngRes_offaxis", r68)
  785 + neng = r68.GetXaxis().GetNbins()
  786 + noffset = r68.GetYaxis().GetNbins()
  787 +
  788 + # Compute scale histogram
  789 + scale = r68.Clone()
  790 + for ioff in range(noffset):
  791 + for ieng in range(neng):
  792 + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1)
  793 + if integral > 0.0:
  794 + value = 1.0 / integral
  795 + else:
  796 + value = 0.0
  797 + scale.SetBinContent(ieng+1,ioff+1,value)
  798 +
  799 + # Set zero histogram
  800 + zero = r68.Clone()
  801 + for ioff in range(noffset):
  802 + for ieng in range(neng):
  803 + zero.SetBinContent(ieng+1,ioff+1,0.0)
  804 +
  805 + # Set boundaries
  806 + bounds = self.make_2D(r68, self.hdu_psf, None, "deg")
537 for b in bounds: 807 for b in bounds:
538 self.psf_bounds.append(b) 808 self.psf_bounds.append(b)
539 self.set_cif_keywords(self.hdu_psf, self.psf_name, \ 809 self.set_cif_keywords(self.hdu_psf, self.psf_name, \
540 self.psf_bounds, self.psf_desc) 810 self.psf_bounds, self.psf_desc)
541 811
542 - # Make 2D 80% point spread function  
543 - if self.hdu_psf != None:  
544 - file.GetObject("AngRes80_offaxis", array)  
545 - self.make_2D(array, self.hdu_psf, "R80", "deg") 812 + # SCALE
  813 + self.make_2D(scale, self.hdu_psf, "SCALE", "")
  814 +
  815 + # SIGMA_1
  816 + self.make_2D(r68, self.hdu_psf, "SIGMA_1", "deg")
546 817
547 - # Make 2D energy resolution 818 + # AMPL_2
  819 + self.make_2D(zero, self.hdu_psf, "AMPL_2", "")
  820 +
  821 + # SIGMA_2
  822 + self.make_2D(zero, self.hdu_psf, "SIGMA_2", "deg")
  823 +
  824 + # AMPL_3
  825 + self.make_2D(zero, self.hdu_psf, "AMPL_3", "")
  826 +
  827 + # SIGMA_3
  828 + self.make_2D(zero, self.hdu_psf, "SIGMA_3", "deg")
  829 +
  830 + # Return
  831 + return
  832 +
  833 + def root2edisp(self, file):
  834 + """
  835 + Translate ROOT to CALDB energy dispersion extension. The following ROOT
  836 + histograms are used:
  837 +
  838 + ERes_offaxis -> ERES
  839 + Ebias_offaxis -> EBIAS
  840 +
  841 + Parameters:
  842 + file - ROOT file.
  843 + Keywords:
  844 + None
  845 + """
  846 + # Continue only if energy dispersion HDU has been opened
548 if self.hdu_edisp != None: 847 if self.hdu_edisp != None:
549 - file.GetObject("ERes", array)  
550 - bounds = self.make_2D(array, self.hdu_edisp, "ERES68", "") 848 +
  849 + # Allocate ROOT 2D array
  850 + eres = TH2F()
  851 + ebias = TH2F()
  852 + file.GetObject("ERes_offaxis", eres)
  853 + file.GetObject("Ebias_offaxis", ebias)
  854 +
  855 + # Set boundaries
  856 + bounds = self.make_2D(eres, self.hdu_edisp, None, "deg")
551 for b in bounds: 857 for b in bounds:
552 self.edisp_bounds.append(b) 858 self.edisp_bounds.append(b)
553 self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \ 859 self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \
554 self.edisp_bounds, self.edisp_desc) 860 self.edisp_bounds, self.edisp_desc)
555 - 861 +
  862 + # ERES
  863 + self.make_2D(eres, self.hdu_edisp, "ERES", "Ereco/Etrue")
  864 +
  865 + # EBIAS
  866 + self.make_2D(ebias, self.hdu_edisp, "EBIAS", "Ereco/Etrue")
  867 +
  868 + # Return
  869 + return
  870 +
  871 + def root2bgd(self, file):
  872 + """
  873 + Translate ROOT to CALDB background extension. The following ROOT
  874 + histograms are used:
  875 +
  876 + BGRatePerSqDeg_offaxis -> BGD
  877 + BGRatePerSqDeg_offaxis -> BGD_RECO
  878 +
  879 + Parameters:
  880 + file - ROOT file.
  881 + Keywords:
  882 + None
  883 + """
  884 + # Continue only if background HDU has been opened
  885 + if self.hdu_bgd != None:
  886 +
  887 + # Allocate ROOT 2D array
  888 + array = TH2F()
  889 + file.GetObject("BGRatePerSqDeg_offaxis", array)
  890 +
  891 + # Set boundaries
  892 + bounds = self.make_2D(array, self.hdu_bgd, None, "deg")
  893 + for b in bounds:
  894 + self.bgd_bounds.append(b)
  895 + self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \
  896 + self.bgd_bounds, self.bgd_desc)
  897 +
  898 + # BGD
  899 + self.make_2D(array, self.hdu_bgd, "BGD", "")
  900 +
  901 + # BGD_RECO
  902 + self.make_2D(array, self.hdu_bgd, "BGD_RECO", "")
  903 +
556 # Return 904 # Return
557 return 905 return
558 906
@@ -589,22 +937,41 @@ class caldb(): @@ -589,22 +937,41 @@ class caldb():
589 #==========================# 937 #==========================#
590 if __name__ == '__main__': 938 if __name__ == '__main__':
591 """ 939 """
592 - This script translates a ROOT 2D performance file in a CALDB compliant  
593 - CTA response. 940 + This script translates a CTA ROOT 2D performance files in a CALDB compliant
  941 + CTA response. The script will create the response for 3 observations to
  942 + illustrate the run-wise file structure.
  943 +
  944 + The CTA ROOT 2D performance files contain the following histograms:
  945 + - DiffSens_offaxis (2D) - Differential sensitivity
  946 + - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?)
  947 + - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment
  948 + - AngRes_offaxis (2D) - Angular resolution 68% containment
  949 + - AngRes80_offaxis (2D) - Angular resolution 80% containment
  950 + - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree
  951 + - BGRate_offaxis (2D) - Background rate
  952 + - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins)
  953 + - MigMatrix_offaxis (3D)
  954 + - EestOverEtrue_offaxis (3D)
  955 + - ERes_offaxis (2D) - Energy Resolution
  956 + - Ebias_offaxis (2D) - Energy bias
594 """ 957 """
595 # Set ROOT filename 958 # Set ROOT filename
596 - #path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxis"  
597 - #irf = "SubarrayE_offaxis.root"  
598 - path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_Nov2011"  
599 - irf = "SubarrayE_IFAE_50hours_20111121_offaxis.root" 959 + path = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEOffaxisPerformanceBEI_May2012"
  960 + irf = "SubarrayE_IFAE_50hours_20120510_offaxis.root"
600 filename = path+"/"+irf 961 filename = path+"/"+irf
601 962
602 - # Allocate caldb  
603 - db = caldb("E") 963 + # Set observation identifiers
  964 + obsids = ["000001", "000002", "000003"]
  965 +
  966 + # Generate CALDB entries, one for each observation (or run)
  967 + for obsid in obsids:
  968 +
  969 + # Allocate caldb
  970 + db = caldb("E", obsid)
604 971
605 - # Open calibration  
606 - db.open("test") 972 + # Open calibration
  973 + db.open("test")
607 974
608 - # Translate ROOT to CALDB information  
609 - db.root2caldb(filename) 975 + # Translate ROOT to CALDB information
  976 + db.root2caldb(filename)
610 977
611 \ No newline at end of file 978 \ No newline at end of file