Leonardo Di Venere / ctools

Commit b1668ffb741935eed09dea7c7c3c734b8ed6fd7e

Authored by Jürgen Knödlseder
1 parent fc6ba477

Implement FFT computation for RING background (#2309) 

A FFT-based computation was implemented in ctskymap for the RING background, which is faster but less accurate than the direct computation.
Inline Side-by-side
Showing 4 changed files with 307 additions and 46 deletions
doc/source/users/reference_manual/ctskymap.rst
  View file @b1668ff
... ... @@ -93,6 +93,13 @@ Background subtraction configuration parameters
93 93 ``(inexclusion = NONE) [file]``
94 94 Exclusion region file as either a FITS map or DS9 region file.
95 95  
  96 +``(usefft = yes) [boolean]``
  97 + Specifies whether a FFT should be used for ring computation. If set to
  98 + ``no`` the ring background will be computed directly. Direct computation is
  99 + more accurate since the precise angular distances between pixels are taken
  100 + into account, but FFT is faster, but less accurate since it assumes
  101 + Euclidean distances between pixels.
  102 +
96 103  
97 104 Standard parameters
98 105 -------------------
... ...
src/ctskymap/ctskymap.cpp
  View file @b1668ff
... ... @@ -35,6 +35,7 @@
35 35 #define G_GET_PARAMETERS "ctskymap::get_parameter()"
36 36 #define G_MAP_EVENTS "ctskymap::map_events(GCTAObservation*)"
37 37 #define G_MAP_BACKGROUND_IRF "ctskymap::map_background_irf(GCTAObservation*)"
  38 +#define G_RING_KERNEL "ctskymap::ring_kernel(double&, double&)"
38 39  
39 40 /* __ Debug definitions __________________________________________________ */
40 41  
... ... @@ -407,6 +408,7 @@ void ctskymap::init_members(void)
407 408 m_roiradius = 0.0;
408 409 m_inradius = 0.0;
409 410 m_outradius = 0.0;
  411 + m_usefft = true;
410 412 m_publish = false;
411 413 m_chatter = static_cast<GChatter>(2);
412 414  
... ... @@ -445,6 +447,7 @@ void ctskymap::copy_members(const ctskymap&amp; app)
445 447 m_roiradius = app.m_roiradius;
446 448 m_inradius = app.m_inradius;
447 449 m_outradius = app.m_outradius;
  450 + m_usefft = app.m_usefft;
448 451 m_publish = app.m_publish;
449 452 m_chatter = app.m_chatter;
450 453  
... ... @@ -508,6 +511,7 @@ void ctskymap::get_parameters(void)
508 511 m_inradius = (*this)["inradius"].real();
509 512 m_outradius = (*this)["outradius"].real();
510 513 m_inexclusion = (*this)["inexclusion"].filename();
  514 + m_usefft = (*this)["usefft"].boolean();
511 515  
512 516 // Make sure that (roiradius < inradius < outradius)
513 517 if (m_roiradius > m_inradius) {
... ... @@ -911,19 +915,42 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
911 915  
912 916  
913 917 /***********************************************************************//**
914   - * @brief Computes the bin-by-bin significance if background subtraction
  918 + * @brief Compute the pixel significance if background subtraction is
915 919 * requested
916 920 *
917   - * Method computes the bin-by bin significance if a background subtraction
918   - * method is specified. If method is "IRF", Poisson statistics (in the
919   - * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance
920   - * is computed for each bin.
  921 + * Computes the pixel significance if a background subtraction method is
  922 + * specified.
  923 + *
  924 + * For the @c IRF method Poisson statistic (in the Gaussian limit) is
  925 + * assumed, and the significance \f$\sigma_i\f$ is computed using
  926 + *
  927 + * \f[\sigma_i = \frac{N_i - B_i}{\sqrt{N_i}}\f]
  928 + *
  929 + * where
  930 + * \f$N_i\f$ is the number of observed counts and
  931 + * \f$B_i\f$ is the estimated number of background counts.
  932 + *
  933 + * For the @c RING method the Li & Ma significance is computed for each
  934 + * pixel, using the @c roiradius, @c inradius and @c outradius parameters
  935 + * to specify the radius of the On region and the Off background ring,
  936 + * respectively. Depending on the @c usefft parameter, either an FFT is used
  937 + * to compute the number of events and the acceptance of the On and Off
  938 + * regions, or a direct computation is performed. While the former is faster
  939 + * it is less accurate since it assumes Euclidean distances in the sky map.
  940 + * The latter is slower but compute exact distance between pixels of the sky
  941 + * map. Use FFT if the sky map is close to a cartesian grid, and use the
  942 + * direct method if the sky map shows important distortions.
921 943 ***************************************************************************/
922 944 void ctskymap::map_significance(void)
923 945 {
924 946 // Compute significance from "RING" method (Li & Ma eq. 17)
925 947 if (m_bkgsubtract == "RING") {
926   - map_significance_ring();
  948 + if (m_usefft) {
  949 + map_significance_ring_fft();
  950 + }
  951 + else {
  952 + map_significance_ring_direct();
  953 + }
927 954 }
928 955  
929 956 // Compute significance using Poisson statistics in the Gaussian limit
... ... @@ -937,15 +964,16 @@ void ctskymap::map_significance(void)
937 964  
938 965  
939 966 /***********************************************************************//**
940   - * @brief Computes the bin-by-bin significance for RING background
  967 + * @brief Compute the pixel significance for RING background
941 968 *
942   - * Computes the a Li & Ma significance for each sky map bin and replaces
943   - * the sky and background maps by the On- and Off-count maps.
  969 + * Computes the Li & Ma significance for each sky map pixel and replaces
  970 + * the sky and background maps by the On- and Off-count maps. The computation
  971 + * is done by computing the exact distances between pixels in the sky map.
944 972 ***************************************************************************/
945   -void ctskymap::map_significance_ring(void)
  973 +void ctskymap::map_significance_ring_direct(void)
946 974 {
947 975 // Log message about what is being done
948   - log_header1(NORMAL, "Computing Ring background map");
  976 + log_header1(NORMAL, "Computing ring background map (direct method)");
949 977 log_value(NORMAL, "Total pixels to process", m_skymap.npix());
950 978  
951 979 // Initialise On- and Off-count maps
... ... @@ -990,21 +1018,15 @@ void ctskymap::map_significance_ring(void)
990 1018 // If alpha is zero then increment the bad alpha counter
991 1019 if (alpha == 0.0) {
992 1020 num_bad_alpha++;
  1021 + #pragma omp critical(ctskymap_map_significance_ring)
  1022 + m_sigmap(i) = 0.0;
993 1023 }
994 1024  
995 1025 // ... otherwise store the results
996 1026 else {
997 1027  
998 1028 // Compute and store significance (Li & Ma eq. 17)
999   - double sigma;
1000   - if (n_on == 0.0) {
1001   - sigma = -2.0 * n_off * std::log(1.0+alpha);
1002   - }
1003   - else {
1004   - sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
1005   - (n_on * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) +
1006   - n_off * std::log((1.0+alpha) * n_off / (n_on+n_off)));
1007   - }
  1029 + double sigma = sigma_li_ma(n_on, n_off, alpha);
1008 1030 #pragma omp critical(ctskymap_map_significance_ring)
1009 1031 m_sigmap(i) = sigma;
1010 1032  
... ... @@ -1013,7 +1035,7 @@ void ctskymap::map_significance_ring(void)
1013 1035 } // endfor: looped over all pixels
1014 1036  
1015 1037 // Log the number of bad-alpha bins
1016   - log_value(NORMAL, "Bins with alpha=0", num_bad_alpha);
  1038 + log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha);
1017 1039  
1018 1040 // Now take the square root. Since some bins can be negative, first
1019 1041 // take the absolute value of the map, square root that, then multiply
... ... @@ -1031,6 +1053,90 @@ void ctskymap::map_significance_ring(void)
1031 1053  
1032 1054  
1033 1055 /***********************************************************************//**
  1056 + * @brief Compute the pixel significance for RING background
  1057 + *
  1058 + * Computes the Li & Ma significance for each sky map pixel and replaces
  1059 + * the sky and background maps by the On- and Off-count maps. The computation
  1060 + * is done using a FFT.
  1061 + ***************************************************************************/
  1062 +void ctskymap::map_significance_ring_fft(void)
  1063 +{
  1064 + // Log message about what is being done
  1065 + log_header1(NORMAL, "Computing ring background map (FFT method)");
  1066 +
  1067 + // Initialise statistics
  1068 + int num_bad_alpha = 0;
  1069 +
  1070 + // Set-up sky and acceptance maps with excluded pixels according to the
  1071 + // exclusion map
  1072 + GSkyMap skymap = m_skymap;
  1073 + GSkyMap acceptance = m_bkgmap;
  1074 + for (int i = 0; i < m_skymap.npix(); ++i) {
  1075 + if (m_exclmap(i) != 0.0) {
  1076 + skymap(i) = 0.0;
  1077 + acceptance(i) = 0.0;
  1078 + }
  1079 + }
  1080 +
  1081 + // Convolve maps for On-region and Off-region
  1082 + GSkyMap on_counts = ring_convolve(skymap, 0.0, m_roiradius);
  1083 + GSkyMap on_alpha = ring_convolve(acceptance, 0.0, m_roiradius);
  1084 + GSkyMap off_counts = ring_convolve(skymap, m_inradius, m_outradius);
  1085 + GSkyMap off_alpha = ring_convolve(acceptance, m_inradius, m_outradius);
  1086 +
  1087 + // Loop over sky map pixels
  1088 + for (int i = 0; i < m_skymap.npix(); ++i) {
  1089 +
  1090 + // Get On-counts, Off-counts and alpha value for this bin
  1091 + double n_on = on_counts(i);
  1092 + double n_off = off_counts(i);
  1093 + double alpha;
  1094 + if (off_alpha(i) == 0.0) {
  1095 + alpha = 0.0;
  1096 + n_on = 0.0;
  1097 + n_off = 0.0;
  1098 + }
  1099 + else {
  1100 + alpha = on_alpha(i) / off_alpha(i);
  1101 + }
  1102 +
  1103 + // Store the On and alpha-weighted Off-counts
  1104 + on_counts(i) = n_on;
  1105 + off_counts(i) = alpha * n_off;
  1106 +
  1107 + // If alpha is zero then increment the bad alpha counter
  1108 + if (alpha == 0.0) {
  1109 + m_sigmap(i) = 0.0;
  1110 + num_bad_alpha++;
  1111 + }
  1112 +
  1113 + // ... otherwise compute and store the results
  1114 + else {
  1115 + m_sigmap(i) = sigma_li_ma(n_on, n_off, alpha);
  1116 + }
  1117 +
  1118 + } // endfor: looped over all pixels
  1119 +
  1120 + // Log the number of bad-alpha bins
  1121 + log_value(NORMAL, "Total number of pixels", m_skymap.npix());
  1122 + log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha);
  1123 +
  1124 + // Now take the square root. Since some bins can be negative, first
  1125 + // take the absolute value of the map, square root that, then multiply
  1126 + // each bin by its sign to preserve the +/- significance.
  1127 + m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap));
  1128 +
  1129 + // Terminate by storing the on-counts map as sky map and the off-counts
  1130 + // map as background map
  1131 + m_skymap = on_counts;
  1132 + m_bkgmap = off_counts;
  1133 +
  1134 + // Return
  1135 + return;
  1136 +}
  1137 +
  1138 +
  1139 +/***********************************************************************//**
1034 1140 * @brief Computes Non, Noff and alpha for a counts map and sensitivity map
1035 1141 *
1036 1142 * @param[in] ipixel Sky map pixel to consider.
... ... @@ -1172,7 +1278,7 @@ void ctskymap::compute_ring_values(const int&amp; ipixel,
1172 1278 * ny is the number of y pixels in the sky map.
1173 1279 ***************************************************************************/
1174 1280 void ctskymap::ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
1175   - int& iy1, int& iy2)
  1281 + int& iy1, int& iy2) const
