Modify script to apply raft-based adjustments for "questionable" rafts.

python/lsst/obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py

@@ -24,9 +24,118 @@
 import os
 import dateutil.parser
 import numpy as np
+from scipy.interpolate import interp1d
+from scipy.optimize import leastsq
 
 import lsst.utils
 from lsst.meas.algorithms.simple_curve import AmpCurve
+import lsst.afw.math
+
+
+class SplineFitter:
+    """Simple spline fitter to adjust rafts.
+
+    Parameters
+    ----------
+    nodes : `np.ndarray`
+        Node wavelengths.
+    wavelengths : `np.ndarray`
+        Wavelengths (nm).
+    throughput_obs : `np.ndarray`
+        Observed throughput (questionable raft median).
+    throughput_ref : `np.ndarray`
+        Reference throughput (good raft median).
+    """
+    def __init__(self, nodes, wavelengths, throughput_obs, throughput_ref):
+        self._nodes = nodes
+        self._wavelengths = wavelengths
+        self._throughput_obs = throughput_obs
+        self._throughput_ref = throughput_ref
+
+    @staticmethod
+    def compute_ratio_model(nodes, pars, wls, tput_obs, tput_ref, return_spline=False):
+        """Compute the ratio between model and observed.
+
+        Parameters
+        ----------
+        nodes : `np.ndarray`
+            Spline nodes.
+        pars : `np.ndarray`
+            Spline parameters.
+        wls : `np.ndarray`
+            Wavelengths.
+        tput_obs : `np.ndarray`
+            Observed throughput.
+        tput_ref : `np.ndarray`
+            Reference throughput.
+        return_spline : `bool`, optional
+            Return spline interpolation object?
+
+        Returns
+        -------
+        ratio_model : `np.ndarray`
+            Ratio between model and observed.
+        spl : `lsst.afw.math.thing`
+            Spline interpolator (returns if return_spline=True).
+        """
+        spl = lsst.afw.math.makeInterpolate(
+            nodes,
+            pars,
+            lsst.afw.math.stringToInterpStyle("AKIMA_SPLINE"),
+        )
+
+        model = spl.interpolate(wls)
+        ratio = (tput_obs * model) / tput_ref
+
+        if return_spline:
+            return ratio, spl
+        else:
+            return ratio
+
+    def fit(self, p0):
+        """Fit the spline function.
+
+        Parameters
+        ----------
+        p0 : `np.ndarray`
+            Array of starting parameters.
+
+        Returns
+        -------
+        pars : `np.ndarray`
+            Best fit spline parameters.
+        """
+        params, cov_params, _, msg, ierr = leastsq(
+            self,
+            p0,
+            full_output=True,
+            ftol=1e-5,
+            maxfev=12000,
+        )
+
+        return params
+
+    def __call__(self, pars):
+        """Compute the residuals for leastsq.
+
+        Parameters
+        ----------
+        pars : `np.ndarray`
+            Spline parameters.
+
+        Returns
+        -------
+        residuals : `np.ndarray`
+            Fit residuals.
+        """
+        ratio_model = self.compute_ratio_model(
+            self._nodes,
+            pars,
+            self._wavelengths,
+            self._throughput_obs,
+            self._throughput_ref,
+        )
+        return ratio_model - 1.0
 
 
 data_path = lsst.utils.getPackageDir("obs_lsst_data")
@@ -39,18 +148,133 @@
 
 data = Table.read(parquet_file)
 
