Settle energy dispersion kludge (#2403)

Jürgen Knödlseder
1 parent 2b15414b
Showing 5 changed files with 259 additions and 458 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.5.1.rst
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/src/GCTAEdisp2D.cpp
-2018-02-23
+2018-03-19
  
         * Bug fix version 1.5.1 released
         ================================
  
-        Add kludge to handle noisy CTA 2D energy dispersion (#2342).
+        Add kludge to handle noisy CTA 2D energy dispersion (#2342, #2403).
         Prevent unintentional deallocation of arguments (#2340).
         Correct sky map smoothing.
  
 New Features and Important Changes in GammaLib 1.5.1
  
-23 February 2018
+19 Mars 2018
  
  
 1. Introduction
@@ -149,9 +149,7 @@ None
 20. CTA interface
 -----------------
 Add a kludge to handle the noisy energy dispersion that is shipped with the 1DC
-data and that is included in the Prod3 release. Noisy pixel values are skipped
-and a slight smoothing is applied in true energy to reduce the noise level
-(#2342).
+data and that is included in the Prod3 release (#2342, #2403).
  
 The GCTAObservation::response() and GCTAOnOffObservation::response() methods
 now only deallocate the previous content if the reference to the argument
@@ -10,6 +10,8 @@ GammaLib 1.5.1 is a bug fix release for GammaLib release 1.5.
  
 The following bugs were fixed:
  
+* [`2403 <https://cta-redmine.irap.omp.eu/issues/2403>`_] -
+  Refine kludge to handle noisy CTA 2D energy dispersion
 * [`2342 <https://cta-redmine.irap.omp.eu/issues/2342>`_] -
   Add kludge to handle noisy CTA 2D energy dispersion
 * [`2340 <https://cta-redmine.irap.omp.eu/issues/2340>`_] -
@@ -126,15 +126,23 @@ private:
     void normalize_table(void);
  
     // Kludge
-    void     clean_table(const int& threshold);
-    void     clip_table(const double& threshold);
-    void     denoise_table(const double& sigma);
-    GNdarray gaussian_array(const double& mean, const double& rms) const;
-    GNdarray get_gaussian_pars(const double& sigma) const;
-    GFft     fft_smooth_kernel(const int&    nbins,
-                               const double& sigma) const;
+    void     smooth_table(void);
+    GNdarray smooth_array(const GNdarray& array,
+                          const int&      nstart,
+                          const int&      nstop,
+                          const double&   limit) const;
+    double   smoothed_array_value(const int&      inx,
+                                  const GNdarray& array,
+                                  const int&      nodes,
+                                  const double&   limit) const;
+    void     get_moments(const int& itheta,
+                         GNdarray*  mean,
+                         GNdarray*  rms,
+                         GNdarray*  total) const;
     void     get_mean_rms(const GNdarray& array, double* mean, double* rms) const;
-    double   get_single_event_value(const GNdarray& array) const;
+    GNdarray gaussian_array(const double& mean,
+                            const double& rms,
+                            const double& total) const;
  
     // Protected classes
     class edisp_kern : public GFunction {
@@ -1200,9 +1200,7 @@ void GCTAEdisp2D::set_table(void)
  
     // Smooth energy dispersion table
     #if defined(G_SMOOTH_EDISP_KLUDGE)
-    denoise_table(10.0);
-    clean_table(15);
-    clip_table(0.001);
+    smooth_table();
     #endif
  
     // Normalize energy dispersion table
@@ -1362,166 +1360,80 @@ void GCTAEdisp2D::normalize_table(void)
 }
  
  
+
 /***********************************************************************//**
- * @brief Clean energy dispersion table
- *
- * @param[in] threshold Minimum number of consecutive non-zeros
- *
- * Removes noise from the energy dispersion table.
- *
- * For each offset angle and true energy, the method determines the maximum
- * value of the energy dispersion table and then walks towards smaller and
- * larger migration value to determine the first pixel for which the energy
- * dispersion becomes zero. All pixels beyond these pixels are then clipped
- * to zero. This removes noise that is at small or large migration values.
- *
- * The method then computes the number of remaining contiguos non-zero pixels
- * for each offset angle and true energy, and if this number is smaller than
- * a threshold, all pixels are set to zero. This removes true energy bins
- * that are sparesely filled.
+ * @brief Smoothed energy dispersion table
  ***************************************************************************/
-void GCTAEdisp2D::clean_table(const int& threshold)
+void GCTAEdisp2D::smooth_table(void)
 {
+    // Log entrance
+    //std::cout << "GCTAEdisp2D::smooth_table in" << std::endl;
+
     // Get axes dimensions
     int netrue = m_edisp.axis_bins(m_inx_etrue);
     int nmigra = m_edisp.axis_bins(m_inx_migra);
     int ntheta = m_edisp.axis_bins(m_inx_theta);
+    int npix   = netrue * nmigra;
  
     // Loop over all offset angles
     for (int itheta = 0; itheta < ntheta; ++itheta) {
  
-        // Loop ober all true energies
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
+        // Compute means, rms and total as function of etrue
+        GNdarray mean(netrue);
+        GNdarray rms(netrue);
+        GNdarray total(netrue);
+        get_moments(itheta, &mean, &rms, &total);
  
-            // Compute base index
-            int ibase = ietrue + itheta * netrue * nmigra;
-
-            // Determine index of maximum value
-            int    imax = 0;
-            double max  = m_edisp(m_inx_matrix,ibase+ietrue);
-            for (int imigra = 1; imigra < nmigra; ++imigra) {
-                int    inx   = ibase + imigra*netrue;
-                double value = m_edisp(m_inx_matrix, inx);
-                if (value > max) {
-                    imax = imigra;
-                    max  = value;
-                }
-            }
+        // Smooth all three
+        mean  = smooth_array(mean,  30, 30, 0.5);
+        rms   = smooth_array(rms,   30, 30, 0.03);
+        total = smooth_array(total, 30, 30, 0.0);
  
-            // Do nothing if all pixels are zero
-            if (max == 0.0) {
-                continue;
+        // Make sure that total is not negative
+        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
+            if (total(ietrue) < 0.0) {
+                total(ietrue) = 0.0;
             }
+        }
+        
+        // Replace matrix by Gaussians
+        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
  
-            // Walk from maximum to the first pixel
-            bool clip = false;
-            for (int imigra = imax; imigra >= 0; imigra--) {
-                int inx = ibase + imigra*netrue;
-                if (clip) {
-                    m_edisp(m_inx_matrix, inx) = 0.0;
-                }
-                else {
-                    double value = m_edisp(m_inx_matrix, inx);
-                    if (value == 0.0) {
-                        clip = true;
-                    }
-                }
-            }
+            // Continue only if there is information
+            if (total(ietrue) > 0.0) {
  
-            // Walk from maximum to the last pixel
-            clip = false;
-            for (int imigra = imax; imigra < nmigra; ++imigra) {
-                int inx = ibase + imigra*netrue;
-                if (clip) {
-                    m_edisp(m_inx_matrix, inx) = 0.0;
-                }
-                else {
-                    double value = m_edisp(m_inx_matrix, inx);
-                    if (value == 0.0) {
-                        clip = true;
-                    }
-                }
-            }
+                // Get Gaussian
+                GNdarray gauss = gaussian_array(mean(ietrue), rms(ietrue), total(ietrue));
  
-            // Determine number of contiguous non-zero pixels
-            int nonzero            = 0;
-            int contiguous_nonzero = 0;
-            for (int imigra = 0; imigra < nmigra; ++imigra) {
-                int    inx   = ibase + imigra*netrue;
-                double value = m_edisp(m_inx_matrix, inx);
-                if (value != 0.0) {
-                    nonzero++;
-                    if (nonzero > contiguous_nonzero) {
-                        contiguous_nonzero = nonzero;
+                // Compute base index
+                int ibase = itheta * npix + ietrue;
+                if ((mean(ietrue) > 0.0) && (rms(ietrue))) {
+                    for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
+                        m_edisp(m_inx_matrix, i) = gauss(imigra);
                     }
                 }
                 else {
-                    nonzero = 0;
-                }
-            }
-
-            // If there are less than threshold contiguous non-zero pixels we
-            // only have noise, hence all pixels are set to zero
-            if (contiguous_nonzero < threshold) {
-                for (int imigra = 0; imigra < nmigra; ++imigra) {
-                    int inx = ibase + imigra*netrue;
-                    m_edisp(m_inx_matrix, inx) = 0.0;
-                }
-            }
-
-        } // endfor: looped over all true energies
-
-    } // endfor: looped over all offset angles
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Clip energy dispersion table
- *
- * @param[in] threshold Clipping threshold as fraction of maximum value.
- *
- * Clips all migration vectors at a threshold that is a fraction of the
- * maximum value in the migration vector.
- ***************************************************************************/
-void GCTAEdisp2D::clip_table(const double& threshold)
-{
-    // Get axes dimensions
-    int netrue = m_edisp.axis_bins(m_inx_etrue);
-    int nmigra = m_edisp.axis_bins(m_inx_migra);
-    int ntheta = m_edisp.axis_bins(m_inx_theta);
-    int npix   = netrue * nmigra;
-
-    // Loop over all offset angles
-    for (int itheta = 0; itheta < ntheta; ++itheta) {
-
-        // Loop over all true energies
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
-
-            // Compute base index
-            int ibase = itheta * npix + ietrue;
-
-            // Compute maximum value
-            double max = 0.0;
-            for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
-                if (m_edisp(m_inx_matrix, i) > max) {
-                    max = m_edisp(m_inx_matrix, i);
+                    for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
+                        m_edisp(m_inx_matrix, i) = 0.0;
+                    }
                 }
-            }
+                
+            } // endif: there was information
  
-            // Clip values below threshold
-            double clip_value = max * threshold;
-            for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
-                if (m_edisp(m_inx_matrix, i) < clip_value) {
+            // ... otherwise reset the matrix to zero
+            else {
+                int ibase = itheta * npix + ietrue;
+                for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
                     m_edisp(m_inx_matrix, i) = 0.0;
                 }
             }
  
         } // endfor: looped over true energies
  
-    } // endfor: looped over all offset angles
+    } // endfor: looped over offset angles
+
+    // Log exit
+    //std::cout << "GCTAEdisp2D::smooth_table out" << std::endl;
  
     // Return
     return;
@@ -1529,334 +1441,198 @@ void GCTAEdisp2D::clip_table(const double&amp; threshold)
  
  
 /***********************************************************************//**
- * @brief Denoise energy dispersion table
- *
- * @param[in] sigma Gaussian smoothing single in true energy.
+ * @brief Smoothed array
  *
- * Denoise the energy dispersion table by replacing the energy dispersion
- * for each offset angle and true energy by a Gaussian function. The mean
- * and rms of the Gaussian is smoothed as function of true energy using a
- * Gaussian smoothing kernel with a width of @p sigma pixels.
+ * @param[in] array Array.
+ * @param[in] nstart Number of nodes used for regression at first array value.
+ * @param[in] nstop Number of nodes used for regression at last array value.
+ * @param[in] limit Limit for array element exclusion.
+ * @return Smoothed array.
  *
- * True energy values with no energy dispersion information are interpolated
- * from neighbouring energy dispersion values.
+ * Returns a smoothed array that is computed by locally adjusting a linear
+ * regression law to the array elements.
  ***************************************************************************/
-void GCTAEdisp2D::denoise_table(const double& sigma)
+GNdarray GCTAEdisp2D::smooth_array(const GNdarray& array,
+                                   const int&      nstart,
+                                   const int&      nstop,
+                                   const double&   limit) const
 {
-    // Log entrance
-    //std::cout << "GCTAEdisp2D::denoise_table in" << std::endl;
-
-    // Get axes dimensions
-    int netrue = m_edisp.axis_bins(m_inx_etrue);
-    int nmigra = m_edisp.axis_bins(m_inx_migra);
-    int ntheta = m_edisp.axis_bins(m_inx_theta);
-    int npix   = netrue * nmigra;
-
-    // Get Gaussian parameters for all offset angles and true energies
-    GNdarray gpars = get_gaussian_pars(sigma);
-
-    // Loop over all offset angles
-    for (int itheta = 0; itheta < ntheta; ++itheta) {
-
-        // Initialise sums
-        std::vector<double> sum_array;
-
-        // Loop over all true energies
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
-
-            // Compute base index
-            int ibase = itheta * npix + ietrue;
-
-            // Extract array
-            GNdarray array(nmigra);
-            for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
-                array(imigra) = m_edisp(m_inx_matrix, i);
-            }
-
-            // Compute single event value
-            double event_value = get_single_event_value(array);
-
-            // Smooth array
-            GNdarray smoothed_array(nmigra);
-            if (max(array) >= 2.0*event_value) {
-                smoothed_array = gaussian_array(gpars(itheta, ietrue, 0),
-                                                gpars(itheta, ietrue, 1));
-            }
-
-            // Compute total of smoothed array
-            double total = sum(smoothed_array);
-
-            // Store sum
-            sum_array.push_back(total);
-
-            // Put back smoothed array values in matrix
-            if (total > 0.0) {
-                for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
-                    m_edisp(m_inx_matrix, i) = smoothed_array(imigra) / total;
-                }
-            }
-
-        } // endfor: looped over all true energies
+    // Initialise smoothed array
+    GNdarray smoothed_array(array.size());
  
-        // Interpolate gaps
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
-
-            // Fall through if array is not empty
-            if (sum_array[ietrue] > 0.0) {
-                continue;
-            }
+    // Compute node step size
+    double nstep = double(nstop - nstart)/double(array.size());
  
-            // Find enclosing indices of non-empty arrays
-            int ileft  = ietrue;
-            int iright = ietrue;
-            for (; ileft >= 0; ileft--) {
-                if (sum_array[ileft] > 0.0) {
-                    break;
-                }
-            }
-            for (; iright < netrue; ++iright) {
-                if (sum_array[iright] > 0.0) {
-                    break;
-                }
-            }
-
-            // If encosing indices are valid then interpolate array
-            if ((ileft >= 0) && (iright < netrue)) {
-
-                // Compute weight factors
-                double norm   = 1.0 / double(iright-ileft);
-                double wleft  = double(iright - ietrue) * norm;
-                double wright = double(ietrue - ileft)  * norm;
-
-                // Fill array by linear interpolation
-                int idst       = itheta * npix + ietrue;
-                int isrc_left  = itheta * npix + ileft;
-                int isrc_right = itheta * npix + iright;
-                for (int imigra = 0; imigra < nmigra; ++imigra) {
-                    m_edisp(m_inx_matrix, idst) =
-                            wleft  * m_edisp(m_inx_matrix, isrc_left) +
-                            wright * m_edisp(m_inx_matrix, isrc_right);
-                    idst       += netrue;
-                    isrc_left  += netrue;
-                    isrc_right += netrue;
-                }
-
-            } // endif: there were valid enclosing indices
-        
-        } // endfor: looped over all true energies
-
-    } // endfor: looped over all offset angles
-
-    // Log exit
-    //std::cout << "GCTAEdisp2D::denoise_table out" << std::endl;
+    // Loop over all array elements and computed the smoothed value
+    for (int i = 0; i < array.size(); ++i) {
+        int nodes         = nstart + int(i*nstep);
+        smoothed_array(i) = smoothed_array_value(i, array, nodes, limit);
+    }
  
-    // Return
-    return;
+    // Return smoothed array
+    return smoothed_array;
 }
  
  
 /***********************************************************************//**
- * @brief Return Gaussian approximation of energy dispersion array
+ * @brief Get smoothed array value
  *
- * @param[in] mean Gaussian mean.
- * @param[in] rms Gaussian rms.
- * @return Gaussian approximation of energy dispersion array.
+ * @param[in] inx Index of array element.
+ * @param[in] array Array.
+ * @param[in] nodes Number of nodes used for regression.
+ * @param[in] limit Limit for array element exclusion.
+ * @return Smoothed array value.
  *
- * Returns a Gaussian approximation of the energy dispersion array by
- * computing the mean migration value and its root mean square and by using
- * these values as the centre and the width of a Gaussian function. The
- * Gaussian function is normalized so that the sum of the output array is
- * unity.
+ * Returns a smoothed array value that is derived from adjusting a linear
+ * law using regression to @p nodes adjacent array elements. Array elements
+ * with values below @p limit are excluded from the regression.
  ***************************************************************************/
