Simon Bonnefoy / gammalib

Commit c23e96356e95d069fee3efb7820b5523486f1c69

Authored by Jürgen Knödlseder
1 parent 868ff9df

Adapt unit tests to revision of energy dispersion interface (#2704)

Inline Side-by-side
Showing 3 changed files with 67 additions and 84 deletions
inst/cta/test/test_CTA.cpp
  View file @c23e963
... ... @@ -901,7 +901,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
901 901 test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
902 902  
903 903 // Test that diagonal element is non-zero
904   - test_value(edisp2(std::log10(30.0),std::log10(30.0)), 14.7025305813298);
  904 + test_value(edisp2(std::log10(30.0),std::log10(30.0)), 2.12840930049511e-07);
905 905  
906 906 // Test if non-diagonal element (above diagonal) is zero
907 907 test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
... ... @@ -911,7 +911,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
911 911 GEnergy etrue(1.0, "TeV");
912 912 GEnergy emin(0.1, "TeV");
913 913 GEnergy emax(10.0, "TeV");
914   - test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.0);
  914 + test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.0);
915 915 test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.5);
916 916 test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.5);
917 917  
... ... @@ -958,7 +958,7 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
958 958 test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
959 959  
960 960 // Test that diagonal element is non-zero
961   - test_value(edisp2(std::log10(30.0),std::log10(30.0)), 7.4480009344395);
  961 + test_value(edisp2(std::log10(30.0),std::log10(30.0)), 1.0785985662216e-07);
962 962  
963 963 // Test if non-diagonal element (above diagonal) is zero
964 964 test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
... ... @@ -967,13 +967,13 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
967 967 GEnergy etrue(1.0, "TeV");
968 968 GEnergy emin(0.1, "TeV");
969 969 GEnergy emax(10.0, "TeV");
970   - test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 0.999664174482241);
971   - test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.494679101633454);
972   - test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.504985072848788);
  970 + test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.00002576227557);
  971 + test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.42065218271377);
  972 + test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.43096191010834);
973 973  
974 974 // Test mc method
975 975 GRan ran;
976   - test_value(edisp2.mc(ran, 0.0).TeV(), 0.942937281664539);
  976 + test_value(edisp2.mc(ran, 0.0).TeV(), 0.930126476703161);
977 977  
978 978 // Test normalisation
979 979 test_edisp_integration(edisp2, 0.5, 10.0);
... ... @@ -1006,14 +1006,14 @@ void TestGCTAResponse::test_response_edisp_2D(void)
1006 1006 test_value(edisp2.ebounds_obs(1.0).emax().TeV(), 13.2011115155942);
1007 1007  
1008 1008 // Test ebounds_src method
1009   - test_value(edisp2.ebounds_src(1.0).emin().TeV(), 7.72198100249182);
1010   - test_value(edisp2.ebounds_src(1.0).emax().TeV(), 19.8348918834691);
  1009 + test_value(edisp2.ebounds_src(1.0).emin().TeV(), 7.7813808018878);
  1010 + test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8059295273666);
1011 1011  
1012 1012 // Test if non-diagonal element (below diagonal) is zero
1013 1013 test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
1014 1014  
1015 1015 // Test that diagonal element is non-zero
1016   - test_value(edisp2(std::log10(30.0),std::log10(30.0)), 11.7884028286084);
  1016 + test_value(edisp2(std::log10(30.0),std::log10(30.0)), 1.74200329848547e-07);
1017 1017  
1018 1018 // Test if non-diagonal element (above diagonal) is zero
1019 1019 test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
... ... @@ -1022,13 +1022,13 @@ void TestGCTAResponse::test_response_edisp_2D(void)
1022 1022 GEnergy etrue(1.0, "TeV");
1023 1023 GEnergy emin(0.1, "TeV");
1024 1024 GEnergy emax(10.0, "TeV");
1025   - test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.00008878464952);
1026   - test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.448839360815326);
1027   - test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.551249618415791);
  1025 + test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.00002101900895);
  1026 + test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.480451962185927);
  1027 + test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.519569056823021);
1028 1028  
1029 1029 // Test mc method
1030 1030 GRan ran;
1031   - test_value(edisp2.mc(ran, 0.0).TeV(), 0.97532927330811);
  1031 + test_value(edisp2.mc(ran, 0.0).TeV(), 0.872387491302767);
1032 1032  
1033 1033 // Test normalisation
1034 1034 test_edisp_integration(edisp2);
... ... @@ -1411,9 +1411,9 @@ void TestGCTAResponse::test_response_edispcube(void)
1411 1411 /***********************************************************************//**
1412 1412 * @brief Utility function for energy dispersion tests
1413 1413 *
1414   - * @param[in] rsp Response
1415   - * @param[in] e_src_min Minimum true energy (TeV) (default: 0.1)
1416   - * @param[in] e_src_max Maximum true energy (TeV) (default: 10.0)
  1414 + * @param[in] rsp Response.
  1415 + * @param[in] e_src_min Minimum true energy (TeV).
  1416 + * @param[in] e_src_max Maximum true energy (TeV).
1417 1417 ***************************************************************************/
1418 1418 void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf& rsp,
1419 1419 const double& e_src_min,
... ... @@ -1429,40 +1429,29 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf& rs
1429 1429 double log10_e_src = std::log10(e_src);
1430 1430  
1431 1431 // Retrieve boundaries in observed energy
1432   - GEbounds ebounds = rsp.edisp()->ebounds_obs(log10_e_src);
1433   - GEnergy emin = ebounds.emin();
1434   - GEnergy emax = ebounds.emax();
1435   - double logE_min = std::log10(emin.TeV());
1436   - double logE_max = std::log10(emax.TeV());
1437   -
1438   - // Compute step size for numerical integration
1439   - const int steps = 1000;
1440   - double dlogE = (logE_max-logE_min)/steps;
  1432 + GEbounds ebounds = rsp.edisp()->ebounds_obs(log10_e_src);
  1433 + double emin = ebounds.emin().MeV();
  1434 + double emax = ebounds.emax().MeV();
1441 1435  
1442 1436 // Perform numerical integration by summing
1443   - double sum = 0.0;
1444   - double logE_obs = logE_min;
1445   - for (int i = 0; i < steps; ++i) {
1446   - double dp_dlogE = (*rsp.edisp())(logE_obs, log10_e_src);
1447   - sum += dp_dlogE * dlogE;
1448   - logE_obs += dlogE;
1449   - }
1450   - GEnergy eng(e_src, "TeV");
1451   - test_value(sum, 1.0, 0.001, "Energy Dispersion integration for "+eng.print());
1452   -
1453   - // And now in linear energies, test GCTAResponse::edisp() method
1454 1437 const int nE = 10000;
1455   - double E_min = emin.MeV();
1456   - double E_max = emax.MeV();
1457   - double dE = (E_max-E_min)/nE;
  1438 + double dE = (emax-emin)/nE;
  1439 + double sum1 = 0.0;
1458 1440 double sum2 = 0.0;
1459   - double E_obs = E_min;
1460   - for (int i = 0; i < nE; ++i) {
1461   - double dp_dE = rsp.edisp(GEnergy(E_obs,"MeV"), 0.0, 0.0, 0.0, 0.0, log10_e_src);
1462   - sum2 += dp_dE * dE;
1463   - E_obs += dE;
  1441 + double E_obs = emin;
  1442 + for (int i = 0; i < nE; ++i, E_obs += dE) {
  1443 + GEnergy ereco(E_obs,"MeV");
  1444 + double dp1 = (*rsp.edisp())(ereco.log10TeV(), log10_e_src);
  1445 + double dp2 = rsp.edisp(ereco, 0.0, 0.0, 0.0, 0.0, log10_e_src);
  1446 + sum1 += dp1 * dE;
  1447 + sum2 += dp2 * dE;
1464 1448 }
1465   - test_value(sum2, 1.0, 0.001, "Energy Dispersion integration for "+eng.print());
  1449 + GEnergy eng(e_src, "TeV");
  1450 + test_value(sum1, 1.0, 0.001, "Energy dispersion integration using "+
  1451 + rsp.edisp()->classname()+"::operator() "
  1452 + "for "+eng.print());
  1453 + test_value(sum2, 1.0, 0.001, "Energy dispersion integration using "
  1454 + "GCTAResponseIrf::edisp() for "+eng.print());
1466 1455 }
1467 1456  
1468 1457 } // endif: energy dispersion was available
... ... @@ -1475,15 +1464,17 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf&amp; rs
1475 1464 /***********************************************************************//**
1476 1465 * @brief Utility function for energy dispersion integration tests
1477 1466 *
1478   - * @param[in] edisp Energy dispersion
1479   - * @param[in] e_src_min Minimum true energy (TeV) (default: 0.1)
1480   - * @param[in] e_src_max Maximum true energy (TeV) (default: 10.0)
  1467 + * @param[in] edisp Energy dispersion instance.
  1468 + * @param[in] e_src_min Minimum true energy (TeV).
  1469 + * @param[in] e_src_max Maximum true energy (TeV).
  1470 + *
  1471 + * Tests whether the energy dispersion integrated over observed energies
  1472 + * gives unity.
1481 1473 ***************************************************************************/
1482   -void TestGCTAResponse::test_edisp_integration(const GCTAEdisp& edisp,
1483   - const double& e_src_min,
1484   - const double& e_src_max)
  1474 +void TestGCTAResponse::test_edisp_integration(const GCTAEdisp& edisp,
  1475 + const double& e_src_min,
  1476 + const double& e_src_max)
