GModelSpatialDiffuseCube.cpp

Go to the documentation of this file.

/***************************************************************************
 *       GModelSpatialDiffuseCube.cpp - Spatial map cube model class       *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2010-2021 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 <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GModelSpatialDiffuseCube.cpp
 * @brief Spatial map cube model class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#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);
#if defined(G_LEGACY_XML_FORMAT)
const GModelSpatialDiffuseCube g_spatial_cube_legacy_seed(true, "MapCubeFunction");
const GModelSpatialRegistry    g_spatial_cube_legacy_registry(&g_spatial_cube_legacy_seed);
#endif
 
/* __ 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&)"
#define G_READ_FITS                  "GModelSpatialDiffuseCube::read(GFits&)"
 
/* __ 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 Model type constructor
 *
 * @param[in] dummy Dummy flag.
 * @param[in] type Model type.
 *
 * Constructs empty map cube model by specifying a model @p type.
 ***************************************************************************/
GModelSpatialDiffuseCube::GModelSpatialDiffuseCube(const bool&        dummy,
                                                   const std::string& type) :
                          GModelSpatialDiffuse()
{
    // Initialise members
    init_members();
 
    // Set model type
    m_type = type;
 
    // 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.
 * @param[in] gradients Compute gradients?
 * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$)
 *
 * Computes the spatial diffuse model as function of photon parameters.
 *
 * If the @p gradients flag is true the method will also evaluate the partial
 * derivatives of the model.
 ***************************************************************************/
double GModelSpatialDiffuseCube::eval(const GPhoton& photon,
                                      const bool&    gradients) 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();
 
    // Optionally compute partial derivatives
    if (gradients) {
 
        // Compute partial derivatives of the parameter value. In case that
        // the value is negative set the gradient to zero.
        double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0;
        if (value < 0.0) {
            g_value = 0.0;
        }
 
        // Set gradient
        m_value.factor_gradient(g_value);
 
    } // endif: computed partial derivatives
 
    // Make sure that value is not negative
    if (value < 0.0) {
        value = 0.0;
    }
 
    // 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.
 * @exception GException::invalid_return_value
 *            Simulation cone not defined, does not overlap with map cube or
 *            map cube is empty for the specified energy.
 *
 * 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. If no sky direction could be
 * determined, the method throws an GException::invalid_return_value
 * exception.
 *
 * @todo Make sure that spatial model value is taken into account.
 ***************************************************************************/
GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy& energy,
                                     const GTime&   time,
                                     GRan&          ran) const
{
    // Fetch cube
    fetch_cube();
 
    // Determine number of skymap pixels
    int npix = cube().npix();
 
    // 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 happen if the simulation cone has not been defined or if there is
    // no overlap with the sky map or if the sky map is empty
    if (max <= 0.0) {
        std::string msg;
        if (m_mc_max.size() == 0) {
            msg = "Simulation cone has not been defined. Please specify a "
                  "valid simulation cone before calling the method.";
        }
        else {
            msg = "The map cube is empty at "+energy.print()+" within the "
                  "simulation cone. Please specify a valid map cube or "
                  "restrict the simulated energies to values where the map "
                  "cube is non-zero within the simulation cone.";
        }
        throw GException::invalid_return_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_cone.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.
 ***************************************************************************/
bool GModelSpatialDiffuseCube::contains(const GSkyDir& dir,
                                        const double&  margin) const
{
    // Make sure that cube in online
    fetch_cube();
 
    // Initialise containment flag
    bool contains = false;
 
    // Continue only if radius is positive
    if (m_region.radius() > 0.0) {
 
        // Compute distance to centre
        double distance = m_region.centre().dist_deg(dir);
 
        // If distance is smaller than radius plus margin we consider
        // the position to be contained within the bounding circle
        if (distance < m_region.radius() + margin) {
            contains = true;
        }
 
    }
 
    // Return containment
    return (contains);
}
 
 
/***********************************************************************//**
 * @brief Read model from XML element
 *
 * @param[in] xml XML element.
 *
 * @exception GException::invalid_value
 *            Model parameters not found in XML element.
 *
 * Read the map cube information from an XML element. The XML element should
 * have the format
 *
 *     <spatialModel type="DiffuseMapCube" file="test_file.fits">
 *       <parameter name="Normalization" ../>
 *     </spatialModel>
 ***************************************************************************/
void GModelSpatialDiffuseCube::read(const GXmlElement& xml)
{
    // Verify number of model parameters
    gammalib::xml_check_parnum(G_READ, xml, 1);
 
    // If "Normalization" parameter exists then read parameter from this
    // XML element
    if (gammalib::xml_has_par(xml, "Normalization")) {
        const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Normalization");
        m_value.read(*par);
    }
 
    // ... otherwise try reading parameter from "Value" parameter
    #if defined(G_LEGACY_XML_FORMAT)
    else if (gammalib::xml_has_par(xml, "Value")) {
        const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Value");
        m_value.read(*par);
    }
    #endif
 
    // ... otherwise throw an exception
    else {
        #if defined(G_LEGACY_XML_FORMAT)
        std::string msg = "Diffuse map cube model requires either "
                          "\"Normalization\" or \"Value\" parameter.";
        #else
        std::string msg = "Diffuse map cube model requires \"Normalization\" "
                          "parameter.";
        #endif
        throw GException::invalid_value(G_READ, msg);
    }
 
    // Save filename
    m_filename = gammalib::xml_file_expand(xml, xml.attribute("file"));
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Write model into XML element
 *
 * @param[in] xml XML element.
 *
 * Write the map cube information into an XML element. The XML element will
 * have either the format
 *
 *     <spatialModel type="DiffuseMapCube" file="test_file.fits">
 *       <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
 *     </spatialModel>
 *
 * or alternatively
 *
 *     <spatialModel type="DiffuseMapCube" file="test_file.fits">
 *       <parameter name="Value" scale="1" value="1" min="0.1" max="10" free="0"/>
 *     </spatialModel>
 *
 * The latter format is the default for newly written XML elements. 
 ***************************************************************************/
void GModelSpatialDiffuseCube::write(GXmlElement& xml) const
{
    // Verify model type
    gammalib::xml_check_type(G_WRITE, xml, type());
 
    // Get or create parameter
    GXmlElement* value = gammalib::xml_need_par(G_WRITE, xml, m_value.name());
 
    // Write parameter
    m_value.write(*value);
 
    // Set filename
    xml.attribute("file", gammalib::xml_file_reduce(xml, 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] cone Monte Carlo simulation cone.
 *
 * Sets the simulation cone 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::mc_cone(const GSkyRegionCircle& cone) const
{
    // Continue only if the simulation cone has changed
    if (cone != m_mc_cone) {
 
        // Save simulation cone definition
        m_mc_cone = cone;
 
        // Pre-compute 1 - cosine of radius
        m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_cone.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  = cube().npix();
        int nmaps = cube().nmaps();
 
        // Continue only if there are pixels and maps
        if (npix > 0 && nmaps > 0) {
 
            // Initial vectors for total flux and maximum intensity
            std::vector<double> total_flux(nmaps, 0.0);
            std::vector<double> max_intensity(nmaps, 0.0);
 
            // Loop over all map pixels
            for (int k = 0; k < npix; ++k) {
 
                // Continue only if pixel is within MC cone
                double distance = m_mc_cone.centre().dist_deg(m_cube.pix2dir(k));
                if (distance <= m_mc_cone.radius()) {
 
                    // Loop over all maps
                    for (int i = 0; i < nmaps; ++i) {
 
                        // Update total flux
                        double flux = m_cube.flux(k,i);
                        if (flux > 0.0) {
                            total_flux[i] += flux;
                        }
 
                        // Update maximum intensity
                        double value = m_cube(k,i);
                        if (value > max_intensity[i]) {
                            max_intensity[i] = value;
                        }
 
                    } // endfor: looped over all maps
 
                } // endif: pixel within MC cone
 
            } // endfor: looped over all pixels
 
            // Reserve space in cache
            m_mc_max.reserve(nmaps);
            m_mc_spectrum.reserve(nmaps);
 
            // Set maximum intensity and spectrum vectors
            for (int i = 0; i < nmaps; ++i) {
 
                // Store maximum map intensity
                m_mc_max.push_back(max_intensity[i]);
 
                // 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[i] > 0.0) {
                        m_mc_spectrum.append(energy, total_flux[i]);
                    }
 
                }
 
            } // endfor: looped over all maps
 
            // Log maximum intensity and total flux for debugging
            #if defined(G_DEBUG_MC_CACHE)
            std::cout << "GModelSpatialDiffuseCube::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 Read diffuse cube from FITS file
 *
 * @param[in] fits FITS file.
 *
 * @exception GException::invalid_value
 *            ENERGIES extension incompatible with map cube.
 *
 * Reads the diffuse cube sky map from the first extension in the @p fits
 * file and the energy extension from the "ENERGIES" extension.
 ***************************************************************************/
void GModelSpatialDiffuseCube::read(const GFits& fits)
{
    // Initialise skymap
    m_cube.clear();
    m_logE.clear();
 
    // Read cube from the first valid WCS image
    for (int extno = 0; extno < fits.size(); ++extno) {
        const GFitsHDU& hdu = *fits.at(extno);
        if ((hdu.exttype() == GFitsHDU::HT_IMAGE) &&
            (hdu.has_card("NAXIS") && (hdu.integer("NAXIS") >= 2))) {
            m_cube.read(hdu);
            break;
        }
    }
 
    // Read energies
    GEnergies energies;
    energies.read(*(fits.table(gammalib::extname_energies)));
 
    // 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. The \"ENERGIES\" extension table shall "
                          "provide one enegy value for each map in the cube.";
        throw GException::invalid_value(G_READ_FITS, 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 Write diffuse cube into FITS file
 *
 * @param[in] fits FITS file.
 *
 * Writes the diffuse cube sky map and the energy intervals into @p fits
 * file.
 ***************************************************************************/
void GModelSpatialDiffuseCube::write(GFits& fits) const
{
    // Make sure that cube is fetched
    fetch_cube();
 
    // Write cube
    m_cube.write(fits);
 
    // 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(fits);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Returns diffuse cube flux integrated in sky region
 *
 * @param[in] region Sky region.
 * @param[in] srcEng Energy.
 * @param[in] srcTime Time.
 * @return Flux (adimensional or ph/cm2/s).
 *
 * Returns diffuse cube flux within a sky region. The flux is computed by
 * log-log interpolation of the flux of the sky maps that encompass the
 * specified @p srcEng energy. Specifically
 *
 * \f[
 *    F = {\tt value} \times \left \{
 *    \begin{array}{l l}
 *      0 & \mbox{if $F_l=0$ and $F_r=0$} \\
 *      F_l & \mbox{if $F_l\ne0$ and $F_r=0$} \\
 *      F_r & \mbox{if $F_l=0$ and $F_r\ne0$} \\
 *      \exp \left( w_l \log F_l + w_r \log F_r \right) & \mbox{else}
 *    \end{array}
 *    \right .
 * \f]
 *
 * where
 * - \f$F\f$ is the diffuse cube flux
 * - \f${\tt value}\f$ is the cube normalisation that is returned by the
 *   value() method,
 * - \f$w_l\f$ and \f$w_r\f$ are the weighting factors from the logarithmic
 *   interpolation in @p srcEng, with \f$w_l+w_r=1\f$, and
 * - \f$F_l\f$ and \f$F_r\f$ are the fluxes in the sky maps that encompass
 *   the @p srcEng energy.
 *
 * This method fetches the diffuse cube if it was not yet fetched before.
 ***************************************************************************/
double GModelSpatialDiffuseCube::flux(const GSkyRegion& region,
                                      const GEnergy&    srcEng,
                                      const GTime&      srcTime) const
{
    // Fetch cube
    fetch_cube();
 
    // Initialise flux value
    double flux = 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(srcEng.log10MeV());
 
        // Compute map cube flux for the left and right map
        double left_flux  = m_cube.flux(region, m_logE.inx_left());
        double right_flux = m_cube.flux(region, m_logE.inx_right());
 
        // Perform log-log interpolation
        if (left_flux > 0.0 && right_flux > 0.0) {
            double log_flux = m_logE.wgt_left()  * std::log(left_flux) +
                              m_logE.wgt_right() * std::log(right_flux);
            flux = std::exp(log_flux);
        }
        else if (left_flux > 0.0) {
            flux = left_flux;
        }
        else if (right_flux > 0.0) {
            flux = right_flux;
        }
 
        // Multiply by cube normalisation
        flux *= m_value.value();
 
    } // endif: energy information was available
 
    // Return flux
    return flux;
}
 
 
/***********************************************************************//**
 * @brief Print map cube information
 *
 * @param[in] chatter Chattiness.
 * @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 model type
    m_type = "DiffuseMapCube";
 
    // 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_cone.clear();
    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_type     = model.m_type;   // Needed to conserve model type
    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_cone             = model.m_mc_cone;
    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<GModelSpatialDiffuseCube*>(this)->load_cube(m_filename);
        } // end of pragma
 
    } // endif: file not
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Load map cube
 *
 * @param[in] filename Map cube file.
 *
 * @exception GException::invalid_value
 *            Number of maps in cube mismatches number of energy bins.
 *
 * Load diffuse map cube.
 ***************************************************************************/
void GModelSpatialDiffuseCube::load_cube(const GFilename& filename)
{
    // Initialise sky map cube, region and log(energy) nodes
    m_cube.clear();
    m_region.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. 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();
 
    // Get region circle
    m_region = m_cube.region_circle();
 
    // Set simulation cone
    mc_cone(m_region);
 
    // 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. See the
 * mc_cone() for more information.
 ***************************************************************************/
void GModelSpatialDiffuseCube::update_mc_cache(void)
{
    // Initialise sky region to all sky (dummy)
    GSkyRegionCircle cone(0.0, 0.0, 360.0);
 
    // Set cone
    mc_cone(cone);
 
    // 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;
}
 
 
/***********************************************************************//**
 * @brief Set boundary sky region
 ***************************************************************************/
void GModelSpatialDiffuseCube::set_region(void) const
{
    // Make sure that cube in online
    fetch_cube();
 
    // Return
    return;
}