Throughput file updates. Each detector has its own throughput table

romanisim/bandpass.py

@@ -10,6 +10,7 @@
 from scipy import integrate
 from astropy import constants
 from astropy.table import Table
+from astropy.io import ascii
 from astropy import units as u
 import functools
 from romanisim import parameters
@@ -38,11 +39,11 @@
 ABZeroSpFluxDens = 3631e-23 * u.erg / (u.s * u.cm**2 * u.hertz)
 
 
-def read_gsfc_effarea(filename=None):
+def read_gsfc_effarea(sca, filename=None):
     """Read an effective area file from Roman.
 
     This just puts together the right invocation to get an Excel-converted
-    CSV file into memory.
+    ECSV file into memory.
 
     Parameters
     ----------
@@ -55,19 +56,17 @@ def read_gsfc_effarea(filename=None):
         table with effective areas for different Roman bandpasses.
     """
 
-    reader = functools.partial(Table.read, format='csv', delimiter=',',
-                               header_start=1, data_start=2)
+    # The throughput files come as ECSV from March 2024 so the invocation changed to properly read the ECSVs. 
     if filename is None:
-        ref = (resources.files('romanisim') / 'data'
-               / 'Roman_effarea_20201130.csv')
-        with resources.as_file(ref) as filename:
-            out = reader(filename)
+        filename = str(resources.files('romanisim') / 'data' / 'Roman_effarea_tables_20240327' 
+               / f'Roman_effarea_v8_SCA{sca:02}_20240301.ecsv')
+        out = ascii.read(filename)
     else:
-        out = reader(filename)
+        out = ascii.read(filename)
     return out
 
 
-def compute_abflux(effarea=None):
+def compute_abflux(sca, effarea=None):
     """Compute the AB zero point fluxes for each filter.
 
     How many electrons would a zeroth magnitude AB star deposit in
@@ -85,7 +84,7 @@ def compute_abflux(effarea=None):
     """
 
     if effarea is None:
-        effarea = read_gsfc_effarea()
+        effarea = read_gsfc_effarea(sca)
 
     # get the filter names.  This is a bit ugly since there's also
     # a wavelength column 'Wave', and Excel appends a column to each line
@@ -95,11 +94,11 @@ def compute_abflux(effarea=None):
     abfv = ABZeroSpFluxDens
     out = dict()
     for bandpass in filter_names:
-        out[bandpass] = compute_count_rate(flux=abfv, bandpass=bandpass, effarea=effarea)
+        out[bandpass] = compute_count_rate(flux=abfv, bandpass=bandpass, sca=sca, effarea=effarea)
     return out
 
 
-def get_abflux(bandpass):
+def get_abflux(bandpass, sca):
     """Get the zero point flux for a particular bandpass.
 
     This is a simple wrapper for compute_abflux, caching the results.
@@ -116,16 +115,16 @@ def get_abflux(bandpass):
     """
     bandpass = galsim2roman_bandpass.get(bandpass, bandpass)
 
-    # If abflux for this bandpass has been calculated, return the calculated
+    # If abflux for this bandpass for the specified SCA has been calculated, return the calculated
     # value instead of rerunning an expensive calculation
     abflux = getattr(get_abflux, 'abflux', None)
-    if abflux is None:
-        abflux = compute_abflux()
+    if (abflux is None) or (f'SCA{sca:02}' not in abflux):
+        abflux = compute_abflux(sca)
         get_abflux.abflux = abflux
-    return abflux[bandpass]
+    return abflux[f'SCA{sca:02}'][bandpass]
 
 
-def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=None):
+def compute_count_rate(flux, bandpass, sca, filename=None, effarea=None, wavedist=None):
     """Compute the count rate in a given filter, for a specified SED.
 
     How many electrons would an object with SED given by
@@ -151,7 +150,7 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non
     """
     # Read in default Roman effective areas from Goddard, if areas not supplied
     if effarea is None:
-        effarea = read_gsfc_effarea(filename)
+        effarea = read_gsfc_effarea(sca, filename)
 
     # If wavelength distribution is supplied, interpolate flux and area
     # over it and the effective area table layout
@@ -184,7 +183,7 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non
     return zpflux
 
 
-def etomjysr(bandpass):
+def etomjysr(bandpas, sca):
     """Compute factor converting e/s/pix to MJy/sr.
 
     Assumes a pixel scale of 0.11" (romanisim.parameters.pixel_scale)
@@ -200,7 +199,7 @@ def etomjysr(bandpass):
         the factor F such that MJy / sr = F * DN/s
     """
 
-    abflux = get_abflux(bandpass)  # e/s corresponding to 3631 Jy
+    abflux = get_abflux(bandpass, sca)  # e/s corresponding to 3631 Jy
     srpix = ((parameters.pixel_scale * u.arcsec) ** 2).to(u.sr).value
     mjysr = 1 / abflux * 3631 / 10 ** 6 / srpix  # should be ~0.3
     return mjysr