Josh Cardenzana / gammalib

Commit 8b30ae91a7ef17e8989bfc41a869581112844776

Authored by Jürgen Knödlseder
1 parent 54700f26

Exclude upper edge boundary in GSkyMap::contains() method (#3012)

Inline Side-by-side
Showing 6 changed files with 41 additions and 18 deletions
ChangeLog
  View file @8b30ae9
... ... @@ -4,6 +4,7 @@
4 4 ========================
5 5  
6 6 Add GObservations::npred() methods (#3062)
  7 + Exclude upper edge boundary in GSkyMap::contains() method (#3012)
7 8 Add GSkyMap::extract() methods (#2990)
8 9 Add GCOMDri::cone_content() method (#2973)
9 10 Add GCOMObservation::drm(GModels&) method (#2973)
... ...
... ... @@ -115,6 +115,12 @@ None
115 115  
116 116 12. Sky module
117 117 --------------
  118 +Exclude the upper sky map boundary in GSkyMap::contains() method. Originally the
  119 +boundary was excluded, but this was changed in commit 0ecb6a07 when the
  120 +GSkyMap::flux() method was added since the flux method evaluates intensities at
  121 +the upper pixel boundaries. Now the GSkyMap::flux() method evaluates intensities
  122 +just below the upper pixel boundaries (#3012).
  123 +
118 124 Two GSkyMap::extract() methods were added. The first method extracts a patch of
119 125 pixels, defined by the start and stop pixel indices in X and Y direction. The second
120 126 method extracts all patch of pixels that is defined by inclusion sky regions (#2990).
... ...
doc/source/admin/release_history/1.7.rst
  View file @8b30ae9
... ... @@ -17,6 +17,8 @@ In particular, this release provides:
17 17 Bug fixes
18 18 ---------
19 19  
  20 +* [`3012 <https://cta-redmine.irap.omp.eu/issues/3012>`_] -
  21 + Exclude upper edge boundary in ``GSkyMap::contains()`` method.
20 22 * Correctly write model attributes in ``GCOMModelDRBFitting`` class.
21 23  
22 24  
... ...
inst/cta/test/test_CTA.cpp
  View file @8b30ae9
... ... @@ -705,12 +705,13 @@ void TestGCTAResponse::test_response_irf_diffuse(void)
705 705 // Set reference value
706 706 //double ref = 13803.800313356;
707 707 //const double ref = 13803.5186374;
708   - //const double ref = 13803.6932774; // after GWcs::solidangle improvement
709   - //const double ref = 13803.3453994; // after Irf computation improvements
710   - //const double ref = 13803.3461768; // using skymap flux() method
711   - //const double ref = 13803.5251338; // after correcting contains() method
712   - //const double ref = 13803.058889257; // after increasing integration accuracy
713   - const double ref = 13803.4194464338; // after increasing integration accuracy
  708 + //const double ref = 13803.6932774; // after GWcs::solidangle improvement
  709 + //const double ref = 13803.3453994; // after Irf computation improvements
  710 + //const double ref = 13803.3461768; // using skymap flux() method
  711 + //const double ref = 13803.5251338; // after correcting contains() method
  712 + //const double ref = 13803.058889257; // after increasing integration accuracy
  713 + //const double ref = 13803.4194464338; // after increasing integration accuracy
  714 + const double ref = 13803.1566304811; // after reverting contains() method
714 715  
715 716 // Set parameters
716 717 double src_ra = 201.3651;
... ... @@ -795,7 +796,8 @@ void TestGCTAResponse::test_response_npred_diffuse(void)
795 796 //const double ref = 11212.7161901; // npred_diffuse precision to iter=8
796 797 //const double ref = 11238.6928076; // npred_diffuse precision to iter=9
797 798 //const double ref = 11238.6933328; // using skymap flux() method
798   - const double ref = 11238.8390288; // after correcting contains() method
  799 + //const double ref = 11238.8390288; // after correcting contains() method
  800 + const double ref = 11238.6933449; // after reverting contains() method
799 801  
800 802 // Set parameters
801 803 double src_ra = 201.3651;
... ...
src/sky/GSkyMap.cpp
  View file @8b30ae9
... ... @@ -1567,14 +1567,18 @@ double GSkyMap::flux(const GSkyPixel&amp; pixel, const int&amp; map) const
1567 1567 // Get centre
1568 1568 GSkyDir centre = pix2dir(pixel);
1569 1569  
  1570 + // Set upper pixel edge just beyond the upper boundary since the
  1571 + // upper boundary is excluded
  1572 + const double up = 0.4999999;
  1573 +
1570 1574 // Get boundaries
1571 1575 GSkyDir boundary1 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()-0.5));
1572 1576 GSkyDir boundary2 = pix2dir(GSkyPixel(pixel.x(), pixel.y()-0.5));
1573   - GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()-0.5));
1574   - GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()));
1575   - GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+0.5, pixel.y()+0.5));
1576   - GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+0.5));
1577   - GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+0.5));
  1577 + GSkyDir boundary3 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()-0.5));
  1578 + GSkyDir boundary4 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()));
  1579 + GSkyDir boundary5 = pix2dir(GSkyPixel(pixel.x()+up, pixel.y()+up));
  1580 + GSkyDir boundary6 = pix2dir(GSkyPixel(pixel.x(), pixel.y()+up));
  1581 + GSkyDir boundary7 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()+up));