-GNdarray GCTAEdisp2D::gaussian_array(const double& mean, const double& rms) const
+double GCTAEdisp2D::smoothed_array_value(const int&      inx,
+                                         const GNdarray& array,
+                                         const int&      nodes,
+                                         const double&   limit) const
 {
-    // Initialise empty Gaussian array
-    int      nmigra = m_edisp.axis_bins(m_inx_migra);
-    GNdarray gaussian(nmigra);
-
-    // If the array contains information then compute a Gaussian
-    // approximation
-    if (rms > 0.0) {
-
-        // Get reference to migration values
-        const GNodeArray& migras = m_edisp.axis_nodes(m_inx_migra);
-        int               nmigra = migras.size();
+    // Initialise variables
+    double mean_x = 0.0;
+    double mean_y = 0.0;
+    int    ileft  = inx - 1;
+    int    iright = inx + 1;
+    int    nleft  = 0;
+    int    nright = 0;
+
+    // Allocate vector of elements used for regression
+    std::vector<double> x_vals;
+    std::vector<double> y_vals;
+
+    // Add nodes on the left of the element of interest
+    while (nleft < nodes) {
+        if (ileft < 0) {
+            break;
+        }
+        if (array(ileft) > limit) {
+            double x = double(ileft);
+            double y = array(ileft);
+            x_vals.push_back(x);
+            y_vals.push_back(y);
+            mean_x += x;
+            mean_y += y;
+            nleft++;
+        }
+        ileft--;
+    }
  
-        // Compute Gaussian
-        double total_gauss = 0.0;
-        for (int imigra = 0; imigra < nmigra; ++imigra) {
-            double arg       = (migras[imigra] - mean) / rms;
-            double value     = std::exp(-0.5*arg*arg);
-            gaussian(imigra) = value;
-            total_gauss     += value;
+    // Add remaining nodes on the right of the element of interest
+    while (nright < 2*nodes-nleft) {
+        if (iright >= array.size()) {
+            break;
+        }
+        if (array(iright) > limit) {
+            double x = double(iright);
+            double y = array(iright);
+            x_vals.push_back(x);
+            y_vals.push_back(y);
+            mean_x += x;
+            mean_y += y;
+            nright++;
         }
+        iright++;
+    }
  
-        // Normalise Gaussian
-        if (total_gauss > 0.0) {
-            gaussian *= 1.0 / total_gauss;
+    // Add remaining nodes on the left if all right nodes were not exhausted
+    while (nleft < 2*nodes-nright) {
+        if (ileft < 0) {
+            break;
+        }
+        if (array(ileft) > limit) {
+            double x = double(ileft);
+            double y = array(ileft);
+            x_vals.push_back(x);
+            y_vals.push_back(y);
+            mean_x += x;
+            mean_y += y;
+            nleft++;
         }
+        ileft--;
+    }
  
-    } // endif: computed Gaussian approximation
+    // Compute mean x and y values
+    double total = double(nleft+nright);
+    if (total > 0.0) {
+        mean_x /= total;
+        mean_y /= total;
+    }
  
-    // Return Gaussian array
-    return gaussian;
+    // Compute regression line slope
+    double beta_nom   = 0.0;
+    double beta_denom = 0.0;
+    for (int i = 0; i < x_vals.size(); ++i) {
+        double x    = x_vals[i];
+        double y    = y_vals[i];
+        beta_nom   += (x - mean_x) * (y - mean_y);
+        beta_denom += (x - mean_x) * (x - mean_x);
+    }
+    double beta = beta_nom / beta_denom;
+
+    // Compute regression line offset
+    double alpha = mean_y - beta * mean_x;
+
+    // Compute value
+    double value = alpha + beta * double(inx);
+
+    // Return value
+    return value;
 }
  
  
 /***********************************************************************//**
- * @brief Compute mean and root mean square of migration array
+ * @brief Compute moments
  *
- * @param[in] sigma Smoothing parameter in true energy.
- * @return Array of mean and rms values for each offset angle and true energy.
+ * @param[in] itheta Offset angle index.
+ * @param[out] mean Pointer to mean array.
+ * @param[out] rms Pointer to rms array.
+ * @param[out] total Pointer to total array.
  *
- * Computes the mean and the root mean square of the migration array for each
- * offset angle and true energy and applies a Gaussian smoothing in true
- * energy.
+ * Computes the mean and root mean square values as function of true energy
+ * for a given offset angle.
  ***************************************************************************/
-GNdarray GCTAEdisp2D::get_gaussian_pars(const double& sigma) const
+void GCTAEdisp2D::get_moments(const int& itheta,
+                              GNdarray*  mean,
+                              GNdarray*  rms,
+                              GNdarray*  total) const
 {
     // Get axes dimensions
-    int netrue        = m_edisp.axis_bins(m_inx_etrue);
-    int nmigra        = m_edisp.axis_bins(m_inx_migra);
-    int ntheta        = m_edisp.axis_bins(m_inx_theta);
-    int npix          = netrue * nmigra;
+    int netrue = m_edisp.axis_bins(m_inx_etrue);
+    int nmigra = m_edisp.axis_bins(m_inx_migra);
+    int npix   = netrue * nmigra;
  
-    // Initialise result
-    GNdarray result(ntheta, netrue, 2);
+    // Loop over all true energies
+    for (int ietrue = 0; ietrue < netrue; ++ietrue) {
  
-    // Setup smoothing kernel
-    GFft fft_kernel = fft_smooth_kernel(netrue, sigma);
+        // Compute base index
+        int ibase = itheta * npix + ietrue;
  
-    // Loop over all offset angles
-    for (int itheta = 0; itheta < ntheta; ++itheta) {
+        // Extract array
+        GNdarray array(nmigra);
+        for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
+            array(imigra) = m_edisp(m_inx_matrix, i);
+        }
  
-        // Allocate working arrays
-        GNdarray work_mean(netrue);
-        GNdarray work_rms(netrue);
-        GNdarray work_sum(netrue);
+        // Store sum
+        (*total)(ietrue) = sum(array);
  
-        // Initialise mean and rms values
+        // Get Gaussian mean and rms
         double mean_value = 0.0;
         double rms_value  = 0.0;
+        get_mean_rms(array, &mean_value, &rms_value);
  
-        // Loop over all true energies
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
-
-            // Compute base index
-            int ibase = itheta * npix + ietrue;
-
-            // Extract array
-            GNdarray array(nmigra);
-            for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
-                array(imigra) = m_edisp(m_inx_matrix, i);
-            }
-
-            // Store sum (only results with positive sum will be kept)
-            work_sum(ietrue) = sum(array);
-
-            // Compute single event value
-            double event_value = get_single_event_value(array);
-
-            // Get Gaussian mean and rms. Use last mean and rms value in case
-            // that the array is empty
-            if (work_sum(ietrue) > 0) {
-                mean_value = 0.0;
-                rms_value  = 0.0;
-                get_mean_rms(array, &mean_value, &rms_value);
-            }
-
-            // Store mean and rms
-            work_mean(ietrue) = mean_value;
-            work_rms(ietrue)  = rms_value;
-
-        } // endfor: looped over all true energies
-
-        // Pad front
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
-            if (work_mean(ietrue) > 0.0) {
-                for (int i = 0; i < ietrue; ++i) {
-                    work_mean(i) = work_mean(ietrue);
-                    work_rms(i)  = work_rms(ietrue);
-                }
-                break;
-            }
-        }
-
-        // Pad back
-        for (int ietrue = netrue-1; ietrue >= 0; ietrue--) {
-            if (work_mean(ietrue) > 0.0) {
-                for (int i = ietrue; i < netrue; ++i) {
-                    work_mean(i) = work_mean(ietrue);
-                    work_rms(i)  = work_rms(ietrue);
-                }
-                break;
-            }
-        }
-
-        // Smooth mean
-        GFft fft_mean(work_mean);
-        GFft fft_mean_smooth = fft_mean * fft_kernel;
-        work_mean = fft_mean_smooth.backward();
+        // Store mean and rms
+        (*mean)(ietrue) = mean_value;
+        (*rms)(ietrue)  = rms_value;
  
-        // Smooth rms
-        GFft fft_rms(work_rms);
-        GFft fft_rms_smooth = fft_rms * fft_kernel;
-        work_rms = fft_rms_smooth.backward();
+    } // endfor: looped over all true energies
  
-        // Store result
-        for (int ietrue = 0; ietrue < netrue; ++ietrue) {
-            if (work_sum(ietrue) > 0.0) {
-                result(itheta, ietrue, 0) = work_mean(ietrue);
-                result(itheta, ietrue, 1) = work_rms(ietrue);
-            }
-        }
-
-    } // endif: looped over all offset angles
-
-    // Return result
-    return result;
-}
-
-
-/***********************************************************************//**
- * @brief Get FFT of smoothing kernel
- *
- * @param[in] nbins Number of kernel bins.
- * @param[in] sigma Gaussian kernel sigma in bins.
- * @return FFT of smoothing kernel
- *
- * Returns the fast-fourrier transform of a Gaussian smoothing kernel.
- ***************************************************************************/
-GFft GCTAEdisp2D::fft_smooth_kernel(const int&    nbins,
-                                    const double& sigma) const
-{
-    // Allocate kernel
-    GNdarray kernel(nbins);
-
-    // Initialise sum and compute Gaussian normalisation
-    double sum  =  0.0;
-    double norm = -0.5 / (sigma * sigma);
-
-    // Set Gaussian kernel
-    for (int i = 0; i < nbins; ++i) {
-        double value = std::exp(norm*double(i*i));
-        kernel(i) += value;
-        sum       += value;
-        if (i > 0) {
-            kernel(nbins-i) += value;
-            sum             += value;
-        }
-    }
-
-    // Normalize kernel
-    if (sum > 0.0) {
-        for (int i = 0; i < nbins; ++i) {
-            kernel(i) /= sum;
-        }
-    }
-
-    // Return FFT of kernel
-    return (GFft(kernel));
+    // Return
+    return;
 }
  
  
