Commit d04d0a7a690f92a05f881eda326215e9ea530475
1 parent
21212a24
Replace cspull ASCII output files by FITS files (#1707)
Showing
13 changed files
with
379 additions
and
65 deletions
ChangeLog
NEWS
| 1 | 1 | New Features and Important Changes in ctools 2.0.0 |
| 2 | 2 | |
| 3 | -11 May 2021 | |
| 3 | +19 May 2021 | |
| 4 | 4 | |
| 5 | 5 | |
| 6 | 6 | Introduction |
| ... | ... | @@ -255,7 +255,7 @@ the width, height and position angle of the rectangular sky regions (#3248). |
| 255 | 255 | |
| 256 | 256 | cspull - Generate pull distributions |
| 257 | 257 | ------------------------------------ |
| 258 | -None | |
| 258 | +The ASCII output file was replaced by a FITS file (#1707). | |
| 259 | 259 | |
| 260 | 260 | |
| 261 | 261 | csresmap - Generate residual map |
| ... | ... | @@ -361,8 +361,9 @@ regions (#3248). |
| 361 | 361 | |
| 362 | 362 | Examples |
| 363 | 363 | -------- |
| 364 | -The show_butterfly and show_sensitivity.py scripts were adapted to deal also | |
| 365 | -with FITS output files (#1707). | |
| 364 | +The show_butterfly, show_pull_evolution.py, show_pull_histogram.py and | |
| 365 | +show_sensitivity.py scripts were adapted to deal also with FITS output files | |
| 366 | +(#1707). | |
| 366 | 367 | |
| 367 | 368 | |
| 368 | 369 | caldb | ... | ... |
README.md
| 1 | 1 | ctools information |
| 2 | 2 | ================== |
| 3 | -* Version: 2.0.0.dev (11 May 2021) | |
| 3 | +* Version: 2.0.0.dev (19 May 2021) | |
| 4 | 4 | * GammaLib dependency: 2.0.0.dev |
| 5 | 5 | |
| 6 | 6 | [](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/) | ... | ... |
cscripts/Makefile.am
| 1 | 1 | ############################################################################# |
| 2 | 2 | # Makefile for cscripts scripts # |
| 3 | 3 | # ------------------------------------------------------------------------- # |
| 4 | -# Copyright (C) 2016-2020 Juergen Knoedlseder # | |
| 4 | +# Copyright (C) 2016-2021 Juergen Knoedlseder # | |
| 5 | 5 | # ------------------------------------------------------------------------- # |
| 6 | 6 | # # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify # |
| ... | ... | @@ -203,6 +203,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \ |
| 203 | 203 | $(top_srcdir)/test/data/pha.fits \ |
| 204 | 204 | $(top_srcdir)/test/data/pha_channels.fits \ |
| 205 | 205 | $(top_srcdir)/test/data/pull.dat \ |
| 206 | + $(top_srcdir)/test/data/pull.fits \ | |
| 206 | 207 | $(top_srcdir)/test/data/sensitivity.dat \ |
| 207 | 208 | $(top_srcdir)/test/data/sensitivity_integral.dat \ |
| 208 | 209 | $(top_srcdir)/test/data/sensitivity.fits \ | ... | ... |
cscripts/cspull.par
| ... | ... | @@ -24,7 +24,7 @@ caldb, s, a, prod2,,, "Calibration database" |
| 24 | 24 | irf, s, a, South_0.5h,,, "Instrument response function" |
| 25 | 25 | edisp, b, h, no,,, "Apply energy dispersion?" |
| 26 | 26 | deadc, r, h, 0.98,0,1, "Deadtime correction factor" |
| 27 | -outfile, f, a, pull.dat,,, "Output pull distribution file" | |
| 27 | +outfile, f, a, pull.fits,,, "Output pull distribution FITS file" | |
| 28 | 28 | |
| 29 | 29 | # |
| 30 | 30 | # Script parameters | ... | ... |
cscripts/cspull.py
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | # ========================================================================== |
| 3 | 3 | # This script generates the pull distribution for all model parameters. |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2011-2018 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2011-2021 Juergen Knoedlseder | |
| 6 | 6 | # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify |
| 8 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -41,7 +41,8 @@ class cspull(ctools.csobservation): |
| 41 | 41 | # Initialise application by calling the appropriate class constructor |
| 42 | 42 | self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) |
| 43 | 43 | |
| 44 | - # Initialise some members | |
| 44 | + # Initialise class members | |
| 45 | + self._fits = None | |
| 45 | 46 | self._nthreads = 0 |
| 46 | 47 | |
| 47 | 48 | # Return |
| ... | ... | @@ -59,6 +60,7 @@ class cspull(ctools.csobservation): |
| 59 | 60 | """ |
| 60 | 61 | # Set pickled dictionary |
| 61 | 62 | state = {'base' : ctools.csobservation.__getstate__(self), |
| 63 | + 'fits' : self._fits, | |
| 62 | 64 | 'nthreads' : self._nthreads} |
| 63 | 65 | |
| 64 | 66 | # Return pickled dictionary |
| ... | ... | @@ -74,11 +76,13 @@ class cspull(ctools.csobservation): |
| 74 | 76 | Pickled instance |
| 75 | 77 | """ |
| 76 | 78 | ctools.csobservation.__setstate__(self, state['base']) |
| 79 | + self._fits = state['fits'] | |
| 77 | 80 | self._nthreads = state['nthreads'] |
| 78 | 81 | |
| 79 | 82 | # Return |
| 80 | 83 | return |
| 81 | 84 | |
| 85 | + | |
| 82 | 86 | # Private methods |
| 83 | 87 | def _get_parameters(self): |
| 84 | 88 | """ |
| ... | ... | @@ -121,9 +125,12 @@ class cspull(ctools.csobservation): |
| 121 | 125 | self['chatter'].integer() |
| 122 | 126 | |
| 123 | 127 | # Query some parameters |
| 124 | - self['outfile'].filename() | |
| 125 | 128 | self['profile'].boolean() |
| 126 | 129 | |
| 130 | + # Read ahead output parameters | |
| 131 | + if self._read_ahead(): | |
| 132 | + self['outfile'].filename() | |
| 133 | + | |
| 127 | 134 | # Write input parameters into logger |
| 128 | 135 | self._log_parameters(gammalib.TERSE) |
| 129 | 136 | |
| ... | ... | @@ -297,13 +304,18 @@ class cspull(ctools.csobservation): |
| 297 | 304 | # Set name as a combination of model name and parameter |
| 298 | 305 | # name separated by an underscore. In that way each |
| 299 | 306 | # parameter has a unique name. |
| 300 | - name = model.name()+'_'+par.name() | |
| 307 | + name = model.name() + '_' + par.name() | |
| 308 | + | |
| 309 | + # Set "Parameter", "Parameter_error" and "Parameter_Pull" | |
| 310 | + # column names | |
| 311 | + col_par = name | |
| 312 | + col_par_error = name+'_error' | |
| 313 | + col_par_pull = name+'_pull' | |
| 301 | 314 | |
| 302 | - # Append parameter, Pull_parameter and e_parameter column | |
| 303 | - # names | |
| 304 | - colnames.append(name) | |
| 305 | - colnames.append('Pull_'+name) | |
| 306 | - colnames.append('e_'+name) | |
| 315 | + # Append column names | |
| 316 | + colnames.append(col_par) | |
| 317 | + colnames.append(col_par_error) | |
| 318 | + colnames.append(col_par_pull) | |
| 307 | 319 | |
| 308 | 320 | # Compute pull for this parameter as the difference |
| 309 | 321 | # (fitted - true) / error |
| ... | ... | @@ -317,9 +329,9 @@ class cspull(ctools.csobservation): |
| 317 | 329 | pull = 99.0 |
| 318 | 330 | |
| 319 | 331 | # Store results in dictionary |
| 320 | - values[name] = fitted_value | |
| 321 | - values['Pull_'+name] = pull | |
| 322 | - values['e_'+name] = error | |
| 332 | + values[col_par] = fitted_value | |
| 333 | + values[col_par_error] = error | |
| 334 | + values[col_par_pull] = pull | |
| 323 | 335 | |
| 324 | 336 | # Write results into logger |
| 325 | 337 | value = '%.4f (%e +/- %e)' % (pull, fitted_value, error) |
| ... | ... | @@ -331,6 +343,100 @@ class cspull(ctools.csobservation): |
| 331 | 343 | # Return |
| 332 | 344 | return result |
| 333 | 345 | |
| 346 | + def _create_fits(self, results): | |
| 347 | + """ | |
| 348 | + Create FITS file from results | |
| 349 | + | |
| 350 | + Parameters | |
| 351 | + ---------- | |
| 352 | + results : list of dict | |
| 353 | + List of result dictionaries | |
| 354 | + """ | |
| 355 | + # Gather headers for parameter columns | |
| 356 | + headers = [] | |
| 357 | + for colname in results[0]['colnames']: | |
| 358 | + if colname != 'LogL' and colname != 'Sim_Events' and \ | |
| 359 | + colname != 'Npred_Events': | |
| 360 | + headers.append(colname) | |
| 361 | + | |
| 362 | + # Create FITS table columns | |
| 363 | + nrows = len(results) | |
| 364 | + logl = gammalib.GFitsTableDoubleCol('LOGL', nrows) | |
| 365 | + nevents = gammalib.GFitsTableDoubleCol('NEVENTS', nrows) | |
| 366 | + npred = gammalib.GFitsTableDoubleCol('NPRED', nrows) | |
| 367 | + logl.unit('') | |
| 368 | + nevents.unit('counts') | |
| 369 | + npred.unit('counts') | |
| 370 | + columns = [] | |
| 371 | + for header in headers: | |
| 372 | + name = gammalib.toupper(header) | |
| 373 | + column = gammalib.GFitsTableDoubleCol(name, nrows) | |
| 374 | + column.unit('') | |
| 375 | + columns.append(column) | |
| 376 | + | |
| 377 | + # Fill FITS table columns | |
| 378 | + for i, result in enumerate(results): | |
| 379 | + logl[i] = result['values']['LogL'] | |
| 380 | + nevents[i] = result['values']['Sim_Events'] | |
| 381 | + npred[i] = result['values']['Npred_Events'] | |
| 382 | + for k, column in enumerate(columns): | |
| 383 | + column[i] = result['values'][headers[k]] | |
| 384 | + | |
| 385 | + # Initialise FITS Table with extension "PULL_DISTRIBUTION" | |
| 386 | + table = gammalib.GFitsBinTable(nrows) | |
| 387 | + table.extname('PULL_DISTRIBUTION') | |
| 388 | + | |
| 389 | + # Add keywors for compatibility with gammalib.GMWLSpectrum | |
| 390 | + table.card('INSTRUME', 'CTA', 'Name of Instrument') | |
| 391 | + table.card('TELESCOP', 'CTA', 'Name of Telescope') | |
| 392 | + | |
| 393 | + # Stamp header | |
| 394 | + self._stamp(table) | |
| 395 | + | |
| 396 | + # Set optional keyword values | |
| 397 | + if gammalib.toupper(self['onsrc'].string()) != 'NONE': | |
| 398 | + onrad = self['onrad'].real() | |
| 399 | + else: | |
| 400 | + onrad = 0.0 | |
| 401 | + if self['enumbins'].integer() > 0: | |
| 402 | + npix = self['npix'].integer() | |
| 403 | + binsz = self['binsz'].real() | |
| 404 | + coordsys = self['coordsys'].string() | |
| 405 | + proj = self['proj'].string() | |
| 406 | + else: | |
| 407 | + npix = 0 | |
| 408 | + binsz = 0.0 | |
| 409 | + coordsys = '' | |
| 410 | + proj = '' | |
| 411 | + | |
| 412 | + # Add script keywords | |
| 413 | + table.card('NPULLS', self['ntrials'].integer(), 'Number of pulls') | |
| 414 | + table.card('SEED', self['seed'].integer(), 'Seed value for pulls') | |
| 415 | + table.card('ONSRC', self['onsrc'].string(), 'Name of On surce for On/Off analysis') | |
| 416 | + table.card('ONRAD', onrad, '[deg] Radius of On region') | |
| 417 | + table.card('ENUMBINS', self['enumbins'].integer(), 'Number of energy bins') | |
| 418 | + table.card('NPIX', npix, 'Number of pixels for binned analysis') | |
| 419 | + table.card('BINSZ', binsz, 'Pixel size for binned analysis') | |
| 420 | + table.card('COORDSYS', coordsys, 'Coordinate system for binned analysis') | |
| 421 | + table.card('PROJ', proj, 'Projection for binned analysis') | |
| 422 | + table.card('STAT', self['statistic'].string(), 'Optimization statistic') | |
| 423 | + table.card('EDISP', self['edisp'].boolean(), 'Use energy dispersion?') | |
| 424 | + table.card('PROFILE', self['profile'].boolean(), 'Use likelihood profile method for errors?') | |
| 425 | + | |
| 426 | + # Append filled columns to fits table | |
| 427 | + table.append(logl) | |
| 428 | + table.append(nevents) | |
| 429 | + table.append(npred) | |
| 430 | + for column in columns: | |
| 431 | + table.append(column) | |
| 432 | + | |
| 433 | + # Create the FITS file now | |
| 434 | + self._fits = gammalib.GFits() | |
| 435 | + self._fits.append(table) | |
| 436 | + | |
| 437 | + # Return | |
| 438 | + return | |
| 439 | + | |
| 334 | 440 | |
| 335 | 441 | # Public methods |
| 336 | 442 | def run(self): |
| ... | ... | @@ -356,6 +462,9 @@ class cspull(ctools.csobservation): |
| 356 | 462 | # Get seed value |
| 357 | 463 | seed = self['seed'].integer() |
| 358 | 464 | |
| 465 | + # Initialise results | |
| 466 | + results = [] | |
| 467 | + | |
| 359 | 468 | # If more than a single thread is requested then use multiprocessing |
| 360 | 469 | if self._nthreads > 1: |
| 361 | 470 | args = [(self, '_trial', i + seed) for i in range(ntrials)] |
| ... | ... | @@ -367,16 +476,21 @@ class cspull(ctools.csobservation): |
| 367 | 476 | # If multiprocessing was used then recover results and put them |
| 368 | 477 | # into the log file |
| 369 | 478 | if self._nthreads > 1: |
| 370 | - result = poolresults[i][0] | |
| 479 | + #result = poolresults[i][0] | |
| 480 | + results.append(poolresults[i][0]) | |
| 371 | 481 | self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False) |
| 372 | 482 | |
| 373 | 483 | # ... otherwise make a trial now |
| 374 | 484 | else: |
| 485 | + | |
| 486 | + # Run trial | |
| 375 | 487 | result = self._trial(i + seed) |
| 376 | 488 | |
| 377 | - # Write out trial result | |
| 378 | - ioutils.write_csv_row(self['outfile'].filename().url(), i, | |
| 379 | - result['colnames'], result['values']) | |
| 489 | + # Append results | |
| 490 | + results.append(result) | |
| 491 | + | |
| 492 | + # Create FITS file | |
| 493 | + self._create_fits(results) | |
| 380 | 494 | |
| 381 | 495 | # Return |
| 382 | 496 | return |
| ... | ... | @@ -396,6 +510,38 @@ class cspull(ctools.csobservation): |
| 396 | 510 | # Return |
| 397 | 511 | return |
| 398 | 512 | |
| 513 | + def save(self): | |
| 514 | + """ | |
| 515 | + Save pull FITS file | |
| 516 | + """ | |
| 517 | + # Write header | |
| 518 | + self._log_header1(gammalib.TERSE, 'Save pull distribution') | |
| 519 | + | |
| 520 | + # Continue only if FITS file is valid | |
| 521 | + if self._fits != None: | |
| 522 | + | |
| 523 | + # Get outmap parameter | |
| 524 | + outfile = self['outfile'].filename() | |
| 525 | + | |
| 526 | + # Log file name | |
| 527 | + self._log_value(gammalib.NORMAL, 'Pull distribution file', outfile.url()) | |
| 528 | + | |
| 529 | + # Save pull distribution | |
| 530 | + self._fits.saveto(outfile, self['clobber'].boolean()) | |
| 531 | + | |
| 532 | + # Return | |
| 533 | + return | |
| 534 | + | |
| 535 | + def pull_distribution(self): | |
| 536 | + """ | |
| 537 | + Return pull distribution FITS file | |
| 538 | + | |
| 539 | + Returns: | |
| 540 | + FITS file containing pull distribution | |
| 541 | + """ | |
| 542 | + # Return | |
| 543 | + return self._fits | |
| 544 | + | |
| 399 | 545 | |
| 400 | 546 | # ======================== # |
| 401 | 547 | # Main routine entry point # | ... | ... |
doc/source/admin/release_history/2.0.rst
doc/source/users/reference_manual/cspull.rst
| ... | ... | @@ -34,12 +34,10 @@ be performed in parallel over as many processes as the number of CPUs available |
| 34 | 34 | machine. The maximum number of parallel processes can be set by the user through the |
| 35 | 35 | ``nthreads`` hidden parameter. |
| 36 | 36 | |
| 37 | -:ref:`cspull` will generate an ASCII file in comma-separated value (CSV) format, | |
| 38 | -containing one row per pull. The first row is a header row providing the | |
| 39 | -column names. The following rows give the pull results, one row per pull. | |
| 40 | -This includes for each parameter the parameter value, error and pull. The | |
| 41 | -maximum likelihood value and the observed and the estimated number of counts | |
| 42 | -are also given. | |
| 37 | +:ref:`cspull` will generate a FITS file containing a binary table with one row per | |
| 38 | +pull. FITS table columns include the maximum log-likelihood value of the pull, | |
| 39 | +the number of simulated and fitted events, and for each parameter the parameter | |
| 40 | +value, error and pull. | |
| 43 | 41 | |
| 44 | 42 | From the output file, pull distribution plots can be generated using for |
| 45 | 43 | example the ``show_pull_histogram.py`` script in the examples folder. The | ... | ... |
examples/show_pull_evolution.py
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | # ========================================================================== |
| 3 | 3 | # Shows the evolution of the mean and rms pull |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2015-2020 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2015-2021 Juergen Knoedlseder | |
| 6 | 6 | # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify |
| 8 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -31,9 +31,80 @@ try: |
| 31 | 31 | except ImportError: |
| 32 | 32 | print('This script needs the "numpy" module') |
| 33 | 33 | sys.exit() |
| 34 | +import gammalib | |
| 34 | 35 | import cscripts |
| 35 | 36 | |
| 36 | 37 | |
| 38 | +# ====================== # | |
| 39 | +# Read pull distribution # | |
| 40 | +# ====================== # | |
| 41 | +def read_pull_distribution(filename, parname): | |
| 42 | + """ | |
| 43 | + Read pull distribution | |
| 44 | + | |
| 45 | + Parameters | |
| 46 | + ---------- | |
| 47 | + filename : str | |
| 48 | + Name of FITS or CSV file | |
| 49 | + parname : str | |
| 50 | + Parameter name | |
| 51 | + | |
| 52 | + Returns | |
| 53 | + ------- | |
| 54 | + pulls : dict | |
| 55 | + Dictionary with pull distribution | |
| 56 | + """ | |
| 57 | + # Read pull distribution dependent on file type | |
| 58 | + fname = gammalib.GFilename(filename) | |
| 59 | + if fname.is_fits(): | |
| 60 | + pulls = read_pull_distribution_fits(filename, parname) | |
| 61 | + else: | |
| 62 | + pulls = cscripts.ioutils.read_pull_values(filename, parname) | |
| 63 | + | |
| 64 | + # Return pulls | |
| 65 | + return pulls | |
| 66 | + | |
| 67 | + | |
| 68 | +# ===================================== # | |
| 69 | +# Read pull distribution from FITS file # | |
| 70 | +# ===================================== # | |
| 71 | +def read_pull_distribution_fits(filename, parname): | |
| 72 | + """ | |
| 73 | + Read pull distribution from FITS file | |
| 74 | + | |
| 75 | + Parameters | |
| 76 | + ---------- | |
| 77 | + filename : str | |
| 78 | + Name of FITS file | |
| 79 | + parname : str | |
| 80 | + Parameter name | |
| 81 | + | |
| 82 | + Returns | |
| 83 | + ------- | |
| 84 | + pulls : dict | |
| 85 | + Dictionary with pull distribution | |
| 86 | + """ | |
| 87 | + # Open FITS file | |
| 88 | + fits = gammalib.GFits(filename) | |
| 89 | + | |
| 90 | + # Get pull distribution table | |
| 91 | + table = fits.table('PULL_DISTRIBUTION') | |
| 92 | + | |
| 93 | + # Get relevant column | |
| 94 | + column = table[parname] | |
| 95 | + | |
| 96 | + # Initialise vector | |
| 97 | + pulls = [] | |
| 98 | + | |
| 99 | + # Fill vectors | |
| 100 | + nrows = table.nrows() | |
| 101 | + for row in range(nrows): | |
| 102 | + pulls.append(column[row]) | |
| 103 | + | |
| 104 | + # Return | |
| 105 | + return pulls | |
| 106 | + | |
| 107 | + | |
| 37 | 108 | # =================== # |
| 38 | 109 | # Plot pull evolution # |
| 39 | 110 | # =================== # |
| ... | ... | @@ -50,8 +121,11 @@ def plot_pull_evolution(filename, parname, plotfile): |
| 50 | 121 | plotfile : str |
| 51 | 122 | Plot filename |
| 52 | 123 | """ |
| 53 | - # Read values from CSV file | |
| 54 | - values = np.array(cscripts.ioutils.read_pull_values(filename, parname)) | |
| 124 | + # Read pull distribution from file | |
| 125 | + pulls = read_pull_distribution(filename, parname) | |
| 126 | + | |
| 127 | + # Create Numpy array | |
| 128 | + values = np.array(pulls) | |
| 55 | 129 | |
| 56 | 130 | # Compute mean and rms |
| 57 | 131 | nsamples = len(values) |
| ... | ... | @@ -97,7 +171,7 @@ def show_pull_evolution(): |
| 97 | 171 | usage = 'show_pull_evolution.py [-p plotfile] file parameter' |
| 98 | 172 | |
| 99 | 173 | # Set default options |
| 100 | - options = [{'option': '-p', 'value': ''}] | |
| 174 | + options = [{'option': '-p', 'value': ''}] | |
| 101 | 175 | |
| 102 | 176 | # Get arguments and options from command line arguments |
| 103 | 177 | args, options = cscripts.ioutils.get_args_options(options, usage) | ... | ... |
examples/show_pull_histogram.py
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | # ========================================================================== |
| 3 | 3 | # Shows the pull histogram |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2011-2020 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2011-2021 Juergen Knoedlseder | |
| 6 | 6 | # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify |
| 8 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -31,9 +31,80 @@ try: |
| 31 | 31 | except ImportError: |
| 32 | 32 | print('This script needs the "numpy" module') |
| 33 | 33 | sys.exit() |
| 34 | +import gammalib | |
| 34 | 35 | import cscripts |
| 35 | 36 | |
| 36 | 37 | |
| 38 | +# ====================== # | |
| 39 | +# Read pull distribution # | |
| 40 | +# ====================== # | |
| 41 | +def read_pull_distribution(filename, parname): | |
| 42 | + """ | |
| 43 | + Read pull distribution | |
| 44 | + | |
| 45 | + Parameters | |
| 46 | + ---------- | |
| 47 | + filename : str | |
| 48 | + Name of FITS or CSV file | |
| 49 | + parname : str | |
| 50 | + Parameter name | |
| 51 | + | |
| 52 | + Returns | |
| 53 | + ------- | |
| 54 | + pulls : dict | |
| 55 | + Dictionary with pull distribution | |
| 56 | + """ | |
| 57 | + # Read pull distribution dependent on file type | |
| 58 | + fname = gammalib.GFilename(filename) | |
| 59 | + if fname.is_fits(): | |
| 60 | + pulls = read_pull_distribution_fits(filename, parname) | |
| 61 | + else: | |
| 62 | + pulls = cscripts.ioutils.read_pull_values(filename, parname) | |
| 63 | + | |
| 64 | + # Return pulls | |
| 65 | + return pulls | |
| 66 | + | |
| 67 | + | |
| 68 | +# ===================================== # | |
| 69 | +# Read pull distribution from FITS file # | |
| 70 | +# ===================================== # | |
| 71 | +def read_pull_distribution_fits(filename, parname): | |
| 72 | + """ | |
| 73 | + Read pull distribution from FITS file | |
| 74 | + | |
| 75 | + Parameters | |
| 76 | + ---------- | |
| 77 | + filename : str | |
| 78 | + Name of FITS file | |
| 79 | + parname : str | |
| 80 | + Parameter name | |
| 81 | + | |
| 82 | + Returns | |
| 83 | + ------- | |
| 84 | + pulls : dict | |
| 85 | + Dictionary with pull distribution | |
| 86 | + """ | |
| 87 | + # Open FITS file | |
| 88 | + fits = gammalib.GFits(filename) | |
| 89 | + | |
| 90 | + # Get pull distribution table | |
| 91 | + table = fits.table('PULL_DISTRIBUTION') | |
| 92 | + | |
| 93 | + # Get relevant column | |
| 94 | + column = table[parname] | |
| 95 | + | |
| 96 | + # Initialise vector | |
| 97 | + pulls = [] | |
| 98 | + | |
| 99 | + # Fill vectors | |
| 100 | + nrows = table.nrows() | |
| 101 | + for row in range(nrows): | |
| 102 | + pulls.append(column[row]) | |
| 103 | + | |
| 104 | + # Return | |
| 105 | + return pulls | |
| 106 | + | |
| 107 | + | |
| 37 | 108 | # =================== # |
| 38 | 109 | # Plot pull histogram # |
| 39 | 110 | # =================== # |
| ... | ... | @@ -52,8 +123,11 @@ def plot_pull_histogram(filename, parname, nbins, plotfile): |
| 52 | 123 | plotfile : str |
| 53 | 124 | Plot filename |
| 54 | 125 | """ |
| 55 | - # Read values from CSV file | |
| 56 | - values = np.array(cscripts.ioutils.read_pull_values(filename, parname)) | |
| 126 | + # Read pull distribution from file | |
| 127 | + pulls = read_pull_distribution(filename, parname) | |
| 128 | + | |
| 129 | + # Create Numpy array | |
| 130 | + values = np.array(pulls) | |
| 57 | 131 | |
| 58 | 132 | # Create histogram |
| 59 | 133 | _, bins, _ = plt.hist(values, nbins, range=[-4.0,4.0], | ... | ... |
test/data/pull.fits
0 → 100644
No preview for this file type
test/test_cspull.py
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | # ========================================================================== |
| 3 | 3 | # This scripts performs unit tests for the cspull script. |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2016-2018 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2016-2021 Juergen Knoedlseder | |
| 6 | 6 | # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify |
| 8 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -77,7 +77,7 @@ class Test(test): |
| 77 | 77 | |
| 78 | 78 | # Setup cspull command |
| 79 | 79 | cmd = cspull+' inmodel="'+self._model+'" onsrc="NONE"'+ \ |
| 80 | - ' outfile="cspull_cmd1.dat"'+ \ | |
| 80 | + ' outfile="cspull_cmd1.fits"'+ \ | |
| 81 | 81 | ' ntrials=2'+ \ |
| 82 | 82 | ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ |
| 83 | 83 | ' ra=83.6331 dec=22.0145 emin=1.0 emax=100.0'+ \ |
| ... | ... | @@ -90,11 +90,11 @@ class Test(test): |
| 90 | 90 | 'Check successful execution from command line') |
| 91 | 91 | |
| 92 | 92 | # Check pull distribution file |
| 93 | - self._check_pull_file('cspull_cmd1.dat') | |
| 93 | + self._check_pull_file('cspull_cmd1.fits') | |
| 94 | 94 | |
| 95 | 95 | # Setup cspull command |
| 96 | 96 | cmd = cspull+' inmodel="model_that_does_not_exist.xml" onsrc="NONE"'+ \ |
| 97 | - ' outfile="cspull_cmd1.dat"'+ \ | |
| 97 | + ' outfile="cspull_cmd1.fits"'+ \ | |
| 98 | 98 | ' ntrials=2'+ \ |
| 99 | 99 | ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ |
| 100 | 100 | ' ra=83.6331 dec=22.0145 emin=1.0 emax=100.0'+ \ |
| ... | ... | @@ -121,7 +121,7 @@ class Test(test): |
| 121 | 121 | pull = cscripts.cspull() |
| 122 | 122 | pull['inmodel'] = self._model |
| 123 | 123 | pull['onsrc'] = 'NONE' |
| 124 | - pull['outfile'] = 'cspull_py1.dat' | |
| 124 | + pull['outfile'] = 'cspull_py1.fits' | |
| 125 | 125 | pull['ntrials'] = 2 |
| 126 | 126 | pull['caldb'] = self._caldb |
| 127 | 127 | pull['irf'] = self._irf |
| ... | ... | @@ -140,15 +140,16 @@ class Test(test): |
| 140 | 140 | # Run cspull script |
| 141 | 141 | pull.logFileOpen() # Make sure we get a log file |
| 142 | 142 | pull.run() |
| 143 | + pull.save() | |
| 143 | 144 | |
| 144 | 145 | # Check pull distribution file |
| 145 | - self._check_pull_file('cspull_py1.dat') | |
| 146 | + self._check_pull_file('cspull_py1.fits') | |
| 146 | 147 | |
| 147 | 148 | # Set-up binned cspull |
| 148 | 149 | pull = cscripts.cspull() |
| 149 | 150 | pull['inmodel'] = self._model |
| 150 | 151 | pull['onsrc'] = 'NONE' |
| 151 | - pull['outfile'] = 'cspull_py2.dat' | |
| 152 | + pull['outfile'] = 'cspull_py2.fits' | |
| 152 | 153 | pull['ntrials'] = 1 |
| 153 | 154 | pull['caldb'] = self._caldb |
| 154 | 155 | pull['irf'] = self._irf |
| ... | ... | @@ -172,14 +173,14 @@ class Test(test): |
| 172 | 173 | pull.execute() |
| 173 | 174 | |
| 174 | 175 | # Check pull distribution file |
| 175 | - self._check_pull_file('cspull_py2.dat', rows=2) | |
| 176 | + self._check_pull_file('cspull_py2.fits', rows=1) | |
| 176 | 177 | |
| 177 | 178 | # Set-up cspull from event list |
| 178 | 179 | pull = cscripts.cspull() |
| 179 | 180 | pull['inobs'] = self._events |
| 180 | 181 | pull['onsrc'] = 'NONE' |
| 181 | 182 | pull['inmodel'] = self._model |
| 182 | - pull['outfile'] = 'cspull_py3.dat' | |
| 183 | + pull['outfile'] = 'cspull_py3.fits' | |
| 183 | 184 | pull['ntrials'] = 1 |
| 184 | 185 | pull['caldb'] = self._caldb |
| 185 | 186 | pull['irf'] = self._irf |
| ... | ... | @@ -191,7 +192,7 @@ class Test(test): |
| 191 | 192 | pull.execute() |
| 192 | 193 | |
| 193 | 194 | # Check pull distribution file |
| 194 | - self._check_pull_file('cspull_py3.dat', rows=2) | |
| 195 | + self._check_pull_file('cspull_py3.fits', rows=1) | |
| 195 | 196 | |
| 196 | 197 | # Build observation container with unbinned observation |
| 197 | 198 | cta = gammalib.GCTAObservation(self._events) |
| ... | ... | @@ -202,7 +203,7 @@ class Test(test): |
| 202 | 203 | pull = cscripts.cspull(obs) |
| 203 | 204 | pull['inmodel'] = self._model |
| 204 | 205 | pull['onsrc'] = 'NONE' |
| 205 | - pull['outfile'] = 'cspull_py4.dat' | |
| 206 | + pull['outfile'] = 'cspull_py4.fits' | |
| 206 | 207 | pull['ntrials'] = 1 |
| 207 | 208 | pull['caldb'] = self._caldb |
| 208 | 209 | pull['irf'] = self._irf |
| ... | ... | @@ -214,7 +215,7 @@ class Test(test): |
| 214 | 215 | pull.execute() |
| 215 | 216 | |
| 216 | 217 | # Check pull distribution file |
| 217 | - self._check_pull_file('cspull_py4.dat', rows=2) | |
| 218 | + self._check_pull_file('cspull_py4.fits', rows=1) | |
| 218 | 219 | |
| 219 | 220 | # Set-up stacked cspull with one observation where response cubes |
| 220 | 221 | # are specified in the input observation |
| ... | ... | @@ -222,7 +223,7 @@ class Test(test): |
| 222 | 223 | pull['inobs'] = self._inobs |
| 223 | 224 | pull['inmodel'] = self._stacked_model |
| 224 | 225 | pull['onsrc'] = 'NONE' |
| 225 | - pull['outfile'] = 'cspull_py5.dat' | |
| 226 | + pull['outfile'] = 'cspull_py5.fits' | |
| 226 | 227 | pull['ntrials'] = 1 |
| 227 | 228 | pull['enumbins'] = 0 |
| 228 | 229 | pull['logfile'] = 'cspull_py5.log' |
| ... | ... | @@ -232,7 +233,7 @@ class Test(test): |
| 232 | 233 | pull.execute() |
| 233 | 234 | |
| 234 | 235 | # Check pull distribution file |
| 235 | - self._check_pull_file('cspull_py5.dat', rows=2) | |
| 236 | + self._check_pull_file('cspull_py5.fits', rows=1) | |
| 236 | 237 | |
| 237 | 238 | # Set-up stacked cspull with two observations from which response |
| 238 | 239 | # cubes will be computed internally. The IRFs are also specified |
| ... | ... | @@ -242,7 +243,7 @@ class Test(test): |
| 242 | 243 | pull['inobs'] = self._inobs_two |
| 243 | 244 | pull['inmodel'] = self._model |
| 244 | 245 | pull['onsrc'] = 'NONE' |
| 245 | - pull['outfile'] = 'cspull_py6.dat' | |
| 246 | + pull['outfile'] = 'cspull_py6.fits' | |
| 246 | 247 | pull['ntrials'] = 1 |
| 247 | 248 | pull['emin'] = 1.0 |
| 248 | 249 | pull['emax'] = 100.0 |
| ... | ... | @@ -258,13 +259,13 @@ class Test(test): |
| 258 | 259 | pull.execute() |
| 259 | 260 | |
| 260 | 261 | # Check pull distribution file |
| 261 | - self._check_pull_file('cspull_py6.dat', rows=2) | |
| 262 | + self._check_pull_file('cspull_py6.fits', rows=1) | |
| 262 | 263 | |
| 263 | 264 | # Set-up cspull without multiprocessing |
| 264 | 265 | pull = cscripts.cspull() |
| 265 | 266 | pull['inmodel'] = self._model |
| 266 | 267 | pull['onsrc'] = 'NONE' |
| 267 | - pull['outfile'] = 'cspull_py7.dat' | |
| 268 | + pull['outfile'] = 'cspull_py7.fits' | |
| 268 | 269 | pull['ntrials'] = 2 |
| 269 | 270 | pull['caldb'] = self._caldb |
| 270 | 271 | pull['irf'] = self._irf |
| ... | ... | @@ -284,9 +285,10 @@ class Test(test): |
| 284 | 285 | # Run cspull script |
| 285 | 286 | pull.logFileOpen() # Make sure we get a log file |
| 286 | 287 | pull.run() |
| 288 | + pull.save() | |
| 287 | 289 | |
| 288 | 290 | # Check pull distribution file |
| 289 | - self._check_pull_file('cspull_py7.dat') | |
| 291 | + self._check_pull_file('cspull_py7.fits') | |
| 290 | 292 | |
| 291 | 293 | # Return |
| 292 | 294 | return |
| ... | ... | @@ -316,7 +318,7 @@ class Test(test): |
| 316 | 318 | pull['deadc'] = 0.98 |
| 317 | 319 | pull['rad'] = 5.0 |
| 318 | 320 | pull['logfile'] = 'cspull_py1_pickle.log' |
| 319 | - pull['outfile'] = 'cspull_py1_pickle.dat' | |
| 321 | + pull['outfile'] = 'cspull_py1_pickle.fits' | |
| 320 | 322 | pull['chatter'] = 2 |
| 321 | 323 | |
| 322 | 324 | # Perform pickeling tests of filled class |
| ... | ... | @@ -325,23 +327,25 @@ class Test(test): |
| 325 | 327 | # Run csphasecrv script and save light curve |
| 326 | 328 | obj.logFileOpen() # Make sure we get a log file |
| 327 | 329 | obj.run() |
| 330 | + obj.save() | |
| 328 | 331 | |
| 329 | 332 | # Check pull distribution file |
| 330 | - self._check_pull_file('cspull_py1_pickle.dat') | |
| 333 | + self._check_pull_file('cspull_py1_pickle.fits') | |
| 331 | 334 | |
| 332 | 335 | # Return |
| 333 | 336 | return |
| 334 | 337 | |
| 335 | 338 | # Check pull file |
| 336 | - def _check_pull_file(self, filename, rows=3): | |
| 339 | + def _check_pull_file(self, filename, rows=2): | |
| 337 | 340 | """ |
| 338 | 341 | Check pull file |
| 339 | 342 | """ |
| 340 | - # Open pull file as CSV file | |
| 341 | - pulls = gammalib.GCsv(filename) | |
| 343 | + # Open pull file as FITS file | |
| 344 | + fits = gammalib.GFits(filename) | |
| 345 | + table = fits.table('PULL_DISTRIBUTION') | |
| 342 | 346 | |
| 343 | 347 | # Check dimensions |
| 344 | - self.test_value(pulls.nrows(), rows, | |
| 348 | + self.test_value(table.nrows(), rows, | |
| 345 | 349 | 'Check for number of rows in pull file') |
| 346 | 350 | |
| 347 | 351 | # Return | ... | ... |
test/test_examples.py.in
| ... | ... | @@ -295,12 +295,19 @@ class Test(gammalib.GPythonTestSuite): |
| 295 | 295 | Test show_pull_evolution |
| 296 | 296 | """ |
| 297 | 297 | # Set arguments |
| 298 | - args = ('-p example_show_pull_evolution.eps %s/pull.dat ' | |
| 298 | + args = ('-p example_show_pull_evolution_csv.eps %s/pull.dat ' | |
| 299 | 299 | 'Pull_Crab_Prefactor' % self._datadir) |
| 300 | 300 | |
| 301 | 301 | # Execute script |
| 302 | 302 | self._execute_python('show_pull_evolution', args=args) |
| 303 | 303 | |
| 304 | + # Set arguments | |
| 305 | + args = ('-p example_show_pull_evolution_fits.eps %s/pull.fits ' | |
| 306 | + 'CRAB_PREFACTOR_PULL' % self._datadir) | |
| 307 | + | |
| 308 | + # Execute script | |
| 309 | + self._execute_python('show_pull_evolution', args=args) | |
| 310 | + | |
| 304 | 311 | # Return |
| 305 | 312 | return |
| 306 | 313 | |
| ... | ... | @@ -310,19 +317,26 @@ class Test(gammalib.GPythonTestSuite): |
| 310 | 317 | Test show_pull_histogram |
| 311 | 318 | """ |
| 312 | 319 | # Set arguments |
| 313 | - args = ('-p example_show_pull_histogram1.eps %s/pull.dat ' | |
| 320 | + args = ('-p example_show_pull_histogram_csv1.eps %s/pull.dat ' | |
| 314 | 321 | 'Pull_Crab_Prefactor' % self._datadir) |
| 315 | 322 | |
| 316 | 323 | # Execute script |
| 317 | 324 | self._execute_python('show_pull_histogram', args=args) |
| 318 | 325 | |
| 319 | 326 | # Set arguments |
| 320 | - args = ('-n 20 -p example_show_pull_histogram2.eps %s/pull.dat ' | |
| 327 | + args = ('-n 20 -p example_show_pull_histogram_csv2.eps %s/pull.dat ' | |
| 321 | 328 | 'Pull_Crab_Prefactor' % self._datadir) |
| 322 | 329 | |
| 323 | 330 | # Execute script |
| 324 | 331 | self._execute_python('show_pull_histogram', args=args) |
| 325 | 332 | |
| 333 | + # Set arguments | |
| 334 | + args = ('-p example_show_pull_histogram_fits.eps %s/pull.fits ' | |
| 335 | + 'CRAB_PREFACTOR_PULL' % self._datadir) | |
| 336 | + | |
| 337 | + # Execute script | |
| 338 | + self._execute_python('show_pull_histogram', args=args) | |
| 339 | + | |
| 326 | 340 | # Return |
| 327 | 341 | return |
| 328 | 342 | ... | ... |