implement effective area access

knodlseder
1 parent fe605e26
Showing 4 changed files with 125 additions and 12 deletions
include/GLATResponse.hpp
include/GLATResponseTable.hpp
pyext/GLATResponse.i
src/rsp/GLATAeff.cpp
@@ -58,6 +58,9 @@ public:
                 const GSkyDir& srcDir, const double& srcEng,
                 const GSkyDir& instPntDir, const double& instPosAng,
                 const double& time);
+    double aeff(const double& logE, const double& ctheta);
+    void   aeff_ctheta_min(const double& ctheta);
+    double aeff_ctheta_min(void) const;
     // PSF Methods
     double  psf(const GSkyDir& obsDir, const double& obsEng,
@@ -112,7 +115,9 @@ private:
     GVector get_fits_vector(const GFitsHDU* hdu, const std::string& colname, int row = 0);
     // Private Aeff data
-    GLATResponseTable m_aeff_bins;       //!< Aeff energy and cos theta binning
+    GLATResponseTable m_aeff_bins;        //!< Aeff energy and cos theta binning
+    double*           m_aeff;             //!< Aeff array
+    double            m_aeff_ctheta_min;  //!< Minimum valid cos(theta)
     // Private PSF data
     GLATResponseTable m_psf_bins;        //!< PSF energy and cos theta binning
@@ -56,6 +56,10 @@ public:
     double energy(const int& ie) const;
     int    num_energy(void) const { return m_energy_num; }
     int    num_ctheta(void) const { return m_ctheta_num; }
+    double energy_lo(const int& inx) const { return m_energy_lo[inx]; }
+    double energy_hi(const int& inx) const { return m_energy_hi[inx]; }
+    double ctheta_lo(const int& inx) const { return m_ctheta_lo[inx]; }
+    double ctheta_hi(const int& inx) const { return m_ctheta_hi[inx]; }
 private:
     // Methods
@@ -34,26 +34,36 @@ public:
     GLATResponse(const GLATResponse& rsp);
     ~GLATResponse();
-    // Methods
+    // IRF Methods
     double irf(const GSkyDir& obsDir, const double& obsEng,
                const GSkyDir& srcDir, const double& srcEng,
                const GSkyDir& instPntDir, const double& instPosAng,
                const double& time);
-    double  psf(const GSkyDir& obsDir, const double& obsEng,
+    // Aeff Methods
+    double aeff(const GSkyDir& obsDir, const double& obsEng,
                 const GSkyDir& srcDir, const double& srcEng,
                 const GSkyDir& instPntDir, const double& instPosAng,
                 const double& time);
-    double  psf(const double& delta, const double& logE, const double& ctheta);
-    GVector psf(const GVector& delta, const double& logE, const double& ctheta);
-    double aeff(const GSkyDir& obsDir, const double& obsEng,
+    double aeff(const double& logE, const double& ctheta);
+    void   aeff_ctheta_min(const double& ctheta);
+    double aeff_ctheta_min(void) const;
+
+    // PSF Methods
+    double  psf(const GSkyDir& obsDir, const double& obsEng,
                 const GSkyDir& srcDir, const double& srcEng,
                 const GSkyDir& instPntDir, const double& instPosAng,
                 const double& time);
+    double  psf(const double& delta, const double& logE, const double& ctheta);
+    GVector psf(const GVector& delta, const double& logE, const double& ctheta);
+
+    // Edisp Methods
     double edisp(const GSkyDir& obsDir, const double& obsEng,
                  const GSkyDir& srcDir, const double& srcEng,
                  const GSkyDir& instPntDir, const double& instPosAng,
                  const double& time);
-    void   set_caldb(const std::string& caldb);
-    void   load(const std::string& rspname, const std::string& rsptype);
-    void   save(const std::string& rspname) const;
+
+    // Other Methods
+    void set_caldb(const std::string& caldb);
+    void load(const std::string& rspname, const std::string& rsptype);
+    void save(const std::string& rspname) const;
 };
@@ -66,6 +66,50 @@ double GLATResponse::aeff(const GSkyDir&amp; obsDir, const double&amp; obsEng,
 }
+/***********************************************************************//**
+ * @brief Return effective area in units of cm2
+ *
+ * @param[in] logE Logarithm of the true photon energy (MeV).
+ * @param[in] ctheta Cosine of the zenith angle with respect to the pointing 
+ *            axis.
+ ***************************************************************************/
+double GLATResponse::aeff(const double& logE, const double& ctheta)
+{
+    // Get effective area value
+    double aeff = (ctheta >= m_aeff_ctheta_min) 
+                  ? m_aeff_bins.interpolate(logE, ctheta, m_aeff) : 0.0;
+
+    // Return effective area value
+    return aeff;
+}
+
+
+/***********************************************************************//**
+ * @brief Set minimum cos(theta) angle for effective area access
+ *
+ * @param[in] ctheta Cosine of the maximum zenith angle for which effective
+ *                   areas will be returned (0 is returned for larger values)
+ ***************************************************************************/
+void GLATResponse::aeff_ctheta_min(const double& ctheta)
+{
+    // Set minimum cos(theta) value
+    m_aeff_ctheta_min = ctheta;
+    
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get minimum cos(theta) angle for effective area access
+ ***************************************************************************/
+double GLATResponse::aeff_ctheta_min(void) const
+{
+    // Return minimum cos(theta) value
+    return m_aeff_ctheta_min;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                            Private methods                              =
@@ -74,6 +118,8 @@ double GLATResponse::aeff(const GSkyDir&amp; obsDir, const double&amp; obsEng,
 /***********************************************************************//**
  * @brief Initialise effective area
+ *
+ * The effective area is converted into units of cm2.
  ***************************************************************************/
 void GLATResponse::aeff_init(void)
 {
@@ -92,8 +138,23 @@ void GLATResponse::aeff_init(void)
     // Get the energy and cos(theta) bins
     m_aeff_bins.load(hdu);
-    // Get the data
-    GVector effarea = get_fits_vector(hdu, "EFFAREA");
+    // Allocate memory for effective area
+    int size = m_aeff_bins.num_energy() * m_aeff_bins.num_ctheta();
+    if (size > 0)
+        m_aeff = new double[size];
+
+    // Get pointer to effective area column
+    GFitsTableCol* ptr = hdu->column("EFFAREA");
+    if (ptr == NULL)
+        throw GException::fits_column_not_found(G_INIT_AEFF, "EFFAREA");
+
+    // Copy data and convert from m2 into cm2
+    int num = ptr->number();
+    for (int i = 0; i < num; ++i)
+        m_aeff[i] = ptr->real(0,i) * 1.0e4;
+
+    // Set minimum cos(theta)
+    m_aeff_ctheta_min = m_aeff_bins.ctheta_lo(0);
     // Return
     return;
@@ -116,8 +177,20 @@ void GLATResponse::aeff_append(GFits&amp; file) const
     m_aeff_bins.save(&hdu_bounds);
     hdu_bounds.extname("ABOUNDS");
+    // Build effective area image
+    int naxes[] = {m_aeff_bins.num_energy(), m_aeff_bins.num_ctheta()};
+    GFitsDblImage image_aeff(2, naxes, m_aeff);
+
+    // Construct effective area HDU
+    GFitsHDU hdu_aeff(image_aeff);
+    hdu_aeff.extname("AEFF");
+    hdu_aeff.header()->update(GFitsHeaderCard("CTYPE1", "Energy", "Energy binning"));
+    hdu_aeff.header()->update(GFitsHeaderCard("CTYPE2", "cos(theta)", "cos(theta) binning"));
+    hdu_aeff.header()->update(GFitsHeaderCard("BUNIT", "cm2", "Pixel unit"));
+
     // Append HDUs to FITS file
     file.append_hdu(hdu_bounds);
+    file.append_hdu(hdu_aeff);
     // Return
     return;
@@ -130,7 +203,9 @@ void GLATResponse::aeff_append(GFits&amp; file) const
 void GLATResponse::aeff_init_members(void)
 {
     // Initialise Aeff members
-
+    m_aeff_ctheta_min = 0.0;
+    m_aeff            = NULL;
+    
     // Return
     return;
 }
@@ -141,6 +216,19 @@ void GLATResponse::aeff_init_members(void)
  ***************************************************************************/
 void GLATResponse::aeff_copy_members(const GLATResponse& rsp)
 {
+    // Copy attributes
+    m_aeff_bins       = rsp.m_aeff_bins;
+    m_aeff_ctheta_min = rsp.m_aeff_ctheta_min;
+
+    // Copy effective area array
+    if (rsp.m_aeff != NULL) {
+        int size = m_aeff_bins.num_energy() * m_aeff_bins.num_ctheta();
+        if (size > 0) {
+            m_aeff = new double[size];
+            memcpy(m_aeff, rsp.m_aeff, size*sizeof(double));
+        }
+    }
+
     // Return
     return;
 }
@@ -151,6 +239,12 @@ void GLATResponse::aeff_copy_members(const GLATResponse&amp; rsp)
 ***************************************************************************/
 void GLATResponse::aeff_free_members(void)
 {
+    // Free effective area memory
+    if (m_aeff != NULL) delete [] m_aeff;
+
+    // Signal that effective area memory is free
+    m_aeff = NULL;
+
     // Return
     return;
 }