From 74c93801ddd5758252381446308f64eb9040df15 Mon Sep 17 00:00:00 2001 From: Josh Cardenzana Date: Tue, 10 Dec 2019 09:16:26 +0100 Subject: [PATCH] (#2938) Adjust column name formatting in csspec fits files --- cscripts/csspec.py | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------- examples/show_spectrum.py | 46 ++++++++++++++++++++++++++-------------------- 2 files changed, 78 insertions(+), 58 deletions(-) diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 511841d..8807dba 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -645,6 +645,7 @@ class csspec(ctools.csobservation): emin = self._ebounds.emin(i) emax = self._ebounds.emax(i) elogmean = self._ebounds.elogmean(i) + e_scale = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg # Select observations for energy bin obs = self._select_obs(emin, emax) @@ -654,6 +655,7 @@ class csspec(ctools.csobservation): 'energy_low': (elogmean - emin).TeV(), 'energy_high': (emax - elogmean).TeV(), 'flux': 0.0, + 'e2dnde': 0.0, 'flux_err': 0.0, 'TS': 0.0, 'ulimit': 0.0, @@ -736,8 +738,7 @@ class csspec(ctools.csobservation): # Compute upper limit if ulimit_value > 0.0: - result['ulimit'] = ulimit_value * elogmean.MeV() * \ - elogmean.MeV() * gammalib.MeV2erg + result['ulimit'] = ulimit_value # Compute differential flux and flux error fitted_flux = source.spectral().eval(elogmean) @@ -758,14 +759,14 @@ class csspec(ctools.csobservation): e_flux *= norm # Convert differential flux and flux error to nuFnu - elogmean2 = elogmean.MeV() * elogmean.MeV() - result['flux'] = fitted_flux * elogmean2 * gammalib.MeV2erg - result['flux_err'] = e_flux * elogmean2 * gammalib.MeV2erg + result['flux'] = fitted_flux + result['e2dnde'] = fitted_flux * e_scale + result['flux_err'] = e_flux # Log information - value = '%e +/- %e' % (result['flux'], result['flux_err']) + value = '%e +/- %e' % (result['e2dnde'], result['flux_err']*e_scale) if self['calc_ulim'].boolean() and result['ulimit'] > 0.0: - value += ' [< %e]' % (result['ulimit']) + value += ' [< %e]' % (result['ulimit']*e_scale) value += ' erg/cm2/s' if self['calc_ts'].boolean() and result['TS'] > 0.0: value += ' (TS = %.3f)' % (result['TS']) @@ -866,6 +867,7 @@ class csspec(ctools.csobservation): 'energy_low': (elogmean - emin).TeV(), 'energy_high': (emax - elogmean).TeV(), 'flux': 0.0, + 'e2dnde': 0.0, 'flux_err': 0.0, 'TS': 0.0, 'ulimit': 0.0, @@ -876,8 +878,9 @@ class csspec(ctools.csobservation): # Convert differential flux and flux error to nuFnu norm = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg - result['flux'] = spectrum[i*2+1].value() * norm - result['flux_err'] = spectrum[i*2+1].error() * norm + result['flux'] = spectrum[i*2+1].value() + result['e2dnde'] = spectrum[i*2+1].value() * norm + result['flux_err'] = spectrum[i*2+1].error() # Compute upper flux limit ulimit_value = -1.0 @@ -925,8 +928,7 @@ class csspec(ctools.csobservation): # Compute upper limit if ulimit_value > 0.0: - result['ulimit'] = ulimit_value * elogmean.MeV() * \ - elogmean.MeV() * gammalib.MeV2erg + result['ulimit'] = ulimit_value # Compute TS if self['calc_ts'].boolean(): @@ -954,9 +956,9 @@ class csspec(ctools.csobservation): result['TS'] = 2.0 * (logL1 - logL0) # Log information - value = '%e +/- %e' % (result['flux'], result['flux_err']) + value = '%e +/- %e' % (result['e2dnde'], result['flux_err']*norm) if self['calc_ulim'].boolean() and result['ulimit'] > 0.0: - value += ' [< %e]' % (result['ulimit']) + value += ' [< %e]' % (result['ulimit']*norm) value += ' erg/cm2/s' if self['calc_ts'].boolean() and result['TS'] > 0.0: value += ' (TS = %.3f)' % (result['TS']) @@ -994,7 +996,6 @@ class csspec(ctools.csobservation): dll_sigmax = self['dll_sigmax'].real() dll_sigstep = self['dll_sigstep'].real() sigsteps = int(2 * (dll_sigmax/dll_sigstep) + 1) - scale = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg # Setup the source model for fitting source = like.obs().models()[self['srcname'].string()] @@ -1016,6 +1017,9 @@ class csspec(ctools.csobservation): # Store the resulting log-likelihood norm_scan = [] dloglike = [] + ref_norm = norm + if self._method == 'SLICE': + ref_norm = spectral.eval(elogmean) # Compute the flux values to evaluate loglike at log_norm = math.log10(norm) @@ -1042,9 +1046,9 @@ class csspec(ctools.csobservation): # Store dlikelihood & norm dloglike.append(loglike - like.obs().logL()) if self._method == 'SLICE': - norm_scan.append(spectral.eval(elogmean) * scale) + norm_scan.append(spectral.eval(elogmean)) else: - norm_scan.append(new_norm * scale) + norm_scan.append(new_norm) return (norm_scan, dloglike, -loglike) @@ -1063,22 +1067,27 @@ class csspec(ctools.csobservation): energy = gammalib.GFitsTableDoubleCol('e_ref', nrows) energy_low = gammalib.GFitsTableDoubleCol('e_min', nrows) energy_high = gammalib.GFitsTableDoubleCol('e_max', nrows) - flux = gammalib.GFitsTableDoubleCol('e2dnde', nrows) - flux_err = gammalib.GFitsTableDoubleCol('e2dnde_err', nrows) + norm = gammalib.GFitsTableDoubleCol('norm', nrows) + norm_err = gammalib.GFitsTableDoubleCol('norm_err', nrows) + ulim_values = gammalib.GFitsTableDoubleCol('norm_ul', nrows) + e2dnde = gammalib.GFitsTableDoubleCol('ref_e2dnde', nrows) + flux = gammalib.GFitsTableDoubleCol('ref_flux', nrows) + Npred_values = gammalib.GFitsTableDoubleCol('ref_npred', nrows) TSvalues = gammalib.GFitsTableDoubleCol('ts', nrows) - ulim_values = gammalib.GFitsTableDoubleCol('e2dnde_ul', nrows) - Npred_values = gammalib.GFitsTableDoubleCol('npred', nrows) loglike = gammalib.GFitsTableDoubleCol('loglike', nrows) norm_scan = gammalib.GFitsTableDoubleCol('norm_scan', nrows, ncols) dloglike_scan= gammalib.GFitsTableDoubleCol('dloglike_scan', nrows, ncols) energy.unit('TeV') energy_low.unit('TeV') energy_high.unit('TeV') - flux.unit('erg/cm2/s') - flux_err.unit('erg/cm2/s') - ulim_values.unit('erg/cm2/s') + norm.unit('') + norm_err.unit('') + ulim_values.unit('') + e2dnde.unit('erg/cm2/s') + flux.unit('counts/MeV/cm2/s') + Npred_values.unit('counts') loglike.unit('') - norm_scan.unit('erg/cm2/s') + norm_scan.unit('') dloglike_scan.unit('') # File FITS table columns @@ -1086,17 +1095,19 @@ class csspec(ctools.csobservation): energy[i] = result['energy'] energy_low[i] = result['energy_low'] energy_high[i] = result['energy_high'] + norm[i] = 1.0 + norm_err[i] = result['flux_err'] / result['flux'] + ulim_values[i] = result['ulimit'] / result['flux'] + e2dnde[i] = result['e2dnde'] flux[i] = result['flux'] - flux_err[i] = result['flux_err'] - TSvalues[i] = result['TS'] - ulim_values[i] = result['ulimit'] Npred_values[i] = result['Npred'] + TSvalues[i] = result['TS'] loglike[i] = result['logL'] # Add likelihood scan values for fbin in range(ncols): dloglike_scan[i,fbin] = result['dloglike'][fbin] - norm_scan[i,fbin] = result['norm_scan'][fbin] + norm_scan[i,fbin] = result['norm_scan'][fbin] / result['flux'] # Initialise FITS Table with extension "SPECTRUM" table = gammalib.GFitsBinTable(nrows) @@ -1110,26 +1121,29 @@ class csspec(ctools.csobservation): table.append(energy) table.append(energy_low) table.append(energy_high) + table.append(norm) + table.append(norm_err) + table.append(ulim_values) + table.append(e2dnde) table.append(flux) - table.append(flux_err) table.append(TSvalues) - table.append(ulim_values) table.append(Npred_values) - + # Define the SED type - sed_type = 'e2dnde' - if ncols > 0: - sed_type = 'likelihood' - table.append(loglike) - table.append(norm_scan) - table.append(dloglike_scan) - table.card('SED_TYPE', sed_type, 'SED type') + table.card('SED_TYPE', 'norm,e2dnde,flux,npred', 'SED type') # Define the upper limit confidence level ulimit = ctools.ctulimit() table.card('UL_CONF', ulimit['confidence'].real(), 'Confidence level of upper limits') + # Add the likelihood data + if ncols > 0: + table.card('SED_TYPE', 'likelihood,e2dnde,flux,npred', 'Likelihood scan data') + table.append(loglike) + table.append(norm_scan) + table.append(dloglike_scan) + # Create the FITS file now self._fits = gammalib.GFits() self._fits.append(table) diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py index 5f90b51..c77b33c 100755 --- a/examples/show_spectrum.py +++ b/examples/show_spectrum.py @@ -82,10 +82,11 @@ def get_spectrum_fits(fits): c_energy = table['e_ref'] c_ed = table['e_min'] c_eu = table['e_max'] - c_flux = table['e2dnde'] - c_eflux = table['e2dnde_err'] + c_flux = table['ref_e2dnde'] + c_norm = table['norm'] + c_eflux = table['norm_err'] + c_upper = table['norm_ul'] c_ts = table['ts'] - c_upper = table['e2dnde_ul'] # Initialise arrays to be filled spec = { @@ -101,18 +102,22 @@ def get_spectrum_fits(fits): 'yerr' : [], 'loglike' : [], 'engs_scan' : [], - 'norm_scan' : [], + 'e2dnde_scan' : [], 'dll_scan' : [] } # Determine if we can load the delta-log-likelihood profiles + has_sedtype = table.has_card('SED_TYPE') load_dll = False - if (fits.table(1).has_card('SED_TYPE') and - fits.table(1).card('SED_TYPE').string() == 'likelihood'): - c_loglike = table['loglike'] - c_norm_scan = table['norm_scan'] - c_dll_scan = table['dloglike_scan'] - load_dll = True + if has_sedtype: + seds = table.card('SED_TYPE').string().split(',') + + # Load likelhood columns if present + if seds[0] == 'likelihood': + load_dll = True + c_loglike = table['loglike'] + c_norm_scan = table['norm_scan'] + c_dll_scan = table['dloglike_scan'] # Loop over rows of the file nrows = table.nrows() @@ -120,36 +125,37 @@ def get_spectrum_fits(fits): # Get Test Statistic, flux and flux error ts = c_ts.real(row) - flx = c_flux.real(row) - e_flx = c_eflux.real(row) + norm = c_norm.real(row) + flx = norm * c_flux.real(row) + e_flx = flx * c_eflux.real(row) # If Test Statistic is larger than 9 and flux error is smaller than # flux then append flux plots ... if ts > 9.0 and e_flx < flx: spec['energies'].append(c_energy.real(row)) - spec['flux'].append(c_flux.real(row)) + spec['flux'].append(flx) spec['ed_engs'].append(c_ed.real(row)) spec['eu_engs'].append(c_eu.real(row)) - spec['e_flux'].append(c_eflux.real(row)) + spec['e_flux'].append(e_flx) # ... otherwise append upper limit else: spec['ul_energies'].append(c_energy.real(row)) - spec['ul_flux'].append(c_upper.real(row)) + spec['ul_flux'].append(flx*c_upper.real(row)) spec['ul_ed_engs'].append(c_ed.real(row)) spec['ul_eu_engs'].append(c_eu.real(row)) # Load the delta log-likelihood values if load_dll: - spec['norm_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) - spec['dll_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) + spec['e2dnde_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) + spec['dll_scan'] = np.zeros((nrows, c_norm_scan.elements(0))) for row in range(nrows): loglike = c_loglike.real(row) spec['loglike'].append(loglike) spec['engs_scan'].append(c_energy.real(row) - c_ed.real(row)) for col in range(c_norm_scan.elements(row)): - spec['norm_scan'][row,col] = c_norm_scan.real(row,col) + spec['e2dnde_scan'][row,col] = c_norm_scan.real(row,col) * c_flux.real(row) spec['dll_scan'][row,col] = c_dll_scan.real(row,col) spec['engs_scan'].append(c_energy.real(nrows-1)+c_eu.real(nrows-1)) @@ -183,10 +189,10 @@ def plot_dloglike(spec): y_bounds = 10 ** y_bounds dll_hist = [] - for ebin in range(len(spec['norm_scan'])): + for ebin in range(len(spec['e2dnde_scan'])): # Interpolate the dlogLike values - dll_interp = np.interp(fluxpnts, spec['norm_scan'][ebin], + dll_interp = np.interp(fluxpnts, spec['e2dnde_scan'][ebin], spec['dll_scan'][ebin]) dll_hist.append(dll_interp) -- libgit2 0.21.4