From 217ed01922444c6dfc06a57bc2e6e7ae0c864da1 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 10 Feb 2018 23:08:34 +0100 Subject: [PATCH] Speed up IRF computation by using GSkyDir::cos_dist() to select relevant pixels. Also implement OpenMP support for multiple observations. --- doc/source/users/reference_manual/ctskymap.rst | 6 +++--- src/ctskymap/ctskymap.cpp | 669 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctskymap/ctskymap.hpp | 29 ++++++++--------------------- src/ctskymap/ctskymap.par | 4 ++-- 4 files changed, 343 insertions(+), 365 deletions(-) diff --git a/doc/source/users/reference_manual/ctskymap.rst b/doc/source/users/reference_manual/ctskymap.rst index 80627f7..ac0296c 100644 --- a/doc/source/users/reference_manual/ctskymap.rst +++ b/doc/source/users/reference_manual/ctskymap.rst @@ -82,13 +82,13 @@ Background subtraction configuration parameters Background subtraction method. ``roiradius [real]`` - Source region radius for estimating on-counts (in degrees). + Source region radius for ``RING`` subtraction (in degrees). ``inradius [real]`` - Inner ring radius for RING subtraction (in degrees). + Inner background ring radius for ``RING`` subtraction (in degrees). ``outradius [real]`` - Outer ring radius for RING subtraction (in degrees). + Outer background ring radius for ``RING`` subtraction (in degrees). ``(inexclusion = NONE) [file]`` Exclusion region file as either a FITS map or DS9 region file. diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index acc9d33..b057f0c 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -35,8 +35,6 @@ #define G_GET_PARAMETERS "ctskymap::get_parameter()" #define G_MAP_EVENTS "ctskymap::map_events(GCTAObservation*)" #define G_MAP_BACKGROUND_IRF "ctskymap::map_background_irf(GCTAObservation*)" -#define G_MAP_BACKGROUND_RING "ctskymap::map_background_ring("\ - "GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ @@ -221,20 +219,44 @@ void ctskymap::run(void) // Get parameters get_parameters(); + // Setup maps + setup_maps(); + // Write input observation container into logger log_observations(NORMAL, m_obs, "Input observation"); // Write header into logger + log_header1(TERSE, gammalib::number("Find unbinned observation", + m_obs.size())); + + // Find all unbinned CTA observations in m_obs + std::vector obs_list(0); + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { + + // Push observation into list + obs_list.push_back(obs); + + // Write message + std::string msg = " Including unbinned "+obs->instrument()+ + " observation"; + log_string(NORMAL, msg); + } + + // Write header into logger log_header1(TERSE, gammalib::number("Map observation", m_obs.size())); // Loop over all unbinned CTA observations in the container - for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; - obs = next_unbinned_observation()) { + #pragma omp parallel for + for (int i = 0; i < obs_list.size(); ++i) { + + // Get pointer to observation + GCTAObservation* obs = obs_list[i]; // Map events into sky map map_events(obs); - // Compute background sky map + // Compute background sky map (only if m_bkgsubtract != "NONE") map_background(obs); // Dispose events to free memory @@ -250,7 +272,8 @@ void ctskymap::run(void) // Subtract background map from counts map m_skymap -= m_bkgmap; - } + + } // endif: background subtraction is selected // Optionally publish sky map if (m_publish) { @@ -292,7 +315,7 @@ void ctskymap::save(void) // If background subtraction is requested then write background map // and significance map to FITS file if (m_bkgsubtract != "NONE") { - + // Write background map into FITS file hdu = m_bkgmap.write(fits); @@ -375,12 +398,8 @@ void ctskymap::publish(const std::string& name) void ctskymap::init_members(void) { // Initialise members - m_skymap.clear(); - m_bkgmap.clear(); - m_sigmap.clear(); - m_exclmap.clear(); - m_alphamap.clear(); - m_onmap.clear(); + m_outmap.clear(); + m_inexclusion.clear(); m_emin = 0.0; m_emax = 0.0; m_bkgsubtract = "NONE"; @@ -390,6 +409,12 @@ void ctskymap::init_members(void) m_publish = false; m_chatter = static_cast(2); + // Initialise maps + m_skymap.clear(); + m_bkgmap.clear(); + m_sigmap.clear(); + m_exclmap.clear(); + // Initialise cache m_solidangle.clear(); m_dirs.clear(); @@ -406,13 +431,9 @@ void ctskymap::init_members(void) ***************************************************************************/ void ctskymap::copy_members(const ctskymap& app) { - // Copy attributes - m_skymap = app.m_skymap; - m_bkgmap = app.m_bkgmap; - m_sigmap = app.m_sigmap; - m_exclmap = app.m_exclmap; - m_alphamap = app.m_alphamap; - m_onmap = app.m_onmap; + // Copy members + m_outmap = app.m_outmap; + m_inexclusion = app.m_inexclusion; m_emin = app.m_emin; m_emax = app.m_emax; m_bkgsubtract = app.m_bkgsubtract; @@ -422,9 +443,15 @@ void ctskymap::copy_members(const ctskymap& app) m_publish = app.m_publish; m_chatter = app.m_chatter; + // Copy maps + m_skymap = app.m_skymap; + m_bkgmap = app.m_bkgmap; + m_sigmap = app.m_sigmap; + m_exclmap = app.m_exclmap; + // Copy cache - m_solidangle = app.m_solidangle; - m_dirs = app.m_dirs; + m_solidangle = app.m_solidangle; + m_dirs = app.m_dirs; // Return return; @@ -456,7 +483,7 @@ void ctskymap::get_parameters(void) // information and do not accept counts cubes. setup_observations(m_obs, false, true, false); - // Create sky map based on task parameters + // Create sky map. Note that the method reads the relevant User parameters. m_skymap = create_map(m_obs); // Get further parameters @@ -468,9 +495,10 @@ void ctskymap::get_parameters(void) if (m_bkgsubtract == "RING") { // Get parameters - m_roiradius = (*this)["roiradius"].real(); - m_inradius = (*this)["inradius"].real(); - m_outradius = (*this)["outradius"].real(); + m_roiradius = (*this)["roiradius"].real(); + m_inradius = (*this)["inradius"].real(); + m_outradius = (*this)["outradius"].real(); + m_inexclusion = (*this)["inexclusion"].filename(); // Make sure that (roiradius < inradius < outradius) if (m_roiradius > m_inradius) { @@ -484,8 +512,8 @@ void ctskymap::get_parameters(void) } // endif: read RING background parameters - // If IRF background subtraction is requested then make sure that the - // CTA observations in the observation container have response information + // If background subtraction is requested then make sure that the CTA + // observations in the observation container have response information if (m_bkgsubtract != "NONE") { set_response(m_obs); } @@ -496,19 +524,38 @@ void ctskymap::get_parameters(void) // Read ahead parameters if (read_ahead()) { - m_outmap = (*this)["outmap"].filename(); + m_outmap = (*this)["outmap"].filename(); } + // Write parameters into logger + log_parameters(TERSE); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Setup maps + * + * Setup maps for sky map generation. + ***************************************************************************/ +void ctskymap::setup_maps(void) +{ + // Clear cache + m_solidangle.clear(); + m_dirs.clear(); + // Create background map and significance map if background subtraction // is requested if (m_bkgsubtract != "NONE") { - // Create backgrond and significance maps + // Create background and significance maps m_bkgmap = create_map(m_obs); m_sigmap = create_map(m_obs); - + // Setup the exclusions map - map_exclusions((*this)["inexclusion"].filename()); + map_exclusions(m_inexclusion); // Cache the pixel solid angles and sky directions m_solidangle.reserve(m_bkgmap.npix()); @@ -518,17 +565,107 @@ void ctskymap::get_parameters(void) m_dirs.push_back(m_bkgmap.inx2dir(i)); } - // If doing a ring background subtraction, generate an alpha map and - // make sure to cache the counts for each observation - if (m_bkgsubtract == "RING") { - m_alphamap = create_map(m_obs); - m_onmap = create_map(m_obs); + } // endif: background subtraction selected + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generates map of pixel exclusions + * + * @param[in] filename Exclusion file name. + * + * Generates a sky map of the pixels that are to be excluded from the + * background estimation. Pixels with values different from 0 will be + * excluded. + ***************************************************************************/ +void ctskymap::map_exclusions(const GFilename& filename) +{ + // Create exlusion map + m_exclmap = create_map(m_obs); + + // Set all pixels to 0 (no pixel excluded) + m_exclmap = 0.0; + + // Make sure the exclusions filename is valid + if (is_valid_filename(filename)) { + + // Fill the exclusions based on the regions supplied + if (filename.is_fits()) { + map_exclusions_fits(filename); + } + else { + map_exclusions_region(filename); } - } // endif: background subtraction selected + } // endif: filename was valid + + // Return + return; +} - // Write parameters into logger - log_parameters(TERSE); + +/***********************************************************************//** + * @brief Fills exclusions map from FITS image + * + * @param[in] filename FITS image file name. + * + * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in + * the exclusion sky map FITS file. + ***************************************************************************/ +void ctskymap::map_exclusions_fits(const GFilename& filename) +{ + // Load the fits image + GSkyMap inmap(filename); + + // Loop through the individual pixels in the exclusion map + for (int i = 0; i < m_exclmap.npix(); ++i) { + + // Get the pixel direction + GSkyDir dir = m_exclmap.pix2dir(i); + + // Check this sky position in the fits map + if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) { + + // Set the pixel to 1 + m_exclmap(i) = 1.0; + + } // endif: pixel,region overlap check + + } // endfor: looped over exclusion map pixels + + // Return + return; +} + + +/***********************************************************************//** + * @brief Fills exclusions map from DS9 region file + * + * @param[in] filename DS9 region file name. + * + * Sets all exclusion map pixels to 1 that are contained in any of the DS9 + * regions. + ***************************************************************************/ +void ctskymap::map_exclusions_region(const GFilename& filename) +{ + // Load the exclusion regions + GSkyRegions regions(filename); + + // Loop through the individual pixels in the exclusion map + for (int i = 0; i < m_exclmap.npix(); ++i) { + + // Get the pixel position + GSkyDir dir = m_exclmap.pix2dir(i); + + // If pixel position overlaps with the regions + if (regions.contains(dir)) { + m_exclmap(i) = 1.0; + } + + } // endfor: looped over exclusion map pixels // Return return; @@ -604,21 +741,26 @@ void ctskymap::map_events(GCTAObservation* obs) } // Fill event in skymap - m_skymap(pixel, 0) += 1.0; + #pragma omp critical(ctskymap_map_events) + m_skymap(pixel) += 1.0; num_in_map++; - + } // endfor: looped over all events // Log binning results - log_value(NORMAL, "Events in list", obs->events()->size()); - log_value(NORMAL, "Events in map", num_in_map); - log_value(NORMAL, "Events outside RoI", num_outside_roi); - log_value(NORMAL, "Events outside map area", num_outside_map); - log_value(NORMAL, "Events outside energies", num_outside_erange); - - // Write sky map into header - log_header1(EXPLICIT, "Sky map"); - log_string(EXPLICIT, m_skymap.print(m_chatter)); + #pragma omp critical(ctskymap_map_events) + { + log_header3(TERSE, get_obs_header(obs)); + log_value(NORMAL, "Events in list", obs->events()->size()); + log_value(NORMAL, "Events in map", num_in_map); + log_value(NORMAL, "Events outside RoI", num_outside_roi); + log_value(NORMAL, "Events outside map area", num_outside_map); + log_value(NORMAL, "Events outside energies", num_outside_erange); + + // Write sky map into header + log_header1(EXPLICIT, "Sky map"); + log_string(EXPLICIT, m_skymap.print(m_chatter)); + } // Return return; @@ -646,7 +788,7 @@ void ctskymap::map_background(GCTAObservation* obs) map_background_irf(obs); } else if (m_bkgsubtract == "RING") { - map_background_ring(obs); + map_background_irf(obs); } // Return @@ -670,7 +812,8 @@ void ctskymap::map_background(GCTAObservation* obs) void ctskymap::map_background_irf(GCTAObservation* obs) { // Get IRF response - const GCTAResponseIrf* rsp = dynamic_cast(obs->response()); + const GCTAResponseIrf* rsp = dynamic_cast + (obs->response()); // Throw an exception if observation has no instrument response function if (rsp == NULL) { @@ -697,8 +840,16 @@ void ctskymap::map_background_irf(GCTAObservation* obs) GEnergy emin(m_emin, "TeV"); GEnergy emax(m_emax, "TeV"); - // Extract region of interest from observation - GCTARoi roi = obs->roi(); + // Extract region of interest center and radius from observation and + // compute the cosine of the RoI radius. This allows to mimimize as + // much as possible the trigonometric computations. Note that the cosine + // of an angle is maximal for an angle of zero. This explains later why + // we chose ">=" to test whether an angle is smaller than a given radius. + GSkyDir roi_centre = obs->roi().centre().dir(); + double cos_roi_radius = std::cos(obs->roi().radius() * gammalib::deg2rad); + + // Extract exposure + double exposure = obs->livetime(); // Initialise statistics double total = 0.0; @@ -707,25 +858,25 @@ void ctskymap::map_background_irf(GCTAObservation* obs) for (int i = 0; i < m_bkgmap.npix(); ++i) { // Get sky direction of pixel - GSkyDir skydir = m_bkgmap.inx2dir(i); - - // Convert sky direction in instrument direction - GCTAInstDir instdir = obs->pointing().instdir(skydir); + GSkyDir& skydir = m_dirs[i]; - // If RoI is valid and instrument direction is not within RoI then - // skip pixel - if (roi.is_valid() && !roi.contains(instdir)) { + // Skip pixel if it is not within the RoI + if (roi_centre.cos_dist(skydir) < cos_roi_radius) { continue; } + // Convert sky direction to instrument direction + GCTAInstDir instdir = obs->pointing().instdir(skydir); + // Compute background value double value = bkg->rate_ebin(instdir, emin, emax); // Multiply background rate with livetime and solid angle - value *= obs->livetime() * m_solidangle[i]; + value *= exposure * m_solidangle[i]; // Add number of background events to background map - m_bkgmap(i,0) += value; + #pragma omp critical(ctskymap_map_background_irf) + m_bkgmap(i) += value; // Update total number of background events total += value; @@ -733,7 +884,10 @@ void ctskymap::map_background_irf(GCTAObservation* obs) } // endfor: looped over background map pixels // Log background subtraction results - log_value(NORMAL, "Events in background", int(total+0.5)); + #pragma omp critical(ctskymap_map_events) + { + log_value(NORMAL, "Events in background", int(total+0.5)); + } // Return return; @@ -741,196 +895,122 @@ void ctskymap::map_background_irf(GCTAObservation* obs) /***********************************************************************//** - * @brief Estimates the background in sky map based on the ring background - * method - * - * @param[in] obs CTA observation. - * - * @exception GException::invalid_value - * No response information available for observation. - * No background template available in instrument response function. + * @brief Computes the bin-by-bin significance if background subtraction + * requested * - * Estimates the background in the sky map by summing the events within a - * ring centered at a given pixel's position. The pixels in the ring are - * weighted also by the background IRF value. + * Method computes the bin-by bin significance if a background subtraction + * method is specified. If method is "IRF", Poisson statistics (in the + * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance + * is computed for each bin. ***************************************************************************/ -void ctskymap::map_background_ring(GCTAObservation* obs) +void ctskymap::map_significance(void) { - // Get IRF response (to scale background counts) - const GCTAResponseIrf* rsp = dynamic_cast(obs->response()); - - // Throw an exception if observation has no instrument response function - if (rsp == NULL) { - std::string msg = "No response information available for "+ - get_obs_header(obs)+" to compute IRF background. " - "Please specify response information or use " - "another background subtraction method."; - throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg); + // Compute significance from "RING" method (Li & Ma eq. 17) + if (m_bkgsubtract == "RING") { + map_significance_ring(); } - // Get IRF background template - const GCTABackground* bkg = rsp->background(); - - // Throw an exception if observation has no IRF background template - if (bkg == NULL) { - std::string msg = "No IRF background template found in instrument " - "response function for "+ - get_obs_header(obs)+". Please specify an instrument " - "response function containing a background template."; - throw GException::invalid_value(G_MAP_BACKGROUND_RING, msg); + // Compute significance using Poisson statistics in the Gaussian limit + else { + m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap); } - // Set minimum and maximum energy as GEnergy instances - GEnergy emin(m_emin, "TeV"); - GEnergy emax(m_emax, "TeV"); - - // Extract region of interest from observation - GCTARoi roi = obs->roi(); - - // Extract exposure - double exposure = obs->livetime(); - - // Loop over all map pixels - for (int i = 0; i < m_bkgmap.npix(); ++i) { - - // Get sky direction of pixel - GSkyDir& skydir = m_dirs[i]; - - // Convert sky direction in instrument direction - GCTAInstDir instdir = obs->pointing().instdir(skydir); - - // If RoI is valid and instrument direction is not within RoI then - // skip pixel - if (roi.is_valid() && !roi.contains(instdir)) { - continue; - } - - // Compute background value - m_alphamap(i,0) += bkg->rate_ebin(instdir, emin, emax) * - m_solidangle[i] * exposure; - - // Store the counts in this bin - m_onmap(i,0) += m_skymap(i,0); - - } // endfor: Loop for caching bkg IRF sensitivity - - // Zero out the counts map to prepare it for the next observation - m_skymap = 0.0; - // Return return; } /***********************************************************************//** - * @brief Write keywords in FITS HDU - * - * @param[in,out] hdu Pointer to FITS HDU. + * @brief Computes the bin-by-bin significance for RING background * - * Writes keywords in FITS HDU. + * Computes the a Li & Ma significance for each sky map bin and replaces + * the sky and background maps by the On- and Off-count maps. ***************************************************************************/ -void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const +void ctskymap::map_significance_ring(void) { - // Continue only if pointer is valid - if (hdu != NULL) { + // Log message about what is being done + log_header1(NORMAL, "Computing Ring background map"); + log_value(NORMAL, "Total pixels to process", m_skymap.npix()); - // Set keywords - hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method"); - hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary"); - hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary"); - hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX"); - - } // endif: pointer was valid + // Initialise On- and Off-count maps + GSkyMap onmap = m_skymap; + GSkyMap offmap = m_skymap; + onmap = 0.0; + offmap = 0.0; - // Return - return; -} + // Initialise statistics + int num_bad_alpha = 0; + + // Set progress logging frequency. By default the user is updated on the + // progress when another 10% of pixels is processed, however, if there + // are many pixels the user is updated each 100000 pixels. + int n_logpix = m_skymap.npix() / 10; + if (n_logpix > 100000) { + n_logpix = 100000; + } + // Loop over sky map pixels + #pragma omp parallel for + for (int i = 0; i < m_skymap.npix(); ++i) { -/***********************************************************************//** - * @brief Computes the bin-by-bin significance if background subtraction - * requested - * - * Method computes the bin-by bin significance if a background subtraction - * method is specified. If method is "IRF", Poisson statistics (in the - * Gaussian limit) are assumed. If method is "RING" a Li & Ma significance - * is computed for each bin (as well as a map of the computed alpha values). - ***************************************************************************/ -void ctskymap::map_significance(void) -{ - // Compute significance from "RING" method (Li & Ma eq. 17) - if (m_bkgsubtract == "RING") { - - // Log message about what is being done - log_header1(NORMAL, "Computing Ring background map"); - log_value(NORMAL, "Total pixels to process", m_onmap.npix()); + // Initialise the on/off counts and alpha for this bin + double n_on = 0.0; + double n_off = 0.0; + double alpha = 0.0; - // Store the number of bins with an inappropriate alpha parameter - int num_bad_alpha = 0; + // Update user about progress + if (i % n_logpix == 0) { + #pragma omp critical(ctskymap_map_significance_ring) + log_value(NORMAL, "Pixels remaining", m_skymap.npix()-i); + } - // Loop through each bin in the on-counts map - for (int i = 0; i < m_onmap.npix(); ++i) { + // Get bin coordinates + GSkyDir& skydir = m_dirs[i]; - // Initialise the on/off counts and alpha for this bin - double n_on = 0.0; - double n_off = 0.0; - double alpha = 0.0; + // Compute the alpha and counts for this bin + compute_ring_values(m_skymap, m_bkgmap, skydir, n_on, n_off, alpha); - // Since this can take a long time, keep the user updated on - // the progress when another 10% of pixels is processed - if (i % (m_onmap.npix()/10) == 0) { - log_value(NORMAL, "Pixels remaining", m_onmap.npix()-i); - } + // Store the On and alpha-weighted Off-counts + onmap(i) = n_on; + offmap(i) = alpha * n_off; - // Get bin coordinates - GSkyDir& skydir = m_dirs[i]; + // If alpha is zero then increment the bad alpha counter + if (alpha == 0.0) { + num_bad_alpha++; + } - // Compute the alpha and counts for this bin - compute_ring_values(m_onmap, m_alphamap, skydir, - n_on, n_off, alpha); + // ... otherwise store the results + else { - // If alpha is zero then increment the bad alpha counter - if (alpha == 0.0) { - num_bad_alpha++; + // Compute and store significance (Li & Ma eq. 17) + double sigma; + if (n_on == 0.0) { + sigma = -2.0 * n_off * std::log(1.0+alpha); } - - // ... otherwise store the results else { + sigma = (n_on < (alpha*n_off) ? -2.0 : 2.0) * + (n_on * std::log((1.0+alpha) * n_on / (alpha * (n_on+n_off))) + + n_off * std::log((1.0+alpha) * n_off / (n_on+n_off))); + } + #pragma omp critical(ctskymap_map_significance_ring) + m_sigmap(i) = sigma; - // Store the on-counts and alpha-weighted off-counts - m_skymap(i,0) = n_on; - m_bkgmap(i,0) = alpha * n_off; - - // Compute and store significance (Li & Ma eq. 17) - if (n_on == 0.0) { - m_sigmap(i,0) = -2.0 * n_off * std::log(1.0+alpha); - } - else { - m_sigmap(i,0) = (n_on < (alpha*n_off) ? -2.0 : 2.0) * - (n_on * std::log((1.0+alpha) * - n_on / (alpha * (n_on+n_off))) + - n_off * std::log((1.0+alpha) * - n_off / (n_on+n_off))); - } - - } // endelse: alpha was non-zero - - } // endfor: looped over all pixels + } // endelse: alpha was non-zero - // Log the number of bad-alpha bins - log_value(NORMAL, "Bins with alpha=0", num_bad_alpha); + } // endfor: looped over all pixels - // Now take the square root. Since some bins can be negative, first - // take the absolute value of the map, square root that, then multiply - // each bin by its sign to preserve the +/- significance. - m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap)); + // Log the number of bad-alpha bins + log_value(NORMAL, "Bins with alpha=0", num_bad_alpha); - } // endif: ring background selected + // Now take the square root. Since some bins can be negative, first + // take the absolute value of the map, square root that, then multiply + // each bin by its sign to preserve the +/- significance. + m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap)); - // Compute significance using Poisson statistics in the Gaussian limit - else { - m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap); - } + // Terminate by storing the on-counts map as sky map and the off-counts + // map as background map + m_skymap = onmap; + m_bkgmap = offmap; // Return return; @@ -940,18 +1020,18 @@ void ctskymap::map_significance(void) /***********************************************************************//** * @brief Computes Non, Noff and alpha for a counts map and sensitivity map * - * @param[in] counts Counts map - * @param[in] sensitivity Sensitivity map - * @param[in] position Position of the ring & roi centers - * @param[out] non Returned estimate of ON counts - * @param[out] noff Returned estimate of OFF counts - * @param[out] alpha Returned estimate of alpha + * @param[in] counts Counts map + * @param[in] background Background map + * @param[in] position Position of the ring & roi centers + * @param[out] non Returned estimate of ON counts + * @param[out] noff Returned estimate of OFF counts + * @param[out] alpha Returned estimate of alpha * * This method computes the Non counts Noff values at a given position. It * also computes the alpha parameter from the passed sensitivity map. ***************************************************************************/ void ctskymap::compute_ring_values(const GSkyMap& counts, - const GSkyMap& sensitivity, + const GSkyMap& background, const GSkyDir& position, double& non, double& noff, @@ -984,13 +1064,13 @@ void ctskymap::compute_ring_values(const GSkyMap& counts, if (position.cos_dist(skydir) >= cos_outradius) { // Check if pixel is inside the background region - if ((m_exclmap(j,0) == 0.0) && position.cos_dist(skydir) < cos_inradius) { + if ((m_exclmap(j) == 0.0) && position.cos_dist(skydir) < cos_inradius) { // Update n_off - noff += counts(j,0); + noff += counts(j); // Update alpha_off - alpha_off += sensitivity(j,0); + alpha_off += background(j); } @@ -998,10 +1078,10 @@ void ctskymap::compute_ring_values(const GSkyMap& counts, else if (position.cos_dist(skydir) >= cos_roiradius) { // Update n_on for significance computation - non += counts(j,0); + non += counts(j); // Update alpha_on - alpha_on += sensitivity(j,0); + alpha_on += background(j); } // endif: source and background region check @@ -1026,121 +1106,32 @@ void ctskymap::compute_ring_values(const GSkyMap& counts, /***********************************************************************//** - * @brief Generates map of pixel exclusions - * - * @param[in] filename Exclusion file name. - * - * Generates a sky map of the pixels that are to be excluded from the - * background estimation. Pixels with values different from 0 will be - * excluded. - ***************************************************************************/ -void ctskymap::map_exclusions(const GFilename& filename) -{ - // Create exlusion map - m_exclmap = create_map(m_obs); - - // Set all pixels to 0 (no pixel excluded) - m_exclmap = 0.0; - - // Make sure the exclusions filename is valid - if (is_valid_filename(filename)) { - - // Fill the exclusions based on the regions supplied - if (filename.is_fits()) { - map_exclusions_fits(filename); - } - else { - map_exclusions_reg(filename); - } - - } // endif: filename was valid - - // Return - return; -} - - -/***********************************************************************//** - * @brief Fills exclusions map from FITS image - * - * @param[in] filename FITS image file name. - * - * Sets all exclusion map pixels to 1 that correspond to non-zero pixels in - * the exclusion sky map FITS file. - ***************************************************************************/ -void ctskymap::map_exclusions_fits(const GFilename& filename) -{ - // Load the fits image - GSkyMap inmap(filename); - - // Loop through the individual pixels in the exclusion map - for (int i = 0; i < m_exclmap.npix(); ++i) { - - // Get the pixel direction - GSkyDir dir = m_exclmap.pix2dir(i); - - // Check this sky position in the fits map - if (inmap.contains(dir) && (inmap(inmap.dir2inx(dir)) != 0.0)) { - - // Set the pixel to 1 - m_exclmap(i) = 1.0; - - } // endif: pixel,region overlap check - - } // endfor: looped over exclusion map pixels - - // Return - return; -} - - -/***********************************************************************//** - * @brief Fills exclusions map from DS9 region file + * @brief Write keywords in FITS HDU * - * @param[in] filename DS9 region file name. + * @param[in,out] hdu Pointer to FITS HDU. * - * Sets all exclusion map pixels to 1 that are contained in any of the DS9 - * regions. + * Writes keywords in FITS HDU. ***************************************************************************/ -void ctskymap::map_exclusions_reg(const GFilename& filename) +void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const { - // Load the exclusion regions - GSkyRegions regions(filename); - - // Loop through the individual pixels in the exclusion map - for (int i = 0; i < m_exclmap.npix(); ++i) { + // Continue only if pointer is valid + if (hdu != NULL) { - // Get the pixel position - GSkyDir dir = m_exclmap.pix2dir(i); + // Set keywords + hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method"); + hdu->card("E_MIN", m_emin, "[TeV] Lower energy boundary"); + hdu->card("E_MAX", m_emax, "[TeV] Upper energy boundary"); + hdu->card("EUNIT", "TeV", "Units for E_MIN and E_MAX"); - // If pixel position overlaps with the regions - if (regions.contains(dir)) { - m_exclmap(i) = 1.0; + // Write RING background keywords + if (m_bkgsubtract == "RING") { + hdu->card("ROIRAD", m_roiradius, "[deg] Source region radius"); + hdu->card("INRAD", m_inradius, "[deg] Inner background ring radius"); + hdu->card("OUTRAD", m_outradius, "[deg] Outer background ring radius"); } - } // endfor: looped over exclusion map pixels + } // endif: pointer was valid // Return return; } - - -/***********************************************************************//** - * @brief Estimates the background in sky map based on IRF template - * - * @param[in] lnE Natural logarithm of energy in MeV. - ***************************************************************************/ -double ctskymap::irf_kern::eval(const double& lnE) -{ - // Get log10 of energy in TeV - double logE = lnE * gammalib::inv_ln10 - 6; - - // Get function value - double value = (*m_bgd)(logE, m_dir->detx(), m_dir->dety()); - - // Correct for variable substitution - value *= std::exp(lnE); - - // Return value - return value; -} diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index c3a1300..4fe1ba2 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -68,35 +68,24 @@ protected: void copy_members(const ctskymap& app); void free_members(void); void get_parameters(void); + void setup_maps(void); + void map_exclusions(const GFilename& filename); + void map_exclusions_fits(const GFilename& filename); + void map_exclusions_region(const GFilename& filename); void map_events(GCTAObservation* obs); void map_background(GCTAObservation* obs); void map_background_irf(GCTAObservation* obs); - void map_background_ring(GCTAObservation* obs); void map_significance(void); - void map_exclusions(const GFilename& filename); - void map_exclusions_fits(const GFilename& filename); - void map_exclusions_reg(const GFilename& filename); - void compute_ring_values(const GSkyMap& counts, - const GSkyMap& sensitivity, + void map_significance_ring(void); + void compute_ring_values(const GSkyMap& counts, + const GSkyMap& background, const GSkyDir& position, double& non, double& noff, double& alpha); void write_hdu_keywords(GFitsHDU* hdu) const; - // Background integration kernel - class irf_kern : public GFunction { - public: - irf_kern(const GCTABackground* bgd, - const GCTAInstDir* dir) : - m_bgd(bgd), - m_dir(dir) { } - double eval(const double& lnE); - protected: - const GCTABackground* m_bgd; //!< Pointer to background - const GCTAInstDir* m_dir; //!< Pointer to instrument direction - }; - // User parameters GFilename m_outmap; //!< Output file name + GFilename m_inexclusion; //!< Exclusion map file name double m_emin; //!< Minimum energy (TeV) double m_emax; //!< Maximum energy (TeV) std::string m_bkgsubtract; //!< Background subtraction method @@ -111,8 +100,6 @@ protected: GSkyMap m_bkgmap; //!< Background map GSkyMap m_sigmap; //!< Significance map GSkyMap m_exclmap; //!< Exclusion map for RING background - GSkyMap m_alphamap; //!< Alpha map values for RING background - GSkyMap m_onmap; //!< On counts map for RING background // Caching variables to prevent multiple computations of the same thing std::vector m_solidangle; //!< Cached pixel solid angles diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par index 86e8398..88b6241 100644 --- a/src/ctskymap/ctskymap.par +++ b/src/ctskymap/ctskymap.par @@ -44,8 +44,8 @@ yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC #================================== bkgsubtract, s, a, NONE, NONE|IRF|RING,, "Background subtraction method" roiradius, r, a, 0.2,,, "Source region radius for estimating on-counts (degrees)" -inradius, r, a, 0.6,,, "Inner ring radius (degrees)" -outradius, r, a, 0.8,,, "Outer ring radius (degrees)" +inradius, r, a, 0.6,,, "Inner background ring radius (degrees)" +outradius, r, a, 0.8,,, "Outer background ring radius (degrees)" inexclusion, f, h, NONE,,, "Exclusion map or DS9 region file name" # -- libgit2 0.21.4