In [1]:
%pylab inline
In [2]:
import os
import matplotlib.gridspec as gridspec
from scipy.interpolate import interp1d
from astropy.io import fits
In [3]:
# top level dir containg boss data, this notebook assumes a directory layout organized as follows:
# '<boss_root>/deltas/delta-<thing_id>.fits'
# '<boss_root>/v5_7_0/spectra/lite/<plate>/spec-<plate>-<mjd>-<fiber>.fits'
boss_root = '/Users/daniel/data/boss'
In [4]:
class BOSSSpectrum:
"""
Represents a BOSS spectrum
"""
def __init__(self, loglam, flux, ivar):
self.loglam = loglam
self.flux = flux
self.ivar = ivar
def combine_pixels(self, num_combine, offset=None, ivar_weighted=True):
"""
Combines groups of consecutive pixels.
Args:
num_combine (int): number of consecutive pixels to combine
offset (None, int): pixel index to start groupings at (useful for aligning wavelength grids)
ivar_weighted (bool): indicate whether or not to combine flux using an ivar weighted average
Returns:
A new BOSSSpectrum object with combined pixels of self.
"""
# Find offset relative to fiducial wavelength grid
if offset is None:
fiducial_offset = self.get_fiducial_pixel_index_offset(self.loglam[0])
offset = fiducial_offset % num_combine
# Trim pixel arrays to multiple of num_combine
num_pixels = len(self.loglam)
num_leftover = num_pixels % num_combine
if offset > 0 or num_leftover > 0:
end_pixel = num_pixels - ((num_pixels - offset) % num_combine)
trimmed = self.select_pixels(slice(offset, end_pixel))
else:
trimmed = self
# Perform combination as a "reduction" of reshaped array dimension
loglam = trimmed.loglam.reshape((-1, num_combine))
flux = trimmed.flux.reshape((-1, num_combine))
ivar = trimmed.ivar.reshape((-1, num_combine))
loglam = np.average(loglam, axis=1)
if ivar_weighted:
flux = np.average(flux, weights=ivar, axis=1)
else:
flux = np.average(flux, axis=1)
ivar = np.sum(ivar, axis=1)
return BOSSSpectrum(loglam, flux, ivar)
def busca_combine(self, num_combine, z, forest_min=1040, forest_max=1200,
observed_min=3600, observed_max=5500, loglam_spacing=1e-4):
loglam_min_rest = np.log10(forest_min)
loglam_max_rest = np.log10(forest_max)
loglam_min = np.log10(observed_min)
loglam_max = np.log10(observed_max)
dloglam = num_combine*loglam_spacing
# calculate start and end wavelengths
loglam_start = loglam_min_rest + np.log10(1+z)
if loglam_start < loglam_min:
loglam_start = loglam_min
loglam_stop = loglam_max_rest + np.log10(1+z)
if loglam_stop > loglam_max:
loglam_stop = loglam_max
# shift to coadd grid
combined_bin_start = np.around((loglam_start-self.loglam[0])/dloglam).astype(int)
loglam_start = self.loglam[0] + combined_bin_start*dloglam
combined_bin_stop = np.around((loglam_stop-self.loglam[0])/dloglam).astype(int)
loglam_stop = self.loglam[0] + combined_bin_stop*dloglam
# calculate number of bins on combined grid
num_combined_bins = np.around((loglam_stop-loglam_start)/dloglam).astype(int)+1
#loglam = loglam_start + np.arange(num_combined_bins)*num_combined_bins
# calculate new combined bin indices of uncombinned pixels
bin_indices = np.around((self.loglam-loglam_start)/dloglam).astype(int);
# trim spectrum to pixels that will go into combined bins
bin_min = np.argmax(bin_indices >= 0)
bin_max = np.argmax(bin_indices >= num_combined_bins)
trimmed = self.select_pixels(slice(bin_min, bin_max))
# combine pixels
return trimmed.combine_pixels(num_combine, offset=0)
def select_pixels(self, selection):
"""
Returns a copy of self after applying the selection provided.
"""
return BOSSSpectrum(self.loglam[selection], self.flux[selection], self.ivar[selection])
def get_forest_spectrum(self, z, forest_min=1040.0, forest_max=1200.0):
"""
Returns a copy of self after selecting pixels in the lyman alpha forest
"""
forest_min *= 1+z
forest_max *= 1+z
forest_pixels = np.where( \
(self.loglam >= np.log10(forest_min)) & \
(self.loglam <= np.log10(forest_max)) )[0]
return self.select_pixels(forest_pixels)
@staticmethod
def get_fiducial_pixel_index_offset(coeff0, coeff1=1e-4, fid_lambda=3500.26):
"""
Returns the pixel index offset from the start of the BOSS co-add fiducial wavelength grid.
If more than 1% away from a fiducial pixel index, returns the exact offset (in pixels).
Args:
coeff0 (float): central wavelength (log10) of first pixel
coeff1 (float, optional): log10 spacing between pixel centers
Returns:
pixel index offset from the start of the fiducial wavelength grid.
"""
delta = (np.log10(fid_lambda)-coeff0)/coeff1
# round to nearest pixel index
offset = np.floor(delta+0.5).astype(int)
# if more than %1 away from pixel index, return exact offset (in pixels)
if np.any(np.fabs(delta-offset) > 0.01):
return -delta
# otherwise return nearest pixel index
return -offset
In [11]:
class DeltaLOS():
"""
Represents a delta file (delta field along a line of sight)
"""
def __init__(self, thing_id):
self.thing_id = thing_id
hdulist = fits.open(os.path.join(boss_root,'deltas','delta-%s.fits' % thing_id))
data = hdulist[1].data
self.loglam = data['loglam']
self.delta = data['delta']
self.weight = data['weight']
self.r_comov = data['r_comov']
self.cont = data['cont']
self.msha = data['msha']
self.mabs = data['mabs']
self.ivar = data['ivar']
self.flux = (1+self.delta)*self.cont*self.msha*self.mabs
header = hdulist[1].header
self.z = header['z']
self.plate = header['plate']
self.mjd = header['mjd']
self.fiber = header['fiberid']
self.ra = header['ra']
self.dec = header['dec']
In [18]:
class LiteSpec():
"""
Represents a lite spec file (single target's observation, co-add only)
"""
def __init__(self, plate, mjd, fiber):
self.plate = plate
self.mjd = mjd
self.fiber = fiber
filename = 'spec-%s-%s-%04d.fits' % (plate, mjd, fiber)
hdulist = fits.open(os.path.join(boss_root,'v5_7_0','spectra','lite',str(plate),filename))
data = hdulist[1].data
self.loglam = data['loglam']
self.flux = data['flux']
self.ivar = data['ivar']
self.ra = hdulist[2].data['ra']
self.dec = hdulist[2].data['dec']
In [19]:
def compare_binning(thing_id, num_combine=3, offset=None):
# Load delta los data and raw BOSS spectrum
delta_los = DeltaLOS(thing_id)
print 'RA, Dec: %f, %f' % (delta_los.ra, delta_los.dec)
lite_spec = LiteSpec(delta_los.plate, delta_los.mjd, delta_los.fiber)
print 'RA, Dec: %f, %f' % (lite_spec.ra, lite_spec.dec)
boss_spectrum = BOSSSpectrum(lite_spec.loglam, lite_spec.flux, lite_spec.ivar)
print 'Plate, mjd, fiber:', delta_los.plate, delta_los.mjd, delta_los.fiber
# combine pixels from raw spectrum into "analysis" pixels
analysis_spectrum = boss_spectrum.combine_pixels(num_combine, offset=offset, ivar_weighted=True)
forest_spectrum = analysis_spectrum.get_forest_spectrum(delta_los.z)
busca_spectrum = boss_spectrum.busca_combine(num_combine, delta_los.z)
test_spectrum = busca_spectrum#.get_forest_spectrum(delta_los.z)
# plot results
# plot forest/spec limits
z = delta_los.z
forest_min = 1040*(1+z)
forest_max = 1200*(1+z)
observed_min = 10**boss_spectrum.loglam[0]
observed_max = 10**boss_spectrum.loglam[-1]
fig = plt.figure(figsize=(14, 6))
fig.suptitle('THING_ID = %s, Z_VI = %.2f' % (thing_id, z))
gs = gridspec.GridSpec(2, 2, height_ratios=[2,1], width_ratios=[1,1], hspace=0.1, wspace=0.05)
ax1 = plt.subplot(gs[0])
# plot the raw boss spectrum
plt.plot(10**boss_spectrum.loglam, boss_spectrum.flux,
color='gray', label='Raw Spectrum', lw=0, marker='+')
# plot busca forest
s = test_spectrum.ivar > 0
plt.scatter(10**test_spectrum.loglam[s], test_spectrum.flux[s],
color='purple', marker='o', facecolors='none', label='Busca Pixels')
# plot the delta los spectrum
s = delta_los.ivar > 0
plt.errorbar(10**delta_los.loglam[s], delta_los.flux[s], yerr=1/np.sqrt(delta_los.ivar)[s],
color='blue', marker='x', fmt='o', alpha=.7, label='FPG Pixels')
# plot the combined and trimmed forest
s = forest_spectrum.ivar > 0
plt.errorbar(10**forest_spectrum.loglam[s], forest_spectrum.flux[s], yerr=1/np.sqrt(forest_spectrum.ivar[s]),
color='red', marker='x', fmt='o', alpha=.7, label='IVAR Weighted Pixels')
plt.axvline(forest_min, color='black', linestyle='--')
plt.axvline(forest_max, color='black', linestyle='--')
plt.axvline(observed_min, color='black', linestyle='-.')
plt.grid()
plt.ylabel('Observed Flux')
# hide x-axis ticklabels
plt.setp(ax1.get_xticklabels(), visible=False)
ax2 = plt.subplot(gs[1], sharey=ax1)
# plot the raw boss spectrum
plt.plot(10**boss_spectrum.loglam, boss_spectrum.flux,
color='gray', label='Raw Spectrum', lw=0, marker='+')
# plot busca forest
s = test_spectrum.ivar > 0
plt.scatter(10**test_spectrum.loglam[s], test_spectrum.flux[s],
color='purple', marker='o', facecolors='none', label='Busca Pixels')
# plot the delta los spectrum
s = delta_los.ivar > 0
plt.errorbar(10**delta_los.loglam[s], delta_los.flux[s], yerr=1/np.sqrt(delta_los.ivar)[s],
color='blue', marker='x', fmt='o', alpha=.7, label='FPG Pixels')
# plot the combined and trimmed forest
s = forest_spectrum.ivar > 0
plt.errorbar(10**forest_spectrum.loglam[s], forest_spectrum.flux[s], yerr=1/np.sqrt(forest_spectrum.ivar[s]),
color='red', marker='x', fmt='o', alpha=.7, label='IVAR Weighted Pixels')
plt.axvline(forest_max, color='black', linestyle='--')
plt.grid()
plt.ylim(np.min(forest_spectrum.flux-1/np.sqrt(forest_spectrum.ivar))*0.9,
np.max(forest_spectrum.flux+1/np.sqrt(forest_spectrum.ivar))*1.1)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax2.get_yticklabels(), visible=False)
# switch to lower subplot, share x axis
ax3 = plt.subplot(gs[2], sharex=ax1)
# create interpolation of delta field fluxes so we can evaulate using grid of our combined pixels
delta_los_interp = interp1d(delta_los.loglam, delta_los.flux, bounds_error=False)
residual_flux = forest_spectrum.flux - delta_los_interp(forest_spectrum.loglam)
# plot residual
plt.step(10**forest_spectrum.loglam, residual_flux,
color='black', where='mid', lw=.5)
plt.step(10**test_spectrum.loglam, test_spectrum.flux - delta_los_interp(test_spectrum.loglam),
color='purple', where='mid', lw=.5)
# plot forest/spec limits
plt.axvline(forest_min, color='black', linestyle='--')
plt.axvline(forest_max, color='black', linestyle='--')
plt.axvline(observed_min, color='black', linestyle='-.')
plt.grid()
plot_min = max(forest_min, observed_min)-10
plot_max = min(plot_min+80, forest_max+20)
plt.xlim(plot_min, plot_max)
plt.xlabel('Observed Wavelength ($\AA$)')
plt.ylabel('Residual Flux')
ax4 = plt.subplot(gs[3], sharex=ax2, sharey=ax3)
# plot residual
plt.step(10**forest_spectrum.loglam, residual_flux,
color='black', where='mid', lw=.5)
plt.step(10**test_spectrum.loglam, test_spectrum.flux - delta_los_interp(test_spectrum.loglam),
color='purple', where='mid', lw=.5)
# plot forest/spec limits
plt.axvline(forest_min, color='black', linestyle='--')
plt.axvline(forest_max, color='black', linestyle='--')
plt.axvline(observed_min, color='black', linestyle='-.')
plt.grid()
plt.xlim(forest_max+10 - (plot_max - plot_min), forest_max+10)
plt.xlabel('Observed Wavelength ($\AA$)')
plt.setp(ax4.get_yticklabels(), visible=False)
In [20]:
compare_binning(66692618, num_combine=3)
In [17]:
compare_binning(401440280, num_combine=3)
In [10]:
compare_binning(401440280, num_combine=3, offset=2)
In [10]: