GammaLib

It turned out that the GCTAPsf2D::mc method did not take into account the special case that the sigma and/or normalization of the second and third Gaussian components are zero. The following code solves this issue

double GCTAPsf2D::mc(GRan&         ran,
                     const double& logE,
                     const double& theta,
                     const double& phi,
                     const double& zenith,
                     const double& azimuth,
                     const bool&   etrue) const
{
    // Update the parameter cache
    update(logE, theta);

    // Select in which Gaussian we are
    double sigma = m_sigma1;
    double sum1  = m_sigma1;
    double sum2  = m_sigma2 * m_norm2;
    double sum3  = m_sigma3 * m_norm3;
    double sum   = sum1 + sum2 + sum3;
    double u     = ran.uniform() * sum;
    if (sum2 > 0.0 && u >= sum2) {
        sigma = m_sigma3;
    }
    else if (sum1 > 0.0 && u >= sum1) {
        sigma = m_sigma2;
    }

    // Now draw from the selected Gaussian
    double delta = sigma * ran.chisq2();

    // Return PSF offset
    return delta;
}