diff --git a/pyext/ctexpcube.i b/pyext/ctexpcube.i
new file mode 100644
index 0000000..ffb2fcd
--- /dev/null
+++ b/pyext/ctexpcube.i
@@ -0,0 +1,72 @@
+/***************************************************************************
+ * ctexpcube - CTA exposure cube tool *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2014 by 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 . *
+ * *
+ ***************************************************************************/
+/**
+ * @file ctexpcube.i
+ * @brief CTA exposure cube tool definition
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctexpcube.hpp"
+#include "GTools.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctexpcube
+ *
+ * @brief CTA exposure cube tool
+ ***************************************************************************/
+class ctexpcube : public GApplication {
+
+public:
+ // Constructors and destructors
+ ctexpcube(void);
+ explicit ctexpcube(GObservations obs);
+ ctexpcube(int argc, char *argv[]);
+ ctexpcube(const ctexpcube& app);
+ virtual ~ctexpcube(void);
+
+ // Operators
+ ctexpcube& operator=(const ctexpcube& app);
+
+ // Methods
+ void clear(void);
+ void execute(void);
+ void run(void);
+ void save(void);
+ GObservations& obs(void);
+ 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 exposure map tool Python extension
+ ***************************************************************************/
+%extend ctexpcube {
+ ctexpcube copy() {
+ return (*self);
+ }
+}
diff --git a/pyext/ctools/ctools.i b/pyext/ctools/ctools.i
index a72ed96..5b762e1 100755
--- a/pyext/ctools/ctools.i
+++ b/pyext/ctools/ctools.i
@@ -1,7 +1,7 @@
/***************************************************************************
* ctools - SWIG file *
* ----------------------------------------------------------------------- *
- * copyright (C) 2010-2013 by Juergen Knoedlseder *
+ * copyright (C) 2010-2014 by Juergen Knoedlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -51,3 +51,4 @@
%include "ctmodel.i"
%include "ctselect.i"
%include "ctskymap.i"
+%include "ctexpcube.i"
diff --git a/src/ctexpcube/ctexpcube b/src/ctexpcube/ctexpcube
index 0e16021..ff4c7de 100755
--- a/src/ctexpcube/ctexpcube
+++ b/src/ctexpcube/ctexpcube
diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp
index 65fd767..0881abb 100644
--- a/src/ctexpcube/ctexpcube.cpp
+++ b/src/ctexpcube/ctexpcube.cpp
@@ -31,13 +31,12 @@
#include
#include "ctexpcube.hpp"
#include "GTools.hpp"
+#include "GWcs.hpp"
/* __ Method name definitions ____________________________________________ */
#define G_GET_EBOUNDS "ctexpcube::get_ebounds()"
-
-
-#define G_SETUP_OBS "ctexpcube::setup_obs()"
-#define G_MODEL_MAP "ctexpcube::model_map(GCTAObservation*)"
+#define G_GET_PARAMETERS "ctexpcube::get_parameters()"
+#define G_SET_FROM_CNTMAP "ctexpcube::set_from_cntmap(std::string&)"
/* __ Debug definitions __________________________________________________ */
@@ -276,37 +275,8 @@ void ctexpcube::run(void)
log.header1("Generate exposure 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");
- }
- }
-
- // Increment number of observations
- n_observations++;
-
- // Generate exposure cube
- //model_map(obs, m_obs.models());
-
- } // endif: CTA observation found
-
- } // endfor: looped over observations
+ // Fill exposure
+ m_expcube.fill(m_obs);
// Restore energy dispersion flag for all CTA observations
for (int i = 0; i < m_obs.size(); ++i) {
@@ -349,24 +319,27 @@ void ctexpcube::save(void)
*
* Get all task parameters from parameter file or (if required) by querying
* the user. The parameters are read in the correct order.
+ *
+ * @todo Setup exposure cube from counts map
***************************************************************************/
void ctexpcube::get_parameters(void)
{
// If we do not have any observations in the container then get an
// input file name or observation descriptor file
if (m_obs.size() == 0) {
- m_infile = (*this)["infile"].filename();
+ get_obs();
}
- // Read response parameters
- m_caldb = (*this)["caldb"].string();
- m_irf = (*this)["irf"].string();
+ // 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 exposure cube from
// the user parameters
- m_cntmap = (*this)["cntmap"].filename();
- if ((m_cntmap == "NONE") || (gammalib::strip_whitespace(m_cntmap) == "")) {
+ 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();
@@ -388,6 +361,10 @@ void ctexpcube::get_parameters(void)
// ... otherwise setup the exposure cube from the counts map
else {
+
+ // Set exposure cube from counts map
+ set_from_cntmap(cntmap);
+
}
// Read output filename (if needed)
@@ -401,8 +378,98 @@ void ctexpcube::get_parameters(void)
/***********************************************************************//**
+ * @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 {
+
+ // 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 ctexpcube::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 it isn't set already
+ if (obs->response().aeff() == NULL) {
+
+ // 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 ctexpcube::get_ebounds(void)
@@ -467,6 +534,74 @@ void ctexpcube::get_ebounds(void)
}
+/***********************************************************************//**
+ * @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);
+ }
+
+ } // 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 =
@@ -479,10 +614,7 @@ void ctexpcube::get_ebounds(void)
void ctexpcube::init_members(void)
{
// Initialise members
- m_infile.clear();
m_outfile.clear();
- m_caldb.clear();
- m_irf.clear();
m_apply_edisp = false;
// Initialise protected members
@@ -507,10 +639,7 @@ void ctexpcube::init_members(void)
void ctexpcube::copy_members(const ctexpcube& app)
{
// Copy attributes
- m_infile = app.m_infile;
m_outfile = app.m_outfile;
- m_caldb = app.m_caldb;
- m_irf = app.m_irf;
m_apply_edisp = app.m_apply_edisp;
// Copy protected members
diff --git a/src/ctexpcube/ctexpcube.hpp b/src/ctexpcube/ctexpcube.hpp
index 4654836..fcd270f 100644
--- a/src/ctexpcube/ctexpcube.hpp
+++ b/src/ctexpcube/ctexpcube.hpp
@@ -61,7 +61,10 @@ public:
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
@@ -70,11 +73,7 @@ protected:
void free_members(void);
// User parameters
- std::string m_infile; //!< Input counts map or XML file
- std::string m_cntmap; //!< Counts map for exposure cube definition
std::string m_outfile; //!< Output exposure cube file
- std::string m_caldb; //!< Calibration database
- std::string m_irf; //!< Instrument response functions
bool m_apply_edisp; //!< Apply energy dispersion?
// Protected members
diff --git a/src/ctexpcube/ctexpcube.par b/src/ctexpcube/ctexpcube.par
index b1b98a0..4560351 100644
--- a/src/ctexpcube/ctexpcube.par
+++ b/src/ctexpcube/ctexpcube.par
@@ -17,7 +17,7 @@
#
# General parameters
#===================
-infile, f, a, "NONE",,, "Input event list, counts map or observation definition file"
+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"
diff --git a/test/test_ctools.sh b/test/test_ctools.sh
index 9f0ea78..72bb3f8 100755
--- a/test/test_ctools.sh
+++ b/test/test_ctools.sh
@@ -2,7 +2,7 @@
# ==========================================================================
# This script tests all ctools
#
-# Copyright (C) 2011-2012 Juergen Knoedlseder
+# Copyright (C) 2011-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
@@ -40,6 +40,7 @@ ctmodel=../src/ctmodel/ctmodel
ctobssim=../src/ctobssim/ctobssim
ctselect=../src/ctselect/ctselect
ctskymap=../src/ctskymap/ctskymap
+ctexpcube=../src/ctexpcube/ctexpcube
#
@@ -340,4 +341,49 @@ else
exit 1
fi
$ECHO " ok"
-#sleep 100
+
+
+#
+# Test ctexpcube
+# ==============
+$ECHO -n "Test ctexpcube: "
+#
+# Run 1
+$ctexpcube infile="data/crab_events.fits" \
+ cntmap="NONE" \
+ outfile="expcube1.fits" \
+ caldb="irf" \
+ irf="cta_dummy_irf" \
+ emin=0.1 \
+ emax=100.0 \
+ enumbins=20 \
+ nxpix=200 \
+ nypix=200 \
+ binsz=0.02 \
+ coordsys="CEL" \
+ xref=83.63 \
+ yref=22.01 \
+ proj="CAR"
+$ECHO -n "."
+if [ -s "expcube1.fits" ]
+then
+ $ECHO -n "."
+else
+ $ECHO " expcube1.fits file is not valid"
+ exit 1
+fi
+#
+# Run 2
+$ctexpcube infile="data/crab_events.fits" \
+ cntmap="data/crab_cntmap.fits" \
+ outfile="expcube2.fits" \
+ caldb="irf" \
+ irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "expcube2.fits" ]
+then
+ $ECHO -n "."
+else
+ $ECHO " expcube2.fits file is not valid"
+ exit 1
+fi