diff --git a/pyext/ctool.i b/pyext/ctool.i index 9923689..c98f8d8 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Juergen Knoedlseder * + * copyright (C) 2014-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -25,6 +25,7 @@ */ %{ /* Put headers and other declarations here that are needed for compilation */ +#define SWIG #include "ctool.hpp" %} @@ -51,4 +52,26 @@ public: // Public methods virtual void execute(void); + + // Protected methods + void init_members(void); + void copy_members(const ctool& app); + void free_members(void); + const bool& read_ahead(void) const; + GObservations get_observations(const bool& get_response = true); + + // Protected methods that create objects from user parameters + GEbounds create_ebounds(void); + GSkymap create_map(const GObservations& obs); + GCTAEventCube create_cube(const GObservations& obs); + GCTAObservation create_cta_obs(void); + + // Protected methods that check user parameters + void require_inobs(const std::string& method); + + // Protected support methods + void set_response(GObservations& obs); + void set_obs_response(GCTAObservation* obs); + GSkyDir get_mean_pointing(const GObservations& obs); + size_t get_current_rss(void); }; diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py index e5659c8..cd0998d 100755 --- a/pyext/ctools/obsutils.py +++ b/pyext/ctools/obsutils.py @@ -42,6 +42,7 @@ def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, binsz - Pixel size for binned simulation (deg/pixel) npix - Number of pixels in X and Y for binned simulation """ + # Allocate ctobssim application and set parameters sim = ctools.ctobssim(obs) sim["seed"].integer(seed) diff --git a/scripts/cspull.py b/scripts/cspull.py index 66fdd5d..114512d 100755 --- a/scripts/cspull.py +++ b/scripts/cspull.py @@ -2,7 +2,7 @@ # ========================================================================== # This script generates the pull distribution for all model parameters. # -# Copyright (C) 2011-2014 Juergen Knoedlseder +# Copyright (C) 2011-2015 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -22,7 +22,6 @@ import gammalib from ctools import obsutils import sys import csv -import math # ============ # @@ -41,12 +40,12 @@ class cspull(gammalib.GApplication): """ # Set name self.name = "cspull" - self.version = "0.3.0" + self.version = "0.4.0" # Initialise some members - self.obs = None - self.model = None - self.m_srcmdl = None + self.obs = None + self.model = None + self.m_inmodel = None # Make sure that parfile exists file = self.parfile() @@ -93,7 +92,7 @@ class cspull(gammalib.GApplication): # Create default parfile pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model")) + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model")) pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name")) pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database")) @@ -123,7 +122,7 @@ class cspull(gammalib.GApplication): """ # Get parameters if self.model == None: - self.m_srcmdl = self["srcmdl"].filename() + self.m_inmodel = self["inmodel"].filename() self.m_outfile = self["outfile"].filename() self.m_ntrials = self["ntrials"].integer() self.m_caldb = self["caldb"].string() @@ -144,15 +143,15 @@ class cspull(gammalib.GApplication): # Set some fixed parameters self.m_log = False # Logging in client tools - self.m_debug = False # Debugging in client tools + self.m_debug = True # Debugging in client tools # Setup observations self.obs = self.set_obs() - + # Load source model - if self.m_srcmdl != None: - self.model = gammalib.GModels(self.m_srcmdl) - + if self.m_inmodel != None: + self.model = gammalib.GModels(self.m_inmodel) + # Append source model to observation self.obs.models(self.model) diff --git a/scripts/csresmap.py b/scripts/csresmap.py index b4bb983..dcb8564 100755 --- a/scripts/csresmap.py +++ b/scripts/csresmap.py @@ -21,10 +21,7 @@ import ctools import gammalib import sys -import tempfile import os -import math - # ============== # # cstsdist class # @@ -42,7 +39,7 @@ class csresmap(gammalib.GApplication): """ # Set name self.name = "csresmap" - self.version = "0.1.0" + self.version = "0.2.0" # Initialise some members self.obs = None @@ -103,9 +100,9 @@ class csresmap(gammalib.GApplication): # Create default parfile pars = gammalib.GApplicationPars() - pars.append(gammalib.GApplicationPar("infile","f","a","obs.xml","","","Observation definition")) - pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab.xml","","","model file name")) - pars.append(gammalib.GApplicationPar("outfile","f","a","resmap.fits","","","Output file name")) + pars.append(gammalib.GApplicationPar("inobs","f","a","obs.xml","","","Observation definition")) + pars.append(gammalib.GApplicationPar("inmodel","f","a","crab.xml","","","model file name")) + pars.append(gammalib.GApplicationPar("outmap","f","a","resmap.fits","","","Output file name")) pars.append(gammalib.GApplicationPar("binfile","f","h","binned.fits","","","Output binned file name")) pars.append(gammalib.GApplicationPar("modfile","f","h","model.fits","","","Output model file name")) pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) @@ -137,7 +134,7 @@ class csresmap(gammalib.GApplication): # Set observation if not done before if self.obs.size() == 0: - obsfile = self["infile"].filename() + obsfile = self["inobs"].filename() try: self.obs = gammalib.GObservations(obsfile) @@ -163,8 +160,8 @@ class csresmap(gammalib.GApplication): # Note that the models should be optimised for # meaningful residual map if self.obs.models().size() == 0: - self.obs.models(self["srcmdl"].filename()) - self.m_outfile = self["outfile"].filename() + self.obs.models(self["inmodel"].filename()) + self.m_outfile = self["outmap"].filename() self.m_modfile = self["modfile"].filename() self.m_binfile = self["binfile"].filename() self.m_xref = self["xref"].real() diff --git a/scripts/cssens.py b/scripts/cssens.py index 6798198..abd82a0 100755 --- a/scripts/cssens.py +++ b/scripts/cssens.py @@ -17,7 +17,7 @@ # # The source is modelled as a point source # -# Copyright (C) 2011-2014 Juergen Knoedlseder +# Copyright (C) 2011-2015 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -33,8 +33,8 @@ # along with this program. If not, see . # # ========================================================================== -from ctools import * -from gammalib import * +import gammalib +from ctools import obsutils import sys import csv import math @@ -43,7 +43,7 @@ import math # ============ # # cssens class # # ============ # -class cssens(GApplication): +class cssens(gammalib.GApplication): """ This class implements the sensitivity computation script. It derives from the GammaLib::GApplication class which provides support for parameter @@ -63,9 +63,9 @@ class cssens(GApplication): # Initialise application if len(argv) == 0: - GApplication.__init__(self, self.name, self.version) + gammalib.GApplication.__init__(self, self.name, self.version) elif len(argv) ==1: - GApplication.__init__(self, self.name, self.version, *argv) + gammalib.GApplication.__init__(self, self.name, self.version, *argv) else: raise TypeError("Invalid number of arguments given.") @@ -97,31 +97,31 @@ class cssens(GApplication): parfile = self.name+".par" try: - pars = GApplicationPars(parfile) + 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 = GApplicationPars() - pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name")) - pars.append(GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database")) - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) - pars.append(GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) - pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)")) - pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)")) - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) - pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold")) - pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?")) - pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index")) - pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius")) - pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width")) - pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations")) - pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average")) + pars = gammalib.GApplicationPars() + pars.append(gammalib.GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name")) + pars.append(gammalib.GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database")) + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) + pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) + pars.append(gammalib.GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) + pars.append(gammalib.GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) + 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("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("sigma","r","h","5.0","","","Significance threshold")) + pars.append(gammalib.GApplicationPar("ts_use","b","h","yes","","","Use TS?")) + pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Assumed spectral index")) + pars.append(gammalib.GApplicationPar("radius","r","h","0.1","","","Extended source model radius")) + pars.append(gammalib.GApplicationPar("width","r","h","0.05","","","Extended source model width")) + 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.save(parfile) @@ -192,8 +192,8 @@ class cssens(GApplication): # Initialise models. Note that we centre the point source at the Galactic # center as our observation is also centred at the Galactic centre, so # we're onaxis. - bkg_model = GModels() - full_model = GModels() + bkg_model = gammalib.GModels() + full_model = gammalib.GModels() bkg_model.append(self.set_bkg_model()) full_model.append(self.set_bkg_model()) full_model.append(self.set_src_model(0.0, \ @@ -263,10 +263,10 @@ class cssens(GApplication): emax - Maximum energy [TeV] (default: 100.0) """ # Allocate observation container - obs = GObservations() + obs = gammalib.GObservations() # Set single pointing - pntdir = GSkyDir() + pntdir = gammalib.GSkyDir() pntdir.lb_deg(lpnt, bpnt) # Create CTA observation @@ -289,21 +289,21 @@ class cssens(GApplication): fitidx - Fit spectral index (default: False) """ # Define radial component - radial = GCTAModelRadialGauss(3.0) + radial = gammalib.GCTAModelRadialGauss(3.0) # Define spectral component. If a background model is given then # use a file function. Otherwise use a static power law. if len(self.m_bkg) > 0: - spectrum = GModelSpectralFunc(self.m_bkg, 1.0) + spectrum = gammalib.GModelSpectralFunc(self.m_bkg, 1.0) else: - spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV")) + spectrum = gammalib.GModelSpectralPlaw(61.8e-6, -1.85, gammalib.GEnergy(1.0, "TeV")) if fitidx: spectrum["Index"].free() else: spectrum["Index"].fix() # Create background model - model = GCTAModelRadialAcceptance(radial, spectrum) + model = gammalib.GCTAModelRadialAcceptance(radial, spectrum) model.name("Background") model.instruments("CTA") @@ -328,11 +328,11 @@ class cssens(GApplication): fitidx - Fit index? (default: False) """ # Set source location - location = GSkyDir() + location = gammalib.GSkyDir() location.lb_deg(l, b) # Set source spectrum - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")) + spectrum = gammalib.GModelSpectralPlaw(flux*5.7e-16, index, gammalib.GEnergy(0.3, "TeV")) if fitidx: spectrum["Index"].free() else: @@ -340,7 +340,7 @@ class cssens(GApplication): # Set source if type == "point": - spatial = GModelSpatialPointSource(location) + spatial = gammalib.GModelSpatialPointSource(location) if fitpos: spatial["RA"].free() spatial["DEC"].free() @@ -348,7 +348,7 @@ class cssens(GApplication): spatial["RA"].fix() spatial["DEC"].fix() elif type == "gauss": - spatial = GModelSpatialRadialGauss(location, self.m_radius) + spatial = gammalib.GModelSpatialRadialGauss(location, self.m_radius) if fitpos: spatial["RA"].free() spatial["DEC"].free() @@ -358,7 +358,7 @@ class cssens(GApplication): spatial["DEC"].fix() spatial["Sigma"].fix() elif type == "disk": - spatial = GModelSpatialRadialDisk(location, self.m_radius) + spatial = gammalib.GModelSpatialRadialDisk(location, self.m_radius) if fitpos: spatial["RA"].free() spatial["DEC"].free() @@ -368,7 +368,7 @@ class cssens(GApplication): spatial["DEC"].fix() spatial["Radius"].fix() elif type == "shell": - spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width) + spatial = gammalib.GModelSpatialRadialShell(location, self.m_radius, self.m_width) if fitpos: spatial["RA"].free() spatial["DEC"].free() @@ -382,7 +382,7 @@ class cssens(GApplication): else: self.log("ERROR: Unknown source type '"+type+"'.\n") return None - source = GModelSky(spatial, spectrum) + source = gammalib.GModelSky(spatial, spectrum) # Set source name source.name("Test") diff --git a/scripts/cstsdist.py b/scripts/cstsdist.py index 7e8afb8..67dd56b 100755 --- a/scripts/cstsdist.py +++ b/scripts/cstsdist.py @@ -3,7 +3,7 @@ # This script generates the TS distribution for a particular model based # on Monte-Carlo simulations. # -# Copyright (C) 2011-2014 Juergen Knoedlseder +# Copyright (C) 2011-2015 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -19,8 +19,8 @@ # along with this program. If not, see . # # ========================================================================== -from ctools import * -from gammalib import * +import gammalib +from ctools import obsutils import sys import csv import math @@ -29,7 +29,7 @@ import math # ============== # # cstsdist class # # ============== # -class cstsdist(GApplication): +class cstsdist(gammalib.GApplication): """ This class implements the TS distribution generation script. It derives from the GammaLib::GApplication class which provides support for parameter @@ -42,7 +42,7 @@ class cstsdist(GApplication): """ # Set name self.name = "cstsdist" - self.version = "0.1.0" + self.version = "0.2.0" # Initialise some members self.obs = None @@ -54,9 +54,9 @@ class cstsdist(GApplication): # Initialise application if len(argv) == 0: - GApplication.__init__(self, self.name, self.version) + gammalib.GApplication.__init__(self, self.name, self.version) elif len(argv) ==1: - GApplication.__init__(self, self.name, self.version, *argv) + gammalib.GApplication.__init__(self, self.name, self.version, *argv) else: raise TypeError("Invalid number of arguments given.") @@ -87,28 +87,28 @@ class cstsdist(GApplication): parfile = self.name+".par" try: - pars = GApplicationPars(parfile) + 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 = GApplicationPars() - pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name")) - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials")) - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) - pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index")) - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)")) - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)")) - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)")) - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) + pars = gammalib.GApplicationPars() + pars.append(gammalib.GApplicationPar("outfile","f","h","ts.dat","","","Output file name")) + pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) + pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) + pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) + pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Spectral index")) + pars.append(gammalib.GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) + pars.append(gammalib.GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) + pars.append(gammalib.GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)")) + pars.append(gammalib.GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)")) + pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) + pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) + pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)")) + 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_standard() pars.save(parfile) @@ -146,8 +146,8 @@ class cstsdist(GApplication): # Initialise models. Note that we centre the point source at the Galactic # center as our observation is also centred at the Galactic centre, so # we're onaxis. - self.bkg_model = GModels() - self.full_model = GModels() + self.bkg_model = gammalib.GModels() + self.full_model = gammalib.GModels() self.bkg_model.append(self.set_bkg_model()) self.full_model.append(self.set_bkg_model()) self.full_model.append(self.set_src_model(0.0, self.m_offset, \ @@ -247,10 +247,10 @@ class cstsdist(GApplication): emax - Maximum energy [TeV] (default: 100.0) """ # Allocate observation container - obs = GObservations() + obs = gammalib.GObservations() # Set single pointing - pntdir = GSkyDir() + pntdir = gammalib.GSkyDir() pntdir.lb_deg(lpnt, bpnt) # Create CTA observation @@ -269,17 +269,17 @@ class cstsdist(GApplication): Setup CTA background model. """ # Define radial component - radial = GCTAModelRadialGauss(3.0) + radial = gammalib.GCTAModelRadialGauss(3.0) if fitsigma: radial["Sigma"].free() else: radial["Sigma"].fix() # Define spectral component - spectrum = GModelSpectralFunc(self.m_bkg, 1.0) + spectrum = gammalib.GModelSpectralFunc(self.m_bkg, 1.0) # Create background model - model = GCTAModelRadialAcceptance(radial, spectrum) + model = gammalib.GCTAModelRadialAcceptance(radial, spectrum) model.name("Background") model.instruments("CTA") @@ -308,11 +308,11 @@ class cstsdist(GApplication): fitidx - Fit index? (default: True) """ # Set source location - location = GSkyDir() + location = gammalib.GSkyDir() location.lb_deg(l, b) # Set source spectrum - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")) + spectrum = gammalib.GModelSpectralPlaw(flux*5.7e-16, index, gammalib.GEnergy(0.3, "TeV")) if fitidx: spectrum["Index"].free() else: @@ -320,32 +320,33 @@ class cstsdist(GApplication): # Set source if type == "point": - spatial = GModelSpatialPointSource(location) + spatial = gammalib.GModelSpatialPointSource(location) if fitpos: spatial[0].free() spatial[1].free() elif type == "gauss": - spatial = GModelSpatialRadialGauss(location, sigma) + spatial = gammalib.GModelSpatialRadialGauss(location, sigma) if fitpos: spatial[0].free() spatial[1].free() elif type == "disk": - spatial = GModelSpatialRadialDisk(location, radius) + spatial = gammalib.GModelSpatialRadialDisk(location, radius) if fitpos: spatial[0].free() spatial[1].free() elif type == "shell": - spatial = GModelSpatialRadialShell(location, radius, width) + spatial = gammalib.GModelSpatialRadialShell(location, radius, width) if fitpos: spatial[0].free() spatial[1].free() else: self.log("ERROR: Unknown source type '"+type+"'.\n") return None - source = GModelSky(spatial, spectrum) + source = gammalib.GModelSky(spatial, spectrum) # Set source name source.name("Test") + source.tscalc(True) # Return source return source @@ -368,7 +369,7 @@ class cstsdist(GApplication): binsz=self.m_binsz, \ npix=self.m_npix, \ log=self.m_log, debug=self.m_debug) - + # Determine number of events in simulation nevents = 0.0 for run in sim: @@ -380,7 +381,7 @@ class cstsdist(GApplication): self.log.parformat("Number of simulated events") self.log(nevents) self.log("\n") - + # Fit background only sim.models(self.bkg_model) like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index a360bb9..e017889 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -33,8 +33,8 @@ #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ -#define G_INIT_CUBE "ctbin::init_cube()" #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)" +#define G_GET_PARAMETERS "ctbin::get_parameters()" /* __ Debug definitions __________________________________________________ */ @@ -238,9 +238,6 @@ void ctbin::run(void) } } - // Initialise counts cube - init_cube(); - // Loop over all observations in the container for (int i = 0; i < m_obs.size(); ++i) { @@ -309,14 +306,14 @@ void ctbin::save(void) } // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); // Get CTA observation from observation container GCTAObservation* obs = dynamic_cast(m_obs[0]); // Save only if observation is valid if (obs != NULL) { - obs->save(m_outfile, clobber()); + obs->save(m_outcube, clobber()); } // Return @@ -336,16 +333,8 @@ void ctbin::save(void) void ctbin::init_members(void) { // Initialise members - m_evfile.clear(); - m_outfile.clear(); - m_proj.clear(); - m_coordsys.clear(); - m_usepnt = false; - m_xref = 0.0; - m_yref = 0.0; - m_binsz = 0.0; - m_nxpix = 0; - m_nypix = 0; + m_outcube.clear(); + m_usepnt = false; // Initialise protected members m_obs.clear(); @@ -368,24 +357,16 @@ void ctbin::init_members(void) void ctbin::copy_members(const ctbin& app) { // Copy attributes - m_evfile = app.m_evfile; - m_outfile = app.m_outfile; - m_proj = app.m_proj; - m_coordsys = app.m_coordsys; - m_usepnt = app.m_usepnt; - m_xref = app.m_xref; - m_yref = app.m_yref; - m_binsz = app.m_binsz; - m_nxpix = app.m_nxpix; - m_nypix = app.m_nypix; + m_outcube = app.m_outcube; + m_usepnt = app.m_usepnt; // Copy protected members - m_obs = app.m_obs; - m_cube = app.m_cube; - m_ebounds = app.m_ebounds; - m_gti = app.m_gti; - m_ontime = app.m_ontime; - m_livetime = app.m_livetime; + m_obs = app.m_obs; + m_cube = app.m_cube; + m_ebounds = app.m_ebounds; + m_gti = app.m_gti; + m_ontime = app.m_ontime; + m_livetime = app.m_livetime; // Return return; @@ -413,59 +394,32 @@ void ctbin::free_members(void) ***************************************************************************/ void ctbin::get_parameters(void) { - // Initialize energy boundaries - m_ebounds.clear(); - - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Get name of CTA events file - m_evfile = (*this)["evfile"].filename(); - - // Allocate CTA observation - GCTAObservation obs; - - // Try first to open as FITS file - try { + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Load event list in CTA observation - obs.load(m_evfile); + // Get observation container without response (not needed) + m_obs = get_observations(false); - // Append CTA observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(m_evfile); + } // endif: there was no observation in the container - } + // Create an event cube based on task parameters + GCTAEventCube cube = create_cube(m_obs); - } // endif: there was no observation in the container + // Get the skymap from the cube and initialise all pixels to zero + m_cube = cube.map(); + m_cube = 0.0; - // Get parameters - m_usepnt = (*this)["usepnt"].boolean(); - if (!m_usepnt) { - m_xref = (*this)["xref"].real(); - m_yref = (*this)["yref"].real(); - } - m_proj = (*this)["proj"].string(); - m_coordsys = (*this)["coordsys"].string(); - m_binsz = (*this)["binsz"].real(); - m_nxpix = (*this)["nxpix"].integer(); - m_nypix = (*this)["nypix"].integer(); - - // Set energy boundaries - m_ebounds = get_ebounds(); + // Get energy boundaries + m_ebounds = cube.ebounds(); // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); } // Return @@ -474,71 +428,6 @@ void ctbin::get_parameters(void) /***********************************************************************//** - * @brief Initialise counts cube information - * - * @exception GException::invalid_value - * No valid CTA observation found to derive the counts cube - * map centre. - * - * Initialises the skymap, energy boundaries and GTI for a counts cube. - ***************************************************************************/ -void ctbin::init_cube(void) -{ - // Initialse cube information - m_ontime = 0.0; - m_livetime = 0.0; - m_cube.clear(); - m_gti.clear(); - - // Set event cube centre, either from the user parameters or from the - // pointing - double xref = m_xref; - double yref = m_yref; - if (m_usepnt) { - - // Dummy: get pointing from first observation. Ultimately, we want - // to get the pointing from a kind of average - bool found = false; - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - const GCTAPointing& pnt = obs->pointing(); - if (gammalib::toupper(m_coordsys) == "GAL") { - xref = pnt.dir().l_deg(); - yref = pnt.dir().b_deg(); - } - else { - xref = pnt.dir().ra_deg(); - yref = pnt.dir().dec_deg(); - } - found = true; - break; - } - } - - // Signal if no pointing is found - if (!found) { - std::string msg = "No valid CTA observation has been found in " - "observation list, hence no pointing information " - "could be extracted. Use the \"usepnt=no\" " - "option and specify pointing explicitly when " - "running ctbin."; - throw GException::invalid_value(G_INIT_CUBE, msg); - } - - } // endif: used pointing - - // Create skymap - m_cube = GSkymap(m_proj, m_coordsys, - xref, yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, m_ebounds.size()); - - // Return - return; -} - - -/***********************************************************************//** * @brief Fill events into counts cube * * @param[in] obs CTA observation. @@ -590,8 +479,8 @@ void ctbin::fill_cube(GCTAObservation* obs) } // Skip if pixel is out of range - if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) || - pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) { + if (pixel.x() < -0.5 || pixel.x() > (m_cube.nx()-0.5) || + pixel.y() < -0.5 || pixel.y() > (m_cube.ny()-0.5)) { num_outside_map++; continue; } diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index ca72d48..50f46dd 100644 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTBIN_NAME "ctbin" -#define CTBIN_VERSION "00-05-00" +#define CTBIN_VERSION "00-06-00" /***********************************************************************//** @@ -87,21 +87,12 @@ protected: void copy_members(const ctbin& app); void free_members(void); void get_parameters(void); - void init_cube(void); void fill_cube(GCTAObservation* obs); void obs_cube(void); // User parameters - std::string m_evfile; //!< Input event list or XML file - std::string m_outfile; //!< Output counts map or XML file + std::string m_outcube; //!< Output counts map file bool m_usepnt; //!< Use pointing instead of xref/yref parameters - std::string m_proj; //!< WCS projection - std::string m_coordsys; //!< Coordinate system - double m_xref; //!< Longitude reference coordinate - double m_yref; //!< Latitude reference coordinate - double m_binsz; //!< Pixel size - int m_nxpix; //!< Number of pixels in longitude - int m_nypix; //!< Number of pixels in latitude // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctbin/ctbin.par b/src/ctbin/ctbin.par index 4fef2c6..bdcb8d7 100644 --- a/src/ctbin/ctbin.par +++ b/src/ctbin/ctbin.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctbin.par # -# Version: 00-05-00 # +# Version: 00-06-00 # # # # Author: Jurgen Knodlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -17,8 +17,8 @@ # # General parameters #=================== -evfile, f, a, "events.fits",,,"Input event list or observation definition file" -outfile, f, a, "cntmap.fits",,, "Output counts cube" +inobs, f, a, "events.fits",,,"Event list or observation definition file" +outcube, f, a, "cntmap.fits",,, "Output counts cube" # # Energy binning parameters #========================== diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index c7d722b..56d8647 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbkgcube - Background cube generation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Chia-Chun Lu * + * copyright (C) 2014-2015 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -34,8 +34,6 @@ /* __ Method name definitions ____________________________________________ */ #define G_GET_PARAMETERS "ctbkgcube::get_parameters()" -#define G_GET_OBS "ctbkgcube::get_obs()" -#define G_SET_FROM_CNTMAP "ctbkgcube::set_from_cntmap(std::string&)" #define G_FILL_CUBE "ctbkgcube::fill_cube(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ @@ -258,6 +256,7 @@ void ctbkgcube::run(void) // container and put all removed components in the output // container int num = m_bkgmdl.size(); + for (int i = num-1; i >= 0; --i) { // Flag removal @@ -407,7 +406,7 @@ void ctbkgcube::save(void) } // Get output filenames - std::string outfile = (*this)["outfile"].filename(); + std::string outfile = (*this)["outcube"].filename(); std::string outmodel = (*this)["outmodel"].filename(); // Create empty FITS file @@ -460,6 +459,8 @@ void ctbkgcube::save(void) void ctbkgcube::init_members(void) { // Initialise members + m_outcube.clear(); + m_outmodel.clear(); m_obs.clear(); m_bkgcube.clear(); m_bkgmdl.clear(); @@ -481,6 +482,8 @@ void ctbkgcube::init_members(void) void ctbkgcube::copy_members(const ctbkgcube& app) { // Copy members + m_outmodel = app.m_outmodel; + m_outcube = app.m_outcube; m_obs = app.m_obs; m_bkgcube = app.m_bkgcube; m_bkgmdl = app.m_bkgmdl; @@ -515,11 +518,18 @@ void ctbkgcube::free_members(void) ***************************************************************************/ void ctbkgcube::get_parameters(void) { - // If we do not have any observations in the container then load them - // using the "infile" parameter + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - get_obs(); - } + + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); + + // Build observation container + m_obs = get_observations(); + + } // endif: there was no observation in the container + // If there are no models associated with the observations then load now // the model definition from the XML file @@ -536,176 +546,47 @@ void ctbkgcube::get_parameters(void) m_obs.models(models); } - // If no counts map is specified then setup the background cube from - // the user parameters - std::string cntmap = (*this)["cntmap"].filename(); - if ((gammalib::toupper(cntmap) == "NONE") || - (gammalib::strip_whitespace(cntmap) == "")) { - - // Get user parameters for counts map definition - std::string proj = (*this)["proj"].string(); - std::string coordsys = (*this)["coordsys"].string(); - double xref = (*this)["xref"].real(); - double yref = (*this)["yref"].real(); - double binsz = (*this)["binsz"].real(); - int nxpix = (*this)["nxpix"].integer(); - int nypix = (*this)["nypix"].integer(); - - // Get energy definition - m_ebounds = get_ebounds(); - - // Define skymap for background cube - GSkymap map(proj, coordsys, xref, yref, - -binsz, binsz, nxpix, nypix, - m_ebounds.size()); - - // Allocate background cube using a dummy GTI - GGti gti; - gti.append(GTime(0.0), GTime(1.0)); // Dummy GTI - m_bkgcube = GCTAEventCube(map, m_ebounds, gti); + // Get the incube filename + std::string incube = (*this)["incube"].filename(); - } + // Check for filename validity + if ((gammalib::toupper(incube) == "NONE") || + (gammalib::strip_whitespace(incube) == "")) { + + // Create an event cube based on task parameters + m_bkgcube = create_cube(m_obs); + + } // endif: filename was not valid // ... otherwise setup the background cube from the counts map else { - // Set background cube from counts map - set_from_cntmap(cntmap); - - } - - // Read output filenames (if needed) - if (read_ahead()) { - (*this)["outfile"].filename(); - (*this)["outmodel"].filename(); - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Get observation definition - * - * @exception GException::invalid_value - * No event list, counts cube or observation definition file - * specified. - * - * Get observation definition from the user parameters. - ***************************************************************************/ -void ctbkgcube::get_obs(void) -{ - // Get input filename - std::string filename = (*this)["infile"].filename(); - - // Check whether infile parameter is undefined - if (gammalib::toupper(filename) == "NONE") { - std::string msg = "No event list, counts cube or observation definition " - "file specified. Please set the \"infile\" parameter " - "to either an event list filename, a counts cube " - "filename or an observation definition filename."; - throw GException::invalid_value(G_GET_OBS, msg); - } - - // Try first to open as FITS file - try { - - // Allocate CTA observation - GCTAObservation obs; - - // Load input file in CTA observation - obs.load(filename); - - // Append CTA observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(filename); - - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Set background cube definition from counts map - * - * @exception GException::invalid_value - * Invalid counts map projection or invalid events encountered. - * - * Set background cube definition from counts map. - ***************************************************************************/ -void ctbkgcube::set_from_cntmap(const std::string& filename) -{ - // Allocate CTA observation - GCTAObservation obs; - - // Load counts map in CTA observation - obs.load(filename); - - // Set background cube from counts map - const GCTAEventCube* cube = dynamic_cast(obs.events()); + // Load event cube from filename + m_bkgcube.load(incube); + // Loop over all bins in background cube + for (int i = 0; i < m_bkgcube.size(); ++i) { - // Continue only if cube is valid - if (cube != NULL) { + // Get event bin + GCTAEventBin* bin = m_bkgcube[i]; + bin->counts(0.0); - // Get sky map projection - const GWcs* wcs = dynamic_cast(cube->map().projection()); - - // Continue only if projection is valid - if (wcs != NULL) { - - // Get user parameters for counts map definition - std::string proj = wcs->code(); - std::string coordsys = wcs->coordsys(); - double xref = wcs->crval(0); - double yref = wcs->crval(1); - double dx = wcs->cdelt(0); - double dy = wcs->cdelt(1); - int nx = cube->map().nx(); - int ny = cube->map().ny(); - - // Get energy definition - m_ebounds = cube->ebounds(); - - // Define skymap for background cube - GSkymap map(proj, coordsys, xref, yref, - dx, dy, nx, ny, - m_ebounds.size()); - - // Allocate background cube - m_bkgcube = GCTAEventCube(map, m_ebounds, obs.events()->gti()); - - } // endif: WCS projection was valid - - // ... projection is not of WCS type - else { - std::string msg = "Counts map project is not of WCS type."; - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); } + + } // endelse: cube was loaded from file - } // endif: observation contained an events cube + // Get energy definition + m_ebounds = m_bkgcube.ebounds(); - // ... there is not events cube - else { - std::string msg = "No events cube found in file \"" - ""+filename+"\"."; - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); + // Read output filenames (if needed) + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); + m_outmodel = (*this)["outmodel"].filename(); } // Return return; } - /***********************************************************************//** * @brief Generate background cube * @@ -717,6 +598,7 @@ void ctbkgcube::set_from_cntmap(const std::string& filename) ***************************************************************************/ void ctbkgcube::fill_cube(GCTAObservation* obs) { + // Continue only if observation pointer is valid if (obs != NULL) { diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp index 1607e8a..aa9d065 100644 --- a/src/ctbkgcube/ctbkgcube.hpp +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbkgcube - Background cube generation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Chia-Chun Lu * + * copyright (C) 2014-2015 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTBKGCUBE_NAME "ctbkgcube" -#define CTBKGCUBE_VERSION "00-01-00" +#define CTBKGCUBE_VERSION "00-02-00" /***********************************************************************//** @@ -68,11 +68,11 @@ protected: void copy_members(const ctbkgcube& app); void free_members(void); void get_parameters(void); - void get_obs(void); - void set_from_cntmap(const std::string& filename); void fill_cube(GCTAObservation* obs); // Protected members + std::string m_outcube; //!< Filename of output cube + std::string m_outmodel; //!< Filename of output XML model GObservations m_obs; //!< Observation container GCTAEventCube m_bkgcube; //!< Background cube GModels m_bkgmdl; //!< CTA background models diff --git a/src/ctbkgcube/ctbkgcube.par b/src/ctbkgcube/ctbkgcube.par index 918e9ed..acfcee3 100644 --- a/src/ctbkgcube/ctbkgcube.par +++ b/src/ctbkgcube/ctbkgcube.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctbkgcube.par # -# Version: 00-01-00 # +# Version: 00-02-00 # # # # Author: Chia-Chun Lu # # c.lu@mpi-hd.mpg.de # @@ -15,13 +15,15 @@ # # ############################################################################# # -# General parameters +# Observational parameters #=================== -infile, f, a, "NONE",,, "Event list or observation definition file" +inobs, f, a, "NONE",,, "Input event list, cube or observation definition file" inmodel, f, a, "NONE",,, "Input (background) model XML file" -cntmap, f, a, "NONE",,, "Counts cube for background cube definition" -outfile, f, a, "bkgcube.fits",,, "Output background cube file" +incube, f, a, "NONE",,, "Counts cube for background cube definition" +outcube, f, a, "bkgcube.fits",,, "Output background cube file" outmodel, f, a, "NONE",,,"Output (background) model XML file" +caldb, s, a, "dummy",,, "Calibration database" +irf, s, a, "cta_dummy_irf",,, "Instrument response function" # # Energy binning parameters #========================== @@ -33,6 +35,7 @@ ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" # # Spatial binning parameters #=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" diff --git a/src/ctbutterfly/ctbutterfly.cpp b/src/ctbutterfly/ctbutterfly.cpp index 92c2c43..396aaa3 100644 --- a/src/ctbutterfly/ctbutterfly.cpp +++ b/src/ctbutterfly/ctbutterfly.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbutterfly - butterfly calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Michael Mayer * + * copyright (C) 2014-2015 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -399,7 +399,6 @@ void ctbutterfly::save(void) void ctbutterfly::init_members(void) { // Initialise members - m_infile.clear(); m_srcname.clear(); m_outfile.clear(); m_ebounds.clear(); @@ -424,7 +423,6 @@ void ctbutterfly::init_members(void) void ctbutterfly::copy_members(const ctbutterfly& app) { // Copy attributes - m_infile = app.m_infile; m_srcname = app.m_srcname; m_outfile = app.m_outfile; m_ebounds = app.m_ebounds; @@ -466,51 +464,25 @@ void ctbutterfly::free_members(void) ***************************************************************************/ void ctbutterfly::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Allocate CTA observation - GCTAObservation obs; + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Get event file name - std::string filename = (*this)["infile"].filename(); - - // Try first to open as FITS file - try { - - // Load data - obs.load(filename); - - // Set response - set_obs_response(&obs); - - // Append observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(filename); - - // Check if all observations have response information. If - // not, get the calibration database parameters and set - // the response properly - set_response(m_obs); - - } // endcatch: file was an XML file + // Build observation container + m_obs = get_observations(); } // endif: there was no observation in the container + // If there is are no models associated with the observations then // load now the model definition if (m_obs.models().size() == 0) { // Get models XML filename - std::string filename = (*this)["srcmdl"].filename(); + std::string filename = (*this)["inmodel"].filename(); // Setup models for optimizing. m_obs.models(GModels(filename)); @@ -526,15 +498,15 @@ void ctbutterfly::get_parameters(void) throw GException::invalid_value(G_GET_PARAMETERS, msg); } - // Get energy binning information - m_ebounds = get_ebounds(); + // Create energy boundaries from user parameters + m_ebounds = create_ebounds(); // Get matrix file name and load if possible std::string matrixfilename = (*this)["matrix"].filename(); if (matrixfilename != "NONE") { std::string msg = "Loading of matrix from file not implemented yet. " - "Use filename = \"NONE\" to induce a recomputation of " - "the matrix internally."; + "Use filename = \"NONE\" to induce a recomputation of " + "the matrix internally."; throw GException::feature_not_implemented(G_GET_PARAMETERS, msg); // m_covariance.load(matrixfilename); } diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp index 329175f..263c6d2 100644 --- a/src/ctbutterfly/ctbutterfly.hpp +++ b/src/ctbutterfly/ctbutterfly.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbutterfly - butterfly calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Michael Mayer * + * copyright (C) 2014-2015 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTBUTTERFLY_NAME "ctbutterfly" -#define CTBUTTERFLY_VERSION "00-01-00" +#define CTBUTTERFLY_VERSION "00-02-00" /***********************************************************************//** @@ -80,12 +80,11 @@ protected: void get_parameters(void); // User parameters - std::string m_infile; //!< Input file std::string m_srcname; //!< Name of source to compute butterfly std::string m_outfile; //!< Output ascii file GEbounds m_ebounds; //!< Energy binning definition - // Protected members + // Protected members GObservations m_obs; //!< Observation container GMatrixSparse m_covariance; //!< Covariance matrix std::vector m_energies; //!< Energy values for storage diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index 7cc0c48..751a077 100644 --- a/src/ctbutterfly/ctbutterfly.par +++ b/src/ctbutterfly/ctbutterfly.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctbutterfly # -# Version: 00-01-00 # +# Version: 00-02-00 # # # # Author: Michael Mayer # # michael.mayer@physik.hu-berlin.de # @@ -17,8 +17,8 @@ # # General parameters #=================== -infile, f, a, "events.fits",,,"Input event list, cube or observation definition file" -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" +inobs, f, a, "events.fits",,,"Input event list, cube or observation definition file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" srcname, s, a, "Crab",,,"Source of interest" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" diff --git a/src/ctcubemask/ctcubemask.cpp b/src/ctcubemask/ctcubemask.cpp index fd57519..3bf4732 100644 --- a/src/ctcubemask/ctcubemask.cpp +++ b/src/ctcubemask/ctcubemask.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctcubemask - Cube filter tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Chia-Chun Lu * + * copyright (C) 2014-2015 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -36,6 +36,7 @@ /* __ Method name definitions ____________________________________________ */ #define G_RUN "ctcubemask::run()" #define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*)" +#define G_GET_PARAMETERS "ctcubemask::get_parameters()" /* __ Debug definitions __________________________________________________ */ @@ -336,6 +337,9 @@ void ctcubemask::save(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctcubemask. + * * Get all task parameters from parameter file or (if required) by querying * the user. * @@ -345,46 +349,29 @@ void ctcubemask::save(void) ***************************************************************************/ void ctcubemask::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Get CTA event list file name - m_infile = (*this)["infile"].filename(); - - // Allocate CTA observation - GCTAObservation obs; - - // Try first to open as FITS file - try { - - // Load event cube in CTA observation - obs.load(m_infile); - - // Append CTA observation to container - m_obs.append(obs); - - // Signal that no XML file should be used for storage - m_use_xml = false; - + // Throw exception if no infile is given + if ((*this)["inobs"].filename() == "NONE" || + (*this)["inobs"].filename() == "") { + std::string msg = "A valid file needs to be specified for the " + "\"inobs\" parameter, yet \""+ + (*this)["inobs"].filename()+"\" was given." + " Specify a vaild observation definition or " + "event list FITS file to proceed."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(m_infile); - // Signal that XML file should be used for storage - m_use_xml = true; - - } + // Build observation container without response (not needed) + m_obs = get_observations(false); } // endif: there was no observation in the container // Get parameters m_regfile = (*this)["regfile"].filename(); - m_usepnt = (*this)["usepnt"].boolean(); + m_usepnt = (*this)["usepnt"].boolean(); if (!m_usepnt) { m_ra = (*this)["ra"].real(); m_dec = (*this)["dec"].real(); @@ -396,7 +383,7 @@ void ctcubemask::get_parameters(void) // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); } // Return @@ -538,9 +525,8 @@ void ctcubemask::apply_mask(GCTAObservation* obs) void ctcubemask::init_members(void) { // Initialise parameters - m_infile.clear(); m_regfile.clear(); - m_outfile.clear(); + m_outcube.clear(); m_prefix.clear(); m_usepnt = false; m_ra = 0.0; @@ -552,7 +538,6 @@ void ctcubemask::init_members(void) // Initialise protected members m_obs.clear(); m_infiles.clear(); - m_use_xml = false; // Return return; @@ -567,9 +552,8 @@ void ctcubemask::init_members(void) void ctcubemask::copy_members(const ctcubemask& app) { // Copy parameters - m_infile = app.m_infile; m_regfile = app.m_regfile; - m_outfile = app.m_outfile; + m_outcube = app.m_outcube; m_prefix = app.m_prefix; m_usepnt = app.m_usepnt; m_ra = app.m_ra; @@ -581,7 +565,6 @@ void ctcubemask::copy_members(const ctcubemask& app) // Copy protected members m_obs = app.m_obs; m_infiles = app.m_infiles; - m_use_xml = app.m_use_xml; // Return return; @@ -622,7 +605,7 @@ std::string ctcubemask::check_infile(const std::string& filename) const } catch (GException::fits_hdu_not_found& e) { message = "No \"IMAGE\" extension found in input file \"" - + m_outfile + "\"."; + + m_outcube + "\"."; } // Return @@ -661,13 +644,13 @@ std::string ctcubemask::set_outfile_name(const std::string& filename) const void ctcubemask::save_fits(void) { // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); // Get CTA observation from observation container GCTAObservation* obs = dynamic_cast(m_obs[0]); // Save event list - save_counts_map(obs, m_outfile); + save_counts_map(obs, m_outcube); // Return return; @@ -688,14 +671,14 @@ void ctcubemask::save_fits(void) void ctcubemask::save_xml(void) { // Get output filename and prefix - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["m_outcube"].filename(); m_prefix = (*this)["prefix"].string(); // Issue warning if output filename has no .xml suffix - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); + std::string suffix = gammalib::tolower(m_outcube.substr(m_outcube.length()-4,4)); if (suffix != ".xml") { log << "*** WARNING: Name of observation definition output file \""+ - m_outfile+"\"" << std::endl; + m_outcube+"\"" << std::endl; log << "*** WARNING: does not terminate with \".xml\"." << std::endl; log << "*** WARNING: This is not an error, but might be misleading." " It is recommended" << std::endl; @@ -726,7 +709,7 @@ void ctcubemask::save_xml(void) } // endfor: looped over observations // Save observations in XML file - m_obs.save(m_outfile); + m_obs.save(m_outcube); // Return return; diff --git a/src/ctcubemask/ctcubemask.hpp b/src/ctcubemask/ctcubemask.hpp index 5b2f961..f7a5756 100644 --- a/src/ctcubemask/ctcubemask.hpp +++ b/src/ctcubemask/ctcubemask.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctcubemask - Cube filter tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Chia-Chun Lu * + * copyright (C) 2014-2015 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -36,7 +36,7 @@ /* __Definitions _________________________________________________________ */ #define CTCUBEMASK_NAME "ctcubemask" -#define CTCUBEMASK_VERSION "00-01-00" +#define CTCUBEMASK_VERSION "00-02-00" /***********************************************************************//** @@ -78,9 +78,9 @@ protected: const std::string& outfile) const; // User parameters - std::string m_infile; //!< Input event list or XML file + //std::string m_infile; //!< Input event list or XML file std::string m_regfile; //!< ds9 region file - std::string m_outfile; //!< Output event list or XML file + std::string m_outcube; //!< Output event list or XML file std::string m_prefix; //!< Prefix for multiple counts maps bool m_usepnt; //!< Use pointing instead of RA/DEC parameters double m_ra; //!< RA of ROI centre @@ -92,7 +92,6 @@ protected: // Protected members GObservations m_obs; //!< Observations container std::vector m_infiles; //!< Input event filenames - bool m_use_xml; //!< Use XML file instead of FITS file }; diff --git a/src/ctcubemask/ctcubemask.par b/src/ctcubemask/ctcubemask.par index 66e6817..d8021a3 100644 --- a/src/ctcubemask/ctcubemask.par +++ b/src/ctcubemask/ctcubemask.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctcubemask.par # -# Version: 00-01-00 # +# Version: 00-02-00 # # # # Author: Chia-Chun Lu # # c.lu@mpi-hd.mpg.de # @@ -17,9 +17,9 @@ # # File information #================= -infile, f, a, "cntmap.fits",,,"Input counts cube or observation definition file" +inobs, f, a, "cntmap.fits",,,"Input counts cube or observation definition file" regfile, f, a, "NONE",,,"Exclusion region file in ds9 format" -outfile, f, a, "filtered_cube.fits",,,"Output event cube or observation definition file" +outcube, f, a, "filtered_cube.fits",,,"Output event cube or observation definition file" prefix, s, h, "filtered_",,,"Prefix for counts cube in observation definition file" # # Selection parameters diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp index bf69e38..fa4c134 100644 --- a/src/ctexpcube/ctexpcube.cpp +++ b/src/ctexpcube/ctexpcube.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctexpcube - Exposure cube generation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Juergen Knoedlseder * + * copyright (C) 2014-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -284,10 +284,10 @@ void ctexpcube::save(void) } // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); // Save exposure cube - m_expcube.save(m_outfile, clobber()); + m_expcube.save(m_outcube, clobber()); // Return return; @@ -306,13 +306,12 @@ void ctexpcube::save(void) void ctexpcube::init_members(void) { // Initialise members - m_outfile.clear(); + m_outcube.clear(); m_apply_edisp = false; // Initialise protected members m_obs.clear(); m_expcube.clear(); - m_ebounds.clear(); // Return return; @@ -327,13 +326,12 @@ void ctexpcube::init_members(void) void ctexpcube::copy_members(const ctexpcube& app) { // Copy attributes - m_outfile = app.m_outfile; + m_outcube = app.m_outcube; m_apply_edisp = app.m_apply_edisp; // Copy protected members - m_obs = app.m_obs; - m_expcube = app.m_expcube; - m_ebounds = app.m_ebounds; + m_obs = app.m_obs; + m_expcube = app.m_expcube; // Return return; @@ -355,164 +353,52 @@ void ctexpcube::free_members(void) * * Get all task parameters from parameter file or (if required) by querying * the user. The parameters are read in the correct order. - * - * @todo Setup exposure cube from counts map ***************************************************************************/ void ctexpcube::get_parameters(void) { - // If we do not have any observations in the container then get an - // input file name or observation descriptor file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - get_obs(); - } - - // For all observations that have no response, set the response - // from the task parameters - set_response(m_obs); - - // Read energy dispersion flag - m_apply_edisp = (*this)["edisp"].boolean(); - // If no counts map is specified then setup the exposure cube from - // the user parameters - std::string cntmap = (*this)["cntmap"].filename(); - if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) { - - // Get user parameters for counts map definition - std::string wcs = (*this)["proj"].string(); - std::string coordsys = (*this)["coordsys"].string(); - double xref = (*this)["xref"].real(); - double yref = (*this)["yref"].real(); - double binsz = (*this)["binsz"].real(); - int nxpix = (*this)["nxpix"].integer(); - int nypix = (*this)["nypix"].integer(); - - // Get energy definition - m_ebounds = get_ebounds(); + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Define exposure cube - m_expcube = GCTAExposure(wcs, coordsys, xref, yref, - -binsz, binsz, nxpix, nypix, - m_ebounds); - } + // Build observation container + m_obs = get_observations(); - // ... otherwise setup the exposure cube from the counts map - else { - - // Set exposure cube from counts map - set_from_cntmap(cntmap); - - } + } // endif: there was no observation in the container - // Read output filename (if needed) - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Get observation definition - * - * Get observation definition from the user parameters. - ***************************************************************************/ -void ctexpcube::get_obs(void) -{ - // Get input filename - std::string filename = (*this)["infile"].filename(); - - // Try first to open as FITS file - try { + // Read energy dispersion flag + m_apply_edisp = (*this)["edisp"].boolean(); - // Allocate CTA observation - GCTAObservation obs; - - // Load input file in CTA observation - obs.load(filename); + // Get the incube filename + std::string incube = (*this)["incube"].filename(); - // Append CTA observation to container - m_obs.append(obs); + // Check for filename validity + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { + // Create an event cube based on task parameters + GCTAEventCube cube = create_cube(m_obs); - // Load observations from XML file - m_obs.load(filename); + // Define exposure cube + m_expcube = GCTAExposure(cube); - } + } // endif: filename was not valid - // Return - return; -} + // ... otherwise setup the exposure cube from the counts map + else { + // Load event cube from filename + GCTAEventCube cube(incube); -/***********************************************************************//** - * @brief Set exposure cube definition from counts map - * - * @exception GException::invalid_value - * Invalid counts map projection or invalid events encountered. - * - * Set exposure cube definition from counts map. - ***************************************************************************/ -void ctexpcube::set_from_cntmap(const std::string& filename) -{ - // Allocate CTA observation - GCTAObservation obs; - - // Load counts map in CTA observation - obs.load(filename); - - // Set exposure cube from counts map - const GCTAEventCube* cube = dynamic_cast(obs.events()); - - // Continue only if cube is valid - if (cube != NULL) { - - // Get sky map projection - const GWcs* wcs = dynamic_cast(cube->map().projection()); - - // Continue only if projection is valid - if (wcs != NULL) { - - // Get user parameters for counts map definition - std::string proj = wcs->code(); - std::string coordsys = wcs->coordsys(); - double xref = wcs->crval(0); - double yref = wcs->crval(1); - double dx = wcs->cdelt(0); - double dy = wcs->cdelt(1); - int nx = cube->map().nx(); - int ny = cube->map().ny(); - - // Get energy definition - m_ebounds = cube->ebounds(); - - // Define exposure cube - m_expcube = GCTAExposure(proj, coordsys, xref, yref, - dx, dy, nx, ny, - m_ebounds); - - } // endif: WCS projection was valid - - // ... projection is not of WCS type - else { - std::string msg = "Counts map project is not of WCS type."; - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); - } + // Define exposure cube + m_expcube = GCTAExposure(cube); - } // endif: observation contained an events cube + } // endelse: cube loaded from file - // ... there is not events cube - else { - std::string msg = "No events cube found in file \"" - ""+filename+"\"."; - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); + // Read output filename (if needed) + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); } // Return diff --git a/src/ctexpcube/ctexpcube.hpp b/src/ctexpcube/ctexpcube.hpp index 6a648d4..cd43db6 100644 --- a/src/ctexpcube/ctexpcube.hpp +++ b/src/ctexpcube/ctexpcube.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTEXPCUBE_NAME "ctexpcube" -#define CTEXPCUBE_VERSION "00-01-00" +#define CTEXPCUBE_VERSION "00-02-00" /***********************************************************************//** @@ -67,17 +67,14 @@ protected: void copy_members(const ctexpcube& app); void free_members(void); void get_parameters(void); - void get_obs(void); - void set_from_cntmap(const std::string& filename); // User parameters - std::string m_outfile; //!< Output exposure cube file + std::string m_outcube; //!< Output exposure cube file bool m_apply_edisp; //!< Apply energy dispersion? // Protected members GObservations m_obs; //!< Observation container GCTAExposure m_expcube; //!< Exposure cube - GEbounds m_ebounds; //!< Energy boundaries }; diff --git a/src/ctexpcube/ctexpcube.par b/src/ctexpcube/ctexpcube.par index c1f02fa..bf9da56 100644 --- a/src/ctexpcube/ctexpcube.par +++ b/src/ctexpcube/ctexpcube.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctexpcube.par # -# Version: 00-01-00 # +# Version: 00-02-00 # # # # Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -15,11 +15,11 @@ # # ############################################################################# # -# General parameters +# Observational parameters #=================== -infile, f, a, "NONE",,, "Event list or observation definition file" -cntmap, f, a, "NONE",,, "Counts map for exposure cube definition" -outfile, f, a, "expcube.fits",,, "Output exposure cube file" +inobs, f, a, "NONE",,, "Event list or observation definition file" +incube, f, a, "NONE",,, "Counts map for exposure cube definition" +outcube, f, a, "expcube.fits",,, "Output exposure cube file" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" @@ -34,6 +34,7 @@ ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" # # Spatial binning parameters #=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index c0092d0..5766de3 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctlike - Maximum likelihood fitting tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -326,11 +326,11 @@ void ctlike::save(void) } // Get output filename - m_outmdl = (*this)["outmdl"].filename(); + m_outmodel = (*this)["outmodel"].filename(); // Write results out as XML model - if (gammalib::toupper(m_outmdl) != "NONE") { - m_obs.models().save(m_outmdl); + if (gammalib::toupper(m_outmodel) != "NONE") { + m_obs.models().save(m_outmodel); } // Return @@ -357,56 +357,34 @@ void ctlike::save(void) ***************************************************************************/ void ctlike::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Allocate CTA observation - GCTAObservation obs; + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Get event file name - std::string filename = (*this)["infile"].filename(); + // Get observation container + m_obs = get_observations(); - // Try first to open as FITS file - try { - - // Load data - obs.load(filename); - - // Get other task parameters - std::string statistics = gammalib::toupper((*this)["stat"].string()); - - // Set statistics - obs.statistics(statistics); - - // Set response - set_obs_response(&obs); - - // Append observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(filename); + } // endif: there was no observation in the container - // For all observations that have no response, set the response - // from the task parameters - set_response(m_obs); + // If only single observation is used, read statistics parameter + if (!m_use_xml) { - } // endcatch: file was an XML file + // Get other task parameters + std::string statistics = gammalib::toupper((*this)["stat"].string()); - } // endif: there was no observation in the container + // Set statistics + (*m_obs[0]).statistics(statistics); + } // If there is are no models associated with the observations then // load now the model definition if (m_obs.models().size() == 0) { // Get models XML filename - std::string filename = (*this)["srcmdl"].filename(); + std::string filename = (*this)["inmodel"].filename(); // Setup models for optimizing. m_obs.models(GModels(filename)); @@ -420,7 +398,7 @@ void ctlike::get_parameters(void) // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { - m_outmdl = (*this)["outmdl"].filename(); + m_outmodel = (*this)["outmodel"].filename(); } // Return @@ -547,7 +525,7 @@ double ctlike::reoptimize_lm(void) void ctlike::init_members(void) { // Initialise members - m_outmdl.clear(); + m_outmodel.clear(); m_obs.clear(); m_refit = false; m_max_iter = 100; // Set maximum number of iterations @@ -573,7 +551,7 @@ void ctlike::copy_members(const ctlike& app) { // Copy attributes m_refit = app.m_refit; - m_outmdl = app.m_outmdl; + m_outmodel = app.m_outmodel; m_obs = app.m_obs; m_max_iter = app.m_max_iter; m_max_stall = app.m_max_stall; diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp index c71d273..6e7ac22 100644 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctlike - Maximum likelihood fitting tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -35,7 +35,7 @@ /* __Definitions _________________________________________________________ */ #define CTLIKE_NAME "ctlike" -#define CTLIKE_VERSION "00-06-00" +#define CTLIKE_VERSION "00-07-00" /***********************************************************************//** @@ -74,7 +74,7 @@ protected: // User parameters bool m_refit; //!< Refitting - std::string m_outmdl; //!< Source model output XML file + std::string m_outmodel; //!< Source model output XML file bool m_apply_edisp; //!< Apply energy dispersion? // Members diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index fc63e6c..bb766d2 100644 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctlike.par # -# Version: 00-05-00 # +# Version: 00-06-00 # # # # Author: Jurgen Knodlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -17,14 +17,14 @@ # # General parameters #=================== -infile, f, a, "events.fits",,, "Event list, counts cube or observation definition file" -stat, s, h, "POISSON",,, "Optimization statistics" -edisp, b, h, no,,, "Apply energy dispersion?" -refit, b, h, no,,, "Do refitting?" -caldb, s, a, "dummy",,, "Calibration database" -irf, s, a, "cta_dummy_irf",,, "Instrument response function" -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" -outmdl, f, a, "crab_results.xml",,, "Source model output file" +inobs, f, a, "events.fits",,, "Event list, counts cube or observation definition file" +stat, s, h, "POISSON",,, "Optimization statistics" +edisp, b, h, no,,, "Apply energy dispersion?" +refit, b, h, no,,, "Do refitting?" +caldb, s, a, "dummy",,, "Calibration database" +irf, s, a, "cta_dummy_irf",,, "Instrument response function" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" +outmodel, f, a, "crab_results.xml",,, "Source model output file" # # Standard parameters #==================== diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 1d77402..2532cba 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -33,6 +33,7 @@ #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_GET_PARAMETERS "ctmodel::get_parameters()" /* __ Debug definitions __________________________________________________ */ @@ -207,9 +208,6 @@ void ctmodel::run(void) // Get task parameters get_parameters(); - // Setup observation container - setup_obs(); - // Write parameters into logger if (logTerse()) { log_parameters(); @@ -246,11 +244,6 @@ void ctmodel::run(void) log.header1("Generate model cube"); } - // Initialise counts cube if cube has not yet been set - if (!has_cube()) { - init_cube(); - } - // Loop over all observations in the container for (int i = 0; i < m_obs.size(); ++i) { @@ -318,10 +311,10 @@ void ctmodel::save(void) } // Make sure we have the FITS filename - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); // Save model cube into FITS file - m_cube.save(m_outfile, clobber()); + m_cube.save(m_outcube, clobber()); // Return return; @@ -365,35 +358,12 @@ void ctmodel::cube(const GCTAEventCube& cube) void ctmodel::init_members(void) { // Initialise members - m_infile.clear(); - m_obsfile.clear(); - m_outfile.clear(); - m_caldb.clear(); - m_irf.clear(); - m_srcmdl.clear(); - m_proj.clear(); - m_coordsys.clear(); - m_ebinalg.clear(); - m_ebinfile.clear(); - m_ra = 0.0; - m_dec = 0.0; - m_deadc = 1.0; - m_tmin = 0.0; - m_tmax = 0.0; - m_emin = 0.0; - m_emax = 0.0; - m_enumbins = 0; - m_xref = 0.0; - m_yref = 0.0; - m_binsz = 0.0; - m_nxpix = 0; - m_nypix = 0; + m_outcube.clear(); m_apply_edisp = false; // Initialise protected members m_obs.clear(); m_cube.clear(); - m_ebounds.clear(); m_gti.clear(); m_has_cube = false; @@ -410,35 +380,12 @@ void ctmodel::init_members(void) void ctmodel::copy_members(const ctmodel& app) { // Copy attributes - m_infile = app.m_infile; - m_obsfile = app.m_obsfile; - m_outfile = app.m_outfile; - m_caldb = app.m_caldb; - m_irf = app.m_irf; - m_srcmdl = app.m_srcmdl; - m_proj = app.m_proj; - m_coordsys = app.m_coordsys; - m_ebinalg = app.m_ebinalg; - m_ebinfile = app.m_ebinfile; - m_ra = app.m_ra; - m_dec = app.m_dec; - m_deadc = app.m_deadc; - m_tmin = app.m_tmin; - m_tmax = app.m_tmax; - m_emin = app.m_emin; - m_emax = app.m_emax; - m_enumbins = app.m_enumbins; - m_xref = app.m_xref; - m_yref = app.m_yref; - m_binsz = app.m_binsz; - m_nxpix = app.m_nxpix; - m_nypix = app.m_nypix; + m_outcube = app.m_outcube; m_apply_edisp = app.m_apply_edisp; // Copy protected members m_obs = app.m_obs; m_cube = app.m_cube; - m_ebounds = app.m_ebounds; m_gti = app.m_gti; m_has_cube = app.m_has_cube; @@ -465,255 +412,65 @@ void ctmodel::free_members(void) ***************************************************************************/ void ctmodel::get_parameters(void) { - // Read input and optionally output cube filenames - if (!has_cube()) { - m_infile = (*this)["infile"].filename(); - } - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); + // If there are no observations in container then load them via user + // parameters + if (m_obs.size() == 0) { + m_obs = get_observations(); } - // Read model definiton file if required + // Read model definition file if required if (m_obs.models().size() == 0) { - m_srcmdl = (*this)["srcmdl"].filename(); + std::string inmodel = (*this)["inmodel"].filename(); + + // Load models from file + m_obs.models(inmodel); + + } // endif: there were no models + + // Read optionally output cube filenames + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); } // Get energy dispersion flag parameters m_apply_edisp = (*this)["edisp"].boolean(); - // If the observation container is empty then read the parameters that - // are necessary to fill it - if (m_obs.size() == 0) { - - // Read observation definition filename and response parameters - m_obsfile = (*this)["obsfile"].filename(); - - // If no observation definition file has been specified then read all - // parameters that are necessary to create an observation from scratch - // (see method setup_obs) - if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) { - m_ra = (*this)["ra"].real(); - m_dec = (*this)["dec"].real(); - m_deadc = (*this)["deadc"].real(); - m_tmin = (*this)["tmin"].real(); - m_tmax = (*this)["tmax"].real(); - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - } + if (!m_has_cube) { - } + // Read cube definition file + std::string incube = (*this)["incube"].filename(); - // Check if we need response information - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - if (!obs->has_response()) { - m_caldb = (*this)["caldb"].string(); - m_irf = (*this)["irf"].string(); - break; - } - } - } + // If no cube file has been specified then create a cube from + // the task parameters ... + if ((incube == "NONE") || + (gammalib::strip_whitespace(incube) == "")) { + + // Create cube from scratch + m_cube = create_cube(m_obs); - // If no cube file has been specified then read all parameters that - // are necessary to create the cube from scratch (see method init_cube) - if ((!has_cube()) && - ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == ""))) { - m_ebinalg = (*this)["ebinalg"].string(); - if (m_ebinalg == "FILE") { - m_ebinfile = (*this)["ebinfile"].filename(); } - else { - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - m_enumbins = (*this)["enumbins"].integer(); - } - m_proj = (*this)["proj"].string(); - m_coordsys = (*this)["coordsys"].string(); - m_xref = (*this)["xref"].real(); - m_yref = (*this)["yref"].real(); - m_binsz = (*this)["binsz"].real(); - m_nxpix = (*this)["nxpix"].integer(); - m_nypix = (*this)["nypix"].integer(); - } - - // Return - return; -} - -/***********************************************************************//** - * @brief Setup observation container - * - * This method makes sure that we have an observation container with an - * associated model and that all CTA observations have a valid response - * functions. - * - * If on input the observation container is empty, the method sets up an - * observation container by either loading the information from an - * observation definition XML file, by loading an event list or a counts - * cube as a single observation, or by constructing an observation from - * the scratch from the task parameters. - * - * If no model exists yet in the observation container, the method will load - * the model from the file given by the task parameter "srcmdl". - * - * The method then loops over all observations in the container and makes - * sure that all CTA observations have a valid response function. If no - * response function exists so far, the response function will be set - * using the "caldb" and "irf" task parameters. - ***************************************************************************/ -void ctmodel::setup_obs(void) -{ - // If there are no observations in the container then get them either - // from the observation definition file or from the task parameters - if (m_obs.size() == 0) { - - // If no observation definition file has been specified, then create - // an observation container from the task parameters - if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) { - - // Set pointing direction - GCTAPointing pnt; - GSkyDir skydir; - skydir.radec_deg(m_ra, m_dec); - pnt.dir(skydir); - - // Setup time interval covered by model - GGti gti; - GTime tmin(m_tmin); - GTime tmax(m_tmax); - gti.append(tmin, tmax); - - // Setup energy range covered by model - GEnergy emin(m_emin, "TeV"); - GEnergy emax(m_emax, "TeV"); - GEbounds ebounds(1, emin, emax); - - // Allocate CTA observation and empty event list - GCTAObservation obs; - GCTAEventList list; - - // Set event list GTI and energy boundaries and append event list - // to observation (we need in fact only the GTI for the - // computations) - list.gti(gti); - list.ebounds(ebounds); - obs.events(list); - - // Set CTA observation attributes - obs.pointing(pnt); - obs.ontime(gti.ontime()); - obs.livetime(gti.ontime()*m_deadc); - obs.deadc(m_deadc); - - // Append CTA observation to container - m_obs.append(obs); - - } // endif: created single CTA observation from task parameters - - // ... otherwise try to load information from the file + // ... otherwise load the cube from file and reset all bins + // to zero else { - // First try to open the file as an event list or counts map - try { - - // Allocate CTA observation - GCTAObservation obs; - - // Load CTA observation from FITS file - obs.load(m_obsfile); - - // Append CTA observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file. This will throw - // an exception if it fails. - m_obs.load(m_obsfile); + // Load cube from given file + m_cube.load(incube); + // Set all cube bins to zero + for (int i = 0; i < m_cube.size(); ++i) { + m_cube[i]->counts(0.0); } - - } // endelse: loaded information from input file - - } // endif: there was no observation in the container - - // If there are no models associated with the observations then - // load now the model definition from the XML file - if (m_obs.models().size() == 0) { - m_obs.models(GModels(m_srcmdl)); - } - - // For all observations that have no response, set the response - // from the task parameters - set_response(m_obs); - - // Return - return; -} - - -/***********************************************************************//** - * @brief Initialise model cube - * - * Initialises model cube by either reading the cube definition from an - * existing cube file specified by the "infile" parameter, or by constructing - * a cube from the energy and spatial binning information provided by the - * task parameters. The GTIs of the model cube will be undefined on - * exit. - ***************************************************************************/ -void ctmodel::init_cube(void) -{ - // Clear model cube and GTIs - m_cube.clear(); - m_gti.clear(); - - // If no cube file has been specified then create a cube from the task - // parameters - if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) { - - // Set dummy GTI that is needed for event cube creation - GGti gti; - gti.append(GTime(0.0), GTime(0.1234)); - - // Set energy boundaries - m_ebounds = get_ebounds(); - - // Setup skymap - GSkymap map = GSkymap(m_proj, m_coordsys, - m_xref, m_yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, m_enumbins); - - // Create model cube - m_cube = GCTAEventCube(map, m_ebounds, gti); - - } - - // ... otherwise load cube from file - else { - - // Load cube - m_cube.load(m_infile); - - // Set all cube bins to zero - for (int i = 0; i < m_cube.size(); ++i) { - m_cube[i]->counts(0.0); } + // Signal that cube has been set + m_has_cube = true; } - // Signal that cube has been set - m_has_cube = true; - // Return return; } - /***********************************************************************//** * @brief Fill model into model cube * diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index c4520a3..90b8cc4 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTMODEL_NAME "ctmodel" -#define CTMODEL_VERSION "00-03-00" +#define CTMODEL_VERSION "00-04-00" /***********************************************************************//** @@ -78,41 +78,16 @@ protected: void copy_members(const ctmodel& app); void free_members(void); void get_parameters(void); - void setup_obs(void); - void init_cube(void); void fill_cube(const GCTAObservation* obs); bool has_cube(void) const; // User parameters - std::string m_infile; //!< Input counts cube - std::string m_obsfile; //!< Event list, counts map or observation definition XML file - std::string m_outfile; //!< Output model cube - std::string m_caldb; //!< Calibration database - std::string m_irf; //!< Instrument response functions - std::string m_srcmdl; //!< Source model - double m_ra; //!< RA of pointing direction - double m_dec; //!< DEC of pointing direction - double m_deadc; //!< Deadtime correction factor - double m_tmin; //!< Start time - double m_tmax; //!< Stop time - std::string m_ebinalg; //!< Algorithm for energy binning - std::string m_ebinfile; //!< FITS-file containing energy binning - double m_emin; //!< Lower energy - double m_emax; //!< Upper energy - int m_enumbins; //!< Number of energy bins - std::string m_proj; //!< WCS projection - std::string m_coordsys; //!< Coordinate system - double m_xref; //!< Longitude reference coordinate - double m_yref; //!< Latitude reference coordinate - double m_binsz; //!< Pixel size - int m_nxpix; //!< Number of pixels in longitude - int m_nypix; //!< Number of pixels in latitude + std::string m_outcube; //!< Output model cube bool m_apply_edisp; //!< Apply energy dispersion? // Protected members GObservations m_obs; //!< Observation container GCTAEventCube m_cube; //!< Model cube - GEbounds m_ebounds; //!< Energy boundaries GGti m_gti; //!< Model cube GTIs bool m_has_cube; //!< Signal if cube has been set or loaded }; diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par index 77e510d..4593962 100644 --- a/src/ctmodel/ctmodel.par +++ b/src/ctmodel/ctmodel.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctmodel.par # -# Version: 00-03-00 # +# Version: 00-04-00 # # # # Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -17,18 +17,19 @@ # # General parameters #=================== -infile, f, a, "NONE",,, "Input counts map or observation definition file" -outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file" -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" +incube, f, a, "NONE",,, "Input counts map or observation definition file" +outcube, f, a, "modmap.fits",,, "Output counts map or observation definition file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" edisp, b, h, no,,, "Apply energy dispersion?" # # Observation parameters #======================= -obsfile, f, a, "NONE",,, "Observation definition file" +inobs, f, a, "NONE",,, "Input event list, cube or observation definition file" caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" ra, r, a, 83.63,0,360, "RA of pointing (degrees)" dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)" +rad, r, a, 5.0, 0,180, "Radius of FOV (degrees)" tmin, r, a, 0.0,0,, "Start time (MET in s)" tmax, r, a, 1800.0,0,, "End time (MET in s)" deadc, r, h, 0.95,0,1, "Average deadtime correction factor" @@ -43,6 +44,7 @@ ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" # # Spatial binning parameters #=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 4da2029..6a44f21 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctobssim - Observation simulator tool * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2014 by Juergen Knoedlseder * + * copyright (C) 2011-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -377,7 +377,7 @@ void ctobssim::run(void) outfile = m_prefix + gammalib::str(i) + ".fits"; } else { - outfile = (*this)["outfile"].filename(); + outfile = (*this)["outevents"].filename(); } // Store output file name in original observation @@ -490,35 +490,20 @@ void ctobssim::save(void) void ctobssim::init_members(void) { // Initialise user parameters - m_infile.clear(); - m_outfile.clear(); + m_outevents.clear(); m_prefix.clear(); - m_seed = 1; - m_ra = 0.0; - m_dec = 0.0; - m_rad = 0.0; - m_tmin = 0.0; - m_tmax = 0.0; - m_emin = 0.0; - m_emax = 0.0; + m_seed = 1; m_apply_edisp = false; - m_deadc = 1.0; // Initialise protected members m_rans.clear(); m_obs.clear(); - m_use_xml = false; m_save_and_dispose = false; // Set fixed parameters m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2 m_max_photons = 1000000; //!< Maximum number of photons / time slice - // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD, - // which is defined in GCTALib.hpp. This is somehow a kluge. We need - // a better mechanism to implement the CTA reference MJD. - m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL"); - // Initialise first event identifier m_event_id = 1; @@ -535,28 +520,17 @@ void ctobssim::init_members(void) void ctobssim::copy_members(const ctobssim& app) { // Copy user parameters - m_infile = app.m_infile; - m_outfile = app.m_outfile; + m_outevents = app.m_outevents; m_prefix = app.m_prefix; m_seed = app.m_seed; - m_ra = app.m_ra; - m_dec = app.m_dec; - m_rad = app.m_rad; - m_tmin = app.m_tmin; - m_tmax = app.m_tmax; - m_emin = app.m_emin; - m_emax = app.m_emax; m_apply_edisp = app.m_apply_edisp; - m_deadc = app.m_deadc; // Copy protected members m_area = app.m_area; m_max_photons = app.m_max_photons; m_rans = app.m_rans; m_obs = app.m_obs; - m_use_xml = app.m_use_xml; m_save_and_dispose = app.m_save_and_dispose; - m_cta_ref = app.m_cta_ref; m_event_id = app.m_event_id; // Return @@ -585,47 +559,16 @@ void ctobssim::free_members(void) ***************************************************************************/ void ctobssim::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { + m_obs = get_observations(); + } - // Get CTA observation parameters - m_infile = (*this)["infile"].filename(); - m_ra = (*this)["ra"].real(); - m_dec = (*this)["dec"].real(); - - // Set pointing direction - GCTAPointing pnt; - GSkyDir skydir; - skydir.radec_deg(m_ra, m_dec); - pnt.dir(skydir); - - // Allocate CTA observation - GCTAObservation obs; - - // Set response - set_obs_response(&obs); - - // Set pointing - obs.pointing(pnt); - - // Set event list (queries remaining parameters) - set_list(&obs); - - // Append CTA observation to container - m_obs.append(obs); - - // Load models into container - m_obs.models(m_infile); - - // Signal that no XML file should be used for storage - m_use_xml = false; - - } // endif: there was no observation in the container - - // ... otherwise we signal that an XML file should be written - else { - m_use_xml = true; + // Read model definition file if required + if (m_obs.models().size() == 0) { + std::string inmodel = (*this)["inmodel"].filename(); + m_obs.models(inmodel); } // Get other parameters @@ -635,18 +578,8 @@ void ctobssim::get_parameters(void) // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - m_prefix = (*this)["prefix"].string(); - } - - // Make sure that all observations hold a CTA event list - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - if (dynamic_cast(obs->events()) == NULL) { - set_list(obs); - } - } + m_outevents = (*this)["outevents"].filename(); + m_prefix = (*this)["prefix"].string(); } // Initialise random number generators. We initialise here one random @@ -700,81 +633,6 @@ void ctobssim::get_parameters(void) /***********************************************************************//** - * @brief Set empty CTA event list - * - * @param[in] obs CTA observation. - * - * Attaches an empty event list to CTA observation. The method also sets the - * pointing direction using the m_ra and m_dec members, the ROI based on - * m_ra, m_dec and m_rad, a single GTI based on m_tmin and m_tmax, and a - * single energy boundary based on m_emin and m_emax. The method furthermore - * sets the ontime, livetime and deadtime correction factor. - ***************************************************************************/ -void ctobssim::set_list(GCTAObservation* obs) -{ - // Continue only if observation is valid - if (obs != NULL) { - - // Get CTA observation parameters - m_ra = (*this)["ra"].real(); - m_dec = (*this)["dec"].real(); - m_rad = (*this)["rad"].real(); - m_tmin = (*this)["tmin"].real(); - m_tmax = (*this)["tmax"].real(); - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - m_deadc = (*this)["deadc"].real(); - - // Allocate CTA event list - GCTAEventList events; - - // Set pointing direction - GCTAPointing pnt; - GSkyDir skydir; - skydir.radec_deg(m_ra, m_dec); - pnt.dir(skydir); - - // Set ROI - GCTAInstDir instdir(skydir); - GCTARoi roi(instdir, m_rad); - - // Set GTI - GGti gti(m_cta_ref); - GTime tstart; - GTime tstop; - tstart.set(m_tmin, m_cta_ref); - tstop.set(m_tmax, m_cta_ref); - gti.append(tstart, tstop); - - // Set energy boundaries - GEbounds ebounds; - GEnergy emin; - GEnergy emax; - emin.TeV(m_emin); - emax.TeV(m_emax); - ebounds.append(emin, emax); - - // Set CTA event list attributes - events.roi(roi); - events.gti(gti); - events.ebounds(ebounds); - - // Attach event list to CTA observation - obs->events(events); - - // Set observation ontime, livetime and deadtime correction factor - obs->ontime(gti.ontime()); - obs->livetime(gti.ontime()*m_deadc); - obs->deadc(m_deadc); - - } // endif: oberservation was valid - - // Return - return; -} - - -/***********************************************************************//** * @brief Simulate source events from photon list * * @param[in] obs Pointer on CTA observation. @@ -1163,13 +1021,13 @@ void ctobssim::save_fits(void) if (!m_save_and_dispose) { // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outevents = (*this)["outevents"].filename(); // Get CTA observation from observation container GCTAObservation* obs = dynamic_cast(m_obs[0]); // Save observation into FITS file - obs->save(m_outfile, clobber()); + obs->save(m_outevents, clobber()); } // endif: event list has not yet been saved and disposed @@ -1191,14 +1049,14 @@ void ctobssim::save_fits(void) void ctobssim::save_xml(void) { // Get output filename and prefix - m_outfile = (*this)["outfile"].filename(); + m_outevents = (*this)["outevents"].filename(); m_prefix = (*this)["prefix"].string(); // Issue warning if output filename has no .xml suffix - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); + std::string suffix = gammalib::tolower(m_outevents.substr(m_outevents.length()-4,4)); if (suffix != ".xml") { log << "*** WARNING: Name of observation definition output file \""+ - m_outfile+"\"" << std::endl; + m_outevents+"\"" << std::endl; log << "*** WARNING: does not terminate with \".xml\"." << std::endl; log << "*** WARNING: This is not an error, but might be misleading." " It is recommended" << std::endl; @@ -1239,7 +1097,7 @@ void ctobssim::save_xml(void) } // endif: event list has not yet been saved and disposed // Save observations in XML file - m_obs.save(m_outfile); + m_obs.save(m_outevents); // Return return; diff --git a/src/ctobssim/ctobssim.hpp b/src/ctobssim/ctobssim.hpp index 6d0a901..e58276e 100644 --- a/src/ctobssim/ctobssim.hpp +++ b/src/ctobssim/ctobssim.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctobssim - Observation simulator tool * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2014 by Juergen Knoedlseder * + * copyright (C) 2011-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTOBSSIM_NAME "ctobssim" -#define CTOBSSIM_VERSION "00-09-00" +#define CTOBSSIM_VERSION "00-10-00" /***********************************************************************//** @@ -73,7 +73,6 @@ protected: void copy_members(const ctobssim& app); void free_members(void); void get_parameters(void); - void set_list(GCTAObservation* obs); void simulate_source(GCTAObservation* obs, const GModels& models, GRan& ran, @@ -86,19 +85,10 @@ protected: void save_xml(void); // User parameters - std::string m_infile; //!< Input model - std::string m_outfile; //!< Output events file + std::string m_outevents; //!< Output events file std::string m_prefix; //!< Prefix for multiple event lists int m_seed; //!< Random number generator seed - double m_ra; //!< RA of pointing direction - double m_dec; //!< DEC of pointing direction - double m_rad; //!< FOV radius (degrees) - double m_tmin; //!< Start time (MET) - double m_tmax; //!< Stop time (MET) - double m_emin; //!< Lower energy (TeV) - double m_emax; //!< Upper energy (TeV) bool m_apply_edisp; //!< Apply energy dispersion? - double m_deadc; //!< Average deadtime correction // Protected members mutable GObservations m_obs; //!< Observation container @@ -106,8 +96,6 @@ protected: double m_area; //!< Surface area for simulation (cm2) int m_max_photons; //!< Maximum number of photons/slice std::vector m_rans; //!< Random number generators - bool m_use_xml; //!< Use XML file instead of FITS file - GTimeReference m_cta_ref; //!< CTA time reference int m_event_id; //!< Event identifier }; @@ -120,7 +108,6 @@ protected: inline const GObservations& ctobssim::obs(void) const { - // Return observation container return m_obs; } diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par index 19e9e49..073e6a1 100644 --- a/src/ctobssim/ctobssim.par +++ b/src/ctobssim/ctobssim.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctobssim.par # -# Version: 00-08-00 # +# Version: 00-10-00 # # # # Author: Jurgen Knodlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -17,9 +17,10 @@ # # File information #================= -infile, f, a, "$CTOOLS/share/models/crab.xml",,, "Model" -outfile,f, a, "events.fits",,, "Output event data file or observation definition file" -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file" +inobs, f, h, "NONE",,,"Input event list or observation definition file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Model" +outevents, f, a, "events.fits",,,"Output event data file or observation definition file" +prefix, s, h, "sim_events_",,,"Prefix for event lists in observation definition file" # # Observation parameters #======================= diff --git a/src/ctpsfcube/ctpsfcube.cpp b/src/ctpsfcube/ctpsfcube.cpp index 64a6097..0e88fe4 100644 --- a/src/ctpsfcube/ctpsfcube.cpp +++ b/src/ctpsfcube/ctpsfcube.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctpsfcube - PSF cube generation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Chia-Chun Lu * + * copyright (C) 2014-2015 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -35,7 +35,6 @@ /* __ Method name definitions ____________________________________________ */ #define G_GET_PARAMETERS "ctpsfcube::get_parameters()" -#define G_SET_FROM_CNTMAP "ctpsfcube::set_from_cntmap(std::string&)" /* __ Debug definitions __________________________________________________ */ @@ -275,10 +274,10 @@ void ctpsfcube::save(void) } // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outcube = (*this)["outcube"].filename(); // Save PSF cube - m_psfcube.save(m_outfile, clobber()); + m_psfcube.save(m_outcube, clobber()); // Return return; @@ -297,13 +296,12 @@ void ctpsfcube::save(void) void ctpsfcube::init_members(void) { // Initialise members - m_outfile.clear(); + m_outcube.clear(); m_apply_edisp = false; // Initialise protected members m_obs.clear(); m_psfcube.clear(); - m_ebounds.clear(); // Return return; @@ -318,13 +316,12 @@ void ctpsfcube::init_members(void) void ctpsfcube::copy_members(const ctpsfcube& app) { // Copy attributes - m_outfile = app.m_outfile; + m_outcube = app.m_outcube; m_apply_edisp = app.m_apply_edisp; // Copy protected members m_obs = app.m_obs; m_psfcube = app.m_psfcube; - m_ebounds = app.m_ebounds; // Return return; @@ -346,167 +343,56 @@ void ctpsfcube::free_members(void) * * Get all task parameters from parameter file or (if required) by querying * the user. The parameters are read in the correct order. - * - * @todo Setup PSF cube from counts map ***************************************************************************/ void ctpsfcube::get_parameters(void) { - // If we do not have any observations in the container then get an - // input file name or observation descriptor file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - get_obs(); - } - - // Make sure that response is set - set_response(m_obs); - // Read energy dispersion flag - m_apply_edisp = (*this)["edisp"].boolean(); + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // If no counts map is specified then setup the PSF cube from - // the user parameters - std::string cntmap = (*this)["cntmap"].filename(); - if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) { - - // Get user parameters for counts map definition - std::string wcs = (*this)["proj"].string(); - std::string coordsys = (*this)["coordsys"].string(); - double xref = (*this)["xref"].real(); - double yref = (*this)["yref"].real(); - double binsz = (*this)["binsz"].real(); - int nxpix = (*this)["nxpix"].integer(); - int nypix = (*this)["nypix"].integer(); - double dmax = (*this)["amax"].real(); - int ndbins = (*this)["anumbins"].integer(); - - // Get energy definition - m_ebounds = get_ebounds(); - - // Define PSF cube - m_psfcube = GCTAMeanPsf(wcs, coordsys, xref, yref, - -binsz, binsz, nxpix, nypix, - m_ebounds, dmax, ndbins); - - } // endif: PSF cube set from user parameters - - // ... otherwise setup the PSF cube from the counts map - else { - - // Set PSF cube from counts map - set_from_cntmap(cntmap); - - } - - // Read output filename (if needed) - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - } + // Build observation container + m_obs = get_observations(); - // Return - return; -} + } // endif: there was no observation in the container + // Get the incube filename + std::string incube = (*this)["incube"].filename(); -/***********************************************************************//** - * @brief Get observation definition - * - * Get observation definition from the user parameters. - ***************************************************************************/ -void ctpsfcube::get_obs(void) -{ - // Get input filename - std::string filename = (*this)["infile"].filename(); + // Get additional binning parameters + double amax = (*this)["amax"].real(); + int anumbins = (*this)["anumbins"].integer(); - // Try first to open as FITS file - try { + // Check for filename validity + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { - // Allocate CTA observation - GCTAObservation obs; - - // Load input file in CTA observation - obs.load(filename); + // Create an event cube based on task parameters + GCTAEventCube cube = create_cube(m_obs); - // Append CTA observation to container - m_obs.append(obs); + // Define psf cube + m_psfcube = GCTAMeanPsf(cube, amax, anumbins); } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - // Load observations from XML file - m_obs.load(filename); + // ... otherwise setup the exposure cube from the counts map + else { - } + // Load event cube from filename + GCTAEventCube cube(incube); - // Return - return; -} + // Define psf cube + m_psfcube = GCTAMeanPsf(cube, amax, anumbins); + } // endelse: cube loaded from file -/***********************************************************************//** - * @brief Set PSF cube definition from counts map - * - * @exception GException::invalid_value - * Invalid counts map projection or invalid events encountered. - * - * Set PSF cube definition from counts map. - ***************************************************************************/ -void ctpsfcube::set_from_cntmap(const std::string& filename) -{ - // Allocate CTA observation - GCTAObservation obs; - - // Load counts map in CTA observation - obs.load(filename); - - // Set PSF cube from counts map - const GCTAEventCube* cube = dynamic_cast(obs.events()); - - // Continue only if cube is valid - if (cube != NULL) { - - // Get sky map projection - const GWcs* wcs = dynamic_cast(cube->map().projection()); - - // Continue only if projection is valid - if (wcs != NULL) { - - // Get user parameters for counts map definition - std::string proj = wcs->code(); - std::string coordsys = wcs->coordsys(); - double xref = wcs->crval(0); - double yref = wcs->crval(1); - double dx = wcs->cdelt(0); - double dy = wcs->cdelt(1); - int nx = cube->map().nx(); - int ny = cube->map().ny(); - double amax = (*this)["amax"].real(); - int anumbins = (*this)["anumbins"].integer(); - - // Get energy definition - m_ebounds = cube->ebounds(); - - // Define PSF cube - m_psfcube = GCTAMeanPsf(proj, coordsys, xref, yref, - dx, dy, nx, ny, - m_ebounds, amax, anumbins); - - } // endif: WCS projection was valid - - // ... projection is not of WCS type - else { - std::string msg = "Counts map project is not of WCS type."; - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); - } - - } // endif: observation contained an events cube + // Read energy dispersion flag + m_apply_edisp = (*this)["edisp"].boolean(); - // ... there is not events cube - else { - std::string msg = "No events cube found in file \"" - ""+filename+"\"."; - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); + // Read output filename (if needed) + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); } // Return diff --git a/src/ctpsfcube/ctpsfcube.hpp b/src/ctpsfcube/ctpsfcube.hpp index f626aac..2cac0b7 100644 --- a/src/ctpsfcube/ctpsfcube.hpp +++ b/src/ctpsfcube/ctpsfcube.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTPSFCUBE_NAME "ctpsfcube" -#define CTPSFCUBE_VERSION "00-01-00" +#define CTPSFCUBE_VERSION "00-02-00" /***********************************************************************//** @@ -67,17 +67,14 @@ protected: void copy_members(const ctpsfcube& app); void free_members(void); void get_parameters(void); - void get_obs(void); - void set_from_cntmap(const std::string& filename); // User parameters - std::string m_outfile; //!< Output PSF cube file + std::string m_outcube; //!< Output PSF cube file bool m_apply_edisp; //!< Apply energy dispersion? // Protected members GObservations m_obs; //!< Observation container GCTAMeanPsf m_psfcube; //!< PSF cube - GEbounds m_ebounds; //!< Energy boundaries }; diff --git a/src/ctpsfcube/ctpsfcube.par b/src/ctpsfcube/ctpsfcube.par index a9520fd..91a687a 100644 --- a/src/ctpsfcube/ctpsfcube.par +++ b/src/ctpsfcube/ctpsfcube.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctpsfcube.par # -# Version: 00-01-00 # +# Version: 00-02-00 # # # # Author: Chia-Chun Lu # # c.lu@mpi-hd.mpg.de # @@ -15,11 +15,11 @@ # # ############################################################################# # -# General parameters +# Observational parameters #=================== -infile, f, a, "NONE",,, "Event list or observation definition file" -cntmap, f, a, "NONE",,, "Counts map for exposure cube definition" -outfile, f, a, "psfcube.fits",,, "Output exposure cube file" +inobs, f, a, "NONE",,, "Event list or observation definition file" +incube, f, a, "NONE",,, "Counts map for psf cube definition" +outcube, f, a, "psfcube.fits",,, "Output psf cube file" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" @@ -34,6 +34,7 @@ ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" # # Spatial binning parameters #=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 4dd0400..08646be 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctselect - Data selection tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -387,8 +387,7 @@ void ctselect::save(void) void ctselect::init_members(void) { // Initialise parameters - m_infile.clear(); - m_outfile.clear(); + m_outobs.clear(); m_prefix.clear(); m_usepnt = false; m_ra = -1.0; @@ -406,12 +405,6 @@ void ctselect::init_members(void) m_infiles.clear(); m_timemin.clear(); m_timemax.clear(); - m_use_xml = false; - - // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD, - // which is defined in GCTALib.hpp. This is somehow a kluge. We need - // a better mechanism to implement the CTA reference MJD. - m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL"); // Return return; @@ -426,8 +419,7 @@ void ctselect::init_members(void) void ctselect::copy_members(const ctselect& app) { // Copy parameters - m_infile = app.m_infile; - m_outfile = app.m_outfile; + m_outobs = app.m_outobs; m_prefix = app.m_prefix; m_usepnt = app.m_usepnt; m_ra = app.m_ra; @@ -445,8 +437,6 @@ void ctselect::copy_members(const ctselect& app) m_infiles = app.m_infiles; m_timemin = app.m_timemin; m_timemax = app.m_timemax; - m_use_xml = app.m_use_xml; - m_cta_ref = app.m_cta_ref; // Return return; @@ -466,7 +456,7 @@ void ctselect::free_members(void) /***********************************************************************//** * @brief Get application parameters * - * Get all task parameters from parameter file or (if required) by querying + * Get all user parameters from parameter file or (if required) by querying * the user. Times are assumed to be in the native CTA MJD format. * * This method also loads observations if no observations are yet allocated. @@ -475,40 +465,15 @@ void ctselect::free_members(void) ***************************************************************************/ void ctselect::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Get CTA event list file name - m_infile = (*this)["infile"].filename(); - - // Allocate CTA observation - GCTAObservation obs; - - // Try first to open as FITS file - try { - - // Load event list in CTA observation - obs.load(m_infile); - - // Append CTA observation to container - m_obs.append(obs); - - // Signal that no XML file should be used for storage - m_use_xml = false; - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Load observations from XML file - m_obs.load(m_infile); - - // Signal that XML file should be used for storage - m_use_xml = true; - - } + // Get observation container without response (not needed) + m_obs = get_observations(false); } // endif: there was no observation in the container @@ -537,8 +502,8 @@ void ctselect::get_parameters(void) // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - m_prefix = (*this)["prefix"].string(); + m_outobs = (*this)["outobs"].filename(); + m_prefix = (*this)["prefix"].string(); } // Set time interval with input times given in CTA reference @@ -1048,7 +1013,7 @@ std::string ctselect::check_infile(const std::string& filename) const if (missing.size() > 0) { message = "The following columns are missing in the" " \"EVENTS\" extension of input file \"" - + m_outfile + "\": "; + + m_outobs + "\": "; for (int i = 0; i < missing.size(); ++i) { message += "\"" + missing[i] + "\""; if (i < missing.size()-1) { @@ -1060,7 +1025,7 @@ std::string ctselect::check_infile(const std::string& filename) const } catch (GException::fits_hdu_not_found& e) { message = "No \"EVENTS\" extension found in input file \"" - + m_outfile + "\"."; + + m_outobs + "\"."; } // Return @@ -1103,13 +1068,13 @@ std::string ctselect::set_outfile_name(const std::string& filename) const void ctselect::save_fits(void) { // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outobs = (*this)["outobs"].filename(); // Get CTA observation from observation container GCTAObservation* obs = dynamic_cast(m_obs[0]); // Save event list - save_event_list(obs, m_infiles[0], m_outfile); + save_event_list(obs, m_infiles[0], m_outobs); // Return return; @@ -1129,14 +1094,14 @@ void ctselect::save_fits(void) void ctselect::save_xml(void) { // Get output filename and prefix - m_outfile = (*this)["outfile"].filename(); + m_outobs = (*this)["outobs"].filename(); m_prefix = (*this)["prefix"].string(); // Issue warning if output filename has no .xml suffix - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4)); + std::string suffix = gammalib::tolower(m_outobs.substr(m_outobs.length()-4,4)); if (suffix != ".xml") { log << "*** WARNING: Name of observation definition output file \""+ - m_outfile+"\"" << std::endl; + m_outobs+"\"" << std::endl; log << "*** WARNING: does not terminate with \".xml\"." << std::endl; log << "*** WARNING: This is not an error, but might be misleading." " It is recommended" << std::endl; @@ -1167,7 +1132,7 @@ void ctselect::save_xml(void) } // endfor: looped over observations // Save observations in XML file - m_obs.save(m_outfile); + m_obs.save(m_outobs); // Return return; diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index 586bf88..5932eab 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctselect - Data selection tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -36,7 +36,7 @@ /* __Definitions _________________________________________________________ */ #define CTSELECT_NAME "ctselect" -#define CTSELECT_VERSION "00-06-00" +#define CTSELECT_VERSION "00-07-00" /***********************************************************************//** @@ -82,8 +82,7 @@ protected: const std::string& outfile) const; // User parameters - std::string m_infile; //!< Input event list or XML file - std::string m_outfile; //!< Output event list or XML file + std::string m_outobs; //!< Output event list or XML file std::string m_prefix; //!< Prefix for multiple event lists bool m_usepnt; //!< Use pointing instead of RA/DEC parameters double m_ra; //!< RA of ROI centre @@ -101,8 +100,7 @@ protected: std::vector m_infiles; //!< Input event filenames GTime m_timemin; //!< Earliest time GTime m_timemax; //!< Latest time - bool m_use_xml; //!< Use XML file instead of FITS file - GTimeReference m_cta_ref; //!< CTA time reference + }; diff --git a/src/ctselect/ctselect.par b/src/ctselect/ctselect.par index b9a50d0..adb8e22 100644 --- a/src/ctselect/ctselect.par +++ b/src/ctselect/ctselect.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctselect.par # -# Version: 00-05-00 # +# Version: 00-06-00 # # # # Author: Jurgen Knodlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -17,9 +17,11 @@ # # File information #================= -infile, f, a, "events.fits",,,"Input event list or observation definition file" -outfile,f, a, "selected_events.fits",,,"Output event list or observation definition file" +inobs, f, a, "events.fits",,,"Input event list or observation definition file" +outobs, f, a, "selected_events.fits",,,"Output event list or observation definition file" prefix, s, h, "selected_",,,"Prefix for event lists in observation definition file" +caldb, s, h, "dummy",,, "Calibration database" +irf, s, h, "cta_dummy_irf",,, "Instrument response function" # # Selection parameters #===================== diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 898adc6..69bd7d2 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctskymap - Sky mapping tool * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2014 by Juergen Knoedlseder * + * copyright (C) 2011-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -35,6 +35,7 @@ /* __ Method name definitions ____________________________________________ */ #define G_INIT_MAP "ctskymap::init_map(GCTAObservation* obs)" #define G_BIN_EVENTS "ctskymap::bin_events(GCTAObservation*)" +#define G_GET_PARAMETERS "ctskymap::get_parameters()" /* __ Debug definitions __________________________________________________ */ @@ -262,12 +263,6 @@ void ctskymap::run(void) } } - // If this is the first valid observation we're working on - // then initialise the sky maps - if (num_obs == 0) { - init_map(obs); - } - // Map events into sky map map_events(obs); @@ -311,10 +306,10 @@ void ctskymap::save(void) } // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outmap = (*this)["outmap"].filename(); // Save sky map - m_skymap.save(m_outfile, clobber()); + m_skymap.save(m_outmap, clobber()); // Return return; @@ -332,39 +327,28 @@ void ctskymap::save(void) ***************************************************************************/ void ctskymap::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Get name of CTA events file - m_evfile = (*this)["evfile"].filename(); + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Load unbinned CTA observation - GCTAObservation obs; - obs.load(m_evfile); - - // Append CTA observation to container - m_obs.append(obs); - - // Use the xref and yref parameters for binning (otherwise the - // pointing direction(s) is/are used) - m_xref = (*this)["xref"].real(); - m_yref = (*this)["yref"].real(); + // Get observation container without response (not needed) + m_obs = get_observations(false); } // endif: there was no observation in the container + // Create sky map based on task parameters + m_skymap = create_map(m_obs); + // Get remaining parameters - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - m_proj = (*this)["proj"].string(); - m_coordsys = (*this)["coordsys"].string(); - m_binsz = (*this)["binsz"].real(); - m_nxpix = (*this)["nxpix"].integer(); - m_nypix = (*this)["nypix"].integer(); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); // Read ahead parameters if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); + m_outmap = (*this)["outmap"].filename(); } // Return @@ -373,53 +357,6 @@ void ctskymap::get_parameters(void) /***********************************************************************//** - * @brief Initialise sky map - * - * @param[in] obs CTA observation (no NULL pointer allowed). - * - * This method initialises the sky map. - ***************************************************************************/ -void ctskymap::init_map(GCTAObservation* obs) -{ - // Clear any existing sky map - m_skymap.clear(); - - // Get map centre. If no centre is specified then extract the map centre - // from the pointing direction. This obviously only works for - double xref; - double yref; - if (m_xref != 9999.0 && m_yref != 9999.0) { - xref = m_xref; - yref = m_yref; - } - else { - - - // Get pointer on CTA pointing - const GCTAPointing& pnt = obs->pointing(); - - // Set reference point to pointing - if (gammalib::toupper(m_coordsys) == "GAL") { - xref = pnt.dir().l_deg(); - yref = pnt.dir().b_deg(); - } - else { - xref = pnt.dir().ra_deg(); - yref = pnt.dir().dec_deg(); - } - - } // endelse: map centre set to pointing - - // Create skymap - m_skymap = GSkymap(m_proj, m_coordsys, xref, yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, 1); - - // Return - return; -} - - -/***********************************************************************//** * @brief Map events into a sky map * * @param[in] obs CTA observation. @@ -473,8 +410,8 @@ void ctskymap::map_events(GCTAObservation* obs) GSkyPixel pixel = m_skymap.dir2pix(dir); // Skip if pixel is out of range - if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) || - pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) { + if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) || + pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) { num_outside_map++; continue; } @@ -525,19 +462,10 @@ void ctskymap::map_events(GCTAObservation* obs) void ctskymap::init_members(void) { // Initialise members - m_evfile.clear(); - m_outfile.clear(); - m_proj.clear(); - m_coordsys.clear(); m_obs.clear(); m_skymap.clear(); - m_emin = 0.0; - m_emax = 0.0; - m_xref = 9999.0; // Flags unset (use pointing direction) - m_yref = 9999.0; // Flags unset (use pointing direction) - m_binsz = 0.0; - m_nxpix = 0; - m_nypix = 0; + m_emin = 0.0; + m_emax = 0.0; // Return return; @@ -552,19 +480,10 @@ void ctskymap::init_members(void) void ctskymap::copy_members(const ctskymap& app) { // Copy attributes - m_evfile = app.m_evfile; - m_outfile = app.m_outfile; - m_proj = app.m_proj; - m_coordsys = app.m_coordsys; - m_obs = app.m_obs; - m_skymap = app.m_skymap; - m_emin = app.m_emin; - m_emax = app.m_emax; - m_xref = app.m_xref; - m_yref = app.m_yref; - m_binsz = app.m_binsz; - m_nxpix = app.m_nxpix; - m_nypix = app.m_nypix; + m_obs = app.m_obs; + m_skymap = app.m_skymap; + m_emin = app.m_emin; + m_emax = app.m_emax; // Return return; diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index de00c4b..b52fcff 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctskymap - Sky mapping tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -69,21 +69,13 @@ protected: void copy_members(const ctskymap& app); void free_members(void); void get_parameters(void); - void init_map(GCTAObservation* obs); void map_events(GCTAObservation* obs); // User parameters - std::string m_evfile; //!< Input event list - std::string m_outfile; //!< Output counts map + std::string m_outmap; //!< Output counts map double m_emin; //!< Minimum energy (TeV) double m_emax; //!< Maximum energy (TeV) - std::string m_proj; //!< WCS projection - std::string m_coordsys; //!< Coordinate system - double m_xref; //!< Longitude reference coordinate - double m_yref; //!< Latitude reference coordinate - double m_binsz; //!< Pixel size - int m_nxpix; //!< Number of pixels in longitude - int m_nypix; //!< Number of pixels in latitude + // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par index 39ebe13..8c7adb9 100644 --- a/src/ctskymap/ctskymap.par +++ b/src/ctskymap/ctskymap.par @@ -17,8 +17,8 @@ # # General parameters #=================== -evfile, f, a, "events.fits",,,"Event data file name" -outfile, f, a, "skymap.fits",,, "Output file name" +inobs, f, a, "events.fits",,,"Event data file name" +outmap, f, a, "skymap.fits",,, "Output file name" # # Event selection parameters #=========================== @@ -27,6 +27,7 @@ emax, r, a, 100.0,,,"Maximum energy in TeV" # # Spatial binning parameters #=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index 1248334..4a17dc4 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * cttsmap - TS map calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Michael Mayer * + * copyright (C) 2014-2015 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -233,8 +233,6 @@ void cttsmap::run(void) log.header1("Initialise TS map"); } - // Initialise maps - init_maps(); // Get bins to be computed // To split the computation on several jobs, @@ -369,7 +367,7 @@ void cttsmap::save(void) } // Get output filename - m_outfile = (*this)["outfile"].filename(); + m_outmap = (*this)["outmap"].filename(); // Create fits file instance GFits fitsfile = GFits(); @@ -392,7 +390,7 @@ void cttsmap::save(void) } // Save FITS file - fitsfile.saveto(m_outfile,true); + fitsfile.saveto(m_outmap,true); // Return return; @@ -411,15 +409,8 @@ void cttsmap::save(void) void cttsmap::init_members(void) { // Initialise members - m_infile.clear(); - m_outfile.clear(); - m_proj.clear(); - m_coordsys.clear(); - m_xref = 0.0; - m_yref = 0.0; - m_binsz = 0.0; - m_nxpix = 0; - m_nypix = 0; + m_srcname.clear(); + m_outmap.clear(); // Initialise protected members m_obs.clear(); @@ -445,25 +436,18 @@ void cttsmap::init_members(void) void cttsmap::copy_members(const cttsmap& app) { // Copy attributes - m_infile = app.m_infile; - m_outfile = app.m_outfile; - m_proj = app.m_proj; - m_coordsys = app.m_coordsys; - m_xref = app.m_xref; - m_yref = app.m_yref; - m_binsz = app.m_binsz; - m_nxpix = app.m_nxpix; - m_nypix = app.m_nypix; + m_srcname = app.m_srcname; + m_outmap = app.m_outmap; // Copy protected members - m_binmin = app.m_binmin; - m_binmax = app.m_binmax; - m_logL0 = app.m_logL0; - m_obs = app.m_obs; - m_tsmap = app.m_tsmap; - m_mapnames = app.m_mapnames; - m_maps = app.m_maps; - m_statusmap = app.m_statusmap; + m_binmin = app.m_binmin; + m_binmax = app.m_binmax; + m_logL0 = app.m_logL0; + m_obs = app.m_obs; + m_tsmap = app.m_tsmap; + m_mapnames = app.m_mapnames; + m_maps = app.m_maps; + m_statusmap = app.m_statusmap; // Clone protected members m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL; @@ -498,42 +482,15 @@ void cttsmap::free_members(void) ***************************************************************************/ void cttsmap::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Allocate CTA observation - GCTAObservation obs; + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - // Get event file name - std::string filename = (*this)["infile"].filename(); - - // Try first to open as FITS file - try { - - // Load data - obs.load(filename); - - // Set response - set_obs_response(&obs); - - // Append observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(filename); - - // Check if all observations have response information. If - // not, get the calibration database parameters and set - // the response properly - set_response(m_obs); - - } // endcatch: file was an XML file + // Get observation container + m_obs = get_observations(); } // endif: there was no observation in the container @@ -542,7 +499,7 @@ void cttsmap::get_parameters(void) if (m_obs.models().size() == 0) { // Get models XML filename - std::string filename = (*this)["srcmdl"].filename(); + std::string filename = (*this)["inmodel"].filename(); // Setup models for optimizing. m_obs.models(GModels(filename)); @@ -579,14 +536,11 @@ void cttsmap::get_parameters(void) (*m_testsource)["DEC"].fix(); } - // Get map parameters - m_xref = (*this)["xref"].real(); - m_yref = (*this)["yref"].real(); - m_proj = (*this)["proj"].string(); - m_coordsys = (*this)["coordsys"].string(); - m_binsz = (*this)["binsz"].real(); - m_nxpix = (*this)["nxpix"].integer(); - m_nypix = (*this)["nypix"].integer(); + // Create sky map based on task parameters + GSkymap map = create_map(m_obs); + + // Initialise maps from user parameters + init_maps(map); // Get optional splitting parameters m_binmin = (*this)["binmin"].integer(); @@ -596,7 +550,7 @@ void cttsmap::get_parameters(void) // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); + m_outmap = (*this)["outmap"].filename(); } // Return @@ -609,23 +563,19 @@ void cttsmap::get_parameters(void) * * Initialises skymaps that will contain map information. ***************************************************************************/ -void cttsmap::init_maps(void) +void cttsmap::init_maps(const GSkymap& map) { // Initialise map information m_tsmap.clear(); // Create skymap - m_tsmap = GSkymap(m_proj, m_coordsys, - m_xref, m_yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, 1); + m_tsmap = GSkymap(map); // Initialise map information m_statusmap.clear(); // Create status map - m_statusmap = GSkymap(m_proj, m_coordsys, - m_xref, m_yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, 1); + m_statusmap = GSkymap(map); // Initialise maps of free parameters if (m_testsource != NULL) { @@ -640,7 +590,7 @@ void cttsmap::init_maps(void) // Add sky map for free parameters. Note that push_back will // create a copy. - m_maps.push_back(m_tsmap); + m_maps.push_back(map); // Store parameter name m_mapnames.push_back((*m_testsource)[i].name()); diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp index f686a24..369c565 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * cttsmap - TS map calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Michael Mayer * + * copyright (C) 2014-2015 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTTSMAP_NAME "cttsmap" -#define CTTSMAP_VERSION "00-01-00" +#define CTTSMAP_VERSION "00-02-00" /***********************************************************************//** @@ -81,24 +81,16 @@ protected: void copy_members(const cttsmap& app); void free_members(void); void get_parameters(void); - void init_maps(void); + void init_maps(const GSkymap& map); // User parameters - std::string m_infile; //!< Input file std::string m_srcname; //!< Name of source which is moved around - std::string m_outfile; //!< Output counts map or XML file - std::string m_proj; //!< WCS projection - std::string m_coordsys; //!< Coordinate system - double m_xref; //!< Longitude reference coordinate - double m_yref; //!< Latitude reference coordinate - double m_binsz; //!< Pixel size - int m_nxpix; //!< Number of pixels in longitude - int m_nypix; //!< Number of pixels in latitude + std::string m_outmap; //!< Output counts map or XML file // Parameters to control speed and job splitting - int m_binmin; //!< Map bin number from which computation should start - int m_binmax; //!< Map bin number where map computation should end - double m_logL0; //!< Likelihood value of null hypothesis + int m_binmin; //!< Map bin number from which computation should start + int m_binmax; //!< Map bin number where map computation should end + double m_logL0; //!< Likelihood value of null hypothesis // Protected members GObservations m_obs; //!< Observation container diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 09c5628..eb29ce5 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: cttsmap.par # -# Version: 00-01-00 # +# Version: 00-02-00 # # # # Author: Michael Mayer # # michael.mayer@desy.de # @@ -17,16 +17,17 @@ # # General parameters #=================== -infile, f, a, "events.fits",,,"Input event list, cube or observation definition file" -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" -srcname, s, a, "TestSource",,,"Test source" -caldb, s, a, "dummy",,, "Calibration database" -irf, s, a, "cta_dummy_irf",,, "Instrument response function" -outfile, f, a, "tsmap.fits",,, "Output counts cube" +inobs, f, a, "events.fits",,, "Input event list, cube or observation definition file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" +srcname, s, a, "TestSource",,, "Test source" +caldb, s, a, "dummy",,, "Calibration database" +irf, s, a, "cta_dummy_irf",,, "Instrument response function" +outmap, f, a, "tsmap.fits",,, "Output counts cube" # # Spatial binning parameters #=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index d025f9a..c4b211b 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Juergen Knoedlseder * + * copyright (C) 2014-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -47,7 +47,8 @@ /* __ Method name definitions ____________________________________________ */ -#define G_GET_EBOUNDS "ctool::get_ebounds()" +#define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)" +#define G_CREATE_EBOUNDS "ctool::create_ebounds()" /* __ Debug definitions __________________________________________________ */ @@ -226,6 +227,12 @@ void ctool::init_members(void) { // Initialise members m_read_ahead = false; + m_use_xml = false; + + // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD, + // which is defined in GCTALib.hpp. This is somehow a kluge. We need + // a better mechanism to implement the CTA reference MJD. + m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL"); // Set logger properties log.date(true); @@ -244,6 +251,8 @@ void ctool::copy_members(const ctool& app) { // Copy members m_read_ahead = app.m_read_ahead; + m_use_xml = app.m_use_xml; + m_cta_ref = app.m_cta_ref; // Return return; @@ -266,7 +275,102 @@ void ctool::free_members(void) /***********************************************************************//** - * @brief Get the energy boundaries + * @brief Get observation container + * + * @param[in] get_response Indicates whether response information should + * be loaded (default: true) + * + * Get an observation container according to the user parameters. The method + * supports loading of a individual FITS file or an observation definition + * file in XML format. If the input filename is empty, parameters are read + * to build a CTA observation from scratch. + ***************************************************************************/ +GObservations ctool::get_observations(const bool& get_response) +{ + // Initialise empty observation container + GObservations obs; + + // Get the filename from the input parameters + std::string filename = (*this)["inobs"].filename(); + + // If no observation definition file has been specified then read all + // parameters that are necessary to create an observation from scratch + if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) { + + // Setup a new CTA observation + GCTAObservation cta_obs = create_cta_obs(); + + // Get response if required + if (get_response) { + + // Set response + set_obs_response(&cta_obs); + + } // endif: response was required + + // Append observation to container + obs.append(cta_obs); + + } // endif: filename was "NONE" or "" + + // ... otherwise we have a file name + else { + + // Try first to open as FITS file + try { + + // Allocate empyt CTA observation + GCTAObservation cta_obs; + + // Load data + cta_obs.load(filename); + + // Get response if required + if (get_response) { + + // Set response + set_obs_response(&cta_obs); + + } // endif: response was required + + // Append observation to container + obs.append(cta_obs); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } + + // ... otherwise try to open as XML file + catch (GException::fits_open_error &e) { + + // Load observations from XML file + obs.load(filename); + + // Get response if required + if (get_response) { + + // For all observations that have no response, set the response + // from the task parameters + set_response(obs); + + } // endif: response was required + + // Signal that XML file should be used for storage + m_use_xml = true; + + } // endcatch: file was an XML file + + } + + // Return observation container + return obs; + +} + + +/***********************************************************************//** + * @brief Create energy boundaries from user parameters * * @exception GException::invalid_value * No valid energy boundary extension found. @@ -275,8 +379,14 @@ void ctool::free_members(void) * supports loading of energy boundary information from the EBOUNDS or * ENERGYBINS extension, or setting energy boundaries using a linear * or logarithmical spacing. + * + * The following parameters are read: + * ebinalg - Energy binning algorithm + * emin - Minimum energy (if ebinalg != FILE) + * emax - Maximum energy (if ebinalg != FILE) + * enumbins - Number of energy bins (if ebinalg != FILE) ***************************************************************************/ -GEbounds ctool::get_ebounds(void) +GEbounds ctool::create_ebounds(void) { // Allocate energy boundaries GEbounds ebounds; @@ -310,7 +420,7 @@ GEbounds ctool::get_ebounds(void) "An \"EBOUNDS\" or \"ENERGYBINS\" extension" " is required if the parameter \"ebinalg\"" " is set to \"FILE\"."; - throw GException::invalid_value(G_GET_EBOUNDS, msg); + throw GException::invalid_value(G_CREATE_EBOUNDS, msg); } } // endif: ebinalg was "FILE" @@ -343,12 +453,229 @@ GEbounds ctool::get_ebounds(void) /***********************************************************************//** + * @brief Create a skymap from user parameters + * + * @param[in] obs Observation container. + * + * Creates an empty skymap from input user parameters. The reference + * coordinate is extracted from the mean pointing of the CTA observations + * in the observation container (see get_mean_pointing()). + * + * The following parameters are read: + * usepnt - Use pointing for reference coordinate? + * xref - Reference in Right Ascension or longitude (if usepnt=no) + * yref - Reference in Declination or latitude (if usepnt=no) + * proj - Coordinate projection + * coordsys - Coordinate system + * binsz - Bin size (deg) + * nxpix - Number of pixels in Right Ascension or longitude + * nypix - Number of pixels in Declination or latitude + ***************************************************************************/ +GSkymap ctool::create_map(const GObservations& obs) +{ + // Read task parameters + double xref = 0.0; + double yref = 0.0; + bool usepnt = (*this)["usepnt"].boolean(); + if (!usepnt) { + xref = (*this)["xref"].real(); + yref = (*this)["yref"].real(); + } + std::string proj = (*this)["proj"].string(); + std::string coordsys = (*this)["coordsys"].string(); + double binsz = (*this)["binsz"].real(); + int nxpix = (*this)["nxpix"].integer(); + int nypix = (*this)["nypix"].integer(); + + // If requested, get pointing from observations + if (usepnt) { + + // Get mean pointing + GSkyDir pnt = get_mean_pointing(obs); + + // Use correct coordinate system + if (gammalib::toupper(coordsys) == "GAL") { + xref = pnt.l_deg(); + yref = pnt.b_deg(); + } + else { + xref = pnt.ra_deg(); + yref = pnt.dec_deg(); + } + + } // endif: got mean pointing as reference + + // Initialise sky map + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, + nxpix, nypix, 1); + + // Return sky map + return map; +} + + +/***********************************************************************//** + * @brief Create a CTA event cube from user parameters + * + * @param[in] obs Observation container. + * + * Creates an empty CTA event cube from input user parameters. The reference + * coordinate is extracted from the mean pointing of the CTA observations + * in the observation container (see get_mean_pointing()). This method will + * be used if a ctool has the input parameter usepnt=true. The method + * appends a dummy GTI to the event cube. + * + * The following parameters are read: + * proj - Coordinate projection + * coordsys - Coordinate system + * binsz - Bin size (deg) + * nxpix - Number of pixels in Right Ascension or longitude + * nypix - Number of pixels in Declination or latitude + * ebinalg - Energy binning algorithm + * emin - Minimum energy (if ebinalg != FILE) + * emax - Maximum energy (if ebinalg != FILE) + * enumbins - Number of energy bins (if ebinalg != FILE) + ***************************************************************************/ +GCTAEventCube ctool::create_cube(const GObservations& obs) +{ + // Get skymap + GSkymap map = create_map(obs); + + // Set energy boundaries + GEbounds ebounds = create_ebounds(); + + // Set dummy GTI that is needed for event cube creation + GGti gti; + gti.append(GTime(0.0), GTime(0.1234)); + + // Extend skymap to the requested number of maps + map.nmaps(ebounds.size()); + + // Initialise event cube from all parameters + GCTAEventCube cube = GCTAEventCube(map, ebounds, gti); + + // Return event cube + return cube; +} + + +/***********************************************************************//** + * @brief Create a CTA observation from user parameters + * + * Creates an empty CTA observation from user parameters. An empty event list + * including RoI, GTI and energy boundary information is attached to the + * observation. The method also sets the pointing direction using the ra and + * dec parameter, the ROI based on ra, dec and rad, a single GTI based on + * tmin and tmax, and a single energy boundary based on emin and emax. The + * method furthermore sets the ontime, livetime and deadtime correction + * factor. + * + * The following parameters are read: + * ra - Right Ascension of pointing and RoI centre (deg) + * dec - Declination of pointing and RoI centre (deg) + * rad - Radius of RoI (deg) + * deadc - Deadtime correction factor + * tmin - Start time + * tmax - Stop time + * emin - Minimum energy (TeV) + * emax - Maximum energy (TeV) + ***************************************************************************/ +GCTAObservation ctool::create_cta_obs(void) +{ + // Get CTA observation parameters + double ra = (*this)["ra"].real(); + double dec = (*this)["dec"].real(); + double rad = (*this)["rad"].real(); + double deadc = (*this)["deadc"].real(); + double t_min = (*this)["tmin"].real(); + double t_max = (*this)["tmax"].real(); + double e_min = (*this)["emin"].real(); + double e_max = (*this)["emax"].real(); + + // Allocate CTA observation and empty event list + GCTAObservation obs; + GCTAEventList list; + + // Set pointing direction + GCTAPointing pnt; + GSkyDir skydir; + skydir.radec_deg(ra, dec); + pnt.dir(skydir); + + // Set ROI + GCTAInstDir instdir(skydir); + GCTARoi roi(instdir, rad); + + // Set GTI + GGti gti(m_cta_ref); + GTime tstart; + GTime tstop; + tstart.set(t_min, m_cta_ref); + tstop.set(t_max, m_cta_ref); + gti.append(tstart, tstop); + + // Set energy boundaries + GEbounds ebounds; + GEnergy emin; + GEnergy emax; + emin.TeV(e_min); + emax.TeV(e_max); + ebounds.append(emin, emax); + + // Set CTA event list attributes + list.roi(roi); + list.gti(gti); + list.ebounds(ebounds); + + // Attach event list to CTA observation + obs.events(list); + + // Set observation ontime, livetime and deadtime correction factor + obs.pointing(pnt); + obs.ontime(gti.ontime()); + obs.livetime(gti.ontime()*deadc); + obs.deadc(deadc); + + // Return CTA observation + return obs; + +} + + +/***********************************************************************//** + * @brief Throws exception in inobs parameter is not valid + * + * @param[in] method Method name. + ***************************************************************************/ +void ctool::require_inobs(const std::string& method) +{ + // Throw exception if no infile is given + if ((*this)["inobs"].filename() == "NONE" || + (*this)["inobs"].filename() == "") { + std::string msg = "A valid file needs to be specified for the " + "\"inobs\" parameter, yet \""+ + (*this)["inobs"].filename()+"\" was given." + " Specify a valid observation definition or " + "FITS file to proceed."; + throw GException::invalid_value(method, msg); + } + + // Return + return; +} + + +/***********************************************************************//** * @brief Set response for all CTA observations in container * * @param[in,out] obs Observation container * * Set the response for a CTA observations in the container that so far have - * no response using the "database" and "irf" task parameters. + * no response using the "caldb" and "irf" task parameters. + * + * The following parameters are read: + * caldb - Calibration database + * irf - Instrument response function ***************************************************************************/ void ctool::set_response(GObservations& obs) { @@ -380,8 +707,12 @@ void ctool::set_response(GObservations& obs) * * @param[in,out] obs CTA observation * - * Set the response for a CTA observation using the "database" and "irf" + * Set the response for a CTA observation using the "caldb" and "irf" * task parameters. + * + * The following parameters are read: + * caldb - Calibration database + * irf - Instrument response function ***************************************************************************/ void ctool::set_obs_response(GCTAObservation* obs) { @@ -409,6 +740,67 @@ void ctool::set_obs_response(GCTAObservation* obs) /***********************************************************************//** + * @brief Derives mean pointing from CTA observations + * + * @param[in] obs Observation container. + * @return Pointing direction. + * + * Computes the mean pointing direction from all CTA observations in the + * observation container. + * + * @todo This method does not work properly for wrap arounds in Right + * Ascension. It certainly will also lead to bad results near the + * celestial pole. + ***************************************************************************/ +GSkyDir ctool::get_mean_pointing(const GObservations& obs) +{ + // Initialise values + double xref = 0.0; + double yref = 0.0; + int n_pnt = 0; + + // Loop over all observations to compute mean pointings + for (int i = 0; i < obs.size(); ++i) { + + // If we have a CTA observation then extract the pointing direction + const GCTAObservation* cta_obs = dynamic_cast(obs[i]); + if (cta_obs != NULL) { + + // Get pointing + const GCTAPointing& pnt = cta_obs->pointing(); + + // Add to coordinates + xref += pnt.dir().ra_deg(); + yref += pnt.dir().dec_deg(); + n_pnt += 1; + + } // endif: cta observation was valid + + } // endfor: loop over all observations + + // Signal if no pointing is found + if (n_pnt < 1) { + std::string msg = "No valid CTA observation has been found in " + "observation list, hence no pointing information " + "could be extracted. Use the \"usepnt=no\" " + "option and specify the pointing explicitly."; + throw GException::invalid_value(G_GET_MEAN_POINTING, msg); + } + + // Calculate mean coordinates + double ra = xref / double(n_pnt); + double dec = yref / double(n_pnt); + + // Initialise sky dir + GSkyDir dir = GSkyDir(); + dir.radec_deg(ra,dec); + + // Return sky direction + return dir; +} + + +/***********************************************************************//** * @brief Get current resident set size (physical memory use) in Bytes * * @return Physical memory use in Bytes. diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 39acdec..ad20618 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Juergen Knoedlseder * + * copyright (C) 2014-2015 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -60,19 +60,37 @@ public: // Public methods virtual void execute(void); +#ifdef SWIG +public: +#else protected: +#endif // Protected methods - void init_members(void); - void copy_members(const ctool& app); - void free_members(void); - const bool& read_ahead(void) const; - GEbounds get_ebounds(void); - void set_response(GObservations& obs); - void set_obs_response(GCTAObservation* obs); - size_t get_current_rss(void); + void init_members(void); + void copy_members(const ctool& app); + void free_members(void); + const bool& read_ahead(void) const; + GObservations get_observations(const bool& get_response = true); + + // Protected methods that create objects from user parameters + GEbounds create_ebounds(void); + GSkymap create_map(const GObservations& obs); + GCTAEventCube create_cube(const GObservations& obs); + GCTAObservation create_cta_obs(void); + + // Protected methods that check user parameters + void require_inobs(const std::string& method); + + // Protected support methods + void set_response(GObservations& obs); + void set_obs_response(GCTAObservation* obs); + GSkyDir get_mean_pointing(const GObservations& obs); + size_t get_current_rss(void); // Protected members - bool m_read_ahead; //!< Read ahead parameters + bool m_read_ahead; //!< Read ahead parameters + bool m_use_xml; //!< Use XML file instead of FITS file for observations + GTimeReference m_cta_ref; //!< CTA time reference }; diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index 9b257fc..c2db50e 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -52,7 +52,7 @@ export PFILES=pfiles $ECHO @ECHO_N@ "Test cspull: " chmod +x @abs_top_srcdir@/scripts/cspull.py @abs_top_srcdir@/scripts/cspull.py \ - srcmdl="data/crab.xml" \ + inmodel="data/crab.xml" \ outfile="pull.dat" \ ntrials="3" \ caldb="irf" \ @@ -138,9 +138,9 @@ $ECHO " ok" $ECHO @ECHO_N@ "Test csresmap: " chmod +x @abs_top_srcdir@/scripts/csresmap.py @abs_top_srcdir@/scripts/csresmap.py \ - infile="data/crab_events.fits" \ - outfile="resmap.fits" \ - srcmdl="data/crab.xml" \ + inobs="data/crab_events.fits" \ + outmap="resmap.fits" \ + inmodel="data/crab.xml" \ caldb="irf" \ irf="cta_dummy_irf" \ nxpix=50 \ diff --git a/test/test_ctbin.py b/test/test_ctbin.py index 3027938..0e15eef 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctbin bin = ctools.ctbin() - bin["evfile"].filename(self.events_name) - bin["outfile"].filename("cntmap.fits") + bin["inobs"].filename(self.events_name) + bin["outcube"].filename("cntmap.fits") bin["ebinalg"].string("LOG") bin["emin"].real(0.1) bin["emax"].real(100.0) diff --git a/test/test_ctbkgcube.py b/test/test_ctbkgcube.py index f17a87d..b2b4d88 100755 --- a/test/test_ctbkgcube.py +++ b/test/test_ctbkgcube.py @@ -40,6 +40,8 @@ class Test(gammalib.GPythonTestSuite): # Set members self.events_name = "data/crab_events.fits" self.bkg_model = "data/crab.xml" + self.caldb = "irf" + self.irf = "cta_dummy_irf" # Return return @@ -65,11 +67,13 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctbkgcube bkgcube = ctools.ctbkgcube() - bkgcube["infile"].filename(self.events_name) + bkgcube["inobs"].filename(self.events_name) bkgcube["inmodel"].filename(self.bkg_model) - bkgcube["cntmap"].filename("NONE") - bkgcube["outfile"].filename("bkgcube.fits") + bkgcube["incube"].filename("NONE") + bkgcube["outcube"].filename("bkgcube.fits") bkgcube["outmodel"].filename("bkgcube.xml") + bkgcube["caldb"].string(self.caldb) + bkgcube["irf"].string(self.irf) bkgcube["ebinalg"].string("LOG") bkgcube["emin"].real(0.1) bkgcube["emax"].real(100.0) diff --git a/test/test_ctbutterfly.py b/test/test_ctbutterfly.py index a3e160b..34aa9c8 100755 --- a/test/test_ctbutterfly.py +++ b/test/test_ctbutterfly.py @@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctbutterfly butterfly = ctools.ctbutterfly() - butterfly["infile"].filename(self.events_name) - butterfly["srcmdl"].filename(self.model_name) + butterfly["inobs"].filename(self.events_name) + butterfly["inmodel"].filename(self.model_name) butterfly["srcname"].string("Crab") butterfly["caldb"].string(self.caldb) butterfly["irf"].string(self.irf) diff --git a/test/test_ctcubemask.py b/test/test_ctcubemask.py index 300ab5b..01c426c 100755 --- a/test/test_ctcubemask.py +++ b/test/test_ctcubemask.py @@ -65,9 +65,9 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctcubemask mask = ctools.ctcubemask() - mask["infile"].filename(self.cntmap_name) + mask["inobs"].filename(self.cntmap_name) mask["regfile"].filename(self.regfile) - mask["outfile"].filename("filtered_cntmap.fits") + mask["outcube"].filename("filtered_cntmap.fits") mask["ra"].real(83.63) mask["dec"].real(22.01) mask["rad"].real(2.0) diff --git a/test/test_ctexpcube.py b/test/test_ctexpcube.py index d2a55b1..6c376de 100755 --- a/test/test_ctexpcube.py +++ b/test/test_ctexpcube.py @@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctexpcube expcube = ctools.ctexpcube() - expcube["infile"].filename(self.events_name) - expcube["cntmap"].filename("NONE") - expcube["outfile"].filename("expcube.fits") + expcube["inobs"].filename(self.events_name) + expcube["incube"].filename("NONE") + expcube["outcube"].filename("expcube.fits") expcube["caldb"].string(self.caldb) expcube["irf"].string(self.irf) expcube["ebinalg"].string("LOG") diff --git a/test/test_ctlike.py b/test/test_ctlike.py index 1f743fa..8618380 100755 --- a/test/test_ctlike.py +++ b/test/test_ctlike.py @@ -67,11 +67,11 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctlike like = ctools.ctlike() - like["infile"].filename(self.events_name) - like["srcmdl"].filename(self.model_name) + like["inobs"].filename(self.events_name) + like["inmodel"].filename(self.model_name) like["caldb"].string(self.caldb) like["irf"].string(self.irf) - like["outmdl"].filename("result.xml") + like["outmodel"].filename("result.xml") # Run tool self.test_try("Run ctlike") diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 2b7ba3d..41ce959 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -66,12 +66,13 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctmodel from scratch model = ctools.ctmodel() - model["infile"].filename("NONE") - model["outfile"].filename("modmap.fits") - model["srcmdl"].filename(self.model_name) - model["obsfile"].filename("NONE") + model["incube"].filename("NONE") + model["outcube"].filename("modmap.fits") + model["inmodel"].filename(self.model_name) + model["inobs"].filename("NONE") model["caldb"].string(self.caldb) model["irf"].string(self.irf) + model["rad"].real(5.0) model["ra"].real(83.63) model["dec"].real(22.01) model["tmin"].real(0.0) diff --git a/test/test_ctobssim.py b/test/test_ctobssim.py index fe092b1..43c47de 100755 --- a/test/test_ctobssim.py +++ b/test/test_ctobssim.py @@ -67,8 +67,8 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctobssim sim = ctools.ctobssim() - sim["infile"].filename(self.model_name) - sim["outfile"].filename("events.fits") + sim["inmodel"].filename(self.model_name) + sim["outevents"].filename("events.fits") sim["caldb"].string(self.caldb) sim["irf"].string(self.irf) sim["ra"].real(83.63) @@ -117,8 +117,8 @@ class Test(gammalib.GPythonTestSuite): # Set-up ctobssim sim = ctools.ctobssim(obs) - sim["outfile"].filename("sim_events.xml") - + sim["outevents"].filename("sim_events.xml") + sim["inmodel"].filename(self.model_name) # Run tool self.test_try("Run ctobssim") try: diff --git a/test/test_ctools.sh b/test/test_ctools.sh index e7c4e99..e36b054 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -67,10 +67,11 @@ export PFILES=pfiles # Test ctobssim # ============= $ECHO -n "Test ctobssim: " + # # Run 1 -$ctobssim infile="data/crab.xml" \ - outfile="events.fits" \ +$ctobssim inmodel="data/crab.xml" \ + outevents="events.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ ra=83.63 \ @@ -90,8 +91,8 @@ else fi # # Run 2 -$ctobssim infile="data/model_background_cube.xml" \ - outfile="events2.fits" \ +$ctobssim inmodel="data/model_background_cube.xml" \ + outevents="events2.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ ra=84.17263 \ @@ -116,8 +117,8 @@ $ECHO " ok" # Test ctskymap # ============= $ECHO -n "Test ctskymap: " -$ctskymap evfile="events.fits" \ - outfile="skymap.fits" \ +$ctskymap inobs="events.fits" \ + outmap="skymap.fits" \ emin=0.1 \ emax=100.0 \ nxpix=200 \ @@ -144,8 +145,8 @@ $ECHO " ok" $ECHO -n "Test ctbin: " # # Run 1 -$ctbin evfile="events.fits" \ - outfile="cntmap.fits" \ +$ctbin inobs="events.fits" \ + outcube="cntmap.fits" \ emin=0.1 \ emax=100.0 \ enumbins=20 \ @@ -167,8 +168,8 @@ else fi # # Run 2 -$ctbin evfile="events2.fits" \ - outfile="cntmap2.fits" \ +$ctbin inobs="events2.fits" \ + outcube="cntmap2.fits" \ emin=0.3 \ emax=8.0 \ enumbins=10 \ @@ -191,18 +192,19 @@ fi $ECHO " ok" + # # Test ctmodel # ============ $ECHO -n "Test ctmodel: " # # Run 1 -$ctmodel infile="cntmap.fits" \ - obsfile="cntmap.fits" \ - outfile="modmap1.fits" \ +$ctmodel inobs="cntmap.fits" \ + incube="cntmap.fits" \ + outcube="modmap1.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ - srcmdl="data/crab.xml" + inmodel="data/crab.xml" $ECHO -n "." if [ -s "modmap1.fits" ] then @@ -213,14 +215,15 @@ else fi # # Run 3 -$ctmodel infile="NONE" \ - obsfile="NONE" \ - outfile="modmap3.fits" \ +$ctmodel inobs="NONE" \ + incube="NONE" \ + outcube="modmap3.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ - srcmdl="data/crab.xml" \ + inmodel="data/crab.xml" \ ra=83.63 \ dec=22.01 \ + rad=10.0 \ tmin=0.0 \ tmax=1800.0 \ emin=0.1 \ @@ -243,12 +246,12 @@ else fi # # Run 2 -$ctmodel infile="cntmap2.fits" \ - obsfile="cntmap2.fits" \ - outfile="modmap2.fits" \ +$ctmodel inobs="cntmap2.fits" \ + incube="cntmap2.fits" \ + outcube="modmap2.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ - srcmdl="data/model_background_cube.xml" + inmodel="data/model_background_cube.xml" $ECHO -n "." if [ -s "modmap2.fits" ] then @@ -264,8 +267,8 @@ $ECHO " ok" # Test ctselect # ============= $ECHO -n "Test ctselect: " -$ctselect infile="events.fits" \ - outfile="selected_events.fits" \ +$ctselect inobs="events.fits" \ + outobs="selected_events.fits" \ ra=83.63 \ dec=22.01 \ rad=10.0 \ @@ -290,9 +293,9 @@ $ECHO " ok" $ECHO -n "Test ctlike: " # # Run 1 -$ctlike infile="cntmap.fits" \ - srcmdl="data/crab.xml" \ - outmdl="results_binned.xml" \ +$ctlike inobs="cntmap.fits" \ + inmodel="data/crab.xml" \ + outmodel="results_binned.xml" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." @@ -305,9 +308,9 @@ else fi # # Run 2 -$ctlike infile="selected_events.fits" \ - srcmdl="data/crab.xml" \ - outmdl="results_unbinned.xml" \ +$ctlike inobs="selected_events.fits" \ + inmodel="data/crab.xml" \ + outmodel="results_unbinned.xml" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." @@ -320,9 +323,9 @@ else fi # # Run 3 -$ctlike infile="events2.fits" \ - srcmdl="data/model_background_cube.xml" \ - outmdl="results_unbinned_background_cube.xml" \ +$ctlike inobs="events2.fits" \ + inmodel="data/model_background_cube.xml" \ + outmodel="results_unbinned_background_cube.xml" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." @@ -335,9 +338,9 @@ else fi # # Run 4 -$ctlike infile="cntmap2.fits" \ - srcmdl="data/model_background_cube.xml" \ - outmdl="results_binned_background_cube.xml" \ +$ctlike inobs="cntmap2.fits" \ + inmodel="data/model_background_cube.xml" \ + outmodel="results_binned_background_cube.xml" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." @@ -355,12 +358,12 @@ $ECHO " ok" # Test cttsmap # ============ $ECHO -n "Test cttsmap: " -$cttsmap infile="selected_events.fits" \ - srcmdl="data/crab.xml" \ +$cttsmap inobs="selected_events.fits" \ + inmodel="data/crab.xml" \ srcname="Crab" \ caldb="irf" \ irf="cta_dummy_irf" \ - outfile="tsmap.fits" \ + outmap="tsmap.fits" \ nxpix=5 \ nypix=5 \ binsz=0.05 \ @@ -385,9 +388,9 @@ $ECHO " ok" $ECHO -n "Test ctexpcube: " # # Run 1 -$ctexpcube infile="data/crab_events.fits" \ - cntmap="NONE" \ - outfile="expcube1.fits" \ +$ctexpcube inobs="data/crab_events.fits" \ + incube="NONE" \ + outcube="expcube1.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ emin=0.1 \ @@ -410,9 +413,9 @@ else fi # # Run 2 -$ctexpcube infile="data/crab_events.fits" \ - cntmap="data/crab_cntmap.fits" \ - outfile="expcube2.fits" \ +$ctexpcube inobs="data/crab_events.fits" \ + incube="data/crab_cntmap.fits" \ + outcube="expcube2.fits" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." @@ -432,9 +435,9 @@ $ECHO " ok" $ECHO -n "Test ctpsfcube: " # # Run 1 -$ctpsfcube infile="data/crab_events.fits" \ - cntmap="NONE" \ - outfile="psfcube1.fits" \ +$ctpsfcube inobs="data/crab_events.fits" \ + incube="NONE" \ + outcube="psfcube1.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ emin=0.1 \ @@ -459,9 +462,9 @@ else fi # # Run 2 -$ctpsfcube infile="data/crab_events.fits" \ - cntmap="data/crab_cntmap.fits" \ - outfile="psfcube2.fits" \ +$ctpsfcube inobs="data/crab_events.fits" \ + incube="data/crab_cntmap.fits" \ + outcube="psfcube2.fits" \ caldb="irf" \ irf="cta_dummy_irf" \ amax=0.3 \ @@ -483,11 +486,13 @@ $ECHO " ok" $ECHO -n "Test ctbkgcube: " # # Run 1 -$ctbkgcube infile="data/crab_events.fits" \ +$ctbkgcube inobs="data/crab_events.fits" \ inmodel="data/crab.xml" \ - cntmap="NONE" \ - outfile="bkgcube1.fits" \ + incube="NONE" \ + outcube="bkgcube1.fits" \ outmodel="bkgcube1.xml" \ + caldb="irf" \ + irf="cta_dummy_irf" \ emin=0.1 \ emax=100.0 \ enumbins=20 \ @@ -515,10 +520,12 @@ else fi # # Run 2 -$ctbkgcube infile="data/crab_events.fits" \ +$ctbkgcube inobs="data/crab_events.fits" \ inmodel="data/crab.xml" \ - cntmap="data/crab_cntmap.fits" \ - outfile="bkgcube2.fits" \ + incube="data/crab_cntmap.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + outcube="bkgcube2.fits" \ outmodel="bkgcube2.xml" $ECHO -n "." if [ -s "bkgcube2.fits" ] @@ -542,9 +549,9 @@ $ECHO " ok" # Test ctcubemask # =============== $ECHO -n "Test ctcubemask: " -$ctcubemask infile="data/crab_cntmap.fits" \ +$ctcubemask inobs="data/crab_cntmap.fits" \ regfile="data/exclusion.reg" \ - outfile="filtered_cube.fits" \ + outcube="filtered_cube.fits" \ ra=83.63 \ dec=22.01 \ rad=2.0 \ @@ -565,8 +572,8 @@ $ECHO " ok" # Test ctbutterfly # =============== $ECHO -n "Test ctbutterfly: " -$ctbutterfly infile="data/crab_events.fits.gz" \ - srcmdl="data/crab.xml" \ +$ctbutterfly inobs="data/crab_events.fits.gz" \ + inmodel="data/crab.xml" \ srcname="Crab" \ outfile="butterfly.txt" \ caldb="irf" \ diff --git a/test/test_ctpsfcube.py b/test/test_ctpsfcube.py index 5870619..1cd3b89 100755 --- a/test/test_ctpsfcube.py +++ b/test/test_ctpsfcube.py @@ -66,9 +66,9 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctpsfcube psfcube = ctools.ctpsfcube() - psfcube["infile"].filename(self.events_name) - psfcube["cntmap"].filename("NONE") - psfcube["outfile"].filename("psfcube.fits") + psfcube["inobs"].filename(self.events_name) + psfcube["incube"].filename("NONE") + psfcube["outcube"].filename("psfcube.fits") psfcube["caldb"].string(self.caldb) psfcube["irf"].string(self.irf) psfcube["ebinalg"].string("LOG") diff --git a/test/test_ctselect.py b/test/test_ctselect.py index 3fbd776..8f24a85 100755 --- a/test/test_ctselect.py +++ b/test/test_ctselect.py @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctselect select = ctools.ctselect() - select["infile"].filename(self.events_name) - select["outfile"].filename("selected_events.fits") + select["inobs"].filename(self.events_name) + select["outobs"].filename("selected_events.fits") select["ra"].real(83.63) select["dec"].real(22.01) select["rad"].real(3.0) diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py index 5dee2cd..5226ed8 100755 --- a/test/test_ctskymap.py +++ b/test/test_ctskymap.py @@ -64,8 +64,8 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up ctskymap skymap = ctools.ctskymap() - skymap["evfile"].filename(self.events_name) - skymap["outfile"].filename("skymap.fits") + skymap["inobs"].filename(self.events_name) + skymap["outmap"].filename("skymap.fits") skymap["emin"].real(0.1) skymap["emax"].real(100.0) skymap["nxpix"].integer(200) diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index ef731a9..1fcba9e 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -67,12 +67,12 @@ class Test(gammalib.GPythonTestSuite): """ # Set-up cttsmap tsmap = ctools.cttsmap() - tsmap["infile"].filename(self.events_name) - tsmap["srcmdl"].filename(self.model_name) + tsmap["inobs"].filename(self.events_name) + tsmap["inmodel"].filename(self.model_name) tsmap["srcname"].string("Crab") tsmap["caldb"].string(self.caldb) tsmap["irf"].string(self.irf) - tsmap["outfile"].filename("tsmap.fits") + tsmap["outmap"].filename("tsmap.fits") tsmap["nxpix"].integer(5) tsmap["nypix"].integer(5) tsmap["binsz"].real(0.02) diff --git a/test/test_pipelines.py b/test/test_pipelines.py index 9af009e..f61692e 100755 --- a/test/test_pipelines.py +++ b/test/test_pipelines.py @@ -78,8 +78,8 @@ class Test(gammalib.GPythonTestSuite): # Simulate events sim = ctools.ctobssim() - sim["infile"].filename(model_name) - sim["outfile"].filename(events_name) + sim["inmodel"].filename(model_name) + sim["outevents"].filename(events_name) sim["caldb"].string(caldb) sim["irf"].string(irf) sim["ra"].real(ra) @@ -98,8 +98,8 @@ class Test(gammalib.GPythonTestSuite): # Select events select = ctools.ctselect() - select["infile"].filename(events_name) - select["outfile"].filename(selected_events_name) + select["inobs"].filename(events_name) + select["outobs"].filename(selected_events_name) select["ra"].real(ra) select["dec"].real(dec) select["rad"].real(rad_select) @@ -116,9 +116,9 @@ class Test(gammalib.GPythonTestSuite): # Perform maximum likelihood fitting like = ctools.ctlike() - like["infile"].filename(selected_events_name) - like["srcmdl"].filename(model_name) - like["outmdl"].filename(result_name) + like["inobs"].filename(selected_events_name) + like["inmodel"].filename(model_name) + like["outmodel"].filename(result_name) like["caldb"].string(caldb) like["irf"].string(irf) self.test_try("Execute ctlike") @@ -149,7 +149,7 @@ class Test(gammalib.GPythonTestSuite): # Simulate events sim = ctools.ctobssim() - sim["infile"].filename(model_name) + sim["inmodel"].filename(model_name) sim["caldb"].string(caldb) sim["irf"].string(irf) sim["ra"].real(ra) diff --git a/test/test_python.py b/test/test_python.py index e63f915..6dd5aaa 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -102,6 +102,7 @@ if __name__ == '__main__': # Run test suite success = suites.run() + # Save test results suites.save("reports/ctools.xml")