Implement response computation for extended sources

Jürgen Knödlseder
1 parent 75b5e302
Showing 2 changed files with 311 additions and 21 deletions
inst/com/include/GCOMResponse.hpp
inst/com/src/GCOMResponse.cpp
@@ -111,6 +111,13 @@ private:
                                 const GObservation& obs,
                                 GMatrix*            gradients = NULL) const;
  
+    // Support methods
+    GVector irf_extended(const GModelSky&    model,
+                         const GObservation& obs,
+                         const GSkyDir&      model_dir,
+                         const double&       theta_max,
+                         GMatrix*            gradients = NULL) const;
+
     // Private data members
     GCaldb              m_caldb;             //!< Calibration database
     std::string         m_rspname;           //!< Response name
@@ -32,6 +32,8 @@
 #include <typeinfo>
 #include "GTools.hpp"
 #include "GMath.hpp"
+#include "GVector.hpp"
+#include "GMatrix.hpp"
 #include "GFits.hpp"
 #include "GCaldb.hpp"
 #include "GEvent.hpp"
@@ -44,6 +46,8 @@
 #include "GFitsImageFloat.hpp"
 #include "GModelSky.hpp"
 #include "GModelSpatialPointSource.hpp"
+#include "GModelSpatialRadial.hpp"
+#include "GModelSpatialElliptical.hpp"
 #include "GModelSpectralConst.hpp"
 #include "GCOMResponse.hpp"
 #include "GCOMObservation.hpp"
@@ -66,6 +70,8 @@
                                                   " GObservation&, GMatrix*)"
 #define G_IRF_DIFFUSE "GCOMResponse::irf_diffuse(GModelSky&, GObservation&, "\
                                                                   "GMatrix*)"
