Add King profile PSF unit test, correctly normalize fixed delta_max PSF, and aut…

…omatically detected King profile PSF upon IRF loading. A unit test was added to verify that the King profile PSF is normalized correctly. For this purpose a King profile PSF calibration was inserted in slot 2 (bcf/000002) of the calibration database. As the delta_max has been fixed to 0.7 degrees, the King profile PSF now returns 0 for delta angles larger than 0.7 deg. Furthermore, the PSF is normalized so that the integral from 0 to 0.7 deg gives 1. In that way we assure that the PSF integral is always 1. The GCTAResponse::load_psf method now automatically detected the King profile PSF. Also, the load method now allow specification of a FITS file as the IRF name.

Add King profile PSF unit test, correctly normalize fixed delta_max PSF, and aut…
…omatically detected King profile PSF upon IRF loading. A unit test was added to verify that the King profile PSF is normalized correctly. For this purpose a King profile PSF calibration was inserted in slot 2 (bcf/000002) of the calibration database. As the delta_max has been fixed to 0.7 degrees, the King profile PSF now returns 0 for delta angles larger than 0.7 deg. Furthermore, the PSF is normalized so that the integral from 0 to 0.7 deg gives 1. In that way we assure that the PSF integral is always 1. The GCTAResponse::load_psf method now automatically detected the King profile PSF. Also, the load method now allow specification of a FITS file as the IRF name.
Jürgen Knödlseder
1 parent eda368fe
Showing 7 changed files with 269 additions and 164 deletions
inst/cta/caldb/data/cta/e/bcf/000002/irf_test.fits
inst/cta/include/GCTAPsfKing.hpp
inst/cta/pyext/GCTAPsfKing.i
inst/cta/src/GCTAPsfKing.cpp
inst/cta/src/GCTAResponse.cpp
inst/cta/test/test_CTA.cpp
inst/cta/test/test_CTA.hpp
 /***************************************************************************
- *       GCTAPsfKing.hpp - CTA point spread function vector class        *
+ *      GCTAPsfKing.hpp - King profile CTA point spread function class     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012 by Juergen Knoedlseder                              *
+ *  copyright (C) 2013 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file GCTAPsfKing.hpp
- * @brief CTA point spread function class definition with a King profile
+ * @brief King profile CTA point spread function class definition
  * @author Michael Mayer
  */
@@ -48,7 +48,7 @@ class GCTAPsfKing : public GCTAPsf {
 public:
     // Constructors and destructors
     GCTAPsfKing(void);
-    GCTAPsfKing(const std::string& filename);
+    explicit GCTAPsfKing(const std::string& filename);
     GCTAPsfKing(const GCTAPsfKing& psf);
     virtual ~GCTAPsfKing(void);
@@ -63,24 +63,24 @@ public:
                       const bool&   etrue = true) const;
     // Implemented pure virtual methods
-    void           clear(void);
+    void         clear(void);
     GCTAPsfKing* clone(void) const;
-    void           load(const std::string& filename);
-    std::string    filename(void) const;
-    double         mc(GRan&         ran,
-                      const double& logE, 
-                      const double& theta = 0.0, 
-                      const double& phi = 0.0,
-                      const double& zenith = 0.0,
-                      const double& azimuth = 0.0,
-                      const bool&   etrue = true) const;
-    double         delta_max(const double& logE, 
-                             const double& theta = 0.0, 
-                             const double& phi = 0.0,
-                             const double& zenith = 0.0,
-                             const double& azimuth = 0.0,
-                             const bool&   etrue = true) const;
-    std::string print(const GChatter& chatter = NORMAL) const;
+    void         load(const std::string& filename);
+    std::string  filename(void) const;
+    double       mc(GRan&         ran,
+                    const double& logE, 
+                    const double& theta = 0.0, 
+                    const double& phi = 0.0,
+                    const double& zenith = 0.0,
+                    const double& azimuth = 0.0,
+                    const bool&   etrue = true) const;
+    double       delta_max(const double& logE, 
+                           const double& theta = 0.0, 
+                           const double& phi = 0.0,
+                           const double& zenith = 0.0,
+                           const double& azimuth = 0.0,
+                           const bool&   etrue = true) const;
+    std::string  print(const GChatter& chatter = NORMAL) const;
 private:
