Domenico Tiziani / ctools

Commit a266b5659e5c9469de5450b7bc21bc26e2e8c14b

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

Revise csroot2caldb script (#1828)

Inline Side-by-side
Showing 2 changed files with 723 additions and 758 deletions
cscripts/csroot2caldb.par
  View file @a266b56
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: csroot2caldb.par #
8   -# Version: 1.0.0 #
  8 +# Version: 1.2.0 #
9 9 # #
10 10 # Author: Juergen Knoedlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -24,6 +24,7 @@ id, s, a, South_50h,,, "Instrument Response Function identifier"
24 24 version, s, a, prod3,,, "Instrument Response Function version"
25 25 analysis, s, a, DESY,,, "Instrument Response Function analysis"
26 26 psftype, s, h, Gauss,Gauss|King,, "Point spread function type"
  27 +split, b, h, no,,, "Split IRF components into different files"
27 28 norm1d, b, h, no,,, "Normalize on 1D histograms"
28 29 rebin, b, h, no,,, "Rebin"
29 30 eascale, r, h, 1.0,,, "Scaling factor for effective areas"
... ...
cscripts/csroot2caldb.py
  View file @a266b56
... ... @@ -39,7 +39,7 @@ except:
39 39 # ================== #
40 40 class csroot2caldb(ctools.cscript):
41 41 """
42   - Generates IRFs in CALDB from a ROOT offaxis performance file.
  42 + Generates IRFs in CALDB from a ROOT offaxis performance file
43 43 """
44 44  
45 45 # Constructor
... ... @@ -51,9 +51,8 @@ class csroot2caldb(ctools.cscript):
51 51 self._name = 'csroot2caldb'
52 52 self._version = '1.2.0'
53 53  
54   - # Initialise some members
  54 + # Initialise file members
55 55 self._cif = None
56   - self._irf = None
57 56 self._ea = None
58 57 self._psf = None
59 58 self._edisp = None
... ... @@ -80,6 +79,7 @@ class csroot2caldb(ctools.cscript):
80 79 self['analysis'].string()
81 80 self['version'].string()
82 81 self['psftype'].string()
  82 + self['split'].boolean()
83 83 self['norm1d'].boolean()
84 84 self['rebin'].boolean()
85 85 self['eascale'].real()
... ... @@ -90,64 +90,66 @@ class csroot2caldb(ctools.cscript):
90 90 # Write input parameters into logger
91 91 if self._logTerse():
92 92 self._log_parameters()
93   - self._log("\n")
  93 + self._log('\n')
94 94  
95 95 # Return
96 96 return
97 97  
98   - def _initialise(self):
99   -
100   - # Set CALDB information
101   - self._cal_tel = "CTA"
102   - self._cal_inst = self["inst"].string().upper()
103   - self._cal_obsid = self["id"].string()
104   - self._cal_det = "NONE"
105   - self._cal_flt = "NONE"
106   - self._cal_class = "BCF"
107   - self._cal_type = "DATA"
108   - self._cal_qual = 0 # 0=good, 1=bad, 2=dubious, ...
109   - self._cal_date = "14/01/30"
110   - self._val_date = "2014-01-30"
111   - self._val_time = "00:00:00"
112   - self._ref_time = 51544.0
113   - self._cal_name = "NAME("+self["id"].string()+")"
114   - self._cal_version = "VERSION("+self["version"].string()+")"
115   - self._cal_cut = "CLASS(BEST)"
116   - self._cal_analysis = "ANALYSIS("+self["analysis"].string()+")"
117   - #
118   - self._ea_name = "EFF_AREA"
119   - self._ea_doc = "CAL/GEN/92-019"
120   - self._ea_bounds = [self._cal_name, self._cal_version,
121   - self._cal_cut, self._cal_analysis]
122   - self._ea_desc = "CTA effective area"
123   - #
124   - self._psf_name = "RPSF"
125   - self._psf_doc = "CAL/GEN/92-020"
126   - self._psf_bounds = [self._cal_name, self._cal_version,
127   - self._cal_cut, self._cal_analysis]
128   - self._psf_desc = "CTA point spread function"
129   - #
130   - self._edisp_name = "EDISP"
131   - self._edisp_doc = "???"
132   - self._edisp_bounds = [self._cal_name, self._cal_version,
133   - self._cal_cut, self._cal_analysis]
134   - self._edisp_desc = "CTA energy dispersion"
135   - #
136   - self._bgd_name = "BGD"
137   - self._bgd_doc = "???"
138   - self._bgd_bounds = [self._cal_name, self._cal_version,
139   - self._cal_cut, self._cal_analysis]
140   - self._bgd_desc = "CTA background"
141   -
142   - # Create directory structure
143   - self._make_dirs()
144   -
145   - # Return
146   - return
  98 + def _init_metadata(self):
  99 + """
  100 + Initialise dictionary containing the IRF metadata
147 101  
148   - def _make_dirs(self):
  102 + Returns
  103 + -------
  104 + irf : dict
  105 + IRF metadata dictionary
149 106 """
150   - Generate CALDB directory structure for one observation identifier.
  107 + # Initialise IRF metadata dictionary
  108 + irf = {}
  109 +
  110 + # Set calibration information
  111 + cal_name = 'NAME('+self['id'].string()+')'
  112 + cal_version = 'VERSION('+self['version'].string()+')'
  113 + cal_cut = 'CLASS(BEST)'
  114 + cal_analysis = 'ANALYSIS('+self['analysis'].string()+')'
  115 +
  116 + # Set IRF information
  117 + irf['CAL_TEL'] = 'CTA'
  118 + irf['CAL_INST'] = self['inst'].string().upper()
  119 + irf['CAL_OBSID'] = self['id'].string()
  120 + irf['CAL_DET'] = 'NONE'
  121 + irf['CAL_FLT'] = 'NONE'
  122 + irf['CAL_CLASS'] = 'BCF'
  123 + irf['CAL_TYPE'] = 'DATA'
  124 + irf['CAL_QUAL'] = 0 # 0=good, 1=bad, 2=dubious, ...
  125 + irf['CAL_DATE'] = '14/01/30'
  126 + irf['VAL_DATE'] = '2014-01-30'
  127 + irf['VAL_TIME'] = '00:00:00'
  128 + irf['REF_TIME'] = 51544.0
  129 + irf['EA_NAME'] = 'EFF_AREA'
  130 + irf['EA_DOC'] = 'CAL/GEN/92-019'
  131 + irf['EA_BOUNDS'] = [cal_name, cal_version, cal_cut, cal_analysis]
  132 + irf['EA_DESC'] = 'CTA effective area'
  133 + irf['PSF_NAME'] = 'RPSF'
  134 + irf['PSF_DOC'] = 'CAL/GEN/92-020'
  135 + irf['PSF_BOUNDS'] = [cal_name, cal_version, cal_cut, cal_analysis]
  136 + irf['PSF_DESC'] = 'CTA point spread function'
  137 + irf['EDISP_NAME'] = 'EDISP'
  138 + irf['EDISP_DOC'] = '???'
  139 + irf['EDISP_BOUNDS'] = [cal_name, cal_version, cal_cut, cal_analysis]
  140 + irf['EDISP_DESC'] = 'CTA energy dispersion'
  141 + irf['BGD_NAME'] = 'BGD'
  142 + irf['BGD_DOC'] = '???'
  143 + irf['BGD_BOUNDS'] = [cal_name, cal_version, cal_cut, cal_analysis]
  144 + irf['BGD_DESC'] = 'CTA background'
  145 +
  146 + # Return metadata
  147 + return irf
  148 +
  149 + def _make_dirs(self, version, irf):
  150 + """
  151 + Generate CALDB directory structure for one observation identifier
  152 +
151 153 The structure is given by
152 154  
153 155 data/<tel>/<inst>/bcf/<obsid>
... ... @@ -155,172 +157,193 @@ class csroot2caldb(ctools.cscript):
155 157 where <tel> is "cta" and <inst> is the instrument specified in the
156 158 CALDB constructor (the instrument may be used for different array
157 159 configurations).
  160 +
  161 + Parameters
  162 + ----------
  163 + version : str
  164 + Version string
  165 + irf : dict
  166 + IRF metadata dictionary
  167 +
  168 + Returns
  169 + -------
  170 + dir : dict
  171 + Directory structure information
158 172 """
  173 + # Initialise directory structure dictionary
  174 + dir = {}
  175 +
159 176 # Set base directory
160   - base_dir = "data"
161   - base_dir += "/"+self._cal_tel.lower()
162   - base_dir += "/"+self._cal_inst.lower()
  177 + base_dir = 'data'
  178 + base_dir += '/'+irf['CAL_TEL'].lower()
  179 + base_dir += '/'+irf['CAL_INST'].lower()
163 180  
164 181 # Set directory names
165   - obs_dir = base_dir+"/bcf/"+self._cal_obsid
166   - self._ea_dir = obs_dir
167   - self._psf_dir = obs_dir
168   - self._edisp_dir = obs_dir
169   - self._bgd_dir = obs_dir
170   -
171   - # Optionally prefix path
172   - outdir = self["outdir"].string()
  182 + obs_dir = base_dir+'/bcf/'+irf['CAL_OBSID']
  183 + dir['EA_DIR'] = obs_dir
  184 + dir['PSF_DIR'] = obs_dir
  185 + dir['EDISP_DIR'] = obs_dir
  186 + dir['BGD_DIR'] = obs_dir
  187 +
  188 + # Set path to response components. If the "outdir" parameter is set
  189 + # then prefix the value to the path of the response components
  190 + outdir = self['outdir'].string()
173 191 if len(outdir) > 0:
174   - self._base_path = outdir+"/"+base_dir
175   - self._ea_path = outdir+"/"+self._ea_dir
176   - self._psf_path = outdir+"/"+self._psf_dir
177   - self._edisp_path = outdir+"/"+self._edisp_dir
178   - self._bgd_path = outdir+"/"+self._bgd_dir
  192 + dir['BASE_PATH'] = outdir+'/'+base_dir
  193 + dir['EA_PATH'] = outdir+'/'+dir['EA_DIR']
  194 + dir['PSF_PATH'] = outdir+'/'+dir['PSF_DIR']
  195 + dir['EDISP_PATH'] = outdir+'/'+dir['EDISP_DIR']
  196 + dir['BGD_PATH'] = outdir+'/'+dir['BGD_DIR']
  197 + else:
  198 + dir['BASE_PATH'] = base_dir
  199 + dir['EA_PATH'] = dir['EA_DIR']
  200 + dir['PSF_PATH'] = dir['PSF_DIR']
  201 + dir['EDISP_PATH'] = dir['EDISP_DIR']
  202 + dir['BGD_PATH'] = dir['BGD_DIR']
  203 +
  204 + # Set IRF component file names. If the "split" parameter is "yes" then
  205 + # split the IRF components over several files.
  206 + if self['split'].boolean():
  207 + dir['EA_FILE'] = 'ea_'+version+'.fits'
  208 + dir['PSF_FILE'] = 'psf_'+version+'.fits'
  209 + dir['EDISP_FILE'] = 'edisp_'+version+'.fits'
  210 + dir['BGD_FILE'] = 'bgd_'+version+'.fits'
179 211 else:
180   - self._base_path = base_dir
181   - self._ea_path = self._ea_dir
182   - self._psf_path = self._psf_dir
183   - self._edisp_path = self._edisp_dir
184   - self._bgd_path = self._bgd_dir
  212 + dir['EA_FILE'] = 'irf_'+version+'.fits'
  213 + dir['PSF_FILE'] = 'irf_'+version+'.fits'
  214 + dir['EDISP_FILE'] = 'irf_'+version+'.fits'
  215 + dir['BGD_FILE'] = 'irf_'+version+'.fits'
185 216  
186 217 # Create directory structure
187   - if not os.path.isdir(self._ea_path):
188   - os.makedirs(self._ea_path)
189   - if not os.path.isdir(self._psf_path):
190   - os.makedirs(self._psf_path)
191   - if not os.path.isdir(self._edisp_path):
192   - os.makedirs(self._edisp_path)
193   - if not os.path.isdir(self._bgd_path):
194   - os.makedirs(self._bgd_path)
195   -
196   - # Return
197   - return
198   -
199   - def _open(self, version, split=False, clobber=True):
  218 + if not os.path.isdir(dir['EA_PATH']):
  219 + os.makedirs(dir['EA_PATH'])
  220 + if not os.path.isdir(dir['PSF_PATH']):
  221 + os.makedirs(dir['PSF_PATH'])
  222 + if not os.path.isdir(dir['EDISP_PATH']):
  223 + os.makedirs(dir['EDISP_PATH'])
  224 + if not os.path.isdir(dir['BGD_PATH']):
  225 + os.makedirs(dir['BGD_PATH'])
  226 +
  227 + # Return dir
  228 + return dir
  229 +
  230 + def _open(self, irf, dir):
200 231 """
201   - Open existing or create new calibration. The actual version will put
202   - all calibrations in the same file, although each part of the response
203   - function will have its own logical name. We can thus easily modify
204   - the script to put each calibration information in a separate file.
205   -
206   - Parameters:
207   - version - Filename version
208   - Keywords:
209   - split - Split IRF over several files?
210   - clobber - Overwrite existing files?
211   - """
212   - # Set calibrate file names
213   - if split:
214   - self._ea_file = "ea_"+version+".fits"
215   - self._psf_file = "psf_"+version+".fits"
216   - self._edisp_file = "edisp_"+version+".fits"
217   - self._bgd_file = "bgd_"+version+".fits"
218   - else:
219   - self._ea_file = "irf_"+version+".fits"
220   - self._psf_file = "irf_"+version+".fits"
221   - self._edisp_file = "irf_"+version+".fits"
222   - self._bgd_file = "irf_"+version+".fits"
  232 + Open existing or create new calibration
  233 +
  234 + The actual version will put all calibrations in the same file, although
  235 + each part of the response function will have its own logical name. We
  236 + can thus easily modify the script to put each calibration information
  237 + in a separate file
223 238  
  239 + Parameters
  240 + ----------
  241 + irf : dict
  242 + IRF metadata dictionary
  243 + dir : dir
  244 + Directory structure dictionary
  245 + """
224 246 # Open calibration database index
225   - self._cif = gammalib.GFits(self._base_path+"/caldb.indx", True)
  247 + self._cif = gammalib.GFits(dir['BASE_PATH']+'/caldb.indx', True)
226 248  
227 249 # If file has no CIF extension than create it now
228 250 try:
229   - self._hdu_cif = self._cif.table("CIF")
  251 + self._hdu_cif = self._cif.table('CIF')
230 252 except:
231   - self._create_cif()
232   - self._hdu_cif = self._cif.table("CIF")
  253 + self._cif.append(calutils.create_cif_table())
  254 + self._hdu_cif = self._cif.table('CIF')
233 255  
234   - # Set filenames
235   - ea_filename = self._ea_path+"/"+self._ea_file
236   - psf_filename = self._psf_path+"/"+self._psf_file
237   - edisp_filename = self._edisp_path+"/"+self._edisp_file
238   - bgd_filename = self._bgd_path+"/"+self._bgd_file
  256 + # Set IRF component filenames
  257 + ea_filename = dir['EA_PATH']+'/'+dir['EA_FILE']
  258 + psf_filename = dir['PSF_PATH']+'/'+dir['PSF_FILE']
  259 + edisp_filename = dir['EDISP_PATH']+'/'+dir['EDISP_FILE']
  260 + bgd_filename = dir['BGD_PATH']+'/'+dir['BGD_FILE']
239 261  
240   - # Open files
241   - if split:
  262 + # Open IRF component files
  263 + if self['split'].boolean():
242 264 self._ea = gammalib.GFits(ea_filename, True)
243 265 self._psf = gammalib.GFits(psf_filename, True)
244 266 self._edisp = gammalib.GFits(edisp_filename, True)
245 267 self._bgd = gammalib.GFits(bgd_filename, True)
246 268 else:
247   - self._irf = gammalib.GFits(ea_filename, True)
248   - self._ea = self._irf
249   - self._psf = self._irf
250   - self._edisp = self._irf
251   - self._bgd = self._irf
  269 + self._ea = gammalib.GFits(ea_filename, True)
  270 + self._psf = self._ea
  271 + self._edisp = self._ea
  272 + self._bgd = self._ea
252 273  
253 274 # Open HDUs
254 275 self._hdu_ea = self._open_hdu(self._ea, "EFFECTIVE AREA",
255   - self._ea_name, self._ea_doc)
  276 + irf['EA_NAME'], irf['EA_DOC'], irf)
256 277 self._hdu_psf = self._open_hdu(self._psf, "POINT SPREAD FUNCTION",
257   - self._psf_name, self._psf_doc)
  278 + irf['PSF_NAME'], irf['PSF_DOC'], irf)
258 279 self._hdu_edisp = self._open_hdu(self._edisp, "ENERGY DISPERSION",
259   - self._edisp_name, self._edisp_doc)
  280 + irf['EDISP_NAME'], irf['EDISP_DOC'], irf)
260 281 self._hdu_bgd = self._open_hdu(self._bgd, "BACKGROUND",
261   - self._bgd_name, self._bgd_doc)
  282 + irf['BGD_NAME'], irf['BGD_DOC'], irf)
262 283  
263 284 # Return
264 285 return
265 286  
266   - def _close(self):
267   - """
268   - Close calibration FITS files.
  287 + def _close(self, irf, dir):
269 288 """
270   - # Add information to CIF. We do this now as before we did not
271   - # necessarily have all the information at hand (in particular
272   - # about the boundaries)
273   - self._add_cif_info()
274   -
275   - # Close CIF
276   - self._close_file(self._cif)
277   -
278   - # If IRF file exists then close it now. All file pointers will
279   - # be set to None
280   - if self._irf != None:
281   - self._close_file(self._irf)
282   - self._ea = None
283   - self._psf = None
284   - self._edisp = None
285   - self._bgd = None
286   -
287   - # ... otherwise we have split files, so we have to close them all
288   - else:
289   - self._close_file(self._ea)
290   - self._close_file(self._psf)
291   - self._close_file(self._edisp)
292   - self._close_file(self._bgd)
  289 + Close calibration FITS files
293 290  
294   - # Return
295   - return
296   -
297   - def _close_file(self, file):
298   - """
299   - Close one file.
  291 + Parameters
  292 + ----------
  293 + irf : dict
  294 + IRF metadata dictionary
  295 + dir : dict
  296 + Directory structure dictionary
300 297 """
301   - # If file is open then close it
302   - if file != None:
303   - file.save(True)
304   - file.close()
305   - file = None
306   -
  298 + # Add information to CIF. We do this now as before we did not
  299 + # necessarily have all the information at hand (in particular about
  300 + # the boundaries)
  301 + self._add_cif_info(irf, dir)
  302 +
  303 + # Save and close CIF
  304 + self._cif.save(True)
  305 + self._cif.close()
  306 +
  307 + # Close all IRF components. If the files have been split all components
  308 + # need to be closed, otherwise only the effective area component needs
  309 + # to be closed as representative for all components.
  310 + self._ea.save(True)
  311 + self._ea.close()
  312 + if self['split'].boolean():
  313 + self._psf.save(True)
  314 + self._edisp.save(True)
  315 + self._bgd.save(True)
  316 + self._psf.close()
  317 + self._edisp.close()
  318 + self._bgd.close()
  319 +
307 320 # Return
308 321 return
309 322  
310   - def _open_hdu(self, fits, extname, name, doc):
  323 + def _open_hdu(self, fits, extname, name, doc, irf):
311 324 """
312 325 Open HDU
313 326  
314   - If HDU does not exist then create one.
315   -
316   - Parameters:
317   - fits: FITS file
318   - extname: Extension name
319   - name: Name string
320   - doc: Document string
  327 + Opens a FITS binary table with given "extname". If HDU does not exist
  328 + in the FITS file it will be created and appended to the FITS file.
321 329  
322   - Returns:
323   - FITS table
  330 + Parameters
  331 + ----------
  332 + fits : `~gammalib.GFits`
  333 + FITS file
  334 + extname : str
  335 + Extension name
  336 + name : str
  337 + Name string
  338 + doc : str
  339 + Document string
  340 + irf : dict
  341 + IRF metadata dictionary
  342 +
  343 + Returns
  344 + -------
  345 + table : `~gammalib.GFitsBinTable`
  346 + FITS binary table
324 347 """
325 348 # Create table if it does not yet exist
326 349 if not fits.contains(extname):
... ... @@ -332,7 +355,7 @@ class csroot2caldb(ctools.cscript):
332 355 table.extname(extname)
333 356  
334 357 # Set OGIP keywords
335   - self._set_ogip_keywords(table, doc, ["RESPONSE", name])
  358 + self._set_ogip_keywords(table, doc, ['RESPONSE', name], irf)
336 359  
337 360 # Append table to FITS file
338 361 fits.append(table)
... ... @@ -340,650 +363,593 @@ class csroot2caldb(ctools.cscript):
340 363 # Return FITS table
341 364 return fits.table(extname)
342 365  
343   - def _set_ogip_keywords(self, hdu, hdudoc, hduclas):
  366 + def _set_ogip_keywords(self, hdu, hdudoc, hduclas, irf):
344 367 """
345   - Set standard OGIP keywords for extension.
  368 + Set standard OGIP keywords for extension
346 369  
347   - Parameters:
348   - hdu: Header Data Unit.
349   - hdudoc: Documentation reference string
350   - hduclas: Array of HDUCLAS fields
  370 + Parameters
  371 + ----------
  372 + hdu : `~gammalig.GFitsHDU`
  373 + Header Data Unit
  374 + hdudoc : str
  375 + Documentation reference string
  376 + hduclas : list of str
  377 + List of HDUCLAS fields
  378 + irf : dict
  379 + IRF dictonary
351 380 """
352 381 # Set UTC date of file creation
353 382 utc = datetime.utcnow().isoformat()[:19]
354 383  
355 384 # Set keywords
356   - hdu.card("ORIGIN", "IRAP", "Name of organization making this file")
357   - hdu.card("DATE", utc, "File creation date (YYYY-MM-DDThh:mm:ss UTC)")
358   - hdu.card("TELESCOP", self._cal_tel, "Name of telescope")
359   - hdu.card("INSTRUME", self._cal_inst, "Name of instrument")
360   - hdu.card("DETNAM", self._cal_det, "Name of detector")
361   - hdu.card("HDUCLASS", "OGIP", "HDU class")
362   - hdu.card("HDUDOC", hdudoc, "HDU documentation")
  385 + hdu.card('ORIGIN', 'IRAP', 'Name of organization making this file')
  386 + hdu.card('DATE', utc, 'File creation date (YYYY-MM-DDThh:mm:ss UTC)')
  387 + hdu.card('TELESCOP', irf['CAL_TEL'], 'Name of telescope')
  388 + hdu.card('INSTRUME', irf['CAL_INST'], 'Name of instrument')
  389 + hdu.card('DETNAM', irf['CAL_DET'], 'Name of detector')
  390 + hdu.card('HDUCLASS', 'OGIP', 'HDU class')
  391 + hdu.card('HDUDOC', hdudoc, 'HDU documentation')
363 392 for i, item in enumerate(hduclas):
364   - key = "HDUCLAS%d" % (i+1)
365   - hdu.card(key, item, "HDU class")
366   - hdu.card("HDUVERS", "1.0.0", "HDU version")
  393 + key = 'HDUCLAS%d' % (i+1)
  394 + hdu.card(key, item, 'HDU class')
  395 + hdu.card('HDUVERS', '1.0.0', 'HDU version')
367 396  
368 397 # Return
369 398 return
370 399  
371   - def _create_cif(self):
  400 + def _add_cif_info(self, irf, dir):
372 401 """
373   - Create Calibration Index File (CIF) extension in CIF FITS file.
  402 + Add information to CIF extension
374 403  
375   - Parameters:
376   - None
377   - Keywords:
378   - None
  404 + Parameters
  405 + ----------
  406 + irf : dict
  407 + IRF metadata dictionary
  408 + dir : dict
  409 + Directory structure dictionary
379 410 """
380   - # Attach table to FITS file. Note that at this moment the FITS table
381   - # gets connected to the FITS file. Yet since nothing was yet written
382   - # to the FITS file, we cannot read anything from it.
383   - self._cif.append(calutils.create_cif_table())
384   -
385   - # Return
386   - return
  411 + # Set list of IRF component names
  412 + names = ['EA', 'PSF', 'EDISP', 'BGD']
  413 +
  414 + # Initialise CIF row index
  415 + row = self._hdu_cif.nrows()
  416 +
  417 + # Append rows for all components to CIF extension
  418 + self._hdu_cif.append_rows(len(names))
  419 +
  420 + # Add information for all components
  421 + for name in names:
  422 +
  423 + # Set dictionary keys for this component
  424 + key_dir = '%s_DIR' % name
  425 + key_file = '%s_FILE' % name
  426 + key_name = '%s_NAME' % name
  427 + key_desc = '%s_DESC' % name
  428 + key_bounds = '%s_BOUNDS' % name
  429 +
  430 + # Set generic information
  431 + self._hdu_cif['TELESCOP'][row] = irf['CAL_TEL']
  432 + self._hdu_cif['INSTRUME'][row] = irf['CAL_INST']
  433 + self._hdu_cif['DETNAM'][row] = irf['CAL_DET']
  434 + self._hdu_cif['FILTER'][row] = irf['CAL_FLT']
  435 + self._hdu_cif['CAL_DEV'][row] = 'ONLINE'
  436 + self._hdu_cif['CAL_CLAS'][row] = irf['CAL_CLASS']
  437 + self._hdu_cif['CAL_DTYP'][row] = irf['CAL_TYPE']
  438 + self._hdu_cif['CAL_VSD'][row] = irf['VAL_DATE']
  439 + self._hdu_cif['CAL_VST'][row] = irf['VAL_TIME']
  440 + self._hdu_cif['REF_TIME'][row] = irf['REF_TIME']
  441 + self._hdu_cif['CAL_QUAL'][row] = irf['CAL_QUAL']
  442 + self._hdu_cif['CAL_DATE'][row] = irf['CAL_DATE']
  443 +
  444 + # Set component specific information
  445 + self._hdu_cif['CAL_DIR'][row] = dir[key_dir]
  446 + self._hdu_cif['CAL_FILE'][row] = dir[key_file]
  447 + self._hdu_cif['CAL_CNAM'][row] = irf[key_name]
  448 + self._hdu_cif['CAL_DESC'][row] = irf[key_desc]
  449 + self._hdu_cif['CAL_XNO'][row] = 1
  450 + for i in range(9):
  451 + if i >= len(irf[key_bounds]):
  452 + self._hdu_cif['CAL_CBD'][row,i] = 'NONE'
  453 + else:
  454 + self._hdu_cif['CAL_CBD'][row,i] = irf[key_bounds][i]
387 455  
388   - def _add_cif_info(self):
389   - """
390   - Add information to CIF extension.
391   - """
392   - # Append 4 rows to CIF extension
393   - self._hdu_cif.append_rows(4)
394   -
395   - # Add generic information for these 4 rows
396   - for i in range(4):
397   -
398   - # Set row number
399   - row = i + self._hdu_cif.nrows() - 4
400   -
401   - # Set element
402   - self._hdu_cif["TELESCOP"][row] = self._cal_tel
403   - self._hdu_cif["INSTRUME"][row] = self._cal_inst
404   - self._hdu_cif["DETNAM"][row] = self._cal_det
405   - self._hdu_cif["FILTER"][row] = self._cal_flt
406   - self._hdu_cif["CAL_DEV"][row] = "ONLINE"
407   - self._hdu_cif["CAL_CLAS"][row] = self._cal_class
408   - self._hdu_cif["CAL_DTYP"][row] = self._cal_type
409   - self._hdu_cif["CAL_VSD"][row] = self._val_date
410   - self._hdu_cif["CAL_VST"][row] = self._val_time
411   - self._hdu_cif["REF_TIME"][row] = self._ref_time
412   - self._hdu_cif["CAL_QUAL"][row] = self._cal_qual
413   - self._hdu_cif["CAL_DATE"][row] = self._cal_date
414   -
415   - # Add effective area information
416   - row = self._hdu_cif.nrows() - 4
417   - self._hdu_cif["CAL_DIR"][row] = self._ea_dir
418   - self._hdu_cif["CAL_FILE"][row] = self._ea_file
419   - self._hdu_cif["CAL_CNAM"][row] = self._ea_name
420   - self._hdu_cif["CAL_DESC"][row] = self._ea_desc
421   - self._hdu_cif["CAL_XNO"][row] = 1
422   - for n in range(9):
423   - if n >= len(self._ea_bounds):
424   - self._hdu_cif["CAL_CBD"][row, n] = "NONE"
425   - else:
426   - self._hdu_cif["CAL_CBD"][row, n] = self._ea_bounds[n]
427   -
428   - # Add point spread function information
429   - row = self._hdu_cif.nrows() - 3
430   - self._hdu_cif["CAL_DIR"][row] = self._psf_dir
431   - self._hdu_cif["CAL_FILE"][row] = self._psf_file
432   - self._hdu_cif["CAL_CNAM"][row] = self._psf_name
433   - self._hdu_cif["CAL_DESC"][row] = self._psf_desc
434   - self._hdu_cif["CAL_XNO"][row] = 1
435   - for n in range(9):
436   - if n >= len(self._psf_bounds):
437   - self._hdu_cif["CAL_CBD"][row, n] = "NONE"
438   - else:
439   - self._hdu_cif["CAL_CBD"][row, n] = self._psf_bounds[n]
440   -
441   - # Add energy dispersion information
442   - row = self._hdu_cif.nrows() - 2
443   - self._hdu_cif["CAL_DIR"][row] = self._edisp_dir
444   - self._hdu_cif["CAL_FILE"][row] = self._edisp_file
445   - self._hdu_cif["CAL_CNAM"][row] = self._edisp_name
446   - self._hdu_cif["CAL_DESC"][row] = self._edisp_desc
447   - self._hdu_cif["CAL_XNO"][row] = 1
448   - for n in range(9):
449   - if n >= len(self._edisp_bounds):
450   - self._hdu_cif["CAL_CBD"][row, n] = "NONE"
451   - else:
452   - self._hdu_cif["CAL_CBD"][row, n] = self._edisp_bounds[n]
453   -
454   - # Add background information
455   - row = self._hdu_cif.nrows() - 1
456   - self._hdu_cif["CAL_DIR"][row] = self._bgd_dir
457   - self._hdu_cif["CAL_FILE"][row] = self._bgd_file
458   - self._hdu_cif["CAL_CNAM"][row] = self._bgd_name
459   - self._hdu_cif["CAL_DESC"][row] = self._bgd_desc
460   - self._hdu_cif["CAL_XNO"][row] = 1
461   - for n in range(9):
462   - if n >= len(self._bgd_bounds):
463   - self._hdu_cif["CAL_CBD"][row, n] = "NONE"
464   - else:
465   - self._hdu_cif["CAL_CBD"][row, n] = self._bgd_bounds[n]
  456 + # Increment row index
  457 + row += 1
466 458  
467 459 # Return
468 460 return
469 461  
470   - def _set_cif_keywords(self, hdu, name, bounds, desc):
  462 + def _set_cif_keywords(self, hdu, name, bounds, desc, irf):
471 463 """
472   - Set standard CIF keywords for extension.
  464 + Set standard CIF keywords for extension
  465 +
  466 + Parameters
  467 + ----------
  468 + hdu : `~gammalib.GFitsHDU`
  469 + FITS HDU
  470 + name : str
  471 + Calibration name
  472 + bounds : list of str
  473 + Calibration boundaries
  474 + desc : str
  475 + Calibration description
  476 + irf : dict
  477 + IRF dictionary
473 478 """
474 479 # Set standard CIF keywords
475   - hdu.card("CSYSNAME", "XMA_POL", "")
476   - hdu.card("CCLS0001", self._cal_class, "Calibration class")
477   - hdu.card("CDTP0001", self._cal_type, "Calibration type")
478   - hdu.card("CCNM0001", name, "Calibration name")
  480 + hdu.card('CSYSNAME', 'XMA_POL', '')
  481 + hdu.card('CCLS0001', irf['CAL_CLASS'], 'Calibration class')
  482 + hdu.card('CDTP0001', irf['CAL_TYPE'], 'Calibration type')
  483 + hdu.card('CCNM0001', name, 'Calibration name')
479 484  
480 485 # Set boundary keywords
481   - for n in range(9):
482   - keyname = "CBD%d0001" % (n+1)
483   - if n >= len(bounds):
484   - value = "NONE"
  486 + for i in range(9):
  487 + keyname = 'CBD%d0001' % (i+1)
  488 + if i >= len(bounds):
  489 + value = 'NONE'
485 490 else:
486   - value = bounds[n]
487   - hdu.card(keyname, value, "Calibration boundary")
  491 + value = bounds[i]
  492 + hdu.card(keyname, value, 'Calibration boundary')
488 493  
489 494 # Set validity keywords
490   - hdu.card("CVSD0001", self._val_date, "Calibration validity start date (UTC)")
491   - hdu.card("CVST0001", self._val_time, "Calibration validity start time (UTC)")
492   - hdu.card("CDEC0001", desc, "Calibration description")
  495 + hdu.card('CVSD0001', irf['VAL_DATE'],
  496 + 'Calibration validity start date (UTC)')
  497 + hdu.card('CVST0001', irf['VAL_TIME'],
  498 + 'Calibration validity start time (UTC)')
  499 + hdu.card('CDEC0001', desc, 'Calibration description')
493 500  
494 501 # Return
495 502 return
496 503  
497 504  
498 505 # ROOT specific private members
499   - def _root2caldb(self):
  506 + def _root2caldb(self, irf):
500 507 """
501   - Translate ROOT to CALDB information.
  508 + Translate ROOT to CALDB information
502 509  
503   - Parameters:
504   - filename - ROOT 2D performance filename.
  510 + Parameters
  511 + ----------
  512 + irf : dict
  513 + IRF dictionary
505 514 """
506 515 # Open ROOT performance file
507   - file = TFile(self["infile"].filename().url())
  516 + file = TFile(self['infile'].filename().url())
508 517  
509 518 # Create effective area
510   - self._root2ea(file, rebin=self["rebin"].boolean(),
511   - eascale=self["eascale"].real())
  519 + self._root2ea(file, irf)
512 520  
513 521 # Create point spread function
514   - self._root2psf(file, self["psftype"].string())
  522 + self._root2psf(file, irf)
515 523  
516 524 # Create energy dispersion
517   - self._root2edisp(file)
  525 + self._root2edisp(file, irf)
518 526  
519 527 # Create background
520   - #self._root2bgd(file)
521   - self._root2bgd3D(file)
  528 + self._root2bgd(file, irf)
522 529  
523 530 # Return
524 531 return
525 532  
526   - def _root2ea(self, file, rebin=False, eascale=1.0):
  533 + def _root2ea(self, file, irf):
527 534 """
528   - Translate ROOT to CALDB effective area extension. The following ROOT
529   - histograms are used:
530   -
531   - EffectiveAreaEtrue_offaxis -> EFFAREA
532   - EffectiveArea_offaxis -> EFFAREA_RECO
533   -
534   - Parameters:
535   - file: ROOT file.
  535 + Translate ROOT to CALDB effective area extension
  536 +
  537 + The following ROOT histograms are used:
  538 + - EffectiveAreaEtrue_offaxis -> EFFAREA
  539 + - EffectiveArea_offaxis -> EFFAREA_RECO
536 540  
537   - Keywords:
538   - rebin: Rebin Etrue histogram (useful for Prod1 IFAE runs)
539   - eascale: Effective area scaling factor
  541 + Parameters
  542 + ----------
  543 + file : `~ROOT.TFile`
  544 + ROOT file
  545 + irf : dict
  546 + IRF metadata dictionary
540 547 """
541   - # Continue only if effective area HDU has been opened
542   - if self._hdu_ea != None:
  548 + # Write header
  549 + if self._logTerse():
  550 + self._log('\n')
  551 + self._log.header1('Generate effective area extension')
  552 +
  553 + # Get relevant ROOT histograms
  554 + etrue = file.Get('EffectiveAreaEtrue_offaxis')
  555 + ereco = file.Get('EffectiveArea_offaxis')
  556 +
  557 + # If requested then normalize the 2D histogram on the on-axis 1D
  558 + # histogram. This assures that the 2D histogram has the same on-axis
  559 + # effective area dependence as the 1D histogram.
  560 + if self['norm1d'].boolean():
  561 + etrue_1D = file.Get('EffectiveAreaEtrue')
  562 + ereco_1D = file.Get('EffectiveArea')
  563 + self._renorm_onaxis(etrue, etrue_1D)
  564 + self._renorm_onaxis(ereco, ereco_1D)
  565 +
  566 + # If requested then rebin the effective area Etrue histogram. This
  567 + # increases the number of bins by a factor 10.
  568 + if self['rebin'].boolean():
  569 + etrue.RebinX(10)
  570 + neng = etrue.GetXaxis().GetNbins()
  571 + noffset = etrue.GetYaxis().GetNbins()
  572 + for ioff in range(noffset):
  573 + for ieng in range(neng):
  574 + value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
  575 + etrue.SetBinContent(ieng+1,ioff+1,value)
543 576  
544   - # Write header
545   - if self._logTerse():
546   - self._log('\n')
547   - self._log.header1('Generate effective area extension')
548   -
549   - # Get User parameters
550   - norm1d = self['norm1d'].boolean()
551   -
552   - # Allocate ROOT 2D arrays
553   - etrue = TH2F()
554   - ereco = TH2F()
555   - file.GetObject('EffectiveAreaEtrue_offaxis', etrue)
556   - file.GetObject('EffectiveArea_offaxis', ereco)
557   -
558   - # Normalize 2D histogram on on-axis 1D histogram
559   - if norm1d:
560   - etrue_1D = TH1F()
561   - ereco_1D = TH1F()
562   - file.GetObject('EffectiveAreaEtrue', etrue_1D)
563   - file.GetObject('EffectiveArea', ereco_1D)
564   - self._renorm_onaxis(etrue, etrue_1D)
565   - self._renorm_onaxis(ereco, ereco_1D)
566   -
567   - # Rebin etrue histogram
568   - if rebin:
569   - etrue.RebinX(10)
570   - neng = etrue.GetXaxis().GetNbins()
571   - noffset = etrue.GetYaxis().GetNbins()
572   - for ioff in range(noffset):
573   - for ieng in range(neng):
574   - value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
575   - etrue.SetBinContent(ieng+1,ioff+1,value)
  577 + # Get effective area multiplicator. This allows renormalising the
  578 + # effective area.
  579 + eascale = self['eascale'].real()
576 580  
577   - # Scale histograms (if needed)
578   - if eascale != 1.0:
579   - neng = etrue.GetXaxis().GetNbins()
580   - noffset = etrue.GetYaxis().GetNbins()
581   - for ioff in range(noffset):
582   - for ieng in range(neng):
583   - value = etrue.GetBinContent(ieng+1,ioff+1) * eascale
584   - etrue.SetBinContent(ieng+1,ioff+1,value)
585   - neng = ereco.GetXaxis().GetNbins()
586   - noffset = ereco.GetYaxis().GetNbins()
587   - for ioff in range(noffset):
588   - for ieng in range(neng):
589   - value = ereco.GetBinContent(ieng+1,ioff+1) * eascale
590   - ereco.SetBinContent(ieng+1,ioff+1,value)
  581 + # If renormalisation has been requested then do it now
  582 + if eascale != 1.0:
  583 + neng = etrue.GetXaxis().GetNbins()
  584 + noffset = etrue.GetYaxis().GetNbins()
  585 + for ioff in range(noffset):
  586 + for ieng in range(neng):
  587 + value = etrue.GetBinContent(ieng+1,ioff+1) * eascale
  588 + etrue.SetBinContent(ieng+1,ioff+1,value)
  589 + neng = ereco.GetXaxis().GetNbins()
  590 + noffset = ereco.GetYaxis().GetNbins()
  591 + for ioff in range(noffset):
  592 + for ieng in range(neng):
  593 + value = ereco.GetBinContent(ieng+1,ioff+1) * eascale
  594 + ereco.SetBinContent(ieng+1,ioff+1,value)
  595 +
  596 + # Set boundaries (use Ereco boundaries)
  597 + bounds = self._make_2D(ereco, self._hdu_ea, None, 'm2')
  598 + for b in bounds:
  599 + irf['EA_BOUNDS'].append(b)
591 600  
592   - # Write boundaries (use Ereco boundaries)
593   - bounds = self._make_2D(ereco, self._hdu_ea, None, "m2")
594   - for b in bounds:
595   - self._ea_bounds.append(b)
596   - self._set_cif_keywords(self._hdu_ea, self._ea_name, \
597   - self._ea_bounds, self._ea_desc)
  601 + # Write boundary keywords
  602 + self._set_cif_keywords(self._hdu_ea, irf['EA_NAME'],
  603 + irf['EA_BOUNDS'], irf['EA_DESC'], irf)
598 604  
599   - # EFFAREA
600   - self._make_2D(etrue, self._hdu_ea, "EFFAREA", "m2")
  605 + # Create "EFFAREA" data column
  606 + self._make_2D(etrue, self._hdu_ea, 'EFFAREA', 'm2')
601 607  
602   - # EFFAREA_RECO
603   - self._make_2D(ereco, self._hdu_ea, "EFFAREA_RECO", "m2")
  608 + # Create "EFFAREA_RECO" data column
  609 + self._make_2D(ereco, self._hdu_ea, 'EFFAREA_RECO', 'm2')
604 610  
605 611 # Return
606 612 return
607 613  
608   - def _root2psf(self, file, psftype):
  614 + def _root2psf(self, file, irf):
609 615 """
610   - Translate ROOT to CALDB point spread function extension.
611   - Parameters:
612   - file - ROOT file.
613   - psftype - PSF type (Gauss, King)
614   - Keywords:
615   - None
  616 + Translate ROOT to CALDB point spread function extension
  617 +
  618 + Parameters
  619 + ----------
  620 + file : `~ROOT.TFile`
  621 + ROOT file
  622 + irf : dict
  623 + IRF metadata dictionary
616 624 """
617   - # Continue only if point spread function HDU has been opened
618   - if self._hdu_psf != None:
  625 + # Write header
  626 + if self._logTerse():
  627 + self._log('\n')
  628 + self._log.header1('Generate point spread function extension')
619 629  
620   - # Write header
621   - if self._logTerse():
622   - self._log('\n')
623   - self._log.header1('Generate point spread function extension')
  630 + # King profile PSF
  631 + if self['psftype'].string() == 'King':
  632 + self._root2psf_king(file, irf)
624 633  
625   - # King profile PSF
626   - if psftype == 'King':
627   - self._root2psf_king(file)
628   - else:
629   - self._root2psf_gauss(file)
  634 + # ... otherwise use Gaussian profile PSF
  635 + else:
  636 + self._root2psf_gauss(file, irf)
630 637  
631 638 # Return
632 639 return
633 640  
634   - def _root2psf_gauss(self, file):
  641 + def _root2psf_gauss(self, file, irf):
635 642 """
636   - Translate ROOT to CALDB point spread function extension. The following
637   - ROOT histograms are used:
638   -
639   - 1/(2*pi*SIGMA_1) -> SCALE
640   - AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8)
641   - 0.0 -> AMPL_2
642   - 0.0 -> SIGMA_2
643   - 0.0 -> AMPL_3
644   - 0.0 -> SIGMA_3
645   -
646   - Parameters:
647   - file - ROOT file.
648   - Keywords:
649   - None
  643 + Translate ROOT to CALDB point spread function extension
  644 +
  645 + The following ROOT histograms are used:
  646 + - 1/(2*pi*SIGMA_1) -> SCALE
  647 + - AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8)
  648 + - 0.0 -> AMPL_2
  649 + - 0.0 -> SIGMA_2
  650 + - 0.0 -> AMPL_3
  651 + - 0.0 -> SIGMA_3
  652 +
  653 + Parameters
  654 + ----------
  655 + file : `~ROOT.TFile`
  656 + ROOT file
  657 + irf : dict
  658 + IRF metadata dictionary
