Lili Yang / gammalib

Commit f8c9f87fba951335f67e876df163012720e1b301

Authored by Jürgen Knödlseder
1 parent b86bdaee

Implement internal linear delta binning. 

Linear delta binning accelerates the interpolation scheme.
Inline Side-by-side
Showing 2 changed files with 93 additions and 24 deletions
inst/cta/include/GCTAMeanPsf.hpp
  View file @f8c9f87
... ... @@ -100,15 +100,18 @@ protected:
100 100 void free_members(void);
101 101 void clear_cube(void);
102 102 void update(const double& delta, const double& logE) const;
  103 + void set_delta_axis(void);
103 104 void set_eng_axis(void);
104 105 void set_to_smooth(void);
105 106  
106 107 // Data
107   - mutable std::string m_filename; //!< Filename
108   - GSkymap m_cube; //!< PSF cube
109   - GEbounds m_ebounds; //!< Energy bounds for the PSF cube
110   - GNodeArray m_elogmeans; //!< Mean log10TeV energy for the PSF cube
111   - GNodeArray m_deltas; //!< Delta bins (deg) for the PSF cube
  108 + mutable std::string m_filename; //!< Filename
  109 + GSkymap m_cube; //!< PSF cube
  110 + GEbounds m_ebounds; //!< Energy bounds for the PSF cube
  111 + GNodeArray m_elogmeans; //!< Mean log10TeV energy for the PSF cube
  112 + GNodeArray m_deltas; //!< Delta bins (deg) for the PSF cube
  113 + GNodeArray m_deltas_cache; //!< Internal delta bins (rad)
  114 + bool m_quadratic_binning; //!< Internal binning is linear
112 115  
113 116 private:
114 117 // Response table computation cache for 2D access
... ...
inst/cta/src/GCTAMeanPsf.cpp
  View file @f8c9f87
... ... @@ -42,8 +42,8 @@
42 42 /* __ Macros _____________________________________________________________ */
43 43  
44 44 /* __ Coding definitions _________________________________________________ */
45   -#define G_SMOOTH_PSF //!< Guarantee no singularities
46   -#define G_SQRT_BINNING //!< Use sqrt delta spacing (samples better the peak)
  45 +#define G_SMOOTH_PSF //!< Guarantee no singularities
  46 +#define G_QUADRATIC_BINNING //!< Use quadratic delta spacing
