From 1ae952bd36fc6a00890fa42dea5033ee7fd5f9e1 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 20 Dec 2015 20:44:08 +0100 Subject: [PATCH] Re-enable ctmodel model computation for stacked analysis. --- ChangeLog | 2 ++ NEWS | 3 ++- src/ctmodel/ctmodel.cpp | 134 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- src/ctmodel/ctmodel.hpp | 3 ++- src/ctmodel/ctmodel.par | 2 +- 5 files changed, 138 insertions(+), 6 deletions(-) diff --git a/ChangeLog b/ChangeLog index 3f48893..a3e0a10 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,8 @@ * Bug fix version 1.0.1 released ================================ + Re-enable ctmodel computation for stacked analysis (#1609) + 2015-12-18 Juergen Knoedlseder diff --git a/NEWS b/NEWS index 4da5b00..696dbc1 100644 --- a/NEWS +++ b/NEWS @@ -59,7 +59,8 @@ None 11. ctmodel - CTA counts cube modelling --------------------------------------- -None +The get_obs() method, that has been accidently removed, has been put back +which allows computation of a model cube for a stacked analysis. 12. ctexpcube - CTA exposure cube generation diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 9aa8bda..8351883 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -411,7 +411,7 @@ void ctmodel::init_members(void) m_gti.clear(); m_has_cube = false; m_append_cube = false; - m_binned = false; + m_binned = false; // Return return; @@ -435,7 +435,7 @@ void ctmodel::copy_members(const ctmodel& app) m_gti = app.m_gti; m_has_cube = app.m_has_cube; m_append_cube = app.m_append_cube; - m_binned = app.m_binned; + m_binned = app.m_binned; // Return return; @@ -466,7 +466,7 @@ void ctmodel::get_parameters(void) // If there are no observations in container then load them via user // parameters if (m_obs.size() == 0) { - m_obs = get_observations(); + get_obs(); } // Check if we got excactly one binned CTA observation @@ -574,6 +574,134 @@ void ctmodel::get_parameters(void) /***********************************************************************//** + * @brief Get observation container + * + * Get an observation container according to the user parameters. The method + * supports loading of a individual FITS file or an observation definition + * file in XML format. + * + * If the input filename is empty, the method checks for the existence of the + * "expcube", "psfcube" and "bkgcube" parameters. If file names have been + * specified, the method loads the files and creates a dummy events cube that + * is appended to the observation container. + * + * If no file names are specified for the "expcube", "psfcube" or "bkgcube" + * parameters, the method reads the necessary parameters to build a CTA + * observation from scratch. + ***************************************************************************/ +void ctmodel::get_obs(void) +{ + // Get the filename from the input parameters + std::string filename = (*this)["inobs"].filename(); + + // If no observation definition file has been specified then read all + // parameters that are necessary to create an observation from scratch + if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) { + + // Get response cube filenames + std::string expcube = (*this)["expcube"].filename(); + std::string psfcube = (*this)["psfcube"].filename(); + std::string bkgcube = (*this)["bkgcube"].filename(); + + // If the filenames are valid then build an observation from cube + // response information + if ((expcube != "NONE") && (psfcube != "NONE") && (bkgcube != "NONE") && + (gammalib::strip_whitespace(expcube) != "") && + (gammalib::strip_whitespace(psfcube) != "") && + (gammalib::strip_whitespace(bkgcube) != "")) { + + // Get exposure, PSF and background cubes + GCTACubeExposure exposure(expcube); + GCTACubePsf psf(psfcube); + GCTACubeBackground background(bkgcube); + + // Create energy boundaries + GEbounds ebounds = create_ebounds(); + + // Create dummy sky map cube + GSkyMap map("CAR","GAL",0.0,0.0,1.0,1.0,1,1,ebounds.size()); + m_append_cube = true; + + // Create event cube + GCTAEventCube cube(map, ebounds, exposure.gti()); + + // Create CTA observation + GCTAObservation cta; + cta.events(cube); + cta.response(exposure, psf, background); + + // Append observation to container + m_obs.append(cta); + + } // endif: cube response information was available + + // ... otherwise build an observation from IRF response information + else { + + // Create CTA observation + GCTAObservation cta = create_cta_obs(); + + // Set response + set_obs_response(&cta); + + // Append observation to container + m_obs.append(cta); + + } + + } // endif: filename was "NONE" or "" + + // ... otherwise we have a file name + else { + + // If file is a FITS file then create an empty CTA observation + // and load file into observation + if (gammalib::is_fits(filename)) { + + // Allocate empty CTA observation + GCTAObservation cta; + + // Load data + cta.load(filename); + + // Set response + set_obs_response(&cta); + + // Append observation to container + m_obs.append(cta); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } + + // ... otherwise load file into observation container + else { + + // Load observations from XML file + m_obs.load(filename); + + // For all observations that have no response, set the response + // from the task parameters + set_response(m_obs); + + // Set observation boundary parameters (emin, emax, rad) + set_obs_bounds(m_obs); + + // Signal that XML file should be used for storage + m_use_xml = true; + + } // endelse: file was an XML file + + } + + // Return + return; + +} + + +/***********************************************************************//** * @brief Fill model into model cube * * @param[in] obs CTA observation. diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index 39a30a8..17bdd5d 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -34,7 +34,7 @@ /* __Definitions _________________________________________________________ */ #define CTMODEL_NAME "ctmodel" -#define CTMODEL_VERSION "1.0.0" +#define CTMODEL_VERSION "1.0.1" /***********************************************************************//** @@ -79,6 +79,7 @@ protected: void copy_members(const ctmodel& app); void free_members(void); void get_parameters(void); + void get_obs(void); void fill_cube(const GCTAObservation* obs); bool has_cube(void) const; diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par index 8534e95..39ae83b 100644 --- a/src/ctmodel/ctmodel.par +++ b/src/ctmodel/ctmodel.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctmodel.par # -# Version: 1.0.0 # +# Version: 1.0.1 # # # # Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # -- libgit2 0.21.4