diff --git a/python/lsst/obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py b/python/lsst/obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py new file mode 100644 index 000000000..aee8fc60e --- /dev/null +++ b/python/lsst/obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py @@ -0,0 +1,359 @@ +# This file is part of obs_lsst. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# This script is a record of a one-time run for DM-40164. +# See obs_lsst_data lsstCam/transmission_sensor/README.md +# for details. + + +import astropy.units as u +from astropy.table import Table, QTable +import re +import os +import dateutil.parser +import numpy as np +from scipy.interpolate import interp1d +from scipy.optimize import leastsq +import copy + +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 + + +debug_display = False + +data_path = lsst.utils.getPackageDir("obs_lsst_data") +transmission_path = os.path.join(data_path, "lsstCam", "transmission_sensor") +parquet_file = os.path.join(transmission_path, "qe_raft_allvalues_nircorrected_20230725.parquet") +parquet_file_update = os.path.join( + transmission_path, + "qe_raft_allvalues_nircorrected_20230725_adjust.parquet", +) + +valid_start = "1970-01-01T00:00:00" +valid_date = dateutil.parser.parse(valid_start) +datestr = ''.join(re.split(r'[:-]', valid_date.isoformat())) + +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. +# The number of nodes is chosen to be large enough to make the +# corrections "well-matched" to the template. +nodes = np.linspace(320.0, 1099.0, 40) + +# 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, + ) + + if debug_display: + import matplotlib.pyplot as plt + + corrected = np.array(spl.interpolate(wavelengths) * questionable_throughputs[raft]) + + plt.clf() + plt.plot(wavelengths, questionable_throughputs[raft], 'r-', label="Questionable") + plt.plot(wavelengths, throughput_ref, "b-", label="Reference") + plt.plot(wavelengths, corrected, "m-", label="Corrected") + plt.legend() + plt.title(f"Raft = {raft}") + plt.show() + + questionable_spline_correctors[raft] = spl + +det_nums = np.unique(data["idet"]) + +data_update = copy.copy(data) + +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) + + data_update["qecorr"][det_use][:] = efficiency + + curve_table = QTable( + { + "amp_name": np.array(data["seg"][det_use]), + "wavelength": wavelength * u.nanometer, + "efficiency": efficiency * u.percent, + } + ) + curve = AmpCurve.fromTable(curve_table) + + out_path = os.path.join(transmission_path, bay.lower() + "_" + slot.lower()) + os.makedirs(out_path, exist_ok=True) + + out_file = os.path.join(out_path, datestr + ".ecsv") + + curve_table.meta.update( + { + "CALIBDATE": valid_start, + "INSTRUME": "LSSTCAM", + "OBSTYPE": "transmission_sensor", + "TYPE": "transmission_sensor", + "DETECTOR": det_num, + "PARQUETFILE": os.path.basename(parquet_file), + "CALIBCLS": "lsst.ip.isr.IntermediateSensorTransmissionCurve", + } + ) + curve_table.meta["CALIB_ID"] = ( + f"raftName={bay} detectorName={slot} " + f"detector={det_num} calibDate={valid_start} " + f"ccd={det_num} ccdnum={det_num} filter=None" + ) + + # We need to remove any previous file if it is there. + if os.path.isfile(out_file): + os.remove(out_file) + + curve.writeText(out_file) + + # And update the average. + det_use_ave, = np.where((data["idet"] == det_num) & (data["seg"] == "Ave")) + wavelength_ave = np.array(data["wl"][det_use_ave]) + efficiency_ave = np.array(data["qecorr"][det_use_ave]) + if bay in questionable_rafts: + spl = questionable_spline_correctors[bay] + efficiency_ave *= spl.interpolate(wavelength_ave) + + data_update["qecorr"][det_use_ave][:] = efficiency_ave + +data_update.write(parquet_file_update, overwrite=True)