Implement background model subtraction in ctskymap (#1882)

Jürgen Knödlseder
1 parent c1b24c0f
Showing 9 changed files with 232 additions and 52 deletions
ChangeLog
NEWS
doc/source/users/reference_manual/ctskymap.rst
doc/source/users/tutorials/1dc/first_skymap.png
doc/source/users/tutorials/1dc/first_skymap.rst
doc/source/users/tutorials/1dc/first_skymap_bkgsubtract.png
src/ctskymap/ctskymap.cpp
src/ctskymap/ctskymap.hpp
src/ctskymap/ctskymap.par
-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)
 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
 -----------------------------------
@@ -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 <CEL|GAL> [string]``
-    Coordinate system (CEL - celestial, GAL - galactic).
+    Coordinate system (CEL - celestial, GAL - galactic)
+``proj <AIT|AZP|CAR|MER|STG|TAN> [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 <AIT|AZP|CAR|MER|STG|TAN> [string]``
-    Projection method.
+``bkgsubtract <NONE|IRF> [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
@@ -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 <glossary_irf>`.
+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*
@@ -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<GChatter>(2);
+    m_emin        = 0.0;
+    m_emax        = 0.0;
+    m_bkgsubtract = "NONE";
+    m_publish     = false;
+    m_chatter     = static_cast<GChatter>(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<GChatter>((*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<GChatter>((*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<GCTAEventList*>
+                            (const_cast<GEvents*>(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<const GCTAEventList*>(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<GCTAEventList*>
-                                       (const_cast<GEvents*>(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<const GCTAResponseIrf*>(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<double> 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;
+}
@@ -28,6 +28,7 @@
 #define CTSKYMAP_HPP
 /* __ Includes ___________________________________________________________ */
+#include <string>
 #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
@@ -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,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 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