Merge branch '1323-write-tsmap' into integration

Jürgen Knödlseder
2 parents 35cf1869 f3df2070
Showing 15 changed files with 1275 additions and 17 deletions
.gitignore
ChangeLog
NEWS
configure.ac
pyext/ctools/ctools.i
pyext/cttsmap.i
src/Makefile.am
src/cttsmap/Makefile.am
src/cttsmap/cttsmap.cpp
src/cttsmap/cttsmap.hpp
src/cttsmap/cttsmap.par
src/cttsmap/main.cpp
test/test_ctools.sh
test/test_cttsmap.py
test/test_python.py
@@ -48,6 +48,7 @@ pyext/ctools/ctools.py
 pyext/ctools/ctools_wrap.cpp
 src/ctbin/ctbin
 src/ctlike/ctlike
+src/cttsmap/cttsmap
 src/ctmodel/ctmodel
 src/ctobssim/ctobssim
 src/ctselect/ctselect
-2014-07-24  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-10-06  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.8.0 released
         ========================
@@ -10,6 +10,7 @@
         Add ctpsfcube tool (generates PSF cube)
         Add ctbkgcube tool (generates background cube)
         Add ctcubemask tool (masks event bins in event cube)
+        Add cttsmap tool (generates TS map)
         Support energy dispersion in ctobssim, ctlike and ctmodel	
         Correct ctobssim for thread safety
         Move obsutils to ctools namespace
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
 Author(s) : Juergen Knoedlseder
-24 July 2014
+6 October 2014
 1. Introduction
@@ -64,8 +64,8 @@ case of binning a single event list and using the usepnt=yes option
 (or specifying the pointing manually as the xref and yref parameters).
-7. ctlike - CTA maximum likelihood fitting
-------------------------------------------
+7. ctlike - Maximum likelihood fitting
+--------------------------------------
 ctlike now has full access to the calibration database, allowing
 the use for 2D response functions in FITS files.
@@ -73,13 +73,28 @@ ctlike now also supports the handling of energy dispersion. Energy
 dispersion can be switched on or off using the hidden edisp parameter.
-8. ctskymap - CTA sky mapping
+8. cttsmap - Generation of Test Statistics map
+----------------------------------------------
+The cttsmap tool has been added that provides a test statistics map
+for a given set of observations. This tool is not restricted to CTA
+data and will work for all types of instruments.
+
+cttsmap will move a test sources over a specified grid of positions
+and compute the likelihood difference between having the source at
+a given position and having not the source. Twice the likelihood
+difference is called the Test Statistic. The computation is done for
+a grid of positions, and the results is saved into a FITS file as
+TS map. Also the value of the fitted parameters are stored as
+extension in the FITS file.
+
+
+9. ctskymap - CTA sky mapping
 -----------------------------
 None
-9. ctmodel - CTA counts map modelling
--------------------------------------
+10. ctmodel - CTA counts map modelling
+--------------------------------------
 ctmodel now has full access to the calibration database, allowing
 the use for 2D response functions in FITS files.
@@ -87,7 +102,7 @@ ctmodel now also supports the handling of energy dispersion. Energy
 dispersion can be switched on or off using the hidden edisp parameter.
-10. ctexpcube - CTA exposure cube generation
+11. ctexpcube - CTA exposure cube generation
 --------------------------------------------
 The first version of the ctexpcube tool is included in the release.
 ctexpcube computes the exposure cube for an observation definition
@@ -95,7 +110,7 @@ file. The cube definition may be taken from a counts cube or may be
 specified by user parameters.
-11. ctpsfcube - CTA point spread function cube generation
+12. ctpsfcube - CTA point spread function cube generation
 ---------------------------------------------------------
 The first version of the ctpsfcube tool is included in the release.
 ctpsfcube computes the point spread function cube for an observation
@@ -103,7 +118,7 @@ definition file. The cube definition may be taken from a counts cube
 or may be specified by user parameters.
-12. ctbkgcube - CTA background cube generation
+13. 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 
@@ -111,7 +126,7 @@ file. The cube definition may be taken from a counts cube or may be
 specified by user parameters.
-13. ctcubemask - CTA cube filtering
+14. 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
@@ -119,12 +134,12 @@ bin content to -1. Filtering can be done by energy, for a circular
 Region of Interest (ROI) or using ds9 exclusion regions.
-14. Examples
+15. Examples
 ------------
 None
-15. Scripts
+16. Scripts
 -----------
 Add logL and event statistics to cspull output file.
@@ -144,11 +159,11 @@ The following functions have been added to obsutils:
   specified energy bins.
-16. caldb
+17. caldb
 ---------
 None
-17. models
+18. models
 ----------
 None
@@ -886,6 +886,7 @@ AC_CONFIG_FILES([Makefile
                  src/ctpsfcube/Makefile
                  src/ctbkgcube/Makefile
                  src/ctcubemask/Makefile
+		 src/cttsmap/Makefile
                  src/ctools-setup
                  test/Makefile
                  test/test_cscripts.sh
@@ -55,3 +55,4 @@
 %include "ctpsfcube.i"
 %include "ctbkgcube.i"
 %include "ctcubemask.i"
+%include "cttsmap.i"
+/***************************************************************************
+ *                      cttsmap - TS map calculation tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file cttsmap.i
+ * @brief TS map calculation tool interface definition
+ * @author Michael Mayer
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "cttsmap.hpp"
+#include "GTools.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class cttsmap
+ *
+ * @brief TS map calculation
+ ***************************************************************************/
+class cttsmap : public GApplication  {
+public:
+    // Constructors and destructors
+    cttsmap(void);
+    explicit cttsmap(const GObservations& obs);
+    cttsmap(int argc, char *argv[]);
+    cttsmap(const cttsmap& app);
+    virtual ~cttsmap(void);
+    
+    // Methods
+    void                 clear(void);
+    void                 execute(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    const GSkymap&       tsmap(void) const;
+};
+
+
+/***********************************************************************//**
+ * @brief TS map calculation tool Python extensions
+ ***************************************************************************/
+%extend cttsmap {
+    cttsmap copy() {
+        return (*self);
+    }
+}
@@ -14,7 +14,8 @@ SUBDIRS = ctbin \
           ctexpcube \
           ctpsfcube \
           ctbkgcube \
-          ctcubemask 
+          ctcubemask \
+          cttsmap
 # Scripts to be distributed
 scripts = ctools-setup \
@@ -35,7 +36,8 @@ libctools_la_LIBADD = ctbin/libctbin.la \
                       ctexpcube/libctexpcube.la \
                       ctpsfcube/libctpsfcube.la \
                       ctbkgcube/libctbkgcube.la \
-                      ctcubemask/libctcubemask.la
+                      ctcubemask/libctcubemask.la \
+                      cttsmap/libcttsmap.la
 # Distribute scripts
 EXTRA_DIST   = $(scripts)
+## Process this file with automake to produce Makefile.in
+
+# Program name
+bin_PROGRAMS = cttsmap
+
+# Sources
+sources = cttsmap.cpp
+
+# Libtool library
+noinst_LTLIBRARIES = libcttsmap.la
+libcttsmap_la_SOURCES = $(sources)
+
+# Program definition
+cttsmap_SOURCES = main.cpp
+cttsmap_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
+cttsmap_CPPFLAGS = @CPPFLAGS@
+cttsmap_LDADD = libcttsmap.la
+
+# Parameter files
+pfilesdir = $(prefix)/syspfiles
+pfiles_DATA = cttsmap.par
+/***************************************************************************
+ *                      cttsmap - TS map calculation tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file cttsmap.hpp
+ * @brief TS map calculation tool interface implementation
+ * @author Michael Mayer
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "cttsmap.hpp"
+#include "GTools.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_GET_PARAMETERS                          "cttsmap::get_parameters()"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+cttsmap::cttsmap(void) : GApplication(CTTSMAP_NAME, CTTSMAP_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.
+ ***************************************************************************/
+cttsmap::cttsmap(const GObservations& obs) :
+         GApplication(CTTSMAP_NAME, CTTSMAP_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.
+ ***************************************************************************/
+cttsmap::cttsmap(int argc, char *argv[]) :
+         GApplication(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
+{
+    // Initialise members
+    init_members();
+
+    // Write header into logger
+    log_header();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+cttsmap::cttsmap(const cttsmap& app) : GApplication(app)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(app);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+cttsmap::~cttsmap(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+cttsmap& cttsmap::operator=(const cttsmap& 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 cttsmap::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 cttsmap class. It is invoked when
+ * the executable is called from command line.
+ *
+ * The method reads the task parameters, computes the TS values, stores the
+ * results into map(s), and writes the maps into a FITS file on disk.
+ ***************************************************************************/
+void cttsmap::execute(void)
+{
+    // Signal that some parameters should be read ahead
+    m_read_ahead = true;
+
+    // Compute the maps
+    run();
+
+    // Save the maps into FITS file
+    save();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes the TS maps
+ *
+ * This method moves a source along a grid of coordinates and refits the
+ * given models including this additional source. The fit parameters of the
+ * source and its TS values are filled into respective position in the map.
+ ***************************************************************************/
+void cttsmap::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("Initialise TS map");
+    }
+
+    // Initialise maps
+    init_maps();
+
+    // Get bins to be computed
+    // To split the computation on several jobs,
+    // the user is able to provide a subset of bins
+    // which should be computed
+    int binmin = (m_binmin == -1) ? 0 : m_binmin;
+    int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax;
+
+    // Initialise optimizer
+    GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(log)
+        	                            : new GOptimizerLM();
+
+    // Store initial models
+    GModels models_orig = m_obs.models();
+
+    // Get model instance for further computations
+    GModels models = m_obs.models();
+
+    // Remove test source
+    models.remove(m_srcname);
+
+    // Get likelihood of null hypothesis if not given before
+	if (m_logL0 == 0.0) {
+
+        // Write header
+        if (logTerse()) {
+            log << std::endl;
+            log.header1("Compute NULL Hypothesis for TS computation");
+        }
+
+		// Compute likelihood without the test source
+		m_obs.models(models);
+		m_obs.optimize(*opt);
+		m_logL0 = -(opt->value());
+
+	}
+
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Generate TS map");
+    }
+
+    // Loop over grid positions
+    for (int i = binmin; i < binmax; ++i) {
+
+    	// Get the coordinate of current bin
+    	GSkyDir bincentre = m_tsmap.inx2dir(i);
+
+    	// Header for verbose logging
+    	if (logExplicit()) {
+            std::string msg = "Computing TS for bin number "+gammalib::str(i)+
+                              " at "+bincentre.print();
+            log << std::endl;
+            log.header2(msg);
+    	}
+
+    	// Add test source at current bin position
+        if (m_testsource != NULL) {
+            (*m_testsource)["RA"].value(bincentre.ra_deg());
+            (*m_testsource)["DEC"].value(bincentre.dec_deg());
+            models.append(*m_testsource);
+        }
+
+    	// Assign models to observations
+    	m_obs.models(models);
+
+    	// Optimize observation container
+    	m_obs.optimize(*opt);
+
+    	// Retrieve the Likelihood value
+    	double logL1 = -(opt->value());
+
+    	// Compute TS value
+    	double ts = 2.0 * (logL1 - m_logL0);
+
+    	// Log information
+    	if (logExplicit()) {
+    		log << " TS value ....: ";
+            log << ts << std::endl;
+    	}
+        else if (logTerse()) {
+    		log << "TS for bin number ";
+            log << i;
+            log << " at ";
+            log << bincentre.print();
+            log << ": ";
+            log << ts << std::endl;
+        }
+
+    	// Get test source model instance
+    	GModels best_fit_model = m_obs.models();
+    	GModel* testsource     = best_fit_model[m_srcname];
+
+    	// Assign values to the maps
+    	m_tsmap(i) = ts;
+
+    	// Extract fitted test source parameters
+        if (m_testsource != NULL) {
+            for (int j = 0; j < m_mapnames.size(); ++j) {
+                m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+            }
+    	}
+
+    	// Set status of bin to true
+    	m_statusmap(i) = 1.0;
+
+    	// Remove model from container
+    	models.remove(m_srcname);
+
+    } // endfor: looped over grid positions
+
+    // Bring models to initial state
+    m_obs.models(models_orig);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save maps
+ *
+ * This method saves the TS map and additional maps into a FITS file.
+ ***************************************************************************/
+void cttsmap::save(void)
+{
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Save TS map");
+    }
+
+    // Get output filename
+    m_outfile = (*this)["outfile"].filename();
+
+    // Create fits file instance
+    GFits fitsfile = GFits();
+
+    // Write the sky maps to the FITS file
+    m_tsmap.write(fitsfile);
+    for (int i = 0; i < m_mapnames.size(); i++) {
+    	m_maps[i].write(fitsfile);
+    }
+
+    // Set extension name for all maps
+    for (int i = 0; i < m_mapnames.size(); i++) {
+    	fitsfile[i+1]->extname(m_mapnames[i]);
+    }
+
+	// Add computation log if not all bins are computed
+    if (m_binmin != -1 || m_binmax != -1) {
+    	m_statusmap.write(fitsfile);
+    	fitsfile[m_mapnames.size()+1]->extname("STATUS MAP");
+    }
+
+    // Save FITS file
+    fitsfile.saveto(m_outfile,true);
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void cttsmap::init_members(void)
+{
+    // Initialise members
+    m_infile.clear();
+    m_caldb.clear();
+    m_irf.clear();
+    m_outfile.clear();
+    m_proj.clear();
+    m_coordsys.clear();
+    m_xref     = 0.0;
+    m_yref     = 0.0;
+    m_binsz    = 0.0;
+    m_nxpix    = 0;
+    m_nypix    = 0;
+
+    // Initialise protected members
+    m_obs.clear();
+    m_read_ahead = false;
+    m_binmin     = -1;
+    m_binmax     = -1;
+    m_logL0      = 0.0;
+    m_tsmap.clear();
+    m_statusmap.clear();
+    m_mapnames.clear();
+    m_maps.clear();
+    m_testsource = NULL;
+
+    // Set logger properties
+    log.date(true);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void cttsmap::copy_members(const cttsmap& app)
+{
+    // Copy attributes
+    m_infile   = app.m_infile;
+    m_caldb    = app.m_caldb;
+    m_irf      = app.m_irf;
+    m_outfile  = app.m_outfile;
+    m_proj     = app.m_proj;
+    m_coordsys = app.m_coordsys;
+    m_xref     = app.m_xref;
+    m_yref     = app.m_yref;
+    m_binsz    = app.m_binsz;
+    m_nxpix    = app.m_nxpix;
+    m_nypix    = app.m_nypix;
+
+    // Copy protected members
+    m_binmin     = app.m_binmin;
+    m_binmax     = app.m_binmax;
+    m_logL0      = app.m_logL0;
+    m_obs        = app.m_obs;
+    m_read_ahead = app.m_read_ahead;
+    m_tsmap      = app.m_tsmap;
+    m_mapnames   = app.m_mapnames;
+    m_maps       = app.m_maps;
+    m_statusmap  = app.m_statusmap;
+
+    // Clone protected members
+    m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void cttsmap::free_members(void)
+{
+    // Write separator into logger
+    if (logTerse()) {
+        log << std::endl;
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * @exception GException::invalid_value
+ *            Test source not found or no RA/DEC parameters found for test
+ *            source.
+ *
+ * Get all task parameters from parameter file or (if required) by querying
+ * the user. Most parameters are only required if no observation exists so
+ * far in the observation container. In this case, a single CTA observation
+ * will be added to the container, using the definition provided in the
+ * parameter file.
+ ***************************************************************************/
+void cttsmap::get_parameters(void)
+{
+    // If there are no observations in container then add a single CTA
+    // observation using the parameters from the parameter file
+    if (m_obs.size() == 0) {
+
+        // Allocate CTA observation
+        GCTAObservation obs;
+
+        // Get event file name
+        std::string filename = (*this)["infile"].filename();
+
+        // Try first to open as FITS file
+        try {
+
+            // Load data
+            obs.load(filename);
+
+            // Get other task parameters
+            m_caldb = (*this)["caldb"].string();
+            m_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(m_caldb)) {
+                caldb.rootdir(m_caldb);
+            }
+            else {
+                caldb.open("cta", m_caldb);
+            }
+
+            // Set reponse
+            obs.response(m_irf, caldb);
+
+            // Append observation to container
+            m_obs.append(obs);
+
+        }
+
+        // ... otherwise try to open as XML file
+        catch (GException::fits_open_error &e) {
+
+            // Load observations from XML file
+            m_obs.load(filename);
+
+            // Check if all observations have response information. If
+            // not, get the calibration database parameters and set
+            // the response properly
+            bool asked_for_response = false;
+            for (int i = 0; i < m_obs.size(); ++i) {
+
+                // Get CTA observation
+                GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
+
+                // Continue only if observation is a CTA observation
+                if (obs != NULL) {
+
+                    // If response is not valid then set response from
+                    // task parameters
+                    if (!obs->has_response()) {
+
+                        // Get calibration parameters
+                        if (!asked_for_response) {
+                            m_caldb = (*this)["caldb"].string();
+                            m_irf   = (*this)["irf"].string();
+                            asked_for_response = true;
+                        }
+
+                        // Set response
+                        GCaldb caldb;
+                        if (gammalib::dir_exists(m_caldb)) {
+                            caldb.rootdir(m_caldb);
+                        }
+                        else {
+                            caldb.open("cta", m_caldb);
+                        }
+                        obs->response(m_irf, caldb);
+
+                    } // endif: there was no response
+                } // endif: observation was a CTA observation
+
+            } // endfor: looped over observations
+
+        } // endcatch: file was an XML file
+
+    } // endif: there was no observation in the container
+
+    // If there is are no models associated with the observations then
+    // load now the model definition
+    if (m_obs.models().size() == 0) {
+
+        // Get models XML filename
+        std::string filename = (*this)["srcmdl"].filename();
+
+        // Setup models for optimizing.
+        m_obs.models(GModels(filename));
+
+    } // endif: no models were associated with observations
+
+    // Get name of test source and check container for this name
+    m_srcname = (*this)["srcname"].string();
+    if (!m_obs.models().contains(m_srcname)) {
+        std::string msg = "Source \""+m_srcname+"\" not found in model "
+                          "container. Please add a source with that name "
+                          "or check for a possible typos.";
+    	throw GException::invalid_value(G_GET_PARAMETERS, msg);
+    }
+
+    // Get Model and store it as protected member
+    m_testsource = m_obs.models()[m_srcname]->clone();
+
+    // Check if model has RA and DEC parameters which are necessary to
+    // continue
+    if (m_testsource != NULL) {
+        if (!m_testsource->has_par("RA") || !m_testsource->has_par("DEC")) {
+            std::string msg = "Source \""+m_srcname+"\" has no \"RA\" and "
+                              "\"DEC\" parameters. Only sources with \"RA\" "
+                              " and \"DEC\" parameters can be used as test "
+                              "sources.";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
+        }
+    }
+
+    // Fix the spatial parameters of the test source
+    if (m_testsource != NULL) {
+        (*m_testsource)["RA"].fix();
+        (*m_testsource)["DEC"].fix();
+    }
+
+    // Get map parameters
+    m_xref     = (*this)["xref"].real();
+    m_yref     = (*this)["yref"].real();
+    m_proj     = (*this)["proj"].string();
+    m_coordsys = (*this)["coordsys"].string();
+    m_binsz    = (*this)["binsz"].real();
+    m_nxpix    = (*this)["nxpix"].integer();
+    m_nypix    = (*this)["nypix"].integer();
+
+    // Get optional splitting parameters
+    m_binmin = (*this)["binmin"].integer();
+    m_binmax = (*this)["binmax"].integer();
+    m_logL0  = (*this)["logL0"].real();
+
+    // Optionally read ahead parameters so that they get correctly
+    // dumped into the log file
+    if (m_read_ahead) {
+        m_outfile = (*this)["outfile"].filename();
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Initialise skymaps
+ *
+ * Initialises skymaps that will contain map information.
+ ***************************************************************************/
+void cttsmap::init_maps(void)
+{
+    // Initialise map information
+    m_tsmap.clear();
+
+    // Create skymap
+    m_tsmap = GSkymap(m_proj, m_coordsys,
+                      m_xref, m_yref, -m_binsz, m_binsz,
+                      m_nxpix, m_nypix, 1);
+
+    // Initialise map information
+	m_statusmap.clear();
+
+	// Create status map
+	m_statusmap = GSkymap(m_proj, m_coordsys,
+					      m_xref, m_yref, -m_binsz, m_binsz,
+					      m_nxpix, m_nypix, 1);
+
+	// Initialise maps of free parameters
+    if (m_testsource != NULL) {
+
+        // Loop over all model parameters
+        for (int i = 0; i < m_testsource->size(); i++) {
+        
+            // Skip fixed parameters
+            if ((*m_testsource)[i].is_fixed()) {
+                continue;
+            }
+
+            // Add sky map for free parameters. Note that push_back will
+            // create a copy.
+            m_maps.push_back(m_tsmap);
+
+            // Store parameter name
+            m_mapnames.push_back((*m_testsource)[i].name());
+        
+        } // endfor: looped over all model parameters
+
+	} // endif: test source was valid
+
+    // Return
+    return;
+}
+/***************************************************************************
+ *                      cttsmap - TS map calculation tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file cttsmap.hpp
+ * @brief TS map calculation tool interface definition
+ * @author Michael Mayer
+ */
+
+#ifndef CTTSMAP_HPP
+#define CTTSMAP_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+
+/* __Definitions _________________________________________________________ */
+#define CTTSMAP_NAME    "cttsmap"
+#define CTTSMAP_VERSION "00-01-00"
+
+
+/***********************************************************************//**
+ * @class cttsmap
+ *
+ * @brief TS map calculation
+ *
+ * This class computes a set of maps from any kind of observations:
+ *
+ *      - TS map
+ *      - Flux map
+ *      - Map of spectral index
+ *
+ * The class operates on predefined observation containers, an individual
+ * event list or an observation definition XML file.
+ *
+ * During the computation a putative point-like source is moved along a
+ * grid of coordinates. The best fit results (TS, flux, index) are stored in
+ * maps which are saved in the output FITS files.
+ ***************************************************************************/
+class cttsmap : public GApplication  {
+public:
+    // Constructors and destructors
+	cttsmap(void);
+    explicit cttsmap(const GObservations& obs);
+    cttsmap(int argc, char *argv[]);
+    cttsmap(const cttsmap& app);
+    virtual ~cttsmap(void);
+
+    // Operators
+    cttsmap& operator=(const cttsmap& app);
+
+    // Methods
+    void                 clear(void);
+    void                 execute(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    const GSkymap&       tsmap(void) const;
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const cttsmap& app);
+    void free_members(void);
+    void get_parameters(void);
+    void init_maps(void);
+
+    // User parameters
+    std::string              m_infile;     //!< Input file
+    std::string              m_srcname;    //!< Name of source which is moved around
+    std::string              m_caldb;      //!< Calibration database
+    std::string              m_irf;        //!< Instrument response functions
+    std::string              m_outfile;    //!< Output counts map or XML file
+    std::string              m_proj;       //!< WCS projection
+    std::string              m_coordsys;   //!< Coordinate system
+    double                   m_xref;       //!< Longitude reference coordinate
+    double                   m_yref;       //!< Latitude reference coordinate
+    double                   m_binsz;      //!< Pixel size
+    int                      m_nxpix;      //!< Number of pixels in longitude
+    int                      m_nypix;      //!< Number of pixels in latitude
+
+    // Parameters to control speed and job splitting
+    int                      m_binmin;  //!< Map bin number from which computation should start
+    int                      m_binmax;  //!< Map bin number where map computation should end
+    double                   m_logL0;   //!< Likelihood value of null hypothesis
+
+    // Protected members
+    GObservations            m_obs;        //!< Observation container
+    bool                     m_read_ahead; //!< Read ahead parameters
+    GSkymap                  m_tsmap;      //!< TS map
+    GSkymap                  m_statusmap;  //!< Map of computed bins
+    std::vector<std::string> m_mapnames;   //!< Names of free parameters
+    std::vector<GSkymap>     m_maps;       //!< Sky maps for each free parameter
+    GModel*                  m_testsource; //!< Pointer to test source for TS computation
+};
+
+
+/***********************************************************************//**
+ * @brief Return observation container
+ *
+ * @return Reference to observation container
+ ***************************************************************************/
+inline
+const GObservations& cttsmap::obs(void) const
+{
+    return m_obs;
+}
+
+
+/***********************************************************************//**
+ * @brief Return TS skymap
+ *
+ * @return Reference to TS skymap
+ ***************************************************************************/
+inline
+const GSkymap& cttsmap::tsmap(void) const
+{
+    return m_tsmap;
+}
+
+#endif /* CTTSMAP_HPP */
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      cttsmap.par                                                   #
+#  Version:   00-01-00                                                      #
+#                                                                           #
+#  Author:    Michael Mayer                                                 #
+#             michael.mayer@desy.de                                         #
+#             DESY                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the cttsmap application                    #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+infile,   f, a, "events.fits",,,"Input event list, cube or observation definition file"
+srcmdl,   f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+srcname,  s, a, "TestSource",,,"Test source"
+caldb,    s, a, "dummy",,, "Calibration database"
+irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
+outfile,  f, a, "tsmap.fits",,, "Output counts cube"
+
+#
+# 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)"
+proj,     s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"
+
+#
+# Parameters for splitting and speed purpose
+#===========================
+binmin,   i, h,   -1,,, "First bin to compute"
+binmax,   i, h,   -1,,, "Last bin to compute"
+logL0,    r, h,  0.0,,, "LogLikelihood value of null hypothesis"
+
+
+#
+# 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, "cttsmap.log",,, "Log filename"
+/***************************************************************************
+ *                      cttsmap - TS map calculation tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file cttsmap/main.cpp
+ * @brief CTA TS map generation
+ * @author Michael Mayer
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "cttsmap.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 cttsmap application. It allocates a
+ * cttsmap object and runs the application.
+ ***************************************************************************/
+int main (int argc, char *argv[])
+{
+    // Initialise return code
+    int rc = 1;
+
+    // Create instance of application
+    cttsmap 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"
+                              " cttsmap. 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;
+}
@@ -36,6 +36,7 @@ $ECHO &quot;***************&quot;
 # ==================
 ctbin=../src/ctbin/ctbin
 ctlike=../src/ctlike/ctlike
+cttsmap=../src/cttsmap/cttsmap
 ctmodel=../src/ctmodel/ctmodel
 ctobssim=../src/ctobssim/ctobssim
 ctselect=../src/ctselect/ctselect
@@ -347,6 +348,34 @@ $ECHO &quot; ok&quot;
 #
+# Test cttsmap
+# ============
+$ECHO -n "Test cttsmap: "
+$cttsmap infile="selected_events.fits" \
+         srcmdl="data/crab.xml" \
+         srcname="Crab" \
+         caldb="irf" \
+         irf="cta_dummy_irf" \
+         outfile="tsmap.fits" \
+         nxpix=5 \
+         nypix=5 \
+         binsz=0.05 \
+         coordsys="CEL" \
+         xref=83.63 \
+         yref=22.01 \
+         proj="CAR"
+$ECHO -n "."
+if [ -s "tsmap.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " tsmap.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+
+
+#
 # Test ctexpcube
 # ==============
 $ECHO -n "Test ctexpcube: "
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cttsmap 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import ctools
+
+
+# ========================== #
+# Test class for cttsmap tool #
+# ========================== #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for cttsmap tool.
+    """
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Set members
+        self.events_name = "data/crab_events.fits"
+        self.model_name  = "data/crab.xml"
+        self.caldb       = "irf"
+        self.irf         = "cta_dummy_irf"
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name("cttsmap")
+
+        # Append tests
+        self.append(self.test_functional, "Test cttsmap functionality")
+
+        # Return
+        return
+
+    # Test cttsmap functionnality
+    def test_functional(self):
+        """
+        Test cttsmap functionnality.
+        """
+        # Set-up cttsmap
+        tsmap = ctools.cttsmap()
+        tsmap["infile"].filename(self.events_name)
+        tsmap["srcmdl"].filename(self.model_name)
+        tsmap["srcname"].string("Crab")
+        tsmap["caldb"].string(self.caldb)
+        tsmap["irf"].string(self.irf)
+        tsmap["outfile"].filename("tsmap.fits")
+        tsmap["nxpix"].integer(5)
+        tsmap["nypix"].integer(5)
+        tsmap["binsz"].real(0.02)
+        tsmap["coordsys"].string("CEL")
+        tsmap["proj"].string("CAR")
+        tsmap["xref"].real(83.63)
+        tsmap["yref"].real(22.01)
+
+        # Run tool
+        self.test_try("Run cttsmap")
+        try:
+            tsmap.run()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in cttsmap.")
+
+        # Save results
+        self.test_try("Save results")
+        try:
+            tsmap.save()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in saving results.")
+        
+        # Return
+        return
@@ -26,6 +26,7 @@ import test_ctobssim
 import test_ctselect
 import test_ctbin
 import test_ctlike
+import test_cttsmap
 import test_ctmodel
 import test_ctskymap
 import test_ctexpcube
@@ -50,6 +51,7 @@ if __name__ == &#39;__main__&#39;:
     suite_ctselect   = test_ctselect.Test()
     suite_ctbin      = test_ctbin.Test()
     suite_ctlike     = test_ctlike.Test()
+    suite_cttsmap    = test_cttsmap.Test()
     suite_ctmodel    = test_ctmodel.Test()
     suite_ctskymap   = test_ctskymap.Test()
     suite_ctexpcube  = test_ctexpcube.Test()
@@ -63,6 +65,7 @@ if __name__ == &#39;__main__&#39;:
     suite_ctselect.set()
     suite_ctbin.set()
     suite_ctlike.set()
+    suite_cttsmap.set()
     suite_ctmodel.set()
     suite_ctskymap.set()
     suite_ctexpcube.set()
@@ -76,6 +79,7 @@ if __name__ == &#39;__main__&#39;:
     suites.append(suite_ctselect)
     suites.append(suite_ctbin)
     suites.append(suite_ctlike)
+    suites.append(suite_cttsmap)
     suites.append(suite_ctmodel)
     suites.append(suite_ctskymap)
     suites.append(suite_ctexpcube)