From 5fd8c80d3c9a78c9a0a365d0e95fd81c28fa3881 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 4 Feb 2014 11:22:32 +0100 Subject: [PATCH] Correct mc() method for the special case that the second and/or third Gaussian are zero. --- ChangeLog | 3 ++- NEWS | 6 +++++- inst/cta/src/GCTAPsf2D.cpp | 21 ++++++++++++++++++--- 3 files changed, 25 insertions(+), 5 deletions(-) diff --git a/ChangeLog b/ChangeLog index 4aea96e..4d8d21f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-01-31 Juergen Knoedlseder +2014-02-04 Juergen Knoedlseder * Bug fix Version 0.8.1 released ================================ @@ -11,6 +11,7 @@ Added vo module and GVOClient class to Python interface Track results of Python unit tests (#1119) Clean-up calibration database interface (#1102) + Correct 2D PSF Monte-Carlo simulations (#1128) 2014-01-18 Juergen Knoedlseder diff --git a/NEWS b/NEWS index 6694c2c..4908547 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-08-01 Author(s) : Juergen Knoedlseder -31 January 2014 +4 February 2014 1. Introduction @@ -130,6 +130,10 @@ the id attribute to the XML file in case that identifiers are specified. Furthermore, GCTAModelRadialAcceptance no longer write the temporal model component. +The GCTAPsf2D::mc() method did not sample the PSF correctly in case that +the sigma and normalization parameters were zero for the second and/or +third component. This has now been corrected for. + 21. COMPTEL interface --------------------- diff --git a/inst/cta/src/GCTAPsf2D.cpp b/inst/cta/src/GCTAPsf2D.cpp index 2f71b8d..3b080c9 100644 --- a/inst/cta/src/GCTAPsf2D.cpp +++ b/inst/cta/src/GCTAPsf2D.cpp @@ -317,10 +317,10 @@ double GCTAPsf2D::mc(GRan& ran, double sum3 = m_sigma3 * m_norm3; double sum = sum1 + sum2 + sum3; double u = ran.uniform() * sum; - if (u >= sum2) { + if (sum2 > 0.0 && u >= sum2) { sigma = m_sigma3; } - else if (u >= sum1) { + else if (sum1 > 0.0 && u >= sum1) { sigma = m_sigma2; } @@ -383,12 +383,27 @@ std::string GCTAPsf2D::print(const GChatter& chatter) const // Continue only if chatter is not silent if (chatter != SILENT) { + // Compute energy boundaries in TeV + double emin = m_psf.axis_lo(0,0); + double emax = m_psf.axis_hi(0,m_psf.axis(0)-1); + + // Compute offset angle boundaries in deg + double omin = m_psf.axis_lo(1,0); + double omax = m_psf.axis_hi(1,m_psf.axis(1)-1); + // Append header result.append("=== GCTAPsf2D ==="); // Append information result.append("\n"+gammalib::parformat("Filename")+m_filename); - result.append("\n"+m_psf.print(chatter)); + result.append("\n"+gammalib::parformat("Number of energy bins") + + gammalib::str(m_psf.axis(0))); + result.append("\n"+gammalib::parformat("Number of offset bins") + + gammalib::str(m_psf.axis(1))); + result.append("\n"+gammalib::parformat("Log10(Energy) range")); + result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV"); + result.append("\n"+gammalib::parformat("Offset angle range")); + result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg"); } // endif: chatter was not silent -- libgit2 0.21.4