@@ -91,15 +91,30 @@ private:
     void update(const double& logE, const double& theta) const;
     // Members
-   std::string       m_filename;   //!< Name of Aeff response file
-   GCTAResponseTable m_psf;        //!< PSF response table
-
-   mutable double m_par_logE; //!< Cache energy
-   mutable double m_par_theta; //!< Cache offset angle
-   mutable double m_par_norm; //!< King profile normalization
-   mutable double m_par_sigma; //!< King profile sigma (radians)
-   mutable double m_par_gamma; //!< King profile gamma parameter
+    std::string       m_filename;   //!< Name of Aeff response file
+    GCTAResponseTable m_psf;        //!< PSF response table
+    // Evaluation cache
+    mutable double m_par_logE;   //!< Cache energy
+    mutable double m_par_theta;  //!< Cache offset angle
+    mutable double m_par_norm;   //!< King profile normalization
+    mutable double m_par_sigma;  //!< King profile sigma (radians)
+    mutable double m_par_sigma2; //!< King profile sigma squared
+    mutable double m_par_gamma;  //!< King profile gamma parameter
 };
+
+/***********************************************************************//**
+ * @brief Return filename
+ *
+ * @return Filename from which point spread function was loaded.
+ *
+ * Returns filename from which point spread function was loaded.
+ ***************************************************************************/
+inline
+std::string GCTAPsfKing::filename(void) const
+{
+    return m_filename;
+}
+
 #endif /* GCTAPsfKing_HPP */
 /***************************************************************************
- *        GCTAPsfKing.i - CTA point spread function  class         *
+ *       GCTAPsfKing.i - King profile CTA point spread function class      *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012 by Michael Mayer                                    *
+ *  copyright (C) 2013 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file GCTAPsfKing.i
- * @brief CTA point spread function with a King Profile
+ * @brief King profile CTA point spread function class definition
  * @author Michael Mayer
  */
 %{
@@ -39,11 +39,10 @@
  * of energy as determined from a FITS table.
  ***************************************************************************/
 class GCTAPsfKing : public GCTAPsf {
-
 public:
     // Constructors and destructors
     GCTAPsfKing(void);
-    GCTAPsfKing(const std::string& filename);
+    explicit GCTAPsfKing(const std::string& filename);
     GCTAPsfKing(const GCTAPsfKing& psf);
     virtual ~GCTAPsfKing(void);
@@ -57,24 +56,23 @@ public:
                       const bool&   etrue = true) const;
     // Implemented pure virtual methods
-    void           clear(void);
+    void         clear(void);
     GCTAPsfKing* clone(void) const;
-    void           load(const std::string& filename);
-    std::string    filename(void) const;
-    double         mc(GRan&         ran,
-                      const double& logE, 
-                      const double& theta = 0.0, 
-                      const double& phi = 0.0,
-                      const double& zenith = 0.0,
-                      const double& azimuth = 0.0,
-                      const bool&   etrue = true) const;
-    double         delta_max(const double& logE, 
-                             const double& theta = 0.0, 
-                             const double& phi = 0.0,
-                             const double& zenith = 0.0,
-                             const double& azimuth = 0.0,
-                             const bool&   etrue = true) const;
-
+    void         load(const std::string& filename);
+    std::string  filename(void) const;
+    double       mc(GRan&         ran,
+                    const double& logE, 
+                    const double& theta = 0.0, 
+                    const double& phi = 0.0,
+                    const double& zenith = 0.0,
+                    const double& azimuth = 0.0,
+                    const bool&   etrue = true) const;
+    double       delta_max(const double& logE, 
+                           const double& theta = 0.0, 
+                           const double& phi = 0.0,
+                           const double& zenith = 0.0,
+                           const double& azimuth = 0.0,
+                           const bool&   etrue = true) const;
 };
