Merge pull request #1861 from pypeit/tilts_fwhm

Allow `fwhm_fromlines` for tilts calibration

doc/releases/1.16.1dev.rst

@@ -42,6 +42,7 @@ Functionality/Performance Improvements and Additions
 - Added a functionality that allows, when multiple frames are combined, to scale each
   frame to have the same mean value before combining. To use this
   functionality, the new parameter ``scale_mean`` should be set to ``True``.
+- Added the possibility to use the parameter ``fwhm_fromlines`` also for the tilts calibration.
 
 Instrument-specific Updates
 ---------------------------

pypeit/calibrations.py

@@ -1176,11 +1176,14 @@ def get_tilts(self):
         _spat_flexure = self.mstilt.spat_flexure \
             if self.par['tiltframe']['process']['spat_flexure_correct'] else None
 
+        # get measured fwhm from wv_calib
+        measured_fwhms = [wvfit.fwhm for wvfit in self.wv_calib.wv_fits]
+
         # Build
         buildwaveTilts = wavetilts.BuildWaveTilts(
             self.mstilt, self.slits, self.spectrograph, self.par['tilts'],
             self.par['wavelengths'], det=self.det, qa_path=self.qa_path,
-            spat_flexure=_spat_flexure)
+            spat_flexure=_spat_flexure, measured_fwhms=measured_fwhms)
 
         # TODO still need to deal with syntax for LRIS ghosts. Maybe we don't need it
         self.wavetilts = buildwaveTilts.run(doqa=self.write_qa, show=self.show)

pypeit/par/pypeitpar.py

@@ -2978,15 +2978,16 @@ def __init__(self, reference=None, method=None, echelle=None, ech_nspec_coeff=No
         defaults['fwhm'] = 4.
         dtypes['fwhm'] = [int, float]
         descr['fwhm'] = 'Spectral sampling of the arc lines. This is the FWHM of an arcline in ' \
-                        'binned pixels of the input arc image'
+                        'binned pixels of the input arc image. Note that this is used also in the ' \
+                        'wave tilts calibration.'
 
         defaults['fwhm_fromlines'] = True
         dtypes['fwhm_fromlines'] = bool
         descr['fwhm_fromlines'] = 'Estimate spectral resolution in each slit using the arc lines. '\
-                                  'If True, the estimated FWHM will override ``fwhm`` only in '\
+                                  'If True, the estimated FWHM will override ``fwhm`` in '\
                                   'the determination of the wavelength solution (including the ' \
                                   'calculation of the threshold for the solution RMS, see ' \
-                                  '``rms_thresh_frac_fwhm``), but not for the wave tilts calibration.' \
+                                  '``rms_thresh_frac_fwhm``), and ALSO for the wave tilts calibration.' \
 
         defaults['fwhm_spat_order'] = 0
         dtypes['fwhm_spat_order'] = int

pypeit/scripts/chk_tilts.py

@@ -32,14 +32,18 @@ def get_parser(cls, width=None):
 
     @staticmethod
     def main(args):
+        from pathlib import Path
         from pypeit import wavetilts
 
         chk_version = not args.try_old
 
+        # tilts file path
+        file = Path(args.file).absolute()
+
         # Load
-        tilts = wavetilts.WaveTilts.from_file(args.file, chk_version=chk_version)
+        tilts = wavetilts.WaveTilts.from_file(file, chk_version=chk_version)
         tilts.show(in_ginga=np.logical_not(args.mpl), show_traces=args.show_traces,
-                   chk_version=chk_version)
+                   calib_dir=file.parent, chk_version=chk_version)
 
 
 

pypeit/wavetilts.py

@@ -25,6 +25,7 @@
 from pypeit.display import display
 from pypeit.core import arc
 from pypeit.core import tracewave
+from pypeit.core.wavecal import autoid
 from pypeit.images import buildimage
 
 
