diff --git a/doc/source/users/reference_manual/ctskymap.rst b/doc/source/users/reference_manual/ctskymap.rst index e9a4ebc..f4c3591 100644 --- a/doc/source/users/reference_manual/ctskymap.rst +++ b/doc/source/users/reference_manual/ctskymap.rst @@ -21,56 +21,68 @@ If ``IRF`` is selected, the background template that are shipped with the Instrument Response Functions will be used for background subtraction. ctskymap generates a FITS file comprising a sky map as primary extension. +If a background subtraction method was selected, the FITS file will contain +the additional extensions ``BACKGRAOUND`` and ``SIGNIFICANCE`` that contain +the background map and a significance map, respectively. For the latter, the +significance of the signal is computed for each sky map pixel. General parameters ------------------ ``inobs [file]`` - Input event list or observation definition XML file - + Input event list or observation definition XML file. + +``caldb [string]`` + Calibration database (only required for IRF background subtraction if no + response information is provided by ``inobs``). + +``irf [string]`` + Instrument response function (only required for IRF background subtraction + if no response information is provided by ``inobs``). + ``outmap [file]`` - Output sky map file + Output sky map file. ``emin [real]`` - Minimum energy in map (in TeV) + Minimum energy in map (in TeV). ``emax [real]`` - Maximum energy in map (in TeV) + Maximum energy in map (in TeV). ``(usepnt = no) [boolean]`` Use pointing instead of xref/yref parameters? ``nxpix [integer]`` - Size of the Right Ascension / Galactic longitude axis (in pixels) + Size of the Right Ascension / Galactic longitude axis (in pixels). ``nypix [integer]`` - Size of the Declination / Galactic latitude axis (in pixels) + Size of the Declination / Galactic latitude axis (in pixels). ``binsz [real]`` Pixel size (in degrees/pixel). ``coordsys [string]`` - Coordinate system (CEL - celestial, GAL - galactic) + Coordinate system (CEL - celestial, GAL - galactic). ``proj [string]`` - Projection method + Projection method. ``xref [real]`` - Right Ascension / Galactic longitude of image centre (J2000, in degrees) + Right Ascension / Galactic longitude of image centre (J2000, in degrees). ``yref [real]`` - Declination / Galactic latitude of image centre (J2000, in degrees) + Declination / Galactic latitude of image centre (J2000, in degrees). ``bkgsubtract [string]`` - Background subtraction method + Background subtraction method. Standard parameters ------------------- ``(publish = no) [boolean]`` - Specifies whether the sky map should be published on VO Hub + Specifies whether the sky map should be published on VO Hub. ``(chatter = 2) [integer]`` Verbosity of the executable: @@ -85,16 +97,17 @@ Standard parameters chatter = 4: detailed report about the task execution ``(clobber = yes) [boolean]`` - Specifies whether an existing output sky map file should be overwritten + Specifies whether an existing output sky map file should be overwritten. ``(debug = no) [boolean]`` - Enables debug mode. In debug mode the executable will dump any log file output to the console. + Enables debug mode. In debug mode the executable will dump any log file + output to the console. ``(mode = ql) [string]`` - Mode of automatic parameters (default is "ql", i.e. "query and learn") + Mode of automatic parameters (default is "ql", i.e. "query and learn"). ``(logfile = ctskymap.log) [string]`` - Name of log file + Name of log file. Related tools or scripts diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 8605fcd..36cb64d 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -233,11 +233,23 @@ void ctskymap::run(void) // Map events into sky map map_events(obs); - // Subtract background from sky map - bkg_subtract(obs); + // Compute background sky map + map_background(obs); } // endfor: looped over observations + // If background subtract was required then compute significance map + // and subtract background map from counts map + if (m_bkgsubtract != "NONE") { + + // Compute significance map + m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap); + + // Subtract background map from counts map + m_skymap -= m_bkgmap; + + } // endif: background subtraction was requested + // Optionally publish sky map if (m_publish) { publish(); @@ -271,6 +283,24 @@ void ctskymap::save(void) // Write sky map into FITS file m_skymap.write(fits); + // 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 + m_bkgmap.write(fits); + + // Set background map extension name + fits[fits.size()-1]->extname("BACKGROUND"); + + // Write significance map into FITS file + m_sigmap.write(fits); + + // Set significance map extension name + fits[fits.size()-1]->extname("SIGNIFICANCE"); + + } // endif: background subtraction was requested + // Save FITS file to disk fits.saveto(m_outmap, clobber()); @@ -328,6 +358,8 @@ void ctskymap::init_members(void) { // Initialise members m_skymap.clear(); + m_bkgmap.clear(); + m_sigmap.clear(); m_emin = 0.0; m_emax = 0.0; m_bkgsubtract = "NONE"; @@ -348,6 +380,8 @@ 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_emin = app.m_emin; m_emax = app.m_emax; m_bkgsubtract = app.m_bkgsubtract; @@ -387,10 +421,18 @@ void ctskymap::get_parameters(void) // Create sky map based on task parameters m_skymap = create_map(m_obs); - // Get remaining parameters + // Get further parameters m_emin = (*this)["emin"].real(); m_emax = (*this)["emax"].real(); m_bkgsubtract = (*this)["bkgsubtract"].string(); + + // If IRF background subtraction is requested then make sure that the + // CTA observations in the observation container have response information + if (m_bkgsubtract == "IRF") { + set_response(m_obs); + } + + // Get remaining parameters m_publish = (*this)["publish"].boolean(); m_chatter = static_cast((*this)["chatter"].integer()); @@ -399,6 +441,13 @@ void ctskymap::get_parameters(void) m_outmap = (*this)["outmap"].filename(); } + // Create background map and significance map if background subtraction + // is requested + if (m_bkgsubtract != "NONE") { + m_bkgmap = create_map(m_obs); + m_sigmap = create_map(m_obs); + } + // Write parameters into logger log_parameters(TERSE); @@ -487,22 +536,23 @@ void ctskymap::map_events(GCTAObservation* obs) /***********************************************************************//** - * @brief Subtract background events from skymap + * @brief Estimates the background in sky map for observation * * @param[in] obs CTA observation. * - * Subtracts background events from the sky map using various methods. The - * background subtraction method is specified by the "bkgsubtract" parameter - * which can take the following values: + * 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 subtraction - * IRF - Background subtraction based on IRF template + * NONE - No background estimation + * IRF - Background estimation based on IRF template ***************************************************************************/ -void ctskymap::bkg_subtract(GCTAObservation* obs) +void ctskymap::map_background(GCTAObservation* obs) { - // Dispatch to apporpriate background subtraction method + // Dispatch to appropriate background estimation method if (m_bkgsubtract == "IRF") { - bkg_subtract_irf(obs); + map_background_irf(obs); } // Return @@ -511,7 +561,7 @@ void ctskymap::bkg_subtract(GCTAObservation* obs) /***********************************************************************//** - * @brief Subtract background events based on IRF template from sky map + * @brief Estimates the background in sky map based on IRF template * * @param[in] obs CTA observation. * @@ -519,9 +569,10 @@ void ctskymap::bkg_subtract(GCTAObservation* obs) * No response information available for observation. * No background template available in instrument response function. * - * Subtracts background events from the sky map using the IRF template. + * Estimates the background in the sky map using the IRF template for a given + * observation and adds this estimate to the background sky map. ***************************************************************************/ -void ctskymap::bkg_subtract_irf(GCTAObservation* obs) +void ctskymap::map_background_irf(GCTAObservation* obs) { // Get IRF response const GCTAResponseIrf* rsp = dynamic_cast(obs->response()); @@ -563,11 +614,11 @@ void ctskymap::bkg_subtract_irf(GCTAObservation* obs) // Setup energy nodes for computation GEnergies energies(nebin, emin, emax); // More complicated nbins later - // Loop over all sky map pixels - for (int i = 0; i < m_skymap.npix(); ++i) { + // Loop over all background map pixels + for (int i = 0; i < m_bkgmap.npix(); ++i) { // Get sky direction of pixel - GSkyDir skydir = m_skymap.inx2dir(i); + GSkyDir skydir = m_bkgmap.inx2dir(i); // Convert sky direction in instrument direction GCTAInstDir instdir = obs->pointing().instdir(skydir); @@ -587,12 +638,12 @@ void ctskymap::bkg_subtract_irf(GCTAObservation* obs) } // Multiply background rate with livetime and solid angle - value *= obs->livetime() * m_skymap.solidangle(i); + value *= obs->livetime() * m_bkgmap.solidangle(i); - // Subtract background value from sky map - m_skymap(i,0) -= value; + // Add number of background events to background map + m_bkgmap(i,0) += value; - } // endfor: looped over sky map pixels + } // endfor: looped over background map pixels // Return return; diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index e62b3ae..01aa6a1 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -70,8 +70,8 @@ protected: void free_members(void); void get_parameters(void); void map_events(GCTAObservation* obs); - void bkg_subtract(GCTAObservation* obs); - void bkg_subtract_irf(GCTAObservation* obs); + void map_background(GCTAObservation* obs); + void map_background_irf(GCTAObservation* obs); // User parameters GFilename m_outmap; //!< Output file name @@ -83,6 +83,8 @@ protected: // Protected members GSkyMap m_skymap; //!< Sky map + GSkyMap m_bkgmap; //!< Background map + GSkyMap m_sigmap; //!< Significance map }; diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par index a8dfbaa..6608657 100644 --- a/src/ctskymap/ctskymap.par +++ b/src/ctskymap/ctskymap.par @@ -17,6 +17,8 @@ # General parameters #=================== 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" outmap, f, a, "skymap.fits",,, "Output skymap file" # diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py index 862bd34..38db35a 100755 --- a/test/test_ctskymap.py +++ b/test/test_ctskymap.py @@ -70,7 +70,7 @@ class Test(test): ' outmap="ctskymap_cmd1.fits"'+ \ ' emin=0.1 emax=100.0 nxpix=200 nypix=200'+ \ ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \ - ' xref=83.63 yref=22.01'+ \ + ' xref=83.63 yref=22.01 bkgsubtract="NONE"'+ \ ' logfile="ctskymap_cmd1.log" chatter=1' # Check if execution of wrong command fails @@ -89,7 +89,7 @@ class Test(test): ' outmap="ctskymap_cmd2.fits"'+ \ ' emin=0.1 emax=100.0 nxpix=200 nypix=200'+ \ ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \ - ' xref=83.63 yref=22.01'+ \ + ' xref=83.63 yref=22.01 bkgsubtract="NONE"'+ \ ' logfile="ctskymap_cmd2.log" chatter=2' # Check if execution failed @@ -134,19 +134,20 @@ class Test(test): skymap.clear() # Now set ctskymap parameters - skymap['inobs'] = self._events - skymap['emin'] = 0.1 - skymap['emax'] = 100 - skymap['nxpix'] = 200 - skymap['nypix'] = 200 - skymap['binsz'] = 0.02 - skymap['coordsys'] = 'CEL' - skymap['proj'] = 'CAR' - skymap['xref'] = 83.63 - skymap['yref'] = 22.01 - skymap['outmap'] = 'ctskymap_py1.fits' - skymap['logfile'] = 'ctskymap_py1.log' - skymap['chatter'] = 2 + skymap['inobs'] = self._events + skymap['emin'] = 0.1 + skymap['emax'] = 100 + skymap['nxpix'] = 200 + skymap['nypix'] = 200 + skymap['binsz'] = 0.02 + skymap['coordsys'] = 'CEL' + skymap['proj'] = 'CAR' + skymap['xref'] = 83.63 + skymap['yref'] = 22.01 + skymap['bkgsubtract'] = 'NONE' + skymap['outmap'] = 'ctskymap_py1.fits' + skymap['logfile'] = 'ctskymap_py1.log' + skymap['chatter'] = 2 # Run ctskymap tool skymap.logFileOpen() # Make sure we get a log file @@ -190,18 +191,19 @@ class Test(test): # Allocate ctskymap tool from observation container skymap = ctools.ctskymap(obs) - skymap['emin'] = 0.1 - skymap['emax'] = 100 - skymap['nxpix'] = 200 - skymap['nypix'] = 200 - skymap['binsz'] = 0.02 - skymap['coordsys'] = 'GAL' - skymap['proj'] = 'CAR' - skymap['xref'] = 184.5575 - skymap['yref'] = -5.7844 - skymap['outmap'] = 'ctskymap_py3.fits' - skymap['logfile'] = 'ctskymap_py3.log' - skymap['chatter'] = 4 + skymap['emin'] = 0.1 + skymap['emax'] = 100 + skymap['nxpix'] = 200 + skymap['nypix'] = 200 + skymap['binsz'] = 0.02 + skymap['coordsys'] = 'GAL' + skymap['proj'] = 'CAR' + skymap['xref'] = 184.5575 + skymap['yref'] = -5.7844 + skymap['bkgsubtract'] = 'NONE' + skymap['outmap'] = 'ctskymap_py3.fits' + skymap['logfile'] = 'ctskymap_py3.log' + skymap['chatter'] = 4 # Execute tool skymap.logFileOpen() # Needed to get a new log file @@ -213,6 +215,31 @@ class Test(test): # Publish with name skymap.publish('My sky map') + # Allocate ctskymap tool from observation container + skymap = ctools.ctskymap(obs) + skymap['emin'] = 0.1 + skymap['emax'] = 100 + skymap['nxpix'] = 200 + skymap['nypix'] = 200 + skymap['binsz'] = 0.02 + skymap['coordsys'] = 'GAL' + skymap['proj'] = 'CAR' + skymap['xref'] = 184.5575 + skymap['yref'] = -5.7844 + skymap['bkgsubtract'] = 'IRF' + skymap['caldb'] = self._caldb + skymap['irf'] = self._irf + skymap['outmap'] = 'ctskymap_py4.fits' + skymap['logfile'] = 'ctskymap_py4.log' + skymap['chatter'] = 4 + + # Execute tool + skymap.logFileOpen() # Needed to get a new log file + skymap.execute() + + # Check result file + self._check_result_file('ctskymap_py3.fits') + # Return return