Skip to content

Commit

Permalink
Merge pull request #172 from spacetelescope/calcos_3.3.10_rc1
Browse files Browse the repository at this point in the history 
Calcos 3.3.10 (#171)
  • Loading branch information
@rendinam
rendinam committed Oct 12, 2020
2 parents 0c0f43c + 795738f commit 1731074
Show file tree
Hide file tree
Showing 5 changed files with 168 additions and 49 deletions.
3 changes: 2 additions & 1 deletion calcos/concurrent.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -602,7 +602,8 @@ def findShifts(self):
snr_ff = 0. # ignore error array
axis = 1 # dispersion is along X axis
dummy_exptime = 1.
(N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i,
(N_i, ERR_i, ERR_LOW_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i,
DQ_ALL_i, LOWER_OUTER_INDEX_i, UPPER_OUTER_INDEX_i,
LOWER_INNER_INDEX_i, UPPER_INNER_INDEX_i,
ENCLOSED_FRACTION_i, AV_E_BKG_i,
Expand Down
37 changes: 32 additions & 5 deletions calcos/cosutil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
import numpy as np
import numpy.linalg as LA
import astropy.io.fits as fits
from astropy.stats import poisson_conf_interval
from . import ccos
from .calcosparam import * # parameter definitions

Expand Down Expand Up @@ -3733,24 +3734,50 @@ def centerOfQuartic(x, coeff):

return x_min

def errGehrels(a):
def errGehrels(counts):
"""Compute error estimate.
The error estimate is computed using the Gehrels approximation for the
upper confidence limit.
Parameters
----------
a: array_like or float
counts: array_like or float
Number of counts (not necessarily integer values).
Returns
-------
array_like or float
The error estimate for a counts.
tuple of 2 array_like or float
(The lower error estimate for counts,
the upper error estimate for counts)
"""
icounts = (counts + .5).astype(np.int)
upper = (1. + np.sqrt(icounts + 0.75))
lower = np.where(icounts > 0., Gehrels_lower(icounts), 0.)
return (lower.astype(np.float32), upper.astype(np.float32))

return (1. + np.sqrt(a + 0.75))
def Gehrels_lower(counts):
return counts - counts * (1.0 - 1.0 / (9.0 * counts) - 1.0 / (3.0 * np.sqrt(counts)))**3

def errFrequentist(counts):
"""Compute errors using the 'frequentist-confidence' option of astropy's poisson_conf_interval
Parameters
----------
counts: array-like or float
Number of counts (not necessarily integer values).
Returns
-------
tuple of 2 array-like or float
(The lower error estimate for counts,
the upper error estimate for counts)
"""

lower, upper = poisson_conf_interval(counts, interval='frequentist-confidence')
err_lower = counts - lower
err_upper = upper - counts
return (err_lower.astype(np.float32), err_upper.astype(np.float32))

def precess(t, target):
"""Precess target to the time of observation.
Expand Down
121 changes: 83 additions & 38 deletions calcos/extract.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import os
import numpy as np
import astropy.io.fits as fits
from astropy.stats import poisson_conf_interval
from . import cosutil
from . import ccos
from . import dispersion
Expand Down Expand Up @@ -163,6 +164,11 @@ def extract1D(input, incounts=None, output=None,
unit="erg /s /cm**2 /angstrom"))
col.append(fits.Column(name="ERROR", format=rpt+"E",
unit="erg /s /cm**2 /angstrom"))
col.append(fits.Column(name="ERROR_LOWER", format=rpt+"E",
unit="erg /s /cm**2 /angstrom"))
col.append(fits.Column(name="VARIANCE_FLAT", format=rpt+"E"))
col.append(fits.Column(name="VARIANCE_COUNTS", format=rpt+"E"))
col.append(fits.Column(name="VARIANCE_BKG", format=rpt+"E"))
col.append(fits.Column(name="GROSS", format=rpt+"E",
unit="count /s"))
col.append(fits.Column(name="GCOUNTS", format=rpt+"E",
Expand Down Expand Up @@ -601,7 +607,8 @@ def doExtract(ifd_e, ifd_c, ofd, nelem,
if is_corrtag:
key = "shift1" + segment[-1]
shift1 = ofd[1].header.get(key, 0.)
(N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i, DQ_ALL_i,
(N_i, ERROR_i, ERROR_LOWER_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i, DQ_ALL_i,
LOWER_OUTER_i, UPPER_OUTER_i, LOWER_INNER_i, UPPER_INNER_i,
ENCLOSED_FRACTION_i, BACKGROUND_PER_ROW_i, EE_LOWER_OUTER_i,
EE_LOWER_INNER_i, EE_UPPER_INNER_i, EE_UPPER_OUTER_i
Expand All @@ -615,7 +622,9 @@ def doExtract(ifd_e, ifd_c, ofd, nelem,
local_find_targ)
else:
if info["xtrctalg"] == 'BOXCAR':
(N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i, DQ_ALL_i,
(N_i, ERROR_i, ERROR_LOWER_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i,
DQ_i, DQ_WGT_i, DQ_ALL_i,
LOWER_OUTER_i, UPPER_OUTER_i, LOWER_INNER_i, UPPER_INNER_i,
ENCLOSED_FRACTION_i, BACKGROUND_PER_ROW_i, EE_LOWER_OUTER_i,
EE_LOWER_INNER_i, EE_UPPER_INNER_i, EE_UPPER_OUTER_i
Expand All @@ -630,7 +639,9 @@ def doExtract(ifd_e, ifd_c, ofd, nelem,
user_xdisp_locn, user_xdisp_size,
local_find_targ)
elif info["xtrctalg"] == 'TWOZONE':
(N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i, DQ_ALL_i,
(N_i, ERROR_i, ERROR_LOWER_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i, DQ_i,
DQ_WGT_i, DQ_ALL_i,
LOWER_OUTER_i, UPPER_OUTER_i, LOWER_INNER_i, UPPER_INNER_i,
ENCLOSED_FRACTION_i, BACKGROUND_PER_ROW_i, EE_LOWER_OUTER_i,
EE_LOWER_INNER_i, EE_UPPER_INNER_i, EE_UPPER_OUTER_i
Expand All @@ -648,7 +659,9 @@ def doExtract(ifd_e, ifd_c, ofd, nelem,
else:
cosutil.printMsg("Unknown extraction method, defaulting to", \
" BOXCAR")
(N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i, DQ_ALL_i,
(N_i, ERROR_i, ERROR_LOWER_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i, DQ_i,
DQ_WGT_i, DQ_ALL_i,
LOWER_OUTER_i, UPPER_OUTER_i, LOWER_INNER_i, UPPER_INNER_i,
ENCLOSED_FRACTION_i, BACKGROUND_PER_ROW_i, EE_LOWER_OUTER_i,
EE_LOWER_INNER_i, EE_UPPER_INNER_i, EE_UPPER_OUTER_i
Expand All @@ -668,9 +681,13 @@ def doExtract(ifd_e, ifd_c, ofd, nelem,
outdata.field("EXPTIME")[row] = exptime
outdata.field("WAVELENGTH")[row][:] = wavelength.copy()
outdata.field("FLUX")[row][:] = 0.
outdata.field("ERROR")[row][:] = ERR_i.copy()
outdata.field("ERROR")[row][:] = ERROR_i.copy()
outdata.field("ERROR_LOWER")[row][:] = ERROR_LOWER_i.copy()
outdata.field("VARIANCE_FLAT")[row][:] = VARIANCE_FLAT_i.copy()
outdata.field("VARIANCE_BKG")[row][:] = VARIANCE_BKG_i.copy()
outdata.field("GROSS")[row][:] = GC_i.copy()
outdata.field("GCOUNTS")[row][:] = GCOUNTS_i.copy()
outdata.field("VARIANCE_COUNTS")[row][:] = VARIANCE_COUNTS_i.copy()
outdata.field("NET")[row][:] = N_i.copy()
# outdata.field("NET_ERROR")[row][:] = ERR_i.copy() xxx
outdata.field("BACKGROUND")[row][:] = BK_i.copy()
Expand Down Expand Up @@ -1128,20 +1145,28 @@ def extractSegmentBoxcar(e_data, c_data, e_dq_data, ofd_header, segment,
del temp_bk
else:
BK_i = np.zeros(axis_length, dtype=np.float32)

# The error in the counts is the sum in quadrature of 3 terms
N_i = eps_i * (GC_i - BK_i)

if snr_ff > 0.:
term1_i = (N_i * exptime / (extr_height * snr_ff))**2
VARIANCE_FLAT_i = (N_i * exptime / (extr_height * snr_ff))**2
else:
term1_i = 0.
term2_i = eps_i**2 * exptime * \
(GC_i + BK_i * (bkg_norm / float(bkg_smooth)))
VARIANCE_FLAT_i = N_i * 0.
VARIANCE_COUNTS_i = eps_i**2 * GC_i * exptime
temp_val = BK_i * exptime * (bkg_norm / float(bkg_smooth))
VARIANCE_BKG_i = eps_i * eps_i * temp_val

variance_i = VARIANCE_FLAT_i + VARIANCE_COUNTS_i + VARIANCE_BKG_i
# Use the frequentist option of the astropy Poisson confidence interval function
# to calculate errors
ERROR_LOWER_i, ERROR_i = cosutil.errFrequentist(variance_i)
# ERR_i is the error in the count RATE
if exptime > 0.:
# Use the Gehrels variance function.
ERR_i = cosutil.errGehrels(term1_i + term2_i) / exptime
ERROR_LOWER_i /= exptime
ERROR_i /= exptime
else:
ERR_i = N_i * 0.
ERROR_LOWER_i = N_i * 0.
ERROR_i = N_i * 0.

updateExtractionKeywords(ofd_header, segment,
slope, extr_height,
Expand All @@ -1165,7 +1190,8 @@ def extractSegmentBoxcar(e_data, c_data, e_dq_data, ofd_header, segment,
UPPER_INNER_VALUE_i = N_i*0.0 + 1.0
UPPER_OUTER_VALUE_i = N_i*0.0 + 1.0

return (N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i,
return (N_i, ERROR_i, ERROR_LOWER_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i,
DQ_ALL_i, LOWER_OUTER_INDEX_i, UPPER_OUTER_INDEX_i,
LOWER_INNER_INDEX_i, UPPER_INNER_INDEX_i,
ENCLOSED_FRACTION_i, AV_E_BKG_i,
Expand Down Expand Up @@ -1540,26 +1566,35 @@ def extractSegmentTwozone(e_data, c_data, e_dq_data, ofd_header, segment,
DQ_ALL_i[column] = DQ_PIXEL_OUT_OF_BOUNDS
DQ_WGT_i[column] = 0.0
N_i = total_ecounts
goodcolumns = np.where(nrows_c_bkg_i > 0)
DQ_WGT_i = np.where(nrows_c_bkg_i > 0, DQ_WGT_i, 0.0)
flat_correction = total_ecounts / np.where(total_ccounts <= 0.0, 1.0,
total_ccounts)
flat_correction = np.where(flat_correction == 0.0, 1.0, flat_correction)
#
# Now calculate the error
VARIANCE_FLAT_i = np.zeros(ncols, dtype=np.float32)
if snr_ff > 0.0:
term1_i = (N_i * exptime / (extr_height_i * snr_ff))**2
VARIANCE_FLAT_i[goodcolumns] = (N_i[goodcolumns] * exptime / (extr_height_i[goodcolumns] * snr_ff))**2
else:
term1_i = 0.0
term2_i = np.zeros(ncols, dtype=np.float32)
goodcolumns = np.where(nrows_c_bkg_i > 0)
term2_i[goodcolumns] = (flat_correction[goodcolumns])**2 * exptime \
* (total_ccounts[goodcolumns] + \
av_c_bkg_i[goodcolumns] * (extr_height_i[goodcolumns])**2 \
/ (nrows_c_bkg_i[goodcolumns] * bkg_smooth))
VARIANCE_FLAT_i = N_i * 0.0
VARIANCE_COUNTS_i = np.zeros(ncols, dtype=np.float32)
VARIANCE_BKG_i = np.zeros(ncols, dtype=np.float32)
VARIANCE_COUNTS_i[goodcolumns] = (flat_correction[goodcolumns])**2 * exptime \
* total_ccounts[goodcolumns]
VARIANCE_BKG_i[goodcolumns] = (flat_correction[goodcolumns])**2 * exptime \
* av_c_bkg_i[goodcolumns] * (extr_height_i[goodcolumns])**2 \
/ (nrows_c_bkg_i[goodcolumns] * bkg_smooth)
if exptime > 0.0:
# Use the Gehrels variance function
ERR_i = cosutil.errGehrels(term1_i + term2_i) / exptime
# Use the frequentist option of the astropy Poisson confidence interval function
# to calculate errors
VARIANCE_i = VARIANCE_FLAT_i + VARIANCE_COUNTS_i + VARIANCE_BKG_i
ERROR_LOWER_i, ERROR_i = cosutil.errFrequentist(VARIANCE_i)
ERROR_LOWER_i /= exptime
ERROR_i /= exptime
else:
ERR_i = N_i * 0.0
ERROR_LOWER_i = N_i * 0.0
ERROR_i = N_i * 0.0
SUMMED_BACKGROUND_i = AV_E_BKG_i * extr_height_i
GC_i = total_ccounts
GCOUNTS_i = GC_i * exptime
Expand Down Expand Up @@ -1593,12 +1628,14 @@ def extractSegmentTwozone(e_data, c_data, e_dq_data, ofd_header, segment,
xd_nominal, centroid, cent_err, offset,
b_bkg1, b_bkg2,
bkg_height1, bkg_height2)
return (N_i, ERR_i, GC_i, GCOUNTS_i, SUMMED_BACKGROUND_i, DQ_i, DQ_WGT_i,
return (N_i, ERROR_i, ERROR_LOWER_i, VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, SUMMED_BACKGROUND_i, DQ_i, DQ_WGT_i,
DQ_ALL_i, LOWER_OUTER_INDEX_i, UPPER_OUTER_INDEX_i,
LOWER_INNER_INDEX_i, UPPER_INNER_INDEX_i,
ENCLOSED_FRACTION_i, AV_E_BKG_i,
LOWER_OUTER_VALUE_i, LOWER_INNER_VALUE_i,
UPPER_INNER_VALUE_i, UPPER_OUTER_VALUE_i)
UPPER_INNER_VALUE_i, UPPER_OUTER_VALUE_i
)

def getProfileCentroid(phdr, segment):
key = "SP_LOC_" + segment[-1]
Expand Down Expand Up @@ -2070,18 +2107,22 @@ def extractCorrtag(xi, eta, dq, epsilon, dq_array,
N_i = eps_i * (GC_i - BK_i)

if snr_ff > 0.:
term1_i = (N_i * exptime / (extr_height * snr_ff))**2
VARIANCE_FLAT_i = (N_i * exptime / (extr_height * snr_ff))**2
else:
term1_i = 0.
term2_i = eps_i**2 * exptime * \
(GC_i + BK_i * (bkg_norm / float(bkg_smooth)))
VARIANCE_FLAT_i = N_i * 0.
VARIANCE_COUNTS_i = eps_i**2 * exptime * GC_i
VARIANCE_BKG_i = eps_i**2 * exptime * BK_i * (bkg_norm / float(bkg_smooth))
if exptime > 0.:
sum_terms = term1_i + term2_i
sum_terms = np.where(sum_terms > 0, sum_terms, 0.)
# Use the Gehrels variance function.
ERR_i = cosutil.errGehrels(sum_terms) / exptime
VARIANCE_i = VARIANCE_FLAT_i + VARIANCE_COUNTS_i + VARIANCE_BKG_i
VARIANCE_i = np.where(VARIANCE_i > 0, VARIANCE_i, 0.)
# Use the frequentist option of astropy Poisson Confidence Interval function
# to calculate errors
ERROR_LOWER_i, ERROR_i = cosutil.errFrequentist(VARIANCE_i)
ERROR_LOWER_i /= exptime
ERROR_i /= exptime
else:
ERR_i = N_i * 0.
ERROR_LOWER_i = N_i * 0.
ERROR_i = N_i * 0.
if ofd_header is not None:
xd_offset = -999. # not implemented yet
updateExtractionKeywords(ofd_header, segment,
Expand All @@ -2100,7 +2141,9 @@ def extractCorrtag(xi, eta, dq, epsilon, dq_array,
UPPER_INNER_VALUE_i = N_i*0.0 + 1.0
UPPER_OUTER_VALUE_i = N_i*0.0 + 1.0

return (N_i, ERR_i, GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i,
return (N_i, ERROR_i, ERROR_LOWER_i,
VARIANCE_FLAT_i, VARIANCE_COUNTS_i, VARIANCE_BKG_i,
GC_i, GCOUNTS_i, BK_i, DQ_i, DQ_WGT_i,
DQ_ALL_i, LOWER_OUTER_INDEX_i, UPPER_OUTER_INDEX_i,
LOWER_INNER_INDEX_i, UPPER_INNER_INDEX_i,
ENCLOSED_FRACTION_i, AV_E_BKG_i,
Expand Down Expand Up @@ -2139,7 +2182,7 @@ def doFluxCorr(ofd, info, reffiles, tdscorr):
net = outdata.field("NET")
flux = outdata.field("FLUX")
error = outdata.field("ERROR")

error_lower = outdata.field("ERROR_LOWER")
fluxtab = reffiles["fluxtab"]

# segment will be added to filter in the loop
Expand All @@ -2165,6 +2208,7 @@ def doFluxCorr(ofd, info, reffiles, tdscorr):
factor = np.where(factor <= 0., 1., factor)
flux[row][:] = net[row] / factor
error[row][:] = error[row] / factor
error_lower[row][:] = error_lower[row] / factor
ofd[0].header["fluxcorr"] = "COMPLETE"

# Compute an array of time-dependent correction factors (a potentially
Expand Down Expand Up @@ -2216,6 +2260,7 @@ def doFluxCorr(ofd, info, reffiles, tdscorr):
ccos.interp1d(wl_tds, factor_tds, wavelength[row], factor)
flux[row][:] /= factor
error[row][:] /= factor
error_lower[row][:] /= factor
if extrapolate and not printed:
cosutil.printWarning("TDS correction was extrapolated.")
printed = True
Expand Down

0 comments on commit 1731074

Please sign in to comment.