#!/usr/bin/env python # ========================================================================== # Dump information about observation into log file # # Copyright (C) 2015-2016 Michael Mayer # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ========================================================================== import sys import gammalib import ctools # =============== # # csobsinfo class # # =============== # class csobsinfo(ctools.cscript): """ Shows the content of an observation container """ # Constructor def __init__(self, *argv): """ Constructor. """ # Set name self._name = 'csobsinfo' self._version = '1.1.0' # Initialise class members self._obj_dir = None self._compute_offset = False self._offsets = [] self._zeniths = [] self._azimuths = [] self._pnt_ra = [] self._pnt_dec = [] self._ebounds = gammalib.GEbounds() self._gti = gammalib.GGti() # Initialise observation container from constructor arguments. self._obs, argv = self._set_input_obs(argv) # Initialise application by calling the appropriate class # constructor. self._init_cscript(argv) # Return return # Private methods def _get_parameters(self): """ Get parameters from parfile and setup the observation """ # Get parameters if self._obs.size() == 0: self._require_inobs('csobsinfo::get_parameters') self._obs = self._get_observations(False) # Initialise object position self._obj_dir = gammalib.GSkyDir() # Get (optional) offset parameters self._compute_offset = self['offset'].boolean() if self._compute_offset: ra = self['ra'].real() dec = self['dec'].real() self._obj_dir.radec_deg(ra,dec) # Read ahead DS9 filename if self._read_ahead(): self['ds9file'].filename() # Write input parameters into logger if self._logTerse(): self._log_parameters() self._log('\n') # Return return # Public methods def run(self): """ Run the script """ # Switch screen logging on in debug mode if self._logDebug(): self._log.cout(True) # Get parameters self._get_parameters() # Initialise arrays to store certain values for reuse # Todo, think about using a python dictionary self._offsets = [] self._zeniths = [] self._azimuths = [] self._pnt_ra = [] self._pnt_dec = [] self._ebounds = gammalib.GEbounds() self._gti = gammalib.GGti() obs_names = [] # Initialise output to be filled ontime = 0.0 livetime = 0.0 n_events = 0 n_eventbins = 0 n_obs_binned = 0 n_obs_unbinned = 0 # Write header if self._logTerse(): self._log('\n') self._log.header1(gammalib.number('Observation', self._obs.size())) # Loop over observations for obs in self._obs: # Skip non-CTA observations if not obs.classname() == 'GCTAObservation': self._log('Skipping '+obs.instrument()+' observation\n') continue # Use observed object as observation name if name is not given obs_name = obs.name() if obs_name == '': obs_name = obs.object() # Logging if self._logTerse(): self._log.header2(obs_name) # Retrieve observation name obs_names.append(obs_name) # Retrieve energy boundaries obs_bounds = obs.events().ebounds() # Retrieve time interval obs_gti = obs.events().gti() # Compute mean time and dead time fraction in percent deadfrac = (1.0-obs.deadc())*100.0 # Retrieve pointing and store Ra,Dec pnt_dir = obs.pointing().dir() self._pnt_ra.append(pnt_dir.ra_deg()) self._pnt_dec.append(pnt_dir.dec_deg()) # If avaliable append energy boundaries if obs_bounds.size() > 0 : self._ebounds.append(obs_bounds.emin(),obs_bounds.emax()) # Append time interval self._gti.append(obs_gti.tstart(), obs_gti.tstop()) # Increment global livetime and ontime ontime += obs.ontime() livetime += obs.livetime() # Bookkeeping if obs.eventtype() == 'CountsCube': n_eventbins += obs.events().size() n_obs_binned += 1 if self._logTerse(): self._log.parformat('Binned') self._log('yes\n') self._log.parformat('Number of bins') self._log(str(obs.events().size())) self._log('\n') else: n_events += obs.events().size() n_obs_unbinned += 1 if self._logTerse(): self._log.parformat('Binned') self._log('no\n') self._log.parformat('Number of events') self._log(str(obs.events().size())) self._log('\n') # Retrieve zenith and azimuth and store for later use zenith = obs.pointing().zenith() azimuth = obs.pointing().azimuth() self._zeniths.append(zenith) self._azimuths.append(azimuth) # Optionally compute offset with respect to target direction if self._compute_offset: offset = pnt_dir.dist_deg(self._obj_dir) self._offsets.append(offset) else: self._offsets.append(-1.0) # Optionally log details if self._logTerse(): # Log the observation energy range (if available) self._log.parformat('Energy range') if obs_bounds.size() == 0: self._log('undefined') else: self._log(str(obs_bounds.emin())) self._log(' - ') self._log(str(obs_bounds.emax())) self._log('\n') # Log observation time interval self._log.parformat('Time range (MJD)') if obs_gti.size() == 0: self._log('undefined') else: self._log(str(obs_gti.tstart().mjd())) self._log(' - ') self._log(str(obs_gti.tstop().mjd())) self._log('\n') # Log ontime self._log.parformat('Ontime') self._log(str(obs.ontime())) self._log(' s\n') # Log livetime self._log.parformat('Livetime') self._log(str(obs.livetime())) self._log(' s\n') # Log dead time fraction self._log.parformat('Deadtime fraction (%)') self._log('%.3f' % (deadfrac)) self._log('\n') # Log pointing direction self._log.parformat('Pointing') self._log(str(pnt_dir)) self._log('\n') # Optionally log offset with respect to target direction if self._compute_offset: self._log.parformat('Offset from target') self._log('%.2f' % (offset)) self._log(' deg\n') # Log Zenith and Azimuth if required self._log.parformat('Zenith angle') self._log('%.2f' % (zenith)) self._log(' deg\n') self._log.parformat('Azimuth angle') self._log('%.2f' % (azimuth)) self._log(' deg\n') # Log summary if self._logTerse(): # Write header self._log('\n') self._log.header1('Summary') # Log general summary self._log.header3('Observations') self._log.parformat('Unbinned observations') self._log(str(n_obs_unbinned)) self._log('\n') self._log.parformat('Binned observations') self._log(str(n_obs_binned)) self._log('\n') self._log.header3('Events') self._log.parformat('Number of events') self._log(str(n_events)) self._log('\n') self._log.parformat('Number of bins') self._log(str(n_eventbins)) self._log('\n') # Log pointing summary self._log.header3('Pointings') # Log mean offset if possible if self._compute_offset: self._log.parformat('Mean offset angle') self._log('%.2f' % (sum(self._offsets)/len(self._offsets))) self._log(' deg\n') # Log mean azimuth and zenith angle self._log.parformat('Mean zenith angle') self._log('%.2f' % (sum(self._zeniths)/len(self._zeniths))) self._log(' deg\n') self._log.parformat('Mean azimuth angle') self._log('%.2f' % (sum(self._azimuths)/len(self._azimuths))) self._log(' deg\n') # Optionally log names of observations. Note that the set class # is used to extract all different observation names from the # list of observation names, and the set class is only available # from Python 2.4 on. if self._logExplicit() and sys.version_info >= (2,4): obs_set = set(obs_names) for name in obs_set: self._log.parformat('"'+name+'"') self._log(str(obs_names.count(name))) self._log('\n') self._log('\n') # Log energy range self._log.header3('Energy range') self._log.parformat('Minimum energy') if self._ebounds.size() == 0: self._log('undefined') else: self._log(str(self._ebounds.emin())) self._log('\n') self._log.parformat('Maximum energy') if self._ebounds.size() == 0: self._log('undefined') else: self._log(str(self._ebounds.emax())) self._log('\n') # Log time range self._log.header3('Time range') self._log.parformat('Start (MJD)') self._log(str(self._gti.tstart().mjd())) self._log('\n') self._log.parformat('Stop (MJD)') self._log(str(self._gti.tstop().mjd())) self._log('\n') # Log ontime and livetime in different units self._log.parformat('Total ontime') self._log('%.2f s = %.2f min = %.2f h' % (ontime, ontime/60., ontime/3600.)) self._log('\n') self._log.parformat('Total livetime') self._log('%.2f s = %.2f min = %.2f h' % (livetime, livetime/60.0, livetime/3600.)) self._log('\n') # Return return def save(self): """ Save pointings into DS9 region file This method saves all pointing directions that are found in the observation container into a DS9 region file. If "NONE" is specified for the "ds9file" parameter the method does nothing. """ # Get output filename in case it was not read ahead ds9file = self['ds9file'].filename() # Check if DS9 file is valid if ds9file.url() != 'NONE': # Write header if self._logTerse(): self._log('\n') self._log.header1('Save pointings in DS9 file') # Log filename if self._logTerse(): self._log(gammalib.parformat('DS9 filename')) self._log(ds9file.url()) self._log('\n') # Open file f = open(ds9file.url(),'w') # Write coordinate system f.write('fk5\n') # Loop over pointings for i in range(len(self._pnt_ra)): # Create string line = 'point(' line += str(self._pnt_ra[i])+','+str(self._pnt_dec[i])+')' line += ' # point=cross 20 width=3\n' # Write to file f.write(line) # Close file f.close() # Return return def zeniths(self): """ Return zenith angles """ return self._zeniths def azimuths(self): """ Return azimuth angles """ return self._azimuths def offsets(self): """ Return offset angles """ return self._offsets def ebounds(self): """ Return energy boundaries """ return self._ebounds def gti(self): """ Return good time intervals """ return self._gti def execute(self): """ Execute the script. """ # Open logfile self.logFileOpen() # Read ahead output parameters self._read_ahead(True) # Run the script self.run() # Save ds9 file if required self.save() # Return return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': # Create instance of application app = csobsinfo(sys.argv) # Execute application app.execute()