Monte Carlo for empty maps or map cubes no longer throws an exception (#1915)

Jürgen Knödlseder
1 parent f2a928d4
Showing 7 changed files with 59 additions and 13 deletions
ChangeLog
NEWS
include/GException.hpp
src/model/GModelSky.cpp
src/model/GModelSpatialDiffuseCube.cpp
src/model/GModelSpatialDiffuseMap.cpp
src/support/GException.cpp
@@ -3,7 +3,7 @@
         * Version 1.2.0 released
         =========================
  
-        Improve exception message in GModelSpatialDiffuseCube::mc (#1915)
+        Monte Carlo for empty maps or map cubes no longer throws an exception (#1915)
         Add GModelTemporalFunc call (#1914)
         Add GSkyMap FITS HDU constructor
         Add leap second for 1-1-2017 (#1887)
@@ -175,8 +175,10 @@ function for the simulation of transient sources. A unit test has been
 added and the mc() and eval() methods have been validated using a dedicated
 test script (#1914).
  
-Improve the exception message for empty map cubes in GModelSpatialDiffuseCube::mc
-(#1915).
+The GModelSpatialDiffuseMap::mc and GModelSpatialDiffuseCube::mc methods
+now throw an GException::invalid_return_value exception that is catched
+in the GModelSky::mc method and leads to the skipping of the concerned
+photon. This allows handling of empty diffuse maps or map cubes (#1915).
  
  
 9. Numerics module
 /***************************************************************************
  *                    GException.hpp - Exception handler                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2006-2015 by Juergen Knoedlseder                         *
+ *  copyright (C) 2006-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -82,6 +82,13 @@ public:
                          const std::string& message);
     };
  
+    // Invalid return value
+    class invalid_return_value : public GExceptionHandler {
+    public:
+        invalid_return_value(const std::string& origin,
+                             const std::string& message);
+    };
+
     // Out of range
     class out_of_range : public GExceptionHandler {
     public:
@@ -942,8 +942,15 @@ GPhotons GModelSky::mc(const double&amp; area,
                 std::cout << "    Energy=" << photon.energy() << std::endl;
                 #endif
  
-                // Set incident photon direction
-                photon.dir(m_spatial->mc(photon.energy(), photon.time(), ran));
+                // Set incident photon direction. If an invalid_return_value
+                // exception occurs the sky direction returned by the spatial
+                // Monte Carlo method is invalid and the photon is skipped.
+                try {
+                    photon.dir(m_spatial->mc(photon.energy(), photon.time(), ran));
+                }
+                catch (GException::invalid_return_value) {
+                    continue;
+                }
  
                 // Debug option: dump direction
                 #if defined(G_DUMP_MC_DETAIL)
@@ -363,11 +363,15 @@ double GModelSpatialDiffuseCube::eval(const GPhoton&amp; photon,
  * @exception GException::invalid_value
  *            No map cube defined.
  *            No energy boundaries defined.
- *            Simulation cone not defined or does not overlap with map cube.
+ * @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.
+ * method to determine the sky direction. If no sky direction could be
+ * determined, the method throws an GException::invalid_return_value
+ * exception.
  ***************************************************************************/
 GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy& energy,
                                      const GTime&   time,
@@ -430,7 +434,7 @@ GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy&amp; energy,
                   "restrict the simulated energies to values where the map "
                   "cube is non-zero within the simulation cone.";
         }
-        throw GException::invalid_value(G_MC, msg);
+        throw GException::invalid_return_value(G_MC, msg);
     }
  
     // Allocate sky direction
@@ -350,12 +350,15 @@ double GModelSpatialDiffuseMap::eval(const GPhoton&amp; photon,
  * @return Sky direction.
  *
  * @exception GException::invalid_value
- *            Simulation cone not defined or does not overlap with sky map.
  *            No sky map defined.
+ * @exception GException::invalid_return_value
+ *            Simulation cone not defined, does not overlap with map or map
+ *            is empty.
  *
  * Returns a random sky direction according to the intensity distribution of
  * the model sky map. The method uses a rejection method to determine the sky
- * direction.
+ * direction. If no sky direction could be determined, the method throws an
+ * GException::invalid_return_value exception.
  ***************************************************************************/
 GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy& energy,
                                     const GTime&   time,
@@ -368,7 +371,7 @@ GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy&amp; energy,
         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);
+        throw GException::invalid_return_value(G_MC, msg);
     }
  
     // Determine number of skymap pixels
 /***************************************************************************
  *                   GException.cpp  -  exception handlers                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2006-2015 by Juergen Knoedlseder                         *
+ *  copyright (C) 2006-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -89,6 +89,29 @@ GException::invalid_argument::invalid_argument(const std::string&amp; origin,
  
  
 /***********************************************************************//**
+ * @brief Invalid return value
+ *
+ * @param[in] origin Name of method that has thrown the exception.
+ * @param[in] message Error message.
+ ***************************************************************************/
+GException::invalid_return_value::invalid_return_value(const std::string& origin,
+                                                       const std::string& message)
+{
+    // Set origin
+    m_origin  = origin;
+
+    // Set message string
+    m_message = "Invalid return value.";
+    if (message.length() > 0) {
+        m_message += (" " + message);
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Index is out of range [0,elements-1]
  *
  * @param[in] origin Method throwing the exception.