diff --git a/.gitignore b/.gitignore index e4ac738..6971cc5 100644 --- a/.gitignore +++ b/.gitignore @@ -54,6 +54,7 @@ src/ctselect/ctselect src/ctskymap/ctskymap src/ctexpcube/ctexpcube src/ctpsfcube/ctpsfcube +src/ctbkgcube/ctbkgcube src/ctcubemask/ctcubemask test/pfiles/ test/*.fits diff --git a/ChangeLog b/ChangeLog index 8595f65..ec35376 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-07-20 Juergen Knoedlseder +2014-07-22 Juergen Knoedlseder * Version 0.8.0 released ======================== @@ -6,7 +6,8 @@ Adapt to GammaLib 0.9.0 interface ctbin now always merges event lists into a single counts cube Add ctexpcube tool (generates exposure cube) - Add ctpsfcube tool (generated PSF cube) + Add ctpsfcube tool (generates PSF cube) + Add ctbkgcube tool (generates background cube) Add ctcubemask tool (masks event bins in event cube) Support energy dispersion in ctobssim, ctlike and ctmodel Correct ctobssim for thread safety diff --git a/NEWS b/NEWS index 4da0e53..8c66f8a 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 Author(s) : Juergen Knoedlseder -20 July 2014 +22 July 2014 1. Introduction @@ -97,7 +97,15 @@ definition file. The cube definition may be taken from a counts cube or may be specified by user parameters. -12. ctcubemask - CTA cube filtering +12. ctbkgcube - CTA background cube generation +---------------------------------------------- +The first version of the ctbkgcube tool is included in the release. +ctbkgcube computes the background cube for an observation definition +file. The cube definition may be taken from a counts cube or may be +specified by user parameters. + + +13. ctcubemask - CTA cube filtering ----------------------------------- The first version of the ctcubemask tool is included in the release. ctcubemask allows filtering of bins in a counts cube by setting the @@ -105,12 +113,12 @@ bin content to -1. Filtering can be done by energy, for a circular Region of Interest (ROI) or using ds9 exclusion regions. -13. Examples +14. Examples ------------ None -14. Scripts +15. Scripts ----------- Add logL and event statistics to cspull output file. @@ -119,11 +127,11 @@ Convert all tabs into spaces. The obsutils script has now been moved into the ctools module namespace. -15. caldb +16. caldb --------- None -16. models +17. models ---------- None diff --git a/configure.ac b/configure.ac index 706bf89..187975d 100644 --- a/configure.ac +++ b/configure.ac @@ -884,6 +884,7 @@ AC_CONFIG_FILES([Makefile src/ctselect/Makefile src/ctexpcube/Makefile src/ctpsfcube/Makefile + src/ctbkgcube/Makefile src/ctcubemask/Makefile src/ctools-setup test/Makefile diff --git a/pyext/ctbkgcube.i b/pyext/ctbkgcube.i new file mode 100644 index 0000000..d93b533 --- /dev/null +++ b/pyext/ctbkgcube.i @@ -0,0 +1,68 @@ +/*************************************************************************** + * ctbkgcube - CTA background cube tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014 by Chia-Chun Lu * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctbkgcube.i + * @brief CTA background cube tool definition + * @author Chia-Chun Lu + */ + +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctbkgcube.hpp" +#include "GTools.hpp" +%} + +/***********************************************************************//** + * @class ctbkgcube + * + * @brief CTA background cube tool + ***************************************************************************/ +class ctbkgcube : public GApplication { + +public: + // Constructors and destructors + ctbkgcube(void); + explicit ctbkgcube(const GObservations& obs); + ctbkgcube(int argc, char *argv[]); + ctbkgcube(const ctbkgcube& app); + virtual ~ctbkgcube(void); + + // Methods + void clear(void); + void execute(void); + void run(void); + void save(void); + const GObservations& obs(void) const; + void get_parameters(void); + void get_obs(void); + void get_ebounds(void); + void set_from_cntmap(const std::string& filename); + void fill_cube(GCTAObservation* obs); +}; + +/***********************************************************************//** + * @brief CTA background cube tool Python extension + ***************************************************************************/ +%extend ctbkgcube { + ctbkgcube copy() { + return (*self); + } +} diff --git a/pyext/ctools/ctools.i b/pyext/ctools/ctools.i index a51c413..956fa6d 100755 --- a/pyext/ctools/ctools.i +++ b/pyext/ctools/ctools.i @@ -53,4 +53,5 @@ %include "ctskymap.i" %include "ctexpcube.i" %include "ctpsfcube.i" +%include "ctbkgcube.i" %include "ctcubemask.i" diff --git a/src/Makefile.am b/src/Makefile.am index 75dcfd4..60dfc1e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -13,6 +13,7 @@ SUBDIRS = ctbin \ ctskymap \ ctexpcube \ ctpsfcube \ + ctbkgcube \ ctcubemask # Scripts to be distributed @@ -33,6 +34,7 @@ libctools_la_LIBADD = ctbin/libctbin.la \ ctskymap/libctskymap.la \ ctexpcube/libctexpcube.la \ ctpsfcube/libctpsfcube.la \ + ctbkgcube/libctbkgcube.la \ ctcubemask/libctcubemask.la # Distribute scripts diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index 2cd91fd..cf4abdf 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -777,7 +777,7 @@ void ctbin::fill_cube(GCTAObservation* obs) // Make sure that the observation holds a CTA event list. If this // is not the case then throw an exception. if (dynamic_cast(obs->events()) == NULL) { - throw GException::no_list(G_BIN_EVENTS); + throw GException::no_list(G_FILL_CUBE); } // Initialise binning statistics diff --git a/src/ctbkgcube/Makefile.am b/src/ctbkgcube/Makefile.am new file mode 100644 index 0000000..fa55a8a --- /dev/null +++ b/src/ctbkgcube/Makefile.am @@ -0,0 +1,21 @@ +## Process this file with automake to produce Makefile.in + +# Program name +bin_PROGRAMS = ctbkgcube + +# Sources +sources = ctbkgcube.cpp + +# Libtool library +noinst_LTLIBRARIES = libctbkgcube.la +libctbkgcube_la_SOURCES = $(sources) + +# Program definition +ctbkgcube_SOURCES = main.cpp +ctbkgcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctbkgcube_CPPFLAGS = @CPPFLAGS@ +ctbkgcube_LDADD = libctbkgcube.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctbkgcube.par diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp new file mode 100644 index 0000000..50ce350 --- /dev/null +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -0,0 +1,795 @@ +/*************************************************************************** + * ctbkgcube - CTA background cube tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014 by Chia-Chun Lu * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctbkgcube.cpp + * @brief CTA background cube tool implementation + * @author Chia-Chun Lu + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctbkgcube.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_GET_PARAMETERS "ctbkgcube::get_parameters()" +#define G_GET_OBS "ctbkgcube::get_obs()" +#define G_GET_EBOUNDS "ctbkgcube::get_ebounds()" +#define G_SET_FROM_CNTMAP "ctbkgcube::set_from_cntmap(std::string&)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctbkgcube::ctbkgcube(void) : GApplication(CTBKGCUBE_NAME, CTBKGCUBE_VERSION) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * param[in] obs Observation container. + * + * This method creates an instance of the class by copying an existing + * observations container. + ***************************************************************************/ +ctbkgcube::ctbkgcube(const GObservations& obs) : + GApplication(CTBKGCUBE_NAME, CTBKGCUBE_VERSION) +{ + // Initialise members + init_members(); + + // Set observations + m_obs = obs; + + // Write header into logger + log_header(); + + // Return + return; +} + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctbkgcube::ctbkgcube(int argc, char *argv[]) : + GApplication(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app ctbkgcube application. + ***************************************************************************/ +ctbkgcube::ctbkgcube(const ctbkgcube& app) : GApplication(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctbkgcube::~ctbkgcube(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app ctbkgcube application. + * @return Returns ctbkgcube application. + ***************************************************************************/ +ctbkgcube& ctbkgcube::operator=(const ctbkgcube& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->GApplication::operator=(app); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear instance + ***************************************************************************/ +void ctbkgcube::clear(void) +{ + // Free members + free_members(); + this->GApplication::free_members(); + + // Initialise members + this->GApplication::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Execute application + * + * This is the main execution method of the ctbkgcube class. It is invoked + * when the executable is called from command line. The method generates + * the background cube and saves the result. + ***************************************************************************/ +void ctbkgcube::execute(void) +{ + // Signal that some parameters should be read ahead + m_read_ahead = true; + + // Create the background cube + run(); + + // Save the background cube into the output FITS file + save(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate the background cube(s). + * + * This method reads the task parameters from the parfile, sets up the + * observation container, loops over all CTA observations in the container + * and generates a background cube from the CTA observations. + ***************************************************************************/ +void ctbkgcube::run(void) +{ + // If we're in debug mode then all output is also dumped on the screen + if (logDebug()) { + log.cout(true); + } + + // Get task parameters + get_parameters(); + + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Observations"); + } + else { + log.header1("Observation"); + } + log << m_obs << std::endl; + } + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Prepare model"); + } + + // Copy models from observation container + m_bkgmdl = m_obs.models(); + + // Remove all model that are not CTA background models from the + // container + int num = m_bkgmdl.size(); + for (int i = num-1; i >= 0; --i) { + + // Flag removal + bool remove = true; + + // Do we have a CTA, HESS or VERITAS specific model? + if (m_bkgmdl[i]->is_valid("CTA", "") || + m_bkgmdl[i]->is_valid("HESS", "") || + m_bkgmdl[i]->is_valid("VERITAS", "")) { + + // Do we have a background model? + if (dynamic_cast(m_bkgmdl[i]) != NULL) { + + // ... then keep model + remove = false; + + } // endif: had a background model + } // endif: had a CTA, HESS or VERITAS model + + // Log results + if (logTerse()) { + if (remove) { + log << gammalib::parformat("Remove model"); + } + else { + log << gammalib::parformat("Keep model"); + } + log << m_bkgmdl[i]->name(); + log << " "; + log << m_bkgmdl[i]->type(); + log << "("; + log << m_bkgmdl[i]->instruments(); + log << ")"; + log << std::endl; + } + + // Remove model if requested + if (remove) { + m_bkgmdl.remove(i); + } + + } // endfor: looped over all background models + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Generate background cube"); + } + + // Initialise observation counter + int n_observations = 0; + + // Loop over all observations in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Continue only if observation is a CTA observation + if (obs != NULL) { + + // Write header for observation + if (logTerse()) { + if (obs->name().length() > 1) { + log.header3("Observation "+obs->name()); + } + else { + log.header3("Observation "+gammalib::str(i)); + } + } + + // Increment number of observations + n_observations++; + + // Fill the cube + fill_cube(obs); + + } // endif: CTA observation found + + } // endfor: looped over observations + + // Log results + if (logTerse()) { + log << gammalib::parformat("Number of observations"); + log << n_observations << std::endl; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate background cube + * + * @param[in] obs Pointer to CTA observation. + * + * Fills + ***************************************************************************/ +void ctbkgcube::fill_cube(GCTAObservation* obs) +{ + // Continue only if observation pointer is valid + if (obs != NULL) { + + // Initialise statistics + double sum = 0.0; + + // Set GTI of actual observations as the GTI of the event cube + m_bkgcube.gti(obs->events()->gti()); + + // Loop over all bins in background cube + for (int i = 0; i < m_bkgcube.size(); ++i) { + + // Get event bin + GCTAEventBin* bin = m_bkgcube[i]; + + // Compute model value for event bin + double model = + m_bkgmdl.eval(*(const_cast(bin)), *obs) * + bin->size(); + + // Add existing number of counts + model += bin->counts(); + + // Store cumulated value + bin->counts(model); + + // Sum all events + sum += model; + } + + // Log results + if (logTerse()) { + log << gammalib::parformat("Background events in cube"); + log << sum << std::endl; + } + + } // endif: observation pointer was not valid + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save background cube + * + * Save the background cube into the file specified by the outfile parameter. + ***************************************************************************/ +void ctbkgcube::save(void) +{ + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Save background cube"); + } + + // Create energies container from energy boundaries + GEnergies energies; + for (int i = 0; i < m_ebounds.size(); ++i) { + energies.append(m_ebounds.elogmean(i)); + } + + // Get output filename + std::string filename = (*this)["outfile"].filename(); + + // Create empty FITS file + GFits fits; + + // Write background cube + m_bkgcube.map().write(fits); + + // Write energies + energies.write(fits); + + // Save FITS file + fits.saveto(filename, clobber()); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * @exception GException::invalid_value + * No background model definition XML file specified. + * + * Get all task parameters from parameter file or (if required) by querying + * the user. The parameters are read in the correct order. + ***************************************************************************/ +void ctbkgcube::get_parameters(void) +{ + // If we do not have any observations in the container then load them + // using the "infile" parameter + if (m_obs.size() == 0) { + get_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) { + if ((*this)["bkgmdl"].is_undefined()) { + std::string msg = "No background model definition XML file " + "specified. Please set the \"bkgmdl\" parameter " + "to the background model definition XML " + "filename."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + std::string bkgmdl = (*this)["bkgmdl"].filename(); + GModels models(bkgmdl); + 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 + 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); + + } + + // ... otherwise setup the background cube from the counts map + else { + + // Set background cube from counts map + set_from_cntmap(cntmap); + + } + + // Read output filename (if needed) + if (m_read_ahead) { + std::string filename = (*this)["outfile"].filename(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get observation definition + * + * @exception GException::invalid_value + * No event list, counts cube or observation definition file + * specified. + * + * Get observation definition from the user parameters. + ***************************************************************************/ +void ctbkgcube::get_obs(void) +{ + // Get input filename + std::string filename = (*this)["infile"].filename(); + + // Check whether infile parameter is undefined + if (gammalib::toupper(filename) == "NONE") { + std::string msg = "No event list, counts cube or observation definition " + "file specified. Please set the \"infile\" parameter " + "to either an event list filename, a counts cube " + "filename or an observation definition filename."; + throw GException::invalid_value(G_GET_OBS, msg); + } + + // Try first to open as FITS file + try { + + // Allocate CTA observation + GCTAObservation obs; + + // Load input file in CTA observation + obs.load(filename); + + // Append CTA observation to container + m_obs.append(obs); + + } + + // ... otherwise try to open as XML file + catch (GException::fits_open_error &e) { + + // Load observations from XML file + m_obs.load(filename); + + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get energy boundaries from parameters + * + * @exception GException::invalid_value + * Energy bin definition file not defined or invalid extension + * name encountered. + * + * Get the energy boundaries from the user parameters. + ***************************************************************************/ +void ctbkgcube::get_ebounds(void) +{ + // Determine the energy binning alogrithm + std::string ebinalg = (*this)["ebinalg"].string(); + + // If we have the binning given by a file then try to get the boundaries + // from that file + if (ebinalg == "FILE") { + + // Get filename + std::string filename = (*this)["ebinfile"].filename(); + + // Check whether infile parameter is undefined + if (gammalib::toupper(filename) == "NONE") { + std::string msg = "No energy bin definition file specified. " + "Please set the \"ebinfile\" parameter to a " + "file containing the energy bin definition."; + throw GException::invalid_value(G_GET_EBOUNDS, msg); + } + + // Open fits file to check which extension is given + GFits file(filename); + + // Check first for EBOUNDS extension + if (file.contains("EBOUNDS")) { + file.close(); + m_ebounds.load(filename, "EBOUNDS"); + } + + // ... then check for ENERGYBINS extension + else if (file.contains("ENERGYBINS")) { + file.close(); + m_ebounds.load(filename, "ENERGYBINS"); + } + + // ... otherwise throw an exception + else { + file.close(); + std::string msg = "No extension with name \"EBOUNDS\" or" + " \"ENERGYBINS\" found in FITS file" + " \""+filename+"\".\n" + "An \"EBOUNDS\" or \"ENERGYBINS\" extension" + " is required if the parameter \"ebinalg\"" + " is set to \"FILE\"."; + throw GException::invalid_value(G_GET_EBOUNDS, msg); + } + } + + // ... otherwise read emin, emax and nebins + else { + + // Get the relevant parameters + double emin = (*this)["emin"].real(); + double emax = (*this)["emax"].real(); + int enumbins = (*this)["enumbins"].integer(); + bool log = ((*this)["ebinalg"].string() == "LIN") ? false : true; + + // Create energy boundaries + m_ebounds = GEbounds(enumbins, + GEnergy(emin, "TeV"), + GEnergy(emax, "TeV"), + log); + + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set background cube definition from counts map + * + * @exception GException::invalid_value + * Invalid counts map projection or invalid events encountered. + * + * Set background cube definition from counts map. + ***************************************************************************/ +void ctbkgcube::set_from_cntmap(const std::string& filename) +{ + // Allocate CTA observation + GCTAObservation obs; + + // Load counts map in CTA observation + obs.load(filename); + + // Set background cube from counts map + const GCTAEventCube* cube = dynamic_cast(obs.events()); + + // 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 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); + } + + } // endif: observation contained an events cube + + // ... 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); + } + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctbkgcube::init_members(void) +{ + // Initialise members + m_read_ahead = false; + m_obs.clear(); + m_bkgcube.clear(); + m_bkgmdl.clear(); + m_ebounds.clear(); + + // Set logger properties + log.date(true); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctbkgcube::copy_members(const ctbkgcube& app) +{ + // Copy members + m_read_ahead = app.m_read_ahead; + m_obs = app.m_obs; + m_bkgcube = app.m_bkgcube; + m_bkgmdl = app.m_bkgmdl; + m_ebounds = app.m_ebounds; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctbkgcube::free_members(void) +{ + // Write separator into logger + if (logTerse()) { + log << std::endl; + } + + // Return + return; +} diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp new file mode 100644 index 0000000..e488e19 --- /dev/null +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -0,0 +1,96 @@ +/*************************************************************************** + * ctbkgcube - CTA background cube tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014 by Chia-Chun Lu * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctbkgcube.hpp + * @brief CTA background cube tool definition + * @author Chia-Chun Lu + */ + +#ifndef CTBKGCUBE_HPP +#define CTBKGCUBE_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTBKGCUBE_NAME "ctbkgcube" +#define CTBKGCUBE_VERSION "00-01-00" + + +/***********************************************************************//** + * @class ctbkgcube + * + * @brief CTA background cube tool + ***************************************************************************/ +class ctbkgcube : public GApplication { + +public: + // Constructors and destructors + ctbkgcube(void); + explicit ctbkgcube(const GObservations& obs); + ctbkgcube(int argc, char *argv[]); + ctbkgcube(const ctbkgcube& app); + virtual ~ctbkgcube(void); + + // Operators + ctbkgcube& operator=(const ctbkgcube& app); + + // Methods + void clear(void); + void execute(void); + void run(void); + void save(void); + const GObservations& obs(void) const; + void get_parameters(void); + void get_obs(void); + void get_ebounds(void); + void set_from_cntmap(const std::string& filename); + void fill_cube(GCTAObservation* obs); + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctbkgcube& app); + void free_members(void); + + // Protected members + bool m_read_ahead; //!< Read ahead parameters + GObservations m_obs; //!< Observation container + GCTAEventCube m_bkgcube; //!< Background cube + GModels m_bkgmdl; //!< CTA background models + GEbounds m_ebounds; //!< Energy boundaries +}; + + +/***********************************************************************//** + * @brief Return observation container + * + * @return Observation container. + ***************************************************************************/ +inline +const GObservations& ctbkgcube::obs(void) const +{ + // Return observation container + return m_obs; +} + +#endif /* CTBKGCUBE_HPP */ diff --git a/src/ctbkgcube/ctbkgcube.par b/src/ctbkgcube/ctbkgcube.par new file mode 100644 index 0000000..1f89a76 --- /dev/null +++ b/src/ctbkgcube/ctbkgcube.par @@ -0,0 +1,50 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctbkgcube.par # +# Version: 00-01-00 # +# # +# Author: Chia-Chun Lu # +# c.lu@mpi-hd.mpg.de # +# MPIK # +# # +# Purpose: Parameter file for the ctbkgcube application # +# # +############################################################################# +# +# General parameters +#=================== +infile, f, a, "NONE",,, "Event list, counts cube or observation definition file" +bkgmdl, f, a, "NONE",,, "Background model XML file" +cntmap, f, a, "NONE",,, "Counts cube for background cube definition" +outfile, f, a, "bkgcube.fits",,, "Output background cube file" +# +# Energy binning parameters +#========================== +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins" +emin, r, a, 0.1,,,"Start value for first energy bin in TeV" +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV" +enumbins, i, a, 20,,,"Number of energy bins" +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" +# +# Spatial binning 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)" +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)" +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees" +proj, s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN" +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, "ql",,, "Mode of automatic parameters" +logfile, f, h, "ctbkgcube.log",,, "Log filename" diff --git a/src/ctbkgcube/main.cpp b/src/ctbkgcube/main.cpp new file mode 100644 index 0000000..93c4277 --- /dev/null +++ b/src/ctbkgcube/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctbkgcube - CTA background cube tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014 by Chia-Chun Lu * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctbkgcube/main.cpp + * @brief CTA background cube tool main code + * @author Chia-Chun Lu + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctbkgcube.hpp" + + +/***********************************************************************//** + * @brief Main entry point of application + * + * @param[in] argc Number of command line arguments. + * @param[in] argv Command line arguments. + * + * This is the main entry point of the ctbkgcube application. It allocates a + * ctbkgcube object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctbkgcube application(argc, argv); + + // Run application + try { + // Execute application + application.execute(); + + // Signal success + rc = 0; + } + catch (std::exception &e) { + + // Extract error message + std::string message = e.what(); + std::string signal = "*** ERROR encounterted in the execution of" + " ctbkgcube. Run aborted ..."; + + // Write error in logger + application.log << signal << std::endl; + application.log << message << std::endl; + + // Write error on standard output + std::cout << signal << std::endl; + std::cout << message << std::endl; + + } // endcatch: catched any application error + + // Return + return rc; +} diff --git a/src/ctexpcube/ctexpcube.par b/src/ctexpcube/ctexpcube.par index 4560351..76fd8a6 100644 --- a/src/ctexpcube/ctexpcube.par +++ b/src/ctexpcube/ctexpcube.par @@ -17,7 +17,7 @@ # # General parameters #=================== -infile, f, a, "NONE",,, "Observation definition file" +infile, f, a, "NONE",,, "Event list, counts cube or observation definition file" cntmap, f, a, "NONE",,, "Counts map for exposure cube definition" outfile, f, a, "expcube.fits",,, "Output exposure cube file" caldb, s, a, "dummy",,, "Calibration database" diff --git a/src/ctpsfcube/ctpsfcube.par b/src/ctpsfcube/ctpsfcube.par index 8a899f8..22cb4ef 100644 --- a/src/ctpsfcube/ctpsfcube.par +++ b/src/ctpsfcube/ctpsfcube.par @@ -17,7 +17,7 @@ # # General parameters #=================== -infile, f, a, "NONE",,, "Observation definition file" +infile, f, a, "NONE",,, "Event list, counts cube or observation definition file" cntmap, f, a, "NONE",,, "Counts map for exposure cube definition" outfile, f, a, "psfcube.fits",,, "Output exposure cube file" caldb, s, a, "dummy",,, "Calibration database" diff --git a/test/test_ctbin.py b/test/test_ctbin.py index e91d2dd..3027938 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ========================= # +# Test class for ctbin tool # +# ========================= # class Test(gammalib.GPythonTestSuite): """ Test class for ctbin tool. diff --git a/test/test_ctbkgcube.py b/test/test_ctbkgcube.py new file mode 100755 index 0000000..6c9b27d --- /dev/null +++ b/test/test_ctbkgcube.py @@ -0,0 +1,101 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the ctbkgcube tool. +# +# Copyright (C) 2014 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools + + +# ============================= # +# Test class for ctbkgcube tool # +# ============================= # +class Test(gammalib.GPythonTestSuite): + """ + Test class for ctbkgcube tool. + """ + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + gammalib.GPythonTestSuite.__init__(self) + + # Set members + self.events_name = "data/crab_events.fits" + self.bkg_model = "data/crab.xml" + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name("ctbkgcube") + + # Append tests + self.append(self.test_functional, "Test ctbkgcube functionality") + + # Return + return + + # Test ctbkgcube functionnality + def test_functional(self): + """ + Test ctbkgcube functionnality. + """ + # Set-up ctbkgcube + bkgcube = ctools.ctbkgcube() + bkgcube["infile"].filename(self.events_name) + bkgcube["bkgmdl"].filename(self.bkg_model) + bkgcube["cntmap"].filename("NONE") + bkgcube["outfile"].filename("psfcube.fits") + bkgcube["ebinalg"].string("LOG") + bkgcube["emin"].real(0.1) + bkgcube["emax"].real(100.0) + bkgcube["enumbins"].integer(20) + bkgcube["nxpix"].integer(10) + bkgcube["nypix"].integer(10) + bkgcube["binsz"].real(0.4) + bkgcube["coordsys"].string("CEL") + bkgcube["proj"].string("CAR") + bkgcube["xref"].real(83.63) + bkgcube["yref"].real(22.01) + + # Run tool + self.test_try("Run ctbkgcube") + try: + bkgcube.run() + self.test_try_success() + except: + self.test_try_failure("Exception occured in ctbkgcube.") + + # Save background cube + self.test_try("Save background cube") + try: + bkgcube.save() + self.test_try_success() + except: + self.test_try_failure("Exception occured in saving background cube.") + + # Return + return diff --git a/test/test_ctcubemask.py b/test/test_ctcubemask.py index fdcea06..300ab5b 100755 --- a/test/test_ctcubemask.py +++ b/test/test_ctcubemask.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ============================== # +# Test class for ctcubemask tool # +# ============================== # class Test(gammalib.GPythonTestSuite): """ Test class for ctcubemask tool. diff --git a/test/test_ctexpcube.py b/test/test_ctexpcube.py index 9a43fab..d2a55b1 100755 --- a/test/test_ctexpcube.py +++ b/test/test_ctexpcube.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ============================= # +# Test class for ctexpcube tool # +# ============================= # class Test(gammalib.GPythonTestSuite): """ Test class for ctexpcube tool. diff --git a/test/test_ctlike.py b/test/test_ctlike.py index 0e48c3b..1f743fa 100755 --- a/test/test_ctlike.py +++ b/test/test_ctlike.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ========================== # +# Test class for ctlike tool # +# ========================== # class Test(gammalib.GPythonTestSuite): """ Test class for ctlike tool. diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 17d43c4..996226c 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# =========================== # +# Test class for ctmodel tool # +# =========================== # class Test(gammalib.GPythonTestSuite): """ Test class for ctmodel tool. diff --git a/test/test_ctools.sh b/test/test_ctools.sh index 8ea3d16..ed940f1 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -42,6 +42,7 @@ ctselect=../src/ctselect/ctselect ctskymap=../src/ctskymap/ctskymap ctexpcube=../src/ctexpcube/ctexpcube ctpsfcube=../src/ctpsfcube/ctpsfcube +ctbkgcube=../src/ctbkgcube/ctbkgcube ctcubemask=../src/ctcubemask/ctcubemask @@ -444,6 +445,51 @@ $ECHO " ok" # +# Test ctbkgcube +# ============== +$ECHO -n "Test ctbkgcube: " +# +# Run 1 +$ctbkgcube infile="data/crab_events.fits" \ + bkgmdl="data/crab.xml" \ + cntmap="NONE" \ + outfile="bkgcube1.fits" \ + emin=0.1 \ + emax=100.0 \ + enumbins=20 \ + nxpix=10 \ + nypix=10 \ + binsz=0.4 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" +$ECHO -n "." +if [ -s "bkgcube1.fits" ] +then + $ECHO -n "." +else + $ECHO " bkgcube1.fits file is not valid" + exit 1 +fi +# +# Run 2 +$ctbkgcube infile="data/crab_events.fits" \ + bkgmdl="data/crab.xml" \ + cntmap="data/crab_cntmap.fits" \ + outfile="bkgcube2.fits" +$ECHO -n "." +if [ -s "bkgcube2.fits" ] +then + $ECHO -n "." +else + $ECHO " bkgcube2.fits file is not valid" + exit 1 +fi +$ECHO " ok" + + +# # Test ctcubemask # =============== $ECHO -n "Test ctcubemask: " diff --git a/test/test_ctpsfcube.py b/test/test_ctpsfcube.py index fa71a04..5870619 100755 --- a/test/test_ctpsfcube.py +++ b/test/test_ctpsfcube.py @@ -22,9 +22,9 @@ import gammalib import ctools -# ============================ # -# Test class for ctobssim tool # -# ============================ # +# ============================= # +# Test class for ctpsfcube tool # +# ============================= # class Test(gammalib.GPythonTestSuite): """ Test class for ctpsfcube tool. diff --git a/test/test_ctselect.py b/test/test_ctselect.py index ad3033e..3fbd776 100755 --- a/test/test_ctselect.py +++ b/test/test_ctselect.py @@ -23,7 +23,7 @@ import ctools # ============================ # -# Test class for ctobssim tool # +# Test class for ctselect tool # # ============================ # class Test(gammalib.GPythonTestSuite): """ diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py index c4236ed..5dee2cd 100755 --- a/test/test_ctskymap.py +++ b/test/test_ctskymap.py @@ -23,7 +23,7 @@ import ctools # ============================ # -# Test class for ctobssim tool # +# Test class for ctskymap tool # # ============================ # class Test(gammalib.GPythonTestSuite): """ diff --git a/test/test_python.py b/test/test_python.py index 4d1de1c..8ab4a76 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -28,9 +28,10 @@ import test_ctbin import test_ctlike import test_ctmodel import test_ctskymap -import test_ctcubemask import test_ctexpcube import test_ctpsfcube +import test_ctbkgcube +import test_ctcubemask import test_pipelines @@ -51,9 +52,10 @@ if __name__ == '__main__': suite_ctlike = test_ctlike.Test() suite_ctmodel = test_ctmodel.Test() suite_ctskymap = test_ctskymap.Test() - suite_ctcubemask = test_ctcubemask.Test() suite_ctexpcube = test_ctexpcube.Test() suite_ctpsfcube = test_ctpsfcube.Test() + suite_ctbkgcube = test_ctbkgcube.Test() + suite_ctcubemask = test_ctcubemask.Test() suite_pipelines = test_pipelines.Test() # Setup unit tests @@ -63,9 +65,10 @@ if __name__ == '__main__': suite_ctlike.set() suite_ctmodel.set() suite_ctskymap.set() - suite_ctcubemask.set() suite_ctexpcube.set() suite_ctpsfcube.set() + suite_ctbkgcube.set() + suite_ctcubemask.set() suite_pipelines.set() # Append tests to container @@ -75,9 +78,10 @@ if __name__ == '__main__': suites.append(suite_ctlike) suites.append(suite_ctmodel) suites.append(suite_ctskymap) - suites.append(suite_ctcubemask) suites.append(suite_ctexpcube) suites.append(suite_ctpsfcube) + suites.append(suite_ctbkgcube) + suites.append(suite_ctcubemask) suites.append(suite_pipelines) # Set PFILES environment variable