diff --git a/.gitignore b/.gitignore index be454c1..51710b2 100644 --- a/.gitignore +++ b/.gitignore @@ -70,6 +70,8 @@ test/*.trs test/test_cscripts.sh test/reports/*.xml test/butterfly.txt +test/model.reg +test/obs.reg examples/*.fits examples/*.xml examples/*.dat diff --git a/ChangeLog b/ChangeLog index 25c6cb0..1a32fc4 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2015-11-24 Juergen Knoedlseder + + * Version 1.1.0 released + ======================== + + Add csobsinfo, csmodelinfo, csmodelmerge, cstsmapmerge (#1578) + + 2015-11-19 Juergen Knoedlseder * Version 1.0.0 released diff --git a/NEWS b/NEWS index 679749a..eefb17a 100644 --- a/NEWS +++ b/NEWS @@ -1,15 +1,15 @@ -New Features and Important Changes in ctools 1.0.0 +New Features and Important Changes in ctools 1.1.0 Author(s) : Juergen Knoedlseder -19 November 2015 +24 November 2015 1. Introduction --------------- This document briefly summarizes executable by executable the most important changes made since the last release of ctools. The code complies to GammaLib -version 1.0.0 or above. +version 1.1.0 or above. 2. Configuration @@ -82,8 +82,7 @@ None 15. ctcubemask - CTA cube filtering ----------------------------------- -Allow to ignore the spatial filtering by specifying NONE parameters for "ra", -"dec" or "rad". +None 16. ctbutterfly - Generate butterfly diagrams @@ -103,10 +102,11 @@ None 19. cscripts ------------ -The csspec script now handled correctly stacked data. Dependent on the -requested energy binning, it selects for each spectral energy bin the -relevant energy bins from the counts cube and performs a maximum likelihood -analysis. +The following scripts have been added: +- csobsinfo: dump information in observation container +- csmodelinfo: dump information in model container +- csmodelmerge: merge model containers into single container +- cstsmapmerge: merge sliced TS map into single map 20. Examples diff --git a/README b/README index acd9ba5..f7d9f7f 100644 --- a/README +++ b/README @@ -1,8 +1,8 @@ ctools information ================== -Version: 1.0.0 (21 November 2015) +Version: 1.1.0 (24 November 2015) Author: Juergen Knodlseder (jurgen.knodlseder@irap.omp.eu) -GammaLib dependency: 1.0.0 +GammaLib dependency: 1.1.0 License information diff --git a/configure.ac b/configure.ac index 5424b7a..c75b2d9 100644 --- a/configure.ac +++ b/configure.ac @@ -29,7 +29,7 @@ ############################################################################# # Initialisation # ############################################################################# -AC_INIT([ctools], [1.0.0], [jurgen.knodlseder@irap.omp.eu], [ctools]) +AC_INIT([ctools], [1.1.0.dev], [jurgen.knodlseder@irap.omp.eu], [ctools]) AC_CONFIG_SRCDIR([src/ctselect/ctselect.cpp]) AC_CONFIG_MACRO_DIR([m4]) #AM_INIT_AUTOMAKE([gnu no-dependencies]) diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 8140284..f2a8d6a 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -1,19 +1,16 @@ ## Process this file with automake to produce Makefile.in -# Scripts to be distributed in bindir -binscripts = tsmerge.py - # Scripts to be distributed in $(pyexecdir)/cscripts pythonscripts = __init__.py \ cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \ - cshessobs.py csobsdef.py csspec.py cslightcrv.py \ - csinfo.py \ + cshessobs.py csobsdef.py csspec.py cslightcrv.py csobsinfo.py \ + csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \ obsutils.py # Scripts for which symbolic links will be created in bindir linkscripts = cspull cssens cstsdist cscaldb csresmap \ - cshessobs csobsdef csspec cslightcrv \ - csinfo + cshessobs csobsdef csspec cslightcrv csobsinfo \ + csinfo csmodelinfo csmodelmerge cstsmapmerge # Distribute scripts EXTRA_DIST = $(binscripts) $(pythonscripts) diff --git a/cscripts/csmodelinfo.py b/cscripts/csmodelinfo.py new file mode 100755 index 0000000..2356984 --- /dev/null +++ b/cscripts/csmodelinfo.py @@ -0,0 +1,511 @@ +#!/usr/bin/env python +# ========================================================================== +# Dump information about models into log file +# +# Copyright (C) 2015 Michael Mayer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools +import sys + + +# ================= # +# csmodelinfo class # +# ================= # +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. + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "csmodelinfo" + self.version = "1.1.0" + + # Make sure that parfile exists + file = self.parfile() + + # Initialise application + if len(argv) == 0: + ctools.cscript.__init__(self, self.name, self.version) + elif len(argv) == 1: + ctools.cscript.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + + # Set logger properties + self.log_header() + self.log.date(True) + + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Return + return + + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a + parfile. + """ + + # Set parfile name + parfile = self.name+".par" + + try: + pars = gammalib.GApplicationPars(parfile) + except: + # Signal 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): + """ + 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() + + # 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() + + # 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") + + if self.logVerbose(): + self.log("\n") + self.log.header1("Models in detail") + self.log("\n") + self.log(str(self.models)) + self.log("\n") + + # Initialise output values + types = {} + instruments = {} + ts = {} + pars_at_limit = {} + free_src_pars = {} + n_par_total = 0 + n_par_free = 0 + n_par_fixed = 0 + n_par_at_limit = 0 + n_par_free_bkg = 0 + n_par_free_src = 0 + n_par_free_spec = 0 + n_par_free_spat = 0 + n_par_free_temp = 0 + + # Loop over models and initialise dictionaries with model names + 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: + + # Model types + types[model.type()] += 1 + + # Model instruments + instruments[model.instruments()] += 1 + + # TS + if model.has_ts(): + ts[model.name()] = model.ts() + + # Skymodel? + skymodel = (model.classname() == "GModelSky") + if skymodel: + free_src_pars[model.name()] = [] + + # Loop over model parameters + for par in model: + + # Increment total parameters + n_par_total += 1 + + # Check if parameter is at limit + if par.has_range(): + if par.value() == par.min() or par.value() == par.max(): + pars_at_limit[model.name()].append(par.name()) + n_par_at_limit += 1 + + # Check if model parameter is free/fixed and + # increment respective variable + if par.is_free(): + n_par_free += 1 + if skymodel: + free_src_pars[model.name()].append(par) + n_par_free_src += 1 + else: + n_par_free_bkg += 1 + else: + n_par_fixed += 1 + + # If sky model, divide free parameters into spatial, spectral and temporal + if skymodel: + for par in model.spectral(): + if par.is_free(): + n_par_free_spec += 1 + for par in model.spatial(): + if par.is_free(): + n_par_free_spat += 1 + for par in model.temporal(): + if par.is_free(): + n_par_free_temp += 1 + + # Log summary + self.log("\n") + self.log.header1("Summary") + self.log("\n") + + # Log instruments + self.log.header3("Instrument specific models") + for inst, n_inst in instruments.iteritems(): + if inst == "": + inst = "All" + self.log.parformat(inst) + self.log(str(n_inst)) + self.log("\n") + + # Log model types + self.log("\n") + self.log.header3("Model types") + for modeltype, n_types in types.iteritems(): + self.log.parformat(modeltype) + self.log(str(n_types)) + self.log("\n") + + # Log parameter information + self.log("\n") + self.log.header2("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") + 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.iteritems(): + 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.iteritems(): + 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.iteritems(): + self.log.parformat(source) + self.log("%.2f"%tsvalue) + self.log("\n") + + # Return + return + + def save(self): + """ + Save models to ds9 region file if required + """ + + # Check if ds9 file is valid + if not self.ds9file == "NONE": + + self.log("\n") + self.log.header1("Save region file") + self.log("\n") + + # Open file + f = open(self.ds9file,"w") + + # Write coordinate system + f.write("fk5\n") + + # Loop over 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 len(line): + f.write(line+"\n") + + # Logging for diffuse components + elif model.type() == "DiffuseSource": + self.log("Skipping diffuse model "+model.name()+"\n") + + # Logging for background components + else: + if self.logExplicit(): + self.log("Skipping background model "+model.name()+"\n") + + # Close file + f.close() + + # 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. + """ + # Run the script + self.run() + + # Save ds9 file if required + self.save() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + """ + Generates model information output. + """ + # Create instance of application + app = csmodelinfo(sys.argv) + + # Open logfile + app.logFileOpen() + + # Execute application + app.execute() diff --git a/cscripts/csmodelmerge.py b/cscripts/csmodelmerge.py new file mode 100755 index 0000000..ff67612 --- /dev/null +++ b/cscripts/csmodelmerge.py @@ -0,0 +1,218 @@ +#!/usr/bin/env python +# ========================================================================== +# Merges model XML files into one file +# +# Copyright (C) 2015 Michael Mayer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools +import sys +import glob + + +# ================== # +# csmodelmerge class # +# ================== # +class csmodelmerge(ctools.cscript): + """ + This class merges an arbitray number of model container into one output + model XML file. + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "csmodelmerge" + self.version = "1.1.0" + + # Make sure that parfile exists + file = self.parfile() + + # Initialise application + if len(argv) == 0: + ctools.cscript.__init__(self, self.name, self.version) + elif len(argv) == 1: + ctools.cscript.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + + # Set logger properties + self.log_header() + self.log.date(True) + + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Return + return + + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a + parfile. + """ + + # Set parfile name + parfile = self.name+".par" + + try: + pars = gammalib.GApplicationPars(parfile) + except: + # Signal 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): + """ + Get parameters from parfile and setup the observation. + """ + # Get input models + self.inmodels = self["inmodels"].string() + + # Interpet inmodels string + # Handle ascii files + if "@" == self.inmodels[0]: + self.files = open(self.inmodels.replace("@","")).read().splitlines() + + # Handle wild card strings + elif "*" in self.inmodels: + self.files = glob.glob(self.inmodels) + + # Handle space-separated list + elif " " in self.inmodels: + self.files = self.inmodels.split(" ") + + # Handle comma-separated list + elif "," in self.inmodels: + self.files = self.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" + raise gammalib.GException.invalid_argument("csmodelmerge::get_parameters", msg) + + # Get output file + self.outmodel = self["outmodel"].filename() + + # 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") + + # Dump header + if self.logTerse(): + self.log("\n") + self.log.header1("Merging Models") + self.log("\n") + + # Initialise output model container + self.models = gammalib.GModels() + + # Loop over model files + for modelfile in self.files: + + # Logging + self.log("Adding "+modelfile) + self.log("\n") + + # Open models + model_add = gammalib.GModels(modelfile) + + # Append models to container + for model in model_add: + self.models.append(model) + + # Return + return + + def save(self): + """ + Save pointings to ds9 region file if required + """ + # Logging + if self.logExplicit(): + self.log("\n") + self.log.header1("Saving model file") + self.log("\n") + + # Save models + self.models.save(self.outmodel) + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Save ds9 file if required + self.save() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + """ + Merges model XML files + """ + # Create instance of application + app = csmodelmerge(sys.argv) + + # Open logfile + app.logFileOpen() + + # Execute application + app.execute() diff --git a/cscripts/csobsinfo.py b/cscripts/csobsinfo.py new file mode 100755 index 0000000..c438563 --- /dev/null +++ b/cscripts/csobsinfo.py @@ -0,0 +1,465 @@ +#!/usr/bin/env python +# ========================================================================== +# Dump information about observation into log file +# +# Copyright (C) 2015 Michael Mayer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools +import sys + + +# =============== # +# csobsinfo class # +# =============== # +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. + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "csobsinfo" + self.version = "1.1.0" + + # Make sure that parfile exists + file = self.parfile() + + # Initialise application + if len(argv) == 0: + ctools.cscript.__init__(self, self.name, self.version) + elif len(argv) == 1: + ctools.cscript.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + + # Set logger properties + self.log_header() + self.log.date(True) + + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Return + return + + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a + parfile. + """ + + # Set parfile name + parfile = self.name+".par" + + try: + pars = gammalib.GApplicationPars(parfile) + except: + # Signal 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): + """ + Get parameters from parfile and setup the observation. + """ + # Get parameters + self.require_inobs("csobsinfo::get_parameters") + self.obs = self.get_observations(False) + + # Initialise object position + self.obj_dir = gammalib.GSkyDir() + + # Get (optional) offset parameters + self.m_offset = self["offset"].boolean() + if self.m_offset: + ra = self["ra"].real() + dec = self["dec"].real() + self.obj_dir.radec_deg(ra,dec) + + # Get (optional) DS9 filename + self.ds9file = self["ds9file"].filename() + + # 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 arrays to store certain values for reuse + # Todo, think about using a python dictionary + self.offsets = [] + self.zeniths = [] + self.azimuths = [] + self.pnt_ra = [] + self.pnt_dec = [] + obs_names = [] + + # Initialise output to be filled + self.ebounds = gammalib.GEbounds() + self.gti = gammalib.GGti() + ontime = 0.0 + livetime = 0.0 + n_events = 0 + n_eventbins = 0 + n_obs_binned = 0 + n_obs_unbinned = 0 + + # Logging + self.log.header1("Looping over observations") + + # Loop over observations + for obs in self.obs: + + + # Skip non-CTA observations + if not obs.classname() == "GCTAObservation": + self.log("Skipping "+obs.instrument()+" observation\n") + continue + + # Use observed object as observation name + # if name is not given + obs_name = obs.name() + if obs_name == "": + obs_name = obs.object() + + # Logging + if self.logTerse(): + self.log.header2(obs_name) + + # Retrieve observation name + obs_names.append(obs_name) + + # Retrieve energy boundaries + obs_bounds = obs.events().ebounds() + + # Retrieve time interval + obs_gti = obs.events().gti() + + # Compute mean time and dead time fraction in percent + tmean = (obs_gti.tstart() + obs_gti.tstop()) + tmean *= 0.5 + deadfrac = (1.0-obs.deadc(tmean))*100.0 + + # Retrieve pointing and store Ra,Dec + pnt_dir = obs.pointing().dir() + self.pnt_ra.append(pnt_dir.ra_deg()) + self.pnt_dec.append(pnt_dir.dec_deg()) + + # If avaliable append energy boundaries + if obs_bounds.size() > 0 : + self.ebounds.merge(obs_bounds.emin(),obs_bounds.emax()) + + # Append time interval + self.gti.merge(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") + + 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") + 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") + + # Log option + if self.logTerse(): + + # Log the observation energy range (if available) + self.log.parformat("Energy range") + if obs_bounds.size() == 0: + self.log("undefined") + else: + self.log(str(obs_bounds.emin())+" - "+str(obs_bounds.emax())) + self.log("\n") + + # Log observation time interval + self.log.parformat("Time range (MJD)") + if obs_gti.size() == 0: + self.log("undefined") + else: + self.log(str(obs_gti.tstart().mjd())+" - "+str(obs_gti.tstop().mjd())) + self.log("\n") + + # Log ontime + self.log.parformat("Ontime [s]") + self.log(str(obs.ontime())) + self.log("\n") + + # Log livetime + self.log.parformat("Livetime [s]") + self.log(str(obs.livetime())) + self.log("\n") + + # Log dead time fraction + self.log.parformat("Dead-time fraction (%)") + self.log("%.3f" % (deadfrac)) + self.log("\n") + + # Log pointing direction + 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.offsets.append(offset) + self.log.parformat("Offset") + self.log("%.2f" % (offset)) + self.log("\n") + else: + self.offsets.append(-1.0) + + # 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) + + # 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") + + # Log summary + 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 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.offsets)/len(self.offsets))) + 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("\n") + self.log.parformat("Mean azimuth angle") + self.log("%.2f" % (sum(self.azimuths)/len(self.azimuths))) + self.log("\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") + + # Log energy range + self.log.header3("Energy range") + self.log.parformat("Emin") + if self.ebounds.size() == 0: + self.log("undefined") + else: + self.log(str(self.ebounds.emin())) + self.log("\n") + self.log.parformat("Emax") + if self.ebounds.size() == 0: + self.log("undefined") + else: + self.log(str(self.ebounds.emax())) + self.log("\n\n") + + # Log time range + self.log.header3("Time range") + self.log.parformat("Start [MJD]") + self.log(str(self.gti.tstart().mjd())) + self.log("\n") + self.log.parformat("Stop [MJD]") + self.log(str(self.gti.tstop().mjd())) + self.log("\n") + + # Log ontime and livetime in different units + self.log.parformat("Total ontime") + self.log("%.2f s = %.2f min = %.2f h" % (ontime, ontime/60., ontime/3600.)) + self.log("\n") + self.log.parformat("Total livetime") + self.log("%.2f s = %.2f min = %.2f h" % (livetime, livetime/60.0, livetime/3600.)) + self.log("\n") + + # Return + return + + def save(self): + """ + Save pointings to ds9 region file if required + """ + + # Check if ds9 file is valid + if not self.ds9file == "NONE": + + # Open file + f = open(self.ds9file,"w") + + # Write coordinate system + f.write("fk5\n") + + # Loop over pointings + for i in range(len(self.pnt_ra)): + + # Create string + line = "point(" + line+=str(self.pnt_ra[i])+","+str(self.pnt_dec[i])+")" + line+=" # point=cross 20 width=3 \n" + + # Write to file + f.write(line) + + # Close file + f.close() + + # Return + return + + + def zeniths(self): + """ + Return zenith angles + """ + return self.zeniths + + def azimuths(self): + """ + Return azimuth angles + """ + return self.azimuths + + def offsets(self): + """ + Return offset angles + """ + return self.offsets + + def ebounds(self): + """ + Return energy boundaries + """ + return self.ebounds + + def gti(self): + """ + Return good time intervals + """ + return self.gti + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Save ds9 file if required + self.save() + + # Return + return + + +# ======================== # +# 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() + + # Execute application + app.execute() diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py new file mode 100755 index 0000000..df444d9 --- /dev/null +++ b/cscripts/cstsmapmerge.py @@ -0,0 +1,407 @@ +#!/usr/bin/env python +# ========================================================================== +# Merges sliced TS map files into one map +# +# Copyright (C) 2015 Michael Mayer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools +import sys +import glob +import os + + +# ================== # +# cstsmapmerge class # +# ================== # +class cstsmapmerge(ctools.cscript): + """ + This class merges TS map files into one TS map + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "cstsmapmerge" + self.version = "1.1.0" + + # Make sure that parfile exists + file = self.parfile() + + # Initialise application + if len(argv) == 0: + ctools.cscript.__init__(self, self.name, self.version) + elif len(argv) == 1: + ctools.cscript.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + + # Set logger properties + self.log_header() + self.log.date(True) + + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Return + return + + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a + parfile. + """ + + # Set parfile name + parfile = self.name+".par" + + try: + pars = gammalib.GApplicationPars(parfile) + except: + # Signal 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): + """ + Get parameters from parfile and setup the observation. + """ + # Get input models + self.inmodels = self["inmaps"].string() + + # Interpet inmodels string + # Handle ascii files + if "@" == self.inmodels[0]: + self.files = open(self.inmodels.replace("@","")).read().splitlines() + + # Handle wild card strings + elif "*" in self.inmodels: + self.files = glob.glob(self.inmodels) + + # Handle space-separated list + elif " " in self.inmodels: + self.files = self.inmodels.split(" ") + + # Handle comma-separated list + elif "," in self.inmodels: + self.files = self.inmodels.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" + raise gammalib.GException.invalid_argument("cstsmapmerge::get_parameters", msg) + + # Check number of files. We need at least two files. + if len(self.files) <= 1: + msg = "Need at least two files to start merging, "+\ + str(len(self.files))+" files given." + raise gammalib.GException.invalid_argument("cstsmapmerge::get_parameters", msg) + + # Get output file + self.outmap = self["outmap"].filename() + + # Get other parameters + self.overwrite = self["overwrite"].boolean() + self.delete = self["delete"].boolean() + + # Return + return + + def init_map(self, fitsfile): + """ + """ + # Set filename + 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"]) + + # Get other maps + self.maps = [] + self.mapnames = [] + + # Loop over extensions + for hdu in fits: + + # Leave out primary and status extension + if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP": + + # Add present maps + skymap = gammalib.GSkyMap() + skymap.read(hdu) + self.maps.append(skymap) + self.mapnames.append(hdu.extname()) + + # Close FITS file + fits.close() + + # 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) + add_tsmap = gammalib.GSkyMap() + add_tsmap.read(fits[0]) + add_statusmap = gammalib.GSkyMap() + add_statusmap.read(fits["STATUS MAP"]) + + # Get other maps + add_maps = [] + + # Loop over extensions + for hdu in fits: + + # Leave out primary and status extension + if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP": + + # Add present maps + skymap = gammalib.GSkyMap() + skymap.read(hdu) + add_maps.append(skymap) + + # 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+"\"") + + # Loop over bins + for i in range(self.tsmap.npix()): + 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 gammalib.GException.invalid_value("cstsmapmerge::add", msg) + + # Copy TS values + self.tsmap[i] = add_tsmap[i] + + # Copy status + 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] + + # 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") + + if self.logTerse(): + self.log("\n") + self.log.header1("Merging sliced TS map files") + self.log("\n") + + # Initialise file to start with + file0 = "" + + self.added_files = [] + + # Test files for the entry status map + # use the first one to appear useful + for fitsfile in self.files: + if not gammalib.is_fits(fitsfile): + self.log("Skipping file \""+fitsfile+" (not a FITS file)\n") + fits = gammalib.GFits(fitsfile) + + if fits.contains("STATUS MAP"): + self.log("Using file \""+fitsfile+"\" as initialising map\n") + file0 = fitsfile + break + else: + if self.logExplicit(): + self.log("File \""+fitsfile+"\" has not \"STATUS MAP\" extension\n") + + # Print an error message if no suitable file was provided + 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") + + # Otherwise merge files + else: + + # 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) + + # Append to added files + self.added_files.append(file0) + + # Loop over files + for fitsfile in workfiles: + + # Skip if file is not FITS + if not gammalib.is_fits(fitsfile): + if self.logTerse(): + self.log("Skipping file \""+fitsfile+" (not a FITS file)\n") + continue + + # 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") + continue + + # Logging + if self.logTerse(): + self.log("Adding TS map slice \""+fitsfile+"\"\n") + + # Add file + self.add(fitsfile) + + # Append to added files + self.added_files.append(fitsfile) + + # Return + return + + def save(self): + """ + Save merged tsmap file and remove slices if requested + """ + + # Logging + if self.logExplicit(): + self.log("\n") + self.log.header1("Saving merged TS map") + self.log("\n") + + # Create FITS file + fits = gammalib.GFits() + + # Write TS map into primary + self.tsmap.write(fits) + + # Loop over maps and write them to 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]) + + # Check if map is fully done + done = True + for pix in self.statusmap: + if pix < 0.5: + done = False + break + + # Write Status map if we are not done yet + if not done: + 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) + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Save ds9 file if required + self.save() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + """ + Merges ts map slices. + """ + # Create instance of application + app = cstsmapmerge(sys.argv) + + # Open logfile + app.logFileOpen() + + # Execute application + app.execute() diff --git a/cscripts/tsmerge.py b/cscripts/tsmerge.py deleted file mode 100644 index b09f099..0000000 --- a/cscripts/tsmerge.py +++ /dev/null @@ -1,181 +0,0 @@ -#!/usr/bin/env python - -import gammalib -import sys -import os -import glob -import json - -def print_help(): - print("") - print("tsmerge.py expects two keyword argumens: use via command line as follows") - print("tsmerge.py files= outfile=") - print("-------------------------------------------------------------------------") - print("*** \"outfile\" should be a path to a fits file") - print("*** \"files\" should be either be:") - print(" - an ascii-file listing the files your want to merge") - print(" OR - a wildcard string describing the files for merging") - print("-------------------------------------------------------------------------") - print("") - print("Example: python tsmerge,py files=tsmap*.fits outfile=tsmap.fits") - print("") - sys.exit(-1) - - -# Maps class holding the different sky maps -# created by cttsmap -class maps: - - # Constructor - def __init__(self,fitsfile): - - self.filename = fitsfile - if not os.path.isfile(fitsfile): - self.error("Maps instance cannot be constructed: File \""+self.filename+"\" does not exist") - - # open FITS file - fits = gammalib.GFits(fitsfile) - - # get TS map - self.tsmap = gammalib.GSkyMap() - self.tsmap.read(fits[0]) - - # Get status map - if not fits.contains("STATUS MAP"): - self.error("No extension \"STATUS MAP\" found in \""+self.filename+"\"") - self.statusmap = gammalib.GSkyMap() - self.statusmap.read(fits["STATUS MAP"]) - - # Get other maps - self.maps = [] - self.mapnames = [] - for hdu in fits: - if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP": - skymap = gammalib.GSkyMap() - skymap.read(hdu) - self.maps.append(skymap) - self.mapnames.append(hdu.extname()) - - # len operator - def __len__(self): - return len(self.maps) - - # operator [] - def __getitem__(self,index): - return self.maps[index] - - # Check the status map for completeness - def check_status(self): - check = True - for pix in self.statusmap: - if pix < 0.5: - check = False - break - return check - - # Save maps to outfile - def save(self,outfile): - - fits = gammalib.GFits() - self.tsmap.write(fits) - - for i in range(len(self.maps)): - self.maps[i].write(fits) - - for i in range(len(self.mapnames)): - fits[i+1].extname(self.mapnames[i]) - - if not self.check_status(): - self.statusmap.write(fits) - fits[fits.size()-1].extname("STATUS MAP") - fits.saveto(outfile,True) - - # Add map instances to each other - # in this way the TS maps get merged - def add(self,addmaps): - if not len(addmaps) == len(self.maps): - self.error("Maps cannot be added: Maps from \""+self.filename+"\" do not match from \""+addmaps.filename+"\"") - - # Loop over bins - for i in range(self.tsmap.npix()): - if addmaps.statusmap[i] > 0.5: - self.tsmap[i] = addmaps.tsmap[i] - self.statusmap[i] = addmaps.statusmap[i] - for j in range(len(self.maps)): - self.maps[j][i] = addmaps[j][i] - # Unify error messages - def error(self,message): - print("*** ERROR ***: "+message) - sys.exit(-1) - - -def pars_from_argv(argv): - - if argv is None: - argv = sys.argv - - args = [] - kwargs = {} - for s in argv[1:]: - if s.count('=') == 1: - key, value = s.split('=', 1) - else: - key, value = None, s - try: - value = json.loads(value) - except ValueError: - pass - if key: - kwargs[key] = value - else: - args.append(value) - return args, kwargs - - -def get_kwargs(argv): - kwargs = {} - - for s in argv: - if s.count('=') == 1: - key, value = s.split('=', 1) - kwargs[key] = value - - else: - print_help() - - return kwargs - -if __name__ == "__main__": - # Initialise list of files - - kwargs = get_kwargs(sys.argv[1:]) - if "outfile" not in kwargs or "files" not in kwargs: - print_help() - outfile = kwargs["outfile"] - files = kwargs["files"] - allfiles = [] - if os.path.isfile(files): - allfiles = open(files).read().splitlines() - else: - allfiles = glob.glob(files) - - # Initialise files to be merged - mergefiles = [] - - # Test files for the entry status map - # and add them to the file list - for fitsfile in allfiles: - fits = gammalib.GFits(fitsfile) - if fits.contains("STATUS MAP"): - mergefiles.append(fitsfile) - fits.close() - - # create a maps instance from the first file - finalmaps = maps(mergefiles[0]) - - # add other files to the maps instance - for fitsfile in mergefiles[1:]: - finalmaps.add(maps(fitsfile)) - - # save the filled maps to file - finalmaps.save(outfile=outfile) diff --git a/ctools.pc.in b/ctools.pc.in index f150c20..340a402 100644 --- a/ctools.pc.in +++ b/ctools.pc.in @@ -6,9 +6,9 @@ includedir=${prefix}/include Name: ctools Description: CTA science analysis tools URL: http://cta.irap.omp.eu/ctools -Version: 1.0.0 +Version: 1.1.0 Requires: -Requires.private: gammalib >= 1.0.0 +Requires.private: gammalib >= 1.1.0 Conflicts: Cflags: -I${includedir}/ctools @OPENMP_CXXFLAGS@ Libs: -L${libdir} -lctools diff --git a/doc/Doxyfile b/doc/Doxyfile index 5e00812..9f07c4f 100644 --- a/doc/Doxyfile +++ b/doc/Doxyfile @@ -32,7 +32,7 @@ PROJECT_NAME = ctools # This could be handy for archiving the generated documentation or # if some version control system is used. -PROJECT_NUMBER = 1.0.0 +PROJECT_NUMBER = 1.1.0.dev # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer diff --git a/doc/source/conf.py b/doc/source/conf.py index 0288704..f4a4890 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -66,10 +66,10 @@ copyright = u'2014-2015, ctools developers' # built documents. # # The short X.Y version. -version = '1.0' +version = '1.1' # The full version, including alpha/beta/rc tags. -release = '1.0.0' +release = '1.1.0.dev' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/reference_manual/cscaldb.rst b/doc/source/reference_manual/cscaldb.rst index 536d98c..bfbe18e 100644 --- a/doc/source/reference_manual/cscaldb.rst +++ b/doc/source/reference_manual/cscaldb.rst @@ -1 +1 @@ -.. _cscaldb: cscaldb ======= Dumps all available instrument response functions into log fil. Synopsis -------- This scripts dumps all instrument response functions that are available in the calibration database into the log file. If you want the information on the screen, add the ``debug=yes`` parameter. General parameters ------------------ None Standard parameters ------------------- ``(chatter = 2) [integer]`` Verbosity of the executable: chatter = 0: no information will be logged chatter = 1: only errors will be logged chatter = 2: errors and actions will be logged chatter = 3: report about the task execution chatter = 4: detailed report about the task execution ``(clobber = yes) [boolean]`` Specifies whether an existing output counts cube should be overwritten. ``(debug = no) [boolean]`` Enables debug mode. In debug mode the executable will dump any log file output to the console. ``(mode = ql) [string]`` Mode of automatic parameters (default is "ql", i.e. "query and learn"). Related tools ------------- None \ No newline at end of file +.. _cscaldb: cscaldb ======= Dumps all available instrument response functions into log file. Synopsis -------- This scripts dumps all instrument response functions that are available in the calibration database into the log file. If you want the information on the screen, add the ``debug=yes`` parameter. General parameters ------------------ None Standard parameters ------------------- ``(chatter = 2) [integer]`` Verbosity of the executable: chatter = 0: no information will be logged chatter = 1: only errors will be logged chatter = 2: errors and actions will be logged chatter = 3: report about the task execution chatter = 4: detailed report about the task execution ``(clobber = yes) [boolean]`` Specifies whether an existing output counts cube should be overwritten. ``(debug = no) [boolean]`` Enables debug mode. In debug mode the executable will dump any log file output to the console. ``(mode = ql) [string]`` Mode of automatic parameters (default is "ql", i.e. "query and learn"). Related tools ------------- None \ No newline at end of file diff --git a/doc/source/reference_manual/csmodelinfo.rst b/doc/source/reference_manual/csmodelinfo.rst new file mode 100644 index 0000000..7baa7e9 --- /dev/null +++ b/doc/source/reference_manual/csmodelinfo.rst @@ -0,0 +1,101 @@ +.. _csmodelinfo: + +csmodelinfo +=========== + +Dumps information about a model container into log file or on screen. + + +Synopsis +-------- + +This script provides detailled information about a given model container. +It presents the number of parameters, the number of free parameters (sub- +divided by spectral, spatial and temporal type). The tool further creates +information about any parameter that is stuck at the parameter boundaries. +This might be quite useful for debugging a fit. Finally the tool also list +the TS values of each source that was fitted, in order have a quick overview +on these quantities on the screen. + +In addition the tool is also capable of exporting the input model XML container +to a region file which is compatible to be read with ds9 or other displying +tools. To create a region file, the hidden parameter "ds9file" needs to be +specified. There are quite a few further hidden parameters to steer e.g. +color, font, etc of the region file output. + + +General parameters +------------------ + +``inmodel [file]`` + Event list, counts cube, or observation definition file + +``(ds9file = NONE) [file]`` + Output DS9 region file containing models + +``(pnt_type=cross) [string]`` + Marker type for point sources + +``pnt_mark_size [integer]`` + Marker size for point sources + +``(show_labels=yes) [boolean]`` + Add source labels? + +``(width=2) [integer]`` + Line width for regions + + ``(font=helvetica) [string]`` + Font for source labels + +``(fontsize=12) [integer]`` + Font size for source labels + +``(fonttype=normal) [string]`` + Use normal or bold font? + +``(fonttype2=roman) [string]`` + Use roman or italic font? + +``(show_ext_type=yes) [boolean]`` + Show type of extended model in source name? + +``(free_color=green) [string]`` + Color for sources with free parameters (any ds9 color or hex code) + +``(fixed_color=magenta) [string]`` + Color for sources without free parameters (any ds9 color or hex code) + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output counts cube should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = csmodelinfo.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +None \ No newline at end of file diff --git a/doc/source/reference_manual/csmodelmerge.rst b/doc/source/reference_manual/csmodelmerge.rst new file mode 100644 index 0000000..cb92c09 --- /dev/null +++ b/doc/source/reference_manual/csmodelmerge.rst @@ -0,0 +1,67 @@ +.. _csmodelmerge: + +csmodelmerge +============ + +Merges an arbitrary number of model XML files into a single file. + + +Synopsis +-------- + +This script merges several model XML file into a single model XML file. The +files are passed via the ``inmodels`` parameter. There are 4 options to specify +a list of files that should be merged: +- a space-separated list of file names (e.g. model1.xml model2.xml) +- a comma-separated list of file names (e.g. model1.xml,model2.xml) +- a wildcard string (e.g. mymodel*.xml) +- an ASCII file containing the file names, one file per line (e.g. @mymodels.txt) + +Note that the "@" needs to be specified in case of an ASCII file name. It is +also important to note that models with the same name cannot be appended. The +tool will throw and exception if the provided file contain models with the same +name. + + +General parameters +------------------ + +``inmodels [string]`` + Input model XML files + +``outmodel [file]`` + Output model file + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output counts cube should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = csmodelmerge.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +None \ No newline at end of file diff --git a/doc/source/reference_manual/csobsinfo.rst b/doc/source/reference_manual/csobsinfo.rst new file mode 100644 index 0000000..6f47bd3 --- /dev/null +++ b/doc/source/reference_manual/csobsinfo.rst @@ -0,0 +1,88 @@ +.. _csobsinfo: + +csobsinfo +========= + +Dumps information about an observation container into log file or on +screen. + + +Synopsis +-------- + +This script provides detailed information about a given observation container. +For IACT observations usually several observations with different configurations +are combined. The script computes e.g. the overall energy range which could be +an important input for :doc:`ctbin` to create a suitable energy range. In +addition the script prints the pointing attributes like the average +zenith/azimuth angle. If the hidden boolean parameter "offset" is specified to +"yes", the script queries for the target position and computes the run-by-run +offset as well as its average. In the observation summary, the time range of +the observations including total ontime and (dead-time corrected) livetime is +shown. + +The script further contains the hidden parameter "ds9file". If specified with +a file name, the tool writes out a ds9 region file including a cross for each +pointing position in the sky. + +When executed from within python, the script provides methods to return arrays +containing observation-wise parameters. This might be useful for graphical +display within a script. The following methods are implemented: +- zeniths (returning all zenith angles) +- azimuths (returning all azimuth angles) +- offsets (returning all offset angles, provided the hidden parameter offset=yes) +- ebounds (returning energy boundaries) +- gtis (returning good time intervals) + + +General parameters +------------------ + +``inobs [file]`` + Event list, counts cube, or observation definition file + +``(offset = no) [file]`` + Compute offset from target to pointing positions + +``ra [real]`` + Target right ascension for offset computation (only queried if offset=yes) + +``dec [real]`` + Target declination for offset computation (only queried if offset=yes) + +``(ds9file = NONE) [file]`` + DS9 region file containing pointings + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output counts cube should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = csobsinfo.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +None \ No newline at end of file diff --git a/doc/source/reference_manual/cstsmapmerge.rst b/doc/source/reference_manual/cstsmapmerge.rst new file mode 100644 index 0000000..1370f10 --- /dev/null +++ b/doc/source/reference_manual/cstsmapmerge.rst @@ -0,0 +1,72 @@ +.. _cstsmapmerge: + +cstsmapmerge +============ + +Merges slices of splitted tsmap computations. + + +Synopsis +-------- + +In case a TS map has been computed using the option to split the calculations +into several jobs, this tool is required to merge the output into one final +map again. Running :doc:`cttsmap` with the hidden parameters ``binmin`` and +``binmax`` specified will result in several sliced TS maps. This script merges +them into one FITS file. The input slices can be passed via the ``inmaps`` +parameter. There are 4 options to specify a list of files that should be merged: +- a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits) +- a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits) +- a wildcard string (e.g. tsmap*.xml) +- an ASCII file containing the file names, one file per line (e.g. @mymaps.txt) + +Note that the "@" needs to be specified in case of an ASCII file name. + +General parameters +------------------ + +``inmaps [string]`` + Input TS map FITS files + +``outmap [file]`` + Output TS map FITS file + +``(overwrite=yes) [boolean]`` + Overwrite previously filled bins? + +``(delete=no) [boolean]`` + Delete merged files? + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output counts cube should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = cstsmapmerge.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +:ref:`cstsmap` diff --git a/doc/source/reference_manual/cttsmap.rst b/doc/source/reference_manual/cttsmap.rst index 9789cfc..6c077e8 100644 --- a/doc/source/reference_manual/cttsmap.rst +++ b/doc/source/reference_manual/cttsmap.rst @@ -108,7 +108,8 @@ Standard parameters Name of log file. -Related tools -------------- +Related tools or scripts +------------------------ + +:ref:`cstsmapmerge` -None diff --git a/doc/source/reference_manual/index.rst b/doc/source/reference_manual/index.rst index 90aacd6..407b27f 100644 --- a/doc/source/reference_manual/index.rst +++ b/doc/source/reference_manual/index.rst @@ -45,3 +45,7 @@ cscripts csspec --- Computes spectral points csresmap --- Generates residual map cstsdist --- Generates TS distribution + csobsinfo --- Shows observation container content + csmodelinfo --- Shows model container content + csmodelmerge --- Merges several model containers into one file + cstsmapmerge --- Merges slices from ts map computations diff --git a/test/data/model_cube_background2.xml b/test/data/model_cube_background2.xml new file mode 100644 index 0000000..82e839b --- /dev/null +++ b/test/data/model_cube_background2.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/test/data/tsmap_0.fits b/test/data/tsmap_0.fits new file mode 100644 index 0000000..57f530b --- /dev/null +++ b/test/data/tsmap_0.fits diff --git a/test/data/tsmap_1.fits b/test/data/tsmap_1.fits new file mode 100644 index 0000000..e905c5b --- /dev/null +++ b/test/data/tsmap_1.fits diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index 628c16e..48691bf 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -36,7 +36,7 @@ $ECHO "*****************" # # Remove any existing result files # ================================ -rm -rf *.fits *.xml *.dat +rm -rf *.fits *.xml *.dat *.reg # @@ -74,6 +74,203 @@ $ECHO " ok" # +# Test csobsinfo +# ============== +$ECHO -n "Test csobsinfo: " +chmod +x @abs_top_srcdir@/cscripts/csobsinfo.py +@abs_top_srcdir@/cscripts/csobsinfo.py \ + inobs="obs.xml" \ + ds9file="obs.reg" \ + offset=yes \ + ra="83.6331" \ + dec="22.01" > csobsinfo.log +$ECHO -n "." +if [ -s "obs.reg" ] +then + $ECHO -n "." +else + $ECHO " obs.reg file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test csmodelinfo +# ================ +$ECHO -n "Test csmodelinfo: " +chmod +x @abs_top_srcdir@/cscripts/csmodelinfo.py +@abs_top_srcdir@/cscripts/csmodelinfo.py \ + inmodel="data/crab.xml" \ + ds9file="model.reg" > csmodelinfo.log +$ECHO -n "." +if [ -s "model.reg" ] +then + $ECHO -n "." +else + $ECHO " model.reg file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test csmodelmerge (space-separated) +# =================================== +$ECHO -n "Test csmodelmerge (space-separated): " +chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py +@abs_top_srcdir@/cscripts/csmodelmerge.py \ + inmodels="data/crab.xml data/source1.xml" \ + outmodel="merged_models.xml" > csmodelmerge.log +$ECHO -n "." +if [ -s "merged_models.xml" ] +then + $ECHO -n "." +else + $ECHO " merged_models.xml file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test csmodelmerge (comma-separated) +# =================================== +$ECHO -n "Test csmodelmerge (comma-separated): " +chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py +@abs_top_srcdir@/cscripts/csmodelmerge.py \ + inmodels="data/crab.xml,data/source1.xml" \ + outmodel="merged_models2.xml" > csmodelmerge.log +$ECHO -n "." +if [ -s "merged_models2.xml" ] +then + $ECHO -n "." +else + $ECHO " merged_models2.xml file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test csmodelmerge (wildcard string) +# =================================== +$ECHO -n "Test csmodelmerge (wildcard string): " +chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py +@abs_top_srcdir@/cscripts/csmodelmerge.py \ + inmodels="data/model_cube_background*.xml" \ + outmodel="merged_models3.xml" > csmodelmerge.log +$ECHO -n "." +if [ -s "merged_models3.xml" ] +then + $ECHO -n "." +else + $ECHO " merged_models3.xml file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test csmodelmerge (ascii-file) +# ============================== +ls data/model_cube_background* > models.txt +$ECHO -n "Test csmodelmerge (ascii-file): " +chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py +@abs_top_srcdir@/cscripts/csmodelmerge.py \ + inmodels="@models.txt" \ + outmodel="merged_models4.xml" > csmodelmerge.log +$ECHO -n "." +if [ -s "merged_models4.xml" ] +then + $ECHO -n "." +else + $ECHO " merged_models4.xml file is not valid" + exit 1 +fi +$ECHO " ok" +rm models.txt + + +# +# Test cstsmapmerge (comma-separated) +# =================================== +$ECHO -n "Test cstsmapmerge (comma-separated): " +chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py +@abs_top_srcdir@/cscripts/cstsmapmerge.py \ + inmaps="data/tsmap_0.fits,data/tsmap_1.fits" \ + outmap="merged_tsmap1.fits" > cstsmapmerge.log +$ECHO -n "." +if [ -s "merged_tsmap1.fits" ] +then + $ECHO -n "." +else + $ECHO " merged_tsmap1.fits file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test cstsmapmerge (space-separated) +# =================================== +$ECHO -n "Test cstsmapmerge (space-separated): " +chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py +@abs_top_srcdir@/cscripts/cstsmapmerge.py \ + inmaps="data/tsmap_0.fits data/tsmap_1.fits" \ + outmap="merged_tsmap2.fits" > cstsmapmerge.log +$ECHO -n "." +if [ -s "merged_tsmap2.fits" ] +then + $ECHO -n "." +else + $ECHO " merged_tsmap2.fits file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test cstsmapmerge (wildcard string) +# =================================== +$ECHO -n "Test cstsmapmerge (wildcard string): " +chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py +@abs_top_srcdir@/cscripts/cstsmapmerge.py \ + inmaps="data/tsmap_*.fits" \ + outmap="merged_tsmap3.fits" > cstsmapmerge.log +$ECHO -n "." +if [ -s "merged_tsmap3.fits" ] +then + $ECHO -n "." +else + $ECHO " merged_tsmap3.fits file is not valid" + exit 1 +fi +$ECHO " ok" + + +# +# Test cstsmapmerge (ascii-file) +# ============================== +ls data/tsmap_*.fits > tsmaps.txt +$ECHO -n "Test cstsmapmerge (ascii-file): " +chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py +@abs_top_srcdir@/cscripts/cstsmapmerge.py \ + inmaps="@tsmaps.txt" \ + outmap="merged_tsmap4.fits" > cstsmapmerge.log +$ECHO -n "." +if [ -s "merged_tsmap4.fits" ] +then + $ECHO -n "." +else + $ECHO " merged_tsmap4.fits file is not valid" + exit 1 +fi +$ECHO " ok" +rm tsmaps.txt + + +# # Test cspull # =========== $ECHO -n "Test cspull: "