diff --git a/ChangeLog b/ChangeLog index b4ef0c8..c55a3e2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-11-23 +2016-11-25 * Version 1.2.0 released ======================== + Add background subtraction to ctskymap tool (#1882) Add csobsselect script (#1881) Add csmodelselect script (#1879) Add startindex parameter to ctobssim (#1878) diff --git a/NEWS b/NEWS index 2fcb748..e556c2c 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.2.0 -23 November 2016 +25 November 2016 Introduction @@ -123,6 +123,11 @@ of the tool (#1835). The ctskymap::map() method has been renamed to ctskymap::skymap() to avoid a name conflict with Python (#1828). +The tool now also supports subtraction of background models that can be enabled +using the "bkgsubtract" parameter. If this parameter is set to "IRF" the +background templates shipped with the Instrument Response Function will be +subtracted from the sky map (#1882). + ctmodel - CTA counts cube modelling ----------------------------------- diff --git a/doc/source/users/reference_manual/ctskymap.rst b/doc/source/users/reference_manual/ctskymap.rst index e7abc02..e9a4ebc 100644 --- a/doc/source/users/reference_manual/ctskymap.rst +++ b/doc/source/users/reference_manual/ctskymap.rst @@ -14,6 +14,12 @@ provided in an observation definition file. The tool will loop over all event lists that are provided and fill all events into a single sky map. Only events within an energy interval spanned by ``emin`` and ``emax`` are considered. +Optionally, the tool will subtract a background model from the sky map. The +background subtraction method can be selected using the ``bkgsubtract`` +parameter. By default, no background model is subtracted (method ``NONE``). +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. @@ -21,47 +27,50 @@ General parameters ------------------ ``inobs [file]`` - Input event list or observation definition XML file. + Input event list or observation definition XML file ``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 + ``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) -``proj [string]`` - Projection method. +``bkgsubtract [string]`` + 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: @@ -76,16 +85,16 @@ 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. ``(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/doc/source/users/tutorials/1dc/first_skymap.png b/doc/source/users/tutorials/1dc/first_skymap.png index 2520697..e421cbe 100644 --- a/doc/source/users/tutorials/1dc/first_skymap.png +++ b/doc/source/users/tutorials/1dc/first_skymap.png diff --git a/doc/source/users/tutorials/1dc/first_skymap.rst b/doc/source/users/tutorials/1dc/first_skymap.rst index 5976cad..06c074f 100644 --- a/doc/source/users/tutorials/1dc/first_skymap.rst +++ b/doc/source/users/tutorials/1dc/first_skymap.rst @@ -21,14 +21,15 @@ You do this with the :ref:`ctskymap` tool by typing: Size of the Y axis in pixels [200] 250 Lower energy limit (TeV) [0.1] Upper energy limit (TeV) [100.0] + Background subtraction method (NONE|IRF) [NONE] Output skymap file [skymap.fits] -We generate here a sky map centred on the Galactic Centre in Galactic +We generated here a sky map centred on the Galactic Centre in Galactic coordinates using a cartesian projection. The sky map is 10° wide and 5° high, with an image scale of 0.02° per pixel. -All events between 0.1 and 100 TeV have been collected in the skymap. -The :ref:`ctskymap` tool writes the result into the FITS file ``skymap.fits`` -that will be created in the ``my_first_analysis`` folder. +All events between 0.1 and 100 TeV were collected in the skymap. +The :ref:`ctskymap` tool wrote the sky map into the FITS file ``skymap.fits`` +that was created in the ``my_first_analysis`` folder. The FITS file can be displayed using for example ds9 (see below): .. figure:: first_skymap.png @@ -37,4 +38,40 @@ The FITS file can be displayed using for example ds9 (see below): *Sky map of the events recorded during the Galactic Plane Survey around the Galactic Centre* -Bla +The sky map shows a wide-spread distribution of events with a number of sources +superimposed. +Most of the events are due to irreducable background that hampers the +recognition of the gamma-ray sources. +To describe the irreducable background in the CTA data, templates of the +background event distribution are shipped together with the +:ref:`Instrument Response Functions `. +These template can be used by :ref:`ctskymap` to subtract the background +contribution. +To enable the background subtraction you should re-run :ref:`ctskymap` with +the background subtraction method set to ``IRF`` as follows: + +.. code-block:: bash + + $ ctskymap + Input event list or observation definition XML file [events.fits] obs.xml + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] 0.0 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] 0.0 + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] GAL + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] 500 + Size of the Y axis in pixels [200] 250 + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Background subtraction method (NONE|IRF) [NONE] IRF + Output skymap file [skymap.fits] + +The figure below shows the resulting sky map. +The irreducable background has been subtracted and the sources of gamma-ray +emission are now clearly descernable. + +.. figure:: first_skymap_bkgsubtract.png + :width: 600px + :align: center + + *Background subtracted sky map of the events recorded during the Galactic Plane Survey around the Galactic Centre* diff --git a/doc/source/users/tutorials/1dc/first_skymap_bkgsubtract.png b/doc/source/users/tutorials/1dc/first_skymap_bkgsubtract.png new file mode 100644 index 0000000..9473e70 --- /dev/null +++ b/doc/source/users/tutorials/1dc/first_skymap_bkgsubtract.png diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index a756c18..8605fcd 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -35,6 +35,7 @@ /* __ Method name definitions ____________________________________________ */ #define G_INIT_MAP "ctskymap::init_map(GCTAObservation* obs)" #define G_BIN_EVENTS "ctskymap::bin_events(GCTAObservation*)" +#define G_BKG_SUBTRACT_IRF "ctskymap::bkg_subtract_irf(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ @@ -231,7 +232,10 @@ void ctskymap::run(void) // Map events into sky map map_events(obs); - + + // Subtract background from sky map + bkg_subtract(obs); + } // endfor: looped over observations // Optionally publish sky map @@ -324,10 +328,11 @@ void ctskymap::init_members(void) { // Initialise members m_skymap.clear(); - m_emin = 0.0; - m_emax = 0.0; - m_publish = false; - m_chatter = static_cast(2); + m_emin = 0.0; + m_emax = 0.0; + m_bkgsubtract = "NONE"; + m_publish = false; + m_chatter = static_cast(2); // Return return; @@ -342,11 +347,12 @@ void ctskymap::init_members(void) void ctskymap::copy_members(const ctskymap& app) { // Copy attributes - m_skymap = app.m_skymap; - m_emin = app.m_emin; - m_emax = app.m_emax; - m_publish = app.m_publish; - m_chatter = app.m_chatter; + m_skymap = app.m_skymap; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_bkgsubtract = app.m_bkgsubtract; + m_publish = app.m_publish; + m_chatter = app.m_chatter; // Return return; @@ -382,10 +388,11 @@ void ctskymap::get_parameters(void) m_skymap = create_map(m_obs); // Get remaining parameters - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - m_publish = (*this)["publish"].boolean(); - m_chatter = static_cast((*this)["chatter"].integer()); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_bkgsubtract = (*this)["bkgsubtract"].string(); + m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read ahead parameters if (read_ahead()) { @@ -407,16 +414,18 @@ void ctskymap::get_parameters(void) * * @exception GException::no_list * No event list found in observation. - * @exception GCTAException::no_pointing - * No valid CTA pointing found. * * This method maps the events found in a CTA events list into a sky map. ***************************************************************************/ void ctskymap::map_events(GCTAObservation* obs) { + // Get non-const pointer on a CTA event list + GCTAEventList* events = dynamic_cast + (const_cast(obs->events())); + // Make sure that the observation holds a CTA event list. If this // is not the case then throw an exception. - if (dynamic_cast(obs->events()) == NULL) { + if (events == NULL) { throw GException::no_list(G_BIN_EVENTS); } @@ -433,8 +442,6 @@ void ctskymap::map_events(GCTAObservation* obs) int num_in_map = 0; // Fill sky map - GCTAEventList* events = static_cast - (const_cast(obs->events())); for (int i = 0; i < events->size(); ++i) { // Get event @@ -477,3 +484,116 @@ void ctskymap::map_events(GCTAObservation* obs) // Return return; } + + +/***********************************************************************//** + * @brief Subtract background events from skymap + * + * @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: + * + * NONE - No background subtraction + * IRF - Background subtraction based on IRF template + ***************************************************************************/ +void ctskymap::bkg_subtract(GCTAObservation* obs) +{ + // Dispatch to apporpriate background subtraction method + if (m_bkgsubtract == "IRF") { + bkg_subtract_irf(obs); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Subtract background events based on IRF template from sky map + * + * @param[in] obs CTA observation. + * + * @exception GException::invalid_value + * 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. + ***************************************************************************/ +void ctskymap::bkg_subtract_irf(GCTAObservation* obs) +{ + // Get IRF response + 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_BKG_SUBTRACT_IRF, msg); + } + + // 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_BKG_SUBTRACT_IRF, msg); + } + + // Setup energy range covered by data + GEnergy emin; + GEnergy emax; + GEbounds ebds; + emin.TeV(m_emin); + emax.TeV(m_emax); + + // Determine number of energy nodes, 20 per decade, at least 10 bins + int nebin = int((emax.log10TeV() - emin.log10TeV()) * 20.0); + if (nebin < 10) { + nebin = 10; + } + + // 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) { + + // Get sky direction of pixel + GSkyDir skydir = m_skymap.inx2dir(i); + + // Convert sky direction in instrument direction + GCTAInstDir instdir = obs->pointing().instdir(skydir); + + // Compute background value for each energy node + std::vector values(energies.size(), 0.0); + for (int k = 0; k < energies.size(); ++k) { + values[k] = (*bkg)(energies[k].log10TeV(), instdir.detx(), + instdir.dety()); + } + + // Integrate background values according to Trapezoid rule + double value = 0; + for (int k = 1; k < energies.size(); ++k) { + value += 0.5 * (values[k] + values[k-1]) * + (energies[k].MeV() - energies[k-1].MeV()); + } + + // Multiply background rate with livetime and solid angle + value *= obs->livetime() * m_skymap.solidangle(i); + + // Subtract background value from sky map + m_skymap(i,0) -= value; + + } // endfor: looped over sky map pixels + + // Return + return; +} diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index ecfd07b..e62b3ae 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -28,6 +28,7 @@ #define CTSKYMAP_HPP /* __ Includes ___________________________________________________________ */ +#include #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ @@ -69,13 +70,16 @@ 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); // User parameters - GFilename m_outmap; //!< Output file name - double m_emin; //!< Minimum energy (TeV) - double m_emax; //!< Maximum energy (TeV) - bool m_publish; //!< Publish sky map? - GChatter m_chatter; //!< Chattiness + GFilename m_outmap; //!< Output file name + double m_emin; //!< Minimum energy (TeV) + double m_emax; //!< Maximum energy (TeV) + std::string m_bkgsubtract; //!< Background subtraction method + bool m_publish; //!< Publish sky map? + GChatter m_chatter; //!< Chattiness // Protected members GSkyMap m_skymap; //!< Sky map diff --git a/src/ctskymap/ctskymap.par b/src/ctskymap/ctskymap.par index 1604c3d..a8dfbaa 100644 --- a/src/ctskymap/ctskymap.par +++ b/src/ctskymap/ctskymap.par @@ -5,11 +5,10 @@ # ------------------------------------------------------------------------- # # # # File: ctskymap.par # -# Version: 1.1.0 # # # -# Author: Jurgen Knodlseder # -# knodlseder@cesr.fr # -# IRAP # +# Author: Juergen Knoedlseder # +# jknodlseder@irap.omp.eu # +# IRAP, Toulouse, France # # # # Purpose: Parameter file for the ctskymap application # # # @@ -33,10 +32,15 @@ usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" -coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +proj, s, a, CAR, AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" -proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" + +# +# Background subtraction parameters +#================================== +bkgsubtract, s, a, NONE, NONE|IRF,, "Background subtraction method" # # Standard parameters