From c4ca5cebf186ac6c4eb178cd6c9cb39049157f38 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 18 Feb 2016 15:47:03 +0100 Subject: [PATCH] Restrict Monte Carlo simulations of GModelSpatialDiffuseConst to simulation cone (#1674). --- ChangeLog | 1 + NEWS | 5 +++++ include/GModelSpatialDiffuseConst.hpp | 27 +++++---------------------- src/model/GModelSpatialDiffuseConst.cpp | 216 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------- 4 files changed, 150 insertions(+), 99 deletions(-) diff --git a/ChangeLog b/ChangeLog index 3152661..43eed90 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.1.0 released ========================= + Restrict Monte Carlo simulations in GModelSpatialDiffuseConst to simulation cone (#1674) Use default GTime in GModelSpectral eval() and eval_gradients() methods (#1672) Optimise GModelSpatialDiffuseMap::set_mc_cone() method (#1673) Optimise GModelSpatialDiffuseCube::set_mc_cone() method (#1673) diff --git a/NEWS b/NEWS index eb53d31..98eacba 100644 --- a/NEWS +++ b/NEWS @@ -111,6 +111,11 @@ None 8. Model module --------------- +Sky directions drawn by GModelSpatialDiffuseConst::mc() are now restricted +to the simulation cone that has previously been defined using the mc_norm() +method. This speeds up considerbly the Monte Carlo simulations for an +isotropic diffuse model (#1674). + The GModelSpectral eval() and eval_gradients() methods use an empty GTim() object as default time argument since all spectral models are so far time- independent. This avoids that a user needs to specify a dummy time argument diff --git a/include/GModelSpatialDiffuseConst.hpp b/include/GModelSpatialDiffuseConst.hpp index 55f5722..3351def 100644 --- a/include/GModelSpatialDiffuseConst.hpp +++ b/include/GModelSpatialDiffuseConst.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialDiffuseConst.hpp - Spatial isotropic model class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2015 by Juergen Knoedlseder * + * copyright (C) 2010-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -86,7 +86,9 @@ protected: void free_members(void); // Protected members - GModelPar m_value; //!< Value + GModelPar m_value; //!< Value + mutable GSkyDir m_mc_centre; //!< Simulation cone centre + mutable double m_mc_cos_radius; //!< Cosine of simulation cone radius }; @@ -146,29 +148,10 @@ void GModelSpatialDiffuseConst::value(const double& value) /***********************************************************************//** - * @brief Return normalization of constant diffuse source for Monte Carlo - * simulations - * - * @param[in] dir Centre of simulation cone. - * @param[in] radius Radius of simulation cone (degrees). - * @return Normalization. - * - * Returns the normalization of a constant diffuse source. The normalization - * is given by the model value and does not depend on @p dir and @p radius. - ***************************************************************************/ -inline -double GModelSpatialDiffuseConst::mc_norm(const GSkyDir& dir, - const double& radius) const -{ - return (value()); -} - - -/***********************************************************************//** * @brief Signals whether model contains sky direction * * @param[in] dir Sky direction. - * @param[in] margin Margin to be added to sky direction (degrees) + * @param[in] margin Margin to be added to sky direction (deg, default: 0.0). * @return True. * * Signals whether a sky direction falls within the bounding circle of diff --git a/src/model/GModelSpatialDiffuseConst.cpp b/src/model/GModelSpatialDiffuseConst.cpp index 18b671b..0e1e7d3 100644 --- a/src/model/GModelSpatialDiffuseConst.cpp +++ b/src/model/GModelSpatialDiffuseConst.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialDiffuseConst.cpp - Spatial isotropic model class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2015 by Juergen Knoedlseder * + * copyright (C) 2010-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -41,7 +41,7 @@ const GModelSpatialDiffuseConst g_spatial_const_seed; const GModelSpatialRegistry g_spatial_const_registry(&g_spatial_const_seed); /* __ Method name definitions ____________________________________________ */ -#define G_MC "GModelSpatialDiffuseConst::mc(GRan&)" +#define G_MC_NORM "GModelSpatialDiffuseConst::mc_norm(GSkyDir&, double&)" #define G_READ "GModelSpatialDiffuseConst::read(GXmlElement&)" #define G_WRITE "GModelSpatialDiffuseConst::write(GXmlElement&)" @@ -60,6 +60,8 @@ const GModelSpatialRegistry g_spatial_const_registry(&g_spatial_const_seed); /***********************************************************************//** * @brief Void constructor + * + * Construct empty isotropic spatial model. ***************************************************************************/ GModelSpatialDiffuseConst::GModelSpatialDiffuseConst(void) : GModelSpatialDiffuse() @@ -98,11 +100,11 @@ GModelSpatialDiffuseConst::GModelSpatialDiffuseConst(const GXmlElement& xml) : /***********************************************************************//** * @brief Value constructor * - * @param[in] value Isotropic intensity value (units sr^-1). + * @param[in] value Isotropic intensity value (\f$sr^{-1}\f$). * - * Constructs isotropic spatial model by assigning the intensity @p value of - * the diffuse emission. This constructor explicitly sets the m_value - * parameter of the model. + * Constructs isotropic spatial model by assigning an intensity of @p value + * in units of \f$sr^{-1}\f$ to the diffuse emission. This constructor + * explicitly sets the @a m_value member of the model. ***************************************************************************/ GModelSpatialDiffuseConst::GModelSpatialDiffuseConst(const double& value) : GModelSpatialDiffuse() @@ -125,6 +127,9 @@ GModelSpatialDiffuseConst::GModelSpatialDiffuseConst(const double& value) : * @brief Copy constructor * * @param[in] model Isotropic spatial model. + * + * Constructs isotropic spatial model by copying another isotropic spatial + * model. ***************************************************************************/ GModelSpatialDiffuseConst::GModelSpatialDiffuseConst(const GModelSpatialDiffuseConst& model) : GModelSpatialDiffuse(model) @@ -142,6 +147,8 @@ GModelSpatialDiffuseConst::GModelSpatialDiffuseConst(const GModelSpatialDiffuseC /***********************************************************************//** * @brief Destructor + * + * Destructs isotropic spatial model. ***************************************************************************/ GModelSpatialDiffuseConst::~GModelSpatialDiffuseConst(void) { @@ -164,6 +171,8 @@ GModelSpatialDiffuseConst::~GModelSpatialDiffuseConst(void) * * @param[in] model Isotropic spatial model. * @return Isotropic spatial model. + * + * Assigns an isotropic spatial model. ***************************************************************************/ GModelSpatialDiffuseConst& GModelSpatialDiffuseConst::operator=(const GModelSpatialDiffuseConst& model) { @@ -197,10 +206,13 @@ GModelSpatialDiffuseConst& GModelSpatialDiffuseConst::operator=(const GModelSpat /***********************************************************************//** * @brief Clear isotropic spatial model + * + * Clears the isotropic spatial model. This method is equivalent to creating + * an isotropic spatial model using the void constructor. ***************************************************************************/ void GModelSpatialDiffuseConst::clear(void) { - // Free class members (base and derived classes, derived class first) + // Free class members free_members(); this->GModelSpatialDiffuse::free_members(); this->GModelSpatial::free_members(); @@ -219,6 +231,8 @@ void GModelSpatialDiffuseConst::clear(void) * @brief Clone isotropic spatial model * * @return Pointer to deep copy of isotropic spatial model. + * + * Returns a pointer to a deep copy of the isotropic spatial model. ***************************************************************************/ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const { @@ -228,15 +242,22 @@ GModelSpatialDiffuseConst* GModelSpatialDiffuseConst::clone(void) const /***********************************************************************//** - * @brief Evaluate function (in units of sr^-1) + * @brief Evaluate isotropic spatial model value + * + * @param[in] photon Incident photon. + * @return Model value (\f$sr^{-1}\f$). + * + * Evaluates the isotropic spatial model for a given @p photon, characterised + * by a sky direction, energy and arrival time \f$(\vec{p},E,t)\f$. By + * definition, the isotropic spatial model and gradient are independent of + * sky direction, energy and time. The model value is given by * - * @param[in] photon Incident photon (ignored). - * @return Model value (units of sr^-1). + * \f[ + * M_\mathrm{S}(\vec{p}|E,t) = \frac{v}{4\pi} + * \f] * - * Evaluates the spatial part for an isotropic source model. By definition - * the model value is independent from the sky direction. The model value is - * given by the value parameter divided by the solid angle of the sphere - * \f$4 \pi\f$. + * where \f$v\f$ is the value() parameter, which is divided by the solid + * angle \f$4\pi\f$ of the celestial sphere. ***************************************************************************/ double GModelSpatialDiffuseConst::eval(const GPhoton& photon) const { @@ -252,14 +273,32 @@ double GModelSpatialDiffuseConst::eval(const GPhoton& photon) const /***********************************************************************//** - * @brief Evaluate function and gradients (in units of sr^-1) + * @brief Evaluate isotropic spatial model value and parameter gradient + * + * @param[in] photon Incident photon. + * @return Model value (\f$sr^{-1}\f$). + * + * Evaluates the isotropic spatial model for a given @p photon, characterised + * by a sky direction, energy and arrival time \f$(\vec{p},E,t)\f$, and + * computes the gradient of the value() parameter. By definition, the + * isotropic spatial model and gradient are independent of sky direction, + * energy and time. The model value is given by * - * @param[in] photon Incident photon (ignored). - * @return Model value (units of sr^-1). + * \f[ + * M_\mathrm{S}(\vec{p}|E,t) = \frac{v}{4\pi} + * \f] * - * Evaluates the spatial part for an isotropic source model and sets the - * parameter gradient. By definition, the model value and gradient are - * independent from the sky direction. + * where \f$v\f$ is the value() parameter, which is divided by the solid + * angle \f$4\pi\f$ of the celestial sphere. The value() gradient is given + * by + * + * \f[ + * \frac{\partial M_\mathrm{S}(\vec{p}|E,t)}{\partial v_v} = + * \frac{v_s}{4\pi} + * \f] + * + * where \f$v=v_v v_s\f$ is the factorisation of the value() parameter into + * a value \f$v_v\f$ and scale \f$v_s\f$ term. ***************************************************************************/ double GModelSpatialDiffuseConst::eval_gradients(const GPhoton& photon) const { @@ -283,41 +322,91 @@ double GModelSpatialDiffuseConst::eval_gradients(const GPhoton& photon) const /***********************************************************************//** * @brief Return MC sky direction * - * @param[in] energy Photon energy (ignored). - * @param[in] time Photon arrival time (ignored). + * @param[in] energy Photon energy. + * @param[in] time Photon arrival time. * @param[in,out] ran Random number generator. * @return Sky direction. * - * Returns an arbitrary position on the celestial sphere. This position is - * independent of event @p energy and arrival @p time. + * Returns an arbitrary direction on the celestial sphere within a simulation + * cone region that has been defined by a previous call to the mc_norm() + * method (if no such call was made the entire sky will be assumed as the + * simulation cone). The sky direction is independent of event @p energy and + * arrival @p time. ***************************************************************************/ GSkyDir GModelSpatialDiffuseConst::mc(const GEnergy& energy, const GTime& time, GRan& ran) const { - // Simulate Right Ascension and Declination - double ra = gammalib::twopi * ran.uniform(); - double dec = std::asin(2.0 * ran.uniform() - 1.0); + // Simulate offset from simulation cone centre + double theta = std::acos(1.0 - ran.uniform() * (1.0 - m_mc_cos_radius)) * + gammalib::rad2deg; + double phi = 360.0 * ran.uniform(); - // Set sky direction - GSkyDir dir; - dir.radec(ra, dec); + // Initialise sky direction with simulation cone centre + GSkyDir dir(m_mc_centre); + // Rotate sky direction by offset angles (phi, theta) + dir.rotate_deg(phi, theta); + // Return sky direction return dir; } /***********************************************************************//** + * @brief Return normalization of constant diffuse source for Monte Carlo + * simulations + * + * @param[in] dir Centre of simulation cone. + * @param[in] radius Radius of simulation cone (deg). + * @return Model normalization. + * + * @exception GException::invalid_argument + * Radius of simulation cone not coprised in [0,180] degrees. + * + * Returns the normalization of a constant diffuse source. The normalization + * is given by the product between the fraction \f$f\f$ of the sky covered + * by the simulation cone + * + * \f[ + * f = \frac{2\pi \left( 1-\cos(r) \right)}{4\pi} + * \f] + * + * (where \f$r\f$ is the radius of the simulation cone) and the model + * parameter value(). The normalization only depends of the radius of the + * simulation cone and is invariant to its centre. + ***************************************************************************/ +double GModelSpatialDiffuseConst::mc_norm(const GSkyDir& dir, + const double& radius) const +{ + // Throw exception if the radius is invalid + if (radius < 0.0 || radius > 180.0) { + std::string msg = "Invalid simulation cone radius "+ + gammalib::str(radius)+" specified. Please provide " + "a value comprised within 0 and 180 degrees."; + throw GException::invalid_argument(G_MC_NORM, msg); + } + + // Store cosine of simulation cone radius and cone centre + m_mc_centre = dir; + m_mc_cos_radius = std::cos(radius * gammalib::deg2rad); + + // Compute fraction of sky covered by the simulation cone + double fraction = 0.5 * (1.0 - m_mc_cos_radius); + + // Multiply by normalization value of model + double norm = fraction * value(); + + // Return normalization + return (norm); +} + + +/***********************************************************************//** * @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 isotropic source model information from an XML element. The XML * element is expected to have the following format: * @@ -331,23 +420,11 @@ void GModelSpatialDiffuseConst::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, - "Isotropic source model requires exactly 1 parameter."); - } - - // Get pointer on model parameter - const GXmlElement* par = xml.element("parameter", 0); + // Get value parameter + const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Value"); - // Get value - if (par->attribute("name") == "Value") { - m_value.read(*par); - } - else { - throw GException::model_invalid_parnames(G_READ, xml, - "Isotropic source model requires \"Value\" parameter."); - } + // Read value parameter + m_value.read(*par); // Return return; @@ -361,10 +438,6 @@ void GModelSpatialDiffuseConst::read(const GXmlElement& xml) * * @exception GException::model_invalid_spatial * Existing XML element is not of type "ConstantValue" - * @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 isotropic source model information into an XML element. The XML * element will have the following format: @@ -387,28 +460,11 @@ void GModelSpatialDiffuseConst::write(GXmlElement& xml) const "Spatial model is not of type \"ConstantValue\"."); } - // If XML element has 0 nodes then append 1 parameter node - if (xml.elements() == 0) { - xml.append(GXmlElement("parameter name=\"Value\"")); - } - - // Verify that XML element has exactly 1 parameter - if (xml.elements() != 1 || xml.elements("parameter") != 1) { - throw GException::model_invalid_parnum(G_WRITE, xml, - "Isotropic source model requires exactly 1 parameter."); - } + // Get or create Value parameter + GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "Value"); - // Get pointers on both model parameters - GXmlElement* par = xml.element("parameter", 0); - - // Set or update parameter - if (par->attribute("name") == "Value") { - m_value.write(*par); - } - else { - throw GException::model_invalid_parnames(G_WRITE, xml, - "Isotropic source model requires \"Value\" parameter."); - } + // Write parameter + m_value.write(*par); // Return return; @@ -467,6 +523,10 @@ void GModelSpatialDiffuseConst::init_members(void) m_value.gradient(0.0); m_value.has_grad(true); + // Initialise other members + m_mc_centre.clear(); + m_mc_cos_radius = -1.0; // Set simulation cone to allsky + // Set parameter pointer(s) m_pars.clear(); m_pars.push_back(&m_value); @@ -484,7 +544,9 @@ void GModelSpatialDiffuseConst::init_members(void) void GModelSpatialDiffuseConst::copy_members(const GModelSpatialDiffuseConst& model) { // Copy members - m_value = model.m_value; + m_value = model.m_value; + m_mc_centre = model.m_mc_centre; + m_mc_cos_radius = model.m_mc_cos_radius; // Set parameter pointer(s) m_pars.clear(); -- libgit2 0.21.4