+# Code to do adjustments of questionable rafts.
+questionable_rafts = ["R03", "R11", "R21", "R32", "R42"]
+
+e2v_rafts = ["R11", "R12", "R13", "R14",
+             "R21", "R22", "R23", "R24",
+             "R30", "R31", "R32", "R33", "R34"]
+itl_rafts = ["R01", "R02", "R03",
+             "R10",
+             "R20",
+             "R41", "R42", "R43"]
+
+n_amp_per_det = 16
+n_det_per_raft = 9
+
+# Nodes chosen to cover the wavelength range of the QE curve data.
+nodes = np.linspace(320.0, 1099.0, 20)
+
+# We will do all fitting at a standardized set of wavelengths.
+wavelengths = np.linspace(np.min(nodes), np.max(nodes), 1000)
+
+for det_type in ["e2v", "itl"]:
+    if det_type == "e2v":
+        rafts = e2v_rafts
+    else:
+        rafts = itl_rafts
+
+    tput_amps = np.zeros((len(wavelengths), n_amp_per_det*n_det_per_raft*len(rafts)))
+    tput_amps[:, :] = np.nan
+
+    counter = 0
+
+    for raft in rafts:
+        if raft in questionable_rafts:
+            continue
+
+        raft_use, = np.where((data["bay"] == raft) & (data["seg"] != "Ave"))
+
+        det_amps = []
+        for row in data[raft_use]:
+            det_amps.append(row["bay"] + row["slot"] + row["seg"])
+        det_amps = np.array(det_amps)
+        unique_det_amps = np.unique(det_amps)
+
+        for i, det_amp in enumerate(unique_det_amps):
+            amp_use, = np.where(det_amps == det_amp)
+
+            interp = interp1d(
+                data["wl"][raft_use][amp_use],
+                data["qecorr"][raft_use][amp_use],
+                bounds_error=False,
+                fill_value=0.0,
+            )
+            tput_amps[:, counter] = interp(wavelengths)
+
+            counter += 1
+
+    if det_type == "e2v":
+        tput_e2v_median = np.nanmedian(tput_amps, axis=1)
+    else:
+        tput_itl_median = np.nanmedian(tput_amps, axis=1)
+
+# Compute the median for each questionable raft.
+questionable_throughputs = {}
+
+for raft in questionable_rafts:
+    counter = 0
+
+    tput_amps = np.zeros((len(wavelengths), n_amp_per_det*n_det_per_raft))
+
+    raft_use, = np.where((data["bay"] == raft) & (data["seg"] != "Ave"))
+
+    det_amps = []
+    for row in data[raft_use]:
+        det_amps.append(row["bay"] + row["slot"] + row["seg"])
+    det_amps = np.array(det_amps)
+    unique_det_amps = np.unique(det_amps)
+
+    for i, det_amp in enumerate(unique_det_amps):
+        amp_use, = np.where(det_amps == det_amp)
+
+        interp = interp1d(
+            data["wl"][raft_use][amp_use],
+            data["qecorr"][raft_use][amp_use],
+            bounds_error=False,
+            fill_value=0.0,
+        )
+        tput_amps[:, counter] = interp(wavelengths)
+
+        counter += 1
+
+    questionable_throughputs[raft] = np.nanmedian(tput_amps, axis=1)
+
+# For each questionable raft, we want to fit some spline nodes.
+questionable_spline_correctors = {}
+for raft in questionable_rafts:
+
+    if raft in e2v_rafts:
+        throughput_ref = tput_e2v_median
+    else:
+        throughput_ref = tput_itl_median
+
+    fitter = SplineFitter(nodes, wavelengths, questionable_throughputs[raft], throughput_ref)
+    pars = fitter.fit(np.ones(len(nodes)))
+
+    _, spl = fitter.compute_ratio_model(nodes, pars, wavelengths, questionable_throughputs[raft], throughput_ref, return_spline=True)
+
+    questionable_spline_correctors[raft] = spl
+
 det_nums = np.unique(data["idet"])
 
 for det_num in det_nums:
     det_use, = np.where((data["idet"] == det_num) & (data["seg"] != "Ave"))
     slot = data["slot"][det_use[0]]
     bay = data["bay"][det_use[0]]
 
+    wavelength = np.array(data["wl"][det_use])
+    efficiency = np.array(data["qecorr"][det_use])
+    if bay in questionable_rafts:
+        # Fix this up with spline.
+        spl = questionable_spline_correctors[bay]
+        efficiency *= spl.interpolate(wavelength)
+
     curve_table = QTable(
         {
             "amp_name": np.array(data["seg"][det_use]),
-            "wavelength": np.array(data["wl"][det_use]) * u.nanometer,
-            "efficiency": np.array(data["qecorr"][det_use]) * u.percent,
+            "wavelength": wavelength * u.nanometer,
+            "efficiency": efficiency * u.percent,
         }
     )
     curve = AmpCurve.fromTable(curve_table)