650 659 """
651   - # Continue only if point spread function HDU has been opened
652   - if self._hdu_psf != None:
  660 + # Get relevant ROOT histograms
  661 + r68 = file.Get('AngRes_offaxis')
653 662  
654   - # Allocate ROOT 2D array
655   - r68 = TH2F()
  663 + # Extract number of bins in histogram
  664 + neng = r68.GetXaxis().GetNbins()
  665 + noffset = r68.GetYaxis().GetNbins()
656 666  
657   - # Read 68% containment histogram
658   - file.GetObject("AngRes_offaxis", r68)
659   - neng = r68.GetXaxis().GetNbins()
660   - noffset = r68.GetYaxis().GetNbins()
  667 + # Converts 68% -> 1 sigma
  668 + r68_to_sigma = 0.6624305
  669 + for ioff in range(noffset):
  670 + for ieng in range(neng):
  671 + sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma
  672 + r68.SetBinContent(ieng+1,ioff+1,sigma)
661 673  
662   - # Converts 68% -> 1 sigma
663   - r68_to_sigma = 0.6624305
664   - for ioff in range(noffset):
665   - for ieng in range(neng):
666   - sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma
667   - r68.SetBinContent(ieng+1,ioff+1,sigma)
  674 + # Compute scale histogram
  675 + scale = r68.Clone()
  676 + for ioff in range(noffset):
  677 + for ieng in range(neng):
  678 + integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1)
  679 + if integral > 0.0:
  680 + value = 1.0 / integral
  681 + else:
  682 + value = 0.0
  683 + scale.SetBinContent(ieng+1,ioff+1,value)
668 684  
669   - # Compute scale histogram
670   - scale = r68.Clone()
671   - for ioff in range(noffset):
672   - for ieng in range(neng):
673   - integral = 2.0 * math.pi * r68.GetBinContent(ieng+1,ioff+1)
674   - if integral > 0.0:
675   - value = 1.0 / integral
676   - else:
677   - value = 0.0
678   - scale.SetBinContent(ieng+1,ioff+1,value)
  685 + # Set zero histogram
  686 + zero = r68.Clone()
  687 + for ioff in range(noffset):
  688 + for ieng in range(neng):
  689 + zero.SetBinContent(ieng+1,ioff+1,0.0)
679 690  
680   - # Set zero histogram
681   - zero = r68.Clone()
682   - for ioff in range(noffset):
683   - for ieng in range(neng):
684   - zero.SetBinContent(ieng+1,ioff+1,0.0)
  691 + # Set boundaries
  692 + bounds = self._make_2D(r68, self._hdu_psf, None, 'deg')
  693 + for b in bounds:
  694 + irf['PSF_BOUNDS'].append(b)
  695 + irf['PSF_BOUNDS'].append('PSF(GAUSS)')
685 696  
686   - # Set boundaries
687   - bounds = self._make_2D(r68, self._hdu_psf, None, "deg")
688   - for b in bounds:
689   - self._psf_bounds.append(b)
690   - self._psf_bounds.append("PSF(GAUSS)")
691   - self._set_cif_keywords(self._hdu_psf, self._psf_name, \
692   - self._psf_bounds, self._psf_desc)
  697 + # Write boundary keywords
  698 + self._set_cif_keywords(self._hdu_psf, irf['PSF_NAME'],
  699 + irf['PSF_BOUNDS'], irf['PSF_DESC'], irf)
693 700  
694   - # SCALE
695   - self._make_2D(scale, self._hdu_psf, "SCALE", "")
  701 + # Create "SCALE" data column
  702 + self._make_2D(scale, self._hdu_psf, 'SCALE', '')
696 703  
697   - # SIGMA_1
698   - self._make_2D(r68, self._hdu_psf, "SIGMA_1", "deg")
  704 + # Create "SIGMA_1" data column
  705 + self._make_2D(r68, self._hdu_psf, 'SIGMA_1', 'deg')
699 706  
700   - # AMPL_2
701   - self._make_2D(zero, self._hdu_psf, "AMPL_2", "")
  707 + # Create "AMPL_2" data column
  708 + self._make_2D(zero, self._hdu_psf, 'AMPL_2', '')
702 709  
703   - # SIGMA_2
704   - self._make_2D(zero, self._hdu_psf, "SIGMA_2", "deg")
  710 + # Create "SIGMA_2" data column
  711 + self._make_2D(zero, self._hdu_psf, 'SIGMA_2', 'deg')
705 712  
706   - # AMPL_3
707   - self._make_2D(zero, self._hdu_psf, "AMPL_3", "")
  713 + # Create "AMPL_3" data column
  714 + self._make_2D(zero, self._hdu_psf, 'AMPL_3', '')
708 715  
709   - # SIGMA_3
710   - self._make_2D(zero, self._hdu_psf, "SIGMA_3", "deg")
  716 + # Create "SIGMA_3" data column
  717 + self._make_2D(zero, self._hdu_psf, 'SIGMA_3', 'deg')
711 718  
712 719 # Return
713 720 return
714 721  
715   - def _root2psf_king(self, file):
  722 + def _root2psf_king(self, file, irf):
716 723 """
717   - Translate ROOT to CALDB point spread function extension. The following
718   - ROOT histograms are used:
719   -
720   - AngRes_offaxis
721   - AngRes80_offaxis
  724 + Translate ROOT to CALDB point spread function extension
  725 +
  726 + The following ROOT histograms are used:
  727 + - AngRes_offaxis
  728 + - AngRes80_offaxis
722 729  
723   - Parameters:
724   - file - ROOT file.
725   - Keywords:
726   - None
  730 + Parameters
  731 + ----------
  732 + file : `~ROOT.TFile`
  733 + ROOT file
  734 + irf : dict
  735 + IRF metadata dictionary
727 736 """
728   - # Continue only if point spread function HDU has been opened
729   - if self._hdu_psf != None:
730   -
731   - # Allocate ROOT 2D arrays
732   - r68 = TH2F()
733   - r80 = TH2F()
  737 + # Get relevant ROOT histograms
  738 + r68 = file.Get('AngRes_offaxis')
  739 + r80 = file.Get('AngRes80_offaxis')
734 740  
735   - # Read 68% and 80% containment histograms
736   - file.GetObject("AngRes_offaxis", r68)
737   - file.GetObject("AngRes80_offaxis", r80)
738   - neng = r68.GetXaxis().GetNbins()
739   - noffset = r68.GetYaxis().GetNbins()
  741 + # Extract number of bins in histogram
  742 + neng = r68.GetXaxis().GetNbins()
  743 + noffset = r68.GetYaxis().GetNbins()
740 744  
741   - # Initialise parameter maps by cloning the r68 2D map
742   - gamma2D = r68.Clone()
743   - sigma2D = r68.Clone()
  745 + # Initialise parameter maps by cloning the r68 2D map
  746 + gamma2D = r68.Clone()
  747 + sigma2D = r68.Clone()
744 748  
745   - # Compute gamma and sigma values
746   - for ioff in range(noffset):
  749 + # Compute gamma and sigma values
  750 + for ioff in range(noffset):
747 751  
748   - # Initialise last results
749   - last_gamma = 0.0
750   - last_sigma = 0.0
  752 + # Initialise last results
  753 + last_gamma = 0.0
  754 + last_sigma = 0.0
751 755  
752   - for ieng in range(neng):
  756 + # Loop over all energies
  757 + for ieng in range(neng):
753 758  
754   - # Extract radii
755   - r_68 = r68.GetBinContent(ieng+1,ioff+1)
756   - r_80 = r80.GetBinContent(ieng+1,ioff+1)
757   -
758   - # Initialise results
759   - gamma = 0.0
760   - sigma = 0.0
761   -
762   - # Continue only if both radii are positive
763   - if r_68 > 0 and r_80 > 0:
764   -
765   - # Derive constants for equation to solve
766   - a = 1.0 - 0.68
767   - b = 1.0 - 0.80
768   - c = r_68*r_68/(r_80*r_80)
769   -
770   - # Solve equation (a^x-1)/(b^x-1)=c for x using secant
771   - # method. Stop when we are better than 1e-6.
772   - x1 = -0.5
773   - x2 = -1.0
774   - f1 = (math.pow(a,x1) - 1.0)/(math.pow(b,x1) - 1.0) - c
775   - f2 = (math.pow(a,x2) - 1.0)/(math.pow(b,x2) - 1.0) - c
776   - iter = 0
777   - while True:
778   - x = x1 - f1 * (x1-x2)/(f1-f2)
779   - f = (math.pow(a,x) - 1.0)/(math.pow(b,x) - 1.0) - c
780   - iter += 1
781   - if abs(f) < 1.0e-6:
782   - break
783   - else:
784   - f2 = f1
785   - x2 = x1
786   - f1 = f
787   - x1 = x
788   -
789   - # Compute gamma.
790   - if x < 0.0:
791   - gamma = 1.0 - 1.0/x
  759 + # Extract radii
  760 + r_68 = r68.GetBinContent(ieng+1,ioff+1)
  761 + r_80 = r80.GetBinContent(ieng+1,ioff+1)
  762 +
  763 + # Initialise results
  764 + gamma = 0.0
  765 + sigma = 0.0
  766 +
  767 + # Continue only if both radii are positive
  768 + if r_68 > 0 and r_80 > 0:
  769 +
  770 + # Derive constants for equation to solve
  771 + a = 1.0 - 0.68
  772 + b = 1.0 - 0.80
  773 + c = r_68*r_68/(r_80*r_80)
  774 +
  775 + # Solve equation (a^x-1)/(b^x-1)=c for x using secant
  776 + # method. Stop when we are better than 1e-6.
  777 + x1 = -0.5
  778 + x2 = -1.0
  779 + f1 = (math.pow(a,x1) - 1.0)/(math.pow(b,x1) - 1.0) - c
  780 + f2 = (math.pow(a,x2) - 1.0)/(math.pow(b,x2) - 1.0) - c
  781 + iter = 0
  782 + while True:
  783 + x = x1 - f1 * (x1-x2)/(f1-f2)
  784 + f = (math.pow(a,x) - 1.0)/(math.pow(b,x) - 1.0) - c
  785 + iter += 1
  786 + if abs(f) < 1.0e-6:
  787 + break
792 788 else:
793   - gamma = 1.0
  789 + f2 = f1
  790 + x2 = x1
  791 + f1 = f
  792 + x1 = x
  793 +
  794 + # Compute gamma.
  795 + if x < 0.0:
  796 + gamma = 1.0 - 1.0/x
  797 + else:
  798 + gamma = 1.0
794 799  
795   - # Compute sigma
796   - denominator = 2.0 * gamma * (math.pow(a, x) - 1.0)
  800 + # Compute sigma
  801 + denominator = 2.0 * gamma * (math.pow(a, x) - 1.0)
  802 + if denominator > 0.0:
  803 + sigma = r_68 * math.sqrt(1.0/denominator)
  804 + else:
  805 + denominator = 2.0 * gamma * (math.pow(b, x) - 1.0)
797 806 if denominator > 0.0:
798   - sigma = r_68 * math.sqrt(1.0/denominator)
  807 + sigma = r_80 * math.sqrt(1.0/denominator)
799 808 else:
800   - denominator = 2.0 * gamma * (math.pow(b, x) - 1.0)
801   - if denominator > 0.0:
802   - sigma = r_80 * math.sqrt(1.0/denominator)
803   - else:
804   - gamma = 0.0
805   - sigma = 0.0
806   -
807   - # Handle special case that no results were found.
808   - # This takes care of pixels that are ill defined
809   - # in the MC file.
810   - if gamma == 0.0 and sigma == 0.0:
811   - gamma = last_gamma
812   - sigma = last_sigma
813   -
814   - # Store surrent result as last result
815   - last_gamma = gamma
816   - last_sigma = sigma
817   -
818   - # Show results on console
819   - #print(ioff,ieng,gamma,sigma,x,f,iter)
820   -
821   - # Set bin contents
822   - gamma2D.SetBinContent(ieng+1,ioff+1,gamma)
823   - sigma2D.SetBinContent(ieng+1,ioff+1,sigma)
824   -
825   - # Set boundaries
826   - bounds = self._make_2D(r68, self._hdu_psf, None, "deg")
827   - for b in bounds:
828   - self._psf_bounds.append(b)
829   - self._psf_bounds.append("PSF(KING)")
830   - self._set_cif_keywords(self._hdu_psf, self._psf_name, \
831   - self._psf_bounds, self._psf_desc)
832   -
833   - # GAMMA
834   - self._make_2D(gamma2D, self._hdu_psf, "GAMMA", "")
835   -
836   - # SIGMA
837   - self._make_2D(sigma2D, self._hdu_psf, "SIGMA", "deg")
  809 + gamma = 0.0
  810 + sigma = 0.0
838 811  
839   - # Return
840   - return
  812 + # Handle special case that no results were found.
  813 + # This takes care of pixels that are ill defined
  814 + # in the MC file.
  815 + if gamma == 0.0 and sigma == 0.0:
  816 + gamma = last_gamma
  817 + sigma = last_sigma
841 818  
842   - def _root2edisp(self, file):
843   - """
844   - Translate ROOT to CALDB energy dispersion extension. The following ROOT
845   - histograms are used:
  819 + # Store surrent result as last result
  820 + last_gamma = gamma
  821 + last_sigma = sigma
