From ba9b66ded85bdbf2a8d77fb8977dbc9a814d8b72 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 18 Jul 2016 16:17:51 +0200 Subject: [PATCH] Fix consistency problem in GCTABackground3D class --- ChangeLog | 9 +++++---- NEWS | 6 +++++- inst/cta/include/GCTABackground3D.hpp | 40 +++++++++++++++++++++++----------------- inst/cta/src/GCTABackground3D.cpp | 654 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- inst/cta/src/GCTAModelIrfBackground.cpp | 2 +- inst/cta/src/GCTAResponseIrf.cpp | 2 +- src/sky/GSkyMap.cpp | 4 ++-- 7 files changed, 453 insertions(+), 264 deletions(-) diff --git a/ChangeLog b/ChangeLog index 10d14e1..c8fe995 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,12 +1,13 @@ -2016-07-12 Juergen Knoedlseder +2016-07-18 Juergen Knoedlseder * Version 1.1.0 released ========================= + Fix consistency problem in GCTABackground3D (#1815) Models accessing file now support relative file names (#1808) Add GXmlElement and GXmlDocument filename() methods (#1808) Move node parent from GXmlElement to GXmlNode class (#1808) - Correct Makefile for example executables (#1806) + Fix Makefile for example executables (#1806) Include examples in unit tests (#1803) Handle additional columns in GArf class GCTAAeff2D no longer extrapolates beyond response table (#1784) @@ -56,8 +57,8 @@ Add GFilename::fullname() method Add GCTAEventList::remove() method Handle also additional columns in GCTAEventList::append() method - Correct GFitsTableCol::insert() method for empty columns - Correct GCTAEdisp2D::compute_ebounds_xxx() methods (#1662) + Fix GFitsTableCol::insert() method for empty columns + Fix GCTAEdisp2D::compute_ebounds_xxx() methods (#1662) Use GFilename in GFits to store filename (#1649) Load CTA GTI extension from the HDU name provided by the DSS keyword (#1598) Support Fermi LAT Pass 8 response functions (#1517) diff --git a/NEWS b/NEWS index d6418c3..5e21137 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0 Author(s) : Juergen Knoedlseder -12 July 2016 +18 July 2016 1. Introduction @@ -440,6 +440,10 @@ handled by the weight array (#1789). The GCTAAeff2D class no longer extrapolates the effective area beyond the boundaries of the response table, but sets the effective area values to zero (#1784). +The GCTABackground3D class now assures that the mc() method covers the same DETX and DETY +range than the access operator. This fixes a consistency problem between Monte Carlo +sampling and function evaluation (#1815). + 21. COMPTEL interface --------------------- diff --git a/inst/cta/include/GCTABackground3D.hpp b/inst/cta/include/GCTABackground3D.hpp index 08e5df1..ca434ba 100644 --- a/inst/cta/include/GCTABackground3D.hpp +++ b/inst/cta/include/GCTABackground3D.hpp @@ -83,33 +83,39 @@ public: private: // Methods - void init_members(void); - void copy_members(const GCTABackground3D& bgd); - void free_members(void); - void init_mc_cache(void) const; + void init_members(void); + void copy_members(const GCTABackground3D& bgd); + void free_members(void); + void set_limits(void); + int index(const int& idetx, const int& idety, const int& iebin) const; + void init_mc_cache(void) const; + void init_mc_max_rate(void) const; + double solid_angle(const double& detx1, const double& dety1, + const double& detx2, const double& dety2, + const double& detx3, const double& dety3) const; // Members GFilename m_filename; //!< Name of background response file GCTAResponseTable m_background; //!< Background response table - double m_mc_max_bin; //!< Maximum spatial binsize for MC - double m_mc_max_logE; //!< Maximum log energy binsize for MC int m_inx_detx; //!< DETX index int m_inx_dety; //!< DETY index int m_inx_energy; //!< Energy index int m_inx_bgd; //!< Background index + int m_num_detx; //!< Number of DETX bins + int m_num_dety; //!< Number of DETY bins + int m_num_energy; //!< Number of energy bins + int m_num[3]; //!< Array of number of bins + double m_detx_min; //!< DETX minimum (radians) + double m_detx_max; //!< DETX maximum (radians) + double m_dety_min; //!< DETY minimum (radians) + double m_dety_max; //!< DETY maximum (radians) + double m_logE_min; //!< Log10(E/TeV) minimum + double m_logE_max; //!< Log10(E/TeV) maximum // Monte Carlo cache - mutable std::vector m_mc_max; //!< Maximum rate - mutable GModelSpectralNodes m_mc_spectrum; //!< Response cube spectrum - mutable double m_mc_detx_min; //!< DETX minimum (deg) - mutable double m_mc_detx_max; //!< DETX maximum (deg) - mutable double m_mc_detx_bin; //!< DETX binsize (deg) - mutable double m_mc_dety_min; //!< DETY minimum (deg) - mutable double m_mc_dety_max; //!< DETY maximum (deg) - mutable double m_mc_dety_bin; //!< DETY binsize (deg) - mutable double m_mc_logE_min; //!< log10 energy minimum (TeV) - mutable double m_mc_logE_max; //!< log10 energy maximum (TeV) - mutable double m_mc_logE_bin; //!< log10 energy binsize (TeV) + mutable std::vector m_mc_max; //!< Maximum background rate + mutable GModelSpectralNodes m_mc_spectrum; //!< Response cube spectrum + mutable double m_mc_one_minus_costheta; //!< 1-cos(theta_max) }; diff --git a/inst/cta/src/GCTABackground3D.cpp b/inst/cta/src/GCTABackground3D.cpp index 0deda46..d6098d1 100644 --- a/inst/cta/src/GCTABackground3D.cpp +++ b/inst/cta/src/GCTABackground3D.cpp @@ -174,27 +174,28 @@ GCTABackground3D& GCTABackground3D::operator=(const GCTABackground3D& bgd) * @brief Return background rate in units of events/s/MeV/sr * * @param[in] logE Log10 of the true photon energy (TeV). - * @param[in] detx Tangential coord in nominal sys (radians). - * @param[in] dety Tangential coord in nominal sys (radians). + * @param[in] detx Tangential X coordinate in nominal system (radians). + * @param[in] dety Tangential Y coordinate in nominal system (radians). * * Returns the background rate in units of events/s/MeV/sr for a given energy - * and detector coordinates. The method assures that the background rate - * never becomes negative. + * and detector coordinates. The operator interpolates linearly in DETX and + * DETY, and logarithmically in energy. * - * The method interpolates linearly in DETX and DETY, and logarithmically - * in energy. - * - * If the background model is not valid the operator returns 0. + * The operator assures that the background rate never becomes negative. For + * invalid background models or detector coordinates outside the range + * covered by the response cube, the operator returns a rate of zero. ***************************************************************************/ double GCTABackground3D::operator()(const double& logE, const double& detx, const double& dety) const { - // Initialise rate - double rate(0.0); + // Initialise background rate + double rate = 0.0; - // Continue only if background model is valid - if (is_valid()) { + // Continue only if background model is valid and if the DETX and DETY + // values are in validity range + if (is_valid() && (detx >= m_detx_min) && (detx <= m_detx_max) && + (dety >= m_dety_min) && (dety <= m_dety_max)) { // Retrieve references to node arrays const GNodeArray& detx_nodes = m_background.axis_nodes(m_inx_detx); @@ -206,46 +207,43 @@ double GCTABackground3D::operator()(const double& logE, dety_nodes.set_value(dety); energy_nodes.set_value(logE); - // Compute offsets of DETY in DETX-DETY plane - int size1 = m_background.axis_bins(m_inx_detx); - int offset_left = dety_nodes.inx_left() * size1; - int offset_right = dety_nodes.inx_right() * size1; - - // Set indices for bi-linear interpolation in DETX-DETY plane - int inx_ll = detx_nodes.inx_left() + offset_left; - int inx_lr = detx_nodes.inx_left() + offset_right; - int inx_rl = detx_nodes.inx_right() + offset_left; - int inx_rr = detx_nodes.inx_right() + offset_right; - // Set weighting factors for bi-linear interpolation in DETX-DETY plane double wgt_ll = detx_nodes.wgt_left() * dety_nodes.wgt_left(); double wgt_lr = detx_nodes.wgt_left() * dety_nodes.wgt_right(); double wgt_rl = detx_nodes.wgt_right() * dety_nodes.wgt_left(); double wgt_rr = detx_nodes.wgt_right() * dety_nodes.wgt_right(); - // Set indices for energy interpolation - int inx_emin = energy_nodes.inx_left(); - int inx_emax = energy_nodes.inx_right(); - // Set weighting factors for energy interpolation double wgt_emin = energy_nodes.wgt_left(); double wgt_emax = energy_nodes.wgt_right(); - // Compute offsets in energy dimension - int npixels = m_background.axis_bins(m_inx_detx) * - m_background.axis_bins(m_inx_dety); - int offset_emin = inx_emin * npixels; - int offset_emax = inx_emax * npixels; - - // Bi-linear interpolation the rates in both energy layers - double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll + offset_emin) + - wgt_lr * m_background(m_inx_bgd, inx_lr + offset_emin) + - wgt_rl * m_background(m_inx_bgd, inx_rl + offset_emin) + - wgt_rr * m_background(m_inx_bgd, inx_rr + offset_emin); - double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll + offset_emax) + - wgt_lr * m_background(m_inx_bgd, inx_lr + offset_emax) + - wgt_rl * m_background(m_inx_bgd, inx_rl + offset_emax) + - wgt_rr * m_background(m_inx_bgd, inx_rr + offset_emax); + // Get bin indices + int inx_ll_emin = index(detx_nodes.inx_left(), dety_nodes.inx_left(), + energy_nodes.inx_left()); + int inx_lr_emin = index(detx_nodes.inx_left(), dety_nodes.inx_right(), + energy_nodes.inx_left()); + int inx_rl_emin = index(detx_nodes.inx_right(), dety_nodes.inx_left(), + energy_nodes.inx_left()); + int inx_rr_emin = index(detx_nodes.inx_right(), dety_nodes.inx_right(), + energy_nodes.inx_left()); + int inx_ll_emax = index(detx_nodes.inx_left(), dety_nodes.inx_left(), + energy_nodes.inx_right()); + int inx_lr_emax = index(detx_nodes.inx_left(), dety_nodes.inx_right(), + energy_nodes.inx_right()); + int inx_rl_emax = index(detx_nodes.inx_right(), dety_nodes.inx_left(), + energy_nodes.inx_right()); + int inx_rr_emax = index(detx_nodes.inx_right(), dety_nodes.inx_right(), + energy_nodes.inx_right()); + + // Bi-linearly interpolate the rates in both energy layers + double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll_emin) + + wgt_lr * m_background(m_inx_bgd, inx_lr_emin) + + wgt_rl * m_background(m_inx_bgd, inx_rl_emin) + + wgt_rr * m_background(m_inx_bgd, inx_rr_emin); + double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll_emax) + + wgt_lr * m_background(m_inx_bgd, inx_lr_emax) + + wgt_rl * m_background(m_inx_bgd, inx_rl_emax) + + wgt_rr * m_background(m_inx_bgd, inx_rr_emax); // If both rates are positive then perform a logarithmic interpolation // in energy @@ -338,12 +336,6 @@ void GCTABackground3D::read(const GFitsTable& table) // Read background table m_background.read(table); - // Get mandatory indices (throw exception if not found) - m_inx_detx = m_background.axis("DETX"); - m_inx_dety = m_background.axis("DETY"); - m_inx_energy = m_background.axis("ENERG"); - m_inx_bgd = m_background.table("BGD"); - // Throw an exception if the table is not three-dimensional if (m_background.axes() != 3) { std::string msg = "Expected three-dimensional background " @@ -354,6 +346,22 @@ void GCTABackground3D::read(const GFitsTable& table) throw GException::invalid_value(G_READ, msg); } + // Get mandatory indices (throw exception if not found) + m_inx_detx = m_background.axis("DETX"); + m_inx_dety = m_background.axis("DETY"); + m_inx_energy = m_background.axis("ENERG"); + m_inx_bgd = m_background.table("BGD"); + + // Get axes dimensions + m_num_detx = m_background.axis_bins(m_inx_detx); + m_num_dety = m_background.axis_bins(m_inx_dety); + m_num_energy = m_background.axis_bins(m_inx_energy); + + // Set dimension array + m_num[m_inx_detx] = m_num_detx; + m_num[m_inx_dety] = m_num_dety; + m_num[m_inx_energy] = m_num_energy; + // Set DETX and DETY axis to radians m_background.axis_radians(m_inx_detx); m_background.axis_radians(m_inx_dety); @@ -361,6 +369,9 @@ void GCTABackground3D::read(const GFitsTable& table) // Set energy axis to logarithmic scale m_background.axis_log10(m_inx_energy); + // Set DETX and DETY limits + set_limits(); + // Return return; } @@ -484,7 +495,7 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy& energy, const GTime& time, GRan& ran) const { - // Initialise Monte Carlo Cache + // If Monte Carlo has not been initialised then initialise the cache now if (m_mc_max.empty()) { init_mc_cache(); } @@ -508,35 +519,35 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy& energy, int inx_left = energy_nodes.inx_left(); int inx_right = energy_nodes.inx_right(); - // Get maximum background rates - double max_rate_left = m_mc_max[inx_left]; - double max_rate_right = m_mc_max[inx_right]; - - // Compute maximum background rate - double max_rate = 0.0; - if (max_rate_left > 0.0 && max_rate_right > 0.0) { - double wgt_left = energy_nodes.wgt_left(); - double wgt_right = energy_nodes.wgt_right(); - max_rate = std::exp(wgt_left * std::log(max_rate_left) + - wgt_right * std::log(max_rate_right)); - } - else { - max_rate = 0.0; + // Get maximum background rate as the maximum of the left and the right + // energy node + double max_rate = m_mc_max[inx_left]; + if (m_mc_max[inx_right] > max_rate) { + max_rate = m_mc_max[inx_right]; } - // Get instrument direction + // Get instrument direction using a rejection method. This assures + // that the Monte Carlo sample follows the model distribution while (true) { - // Get randomized detx and dety (radians) - double detx = (m_mc_detx_min + - ran.uniform() * (m_mc_detx_max - m_mc_detx_min)) * - gammalib::deg2rad; - double dety = (m_mc_dety_min + - ran.uniform() * (m_mc_dety_max - m_mc_dety_min)) * - gammalib::deg2rad; + // Get randomized detx and dety in radians. If the direction falls + // outside the validity range then get another pair of detx and + // dety values. + double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_costheta); + double phi = gammalib::twopi * ran.uniform(); + double detx = theta * std::cos(phi); + double dety = theta * std::sin(phi); + if ((detx < m_detx_min) || (detx > m_detx_max) || + (dety < m_dety_min) || (dety > m_dety_max)) { + continue; + } - // Get background rate for these coordinates + // Get background rate for these coordinates. If the background + // rate is zero then get other coordinates. double value = this->operator()(logE, detx, dety); + if (value <= 0.0) { + continue; + } // Get uniform random number double uniform = ran.uniform() * max_rate; @@ -581,37 +592,30 @@ std::string GCTABackground3D::print(const GChatter& chatter) const if (is_valid()) { // Compute DETX boundaries in deg - double detx_min = m_background.axis_lo(m_inx_detx,0); - double detx_max = m_background.axis_hi(m_inx_detx, - m_background.axis_bins(m_inx_detx)-1); + double detx_min = m_background.axis_lo(m_inx_detx, 0); + double detx_max = m_background.axis_hi(m_inx_detx, m_num_detx-1); // Compute DETY boundaries in deg double dety_min = m_background.axis_lo(m_inx_dety,0); - double dety_max = m_background.axis_hi(m_inx_dety, - m_background.axis_bins(m_inx_dety)-1); + double dety_max = m_background.axis_hi(m_inx_dety, m_num_dety-1); // Compute energy boundaries in TeV - double emin = m_background.axis_lo(m_inx_energy,0); - double emax = m_background.axis_hi(m_inx_energy, - m_background.axis_bins(m_inx_energy)-1); + double emin = m_background.axis_lo(m_inx_energy, 0); + double emax = m_background.axis_hi(m_inx_energy, m_num_energy-1); // Append information result.append("\n"+gammalib::parformat("Number of DETX bins") + - gammalib::str(m_background.axis_bins(m_inx_detx))); + gammalib::str(m_num_detx)); result.append("\n"+gammalib::parformat("Number of DETY bins") + - gammalib::str(m_background.axis_bins(m_inx_dety))); + gammalib::str(m_num_dety)); result.append("\n"+gammalib::parformat("Number of energy bins") + - gammalib::str(m_background.axis_bins(m_inx_energy))); + gammalib::str(m_num_energy)); result.append("\n"+gammalib::parformat("DETX range")); result.append(gammalib::str(detx_min)+" - "+gammalib::str(detx_max)+" deg"); result.append("\n"+gammalib::parformat("DETX range")); result.append(gammalib::str(dety_min)+" - "+gammalib::str(dety_max)+" deg"); result.append("\n"+gammalib::parformat("Energy range")); result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV"); - result.append("\n"+gammalib::parformat("Maximum bin size for MC")); - result.append(gammalib::str(m_mc_max_bin)+" deg"); - result.append("\n"+gammalib::parformat("Maximum logE step for MC")); - result.append(gammalib::str(m_mc_max_logE)+"^10 TeV"); } // endif: there were 3 axis @@ -646,25 +650,27 @@ void GCTABackground3D::init_members(void) // Initialise members m_filename.clear(); m_background.clear(); - m_mc_max_bin = 0.05; //!< Spatial binning not worse than 0.05 deg - m_mc_max_logE = 0.02; //!< Spectral binning not worse than 0.02^10 TeV - m_inx_detx = 0; - m_inx_dety = 1; - m_inx_energy = 2; - m_inx_bgd = 0; + m_inx_detx = 0; + m_inx_dety = 1; + m_inx_energy = 2; + m_inx_bgd = 0; + m_num_detx = 0; + m_num_dety = 0; + m_num_energy = 0; + m_num[0] = 0; + m_num[1] = 0; + m_num[2] = 0; + m_detx_min = 0.0; + m_detx_max = 0.0; + m_dety_min = 0.0; + m_dety_max = 0.0; + m_logE_min = 0.0; + m_logE_max = 0.0; // Initialise MC cache m_mc_max.clear(); m_mc_spectrum.clear(); - m_mc_detx_min = 0.0; - m_mc_detx_max = 0.0; - m_mc_detx_bin = 0.0; - m_mc_dety_min = 0.0; - m_mc_dety_max = 0.0; - m_mc_dety_bin = 0.0; - m_mc_logE_min = 0.0; - m_mc_logE_max = 0.0; - m_mc_logE_bin = 0.0; + m_mc_one_minus_costheta = 0.0; // Return return; @@ -679,27 +685,29 @@ void GCTABackground3D::init_members(void) void GCTABackground3D::copy_members(const GCTABackground3D& bgd) { // Copy members - m_filename = bgd.m_filename; - m_background = bgd.m_background; - m_mc_max_bin = bgd.m_mc_max_bin; - m_mc_max_logE = bgd.m_mc_max_logE; - m_inx_detx = bgd.m_inx_detx; - m_inx_dety = bgd.m_inx_dety; - m_inx_energy = bgd.m_inx_energy; - m_inx_bgd = bgd.m_inx_bgd; + m_filename = bgd.m_filename; + m_background = bgd.m_background; + m_inx_detx = bgd.m_inx_detx; + m_inx_dety = bgd.m_inx_dety; + m_inx_energy = bgd.m_inx_energy; + m_inx_bgd = bgd.m_inx_bgd; + m_num_detx = bgd.m_num_detx; + m_num_dety = bgd.m_num_dety; + m_num_energy = bgd.m_num_energy; + m_num[0] = bgd.m_num[0]; + m_num[1] = bgd.m_num[1]; + m_num[2] = bgd.m_num[2]; + m_detx_min = bgd.m_detx_min; + m_detx_max = bgd.m_detx_max; + m_dety_min = bgd.m_dety_min; + m_dety_max = bgd.m_dety_max; + m_logE_min = bgd.m_logE_min; + m_logE_max = bgd.m_logE_max; // Copy MC cache - m_mc_max = bgd.m_mc_max; - m_mc_spectrum = bgd.m_mc_spectrum; - m_mc_detx_min = bgd.m_mc_detx_min; - m_mc_detx_max = bgd.m_mc_detx_max; - m_mc_detx_bin = bgd.m_mc_detx_bin; - m_mc_dety_min = bgd.m_mc_dety_min; - m_mc_dety_max = bgd.m_mc_dety_max; - m_mc_dety_bin = bgd.m_mc_dety_bin; - m_mc_logE_min = bgd.m_mc_logE_min; - m_mc_logE_max = bgd.m_mc_logE_max; - m_mc_logE_bin = bgd.m_mc_logE_bin; + m_mc_max = bgd.m_mc_max; + m_mc_spectrum = bgd.m_mc_spectrum; + m_mc_one_minus_costheta = bgd.m_mc_one_minus_costheta; // Return return; @@ -707,6 +715,66 @@ void GCTABackground3D::copy_members(const GCTABackground3D& bgd) /***********************************************************************//** + * @brief Set DETX, DETY and log10(E/TeV) limits + * + * Set the DETX, DETY and log10(E/TeV) limits. + ***************************************************************************/ +void GCTABackground3D::set_limits(void) +{ + // Compute DETX boundaries in radians + m_detx_min = m_background.axis_lo(m_inx_detx, 0) * + gammalib::deg2rad; + m_detx_max = m_background.axis_hi(m_inx_detx, m_num_detx-1) * + gammalib::deg2rad; + + // Compute DETY boundaries in radians + m_dety_min = m_background.axis_lo(m_inx_dety, 0) * + gammalib::deg2rad; + m_dety_max = m_background.axis_hi(m_inx_dety, m_num_dety-1) * + gammalib::deg2rad; + + // Compute energy boundaries in log10(TeV) + GEnergy emin(m_background.axis_lo(m_inx_energy, 0), + m_background.axis_lo_unit(m_inx_energy)); + GEnergy emax(m_background.axis_hi(m_inx_energy, m_num_energy-1), + m_background.axis_hi_unit(m_inx_energy)); + m_logE_min = emin.log10TeV(); + m_logE_max = emax.log10TeV(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return background rate bin index + * + * @param[in] idetx DETX index. + * @param[in] idety DETY index. + * @param[in] iebin Energy index. + * @return Background rate bin index. + * + * Returns the background rate bin index independent of the ordering. + ***************************************************************************/ +int GCTABackground3D::index(const int& idetx, + const int& idety, + const int& iebin) const +{ + // Set index arrays + int inx[3]; + inx[m_inx_detx] = idetx; + inx[m_inx_dety] = idety; + inx[m_inx_energy] = iebin; + + // Compute index + int index = inx[0] + (inx[1] + inx[2] * m_num[1]) * m_num[0]; + + // Return index + return index; +} + + +/***********************************************************************//** * @brief Delete class members ***************************************************************************/ void GCTABackground3D::free_members(void) @@ -719,15 +787,7 @@ void GCTABackground3D::free_members(void) /***********************************************************************//** * @brief Initialise Monte Carlo cache * - * @exception GException::invalid_value - * No valid background model defined. - * * Initialises the cache for Monte Carlo sampling. - * - * @todo Verify assumption made about the solid angles of the response table - * elements. - * @todo Implement method that computes the "flux" in a background map - * pixel (inspired from GSkyMap::flux()) ***************************************************************************/ void GCTABackground3D::init_mc_cache(void) const { @@ -735,125 +795,243 @@ void GCTABackground3D::init_mc_cache(void) const m_mc_max.clear(); m_mc_spectrum.clear(); - // Determine number of response cube pixels and maps - int nx = m_background.axis_bins(m_inx_detx); - int ny = m_background.axis_bins(m_inx_dety); - int npix = nx * ny; - int nmaps = m_background.axis_bins(m_inx_energy); - - // Compute DETX boundaries and binsize in degress - m_mc_detx_min = m_background.axis_lo(m_inx_detx,0); - m_mc_detx_max = m_background.axis_hi(m_inx_detx,nx-1); - m_mc_detx_bin = (m_mc_detx_max - m_mc_detx_min) / nx; + // Initialise maximum rate array + init_mc_max_rate(); + + // Compute 1-cos(theta_max) + double theta1 = std::sqrt(m_detx_min * m_detx_min + m_dety_min * m_dety_min); + double theta2 = std::sqrt(m_detx_min * m_detx_min + m_dety_max * m_dety_max); + double theta3 = std::sqrt(m_detx_max * m_detx_max + m_dety_min * m_dety_min); + double theta4 = std::sqrt(m_detx_max * m_detx_max + m_dety_max * m_dety_max); + double theta = theta1; + if (theta2 > theta) { + theta = theta2; + } + if (theta3 > theta) { + theta = theta3; + } + if (theta4 > theta) { + theta = theta4; + } + m_mc_one_minus_costheta = 1.0 - std::cos(theta); - // Compute DETY boundaries and binsize in degress - m_mc_dety_min = m_background.axis_lo(m_inx_dety,0); - m_mc_dety_max = m_background.axis_hi(m_inx_dety,ny-1); - m_mc_dety_bin = (m_mc_dety_max - m_mc_dety_min) / ny; + // Set number of energy bins. Make sure that we compute at least 100 energy + // bins + int nbins = m_num_energy * 10; + if (nbins < 100) { + nbins = 100; + } - // Compute energy boundaries and binsize in log10(TeV) - GEnergy emin(m_background.axis_lo(m_inx_energy,0), - m_background.axis_lo_unit(m_inx_energy)); - GEnergy emax(m_background.axis_hi(m_inx_energy,nmaps-1), - m_background.axis_hi_unit(m_inx_energy)); - m_mc_logE_min = emin.log10TeV(); - m_mc_logE_max = emax.log10TeV(); - m_mc_logE_bin = (m_mc_logE_max - m_mc_logE_min) / nmaps; - - // Determine solid angle of the pixel. We assume here simply that we - // have square pixels of identical solid angle. It needs to be checked - // whether this is a valid assumption. - double solidangle = m_mc_detx_bin * gammalib::deg2rad * - m_mc_dety_bin * gammalib::deg2rad; - - // Compute flux normalisation factor. As we add-up 3 intensity we have - // do divide by 3, and as we add up 8 wedges we have to divide the - // solid angle by 8. - double flux_norm = gammalib::onethird * solidangle / 8.0; - - // Continue only if there are pixels and maps - if (npix > 0 && nmaps > 0) { - - // Loop over all maps - for (int i = 0; i < nmaps; ++i) { - - // Get logE - double logE = m_mc_logE_min + (i+0.5) * m_mc_logE_bin; - - // Initialise cache with cumulative pixel fluxes and compute - // total flux in response table for normalization. Negative - // pixels are excluded from the cumulative map. - double total_flux = 0.0; // units: events/s/MeV - double max_rate = 0.0; // units: events/s/sr/MeV - double xmin = m_mc_detx_min * gammalib::deg2rad; - double xbin = m_mc_detx_bin * gammalib::deg2rad; - double ymin = m_mc_dety_min * gammalib::deg2rad; - double ybin = m_mc_dety_bin * gammalib::deg2rad; - double dx = 0.5 * xbin; - double dy = 0.5 * ybin; - double detx = xmin + 0.5 * xbin; - for (int ix = 0; ix < nx; ++ix, detx += xbin) { - double dety = ymin + 0.5 * ybin; - for (int iy = 0; iy < ny; ++iy, dety += ybin) { - - // Compute flux by summing the flux in 8 pixel wedges - double i0 = (*this)(logE, detx, dety); - double i1 = (*this)(logE, detx-dx, dety-dy); - double i2 = (*this)(logE, detx, dety-dy); - double i3 = (*this)(logE, detx+dx, dety-dy); - double i4 = (*this)(logE, detx+dx, dety); - double i5 = (*this)(logE, detx+dx, dety+dy); - double i6 = (*this)(logE, detx, dety+dy); - double i7 = (*this)(logE, detx-dx, dety+dy); - double i8 = (*this)(logE, detx-dx, dety); - double flux1 = flux_norm * (i1 + i2 + i0); - double flux2 = flux_norm * (i2 + i3 + i0); - double flux3 = flux_norm * (i3 + i4 + i0); - double flux4 = flux_norm * (i4 + i5 + i0); - double flux5 = flux_norm * (i5 + i6 + i0); - double flux6 = flux_norm * (i6 + i7 + i0); - double flux7 = flux_norm * (i7 + i8 + i0); - double flux8 = flux_norm * (i8 + i1 + i0); - double flux = (flux1 + flux2 + flux3 + flux4 + - flux5 + flux6 + flux7 + flux8); - - // Sum flux - if (flux > 0.0) { - total_flux += flux; - } - - // Get maximum rate - if (i0 > max_rate) { - max_rate = i0; - } + // Compute DETX and DETY binsize in radians and energy binsize in log10(TeV) + double detx_bin = (m_detx_max - m_detx_min) / m_num_detx; + double dety_bin = (m_dety_max - m_dety_min) / m_num_dety; + double logE_bin = (m_logE_max - m_logE_min) / (nbins+1); + + // Loop over all energy bins + for (int i = 0; i < nbins; ++i) { + + // Get logE of energy bin + double logE = m_logE_min + i * logE_bin; + + // Initialise cache with cumulative pixel fluxes and compute total flux + // in response table for normalization. Negative pixels are excluded + // from the cumulative map. + double total_flux = 0.0; // units: events/s/MeV + double dx = 0.5 * detx_bin; + double dy = 0.5 * dety_bin; + double detx = m_detx_min + dx; + for (int ix = 0; ix < m_num_detx; ++ix, detx += detx_bin) { + double dety = m_dety_min + dy; + for (int iy = 0; iy < m_num_dety; ++iy, dety += dety_bin) { + + // Compute intensities + double i0 = (*this)(logE, detx, dety); + double i1 = (*this)(logE, detx-dx, dety-dy); + double i2 = (*this)(logE, detx, dety-dy); + double i3 = (*this)(logE, detx+dx, dety-dy); + double i4 = (*this)(logE, detx+dx, dety); + double i5 = (*this)(logE, detx+dx, dety+dy); + double i6 = (*this)(logE, detx, dety+dy); + double i7 = (*this)(logE, detx-dx, dety+dy); + double i8 = (*this)(logE, detx-dx, dety); + + // Compute solid angle of the 8 pixel wedges + double s1 = solid_angle(detx, dety, detx-dx, dety-dy, detx, dety-dy); + double s2 = solid_angle(detx, dety, detx, dety-dy, detx+dx, dety-dy); + double s3 = solid_angle(detx, dety, detx+dx, dety-dy, detx+dx, dety); + double s4 = solid_angle(detx, dety, detx+dx, dety, detx+dx, dety+dy); + double s5 = solid_angle(detx, dety, detx+dx, dety+dy, detx, dety+dy); + double s6 = solid_angle(detx, dety, detx, dety+dy, detx-dx, dety+dy); + double s7 = solid_angle(detx, dety, detx-dx, dety+dy, detx-dx, dety); + double s8 = solid_angle(detx, dety, detx-dx, dety, detx-dx, dety-dy); + + // Compute flux by summing the flux in 8 pixel wedges + double flux1 = s1 * (i1 + i2 + i0); + double flux2 = s2 * (i2 + i3 + i0); + double flux3 = s3 * (i3 + i4 + i0); + double flux4 = s4 * (i4 + i5 + i0); + double flux5 = s5 * (i5 + i6 + i0); + double flux6 = s6 * (i6 + i7 + i0); + double flux7 = s7 * (i7 + i8 + i0); + double flux8 = s8 * (i8 + i1 + i0); + double flux = (flux1 + flux2 + flux3 + flux4 + + flux5 + flux6 + flux7 + flux8) / 3.0; + + // Sum flux + if (flux > 0.0) { + total_flux += flux; } - } - // Append maximum rate - m_mc_max.push_back(max_rate); + } // endfor: loop over DETY - // Set energy value (unit independent) - GEnergy energy; - energy.log10(logE, m_background.axis_lo_unit(m_inx_energy)); + } // endfor: loop over DETX - // Only append node if flux is positive - if (total_flux > 0.0) { - m_mc_spectrum.append(energy, total_flux); - } + // Set energy value (unit independent) + GEnergy energy; + energy.log10(logE, m_background.axis_lo_unit(m_inx_energy)); - // Dump spectrum for debugging - #if defined(G_DEBUG_MC_INIT) - std::cout << "Energy=" << energy; - std::cout << " Rate_max=" << max_rate; - std::cout << " events/s/sr/MeV"; - std::cout << " Rate_total=" << total_flux; - std::cout << " events/s/MeV" << std::endl; - #endif + // Only append node if flux is positive + if (total_flux > 0.0) { + m_mc_spectrum.append(energy, total_flux); + } - } // endfor: looped over all maps + // Dump spectrum for debugging + #if defined(G_DEBUG_MC_INIT) + std::cout << "Energy=" << energy; + std::cout << " Rate_total=" << total_flux; + std::cout << " events/s/MeV" << std::endl; + #endif + + } // endfor: looped over energy bins + + // Dump maximum rate + #if defined(G_DEBUG_MC_INIT) + for (int i = 0; i < m_mc_max.size(); ++i) { + double logE = m_logE_min + (i+0.5) * logE_bin; + GEnergy energy; + energy.log10(logE, m_background.axis_lo_unit(m_inx_energy)); + std::cout << "Energy=" << energy; + std::cout << " Rate_max=" << max_rate; + std::cout << " events/s/sr/MeV"; + } + #endif - } // endif: there were cube pixels and maps + // Return + return; +} + + +/***********************************************************************//** + * @brief Initialise array of maximum background rate + * + * Initialises the array m_mc_max that holds the maximum background rate in + * each DETX-DETY plane of the background model. + ***************************************************************************/ +void GCTABackground3D::init_mc_max_rate(void) const +{ + // Initialise cache + m_mc_max.clear(); + + // Loop over all maps to determine the maximum map intensity in units + // of events/s/sr/MeV + for (int iebin = 0; iebin < m_num_energy; ++iebin) { + + // Initialise maximum rate + double max_rate = 0.0; + + // Loop over DETX-DETY plane + for (int ix = 0; ix < m_num_detx; ++ix) { + for (int iy = 0; iy < m_num_dety; ++iy) { + + // Get bin index + int inx = index(ix, iy, iebin); + + // Get background rate + double rate = m_background(m_inx_bgd, inx); + + // If background rate is larger than maximum then store the + // background rate as new maximum + if (rate > max_rate) { + max_rate = rate; + } + + } // endfor: looped over DETY + } // endfor: looped over DETX + + // Append maximum rate + m_mc_max.push_back(max_rate); + + } // endif: looped over all maps // Return return; } + + +/***********************************************************************//** + * @brief Compute solid angle of pixel wedge + * + * @param[in] detx1 DETX of first edge. + * @param[in] dety1 DETY of first edge. + * @param[in] detx2 DETX of second edge. + * @param[in] dety2 DETY of second edge. + * @param[in] detx3 DETX of third edge. + * @param[in] dety3 DETY of third edge. + * @return Solid angle (steradians). + * + * Estimate the solid angle subtended by 3 coordinates using Huilier's + * theorem. + * + * Below, the definiton of the pixel cornes and sides are shown as used + * within the code. + * + * a12 + * 1---------2 + * | / + * | / + * | / + * | / + * a13| /a23 + * | / + * | / + * | / + * |/ + * 3 + * + ***************************************************************************/ +double GCTABackground3D::solid_angle(const double& detx1, const double& dety1, + const double& detx2, const double& dety2, + const double& detx3, const double& dety3) const +{ + // Initialise solid angle + double solidangle = 0.0; + + // Set sky directions + double theta1 = std::sqrt(detx1 * detx1 + dety1 * dety1); + double phi1 = std::atan2(dety1, detx1); + double theta2 = std::sqrt(detx2 * detx2 + dety2 * dety2); + double phi2 = std::atan2(dety2, detx2); + double theta3 = std::sqrt(detx3 * detx3 + dety3 * dety3); + double phi3 = std::atan2(dety3, detx3); + GSkyDir dir1; + GSkyDir dir2; + GSkyDir dir3; + dir1.radec(phi1, gammalib::pihalf - theta1); + dir2.radec(phi2, gammalib::pihalf - theta2); + dir3.radec(phi3, gammalib::pihalf - theta3); + + // Compute angular distances between pixel corners + double a12 = dir1.dist(dir2); + double a13 = dir1.dist(dir3); + double a23 = dir2.dist(dir3); + + // Compute solid angle + double s = 0.5 * (a12 + a23 + a13); + solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) * + std::tan(0.5*(s-a12)) * + std::tan(0.5*(s-a23)) * + std::tan(0.5*(s-a13)))); + + // Return solid angle + return solidangle; +} diff --git a/inst/cta/src/GCTAModelIrfBackground.cpp b/inst/cta/src/GCTAModelIrfBackground.cpp index 819abe4..fd6737c 100644 --- a/inst/cta/src/GCTAModelIrfBackground.cpp +++ b/inst/cta/src/GCTAModelIrfBackground.cpp @@ -659,7 +659,7 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation& obs, GRan& ran) co for (int i = 0; i < spectral.nodes(); ++i) { GEnergy energy = spectral.energy(i); double intensity = spectral.intensity(i); - double norm = m_spectral->eval(energy, events->tstart()); + double norm = m_spectral->eval(energy); spectral.intensity(i, norm*intensity); } diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index 544eba4..c9a7ce1 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -634,7 +634,7 @@ GCTAEventAtom* GCTAResponseIrf::mc(const double& area, const GPhoton& photon, } // Continue only if event is detected - if (ran.uniform() <= ulimite) { + if ((ulimite > 0.0) && (ran.uniform() <= ulimite)) { // Apply deadtime correction double deadc = obs.deadc(photon.time()); diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 86033ae..ca8c42d 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -2802,7 +2802,7 @@ GFitsImageDouble* GSkyMap::create_wcs_hdu(void) const * @return Solid angle (steradians). * * Estimate the solid angle subtended by 4 sky directions using Huilier's - * theroem. + * theorem. * * Below, the definiton of the pixel cornes and sides are shown as used * within the code. @@ -2891,7 +2891,7 @@ double GSkyMap::solidangle(const GSkyDir& dir1, const GSkyDir& dir2, * @return Solid angle (steradians). * * Estimate the solid angle subtended by 3 sky directions using Huilier's - * theroem. + * theorem. * * Below, the definiton of the pixel cornes and sides are shown as used * within the code. -- libgit2 0.21.4