Commit 8590d4867d39b8803424c7528d8d83cfc2728b3a
1 parent
ce0042b4
Complete adding of unit tests for all example scripts
Showing
20 changed files
with
611 additions
and
237 deletions
ChangeLog
| 1 | -2016-06-21 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu> | |
| 1 | +2016-06-25 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu> | |
| 2 | 2 | |
| 3 | 3 | * Version 1.1.0 released |
| 4 | 4 | ======================== |
| 5 | 5 | |
| 6 | + Add unit tests for all example scripts | |
| 6 | 7 | Handle zero parameter errors in ctulimit (#1795) |
| 7 | 8 | Compute counts cube weights in ctbin (#1789) |
| 8 | 9 | Adapt to new stacked cube interface (#1788) | ... | ... |
NEWS
README.md
examples/README.md
| 1 | 1 | ctools examples |
| 2 | 2 | =============== |
| 3 | 3 | |
| 4 | -This folder contains Python example scripts that illustrate the | |
| 5 | -usage of the ctools. | |
| 6 | - | |
| 7 | -`benchmark_cta_analysis.py` | |
| 8 | -> This script peforms a benchmark of the various CTA analysis types | |
| 9 | -> (unbinned, binned, stacked), and if matplotlib is installed, creates | |
| 10 | -> a plot of the benchmark results. | |
| 11 | - | |
| 12 | -`check_models.py` | |
| 13 | -> Checks a variety of models using an unbinned analysis. | |
| 4 | +This folder contains Python example scripts that illustrate the usage of | |
| 5 | +the ctools. | |
| 14 | 6 | |
| 15 | 7 | `generate_prod3_irfs/py` |
| 16 | -> Generate Prod3 Instrument Response Functions from FITS files. | |
| 8 | +> Generates Prod3 Instrument Response Functions from FITS files. | |
| 17 | 9 | |
| 18 | 10 | `make_pointings.py` |
| 19 | -> This script creates pointings for CTA surveys. For the moment the | |
| 20 | -> Galactic Plane Survey (gps), the Extragalactic Survey (extgal), | |
| 21 | -> the Galactic Centre Survey (gc) and the LMC survey (lmc) are | |
| 22 | -> implemented. The script outputs and ASCII file that can be used by | |
| 23 | -> csobsdef. | |
| 24 | - | |
| 25 | -`make_pull_at_sensitivity_limit.py` | |
| 26 | -> This script creates pull distributions as function of energy | |
| 27 | -> for a test source that is at the 5 sigma sensitivity limit for | |
| 28 | -> 50h of effective observing time. It takes as input file the | |
| 29 | -> output file generated using cspull. | |
| 30 | -> If the Python processing module is installed the script can make | |
| 31 | -> use of multiple cores and runs in parallel. | |
| 11 | +> Creates pointings for CTA surveys. For the moment the Galactic Plane | |
| 12 | +> Survey (gps), the Extragalactic Survey (extgal), the Galactic Centre | |
| 13 | +> Survey (gc) and the LMC survey (lmc) are implemented. The script outputs | |
| 14 | +> and ASCII file that can be used by csobsdef. | |
| 32 | 15 | |
| 33 | 16 | `make_spectrum.py` |
| 34 | -> This script shows how to make a spectrum using obsutils. | |
| 35 | - | |
| 36 | -`make_survey.py` | |
| 37 | -> This script simulates a CTA galactic plane survey. | |
| 17 | +> Simulates events and generates a spectrum from the using csspec. | |
| 38 | 18 | |
| 39 | 19 | `make_ts_distributions.py` |
| 40 | -> This script creates TS distributions as function of energy for | |
| 41 | -> an off observation of 50h effective observing time. | |
| 42 | -> If the Python processing module is installed the script can make | |
| 43 | -> use of multiple cores and runs in parallel. | |
| 20 | +> Creates TS distributions as function of energy for an OFF observaition | |
| 21 | +> of 50 hours effective observing time. If the Python processing module is | |
| 22 | +> installed the script can make use of multiple cores and runs in parallel. | |
| 44 | 23 | |
| 45 | 24 | `pipeline_binned_disk.py` |
| 46 | -> Example pipeline for simulating events and performing a binned | |
| 47 | -> analysis by saving intermediate results to disk. | |
| 25 | +> Example pipeline for simulating events and performing a binned analysis | |
| 26 | +> by saving intermediate results to disk. | |
| 48 | 27 | |
| 49 | 28 | `pipeline_binned_mem.py` |
| 50 | -> Example pipeline for simulating events and performing a binned | |
| 51 | -> analysis in memory. No files will be written to disk. | |
| 29 | +> Example pipeline for simulating events and performing a binned analysis | |
| 30 | +> in memory. No files will be written to disk. | |
| 52 | 31 | |
| 53 | 32 | `pipeline_stacked_disk.py` |
| 54 | -> Example pipeline for simulating events and performing a stacked | |
| 55 | -> analysis by saving intermediate results to disk. | |
| 33 | +> Example pipeline for simulating events and performing a stacked analysis | |
| 34 | +> by saving intermediate results to disk. | |
| 56 | 35 | |
| 57 | 36 | `pipeline_stacked_mem.py` |
| 58 | -> Example pipeline for simulating events and performing a stacked | |
| 59 | -> analysis in memory. No files will be written to disk. | |
| 37 | +> Example pipeline for simulating events and performing a stacked analysis | |
| 38 | +> in memory. No files will be written to disk. | |
| 60 | 39 | |
| 61 | 40 | `pipeline_unbinned_disk.py` |
| 62 | -> Example pipeline for simulating events and performing an unbinned | |
| 63 | -> analysis by saving intermediate results to disk. | |
| 41 | +> Example pipeline for simulating events and performing an unbinned analysis | |
| 42 | +> by saving intermediate results to disk. | |
| 64 | 43 | |
| 65 | 44 | `pipeline_unbinned_mem.py` |
| 66 | -> Example pipeline for simulating events and performing an unbinned | |
| 67 | -> analysis in memory. No files will be written to disk. | |
| 45 | +> Example pipeline for simulating events and performing an unbinned analysis | |
| 46 | +> in memory. No files will be written to disk. | |
| 68 | 47 | |
| 69 | 48 | `show_butterfly.py` |
| 70 | 49 | > Displays a butterfly generated using the ctbutterfly tool. |
| ... | ... | @@ -86,14 +65,10 @@ usage of the ctools. |
| 86 | 65 | > Displays one or several CTA response functions. |
| 87 | 66 | |
| 88 | 67 | `show_sensitivity.py` |
| 89 | -> Displays sensitivity information generated by the cssens script. | |
| 68 | +> Displays the sensitivity generated by the cssens script. | |
| 90 | 69 | |
| 91 | 70 | `show_spectrum.py` |
| 92 | 71 | > Displays spectrum generated by the csspec script. |
| 93 | 72 | |
| 94 | 73 | `show_ts_distribution.py` |
| 95 | -> Displays Test Statistic histograms generated using the cstsdist | |
| 96 | -> script. | |
| 97 | - | |
| 98 | -`simulate_events.py` | |
| 99 | -> Simulates events and compares the simulation to the model. | |
| 74 | +> Displays Test Statistic histograms generated using the cstsdist script. | ... | ... |
examples/make_ts_distributions.py
| 1 | 1 | #! /usr/bin/env python |
| 2 | 2 | # ========================================================================== |
| 3 | -# Generate TS distribution | |
| 3 | +# Generate TS distributions as function of energy for an OFF observation | |
| 4 | 4 | # |
| 5 | 5 | # Copyright (C) 2011-2016 Jurgen Knodlseder |
| 6 | 6 | # |
| ... | ... | @@ -65,7 +65,7 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, |
| 65 | 65 | tsdist['outfile'] = outfile |
| 66 | 66 | tsdist['ntrials'] = ntrials |
| 67 | 67 | tsdist['caldb'] = 'prod2' |
| 68 | - tsdist['irf'] = 'South_0.5h' | |
| 68 | + tsdist['irf'] = 'South_50h' | |
| 69 | 69 | tsdist['ra'] = 83.63 |
| 70 | 70 | tsdist['dec'] = 22.01 |
| 71 | 71 | tsdist['emin'] = emin |
| ... | ... | @@ -76,7 +76,6 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, |
| 76 | 76 | tsdist['rad'] = 5.0 |
| 77 | 77 | tsdist['npix'] = 200 |
| 78 | 78 | tsdist['binsz'] = 0.05 |
| 79 | - #tsdist['debug'] = True | |
| 80 | 79 | |
| 81 | 80 | # Optionally open the log file |
| 82 | 81 | if log: |
| ... | ... | @@ -95,7 +94,7 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, |
| 95 | 94 | if __name__ == '__main__': |
| 96 | 95 | |
| 97 | 96 | # Get input arguments |
| 98 | - usage = 'make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]' | |
| 97 | + usage = 'make_ts_distributions.py [-n ntrials] [-e enumbins] [-m max_threads]' | |
| 99 | 98 | if len(sys.argv) < 1: |
| 100 | 99 | print(usage) |
| 101 | 100 | sys.exit() | ... | ... |
examples/show_butterfly.py
| ... | ... | @@ -32,7 +32,7 @@ except: |
| 32 | 32 | # ============ # |
| 33 | 33 | # Script entry # |
| 34 | 34 | # ============ # |
| 35 | -if __name__ == "__main__": | |
| 35 | +if __name__ == '__main__': | |
| 36 | 36 | |
| 37 | 37 | # Check for given butterfly file |
| 38 | 38 | if len(sys.argv) < 2: |
| ... | ... | @@ -83,7 +83,7 @@ if __name__ == "__main__": |
| 83 | 83 | plt.plot(line_x,line_y,color='black',ls='-') |
| 84 | 84 | plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5) |
| 85 | 85 | plt.xlabel('Energy (MeV)') |
| 86 | - plt.ylabel('E dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)') | |
| 86 | + plt.ylabel(r'E $\times$ dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)') | |
| 87 | 87 | |
| 88 | 88 | # Show spectrum or save it into file |
| 89 | 89 | if len(plotfile) > 0: | ... | ... |
examples/show_model.py
| ... | ... | @@ -103,7 +103,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100, plotfile=''): |
| 103 | 103 | # ============= # |
| 104 | 104 | # Script entry # |
| 105 | 105 | # ============= # |
| 106 | -if __name__ == "__main__": | |
| 106 | +if __name__ == '__main__': | |
| 107 | 107 | |
| 108 | 108 | # Check for model filename |
| 109 | 109 | if len(sys.argv) < 2: | ... | ... |
examples/show_pha.py
| ... | ... | @@ -95,7 +95,7 @@ def show_pha(pha, plotfile=''): |
| 95 | 95 | if __name__ == '__main__': |
| 96 | 96 | |
| 97 | 97 | # Print usage information |
| 98 | - usage = 'Usage: show_pha filename [file]' | |
| 98 | + usage = 'Usage: show_pha.py filename [file]' | |
| 99 | 99 | if len(sys.argv) < 2: |
| 100 | 100 | print(usage) |
| 101 | 101 | sys.exit() | ... | ... |
examples/show_pull_evolution.py
| ... | ... | @@ -89,7 +89,7 @@ def read_pull(filename, parname): |
| 89 | 89 | if __name__ == '__main__': |
| 90 | 90 | |
| 91 | 91 | # Print usage information |
| 92 | - usage = 'Usage: show_pull_evolution filename parname [file]' | |
| 92 | + usage = 'Usage: show_pull_evolution.py filename parname [file]' | |
| 93 | 93 | if len(sys.argv) < 3: |
| 94 | 94 | print(usage) |
| 95 | 95 | sys.exit() | ... | ... |
examples/show_pull_histogram.py
| ... | ... | @@ -90,7 +90,7 @@ def read_pull(filename, parname): |
| 90 | 90 | if __name__ == '__main__': |
| 91 | 91 | |
| 92 | 92 | # Print usage information |
| 93 | - usage = 'Usage: show_pull_histogram filename parname [bins] [file]' | |
| 93 | + usage = 'Usage: show_pull_histogram.py filename parname [bins] [file]' | |
| 94 | 94 | if len(sys.argv) < 3: |
| 95 | 95 | print(usage) |
| 96 | 96 | sys.exit() | ... | ... |
examples/show_response.py
| 1 | 1 | #! /usr/bin/env python |
| 2 | 2 | # ========================================================================== |
| 3 | -# This script displays one or several CTA response functions. | |
| 3 | +# Shows one or several CTA response functions. | |
| 4 | 4 | # |
| 5 | -# Required 3rd party modules: | |
| 6 | -# - matplotlib | |
| 7 | -# - numpy | |
| 8 | -# | |
| 9 | -# Copyright (C) 2014-2015 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2014-2016 Juergen Knoedlseder | |
| 10 | 6 | # |
| 11 | 7 | # This program is free software: you can redistribute it and/or modify |
| 12 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -25,7 +21,13 @@ |
| 25 | 21 | import sys |
| 26 | 22 | import math |
| 27 | 23 | import gammalib |
| 28 | -import matplotlib.pyplot as plt | |
| 24 | +try: | |
| 25 | + import matplotlib.pyplot as plt | |
| 26 | + plt.figure() | |
| 27 | + plt.close() | |
| 28 | +except: | |
| 29 | + print('This script needs the "matplotlib" module') | |
| 30 | + sys.exit() | |
| 29 | 31 | |
| 30 | 32 | |
| 31 | 33 | # ================================== # |
| ... | ... | @@ -35,11 +37,19 @@ def sigma_lima(Non, Noff, alpha=0.2): |
| 35 | 37 | """ |
| 36 | 38 | Compute Eq. (17) of Li & Ma (1983). |
| 37 | 39 | |
| 38 | - Parameters: | |
| 39 | - Non - Number of on counts | |
| 40 | - Noff - Number of off counts | |
| 41 | - Keywords: | |
| 42 | - alpha - Ratio of on-to-off exposure | |
| 40 | + Parameters | |
| 41 | + ---------- | |
| 42 | + Non : float | |
| 43 | + Number of On counts | |
| 44 | + Noff : float | |
| 45 | + Number of Off counts | |
| 46 | + alpha : float, optional | |
| 47 | + Ratio of on-to-off exposure | |
| 48 | + | |
| 49 | + Returns | |
| 50 | + ------- | |
| 51 | + sigma : float | |
| 52 | + Detection significance in Gaussian sigma | |
| 43 | 53 | """ |
| 44 | 54 | # Compute sensitivity |
| 45 | 55 | alpha1 = alpha + 1.0 |
| ... | ... | @@ -61,11 +71,19 @@ def Non_lima(sigma, Noff, alpha=0.2): |
| 61 | 71 | """ |
| 62 | 72 | Solve Eq. (17) of Li & Ma (1983) for Non |
| 63 | 73 | |
| 64 | - Parameters: | |
| 65 | - sigma - Required significance | |
| 66 | - Noff - Number of off counts | |
| 67 | - Keywords: | |
| 68 | - alpha - Ratio of on-to-off exposure | |
| 74 | + Parameters | |
| 75 | + ---------- | |
| 76 | + sigma : float | |
| 77 | + Required detection significance in Gaussian sigma | |
| 78 | + Noff : float | |
| 79 | + Number of Off counts | |
| 80 | + alpha : float, optional | |
| 81 | + Ratio of on-to-off exposure | |
| 82 | + | |
| 83 | + Returns | |
| 84 | + ------- | |
| 85 | + Non : float | |
| 86 | + Number of On counts | |
| 69 | 87 | """ |
| 70 | 88 | # Set initial guess for Non |
| 71 | 89 | S = 25.0 |
| ... | ... | @@ -89,9 +107,18 @@ def Non_lima(sigma, Noff, alpha=0.2): |
| 89 | 107 | # ======================= # |
| 90 | 108 | # Show one effective area # |
| 91 | 109 | # ======================= # |
| 92 | -def show_one_effective_area(rsp, name, color="r"): | |
| 110 | +def show_one_effective_area(rsp, name, color='r'): | |
| 93 | 111 | """ |
| 94 | - Show one effective area. | |
| 112 | + Show one effective area | |
| 113 | + | |
| 114 | + Parameters | |
| 115 | + ---------- | |
| 116 | + rsp : `~gammalib.GCTAResponse` | |
| 117 | + Response function | |
| 118 | + name : str | |
| 119 | + Name of the response function | |
| 120 | + color : str, optional | |
| 121 | + Color for plot | |
| 95 | 122 | """ |
| 96 | 123 | # Set to figure 1 |
| 97 | 124 | plt.figure(1) |
| ... | ... | @@ -111,11 +138,11 @@ def show_one_effective_area(rsp, name, color="r"): |
| 111 | 138 | plt.loglog(E, aeff, color+'o') |
| 112 | 139 | |
| 113 | 140 | # Set axes |
| 114 | - plt.xlabel("Energy [TeV]") | |
| 115 | - plt.ylabel("Effective area [cm2]") | |
| 141 | + plt.xlabel('Energy (TeV)') | |
| 142 | + plt.ylabel(r'Effective area (cm$^{2}$)') | |
| 116 | 143 | |
| 117 | 144 | # Set legend |
| 118 | - plt.legend(loc="lower right") | |
| 145 | + plt.legend(loc='lower right') | |
| 119 | 146 | |
| 120 | 147 | # Return |
| 121 | 148 | return |
| ... | ... | @@ -124,9 +151,18 @@ def show_one_effective_area(rsp, name, color="r"): |
| 124 | 151 | # ======================== # |
| 125 | 152 | # Show one background rate # |
| 126 | 153 | # ======================== # |
| 127 | -def show_one_background_rate(rsp, name, color="r"): | |
| 154 | +def show_one_background_rate(rsp, name, color='r'): | |
| 128 | 155 | """ |
| 129 | - Show one background rate. | |
| 156 | + Show one background rate | |
| 157 | + | |
| 158 | + Parameters | |
| 159 | + ---------- | |
| 160 | + rsp : `~gammalib.GCTAResponse` | |
| 161 | + Response function | |
| 162 | + name : str | |
| 163 | + Name of the response function | |
| 164 | + color : str, optional | |
| 165 | + Color for plot | |
| 130 | 166 | """ |
| 131 | 167 | # Set to figure 2 |
| 132 | 168 | plt.figure(2) |
| ... | ... | @@ -148,11 +184,11 @@ def show_one_background_rate(rsp, name, color="r"): |
| 148 | 184 | plt.loglog(E, bgrate, color+'o') |
| 149 | 185 | |
| 150 | 186 | # Set axes |
| 151 | - plt.xlabel("Energy [TeV]") | |
| 152 | - plt.ylabel("Background rate [events/s/MeV/sr]") | |
| 187 | + plt.xlabel('Energy (TeV)') | |
| 188 | + plt.ylabel(r'Background rate (counts s$^{-1}$ MeV$^{-1}$ sr$^{-1}$)') | |
| 153 | 189 | |
| 154 | 190 | # Set legend |
| 155 | - plt.legend(loc="lower left") | |
| 191 | + plt.legend(loc='lower left') | |
| 156 | 192 | |
| 157 | 193 | # Return |
| 158 | 194 | return |
| ... | ... | @@ -161,9 +197,25 @@ def show_one_background_rate(rsp, name, color="r"): |
| 161 | 197 | # ==================== # |
| 162 | 198 | # Show one sensitivity # |
| 163 | 199 | # ==================== # |
| 164 | -def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sigma=5.0): | |
| 200 | +def show_one_sensitivity(rsp, name, color='r', duration=180000.0, alpha=0.2, | |
| 201 | + sigma=5.0): | |
| 165 | 202 | """ |
| 166 | - Show one sensitivity. | |
| 203 | + Show one sensitivity | |
| 204 | + | |
| 205 | + Parameters | |
| 206 | + ---------- | |
| 207 | + rsp : `~gammalib.GCTAResponse` | |
| 208 | + Response function | |
| 209 | + name : str | |
| 210 | + Name of the response function | |
| 211 | + color : str, optional | |
| 212 | + Color for plot | |
| 213 | + duration : float, optional | |
| 214 | + Duration of observation (s) | |
| 215 | + alpha : float, optional | |
| 216 | + Ratio of on-to-off exposure | |
| 217 | + sigma : float, optional | |
| 218 | + Required detection significance in Gaussian sigma | |
| 167 | 219 | """ |
| 168 | 220 | # Set constants |
| 169 | 221 | r68_to_sigma = 0.6624305 |
| ... | ... | @@ -194,9 +246,9 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig |
| 194 | 246 | src_counts = 0.05*bgd_counts |
| 195 | 247 | if src_counts < 10: |
| 196 | 248 | src_counts = 10.0 |
| 197 | - emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), "TeV") | |
| 198 | - emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), "TeV") | |
| 199 | - epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), "TeV") | |
| 249 | + emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), 'TeV') | |
| 250 | + emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), 'TeV') | |
| 251 | + epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), 'TeV') | |
| 200 | 252 | plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot) |
| 201 | 253 | conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax) |
| 202 | 254 | 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 |
| 206 | 258 | plt.loglog(E, flux, color+'o') |
| 207 | 259 | |
| 208 | 260 | # Set axes |
| 209 | - plt.xlabel("Energy [TeV]") | |
| 210 | - plt.ylabel("Sensitivity [erg/cm2/s]") | |
| 261 | + plt.xlabel('Energy (TeV)') | |
| 262 | + plt.ylabel(r'Sensitivity (erg cm$^{-2}$ s$^{-1}$)') | |
| 211 | 263 | |
| 212 | 264 | # Set legend |
| 213 | - plt.legend(loc="upper right") | |
| 265 | + plt.legend(loc='upper right') | |
| 214 | 266 | |
| 215 | 267 | # Return |
| 216 | 268 | return |
| ... | ... | @@ -219,9 +271,22 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig |
| 219 | 271 | # ==================== # |
| 220 | 272 | # Show one sensitivity # |
| 221 | 273 | # ==================== # |
| 222 | -def show_one_response(rspname, dbname, name, rootdir=None, color="r"): | |
| 274 | +def show_one_response(rspname, dbname, name, rootdir=None, color='r'): | |
| 223 | 275 | """ |
| 224 | 276 | Show one response. |
| 277 | + | |
| 278 | + Parameters | |
| 279 | + ---------- | |
| 280 | + rspname : str | |
| 281 | + Response name | |
| 282 | + dbname : str | |
| 283 | + Database name | |
| 284 | + name : str | |
| 285 | + Name of the response function | |
| 286 | + rootdir : str, optional | |
| 287 | + Response root directory | |
| 288 | + color : str, optional | |
| 289 | + Color for plot | |
| 225 | 290 | """ |
| 226 | 291 | # Set-up calibration database |
| 227 | 292 | caldb = gammalib.GCaldb() |
| ... | ... | @@ -230,7 +295,7 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): |
| 230 | 295 | if gammalib.dir_exists(dbname): |
| 231 | 296 | caldb.rootdir(dbname) |
| 232 | 297 | else: |
| 233 | - caldb.open("cta", dbname) | |
| 298 | + caldb.open('cta', dbname) | |
| 234 | 299 | |
| 235 | 300 | # Load response function |
| 236 | 301 | rsp = gammalib.GCTAResponseIrf(rspname, caldb) |
| ... | ... | @@ -253,28 +318,35 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): |
| 253 | 318 | # ======================== # |
| 254 | 319 | if __name__ == '__main__': |
| 255 | 320 | |
| 321 | + # Print usage information | |
| 322 | + usage = 'Usage: show_response.py [file]' | |
| 323 | + if len(sys.argv) < 1: | |
| 324 | + print(usage) | |
| 325 | + sys.exit() | |
| 326 | + | |
| 327 | + # Check if plotting in file is requested | |
| 328 | + plotfile = '' | |
| 329 | + if len(sys.argv) == 2: | |
| 330 | + plotfile = sys.argv[1] | |
| 331 | + | |
| 256 | 332 | # Create figures |
| 257 | 333 | plt.figure(1) |
| 258 | - plt.title("Effective area") | |
| 334 | + plt.title('Effective area') | |
| 259 | 335 | plt.figure(2) |
| 260 | - plt.title("Background rate") | |
| 336 | + plt.title('Background rate') | |
| 261 | 337 | plt.figure(3) |
| 262 | - plt.title("Sensitivity") | |
| 338 | + plt.title('Sensitivity') | |
| 263 | 339 | |
| 264 | 340 | # Set response dictionary |
| 265 | 341 | rsps = [ |
| 266 | - {'dbname': "dummy", | |
| 267 | - 'rspname': "cta_dummy_irf", | |
| 268 | - 'name': "Prod1 (Array E, MPIK)", | |
| 269 | - 'color': "b"}, | |
| 270 | - {'dbname': "e", | |
| 271 | - 'rspname': "IFAE20120510_50h", | |
| 272 | - 'name': "Prod1 (Array E, IFAE)", | |
| 273 | - 'color': "g"}, | |
| 274 | - {'dbname': "aar", | |
| 275 | - 'rspname': "DESY20140105_50h", | |
| 276 | - 'name': "Prod2 (Aar, DESY)", | |
| 277 | - 'color': "r"} | |
| 342 | + {'dbname': 'irf', | |
| 343 | + 'rspname': 'cta_dummy_irf.dat', | |
| 344 | + 'name': 'Prod1 (Array E, MPIK)', | |
| 345 | + 'color': 'b'}, | |
| 346 | + {'dbname': 'irf', | |
| 347 | + 'rspname': 'irf_file.fits', | |
| 348 | + 'name': 'Prod2 (Aar, IFAE)', | |
| 349 | + 'color': 'r'} | |
| 278 | 350 | ] |
| 279 | 351 | |
| 280 | 352 | # Loop over all responses |
| ... | ... | @@ -293,5 +365,8 @@ if __name__ == '__main__': |
| 293 | 365 | # Show response |
| 294 | 366 | show_one_response(rspname, dbname, name, rootdir=rootdir, color=color) |
| 295 | 367 | |
| 296 | - # Show plot | |
| 297 | - plt.show() | |
| 368 | + # Show figure | |
| 369 | + if len(plotfile) > 0: | |
| 370 | + plt.savefig(plotfile) | |
| 371 | + else: | |
| 372 | + plt.show() | ... | ... |
examples/show_sensitivity.py
| 1 | 1 | #! /usr/bin/env python |
| 2 | 2 | # ========================================================================== |
| 3 | -# This script displays a differential sensitivity plot generated by cssens. | |
| 3 | +# Display differential sensitivity plot generated by cssens | |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2015 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2015-2016 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 |
| ... | ... | @@ -23,8 +23,11 @@ import math |
| 23 | 23 | import gammalib |
| 24 | 24 | try: |
| 25 | 25 | import matplotlib.pyplot as plt |
| 26 | + plt.figure() | |
| 27 | + plt.close() | |
| 26 | 28 | except: |
| 27 | - sys.exit("This script needs matplotlib") | |
| 29 | + print('This script needs the "matplotlib" module') | |
| 30 | + sys.exit() | |
| 28 | 31 | |
| 29 | 32 | |
| 30 | 33 | # ===================== # |
| ... | ... | @@ -32,7 +35,17 @@ except: |
| 32 | 35 | # ===================== # |
| 33 | 36 | def read_sensitivity(filename): |
| 34 | 37 | """ |
| 35 | - Read sensitivity information from CSV file. | |
| 38 | + Read sensitivity information from CSV file | |
| 39 | + | |
| 40 | + Parameters | |
| 41 | + ---------- | |
| 42 | + filename : str | |
| 43 | + Name of CSV file | |
| 44 | + | |
| 45 | + Returns | |
| 46 | + ------- | |
| 47 | + sensitivity : dict | |
| 48 | + Dictionary with sensitivity information | |
| 36 | 49 | """ |
| 37 | 50 | # Read filename |
| 38 | 51 | csv = gammalib.GCsv(filename,',') |
| ... | ... | @@ -40,39 +53,39 @@ def read_sensitivity(filename): |
| 40 | 53 | # Create dictionnary |
| 41 | 54 | sensitivity = {} |
| 42 | 55 | for column in range(csv.ncols()): |
| 43 | - name = csv[0,column] | |
| 56 | + name = csv[0,column].rstrip('\r') | |
| 44 | 57 | values = [] |
| 45 | 58 | for row in range(csv.nrows()-1): |
| 46 | 59 | values.append(float(csv[row+1,column])) |
| 47 | 60 | sensitivity[name] = values |
| 48 | 61 | |
| 49 | 62 | # Check where file contains differential or integral sensitivity |
| 50 | - mode = "Integral" | |
| 51 | - if sensitivity.has_key("emax"): | |
| 63 | + mode = 'Integral' | |
| 64 | + if sensitivity.has_key('emax'): | |
| 52 | 65 | emax_ref = -1.0 |
| 53 | - for value in sensitivity["emax"]: | |
| 66 | + for value in sensitivity['emax']: | |
| 54 | 67 | if emax_ref < 0.0: |
| 55 | 68 | emax_ref = value |
| 56 | 69 | elif emax_ref != value: |
| 57 | - mode = "Differential" | |
| 70 | + mode = 'Differential' | |
| 58 | 71 | break |
| 59 | 72 | |
| 60 | 73 | # Add mode |
| 61 | - sensitivity["mode"] = mode | |
| 74 | + sensitivity['mode'] = mode | |
| 62 | 75 | |
| 63 | 76 | # Add linear energy values |
| 64 | - if mode == "Differential": | |
| 65 | - if sensitivity.has_key("loge"): | |
| 77 | + if mode == 'Differential': | |
| 78 | + if sensitivity.has_key('loge'): | |
| 66 | 79 | name = 'energy' |
| 67 | 80 | values = [] |
| 68 | - for value in sensitivity["loge"]: | |
| 81 | + for value in sensitivity['loge']: | |
| 69 | 82 | values.append(math.pow(10.0, value)) |
| 70 | 83 | sensitivity[name] = values |
| 71 | 84 | else: |
| 72 | - if sensitivity.has_key("emin"): | |
| 85 | + if sensitivity.has_key('emin'): | |
| 73 | 86 | name = 'energy' |
| 74 | 87 | values = [] |
| 75 | - for value in sensitivity["emin"]: | |
| 88 | + for value in sensitivity['emin']: | |
| 76 | 89 | values.append(value) |
| 77 | 90 | sensitivity[name] = values |
| 78 | 91 | |
| ... | ... | @@ -83,9 +96,18 @@ def read_sensitivity(filename): |
| 83 | 96 | # ===================== # |
| 84 | 97 | # Show sensitivity data # |
| 85 | 98 | # ===================== # |
| 86 | -def show_sensitivity(sensitivity, filename): | |
| 99 | +def show_sensitivity(sensitivity, filename, plotfile=''): | |
| 87 | 100 | """ |
| 88 | - Plot sensitivity. | |
| 101 | + Plot sensitivity data | |
| 102 | + | |
| 103 | + Parameters | |
| 104 | + ---------- | |
| 105 | + sensitivity : dict | |
| 106 | + Dictionary with sensitivity information | |
| 107 | + filename : str | |
| 108 | + Name of CSV file | |
| 109 | + plotfile : str, optional | |
| 110 | + Plot file name | |
| 89 | 111 | """ |
| 90 | 112 | # Build title |
| 91 | 113 | title = sensitivity['mode'] + ' sensitivity ('+filename+')' |
| ... | ... | @@ -102,14 +124,17 @@ def show_sensitivity(sensitivity, filename): |
| 102 | 124 | plt.plot(sensitivity['energy'], sensitivity['sensitivity'], 'ro-') |
| 103 | 125 | |
| 104 | 126 | # Set labels |
| 105 | - plt.xlabel("Energy (TeV)") | |
| 106 | - if sensitivity['mode'] == "Differential": | |
| 107 | - plt.ylabel(r"E$\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)") | |
| 127 | + plt.xlabel('Energy (TeV)') | |
| 128 | + if sensitivity['mode'] == 'Differential': | |
| 129 | + plt.ylabel(r'E $\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)') | |
| 108 | 130 | else: |
| 109 | - plt.ylabel(r"E$\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)") | |
| 131 | + plt.ylabel(r'E $\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)') | |
| 110 | 132 | |
| 111 | - # Show plot | |
| 112 | - plt.show() | |
| 133 | + # Show figure | |
| 134 | + if len(plotfile) > 0: | |
| 135 | + plt.savefig(plotfile) | |
| 136 | + else: | |
| 137 | + plt.show() | |
| 113 | 138 | |
| 114 | 139 | # Return |
| 115 | 140 | return |
| ... | ... | @@ -121,9 +146,9 @@ def show_sensitivity(sensitivity, filename): |
| 121 | 146 | if __name__ == '__main__': |
| 122 | 147 | |
| 123 | 148 | # Print usage information |
| 124 | - usage = "Usage: show_sensitivity [filename]" | |
| 149 | + usage = 'Usage: show_sensitivity.py [filename] [file]' | |
| 125 | 150 | if len(sys.argv) < 1: |
| 126 | - sys.stdout.write(usage+"\n") | |
| 151 | + print(usage) | |
| 127 | 152 | sys.exit() |
| 128 | 153 | |
| 129 | 154 | # Extract parameters |
| ... | ... | @@ -131,9 +156,12 @@ if __name__ == '__main__': |
| 131 | 156 | filename = 'sensitivity.dat' |
| 132 | 157 | else: |
| 133 | 158 | filename = sys.argv[1] |
| 159 | + plotfile = '' | |
| 160 | + if len(sys.argv) == 3: | |
| 161 | + plotfile = sys.argv[2] | |
| 134 | 162 | |
| 135 | 163 | # Read sensitivity data |
| 136 | 164 | sensitivity = read_sensitivity(filename) |
| 137 | 165 | |
| 138 | 166 | # Show sensitivity data |
| 139 | - show_sensitivity(sensitivity, filename) | |
| 167 | + show_sensitivity(sensitivity, filename, plotfile=plotfile) | ... | ... |
examples/show_spectrum.py
| 1 | 1 | #! /usr/bin/env python |
| 2 | 2 | # ========================================================================== |
| 3 | -# This script shows a spectrum created with csspec | |
| 3 | +# Display spectrum generated by csspec | |
| 4 | 4 | # |
| 5 | 5 | # Copyright (C) 2015-2016 Juergen Knoedlseder |
| 6 | 6 | # |
| ... | ... | @@ -18,35 +18,41 @@ |
| 18 | 18 | # along with this program. If not, see <http://www.gnu.org/licenses/>. |
| 19 | 19 | # |
| 20 | 20 | # ========================================================================== |
| 21 | - | |
| 22 | -import gammalib | |
| 23 | 21 | import sys |
| 22 | +import gammalib | |
| 24 | 23 | try: |
| 25 | 24 | import matplotlib.pyplot as plt |
| 25 | + plt.figure() | |
| 26 | + plt.close() | |
| 26 | 27 | except: |
| 27 | - sys.exit("This script needs matplotlib") | |
| 28 | + print('This script needs the "matplotlib" module') | |
| 29 | + sys.exit() | |
| 28 | 30 | |
| 29 | 31 | |
| 30 | -# ============ # | |
| 31 | -# Script entry # | |
| 32 | -# ============ # | |
| 33 | -if __name__ == "__main__": | |
| 34 | - | |
| 35 | - # Check for given spectrum file | |
| 36 | - if not len(sys.argv) == 2: | |
| 37 | - sys.exit("Usage: show_spectrum.py spectrum.fits") | |
| 32 | +# ============= # | |
| 33 | +# Show spectrum # | |
| 34 | +# ============= # | |
| 35 | +def show_spectrum(filename, plotfile=''): | |
| 36 | + """ | |
| 37 | + Plot sensitivity data | |
| 38 | 38 | |
| 39 | + Parameters | |
| 40 | + ---------- | |
| 41 | + filename : str | |
| 42 | + Name of spectrum FITS file | |
| 43 | + plotfile : str, optional | |
| 44 | + Plot file name | |
| 45 | + """ | |
| 39 | 46 | # Read spectrum file |
| 40 | - filename = sys.argv[1] | |
| 41 | 47 | fits = gammalib.GFits(filename) |
| 42 | 48 | table = fits.table(1) |
| 43 | - c_energy = table["Energy"] | |
| 44 | - c_ed = table["ed_Energy"] | |
| 45 | - c_eu = table["eu_Energy"] | |
| 46 | - c_flux = table["Flux"] | |
| 47 | - c_eflux = table["e_Flux"] | |
| 48 | - c_ts = table["TS"] | |
| 49 | - c_upper = table["UpperLimit"] | |
| 49 | + c_energy = table['Energy'] | |
| 50 | + c_ed = table['ed_Energy'] | |
| 51 | + c_eu = table['eu_Energy'] | |
| 52 | + c_flux = table['Flux'] | |
| 53 | + c_eflux = table['e_Flux'] | |
| 54 | + c_ts = table['TS'] | |
| 55 | + c_upper = table['UpperLimit'] | |
| 50 | 56 | |
| 51 | 57 | # Initialise arrays to be filled |
| 52 | 58 | energies = [] |
| ... | ... | @@ -63,25 +69,22 @@ if __name__ == "__main__": |
| 63 | 69 | nrows = table.nrows() |
| 64 | 70 | for row in range(nrows): |
| 65 | 71 | |
| 66 | - # Get TS | |
| 72 | + # Get Test Statistic, flux and flux error | |
| 67 | 73 | ts = c_ts.real(row) |
| 68 | 74 | flx = c_flux.real(row) |
| 69 | 75 | e_flx = c_eflux.real(row) |
| 70 | 76 | |
| 71 | - # Switch | |
| 77 | + # If Test Statistic is larger than 9 and flux error is smaller than | |
| 78 | + # flux then append flux plots ... | |
| 72 | 79 | if ts > 9.0 and e_flx < flx: |
| 73 | - | |
| 74 | - # Add information | |
| 75 | 80 | energies.append(c_energy.real(row)) |
| 76 | 81 | flux.append(c_flux.real(row)) |
| 77 | 82 | ed_engs.append(c_ed.real(row)) |
| 78 | 83 | eu_engs.append(c_eu.real(row)) |
| 79 | 84 | e_flux.append(c_eflux.real(row)) |
| 80 | 85 | |
| 81 | - # | |
| 86 | + # ... otherwise append upper limit | |
| 82 | 87 | else: |
| 83 | - | |
| 84 | - # Add information | |
| 85 | 88 | ul_energies.append(c_energy.real(row)) |
| 86 | 89 | ul_flux.append(c_upper.real(row)) |
| 87 | 90 | ul_ed_engs.append(c_ed.real(row)) |
| ... | ... | @@ -91,8 +94,39 @@ if __name__ == "__main__": |
| 91 | 94 | plt.figure() |
| 92 | 95 | plt.loglog() |
| 93 | 96 | plt.grid() |
| 94 | - plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], fmt='ro') | |
| 95 | - plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], yerr=1.0e-11, uplims=True, fmt='ro') | |
| 96 | - plt.xlabel("Energy (TeV)") | |
| 97 | - plt.ylabel(r"E dN/dE (erg cm$^{-2}$ s$^{-1}$)") | |
| 98 | - plt.show() | |
| 97 | + plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], | |
| 98 | + fmt='ro') | |
| 99 | + plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], | |
| 100 | + yerr=1.0e-11, uplims=True, fmt='ro') | |
| 101 | + plt.xlabel('Energy (TeV)') | |
| 102 | + plt.ylabel(r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)') | |
| 103 | + | |
| 104 | + # Show figure | |
| 105 | + if len(plotfile) > 0: | |
| 106 | + plt.savefig(plotfile) | |
| 107 | + else: | |
| 108 | + plt.show() | |
| 109 | + | |
| 110 | + # Return | |
| 111 | + return | |
| 112 | + | |
| 113 | + | |
| 114 | +# ======================== # | |
| 115 | +# Main routine entry point # | |
| 116 | +# ======================== # | |
| 117 | +if __name__ == '__main__': | |
| 118 | + | |
| 119 | + # Print usage information | |
| 120 | + usage = 'Usage: show_spectrum.py filename [file]' | |
| 121 | + if len(sys.argv) < 2: | |
| 122 | + print(usage) | |
| 123 | + sys.exit() | |
| 124 | + | |
| 125 | + # Extract parameters | |
| 126 | + filename = sys.argv[1] | |
| 127 | + plotfile = '' | |
| 128 | + if len(sys.argv) == 3: | |
| 129 | + plotfile = sys.argv[2] | |
| 130 | + | |
| 131 | + # Show spectrum | |
| 132 | + show_spectrum(filename, plotfile=plotfile) | ... | ... |
examples/show_ts_distribution.py
| 1 | 1 | #! /usr/bin/env python |
| 2 | 2 | # ========================================================================== |
| 3 | -# This script displays the TS distribution values generated using cstsdist. | |
| 3 | +# Display TS distribution generated by cstsdist | |
| 4 | 4 | # |
| 5 | -# Required 3rd party modules: | |
| 6 | -# - matplotlib | |
| 7 | -# - numpy | |
| 8 | -# | |
| 9 | -# Copyright (C) 2011 Jurgen Knodlseder | |
| 5 | +# Copyright (C) 2011-2016 Jurgen Knodlseder | |
| 10 | 6 | # |
| 11 | 7 | # This program is free software: you can redistribute it and/or modify |
| 12 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -22,20 +18,42 @@ |
| 22 | 18 | # along with this program. If not, see <http://www.gnu.org/licenses/>. |
| 23 | 19 | # |
| 24 | 20 | # ========================================================================== |
| 25 | -import matplotlib.mlab as mlab | |
| 26 | -import matplotlib.pyplot as plt | |
| 27 | -import numpy as np | |
| 28 | 21 | import sys |
| 29 | 22 | import csv |
| 30 | 23 | import math |
| 24 | +try: | |
| 25 | + import matplotlib.pyplot as plt | |
| 26 | + import matplotlib.mlab as mlab | |
| 27 | + plt.figure() | |
| 28 | + plt.close() | |
| 29 | +except: | |
| 30 | + print('This script needs the "matplotlib" module') | |
| 31 | + sys.exit() | |
| 32 | +try: | |
| 33 | + import numpy as np | |
| 34 | +except: | |
| 35 | + print('This script needs the "numpy" module') | |
| 36 | + sys.exit() | |
| 31 | 37 | |
| 32 | 38 | |
| 33 | 39 | # ==================== # |
| 34 | 40 | # Read TS distribution # |
| 35 | 41 | # ==================== # |
| 36 | -def read_ts(filename, tsname="TS"): | |
| 42 | +def read_ts(filename, tsname='TS'): | |
| 37 | 43 | """ |
| 38 | - Read TS distribution values from CSV file. | |
| 44 | + Read TS distribution values from CSV file | |
| 45 | + | |
| 46 | + Parameters | |
| 47 | + ---------- | |
| 48 | + filename : str | |
| 49 | + Name of CSV file | |
| 50 | + tsname : str, optional | |
| 51 | + Column name of TS | |
| 52 | + | |
| 53 | + Returns | |
| 54 | + ------- | |
| 55 | + a : `numpy.array` | |
| 56 | + Array of TS values | |
| 39 | 57 | """ |
| 40 | 58 | # Initialise list |
| 41 | 59 | values = [] |
| ... | ... | @@ -53,9 +71,9 @@ def read_ts(filename, tsname="TS"): |
| 53 | 71 | try: |
| 54 | 72 | index = row.index(tsname) |
| 55 | 73 | except: |
| 56 | - sys.stdout.write('ERROR: Column "TS" not found in file\n') | |
| 57 | - sys.stdout.write(row+"\n") | |
| 58 | - raise NameError("TS") | |
| 74 | + print('ERROR: Column "'+tsname+'" not found in file') | |
| 75 | + print(row) | |
| 76 | + raise NameError(tsname) | |
| 59 | 77 | |
| 60 | 78 | # Handle data rows |
| 61 | 79 | else: |
| ... | ... | @@ -76,7 +94,17 @@ def read_ts(filename, tsname="TS"): |
| 76 | 94 | # =========== # |
| 77 | 95 | def erf(x): |
| 78 | 96 | """ |
| 79 | - Compute error function. | |
| 97 | + Compute error function | |
| 98 | + | |
| 99 | + Parameters | |
| 100 | + ---------- | |
| 101 | + x : float | |
| 102 | + Value | |
| 103 | + | |
| 104 | + Returns | |
| 105 | + ------- | |
| 106 | + erf(-x) : float | |
| 107 | + Error function | |
| 80 | 108 | """ |
| 81 | 109 | # save the sign of x |
| 82 | 110 | sign = 1 |
| ... | ... | @@ -103,9 +131,20 @@ def erf(x): |
| 103 | 131 | # =============================== # |
| 104 | 132 | # Display cumulative distribution # |
| 105 | 133 | # =============================== # |
| 106 | -def dist_cdf(values, nbins, title): | |
| 134 | +def dist_cdf(values, nbins, title, expect='', plotfile=''): | |
| 107 | 135 | """ |
| 108 | - Displays cumulative TS distribution and overplots expectation. | |
| 136 | + Display cumulative TS distribution and overplots expectation | |
| 137 | + | |
| 138 | + Parameters | |
| 139 | + ---------- | |
| 140 | + values : float | |
| 141 | + TS values | |
| 142 | + nbins : int | |
| 143 | + Number of bins | |
| 144 | + title : str | |
| 145 | + Plot title | |
| 146 | + plotfile : str, optional | |
| 147 | + Plot file name | |
| 109 | 148 | """ |
| 110 | 149 | # Set range and adapt number of bins. We make sure that we have a |
| 111 | 150 | # histogram bin centred on 0 that will capture the TS=0 part. |
| ... | ... | @@ -142,19 +181,23 @@ def dist_cdf(values, nbins, title): |
| 142 | 181 | hist_expect.append(norm) |
| 143 | 182 | |
| 144 | 183 | # Show expected distribution (half) |
| 145 | - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations") | |
| 146 | - l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation") | |
| 147 | - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) | |
| 184 | + plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label='simulations') | |
| 185 | + if expect == 'show': | |
| 186 | + plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label='expectation') | |
| 187 | + plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) | |
| 148 | 188 | |
| 149 | 189 | # Set plot |
| 150 | 190 | plt.xlabel('TS') |
| 151 | 191 | plt.ylabel('Number of trials > TS') |
| 152 | 192 | plt.title(title) |
| 153 | 193 | plt.grid(True) |
| 154 | - plt.legend(loc="upper right") | |
| 194 | + plt.legend(loc='upper right') | |
| 155 | 195 | |
| 156 | - # Show histogram | |
| 157 | - plt.show() | |
| 196 | + # Show figure | |
| 197 | + if len(plotfile) > 0: | |
| 198 | + plt.savefig(plotfile) | |
| 199 | + else: | |
| 200 | + plt.show() | |
| 158 | 201 | |
| 159 | 202 | # Return |
| 160 | 203 | return |
| ... | ... | @@ -163,10 +206,20 @@ def dist_cdf(values, nbins, title): |
| 163 | 206 | # ======================================== # |
| 164 | 207 | # Display probability density distribution # |
| 165 | 208 | # ======================================== # |
| 166 | -def dist_pdf(values, nbins, title): | |
| 209 | +def dist_pdf(values, nbins, title, expect='', plotfile=''): | |
| 167 | 210 | """ |
| 168 | - Display probability density TS distribution and overplots | |
| 169 | - expectation. | |
| 211 | + Display probability density TS distribution and overplots expectation | |
| 212 | + | |
| 213 | + Parameters | |
| 214 | + ---------- | |
| 215 | + values : float | |
| 216 | + TS values | |
| 217 | + nbins : int | |
| 218 | + Number of bins | |
| 219 | + title : str | |
| 220 | + Plot title | |
| 221 | + plotfile : str, optional | |
| 222 | + Plot file name | |
| 170 | 223 | """ |
| 171 | 224 | # Set range and adapt number of bins. We make sure that we have a |
| 172 | 225 | # histogram bin centred on 0 that will capture the TS=0 part. |
| ... | ... | @@ -180,9 +233,9 @@ def dist_pdf(values, nbins, title): |
| 180 | 233 | # Create histogram |
| 181 | 234 | n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], |
| 182 | 235 | align='mid', facecolor='green', |
| 183 | - log=True, label="simulations") | |
| 236 | + log=True, label='simulations') | |
| 184 | 237 | |
| 185 | - # Create expected distribution (onyl for positive TS). We compute the | |
| 238 | + # Create expected distribution (only for positive TS). We compute the | |
| 186 | 239 | # full and half of the distribution, as for positively constrained |
| 187 | 240 | # quantities half of the values have TS=0 |
| 188 | 241 | x = [] |
| ... | ... | @@ -202,17 +255,24 @@ def dist_pdf(values, nbins, title): |
| 202 | 255 | y_half.append(0.5*y) |
| 203 | 256 | |
| 204 | 257 | # Show expected distribution (half) |
| 205 | - l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation") | |
| 258 | + if expect == 'show': | |
| 259 | + plt.semilogy(x, y_half, 'ro', linewidth=2, label='expectation') | |
| 206 | 260 | |
| 207 | 261 | # Set plot |
| 208 | 262 | plt.xlabel('TS') |
| 209 | 263 | plt.ylabel('Number of trials') |
| 210 | 264 | plt.title(title) |
| 211 | 265 | plt.grid(True) |
| 212 | - plt.legend(loc="upper right") | |
| 266 | + plt.legend(loc='upper right') | |
| 213 | 267 | |
| 214 | - # Show histogram | |
| 215 | - plt.show() | |
| 268 | + # Show figure | |
| 269 | + if len(plotfile) > 0: | |
| 270 | + plt.savefig(plotfile) | |
| 271 | + else: | |
| 272 | + plt.show() | |
| 273 | + | |
| 274 | + # Return | |
| 275 | + return | |
| 216 | 276 | |
| 217 | 277 | |
| 218 | 278 | # ======================== # |
| ... | ... | @@ -221,26 +281,22 @@ def dist_pdf(values, nbins, title): |
| 221 | 281 | if __name__ == '__main__': |
| 222 | 282 | |
| 223 | 283 | # Print usage information |
| 224 | - usage = "Usage: show_ts_distribution filename" \ | |
| 225 | - " [-n bins] [-c column] [-t title] [-p plot]" | |
| 284 | + usage = 'Usage: show_ts_distribution.py filename' \ | |
| 285 | + ' [-n bins] [-c column] [-t title] [-p plot] [-f file]' | |
| 226 | 286 | if len(sys.argv) < 2: |
| 227 | - sys.stdout.write(usage+"\n") | |
| 287 | + print(usage) | |
| 228 | 288 | sys.exit() |
| 229 | 289 | |
| 230 | 290 | # Extract parameters |
| 231 | 291 | filename = sys.argv[1] |
| 232 | 292 | |
| 233 | - # Set default parameters | |
| 234 | - nbins = 30 | |
| 235 | - tsname = "TS" | |
| 236 | - title = "TS distribution" | |
| 237 | - plot = "pdf" | |
| 238 | - | |
| 239 | 293 | # Parameter dictionnary |
| 240 | - pars = [{'option': '-n', 'value': nbins}, | |
| 241 | - {'option': '-c', 'value': tsname}, | |
| 242 | - {'option': '-t', 'value': title}, | |
| 243 | - {'option': '-p', 'value': plot}] | |
| 294 | + pars = [{'option': '-n', 'value': 30}, | |
| 295 | + {'option': '-c', 'value': 'TS'}, | |
| 296 | + {'option': '-t', 'value': 'TS distribution'}, | |
| 297 | + {'option': '-p', 'value': 'pdf'}, | |
| 298 | + {'option': '-e', 'value': ''}, | |
| 299 | + {'option': '-f', 'value': ''}] | |
| 244 | 300 | |
| 245 | 301 | # Gather parameters from command line |
| 246 | 302 | i = 2 |
| ... | ... | @@ -254,27 +310,29 @@ if __name__ == '__main__': |
| 254 | 310 | try: |
| 255 | 311 | par['value'] = sys.argv[i] |
| 256 | 312 | except: |
| 257 | - sys.stdout.write(usage+"\n") | |
| 313 | + print(usage) | |
| 258 | 314 | sys.exit() |
| 259 | 315 | else: |
| 260 | - sys.stdout.write(usage+"\n") | |
| 316 | + print(usage) | |
| 261 | 317 | sys.exit() |
| 262 | 318 | |
| 263 | 319 | # Next item |
| 264 | 320 | i += 1 |
| 265 | 321 | |
| 266 | 322 | # Recover parameters |
| 267 | - nbins = int(pars[0]['value']) | |
| 268 | - tsname = pars[1]['value'] | |
| 269 | - title = pars[2]['value'] | |
| 270 | - plot = pars[3]['value'] | |
| 323 | + nbins = int(pars[0]['value']) | |
| 324 | + tsname = pars[1]['value'] | |
| 325 | + title = pars[2]['value'] | |
| 326 | + plot = pars[3]['value'] | |
| 327 | + expect = pars[4]['value'] | |
| 328 | + plotfile = pars[5]['value'] | |
| 271 | 329 | |
| 272 | 330 | # Read values from CSV file |
| 273 | 331 | values = read_ts(filename, tsname=tsname) |
| 274 | - sys.stdout.write(len(values), "values read.\n") | |
| 332 | + print('%d values read.' % len(values)) | |
| 275 | 333 | |
| 276 | 334 | # Show histogram |
| 277 | - if plot == "pdf": | |
| 278 | - dist_pdf(values, nbins, title) | |
| 335 | + if plot == 'pdf': | |
| 336 | + dist_pdf(values, nbins, title, expect=expect, plotfile=plotfile) | |
| 279 | 337 | else: |
| 280 | - dist_cdf(values, nbins, title) | |
| 338 | + dist_cdf(values, nbins, title, expect=expect, plotfile=plotfile) | ... | ... |
test/README.md
test/data/sensitivity.dat
0 → 100644
| 1 | +loge,emin,emax,crab_flux,photon_flux,energy_flux,sensitivity | |
| 2 | +-1.51402150412,0.02,0.0468734582318,0.144989163852,6.60543533444e-10,3.05652923293e-11,3.56378721075e-11 | |
| 3 | +-1.14412450369,0.0468734582318,0.109856054331,0.0561134683318,7.24745374549e-11,7.85976500012e-12,9.16416230711e-12 | |
| 4 | +-0.774227503252,0.109856054331,0.257466658709,0.0315037806261,1.15354211075e-11,2.93193651747e-12,3.41851723554e-12 | |
| 5 | +-0.404330502818,0.257466658709,0.603417633655,0.0244879762385,2.54200331885e-12,1.51423828719e-12,1.76553948308e-12 | |
| 6 | +-0.0344335023848,0.603417633655,1.41421356237,0.022256802649,6.54995675049e-13,9.14436478628e-13,1.0661952755e-12 | |
| 7 | +0.335463498049,1.41421356237,3.31445401734,0.0250217133244,2.08758967789e-13,6.83057479879e-13,7.96416891675e-13 | |
| 8 | +0.705360498482,3.31445401734,7.76799609716,0.042333217441,1.00129446033e-13,7.67840339358e-13,8.95270214276e-13 | |
| 9 | +1.07525749892,7.76799609716,18.2056420303,0.0929194773185,6.23074149979e-14,1.11981387734e-12,1.30565686449e-12 | |
| 10 | +1.44515449935,18.2056420303,42.6680700645,0.14972962577,2.8463818009e-14,1.19893673197e-12,1.39791085453e-12 | |
| 11 | +1.81505149978,42.6680700645,100.0,0.382388032371,2.06082989569e-14,2.03442769003e-12,2.37205890421e-12 | ... | ... |
test/data/spectrum.fits
0 → 100644
No preview for this file type
test/data/ts.dat
0 → 100644
| 1 | +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 | |
| 2 | +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 | |
| 3 | +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 | |
| 4 | +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 | |
| 5 | +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 | |
| 6 | +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 | |
| 7 | +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 | |
| 8 | +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 | |
| 9 | +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 | |
| 10 | +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 | |
| 11 | +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 | |
| 12 | +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 | |
| 13 | +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 | |
| 14 | +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 | |
| 15 | +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 | |
| 16 | +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 | |
| 17 | +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 | |
| 18 | +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 | |
| 19 | +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 | |
| 20 | +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 | |
| 21 | +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 | |
| 22 | +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 | |
| 23 | +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 | |
| 24 | +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 | |
| 25 | +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 | |
| 26 | +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 | |
| 27 | +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 | |
| 28 | +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 | |
| 29 | +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 | |
| 30 | +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 | |
| 31 | +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 | |
| 32 | +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 | |
| 33 | +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 | |
| 34 | +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 | |
| 35 | +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 | |
| 36 | +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 | |
| 37 | +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 | |
| 38 | +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 | |
| 39 | +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 | |
| 40 | +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 | |
| 41 | +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 | |
| 42 | +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 | |
| 43 | +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 | |
| 44 | +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 | |
| 45 | +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 | |
| 46 | +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 | |
| 47 | +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 | |
| 48 | +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 | |
| 49 | +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 | |
| 50 | +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 | |
| 51 | +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 | |
| 52 | +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 | |
| 53 | +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 | |
| 54 | +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 | |
| 55 | +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 | |
| 56 | +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 | |
| 57 | +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 | |
| 58 | +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 | |
| 59 | +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 | |
| 60 | +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 | |
| 61 | +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 | |
| 62 | +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 | |
| 63 | +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 | |
| 64 | +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 | |
| 65 | +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 | |
| 66 | +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 | |
| 67 | +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 | |
| 68 | +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 | |
| 69 | +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 | |
| 70 | +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 | |
| 71 | +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 | |
| 72 | +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 | |
| 73 | +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 | |
| 74 | +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 | |
| 75 | +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 | |
| 76 | +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 | |
| 77 | +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 | |
| 78 | +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 | |
| 79 | +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 | |
| 80 | +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 | |
| 81 | +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 | |
| 82 | +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 | |
| 83 | +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 | |
| 84 | +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 | |
| 85 | +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 | |
| 86 | +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 | |
| 87 | +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 | |
| 88 | +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 | |
| 89 | +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 | |
| 90 | +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 | |
| 91 | +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 | |
| 92 | +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 | |
| 93 | +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 | |
| 94 | +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 | |
| 95 | +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 | |
| 96 | +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 | |
| 97 | +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 | |
| 98 | +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 | |
| 99 | +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 | |
| 100 | +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 | |
| 101 | +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 | ... | ... |
test/dev/README.md
0 → 100644
| 1 | +Developer test scripts | |
| 2 | +====================== | |
| 3 | + | |
| 4 | +This folder contains scripts that are used for various testing purposes | |
| 5 | +during code development. None of the scripts actually needs to work, they | |
| 6 | +may actually be outdated, but could still be useful for further testing | |
| 7 | +if needed. | |
| 8 | + | |
| 9 | +`benchmark_cta_analysis.py` | |
| 10 | +> This script peforms a benchmark of the various CTA analysis types | |
| 11 | +> (unbinned, binned, stacked), and if matplotlib is installed, creates | |
| 12 | +> a plot of the benchmark results. | |
| 13 | + | |
| 14 | +`check_models.py` | |
| 15 | +> Checks a variety of models using an unbinned analysis. | |
| 16 | + | |
| 17 | +`make_pull_at_sensitivity_limit.py` | |
| 18 | +> This script creates pull distributions as function of energy | |
| 19 | +> for a test source that is at the 5 sigma sensitivity limit for | |
| 20 | +> 50h of effective observing time. It takes as input file the | |
| 21 | +> output file generated using cspull. | |
| 22 | +> If the Python processing module is installed the script can make | |
| 23 | +> use of multiple cores and runs in parallel. | |
| 24 | + | |
| 25 | +`make_survey.py` | |
| 26 | +> This script simulates a CTA galactic plane survey. | |
| 27 | + | |
| 28 | +`simulate_events.py` | |
| 29 | +> Simulates events and compares the simulation to the model. | |
| 30 | + | ... | ... |
test/test_examples.py
| ... | ... | @@ -104,6 +104,10 @@ class Test(gammalib.GPythonTestSuite): |
| 104 | 104 | self.append(self.test_show_pha, 'Test show_pha.py') |
| 105 | 105 | self.append(self.test_show_pull_evolution, 'Test show_pull_evolution.py') |
| 106 | 106 | self.append(self.test_show_pull_histogram, 'Test show_pull_histogram.py') |
| 107 | + self.append(self.test_show_response, 'Test test_show_response.py') | |
| 108 | + self.append(self.test_show_sensitivity, 'Test test_show_sensitivity.py') | |
| 109 | + self.append(self.test_show_spectrum, 'Test test_show_spectrum.py') | |
| 110 | + self.append(self.test_show_ts_distribution, 'Test show_ts_distribution.py') | |
| 107 | 111 | |
| 108 | 112 | # Return |
| 109 | 113 | return |
| ... | ... | @@ -253,7 +257,8 @@ class Test(gammalib.GPythonTestSuite): |
| 253 | 257 | """ |
| 254 | 258 | # Execute script |
| 255 | 259 | self._execute_python('show_butterfly', |
| 256 | - args='data/butterfly.txt example_butterfly.eps') | |
| 260 | + args='data/butterfly.txt ' | |
| 261 | + 'example_show_butterfly.eps') | |
| 257 | 262 | |
| 258 | 263 | #TODO: Do any testing |
| 259 | 264 | |
| ... | ... | @@ -267,7 +272,7 @@ class Test(gammalib.GPythonTestSuite): |
| 267 | 272 | """ |
| 268 | 273 | # Execute script |
| 269 | 274 | self._execute_python('show_model', |
| 270 | - args='data/crab.xml Crab example_model.eps') | |
| 275 | + args='data/crab.xml Crab example_show_model.eps') | |
| 271 | 276 | |
| 272 | 277 | #TODO: Do any testing |
| 273 | 278 | |
| ... | ... | @@ -280,7 +285,8 @@ class Test(gammalib.GPythonTestSuite): |
| 280 | 285 | Test show_pha |
| 281 | 286 | """ |
| 282 | 287 | # Execute script |
| 283 | - self._execute_python('show_pha', args='data/pha.fits example_pha.eps') | |
| 288 | + self._execute_python('show_pha', | |
| 289 | + args='data/pha.fits example_show_pha.eps') | |
| 284 | 290 | |
| 285 | 291 | #TODO: Do any testing |
| 286 | 292 | |
| ... | ... | @@ -295,7 +301,7 @@ class Test(gammalib.GPythonTestSuite): |
| 295 | 301 | # Execute script |
| 296 | 302 | self._execute_python('show_pull_evolution', |
| 297 | 303 | args='data/pull.dat Pull_Crab_Prefactor ' |
| 298 | - 'example_pull_evolution.eps') | |
| 304 | + 'example_show_pull_evolution.eps') | |
| 299 | 305 | |
| 300 | 306 | #TODO: Do any testing |
| 301 | 307 | |
| ... | ... | @@ -310,7 +316,65 @@ class Test(gammalib.GPythonTestSuite): |
| 310 | 316 | # Execute script |
| 311 | 317 | self._execute_python('show_pull_histogram', |
| 312 | 318 | args='data/pull.dat Pull_Crab_Prefactor 50 ' |
| 313 | - 'example_pull_histogram.eps') | |
| 319 | + 'example_show_pull_histogram.eps') | |
| 320 | + | |
| 321 | + #TODO: Do any testing | |
| 322 | + | |
| 323 | + # Return | |
| 324 | + return | |
| 325 | + | |
| 326 | + # Test show_response | |
| 327 | + def test_show_response(self): | |
| 328 | + """ | |
| 329 | + Test show_response | |
| 330 | + """ | |
| 331 | + # Execute script | |
| 332 | + self._execute_python('show_response', args='example_show_response.eps') | |
| 333 | + | |
| 334 | + #TODO: Do any testing | |
| 335 | + | |
| 336 | + # Return | |
| 337 | + return | |
| 338 | + | |
| 339 | + # Test show_sensitivity | |
| 340 | + def test_show_sensitivity(self): | |
| 341 | + """ | |
| 342 | + Test show_sensitivity | |
| 343 | + """ | |
| 344 | + # Execute script | |
| 345 | + self._execute_python('show_sensitivity', | |
| 346 | + args='data/sensitivity.dat ' | |
| 347 | + 'example_show_sensitivity.eps') | |
| 348 | + | |
| 349 | + #TODO: Do any testing | |
| 350 | + | |
| 351 | + # Return | |
| 352 | + return | |
| 353 | + | |
| 354 | + # Test show_spectrum | |
| 355 | + def test_show_spectrum(self): | |
| 356 | + """ | |
| 357 | + Test show_spectrum | |
| 358 | + """ | |
| 359 | + # Execute script | |
| 360 | + self._execute_python('show_spectrum', | |
| 361 | + args='data/spectrum.fits ' | |
| 362 | + 'example_show_spectrum.eps') | |
| 363 | + | |
| 364 | + #TODO: Do any testing | |
| 365 | + | |
| 366 | + # Return | |
| 367 | + return | |
| 368 | + | |
| 369 | + # Test show_ts_distribution | |
| 370 | + def test_show_ts_distribution(self): | |
| 371 | + """ | |
| 372 | + Test show_ts_distribution | |
| 373 | + """ | |
| 374 | + # Execute script | |
| 375 | + self._execute_python('show_ts_distribution', | |
| 376 | + args='data/ts.dat -n 100 ' | |
| 377 | + '-f example_show_ts_distribution.eps') | |
| 314 | 378 | |
| 315 | 379 | #TODO: Do any testing |
| 316 | 380 | ... | ... |