@@ -174,7 +175,7 @@ def spatid_to_zero(self, spat_id):
         return np.where(mtch)[0][0]
 
     def show(self, waveimg=None, wcs_match=True, in_ginga=True, show_traces=False,
-             chk_version=True):
+             calib_dir=None, chk_version=True):
         """
         Show in ginga or mpl Tiltimg with the tilts traced and fitted overlaid
 
@@ -188,22 +189,32 @@ def show(self, waveimg=None, wcs_match=True, in_ginga=True, show_traces=False,
                 If True, show the image in ginga. Otherwise, use matplotlib.
             show_traces (bool, optional):
                 If True, show the traces of the tilts on the image.
+            calib_dir (`Path`_):
+                Path to the calibration directory.  If None, the path is taken from the
+                WaveTilts object.
             chk_version (:obj:`bool`, optional):
                 When reading in existing files written by PypeIt, perform strict
                 version checking to ensure a valid file.  If False, the code
                 will try to keep going, but this may lead to faults and quiet
                 failures.  User beware!
         """
         # get tilt_img_dict
-        cal_file = Path(self.calib_dir).absolute() / self.tiltimg_filename
+        _calib_dir = self.calib_dir
+        if calib_dir is not None and calib_dir.exists():
+            _calib_dir = calib_dir
+            msgs.info(f'Searching for other calibration files in {str(_calib_dir)}')
+        else:
+            msgs.info(f'Searching for other calibration files in the default directory {str(_calib_dir)}')
+
+        cal_file = Path(_calib_dir).absolute() / self.tiltimg_filename
         if cal_file.exists():
             tilt_img_dict = buildimage.TiltImage.from_file(cal_file, chk_version=chk_version)
         else:
             msgs.error(f'Tilt image {str(cal_file)} NOT FOUND.')
 
         # get slits
         slitmask = None
-        cal_file = Path(self.calib_dir).absolute() / self.slits_filename
+        cal_file = Path(_calib_dir).absolute() / self.slits_filename
         if cal_file.exists():
             slits = slittrace.SlitTraceSet.from_file(cal_file, chk_version=chk_version)
             _slitmask = slits.slit_img(initial=True, flexure=self.spat_flexure)
@@ -220,7 +231,7 @@ def show(self, waveimg=None, wcs_match=True, in_ginga=True, show_traces=False,
         # get waveimg
         same_size = (slits.nspec, slits.nspat) == tilt_img_dict.image.shape
         if waveimg is None and slits is not None and same_size and in_ginga:
-            wv_calib_name = wavecalib.WaveCalib.construct_file_name(self.calib_key, calib_dir=self.calib_dir)
+            wv_calib_name = wavecalib.WaveCalib.construct_file_name(self.calib_key, calib_dir=_calib_dir)
             if Path(wv_calib_name).absolute().exists():
                 wv_calib = wavecalib.WaveCalib.from_file(wv_calib_name, chk_version=chk_version)
                 tilts = self.fit2tiltimg(slitmask, flexure=self.spat_flexure)
@@ -268,24 +279,32 @@ class BuildWaveTilts:
         slits (:class:`~pypeit.slittrace.SlitTraceSet`):
             Slit edges
         spectrograph (:class:`~pypeit.spectrographs.spectrograph.Spectrograph`):
-            Spectrograph object
+            The `Spectrograph` instance that sets the instrument used.  Used to set
+            :attr:`spectrograph`.
         par (:class:`~pypeit.par.pypeitpar.WaveTiltsPar` or None):
-            The parameters used to fuss with the tilts
+            The parameters used for the tilt calibration.
+            Uses ``['calibrations']['tilts']``.
         wavepar (:class:`~pypeit.par.pypeitpar.WavelengthSolutionPar` or None):
-            The parameters used for the wavelength solution
-        det (int): Detector index
+            The parameters used for the wavelength solution.
+            Uses ``['calibrations']['wavelengths']``.
+        det (int):
+            Detector index
         qa_path (:obj:`str`, optional):
             Directory for QA output.
         spat_flexure (float, optional):
             If input, the slitmask and slit edges are shifted prior
             to tilt analysis.
+        measured_fwhms (`numpy.ndarray`_, optional):
+            FWHM of the arc lines measured during wavecalib.  If provided, this
+            will be used for arc/sky line detection.
+
 
 
     Attributes:
         spectrograph (:class:`~pypeit.spectrographs.spectrograph.Spectrograph`):
-            ??
+            The `Spectrograph` instance that sets the instrument used.
         steps (list):
-            ??
+            List of the processing steps performed
         mask (`numpy.ndarray`_):
             boolean array; True = Ignore this slit
         all_trcdict (list):
@@ -304,7 +323,7 @@ class BuildWaveTilts:
 
     # TODO This needs to be modified to take an inmask
     def __init__(self, mstilt, slits, spectrograph, par, wavepar, det=1, qa_path=None,
-                 spat_flexure=None):
+                 spat_flexure=None, measured_fwhms=None):
 
         # TODO: Perform type checking
         self.spectrograph = spectrograph
