diff --git a/include/GLATResponse.hpp b/include/GLATResponse.hpp index 08f7c7e..070b862 100644 --- a/include/GLATResponse.hpp +++ b/include/GLATResponse.hpp @@ -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 diff --git a/include/GLATResponseTable.hpp b/include/GLATResponseTable.hpp index fb3e25a..19116ed 100644 --- a/include/GLATResponseTable.hpp +++ b/include/GLATResponseTable.hpp @@ -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 diff --git a/pyext/GLATResponse.i b/pyext/GLATResponse.i index 8117d0f..a5fd7d8 100644 --- a/pyext/GLATResponse.i +++ b/pyext/GLATResponse.i @@ -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; }; diff --git a/src/rsp/GLATAeff.cpp b/src/rsp/GLATAeff.cpp index b8e17bb..a8a976e 100644 --- a/src/rsp/GLATAeff.cpp +++ b/src/rsp/GLATAeff.cpp @@ -66,6 +66,50 @@ double GLATResponse::aeff(const GSkyDir& obsDir, const double& 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& obsDir, const double& 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& 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& 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& 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; }