@@ -82,4 +80,7 @@ public:
  * @brief GCTAPsfKing class extension
  ***************************************************************************/
 %extend GCTAPsfKing {
+    GCTAPsfKing copy() {
+        return (*self);
+    }
 };
 /***************************************************************************
- *       GCTAPsfKing.cpp - CTA point spread function vector class        *
+ *      GCTAPsfKing.cpp - King profile CTA point spread function class     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012 by Juergen Knoedlseder                              *
+ *  copyright (C) 2013 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -19,8 +19,8 @@
  *                                                                         *
  ***************************************************************************/
 /**
- * @file GCTAPsfKing.cpp
- * @brief CTA point spread function using a King profile
+ * @file GCTAPsfKing.hpp
+ * @brief King profile CTA point spread function class definition
  * @author Michael Mayer
  */
@@ -35,15 +35,18 @@
 #include "GCTAException.hpp"
 /* __ Method name definitions ____________________________________________ */
-#define G_LOAD                            "GCTAPsfKing::load(std::string&)"
+#define G_UPDATE                      "GCTAPsfKing::update(double&, double&)"
 /* __ Macros _____________________________________________________________ */
 /* __ Coding definitions _________________________________________________ */
+#define G_FIX_DELTA_MAX
 /* __ Debug definitions __________________________________________________ */
 /* __ Constants __________________________________________________________ */
