Improve elliptical Irf and Npred computation.

The binned analysis looks okay, but unbinned shows a strange behaviour for the semiminor axis parameter which results in a strong dependency of Npred on the semiminor axis value. This seems to prevent a nice convergence in unbinned mode. There is also an open issue to be solved related to the possibility that the minor axis is larger than the major axis. The actual code assumes that minor axis < major axis. If this is not fulfilled, the omega arc computations will be wrong.

Improve elliptical Irf and Npred computation.
The binned analysis looks okay, but unbinned shows a strange behaviour for the semiminor axis parameter which results in a strong dependency of Npred on the semiminor axis value. This seems to prevent a nice convergence in unbinned mode. There is also an open issue to be solved related to the possibility that the minor axis is larger than the major axis. The actual code assumes that minor axis < major axis. If this is not fulfilled, the omega arc computations will be wrong.
Jürgen Knödlseder
1 parent 1ee6318e
Showing 7 changed files with 172 additions and 57 deletions
include/GModelSpatialEllipticalDisk.hpp
inst/cta/src/GCTAResponseIrf.cpp
inst/cta/src/GCTAResponse_helpers.cpp
inst/cta/src/GCTAResponse_helpers.hpp
inst/cta/test/test_irf_elliptical.py
inst/cta/test/test_npred_elliptical.py
src/model/GModelSpatialEllipticalDisk.cpp
 /***************************************************************************
  *   GModelSpatialEllipticalDisk.hpp - Elliptical disk source model class  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013 by Michael Mayer                                    *
+ *  copyright (C) 2013-2014 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -1755,10 +1755,21 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp;       event,
     double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
  
     // Get semimajor and semiminor ellipse boundaries [radians]
-    double aspect_ratio = (model->semiminor() > 0) ?
-                          model->semimajor() / model->semiminor() : 0.0;
+    /*
+    double aspect_ratio = (model->semimajor() > 0) ?
+                           model->semiminor() / model->semimajor() : 0.0;
     double semimajor    = model->theta_max();
     double semiminor    = semimajor * aspect_ratio;
+    */
+    double semimajor    = model->semimajor() * gammalib::deg2rad;
+    double semiminor    = model->semiminor() * gammalib::deg2rad;
+
+    //TEST TEST TEST
+    if (semimajor < semiminor) {
+        std::cout << "Warning:";
+        std::cout << " semiminor=" << semiminor * gammalib::rad2deg;
+        std::cout << " semimajor=" << semimajor * gammalib::rad2deg;
+    }
  
     // Set zenith angle integration range for elliptical model
     double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
@@ -2326,10 +2337,21 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
     double rho_roi = roi.centre().dir().dist(centre);
  
     // Get semimajor and semiminor ellipse boundaries [radians]
-    double aspect_ratio = (model->semiminor() > 0) ?
-                          model->semimajor() / model->semiminor() : 0.0;
+    /*
+    double aspect_ratio = (model->semimajor() > 0) ?
+                           model->semiminor() / model->semimajor() : 0.0;
     double semimajor    = model->theta_max();
     double semiminor    = semimajor * aspect_ratio;
+    */
+    double semimajor    = model->semimajor() * gammalib::deg2rad;
+    double semiminor    = model->semiminor() * gammalib::deg2rad;
+
+    //TEST TEST TEST
+    if (semimajor < semiminor) {
+        std::cout << "Warning:";
+        std::cout << " semiminor=" << semiminor * gammalib::rad2deg;
+        std::cout << " semimajor=" << semimajor * gammalib::rad2deg;
+    }
  
     // Set offset angle integration range. We take here the ROI+PSF into
     // account to make no integrations beyond the point where the
@@ -2361,8 +2383,8 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
                                                 *model,
                                                 semimajor,
                                                 semiminor,
-                                                source.energy(),
-                                                source.time(),
+                                                srcEng,
+                                                srcTime,
                                                 cta,
                                                 rot,
                                                 rho_roi,
@@ -46,6 +46,7 @@
  
 /* __ Debug definitions __________________________________________________ */
 //#define G_DEBUG_INTEGRAL                             //!< Debug integration
