diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py index 7dffc26..fb35361 100755 --- a/cscripts/cstsmapmerge.py +++ b/cscripts/cstsmapmerge.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()) diff --git a/cscripts/cstsmapsplit.par b/cscripts/cstsmapsplit.par index 4736f66..a54e0ec 100644 --- a/cscripts/cstsmapsplit.par +++ b/cscripts/cstsmapsplit.par @@ -27,6 +27,7 @@ bkgcube, f, a, NONE,,, "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?" diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index 27c3ff2..d772d79 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -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& 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& 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& 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 diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp index 49ff6cf..a17329f 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -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 m_mapnames; //!< Names of free parameters std::vector m_maps; //!< Sky maps for each free parameter GModel* m_testsource; //!< Pointer to test source for TS computation diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 6c23137..d3450d4 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -27,6 +27,7 @@ bkgcube, f, a, "NONE",,, "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" # diff --git a/test/data/tsmap_0.fits b/test/data/tsmap_0.fits index 57f530b..36f649a 100644 --- a/test/data/tsmap_0.fits +++ b/test/data/tsmap_0.fits diff --git a/test/data/tsmap_1.fits b/test/data/tsmap_1.fits index e905c5b..5359f0b 100644 --- a/test/data/tsmap_1.fits +++ b/test/data/tsmap_1.fits diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index 925694c..327410d 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -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