diff --git a/configure.ac b/configure.ac
index 9023bb7..706bf89 100644
--- a/configure.ac
+++ b/configure.ac
@@ -883,6 +883,7 @@ AC_CONFIG_FILES([Makefile
src/ctskymap/Makefile
src/ctselect/Makefile
src/ctexpcube/Makefile
+ src/ctpsfcube/Makefile
src/ctcubemask/Makefile
src/ctools-setup
test/Makefile
diff --git a/pyext/ctexpcube.i b/pyext/ctexpcube.i
index 1d435a0..446ec5e 100644
--- a/pyext/ctexpcube.i
+++ b/pyext/ctexpcube.i
@@ -60,7 +60,7 @@ public:
/***********************************************************************//**
- * @brief CTA exposure map tool Python extension
+ * @brief CTA exposure cube tool Python extension
***************************************************************************/
%extend ctexpcube {
ctexpcube copy() {
diff --git a/pyext/ctpsfcube.i b/pyext/ctpsfcube.i
new file mode 100644
index 0000000..a878643
--- /dev/null
+++ b/pyext/ctpsfcube.i
@@ -0,0 +1,71 @@
+/***************************************************************************
+ * ctpsfcube - CTA PSF 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 ctpsfcube.i
+ * @brief CTA PSF cube tool definition
+ * @author Chia-Chun Lu
+ */
+
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctpsfcube.hpp"
+#include "GTools.hpp"
+%}
+
+/***********************************************************************//**
+ * @class ctpsfcube
+ *
+ * @brief CTA PSF cube tool
+ ***************************************************************************/
+class ctpsfcube : public GApplication {
+
+public:
+ // Constructors and destructors
+ ctpsfcube(void);
+ explicit ctpsfcube(GObservations obs);
+ ctpsfcube(int argc, char *argv[]);
+ ctpsfcube(const ctpsfcube& app);
+ virtual ~ctpsfcube(void);
+
+ // Operators
+ ctpsfcube& operator=(const ctpsfcube& app);
+
+ // Methods
+ void clear(void);
+ void execute(void);
+ void run(void);
+ void save(void);
+ GObservations& obs(void) { return m_obs; }
+ void get_parameters(void);
+ void get_obs(void);
+ void set_response(void);
+ void get_ebounds(void);
+ void set_from_cntmap(const std::string& filename);
+};
+
+/***********************************************************************//**
+ * @brief CTA PSF cube tool Python extension
+ ***************************************************************************/
+%extend ctpsfcube {
+ ctpsfcube copy() {
+ return (*self);
+ }
+}
diff --git a/src/Makefile.am b/src/Makefile.am
index 7d789f6..75dcfd4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -12,7 +12,8 @@ SUBDIRS = ctbin \
ctselect \
ctskymap \
ctexpcube \
- ctcubemask
+ ctpsfcube \
+ ctcubemask
# Scripts to be distributed
scripts = ctools-setup \
@@ -31,6 +32,7 @@ libctools_la_LIBADD = ctbin/libctbin.la \
ctselect/libctselect.la \
ctskymap/libctskymap.la \
ctexpcube/libctexpcube.la \
+ ctpsfcube/libctpsfcube.la \
ctcubemask/libctcubemask.la
# Distribute scripts
diff --git a/src/ctcubemask/ctcubemask.cpp b/src/ctcubemask/ctcubemask.cpp
index 75587ad..cdf2da9 100644
--- a/src/ctcubemask/ctcubemask.cpp
+++ b/src/ctcubemask/ctcubemask.cpp
@@ -201,7 +201,7 @@ void ctcubemask::clear(void)
/***********************************************************************//**
* @brief Execute application
*
- * This method performs the event selection and saves the result
+ * This method applies the cube mask and saves the result.
***************************************************************************/
void ctcubemask::execute(void)
{
@@ -286,7 +286,7 @@ void ctcubemask::run(void)
// Save event file name (for possible saving)
m_infiles[i] = obs->eventfile();
- // Load observation from temporary file, including event selection
+ // Apply mask on the event cube
apply_mask(obs);
} // endif: had a CTA observation
@@ -302,7 +302,7 @@ void ctcubemask::run(void)
// Write observation(s) into logger
if (logTerse()) {
log << std::endl;
- log.header1("Observations after selection");
+ log.header1("Observations after masked");
log << m_obs << std::endl;
}
@@ -312,22 +312,22 @@ void ctcubemask::run(void)
/***********************************************************************//**
- * @brief Save the selected event list(s)
+ * @brief Save the masked event cube(s)
*
- * This method saves the selected event list(s) into FITS file(s). There are
+ * This method saves the masked event cube(s) into FITS file(s). There are
* two modes, depending on the m_use_xml flag.
*
- * If m_use_xml is true, all selected event list(s) will be saved into FITS
+ * If m_use_xml is true, all masked event cube(s) will be saved into FITS
* files, where the output filenames are constructued from the input
* filenames by prepending the m_prefix string to name. Any path information
- * will be stripped form the input name, hence event files will be written
+ * will be stripped form the input name, hence event cube files will be written
* into the local working directory (unless some path information is present
* in the prefix). In addition, an XML file will be created that gathers
- * the filename information for the selected event list(s). If an XML file
+ * the filename information for the masked event cube(s). If an XML file
* was present on input, all metadata information will be copied from this
* input file.
*
- * If m_use_xml is false, the selected event list will be saved into a FITS
+ * If m_use_xml is false, the masked event cubes will be saved into a FITS
* file.
***************************************************************************/
void ctcubemask::save(void)
@@ -362,10 +362,10 @@ void ctcubemask::save(void)
* @brief Get application parameters
*
* 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.
+ * the user.
*
* This method also loads observations if no observations are yet allocated.
- * Observations are either loaded from a single CTA even list, or from a
+ * Observations are either loaded from a single CTA even cube, or from a
* XML file using the metadata information that is stored in that file.
***************************************************************************/
void ctcubemask::get_parameters(void)
@@ -383,7 +383,7 @@ void ctcubemask::get_parameters(void)
// Try first to open as FITS file
try {
- // Load event list in CTA observation
+ // Load event cube in CTA observation
obs.load(m_infile);
// Append CTA observation to container
@@ -430,11 +430,13 @@ void ctcubemask::get_parameters(void)
/***********************************************************************//**
- * @brief Apply mask to counts cube
+ * @brief Mask event cube
*
* @param[in] obs CTA observation.
*
- * Masks out event bins by setting their content to -1.
+ * Mask an event cube. The mask sets pixel values outside the region of interest
+ * or the energy of interest to -1. These pixels will be ignored in likelihood
+ * fitting.
***************************************************************************/
void ctcubemask::apply_mask(GCTAObservation* obs)
{
diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp
index 39ff33e..89fd687 100644
--- a/src/ctexpcube/ctexpcube.cpp
+++ b/src/ctexpcube/ctexpcube.cpp
@@ -223,7 +223,7 @@ void ctexpcube::execute(void)
/***********************************************************************//**
- * @brief Generate the model map(s)
+ * @brief Generate the exposure cube(s).
*
* This method reads the task parameters from the parfile, sets up the
* observation container, loops over all CTA observations in the container
diff --git a/src/ctpsfcube/Makefile.am b/src/ctpsfcube/Makefile.am
new file mode 100644
index 0000000..3c402a0
--- /dev/null
+++ b/src/ctpsfcube/Makefile.am
@@ -0,0 +1,21 @@
+## Process this file with automake to produce Makefile.in
+
+# Program name
+bin_PROGRAMS = ctpsfcube
+
+# Sources
+sources = ctpsfcube.cpp
+
+# Libtool library
+noinst_LTLIBRARIES = libctpsfcube.la
+libctpsfcube_la_SOURCES = $(sources)
+
+# Program definition
+ctpsfcube_SOURCES = main.cpp
+ctpsfcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
+ctpsfcube_CPPFLAGS = @CPPFLAGS@
+ctpsfcube_LDADD = libctpsfcube.la
+
+# Parameter files
+pfilesdir = $(prefix)/syspfiles
+pfiles_DATA = ctpsfcube.par
diff --git a/src/ctpsfcube/ctpsfcube b/src/ctpsfcube/ctpsfcube
new file mode 100755
index 0000000..c65d28d
--- /dev/null
+++ b/src/ctpsfcube/ctpsfcube
diff --git a/src/ctpsfcube/ctpsfcube.cpp b/src/ctpsfcube/ctpsfcube.cpp
new file mode 100644
index 0000000..c260a38
--- /dev/null
+++ b/src/ctpsfcube/ctpsfcube.cpp
@@ -0,0 +1,672 @@
+/***************************************************************************
+ * ctpsfcube - CTA PSF 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 ctpsfcube.cpp
+ * @brief CTA PSF cube tool implementation
+ * @author Chia-Chun Lu
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include
+#endif
+#include
+#include "ctpsfcube.hpp"
+#include "GTools.hpp"
+#include "GWcs.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_GET_EBOUNDS "ctpsfcube::get_ebounds()"
+#define G_GET_PARAMETERS "ctpsfcube::get_parameters()"
+#define G_SET_FROM_CNTMAP "ctpsfcube::set_from_cntmap(std::string&)"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+
+/*==========================================================================
+ = =
+ = Constructors/destructors =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+ctpsfcube::ctpsfcube(void) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
+{
+ // Initialise members
+ init_members();
+
+ // Write header into logger
+ log_header();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Observations constructor
+ *
+ * This method creates an instance of the class by copying an existing
+ * observations container.
+ ***************************************************************************/
+ctpsfcube::ctpsfcube(GObservations obs) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_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.
+ ***************************************************************************/
+ctpsfcube::ctpsfcube(int argc, char *argv[]) :
+ GApplication(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, argc, argv)
+{
+ // Initialise members
+ init_members();
+
+ // Write header into logger
+ log_header();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app ctpsfcube application.
+ ***************************************************************************/
+ctpsfcube::ctpsfcube(const ctpsfcube& app) : GApplication(app)
+{
+ // Initialise members
+ init_members();
+
+ // Copy members
+ copy_members(app);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+ctpsfcube::~ctpsfcube(void)
+{
+ // Free members
+ free_members();
+
+ // Return
+ return;
+}
+
+
+/*==========================================================================
+ = =
+ = Operators =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app ctpsfcube application.
+ * @return Returns ctpsfcube application.
+ ***************************************************************************/
+ctpsfcube& ctpsfcube::operator=(const ctpsfcube& 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 ctpsfcube::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 ctpsfcube class. It is invoked
+ * when the executable is called from command line. The method generates
+ * the PSF cube and saves the result.
+ ***************************************************************************/
+void ctpsfcube::execute(void)
+{
+ // Signal that some parameters should be read ahead
+ m_read_ahead = true;
+
+ // Create the PSF cube
+ run();
+
+ // Save the PSF cube into the output FITS file
+ save();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Generate the model map(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 PSF cube from the CTA observations.
+ ***************************************************************************/
+void ctpsfcube::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;
+ }
+
+ // Set energy dispersion flag for all CTA observations and save old
+ // values in save_edisp vector
+ std::vector save_edisp;
+ save_edisp.assign(m_obs.size(), false);
+ for (int i = 0; i < m_obs.size(); ++i) {
+ GCTAObservation* obs = dynamic_cast(m_obs[i]);
+ if (obs != NULL) {
+ save_edisp[i] = obs->response()->apply_edisp();
+ obs->response()->apply_edisp(m_apply_edisp);
+ }
+ }
+
+ // 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("Generate PSF cube");
+ }
+
+ // Fill PSF
+ m_psfcube.fill(m_obs);
+
+ // Restore energy dispersion flag for all CTA observations
+ for (int i = 0; i < m_obs.size(); ++i) {
+ GCTAObservation* obs = dynamic_cast(m_obs[i]);
+ if (obs != NULL) {
+ obs->response()->apply_edisp(save_edisp[i]);
+ }
+ }
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save PSF cube
+ *
+ ***************************************************************************/
+void ctpsfcube::save(void)
+{
+ // Write header
+ if (logTerse()) {
+ log << std::endl;
+ log.header1("Save PSF cube");
+ }
+
+ // Get output filename
+ m_outfile = (*this)["outfile"].filename();
+
+ // Save PSF cube
+ m_psfcube.save(m_outfile, clobber());
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * Get all task parameters from parameter file or (if required) by querying
+ * the user. The parameters are read in the correct order.
+ *
+ * @todo Setup PSF cube from counts map
+ ***************************************************************************/
+void ctpsfcube::get_parameters(void)
+{
+ // If we do not have any observations in the container then get an
+ // input file name or observation descriptor file
+ if (m_obs.size() == 0) {
+ get_obs();
+ }
+
+ // Make sure that response is set
+ set_response();
+
+ // Read energy dispersion flag
+ m_apply_edisp = (*this)["edisp"].boolean();
+
+ // 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)["dmax"].real();
+ int ndbins = (*this)["ndbins"].integer();
+
+ // Get energy definition
+ get_ebounds();
+
+ // Define exposure cube
+ m_psfcube = GCTAMeanPsf(wcs, coordsys, xref, yref,
+ -binsz, binsz, nxpix, nypix,
+ m_ebounds, dmax, ndbins);
+ }
+
+ // ... otherwise setup the exposure cube from the counts map
+ else {
+
+ // Set PSF cube from counts map
+ set_from_cntmap(cntmap);
+
+ }
+
+ // Read output filename (if needed)
+ if (m_read_ahead) {
+ m_outfile = (*this)["outfile"].filename();
+ }
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @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();
+
+ // Try first to open as FITS file
+ try {
+
+ // Allocate CTA observation
+ GCTAObservation obs;
+
+ // Load input file in CTA observation
+ obs.load(filename);
+
+ // Append CTA observation to container
+ m_obs.append(obs);
+
+
+ }
+
+ // ... otherwise try to open as XML file
+ catch (GException::fits_open_error &e) {
+
+ // Load observations from XML file
+ m_obs.load(filename);
+
+ }
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set observation response
+ *
+ * Set response for all observations that so have no response.
+ ***************************************************************************/
+void ctpsfcube::set_response(void)
+{
+ // Loop over all observations
+ for (int i = 0; i < m_obs.size(); ++i) {
+
+ // Is this observation a CTA observation?
+ GCTAObservation* obs = dynamic_cast(m_obs[i]);
+
+ // Yes ...
+ if (obs != NULL) {
+
+ // Set response if we don't have one
+ if (!obs->hasresponse()) {
+
+ // Load response information
+ std::string database = (*this)["caldb"].string();
+ std::string irf = (*this)["irf"].string();
+
+ // Set calibration database. If specified parameter is a
+ // directory then use this as the pathname to the calibration
+ // database. Otherwise interpret this as the instrument name,
+ // the mission being "cta"
+ GCaldb caldb;
+ if (gammalib::dir_exists(database)) {
+ caldb.rootdir(database);
+ }
+ else {
+ caldb.open("cta", database);
+ }
+
+ // Set reponse
+ obs->response(irf, caldb);
+
+ } // endif: observation already has a response
+
+ } // endif: observation was a CTA observation
+
+ } // endfor: looped over all observations
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get energy boundaries from parameters
+ *
+ * @exception GException::invalid_value
+ * Invalid extension name encountered.
+ *
+ * Get the energy boundaries from the user parameters.
+ ***************************************************************************/
+void ctpsfcube::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();
+
+ // 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 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 ctpsfcube::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();
+ double dmax = (*this)["dmax"].real();
+ int ndbins = (*this)["ndbins"].integer();
+
+ // Get energy definition
+ m_ebounds = cube->ebounds();
+
+ // Define exposure cube
+ m_psfcube = GCTAMeanPsf(proj, coordsys, xref, yref,
+ dx, dy, nx, ny,
+ m_ebounds, dmax, ndbins);
+
+ } // 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 ctpsfcube::init_members(void)
+{
+ // Initialise members
+ m_outfile.clear();
+ m_apply_edisp = false;
+
+ // Initialise protected members
+ m_read_ahead = false;
+ m_obs.clear();
+ m_psfcube.clear();
+ m_ebounds.clear();
+
+ // Set logger properties
+ log.date(true);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctpsfcube::copy_members(const ctpsfcube& app)
+{
+ // Copy attributes
+ m_outfile = app.m_outfile;
+ m_apply_edisp = app.m_apply_edisp;
+
+ // Copy protected members
+ m_read_ahead = app.m_read_ahead;
+ m_obs = app.m_obs;
+ m_psfcube = app.m_psfcube;
+ m_ebounds = app.m_ebounds;
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctpsfcube::free_members(void)
+{
+ // Write separator into logger
+ if (logTerse()) {
+ log << std::endl;
+ }
+
+ // Return
+ return;
+}
diff --git a/src/ctpsfcube/ctpsfcube.hpp b/src/ctpsfcube/ctpsfcube.hpp
new file mode 100644
index 0000000..b34c48c
--- /dev/null
+++ b/src/ctpsfcube/ctpsfcube.hpp
@@ -0,0 +1,86 @@
+/***************************************************************************
+ * ctpsfcube - CTA PSF 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 ctpsfcube.hpp
+ * @brief CTA PSF cube tool definition
+ * @author Chia-Chun Lu
+ */
+
+#ifndef CTPSFCUBE_HPP
+#define CTPSFCUBE_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+
+/* __Definitions _________________________________________________________ */
+#define CTPSFCUBE_NAME "ctpsfcube"
+#define CTPSFCUBE_VERSION "00-01-00"
+
+
+/***********************************************************************//**
+ * @class ctpsfcube
+ *
+ * @brief CTA PSF cube tool
+ ***************************************************************************/
+class ctpsfcube : public GApplication {
+
+public:
+ // Constructors and destructors
+ ctpsfcube(void);
+ explicit ctpsfcube(GObservations obs);
+ ctpsfcube(int argc, char *argv[]);
+ ctpsfcube(const ctpsfcube& app);
+ virtual ~ctpsfcube(void);
+
+ // Operators
+ ctpsfcube& operator=(const ctpsfcube& app);
+
+ // Methods
+ void clear(void);
+ void execute(void);
+ void run(void);
+ void save(void);
+ GObservations& obs(void) { return m_obs; }
+ void get_parameters(void);
+ void get_obs(void);
+ void set_response(void);
+ void get_ebounds(void);
+ void set_from_cntmap(const std::string& filename);
+
+protected:
+ // Protected methods
+ void init_members(void);
+ void copy_members(const ctpsfcube& app);
+ void free_members(void);
+
+ // User parameters
+ std::string m_outfile; //!< Output PSF cube file
+ bool m_apply_edisp; //!< Apply energy dispersion?
+
+ // Protected members
+ bool m_read_ahead; //!< Read ahead parameters
+ GObservations m_obs; //!< Observation container
+ GCTAMeanPsf m_psfcube; //!< PSF cube
+ GEbounds m_ebounds; //!< Energy boundaries
+};
+
+#endif /* CTPSFCUBE_HPP */
diff --git a/src/ctpsfcube/ctpsfcube.par b/src/ctpsfcube/ctpsfcube.par
new file mode 100644
index 0000000..1c3d219
--- /dev/null
+++ b/src/ctpsfcube/ctpsfcube.par
@@ -0,0 +1,57 @@
+#############################################################################
+# #
+# CTA Science Analysis Tools #
+# #
+# ------------------------------------------------------------------------- #
+# #
+# File: ctexpcube.par #
+# Version: 00-01-00 #
+# #
+# Author: Juergen Knoedlseder #
+# Jurgen.Knodlseder@irap.omp.eu #
+# IRAP #
+# #
+# Purpose: Parameter file for the ctexpcube application #
+# #
+#############################################################################
+#
+# General parameters
+#===================
+infile, f, a, "NONE",,, "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"
+irf, s, a, "cta_dummy_irf",,, "Instrument response function"
+edisp, b, h, no,,, "Apply energy dispersion?"
+#
+# 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"
+#
+# Histogram binning parameters
+#=============================
+dmax, r, a, 0.3,,, "Upper bound of delta"
+ndbins, i, a, 60,,, "Number of delta bins"
+#
+# 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, "ctexpcube.log",,, "Log filename"
diff --git a/src/ctpsfcube/main.cpp b/src/ctpsfcube/main.cpp
new file mode 100644
index 0000000..bb7c94c
--- /dev/null
+++ b/src/ctpsfcube/main.cpp
@@ -0,0 +1,79 @@
+/***************************************************************************
+ * ctpsfcube - CTA PSF 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 ctpsfcube/main.cpp
+ * @brief CTA PSF cube tool main code
+ * @author Chia-Chun Lu
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include
+#endif
+#include
+#include "ctpsfcube.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 ctpsfcube application. It allocates a
+ * ctpsfcube object and runs the application.
+ ***************************************************************************/
+int main (int argc, char *argv[])
+{
+ // Initialise return code
+ int rc = 1;
+
+ // Create instance of application
+ ctpsfcube 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"
+ " ctpsfcube. 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;
+}