diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 6826fc9..7fbfd42 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -11,9 +11,10 @@ pythonscripts = __init__.py \ # Parfiles to be distributed in $(prefix)/syspfiles parfiles = cspull.par cssens.par cstsdist.par cscaldb.par csresmap.par \ - csobsdef.par csspec.par csobsinfo.par \ + csobsdef.par csspec.par cslightcrv.par csobsinfo.par \ csmodelinfo.par csmodelmerge.par cstsmapmerge.par \ - csobs2caldb.par + csobs2caldb.par \ + csfindobs.par csiactcopy.par csiactdata.par csiactobs.par # Scripts for which symbolic links will be created in bindir linkscripts = cspull cssens cstsdist cscaldb csresmap \ diff --git a/cscripts/cscaldb.py b/cscripts/cscaldb.py index 04f3683..5f3f589 100755 --- a/cscripts/cscaldb.py +++ b/cscripts/cscaldb.py @@ -37,7 +37,7 @@ class cscaldb(ctools.cscript): into the "cscaldb.log" log file. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -46,28 +46,13 @@ class cscaldb(ctools.cscript): self._name = "cscaldb" self._version = "1.1.0" - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - # Private methods def _get_parameters(self): @@ -164,16 +149,6 @@ class cscaldb(ctools.cscript): # Public methods - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Return - return - def run(self): """ Run the script. @@ -238,6 +213,19 @@ class cscaldb(ctools.cscript): # Return return + def execute(self): + """ + Execute the script. + """ + # Open logfile + self._logFileOpen() + + # Run the script + self.run() + + # Return + return + # ======================== # # Main routine entry point # @@ -247,8 +235,5 @@ if __name__ == '__main__': # Create instance of application app = cscaldb(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csfindobs.par b/cscripts/csfindobs.par new file mode 100644 index 0000000..a4767d6 --- /dev/null +++ b/cscripts/csfindobs.par @@ -0,0 +1,44 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csfindobs.par # +# Version: 1.1.0 # +# # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # +# # +# Purpose: Parameter file for the csfindobs script # +# # +############################################################################# +# +# Input and output parameters +#============================ +datapath, s, a, $VHEFITS,,, "Path were data are located" +prodname, s, a, prod-name,,, "Name of FITS production (Run csiactdata to view your options)" +outfile, f, a, runlist.lis,,, "Output runlist file" + +# +# Script parameters +#================== +ra, r, a, 83.63,,, "Right ascension of selection region centre (deg)" +dec, r, a, 22.01,,, "Declination of selection region centre (deg)" +rad, r, a, 2.5,,, "Search radius of selection region (deg)" +min_qual, i, h, 0,0|1|2,, "Minimum data quality (0=perfect, 1=ok, 2=bad)" +expression, s, h, NONE,,, "Additional expression" +master_indx, s, h, master.json,,, "Name of master index file" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, csfindobs.log,,, "Log filename" + + + diff --git a/cscripts/csfindobs.py b/cscripts/csfindobs.py index e163020..086774d 100755 --- a/cscripts/csfindobs.py +++ b/cscripts/csfindobs.py @@ -33,7 +33,7 @@ class csfindobs(ctools.cscript): Find observations from an IACT data store. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -44,68 +44,15 @@ class csfindobs(ctools.cscript): self._verbose = False self._datapath = os.getenv("VHEFITS","") - # Make sure that parfile exists - self._parfile() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Write separator into logger - if self._logTerse(): - self._log("\n") - - # Return - return - def _parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a - parfile. - """ - # Set parfile name - parfile = self._name+".par" - - try: - pars = gammalib.GApplicationPars(parfile) - except: - # Signal if parfile was not found - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("datapath","s","a",self._datapath,"","","Path were data are located")) - pars.append(gammalib.GApplicationPar("prodname","s","a","prod-name","","","Name of FITS production (Run csiactdata to view your options)")) - pars.append(gammalib.GApplicationPar("outfile","f","a","runlist.lis","","","Runlist outfile")) - pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","","","Right ascension")) - pars.append(gammalib.GApplicationPar("dec","r","a","22.01","","","Declination")) - pars.append(gammalib.GApplicationPar("rad","r","a","2.5","","","Search radius")) - pars.append(gammalib.GApplicationPar("min_qual","i","h","0","0|1|2","","Minimum data quality (0=perfect, 1=ok, 2=bad)")) - pars.append(gammalib.GApplicationPar("expression","s","h","NONE","","","Additional expression")) - pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csfindobs.log","","","Log filename")) - pars.save(parfile) - - # Return - return - + # Private methods def _get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -209,27 +156,7 @@ class csfindobs(ctools.cscript): # Return return - def obs_ids(self): - """ - Return OBS IDs - """ - - # Return - return self.runs - - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Save residual map - self.save(self._outfile) - - # Return - return - + # Public methods def run(self): """ Run the script. @@ -304,6 +231,22 @@ class csfindobs(ctools.cscript): # Return return + def execute(self): + """ + Execute the script. + """ + # Open logfile + self._logFileOpen() + + # Run the script + self.run() + + # Save residual map + self.save(self._outfile) + + # Return + return + def save(self, outfile): """ Save. @@ -323,6 +266,14 @@ class csfindobs(ctools.cscript): # Return return + def obs_ids(self): + """ + Return OBS IDs + """ + + # Return + return self._runs + # ======================== # # Main routine entry point # @@ -332,8 +283,5 @@ if __name__ == '__main__': # Create instance of application app = csfindobs(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csiactcopy.par b/cscripts/csiactcopy.par new file mode 100644 index 0000000..a9257a4 --- /dev/null +++ b/cscripts/csiactcopy.par @@ -0,0 +1,32 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csiactcopy.par # +# Version: 1.1.0 # +# # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # +# # +# Purpose: Parameter file for the csiactcopy script # +# # +############################################################################# +# +# Input and output parameters +#============================ +remote_master, f, a, /Volumes/mountpoint/master.json,,, "Location of remote master file" +prodname, s, a, iact-fits,,, "Name of FITS production to download" +runlist, f, h, NONE,,, "Input list of observation IDs" +outpath, s, a, /path/to/fits,,, "Destination path of FITS data" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, csiactcopy.log,,, "Log filename" diff --git a/cscripts/csiactcopy.py b/cscripts/csiactcopy.py index 5ea41de..327f5e0 100755 --- a/cscripts/csiactcopy.py +++ b/cscripts/csiactcopy.py @@ -34,6 +34,8 @@ class csiactcopy(ctools.cscript): It can take a list of observation IDs to allow the download specific observations. Index files get merged and updated accordingly. """ + + # Constructor def __init__(self, *argv): """ Constructor. @@ -43,60 +45,15 @@ class csiactcopy(ctools.cscript): self._version = "1.1.0" self._datapath = os.getenv("VHEFITS","") - # Make sure that parfile exists - self._parfile() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - - # Return - return - - def _parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a - parfile. - """ - # Set parfile name - parfile = self._name+".par" - - try: - pars = gammalib.GApplicationPars(parfile) - except: - # Signal if parfile was not found - print("Parfile "+parfile+" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("remote_master","f","a","/Volumes/mountpoint/master.json","","","Location of remote master file")) - pars.append(gammalib.GApplicationPar("prodname","s","a","iact-fits","","","Name of FITS production to download")) - pars.append(gammalib.GApplicationPar("outpath","s","a","/path/to/fits","","","Destination path of FITS data")) - pars.append(gammalib.GApplicationPar("runlist","f","h","NONE","","","List of observation IDs")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csiactcopy.log","","","Log filename")) - pars.save(parfile) - - # Return - return + # Private methods def _get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -114,20 +71,233 @@ class csiactcopy(ctools.cscript): self._prodname = self["prodname"].string() self._outpath = gammalib.expand_env(self["outpath"].string()) self._runlist = self["runlist"].filename() - + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") + # Return return - def execute(self): + def copy(self, source, clobber): """ - Execute the script. + copy file to _outpath """ - # Run the script - self.run() + # Get file destination + destination = os.path.join(self._outpath, + os.path.relpath(source, self._remote_base)) + + # Initialise return value + filesize = 0.0 + + # Logging + if self._logVerbose(): + self._log("\n") + self._log.header3("Copy file") + self._log.parformat("Source") + self._log(str(source)) + self._log("\n") + self._log.parformat("Destination") + self._log(str(destination)) + self._log("\n") + self._log.parformat("Already exists") + self._log(str(os.path.isfile(destination))) + self._log("\n") + self._log.parformat("Overwrite") + self._log(str(clobber)) + self._log("\n") + self._log.parformat("Copying") + + # Flag if file destination is already available + is_file = os.path.isfile(destination) + + # check if file could be skipped because clobber=no + if is_file and clobber == False: + if self._logVerbose(): + self._log("Skip (clobber=no)") + self._log("\n") + + # check if file could be skipped because it is the same file + elif is_file and os.path.samefile(destination, source): + if self._logVerbose(): + self._log("Skip (same file)") + self._log("\n") + + else: + # Create directory if not existent + dest_dir = os.path.dirname(destination) + if not os.path.isdir(dest_dir): + os.makedirs(dest_dir) + + # Copy file (if source not identical to destination) + shutil.copy2(source, destination) + + # Get filesize + filesize = os.stat(destination).st_size + + # Logging + if self._logVerbose(): + self._log("Done!") + self._log("\n") + + # Return + return filesize + + + def merge(self, localfits, remotefits, hduname, clobber): + """ + merge remote and local fits file + If local fits file not present, a new one is created + """ + + if self._logTerse(): + self._log.parformat("Local file") + self._log(str(localfits)) + self._log("\n") + self._log.parformat("Remote file") + self._log(str(remotefits)) + self._log("\n") + self._log.parformat("HDU name") + self._log(str(hduname)) + self._log("\n") + self._log.parformat("Conflict preference") + if clobber: + self._log("Remote") + else: + self._log("Local") + self._log("\n") + + # Check if local fits file is available + if os.path.isfile(localfits): + local = gammalib.GFits(str(localfits)) + local_hdu = local[hduname] + else: + # Otherwise load remote fits file and delete rows + local = gammalib.GFits(str(remotefits)) + local_hdu = local[hduname] + local_hdu.remove_rows(0,local_hdu.nrows()) + + # find local obs_id + local_obs = [] + lobs_col = local_hdu["OBS_ID"] + for obs_id in lobs_col: + local_obs.append(obs_id) + + # Create selection expression for opening remote fits file + selection = "" + if len(self._runs): + selection += "["+hduname+"][" + for run in self._runs: + selection += "OBS_ID=="+str(run) + selection += "||" + selection = selection[:-2] + selection += "]" + + # Create file name + remotefile = remotefits + selection + + # open remote fits file + remote = gammalib.GFits(str(remotefile)) + remote_hdu = remote[hduname] + + # Get remote obs_id + remote_obs = [] + robs_col = remote_hdu["OBS_ID"] + for obs_id in robs_col: + remote_obs.append(obs_id) + + if self._logTerse(): + self._log.parformat("Remote entries") + self._log(str(len(remote_obs))) + self._log("\n") + self._log.parformat("Local entries") + self._log(str(len(local_obs))) + self._log("\n") + + # initialise counter for logging + removed_rows = 0 + + # If clobber=yes, we use the remote content to overwrite + # corresponding local content + if clobber: + + # Loop over remote obs_ids + for remote_obs_id in remote_obs: + + # Check if remote obs_id is present locally + if remote_obs_id in local_obs: + + # Remove local obs_id entries + table_has_obsid = True + while table_has_obsid: + for i in range(local_hdu.nrows()): + if remote_obs_id == local_hdu["OBS_ID"][i]: + local_hdu.remove_rows(i,1) + removed_rows += 1 + break + table_has_obsid = False + for i in range(local_hdu.nrows()): + if remote_obs_id == local_hdu["OBS_ID"][i]: + table_has_obsid = True + break + + # If clobber=false, we keep the local content unchanged and just + # expand by remote content + else: + # Loop over remote obs_ids + for local_obs_id in local_obs: + + # Check if local obs_id is present remotely + if local_obs_id in remote_obs: + + # Remove remote obs_id entries + table_has_obsid = True + while table_has_obsid: + for i in range(remote_hdu.nrows()): + if local_obs_id == remote_hdu["OBS_ID"][i]: + remote_hdu.remove_rows(i,1) + removed_rows += 1 + break + table_has_obsid = False + for i in range(remote_hdu.nrows()): + if local_obs_id == remote_hdu["OBS_ID"][i]: + table_has_obsid = True + break + + # Store number of local rows + old_local_rows = local_hdu.nrows() + + # Add rmeote HDUs + local_hdu.insert_rows(old_local_rows, remote_hdu.nrows()) + + if self._logTerse(): + if clobber: + self._log.parformat("Removed local rows") + else: + self._log.parformat("Skipped remote rows") + self._log(str(removed_rows)) + self._log("\n") + # Loop over columns + for i in range(local_hdu.ncols()): + + # Get local and remote columns + local_col = local_hdu[i] + remote_col = remote_hdu[i] + + # Loop over entries and merge + for j in range(remote_col.length()): + local_col[j+old_local_rows] = remote_col[j] + + # Save local fits file + local.saveto(str(localfits), True) + # Return return - + + + # Public methods def run(self): """ Run the script. @@ -138,11 +308,6 @@ class csiactcopy(ctools.cscript): # Get parameters self._get_parameters() - # Write input parameters into logger - if self._logTerse(): - self._log_parameters() - self._log("\n") - # Make destination directory if not available if not os.path.isdir(self._outpath): os.makedirs(self._outpath) @@ -516,221 +681,19 @@ class csiactcopy(ctools.cscript): # Return return - def copy(self, source, clobber): - """ - copy file to m_outpath + def execute(self): """ - # Get file destination - destination = os.path.join(self._outpath, - os.path.relpath(source, self._remote_base)) - - # Initialise return value - filesize = 0.0 - - # Logging - if self._logVerbose(): - self._log("\n") - self._log.header3("Copy file") - self._log.parformat("Source") - self._log(str(source)) - self._log("\n") - self._log.parformat("Destination") - self._log(str(destination)) - self._log("\n") - self._log.parformat("Already exists") - self._log(str(os.path.isfile(destination))) - self._log("\n") - self._log.parformat("Overwrite") - self._log(str(clobber)) - self._log("\n") - self._log.parformat("Copying") - - # Flag if file destination is already available - is_file = os.path.isfile(destination) - - # check if file could be skipped because clobber=no - if is_file and clobber == False: - if self._logVerbose(): - self._log("Skip (clobber=no)") - self._log("\n") - - # check if file could be skipped because it is the same file - elif is_file and os.path.samefile(destination, source): - if self._logVerbose(): - self._log("Skip (same file)") - self._log("\n") - - else: - # Create directory if not existent - dest_dir = os.path.dirname(destination) - if not os.path.isdir(dest_dir): - os.makedirs(dest_dir) - - # Copy file (if source not identical to destination) - shutil.copy2(source, destination) - - # Get filesize - filesize = os.stat(destination).st_size - - # Logging - if self._logVerbose(): - self._log("Done!") - self._log("\n") - - # Return - return filesize - - - def merge(self, localfits, remotefits, hduname, clobber): + Execute the script. """ - merge remote and local fits file - If local fits file not present, a new one is created - """ - - if self._logTerse(): - self._log.parformat("Local file") - self._log(str(localfits)) - self._log("\n") - self._log.parformat("Remote file") - self._log(str(remotefits)) - self._log("\n") - self._log.parformat("HDU name") - self._log(str(hduname)) - self._log("\n") - self._log.parformat("Conflict preference") - if clobber: - self._log("Remote") - else: - self._log("Local") - self._log("\n") - - # Check if local fits file is available - if os.path.isfile(localfits): - local = gammalib.GFits(str(localfits)) - local_hdu = local[hduname] - else: - # Otherwise load remote fits file and delete rows - local = gammalib.GFits(str(remotefits)) - local_hdu = local[hduname] - local_hdu.remove_rows(0,local_hdu.nrows()) - - # find local obs_id - local_obs = [] - lobs_col = local_hdu["OBS_ID"] - for obs_id in lobs_col: - local_obs.append(obs_id) - - # Create selection expression for opening remote fits file - selection = "" - if len(self._runs): - selection += "["+hduname+"][" - for run in self._runs: - selection += "OBS_ID=="+str(run) - selection += "||" - selection = selection[:-2] - selection += "]" - - # Create file name - remotefile = remotefits + selection - - # open remote fits file - remote = gammalib.GFits(str(remotefile)) - remote_hdu = remote[hduname] - - # Get remote obs_id - remote_obs = [] - robs_col = remote_hdu["OBS_ID"] - for obs_id in robs_col: - remote_obs.append(obs_id) - - if self._logTerse(): - self._log.parformat("Remote entries") - self._log(str(len(remote_obs))) - self._log("\n") - self._log.parformat("Local entries") - self._log(str(len(local_obs))) - self._log("\n") - - # initialise counter for logging - removed_rows = 0 - - # If clobber=yes, we use the remote content to overwrite - # corresponding local content - if clobber: - - # Loop over remote obs_ids - for remote_obs_id in remote_obs: - - # Check if remote obs_id is present locally - if remote_obs_id in local_obs: - - # Remove local obs_id entries - table_has_obsid = True - while table_has_obsid: - for i in range(local_hdu.nrows()): - if remote_obs_id == local_hdu["OBS_ID"][i]: - local_hdu.remove_rows(i,1) - removed_rows += 1 - break - table_has_obsid = False - for i in range(local_hdu.nrows()): - if remote_obs_id == local_hdu["OBS_ID"][i]: - table_has_obsid = True - break - - # If clobber=false, we keep the local content unchanged and just - # expand by remote content - else: - # Loop over remote obs_ids - for local_obs_id in local_obs: - - # Check if local obs_id is present remotely - if local_obs_id in remote_obs: - - # Remove remote obs_id entries - table_has_obsid = True - while table_has_obsid: - for i in range(remote_hdu.nrows()): - if local_obs_id == remote_hdu["OBS_ID"][i]: - remote_hdu.remove_rows(i,1) - removed_rows += 1 - break - table_has_obsid = False - for i in range(remote_hdu.nrows()): - if local_obs_id == remote_hdu["OBS_ID"][i]: - table_has_obsid = True - break - - # Store number of local rows - old_local_rows = local_hdu.nrows() - - # Add rmeote HDUs - local_hdu.insert_rows(old_local_rows, remote_hdu.nrows()) - - if self._logTerse(): - if clobber: - self._log.parformat("Removed local rows") - else: - self._log.parformat("Skipped remote rows") - self._log(str(removed_rows)) - self._log("\n") + # Open logfile + self._logFileOpen() + + # Run the script + self.run() - # Loop over columns - for i in range(local_hdu.ncols()): - - # Get local and remote columns - local_col = local_hdu[i] - remote_col = remote_hdu[i] - - # Loop over entries and merge - for j in range(remote_col.length()): - local_col[j+old_local_rows] = remote_col[j] - - # Save local fits file - local.saveto(str(localfits), True) - # Return return + # ======================== # # Main routine entry point # @@ -740,8 +703,5 @@ if __name__ == '__main__': # Create instance of application app = csiactcopy(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csiactdata.par b/cscripts/csiactdata.par new file mode 100644 index 0000000..b13434f --- /dev/null +++ b/cscripts/csiactdata.par @@ -0,0 +1,31 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csiactdata.par # +# Version: 1.1.0 # +# # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # +# # +# Purpose: Parameter file for the csiactdata script # +# # +############################################################################# +# +# Input and output parameters +#============================ +datapath, s, a, $VHEFITS,,, "Path were data is located" +master_indx, s, h, master.json,,, "Name of master index file" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, csiactdata.log,,, "Log filename" + diff --git a/cscripts/csiactdata.py b/cscripts/csiactdata.py index c8524e3..3ebeb18 100755 --- a/cscripts/csiactdata.py +++ b/cscripts/csiactdata.py @@ -2,7 +2,7 @@ # ========================================================================== # Inspection of IACT data storage # -# Copyright (C) 2015 Michael Mayer +# 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 @@ -35,6 +35,8 @@ class csiactdata(ctools.cscript): file or prints it on the screen. These information can be used as input for 'csiactobs'. """ + + # Constructor def __init__(self, *argv): """ Constructor. @@ -45,57 +47,15 @@ class csiactdata(ctools.cscript): self._datapath = os.getenv("VHEFITS","") self._prodnames = [] - # Make sure that parfile exists - file = self._parfile() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - - def _parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a - parfile. - """ - # Set parfile name - parfile = self._name+".par" - - try: - pars = gammalib.GApplicationPars(parfile) - except: - # Signal if parfile was not found - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("datapath","s","a",self._datapath,"","","Path were data is located")) - pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csiactdata.log","","","Log filename")) - pars.save(parfile) - - # Return - return + # Private methods def get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -119,16 +79,7 @@ class csiactdata(ctools.cscript): # Return return - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Return - return - + # Public methods def run(self): """ Run the script. @@ -222,6 +173,20 @@ class csiactdata(ctools.cscript): # Return return + def execute(self): + """ + Execute the script. + """ + # Open logfile + self._logFileOpen() + + # Run the script + self.run() + + # Return + return + + def names(self): """ Return available FITS production names @@ -239,8 +204,5 @@ if __name__ == '__main__': # Create instance of application app = csiactdata(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csiactobs.par b/cscripts/csiactobs.par new file mode 100644 index 0000000..d6e2aa5 --- /dev/null +++ b/cscripts/csiactobs.par @@ -0,0 +1,54 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csiactobs.par # +# Version: 1.1.0 # +# # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # +# # +# Purpose: Parameter file for the csiactobs script # +# # +############################################################################# +# +# Input and output parameters +#============================ +datapath, s, a, $VHEFITS,,, "Path were data is located" +prodname, s, a, ,,, "Data storage name" +infile, f, a, runlist.lis,,, "Input runlist file" +inmodel, f, h, NONE,,, "Input model definition XML file (optional)" +master_indx, s, h, master.json,,, "Input master index file" +outobs, f, a, obs.xml,,, "Output observation definition XML file" +outmodel, f, a, bgmodels.xml,,, "Output model definition XML file" + +# +# Script parameters +#================== +bkgpars, i, a, 1,,, "Number of free parameters per background model" +bkg_scale, b, h, yes,,, "Apply background scaling factor from observation index file if available?" +ev_hiera, s, h, events,,, "Hierarchy of event formats" +aeff_hiera, s, h, aeff_2d,,, "Hierarchy of effective area formats" +psf_hiera, s, h, psf_king|psf_3gauss,,, "Hierarchy of psf formats" +edisp_hiera, s, h, edisp_2d,,, "Hierarchy of energy dispersion formats" +bkg_hiera, s, h, bkg_3d,,, "Hierarchy of background formats" +bkg_mod_hiera, s, h, irf|aeff|gauss,,, "Hierarchy of background models" +bkg_gauss_norm, r, h, 1e-8,,, "Input normalisation for Gaussian background" +bkg_gauss_index, r, h, -2.0,,, "Input spectral index for Gaussian background" +bkg_gauss_sigma, r, h, 2.5,,, "Input sigma for Gaussian background" +bkg_aeff_norm, r, h, 1e-14,,, "Input normalisation for effective area background" +bkg_aeff_index, r, h, -2.0,,, "Input spectral index for effective area background" +bkg_range_factor, r, h, 100.0,,, "Factor to determine range of background normalisation" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, csiactobs.log,,, "Log filename" + diff --git a/cscripts/csiactobs.py b/cscripts/csiactobs.py index f4ae29c..3c74c5c 100755 --- a/cscripts/csiactobs.py +++ b/cscripts/csiactobs.py @@ -37,6 +37,8 @@ class csiactobs(ctools.cscript): described here: http://gamma-astro-data-formats.readthedocs.org/en/latest/ """ + + # Constructor def __init__(self, *argv): """ Constructor. @@ -53,75 +55,15 @@ class csiactobs(ctools.cscript): self._models = gammalib.GModels() self._runlist = [] - # Make sure that parfile exists - file = self._parfile() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) - - # Return - return + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - def __del__(self): - """ - Destructor. - """ # Return return - def _parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a parfile. - """ - # Set parfile name - parfile = self._name+".par" - - try: - pars = gammalib.GApplicationPars(parfile) - except: - # Signal if parfile was not found - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("datapath","s","a",self._datapath,"","","Path were data is located")) - pars.append(gammalib.GApplicationPar("prodname","s","a","","","","Data storage name (Run csiactdata to view your options)")) - pars.append(gammalib.GApplicationPar("infile","f","a","runlist.lis","","","Runlist file")) - pars.append(gammalib.GApplicationPar("inmodel","f","h","NONE","","","Input model XML file (optional)")) - pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Observation XML outfile")) - pars.append(gammalib.GApplicationPar("outmodel","f","a","bgmodels.xml","","","Output model XML file")) - pars.append(gammalib.GApplicationPar("bkgpars","i","a","1","","","Number of free parameters per background model")) - pars.append(gammalib.GApplicationPar("master_indx","s","h","master.json","","","Name of master index file")) - pars.append(gammalib.GApplicationPar("bkg_scale","b","h","yes","","","Specifies whether the background scaling factor from the observation index file should be applied if available. ")) - pars.append(gammalib.GApplicationPar("ev_hiera","s","h","events","","","Hierarchy of event formats")) - pars.append(gammalib.GApplicationPar("aeff_hiera","s","h","aeff_2d","","","Hierarchy of effective area formats")) - pars.append(gammalib.GApplicationPar("psf_hiera","s","h","psf_king|psf_3gauss","","","Hierarchy of psf formats")) - pars.append(gammalib.GApplicationPar("edisp_hiera","s","h","edisp_2d","","","Hierarchy of energy dispersion formats")) - pars.append(gammalib.GApplicationPar("bkg_hiera","s","h","bkg_3d","","","Hierarchy of background formats")) - pars.append(gammalib.GApplicationPar("bkg_mod_hiera","s","h","irf|aeff|gauss","","","Hierarchy of background models")) - pars.append(gammalib.GApplicationPar("bkg_gauss_norm","r","h","1e-8","","","Input normalisation for Gaussian background")) - pars.append(gammalib.GApplicationPar("bkg_gauss_index","r","h","-2.0","","","Input spectral index for Gaussian background")) - pars.append(gammalib.GApplicationPar("bkg_gauss_sigma","r","h","2.5","","","Input sigma for Gaussian background")) - pars.append(gammalib.GApplicationPar("bkg_aeff_norm","r","h","1e-14","","","Input normalisation for effective area background")) - pars.append(gammalib.GApplicationPar("bkg_aeff_index","r","h","-2.0","","","Input spectral index for effective area background")) - pars.append(gammalib.GApplicationPar("bkg_range_factor","r","h","100.0","","","Factor to determine range of background normalisation")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csiactobs.log","","","Log filename")) - pars.save(parfile) - - # Return - return + # Private methods def _get_parameters(self): """ Get parameters from parfile and setup the observation. @@ -366,43 +308,8 @@ class csiactobs(ctools.cscript): # Return model return model - def ebounds(self): - """ - Returns runlist energy range - """ - return self._ebounds - - def obs(self): - """ - Returns GObservations object - """ - # Initialise observations - obs = gammalib.GObservations() - - # Read observations if XML is filled - if self._xml.size(): - obs.read(self._xml) - - # Assign models - obs.models(self._models) - - # Return observations - return obs - - - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Save residual map - self.save(self._outobs, self._clobber) - - # Return - return + # Public methods def run(self): """ Run the script. @@ -665,6 +572,22 @@ class csiactobs(ctools.cscript): # Return return + def execute(self): + """ + Execute the script. + """ + # Open logfile + self._logFileOpen() + + # Run the script + self.run() + + # Save residual map + self.save(self._outobs, self._clobber) + + # Return + return + def save(self,outfile,clobber): # Save observation XML file @@ -676,6 +599,29 @@ class csiactobs(ctools.cscript): # Return return + def ebounds(self): + """ + Returns runlist energy range + """ + return self._ebounds + + def obs(self): + """ + Returns GObservations object + """ + # Initialise observations + obs = gammalib.GObservations() + + # Read observations if XML is filled + if self._xml.size(): + obs.read(self._xml) + + # Assign models + obs.models(self._models) + + # Return observations + return obs + # ======================== # # Main routine entry point # @@ -685,8 +631,5 @@ if __name__ == '__main__': # Create instance of application app = csiactobs(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/cslightcrv.par b/cscripts/cslightcrv.par new file mode 100644 index 0000000..7a3e8e9 --- /dev/null +++ b/cscripts/cslightcrv.par @@ -0,0 +1,65 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cslightcrv.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the cslightcrv script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inobs, f, h, events.fits,,, "Input event list, counts cube, or observation definition XML file" +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file" +srcname, s, a, Crab,,, "Source name" +expcube, f, a, NONE,,, "Input exposure cube file (only needed for stacked analysis)" +psfcube, f, a, NONE,,, "Input PSF cube file (only needed for stacked analysis)" +edispcube, f, a, NONE,,, "Input energy dispersion cube file (only needed for stacked analysis)" +bkgcube, f, a, NONE,,, "Input background cube file (only needed for stacked analysis)" +caldb, s, a, prod2,,, "Calibration database" +irf, s, a, South_0.5h,,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" +outfile, f, a, lightcurve.fits,,, "Output light curve file" + +# +# Script parameters +#================== +tbinalg, s, a, GTI,FILE|LIN|GTI,, "Algorithm for defining time bins" +tmin, r, a, 51544.5,,,"Lightcurve start time (MJD])" +tmax, r, a, 51544.6,,,"Lightcurve stop time (MJD)" +tbins, i, a, 5,1,10000, "Number of time bins" +tbinfile, f, a, NONE,,, "File defining the time binning" +emin, r, a, 0.1,,, "Lower energy limit of events (TeV)" +emax, r, a, 100.0,,, "Upper energy limit of events (TeV)" +enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)" +coordsys, s, a, CEL,CEL|GAL,, "Coordinate System" +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" +nxpix, i, a, 200,,, "Size of the X axis in pixels" +nypix, i, a, 200,,, "Size of the Y axis in pixels" +binsz, r, a, 0.02,,, "Pixel size (deg/pixel)" +calc_ts, b, h, yes,,, "Compute TS value for each bin?" +calc_ulim, b, h, yes,,, "Compute upper limit for each bin?" +fix_srcs, b, h, yes,,, "Fix other sky model parameters?" +fix_bkg, b, h, no,,, "Fix background model parameters?" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, cslightcrv.log,,, "Log filename" + + + diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index b90268c..313bc55 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -30,159 +30,70 @@ class cslightcrv(ctools.cscript): """ Generates a lightcurve. - This class implements the creation of a light curve. It derives from - the ctools.cscript class which provides support for parameter files, - command line arguments, and logging. In that way the Python script - behaves just as a regular ctool. + The cslightcrv class generates a light curve for Imaging Air Cherenkov + Telescope event data by performing a maximum likelihood fit using + ctlike in a series of time bins. The time bins can be either + specified in an ASCII file, as an interval divided into equally + sized time bins, or can be taken from the Good Time Intervals of the + observation(s). + + The format of the ASCII file is one row per time bin, each specifying + the start of stop value of the bin, separated by a whitespace. The + times are given in Modified Julian Days (MJD). + + Examples: + >>> lcrv = cslightcrv() + >>> lcrv.run() + >>> ... (querying for parameters) ... + >>> fits = lcrv.lightcurve() + Generates a light curve and retrieves the results in + a FITS file. + + >>> lcrv = cslightcrv() + >>> lcrv.execute() + >>> ... (querying for parameters) ... + Generates a light curve and saves results in a FITS file. + + >>> lcrv = cslightcrv(obs) + >>> lcrv.execute() + >>> ... (querying for parameters) ... + Generates a light curve from the observations in an + observation container and saves results in a FITS file. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. """ - # Set name self._name = "cslightcrv" self._version = "1.1.0" # Initialise some members - self._outfile = "" - - # Initialise observation - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - self._obs = argv[0] - argv = argv[1:] - else: - self._obs = gammalib.GObservations() - - # Make sure that parfile exists - self._parfile() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + self._srcname = "" + self._tbins = gammalib.GGti() + self._stacked = False - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise observation container from constructor arguments. + self._obs = self._set_input_obs(argv) - # Return - return + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - def __del__(self): - """ - Destructor. - """ # Return return # Private methods - def _parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a parfile. - """ - # Set parfile name - parfile = self._name+".par" - - try: - pars = gammalib.GApplicationPars(parfile) - except: - # Signal if parfile was not found - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Input event list, counts cube, or observation definition XML file")) - pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Input model XML file")) - pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name")) - pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","Input PSF cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("edispcube","f","a","NONE","","","Input energy dispersion cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Input background cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) - pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) - pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) - pars.append(gammalib.GApplicationPar("outfile","f","a","lightcurve.fits","","","Output light curve file")) - pars.append(gammalib.GApplicationPar("tbinalg","s","a","GTI","FILE|LIN|GTI","", "Algorithm for defining time bins")) - pars.append(gammalib.GApplicationPar("tmin","r","a","51544.5","","days","Lightcurve start time [MJD]")) - pars.append(gammalib.GApplicationPar("tmax","r","a","51544.6","","days","Lightcurve stop time [MJD]")) - pars.append(gammalib.GApplicationPar("tbins","i","a","5","","","Number of time bins")) - pars.append(gammalib.GApplicationPar("tbinfile","f","a","tbins.fits","","", "File defining the time binning")) - pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins per light curve bin (0=unbinned)")) - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit of events (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit of events (TeV)")) - pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System")) - pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method")) - pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)")) - pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)")) - pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels")) - pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels")) - pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)")) - pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value for each bin?")) - pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit for each bin?")) - pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other sky model parameters?")) - pars.append(gammalib.GApplicationPar("fix_bkg","b","h","no","yes|no","","Fix background model parameters?")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","cslightcrv.log","","","Log filename")) - pars.save(parfile) - - # Return - return - - def _create_tbounds(self): - - # Create time bin container - self._tbins = gammalib.GGti() - - # check for binning algorithm - if self["tbinalg"].string() == "FILE": - try: - self._tbins.load(self["tbinfile"].filename()) - except: - csv = gammalib.GCsv(self["tbinfile"].filename()) - for i in range(csv.nrows()): - tmin = gammalib.GTime() - tmax = gammalib.GTime() - tmin.mjd(csv.real(i,0)) - tmax.mjd(csv.real(i,1)) - self._tbins.append(tmin,tmax) - elif self["tbinalg"].string() == "LIN": - # Use linear time binning - time_min = self["tmin"].real() - time_max = self["tmax"].real() - nbins = self["tbins"].integer() - step = (time_max - time_min) / float(nbins) - for i in range(nbins): - tmin = gammalib.GTime() - tmin.mjd(time_min + i*step)# ref) - tmax = gammalib.GTime() - tmax.mjd(time_min + (i+1)*step)#ref) - self._tbins.append(tmin,tmax) - elif self["tbinalg"].string() == "GTI": - # Use GTIs of observations - for obs in self._obs: - for i in range(obs.events().gti().size()): - self._tbins.append(obs.events().gti().tstart(i),obs.events().gti().tstop(i)) - else: - raise AttributeError("tbinalg=\""+self["tbinalg"].string()+"\" unkown. Must be one of \"FILE\", \"LIN\" or \"GTI\"") - - # Return - return - def _get_parameters(self): """ - Get parameters from parfile and setup the observation. + Get parameters from parfile. """ # Set observation if not done before if self._obs == None or self._obs.size() == 0: - self._require_inobs("cslightcrv::get_parameters()") + self._require_inobs("cslightcrv::_get_parameters()") self._obs = self._get_observations() # Set models if we have none @@ -192,46 +103,35 @@ class cslightcrv(ctools.cscript): # Get source name self._srcname = self["srcname"].string() - # Get energy dispersion flag - self._edisp = self["edisp"].boolean() - # Get time boundaries - self._create_tbounds() + self._tbins = self._create_tbounds() + + # Set stacked analysis flag if the requested number of energy bins + # is not zero; otherwise an unbinned analysis will be done. + self._stacked = False if self["enumbins"].integer() == 0 else True + + # Make sure that other parameters are queried now and not later + self["emin"].real() + self["emax"].real() + if self._stacked: + self["coordsys"].string() + self["proj"].string() + self["xref"].real() + self["yref"].real() + self["nxpix"].integer() + self["nypix"].integer() + self["binsz"].real() + + # Do the same for the hidden parameters, just in case + self["edisp"].boolean() + self["calc_ulim"].boolean() + self["calc_ts"].boolean() + self["fix_bkg"].boolean() + self["fix_srcs"].boolean() - # Unbinned or binned analysis? - self._ebins = self["enumbins"].integer() - if self._ebins == 0: - self._binned = False - else: - self._binned = True - - # Get energy range - self._emin = self["emin"].real() - self._emax = self["emax"].real() - - # Get binning flag - if self._binned: - self._coordsys = self["coordsys"].string() - self._proj = self["proj"].string() - self._xref = self["xref"].real() - self._yref = self["yref"].real() - self._nxpix = self["nxpix"].integer() - self._nypix = self["nypix"].integer() - self._binsz = self["binsz"].real() - - # Read other parameters - self._outfile = self["outfile"].filename() - - # Get other parameeters - self._calc_ulimit = self["calc_ulim"].boolean() - self._calc_ts = self["calc_ts"].boolean() - self._fix_bkg = self["fix_bkg"].boolean() - self._fix_srcs = self["fix_srcs"].boolean() - - # Set some fixed parameters - self._chatter = self["chatter"].integer() - self._clobber = self["clobber"].boolean() - self._debug = self["debug"].boolean() + # Read ahead output parameters + if self._read_ahead(): + self["outfile"].filename() # Write input parameters into logger if self._logTerse(): @@ -241,26 +141,101 @@ class cslightcrv(ctools.cscript): # Return return + def _create_tbounds(self): + """ + Creates light curve time bins. - # Public methods - def run(self): + The method reads the following user parameters: + tbinalg: Time binning algorithm + tbinfile: Time binning file (FITS or ASCII) + tmin: Start time (MJD) + tmax: Stop time (MJD) + tbins: Number of time bins + + Returns: + Light curve bins in form of a GTI. """ - Run the script. + # Initialise Good Time Intervals + gti = gammalib.GGti() + + # Get algorithm to use for defining the time intervals + algorithm = self["tbinalg"].string() + + # Handle a FITS or a ASCII file for time bin definition + if algorithm == "FILE": + + # Get the filename + filename = self["tbinfile"].filename() + + # If the file a FITS file then load GTIs + if filename.is_fits(): + gti.load(filename) + + # ... otherwise load file as CSV ASCII file + csv = gammalib.GCsv(filename) + for i in range(csv.nrows()): + tmin = gammalib.GTime() + tmax = gammalib.GTime() + tmin.mjd(csv.real(i,0)) + tmax.mjd(csv.real(i,1)) + gti.append(tmin,tmax) + + # Handle linear time binning + elif algorithm == "LIN": + + # Get start and stop time and number of time bins + time_min = self["tmin"].real() + time_max = self["tmax"].real() + nbins = self["tbins"].integer() + + # Compute time step and setup time intervals + time_step = (time_max - time_min) / float(nbins) + for i in range(nbins): + tmin = gammalib.GTime() + tmax = gammalib.GTime() + tmin.mjd(time_min + i *time_step) + tmax.mjd(time_min + (i+1)*time_step) + gti.append(tmin,tmax) + + # Handle usage of observation GTIs + elif algorithm == "GTI": + + # Append the GTIs of all observations + for obs in self._obs: + for i in range(obs.events().gti().size()): + gti.append(obs.events().gti().tstart(i), + obs.events().gti().tstop(i)) + + # ... otherwise raise an exception (this should never occur) + else: + raise AttributeError('Paramter tbinalg="'+algorithm+'" unknown. ' + 'Must be one of "FILE", "LIN" or "GTI".') + + # Return Good Time Intervals + return gti + + def _get_free_par_names(self): """ - # Switch screen logging on in debug mode - if self._logDebug(): - self._log.cout(True) + Return list of free parameter names. + + Returns: + List of free parameter names. + """ + # Initialise list of free parameter names + names = [] - # Get parameters - self._get_parameters() + # Collect list of free parameter names + for par in self._obs.models()[self._srcname]: + if par.is_free(): + names.append(par.name()) - # Write observation into logger - if self._logTerse(): - self._log("\n") - self._log.header1("Observation") - self._log(str(self._obs)) - self._log("\n") + # Return names + return names + def _adjust_model_pars(self): + """ + Adjust model parameters dependent on user parameters. + """ # Write header if self._logTerse(): self._log("\n") @@ -269,85 +244,341 @@ class cslightcrv(ctools.cscript): # Adjust model parameters dependent on input user parameters for model in self._obs.models(): - # Set TS flag for all models to false. - # Source of interest will be set to true later + # Set TS flag for all models to false. The source of interest + # will be set to true later model.tscalc(False) # Log model name - if self._logExplicit(): + if self._logNormal(): self._log.header3(model.name()) # Deal with the source of interest if model.name() == self._srcname: - if self._calc_ts: + if self["calc_ts"].boolean(): model.tscalc(True) - elif self._fix_bkg and not model.classname() == "GModelSky": + elif (self["fix_bkg"].boolean() and + model.classname() != "GModelSky"): for par in model: - if par.is_free() and self._logExplicit(): - self._log(" Fixing \""+par.name()+"\"\n") - par.fix() - - elif self._fix_srcs and model.classname() == "GModelSky": + if par.is_free(): + par.fix() + if self._logNormal(): + self._log(gammalib.parformat(par.name())) + self._log("fixed\n") + + elif (self["fix_srcs"].boolean() and + model.classname() == "GModelSky"): for par in model: - if par.is_free() and self._logExplicit(): - self._log(" Fixing \""+par.name()+"\"\n") - par.fix() + if par.is_free(): + par.fix() + if self._logNormal(): + self._log(gammalib.parformat(par.name())) + self._log("fixed\n") - # Write header - if self._logTerse(): - self._log("\n") - self._log.header1("Generate lightcurve") + # Return + return - # Initialise FITS Table with extension "LIGHTCURVE" - table = gammalib.GFitsBinTable(self._tbins.size()) - table.extname("LIGHTCURVE") + def _create_fits_table(self, results): + """ + Creates FITS binary table containing light curve results. + + Args: + results: List of result dictionaries. - # Add Header for compatibility with gammalib.GMWLSpectrum - table.card("INSTRUME", "CTA", "Name of Instrument") - table.card("TELESCOP", "CTA", "Name of Telescope") + Returns: + FITS binary table containing light curve. + """ + # Determine number of rows in FITS table + nrows = len(results) - # Create FITS table columns - MJD = gammalib.GFitsTableDoubleCol("MJD", self._tbins.size()) - MJD.unit("days") - e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self._tbins.size()) - e_MJD.unit("days") + # Create FITS Table with extension "LIGHTCURVE" + table = gammalib.GFitsBinTable(nrows) + table.extname("LIGHTCURVE") - # Create a FITS column for every free parameter - columns = [] + # Append time columns + mjd = gammalib.GFitsTableDoubleCol("MJD", nrows) + e_mjd = gammalib.GFitsTableDoubleCol("e_MJD", nrows) + mjd.unit("days") + e_mjd.unit("days") + for i, result in enumerate(results): + mjd[i] = result["mjd"] + e_mjd[i] = result["e_mjd"] + table.append(mjd) + table.append(e_mjd) + + # Create parameter columns for par in self._obs.models()[self._srcname]: if par.is_free(): - col = gammalib.GFitsTableDoubleCol(par.name(), self._tbins.size()) + name = par.name() + e_name = "e_"+par.name() + col = gammalib.GFitsTableDoubleCol(name, nrows) + e_col = gammalib.GFitsTableDoubleCol(e_name, nrows) col.unit(par.unit()) - columns.append(col) - e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self._tbins.size()) e_col.unit(par.unit()) - columns.append(e_col) + for i, result in enumerate(results): + col[i] = result['values'][name] + e_col[i] = result['values'][e_name] + table.append(col) + table.append(e_col) + + # Append Test Statistic column + ts = gammalib.GFitsTableDoubleCol("TS", nrows) + for i, result in enumerate(results): + ts[i] = result["ts"] + table.append(ts) + + # Append upper limit column + ulimit = gammalib.GFitsTableDoubleCol("UpperLimit", nrows) + ulimit.unit("ph/cm2/s") + for i, result in enumerate(results): + ulimit[i] = result['ulimit'] + table.append(ulimit) + + # Return table + return table + + def _compute_ulimit(self, obs): + """ + Computes upper flux limit. - # Create TS and upper limit columns - TSvalues = gammalib.GFitsTableDoubleCol("TS", self._tbins.size()) - ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self._tbins.size()) - ulim_values.unit("ph/cm2/s") + Args: + obs: Observation container. - # Loop over energy bins - for i in range(self._tbins.size()): + Returns: + Upper flux limit (-1 of not computed). + """ + # Initialise upper flux limit + ulimit_value = -1.0 + + # Perform computation only if requested + if self["calc_ulim"].boolean(): - # Log information - if self._logTerse(): - self._log("\n") - self._log.header2("Time bin "+str(i)) + # Write header in logger + if self._logExplicit(): + self._log.header3("Computing upper flux limit") + + # Create upper limit object + ulimit = ctools.ctulimit(obs) + ulimit["srcname"] = self._srcname + ulimit["eref"] = 1.0 + + # Try to run upper limit and catch exceptions + try: + ulimit.run() + ulimit_value = ulimit.flux_ulimit() + except: + if self._logTerse(): + self._log("Upper limit flux calculation failed.s\n") + ulimit_value = -1.0 + + # Return upper limit + return ulimit_value + + def _bin_observation(self, obs): + """ + Bin an observation is a binned analysis was requested. + + Args: + obs: Observation container. + + Returns: + Observation container with a binned. + """ + # Header + if self._logExplicit(): + self._log.header3("Binning events") + + # Bin events + bin = ctools.ctbin(obs) + bin["usepnt"] = False + bin["ebinalg"] = "LOG" + bin["xref"] = self["xref"].real() + bin["yref"] = self["yref"].real() + bin["binsz"] = self["binsz"].real() + bin["nxpix"] = self["nxpix"].integer() + bin["nypix"] = self["nypix"].integer() + bin["enumbins"] = self["enumbins"].integer() + bin["emin"] = self["emin"].real() + bin["emax"] = self["emax"].real() + bin["coordsys"] = self["coordsys"].string() + bin["proj"] = self["proj"].string() + bin.run() + + # Header + if self._logExplicit(): + self._log.header3("Creating exposure cube") + + # Create exposure cube + expcube = ctools.ctexpcube(obs) + expcube["incube"] = "NONE" + expcube["usepnt"] = False + expcube["ebinalg"] = "LOG" + expcube["xref"] = self["xref"].real() + expcube["yref"] = self["yref"].real() + expcube["binsz"] = self["binsz"].real() + expcube["nxpix"] = self["nxpix"].integer() + expcube["nypix"] = self["nypix"].integer() + expcube["enumbins"] = self["enumbins"].integer() + expcube["emin"] = self["emin"].real() + expcube["emax"] = self["emax"].real() + expcube["coordsys"] = self["coordsys"].string() + expcube["proj"] = self["proj"].string() + expcube.run() + + # Header + if self._logExplicit(): + self._log.header3("Creating point spread function cube") + + # Create psf cube + psfcube = ctools.ctpsfcube(obs) + psfcube["incube"] = "NONE" + psfcube["usepnt"] = False + psfcube["ebinalg"] = "LOG" + psfcube["xref"] = self["xref"].real() + psfcube["yref"] = self["yref"].real() + psfcube["binsz"] = self["binsz"].real() + psfcube["nxpix"] = self["nxpix"].integer() + psfcube["nypix"] = self["nypix"].integer() + psfcube["enumbins"] = self["enumbins"].integer() + psfcube["emin"] = self["emin"].real() + psfcube["emax"] = self["emax"].real() + psfcube["coordsys"] = self["coordsys"].string() + psfcube["proj"] = self["proj"].string() + psfcube.run() + + # Check if we need to include energy dispersion + if self["edisp"].boolean(): + + # Header + if self._logExplicit(): + self._log.header3("Creating energy dispersion cube") + + # Create edisp cube + edispcube = ctools.ctedispcube(obs) + edispcube["incube"] = "NONE" + edispcube["usepnt"] = False + edispcube["ebinalg"] = "LOG" + edispcube["xref"] = self["xref"].real() + edispcube["yref"] = self["yref"].real() + edispcube["binsz"] = self["binsz"].real() + edispcube["nxpix"] = self["nxpix"].integer() + edispcube["nypix"] = self["nypix"].integer() + edispcube["enumbins"] = self["enumbins"].integer() + edispcube["emin"] = self["emin"].real() + edispcube["emax"] = self["emax"].real() + edispcube["coordsys"] = self["coordsys"].string() + edispcube["proj"] = self["proj"].string() + edispcube.run() + + # Header + if self._logExplicit(): + self._log.header3("Creating background cube") + + # Create background cube + bkgcube = ctools.ctbkgcube(obs) + bkgcube["incube"] = "NONE" + bkgcube["usepnt"] = False + bkgcube["ebinalg"] = "LOG" + bkgcube["xref"] = self["xref"].real() + bkgcube["yref"] = self["yref"].real() + bkgcube["binsz"] = self["binsz"].real() + bkgcube["nxpix"] = self["nxpix"].integer() + bkgcube["nypix"] = self["nypix"].integer() + bkgcube["enumbins"] = self["enumbins"].integer() + bkgcube["emin"] = self["emin"].real() + bkgcube["emax"] = self["emax"].real() + bkgcube["coordsys"] = self["coordsys"].string() + bkgcube["proj"] = self["proj"].string() + bkgcube.run() + + # Retrieve a new oberservation container + new_obs = bin.obs().copy() + + # Get new models + models = bkgcube.models() + + # Set stacked response + if self["edisp"].boolean(): + new_obs[0].response(expcube.expcube(), + psfcube.psfcube(), + edispcube.edispcube(), + bkgcube.bkgcube()) + else: + new_obs[0].response(expcube.expcube(), + psfcube.psfcube(), + bkgcube.bkgcube()) + + # Fix background models if required + if self["fix_bkg"].boolean(): + for model in models: + if model.classname() != "GModelSky": + for par in model: + par.fix() + + # Set models for new oberservation container + new_obs.models(models) + + # Return new oberservation container + return new_obs + + + # 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() + + # Write observation into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Observation") + self._log(str(self._obs)) + self._log("\n") + + # Adjust model parameters dependent on user parameters + self._adjust_model_pars() + + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Generate lightcurve") + + # Initialise list of result dictionaries + results = [] + + # Get source parameters + pars = self._get_free_par_names() + + # Loop over time bins + for i in range(self._tbins.size()): # Get time boundaries tmin = self._tbins.tstart(i) tmax = self._tbins.tstop(i) + # Write time bin into header + if self._logTerse(): + self._log.header2("MJD "+ + str(tmin.mjd())+"-"+ + str(tmax.mjd())) + # Compute time bin center and time width twidth = 0.5 * (tmax - tmin) # in seconds tmean = tmin + twidth - # Store time as MJD - MJD[i] = tmean.mjd() - e_MJD[i] = twidth / gammalib.sec_in_day # in days + # Initialise result dictionary + result = {'mjd': tmean.mjd(), + 'e_mjd': twidth / gammalib.sec_in_day, + 'ts': 0.0, + 'ulimit': 0.0, + 'pars': pars, + 'values': {}} # Log information if self._logExplicit(): @@ -355,8 +586,8 @@ class cslightcrv(ctools.cscript): # Select events select = ctools.ctselect(self._obs) - select["emin"] = self._emin - select["emax"] = self._emax + select["emin"] = self["emin"].real() + select["emax"] = self["emax"].real() select["tmin"] = tmin.convert(select._time_reference()) select["tmax"] = tmax.convert(select._time_reference()) select["rad"] = "UNDEFINED" @@ -367,233 +598,86 @@ class cslightcrv(ctools.cscript): # Retrieve observation obs = select.obs() - # Binned analysis - if self._binned: - - # Header - if self._logTerse(): - self._log.header3("Binning events") - - # Bin events - bin = ctools.ctbin(select.obs()) - bin["usepnt"] = False - bin["ebinalg"] = "LOG" - bin["xref"] = self._xref - bin["yref"] = self._yref - bin["binsz"] = self._binsz - bin["nxpix"] = self._nxpix - bin["nypix"] = self._nypix - bin["enumbins"] = self._ebins - bin["emin"] = self._emin - bin["emax"] = self._emax - bin["coordsys"] = self._coordsys - bin["proj"] = self._proj - bin.run() - - # Header - if self._logTerse(): - self._log.header3("Creating exposure cube") - - # Create exposure cube - expcube = ctools.ctexpcube(select.obs()) - expcube["incube"] = "NONE" - expcube["usepnt"] = False - expcube["ebinalg"] = "LOG" - expcube["xref"] = self._xref - expcube["yref"] = self._yref - expcube["binsz"] = self._binsz - expcube["nxpix"] = self._nxpix - expcube["nypix"] = self._nypix - expcube["enumbins"] = self._ebins - expcube["emin"] = self._emin - expcube["emax"] = self._emax - expcube["coordsys"] = self._coordsys - expcube["proj"] = self._proj - expcube.run() - - # Header - if self._logTerse(): - self._log.header3("Creating point spread function cube") - - # Create psf cube - psfcube = ctools.ctpsfcube(select.obs()) - psfcube["incube"] = "NONE" - psfcube["usepnt"] = False - psfcube["ebinalg"] = "LOG" - psfcube["xref"] = self._xref - psfcube["yref"] = self._yref - psfcube["binsz"] = self._binsz - psfcube["nxpix"] = self._nxpix - psfcube["nypix"] = self._nypix - psfcube["enumbins"] = self._ebins - psfcube["emin"] = self._emin - psfcube["emax"] = self._emax - psfcube["coordsys"] = self._coordsys - psfcube["proj"] = self._proj - psfcube.run() - - # Check if we need to include energy dispersion - if self._edisp: - - # Header - if self._logTerse(): - self._log.header3("Creating energy dispersion cube") - - # Create edisp cube - edispcube = ctools.ctedispcube(select.obs()) - edispcube["incube"] = "NONE" - edispcube["usepnt"] = False - edispcube["ebinalg"] = "LOG" - edispcube["xref"] = self._xref - edispcube["yref"] = self._yref - edispcube["binsz"] = self._binsz - edispcube["nxpix"] = self._nxpix - edispcube["nypix"] = self._nypix - edispcube["enumbins"] = self._ebins - edispcube["emin"] = self._emin - edispcube["emax"] = self._emax - edispcube["coordsys"] = self._coordsys - edispcube["proj"] = self._proj - edispcube.run() - - # Header - if self._logTerse(): - self._log.header3("Creating background cube") - - # Create background cube - bkgcube = ctools.ctbkgcube(select.obs()) - bkgcube["incube"] = "NONE" - bkgcube["usepnt"] = False - bkgcube["ebinalg"] = "LOG" - bkgcube["xref"] = self._xref - bkgcube["yref"] = self._yref - bkgcube["binsz"] = self._binsz - bkgcube["nxpix"] = self._nxpix - bkgcube["nypix"] = self._nypix - bkgcube["enumbins"] = self._ebins - bkgcube["emin"] = self._emin - bkgcube["emax"] = self._emax - bkgcube["coordsys"] = self._coordsys - bkgcube["proj"] = self._proj - bkgcube.run() - - # Set new binned observation - obs = bin.obs() - - # Get new models - models = bkgcube.models() - - # Set precomputed binned response - if self._edisp: - obs[0].response(expcube.expcube(), psfcube.psfcube(), - edispcube.edispcube(), bkgcube.bkgcube()) - else: - obs[0].response(expcube.expcube(), psfcube.psfcube(), - bkgcube.bkgcube()) - - # Fix background models if required - if self._fix_bkg: - for model in models: - if not model.classname() == "GModelSky": - for par in model: - par.fix() - - # Set new models to binned observation - obs.models(models) + # If a stacked analysis is requested then bin the events + # and compute the stacked response functions and setup + # an observation container with a single stacked observation. + if self._stacked: + obs = self._bin_observation(obs) # Header - if self._logTerse(): - self._log.header3("Performing fit") + if self._logExplicit(): + self._log.header3("Fitting the data") - # Likelihood + # Do maximum likelihood model fitting like = ctools.ctlike(obs) - like["edisp"] = self._edisp + like["edisp"] = self["edisp"].boolean() like.run() # Skip bin if no event was present if like.obs().logL() == 0.0: - # Log information + # Signal skipping of bin if self._logTerse(): - self._log("No event in this time bin. Bin is skipped\n") - - # Set all values to 0 - for col in columns: - col[i] = 0.0 - TSvalues[i] = 0.0 - ulim_values[i] = 0.0 - continue + self._log(gammalib.parformat("Warning")) + self._log("No event in this time bin, skip bin.\n") - # Get results - fitted_models = like.obs().models() - source = fitted_models[self._srcname] + # Set all results to 0 + for par in pars: + result['values'][par] = 0.0 + result['values']["e_"+par] = 0.0 - # Calculate Upper Limit - ulimit_value = -1.0 - if self._calc_ulimit: + # Append result + results.append(result) - # Logging information - if self._logTerse(): - self._log.header3("Computing upper limit") - - # Create upper limit object - ulimit = ctools.ctulimit(like.obs()) - ulimit["srcname"] = self._srcname - ulimit["eref"] = 1.0 - - # Try to run upper limit and catch exceptions - try: - ulimit.run() - ulimit_value = ulimit.flux_ulimit() - except: - if self._logTerse(): - self._log("Upper limit calculation failed\n") - ulimit_value = -1.0 - - # Get TS value - TS = -1.0 - if self._calc_ts: - TS = source.ts() - - # Set values for storage - TSvalues[i] = TS - - # Set FITS column values - for col in columns: - if "e_" == col.name()[:2]: - col[i] = source.spectral()[col.name()[2:]].error() - else: - col[i] = source.spectral()[col.name()].value() - - # Store upper limit value if available - if ulimit_value > 0.0: - ulim_values[i] = ulimit_value + # Continue with next time bin + continue - # Log information - if self._logExplicit(): - self._log.header3("Results of bin "+str(i)+": MJD "+str(tmin.mjd())+"-"+str(tmax.mjd())) - for col in columns: - if "e_" == col.name()[:2]: - continue - value = source.spectral()[col.name()].value() - error = source.spectral()[col.name()].error() - unit = source.spectral()[col.name()].unit() - self._log(" > "+col.name()+": "+str(value)+" +- "+str(error)+" "+unit+"\n") - if self._calc_ts and TSvalues[i] > 0.0: - self._log(" > TS = "+str(TS)+" \n") - if self._calc_ulimit and ulim_values[i] > 0.0: - self._log(" > UL = "+str(ulim_values[i])+" [ph/cm2/s]") - self._log("\n") + # Retrieve model fitting results for source of interest + source = like.obs().models()[self._srcname] - # Append filles columns to fits table - table.append(MJD) - table.append(e_MJD) - for col in columns: - table.append(col) - table.append(TSvalues) - table.append(ulim_values) + # Extract parameter values + for par in pars: + result['values'][par] = source.spectral()[par].value() + result['values']["e_"+par] = source.spectral()[par].error() - # Create the FITS file now + # Calculate upper limit (-1 if not computed) + ulimit_value = self._compute_ulimit(like.obs()) + if ulimit_value > 0.0: + result['ulimit'] = ulimit_value + + # Extract Test Statistic value + if self["calc_ts"].boolean(): + result['ts'] = source.ts() + + # Append result to list of dictionaries + results.append(result) + + # Log results for this time bin + if self._logNormal(): + self._log.header3("Results") + pars = self._get_free_par_names() + for par in pars: + value = source.spectral()[par].value() + error = source.spectral()[par].error() + unit = source.spectral()[par].unit() + self._log(gammalib.parformat(par)) + self._log(str(value)) + self._log(" +/- ") + self._log(str(error)) + self._log(" ") + self._log(unit) + self._log("\n") + if result['ulimit'] > 0.0: + self._log(gammalib.parformat("Upper flux limit")) + self._log(str(result['ulimit'])+" ph/cm2/s\n") + if self["calc_ts"].boolean(): + self._log(gammalib.parformat("Test Statistic")) + self._log(str(result['ts'])+"\n") + + # Create FITS table from results + table = self._create_fits_table(results) + + # Create FITS file and append FITS table to FITS file self._fits = gammalib.GFits() self._fits.append(table) @@ -604,6 +688,9 @@ class cslightcrv(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Read ahead output parameters self._read_ahead(True) @@ -618,35 +705,37 @@ class cslightcrv(ctools.cscript): def save(self): """ - Save lightcurve. + Save light curve. """ # Write header if self._logTerse(): - self._log("\n") - self._log.header1("Save lightcurve") + self._log('\n') + self._log.header1('Save light curve') - # Get outmap parameter - self._outfile = self["outfile"].filename() + # Get light curve filename + outfile = self["outfile"].filename() # Continue only filename and residual map are valid - if self._outfile != "NONE" and self._fits != None: + if self._fits != None: # Log file name if self._logTerse(): - self._log("Save lightcurve into \""+self._outfile+"\".\n") + self._log(gammalib.parformat("Light curve file")) + self._log(outfile.url()) + self._log("\n") # Save spectrum - self._fits.saveto(self._outfile, self._clobber) + self._fits.saveto(outfile, self._clobber()) # Return return def lightcurve(self): """ - Return lightcurve FITS file. + Return light curve FITS file. Returns: - FITS file containing lightcurve. + FITS file containing light curve. """ # Return return self._fits @@ -670,8 +759,5 @@ if __name__ == '__main__': # Create instance of application app = cslightcrv(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csmodelinfo.par b/cscripts/csmodelinfo.par index a2516e5..5506773 100644 --- a/cscripts/csmodelinfo.par +++ b/cscripts/csmodelinfo.par @@ -7,9 +7,9 @@ # File: csmodelinfo.par # # Version: 1.1.0 # # # -# Author: Juergen Knoedlseder # -# Jurgen.Knodlseder@irap.omp.eu # -# IRAP # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the csmodelinfo script # # # diff --git a/cscripts/csmodelinfo.py b/cscripts/csmodelinfo.py index 6f6fc1e..475d079 100755 --- a/cscripts/csmodelinfo.py +++ b/cscripts/csmodelinfo.py @@ -36,7 +36,7 @@ class csmodelinfo(ctools.cscript): with DS9. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -47,7 +47,7 @@ class csmodelinfo(ctools.cscript): # Initialise class members self._models = gammalib.GModels() - self._ds9file = "NONE" + self._ds9file = gammalib.GFilename("NONE") self._pnt_type = "" self._pnt_mark_size = 12 self._show_labels = True @@ -59,27 +59,11 @@ class csmodelinfo(ctools.cscript): self._show_ext_type = True self._free_color = "green" self._fixed_color = "magenta" - self._read_ahead = False - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - # Set logger properties - self._log_header() - self._log.date(True) - - # Return - return - - def __del__(self): - """ - Destructor. - """ # Return return @@ -106,7 +90,7 @@ class csmodelinfo(ctools.cscript): self._fixed_color = self["fixed_color"].string() # Read ahead DS9 filename - if self._read_ahead: + if self._read_ahead(): self._ds9file = self["ds9file"].filename() # Write input parameters into logger @@ -500,8 +484,11 @@ class csmodelinfo(ctools.cscript): """ Execute the script. """ - # Set read ahead flag - self._read_ahead = True + # Open logfile + self._logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) # Run the script self.run() @@ -521,8 +508,5 @@ if __name__ == '__main__': # Create instance of application app = csmodelinfo(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csmodelmerge.par b/cscripts/csmodelmerge.par index 8f36965..c52ea90 100644 --- a/cscripts/csmodelmerge.par +++ b/cscripts/csmodelmerge.par @@ -7,9 +7,9 @@ # File: csmodelmerge.par # # Version: 1.1.0 # # # -# Author: Juergen Knoedlseder # -# Jurgen.Knodlseder@irap.omp.eu # -# IRAP # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the csmodelmerge script # # # diff --git a/cscripts/csmodelmerge.py b/cscripts/csmodelmerge.py index cfc20a0..3d998d1 100755 --- a/cscripts/csmodelmerge.py +++ b/cscripts/csmodelmerge.py @@ -35,7 +35,7 @@ class csmodelmerge(ctools.cscript): a single model definition XML file. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -47,28 +47,11 @@ class csmodelmerge(ctools.cscript): # Initialise class members self._files = None self._models = gammalib.GModels() - self._outmodel = "NONE" - self._read_ahead = False - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) - # Return - return + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - def __del__(self): - """ - Destructor. - """ # Return return @@ -105,8 +88,8 @@ class csmodelmerge(ctools.cscript): raise RuntimeError(msg) # Read ahead output filename - if self._read_ahead: - self.outmodel = self["outmodel"].filename() + if self._read_ahead(): + self["outmodel"].filename() # Write input parameters into logger if self._logTerse(): @@ -178,16 +161,16 @@ class csmodelmerge(ctools.cscript): self._log.header1("Save models") # Get output filename in case it was not read ahead - self._outmodel = self["outmodel"].filename() + outmodel = self["outmodel"].filename() # Log filename if self._logTerse(): self._log(gammalib.parformat("Model definition XML file")) - self._log(self._outmodel.url()) + self._log(outmodel.url()) self._log("\n") # Save models - self._models.save(self._outmodel) + self._models.save(outmodel) # Return return @@ -196,8 +179,11 @@ class csmodelmerge(ctools.cscript): """ Execute the script. """ - # Set read ahead flag - self._read_ahead = True + # Open logfile + self._logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) # Run the script self.run() @@ -217,8 +203,5 @@ if __name__ == '__main__': # Create instance of application app = csmodelmerge(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csobs2caldb.par b/cscripts/csobs2caldb.par index 9f8c8a0..21d2316 100644 --- a/cscripts/csobs2caldb.par +++ b/cscripts/csobs2caldb.par @@ -7,9 +7,9 @@ # File: csobs2caldb.par # # Version: 1.1.0 # # # -# Author: Juergen Knoedlseder # -# Jurgen.Knodlseder@irap.omp.eu # -# IRAP # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the csobs2caldb script # # # diff --git a/cscripts/csobs2caldb.py b/cscripts/csobs2caldb.py index 1dd5076..664bac7 100755 --- a/cscripts/csobs2caldb.py +++ b/cscripts/csobs2caldb.py @@ -40,7 +40,7 @@ class csobs2caldb(ctools.cscript): index of any observation using the hidden "index" parameter. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -59,28 +59,13 @@ class csobs2caldb(ctools.cscript): self._base_dir = "" self._cal_dir = "" - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - # Private methods def _get_parameters(self): @@ -395,6 +380,9 @@ class csobs2caldb(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Run the script self.run() @@ -413,8 +401,5 @@ if __name__ == '__main__': # Create instance of application app = csobs2caldb(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csobsdef.py b/cscripts/csobsdef.py index 67f4e0c..ebe2592 100755 --- a/cscripts/csobsdef.py +++ b/cscripts/csobsdef.py @@ -72,7 +72,7 @@ class csobsdef(ctools.cscript): IRF in the "dummy" calibration database. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -85,31 +85,14 @@ class csobsdef(ctools.cscript): self._obs = gammalib.GObservations() self._inpnt = "" self._tmin = 0.0 - self._outobs = "NONE" self._chatter = 2 self._clobber = True self._debug = False - self._read_ahead = False - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - # Set logger properties - self._log_header() - self._log.date(True) - - # Return - return - - def __del__(self): - """ - Destructor. - """ # Return return @@ -123,8 +106,8 @@ class csobsdef(ctools.cscript): self._inpnt = self["inpnt"].filename() # Read ahead parameters - if self._read_ahead: - self._outobs = self["outobs"].filename() + if self._read_ahead(): + self["outobs"].filename() # Set some fixed parameters self._chatter = self["chatter"].integer() @@ -347,19 +330,19 @@ class csobsdef(ctools.cscript): self._log.header1("Save observation definition XML file") # Get output filename in case it was not read ahead - self._outobs = self["outobs"].filename() + outobs = self["outobs"].filename() # Check if observation definition XML file is valid - if self._outobs != "NONE": + if outobs.url() != "NONE": # Log filename if self._logTerse(): self._log(gammalib.parformat("Observation XML file")) - self._log(self._outobs.url()) + self._log(outobs.url()) self._log("\n") # Save observation definition XML file - self._obs.save(self._outobs) + self._obs.save(outobs) # Return return @@ -368,8 +351,11 @@ class csobsdef(ctools.cscript): """ Execute the script. """ - # Set read ahead flag - self._read_ahead = True + # Open logfile + self._logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) # Run the script self.run() @@ -389,8 +375,5 @@ if __name__ == '__main__': # Create instance of application app = csobsdef(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csobsinfo.par b/cscripts/csobsinfo.par index 23a3ab5..ac66a73 100644 --- a/cscripts/csobsinfo.par +++ b/cscripts/csobsinfo.par @@ -7,9 +7,9 @@ # File: csobsinfo.par # # Version: 1.1.0 # # # -# Author: Juergen Knoedlseder # -# Jurgen.Knodlseder@irap.omp.eu # -# IRAP # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the csobsinfo script # # # diff --git a/cscripts/csobsinfo.py b/cscripts/csobsinfo.py index fe6418d..10ef2cc 100755 --- a/cscripts/csobsinfo.py +++ b/cscripts/csobsinfo.py @@ -31,7 +31,7 @@ class csobsinfo(ctools.cscript): Shows the content of an observation container. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -43,7 +43,6 @@ class csobsinfo(ctools.cscript): # Initialise class members self._obj_dir = None self._compute_offset = False - self._ds9file = "" self._offsets = [] self._zeniths = [] self._azimuths = [] @@ -51,36 +50,17 @@ class csobsinfo(ctools.cscript): self._pnt_dec = [] self._ebounds = gammalib.GEbounds() self._gti = gammalib.GGti() - self._read_ahead = False - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - self._obs = argv[0] - argv = argv[1:] - else: - self._obs = gammalib.GObservations() - self._obs.clear() - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + # Initialise observation container from constructor arguments. + self._obs = self._set_input_obs(argv) - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - # Private methods def _get_parameters(self): @@ -103,8 +83,8 @@ class csobsinfo(ctools.cscript): self._obj_dir.radec_deg(ra,dec) # Read ahead DS9 filename - if self._read_ahead: - self._ds9file = self["ds9file"].filename() + if self._read_ahead(): + self["ds9file"].filename() # Write input parameters into logger if self._logTerse(): @@ -387,19 +367,19 @@ class csobsinfo(ctools.cscript): self._log.header1("Save pointings in DS9 file") # Get output filename in case it was not read ahead - self._ds9file = self["ds9file"].filename() + ds9file = self["ds9file"].filename() # Check if DS9 file is valid - if self._ds9file != "NONE": + if ds9file.url() != "NONE": # Log filename if self._logTerse(): self._log(gammalib.parformat("DS9 filename")) - self._log(self._ds9file.url()) + self._log(ds9file.url()) self._log("\n") # Open file - f = open(self._ds9file.url(),"w") + f = open(ds9file.url(),"w") # Write coordinate system f.write("fk5\n") @@ -455,8 +435,11 @@ class csobsinfo(ctools.cscript): """ Execute the script. """ - # Set read ahead flag - self._read_ahead = True + # Open logfile + self._logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) # Run the script self.run() @@ -476,8 +459,5 @@ if __name__ == '__main__': # Create instance of application app = csobsinfo(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/cspull.py b/cscripts/cspull.py index 1e05f3f..ec659ef 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -38,7 +38,7 @@ class cspull(ctools.cscript): script behaves just as a regular ctool. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -48,7 +48,6 @@ class cspull(ctools.cscript): self._version = "1.1.0" # Initialise some members - self._obs = None self._model = None self._inmodel = None self._edisp = False @@ -62,25 +61,13 @@ class cspull(ctools.cscript): self._proj = "TAN" self._log_clients = False - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise observation container from constructor arguments. + self._obs = self._set_input_obs(argv) - # Return - return + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - def __del__(self): - """ - Destructor. - """ # Return return @@ -310,10 +297,10 @@ class cspull(ctools.cscript): # Fit model if self._profile: - models = _obs.models() + models = self._obs.models() for i in range(models.size()): model_name = models[i].name() - like = obsutils.cterror(_obs, model_name, + like = obsutils.cterror(self._obs, model_name, log=self._log_clients, debug=self._debug, chatter=self._chatter) @@ -442,6 +429,9 @@ class cspull(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Run the script self.run() @@ -467,8 +457,5 @@ if __name__ == '__main__': # Create instance of application app = cspull(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/csresmap.par b/cscripts/csresmap.par index 533ea76..e529caf 100644 --- a/cscripts/csresmap.par +++ b/cscripts/csresmap.par @@ -7,9 +7,9 @@ # File: csresmap.par # # Version: 1.1.0 # # # -# Author: Juergen Knoedlseder # -# Jurgen.Knodlseder@irap.omp.eu # -# IRAP # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the csresmap script # # # diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index 9db3ea5..8559c87 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -36,7 +36,7 @@ class csresmap(ctools.cscript): behaves just as a regular ctool. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -49,7 +49,6 @@ class csresmap(ctools.cscript): self._algorithm = "SUB" self._resmap = None self._modcube = "NONE" - self._outmap = "NONE" self._use_maps = False self._skip_binning = False self._publish = False @@ -74,35 +73,16 @@ class csresmap(ctools.cscript): self._debug = False self._log_clients = False - # Initialise observation - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - self._obs = argv[0] - argv = argv[1:] - else: - self._obs = gammalib.GObservations() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + # Initialise observation container from constructor arguments. + self._obs = self._set_input_obs(argv) - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - # Private methods def _get_parameters(self): @@ -177,7 +157,7 @@ class csresmap(ctools.cscript): # Read ahead output parameters if (self._read_ahead()): - self._outmap = self["outmap"].filename() + self["outmap"].filename() # Write input parameters into logger if self._logTerse(): @@ -318,6 +298,9 @@ class csresmap(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Read ahead output parameters self._read_ahead(True) @@ -340,17 +323,19 @@ class csresmap(ctools.cscript): self._log.header1("Save residual map") # Get outmap parameter - self._outmap = self["outmap"].filename() + outmap = self["outmap"].filename() # Continue only filename and residual map are valid - if self._outmap != "NONE" and self._resmap != None: + if self._resmap != None: # Log file name if self._logTerse(): - self._log("Save residual map into \""+self._outmap+"\".\n") + self._log(gammalib.parformat("Residual map file")) + self._log(outmap.url()) + self._log("\n") # Save residual map - self._resmap.save(self._outmap, self._clobber) + self._resmap.save(outmap, self._clobber) # Return return @@ -405,8 +390,5 @@ if __name__ == '__main__': # Create instance of application app = csresmap(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 7d5eef0..e526b67 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -41,7 +41,7 @@ class cssens(ctools.cscript): the likelihood difference between fitting with and without the test source. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -54,7 +54,7 @@ class cssens(ctools.cscript): self._obs = gammalib.GObservations() self._ebounds = gammalib.GEbounds() self._srcname = "" - self._outfile = "" + self._outfile = gammalib.GFilename() self._ra = None self._dec = None self._edisp = False @@ -71,28 +71,13 @@ class cssens(ctools.cscript): self._log_clients = False self._debug = False - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - # Private methods def _get_parameters(self): @@ -734,6 +719,9 @@ class cssens(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Run the script self.run() @@ -750,8 +738,5 @@ if __name__ == '__main__': # Create instance of application app = cssens(sys.argv) - # Open logfile - app._logFileOpen() - # Run application app.execute() diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 327e029..f011e20 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -51,7 +51,6 @@ class csspec(ctools.cscript): self._binned_mode = False self._srcname = "" self._edisp = False - self._outfile = "NONE" self._calc_ulimit = True self._calc_ts = True self._fix_bkg = False @@ -60,24 +59,12 @@ class csspec(ctools.cscript): self._clobber = True self._debug = False - # Initialise observation - if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - self._obs = argv[0] - argv = argv[1:] - else: - self._obs = gammalib.GObservations() - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + # Initialise observation container from constructor arguments. + self._obs = self._set_input_obs(argv) - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return @@ -132,8 +119,8 @@ class csspec(ctools.cscript): self._debug = self["debug"].boolean() # Read ahead output parameters - if (self._read_ahead()): - self._outfile = self["outfile"].filename() + if self._read_ahead(): + self["outfile"].filename() # Write input parameters into logger if self._logTerse(): @@ -494,6 +481,9 @@ class csspec(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Read ahead output parameters self._read_ahead(True) @@ -516,17 +506,19 @@ class csspec(ctools.cscript): self._log.header1("Save spectrum") # Get outmap parameter - self._outfile = self["outfile"].filename() + outfile = self["outfile"].filename() # Continue only filename and residual map are valid - if self._outfile != "NONE" and self._fits != None: + if self._fits != None: # Log file name if self._logTerse(): - self._log("Save spectrum into \""+self._outfile+"\".\n") + self._log(gammalib.parformat("Spectrum file")) + self._log(outfile.url()) + self._log("\n") # Save spectrum - self._fits.saveto(self._outfile, self._clobber) + self._fits.saveto(outfile, self._clobber) # Return return @@ -560,8 +552,5 @@ if __name__ == '__main__': # Create instance of application app = csspec(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index 0af2a0a..137671b 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -39,7 +39,7 @@ class cstsdist(ctools.cscript): script behaves just as a regular ctool. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -55,34 +55,19 @@ class cstsdist(ctools.cscript): self._enumbins = 0 self._npix = 0 self._binsz = 0.0 - self._outfile = "ts.dat" + self._outfile = gammalib.GFilename("ts.dat") self._ntrials = 10 self._edisp = False self._debug = False self._log_clients = False - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) ==1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) # Return return - def __del__(self): - """ - Destructor. - """ - # Return - return - # Private methods def _get_parameters(self): @@ -463,6 +448,9 @@ class cstsdist(ctools.cscript): """ Execute the script. """ + # Open logfile + self._logFileOpen() + # Run the script self.run() @@ -478,8 +466,5 @@ if __name__ == '__main__': # Create instance of application app = cstsdist(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/cstsmapmerge.par b/cscripts/cstsmapmerge.par index f195ec7..ce531b8 100644 --- a/cscripts/cstsmapmerge.par +++ b/cscripts/cstsmapmerge.par @@ -7,9 +7,9 @@ # File: cstsmapmerge.par # # Version: 1.1.0 # # # -# Author: Juergen Knoedlseder # -# Jurgen.Knodlseder@irap.omp.eu # -# IRAP # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the cstsmapmerge script # # # diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py index af41ea8..4eccec7 100755 --- a/cscripts/cstsmapmerge.py +++ b/cscripts/cstsmapmerge.py @@ -33,7 +33,7 @@ class cstsmapmerge(ctools.cscript): Merge Test Statistic maps. """ - # Constructors and destructors + # Constructor def __init__(self, *argv): """ Constructor. @@ -52,28 +52,11 @@ class cstsmapmerge(ctools.cscript): self._merged_files = [] self._overwrite = True self._delete = False - self._outmap = "NONE" - self._read_ahead = False - - # Initialise application - if len(argv) == 0: - ctools.cscript.__init__(self, self._name, self._version) - elif len(argv) == 1: - ctools.cscript.__init__(self, self._name, self._version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self._log_header() - self._log.date(True) - # Return - return + # Initialise application by calling the appropriate class + # constructor. + self._init_cscript(argv) - def __del__(self): - """ - Destructor. - """ # Return return @@ -120,8 +103,8 @@ class cstsmapmerge(ctools.cscript): self._delete = self["delete"].boolean() # Read ahead output filename - if self._read_ahead: - self._outmap = self["outmap"].filename() + if self._read_ahead(): + self["outmap"].filename() # Write input parameters into logger if self._logTerse(): @@ -391,12 +374,12 @@ class cstsmapmerge(ctools.cscript): self._log.header1("Save TS map") # Get output filename in case it was not read ahead - self._outmap = self["outmap"].filename() + outmap = self["outmap"].filename() # Log filename if self._logTerse(): self._log(gammalib.parformat("TS map file")) - self._log(self._outmap.url()) + self._log(outmap.url()) self._log("\n") # Create FITS file @@ -426,7 +409,7 @@ class cstsmapmerge(ctools.cscript): fits[fits.size()-1].extname("STATUS MAP") # Save FITS file - fits.saveto(self._outmap, self._clobber()) + fits.saveto(outmap, self._clobber()) # Delete TS input maps if requested if self._delete: @@ -444,8 +427,11 @@ class cstsmapmerge(ctools.cscript): """ Execute the script. """ - # Set read ahead flag - self._read_ahead = True + # Open logfile + self._logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) # Run the script self.run() @@ -465,8 +451,5 @@ if __name__ == '__main__': # Create instance of application app = cstsmapmerge(sys.argv) - # Open logfile - app._logFileOpen() - # Execute application app.execute() diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index fc07d24..d2871b1 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -1,8 +1,7 @@ # ========================================================================== -# This script provides a number of functions that are useful for handling -# CTA observations. +# CTA observation handling support functions. # -# Copyright (C) 2011-2015 Juergen Knoedlseder +# Copyright (C) 2011-2016 Juergen Knoedlseder # # 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 @@ -31,21 +30,28 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, """ Simulate events for all observations in the container. - Parameters: - obs - Observation container - Keywords: - log - Create log file(s) - debug - Create console dump? - edisp - Apply energy dispersion? - seed - Seed value for simulations (default: 0) - nbins - Number of energy bins (default: 0=unbinned) - binsz - Pixel size for binned simulation (deg/pixel) - npix - Number of pixels in X and Y for binned simulation + Args: + obs: Observation container without events + + Kwargs: + log: Create log file(s) + debug: Create console dump? + chatter: Chatter level + edisp: Apply energy dispersion? + seed: Seed value for simulations (default: 0) + nbins: Number of energy bins (default: 0=unbinned) + binsz: Pixel size for binned simulation (deg/pixel) + npix: Number of pixels in X and Y for binned simulation + proj: Projection for binned simulation + coord: Coordinate system for binned simulation + + Returns: + Observation container, filled with simulated events. """ # Allocate ctobssim application and set parameters sim = ctools.ctobssim(obs) - sim["seed"] = seed + sim["seed"] = seed sim["edisp"] = edisp # Optionally open the log file @@ -85,16 +91,16 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0, # Allocate ctbin application and set parameters bin = ctools.ctbin(sim.obs()) - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax + bin["ebinalg"] = "LOG" + bin["emin"] = emin + bin["emax"] = emax bin["enumbins"] = nbins - bin["usepnt"] = True # Use pointing for map centre - bin["nxpix"] = npix - bin["nypix"] = npix - bin["binsz"] = binsz + bin["usepnt"] = True # Use pointing for map centre + bin["nxpix"] = npix + bin["nypix"] = npix + bin["binsz"] = binsz bin["coordsys"] = coord - bin["proj"] = proj + bin["proj"] = proj # Optionally open the log file if log: @@ -350,28 +356,30 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ # ======================= # def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ emin=0.1, emax=100.0, rad=5.0, \ - irf="South_50h", caldb="prod2", id="000000", instrument="CTA"): + irf="South_50h", caldb="prod2", id="000000"): """ - Returns a single CTA observation containing an empty CTA event list. - By looping over this function you can add CTA observations to the - observation container. - - Parameters: - pntdir - Pointing direction [GSkyDir] - Keywords: - tstart - Start time [seconds] (default: 0.0) - duration - Duration of observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - emin - Minimum event energy [TeV] (default: 0.1) - emax - Maximum event energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - irf - Instrument response function (default: cta_dummy_irf) - caldb - Calibration database path (default: "dummy") - id - Run identifier (default: "000000") - instrument - Intrument (default: "CTA") + Set a single CTA observation. + + The function sets a single CTA observation containing an empty CTA + event list. By looping over this function you can add CTA observations + to the observation container. + + Args: + pntdir: Pointing direction [GSkyDir] + + Kwargs: + tstart: Start time (seconds) (default: 0.0) + duration: Duration of observation (seconds) (default: 1800.0) + deadc: Deadtime correction factor (default: 0.95) + emin: Minimum event energy (TeV) (default: 0.1) + emax: Maximum event energy (TeV) (default: 100.0) + rad: ROI radius used for analysis (deg) (default: 5.0) + irf: Instrument response function (default: "South_50h") + caldb: Calibration database path (default: "prod2") + id: Run identifier (default: "000000") """ # Allocate CTA observation - obs_cta = gammalib.GCTAObservation(instrument) + obs_cta = gammalib.GCTAObservation() # Set calibration database db = gammalib.GCaldb() @@ -397,7 +405,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) # Set energy boundaries - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), gammalib.GEnergy(emax, "TeV")) # Allocate event list diff --git a/pyext/ctedispcube.i b/pyext/ctedispcube.i new file mode 100644 index 0000000..3fd8d20 --- /dev/null +++ b/pyext/ctedispcube.i @@ -0,0 +1,61 @@ +/*************************************************************************** + * ctedispcube - Energy dispersion cube generation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Maria Haupt * + * ----------------------------------------------------------------------- * + * * + * 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 . * + * * + ***************************************************************************/ +/** + * @file ctedispcube.i + * @brief Energy dispersion cube generation tool definition + * @author Chia-Chun Lu + */ + +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctedispcube.hpp" +%} + +/***********************************************************************//** + * @class ctedispcube + * + * @brief Energy dispersion cube generation tool + ***************************************************************************/ +class ctedispcube : public ctool { + +public: + // Constructors and destructors + ctedispcube(void); + explicit ctedispcube(const GObservations& obs); + ctedispcube(int argc, char *argv[]); + ctedispcube(const ctedispcube& app); + virtual ~ctedispcube(void); + + // Methods + void clear(void); + void run(void); + void save(void); + const GCTACubeEdisp& edispcube(void) const; +}; + +/***********************************************************************//** + * @brief Energy dispersion cube generation tool Python extensions + ***************************************************************************/ +%extend ctedispcube { + ctedispcube copy() { + return (*self); + } +} diff --git a/pyext/ctool.i b/pyext/ctool.i index 00eaff8..89234ec 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -99,6 +99,7 @@ public: %rename(_set_obs_bounds) set_obs_bounds(GObservations& obs); %rename(_get_mean_pointing) get_mean_pointing(const GObservations& obs); %rename(_get_current_rss) get_current_rss(); + %rename(_get_obs_header) get_obs_header(); // Protected methods const bool& read_ahead(void) const; @@ -122,6 +123,7 @@ public: void set_obs_bounds(GObservations& obs); GSkyDir get_mean_pointing(const GObservations& obs); size_t get_current_rss(void); + std::string get_obs_header(const GObservation* obs); }; @@ -152,7 +154,7 @@ public: /***********************************************************************//** - * @brief ctool base class Python extension + * @brief ctool base class Python extensions ***************************************************************************/ %extend ctool { void _read_ahead(const bool& flag) { @@ -160,3 +162,35 @@ public: return; } } + +/***********************************************************************//** + * @brief cscript base class Python extensions + ***************************************************************************/ +%pythoncode %{ + +# Initialise observation container from constructor arguments. In case that +# an observation container is provided as argument, this container will be +# used to initialise the class. +def _set_input_obs(self, argv): + if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): + obs = argv[0] + argv = argv[1:] + else: + obs = gammalib.GObservations() + return obs +cscript._set_input_obs = _set_input_obs + +# Initialise application by calling the appropriate class constructor. +def _init_cscript(self, argv): + if len(argv) == 0: + cscript.__init__(self, self._name, self._version) + elif len(argv) == 1: + cscript.__init__(self, self._name, self._version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + # Set logger properties + self._log_header() + self._log.date(True) +cscript._init_cscript = _init_cscript + +%} diff --git a/pyext/ctools/tools.i b/pyext/ctools/tools.i index 2cbeaf4..d5ffe3a 100755 --- a/pyext/ctools/tools.i +++ b/pyext/ctools/tools.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctools - SWIG file * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -54,6 +54,7 @@ %include "ctskymap.i" %include "ctexpcube.i" %include "ctpsfcube.i" +%include "ctedispcube.i" %include "ctbkgcube.i" %include "ctcubemask.i" %include "cttsmap.i" @@ -66,4 +67,4 @@ def test(): from ctools.tests import test_python_ctools test_python_ctools.test(installed=True) -%} \ No newline at end of file +%} diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index b169956..4da9f0e 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -52,6 +52,8 @@ const GEnergy g_energy_margin(1.0e-12, "TeV"); /***********************************************************************//** * @brief Void constructor + * + * Constructs an empty ctbin tool. ***************************************************************************/ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) { @@ -68,8 +70,7 @@ ctbin::ctbin(void) : ctool(CTBIN_NAME, CTBIN_VERSION) * * param[in] obs Observation container. * - * This method creates an instance of the class by copying an existing - * observations container. + * Constructs ctbin tool from an observation container. ***************************************************************************/ ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) { @@ -90,6 +91,9 @@ ctbin::ctbin(const GObservations& obs) : ctool(CTBIN_NAME, CTBIN_VERSION) * * @param[in] argc Number of arguments in command line. * @param[in] argv Array of command line arguments. + * + * Constructs ctbin tool using command line arguments for user parameter + * setting. ***************************************************************************/ ctbin::ctbin(int argc, char *argv[]) : ctool(CTBIN_NAME, CTBIN_VERSION, argc, argv) @@ -106,6 +110,8 @@ ctbin::ctbin(int argc, char *argv[]) : * @brief Copy constructor * * @param[in] app Application. + * + * Constructs ctbin tool from another ctbin instance. ***************************************************************************/ ctbin::ctbin(const ctbin& app) : ctool(app) { @@ -122,6 +128,8 @@ ctbin::ctbin(const ctbin& app) : ctool(app) /***********************************************************************//** * @brief Destructor + * + * Destructs ctbin tool. ***************************************************************************/ ctbin::~ctbin(void) { @@ -142,7 +150,10 @@ ctbin::~ctbin(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] app Application. + * @param[in] app ctbin tool. + * @return ctbin tool. + * + * Assigns ctbin tool. ***************************************************************************/ ctbin& ctbin::operator=(const ctbin& app) { @@ -175,7 +186,9 @@ ctbin& ctbin::operator=(const ctbin& app) ==========================================================================*/ /***********************************************************************//** - * @brief Clear instance + * @brief Clear ctbin tool + * + * Clears ctbin tool. ***************************************************************************/ void ctbin::clear(void) { @@ -195,12 +208,10 @@ void ctbin::clear(void) /***********************************************************************//** - * @brief Bin the event data + * @brief Run the ctbin tool * - * This method loops over all observations found in the observation conatiner - * and bins all events from the event list(s) into counts map(s). Note that - * each event list is binned in a separate counts map, hence no summing of - * events is done. + * Gets the user parameters and loops over all CTA observations in the + * observation container to bin the events into a single counts cube. ***************************************************************************/ void ctbin::run(void) { @@ -212,33 +223,17 @@ void ctbin::run(void) // Get task parameters get_parameters(); - // Write parameters into logger - if (logTerse()) { - log_parameters(); - log << std::endl; - } - // Write observation(s) into logger if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } + log.header1(gammalib::number("Observation", m_obs.size())); log << m_obs << std::endl; } // Write header if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Bin observations"); - } - else { - log.header1("Bin observation"); - } + log.header1(gammalib::number("Bin observation", m_obs.size())); } // Loop over all observations in the container @@ -246,14 +241,7 @@ void ctbin::run(void) // Write header for observation if (logTerse()) { - std::string header = m_obs[i]->instrument() + " observation"; - if (m_obs[i]->name().length() > 1) { - header += " \"" + m_obs[i]->name() + "\""; - } - if (m_obs[i]->id().length() > 1) { - header += " (id=" + m_obs[i]->id() +")"; - } - log.header3(header); + log.header3(get_obs_header(m_obs[i])); } // Get CTA observation @@ -293,12 +281,7 @@ void ctbin::run(void) // Write observation(s) into logger if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Binned observations"); - } - else { - log.header1("Binned observation"); - } + log.header1(gammalib::number("Binned observation", m_obs.size())); log << m_obs << std::endl; } @@ -339,8 +322,8 @@ void ctbin::save(void) // Log filename if (logTerse()) { - log << "Save counts cube into file \""+m_outcube+"\"."; - log << std::endl; + log << gammalib::parformat("Counts cube file"); + log << m_outcube.url() << std::endl; } // Save cube @@ -376,7 +359,8 @@ void ctbin::publish(const std::string& name) // Log filename if (logTerse()) { - log << "Publish \""+user_name+"\" counts cube." << std::endl; + log << gammalib::parformat("Counts cube name"); + log << user_name << std::endl; } // Publish counts cube @@ -496,6 +480,12 @@ void ctbin::get_parameters(void) m_outcube = (*this)["outcube"].filename(); } + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + // Return return; } @@ -521,8 +511,8 @@ void ctbin::fill_cube(GCTAObservation* obs) const GCTAEventList* events = dynamic_cast(obs->events()); if (events == NULL) { std::string msg = "CTA Observation does not contain an event " - "list. Event list information is needed to " - "fill the counts map."; + "list. An event list is needed to fill the " + "counts cube."; throw GException::invalid_value(G_FILL_CUBE, msg); } @@ -534,7 +524,7 @@ void ctbin::fill_cube(GCTAObservation* obs) // Check for RoI sanity if (!roi.is_valid()) { - std::string msg = "No RoI information found in input observation " + std::string msg = "No valid RoI found in input observation " "\""+obs->name()+"\". Run ctselect to specify " "an RoI for this observation before running " "ctbin."; @@ -624,10 +614,12 @@ void ctbin::fill_cube(GCTAObservation* obs) * @brief Create output observation container. * * Creates an output observation container that combines all input CTA - * observation into a single cube-style observation. All non CTA observations - * present in the observation container are kept. The method furthermore - * conserves any response information in case that a single CTA observation - * is provided. This supports the original binned analysis. + * observation into a single stacked observation. All non-CTA observations + * and all binned CTA observations that were present in the observation + * container are append to the observation container so that they can + * be used by other tools. The method furthermore conserves any response + * information in case that a single CTA observation is provided to support + * binned analysis. ***************************************************************************/ void ctbin::obs_cube(void) { @@ -641,26 +633,37 @@ void ctbin::obs_cube(void) GCTAObservation* obs = dynamic_cast(m_obs[0]); if (obs != NULL) { - // Only change the event type if we had an unbinned observation + // Change the event type if we had an unbinned observation if (obs->eventtype() == "EventList") { + // Assign cube to the observation obs->events(this->cube()); - obs->eventfile(""); - } // endif: observation was unbinned - else { // input observation was binned and skipped - // Create new and empty cube - GCTAEventCube cube = GCTAEventCube(m_cube, m_ebounds, obs->gti()); + } + + // ... otherwise the input observation was binned and hence + // skipped. In that case we simply append an empty counts + // cube + else { - // Assign empty cube to have a new binned observation + // Create empty counts cube + GCTAEventCube cube(m_cube, m_ebounds, obs->gti()); + + // Assign empty cube to observation obs->events(cube); + } + // Reset file name + obs->eventfile(""); + } // endif: obervation was valid } // endif: we only had one observation in the container - // ... otherwise put a single CTA observation in container + // ... otherwise put a single CTA observation in container and + // append all observations that have not been used in the binning + // to the container else { // Allocate observation container @@ -672,10 +675,25 @@ void ctbin::obs_cube(void) // Attach event cube to CTA observation obs.events(this->cube()); - // Set map centre as pointing - GSkyPixel pixel(0.5*double(m_cube.nx()), 0.5*double(m_cube.ny())); - GSkyDir centre = m_cube.pix2dir(pixel); - GCTAPointing pointing(centre); + // Compute average pointing direction for all CTA event lists + double ra = 0.0; + double dec = 0.0; + double number = 0.0; + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* cta = dynamic_cast(m_obs[i]); + if ((cta != NULL) && (cta->eventtype() == "EventList")) { + ra += cta->pointing().dir().ra(); + dec += cta->pointing().dir().dec(); + number += 1.0; + } + } + if (number > 0.0) { + ra /= number; + dec /= number; + } + GSkyDir dir; + dir.radec(ra, dec); + GCTAPointing pointing(dir); // Compute deadtime correction double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0; @@ -683,8 +701,8 @@ void ctbin::obs_cube(void) // Set CTA observation attributes obs.pointing(pointing); obs.obs_id(0); - obs.ra_obj(centre.ra_deg()); //!< Dummy - obs.dec_obj(centre.dec_deg()); //!< Dummy + obs.ra_obj(dir.ra_deg()); //!< Dummy + obs.dec_obj(dir.dec_deg()); //!< Dummy obs.ontime(m_ontime); obs.livetime(m_livetime); obs.deadc(deadc); @@ -701,6 +719,9 @@ void ctbin::obs_cube(void) if (obs == NULL) { container.append(*m_obs[i]); } + else if (obs->eventtype() != "EventList") { + container.append(*m_obs[i]); + } } // Set observation container diff --git a/src/ctbin/ctbin.par b/src/ctbin/ctbin.par index 8ce64dc..a55b087 100644 --- a/src/ctbin/ctbin.par +++ b/src/ctbin/ctbin.par @@ -7,27 +7,27 @@ # File: ctbin.par # # Version: 1.1.0 # # # -# Author: Jurgen Knodlseder # +# Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # # IRAP # # # -# Purpose: Parameter file for the ctbin application # +# Purpose: Parameter file for the ctbin tool # # # ############################################################################# # # General parameters #=================== -inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" -outcube, f, a, "cntcube.fits",,, "Output counts cube file" +inobs, f, a, events.fits,,, "Input event list or observation definition XML file" +outcube, f, a, cntcube.fits,,, "Output counts cube file" # # Energy binning parameters #========================== -ebinalg, s, a, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins" +ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins" emin, r, a, 0.1,,, "Start value for first energy bin in TeV" emax, r, a, 100.0,,, "Stop value for last energy bin in TeV" -enumbins, i, a, 20,,, "Number of energy bins" -ebinfile, f, a, "NONE",,, "Name of the file containing the energy bin definition" +enumbins, i, a, 20,1,200, "Number of energy bins" +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" # # Spatial binning parameters @@ -36,10 +36,10 @@ usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" -coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" -proj, s, a, "CAR",AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" # # Standard parameters @@ -48,5 +48,5 @@ publish, b, h, no,,, "Publish counts cube on VO Hub?" chatter, i, h, 2,0,4, "Chattiness of output" clobber, b, h, yes,,, "Overwrite existing output files with new output files?" debug, b, h, no,,, "Debugging mode activated" -mode, s, h, "ql",,, "Mode of automatic parameters" -logfile, f, h, "ctbin.log",,, "Log filename" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, ctbin.log,,, "Log filename" diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 68b664e..4486bea 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -60,6 +60,8 @@ const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees) /***********************************************************************//** * @brief Void constructor + * + * Constructs an empty ctobssim tool. ***************************************************************************/ ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION) { @@ -76,7 +78,7 @@ ctobssim::ctobssim(void) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION) * * @param[in] obs Observation container. * - * Constructs application by initialising from an observation container. + * Constructs ctobssim tool from an observation container. ***************************************************************************/ ctobssim::ctobssim(const GObservations& obs) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION) @@ -98,7 +100,8 @@ ctobssim::ctobssim(const GObservations& obs) : * @param[in] argc Number of arguments in command line. * @param[in] argv Array of command line arguments. * - * Constructs application using command line arguments for parameter setting. + * Constructs ctobssim tool using command line arguments for user parameter + * setting. ***************************************************************************/ ctobssim::ctobssim(int argc, char *argv[]) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv) @@ -115,6 +118,8 @@ ctobssim::ctobssim(int argc, char *argv[]) : * @brief Copy constructor * * @param[in] app Application. + * + * Constructs ctobssim tool from another ctobssim instance. ***************************************************************************/ ctobssim::ctobssim(const ctobssim& app) : ctool(app) { @@ -131,6 +136,8 @@ ctobssim::ctobssim(const ctobssim& app) : ctool(app) /***********************************************************************//** * @brief Destructor + * + * Destructs ctobssim tool. ***************************************************************************/ ctobssim::~ctobssim(void) { @@ -151,8 +158,10 @@ ctobssim::~ctobssim(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] app Application. - * @return Application. + * @param[in] app ctobssim tool. + * @return ctobssim tool. + * + * Assigns ctobssim tool. ***************************************************************************/ ctobssim& ctobssim::operator=(const ctobssim& app) { @@ -185,7 +194,9 @@ ctobssim& ctobssim::operator=(const ctobssim& app) ==========================================================================*/ /***********************************************************************//** - * @brief Clear application + * @brief Clear ctobssim tool + * + * Clears ctobssim tool. ***************************************************************************/ void ctobssim::clear(void) { @@ -205,10 +216,10 @@ void ctobssim::clear(void) /***********************************************************************//** - * @brief Simulate event data + * @brief Run the ctobssim tool. * - * This method runs the simulation. Results are not saved by this method. - * Invoke "save" to save the results. + * Gets the user parameters, loops over all CTA observations in the + * observation container, and simulate events for each observation. ***************************************************************************/ void ctobssim::run(void) { @@ -286,24 +297,14 @@ void ctobssim::run(void) // Write observation(s) into logger if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } + log.header1(gammalib::number("Observation", m_obs.size())); log << m_obs << std::endl; } // Write header if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Simulate observations"); - } - else { - log.header1("Simulate observation"); - } + log.header1(gammalib::number("Simulate observation", m_obs.size())); } // From here on the code can be parallelized if OpenMP support @@ -334,14 +335,7 @@ void ctobssim::run(void) // Write header for observation if (logTerse()) { - std::string header = m_obs[i]->instrument() + " observation"; - if (m_obs[i]->name().length() > 1) { - header += " \"" + m_obs[i]->name() + "\""; - } - if (m_obs[i]->id().length() > 1) { - header += " (id=" + m_obs[i]->id() +")"; - } - wrklog.header3(header); + wrklog.header3(get_obs_header(m_obs[i])); } // Get pointer on CTA observation @@ -1455,8 +1449,8 @@ void ctobssim::save_fits(void) // Log filename if (logTerse()) { - log << "Save event list into file \""+m_outevents+"\"."; - log << std::endl; + log << gammalib::parformat("Event list file"); + log << m_outevents << std::endl; } // Save observation into FITS file @@ -1483,7 +1477,7 @@ void ctobssim::save_xml(void) { // Get output filename and prefix m_outevents = (*this)["outevents"].filename(); - m_prefix = (*this)["prefix"].string(); + m_prefix = (*this)["prefix"].string(); // Issue warning if output filename has no .xml suffix std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4)); @@ -1513,15 +1507,16 @@ void ctobssim::save_xml(void) if (obs->events()->size() != 0) { // Set event output file name - std::string outfile = m_prefix + gammalib::str(i) + ".fits"; + std::string outfile = m_prefix + gammalib::str(i) + + ".fits"; // Store output file name in observation obs->eventfile(outfile); // Log filename if (logTerse()) { - log << "Save event list into file \""+outfile+"\"."; - log << std::endl; + log << gammalib::parformat("Event list file"); + log << outfile << std::endl; } // Save observation into FITS file diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par index e2a8101..f6bfb2e 100644 --- a/src/ctobssim/ctobssim.par +++ b/src/ctobssim/ctobssim.par @@ -7,7 +7,7 @@ # File: ctobssim.par # # Version: 1.1.0 # # # -# Author: Jurgen Knodlseder # +# Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # # IRAP # # # @@ -17,13 +17,13 @@ # # General parameters #=================== -inobs, f, h, "NONE",,, "Input event list or observation definition XML file" -inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file" -caldb, s, a, "prod2",,, "Calibration database" -irf, s, a, "South_0.5h",,, "Instrument response function" +inobs, f, h, NONE,,, "Input event list or observation definition XML file" +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file" +caldb, s, a, prod2,,, "Calibration database" +irf, s, a, South_0.5h,,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" -outevents, f, a, "events.fits",,, "Output event data file or observation definition XML file" -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition XML file" +outevents, f, a, events.fits,,, "Output event data file or observation definition XML file" +prefix, s, h, sim_events_,,, "Prefix for event lists in observation definition XML file" # # Simulation parameters @@ -38,7 +38,7 @@ emin, r, a, 0.1,,, "Lower energy limit (TeV)" emax, r, a, 100.0,,, "Upper energy limit (TeV)" deadc, r, h, 0.95,0,1, "Average deadtime correction factor" maxrate, r, h, 1.0e6,,, "Maximum photon rate" -eslices, i, h, 10,1,100,"Number of energy slices" +eslices, i, h, 10,1,100, "Number of energy slices" # # Standard parameters @@ -46,5 +46,5 @@ eslices, i, h, 10,1,100,"Number of energy slices" chatter, i, h, 2,0,4, "Chattiness of output" clobber, b, h, yes,,, "Overwrite existing output files with new output files?" debug, b, h, no,,, "Debugging mode activated" -mode, s, h, "ql",,, "Mode of automatic parameters" -logfile, f, h, "ctobssim.log",,, "Log filename" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, ctobssim.log,,, "Log filename" diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index 5cc444d..dd9608a 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1108,6 +1108,39 @@ size_t ctool::get_current_rss(void) /***********************************************************************//** + * @brief Return observation header string + * + * @param[in] obs Pointer to observation. + * @return String containing observation information. + * + * Returns a string that contains observation information, including the + * instrument name, the observation name, and the observation ID. The format + * of the string is + * + * XXX observation "name" (id=YYYYY) + * + ***************************************************************************/ +std::string ctool::get_obs_header(const GObservation* obs) +{ + // Initialise header string + std::string header = obs->instrument() + " observation"; + + // If observation name is not empty then add name + if (!obs->name().empty()) { + header += " \"" + obs->name() + "\""; + } + + // If observation ID is not empty then add ID + if (!obs->id().empty()) { + header += " (id=" + obs->id() +")"; + } + + // Return header + return header; +} + + +/***********************************************************************//** * @brief Dumps help text in the console ***************************************************************************/ void ctool::provide_help(void) const diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 5355efa..e26a008 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2015 by Juergen Knoedlseder * + * copyright (C) 2014-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -90,17 +90,18 @@ protected: void set_obs_bounds(GObservations& obs); GSkyDir get_mean_pointing(const GObservations& obs); size_t get_current_rss(void); + std::string get_obs_header(const GObservation* obs); // Protected members - bool m_read_ahead; //!< Read ahead parameters + bool m_read_ahead; //!< Read ahead output parameters protected: // Protected methods - void provide_help(void) const; + void provide_help(void) const; // Protected members - bool m_use_xml; //!< Use XML file instead of FITS file for observations - GTimeReference m_cta_ref; //!< CTA time reference + bool m_use_xml; //!< Use XML file instead of FITS file for observations + GTimeReference m_cta_ref; //!< CTA time reference }; diff --git a/test/test_ctbin.py b/test/test_ctbin.py index 5368618..8884f0a 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the ctbin tool. +# This scripts performs unit tests for the ctbin module. # # Copyright (C) 2014-2016 Juergen Knoedlseder # @@ -22,13 +22,14 @@ import gammalib import ctools -# ========================= # -# Test class for ctbin tool # -# ========================= # +# =========================== # +# Test class for ctbin module # +# =========================== # class Test(gammalib.GPythonTestSuite): """ - Test class for ctbin tool. + Test class for ctbin module. """ + # Constructor def __init__(self): """ @@ -38,7 +39,7 @@ class Test(gammalib.GPythonTestSuite): gammalib.GPythonTestSuite.__init__(self) # Set members - self.events_name = "data/crab_events.fits" + self._events_name = "data/crab_events.fits" # Return return @@ -52,19 +53,19 @@ class Test(gammalib.GPythonTestSuite): self.name("ctbin") # Append tests - self.append(self.test_functional, "Test ctbin functionality") + self.append(self._test_ctbin, "Test ctbin") # Return return - # Test ctbin functionnality - def test_functional(self): + # Test ctbin + def _test_ctbin(self): """ - Test ctbin functionnality. + Test ctbin. """ # Set-up ctbin bin = ctools.ctbin() - bin["inobs"] = self.events_name + bin["inobs"] = self._events_name bin["outcube"] = "cntmap.fits" bin["ebinalg"] = "LOG" bin["emin"] = 0.1 @@ -78,27 +79,154 @@ class Test(gammalib.GPythonTestSuite): bin["xref"] = 83.63 bin["yref"] = 22.01 - # Run tool - self.test_try("Run ctbin") - try: - bin.run() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in ctbin: %s." % (e,) - except: - msg = "Exception occured in ctbin." - self.test_try_failure(msg) + # Run ctbin tool + bin.run() + + # Check content of observation and cube + self._test_observation(bin, 5542) + self._test_cube(bin.cube(), 5542) + + # Test copy constructor + cpy_bin = bin.copy() + + # Check content of observation and cube + self._test_observation(cpy_bin, 5542) + self._test_cube(cpy_bin.cube(), 5542) + + # Run copy of ctbin tool again + cpy_bin.run() + + # Check content of observation and cube. We expect now an empty + # event cube as on input the observation is binned, and any binned + # observation will be skipped, hence the counts cube should be + # empty. + self._test_observation(cpy_bin, 0) + self._test_cube(cpy_bin.cube(), 0) + + # Clear and run copy of ctbin again + #cpy_bin.clear() + #cpy_bin.run() +#RuntimeError: *** ERROR in GXmlAttribute::value(std::string): Invalid XML attribute value: *** ERROR in GPythonTestSuite::test: +#*** ERROR in GApplicationPars::operator[](std::string&): Invalid argument. Parameter "debug" has not been found in parameter file. +#Please specify a valid parameter name. # Save counts cube - self.test_try("Save counts cube") + bin.save() + + # Load counts cube and check content. + evt = gammalib.GCTAEventCube("cntmap.fits") + self._test_cube(evt, 5542) + + # Prepare observation container for stacked analysis + cta = gammalib.GCTAObservation(self._events_name) + obs = gammalib.GObservations() + cta.id("0001") + obs.append(cta) + cta.id("0002") + obs.append(cta) + cta.id("0003") + obs.append(cta) + + # Set-up ctbin using an observation container + bin = ctools.ctbin(obs) + bin["outcube"] = "cntmap.fits" + bin["ebinalg"] = "LOG" + bin["emin"] = 0.1 + bin["emax"] = 100.0 + bin["enumbins"] = 20 + bin["nxpix"] = 200 + bin["nypix"] = 200 + bin["binsz"] = 0.02 + bin["coordsys"] = "CEL" + bin["proj"] = "CAR" + bin["xref"] = 83.63 + bin["yref"] = 22.01 + + # Run ctbin tool + bin.run() + + # Check content of observation and cube (need multiplier=3 since + # three identical observations have been appended) + self._test_observation(bin, 5542, multiplier=3) + self._test_cube(bin.cube(), 5542, multiplier=3) + + # Set-up ctbin with invalid event file + bin = ctools.ctbin() + bin["inobs"] = "event_file_that_does_not_exist.fits" + bin["outcube"] = "cntmap.fits" + bin["ebinalg"] = "LOG" + bin["emin"] = 0.1 + bin["emax"] = 100.0 + bin["enumbins"] = 20 + bin["nxpix"] = 200 + bin["nypix"] = 200 + bin["binsz"] = 0.02 + bin["coordsys"] = "CEL" + bin["proj"] = "CAR" + bin["xref"] = 83.63 + bin["yref"] = 22.01 + + # Run ctbin tool + self.test_try("Run ctbin with invalid event file") try: - bin.save() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in saving counts cube: %s." % (e,) + bin.run() + self.test_try_failure() except: - msg = "Exception occured in saving counts cube." - self.test_try_failure(msg) + self.test_try_success() + + # Return + return + + # Check observation + def _test_observation(self, ctbin, nevents, multiplier=1): + """ + Test content of an observation. + + Args: + ctbin: ctbin instance + nevents: Expected number of events + + Kwargs: + multiplier: Observation multiplier + """ + # Test observation container + obs = gammalib.GCTAObservation(ctbin.obs()[0]) + pnt = obs.pointing() + self.test_value(ctbin.obs().size(), 1, + "There is one observation") + self.test_assert(obs.instrument() == "CTA", + "Observation is CTA observation") + self.test_value(obs.ontime(), 1800.0*multiplier, 1.0e-6, + "Ontime is 1800 sec") + self.test_value(obs.livetime(), 1710.0*multiplier, 1.0e-6, + "Livetime is 1710 sec") + self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, + "Pointing Right Ascension is 83.63 deg") + self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, + "Pointing Declination is 22.01 deg") + + # Test event cube + self._test_cube(obs.events(), nevents, multiplier=multiplier) + + # Return + return + + # Check event cube + def _test_cube(self, cube, nevents, multiplier=1): + """ + Test content of event cube." + + Args: + cube: Event cube + nevents: Expected number of events. + + Kwargs: + multiplier: Event number multiplier + """ + # Test event cube + self.test_value(cube.size(), 800000, "800000 event bins") + self.test_value(cube.number(), nevents*multiplier, + str(nevents*multiplier)+" events") # Return return diff --git a/test/test_ctobssim.py b/test/test_ctobssim.py index 1b36ae2..5816c06 100755 --- a/test/test_ctobssim.py +++ b/test/test_ctobssim.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the ctobssim tool. +# This scripts performs unit tests for the ctobssim module. # # Copyright (C) 2014-2016 Juergen Knoedlseder # @@ -20,15 +20,17 @@ # ========================================================================== import gammalib import ctools +from cscripts import obsutils -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ============================== # +# Test class for ctobssim module # +# ============================== # class Test(gammalib.GPythonTestSuite): """ - Test class for ctobssim tool. + Test class for ctobssim module. """ + # Constructor def __init__(self): """ @@ -38,9 +40,9 @@ class Test(gammalib.GPythonTestSuite): gammalib.GPythonTestSuite.__init__(self) # Set members - self.model_name = "data/crab.xml" - self.caldb = "irf" - self.irf = "cta_dummy_irf" + self._model_name = "data/crab.xml" + self._caldb = "irf" + self._irf = "cta_dummy_irf" # Return return @@ -54,23 +56,23 @@ class Test(gammalib.GPythonTestSuite): self.name("ctobssim") # Append tests - self.append(self.test_functional, "Test ctobssim functionality") - self.append(self.test_container, "Test ctobssim on observation container") + self.append(self._test_ctobssim, "Test ctobssim") + #self.append(self._test_container, "Test ctobssim on observation container") # Return return - # Test ctobssim functionnality - def test_functional(self): + # Test ctobssim + def _test_ctobssim(self): """ - Test ctobssim functionnality. + Test ctobssim. """ # Set-up ctobssim sim = ctools.ctobssim() - sim["inmodel"] = self.model_name + sim["inmodel"] = self._model_name sim["outevents"] = "events.fits" - sim["caldb"] = self.caldb - sim["irf"] = self.irf + sim["caldb"] = self._caldb + sim["irf"] = self._irf sim["ra"] = 83.63 sim["dec"] = 22.01 sim["rad"] = 10.0 @@ -80,154 +82,123 @@ class Test(gammalib.GPythonTestSuite): sim["emax"] = 100.0 # Run tool - self.test_try("Run ctobssim") - try: - sim.run() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in ctobssim: %s." % (e,) - except: - msg = "Exception occured in ctobssim." - self.test_try_failure(msg) + sim.run() - # Retrieve observation and check content - obs = gammalib.GCTAObservation(sim.obs()[0]) - pnt = obs.pointing() - evt = obs.events() - self.test_assert(obs.instrument() == "CTA", "Observation not a CTA observation") - self.test_value(sim.obs().size(), 1, "There is not a single observation") - self.test_value(obs.ontime(), 1800.0, 1.0e-6, "Ontime is not 1800 sec") - self.test_value(obs.livetime(), 1710.0, 1.0e-6, "Livetime is not 1710 sec") - self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "ROI Right Ascension is not 83.63 deg") - self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "ROI Declination is not 22.01 deg") - self.test_value(evt.size(), 4119, "Number of events is not 4119") + # Check content of observation + self._test_observation(sim) + self._test_list(sim.obs()[0].events(), 4119) # Save events - self.test_try("Save events") - try: - sim.save() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in saving events: %s." % (e,) - except: - msg = "Exception occured in saving events." - self.test_try_failure(msg) + sim.save() + + # Load counts cube and check content. + evt = gammalib.GCTAEventList("events.fits") + self._test_list(evt, 4119) - # Test ctobssim on observation container - def test_container(self): - """ - Test ctobssim on observation container. - """ # Set-up observation container - obs = self.set_obs(4) + pnts = [{'ra': 83.63, 'dec': 21.01}, + {'ra': 84.63, 'dec': 22.01}, + {'ra': 83.63, 'dec': 23.01}, + {'ra': 82.63, 'dec': 22.01}] + obs = obsutils.set_obs_list(pnts, caldb=self._caldb, irf=self._irf) - # Set-up ctobssim + # Set-up ctobssim from observation container sim = ctools.ctobssim(obs) sim["outevents"] = "sim_events.xml" - sim["inmodel"] = self.model_name + sim["inmodel"] = self._model_name + # Run tool - self.test_try("Run ctobssim") - try: - sim.run() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in ctobssim: %s." % (e,) - except: - msg = "Exception occured in ctobssim." - self.test_try_failure(msg) + sim.run() # Retrieve observation and check content - self.test_value(sim.obs().size(), 4, "There are not 4 observations") + self._test_observation(sim, nobs=4, pnts=pnts) + self._test_list(sim.obs()[0].events(), 4063) + self._test_list(sim.obs()[1].events(), 4059) + self._test_list(sim.obs()[2].events(), 4010) + self._test_list(sim.obs()[3].events(), 4138) # Save events - self.test_try("Save events") + sim.save() + + # Load events + obs = gammalib.GObservations("sim_events.xml") + + # Retrieve observation and check content + self._test_list(obs[0].events(), 4063) + self._test_list(obs[1].events(), 4059) + self._test_list(obs[2].events(), 4010) + self._test_list(obs[3].events(), 4138) + + # Set-up ctobssim with invalid event file + sim = ctools.ctobssim() + sim["inmodel"] = "model_file_that_does_not_exist.xml" + sim["outevents"] = "events.fits" + sim["caldb"] = self._caldb + sim["irf"] = self._irf + sim["ra"] = 83.63 + sim["dec"] = 22.01 + sim["rad"] = 10.0 + sim["tmin"] = 0.0 + sim["tmax"] = 1800.0 + sim["emin"] = 0.1 + sim["emax"] = 100.0 + + # Run ctbin tool + self.test_try("Run ctobssim with invalid model file") try: - sim.save() - self.test_try_success() - #except Exception as e: - # msg = "Exception occured in saving events: %s." % (e,) + sim.run() + self.test_try_failure() except: - msg = "Exception occured in saving events." - self.test_try_failure(msg) + self.test_try_success() + + # Return + return - # Setup observation container - def set_obs(self, number): + # Check observation + def _test_observation(self, ctobssim, nobs=1, + pnts=[{'ra': 83.63, 'dec': 22.01}]): """ - Setup observation container. + Test content of an observation. + + Args: + ctobssim: ctobssim instance + + Kwargs: + nobs: Number of observations + pnts: List of pointing dictionaries """ - # Initialise empty observation container - obs = gammalib.GObservations() - - # Initialise first time and identifier - tstart = 0.0 - offset = -float(number)/2.0 - - # Loop over number of observations - for i in range(number): - id = "%6.6d" % i - obs_cta = self.set_one_obs(id, offset) - offset += 1.0 - tstart += 1830.0 - obs.append(obs_cta) - - # Return observation container - return obs - - # Setup one observation - def set_one_obs(self, id, offset, \ - tstart=0.0, duration=1800.0, deadc=0.95, \ - emin=0.1, emax=100.0, rad=5.0, \ - irf="cta_dummy_irf", caldb="dummy"): + # Test observation container + self.test_value(ctobssim.obs().size(), nobs, + "There is one observation") + for i in range(ctobssim.obs().size()): + obs = gammalib.GCTAObservation(ctobssim.obs()[i]) + pnt = obs.pointing() + self.test_assert(obs.instrument() == "CTA", + "Observation is CTA observation") + self.test_value(obs.ontime(), 1800.0, 1.0e-6, + "Ontime is 1800 sec") + self.test_value(obs.livetime(), 1710.0, 1.0e-6, + "Livetime is 1710 sec") + self.test_value(pnt.dir().ra_deg(), pnts[i]['ra'], 1.0e-6, + "Pointing Right Ascension is "+str(pnts[i]['ra'])+" deg") + self.test_value(pnt.dir().dec_deg(), pnts[i]['dec'], 1.0e-6, + "Pointing Declination is "+str(pnts[i]['dec'])+" deg") + + # Return + return + + # Check event list + def _test_list(self, list, nevents): """ - Setup one observation for test purposes. + Test content of event list." + + Args: + list: Event list + nevents: Expected number of events. """ - # Allocate CTA observation - obs_cta = gammalib.GCTAObservation() - - # Set calibration database - db = gammalib.GCaldb() - if (gammalib.dir_exists(self.caldb)): - db.rootdir(self.caldb) - else: - db.open("cta", self.caldb) - - # Set pointing direction - pntdir = gammalib.GSkyDir() - pntdir.radec_deg(83.63, 22.01+offset) - pnt = gammalib.GCTAPointing() - pnt.dir(pntdir) - obs_cta.pointing(pnt) - - # Set ROI - roi = gammalib.GCTARoi() - instdir = gammalib.GCTAInstDir() - instdir.dir(pntdir) - roi.centre(instdir) - roi.radius(rad) - - # Set GTI - gti = gammalib.GGti() - gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) - - # Set energy boundaries - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ - gammalib.GEnergy(emax, "TeV")) - - # Allocate event list - events = gammalib.GCTAEventList() - events.roi(roi) - events.gti(gti) - events.ebounds(ebounds) - obs_cta.events(events) - - # Set instrument response - obs_cta.response(self.irf, db) - - # Set ontime, livetime, and deadtime correction factor - obs_cta.ontime(duration) - obs_cta.livetime(duration*deadc) - obs_cta.deadc(deadc) - obs_cta.id(id) - - # Return CTA observation - return obs_cta + # Test event list + self.test_value(list.size(), nevents, str(nevents)+" elements") + self.test_value(list.number(), nevents, str(nevents)+" events") + + # Return + return diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py index 9c57f1a..c3266af 100755 --- a/test/test_python_cscripts.py +++ b/test/test_python_cscripts.py @@ -34,7 +34,7 @@ def test(installed=False): Perform unit testing for Python interface. """ # Allocate test suite container - suites = gammalib.GTestSuites("cscripts unit testing") + suites = gammalib.GTestSuites("cscripts Python module unit testing") # Allocate test suites and append them to the container suite_cscaldb = test_cscaldb.Test() diff --git a/test/test_python_ctools.py b/test/test_python_ctools.py index 3d4aacd..1f72500 100755 --- a/test/test_python_ctools.py +++ b/test/test_python_ctools.py @@ -48,7 +48,7 @@ def test(installed=False): Perform unit testing for Python interface. """ # Allocate test suite container - suites = gammalib.GTestSuites("ctools unit testing") + suites = gammalib.GTestSuites("ctools Python module unit testing") # Allocate test suites and append them to the container suite_ctobssim = test_ctobssim.Test()