+#define G_DEBUG_MODEL_ZERO           //!< Debug check for zero model values
  
 /* __ Constants __________________________________________________________ */
  
@@ -623,6 +624,14 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
             double cos_ph  = cos_rho * m_cos_rho_pnt;
             double sin_ph  = sin_rho * m_sin_rho_pnt;
  
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model (e.g. an elliptical
+            // disk model)
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
+
             // Setup integration kernel
             cta_irf_elliptical_kern_omega integrand(m_rsp,
                                                     m_model,
@@ -634,7 +643,7 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
                                                     m_obsEng,
                                                     m_posangle_obs,
                                                     m_omega_pnt,
-                                                    rho,
+                                                    rho_kluge,
                                                     cos_psf,
                                                     sin_psf,
                                                     cos_ph,
@@ -790,8 +799,32 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
     // Initialise IRF value
     double irf = 0.0;
  
+    // Compute azimuth angle in model coordinate system (radians)
+    double omega_model = omega + m_posangle_obs;
+
     // Evaluate sky model
-    double model = m_model.eval(m_rho, omega + m_posangle_obs, m_srcEng, m_srcTime);
+    double model = m_model.eval(m_rho, omega_model, m_srcEng, m_srcTime);
+
+    // Debug: test if model is non positive
+    #if defined(G_DEBUG_MODEL_ZERO)
+    if (model <= 0.0) {
+        double m_semiminor_rad = m_model.semiminor() * gammalib::deg2rad;
+        double m_semimajor_rad = m_model.semimajor() * gammalib::deg2rad;
+        double diff_angle      = omega_model - m_model.posangle() * gammalib::deg2rad;
+        double cosinus         = std::cos(diff_angle);
+        double sinus           = std::sin(diff_angle);
+        double arg1            = m_semiminor_rad * cosinus;
+        double arg2            = m_semimajor_rad * sinus;
+        double r_ellipse       = m_semiminor_rad * m_semimajor_rad /
+                                 std::sqrt(arg1*arg1 + arg2*arg2);
+        std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval";
+        std::cout << " zero model for (rho,omega)=";
+        std::cout << m_rho*gammalib::rad2deg << ",";
+        std::cout << omega*gammalib::rad2deg << ")";
+        std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
+        std::cout << std::endl;
+    }
+    #endif
  
     // Continue only if model is positive
     if (model > 0.0) {
@@ -851,10 +884,20 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
  *                     N_{\rm pred}(\rho,\omega) d\omega
  * \f]
  *
- * The azimuth angle integration range 
- * \f$[\omega_{\rm min}, \omega_{\rm max}\f$
- * is limited to an arc around the vector connecting the model centre to
- * the ROI centre. This limitation assures proper converges properly.
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the ROI centre, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
+ * for a given true photon energy and photon arrival time,
+ * \f$N_{\rm pred}(\rho,\omega)\f$ is the data space integral over the
+ * response function.
+ *
+ * The method performs the required coordinate transformations from the
+ * model system, spanned by \f$(\rho, \omega)\f$, to the system needed for
+ * Npred computations. Furthermore, the method limits the integration range
+ * to area where the ellipse intersects the ROI.
  ***************************************************************************/
 double cta_npred_elliptical_kern_rho::eval(const double& rho)
 {
@@ -875,14 +918,18 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
         // Continue only if arc length is positive
         if (domega > 0.0) {
  
-            // Compute phi integration range
-            double omega_min = m_posangle_roi - domega;
-            double omega_max = m_posangle_roi + domega;
-
-            // Compute sine and cosine of offset angle
+            // Compute sine and cosine terms for azimuthal integration
             double sin_rho = std::sin(rho);
             double cos_rho = std::cos(rho);
  
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model (e.g. an elliptical
+            // disk model)
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
+
             // Setup phi integration kernel
             cta_npred_elliptical_kern_omega integrand(m_rsp,
                                                       m_model,
@@ -890,8 +937,10 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
                                                       m_srcTime,
                                                       m_obs,
                                                       m_rot,
+                                                      rho_kluge,
                                                       sin_rho,
-                                                      cos_rho);
+                                                      cos_rho,
+                                                      m_posangle_roi);
  
             // Setup integrator
             GIntegral integral(&integrand);
@@ -928,11 +977,10 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
                 if (omega_width > 0.0) {
  
                     // Compute azimuth angle difference between ellipse
-                    // position angle and position angle of observed
-                    // photon in the model system. This angle will define
-                    // the reference point around which the circle arc. Make
-                    // sure that omega_0 is within [-pi,pi] thus that the omega
-                    // intervals are between [-2pi,2pi]
+                    // position angle and position angle of ROI in the model
+                    // system. This angle will define the reference point for
+                    // the circle arcs. Make sure that omega_0 is within [-pi,pi]
+                    // thus that the omega intervals are between [-2pi,2pi]
                     double omega_0 = m_model.posangle() * gammalib::deg2rad - 
                                      m_posangle_roi;
                     if (omega_0 > gammalib::pi) {
@@ -1005,34 +1053,68 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
  * Computes
  *
  * \f[
- *    S_{\rm p}(\rho,\omega | E, t) \, N_{\rm pred}(\rho,\omega)
+ *    S_{\rm p}(\omega | \rho, E, t) \, N_{\rm pred}(\omega | \rho)
  * \f]
+ *
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the centre of the Region of Interest (ROI), and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ *
+ * @todo Npred computation goes over sky coordinates. This can maybe be
+ *       optimized to reduce the number of coordinate transformations.
+ * @todo Check whether the Npred omega argument is the right one.
  ***************************************************************************/
 double cta_npred_elliptical_kern_omega::eval(const double& omega)
 {
     // Initialise Npred value
     double npred = 0.0;
  
-    // Compute sky direction vector in native coordinates
-    double  cos_omega = std::cos(omega);
-    double  sin_omega = std::sin(omega);
-    GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
-
-    // Rotate from native into celestial system
-    GVector cel = m_rot * native;
+    // Compute azimuth angle in model coordinate system (radians)
+    double omega_model = omega + m_posangle_roi;
  
-    // Set sky direction
-    GSkyDir srcDir;
-    srcDir.celvector(cel);
+    // Evaluate sky model
+    double model = m_model.eval(m_rho, omega_model, m_srcEng, m_srcTime);
+
+    // Debug: test if model is non positive
+    #if defined(G_DEBUG_MODEL_ZERO)
+    if (model <= 0.0) {
+        double m_semiminor_rad = m_model.semiminor() * gammalib::deg2rad;
+        double m_semimajor_rad = m_model.semimajor() * gammalib::deg2rad;
+        double diff_angle      = omega_model - m_model.posangle() * gammalib::deg2rad;
+        double cosinus         = std::cos(diff_angle);
+        double sinus           = std::sin(diff_angle);
+        double arg1            = m_semiminor_rad * cosinus;
+        double arg2            = m_semimajor_rad * sinus;
+        double r_ellipse       = m_semiminor_rad * m_semimajor_rad /
+                                 std::sqrt(arg1*arg1 + arg2*arg2);
+        std::cout << "*** WARNING: cta_npred_elliptical_kern_omega::eval";
+        std::cout << " zero model for (rho,omega)=";
+        std::cout << m_rho*gammalib::rad2deg << ",";
+        std::cout << omega*gammalib::rad2deg << ")";
+        std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
+        std::cout << std::endl;
+    }
+    #endif
+    
+    // Continue only if model is positive
+    if (model > 0.0) {
+    
+        // Compute sky direction vector in native coordinates
+        double  cos_omega = std::cos(omega_model);
+        double  sin_omega = std::sin(omega_model);
+        GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
  
-    // Set Photon
-    GPhoton photon(srcDir, m_srcEng, m_srcTime);
+        // Rotate from native into celestial system
+        GVector cel = m_rot * native;
  
-    // Get model value for this photon
-    double model = m_model.eval(photon);
+        // Set sky direction
+        GSkyDir srcDir;
+        srcDir.celvector(cel);
  
-    // Continue only if sky intensity is positive
-    if (model > 0.0) {
+        // Set Photon
+        GPhoton photon(srcDir, m_srcEng, m_srcTime);
  
         // Compute Npred for this sky direction
         npred = m_rsp.npred(photon, m_obs) * model;
@@ -687,26 +687,32 @@ public:
                                     const GTime&                   srcTime,
                                     const GCTAObservation&         obs,
                                     const GMatrix&                 rot,
+                                    const double&                  rho,
                                     const double&                  sin_rho,
-                                    const double&                  cos_rho) :
+                                    const double&                  cos_rho,
+                                    const double&                  posangle_roi) :
                                     m_rsp(rsp),
                                     m_model(model),
                                     m_srcEng(srcEng),
                                     m_srcTime(srcTime),
                                     m_obs(obs),
                                     m_rot(rot),
+                                    m_rho(rho),
                                     m_sin_rho(sin_rho),
-                                    m_cos_rho(cos_rho) { }
+                                    m_cos_rho(cos_rho),
+                                    m_posangle_roi(posangle_roi) { }
     double eval(const double& omega);
 protected:
-    const GCTAResponseIrf&         m_rsp;      //!< CTA response
-    const GModelSpatialElliptical& m_model;    //!< Model
-    const GEnergy&                 m_srcEng;   //!< True photon energy
-    const GTime&                   m_srcTime;  //!< True photon arrival time
-    const GCTAObservation&         m_obs;      //!< Pointer to observation
-    const GMatrix&                 m_rot;      //!< Rotation matrix
-    const double&                  m_sin_rho;  //!< Sine of offset angle
-    const double&                  m_cos_rho;  //!< Cosine of offset angle
+    const GCTAResponseIrf&         m_rsp;          //!< CTA response
+    const GModelSpatialElliptical& m_model;        //!< Model
+    const GEnergy&                 m_srcEng;       //!< True photon energy
+    const GTime&                   m_srcTime;      //!< True photon arrival time
+    const GCTAObservation&         m_obs;          //!< Pointer to observation
+    const GMatrix&                 m_rot;          //!< Rotation matrix
+    const double&                  m_rho;
+    const double&                  m_sin_rho;      //!< Sine of offset angle
+    const double&                  m_cos_rho;      //!< Cosine of offset angle
+    const double&                  m_posangle_roi; //!< Position angle of ROI
 };
  
  
@@ -413,8 +413,10 @@ if __name__ == &#39;__main__&#39;:
     # Set model parameters
     centre = gammalib.GSkyDir()
     centre.radec_deg(0.0, 0.0)
-    semimajor = 1.0
-    semiminor = 0.50
+    #semimajor = 1.0
+    #semiminor = 0.50
+    semimajor = 0.5
+    semiminor = 1.0
     posang    = 120.0
  
     # Set observed photon direction
@@ -180,9 +180,12 @@ def show_response(centre, semimajor, semiminor, posang, roiDir, pntDir, roi_radi
  
     # Rho limits
     rho_min = 0.0
-    rho_max = theta_max
+    #rho_max = theta_max
+    rho_max = rho_roi + radius_roi
     if (rho_roi > radius_roi):
         rho_min = rho_roi - radius_roi
+    if (rho_max > theta_max):
+        rho_max = theta_max
  
     # Setup ellipse as vector
     x = []
@@ -370,8 +373,8 @@ if __name__ == &#39;__main__&#39;:
         sys.exit()
  
     # Set Parameters
-    roi_radius = 0.5
-    delta_max  = 0.2
+    roi_radius = 1.0
+    delta_max  = 0.4
  
     # Set model parameters
     centre = gammalib.GSkyDir()
@@ -692,7 +692,7 @@ void GModelSpatialEllipticalDisk::free_members(void)
  ***************************************************************************/
 void GModelSpatialEllipticalDisk::update() const
 {
-    // Update if one axis has changed has changed
+    // Update if one axis has changed
     if (m_last_semiminor != semiminor() || m_last_semimajor != semimajor()) {
  
         // Store last values