/***************************************************************************
* ctbkgcube - Background cube generation tool *
* ----------------------------------------------------------------------- *
* copyright (C) 2014-2016 by Chia-Chun Lu *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see . *
* *
***************************************************************************/
/**
* @file ctbkgcube.cpp
* @brief Background cube generation tool definition
* @author Chia-Chun Lu
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include
#endif
#include
#include "ctbkgcube.hpp"
#include "GTools.hpp"
/* __ Method name definitions ____________________________________________ */
#define G_GET_PARAMETERS "ctbkgcube::get_parameters()"
/* __ Debug definitions __________________________________________________ */
/* __ Coding definitions _________________________________________________ */
/* __ Constants __________________________________________________________ */
/*==========================================================================
= =
= Constructors/destructors =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Void constructor
***************************************************************************/
ctbkgcube::ctbkgcube(void) : ctool(CTBKGCUBE_NAME, CTBKGCUBE_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.
***************************************************************************/
ctbkgcube::ctbkgcube(const GObservations& obs) :
ctool(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
{
// Initialise members
init_members();
// Set observations
m_obs = obs;
// Return
return;
}
/***********************************************************************//**
* @brief Command line constructor
*
* @param[in] argc Number of arguments in command line.
* @param[in] argv Array of command line arguments.
***************************************************************************/
ctbkgcube::ctbkgcube(int argc, char *argv[]) :
ctool(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, argc, argv)
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Copy constructor
*
* @param[in] app Application.
***************************************************************************/
ctbkgcube::ctbkgcube(const ctbkgcube& app) : ctool(app)
{
// Initialise members
init_members();
// Copy members
copy_members(app);
// Return
return;
}
/***********************************************************************//**
* @brief Destructor
***************************************************************************/
ctbkgcube::~ctbkgcube(void)
{
// Free members
free_members();
// Return
return;
}
/*==========================================================================
= =
= Operators =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Assignment operator
*
* @param[in] app Application.
* @return Returns application.
***************************************************************************/
ctbkgcube& ctbkgcube::operator=(const ctbkgcube& app)
{
// Execute only if object is not identical
if (this != &app) {
// Copy base class members
this->ctool::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 ctbkgcube::clear(void)
{
// Free members
free_members();
this->ctool::free_members();
this->GApplication::free_members();
// Initialise members
this->GApplication::init_members();
this->ctool::init_members();
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Generate the background cube(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 background cube from the CTA observations.
***************************************************************************/
void ctbkgcube::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;
}
// Log input model
if (logTerse()) {
log << std::endl;
log.header1("Input model");
log << m_obs.models() << std::endl;
}
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Prepare model");
}
// Copy models from observation container and reset output model
// container
GModels models_orig = m_obs.models();
m_bkgmdl = m_obs.models();
m_outmdl.clear();
// Remove all models that are not CTA background models from the
// container and put all removed components in the output
// container
int num = m_bkgmdl.size();
for (int i = num-1; i >= 0; --i) {
// Flag removal
bool remove = true;
// Do we have a CTA, HESS, MAGIC or VERITAS specific model?
if (m_bkgmdl[i]->is_valid("CTA", "") ||
m_bkgmdl[i]->is_valid("HESS", "") ||
m_bkgmdl[i]->is_valid("MAGIC", "") ||
m_bkgmdl[i]->is_valid("VERITAS", "")) {
// Do we have a background model?
if (dynamic_cast(m_bkgmdl[i]) != NULL) {
// ... then keep model
remove = false;
} // endif: had a background model
} // endif: had a CTA, HESS, MAGIC or VERITAS model
// Log results
if (logTerse()) {
if (remove) {
log << gammalib::parformat("Remove model");
}
else {
log << gammalib::parformat("Keep model");
}
log << m_bkgmdl[i]->name();
log << " ";
log << m_bkgmdl[i]->type();
log << "(";
log << m_bkgmdl[i]->instruments();
log << ")";
log << std::endl;
}
// If removal is requested, append model to output container and
// remove it from the background model container ...
if (remove) {
m_outmdl.append(*(m_bkgmdl[i]));
m_bkgmdl.remove(i);
}
} // endfor: looped over all background models
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Generate background cube");
}
// Assign background models to container
m_obs.models(m_bkgmdl);
// Fill background cube from observations
m_background.fill(m_obs, &log);
// Create a background model for the output background cube and append
// that model to the input model in place of the original
// background models
// TODO: We might think of creating the spectral model via user parameter
GModelSpectralPlaw spectral(1.0, 0.0, GEnergy(1.0, "TeV"));
spectral["Prefactor"].range(0.01, 100.0);
spectral["Index"].range(-5.0, 5.0);
GCTAModelCubeBackground model(spectral);
// Set model name
model.name("BackgroundModel");
// Set model instrument
// TODO: Account for possibility to have observations from different
// IACTs in the same container
model.instruments("CTA,HESS,MAGIC,VERITAS");
// Append model to output container
m_outmdl.append(model);
// Log output model
if (logTerse()) {
log << std::endl;
log.header1("Output model");
log << m_outmdl << std::endl;
}
// Recover original models
m_obs.models(models_orig);
// Optionally publish background cube
if (m_publish) {
publish();
}
// Return
return;
}
/***********************************************************************//**
* @brief Save background cube
*
* Save the background cube into a FITS file and the output model into an
* XML file.
***************************************************************************/
void ctbkgcube::save(void)
{
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Save background cube");
}
// Get background cube and model file names
m_outcube = (*this)["outcube"].filename();
m_outmodel = (*this)["outmodel"].filename();
// Save only if filename is non-empty
if (!m_outcube.is_empty()) {
// Log filename
if (logTerse()) {
log << "Save background cube into file \""+m_outcube+"\".";
log << std::endl;
}
// Save background cube
m_background.save(m_outcube, clobber());
}
// Write output models if filename is valid
if ((!m_outmodel.is_empty()) &&
(gammalib::tolower(m_outmodel.url()) != "none")) {
// Log filename
if (logTerse()) {
log << "Save model into file \""+m_outmodel.url()+"\".";
log << std::endl;
}
// Save output model for stacked analyses
m_outmdl.save(m_outmodel.url());
}
// Return
return;
}
/***********************************************************************//**
* @brief Publish background cube
*
* @param[in] name Background cube name.
***************************************************************************/
void ctbkgcube::publish(const std::string& name)
{
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Publish background cube");
}
// Set default name is user name is empty
std::string user_name(name);
if (user_name.empty()) {
user_name = CTBKGCUBE_NAME;
}
// Log filename
if (logTerse()) {
log << "Publish \""+user_name+"\" background cube." << std::endl;
}
// Publish exposure cube
m_background.cube().publish(user_name);
// Return
return;
}
/*==========================================================================
= =
= Private methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Initialise class members
***************************************************************************/
void ctbkgcube::init_members(void)
{
// Initialise members
m_outcube.clear();
m_outmodel.clear();
m_publish = false;
m_obs.clear();
m_background.clear();
m_bkgmdl.clear();
m_outmdl.clear();
// Return
return;
}
/***********************************************************************//**
* @brief Copy class members
*
* @param[in] app Application.
***************************************************************************/
void ctbkgcube::copy_members(const ctbkgcube& app)
{
// Copy members
m_outmodel = app.m_outmodel;
m_outcube = app.m_outcube;
m_publish = app.m_publish;
m_obs = app.m_obs;
m_background = app.m_background;
m_bkgmdl = app.m_bkgmdl;
m_outmdl = app.m_outmdl;
// Return
return;
}
/***********************************************************************//**
* @brief Delete class members
***************************************************************************/
void ctbkgcube::free_members(void)
{
// Return
return;
}
/***********************************************************************//**
* @brief Get application parameters
*
* @exception GException::invalid_value
* No background model definition XML file specified.
*
* Get all task parameters from parameter file or (if required) by querying
* the user. The parameters are read in the correct order.
***************************************************************************/
void ctbkgcube::get_parameters(void)
{
// If there are no observations in container then load them via user
// parameters
if (m_obs.size() == 0) {
// Throw exception if no input observation file is given
require_inobs(G_GET_PARAMETERS);
// Throw exception if counts cube is given
require_inobs_nocube(G_GET_PARAMETERS);
// Build observation container
m_obs = get_observations();
} // endif: there was no observation in the container
// Get the incube filename
std::string incube = (*this)["incube"].filename();
// Check for filename validity
if ((gammalib::toupper(incube) == "NONE") ||
(gammalib::strip_whitespace(incube) == "")) {
// Create an event cube based on task parameters
GCTAEventCube cube = create_cube(m_obs);
// Define background cube
m_background = GCTACubeBackground(cube);
} // endif: filename was not valid
// ... otherwise setup the background cube from the counts map
else {
// Load event cube from filename
GCTAEventCube cube(incube);
// Define background cube from file
m_background = GCTACubeBackground(cube);
} // endelse: cube was loaded from file
// If there are no models associated with the observations then load now
// the model definition from the XML file
if (m_obs.models().size() == 0) {
if ((*this)["inmodel"].is_undefined()) {
std::string msg = "No model definition XML file specified. "
"Please set the \"inmodel\" parameter to the "
"XML file that contains the background model "
"definition.";
throw GException::invalid_value(G_GET_PARAMETERS, msg);
}
std::string inmodel = (*this)["inmodel"].filename();
GModels models(inmodel);
m_obs.models(models);
}
// Get remaining parameters
m_publish = (*this)["publish"].boolean();
// Read output filenames (if needed)
if (read_ahead()) {
m_outcube = (*this)["outcube"].filename();
m_outmodel = (*this)["outmodel"].filename();
}
// Return
return;
}