diff --git a/cscripts/csfindobs.py b/cscripts/csfindobs.py index d422279..e163020 100755 --- a/cscripts/csfindobs.py +++ b/cscripts/csfindobs.py @@ -30,32 +30,34 @@ import json # =============== # class csfindobs(ctools.cscript): """ - This class can be used to find observations from an IACT data store. + Find observations from an IACT data store. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "csfindobs" - self.version = "1.0.0" - self.verbose = False - self.datapath = os.getenv("VHEFITS","") + self._name = "csfindobs" + self._version = "1.1.0" + self._verbose = False + self._datapath = os.getenv("VHEFITS","") # Make sure that parfile exists - file = self.parfile() + self._parfile() # Initialise application if len(argv) == 0: - ctools.cscript.__init__(self, self.name, self.version) + ctools.cscript.__init__(self, self._name, self._version) elif len(argv) ==1: - ctools.cscript.__init__(self, self.name, self.version, *argv) + 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) + self._log_header() + self._log.date(True) # Return return @@ -65,19 +67,20 @@ class csfindobs(ctools.cscript): Destructor. """ # Write separator into logger - if self.logTerse(): - self.log("\n") + if self._logTerse(): + self._log("\n") # Return return - def parfile(self): + 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. + default parfile. This kluge avoids shipping the cscript with a + parfile. """ # Set parfile name - parfile = self.name+".par" + parfile = self._name+".par" try: pars = gammalib.GApplicationPars(parfile) @@ -87,7 +90,7 @@ class csfindobs(ctools.cscript): # Create default parfile pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("datapath","s","a",self.datapath,"","","Path were data are located")) + 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")) @@ -103,89 +106,105 @@ class csfindobs(ctools.cscript): # Return return - def get_parameters(self): + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get parameters - if self.datapath == "": - self.datapath = self["datapath"].string() + if self._datapath == "": + self._datapath = self["datapath"].string() # Expand environment - self.datapath = gammalib.expand_env(self.datapath) + self._datapath = gammalib.expand_env(self._datapath) # Get production name - self.m_prodname = self["prodname"].string() + self._prodname = self["prodname"].string() # Master index file name - self.m_master_indx = self["master_indx"].string() + self._master_indx = self["master_indx"].string() # Intiialise flag if spatial selection is required - self.select_radec = True + self._select_radec = True # Initialise invalid radius - self.m_radius = 0.0 + self._radius = 0.0 # Check for vailidity of spatial parameters - if self["ra"].is_valid() and self["dec"].is_valid() and self["rad"].is_valid(): + if (self["ra"].is_valid() and + self["dec"].is_valid() and + self["rad"].is_valid()): # Read spatial parameters - self.m_ra = self["ra"].real() - self.m_dec = self["dec"].real() - self.m_radius = self["rad"].real() + self._ra = self["ra"].real() + self._dec = self["dec"].real() + self._radius = self["rad"].real() else: - self.select_radec = False + self._select_radec = False # Check Radius for validity - if self.m_radius <= 0.0: - self.select_radec = False + if self._radius <= 0.0: + self._select_radec = False # Read other parameters - self.m_quality = self["min_qual"].integer() - self.m_expression = self["expression"].string() + self._quality = self["min_qual"].integer() + self._expression = self["expression"].string() # Check for validity of expression - if self.m_expression == "NONE" or self.m_expression == "INDEF": - self.m_expression = "" + if self._expression == "NONE" or self._expression == "INDEF": + self._expression = "" # Read outfile - self.m_outfile = self["outfile"].filename() + self._outfile = self["outfile"].filename() # Open master index file and look for prodname - master_file = os.path.join(self.datapath, self.m_master_indx) + master_file = os.path.join(self._datapath, self._master_indx) if not os.path.isfile(master_file): - raise RuntimeError("FITS data store not available. No master index file found at \""+master_file+"\". Make sure the file is copied from the server and your datapath is set correctly.") + raise RuntimeError("FITS data store not available. "+ + "No master index file found at \""+ + master_file+"\". Make sure the file "+ + "is copied from the server and your "+ + "datapath is set correctly.") # Open and load JSON file json_data = open(master_file).read() data = json.loads(json_data) if not "datasets" in data: - raise RuntimeError("Key \"datasets\" not available in master index file.") + raise RuntimeError("Key \"datasets\" not available in "+ + "master index file.") # Get configurations configs = data["datasets"] # Initialise obs index file - self.m_obs_index = "" + self._obs_index = "" # Get HDUs for config in configs: - if self.m_prodname == config["name"]: - self.m_obs_index = str(os.path.join(self.datapath, config["obsindx"])) + if self._prodname == config["name"]: + self._obs_index = str(os.path.join(self._datapath, + config["obsindx"])) break # Check HDUs - if self.m_obs_index == "": - raise RuntimeError("FITS data store \""+self.m_prodname+"\" not available. Run csiactdata to get a list of available storage names") - filename = gammalib.GFilename(self.m_obs_index+"[OBS_INDEX]") + if self._obs_index == "": + raise RuntimeError("FITS data store \""+self._prodname+ + "\" not available. Run csiactdata to get "+ + "a list of available storage names") + filename = gammalib.GFilename(self._obs_index+"[OBS_INDEX]") if not filename.is_fits(): - raise RuntimeError("Observation index file \""+self.m_obs_index+"[OBS_INDEX]\" for FITS data store \""+self.m_prodname+"\" not available. Check your master index file or run csiactdata to get a list of available storage names.") + raise RuntimeError("Observation index file \""+ + self._obs_index+ + "[OBS_INDEX]\" for FITS data store \""+ + self._prodname+ + "\" not available. "+ + "Check your master index file or run "+ + "csiactdata to get a list of available "+ + "storage names.") # Set other members - self.m_log = False # Logging in client tools - self.m_debug = False # Debugging in client tools - self.m_clobber = self["clobber"].boolean() + self._debug = False # Debugging in client tools + self._clobber = self["clobber"].boolean() # Return return @@ -206,7 +225,7 @@ class csfindobs(ctools.cscript): self.run() # Save residual map - self.save(self.m_outfile) + self.save(self._outfile) # Return return @@ -216,40 +235,41 @@ class csfindobs(ctools.cscript): Run the script. """ # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + if self._logDebug(): + self._log.cout(True) # Get parameters - self.get_parameters() + self._get_parameters() # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + if self._logTerse(): + self._log_parameters() + self._log("\n") # Initialise run list - self.runs = [] + self._runs = [] # Initialise selction expression expr = "" # Add spatial expression if possible - if self.select_radec: - expr += "ANGSEP("+str(self.m_ra)+","+str(self.m_dec)+",RA_PNT,DEC_PNT)<="+str(self.m_radius) + if self._select_radec: + expr += "ANGSEP("+str(self._ra)+","+str(self._dec)+ \ + ",RA_PNT,DEC_PNT)<="+str(self._radius) # Add connector if expression is empty if len(expr): expr += "&&" # Add quality expression - expr += "QUALITY<="+str(self.m_quality) + expr += "QUALITY<="+str(self._quality) # Add user expression - if self.m_expression == "NONE" or len(self.m_expression) > 0: - expr += "&&"+self.m_expression + if self._expression == "NONE" or len(self._expression) > 0: + expr += "&&"+self._expression # Initialise filename - openfile = self.m_obs_index+"[OBS_INDEX]["+expr+"]" + openfile = self._obs_index+"[OBS_INDEX]["+expr+"]" # Open file fits = gammalib.GFits(openfile) @@ -262,21 +282,21 @@ class csfindobs(ctools.cscript): # Read observation IDs into array for i in range(obs_indx.nrows()): - self.runs.append(obs_indx["OBS_ID"][i]) + self._runs.append(obs_indx["OBS_ID"][i]) # Dump expression - if self.logNormal(): - self.log("\n") - self.log.header3("Expression") - self.log(expr) - self.log("\n") + if self._logNormal(): + self._log("\n") + self._log.header3("Expression") + self._log(expr) + self._log("\n") # Dump number of observations - if self.logTerse(): - self.log("\n") - self.log.header3("Observations") - self.log("Found "+str(len(self.runs))+" observations") - self.log("\n") + if self._logTerse(): + self._log("\n") + self._log.header3("Observations") + self._log("Found "+str(len(self._runs))+" observations") + self._log("\n") # Close file fits.close() @@ -285,34 +305,35 @@ class csfindobs(ctools.cscript): return def save(self, outfile): - + """ + Save. + """ # Check for clobber - if outfile.exists() and not self.clobber(): - self.log("File "+outfile+" already exists. Set clobber=yes to overwrite file") - self.log("\n") + if outfile.exists() and not self._clobber(): + self._log("File \""+outfile+"\" already exists. ") + self._log("Set clobber=yes to overwrite file.") + self._log("\n") else: # Write runs to file f = open(repr(outfile),"w") - for run in self.runs: + for run in self._runs: f.write(str(run)+" \n") f.close() # Return return - + + # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Find IACT observations. - """ + # Create instance of application app = csfindobs(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() - diff --git a/cscripts/csiactcopy.py b/cscripts/csiactcopy.py index 2db33bf..5ea41de 100755 --- a/cscripts/csiactcopy.py +++ b/cscripts/csiactcopy.py @@ -39,24 +39,24 @@ class csiactcopy(ctools.cscript): Constructor. """ # Set name - self.name = "csiactcopy" - self.version = "1.1.0" - self.datapath = os.getenv("VHEFITS","") + self._name = "csiactcopy" + self._version = "1.1.0" + self._datapath = os.getenv("VHEFITS","") # Make sure that parfile exists - file = self.parfile() + self._parfile() # Initialise application if len(argv) == 0: - ctools.cscript.__init__(self, self.name, self.version) + ctools.cscript.__init__(self, self._name, self._version) elif len(argv) == 1: - ctools.cscript.__init__(self, self.name, self.version, *argv) + 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) + self._log_header() + self._log.date(True) # Return return @@ -69,14 +69,14 @@ class csiactcopy(ctools.cscript): # Return return - def parfile(self): + 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" + parfile = self._name+".par" try: pars = gammalib.GApplicationPars(parfile) @@ -97,23 +97,23 @@ class csiactcopy(ctools.cscript): # Return return - def get_parameters(self): + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get Parameters - self.m_remote_master = self["remote_master"].filename() - if not self.m_remote_master.exists(): + self._remote_master = self["remote_master"].filename() + if not self._remote_master.exists(): raise RuntimeError("*** ERROR: FITS data not available. No " "master index file found in \""+ - self.m_remote_master+"\". Make sure remote " + self._remote_master+"\". Make sure remote " "file system is properly mounted.") # Get parameters - self.m_prodname = self["prodname"].string() - self.m_outpath = gammalib.expand_env(self["outpath"].string()) - self.m_runlist = self["runlist"].filename() + self._prodname = self["prodname"].string() + self._outpath = gammalib.expand_env(self["outpath"].string()) + self._runlist = self["runlist"].filename() # Return return @@ -133,76 +133,76 @@ class csiactcopy(ctools.cscript): Run the script. """ # Switch screen logging on - self.log.cout(True) + self._log.cout(True) # Get parameters - self.get_parameters() + self._get_parameters() # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + if self._logTerse(): + self._log_parameters() + self._log("\n") # Make destination directory if not available - if not os.path.isdir(self.m_outpath): - os.makedirs(self.m_outpath) + if not os.path.isdir(self._outpath): + os.makedirs(self._outpath) # Log output - if self.logTerse(): - self.log("\n") - self.log.header2("Runlist information") + if self._logTerse(): + self._log("\n") + self._log.header2("Runlist information") # Initialise runlist array - self.runlist = [] + self._runs = [] # Check if runlist file is available - if self.m_runlist.exists(): - runfile = open(self.m_runlist) + if self._runlist.exists(): + runfile = open(self._runlist) for line in runfile.readlines(): if len(line) == 0: continue if line[0] == '#': continue if len(line.split()) > 0: - self.runlist.append(int(line.split()[0])) + self._runs.append(int(line.split()[0])) runfile.close() # Logging - if self.logTerse(): - self.log.parformat("Number of observations") - self.log(str(len(self.runlist))) - self.log("\n") + if self._logTerse(): + self._log.parformat("Number of observations") + self._log(str(len(self._runs))) + self._log("\n") # Copy all data if runlist file is "NONE" - elif self.m_runlist == "NONE": - if self.logTerse(): - self.log("Copy all available data") - self.log("\n") + elif self._runlist == "NONE": + if self._logTerse(): + self._log("Copy all available data") + self._log("\n") # Raise exception if file not valid else: raise RuntimeError("*** ERROR: Runlist file \""+ - self.m_runlist+"\" not available.") + self._runlist+"\" not available.") # Check for availability of remote master file - if not self.m_remote_master.exists(): + if not self._remote_master.exists(): raise RuntimeError("*** ERROR: Remote master file \""+ - self.m_remote_master+"\" does not exist.") + self._remote_master+"\" does not exist.") # Retrieve json data from remote master - json_data = open(self.m_remote_master.url()).read() + json_data = open(self._remote_master.url()).read() data = json.loads(json_data) if not "datasets" in data: raise RuntimeError("*** ERROR: Key \"datasets\" not available " "in remote master index file \""+ - self.m_remote_master+"\".") + self._remote_master+"\".") # Get array of configurations configs = data["datasets"] # Get remote paths - #self.remote_base = os.path.dirname(self.m_remote_master.url()) - self.remote_base = self.m_remote_master.path() + #self._remote_base = os.path.dirname(self._remote_master.url()) + self._remote_base = self._remote_master.path() # Initialise flag if prod has been found has_prod = False @@ -211,36 +211,36 @@ class csiactcopy(ctools.cscript): files = set() # Logging - if self.logTerse(): - self.log("\n") - self.log.header2("Loop over remote configs") + if self._logTerse(): + self._log("\n") + self._log.header2("Loop over remote configs") # Loop over configs for config in configs: # Logging - if self.logVerbose(): - self.log.header2(str(config["name"])) + if self._logVerbose(): + self._log.header2(str(config["name"])) # Check if prodname was found - if config["name"] == self.m_prodname: + if config["name"] == self._prodname: # Inidicate that prodname has been found has_prod = True # Build path of index files - remote_hdu = os.path.join(self.remote_base, config["hduindx"]) - remote_obs = os.path.join(self.remote_base, config["obsindx"]) + remote_hdu = os.path.join(self._remote_base, config["hduindx"]) + remote_obs = os.path.join(self._remote_base, config["obsindx"]) # Log information - if self.logTerse(): - self.log.header3("Remote config \""+self.m_prodname+"\"") - self.log.parformat("HDU index") - self.log(str(remote_hdu)) - self.log("\n") - self.log.parformat("Observation index") - self.log(str(remote_obs)) - self.log("\n") + if self._logTerse(): + self._log.header3("Remote config \""+self._prodname+"\"") + self._log.parformat("HDU index") + self._log(str(remote_hdu)) + self._log("\n") + self._log.parformat("Observation index") + self._log(str(remote_obs)) + self._log("\n") # Open remote HDU index file fits = gammalib.GFits(str(remote_hdu)) @@ -259,12 +259,12 @@ class csiactcopy(ctools.cscript): remote_ids.add(table["OBS_ID"][row]) # Log runs that are not available remotely - for run in self.runlist: + for run in self._runs: if not run in remote_ids: - if self.logNormal(): - self.log("Skip observation "+str(run)+": ") - self.log("ID not available remotely") - self.log("\n") + if self._logNormal(): + self._log("Skip observation "+str(run)+": ") + self._log("ID not available remotely") + self._log("\n") # Loop over remote HDU index file for row in range(table.nrows()): @@ -279,10 +279,10 @@ class csiactcopy(ctools.cscript): continue # Check if we need to consider an input runlist - if len(self.runlist): + if len(self._runs): # Check if obs id is in runlist - if obs_id in self.runlist: + if obs_id in self._runs: # Get filename fname = os.path.join(os.path.dirname(remote_hdu), file_dir, file_name) @@ -311,44 +311,44 @@ class csiactcopy(ctools.cscript): cp_size += table["SIZE"][row] # Log file information - if self.logTerse(): - self.log("\n") - self.log.header2("File information") - self.log.parformat("Number of files") - self.log(str(len(files))) - self.log("\n") + if self._logTerse(): + self._log("\n") + self._log.header2("File information") + self._log.parformat("Number of files") + self._log(str(len(files))) + self._log("\n") if has_size: - self.log.parformat("Size") - self.log("%.2f"%(float(cp_size)*1e-6)+" MB") - self.log("\n") - if self.logVerbose(): - self.log("\n") - self.log.header3("File names") + self._log.parformat("Size") + self._log("%.2f"%(float(cp_size)*1e-6)+" MB") + self._log("\n") + if self._logVerbose(): + self._log("\n") + self._log.header3("File names") for filename in files: - self.log(str(filename)+"\n") + self._log(str(filename)+"\n") # Close HDU index file fits.close() # If prodname is not found just log that we skip the config else: - if self.logExplicit(): - self.log("\n") - self.log.header3("Skipping config \""+str(config["name"])+"\"") - self.log("\n") + if self._logExplicit(): + self._log("\n") + self._log.header3("Skipping config \""+str(config["name"])+"\"") + self._log("\n") # Raise Exception if prodname was not found if not has_prod: - msg = "*** ERROR: FITS production \""+self.m_prodname+\ + msg = "*** ERROR: FITS production \""+self._prodname+\ "\" not available. Available productions are:\n" for config in configs: msg += " - "+config["name"]+"\n" raise RuntimeError(msg) # Logging - if self.logNormal(): - self.log("\n") - self.log.header1("Copying files") + if self._logNormal(): + self._log("\n") + self._log.header1("Copying files") # intialise counter k = 0 @@ -356,11 +356,11 @@ class csiactcopy(ctools.cscript): # initialise values for logging last_fraction = 0.0 fraction_increment = 20.0 - if self.logNormal(): + if self._logNormal(): fraction_increment = 10.0 - if self.logTerse(): + if self._logTerse(): fraction_increment = 5.0 - if self.logExplicit(): + if self._logExplicit(): fraction_increment = 2.0 # initialise logging properties @@ -373,14 +373,14 @@ class csiactcopy(ctools.cscript): # Log progress fraction = float(k) / float(len(files)) * 100.0 while fraction > last_fraction: - if self.logNormal() and not self.logVerbose(): - self.log.parformat("Status") - self.log(str(int(last_fraction))+"%") - self.log("\n") + if self._logNormal() and not self._logVerbose(): + self._log.parformat("Status") + self._log(str(int(last_fraction))+"%") + self._log("\n") last_fraction += fraction_increment # Copy file - filesize = self.copy(filename, self.clobber()) + filesize = self.copy(filename, self._clobber()) if filesize > 0.0: total_size += filesize @@ -390,51 +390,51 @@ class csiactcopy(ctools.cscript): k += 1 # Logging - if self.logNormal() and not self.logVerbose(): - self.log.parformat("Status") - self.log("Finished") - self.log("\n") - if self.logTerse(): - self.log("\n") - self.log.header1("Updating index files") + if self._logNormal() and not self._logVerbose(): + self._log.parformat("Status") + self._log("Finished") + self._log("\n") + if self._logTerse(): + self._log("\n") + self._log.header1("Updating index files") # Build local index file names - local_hdu = os.path.join(self.m_outpath, os.path.relpath(remote_hdu, self.remote_base)) - local_obs = os.path.join(self.m_outpath, os.path.relpath(remote_obs, self.remote_base)) + local_hdu = os.path.join(self._outpath, os.path.relpath(remote_hdu, self._remote_base)) + local_obs = os.path.join(self._outpath, os.path.relpath(remote_obs, self._remote_base)) # If we have a runlist then merge index file - if len(self.runlist): + if len(self._runs): # Logging - if self.logTerse(): - self.log("\n") - self.log.header3("HDU index") + if self._logTerse(): + self._log("\n") + self._log.header3("HDU index") # Merge remote index files with local files - self.merge(local_hdu, remote_hdu, "HDU_INDEX", self.clobber()) + self.merge(local_hdu, remote_hdu, "HDU_INDEX", self._clobber()) - if self.logTerse(): - self.log("\n") - self.log.header3("OBS index") + if self._logTerse(): + self._log("\n") + self._log.header3("OBS index") # Merge remote index files with local files - self.merge(local_obs, remote_obs, "OBS_INDEX", self.clobber()) + self.merge(local_obs, remote_obs, "OBS_INDEX", self._clobber()) else: # If all files were copied, just copy index files too - self.copy(remote_hdu, self.clobber()) - self.copy(remote_obs, self.clobber()) + self.copy(remote_hdu, self._clobber()) + self.copy(remote_obs, self._clobber()) - if self.logTerse(): - self.log("\n") - self.log.header3("Master index file") + if self._logTerse(): + self._log("\n") + self._log.header3("Master index file") # Adding prodname to local master - localmaster = os.path.join(self.m_outpath, "master.json") + localmaster = os.path.join(self._outpath, "master.json") # If no local master is found, copy master over first if not os.path.isfile(localmaster): - self.copy(self.m_remote_master.url(), self.clobber()) + self.copy(self._remote_master.url(), self._clobber()) # Load local master json_data = open(localmaster).read() @@ -452,37 +452,37 @@ class csiactcopy(ctools.cscript): for config in configs: # Get hdu and obs index files - hdu = os.path.join(self.m_outpath, config["hduindx"]) - obs = os.path.join(self.m_outpath, config["obsindx"]) + hdu = os.path.join(self._outpath, config["hduindx"]) + obs = os.path.join(self._outpath, config["obsindx"]) # Check if index files are available if not (gammalib.GFilename(str(hdu)).is_fits() and gammalib.GFilename(str(obs)).is_fits()): - if self.logTerse(): - self.log("Removing \""+str(config["name"])+"\" ") - self.log("(not available)") - self.log("\n") + if self._logTerse(): + self._log("Removing \""+str(config["name"])+"\" ") + self._log("(not available)") + self._log("\n") else: # Append config if valid newconfigs.append(config) - if self.logTerse(): - self.log("Keeping \""+str(config["name"])+"\"") - self.log("\n") + if self._logTerse(): + self._log("Keeping \""+str(config["name"])+"\"") + self._log("\n") # Signals that downloaded config is available - if config["name"] == self.m_prodname: + if config["name"] == self._prodname: has_config = True # Create new entry if config was not available if not has_config: newdict = dict.fromkeys(["name","hduindx","obsindx"]) - newdict["name"] = self.m_prodname - newdict["hduindx"] = os.path.relpath(local_hdu, self.m_outpath) - newdict["obsindx"] = os.path.relpath(local_obs, self.m_outpath) + newdict["name"] = self._prodname + newdict["hduindx"] = os.path.relpath(local_hdu, self._outpath) + newdict["obsindx"] = os.path.relpath(local_obs, self._outpath) newconfigs.append(newdict) - if self.logTerse(): - self.log("Adding \""+str(newdict["name"])+"\"") - self.log("\n") + if self._logTerse(): + self._log("Adding \""+str(newdict["name"])+"\"") + self._log("\n") # Write new json master file f = open(localmaster, "w") @@ -491,27 +491,27 @@ class csiactcopy(ctools.cscript): f.close() # Log summary - if self.logNormal(): - self.log("\n") - self.log.header1("Summary") - self.log.parformat("Data files found") - self.log(str(len(files))) - self.log("\n") - self.log.parformat("Data files copied") - self.log(str(n_copied)) - self.log("\n") - self.log.parformat("Size") - self.log("%.2f"%(float(total_size)*1e-6)+" MB") - self.log("\n") - self.log.parformat("Local configs") - self.log(str(len(newconfigs))) - self.log("\n") - if self.logTerse(): - self.log("\n") - self.log.header3("Content of master index") + if self._logNormal(): + self._log("\n") + self._log.header1("Summary") + self._log.parformat("Data files found") + self._log(str(len(files))) + self._log("\n") + self._log.parformat("Data files copied") + self._log(str(n_copied)) + self._log("\n") + self._log.parformat("Size") + self._log("%.2f"%(float(total_size)*1e-6)+" MB") + self._log("\n") + self._log.parformat("Local configs") + self._log(str(len(newconfigs))) + self._log("\n") + if self._logTerse(): + self._log("\n") + self._log.header3("Content of master index") for config in newconfigs: - self.log(str(config["name"])) - self.log("\n") + self._log(str(config["name"])) + self._log("\n") # Return return @@ -521,43 +521,44 @@ class csiactcopy(ctools.cscript): copy file to m_outpath """ # Get file destination - destination = os.path.join(self.m_outpath,os.path.relpath(source, self.remote_base)) + 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") + 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") + 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") + if self._logVerbose(): + self._log("Skip (same file)") + self._log("\n") else: # Create directory if not existent @@ -572,9 +573,9 @@ class csiactcopy(ctools.cscript): filesize = os.stat(destination).st_size # Logging - if self.logVerbose(): - self.log("Done!") - self.log("\n") + if self._logVerbose(): + self._log("Done!") + self._log("\n") # Return return filesize @@ -586,22 +587,22 @@ class csiactcopy(ctools.cscript): 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 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") + self._log("Remote") else: - self.log("Local") - self.log("\n") + self._log("Local") + self._log("\n") # Check if local fits file is available if os.path.isfile(localfits): @@ -621,9 +622,9 @@ class csiactcopy(ctools.cscript): # Create selection expression for opening remote fits file selection = "" - if len(self.runlist): + if len(self._runs): selection += "["+hduname+"][" - for run in self.runlist: + for run in self._runs: selection += "OBS_ID=="+str(run) selection += "||" selection = selection[:-2] @@ -642,13 +643,13 @@ class csiactcopy(ctools.cscript): 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") + 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 @@ -706,13 +707,13 @@ class csiactcopy(ctools.cscript): # Add rmeote HDUs local_hdu.insert_rows(old_local_rows, remote_hdu.nrows()) - if self.logTerse(): + if self._logTerse(): if clobber: - self.log.parformat("Removed local rows") + self._log.parformat("Removed local rows") else: - self.log.parformat("Skipped remote rows") - self.log(str(removed_rows)) - self.log("\n") + self._log.parformat("Skipped remote rows") + self._log(str(removed_rows)) + self._log("\n") # Loop over columns for i in range(local_hdu.ncols()): @@ -735,14 +736,12 @@ class csiactcopy(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Copy IACT data - """ + # Create instance of application app = csiactcopy(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csiactdata.py b/cscripts/csiactdata.py index c396292..c8524e3 100755 --- a/cscripts/csiactdata.py +++ b/cscripts/csiactdata.py @@ -40,25 +40,25 @@ class csiactdata(ctools.cscript): Constructor. """ # Set name - self.name = "csiactdata" - self.version = "1.0.0" - self.datapath = os.getenv("VHEFITS","") - self.prodnames = [] + self._name = "csiactdata" + self._version = "1.1.0" + self._datapath = os.getenv("VHEFITS","") + self._prodnames = [] # Make sure that parfile exists - file = self.parfile() + file = self._parfile() # Initialise application if len(argv) == 0: - ctools.cscript.__init__(self, self.name, self.version) + ctools.cscript.__init__(self, self._name, self._version) elif len(argv) == 1: - ctools.cscript.__init__(self, self.name, self.version, *argv) + 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) + self._log_header() + self._log.date(True) # Return return @@ -70,14 +70,14 @@ class csiactdata(ctools.cscript): # Return return - def parfile(self): + 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" + parfile = self._name+".par" try: pars = gammalib.GApplicationPars(parfile) @@ -87,7 +87,7 @@ class csiactdata(ctools.cscript): # Create default parfile pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("datapath","s","a",self.datapath,"","","Path were data is located")) + 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")) @@ -102,19 +102,19 @@ class csiactdata(ctools.cscript): """ # Get datapath if not already set - if self.datapath == "": - self.datapath = self["datapath"].string() + if self._datapath == "": + self._datapath = self["datapath"].string() # Expand environment - self.datapath = gammalib.expand_env(self.datapath) + self._datapath = gammalib.expand_env(self._datapath) # Get filename of master index file - self.m_master_indx = self["master_indx"].string() - self.m_master_file = os.path.join(self.datapath, self.m_master_indx) + self._master_indx = self["master_indx"].string() + self._master_file = os.path.join(self._datapath, self._master_indx) # Check for presence of master index file - if not os.path.isfile(self.m_master_file): - raise RuntimeError("Master index file \""+self.m_master_file+"\" not found. Use hidden parameter \"master_indx\" to specifiy a different filename.") + if not os.path.isfile(self._master_file): + raise RuntimeError("Master index file \""+self._master_file+"\" not found. Use hidden parameter \"master_indx\" to specifiy a different filename.") # Return return @@ -134,58 +134,58 @@ class csiactdata(ctools.cscript): Run the script. """ # Switch screen logging on - self.log.cout(True) + self._log.cout(True) # Get parameters - self.get_parameters() + self._get_parameters() # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + if self._logTerse(): + self._log_parameters() + self._log("\n") - self.log("\n") - self.log.header1("Data storage entries") - self.log("\n") - self.log.parformat("Master index file") - self.log(self.m_master_file) - self.log("\n") - self.log("\n") + self._log("\n") + self._log.header1("Data storage entries") + self._log("\n") + self._log.parformat("Master index file") + self._log(self._master_file) + self._log("\n") + self._log("\n") # Open and load JSON file - json_data = open(self.m_master_file).read() + json_data = open(self._master_file).read() data = json.loads(json_data) configs = data["datasets"] # Initialise array for available names - self.prodnames = [] + self._prodnames = [] # Loop over configs and display unavailable storage first for config in configs: # Create hdu and obs index files - hdu = os.path.join(self.datapath, config["hduindx"]) - obs = os.path.join(self.datapath, config["obsindx"]) + hdu = os.path.join(self._datapath, config["hduindx"]) + obs = os.path.join(self._datapath, config["obsindx"]) filename_hdu = gammalib.GFilename(str(hdu)+"[HDU_INDEX]") filename_obs = gammalib.GFilename(str(obs)+"[OBS_INDEX]") # Check if index files are available if not (filename_hdu.is_fits() and filename_obs.is_fits()): - self.log.parformat(str(config["name"])) - self.log("Not available") - self.log("\n") + self._log.parformat(str(config["name"])) + self._log("Not available") + self._log("\n") - self.log("\n") - self.log.header2("Available data configs") + self._log("\n") + self._log.header2("Available data configs") # Loop over configs and log available configs for config in configs: # Create hdu and obs index files - hdu = os.path.join(self.datapath, config["hduindx"]) - obs = os.path.join(self.datapath, config["obsindx"]) + hdu = os.path.join(self._datapath, config["hduindx"]) + obs = os.path.join(self._datapath, config["obsindx"]) filename_hdu = gammalib.GFilename(str(hdu)+"[HDU_INDEX]") filename_obs = gammalib.GFilename(str(obs)+"[OBS_INDEX]") @@ -194,30 +194,30 @@ class csiactdata(ctools.cscript): if (filename_hdu.is_fits() and filename_obs.is_fits()): # Log data information - self.log("\n") - self.log.header3(str(config["name"])) + self._log("\n") + self._log.header3(str(config["name"])) # Print important informations first - self.log.parformat("Name") - self.log(str(config["name"])) - self.log("\n") - self.log.parformat("Observation index") - self.log(str(config["obsindx"])) - self.log("\n") - self.log.parformat("HDU index") - self.log(str(config["hduindx"])) - self.log("\n") + self._log.parformat("Name") + self._log(str(config["name"])) + self._log("\n") + self._log.parformat("Observation index") + self._log(str(config["obsindx"])) + self._log("\n") + self._log.parformat("HDU index") + self._log(str(config["hduindx"])) + self._log("\n") # Append to available names - self.prodnames.append(str(config["name"])) + self._prodnames.append(str(config["name"])) # Print additional information for key in config: if key in ["name","hduindx","obsindx"]: continue - self.log.parformat(str(key)) - self.log(str(config[key])) - self.log("\n") + self._log.parformat(str(key)) + self._log(str(config[key])) + self._log("\n") # Return return @@ -228,20 +228,19 @@ class csiactdata(ctools.cscript): """ # Return - return self.prodnames + return self._prodnames + # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Shows data storage information - """ + # Create instance of application app = csiactdata(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csiactobs.py b/cscripts/csiactobs.py index af1805a..f4ae29c 100755 --- a/cscripts/csiactobs.py +++ b/cscripts/csiactobs.py @@ -42,31 +42,31 @@ class csiactobs(ctools.cscript): Constructor. """ # Set name and version - self.name = "csiactobs" - self.version = "1.1.0" + self._name = "csiactobs" + self._version = "1.1.0" # Initialise some members - self.m_ebounds = gammalib.GEbounds() - self.datapath = os.getenv("VHEFITS","") - self.inmodels = None - self.xml = gammalib.GXml() - self.models = gammalib.GModels() - self.m_runlist = [] + self._ebounds = gammalib.GEbounds() + self._datapath = os.getenv("VHEFITS","") + self._inmodels = None + self._xml = gammalib.GXml() + self._models = gammalib.GModels() + self._runlist = [] # Make sure that parfile exists - file = self.parfile() + file = self._parfile() # Initialise application if len(argv) == 0: - ctools.cscript.__init__(self, self.name, self.version) + ctools.cscript.__init__(self, self._name, self._version) elif len(argv) == 1: - ctools.cscript.__init__(self, self.name, self.version, *argv) + 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) + self._log_header() + self._log.date(True) # Return return @@ -78,13 +78,13 @@ class csiactobs(ctools.cscript): # Return return - def parfile(self): + 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" + parfile = self._name+".par" try: pars = gammalib.GApplicationPars(parfile) @@ -94,7 +94,7 @@ class csiactobs(ctools.cscript): # Create default parfile pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("datapath","s","a",self.datapath,"","","Path were data is located")) + 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)")) @@ -122,74 +122,74 @@ class csiactobs(ctools.cscript): # Return return - def get_parameters(self): + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ - if self.datapath == "": - self.datapath = self["datapath"].string() + if self._datapath == "": + self._datapath = self["datapath"].string() # Check for input models if not self["inmodel"].filename() == "NONE": - self.inmodels = gammalib.GModels(self["inmodel"].filename()) + self._inmodels = gammalib.GModels(self["inmodel"].filename()) # Expand environment - self.datapath = gammalib.expand_env(self.datapath) + self._datapath = gammalib.expand_env(self._datapath) # Read FITS production - self.m_prodname = self["prodname"].string() + self._prodname = self["prodname"].string() # Read runlist file if list not already filled - if len(self.m_runlist) == 0: + if len(self._runlist) == 0: # Get file name - self.m_runlistfile = self["infile"].filename() + self._runlistfile = self["infile"].filename() # Read runlist from file - runfile = open(self.m_runlistfile.url()) + runfile = open(self._runlistfile.url()) for line in runfile.readlines(): if len(line) == 0: continue if line[0] == '#': continue if len(line.split()) > 0: - self.m_runlist.append(line.split()[0]) + self._runlist.append(line.split()[0]) runfile.close() # Read number of background parameters - self.m_bkgpars = self["bkgpars"].integer() + self._bkgpars = self["bkgpars"].integer() # Output model file - self.m_outmodel = self["outmodel"].filename() + self._outmodel = self["outmodel"].filename() # Observation outfile - self.m_outobs = self["outobs"].filename() + self._outobs = self["outobs"].filename() # Master index file name - self.m_master_indx = self["master_indx"].string() + self._master_indx = self["master_indx"].string() # Read flag for background scaling factor - self.m_use_bkg_scale = self["bkg_scale"].boolean() + self._use_bkg_scale = self["bkg_scale"].boolean() # Read hierarchy of file loading - self.m_ev_hiera = self["ev_hiera"].string().split("|") - self.m_aeff_hiera = self["aeff_hiera"].string().split("|") - self.m_psf_hiera = self["psf_hiera"].string().split("|") - self.m_bkg_hiera = self["bkg_hiera"].string().split("|") - self.m_edisp_hiera = self["edisp_hiera"].string().split("|") - self.m_bkg_mod_hiera = self["bkg_mod_hiera"].string().split("|") + self._ev_hiera = self["ev_hiera"].string().split("|") + self._aeff_hiera = self["aeff_hiera"].string().split("|") + self._psf_hiera = self["psf_hiera"].string().split("|") + self._bkg_hiera = self["bkg_hiera"].string().split("|") + self._edisp_hiera = self["edisp_hiera"].string().split("|") + self._bkg_mod_hiera = self["bkg_mod_hiera"].string().split("|") # Read hidden background parameters - self.m_bkg_gauss_norm = self["bkg_gauss_norm"].real() - self.m_bkg_gauss_index = self["bkg_gauss_index"].real() - self.m_bkg_gauss_sigma = self["bkg_gauss_sigma"].real() - self.m_bkg_aeff_index = self["bkg_aeff_index"].real() - self.m_bkg_aeff_norm = self["bkg_aeff_norm"].real() - self.m_bkg_range_factor = self["bkg_range_factor"].real() + self._bkg_gauss_norm = self["bkg_gauss_norm"].real() + self._bkg_gauss_index = self["bkg_gauss_index"].real() + self._bkg_gauss_sigma = self["bkg_gauss_sigma"].real() + self._bkg_aeff_index = self["bkg_aeff_index"].real() + self._bkg_aeff_norm = self["bkg_aeff_norm"].real() + self._bkg_range_factor = self["bkg_range_factor"].real() # Open master index file and look for prodname - master_file = os.path.join(self.datapath, self.m_master_indx) + master_file = os.path.join(self._datapath, self._master_indx) if not os.path.isfile(master_file): raise RuntimeError("FITS data store not available. No master index file found. Make sure the file is copied from the server and your datapath is set correctly.") @@ -203,38 +203,37 @@ class csiactobs(ctools.cscript): configs = data["datasets"] # Initialise HDUs - self.m_hdu_index = self.m_obs_index = "" + self._hdu_index = self._obs_index = "" # Get HDUs for config in configs: - if self.m_prodname == config["name"]: - self.m_hdu_index = str(os.path.join(self.datapath, config["hduindx"])) - self.m_obs_index = str(os.path.join(self.datapath, config["obsindx"])) + if self._prodname == config["name"]: + self._hdu_index = str(os.path.join(self._datapath, config["hduindx"])) + self._obs_index = str(os.path.join(self._datapath, config["obsindx"])) break # Check HDUs - if self.m_hdu_index == "" or self.m_obs_index == "": - raise RuntimeError("*** ERROR: FITS data store \""+self.m_prodname+"\" not available. Run csiactdata to get a list of available storage names") - filename = gammalib.GFilename(self.m_hdu_index+"[HDU_INDEX]") + if self._hdu_index == "" or self._obs_index == "": + raise RuntimeError("*** ERROR: FITS data store \""+self._prodname+"\" not available. Run csiactdata to get a list of available storage names") + filename = gammalib.GFilename(self._hdu_index+"[HDU_INDEX]") if not filename.is_fits(): - raise RuntimeError("*** ERROR: HDU index file \""+self.m_hdu_index+"[HDU_INDEX]\" for FITS data store \""+self.m_prodname+"\" not available. Check your master index file or run csiactdata to get a list of available storage names.") + raise RuntimeError("*** ERROR: HDU index file \""+self._hdu_index+"[HDU_INDEX]\" for FITS data store \""+self._prodname+"\" not available. Check your master index file or run csiactdata to get a list of available storage names.") # Check for existence of "BKG_SCALE" in the observation index file if required - if self.m_use_bkg_scale: - filename = gammalib.GFilename(self.m_obs_index+"[OBS_INDEX]") + if self._use_bkg_scale: + filename = gammalib.GFilename(self._obs_index+"[OBS_INDEX]") if filename.is_fits(): - fits = gammalib.GFits(self.m_obs_index) + fits = gammalib.GFits(self._obs_index) if not fits["OBS_INDEX"].contains("BKG_SCALE"): - self.m_use_bkg_scale = False + self._use_bkg_scale = False fits.close() else: - self.m_use_bkg_scale = False + self._use_bkg_scale = False # Create base data directory from hdu index file location - self.subdir = os.path.dirname(self.m_hdu_index) - self.m_log = False # Logging in client tools - self.m_debug = False # Debugging in client tools - self.m_clobber = self["clobber"].boolean() + self._subdir = os.path.dirname(self._hdu_index) + self._debug = False # Debugging in client tools + self._clobber = self["clobber"].boolean() # Return return @@ -248,7 +247,7 @@ class csiactobs(ctools.cscript): try: # Store a copy of input argument - self.m_runlist = list(runlist) + self._runlist = list(runlist) except: @@ -258,15 +257,15 @@ class csiactobs(ctools.cscript): # Return return - def background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0): + def _background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0): # Handle constant spectral model - if index == 0.0 and self.m_bkgpars <= 1: + if index == 0.0 and self._bkgpars <= 1: spec = gammalib.GModelSpectralConst() - spec["Value"].min(prefactor / self.m_bkg_range_factor) - spec["Value"].max(prefactor * self.m_bkg_range_factor) + spec["Value"].min(prefactor / self._bkg_range_factor) + spec["Value"].max(prefactor * self._bkg_range_factor) spec["Value"].value(prefactor) - if self.m_bkgpars == 0: + if self._bkgpars == 0: spec["Value"].fix() else: spec["Value"].free() @@ -274,7 +273,7 @@ class csiactobs(ctools.cscript): else: # Create power law model - if self.m_bkgpars <= 2: + if self._bkgpars <= 2: e = gammalib.GEnergy(1.0,"TeV") spec = gammalib.GModelSpectralPlaw(prefactor, index, e) @@ -286,10 +285,10 @@ class csiactobs(ctools.cscript): spec[1].max(5.0) # Set number of free parameters - if self.m_bkgpars == 0: + if self._bkgpars == 0: spec[0].fix() spec[1].fix() - elif self.m_bkgpars == 1: + elif self._bkgpars == 1: spec[0].free() spec[1].fix() else: @@ -303,7 +302,7 @@ class csiactobs(ctools.cscript): # Create spectral model and energy values spec = gammalib.GModelSpectralNodes() - bounds = gammalib.GEbounds(self.m_bkgpars,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"), True) + bounds = gammalib.GEbounds(self._bkgpars,gammalib.GEnergy(emin,"TeV"),gammalib.GEnergy(emax,"TeV"), True) for i in range(bounds.size()): energy = bounds.elogmean(i) value = plaw.eval(energy, gammalib.GTime()) @@ -313,41 +312,41 @@ class csiactobs(ctools.cscript): if "Energy" in par.name(): par.fix() elif "Intensity" in par.name(): - par.min(value / self.m_bkg_range_factor) - par.max(value * self.m_bkg_range_factor) + par.min(value / self._bkg_range_factor) + par.max(value * self._bkg_range_factor) # Return spectrum return spec - def iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100): + def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100): # handle IrfBackground if bkgtype == "irf": prefactor = 1.0 index = 0.0 prefactor *= bkg_scale - spec = self.background_spectrum(obs_id, prefactor, index, emin, emax) + spec = self._background_spectrum(obs_id, prefactor, index, emin, emax) # Create background model instance bck = gammalib.GCTAModelIrfBackground(spec) # Set AeffBackground elif bkgtype == "aeff": - prefactor = bkg_scale * self.m_bkg_aeff_norm - spec = self.background_spectrum(obs_id, prefactor, self.m_bkg_aeff_index, emin, emax) + prefactor = bkg_scale * self._bkg_aeff_norm + spec = self._background_spectrum(obs_id, prefactor, self._bkg_aeff_index, emin, emax) # Create background model instance bck = gammalib.GCTAModelAeffBackground(spec) # Set Gaussian Background elif bkgtype == "gauss": - prefactor = bkg_scale * self.m_bkg_gauss_norm - spec = self.background_spectrum(obs_id, prefactor, self.m_bkg_gauss_index, emin, emax) - radial = gammalib.GCTAModelRadialGauss(self.m_bkg_gauss_sigma) + prefactor = bkg_scale * self._bkg_gauss_norm + spec = self._background_spectrum(obs_id, prefactor, self._bkg_gauss_index, emin, emax) + radial = gammalib.GCTAModelRadialGauss(self._bkg_gauss_sigma) bck = gammalib.GCTAModelRadialAcceptance(radial, spec) else: - sys.exit("Background type \""+self.bkgtype+"\" unsupported") + sys.exit("Background type \""+self._bkgtype+"\" unsupported") # Copy model model = bck.clone() @@ -371,7 +370,7 @@ class csiactobs(ctools.cscript): """ Returns runlist energy range """ - return self.m_ebounds + return self._ebounds def obs(self): """ @@ -381,11 +380,11 @@ class csiactobs(ctools.cscript): obs = gammalib.GObservations() # Read observations if XML is filled - if self.xml.size(): - obs.read(self.xml) + if self._xml.size(): + obs.read(self._xml) # Assign models - obs.models(self.models) + obs.models(self._models) # Return observations return obs @@ -399,7 +398,7 @@ class csiactobs(ctools.cscript): self.run() # Save residual map - self.save(self.m_outobs, self.m_clobber) + self.save(self._outobs, self._clobber) # Return return @@ -409,41 +408,41 @@ class csiactobs(ctools.cscript): Run the script. """ # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + if self._logDebug(): + self._log.cout(True) # Get parameters - self.get_parameters() + self._get_parameters() # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + if self._logTerse(): + self._log_parameters() + self._log("\n") # Clear models - self.models.clear() + self._models.clear() # Clear xml file and append an observation list - self.xml.clear() - self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\"")) - lib = self.xml.element("observation_list", 0) + self._xml.clear() + self._xml.append(gammalib.GXmlElement("observation_list title=\"observation list\"")) + lib = self._xml.element("observation_list", 0) # Log output - if self.logTerse(): - self.log("\n") - self.log.header1("Looping over "+str(len(self.m_runlist))+" runs") + if self._logTerse(): + self._log("\n") + self._log.header1("Looping over "+str(len(self._runlist))+" runs") # Initialise energy range values for logging - self.m_ebounds.clear() + self._ebounds.clear() # Loop over runs - for obs_id in self.m_runlist: + for obs_id in self._runlist: # Create selection string obs_selection = "[OBS_ID=="+str(obs_id)+"]" # Open HDU index file - hduindx = gammalib.GFits(self.m_hdu_index+"[HDU_INDEX]"+obs_selection) + hduindx = gammalib.GFits(self._hdu_index+"[HDU_INDEX]"+obs_selection) hduindx_hdu = hduindx["HDU_INDEX"] # Initialise files and hdu names @@ -458,80 +457,80 @@ class csiactobs(ctools.cscript): # Handle events index = -1 - for version in self.m_ev_hiera: + for version in self._ev_hiera: n = formats.count(version) if n > 0: index = formats.index(version) break if index >= 0: - eventfile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) + eventfile = os.path.join(self._subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) eventhdu = hduindx_hdu["HDU_NAME"][index] eventfile += "["+eventhdu+"]" if not gammalib.GFilename(eventfile).is_fits(): - if self.logTerse(): - self.log("Skipping observation "+str(obs_id)+": eventfile \""+eventfile+"\" not found\n") + if self._logTerse(): + self._log("Skipping observation "+str(obs_id)+": eventfile \""+eventfile+"\" not found\n") continue # Handle Aeff index = -1 - for version in self.m_aeff_hiera: + for version in self._aeff_hiera: n = formats.count(version) if n > 0: index = formats.index(version) break if index >= 0: - aefffile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) + aefffile = os.path.join(self._subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) aeffhdu = hduindx_hdu["HDU_NAME"][index] aefffile += "["+aeffhdu+"]" if not gammalib.GFilename(aefffile).is_fits(): - if self.logTerse(): - self.log("Skipping observation "+str(obs_id)+": effective area \""+aefffile+"\" not found\n") + if self._logTerse(): + self._log("Skipping observation "+str(obs_id)+": effective area \""+aefffile+"\" not found\n") continue # Handle psf index = -1 - for version in self.m_psf_hiera: + for version in self._psf_hiera: n = formats.count(version) if n > 0: index = formats.index(version) break if index >= 0: - psffile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) + psffile = os.path.join(self._subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) psffile += "["+hduindx_hdu["HDU_NAME"][index]+"]" if not gammalib.GFilename(psffile).is_fits(): - if self.logTerse(): - self.log("Skipping observation "+str(obs_id)+": point spread function \""+psffile+"\" not found\n") + if self._logTerse(): + self._log("Skipping observation "+str(obs_id)+": point spread function \""+psffile+"\" not found\n") continue # Handle edisp index = -1 - for version in self.m_edisp_hiera: + for version in self._edisp_hiera: n = formats.count(version) if n > 0: index = formats.index(version) break if index >= 0: - edispfile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) + edispfile = os.path.join(self._subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) edispfile += "["+hduindx_hdu["HDU_NAME"][index]+"]" if not gammalib.GFilename(edispfile).is_fits(): - if self.logTerse(): - self.log("Warning: observation "+str(obs_id)+" has no energy dispersion \""+edispfile+"\" information\n") + if self._logTerse(): + self._log("Warning: observation "+str(obs_id)+" has no energy dispersion \""+edispfile+"\" information\n") edispfile = "" # Handle background index = -1 - bkg_mod_hierarchy = list(self.m_bkg_mod_hiera) - for version in self.m_bkg_hiera: + bkg_mod_hierarchy = list(self._bkg_mod_hiera) + for version in self._bkg_hiera: n = formats.count(version) if n > 0: index = formats.index(version) break if index >= 0: - bkgfile = os.path.join(self.subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) + bkgfile = os.path.join(self._subdir, hduindx_hdu["FILE_DIR"][index], hduindx_hdu["FILE_NAME"][index]) bkgfile += "["+hduindx_hdu["HDU_NAME"][index]+"]" if not gammalib.GFilename(bkgfile).is_fits(): - if self.logTerse(): - self.log("Warning: observation "+str(obs_id)+" has no background information (file=\""+bkgfile+"\"). IRF background cannot be used\n") + if self._logTerse(): + self._log("Warning: observation "+str(obs_id)+" has no background information (file=\""+bkgfile+"\"). IRF background cannot be used\n") bkgfile = "" if "irf" in bkg_mod_hierarchy: bkg_mod_hierarchy.remove("irf") @@ -541,8 +540,8 @@ class csiactobs(ctools.cscript): # Handle background scale information if available bkg_scale = 1.0 - if self.m_use_bkg_scale: - obsindx = gammalib.GFits(self.m_obs_index+"[OBS_INDEX]"+obs_selection) + if self._use_bkg_scale: + obsindx = gammalib.GFits(self._obs_index+"[OBS_INDEX]"+obs_selection) bkg_scale = obsindx["OBS_INDEX"]["BKG_SCALE"][0] obsindx.close() @@ -567,7 +566,7 @@ class csiactobs(ctools.cscript): run_emax = gammalib.GEnergy(aeff_table.real("HI_THRES"), "TeV") # Append to ebounds - self.m_ebounds.append(run_emin, run_emax) + self._ebounds.append(run_emin, run_emax) else: # Set default values for energy range @@ -578,21 +577,21 @@ class csiactobs(ctools.cscript): aeff_fits.close() # Append instrumental background model - self.models.append(self.iact_background(telescope, obs_id, bkg_scale, bkg_mod_hierarchy[0], run_emin.TeV(), run_emax.TeV())) + self._models.append(self._iact_background(telescope, obs_id, bkg_scale, bkg_mod_hierarchy[0], run_emin.TeV(), run_emax.TeV())) # Logging - if self.logTerse(): - self.log("Adding observation "+str(obs_id)+" (\""+object_name+"\")\n") - - if self.logExplicit(): - self.log(" Event file: "+eventfile+"\n") - self.log(" Effective area: "+aefffile+"\n") - self.log(" Point spread function: "+psffile+"\n") - self.log(" Energy dispersion: "+edispfile+"\n") - self.log(" Background file: "+bkgfile+"\n") - if self.m_use_bkg_scale: - self.log(" Background scale: "+str(bkg_scale)+"\n") - self.log("\n") + if self._logTerse(): + self._log("Adding observation "+str(obs_id)+" (\""+object_name+"\")\n") + + if self._logExplicit(): + self._log(" Event file: "+eventfile+"\n") + self._log(" Effective area: "+aefffile+"\n") + self._log(" Point spread function: "+psffile+"\n") + self._log(" Energy dispersion: "+edispfile+"\n") + self._log(" Background file: "+bkgfile+"\n") + if self._use_bkg_scale: + self._log(" Background scale: "+str(bkg_scale)+"\n") + self._log("\n") # Append observation to XML and set attributes obs = lib.append("observation"); @@ -631,37 +630,37 @@ class csiactobs(ctools.cscript): if lib.size(): # Log header of energy range - if self.logTerse(): - self.log("\n") - self.log.header3("Energy range of obervation list") + if self._logTerse(): + self._log("\n") + self._log.header3("Energy range of obervation list") # Logging if energy range is available - if self.m_ebounds.size() > 0: + if self._ebounds.size() > 0: # Write energy range into log file - self.log(str(self.m_ebounds.emin())+" - "+str(self.m_ebounds.emax())) - self.log("\n") + self._log(str(self._ebounds.emin())+" - "+str(self._ebounds.emax())) + self._log("\n") else: # Write 'not available' into log file - self.log("not available") - self.log("\n") + self._log("not available") + self._log("\n") else: - self.log.header2("WARNING: No observation from given runlist available") + self._log.header2("WARNING: No observation from given runlist available") # Append models provided by 'inmodels' if necessary - if not self.inmodels == None: - if self.logTerse(): - self.log("\n") - self.log.header1("Appending models") + if not self._inmodels == None: + if self._logTerse(): + self._log("\n") + self._log.header1("Appending models") # Loop over input models - for model in self.inmodels: - if self.logTerse(): - self.log.header3("Adding model \""+model.name()+"\"") - self.models.append(model) + for model in self._inmodels: + if self._logTerse(): + self._log.header3("Adding model \""+model.name()+"\"") + self._models.append(model) # Return return @@ -669,10 +668,10 @@ class csiactobs(ctools.cscript): def save(self,outfile,clobber): # Save observation XML file - self.xml.save(outfile) + self._xml.save(outfile) # Save model XML file - self.models.save(self.m_outmodel) + self._models.save(self._outmodel) # Return return @@ -682,14 +681,12 @@ class csiactobs(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates observation file. - """ + # Create instance of application app = csiactobs(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index 60f177d..b90268c 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# Light curve generation script. +# Generates a lightcurve. # # Copyright (C) 2014-2016 Michael Mayer # @@ -23,51 +23,53 @@ import ctools import sys -# ============ # -# csspec class # -# ============ # +# ================ # +# cslightcrv class # +# ================ # 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. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "cslightcrv" - self.version = "1.1.0" + self._name = "cslightcrv" + self._version = "1.1.0" # Initialise some members - self.obs = None + self._outfile = "" - # Initialise some members + # Initialise observation if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - self.obs = argv[0] - argv = argv[1:] + self._obs = argv[0] + argv = argv[1:] else: - self.obs = gammalib.GObservations() - self.obs.clear() - self.m_outfile = "" + self._obs = gammalib.GObservations() # Make sure that parfile exists - file = self.parfile() + self._parfile() # Initialise application if len(argv) == 0: - ctools.cscript.__init__(self, self.name, self.version) + ctools.cscript.__init__(self, self._name, self._version) elif len(argv) ==1: - ctools.cscript.__init__(self, self.name, self.version, *argv) + 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) + self._log_header() + self._log.date(True) # Return return @@ -76,20 +78,18 @@ class cslightcrv(ctools.cscript): """ Destructor. """ - # Write separator into logger - if self.logTerse(): - self.log("\n") - # Return return - def parfile(self): + + # 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" + parfile = self._name+".par" try: pars = gammalib.GApplicationPars(parfile) @@ -136,15 +136,15 @@ class cslightcrv(ctools.cscript): # Return return - def create_tbounds(self): + def _create_tbounds(self): # Create time bin container - self.m_tbins = gammalib.GGti() + self._tbins = gammalib.GGti() # check for binning algorithm if self["tbinalg"].string() == "FILE": try: - self.m_tbins.load(self["tbinfile"].filename()) + self._tbins.load(self["tbinfile"].filename()) except: csv = gammalib.GCsv(self["tbinfile"].filename()) for i in range(csv.nrows()): @@ -152,7 +152,7 @@ class cslightcrv(ctools.cscript): tmax = gammalib.GTime() tmin.mjd(csv.real(i,0)) tmax.mjd(csv.real(i,1)) - self.m_tbins.append(tmin,tmax) + self._tbins.append(tmin,tmax) elif self["tbinalg"].string() == "LIN": # Use linear time binning time_min = self["tmin"].real() @@ -164,165 +164,143 @@ class cslightcrv(ctools.cscript): tmin.mjd(time_min + i*step)# ref) tmax = gammalib.GTime() tmax.mjd(time_min + (i+1)*step)#ref) - self.m_tbins.append(tmin,tmax) + self._tbins.append(tmin,tmax) elif self["tbinalg"].string() == "GTI": # Use GTIs of observations - for obs in self.obs: + for obs in self._obs: for i in range(obs.events().gti().size()): - self.m_tbins.append(obs.events().gti().tstart(i),obs.events().gti().tstop(i)) + 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): + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Set observation if not done before - if self.obs == None or self.obs.size() == 0: + if self._obs == None or self._obs.size() == 0: self._require_inobs("cslightcrv::get_parameters()") - self.obs = self._get_observations() + self._obs = self._get_observations() # Set models if we have none - if self.obs.models().size() == 0: - self.obs.models(self["inmodel"].filename()) + if self._obs.models().size() == 0: + self._obs.models(self["inmodel"].filename()) # Get source name - self.m_srcname = self["srcname"].string() + self._srcname = self["srcname"].string() # Get energy dispersion flag - self.m_edisp = self["edisp"].boolean() + self._edisp = self["edisp"].boolean() # Get time boundaries - self.create_tbounds() + self._create_tbounds() # Unbinned or binned analysis? - self.m_ebins = self["enumbins"].integer() - if self.m_ebins == 0: - self.m_binned = False + self._ebins = self["enumbins"].integer() + if self._ebins == 0: + self._binned = False else: - self.m_binned = True + self._binned = True # Get energy range - self.m_emin = self["emin"].real() - self.m_emax = self["emax"].real() + self._emin = self["emin"].real() + self._emax = self["emax"].real() # Get binning flag - if self.m_binned: - self.m_coordsys = self["coordsys"].string() - self.m_proj = self["proj"].string() - self.m_xref = self["xref"].real() - self.m_yref = self["yref"].real() - self.m_nxpix = self["nxpix"].integer() - self.m_nypix = self["nypix"].integer() - self.m_binsz = self["binsz"].real() + 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.m_outfile = self["outfile"].filename() + self._outfile = self["outfile"].filename() # Get other parameeters - self.m_calc_ulimit = self["calc_ulim"].boolean() - self.m_calc_ts = self["calc_ts"].boolean() - self.m_fix_bkg = self["fix_bkg"].boolean() - self.m_fix_srcs = self["fix_srcs"].boolean() + 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.m_log = False # Logging in client tools - self.m_chatter = self["chatter"].integer() - self.m_clobber = self["clobber"].boolean() - self.m_debug = self["debug"].boolean() + self._chatter = self["chatter"].integer() + self._clobber = self["clobber"].boolean() + self._debug = self["debug"].boolean() - # Return - return - - def models(self, models): - """ - Set model. - """ - # Copy models - self.obs.models(models.clone()) + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Save residual map - self.fits.saveto(self.m_outfile, self.m_clobber) - - # Return - return + # Public methods def run(self): """ Run the script. """ # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + if self._logDebug(): + self._log.cout(True) # Get parameters - self.get_parameters() - - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + 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") + if self._logTerse(): + self._log("\n") + self._log.header1("Observation") + self._log(str(self._obs)) + self._log("\n") # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Adjust model parameters") + if self._logTerse(): + self._log("\n") + self._log.header1("Adjust model parameters") # Adjust model parameters dependent on input user parameters - for model in self.obs.models(): + for model in self._obs.models(): # Set TS flag for all models to false. # Source of interest will be set to true later model.tscalc(False) # Log model name - if self.logExplicit(): - self.log.header3(model.name()) + if self._logExplicit(): + self._log.header3(model.name()) # Deal with the source of interest - if model.name() == self.m_srcname: - if self.m_calc_ts: + if model.name() == self._srcname: + if self._calc_ts: model.tscalc(True) - elif self.m_fix_bkg and not model.classname() == "GModelSky": + elif self._fix_bkg and not model.classname() == "GModelSky": for par in model: - if par.is_free() and self.logExplicit(): - self.log(" Fixing \""+par.name()+"\"\n") + if par.is_free() and self._logExplicit(): + self._log(" Fixing \""+par.name()+"\"\n") par.fix() - elif self.m_fix_srcs and model.classname() == "GModelSky": + elif self._fix_srcs and model.classname() == "GModelSky": for par in model: - if par.is_free() and self.logExplicit(): - self.log(" Fixing \""+par.name()+"\"\n") + if par.is_free() and self._logExplicit(): + self._log(" Fixing \""+par.name()+"\"\n") par.fix() # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate lightcurve") + if self._logTerse(): + self._log("\n") + self._log.header1("Generate lightcurve") # Initialise FITS Table with extension "LIGHTCURVE" - table = gammalib.GFitsBinTable(self.m_tbins.size()) + table = gammalib.GFitsBinTable(self._tbins.size()) table.extname("LIGHTCURVE") # Add Header for compatibility with gammalib.GMWLSpectrum @@ -330,38 +308,38 @@ class cslightcrv(ctools.cscript): table.card("TELESCOP", "CTA", "Name of Telescope") # Create FITS table columns - MJD = gammalib.GFitsTableDoubleCol("MJD", self.m_tbins.size()) + MJD = gammalib.GFitsTableDoubleCol("MJD", self._tbins.size()) MJD.unit("days") - e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self.m_tbins.size()) + e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self._tbins.size()) e_MJD.unit("days") # Create a FITS column for every free parameter columns = [] - for par in self.obs.models()[self.m_srcname]: + for par in self._obs.models()[self._srcname]: if par.is_free(): - col = gammalib.GFitsTableDoubleCol(par.name(), self.m_tbins.size()) + col = gammalib.GFitsTableDoubleCol(par.name(), self._tbins.size()) col.unit(par.unit()) columns.append(col) - e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self.m_tbins.size()) + e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self._tbins.size()) e_col.unit(par.unit()) columns.append(e_col) # Create TS and upper limit columns - TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_tbins.size()) - ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_tbins.size()) + TSvalues = gammalib.GFitsTableDoubleCol("TS", self._tbins.size()) + ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self._tbins.size()) ulim_values.unit("ph/cm2/s") # Loop over energy bins - for i in range(self.m_tbins.size()): + for i in range(self._tbins.size()): # Log information - if self.logTerse(): - self.log("\n") - self.log.header2("Time bin "+str(i)) + if self._logTerse(): + self._log("\n") + self._log.header2("Time bin "+str(i)) # Get time boundaries - tmin = self.m_tbins.tstart(i) - tmax = self.m_tbins.tstop(i) + tmin = self._tbins.tstart(i) + tmax = self._tbins.tstop(i) # Compute time bin center and time width twidth = 0.5 * (tmax - tmin) # in seconds @@ -372,13 +350,13 @@ class cslightcrv(ctools.cscript): e_MJD[i] = twidth / gammalib.sec_in_day # in days # Log information - if self.logExplicit(): - self.log.header3("Selecting events") + if self._logExplicit(): + self._log.header3("Selecting events") # Select events - select = ctools.ctselect(self.obs) - select["emin"] = self.m_emin - select["emax"] = self.m_emax + select = ctools.ctselect(self._obs) + select["emin"] = self._emin + select["emax"] = self._emax select["tmin"] = tmin.convert(select._time_reference()) select["tmax"] = tmax.convert(select._time_reference()) select["rad"] = "UNDEFINED" @@ -390,113 +368,113 @@ class cslightcrv(ctools.cscript): obs = select.obs() # Binned analysis - if self.m_binned: + if self._binned: # Header - if self.logTerse(): - self.log.header3("Binning events") + if self._logTerse(): + self._log.header3("Binning events") # Bin events bin = ctools.ctbin(select.obs()) bin["usepnt"] = False bin["ebinalg"] = "LOG" - bin["xref"] = self.m_xref - bin["yref"] = self.m_yref - bin["binsz"] = self.m_binsz - bin["nxpix"] = self.m_nxpix - bin["nypix"] = self.m_nypix - bin["enumbins"] = self.m_ebins - bin["emin"] = self.m_emin - bin["emax"] = self.m_emax - bin["coordsys"] = self.m_coordsys - bin["proj"] = self.m_proj + 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") + 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.m_xref - expcube["yref"] = self.m_yref - expcube["binsz"] = self.m_binsz - expcube["nxpix"] = self.m_nxpix - expcube["nypix"] = self.m_nypix - expcube["enumbins"] = self.m_ebins - expcube["emin"] = self.m_emin - expcube["emax"] = self.m_emax - expcube["coordsys"] = self.m_coordsys - expcube["proj"] = self.m_proj + 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") + 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.m_xref - psfcube["yref"] = self.m_yref - psfcube["binsz"] = self.m_binsz - psfcube["nxpix"] = self.m_nxpix - psfcube["nypix"] = self.m_nypix - psfcube["enumbins"] = self.m_ebins - psfcube["emin"] = self.m_emin - psfcube["emax"] = self.m_emax - psfcube["coordsys"] = self.m_coordsys - psfcube["proj"] = self.m_proj + 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.m_edisp: + if self._edisp: # Header - if self.logTerse(): - self.log.header3("Creating energy dispersion cube") + 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.m_xref - edispcube["yref"] = self.m_yref - edispcube["binsz"] = self.m_binsz - edispcube["nxpix"] = self.m_nxpix - edispcube["nypix"] = self.m_nypix - edispcube["enumbins"] = self.m_ebins - edispcube["emin"] = self.m_emin - edispcube["emax"] = self.m_emax - edispcube["coordsys"] = self.m_coordsys - edispcube["proj"] = self.m_proj + 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") + 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.m_xref - bkgcube["yref"] = self.m_yref - bkgcube["binsz"] = self.m_binsz - bkgcube["nxpix"] = self.m_nxpix - bkgcube["nypix"] = self.m_nypix - bkgcube["enumbins"] = self.m_ebins - bkgcube["emin"] = self.m_emin - bkgcube["emax"] = self.m_emax - bkgcube["coordsys"] = self.m_coordsys - bkgcube["proj"] = self.m_proj + 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 @@ -506,7 +484,7 @@ class cslightcrv(ctools.cscript): models = bkgcube.models() # Set precomputed binned response - if self.m_edisp: + if self._edisp: obs[0].response(expcube.expcube(), psfcube.psfcube(), edispcube.edispcube(), bkgcube.bkgcube()) else: @@ -514,7 +492,7 @@ class cslightcrv(ctools.cscript): bkgcube.bkgcube()) # Fix background models if required - if self.m_fix_bkg: + if self._fix_bkg: for model in models: if not model.classname() == "GModelSky": for par in model: @@ -524,20 +502,20 @@ class cslightcrv(ctools.cscript): obs.models(models) # Header - if self.logTerse(): - self.log.header3("Performing fit") + if self._logTerse(): + self._log.header3("Performing fit") # Likelihood like = ctools.ctlike(obs) - like["edisp"] = self.m_edisp + like["edisp"] = self._edisp like.run() # Skip bin if no event was present if like.obs().logL() == 0.0: # Log information - if self.logTerse(): - self.log("No event in this time bin. Bin is skipped\n") + if self._logTerse(): + self._log("No event in this time bin. Bin is skipped\n") # Set all values to 0 for col in columns: @@ -548,19 +526,19 @@ class cslightcrv(ctools.cscript): # Get results fitted_models = like.obs().models() - source = fitted_models[self.m_srcname] + source = fitted_models[self._srcname] # Calculate Upper Limit ulimit_value = -1.0 - if self.m_calc_ulimit: + if self._calc_ulimit: # Logging information - if self.logTerse(): - self.log.header3("Computing upper limit") + if self._logTerse(): + self._log.header3("Computing upper limit") # Create upper limit object ulimit = ctools.ctulimit(like.obs()) - ulimit["srcname"] = self.m_srcname + ulimit["srcname"] = self._srcname ulimit["eref"] = 1.0 # Try to run upper limit and catch exceptions @@ -568,13 +546,13 @@ class cslightcrv(ctools.cscript): ulimit.run() ulimit_value = ulimit.flux_ulimit() except: - if self.logTerse(): - self.log("Upper limit calculation failed\n") + if self._logTerse(): + self._log("Upper limit calculation failed\n") ulimit_value = -1.0 # Get TS value TS = -1.0 - if self.m_calc_ts: + if self._calc_ts: TS = source.ts() # Set values for storage @@ -592,20 +570,20 @@ class cslightcrv(ctools.cscript): ulim_values[i] = ulimit_value # Log information - if self.logExplicit(): - self.log.header3("Results of bin "+str(i)+": MJD "+str(tmin.mjd())+"-"+str(tmax.mjd())) + 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.m_calc_ts and TSvalues[i] > 0.0: - self.log(" > TS = "+str(TS)+" \n") - if self.m_calc_ulimit and ulim_values[i] > 0.0: - self.log(" > UL = "+str(ulim_values[i])+" [ph/cm2/s]") - self.log("\n") + 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") # Append filles columns to fits table table.append(MJD) @@ -616,8 +594,69 @@ class cslightcrv(ctools.cscript): table.append(ulim_values) # Create the FITS file now - self.fits = gammalib.GFits() - self.fits.append(table) + self._fits = gammalib.GFits() + self._fits.append(table) + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Read ahead output parameters + self._read_ahead(True) + + # Run the script + self.run() + + # Save lightcurve + self.save() + + # Return + return + + def save(self): + """ + Save lightcurve. + """ + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save lightcurve") + + # Get outmap parameter + self._outfile = self["outfile"].filename() + + # Continue only filename and residual map are valid + if self._outfile != "NONE" and self._fits != None: + + # Log file name + if self._logTerse(): + self._log("Save lightcurve into \""+self._outfile+"\".\n") + + # Save spectrum + self._fits.saveto(self._outfile, self._clobber) + + # Return + return + + def lightcurve(self): + """ + Return lightcurve FITS file. + + Returns: + FITS file containing lightcurve. + """ + # Return + return self._fits + + def models(self, models): + """ + Set model. + """ + # Copy models + self._obs.models(models.clone()) # Return return @@ -627,14 +666,12 @@ class cslightcrv(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates light curve. - """ + # Create instance of application app = cslightcrv(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 0b82fd4..327e029 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -59,7 +59,6 @@ class csspec(ctools.cscript): self._chatter = 2 self._clobber = True self._debug = False - self._log_clients = False # Initialise observation if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): @@ -206,6 +205,8 @@ class csspec(ctools.cscript): # Return return + + # Public methods def run(self): """ Run the script.