Leonardo Di Venere / gammalib

Commit f02bf3f306d209454f202b4579c5073388bb0823

Authored by Jürgen Knödlseder
1 parent dabe9da8

Revise GCTAEdispCube::operator() interface (#2721)

Inline Side-by-side
Showing 9 changed files with 148 additions and 88 deletions
ChangeLog
  View file @f02bf3f
1   -2018-11-06
  1 +2018-11-08
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Revise GCTAEdispCube::operator() interface (#2721)
6 7 Rename GCTAResponseIrf::lo_save_thres() and GCTAResponseIrf::hi_save_thres()
7 8 to GCTAResponseIrf::lo_safe_thres() and GCTAResponseIrf::hi_safe_thres()
8 9 (#2716)
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -6 November 2018
  3 +8 November 2018
4 4  
5 5  
6 6 1. Introduction
... ... @@ -151,12 +151,12 @@ The arguments for the following methods have been changed:
151 151 - GCTAModelRadial::mc(GRan&)
152 152 - GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&)
153 153 - GCTAEdisp::mc(GRan&, GEnergy&, double&, double&, double&, double&)
  154 +- GCTACubeEdisp::operator()(GEnergy&, GEnergy&, GSkyDir&)
154 155 - GCTAResponseIrf::edisp(GEnergy&, GEnergy&, double&, double&, double&, double&)
155 156 - GLATAeff(GFits&, std::string&)
156 157 - GLATEdisp(GFits&, std::string&)
157 158 - GLATPsf(GFits&, std::string&)
158 159  
159   -
160 160 The return value of the following methods has been changed:
161 161 - none
162 162  
... ... @@ -314,6 +314,10 @@ None
314 314  
315 315 20. CTA interface
316 316 -----------------
  317 +The arguments of the GCTACubeEdisp::operator() were changed from, passing
  318 +now reconstructed and true energy objects, followed by a true sky direction.
  319 +This is more consistent with the GCTAEdisp interface (#2721).
  320 +
317 321 Rename GCTAResponseIrf::lo_save_thres() and GCTAResponseIrf::hi_save_thres()
318 322 to GCTAResponseIrf::lo_safe_thres() and GCTAResponseIrf::hi_safe_thres()
319 323 (#2716).
... ...
doc/source/admin/release_history/1.6.rst
  View file @f02bf3f
... ... @@ -47,6 +47,8 @@ Bug fixes
47 47 Improvements
48 48 ------------
49 49  
  50 +* [`2721 <https://cta-redmine.irap.omp.eu/issues/2721>`_] -
  51 + Revise ``GCTAEdispCube::operator()`` interface
50 52 * [`2715 <https://cta-redmine.irap.omp.eu/issues/2715>`_] -
51 53 Add generic CTA response cache
52 54 * [`2704 <https://cta-redmine.irap.omp.eu/issues/2704>`_] -
... ...
inst/cta/include/GCTACubeEdisp.hpp
  View file @f02bf3f
1 1 /***************************************************************************
2 2 * GCTACubeEdisp.hpp - CTA cube analysis energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2017 by Michael Mayer *
  4 + * copyright (C) 2016-2018 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -84,9 +84,9 @@ public:
84 84  
85 85 // Operators
86 86 GCTACubeEdisp& operator=(const GCTACubeEdisp& cube);
87   - double operator()(const GSkyDir& dir,
88   - const double& migra,
89   - const GEnergy& energy) const;
  87 + double operator()(const GEnergy& ereco,
  88 + const GEnergy& etrue,
  89 + const GSkyDir& dir) const;
90 90  
91 91 // Methods
92 92 void clear(void);
... ... @@ -116,7 +116,7 @@ protected:
116 116 void clear_cube(void);
117 117 void fill_cube(const GCTAObservation& obs, GSkyMap* exposure = NULL,
118 118 GLog* log = NULL);
119   - void update(const double& migra, const double& logEsrc) const;
  119 + void update(const GEnergy& ereco, const GEnergy& etrue) const;
120 120 void set_eng_axis(void);
121 121 void compute_ebounds(void) const;
122 122  
... ...
inst/cta/pyext/GCTACubeEdisp.i
  View file @f02bf3f
... ... @@ -64,9 +64,9 @@ public:
64 64 virtual ~GCTACubeEdisp(void);
65 65  
66 66 // Interpolation operator
67   - double operator()(const GSkyDir& dir,
68   - const double& migra,
69   - const GEnergy& energy) const;
  67 + double operator()(const GEnergy& ereco,
  68 + const GEnergy& etrue,
  69 + const GSkyDir& dir) const;
70 70  
71 71 // Methods
72 72 void clear(void);
... ...
inst/cta/src/GCTACubeEdisp.cpp
  View file @f02bf3f
... ... @@ -323,23 +323,23 @@ GCTACubeEdisp&amp; GCTACubeEdisp::operator=(const GCTACubeEdisp&amp; cube)
323 323  
324 324  
325 325 /***********************************************************************//**
326   - * @brief Return energy dispersion (in units of MeV^-1)
  326 + * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
327 327 *
  328 + * @param[in] ereco Reconstructed event energy.
  329 + * @param[in] etrue True photon energy.
328 330 * @param[in] dir Coordinate of the true photon position.
329   - * @param[in] migra Energy migration (reconstructed over true energy).
330   - * @param[in] srcEng True photon energy.
331   - * @return Energy dispersion (in units of MeV^-1)
  331 + * @return Energy dispersion (MeV\f$^{-1}\f$)
332 332 *
333 333 * Returns the energy dispersion for a given energy migration (reconstructed
334 334 * over true energy), true photon energy, and sky direction in units of
335   - * MeV^-1.
  335 + * MeV\f$^{-1}\f$.
336 336 ***************************************************************************/
337   -double GCTACubeEdisp::operator()(const GSkyDir& dir,
338   - const double& migra,
339   - const GEnergy& srcEng) const
  337 +double GCTACubeEdisp::operator()(const GEnergy& ereco,
  338 + const GEnergy& etrue,
  339 + const GSkyDir& dir) const
340 340 {
341 341 // Update indices and weighting factors for interpolation
342   - update(migra, srcEng.log10TeV());
  342 + update(ereco, etrue);
343 343  
344 344 // Perform bi-linear interpolation
345 345 double edisp = m_wgt1 * m_cube(dir, m_inx1) +
... ... @@ -837,12 +837,21 @@ void GCTACubeEdisp::clear_cube(void)
837 837 /***********************************************************************//**
838 838 * @brief Fill energy dispersion cube from observation container
839 839 *
840   - * @param[in] obs Observation.
  840 + * @param[in] obs CTA observation.
841 841 * @param[in] exposure Pointer towards exposure map.
842 842 * @param[in] log Pointer towards logger.
843 843 *
844 844 * @exception GException::invalid_value
845 845 * No RoI or response found in CTA observation.
  846 + *
  847 + * Fills the energy dispersion cube from the observation container.
  848 + *
  849 + * If the @p exposure argument is not NULL, the energy dispersion is weighted
  850 + * by the exposure, and the exposure weighting is added to the exposure map.
  851 + *
  852 + * If the @p log argument is not NULL, the method puts information about
  853 + * exclusion and inclusion of a CTA observation into the energy dispersion
  854 + * cube.
846 855 ***************************************************************************/
847 856 void GCTACubeEdisp::fill_cube(const GCTAObservation& obs,
848 857 GSkyMap* exposure,
... ... @@ -965,12 +974,14 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
965 974 // Get log10 of true energy in TeV
966 975 double logEsrc = m_energies[iebin].log10TeV();
967 976  
968   - // Compute exposure weight
969   - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
970   - obs.livetime();
971   -
972   - // If available, accumulate weights
  977 + // Compute exposure weight. If no exposure map is
  978 + // specified, the weight is one. Otherwise, the weight
  979 + // is the effective area for this offset angle and
  980 + // source energy times the livetime of the observation.
  981 + double weight = 1.0;
973 982 if (exposure != NULL) {
  983 + weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
  984 + obs.livetime();
974 985 (*exposure)(pixel, iebin) += weight;
975 986 }
976 987  
... ... @@ -985,9 +996,8 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
985 996 continue;
986 997 }
987 998  
988   - // Compute log10 of reconstructed energy in TeV
989   - double logEobs = std::log10(migra) + logEsrc;
990   - GEnergy Eobs(std::pow(10.0, logEobs), "TeV");
  999 + // Compute reconstructed energy
  1000 + GEnergy Eobs = migra * m_energies[iebin];
991 1001  
992 1002 // Set map index
993 1003 int imap = offset(imigra, iebin);
... ... @@ -1016,19 +1026,22 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
1016 1026 /***********************************************************************//**
1017 1027 * @brief Update energy dispersion cube indices and weights cache
1018 1028 *
1019   - * @param[in] migra Reconstructed over true energy.
1020   - * @param[in] logEsrc Log10 of true photon energy in TeV.
  1029 + * @param[in] ereco Reconstructed event energy.
  1030 + * @param[in] etrue True photon energy.
1021 1031 *
1022 1032 * Updates the energy dispersion cube indices, stored in the members m_inx1,
1023 1033 * m_inx2, m_inx3, and m_inx4, and weights cache, stored in m_wgt1, m_wgt2,
1024 1034 * m_wgt3, and m_wgt4.
1025 1035 ***************************************************************************/
1026   -void GCTACubeEdisp::update(const double& migra, const double& logEsrc) const
  1036 +void GCTACubeEdisp::update(const GEnergy& ereco, const GEnergy& etrue) const
1027 1037 {
1028 1038 // Set node array for energy interpolation
1029   - m_elogmeans.set_value(logEsrc);
  1039 + m_elogmeans.set_value(etrue.log10TeV());
  1040 +
  1041 + // Compute migration value
  1042 + double migra = ereco/etrue;
1030 1043  
1031   - // Set node array for delta interpolation
  1044 + // Set node array for migra interpolation
1032 1045 m_migras.set_value(migra);
1033 1046  
1034 1047 // Set indices for bi-linear interpolation
... ...
inst/cta/src/GCTAResponseCube.cpp
  View file @f02bf3f
... ... @@ -354,14 +354,8 @@ double GCTAResponseCube::irf(const GEvent&amp; event,
354 354  
355 355 // Multiply-in energy dispersion
356 356 if (use_edisp() && irf > 0.0) {
357   -
358   - // Determine fraction of reconstructed and true energy
359   - double migra = obsEng.TeV() / srcEng.TeV();
360   -
361   - // Multiply-in energy dispersion
362   - irf *= edisp()(srcDir, migra, srcEng);
363   -
364   - } // endif: energy dispersion was available and psf was non-zero
  357 + irf *= edisp()(obsEng, srcEng, srcDir);
  358 + }
365 359  
366 360 // Divide by livetime
367 361 irf /= livetime;
... ... @@ -1135,9 +1129,6 @@ double GCTAResponseCube::irf_ptsrc(const GEvent&amp; event,
1135 1129 // direction in radians
1136 1130 double delta = bin->dir().dir().dist(srcDir);
1137 1131  
1138   - // Determine fraction of reconstructed and true energy
1139   - double migra = bin->energy().TeV() / srcEng.TeV();
1140   -
1141 1132 // Get maximum angular separation for PSF (in radians)
1142 1133 double delta_max = psf().delta_max();
1143 1134  
... ... @@ -1162,11 +1153,8 @@ double GCTAResponseCube::irf_ptsrc(const GEvent&amp; event,
1162 1153  
1163 1154 // Multiply-in energy dispersion
1164 1155 if (use_edisp() && irf > 0.0) {
1165   -
1166   - // Multiply-in energy dispersion
1167   - irf *= edisp()(srcDir, migra, srcEng);
1168   -
1169   - } // endif: energy dispersion was available and psf was non-zero
  1156 + irf *= edisp()(bin->energy(), srcEng, srcDir);
  1157 + }
1170 1158  
1171 1159 // Apply deadtime correction
1172 1160 irf *= exposure().deadc();
... ... @@ -1236,9 +1224,6 @@ double GCTAResponseCube::irf_radial(const GEvent&amp; event,
1236 1224 // (radians)
1237 1225 double rho_obs = model->dir().dist(obsDir);
1238 1226  
1239   - // Determine fraction of reconstructed and true energy
1240   - double migra = bin->energy().TeV() / srcEng.TeV();
1241   -
1242 1227 // Get livetime (in seconds)
1243 1228 double livetime = exposure().livetime();
1244 1229  
... ... @@ -1264,17 +1249,14 @@ double GCTAResponseCube::irf_radial(const GEvent&amp; event,
1264 1249 irf *= psf_radial(model, rho_obs, obsDir, srcEng, obsTime);
1265 1250  
1266 1251 // Multiply-in energy dispersion
  1252 + //
  1253 + // The current code assumes that the energy dispersion at the
  1254 + // observed and true event direction does not vary
  1255 + // significantly. In other words, the code assumes that the
  1256 + // energy dispersion is constant over the size of the PSF.
1267 1257 if (use_edisp() && irf > 0.0) {
1268   -
1269   - // Multiply-in energy dispersion
1270   - //
1271   - // The current code assumes that the energy dispersion at the
1272   - // observed and true event direction does not vary
1273   - // significantly. In other words, the code assumes that the
1274   - // energy dispersion is constant over the size of the PSF.
1275   - irf *= edisp()(obsDir, migra, srcEng);
1276   -
1277   - } // endif: energy dispersion was available and psf was non-zero
  1258 + irf *= edisp()(bin->energy(), srcEng, obsDir);
  1259 + }
1278 1260  
1279 1261 // Apply deadtime correction
1280 1262 irf *= exposure().deadc();
... ... @@ -1344,9 +1326,6 @@ double GCTAResponseCube::irf_elliptical(const GEvent&amp; event,
1344 1326 double rho_obs = model->dir().dist(obsDir);
1345 1327 double posangle_obs = model->dir().posang(obsDir);
1346 1328  
1347   - // Determine fraction of reconstructed and true energy
1348   - double migra = bin->energy().TeV() / srcEng.TeV();
1349   -
1350 1329 // Get livetime (in seconds)
1351 1330 double livetime = exposure().livetime();
1352 1331  
... ... @@ -1372,23 +1351,20 @@ double GCTAResponseCube::irf_elliptical(const GEvent&amp; event,
1372 1351 irf *= psf_elliptical(model, rho_obs, posangle_obs, obsDir, srcEng, obsTime);
1373 1352  
1374 1353 // Multiply-in energy dispersion
  1354 + //
  1355 + // The current code assumes that the energy dispersion at the
  1356 + // observed and true event direction does not vary
  1357 + // significantly. In other words, the code assumes that the
  1358 + // energy dispersion is constant over the size of the PSF.
1375 1359 if (use_edisp() && irf > 0.0) {
1376   -
1377   - // Multiply-in energy dispersion
1378   - //
1379   - // The current code assumes that the energy dispersion at the
1380   - // observed and true event direction does not vary
1381   - // significantly. In other words, the code assumes that the
1382   - // energy dispersion is constant over the size of the PSF.
1383   - irf *= edisp()(obsDir, migra, srcEng);
1384   -
1385   - } // endif: energy dispersion was available and psf was non-zero
  1360 + irf *= edisp()(bin->energy(), srcEng, obsDir);
  1361 + }
1386 1362  
1387 1363 // Apply deadtime correction
1388 1364 irf *= exposure().deadc();
1389 1365  
1390 1366 } // endif: exposure was positive
1391   -
  1367 +
1392 1368 } // endif: we were sufficiently close and livetime >0
1393 1369  
1394 1370 // Compile option: Check for NaN/Inf
... ... @@ -1471,8 +1447,7 @@ double GCTAResponseCube::irf_diffuse(const GEvent&amp; event,
1471 1447 if (use_edisp()) {
1472 1448 irf = cache->irf(bin->ipix(), source.energy());
1473 1449 if (irf > 0.0) {
1474   - irf *= edisp()(bin->dir().dir(), bin->energy()/source.energy(),
1475   - source.energy());
  1450 + irf *= edisp()(bin->energy(), source.energy(), bin->dir().dir());
1476 1451 }
1477 1452 }
1478 1453  
... ...
inst/cta/test/test_CTA.cpp
  View file @f02bf3f
... ... @@ -1364,8 +1364,8 @@ void TestGCTAResponse::test_response_bkgcube(void)
1364 1364 ***************************************************************************/
1365 1365 void TestGCTAResponse::test_response_edispcube(void)
1366 1366 {
1367   - // Test PSF cube constructors
1368   - test_try("CTA Edisp cube void constructor");
  1367 + // Test energy dispersion cube constructors
  1368 + test_try("CTA energy dispersion cube void constructor");
1369 1369 try {
1370 1370 GCTACubeEdisp cube;
1371 1371 test_try_success();
... ... @@ -1373,7 +1373,7 @@ void TestGCTAResponse::test_response_edispcube(void)
1373 1373 catch (std::exception &e) {
1374 1374 test_try_failure(e);
1375 1375 }
1376   - test_try("CTA Edisp cube map constructor");
  1376 + test_try("CTA energy dispersion cube map constructor");
1377 1377 try {
1378 1378 GEnergies energies(20, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV"));
1379 1379 GCTACubeEdisp cube("CAR", "CEL", 83.63, 22.01, 0.4, 0.4, 10, 10,
... ... @@ -1384,24 +1384,41 @@ void TestGCTAResponse::test_response_edispcube(void)
1384 1384 test_try_failure(e);
1385 1385 }
1386 1386  
1387   - // Test set method
  1387 + // Load IRF for energy dispersion computation
1388 1388 GCTAObservation obs_cta;
1389 1389 obs_cta.load(cta_events);
1390 1390 obs_cta.response(cta_irf, GCaldb(cta_caldb));
1391   - GEnergies energies(20, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV"));
  1391 +
  1392 + // Initialise energy dispersion cube
  1393 + GEnergies energies(100, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV"));
1392 1394 GCTACubeEdisp cube("CAR", "CEL", 83.63, 22.01, 0.4, 0.4, 10, 10,
1393   - energies, 2.0, 20);
  1395 + energies, 3.0, 100);
  1396 +
  1397 + // Test setting of observation
1394 1398 cube.set(obs_cta);
1395   - cube.save("test_cta_edispcube_one.fits", true);
1396 1399  
1397   - // Test fill method
  1400 + // Test energy dispersion integration
  1401 + test_edispcube_integration(cube, 0.1, 10.0);
  1402 +
  1403 + // Test filling of single observation
1398 1404 GObservations obs;
1399 1405 obs_cta.id("000001");
1400 1406 obs.append(obs_cta);
  1407 + cube.fill(obs);
  1408 +
  1409 + // Test energy dispersion integration
  1410 + test_edispcube_integration(cube, 0.1, 10.0);
  1411 +
  1412 + // Test fill of two observations
1401 1413 obs_cta.id("000002");
1402 1414 obs.append(obs_cta);
1403 1415 cube.fill(obs);
1404   - cube.save("test_cta_edispcube_two.fits", true);
  1416 +
  1417 + // Test energy dispersion integration
  1418 + test_edispcube_integration(cube, 0.1, 10.0);
  1419 +
  1420 + // Test saving
  1421 + cube.save("test_cta_edispcube.fits", true);
1405 1422  
1406 1423 // Return
1407 1424 return;
... ... @@ -1688,6 +1705,51 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1688 1705  
1689 1706  
1690 1707 /***********************************************************************//**
  1708 + * @brief Utility function for energy dispersion cube integration tests
  1709 + *
  1710 + * @param[in] edisp Energy dispersion cube instance.
  1711 + * @param[in] e_src_min Minimum true energy (TeV).
  1712 + * @param[in] e_src_max Maximum true energy (TeV).
  1713 + *
  1714 + * Tests whether the energy dispersion cube integrated over observed
  1715 + * energies gives unity.
  1716 + ***************************************************************************/
  1717 +void TestGCTAResponse::test_edispcube_integration(const GCTACubeEdisp& edisp,
  1718 + const double& e_src_min,
  1719 + const double& e_src_max)
  1720 +{
  1721 + // Set sky direction
  1722 + GSkyDir dir;
  1723 + dir.radec_deg(83.63, 22.01);
  1724 +
  1725 + // Test energy dispersion cube integration
  1726 + for (double e_src = e_src_min; e_src <= e_src_max; e_src *= 2.0) {
  1727 +
  1728 + // Compute log10 of true energy
  1729 + GEnergy etrue(e_src, "TeV");
  1730 +
  1731 + // Perform numerical integration by summing
  1732 + const int nE = 10000;
  1733 + double emin = 1000.0; //!< 1 GeV
  1734 + double emax = 3.0 * etrue.MeV(); //!< Migra_max = 3.0
  1735 + double dE = (emax-emin)/nE;
  1736 + double sum = 0.0;
  1737 + double E_obs = emin;
  1738 + for (int i = 0; i < nE; ++i, E_obs += dE) {
  1739 + GEnergy ereco(E_obs,"MeV");
  1740 + double dp = edisp(ereco, etrue, dir);
  1741 + sum += dp * dE;
  1742 + }
  1743 + test_value(sum, 1.0, 0.001, "Response cube energy dispersion integration "
  1744 + " for "+etrue.print());
  1745 + }
  1746 +
  1747 + // Return
  1748 + return;
  1749 +}
  1750 +
  1751 +
  1752 +/***********************************************************************//**
1691 1753 * @brief Test CTA background
1692 1754 ***************************************************************************/
1693 1755 void TestGCTAModel::test_model_bgd(void)
... ...
inst/cta/test/test_CTA.hpp
  View file @f02bf3f
... ... @@ -99,6 +99,9 @@ public:
99 99 void test_edisp_integration(const GCTAEdisp& edisp,
100 100 const double& e_src_min = 0.1,
101 101 const double& e_src_max = 10.0);
  102 + void test_edispcube_integration(const GCTACubeEdisp& edisp,
  103 + const double& e_src_min = 0.1,
  104 + const double& e_src_max = 10.0);
102 105 };
103 106  
104 107  
... ...