@@ -1910,35 +1686,52 @@ void GCTAEdisp2D::get_mean_rms(const GNdarray&amp; array,
  
  
 /***********************************************************************//**
- * @brief Estimate the value of a single event
+ * @brief Return Gaussian approximation of energy dispersion array
  *
- * @param[in] array Energy dispersion array.
- * @return Value of a single event.
+ * @param[in] mean Gaussian mean.
+ * @param[in] rms Gaussian rms.
+ * @param[in] total Gaussian total.
+ * @return Gaussian approximation of energy dispersion array.
  *
- * Returns the minimum non-zero value in an array as an estimate of the value
- * of a single event. If the array is empty the method returns zero.
+ * Returns a Gaussian approximation of the energy dispersion array by
+ * computing the mean migration value and its root mean square and by using
+ * these values as the centre and the width of a Gaussian function. The
+ * Gaussian function is normalized so that the sum of the output array is
+ * unity.
  ***************************************************************************/
-double GCTAEdisp2D::get_single_event_value(const GNdarray& array) const
+GNdarray GCTAEdisp2D::gaussian_array(const double& mean,
+                                     const double& rms,
+                                     const double& total) const
 {
-    // Initialise minimum value
-    double event_value = 1.0e30;
+    // Initialise empty Gaussian array
+    int      nmigra = m_edisp.axis_bins(m_inx_migra);
+    GNdarray gaussian(nmigra);
  
-    // Get minimum non-zero value
-    for (int i = 0; i < array.size(); ++i) {
-        double value = array(i);
-        if ((value > 0.0) && (value < event_value)) {
-            event_value = value;
+    // If the rms is valid then compute Gaussian
+    if (rms > 0.0) {
+
+        // Get reference to migration values
+        const GNodeArray& migras = m_edisp.axis_nodes(m_inx_migra);
+        int               nmigra = migras.size();
+
+        // Compute Gaussian
+        double total_gauss = 0.0;
+        for (int imigra = 0; imigra < nmigra; ++imigra) {
+            double arg       = (migras[imigra] - mean) / rms;
+            double value     = std::exp(-0.5*arg*arg);
+            gaussian(imigra) = value;
+            total_gauss     += value;
         }
-    }
  
-    // If no minimum value was encountered then set single event value to
-    // zero
-    if (event_value == 1.0e30) {
-        event_value = 0.0;
-    }
+        // Normalise Gaussian
+        if (total_gauss > 0.0) {
+            gaussian *= total / total_gauss;
+        }
  
-    // Return
-    return event_value;
+    } // endif: computed Gaussian
+
+    // Return Gaussian array
+    return gaussian;
 }