diff --git a/ChangeLog b/ChangeLog index bddbb13..a161fc6 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-06-21 Juergen Knoedlseder +2016-06-25 Juergen Knoedlseder * Version 1.1.0 released ======================== + Add unit tests for all example scripts Handle zero parameter errors in ctulimit (#1795) Compute counts cube weights in ctbin (#1789) Adapt to new stacked cube interface (#1788) diff --git a/NEWS b/NEWS index 23da3f5..db7bb1f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.1.0 -21 June 2016 +25 June 2016 Introduction diff --git a/README.md b/README.md index 41cc4f9..dee57cc 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ctools information ================== -* Version: 1.1.0.dev1 (21 June 2016) +* Version: 1.1.0.dev1 (25 June 2016) * Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) * GammaLib dependency: 1.1.0.dev1 diff --git a/examples/README.md b/examples/README.md index 291918c..5aabb5b 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1,70 +1,49 @@ ctools examples =============== -This folder contains Python example scripts that illustrate the -usage of the ctools. - -`benchmark_cta_analysis.py` -> This script peforms a benchmark of the various CTA analysis types -> (unbinned, binned, stacked), and if matplotlib is installed, creates -> a plot of the benchmark results. - -`check_models.py` -> Checks a variety of models using an unbinned analysis. +This folder contains Python example scripts that illustrate the usage of +the ctools. `generate_prod3_irfs/py` -> Generate Prod3 Instrument Response Functions from FITS files. +> Generates Prod3 Instrument Response Functions from FITS files. `make_pointings.py` -> This script creates pointings for CTA surveys. For the moment the -> Galactic Plane Survey (gps), the Extragalactic Survey (extgal), -> the Galactic Centre Survey (gc) and the LMC survey (lmc) are -> implemented. The script outputs and ASCII file that can be used by -> csobsdef. - -`make_pull_at_sensitivity_limit.py` -> This script creates pull distributions as function of energy -> for a test source that is at the 5 sigma sensitivity limit for -> 50h of effective observing time. It takes as input file the -> output file generated using cspull. -> If the Python processing module is installed the script can make -> use of multiple cores and runs in parallel. +> Creates pointings for CTA surveys. For the moment the Galactic Plane +> Survey (gps), the Extragalactic Survey (extgal), the Galactic Centre +> Survey (gc) and the LMC survey (lmc) are implemented. The script outputs +> and ASCII file that can be used by csobsdef. `make_spectrum.py` -> This script shows how to make a spectrum using obsutils. - -`make_survey.py` -> This script simulates a CTA galactic plane survey. +> Simulates events and generates a spectrum from the using csspec. `make_ts_distributions.py` -> This script creates TS distributions as function of energy for -> an off observation of 50h effective observing time. -> If the Python processing module is installed the script can make -> use of multiple cores and runs in parallel. +> Creates TS distributions as function of energy for an OFF observaition +> of 50 hours effective observing time. If the Python processing module is +> installed the script can make use of multiple cores and runs in parallel. `pipeline_binned_disk.py` -> Example pipeline for simulating events and performing a binned -> analysis by saving intermediate results to disk. +> Example pipeline for simulating events and performing a binned analysis +> by saving intermediate results to disk. `pipeline_binned_mem.py` -> Example pipeline for simulating events and performing a binned -> analysis in memory. No files will be written to disk. +> Example pipeline for simulating events and performing a binned analysis +> in memory. No files will be written to disk. `pipeline_stacked_disk.py` -> Example pipeline for simulating events and performing a stacked -> analysis by saving intermediate results to disk. +> Example pipeline for simulating events and performing a stacked analysis +> by saving intermediate results to disk. `pipeline_stacked_mem.py` -> Example pipeline for simulating events and performing a stacked -> analysis in memory. No files will be written to disk. +> Example pipeline for simulating events and performing a stacked analysis +> in memory. No files will be written to disk. `pipeline_unbinned_disk.py` -> Example pipeline for simulating events and performing an unbinned -> analysis by saving intermediate results to disk. +> Example pipeline for simulating events and performing an unbinned analysis +> by saving intermediate results to disk. `pipeline_unbinned_mem.py` -> Example pipeline for simulating events and performing an unbinned -> analysis in memory. No files will be written to disk. +> Example pipeline for simulating events and performing an unbinned analysis +> in memory. No files will be written to disk. `show_butterfly.py` > Displays a butterfly generated using the ctbutterfly tool. @@ -86,14 +65,10 @@ usage of the ctools. > Displays one or several CTA response functions. `show_sensitivity.py` -> Displays sensitivity information generated by the cssens script. +> Displays the sensitivity generated by the cssens script. `show_spectrum.py` > Displays spectrum generated by the csspec script. `show_ts_distribution.py` -> Displays Test Statistic histograms generated using the cstsdist -> script. - -`simulate_events.py` -> Simulates events and compares the simulation to the model. +> Displays Test Statistic histograms generated using the cstsdist script. diff --git a/examples/make_ts_distributions.py b/examples/make_ts_distributions.py index c362279..a1412ca 100755 --- a/examples/make_ts_distributions.py +++ b/examples/make_ts_distributions.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# Generate TS distribution +# Generate TS distributions as function of energy for an OFF observation # # Copyright (C) 2011-2016 Jurgen Knodlseder # @@ -65,7 +65,7 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, tsdist['outfile'] = outfile tsdist['ntrials'] = ntrials tsdist['caldb'] = 'prod2' - tsdist['irf'] = 'South_0.5h' + tsdist['irf'] = 'South_50h' tsdist['ra'] = 83.63 tsdist['dec'] = 22.01 tsdist['emin'] = emin @@ -76,7 +76,6 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, tsdist['rad'] = 5.0 tsdist['npix'] = 200 tsdist['binsz'] = 0.05 - #tsdist['debug'] = True # Optionally open the log file if log: @@ -95,7 +94,7 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, if __name__ == '__main__': # Get input arguments - usage = 'make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]' + usage = 'make_ts_distributions.py [-n ntrials] [-e enumbins] [-m max_threads]' if len(sys.argv) < 1: print(usage) sys.exit() diff --git a/examples/show_butterfly.py b/examples/show_butterfly.py index 1704df5..a3539e7 100755 --- a/examples/show_butterfly.py +++ b/examples/show_butterfly.py @@ -32,7 +32,7 @@ except: # ============ # # Script entry # # ============ # -if __name__ == "__main__": +if __name__ == '__main__': # Check for given butterfly file if len(sys.argv) < 2: @@ -83,7 +83,7 @@ if __name__ == "__main__": plt.plot(line_x,line_y,color='black',ls='-') plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5) plt.xlabel('Energy (MeV)') - plt.ylabel('E dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)') + plt.ylabel(r'E $\times$ dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)') # Show spectrum or save it into file if len(plotfile) > 0: diff --git a/examples/show_model.py b/examples/show_model.py index d8e0792..a160016 100755 --- a/examples/show_model.py +++ b/examples/show_model.py @@ -103,7 +103,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100, plotfile=''): # ============= # # Script entry # # ============= # -if __name__ == "__main__": +if __name__ == '__main__': # Check for model filename if len(sys.argv) < 2: diff --git a/examples/show_pha.py b/examples/show_pha.py index d760e5e..904ea9c 100755 --- a/examples/show_pha.py +++ b/examples/show_pha.py @@ -95,7 +95,7 @@ def show_pha(pha, plotfile=''): if __name__ == '__main__': # Print usage information - usage = 'Usage: show_pha filename [file]' + usage = 'Usage: show_pha.py filename [file]' if len(sys.argv) < 2: print(usage) sys.exit() diff --git a/examples/show_pull_evolution.py b/examples/show_pull_evolution.py index 5066e9d..06bace2 100755 --- a/examples/show_pull_evolution.py +++ b/examples/show_pull_evolution.py @@ -89,7 +89,7 @@ def read_pull(filename, parname): if __name__ == '__main__': # Print usage information - usage = 'Usage: show_pull_evolution filename parname [file]' + usage = 'Usage: show_pull_evolution.py filename parname [file]' if len(sys.argv) < 3: print(usage) sys.exit() diff --git a/examples/show_pull_histogram.py b/examples/show_pull_histogram.py index ccd2311..7a8521e 100755 --- a/examples/show_pull_histogram.py +++ b/examples/show_pull_histogram.py @@ -90,7 +90,7 @@ def read_pull(filename, parname): if __name__ == '__main__': # Print usage information - usage = 'Usage: show_pull_histogram filename parname [bins] [file]' + usage = 'Usage: show_pull_histogram.py filename parname [bins] [file]' if len(sys.argv) < 3: print(usage) sys.exit() diff --git a/examples/show_response.py b/examples/show_response.py index 1e8900d..468f736 100755 --- a/examples/show_response.py +++ b/examples/show_response.py @@ -1,12 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays one or several CTA response functions. +# Shows one or several CTA response functions. # -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2014-2015 Juergen Knoedlseder +# Copyright (C) 2014-2016 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -25,7 +21,13 @@ import sys import math import gammalib -import matplotlib.pyplot as plt +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() # ================================== # @@ -35,11 +37,19 @@ def sigma_lima(Non, Noff, alpha=0.2): """ Compute Eq. (17) of Li & Ma (1983). - Parameters: - Non - Number of on counts - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure + Parameters + ---------- + Non : float + Number of On counts + Noff : float + Number of Off counts + alpha : float, optional + Ratio of on-to-off exposure + + Returns + ------- + sigma : float + Detection significance in Gaussian sigma """ # Compute sensitivity alpha1 = alpha + 1.0 @@ -61,11 +71,19 @@ def Non_lima(sigma, Noff, alpha=0.2): """ Solve Eq. (17) of Li & Ma (1983) for Non - Parameters: - sigma - Required significance - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure + Parameters + ---------- + sigma : float + Required detection significance in Gaussian sigma + Noff : float + Number of Off counts + alpha : float, optional + Ratio of on-to-off exposure + + Returns + ------- + Non : float + Number of On counts """ # Set initial guess for Non S = 25.0 @@ -89,9 +107,18 @@ def Non_lima(sigma, Noff, alpha=0.2): # ======================= # # Show one effective area # # ======================= # -def show_one_effective_area(rsp, name, color="r"): +def show_one_effective_area(rsp, name, color='r'): """ - Show one effective area. + Show one effective area + + Parameters + ---------- + rsp : `~gammalib.GCTAResponse` + Response function + name : str + Name of the response function + color : str, optional + Color for plot """ # Set to figure 1 plt.figure(1) @@ -111,11 +138,11 @@ def show_one_effective_area(rsp, name, color="r"): plt.loglog(E, aeff, color+'o') # Set axes - plt.xlabel("Energy [TeV]") - plt.ylabel("Effective area [cm2]") + plt.xlabel('Energy (TeV)') + plt.ylabel(r'Effective area (cm$^{2}$)') # Set legend - plt.legend(loc="lower right") + plt.legend(loc='lower right') # Return return @@ -124,9 +151,18 @@ def show_one_effective_area(rsp, name, color="r"): # ======================== # # Show one background rate # # ======================== # -def show_one_background_rate(rsp, name, color="r"): +def show_one_background_rate(rsp, name, color='r'): """ - Show one background rate. + Show one background rate + + Parameters + ---------- + rsp : `~gammalib.GCTAResponse` + Response function + name : str + Name of the response function + color : str, optional + Color for plot """ # Set to figure 2 plt.figure(2) @@ -148,11 +184,11 @@ def show_one_background_rate(rsp, name, color="r"): plt.loglog(E, bgrate, color+'o') # Set axes - plt.xlabel("Energy [TeV]") - plt.ylabel("Background rate [events/s/MeV/sr]") + plt.xlabel('Energy (TeV)') + plt.ylabel(r'Background rate (counts s$^{-1}$ MeV$^{-1}$ sr$^{-1}$)') # Set legend - plt.legend(loc="lower left") + plt.legend(loc='lower left') # Return return @@ -161,9 +197,25 @@ def show_one_background_rate(rsp, name, color="r"): # ==================== # # Show one sensitivity # # ==================== # -def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sigma=5.0): +def show_one_sensitivity(rsp, name, color='r', duration=180000.0, alpha=0.2, + sigma=5.0): """ - Show one sensitivity. + Show one sensitivity + + Parameters + ---------- + rsp : `~gammalib.GCTAResponse` + Response function + name : str + Name of the response function + color : str, optional + Color for plot + duration : float, optional + Duration of observation (s) + alpha : float, optional + Ratio of on-to-off exposure + sigma : float, optional + Required detection significance in Gaussian sigma """ # Set constants r68_to_sigma = 0.6624305 @@ -194,9 +246,9 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig src_counts = 0.05*bgd_counts if src_counts < 10: src_counts = 10.0 - emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), "TeV") - emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), "TeV") - epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), "TeV") + emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), 'TeV') + emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), 'TeV') + epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), 'TeV') plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot) conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax) flux[i] = conv * src_counts / (duration * aeff*0.68) @@ -206,11 +258,11 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig plt.loglog(E, flux, color+'o') # Set axes - plt.xlabel("Energy [TeV]") - plt.ylabel("Sensitivity [erg/cm2/s]") + plt.xlabel('Energy (TeV)') + plt.ylabel(r'Sensitivity (erg cm$^{-2}$ s$^{-1}$)') # Set legend - plt.legend(loc="upper right") + plt.legend(loc='upper right') # Return return @@ -219,9 +271,22 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig # ==================== # # Show one sensitivity # # ==================== # -def show_one_response(rspname, dbname, name, rootdir=None, color="r"): +def show_one_response(rspname, dbname, name, rootdir=None, color='r'): """ Show one response. + + Parameters + ---------- + rspname : str + Response name + dbname : str + Database name + name : str + Name of the response function + rootdir : str, optional + Response root directory + color : str, optional + Color for plot """ # Set-up calibration database caldb = gammalib.GCaldb() @@ -230,7 +295,7 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): if gammalib.dir_exists(dbname): caldb.rootdir(dbname) else: - caldb.open("cta", dbname) + caldb.open('cta', dbname) # Load response function rsp = gammalib.GCTAResponseIrf(rspname, caldb) @@ -253,28 +318,35 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): # ======================== # if __name__ == '__main__': + # Print usage information + usage = 'Usage: show_response.py [file]' + if len(sys.argv) < 1: + print(usage) + sys.exit() + + # Check if plotting in file is requested + plotfile = '' + if len(sys.argv) == 2: + plotfile = sys.argv[1] + # Create figures plt.figure(1) - plt.title("Effective area") + plt.title('Effective area') plt.figure(2) - plt.title("Background rate") + plt.title('Background rate') plt.figure(3) - plt.title("Sensitivity") + plt.title('Sensitivity') # Set response dictionary rsps = [ - {'dbname': "dummy", - 'rspname': "cta_dummy_irf", - 'name': "Prod1 (Array E, MPIK)", - 'color': "b"}, - {'dbname': "e", - 'rspname': "IFAE20120510_50h", - 'name': "Prod1 (Array E, IFAE)", - 'color': "g"}, - {'dbname': "aar", - 'rspname': "DESY20140105_50h", - 'name': "Prod2 (Aar, DESY)", - 'color': "r"} + {'dbname': 'irf', + 'rspname': 'cta_dummy_irf.dat', + 'name': 'Prod1 (Array E, MPIK)', + 'color': 'b'}, + {'dbname': 'irf', + 'rspname': 'irf_file.fits', + 'name': 'Prod2 (Aar, IFAE)', + 'color': 'r'} ] # Loop over all responses @@ -293,5 +365,8 @@ if __name__ == '__main__': # Show response show_one_response(rspname, dbname, name, rootdir=rootdir, color=color) - # Show plot - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() diff --git a/examples/show_sensitivity.py b/examples/show_sensitivity.py index bedcf27..f72dcd3 100755 --- a/examples/show_sensitivity.py +++ b/examples/show_sensitivity.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays a differential sensitivity plot generated by cssens. +# Display differential sensitivity plot generated by cssens # -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -23,8 +23,11 @@ import math import gammalib try: import matplotlib.pyplot as plt + plt.figure() + plt.close() except: - sys.exit("This script needs matplotlib") + print('This script needs the "matplotlib" module') + sys.exit() # ===================== # @@ -32,7 +35,17 @@ except: # ===================== # def read_sensitivity(filename): """ - Read sensitivity information from CSV file. + Read sensitivity information from CSV file + + Parameters + ---------- + filename : str + Name of CSV file + + Returns + ------- + sensitivity : dict + Dictionary with sensitivity information """ # Read filename csv = gammalib.GCsv(filename,',') @@ -40,39 +53,39 @@ def read_sensitivity(filename): # Create dictionnary sensitivity = {} for column in range(csv.ncols()): - name = csv[0,column] + name = csv[0,column].rstrip('\r') values = [] for row in range(csv.nrows()-1): values.append(float(csv[row+1,column])) sensitivity[name] = values # Check where file contains differential or integral sensitivity - mode = "Integral" - if sensitivity.has_key("emax"): + mode = 'Integral' + if sensitivity.has_key('emax'): emax_ref = -1.0 - for value in sensitivity["emax"]: + for value in sensitivity['emax']: if emax_ref < 0.0: emax_ref = value elif emax_ref != value: - mode = "Differential" + mode = 'Differential' break # Add mode - sensitivity["mode"] = mode + sensitivity['mode'] = mode # Add linear energy values - if mode == "Differential": - if sensitivity.has_key("loge"): + if mode == 'Differential': + if sensitivity.has_key('loge'): name = 'energy' values = [] - for value in sensitivity["loge"]: + for value in sensitivity['loge']: values.append(math.pow(10.0, value)) sensitivity[name] = values else: - if sensitivity.has_key("emin"): + if sensitivity.has_key('emin'): name = 'energy' values = [] - for value in sensitivity["emin"]: + for value in sensitivity['emin']: values.append(value) sensitivity[name] = values @@ -83,9 +96,18 @@ def read_sensitivity(filename): # ===================== # # Show sensitivity data # # ===================== # -def show_sensitivity(sensitivity, filename): +def show_sensitivity(sensitivity, filename, plotfile=''): """ - Plot sensitivity. + Plot sensitivity data + + Parameters + ---------- + sensitivity : dict + Dictionary with sensitivity information + filename : str + Name of CSV file + plotfile : str, optional + Plot file name """ # Build title title = sensitivity['mode'] + ' sensitivity ('+filename+')' @@ -102,14 +124,17 @@ def show_sensitivity(sensitivity, filename): plt.plot(sensitivity['energy'], sensitivity['sensitivity'], 'ro-') # Set labels - plt.xlabel("Energy (TeV)") - if sensitivity['mode'] == "Differential": - plt.ylabel(r"E$\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)") + plt.xlabel('Energy (TeV)') + if sensitivity['mode'] == 'Differential': + plt.ylabel(r'E $\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)') else: - plt.ylabel(r"E$\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)") + plt.ylabel(r'E $\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)') - # Show plot - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() # Return return @@ -121,9 +146,9 @@ def show_sensitivity(sensitivity, filename): if __name__ == '__main__': # Print usage information - usage = "Usage: show_sensitivity [filename]" + usage = 'Usage: show_sensitivity.py [filename] [file]' if len(sys.argv) < 1: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Extract parameters @@ -131,9 +156,12 @@ if __name__ == '__main__': filename = 'sensitivity.dat' else: filename = sys.argv[1] + plotfile = '' + if len(sys.argv) == 3: + plotfile = sys.argv[2] # Read sensitivity data sensitivity = read_sensitivity(filename) # Show sensitivity data - show_sensitivity(sensitivity, filename) + show_sensitivity(sensitivity, filename, plotfile=plotfile) diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py index 2455946..970bdc1 100755 --- a/examples/show_spectrum.py +++ b/examples/show_spectrum.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This script shows a spectrum created with csspec +# Display spectrum generated by csspec # # Copyright (C) 2015-2016 Juergen Knoedlseder # @@ -18,35 +18,41 @@ # along with this program. If not, see . # # ========================================================================== - -import gammalib import sys +import gammalib try: import matplotlib.pyplot as plt + plt.figure() + plt.close() except: - sys.exit("This script needs matplotlib") + print('This script needs the "matplotlib" module') + sys.exit() -# ============ # -# Script entry # -# ============ # -if __name__ == "__main__": - - # Check for given spectrum file - if not len(sys.argv) == 2: - sys.exit("Usage: show_spectrum.py spectrum.fits") +# ============= # +# Show spectrum # +# ============= # +def show_spectrum(filename, plotfile=''): + """ + Plot sensitivity data + Parameters + ---------- + filename : str + Name of spectrum FITS file + plotfile : str, optional + Plot file name + """ # Read spectrum file - filename = sys.argv[1] fits = gammalib.GFits(filename) table = fits.table(1) - c_energy = table["Energy"] - c_ed = table["ed_Energy"] - c_eu = table["eu_Energy"] - c_flux = table["Flux"] - c_eflux = table["e_Flux"] - c_ts = table["TS"] - c_upper = table["UpperLimit"] + c_energy = table['Energy'] + c_ed = table['ed_Energy'] + c_eu = table['eu_Energy'] + c_flux = table['Flux'] + c_eflux = table['e_Flux'] + c_ts = table['TS'] + c_upper = table['UpperLimit'] # Initialise arrays to be filled energies = [] @@ -63,25 +69,22 @@ if __name__ == "__main__": nrows = table.nrows() for row in range(nrows): - # Get TS + # Get Test Statistic, flux and flux error ts = c_ts.real(row) flx = c_flux.real(row) e_flx = c_eflux.real(row) - # Switch + # If Test Statistic is larger than 9 and flux error is smaller than + # flux then append flux plots ... if ts > 9.0 and e_flx < flx: - - # Add information energies.append(c_energy.real(row)) flux.append(c_flux.real(row)) ed_engs.append(c_ed.real(row)) eu_engs.append(c_eu.real(row)) e_flux.append(c_eflux.real(row)) - # + # ... otherwise append upper limit else: - - # Add information ul_energies.append(c_energy.real(row)) ul_flux.append(c_upper.real(row)) ul_ed_engs.append(c_ed.real(row)) @@ -91,8 +94,39 @@ if __name__ == "__main__": plt.figure() plt.loglog() plt.grid() - plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], fmt='ro') - plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], yerr=1.0e-11, uplims=True, fmt='ro') - plt.xlabel("Energy (TeV)") - plt.ylabel(r"E dN/dE (erg cm$^{-2}$ s$^{-1}$)") - plt.show() + plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], + fmt='ro') + plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], + yerr=1.0e-11, uplims=True, fmt='ro') + plt.xlabel('Energy (TeV)') + plt.ylabel(r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)') + + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Print usage information + usage = 'Usage: show_spectrum.py filename [file]' + if len(sys.argv) < 2: + print(usage) + sys.exit() + + # Extract parameters + filename = sys.argv[1] + plotfile = '' + if len(sys.argv) == 3: + plotfile = sys.argv[2] + + # Show spectrum + show_spectrum(filename, plotfile=plotfile) diff --git a/examples/show_ts_distribution.py b/examples/show_ts_distribution.py index ef07b40..13f580f 100755 --- a/examples/show_ts_distribution.py +++ b/examples/show_ts_distribution.py @@ -1,12 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays the TS distribution values generated using cstsdist. +# Display TS distribution generated by cstsdist # -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2016 Jurgen Knodlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -22,20 +18,42 @@ # along with this program. If not, see . # # ========================================================================== -import matplotlib.mlab as mlab -import matplotlib.pyplot as plt -import numpy as np import sys import csv import math +try: + import matplotlib.pyplot as plt + import matplotlib.mlab as mlab + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() +try: + import numpy as np +except: + print('This script needs the "numpy" module') + sys.exit() # ==================== # # Read TS distribution # # ==================== # -def read_ts(filename, tsname="TS"): +def read_ts(filename, tsname='TS'): """ - Read TS distribution values from CSV file. + Read TS distribution values from CSV file + + Parameters + ---------- + filename : str + Name of CSV file + tsname : str, optional + Column name of TS + + Returns + ------- + a : `numpy.array` + Array of TS values """ # Initialise list values = [] @@ -53,9 +71,9 @@ def read_ts(filename, tsname="TS"): try: index = row.index(tsname) except: - sys.stdout.write('ERROR: Column "TS" not found in file\n') - sys.stdout.write(row+"\n") - raise NameError("TS") + print('ERROR: Column "'+tsname+'" not found in file') + print(row) + raise NameError(tsname) # Handle data rows else: @@ -76,7 +94,17 @@ def read_ts(filename, tsname="TS"): # =========== # def erf(x): """ - Compute error function. + Compute error function + + Parameters + ---------- + x : float + Value + + Returns + ------- + erf(-x) : float + Error function """ # save the sign of x sign = 1 @@ -103,9 +131,20 @@ def erf(x): # =============================== # # Display cumulative distribution # # =============================== # -def dist_cdf(values, nbins, title): +def dist_cdf(values, nbins, title, expect='', plotfile=''): """ - Displays cumulative TS distribution and overplots expectation. + Display cumulative TS distribution and overplots expectation + + Parameters + ---------- + values : float + TS values + nbins : int + Number of bins + title : str + Plot title + plotfile : str, optional + Plot file name """ # Set range and adapt number of bins. We make sure that we have a # histogram bin centred on 0 that will capture the TS=0 part. @@ -142,19 +181,23 @@ def dist_cdf(values, nbins, title): hist_expect.append(norm) # Show expected distribution (half) - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations") - l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation") - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) + plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label='simulations') + if expect == 'show': + plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label='expectation') + plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) # Set plot plt.xlabel('TS') plt.ylabel('Number of trials > TS') plt.title(title) plt.grid(True) - plt.legend(loc="upper right") + plt.legend(loc='upper right') - # Show histogram - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() # Return return @@ -163,10 +206,20 @@ def dist_cdf(values, nbins, title): # ======================================== # # Display probability density distribution # # ======================================== # -def dist_pdf(values, nbins, title): +def dist_pdf(values, nbins, title, expect='', plotfile=''): """ - Display probability density TS distribution and overplots - expectation. + Display probability density TS distribution and overplots expectation + + Parameters + ---------- + values : float + TS values + nbins : int + Number of bins + title : str + Plot title + plotfile : str, optional + Plot file name """ # Set range and adapt number of bins. We make sure that we have a # histogram bin centred on 0 that will capture the TS=0 part. @@ -180,9 +233,9 @@ def dist_pdf(values, nbins, title): # Create histogram n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], align='mid', facecolor='green', - log=True, label="simulations") + log=True, label='simulations') - # Create expected distribution (onyl for positive TS). We compute the + # Create expected distribution (only for positive TS). We compute the # full and half of the distribution, as for positively constrained # quantities half of the values have TS=0 x = [] @@ -202,17 +255,24 @@ def dist_pdf(values, nbins, title): y_half.append(0.5*y) # Show expected distribution (half) - l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation") + if expect == 'show': + plt.semilogy(x, y_half, 'ro', linewidth=2, label='expectation') # Set plot plt.xlabel('TS') plt.ylabel('Number of trials') plt.title(title) plt.grid(True) - plt.legend(loc="upper right") + plt.legend(loc='upper right') - # Show histogram - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return # ======================== # @@ -221,26 +281,22 @@ def dist_pdf(values, nbins, title): if __name__ == '__main__': # Print usage information - usage = "Usage: show_ts_distribution filename" \ - " [-n bins] [-c column] [-t title] [-p plot]" + usage = 'Usage: show_ts_distribution.py filename' \ + ' [-n bins] [-c column] [-t title] [-p plot] [-f file]' if len(sys.argv) < 2: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Extract parameters filename = sys.argv[1] - # Set default parameters - nbins = 30 - tsname = "TS" - title = "TS distribution" - plot = "pdf" - # Parameter dictionnary - pars = [{'option': '-n', 'value': nbins}, - {'option': '-c', 'value': tsname}, - {'option': '-t', 'value': title}, - {'option': '-p', 'value': plot}] + pars = [{'option': '-n', 'value': 30}, + {'option': '-c', 'value': 'TS'}, + {'option': '-t', 'value': 'TS distribution'}, + {'option': '-p', 'value': 'pdf'}, + {'option': '-e', 'value': ''}, + {'option': '-f', 'value': ''}] # Gather parameters from command line i = 2 @@ -254,27 +310,29 @@ if __name__ == '__main__': try: par['value'] = sys.argv[i] except: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() else: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Next item i += 1 # Recover parameters - nbins = int(pars[0]['value']) - tsname = pars[1]['value'] - title = pars[2]['value'] - plot = pars[3]['value'] + nbins = int(pars[0]['value']) + tsname = pars[1]['value'] + title = pars[2]['value'] + plot = pars[3]['value'] + expect = pars[4]['value'] + plotfile = pars[5]['value'] # Read values from CSV file values = read_ts(filename, tsname=tsname) - sys.stdout.write(len(values), "values read.\n") + print('%d values read.' % len(values)) # Show histogram - if plot == "pdf": - dist_pdf(values, nbins, title) + if plot == 'pdf': + dist_pdf(values, nbins, title, expect=expect, plotfile=plotfile) else: - dist_cdf(values, nbins, title) + dist_cdf(values, nbins, title, expect=expect, plotfile=plotfile) diff --git a/test/README.md b/test/README.md index ca1a5ea..c93c30e 100644 --- a/test/README.md +++ b/test/README.md @@ -1,6 +1,4 @@ - ctools unit test ================ This folder contains files used for ctools unit testing. -test_cscripts.sh diff --git a/test/data/sensitivity.dat b/test/data/sensitivity.dat new file mode 100644 index 0000000..8872926 --- /dev/null +++ b/test/data/sensitivity.dat @@ -0,0 +1,11 @@ +loge,emin,emax,crab_flux,photon_flux,energy_flux,sensitivity +-1.51402150412,0.02,0.0468734582318,0.144989163852,6.60543533444e-10,3.05652923293e-11,3.56378721075e-11 +-1.14412450369,0.0468734582318,0.109856054331,0.0561134683318,7.24745374549e-11,7.85976500012e-12,9.16416230711e-12 +-0.774227503252,0.109856054331,0.257466658709,0.0315037806261,1.15354211075e-11,2.93193651747e-12,3.41851723554e-12 +-0.404330502818,0.257466658709,0.603417633655,0.0244879762385,2.54200331885e-12,1.51423828719e-12,1.76553948308e-12 +-0.0344335023848,0.603417633655,1.41421356237,0.022256802649,6.54995675049e-13,9.14436478628e-13,1.0661952755e-12 +0.335463498049,1.41421356237,3.31445401734,0.0250217133244,2.08758967789e-13,6.83057479879e-13,7.96416891675e-13 +0.705360498482,3.31445401734,7.76799609716,0.042333217441,1.00129446033e-13,7.67840339358e-13,8.95270214276e-13 +1.07525749892,7.76799609716,18.2056420303,0.0929194773185,6.23074149979e-14,1.11981387734e-12,1.30565686449e-12 +1.44515449935,18.2056420303,42.6680700645,0.14972962577,2.8463818009e-14,1.19893673197e-12,1.39791085453e-12 +1.81505149978,42.6680700645,100.0,0.382388032371,2.06082989569e-14,2.03442769003e-12,2.37205890421e-12 diff --git a/test/data/spectrum.fits b/test/data/spectrum.fits new file mode 100644 index 0000000..80da2e5 --- /dev/null +++ b/test/data/spectrum.fits diff --git a/test/data/ts.dat b/test/data/ts.dat new file mode 100644 index 0000000..573ca7f --- /dev/null +++ b/test/data/ts.dat @@ -0,0 +1,101 @@ +TS,LogL_bgm,LogL_all,Nevents,Npred_bkg,Npred_all,Crab_Prefactor,Unc_Crab_Prefactor,CTABackgroundModel_Prefactor,Unc_CTABackgroundModel_Prefactor,CTABackgroundModel_Index,Unc_CTABackgroundModel_Index +21858.52363364247,164117.78313692164,153188.5213201004,22947.0,22946.99999940957,22946.999182570653,5.793209910482599e-16,9.843787551860997e-18,0.9961007850466849,0.01197837072408249,-0.002988471615299611,0.007462197925324776 +22254.447332882788,165680.6922036114,154553.46853717,23099.0,23098.999999613512,23098.998093109363,5.834991758798705e-16,9.870580046002392e-18,1.0126678419515491,0.011967610062684665,0.004782477384744285,0.007316347479857281 +21669.561687524547,165499.35987704716,154664.57903328488,23134.0,23133.99999969615,23133.999315025623,5.759131518978944e-16,9.8078414532141e-18,1.0105363690914626,0.011967661087729033,-0.00017543331960004447,0.0073220367405137776 +20586.996292978816,164330.2704312009,154036.77228471148,22897.0,22896.999999836524,22896.997968569693,5.555266953291272e-16,9.65303499566792e-18,1.0126646529020689,0.011972560615151874,0.005782431464096715,0.007323711254845693 +21217.911190702056,163471.7261244653,152862.77052911426,22869.0,22868.999999925014,22868.99897334564,5.70225985492044e-16,9.781422622377391e-18,0.9910189865469351,0.011911721710682954,-0.006124951668038474,0.007435839123654169 +21499.99581397162,164613.52000714888,153863.52210016307,22972.0,22971.999999709296,22971.999752804753,5.6910222413269675e-16,9.758624936302264e-18,0.996047801086967,0.011904474629350275,-0.006618987374187345,0.007377903239793396 +21218.138624030922,164860.51429341873,154251.44498140327,23028.0,23027.999999721134,23027.999861910866,5.631722425238875e-16,9.709038606307067e-18,1.0013335385708708,0.011894934859389805,-0.006251748911501181,0.007319737866598979 +21516.115901010286,164574.57805417257,153816.52010366743,22977.0,22976.999999663327,22976.99992616064,5.737703984639025e-16,9.8058501757969e-18,1.0017304453830707,0.011926394544762287,-0.0012345022004066025,0.007362686667357548 +21814.264059167763,163757.54297349413,152850.41094391025,22865.0,22864.999999773583,22864.999780063565,5.781984731779256e-16,9.833138684097321e-18,0.9972279507003987,0.01184624889933737,0.0009075599364746032,0.007331742276660691 +21722.886111398868,165089.47672604493,154228.0336703455,23109.0,23108.999999674452,23108.9997904007,5.709992445993981e-16,9.76964807969154e-18,0.9967680654306393,0.011876107200647554,-0.011021721057770708,0.007332464028069347 +21598.863410390506,165201.15784713515,154401.7261419399,23121.0,23120.999999727013,23120.99956818323,5.796621137864304e-16,9.85685177987898e-18,1.0005077278386079,0.011874257987756785,-0.006437475225365512,0.007308543239518941 +21794.114624339214,164895.88104651406,153998.82373434445,23001.0,23000.99999958087,23000.999482800085,5.798901936139489e-16,9.862163118374043e-18,1.005069551247296,0.011900250356816858,0.0019410618022905287,0.0073140103821301514 +21227.425172272255,166038.31113511338,155424.59854897726,23106.0,23105.999999409567,23105.996163342123,5.683306232669589e-16,9.767724432554512e-18,1.0208070375514615,0.012000389287969418,0.006896248020332655,0.007278861707651359 +21573.883430126647,163893.0064978525,153106.06478278918,22893.0,22892.999999855037,22892.997657277247,5.746330138649936e-16,9.809856798510093e-18,0.9887958563489442,0.011919769035438853,-0.007702982005473525,0.007460539914015142 +21216.534766604716,164024.07526321124,153415.80787990888,22891.0,22890.99999976823,22890.999903703054,5.634003761042508e-16,9.713723962581571e-18,0.9988843468354479,0.011884925398741884,-0.002623626233923857,0.007345635611935417 +21597.08751892927,165975.3560595174,155176.81230005276,23150.0,23149.99999944865,23149.998888937902,5.718580414582826e-16,9.78000524693473e-18,1.0132558638441092,0.011992693047890666,0.0002528999590185808,0.007321935314776597 +21387.703713911586,163364.79937473897,152670.94751778318,22794.0,22793.999999815976,22793.999973829385,5.6892940405553605e-16,9.758446104823037e-18,0.9943975957184543,0.011824486448349799,-0.0008002832127617221,0.007333535186744982 +20912.907039812184,162038.81793470908,151582.36441480299,22666.0,22665.99999997628,22665.990260466566,5.601818185889693e-16,9.694320601561969e-18,0.9751815451597996,0.011793838389151424,-0.01294140924200319,0.007455055791860874 +21911.95291238872,163189.5450744877,152233.56861829333,22810.0,22809.999999897966,22809.998057095352,5.791145848008091e-16,9.84239549367057e-18,0.9854579489782141,0.011857156203560955,-0.005834836728800691,0.007437257754506665 +21121.265116654686,164055.15655683712,153494.52399850977,22912.0,22911.999999787822,22911.99984896773,5.654325050779409e-16,9.741375490234794e-18,0.9999753099585894,0.01194854296456228,-0.0021040061041819437,0.00739916979858624 +21195.306229041773,164281.77518035137,153684.12206583048,23034.0,23033.999999960743,23033.998176414847,5.70433524603058e-16,9.77671925031924e-18,0.9886794379131255,0.011878644954228475,-0.01446464700955386,0.007402981305713397 +20846.777400515275,162770.057188395,152346.66848813737,22731.0,22730.999999907664,22730.99981226707,5.624363678525482e-16,9.712847646970035e-18,0.9909653476365609,0.011763408213370439,-0.002595590722067562,0.007301229915792391 +20897.083739825,164019.2197518424,153570.6778819299,23003.0,23002.999999907483,23002.999556953055,5.649662099748417e-16,9.736547752010038e-18,0.9935802602295315,0.01188380653681554,-0.010805113941356031,0.007371306800476196 +20602.686458923854,164911.49560107855,154610.15237161663,22935.0,22934.99999955501,22934.99341642103,5.530693214901318e-16,9.625543333068213e-18,1.0206021943272432,0.011985188441572748,0.00972446272247772,0.00727243068453318 +21836.23271708103,163824.30785453672,152906.1914959962,22755.0,22754.999999495434,22754.999710735854,5.75568768540384e-16,9.811525606728942e-18,1.0043798388913345,0.01188200956164694,0.010276896104008231,0.007325165400710747 +21420.40131879633,164409.64030012718,153699.43964072902,23001.0,23000.999999780248,23000.99928771896,5.684570967501337e-16,9.762208897347998e-18,0.9913746713865513,0.011878388317921098,-0.011560419259827103,0.007384451973090443 +21727.81684923329,164232.132445953,153368.22402133635,22985.0,22984.99999976957,22984.998900142164,5.738813127112364e-16,9.809036890551794e-18,0.9911339994790322,0.01189255974895351,-0.009747052205701144,0.00740674910573882 +20868.72095527401,164301.13159811447,153866.77112047747,22992.0,22991.999999802247,22991.99969008748,5.605800164321631e-16,9.703156512461817e-18,0.9973271760042328,0.011871321613801066,-0.008569072284481264,0.0073309043881014035 +22492.009746829106,164682.28638778496,153436.2815143704,22951.0,22950.9999996003,22950.99665538983,5.879507535878784e-16,9.910190021652269e-18,1.001162832181959,0.011846315465073663,0.002987242954331668,0.0073024800531300005 +22653.864335719903,165519.4792556914,154192.54708783145,23076.0,23075.999999249838,23075.993670520074,5.912861048051989e-16,9.940531531808007e-18,1.0024002934241467,0.011941973255326052,-0.00023952533210942385,0.007373252448619404 +21449.306879230775,166768.10502219904,156043.45158258366,23297.0,23296.99999917116,23296.992123045784,5.7341518582148665e-16,9.807974620497261e-18,1.024760753033504,0.012028256916479971,0.0035902945030410723,0.007258563117465229 +21486.346823231317,164876.94234162793,154133.76893001227,23003.0,23002.999999741518,23002.999859281495,5.757455861467064e-16,9.822563593405131e-18,1.004168442832293,0.011942314003147688,0.00011042172295081477,0.0073583821510989885 +20652.239219280833,163809.36895235302,153483.2493427126,22957.0,22956.999999943808,22956.995521902107,5.528577248898855e-16,9.622644416274187e-18,0.9853639781514777,0.011834892876956514,-0.018391417942683855,0.0073803794078678014 +21787.521341275133,164528.73527216076,153634.9746015232,22991.0,22990.99999973,22990.998358757795,5.720325280614957e-16,9.785480546850936e-18,0.9855686593361147,0.011839563015884493,-0.014781546186952868,0.007392004202286942 +21402.38363453705,164305.6816226261,153604.4898053576,22994.0,22993.99999989152,22993.99965618869,5.712540995621592e-16,9.778682761426494e-18,0.9949474763637698,0.011898956019282798,-0.007806261712544803,0.007380964683300804 +20141.589426982973,163856.60168552134,153785.80697202985,22929.0,22928.999999913434,22928.999999743348,5.438078945642788e-16,9.558024400005279e-18,0.9932852802680611,0.011875856351872056,-0.01340284403577944,0.007359275829012058 +21943.395299219526,163065.74164094013,152094.04399133037,22784.0,22783.999999305994,22783.99675835668,5.79676765978459e-16,9.848370996477712e-18,0.9837834848183966,0.01187815045838842,-0.005954920583608014,0.007472494212571003 +20834.792521972733,164064.23789739315,153646.84163640678,22895.0,22894.99999975629,22894.999530192188,5.589302026667169e-16,9.677416099688684e-18,0.9986526033294748,0.011854084906733467,-0.004096128598782811,0.007311666493293262 +21305.324882587476,163265.06469459928,152612.40225330554,22828.0,22827.999999891683,22827.99808304598,5.713178115397861e-16,9.79512015984894e-18,0.987396852431259,0.011894410275773492,-0.007029367080902169,0.007451877733083311 +20847.034227716562,163814.44659629313,153390.92948243485,22897.0,22896.9999996484,22896.99953025752,5.607597257232211e-16,9.688770805376503e-18,1.0047595349843508,0.011799861066908465,0.0010146013322702991,0.007215496921685529 +21216.92292083666,164176.47242941064,153568.0109689923,22999.0,22998.999999831543,22998.9975955996,5.642001984733995e-16,9.72351332200732e-18,0.9899066275912457,0.011903062056425938,-0.013687904197762802,0.007417248945029567 +20816.64541355765,163181.06644880498,152772.74374202616,22777.0,22776.99999988913,22776.994775080595,5.516451209091623e-16,9.618999946029458e-18,0.9869057530226409,0.011856380365721023,-0.010358640899461943,0.007409856539614678 +21679.522176688828,163799.51890170656,152959.75781336214,22837.0,22836.99999978841,22836.99973823716,5.717636386575965e-16,9.776297868466177e-18,0.9908230482181739,0.011847216840415472,-0.0046007921454711475,0.007378504976736998 +21106.35513021017,164821.97117069492,154268.79360558983,22917.0,22916.999999453958,22916.99841731256,5.647103860714512e-16,9.727705266059672e-18,1.01541521494636,0.011955291437595085,0.00943775316656054,0.007291252547122076 +20837.813223739155,163866.034636364,153447.12802449442,22791.0,22790.999999713706,22790.998309820432,5.572360798379893e-16,9.654736904973028e-18,1.0091814893739082,0.011893607890048716,0.007811373787627806,0.007283734723217876 +20780.997344649804,164108.8221487588,153718.3234764339,22982.0,22981.999999861746,22981.99986802358,5.630081040612143e-16,9.720705407249754e-18,0.9960740698856181,0.011864470017298944,-0.008531796083847434,0.007334719933560889 +21835.693036562065,164481.7513907408,153563.90487245977,22863.0,22862.999999581152,22862.999821485086,5.719539799575279e-16,9.779590491901723e-18,1.0043623557350347,0.011875625884065816,0.0049348547201782614,0.007305172253576122 +21653.10258389666,165542.21112609853,154715.6598341502,23262.0,23261.999999841548,23261.99499680352,5.749522283549226e-16,9.819619041947758e-18,0.9902674093404891,0.011918436489673211,-0.02105969680201884,0.007407329162175303 +21846.154768515087,164624.27857300622,153701.20118874867,23012.0,23011.99999979313,23011.998528939348,5.833817292347365e-16,9.884966321550659e-18,0.9935369248121374,0.011882112655646851,-0.006556706861582095,0.007382216624463548 +22013.17601596308,166235.52968944985,155228.9416814683,23227.0,23226.999999677897,23226.998673585054,5.797707114245376e-16,9.848651501252573e-18,1.006518546940583,0.011988963639015245,-0.00598127326508908,0.007364840449676914 +21651.91926104814,163154.12158821282,152328.16195768875,22773.0,22772.99999985566,22772.998844177993,5.714396502849746e-16,9.769436981298296e-18,0.9851000302656024,0.011819450107849325,-0.006582743622206519,0.007400947537603796 +21989.977012775722,165301.46496142616,154306.4764550383,23054.0,23053.999999486277,23053.99866206013,5.831451508774712e-16,9.873936836478257e-18,1.0070990373467372,0.011907048508928606,0.0022074936263132647,0.007300838328239046 +21438.695088208653,163698.03314021966,152978.68559611533,22850.0,22849.999999890217,22849.999942469298,5.723101116737367e-16,9.79415042715683e-18,0.9959990059596918,0.011909255151422262,-0.0009492657684907494,0.007402029667490372 +20838.374803555547,163909.36890479646,153490.18150301868,22825.0,22824.999999715288,22824.99939814383,5.603315148685522e-16,9.69540670910627e-18,1.0087802770334235,0.011925882318992111,0.006917357818017513,0.007318975346866552 +21270.20557162707,165625.21774982946,154990.11496401593,23098.0,23097.999999612013,23097.99846084131,5.695134364689924e-16,9.766410999308169e-18,1.0165267626116348,0.012010348326125827,0.0042397629393144945,0.007319288376668166 +21527.705370273616,164243.86907532276,153480.01639018595,22965.0,22964.999999743457,22964.99993603194,5.701748498974925e-16,9.771290624802224e-18,0.995702026500132,0.011874291328625755,-0.006337078712909612,0.007353017018362593 +21055.645831402973,164272.05201691244,153744.22910121095,22913.0,22912.99999970991,22912.999921720784,5.630081490231017e-16,9.7255367108198e-18,0.9979621989156474,0.011868129706833142,-0.0043189322848428005,0.0073321006440042955 +21156.85824884911,162399.58053200433,151821.15140757977,22706.0,22705.99999991407,22705.997403443518,5.654797040066863e-16,9.740490054150095e-18,0.9784619159402815,0.011744743142957838,-0.010617863436450973,0.007382073296970704 +21201.688190487737,164789.7753571481,154188.93126190422,23026.0,23025.99999906124,23025.99972065303,5.667600645790645e-16,9.742584040372483e-18,1.0120410866982648,0.01207095002314141,0.002979507608702241,0.007418035938329006 +21579.612940530176,164202.44370278623,153412.63723252114,22874.0,22873.999999804666,22874.00001950805,5.716335760382397e-16,9.78228682145465e-18,1.002323437487339,0.011916380859040334,0.0028258681377384813,0.007357128369306409 +21460.62560378993,164011.29522604778,153280.9824241528,22951.0,22950.99999985298,22950.999919398004,5.718007013749168e-16,9.78668674780352e-18,0.9941518980348495,0.011860449976759028,-0.006570586126134151,0.007354999480418171 +20888.288302358415,164286.87337022196,153842.72921904275,22976.0,22975.99999977232,22975.998047214835,5.544737396797197e-16,9.641570891212174e-18,1.0007779663248277,0.011916634775209257,-0.006801761329085501,0.0073460207848977796 +20408.074956183205,164422.41163892217,154218.37416083057,22955.0,22954.999999835327,22954.995893393418,5.513465666103151e-16,9.625616231421954e-18,1.0058950024352238,0.011974149384797315,-0.0028091604581978514,0.007363395502675995 +21802.37343549641,164500.592398738,153599.4056809898,22961.0,22960.999999757954,22960.9989212055,5.832468881001641e-16,9.890309949834711e-18,1.0034751072565908,0.011921220309959827,0.0031744106933537897,0.0073533366536048545 +20876.450997727632,164582.49714246154,154144.27164359772,22993.0,22992.999999842265,22992.9998257974,5.648090586661093e-16,9.736675647490638e-18,1.0025225514093496,0.011933600093706552,-0.0035293846222004244,0.007354639908367345 +21706.519023373083,164473.15900975114,153619.8994980646,22978.0,22977.999999763233,22977.999793557632,5.759940859191368e-16,9.813811252134857e-18,0.9984861958740264,0.011920224946766407,-0.0032243711331649107,0.007381393296644344 +21656.740131683822,165712.62332411483,154884.25325827292,23262.0,23261.99999977213,23261.998757088804,5.797595166213485e-16,9.856541718951383e-18,0.9985566042942451,0.011918622406145106,-0.013482108055460828,0.007351626660765964 +21108.65355328942,164601.01057088916,154046.68379424445,23016.0,23015.99999983359,23015.999758528997,5.653188521209414e-16,9.736939795430797e-18,0.9959645629245217,0.011889610707927175,-0.00938547147416546,0.007358851980664596 +20647.59149535757,162669.5612982715,152345.7655505927,22734.0,22733.99999995898,22733.995853601802,5.5190264735296685e-16,9.620164290713148e-18,0.9903309324372139,0.01187206487521323,-0.005898045252815988,0.007403049951329795 +22290.401570641203,163883.56355593016,152738.36277060956,22782.0,22781.999999546068,22781.99916149079,5.824528487203314e-16,9.864372038097642e-18,0.9966719584506485,0.011871832280126494,0.004939494731895152,0.007375380023497896 +21624.4816204338,164856.28447951959,154044.0436693027,22996.0,22995.99999963064,22995.999887018945,5.729990754258132e-16,9.793654933564995e-18,1.0001246408759745,0.011920784556020848,-0.003433668355429103,0.007365033404665828 +21599.218091811636,162947.49765919705,152147.88861329123,22666.0,22665.999999721185,22665.9995259345,5.76124788921756e-16,9.820444219997191e-18,0.9997434910652248,0.011837575389850108,0.010453639932065167,0.0073277935176205644 +21352.83466984972,163983.0417331816,153306.62439825674,22833.0,22832.999999748183,22832.999677072614,5.611046191249927e-16,9.692013011636335e-18,0.9935935257673363,0.011855675126920649,-0.004991283544181841,0.007359215729362822 +22138.833458752255,164579.4216104209,153510.00488104476,22972.0,22971.99999954863,22971.999029847222,5.82889920739159e-16,9.860552383071182e-18,1.0016245840621145,0.011965267887535018,0.0012043877283222457,0.007408104140671816 +21460.73646459746,164025.76170268297,153295.39347038424,22919.0,22918.99999972083,22918.99998951889,5.739817859219702e-16,9.810075323312415e-18,0.9966380414703272,0.01185241980257522,-0.0028036440323137035,0.007332677333553906 +20272.70840003481,163913.44129798363,153777.08709796623,22905.0,22904.99999986819,22904.993807344785,5.491116393214731e-16,9.603700699255587e-18,0.9993535632641539,0.01192103384820017,-0.006429848120336535,0.007364414061681116 +22040.583799355198,164196.0694411831,153175.7775415055,22948.0,22947.999999762113,22947.998520762158,5.830729768890433e-16,9.878618700156061e-18,0.9889198731570649,0.01183479859649009,-0.00766945927540272,0.00737644490898403 +22189.423341836373,162999.79537715373,151905.08370623554,22762.0,22761.999999802018,22761.99812426319,5.840469854404995e-16,9.88689591899704e-18,0.9874701389124014,0.011840024598898219,-0.001045075500813165,0.007415569662489257 +21404.67148317158,163675.38150747615,152973.04576589036,22873.0,22872.99999985997,22872.999797008175,5.707784704943953e-16,9.776315626974606e-18,0.9932817805516895,0.011877174992673342,-0.004381770811720796,0.007385324275316955 +21680.787235579046,164409.20965859844,153568.81604080892,23000.0,22999.999999595042,22999.996041215003,5.756363373649316e-16,9.818959282992611e-18,0.9911211088710422,0.011963680831576883,-0.009911884410269314,0.007474654715810191 +21579.654693948105,164174.39769977125,153384.5703527972,22910.0,22909.999999767482,22909.999826840383,5.71987722050035e-16,9.791834857909584e-18,0.9969717770530057,0.011930237932024544,-0.0026897435561697674,0.007407652576797861 +21297.579771767836,163385.8273468829,152737.03746099898,22844.0,22843.999999879547,22843.9991035317,5.688808155458493e-16,9.7608380955326e-18,0.987985497494533,0.011852158347747057,-0.007831002182125445,0.007402005571818617 +22413.117900692392,165160.81083219431,153954.25188184812,23075.0,23074.99999963741,23074.9980466868,5.849262030225733e-16,9.88285352545655e-18,0.9941410946295832,0.011801812232478627,-0.008211194148192122,0.007291859455219944 +21122.470503843215,163334.07738528133,152772.84213335972,22879.0,22878.999999898744,22878.99709745025,5.633723709743721e-16,9.722100961324131e-18,0.984984844375672,0.011857253919565828,-0.01296529513990941,0.007421939723453924 +21282.707195805793,163941.15171670247,153299.79811879958,22909.0,22908.9999998743,22908.99996153768,5.729219609420647e-16,9.797812724663524e-18,0.9954654665694013,0.011861939156832611,-0.003583079120707077,0.007352470265065281 +20724.34760198067,163318.72421718622,152956.55041619588,22790.0,22789.999999879707,22789.99813018993,5.550590116159096e-16,9.645783754938671e-18,0.9933982950272806,0.01187452330418287,-0.00495332364849543,0.0073808910696461375 +21757.13207008707,165179.7837367039,154301.21770166035,23030.0,23029.999999515534,23029.999095339273,5.763186525066147e-16,9.828374022376591e-18,1.0136834107892345,0.011971605614733427,0.006510480574850554,0.007314863606454032 +21671.1096591925,166527.63624669318,155692.08141709692,23264.0,23263.99999942336,23263.992478841574,5.84129987251449e-16,9.896922008537031e-18,1.0234213270128842,0.012014512018187535,0.006575992519346971,0.007267333319003582 +21221.41757199116,163903.57706571487,153292.8682797193,22906.0,22905.999999861677,22905.998783516217,5.646919165369301e-16,9.724700045763518e-18,0.9900591679883312,0.011885002043500589,-0.009753133051711622,0.007408431197796472 +21681.41252249159,165556.59745661478,154715.89119536898,23196.0,23195.999999720592,23195.99929739878,5.773631747722166e-16,9.836520659176885e-18,0.9977101559045379,0.01186422662665964,-0.012136729410717443,0.0073081040770472065 +21153.477707419137,166277.3899099234,155700.65105621383,23255.0,23254.999999606127,23254.99604384548,5.653184419662256e-16,9.730321230819675e-18,1.019156144918438,0.01201068815115771,-0.0010363684582943512,0.007280823299749718 +22226.132660409086,164336.5451073256,153223.47877712105,22935.0,22934.999999669544,22934.9976589201,5.8647466247030845e-16,9.907864241886007e-18,0.9988210866347508,0.011904203335735244,0.001435599985409938,0.0073741953380605406 +21424.358322370157,164683.5971306915,153971.4179695064,22940.0,22939.99999953921,22939.999271899742,5.679278962446394e-16,9.75563702408002e-18,1.0098384925060064,0.011899202395887487,0.0050329170142471,0.007276269049906902 +21455.353563566052,164947.2565606395,154219.57977885648,22987.0,22986.99999956718,22986.998933887335,5.706645322939917e-16,9.770262988270668e-18,1.0142244697962193,0.011960235218104793,0.007217238595620098,0.007300805026314625 +21452.114580445574,164273.43465029492,153547.37736007213,22845.0,22844.9999995031,22844.999775048265,5.677909293087916e-16,9.751768271419505e-18,1.0082800787052053,0.011922505792565435,0.007632916428577924,0.007322032069579832 +21019.69961927872,162935.32861177463,152425.47880213527,22657.0,22656.999999781452,22656.999043407268,5.573701521579324e-16,9.661153343333429e-18,0.994168966337365,0.011892058320532648,0.0016080034919236333,0.007408684998744538 +21441.518593108456,164053.07922474376,153332.31992818954,22960.0,22959.99999984039,22959.999974187285,5.727077800111475e-16,9.793315534540368e-18,0.9933895462403782,0.011810199688400778,-0.0072931345463456265,0.007310683986440049 +21781.365669648803,164462.0895137939,153571.4066789695,22877.0,22876.999999410484,22876.999794261017,5.738537591812932e-16,9.802523721098818e-18,1.0067930568779322,0.0119066595609444,0.006739517264159702,0.007318677995153281 +21013.571768569178,162221.01090832267,151714.22502403808,22714.0,22713.999999970183,22714.000000122807,5.578526953798294e-16,9.670657005596775e-18,0.96868333949408,0.011813808188693846,-0.02057929881999145,0.007517120515436912 +20965.61130349373,164112.7074564671,153629.90180472023,22893.0,22892.999999919695,22892.99927046633,5.5999754133615e-16,9.687391521031491e-18,0.9954549917387471,0.011910503281561615,-0.0062261048953279885,0.007392491292287055 diff --git a/test/dev/README.md b/test/dev/README.md new file mode 100644 index 0000000..0de5415 --- /dev/null +++ b/test/dev/README.md @@ -0,0 +1,30 @@ +Developer test scripts +====================== + +This folder contains scripts that are used for various testing purposes +during code development. None of the scripts actually needs to work, they +may actually be outdated, but could still be useful for further testing +if needed. + +`benchmark_cta_analysis.py` +> This script peforms a benchmark of the various CTA analysis types +> (unbinned, binned, stacked), and if matplotlib is installed, creates +> a plot of the benchmark results. + +`check_models.py` +> Checks a variety of models using an unbinned analysis. + +`make_pull_at_sensitivity_limit.py` +> This script creates pull distributions as function of energy +> for a test source that is at the 5 sigma sensitivity limit for +> 50h of effective observing time. It takes as input file the +> output file generated using cspull. +> If the Python processing module is installed the script can make +> use of multiple cores and runs in parallel. + +`make_survey.py` +> This script simulates a CTA galactic plane survey. + +`simulate_events.py` +> Simulates events and compares the simulation to the model. + diff --git a/test/test_examples.py b/test/test_examples.py index 4418ac8..21032ee 100755 --- a/test/test_examples.py +++ b/test/test_examples.py @@ -104,6 +104,10 @@ class Test(gammalib.GPythonTestSuite): self.append(self.test_show_pha, 'Test show_pha.py') self.append(self.test_show_pull_evolution, 'Test show_pull_evolution.py') self.append(self.test_show_pull_histogram, 'Test show_pull_histogram.py') + self.append(self.test_show_response, 'Test test_show_response.py') + self.append(self.test_show_sensitivity, 'Test test_show_sensitivity.py') + self.append(self.test_show_spectrum, 'Test test_show_spectrum.py') + self.append(self.test_show_ts_distribution, 'Test show_ts_distribution.py') # Return return @@ -253,7 +257,8 @@ class Test(gammalib.GPythonTestSuite): """ # Execute script self._execute_python('show_butterfly', - args='data/butterfly.txt example_butterfly.eps') + args='data/butterfly.txt ' + 'example_show_butterfly.eps') #TODO: Do any testing @@ -267,7 +272,7 @@ class Test(gammalib.GPythonTestSuite): """ # Execute script self._execute_python('show_model', - args='data/crab.xml Crab example_model.eps') + args='data/crab.xml Crab example_show_model.eps') #TODO: Do any testing @@ -280,7 +285,8 @@ class Test(gammalib.GPythonTestSuite): Test show_pha """ # Execute script - self._execute_python('show_pha', args='data/pha.fits example_pha.eps') + self._execute_python('show_pha', + args='data/pha.fits example_show_pha.eps') #TODO: Do any testing @@ -295,7 +301,7 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('show_pull_evolution', args='data/pull.dat Pull_Crab_Prefactor ' - 'example_pull_evolution.eps') + 'example_show_pull_evolution.eps') #TODO: Do any testing @@ -310,7 +316,65 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('show_pull_histogram', args='data/pull.dat Pull_Crab_Prefactor 50 ' - 'example_pull_histogram.eps') + 'example_show_pull_histogram.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_response + def test_show_response(self): + """ + Test show_response + """ + # Execute script + self._execute_python('show_response', args='example_show_response.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_sensitivity + def test_show_sensitivity(self): + """ + Test show_sensitivity + """ + # Execute script + self._execute_python('show_sensitivity', + args='data/sensitivity.dat ' + 'example_show_sensitivity.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_spectrum + def test_show_spectrum(self): + """ + Test show_spectrum + """ + # Execute script + self._execute_python('show_spectrum', + args='data/spectrum.fits ' + 'example_show_spectrum.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_ts_distribution + def test_show_ts_distribution(self): + """ + Test show_ts_distribution + """ + # Execute script + self._execute_python('show_ts_distribution', + args='data/ts.dat -n 100 ' + '-f example_show_ts_distribution.eps') #TODO: Do any testing