From 997ffff511b58a071bdf5d64ab3bd5608d813719 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 13 Feb 2018 16:00:35 +0100 Subject: [PATCH] Allow specifying of input sky map in ctskymap (#2311) --- ChangeLog | 3 ++- NEWS | 19 +++++++++++++++++-- doc/source/admin/release_history/1.6.rst | 2 ++ doc/source/users/reference_manual/ctskymap.rst | 18 ++++++++++++++++++ pyext/ctobservation.i | 1 + src/ctskymap/ctskymap.cpp | 683 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctskymap/ctskymap.hpp | 27 +++++++++++++++++---------- src/ctskymap/ctskymap.par | 5 ++++- src/support/ctobservation.cpp | 42 ++++++++++++++++++++++++++++++++++++++++++ src/support/ctobservation.hpp | 3 ++- src/support/ctool.cpp | 16 ++++++++++------ 11 files changed, 519 insertions(+), 300 deletions(-) diff --git a/ChangeLog b/ChangeLog index 00bbb77..b534a1c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2018-02-10 +2018-02-13 * Version 1.6.0 released ======================== + Allow specifying of input sky map in ctskymap (#2311) Accelerate IRF and RING background methods in ctskymap (#2309) diff --git a/NEWS b/NEWS index 18177e3..2841905 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.6.0 -10 February 2018 +13 February 2018 Introduction @@ -106,7 +106,22 @@ None ctskymap - CTA sky mapping -------------------------- -Accelerate the IRF and RING background methods (#2309). +Add the "inmap", "iterations" and "threshold" parameters. The hidden +"inmap" parameter allows specifying a sky map file containing a "COUNTS" +and "ACCEPTANCE" extension which avoids looping over event lists in an +observation container. The "iterations" parameter allows the iterative +determination of an exclusion map from the resulting significance map for +the RING background method. The exclusion map will flag all pixels with +significances above the "threshold" value so that they will be excluded +in the subsequent RING background determination. Since this will impact +the pixel significances, the procedure needs to be repeated a couple of +times to stabilise (#2311). + +Accelerate the IRF and RING background methods. For the RING background +method the computation is either performed using a FFT or using a direct +computation. The FFT computation is faster, but assumed Euclidian +distances for the sky map grid, the direct computation is slower, but +accurately computes the distances between sky map pixels (#2309). cttsmap - Generation of Test Statistic map diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 96c92a2..6bd9a20 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -22,6 +22,8 @@ Bug fixes Improvements ------------ +* [`2311 `_] - + Allow specifying of input skymap to :ref:`ctskymap` tool * [`2309 `_] - Accelerate `IRF` and `RING` background methods in :ref:`ctskymap` diff --git a/doc/source/users/reference_manual/ctskymap.rst b/doc/source/users/reference_manual/ctskymap.rst index d0bfc2b..6af9c7e 100644 --- a/doc/source/users/reference_manual/ctskymap.rst +++ b/doc/source/users/reference_manual/ctskymap.rst @@ -41,6 +41,11 @@ General parameters Instrument response function (only required for IRF background subtraction if no response information is provided by ``inobs``). +``(inmap = NONE) [file]`` + Input sky map file containing a ``COUNTS`` and an ``ACCEPTANCE`` extension. + Such extensions are produce by a previous run of the :ref:`ctskymap` tool with + the ``IRF`` or ``RING`` background method. + ``outmap [file]`` Output sky map file. @@ -90,6 +95,19 @@ Background subtraction configuration parameters ``outradius [real]`` Outer background ring radius for ``RING`` subtraction (in degrees). +``iterations [integer]`` + Number of iterations for the automatic computation of exclusion regions. + Sky map pixels with detection significance above the value specified by the + ``threshold`` parameter will be defined as exclusion regions, and the ``RING`` + background will be recomputed based on these exclusion regions. Since this + will likely change the significance of the pixels, the procedure can be + iteratively repeated, and the ``iterations`` parameter specifies how often + the procedure should be repeated. + +``threshold [real]`` + Significance threshold above which pixels should be excluded for the + background computation in the ``RING`` background method. + ``(inexclusion = NONE) [file]`` Exclusion region file as either a FITS map or DS9 region file. diff --git a/pyext/ctobservation.i b/pyext/ctobservation.i index 9270a21..38f221f 100644 --- a/pyext/ctobservation.i +++ b/pyext/ctobservation.i @@ -94,6 +94,7 @@ public: // Make methods private in Python by prepending an underscore %rename(_first_unbinned_observation) first_unbinned_observation; %rename(_next_unbinned_observation) next_unbinned_observation; + %rename(_read_ogip_keywords) read_ogip_keywords; %rename(_write_ogip_keywords) write_ogip_keywords; %rename(_set_obs_statistic) set_obs_statistic; %rename(_save_events_fits) save_events_fits; diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index ff17af5..a75fe29 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -35,6 +35,8 @@ #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_RING_BOUNDING_BOX "ctskymap::ring_bounding_box(int&, int&, int&, "\ + "int&, int&)" #define G_RING_KERNEL "ctskymap::ring_kernel(double&, double&)" /* __ Debug definitions __________________________________________________ */ @@ -223,58 +225,14 @@ void ctskymap::run(void) // 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); + // If no input map exists then fill counts and acceptance maps from + // observation container + if (!m_has_inmap) { + fill_maps(); } - // Write header into logger - log_header1(TERSE, gammalib::number("Map observation", m_obs.size())); - - // Loop over all unbinned CTA observations in the container - #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 (only if m_bkgsubtract != "NONE") - map_background(obs); - - // Dispose events to free memory - obs->dispose_events(); - - } // endfor: looped over observations - - // If background subtraction is selected then compute significance map - if (m_bkgsubtract != "NONE") { - - // Compute the significance map - map_significance(); - - // Subtract background map from counts map - m_skymap -= m_bkgmap; - - } // endif: background subtraction is selected + // Compute maps + compute_maps(); // Optionally publish sky map if (m_publish) { @@ -289,8 +247,8 @@ void ctskymap::run(void) /***********************************************************************//** * @brief Save sky map * - * Saves the sky map into a FITS file. The FITS file name is specified by the - * @p outname parameter. + * Saves the sky map into a FITS file. The FITS file name is specified by + * the @p outname parameter. ***************************************************************************/ void ctskymap::save(void) { @@ -307,39 +265,31 @@ void ctskymap::save(void) GFits fits; // Write sky map into FITS file - GFitsHDU* hdu = m_skymap.write(fits); + write_map(fits, m_skymap, "SKYMAP"); - // Write keywords into sky map extension - write_ogip_keywords(hdu); - write_hdu_keywords(hdu); - - // If background subtraction is requested then write background map - // and significance map to FITS file + // If background subtraction is requested then write the background map, + // the significance map, the counts map and the acceptance map to the + // FITS file if (m_bkgsubtract != "NONE") { // Write background map into FITS file - hdu = m_bkgmap.write(fits); - - // Set background map extension name - if (hdu != NULL) { - hdu->extname("BACKGROUND"); - } + write_map(fits, m_bkgmap, "BACKGROUND"); - // Write keywords into background extension - write_ogip_keywords(hdu); - write_hdu_keywords(hdu); // Write significance map into FITS file - hdu = m_sigmap.write(fits); + write_map(fits, m_sigmap, "SIGNIFICANCE"); - // Set significance map extension name - if (hdu != NULL) { - hdu->extname("SIGNIFICANCE"); - } + // Write counts map into FITS file + write_map(fits, m_counts, "COUNTS"); - // Write keywords into significance extension - write_ogip_keywords(hdu); - write_hdu_keywords(hdu); + // Write acceptance map into FITS file + write_map(fits, m_acceptance, "ACCEPTANCE"); + + // If RING background subtraction is requested then write also + // the exclusion map to the FITS file + if (m_bkgsubtract == "RING") { + write_map(fits, m_exclmap, "EXCLUSION"); + } } // endif: background subtraction was requested @@ -401,21 +351,25 @@ void ctskymap::init_members(void) // Initialise User parameters m_outmap.clear(); m_inexclusion.clear(); - m_binsz = 0.0; m_emin = 0.0; m_emax = 0.0; m_bkgsubtract = "NONE"; m_roiradius = 0.0; m_inradius = 0.0; m_outradius = 0.0; + m_iterations = 0; + m_threshold = 5.0; m_usefft = true; m_publish = false; m_chatter = static_cast(2); // Initialise members + m_has_inmap = false; m_skymap.clear(); m_bkgmap.clear(); m_sigmap.clear(); + m_counts.clear(); + m_acceptance.clear(); m_exclmap.clear(); // Initialise cache @@ -440,22 +394,26 @@ void ctskymap::copy_members(const ctskymap& app) // Copy User parameters m_outmap = app.m_outmap; m_inexclusion = app.m_inexclusion; - m_binsz = app.m_binsz; m_emin = app.m_emin; m_emax = app.m_emax; m_bkgsubtract = app.m_bkgsubtract; m_roiradius = app.m_roiradius; m_inradius = app.m_inradius; m_outradius = app.m_outradius; + m_iterations = app.m_iterations; + m_threshold = app.m_threshold; m_usefft = app.m_usefft; m_publish = app.m_publish; m_chatter = app.m_chatter; // Copy members - m_skymap = app.m_skymap; - m_bkgmap = app.m_bkgmap; - m_sigmap = app.m_sigmap; - m_exclmap = app.m_exclmap; + m_has_inmap = app.m_has_inmap; + m_skymap = app.m_skymap; + m_bkgmap = app.m_bkgmap; + m_sigmap = app.m_sigmap; + m_counts = app.m_counts; + m_acceptance = app.m_acceptance; + m_exclmap = app.m_exclmap; // Copy cache m_solidangle = app.m_solidangle; @@ -490,17 +448,61 @@ void ctskymap::free_members(void) ***************************************************************************/ void ctskymap::get_parameters(void) { - // Setup observations from "inobs" parameter. Do not request response - // information and do not accept counts cubes. - setup_observations(m_obs, false, true, false); + // So far we have no valid input map + m_has_inmap = false; + + // If an input map was specified then extract try loading the COUNTS + // and ACCEPTANCE extensions from that file + GFilename inmap = (*this)["inmap"].filename(); + if (inmap != "NONE") { + GFits fits(inmap); + if (fits.contains("COUNTS") && fits.contains("ACCEPTANCE")) { + + // Get HDUs + GFitsHDU &counts = *fits["COUNTS"]; + GFitsHDU &acceptance = *fits["ACCEPTANCE"]; - // Create sky map. Note that the method reads the relevant User parameters. - m_skymap = create_map(m_obs); + // Read counts and acceptance map + m_counts.read(counts); + m_acceptance.read(acceptance); - // Get further parameters - m_binsz = (*this)["binsz"].real(); - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); + // Extract OGIP parameters + read_ogip_keywords(&counts); + + // Extract energy range + m_emin = (counts.has_card("E_MIN")) ? counts.real("E_MIN") : 0.0; + m_emax = (counts.has_card("E_MAX")) ? counts.real("E_MAX") : 0.0; + + // Signal availability of sky map + m_has_inmap = true; + } + else { + std::string msg = "Input sky map \""+inmap.url()+"\" does not " + "contain \"COUNTS\" and \"ACCEPTANCE\" " + "extensions."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + fits.close(); + } + + // If no valid input map exists then request observations + if (!m_has_inmap) { + + // Setup observations from "inobs" parameter. Do not request response + // information and do not accept counts cubes. + setup_observations(m_obs, false, true, false); + + // Create counts and acceptance maps from the User parameters + m_counts = create_map(m_obs); + m_acceptance = m_counts; + + // Get further parameters + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + + } // endif: no valid input map existed + + // Get background method m_bkgsubtract = (*this)["bkgsubtract"].string(); // Get RING background parameters @@ -510,16 +512,20 @@ void ctskymap::get_parameters(void) m_roiradius = (*this)["roiradius"].real(); m_inradius = (*this)["inradius"].real(); m_outradius = (*this)["outradius"].real(); + m_iterations = (*this)["iterations"].integer(); + if (m_iterations > 0) { + m_threshold = (*this)["threshold"].real(); + } m_inexclusion = (*this)["inexclusion"].filename(); m_usefft = (*this)["usefft"].boolean(); // Make sure that (roiradius < inradius < outradius) if (m_roiradius > m_inradius) { - std::string msg("'roiradius' must be smaller than 'inradius'"); + std::string msg("\"roiradius\" must be smaller than \"inradius\""); throw GException::invalid_value(G_GET_PARAMETERS, msg); } else if (m_inradius > m_outradius) { - std::string msg("'inradius' must be smaller than 'outradius'"); + std::string msg("\"inradius\" must be smaller than \"outradius\""); throw GException::invalid_value(G_GET_PARAMETERS, msg); } @@ -527,7 +533,7 @@ void ctskymap::get_parameters(void) // If background subtraction is requested then make sure that the CTA // observations in the observation container have response information - if (m_bkgsubtract != "NONE") { + if (m_bkgsubtract != "NONE" && !m_has_inmap) { set_response(m_obs); } @@ -563,19 +569,15 @@ void ctskymap::setup_maps(void) // is requested if (m_bkgsubtract != "NONE") { - // Create background and significance maps - m_bkgmap = create_map(m_obs); - m_sigmap = create_map(m_obs); - // Setup the exclusions map - map_exclusions(m_inexclusion); + setup_exclusion_map(m_inexclusion); // Cache the pixel solid angles and sky directions - m_solidangle.reserve(m_bkgmap.npix()); - m_dirs.reserve(m_bkgmap.npix()); - for (int i = 0; i < m_bkgmap.npix(); ++i) { - m_solidangle.push_back(m_bkgmap.solidangle(i)); - m_dirs.push_back(m_bkgmap.inx2dir(i)); + m_solidangle.reserve(m_counts.npix()); + m_dirs.reserve(m_counts.npix()); + for (int i = 0; i < m_counts.npix(); ++i) { + m_solidangle.push_back(m_counts.solidangle(i)); + m_dirs.push_back(m_counts.inx2dir(i)); } } // endif: background subtraction selected @@ -601,10 +603,10 @@ void ctskymap::setup_maps(void) * background estimation. Pixels with values different from 0 will be * excluded. ***************************************************************************/ -void ctskymap::map_exclusions(const GFilename& filename) +void ctskymap::setup_exclusion_map(const GFilename& filename) { // Create exlusion map - m_exclmap = create_map(m_obs); + m_exclmap = m_counts; // Set all pixels to 0 (no pixel excluded) m_exclmap = 0.0; @@ -614,10 +616,10 @@ void ctskymap::map_exclusions(const GFilename& filename) // Fill the exclusions based on the regions supplied if (filename.is_fits()) { - map_exclusions_fits(filename); + setup_exclusion_map_fits(filename); } else { - map_exclusions_region(filename); + setup_exclusion_map_region(filename); } } // endif: filename was valid @@ -635,7 +637,7 @@ void ctskymap::map_exclusions(const GFilename& filename) * 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) +void ctskymap::setup_exclusion_map_fits(const GFilename& filename) { // Load the fits image GSkyMap inmap(filename); @@ -669,7 +671,7 @@ void ctskymap::map_exclusions_fits(const GFilename& filename) * Sets all exclusion map pixels to 1 that are contained in any of the DS9 * regions. ***************************************************************************/ -void ctskymap::map_exclusions_region(const GFilename& filename) +void ctskymap::setup_exclusion_map_region(const GFilename& filename) { // Load the exclusion regions GSkyRegions regions(filename); @@ -693,16 +695,75 @@ void ctskymap::map_exclusions_region(const GFilename& filename) /***********************************************************************//** - * @brief Map events into a sky map + * @brief Fill maps from observation container + * + * Fill counts and optionally acceptance maps from data in the observation + * container. + ***************************************************************************/ +void ctskymap::fill_maps(void) +{ + // 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("Fill map from observation", + m_obs.size())); + + // Loop over all unbinned CTA observations in the container + #pragma omp parallel for + for (int i = 0; i < obs_list.size(); ++i) { + + // Get pointer to observation + GCTAObservation* obs = obs_list[i]; + + // Fill events into counts map + fill_maps_counts(obs); + + // Optionally compute acceptance sky map + if (m_bkgsubtract != "NONE") { + fill_maps_acceptance(obs); + } + + // Dispose events to free memory + obs->dispose_events(); + + } // endfor: looped over observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Fill events into counts map * * @param[in] obs CTA observation. * * @exception GException::no_list * No event list found in observation. * - * This method maps the events found in a CTA events list into a sky map. + * Fills the events found in a CTA events list into a sky map. The method + * adds the events to the m_counts member. ***************************************************************************/ -void ctskymap::map_events(GCTAObservation* obs) +void ctskymap::fill_maps_counts(GCTAObservation* obs) { // Get non-const pointer on a CTA event list GCTAEventList* events = dynamic_cast @@ -745,11 +806,11 @@ void ctskymap::map_events(GCTAObservation* obs) // Determine sky pixel GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); GSkyDir dir = inst->dir(); - GSkyPixel pixel = m_skymap.dir2pix(dir); + GSkyPixel pixel = m_counts.dir2pix(dir); // Skip if pixel is out of range - if (pixel.x() < -0.5 || pixel.x() > (m_skymap.nx() - 0.5) || - pixel.y() < -0.5 || pixel.y() > (m_skymap.ny() - 0.5)) { + if (pixel.x() < -0.5 || pixel.x() > (m_counts.nx() - 0.5) || + pixel.y() < -0.5 || pixel.y() > (m_counts.ny() - 0.5)) { num_outside_map++; continue; } @@ -762,7 +823,7 @@ void ctskymap::map_events(GCTAObservation* obs) // Fill event in skymap #pragma omp critical(ctskymap_map_events) - m_skymap(pixel) += 1.0; + m_counts(pixel) += 1.0; num_in_map++; } // endfor: looped over all events @@ -779,7 +840,7 @@ void ctskymap::map_events(GCTAObservation* obs) // Write sky map into header log_header1(EXPLICIT, "Sky map"); - log_string(EXPLICIT, m_skymap.print(m_chatter)); + log_string(EXPLICIT, m_counts.print(m_chatter)); } // Return @@ -788,36 +849,7 @@ void ctskymap::map_events(GCTAObservation* obs) /***********************************************************************//** - * @brief Estimates the background in sky map for observation - * - * @param[in] obs CTA observation. - * - * Estimates the background in the sky map for a given observation and adds - * this estimate to the background sky map. The background estimation method - * is specified by the "bkgsubtract" parameter which can take the following - * values: - * - * NONE - No background estimation - * IRF - Background estimation based on IRF template - * RING - Ring background estimation - ***************************************************************************/ -void ctskymap::map_background(GCTAObservation* obs) -{ - // Dispatch to appropriate background estimation method - if (m_bkgsubtract == "IRF") { - map_background_irf(obs); - } - else if (m_bkgsubtract == "RING") { - map_background_irf(obs); - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Estimates the background in sky map based on IRF template + * @brief Compute background acceptance sky map based on IRF template * * @param[in] obs CTA observation. * @@ -825,11 +857,10 @@ void ctskymap::map_background(GCTAObservation* obs) * No response information available for observation. * No background template available in instrument response function. * - * Estimates the background in the sky map using the IRF template for a given - * observation and adds this estimate to the background sky map. The IRF - * template is integrated numerically in energy. + * Computes a background acceptance sky map using the IRF template for a + * given observation and adds the estimate to the acceptance sky map. ***************************************************************************/ -void ctskymap::map_background_irf(GCTAObservation* obs) +void ctskymap::fill_maps_acceptance(GCTAObservation* obs) { // Get IRF response const GCTAResponseIrf* rsp = dynamic_cast @@ -874,8 +905,8 @@ void ctskymap::map_background_irf(GCTAObservation* obs) // Initialise statistics double total = 0.0; - // Loop over all background map pixels - for (int i = 0; i < m_bkgmap.npix(); ++i) { + // Loop over all acceptance map pixels + for (int i = 0; i < m_acceptance.npix(); ++i) { // Get sky direction of pixel GSkyDir& skydir = m_dirs[i]; @@ -894,16 +925,16 @@ void ctskymap::map_background_irf(GCTAObservation* obs) // Multiply background rate with livetime and solid angle value *= exposure * m_solidangle[i]; - // Add number of background events to background map - #pragma omp critical(ctskymap_map_background_irf) - m_bkgmap(i) += value; + // Add number of background events to acceptance map + #pragma omp critical(ctskymap_map_acceptance) + m_acceptance(i) += value; // Update total number of background events total += value; - } // endfor: looped over background map pixels + } // endfor: looped over acceptance map pixels - // Log background subtraction results + // Log acceptance results #pragma omp critical(ctskymap_map_events) { log_value(NORMAL, "Events in background", int(total+0.5)); @@ -915,14 +946,24 @@ void ctskymap::map_background_irf(GCTAObservation* obs) /***********************************************************************//** - * @brief Compute the pixel significance if background subtraction is - * requested + * @brief Compute sky map, background map and significance map + * + * Computes the sky map, background map and significance map. The following + * background subtraction methods are supported: + * + * NONE + * ==== * - * Computes the pixel significance if a background subtraction method is - * specified. + * The sky map is simply the binned counts map, no background and + * significance maps are computed. * - * For the @c IRF method Poisson statistic (in the Gaussian limit) is - * assumed, and the significance \f$\sigma_i\f$ is computed using + * IRF + * === + * + * The background map is the acceptance map, and this background map is + * subtracted from the counts map and assigned to the sky map. The + * significance \f$\sigma_i\f$ is computed assuming Poisson statistic in + * the Gaussian limit, using * * \f[\sigma_i = \frac{N_i - B_i}{\sqrt{N_i}}\f] * @@ -930,6 +971,9 @@ void ctskymap::map_background_irf(GCTAObservation* obs) * \f$N_i\f$ is the number of observed counts and * \f$B_i\f$ is the estimated number of background counts. * + * RING + * ==== + * * For the @c RING method the Li & Ma significance is computed for each * pixel, using the @c roiradius, @c inradius and @c outradius parameters * to specify the radius of the On region and the Off background ring, @@ -941,21 +985,59 @@ void ctskymap::map_background_irf(GCTAObservation* obs) * map. Use FFT if the sky map is close to a cartesian grid, and use the * direct method if the sky map shows important distortions. ***************************************************************************/ -void ctskymap::map_significance(void) +void ctskymap::compute_maps(void) { - // Compute significance from "RING" method (Li & Ma eq. 17) - if (m_bkgsubtract == "RING") { + // Write header into logger + log_header1(TERSE, "Compute maps"); + + // If no background subtraction is requested then simply assign the + // counts map to the skymap + if (m_bkgsubtract == "NONE") { + m_skymap = m_counts; + } + + // else if IRF background is selected then use the acceptance as background + else if (m_bkgsubtract == "IRF") { + m_skymap = m_counts - m_acceptance; + m_bkgmap = m_acceptance; + m_sigmap = (m_counts - m_acceptance) / sqrt(m_counts); + } + + // else if RING background is selected then compute the On + else if (m_bkgsubtract == "RING") { + + // Compute initial RING background if (m_usefft) { - map_significance_ring_fft(); + compute_maps_ring_fft(); } else { - map_significance_ring_direct(); + compute_maps_ring_direct(); } - } - // Compute significance using Poisson statistics in the Gaussian limit - else { - m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap); + // Store copy of exclusion regions + GSkyMap exclmap = m_exclmap; + + // Iterative computation of exclusion regions + for (int iter = 0; iter < m_iterations; ++iter) { + + // Set exclusion region from significance map + m_exclmap = exclmap; + for (int i = 0; i < m_exclmap.npix(); ++i) { + if (m_sigmap(i) > m_threshold) { + m_exclmap(i) = 1.0; + } + } + + // Re-compute RING background + if (m_usefft) { + compute_maps_ring_fft(); + } + else { + compute_maps_ring_direct(); + } + + } // endfor: looped over iterations + } // Return @@ -964,77 +1046,79 @@ void ctskymap::map_significance(void) /***********************************************************************//** - * @brief Compute the pixel significance for RING background + * @brief Compute the maps for RING background using a FFT * * Computes the Li & Ma significance for each sky map pixel and replaces * the sky and background maps by the On- and Off-count maps. The computation - * is done by computing the exact distances between pixels in the sky map. + * is done using a FFT. ***************************************************************************/ -void ctskymap::map_significance_ring_direct(void) +void ctskymap::compute_maps_ring_fft(void) { // Log message about what is being done - log_header1(NORMAL, "Computing ring background map (direct method)"); - log_value(NORMAL, "Total pixels to process", m_skymap.npix()); + log_header3(NORMAL, "Computing ring background map (FFT method)"); - // Initialise On- and Off-count maps - GSkyMap onmap = m_skymap; - GSkyMap offmap = m_skymap; - onmap = 0.0; - offmap = 0.0; + // Initialise maps + m_skymap = m_counts; + m_bkgmap = m_counts; + m_sigmap = m_counts; + m_skymap = 0.0; + m_bkgmap = 0.0; + m_sigmap = 0.0; // 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; + // Set-up sky and acceptance maps with excluded pixels according to the + // exclusion map + GSkyMap excl_counts = m_counts; + GSkyMap excl_acceptance = m_acceptance; + for (int i = 0; i < m_counts.npix(); ++i) { + if (m_exclmap(i) != 0.0) { + excl_counts(i) = 0.0; + excl_acceptance(i) = 0.0; + } } - // Loop over sky map pixels - #pragma omp parallel for - for (int i = 0; i < m_skymap.npix(); ++i) { + // Convolve maps for On-region and Off-region + GSkyMap on_counts = ring_convolve(m_counts, 0.0, m_roiradius); + GSkyMap on_alpha = ring_convolve(m_acceptance, 0.0, m_roiradius); + GSkyMap off_counts = ring_convolve(excl_counts, m_inradius, m_outradius); + GSkyMap off_alpha = ring_convolve(excl_acceptance, m_inradius, m_outradius); - // 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 maps by loop over sky map pixels + for (int i = 0; i < m_sigmap.npix(); ++i) { - // Update user about progress (only if there are more than 10000 pixels) - if (i % n_logpix == 0 && n_logpix > 10000) { - #pragma omp critical(ctskymap_map_significance_ring) - log_value(NORMAL, "Pixels remaining", m_skymap.npix()-i); + // Get On-counts, Off-counts and alpha value for this bin + double n_on = on_counts(i); + double n_off = off_counts(i); + double alpha; + if (off_alpha(i) == 0.0) { + alpha = 0.0; + n_on = 0.0; + n_off = 0.0; + } + else { + alpha = on_alpha(i) / off_alpha(i); } - - // Compute the alpha and counts for this bin - compute_ring_values(i, m_skymap, m_bkgmap, n_on, n_off, alpha); // Store the On and alpha-weighted Off-counts - onmap(i) = n_on; - offmap(i) = alpha * n_off; + m_bkgmap(i) = alpha * n_off; + m_skymap(i) = n_on - m_bkgmap(i); // If alpha is zero then increment the bad alpha counter if (alpha == 0.0) { num_bad_alpha++; - #pragma omp critical(ctskymap_map_significance_ring) - m_sigmap(i) = 0.0; } - // ... otherwise store the results + // ... otherwise compute and store the results else { - - // Compute and store significance (Li & Ma eq. 17) - double sigma = sigma_li_ma(n_on, n_off, alpha); - #pragma omp critical(ctskymap_map_significance_ring) - m_sigmap(i) = sigma; - - } // endelse: alpha was non-zero + m_sigmap(i) = sigma_li_ma(n_on, n_off, alpha); + } } // endfor: looped over all pixels // Log the number of bad-alpha bins + log_value(NORMAL, "Total number of pixels", m_sigmap.npix()); log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha); // Now take the square root. Since some bins can be negative, first @@ -1042,11 +1126,6 @@ void ctskymap::map_significance_ring_direct(void) // each bin by its sign to preserve the +/- significance. m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap)); - // 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; } @@ -1057,68 +1136,75 @@ void ctskymap::map_significance_ring_direct(void) * * Computes the Li & Ma significance for each sky map pixel and replaces * the sky and background maps by the On- and Off-count maps. The computation - * is done using a FFT. + * is done by computing the exact distances between pixels in the sky map. ***************************************************************************/ -void ctskymap::map_significance_ring_fft(void) +void ctskymap::compute_maps_ring_direct(void) { // Log message about what is being done - log_header1(NORMAL, "Computing ring background map (FFT method)"); + log_header3(NORMAL, "Computing ring background map (direct method)"); + log_value(NORMAL, "Total pixels to process", m_counts.npix()); + + // Initialise maps + m_skymap = m_counts; + m_bkgmap = m_counts; + m_sigmap = m_counts; + m_skymap = 0.0; + m_bkgmap = 0.0; + m_sigmap = 0.0; // Initialise statistics int num_bad_alpha = 0; - // Set-up sky and acceptance maps with excluded pixels according to the - // exclusion map - GSkyMap skymap = m_skymap; - GSkyMap acceptance = m_bkgmap; - for (int i = 0; i < m_skymap.npix(); ++i) { - if (m_exclmap(i) != 0.0) { - skymap(i) = 0.0; - acceptance(i) = 0.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_counts.npix() / 10; + if (n_logpix > 100000) { + n_logpix = 100000; } - // Convolve maps for On-region and Off-region - GSkyMap on_counts = ring_convolve(skymap, 0.0, m_roiradius); - GSkyMap on_alpha = ring_convolve(acceptance, 0.0, m_roiradius); - GSkyMap off_counts = ring_convolve(skymap, m_inradius, m_outradius); - GSkyMap off_alpha = ring_convolve(acceptance, m_inradius, m_outradius); + // Compute maps by loop over sky map pixels + #pragma omp parallel for + for (int i = 0; i < m_sigmap.npix(); ++i) { - // Loop over sky map pixels - for (int i = 0; i < m_skymap.npix(); ++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; - // Get On-counts, Off-counts and alpha value for this bin - double n_on = on_counts(i); - double n_off = off_counts(i); - double alpha; - if (off_alpha(i) == 0.0) { - alpha = 0.0; - n_on = 0.0; - n_off = 0.0; - } - else { - alpha = on_alpha(i) / off_alpha(i); + // Update user about progress (only if there are more than 10000 pixels) + if (i % n_logpix == 0 && n_logpix > 10000) { + #pragma omp critical(ctskymap_map_significance_ring) + log_value(NORMAL, "Pixels remaining", m_sigmap.npix()-i); } + // Compute the alpha and counts for this bin + compute_ring_values(i, m_counts, m_acceptance, n_on, n_off, alpha); + // Store the On and alpha-weighted Off-counts - on_counts(i) = n_on; - off_counts(i) = alpha * n_off; + m_bkgmap(i) = alpha * n_off; + m_skymap(i) = n_on - m_bkgmap(i); // If alpha is zero then increment the bad alpha counter if (alpha == 0.0) { - m_sigmap(i) = 0.0; num_bad_alpha++; + #pragma omp critical(ctskymap_compute_maps_ring_direct) + m_sigmap(i) = 0.0; } - // ... otherwise compute and store the results + // ... otherwise store the results else { - m_sigmap(i) = sigma_li_ma(n_on, n_off, alpha); - } + + // Compute and store significance (Li & Ma eq. 17) + double sigma = sigma_li_ma(n_on, n_off, alpha); + #pragma omp critical(ctskymap_compute_maps_ring_direct) + m_sigmap(i) = sigma; + + } // endelse: alpha was non-zero } // endfor: looped over all pixels // Log the number of bad-alpha bins - log_value(NORMAL, "Total number of pixels", m_skymap.npix()); log_value(NORMAL, "Pixels with alpha=0", num_bad_alpha); // Now take the square root. Since some bins can be negative, first @@ -1126,11 +1212,6 @@ void ctskymap::map_significance_ring_fft(void) // each bin by its sign to preserve the +/- significance. m_sigmap = sign(m_sigmap) * sqrt(abs(m_sigmap)); - // Terminate by storing the on-counts map as sky map and the off-counts - // map as background map - m_skymap = on_counts; - m_bkgmap = off_counts; - // Return return; } @@ -1281,21 +1362,33 @@ void ctskymap::ring_bounding_box(const int& ipixel, int& ix1, int& ix2, int& iy1, int& iy2) const { // Get number of pixels in x and y direction - int nx = m_skymap.nx(); - int ny = m_skymap.ny(); + int nx = m_counts.nx(); + int ny = m_counts.ny(); // Get x and y index of pixel to consider int ix0 = ipixel % nx; int iy0 = ipixel / nx; + // Get pointer on WCS projection + const GWcs* wcs = dynamic_cast(m_counts.projection()); + if (wcs == NULL) { + std::string msg = "Sky map is not a WCS projection. Method is only " + "valid for WCS projections."; + throw GException::invalid_value(G_RING_BOUNDING_BOX, msg); + } + + // Get X and Y step size + double dx_binsz = wcs->cdelt(0); + double dy_binsz = wcs->cdelt(1); + // Compute pixel increment in x-direction double dx = 0.0; int ix = (ix0 > 0) ? ix0 - 1 : ix0 + 1; if (ix < nx) { dx = m_dirs[ipixel].dist_deg(m_dirs[ix+iy0*nx]); } - if (dx < m_binsz) { - dx = m_binsz; + if (dx < dx_binsz) { + dx = dx_binsz; } // Compute pixel increment in y-direction @@ -1304,8 +1397,8 @@ void ctskymap::ring_bounding_box(const int& ipixel, int& ix1, int& ix2, if (iy < ny) { dx = m_dirs[ipixel].dist_deg(m_dirs[ix0+iy*nx]); } - if (dy < m_binsz) { - dy = m_binsz; + if (dy < dy_binsz) { + dy = dy_binsz; } // Compute bounding box half size in x and y. The outer radius is @@ -1410,7 +1503,7 @@ GSkyMap ctskymap::ring_convolve(const GSkyMap& map, const double& rmin, GNdarray ctskymap::ring_kernel(const double& rmin, const double& rmax) const { // Get pointer on WCS projection - const GWcs* wcs = dynamic_cast(m_skymap.projection()); + const GWcs* wcs = dynamic_cast(m_counts.projection()); if (wcs == NULL) { std::string msg = "Sky map is not a WCS projection. Method is only " "valid for WCS projections."; @@ -1422,24 +1515,24 @@ GNdarray ctskymap::ring_kernel(const double& rmin, const double& rmax) const double dy = wcs->cdelt(1); // Initialise kernel - GNdarray kern(m_skymap.nx(), m_skymap.ny()); + GNdarray kern(m_counts.nx(), m_counts.ny()); // Fill kernel - for (int ix1 = 0, ix2 = m_skymap.nx(); ix1 < m_skymap.nx(); ++ix1, --ix2) { + for (int ix1 = 0, ix2 = m_counts.nx(); ix1 < m_counts.nx(); ++ix1, --ix2) { double x = ix1 * dx; double xqs = x * x; - for (int iy1 = 0, iy2 = m_skymap.ny(); iy1 < m_skymap.ny(); ++iy1, --iy2) { + for (int iy1 = 0, iy2 = m_counts.ny(); iy1 < m_counts.ny(); ++iy1, --iy2) { double y = iy1 * dy; double r = std::sqrt(xqs + y*y); if ((r >= rmin) && (r <= rmax)) { kern(ix1,iy1) += 1.0; - if (ix2 < m_skymap.nx()) { + if (ix2 < m_counts.nx()) { kern(ix2,iy1) += 1.0; } - if (iy2 < m_skymap.ny()) { + if (iy2 < m_counts.ny()) { kern(ix1,iy2) += 1.0; } - if ((ix2 < m_skymap.nx()) && (iy2 < m_skymap.ny())) { + if ((ix2 < m_counts.nx()) && (iy2 < m_counts.ny())) { kern(ix2,iy2) += 1.0; } } @@ -1486,6 +1579,34 @@ double ctskymap::sigma_li_ma(const double& n_on, /***********************************************************************//** + * @brief Write sky map into FITS file + * + * @param[in,out] fits FITS file. + * @param[in] map Sky map. + * @param[in] extname Extension name. + * + * Write one sky map with the extension name and keywords into the FITS file. + ***************************************************************************/ +void ctskymap::write_map(GFits& fits, const GSkyMap& map, const std::string& extname) const +{ + // Write map into FITS file + GFitsHDU* hdu = map.write(fits); + + // Set map extension name + if (hdu != NULL) { + hdu->extname(extname); + } + + // Write keywords into map extension + write_ogip_keywords(hdu); + write_hdu_keywords(hdu); + + // Return + return; +} + + +/***********************************************************************//** * @brief Write keywords in FITS HDU * * @param[in,out] hdu Pointer to FITS HDU. @@ -1505,9 +1626,13 @@ void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const // 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"); + hdu->card("USEFFT", m_usefft, "Use FFT for RING background"); + 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"); + hdu->card("ITER", m_iterations, "Exclusion map iterations"); + hdu->card("THRES", m_threshold, "Exclusion map threshold"); + hdu->card("EXCLMAP", m_inexclusion.url(), "Exclusion map name"); } } // endif: pointer was valid diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index a72a1e7..f23e292 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -69,15 +69,15 @@ protected: 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_significance(void); - void map_significance_ring_direct(void); - void map_significance_ring_fft(void); + void setup_exclusion_map(const GFilename& filename); + void setup_exclusion_map_fits(const GFilename& filename); + void setup_exclusion_map_region(const GFilename& filename); + void fill_maps(void); + void fill_maps_counts(GCTAObservation* obs); + void fill_maps_acceptance(GCTAObservation* obs); + void compute_maps(void); + void compute_maps_ring_fft(void); + void compute_maps_ring_direct(void); void compute_ring_values(const int& ipixel, const GSkyMap& counts, const GSkyMap& background, @@ -92,26 +92,33 @@ protected: double sigma_li_ma(const double& n_on, const double& n_off, const double& alpha) const; + void write_map(GFits& fits, + const GSkyMap& map, + const std::string& extname) const; void write_hdu_keywords(GFitsHDU* hdu) const; // User parameters GFilename m_outmap; //!< Output file name GFilename m_inexclusion; //!< Exclusion map file name - double m_binsz; //!< Sky map binsize double m_emin; //!< Minimum energy (TeV) double m_emax; //!< Maximum energy (TeV) std::string m_bkgsubtract; //!< Background subtraction method double m_roiradius; //!< Region of interest radius for RING bkg. double m_inradius; //!< Inner ring radius for RING background double m_outradius; //!< Outer ring radius for RING background + int m_iterations; //!< Number of iterations for RING background + double m_threshold; //!< Threshold for RING background bool m_usefft; //!< Use FFT for RING background bool m_publish; //!< Publish sky map? GChatter m_chatter; //!< Chattiness // Protected members + bool m_has_inmap; //!< Has valid input map GSkyMap m_skymap; //!< Sky map GSkyMap m_bkgmap; //!< Background map GSkyMap m_sigmap; //!< Significance map + GSkyMap m_counts; //!< Counts map + GSkyMap m_acceptance; //!< Acceptance map GSkyMap m_exclmap; //!< Exclusion map for RING background // Caching variables to prevent multiple computations of the same thing diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par index ef3357f..022967e 100644 --- a/src/ctskymap/ctskymap.par +++ b/src/ctskymap/ctskymap.par @@ -19,6 +19,7 @@ inobs, f, a, events.fits,,, "Input event list or observation definition XML file" caldb, s, a, prod2,,, "Calibration database" irf, s, a, South_0.5h,,, "Instrument response function" +inmap, f, h, NONE,,, "Input skymap file" outmap, f, a, skymap.fits,,, "Output skymap file" # @@ -46,8 +47,10 @@ 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 background ring radius (degrees)" outradius, r, a, 0.8,,, "Outer background ring radius (degrees)" +iterations, i, a, 0,0,100, "Number of iterations for exclusion regions computation" +threshold, r, a, 5.0,,, "Significance threshold for exclusion regions computation" inexclusion, f, h, NONE,,, "Exclusion map or DS9 region file name" -usefft, b, h, yes,,, "Use FFT for ring computation?" +usefft, b, h, yes,,, "Use FFT for ring computation?" # # Standard parameters diff --git a/src/support/ctobservation.cpp b/src/support/ctobservation.cpp index d15ee6f..e6511ec 100644 --- a/src/support/ctobservation.cpp +++ b/src/support/ctobservation.cpp @@ -319,6 +319,48 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const /***********************************************************************//** + * @brief Read OGIP keywords from FITS HDU + * + * @param[in,out] hdu Pointer to FITS HDU. + * + * Read OGIP keywords from FITS HDU. + ***************************************************************************/ +void ctobservation::read_ogip_keywords(GFitsHDU* hdu) const +{ + // Continue only if pointer is valid + if (hdu != NULL) { + + // Get observation information + m_ogip_telescope = (hdu->has_card("TELESCOP")) ? hdu->string("TELESCOP") : ""; + + // Get observation time information + std::string date_obs = (hdu->has_card("DATE-OBS")) ? hdu->string("DATE-OBS") : ""; + std::string time_obs = (hdu->has_card("TIME-OBS")) ? hdu->string("TIME-OBS") : ""; + std::string date_end = (hdu->has_card("DATE-END")) ? hdu->string("DATE-END") : ""; + std::string time_end = (hdu->has_card("TIME-END")) ? hdu->string("TIME-END") : ""; + m_ogip_telapse = (hdu->has_card("TELAPSE")) ? hdu->real("TELAPSE") : 0.0; + m_ogip_ontime = (hdu->has_card("ONTIME")) ? hdu->real("ONTIME") : 0.0; + m_ogip_livetime = (hdu->has_card("LIVETIME")) ? hdu->real("LIVETIME") : 0.0; + m_ogip_exposure = (hdu->has_card("EXPOSURE")) ? hdu->real("EXPOSURE") : 0.0; + + // Get OGIP start and stop time + std::string tstart = date_obs + "T" + time_obs; + std::string tstop = date_end + "T" + time_end; + if (tstart != "T") { + m_ogip_tstart.utc(tstart); + } + if (tstop != "T") { + m_ogip_tstop.utc(tstop); + } + + } // endif: pointer was valid + + // Return + return; +} + + +/***********************************************************************//** * @brief Write OGIP keywords in FITS HDU * * @param[in,out] hdu Pointer to FITS HDU. diff --git a/src/support/ctobservation.hpp b/src/support/ctobservation.hpp index ea823e6..797ac92 100644 --- a/src/support/ctobservation.hpp +++ b/src/support/ctobservation.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctobservation - Base class for observation tools * * ----------------------------------------------------------------------- * - * copyright (C) 2016-2017 by Juergen Knoedlseder * + * copyright (C) 2016-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -73,6 +73,7 @@ protected: GCTAObservation* next_unbinned_observation(void); const GCTAObservation* first_unbinned_observation(void) const; const GCTAObservation* next_unbinned_observation(void) const; + void read_ogip_keywords(GFitsHDU* hdu) const; void write_ogip_keywords(GFitsHDU* hdu) const; void set_obs_statistic(const std::string& statistic); void save_events_fits(void); diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index 119bead..faf9125 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -650,7 +650,11 @@ GEbounds ctool::create_ebounds(void) ***************************************************************************/ GSkyMap ctool::create_map(const GObservations& obs) { - // Read task parameters + // Read coordinate system and projection + std::string coordsys = (*this)["coordsys"].string(); + std::string proj = (*this)["proj"].string(); + + // Read sky map reference double xref = 0.0; double yref = 0.0; bool usepnt = (*this)["usepnt"].boolean(); @@ -658,11 +662,11 @@ GSkyMap ctool::create_map(const GObservations& obs) xref = (*this)["xref"].real(); yref = (*this)["yref"].real(); } - std::string proj = (*this)["proj"].string(); - std::string coordsys = (*this)["coordsys"].string(); - double binsz = (*this)["binsz"].real(); - int nxpix = (*this)["nxpix"].integer(); - int nypix = (*this)["nypix"].integer(); + + // Read sky map pixel size and number of bins + double binsz = (*this)["binsz"].real(); + int nxpix = (*this)["nxpix"].integer(); + int nypix = (*this)["nypix"].integer(); // If requested, get pointing from observations if (usepnt) { -- libgit2 0.21.4