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

Jump to bottom

Add fix for HASP data sometimes creating a product with no good pixels #47

Merged
merged 3 commits into from
Aug 25, 2023
Merged

Add fix for HASP data sometimes creating a product with no good pixels #47

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d2bd345
Add fix for HASP data sometimes creating a product with no good pixels,
stscirij May 8, 2023
597d42e
Added fix for STIS and CCD data too
stscirij Jun 5, 2023
85f2ac4
Fix a couple of issues that arose from HASP testing:
stscirij Jul 28, 2023
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
87 changes: 55 additions & 32 deletions ullyses/coadd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,8 +191,11 @@ def coadd(self):
weight = flux_weight[goodpixels]
net_counts = segment.data['net'][goodpixels] * segment.exptime
if self.instrument == 'COS':
variances = segment.data['variance_counts'][goodpixels] + segment.data['variance_bkg'][goodpixels] + segment.data['variance_flat'][goodpixels]

try:
variances = segment.data['variance_counts'][goodpixels] + segment.data['variance_bkg'][goodpixels] + segment.data['variance_flat'][goodpixels]
except KeyError:
print('WARNING: Cant get variance keywords for COS data, using gross counts instead')
variances = gcounts
flux = segment.data['flux'][goodpixels]
for i in range(npts):
self.sumnetcounts[indices[i]] = self.sumnetcounts[indices[i]] + net_counts[i]
Expand All @@ -203,30 +206,35 @@ def coadd(self):
if self.instrument == 'COS':
self.output_varsum[indices[i]] = self.output_varsum[indices[i]] + variances[i]
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
nonzeros = np.where(self.output_sumweight != 0)
if self.instrument == 'COS':
# Using the variances (which only COS has) gives spikes in the error when the flux goes negative.
self.output_varsum[nonzeros] = np.where(self.output_varsum[nonzeros] < 0.5, 0.5, self.output_varsum[nonzeros])
self.output_errors[nonzeros] = np.sqrt(self.output_varsum[nonzeros])
if len(good_dq[0]) == 0:
self.first_good_wavelength = None
self.last_good_wavelength = None
print("Warning: No good pixels in product")
else:
# For the moment calculate errors from the gross counts
self.output_errors[nonzeros] = np.sqrt(self.output_sumgcounts[nonzeros])
self.output_flux[nonzeros] = self.output_sumflux[nonzeros] / self.output_sumweight[nonzeros]
# Calculate the conversion from counts to flux units for errors
conversion = self.output_flux[nonzeros] / self.sumnetcounts[nonzeros]
# Clean out NaNs from where flux and net are zero
good = np.where(~np.isnan(conversion))
bad = np.where(np.isnan(conversion))
lenwl = len(wavelength)
interpolated_values = np.interp(self.output_wavelength[bad],
self.output_wavelength[good],
conversion[good])
conversion[bad] = interpolated_values
# Use conversion to calculate error in flux units (erg/cm^2/s/Ang)
self.output_errors[nonzeros] = self.output_errors[nonzeros] * conversion
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
nonzeros = np.where(self.output_sumweight != 0)
if self.instrument == 'COS':
# Using the variances (which only COS has) gives spikes in the error when the flux goes negative.
self.output_varsum[nonzeros] = np.where(self.output_varsum[nonzeros] < 0.5, 0.5, self.output_varsum[nonzeros])
self.output_errors[nonzeros] = np.sqrt(self.output_varsum[nonzeros])
else:
# For the moment calculate errors from the gross counts
self.output_errors[nonzeros] = np.sqrt(self.output_sumgcounts[nonzeros])
self.output_flux[nonzeros] = self.output_sumflux[nonzeros] / self.output_sumweight[nonzeros]
# Calculate the conversion from counts to flux units for errors
conversion = self.output_flux[nonzeros] / self.sumnetcounts[nonzeros]
# Clean out NaNs from where flux and net are zero
good = np.where(~np.isnan(conversion))
bad = np.where(np.isnan(conversion))
lenwl = len(wavelength)
interpolated_values = np.interp(self.output_wavelength[bad],
self.output_wavelength[good],
conversion[good])
conversion[bad] = interpolated_values
# Use conversion to calculate error in flux units (erg/cm^2/s/Ang)
self.output_errors[nonzeros] = self.output_errors[nonzeros] * conversion
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
return

def get_targname(self):
Expand Down Expand Up @@ -339,8 +347,13 @@ def coadd(self):
nonzeros = np.where(self.output_errors != 0.0)
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
if len(good_dq[0]) > 0:
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
else:
self.first_good_wavelength = None
self.last_good_wavelength = None
print('No good data in product')
return

class CCDSegmentList(SegmentList):
Expand Down Expand Up @@ -395,7 +408,7 @@ def create_output_wavelength_grid(self):

self.delta_wavelength = max_delta_wavelength

wavegrid = np.arange(self.min_wavelength, self.max_wavelength, self.delta_wavelength)
wavegrid = np.arange(self.min_wavelength, self.max_wavelength+self.delta_wavelength, self.delta_wavelength)

self.output_wavelength = wavegrid

Expand Down Expand Up @@ -425,8 +438,13 @@ def coadd(self):
nonzeros = np.where(self.output_errors != 0.0)
self.signal_to_noise[nonzeros] = self.output_flux[nonzeros] / self.output_errors[nonzeros]
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
if len(good_dq[0]) > 0:
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
else:
self.first_good_wavelength = None
self.last_good_wavelength = None
print('No good data in product')
else:
nelements = len(self.output_wavelength)
self.output_sumsqerrors = np.zeros(nelements)
Expand All @@ -444,8 +462,13 @@ def coadd(self):
self.output_sumsqerrors[indices[i]] = self.output_sumsqerrors[indices[i]] + (err[i] * weight[i])**2
self.output_exptime[indices[i]] = self.output_exptime[indices[i]] + segment.exptime
good_dq = np.where(self.output_exptime > 0.)
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
if len(good_dq[0]) > 0:
self.first_good_wavelength = self.output_wavelength[good_dq][0]
self.last_good_wavelength = self.output_wavelength[good_dq][-1]
else:
self.first_good_wavelength = None
self.last_good_wavelength = None
print('No good data in product')
nonzeros = np.where(self.output_sumweight != 0)
self.output_flux[nonzeros] = self.output_sumflux[nonzeros] / self.output_sumweight[nonzeros]
self.output_errors[nonzeros] = np.sqrt(self.output_sumsqerrors[nonzeros]) / self.output_sumweight[nonzeros]
Expand Down
Loading