Add caldb and irf parameters to ctskymap, and store background and significance map

Jürgen Knödlseder
1 parent 05643325
Showing 5 changed files with 163 additions and 68 deletions
doc/source/users/reference_manual/ctskymap.rst
src/ctskymap/ctskymap.cpp
src/ctskymap/ctskymap.hpp
src/ctskymap/ctskymap.par
test/test_ctskymap.py
@@ -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 <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
+    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 <NONE|IRF> [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
@@ -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&amp; 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<GChatter>((*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<const GCTAResponseIrf*>(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;
@@ -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
 };
  
  
@@ -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"
  
 #
@@ -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