diff --git a/ChangeLog b/ChangeLog index 6bbb12c..8189ab4 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,9 +1,9 @@ -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. diff --git a/NEWS b/NEWS index 1d21f40..801774f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ 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 diff --git a/doc/source/admin/release_history/1.5.1.rst b/doc/source/admin/release_history/1.5.1.rst index 33b696f..7456838 100644 --- a/doc/source/admin/release_history/1.5.1.rst +++ b/doc/source/admin/release_history/1.5.1.rst @@ -10,6 +10,8 @@ GammaLib 1.5.1 is a bug fix release for GammaLib release 1.5. The following bugs were fixed: +* [`2403 `_] - + Refine kludge to handle noisy CTA 2D energy dispersion * [`2342 `_] - Add kludge to handle noisy CTA 2D energy dispersion * [`2340 `_] - diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp index f3bf420..3f57195 100644 --- a/inst/cta/include/GCTAEdisp2D.hpp +++ b/inst/cta/include/GCTAEdisp2D.hpp @@ -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 { diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp index 7ce4cb4..6dd72bb 100644 --- a/inst/cta/src/GCTAEdisp2D.cpp +++ b/inst/cta/src/GCTAEdisp2D.cpp @@ -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& 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 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 x_vals; + std::vector 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& 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; }