/***************************************************************************
* ctmodel - Model cube generation tool *
* ----------------------------------------------------------------------- *
* copyright (C) 2012-2017 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 ctmodel.cpp
* @brief Model cube generation tool implementation
* @author Juergen Knoedlseder
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include
#endif
#include
#include "ctmodel.hpp"
#include "GTools.hpp"
/* __ Method name definitions ____________________________________________ */
#define G_GET_PARAMETERS "ctmodel::get_parameters()"
#define G_FILL_CUBE "ctmodel::fill_cube(GCTAObservation*)"
/* __ Debug definitions __________________________________________________ */
/* __ Coding definitions _________________________________________________ */
/* __ Constants __________________________________________________________ */
const GEnergy g_energy_margin(1.0e-12, "TeV");
/*==========================================================================
= =
= Constructors/destructors =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Void constructor
***************************************************************************/
ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, VERSION)
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Observations constructor
*
* param[in] obs Observation container.
*
* This method creates an instance of the class by copying an existing
* observations container.
***************************************************************************/
ctmodel::ctmodel(const GObservations& obs) :
ctobservation(CTMODEL_NAME, VERSION, obs)
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Command line constructor
*
* @param[in] argc Number of arguments in command line.
* @param[in] argv Array of command line arguments.
***************************************************************************/
ctmodel::ctmodel(int argc, char *argv[]) :
ctobservation(CTMODEL_NAME, VERSION, argc, argv)
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Copy constructor
*
* @param[in] app Application.
***************************************************************************/
ctmodel::ctmodel(const ctmodel& app) : ctobservation(app)
{
// Initialise members
init_members();
// Copy members
copy_members(app);
// Return
return;
}
/***********************************************************************//**
* @brief Destructor
***************************************************************************/
ctmodel::~ctmodel(void)
{
// Free members
free_members();
// Return
return;
}
/*==========================================================================
= =
= Operators =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Assignment operator
*
* @param[in] app Application.
* @return Application.
***************************************************************************/
ctmodel& ctmodel::operator=(const ctmodel& app)
{
// Execute only if object is not identical
if (this != &app) {
// Copy base class members
this->ctobservation::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 ctmodel tool
*
* Clears ctmodel tool.
***************************************************************************/
void ctmodel::clear(void)
{
// Free members
free_members();
this->ctobservation::free_members();
this->ctool::free_members();
// Clear base class (needed to conserve tool name and version)
this->GApplication::clear();
// Initialise members
this->ctool::init_members();
this->ctobservation::init_members();
init_members();
// Write header into logger
log_header();
// 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 model map for each CTA observation.
***************************************************************************/
void ctmodel::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();
// Extract cube properties
extract_cube_properties();
// Set energy dispersion flags of all CTA observations and save old
// values in save_edisp vector
std::vector save_edisp = set_edisp(m_obs, m_apply_edisp);
// Write input observation container into logger
log_observations(NORMAL, m_obs, "Input observation");
// Write input model container into logger
log_models(NORMAL, m_obs.models(), "Input model");
// Write header
log_header1(TERSE, "Generate model cube");
// Loop over all observations in the container
#pragma omp parallel
{
// Create a models object for this core
GModels models(m_obs.models());
#pragma omp for schedule(static,1)
for (int i = 0; i < m_obs.size(); ++i) {
// Get CTA observation
GCTAObservation* obs = dynamic_cast(m_obs[i]);
// Skip observation if it's not CTA
if (obs == NULL) {
log_header3(TERSE, get_obs_header(m_obs[i]));
std::string msg = " Skipping "+m_obs[i]->instrument()+
" observation";
log_string(NORMAL, msg);
}
// Fill cube and leave loop if we are binned mode (meaning we
// only have one binned observation)
else if (m_binned) {
fill_cube(obs, models);
i = m_obs.size();
}
// Skip observation if we have a binned observation
else if (obs->eventtype() == "CountsCube") {
log_header3(TERSE, get_obs_header(m_obs[i]));
std::string msg = " Skipping binned "+m_obs[i]->instrument()+
" observation";
log_string(NORMAL, msg);
}
// Otherwise, everything seems to be fine, so fill the cube from obs
else {
// Fill the cube
fill_cube(obs, models);
// Dispose events to free memory if event file exists on disk
if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
obs->dispose_events();
}
}
} // endfor: looped over observations
} // end openmp parallel section
// Write model cube into header
log_header1(NORMAL, "Model cube");
log_string(NORMAL, m_cube.print());
// Restore energy dispersion flags of all CTA observations
restore_edisp(m_obs, save_edisp);
// Optionally publish model cube
if (m_publish) {
publish();
}
// Return
return;
}
/***********************************************************************//**
* @brief Save model cube
*
* Saves the model cube into a FITS file specified using the "outfile"
* task parameter.
***************************************************************************/
void ctmodel::save(void)
{
// Write header
log_header1(TERSE, "Save model cube");
// Get model cube filename
m_outcube = (*this)["outcube"].filename();
// Save only if filename is non-empty and cube has some size
if (!m_outcube.is_empty() && m_cube.size() > 0) {
m_cube.save(m_outcube, clobber());
}
// Write into logger what has been done
std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url();
if (m_cube.size() == 0) {
fname.append(" (cube is empty, no file created)");
}
log_value(NORMAL, "Model cube file", fname);
// Return
return;
}
/***********************************************************************//**
* @brief Publish model cube
*
* @param[in] name Model cube name.
***************************************************************************/
void ctmodel::publish(const std::string& name)
{
// Write header into logger
log_header1(TERSE, "Publish model cube");
// Set default name is user name is empty
std::string user_name(name);
if (user_name.empty()) {
user_name = CTMODEL_NAME;
}
// Log filename
log_value(NORMAL, "Model cube name", user_name);
// Publish model cube
m_cube.counts().publish(user_name);
// Return
return;
}
/***********************************************************************//**
* @brief Set model cube
*
* @param[in] cube Model cube.
*
* Set model cube and set all cube bins to zero.
***************************************************************************/
void ctmodel::cube(const GCTAEventCube& cube)
{
// Set cube
m_cube = cube;
// Set all cube bins to zero
for (int i = 0; i < m_cube.size(); ++i) {
m_cube[i]->counts(0.0);
}
// Signal that cube has been set
m_has_cube = true;
// Return
return;
}
/*==========================================================================
= =
= Private methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Initialise class members
***************************************************************************/
void ctmodel::init_members(void)
{
// Initialise members
m_outcube.clear();
m_apply_edisp = false;
m_publish = false;
m_chatter = static_cast(2);
// Initialise protected members
m_cube.clear();
m_gti.clear();
m_has_cube = false;
m_append_cube = false;
m_binned = false;
m_dir.clear();
m_solidangle.clear();
m_energy.clear();
m_ewidth.clear();
m_time.clear();
// Return
return;
}
/***********************************************************************//**
* @brief Copy class members
*
* @param[in] app Application.
***************************************************************************/
void ctmodel::copy_members(const ctmodel& app)
{
// Copy attributes
m_outcube = app.m_outcube;
m_apply_edisp = app.m_apply_edisp;
m_publish = app.m_publish;
m_chatter = app.m_chatter;
// Copy protected members
m_cube = app.m_cube;
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_dir = app.m_dir;
m_solidangle = app.m_solidangle;
m_energy = app.m_energy;
m_ewidth = app.m_ewidth;
m_time = app.m_time;
// Return
return;
}
/***********************************************************************//**
* @brief Delete class members
***************************************************************************/
void ctmodel::free_members(void)
{
// 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.
***************************************************************************/
void ctmodel::get_parameters(void)
{
// Reset cube append flag
m_append_cube = false;
// If there are no observations in container then load them via user
// parameters.
if (m_obs.size() == 0) {
get_obs();
}
// ... otherwise add response information and energy boundaries in case
// that they are missing
else {
setup_observations(m_obs);
}
// If we have now excactly one CTA observation (but no cube has yet been
// appended to the observation) then check whether this observation
// is a binned observation, and if yes, extract the counts cube for
// model generation
if ((m_obs.size() == 1) && (m_append_cube == false)) {
// Get CTA observation
GCTAObservation* obs = dynamic_cast(m_obs[0]);
// Continue only if observation is a CTA observation
if (obs != NULL) {
// Check for binned observation
if (obs->eventtype() == "CountsCube") {
// Set cube from binned observation
GCTAEventCube* evtcube = dynamic_cast
(const_cast(obs->events()));
cube(*evtcube);
// Signal that cube has been set
m_has_cube = true;
// Signal that we are in binned mode
m_binned = true;
} // endif: observation was binned
} // endif: observation was CTA
} // endif: had exactly one observation
// Read model definition file if required
if (m_obs.models().size() == 0) {
// Get model filename
std::string inmodel = (*this)["inmodel"].filename();
// Load models from file
m_obs.models(inmodel);
} // endif: there were no models
// Get energy dispersion flag parameters
m_apply_edisp = (*this)["edisp"].boolean();
// If we do not have yet a counts cube for model computation then check
// whether we should read it from the "incube" parameter or whether we
// should create it from scratch using the task parameters
if (!m_has_cube) {
// Read cube definition file
std::string incube = (*this)["incube"].filename();
// If the cube filename is valid the load the cube and set all cube
// bins to zero
if (is_valid_filename(incube)) {
// Load cube from given file
m_cube.load(incube);
// Set all cube bins to zero
for (int i = 0; i < m_cube.size(); ++i) {
m_cube[i]->counts(0.0);
}
} // endif: cube filename was valid
// ... otherwise create a cube from the user parameters
else {
m_cube = create_cube(m_obs);
}
// Signal that cube has been set
m_has_cube = true;
} // endif: we had no cube yet
// Get remaining parameters
m_publish = (*this)["publish"].boolean();
m_chatter = static_cast((*this)["chatter"].integer());
// Read optionally output cube filenames
if (read_ahead()) {
m_outcube = (*this)["outcube"].filename();
}
// If cube should be appended to first observation then do that now.
// This is a kluge that makes sure that the cube is passed as part
// of the observation in case that a cube response is used. The kluge
// is needed because the GCTACubeSourceDiffuse::set method needs to
// get the full event cube from the observation. It is also at this
// step that the GTI, which may just be a dummy GTI when create_cube()
// has been used, will be set.
if (m_append_cube) {
//TODO: Check that energy boundaries are compatible
// Attach GTI of observations to model cube
m_cube.gti(m_obs[0]->events()->gti());
// Attach model cube to observations
m_obs[0]->events(m_cube);
} // endif: cube was scheduled for appending
// Write parameters into logger
log_parameters(TERSE);
// Return
return;
}
/***********************************************************************//**
* @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.
*
* The method sets m_append_cube = true and m_binned = true in case that
* a stacked observation is requested (as detected by the presence of the
* "expcube", "psfcube", and "bkgcube" parameters). In that case, it appended
* a dummy event cube to the observation.
*
* @todo Support stacked energy dispersion
***************************************************************************/
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 (!is_valid_filename(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 (is_valid_filename(expcube) && is_valid_filename(psfcube) &&
is_valid_filename(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());
// 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);
// Signal that we are in binned mode
m_binned = true;
// Signal that we appended a cube
m_append_cube = true;
} // 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 not valid
// ... 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 (GFilename(filename).is_fits()) {
// 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 Extract cube properties in data members
***************************************************************************/
void ctmodel::extract_cube_properties(void)
{
// Clear cube properties
m_dir.clear();
m_solidangle.clear();
m_energy.clear();
m_ewidth.clear();
// Get number of spatial pixels and energy layers in counts cube
int npix = m_cube.npix();
int nebins = m_cube.ebins();
// Reserve capacity in vectors
m_dir.reserve(npix);
m_solidangle.reserve(npix);
m_energy.reserve(nebins);
m_ewidth.reserve(nebins);
// Loop over all the spatial bins and extract spatial cube properties
for (int i = 0; i < npix; ++i) {
// Get cube bin
GCTAEventBin* bin = m_cube[i];
// Extract sky direction and solid angle
m_dir.push_back(bin->dir());
m_solidangle.push_back(bin->solidangle());
} // endfor: looped over spatial bins
// Loop over all the spectral bins and extract spectral cube properties
for (int iebin = 0, ibin = 0; iebin < nebins; ++iebin, ibin += npix) {
// Get cube bin
GCTAEventBin* bin = m_cube[ibin];
// Extract energy and energy width
m_energy.push_back(bin->energy());
m_ewidth.push_back(bin->ewidth());
} // endfor: looped over spectral bins
// Set time
m_time = m_cube[0]->time();
// Return
return;
}
/***********************************************************************//**
* @brief Fill model into model cube
*
* @param[in] obs CTA observation.
*
* Adds the expected number of events for a given observation to the events
* that are already found in the model cube. The method also updates the
* GTI of the model cube so that cube GTI is a list of the GTIs of all
* observations that were used to generate the model cube.
***************************************************************************/
void ctmodel::fill_cube(const GCTAObservation* obs, GModels& models)
{
// Get number of spatial and spectral bins in counts cube
int npix = m_cube.npix();
int nebins = m_cube.ebins();
// Get references to GTI and energy boundaries for the event list
const GGti& gti = obs->events()->gti();
const GEbounds& obs_ebounds = obs->ebounds();
// Get cube energy boundaries
const GEbounds& cube_ebounds = m_cube.ebounds();
// Get counts cube usage flags
std::vector usage = cube_layer_usage(cube_ebounds, obs_ebounds);
// Initialise empty, invalid RoI
GCTARoi roi;
// Retrieve RoI in case we have an unbinned observation
if (obs->eventtype() == "EventList") {
roi = obs->roi();
}
// Initialise statistics
double sum = 0.0;
int num_outside_ebds = 0;
int num_outside_roi = 0;
// Extract copy of models from observation container
GModels ex_models = trim_models(models, roi);
// Get pointer to event cube pixels
double* pixels = const_cast(m_cube.counts().pixels());
// Initialise event bin
GCTAEventBin bin;
// Set event bin attributes that are constant for a counts cube
bin.time(m_time);
bin.ontime(obs->ontime());
bin.weight(1.0);
// Loop over all the spatial bins
for (int i = 0; i < npix; ++i) {
// If RoI is valid then skip bin if it is outside the RoI of the
// observation
if (roi.is_valid() && !roi.contains(m_dir[i])) {
num_outside_roi += nebins;
continue;
}
// Set instrument direction, solid angle and pixel index of bin
bin.dir(m_dir[i]);
bin.solidangle(m_solidangle[i]);
bin.ipix(i);
// Loop over all of the energy bins of the cube
for (int iebin = 0, ibin = i; iebin < nebins; ++iebin, ibin += npix) {
// Skip bin if the corresponding counts cube energy bin is not fully
// contained in the event list energy range. This avoids having
// partially filled bins.
if (!usage[iebin]) {
num_outside_ebds++;
continue;
}
// Set energy, energy width and energy index of bin
bin.energy(m_energy[iebin]);
bin.ewidth(m_ewidth[iebin]);
bin.ieng(iebin);
// Compute model value for cube bin
double model = models.eval(bin, *obs) * bin.size();
// Sum model
sum += model;
// Store value
#pragma omp critical(ctmodel_fill_cube)
pixels[ibin] += model;
} // endfor: looped over all spatial bins
} // endfor: looped over all cube bins
// Re-append the models to the main model
models.extend(ex_models);
// Update GTIs
#pragma omp critical(ctmodel_fill_cube)
{
// Append GTIs of observation to list of GTIs
m_gti.extend(gti);
// Update GTIs
m_cube.gti(m_gti);
}
// Log filling results
#pragma omp critical(ctmodel_fill_cube)
{
log_header3(TERSE, get_obs_header(obs));
log_value(NORMAL, "Model events in cube", sum);
log_value(NORMAL, "Bins outside energy range", num_outside_ebds);
log_value(NORMAL, "Bins outside RoI", num_outside_roi);
}
// Write model cube into header
if (m_chatter >= EXPLICIT) {
log_header2(EXPLICIT, "Model cube");
log_string(EXPLICIT, m_cube.print(m_chatter));
}
// Return
return;
}
/***********************************************************************//**
* @brief Find the models falling inside a defined region of interest.
*
* @param[in,out] all_models Model container to be trimmed.
* @param[in] roi Observation region of interest.
* @return New model container containing only models inside RoI.
*
* Note that a buffer is added to the observation region to ensure that
* point sources (which have radius=0) are appropriately included if they
* fall near the edge of the observation.
***************************************************************************/
GModels ctmodel::trim_models(GModels& all_models, const GCTARoi& roi)
{
// Create the model container containing only models in RoI
GModels excluded_models;
// Do nothing if roi is not valid. This will be the case for binned data
// sets as they have the default roi radius of 0. In this case we would
// remove ALL models, which is obviously incorrect. So we keep all models
// when filling based on a counts cube.
if (!roi.is_valid()) {
return excluded_models;
}
// Remove all models that dont overlap with the region of interest. Note
// that an extra factor is used since point sources have regions of
// radius 0, which is a problem when they fall just barely outside the
// RoI. This ensures point sources are appropriately included.
GSkyRegionCircle obsreg(roi.centre().dir(), roi.radius()+0.5);
// Loop over the models in the passed container
for (int i = 0; i < all_models.size(); ++i) {
// Cast the model to a GModelSky object
GModelSky* model = dynamic_cast(all_models[i]);
// If model is not a GModelSky model than skip
if (model == NULL) {
continue;
}
// Otherwise, exclude the model if it isn't a diffuse cube and it
// doesn't overlap with the observation.
if (!(model->spatial()->code() == GMODEL_SPATIAL_DIFFUSE) &&
(!model->spatial()->region()->overlaps(obsreg))) {
// Append model to the excluded models
excluded_models.append(*model);
// Remove model
all_models.remove(i);
// Decrement to prevent skipping models
i--;
}
} // endfor: looped over all models
// Return excluded models
return excluded_models;
}