Alexander Ziegler / ctools

Commit 523323c92ca0e3e723116d51564677176a3a2ff4

Authored by Jürgen Knödlseder
2 parents a446d51c add41cca

Merge branch '1828-quality-of-iact-cscripts' into integration 

# Conflicts:
#	cscripts/csiactcopy.py
#	cscripts/csiactdata.py
Inline Side-by-side
Showing 9 changed files with 407 additions and 113 deletions
cscripts/Makefile.am
  View file @523323c
... ... @@ -152,7 +152,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \
152 152 $(top_srcdir)/test/data/models_spectral.xml \
153 153 $(top_srcdir)/test/data/model_filefunction.txt \
154 154 $(top_srcdir)/test/data/model_map.fits.gz \
155   - $(top_srcdir)/test/data/model_mapcube.fits.gz
  155 + $(top_srcdir)/test/data/model_mapcube.fits.gz \
  156 + $(top_srcdir)/test/data/iact_runlist.dat
156 157  
157 158 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata
158 159 testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
... ...
cscripts/csiactcopy.py
  View file @523323c
... ... @@ -91,7 +91,18 @@ class csiactcopy(ctools.cscript):
91 91 def _copy(self, source, clobber):
92 92 """
93 93 Copy file to outpath
  94 +
  95 + Parameters
  96 + ----------
  97 + source : string
  98 + Path of file to be copied
  99 +
  100 + Returns
  101 + -------
  102 + filesize : float
  103 + Size of the file that was copied
94 104 """
  105 +
95 106 # Get file destination
96 107 destination = os.path.join(self._outpath,
97 108 os.path.relpath(source, self._remote_base))
... ... @@ -147,10 +158,21 @@ class csiactcopy(ctools.cscript):
147 158 def _merge(self, localfits, remotefits, hduname, clobber):
148 159 """
149 160 Merge remote and local fits files
150   -
151 161 If the local fits file is not present, a new one is created.
  162 +
  163 + Parameters
  164 + ----------
  165 + localfits : string
  166 + Path of local index FITS file
  167 + remotefits : string
  168 + Path of remote index FITS file
  169 + hduname : string
  170 + Name of HDU extension to be merged
  171 + clobber : bool
  172 + Flag if remote content should overwrite local content
152 173 """
153   - # ...
  174 +
  175 + # Logging
154 176 if self._logTerse():
155 177 self._log_value('Local file', localfits)
156 178 self._log_value('Remote file', remotefits)
... ... @@ -214,7 +236,7 @@ class csiactcopy(ctools.cscript):
214 236 tmp_hdu = remote[hduname]
215 237  
216 238 # Get old size of remote HDU
217   - size = remote_hdu.size()
  239 + size = remote_hdu.nrows()
218 240  
219 241 # Extend remote HDU
220 242 remote_hdu.append_rows(tmp_hdu.nrows())
... ... @@ -229,8 +251,8 @@ class csiactcopy(ctools.cscript):
229 251 robs_col = remote_hdu['OBS_ID']
230 252 for obs_id in robs_col:
231 253 remote_obs.append(obs_id)
232   -
233   - # Dump the number of remote and local entries
  254 +
  255 + # Log entries of index files
234 256 if self._logTerse():
235 257 self._log_value('Remote entries', len(remote_obs))
236 258 self._log_value('Local entries', len(local_obs))
... ... @@ -250,12 +272,18 @@ class csiactcopy(ctools.cscript):
250 272  
251 273 # Remove local obs_id entries
252 274 table_has_obsid = True
  275 +
  276 + # Loop over obs ids and overwrite with remote content
253 277 while table_has_obsid:
  278 +
  279 + # First remove old rows
254 280 for i in range(local_hdu.nrows()):
255 281 if remote_obs_id == local_hdu['OBS_ID'][i]:
256 282 local_hdu.remove_rows(i,1)
257 283 removed_rows += 1
258 284 break
  285 +
  286 + # Replace with remote content
259 287 table_has_obsid = False
260 288 for i in range(local_hdu.nrows()):
261 289 if remote_obs_id == local_hdu['OBS_ID'][i]:
... ... @@ -290,8 +318,8 @@ class csiactcopy(ctools.cscript):
290 318  
291 319 # Add remote HDUs
292 320 local_hdu.insert_rows(old_local_rows, remote_hdu.nrows())
293   -
294   - # Dump ...
  321 +
  322 + # Log actions
295 323 if self._logTerse():
296 324 if clobber:
297 325 what = 'Removed local rows'
... ... @@ -344,6 +372,8 @@ class csiactcopy(ctools.cscript):
344 372 # Check if runlist file is available
345 373 if self._runlist.exists():
346 374 runfile = open(self._runlist.url())
  375 +
  376 + # Read runlist
347 377 for line in runfile.readlines():
348 378 if len(line) == 0:
349 379 continue
... ... @@ -351,6 +381,8 @@ class csiactcopy(ctools.cscript):
351 381 continue
352 382 if len(line.split()) > 0:
353 383 self._runs.append(int(line.split()[0]))
  384 +
  385 + # Close runlist file
354 386 runfile.close()
355 387  
356 388 # Logging
... ... @@ -428,8 +460,8 @@ class csiactcopy(ctools.cscript):
428 460 # Initialise flag if SIZE column is present
429 461 has_size = table.contains('SIZE')
430 462  
  463 + # Initialise file size
