Fix consistency problem in GCTABackground3D class

Jürgen Knödlseder
1 parent c31806e6
Showing 7 changed files with 453 additions and 264 deletions
ChangeLog
NEWS
inst/cta/include/GCTABackground3D.hpp
inst/cta/src/GCTABackground3D.cpp
inst/cta/src/GCTAModelIrfBackground.cpp
inst/cta/src/GCTAResponseIrf.cpp
src/sky/GSkyMap.cpp
-2016-07-12  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2016-07-18  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 1.1.0 released
         =========================
+        Fix consistency problem in GCTABackground3D (#1815)
         Models accessing file now support relative file names (#1808)
         Add GXmlElement and GXmlDocument filename() methods (#1808)
         Move node parent from GXmlElement to GXmlNode class (#1808)
-        Correct Makefile for example executables (#1806)
+        Fix Makefile for example executables (#1806)
         Include examples in unit tests (#1803)
         Handle additional columns in GArf class
         GCTAAeff2D no longer extrapolates beyond response table (#1784)
@@ -56,8 +57,8 @@
         Add GFilename::fullname() method
         Add GCTAEventList::remove() method
         Handle also additional columns in GCTAEventList::append() method
-        Correct GFitsTableCol::insert() method for empty columns
-        Correct GCTAEdisp2D::compute_ebounds_xxx() methods (#1662)
+        Fix GFitsTableCol::insert() method for empty columns
+        Fix GCTAEdisp2D::compute_ebounds_xxx() methods (#1662)
         Use GFilename in GFits to store filename (#1649)
         Load CTA GTI extension from the HDU name provided by the DSS keyword (#1598)
         Support Fermi LAT Pass 8 response functions (#1517)
@@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
 Author(s) : Juergen Knoedlseder
-12 July 2016
+18 July 2016
 1. Introduction
@@ -440,6 +440,10 @@ handled by the weight array (#1789).
 The GCTAAeff2D class no longer extrapolates the effective area beyond the boundaries
 of the response table, but sets the effective area values to zero (#1784).
+The GCTABackground3D class now assures that the mc() method covers the same DETX and DETY
+range than the access operator. This fixes a consistency problem between Monte Carlo
+sampling and function evaluation (#1815).
+
 21. COMPTEL interface
 ---------------------
@@ -83,33 +83,39 @@ public:
 private:
     // Methods
-    void init_members(void);
-    void copy_members(const GCTABackground3D& bgd);
-    void free_members(void);
-    void init_mc_cache(void) const;
+    void   init_members(void);
+    void   copy_members(const GCTABackground3D& bgd);
+    void   free_members(void);
+    void   set_limits(void);
+    int    index(const int& idetx, const int& idety, const int& iebin) const;
+    void   init_mc_cache(void) const;
+    void   init_mc_max_rate(void) const;
+    double solid_angle(const double& detx1, const double& dety1,
+                       const double& detx2, const double& dety2,
+                       const double& detx3, const double& dety3) const;
     // Members
     GFilename         m_filename;    //!< Name of background response file
     GCTAResponseTable m_background;  //!< Background response table
-    double            m_mc_max_bin;  //!< Maximum spatial binsize for MC
-    double            m_mc_max_logE; //!< Maximum log energy binsize for MC
     int               m_inx_detx;    //!< DETX index
     int               m_inx_dety;    //!< DETY index
     int               m_inx_energy;  //!< Energy index
     int               m_inx_bgd;     //!< Background index
+    int               m_num_detx;    //!< Number of DETX bins
+    int               m_num_dety;    //!< Number of DETY bins
+    int               m_num_energy;  //!< Number of energy bins
+    int               m_num[3];      //!< Array of number of bins
+    double            m_detx_min;    //!< DETX minimum (radians)
+    double            m_detx_max;    //!< DETX maximum (radians)
+    double            m_dety_min;    //!< DETY minimum (radians)
+    double            m_dety_max;    //!< DETY maximum (radians)
+    double            m_logE_min;    //!< Log10(E/TeV) minimum
+    double            m_logE_max;    //!< Log10(E/TeV) maximum
     // Monte Carlo cache
-    mutable std::vector<double> m_mc_max;      //!< Maximum rate
-    mutable GModelSpectralNodes m_mc_spectrum; //!< Response cube spectrum
-    mutable double              m_mc_detx_min; //!< DETX minimum (deg)
-    mutable double              m_mc_detx_max; //!< DETX maximum (deg)
-    mutable double              m_mc_detx_bin; //!< DETX binsize (deg)
-    mutable double              m_mc_dety_min; //!< DETY minimum (deg)
-    mutable double              m_mc_dety_max; //!< DETY maximum (deg)
-    mutable double              m_mc_dety_bin; //!< DETY binsize (deg)
-    mutable double              m_mc_logE_min; //!< log10 energy minimum (TeV)
-    mutable double              m_mc_logE_max; //!< log10 energy maximum (TeV)
-    mutable double              m_mc_logE_bin; //!< log10 energy binsize (TeV)
+    mutable std::vector<double> m_mc_max;                //!< Maximum background rate
+    mutable GModelSpectralNodes m_mc_spectrum;           //!< Response cube spectrum
+    mutable double              m_mc_one_minus_costheta; //!< 1-cos(theta_max)
 };
@@ -174,27 +174,28 @@ GCTABackground3D&amp; GCTABackground3D::operator=(const GCTABackground3D&amp; bgd)
  * @brief Return background rate in units of events/s/MeV/sr
  *
  * @param[in] logE Log10 of the true photon energy (TeV).
- * @param[in] detx Tangential coord in nominal sys (radians).
- * @param[in] dety Tangential coord in nominal sys (radians).
+ * @param[in] detx Tangential X coordinate in nominal system (radians).
+ * @param[in] dety Tangential Y coordinate in nominal system (radians).
  *
  * Returns the background rate in units of events/s/MeV/sr for a given energy
- * and detector coordinates. The method assures that the background rate
- * never becomes negative.
+ * and detector coordinates. The operator interpolates linearly in DETX and
+ * DETY, and logarithmically in energy.
  *
- * The method interpolates linearly in DETX and DETY, and logarithmically
- * in energy.
- *
- * If the background model is not valid the operator returns 0.
+ * The operator assures that the background rate never becomes negative. For
+ * invalid background models or detector coordinates outside the range
+ * covered by the response cube, the operator returns a rate of zero.
  ***************************************************************************/
 double GCTABackground3D::operator()(const double& logE, 
                                     const double& detx, 
                                     const double& dety) const
 {
-    // Initialise rate
-    double rate(0.0);
+    // Initialise background rate
+    double rate = 0.0;
-    // Continue only if background model is valid
-    if (is_valid()) {
+    // Continue only if background model is valid and if the DETX and DETY
+    // values are in validity range
+    if (is_valid() && (detx >= m_detx_min) && (detx <= m_detx_max) &&
+                      (dety >= m_dety_min) && (dety <= m_dety_max)) {
         // Retrieve references to node arrays
         const GNodeArray& detx_nodes   = m_background.axis_nodes(m_inx_detx);
@@ -206,46 +207,43 @@ double GCTABackground3D::operator()(const double&amp; logE,
         dety_nodes.set_value(dety);
         energy_nodes.set_value(logE);
-        // Compute offsets of DETY in DETX-DETY plane
-        int size1        = m_background.axis_bins(m_inx_detx);
-        int offset_left  = dety_nodes.inx_left()  * size1;
-        int offset_right = dety_nodes.inx_right() * size1;
-
-        // Set indices for bi-linear interpolation in DETX-DETY plane
-        int inx_ll = detx_nodes.inx_left()  + offset_left;
-        int inx_lr = detx_nodes.inx_left()  + offset_right;
-        int inx_rl = detx_nodes.inx_right() + offset_left;
-        int inx_rr = detx_nodes.inx_right() + offset_right;
-
         // Set weighting factors for bi-linear interpolation in DETX-DETY plane
         double wgt_ll = detx_nodes.wgt_left()  * dety_nodes.wgt_left();
         double wgt_lr = detx_nodes.wgt_left()  * dety_nodes.wgt_right();
         double wgt_rl = detx_nodes.wgt_right() * dety_nodes.wgt_left();
         double wgt_rr = detx_nodes.wgt_right() * dety_nodes.wgt_right();
-        // Set indices for energy interpolation
-        int inx_emin = energy_nodes.inx_left();
-        int inx_emax = energy_nodes.inx_right();
-
         // Set weighting factors for energy interpolation
         double wgt_emin = energy_nodes.wgt_left();
         double wgt_emax = energy_nodes.wgt_right();
-        // Compute offsets in energy dimension
-        int npixels     = m_background.axis_bins(m_inx_detx) *
-                          m_background.axis_bins(m_inx_dety);
-        int offset_emin = inx_emin * npixels;
-        int offset_emax = inx_emax * npixels;
-
-        // Bi-linear interpolation the rates in both energy layers
-        double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll + offset_emin) +
-                           wgt_lr * m_background(m_inx_bgd, inx_lr + offset_emin) +
-                           wgt_rl * m_background(m_inx_bgd, inx_rl + offset_emin) +
-                           wgt_rr * m_background(m_inx_bgd, inx_rr + offset_emin);
-        double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll + offset_emax) +
-                           wgt_lr * m_background(m_inx_bgd, inx_lr + offset_emax) +
-                           wgt_rl * m_background(m_inx_bgd, inx_rl + offset_emax) +
-                           wgt_rr * m_background(m_inx_bgd, inx_rr + offset_emax);
+        // Get bin indices
+        int inx_ll_emin = index(detx_nodes.inx_left(),  dety_nodes.inx_left(),
+                                energy_nodes.inx_left());
+        int inx_lr_emin = index(detx_nodes.inx_left(),  dety_nodes.inx_right(),
+                                energy_nodes.inx_left());
+        int inx_rl_emin = index(detx_nodes.inx_right(), dety_nodes.inx_left(),
+                                energy_nodes.inx_left());
+        int inx_rr_emin = index(detx_nodes.inx_right(), dety_nodes.inx_right(),
+                                energy_nodes.inx_left());
+        int inx_ll_emax = index(detx_nodes.inx_left(),  dety_nodes.inx_left(),
+                                energy_nodes.inx_right());
+        int inx_lr_emax = index(detx_nodes.inx_left(),  dety_nodes.inx_right(),
+                                energy_nodes.inx_right());
+        int inx_rl_emax = index(detx_nodes.inx_right(), dety_nodes.inx_left(),
+                                energy_nodes.inx_right());
+        int inx_rr_emax = index(detx_nodes.inx_right(), dety_nodes.inx_right(),
+                                energy_nodes.inx_right());
+
+        // Bi-linearly interpolate the rates in both energy layers
+        double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll_emin) +
+                           wgt_lr * m_background(m_inx_bgd, inx_lr_emin) +
+                           wgt_rl * m_background(m_inx_bgd, inx_rl_emin) +
+                           wgt_rr * m_background(m_inx_bgd, inx_rr_emin);
+        double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll_emax) +
+                           wgt_lr * m_background(m_inx_bgd, inx_lr_emax) +
+                           wgt_rl * m_background(m_inx_bgd, inx_rl_emax) +
+                           wgt_rr * m_background(m_inx_bgd, inx_rr_emax);
         // If both rates are positive then perform a logarithmic interpolation
         // in energy
@@ -338,12 +336,6 @@ void GCTABackground3D::read(const GFitsTable&amp; table)
     // Read background table
     m_background.read(table);
-    // Get mandatory indices (throw exception if not found)
-    m_inx_detx   = m_background.axis("DETX");
-    m_inx_dety   = m_background.axis("DETY");
-    m_inx_energy = m_background.axis("ENERG");
-    m_inx_bgd    = m_background.table("BGD");
-
     // Throw an exception if the table is not three-dimensional
     if (m_background.axes() != 3) {
         std::string msg = "Expected three-dimensional background "
@@ -354,6 +346,22 @@ void GCTABackground3D::read(const GFitsTable&amp; table)
         throw GException::invalid_value(G_READ, msg);
     }
+    // Get mandatory indices (throw exception if not found)
+    m_inx_detx   = m_background.axis("DETX");
+    m_inx_dety   = m_background.axis("DETY");
+    m_inx_energy = m_background.axis("ENERG");
+    m_inx_bgd    = m_background.table("BGD");
+
+    // Get axes dimensions
+    m_num_detx   = m_background.axis_bins(m_inx_detx);
+    m_num_dety   = m_background.axis_bins(m_inx_dety);
+    m_num_energy = m_background.axis_bins(m_inx_energy);
+
+    // Set dimension array
+    m_num[m_inx_detx]   = m_num_detx;
+    m_num[m_inx_dety]   = m_num_dety;
+    m_num[m_inx_energy] = m_num_energy;
+
     // Set DETX and DETY axis to radians
     m_background.axis_radians(m_inx_detx);
     m_background.axis_radians(m_inx_dety);
@@ -361,6 +369,9 @@ void GCTABackground3D::read(const GFitsTable&amp; table)
     // Set energy axis to logarithmic scale
     m_background.axis_log10(m_inx_energy);
+    // Set DETX and DETY limits
+    set_limits();
+
     // Return
     return;
 }
@@ -484,7 +495,7 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy&amp; energy,
                                  const GTime&   time,
                                  GRan&          ran) const
 {
-    // Initialise Monte Carlo Cache
+    // If Monte Carlo has not been initialised then initialise the cache now
     if (m_mc_max.empty()) {
         init_mc_cache();
     }
@@ -508,35 +519,35 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy&amp; energy,
         int inx_left  = energy_nodes.inx_left();
         int inx_right = energy_nodes.inx_right();
-        // Get maximum background rates
-        double max_rate_left  = m_mc_max[inx_left];
-        double max_rate_right = m_mc_max[inx_right];
-
-        // Compute maximum background rate
-        double max_rate = 0.0;
-        if (max_rate_left > 0.0 && max_rate_right > 0.0) {
-            double wgt_left  = energy_nodes.wgt_left();
-            double wgt_right = energy_nodes.wgt_right();
-            max_rate         = std::exp(wgt_left  * std::log(max_rate_left) +
-                                        wgt_right * std::log(max_rate_right));
-        }
-        else {
-            max_rate = 0.0;
+        // Get maximum background rate as the maximum of the left and the right
+        // energy node
+        double max_rate = m_mc_max[inx_left];
+        if (m_mc_max[inx_right] > max_rate) {
+            max_rate = m_mc_max[inx_right];
         }
-        // Get instrument direction
+        // Get instrument direction using a rejection method. This assures
+        // that the Monte Carlo sample follows the model distribution
         while (true) {
-            // Get randomized detx and dety (radians)
-            double detx = (m_mc_detx_min +
-                           ran.uniform() * (m_mc_detx_max - m_mc_detx_min)) *
-                          gammalib::deg2rad;
-            double dety = (m_mc_dety_min +
-                           ran.uniform() * (m_mc_dety_max - m_mc_dety_min)) *
-                          gammalib::deg2rad;
+            // Get randomized detx and dety in radians. If the direction falls
+            // outside the validity range then get another pair of detx and
+            // dety values.
+            double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_costheta);
+            double phi   = gammalib::twopi * ran.uniform();
+            double detx  = theta * std::cos(phi);
+            double dety  = theta * std::sin(phi);
+            if ((detx < m_detx_min) || (detx > m_detx_max) ||
+                (dety < m_dety_min) || (dety > m_dety_max)) {
+                continue;
+            }
-            // Get background rate for these coordinates
+            // Get background rate for these coordinates. If the background
+            // rate is zero then get other coordinates.
             double value = this->operator()(logE, detx, dety);
+            if (value <= 0.0) {
+                continue;
+            }
             // Get uniform random number
             double uniform = ran.uniform() * max_rate;
@@ -581,37 +592,30 @@ std::string GCTABackground3D::print(const GChatter&amp; chatter) const
         if (is_valid()) {
             // Compute DETX boundaries in deg
-            double detx_min = m_background.axis_lo(m_inx_detx,0);
-            double detx_max = m_background.axis_hi(m_inx_detx,
-                              m_background.axis_bins(m_inx_detx)-1);
+            double detx_min = m_background.axis_lo(m_inx_detx, 0);
+            double detx_max = m_background.axis_hi(m_inx_detx, m_num_detx-1);
             // Compute DETY boundaries in deg
             double dety_min = m_background.axis_lo(m_inx_dety,0);
-            double dety_max = m_background.axis_hi(m_inx_dety,
-                              m_background.axis_bins(m_inx_dety)-1);
+            double dety_max = m_background.axis_hi(m_inx_dety, m_num_dety-1);
             // Compute energy boundaries in TeV
-            double emin = m_background.axis_lo(m_inx_energy,0);
-            double emax = m_background.axis_hi(m_inx_energy,
-                          m_background.axis_bins(m_inx_energy)-1);
+            double emin = m_background.axis_lo(m_inx_energy, 0);
+            double emax = m_background.axis_hi(m_inx_energy, m_num_energy-1);
             // Append information
             result.append("\n"+gammalib::parformat("Number of DETX bins") +
-                          gammalib::str(m_background.axis_bins(m_inx_detx)));
+                          gammalib::str(m_num_detx));
             result.append("\n"+gammalib::parformat("Number of DETY bins") +
-                          gammalib::str(m_background.axis_bins(m_inx_dety)));
+                          gammalib::str(m_num_dety));
             result.append("\n"+gammalib::parformat("Number of energy bins") +
-                          gammalib::str(m_background.axis_bins(m_inx_energy)));
+                          gammalib::str(m_num_energy));
             result.append("\n"+gammalib::parformat("DETX range"));
             result.append(gammalib::str(detx_min)+" - "+gammalib::str(detx_max)+" deg");
             result.append("\n"+gammalib::parformat("DETX range"));
             result.append(gammalib::str(dety_min)+" - "+gammalib::str(dety_max)+" deg");
             result.append("\n"+gammalib::parformat("Energy range"));
             result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
-            result.append("\n"+gammalib::parformat("Maximum bin size for MC"));
-            result.append(gammalib::str(m_mc_max_bin)+" deg");
-            result.append("\n"+gammalib::parformat("Maximum logE step for MC"));
-            result.append(gammalib::str(m_mc_max_logE)+"^10 TeV");
         } // endif: there were 3 axis
@@ -646,25 +650,27 @@ void GCTABackground3D::init_members(void)
     // Initialise members
     m_filename.clear();
     m_background.clear();
-    m_mc_max_bin  = 0.05;  //!< Spatial binning not worse than 0.05 deg
-    m_mc_max_logE = 0.02;  //!< Spectral binning not worse than 0.02^10 TeV
-    m_inx_detx    = 0;
-    m_inx_dety    = 1;
-    m_inx_energy  = 2;
-    m_inx_bgd     = 0;
+    m_inx_detx   = 0;
+    m_inx_dety   = 1;
+    m_inx_energy = 2;
+    m_inx_bgd    = 0;
+    m_num_detx   = 0;
+    m_num_dety   = 0;
+    m_num_energy = 0;
+    m_num[0]     = 0;
+    m_num[1]     = 0;
+    m_num[2]     = 0;
+    m_detx_min   = 0.0;
+    m_detx_max   = 0.0;
+    m_dety_min   = 0.0;
+    m_dety_max   = 0.0;
+    m_logE_min   = 0.0;
+    m_logE_max   = 0.0;
     // Initialise MC cache
     m_mc_max.clear();
     m_mc_spectrum.clear();
-    m_mc_detx_min = 0.0;
-    m_mc_detx_max = 0.0;
-    m_mc_detx_bin = 0.0;
-    m_mc_dety_min = 0.0;
-    m_mc_dety_max = 0.0;
-    m_mc_dety_bin = 0.0;
-    m_mc_logE_min = 0.0;
-    m_mc_logE_max = 0.0;
-    m_mc_logE_bin = 0.0;
+    m_mc_one_minus_costheta = 0.0;
     // Return
     return;
@@ -679,27 +685,29 @@ void GCTABackground3D::init_members(void)
 void GCTABackground3D::copy_members(const GCTABackground3D& bgd)
 {
     // Copy members
-    m_filename    = bgd.m_filename;
-    m_background  = bgd.m_background;
-    m_mc_max_bin  = bgd.m_mc_max_bin;
-    m_mc_max_logE = bgd.m_mc_max_logE;
-    m_inx_detx    = bgd.m_inx_detx;
-    m_inx_dety    = bgd.m_inx_dety;
-    m_inx_energy  = bgd.m_inx_energy;
-    m_inx_bgd     = bgd.m_inx_bgd;
+    m_filename   = bgd.m_filename;
+    m_background = bgd.m_background;
+    m_inx_detx   = bgd.m_inx_detx;
+    m_inx_dety   = bgd.m_inx_dety;
+    m_inx_energy = bgd.m_inx_energy;
+    m_inx_bgd    = bgd.m_inx_bgd;
+    m_num_detx   = bgd.m_num_detx;
+    m_num_dety   = bgd.m_num_dety;
+    m_num_energy = bgd.m_num_energy;
+    m_num[0]     = bgd.m_num[0];
+    m_num[1]     = bgd.m_num[1];
+    m_num[2]     = bgd.m_num[2];
+    m_detx_min   = bgd.m_detx_min;
+    m_detx_max   = bgd.m_detx_max;
+    m_dety_min   = bgd.m_dety_min;
+    m_dety_max   = bgd.m_dety_max;
+    m_logE_min   = bgd.m_logE_min;
+    m_logE_max   = bgd.m_logE_max;
     // Copy MC cache
-    m_mc_max      = bgd.m_mc_max;
-    m_mc_spectrum = bgd.m_mc_spectrum;
-    m_mc_detx_min = bgd.m_mc_detx_min;
-    m_mc_detx_max = bgd.m_mc_detx_max;
-    m_mc_detx_bin = bgd.m_mc_detx_bin;
-    m_mc_dety_min = bgd.m_mc_dety_min;
-    m_mc_dety_max = bgd.m_mc_dety_max;
-    m_mc_dety_bin = bgd.m_mc_dety_bin;
-    m_mc_logE_min = bgd.m_mc_logE_min;
-    m_mc_logE_max = bgd.m_mc_logE_max;
-    m_mc_logE_bin = bgd.m_mc_logE_bin;
+    m_mc_max                = bgd.m_mc_max;
+    m_mc_spectrum           = bgd.m_mc_spectrum;
+    m_mc_one_minus_costheta = bgd.m_mc_one_minus_costheta;
     // Return
     return;
@@ -707,6 +715,66 @@ void GCTABackground3D::copy_members(const GCTABackground3D&amp; bgd)
 /***********************************************************************//**
+ * @brief Set DETX, DETY and log10(E/TeV) limits
+ *
+ * Set the DETX, DETY and log10(E/TeV) limits.
+ ***************************************************************************/
+void GCTABackground3D::set_limits(void)
+{
+    // Compute DETX boundaries in radians
+    m_detx_min = m_background.axis_lo(m_inx_detx, 0) *
+                 gammalib::deg2rad;
+    m_detx_max = m_background.axis_hi(m_inx_detx, m_num_detx-1) *
+                 gammalib::deg2rad;
+
+    // Compute DETY boundaries in radians
+    m_dety_min = m_background.axis_lo(m_inx_dety, 0) *
+                 gammalib::deg2rad;
+    m_dety_max = m_background.axis_hi(m_inx_dety, m_num_dety-1) *
+                 gammalib::deg2rad;
+
+    // Compute energy boundaries in log10(TeV)
+    GEnergy emin(m_background.axis_lo(m_inx_energy, 0),
+                 m_background.axis_lo_unit(m_inx_energy));
+    GEnergy emax(m_background.axis_hi(m_inx_energy, m_num_energy-1),
+                 m_background.axis_hi_unit(m_inx_energy));
+    m_logE_min = emin.log10TeV();
+    m_logE_max = emax.log10TeV();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return background rate bin index
+ *
+ * @param[in] idetx DETX index.
+ * @param[in] idety DETY index.
+ * @param[in] iebin Energy index.
+ * @return Background rate bin index.
+ *
+ * Returns the background rate bin index independent of the ordering.
+ ***************************************************************************/
+int GCTABackground3D::index(const int& idetx,
+                            const int& idety,
+                            const int& iebin) const
+{
+    // Set index arrays
+    int inx[3];
+    inx[m_inx_detx]   = idetx;
+    inx[m_inx_dety]   = idety;
+    inx[m_inx_energy] = iebin;
+
+    // Compute index
+    int index = inx[0] + (inx[1] + inx[2] * m_num[1]) * m_num[0];
+
+    // Return index
+    return index;
+}
+
+
+/***********************************************************************//**
  * @brief Delete class members
  ***************************************************************************/
 void GCTABackground3D::free_members(void)
@@ -719,15 +787,7 @@ void GCTABackground3D::free_members(void)
 /***********************************************************************//**
  * @brief Initialise Monte Carlo cache
  *
- * @exception GException::invalid_value
- *            No valid background model defined.
- *
  * Initialises the cache for Monte Carlo sampling.
- * 
- * @todo Verify assumption made about the solid angles of the response table
- *       elements.
- * @todo Implement method that computes the "flux" in a background map
- *       pixel (inspired from GSkyMap::flux())
  ***************************************************************************/
 void GCTABackground3D::init_mc_cache(void) const
 {
@@ -735,125 +795,243 @@ void GCTABackground3D::init_mc_cache(void) const
     m_mc_max.clear();
     m_mc_spectrum.clear();
-    // Determine number of response cube pixels and maps
-    int nx    = m_background.axis_bins(m_inx_detx);
-    int ny    = m_background.axis_bins(m_inx_dety);
-    int npix  = nx * ny;
-    int nmaps = m_background.axis_bins(m_inx_energy);
-
-    // Compute DETX boundaries and binsize in degress
-    m_mc_detx_min = m_background.axis_lo(m_inx_detx,0);
-    m_mc_detx_max = m_background.axis_hi(m_inx_detx,nx-1);
-    m_mc_detx_bin = (m_mc_detx_max - m_mc_detx_min) / nx;
+    // Initialise maximum rate array
+    init_mc_max_rate();
+
+    // Compute 1-cos(theta_max)
+    double theta1 = std::sqrt(m_detx_min * m_detx_min + m_dety_min * m_dety_min);
+    double theta2 = std::sqrt(m_detx_min * m_detx_min + m_dety_max * m_dety_max);
+    double theta3 = std::sqrt(m_detx_max * m_detx_max + m_dety_min * m_dety_min);
+    double theta4 = std::sqrt(m_detx_max * m_detx_max + m_dety_max * m_dety_max);
+    double theta = theta1;
+    if (theta2 > theta) {
+        theta = theta2;
+    }
+    if (theta3 > theta) {
+        theta = theta3;
+    }
+    if (theta4 > theta) {
+        theta = theta4;
+    }
+    m_mc_one_minus_costheta = 1.0 - std::cos(theta);
-    // Compute DETY boundaries and binsize in degress
-    m_mc_dety_min = m_background.axis_lo(m_inx_dety,0);
-    m_mc_dety_max = m_background.axis_hi(m_inx_dety,ny-1);
-    m_mc_dety_bin = (m_mc_dety_max - m_mc_dety_min) / ny;
+    // Set number of energy bins. Make sure that we compute at least 100 energy
+    // bins
+    int nbins = m_num_energy * 10;
+    if (nbins < 100) {
+        nbins = 100;
+    }
-    // Compute energy boundaries and binsize in log10(TeV)
-    GEnergy emin(m_background.axis_lo(m_inx_energy,0),
-                 m_background.axis_lo_unit(m_inx_energy));
-    GEnergy emax(m_background.axis_hi(m_inx_energy,nmaps-1),
-                 m_background.axis_hi_unit(m_inx_energy));
-    m_mc_logE_min = emin.log10TeV();
-    m_mc_logE_max = emax.log10TeV();
-    m_mc_logE_bin = (m_mc_logE_max - m_mc_logE_min) / nmaps;
-
-    // Determine solid angle of the pixel. We assume here simply that we
-    // have square pixels of identical solid angle. It needs to be checked
-    // whether this is a valid assumption.
-    double solidangle = m_mc_detx_bin * gammalib::deg2rad *
-                        m_mc_dety_bin * gammalib::deg2rad;
-
-    // Compute flux normalisation factor. As we add-up 3 intensity we have
-    // do divide by 3, and as we add up 8 wedges we have to divide the
-    // solid angle by 8.
-    double flux_norm  = gammalib::onethird * solidangle / 8.0;
-
-    // Continue only if there are pixels and maps
-    if (npix > 0 && nmaps > 0) {
-
-        // Loop over all maps
-        for (int i = 0; i < nmaps; ++i) {
-
-            // Get logE
-            double logE = m_mc_logE_min + (i+0.5) * m_mc_logE_bin;
-
-            // Initialise cache with cumulative pixel fluxes and compute
-            // total flux in response table for normalization. Negative
-            // pixels are excluded from the cumulative map.
-            double total_flux = 0.0; // units: events/s/MeV
-            double max_rate   = 0.0; // units: events/s/sr/MeV
-            double xmin       = m_mc_detx_min * gammalib::deg2rad;
-            double xbin       = m_mc_detx_bin * gammalib::deg2rad;
-            double ymin       = m_mc_dety_min * gammalib::deg2rad;
-            double ybin       = m_mc_dety_bin * gammalib::deg2rad;
-            double dx         = 0.5 * xbin;
-            double dy         = 0.5 * ybin;
-            double detx       = xmin + 0.5 * xbin;
-            for (int ix = 0; ix < nx; ++ix, detx += xbin) {
-                double dety = ymin + 0.5 * ybin;
-                for (int iy = 0; iy < ny; ++iy, dety += ybin) {
-
-                    // Compute flux by summing the flux in 8 pixel wedges
-                    double i0    = (*this)(logE, detx,    dety);
-                    double i1    = (*this)(logE, detx-dx, dety-dy);
-                    double i2    = (*this)(logE, detx,    dety-dy);
-                    double i3    = (*this)(logE, detx+dx, dety-dy);
-                    double i4    = (*this)(logE, detx+dx, dety);
-                    double i5    = (*this)(logE, detx+dx, dety+dy);
-                    double i6    = (*this)(logE, detx,    dety+dy);
-                    double i7    = (*this)(logE, detx-dx, dety+dy);
-                    double i8    = (*this)(logE, detx-dx, dety);
-                    double flux1 = flux_norm * (i1 + i2 + i0);
-                    double flux2 = flux_norm * (i2 + i3 + i0);
-                    double flux3 = flux_norm * (i3 + i4 + i0);
-                    double flux4 = flux_norm * (i4 + i5 + i0);
-                    double flux5 = flux_norm * (i5 + i6 + i0);
-                    double flux6 = flux_norm * (i6 + i7 + i0);
-                    double flux7 = flux_norm * (i7 + i8 + i0);
-                    double flux8 = flux_norm * (i8 + i1 + i0);
-                    double flux  = (flux1 + flux2  + flux3  + flux4 +
-                                    flux5 + flux6  + flux7  + flux8);
-
-                    // Sum flux
-                    if (flux > 0.0) {
-                        total_flux += flux;
-                    }
-
-                    // Get maximum rate
-                    if (i0 > max_rate) {
-                        max_rate = i0;
-                    }
+    // Compute DETX and DETY binsize in radians and energy binsize in log10(TeV)
+    double detx_bin = (m_detx_max - m_detx_min) / m_num_detx;
+    double dety_bin = (m_dety_max - m_dety_min) / m_num_dety;
+    double logE_bin = (m_logE_max - m_logE_min) / (nbins+1);
+
+    // Loop over all energy bins
+    for (int i = 0; i < nbins; ++i) {
+
+        // Get logE of energy bin
+        double logE = m_logE_min + i * logE_bin;
+
+        // Initialise cache with cumulative pixel fluxes and compute total flux
+        // in response table for normalization. Negative pixels are excluded
+        // from the cumulative map.
+        double total_flux = 0.0; // units: events/s/MeV
+        double dx         = 0.5 * detx_bin;
+        double dy         = 0.5 * dety_bin;
+        double detx       = m_detx_min + dx;
+        for (int ix = 0; ix < m_num_detx; ++ix, detx += detx_bin) {
+            double dety = m_dety_min + dy;
+            for (int iy = 0; iy < m_num_dety; ++iy, dety += dety_bin) {
+
+                // Compute intensities
+                double i0 = (*this)(logE, detx,    dety);
+                double i1 = (*this)(logE, detx-dx, dety-dy);
+                double i2 = (*this)(logE, detx,    dety-dy);
+                double i3 = (*this)(logE, detx+dx, dety-dy);
+                double i4 = (*this)(logE, detx+dx, dety);
+                double i5 = (*this)(logE, detx+dx, dety+dy);
+                double i6 = (*this)(logE, detx,    dety+dy);
+                double i7 = (*this)(logE, detx-dx, dety+dy);
+                double i8 = (*this)(logE, detx-dx, dety);
+
+                // Compute solid angle of the 8 pixel wedges
+                double s1 = solid_angle(detx, dety, detx-dx, dety-dy, detx,    dety-dy);
+                double s2 = solid_angle(detx, dety, detx,    dety-dy, detx+dx, dety-dy);
+                double s3 = solid_angle(detx, dety, detx+dx, dety-dy, detx+dx, dety);
+                double s4 = solid_angle(detx, dety, detx+dx, dety,    detx+dx, dety+dy);
+                double s5 = solid_angle(detx, dety, detx+dx, dety+dy, detx,    dety+dy);
+                double s6 = solid_angle(detx, dety, detx,    dety+dy, detx-dx, dety+dy);
+                double s7 = solid_angle(detx, dety, detx-dx, dety+dy, detx-dx, dety);
+                double s8 = solid_angle(detx, dety, detx-dx, dety,    detx-dx, dety-dy);
+
+                // Compute flux by summing the flux in 8 pixel wedges
+                double flux1 = s1 * (i1 + i2 + i0);
+                double flux2 = s2 * (i2 + i3 + i0);
+                double flux3 = s3 * (i3 + i4 + i0);
+                double flux4 = s4 * (i4 + i5 + i0);
+                double flux5 = s5 * (i5 + i6 + i0);
+                double flux6 = s6 * (i6 + i7 + i0);
+                double flux7 = s7 * (i7 + i8 + i0);
+                double flux8 = s8 * (i8 + i1 + i0);
+                double flux  = (flux1 + flux2  + flux3  + flux4 +
+                                flux5 + flux6  + flux7  + flux8) / 3.0;
+
+                // Sum flux
+                if (flux > 0.0) {
+                    total_flux += flux;
                 }
-            }
-            // Append maximum rate
-            m_mc_max.push_back(max_rate);
+            } // endfor: loop over DETY
-            // Set energy value (unit independent)
-            GEnergy energy;
-            energy.log10(logE, m_background.axis_lo_unit(m_inx_energy));
+        } // endfor: loop over DETX
-            // Only append node if flux is positive
-            if (total_flux > 0.0) {
-                m_mc_spectrum.append(energy, total_flux);
-            }
+        // Set energy value (unit independent)
+        GEnergy energy;
+        energy.log10(logE, m_background.axis_lo_unit(m_inx_energy));
-            // Dump spectrum for debugging
-            #if defined(G_DEBUG_MC_INIT)
-            std::cout << "Energy=" << energy;
-            std::cout << " Rate_max=" << max_rate;
-            std::cout << " events/s/sr/MeV";
-            std::cout << " Rate_total=" << total_flux;
-            std::cout << " events/s/MeV" << std::endl;
-            #endif
+        // Only append node if flux is positive
+        if (total_flux > 0.0) {
+            m_mc_spectrum.append(energy, total_flux);
+        }
-        } // endfor: looped over all maps
+        // Dump spectrum for debugging
+        #if defined(G_DEBUG_MC_INIT)
+        std::cout << "Energy=" << energy;
+        std::cout << " Rate_total=" << total_flux;
+        std::cout << " events/s/MeV" << std::endl;
+        #endif
+
+    } // endfor: looped over energy bins
+
+    // Dump maximum rate
+    #if defined(G_DEBUG_MC_INIT)
+    for (int i = 0; i < m_mc_max.size(); ++i) {
+        double logE = m_logE_min + (i+0.5) * logE_bin;
+        GEnergy energy;
+        energy.log10(logE, m_background.axis_lo_unit(m_inx_energy));
+        std::cout << "Energy=" << energy;
+        std::cout << " Rate_max=" << max_rate;
+        std::cout << " events/s/sr/MeV";
+    }
+    #endif
-    } // endif: there were cube pixels and maps
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Initialise array of maximum background rate
+ *
+ * Initialises the array m_mc_max that holds the maximum background rate in
+ * each DETX-DETY plane of the background model.
+ ***************************************************************************/
+void GCTABackground3D::init_mc_max_rate(void) const
+{
+    // Initialise cache
+    m_mc_max.clear();
+
+    // Loop over all maps to determine the maximum map intensity in units
+    // of events/s/sr/MeV
+    for (int iebin = 0; iebin < m_num_energy; ++iebin) {
+
+        // Initialise maximum rate
+        double max_rate = 0.0;
+
+        // Loop over DETX-DETY plane
+        for (int ix = 0; ix < m_num_detx; ++ix) {
+            for (int iy = 0; iy < m_num_dety; ++iy) {
+
+                // Get bin index
+                int inx = index(ix, iy, iebin);
+
+                // Get background rate
+                double rate = m_background(m_inx_bgd, inx);
+                
+                // If background rate is larger than maximum then store the
+                // background rate as new maximum
+                if (rate > max_rate) {
+                    max_rate = rate;
+                }
+
+            } // endfor: looped over DETY
+        } // endfor: looped over DETX
+
+        // Append maximum rate
+        m_mc_max.push_back(max_rate);
+
+    } // endif: looped over all maps
     // Return
     return;
 }
+
+
+/***********************************************************************//**
+ * @brief Compute solid angle of pixel wedge
+ *
+ * @param[in] detx1 DETX of first edge.
+ * @param[in] dety1 DETY of first edge.
+ * @param[in] detx2 DETX of second edge.
+ * @param[in] dety2 DETY of second edge.
+ * @param[in] detx3 DETX of third edge.
+ * @param[in] dety3 DETY of third edge.
+ * @return Solid angle (steradians).
+ *
+ * Estimate the solid angle subtended by 3 coordinates using Huilier's
+ * theorem.
+ *
+ * Below, the definiton of the pixel cornes and sides are shown as used
+ * within the code.
+ *
+ *             a12
+ *         1---------2
+ *         |        / 
+ *         |       /  
+ *         |      /   
+ *         |     /    
+ *      a13|    /a23
+ *         |   / 
+ *         |  / 
+ *         | /
+ *         |/
+ *         3
+ *
+ ***************************************************************************/
+double GCTABackground3D::solid_angle(const double& detx1, const double& dety1,
+                                     const double& detx2, const double& dety2,
+                                     const double& detx3, const double& dety3) const
+{
+    // Initialise solid angle
+    double solidangle = 0.0;
+
+    // Set sky directions
+	double  theta1 = std::sqrt(detx1 * detx1 + dety1 * dety1);
+	double  phi1   = std::atan2(dety1, detx1);
+	double  theta2 = std::sqrt(detx2 * detx2 + dety2 * dety2);
+	double  phi2   = std::atan2(dety2, detx2);
+	double  theta3 = std::sqrt(detx3 * detx3 + dety3 * dety3);
+	double  phi3   = std::atan2(dety3, detx3);
+    GSkyDir dir1;
+    GSkyDir dir2;
+    GSkyDir dir3;
+    dir1.radec(phi1, gammalib::pihalf - theta1);
+    dir2.radec(phi2, gammalib::pihalf - theta2);
+    dir3.radec(phi3, gammalib::pihalf - theta3);
+
+    // Compute angular distances between pixel corners
+    double a12 = dir1.dist(dir2);
+    double a13 = dir1.dist(dir3);
+    double a23 = dir2.dist(dir3);
+
+    // Compute solid angle
+    double s   = 0.5 * (a12 + a23 + a13);
+    solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
+                                           std::tan(0.5*(s-a12)) *
+                                           std::tan(0.5*(s-a23)) *
+                                           std::tan(0.5*(s-a13))));
+
+    // Return solid angle
+    return solidangle;
+}
@@ -659,7 +659,7 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
         for (int i = 0; i < spectral.nodes(); ++i) {
             GEnergy energy    = spectral.energy(i);
             double  intensity = spectral.intensity(i);
-            double  norm      = m_spectral->eval(energy, events->tstart());
+            double  norm      = m_spectral->eval(energy);
             spectral.intensity(i, norm*intensity);
         }
@@ -634,7 +634,7 @@ GCTAEventAtom* GCTAResponseIrf::mc(const double&amp; area, const GPhoton&amp; photon,
     }
     // Continue only if event is detected
-    if (ran.uniform() <= ulimite) {
+    if ((ulimite > 0.0) && (ran.uniform() <= ulimite)) {
         // Apply deadtime correction
         double deadc = obs.deadc(photon.time());
@@ -2802,7 +2802,7 @@ GFitsImageDouble* GSkyMap::create_wcs_hdu(void) const
  * @return Solid angle (steradians).
  *
  * Estimate the solid angle subtended by 4 sky directions using Huilier's
- * theroem.
+ * theorem.
  *
  * Below, the definiton of the pixel cornes and sides are shown as used
  * within the code.
@@ -2891,7 +2891,7 @@ double GSkyMap::solidangle(const GSkyDir&amp; dir1, const GSkyDir&amp; dir2,
  * @return Solid angle (steradians).
  *
  * Estimate the solid angle subtended by 3 sky directions using Huilier's
- * theroem.
+ * theorem.
  *
  * Below, the definiton of the pixel cornes and sides are shown as used
  * within the code.