From 8a0cdb465d9c401a249645a769f9921eb9c41a7c Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Wed, 2 Aug 2023 09:58:50 -0700 Subject: [PATCH 1/5] Add script to convert parquet QE curves for LSSTCam to ecsv files. --- .../rewrite_lsstcam_qe_files_DM-40164.py | 84 +++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 python/lsst/obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py 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..b9b81c6d1 --- /dev/null +++ b/python/lsst/obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py @@ -0,0 +1,84 @@ +# 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 . +import astropy.units as u +from astropy.table import Table, QTable +import re +import os +import dateutil.parser +import numpy as np + +import lsst.utils +from lsst.meas.algorithms.simple_curve import AmpCurve + + +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") + +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) + +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]] + + 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, + } + ) + 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) From 66acc116b771bf5a77a4383c89f78537be67cbc3 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Mon, 21 Aug 2023 19:29:52 -0700 Subject: [PATCH 2/5] Modify script to apply raft-based adjustments for "questionable" rafts. --- .../rewrite_lsstcam_qe_files_DM-40164.py | 228 +++++++++++++++++- 1 file changed, 226 insertions(+), 2 deletions(-) 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 index b9b81c6d1..53b1aa1ee 100644 --- 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 @@ -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,6 +148,114 @@ 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: @@ -46,11 +263,18 @@ 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) From 6cf10989cfcc59a41dffe0949b896468e060b2ca Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Tue, 22 Aug 2023 09:43:44 -0700 Subject: [PATCH 3/5] Add output of adjusted parquet file. --- .../rewrite_lsstcam_qe_files_DM-40164.py | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) 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 index 53b1aa1ee..127c5d68d 100644 --- 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 @@ -26,6 +26,7 @@ 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 @@ -141,6 +142,10 @@ def __call__(self, pars): 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) @@ -252,12 +257,21 @@ def __call__(self, pars): 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) + _, 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"]) +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]] @@ -270,6 +284,8 @@ def __call__(self, pars): 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]), @@ -306,3 +322,15 @@ def __call__(self, pars): 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) From 3553b6785929781fbb84a12af70eed8f6ab8098c Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 18 Jan 2024 14:42:58 -0800 Subject: [PATCH 4/5] Increase number of nodes and add debug display feature. --- .../rewrite_lsstcam_qe_files_DM-40164.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) 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 index 127c5d68d..f71a20550 100644 --- 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 @@ -139,6 +139,8 @@ def __call__(self, pars): 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") @@ -168,7 +170,9 @@ def __call__(self, pars): 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) +# 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) @@ -266,6 +270,19 @@ def __call__(self, pars): 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"]) From 6e172191d517211474ef40ba78b5d91e9e384100 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Wed, 31 Jan 2024 13:08:22 -0800 Subject: [PATCH 5/5] Add comment to top of one-off script. --- .../obs/lsst/script/rewrite_lsstcam_qe_files_DM-40164.py | 6 ++++++ 1 file changed, 6 insertions(+) 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 index f71a20550..aee8fc60e 100644 --- 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 @@ -18,6 +18,12 @@ # # 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