846 822  
847   - EestOverEtrue_offaxis -> MATRIX
  823 + # Show results on console
  824 + #print(ioff,ieng,gamma,sigma,x,f,iter)
848 825  
849   - Parameters:
850   - file - ROOT file.
851   - Keywords:
852   - None
853   - """
854   - # Continue only if energy dispersion HDU has been opened
855   - if self._hdu_edisp != None:
  826 + # Set bin contents
  827 + gamma2D.SetBinContent(ieng+1,ioff+1,gamma)
  828 + sigma2D.SetBinContent(ieng+1,ioff+1,sigma)
856 829  
857   - # Write header
858   - if self._logTerse():
859   - self._log('\n')
860   - self._log.header1('Generate energy dispersion extension')
  830 + # Set boundaries
  831 + bounds = self._make_2D(r68, self._hdu_psf, None, 'deg')
  832 + for b in bounds:
  833 + irf['PSF_BOUNDS'].append(b)
  834 + irf['PSF_BOUNDS'].append('PSF(KING)')
861 835  
862   - # Allocate ROOT 3D array
863   - matrix = TH3F()
864   - file.GetObject('EestOverEtrue_offaxis', matrix)
  836 + # Write boundary keywords
  837 + self._set_cif_keywords(self._hdu_psf, irf['PSF_NAME'],
  838 + irf['PSF_BOUNDS'], irf['PSF_DESC'], irf)
865 839  
866   - # Set boundaries
867   - bounds = self._make_3D_migra(matrix, self._hdu_edisp, None, '')
868   - for b in bounds:
869   - self._edisp_bounds.append(b)
870   - self._set_cif_keywords(self._hdu_edisp, self._edisp_name, \
871   - self._edisp_bounds, self._edisp_desc)
  840 + # Create "GAMMA" data column
  841 + self._make_2D(gamma2D, self._hdu_psf, 'GAMMA', '')
872 842  
873   - # MATRIX
874   - self._make_3D_migra(matrix, self._hdu_edisp, 'MATRIX', '')
  843 + # Create "SIGMA" data column
  844 + self._make_2D(sigma2D, self._hdu_psf, 'SIGMA', 'deg')
875 845  
876 846 # Return
877 847 return
878 848  
879   - def _root2bgd(self, file):
  849 + def _root2edisp(self, file, irf):
880 850 """
881   - Translate ROOT to CALDB background extension. The following ROOT
882   - histograms are used:
883   -
884   - BGRatePerSqDeg_offaxis -> BGD
885   - BGRatePerSqDeg_offaxis -> BGD_RECO
  851 + Translate ROOT to CALDB energy dispersion extension
  852 +
  853 + The following ROOT histograms are used:
  854 + - EestOverEtrue_offaxis -> MATRIX