431 464 if has_size:
432   - # Initialise file size
433 465 cp_size = 0
434 466  
435 467 # Initialise remote observation IDs
... ... @@ -439,6 +471,8 @@ class csiactcopy(ctools.cscript):
439 471  
440 472 # Log runs that are not available remotely
441 473 for run in self._runs:
  474 +
  475 + # Check for run not in remote data store
442 476 if not run in remote_ids:
443 477 if self._logNormal():
444 478 self._log('Skip observation '+str(run)+': ')
... ... @@ -529,16 +563,22 @@ class csiactcopy(ctools.cscript):
529 563 self._log('\n')
530 564 self._log.header1('Copying files')
531 565  
532   - # intialise counter
  566 + # Intialise counter
533 567 k = 0
534 568  
535   - # initialise values for logging
  569 + # Initialise values for logging
536 570 last_fraction = 0.0
537 571 fraction_increment = 20.0
  572 +
  573 + # Use 10% step increase
538 574 if self._logNormal():
539 575 fraction_increment = 10.0
  576 +
  577 + # Use 5% step increase
540 578 if self._logTerse():
541 579 fraction_increment = 5.0
  580 +
  581 + # Use 2% step increase
542 582 if self._logExplicit():
543 583 fraction_increment = 2.0
544 584  
... ... @@ -552,6 +592,8 @@ class csiactcopy(ctools.cscript):
552 592 # Log progress
553 593 fraction = float(k) / float(len(files)) * 100.0
554 594 while fraction > last_fraction:
  595 +
  596 + # Print status of copying procedure
555 597 if self._logNormal() and not self._logVerbose():
556 598 self._log_value('Status', '%d %%' % int(last_fraction))
557 599 last_fraction += fraction_increment
... ... @@ -565,7 +607,7 @@ class csiactcopy(ctools.cscript):
565 607 total_size += filesize
566 608 n_copied += 1
567 609  
568   - # Increment counter
  610 + # Increment file counter
569 611 k += 1
570 612  
571 613 # Logging
... ... @@ -573,14 +615,17 @@ class csiactcopy(ctools.cscript):
573 615 self._log.parformat('Status')
574 616 self._log('Finished')
575 617 self._log('\n')
  618 +
  619 + # Logging about index files
576 620 if self._logTerse():
577 621 self._log('\n')
578 622 self._log.header1('Updating index files')
579 623  
580   - # Build local index file names
  624 + # Build local hdu index file name
581 625 local_hdu = os.path.join(self._outpath,
582 626 os.path.relpath(remote_hdu,
583 627 self._remote_base))
  628 + # Build local obs index file name
584 629 local_obs = os.path.join(self._outpath,
585 630 os.path.relpath(remote_obs,
586 631 self._remote_base))
... ... @@ -596,6 +641,7 @@ class csiactcopy(ctools.cscript):
596 641 # Merge remote index files with local files
597 642 self._merge(local_hdu, remote_hdu, 'HDU_INDEX', self._clobber())
598 643  
  644 + # Logging
599 645 if self._logTerse():
600 646 self._log('\n')
601 647 self._log.header3('OBS index')
... ... @@ -607,7 +653,8 @@ class csiactcopy(ctools.cscript):
607 653 # If all files were copied, just copy index files too
608 654 self._copy(remote_hdu, self._clobber())
609 655 self._copy(remote_obs, self._clobber())
610   -
  656 +
  657 + # Logging
611 658 if self._logTerse():
612 659 self._log('\n')
613 660 self._log.header3('Master index file')
... ...
cscripts/csiactdata.py
  View file @523323c
... ... @@ -95,6 +95,12 @@ class csiactdata(ctools.cscript):
95 95 def _log_configuration(self, config):
96 96 """
97 97 Log configuration
  98 +
  99 + Parameters
  100 + ----------
  101 + config : dict
  102 + Dictionary containing data configuration
  103 +
98 104 """
99 105 # Log only if the level is at least "Terse"
100 106 if self._logTerse():
... ... @@ -107,6 +113,8 @@ class csiactdata(ctools.cscript):
107 113  
108 114 # Print additional information
109 115 for key in config:
  116 +
  117 + # Skip keys that were treated above
110 118 if key in ['name','hduindx','obsindx']:
111 119 continue
112 120 self._log_value(key, config[key])
... ...
cscripts/csiactobs.py
  View file @523323c
... ... @@ -47,9 +47,10 @@ class csiactobs(ctools.cscript):
47 47 """
48 48 # Set name and version
49 49 self._name = 'csiactobs'
50   - self._version = '1.1.0'
  50 + self._version = '1.2.0'
51 51  
52 52 # Initialise some members
  53 + self._obs = gammalib.GObservations()
53 54 self._ebounds = gammalib.GEbounds()
54 55 self._datapath = os.getenv('VHEFITS','')
55 56 self._inmodels = gammalib.GModels()
... ... @@ -59,8 +60,8 @@ class csiactobs(ctools.cscript):
59 60 self._runlist = []
60 61 self._runlistfile = gammalib.GFilename()
61 62 self._bkgpars = 0
62   - self._outmodel = gammalib.GFilename()
63   - self._outobs = gammalib.GFilename()
  63 + #self._outmodel = gammalib.GFilename()
  64 + #self._outobs = gammalib.GFilename()
64 65 self._master_indx = ''
65 66 self._use_bkg_scale = False
66 67 self._ev_hiera = ['']
... ... @@ -80,8 +81,7 @@ class csiactobs(ctools.cscript):
80 81 self._subdir = ''
81 82 self._debug = False
82 83  
83   - # Initialise application by calling the appropriate class
84   - # constructor.
  84 + # Initialise application by calling the appropriate class constructor
85 85 self._init_cscript(argv)
86 86  
87 87 # Return
... ... @@ -91,7 +91,7 @@ class csiactobs(ctools.cscript):
91 91 # Private methods
92 92 def _get_parameters(self):
93 93 """
94   - Get parameters from parfile and setup the observation.
  94 + Get parameters from parfile and setup the observation
95 95 """
96 96 # Set data path
97 97 if self._datapath == '':
... ... @@ -105,7 +105,7 @@ class csiactobs(ctools.cscript):
105 105 self._datapath = gammalib.expand_env(self._datapath)
106 106  
107 107 # Read FITS production
108   - self._prodname = self['prodname'].string()
  108 + self._prodname = self['prodname'].string()
109 109  
110 110 # Read runlist file if list not already filled
111 111 if len(self._runlist) == 0:
... ... @@ -125,13 +125,12 @@ class csiactobs(ctools.cscript):
125 125 runfile.close()
126 126  
127 127 # Read number of background parameters
128   - self._bkgpars = self['bkgpars'].integer()
129   -
130   - # Output model file
131   - self._outmodel = self['outmodel'].filename()
132   -
133   - # Observation outfile
134   - self._outobs = self['outobs'].filename()
  128 + self._bkgpars = self['bkgpars'].integer()
  129 +
  130 + # Query ahead output parameters
  131 + if self._read_ahead():
  132 + self['outmodel'].filename()
  133 + self['outobs'].filename()
135 134  
136 135 # Master index file name
137 136 self._master_indx = self['master_indx'].string()
... ... @@ -174,27 +173,45 @@ class csiactobs(ctools.cscript):
174 173  
175 174 # Get HDUs
176 175 for config in configs:
  176 +
  177 + # Check if prodname is present
177 178 if self._prodname == config['name']:
178 179 self._hdu_index = str(os.path.join(self._datapath, config['hduindx']))
179 180 self._obs_index = str(os.path.join(self._datapath, config['obsindx']))
  181 +
  182 + # Leave loop if index file names were found
180 183 break
181 184  
182   - # Check HDUs
  185 + # Check index files
183 186 if self._hdu_index == '' or self._obs_index == '':
184 187 raise RuntimeError('*** ERROR: FITS data store "'+self._prodname+'" not available. Run csiactdata to get a list of available storage names')
  188 +
  189 + # Check HDU names
185 190 filename = gammalib.GFilename(self._hdu_index+'[HDU_INDEX]')
186 191 if not filename.is_fits():
187 192 raise RuntimeError('*** ERROR: HDU index file "'+self._hdu_index+'[HDU_INDEX]" for FITS data store "'+self._prodname+'" not available. Check your master index file or run csiactdata to get a list of available storage names.')
188 193  
189 194 # Check for existence of 'BKG_SCALE' in the observation index file if required
190 195 if self._use_bkg_scale:
  196 +
  197 + # Create filename
191 198 filename = gammalib.GFilename(self._obs_index+'[OBS_INDEX]')
  199 +
  200 + # Check if it is a FITS file
192 201 if filename.is_fits():
  202 +
  203 + # Open FITS file
193 204 fits = gammalib.GFits(self._obs_index)
  205 +
  206 + # Check if column "BKG_SCALE" is found and signal its possible usage
194 207 if not fits['OBS_INDEX'].contains('BKG_SCALE'):
195 208 self._use_bkg_scale = False
  209 +
  210 + # Close FITS file
196 211 fits.close()
  212 +
197 213 else:
  214 + # Signal that there is no background scale
198 215 self._use_bkg_scale = False
199 216  
200 217 # Create base data directory from hdu index file location
... ... @@ -209,33 +226,37 @@ class csiactobs(ctools.cscript):
209 226 # Return
210 227 return
211 228  
212   - def runlist(self, runlist):
213   - """
214   - Set observation list
  229 + def _background_spectrum(self, prefactor, index, emin = 0.01, emax = 100.0):
215 230 """
  231 + Create a background spectrum model dependent on user parameters
216 232  
217   - # Check if convertable to a python list
218   - try:
219   -
220   - # Store a copy of input argument
221   - self._runlist = list(runlist)
222   -
223   - except:
224   -
225   - # raise error if wrong type detected
226   - raise RuntimeError('Argument of csiactobs.runlist() must be of type "list"')
  233 + Parameters
  234 + ----------
  235 + prefactor : float
  236 + Prefactor parameter of spectral model
  237 + index : float
  238 + Index parameter of spectral model
  239 + emin : float
  240 + Minimum energy (in case a spectral node function is required)
  241 + emax : float
  242 + Maximum energy (in case a spectral node function is required)
227 243  
228   - # Return
229   - return
230   -
231   - def _background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0):
  244 + Returns
  245 + -------
  246 + spec : `~gammalib.GModelSpectral()`
  247 + Spectral model for the background shape
  248 + """
232 249  
233 250 # Handle constant spectral model
234 251 if index == 0.0 and self._bkgpars <= 1:
235 252 spec = gammalib.GModelSpectralConst()
  253 +
  254 + # Set parameter range
236 255 spec['Normalization'].min(prefactor / self._bkg_range_factor)
237 256 spec['Normalization'].max(prefactor * self._bkg_range_factor)
238 257 spec['Normalization'].value(prefactor)
  258 +
  259 + # Freeze or release normalisation parameter
