Added option to store errors in cttsmap

Michael Mayer
1 parent 2b3d392a
Showing 8 changed files with 110 additions and 46 deletions
cscripts/cstsmapmerge.py
cscripts/cstsmapsplit.par
src/cttsmap/cttsmap.cpp
src/cttsmap/cttsmap.hpp
src/cttsmap/cttsmap.par
test/data/tsmap_0.fits
test/data/tsmap_1.fits
test/test_cttsmap.py
@@ -46,7 +46,7 @@ class cstsmapmerge(ctools.cscript):
         self._files        = None
         self._in_filename  = ''
         self._tsmap        = gammalib.GSkyMap() # Empty sky map
-        self._statusmap    = gammalib.GSkyMap() # Empty sky map
+        self._progressmap  = gammalib.GSkyMap() # Empty sky map
         self._maps         = []
         self._mapnames     = []
         self._merged_files = []
@@ -136,11 +136,11 @@ class cstsmapmerge(ctools.cscript):
         # Open FITS file
         fits = gammalib.GFits(fitsfile)
-        # Read TS and status maps
+        # Read TS and progress maps
         self._tsmap     = gammalib.GSkyMap()
         self._tsmap.read(fits[0])
-        self._statusmap = gammalib.GSkyMap()
-        self._statusmap.read(fits["STATUS MAP"])
+        self._progressmap = gammalib.GSkyMap()
+        self._progressmap.read(fits["PROGRESS MAP"])
         # Get other maps 
         self._maps     = []
@@ -149,8 +149,8 @@ class cstsmapmerge(ctools.cscript):
         # Loop over extensions
         for hdu in fits:
-            # Leave out primary and status extension
-            if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP":
+            # Leave out primary and progress extension
+            if hdu.extname() != "IMAGE" and hdu.extname() != "PROGRESS MAP":
                 # Add present maps
                 skymap = gammalib.GSkyMap()
@@ -174,11 +174,11 @@ class cstsmapmerge(ctools.cscript):
         # Open FITS file
         fits = gammalib.GFits(fitsfile)
-        # Read TS and status maps
+        # Read TS and progress maps
         add_tsmap     = gammalib.GSkyMap()
         add_tsmap.read(fits[0])
-        add_statusmap = gammalib.GSkyMap()
-        add_statusmap.read(fits["STATUS MAP"])
+        add_progressmap = gammalib.GSkyMap()
+        add_progressmap.read(fits["PROGRESS MAP"])
         # Get other maps 
         add_maps = []
@@ -186,8 +186,8 @@ class cstsmapmerge(ctools.cscript):
         # Loop over extensions
         for hdu in fits:
-            # Leave out primary and status extension
-            if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP":
+            # Leave out primary and progress extension
+            if hdu.extname() != "IMAGE" and hdu.extname() != "PROGRESS MAP":
                 # Add present maps
                 skymap = gammalib.GSkyMap()
@@ -208,10 +208,10 @@ class cstsmapmerge(ctools.cscript):
         for i in range(self._tsmap.npix()):
             # Consider only bins that have been computed
-            if add_statusmap[i] > 0.5:
+            if add_progressmap[i] > 0.5:
                 # Raise exception if this bin has already been computed
-                if self._statusmap[i] > 0.5 and not self._overwrite:
+                if self._progressmap[i] > 0.5 and not self._overwrite:
                     msg = "Attempt to merge bin which apparently has "+\
                           "already been merged. File \""+fitsfile+"\" "+\
                           "contains already merged bins. Set hidden "+\
@@ -221,8 +221,8 @@ class cstsmapmerge(ctools.cscript):
                 # Copy TS values
                 self._tsmap[i] = add_tsmap[i]
-                # Copy status 
-                self._statusmap[i] = add_statusmap[i]
+                # Copy progress 
+                self._progressmap[i] = add_progressmap[i]
                 # Loop over maps and copy entries
                 for j in range(len(self._maps)):