+#define G_IRF_EXTENDED              "GCOMResponse::irf_extended(GModelSky&, "\
+                                "GObservation&, GSkyDir&, double&, GMatrix*)"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -972,25 +978,25 @@ GVector GCOMResponse::irf_ptsrc(const GModelSky&amp;    model,
  * @return Instrument response to radial source  for all events in
  *         observation (\f$cm^2\f$).
  *
- * @todo Implement method.
+ * Returns the instrument response to a radial source for all events in
+ * the observations. See GCOMResponse::irf_extended for more information.
  ***************************************************************************/
 GVector GCOMResponse::irf_radial(const GModelSky&    model,
                                  const GObservation& obs,
                                  GMatrix*            gradients) const
 {
-    // Throw exception
-    std::string msg = "Response computation not yet implemented "
-                      "for spatial model type \""+
-                      model.spatial()->type()+"\".";
-    throw GException::feature_not_implemented(G_IRF_RADIAL, msg);
+    // Get pointer to radial model
+    const GModelSpatialRadial* radial =
+          static_cast<const GModelSpatialRadial*>(model.spatial());
  
-    // Get number of events
-    int nevents = obs.events()->size();
+    // Get references to model direction and maximum model radius
+    const GSkyDir& model_dir = radial->dir();
+    const double&  theta_max = radial->theta_max();
  
-    // Initialise result
-    GVector irfs(nevents);
+    // Call support method
+    GVector irfs = irf_extended(model, obs, model_dir, theta_max, gradients);
  
-    // Return IRF value
+    // Return IRF vector
     return irfs;
 }
  
@@ -1004,23 +1010,23 @@ GVector GCOMResponse::irf_radial(const GModelSky&amp;    model,
  * @return Instrument response to elliptical source for all events in
  *         observation (\f$cm^2\f$).
  *
- * @todo Implement method.
+ * Returns the instrument response to an elliptical source for all events in
+ * the observations. See GCOMResponse::irf_extended for more information.
  ***************************************************************************/
 GVector GCOMResponse::irf_elliptical(const GModelSky&    model,
                                      const GObservation& obs,
                                      GMatrix*            gradients) const
 {
-    // Throw exception
-    std::string msg = "Response computation not yet implemented "
-                      "for spatial model type \""+
-                      model.spatial()->type()+"\".";
-    throw GException::feature_not_implemented(G_IRF_ELLIPTICAL, msg);
+    // Get pointer to elliptical model
+    const GModelSpatialElliptical* elliptical =
+          static_cast<const GModelSpatialElliptical*>(model.spatial());
  
-    // Get number of events
-    int nevents = obs.events()->size();
+    // Get references to model direction and maximum model radius
+    const GSkyDir& model_dir = elliptical->dir();
+    const double&  theta_max = elliptical->theta_max();
  
-    // Initialise result
-    GVector irfs(nevents);
+    // Call support method
+    GVector irfs = irf_extended(model, obs, model_dir, theta_max, gradients);
  
     // Return IRF value
     return irfs;
@@ -1209,3 +1215,280 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp;    model,
     // Return IRF vector
     return irfs;
 }
+
+
+/***********************************************************************//**
+ * @brief Return instrument response to extended source
+ *
+ * @param[in] model Sky model.
+ * @param[in] obs Observation.
+ * @param[in] model_dir Centre direction of extended model.
+ * @param[in] theta_max Maximum theta angle (radians).
+ * @param[out] gradients Gradients matrix.
+ * @return Instrument response to radial source  for all events in
+ *         observation (\f$cm^2\f$).
+ *
+ * Returns the instrument response to an extended source for all events in
+ * the observations.
+ *
+ * The computation is done by integrating the extended model for each pixel
+ * in Chi and Psi over a circular region centred on the Chi/Psi pixel with
+ * a radius equal to the maximum Phigeo value. The integration is done in
+ * a spherical system, spanned by Phigeo and azimuth angle. Only coordiates
+ * are considered that plausibly overlap with the extended model. The
+ * integration is done by numerical summation. The integration step is by
+ * default set to 1/10 of the model size, with the constraint that it never
+ * becomes larger than the IAQ Phigeo binsize. This should allow to capture
+ * all extended models with sufficient precision in an acceptable amount of
+ * computing time.
+ *
+ * @p gradients is an optional matrix where the number of rows corresponds
+ * to the number of events in the observation and the number of columns
+ * corresponds to the number of spatial model parameters. Since no gradients
+ * are computed in this method the content of @p gradients if not altered.
+ ***************************************************************************/
+GVector GCOMResponse::irf_extended(const GModelSky&    model,
+                                   const GObservation& obs,
+                                   const GSkyDir&      model_dir,
+                                   const double&       theta_max,
+                                   GMatrix*            gradients) const
+{
+    // Extract COMPTEL observation
+    const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
+    if (obs_ptr == NULL) {
+        std::string cls = std::string(typeid(&obs).name());
+        std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
+                          "observations. Please specify a COMPTEL observation "
+                          "as argument.";
+        throw GException::invalid_argument(G_IRF_EXTENDED, msg);
+    }
+
+    // Extract COMPTEL event cube
+    const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(obs_ptr->events());
+    if (cube == NULL) {
+        std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
+                          "not contain a COMPTEL event cube. Please specify "
+                          "a COMPTEL observation containing and event cube.";
+        throw GException::invalid_argument(G_IRF_EXTENDED, msg);
+    }
+
+    // Throw an exception if COMPTEL response is not set or if
+    if (m_iaq.empty()) {
+        std::string msg = "COMPTEL response is empty. Please initialise the "
+                          "response with an \"IAQ\".";
+        throw GException::invalid_value(G_IRF_EXTENDED, msg);
+    }
+    else if (m_phigeo_bin_size == 0.0) {
+        std::string msg = "COMPTEL response has a zero Phigeo bin size. "
+                          "Please initialise the response with a valid "
+                          "\"IAQ\".";
+        throw GException::invalid_value(G_IRF_EXTENDED, msg);
+    }
+
+    // Get number of Chi/Psi pixels, Phibar layers and event bins
+    int npix    = cube->naxis(0) * cube->naxis(1);
+    int nphibar = cube->naxis(2);
+    int nevents = cube->size();
+
+    // Initialise result
+    GVector irfs(nevents);
+
+    // Initialise some variables
+    double         phigeo_min    = m_phigeo_min      * gammalib::deg2rad;
+    double         phigeo_bin    = m_phigeo_bin_size * gammalib::deg2rad;
+    const GSkyMap& drx           = obs_ptr->drx().map();
+    const double*  drg           = obs_ptr->drg().map().pixels();
+    double         cos_theta_max = std::cos(theta_max);
+
+    // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
+    double iaq_norm = obs_ptr->deadc() /
+                      (obs_ptr->ontime() * cube->dre().tof_correction());
+
+    // Loop over Chi and Psi
+    for (int ipix = 0; ipix < npix; ++ipix) {
+
+        // Get pointer to event bin
+        const GCOMEventBin* bin = (*cube)[ipix];
+
+        // Get reference to instrument direction
+        const GCOMInstDir& obsDir = bin->dir();
+
+        // Get reference to Phigeo sky direction
+        const GSkyDir& phigeoDir = obsDir.dir();
+
+        // Compute angle between model centre and Phigeo sky direction (radians)
+        // and position angle of model with respect to Phigeo sky direction (radians)
+        //double zeta = phigeoDir.dist(radial->dir());
+        double zeta = model_dir.dist(phigeoDir);
+        double phi0 = phigeoDir.posang(model_dir);
+
+        // Setup rotation matrix for faster conversion towards sky coordinates
+        GMatrix ry;
+        GMatrix rz;
+        ry.eulery(phigeoDir.dec_deg() - 90.0);
+        rz.eulerz(-phigeoDir.ra_deg());
+        GMatrix rot = (ry * rz).transpose();
+
+        // Setup azimuth angle range computation
+        double sin_zeta = std::sin(zeta);
+        double cos_zeta = std::cos(zeta);
+
+        // Setup Phigeo integration range
+        double phigeo_low  = (zeta > theta_max) ? zeta - theta_max : 0.0;
+        double phigeo_up   = zeta + theta_max;
+
+        // Setup Phigeo integration step-size and number of integration steps. By
+        // default the integration step size is one tenth of the model diameter, but
+        // not larger than the Phigeo binning of the IAQ. We also make sure that there
+        // are at least two Phigeo integration steps.
+        double phigeo_length = phigeo_up - phigeo_low;
+        double phigeo_step   = 0.1 * phigeo_length;
+        if (phigeo_step > phigeo_bin) {
+            phigeo_step = phigeo_bin;
+        }
+        int nphigeo = int(phigeo_length/phigeo_step + 0.5);
+        if (nphigeo < 2) {
+            nphigeo = 2;
+        }
+        phigeo_step = phigeo_length / double(nphigeo);
+
+        // Compute solid-angle for selected Phigeo integration step size
+        double omega0 = 2.0 * std::sin(0.5 * phigeo_step);
+
+        // Integrate over Phigeo
+        double phigeo = phigeo_low;
+        for (int iphigeo = 0; iphigeo <= nphigeo; ++iphigeo, phigeo += phigeo_step) {
+
+            // Compute half length of the arc (in radians) from a circle with
+            // radius phigeo that intersects with the model, defined as a circle
+            // with maximum radius m_theta_max
+            double dphi = 0.5 * gammalib::roi_arclength(phigeo,
+                                                        zeta,
+                                                        cos_zeta,
+                                                        sin_zeta,
+                                                        theta_max,
+                                                        cos_theta_max);
+
+            // If arc length is not positive then skip this Phigeo value
+            if (dphi <= 0.0) {
+                continue;
+            }
+
+            // Compute sine and cosine of Phigeo
+            double sin_phigeo = std::sin(phigeo);
+            double cos_phigeo = std::cos(phigeo);
+
+            // Determine number of azimuthal integration steps and step size
+            // in radians. By default the integration step size is one tenth of
+            // the arclength, but larger than the Phigeo binning of the IAQ. We
+            // also make sure that there are at least two integration steps.
+            double azimuth_length = 2.0 * dphi * sin_phigeo;
+            double az_step        = 0.1 * azimuth_length;
+            if (az_step > phigeo_bin) {
+                az_step = phigeo_bin;
+            }
+            int naz = int(azimuth_length/az_step + 0.5);
+            if (naz < 2) {
+                naz = 2;
+            }
+            az_step = 2.0 * dphi / double(naz);
+
+            // Computes solid angle of integration bin multiplied by Jaccobian
+            double omega = omega0 * az_step * sin_phigeo;
+
+            // Loop over azimuth angle
+            double az = phi0 - dphi;
+            for (int iaz = 0; iaz <= naz; ++iaz, az += az_step) {
+
+                // Get sky direction
+                double sin_az = std::sin(az);
+                double cos_az = std::cos(az);
+                GVector native(-cos_az*sin_phigeo, sin_az*sin_phigeo, cos_phigeo);
+                GVector dir = rot * native;
+                GSkyDir skyDir;
+                skyDir.celvector(dir);
+
+                // Fall through if sky direction is not contained in DRX
+                if (!drx.contains(skyDir)) {
+                    continue;
+                }
+
+                // Set photon
+                GPhoton photon(skyDir, bin->energy(), bin->time());
+
+                // Get model sky intensity for photon (unit: sr^-1)
+                double intensity = model.spatial()->eval(photon);
+
+                // Fall through if intensity is zero
+                if (intensity == 0.0) {
+                    continue;
+                }
+
+                // Multiply intensity by DRX value (unit: cm^2 s sr^-1)
+                intensity *= drx(skyDir);
+
+                // Multiply intensity by solid angle (unit: cm^2 s)
+                intensity *= omega;
+
+                // Extract IAQ value by linear inter/extrapolation in Phigeo
+                double phirat   = phigeo / phigeo_bin;     // 0.5 at bin centre
+                int    i_phigeo = int(phirat);             // index into which Phigeo falls
+                double eps      = phirat - i_phigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
+                double iaq      = 0.0;
+
+                // Continue only if Phigeo index is valid
+                if (i_phigeo < m_phigeo_bins) {
+
+                    // Loop over Phibar
+                    for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
+
+                        // Get IAQ index
+                        int i = iphibar * m_phigeo_bins + i_phigeo;
+
+                        // Get interpolated IAQ value
+                        if (eps < 0.0) { // interpolate towards left
+                            if (i_phigeo > 0) {
+                                iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
+                            }
+                            else {
+                                iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
+                            }
+                        }
+                        else {           // interpolate towards right
+                            if (i_phigeo < m_phigeo_bins-1) {
+                                iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
+                            }
+                            else {
+                                iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
+                            }
+                        }
+
+                        // Fall through if IAQ is not positive
+                        if (iaq <= 0.0) {
+                            continue;
+                        }
+
+                        // Get DRI index
+                        int idri = ipix + iphibar * npix;
+
+                        // Compute IRF value (unit: cm^2)
+                        double irf = iaq * drg[idri] * iaq_norm * intensity;
+
+                        // Add IRF value if it is positive
+                        if (irf > 0.0) {
+                            irfs[idri] += irf;
+                        }
+
+                    } // endfor: looped over Phibar
+
+                } // endif: Phigeo bin was valid
+
+            } // endfor: looped over azimuth angle around locus of IAQ
+
+        } // endfor: looped over Phigeo angles
+
+    } // endfor: looped over Chi and Psi pixels
+
+    // Return IRF vector
+    return irfs;
+}