239 260 if self._bkgpars == 0:
240 261 spec['Normalization'].fix()
241 262 else:
... ... @@ -245,8 +266,10 @@ class csiactobs(ctools.cscript):
245 266  
246 267 # Create power law model
247 268 if self._bkgpars <= 2:
248   - e = gammalib.GEnergy(1.0,'TeV')
249   - spec = gammalib.GModelSpectralPlaw(prefactor, index, e)
  269 +
  270 + # Set Power Law model with arbitrary pivot energy
  271 + pivot = gammalib.GEnergy(1.0,'TeV')
  272 + spec = gammalib.GModelSpectralPlaw(prefactor, index, pivot)
250 273  
251 274 # Set parameter ranges
252 275 spec[0].min(prefactor / self._bkg_range_factor)
... ... @@ -269,19 +292,31 @@ class csiactobs(ctools.cscript):
269 292 else:
270 293  
271 294 # Create reference powerlaw
272   - plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,'TeV'))
  295 + pivot = gammalib.GEnergy(1.0,'TeV')
  296 + plaw = gammalib.GModelSpectralPlaw(prefactor, index, pivot)
273 297  
274 298 # Create spectral model and energy values
275 299 spec = gammalib.GModelSpectralNodes()
  300 +
  301 + # Create logarithmic energy nodes
276 302 bounds = gammalib.GEbounds(self._bkgpars,gammalib.GEnergy(emin,'TeV'),gammalib.GEnergy(emax,'TeV'), True)
  303 +
  304 + # Loop over bounds and set intensity value
277 305 for i in range(bounds.size()):
278 306 energy = bounds.elogmean(i)
279   - value = plaw.eval(energy, gammalib.GTime())
  307 + value = plaw.eval(energy)
  308 +
  309 + # Append energy, value - tuple to node function
280 310 spec.append(energy, value)
  311 +
  312 + # Loop over parameters
281 313 for par in spec:
282 314  
  315 + # Fix energy nodes
283 316 if 'Energy' in par.name():
284 317 par.fix()
  318 +
  319 + # Release intensity nodes
285 320 elif 'Intensity' in par.name():
286 321 value = par.value()
287 322 par.scale(value)
... ... @@ -292,13 +327,36 @@ class csiactobs(ctools.cscript):
292 327 return spec
293 328  
294 329 def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100):
  330 + """
  331 + Create an IACT background model
  332 +
  333 + Parameters
  334 + ----------
  335 + telescope : string
  336 + Name of telescope
  337 + obs_id : string
  338 + Observation ID
  339 + bkg_scale : float
  340 + Background scaling factor
  341 + bkgtype : string
  342 + Type of background (irf,aeff, or gauss)
295 343  
  344 + Returns
  345 + -------
  346 + model : `~gammalib.GModelData()`
  347 + Background model for IACT observation
  348 + """
  349 +
296 350 # handle IrfBackground
297 351 if bkgtype == 'irf':
  352 +
  353 + # Set parameters to have a constant spectral model
298 354 prefactor = 1.0
299 355 index = 0.0
300 356 prefactor *= bkg_scale
301   - spec = self._background_spectrum(obs_id, prefactor, index, emin, emax)
  357 +
  358 + # Create background spectrum
  359 + spec = self._background_spectrum(prefactor, index, emin, emax)
302 360  
303 361 # Create background model instance
304 362 bck = gammalib.GCTAModelIrfBackground(spec)
... ... @@ -306,7 +364,7 @@ class csiactobs(ctools.cscript):
306 364 # Set AeffBackground
307 365 elif bkgtype == 'aeff':
308 366 prefactor = bkg_scale * self._bkg_aeff_norm
309   - spec = self._background_spectrum(obs_id, prefactor, self._bkg_aeff_index, emin, emax)
  367 + spec = self._background_spectrum(prefactor, self._bkg_aeff_index, emin, emax)
310 368  
311 369 # Create background model instance
312 370 bck = gammalib.GCTAModelAeffBackground(spec)
... ... @@ -314,12 +372,14 @@ class csiactobs(ctools.cscript):
314 372 # Set Gaussian Background
315 373 elif bkgtype == 'gauss':
316 374 prefactor = bkg_scale * self._bkg_gauss_norm
317   - spec = self._background_spectrum(obs_id, prefactor, self._bkg_gauss_index, emin, emax)
  375 +
  376 + # Create background spectrum
  377 + spec = self._background_spectrum(prefactor, self._bkg_gauss_index, emin, emax)
318 378 radial = gammalib.GCTAModelRadialGauss(self._bkg_gauss_sigma)
319 379 bck = gammalib.GCTAModelRadialAcceptance(radial, spec)
320 380  
321 381 else:
322   - msg = 'Background type "'+self._bkgtype+'" unsupported'
  382 + msg = 'Background type "'+bkgtype+'" unsupported'
323 383 raise RuntimeError(msg)
324 384  
325 385 # Copy model
... ... @@ -383,6 +443,7 @@ class csiactobs(ctools.cscript):
383 443 eventfile = aefffile = psffile = edispfile = bkgfile = ''
384 444 eventhdu = aeffhdu = ''
385 445  
  446 + # Initialise types and formats