@@ -241,12 +241,12 @@ class cstsmapmerge(ctools.cscript):
         Returns:
             Number of pixels for which TS has been computed.
         """
-        # Get status map for this file
-        status = gammalib.GSkyMap(fitsfile+"[STATUS MAP]")
+        # Get progress map for this file
+        progress = gammalib.GSkyMap(fitsfile+"[PROGRESS MAP]")
         # Count number of pixels with TS status set
         count = 0
-        for pix in status:
+        for pix in progress:
             if pix > 0.5:
                 count += 1
@@ -275,7 +275,7 @@ class cstsmapmerge(ctools.cscript):
         file0              = ""
         self._merged_files = []
-        # Test files for the entry status map. Use the first one to appear
+        # Test files for the entry progress map. Use the first one to appear
         # useful
         for fitsfile in self._files:
@@ -290,8 +290,8 @@ class cstsmapmerge(ctools.cscript):
             # Open FITS file
             fits = gammalib.GFits(fitsfile)
-            # If file contains a status map then use it
-            if fits.contains("STATUS MAP"):
+            # If file contains a progress map then use it
+            if fits.contains("PROGRESS MAP"):
                 fits.close()
                 file0 = fitsfile
                 if self._logTerse():
@@ -308,14 +308,14 @@ class cstsmapmerge(ctools.cscript):
                 if self._logExplicit():
                     self._log(gammalib.parformat("Skip file"))
                     self._log(fitsfile)
-                    self._log(" (no \"STATUS MAP\" extension)\n")
+                    self._log(" (no \"PROGRESS MAP\" extension)\n")
                 continue
         # Signal if no suitable file was found
         if file0 == "":
             if self._logTerse():
                 self._log("None of the provided files seems to be a sliced "
-                          "TS map file (none has a \"STATUS MAP\" "
+                          "TS map file (none has a \"PROGRESS MAP\" "
                           "extension).\n") 
         # ... otherwise merge files
@@ -347,13 +347,13 @@ class cstsmapmerge(ctools.cscript):
                 # Open FITS file
                 fits = gammalib.GFits(fitsfile)
-                # Skip if file does not contain status map
-                if not fits.contains("STATUS MAP"):
+                # Skip if file does not contain progress map
+                if not fits.contains("PROGRESS MAP"):
                     fits.close()
                     if self._logExplicit():
                         self._log(gammalib.parformat("Skip file"))
                         self._log(fitsfile)
-                        self._log(" (no \"STATUS MAP\" extension)\n")
+                        self._log(" (no \"PROGRESS MAP\" extension)\n")
                     continue
                 # Close FITS file
@@ -410,15 +410,15 @@ class cstsmapmerge(ctools.cscript):
         # Check if map is fully done
         done = True
-        for pix in self._statusmap:
+        for pix in self._progressmap:
             if pix < 0.5:
                 done = False
                 break
-        # Write status map if we are not done yet
+        # Write progress map if we are not done yet
         if not done:
-            self._statusmap.write(fits)
-            fits[fits.size()-1].extname("STATUS MAP")
+            self._progressmap.write(fits)
+            fits[fits.size()-1].extname("PROGRESS MAP")
         # Save FITS file
         fits.saveto(outmap, self._clobber())
@@ -27,6 +27,7 @@ bkgcube,       f, a, NONE,,, &quot;Input background cube file (only needed for stacke
 caldb,         s, a, prod2,,, "Calibration database"
 irf,           s, a, South_0.5h,,, "Instrument response function"
 edisp,         b, h, no,,, "Apply energy dispersion?"
+errors,        b, h, no,,, "Compute/store parameter errors?"
 outmap,        f, a, tsmap.fits,,, "Output Test Statistic map file"
 bins_per_job,  i, a, 5,,, "Number of TS map bins per task"
 compute_null,  b, a, yes,,,"Compute null hypothesis first?"
@@ -315,6 +315,14 @@ void cttsmap::run(void)
     	// Optimize observation container
     	m_obs.optimize(*opt);
+    	// Compute errors if necessary
+    	if (m_errors) {
+    		m_obs.errors(*opt);
+    	}
+
+    	// Get status of optimization
+    	int status = opt->status();
+
     	// Retrieve the Likelihood value
     	double logL1 = -(opt->value());
@@ -326,6 +334,7 @@ void cttsmap::run(void)
     		log << " TS value ....: ";
             log << ts << std::endl;
     	}
+
         else if (logTerse()) {
     		log << "TS for bin number ";
             log << i;
@@ -345,12 +354,26 @@ void cttsmap::run(void)
     	// Extract fitted test source parameters
         if (m_testsource != NULL) {
             for (int j = 0; j < m_mapnames.size(); ++j) {
-                m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+
+            	// Set parameter value and error maps
+            	if (m_mapnames[j].substr(0,2) == "e_") {
+
+            		// Get Parameter name by removing the error prefix
+            		std::string parname = m_mapnames[j].substr(2, m_mapnames[j].size());
+            		m_maps[j](i) = (*testsource)[parname].error();
+
+            	}
+            	else {
+            		m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+            	}
             }
     	}
-    	// Set status of bin to true
-    	m_statusmap(i) = 1.0;
+        // Set Fit status of the bin
+        m_statusmap(i) = status;
+
+    	// Add bin to progress map
+    	m_progressmap(i) = 1.0;
     	// Remove model from container
     	models.remove(m_srcname);
@@ -417,10 +440,14 @@ void cttsmap::save(void)
             fits[i+1]->extname(m_mapnames[i]);
         }
+        // Write Fit status map
+        m_statusmap.write(fits);
+        fits[fits.size()-1]->extname("FITSTATUS");
+
         // Add computation log if not all bins are computed
         if (m_binmin != -1 || m_binmax != -1) {
-            m_statusmap.write(fits);
-            fits[m_mapnames.size()+1]->extname("STATUS MAP");
+            m_progressmap.write(fits);
+            fits[m_mapnames.size()+1]->extname("PROGRESS MAP");
         }
         // Save FITS file
@@ -482,6 +509,7 @@ void cttsmap::init_members(void)
     m_outmap.clear();
     m_apply_edisp = false;
     m_publish     = false;
+    m_errors      = false;
     // Initialise protected members
     m_obs.clear();
@@ -490,6 +518,7 @@ void cttsmap::init_members(void)
     m_logL0      = 0.0;
     m_tsmap.clear();
     m_statusmap.clear();
+    m_progressmap.clear();
     m_mapnames.clear();
     m_maps.clear();
     m_testsource = NULL;
@@ -511,16 +540,18 @@ void cttsmap::copy_members(const cttsmap&amp; app)
     m_outmap      = app.m_outmap;
     m_apply_edisp = app.m_apply_edisp;
     m_publish     = app.m_publish;
+    m_errors      = app.m_errors;
     // Copy protected members
-    m_binmin    = app.m_binmin;
-    m_binmax    = app.m_binmax;
-    m_logL0     = app.m_logL0;
-    m_obs       = app.m_obs;
-    m_tsmap     = app.m_tsmap;
-    m_mapnames  = app.m_mapnames;
-    m_maps      = app.m_maps;
-    m_statusmap = app.m_statusmap;
+    m_binmin      = app.m_binmin;
+    m_binmax      = app.m_binmax;
+    m_logL0       = app.m_logL0;
+    m_obs         = app.m_obs;
+    m_tsmap       = app.m_tsmap;
+    m_mapnames    = app.m_mapnames;
+    m_maps        = app.m_maps;
+    m_statusmap   = app.m_statusmap;
+    m_progressmap = app.m_progressmap;
     // Clone protected members
     m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL;
@@ -612,6 +643,9 @@ void cttsmap::get_parameters(void)
     // Create sky map based on task parameters
     GSkyMap map = create_map(m_obs);
+    // Check if errors should be computed
+    m_errors  = (*this)["errors"].boolean();
+
     // Initialise maps from user parameters
     init_maps(map);
@@ -650,9 +684,11 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
     // Initialise map information
 	m_statusmap.clear();
+	m_progressmap.clear();
-	// Create status map
-	m_statusmap = GSkyMap(map);
+	// Create status map and progress map
+	m_statusmap   = GSkyMap(map);
+	m_progressmap = GSkyMap(map);
 	// Initialise maps of free parameters
     if (m_testsource != NULL) {
@@ -671,6 +707,16 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
             // Store parameter name
             m_mapnames.push_back((*m_testsource)[i].name());
+
+            // Add parameter errors if requested
+            if (m_errors) {
+
+				// Add sky map for fit error of the free parameter
+				m_maps.push_back(map);
+
+				// Store name of parameter error
+				m_mapnames.push_back("e_" + (*m_testsource)[i].name());
+            }
         } // endfor: looped over all model parameters
@@ -89,6 +89,7 @@ protected:
     GFilename                m_outmap;      //!< Output TS map file name
     bool                     m_apply_edisp; //!< Apply energy dispersion?
     bool                     m_publish;     //!< Publish TS map?
+    bool                     m_errors;      //!< Compute and store parameter errors?
     // Parameters to control speed and job splitting
     int                      m_binmin;      //!< Map bin number from which computation should start
@@ -98,7 +99,8 @@ protected:
     // Protected members
     GObservations            m_obs;         //!< Observation container
     GSkyMap                  m_tsmap;       //!< TS map
-    GSkyMap                  m_statusmap;   //!< Map of computed bins
+    GSkyMap                  m_statusmap;   //!< Map of optimizer fit status
+    GSkyMap                  m_progressmap; //!< Map of computed bins
     std::vector<std::string> m_mapnames;    //!< Names of free parameters
     std::vector<GSkyMap>     m_maps;        //!< Sky maps for each free parameter
     GModel*                  m_testsource;  //!< Pointer to test source for TS computation
@@ -27,6 +27,7 @@ bkgcube,   f, a, &quot;NONE&quot;,,, &quot;Input background cube file (only needed for stacked 
 caldb,     s, a, "prod2",,, "Calibration database"
 irf,       s, a, "South_0.5h",,, "Instrument response function"
 edisp,     b, h, no,,, "Apply energy dispersion?"
+errors,    b, h, no,,, "Compute/store parameter errors?"
 outmap,    f, a, "tsmap.fits",,, "Output Test Statistic map file"
 #
@@ -96,6 +96,20 @@ class Test(test):
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
              'Check invalid input file when executed from command line')
+        
+        # Setup cttsmap command
+        cmd = cttsmap+' inobs="'+self._events+'"'+ \
+                      ' inmodel="'+self._model+'" srcname="Crab"'+ \
+                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+                      ' errors=yes'+ \
+                      ' outmap="cttsmap_cmd1.fits"'+ \
+                      ' nxpix=5 nypix=5 binsz=0.02'+ \
+                      ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
+                      ' logfile="cttsmap_cmd1.log" chatter=1'
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
         # Return
         return