fixing stretching factors in wavelength/LSF fitting & improving cross…

… matching

python/lvmdrp/core/fiberrows.py

@@ -870,6 +870,8 @@ def measureArcLines(
         bg = numpy.ones((self._fibers, nlines), dtype=numpy.float32) * numpy.nan
         masked = numpy.zeros((self._fibers, nlines), dtype="bool")
 
+        stretch_min, stretch_max, stretch_steps = 0.998, 1.002, 40
+
         spec = self.getSpec(ref_fiber)
         flux[ref_fiber], cent_wave[ref_fiber], fwhm[ref_fiber], bg[ref_fiber] = spec.fitSepGauss(ref_cent, aperture, fwhm_guess, bg_guess, flux_range, cent_range, fwhm_range, bg_range, axs=axs[ref_fiber][1])
         masked[ref_fiber] = numpy.isnan(flux[ref_fiber])|numpy.isnan(cent_wave[ref_fiber])|numpy.isnan(fwhm[ref_fiber])
@@ -899,15 +901,19 @@ def measureArcLines(
             cc, bhat, mhat = _cross_match_float(
                 ref_spec=last_spec._data,
                 obs_spec=spec._data,
-                stretch_factors=numpy.linspace(0.99,1.01,20),
-                shift_range=[-5, 5],
+                stretch_factors=numpy.linspace(stretch_min, stretch_max, stretch_steps),
+                shift_range=[-10, 10],
                 normalize_spectra=False,
             )
+            if mhat == stretch_min or mhat == stretch_max:
+                log.warning(f"boundary of stretch factors: {mhat = } ({stretch_min, stretch_max = })")
             cent_guess = mhat * last_cent + bhat
             flux[i], cent_wave[i], fwhm[i], bg[i] = spec.fitSepGauss(cent_guess, aperture, fwhm_guess, bg_guess, flux_range, cent_range, fwhm_range, bg_range, axs=axs_fiber)
             masked[i] = numpy.isnan(flux[i])|numpy.isnan(cent_wave[i])|numpy.isnan(fwhm[i])
             if masked[i].any():
                 log.warning(f"some lines were not fitted properly in fiber {i}: ")
+                log.warning(f"  guess = {numpy.round(cent_guess, 3)} ({mhat = :.5f}, {bhat = :.5f})")
+                log.warning(f"   mask = {masked[i]}")
                 log.warning(f"   flux = {numpy.round(flux[i],3)}")
                 log.warning(f"   cent = {numpy.round(cent_wave[i],3)}")
                 log.warning(f"   fwhm = {numpy.round(fwhm[i],3)}")
@@ -939,13 +945,23 @@ def measureArcLines(
             cc, bhat, mhat = _cross_match_float(
                 ref_spec=last_spec._data,
                 obs_spec=spec._data,
-                stretch_factors=numpy.linspace(0.99,1.01,20),
-                shift_range=[-5, 5],
-                normalize_spectra=False
+                stretch_factors=numpy.linspace(stretch_min, stretch_max, stretch_steps),
+                shift_range=[-10, 10],
+                normalize_spectra=False,
             )
+            if mhat == stretch_min or mhat == stretch_max:
+                log.warning(f"boundary of stretch factors: {mhat = } ({stretch_min, stretch_max = })")
             cent_guess = mhat * last_cent + bhat
             flux[i], cent_wave[i], fwhm[i], bg[i] = spec.fitSepGauss(cent_guess, aperture, fwhm_guess, bg_guess, flux_range, cent_range, fwhm_range, bg_range, axs=axs_fiber)
             masked[i] = numpy.isnan(flux[i])|numpy.isnan(cent_wave[i])|numpy.isnan(fwhm[i])
+            if masked[i].any():
+                log.warning(f"some lines were not fitted properly in fiber {i}: ")
+                log.warning(f"  guess = {numpy.round(cent_guess, 3)} ({mhat = :.5f}, {bhat = :.5f})")
+                log.warning(f"   mask = {masked[i]}")
+                log.warning(f"   flux = {numpy.round(flux[i],3)}")
+                log.warning(f"   cent = {numpy.round(cent_wave[i],3)}")
+                log.warning(f"   fwhm = {numpy.round(fwhm[i],3)}")
+                log.warning(f"   bg   = {numpy.round(bg[i],3)}")
 
             last_spec = copy(spec)
             last_cent = copy(cent_wave[i])

python/lvmdrp/core/spectrum1d.py

@@ -290,6 +290,9 @@ def _cross_match_float(
     best_shift = 0
     best_stretch_factor = 1
 
+    # define pixels array centered around middle part of the spectrum
+    # pixels = numpy.arange(ref_spec.size) - ref_spec.size // 2
+
     # normalize the peaks to roughly magnitude 1, so that individual very bright
     # fibers do not dominate the signal
     if normalize_spectra:
@@ -303,6 +306,7 @@ def _cross_match_float(
     for factor in stretch_factors:
         # Stretch the first signal
         stretched_signal1 = zoom(ref_spec_, factor, mode="constant", prefilter=True)
+        # stretched_signal1 = numpy.interp(pixels * factor, pixels, ref_spec_)
 
         # Make the lengths equal
         len_diff = len(obs_spec_) - len(stretched_signal1)
@@ -3496,7 +3500,7 @@ def obtainGaussFluxPeaks(self, pos, sigma, replace_error=1e10, plot=False):
         if bad_pix is not None and bn.nansum(bad_pix) > 0:
             error[bad_pix] = replace_error
 
-        # pyfits.writeto('B.fits', B.toarray(), overwrite=True)    
+        # pyfits.writeto('B.fits', B.toarray(), overwrite=True)
         # if plot:
         #     plt.plot(self._data, "ok")
         #     plt.plot(numpy.dot(A * self._error[:, None], out[0]), "-r")

python/lvmdrp/functions/rssMethod.py

@@ -364,14 +364,16 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf
 
         # fix cc_max_shift
         # cross-match spectrum and pixwav map
+        stretch_min, stretch_max, stretch_steps = 0.95, 1.05, 10000
         cc, bhat, mhat = _cross_match_float(
             ref_spec=pix_spec,
             obs_spec=arc._data[ref_fiber],
-            stretch_factors=numpy.linspace(0.8,1.2,10000),
+            stretch_factors=numpy.linspace(stretch_min, stretch_max, stretch_steps),
             shift_range=[-cc_max_shift, cc_max_shift],
             normalize_spectra=False,
         )
-
+        if mhat == stretch_min or mhat == stretch_max:
+            log.warning(f"boundary of stretch factors: {mhat = } ({stretch_min, stretch_max = })")
         log.info(f"max CC = {cc:.2f} for strech = {mhat:.8f} and shift = {bhat:.8f}")
     else:
         mhat, bhat = 1.0, 0.0