47 47  
48 48 /* __ Debug definitions __________________________________________________ */
49 49  
... ... @@ -149,7 +149,7 @@ GCTAMeanPsf::GCTAMeanPsf(const std::string&amp; wcs,
149 149 // Set delta node array
150 150 m_deltas.clear();
151 151 for (int i = 0; i < ndbins; ++i) {
152   - #if defined(G_SQRT_BINNING)
  152 + #if defined(G_QUADRATIC_BINNING)
153 153 double binsize = std::sqrt(dmax) / double(ndbins);
154 154 double delta = binsize * (double(i) + 0.5); // avoid central singularity
155 155 delta *= delta;
... ... @@ -505,6 +505,9 @@ void GCTAMeanPsf::read(const GFits&amp; fits)
505 505 // Read delta nodes
506 506 m_deltas.read(hdu_deltas);
507 507  
  508 + // Set delta node array for computation
  509 + set_delta_axis();
  510 +
508 511 // Set energy node array
509 512 set_eng_axis();
510 513  
... ... @@ -669,6 +672,8 @@ void GCTAMeanPsf::init_members(void)
669 672 m_ebounds.clear();
670 673 m_elogmeans.clear();
671 674 m_deltas.clear();
  675 + m_deltas_cache.clear();
  676 + m_quadratic_binning = false;
672 677  
673 678 // Initialise cache
674 679 m_inx1 = 0;
... ... @@ -693,11 +698,13 @@ void GCTAMeanPsf::init_members(void)
693 698 void GCTAMeanPsf::copy_members(const GCTAMeanPsf& cube)
694 699 {
695 700 // Copy members
696   - m_filename = cube.m_filename;
697   - m_cube = cube.m_cube;
698   - m_ebounds = cube.m_ebounds;
699   - m_elogmeans = cube.m_elogmeans;
700   - m_deltas = cube.m_deltas;
  701 + m_filename = cube.m_filename;
  702 + m_cube = cube.m_cube;
  703 + m_ebounds = cube.m_ebounds;
  704 + m_elogmeans = cube.m_elogmeans;
  705 + m_deltas = cube.m_deltas;
  706 + m_deltas_cache = cube.m_deltas_cache;
  707 + m_quadratic_binning = cube.m_quadratic_binning;
701 708  
702 709 // Copy cache
703 710 m_inx1 = cube.m_inx1;
... ... @@ -749,28 +756,87 @@ void GCTAMeanPsf::clear_cube(void)
749 756 * @brief Update PSF parameter cache
750 757 *
751 758 * @param[in] delta Angular separation between true and measured photon
752   - * directions (rad).
  759 + * directions (radians).
753 760 * @param[in] logE Log10 true photon energy (TeV).
754 761 *
755 762 * This method updates the PSF parameter cache.
756 763 ***************************************************************************/
757 764 void GCTAMeanPsf::update(const double& delta, const double& logE) const
758 765 {
759   - // Set node array interpolation values
760   - m_deltas.set_value(delta*gammalib::rad2deg);
  766 + // Set node array for delta interpolation
  767 + if (m_quadratic_binning) {
  768 + m_deltas_cache.set_value(std::sqrt(delta));
  769 + }
  770 + else {
  771 + m_deltas_cache.set_value(delta);
  772 + }
  773 +
  774 + // Set node array for energy interpolation
761 775 m_elogmeans.set_value(logE);
762 776  
763 777 // Set indices for bi-linear interpolation
764   - m_inx1 = offset(m_deltas.inx_left(), m_elogmeans.inx_left());
765   - m_inx2 = offset(m_deltas.inx_left(), m_elogmeans.inx_right());
766   - m_inx3 = offset(m_deltas.inx_right(), m_elogmeans.inx_left());
767   - m_inx4 = offset(m_deltas.inx_right(), m_elogmeans.inx_right());
  778 + m_inx1 = offset(m_deltas_cache.inx_left(), m_elogmeans.inx_left());
  779 + m_inx2 = offset(m_deltas_cache.inx_left(), m_elogmeans.inx_right());
  780 + m_inx3 = offset(m_deltas_cache.inx_right(), m_elogmeans.inx_left());
  781 + m_inx4 = offset(m_deltas_cache.inx_right(), m_elogmeans.inx_right());
768 782  
769 783 // Set weighting factors for bi-linear interpolation
770   - m_wgt1 = m_deltas.wgt_left() * m_elogmeans.wgt_left();
771   - m_wgt2 = m_deltas.wgt_left() * m_elogmeans.wgt_right();
772   - m_wgt3 = m_deltas.wgt_right() * m_elogmeans.wgt_left();
773   - m_wgt4 = m_deltas.wgt_right() * m_elogmeans.wgt_right();
  784 + m_wgt1 = m_deltas_cache.wgt_left() * m_elogmeans.wgt_left();
  785 + m_wgt2 = m_deltas_cache.wgt_left() * m_elogmeans.wgt_right();
  786 + m_wgt3 = m_deltas_cache.wgt_right() * m_elogmeans.wgt_left();
  787 + m_wgt4 = m_deltas_cache.wgt_right() * m_elogmeans.wgt_right();
  788 +
  789 + // Return
  790 + return;
  791 +}
  792 +
  793 +
  794 +/***********************************************************************//**
  795 + * @brief Set nodes for delta axis in radians
  796 + *
  797 + * Set the delta axis nodes in radians. If a quadratic binning is detected,
  798 + * the axis is set to sqrt(delta) so that a linear interpolation scheme
  799 + * can be used which is much faster than a quadratic scheme.
  800 + *
  801 + * @todo Check that none of the axis boundaries is non-positive.
  802 + ***************************************************************************/
  803 +void GCTAMeanPsf::set_delta_axis(void)
  804 +{
  805 + // Initialise computation members
  806 + m_deltas_cache.clear();
  807 + m_quadratic_binning = false;
  808 +
  809 + // Set up a linear array that represents the square root of the delta
  810 + // values. If this array corresponds within some numerical precision
  811 + // to the actual delta values, will will use these values for the
  812 + // interpolation (the m_quadratic_binning member will be set to true).
  813 + // Otherwise, the original delta values will be kept and the
  814 + // m_quadratic_binning member will be set to false.
  815 + int ndbins = m_deltas.size();
  816 + if (ndbins > 1) {
  817 + m_quadratic_binning = true;
  818 + double binsize = (std::sqrt(m_deltas[ndbins-1] * gammalib::deg2rad) -
  819 + std::sqrt(m_deltas[0] * gammalib::deg2rad)) /
  820 + double(ndbins-1);
  821 + for (int i = 0; i < ndbins; ++i) {
  822 + double delta = binsize * (double(i) + 0.5);
  823 + double diff = std::abs(delta*delta-m_deltas[i] * gammalib::deg2rad);
  824 + if (diff > 1.0e-6) {
  825 + m_quadratic_binning = false;
  826 + break;
  827 + }
  828 + m_deltas_cache.append(delta);
  829 + }
  830 + }
  831 +
  832 + // If we do not have square root binning then use the original values
  833 + // but convert them to radians
  834 + if (!m_quadratic_binning) {
  835 + m_deltas_cache.clear();
  836 + for (int i = 0; i < m_deltas.size(); ++i) {
  837 + m_deltas_cache.append(m_deltas[i] * gammalib::deg2rad);
  838 + }
  839 + }
774 840  
775 841 // Return
776 842 return;
... ...