886 855  
887   - Parameters:
888   - file - ROOT file.
889   - Keywords:
890   - None
  856 + Parameters
  857 + ----------
  858 + file : `~ROOT.TFile`
  859 + ROOT file
  860 + irf : dict
  861 + IRF metadata dictionary
891 862 """
892   - # Continue only if background HDU has been opened
893   - if self._hdu_bgd != None:
894   -
895   - # Write header
896   - if self._logTerse():
897   - self._log('\n')
898   - self._log.header1('Generate 2D background extension')
  863 + # Write header
  864 + if self._logTerse():
  865 + self._log('\n')
  866 + self._log.header1('Generate energy dispersion extension')
899 867  
900   - # Allocate ROOT 2D array
901   - array = TH2F()
902   - file.GetObject('BGRatePerSqDeg_offaxis', array)
  868 + # Get relevant ROOT histograms
  869 + matrix = file.Get('EestOverEtrue_offaxis')
903 870  
904   - # Set boundaries
905   - bounds = self._make_2D(array, self._hdu_bgd, None, 'deg')
906   - for b in bounds:
907   - self._bgd_bounds.append(b)
908   - self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \
909   - self._bgd_bounds, self._bgd_desc)
  871 + # Set boundaries
  872 + bounds = self._make_3D_migra(matrix, self._hdu_edisp, None, '')
  873 + for b in bounds:
  874 + irf['EDISP_BOUNDS'].append(b)
910 875  
911   - # BGD
912   - self._make_2D(array, self._hdu_bgd, 'BGD', '')
  876 + # Write boundary keywords
  877 + self._set_cif_keywords(self._hdu_edisp, irf['EDISP_NAME'],
  878 + irf['EDISP_BOUNDS'], irf['EDISP_DESC'], irf)
913 879  
914   - # BGD_RECO
915   - self._make_2D(array, self._hdu_bgd, 'BGD_RECO', '')
  880 + # Create "MATRIX" data column
  881 + self._make_3D_migra(matrix, self._hdu_edisp, 'MATRIX', '')
916 882  
917 883 # Return
918 884 return
919 885  
920   - def _root2bgd3D(self, file):
  886 + def _root2bgd(self, file, irf):
921 887 """
922 888 Translate ROOT to CALDB background extension.
923 889  
924 890 The following ROOT histograms are used:
925   - BGRatePerSqDeg_offaxis -> BGD
  891 + - BGRatePerSqDeg_offaxis -> BGD