1176 1282 {
1177 1283 // Get number of pixels in x and y direction
1178 1284 int nx = m_skymap.nx();
... ... @@ -1241,6 +1347,145 @@ void ctskymap::ring_bounding_box(const int&amp; ipixel, int&amp; ix1, int&amp; ix2,
1241 1347  
1242 1348  
1243 1349 /***********************************************************************//**
  1350 + * @brief Return FFT kernel for background ring
  1351 + *
  1352 + * @param[in] map Sky map to be convolved with a background ring.
  1353 + * @param[in] rmin Minimum ring radius (degrees).
  1354 + * @param[in] rmax Maximum ring radius (degrees).
  1355 + * @return FFT kernel.
  1356 + *
  1357 + * Computes the FFT kernel for a backgrund ring.
  1358 + ***************************************************************************/
  1359 +GSkyMap ctskymap::ring_convolve(const GSkyMap& map, const double& rmin,
  1360 + const double& rmax) const
  1361 +{
  1362 + // Copy input map
  1363 + GSkyMap conv_map = map;
  1364 +
  1365 + // Get FFT of ring kernel
  1366 + GFft fft_kernel(ring_kernel(rmin, rmax));
  1367 +
  1368 + // Extract sky map into Ndarray
  1369 + GNdarray array(conv_map.nx(), conv_map.ny());
  1370 + double *dst = array.data();
  1371 + for (int i = 0; i < conv_map.npix(); ++i) {
  1372 + *dst++ = conv_map(i);
  1373 + }
  1374 +
  1375 + // FFT of sky map
  1376 + GFft fft_array = GFft(array);
  1377 +
  1378 + // Multiply FFT of sky map with FFT of kernel
  1379 + GFft fft_smooth = fft_array * fft_kernel;
  1380 +
  1381 + // Backward transform sky map
  1382 + GNdarray smooth = fft_smooth.backward();
  1383 +
  1384 + // Insert sky map
  1385 + double *src = smooth.data();
  1386 + for (int i = 0; i < conv_map.npix(); ++i) {
  1387 + conv_map(i) = *src++;
  1388 + }
  1389 +
  1390 + // Return convolved map
  1391 + return conv_map;
  1392 +}
  1393 +
  1394 +
  1395 +/***********************************************************************//**
  1396 + * @brief Return FFT kernel for background ring
  1397 + *
  1398 + * @param[in] rmin Minimum ring radius (degrees).
  1399 + * @param[in] rmax Maximum ring radius (degrees).
  1400 + * @return FFT kernel.
  1401 + *
  1402 + * @exception GException::invalid_value
  1403 + * Sky map is not a WCS projection.
  1404 + *
  1405 + * Computes the FFT kernel for a background ring. The computation is done
  1406 + * assuming that the sky map represents a cartesian grid, and distances are
  1407 + * computed in that grid using Euclidean distances. The @p rmin and @p rmax
  1408 + * arguments refer to a ring radius in these Euclidean distances.
  1409 + ***************************************************************************/
  1410 +GNdarray ctskymap::ring_kernel(const double& rmin, const double& rmax) const
  1411 +{
  1412 + // Get pointer on WCS projection
  1413 + const GWcs* wcs = dynamic_cast<const GWcs*>(m_skymap.projection());
  1414 + if (wcs == NULL) {
  1415 + std::string msg = "Sky map is not a WCS projection. Method is only "
  1416 + "valid for WCS projections.";
  1417 + throw GException::invalid_value(G_RING_KERNEL, msg);
  1418 + }
  1419 +
  1420 + // Get X and Y step size
  1421 + double dx = wcs->cdelt(0);
  1422 + double dy = wcs->cdelt(1);
  1423 +
  1424 + // Initialise kernel
  1425 + GNdarray kern(m_skymap.nx(), m_skymap.ny());
  1426 +
  1427 + // Fill kernel
  1428 + for (int ix1 = 0, ix2 = m_skymap.nx(); ix1 < m_skymap.nx(); ++ix1, --ix2) {
  1429 + double x = ix1 * dx;
  1430 + double xqs = x * x;
  1431 + for (int iy1 = 0, iy2 = m_skymap.ny(); iy1 < m_skymap.ny(); ++iy1, --iy2) {
  1432 + double y = iy1 * dy;
  1433 + double r = std::sqrt(xqs + y*y);
  1434 + if ((r >= rmin) && (r <= rmax)) {
  1435 + kern(ix1,iy1) += 1.0;
  1436 + if (ix2 < m_skymap.nx()) {
  1437 + kern(ix2,iy1) += 1.0;
  1438 + }
  1439 + if (iy2 < m_skymap.ny()) {
  1440 + kern(ix1,iy2) += 1.0;
  1441 + }
  1442 + if ((ix2 < m_skymap.nx()) && (iy2 < m_skymap.ny())) {
  1443 + kern(ix2,iy2) += 1.0;
  1444 + }
  1445 + }
  1446 + }
  1447 + }
  1448 +
  1449 + // Return kernel
  1450 + return kern;
  1451 +}
  1452 +
  1453 +
  1454 +/***********************************************************************//**
  1455 + * @brief Compute significance following Li & Ma
  1456 + *
  1457 + * @param[in] n_on Number of On-counts.
  1458 + * @param[in] n_off Number of Off-counts.
  1459 + * @param[in] alpha Alpha parameters.
  1460 + * @return Significance.
  1461 + *
  1462 + * Computes the significance following Li & Ma, Equation (17).
  1463 + ***************************************************************************/
  1464 +double ctskymap::sigma_li_ma(const double& n_on,
  1465 + const double& n_off,
  1466 + const double& alpha) const
  1467 +{
  1468 + // Allocate result
  1469 + double sigma;
  1470 +
  1471 + // Handle special case of no On-counts
  1472 + if (n_on == 0.0) {
  1473 + sigma = -2.0 * n_off * std::log(1.0+alpha);
  1474 + }
  1475 +
  1476 + // ... otherwise handle general case
  1477 + else {
  1478 + sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) *
  1479 + (n_on * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) +
  1480 + n_off * std::log((1.0+alpha) * n_off / (n_on+n_off)));
  1481 + }
  1482 +
  1483 + // Return sigma
  1484 + return sigma;
  1485 +}
  1486 +
  1487 +
  1488 +/***********************************************************************//**
1244 1489 * @brief Write keywords in FITS HDU
1245 1490 *
1246 1491 * @param[in,out] hdu Pointer to FITS HDU.
... ...
src/ctskymap/ctskymap.hpp
  View file @b1668ff
... ... @@ -64,28 +64,35 @@ public:
64 64  
65 65 protected:
66 66 // Protected methods
67   - void init_members(void);
68   - void copy_members(const ctskymap& app);
69   - void free_members(void);
70   - void get_parameters(void);
71   - void setup_maps(void);
72   - void map_exclusions(const GFilename& filename);
73   - void map_exclusions_fits(const GFilename& filename);
74   - void map_exclusions_region(const GFilename& filename);
75   - void map_events(GCTAObservation* obs);
76   - void map_background(GCTAObservation* obs);
77   - void map_background_irf(GCTAObservation* obs);
78   - void map_significance(void);
79   - void map_significance_ring(void);
80   - void compute_ring_values(const int& ipixel,
81   - const GSkyMap& counts,
82   - const GSkyMap& background,
83   - double& non,
84   - double& noff,
85   - double& alpha);
86   - void ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
87   - int& iy1, int& iy2);
88   - void write_hdu_keywords(GFitsHDU* hdu) const;
  67 + void init_members(void);
  68 + void copy_members(const ctskymap& app);
  69 + void free_members(void);
  70 + void get_parameters(void);
  71 + void setup_maps(void);
  72 + void map_exclusions(const GFilename& filename);
  73 + void map_exclusions_fits(const GFilename& filename);
  74 + void map_exclusions_region(const GFilename& filename);
  75 + void map_events(GCTAObservation* obs);
  76 + void map_background(GCTAObservation* obs);
  77 + void map_background_irf(GCTAObservation* obs);
  78 + void map_significance(void);
  79 + void map_significance_ring_direct(void);
  80 + void map_significance_ring_fft(void);
  81 + void compute_ring_values(const int& ipixel,
  82 + const GSkyMap& counts,
  83 + const GSkyMap& background,
  84 + double& non,
  85 + double& noff,
  86 + double& alpha);
  87 + void ring_bounding_box(const int& ipixel, int& ix1, int& ix2,
  88 + int& iy1, int& iy2) const;
  89 + GSkyMap ring_convolve(const GSkyMap& map, const double& rmin,
  90 + const double& rmax) const;
  91 + GNdarray ring_kernel(const double& rmin, const double& rmax) const;
  92 + double sigma_li_ma(const double& n_on,
  93 + const double& n_off,
  94 + const double& alpha) const;
  95 + void write_hdu_keywords(GFitsHDU* hdu) const;
89 96  
90 97 // User parameters
91 98 GFilename m_outmap; //!< Output file name
... ... @@ -95,8 +102,9 @@ protected:
95 102 double m_emax; //!< Maximum energy (TeV)
96 103 std::string m_bkgsubtract; //!< Background subtraction method
97 104 double m_roiradius; //!< Region of interest radius for RING bkg.
98   - double m_inradius; //!< Inner ring radius for RING bkg.
99   - double m_outradius; //!< Outer ring radius for RING bkg.
  105 + double m_inradius; //!< Inner ring radius for RING background
  106 + double m_outradius; //!< Outer ring radius for RING background
  107 + bool m_usefft; //!< Use FFT for RING background
100 108 bool m_publish; //!< Publish sky map?
101 109 GChatter m_chatter; //!< Chattiness
102 110  
... ...
src/ctskymap/ctskymap.par
  View file @b1668ff
... ... @@ -47,6 +47,7 @@ roiradius, r, a, 0.2,,, &quot;Source region radius for estimating on-counts (degree
47 47 inradius, r, a, 0.6,,, "Inner background ring radius (degrees)"
48 48 outradius, r, a, 0.8,,, "Outer background ring radius (degrees)"
49 49 inexclusion, f, h, NONE,,, "Exclusion map or DS9 region file name"
  50 +usefft, b, h, yes,,, "Use FFT for ring computation?"
50 51  
51 52 #
52 53 # Standard parameters
... ...