/***************************************************************************
* ctmapcube - Map cube generation tool *
* ----------------------------------------------------------------------- *
* copyright (C) 2016 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 ctmapcube.cpp
* @brief Map cube generation tool implementation
* @author Juergen Knoedlseder
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include
#endif
#include
#include "ctmapcube.hpp"
#include "GTools.hpp"
/* __ Method name definitions ____________________________________________ */
#define G_GET_PARAMETERS "ctmapcube::get_parameters()"
/* __ Debug definitions __________________________________________________ */
/* __ Coding definitions _________________________________________________ */
/* __ Constants __________________________________________________________ */
/*==========================================================================
= =
= Constructors/destructors =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Void constructor
***************************************************************************/
ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION)
{
// 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.
***************************************************************************/
ctmapcube::ctmapcube(int argc, char *argv[]) :
ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION, argc, argv)
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Copy constructor
*
* @param[in] app Application.
***************************************************************************/
ctmapcube::ctmapcube(const ctmapcube& app) : ctool(app)
{
// Initialise members
init_members();
// Copy members
copy_members(app);
// Return
return;
}
/***********************************************************************//**
* @brief Destructor
***************************************************************************/
ctmapcube::~ctmapcube(void)
{
// Free members
free_members();
// Return
return;
}
/*==========================================================================
= =
= Operators =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Assignment operator
*
* @param[in] app Application.
* @return Application.
***************************************************************************/
ctmapcube& ctmapcube::operator=(const ctmapcube& 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 ctmapcube::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 map cube.
*
* This method reads the task parameters from the parfile and generates the
* map cube.
***************************************************************************/
void ctmapcube::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 models into logger
if (logTerse()) {
log << std::endl;
log.header1(gammalib::number("Model", m_models.size()));
log << m_models << std::endl;
}
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Generate map cube");
}
//TODO
// Log cube
if (logTerse()) {
log << std::endl;
log.header1("Map cube");
log << m_cube << std::endl;
}
// Optionally publish map cube
if (m_publish) {
publish();
}
// Return
return;
}
/***********************************************************************//**
* @brief Save map cube
*
* Saves the map cube into a FITS file.
***************************************************************************/
void ctmapcube::save(void)
{
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Save map cube");
}
// Get map cube filename
m_outcube = (*this)["outcube"].filename();
// Save only if filename is non-empty
if (!m_outcube.is_empty()) {
// Log filename
if (logTerse()) {
log << "Save map cube into file \""+m_outcube+"\".";
log << std::endl;
}
// Save map cube into FITS file
m_cube.save(m_outcube, clobber());
}
// Return
return;
}
/***********************************************************************//**
* @brief Publish map cube
*
* @param[in] name Map cube name.
***************************************************************************/
void ctmapcube::publish(const std::string& name)
{
// Write header
if (logTerse()) {
log << std::endl;
log.header1("Publish map cube");
}
// Set default name is user name is empty
std::string user_name(name);
if (user_name.empty()) {
user_name = CTMAPCUBE_NAME;
}
// Log filename
if (logTerse()) {
log << "Publish \""+user_name+"\" map cube." << std::endl;
}
// Publish map cube
m_cube.cube().publish(user_name);
// Return
return;
}
/*==========================================================================
= =
= Private methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Initialise class members
***************************************************************************/
void ctmapcube::init_members(void)
{
// Initialise members
m_outcube.clear();
m_publish = false;
// Initialise protected members
m_models.clear();
m_cube.clear();
// Return
return;
}
/***********************************************************************//**
* @brief Copy class members
*
* @param[in] app Application.
***************************************************************************/
void ctmapcube::copy_members(const ctmapcube& app)
{
// Copy attributes
m_outcube = app.m_outcube;
m_publish = app.m_publish;
// Copy protected members
m_models = app.m_models;
m_cube = app.m_cube;
// Return
return;
}
/***********************************************************************//**
* @brief Delete class members
***************************************************************************/
void ctmapcube::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 ctmapcube::get_parameters(void)
{
// Get energy binning parameters
GEnergies energies = create_energies();
// Get spatial binning parameters
double xref = (*this)["xref"].real();
double yref = (*this)["yref"].real();
std::string proj = (*this)["proj"].string();
std::string coordsys = (*this)["coordsys"].string();
double binsz = (*this)["binsz"].real();
int nxpix = (*this)["nxpix"].integer();
int nypix = (*this)["nypix"].integer();
// Allocate map cube
GSkyMap cube(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix,
energies.size());
// Allocate map cube
m_cube = GModelSpatialDiffuseCube(cube, energies);
// Read model definition file if required
if (m_models.size() == 0) {
// Get model filename
GFilename inmodel = (*this)["inmodel"].filename();
// Load models from file
m_models.load(inmodel);
} // endif: there were no models
// Get remaining parameters
m_publish = (*this)["publish"].boolean();
// Read optionally output cube filenames
if (read_ahead()) {
m_outcube = (*this)["outcube"].filename();
}
// Write parameters into logger
if (logTerse()) {
log_parameters();
log << std::endl;
}
// Return
return;
}
/***********************************************************************//**
* @brief Create energy vector from user parameters
*
* @exception GException::invalid_value
* No valid energy boundary extension found.
*
* Get the energy vector according to the user parameters. The method
* supports loading of energy information from the ENERGIES extension (or
* any other extension specified by the "ebinfile" parameters), or setting
* energies using a linear or logarithmical spacing.
*
* The following parameters are read:
* ebinalg - Energy binning algorithm
* emin - Minimum energy (if ebinalg != FILE)
* emax - Maximum energy (if ebinalg != FILE)
* enumbins - Number of energies (if ebinalg != FILE)
***************************************************************************/
GEnergies ctmapcube::create_energies(void)
{
// Allocate energies
GEnergies energies;
// Get energy binning algorithm
std::string ebinalg = (*this)["ebinalg"].string();
// If energy binning algorithm is of type "FILE" (case sensitive), then
// read energies from FITS file ...
if (ebinalg == "FILE") {
// Get filename
GFilename ebinfile = (*this)["ebinfile"].filename();
// Load energies
energies.load(ebinfile);
} // endif: ebinalg was "FILE"
// ... otherwise use a linear or a logarithmically-spaced energy binning
else {
// Get task parameters
double emin = (*this)["emin"].real();
double emax = (*this)["emax"].real();
int enumbins = (*this)["enumbins"].integer();
// Initialise log mode for ebinning
bool log = true;
// Check if algorithm is linear
if (ebinalg == "LIN") {
log = false;
}
// Setup energies
energies = GEnergies(enumbins, GEnergy(emin, "TeV"),
GEnergy(emax, "TeV"), log);
} // endelse: ebinalg was not "FILE"
// Return energies
return energies;
}