926 892  
927 893 Parameters
928 894 ----------
929 895 file : `~ROOT.TFile`
930   - ROOT file.
  896 + ROOT file
  897 + irf : dict
  898 + IRF metadata dictionary
931 899 """
932   - # Continue only if background HDU has been opened
933   - if self._hdu_bgd != None:
934   -
935   - # Write header
936   - if self._logTerse():
937   - self._log('\n')
938   - self._log.header1('Generate 3D background extension')
  900 + # Write header
  901 + if self._logTerse():
  902 + self._log('\n')
  903 + self._log.header1('Generate 3D background extension')
939 904  
940   - # Get User parameters
941   - norm1d = self['norm1d'].boolean()
942   - bgdethres = self['bgdethres'].real()
943   - bgdinfill = self['bgdinfill'].boolean()
  905 + # Get relevant ROOT histograms
  906 + array = file.Get('BGRatePerSqDeg_offaxis')
944 907  
945   - # Allocate ROOT 2D array
946   - array = TH2F()
947   - file.GetObject('BGRatePerSqDeg_offaxis', array)
  908 + # Replace 2D histogram values by power law extrapolation
  909 + self._plaw_replace(array, self['bgdethres'].real())
948 910  
949   - # Replace 2D histogram values by power law extrapolation
950   - self._plaw_replace(array, bgdethres)
  911 + # If requested then fill-in empty values in 2D histogram
  912 + if self['bgdinfill'].boolean():
  913 + self._infill_bgd(array)
951 914  
952   - # Optionally fill in empty values in 2D histogram
953   - if bgdinfill:
954   - self._infill_bgd(array)
  915 + # If requested then normalize the 2D histogram on the on-axis 1D
  916 + # histogram. This assures that the 2D histogram has the same on-axis
  917 + # effective area dependence as the 1D histogram.
  918 + if self['norm1d'].boolean():
  919 + array_1D = file.Get('BGRatePerSqDeg')
  920 + self._renorm_onaxis(array, array_1D)
955 921  
956   - # Optionally normalize 2D histogram on on-axis 1D histogram
957   - if norm1d:
958   - array_1D = TH1F()
959   - file.GetObject('BGRatePerSqDeg', array_1D)
960   - self._renorm_onaxis(array, array_1D)
  922 + # Set boundaries
  923 + bounds = self._make_3D(array, self._hdu_bgd, None, 'deg')
  924 + for b in bounds:
  925 + irf['BGD_BOUNDS'].append(b)
961 926  
962   - # Set boundaries
963   - bounds = self._make_3D(array, self._hdu_bgd, None, 'deg')
964   - for b in bounds:
965   - self._bgd_bounds.append(b)
966   - self._set_cif_keywords(self._hdu_bgd, self._bgd_name,
967   - self._bgd_bounds, self._bgd_desc)
  927 + # Write boundary keywords
  928 + self._set_cif_keywords(self._hdu_bgd, irf['BGD_NAME'],
  929 + irf['BGD_BOUNDS'], irf['BGD_DESC'], irf)
968 930  
969   - # BGD (reconstructed energy)
970   - self._make_3D(array, self._hdu_bgd, 'BGD', '1/s/MeV/sr')
  931 + # Create "BGD" data column
  932 + self._make_3D(array, self._hdu_bgd, 'BGD', '1/s/MeV/sr')
971 933  
972 934 # Return
973 935 return
974 936  
975 937 def _append_column(self, hdu, name, unit, axis, log=False):
976 938 """
977   - Append column to HDU.
  939 + Append column to HDU
