diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py index 10ef14a..9b12c8c 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._statusmap = gammalib.GSkyMap() # Empty sky map self._maps = [] self._mapnames = [] self._merged_files = [] @@ -148,7 +148,7 @@ class cstsmapmerge(ctools.cscript): for hdu in fits: # Leave out primary and status extension - if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP': + if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP': # Add present maps skymap = gammalib.GSkyMap() @@ -185,7 +185,7 @@ class cstsmapmerge(ctools.cscript): for hdu in fits: # Leave out primary and status extension - if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP': + if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP': # Add present maps skymap = gammalib.GSkyMap() @@ -197,9 +197,9 @@ class cstsmapmerge(ctools.cscript): # Compare size of maps if not len(add_maps) == len(self._maps): - msg = "Cannot merge map \""+fitsfile+"\" into map \""+\ - self._in_filename+"\" since the number of parameters "\ - "between both maps is different." + msg = 'Cannot merge map "'+fitsfile+'" into map "'+\ + self._in_filename+'" since the number of parameters ' +\ + 'between both maps is different.' raise RuntimeError(msg) # Loop over bins @@ -210,10 +210,10 @@ class cstsmapmerge(ctools.cscript): # Raise exception if this bin has already been computed if self._statusmap[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 "+\ - "parameter \"overwrite=yes\" to avoid this error." + msg = 'Attempt to merge bin which apparently has '+\ + 'already been merged. File "'+fitsfile+'" '+\ + 'contains already merged bins. Set hidden '+\ + 'parameter "overwrite=yes" to avoid this error.' raise RuntimeError(msg) # Copy TS values @@ -240,7 +240,7 @@ class cstsmapmerge(ctools.cscript): Number of pixels for which TS has been computed. """ # Get status map for this file - status = gammalib.GSkyMap(fitsfile+"[STATUS MAP]") + status = gammalib.GSkyMap(fitsfile+'[STATUS MAP]') # Count number of pixels with TS status set count = 0 @@ -265,10 +265,10 @@ class cstsmapmerge(ctools.cscript): self._get_parameters() # Write header into logger - self._log_header1(gammalib.TERSE, "Merge TS maps") + self._log_header1(gammalib.TERSE, 'Merge TS maps') # Initialise file to start with - file0 = "" + file0 = '' self._merged_files = [] # Test files for the entry status map. Use the first one to appear @@ -279,7 +279,7 @@ class cstsmapmerge(ctools.cscript): if not gammalib.GFilename(fitsfile).is_fits(): # Log message - self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)") + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)') # Continue continue @@ -297,8 +297,8 @@ class cstsmapmerge(ctools.cscript): file0 = fitsfile # Log message - pix_info = " (%d TS pixels computed)" % self._get_number_of_ts_pixels(fitsfile) - self._log_value(gammalib.TERSE, "Initial TS map file", fitsfile + pix_info) + pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile) + self._log_value(gammalib.TERSE, 'Initial TS map file', fitsfile + pix_info) # Leave loop break @@ -310,16 +310,18 @@ class cstsmapmerge(ctools.cscript): fits.close() # Info message - self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)") + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)') # Continue continue # Signal if no suitable file was found - if file0 == "": - self._log_string(gammalib.TERSE, "None of the provided files seems to be a sliced " - "TS map file (none has a \"STATUS MAP\" " - "extension).\n") + if file0 == '': + msg = 'None of the provided files seems to be a sliced ' + \ + 'TS map file (none has a "STATUS MAP" ' + \ + 'extension)' + raise RuntimeError(msg) + # ... otherwise merge files else: @@ -341,7 +343,7 @@ class cstsmapmerge(ctools.cscript): # Skip if file is not FITS if not gammalib.GFilename(fitsfile).is_fits(): - self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)") + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)') continue # Open FITS file @@ -350,15 +352,15 @@ class cstsmapmerge(ctools.cscript): # Skip if file does not contain status map if not fits.contains('STATUS MAP'): fits.close() - self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)") + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)') continue # Close FITS file fits.close() # Logging - pix_info = " (%d TS pixels computed)" % self._get_number_of_ts_pixels(fitsfile) - self._log_value(gammalib.TERSE, "Merge TS map file", fitsfile + pix_info) + pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile) + self._log_value(gammalib.TERSE, 'Merge TS map file', fitsfile + pix_info) # Merge TS map self._merge_ts_map(fitsfile) @@ -375,7 +377,7 @@ class cstsmapmerge(ctools.cscript): """ # Get output filename in case it was not read ahead - outmap = self["outmap"].filename() + outmap = self['outmap'].filename() # Create FITS file fits = gammalib.GFits() @@ -402,29 +404,29 @@ class cstsmapmerge(ctools.cscript): bins_merged += 1 # Set boolean if map is fully done - done = (bins_merged == self._tsmaps.npix()) + done = (bins_merged == self._tsmap.npix()) # Log summary header - self._log_header2(gammalib.TERSE, "Merging Summary") + self._log_header2(gammalib.TERSE, 'Merging Summary') - # Write status map if we are not done yet + # Write computing status if we are not done yet if not done: # Log bins that were computed and those that were missing - self._log_value(gammalib.TERSE, "TS map bins", self._tsmap.npix()) - self._log_value(gammalib.TERSE, "Bins computed and merged", bins_merged ) - self._log_value(gammalib.TERSE, "Bins missing", self._tsmap.npix() - bins_merged) + self._log_value(gammalib.TERSE, 'TS map bins', self._tsmap.npix()) + self._log_value(gammalib.TERSE, 'Bins computed and merged', bins_merged ) + self._log_value(gammalib.TERSE, 'Bins missing', self._tsmap.npix() - bins_merged) else: # Write success message - self._log_string(gammalib.TERSE, "TS map was fully computed and successfully merged") + self._log_string(gammalib.TERSE, 'TS map was fully computed and successfully merged') # Write header - self._log_header1(gammalib.TERSE, "Save TS map") + self._log_header1(gammalib.TERSE, 'Save TS map') # Log filename - self._log_value(gammalib.TERSE, "TS map file", outmap.url()) + self._log_value(gammalib.TERSE, 'TS map file', outmap.url()) # Save FITS file fits.saveto(outmap, self._clobber()) @@ -433,7 +435,7 @@ class cstsmapmerge(ctools.cscript): if self._delete and done: # Log header - self._log_header2(gammalib.TERSE, "Delete slices TS map files") + self._log_header2(gammalib.TERSE, 'Delete slices TS map files') # Loop over sliced files for filename in self._merged_files: @@ -442,7 +444,7 @@ class cstsmapmerge(ctools.cscript): os.remove(filename) # Log message about file - self._log_value(gammalib.TERSE,"Deleted input file", filename) + self._log_value(gammalib.TERSE,'Deleted input file', filename) # Return return diff --git a/doc/source/users/reference_manual/cttsmap.rst b/doc/source/users/reference_manual/cttsmap.rst index aa26737..401125c 100644 --- a/doc/source/users/reference_manual/cttsmap.rst +++ b/doc/source/users/reference_manual/cttsmap.rst @@ -21,7 +21,10 @@ Gaussian sigma). cttsmap generates a FITS file comprising a sky map of TS values followed by one extension per free parameter that contains sky maps of the fitted -parameter values. +parameter values. To save computation time, the errors are not computed by +default. Specifying the hidden parameter ``errors=yes`` will switch on the +error computation. Accordingly, for each free parameter, the errors are also +stored in a separate sky map. General parameters @@ -56,6 +59,9 @@ General parameters ``outmap [file]`` Output Test Statistic map file. + +``(errors = no) [boolean]`` + Compute and store parameter errors? ``(usepnt = no) [boolean]`` Use CTA pointing direction for map centre instead of xref/yref parameters? diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index 2fb45c0..2e97e02 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -173,26 +173,20 @@ cttsmap& cttsmap::operator=(const cttsmap& app) ==========================================================================*/ /***********************************************************************//** - * @brief Clear cttsmap tool - * - * Clears cttsmap tool. + * @brief Clear instance ***************************************************************************/ void cttsmap::clear(void) { // Free members free_members(); this->ctool::free_members(); - - // Clear base class (needed to conserve tool name and version) - this->GApplication::clear(); + this->GApplication::free_members(); // Initialise members + this->GApplication::init_members(); this->ctool::init_members(); init_members(); - // Write header into logger - log_header(); - // Return return; } @@ -215,12 +209,9 @@ void cttsmap::run(void) // Get task parameters get_parameters(); -<<<<<<< HEAD -======= // Write parameters into logger log_parameters(TERSE); ->>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 // Set energy dispersion flag for all CTA observations and save old // values in save_edisp vector std::vector save_edisp; @@ -234,22 +225,16 @@ void cttsmap::run(void) } // Write observation(s) into logger - if (logTerse()) { - log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } - log << m_obs << std::endl; + if (m_obs.size() > 1) { + log_header1(TERSE, "Observations"); + } + else { + log_header1(TERSE, "Observation"); } + log_string(TERSE, m_obs.print()); // Write header - if (logTerse()) { - log << std::endl; - log.header1("Initialise TS map"); - } + log_header1(TERSE, "Initialise TS map"); // Get bins to be computed @@ -276,10 +261,7 @@ void cttsmap::run(void) if (m_logL0 == 0.0) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Compute NULL Hypothesis for TS computation"); - } + log_header1(TERSE, "Compute NULL Hypothesis for TS computation"); // Compute likelihood without the test source m_obs.models(models); @@ -289,10 +271,7 @@ void cttsmap::run(void) } // Write header - if (logTerse()) { - log << std::endl; - log.header1("Generate TS map"); - } + log_header1(TERSE, "Generate TS map"); // Loop over grid positions for (int i = binmin; i < binmax; ++i) { @@ -301,12 +280,9 @@ void cttsmap::run(void) GSkyDir bincentre = m_tsmap.inx2dir(i); // Header for verbose logging - if (logExplicit()) { - std::string msg = "Computing TS for bin number "+gammalib::str(i)+ - " at "+bincentre.print(); - log << std::endl; - log.header2(msg); - } + std::string msg = "Computing TS for bin number "+gammalib::str(i)+ + " at "+bincentre.print(); + log_header2(EXPLICIT, msg); // Add test source at current bin position if (m_testsource != NULL) { @@ -340,7 +316,6 @@ void cttsmap::run(void) log << " TS value ....: "; log << ts << std::endl; } - else if (logTerse()) { log << "TS for bin number "; log << i; @@ -445,7 +420,7 @@ void cttsmap::save(void) // Write Fit status map m_statusmap.write(fits); - fits[fits.size()-1]->extname("FITSTATUS"); + fits[fits.size()-1]->extname("STATUS MAP"); // Save FITS file fits.saveto(m_outmap, clobber()); @@ -465,10 +440,7 @@ void cttsmap::save(void) void cttsmap::publish(const std::string& name) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Publish TS map"); - } + log_header1(TERSE, "Publish TS map"); // Set default name is user name is empty std::string user_name(name); @@ -477,10 +449,7 @@ void cttsmap::publish(const std::string& name) } // Log filename - if (logTerse()) { - log << gammalib::parformat("Publish TS map"); - log << user_name << std::endl; - } + log_value(TERSE, "Publish TS map", user_name); // Publish TS map m_tsmap.publish(user_name); @@ -659,9 +628,6 @@ void cttsmap::get_parameters(void) m_outmap = (*this)["outmap"].filename(); } - // Write parameters into logger - log_parameters(TERSE); - // Return return; } @@ -716,7 +682,8 @@ void cttsmap::init_maps(const GSkyMap& map) // Store name of parameter error m_mapnames.push_back("e_" + (*m_testsource)[i].name()); - } + + } // endif: error calculation was requested } // endfor: looped over all model parameters diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index d3450d4..e2cfba9 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -27,7 +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?" +errors, b, h, no,,, "Compute and 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 36f649a..8aafb23 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 5359f0b..01141ae 100644 --- a/test/data/tsmap_1.fits +++ b/test/data/tsmap_1.fits diff --git a/test/test_cstsmapmerge.py b/test/test_cstsmapmerge.py index c1dc352..8359cdc 100755 --- a/test/test_cstsmapmerge.py +++ b/test/test_cstsmapmerge.py @@ -21,6 +21,8 @@ import gammalib import cscripts from testing import test +from __builtin__ import False +from Carbon.Aliases import false # ================================== # @@ -98,7 +100,16 @@ class Test(test): # Check if execution failed self.test_assert(self._execute(cmd) != 0, 'Check invalid input file when executed from command line') - + + # Setup cstsmapmerge command + cmd = cstsmapmerge+' inmaps="'+self._events+' '+self._model+'"'+ \ + ' outmap="cstsmapmerge_cmd3.fits"'+ \ + ' logfile="cstsmapmerge_cmd3.log" chatter=4' + + # Check if execution failed + self.test_assert(self._execute(cmd) != 0, + 'Check input file that are no TS maps when executed from command line') + # Return return @@ -160,10 +171,10 @@ class Test(test): # Check merged TS map self._check_result_file('cstsmapmerge_py4.fits') - - # Set-up space-separated cstsmapmerge with invalid files + + # Set-up cstsmapmerge with incomplete file list merge = cscripts.cstsmapmerge() - merge['inmaps'] = self._events+' '+self._model + merge['inmaps'] = self._map0 + ' ' + self._events merge['outmap'] = 'cstsmapmerge_py5.fits' merge['logfile'] = 'cstsmapmerge_py5.log' merge['chatter'] = 4 @@ -172,13 +183,13 @@ class Test(test): merge.execute() # Check merged TS map - #self._check_result_file('cstsmapmerge_py5.fits') + self._check_result_file('cstsmapmerge_py5.fits', False) # Return return # Check result file - def _check_result_file(self, filename): + def _check_result_file(self, filename, test_complete = True): """ Check result file """ @@ -189,6 +200,7 @@ class Test(test): ts = fits['Primary'] prefactor = fits['Prefactor'] index = fits['Index'] + status = fits['STATUS MAP'] # Check dimensions self.test_value(ts.naxis(), 2, 'Check for 2 dimensions') @@ -200,6 +212,23 @@ class Test(test): self.test_value(index.naxis(), 2, 'Check for 2 dimensions') self.test_value(index.naxes(0), 2, 'Check for 2 pixels in X') self.test_value(index.naxes(1), 1, 'Check for 1 pixel in Y') + self.test_value(status.naxis(), 2, 'Check for 2 dimensions') + self.test_value(status.naxes(0), 2, 'Check for 2 pixels in X') + self.test_value(status.naxes(1), 1, 'Check for 1 pixel in Y') + + # Initialise flag if map is complete + done = True + + # Loop over status sky map + for pix in status: + + # Check if pix is larger than threshold + if pix < -0.5: + done = False + break + + # Test for completeness of merged map + self.test_assert(done == test_complete, 'Test if map was merged completely') # Return return