@@ -316,6 +335,7 @@ def __init__(self, mstilt, slits, spectrograph, par, wavepar, det=1, qa_path=Non
         self.det = det
         self.qa_path = qa_path
         self.spat_flexure = spat_flexure
+        self.measured_fwhms = measured_fwhms if measured_fwhms is not None else np.array([None] * slits.nslits)
 
         # Get the non-linear count level
         if self.mstilt.is_mosaic:
@@ -385,23 +405,26 @@ def extract_arcs(self):
 
         return arccen, arccen_bpm
 
-    def find_lines(self, arcspec, slit_cen, slit_idx, bpm=None, debug=False):
+    def find_lines(self, arcspec, slit_cen, slit_idx, fwhm, bpm=None, debug=False):
         """
         Find the lines for tracing
 
         Wrapper to tracewave.tilts_find_lines()
 
         Args:
-            arcspec ():
-                ??
-            slit_cen ():
-                ??
+            arcspec (`numpy.ndarray`_):
+                1D spectrum to be searched for significant detections.
+            slit_cen (`numpy.ndarray`_):
+                Trace down the center of the slit. Must match the shape of arcspec.
             slit_idx (int):
-                Slit index, zero-based
+                Slit index, zero-based.
+            fwhm (float):
+                FWHM of the arc lines.
             bpm (`numpy.ndarray`_, optional):
-                ??
+                Bad-pixel mask for input spectrum. If None, all pixels considered good.
+                If passed, must match the shape of arcspec.
             debug (bool, optional):
-                ??
+                Show a QA plot for the line detection.
 
         Returns:
             tuple:  2 objectcs
@@ -423,7 +446,7 @@ def find_lines(self, arcspec, slit_cen, slit_idx, bpm=None, debug=False):
                                              sig_neigh=self.par['sig_neigh'],
                                              nfwhm_neigh=self.par['nfwhm_neigh'],
                                              only_these_lines=only_these_lines,
-                                             fwhm=self.wavepar['fwhm'],
+                                             fwhm=fwhm,
                                              nonlinear_counts=self.nonlinear_counts,
                                              bpm=bpm, debug_peaks=False, debug_lines=debug)
 
@@ -478,7 +501,7 @@ def fit_tilts(self, trc_tilt_dict, thismask, slit_cen, spat_order, spec_order, s
         self.steps.append(inspect.stack()[0][3])
         return self.all_fit_dict[slit_idx]['coeff2']
 
-    def trace_tilts(self, arcimg, lines_spec, lines_spat, thismask, slit_cen,
+    def trace_tilts(self, arcimg, lines_spec, lines_spat, thismask, slit_cen, fwhm,
                     debug_pca=False, show_tracefits=False):
         """
         Trace the tilts
@@ -500,14 +523,20 @@ def trace_tilts(self, arcimg, lines_spec, lines_spat, thismask, slit_cen,
                is (nspec, nspat) with dtype=bool.
             slit_cen (:obj:`int`):
                 Integer index indicating the slit in question.
+            fwhm (:obj:`float`):
+                FWHM of the arc lines.
+            debug_pca (:obj:`bool`, optional):
+                Show the PCA modeling QA plots.
+            show_tracefits (:obj:`bool`, optional):
+                Show the trace fits.
 
         Returns:
             dict: Dictionary containing information on the traced tilts required
             to fit the filts.
 
         """
         trace_dict = tracewave.trace_tilts(arcimg, lines_spec, lines_spat, thismask, slit_cen,
-                                           inmask=self.gpm, fwhm=self.wavepar['fwhm'],
+                                           inmask=self.gpm, fwhm=fwhm,
                                            spat_order=self.par['spat_order'],
                                            maxdev_tracefit=self.par['maxdev_tracefit'],
                                            sigrej_trace=self.par['sigrej_trace'],
@@ -561,11 +590,13 @@ def model_arc_continuum(self, debug=False):
         for i in range(nslits):
             if self.tilt_bpm[i]:
                 continue
+            # get FWHM for this slit
+            fwhm = autoid.set_fwhm(self.wavepar, measured_fwhm=self.measured_fwhms[i])
             # TODO: What to do with the following iter_continuum parameters?:
             #       sigthresh, sigrej, niter_cont, cont_samp, cont_frac_fwhm
             arc_continuum[:,i], arc_fitmask[:,i] \
                     = arc.iter_continuum(self.arccen[:,i], gpm=np.invert(self.arccen_bpm[:,i]),
-                                         fwhm=self.wavepar['fwhm'])
+                                         fwhm=fwhm)
             # TODO: Original version.  Please leave it for now.
 #            arc_fitmask[:,i], coeff \
 #                    = utils.robust_polyfit_djs(spec, self.arccen[:,i], self.par['cont_order'],
@@ -677,6 +708,7 @@ def run(self, doqa=True, debug=False, show=False):
         # Subtract arc continuum
         _mstilt = self.mstilt.image.copy()
         if self.par['rm_continuum']:
+            msgs.info('Subtracting the continuum')
             continuum = self.model_arc_continuum(debug=debug)
             _mstilt -= continuum
             if debug:
@@ -713,11 +745,13 @@ def run(self, doqa=True, debug=False, show=False):
                 self.slits.mask[slit_idx] = self.slits.bitmask.turn_on(self.slits.mask[slit_idx], 'BADTILTCALIB')
                 continue
             msgs.info(f'Computing tilts for slit/order {self.slits.slitord_id[slit_idx]} ({slit_idx+1}/{self.slits.nslits})')
+            # Get the arc FWHM for this slit
+            fwhm = autoid.set_fwhm(self.wavepar, measured_fwhm=self.measured_fwhms[slit_idx], verbose=True)
             # Identify lines for tracing tilts
             msgs.info('Finding lines for tilt analysis')
             self.lines_spec, self.lines_spat \
                     = self.find_lines(self.arccen[:,slit_idx], self.slitcen[:,slit_idx],
-                                      slit_idx,
+                                      slit_idx, fwhm,
                                       bpm=self.arccen_bpm[:,slit_idx], debug=debug)
 
             if self.lines_spec is None:
@@ -733,7 +767,7 @@ def run(self, doqa=True, debug=False, show=False):
             # each line.
             msgs.info('Trace the tilts')
             self.trace_dict = self.trace_tilts(_mstilt, self.lines_spec, self.lines_spat,
-                                               thismask, self.slitcen[:, slit_idx])
+                                               thismask, self.slitcen[:, slit_idx], fwhm)
             # IF there are < 2 usable arc lines for tilt tracing, PCA fit does not work and the reduction crushes
             # TODO investigate why some slits have <2 usable arc lines
             if np.sum(self.trace_dict['use_tilt']) < 2: