Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor S/N QA #459

Merged
merged 3 commits into from
Dec 10, 2018
Merged

Refactor S/N QA #459

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3666467
refactor underway
profxj Dec 4, 2018
1221be2
full run
profxj Dec 4, 2018
5495961
changes
profxj Dec 4, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion doc/changes.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ desisim change log
0.31.1 (unreleased)
-------------------

* No changes yet
* Refactor S/N qa to load cframes only once (also updates OII for new TRUTH table)

0.31.0 (2018-11-08)
-------------------
Expand Down
26 changes: 18 additions & 8 deletions py/desisim/scripts/qa_s2n.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,17 +22,18 @@ def parse(options=None):
if options is None:
args = parser.parse_args()
else:
args = parser.parse_args(options)
args = parser.parse_args(namespace=options)
return args

def main(args):
import os.path
import numpy as np
import pdb

import desispec.io
from desiutil.log import get_logger
from desisim.spec_qa.s2n import load_s2n_values, obj_s2n_wave, obj_s2n_z
from desisim.spec_qa.s2n import load_all_s2n_values
from desisim.spec_qa.s2n import parse_s2n_values

# Initialize
if args.qaprod_dir is not None:
Expand All @@ -43,8 +44,16 @@ def main(args):
# Grab nights
nights = desispec.io.get_nights()


# Load all s2n (once)
all_s2n_values = []
channels = ['b', 'r', 'z']
for channel in channels:
print("Loading S/N for channel {}".format(channel))
all_s2n_values.append(load_all_s2n_values(nights, channel))

# Loop on channel
for channel in ['b', 'r', 'z']:
for ss, channel in enumerate(channels):
if channel == 'b':
wv_bins = np.arange(3570., 5700., 20.)
elif channel == 'r':
Expand All @@ -53,7 +62,7 @@ def main(args):
wv_bins = np.arange(7500., 9800., 20.)
z_bins = np.linspace(1.0, 1.6, 100) # z camera
else:
pdb.set_trace()
raise IOError("Bad channel value: {}".format(channel))
# Loop on OBJTYPE
for objtype in ['ELG', 'LRG', 'QSO']:
if objtype == 'ELG':
Expand All @@ -63,16 +72,17 @@ def main(args):
flux_bins = np.linspace(16., 22., 6)
elif objtype == 'QSO':
flux_bins = np.linspace(15., 24., 6)
# Load
s2n_values = load_s2n_values(objtype, nights, channel)#, sub_exposures=exposures)
# Parse
fdict = all_s2n_values[ss]
s2n_dict = parse_s2n_values(objtype, fdict)
# Plot
outfile = qaprod_dir+'/QA_s2n_{:s}_{:s}.png'.format(objtype, channel)
desispec.io.util.makepath(outfile)
obj_s2n_wave(s2n_values, wv_bins, flux_bins, objtype, outfile=outfile)
obj_s2n_wave(s2n_dict, wv_bins, flux_bins, objtype, outfile=outfile)
# S/N vs. z for ELG
if (channel == 'z') & (objtype=='ELG'):
outfile = qaprod_dir+'/QA_s2n_{:s}_{:s}_redshift.png'.format(objtype,channel)
desispec.io.util.makepath(outfile)
obj_s2n_z(s2n_values, z_bins, oii_bins, objtype, outfile=outfile)
obj_s2n_z(s2n_dict, z_bins, oii_bins, objtype, outfile=outfile)


126 changes: 125 additions & 1 deletion py/desisim/spec_qa/s2n.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
# matplotlib.use('Agg')

import numpy as np
import sys, os, pdb, glob
import sys, os, glob

from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
Expand All @@ -23,7 +23,129 @@

log = get_logger()


def load_all_s2n_values(nights, channel, sub_exposures=None):
"""
Calculate S/N values for a set of spectra from an input list of nights

Args:
nights: list
channel: str ('b','r','z')
sub_exposures:

Returns:
fdict: dict
Contains all the S/N info for all nights in the given channel

"""
fdict = dict(waves=[], s2n=[], fluxes=[], exptime=[], OII=[], objtype=[])
for night in nights:
if sub_exposures is not None:
exposures = sub_exposures
else:
exposures = get_exposures(night)#, raw=True)
for exposure in exposures:
fibermap_path = findfile(filetype='fibermap', night=night, expid=exposure)
fibermap_data = read_fibermap(fibermap_path)
flavor = fibermap_data.meta['FLAVOR']
if flavor.lower() in ('arc', 'flat', 'bias'):
log.debug('Skipping calibration {} exposure {:08d}'.format(flavor, exposure))
continue
# Load simspec
simspec_file = fibermap_path.replace('fibermap', 'simspec')
log.debug('Getting truth from {}'.format(simspec_file))
sps_hdu = fits.open(simspec_file)
sps_tab = Table(sps_hdu['TRUTH'].data,masked=True)

#- Get OIIFLUX from separate HDU and join
'''
if 'TRUTH_ELG' in sps_hdu:
elg_truth = Table(sps_hdu['TRUTH_ELG'].data)
sps_tab = join(sps_tab, elg_truth['TARGETID', 'OIIFLUX'],
keys='TARGETID', join_type='left')
else:
sps_tab['OIIFLUX'] = 0.0
'''

sps_hdu.close()
#objs = sps_tab['TEMPLATETYPE'] == objtype
#if np.sum(objs) == 0:
# continue

# Load spectra (flux or not fluxed; should not matter)
for ii in range(10):
camera = channel+str(ii)
cframe_path = findfile(filetype='cframe', night=night, expid=exposure, camera=camera)
try:
log.debug('Reading from {}'.format(cframe_path))
cframe = read_frame(cframe_path)
except:
log.warn("Cannot find file: {:s}".format(cframe_path))
continue
# Calculate S/N per Ang
dwave = cframe.wave - np.roll(cframe.wave,1)
dwave[0] = dwave[1]
# Calculate
s2n = cframe.flux * np.sqrt(cframe.ivar) / np.sqrt(dwave)
#s2n = cframe.flux[iobjs,:] * np.sqrt(cframe.ivar[iobjs,:]) / np.sqrt(dwave)
# Save
fdict['objtype'].append(sps_tab['TEMPLATETYPE'].data[cframe.fibers])
fdict['waves'].append(cframe.wave)
fdict['s2n'].append(s2n)
fdict['fluxes'].append(sps_tab['MAG'].data[cframe.fibers])
fdict['OII'].append(sps_tab['OIIFLUX'].data[cframe.fibers])
fdict['exptime'].append(cframe.meta['EXPTIME'])
# Return
return fdict

def parse_s2n_values(objtype, fdict):
"""
Parse the input set of S/N measurements on objtype

Args:
objtype: str
fdict: dict
Contains all the S/N info for all nights in a given channel

Returns:
pdict: dict
Contains all the S/N info for the given objtype

"""
pdict = dict(waves=[], s2n=[], fluxes=[], exptime=[], OII=[], objtype=[])
# Loop on all the entries
for ss, wave in enumerate(fdict['waves']):
objs = fdict['objtype'][ss] == objtype
if np.sum(objs) == 0:
continue
iobjs = np.where(objs)[0]
# Parse/Save
pdict['waves'].append(wave)
pdict['s2n'].append(fdict['s2n'][ss][iobjs,:])
pdict['fluxes'].append(fdict['fluxes'][ss][iobjs])
if objtype == 'ELG':
pdict['OII'].append(fdict['OII'][ss][iobjs])
pdict['exptime'].append(fdict['exptime'][ss])
# Return
return pdict

def load_s2n_values(objtype, nights, channel, sub_exposures=None):
"""
DEPRECATED

Calculate S/N values for a set of spectra

Args:
objtype: str
nights: list
channel: str
sub_exposures:

Returns:
fdict: dict
Contains S/N info

"""
fdict = dict(waves=[], s2n=[], fluxes=[], exptime=[], OII=[])
for night in nights:
if sub_exposures is not None:
Expand All @@ -44,12 +166,14 @@ def load_s2n_values(objtype, nights, channel, sub_exposures=None):
sps_tab = Table(sps_hdu['TRUTH'].data,masked=True)

#- Get OIIFLUX from separate HDU and join
'''
if 'TRUTH_ELG' in sps_hdu:
elg_truth = Table(sps_hdu['TRUTH_ELG'].data)
sps_tab = join(sps_tab, elg_truth['TARGETID', 'OIIFLUX'],
keys='TARGETID', join_type='left')
else:
sps_tab['OIIFLUX'] = 0.0
'''

sps_hdu.close()
objs = sps_tab['TEMPLATETYPE'] == objtype
Expand Down