+const double r_max = 0.7 * gammalib::deg2rad;  // Maximum delta for fixed 
+                                               // delta_max computation
 /*==========================================================================
@@ -169,16 +172,21 @@ GCTAPsfKing&amp; GCTAPsfKing::operator=(const GCTAPsfKing&amp; psf)
  * of sr^-1 for a given energy.
  ***************************************************************************/
 double GCTAPsfKing::operator()(const double& delta,
-                                 const double& logE, 
-                                 const double& theta, 
-                                 const double& phi,
-                                 const double& zenith,
-                                 const double& azimuth,
-                                 const bool&   etrue) const
+                               const double& logE, 
+                               const double& theta, 
+                               const double& phi,
+                               const double& zenith,
+                               const double& azimuth,
+                               const bool&   etrue) const
 {
 	// Initialise PSF value
 	double psf = 0.0;
+    // Compile option: set PSF to zero outside delta_max
+    #if defined(G_FIX_DELTA_MAX)
+    if (delta <= r_max) {
+    #endif
+
     // Update the parameter cache
     update(logE, theta);
@@ -186,10 +194,18 @@ double GCTAPsfKing::operator()(const double&amp; delta,
     if (m_par_norm > 0.0) {
 		// Compute PSF value
-		psf = m_par_norm*pow((1 + 1 / (2 * m_par_gamma) * pow(delta / m_par_sigma, 2)), -m_par_gamma);
+        double arg  = delta / m_par_sigma;
+        double arg2 = arg * arg;
+		psf = m_par_norm * 
+              std::pow((1.0 + 1.0 / (2.0 * m_par_gamma) * arg2), -m_par_gamma);
     } // endif: normalization was positive
+    // Compile option: set PSF to zero outside delta_max
+    #if defined(G_FIX_DELTA_MAX)
+    }
+    #endif
+
     // Return PSF
     return psf;
 }
@@ -249,10 +265,10 @@ void GCTAPsfKing::load(const std::string&amp; filename)
     GFits file(filename);
     // Get PSF table
-    const GFitsTable* table = file.table("POINT SPREAD FUNCTION");
+    const GFitsTable& table = *file.table("POINT SPREAD FUNCTION");
     // Read PSF table
-    m_psf.read(*table);
+    m_psf.read(table);
     // Set energy axis to logarithmic scale
     m_psf.axis_log10(0);
@@ -260,7 +276,7 @@ void GCTAPsfKing::load(const std::string&amp; filename)
     // Set offset angle axis to radians
     m_psf.axis_radians(1);
-    // Convert sigma parameters to radians
+    // Convert sigma parameter to radians
     m_psf.scale(1, gammalib::deg2rad);
     // Close PSF FITS file
@@ -275,18 +291,6 @@ void GCTAPsfKing::load(const std::string&amp; filename)
 /***********************************************************************//**
- * @brief Return filename
- *
- * @return Returns filename from which point spread function was loaded
- ***************************************************************************/
-std::string GCTAPsfKing::filename(void) const
-{
-    // Return filename
-    return m_filename;
-}
-
-
-/***********************************************************************//**
  * @brief Simulate PSF offset (radians)
  *
  * @param[in] ran Random number generator.
@@ -300,12 +304,12 @@ std::string GCTAPsfKing::filename(void) const
  * Draws a random offset for the King Profile
  ***************************************************************************/
 double GCTAPsfKing::mc(GRan&         ran,
-                         const double& logE, 
-                         const double& theta, 
-                         const double& phi,
-                         const double& zenith,
-                         const double& azimuth,
-                         const bool&   etrue) const
+                       const double& logE, 
+                       const double& theta, 
+                       const double& phi,
+                       const double& zenith,
+                       const double& azimuth,
+                       const bool&   etrue) const
 {
 	// Initialise random offset
 	double delta = 0.0;
@@ -317,7 +321,9 @@ double GCTAPsfKing::mc(GRan&amp;         ran,
     double u = ran.uniform();
     // Draw random offset using inversion sampling
-    delta = std::sqrt( ( pow(1.0 - u, 1.0 / (1.0 - m_par_gamma) ) - 1.0 ) * 2.0 * pow(m_par_sigma, 2.0) * m_par_gamma);
+    //delta = std::sqrt( ( pow(1.0 - u, 1.0 / (1.0 - m_par_gamma) ) - 1.0 ) * 2.0 * pow(m_par_sigma, 2.0) * m_par_gamma);
+    delta = std::sqrt( ( std::pow(1.0 - u, 1.0 / (1.0 - m_par_gamma) ) - 1.0 ) * 
+                      2.0 * m_par_sigma2 * m_par_gamma);
     // Return PSF offset
     return delta;
@@ -335,26 +341,30 @@ double GCTAPsfKing::mc(GRan&amp;         ran,
  * @param[in] etrue Use true energy (true/false). Not used.
  *
  * Determine the radius beyond which the PSF becomes negligible. This radius
- * is set by this method to where the containment fraction become 99.995% which equal
- * \f$5 \times \sigma\f$ of a Gaussian width.
+ * is set by this method to where the containment fraction become 99.995% 
+ * which equals \f$5 \times \sigma\f$ of a Gaussian width.
  ***************************************************************************/
 double GCTAPsfKing::delta_max(const double& logE,
-                                   const double& theta, 
-                                   const double& phi,
-                                   const double& zenith,
-                                   const double& azimuth,
-                                   const bool&   etrue) const
+                              const double& theta, 
+                              const double& phi,
+                              const double& zenith,
+                              const double& azimuth,
+                              const bool&   etrue) const
 {
+    // Compile option: use fixed maximum delta
+    #if defined(G_FIX_DELTA_MAX)
+    double radius = r_max;
+    #else
+
     // Update the parameter cache
     update(logE, theta);
     // Compute maximum PSF radius (99.995% containment)
-    double F = 0.99;
-
-    // todo: For now fix maximum radius to 0.7 degree until we found a better method.
-    // The 99% containment radius can become quite large due to the long tails of the
-    // function
-    double radius = 0.7*gammalib::deg2rad; //m_par_sigma * std::sqrt(2 * ( std::pow(1 - F, 1 / (1-m_par_gamma) ) -1 )* m_par_gamma);
+    double F      = 0.99995;
+    double u_max  = (std::pow((1.0 - F), (1.0/(1.0-m_par_gamma))) - 1.0) * 
+                    m_par_gamma;
+    double radius = m_par_sigma * std::sqrt(2.0 * u_max);
+    #endif
     // Return maximum PSF radius
     return radius;
@@ -398,11 +408,12 @@ void GCTAPsfKing::init_members(void)
     // Initialise members
     m_filename.clear();
     m_psf.clear();
-    m_par_logE = -1.0e30;
-    m_par_theta = -1.0;
-    m_par_norm = 0.0;
-    m_par_gamma = 0.0;
-    m_par_sigma = 0.0;
+    m_par_logE   = -1.0e30;
+    m_par_theta  = -1.0;
+    m_par_norm   = 0.0;
+    m_par_gamma  = 0.0;
+    m_par_sigma  = 0.0;
+    m_par_sigma2 = 0.0;
     // Return
     return;
@@ -417,12 +428,14 @@ void GCTAPsfKing::init_members(void)
 void GCTAPsfKing::copy_members(const GCTAPsfKing& psf)
 {
     // Copy members
-    m_filename  = psf.m_filename;
-    m_psf = psf.m_psf;
-    m_par_logE  = psf.m_par_logE;
-    m_par_sigma = psf.m_par_sigma;
-    m_par_norm = psf.m_par_norm;
-    m_par_gamma = psf.m_par_gamma;
+    m_filename   = psf.m_filename;
+    m_psf        = psf.m_psf;
+    m_par_logE   = psf.m_par_logE;
+    m_par_theta  = psf.m_par_theta;
+    m_par_norm   = psf.m_par_norm;
+    m_par_gamma  = psf.m_par_gamma;
+    m_par_sigma  = psf.m_par_sigma;
+    m_par_sigma2 = psf.m_par_sigma2;
     // Return
     return;
@@ -434,9 +447,6 @@ void GCTAPsfKing::copy_members(const GCTAPsfKing&amp; psf)
  ***************************************************************************/
 void GCTAPsfKing::free_members(void)
 {
-	m_filename.clear();
-	m_psf.clear();
-
     // Return
     return;
 }
@@ -447,6 +457,9 @@ void GCTAPsfKing::free_members(void)
  *
  * @param[in] logE Log10 of the true photon energy (TeV).
  *
+ * @exception GException::invalid_value
+ *            No valid point spread function information has been found.
+ *
  * This method updates the PSF parameter cache. As the performance table PSF
  * only depends on energy, the only parameter on which the cache values
  * depend is the energy.
@@ -457,20 +470,42 @@ void GCTAPsfKing::update(const double&amp; logE, const double&amp; theta) const
     if (logE != m_par_logE || theta != m_par_theta) {
         // Save parameters
-        m_par_logE = logE;
+        m_par_logE  = logE;
         m_par_theta = theta;
         // Determine sigma and gamma by interpolating between nodes
         std::vector<double> pars = m_psf(logE,theta);
-        m_par_gamma = pars[0];
-        m_par_sigma = pars[1];
+        // Throw an exception if there are not 2 parameters
+        if (pars.size() != 2) {
+            std::string msg = gammalib::str(pars.size()) + " parameters have"
+                              " been found in the response table of the"
+                              " King profile response function while 2"
+                              " parameters are expected.\n"
+                              "Possibly, the point spread function information"
+                              " has not yet been loaded. Please load the point"
+                              " spread function before using it.";
+            throw GException::invalid_value(G_UPDATE, msg);
+        }
+
+        // Set parameters
+        m_par_gamma  = pars[0];
+        m_par_sigma  = pars[1];
+        m_par_sigma2 = m_par_sigma * m_par_sigma;
         // Determine normalisation for given parameters
-        m_par_norm = 1.0 / gammalib::twopi * pow(m_par_sigma, -2)*(1.0 - 1.0 / m_par_gamma);
+        m_par_norm = 1.0 / gammalib::twopi * (1.0 - 1.0 / m_par_gamma) /
+                     m_par_sigma2;
+
+        // Optionally correct for fixed delta_max
+        #if defined(G_FIX_DELTA_MAX)
+        double u_max = (r_max*r_max) / (2.0 * m_par_sigma2);
+        double norm  = 1.0 - std::pow((1.0 + u_max/m_par_gamma), 1.0-m_par_gamma);
+        m_par_norm /= norm;
+        #endif
+
     }
     // Return
     return;
 }
-
@@ -53,6 +53,7 @@
 #include "GCTAPsf2D.hpp"
 #include "GCTAPsfVector.hpp"
 #include "GCTAPsfPerfTable.hpp"
+#include "GCTAPsfKing.hpp"
 #include "GCTAAeff.hpp"
 #include "GCTAPsf.hpp"
 #include "GCTAEdisp.hpp"
@@ -549,8 +550,13 @@ void GCTAResponse::load(const std::string&amp; irfname)
     // Restore calibration database name
     m_caldb = caldb;
-    // Build filename
-    std::string filename = m_caldb + "/" + irfname + ".dat";
+    // Build filename. If file is not a FITS file then add ".dat" extension
+    // which is typical for a performance table.
+    std::string filename = m_caldb + "/" + irfname;
+    if (irfname.find_last_of(".fits")    != irfname.length()-1 ||
+        irfname.find_last_of(".fits.gz") != irfname.length()-1) {
+        filename += ".dat";
+    }
     // Load effective area
     load_aeff(filename);
@@ -679,10 +685,22 @@ void GCTAResponse::load_psf(const std::string&amp; filename)
         GFits file(filename);
         // If file contains a "POINT SPREAD FUNCTION" extension then load it
-        // as CTA response table
+        // as either a King profile PSF or a 2D PSF
         if (file.contains("POINT SPREAD FUNCTION")) {
-            file.close();
-            m_psf = new GCTAPsf2D(filename);
+            const GFitsTable& table = *file.table("POINT SPREAD FUNCTION");
+            if (table.contains("GAMMA") && table.contains("SIGMA")) {
+                file.close();
+                m_psf = new GCTAPsfKing(filename);
+            }
+            else if (table.contains("SCALE") && table.contains("SIGMA_1") &&
+                     table.contains("AMPL_2") && table.contains("SIGMA_2") &&
+                     table.contains("AMPL_3") && table.contains("SIGMA_3")) {
+                file.close();
+                m_psf = new GCTAPsf2D(filename);
+            }
+            else {
+                file.close();
+            }
         }
         // ... else load it has PSF vector 
@@ -40,16 +40,18 @@
 /* __ Globals ____________________________________________________________ */
 /* __ Constants __________________________________________________________ */
-const std::string datadir       = PACKAGE_SOURCE"/inst/cta/test/data";
-const std::string cta_caldb     = PACKAGE_SOURCE"/inst/cta/caldb";
-const std::string cta_irf       = "kb_E_50h_v3";
-const std::string cta_events    = datadir+"/crab_events.fits.gz";
-const std::string cta_cntmap    = datadir+"/crab_cntmap.fits.gz";
-const std::string cta_bin_xml   = datadir+"/obs_binned.xml";
-const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml";
-const std::string cta_model_xml = datadir+"/crab.xml";
-const std::string cta_rsp_xml   = datadir+"/rsp_models.xml";
+const std::string datadir        = PACKAGE_SOURCE"/inst/cta/test/data";
+const std::string cta_caldb      = PACKAGE_SOURCE"/inst/cta/caldb";
+const std::string cta_irf        = "kb_E_50h_v3";
+const std::string cta_events     = datadir+"/crab_events.fits.gz";
+const std::string cta_cntmap     = datadir+"/crab_cntmap.fits.gz";
+const std::string cta_bin_xml    = datadir+"/obs_binned.xml";
+const std::string cta_unbin_xml  = datadir+"/obs_unbinned.xml";
+const std::string cta_model_xml  = datadir+"/crab.xml";
+const std::string cta_rsp_xml    = datadir+"/rsp_models.xml";
 const std::string cta_modbck_xml = datadir+"/cta_modelbg.xml";
+const std::string cta_caldb_king = PACKAGE_SOURCE"/inst/cta/caldb/data/cta/e/bcf/000002";
+const std::string cta_irf_king   = "irf_test.fits";
 /***********************************************************************//**
@@ -64,6 +66,7 @@ void TestGCTAResponse::set(void)
     append(static_cast<pfunction>(&TestGCTAResponse::test_response), "Test response");
     append(static_cast<pfunction>(&TestGCTAResponse::test_response_aeff), "Test effective area");
     append(static_cast<pfunction>(&TestGCTAResponse::test_response_psf), "Test PSF");
+    append(static_cast<pfunction>(&TestGCTAResponse::test_response_psf_king), "Test King profile PSF");
     append(static_cast<pfunction>(&TestGCTAResponse::test_response_npsf), "Test integrated PSF");
     append(static_cast<pfunction>(&TestGCTAResponse::test_response_irf_diffuse), "Test diffuse IRF");
     append(static_cast<pfunction>(&TestGCTAResponse::test_response_npred_diffuse), "Test diffuse IRF integration");
@@ -236,6 +239,44 @@ void TestGCTAResponse::test_response_psf(void)
 /***********************************************************************//**
+ * @brief Test CTA King profile psf computation
+ *
+ * The Psf computation is tested by integrating numerically the Psf
+ * function. Integration is done in a rather simplistic way, by stepping
+ * radially away from the centre. The integration is done for a set of
+ * energies from 0.1-10 TeV.
+ ***************************************************************************/
+void TestGCTAResponse::test_response_psf_king(void)
+{
+    // Load response
+    GCTAResponse rsp;
+    rsp.caldb(cta_caldb_king);
+    rsp.load(cta_irf_king);
+
+    // Integrate Psf
+    GEnergy eng;
+    for (double e = 0.1; e < 10.0; e *= 2.0) {
+        eng.TeV(e);
+        double r_max = rsp.psf()->delta_max(eng.log10TeV()) * gammalib::rad2deg;
+        double r     = 0.0;
+        double dr    = 0.0001;
+        int    steps = int(1.0/dr);
+        double sum   = 0.0;
+        for (int i = 0; i < steps; ++i) {
+            r   += dr;
+            sum += rsp.psf(r * gammalib::deg2rad, 0.0, 0.0, 0.0, 0.0, eng.log10TeV()) *
+                   gammalib::twopi * std::sin(r * gammalib::deg2rad) * dr *
+                   gammalib::deg2rad;
+        }
+        test_value(sum, 1.0, 0.001, "PSF integration for "+eng.print());
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Test CTA npsf computation
  ***************************************************************************/
 void TestGCTAResponse::test_response_npsf(void)
@@ -427,6 +468,30 @@ void TestGCTAResponse::test_response_npred_diffuse(void)
     return;
 }
+
+/***********************************************************************//**
+ * @brief Test CTA response handling
+ ***************************************************************************/
+void TestGCTAResponse::test_response(void)
+{
+    // Test CTA response loading
+    test_try("Test CTA response loading");
+    try {
+        // Load response
+        GCTAResponse rsp;
+        rsp.caldb(cta_caldb);
+        rsp.load(cta_irf);
+        test_try_success();
+    }
+    catch (std::exception &e) {
+        test_try_failure(e);
+    }
+
+    // Return
+    return;
+}
+
+
 /***********************************************************************//**
  * @brief Test CTA Model Background Npred computation
  *
@@ -511,29 +576,6 @@ void TestGCTAModelBackground::test_modelbg_npred(void)
 /***********************************************************************//**
- * @brief Test CTA response handling
- ***************************************************************************/
-void TestGCTAResponse::test_response(void)
-{
-    // Test CTA response loading
-    test_try("Test CTA response loading");
-    try {
-        // Load response
-        GCTAResponse rsp;
-        rsp.caldb(cta_caldb);
-        rsp.load(cta_irf);
-        test_try_success();
-    }
-    catch (std::exception &e) {
-        test_try_failure(e);
-    }
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
  * @brief Test unbinned observation handling
  ***************************************************************************/
 void TestGCTAObservation::test_unbinned_obs(void)
@@ -49,6 +49,7 @@ public:
     virtual TestGCTAResponse* clone(void) const;
     void                      test_response_aeff(void);
     void                      test_response_psf(void);
+    void                      test_response_psf_king(void);
     void                      test_response_npsf(void);
     void                      test_response_irf_diffuse(void);
     void                      test_response_npred_diffuse(void);