386 447 types = []
387 448 formats = []
388 449 for i in range(hduindx_hdu.nrows()):
... ... @@ -402,7 +463,11 @@ class csiactobs(ctools.cscript):
402 463 hduindx_hdu['FILE_NAME'][index])
403 464 eventhdu = hduindx_hdu['HDU_NAME'][index]
404 465 eventfile += '['+eventhdu+']'
  466 +
  467 + # Check if psf file is available
405 468 if not gammalib.GFilename(eventfile).is_fits():
  469 +
  470 + # Print warning that observation will be skipped since no events are present
406 471 if self._logTerse():
407 472 self._log('Skipping observation '+str(obs_id)+
408 473 ': eventfile "'+eventfile+'" not found\n')
... ... @@ -421,7 +486,11 @@ class csiactobs(ctools.cscript):
421 486 hduindx_hdu['FILE_NAME'][index])
422 487 aeffhdu = hduindx_hdu['HDU_NAME'][index]
423 488 aefffile += '['+aeffhdu+']'
  489 +
  490 + # Check if psf file is available
424 491 if not gammalib.GFilename(aefffile).is_fits():
  492 +
  493 + # Print warning that observation will be skipped since no effective area is present
425 494 if self._logTerse():
426 495 self._log('Skipping observation '+str(obs_id)+
427 496 ': effective area "'+aefffile+'" not found\n')
... ... @@ -439,7 +508,11 @@ class csiactobs(ctools.cscript):
439 508 hduindx_hdu['FILE_DIR'][index],
440 509 hduindx_hdu['FILE_NAME'][index])
441 510 psffile += '['+hduindx_hdu['HDU_NAME'][index]+']'
  511 +
  512 + # Check if psf file is available
442 513 if not gammalib.GFilename(psffile).is_fits():
  514 +
  515 + # Print warning that observation will be skipped since no PSF is present
443 516 if self._logTerse():
444 517 self._log('Skipping observation '+str(obs_id)+
445 518 ': point spread function "'+psffile+
... ... @@ -458,7 +531,11 @@ class csiactobs(ctools.cscript):
458 531 hduindx_hdu['FILE_DIR'][index],
459 532 hduindx_hdu['FILE_NAME'][index])
460 533 edispfile += '['+hduindx_hdu['HDU_NAME'][index]+']'
  534 +
  535 + # Check if edisp file is available
461 536 if not gammalib.GFilename(edispfile).is_fits():
  537 +
  538 + # Print warning that edisp cannot be used for this observation
