Commit 1d30baa73e8507021e6b515dac3499c2f1d816c4

Authored by Jürgen Knödlseder
1 parent 152493f8

Support handling of difference unit in show_lightcurve.py script

Showing 1 changed file with 20 additions and 4 deletions
examples/show_lightcurve.py
@@ -47,12 +47,28 @@ def plot_lightcurve(filename, plotfile): @@ -47,12 +47,28 @@ def plot_lightcurve(filename, plotfile):
47 # Read spectrum file 47 # Read spectrum file
48 fits = gammalib.GFits(filename) 48 fits = gammalib.GFits(filename)
49 table = fits.table(1) 49 table = fits.table(1)
  50 +
  51 + # Extract standard columns
50 c_mjd = table['MJD'] 52 c_mjd = table['MJD']
51 c_emjd = table['e_MJD'] 53 c_emjd = table['e_MJD']
52 - c_flux = table['EnergyFlux']  
53 - c_eflux = table['e_EnergyFlux']  
54 c_ts = table['TS'] 54 c_ts = table['TS']
55 - c_upper = table['EFluxUpperLimit'] 55 +
  56 + # Extract columns dependent on flux type
  57 + if table.contains('EnergyFlux'):
  58 + c_flux = table['EnergyFlux']
  59 + c_eflux = table['e_EnergyFlux']
  60 + c_upper = table['EFluxUpperLimit']
  61 + ylabel = r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)'
  62 + elif table.contains('PhotonFlux'):
  63 + c_flux = table['PhotonFlux']
  64 + c_eflux = table['e_PhotonFlux']
  65 + c_upper = table['FluxUpperLimit']
  66 + ylabel = r'N(E) (ph cm$^{-2}$ s$^{-1}$)'
  67 + else:
  68 + c_flux = table['Prefactor']
  69 + c_eflux = table['e_Prefactor']
  70 + c_upper = table['DiffUpperLimit']
  71 + ylabel = r'dN/dE (cm$^{-2}$ s$^{-1}$ MeV$^{-1}$)'
56 72
57 # Initialise arrays to be filled 73 # Initialise arrays to be filled
58 mjd = [] 74 mjd = []
@@ -97,7 +113,7 @@ def plot_lightcurve(filename, plotfile): @@ -97,7 +113,7 @@ def plot_lightcurve(filename, plotfile):
97 plt.errorbar(ul_mjd, ul_flux, xerr=[ul_e_mjd, ul_e_mjd], 113 plt.errorbar(ul_mjd, ul_flux, xerr=[ul_e_mjd, ul_e_mjd],
98 yerr=ul_e_flux, uplims=True, fmt='ro') 114 yerr=ul_e_flux, uplims=True, fmt='ro')
99 plt.xlabel('MJD (days)') 115 plt.xlabel('MJD (days)')
100 - plt.ylabel(r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)') 116 + plt.ylabel(ylabel)
101 117
102 # Show figure 118 # Show figure
103 if len(plotfile) > 0: 119 if len(plotfile) > 0: