From 25a6a46b78a0f2a94fe8860a7160881863551302 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 12 Mar 2016 10:56:03 +0100 Subject: [PATCH] Adapt to new GApplication interface with private methods and add par file. --- cscripts/Makefile.am | 5 ++++- cscripts/cscaldb.par | 6 +++--- cscripts/csmodelinfo.par | 45 +++++++++++++++++++++++++++++++++++++++++++++ cscripts/csmodelinfo.py | 611 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- cscripts/csmodelmerge.par | 30 ++++++++++++++++++++++++++++++ cscripts/csmodelmerge.py | 196 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------- cscripts/csobs2caldb.par | 32 ++++++++++++++++++++++++++++++++ cscripts/csobs2caldb.py | 325 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------- cscripts/csobsdef.par | 41 +++++++++++++++++++++++++++++++++++++++++ cscripts/csobsdef.py | 262 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------------------------------------------- cscripts/csobsinfo.par | 37 +++++++++++++++++++++++++++++++++++++ cscripts/csobsinfo.py | 465 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- cscripts/cspull.par | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ cscripts/cspull.py | 544 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ cscripts/csresmap.par | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cscripts/csresmap.py | 480 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ cscripts/cssens.par | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ cscripts/cssens.py | 808 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- cscripts/csspec.par | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++++ cscripts/csspec.py | 487 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- cscripts/cstsdist.par | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cscripts/cstsdist.py | 500 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- cscripts/cstsmapmerge.par | 38 ++++++++++++++++++++++++++++++++++++++ cscripts/cstsmapmerge.py | 401 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------ test/test_cscripts.sh.in | 8 ++++---- 25 files changed, 3055 insertions(+), 2542 deletions(-) create mode 100644 cscripts/csmodelinfo.par create mode 100644 cscripts/csmodelmerge.par create mode 100644 cscripts/csobs2caldb.par create mode 100644 cscripts/csobsdef.par create mode 100644 cscripts/csobsinfo.par create mode 100644 cscripts/cspull.par create mode 100644 cscripts/csresmap.par create mode 100644 cscripts/cssens.par create mode 100644 cscripts/csspec.par create mode 100644 cscripts/cstsdist.par create mode 100644 cscripts/cstsmapmerge.par diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index ccd898b..6826fc9 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -10,7 +10,10 @@ pythonscripts = __init__.py \ obsutils.py # Parfiles to be distributed in $(prefix)/syspfiles -parfiles = cscaldb.par +parfiles = cspull.par cssens.par cstsdist.par cscaldb.par csresmap.par \ + csobsdef.par csspec.par csobsinfo.par \ + csmodelinfo.par csmodelmerge.par cstsmapmerge.par \ + csobs2caldb.par # Scripts for which symbolic links will be created in bindir linkscripts = cspull cssens cstsdist cscaldb csresmap \ diff --git a/cscripts/cscaldb.par b/cscripts/cscaldb.par index efc4a71..4e9ebb5 100644 --- a/cscripts/cscaldb.par +++ b/cscripts/cscaldb.par @@ -7,7 +7,7 @@ # File: cscaldb.par # # Version: 1.1.0 # # # -# Author: Jurgen Knodlseder # +# Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # # IRAP # # # @@ -20,5 +20,5 @@ 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, "cscaldb.log",,, "Log filename" +mode, s, h, ql,,, "Mode of automatic parameters" +logfile, f, h, cscaldb.log,,, "Log filename" diff --git a/cscripts/csmodelinfo.par b/cscripts/csmodelinfo.par new file mode 100644 index 0000000..a2516e5 --- /dev/null +++ b/cscripts/csmodelinfo.par @@ -0,0 +1,45 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csmodelinfo.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csmodelinfo script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inmodel, f, a, model.xml,,, "Input model XML file" +ds9file, f, h, NONE,,, "Output DS9 region file containing models" + +# +# Script parameters +#================== +pnt_type, s, h, cross,circle|box|diamond|cross|x|arrow|boxcircle,, "Marker type for point sources" +pnt_mark_size, i, h, 12,,, "Marker size for point sources" +show_labels, b, h, yes,yes|no,, "Add source labels?" +width, i, h, 2,,, "Line width for regions" +fontfamily, s, h, helvetica,helvetica|times|courier,, "Font family for source labels" +fontsize, i, h, 12,,, "Font size for source labels" +fontweight, s, h, normal,normal|bold,, "Font weight for source labels" +fontslant, s, h, roman,roman|italic,, "Font slant for source labels" +show_ext_type, b, h, yes,,, "Show type of extended model in source name?" +free_color, s, h, green,,, "Color for sources with free parameters (any DS9 color or hex code)" +fixed_color, s, h, magenta,,, "Color for source without free parameters (any DS9 color or hex code)" + +# +# 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, csmodelinfo.log,,, "Log filename" diff --git a/cscripts/csmodelinfo.py b/cscripts/csmodelinfo.py index 0e61b6b..6f6fc1e 100755 --- a/cscripts/csmodelinfo.py +++ b/cscripts/csmodelinfo.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # ========================================================================== -# Dump information about models into log file +# Shows the content of a model container # # Copyright (C) 2015-2016 Michael Mayer # @@ -28,34 +28,50 @@ import sys # ================= # class csmodelinfo(ctools.cscript): """ - This class dumps information about a model container into a logfile or - on screen. This might be helpful for quick access to the model container, - showing, e.g., number of sources, background models, free and fixed - parameters etc. The tool can also dump the models into a region file to - display e.g. with ds9. + Shows the content of a model container. + + The script is useful to show, e.g., the number of sources in a model + container, the background models, free and fixed parameters etc. + The script can also dump the models into a region file to for display + with DS9. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "csmodelinfo" - self.version = "1.1.0" + self._name = "csmodelinfo" + self._version = "1.1.0" - # Make sure that parfile exists - file = self.parfile() + # Initialise class members + self._models = gammalib.GModels() + self._ds9file = "NONE" + self._pnt_type = "" + self._pnt_mark_size = 12 + self._show_labels = True + self._width = 2 + self._fontfamily = "helvetica" + self._fontsize = 12 + self._fontweight = "normal" + self._fontslant = "roman" + 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) + 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 @@ -67,92 +83,197 @@ class csmodelinfo(ctools.cscript): # 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 that parfile was not found - print("Parfile \""+parfile+"\" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inmodel","f","a","model.xml","","","Input model XML file")) - pars.append(gammalib.GApplicationPar("ds9file","f","h","NONE","","","DS9 region file containing models")) - pars.append(gammalib.GApplicationPar("pnt_type","s","h","cross","circle|box|diamond|cross|x|arrow|boxcircle","","Marker type for point sources")) - pars.append(gammalib.GApplicationPar("pnt_mark_size","i","h","12","","","Marker size for point sources")) - pars.append(gammalib.GApplicationPar("show_labels","b","h","yes","yes|no","","Add source labels?")) - pars.append(gammalib.GApplicationPar("width","i","h","2","","","Line width for regions")) - pars.append(gammalib.GApplicationPar("font","s","h","helvetica","helvetica|times|courier","","Font for source labels")) - pars.append(gammalib.GApplicationPar("fontsize","i","h","12","","","Font size for source labels")) - pars.append(gammalib.GApplicationPar("fonttype","s","h","normal","normal|bold","","Use normal or bold font?")) - pars.append(gammalib.GApplicationPar("fonttype2","s","h","roman","roman|italic","","Use roman or italic font?")) - pars.append(gammalib.GApplicationPar("show_ext_type","b","h","yes","yes|no","","Show type of extended model in source name")) - pars.append(gammalib.GApplicationPar("free_color","s","h","green","","","Color for sources with free parameters (any ds9 color or hex code)")) - pars.append(gammalib.GApplicationPar("fixed_color","s","h","magenta","","","Color for source without free parameters (any ds9 color or hex code)")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csmodelinfo.log","","","Log filename")) - pars.save(parfile) - # Return - return - - def get_parameters(self): + # Private methods + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get models - self.models = gammalib.GModels(self["inmodel"].filename()) - - # Get ds9 filename - self.ds9file = self["ds9file"].filename() + self._models = gammalib.GModels(self["inmodel"].filename()) # Get hidden parameters for region file - self.pnt_type = self["pnt_type"].string() - self.pnt_mark_size = self["pnt_mark_size"].integer() - self.show_labels = self["show_labels"].boolean() - self.width = self["width"].integer() - self.font = self["font"].string() - self.fontsize = self["fontsize"].integer() - self.fonttype = self["fonttype"].string() - self.fonttype2 = self["fonttype2"].string() - self.show_ext_type = self["show_ext_type"].boolean() - self.free_color = self["free_color"].string() - self.fixed_color = self["fixed_color"].string() + self._pnt_type = self["pnt_type"].string() + self._pnt_mark_size = self["pnt_mark_size"].integer() + self._show_labels = self["show_labels"].boolean() + self._width = self["width"].integer() + self._fontfamily = self["fontfamily"].string() + self._fontsize = self["fontsize"].integer() + self._fontweight = self["fontweight"].string() + self._fontslant = self["fontslant"].string() + self._show_ext_type = self["show_ext_type"].boolean() + self._free_color = self["free_color"].string() + self._fixed_color = self["fixed_color"].string() + + # Read ahead DS9 filename + if self._read_ahead: + self._ds9file = self["ds9file"].filename() + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return + + def _model2ds9string(self, model): + """ + Converts model into a DS9 region string. + + Args: + model: Model. + + Returns: + DS9 region string. + """ + # Initialise DS9 region string + ds9string = "" + + # Determine region color. The color will change between a model + # where all parameters are fixed and a model with at least one + # free parameter. + color = self._fixed_color + for par in model: + if par.is_free(): + color = self._free_color + break + + # Retrieve model sky direction. The model is skipped in case that + # it does not provide a sky direction. + is_valid = True + try: + modelpos = model.spatial().dir() + except AttributeError: + self._log("Skip source \""+model.name()+"\" since no " + "sky direction is associated with that model.\n") + is_valid = False + + # Continue only if sky direction was found + if is_valid: + + # Retrieve model type and name + modeltype = model.type() + modelname = model.name() + + # Handle point source + if modeltype == "PointSource": + + # Append point to DS9 string + ds9string += "point(" + ds9string += str(modelpos.ra_deg()) + ds9string += "," + ds9string += str(modelpos.dec_deg()) + ds9string += ") # " + ds9string += "point="+self._pnt_type+" " + ds9string += str(self._pnt_mark_size)+" " + + # Handle extended sources + elif modeltype == "ExtendedSource": + + # Retrieve spatial model + spatial = model.spatial() + classname = spatial.classname() + + # Handle radial sources + if "Radial" in classname: + + # Retrieve short name of model class (e.g. "Gauss" + # or "Disk) + shorttype = classname.split("Radial")[-1] + + # Handle Disk and Shell model + if (classname == "GModelSpatialRadialDisk" or + classname == "GModelSpatialRadialShell"): + size = spatial.radius() + + # Handle Gauss Model + elif classname == "GModelSpatialRadialGauss": + size = spatial.sigma() + + # Skip if source is unknown + else: + self._log("Skip source \""+model.name()+"\"" + " since the radial model \""+classname+ + "\" is unkown.\n") + is_valid = False + + # Append circle to DS9 string + if is_valid: + ds9string += "circle(" + ds9string += str(modelpos.ra_deg()) + ds9string += "," + ds9string += str(modelpos.dec_deg()) + ds9string += "," + ds9string += str(size*3600.)+"\") # " + + # Handle elliptical sources + elif "Elliptical" in classname: + + # Retrieve short name and source size + shorttype = classname.split("Elliptical")[-1] + size1 = spatial.semimajor() + size2 = spatial.semiminor() + angle = spatial.posangle() + + # Append ellipse to DS9 string + ds9string += "ellipse(" + ds9string += str(modelpos.ra_deg()) + ds9string += "," + ds9string += str(modelpos.dec_deg()) + ds9string += "," + ds9string += str(size1*3600.)+"\"" + ds9string += "," + ds9string += str(size2*3600.)+"\"" + ds9string += "," + ds9string += str(angle+90.0)+") # " + + # Skip if source is neither radial nor elliptical + else: + self._log("Skip source \""+model.name()+"\" " + "since the model \""+classname+"\" is " + "neither a point source, a radial source, " + "nor a elliptical source.\n") + is_valid = False + + # Add short model type to modelname + if self._show_ext_type: + modelname +=" ("+shorttype+")" + + # Add DS9 attributes + if is_valid: + ds9string += "color="+color+" " + ds9string += "width="+str(self._width)+" " + ds9string += "font=\""+self._fontfamily+" " + ds9string += str(self._fontsize)+" " + ds9string += self._fontweight+" " + ds9string += self._fontslant+"\"" + if self._show_labels: + ds9string += " text={"+modelname+"}" + + # Return string + return ds9string + + # 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() + self._get_parameters() - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") - - if self.logVerbose(): - self.log("\n") - self.log.header1("Models in detail") - self.log("\n") - self.log(str(self.models)) - self.log("\n") + # Optionally log models + if self._logVerbose(): + self._log("\n") + self._log.header1("Models") + self._log("\n") + self._log(str(self._models)) + self._log("\n") # Initialise output values types = {} @@ -171,13 +292,13 @@ class csmodelinfo(ctools.cscript): n_par_free_temp = 0 # Loop over models and initialise dictionaries with model names - for model in self.models: + for model in self._models: types[model.type()] = 0 instruments[model.instruments()] = 0 pars_at_limit[model.name()] = [] # Loop over models again and fill values to be dumped - for model in self.models: + for model in self._models: # Model types types[model.type()] += 1 @@ -218,7 +339,8 @@ class csmodelinfo(ctools.cscript): else: n_par_fixed += 1 - # If sky model, divide free parameters into spatial, spectral and temporal + # If sky model, divide free parameters into spatial, spectral + # and temporal if skymodel: for par in model.spectral(): if par.is_free(): @@ -229,90 +351,93 @@ class csmodelinfo(ctools.cscript): for par in model.temporal(): if par.is_free(): n_par_free_temp += 1 + + # Optionally log model information + if self._logTerse(): - # Log summary - self.log("\n") - self.log.header1("Summary") + # Log summary + self._log("\n") + self._log.header1("Summary") - # Log instruments - self.log.header3("Instrument specific models") - for inst, n_inst in instruments.items(): - if inst == "": - inst = "All" - self.log.parformat(inst) - self.log(str(n_inst)) - self.log("\n") + # Log instruments + self._log.header3("Instrument specific models") + for inst, n_inst in instruments.items(): + if inst == "": + inst = "All" + self._log.parformat(inst) + self._log(str(n_inst)) + self._log("\n") - # Log model types - self.log.header3("Model types") - for modeltype, n_types in types.items(): - self.log.parformat(modeltype) - self.log(str(n_types)) - self.log("\n") + # Log model types + self._log.header3("Model types") + for modeltype, n_types in types.items(): + self._log.parformat(modeltype) + self._log(str(n_types)) + self._log("\n") - # Log parameter information - self.log("\n") - self.log.header1("Parameter information") - self.log.parformat("All parameters") - self.log(str(n_par_total)) - self.log("\n") - self.log.parformat("Fixed parameters") - self.log(str(n_par_fixed)) - self.log("\n") - self.log.parformat("Free parameters (total)") - self.log(str(n_par_free)) - self.log("\n") - self.log.parformat("Free background parameters") - self.log(str(n_par_free_bkg)) - self.log("\n") - self.log.parformat("Free source parameters") - self.log(str(n_par_free_src)) - self.log("\n") - if n_par_free_spec > 0: - self.log.parformat("Free spectral parameters") - self.log(str(n_par_free_spec)) - self.log("\n") - if n_par_free_spat > 0: - self.log.parformat("Free spatial parameters") - self.log(str(n_par_free_spat)) - self.log("\n") - if n_par_free_temp > 0: - self.log.parformat("Free temporal parameters") - self.log(str(n_par_free_temp)) - self.log("\n") - self.log.parformat("Parameters at limit") - self.log(str(n_par_at_limit)) - self.log("\n") + # Log model parameter information + self._log("\n") + self._log.header1("Parameter information") + self._log.parformat("All parameters") + self._log(str(n_par_total)) + self._log("\n") + self._log.parformat("Fixed parameters") + self._log(str(n_par_fixed)) + self._log("\n") + self._log.parformat("Free parameters (total)") + self._log(str(n_par_free)) + self._log("\n") + self._log.parformat("Free background parameters") + self._log(str(n_par_free_bkg)) + self._log("\n") + self._log.parformat("Free source parameters") + self._log(str(n_par_free_src)) + self._log("\n") + if n_par_free_spec > 0: + self._log.parformat("Free spectral parameters") + self._log(str(n_par_free_spec)) + self._log("\n") + if n_par_free_spat > 0: + self._log.parformat("Free spatial parameters") + self._log(str(n_par_free_spat)) + self._log("\n") + if n_par_free_temp > 0: + self._log.parformat("Free temporal parameters") + self._log(str(n_par_free_temp)) + self._log("\n") + self._log.parformat("Parameters at limit") + self._log(str(n_par_at_limit)) + self._log("\n") - # Log parameters at limit (if any) - if n_par_at_limit > 0: - self.log.header3("Parameters at limit") - for source, parameter in pars_at_limit.items(): - if len(parameter): - for par in parameter: - self.log.parformat(source) - self.log(par) - self.log("\n") + # Log parameters at limit (if any) + if n_par_at_limit > 0: + self._log.header3("Parameters at limit") + for source, parameter in pars_at_limit.items(): + if len(parameter): + for par in parameter: + self._log.parformat(source) + self._log(par) + self._log("\n") - # Log free parameters if requested - if self.logExplicit() and len(free_src_pars): - self.log("\n") - self.log.header2("Free source parameters") - for source, parameter in free_src_pars.items(): - if len(parameter): - self.log.header3(source) - for par in parameter: - self.log(str(par)) - self.log("\n") - self.log("\n") + # Optionally log free parameters + if self._logExplicit() and len(free_src_pars): + self._log("\n") + self._log.header2("Free source parameters") + for source, parameter in free_src_pars.items(): + if len(parameter): + self._log.header3(source) + for par in parameter: + self._log(str(par)) + self._log("\n") + self._log("\n") - # Log TS values if available - if len(ts): - self.log.header3("Test statistics") - for source,tsvalue in ts.items(): - self.log.parformat(source) - self.log("%.2f"%tsvalue) - self.log("\n") + # Log TS values if available + if len(ts): + self._log.header3("Test statistics") + for source,tsvalue in ts.items(): + self._log.parformat(source) + self._log("%.2f"%tsvalue) + self._log("\n") # Return return @@ -321,37 +446,49 @@ class csmodelinfo(ctools.cscript): """ Save models to ds9 region file if required """ - - # Check if ds9 file is valid - if not self.ds9file == "NONE": + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save models in DS9 file") + + # Get output filename in case it was not read ahead + self._ds9file = self["ds9file"].filename() + + # Check if DS9 file is valid + if self._ds9file != "NONE": - self.log("\n") - self.log.header1("Save region file") + # Log filename + if self._logTerse(): + self._log(gammalib.parformat("DS9 filename")) + self._log(self._ds9file.url()) + self._log("\n") # Open file - f = open(self.ds9file.url(),"w") + f = open(self._ds9file.url(),"w") # Write coordinate system f.write("fk5\n") # Loop over models - for model in self.models: + for model in self._models: # Continue only if point source or extended source model - if model.type() == "PointSource" or model.type() == "ExtendedSource": - - line = self.model2ds9string(model) + if (model.type() == "PointSource" or + model.type() == "ExtendedSource"): + line = self._model2ds9string(model) if len(line): f.write(line+"\n") # Logging for diffuse components elif model.type() == "DiffuseSource": - self.log("Skipping diffuse model "+model.name()+"\n") + self._log("Skipping diffuse model \""+model.name()+ + "\"\n") # Logging for background components else: - if self.logExplicit(): - self.log("Skipping background model "+model.name()+"\n") + if self._logExplicit(): + self._log("Skipping background model \""+ + model.name()+"\"\n") # Close file f.close() @@ -359,127 +496,13 @@ class csmodelinfo(ctools.cscript): # Return return - def model2ds9string(self, model): - """ - Converts model instance to ds9 region string - """ - - # Initialise return string - ds9string = "" - - # Determine color - color = self.fixed_color - for par in model: - if par.is_free(): - color = self.free_color - break - - # Initialise flag ot continue - is_valid = True - - # Retrieve model position - try: - modelpos = model.spatial().dir() - except AttributeError: - self.log("Skipping source "+model.name()+": Could not determine sky position\n") - is_valid = False - - # Continue only if position was found - if is_valid: - - # Retrieve model type and name - modeltype = model.type() - modelname = model.name() - - # Handle point source - if modeltype == "PointSource": - - ds9string += "point(" - ds9string += str(modelpos.ra_deg()) - ds9string += "," - ds9string += str(modelpos.dec_deg()) - ds9string += ") # " - ds9string += "point="+self.pnt_type+" " - ds9string += str(self.pnt_mark_size)+" " - - # Hanlde extended sources - elif modeltype == "ExtendedSource": - - # Retrieve spatial model - spatial = model.spatial() - classname = spatial.classname() - - # Handle radial sources - if "Radial" in classname: - - # Retrieve short name of model class (e.g. "Gauss" or "Disk) - shorttype = classname.split("Radial")[-1] - - # Handle Disk and Shell model - if classname == "GModelSpatialRadialDisk" or classname == "GModelSpatialRadialShell": - size = spatial.radius() - - # Handle Gauss Model - elif classname == "GModelSpatialRadialGauss": - size = spatial.sigma() - - # Skip if source is unknown - else: - self.log("Skipping source "+model.name()+": Spatial model \""+classname+"\" unkown\n") - is_valid = False - - # Append radial ds9 string - if is_valid: - ds9string += "circle(" - ds9string += str(modelpos.ra_deg()) - ds9string += "," - ds9string += str(modelpos.dec_deg()) - ds9string += "," - ds9string += str(size*3600.)+"\") #" - - # Handle elliptical sources - elif "Elliptical" in classname: - - # Retrieve short name and source size - shorttype = classname.split("Elliptical")[-1] - size1 = spatial.semimajor() - size2 = spatial.semiminor() - angle = spatial.posangle() - - # Append elliptical ds9 string - ds9string += "ellipse(" - ds9string += str(modelpos.ra_deg()) - ds9string += "," - ds9string += str(modelpos.dec_deg()) - ds9string += "," - ds9string += str(size1*3600.)+"\"" - ds9string += "," - ds9string += str(size2*3600.)+"\"" - ds9string += "," - ds9string += str(angle+90.0)+") # " - - # Skip if source is neither radial nor elliptical - else: - self.log("Skipping source "+model.name()+": Spatial model neither \"Radial\" nor \"Elliptical\"\n") - - # Add short name to modelname - if self.show_ext_type: - modelname +=" ("+shorttype+")" - - # Add user defined attributes - ds9string += "color="+color+" " - ds9string += "width="+str(self.width)+" " - ds9string += "font=\""+self.font+" "+str(self.fontsize)+" "+self.fonttype+" "+self.fonttype2+"\" " - if self.show_labels: - ds9string += "text={"+modelname+"}" - - # Return string - return ds9string - def execute(self): """ Execute the script. """ + # Set read ahead flag + self._read_ahead = True + # Run the script self.run() @@ -494,14 +517,12 @@ class csmodelinfo(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates model information output. - """ + # Create instance of application app = csmodelinfo(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csmodelmerge.par b/cscripts/csmodelmerge.par new file mode 100644 index 0000000..8f36965 --- /dev/null +++ b/cscripts/csmodelmerge.par @@ -0,0 +1,30 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csmodelmerge.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csmodelmerge script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inmodels, s, a, model.xml,,, "Input model definition XML files" +outmodel, f, a, NONE,,, "Output model definition XML 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, csmodelmerge.log,,, "Log filename" diff --git a/cscripts/csmodelmerge.py b/cscripts/csmodelmerge.py index 1a15a72..cfc20a0 100755 --- a/cscripts/csmodelmerge.py +++ b/cscripts/csmodelmerge.py @@ -1,8 +1,8 @@ #!/usr/bin/env python # ========================================================================== -# Merges model XML files into one file +# Merge model definition XML files # -# 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 @@ -29,31 +29,38 @@ import glob # ================== # class csmodelmerge(ctools.cscript): """ - This class merges an arbitray number of model container into one output - model XML file. + Merge model definition XML files. + + An arbitrary number of model definition XML files will be merged into + a single model definition XML file. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "csmodelmerge" - self.version = "1.1.0" + self._name = "csmodelmerge" + self._version = "1.1.0" - # Make sure that parfile exists - file = self.parfile() + # 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) + 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,124 +72,122 @@ class csmodelmerge(ctools.cscript): # 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 that parfile was not found - print("Parfile \""+parfile+"\" not found. Create default parfile.") - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inmodels","s","a","model.xml","","","Input model XML files")) - pars.append(gammalib.GApplicationPar("outmodel","f","a","NONE","","","Output model file")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csmodelmerge.log","","","Log filename")) - pars.save(parfile) - - # Return - return - - def get_parameters(self): + # Private methods + def _get_parameters(self): """ - Get parameters from parfile and setup the observation. + Get parameters from parfile. """ - # Get input models - self.inmodels = self["inmodels"].string() + # Get input models string + inmodels = self["inmodels"].string() - # Interpet inmodels string # Handle ascii files - if "@" == self.inmodels[0]: - self.files = open(self.inmodels.replace("@","")).read().splitlines() + if "@" == inmodels[0]: + self._files = open(inmodels.replace("@","")).read().splitlines() # Handle wild card strings - elif "*" in self.inmodels: - self.files = glob.glob(self.inmodels) + elif "*" in inmodels: + self._files = glob.glob(inmodels) - # Handle space-separated list - elif " " in self.inmodels: - self.files = self.inmodels.split(" ") + # Handle space separated list + elif " " in inmodels: + self._files = inmodels.split(" ") - # Handle comma-separated list - elif "," in self.inmodels: - self.files = self.inmodels.split(",") + # Handle semi-colon separated list + elif ";" in inmodels: + self._files = inmodels.split(";") # Throw exception if input models cannot be decoded else: - msg = "Parameter \"inmodels\" must contain either an @ASCII"+\ - " file, a comma-separated or whitespace-separated list"+\ - " of files or a wildcard string" + msg = "Parameter \"inmodels\" must contain either an @ASCII "\ + "file, a semi-colon-separated or whitespace-separated "\ + "list of files or a wildcard string." raise RuntimeError(msg) - # Get output file - self.outmodel = self["outmodel"].filename() + # Read ahead output filename + if self._read_ahead: + self.outmodel = self["outmodel"].filename() + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return + + # Public methods def run(self): """ Run the script. """ # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + 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") - - # Dump header - if self.logTerse(): - self.log("\n") - self.log.header1("Merging Models") - self.log("\n") - - # Initialise output model container - self.models = gammalib.GModels() + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Merge models") + + # Initialise model container + self._models = gammalib.GModels() # Loop over model files - for modelfile in self.files: - - # Logging - self.log("Adding "+modelfile) - self.log("\n") + for file in self._files: + + # Construct container from XML file + models = gammalib.GModels(file) - # Open models - model_add = gammalib.GModels(modelfile) + # Log number of models to add + if self._logTerse(): + nmodels = models.size() + if nmodels == 0: + self._log(gammalib.parformat("Add no model from file")) + elif nmodels == 1: + self._log(gammalib.parformat("Add 1 model from file")) + else: + self._log(gammalib.parformat("Add %d models from file" % + nmodels)) + self._log(file) + self._log("\n") - # Append models to container - for model in model_add: - self.models.append(model) + # Extend model container by adding all models in the model file + self._models.extend(models) + # Log total number of models + if self._logTerse(): + self._log(gammalib.parformat("Models after merging")) + self._log(self._models.size()) + self._log("\n") + # Return return def save(self): """ - Save pointings to ds9 region file if required + Save model definition XML file """ - # Logging - if self.logExplicit(): - self.log("\n") - self.log.header1("Saving model file") - self.log("\n") + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save models") + + # Get output filename in case it was not read ahead + self._outmodel = self["outmodel"].filename() + # Log filename + if self._logTerse(): + self._log(gammalib.parformat("Model definition XML file")) + self._log(self._outmodel.url()) + self._log("\n") + # Save models - self.models.save(self.outmodel) + self._models.save(self._outmodel) # Return return @@ -191,6 +196,9 @@ class csmodelmerge(ctools.cscript): """ Execute the script. """ + # Set read ahead flag + self._read_ahead = True + # Run the script self.run() @@ -205,14 +213,12 @@ class csmodelmerge(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Merges model XML files - """ + # Create instance of application app = csmodelmerge(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csobs2caldb.par b/cscripts/csobs2caldb.par new file mode 100644 index 0000000..9f8c8a0 --- /dev/null +++ b/cscripts/csobs2caldb.par @@ -0,0 +1,32 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csobs2caldb.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csobs2caldb script # +# # +############################################################################# +# +# General parameters +#=================== +inobs, f, a, obs.xml,,, "Observation definition XML file" +rspname, s, a, NONE,,, "Response output name (e.g. Zenith50)" +outirfs, f, h, irf_file.fits,,, "Output IRF file name" +index, i, h, 0,,, "Index of observation to be used from XML 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, csobs2caldb.log,,, "Log filename" diff --git a/cscripts/csobs2caldb.py b/cscripts/csobs2caldb.py index 62ef217..1dd5076 100755 --- a/cscripts/csobs2caldb.py +++ b/cscripts/csobs2caldb.py @@ -30,34 +30,46 @@ import datetime # ================= # class csobs2caldb(ctools.cscript): """ - This class implements the creation of a caldb entry for a particular - observation, which might be helpful for running simulations for current - IACTs. It takes an observation definition file as input and uses the - first observation to create a caldb entry which can be used for - simulations and sensitivity studies later on. + Creates a calibration database entry for an observation. + + The creation of a calibration database entry is useful for performing + simulations for current Imaging Air Cherenkov Telescopes (IACTs). + The class takes an observation definition XML file as input and uses + one observation to create a calibration database entry. By default + the first observation will be used, but the user can specify the + index of any observation using the hidden "index" parameter. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ - # Set name - self.name = "csobs2caldb" - self.version = "1.1.0" - - # Make sure that parfile exists - file = self.parfile() - + # Set name and version + self._name = "csobs2caldb" + self._version = "1.1.0" + + # Initialise members + self._obs = [] # Observation container + self._observation = None + self._mission = "CTA" + self._instrument = "" + self._rspname = "" + self._outfile = "" + self._base_dir = "" + self._cal_dir = "" + # 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,87 +81,76 @@ class csobs2caldb(ctools.cscript): # 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 that parfile was not found - print("Parfile \""+parfile+"\" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inobs","f","a","obs.xml","","","Input observation definition file")) - pars.append(gammalib.GApplicationPar("rspname","s","a","NONE","","","Response output name (e.g. Zenith50)")) - pars.append(gammalib.GApplicationPar("outirfs","f","h","irf_file.fits","","","Output IRF file name")) - pars.append(gammalib.GApplicationPar("index","i","h","0","","","Index of observation to be used from XML file")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csobs2caldb.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ - # Get parameters - self.inobs = self["inobs"].filename() - self.obs = gammalib.GObservations(self.inobs) + # Load observation container + self._obs = gammalib.GObservations(self["inobs"].filename()) # Get index of observation to be used (0=default) inx = self["index"].integer() # Get observation - self.observation = self.obs[inx] - - # Set instrument mission to cta - # This is to be compliant with standard gammalib caldb - self.m_mission = "CTA" + self._observation = self._obs[inx] # Get instrument name - self.m_instrument = self.observation.instrument() + self._instrument = self._observation.instrument() # Get config and file name - self.m_rspname = self["rspname"].string() - self.m_outfile = self["outirfs"].filename() + self._rspname = self["rspname"].string() + self._outfile = self["outirfs"].filename() # Make sure we have a cta observation - if not self.observation.classname() == "GCTAObservation": - raise gammalib.Exception.runtime_error("Input observation not of type \"GCTAObservation\"") + if not self._observation.classname() == "GCTAObservation": + raise RuntimeError("Input observation not of type "+ + "\"GCTAObservation\".") # Make sure we have an irf response - if not self.observation.response().classname() == "GCTAResponseIrf": - raise gammalib.Exception.runtime_error("Response of input observation not of type \"GCTAResponseIrf\"") + if not self._observation.response().classname() == "GCTAResponseIrf": + raise RuntimeError("Response of input observation not of "+ + "type \"GCTAResponseIrf\".") # Get other parameters - self.m_clobber = self["clobber"].boolean() + self._clobber = self["clobber"].boolean() # Return return - def make_irf_file(self): - - # retrieve response component - rsp = self.observation.response() - - # Initialise FITS file - fits = gammalib.GFits() + def _make_irf_file(self): + """ + Creates an IRF FITS file. + """ + # Write header into logger + if self._logTerse(): + self._log("\n") + self._log.header2("Creating IRF file") + + # Get response for the observation + rsp = self._observation.response() - # Create GFilename instances - fname_aeff = gammalib.GFilename(rsp.aeff().filename()) - fname_psf = gammalib.GFilename(rsp.psf().filename()) - fname_edisp = gammalib.GFilename(rsp.edisp().filename()) - fname_bkg = gammalib.GFilename(rsp.background().filename()) + # Extract response file names + fname_aeff = rsp.aeff().filename() + fname_psf = rsp.psf().filename() + fname_edisp = rsp.edisp().filename() + fname_bkg = rsp.background().filename() + + # Log filename instances + if self._logTerse(): + self._log(gammalib.parformat("Effective area")) + self._log(repr(fname_aeff)) + self._log("\n") + self._log(gammalib.parformat("Point spread function")) + self._log(fname_psf.url()) + self._log("\n") + self._log(gammalib.parformat("Energy dispersion")) + self._log(fname_edisp.url()) + self._log("\n") + self._log(gammalib.parformat("Background rate")) + self._log(fname_bkg.url()) + self._log("\n") # Open FITS files of response components fits_aeff = gammalib.GFits(fname_aeff) @@ -157,7 +158,10 @@ class csobs2caldb(ctools.cscript): fits_edisp = gammalib.GFits(fname_edisp) fits_bkg = gammalib.GFits(fname_bkg) - # Bundle IRFs into one file + # Create empty FITS file + fits = gammalib.GFits() + + # Append IRF component to FITS file fits.append(fits_aeff[fname_aeff.extname("EFFECTIVE AREA")]) fits.append(fits_psf[fname_psf.extname("POINT SPREAD FUNCTION")]) fits.append(fits_edisp[fname_edisp.extname("ENERGY DISPERSION")]) @@ -166,30 +170,37 @@ class csobs2caldb(ctools.cscript): # Return fits file return fits - def make_caldb_index(self): - + def _make_caldb_index(self): + """ + Creates an IRF index FITS file. + """ + # Write header into logger + if self._logTerse(): + self._log("\n") + self._log.header2("Creating database index file") + # Open calibration database index - indx_file = self.base_dir+"/caldb.indx" + indx_file = self._base_dir+"/caldb.indx" - # Check if indx file exists (create if not existent) - if os.path.isfile(indx_file): - cif = gammalib.GFits(indx_file) - else: - cif = gammalib.GFits(indx_file, True) + # Open index file (or create one if it does not exist) + cif = gammalib.GFits(indx_file, True) # Retrieve "CIF" table if cif.contains("CIF"): table = cif.table("CIF") + + # ... or create binary table if no "CIF" table exists yet else: - # Create binary table if not available + + # Create empty binary table bintable = gammalib.GFitsBinTable() bintable.extname("CIF") bintable.card("CIFVERSN", "1992a", "Version of CIF format") - cif.append(bintable) - - table = cif.table("CIF") + + # Append table to FITS file and recover a reference to the + # appended table + table = cif.append(bintable) - # Copied from cta_root2caldb.py # Attach columns. Reference: CAL/GEN/92-008 table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10)) table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10)) @@ -215,7 +226,7 @@ class csobs2caldb(ctools.cscript): row_index = -1 for caldir in table["CAL_DIR"]: row_index += 1 - if caldir == self.cal_dir: + if caldir == self._cal_dir: has_config = True break @@ -238,8 +249,8 @@ class csobs2caldb(ctools.cscript): today, time = now.split() # Set element - table["TELESCOP"][row] = self.m_mission - table["INSTRUME"][row] = self.m_instrument + table["TELESCOP"][row] = self._mission + table["INSTRUME"][row] = self._instrument table["DETNAM"][row] = "NONE" table["FILTER"][row] = "NONE" table["CAL_DEV"][row] = "ONLINE" @@ -249,115 +260,133 @@ class csobs2caldb(ctools.cscript): table["CAL_VST"][row] = time.split(".")[0] table["REF_TIME"][row] = 51544.0 table["CAL_QUAL"][row] = 0 - table["CAL_CBD"][row] = "NAME("+self.m_rspname+")" + table["CAL_CBD"][row] = "NAME("+self._rspname+")" table["CAL_DATE"][row] = today.replace("-","/")[2:] - table["CAL_DIR"][row] = self.cal_dir - table["CAL_FILE"][row] = self.m_outfile + table["CAL_DIR"][row] = self._cal_dir + table["CAL_FILE"][row] = self._outfile # Add effective area information row = row_index-4 - table["CAL_CNAM"][row] = "EFF_AREA" - table["CAL_DESC"][row] = self.m_instrument+" effective area" + table["CAL_CNAM"][row] = "EFF_AREA" + table["CAL_DESC"][row] = self._instrument+" effective area" # Add point spread function information row = row_index-3 - table["CAL_CNAM"][row] = "RPSF" - table["CAL_DESC"][row] = self.m_instrument+" point spread function" + table["CAL_CNAM"][row] = "RPSF" + table["CAL_DESC"][row] = self._instrument+" point spread function" # Add energy dispersion information row = row_index-2 - table["CAL_CNAM"][row] = "EDISP" - table["CAL_DESC"][row] = self.m_instrument+" energy dispersion" + table["CAL_CNAM"][row] = "EDISP" + table["CAL_DESC"][row] = self._instrument+" energy dispersion" # Add background information row = row_index-1 - table["CAL_CNAM"][row] = "BGD" - table["CAL_DESC"][row] = self.m_instrument+" background" + table["CAL_CNAM"][row] = "BGD" + table["CAL_DESC"][row] = self._instrument+" background" # Return CIF FITS file return cif - def make_dirs(self): + def _make_dirs(self): + """ + Make CALDB directories. + """ + # Write header into logger + if self._logTerse(): + self._log("\n") + self._log.header2("Creating directory structure") - # Create gammalib Caldb + # Create calibration database caldb = gammalib.GCaldb() # Set calibration directory - self.cal_dir = "data" - self.cal_dir += "/"+self.m_mission.lower() - self.cal_dir += "/"+self.m_instrument.lower() - self.cal_dir += "/bcf/"+self.m_rspname + self._cal_dir = "data" + self._cal_dir += "/"+self._mission.lower() + self._cal_dir += "/"+self._instrument.lower() + self._cal_dir += "/bcf/"+self._rspname # Set absolute path - self.base_dir = caldb.rootdir() +"/data" - self.base_dir += "/"+self.m_mission.lower() - self.base_dir += "/"+self.m_instrument.lower() + self._base_dir = caldb.rootdir() +"/data" + self._base_dir += "/"+self._mission.lower() + self._base_dir += "/"+self._instrument.lower() # Set directory for irf file - self.rsp_dir = caldb.rootdir() + "/" + self.cal_dir + self._rsp_dir = caldb.rootdir() + "/" + self._cal_dir # Create directories and log information - if not os.path.isdir(self.rsp_dir): - if self.logExplicit(): - self.log(gammalib.parformat("Directory")) - self.log(self.rsp_dir) - self.log("\n") - os.makedirs(self.rsp_dir) + if not os.path.isdir(self._rsp_dir): + if self._logExplicit(): + self._log(gammalib.parformat("Directory")) + self._log(self._rsp_dir) + self._log("\n") + os.makedirs(self._rsp_dir) else: - if self.logExplicit(): - self.log(gammalib.parformat("Directory (existing)")) - self.log(self.rsp_dir) - self.log("\n") + if self._logExplicit(): + self._log(gammalib.parformat("Directory (existing)")) + self._log(self._rsp_dir) + self._log("\n") # Return return + + # Public methods def run(self): """ Run the script. """ # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + 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\n") - self.log.header1("Creating directory structure") + # Write input parameters and header into logger + if self._logTerse(): + self._log_parameters() + self._log("\n\n") + self._log.header1("Creating CALDB entry") # Create directory structure - self.make_dirs() - - # Logging - if self.logTerse(): - self.log("\n") - self.log.header1("Updating calibration database") + self._make_dirs() # Create/update calibration database - self.caldb_inx = self.make_caldb_index() + self._caldb_inx = self._make_caldb_index() - # Logging - if self.logTerse(): - self.log("\n") - self.log.header1("Creating response file") - - # Create IRF file - self.irf_fits = self.make_irf_file() + # Create response file + self._irf_fits = self._make_irf_file() # Return return def save(self): + """ + Save calibration database FITS file. + """ + # Write header into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Save calibration database") + + # Set response filename + filename = self._rsp_dir + "/" + self._outfile + + # Write filenames into logger + if self._logTerse(): + self._log(gammalib.parformat("CALDB index file")) + self._log(self._caldb_inx.filename().url()) + self._log("\n") + self._log(gammalib.parformat("Response file")) + self._log(filename) + self._log("\n") # Save caldb index file - self.caldb_inx.save(self.m_clobber) + self._caldb_inx.save(self._clobber) # Save response file - self.irf_fits.saveto(self.rsp_dir + "/" + self.m_outfile, self.m_clobber) + self._irf_fits.saveto(filename, self._clobber) # Return return @@ -380,14 +409,12 @@ class csobs2caldb(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates caldb entry. - """ + # Create instance of application app = csobs2caldb(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csobsdef.par b/cscripts/csobsdef.par new file mode 100644 index 0000000..e4994c7 --- /dev/null +++ b/cscripts/csobsdef.par @@ -0,0 +1,41 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csobsdef.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csobsdef script # +# # +############################################################################# +# +# General parameters +#=================== +inpnt, f, a, NONE,,, "Input pointing definition ASCII file" +outobs, f, a, obs.xml,,, "Output observation definition XML file" + +# +# Common pointing parameter values +#================================= +duration, r, a, 1800.0,,, "Pointing duration (seconds)" +caldb, s, h, ,,, "Calibration database" +irf, s, h, ,,, "Instrument response function" +emin, r, h, UNDEF,,, "Lower energy limit (TeV)" +emax, r, h, UNDEF,,, "Upper energy limit (TeV)" +rad, r, h, UNDEF,,, "Radius of ROI (degrees)" +deadc, r, h, 0.95,,, "Deadtime correction factor" + +# +# 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, csobsdef.log,,, "Log filename" diff --git a/cscripts/csobsdef.py b/cscripts/csobsdef.py index cfd0f94..67f4e0c 100755 --- a/cscripts/csobsdef.py +++ b/cscripts/csobsdef.py @@ -28,11 +28,12 @@ import sys # ============== # class csobsdef(ctools.cscript): """ - The csobsdef class generates an observation definition file from a - pointing list. The pointing list is a comma-separated value ASCII - file with header keywords in the first row followed by a list of - pointings (one pointing per row). The following header keywords - are supported (case sensitive, column order irrelevant): + Creates an observation definition file from a pointing list. + + The pointing list is a comma-separated value ASCII file with header + keywords in the first row followed by a list of pointings (one + pointing per row). The following header keywords are supported (case + sensitive, column order irrelevant): name - Observation name string id - Unique observation identifier string @@ -53,53 +54,54 @@ class csobsdef(ctools.cscript): csobsdef as user parameters. The only exception is the "duration" keyword that will automatically be queried. - Here some usage examples: + Examples: - csobsdef - Creates minimal observation definition file. + ./csobsdef + Creates minimal observation definition file. - csobsdef emin=0.1 emax=100.0 - Creates observation definition file with an energy range 100 GeV - 100 TeV. + ./csobsdef emin=0.1 emax=100.0 + Creates observation definition file with an energy range + 100 GeV - 100 TeV. - csobsdef rad=5 - Creates observation definition file with a ROI radius of 5 deg. + ./csobsdef rad=5 + Creates observation definition file with a ROI radius of + 5 degrees. - csobsdef caldb=dummy irf=cta_dummy_irf - Creates observation definition file using the "cta_dummy_irf" IRF in the - "dummy" calibration database. + ./csobsdef caldb=dummy irf=cta_dummy_irf + Creates observation definition file using the "cta_dummy_irf" + IRF in the "dummy" calibration database. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name and version - self.name = "csobsdef" - self.version = "1.0.0" + self._name = "csobsdef" + self._version = "1.1.0" # Initialise class members - self.obs = gammalib.GObservations() - self.inpnt = "" - self.tmin = 0.0 - self.outfile = "" - self.chatter = 2 - self.clobber = True - self.debug = False - self.read_ahead = False - - # Make sure that parfile exists - file = self.parfile() + 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) + 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 @@ -111,92 +113,78 @@ class csobsdef(ctools.cscript): # 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" - - # Test parfile existence, and create parfile if it does not exist - try: - pars = gammalib.GApplicationPars(parfile) - except: - # Signal that parfile was not found - print("Parfile \""+parfile+"\" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inpnt","f","a","NONE","","","Input pointing definition file")) - pars.append(gammalib.GApplicationPar("outobs","f","a","obs.xml","","","Output observation definition XML file")) - pars.append(gammalib.GApplicationPar("duration","r","a","1800.0","","","Pointing duration (seconds)")) - pars.append(gammalib.GApplicationPar("caldb","s","h","","","","Calibration database")) - pars.append(gammalib.GApplicationPar("irf","s","h","","","","Instrument response function")) - pars.append(gammalib.GApplicationPar("emin","r","h","UNDEF","","","Lower energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","h","UNDEF","","","Upper energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("rad","r","h","UNDEF","","","Radius of ROI (degrees)")) - pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csobsdef.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + def _get_parameters(self): """ Get parameters from parfile. """ # Read parameters - self.inpnt = self["inpnt"].filename() + self._inpnt = self["inpnt"].filename() # Read ahead parameters - if self.read_ahead: - self.outfile = self["outobs"].filename() + if self._read_ahead: + self._outobs = self["outobs"].filename() # Set some fixed parameters - self.chatter = self["chatter"].integer() - self.clobber = self["clobber"].boolean() - self.debug = self["debug"].boolean() + self._chatter = self["chatter"].integer() + self._clobber = self["clobber"].boolean() + self._debug = self["debug"].boolean() + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - def execute(self): + def _set_response(self, obs, caldb, irf): """ - Execute the script. + Set response for an observation. + + The method creates an XML element so that that the response XML + writer will write the database and response name into the + observation definition file. """ - # Set read ahead flag - self.read_ahead = True + # Create XML element + xml = gammalib.GXmlElement() - # Run the script - self.run() + # Append parameter + parameter = "parameter name=\"Calibration\" database=\""+\ + caldb+"\" response=\""+irf+"\"" + xml.append(gammalib.GXmlElement(parameter)) - # Save observation definition file - self.obs.save(self.outfile) + # Create CTA response + response = gammalib.GCTAResponseIrf() + response.read(xml) - # Return - return + # Attach response to observation + obs.response(response) + + # Return observation + return obs + + # 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() + self._get_parameters() - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + # Write header into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Creating observation definition XML file") # Load pointing definition file - pntdef = gammalib.GCsv(self.inpnt, ",") + pntdef = gammalib.GCsv(self._inpnt, ",") ncols = pntdef.ncols() npnt = pntdef.nrows()-1 @@ -205,8 +193,8 @@ class csobsdef(ctools.cscript): raise RuntimeError("No header found in pointing definition file.") # Clear observation container - self.obs.clear() - self.id = 1 + self._obs.clear() + self._id = 1 # Extract header from pointing definition file header = [] @@ -233,8 +221,8 @@ class csobsdef(ctools.cscript): if "id" in header: id = pntdef[row, header.index("id")] else: - id = "%6.6d" % self.id - self.id += 1 + id = "%6.6d" % self._id + self._id += 1 obs.id(id) # Set pointing @@ -249,8 +237,8 @@ class csobsdef(ctools.cscript): pntdir = gammalib.GSkyDir() pntdir.lb_deg(lon,lat) else: - raise RuntimeError("No (ra,dec) or (lon,lat) columns found in pointing" - " definition file.") + raise RuntimeError("No (ra,dec) or (lon,lat) columns " + "found in pointing definition file.") obs.pointing(gammalib.GCTAPointing(pntdir)) # Set response function @@ -263,7 +251,7 @@ class csobsdef(ctools.cscript): else: irf = self["irf"].string() if caldb != "" and irf != "": - obs = self.__set_response(obs, caldb, irf) + obs = self._set_response(obs, caldb, irf) # Set deadtime correction factor if "deadc" in header: @@ -277,12 +265,12 @@ class csobsdef(ctools.cscript): duration = float(pntdef[row, header.index("duration")]) else: duration = self["duration"].real() - tmin = self.tmin - tmax = self.tmin + duration - gti = gammalib.GGti(self._time_reference()) - tstart = gammalib.GTime(tmin, self._time_reference()) - tstop = gammalib.GTime(tmax, self._time_reference()) - self.tmin = tmax + tmin = self._tmin + tmax = self._tmin + duration + gti = gammalib.GGti(self._time_reference()) + tstart = gammalib.GTime(tmin, self._time_reference()) + tstop = gammalib.GTime(tmax, self._time_reference()) + self._tmin = tmax gti.append(tstart, tstop) obs.ontime(gti.ontime()) obs.livetime(gti.ontime()*deadc) @@ -334,49 +322,75 @@ class csobsdef(ctools.cscript): # Attach event list to CTA observation obs.events(list) + # Write observation into logger + if self._logExplicit(): + self._log(str(obs)) + self._log("\n") + elif self._logTerse(): + self._log(gammalib.parformat(obs.instrument()+" observation")) + self._log("Name=\""+obs.name()+"\" ") + self._log("ID=\""+obs.id()+"\"\n") + # Append observation - self.obs.append(obs) + self._obs.append(obs) # Return return - def __set_response(self, obs, caldb, irf): + def save(self): """ - Set response for an observation. We create an XML element for that - so that the response XML writer will write the database and response - name into the observation definition file. + Save observation definition XML file. """ - # Create XML element - xml = gammalib.GXmlElement() + # Write header and filename into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Save observation definition XML file") + + # Get output filename in case it was not read ahead + self._outobs = self["outobs"].filename() + + # Check if observation definition XML file is valid + if self._outobs != "NONE": + + # Log filename + if self._logTerse(): + self._log(gammalib.parformat("Observation XML file")) + self._log(self._outobs.url()) + self._log("\n") + + # Save observation definition XML file + self._obs.save(self._outobs) + + # Return + return - # Append parameter - parameter = "parameter name=\"Calibration\" database=\""+\ - caldb+"\" response=\""+irf+"\"" - xml.append(gammalib.GXmlElement(parameter)) + def execute(self): + """ + Execute the script. + """ + # Set read ahead flag + self._read_ahead = True - # Create CTA response - response = gammalib.GCTAResponseIrf() - response.read(xml) + # Run the script + self.run() - # Attach response to observation - obs.response(response) + # Save observation definition file + self.save() - # Return observation - return obs + # Return + return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates observation definition file. - """ + # Create instance of application app = csobsdef(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csobsinfo.par b/cscripts/csobsinfo.par new file mode 100644 index 0000000..23a3ab5 --- /dev/null +++ b/cscripts/csobsinfo.par @@ -0,0 +1,37 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csobsinfo.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csobsinfo script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inobs, f, a, obs.xml,,, "Input event list, counts cube, or observation definition XML file" +ds9file, f, h, NONE,,, "Output DS9 region file containing pointings" + +# +# Script parameters +#================== +offset, b, h, no,,, "Compute offset from target to pointing positions?" +ra, r, a, 0.0,,, "Target Right Ascension (deg)" +dec, r, a, 0.0,,, "Target Declination (deg)") + +# +# 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, csobsinfo.log,,, "Log filename" diff --git a/cscripts/csobsinfo.py b/cscripts/csobsinfo.py index ca0b09b..fe6418d 100755 --- a/cscripts/csobsinfo.py +++ b/cscripts/csobsinfo.py @@ -28,40 +28,48 @@ import sys # =============== # class csobsinfo(ctools.cscript): """ - This class dumps information about an observation container into a - logfile or on screen. This might be helpful for quick access to an - observation container showing, e.g., total lifetime, energy range etc. + Shows the content of an observation container. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "csobsinfo" - self.version = "1.1.0" + self._name = "csobsinfo" + self._version = "1.1.0" - # Initialise some members + # Initialise class members + self._obj_dir = None + self._compute_offset = False + self._ds9file = "" + self._offsets = [] + self._zeniths = [] + self._azimuths = [] + self._pnt_ra = [] + 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:] + self._obs = argv[0] + argv = argv[1:] else: - self.obs = gammalib.GObservations() - self.obs.clear() - - # Make sure that parfile exists - file = self.parfile() + self._obs = gammalib.GObservations() + self._obs.clear() # 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 @@ -73,89 +81,64 @@ class csobsinfo(ctools.cscript): # 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 that parfile was not found - print("Parfile \""+parfile+"\" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inobs","f","a","obs.xml","","","Event list, counts cube, or observation definition file")) - pars.append(gammalib.GApplicationPar("offset","b","h","no","","","Compute offset from target to pointing positions")) - pars.append(gammalib.GApplicationPar("ra","r","a","0.0","","","Target Right Ascension")) - pars.append(gammalib.GApplicationPar("dec","r","a","0.0","","","Target Declination")) - pars.append(gammalib.GApplicationPar("ds9file","f","h","NONE","","","DS9 region file containing pointings")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csobsinfo.log","","","Log filename")) - pars.save(parfile) - # Return - return - - def get_parameters(self): + # Private methods + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Get parameters - if self.obs.size() == 0: + if self._obs.size() == 0: self._require_inobs("csobsinfo::get_parameters") - self.obs = self._get_observations(False) + self._obs = self._get_observations(False) # Initialise object position - self.obj_dir = gammalib.GSkyDir() + self._obj_dir = gammalib.GSkyDir() # Get (optional) offset parameters - self.m_offset = self["offset"].boolean() - if self.m_offset: + self._compute_offset = self["offset"].boolean() + if self._compute_offset: ra = self["ra"].real() dec = self["dec"].real() - self.obj_dir.radec_deg(ra,dec) + self._obj_dir.radec_deg(ra,dec) + + # Read ahead DS9 filename + if self._read_ahead: + self._ds9file = self["ds9file"].filename() - # Get (optional) DS9 filename - self.ds9file = self["ds9file"].filename() + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return + + # Public methods def run(self): """ Run the script. """ # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + 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() # Initialise arrays to store certain values for reuse # Todo, think about using a python dictionary - self.m_offsets = [] - self.m_zeniths = [] - self.m_azimuths = [] - self.pnt_ra = [] - self.pnt_dec = [] - obs_names = [] + self._offsets = [] + self._zeniths = [] + self._azimuths = [] + self._pnt_ra = [] + self._pnt_dec = [] + self._ebounds = gammalib.GEbounds() + self._gti = gammalib.GGti() + obs_names = [] # Initialise output to be filled - self.m_ebounds = gammalib.GEbounds() - self.m_gti = gammalib.GGti() ontime = 0.0 livetime = 0.0 n_events = 0 @@ -163,27 +146,30 @@ class csobsinfo(ctools.cscript): n_obs_binned = 0 n_obs_unbinned = 0 - # Logging - self.log.header1("Looping over observations") + # Write header + if self._logTerse(): + self._log("\n") + if self._obs.size() > 1: + self._log.header1("Observations") + else: + self._log.header1("Observation") # Loop over observations - for obs in self.obs: - + for obs in self._obs: # Skip non-CTA observations if not obs.classname() == "GCTAObservation": - self.log("Skipping "+obs.instrument()+" observation\n") + self._log("Skipping "+obs.instrument()+" observation\n") continue - # Use observed object as observation name - # if name is not given + # Use observed object as observation name if name is not given obs_name = obs.name() if obs_name == "": obs_name = obs.object() # Logging - if self.logTerse(): - self.log.header2(obs_name) + if self._logTerse(): + self._log.header2(obs_name) # Retrieve observation name obs_names.append(obs_name) @@ -199,206 +185,231 @@ class csobsinfo(ctools.cscript): # Retrieve pointing and store Ra,Dec pnt_dir = obs.pointing().dir() - self.pnt_ra.append(pnt_dir.ra_deg()) - self.pnt_dec.append(pnt_dir.dec_deg()) + self._pnt_ra.append(pnt_dir.ra_deg()) + self._pnt_dec.append(pnt_dir.dec_deg()) # If avaliable append energy boundaries if obs_bounds.size() > 0 : - self.m_ebounds.append(obs_bounds.emin(),obs_bounds.emax()) + self._ebounds.append(obs_bounds.emin(),obs_bounds.emax()) # Append time interval - self.m_gti.append(obs_gti.tstart(), obs_gti.tstop()) + self._gti.append(obs_gti.tstart(), obs_gti.tstop()) # Increment global livetime and ontime ontime += obs.ontime() livetime += obs.livetime() - # Log the event type (binned or unbinned) - if self.logTerse(): - self.log.parformat("binned") - + # Bookkeeping if obs.eventtype() == "CountsCube": n_eventbins += obs.events().size() n_obs_binned += 1 - if self.logTerse(): - self.log("yes") - self.log("\n") - self.log.parformat("Number of bins") - self.log(str(obs.events().size())) - self.log("\n") + if self._logTerse(): + self._log.parformat("Binned") + self._log("yes\n") + self._log.parformat("Number of bins") + self._log(str(obs.events().size())) + self._log("\n") else: n_events += obs.events().size() n_obs_unbinned += 1 - if self.logTerse(): - self.log("no") - self.log("\n") - self.log.parformat("Number of events") - self.log(str(obs.events().size())) - self.log("\n") + if self._logTerse(): + self._log.parformat("Binned") + self._log("no\n") + self._log.parformat("Number of events") + self._log(str(obs.events().size())) + self._log("\n") + + # Retrieve zenith and azimuth and store for later use + zenith = obs.pointing().zenith() + azimuth = obs.pointing().azimuth() + self._zeniths.append(zenith) + self._azimuths.append(azimuth) + + # Optionally compute offset with respect to target direction + if self._compute_offset: + offset = pnt_dir.dist_deg(self._obj_dir) + self._offsets.append(offset) + else: + self._offsets.append(-1.0) - # Log option - if self.logTerse(): + # Optionally log details + if self._logTerse(): # Log the observation energy range (if available) - self.log.parformat("Energy range") + self._log.parformat("Energy range") if obs_bounds.size() == 0: - self.log("undefined") + self._log("undefined") else: - self.log(str(obs_bounds.emin())+" - "+str(obs_bounds.emax())) - self.log("\n") + self._log(str(obs_bounds.emin())) + self._log(" - ") + self._log(str(obs_bounds.emax())) + self._log("\n") # Log observation time interval - self.log.parformat("Time range (MJD)") + self._log.parformat("Time range (MJD)") if obs_gti.size() == 0: - self.log("undefined") + self._log("undefined") else: - self.log(str(obs_gti.tstart().mjd())+" - "+str(obs_gti.tstop().mjd())) - self.log("\n") + self._log(str(obs_gti.tstart().mjd())) + self._log(" - ") + self._log(str(obs_gti.tstop().mjd())) + self._log("\n") # Log ontime - self.log.parformat("Ontime [s]") - self.log(str(obs.ontime())) - self.log("\n") + self._log.parformat("Ontime") + self._log(str(obs.ontime())) + self._log(" s\n") # Log livetime - self.log.parformat("Livetime [s]") - self.log(str(obs.livetime())) - self.log("\n") + self._log.parformat("Livetime") + self._log(str(obs.livetime())) + self._log(" s\n") # Log dead time fraction - self.log.parformat("Dead-time fraction (%)") - self.log("%.3f" % (deadfrac)) - self.log("\n") + self._log.parformat("Deadtime fraction (%)") + self._log("%.3f" % (deadfrac)) + self._log("\n") # Log pointing direction - self.log.parformat("Pointing") - self.log(str(pnt_dir)) - self.log("\n") + self._log.parformat("Pointing") + self._log(str(pnt_dir)) + self._log("\n") - # Log offset if possible - if self.m_offset: - offset = pnt_dir.dist_deg(self.obj_dir) - self.m_offsets.append(offset) - self.log.parformat("Offset") - self.log("%.2f" % (offset)) - self.log("\n") - else: - self.m_offsets.append(-1.0) + # Optionally log offset with respect to target direction + if self._compute_offset: + self._log.parformat("Offset from target") + self._log("%.2f" % (offset)) + self._log(" deg\n") - # Retrieve zenith and azimuth and store for later use - zenith = obs.pointing().zenith() - azimuth = obs.pointing().azimuth() - self.m_zeniths.append(zenith) - self.m_azimuths.append(azimuth) - # Log Zenith and Azimuth if required - if self.logExplicit(): - self.log.parformat("Zenith angle") - self.log("%.2f" % (zenith)) - self.log("\n") - self.log.parformat("Azimuth angle") - self.log("%.2f" % (azimuth)) - self.log("\n\n") + self._log.parformat("Zenith angle") + self._log("%.2f" % (zenith)) + self._log(" deg\n") + self._log.parformat("Azimuth angle") + self._log("%.2f" % (azimuth)) + self._log(" deg\n") # Log summary - self.log("\n") - self.log.header1("Summary") + if self._logTerse(): + + # Write header + self._log("\n") + self._log.header1("Summary") - # Log general summary - self.log.parformat("Unbinned observations") - self.log(str(n_obs_unbinned)) - self.log("\n") - self.log.parformat("Number of events") - self.log(str(n_events)) - self.log("\n") - self.log.parformat("Binned observations") - self.log(str(n_obs_binned)) - self.log("\n") - self.log.parformat("Number of bins") - self.log(str(n_eventbins)) - self.log("\n\n") + # Log general summary + self._log.header3("Observations") + self._log.parformat("Unbinned observations") + self._log(str(n_obs_unbinned)) + self._log("\n") + self._log.parformat("Binned observations") + self._log(str(n_obs_binned)) + self._log("\n") + self._log.header3("Events") + self._log.parformat("Number of events") + self._log(str(n_events)) + self._log("\n") + self._log.parformat("Number of bins") + self._log(str(n_eventbins)) + self._log("\n") - # Log pointing summary - self.log.header3("Pointings") + # Log pointing summary + self._log.header3("Pointings") - # Log mean offset if possible - if self.m_offset: - self.log.parformat("Mean offset angle") - self.log("%.2f" % (sum(self.m_offsets)/len(self.m_offsets))) - self.log("\n") + # Log mean offset if possible + if self._compute_offset: + self._log.parformat("Mean offset angle") + self._log("%.2f" % (sum(self._offsets)/len(self._offsets))) + self._log(" deg\n") - # Log mean azimuth and zenith angle - self.log.parformat("Mean zenith angle") - self.log("%.2f" % (sum(self.m_zeniths)/len(self.m_zeniths))) - self.log("\n") - self.log.parformat("Mean azimuth angle") - self.log("%.2f" % (sum(self.m_azimuths)/len(self.m_azimuths))) - self.log("\n") + # Log mean azimuth and zenith angle + self._log.parformat("Mean zenith angle") + self._log("%.2f" % (sum(self._zeniths)/len(self._zeniths))) + self._log(" deg\n") + self._log.parformat("Mean azimuth angle") + self._log("%.2f" % (sum(self._azimuths)/len(self._azimuths))) + self._log(" deg\n") - # Log name of observations if requested - if self.logExplicit(): - obs_set = set(obs_names) - for name in obs_set: - self.log.parformat("\""+name+"\"") - self.log(str(obs_names.count(name))) - self.log("\n") - self.log("\n") + # Optionally log name of observations + if self._logExplicit(): + obs_set = set(obs_names) + for name in obs_set: + self._log.parformat("\""+name+"\"") + self._log(str(obs_names.count(name))) + self._log("\n") + self._log("\n") - # Log energy range - self.log.header3("Energy range") - self.log.parformat("Emin") - if self.m_ebounds.size() == 0: - self.log("undefined") - else: - self.log(str(self.m_ebounds.emin())) - self.log("\n") - self.log.parformat("Emax") - if self.m_ebounds.size() == 0: - self.log("undefined") - else: - self.log(str(self.m_ebounds.emax())) - self.log("\n\n") + # Log energy range + self._log.header3("Energy range") + self._log.parformat("Minimum energy") + if self._ebounds.size() == 0: + self._log("undefined") + else: + self._log(str(self._ebounds.emin())) + self._log("\n") + self._log.parformat("Maximum energy") + if self._ebounds.size() == 0: + self._log("undefined") + else: + self._log(str(self._ebounds.emax())) + self._log("\n") - # Log time range - self.log.header3("Time range") - self.log.parformat("Start [MJD]") - self.log(str(self.m_gti.tstart().mjd())) - self.log("\n") - self.log.parformat("Stop [MJD]") - self.log(str(self.m_gti.tstop().mjd())) - self.log("\n") + # Log time range + self._log.header3("Time range") + self._log.parformat("Start (MJD)") + self._log(str(self._gti.tstart().mjd())) + self._log("\n") + self._log.parformat("Stop (MJD)") + self._log(str(self._gti.tstop().mjd())) + self._log("\n") - # Log ontime and livetime in different units - self.log.parformat("Total ontime") - self.log("%.2f s = %.2f min = %.2f h" % (ontime, ontime/60., ontime/3600.)) - self.log("\n") - self.log.parformat("Total livetime") - self.log("%.2f s = %.2f min = %.2f h" % (livetime, livetime/60.0, livetime/3600.)) - self.log("\n") + # Log ontime and livetime in different units + self._log.parformat("Total ontime") + self._log("%.2f s = %.2f min = %.2f h" % (ontime, ontime/60., ontime/3600.)) + self._log("\n") + self._log.parformat("Total livetime") + self._log("%.2f s = %.2f min = %.2f h" % (livetime, livetime/60.0, livetime/3600.)) + self._log("\n") # Return return def save(self): """ - Save pointings to ds9 region file if required + Save pointings into DS9 region file. + + This method saves all pointing directions that are found in the + observation container into a DS9 region file. If "NONE" is + specified for the "ds9file" parameter the method does nothing. """ - - # Check if ds9 file is valid - if not self.ds9file == "NONE": + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save pointings in DS9 file") + + # Get output filename in case it was not read ahead + self._ds9file = self["ds9file"].filename() + + # Check if DS9 file is valid + if self._ds9file != "NONE": + + # Log filename + if self._logTerse(): + self._log(gammalib.parformat("DS9 filename")) + self._log(self._ds9file.url()) + self._log("\n") # Open file - f = open(self.ds9file.url(),"w") + f = open(self._ds9file.url(),"w") # Write coordinate system f.write("fk5\n") # Loop over pointings - for i in range(len(self.pnt_ra)): + for i in range(len(self._pnt_ra)): # Create string line = "point(" - line+=str(self.pnt_ra[i])+","+str(self.pnt_dec[i])+")" + line+=str(self._pnt_ra[i])+","+str(self._pnt_dec[i])+")" line+=" # point=cross 20 width=3 \n" # Write to file @@ -410,41 +421,43 @@ class csobsinfo(ctools.cscript): # Return return - def zeniths(self): """ Return zenith angles """ - return self.m_zeniths + return self._zeniths def azimuths(self): """ Return azimuth angles """ - return self.m_azimuths + return self._azimuths def offsets(self): """ Return offset angles """ - return self.m_offsets + return self._offsets def ebounds(self): """ Return energy boundaries """ - return self.m_ebounds + return self._ebounds def gti(self): """ Return good time intervals """ - return self.m_gti + return self._gti def execute(self): """ Execute the script. """ + # Set read ahead flag + self._read_ahead = True + # Run the script self.run() @@ -459,14 +472,12 @@ class csobsinfo(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates information about an observation container. - """ + # Create instance of application app = csobsinfo(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/cspull.par b/cscripts/cspull.par new file mode 100644 index 0000000..638886f --- /dev/null +++ b/cscripts/cspull.par @@ -0,0 +1,55 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cspull.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the cspull script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inobs, f, h, NONE,,, "Input event list, counts cube, 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?" +deadc, r, h, 0.95,0,1, "Deadtime correction factor" +outfile, f, a, pull.dat,,, "Output pull distribution file" + +# +# Script parameters +#================== +ntrials, i, a, 10,,, "Number of trials" +ra, r, a, 83.6331,0,360, "RA of pointing (deg)" +dec, r, a, 22.0145,-90,90, "Dec of pointing (deg)" +coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)" +proj, s, a, TAN,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (TeV)" +enumbins, i, a, 0,,, "Number of energy bins (0=unbinned)" +tmin, r, h, 0.0,,, "Start time (MET in s)" +tmax, r, a, 1800.0,,, "Duration (s)" +npix, i, a, 200,,, "Number of pixels for binned" +binsz, r, a, 0.05,,, "Pixel size for binned (deg/pixel)" +rad, r, h, 5.0,0,180, "Radius of ROI (deg)" +pattern, s, h, single,single|four,, "Observation pattern" +offset, r, h, 1.5,,, "Observation pattern offset (deg)" +profile, b, h, no,,, "Use likelihood profile method for errors?" + +# +# 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, cspull.log,,, "Log filename" diff --git a/cscripts/cspull.py b/cscripts/cspull.py index fd0baf2..1e05f3f 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -30,47 +30,49 @@ import csv # ============ # class cspull(ctools.cscript): """ + Generates pull distributions for a model. + This class implements the pull distribution generation script. 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 = "cspull" - self.version = "1.0.0" + self._name = "cspull" + self._version = "1.1.0" # Initialise some members - self.obs = None - self.model = None - self.m_inmodel = None - self.m_edisp = False - self.m_profile = False - self.m_exposure = None - self.m_psfcube = None - self.m_bckcube = None - self.m_edispcube = None - self.m_stackmodels = None - self.m_coordsys = "CEL" - self.m_proj = "TAN" - - # Make sure that parfile exists - file = self.parfile() + self._obs = None + self._model = None + self._inmodel = None + self._edisp = False + self._profile = False + self._exposure = None + self._psfcube = None + self._bckcube = None + self._edispcube = None + self._stackmodels = None + self._coordsys = "CEL" + self._proj = "TAN" + self._log_clients = False # 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 @@ -79,348 +81,220 @@ class cspull(ctools.cscript): """ 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 - print("Parfile "+parfile+" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","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("caldb","s","a","prod2","","","Calibration database")) - pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) - pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) - pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor")) - pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output pull distribution file")) - pars.append(gammalib.GApplicationPar("profile","b","h","no","","","Use likelihood profile method for errors?")) - pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) - pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)")) - pars.append(gammalib.GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)")) - pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate system (CEL - celestial, GAL - galactic)")) - pars.append(gammalib.GApplicationPar("proj","s","a","TAN","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method")) - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) - pars.append(gammalib.GApplicationPar("tmin","r","h","0.0","","","Start time (MET in s)")) - pars.append(gammalib.GApplicationPar("tmax","r","a","1800.0","","","Duration (in s)")) - pars.append(gammalib.GApplicationPar("npix","i","a","200","","","Number of pixels for binned")) - pars.append(gammalib.GApplicationPar("binsz","r","a","0.05","","","Pixel size for binned (deg/pixel)")) - pars.append(gammalib.GApplicationPar("rad","r","h","5.0","0","180","Radius of ROI (deg)")) - pars.append(gammalib.GApplicationPar("pattern","s","h","single","","","Observation pattern (single/four)")) - pars.append(gammalib.GApplicationPar("offset","r","h","1.5","","","Observation pattern offset (deg)")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","cspull.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + def _get_parameters(self): """ - Get parameters from parfile and setup the observation. + Get parameters from parfile. """ # Get parameters # Set observation if not done before - if self.obs == None or self.obs.size() == 0: - self.obs = self._get_observations() + if self._obs == None or self._obs.size() == 0: + self._obs = self._get_observations() # Check for requested pattern and use above # observation parameters to set wobble pattern - self.m_pattern = self["pattern"].string() - if self.m_pattern == "four": - self.obs = self.set_obs() + self._pattern = self["pattern"].string() + if self._pattern == "four": + self._obs = self._set_obs() # Get number of energy bins - self.m_enumbins = self["enumbins"].integer() + self._enumbins = self["enumbins"].integer() # Read parameters for binned if requested - if not self.m_enumbins == 0: - self.m_npix = self["npix"].integer() - self.m_binsz = self["binsz"].real() - self.m_coordsys = self["coordsys"].string() - self.m_proj = self["proj"].string() + if not self._enumbins == 0: + self._npix = self["npix"].integer() + self._binsz = self["binsz"].real() + self._coordsys = self["coordsys"].string() + self._proj = self["proj"].string() else: # Set dummy values (required by obsutils) - self.m_npix = 0 - self.m_binsz = 0.0 + self._npix = 0 + self._binsz = 0.0 # 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()) # Read other parameters - self.m_outfile = self["outfile"].filename() - self.m_ntrials = self["ntrials"].integer() - self.m_edisp = self["edisp"].boolean() - self.m_offset = self["offset"].real() - self.m_profile = self["profile"].boolean() + self._outfile = self["outfile"].filename() + self._ntrials = self["ntrials"].integer() + self._edisp = self["edisp"].boolean() + self._offset = self["offset"].real() + self._profile = self["profile"].boolean() + self._chatter = self["chatter"].integer() + self._debug = self["debug"].boolean() - # Set some fixed parameters - self.m_log = False # Logging in client tools - self.m_chatter = self["chatter"].integer() # Chatter level - self.m_debug = self["debug"].boolean() # Debugging in client tools + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - - def models(self, models): + def _set_obs(self): """ - Set model. - """ - # Copy models - self.obs.models(models.copy()) - - # Return - return - + Set CTA observation(s). - def execute(self): - """ - Execute the script. + Returns: + An observation container. """ - # Run the script - self.run() + # Setup observation definition list + obsdeflist = obsutils.set_obs_patterns(self._pattern, + ra=self["ra"].real(), + dec=self["dec"].real(), + offset=self["offset"].real()) - # Return - return + # Create list of observations + obs = obsutils.set_obs_list(obsdeflist, + tstart=self["tmin"].real(), + duration=self["tmax"].real()-self["tmin"].real(), + deadc=self["deadc"].real(), + emin=self["emin"].real(), + emax=self["emax"].real(), + rad=self["rad"].real(), + irf=self["irf"].string(), + caldb=self["caldb"].string()) + # Return observation container + return obs - def set_stacked_irf(self): + def _set_stacked_irf(self): """ For stacked analysis prepare stacked irfs """ # Write header into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Compute stacked response") + if self._logTerse(): + self._log("\n") + self._log.header1("Compute stacked response") # Get stacked exposure - expcube = ctools.ctexpcube(self.obs) + expcube = ctools.ctexpcube(self._obs) expcube["incube"] = "NONE" expcube["usepnt"] = True expcube["ebinalg"] = "LOG" - expcube["binsz"] = self.m_binsz - expcube["nxpix"] = self.m_npix - expcube["nypix"] = self.m_npix - expcube["enumbins"] = self.m_enumbins + expcube["binsz"] = self._binsz + expcube["nxpix"] = self._npix + expcube["nypix"] = self._npix + expcube["enumbins"] = self._enumbins expcube["emin"] = self["emin"].real() expcube["emax"] = self["emax"].real() - expcube["coordsys"] = self.m_coordsys - expcube["proj"] = self.m_proj + expcube["coordsys"] = self._coordsys + expcube["proj"] = self._proj expcube.run() # Notify exposure computation - if self.logTerse(): - self.log("Computed exposure cube\n") + if self._logTerse(): + self._log("Computed exposure cube\n") # Get stacked Psf - psfcube = ctools.ctpsfcube(self.obs) + psfcube = ctools.ctpsfcube(self._obs) psfcube["incube"] = "NONE" psfcube["usepnt"] = True psfcube["ebinalg"] = "LOG" - psfcube["binsz"] = self.m_binsz*10.0 - psfcube["nxpix"] = self.m_npix/10 - psfcube["nypix"] = self.m_npix/10 - psfcube["enumbins"] = self.m_enumbins + psfcube["binsz"] = self._binsz*10.0 + psfcube["nxpix"] = self._npix/10 + psfcube["nypix"] = self._npix/10 + psfcube["enumbins"] = self._enumbins psfcube["emin"] = self["emin"].real() psfcube["emax"] = self["emax"].real() - psfcube["coordsys"] = self.m_coordsys - psfcube["proj"] = self.m_proj + psfcube["coordsys"] = self._coordsys + psfcube["proj"] = self._proj psfcube.run() # Notify Psf computation - if self.logTerse(): - self.log("Computed Psf cube\n") + if self._logTerse(): + self._log("Computed point spread function cube\n") # Optionally get stacked Edisp - if self.m_edisp: - edispcube = ctools.edispcube(self.obs) + if self._edisp: + edispcube = ctools.edispcube(self._obs) edispcube["incube"] = "NONE" edispcube["usepnt"] = True edispcube["ebinalg"] = "LOG" - edispcube["binsz"] = self.m_binsz*10.0 - edispcube["nxpix"] = self.m_npix/10 - edispcube["nypix"] = self.m_npix/10 - edispcube["enumbins"] = self.m_enumbins + edispcube["binsz"] = self._binsz*10.0 + edispcube["nxpix"] = self._npix/10 + edispcube["nypix"] = self._npix/10 + edispcube["enumbins"] = self._enumbins edispcube["emin"] = self["emin"].real() edispcube["emax"] = self["emax"].real() - edispcube["coordsys"] = self.m_coordsys - edispcube["proj"] = self.m_proj + edispcube["coordsys"] = self._coordsys + edispcube["proj"] = self._proj edispcube.run() # Store result - self.m_edispcube = edispcube.edispcube().copy() + self._edispcube = edispcube.edispcube().copy() # Logging - if self.logTerse(): - self.log("Computed energy dispersion cube\n") + if self._logTerse(): + self._log("Computed energy dispersion cube\n") # Get stacked background - bkgcube = ctools.ctbkgcube(self.obs) + bkgcube = ctools.ctbkgcube(self._obs) bkgcube["incube"] = "NONE" bkgcube["usepnt"] = True bkgcube["ebinalg"] = "LOG" - bkgcube["binsz"] = self.m_binsz - bkgcube["nxpix"] = self.m_npix - bkgcube["nypix"] = self.m_npix - bkgcube["enumbins"] = self.m_enumbins + bkgcube["binsz"] = self._binsz + bkgcube["nxpix"] = self._npix + bkgcube["nypix"] = self._npix + bkgcube["enumbins"] = self._enumbins bkgcube["emin"] = self["emin"].real() bkgcube["emax"] = self["emax"].real() - bkgcube["coordsys"] = self.m_coordsys - bkgcube["proj"] = self.m_proj + bkgcube["coordsys"] = self._coordsys + bkgcube["proj"] = self._proj bkgcube.run() # Notify background cube computation - if self.logTerse(): - self.log("Computed background cube\n") + if self._logTerse(): + self._log("Computed background cube\n") # Store results - self.m_exposure = expcube.expcube().copy() - self.m_psfcube = psfcube.psfcube().copy() - self.m_bckcube = bkgcube.bkgcube().copy() - self.m_stackmodels = bkgcube.models().copy() - - # Return - return - - - 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 input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") - - # Write observation into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Observation(s)") - self.log(str(self.obs)) - self.log("\n") - - # If several observations and binned: prepare stacked irfs - if self.obs.size() > 1 and self.m_enumbins > 0: - self.set_stacked_irf() - - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate pull distribution") - - # Loop over trials - for seed in range(self.m_ntrials): - - # Make a trial - result = self.trial(seed) - - # Write out result immediately - if seed == 0: - file = open(self.m_outfile.url(), 'w') - writer = csv.DictWriter(file, result['colnames']) - headers = {} - for n in result['colnames']: - headers[n] = n - writer.writerow(headers) - else: - file = open(self.m_outfile.url(), 'a') - writer = csv.DictWriter(file, result['colnames']) - writer.writerow(result['values']) - file.close() + self._exposure = expcube.expcube().copy() + self._psfcube = psfcube.psfcube().copy() + self._bckcube = bkgcube.bkgcube().copy() + self._stackmodels = bkgcube.models().copy() # Return return - - def set_obs(self): - """ - Returns an observation container with a set of CTA observations. - - Keywords: + def _trial(self, seed): """ + Compute the pull for a single trial. - # Setup observation definition list - obsdeflist = obsutils.set_obs_patterns(self.m_pattern, - ra=self["ra"].real(), - dec=self["dec"].real(), - offset=self["offset"].real()) - - # Create list of observations - obs = obsutils.set_obs_list(obsdeflist, - tstart=self["tmin"].real(), - duration=self["tmax"].real()-self["tmin"].real(), - deadc=self["deadc"].real(), - emin=self["emin"].real(), - emax=self["emax"].real(), - rad=self["rad"].real(), - irf=self["irf"].string(), - caldb=self["caldb"].string()) - - # Return observation container - return obs - - - def trial(self, seed): + Args: + seed: Random number generator seed. """ - Create the pull for a single trial. - - Parameters: - seed - Random number generator seed - """ - # Write header - if self.logNormal(): - self.log.header2("Trial "+str(seed+1)) + if self._logNormal(): + self._log.header2("Trial "+str(seed+1)) # Simulate events - obs = obsutils.sim(self.obs, - nbins=self.m_enumbins, + obs = obsutils.sim(self._obs, + nbins=self._enumbins, seed=seed, - binsz=self.m_binsz, - npix=self.m_npix, - proj=self.m_proj, - coord=self.m_coordsys, - edisp=self.m_edisp, - log=self.m_log, - debug=self.m_debug, - chatter=self.m_chatter) - - # If stacked, add stacked responses and model - if self.obs.size() > 1 and self.m_enumbins > 0: - if self.m_edisp: - obs[0].response(self.m_exposure, self.m_psfcube, - self.m_edispcube, self.m_bckcube) + binsz=self._binsz, + npix=self._npix, + proj=self._proj, + coord=self._coordsys, + edisp=self._edisp, + log=self._log_clients, + debug=self._debug, + chatter=self._chatter) + + # Set response for a stacked observation + if self._obs.size() > 1 and self._enumbins > 0: + if self._edisp: + obs[0].response(self._exposure, self._psfcube, + self._edispcube, self._bckcube) else: - obs[0].response(self.m_exposure, self.m_psfcube, - self.m_bckcube) - obs.models(self.m_stackmodels) + obs[0].response(self._exposure, self._psfcube, + self._bckcube) + obs.models(self._stackmodels) # Determine number of events in simulation nevents = 0.0 @@ -428,26 +302,26 @@ class cspull(ctools.cscript): nevents += run.events().number() # Write simulation results - if self.logNormal(): - self.log.header3("Simulation") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") + if self._logNormal(): + self._log.header3("Simulation") + self._log.parformat("Number of simulated events") + self._log(nevents) + self._log("\n") # Fit model - if self.m_profile: - models = obs.models() + if self._profile: + models = _obs.models() for i in range(models.size()): model_name = models[i].name() - like = obsutils.cterror(obs, model_name, - log=self.m_log, - debug=self.m_debug, - chatter=self.m_chatter) + like = obsutils.cterror(_obs, model_name, + log=self._log_clients, + debug=self._debug, + chatter=self._chatter) else: - like = obsutils.fit(obs, edisp=self.m_edisp, - log=self.m_log, - debug=self.m_debug, - chatter=self.m_chatter) + like = obsutils.fit(obs, edisp=self._edisp, + log=self._log_clients, + debug=self._debug, + chatter=self._chatter) # Store results logL = like.opt().value() @@ -455,8 +329,8 @@ class cspull(ctools.cscript): models = like.obs().models() # Write result header - if self.logNormal(): - self.log.header3("Pulls") + if self._logNormal(): + self._log.header3("Pulls") # Gather results colnames = [] @@ -484,7 +358,7 @@ class cspull(ctools.cscript): # Compute pull fitted_value = par.value() - real_value = self.obs.models()[i][k].value() + real_value = self._obs.models()[i][k].value() error = par.error() if error != 0.0: pull = (fitted_value - real_value) / error @@ -497,14 +371,14 @@ class cspull(ctools.cscript): values["Unc_"+name] = error # Write result - if self.logNormal(): - self.log.parformat(name) - self.log(pull) - self.log(" (") - self.log(fitted_value) - self.log(" +/- ") - self.log(error) - self.log(")\n") + if self._logNormal(): + self._log.parformat(name) + self._log(pull) + self._log(" (") + self._log(fitted_value) + self._log(" +/- ") + self._log(error) + self._log(")\n") # Bundle together results result = {'colnames': colnames, 'values': values} @@ -513,18 +387,88 @@ class cspull(ctools.cscript): return result + # 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(s)") + self._log(str(self._obs)) + self._log("\n") + + # If several observations and binned: prepare stacked irfs + if self._obs.size() > 1 and self._enumbins > 0: + self._set_stacked_irf() + + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Generate pull distribution") + + # Loop over trials + for seed in range(self._ntrials): + + # Make a trial + result = self._trial(seed) + + # Write out result immediately + if seed == 0: + file = open(self._outfile.url(), 'w') + writer = csv.DictWriter(file, result['colnames']) + headers = {} + for n in result['colnames']: + headers[n] = n + writer.writerow(headers) + else: + file = open(self._outfile.url(), 'a') + writer = csv.DictWriter(file, result['colnames']) + writer.writerow(result['values']) + file.close() + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return + + def models(self, models): + """ + Set model. + """ + # Copy models + self.obs.models(models.copy()) + + # Return + return + + # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates pull distribution for a source model. - """ + # Create instance of application app = cspull(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/csresmap.par b/cscripts/csresmap.par new file mode 100644 index 0000000..533ea76 --- /dev/null +++ b/cscripts/csresmap.par @@ -0,0 +1,57 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csresmap.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csresmap 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" +modcube, f, a, NONE,,, "Input model cube file (generated with ctmodel)" +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?" +outmap, f, a, resmap.fits,,, "Output residual map file" + +# +# Script parameters +#================== +ebinalg, s, h, LOG,LIN|LOG|FILE,, "Binning algorithm" +emin, r, h, 0.1,,, "Lower energy limit (TeV)" +emax, r, h, 100.0,,, "Upper energy limit (TeV)" +enumbins, i, h, 20,,, "Number of energy bins (0=unbinned)" +ebinfile, f, h, NONE,,, "Name of the file containing the energy bin definition" +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)" +algorithm, s, a, SUBDIV,SUB|SUBDIV|SUBDIVSQRT,, "Residual map computation algorithm" + +# +# Standard parameters +#==================== +publish, b, h, no,,, "Publish residual map 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, csresmap.log,,, "Log filename" diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index 036c0e3..9db3ea5 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# Residual map generation script. +# Generates a residual map. # # Copyright (C) 2014-2016 Michael Mayer # @@ -22,57 +22,76 @@ import gammalib import ctools import sys + # ============== # -# cstsdist class # +# csresmap class # # ============== # class csresmap(ctools.cscript): """ + Generates a residual map. + This class implements the creation of a residual map. 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 = "csresmap" - self.version = "1.1.0" - - # Initialise some members - self.obs = None - self.algorithm = "SUB" - self.resmap = None - self.m_modcube = "None" - self.m_outmap = "None" - self.m_use_maps = False - self.m_skip_binning = False - self.m_publish = False - - # Initialise some members + self._name = "csresmap" + self._version = "1.1.0" + + # Initialise class members + self._algorithm = "SUB" + self._resmap = None + self._modcube = "NONE" + self._outmap = "NONE" + self._use_maps = False + self._skip_binning = False + self._publish = False + self._outfile = "" + self._use_maps = False + self._skip_binning = False + self._edisp = False + self._xref = 83.63 + self._yref = 22.01 + self._emin = 0.1 + self._emax = 100.0 + self._enumbins = 0 + self._ebinalg = "LOG" + self._coordsys = "CEL" + self._proj = "CAR" + self._nxpix = 200 + self._nypix = 200 + self._binsz = 0.02 + self._publish = False + 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): - self.obs = argv[0] - argv = argv[1:] + self._obs = argv[0] + argv = argv[1:] else: - self.obs = gammalib.GObservations() - self.obs.clear() - self.outfile = "" - - # Make sure that parfile exists - file = self.parfile() + self._obs = gammalib.GObservations() # 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 @@ -81,271 +100,148 @@ class csresmap(ctools.cscript): """ 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("inobs","f","a","events.fits","","","Input event list, counts cube, or observation definition XML file")) - pars.append(gammalib.GApplicationPar("modcube","f","a","NONE","","","Input model cube file (generated with ctmodel)")) - 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","f","a","NONE","","","Input background cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Input model XML file")) - 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("edisp","b","h","no","","","Apply energy dispersion?")) - pars.append(gammalib.GApplicationPar("outmap","f","a","resmap.fits","","","Output residual map file")) - pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm")) - pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("enumbins","i","h","20","","","Number of energy bins")) - pars.append(gammalib.GApplicationPar("ebinfile","f","h","NONE","","","Name of the file containing the energy bin definition")) - 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("algorithm","s","a","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residual map computation algorithm")) - pars.append(gammalib.GApplicationPar("publish","b","h","no","","","Publish residual map on VO Hub?")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csresmap.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + def _get_parameters(self): """ Get parameters from parfile and setup the observation. """ # Initialise some flags - self.m_use_maps = False - self.m_skip_binning = False + self._use_maps = False + self._skip_binning = False # First check if the inobs parameter is a counts cube - if self.obs.size() == 0 and self["inobs"].filename() != "NONE": + if self._obs.size() == 0 and self["inobs"].filename() != "NONE": filename = gammalib.GFilename(self["inobs"].filename()) if filename.is_fits(): cta = gammalib.GCTAObservation() cta.load(filename) if cta.eventtype() == "CountsCube": - self.m_skip_binning = True + self._skip_binning = True # If we have a counts cube, then ask whether we also have a model - if self.m_skip_binning: - self.m_modcube = self["modcube"].filename() - if self.m_modcube != "NONE": - self.m_use_maps = True + if self._skip_binning: + self._modcube = self["modcube"].filename() + if self._modcube != "NONE": + self._use_maps = True # If not two maps are given, proceed to set up observation - if not self.m_use_maps: + if not self._use_maps: # Set observation if not done before - if self.obs.size() == 0: + if self._obs.size() == 0: self._require_inobs("csresmap.get_parameters()") - self.obs = self._get_observations() + self._obs = self._get_observations() # Check if we have exactly one binned CTA observation - if self.obs.size() == 1: + if self._obs.size() == 1: - if self.obs[0].classname() == "GCTAObservation": - if self.obs[0].eventtype() == "CountsCube": + if self._obs[0].classname() == "GCTAObservation": + if self._obs[0].eventtype() == "CountsCube": # Skip ctbin step later on - self.m_skip_binning = True + self._skip_binning = True # 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()) # Skip query for spatial parameters if a binning is provided in the observation - if not self.m_skip_binning: + if not self._skip_binning: # Read other parameters - self.m_xref = self["xref"].real() - self.m_yref = self["yref"].real() - self.m_emin = self["emin"].real() - self.m_emax = self["emax"].real() - self.m_enumbins = self["enumbins"].integer() - self.m_ebinalg = self["ebinalg"].string() - self.m_coordsys = self["coordsys"].string() - self.m_proj = self["proj"].string() - self.m_nxpix = self["nxpix"].integer() - self.m_nypix = self["nypix"].integer() - self.m_binsz = self["binsz"].real() + self._xref = self["xref"].real() + self._yref = self["yref"].real() + self._emin = self["emin"].real() + self._emax = self["emax"].real() + self._enumbins = self["enumbins"].integer() + self._ebinalg = self["ebinalg"].string() + self._coordsys = self["coordsys"].string() + self._proj = self["proj"].string() + self._nxpix = self["nxpix"].integer() + self._nypix = self["nypix"].integer() + self._binsz = self["binsz"].real() # Read energy dispersion flag - self.m_edisp = self["edisp"].boolean() + self._edisp = self["edisp"].boolean() # Read algorithm - self.m_algorithm = self["algorithm"].string() + self._algorithm = self["algorithm"].string() # Read standard parameters - self.m_publish = self["publish"].boolean() - 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._publish = self["publish"].boolean() + self._chatter = self["chatter"].integer() + self._clobber = self["clobber"].boolean() + self._debug = self["debug"].boolean() # Read ahead output parameters if (self._read_ahead()): - self.m_outmap = self["outmap"].filename() - - # Return - return - - def models(self, models): - """ - Set model. - """ - # Copy models - self.obs.models(models.clone()) - - # Return - return - - def execute(self): - """ - Execute the script. - """ - # Read ahead output parameters - self._read_ahead(True) - - # Run the script - self.run() - - # Save residual map - self.save() + self._outmap = self["outmap"].filename() - # Return - return - - def save(self): - """ - Save residual map. - """ - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Save residual map") - - # Continue only filename and residual map are valid - if self.m_outmap != "None" and self.resmap != None: - - # Log file name - if self.logTerse(): - self.log("Save residual map into \""+self.m_outmap+"\".\n") - - # Save residual map - self.resmap.save(self.m_outmap, self.m_clobber) + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - def publish(self, name=""): - """ - Publish residual map. - """ - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Publish sky map") - - # Continue only if residual map is valid - if self.resmap != None: - - # Set default name is user name is empty - if not name: - user_name = self.name - else: - user_name = name - - # Log map name - if self.logTerse(): - self.log("Publish residual map \""+user_name+"\".\n") - - # Publish map - self.resmap.publish(user_name) - - # 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") # Use input file directly if given - if self.m_use_maps: + if self._use_maps: countmap = gammalib.GSkyMap(self["inobs"].filename()) - modelmap = gammalib.GSkyMap(self.m_modcube) + modelmap = gammalib.GSkyMap(self._modcube) else: - if self.m_skip_binning: - cta_counts_cube = gammalib.GCTAEventCube(self.obs[0].events().clone()) + # ... + if self._skip_binning: + cta_counts_cube = gammalib.GCTAEventCube(self._obs[0].events().clone()) + # ... else: # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate binned map (ctbin)") + if self._logTerse(): + self._log("\n") + self._log.header1("Generate binned map (ctbin)") # Create countsmap - bin = ctools.ctbin(self.obs) - bin["nxpix"].integer(self.m_nxpix) - bin["nypix"].integer(self.m_nypix) - bin["proj"].string(self.m_proj) - bin["coordsys"].string(self.m_coordsys) - bin["xref"].real(self.m_xref) - bin["yref"].real(self.m_yref) - bin["enumbins"].integer(self.m_enumbins) - bin["ebinalg"].string(self.m_ebinalg) - bin["emin"].real(self.m_emin) - bin["emax"].real(self.m_emax) - bin["binsz"].real(self.m_binsz) - bin["chatter"].integer(self.m_chatter) - bin["clobber"].boolean(self.m_clobber) - bin["debug"].boolean(self.m_debug) + bin = ctools.ctbin(self._obs) + bin["nxpix"].integer(self._nxpix) + bin["nypix"].integer(self._nypix) + bin["proj"].string(self._proj) + bin["coordsys"].string(self._coordsys) + bin["xref"].real(self._xref) + bin["yref"].real(self._yref) + bin["enumbins"].integer(self._enumbins) + bin["ebinalg"].string(self._ebinalg) + bin["emin"].real(self._emin) + bin["emax"].real(self._emax) + bin["binsz"].real(self._binsz) + bin["chatter"].integer(self._chatter) + bin["clobber"].boolean(self._clobber) + bin["debug"].boolean(self._debug) bin.run() # Retrieve counts cube @@ -355,17 +251,17 @@ class csresmap(ctools.cscript): countmap = cta_counts_cube.map() # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate model map (ctmodel)") + if self._logTerse(): + self._log("\n") + self._log.header1("Generate model map (ctmodel)") # Create model map - model = ctools.ctmodel(self.obs) + model = ctools.ctmodel(self._obs) model.cube(cta_counts_cube) - model["chatter"].integer(self.m_chatter) - model["clobber"].boolean(self.m_clobber) - model["debug"].boolean(self.m_debug) - model["edisp"].boolean(self.m_edisp) + model["chatter"].integer(self._chatter) + model["clobber"].boolean(self._clobber) + model["debug"].boolean(self._debug) + model["edisp"].boolean(self._edisp) model.run() # Get model map into GSkyMap object @@ -376,61 +272,141 @@ class csresmap(ctools.cscript): # to be a SWIG feature as SWIG creates a new object when # calling bin.cube() #residualmap = bin.cube().map() - self.resmap = countmap.copy() - self.resmap.stack_maps() + self._resmap = countmap.copy() + self._resmap.stack_maps() modelmap.stack_maps() # Continue calculations depending on given algorithm - if self.m_algorithm == "SUB": + if self._algorithm == "SUB": # Subtract maps - self.resmap -= modelmap + self._resmap -= modelmap - elif self.m_algorithm == "SUBDIV": + elif self._algorithm == "SUBDIV": # Subtract and divide by model map - self.resmap -= modelmap - self.resmap /= modelmap + self._resmap -= modelmap + self._resmap /= modelmap #for pixel in modelmap: # if pixel != 0.0: # pixel = 1.0/pixel - #self.resmap *= modelmap + #self._resmap *= modelmap - elif self.m_algorithm == "SUBDIVSQRT": + elif self._algorithm == "SUBDIVSQRT": # subtract and divide by sqrt of model map - self.resmap -= modelmap - self.resmap /= modelmap.sqrt() + self._resmap -= modelmap + self._resmap /= modelmap.sqrt() #for pixel in modelmap: # if pixel != 0.0: # pixel = 1.0/math.sqrt(pixel) - #self.resmap *= modelmap + #self._resmap *= modelmap else: # Raise error if algorithm is unkown - raise TypeError("Algorithm \""+self.m_algorithm+"\" not known") + raise TypeError("Algorithm \""+self._algorithm+"\" not known") # Optionally publish map - if self.m_publish: + if self._publish: self.publish() # Return return + def execute(self): + """ + Execute the script. + """ + # Read ahead output parameters + self._read_ahead(True) + + # Run the script + self.run() + + # Save residual map + self.save() + + # Return + return + + def save(self): + """ + Save residual map. + """ + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save residual map") + + # Get outmap parameter + self._outmap = self["outmap"].filename() + + # Continue only filename and residual map are valid + if self._outmap != "NONE" and self._resmap != None: + + # Log file name + if self._logTerse(): + self._log("Save residual map into \""+self._outmap+"\".\n") + + # Save residual map + self._resmap.save(self._outmap, self._clobber) + + # Return + return + + def publish(self, name=""): + """ + Publish residual map. + + Kwargs: + name: Name of residual map. + """ + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Publish residual map") + + # Continue only if residual map is valid + if self._resmap != None: + + # Set default name is user name is empty + if not name: + user_name = self.name + else: + user_name = name + + # Log map name + if self._logTerse(): + self._log("Publish residual map \""+user_name+"\".\n") + + # Publish map + self._resmap.publish(user_name) + + # Return + return + + def models(self, models): + """ + Set model. + """ + # Copy models + self._obs.models(models.clone()) + + # Return + return + # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates residual count map. - """ + # Create instance of application app = csresmap(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/cssens.par b/cscripts/cssens.par new file mode 100644 index 0000000..90b544e --- /dev/null +++ b/cscripts/cssens.par @@ -0,0 +1,55 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cssens.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the cssens script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inobs, f, h, NONE,,, "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" +caldb, s, a, prod2,,, "Calibration database" +irf, s, a, South_0.5h,,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" +deadc, r, h, 0.95,0,1, "Deadtime correction factor" +outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file" + +# +# Script parameters +#================== +offset, r, h, 0.0,,, "Source offset angle (deg)" +duration, r, a, 180000.0,,, "Effective exposure time (s)" +rad, r, a, 5.0,,, "Radius of ROI (deg)" +emin, r, a, 0.020,,, "Lower energy limit (TeV)" +emax, r, a, 200.0,,, "Upper energy limit (TeV)" +bins, i, a, 21,,, "Number of energy bins for differential sensitivity computation" +enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)" +npix, i, h, 200,,, "Number of pixels for binned analysis" +binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)" +type, s, h, Differential,Differential|Integral,, "Sensitivity type" +sigma, r, h, 5.0,,, "Significance threshold (Gaussian sigma)" +max_iter, i, h, 50,,, "Maximum number of iterations" +num_avg, i, h, 3,,, "Number of iterations for sliding average" + +# +# 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, cspull.log,,, "Log filename" + + diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 378c5dc..7d5eef0 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -31,6 +31,8 @@ import math # ============ # class cssens(ctools.cscript): """ + Computes the CTA sensitivity. + This class computes the CTA sensitivity for a number of energy bins using ctlike. Spectra are fitted in narrow energy bins to simulated data, and the flux level is determined that leads to a particular significance. @@ -38,34 +40,48 @@ class cssens(ctools.cscript): The significance is determined using the Test Statistic, defined as twice the likelihood difference between fitting with and without the test source. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "cssens" - self.version = "1.0.0" - - # Initialise some parameters - self.obs = None - self.m_ra = None - self.m_dec = None - - # Make sure that parfile exists - file = self.parfile() - + self._name = "cssens" + self._version = "1.1.0" + + # Initialise class members + self._obs = gammalib.GObservations() + self._ebounds = gammalib.GEbounds() + self._srcname = "" + self._outfile = "" + self._ra = None + self._dec = None + self._edisp = False + self._emin = 0.020 + self._emax = 200.0 + self._bins = 21 + self._enumbins = 0 + self._npix = 200 + self._binsz = 0.05 + self._type = "Differential" + self._ts_thres = 25.0 + self._max_iter = 50 + self._num_avg = 3 + self._log_clients = False + self._debug = False + # 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.logFileOpen() - self.log_header() - self.log.date(True) + self._log_header() + self._log.date(True) # Return return @@ -74,216 +90,79 @@ class cssens(ctools.cscript): """ 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 - print("Parfile "+parfile+" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","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("caldb","s","a","prod2","","","Calibration database")) - pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) - pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor")) - pars.append(gammalib.GApplicationPar("outfile","f","h","sensitivity.dat","","","Output sensitivity ASCII file")) - pars.append(gammalib.GApplicationPar("offset","r","h","0.0","","","Source offset angle (deg)")) - pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) - pars.append(gammalib.GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)")) - pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) - pars.append(gammalib.GApplicationPar("emin","r","a","0.020","","","Lower energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","a","200.0","","","Upper energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("bins","i","a","21","","","Number of energy bins for differential sensitivity computation")) - pars.append(gammalib.GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)")) - pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) - pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) - pars.append(gammalib.GApplicationPar("type","s","h","Differential","Differential|Integral","","Sensitivity type")) - pars.append(gammalib.GApplicationPar("sigma","r","h","5.0","","","Significance threshold")) - pars.append(gammalib.GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations")) - pars.append(gammalib.GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","cssens.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + def _get_parameters(self): """ Get user parameters from parfile. """ # Set observation if not done before - if self.obs == None or self.obs.size() == 0: - self.obs = self.set_obs() + if self._obs == None or self._obs.size() == 0: + self._obs = self._set_obs() # 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() # Read further parameters - self.m_outfile = self["outfile"].filename() - self.m_emin = self["emin"].real() - self.m_emax = self["emax"].real() - self.m_bins = self["bins"].integer() + self._outfile = self["outfile"].filename() + self._emin = self["emin"].real() + self._emax = self["emax"].real() + self._bins = self["bins"].integer() # Read parameters for binned if requested - self.m_enumbins = self["enumbins"].integer() - if not self.m_enumbins == 0: - self.m_npix = self["npix"].integer() - self.m_binsz = self["binsz"].real() - else: - # Set dummy values (required by obsutils) - self.m_npix = 0 - self.m_binsz = 0.0 + self._enumbins = self["enumbins"].integer() + if not self._enumbins == 0: + self._npix = self["npix"].integer() + self._binsz = self["binsz"].real() # Read remaining parameters - self.m_edisp = self["edisp"].boolean() - self.m_ts_thres = self["sigma"].real()*self["sigma"].real() - self.m_max_iter = self["max_iter"].integer() - self.m_num_avg = self["num_avg"].integer() - self.m_type = self["type"].string() + self._edisp = self["edisp"].boolean() + self._ts_thres = self["sigma"].real()*self["sigma"].real() + self._max_iter = self["max_iter"].integer() + self._num_avg = self["num_avg"].integer() + self._type = self["type"].string() # Set some fixed parameters - self.m_log = False # Logging in client tools - self.m_debug = self["debug"].boolean() # Debugging in client tools + self._debug = self["debug"].boolean() # Debugging in client tools # Derive some parameters - self.m_ebounds = gammalib.GEbounds(self.m_bins, - gammalib.GEnergy(self.m_emin, "TeV"), - gammalib.GEnergy(self.m_emax, "TeV")) + self._ebounds = gammalib.GEbounds(self._bins, + gammalib.GEnergy(self._emin, "TeV"), + gammalib.GEnergy(self._emax, "TeV")) - # Return - return - - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - def run(self): - """ - Run the script. + def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): """ - # Switch screen logging on in debug mode - 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") - - # Initialise script - colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', - 'energy_flux', 'sensitivity'] - results = [] - - # Initialise models - full_model, bkg_model = self.set_models() - - # Write models into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Models") - self.log.header2("Background model") - self.log(str(bkg_model)) - self.log("\n\n") - self.log.header2("Full model") - self.log(str(full_model)) - self.log("\n") - - # Write heder - if self.logTerse(): - self.log("\n") - self.log.header1("Sensitivity determination") - self.log.parformat("Type") - self.log(self.m_type) - self.log("\n") - - # Loop over energy bins - for ieng in range(self.m_ebounds.size()): - - # Set energies - if self.m_type == "Differential": - emin = self.m_ebounds.emin(ieng) - emax = self.m_ebounds.emax(ieng) - elif self.m_type == "Integral": - emin = self.m_ebounds.emin(ieng) - emax = self.m_ebounds.emax() - else: - msg = "Invalid sensitivity type \""+self.m_type+"\" encountered."+ \ - " Either use \"Differential\" or \"Integral\"." - raise RuntimeError(msg) + Set an observation container. - # Determine sensitivity - result = self.get_sensitivity(self.obs, emin, emax, - bkg_model, full_model) + Kwargs: + lpnt: Galactic longitude of pointing [deg] (default: 0.0) + bpnt: Galactic latitude of pointing [deg] (default: 0.0) + emin: Minimum energy [TeV] (default: 0.1) + emax: Maximum energy [TeV] (default: 100.0) - # Write results - if ieng == 0: - f = open(self.m_outfile.url(), 'w') - writer = csv.DictWriter(f, colnames) - headers = {} - for n in colnames: - headers[n] = n - writer.writerow(headers) - #writer.writerow(dict((_,_) for _ in colnames)) - writer.writerow(result) - f.close() - else: - f = open(self.m_outfile.url(), 'a') - writer = csv.DictWriter(f, colnames) - writer.writerow(result) - f.close() - - # Append results - results.append(result) - - # Return - return - - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): + Returns: + Observation container. """ - Returns an observation container. - - Keywords: - lpnt - Galactic longitude of pointing [deg] (default: 0.0) - bpnt - Galactic latitude of pointing [deg] (default: 0.0) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - """ - # If an observation was provided on input then load it from XML file + # If an observation was provided on input then load it from XML + # file filename = self["inobs"].filename() if filename != "NONE" and filename != "": - obs = self.get_observations() + obs = self._get_observations() # ... otherwise allocate a single observation else: @@ -311,56 +190,46 @@ class cssens(ctools.cscript): obs.append(run) # Set source position - offset = self["offset"].real() + offset = self["offset"].real() pntdir.lb_deg(lpnt, bpnt+offset) - self.m_ra = pntdir.ra_deg() - self.m_dec = pntdir.dec_deg() + self._ra = pntdir.ra_deg() + self._dec = pntdir.dec_deg() # Return observation container return obs - def set_obs_ebounds(self, emin, emax): - """ - Set energy boundaries for observation in container. - - Parameters: - emin - Minimum energy - emax - Maximum energy + def _set_models(self, fitspat=False, fitspec=False): """ - # Loop over all observations in container - for obs in self.obs: - - # Set energy boundaries - ebounds = gammalib.GEbounds(emin, emax) - obs.events().ebounds(ebounds) + Set full model and background model. - # Return - return + Kwargs: + fitspat: Fit spatial parameter (default: False). + fitspec: Fit spectral parameters (default: False). - def set_models(self, fitpos=False, fitspec=False): - """ - Set full and background model. + Returns: + Tuple containing full model and background model. """ # Retrieve full model from observation container - full_model = self.obs.models().copy() + full_model = self._obs.models().copy() # Get source model - model = full_model[self.m_srcname] + model = full_model[self._srcname] # Check that model has a Prefactor if not model.has_par("Prefactor"): - msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \ - " Only spectral models with a \"Prefactor\" parameter are supported." + msg = "Model \""+self._srcname+"\" has no parameter "+\ + "\"Prefactor\". Only spectral models with a "+\ + "\"Prefactor\" parameter are supported." raise RuntimeError(msg) # Set source position - if self.m_ra != None and self.m_dec != None: + if self._ra != None and self._dec != None: if model.has_par("RA") and model.has_par("DEC"): - model["RA"].value(self.m_ra) - model["DEC"].value(self.m_dec) + model["RA"].value(self._ra) + model["DEC"].value(self._dec) # Fit or fix spatial parameters - if fitpos: + if fitspat: if model.has_par("RA"): model["RA"].free() if model.has_par("DEC"): @@ -397,27 +266,71 @@ class cssens(ctools.cscript): # Create background model bkg_model = full_model.copy() - bkg_model.remove(self.m_srcname) + bkg_model.remove(self._srcname) # Return models return full_model, bkg_model - def get_sensitivity(self, obs, emin, emax, bkg_model, full_model): + def _set_obs_ebounds(self, emin, emax): + """ + Set energy boundaries for observation container. + + Sets the energy boundaries for all observations in the observation + container. + + Args: + emin: Minimum energy + emax: Maximum energy + """ + # Loop over all observations in container + for obs in self._obs: + + # Set energy boundaries + ebounds = gammalib.GEbounds(emin, emax) + obs.events().ebounds(ebounds) + + # Return + return + + def _get_crab_flux(self, emin, emax): + """ + Return Crab photon flux in a given energy interval. + + Args: + emin: Minimum energy + emax: Maximum energy + + Returns: + Crab photon flux in specified energy interval. + """ + # Set Crab TeV spectral model based on a power law + crab = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, + gammalib.GEnergy(0.3, "TeV")) + + # Determine flux + flux = crab.flux(emin, emax) + + # Return flux + return flux + + def _get_sensitivity(self, emin, emax, bkg_model, full_model): """ Determine sensitivity for given observations. - Parameters: - obs - Observation container - emin - Minimum energy for fitting and flux computation - emax - Maximum energy for fitting and flux computation - bkg_model - Background model - full_model - Source model + Args: + emin: Minimum energy for fitting and flux computation + emax: Maximum energy for fitting and flux computation + bkg_model: Background model + full_model: Source model + + Returns: + Result dictionary """ # Set TeV->erg conversion factor tev2erg = 1.6021764 # Set energy boundaries - self.set_obs_ebounds(emin, emax) + self._set_obs_ebounds(emin, emax) # Determine energy boundaries from first observation in the container loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) @@ -426,23 +339,23 @@ class cssens(ctools.cscript): # Compute Crab unit (this is the factor with which the Prefactor needs # to be multiplied to get 1 Crab - crab_flux = self.get_crab_flux(emin, emax) - src_flux = full_model[self.m_srcname].spectral().flux(emin, emax) + crab_flux = self._get_crab_flux(emin, emax) + src_flux = full_model[self._srcname].spectral().flux(emin, emax) crab_unit = crab_flux/src_flux # Write header - if self.logTerse(): - self.log("\n") - self.log.header2("Energies: "+str(emin)+" - "+str(emax)) - self.log.parformat("Crab flux") - self.log(crab_flux) - self.log(" ph/cm2/s\n") - self.log.parformat("Source model flux") - self.log(src_flux) - self.log(" ph/cm2/s\n") - self.log.parformat("Crab unit factor") - self.log(crab_unit) - self.log("\n") + if self._logTerse(): + self._log("\n") + self._log.header2("Energies: "+str(emin)+" - "+str(emax)) + self._log.parformat("Crab flux") + self._log(crab_flux) + self._log(" ph/cm2/s\n") + self._log.parformat("Source model flux") + self._log(src_flux) + self._log(" ph/cm2/s\n") + self._log.parformat("Crab unit factor") + self._log(crab_unit) + self._log("\n") # Initialise loop crab_flux_value = [] @@ -459,20 +372,23 @@ class cssens(ctools.cscript): iter += 1 # Write header - if self.logExplicit(): - self.log.header2("Iteration "+str(iter)) + if self._logExplicit(): + self._log.header2("Iteration "+str(iter)) # Set source model. crab_prefactor is the Prefactor that # corresponds to 1 Crab src_model = full_model.copy() - crab_prefactor = src_model[self.m_srcname]['Prefactor'].value() * crab_unit - src_model[self.m_srcname]['Prefactor'].value(crab_prefactor * test_crab_flux) - obs.models(src_model) + crab_prefactor = src_model[self._srcname]['Prefactor'].value() * \ + crab_unit + src_model[self._srcname]['Prefactor'].value(crab_prefactor * \ + test_crab_flux) + self._obs.models(src_model) # Simulate events - sim = obsutils.sim(obs, nbins=self.m_enumbins, seed=iter, - binsz=self.m_binsz, npix=self.m_npix, - log=self.m_log, debug=self.m_debug, edisp=self.m_edisp) + sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iter, + binsz=self._binsz, npix=self._npix, + log=self._log_clients, debug=self._debug, + edisp=self._edisp) # Determine number of events in simulation nevents = 0.0 @@ -480,15 +396,16 @@ class cssens(ctools.cscript): nevents += run.events().number() # Write simulation results - if self.logExplicit(): - self.log.header3("Simulation") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") + if self._logExplicit(): + self._log.header3("Simulation") + self._log.parformat("Number of simulated events") + self._log(nevents) + self._log("\n") # Fit background only sim.models(bkg_model) - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug, edisp=self.m_edisp) + like = obsutils.fit(sim, log=self._log_clients, + debug=self._debug, edisp=self._edisp) result_bgm = like.obs().models().copy() LogL_bgm = like.opt().value() npred_bgm = like.obs().npred() @@ -497,36 +414,31 @@ class cssens(ctools.cscript): quality_bgm = npred_bgm-nevents # Write background fit results - if self.logExplicit(): - self.log.header3("Background model fit") - self.log.parformat("log likelihood") - self.log(LogL_bgm) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_bgm) - self.log("\n") - self.log.parformat("Fit quality") - self.log(quality_bgm) - self.log("\n") - - # Start over if the fit quality was bad - #if abs(quality_bgm) > 3.0: - # if self.logExplicit(): - # self.log("Fit quality outside required range. Start over.\n") - # continue + if self._logExplicit(): + self._log.header3("Background model fit") + self._log.parformat("log likelihood") + self._log(LogL_bgm) + self._log("\n") + self._log.parformat("Number of predicted events") + self._log(npred_bgm) + self._log("\n") + self._log.parformat("Fit quality") + self._log(quality_bgm) + self._log("\n") # Write model fit results - if self.logExplicit(): + if self._logExplicit(): for model in result_bgm: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") + self._log(str(par)+"\n") # Fit background and test source sim.models(src_model) - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug, edisp=self.m_edisp) + like = obsutils.fit(sim, log=self._log_clients, + debug=self._debug, edisp=self._edisp) result_all = like.obs().models().copy() LogL_all = like.opt().value() npred_all = like.obs().npred() @@ -536,56 +448,52 @@ class cssens(ctools.cscript): quality_all = npred_all-nevents # Write background and test source fit results - if self.logExplicit(): - self.log.header3("Background and test source model fit") - self.log.parformat("Test statistics") - self.log(ts) - self.log("\n") - self.log.parformat("log likelihood") - self.log(LogL_all) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_all) - self.log("\n") - self.log.parformat("Fit quality") - self.log(quality_all) - self.log("\n") + if self._logExplicit(): + self._log.header3("Background and test source model fit") + self._log.parformat("Test statistics") + self._log(ts) + self._log("\n") + self._log.parformat("log likelihood") + self._log(LogL_all) + self._log("\n") + self._log.parformat("Number of predicted events") + self._log(npred_all) + self._log("\n") + self._log.parformat("Fit quality") + self._log(quality_all) + self._log("\n") # for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") - - # Start over if the fit quality was bad - #if abs(quality_all) > 3.0: - # if self.logExplicit(): - # self.log("Fit quality outside required range. Start over.\n") - # continue + self._log(str(par)+"\n") # Start over if TS was non-positive if ts <= 0.0: - if self.logExplicit(): - self.log("Non positive TS. Start over.\n") + if self._logExplicit(): + self._log("Non positive TS. Start over.\n") continue # Get fitted Crab, photon and energy fluxes - crab_flux = result_all[self.m_srcname]['Prefactor'].value() / crab_prefactor - crab_flux_err = result_all[self.m_srcname]['Prefactor'].error() / crab_prefactor - photon_flux = result_all[self.m_srcname].spectral().flux(emin, emax) - energy_flux = result_all[self.m_srcname].spectral().eflux(emin, emax) + crab_flux = result_all[self._srcname]['Prefactor'].value() / \ + crab_prefactor + crab_flux_err = result_all[self._srcname]['Prefactor'].error() / \ + crab_prefactor + photon_flux = result_all[self._srcname].spectral().flux(emin, emax) + energy_flux = result_all[self._srcname].spectral().eflux(emin, emax) # Compute differential sensitivity in unit erg/cm2/s energy = gammalib.GEnergy(e_mean, "TeV") time = gammalib.GTime() - sensitivity = result_all[self.m_srcname].spectral().eval(energy, time) * \ + sensitivity = result_all[self._srcname].spectral().eval(energy, time) * \ e_mean*erg_mean * 1.0e6 # Compute flux correction factor based on average TS correct = 1.0 if ts > 0: - correct = math.sqrt(self.m_ts_thres/ts) + correct = math.sqrt(self._ts_thres/ts) # Compute extrapolated fluxes crab_flux = correct * crab_flux @@ -598,41 +506,41 @@ class cssens(ctools.cscript): sensitivity_value.append(sensitivity) # Write background and test source fit results - if self.logExplicit(): - self.log.parformat("Photon flux") - self.log(photon_flux) - self.log(" ph/cm2/s\n") - self.log.parformat("Energy flux") - self.log(energy_flux) - self.log(" erg/cm2/s\n") - self.log.parformat("Crab flux") - self.log(crab_flux*1000.0) - self.log(" mCrab\n") - self.log.parformat("Differential sensitivity") - self.log(sensitivity) - self.log(" erg/cm2/s\n") + if self._logExplicit(): + self._log.parformat("Photon flux") + self._log(photon_flux) + self._log(" ph/cm2/s\n") + self._log.parformat("Energy flux") + self._log(energy_flux) + self._log(" erg/cm2/s\n") + self._log.parformat("Crab flux") + self._log(crab_flux*1000.0) + self._log(" mCrab\n") + self._log.parformat("Differential sensitivity") + self._log(sensitivity) + self._log(" erg/cm2/s\n") for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") - elif self.logTerse(): - self.log.parformat("Iteration "+str(iter)) - self.log("TS=") - self.log(ts) - self.log(" ") - self.log("corr=") - self.log(correct) - self.log(" ") - self.log(photon_flux) - self.log(" ph/cm2/s = ") - self.log(energy_flux) - self.log(" erg/cm2/s = ") - self.log(crab_flux*1000.0) - self.log(" mCrab = ") - self.log(sensitivity) - self.log(" erg/cm2/s\n") + self._log(str(par)+"\n") + elif self._logTerse(): + self._log.parformat("Iteration "+str(iter)) + self._log("TS=") + self._log(ts) + self._log(" ") + self._log("corr=") + self._log(correct) + self._log(" ") + self._log(photon_flux) + self._log(" ph/cm2/s = ") + self._log(energy_flux) + self._log(" erg/cm2/s = ") + self._log(crab_flux*1000.0) + self._log(" mCrab = ") + self._log(sensitivity) + self._log(" erg/cm2/s\n") # Compute sliding average of extrapolated fitted prefactor, # photon and energy flux. This damps out fluctuations and @@ -642,7 +550,7 @@ class cssens(ctools.cscript): energy_flux = 0.0 sensitivity = 0.0 num = 0.0 - for k in range(self.m_num_avg): + for k in range(self._num_avg): inx = len(crab_flux_value) - k - 1 if inx >= 0: crab_flux += crab_flux_value[inx] @@ -656,7 +564,7 @@ class cssens(ctools.cscript): sensitivity /= num # Compare average flux to last average - if iter > self.m_num_avg: + if iter > self._num_avg: if test_crab_flux > 0: ratio = crab_flux/test_crab_flux @@ -665,70 +573,71 @@ class cssens(ctools.cscript): # 2. The flux correction factor is small if ratio >= 0.99 and ratio <= 1.01 and \ correct >= 0.9 and correct <= 1.1: - if self.logTerse(): - self.log(" Converged ("+str(ratio)+")\n") + if self._logTerse(): + self._log(" Converged ("+str(ratio)+")\n") break else: - if self.logTerse(): - self.log(" Flux is zero.\n") + if self._logTerse(): + self._log(" Flux is zero.\n") break # Use average for next iteration test_crab_flux = crab_flux # Exit loop if number of trials exhausted - if (iter >= self.m_max_iter): - if self.logTerse(): - self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n") + if (iter >= self._max_iter): + if self._logTerse(): + self._log(" Test ended after "+str(self._max_iter)+ + " iterations.\n") break # Write fit results - if self.logTerse(): - self.log.header3("Fit results") - self.log.parformat("Test statistics") - self.log(ts) - self.log("\n") - self.log.parformat("Photon flux") - self.log(photon_flux) - self.log(" ph/cm2/s\n") - self.log.parformat("Energy flux") - self.log(energy_flux) - self.log(" erg/cm2/s\n") - self.log.parformat("Crab flux") - self.log(crab_flux*1000.0) - self.log(" mCrab\n") - self.log.parformat("Differential sensitivity") - self.log(sensitivity) - self.log(" erg/cm2/s\n") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") - self.log.header3("Background and test source model fitting") - self.log.parformat("log likelihood") - self.log(LogL_all) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_all) - self.log("\n") + if self._logTerse(): + self._log.header3("Fit results") + self._log.parformat("Test statistics") + self._log(ts) + self._log("\n") + self._log.parformat("Photon flux") + self._log(photon_flux) + self._log(" ph/cm2/s\n") + self._log.parformat("Energy flux") + self._log(energy_flux) + self._log(" erg/cm2/s\n") + self._log.parformat("Crab flux") + self._log(crab_flux*1000.0) + self._log(" mCrab\n") + self._log.parformat("Differential sensitivity") + self._log(sensitivity) + self._log(" erg/cm2/s\n") + self._log.parformat("Number of simulated events") + self._log(nevents) + self._log("\n") + self._log.header3("Background and test source model fitting") + self._log.parformat("log likelihood") + self._log(LogL_all) + self._log("\n") + self._log.parformat("Number of predicted events") + self._log(npred_all) + self._log("\n") for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") - self.log.header3("Background model fit") - self.log.parformat("log likelihood") - self.log(LogL_bgm) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_bgm) - self.log("\n") + self._log(str(par)+"\n") + self._log.header3("Background model fit") + self._log.parformat("log likelihood") + self._log(LogL_bgm) + self._log("\n") + self._log.parformat("Number of predicted events") + self._log(npred_bgm) + self._log("\n") for model in result_bgm: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") + self._log(str(par)+"\n") # Store result result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ @@ -739,29 +648,110 @@ class cssens(ctools.cscript): # Return result return result - def get_crab_flux(self, emin, emax): + + # Public methods + def run(self): """ - Determine the Crab photon flux in a given energy interval. + Run the script. """ - # Set Crab spectral model - crab = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, gammalib.GEnergy(0.3, "TeV")) + # Switch screen logging on in debug mode + if self._logDebug(): + self._log.cout(True) - # Determine flux - flux = crab.flux(emin, emax) + # Get parameters + self._get_parameters() + + # Initialise script + colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', + 'energy_flux', 'sensitivity'] + results = [] + + # Initialise models + full_model, bkg_model = self._set_models() + + # Write models into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Models") + self._log.header2("Background model") + self._log(str(bkg_model)) + self._log("\n\n") + self._log.header2("Full model") + self._log(str(full_model)) + self._log("\n") + + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Sensitivity determination") + self._log.parformat("Type") + self._log(self._type) + self._log("\n") + + # Loop over energy bins + for ieng in range(self._ebounds.size()): + + # Set energies + if self._type == "Differential": + emin = self._ebounds.emin(ieng) + emax = self._ebounds.emax(ieng) + elif self._type == "Integral": + emin = self._ebounds.emin(ieng) + emax = self._ebounds.emax() + else: + msg = "Invalid sensitivity type \""+self._type+"\" "+\ + "encountered. Either specify \"Differential\" or "+\ + "\"Integral\"." + raise RuntimeError(msg) + + # Determine sensitivity + result = self._get_sensitivity(emin, emax, + bkg_model, full_model) + + # Write results + if ieng == 0: + f = open(self._outfile.url(), 'w') + writer = csv.DictWriter(f, colnames) + headers = {} + for n in colnames: + headers[n] = n + writer.writerow(headers) + writer.writerow(result) + f.close() + else: + f = open(self._outfile.url(), 'a') + writer = csv.DictWriter(f, colnames) + writer.writerow(result) + f.close() + + # Append results + results.append(result) + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return - # Return flux - return flux # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - CTA ctlike-based sensitivity estimator. - """ + # Create instance of application app = cssens(sys.argv) + # Open logfile + app._logFileOpen() + # Run application app.execute() diff --git a/cscripts/csspec.par b/cscripts/csspec.par new file mode 100644 index 0000000..05b7042 --- /dev/null +++ b/cscripts/csspec.par @@ -0,0 +1,52 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csspec.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the csspec 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, spectrum.fits,,, "Output spectrum file" + +# +# Script parameters +#================== +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (TeV)" +enumbins, i, a, 20,,, "Number of energy bins (0=unbinned)" +ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm" +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" +calc_ts, b, h, yes,,, "Compute TS value in each bin" +calc_ulim, b, h, yes,,, "Compute upper limit in each bin" +fix_srcs, b, h, yes,,, "Fix other skymodel parameters" +fix_bkg, b, h, no,,, "Fix background 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, csspec.log,,, "Log filename" diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 98263a3..0b82fd4 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# Spectral points generation script. +# Generates a spectrum. # # Copyright (C) 2014-2016 Michael Mayer # @@ -28,47 +28,57 @@ import sys # ============ # class csspec(ctools.cscript): """ + Generates a spectrum. + This class implements the creation of spectral points. 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 = "csspec" - self.version = "1.1.0" - - # Initialise some members - self.obs = gammalib.GObservations() - self.fits = None + self._name = "csspec" + self._version = "1.1.0" # Initialise some members + self._ebounds = gammalib.GEbounds() + self._fits = None + self._binned_mode = False + self._srcname = "" + self._edisp = False + self._outfile = "NONE" + self._calc_ulimit = True + self._calc_ts = True + self._fix_bkg = False + self._fix_srcs = True + 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): - 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 = "" - - # Make sure that parfile exists - file = self.parfile() + self._obs = gammalib.GObservations() # 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 @@ -77,121 +87,81 @@ class csspec(ctools.cscript): """ 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("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","f","a","NONE","","","Input background cube file (only needed for stacked analysis)")) - 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("edisp","b","h","no","","","Apply energy dispersion?")) - pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output spectrum file")) - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit for spectral points (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit for spectral points (TeV)")) - pars.append(gammalib.GApplicationPar("enumbins","i","a","20","1","10000","Number of spectral points")) - pars.append(gammalib.GApplicationPar("ebinalg","s","a","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins")) - pars.append(gammalib.GApplicationPar("ebinfile","f","a","NONE","","", "Name of the file containing the energy bin definition")) - pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin")) - pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin")) - pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters")) - pars.append(gammalib.GApplicationPar("fix_bkg","b","h","no","yes|no","","Fix background parameters")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","csspec.log","","","Log filename")) - pars.save(parfile) - # Return - return - - def get_parameters(self): + # Private methods + 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("csspec::get_parameters()") - self.obs = self._get_observations() + self._obs = self._get_observations() # Check if we have one binned cta observation, i.e. if we are in binned mode - self.m_binned_mode = False - if self.obs.size() == 1: - if self.obs[0].classname() == "GCTAObservation": - if self.obs[0].eventtype() == "CountsCube": - self.m_binned_mode = True + self._binned_mode = False + if self._obs.size() == 1: + if self._obs[0].classname() == "GCTAObservation": + if self._obs[0].eventtype() == "CountsCube": + self._binned_mode = True # 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() # Set ebounds - self.set_ebounds() + self._set_ebounds() # Get edisp flag - self.m_edisp = self["edisp"].boolean() - - # Read other parameters - self.m_outfile = self["outfile"].filename() + self._edisp = self["edisp"].boolean() # 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() + + # Read ahead output parameters + if (self._read_ahead()): + self._outfile = self["outfile"].filename() + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - def set_ebounds(self): + def _set_ebounds(self): """ Set energy boundaries. """ # Get ebounds - if self.m_binned_mode: + if self._binned_mode: # Extract cube ebounds - cube_ebounds = self.obs[0].events().ebounds() + cube_ebounds = self._obs[0].events().ebounds() # Read user parameters - self.m_enumbins = self["enumbins"].integer() - self.m_emin = self["emin"].real() - 1e-6 # Rounding tolerance - self.m_emax = self["emax"].real() + 1e-6 # Rounding tolerance + enumbins = self["enumbins"].integer() + emin_val = self["emin"].real() - 1e-6 # Rounding tolerance + emax_val = self["emax"].real() + 1e-6 # Rounding tolerance # Set energy range - emin = gammalib.GEnergy(self.m_emin, "TeV") - emax = gammalib.GEnergy(self.m_emax, "TeV") + emin = gammalib.GEnergy(emin_val, "TeV") + emax = gammalib.GEnergy(emax_val, "TeV") # Determine all cube layers that overlap with the [emin,emax] # energy range. @@ -210,14 +180,14 @@ class csspec(ctools.cscript): raise RuntimeError(msg) # Determine number of layers to use for each spectral bin - n_layers = int(len(use_layers) / self.m_enumbins) + n_layers = int(len(use_layers) / enumbins) # Don't allow zero layers if n_layers == 0: n_layers = 1 # Create energy boundaries - self.m_ebounds = gammalib.GEbounds() + self._ebounds = gammalib.GEbounds() i_start = 0 i_stop = n_layers - 1 while i_start < len(use_layers): @@ -225,51 +195,13 @@ class csspec(ctools.cscript): i_stop = len(use_layers) - 1 ebinmin = cube_ebounds.emin(use_layers[i_start]) ebinmax = cube_ebounds.emax(use_layers[i_stop]) - self.m_ebounds.append(ebinmin, ebinmax) + self._ebounds.append(ebinmin, ebinmax) i_start += n_layers i_stop += n_layers - #for i in range(len(use_layers))[::n_layers]: - # i_start = i - # i_stop = i + n_layers - 1 - # if i_stop >= len(use_layers): - # i_stop = len(use_layers) - 1 - # ebinmin = cube_ebounds.emin(use_layers[i_start]) - # ebinmax = cube_ebounds.emax(use_layers[i_stop]) - # self.m_ebounds.append(ebinmin, ebinmax) # Unbinned mode else: - self.m_ebounds = self._create_ebounds() - - # Return - return - - def models(self, models): - """ - Set model. - """ - # Copy models - self.obs.models(models.clone()) - - # Return - return - - def spectrum(self): - """ - Return spectrum as FITS object. - """ - # Return - return self.fits - - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Save residual map - self.fits.saveto(self.m_outfile, self.m_clobber) + self._ebounds = self._create_ebounds() # Return return @@ -279,91 +211,86 @@ class csspec(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() - - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") + self._get_parameters() # Write spectral binning into header - if self.logTerse(): - self.log("\n") - self.log.header1("Spectral binning") - if self.m_binned_mode: - cube_ebounds = self.obs[0].events().ebounds() - self.log.parformat("Counts cube energy range") - self.log(str(cube_ebounds.emin())) - self.log(" - ") - self.log(str(cube_ebounds.emax())) - self.log("\n") - for i in range(self.m_ebounds.size()): - self.log.parformat("Bin "+str(i+1)) - self.log(str(self.m_ebounds.emin(i))) - self.log(" - ") - self.log(str(self.m_ebounds.emax(i))) - self.log("\n") + if self._logTerse(): + self._log("\n") + self._log.header1("Spectral binning") + if self._binned_mode: + cube_ebounds = self._obs[0].events().ebounds() + self._log.parformat("Counts cube energy range") + self._log(str(cube_ebounds.emin())) + self._log(" - ") + self._log(str(cube_ebounds.emax())) + self._log("\n") + for i in range(self._ebounds.size()): + self._log.parformat("Bin "+str(i+1)) + self._log(str(self._ebounds.emin(i))) + self._log(" - ") + self._log(str(self._ebounds.emax(i))) + self._log("\n") # 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 model.name() == self._srcname: 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() normpar = model.spectral()[0] - if normpar.is_fixed() and self.logExplicit(): - self.log(" Freeing \""+normpar.name()+"\"\n") + if normpar.is_fixed() and self._logExplicit(): + self._log(" Freeing \""+normpar.name()+"\"\n") normpar.free() - if self.m_calc_ts: + 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 spectrum") - self.log(str(self.m_ebounds)) + if self._logTerse(): + self._log("\n") + self._log.header1("Generate spectrum") + self._log(str(self._ebounds)) # Initialise FITS Table with extension "SPECTRUM" - table = gammalib.GFitsBinTable(self.m_ebounds.size()) + table = gammalib.GFitsBinTable(self._ebounds.size()) table.extname("SPECTRUM") # Add Header for compatibility with gammalib.GMWLSpectrum @@ -371,14 +298,15 @@ class csspec(ctools.cscript): table.card("TELESCOP", "CTA", "Name of Telescope") # Create FITS table columns - energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size()) - energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size()) - energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size()) - flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size()) - flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size()) - TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size()) - ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size()) - Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size()) + nrows = self._ebounds.size() + energy = gammalib.GFitsTableDoubleCol("Energy", nrows) + energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", nrows) + energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", nrows) + flux = gammalib.GFitsTableDoubleCol("Flux", nrows) + flux_err = gammalib.GFitsTableDoubleCol("e_Flux", nrows) + TSvalues = gammalib.GFitsTableDoubleCol("TS", nrows) + ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", nrows) + Npred_values = gammalib.GFitsTableDoubleCol("Npred", nrows) energy.unit("TeV") energy_low.unit("TeV") energy_high.unit("TeV") @@ -387,17 +315,17 @@ class csspec(ctools.cscript): ulim_values.unit("erg/cm2/s") # Loop over energy bins - for i in range(self.m_ebounds.size()): + for i in range(nrows): # Log information - if self.logExplicit(): - self.log("\n") - self.log.header2("Energy bin "+str(i+1)) + if self._logExplicit(): + self._log("\n") + self._log.header2("Energy bin "+str(i+1)) # Get energy boundaries - emin = self.m_ebounds.emin(i) - emax = self.m_ebounds.emax(i) - elogmean = self.m_ebounds.elogmean(i) + emin = self._ebounds.emin(i) + emax = self._ebounds.emax(i) + elogmean = self._ebounds.elogmean(i) elogmean2 = elogmean.MeV() * elogmean.MeV() # Store energy as TeV @@ -407,15 +335,15 @@ class csspec(ctools.cscript): energy_low[i] = (elogmean - emin).TeV() energy_high[i] = (emax - elogmean).TeV() - # use ctselect for unbinned analysis - if not self.m_binned_mode: + # Use ctselect for unbinned analysis + if not self._binned_mode: # 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 = ctools.ctselect(self._obs) select["emin"] = emin.TeV() select["emax"] = emax.TeV() select["tmin"] = "UNDEFINED" @@ -428,15 +356,15 @@ class csspec(ctools.cscript): # Retrieve observation obs = select.obs() - # use ctcubemask for binned analysis + # Use ctcubemask for binned analysis else: # Header - if self.logExplicit(): - self.log.header3("Filtering cube") + if self._logExplicit(): + self._log.header3("Filtering cube") # Select layers - cubemask = ctools.ctcubemask(self.obs) + cubemask = ctools.ctcubemask(self._obs) cubemask["regfile"] = "NONE" cubemask["ra"] = "UNDEFINED" cubemask["dec"] = "UNDEFINED" @@ -449,22 +377,22 @@ class csspec(ctools.cscript): obs = cubemask.obs() # Header - if self.logExplicit(): - self.log.header3("Performing fit") + if self._logExplicit(): + 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.logExplicit(): - self.log("No event in this bin. ") - self.log("Likelihood is zero. ") - self.log("Bin is skipped.") + if self._logExplicit(): + self._log("No event in this bin. ") + self._log("Likelihood is zero. ") + self._log("Bin is skipped.") # Set all values to 0 flux[i] = 0.0 @@ -476,19 +404,19 @@ class csspec(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.logExplicit(): - self.log.header3("Computing upper limit") + if self._logExplicit(): + 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"] = elogmean.TeV() # Try to run upper limit and catch exceptions @@ -496,23 +424,23 @@ class csspec(ctools.cscript): ulimit.run() ulimit_value = ulimit.diff_ulimit() except: - if self.logExplicit(): - self.log("Upper limit calculation failed.") + if self._logExplicit(): + self._log("Upper limit calculation failed.") ulimit_value = -1.0 # Get TS value TS = -1.0 - if self.m_calc_ts: + if self._calc_ts: TS = source.ts() # Compute Npred value (only works for unbinned analysis) Npred = 0.0 - if not self.m_binned_mode: + if not self._binned_mode: for observation in like.obs(): Npred += observation.npred(source) # Get differential flux - fitted_flux = source.spectral().eval(elogmean,gammalib.GTime()) + fitted_flux = source.spectral().eval(elogmean) # Compute flux error parvalue = source.spectral()[0].value() @@ -532,17 +460,17 @@ class csspec(ctools.cscript): ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg # Log information - if self.logTerse(): - self.log("\n") - self.log.parformat("Bin "+str(i+1)) - self.log(str(flux[i])) - self.log(" +/- ") - self.log(str(flux_err[i])) - if self.m_calc_ulimit and ulim_values[i] > 0.0: - self.log(" [< "+str(ulim_values[i])+"]") - self.log(" erg/cm2/s") - if self.m_calc_ts and TSvalues[i] > 0.0: - self.log(" (TS = "+str(TS)+")") + if self._logTerse(): + self._log("\n") + self._log.parformat("Bin "+str(i+1)) + self._log(str(flux[i])) + self._log(" +/- ") + self._log(str(flux_err[i])) + if self._calc_ulimit and ulim_values[i] > 0.0: + self._log(" [< "+str(ulim_values[i])+"]") + self._log(" erg/cm2/s") + if self._calc_ts and TSvalues[i] > 0.0: + self._log(" (TS = "+str(TS)+")") # Append filled columns to fits table table.append(energy) @@ -555,8 +483,69 @@ class csspec(ctools.cscript): table.append(Npred_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 spectrum + self.save() + + # Return + return + + def save(self): + """ + Save spectrum. + """ + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save spectrum") + + # 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 spectrum into \""+self._outfile+"\".\n") + + # Save spectrum + self._fits.saveto(self._outfile, self._clobber) + + # Return + return + + def spectrum(self): + """ + Return spectrum FITS file. + + Returns: + FITS file containing spectrum. + """ + # Return + return self._fits + + def models(self, models): + """ + Set model. + """ + # Copy models + self._obs.models(models.clone()) # Return return @@ -566,14 +555,12 @@ class csspec(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates spectral points. - """ + # Create instance of application app = csspec(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/cstsdist.par b/cscripts/cstsdist.par new file mode 100644 index 0000000..c8922c9 --- /dev/null +++ b/cscripts/cstsdist.par @@ -0,0 +1,57 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cstsdist.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the cstsdist script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inobs, f, h, NONE,,, "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?" +deadc, r, h, 0.95,0,1, "Deadtime correction factor" +outfile, f, a, ts.dat,,, "Output Test Statistic distribution file" + +# +# Script parameters +#================== +ntrials, i, a, 10,,, "Number of trials" +ra, r, a, 83.63,0,360, "RA of pointing (deg)" +dec, r, a, 22.01,-90,90, "Dec of pointing (deg)" +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (TeV)" +enumbins, i, a, 0,,, "Number of energy bins (0=unbinned)" +tmin, r, h, 0.0,,, "Start time (MET in s)" +tmax, r, a, 1800.0,,, "Duration (s)" +npix, i, a, 200,,, "Number of pixels for binned" +binsz, r, a, 0.05,,, "Pixel size for binned (deg/pixel)" +rad, r, h, 5.0,0,180, "Radius of ROI (deg)" +pattern, s, h, single,single|four,, "Observation pattern" +offset, r, h, 1.5,,, "Observation pattern offset (deg)" + +# +# 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, cstsdist.log,,, "Log filename" diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index 2b41d49..0af2a0a 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -1,7 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This script generates the TS distribution for a particular model based -# on Monte-Carlo simulations. +# Generates the TS distribution for a particular model # # Copyright (C) 2011-2016 Juergen Knoedlseder # @@ -31,38 +30,48 @@ import csv # ============== # class cstsdist(ctools.cscript): """ + Generates Test Statistic distribution for a model. + + This class implements the TS distribution generation script. 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 = "cstsdist" - self.version = "1.1.0" + self._name = "cstsdist" + self._version = "1.1.0" # Initialise some members - self.obs = None - self.bkg_model = None - self.full_model = None - - # Make sure that parfile exists - file = self.parfile() + self._obs = gammalib.GObservations() + self._pattern = "single" + self._srcname = "" + self._enumbins = 0 + self._npix = 0 + self._binsz = 0.0 + self._outfile = "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) + 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 @@ -71,219 +80,109 @@ class cstsdist(ctools.cscript): """ 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("inobs","f","h","NONE","","","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","s","a","NONE","","","Input exposure cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("psfcube","s","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("caldb","s","a","prod2","","","Calibration database")) - pars.append(gammalib.GApplicationPar("irf","s","a","South_0.5h","","","Instrument response function")) - pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor")) - pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) - pars.append(gammalib.GApplicationPar("outfile","f","a","ts.dat","","","Output Test Statistics distribution file")) - pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) - pars.append(gammalib.GApplicationPar("ra","r","a","83.63","0","360","RA of pointing (deg)")) - pars.append(gammalib.GApplicationPar("dec","r","a","22.01","-90","90","Dec of pointing (deg)")) - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) - pars.append(gammalib.GApplicationPar("tmin","r","h","0.0","","","Start time (MET in s)")) - pars.append(gammalib.GApplicationPar("tmax","r","a","1800.0","","","Observation duration (in s)")) - pars.append(gammalib.GApplicationPar("npix","i","a","200","","","Number of pixels for binned")) - pars.append(gammalib.GApplicationPar("binsz","r","a","0.05","","","Pixel size for binned (deg/pixel)")) - pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)")) - pars.append(gammalib.GApplicationPar("pattern","s","h","single","","","Observation pattern (single/four)")) - pars.append(gammalib.GApplicationPar("offset","r","h","1.5","","","Observation pattern offset (deg)")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","cstsdist.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + 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: - self.obs = self._get_observations() + if self._obs == None or self._obs.size() == 0: + self._obs = self._get_observations() # Check for requested pattern and use above # observation parameters to set wobble pattern - self.m_pattern = self["pattern"].string() - if self.m_pattern == "four": - self.obs = self.set_obs() + self._pattern = self["pattern"].string() + if self._pattern == "four": + self._obs = self._set_obs() # Get source name - self.m_srcname = self["srcname"].string() + self._srcname = self["srcname"].string() # Get number of energy bins - self.m_enumbins = self["enumbins"].integer() + self._enumbins = self["enumbins"].integer() # Read parameters for binned if requested - if not self.m_enumbins == 0: - self.m_npix = self["npix"].integer() - self.m_binsz = self["binsz"].real() - else: - # Set dummy values (required by obsutils) - self.m_npix = 0 - self.m_binsz = 0.0 + if not self._enumbins == 0: + self._npix = self["npix"].integer() + self._binsz = self["binsz"].real() # 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 other parameters - self.m_outfile = self["outfile"].filename() - self.m_ntrials = self["ntrials"].integer() - self.m_edisp = self["edisp"].boolean() - self.m_offset = self["offset"].real() + self._outfile = self["outfile"].filename() + self._ntrials = self["ntrials"].integer() + self._edisp = self["edisp"].boolean() + self._debug = self["debug"].boolean() - # Set some fixed parameters - self.m_log = False # Logging in client tools - self.m_debug = self["debug"].boolean() # Debugging in client tools + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") # Return return - def models(self, models): + def _set_obs(self): """ - Set model. - """ - # Copy models - self.obs.models(models.copy()) - - # Return - return + Set an observation container. - def execute(self): + Returns: + Observation container. """ - Execute the script. - """ - # Run the script - self.run() - - # Return - return - - 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 input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") - - # Initialise models - full_model, bkg_model = self.set_models() - - # Write models into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Models") - self.log.header2("Background model") - self.log(str(bkg_model)) - self.log("\n\n") - self.log.header2("Full model") - self.log(str(full_model)) - self.log("\n") - - # Write observation into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Observation") - self.log(str(self.obs)) - self.log("\n") - # Write models into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Test model") - self.log(str(self.full_model)) - self.log("\n") - - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate TS distribution") + # Setup observation definition list + obsdeflist = obsutils.set_obs_patterns(self._pattern, + ra=self["ra"].real(), + dec=self["dec"].real(), + offset=self["offset"].real()) - # Loop over trials - for seed in range(self.m_ntrials): + # Create list of observations + obs = obsutils.set_obs_list(obsdeflist, + tstart=self["tmin"].real(), + duration=self["tmax"].real()-self["tmin"].real(), + deadc=self["deadc"].real(), + emin=self["emin"].real(), + emax=self["emax"].real(), + rad=self["rad"].real(), + irf=self["irf"].string(), + caldb=self["caldb"].string()) - # Make a trial - result = self.trial(seed, full_model, bkg_model) + # Return observation container + return obs - # Write out result immediately - if seed == 0: - file = open(self.m_outfile.url(), 'w') - writer = csv.DictWriter(file, result['colnames']) - headers = {} - for n in result['colnames']: - headers[n] = n - writer.writerow(headers) - #writer.writerow(dict((_,_) for _ in result['colnames'])) - else: - file = open(self.m_outfile.url(), 'a') - writer = csv.DictWriter(file, result['colnames']) - writer.writerow(result['values']) - file.close() + def _set_models(self, fitspat=False, fitspec=False): + """ + Set full model and background model. - # Return - return + Kwargs: + fitspat: Fit spatial parameter (default: False). + fitspec: Fit spectral parameters (default: False). - def set_models(self, fitpos=False, fitspec=False): - """ - Set full and background model. + Returns: + Tuple containing full model and background model. """ # Retrieve full model from observation container - full_model = self.obs.models().copy() + full_model = self._obs.models().copy() # Get source model - model = full_model[self.m_srcname] + model = full_model[self._srcname] # Check that model has a Prefactor if not model.has_par("Prefactor"): - msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \ - " Only spectral models with a \"Prefactor\" parameter are supported." + msg = "Model \""+self._srcname+"\" has no parameter "+\ + "\"Prefactor\". Only spectral models with a "+\ + "\"Prefactor\" parameter are supported." raise RuntimeError(msg) # Fit or fix spatial parameters - if fitpos: + if fitspat: if model.has_par("RA"): model["RA"].free() if model.has_par("DEC"): @@ -320,55 +219,35 @@ class cstsdist(ctools.cscript): # Create background model bkg_model = full_model.copy() - bkg_model.remove(self.m_srcname) + bkg_model.remove(self._srcname) # Return models return full_model, bkg_model - def set_obs(self): - """ - Returns an observation container with a set of CTA observations. - - Keywords: - """ - - # Setup observation definition list - obsdeflist = obsutils.set_obs_patterns(self.m_pattern, - ra=self["ra"].real(), - dec=self["dec"].real(), - offset=self["offset"].real()) - - # Create list of observations - obs = obsutils.set_obs_list(obsdeflist, - tstart=self["tmin"].real(), - duration=self["tmax"].real()-self["tmin"].real(), - deadc=self["deadc"].real(), - emin=self["emin"].real(), - emax=self["emax"].real(), - rad=self["rad"].real(), - irf=self["irf"].string(), caldb=self["caldb"].string()) - - # Return observation container - return obs - - def trial(self, seed, full_model, bkg_model): + def _trial(self, seed, full_model, bkg_model): """ Create the TS for a single trial. - Parameters: - seed - Random number generator seed + Args: + seed: Random number generator seed + full_model: Full model + bkg_model: Background model + + Returns: + Result dictionary """ # Write header - if self.logExplicit(): - self.log.header2("Trial "+str(seed+1)) + if self._logExplicit(): + self._log.header2("Trial "+str(seed+1)) # Simulate events - sim = obsutils.sim(self.obs, - nbins=self.m_enumbins, + sim = obsutils.sim(self._obs, + nbins=self._enumbins, seed=seed, - binsz=self.m_binsz, - npix=self.m_npix, - log=self.m_log, debug=self.m_debug) + binsz=self._binsz, + npix=self._npix, + log=self._log_clients, + debug=self._debug) # Determine number of events in simulation nevents = 0.0 @@ -376,72 +255,76 @@ class cstsdist(ctools.cscript): nevents += run.events().number() # Write simulation results - if self.logExplicit(): - self.log.header3("Simulation") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") + if self._logExplicit(): + self._log.header3("Simulation") + self._log.parformat("Number of simulated events") + self._log(nevents) + self._log("\n") # Fit background only sim.models(bkg_model) - like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + like_bgm = obsutils.fit(sim, + log=self._log_clients, + debug=self._debug) result_bgm = like_bgm.obs().models().copy() LogL_bgm = like_bgm.opt().value() npred_bgm = like_bgm.obs().npred() # Write background fit results - if self.logExplicit(): - self.log.header3("Background model fit") - self.log.parformat("log likelihood") - self.log(LogL_bgm) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_bgm) - self.log("\n") + if self._logExplicit(): + self._log.header3("Background model fit") + self._log.parformat("log likelihood") + self._log(LogL_bgm) + self._log("\n") + self._log.parformat("Number of predicted events") + self._log(npred_bgm) + self._log("\n") for model in result_bgm: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") + self._log(str(par)+"\n") # Fit background and test source sim.models(full_model) - like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + like_all = obsutils.fit(sim, + log=self._log_clients, + debug=self._debug) result_all = like_all.obs().models().copy() LogL_all = like_all.opt().value() npred_all = like_all.obs().npred() ts = 2.0*(LogL_bgm-LogL_all) # Write background and test source fit results - if self.logExplicit(): - self.log.header3("Background and test source model fit") - self.log.parformat("Test statistics") - self.log(ts) - self.log("\n") - self.log.parformat("log likelihood") - self.log(LogL_all) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_all) - self.log("\n") + if self._logExplicit(): + self._log.header3("Background and test source model fit") + self._log.parformat("Test statistics") + self._log(ts) + self._log("\n") + self._log.parformat("log likelihood") + self._log(LogL_all) + self._log("\n") + self._log.parformat("Number of predicted events") + self._log(npred_all) + self._log("\n") for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") + self._log.parformat("Model") + self._log(model.name()) + self._log("\n") for par in model: - self.log(str(par)+"\n") + self._log(str(par)+"\n") # Write result - elif self.logTerse(): - self.log.parformat("Trial "+str(seed)) - self.log("TS=") - self.log(ts) - self.log(" Prefactor=") - self.log(result_all[self.m_srcname]["Prefactor"].value()) - self.log("+/-") - self.log(result_all[self.m_srcname]["Prefactor"].error()) - self.log("\n") + elif self._logTerse(): + self._log.parformat("Trial "+str(seed)) + self._log("TS=") + self._log(ts) + self._log(" Prefactor=") + self._log(result_all[self._srcname]["Prefactor"].value()) + self._log("+/-") + self._log(result_all[self._srcname]["Prefactor"].error()) + self._log("\n") # Initialise results colnames = [] @@ -498,18 +381,105 @@ class cstsdist(ctools.cscript): return result + # Public methods + def run(self): + """ + Run the script. + """ + # Switch screen logging on in debug mode + if self._logDebug(): + self._log.cout(True) + + # Get parameters + self._get_parameters() + + # Initialise models + full_model, bkg_model = self._set_models() + + # Write models into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Models") + self._log.header2("Background model") + self._log(str(bkg_model)) + self._log("\n\n") + self._log.header2("Full model") + self._log(str(full_model)) + self._log("\n") + + # Write observation into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Observation") + self._log(str(self._obs)) + self._log("\n") + + # Write models into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Test model") + self._log(str(full_model)) + self._log("\n") + + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Generate TS distribution") + + # Loop over trials + for seed in range(self._ntrials): + + # Make a trial + result = self._trial(seed, full_model, bkg_model) + + # Write out result immediately + if seed == 0: + file = open(self._outfile.url(), 'w') + writer = csv.DictWriter(file, result['colnames']) + headers = {} + for n in result['colnames']: + headers[n] = n + writer.writerow(headers) + else: + file = open(self._outfile.url(), 'a') + writer = csv.DictWriter(file, result['colnames']) + writer.writerow(result['values']) + file.close() + + # Return + return + + def models(self, models): + """ + Set model. + """ + # Copy models + self._obs.models(models.copy()) + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return + + # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates TS distribution. - """ + # Create instance of application app = cstsdist(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/cscripts/cstsmapmerge.par b/cscripts/cstsmapmerge.par new file mode 100644 index 0000000..f195ec7 --- /dev/null +++ b/cscripts/cstsmapmerge.par @@ -0,0 +1,38 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cstsmapmerge.par # +# Version: 1.1.0 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the cstsmapmerge script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inmaps, s, a, tsmap.fits,,, "Input TS map FITS files" +outmap, f, a, NONE,,, "Output TS map FITS file" + +# +# Script parameters +#================== +overwrite, b, h, yes,,, "Overwrite previously filled bins?" +delete, b, h, no,,, "Delete merged files?" + +# +# 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, cstsmapmerge.log,,, "Log filename" + + diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py index 2ef812f..af41ea8 100755 --- a/cscripts/cstsmapmerge.py +++ b/cscripts/cstsmapmerge.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # ========================================================================== -# Merges sliced TS map files into one map +# Merge Test Statistic maps # # Copyright (C) 2015-2016 Michael Mayer # @@ -30,30 +30,42 @@ import os # ================== # class cstsmapmerge(ctools.cscript): """ - This class merges TS map files into one TS map + Merge Test Statistic maps. """ + + # Constructors and destructors def __init__(self, *argv): """ Constructor. """ # Set name - self.name = "cstsmapmerge" - self.version = "1.1.0" + self._name = "cstsmapmerge" + self._version = "1.1.0" - # Make sure that parfile exists - file = self.parfile() + # Initialise class members + self._files = None + self._in_filename = "" + self._tsmap = gammalib.GSkyMap() # Empty sky map + self._statusmap = gammalib.GSkyMap() # Empty sky map + self._maps = [] + self._mapnames = [] + 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) + 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,98 +77,79 @@ class cstsmapmerge(ctools.cscript): # 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 that parfile was not found - print("Parfile \""+parfile+"\" not found. Create default parfile.") - - # Create default parfile - pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("inmaps","s","a","tsmap.fits","","","Input TS map FITS files")) - pars.append(gammalib.GApplicationPar("outmap","f","a","NONE","","","Output TS map FITS file")) - pars.append(gammalib.GApplicationPar("overwrite","b","h","yes","","","Overwrite previously filled bins?")) - pars.append(gammalib.GApplicationPar("delete","b","h","no","","","Delete merged files?")) - pars.append_standard() - pars.append(gammalib.GApplicationPar("logfile","f","h","cstsmapmerge.log","","","Log filename")) - pars.save(parfile) - - # Return - return - def get_parameters(self): + # Private methods + def _get_parameters(self): """ - Get parameters from parfile and setup the observation. + Get parameters from parfile. """ - # Get input models - self.inmodels = self["inmaps"].string() + # Get input maps string + inmaps = self["inmaps"].string() - # Interpet inmodels string # Handle ascii files - if "@" == self.inmodels[0]: - self.files = open(self.inmodels.replace("@","")).read().splitlines() + if "@" == inmaps[0]: + self._files = open(inmaps.replace("@","")).read().splitlines() # Handle wild card strings - elif "*" in self.inmodels: - self.files = glob.glob(self.inmodels) + elif "*" in inmaps: + self._files = glob.glob(inmaps) # Handle space-separated list - elif " " in self.inmodels: - self.files = self.inmodels.split(" ") + elif " " in inmaps: + self._files = inmaps.split(" ") - # Handle comma-separated list - elif "," in self.inmodels: - self.files = self.inmodels.split(",") + # Handle semi-colon separated list + elif ";" in inmaps: + self._files = inmaps.split(";") # Throw exception if input models cannot be decoded else: - msg = "Parameter \"inmmaps\" must contain either an @ASCII file,"+\ - " a comma-separated or whitespace-separated list of files"+\ - " or a wildcard string" + msg = "Parameter \"inmaps\" must contain either an @ASCII "\ + "file, a semi-colon-separated or whitespace-separated "\ + "list of files or a wildcard string." raise RuntimeError(msg) # Check number of files. We need at least two files. - if len(self.files) <= 1: + if len(self._files) <= 1: msg = "Need at least two files to start merging, "+\ - str(len(self.files))+" files given." + str(len(self._files))+" file(s) given." raise RuntimeError(msg) - - # Get output file - self.outmap = self["outmap"].filename() - + # Get other parameters - self.overwrite = self["overwrite"].boolean() - self.delete = self["delete"].boolean() + self._overwrite = self["overwrite"].boolean() + self._delete = self["delete"].boolean() + # Read ahead output filename + if self._read_ahead: + self._outmap = self["outmap"].filename() + + # Write input parameters into logger + if self._logTerse(): + self._log_parameters() + self._log("\n") + # Return return - def init_map(self, fitsfile): + def _init_ts_map(self, fitsfile): """ + Initialise Test Statistic map. """ # Set filename - self.in_filename = fitsfile + self._in_filename = fitsfile # Open FITS file - fits = gammalib.GFits(fitsfile) - self.tsmap = gammalib.GSkyMap() - self.tsmap.read(fits[0]) - self.statusmap = gammalib.GSkyMap() - self.statusmap.read(fits["STATUS MAP"]) + fits = gammalib.GFits(fitsfile) + + # Read TS and status maps + self._tsmap = gammalib.GSkyMap() + self._tsmap.read(fits[0]) + self._statusmap = gammalib.GSkyMap() + self._statusmap.read(fits["STATUS MAP"]) # Get other maps - self.maps = [] - self.mapnames = [] + self._maps = [] + self._mapnames = [] # Loop over extensions for hdu in fits: @@ -167,8 +160,8 @@ class cstsmapmerge(ctools.cscript): # Add present maps skymap = gammalib.GSkyMap() skymap.read(hdu) - self.maps.append(skymap) - self.mapnames.append(hdu.extname()) + self._maps.append(skymap) + self._mapnames.append(hdu.extname()) # Close FITS file fits.close() @@ -176,11 +169,17 @@ class cstsmapmerge(ctools.cscript): # Return return - # Add map instances to each other. In this way the TS maps get merged - def add(self, fitsfile): - - # open FITS file - fits = gammalib.GFits(fitsfile) + def _merge_ts_map(self, fitsfile): + """ + Merge TS map from FITS file into output TS map. + + Args: + fitsfile: FITS file to be merged. + """ + # Open FITS file + fits = gammalib.GFits(fitsfile) + + # Read TS and status maps add_tsmap = gammalib.GSkyMap() add_tsmap.read(fits[0]) add_statusmap = gammalib.GSkyMap() @@ -203,175 +202,240 @@ class cstsmapmerge(ctools.cscript): # Close FITS file fits.close() - # Compare size of other maps - if not len(add_maps) == len(self.maps): - self.error("Maps cannot be added: Maps from \""+self.in_filename+"\" do not match from \""+fitsfile+"\"") + # Compare size of maps + if not len(add_maps) == len(self._maps): + msg = "Cannot merge map \""+fitsfile+"\" into map \""+\ + self._in_filename+"\" since the number of parameters "\ + "between both maps is different." + raise RuntimeError(msg) # Loop over bins - for i in range(self.tsmap.npix()): + for i in range(self._tsmap.npix()): + + # Consider only bins that have been computed if add_statusmap[i] > 0.5: - # Throw exception if this bin has already been computed - if self.statusmap[i] > 0.5 and not self.overwrite: - msg = "Attempt to merge bin which apparently has already"+\ - " been merged. File \""+fitsfile+"\" contains"+\ - " already merged bins. Set hidden parameter"+\ - " overwrite=yes to avoid this error." + # Raise exception if this bin has already been computed + if self._statusmap[i] > 0.5 and not self._overwrite: + msg = "Attempt to merge bin which apparently has "+\ + "already been merged. File \""+fitsfile+"\" "+\ + "contains already merged bins. Set hidden "+\ + "parameter \"overwrite=yes\" to avoid this error." raise RuntimeError(msg) # Copy TS values - self.tsmap[i] = add_tsmap[i] + self._tsmap[i] = add_tsmap[i] # Copy status - self.statusmap[i] = add_statusmap[i] + self._statusmap[i] = add_statusmap[i] # Loop over maps and copy entries - for j in range(len(self.maps)): - self.maps[j][i] = add_maps[j][i] + for j in range(len(self._maps)): + self._maps[j][i] = add_maps[j][i] # Return return + def _get_number_of_ts_pixels(self, fitsfile): + """ + Return number of pixels with TS values. + + Args: + fitsfile: FITS file to be merged. + + Returns: + Number of pixels for which TS has been computed. + """ + # Get status map for this file + status = gammalib.GSkyMap(fitsfile+"[STATUS MAP]") + + # Count number of pixels with TS status set + count = 0 + for pix in status: + if pix > 0.5: + count += 1 + + # Return number of pixels + return count + + + # 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() + self._get_parameters() - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") - - if self.logTerse(): - self.log("\n") - self.log.header1("Merging sliced TS map files") - self.log("\n") + # Write header into logger + if self._logTerse(): + self._log("\n") + self._log.header1("Merge TS maps") # Initialise file to start with - file0 = "" - - self.added_files = [] + file0 = "" + self._merged_files = [] - # Test files for the entry status map - # use the first one to appear useful - for fitsfile in self.files: + # Test files for the entry status map. Use the first one to appear + # useful + for fitsfile in self._files: + + # Skip file if it's not a FITS file if not gammalib.GFilename(fitsfile).is_fits(): - self.log("Skipping file \""+fitsfile+" (not a FITS file)\n") + if self._logExplicit(): + self._log(gammalib.parformat("Skip file")) + self._log(fitsfile) + self._log(" (not a FITS file)\n") + continue + + # Open FITS file fits = gammalib.GFits(fitsfile) - + + # If file contains a status map then use it if fits.contains("STATUS MAP"): - self.log("Using file \""+fitsfile+"\" as initialising map\n") + fits.close() file0 = fitsfile + if self._logTerse(): + self._log(gammalib.parformat("Initial TS map file")) + self._log(fitsfile) + count = self._get_number_of_ts_pixels(fitsfile) + self._log(" (%d TS pixels computed)" % count) + self._log("\n") break + + # ... otherwise signal that file is useless else: - if self.logExplicit(): - self.log("File \""+fitsfile+"\" has not \"STATUS MAP\" extension\n") - - # Print an error message if no suitable file was provided + fits.close() + if self._logExplicit(): + self._log(gammalib.parformat("Skip file")) + self._log(fitsfile) + self._log(" (no \"STATUS MAP\" extension)\n") + continue + + # Signal if no suitable file was found if file0 == "": - if self.logTerse(): - self.log("None of the provided file seems to be a sliced ts map file (no \"STATUS MAP\" extension)\n") + if self._logTerse(): + self._log("None of the provided files seems to be a sliced " + "TS map file (none has a \"STATUS MAP\" " + "extension).\n") - # Otherwise merge files + # ... otherwise merge files else: - # copy file list - workfiles = self.files + # Copy file list + workfiles = self._files # Remove entry which will be used to initalise the map workfiles.remove(file0) # Initialise map from first file - self.init_map(file0) + self._init_ts_map(file0) # Append to added files - self.added_files.append(file0) + self._merged_files.append(file0) # Loop over files for fitsfile in workfiles: # Skip if file is not FITS if not gammalib.GFilename(fitsfile).is_fits(): - if self.logTerse(): - self.log("Skipping file \""+fitsfile+" (not a FITS file)\n") + if self._logExplicit(): + self._log(gammalib.parformat("Skip file")) + self._log(fitsfile) + self._log(" (not a FITS file)\n") continue + + # Open FITS file + fits = gammalib.GFits(fitsfile) # Skip if file does not contain status map - fits = gammalib.GFits(fitsfile) if not fits.contains("STATUS MAP"): - if self.logTerse(): - self.log("Skipping file \""+fitsfile+" (not a FITS file)\n") + fits.close() + if self._logExplicit(): + self._log(gammalib.parformat("Skip file")) + self._log(fitsfile) + self._log(" (no \"STATUS MAP\" extension)\n") continue + + # Close FITS file + fits.close() # Logging - if self.logTerse(): - self.log("Adding TS map slice \""+fitsfile+"\"\n") + if self._logTerse(): + self._log(gammalib.parformat("Merge TS map file")) + self._log(fitsfile) + count = self._get_number_of_ts_pixels(fitsfile) + self._log(" (%d TS pixels computed)" % count) + self._log("\n") - # Add file - self.add(fitsfile) + # Merge TS map + self._merge_ts_map(fitsfile) - # Append to added files - self.added_files.append(fitsfile) + # Append FITS file to merged files + self._merged_files.append(fitsfile) # Return return def save(self): """ - Save merged tsmap file and remove slices if requested + Save TS map and remove slices if requested. """ - - # Logging - if self.logExplicit(): - self.log("\n") - self.log.header1("Saving merged TS map") - self.log("\n") + # Write header + if self._logTerse(): + self._log("\n") + self._log.header1("Save TS map") + + # Get output filename in case it was not read ahead + self._outmap = self["outmap"].filename() + + # Log filename + if self._logTerse(): + self._log(gammalib.parformat("TS map file")) + self._log(self._outmap.url()) + self._log("\n") # Create FITS file fits = gammalib.GFits() # Write TS map into primary - self.tsmap.write(fits) + self._tsmap.write(fits) # Loop over maps and write them to fits - for i in range(len(self.maps)): - self.maps[i].write(fits) + for i in range(len(self._maps)): + self._maps[i].write(fits) # Set map names as extensions - for i in range(len(self.mapnames)): - fits[i+1].extname(self.mapnames[i]) + for i in range(len(self._mapnames)): + fits[i+1].extname(self._mapnames[i]) # Check if map is fully done done = True - for pix in self.statusmap: + for pix in self._statusmap: if pix < 0.5: done = False break - # Write Status map if we are not done yet + # Write status map if we are not done yet if not done: - self.statusmap.write(fits) + self._statusmap.write(fits) fits[fits.size()-1].extname("STATUS MAP") # Save FITS file - fits.saveto(self.outmap, self.clobber()) - if self.logExplicit(): - self.log("Saved output map to \""+self.outmap+"\"\n") - if self.delete: - if self.logTerse(): - self.log("\n") - self.log.header1("Deleting slices") - self.log("\n") - for filename in self.added_files: - if self.logExplicit(): - self.log("Deleting file \""+filename+"\"\n") - os.remove(filename) + fits.saveto(self._outmap, self._clobber()) + + # Delete TS input maps if requested + if self._delete: + for filename in self._merged_files: + os.remove(filename) + if self._logTerse(): + self._log(gammalib.parformat("Deleted input file")) + self._log(filename) + self._log("\n") # Return return @@ -380,10 +444,13 @@ class cstsmapmerge(ctools.cscript): """ Execute the script. """ + # Set read ahead flag + self._read_ahead = True + # Run the script self.run() - # Save ds9 file if required + # Save TS map if required self.save() # Return @@ -394,14 +461,12 @@ class cstsmapmerge(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Merges ts map slices. - """ + # Create instance of application app = cstsmapmerge(sys.argv) # Open logfile - app.logFileOpen() + app._logFileOpen() # Execute application app.execute() diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index cbd8fe1..0f66895 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -173,9 +173,9 @@ else $ECHO " merged_models.xml file is not valid" exit 1 fi -# === comma-separated === +# === semi-colon separated === @PYTHON_COVERAGE@@abs_top_srcdir@/cscripts/csmodelmerge.py \ - inmodels="data/crab.xml,data/source1.xml" \ + inmodels="data/crab.xml;data/source1.xml" \ outmodel="merged_models2.xml" > csmodelmerge.log $ECHO -n "." if [ -s "merged_models2.xml" ] @@ -218,10 +218,10 @@ rm models.txt # Test cstsmapmerge # ================= $ECHO -n "Test cstsmapmerge: " -# === comma-separated === +# === semi-colon separated === chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py @PYTHON_COVERAGE@@abs_top_srcdir@/cscripts/cstsmapmerge.py \ - inmaps="data/tsmap_0.fits,data/tsmap_1.fits" \ + inmaps="data/tsmap_0.fits;data/tsmap_1.fits" \ outmap="merged_tsmap1.fits" > cstsmapmerge.log $ECHO -n "." if [ -s "merged_tsmap1.fits" ] -- libgit2 0.21.4