Thierry Louge / gammalib

Commit 555c1e7cd110594b264e6c45f317eef16a02b096

Authored by Jürgen Knödlseder
1 parent 8c9cb5bf

Create new test files for stacked analysis 

The test files for stacked analysis are now based on 3 observations of the Crab performed at different aspect angles, simulated using the Prod2 South_0.5h IRFs.
Inline Side-by-side
Showing 14 changed files with 98 additions and 58 deletions
inst/cta/test/data/bkgcube.fits.gz deleted
View file @8c9cb5b
No preview for this file type
inst/cta/test/data/edispcube.fits.gz deleted
View file @8c9cb5b
No preview for this file type
inst/cta/test/data/expcube.fits.gz deleted
View file @8c9cb5b
No preview for this file type
inst/cta/test/data/obs_cube.xml deleted
View file @8c9cb5b
1   -<observation_list title="observation library">
2   - <observation name="Crab" id="00001" instrument="CTA">
3   - <parameter name="CountsCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/crab_cntmap.fits"/>
4   - <parameter name="ExposureCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/expcube.fits"/>
5   - <parameter name="PsfCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/psfcube.fits"/>
6   - <parameter name="BkgCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/bkgcube.fits"/>
7   - </observation>
8   -</observation_list>
inst/cta/test/data/psfcube.fits.gz deleted
View file @8c9cb5b
No preview for this file type
inst/cta/test/data/stacked_bkgcube.fits.gz 0 → 100644
View file @555c1e7
No preview for this file type
inst/cta/test/data/stacked_cntcube.fits.gz 0 → 100644
View file @555c1e7
No preview for this file type
inst/cta/test/data/stacked_edispcube.fits.gz 0 → 100644
View file @555c1e7
No preview for this file type
inst/cta/test/data/stacked_expcube.fits.gz 0 → 100644
View file @555c1e7
No preview for this file type
inst/cta/test/data/crab_stacked.xml renamed to inst/cta/test/data/stacked_model.xml
  View file @555c1e7
1 1 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
2 2 <source_library title="source library">
3   - <source name="MAGIC Crab (Albert et al. 2008, ApJ, 674, 1037)" type="PointSource">
  3 + <source name="Crab" type="PointSource">
4 4 <spectrum type="PowerLaw">
5 5 <parameter name="Prefactor" value="5.7" error="0" scale="1e-16" min="1e-07" max="1000" free="1" />
6 6 <parameter name="Index" value="2.48" error="0" scale="-1" min="0" max="5" free="1" />
... ...
inst/cta/test/data/stacked_obs.xml 0 → 100644
  View file @555c1e7
  1 +<observation_list title="observation library">
  2 + <observation name="Crab" id="00001" instrument="CTA">
  3 + <parameter name="CountsCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/stacked_cntcube.fits"/>
  4 + <parameter name="ExposureCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/stacked_expcube.fits"/>
  5 + <parameter name="PsfCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/stacked_psfcube.fits"/>
  6 + <parameter name="BkgCube" file="$(PACKAGE_SOURCE)/inst/cta/test/data/stacked_bkgcube.fits"/>
  7 + </observation>
  8 +</observation_list>
... ...
inst/cta/test/data/stacked_obs_input.xml 0 → 100644
  View file @555c1e7
  1 +<?xml version="1.0" encoding="UTF-8" standalone="no"?>
  2 +<observation_list title="observation list">
  3 + <observation name="On Source" id="000001" instrument="CTA">
  4 + <parameter name="Pointing" ra="83.6331" dec="22.0145" />
  5 + <parameter name="EnergyBoundaries" emin="100000.0" emax="100000000.0" />
  6 + <parameter name="GoodTimeIntervals" tmin="0" tmax="1100.0" />
  7 + <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
  8 + <parameter name="RegionOfInterest" ra="83.6331" dec="22.0145" rad="1" />
  9 + <parameter name="Deadtime" deadc="0.95" />
  10 + <parameter name="Calibration" database="prod2" response="South_0.5h" />
  11 + </observation>
  12 + <observation name="Negative Offset" id="000002" instrument="CTA">
  13 + <parameter name="Pointing" ra="83.6331" dec="21.0145" />
  14 + <parameter name="EnergyBoundaries" emin="300000.0" emax="10000000.0" />
  15 + <parameter name="GoodTimeIntervals" tmin="0" tmax="1000.0" />
  16 + <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
  17 + <parameter name="RegionOfInterest" ra="83.6331" dec="21.0145" rad="1.2" />
  18 + <parameter name="Deadtime" deadc="0.95" />
  19 + <parameter name="Calibration" database="prod2" response="South_0.5h" />
  20 + </observation>
  21 + <observation name="Positive Offset" id="000003" instrument="CTA">
  22 + <parameter name="Pointing" ra="84.1331" dec="22.5145" />
  23 + <parameter name="EnergyBoundaries" emin="200000.0" emax="80000000.0" />
  24 + <parameter name="GoodTimeIntervals" tmin="0" tmax="900.0" />
  25 + <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
  26 + <parameter name="RegionOfInterest" ra="84.1331" dec="22.5145" rad="1.5" />
  27 + <parameter name="Deadtime" deadc="0.95" />
  28 + <parameter name="Calibration" database="prod2" response="South_0.5h" />
  29 + </observation>
  30 +</observation_list>
... ...
inst/cta/test/data/stacked_psfcube.fits.gz 0 → 100644
View file @555c1e7
No preview for this file type
inst/cta/test/test_CTA.cpp
  View file @555c1e7
... ... @@ -46,30 +46,36 @@
46 46 /* __ Globals ____________________________________________________________ */
47 47  
48 48 /* __ Constants __________________________________________________________ */
49   -const std::string datadir = PACKAGE_SOURCE"/inst/cta/test/data";
50   -const std::string cta_caldb = PACKAGE_SOURCE"/inst/cta/caldb";
51   -const std::string cta_irf = "cta_dummy_irf";
52   -const std::string cta_events = datadir+"/crab_events.fits";
53   -const std::string cta_events_gti = datadir+"/crab_events_gti.fits[EVENTS2]";
54   -const std::string cta_cntmap = datadir+"/crab_cntmap.fits";
55   -const std::string cta_cube_xml = datadir+"/obs_cube.xml";
56   -const std::string cta_bin_xml = datadir+"/obs_binned.xml";
57   -const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml";
58   -const std::string cta_model_xml = datadir+"/crab.xml";
59   -const std::string cta_stack_model = datadir+"/crab_stacked.xml";
60   -const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
61   -const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
62   -const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
63   -const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
64   -const std::string cta_caldb_king = PACKAGE_SOURCE"/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h_King";
65   -const std::string cta_irf_king = "irf_file.fits";
66   -const std::string cta_psf_table = PACKAGE_SOURCE"/inst/cta/test/caldb/psf_table.fits[PSF_2D_TABLE]";
67   -const std::string cta_edisp_perf = PACKAGE_SOURCE"/inst/cta/test/caldb/cta_dummy_irf.dat";
68   -const std::string cta_edisp_rmf = PACKAGE_SOURCE"/inst/cta/test/caldb/dc1/rmf.fits";
69   -const std::string cta_edisp_2D = PACKAGE_SOURCE"/inst/cta/test/caldb/edisp_matrix.fits";
70   -const std::string cta_bkgcube = datadir+"/bkgcube.fits";
71   -const std::string cta_modbck_fit = datadir+"/bg_test.fits";
72   -const std::string cta_point_table = datadir+"/crab_pointing.fits";
  49 +const std::string datadir = PACKAGE_SOURCE"/inst/cta/test/data";
  50 +const std::string cta_caldb = PACKAGE_SOURCE"/inst/cta/caldb";
  51 +const std::string cta_irf = "cta_dummy_irf";
  52 +const std::string cta_events = datadir+"/crab_events.fits";
  53 +const std::string cta_events_gti = datadir+"/crab_events_gti.fits[EVENTS2]";
  54 +const std::string cta_cntmap = datadir+"/crab_cntmap.fits";
  55 +const std::string cta_bin_xml = datadir+"/obs_binned.xml";
  56 +const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml";
  57 +const std::string cta_model_xml = datadir+"/crab.xml";
  58 +const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
  59 +const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
  60 +const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
  61 +const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
  62 +const std::string cta_caldb_king = PACKAGE_SOURCE"/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h_King";
  63 +const std::string cta_irf_king = "irf_file.fits";
  64 +const std::string cta_psf_table = PACKAGE_SOURCE"/inst/cta/test/caldb/psf_table.fits[PSF_2D_TABLE]";
  65 +const std::string cta_edisp_perf = PACKAGE_SOURCE"/inst/cta/test/caldb/cta_dummy_irf.dat";
  66 +const std::string cta_edisp_rmf = PACKAGE_SOURCE"/inst/cta/test/caldb/dc1/rmf.fits";
  67 +const std::string cta_edisp_2D = PACKAGE_SOURCE"/inst/cta/test/caldb/edisp_matrix.fits";
  68 +const std::string cta_modbck_fit = datadir+"/bg_test.fits";
  69 +const std::string cta_point_table = datadir+"/crab_pointing.fits";
  70 +
  71 +/* __ Test files for stacked analysis (based on Prod2::South_0.5h) _______ */
  72 +const std::string cta_stacked_xml = datadir+"/stacked_obs.xml";
  73 +const std::string cta_stacked_model = datadir+"/stacked_model.xml";
  74 +const std::string cta_stacked_cntcube = datadir+"/stacked_cntcube.fits";
  75 +const std::string cta_stacked_expcube = datadir+"/stacked_expcube.fits";
  76 +const std::string cta_stacked_psfcube = datadir+"/stacked_psfcube.fits";
  77 +const std::string cta_stacked_edispcube = datadir+"/stacked_edispcube.fits";
  78 +const std::string cta_stacked_bkgcube = datadir+"/stacked_bkgcube.fits";
73 79  
74 80  
75 81 /***********************************************************************//**
... ... @@ -1022,7 +1028,7 @@ void TestGCTAResponse::test_response_psfcube(void)
1022 1028 }
1023 1029  
1024 1030 /***********************************************************************//**
1025   - * @brief Test PSF cube handling
  1031 + * @brief Test background cube handling
1026 1032 ***************************************************************************/
1027 1033 void TestGCTAResponse::test_response_bkgcube(void)
1028 1034 {
... ... @@ -1037,7 +1043,7 @@ void TestGCTAResponse::test_response_bkgcube(void)
1037 1043 }
1038 1044  
1039 1045 GCTACubeBackground cube;
1040   - cube.load(cta_bkgcube);
  1046 + cube.load(cta_stacked_bkgcube);
1041 1047 cube.save("test_cta_bkgcube.fits", true);
1042 1048  
1043 1049 // Return
... ... @@ -1045,7 +1051,7 @@ void TestGCTAResponse::test_response_bkgcube(void)
1045 1051 }
1046 1052  
1047 1053 /***********************************************************************//**
1048   - * @brief Test Edisp cube handling
  1054 + * @brief Test energy dispersion cube handling
1049 1055 ***************************************************************************/
1050 1056 void TestGCTAResponse::test_response_edispcube(void)
1051 1057 {
... ... @@ -1745,14 +1751,13 @@ void TestGCTAObservation::test_binned_obs(void)
1745 1751 ***************************************************************************/
1746 1752 void TestGCTAObservation::test_stacked_obs(void)
1747 1753 {
1748   - // Set filenames
1749   - const std::string filename = "test_cta_obs_cube.xml";
  1754 + // Construct stacked observation without energy dispersion
  1755 + GCTAObservation cta1(cta_stacked_cntcube,
  1756 + cta_stacked_expcube,
  1757 + cta_stacked_psfcube,
  1758 + cta_stacked_bkgcube);
1750 1759  
1751   - // Construct stacked observation
1752   - GCTAObservation cta1(cta_cntmap,
1753   - datadir+"/expcube.fits",
1754   - datadir+"/psfcube.fits",
1755   - datadir+"/bkgcube.fits");
  1760 + // Test for presence of response
1756 1761 const GCTAResponseCube* rsp =
1757 1762 dynamic_cast<const GCTAResponseCube*>(cta1.response());
1758 1763 test_assert((rsp != NULL), "Observation contains cube response");
... ... @@ -1762,11 +1767,13 @@ void TestGCTAObservation::test_stacked_obs(void)
1762 1767 }
1763 1768  
1764 1769 // Construct stacked observation with energy dispersion
1765   - GCTAObservation cta2(cta_cntmap,
1766   - datadir+"/expcube.fits",
1767   - datadir+"/psfcube.fits",
1768   - datadir+"/edispcube.fits",
1769   - datadir+"/bkgcube.fits");
  1770 + GCTAObservation cta2(cta_stacked_cntcube,
  1771 + cta_stacked_expcube,
  1772 + cta_stacked_psfcube,
  1773 + cta_stacked_edispcube,
  1774 + cta_stacked_bkgcube);
  1775 +
  1776 + // Test for presence of response
1770 1777 rsp = dynamic_cast<const GCTAResponseCube*>(cta2.response());
1771 1778 test_assert((rsp != NULL), "Observation contains cube response");
1772 1779 if (rsp != NULL) {
... ... @@ -1774,9 +1781,11 @@ void TestGCTAObservation::test_stacked_obs(void)
1774 1781 test_assert(rsp->use_edisp(), "Response has energy dispersion");
1775 1782 }
1776 1783  
1777   - // Construct stacked observation from XML file (w/o energy dispersion)
1778   - GObservations obs(cta_cube_xml);
  1784 + // Construct stacked observation from XML file without energy dispersion
  1785 + GObservations obs(cta_stacked_xml);
1779 1786 test_value(obs.size(), 1, "One observation in container");
  1787 +
  1788 + // Test for presence of response
1780 1789 rsp = dynamic_cast<const GCTAResponseCube*>(obs[0]->response());
1781 1790 test_assert((rsp != NULL), "Observation contains cube response");
1782 1791 if (rsp != NULL) {
... ... @@ -1785,7 +1794,7 @@ void TestGCTAObservation::test_stacked_obs(void)
1785 1794 }
1786 1795  
1787 1796 // Save observation container into XML file
1788   - obs.save(filename);
  1797 + obs.save("test_cta_obs_cube.xml");
1789 1798  
1790 1799 // Return
1791 1800 return;
... ... @@ -1907,27 +1916,28 @@ void TestGCTAOptimize::test_binned_optimizer(void)
1907 1916 /***********************************************************************//**
1908 1917 * @brief Test stacked optimizer
1909 1918 *
1910   - * The stacked response cubes have been computed using the dummy_irf.
  1919 + * The stacked response cubes have been computed using the Prod2::South_0.5h
  1920 + * response.
1911 1921 ***************************************************************************/
1912 1922 void TestGCTAOptimize::test_stacked_optimizer(void)
1913 1923 {
1914 1924 // Set reference result
1915 1925 double fit_results[] = {83.6331, 0,
1916 1926 22.0145, 0,
1917   - 5.98597e-16, 2.00288e-17,
1918   - -2.5012, 0.0250613,
  1927 + 5.98351e-16, 1.09268e-17,
  1928 + -2.49903, 0.0145957,
1919 1929 300000, 0,
1920 1930 1, 0,
1921   - 0.978511, 0.0255432,
1922   - 0.0290683, 0.0163737,
  1931 + 1.07386, 0.0218562,
  1932 + 0.00440413, 0.0145699,
1923 1933 1.0e6, 0,
1924 1934 1, 0};
1925 1935  
1926 1936 // Load stacked CTA observation
1927   - GObservations obs(cta_cube_xml);
  1937 + GObservations obs(cta_stacked_xml);
1928 1938  
1929 1939 // Load models from XML file
1930   - obs.models(cta_stack_model);
  1940 + obs.models(cta_stacked_model);
1931 1941  
1932 1942 // Perform LM optimization
1933 1943 GOptimizerLM opt;
... ...