diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py index e17c24f..9322e43 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 = [] @@ -118,7 +118,11 @@ class cstsmapmerge(ctools.cscript): if self._read_ahead(): self['outmap'].filename() +<<<<<<< HEAD # Write input parameters into logger +======= + # Write input parameters into logger +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 self._log_parameters(gammalib.TERSE) # Return @@ -153,7 +157,11 @@ class cstsmapmerge(ctools.cscript): for hdu in fits: # Leave out primary and status extension +<<<<<<< HEAD if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP': +======= + if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP': +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Add present maps skymap = gammalib.GSkyMap() @@ -192,7 +200,11 @@ class cstsmapmerge(ctools.cscript): for hdu in fits: # Leave out primary and status extension +<<<<<<< HEAD if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP': +======= + if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP': +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Add present maps skymap = gammalib.GSkyMap() @@ -213,14 +225,22 @@ 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_statusmap[i] > -0.5: # Raise exception if this bin has already been computed +<<<<<<< HEAD if self._statusmap[i] > 0.5 and not self._overwrite: msg = ('Attempt to merge bin which apparently has already ' 'been merged. File "%s" contains already merged ' 'bins. Set hidden parameter "overwrite=yes" to ' 'avoid this error.' % fitsfile) +======= + 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." +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 raise RuntimeError(msg) # Copy TS values @@ -256,7 +276,7 @@ class cstsmapmerge(ctools.cscript): # Count number of pixels with TS status set count = 0 for pix in status: - if pix > 0.5: + if pix > -0.5: count += 1 # Return number of pixels @@ -290,8 +310,16 @@ class cstsmapmerge(ctools.cscript): # Skip file if it's not a FITS file if not gammalib.GFilename(fitsfile).is_fits(): +<<<<<<< HEAD self._log_value(gammalib.EXPLICIT, 'Skip file', '%s (not a FITS file)' % fitsfile) +======= + + # Log message + self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)") + + # Continue +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 continue # Open FITS file @@ -299,18 +327,35 @@ class cstsmapmerge(ctools.cscript): # If file contains a status map then use it if fits.contains('STATUS MAP'): +<<<<<<< HEAD +======= + + # Close FITS file +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 fits.close() + + # Set filename to initial filename file0 = fitsfile +<<<<<<< HEAD if self._logTerse(): self._log(gammalib.parformat('Initial TS map file')) self._log(fitsfile) count = self._get_number_of_ts_pixels(fitsfile) self._log(' (%d TS pixels computed)' % count) self._log('\n') +======= + + # 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) + + # Leave loop +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 break # ... otherwise signal that file is useless else: +<<<<<<< HEAD fits.close() self._log_value(gammalib.EXPLICIT, 'Skip file', '%s (no "STATUS MAP" extension)' % fitsfile) @@ -322,6 +367,24 @@ class cstsmapmerge(ctools.cscript): self._log('None of the provided files seems to be a sliced ' 'TS map file (none has a "STATUS MAP" extension).\n') +======= + + # Close fits file + fits.close() + + # Info message + 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") + +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # ... otherwise merge files else: @@ -342,9 +405,14 @@ class cstsmapmerge(ctools.cscript): # Skip if file is not FITS if not gammalib.GFilename(fitsfile).is_fits(): +<<<<<<< HEAD self._log_value(gammalib.EXPLICIT, 'Skip file', '%s (not a FITS file)' % fitsfile) continue +======= + self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)") + continue +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Open FITS file fits = gammalib.GFits(fitsfile) @@ -352,20 +420,29 @@ class cstsmapmerge(ctools.cscript): # Skip if file does not contain status map if not fits.contains('STATUS MAP'): fits.close() +<<<<<<< HEAD self._log_value(gammalib.EXPLICIT, 'Skip file', '%s (no "STATUS MAP" extension)' % fitsfile) +======= + self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)") +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 continue # Close FITS file fits.close() # Logging +<<<<<<< HEAD if self._logTerse(): self._log(gammalib.parformat('Merge TS map file')) self._log(fitsfile) count = self._get_number_of_ts_pixels(fitsfile) self._log(' (%d TS pixels computed)' % count) self._log('\n') +======= + 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) +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Merge TS map self._merge_ts_map(fitsfile) @@ -380,6 +457,7 @@ class cstsmapmerge(ctools.cscript): """ Save TS map and remove slices if requested """ +<<<<<<< HEAD # Write header self._log_header1(gammalib.TERSE, 'Save TS map') @@ -388,6 +466,11 @@ class cstsmapmerge(ctools.cscript): # Log filename self._log_value(gammalib.TERSE, 'TS map file', outmap.url()) +======= + + # Get output filename in case it was not read ahead + outmap = self["outmap"].filename() +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Create FITS file fits = gammalib.GFits() @@ -403,26 +486,67 @@ class cstsmapmerge(ctools.cscript): for i in range(len(self._mapnames)): fits[i+1].extname(self._mapnames[i]) + # Treat status map separately + self._statusmap.write(fits) + fits[fits.size()-1].extname('STATUS MAP') + # Check if map is fully done - done = True + bins_merged = 0 for pix in self._statusmap: - if pix < 0.5: - done = False - break + if pix > -0.5: + bins_merged += 1 + + # Set boolean if map is fully done + done = (bins_merged == self._tsmaps.npix()) + # Log summary header + self._log_header2(gammalib.TERSE, "Merging Summary") + # Write status map if we are not done yet if not done: +<<<<<<< HEAD self._statusmap.write(fits) fits[fits.size()-1].extname('STATUS MAP') +======= + + # 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) + + else: + + # Write success message + self._log_string(gammalib.TERSE, "TS map was fully computed and successfully merged") + + # Write header + self._log_header1(gammalib.TERSE, "Save TS map") + + # Log filename + self._log_value(gammalib.TERSE, "TS map file", outmap.url()) +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Save FITS file - fits.saveto(outmap, self._clobber()) + fits.saveto(outmap, self._clobber()) - # Delete TS input maps if requested - if self._delete: + # Delete TS input maps if requested and map is fully done + if self._delete and done: + + # Log header + self._log_header2(gammalib.TERSE, "Delete slices TS map files") + + # Loop over sliced files for filename in self._merged_files: + + # Remove file os.remove(filename) +<<<<<<< HEAD self._log_value(gammalib.TERSE, 'Deleted input file', filename) +======= + + # Log message about file + self._log_value(gammalib.TERSE,"Deleted input file", filename) +>>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5 # Return return 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 716a8e5..2fb45c0 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -215,6 +215,12 @@ 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; @@ -315,6 +321,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 +340,7 @@ void cttsmap::run(void) log << " TS value ....: "; log << ts << std::endl; } + else if (logTerse()) { log << "TS for bin number "; log << i; @@ -345,12 +360,23 @@ 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; // Remove model from container models.remove(m_srcname); @@ -417,11 +443,9 @@ void cttsmap::save(void) fits[i+1]->extname(m_mapnames[i]); } - // 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"); - } + // Write Fit status map + m_statusmap.write(fits); + fits[fits.size()-1]->extname("FITSTATUS"); // Save FITS file fits.saveto(m_outmap, clobber()); @@ -482,6 +506,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(); @@ -511,16 +536,17 @@ 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; // Clone protected members m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL; @@ -612,6 +638,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); @@ -654,9 +683,13 @@ void cttsmap::init_maps(const GSkyMap& map) // Initialise map information m_statusmap.clear(); - // Create status map + // Create status map and progress map m_statusmap = GSkyMap(map); + // Set status map to invalid value -1 to signal + // that bin wasn't computed yet + m_statusmap = -1.0; + // Initialise maps of free parameters if (m_testsource != NULL) { @@ -674,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..4db5078 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,7 @@ 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 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 e81ad49..ca9ebc3 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