Fixing sensitivity bumps

python/lvmdrp/core/fluxcal.py

@@ -172,9 +172,9 @@ def butter_lowpass_filter(data, cutoff_freq, nyq_freq, order=4):
     return y
 
 
-def filter_channel(w, f, k=3):
+def filter_channel(w, f, k=3, cutoff_freq=0.01):
     c = np.where(np.isfinite(f))
-    s = butter_lowpass_filter(f[c], 0.01, 2)
+    s = butter_lowpass_filter(f[c], cutoff_freq, 2)
     res = s - f[c]
     # plt.plot(w[c], f[c], 'k.')
     # plt.plot(w[c], s, 'b-')

python/lvmdrp/functions/fluxCalMethod.py

@@ -222,13 +222,15 @@ def apply_fluxcal(in_rss: str, out_fframe: str, method: str = 'STD', display_plo
     return fframe
 
 
-def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width=3):
+def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width=3, cutoff_freq=0.01):
     # load the sky masks
     channel = rss._header['CCD']
     w = rss._wave
 
     m = get_sky_mask_uves(w, width=width)
     m2 = None
+    # if channel == "r":
+    #     m2 = ~(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))
     if channel == "z":
         m2 = get_z_continuum_mask(w)
 
@@ -262,6 +264,16 @@ def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width=
 
         spec = (rss._data[fibidx[0],:] - master_sky._data[fibidx[0],:])/exptime  #TODO: check exptime?
 
+        # remove nans
+        # spec[(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))] = np.nan
+        # plt.figure()
+        # plt.plot(w, spec)
+
+        spec = fluxcal.interpolate_mask(w, spec, np.isnan(spec), fill_value="extrapolate")
+
+        # print(">>>>>>>>>>", w[(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))])
+
+
         # interpolate over bright sky lines
         spec = fluxcal.interpolate_mask(w, spec, m, fill_value="extrapolate")
         if channel == "z":