1485 1477 {
1486   -
1487 1478 // Loop over source energies
1488 1479 for (double e_src = e_src_min; e_src <= e_src_max; e_src *= 2.0) {
1489 1480  
... ... @@ -1498,30 +1489,22 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1498 1489 continue;
1499 1490 }
1500 1491  
1501   - // Set measured energy range
1502   - GEnergy emin = ebounds.emin();
1503   - GEnergy emax = ebounds.emax();
1504   - double logE_min = std::log10(emin.TeV());
1505   - double logE_max = std::log10(emax.TeV());
1506   -
1507   - // Compute step size for numerical integration
1508   - const int steps = 1000;
1509   - double dlogE = (logE_max-logE_min)/steps;
1510   -
1511 1492 // Perform numerical integration by summing
1512   - double sum = 0.0;
1513   - double logE_obs = logE_min;
1514   - for (int i = 0; i <= steps; ++i) {
1515   - double dp_dlogE = edisp(logE_obs, log10_e_src);
1516   - sum += dp_dlogE * dlogE;
1517   - logE_obs += dlogE;
  1493 + const int nE = 10000;
  1494 + double emin = ebounds.emin().MeV();
  1495 + double emax = ebounds.emax().MeV();
  1496 + double dE = (emax-emin)/nE;
  1497 + double sum = 0.0;
  1498 + double E_obs = emin;
  1499 + for (int i = 0; i < nE; ++i, E_obs += dE) {
  1500 + GEnergy ereco(E_obs,"MeV");
  1501 + double dp = edisp(ereco.log10TeV(), log10_e_src);
  1502 + sum += dp * dE;
1518 1503 }
1519   -
1520   - // Set message string
1521   - GEnergy eng(e_src, "TeV");
1522   - std::string msg = edisp.classname()+" integration for "+eng.print();
1523   - test_value(sum, 1.0, 0.01, msg);
1524   -
  1504 + GEnergy eng(e_src, "TeV");
  1505 + test_value(sum, 1.0, 0.001, "Energy dispersion integration using "+
  1506 + edisp.classname()+"::operator() for "+
  1507 + eng.print());
1525 1508 }
1526 1509  
1527 1510 // Return
... ...
inst/cta/test/test_CTA.hpp
  View file @c23e963
... ... @@ -95,9 +95,9 @@ public:
95 95 void test_response_edisp_integration(const GCTAResponseIrf& rsp,
96 96 const double& e_src_min = 0.1,
97 97 const double& e_src_max = 10.0);
98   - void test_edisp_integration(const GCTAEdisp& edisp,
99   - const double& e_src_min = 0.1,
100   - const double& e_src_max = 10.0);
  98 + void test_edisp_integration(const GCTAEdisp& edisp,
  99 + const double& e_src_min = 0.1,
  100 + const double& e_src_max = 10.0);
101 101 };
102 102  
103 103  
... ...
inst/cta/test/test_CTA.py
  View file @c23e963
... ... @@ -250,12 +250,12 @@ class Test(gammalib.GPythonTestSuite):
250 250 edisp = gammalib.GCTAEdispPerfTable(filename)
251 251  
252 252 # Test energy dispersion values
253   - self.test_value(edisp(0.0, 0.0), 9.99019627861, 1.0e-6)
254   - self.test_value(edisp(0.001, 0.0), 9.9870644077, 1.0e-6)
255   - self.test_value(edisp(0.01, 0.0), 9.68182, 1.0e-6)
256   - self.test_value(edisp(0.1, 0.0), 0.434382, 1.0e-6)
257   - self.test_value(edisp(1.0, 1.0), 18.064868197, 1.0e-6)
258   - self.test_value(edisp(1.001, 1.0, 0.0), 18.0463571212, 1.0e-6)
  253 + self.test_value(edisp(0.0, 0.0), 4.3386871e-06, 1.0e-6)
  254 + self.test_value(edisp(0.001, 0.0), 4.3273514e-06, 1.0e-6)
  255 + self.test_value(edisp(0.01, 0.0), 4.1090489e-06, 1.0e-6)
  256 + self.test_value(edisp(0.1, 0.0), 1.4984969e-07, 1.0e-6)
  257 + self.test_value(edisp(1.0, 1.0), 7.8454726e-07, 1.0e-6)
  258 + self.test_value(edisp(1.001, 1.0, 0.0), 7.8194077e-07, 1.0e-6)
259 259  
260 260 # Test GCTAResponseIrf file constructor
261 261 db = gammalib.GCaldb(self._caldb)
... ...