Commit c59d9b3616abe89dab6cd151160b8d7be3074612
Merge branch '1816-cttsmap-errors' into integration
# Conflicts: # src/cttsmap/cttsmap.cpp
Showing
10 changed files
with
235 additions
and
123 deletions
cscripts/cstsmapmerge.py
| ... | ... | @@ -30,17 +30,17 @@ import ctools |
| 30 | 30 | # ================== # |
| 31 | 31 | class cstsmapmerge(ctools.cscript): |
| 32 | 32 | """ |
| 33 | - Merge Test Statistic maps | |
| 33 | + Merge Test Statistic maps. | |
| 34 | 34 | """ |
| 35 | 35 | |
| 36 | 36 | # Constructor |
| 37 | 37 | def __init__(self, *argv): |
| 38 | 38 | """ |
| 39 | - Constructor | |
| 39 | + Constructor. | |
| 40 | 40 | """ |
| 41 | 41 | # Set name |
| 42 | 42 | self._name = 'cstsmapmerge' |
| 43 | - self._version = '1.2.0' | |
| 43 | + self._version = '1.1.0' | |
| 44 | 44 | |
| 45 | 45 | # Initialise class members |
| 46 | 46 | self._files = None |
| ... | ... | @@ -64,7 +64,7 @@ class cstsmapmerge(ctools.cscript): |
| 64 | 64 | # Private methods |
| 65 | 65 | def _get_parameters(self): |
| 66 | 66 | """ |
| 67 | - Get parameters from parfile | |
| 67 | + Get parameters from parfile. | |
| 68 | 68 | """ |
| 69 | 69 | # Get input maps string |
| 70 | 70 | inmaps = self['inmaps'].string() |
| ... | ... | @@ -118,7 +118,7 @@ class cstsmapmerge(ctools.cscript): |
| 118 | 118 | if self._read_ahead(): |
| 119 | 119 | self['outmap'].filename() |
| 120 | 120 | |
| 121 | - # Write input parameters into logger | |
| 121 | + # Write input parameters into logger | |
| 122 | 122 | self._log_parameters(gammalib.TERSE) |
| 123 | 123 | |
| 124 | 124 | # Return |
| ... | ... | @@ -126,12 +126,7 @@ class cstsmapmerge(ctools.cscript): |
| 126 | 126 | |
| 127 | 127 | def _init_ts_map(self, fitsfile): |
| 128 | 128 | """ |
| 129 | - Initialise Test Statistic map | |
| 130 | - | |
| 131 | - Parameters | |
| 132 | - ---------- | |
| 133 | - fitsfile : str | |
| 134 | - FITS file to be merged | |
| 129 | + Initialise Test Statistic map. | |
| 135 | 130 | """ |
| 136 | 131 | # Set filename |
| 137 | 132 | self._in_filename = fitsfile |
| ... | ... | @@ -169,12 +164,10 @@ class cstsmapmerge(ctools.cscript): |
| 169 | 164 | |
| 170 | 165 | def _merge_ts_map(self, fitsfile): |
| 171 | 166 | """ |
| 172 | - Merge TS map from FITS file into output TS map | |
| 167 | + Merge TS map from FITS file into output TS map. | |
| 173 | 168 | |
| 174 | - Parameters | |
| 175 | - ---------- | |
| 176 | - fitsfile : str | |
| 177 | - FITS file to be merged | |
| 169 | + Args: | |
| 170 | + fitsfile: FITS file to be merged. | |
| 178 | 171 | """ |
| 179 | 172 | # Open FITS file |
| 180 | 173 | fits = gammalib.GFits(fitsfile) |
| ... | ... | @@ -204,23 +197,23 @@ class cstsmapmerge(ctools.cscript): |
| 204 | 197 | |
| 205 | 198 | # Compare size of maps |
| 206 | 199 | if not len(add_maps) == len(self._maps): |
| 207 | - msg = ('Cannot merge map "%s" into map "%s" since the number of ' | |
| 208 | - 'parameters between both maps is different.' % | |
| 209 | - (fitsfile, self._in_filename)) | |
| 200 | + msg = 'Cannot merge map "'+fitsfile+'" into map "'+\ | |
| 201 | + self._in_filename+'" since the number of parameters ' +\ | |
| 202 | + 'between both maps is different.' | |
| 210 | 203 | raise RuntimeError(msg) |
| 211 | 204 | |
| 212 | 205 | # Loop over bins |
| 213 | 206 | for i in range(self._tsmap.npix()): |
| 214 | 207 | |
| 215 | 208 | # Consider only bins that have been computed |
| 216 | - if add_statusmap[i] > 0.5: | |
| 209 | + if add_statusmap[i] > -0.5: | |
| 217 | 210 | |
| 218 | 211 | # Raise exception if this bin has already been computed |
| 219 | - if self._statusmap[i] > 0.5 and not self._overwrite: | |
| 220 | - msg = ('Attempt to merge bin which apparently has already ' | |
| 221 | - 'been merged. File "%s" contains already merged ' | |
| 222 | - 'bins. Set hidden parameter "overwrite=yes" to ' | |
| 223 | - 'avoid this error.' % fitsfile) | |
| 212 | + if self._statusmap[i] > -0.5 and not self._overwrite: | |
| 213 | + msg = 'Attempt to merge bin which apparently has '+\ | |
| 214 | + 'already been merged. File "'+fitsfile+'" '+\ | |
| 215 | + 'contains already merged bins. Set hidden '+\ | |
| 216 | + 'parameter "overwrite=yes" to avoid this error.' | |
| 224 | 217 | raise RuntimeError(msg) |
| 225 | 218 | |
| 226 | 219 | # Copy TS values |
| ... | ... | @@ -238,17 +231,13 @@ class cstsmapmerge(ctools.cscript): |
| 238 | 231 | |
| 239 | 232 | def _get_number_of_ts_pixels(self, fitsfile): |
| 240 | 233 | """ |
| 241 | - Return number of pixels with TS values | |
| 234 | + Return number of pixels with TS values. | |
| 242 | 235 | |
| 243 | - Parameters | |
| 244 | - ---------- | |
| 245 | - fitsfile : str | |
| 246 | - FITS file to be merged | |
| 236 | + Args: | |
| 237 | + fitsfile: FITS file to be merged. | |
| 247 | 238 | |
| 248 | - Returns | |
| 249 | - ------- | |
| 250 | - counts : int | |
| 251 | - Number of pixels for which TS has been computed | |
| 239 | + Returns: | |
| 240 | + Number of pixels for which TS has been computed. | |
| 252 | 241 | """ |
| 253 | 242 | # Get status map for this file |
| 254 | 243 | status = gammalib.GSkyMap(fitsfile+'[STATUS MAP]') |
| ... | ... | @@ -256,7 +245,7 @@ class cstsmapmerge(ctools.cscript): |
| 256 | 245 | # Count number of pixels with TS status set |
| 257 | 246 | count = 0 |
| 258 | 247 | for pix in status: |
| 259 | - if pix > 0.5: | |
| 248 | + if pix > -0.5: | |
| 260 | 249 | count += 1 |
| 261 | 250 | |
| 262 | 251 | # Return number of pixels |
| ... | ... | @@ -266,7 +255,7 @@ class cstsmapmerge(ctools.cscript): |
| 266 | 255 | # Public methods |
| 267 | 256 | def run(self): |
| 268 | 257 | """ |
| 269 | - Run the script | |
| 258 | + Run the script. | |
| 270 | 259 | """ |
| 271 | 260 | # Switch screen logging on in debug mode |
| 272 | 261 | if self._logDebug(): |
| ... | ... | @@ -276,9 +265,7 @@ class cstsmapmerge(ctools.cscript): |
| 276 | 265 | self._get_parameters() |
| 277 | 266 | |
| 278 | 267 | # Write header into logger |
| 279 | - if self._logTerse(): | |
| 280 | - self._log('\n') | |
| 281 | - self._log.header1('Merge TS maps') | |
| 268 | + self._log_header1(gammalib.TERSE, 'Merge TS maps') | |
| 282 | 269 | |
| 283 | 270 | # Initialise file to start with |
| 284 | 271 | file0 = '' |
| ... | ... | @@ -290,8 +277,11 @@ class cstsmapmerge(ctools.cscript): |
| 290 | 277 | |
| 291 | 278 | # Skip file if it's not a FITS file |
| 292 | 279 | if not gammalib.GFilename(fitsfile).is_fits(): |
| 293 | - self._log_value(gammalib.EXPLICIT, 'Skip file', | |
| 294 | - '%s (not a FITS file)' % fitsfile) | |
| 280 | + | |
| 281 | + # Log message | |
| 282 | + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)') | |
| 283 | + | |
| 284 | + # Continue | |
| 295 | 285 | continue |
| 296 | 286 | |
| 297 | 287 | # Open FITS file |
| ... | ... | @@ -299,29 +289,40 @@ class cstsmapmerge(ctools.cscript): |
| 299 | 289 | |
| 300 | 290 | # If file contains a status map then use it |
| 301 | 291 | if fits.contains('STATUS MAP'): |
| 292 | + | |
| 293 | + # Close FITS file | |
| 302 | 294 | fits.close() |
| 295 | + | |
| 296 | + # Set filename to initial filename | |
| 303 | 297 | file0 = fitsfile |
| 304 | - if self._logTerse(): | |
| 305 | - self._log(gammalib.parformat('Initial TS map file')) | |
| 306 | - self._log(fitsfile) | |
| 307 | - count = self._get_number_of_ts_pixels(fitsfile) | |
| 308 | - self._log(' (%d TS pixels computed)' % count) | |
| 309 | - self._log('\n') | |
| 298 | + | |
| 299 | + # Log message | |
| 300 | + pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile) | |
| 301 | + self._log_value(gammalib.TERSE, 'Initial TS map file', fitsfile + pix_info) | |
| 302 | + | |
| 303 | + # Leave loop | |
| 310 | 304 | break |
| 311 | 305 | |
| 312 | 306 | # ... otherwise signal that file is useless |
| 313 | 307 | else: |
| 314 | - fits.close() | |
| 315 | - self._log_value(gammalib.EXPLICIT, 'Skip file', | |
| 316 | - '%s (no "STATUS MAP" extension)' % fitsfile) | |
| 308 | + | |
| 309 | + # Close fits file | |
| 310 | + fits.close() | |
| 311 | + | |
| 312 | + # Info message | |
| 313 | + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)') | |
| 314 | + | |
| 315 | + # Continue | |
| 317 | 316 | continue |
| 318 | 317 | |
| 319 | 318 | # Signal if no suitable file was found |
| 320 | 319 | if file0 == '': |
| 321 | - if self._logTerse(): | |
| 322 | - self._log('None of the provided files seems to be a sliced ' | |
| 323 | - 'TS map file (none has a "STATUS MAP" extension).\n') | |
| 324 | - | |
| 320 | + msg = 'None of the provided files seems to be a sliced ' + \ | |
| 321 | + 'TS map file (none has a "STATUS MAP" ' + \ | |
| 322 | + 'extension)' | |
| 323 | + raise RuntimeError(msg) | |
| 324 | + | |
| 325 | + | |
| 325 | 326 | # ... otherwise merge files |
| 326 | 327 | else: |
| 327 | 328 | |
| ... | ... | @@ -342,9 +343,8 @@ class cstsmapmerge(ctools.cscript): |
| 342 | 343 | |
| 343 | 344 | # Skip if file is not FITS |
| 344 | 345 | if not gammalib.GFilename(fitsfile).is_fits(): |
| 345 | - self._log_value(gammalib.EXPLICIT, 'Skip file', | |
| 346 | - '%s (not a FITS file)' % fitsfile) | |
| 347 | - continue | |
| 346 | + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)') | |
| 347 | + continue | |
| 348 | 348 | |
| 349 | 349 | # Open FITS file |
| 350 | 350 | fits = gammalib.GFits(fitsfile) |
| ... | ... | @@ -352,20 +352,15 @@ class cstsmapmerge(ctools.cscript): |
| 352 | 352 | # Skip if file does not contain status map |
| 353 | 353 | if not fits.contains('STATUS MAP'): |
| 354 | 354 | fits.close() |
| 355 | - self._log_value(gammalib.EXPLICIT, 'Skip file', | |
| 356 | - '%s (no "STATUS MAP" extension)' % fitsfile) | |
| 355 | + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)') | |
| 357 | 356 | continue |
| 358 | 357 | |
| 359 | 358 | # Close FITS file |
| 360 | 359 | fits.close() |
| 361 | 360 | |
| 362 | 361 | # Logging |
| 363 | - if self._logTerse(): | |
| 364 | - self._log(gammalib.parformat('Merge TS map file')) | |
| 365 | - self._log(fitsfile) | |
| 366 | - count = self._get_number_of_ts_pixels(fitsfile) | |
| 367 | - self._log(' (%d TS pixels computed)' % count) | |
| 368 | - self._log('\n') | |
| 362 | + pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile) | |
| 363 | + self._log_value(gammalib.TERSE, 'Merge TS map file', fitsfile + pix_info) | |
| 369 | 364 | |
| 370 | 365 | # Merge TS map |
| 371 | 366 | self._merge_ts_map(fitsfile) |
| ... | ... | @@ -378,16 +373,11 @@ class cstsmapmerge(ctools.cscript): |
| 378 | 373 | |
| 379 | 374 | def save(self): |
| 380 | 375 | """ |
| 381 | - Save TS map and remove slices if requested | |
| 376 | + Save TS map and remove slices if requested. | |
| 382 | 377 | """ |
| 383 | - # Write header | |
| 384 | - self._log_header1(gammalib.TERSE, 'Save TS map') | |
| 385 | 378 | |
| 386 | 379 | # Get output filename in case it was not read ahead |
| 387 | 380 | outmap = self['outmap'].filename() |
| 388 | - | |
| 389 | - # Log filename | |
| 390 | - self._log_value(gammalib.TERSE, 'TS map file', outmap.url()) | |
| 391 | 381 | |
| 392 | 382 | # Create FITS file |
| 393 | 383 | fits = gammalib.GFits() |
| ... | ... | @@ -403,33 +393,65 @@ class cstsmapmerge(ctools.cscript): |
| 403 | 393 | for i in range(len(self._mapnames)): |
| 404 | 394 | fits[i+1].extname(self._mapnames[i]) |
| 405 | 395 | |
| 396 | + # Treat status map separately | |
| 397 | + self._statusmap.write(fits) | |
| 398 | + fits[fits.size()-1].extname('STATUS MAP') | |
| 399 | + | |
| 406 | 400 | # Check if map is fully done |
| 407 | - done = True | |
| 401 | + bins_merged = 0 | |
| 408 | 402 | for pix in self._statusmap: |
| 409 | - if pix < 0.5: | |
| 410 | - done = False | |
| 411 | - break | |
| 403 | + if pix > -0.5: | |
| 404 | + bins_merged += 1 | |
| 405 | + | |
| 406 | + # Set boolean if map is fully done | |
| 407 | + done = (bins_merged == self._tsmap.npix()) | |
| 412 | 408 | |
| 413 | - # Write status map if we are not done yet | |
| 409 | + # Log summary header | |
| 410 | + self._log_header2(gammalib.TERSE, 'Merging Summary') | |
| 411 | + | |
| 412 | + # Write computing status if we are not done yet | |
| 414 | 413 | if not done: |
| 415 | - self._statusmap.write(fits) | |
| 416 | - fits[fits.size()-1].extname('STATUS MAP') | |
| 414 | + | |
| 415 | + # Log bins that were computed and those that were missing | |
| 416 | + self._log_value(gammalib.TERSE, 'TS map bins', self._tsmap.npix()) | |
| 417 | + self._log_value(gammalib.TERSE, 'Bins computed and merged', bins_merged ) | |
| 418 | + self._log_value(gammalib.TERSE, 'Bins missing', self._tsmap.npix() - bins_merged) | |
| 419 | + | |
| 420 | + else: | |
| 421 | + | |
| 422 | + # Write success message | |
| 423 | + self._log_string(gammalib.TERSE, 'TS map was fully computed and successfully merged') | |
| 424 | + | |
| 425 | + # Write header | |
| 426 | + self._log_header1(gammalib.TERSE, 'Save TS map') | |
| 427 | + | |
| 428 | + # Log filename | |
| 429 | + self._log_value(gammalib.TERSE, 'TS map file', outmap.url()) | |
| 417 | 430 | |
| 418 | 431 | # Save FITS file |
| 419 | - fits.saveto(outmap, self._clobber()) | |
| 432 | + fits.saveto(outmap, self._clobber()) | |
| 420 | 433 | |
| 421 | - # Delete TS input maps if requested | |
| 422 | - if self._delete: | |
| 434 | + # Delete TS input maps if requested and map is fully done | |
| 435 | + if self._delete and done: | |
| 436 | + | |
| 437 | + # Log header | |
| 438 | + self._log_header2(gammalib.TERSE, 'Delete slices TS map files') | |
| 439 | + | |
| 440 | + # Loop over sliced files | |
| 423 | 441 | for filename in self._merged_files: |
| 442 | + | |
| 443 | + # Remove file | |
| 424 | 444 | os.remove(filename) |
| 425 | - self._log_value(gammalib.TERSE, 'Deleted input file', filename) | |
| 445 | + | |
| 446 | + # Log message about file | |
| 447 | + self._log_value(gammalib.TERSE,'Deleted input file', filename) | |
| 426 | 448 | |
| 427 | 449 | # Return |
| 428 | 450 | return |
| 429 | 451 | |
| 430 | 452 | def execute(self): |
| 431 | 453 | """ |
| 432 | - Execute the script | |
| 454 | + Execute the script. | |
| 433 | 455 | """ |
| 434 | 456 | # Open logfile |
| 435 | 457 | self.logFileOpen() | ... | ... |
cscripts/cstsmapsplit.par
| ... | ... | @@ -27,6 +27,7 @@ bkgcube, f, a, NONE,,, "Input background cube file (only needed for stacke |
| 27 | 27 | caldb, s, a, prod2,,, "Calibration database" |
| 28 | 28 | irf, s, a, South_0.5h,,, "Instrument response function" |
| 29 | 29 | edisp, b, h, no,,, "Apply energy dispersion?" |
| 30 | +errors, b, h, no,,, "Compute/store parameter errors?" | |
| 30 | 31 | outmap, f, a, tsmap.fits,,, "Output Test Statistic map file" |
| 31 | 32 | bins_per_job, i, a, 5,,, "Number of TS map bins per task" |
| 32 | 33 | compute_null, b, a, yes,,,"Compute null hypothesis first?" | ... | ... |
doc/source/users/reference_manual/cttsmap.rst
| ... | ... | @@ -21,7 +21,10 @@ Gaussian sigma). |
| 21 | 21 | |
| 22 | 22 | cttsmap generates a FITS file comprising a sky map of TS values followed by |
| 23 | 23 | one extension per free parameter that contains sky maps of the fitted |
| 24 | -parameter values. | |
| 24 | +parameter values. To save computation time, the errors are not computed by | |
| 25 | +default. Specifying the hidden parameter ``errors=yes`` will switch on the | |
| 26 | +error computation. Accordingly, for each free parameter, the errors are also | |
| 27 | +stored in a separate sky map. | |
| 25 | 28 | |
| 26 | 29 | |
| 27 | 30 | General parameters |
| ... | ... | @@ -56,6 +59,9 @@ General parameters |
| 56 | 59 | |
| 57 | 60 | ``outmap [file]`` |
| 58 | 61 | Output Test Statistic map file. |
| 62 | + | |
| 63 | +``(errors = no) [boolean]`` | |
| 64 | + Compute and store parameter errors? | |
| 59 | 65 | |
| 60 | 66 | ``(usepnt = no) [boolean]`` |
| 61 | 67 | Use CTA pointing direction for map centre instead of xref/yref parameters? | ... | ... |
src/cttsmap/cttsmap.cpp
| ... | ... | @@ -173,9 +173,7 @@ cttsmap& cttsmap::operator=(const cttsmap& app) |
| 173 | 173 | ==========================================================================*/ |
| 174 | 174 | |
| 175 | 175 | /***********************************************************************//** |
| 176 | - * @brief Clear cttsmap tool | |
| 177 | - * | |
| 178 | - * Clears cttsmap tool. | |
| 176 | + * @brief Clear instance | |
| 179 | 177 | ***************************************************************************/ |
| 180 | 178 | void cttsmap::clear(void) |
| 181 | 179 | { |
| ... | ... | @@ -194,9 +192,6 @@ void cttsmap::clear(void) |
| 194 | 192 | this->ctlikelihood::init_members(); |
| 195 | 193 | init_members(); |
| 196 | 194 | |
| 197 | - // Write header into logger | |
| 198 | - log_header(); | |
| 199 | - | |
| 200 | 195 | // Return |
| 201 | 196 | return; |
| 202 | 197 | } |
| ... | ... | @@ -290,6 +285,14 @@ void cttsmap::run(void) |
| 290 | 285 | // Optimize observation container |
| 291 | 286 | m_obs.optimize(m_opt); |
| 292 | 287 | |
| 288 | + // Compute errors if necessary | |
| 289 | + if (m_errors) { | |
| 290 | + m_obs.errors(m_opt); | |
| 291 | + } | |
| 292 | + | |
| 293 | + // Get status of optimization | |
| 294 | + int status = m_opt.status(); | |
| 295 | + | |
| 293 | 296 | // Retrieve the Likelihood value |
| 294 | 297 | double logL1 = -(m_opt.value()); |
| 295 | 298 | |
| ... | ... | @@ -314,11 +317,27 @@ void cttsmap::run(void) |
| 314 | 317 | |
| 315 | 318 | // Extract fitted test source parameters |
| 316 | 319 | for (int j = 0; j < m_mapnames.size(); ++j) { |
| 317 | - m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); | |
| 318 | - } | |
| 319 | 320 | |
| 320 | - // Set status of bin to true | |
| 321 | - m_statusmap(i) = 1.0; | |
| 321 | + // If map name start with "e_" then set the parameter error ... | |
| 322 | + if (m_mapnames[j].substr(0,2) == "e_") { | |
| 323 | + | |
| 324 | + // Get parameter name by removing the error prefix | |
| 325 | + std::string parname = m_mapnames[j].substr(2, m_mapnames[j].size()); | |
| 326 | + | |
| 327 | + // Get parameter error | |
| 328 | + m_maps[j](i) = (*testsource)[parname].error(); | |
| 329 | + | |
| 330 | + } | |
| 331 | + | |
| 332 | + // ... otherwise set the parameter value | |
| 333 | + else { | |
| 334 | + m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); | |
| 335 | + } | |
| 336 | + | |
| 337 | + } // endfor: looped over all maps | |
| 338 | + | |
| 339 | + // Set Fit status of the bin | |
| 340 | + m_statusmap(i) = status; | |
| 322 | 341 | |
| 323 | 342 | // Remove model from container |
| 324 | 343 | models.remove(m_srcname); |
| ... | ... | @@ -371,11 +390,9 @@ void cttsmap::save(void) |
| 371 | 390 | fits[i+1]->extname(m_mapnames[i]); |
| 372 | 391 | } |
| 373 | 392 | |
| 374 | - // Add computation log if not all bins are computed | |
| 375 | - if (m_binmin != -1 || m_binmax != -1) { | |
| 376 | - m_statusmap.write(fits); | |
| 377 | - fits[m_mapnames.size()+1]->extname("STATUS MAP"); | |
| 378 | - } | |
| 393 | + // Write Fit status map | |
| 394 | + m_statusmap.write(fits); | |
| 395 | + fits[fits.size()-1]->extname("STATUS MAP"); | |
| 379 | 396 | |
| 380 | 397 | // Save FITS file |
| 381 | 398 | fits.saveto(m_outmap, clobber()); |
| ... | ... | @@ -437,6 +454,7 @@ void cttsmap::init_members(void) |
| 437 | 454 | m_outmap.clear(); |
| 438 | 455 | m_apply_edisp = false; |
| 439 | 456 | m_publish = false; |
| 457 | + m_errors = false; | |
| 440 | 458 | |
| 441 | 459 | // Initialise protected members |
| 442 | 460 | m_binmin = -1; |
| ... | ... | @@ -465,6 +483,7 @@ void cttsmap::copy_members(const cttsmap& app) |
| 465 | 483 | m_outmap = app.m_outmap; |
| 466 | 484 | m_apply_edisp = app.m_apply_edisp; |
| 467 | 485 | m_publish = app.m_publish; |
| 486 | + m_errors = app.m_errors; | |
| 468 | 487 | |
| 469 | 488 | // Copy protected members |
| 470 | 489 | m_binmin = app.m_binmin; |
| ... | ... | @@ -516,7 +535,7 @@ void cttsmap::get_parameters(void) |
| 516 | 535 | { |
| 517 | 536 | // Setup observations from "inobs" parameter |
| 518 | 537 | setup_observations(m_obs); |
| 519 | - | |
| 538 | + | |
| 520 | 539 | // Setup models from "inmodel" parameter |
| 521 | 540 | setup_models(m_obs, (*this)["srcname"].string()); |
| 522 | 541 | |
| ... | ... | @@ -545,6 +564,9 @@ void cttsmap::get_parameters(void) |
| 545 | 564 | // Create sky map based on task parameters |
| 546 | 565 | GSkyMap map = create_map(m_obs); |
| 547 | 566 | |
| 567 | + // Check if errors should be computed | |
| 568 | + m_errors = (*this)["errors"].boolean(); | |
| 569 | + | |
| 548 | 570 | // Initialise maps from user parameters |
| 549 | 571 | init_maps(map); |
| 550 | 572 | |
| ... | ... | @@ -585,14 +607,17 @@ void cttsmap::init_maps(const GSkyMap& map) |
| 585 | 607 | m_tsmap = GSkyMap(map); |
| 586 | 608 | |
| 587 | 609 | // Initialise map information |
| 588 | - m_statusmap.clear(); | |
| 610 | + m_statusmap.clear(); | |
| 611 | + | |
| 612 | + // Create status map and progress map | |
| 613 | + m_statusmap = GSkyMap(map); | |
| 589 | 614 | |
| 590 | - // Create status map | |
| 591 | - m_statusmap = GSkyMap(map); | |
| 615 | + // Set status map to invalid value -1 to signal that bin wasn't computed yet | |
| 616 | + m_statusmap = -1.0; | |
| 592 | 617 | |
| 593 | 618 | // Loop over all model parameters |
| 594 | 619 | for (int i = 0; i < m_testsource->size(); i++) { |
| 595 | - | |
| 620 | + | |
| 596 | 621 | // Skip fixed parameters |
| 597 | 622 | if ((*m_testsource)[i].is_fixed()) { |
| 598 | 623 | continue; |
| ... | ... | @@ -604,7 +629,18 @@ void cttsmap::init_maps(const GSkyMap& map) |
| 604 | 629 | |
| 605 | 630 | // Store parameter name |
| 606 | 631 | m_mapnames.push_back((*m_testsource)[i].name()); |
| 607 | - | |
| 632 | + | |
| 633 | + // Add parameter errors if requested | |
| 634 | + if (m_errors) { | |
| 635 | + | |
| 636 | + // Add sky map for fit error of the free parameter | |
| 637 | + m_maps.push_back(map); | |
| 638 | + | |
| 639 | + // Store name of parameter error | |
| 640 | + m_mapnames.push_back("e_" + (*m_testsource)[i].name()); | |
| 641 | + | |
| 642 | + } // endif: error calculation was requested | |
| 643 | + | |
| 608 | 644 | } // endfor: looped over all model parameters |
| 609 | 645 | |
| 610 | 646 | // Return | ... | ... |
src/cttsmap/cttsmap.hpp
| ... | ... | @@ -82,19 +82,20 @@ protected: |
| 82 | 82 | void init_maps(const GSkyMap& map); |
| 83 | 83 | |
| 84 | 84 | // User parameters |
| 85 | - std::string m_srcname; //!< Name of source which is moved around | |
| 86 | - GFilename m_outmap; //!< Output TS map file name | |
| 87 | - bool m_apply_edisp; //!< Apply energy dispersion? | |
| 88 | - bool m_publish; //!< Publish TS map? | |
| 85 | + std::string m_srcname; //!< Name of source which is moved around | |
| 86 | + GFilename m_outmap; //!< Output TS map file name | |
| 87 | + bool m_apply_edisp; //!< Apply energy dispersion? | |
| 88 | + bool m_publish; //!< Publish TS map? | |
| 89 | + bool m_errors; //!< Compute and store parameter errors? | |
| 89 | 90 | |
| 90 | 91 | // Parameters to control speed and job splitting |
| 91 | - int m_binmin; //!< Map bin number from which computation should start | |
| 92 | - int m_binmax; //!< Map bin number where map computation should end | |
| 93 | - double m_logL0; //!< Likelihood value of null hypothesis | |
| 92 | + int m_binmin; //!< Map bin number from which computation should start | |
| 93 | + int m_binmax; //!< Map bin number where map computation should end | |
| 94 | + double m_logL0; //!< Likelihood value of null hypothesis | |
| 94 | 95 | |
| 95 | 96 | // Protected members |
| 96 | 97 | GSkyMap m_tsmap; //!< TS map |
| 97 | - GSkyMap m_statusmap; //!< Map of computed bins | |
| 98 | + GSkyMap m_statusmap; //!< Map of optimizer fit status | |
| 98 | 99 | std::vector<std::string> m_mapnames; //!< Names of free parameters |
| 99 | 100 | std::vector<GSkyMap> m_maps; //!< Sky maps for each free parameter |
| 100 | 101 | GModel* m_testsource; //!< Pointer to test source for TS computation | ... | ... |
src/cttsmap/cttsmap.par
| ... | ... | @@ -5,11 +5,11 @@ |
| 5 | 5 | # ------------------------------------------------------------------------- # |
| 6 | 6 | # # |
| 7 | 7 | # File: cttsmap.par # |
| 8 | -# Version: 1.1.0 # | |
| 8 | +# Version: 1.2.0 # | |
| 9 | 9 | # # |
| 10 | 10 | # Author: Michael Mayer # |
| 11 | -# michael.mayer@desy.de # | |
| 12 | -# DESY # | |
| 11 | +# michael.mayer@physik.hu-berlin.de # | |
| 12 | +# HU Berlin # | |
| 13 | 13 | # # |
| 14 | 14 | # Purpose: Parameter file for the cttsmap application # |
| 15 | 15 | # # |
| ... | ... | @@ -30,6 +30,11 @@ edisp, b, h, no,,, "Apply energy dispersion?" |
| 30 | 30 | outmap, f, a, "tsmap.fits",,, "Output Test Statistic map file" |
| 31 | 31 | |
| 32 | 32 | # |
| 33 | +# TS map computation parameters | |
| 34 | +#============================== | |
| 35 | +errors, b, h, no,,, "Compute and store parameter errors?" | |
| 36 | + | |
| 37 | +# | |
| 33 | 38 | # Spatial binning parameters |
| 34 | 39 | #=========================== |
| 35 | 40 | usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" |
| ... | ... | @@ -43,7 +48,7 @@ proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" |
| 43 | 48 | |
| 44 | 49 | # |
| 45 | 50 | # Parameters for splitting and speed purpose |
| 46 | -#=========================== | |
| 51 | +#=========================================== | |
| 47 | 52 | binmin, i, h, -1,,, "First bin to compute" |
| 48 | 53 | binmax, i, h, -1,,, "Last bin to compute" |
| 49 | 54 | logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis" | ... | ... |
test/data/tsmap_0.fits
No preview for this file type
test/data/tsmap_1.fits
No preview for this file type
test/test_cstsmapmerge.py
| ... | ... | @@ -98,7 +98,16 @@ class Test(test): |
| 98 | 98 | # Check if execution failed |
| 99 | 99 | self.test_assert(self._execute(cmd) != 0, |
| 100 | 100 | 'Check invalid input file when executed from command line') |
| 101 | - | |
| 101 | + | |
| 102 | + # Setup cstsmapmerge command | |
| 103 | + cmd = cstsmapmerge+' inmaps="'+self._events+' '+self._model+'"'+ \ | |
| 104 | + ' outmap="cstsmapmerge_cmd3.fits"'+ \ | |
| 105 | + ' logfile="cstsmapmerge_cmd3.log" chatter=4' | |
| 106 | + | |
| 107 | + # Check if execution failed | |
| 108 | + self.test_assert(self._execute(cmd) != 0, | |
| 109 | + 'Check input file that are no TS maps when executed from command line') | |
| 110 | + | |
| 102 | 111 | # Return |
| 103 | 112 | return |
| 104 | 113 | |
| ... | ... | @@ -160,10 +169,10 @@ class Test(test): |
| 160 | 169 | |
| 161 | 170 | # Check merged TS map |
| 162 | 171 | self._check_result_file('cstsmapmerge_py4.fits') |
| 163 | - | |
| 164 | - # Set-up space-separated cstsmapmerge with invalid files | |
| 172 | + | |
| 173 | + # Set-up cstsmapmerge with incomplete file list | |
| 165 | 174 | merge = cscripts.cstsmapmerge() |
| 166 | - merge['inmaps'] = self._events+' '+self._model | |
| 175 | + merge['inmaps'] = self._map0 + ' ' + self._events | |
| 167 | 176 | merge['outmap'] = 'cstsmapmerge_py5.fits' |
| 168 | 177 | merge['logfile'] = 'cstsmapmerge_py5.log' |
| 169 | 178 | merge['chatter'] = 4 |
| ... | ... | @@ -172,13 +181,13 @@ class Test(test): |
| 172 | 181 | merge.execute() |
| 173 | 182 | |
| 174 | 183 | # Check merged TS map |
| 175 | - #self._check_result_file('cstsmapmerge_py5.fits') | |
| 184 | + self._check_result_file('cstsmapmerge_py5.fits', False) | |
| 176 | 185 | |
| 177 | 186 | # Return |
| 178 | 187 | return |
| 179 | 188 | |
| 180 | 189 | # Check result file |
| 181 | - def _check_result_file(self, filename): | |
| 190 | + def _check_result_file(self, filename, test_complete = True): | |
| 182 | 191 | """ |
| 183 | 192 | Check result file |
| 184 | 193 | """ |
| ... | ... | @@ -189,6 +198,7 @@ class Test(test): |
| 189 | 198 | ts = fits['Primary'] |
| 190 | 199 | prefactor = fits['Prefactor'] |
| 191 | 200 | index = fits['Index'] |
| 201 | + status = fits['STATUS MAP'] | |
| 192 | 202 | |
| 193 | 203 | # Check dimensions |
| 194 | 204 | self.test_value(ts.naxis(), 2, 'Check for 2 dimensions') |
| ... | ... | @@ -200,6 +210,23 @@ class Test(test): |
| 200 | 210 | self.test_value(index.naxis(), 2, 'Check for 2 dimensions') |
| 201 | 211 | self.test_value(index.naxes(0), 2, 'Check for 2 pixels in X') |
| 202 | 212 | self.test_value(index.naxes(1), 1, 'Check for 1 pixel in Y') |
| 213 | + self.test_value(status.naxis(), 2, 'Check for 2 dimensions') | |
| 214 | + self.test_value(status.naxes(0), 2, 'Check for 2 pixels in X') | |
| 215 | + self.test_value(status.naxes(1), 1, 'Check for 1 pixel in Y') | |
| 216 | + | |
| 217 | + # Initialise flag if map is complete | |
| 218 | + done = True | |
| 219 | + | |
| 220 | + # Loop over status sky map | |
| 221 | + for pix in status: | |
| 222 | + | |
| 223 | + # Check if pix is larger than threshold | |
| 224 | + if pix < -0.5: | |
| 225 | + done = False | |
| 226 | + break | |
| 227 | + | |
| 228 | + # Test for completeness of merged map | |
| 229 | + self.test_assert(done == test_complete, 'Test if map was merged completely') | |
| 203 | 230 | |
| 204 | 231 | # Return |
| 205 | 232 | return | ... | ... |
test/test_cttsmap.py
| ... | ... | @@ -97,6 +97,20 @@ class Test(test): |
| 97 | 97 | # Check if execution failed |
| 98 | 98 | self.test_assert(self._execute(cmd) != 0, |
| 99 | 99 | 'Check invalid input file when executed from command line') |
| 100 | + | |
| 101 | + # Setup cttsmap command | |
| 102 | + cmd = cttsmap+' inobs="'+self._events+'"'+ \ | |
| 103 | + ' inmodel="'+self._model+'" srcname="Crab"'+ \ | |
| 104 | + ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ | |
| 105 | + ' errors=yes'+ \ | |
| 106 | + ' outmap="cttsmap_cmd1.fits"'+ \ | |
| 107 | + ' nxpix=5 nypix=5 binsz=0.02'+ \ | |
| 108 | + ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ | |
| 109 | + ' logfile="cttsmap_cmd1.log" chatter=1' | |
| 110 | + | |
| 111 | + # Check if execution was successful | |
| 112 | + self.test_assert(self._execute(cmd) == 0, | |
| 113 | + 'Check successful execution from command line') | |
| 100 | 114 | |
| 101 | 115 | # Setup cttsmap --help option |
| 102 | 116 | cmd = cttsmap+' --help' | ... | ... |