462 539 if self._logTerse():
463 540 self._log('Warning: observation '+str(obs_id)+
464 541 ' has no energy dispersion "'+edispfile+
... ... @@ -480,17 +557,22 @@ class csiactobs(ctools.cscript):
480 557 bkgfile += '['+hduindx_hdu['HDU_NAME'][index]+']'
481 558 if 'irf' in bkg_mod_hierarchy and not gammalib.GFilename(bkgfile).is_fits():
482 559 bkg_mod_hierarchy.remove('irf')
  560 +
  561 + # Throw warning if IRF background cannot be used
483 562 if self._logTerse():
484 563 self._log('Warning: observation '+str(obs_id)+
485 564 ' has no background information (file="'+
486 565 bkgfile+'"). IRF background cannot be used\n')
487 566 bkgfile = ''
488 567 if len(bkg_mod_hierarchy) == 0:
  568 +
  569 + # Skip observation if no background can be used
489 570 if self._logTerse():
490 571 self._log('Skipping observation '+str(obs_id)+
491 572 ': No background can be used\n')
492 573 continue
493 574 else:
  575 + # Log if we fall back to next background approach
494 576 if self._logTerse():
495 577 self._log('Observation '+str(obs_id)+
496 578 ': Falling back to background "'+
... ... @@ -499,15 +581,21 @@ class csiactobs(ctools.cscript):
499 581 # Close hdu index file
500 582 hduindx.close()
501 583  
502   - # Handle background scale information if available
  584 + # Initialise background scale
503 585 bkg_scale = 1.0
  586 +
  587 + # Check if background scale should be used
504 588 if self._use_bkg_scale:
  589 +
  590 + # Read background scale from fits file
505 591 obsindx = gammalib.GFits(self._obs_index+'[OBS_INDEX]'+obs_selection)
506 592 bkg_scale = obsindx['OBS_INDEX']['BKG_SCALE'][0]
507 593 obsindx.close()
508 594  
509   - # Open fits file to determine the observation name
  595 + # Open event FITS file
510 596 fits = gammalib.GFits(eventfile)
  597 +
  598 + # Get object and telescope strings from event hdu
511 599 events = fits[eventhdu]
512 600 object_name = events.string('OBJECT')
513 601 telescope = events.string('TELESCOP')
... ... @@ -550,6 +638,7 @@ class csiactobs(ctools.cscript):
550 638 self._log('Adding observation '+str(obs_id)+
551 639 ' ("'+object_name+'")\n')
552 640  
  641 + # log more details
553 642 if self._logExplicit():
554 643 self._log(' Event file: '+eventfile+'\n')
555 644 self._log(' Effective area: '+aefffile+'\n')
... ... @@ -615,10 +704,13 @@ class csiactobs(ctools.cscript):
615 704 self._log('\n')
616 705  
617 706 else:
  707 + # Print warning if no observation could be be used from runlist
618 708 self._log.header2('WARNING: No observation from given runlist available')
619 709  
620 710 # Append models provided by 'inmodels' if necessary
621   - if not self._inmodels == None:
  711 + if self._inmodels.size() > 0:
  712 +
  713 + # Logging
622 714 if self._logTerse():
623 715 self._log('\n')
624 716 self._log.header1('Appending models')
... ... @@ -650,17 +742,64 @@ class csiactobs(ctools.cscript):
650 742  
651 743 def save(self):
652 744 """
653   - Save model definition XML file
  745 + Save observation definition and model definition XML file
654 746 """
655   - # Save observation XML file
656   - self._xml.save(self._outobs)
  747 + # Write header
  748 + if self._logTerse():
  749 + self._log('\n')
  750 + self._log.header1('Save XML files')
  751 +
  752 + # Get filenames
  753 + outobs = self['outobs'].filename()
  754 + outmodel = self['outmodel'].filename()
  755 +
  756 + # Log file name
  757 + if self._logTerse():
  758 + self._log(gammalib.parformat('Obs. definition XML file'))
  759 + self._log(outobs.url())
  760 + self._log('\n')
  761 +
  762 + # Save observation definition XML file
  763 + self._xml.save(outobs.url())
  764 +
  765 + # Log file name
  766 + if self._logTerse():
  767 + self._log(gammalib.parformat('Model definiton XML file'))
  768 + self._log(outmodel.url())
  769 + self._log('\n')
657 770  
658 771 # Save model XML file
659   - self._models.save(self._outmodel)
  772 + self._models.save(outmodel.url())
660 773  
661 774 # Return
662 775 return
663 776  
  777 + def runlist(self, runlist):
  778 + """
  779 + Set runlist
  780 +
  781 + Parameters
  782 + ----------
  783 + runlist : list
  784 + List of observation IDs
  785 +
  786 + Raises
  787 + ------
  788 + RuntimeError
  789 + Input runlist is not a Python list
  790 + """
  791 + # If the runlist is convertable to a Python list then store the
  792 + # runlist in the class
  793 + try:
  794 + self._runlist = list(runlist)
  795 +
  796 + # ... otherwise raise an exception
  797 + except:
  798 + raise RuntimeError('Argument is not a Python list')
  799 +
  800 + # Return
  801 + return
  802 +
664 803 def ebounds(self):
665 804 """
666 805 Return runlist energy boundaries
... ... @@ -670,20 +809,20 @@ class csiactobs(ctools.cscript):
670 809  
671 810 def obs(self):
672 811 """
673   - Returns GObservations object
  812 + Return observations container
674 813 """
675   - # Initialise observations
676   - obs = gammalib.GObservations()
  814 + # Clear observations
  815 + self._obs.clear()
677 816  
678 817 # Read observations if XML is filled
679 818 if self._xml.size():
680   - obs.read(self._xml)
  819 + self._obs.read(self._xml)
681 820  
682 821 # Assign models
683   - obs.models(self._models)
  822 + self._obs.models(self._models)
684 823  
685 824 # Return observations
686   - return obs
  825 + return self._obs
687 826  
688 827  
689 828 # ======================== #
... ...
pyext/ctool.i
  View file @523323c
... ... @@ -198,7 +198,7 @@ cscript._init_cscript = _init_cscript
198 198  
199 199 # Log the value of a parameter
200 200 def _log_value(self, name, value):
201   - self._log.parformat(name)
  201 + self._log.parformat(str(name))
202 202 self._log(str(value))
203 203 self._log('\n')
204 204 cscript._log_value = _log_value
... ...
test/data/iact_runlist.dat 0 → 100644
  View file @523323c
  1 +0
... ...
test/data/iactdata/master.json
  View file @523323c
... ... @@ -3,7 +3,14 @@
3 3 {
4 4 "name": "unit-test",
5 5 "obsindx": "obs-index.fits.gz",
6   - "hduindx": "hdu-index.fits.gz"
  6 + "hduindx": "hdu-index.fits.gz",
  7 + "comment": "this is a comment"
  8 + },
  9 + {
  10 + "name": "unit-test-that-doesnt-exist",
  11 + "obsindx": "badfile.fits.gz",
  12 + "hduindx": "unavailable_file.fits.gz",
  13 + "comment": "this is a comment"
7 14 }
8 15 ]
9 16 }
... ...
test/test_csiactdata.py
  View file @523323c
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This scripts performs unit tests for the csiactdata script.
  3 +# This scripts performs unit tests for the csiactdata script
4 4 #
5 5 # Copyright (C) 2016 Juergen Knoedlseder
6 6 #
... ... @@ -108,6 +108,13 @@ class Test(test):
108 108 # Run csiactdata script and save run list
109 109 iactdata.logFileOpen() # Make sure we get a log file
110 110 iactdata.run()
  111 +
  112 + # Get available configs
  113 + names = iactdata.names()
  114 +
  115 + # Check if there is at least one config available
  116 + self.test_assert(len(names) > 0,
  117 + 'Check successful execution from Python')
111 118  
112 119 # Return
113 120 return
... ...
test/test_csiactobs.py
  View file @523323c
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This scripts performs unit tests for the csiactobs script.
  3 +# This scripts performs unit tests for the csiactobs script
4 4 #
5 5 # Copyright (C) 2016 Juergen Knoedlseder
6 6 #
... ... @@ -44,6 +44,7 @@ class Test(test):
44 44  
45 45 # Set data members
46 46 self._datapath = self._datadir + '/iactdata'
  47 + self._runlist = self._datadir + '/iact_runlist.dat'
47 48  
48 49 # Return
49 50 return
... ... @@ -74,10 +75,10 @@ class Test(test):
74 75 # Setup csiactobs command
75 76 cmd = csiactobs+' datapath="'+self._datapath+'"'+ \
76 77 ' prodname="unit-test"'+ \
77   - ' infile="runlist_cmd1.dat"'+ \
  78 + ' infile="'+self._runlist+'"'+ \
78 79 ' bkgpars=1'+\
79   - ' outobs="obs_iactobs_cmd1.xml"'+ \
80   - ' outmodel="bgd_iactobs_cmd1.xml"'+ \
  80 + ' outobs="csiactobs_obs_cmd1.xml"'+ \
  81 + ' outmodel="csiactobs_bgd_cmd1.xml"'+ \
81 82 ' logfile="csiactobs_cmd1.log" chatter=1'
82 83  
83 84 # Check if execution of wrong command fails
... ... @@ -89,23 +90,36 @@ class Test(test):
89 90 'Check successful execution from command line')
90 91  
91 92 # Check observation definition XML file
92   - self._check_obsdef('obs_iactobs_cmd1.xml')
  93 + self._check_obsdef('csiactobs_obs_cmd1.xml')
93 94  
94 95 # Check model definition XML file
95   - self._check_moddef('bgd_iactobs_cmd1.xml')
  96 + self._check_moddef('csiactobs_bgd_cmd1.xml')
96 97  
97 98 # Setup csiactobs command
98 99 cmd = csiactobs+' datapath="data_path_that_does_not_exist"'+ \
99 100 ' prodname="unit-test"'+ \
100   - ' infile="runlist_cmd1.dat"'+ \
  101 + ' infile="'+self._runlist+'"'+ \
101 102 ' bkgpars=1'+\
102   - ' outobs="obs_iactobs_cmd2.xml"'+ \
103   - ' outmodel="bgd_iactobs_cmd2.xml"'+ \
  103 + ' outobs="csiactobs_obs_cmd2.xml"'+ \
  104 + ' outmodel="csiactobs_bgd_cmd2.xml"'+ \
104 105 ' logfile="csiactobs_cmd2.log" chatter=1'
105 106  
106 107 # Check if execution failed
107 108 self.test_assert(self._execute(cmd) != 0,
108   - 'Check invalid input file when executed from command line')
  109 + 'Check invalid input datapath when executed from command line')
  110 +
  111 + # Setup csiactobs command
  112 + cmd = csiactobs+' datapath="'+self._datapath+'"'+ \
  113 + ' prodname="unit-test-doesnt-exist"'+ \
  114 + ' infile="'+self._runlist+'"'+ \
  115 + ' bkgpars=1'+\
  116 + ' outobs="csiactobs_obs_cmd3.xml"'+ \
  117 + ' outmodel="csiactobs_bgd_cmd3.xml"'+ \
  118 + ' logfile="csiactobs_cmd3.log" chatter=1'
  119 +
  120 + # Check if execution failed
  121 + self.test_assert(self._execute(cmd) != 0,
  122 + 'Check invalid input prodname when executed from command line')
109 123  
110 124 # Return
111 125 return
... ... @@ -119,10 +133,10 @@ class Test(test):
119 133 iactobs = cscripts.csiactobs()
120 134 iactobs['datapath'] = self._datapath
121 135 iactobs['prodname'] = 'unit-test'
122   - iactobs['infile'] = 'runlist_py1.dat'
  136 + iactobs['infile'] = self._runlist
123 137 iactobs['bkgpars'] = 1
124   - iactobs['outobs'] = 'obs_iactobs_py1.xml'
125   - iactobs['outmodel'] = 'bgd_iactobs_py1.xml'
  138 + iactobs['outobs'] = 'csiactobs_obs_py1.xml'
  139 + iactobs['outmodel'] = 'csiactobs_bgd_py1.xml'
126 140 iactobs['logfile'] = 'csiactobs_py1.log'
127 141 iactobs['chatter'] = 2
128 142  
... ... @@ -132,10 +146,88 @@ class Test(test):
132 146 iactobs.save()
133 147  
134 148 # Check observation definition XML file
135   - self._check_obsdef('obs_iactobs_py1.xml')
  149 + self._check_obsdef('csiactobs_obs_py1.xml')
136 150  
137 151 # Check model definition XML file
138   - self._check_moddef('bgd_iactobs_py1.xml')
  152 + self._check_moddef('csiactobs_bgd_py1.xml')
  153 +
  154 + # Create test runlist
  155 + runlist = ['0']
  156 +
  157 + # Set-up csiactobs using a runlist with 2 background parameters
  158 + iactobs = cscripts.csiactobs()
  159 + iactobs['datapath'] = self._datapath
  160 + iactobs['prodname'] = 'unit-test'
  161 + iactobs['bkgpars'] = 2
  162 + iactobs['outobs'] = 'csiactobs_obs_py2.xml'
  163 + iactobs['outmodel'] = 'csiactobs_bgd_py2.xml'
  164 + iactobs['logfile'] = 'csiactobs_py2.log'
  165 + iactobs['chatter'] = 3
  166 + iactobs.runlist(runlist)
  167 +
  168 + # Run csiactobs script and save run list
  169 + iactobs.logFileOpen() # Make sure we get a log file
  170 + iactobs.run()
  171 + iactobs.save()
  172 +
  173 + # Test return functions
  174 + self.test_value(iactobs.obs().size(), 1,
  175 + 'Check number of observations in container')
  176 + self.test_value(iactobs.ebounds().size(), 0,
  177 + 'Check number of energy boundaries')
  178 +
  179 + # Check observation definition XML file
  180 + self._check_obsdef('csiactobs_obs_py2.xml')
  181 +
  182 + # Check model definition XML file
  183 + self._check_moddef('csiactobs_bgd_py2.xml')
  184 +
  185 + # Set-up csiactobs with a large number of free parameters and "aeff"
  186 + # background
  187 + iactobs = cscripts.csiactobs()
  188 + iactobs['datapath'] = self._datapath
  189 + iactobs['prodname'] = 'unit-test'
  190 + iactobs['infile'] = self._runlist
  191 + iactobs['bkgpars'] = 8
  192 + iactobs['bkg_mod_hiera'] = 'aeff'
  193 + iactobs['outobs'] = 'csiactobs_obs_py3.xml'
  194 + iactobs['outmodel'] = 'csiactobs_bgd_py3.xml'
  195 + iactobs['logfile'] = 'csiactobs_py3.log'
  196 + iactobs['chatter'] = 4
  197 +
  198 + # Execute csiactobs script
  199 + iactobs.execute()
  200 +
  201 + # Check observation definition XML file
  202 + self._check_obsdef('csiactobs_obs_py3.xml')
  203 +
  204 + # Check model definition XML file
  205 + self._check_moddef('csiactobs_bgd_py3.xml')
  206 +
  207 + # Set-up csiactobs with a "gauss" background and "inmodel" parameter
  208 + iactobs = cscripts.csiactobs()
  209 + iactobs['datapath'] = self._datapath
  210 + iactobs['inmodel'] = self._model
  211 + iactobs['prodname'] = 'unit-test'
  212 + iactobs['infile'] = self._runlist
  213 + iactobs['bkgpars'] = 1
  214 + iactobs['bkg_mod_hiera'] = 'gauss'
  215 + iactobs['outobs'] = 'NONE'
  216 + iactobs['outmodel'] = 'NONE'
  217 + iactobs['logfile'] = 'csiactobs_py4.log'
  218 + iactobs['chatter'] = 4
  219 +
  220 + # Run csiactobs script
  221 + iactobs.logFileOpen() # Make sure we get a log file
  222 + iactobs.run()
  223 +
  224 + # Check number of observations
  225 + self.test_value(iactobs.obs().size(), 1,
  226 + 'Check number of observations in container')
  227 +
  228 + # Check number of models
  229 + self.test_value(iactobs.obs().models().size(), 3,
  230 + 'Check number of models in container')
139 231  
140 232 # Return
141 233 return
... ... @@ -154,34 +246,26 @@ class Test(test):
154 246 # Check number of observations
155 247 self.test_value(obs.size(), 1,
156 248 'Check for single observation in XML file')
157   - self.test_assert(obs[0].eventfile().file() == 'events.fits.gz',
158   - 'Check event file name')
159   - self.test_assert(obs[0].eventfile().extname() == 'EVENTS',
160   - 'Check event extension name')
161   - self.test_assert(rsp.aeff().filename().file() ==
162   - 'irf_file.fits.gz',
163   - 'Check effective area name')
164   - self.test_assert(rsp.aeff().filename().extname() ==
165   - 'EFFECTIVE AREA',
166   - 'Check effective area extension name')
167   - self.test_assert(rsp.psf().filename().file() ==
168   - 'irf_file.fits.gz',
169   - 'Check point spread function name')
170   - self.test_assert(rsp.psf().filename().extname() ==
171   - 'POINT SPREAD FUNCTION',
172   - 'Check point spread function extension name')
173   - self.test_assert(rsp.edisp().filename().file() ==
174   - 'irf_file.fits.gz',
175   - 'Check energy dispersion name')
176   - self.test_assert(rsp.edisp().filename().extname() ==
177   - 'ENERGY DISPERSION',
178   - 'Check energy dispersion extension name')
179   - self.test_assert(rsp.background().filename().file() ==
180   - 'irf_file.fits.gz',
181   - 'Check background name')
182   - self.test_assert(rsp.background().filename().extname() ==
183   - 'BACKGROUND',
184   - 'Check background extension name')
  249 + self.test_value(obs[0].eventfile().file(), 'events.fits.gz',
  250 + 'Check event file name')
  251 + self.test_value(obs[0].eventfile().extname(), 'EVENTS',
  252 + 'Check event extension name')
  253 + self.test_value(rsp.aeff().filename().file(), 'irf_file.fits.gz',
  254 + 'Check effective area file name')
  255 + self.test_value(rsp.aeff().filename().extname(), 'EFFECTIVE AREA',
  256 + 'Check effective area extension name')
  257 + self.test_value(rsp.psf().filename().file(), 'irf_file.fits.gz',
  258 + 'Check point spread function file name')
  259 + self.test_value(rsp.psf().filename().extname(), 'POINT SPREAD FUNCTION',
  260 + 'Check point spread function extension name')
  261 + self.test_value(rsp.edisp().filename().file(), 'irf_file.fits.gz',
  262 + 'Check energy dispersion file name')
  263 + self.test_value(rsp.edisp().filename().extname(), 'ENERGY DISPERSION',
  264 + 'Check energy dispersion extension name')
  265 + self.test_value(rsp.background().filename().file(), 'irf_file.fits.gz',
  266 + 'Check background file name')
  267 + self.test_value(rsp.background().filename().extname(), 'BACKGROUND',
  268 + 'Check background extension name')
185 269  
186 270 # Return
187 271 return
... ...