Sergio Hernández Cadena / gammalib

Commit a96dc02e2bfc5673f6291b5b7d91285f0354e515

Authored by Jürgen Knödlseder
1 parent 07220862

Properly set region for GModelSpatialDiffuseCube class (#3145)

Inline Side-by-side
Showing 6 changed files with 94 additions and 89 deletions
ChangeLog
  View file @a96dc02
1   -2020-06-26
  1 +2020-06-30
2 2  
3 3 * Version 1.7.0 released
4 4 ========================
... ... @@ -9,6 +9,7 @@
9 9 Fix GModelSpatialDiffuseCube pickeling issue (#3179)
10 10 Rename GModelSpatialDiffuseCube::set_mc_cone() to mc_cone() (#3179)
11 11 Rename GModelSpatialDiffuseMap::set_mc_cone() to mc_cone() (#3179)
  12 + Properly set region for GModelSpatialDiffuseCube class (#3145)
12 13 Avoid rounding errors in GSkyMap operators (#3126)
13 14 Fix inverted min-max boundary in GOptimizePar for negative scale factors (#3072)
14 15 Add GObservations::npred() methods (#3062)
... ...
... ... @@ -122,6 +122,9 @@ The set_mc_cone() methods of the GModelSpatialDiffuseCube and GModelSpatialDiffu
122 122 class were renamed to mc_cone() and mc_cone() getter methods were implemented
123 123 for both classes. This solves a problem with pickling both classes (#3179).
124 124  
  125 +Set the GModelSpatialDiffuseCube::m_region member using the m_cube::region_circle()
  126 +method and implemented the GModelSpatialDiffuseCube::contains() method (#3145).
  127 +
125 128 The method GModelSpatial::flux() was added that extracts the flux by integrating
126 129 over a specific sky region. For the moment only circular sky regions are
127 130 implemented (#2961).
... ...
doc/source/admin/release_history/1.7.rst
  View file @a96dc02
... ... @@ -37,6 +37,8 @@ Improvements
37 37 * Add equality and non-equality operators to ``GSkyRegionCircle`` class.
38 38 * [`3202 <https://cta-redmine.irap.omp.eu/issues/3202>`_] -
39 39 Replace ``GCTAResponseCache`` class by ``GResponseCache`` class.
  40 +* [`3145 <https://cta-redmine.irap.omp.eu/issues/3145>`_] -
  41 + Properly set region for ``GModelSpatialDiffuseCube``class.
40 42 * [`3062 <https://cta-redmine.irap.omp.eu/issues/3062>`_] -
41 43 Add ``GObservations::npred()`` methods.
42 44 * [`2990 <https://cta-redmine.irap.omp.eu/issues/2990>`_] -
... ...
src/model/GModelSpatialDiffuseCube.cpp
  View file @a96dc02
... ... @@ -511,13 +511,32 @@ GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy&amp; energy,
511 511 *
512 512 * Signals whether a sky direction falls within the bounding circle of
513 513 * the diffuse cube.
514   - *
515   - * @todo To be implemented.
516 514 ***************************************************************************/
517 515 bool GModelSpatialDiffuseCube::contains(const GSkyDir& dir,
518 516 const double& margin) const
519 517 {
520   - return (true);
  518 + // Make sure that cube in online
  519 + fetch_cube();
  520 +
  521 + // Initialise containment flag
  522 + bool contains = false;
  523 +
  524 + // Continue only if radius is positive
  525 + if (m_region.radius() > 0.0) {
  526 +
  527 + // Compute distance to centre
  528 + double distance = m_region.centre().dist_deg(dir);
  529 +
  530 + // If distance is smaller than radius plus margin we consider
  531 + // the position to be contained within the bounding circle
  532 + if (distance < m_region.radius() + margin) {
  533 + contains = true;
  534 + }
  535 +
  536 + }
  537 +
  538 + // Return containment
  539 + return (contains);
521 540 }
522 541  
523 542  
... ... @@ -1224,8 +1243,9 @@ void GModelSpatialDiffuseCube::fetch_cube(void) const
1224 1243 ***************************************************************************/
1225 1244 void GModelSpatialDiffuseCube::load_cube(const GFilename& filename)
1226 1245 {
1227   - // Initialise skymap
  1246 + // Initialise sky map cube, region and log(energy) nodes
1228 1247 m_cube.clear();
  1248 + m_region.clear();
1229 1249 m_logE.clear();
1230 1250  
1231 1251 // Store filename of cube (for XML writing). Note that we do
... ... @@ -1264,6 +1284,12 @@ void GModelSpatialDiffuseCube::load_cube(const GFilename&amp; filename)
1264 1284 // Set energy boundaries
1265 1285 set_energy_boundaries();
1266 1286  
  1287 + // Get region circle
  1288 + m_region = m_cube.region_circle();
  1289 +
  1290 + // Set simulation cone
  1291 + mc_cone(m_region);
  1292 +
1267 1293 // Return
1268 1294 return;
1269 1295 }
... ... @@ -1381,16 +1407,11 @@ double GModelSpatialDiffuseCube::cube_intensity(const GPhoton&amp; photon) const
1381 1407  
1382 1408 /***********************************************************************//**
1383 1409 * @brief Set boundary sky region
1384   - *
1385   - * @todo Implement determination of the cube boundary circle
1386 1410 ***************************************************************************/
1387 1411 void GModelSpatialDiffuseCube::set_region(void) const
1388 1412 {
1389   - // Set sky region circle (all sky)
1390   - GSkyRegionCircle region(0.0, 0.0, 180.0);
1391   -
1392   - // Set region (circumvent const correctness)
1393   - const_cast<GModelSpatialDiffuseCube*>(this)->m_region = region;
  1413 + // Make sure that cube in online
  1414 + fetch_cube();
1394 1415  
1395 1416 // Return
1396 1417 return;
... ...
src/model/GModelSpatialDiffuseMap.cpp
  View file @a96dc02
... ... @@ -923,9 +923,8 @@ void GModelSpatialDiffuseMap::free_members(void)
923 923 ***************************************************************************/
924 924 void GModelSpatialDiffuseMap::prepare_map(void)
925 925 {
926   - // Initialise centre and radius
  926 + // Initialise region
927 927 m_region.clear();
928   - //m_radius = 0.0;
929 928  
930 929 // Determine number of skymap pixels
931 930 int npix = m_map.npix();
... ... @@ -958,33 +957,8 @@ void GModelSpatialDiffuseMap::prepare_map(void)
958 957 }
959 958 }
960 959  
961   - // If we have a HealPix map then set radius to 180 deg
962   - if (m_map.projection()->code() == "HPX") {
963   - m_region = GSkyRegionCircle(0.0, 0.0, 180.0);
964   - }
965   -
966   - // ... otherwise compute map centre and radius
967   - else {
968   -
969   - // Initialise maximum radius
970   - double max_radius = 0.0;
971   -
972   - // Get map centre
973   - GSkyPixel pixel(m_map.nx()/2.0, m_map.ny()/2.0);
974   - GSkyDir centre = m_map.pix2dir(pixel);
975   -
976   - // Determine map radius
977   - for (int i = 0; i < npix; ++i) {
978   - double radius = m_map.inx2dir(i).dist_deg(centre);
979   - if (radius > max_radius) {
980   - max_radius = radius;
981   - }
982   - }
983   -
984   - // Set sky region
985   - m_region = GSkyRegionCircle(centre, max_radius);
986   -
987   - } // endelse: computed map centre and radius
  960 + // Get region circle
  961 + m_region = m_map.region_circle();
988 962  
989 963 // Set simulation cone
990 964 mc_cone(m_region);
... ... @@ -1001,14 +975,6 @@ void GModelSpatialDiffuseMap::prepare_map(void)
1001 975 ***************************************************************************/
1002 976 void GModelSpatialDiffuseMap::set_region(void) const
1003 977 {
1004   - /*
1005   - // Set sky region circle (all sky)
1006   - GSkyRegionCircle region(m_centre, m_radius);
1007   -
1008   - // Set region (circumvent const correctness)
1009   - const_cast<GModelSpatialDiffuseMap*>(this)->m_region = region;
1010   - */
1011   -
1012 978 // Return
1013 979 return;
1014 980 }
... ...
test/test_GModel.cpp
  View file @a96dc02
... ... @@ -261,7 +261,7 @@ void TestGModel::test_model_par(void)
261 261 test_assert(!par1.has_min(), "Parameter shall have no minimum.");
262 262 test_assert(!par1.has_max(), "Parameter shall have no maximum.");
263 263 test_assert(!par1.has_range(), "Parameter shall have no range.");
264   -
  264 +
265 265 // Test parameter value constructor (value version)
266 266 GModelPar par2("Test parameter", 47.0);
267 267 test_value(par2.value(), 47.0);
... ... @@ -272,7 +272,7 @@ void TestGModel::test_model_par(void)
272 272 test_assert(!par2.has_min(), "Parameter shall have no minimum.");
273 273 test_assert(!par2.has_max(), "Parameter shall have no maximum.");
274 274 test_assert(!par2.has_range(), "Parameter shall have no range.");
275   -
  275 +
276 276 // Test invalid parameter constructor (factor & scale version)
277 277 test_try("Test invalid parameter constructor");
278 278 try {
... ... @@ -540,7 +540,7 @@ void TestGModel::test_sky_model(void)
540 540 sky4.spectral(&spec_plaw);
541 541 test_value(sky4.size(), 5); // 1 spatial, 3 spectral, 1 temporal
542 542 test_value(sky4.spectral()->classname(), "GModelSpectralPlaw");
543   -
  543 +
544 544 // Test XML constructor
545 545 GXml xml1(m_xml_file);
546 546 GXmlElement* element1 = xml1.element(0)->element(0);
... ... @@ -622,7 +622,7 @@ void TestGModel::test_point_source(void)
622 622 GModelSpatialPointSource model3(83.6331, +22.0145);
623 623 test_value(model3.ra(), 83.6331);
624 624 test_value(model3.dec(), +22.0145);
625   -
  625 +
626 626 // Test XML constructor
627 627 GXml xml(m_xml_model_point_plaw);
628 628 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -701,7 +701,7 @@ void TestGModel::test_diffuse_const(void)
701 701 // Test value constructor
702 702 GModelSpatialDiffuseConst model2(3.0);
703 703 test_value(model2.value(), 3.0);
704   -
  704 +
705 705 // Test XML constructor
706 706 GXml xml(m_xml_model_diffuse_const);
707 707 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -759,7 +759,7 @@ void TestGModel::test_diffuse_cube(void)
759 759 test_value(model3.value(), 3.0);
760 760 test_value(model3.filename().url(), "");
761 761 //test_assert(model.cube() == map, "Map cube is not the expected one");
762   -
  762 +
763 763 // Test XML constructor
764 764 GXml xml(m_xml_model_diffuse_cube);
765 765 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -769,20 +769,28 @@ void TestGModel::test_diffuse_cube(void)
769 769 test_value(model4.value(), 1.0);
770 770 test_value(model4.filename().url(), m_cube_file);
771 771  
772   - // Test value method
  772 + // Test value() method
773 773 model4.value(3.9);
774 774 test_value(model4.value(), 3.9);
775 775  
776   - // Test region method
  776 + // Test region() method
777 777 GSkyDir dir1;
778   - dir1.radec_deg(3.9, 3.9);
  778 + dir1.radec_deg(84.2, 22.0);
779 779 test_assert(model4.region()->contains(dir1), "Test region() method (inside)");
  780 + dir1.radec_deg(80.0, 24.0);
  781 + test_assert(!model4.region()->contains(dir1), "Test region() method (outside)");
780 782  
781   - // Test filename method
  783 + // Test contains() method
  784 + dir1.radec_deg(84.2, 22.0);
  785 + test_assert(model4.contains(dir1), "Test contains() method (inside)");
  786 + dir1.radec_deg(80.0, 24.0);
  787 + test_assert(!model4.contains(dir1), "Test contains() method (outside)");
  788 +
  789 + // Test filename() method
782 790 model4.filename("Help me!");
783 791 test_value(model4.filename().url(), "Help me!");
784 792  
785   - // Test cube method
  793 + // Test cube() method
786 794 model4.cube(GSkyMap("GAL", 16, "RING", 10));
787 795 test_value(model4.cube().npix(), 3072);
788 796  
... ... @@ -822,7 +830,7 @@ void TestGModel::test_diffuse_map(void)
822 830 GModelSpatialDiffuseMap model3(map, 3.0);
823 831 test_value(model3.value(), 3.0);
824 832 test_value(model3.filename().url(), "");
825   -
  833 +
826 834 // Test XML constructor and attribute methods
827 835 GXml xml(m_xml_model_diffuse_map);
828 836 GXmlElement* elementp = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -832,22 +840,28 @@ void TestGModel::test_diffuse_map(void)
832 840 test_value(model.value(), 1.0);
833 841 test_value(model.filename().url(), m_map_file);
834 842  
835   - // Test value method
  843 + // Test value() method
836 844 model.value(3.9);
837 845 test_value(model.value(), 3.9);
838 846  
839   - // Test region method
  847 + // Test region() method
840 848 GSkyDir dir1;
841 849 dir1.radec_deg(201.365, -43.019);
842 850 test_assert(model.region()->contains(dir1), "Test region() method (inside)");
843 851 dir1.radec_deg(180.0, -43.019);
844 852 test_assert(!model.region()->contains(dir1), "Test region() method (outside)");
845 853  
846   - // Test load method
  854 + // Test contains() method
  855 + dir1.radec_deg(201.365, -43.019);
  856 + test_assert(model.contains(dir1), "Test region() method (inside)");
  857 + dir1.radec_deg(180.0, -43.019);
  858 + test_assert(!model.contains(dir1), "Test region() method (outside)");
  859 +
  860 + // Test load() method
847 861 model.load(m_map_file);
848 862 test_value(model.filename().url(), m_map_file);
849 863  
850   - // Test map method
  864 + // Test map() method
851 865 model.map(GSkyMap("GAL", 16, "RING", 10));
852 866 test_value(model.map().npix(), 3072);
853 867  
... ... @@ -916,7 +930,7 @@ void TestGModel::test_radial_disk(void)
916 930 test_value(model2.ra(), 83.6331);
917 931 test_value(model2.dec(), 22.0145);
918 932 test_value(model2.radius(), 3.0);
919   -
  933 +
920 934 // Test XML constructor
921 935 GXml xml(m_xml_model_radial_disk);
922 936 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -1067,7 +1081,7 @@ void TestGModel::test_radial_gauss(void)
1067 1081 test_value(model2.ra(), 83.6331);
1068 1082 test_value(model2.dec(), 22.0145);
1069 1083 test_value(model2.sigma(), 3.0);
1070   -
  1084 +
1071 1085 // Test XML constructor
1072 1086 GXml xml(m_xml_model_radial_gauss);
1073 1087 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -1137,7 +1151,7 @@ void TestGModel::test_radial_shell(void)
1137 1151 test_value(model2.dec(), 22.0145);
1138 1152 test_value(model2.radius(), 3.0);
1139 1153 test_value(model2.width(), 1.0);
1140   -
  1154 +
1141 1155 // Test XML constructor
1142 1156 GXml xml(m_xml_model_radial_shell);
1143 1157 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -1213,7 +1227,7 @@ void TestGModel::test_elliptical_disk(void)
1213 1227 test_value(model2.posangle(), 45.0);
1214 1228 test_value(model2.semimajor(), 3.0);
1215 1229 test_value(model2.semiminor(), 2.0);
1216   -
  1230 +
1217 1231 // Test XML constructor
1218 1232 GXml xml(m_xml_model_elliptical_disk);
1219 1233 GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
... ... @@ -1389,7 +1403,7 @@ void TestGModel::test_spatial_composite(void)
1389 1403 GSkyDir dir2;
1390 1404 dir2.radec_deg(83.6331, 22.01);
1391 1405 test_assert(model2.region()->contains(dir2), "Test region() method (inside)");
1392   -
  1406 +
1393 1407 // Test access of individual parameters
1394 1408 test_value(model2["2:RA"].value(), 83.6331);
1395 1409 test_value(model2["2:DEC"].value(), 22.0145);
... ... @@ -1420,7 +1434,7 @@ void TestGModel::test_const(void)
1420 1434 // Test value constructor
1421 1435 GModelSpectralConst model2(3.0);
1422 1436 test_value(model2.value(), 3.0);
1423   -
  1437 +
1424 1438 // Test XML constructor
1425 1439 GXml xml(m_xml_model_point_const);
1426 1440 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1520,7 +1534,7 @@ void TestGModel::test_plaw(void)
1520 1534 test_value(model2.prefactor(), 2.0);
1521 1535 test_value(model2.index(), -2.1);
1522 1536 test_value(model2.pivot().MeV(), 100.0);
1523   -
  1537 +
1524 1538 // Test XML constructor
1525 1539 GXml xml(m_xml_model_point_plaw);
1526 1540 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1641,7 +1655,7 @@ void TestGModel::test_plaw_eflux(void)
1641 1655 test_value(model2.index(), -2.1);
1642 1656 test_value(model2.emin().MeV(), 10.0);
1643 1657 test_value(model2.emax().MeV(), 100.0);
1644   -
  1658 +
1645 1659 // Test XML constructor and value
1646 1660 GXml xml(m_xml_model_point_plaw_eflux);
1647 1661 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1702,7 +1716,7 @@ void TestGModel::test_eplaw(void)
1702 1716 test_value(model2.index(), -2.1);
1703 1717 test_value(model2.pivot().MeV(), 100.0);
1704 1718 test_value(model2.cutoff().GeV(), 1.0);
1705   -
  1719 +
1706 1720 // Test XML constructor
1707 1721 GXml xml(m_xml_model_point_eplaw);
1708 1722 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1841,7 +1855,7 @@ void TestGModel::test_supeplaw(void)
1841 1855 test_value(model2.pivot().MeV(), 100.0);
1842 1856 test_value(model2.cutoff().GeV(), 1.0);
1843 1857 test_value(model2.index2(), 1.1);
1844   -
  1858 +
1845 1859 // Test XML constructor
1846 1860 GXml xml(m_xml_model_point_supeplaw);
1847 1861 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1908,7 +1922,7 @@ void TestGModel::test_bplaw(void)
1908 1922 test_value(model2.index1(), -2.1);
1909 1923 test_value(model2.breakenergy().MeV(), 100.0);
1910 1924 test_value(model2.index2(), -2.8);
1911   -
  1925 +
1912 1926 // Test XML constructor
1913 1927 GXml xml(m_xml_model_point_bplaw);
1914 1928 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1962,7 +1976,7 @@ void TestGModel::test_smoothbplaw(void)
1962 1976 // Test void constructor
1963 1977 GModelSpectralSmoothBrokenPlaw model1;
1964 1978 test_value(model1.type(), "SmoothBrokenPowerLaw", "Check void model type");
1965   -
  1979 +
1966 1980 // Test value constructor
1967 1981 GModelSpectralSmoothBrokenPlaw model2(2.0, -2.1, GEnergy(100.0, "MeV"), -2.8,
1968 1982 GEnergy(1000.0, "MeV"), 2.0);
... ... @@ -1972,8 +1986,7 @@ void TestGModel::test_smoothbplaw(void)
1972 1986 test_value(model2.index2(), -2.8);
1973 1987 test_value(model2.breakenergy().MeV(), 1000.0);
1974 1988 test_value(model2.beta(), 2.0);
1975   -
1976   -
  1989 +
1977 1990 // Test XML constructor
1978 1991 GXml xml(m_xml_model_point_smoothbplaw);
1979 1992 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -1985,20 +1998,19 @@ void TestGModel::test_smoothbplaw(void)
1985 1998 test_value(model3.pivot().TeV(), 1.0);
1986 1999 test_value(model3.index2(), -2.70);
1987 2000 test_value(model3.breakenergy().TeV(), 0.3);
1988   -
1989   -
  2001 +
1990 2002 // Test prefactor method
1991 2003 model3.prefactor(2.3e-16);
1992 2004 test_value(model3.prefactor(), 2.3e-16);
1993   -
  2005 +
1994 2006 // Test index1 method
1995 2007 model3.index1(-2.6);
1996 2008 test_value(model3.index1(), -2.6);
1997   -
  2009 +
1998 2010 // Test pivot method
1999 2011 model3.pivot(GEnergy(0.5, "TeV"));
2000 2012 test_value(model3.pivot().TeV(), 0.5);
2001   -
  2013 +
2002 2014 // Test index2 method
2003 2015 model3.index2(-3.6);
2004 2016 test_value(model3.index2(), -3.6);
... ... @@ -2006,7 +2018,7 @@ void TestGModel::test_smoothbplaw(void)
2006 2018 // Test breakenergy method
2007 2019 model3.breakenergy(GEnergy(0.5, "TeV"));
2008 2020 test_value(model3.breakenergy().TeV(), 0.5);
2009   -
  2021 +
2010 2022 // Test operator access
2011 2023 const char* strarray[] = {"Prefactor", "Index1", "PivotEnergy", "Index2",
2012 2024 "BreakEnergy", "BreakSmoothness"};
... ... @@ -2020,7 +2032,7 @@ void TestGModel::test_smoothbplaw(void)
2020 2032 test_value(model3[keyname].error(), 1.9);
2021 2033 test_value(model3[keyname].gradient(), 0.8);
2022 2034 }
2023   -
  2035 +
2024 2036 // Exit test
2025 2037 return;
2026 2038 }
... ... @@ -2041,7 +2053,7 @@ void TestGModel::test_logparabola(void)
2041 2053 test_value(model2.index(), -2.1);
2042 2054 test_value(model2.pivot().MeV(), 100.0);
2043 2055 test_value(model2.curvature(), -0.2);
2044   -
  2056 +
2045 2057 // Test XML constructor
2046 2058 GXml xml(m_xml_model_point_logparabola);
2047 2059 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -2286,7 +2298,7 @@ void TestGModel::test_filefct(void)
2286 2298 test_value(model2.filename().url(), m_filefct,
2287 2299 "Check file function data file name");
2288 2300 test_value(model2.norm(), 2.0);
2289   -
  2301 +
2290 2302 // Test XML constructor
2291 2303 GXml xml(m_xml_model_point_filefct);
2292 2304 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ... @@ -2338,7 +2350,7 @@ void TestGModel::test_table(void)
2338 2350 test_value(model2.filename().url(), m_table,
2339 2351 "Check file function data file name");
2340 2352 test_value(model2.norm(), 2.0);
2341   -
  2353 +
2342 2354 // Test XML constructor
2343 2355 GXml xml(m_xml_model_point_table);
2344 2356 GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
... ...