/***************************************************************************
* GModelSpatialDiffuseCube.cpp - Spatial map cube model class *
* ----------------------------------------------------------------------- *
* copyright (C) 2010-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 GModelSpatialDiffuseCube.cpp
* @brief Spatial map cube model class implementation
* @author Juergen Knoedlseder
*/
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include
#endif
#include "GException.hpp"
#include "GTools.hpp"
#include "GMath.hpp"
#include "GHealpix.hpp"
#include "GEnergies.hpp"
#include "GFits.hpp"
#include "GFitsTable.hpp"
#include "GFitsTableCol.hpp"
#include "GXmlElement.hpp"
#include "GModelSpatialRegistry.hpp"
#include "GModelSpatialDiffuseCube.hpp"
/* __ Constants __________________________________________________________ */
/* __ Globals ____________________________________________________________ */
const GModelSpatialDiffuseCube g_spatial_cube_seed;
const GModelSpatialRegistry g_spatial_cube_registry(&g_spatial_cube_seed);
/* __ Method name definitions ____________________________________________ */
#define G_MC "GModelSpatialDiffuseCube::mc(GEnergy&, GTime&, GRan&)"
#define G_READ "GModelSpatialDiffuseCube::read(GXmlElement&)"
#define G_WRITE "GModelSpatialDiffuseCube::write(GXmlElement&)"
#define G_ENERGIES "GModelSpatialDiffuseCube::energies(GEnergies&)"
#define G_LOAD_CUBE "GModelSpatialDiffuseCube::load_cube(GFilename&)"
/* __ Macros _____________________________________________________________ */
/* __ Coding definitions _________________________________________________ */
/* __ Debug definitions __________________________________________________ */
//#define G_DEBUG_MC //!< Debug MC method
//#define G_DEBUG_MC_CACHE //!< Debug MC cache
/*==========================================================================
= =
= Constructors/destructors =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Void constructor
*
* Constructs empty map cube model.
***************************************************************************/
GModelSpatialDiffuseCube::GModelSpatialDiffuseCube(void) :
GModelSpatialDiffuse()
{
// Initialise members
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief XML constructor
*
* @param[in] xml XML element.
*
* Constructs map cube model by extracting information from an XML element.
* See the read() method for more information about the expected structure
* of the XML element.
***************************************************************************/
GModelSpatialDiffuseCube::GModelSpatialDiffuseCube(const GXmlElement& xml) :
GModelSpatialDiffuse()
{
// Initialise members
init_members();
// Read information from XML element
read(xml);
// Return
return;
}
/***********************************************************************//**
* @brief Filename constructor
*
* @param[in] filename File name.
* @param[in] value Normalization factor (defaults to 1).
*
* Constructs map cube model by loading a map cube from @p filename and by
* assigning the normalization @p value.
***************************************************************************/
GModelSpatialDiffuseCube::GModelSpatialDiffuseCube(const GFilename& filename,
const double& value) :
GModelSpatialDiffuse()
{
// Initialise members
init_members();
// Set parameter
m_value.value(value);
// Perform autoscaling of parameter
autoscale();
// Store filename
m_filename = filename;
// Return
return;
}
/***********************************************************************//**
* @brief Sky map constructor
*
* @param[in] cube Sky map cube.
* @param[in] energies Sky map energies.
* @param[in] value Normalization factor (defaults to 1).
*
* Constructs map cube model by extracting a @p cube from a sky map. The
* constructor also assigns the energy values for all maps and sets the
* scaling @p value. The filename will remain blank.
***************************************************************************/
GModelSpatialDiffuseCube::GModelSpatialDiffuseCube(const GSkyMap& cube,
const GEnergies& energies,
const double& value) :
GModelSpatialDiffuse()
{
// Initialise members
init_members();
// Set parameter
m_value.value(value);
// Perform autoscaling of parameter
autoscale();
// Set map cube
this->cube(cube);
// Set energies
this->energies(energies);
// Return
return;
}
/***********************************************************************//**
* @brief Copy constructor
*
* @param[in] model Map cube model.
***************************************************************************/
GModelSpatialDiffuseCube::GModelSpatialDiffuseCube(const GModelSpatialDiffuseCube& model) :
GModelSpatialDiffuse(model)
{
// Initialise members
init_members();
// Copy members
copy_members(model);
// Return
return;
}
/***********************************************************************//**
* @brief Destructor
***************************************************************************/
GModelSpatialDiffuseCube::~GModelSpatialDiffuseCube(void)
{
// Free members
free_members();
// Return
return;
}
/*==========================================================================
= =
= Operators =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Assignment operator
*
* @param[in] model Map cube model.
* @return Map cube model.
***************************************************************************/
GModelSpatialDiffuseCube& GModelSpatialDiffuseCube::operator=(const GModelSpatialDiffuseCube& model)
{
// Execute only if object is not identical
if (this != &model) {
// Copy base class members
this->GModelSpatialDiffuse::operator=(model);
// Free members
free_members();
// Initialise members
init_members();
// Copy members
copy_members(model);
} // endif: object was not identical
// Return
return *this;
}
/*==========================================================================
= =
= Public methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Clear map cube model
***************************************************************************/
void GModelSpatialDiffuseCube::clear(void)
{
// Free class members (base and derived classes, derived class first)
free_members();
this->GModelSpatialDiffuse::free_members();
this->GModelSpatial::free_members();
// Initialise members
this->GModelSpatial::init_members();
this->GModelSpatialDiffuse::init_members();
init_members();
// Return
return;
}
/***********************************************************************//**
* @brief Clone map cube model
*
* @return Pointer to deep copy of map cube model.
***************************************************************************/
GModelSpatialDiffuseCube* GModelSpatialDiffuseCube::clone(void) const
{
// Clone map cube model
return new GModelSpatialDiffuseCube(*this);
}
/***********************************************************************//**
* @brief Evaluate function
*
* @param[in] photon Incident photon.
* @return Model value.
*
* Computes the spatial diffuse model as function of photon parameters.
***************************************************************************/
double GModelSpatialDiffuseCube::eval(const GPhoton& photon) const
{
// Get log-log interpolated cube intensity
double intensity = cube_intensity(photon);
// Set the intensity times the scaling factor as model value
double value = intensity * m_value.value();
// Make sure that value is not negative
if (value < 0.0) {
value = 0.0;
}
// Return value
return value;
}
/***********************************************************************//**
* @brief Evaluate function and gradients
*
* @param[in] photon Incident photon.
* @return Model value.
*
* Computes the spatial diffuse model as function of photon parameters and
* sets the value gradient.
***************************************************************************/
double GModelSpatialDiffuseCube::eval_gradients(const GPhoton& photon) const
{
// Get log-log interpolated cube intensity
double intensity = cube_intensity(photon);
// Compute the model value
double value = intensity * m_value.value();
// Compute partial derivatives of the parameter value
double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0;
// Make sure that value is not negative
if (value < 0.0) {
value = 0.0;
g_value = 0.0;
}
// Set gradient to 0 (circumvent const correctness)
const_cast(this)->m_value.factor_gradient(g_value);
// Return value
return value;
}
/***********************************************************************//**
* @brief Returns MC sky direction
*
* @param[in] energy Photon energy.
* @param[in] time Photon arrival time.
* @param[in,out] ran Random number generator.
* @return Sky direction.
*
* @exception GException::invalid_value
* No map cube defined.
* No energy boundaries defined.
* Simulation cone not defined or does not overlap with map cube.
*
* Returns a random sky direction according to the intensity distribution of
* the model sky map and the specified energy. The method uses a rejection
* method to determine the sky direction.
***************************************************************************/
GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy& energy,
const GTime& time,
GRan& ran) const
{
// Fetch cube
fetch_cube();
// Determine number of skymap pixels
int npix = pixels();
// Throw an exception if there are no sky map pixels
if (npix <= 0) {
std::string msg = "No map cube defined. Please specify a valid map cube.";
throw GException::invalid_value(G_MC, msg);
}
// Throw an exception if no energy boundaries are defined
if (m_ebounds.size() < 1) {
std::string msg = "The energy boundaries of the maps in the cube have "
"not been defined. Maybe the map cube file is missing "
"the \"ENERGIES\" extension which defines the energy "
"of each map in the cube. Please provide the energy "
"information.";
throw GException::invalid_value(G_MC, msg);
}
// Set energy for interpolation
m_logE.set_value(energy.log10MeV());
// Compute maximum map value within simulation cone by log-log interpolation
// of the maximum values that are stored for each map
double max = 0.0;
double max_left = m_mc_max[m_logE.inx_left()];
double max_right = m_mc_max[m_logE.inx_right()];
if (max_left > 0.0 && max_right > 0.0) {
double max_log = m_logE.wgt_left() * std::log(max_left) +
m_logE.wgt_right() * std::log(max_right);
max = std::exp(max_log);
}
else if (max_left > 0.0) {
max = max_left;
}
else if (max_right > 0.0) {
max = max_right;
}
// Throw an exception if the maximum MC intensity is not positive. This
// can be the case because the simulation cone has not been defined or
// because it does not overlap with the sky map
if (max <= 0.0) {
std::string msg = "Simulation cone has not been defined or does not "
"overlap with the sky map. Please specify a valid "
"simulation cone.";
throw GException::invalid_value(G_MC, msg);
}
// Allocate sky direction
GSkyDir dir;
// Debug option: initialise counter
#if defined(G_DEBUG_MC)
int num_iterations = 0;
#endif
// Get sky direction
while (true) {
// Debug option: increment counter
#if defined(G_DEBUG_MC)
num_iterations++;
#endif
// Simulate random sky direction within Monte Carlo simulation cone
double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
gammalib::rad2deg;
double phi = 360.0 * ran.uniform();
dir = m_mc_centre;
dir.rotate_deg(phi, theta);
// Get map value at simulated sky direction
double value = 0.0;
double value_left = m_cube(dir, m_logE.inx_left());
double value_right = m_cube(dir, m_logE.inx_right());
if (value_left > 0.0 && value_right > 0.0) {
double value_log = m_logE.wgt_left() * std::log(value_left) +
m_logE.wgt_right() * std::log(value_right);
value = std::exp(value_log);
}
else if (value_left > 0.0) {
value = value_left;
}
else if (value_right > 0.0) {
value = value_right;
}
// If map value is non-positive then simulate a new sky direction
if (value <= 0.0) {
continue;
}
// Get uniform random number
double uniform = ran.uniform() * max;
// Exit loop if the random number is not larger than the map cube value
if (uniform <= value) {
break;
}
} // endwhile: loop until sky direction was accepted
// Debug option: log counter
#if defined(G_DEBUG_MC)
std::cout << num_iterations << " ";
#endif
// Return sky direction
return dir;
}
/***********************************************************************//**
* @brief Signals whether model contains sky direction
*
* @param[in] dir Sky direction.
* @param[in] margin Margin to be added to sky direction (degrees)
* @return True if the model contains the sky direction.
*
* Signals whether a sky direction falls within the bounding circle of
* the diffuse cube.
*
* @todo To be implemented.
***************************************************************************/
bool GModelSpatialDiffuseCube::contains(const GSkyDir& dir,
const double& margin) const
{
return (true);
}
/***********************************************************************//**
* @brief Read model from XML element
*
* @param[in] xml XML element.
*
* @exception GException::model_invalid_parnum
* Invalid number of model parameters found in XML element.
* @exception GException::model_invalid_parnames
* Invalid model parameter names found in XML element.
*
* Read the map cube information from an XML element. The XML element should
* have either the format
*
*
*
*
*
* or alternatively
*
*
*
*
*
***************************************************************************/
void GModelSpatialDiffuseCube::read(const GXmlElement& xml)
{
// Clear model
clear();
// Verify that XML element has exactly 1 parameters
if (xml.elements() != 1 || xml.elements("parameter") != 1) {
throw GException::model_invalid_parnum(G_READ, xml,
"Map cube spatial model requires exactly 1 parameter.");
}
// Get pointer on model parameter
const GXmlElement* par = xml.element("parameter", 0);
// Get value
if (par->attribute("name") == "Normalization" ||
par->attribute("name") == "Value") {
m_value.read(*par);
}
else {
throw GException::model_invalid_parnames(G_READ, xml,
"Map cube spatial model requires either \"Value\" or"
" \"Normalization\" parameter.");
}
// Save filename
m_filename = xml.attribute("file");
// Return
return;
}
/***********************************************************************//**
* @brief Write model into XML element
*
* @param[in] xml XML element.
*
* @exception GException::model_invalid_spatial
* Existing XML element is not of type "MapCubeFunction"
* @exception GException::model_invalid_parnum
* Invalid number of model parameters found in XML element.
* @exception GException::model_invalid_parnames
* Invalid model parameter names found in XML element.
*
* Write the map cube information into an XML element. The XML element will
* have either the format
*
*
*
*
*
* or alternatively
*
*
*
*
*
* The latter format is the default for newly written XML elements.
***************************************************************************/
void GModelSpatialDiffuseCube::write(GXmlElement& xml) const
{
// Set model type
if (xml.attribute("type") == "") {
xml.attribute("type", "MapCubeFunction");
}
// Verify model type
if (xml.attribute("type") != "MapCubeFunction") {
throw GException::model_invalid_spatial(G_WRITE, xml.attribute("type"),
"Spatial model is not of type \"MapCubeFunction\".");
}
// If XML element has 0 nodes then append parameter node. The name
// of the node is "Normalization" as this is the Fermi-LAT standard.
// We thus assure that the XML files will be compatible with
// Fermi-LAT.
if (xml.elements() == 0) {
xml.append(GXmlElement("parameter name=\"Normalization\""));
}
// Verify that XML element has exactly 1 parameter
if (xml.elements() != 1 || xml.elements("parameter") != 1) {
throw GException::model_invalid_parnum(G_WRITE, xml,
"Map cube spatial model requires exactly 1 parameter.");
}
// Get pointers on model parameter
GXmlElement* par = xml.element("parameter", 0);
// Set or update parameter
if (par->attribute("name") == "Normalization" ||
par->attribute("name") == "Value") {
m_value.write(*par);
}
else {
throw GException::model_invalid_parnames(G_WRITE, xml,
"Map cube spatial model requires either \"Value\" or"
" \"Normalization\" parameter.");
}
// Set filename
xml.attribute("file", m_filename);
// Return
return;
}
/***********************************************************************//**
* @brief Set map cube
*
* @param[in] cube Sky map.
*
* Set the map cube of the spatial map cube model.
***************************************************************************/
void GModelSpatialDiffuseCube::cube(const GSkyMap& cube)
{
// Clear filename and signal map as unloaded
m_filename.clear();
m_loaded = false;
// Assign map
m_cube = cube;
// Update Monte Carlo cache
update_mc_cache();
// Return
return;
}
/***********************************************************************//**
* @brief Set energies for map cube
*
* @param[in] energies Sky map energies.
*
* @exception GException::invalid_argument
* Specified sky map energies incompatible with map cube.
*
* Sets the energies for the map cube.
***************************************************************************/
void GModelSpatialDiffuseCube::energies(const GEnergies& energies)
{
// Initialise energies
m_logE.clear();
// Fetch cube
fetch_cube();
// Extract number of energies in vector
int num = energies.size();
// Check if energy binning is consistent with number of maps in the cube
if (num != m_cube.nmaps() ) {
std::string msg = "Number of specified energies ("+gammalib::str(num)+")"
" does not match the number of maps ("
""+gammalib::str(m_cube.nmaps())+" in the map cube.\n"
"The energies argument shall provide a vector of length"
" "+gammalib::str(m_cube.nmaps())+".";
throw GException::invalid_argument(G_ENERGIES, msg);
}
// Set log10(energy) nodes, where energy is in units of MeV
for (int i = 0; i < num; ++i) {
m_logE.append(energies[i].log10MeV());
}
// Set energy boundaries
set_energy_boundaries();
// Update MC cache
update_mc_cache();
// Return
return;
}
/***********************************************************************//**
* @brief Returns map cube energies
*
* @return Map cube energies.
*
* Returns the energies for the map cube in a vector.
***************************************************************************/
GEnergies GModelSpatialDiffuseCube::energies(void)
{
// Initialise energies container
GEnergies energies;
// Fetch cube
fetch_cube();
// Get number of map energies
int num = m_logE.size();
// Continue only if there are maps in the cube
if (num > 0) {
// Reserve space for all energies
energies.reserve(num);
// Set log10(energy) nodes, where energy is in units of MeV
for (int i = 0; i < num; ++i) {
GEnergy energy;
energy.log10MeV(m_logE[i]);
energies.append(energy);
}
} // endif: there were maps in the cube
// Return energies
return energies;
}
/***********************************************************************//**
* @brief Set Monte Carlo simulation cone
*
* @param[in] centre Simulation cone centre.
* @param[in] radius Simulation cone radius (degrees).
*
* Sets the simulation cone centre and radius that defines the directions
* that will be simulated using the mc() method and pre-computes the maximum
* intensity and the spatially integrated flux of each map within the
* simulation cone region.
***************************************************************************/
void GModelSpatialDiffuseCube::set_mc_cone(const GSkyDir& centre,
const double& radius) const
{
// Continue only if the simulation cone has changed
if ((centre != m_mc_centre) || (radius != m_mc_radius)) {
// Save simulation cone definition
m_mc_centre = centre;
m_mc_radius = radius;
// Pre-compute 1 - cosine of radius
m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_radius*gammalib::deg2rad);
// Initialise Monte Carlo cache
m_mc_max.clear();
m_mc_spectrum.clear();
// Fetch cube
fetch_cube();
// Determine number of cube pixels and maps
int npix = pixels();
int nmaps = maps();
// Continue only if there are pixels and maps
if (npix > 0 && nmaps > 0) {
// Reserve space in cache
m_mc_max.reserve(nmaps);
m_mc_spectrum.reserve(nmaps);
// Loop over all maps
for (int i = 0; i < nmaps; ++i) {
// Compute flux and maximum map intensity within the
// simulation cone
double total_flux = 0.0;
double max_intensity = 0.0;
for (int k = 0; k < npix; ++k) {
double distance = centre.dist_deg(m_cube.pix2dir(k));
if (distance <= radius) {
double flux = m_cube.flux(k,i);
if (flux > 0.0) {
total_flux += flux;
}
double value = m_cube(k,i);
if (value > max_intensity) {
max_intensity = value;
}
}
}
// Store maximum map intensity
m_mc_max.push_back(max_intensity);
// Store flux as spectral node
if (m_logE.size() == nmaps) {
// Set map energy
GEnergy energy;
energy.log10MeV(m_logE[i]);
// Only append node if the integrated flux is positive
// (as we do a log-log interpolation)
if (total_flux > 0.0) {
m_mc_spectrum.append(energy, total_flux);
}
}
} // endfor: looped over all maps
// Log maximum intensity and total flux for debugging
#if defined(G_DEBUG_MC_CACHE)
std::cout << "GModelSpatialDiffuseCube::set_mc_cone:" << std::endl;
std::cout << " Maximum map intensity:" << std::endl;
for (int i = 0; i < m_mc_max.size(); ++i) {
GEnergy energy;
energy.log10MeV(m_logE[i]);
std::cout << " " << i;
std::cout << " " << energy;
std::cout << " " << m_mc_max[i] << " ph/cm2/s/sr/MeV";
std::cout << std::endl;
}
std::cout << " Spatially integrated flux:" << std::endl;
for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
std::cout << " " << i;
std::cout << " " << m_mc_spectrum.energy(i);
std::cout << " " << m_mc_spectrum.intensity(i);
std::cout << " ph/cm2/s/MeV" << std::endl;
}
#endif
} // endif: there were cube pixels and maps
} // endif: simulation cone has changed
// Return
return;
}
/***********************************************************************//**
* @brief Load cube into the model class
*
* @param[in] filename cube file.
*
* Loads cube into the model class.
***************************************************************************/
void GModelSpatialDiffuseCube::load(const GFilename& filename)
{
// Load cube
load_cube(filename);
// Update Monte Carlo cache
update_mc_cache();
// Return
return;
}
/***********************************************************************//**
* @brief Save cube into FITS file
*
* @param[in] filename Cube FITS file name.
* @param[in] clobber Overwrite existing FITS file (default: false).
*
* Saves spatial cube model into a FITS file.
***************************************************************************/
void GModelSpatialDiffuseCube::save(const GFilename& filename,
const bool& clobber) const
{
// Create empty FITS file
GFits fits;
// Write event cube
write(fits);
// Save FITS file
fits.saveto(filename, clobber);
// Return
return;
}
/***********************************************************************//**
* @brief Write skymap into FITS file
*
* @param[in] file FITS file pointer.
***************************************************************************/
void GModelSpatialDiffuseCube::write(GFits& file) const
{
// Write cube
m_cube.write(file);
// Create energy intervals
GEnergies energies;
for (int i = 0; i < m_logE.size(); ++i) {
double energy = std::pow(10.0, m_logE[i]);
energies.append(GEnergy(energy, "MeV"));
}
// Write energy intervals
energies.write(file);
// Return
return;
}
/***********************************************************************//**
* @brief Print map cube information
*
* @param[in] chatter Chattiness (defaults to NORMAL).
* @return String containing model information.
***************************************************************************/
std::string GModelSpatialDiffuseCube::print(const GChatter& chatter) const
{
// Initialise result string
std::string result;
// Continue only if chatter is not silent
if (chatter != SILENT) {
// Append header
result.append("=== GModelSpatialDiffuseCube ===");
// Append parameters
result.append("\n"+gammalib::parformat("Map cube file"));
result.append(m_filename);
if (m_loaded) {
result.append(" [loaded]");
}
result.append("\n"+gammalib::parformat("Number of parameters"));
result.append(gammalib::str(size()));
for (int i = 0; i < size(); ++i) {
result.append("\n"+m_pars[i]->print(chatter));
}
// Append detailed information only if a map cube exists
if (m_cube.npix() > 0) {
// NORMAL: Append sky map
if (chatter >= NORMAL) {
result.append("\n"+m_cube.print(chatter));
}
// EXPLICIT: Append energy nodes
if (chatter >= EXPLICIT && m_logE.size() > 0) {
result.append("\n"+gammalib::parformat("Map energy values"));
if (m_logE.size() > 0) {
for (int i = 0; i < m_logE.size(); ++i) {
result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
result.append(gammalib::str(std::pow(10.0, m_logE[i])));
result.append(" MeV (log10E=");
result.append(gammalib::str(m_logE[i]));
result.append(")");
if (m_ebounds.size() == m_logE.size()) {
result.append(" [");
result.append(m_ebounds.emin(i).print());
result.append(", ");
result.append(m_ebounds.emax(i).print());
result.append("]");
}
}
}
else {
result.append("not specified");
}
}
// VERBOSE: Append MC cache
if (chatter >= VERBOSE) {
result.append("\n"+gammalib::parformat("Map flux"));
if (m_mc_spectrum.nodes() > 0) {
for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
result.append(gammalib::str(m_mc_spectrum.intensity(i)));
}
}
else {
result.append("not specified");
}
}
} // endif: map cube exists
} // endif: chatter was not silent
// Return result
return result;
}
/*==========================================================================
= =
= Private methods =
= =
==========================================================================*/
/***********************************************************************//**
* @brief Initialise class members
***************************************************************************/
void GModelSpatialDiffuseCube::init_members(void)
{
// Initialise Value
m_value.clear();
m_value.name("Normalization");
m_value.value(1.0);
m_value.scale(1.0);
m_value.range(0.001, 1000.0);
m_value.gradient(0.0);
m_value.fix();
m_value.has_grad(true);
// Set parameter pointer(s)
m_pars.clear();
m_pars.push_back(&m_value);
// Initialise other members
m_filename.clear();
m_cube.clear();
m_logE.clear();
m_ebounds.clear();
m_loaded = false;
// Initialise MC cache
m_mc_centre.clear();
m_mc_radius = -1.0; //!< Signal that initialisation is needed
m_mc_one_minus_cosrad = 1.0;
m_mc_max.clear();
m_mc_spectrum.clear();
// Return
return;
}
/***********************************************************************//**
* @brief Copy class members
*
* @param[in] model Spatial map cube model.
***************************************************************************/
void GModelSpatialDiffuseCube::copy_members(const GModelSpatialDiffuseCube& model)
{
// Copy members
m_value = model.m_value;
m_filename = model.m_filename;
m_cube = model.m_cube;
m_logE = model.m_logE;
m_ebounds = model.m_ebounds;
m_loaded = model.m_loaded;
// Copy MC cache
m_mc_centre = model.m_mc_centre;
m_mc_radius = model.m_mc_radius;
m_mc_one_minus_cosrad = model.m_mc_one_minus_cosrad;
m_mc_max = model.m_mc_max;
m_mc_spectrum = model.m_mc_spectrum;
// Set parameter pointer(s)
m_pars.clear();
m_pars.push_back(&m_value);
// Return
return;
}
/***********************************************************************//**
* @brief Delete class members
***************************************************************************/
void GModelSpatialDiffuseCube::free_members(void)
{
// Return
return;
}
/***********************************************************************//**
* @brief Fetch cube
*
* Load diffuse cube if it is not yet loaded. The loading is thread save.
***************************************************************************/
void GModelSpatialDiffuseCube::fetch_cube(void) const
{
// Load cube if it is not yet loaded
if (!m_loaded && !m_filename.is_empty()) {
// Put in a OMP critical zone
#pragma omp critical(GModelSpatialDiffuseCube_fetch_cube)
{
const_cast(this)->load_cube(m_filename);
} // end of pragma
} // endif: file not
// Return
return;
}
/***********************************************************************//**
* @brief Load cube
*
* @param[in] filename Cube file.
*
* @exception GException::invalid_value
* Number of maps in cube mismatches number of energy bins.
*
* Load diffuse cube.
***************************************************************************/
void GModelSpatialDiffuseCube::load_cube(const GFilename& filename)
{
// Initialise skymap
m_cube.clear();
m_logE.clear();
// Store filename of cube (for XML writing). Note that we do
// not expand any environment variable at this level, so that
// if we write back the XML element we write the filepath with
// the environment variables
m_filename = filename;
// Load cube
m_cube.load(filename);
// Load energies
GEnergies energies(filename);
// Extract number of energy bins
int num = energies.size();
// Check if energy binning is consistent with primary image hdu
if (num != m_cube.nmaps() ) {
std::string msg = "Number of energies in \"ENERGIES\""
" extension ("+gammalib::str(num)+")"
" does not match the number of maps ("+
gammalib::str(m_cube.nmaps())+" in the"
" map cube.\n"
"The \"ENERGIES\" extension table shall"
" provide one enegy value for each map"
" in the cube.";
throw GException::invalid_value(G_LOAD_CUBE, msg);
}
// Set log10(energy) nodes, where energy is in units of MeV
for (int i = 0; i < num; ++i) {
m_logE.append(energies[i].log10MeV());
}
// Signal that cube has been loaded
m_loaded = true;
// Set energy boundaries
set_energy_boundaries();
// Return
return;
}
/***********************************************************************//**
* @brief Set energy boundaries
*
* Computes energy boundaries from the energy nodes. The boundaries will be
* computed as the centre in log10(energy) between the nodes. If only a
* single map exists, no boundaries will be computed.
***************************************************************************/
void GModelSpatialDiffuseCube::set_energy_boundaries(void)
{
// Initialise energy boundaries
m_ebounds.clear();
// Determine number of energy bins
int num = m_logE.size();
// Continue only if there are at least two energy nodes
if (num > 1) {
// Loop over all nodes
for (int i = 0; i < num; ++i) {
// Calculate minimum energy
double e_min = (i == 0) ? m_logE[i] - 0.5 * (m_logE[i+1] - m_logE[i])
: m_logE[i] - 0.5 * (m_logE[i] - m_logE[i-1]);
// Calculate maximum energy
double e_max = (i < num-1) ? m_logE[i] + 0.5 * (m_logE[i+1] - m_logE[i])
: m_logE[i] + 0.5 * (m_logE[i] - m_logE[i-1]);
// Set energy boundaries
GEnergy emin;
GEnergy emax;
emin.log10MeV(e_min);
emax.log10MeV(e_max);
// Append energy bin to energy boundary arra
m_ebounds.append(emin, emax);
} // endfor: looped over energy nodes
} // endif: there were at least two energy nodes
// Return
return;
}
/***********************************************************************//**
* @brief Update Monte Carlo cache
*
* Initialise the cache for Monte Carlo sampling of the map cube. The Monte
* Carlo cache consists of a linear array that maps a value between 0 and 1
* into the skymap pixel for all maps in the cube.
***************************************************************************/
void GModelSpatialDiffuseCube::update_mc_cache(void)
{
// Set centre and radius to all sky
GSkyDir centre;
double radius = 360.0;
// Compute cache
set_mc_cone(centre, radius);
// Return
return;
}
/***********************************************************************//**
* @brief Compute cube intensity by log-log interpolation
*
* @param[in] photon Incident photon.
* @return Cube intensity value.
***************************************************************************/
double GModelSpatialDiffuseCube::cube_intensity(const GPhoton& photon) const
{
// Fetch cube
fetch_cube();
// Initialise intensity value
double intensity = 0.0;
// Continue only if there is energy information for the map cube
if (m_logE.size() > 0) {
// Set energy for interpolation in log-energy
m_logE.set_value(photon.energy().log10MeV());
// Compute map cube intensity for the left and right map
double left_intensity = m_cube(photon.dir(), m_logE.inx_left());
double right_intensity = m_cube(photon.dir(), m_logE.inx_right());
// Perform log-log interpolation
if (left_intensity > 0.0 && right_intensity > 0.0) {
double log_intensity = m_logE.wgt_left() * std::log(left_intensity) +
m_logE.wgt_right() * std::log(right_intensity);
intensity = std::exp(log_intensity);
}
else if (left_intensity > 0.0) {
intensity = left_intensity;
}
else if (right_intensity > 0.0) {
intensity = right_intensity;
}
} // endif: energy information was available
// Return intensity
return intensity;
}