Eirik Hatlen / ctools

Commit 74c93801ddd5758252381446308f64eb9040df15

Authored by Josh Cardenzana
1 parent 8c7309c7

(#2938) Adjust column name formatting in csspec fits files 

This ensures the SED fits files are consistent with the current suggested data format
Inline Side-by-side
Showing 2 changed files with 78 additions and 58 deletions
cscripts/csspec.py
  View file @74c9380
@@ -645,6 +645,7 @@ class csspec(ctools.csobservation): @@ -645,6 +645,7 @@ class csspec(ctools.csobservation):
645 emin = self._ebounds.emin(i) 645 emin = self._ebounds.emin(i)
646 emax = self._ebounds.emax(i) 646 emax = self._ebounds.emax(i)
647 elogmean = self._ebounds.elogmean(i) 647 elogmean = self._ebounds.elogmean(i)
  648 + e_scale = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg
648 649
649 # Select observations for energy bin 650 # Select observations for energy bin
650 obs = self._select_obs(emin, emax) 651 obs = self._select_obs(emin, emax)
@@ -654,6 +655,7 @@ class csspec(ctools.csobservation): @@ -654,6 +655,7 @@ class csspec(ctools.csobservation):
654 'energy_low': (elogmean - emin).TeV(), 655 'energy_low': (elogmean - emin).TeV(),
655 'energy_high': (emax - elogmean).TeV(), 656 'energy_high': (emax - elogmean).TeV(),
656 'flux': 0.0, 657 'flux': 0.0,
  658 + 'e2dnde': 0.0,
657 'flux_err': 0.0, 659 'flux_err': 0.0,
658 'TS': 0.0, 660 'TS': 0.0,
659 'ulimit': 0.0, 661 'ulimit': 0.0,
@@ -736,8 +738,7 @@ class csspec(ctools.csobservation): @@ -736,8 +738,7 @@ class csspec(ctools.csobservation):
736 738
737 # Compute upper limit 739 # Compute upper limit
738 if ulimit_value > 0.0: 740 if ulimit_value > 0.0:
739 - result['ulimit'] = ulimit_value * elogmean.MeV() * \  
740 - elogmean.MeV() * gammalib.MeV2erg 741 + result['ulimit'] = ulimit_value
741 742
742 # Compute differential flux and flux error 743 # Compute differential flux and flux error
743 fitted_flux = source.spectral().eval(elogmean) 744 fitted_flux = source.spectral().eval(elogmean)
@@ -758,14 +759,14 @@ class csspec(ctools.csobservation): @@ -758,14 +759,14 @@ class csspec(ctools.csobservation):
758 e_flux *= norm 759 e_flux *= norm
759 760
760 # Convert differential flux and flux error to nuFnu 761 # Convert differential flux and flux error to nuFnu
761 - elogmean2 = elogmean.MeV() * elogmean.MeV()  
762 - result['flux'] = fitted_flux * elogmean2 * gammalib.MeV2erg  
763 - result['flux_err'] = e_flux * elogmean2 * gammalib.MeV2erg 762 + result['flux'] = fitted_flux
  763 + result['e2dnde'] = fitted_flux * e_scale
  764 + result['flux_err'] = e_flux
764 765
765 # Log information 766 # Log information
766 - value = '%e +/- %e' % (result['flux'], result['flux_err']) 767 + value = '%e +/- %e' % (result['e2dnde'], result['flux_err']*e_scale)
767 if self['calc_ulim'].boolean() and result['ulimit'] > 0.0: 768 if self['calc_ulim'].boolean() and result['ulimit'] > 0.0:
768 - value += ' [< %e]' % (result['ulimit']) 769 + value += ' [< %e]' % (result['ulimit']*e_scale)
769 value += ' erg/cm2/s' 770 value += ' erg/cm2/s'
770 if self['calc_ts'].boolean() and result['TS'] > 0.0: 771 if self['calc_ts'].boolean() and result['TS'] > 0.0:
771 value += ' (TS = %.3f)' % (result['TS']) 772 value += ' (TS = %.3f)' % (result['TS'])
@@ -866,6 +867,7 @@ class csspec(ctools.csobservation): @@ -866,6 +867,7 @@ class csspec(ctools.csobservation):
866 'energy_low': (elogmean - emin).TeV(), 867 'energy_low': (elogmean - emin).TeV(),
867 'energy_high': (emax - elogmean).TeV(), 868 'energy_high': (emax - elogmean).TeV(),
868 'flux': 0.0, 869 'flux': 0.0,
  870 + 'e2dnde': 0.0,
869 'flux_err': 0.0, 871 'flux_err': 0.0,
870 'TS': 0.0, 872 'TS': 0.0,
871 'ulimit': 0.0, 873 'ulimit': 0.0,
@@ -876,8 +878,9 @@ class csspec(ctools.csobservation): @@ -876,8 +878,9 @@ class csspec(ctools.csobservation):
876 878
877 # Convert differential flux and flux error to nuFnu 879 # Convert differential flux and flux error to nuFnu
878 norm = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg 880 norm = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg
879 - result['flux'] = spectrum[i*2+1].value() * norm  
880 - result['flux_err'] = spectrum[i*2+1].error() * norm 881 + result['flux'] = spectrum[i*2+1].value()
  882 + result['e2dnde'] = spectrum[i*2+1].value() * norm
  883 + result['flux_err'] = spectrum[i*2+1].error()
881 884
882 # Compute upper flux limit 885 # Compute upper flux limit
883 ulimit_value = -1.0 886 ulimit_value = -1.0
@@ -925,8 +928,7 @@ class csspec(ctools.csobservation): @@ -925,8 +928,7 @@ class csspec(ctools.csobservation):
925 928
926 # Compute upper limit 929 # Compute upper limit
927 if ulimit_value > 0.0: 930 if ulimit_value > 0.0:
928 - result['ulimit'] = ulimit_value * elogmean.MeV() * \  
929 - elogmean.MeV() * gammalib.MeV2erg 931 + result['ulimit'] = ulimit_value
930 932
931 # Compute TS 933 # Compute TS
932 if self['calc_ts'].boolean(): 934 if self['calc_ts'].boolean():
@@ -954,9 +956,9 @@ class csspec(ctools.csobservation): @@ -954,9 +956,9 @@ class csspec(ctools.csobservation):
954 result['TS'] = 2.0 * (logL1 - logL0) 956 result['TS'] = 2.0 * (logL1 - logL0)
955 957
956 # Log information 958 # Log information
957 - value = '%e +/- %e' % (result['flux'], result['flux_err']) 959 + value = '%e +/- %e' % (result['e2dnde'], result['flux_err']*norm)
958 if self['calc_ulim'].boolean() and result['ulimit'] > 0.0: 960 if self['calc_ulim'].boolean() and result['ulimit'] > 0.0:
959 - value += ' [< %e]' % (result['ulimit']) 961 + value += ' [< %e]' % (result['ulimit']*norm)
960 value += ' erg/cm2/s' 962 value += ' erg/cm2/s'
961 if self['calc_ts'].boolean() and result['TS'] > 0.0: 963 if self['calc_ts'].boolean() and result['TS'] > 0.0:
962 value += ' (TS = %.3f)' % (result['TS']) 964 value += ' (TS = %.3f)' % (result['TS'])
@@ -994,7 +996,6 @@ class csspec(ctools.csobservation): @@ -994,7 +996,6 @@ class csspec(ctools.csobservation):
994 dll_sigmax = self['dll_sigmax'].real() 996 dll_sigmax = self['dll_sigmax'].real()
995 dll_sigstep = self['dll_sigstep'].real() 997 dll_sigstep = self['dll_sigstep'].real()
996 sigsteps = int(2 * (dll_sigmax/dll_sigstep) + 1) 998 sigsteps = int(2 * (dll_sigmax/dll_sigstep) + 1)
997 - scale = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg  
998 999
999 # Setup the source model for fitting 1000 # Setup the source model for fitting
1000 source = like.obs().models()[self['srcname'].string()] 1001 source = like.obs().models()[self['srcname'].string()]
@@ -1016,6 +1017,9 @@ class csspec(ctools.csobservation): @@ -1016,6 +1017,9 @@ class csspec(ctools.csobservation):
1016 # Store the resulting log-likelihood 1017 # Store the resulting log-likelihood
1017 norm_scan = [] 1018 norm_scan = []
1018 dloglike = [] 1019 dloglike = []
  1020 + ref_norm = norm
  1021 + if self._method == 'SLICE':
  1022 + ref_norm = spectral.eval(elogmean)
1019 1023
1020 # Compute the flux values to evaluate loglike at 1024 # Compute the flux values to evaluate loglike at
1021 log_norm = math.log10(norm) 1025 log_norm = math.log10(norm)
@@ -1042,9 +1046,9 @@ class csspec(ctools.csobservation): @@ -1042,9 +1046,9 @@ class csspec(ctools.csobservation):
1042 # Store dlikelihood & norm 1046 # Store dlikelihood & norm
1043 dloglike.append(loglike - like.obs().logL()) 1047 dloglike.append(loglike - like.obs().logL())
1044 if self._method == 'SLICE': 1048 if self._method == 'SLICE':
1045 - norm_scan.append(spectral.eval(elogmean) * scale) 1049 + norm_scan.append(spectral.eval(elogmean))
1046 else: 1050 else:
1047 - norm_scan.append(new_norm * scale) 1051 + norm_scan.append(new_norm)
1048 1052
1049 return (norm_scan, dloglike, -loglike) 1053 return (norm_scan, dloglike, -loglike)
1050 1054
@@ -1063,22 +1067,27 @@ class csspec(ctools.csobservation): @@ -1063,22 +1067,27 @@ class csspec(ctools.csobservation):
1063 energy = gammalib.GFitsTableDoubleCol('e_ref', nrows) 1067 energy = gammalib.GFitsTableDoubleCol('e_ref', nrows)
1064 energy_low = gammalib.GFitsTableDoubleCol('e_min', nrows) 1068 energy_low = gammalib.GFitsTableDoubleCol('e_min', nrows)
1065 energy_high = gammalib.GFitsTableDoubleCol('e_max', nrows) 1069 energy_high = gammalib.GFitsTableDoubleCol('e_max', nrows)
1066 - flux = gammalib.GFitsTableDoubleCol('e2dnde', nrows)  
1067 - flux_err = gammalib.GFitsTableDoubleCol('e2dnde_err', nrows) 1070 + norm = gammalib.GFitsTableDoubleCol('norm', nrows)
  1071 + norm_err = gammalib.GFitsTableDoubleCol('norm_err', nrows)
  1072 + ulim_values = gammalib.GFitsTableDoubleCol('norm_ul', nrows)
  1073 + e2dnde = gammalib.GFitsTableDoubleCol('ref_e2dnde', nrows)
  1074 + flux = gammalib.GFitsTableDoubleCol('ref_flux', nrows)
  1075 + Npred_values = gammalib.GFitsTableDoubleCol('ref_npred', nrows)
1068 TSvalues = gammalib.GFitsTableDoubleCol('ts', nrows) 1076 TSvalues = gammalib.GFitsTableDoubleCol('ts', nrows)
1069 - ulim_values = gammalib.GFitsTableDoubleCol('e2dnde_ul', nrows)  
1070 - Npred_values = gammalib.GFitsTableDoubleCol('npred', nrows)  
1071 loglike = gammalib.GFitsTableDoubleCol('loglike', nrows) 1077 loglike = gammalib.GFitsTableDoubleCol('loglike', nrows)
1072 norm_scan = gammalib.GFitsTableDoubleCol('norm_scan', nrows, ncols) 1078 norm_scan = gammalib.GFitsTableDoubleCol('norm_scan', nrows, ncols)
1073 dloglike_scan= gammalib.GFitsTableDoubleCol('dloglike_scan', nrows, ncols) 1079 dloglike_scan= gammalib.GFitsTableDoubleCol('dloglike_scan', nrows, ncols)
1074 energy.unit('TeV') 1080 energy.unit('TeV')
1075 energy_low.unit('TeV') 1081 energy_low.unit('TeV')
1076 energy_high.unit('TeV') 1082 energy_high.unit('TeV')
1077 - flux.unit('erg/cm2/s')  
1078 - flux_err.unit('erg/cm2/s')  
1079 - ulim_values.unit('erg/cm2/s') 1083 + norm.unit('')
  1084 + norm_err.unit('')
  1085 + ulim_values.unit('')
  1086 + e2dnde.unit('erg/cm2/s')
  1087 + flux.unit('counts/MeV/cm2/s')
  1088 + Npred_values.unit('counts')
1080 loglike.unit('') 1089 loglike.unit('')
1081 - norm_scan.unit('erg/cm2/s') 1090 + norm_scan.unit('')
1082 dloglike_scan.unit('') 1091 dloglike_scan.unit('')
1083 1092
1084 # File FITS table columns 1093 # File FITS table columns
@@ -1086,17 +1095,19 @@ class csspec(ctools.csobservation): @@ -1086,17 +1095,19 @@ class csspec(ctools.csobservation):
1086 energy[i] = result['energy'] 1095 energy[i] = result['energy']
1087 energy_low[i] = result['energy_low'] 1096 energy_low[i] = result['energy_low']
1088 energy_high[i] = result['energy_high'] 1097 energy_high[i] = result['energy_high']
  1098 + norm[i] = 1.0
  1099 + norm_err[i] = result['flux_err'] / result['flux']
  1100 + ulim_values[i] = result['ulimit'] / result['flux']
  1101 + e2dnde[i] = result['e2dnde']
1089 flux[i] = result['flux'] 1102 flux[i] = result['flux']
1090 - flux_err[i] = result['flux_err']  
1091 - TSvalues[i] = result['TS']  
1092 - ulim_values[i] = result['ulimit']  
1093 Npred_values[i] = result['Npred'] 1103 Npred_values[i] = result['Npred']
  1104 + TSvalues[i] = result['TS']
1094 loglike[i] = result['logL'] 1105 loglike[i] = result['logL']
1095 1106
1096 # Add likelihood scan values 1107 # Add likelihood scan values
1097 for fbin in range(ncols): 1108 for fbin in range(ncols):
1098 dloglike_scan[i,fbin] = result['dloglike'][fbin] 1109 dloglike_scan[i,fbin] = result['dloglike'][fbin]
1099 - norm_scan[i,fbin] = result['norm_scan'][fbin] 1110 + norm_scan[i,fbin] = result['norm_scan'][fbin] / result['flux']
1100 1111
1101 # Initialise FITS Table with extension "SPECTRUM" 1112 # Initialise FITS Table with extension "SPECTRUM"
1102 table = gammalib.GFitsBinTable(nrows) 1113 table = gammalib.GFitsBinTable(nrows)
@@ -1110,26 +1121,29 @@ class csspec(ctools.csobservation): @@ -1110,26 +1121,29 @@ class csspec(ctools.csobservation):
1110 table.append(energy) 1121 table.append(energy)
1111 table.append(energy_low) 1122 table.append(energy_low)
1112 table.append(energy_high) 1123 table.append(energy_high)
  1124 + table.append(norm)
  1125 + table.append(norm_err)
  1126 + table.append(ulim_values)
  1127 + table.append(e2dnde)
1113 table.append(flux) 1128 table.append(flux)
1114 - table.append(flux_err)  
1115 table.append(TSvalues) 1129 table.append(TSvalues)
1116 - table.append(ulim_values)  
1117 table.append(Npred_values) 1130 table.append(Npred_values)
1118 - 1131 +
1119 # Define the SED type 1132 # Define the SED type
1120 - sed_type = 'e2dnde'  
1121 - if ncols > 0:  
1122 - sed_type = 'likelihood'  
1123 - table.append(loglike)  
1124 - table.append(norm_scan)  
1125 - table.append(dloglike_scan)  
1126 - table.card('SED_TYPE', sed_type, 'SED type') 1133 + table.card('SED_TYPE', 'norm,e2dnde,flux,npred', 'SED type')
1127 1134
1128 # Define the upper limit confidence level 1135 # Define the upper limit confidence level
1129 ulimit = ctools.ctulimit() 1136 ulimit = ctools.ctulimit()
1130 table.card('UL_CONF', ulimit['confidence'].real(), 1137 table.card('UL_CONF', ulimit['confidence'].real(),
1131 'Confidence level of upper limits') 1138 'Confidence level of upper limits')
1132 1139
  1140 + # Add the likelihood data
  1141 + if ncols > 0:
  1142 + table.card('SED_TYPE', 'likelihood,e2dnde,flux,npred', 'Likelihood scan data')
  1143 + table.append(loglike)
  1144 + table.append(norm_scan)
  1145 + table.append(dloglike_scan)
  1146 +
1133 # Create the FITS file now 1147 # Create the FITS file now
1134 self._fits = gammalib.GFits() 1148 self._fits = gammalib.GFits()
1135 self._fits.append(table) 1149 self._fits.append(table)
examples/show_spectrum.py
  View file @74c9380
@@ -82,10 +82,11 @@ def get_spectrum_fits(fits): @@ -82,10 +82,11 @@ def get_spectrum_fits(fits):
82 c_energy = table['e_ref'] 82 c_energy = table['e_ref']
83 c_ed = table['e_min'] 83 c_ed = table['e_min']
84 c_eu = table['e_max'] 84 c_eu = table['e_max']
85 - c_flux = table['e2dnde']  
86 - c_eflux = table['e2dnde_err'] 85 + c_flux = table['ref_e2dnde']
  86 + c_norm = table['norm']
  87 + c_eflux = table['norm_err']
  88 + c_upper = table['norm_ul']
87 c_ts = table['ts'] 89 c_ts = table['ts']
88 - c_upper = table['e2dnde_ul']  
89 90
90 # Initialise arrays to be filled 91 # Initialise arrays to be filled
91 spec = { 92 spec = {
@@ -101,18 +102,22 @@ def get_spectrum_fits(fits): @@ -101,18 +102,22 @@ def get_spectrum_fits(fits):
101 'yerr' : [], 102 'yerr' : [],
102 'loglike' : [], 103 'loglike' : [],
103 'engs_scan' : [], 104 'engs_scan' : [],
104 - 'norm_scan' : [], 105 + 'e2dnde_scan' : [],
105 'dll_scan' : [] 106 'dll_scan' : []
106 } 107 }
107 108
108 # Determine if we can load the delta-log-likelihood profiles 109 # Determine if we can load the delta-log-likelihood profiles
  110 + has_sedtype = table.has_card('SED_TYPE')
109 load_dll = False 111 load_dll = False
110 - if (fits.table(1).has_card('SED_TYPE') and  
111 - fits.table(1).card('SED_TYPE').string() == 'likelihood'):  
112 - c_loglike = table['loglike']  
113 - c_norm_scan = table['norm_scan']  
114 - c_dll_scan = table['dloglike_scan']  
115 - load_dll = True 112 + if has_sedtype:
  113 + seds = table.card('SED_TYPE').string().split(',')
  114 +
  115 + # Load likelhood columns if present
  116 + if seds[0] == 'likelihood':
  117 + load_dll = True
  118 + c_loglike = table['loglike']
  119 + c_norm_scan = table['norm_scan']
  120 + c_dll_scan = table['dloglike_scan']
116 121
117 # Loop over rows of the file 122 # Loop over rows of the file
118 nrows = table.nrows() 123 nrows = table.nrows()
@@ -120,36 +125,37 @@ def get_spectrum_fits(fits): @@ -120,36 +125,37 @@ def get_spectrum_fits(fits):
120 125
121 # Get Test Statistic, flux and flux error 126 # Get Test Statistic, flux and flux error
122 ts = c_ts.real(row) 127 ts = c_ts.real(row)
123 - flx = c_flux.real(row)  
124 - e_flx = c_eflux.real(row) 128 + norm = c_norm.real(row)
  129 + flx = norm * c_flux.real(row)
  130 + e_flx = flx * c_eflux.real(row)
125 131
126 # If Test Statistic is larger than 9 and flux error is smaller than 132 # If Test Statistic is larger than 9 and flux error is smaller than
127 # flux then append flux plots ... 133 # flux then append flux plots ...
128 if ts > 9.0 and e_flx < flx: 134 if ts > 9.0 and e_flx < flx:
129 spec['energies'].append(c_energy.real(row)) 135 spec['energies'].append(c_energy.real(row))
130 - spec['flux'].append(c_flux.real(row)) 136 + spec['flux'].append(flx)
131 spec['ed_engs'].append(c_ed.real(row)) 137 spec['ed_engs'].append(c_ed.real(row))
132 spec['eu_engs'].append(c_eu.real(row)) 138 spec['eu_engs'].append(c_eu.real(row))
133 - spec['e_flux'].append(c_eflux.real(row)) 139 + spec['e_flux'].append(e_flx)
134 140
135 # ... otherwise append upper limit 141 # ... otherwise append upper limit
136 else: 142 else:
137 spec['ul_energies'].append(c_energy.real(row)) 143 spec['ul_energies'].append(c_energy.real(row))
138 - spec['ul_flux'].append(c_upper.real(row)) 144 + spec['ul_flux'].append(flx*c_upper.real(row))
139 spec['ul_ed_engs'].append(c_ed.real(row)) 145 spec['ul_ed_engs'].append(c_ed.real(row))
140 spec['ul_eu_engs'].append(c_eu.real(row)) 146 spec['ul_eu_engs'].append(c_eu.real(row))
141 147
142 # Load the delta log-likelihood values 148 # Load the delta log-likelihood values
143 if load_dll: 149 if load_dll:
144 - spec['norm_scan'] = np.zeros((nrows, c_norm_scan.elements(0)))  
145 - spec['dll_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) 150 + spec['e2dnde_scan'] = np.zeros((nrows, c_norm_scan.elements(0)))
  151 + spec['dll_scan'] = np.zeros((nrows, c_norm_scan.elements(0)))
146 152
147 for row in range(nrows): 153 for row in range(nrows):
148 loglike = c_loglike.real(row) 154 loglike = c_loglike.real(row)
149 spec['loglike'].append(loglike) 155 spec['loglike'].append(loglike)
150 spec['engs_scan'].append(c_energy.real(row) - c_ed.real(row)) 156 spec['engs_scan'].append(c_energy.real(row) - c_ed.real(row))
151 for col in range(c_norm_scan.elements(row)): 157 for col in range(c_norm_scan.elements(row)):
152 - spec['norm_scan'][row,col] = c_norm_scan.real(row,col) 158 + spec['e2dnde_scan'][row,col] = c_norm_scan.real(row,col) * c_flux.real(row)
153 spec['dll_scan'][row,col] = c_dll_scan.real(row,col) 159 spec['dll_scan'][row,col] = c_dll_scan.real(row,col)
154 spec['engs_scan'].append(c_energy.real(nrows-1)+c_eu.real(nrows-1)) 160 spec['engs_scan'].append(c_energy.real(nrows-1)+c_eu.real(nrows-1))
155 161
@@ -183,10 +189,10 @@ def plot_dloglike(spec): @@ -183,10 +189,10 @@ def plot_dloglike(spec):
183 y_bounds = 10 ** y_bounds 189 y_bounds = 10 ** y_bounds
184 190
185 dll_hist = [] 191 dll_hist = []
186 - for ebin in range(len(spec['norm_scan'])): 192 + for ebin in range(len(spec['e2dnde_scan'])):
187 193
188 # Interpolate the dlogLike values 194 # Interpolate the dlogLike values
189 - dll_interp = np.interp(fluxpnts, spec['norm_scan'][ebin], 195 + dll_interp = np.interp(fluxpnts, spec['e2dnde_scan'][ebin],
190 spec['dll_scan'][ebin]) 196 spec['dll_scan'][ebin])
191 dll_hist.append(dll_interp) 197 dll_hist.append(dll_interp)
192 198