1578 1582 GSkyDir boundary8 = pix2dir(GSkyPixel(pixel.x()-0.5, pixel.y()));
1579 1583  
1580 1584 // Compute intensities
... ... @@ -1789,7 +1793,15 @@ bool GSkyMap::contains(const GSkyDir&amp; dir) const
1789 1793 * @return Trus if pixels is within map, false otherwise.
1790 1794 *
1791 1795 * Checks whether the specified sky map @p pixel falls within the skymap
1792   - * or not.
  1796 + * or not. A pixel is considered to fall in the sky may if the value is
  1797 + * contained in the interval
  1798 + *
  1799 + * \f[
  1800 + * [-0.5, n-0.5[
  1801 + * \f]
  1802 + *
  1803 + * where \f$n\f$ is the number of sky map pixels. Note that the upper bound
  1804 + * is excluded from the interval.
1793 1805 ***************************************************************************/
1794 1806 bool GSkyMap::contains(const GSkyPixel& pixel) const
1795 1807 {
... ... @@ -1799,7 +1811,7 @@ bool GSkyMap::contains(const GSkyPixel&amp; pixel) const
1799 1811 // Test 1D pixels
1800 1812 if (pixel.is_1D()) {
1801 1813 if (pixel.index() >= -0.5 &&
1802   - pixel.index() <= double(m_num_pixels)-0.5) {
  1814 + pixel.index() < double(m_num_pixels)-0.5) {
1803 1815 inmap = true;
1804 1816 }
1805 1817 }
... ... @@ -1808,8 +1820,8 @@ bool GSkyMap::contains(const GSkyPixel&amp; pixel) const
1808 1820 else if (pixel.is_2D()) {
1809 1821  
1810 1822 // If pixel is in range then set containment flag to true
1811   - if ((pixel.x() >= -0.5 && pixel.x() <= double(m_num_x)-0.5) &&
1812   - (pixel.y() >= -0.5 && pixel.y() <= double(m_num_y)-0.5)) {
  1823 + if ((pixel.x() >= -0.5 && pixel.x() < double(m_num_x)-0.5) &&
  1824 + (pixel.y() >= -0.5 && pixel.y() < double(m_num_y)-0.5)) {
1813 1825 inmap = true;
1814 1826 }
1815 1827  
... ...
test/test_GModel.cpp
  View file @8b30ae9
... ... @@ -834,14 +834,14 @@ void TestGModel::test_diffuse_map(void)
834 834 GModelSpatialDiffuseMap map_norm(m_map_file, 3.0, true);
835 835 test_assert(map_norm.normalize(),
836 836 "Check that model is normalized.");
837   - test_value(map_norm.eval(photon), 13069.7741687, 0.2,
  837 + test_value(map_norm.eval(photon), 13069.6047507, 0.2,
838 838 "Check that model has the correct skymap intensity.");
839 839 GXmlElement element;
840 840 map_norm.write(element);
841 841 map_norm.read(element);
842 842 test_assert(map_norm.normalize(),
843 843 "Check that model is normalized.");
844   - test_value(map_norm.eval(photon), 13069.7741687, 0.2,
  844 + test_value(map_norm.eval(photon), 13069.6047507, 0.2,
845 845 "Check that model has the correct skymap intensity.");
846 846  
847 847 // Test non-normalized map
... ...