978 940  
979   - Parameters:
980   - hdu: HDU
981   - name: Column name
982   - unit: Column unit
983   - axis: ROOT histogram axis
984   -
985   - Keywords:
986   - log: Logarithmic axis
  941 + Parameters
  942 + ----------
  943 + hdu : `~gammalib.GFitsHDU`
  944 + FITS HDU
  945 + name : str
  946 + Column name
  947 + unit : str
  948 + Column unit
  949 + axis : `~ROOT.TAxis`
  950 + ROOT histogram axis
  951 + log : bool, optional
  952 + Axis is logarithmic
987 953 """
988 954 # Continue only if columns does not yet exist
989 955 if not hdu.contains(name):
... ... @@ -1024,9 +990,9 @@ class csroot2caldb(ctools.cscript):
1024 990 Parameters
1025 991 ----------
1026 992 hist2D : `~ROOT.TH2F`
1027   - ROOT 2D histogram.
  993 + ROOT 2D histogram
1028 994 hist1D : `~ROOT.TH1F`
1029   - ROOT 1D histogram.
  995 + ROOT 1D histogram
1030 996 """
1031 997 # Get 2D dimensions
1032 998 neng = hist2D.GetXaxis().GetNbins()
... ... @@ -1101,9 +1067,9 @@ class csroot2caldb(ctools.cscript):
1101 1067 Parameters
1102 1068 ----------
1103 1069 hist2D : `~ROOT.TH2F`
1104   - ROOT 2D histogram.
  1070 + ROOT 2D histogram
1105 1071 ethres : float
1106   - Energy threshold.
  1072 + Energy threshold
1107 1073 """
1108 1074 # Extract energy and offset angle vectors
1109 1075 energies = hist2D.GetXaxis()
... ... @@ -1260,7 +1226,7 @@ class csroot2caldb(ctools.cscript):
1260 1226 Parameters
1261 1227 ----------
1262 1228 hist2D : `~ROOT.TH2F`
1263   - ROOT 2D histogram.
  1229 + ROOT 2D histogram
1264 1230 """
1265 1231 # Extract energy and offset angle vectors
1266 1232 energies = hist2D.GetXaxis()
... ... @@ -1316,19 +1282,24 @@ class csroot2caldb(ctools.cscript):
1316 1282  
1317 1283 def _make_2D(self, array, hdu, name, unit, scale=1.0):
1318 1284 """
1319   - Make 2D matrix as function of energy and offset angle from a
1320   - ROOT 2D histogram. If the HDU has already the energy and
1321   - offset angle columns, this method will simply add another data
1322   - column.
1323   - If name==None, the method will not append any data column.
1324   -
1325   - Parameters:
1326   - array - ROOT 2D histogram.
1327   - hdu - FITS HDU.
1328   - name - Data column name.
1329   - unit - Data unit.
1330   - Keywords:
1331   - scale - Scaling factor for histogram values.
  1285 + Make 2D data column as function of energy and offset angle
  1286 +
  1287 + If the HDU has already the energy and offset angle columns, this method
  1288 + will simply add another data column. If name==None, the method will not
  1289 + append any data column.
  1290 +
  1291 + Parameters
  1292 + ----------
  1293 + array : `~ROOT.TH2F`
  1294 + ROOT 2D histogram
  1295 + hdu : `~gammalib.GFitsHDU`
  1296 + FITS HDU
  1297 + name : str
  1298 + Data column name
  1299 + unit : str
  1300 + Data unit
  1301 + scale : float, optional
  1302 + Scaling factor for histogram values
1332 1303 """
1333 1304 # Extract energy and offset angle vectors
1334 1305 energies = array.GetXaxis()
... ... @@ -1337,10 +1308,10 @@ class csroot2caldb(ctools.cscript):
1337 1308 noffset = offsets.GetNbins()
1338 1309  
1339 1310 # Append axis columns to HDU
1340   - self._append_column(hdu, "ENERG_LO", "TeV", energies, log=True)
1341   - self._append_column(hdu, "ENERG_HI", "TeV", energies, log=True)
1342   - self._append_column(hdu, "THETA_LO", "deg", offsets)
1343   - self._append_column(hdu, "THETA_HI", "deg", offsets)
  1311 + self._append_column(hdu, 'ENERG_LO', 'TeV', energies, log=True)
  1312 + self._append_column(hdu, 'ENERG_HI', 'TeV', energies, log=True)
  1313 + self._append_column(hdu, 'THETA_LO', 'deg', offsets)
  1314 + self._append_column(hdu, 'THETA_HI', 'deg', offsets)