@@ -277,7 +289,13 @@ def standard_sensitivity(stds, rss, GAIA_CACHE_DIR, ext, res, plot=False, width=
 
         # divide to find sensitivity and smooth
         sens = stdflux / spec
-        wgood, sgood = fluxcal.filter_channel(w, sens, 2)
+        if channel == "r":
+            sens[(((w>6563-80)&(w<6563+80))|((w>6850)&(w<6950)))] = np.nan
+        if channel == "b":
+            sens[(((w>4860-50)&(w<4860+50))|((w>4340-50)&(w<4340+30)))] = np.nan
+        wgood, sgood = fluxcal.filter_channel(w, sens, 3, cutoff_freq=cutoff_freq)
+        # wgood = w[np.isfinite(sens)]
+        # sgood = sens[np.isfinite(sens)]
         s = interpolate.make_smoothing_spline(wgood, sgood, lam=1e4)
         res[f"STD{nn}SEN"] = s(w).astype(np.float32) * u.Unit("erg / (ct cm2)")
 

python/lvmdrp/functions/imageMethod.py

@@ -2975,52 +2975,84 @@ def reprojectRSS_drp(
     rep.writeto(f"{out_path}/{out_name}_2d.fits", overwrite=True)
 
 
-def testres_drp(image, trace, fwhm, flux):
-    """
-    Historic task used for debugging of the the extraction routine...
+def validate_extraction(in_image, in_cent, in_width, in_rss, plot_columns=[1000, 2000, 3000], display_plots=False):
+    """Evaluates the extracted flux into the original 2D pixel grid
+
+    This routine will evaluate the extracted flux in the original
+    2D grid and compare the resulting 2D model against the original
+    2D image. Three images are stored as outputs:
+
+        - The residual 2D image: model - data
+        - The ratio 2D image: model / data
+        - The 2D model
+
+    Parameters
+    ----------
+    in_image : str
+        Path to the original 2D image of the extracted flux
+    in_cent : str
+        Path to the fiber centroids trace
+    in_width : str
+        Path to the fiber width (FWHM) trace
+    in_rss : str
+        Path to the extracted flux in RSS format
+    plot_columns : array_like, optional
+        columns to show in plot, by default [1000, 2000, 3000]
+    display_plots : bool, optional
+        whether to display plots to screen or not, by dafult False
     """
+    log.info(f"loading 2D image {in_image}")
     img = Image()
-    # t1 = time.time()
-    img.loadFitsData(image, extension_data=0)
-    trace_mask = TraceMask()
-    trace_mask.loadFitsData(trace, extension_data=0)
-    trace_fwhm = TraceMask()
-    #   trace_fwhm.setData(data=numpy.ones(trace_mask._data.shape)*2.5)
-    trace_fwhm.loadFitsData(fwhm, extension_data=0)
+    img._data = numpy.nan_to_num(img._data)
+    img.loadFitsData(in_image)
 
-    trace_flux = TraceMask()
-    trace_flux.loadFitsData(flux, extension_data=0)
+    log.info(f"loading fiber parameters in {in_cent} and {in_width}")
+    cent = TraceMask.from_file(in_cent)
+    width = TraceMask.from_file(in_width)
+    width._data /= 2.354
+
+    log.info(f"loading extracted flux in {in_rss}")
+    rss = RSS.from_file(in_rss)
+    rss._data = numpy.nan_to_num(rss._data)
+
+    ypix_cor = rss._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"]
+    ypix_ori = img._slitmap[["spectrographid"] == int(img._header["CCD"][1])]["ypix_z"]
+    thermal_shift = ypix_cor - ypix_ori
+    log.info(f"fiber thermal shift in slitmap: {thermal_shift:.4f}")
+    cent._data += thermal_shift
+
+    log.info(f"evaluating extracted flux into 2D pixel grid for {img._dim[1]} columns")
     x = numpy.arange(img._dim[0])
     out = numpy.zeros(img._dim)
     fact = numpy.sqrt(2.0 * numpy.pi)
+
+    fig, axs = create_subplots(to_display=display_plots, nrows=len(plot_columns), ncols=1, figsize=(15,5), sharex=True, layout="constrained")
     for i in range(img._dim[1]):
-        #  print i
-        A = (
-            1.0
-            * numpy.exp(
-                -0.5
-                * (
-                    (x[:, numpy.newaxis] - trace_mask._data[:, i][numpy.newaxis, :])
-                    / abs(trace_fwhm._data[:, i][numpy.newaxis, :] / 2.354)
-                )
-                ** 2
-            )
-            / (fact * abs(trace_fwhm._data[:, i][numpy.newaxis, :] / 2.354))
-        )
-        spec = numpy.dot(A, trace_flux._data[:, i])
+        A = (numpy.exp(-0.5 * ((x[:, None] - cent._data[:, i][None, :]) / abs(width._data[:, i][None, :])) ** 2) / (fact * abs(width._data[:, i][None, :])))
+        spec = numpy.dot(A, rss._data[:, i])
         out[:, i] = spec
-        if i == 1000:
-            plt.plot(spec, "-r")
-            plt.plot(img._data[:, i], "ok")
-            plt.show()
-
+        if i in plot_columns:
+            axs[plot_columns.index(i)].step(x, img._data[:, i], color="k", lw=1, where="mid")
+            axs[plot_columns.index(i)].step(x, spec, color="r", lw=1, where="mid")
+
+    out_path = os.path.dirname(in_image)
+    out_name = os.path.basename(in_image).split(".fits")[0]
+    out_residual = os.path.join(out_path, f"{out_name}_residual.fits")
+    out_2dimage = os.path.join(out_path, f"{out_name}_2dimage.fits")
+    out_ratio = os.path.join(out_path, f"{out_name}_ratio.fits")
+    save_fig(fig, product_path=out_2dimage, to_display=display_plots, figure_path="qa", label="2D_extracted_model")
+
+    log.info(f"writing residual to {out_residual}")
     hdu = pyfits.PrimaryHDU(img._data - out)
-    hdu.writeto("res.fits", overwrite=True)
-    hdu = pyfits.PrimaryHDU(out)
-    hdu.writeto("fit.fits", overwrite=True)
+    hdu.writeto(out_residual, overwrite=True)
 
-    hdu = pyfits.PrimaryHDU((img._data - out) / img._data)
-    hdu.writeto("res_rel.fits", overwrite=True)
+    log.info(f"writing ratio to {out_ratio}")
+    hdu = pyfits.PrimaryHDU(out / img._data)
+    hdu.writeto(out_ratio, overwrite=True)
+
+    log.info(f"writing 2D model {out_2dimage}")
+    hdu = pyfits.PrimaryHDU(out)
+    hdu.writeto(out_2dimage, overwrite=True)
 
 
 # TODO: for arcs take short exposures for bright lines & long exposures for faint lines