Luigi Tibaldo / ctools

Commit b1e044a447312eaf6cf3f3ad5b9dfb97c85c27a2

Authored by Jürgen Knödlseder
1 parent 1cbb5663

show_spectrum.py now adapts to energy scale of data and shows filename in title

Inline Side-by-side
Showing 1 changed file with 39 additions and 17 deletions
examples/show_spectrum.py
  View file @b1e044a
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # Display spectrum generated by csspec
4 4 #
5   -# Copyright (C) 2015-2020 Juergen Knoedlseder
  5 +# Copyright (C) 2015-2021 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -18,6 +18,7 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
  21 +import os
21 22 import sys
22 23 try:
23 24 import matplotlib.pyplot as plt
... ... @@ -70,7 +71,7 @@ def get_spectrum_fits(fits):
70 71 ----------
71 72 fits : `~gammalib.GFits`
72 73 Spectral GFits object
73   -
  74 +
74 75 Returns
75 76 -------
76 77 spec : dict
... ... @@ -168,7 +169,7 @@ def get_spectrum_fits(fits):
168 169 # ========================= #
169 170 # Plot delta log-likelihood #
170 171 # ========================= #
171   -def plot_dloglike(spec):
  172 +def plot_dloglike(ax, spec):
172 173 """
173 174 Plot delta log-likelihood
174 175  
... ... @@ -187,7 +188,7 @@ def plot_dloglike(spec):
187 188  
188 189 # Compute the logarithmic center of each bin
189 190 fluxpnts = 10 ** ((y_bounds[1:]+y_bounds[:-1]) / 2.0)
190   -
  191 +
191 192 # Scale the boundaries on the grid
192 193 y_bounds = 10 ** y_bounds
193 194  
... ... @@ -201,12 +202,12 @@ def plot_dloglike(spec):
201 202 dll_hist = [[row[i] for row in dll_hist] for i in range(steps)]
202 203  
203 204 # Plot the likelihood profiles
204   - plt.pcolormesh(np.array(spec['engs_scan']), y_bounds,
205   - np.array(dll_hist), vmin=-10.0, vmax=0.0,
206   - cmap='Reds', alpha=0.9)
207   - cbar = plt.colorbar()
  205 + c = ax.pcolormesh(np.array(spec['engs_scan']), y_bounds,
  206 + np.array(dll_hist), vmin=-10.0, vmax=0.0,
  207 + cmap='Reds', alpha=0.9)
  208 + cbar = plt.colorbar(c)
208 209 cbar.set_label(r'$\Delta$ log-likelihood')
209   - plt.grid()
  210 + ax.grid()
210 211  
211 212 # Return
212 213 return
... ... @@ -229,27 +230,48 @@ def plot_spectrum(filename, plotfile):
229 230 # Get spectrum parameters
230 231 spec = get_spectrum_file(filename)
231 232  
  233 + # Set energy scale of spectrum
  234 + if max(spec['energies']) < 1.0e-3:
  235 + xlabel = 'Energy (MeV)'
  236 + conv = 1.0e6
  237 + elif max(spec['energies']) < 1.0:
  238 + xlabel = 'Energy (GeV)'
  239 + conv = 1.0e3
  240 + else:
  241 + xlabel = 'Energy (TeV)'
  242 + conv = 1.0
  243 + if conv != 1.0:
  244 + spec['energies'] = [e * 1.0e6 for e in spec['energies']]
  245 + spec['ed_engs'] = [e * 1.0e6 for e in spec['ed_engs']]
  246 + spec['eu_engs'] = [e * 1.0e6 for e in spec['eu_engs']]
  247 + spec['ul_energies'] = [e * 1.0e6 for e in spec['ul_energies']]
  248 + spec['ul_ed_engs'] = [e * 1.0e6 for e in spec['ul_ed_engs']]
  249 + spec['ul_eu_engs'] = [e * 1.0e6 for e in spec['ul_eu_engs']]
  250 +
232 251 # Create the plot
233   - plt.figure()
234   - plt.loglog()
235   - plt.grid()
236   - plt.xlabel('Energy (TeV)')
237   - plt.ylabel(r'E$^2$ $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)')
  252 + fig = plt.figure()
  253 + fig.subplots_adjust(left=0.17, right=0.98, top=0.92, bottom=0.10)
  254 + ax = fig.add_subplot(111)
  255 + ax.loglog()
  256 + ax.grid()
  257 + ax.set_xlabel(xlabel)
  258 + ax.set_ylabel(r'E$^2$ $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)')
  259 + ax.set_title(os.path.basename(filename))
238 260  
239 261 # Plot the spectrum
240   - plt.errorbar(spec['energies'], spec['flux'],
  262 + ax.errorbar(spec['energies'], spec['flux'],
241 263 yerr=spec['e_flux'], xerr=[spec['ed_engs'], spec['eu_engs']],
242 264 fmt='ro')
243 265  
244 266 # Plot upper limits
245 267 if len(spec['ul_energies']) > 0:
246   - plt.errorbar(spec['ul_energies'], spec['ul_flux'], yerr=spec['yerr'],
  268 + ax.errorbar(spec['ul_energies'], spec['ul_flux'], yerr=spec['yerr'],
247 269 xerr=[spec['ul_ed_engs'], spec['ul_eu_engs']],
248 270 uplims=True, fmt='ro')
249 271  
250 272 # Plot log-likelihood profiles
251 273 if len(spec['dll_scan']) > 0:
252   - plot_dloglike(spec)
  274 + plot_dloglike(ax, spec)
253 275  
254 276 # Show figure
255 277 if len(plotfile) > 0:
... ...