1344 1315  
1345 1316 # Append array column
1346 1317 if name != None and not hdu.contains(name):
... ... @@ -1354,13 +1325,13 @@ class csroot2caldb(ctools.cscript):
1354 1325 hdu[name][0,index] = value * scale
1355 1326  
1356 1327 # Collect boundary information
1357   - bd_eng = "ENERG(%.4f-%.2f)TeV" % \
  1328 + bd_eng = 'ENERG(%.4f-%.2f)TeV' % \
1358 1329 (pow(10.0, energies.GetBinLowEdge(1)), \
1359 1330 pow(10.0, energies.GetBinUpEdge(neng)))
1360   - bd_off = "THETA(%.2f-%.2f)deg" % \
  1331 + bd_off = 'THETA(%.2f-%.2f)deg' % \
1361 1332 (offsets.GetBinLowEdge(1), \
1362 1333 offsets.GetBinUpEdge(noffset))
1363   - bd_phi = "PHI(0-360)deg"
  1334 + bd_phi = 'PHI(0-360)deg'
1364 1335 bounds = [bd_eng, bd_off, bd_phi]
1365 1336  
1366 1337 # Return boundary information
... ... @@ -1368,23 +1339,22 @@ class csroot2caldb(ctools.cscript):
1368 1339  
1369 1340 def _make_3D(self, array, hdu, name, unit):
1370 1341 """
1371   - Generate background 3D cube
  1342 + Make 3D data column as function of DETX, DETY and energy
1372 1343  
1373   - Generates 3D cube as function of DETX, DETY and energy from a
1374   - ROOT 2D histogram. If the HDU has already the energy and
1375   - offset angle columns, this method will simply add another data
1376   - column. If name==None, the method will not append any data column.
  1344 + If the HDU has already the energy and offset angle columns, this method
  1345 + will simply add another data column. If name==None, the method will not
  1346 + append any data column.
1377 1347  
1378 1348 Parameters
1379 1349 ----------
1380 1350 array : `~ROOT.TH2F`
1381   - ROOT 2D histogram.
  1351 + ROOT 2D histogram
1382 1352 hdu : `~gammalib.GFitsHDU`
1383   - FITS HDU.
  1353 + FITS HDU
1384 1354 name : str
1385   - Data column name.
  1355 + Data column name
1386 1356 unit : str
1387   - Data unit.
  1357 + Data unit
1388 1358 """
1389 1359 # Get User parameters
1390 1360 scale = self['bgdscale'].real()
... ... @@ -1468,21 +1438,6 @@ class csroot2caldb(ctools.cscript):
1468 1438 value = binsq / deg2sr / (ewidth[ieng])
1469 1439 hdu[name][0,index] = value * scale
1470 1440  
1471   - # Debugging
1472   - if False:
1473   - for ieng in range(neng):
1474   - sum = 0.0
1475   - for ioff in range(noffset):
1476   - index = ieng + ioff * neng
1477   - binsq = array.GetBinContent(ieng+1,ioff+1)
1478   - value = binsq / deg2sr / ewidth[ieng]
1479   - sum += value
1480   - if ioff == 0:
1481   - peak = value
1482   - peaksq = binsq
1483   - print(pow(10.0, energies.GetBinLowEdge(ieng+1))+
1484   - pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq)
1485   -
1486 1441 # Collect boundary information
1487 1442 bd_detx = 'DETX(%.2f-%.2f)deg' % (-det_max, det_max)
1488 1443 bd_dety = 'DETY(%.2f-%.2f)deg' % (-det_max, det_max)
... ... @@ -1496,18 +1451,24 @@ class csroot2caldb(ctools.cscript):
1496 1451  
1497 1452 def _make_3D_migra(self, array, hdu, name, unit, scale=1.0):
1498 1453 """
1499   - Make 3D cube as function of ETRUE, MIGRA and THETA form ROOT 3D
1500   - histogram. If the HDU has already the columns, this method will
1501   - simply add another data column.
1502   - If name==None, the method will not append any data column.
1503   -
1504   - Parameters:
1505   - array - ROOT 3D histogram.
1506   - hdu - FITS HDU.
1507   - name - Data column name.
1508   - unit - Data unit.
1509   - Keywords:
1510   - scale - Scaling factor for histogram values.
  1454 + Make 3D data column as function of ETRUE, MIGRA and THETA
  1455 +
  1456 + If the HDU has already the energy and offset angle columns, this method
  1457 + will simply add another data column. If name==None, the method will not
  1458 + append any data column.
  1459 +
  1460 + Parameters
  1461 + ----------
  1462 + array : `~ROOT.TH3F`
  1463 + ROOT 2D histogram
  1464 + hdu : `~gammalib.GFitsHDU`
  1465 + FITS HDU
  1466 + name : str
  1467 + Data column name
  1468 + unit : str
  1469 + Data unit
  1470 + scale : float, optional
  1471 + Scaling factor for histogram values
1511 1472 """
1512 1473 # Extract Etrue, Eobs/Etrue and offset angle vectors
1513 1474 etrue = array.GetXaxis()
... ... @@ -1522,12 +1483,12 @@ class csroot2caldb(ctools.cscript):
1522 1483 pow(10.0, etrue.GetBinLowEdge(ieng+1)+6.0))
1523 1484  
1524 1485 # Append axis columns to HDU
1525   - self._append_column(hdu, "ETRUE_LO", "TeV", etrue, log=True)
1526   - self._append_column(hdu, "ETRUE_HI", "TeV", etrue, log=True)
1527   - self._append_column(hdu, "MIGRA_LO", "", migra)
1528   - self._append_column(hdu, "MIGRA_HI", "", migra)
1529   - self._append_column(hdu, "THETA_LO", "deg", offsets)
1530   - self._append_column(hdu, "THETA_HI", "deg", offsets)
  1486 + self._append_column(hdu, 'ETRUE_LO', 'TeV', etrue, log=True)
  1487 + self._append_column(hdu, 'ETRUE_HI', 'TeV', etrue, log=True)
  1488 + self._append_column(hdu, 'MIGRA_LO', '', migra)
  1489 + self._append_column(hdu, 'MIGRA_HI', '', migra)
  1490 + self._append_column(hdu, 'THETA_LO', 'deg', offsets)
  1491 + self._append_column(hdu, 'THETA_HI', 'deg', offsets)
1531 1492  
1532 1493 # Append migration matrix to HDU
1533 1494 if name != None and not hdu.contains(name):
... ... @@ -1542,16 +1503,16 @@ class csroot2caldb(ctools.cscript):
1542 1503 hdu[name][0,index] = value * scale
1543 1504  
1544 1505 # Collect boundary information
1545   - bd_eng = "ETRUE(%.4f-%.2f)TeV" % \
  1506 + bd_eng = 'ETRUE(%.4f-%.2f)TeV' % \
1546 1507 (pow(10.0, etrue.GetBinLowEdge(1)),
1547 1508 pow(10.0, etrue.GetBinUpEdge(netrue)))
1548   - bd_migra = "MIGRA(%.3f-%.3f)" % \
  1509 + bd_migra = 'MIGRA(%.3f-%.3f)' % \
1549 1510 (migra.GetBinLowEdge(1),
1550 1511 migra.GetBinUpEdge(nmigra))
1551   - bd_off = "THETA(%.2f-%.2f)deg" % \
  1512 + bd_off = 'THETA(%.2f-%.2f)deg' % \
1552 1513 (offsets.GetBinLowEdge(1),
1553 1514 offsets.GetBinUpEdge(noffset))
1554   - bd_phi = "PHI(0-360)deg"
  1515 + bd_phi = 'PHI(0-360)deg'
1555 1516 bounds = [bd_eng, bd_migra, bd_off, bd_phi]
1556 1517  
1557 1518 # Return boundary information
... ... @@ -1561,7 +1522,7 @@ class csroot2caldb(ctools.cscript):
1561 1522 # Public methods
1562 1523 def run(self):
1563 1524 """
1564   - Run the script.
  1525 + Run the script
1565 1526 """
1566 1527 # Warn if ROOT module is missing
1567 1528 if not _has_root:
... ... @@ -1575,25 +1536,28 @@ class csroot2caldb(ctools.cscript):
1575 1536 # Get parameters
1576 1537 self._get_parameters()
1577 1538  
1578   - # Initialise IRF
1579   - self._initialise()
  1539 + # Initialise IRF metadata
  1540 + irf = self._init_metadata()
  1541 +
  1542 + # Create directory structure
  1543 + dir = self._make_dirs('file', irf)
1580 1544  
1581   - # Open calibration file
1582   - self._open("file")
  1545 + # Open calibration files
  1546 + self._open(irf, dir)
1583 1547  
1584 1548 # Translate ROOT to CALDB information
1585 1549 if _has_root:
1586   - self._root2caldb()
  1550 + self._root2caldb(irf)
1587 1551  
1588   - # Close calibration file
1589   - self._close()
  1552 + # Close calibration files
  1553 + self._close(irf, dir)
1590 1554  
1591 1555 # Return
1592 1556 return
1593 1557  
1594 1558 def execute(self):
1595 1559 """
1596   - Execute the script.
  1560 + Execute the script
1597 1561 """
1598 1562 # Open logfile
1599 1563 self.logFileOpen()
... ...