From be3e4fa07e18a76816e99234bd10930217919f5e Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Fri, 21 Nov 2014 15:28:33 +0100 Subject: [PATCH] initial commit to extend ctool-base class and unify/simplify the get_parameters()-function of individual ctools --- pyext/ctools/obsutils.py | 1 + scripts/cspull.py | 64 +++++++++++++++++++++++++++++++++++++++++----------------------- scripts/csresmap.py | 17 +++++++---------- scripts/cstsdist.py | 9 ++++++++- src/ctbin/ctbin.cpp | 167 ++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------- src/ctbin/ctbin.hpp | 13 ++----------- src/ctbin/ctbin.par | 6 +++--- src/ctbkgcube/ctbkgcube.cpp | 212 +++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctbkgcube/ctbkgcube.hpp | 7 ++++--- src/ctbkgcube/ctbkgcube.par | 13 ++++++++----- src/ctbutterfly/ctbutterfly.cpp | 50 ++++++++++++++++---------------------------------- src/ctbutterfly/ctbutterfly.hpp | 3 +-- src/ctbutterfly/ctbutterfly.par | 6 +++--- src/ctcubemask/ctcubemask.cpp | 77 +++++++++++++++++++++++++++++++---------------------------------------------- src/ctcubemask/ctcubemask.hpp | 8 ++++---- src/ctcubemask/ctcubemask.par | 6 +++--- src/ctexpcube/ctexpcube.cpp | 188 +++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------- src/ctexpcube/ctexpcube.hpp | 7 ++----- src/ctexpcube/ctexpcube.par | 11 ++++++----- src/ctlike/ctlike.cpp | 69 +++++++++++++++++++++++++++------------------------------------------ src/ctlike/ctlike.hpp | 4 ++-- src/ctlike/ctlike.par | 8 ++++---- src/ctmodel/ctmodel.cpp | 326 ++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ src/ctmodel/ctmodel.hpp | 29 ++--------------------------- src/ctmodel/ctmodel.par | 12 +++++++----- src/ctobssim/ctobssim.cpp | 180 +++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctobssim/ctobssim.hpp | 17 ++--------------- src/ctobssim/ctobssim.par | 7 ++++--- src/ctpsfcube/ctpsfcube.cpp | 198 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------------------------------------------------------------- src/ctpsfcube/ctpsfcube.hpp | 7 ++----- src/ctpsfcube/ctpsfcube.par | 11 ++++++----- src/ctselect/ctselect.cpp | 73 +++++++++++++++++++++++-------------------------------------------------- src/ctselect/ctselect.hpp | 8 +++----- src/ctselect/ctselect.par | 8 +++++--- src/ctskymap/ctskymap.cpp | 134 +++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------- src/ctskymap/ctskymap.hpp | 12 ++---------- src/ctskymap/ctskymap.par | 5 +++-- src/cttsmap/cttsmap.cpp | 116 +++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------- src/cttsmap/cttsmap.hpp | 14 +++----------- src/cttsmap/cttsmap.par | 9 +++++---- src/support/ctool.cpp | 416 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/support/ctool.hpp | 11 +++++++++++ test/test_cscripts.sh.in | 8 ++++---- test/test_ctbin.py | 4 ++-- test/test_ctbkgcube.py | 10 +++++++--- test/test_ctbutterfly.py | 4 ++-- test/test_ctcubemask.py | 4 ++-- test/test_ctexpcube.py | 6 +++--- test/test_ctlike.py | 6 +++--- test/test_ctmodel.py | 9 +++++---- test/test_ctobssim.py | 8 ++++---- test/test_ctools.sh | 129 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------- test/test_ctpsfcube.py | 6 +++--- test/test_ctselect.py | 4 ++-- test/test_ctskymap.py | 4 ++-- test/test_cttsmap.py | 6 +++--- test/test_pipelines.py | 16 ++++++++-------- test/test_python.py | 1 + 58 files changed, 1139 insertions(+), 1625 deletions(-) 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..dc6ebf6 100755 --- a/scripts/cspull.py +++ b/scripts/cspull.py @@ -19,10 +19,12 @@ # # ========================================================================== import gammalib +import ctools +#import obsutils from ctools import obsutils import sys import csv -import math +#import math # ============ # @@ -41,12 +43,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.m_inmodel = None # Make sure that parfile exists file = self.parfile() @@ -93,7 +95,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")) @@ -121,9 +123,18 @@ class cspull(gammalib.GApplication): """ Get parameters from parfile and setup the observation. """ - # Get parameters + +# sim = ctools.ctobssim() +# sim.run() +# self.obs = sim.obs() +# sim = ctools.ctobssim() +# sim.run() +# self.obs = sim.obs().copy() +# print self.obs +# +# # 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 +155,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) @@ -260,18 +271,23 @@ class cspull(gammalib.GApplication): if self.logExplicit(): self.log.header2("Trial "+str(seed+1)) +# sim = ctools.ctobssim(self.obs) +# sim["seed"].integer(seed) +# sim.run() + print self.obs # Simulate events - obs = obsutils.sim(self.obs, \ - nbins=self.m_enumbins, \ - seed=seed, \ - binsz=self.m_binsz, \ - npix=self.m_npix, \ - edisp=self.m_edisp, \ - log=self.m_log, debug=self.m_debug) - + sim_obs = obsutils.sim(self.obs, \ + nbins=self.m_enumbins, \ + seed=seed, \ + binsz=self.m_binsz, \ + npix=self.m_npix, \ + edisp=self.m_edisp, \ + log=self.m_log, debug=self.m_debug) + print "SIM" + print sim_obs # Determine number of events in simulation nevents = 0.0 - for run in obs: + for run in sim_obs: nevents += run.events().number() # Write simulation results @@ -282,8 +298,10 @@ class cspull(gammalib.GApplication): self.log("\n") # Fit model - like = obsutils.fit(obs, edisp=self.m_edisp, \ - log=self.m_log, debug=self.m_debug) + like = ctools.ctlike(sim_obs) + like.run() +# like = obsutils.fit(obs, edisp=self.m_edisp, \ +# log=self.m_log, debug=self.m_debug) # Store results logL = like.opt().value() @@ -320,7 +338,7 @@ class cspull(gammalib.GApplication): # Compute pull fitted_value = par.value() - real_value = self.model[i][k].value() + real_value = self.obs.models()[i][k].value() error = par.error() if error != 0.0: pull = (fitted_value - real_value) / error 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/cstsdist.py b/scripts/cstsdist.py index 7e8afb8..477d0fe 100755 --- a/scripts/cstsdist.py +++ b/scripts/cstsdist.py @@ -346,6 +346,7 @@ class cstsdist(GApplication): # 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: @@ -381,6 +382,11 @@ class cstsdist(GApplication): self.log(nevents) self.log("\n") +# sim.models(self.full_model) + +# like = obsutils.fit(sim, log = self.m_log, debug = self.m_debug) +# ts = like.obs().models()["Test"].ts() + # Fit background only sim.models(self.bkg_model) like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) @@ -407,6 +413,7 @@ class cstsdist(GApplication): # Fit background and test source sim.models(self.full_model) like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + print like_all.obs().models() result_all = like_all.obs().models() LogL_all = like_all.opt().value() npred_all = like_all.obs().npred() diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index a360bb9..0684e29 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -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_outcube.clear(); m_usepnt = false; - m_xref = 0.0; - m_yref = 0.0; - m_binsz = 0.0; - m_nxpix = 0; - m_nypix = 0; // Initialise protected members m_obs.clear(); @@ -368,16 +357,8 @@ 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_outcube = app.m_outcube; 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; // Copy protected members m_obs = app.m_obs; @@ -405,6 +386,9 @@ void ctbin::free_members(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctbin. + * * Get all task parameters from parameter file or (if required) by querying * the user. Most parameters are only required if no observation exists so * far in the observation container. In this case, a single CTA observation @@ -413,131 +397,60 @@ 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 { - - // Load event list in CTA observation - obs.load(m_evfile); - - // Append CTA observation to container - m_obs.append(obs); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { + std::string msg = "Parameter \"inobs\" is required to be given in ctbin." + "Specify a vaild observation definition (XML or 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_evfile); - - } + // Build observation container without response (not needed) + m_obs = get_observations(false); } // endif: there was no observation in the container // 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(); - // Optionally read ahead parameters so that they get correctly - // dumped into the log file - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - } + // Initialise event cube + GCTAEventCube cube; - // Return - return; -} + // Check if pointing should be used + if (m_usepnt) { + // buld cube from user parameters and pointing information + cube = get_cube(m_obs); -/***********************************************************************//** - * @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(); + } // endif: pointing used as map center - // Set event cube centre, either from the user parameters or from the - // pointing - double xref = m_xref; - double yref = m_yref; - if (m_usepnt) { + else { - // 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; - } - } + // Build cube from user parameters + cube = get_cube(); - // 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); - } + } // endelse: m_usepnt was false - } // endif: used pointing + // Get the skymap + m_cube = cube.map(); - // Create skymap - m_cube = GSkymap(m_proj, m_coordsys, - xref, yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, m_ebounds.size()); + // 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_outcube = (*this)["outcube"].filename(); + } // Return return; } - /***********************************************************************//** * @brief Fill events into counts cube * @@ -590,8 +503,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..8259f04 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..3c99886 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -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 __________________________________________________ */ @@ -407,7 +405,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 +458,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 +481,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; @@ -507,6 +509,8 @@ void ctbkgcube::free_members(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctbkgcube. * @exception GException::invalid_value * No background model definition XML file specified. * @@ -515,11 +519,23 @@ 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 infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { + + std::string msg = "Parameter \"inobs\" is required to be given in ctbkgcube." + "Specify a vaild observation definition (XML or FITS) file to proceed"; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + + // 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 +552,54 @@ 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) == "")) { - // ... otherwise setup the background cube from the counts map - else { - - // Set background cube from counts map - set_from_cntmap(cntmap); - - } + // Check for usepnt flag + bool usepnt = (*this)["usepnt"].boolean(); - // 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); - } + // Check if pointing should be used + if (usepnt) { - // Try first to open as FITS file - try { + // build cube from user parameters and pointing information + m_bkgcube = get_cube(m_obs); - // 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; -} + } // endif: pointing used as map center + else { -/***********************************************************************//** - * @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); + // Build cube from user parameters + m_bkgcube = get_cube(); - // Set background cube from counts map - const GCTAEventCube* cube = dynamic_cast(obs.events()); + } // endelse: m_usepnt was false - // Continue only if cube is valid - if (cube != NULL) { + } // endif: filename was not valid - // 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); - } + // ... otherwise setup the background cube from the counts map + else { + + // Set background cube from counts map + m_bkgcube = set_from_cntmap(incube); + + } // 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 * diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp index 1607e8a..22b6d5e 100644 --- a/src/ctbkgcube/ctbkgcube.hpp +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTBKGCUBE_NAME "ctbkgcube" -#define CTBKGCUBE_VERSION "00-01-00" +#define CTBKGCUBE_VERSION "00-02-00" /***********************************************************************//** @@ -68,11 +68,12 @@ 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..0089ba6 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..d454189 100644 --- a/src/ctbutterfly/ctbutterfly.cpp +++ b/src/ctbutterfly/ctbutterfly.cpp @@ -399,7 +399,7 @@ void ctbutterfly::save(void) void ctbutterfly::init_members(void) { // Initialise members - m_infile.clear(); + //m_inobs.clear(); m_srcname.clear(); m_outfile.clear(); m_ebounds.clear(); @@ -424,7 +424,7 @@ void ctbutterfly::init_members(void) void ctbutterfly::copy_members(const ctbutterfly& app) { // Copy attributes - m_infile = app.m_infile; + // m_inobs = app.m_inobs; m_srcname = app.m_srcname; m_outfile = app.m_outfile; m_ebounds = app.m_ebounds; @@ -455,6 +455,9 @@ void ctbutterfly::free_members(void) * @brief Get application parameters * * @exception GException::invalid_value + * Parameter "inobs" is required for ctbutterfly. + + * @exception GException::invalid_value * Test source not found or no RA/DEC parameters found for test * source. * @@ -466,51 +469,30 @@ 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; - - // 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); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { + std::string msg = "Parameter \"inobs\" is required to be given in ctbutterfly." + "Specify a vaild observation definition (XML or 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(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)); diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp index 329175f..4c83796 100644 --- a/src/ctbutterfly/ctbutterfly.hpp +++ b/src/ctbutterfly/ctbutterfly.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTBUTTERFLY_NAME "ctbutterfly" -#define CTBUTTERFLY_VERSION "00-01-00" +#define CTBUTTERFLY_VERSION "00-02-00" /***********************************************************************//** @@ -80,7 +80,6 @@ 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 diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index 7cc0c48..0fea881 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..35e6d65 100644 --- a/src/ctcubemask/ctcubemask.cpp +++ b/src/ctcubemask/ctcubemask.cpp @@ -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,42 +349,23 @@ 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 (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); + // If there are no observations in container then load them via user parameters + if (m_obs.size() == 0) { - // Append CTA observation to container - m_obs.append(obs); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { - // 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 - m_obs.load(m_infile); - - // Signal that XML file should be used for storage - m_use_xml = true; + std::string msg = "Parameter \"inobs\" is required to be given in ctcubemask." + "Specify a vaild observation definition (XML or FITS) file to proceed" + "It is not possible to build a binned observation from scratch"; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } - } + // Build observation container without response (not needed) + m_obs = get_observations(false); - } // endif: there was no observation in the container + } // endif: there was no observation in the container // Get parameters m_regfile = (*this)["regfile"].filename(); @@ -396,7 +381,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 +523,9 @@ void ctcubemask::apply_mask(GCTAObservation* obs) void ctcubemask::init_members(void) { // Initialise parameters - m_infile.clear(); + //m_infile.clear(); m_regfile.clear(); - m_outfile.clear(); + m_outcube.clear(); m_prefix.clear(); m_usepnt = false; m_ra = 0.0; @@ -552,7 +537,7 @@ void ctcubemask::init_members(void) // Initialise protected members m_obs.clear(); m_infiles.clear(); - m_use_xml = false; +// m_use_xml = false; // Return return; @@ -567,9 +552,9 @@ void ctcubemask::init_members(void) void ctcubemask::copy_members(const ctcubemask& app) { // Copy parameters - m_infile = app.m_infile; + //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 +566,7 @@ 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; +// m_use_xml = app.m_use_xml; // Return return; @@ -622,7 +607,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 +646,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 +673,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 +711,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..034fd57 100644 --- a/src/ctcubemask/ctcubemask.hpp +++ b/src/ctcubemask/ctcubemask.hpp @@ -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,7 @@ 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 + //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..d68308e 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..ea01fd3 100644 --- a/src/ctexpcube/ctexpcube.cpp +++ b/src/ctexpcube/ctexpcube.cpp @@ -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,13 @@ 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(); + //m_ebounds.clear(); // Return return; @@ -327,13 +327,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; // Return return; @@ -352,6 +351,8 @@ void ctexpcube::free_members(void) /***********************************************************************//** * @brief Get application parameters + * @exception GException::invalid_value + * Parameter "inobs" is required for ctexpcube. * * Get all task parameters from parameter file or (if required) by querying * the user. The parameters are read in the correct order. @@ -360,159 +361,72 @@ void ctexpcube::free_members(void) ***************************************************************************/ 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); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { + + std::string msg = "Parameter \"inobs\" is required to be given in ctexpcube." + "Specify a vaild observation definition (XML or FITS) file to proceed"; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + + // Build observation container + m_obs = get_observations(); + + } // endif: there was no observation in the container // 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(); + // Get the incube filename + std::string incube = (*this)["incube"].filename(); - // Define exposure cube - m_expcube = GCTAExposure(wcs, coordsys, xref, yref, - -binsz, binsz, nxpix, nypix, - m_ebounds); - } - - // ... otherwise setup the exposure cube from the counts map - else { + // Check for filename validity + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { - // Set exposure cube from counts map - set_from_cntmap(cntmap); - - } + // Initialise event cube + GCTAEventCube cube; - // Read output filename (if needed) - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - } + // Check for usepnt flag + bool usepnt = (*this)["usepnt"].boolean(); - // Return - return; -} + // Check if pointing should be used + if (usepnt) { + // build cube from user parameters and pointing information + cube = get_cube(m_obs); -/***********************************************************************//** - * @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 { + } // endif: pointing used as map center - // Allocate CTA observation - GCTAObservation obs; - - // Load input file in CTA observation - obs.load(filename); + else { - // Append CTA observation to container - m_obs.append(obs); + // Build cube from user parameters + cube = get_cube(); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { + } // endelse: m_usepnt was false - // 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 { + // Set exposure cube from counts map + GCTAEventCube cube = set_from_cntmap(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..546ed6e 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..8ddef07 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -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 @@ -341,6 +341,9 @@ void ctlike::save(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctlike. + * * Get all required task parameters from the parameter file or (if specified) * by querying the user. Observation dependent parameters will only be read * if the observation container is actually empty. Observation dependent @@ -357,56 +360,38 @@ 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; - - // Get event file name - std::string filename = (*this)["infile"].filename(); - - // 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); - + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { + std::string msg = "Parameter \"inobs\" is required to be given in ctlike." + "Specify a vaild observation definition (XML or 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) { + // Build observation container + m_obs = get_observations(); - // 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 +405,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 +532,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 +558,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..ad8b852 100644 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -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..7ac742f 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" +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" -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" -outmdl, f, a, "crab_results.xml",,, "Source model output file" +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..c02316c 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -207,9 +207,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 +243,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 +310,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 +357,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 +379,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 +411,83 @@ 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) { - // Read model definiton file if required - if (m_obs.models().size() == 0) { - m_srcmdl = (*this)["srcmdl"].filename(); - } + // Build observation container + m_obs = get_observations(); - // Get energy dispersion flag parameters - m_apply_edisp = (*this)["edisp"].boolean(); + } // endif: there was no observation in the container - // 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(); - } + // Read model definition file if required + if (m_obs.models().size() == 0) { + std::string inmodel = (*this)["inmodel"].filename(); - } + // Load models from file + m_obs.models(inmodel); - // 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; - } - } - } + } // endif: there were no models - // 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(); + // Read optionally output cube filenames + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); } - // 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 - else { - - // First try to open the file as an event list or counts map - try { + // Get energy dispersion flag parameters + m_apply_edisp = (*this)["edisp"].boolean(); - // Allocate CTA observation - GCTAObservation obs; + if (!m_has_cube) { - // Load CTA observation from FITS file - obs.load(m_obsfile); + // Read cube definition file + std::string incube = (*this)["incube"].filename(); - // Append CTA observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { + // If no cube file has been specified then read all parameters that + // are necessary to create the cube from scratch (see method ctool::get_cube()) + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { - // Load observations from XML file. This will throw - // an exception if it fails. - m_obs.load(m_obsfile); + // Get parameters + bool usepnt = (*this)["usepnt"].boolean(); - } + // Initialise event cube + GCTAEventCube cube; - } // endelse: loaded information from input file + // Check if pointing should be used + if (usepnt) { - } // endif: there was no observation in the container + // buld cube from user parameters and pointing information + m_cube = get_cube(m_obs); - // 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)); - } + } // endif: pointing used as map center - // For all observations that have no response, set the response - // from the task parameters - set_response(m_obs); + else { - // Return - return; -} + // Build cube from user parameters + m_cube = get_cube(); + } // endelse: m_usepnt was false -/***********************************************************************//** - * @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); - - } + } // endif: incube file was not vaild - // ... otherwise load cube from file - else { + // ... else: build map information from file + else { - // Load cube - m_cube.load(m_infile); + // 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); + // 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..8e2f072 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..e07e21b 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -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 @@ -427,7 +427,6 @@ void ctobssim::run(void) return; } - /***********************************************************************//** * @brief Save the selected event list(s) * @@ -490,35 +489,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_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 +519,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,48 +558,22 @@ 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 (m_obs.size() == 0) { - - // 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); + // If there are no observations in container then load them via user parameters + if (m_obs.size() == 0) { - // Append CTA observation to container - m_obs.append(obs); + // Build observation container + m_obs = get_observations(); - // Load models into container - m_obs.models(m_infile); + } // endif: there was no observation in the container - // Signal that no XML file should be used for storage - m_use_xml = false; + // Read model definition file if required + if (m_obs.models().size() == 0) { + std::string inmodel = (*this)["inmodel"].filename(); - } // endif: there was no observation in the container + // Load models from file + m_obs.models(inmodel); - // ... otherwise we signal that an XML file should be written - else { - m_use_xml = true; - } + } // endif: there were no models // Get other parameters m_seed = (*this)["seed"].integer(); @@ -635,20 +582,10 @@ 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_outevents = (*this)["outevents"].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); - } - } - } - // Initialise random number generators. We initialise here one random // number generator per observation so that each observation will // get it's own random number generator. This will lead to identical @@ -700,81 +637,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 +1025,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 +1053,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 +1101,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..695f922 100644 --- a/src/ctobssim/ctobssim.hpp +++ b/src/ctobssim/ctobssim.hpp @@ -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..201e142 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,8 +17,9 @@ # # File information #================= -infile, f, a, "$CTOOLS/share/models/crab.xml",,, "Model" -outfile,f, a, "events.fits",,, "Output event data file or 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..c690460 100644 --- a/src/ctpsfcube/ctpsfcube.cpp +++ b/src/ctpsfcube/ctpsfcube.cpp @@ -275,10 +275,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 +297,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 +317,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; @@ -344,6 +342,9 @@ void ctpsfcube::free_members(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctbin. + * * Get all task parameters from parameter file or (if required) by querying * the user. The parameters are read in the correct order. * @@ -351,162 +352,75 @@ void ctpsfcube::free_members(void) ***************************************************************************/ 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); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { - // Read energy dispersion flag - m_apply_edisp = (*this)["edisp"].boolean(); + std::string msg = "Parameter \"inobs\" is required to be given in ctpsfcube." + "Specify a vaild observation definition (XML or FITS) file to proceed"; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } - // 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); - - } + // Build observation container + m_obs = get_observations(); - // Read output filename (if needed) - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - } + } // endif: there was no observation in the container - // Return - return; -} + // Get the incube filename + std::string incube = (*this)["incube"].filename(); + // Get additional binning parameters + double dmax = (*this)["amax"].real(); + int ndbins = (*this)["anumbins"].integer(); -/***********************************************************************//** - * @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(); + // Check for filename validity + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { - // Try first to open as FITS file - try { + // Get parameters + bool usepnt = (*this)["usepnt"].boolean(); - // Allocate CTA observation - GCTAObservation obs; - - // Load input file in CTA observation - obs.load(filename); + // Initialise event cube + GCTAEventCube cube; - // Append CTA observation to container - m_obs.append(obs); + // Check if pointing should be used + if (usepnt) { - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { + // buld cube from user parameters and pointing information + cube = get_cube(m_obs); - // Load observations from XML file - m_obs.load(filename); + } // endif: pointing used as map center + + else { + // Build cube from user parameters + cube = get_cube(); + + } // endelse: m_usepnt was false + + // Define psf cube + m_psfcube = GCTAMeanPsf(cube, dmax, ndbins); } - // Return - return; -} + // ... otherwise setup the exposure cube from the counts map + else { + // Set exposure cube from counts map + GCTAEventCube cube = set_from_cntmap(incube); -/***********************************************************************//** - * @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); - } + // Define psf cube + m_psfcube = GCTAMeanPsf(cube, dmax, ndbins); - } // 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 energy dispersion flag + m_apply_edisp = (*this)["edisp"].boolean(); + + // 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..0583757 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 exposure cube definition" +outcube, f, a, "psfcube.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/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 4dd0400..b66cb32 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -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,6 +456,9 @@ void ctselect::free_members(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctselect. + * * Get all task parameters from parameter file or (if required) by querying * the user. Times are assumed to be in the native CTA MJD format. * @@ -475,40 +468,20 @@ 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); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { - // Signal that no XML file should be used for storage - m_use_xml = false; - + std::string msg = "Parameter \"inobs\" is required to be given in ctselect." + "Specify a vaild observation definition (XML or 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 @@ -537,7 +510,7 @@ 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_outobs = (*this)["outobs"].filename(); m_prefix = (*this)["prefix"].string(); } @@ -1048,7 +1021,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 +1033,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 +1076,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 +1102,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 +1140,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..9a27231 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -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..e1642eb 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..5fac288 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -35,7 +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 __________________________________________________ */ /* __ Coding definitions _________________________________________________ */ @@ -262,12 +262,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 +305,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; @@ -324,6 +318,9 @@ void ctskymap::save(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for ctskymap. + * * Get all task parameters from parameter file or (if required) by querying * the user. Most parameters are only required if no observation exists so * far in the observation container. In this case, a single CTA observation @@ -332,93 +329,54 @@ 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 (m_obs.size() == 0) { + // 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 infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { - // Load unbinned CTA observation - GCTAObservation obs; - obs.load(m_evfile); + std::string msg = "Parameter \"inobs\" is required to be given in ctskymap." + "Specify a vaild observation definition (XML or FITS) file to proceed"; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } - // 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(); + // Build observation container without response (not needed) + m_obs = get_observations(false); - } // endif: there was no observation in the container + } // endif: there was no observation in the container - // 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(); + // Get parameters + bool usepnt = (*this)["usepnt"].boolean(); - // Read ahead parameters - if (read_ahead()) { - m_outfile = (*this)["outfile"].filename(); - } + // Check if pointing should be used + if (usepnt) { - // Return - return; -} + // buld map from user parameters and pointing information + m_skymap = get_map(m_obs); + } // endif: pointing used as map center -/***********************************************************************//** - * @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 { + // Build cube from user parameters + m_skymap = get_map(); - // 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: m_usepnt was false - } // endelse: map centre set to pointing + // Get remaining parameters + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); - // Create skymap - m_skymap = GSkymap(m_proj, m_coordsys, xref, yref, -m_binsz, m_binsz, - m_nxpix, m_nypix, 1); + // Read ahead parameters + if (read_ahead()) { + m_outmap = (*this)["outmap"].filename(); + } // Return return; } - /***********************************************************************//** * @brief Map events into a sky map * @@ -473,8 +431,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 +483,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; // Return return; @@ -552,19 +501,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; // Return return; diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index de00c4b..b5ebb97 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -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..a0416b8 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -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,15 +436,8 @@ 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; @@ -486,6 +470,9 @@ void cttsmap::free_members(void) /***********************************************************************//** * @brief Get application parameters * + * @exception GException::invalid_value + * Parameter "inobs" is required for cttsmap. + * * @exception GException::invalid_value * Test source not found or no RA/DEC parameters found for test * source. @@ -498,42 +485,20 @@ 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; - - // 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); + // Throw exception if no infile is given, since this tool needs an observation + // including events + if ((*this)["inobs"].filename()=="NONE" || (*this)["inobs"].filename() == "") { + std::string msg = "Parameter \"inobs\" is required to be given in cttsmap." + "Specify a vaild observation definition (XML or 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(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 @@ -542,7 +507,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 +544,27 @@ 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(); + + // Get sky map binning from user parameters + // and check if pointing should be used + GSkymap map; + bool usepnt = (*this)["usepnt"].boolean(); + if (usepnt) { + + // build cube from user parameters and pointing information + map = get_map(m_obs); + + } // endif: pointing used as map center + + else { + + // Build cube from user parameters + map = get_map(); + + } // endelse: m_usepnt was false + + // Initialise maps from user parameters + init_maps(map); // Get optional splitting parameters m_binmin = (*this)["binmin"].integer(); @@ -596,7 +574,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 +587,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 +614,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..f17bfdc 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTTSMAP_NAME "cttsmap" -#define CTTSMAP_VERSION "00-01-00" +#define CTTSMAP_VERSION "00-02-00" /***********************************************************************//** @@ -81,19 +81,11 @@ 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 diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 09c5628..d0db7b7 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" +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" -outfile, f, a, "tsmap.fits",,, "Output counts cube" +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..3c83843 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -48,6 +48,8 @@ /* __ Method name definitions ____________________________________________ */ #define G_GET_EBOUNDS "ctool::get_ebounds()" +#define G_SET_FROM_CNTMAP "ctool::set_from_cntmap(std::string&)" +#define G_GET_CUBE "ctool::get_cube(GObservations&)" /* __ Debug definitions __________________________________________________ */ @@ -226,6 +228,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 +252,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; @@ -342,6 +352,412 @@ GEbounds ctool::get_ebounds(void) } + + +/***********************************************************************//** + * @brief Get a GCTAEventCube from file + * + * @param[in] filename FITS file of count cube. + * @exception GException::invalid_value + * Counts map not of WCS type + * + * Gets an cta event cube from file. This function is used in many ctools + * to obtain a skymap definition from file + * + ***************************************************************************/ +GCTAEventCube ctool::set_from_cntmap(const std::string filename) +{ + // Allocate CTA observation + GCTAEventCube cube; + + // Load counts map in CTA observation + cube.load(filename); + + // Check for valid WCS coordinate system + const GWcs* wcs = dynamic_cast(cube.map().projection()); + if (wcs == NULL) { + + std::string msg = "Counts map project is not of WCS type."; + throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); + + } // endif: WCS projection was not valid + + // Return cube + return cube; +} + + + +/***********************************************************************//** + * @brief Get observation container + * + * @param[in] get_response flag whether response information should be loaded + * + * Get an observation container according to the user parameters. The method + * supports loading of a individual FITS file, 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 + // (see method setup_obs) + if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) { + + // Setup a new CTA observation + GCTAObservation cta_obs = setup_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 "" + + 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 Build a CTA observation with an empty event list + * + * 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. + * + ***************************************************************************/ +GCTAObservation ctool::setup_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; + +} + +GSkyDir ctool::get_pointing(const GObservations& obs) +{ + // Initialise values + double xref = 0.0; + double yref = 0.0; + int n_pnt = 0; + bool found = false; + + // Loop over all observations to compute mean pointings + for (int i = 0; i < obs.size(); ++i) { + 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; + + // Signal that pointing has been found + found = true; + + } // endif: cta observation was valid + + } // endfor: loop over alle observations + + // 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"; + throw GException::invalid_value(G_GET_CUBE, 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 Build a GCTAEventCube based on user parameters + * + * Creates an empty CTA event cube from input user parameters. It appends a dummy + * GTi to the event cube. + * + ***************************************************************************/ +GCTAEventCube ctool::get_cube(void) +{ + // Get cube parameters + double xref = (*this)["xref"].real(); + double 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(); + + // Set energy boundaries + GEbounds ebounds = get_ebounds(); + + // Set dummy GTI that is needed for event cube creation + GGti gti; + gti.append(GTime(0.0), GTime(0.1234)); + + // Initialise sky map + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, ebounds.size()); + + // Build event cube from given information + GCTAEventCube cube = GCTAEventCube(map,ebounds,gti); + + // Return cube + return cube; + +} + +/***********************************************************************//** + * @brief Build a GCTAEventCube based on user parameters + * + * @param[in] obs observation container + * + * Creates an empty CTA event cube from input user parameters. It appends a dummy + * GTi to the event cube. The input parameter is used to extract the pointing information + * from the given observations. This function is used if a ctool has the input parameter + * usepnt=true. + * + ***************************************************************************/ +GCTAEventCube ctool::get_cube(const GObservations& obs) +{ + // Read ahead map information + std::string proj = (*this)["proj"].string(); + std::string coordsys = (*this)["coordsys"].string(); + + // Initialise reference coordinates + double xref = 0.0; + double yref = 0.0; + + // Get pointing from observations + GSkyDir pnt = get_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(); + } + + // Read further binning parameters + double binsz = (*this)["binsz"].real(); + int nxpix = (*this)["nxpix"].integer(); + int nypix = (*this)["nypix"].integer(); + + // Set energy boundaries + GEbounds ebounds = get_ebounds(); + + // Set dummy GTI that is needed for event cube creation + GGti gti; + gti.append(GTime(0.0), GTime(0.1234)); + + // Build skymap from user parameters + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, ebounds.size()); + + // Initialise event cube from all parameters + GCTAEventCube cube = GCTAEventCube(map,ebounds,gti); + + // Return event cube + return cube; + +} + +/***********************************************************************//** + * @brief Build a skymap based on user parameters + * + * Creates an empty skymap from input user parameters. + * + ***************************************************************************/ +GSkymap ctool::get_map() +{ + // Read map information + double xref = (*this)["xref"].real(); + double 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(); + + // Initialise sky map + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, 1); + + // Return sky map + return map; +} + +/***********************************************************************//** + * @brief Build a skymap based on user parameters + * + * Creates an empty skymap from input user parameters. + * + ***************************************************************************/ +GSkymap ctool::get_map(const GObservations& obs) +{ + // Read ahead map information + std::string proj = (*this)["proj"].string(); + std::string coordsys = (*this)["coordsys"].string(); + + // Initialise reference coordinates + double xref = 0.0; + double yref = 0.0; + + // Get pointing from observations + GSkyDir pnt = get_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(); + } + + double binsz = (*this)["binsz"].real(); + int nxpix = (*this)["nxpix"].integer(); + int nypix = (*this)["nypix"].integer(); + + // Initialise sky map + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix, 1); + + // Return sky map + return map; +} + /***********************************************************************//** * @brief Set response for all CTA observations in container * diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 39acdec..9035b18 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -67,12 +67,23 @@ protected: void free_members(void); const bool& read_ahead(void) const; GEbounds get_ebounds(void); + GObservations get_observations(const bool& get_response=true); + GCTAObservation setup_obs(void); + GSkymap get_map(void); + GCTAEventCube get_cube(void); + GSkymap get_map(const GObservations& obs); + GCTAEventCube get_cube(const GObservations& obs); + GCTAEventCube set_from_cntmap(const std::string filename); + GSkyDir get_pointing(const GObservations& obs); + void set_response(GObservations& obs); void set_obs_response(GCTAObservation* obs); size_t get_current_rss(void); // Protected members 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") -- libgit2 0.21.4