From 804a451c5c8bc1e9b91b9f7c96977dcb51bdfd51 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 13 Sep 2023 12:44:34 -1000 Subject: [PATCH 01/12] refactored in wavecalib.py and autoid.py --- pypeit/core/wavecal/autoid.py | 86 +++++++++++++++++++++++------------ pypeit/wavecalib.py | 21 +++++---- 2 files changed, 69 insertions(+), 38 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index d4ba5fd91f..311b61b770 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -519,7 +519,7 @@ def reidentify(spec, spec_arxiv_in, wave_soln_arxiv_in, line_list, # Search for lines no matter what to continuum subtract the input arc tcent, ecent, cut_tcent, icut, spec_cont_sub = wvutils.arc_lines_from_spec( - spec, sigdetect=sigdetect,nonlinear_counts=nonlinear_counts, fwhm = fwhm, debug = debug_peaks) + spec, sigdetect=sigdetect, nonlinear_counts=nonlinear_counts, fwhm=fwhm, debug=debug_peaks) # If the detections were not passed in measure them if detections is None: detections = tcent[icut] @@ -531,7 +531,7 @@ def reidentify(spec, spec_arxiv_in, wave_soln_arxiv_in, line_list, det_arxiv1 = {} for iarxiv in range(narxiv): tcent_arxiv, ecent_arxiv, cut_tcent_arxiv, icut_arxiv, spec_cont_sub_now = wvutils.arc_lines_from_spec( - spec_arxiv[:,iarxiv], sigdetect=sigdetect,nonlinear_counts=nonlinear_counts, fwhm = fwhm, debug = debug_peaks) + spec_arxiv[:,iarxiv], sigdetect=sigdetect, nonlinear_counts=nonlinear_counts, fwhm=fwhm, debug=debug_peaks) spec_arxiv_cont_sub[:,iarxiv] = spec_cont_sub_now det_arxiv1[str(iarxiv)] = tcent_arxiv[icut_arxiv] @@ -764,7 +764,7 @@ def match_to_arxiv(lamps:list, spec:np.ndarray, wv_guess:np.ndarray, tcent, ecent, cut_tcent, icut, spec_cont_sub = wvutils.arc_lines_from_spec( spec, sigdetect=sigdetect, nonlinear_counts=nonlinear_counts, - fwhm=fwhm, debug = debug_peaks) + fwhm=fwhm, debug=debug_peaks) # If there are no lines in the input arc, return if tcent.size == 0: return None @@ -772,7 +772,7 @@ def match_to_arxiv(lamps:list, spec:np.ndarray, wv_guess:np.ndarray, # Search for lines in the arxiv arc tcent_arxiv, ecent_arxiv, cut_tcent_arxiv, icut_arxiv, spec_cont_sub_now = wvutils.arc_lines_from_spec( spec_arxiv, sigdetect=sigdetect, - nonlinear_counts=nonlinear_counts, fwhm = fwhm, debug = debug_peaks) + nonlinear_counts=nonlinear_counts, fwhm=fwhm, debug=debug_peaks) # If there are no lines in the arxiv arc, return if tcent_arxiv.size == 0: return None @@ -1114,7 +1114,7 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, ax = plt.gca() # ax.plot(xvals, tspec, label='template') # Template - ax.plot(xvals, np.roll(pspec, int(shift_cc)), 'k', label='input') # Input + ax.plot(xvals, np.roll(pad_spec, int(shift_cc)), 'k', label='input') # Input ax.legend() plt.show() #embed(header='909 autoid') @@ -1195,8 +1195,8 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, # Finish return wvcalib -def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, - ok_mask=None, use_unknowns=True, debug_all=False, +def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, rms_thresh, + ok_mask=None, measured_fwhms=None, use_unknowns=True, debug_all=False, debug_peaks=False, debug_xcorr=False, debug_reid=False, debug_fits=False, nonlinear_counts=1e10, redo_slits:list=None): @@ -1222,9 +1222,13 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, par : :class:`~pypeit.par.pypeitpar.WavelengthSolutionPar` Key parameters that drive the behavior of the wavelength-solution algorithms. + rms_thresh: float + RMS threshold for the wavelength solution fit ok_mask : `numpy.ndarray`, optional Integer array with the list of valid spectra ``spec`` to use. If None, all spectra are used. + measured_fwhms: ndarray, optional + Array of FWHM (in binned pixels) measured from the arc lines. Shape :math:`(N_{\rm orders},)` use_unknowns : bool, default = True, optional If True, arc lines that are known to be present in the spectra, but have not been attributed to an element+ion, will @@ -1305,14 +1309,16 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, msgs.info('Reidentifying and fitting Order = {0:d}, which is {1:d}/{2:d}'.format(orders[iord], iord+1, norders)) sigdetect = wvutils.parse_param(par, 'sigdetect', iord) cc_thresh = wvutils.parse_param(par, 'cc_thresh', iord) - rms_threshold = wvutils.parse_param(par, 'rms_threshold', iord) + rms_thresh_frac_fwhm = par['rms_thresh_frac_fwhm'] msgs.info("Using sigdetect = {}".format(sigdetect)) - msgs.info("Using rms_threshold = {}".format(rms_threshold)) + msgs.info(f"Using rms_threshold = {rms_thresh} ({rms_thresh_frac_fwhm} * median FWHM)") + # Set FWHM for this order + fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[iord]) detections[str(iord)], spec_cont_sub[:, iord], all_patt_dict[str(iord)] = reidentify( spec[:, iord], spec_arxiv[:, iord], wave_arxiv[:, iord], tot_line_list, par['nreid_min'], cc_thresh=cc_thresh, match_toler=par['match_toler'], cc_local_thresh=par['cc_local_thresh'], nlocal_cc=par['nlocal_cc'], - nonlinear_counts=nonlinear_counts, sigdetect=sigdetect, fwhm=par['fwhm'], + nonlinear_counts=nonlinear_counts, sigdetect=sigdetect, fwhm=fwhm, debug_peaks=(debug_peaks or debug_all), debug_xcorr=(debug_xcorr or debug_all), debug_reid=(debug_reid or debug_all)) @@ -1342,7 +1348,7 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, bad_orders = np.append(bad_orders, iord) continue # Is the RMS below the threshold? - if final_fit['rms'] > rms_threshold: + if final_fit['rms'] > rms_thresh: msgs.warn('---------------------------------------------------' + msgs.newline() + 'Reidentify report for slit {0:d}/{1:d}:'.format(iord+1, norders) + msgs.newline() + ' Poor RMS ({0:.3f})! Need to add additional spectra to arxiv to improve fits'.format( @@ -1444,6 +1450,8 @@ class ArchiveReid: par : :class:`~pypeit.par.pypeitpar.WavelengthSolutionPar` Key parameters that drive the behavior of the wavelength-solution algorithms. + rms_thresh: float + RMS threshold for the wavelength solution fit ech_fixed_format: bool Set to True if this is a fixed format echelle spectrograp. The code will then align the archive_arc and the extracted arc for each order for the reidentification. @@ -1451,7 +1459,7 @@ class ArchiveReid: Integer array with the list of valid spectra ``spec`` to use. If None, all spectra are used. measured_fwhms: ndarray, optional - Array of FWHM (in pixels) measured from the arc lines. Shape (nslit,) + Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,) use_unknowns : bool, default = True, optional If True, arc lines that are known to be present in the spectra, but have not been attributed to an element+ion, will @@ -1492,7 +1500,8 @@ class ArchiveReid: """ # TODO: Because we're passing orders directly, we no longer need spectrograph... - def __init__(self, spec, lamps, par, ech_fixed_format=False, ok_mask=None, measured_fwhms=None, use_unknowns=True, debug_all=False, + def __init__(self, spec, lamps, par, rms_thresh, ech_fixed_format=False, ok_mask=None, + measured_fwhms=None, use_unknowns=True, debug_all=False, debug_peaks=False, debug_xcorr=False, debug_reid=False, debug_fits=False, orders=None, nonlinear_counts=1e10): @@ -1616,9 +1625,9 @@ def __init__(self, spec, lamps, par, ech_fixed_format=False, ok_mask=None, meas msgs.info(f'Order: {orders[slit]}') sigdetect = wvutils.parse_param(self.par, 'sigdetect', slit) cc_thresh = wvutils.parse_param(self.par, 'cc_thresh', slit) - rms_threshold = wvutils.parse_param(self.par, 'rms_threshold', slit) + rms_thresh_frac_fwhm = self.par['rms_thresh_frac_fwhm'] msgs.info("Using sigdetect = {}".format(sigdetect)) - msgs.info("Using rms_threshold = {}".format(rms_threshold)) + msgs.info(f"Using rms_threshold = {rms_thresh} ({rms_thresh_frac_fwhm} * median FWHM)") # get FWHM for this slit fwhm = set_fwhm(self.par, measured_fwhm=measured_fwhms[slit]) self.detections[str(slit)], self.spec_cont_sub[:,slit], self.all_patt_dict[str(slit)] = \ @@ -1648,7 +1657,7 @@ def __init__(self, spec, lamps, par, ech_fixed_format=False, ok_mask=None, meas self.bad_slits = np.append(self.bad_slits, slit) continue # Is the RMS below the threshold? - if final_fit['rms'] > rms_threshold: + if final_fit['rms'] > rms_thresh: order_str = '' if orders is None else ', order={}'.format(orders[slit]) msgs.warn('---------------------------------------------------' + msgs.newline() + 'Reidentify report for slit {0:d}/{1:d}'.format(slit, self.nslits-1) + order_str + msgs.newline() + @@ -1711,6 +1720,10 @@ class HolyGrail: islinelist : bool, optional Is lamps a linelist (True), or a list of ions (False) The former is not recommended except by expert users/developers + rms_thresh : float, optional + RMS threshold for the wavelength solution fit + measured_fwhms : ndarray, optional + Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,) outroot : str, optional Name of output file debug : bool, optional @@ -1744,7 +1757,8 @@ class HolyGrail: """ - def __init__(self, spec, lamps, par=None, ok_mask=None, islinelist=False, + def __init__(self, spec, lamps, par=None, ok_mask=None, + islinelist=False, rms_thresh=None, measured_fwhms=None, outroot=None, debug=False, verbose=False, binw=None, bind=None, nstore=1, use_unknowns=True, nonlinear_counts=None, spectrograph=None): @@ -1757,6 +1771,10 @@ def __init__(self, spec, lamps, par=None, ok_mask=None, islinelist=False, self._nstore = nstore self._binw = binw self._bind = bind + self._measured_fwhms = measured_fwhms + self.rms_thresh = \ + round(self._par['rms_thresh_frac_fwhm'] * set_fwhm(self._par, measured_fwhm=np.median(measured_fwhms))) \ + if rms_thresh is None else rms_thresh # Mask info if ok_mask is None: @@ -1851,8 +1869,8 @@ def run_brute_loop(self, slit, tcent_ecent, wavedata=None): idthresh = 0.5 # Criteria for early return (at least this fraction of lines must have # an ID on either side of the spectrum) - rms_threshold = wvutils.parse_param(self._par, 'rms_threshold', slit) - msgs.info("Using rms_threshold = {}".format(rms_threshold)) + rms_thresh_frac_fwhm = self._par['rms_thresh_frac_fwhm'] + msgs.info(f"Using rms_threshold = {self.rms_thresh} ({rms_thresh_frac_fwhm} * median FWHM)") best_patt_dict, best_final_fit = None, None # Loop through parameter space for poly in rng_poly: @@ -1872,7 +1890,7 @@ def run_brute_loop(self, slit, tcent_ecent, wavedata=None): # First time a fit is found best_patt_dict, best_final_fit = copy.deepcopy(patt_dict), copy.deepcopy(final_fit) continue - elif final_fit['rms'] < rms_threshold: + elif final_fit['rms'] < self.rms_thresh: # Has a better fit been identified (i.e. more lines identified)? if len(final_fit['pixel_fit']) > len(best_final_fit['pixel_fit']): best_patt_dict, best_final_fit = copy.deepcopy(patt_dict), copy.deepcopy(final_fit) @@ -1907,10 +1925,14 @@ def run_brute(self, min_nlines=10): # Detect lines, and decide which tcent to use sigdetect = wvutils.parse_param(self._par, 'sigdetect', slit) msgs.info("Using sigdetect = {}".format(sigdetect)) - self._all_tcent, self._all_ecent, self._cut_tcent, self._icut, _ =\ - wvutils.arc_lines_from_spec(self._spec[:, slit].copy(), sigdetect=sigdetect, nonlinear_counts = self._nonlinear_counts) - self._all_tcent_weak, self._all_ecent_weak, self._cut_tcent_weak, self._icut_weak, _ =\ - wvutils.arc_lines_from_spec(self._spec[:, slit].copy(), sigdetect=sigdetect, nonlinear_counts = self._nonlinear_counts) + # get FWHM for this slit + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit]) + self._all_tcent, self._all_ecent, self._cut_tcent, self._icut, _ =\ + wvutils.arc_lines_from_spec(self._spec[:, slit].copy(), sigdetect=sigdetect, fwhm=fwhm, + nonlinear_counts=self._nonlinear_counts) + self._all_tcent_weak, self._all_ecent_weak, self._cut_tcent_weak, self._icut_weak, _ =\ + wvutils.arc_lines_from_spec(self._spec[:, slit].copy(), sigdetect=sigdetect, fwhm=fwhm, + nonlinear_counts =self._nonlinear_counts) # Were there enough lines? This mainly deals with junk slits if self._all_tcent.size < min_nlines: @@ -2013,12 +2035,16 @@ def run_kdtree(self, polygon=4, detsrch=7, lstsrch=10, pixtol=5): if slit not in self._ok_mask: self._all_final_fit[str(slit)] = {} continue + # get FWHM for this slit + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit]) # Detect lines, and decide which tcent to use sigdetect = wvutils.parse_param(self._par, 'sigdetect', slit) self._all_tcent, self._all_ecent, self._cut_tcent, self._icut, _ =\ - wvutils.arc_lines_from_spec(self._spec[:, slit], sigdetect=sigdetect, nonlinear_counts = self._nonlinear_counts) + wvutils.arc_lines_from_spec(self._spec[:, slit], sigdetect=sigdetect, fwhm=fwhm, + nonlinear_counts=self._nonlinear_counts) self._all_tcent_weak, self._all_ecent_weak, self._cut_tcent_weak, self._icut_weak, _ =\ - wvutils.arc_lines_from_spec(self._spec[:, slit], sigdetect=sigdetect, nonlinear_counts = self._nonlinear_counts) + wvutils.arc_lines_from_spec(self._spec[:, slit], sigdetect=sigdetect, fwhm=fwhm, + nonlinear_counts = self._nonlinear_counts) if self._all_tcent.size == 0: msgs.warn("No lines to identify in slit {0:d}!".format(slit+ 1)) continue @@ -2332,7 +2358,7 @@ def cross_match(self, good_fit, detections): new_bad_slits = np.append(new_bad_slits, bs) continue - if final_fit['rms'] > wvutils.parse_param(self._par, 'rms_threshold', bs): + if final_fit['rms'] > self.rms_thresh: msgs.warn('---------------------------------------------------' + msgs.newline() + 'Cross-match report for slit {0:d}/{1:d}:'.format(bs + 1, self._nslit) + msgs.newline() + ' Poor RMS ({0:.3f})! Will try cross matching iteratively'.format(final_fit['rms']) + msgs.newline() + @@ -2509,7 +2535,7 @@ def cross_match_order(self, good_fit): ' Fit was not good enough! Will try cross matching iteratively' + msgs.newline() + '---------------------------------------------------') continue - if final_fit['rms'] > wvutils.parse_param(self._par, 'rms_threshold', slit): + if final_fit['rms'] > self.rms_thresh: msgs.warn('---------------------------------------------------' + msgs.newline() + 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + ' Poor RMS ({0:.3f})! Will try cross matching iteratively'.format(final_fit['rms']) + msgs.newline() + @@ -2884,7 +2910,7 @@ def solve_slit(self, slit, psols, msols, tcent_ecent, nstore=1, nselw=3, nseld=3 # First time a fit is found patt_dict, final_dict = tpatt_dict, tfinal_dict continue - elif tfinal_dict['rms'] < wvutils.parse_param(self._par, 'rms_threshold', slit): + elif tfinal_dict['rms'] < self.rms_thresh: # Has a better fit been identified (i.e. more lines ID)? if len(tfinal_dict['pixel_fit']) > len(final_dict['pixel_fit']): patt_dict, final_dict = copy.deepcopy(tpatt_dict), copy.deepcopy(tfinal_dict) @@ -3014,7 +3040,7 @@ def report_prelim(self, slit, best_patt_dict, best_final_fit): '---------------------------------------------------') self._all_patt_dict[str(slit)] = None self._all_final_fit[str(slit)] = None - elif best_final_fit['rms'] > wvutils.parse_param(self._par, 'rms_threshold', slit): + elif best_final_fit['rms'] > self.rms_thresh: msgs.warn('---------------------------------------------------' + msgs.newline() + 'Preliminary report for slit {0:d}/{1:d}:'.format(slit+1, self._nslit) + msgs.newline() + ' Poor RMS ({0:.3f})! Attempting to cross match.'.format(best_final_fit['rms']) + msgs.newline() + diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 0a96906f34..9735b5111f 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -611,6 +611,9 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Save for redo's self.measured_fwhms = measured_fwhms + # determine rms threshold + fwhm_for_thresh = autoid.set_fwhm(self.par, measured_fwhm=np.median(measured_fwhms)) + self.wave_rms_thresh = round(self.par['rms_thresh_frac_fwhm'] * fwhm_for_thresh,3) # Obtain calibration for all slits if method == 'holy-grail': @@ -636,7 +639,7 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Now preferred # Slit positions arcfitter = autoid.ArchiveReid( - arccen, self.lamps, self.par, + arccen, self.lamps, self.par, self.wave_rms_thresh, ech_fixed_format=self.spectrograph.ech_fixed_format, ok_mask=ok_mask_idx, measured_fwhms=self.measured_fwhms, @@ -680,7 +683,8 @@ def build_wv_calib(self, arccen, method, skip_QA=False, #ok_mask_idx = ok_mask_idx[:-1] patt_dict, final_fit = autoid.echelle_wvcalib( arccen, order_vec, arcspec_arxiv, wave_soln_arxiv, - self.lamps, self.par, ok_mask=ok_mask_idx, + self.lamps, self.par, self.wave_rms_thresh, ok_mask=ok_mask_idx, + measured_fwhms=self.measured_fwhms, nonlinear_counts=self.nonlinear_counts, debug_all=False, redo_slits=np.atleast_1d(self.par['redo_slits']) if self.par['redo_slits'] is not None else None) @@ -702,8 +706,8 @@ def build_wv_calib(self, arccen, method, skip_QA=False, self.wv_calib.arc_spectra = arccen # Save the new fits (if they meet tolerance) for key in final_fit.keys(): - if final_fit[key]['rms'] < self.par['rms_threshold']: - idx = int(key) + idx = int(key) + if final_fit[key]['rms'] < self.wave_rms_thresh: self.wv_calib.wv_fits[idx] = final_fit[key] self.wv_calib.wv_fits[idx].spat_id = self.slits.spat_id[idx] self.wv_calib.wv_fits[idx].fwhm = self.measured_fwhms[idx] @@ -801,7 +805,8 @@ def redo_echelle_orders(self, bad_orders:np.ndarray, dets:np.ndarray, order_dets spec_vec_norm = np.linspace(0,1,nspec) wv_order_mod = self.wv_calib.wv_fit2d[idet].eval(spec_vec_norm, x2=np.ones_like(spec_vec_norm)*order)/order - + # get FWHM for this order + fwhm = autoid.set_fwhm(self.par, measured_fwhm=self.measured_fwhms[iord]) # Link me tcent, spec_cont_sub, patt_dict_slit, tot_llist = autoid.match_to_arxiv( self.lamps, self.arccen[:,iord], wv_order_mod, @@ -810,7 +815,7 @@ def redo_echelle_orders(self, bad_orders:np.ndarray, dets:np.ndarray, order_dets match_toler=self.par['match_toler'], nonlinear_counts=self.nonlinear_counts, sigdetect=wvutils.parse_param(self.par, 'sigdetect', iord), - fwhm=self.par['fwhm']) + fwhm=fwhm) if not patt_dict_slit['acceptable']: msgs.warn(f"Order {order} is still not acceptable after attempt to reidentify.") @@ -833,7 +838,7 @@ def redo_echelle_orders(self, bad_orders:np.ndarray, dets:np.ndarray, order_dets # Keep? # TODO -- Make this a parameter? increase_rms = 1.5 - if final_fit['rms'] < increase_rms*self.par['rms_threshold']* np.median(self.measured_fwhms)/self.par['fwhm']: + if final_fit['rms'] < increase_rms*self.wave_rms_thresh: # TODO -- This is repeated from build_wv_calib() # Would be nice to consolidate # QA @@ -1077,7 +1082,7 @@ def run(self, skip_QA=False, debug=False, # Assess the fits rms = np.array([999. if wvfit.rms is None else wvfit.rms for wvfit in self.wv_calib.wv_fits]) # Test and scale by measured_fwhms - bad_rms = rms > (self.par['rms_threshold'] * np.median(self.measured_fwhms)/self.par['fwhm']) + bad_rms = rms > self.wave_rms_thresh #embed(header='line 975 of wavecalib.py') if np.any(bad_rms): self.wvc_bpm[bad_rms] = True From 7120c9c2a6c4d5205e7a2e36ac7a47516902fa4e Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 13 Sep 2023 13:08:16 -1000 Subject: [PATCH 02/12] edit pypeitpar --- pypeit/par/pypeitpar.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index 6108edcaae..57ac0083c9 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -2507,7 +2507,7 @@ class WavelengthSolutionPar(ParSet): def __init__(self, reference=None, method=None, echelle=None, ech_nspec_coeff=None, ech_norder_coeff=None, ech_sigrej=None, lamps=None, sigdetect=None, fwhm=None, fwhm_fromlines=None, fwhm_spat_order=None, fwhm_spec_order=None, reid_arxiv=None, nreid_min=None, cc_thresh=None, cc_local_thresh=None, nlocal_cc=None, - rms_threshold=None, match_toler=None, func=None, n_first=None, n_final=None, + rms_thresh_frac_fwhm=None, match_toler=None, func=None, n_first=None, n_final=None, sigrej_first=None, sigrej_final=None, numsearch=None, nfitpix=None, refframe=None, nsnippet=None, use_instr_flag=None, wvrng_arxiv=None, @@ -2626,12 +2626,13 @@ def __init__(self, reference=None, method=None, echelle=None, ech_nspec_coeff=No descr['fwhm'] = 'Spectral sampling of the arc lines. This is the FWHM of an arcline in ' \ 'binned pixels of the input arc image' - defaults['fwhm_fromlines'] = False + 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 '\ - 'the determination of the wavelength solution (`i.e.`, not in '\ - 'WaveTilts).' + '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.' \ defaults['fwhm_spat_order'] = 0 dtypes['fwhm_spat_order'] = int @@ -2698,12 +2699,15 @@ def __init__(self, reference=None, method=None, echelle=None, ech_nspec_coeff=No 'be added to it to make it odd.' # These are the parameters used for the iterative fitting of the arc lines - defaults['rms_threshold'] = 0.15 - dtypes['rms_threshold'] = float - descr['rms_threshold'] = 'Maximum RMS (in binned pixels) for keeping a slit/order solution. ' \ - 'Used for echelle spectrographs, the \'reidentify\' method, and when re-analyzing a slit with the redo_slits parameter.' \ - 'In a future PR, we will refactor the code to always scale this threshold off the measured FWHM of the arc lines.' - + defaults['rms_thresh_frac_fwhm'] = 0.15 + dtypes['rms_thresh_frac_fwhm'] = float + descr['rms_thresh_frac_fwhm'] = 'Maximum RMS (expressed as fraction of the FWHM) for keeping ' \ + 'a slit/order solution. If ``fwhm_fromlines`` is True, ' \ + 'FWHM will be computed from the arc lines in each slits ' \ + '(the median value among all the slits is used), otherwise ``fwhm`` ' \ + 'will be used. This parameter is used for the \'holy-grail\', ' \ + '\'reidentify\', and \'echelle\' methods and when re-analyzing ' \ + 'a slit using the ``redo_slits`` parameter. ' defaults['match_toler'] = 2.0 dtypes['match_toler'] = float @@ -2785,7 +2789,7 @@ def from_dict(cls, cfg): 'ech_norder_coeff', 'ech_sigrej', 'ech_separate_2d', 'lamps', 'sigdetect', 'fwhm', 'fwhm_fromlines', 'fwhm_spat_order', 'fwhm_spec_order', 'reid_arxiv', 'nreid_min', 'cc_thresh', 'cc_local_thresh', - 'nlocal_cc', 'rms_threshold', 'match_toler', 'func', 'n_first','n_final', + 'nlocal_cc', 'rms_thresh_frac_fwhm', 'match_toler', 'func', 'n_first','n_final', 'sigrej_first', 'sigrej_final', 'numsearch', 'nfitpix', 'refframe', 'nsnippet', 'use_instr_flag', 'wvrng_arxiv', 'redo_slits', 'qa_log'] From f5c3d29c6f0f8859fd77c60f36b9520a71fe2996 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 13 Sep 2023 13:12:55 -1000 Subject: [PATCH 03/12] minor change --- pypeit/core/wavecal/autoid.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index 9c5c8791ac..b586dffaea 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -1310,9 +1310,8 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, rms_thresh msgs.info('Reidentifying and fitting Order = {0:d}, which is {1:d}/{2:d}'.format(orders[iord], iord+1, norders)) sigdetect = wvutils.parse_param(par, 'sigdetect', iord) cc_thresh = wvutils.parse_param(par, 'cc_thresh', iord) - rms_thresh_frac_fwhm = par['rms_thresh_frac_fwhm'] msgs.info("Using sigdetect = {}".format(sigdetect)) - msgs.info(f"Using rms_threshold = {rms_thresh} ({rms_thresh_frac_fwhm} * median FWHM)") + msgs.info(f"Using rms_threshold = {rms_thresh} ({par['rms_thresh_frac_fwhm']} * median FWHM)") # Set FWHM for this order fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[iord]) detections[str(iord)], spec_cont_sub[:, iord], all_patt_dict[str(iord)] = reidentify( @@ -1626,9 +1625,8 @@ def __init__(self, spec, lamps, par, rms_thresh, ech_fixed_format=False, ok_mask msgs.info(f'Order: {orders[slit]}') sigdetect = wvutils.parse_param(self.par, 'sigdetect', slit) cc_thresh = wvutils.parse_param(self.par, 'cc_thresh', slit) - rms_thresh_frac_fwhm = self.par['rms_thresh_frac_fwhm'] msgs.info("Using sigdetect = {}".format(sigdetect)) - msgs.info(f"Using rms_threshold = {rms_thresh} ({rms_thresh_frac_fwhm} * median FWHM)") + msgs.info(f"Using rms_threshold = {rms_thresh} ({self.par['rms_thresh_frac_fwhm']} * median FWHM)") # get FWHM for this slit fwhm = set_fwhm(self.par, measured_fwhm=measured_fwhms[slit]) self.detections[str(slit)], self.spec_cont_sub[:,slit], self.all_patt_dict[str(slit)] = \ @@ -1870,8 +1868,7 @@ def run_brute_loop(self, slit, tcent_ecent, wavedata=None): idthresh = 0.5 # Criteria for early return (at least this fraction of lines must have # an ID on either side of the spectrum) - rms_thresh_frac_fwhm = self._par['rms_thresh_frac_fwhm'] - msgs.info(f"Using rms_threshold = {self.rms_thresh} ({rms_thresh_frac_fwhm} * median FWHM)") + msgs.info(f"Using rms_threshold = {self.rms_thresh} ({self._par['rms_thresh_frac_fwhm']} * median FWHM)") best_patt_dict, best_final_fit = None, None # Loop through parameter space for poly in rng_poly: From 0f9c8de51fa213c714c89791ab134081baee9b99 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 13 Sep 2023 13:47:28 -1000 Subject: [PATCH 04/12] fix QA --- pypeit/core/arc.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pypeit/core/arc.py b/pypeit/core/arc.py index 7f5ab21b2d..45b9745907 100644 --- a/pypeit/core/arc.py +++ b/pypeit/core/arc.py @@ -179,7 +179,7 @@ def fit2darc_global_qa(pypeitFit, nspec, outfile=None): # Finish if outfile is not None: - plt.savefig(outfile, dpi=800) + plt.savefig(outfile, dpi=300) plt.close() else: plt.show() @@ -291,9 +291,9 @@ def fit2darc_orders_qa(pypeitFit, nspec, outfile=None): ax1.set_ylabel(r'Res. [pix]') - ax0.text(0.1, 0.9, r'RMS={0:.3f} Pixel'.format(rms_order / np.abs(dwl)), ha="left", va="top", + ax0.text(0.1, 0.8, r'RMS={0:.3f} Pixel'.format(rms_order / np.abs(dwl)), ha="left", va="top", transform=ax0.transAxes) - ax0.text(0.1, 0.8, r'$\Delta\lambda$={0:.3f} Pixel/$\AA$'.format(np.abs(dwl)), ha="left", va="top", + ax0.text(0.1, 0.9, r'$\Delta\lambda$={0:.3f} $\AA$/Pixel'.format(np.abs(dwl)), ha="left", va="top", transform=ax0.transAxes) ax0.get_yaxis().set_label_coords(-0.15, 0.5) @@ -311,7 +311,7 @@ def fit2darc_orders_qa(pypeitFit, nspec, outfile=None): # Finish if outfile is not None: - plt.savefig(outfile, dpi=800) + plt.savefig(outfile, dpi=200) plt.close() else: plt.show() From 861ac9704c7a736b469affd7dba2cbfdbd7378e9 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 13 Sep 2023 14:46:38 -1000 Subject: [PATCH 05/12] fix bug --- pypeit/core/wavecal/autoid.py | 5 +++-- pypeit/wavecalib.py | 2 ++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index b586dffaea..e336f3528e 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -1771,9 +1771,10 @@ def __init__(self, spec, lamps, par=None, ok_mask=None, self._binw = binw self._bind = bind self._measured_fwhms = measured_fwhms + _med_fwhm = np.median(measured_fwhms) if measured_fwhms is not None else None self.rms_thresh = \ - round(self._par['rms_thresh_frac_fwhm'] * set_fwhm(self._par, measured_fwhm=np.median(measured_fwhms))) \ - if rms_thresh is None else rms_thresh + round(self._par['rms_thresh_frac_fwhm'] * set_fwhm(self._par, measured_fwhm=_med_fwhm)) \ + if rms_thresh is None else rms_thresh # Mask info if ok_mask is None: diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 4087ea44ac..3144dc2843 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -620,6 +620,8 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Sometimes works, sometimes fails arcfitter = autoid.HolyGrail(arccen, self.lamps, par=self.par, ok_mask=ok_mask_idx, + rms_thresh=self.wave_rms_thresh, + measured_fwhms=self.measured_fwhms, nonlinear_counts=self.nonlinear_counts, spectrograph=self.spectrograph.name) patt_dict, final_fit = arcfitter.get_results() From 71d36f76367007a37458c6f9fda3b9bc5389a678 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Fri, 15 Sep 2023 09:27:47 -1000 Subject: [PATCH 06/12] changes to the spectrographs --- pypeit/core/wavecal/autoid.py | 28 ++++++++++++++---------- pypeit/spectrographs/bok_bc.py | 4 ++-- pypeit/spectrographs/gemini_flamingos.py | 8 +++---- pypeit/spectrographs/gemini_gmos.py | 3 ++- pypeit/spectrographs/gemini_gnirs.py | 5 +++-- pypeit/spectrographs/keck_esi.py | 3 +-- pypeit/spectrographs/keck_hires.py | 2 +- pypeit/spectrographs/keck_lris.py | 4 ++-- pypeit/spectrographs/keck_mosfire.py | 2 +- pypeit/spectrographs/keck_nires.py | 2 +- pypeit/spectrographs/keck_nirspec.py | 4 ++-- pypeit/spectrographs/lbt_luci.py | 8 +++---- pypeit/spectrographs/lbt_mods.py | 8 +++---- pypeit/spectrographs/ldt_deveny.py | 4 ++-- pypeit/spectrographs/magellan_fire.py | 7 +++--- pypeit/spectrographs/magellan_mage.py | 3 +-- pypeit/spectrographs/mmt_binospec.py | 4 ++-- pypeit/spectrographs/mmt_bluechannel.py | 4 ++-- pypeit/spectrographs/mmt_mmirs.py | 4 ++-- pypeit/spectrographs/ntt_efosc2.py | 2 +- pypeit/spectrographs/p200_tspec.py | 3 +-- pypeit/spectrographs/shane_kast.py | 4 ++-- pypeit/spectrographs/soar_goodman.py | 4 ++-- pypeit/spectrographs/vlt_fors.py | 2 +- pypeit/spectrographs/vlt_sinfoni.py | 2 +- pypeit/spectrographs/vlt_xshooter.py | 9 ++++---- pypeit/wavecalib.py | 5 +++-- 27 files changed, 71 insertions(+), 67 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index e336f3528e..138e53fdb6 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -960,7 +960,7 @@ def measure_fwhm(spec, sigdetect=10., fwhm=5.): return measured_fwhm -def set_fwhm(par, measured_fwhm=None): +def set_fwhm(par, measured_fwhm=None, verbose=False): """ Set the value of the arc lines FWHM by choosing between the provided parset and the measured_fwhm @@ -971,6 +971,8 @@ def set_fwhm(par, measured_fwhm=None): wavelength-solution algorithms. measured_fwhm (:obj:`float`): Measured arc lines FWHM in binned pixels of the input arc image + verbose (:obj:`bool`, optional): + Print a message to screen reporting the chosen FWHM Returns: :obj:`float`: Chosen arc lines FWHM in binned pixels of the input arc image @@ -979,13 +981,16 @@ def set_fwhm(par, measured_fwhm=None): # Set FWHM for the methods that follow if par['fwhm_fromlines'] is False: fwhm = par['fwhm'] - msgs.info(f"User-provided arc lines FWHM: {fwhm:.1f} pixels") + if verbose: + msgs.info(f"User-provided arc lines FWHM: {fwhm:.1f} pixels") elif measured_fwhm is None: fwhm = par['fwhm'] - msgs.warn(f"Assumed arc lines FWHM: {fwhm:.1f} pixels") + if verbose: + msgs.warn(f"Assumed arc lines FWHM: {fwhm:.1f} pixels") else: fwhm = measured_fwhm - msgs.info(f"Measured arc lines FWHM: {fwhm:.1f} pixels") + if verbose: + msgs.info(f"Measured arc lines FWHM: {fwhm:.1f} pixels") return fwhm @@ -1083,7 +1088,7 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, slit_ # Grab the observed arc spectrum obs_spec_i = spec[:,slit] # get FWHM for this slit - fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[slit]) + fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[slit], verbose=True) # Find the shift ncomb = temp_spec.size @@ -1313,7 +1318,7 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, rms_thresh msgs.info("Using sigdetect = {}".format(sigdetect)) msgs.info(f"Using rms_threshold = {rms_thresh} ({par['rms_thresh_frac_fwhm']} * median FWHM)") # Set FWHM for this order - fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[iord]) + fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[iord], verbose=True) detections[str(iord)], spec_cont_sub[:, iord], all_patt_dict[str(iord)] = reidentify( spec[:, iord], spec_arxiv[:, iord], wave_arxiv[:, iord], tot_line_list, par['nreid_min'], cc_thresh=cc_thresh, match_toler=par['match_toler'], @@ -1628,7 +1633,7 @@ def __init__(self, spec, lamps, par, rms_thresh, ech_fixed_format=False, ok_mask msgs.info("Using sigdetect = {}".format(sigdetect)) msgs.info(f"Using rms_threshold = {rms_thresh} ({self.par['rms_thresh_frac_fwhm']} * median FWHM)") # get FWHM for this slit - fwhm = set_fwhm(self.par, measured_fwhm=measured_fwhms[slit]) + fwhm = set_fwhm(self.par, measured_fwhm=measured_fwhms[slit], verbose=True) self.detections[str(slit)], self.spec_cont_sub[:,slit], self.all_patt_dict[str(slit)] = \ reidentify(self.spec[:,slit], self.spec_arxiv[:,ind_sp], self.wave_soln_arxiv[:,ind_sp], @@ -1771,9 +1776,10 @@ def __init__(self, spec, lamps, par=None, ok_mask=None, self._binw = binw self._bind = bind self._measured_fwhms = measured_fwhms - _med_fwhm = np.median(measured_fwhms) if measured_fwhms is not None else None + _med_fwhm = np.median(measured_fwhms[measured_fwhms!=0]) \ + if measured_fwhms is not None and np.any(measured_fwhms!=0) else None self.rms_thresh = \ - round(self._par['rms_thresh_frac_fwhm'] * set_fwhm(self._par, measured_fwhm=_med_fwhm)) \ + round(self._par['rms_thresh_frac_fwhm'] * set_fwhm(self._par, measured_fwhm=_med_fwhm),3) \ if rms_thresh is None else rms_thresh # Mask info @@ -1925,7 +1931,7 @@ def run_brute(self, min_nlines=10): sigdetect = wvutils.parse_param(self._par, 'sigdetect', slit) msgs.info("Using sigdetect = {}".format(sigdetect)) # get FWHM for this slit - fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit]) + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) self._all_tcent, self._all_ecent, self._cut_tcent, self._icut, _ =\ wvutils.arc_lines_from_spec(self._spec[:, slit].copy(), sigdetect=sigdetect, fwhm=fwhm, nonlinear_counts=self._nonlinear_counts) @@ -2035,7 +2041,7 @@ def run_kdtree(self, polygon=4, detsrch=7, lstsrch=10, pixtol=5): self._all_final_fit[str(slit)] = {} continue # get FWHM for this slit - fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit]) + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) # Detect lines, and decide which tcent to use sigdetect = wvutils.parse_param(self._par, 'sigdetect', slit) self._all_tcent, self._all_ecent, self._cut_tcent, self._icut, _ =\ diff --git a/pypeit/spectrographs/bok_bc.py b/pypeit/spectrographs/bok_bc.py index 11a9c642df..63f30da387 100644 --- a/pypeit/spectrographs/bok_bc.py +++ b/pypeit/spectrographs/bok_bc.py @@ -231,9 +231,9 @@ def default_pypeit_par(cls): par['calibrations']['wavelengths']['lamps'] = ['NeI', 'ArI', 'ArII', 'HeI'] # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.19 par['calibrations']['wavelengths']['sigdetect'] = 5. - par['calibrations']['wavelengths']['fwhm']= 5.0 + par['calibrations']['wavelengths']['fwhm']= 2.6 #par['calibrations']['wavelengths']['n_first'] = 3 #par['calibrations']['wavelengths']['n_final'] = 5 diff --git a/pypeit/spectrographs/gemini_flamingos.py b/pypeit/spectrographs/gemini_flamingos.py index 52dd3400b3..6873fb4631 100644 --- a/pypeit/spectrographs/gemini_flamingos.py +++ b/pypeit/spectrographs/gemini_flamingos.py @@ -109,9 +109,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with arc lines - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.1 par['calibrations']['wavelengths']['sigdetect']=5 - par['calibrations']['wavelengths']['fwhm'] = 5 + par['calibrations']['wavelengths']['fwhm'] = 5. par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=4 par['calibrations']['wavelengths']['lamps'] = ['OH_NIRES'] @@ -276,9 +276,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with arc lines - par['calibrations']['wavelengths']['rms_threshold'] = 1.0 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.05 # this needs to be updated par['calibrations']['wavelengths']['sigdetect']=3 - par['calibrations']['wavelengths']['fwhm'] = 20 + par['calibrations']['wavelengths']['fwhm'] = 20 # we don't know this value, no dataset in the repo par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=4 par['calibrations']['wavelengths']['lamps'] = ['ArI', 'ArII', 'ThAr', 'NeI'] diff --git a/pypeit/spectrographs/gemini_gmos.py b/pypeit/spectrographs/gemini_gmos.py index e64a3ca707..3012af71d7 100644 --- a/pypeit/spectrographs/gemini_gmos.py +++ b/pypeit/spectrographs/gemini_gmos.py @@ -263,7 +263,8 @@ def default_pypeit_par(cls): par['calibrations']['slitedges']['fit_order'] = 3 # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.40 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.08 # Might be grating dependent.. + par['calibrations']['wavelengths']['fwhm'] = 5. par['calibrations']['wavelengths']['sigdetect'] = 5. # Doesn't work for reddest chip par['calibrations']['wavelengths']['lamps'] = ['CuI', 'ArI', 'ArII'] par['calibrations']['wavelengths']['method'] = 'full_template' diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py index eca3e77cc9..920eed62b5 100644 --- a/pypeit/spectrographs/gemini_gnirs.py +++ b/pypeit/spectrographs/gemini_gnirs.py @@ -322,7 +322,8 @@ def config_specific_par(self, scifile, inp_par=None): par['calibrations']['slitedges']['fit_min_spec_length'] = 0.5 # Wavelengths - par['calibrations']['wavelengths']['rms_threshold'] = 1.0 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.4 # Might be grating dependent.. + par['calibrations']['wavelengths']['fwhm'] = 2.5 par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['lamps'] = ['OH_GNIRS'] # par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] @@ -354,7 +355,7 @@ def config_specific_par(self, scifile, inp_par=None): par['calibrations']['slitedges']['sync_predict'] = 'nearest' # Wavelengths - par['calibrations']['wavelengths']['rms_threshold'] = 1.0 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 1.0 # Might be grating dependent.. par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['lamps'] = ['Ar_IR_GNIRS'] # par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] diff --git a/pypeit/spectrographs/keck_esi.py b/pypeit/spectrographs/keck_esi.py index a8fbd33ea1..c7292ceb5b 100644 --- a/pypeit/spectrographs/keck_esi.py +++ b/pypeit/spectrographs/keck_esi.py @@ -91,9 +91,8 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution # This is for 1x1 - par['calibrations']['wavelengths']['rms_threshold'] = 0.30 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.103 par['calibrations']['wavelengths']['fwhm'] = 2.9 - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['lamps'] = ['CuI', 'ArI', 'NeI', 'HgI', 'XeI', 'ArII'] diff --git a/pypeit/spectrographs/keck_hires.py b/pypeit/spectrographs/keck_hires.py index df4c8eab03..a3d325321f 100644 --- a/pypeit/spectrographs/keck_hires.py +++ b/pypeit/spectrographs/keck_hires.py @@ -115,7 +115,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['ThAr'] # This is for 1x1 binning. TODO GET BINNING SORTED OUT!! - par['calibrations']['wavelengths']['rms_threshold'] = 0.50 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.1 par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['n_final'] = 4 #[3] + 13 * [4] + [3] # This is for 1x1 binning. Needs to be divided by binning for binned data!! diff --git a/pypeit/spectrographs/keck_lris.py b/pypeit/spectrographs/keck_lris.py index 11a4c93810..96c51874a3 100644 --- a/pypeit/spectrographs/keck_lris.py +++ b/pypeit/spectrographs/keck_lris.py @@ -75,7 +75,7 @@ def default_pypeit_par(cls): # Remove slits that are too short par['calibrations']['slitedges']['minimum_slit_length'] = 3. # 1D wavelengths - par['calibrations']['wavelengths']['rms_threshold'] = 0.20 # Might be grism dependent + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.05 # Might be grism dependent # Set the default exposure time ranges for the frame typing par['calibrations']['biasframe']['exprng'] = [None, 0.001] par['calibrations']['darkframe']['exprng'] = [999999, None] # No dark frames @@ -846,7 +846,7 @@ def default_pypeit_par(cls): par['calibrations']['slitedges']['fit_min_spec_length'] = 0.2 # 1D wavelength solution -- Additional parameters are grism dependent - par['calibrations']['wavelengths']['rms_threshold'] = 0.20 # Might be grism dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.06 # Might be grism dependent.. par['calibrations']['wavelengths']['sigdetect'] = 10.0 #par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] diff --git a/pypeit/spectrographs/keck_mosfire.py b/pypeit/spectrographs/keck_mosfire.py index fe78ec576e..177a74134e 100644 --- a/pypeit/spectrographs/keck_mosfire.py +++ b/pypeit/spectrographs/keck_mosfire.py @@ -84,7 +84,7 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.30 #0.20 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.11 #0.20 # Might be grating dependent.. par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 5.0 par['calibrations']['wavelengths']['n_final']= 4 diff --git a/pypeit/spectrographs/keck_nires.py b/pypeit/spectrographs/keck_nires.py index e4d5ef1747..af527ccfb1 100644 --- a/pypeit/spectrographs/keck_nires.py +++ b/pypeit/spectrographs/keck_nires.py @@ -80,7 +80,7 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.30 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.136 par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 2.2 # Measured par['calibrations']['wavelengths']['fwhm_fromlines'] = True diff --git a/pypeit/spectrographs/keck_nirspec.py b/pypeit/spectrographs/keck_nirspec.py index 07ac4fa3d7..eab7cafc91 100644 --- a/pypeit/spectrographs/keck_nirspec.py +++ b/pypeit/spectrographs/keck_nirspec.py @@ -70,9 +70,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.20 #0.20 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.06 #0.20 # Might be grating dependent.. par['calibrations']['wavelengths']['sigdetect']=5.0 - par['calibrations']['wavelengths']['fwhm']= 5.0 + par['calibrations']['wavelengths']['fwhm']= 3.5 par['calibrations']['wavelengths']['n_final']= 4 par['calibrations']['wavelengths']['lamps'] = ['OH_NIRES'] #par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] diff --git a/pypeit/spectrographs/lbt_luci.py b/pypeit/spectrographs/lbt_luci.py index 0ac9d14a04..1406486616 100644 --- a/pypeit/spectrographs/lbt_luci.py +++ b/pypeit/spectrographs/lbt_luci.py @@ -329,10 +329,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths'][ - 'rms_threshold'] = 0.20 # 0.20 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.05 # 0.20 # Might be grating dependent.. par['calibrations']['wavelengths']['sigdetect'] = 5.0 - par['calibrations']['wavelengths']['fwhm'] = 5.0 + par['calibrations']['wavelengths']['fwhm'] = 4.5 par['calibrations']['wavelengths']['n_final'] = 4 par['calibrations']['wavelengths']['lamps'] = ['OH_NIRES'] #par['calibrations']['wavelengths']['nonlinear_counts'] = \ @@ -444,8 +443,7 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths'][ - 'rms_threshold'] = 0.20 # 0.20 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.04 # 0.20 # Might be grating dependent.. par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['fwhm'] = 5.0 par['calibrations']['wavelengths']['n_final'] = 4 diff --git a/pypeit/spectrographs/lbt_mods.py b/pypeit/spectrographs/lbt_mods.py index 026d727e2e..e893f00b7c 100644 --- a/pypeit/spectrographs/lbt_mods.py +++ b/pypeit/spectrographs/lbt_mods.py @@ -327,7 +327,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['sigdetect'] = 5. - par['calibrations']['wavelengths']['rms_threshold'] = 0.4 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.09 par['calibrations']['wavelengths']['fwhm'] = 10. #par['calibrations']['wavelengths']['lamps'] = ['XeI','ArII','ArI','NeI','KrI']] par['calibrations']['wavelengths']['lamps'] = ['ArI','NeI','KrI','XeI'] @@ -491,7 +491,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['sigdetect'] = 10. - par['calibrations']['wavelengths']['rms_threshold'] = 0.4 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.09 par['calibrations']['wavelengths']['lamps'] = ['XeI','KrI','ArI','HgI'] # slit @@ -648,7 +648,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['sigdetect'] = 5. - par['calibrations']['wavelengths']['rms_threshold'] = 1.0 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.22 par['calibrations']['wavelengths']['fwhm'] = 10. #par['calibrations']['wavelengths']['lamps'] = ['XeI','ArII','ArI','NeI','KrI']] par['calibrations']['wavelengths']['lamps'] = ['ArI','NeI','KrI','XeI'] @@ -810,7 +810,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['sigdetect'] = 10. - par['calibrations']['wavelengths']['rms_threshold'] = 0.4 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.09 par['calibrations']['wavelengths']['lamps'] = ['XeI','KrI','ArI','HgI'] # slit diff --git a/pypeit/spectrographs/ldt_deveny.py b/pypeit/spectrographs/ldt_deveny.py index df2b44cc4e..2d5ce4d884 100644 --- a/pypeit/spectrographs/ldt_deveny.py +++ b/pypeit/spectrographs/ldt_deveny.py @@ -492,7 +492,7 @@ def config_specific_par(self, scifile, inp_par=None): # and it's associated tweaks in parameters par['calibrations']['wavelengths']['method'] = 'holy-grail' par['calibrations']['wavelengths']['sigdetect'] = 10.0 # Default: 5.0 - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 # Default: 0.15 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.15 # Default: 0.15 # Start with a lower-order Legendre polymonial for the wavelength fit par['calibrations']['wavelengths']['n_first'] = 2 # Default: 3 # The approximate resolution of this grating @@ -527,7 +527,7 @@ def config_specific_par(self, scifile, inp_par=None): # and it's associated tweaks in parameters par['calibrations']['wavelengths']['method'] = 'holy-grail' par['calibrations']['wavelengths']['sigdetect'] = 10.0 # Default: 5.0 - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 # Default: 0.15 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.25 # Default: 0.15 # Start/end with a lower-order Legendre polymonial for the wavelength fit par['calibrations']['wavelengths']['n_first'] = 2 # Default: 3 par['calibrations']['wavelengths']['n_final'] = 4 # Default: 5 diff --git a/pypeit/spectrographs/magellan_fire.py b/pypeit/spectrographs/magellan_fire.py index 1f103d8cba..7172708e2e 100644 --- a/pypeit/spectrographs/magellan_fire.py +++ b/pypeit/spectrographs/magellan_fire.py @@ -145,7 +145,8 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with OH lines - par['calibrations']['wavelengths']['rms_threshold'] = 1.0 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.2 + par['calibrations']['wavelengths']['fwhm'] = 5. par['calibrations']['wavelengths']['sigdetect']=[5,10,10,10,10,20,30,30,30,30,30,10,30,30,60,30,30,10,20,30,10] par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=[3,3,3,2,4,4,4,3,4,4,4,3,4,4,4,4,4,4,6,6,4] @@ -391,9 +392,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with arc lines - par['calibrations']['wavelengths']['rms_threshold'] = 1.0 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.1 par['calibrations']['wavelengths']['sigdetect']=3 - par['calibrations']['wavelengths']['fwhm'] = 20 + par['calibrations']['wavelengths']['fwhm'] = 10 par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=4 par['calibrations']['wavelengths']['lamps'] = ['ArI', 'ArII', 'ThAr', 'NeI'] diff --git a/pypeit/spectrographs/magellan_mage.py b/pypeit/spectrographs/magellan_mage.py index fc8fdc9d1f..731aa8197f 100644 --- a/pypeit/spectrographs/magellan_mage.py +++ b/pypeit/spectrographs/magellan_mage.py @@ -95,9 +95,8 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution # The following is for 1x1 binning - par['calibrations']['wavelengths']['rms_threshold'] = 0.40 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.133 par['calibrations']['wavelengths']['fwhm'] = 3.0 - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['lamps'] = ['ThAr_MagE'] diff --git a/pypeit/spectrographs/mmt_binospec.py b/pypeit/spectrographs/mmt_binospec.py index 5dfb82155a..7b926fd763 100644 --- a/pypeit/spectrographs/mmt_binospec.py +++ b/pypeit/spectrographs/mmt_binospec.py @@ -176,9 +176,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.125 par['calibrations']['wavelengths']['sigdetect'] = 5. - par['calibrations']['wavelengths']['fwhm']= 5.0 + par['calibrations']['wavelengths']['fwhm']= 4.0 par['calibrations']['wavelengths']['lamps'] = ['ArI', 'ArII'] #par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] par['calibrations']['wavelengths']['method'] = 'full_template' diff --git a/pypeit/spectrographs/mmt_bluechannel.py b/pypeit/spectrographs/mmt_bluechannel.py index f6257758da..0e6d742234 100644 --- a/pypeit/spectrographs/mmt_bluechannel.py +++ b/pypeit/spectrographs/mmt_bluechannel.py @@ -227,9 +227,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.15 + par['calibrations']['wavelengths']['fwhm'] = 3.1 par['calibrations']['wavelengths']['sigdetect'] = 5. - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # Parse the lamps used from the image header. par['calibrations']['wavelengths']['lamps'] = ['use_header'] diff --git a/pypeit/spectrographs/mmt_mmirs.py b/pypeit/spectrographs/mmt_mmirs.py index 36ae51c027..a341fb42f3 100644 --- a/pypeit/spectrographs/mmt_mmirs.py +++ b/pypeit/spectrographs/mmt_mmirs.py @@ -158,9 +158,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with arc lines - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.125 par['calibrations']['wavelengths']['sigdetect']=5 - par['calibrations']['wavelengths']['fwhm'] = 5 + par['calibrations']['wavelengths']['fwhm'] = 4. par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=4 par['calibrations']['wavelengths']['lamps'] = ['OH_NIRES'] diff --git a/pypeit/spectrographs/ntt_efosc2.py b/pypeit/spectrographs/ntt_efosc2.py index 71c425c7db..e6355c3624 100644 --- a/pypeit/spectrographs/ntt_efosc2.py +++ b/pypeit/spectrographs/ntt_efosc2.py @@ -253,7 +253,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['method'] = 'full_template' par['calibrations']['wavelengths']['lamps'] = ['HeI', 'ArI'] - par['calibrations']['wavelengths']['rms_threshold'] = 0.25 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.07 par['calibrations']['wavelengths']['sigdetect'] = 10.0 par['calibrations']['wavelengths']['fwhm'] = 4.0 par['calibrations']['wavelengths']['n_final'] = 4 diff --git a/pypeit/spectrographs/p200_tspec.py b/pypeit/spectrographs/p200_tspec.py index 28da608420..1d4000d047 100644 --- a/pypeit/spectrographs/p200_tspec.py +++ b/pypeit/spectrographs/p200_tspec.py @@ -158,10 +158,9 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.3 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.103 par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 2.9 # As measured in DevSuite - par['calibrations']['wavelengths']['fwhm_fromlines'] = True par['calibrations']['wavelengths']['n_final']= [3,4,4,4,4] par['calibrations']['wavelengths']['lamps'] = ['OH_NIRES'] par['calibrations']['wavelengths']['method'] = 'reidentify' diff --git a/pypeit/spectrographs/shane_kast.py b/pypeit/spectrographs/shane_kast.py index 45d86f1282..e58691757f 100644 --- a/pypeit/spectrographs/shane_kast.py +++ b/pypeit/spectrographs/shane_kast.py @@ -259,7 +259,7 @@ def default_pypeit_par(cls): par['flexure']['spectrum'] = 'sky_kastb_600.fits' # 1D wavelength solution par['calibrations']['wavelengths']['sigdetect'] = 5. - par['calibrations']['wavelengths']['rms_threshold'] = 0.20 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.07 par['calibrations']['wavelengths']['lamps'] = ['CdI','HgI','HeI'] par['calibrations']['wavelengths']['method'] = 'full_template' @@ -681,7 +681,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['NeI', 'HgI', 'HeI', 'ArI'] - par['calibrations']['wavelengths']['rms_threshold'] = 0.20 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.09 par['calibrations']['wavelengths']['sigdetect'] = 5. par['calibrations']['wavelengths']['use_instr_flag'] = True par['sensfunc']['IR']['telgridfile'] = 'TelFit_Lick_3100_11100_R10000.fits' diff --git a/pypeit/spectrographs/soar_goodman.py b/pypeit/spectrographs/soar_goodman.py index 1c523674e5..1f1f9c5bc4 100644 --- a/pypeit/spectrographs/soar_goodman.py +++ b/pypeit/spectrographs/soar_goodman.py @@ -310,7 +310,7 @@ def default_pypeit_par(cls): par['calibrations']['wavelengths']['lamps'] = ['NeI', 'ArI', 'HgI'] # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.17 par['calibrations']['wavelengths']['sigdetect'] = 5. par['calibrations']['wavelengths']['fwhm']= 5.0 par['calibrations']['flatfield']['slit_illum_finecorr'] = False @@ -508,7 +508,7 @@ def default_pypeit_par(cls): par['calibrations']['wavelengths']['lamps'] = ['NeI', 'ArI', 'HgI'] # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.17 par['calibrations']['wavelengths']['sigdetect'] = 5. par['calibrations']['wavelengths']['fwhm'] = 5.0 diff --git a/pypeit/spectrographs/vlt_fors.py b/pypeit/spectrographs/vlt_fors.py index 54e4a9772f..cdbc0d2cc6 100644 --- a/pypeit/spectrographs/vlt_fors.py +++ b/pypeit/spectrographs/vlt_fors.py @@ -57,7 +57,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['HeI', 'ArI'] # Grating dependent - par['calibrations']['wavelengths']['rms_threshold'] = 0.25 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.07 par['calibrations']['wavelengths']['sigdetect'] = 10.0 par['calibrations']['wavelengths']['fwhm'] = 4.0 # Good for 2x binning par['calibrations']['wavelengths']['n_final'] = 4 diff --git a/pypeit/spectrographs/vlt_sinfoni.py b/pypeit/spectrographs/vlt_sinfoni.py index 5bfbcf855c..42ba6cf71c 100644 --- a/pypeit/spectrographs/vlt_sinfoni.py +++ b/pypeit/spectrographs/vlt_sinfoni.py @@ -76,7 +76,7 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_threshold'] = 0.30 #0.20 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.23 #0.20 # Might be grating dependent.. par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 5.0 par['calibrations']['wavelengths']['n_final']= 4 diff --git a/pypeit/spectrographs/vlt_xshooter.py b/pypeit/spectrographs/vlt_xshooter.py index c9f7aa9a42..7ec5db8c21 100644 --- a/pypeit/spectrographs/vlt_xshooter.py +++ b/pypeit/spectrographs/vlt_xshooter.py @@ -295,9 +295,9 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['OH_XSHOOTER'] - par['calibrations']['wavelengths']['rms_threshold'] = 0.60 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.15 par['calibrations']['wavelengths']['sigdetect'] = 10.0 - par['calibrations']['wavelengths']['fwhm'] = 5.0 + par['calibrations']['wavelengths']['fwhm'] = 4. par['calibrations']['wavelengths']['n_final'] = 4 # Reidentification parameters par['calibrations']['wavelengths']['method'] = 'reidentify' @@ -682,9 +682,8 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['ThAr_XSHOOTER_VIS'] # The following is for 1x1 binning. TODO GET BINNING SORTED OUT!! - par['calibrations']['wavelengths']['rms_threshold'] = 1.2 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.15 par['calibrations']['wavelengths']['fwhm'] = 8.0 - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['n_final'] = [3] + 13*[4] + [3] @@ -955,7 +954,7 @@ def default_pypeit_par(cls): par['calibrations']['wavelengths']['lamps'] = ['ThAr_XSHOOTER_UVB'] par['calibrations']['wavelengths']['n_final'] = [3] + 10*[4] # This is for 1x1 - par['calibrations']['wavelengths']['rms_threshold'] = 0.70 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.184 par['calibrations']['wavelengths']['fwhm'] = 3.8 par['calibrations']['wavelengths']['fwhm_fromlines'] = True # diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 3144dc2843..e760108f3c 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -612,7 +612,8 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Save for redo's self.measured_fwhms = measured_fwhms # determine rms threshold - fwhm_for_thresh = autoid.set_fwhm(self.par, measured_fwhm=np.median(measured_fwhms)) + _med_fwhm = np.median(measured_fwhms[measured_fwhms!=0]) if np.any(measured_fwhms!=0) else None + fwhm_for_thresh = autoid.set_fwhm(self.par, measured_fwhm=_med_fwhm) self.wave_rms_thresh = round(self.par['rms_thresh_frac_fwhm'] * fwhm_for_thresh,3) # Obtain calibration for all slits @@ -808,7 +809,7 @@ def redo_echelle_orders(self, bad_orders:np.ndarray, dets:np.ndarray, order_dets wv_order_mod = self.wv_calib.wv_fit2d[idet].eval(spec_vec_norm, x2=np.ones_like(spec_vec_norm)*order)/order # get FWHM for this order - fwhm = autoid.set_fwhm(self.par, measured_fwhm=self.measured_fwhms[iord]) + fwhm = autoid.set_fwhm(self.par, measured_fwhm=self.measured_fwhms[iord], verbose=True) # Link me tcent, spec_cont_sub, patt_dict_slit, tot_llist = autoid.match_to_arxiv( self.lamps, self.arccen[:,iord], wv_order_mod, From 12eb635f169df72976589ab10c272982226a51b3 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Tue, 19 Sep 2023 17:05:06 -0700 Subject: [PATCH 07/12] comments --- pypeit/core/wavecal/autoid.py | 494 ++++++++++++++------------ pypeit/spectrographs/gemini_gnirs.py | 5 +- pypeit/spectrographs/magellan_fire.py | 3 +- pypeit/wavecalib.py | 30 +- 4 files changed, 283 insertions(+), 249 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index 138e53fdb6..fea9b87f00 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -970,7 +970,8 @@ def set_fwhm(par, measured_fwhm=None, verbose=False): Key parameters that drive the behavior of the wavelength-solution algorithms. measured_fwhm (:obj:`float`): - Measured arc lines FWHM in binned pixels of the input arc image + Measured arc lines FWHM in binned pixels of the input arc image. + If None, the value provided by the user in the `fwhm` parset is used. verbose (:obj:`bool`, optional): Print a message to screen reporting the chosen FWHM @@ -1031,7 +1032,8 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, slit_ slit_ids: ndarray, optional Array of slit/order IDs. Shape (nslit,) measured_fwhms : `numpy.ndarray`_, optional - Array of FWHM (in pixels) measured from the arc lines. Shape (nslit,) + Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,). + If None, the value provided by the user in the `fwhm` parset is used. x_percentile : float, optional Passed to reidentify to reduce the dynamic range of arc line amplitudes template_dict : dict, optional @@ -1201,7 +1203,7 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, slit_ # Finish return wvcalib -def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, rms_thresh, +def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, ok_mask=None, measured_fwhms=None, use_unknowns=True, debug_all=False, debug_peaks=False, debug_xcorr=False, debug_reid=False, debug_fits=False, nonlinear_counts=1e10, @@ -1228,13 +1230,12 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, rms_thresh par : :class:`~pypeit.par.pypeitpar.WavelengthSolutionPar` Key parameters that drive the behavior of the wavelength-solution algorithms. - rms_thresh: float - RMS threshold for the wavelength solution fit ok_mask : `numpy.ndarray`, optional Integer array with the list of valid spectra ``spec`` to use. If None, all spectra are used. measured_fwhms: ndarray, optional - Array of FWHM (in binned pixels) measured from the arc lines. Shape :math:`(N_{\rm orders},)` + Array of FWHM (in binned pixels) measured from the arc lines. Shape :math:`(N_{\rm orders},)`. + If None, the value provided by the user in the `fwhm` parset is used. use_unknowns : bool, default = True, optional If True, arc lines that are known to be present in the spectra, but have not been attributed to an element+ion, will @@ -1316,9 +1317,11 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, rms_thresh sigdetect = wvutils.parse_param(par, 'sigdetect', iord) cc_thresh = wvutils.parse_param(par, 'cc_thresh', iord) msgs.info("Using sigdetect = {}".format(sigdetect)) - msgs.info(f"Using rms_threshold = {rms_thresh} ({par['rms_thresh_frac_fwhm']} * median FWHM)") # Set FWHM for this order fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[iord], verbose=True) + # get rms threshold for this slit + rms_thresh = round(par['rms_thresh_frac_fwhm'] * fwhm, 3) + msgs.info(f"Using rms_threshold = {rms_thresh} ({par['rms_thresh_frac_fwhm']} * median FWHM)") detections[str(iord)], spec_cont_sub[:, iord], all_patt_dict[str(iord)] = reidentify( spec[:, iord], spec_arxiv[:, iord], wave_arxiv[:, iord], tot_line_list, par['nreid_min'], cc_thresh=cc_thresh, match_toler=par['match_toler'], @@ -1455,8 +1458,6 @@ class ArchiveReid: par : :class:`~pypeit.par.pypeitpar.WavelengthSolutionPar` Key parameters that drive the behavior of the wavelength-solution algorithms. - rms_thresh: float - RMS threshold for the wavelength solution fit ech_fixed_format: bool Set to True if this is a fixed format echelle spectrograp. The code will then align the archive_arc and the extracted arc for each order for the reidentification. @@ -1464,7 +1465,8 @@ class ArchiveReid: Integer array with the list of valid spectra ``spec`` to use. If None, all spectra are used. measured_fwhms: ndarray, optional - Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,) + Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,). + If None, the value provided by the user in the `fwhm` parset is used. use_unknowns : bool, default = True, optional If True, arc lines that are known to be present in the spectra, but have not been attributed to an element+ion, will @@ -1505,7 +1507,7 @@ class ArchiveReid: """ # TODO: Because we're passing orders directly, we no longer need spectrograph... - def __init__(self, spec, lamps, par, rms_thresh, ech_fixed_format=False, ok_mask=None, + def __init__(self, spec, lamps, par, ech_fixed_format=False, ok_mask=None, measured_fwhms=None, use_unknowns=True, debug_all=False, debug_peaks=False, debug_xcorr=False, debug_reid=False, debug_fits=False, orders=None, nonlinear_counts=1e10): @@ -1631,9 +1633,11 @@ def __init__(self, spec, lamps, par, rms_thresh, ech_fixed_format=False, ok_mask sigdetect = wvutils.parse_param(self.par, 'sigdetect', slit) cc_thresh = wvutils.parse_param(self.par, 'cc_thresh', slit) msgs.info("Using sigdetect = {}".format(sigdetect)) - msgs.info(f"Using rms_threshold = {rms_thresh} ({self.par['rms_thresh_frac_fwhm']} * median FWHM)") # get FWHM for this slit fwhm = set_fwhm(self.par, measured_fwhm=measured_fwhms[slit], verbose=True) + # get rms threshold for this slit + rms_thresh = round(par['rms_thresh_frac_fwhm'] * fwhm, 3) + msgs.info(f"Using rms_threshold = {rms_thresh} ({self.par['rms_thresh_frac_fwhm']} * median FWHM)") self.detections[str(slit)], self.spec_cont_sub[:,slit], self.all_patt_dict[str(slit)] = \ reidentify(self.spec[:,slit], self.spec_arxiv[:,ind_sp], self.wave_soln_arxiv[:,ind_sp], @@ -1724,10 +1728,9 @@ class HolyGrail: islinelist : bool, optional Is lamps a linelist (True), or a list of ions (False) The former is not recommended except by expert users/developers - rms_thresh : float, optional - RMS threshold for the wavelength solution fit measured_fwhms : ndarray, optional - Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,) + Array of FWHM (in binned pixels) measured from the arc lines. Shape (nslit,). + If None, the value provided by the user in the `fwhm` parset is used. outroot : str, optional Name of output file debug : bool, optional @@ -1762,7 +1765,7 @@ class HolyGrail: """ def __init__(self, spec, lamps, par=None, ok_mask=None, - islinelist=False, rms_thresh=None, measured_fwhms=None, + islinelist=False, measured_fwhms=None, outroot=None, debug=False, verbose=False, binw=None, bind=None, nstore=1, use_unknowns=True, nonlinear_counts=None, spectrograph=None): @@ -1776,11 +1779,6 @@ def __init__(self, spec, lamps, par=None, ok_mask=None, self._binw = binw self._bind = bind self._measured_fwhms = measured_fwhms - _med_fwhm = np.median(measured_fwhms[measured_fwhms!=0]) \ - if measured_fwhms is not None and np.any(measured_fwhms!=0) else None - self.rms_thresh = \ - round(self._par['rms_thresh_frac_fwhm'] * set_fwhm(self._par, measured_fwhm=_med_fwhm),3) \ - if rms_thresh is None else rms_thresh # Mask info if ok_mask is None: @@ -1866,7 +1864,18 @@ def set_grids(self, ngridw = 300, ngridd=3000): #ngridw = 200, ngridd=2000): self._ngridd = self._bind.size return - def run_brute_loop(self, slit, tcent_ecent, wavedata=None): + def run_brute_loop(self, slit, tcent_ecent, rms_thresh, wavedata=None): + """ + + Args: + slit (int): + tcent_ecent: + rms_thresh: + wavedata: + + Returns: + + """ # Set the parameter space that gets searched rng_poly = [3, 4] # Range of algorithms to check (only trigons+tetragons are supported) rng_list = range(3, 6) # Number of lines to search over for the linelist @@ -1875,7 +1884,7 @@ def run_brute_loop(self, slit, tcent_ecent, wavedata=None): idthresh = 0.5 # Criteria for early return (at least this fraction of lines must have # an ID on either side of the spectrum) - msgs.info(f"Using rms_threshold = {self.rms_thresh} ({self._par['rms_thresh_frac_fwhm']} * median FWHM)") + msgs.info(f"Using rms_threshold = {rms_thresh} ({self._par['rms_thresh_frac_fwhm']} * median FWHM)") best_patt_dict, best_final_fit = None, None # Loop through parameter space for poly in rng_poly: @@ -1895,7 +1904,7 @@ def run_brute_loop(self, slit, tcent_ecent, wavedata=None): # First time a fit is found best_patt_dict, best_final_fit = copy.deepcopy(patt_dict), copy.deepcopy(final_fit) continue - elif final_fit['rms'] < self.rms_thresh: + elif final_fit['rms'] < rms_thresh: # Has a better fit been identified (i.e. more lines identified)? if len(final_fit['pixel_fit']) > len(best_final_fit['pixel_fit']): best_patt_dict, best_final_fit = copy.deepcopy(patt_dict), copy.deepcopy(final_fit) @@ -1932,6 +1941,8 @@ def run_brute(self, min_nlines=10): msgs.info("Using sigdetect = {}".format(sigdetect)) # get FWHM for this slit fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + # get rms threshold for this slit + rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) self._all_tcent, self._all_ecent, self._cut_tcent, self._icut, _ =\ wvutils.arc_lines_from_spec(self._spec[:, slit].copy(), sigdetect=sigdetect, fwhm=fwhm, nonlinear_counts=self._nonlinear_counts) @@ -1955,7 +1966,7 @@ def run_brute(self, min_nlines=10): self._det_stro[str(slit)] = [self._all_tcent[self._icut].copy(),self._all_ecent[self._icut].copy()] # Run brute force algorithm on the weak lines - best_patt_dict, best_final_fit = self.run_brute_loop(slit, self._det_weak[str(slit)]) + best_patt_dict, best_final_fit = self.run_brute_loop(slit, self._det_weak[str(slit)], rms_thresh) # Print preliminary report good_fit[slit] = self.report_prelim(slit, best_patt_dict, best_final_fit) @@ -2363,7 +2374,12 @@ def cross_match(self, good_fit, detections): new_bad_slits = np.append(new_bad_slits, bs) continue - if final_fit['rms'] > self.rms_thresh: + # get FWHM for this slit + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + # get rms threshold for this slit + rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) + + if final_fit['rms'] > rms_thresh: msgs.warn('---------------------------------------------------' + msgs.newline() + 'Cross-match report for slit {0:d}/{1:d}:'.format(bs + 1, self._nslit) + msgs.newline() + ' Poor RMS ({0:.3f})! Will try cross matching iteratively'.format(final_fit['rms']) + msgs.newline() + @@ -2381,211 +2397,212 @@ def cross_match(self, good_fit, detections): plt.show() return new_bad_slits - def cross_match_order(self, good_fit): - """Using the solutions of all orders, identify the good solutions, and refit the bad ones! - - TODO: This function needs work... The first few lines of code successfully pick up the good orders, - but we need a new routine that (based on an estimated central wavelength and dispersion) can successfully - ID all of the lines. - """ - # DEPRECATED (NOT USED) - - # First determine the central wavelength and dispersion of every slit, using the known good solutions - xplt = np.arange(self._nslit) - yplt, dplt = np.zeros(self._nslit), np.zeros(self._nslit) - imsk = np.ones(self._nslit, dtype=int) - for slit in range(self._nslit): - if good_fit[slit]: - yplt[slit] = self._all_patt_dict[str(slit)]['bwv'] - dplt[slit] = self._all_patt_dict[str(slit)]['bdisp'] - imsk[slit] = 0 - - mask, fit = utils.robust_polyfit(xplt, yplt, 2, function='polynomial', sigma=2, - initialmask=imsk, forceimask=True) - good_fit[mask == 1] = False - wavemodel = utils.func_val(fit, xplt, 'polynomial') - disp = np.median(dplt[good_fit]) - - # TODO: maybe rethink the model at this point? Using the derived - # central wavelength and dispersion identify liens in all orders? - - if self._debug: - plt.subplot(211) - plt.plot(xplt, wavemodel, 'r-') - ww = np.where(mask==0) - plt.plot(xplt[ww], yplt[ww], 'bx') - ww = np.where(mask==1) - plt.plot(xplt[ww], yplt[ww], 'rx') - plt.subplot(212) - plt.plot(xplt, dplt, 'bx') - plt.show() - #embed() - - fact_nl = 1.2 # Non linear factor - new_good_fit = np.zeros(self._nslit, dtype=bool) - for slit in range(self._nslit): - wmin = wavemodel[slit] - fact_nl*disp*self._npix/2 - wmax = wavemodel[slit] + fact_nl*disp*self._npix/2 - ww = np.where((self._wvdata > wmin) & (self._wvdata < wmax)) - wavedata = self._wvdata[ww] - msgs.info('Brute force ID for slit {0:d}/{1:d}'.format(slit+1, self._nslit)) - best_patt_dict, best_final_fit =\ - self.run_brute_loop(slit, arrerr=self._det_weak[str(slit)], wavedata=wavedata) - - self._all_patt_dict[str(slit)] = copy.deepcopy(best_patt_dict) - self._all_final_fit[str(slit)] = copy.deepcopy(best_final_fit) - new_good_fit[slit] = self.report_prelim(slit, best_patt_dict, best_final_fit) - return new_good_fit - - - # Set some fitting parameters - if self._n_final is None: - order = 4 - else: - order = self._n_final - - ofit = [5, 3, 1, 0] - lnpc = len(ofit) - 1 - - # Prepare the fitting coefficients - xv = np.arange(self._npix)/(self._npix-1) - ords = np.arange(self._nslit) - xcen = xv[:, np.newaxis].repeat(self._nslit, axis=1) - extrapord = ~good_fit - maskord = np.where(extrapord)[0] - - coeffs = None - waves = np.zeros(xcen.shape, dtype=float) - for slit in range(self._nslit): - if good_fit[slit]: - func = self._all_final_fit[str(slit)]['function'] - fmin = self._all_final_fit[str(slit)]['fmin'] - fmax = self._all_final_fit[str(slit)]['fmax'] - fitc = self._all_final_fit[str(slit)]['fitc'] - if coeffs is None: - coeffs = np.zeros((fitc.size, self._nslit)) - coeffs[:, slit] = fitc.copy() - waves[:, slit] = utils.func_val(fitc, xv, func, minx=fmin, maxx=fmax) - - msgs.info("Performing a PCA on the order wavelength solutions") - #embed() - pca_wave, outpar = pca.basis(xcen, waves, coeffs, lnpc, ofit, x0in=ords, mask=maskord, skipx0=False, function=func) - - # Report the QA - # TODO: fix setup passing - setup = "BLAH" - pca.pca_plot(setup, outpar, ofit, "wave_cross_match", pcadesc="Wavelength calibration PCA") - - - # Extrapolate the remaining orders requested - #extrap_wave, outpar = pca.extrapolate(outpar, ords) - - # Determine if pixels correlate and anti-correlate with wavelength - signs = np.zeros(self._nslit, dtype=int) - for slit in range(self._nslit): - wvval = pca_wave[:, slit] - if wvval[wvval.size//2] > wvval[wvval.size//2-1]: - signs[slit] = 1 - else: - signs[slit] = -1 - sign = 1 - if np.sum(signs) < 0: - sign = -1 - - new_bad_slits = np.array([], dtype=int) - # Using the first guesses at the wavelength solution, identify lines - for slit in range(self._nslit): - # Get the detections - dets, _ = self.get_use_tcent(sign, self._det_weak[str(slit)]) - lindex = np.array([], dtype=int) - dindex = np.array([], dtype=int) - # Calculate wavelengths for the gsdet detections - wvval = pca_wave[:, slit] - wvcen = wvval[wvval.size//2] - disp = abs(wvval[wvval.size//2] - wvval[wvval.size//2-1]) - for dd in range(dets.size): - pdiff = np.abs(dets[dd] - xv) - bstpx = np.argmin(pdiff) - bstwv = np.abs(self._wvdata - wvval[bstpx]) - if bstwv[np.argmin(bstwv)] > 10.0 * disp: - # This is probably not a good match - continue - lindex = np.append(lindex, np.argmin(bstwv)) - dindex = np.append(dindex, dd) - - # Finalize the best guess of each line - # Initialise the patterns dictionary - patt_dict = dict(acceptable=False, nmatch=0, ibest=-1, bwv=0., - sigdetect=wvutils.parse_param(self._par, 'sigdetect', slit), - mask=np.zeros(dets.size, dtype=bool)) - patt_dict['sign'] = sign - patt_dict['bwv'] = wvcen - patt_dict['bdisp'] = disp - - patterns.solve_triangles(dets, self._wvdata, dindex, lindex, patt_dict) - # Check if a solution was found - if not patt_dict['acceptable']: - new_bad_slits = np.append(new_bad_slits, slit) - msgs.warn('---------------------------------------------------' + msgs.newline() + - 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + - ' Lines could not be identified! Will try cross matching iteratively' + msgs.newline() + - '---------------------------------------------------') - continue - final_fit = self.fit_slit(slit, patt_dict, dets) - if final_fit is None: - # This pattern wasn't good enough - new_bad_slits = np.append(new_bad_slits, slit) - msgs.warn('---------------------------------------------------' + msgs.newline() + - 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + - ' Fit was not good enough! Will try cross matching iteratively' + msgs.newline() + - '---------------------------------------------------') - continue - if final_fit['rms'] > self.rms_thresh: - msgs.warn('---------------------------------------------------' + msgs.newline() + - 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + - ' Poor RMS ({0:.3f})! Will try cross matching iteratively'.format(final_fit['rms']) + msgs.newline() + - '---------------------------------------------------') - # Store this result in new_bad_slits, so the iteration can be performed, - # but make sure to store the result, as this might be the best possible. - new_bad_slits = np.append(new_bad_slits, slit) - self._all_patt_dict[str(slit)] = copy.deepcopy(patt_dict) - self._all_final_fit[str(slit)] = copy.deepcopy(final_fit) - if self._debug: - xplt = np.linspace(0.0, 1.0, self._npix) - yplt = utils.func_val(final_fit['fitc'], xplt, 'legendre', minx=0.0, maxx=1.0) - plt.plot(final_fit['pixel_fit'], final_fit['wave_fit'], 'bx') - plt.plot(xplt, yplt, 'r-') - plt.show() - #embed() - - # debugging - if self._debug: - # First determine the central wavelength and dispersion of every slit, using the known good solutions - xplt = np.arange(self._nslit) - yplt, dplt = np.zeros(self._nslit), np.zeros(self._nslit) - imsk = np.ones(self._nslit, dtype=int) - for slit in range(self._nslit): - if good_fit[slit]: - yplt[slit] = self._all_patt_dict[str(slit)]['bwv'] - dplt[slit] = self._all_patt_dict[str(slit)]['bdisp'] - imsk[slit] = 0 - - mask, fit = utils.robust_polyfit(xplt, yplt, 2, function='polynomial', sigma=2, - initialmask=imsk, forceimask=True) - - ymodel = utils.func_val(fit, xplt, 'polynomial') - plt.subplot(211) - plt.plot(xplt, ymodel, 'r-') - ww = np.where(mask==0) - plt.plot(xplt[ww], yplt[ww], 'bx') - ww = np.where(mask==1) - plt.plot(xplt[ww], yplt[ww], 'rx') - plt.subplot(212) - plt.plot(xplt, dplt, 'bx') - plt.show() - #embed() - - return new_bad_slits + # This routine is commented out because it is not used. + # def cross_match_order(self, good_fit): + # """Using the solutions of all orders, identify the good solutions, and refit the bad ones! + # + # TODO: This function needs work... The first few lines of code successfully pick up the good orders, + # but we need a new routine that (based on an estimated central wavelength and dispersion) can successfully + # ID all of the lines. + # """ + # # DEPRECATED (NOT USED) + # + # # First determine the central wavelength and dispersion of every slit, using the known good solutions + # xplt = np.arange(self._nslit) + # yplt, dplt = np.zeros(self._nslit), np.zeros(self._nslit) + # imsk = np.ones(self._nslit, dtype=int) + # for slit in range(self._nslit): + # if good_fit[slit]: + # yplt[slit] = self._all_patt_dict[str(slit)]['bwv'] + # dplt[slit] = self._all_patt_dict[str(slit)]['bdisp'] + # imsk[slit] = 0 + # + # mask, fit = utils.robust_polyfit(xplt, yplt, 2, function='polynomial', sigma=2, + # initialmask=imsk, forceimask=True) + # good_fit[mask == 1] = False + # wavemodel = utils.func_val(fit, xplt, 'polynomial') + # disp = np.median(dplt[good_fit]) + # + # # TODO: maybe rethink the model at this point? Using the derived + # # central wavelength and dispersion identify liens in all orders? + # + # if self._debug: + # plt.subplot(211) + # plt.plot(xplt, wavemodel, 'r-') + # ww = np.where(mask==0) + # plt.plot(xplt[ww], yplt[ww], 'bx') + # ww = np.where(mask==1) + # plt.plot(xplt[ww], yplt[ww], 'rx') + # plt.subplot(212) + # plt.plot(xplt, dplt, 'bx') + # plt.show() + # #embed() + # + # fact_nl = 1.2 # Non linear factor + # new_good_fit = np.zeros(self._nslit, dtype=bool) + # for slit in range(self._nslit): + # wmin = wavemodel[slit] - fact_nl*disp*self._npix/2 + # wmax = wavemodel[slit] + fact_nl*disp*self._npix/2 + # ww = np.where((self._wvdata > wmin) & (self._wvdata < wmax)) + # wavedata = self._wvdata[ww] + # msgs.info('Brute force ID for slit {0:d}/{1:d}'.format(slit+1, self._nslit)) + # best_patt_dict, best_final_fit =\ + # self.run_brute_loop(slit, arrerr=self._det_weak[str(slit)], wavedata=wavedata) + # + # self._all_patt_dict[str(slit)] = copy.deepcopy(best_patt_dict) + # self._all_final_fit[str(slit)] = copy.deepcopy(best_final_fit) + # new_good_fit[slit] = self.report_prelim(slit, best_patt_dict, best_final_fit) + # return new_good_fit + # + # + # # Set some fitting parameters + # if self._n_final is None: + # order = 4 + # else: + # order = self._n_final + # + # ofit = [5, 3, 1, 0] + # lnpc = len(ofit) - 1 + # + # # Prepare the fitting coefficients + # xv = np.arange(self._npix)/(self._npix-1) + # ords = np.arange(self._nslit) + # xcen = xv[:, np.newaxis].repeat(self._nslit, axis=1) + # extrapord = ~good_fit + # maskord = np.where(extrapord)[0] + # + # coeffs = None + # waves = np.zeros(xcen.shape, dtype=float) + # for slit in range(self._nslit): + # if good_fit[slit]: + # func = self._all_final_fit[str(slit)]['function'] + # fmin = self._all_final_fit[str(slit)]['fmin'] + # fmax = self._all_final_fit[str(slit)]['fmax'] + # fitc = self._all_final_fit[str(slit)]['fitc'] + # if coeffs is None: + # coeffs = np.zeros((fitc.size, self._nslit)) + # coeffs[:, slit] = fitc.copy() + # waves[:, slit] = utils.func_val(fitc, xv, func, minx=fmin, maxx=fmax) + # + # msgs.info("Performing a PCA on the order wavelength solutions") + # #embed() + # pca_wave, outpar = pca.basis(xcen, waves, coeffs, lnpc, ofit, x0in=ords, mask=maskord, skipx0=False, function=func) + # + # # Report the QA + # # TODO: fix setup passing + # setup = "BLAH" + # pca.pca_plot(setup, outpar, ofit, "wave_cross_match", pcadesc="Wavelength calibration PCA") + # + # + # # Extrapolate the remaining orders requested + # #extrap_wave, outpar = pca.extrapolate(outpar, ords) + # + # # Determine if pixels correlate and anti-correlate with wavelength + # signs = np.zeros(self._nslit, dtype=int) + # for slit in range(self._nslit): + # wvval = pca_wave[:, slit] + # if wvval[wvval.size//2] > wvval[wvval.size//2-1]: + # signs[slit] = 1 + # else: + # signs[slit] = -1 + # sign = 1 + # if np.sum(signs) < 0: + # sign = -1 + # + # new_bad_slits = np.array([], dtype=int) + # # Using the first guesses at the wavelength solution, identify lines + # for slit in range(self._nslit): + # # Get the detections + # dets, _ = self.get_use_tcent(sign, self._det_weak[str(slit)]) + # lindex = np.array([], dtype=int) + # dindex = np.array([], dtype=int) + # # Calculate wavelengths for the gsdet detections + # wvval = pca_wave[:, slit] + # wvcen = wvval[wvval.size//2] + # disp = abs(wvval[wvval.size//2] - wvval[wvval.size//2-1]) + # for dd in range(dets.size): + # pdiff = np.abs(dets[dd] - xv) + # bstpx = np.argmin(pdiff) + # bstwv = np.abs(self._wvdata - wvval[bstpx]) + # if bstwv[np.argmin(bstwv)] > 10.0 * disp: + # # This is probably not a good match + # continue + # lindex = np.append(lindex, np.argmin(bstwv)) + # dindex = np.append(dindex, dd) + # + # # Finalize the best guess of each line + # # Initialise the patterns dictionary + # patt_dict = dict(acceptable=False, nmatch=0, ibest=-1, bwv=0., + # sigdetect=wvutils.parse_param(self._par, 'sigdetect', slit), + # mask=np.zeros(dets.size, dtype=bool)) + # patt_dict['sign'] = sign + # patt_dict['bwv'] = wvcen + # patt_dict['bdisp'] = disp + # + # patterns.solve_triangles(dets, self._wvdata, dindex, lindex, patt_dict) + # # Check if a solution was found + # if not patt_dict['acceptable']: + # new_bad_slits = np.append(new_bad_slits, slit) + # msgs.warn('---------------------------------------------------' + msgs.newline() + + # 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + + # ' Lines could not be identified! Will try cross matching iteratively' + msgs.newline() + + # '---------------------------------------------------') + # continue + # final_fit = self.fit_slit(slit, patt_dict, dets) + # if final_fit is None: + # # This pattern wasn't good enough + # new_bad_slits = np.append(new_bad_slits, slit) + # msgs.warn('---------------------------------------------------' + msgs.newline() + + # 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + + # ' Fit was not good enough! Will try cross matching iteratively' + msgs.newline() + + # '---------------------------------------------------') + # continue + # if final_fit['rms'] > rms_thresh: + # msgs.warn('---------------------------------------------------' + msgs.newline() + + # 'Cross-match report for slit {0:d}/{1:d}:'.format(slit, self._nslit-1) + msgs.newline() + + # ' Poor RMS ({0:.3f})! Will try cross matching iteratively'.format(final_fit['rms']) + msgs.newline() + + # '---------------------------------------------------') + # # Store this result in new_bad_slits, so the iteration can be performed, + # # but make sure to store the result, as this might be the best possible. + # new_bad_slits = np.append(new_bad_slits, slit) + # self._all_patt_dict[str(slit)] = copy.deepcopy(patt_dict) + # self._all_final_fit[str(slit)] = copy.deepcopy(final_fit) + # if self._debug: + # xplt = np.linspace(0.0, 1.0, self._npix) + # yplt = utils.func_val(final_fit['fitc'], xplt, 'legendre', minx=0.0, maxx=1.0) + # plt.plot(final_fit['pixel_fit'], final_fit['wave_fit'], 'bx') + # plt.plot(xplt, yplt, 'r-') + # plt.show() + # #embed() + # + # # debugging + # if self._debug: + # # First determine the central wavelength and dispersion of every slit, using the known good solutions + # xplt = np.arange(self._nslit) + # yplt, dplt = np.zeros(self._nslit), np.zeros(self._nslit) + # imsk = np.ones(self._nslit, dtype=int) + # for slit in range(self._nslit): + # if good_fit[slit]: + # yplt[slit] = self._all_patt_dict[str(slit)]['bwv'] + # dplt[slit] = self._all_patt_dict[str(slit)]['bdisp'] + # imsk[slit] = 0 + # + # mask, fit = utils.robust_polyfit(xplt, yplt, 2, function='polynomial', sigma=2, + # initialmask=imsk, forceimask=True) + # + # ymodel = utils.func_val(fit, xplt, 'polynomial') + # plt.subplot(211) + # plt.plot(xplt, ymodel, 'r-') + # ww = np.where(mask==0) + # plt.plot(xplt[ww], yplt[ww], 'bx') + # ww = np.where(mask==1) + # plt.plot(xplt[ww], yplt[ww], 'rx') + # plt.subplot(212) + # plt.plot(xplt, dplt, 'bx') + # plt.show() + # #embed() + # + # return new_bad_slits def get_use_tcent_old(self, corr, cut=True, arr_err=None, weak=False): """ @@ -2907,6 +2924,12 @@ def solve_slit(self, slit, psols, msols, tcent_ecent, nstore=1, nselw=3, nseld=3 use_tcent, _ = self.get_use_tcent(tpatt_dict['sign'], tcent_ecent) tfinal_dict = wv_fitting.fit_slit(self._spec[:, slit], tpatt_dict, use_tcent, self._line_lists) # tfinal_dict = self.fit_slit(slit, tpatt_dict, use_tcent) + + # get FWHM for this slit + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + # get rms threshold for this slit + rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) + if tfinal_dict is None: # This pattern wasn't good enough continue @@ -2915,7 +2938,7 @@ def solve_slit(self, slit, psols, msols, tcent_ecent, nstore=1, nselw=3, nseld=3 # First time a fit is found patt_dict, final_dict = tpatt_dict, tfinal_dict continue - elif tfinal_dict['rms'] < self.rms_thresh: + elif tfinal_dict['rms'] < rms_thresh: # Has a better fit been identified (i.e. more lines ID)? if len(tfinal_dict['pixel_fit']) > len(final_dict['pixel_fit']): patt_dict, final_dict = copy.deepcopy(tpatt_dict), copy.deepcopy(tfinal_dict) @@ -3036,6 +3059,11 @@ def finalize_fit(self, detections): def report_prelim(self, slit, best_patt_dict, best_final_fit): + # get FWHM for this slit + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + # get rms threshold for this slit + rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) + good_fit = False # Report on the best preliminary result if best_final_fit is None: @@ -3045,7 +3073,7 @@ def report_prelim(self, slit, best_patt_dict, best_final_fit): '---------------------------------------------------') self._all_patt_dict[str(slit)] = None self._all_final_fit[str(slit)] = None - elif best_final_fit['rms'] > self.rms_thresh: + elif best_final_fit['rms'] > rms_thresh: msgs.warn('---------------------------------------------------' + msgs.newline() + 'Preliminary report for slit {0:d}/{1:d}:'.format(slit+1, self._nslit) + msgs.newline() + ' Poor RMS ({0:.3f})! Attempting to cross match.'.format(best_final_fit['rms']) + msgs.newline() + diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py index 920eed62b5..b38e731e50 100644 --- a/pypeit/spectrographs/gemini_gnirs.py +++ b/pypeit/spectrographs/gemini_gnirs.py @@ -322,8 +322,7 @@ def config_specific_par(self, scifile, inp_par=None): par['calibrations']['slitedges']['fit_min_spec_length'] = 0.5 # Wavelengths - par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.4 # Might be grating dependent.. - par['calibrations']['wavelengths']['fwhm'] = 2.5 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.4 par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['lamps'] = ['OH_GNIRS'] # par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] @@ -355,7 +354,7 @@ def config_specific_par(self, scifile, inp_par=None): par['calibrations']['slitedges']['sync_predict'] = 'nearest' # Wavelengths - par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 1.0 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.05 par['calibrations']['wavelengths']['sigdetect'] = 5.0 par['calibrations']['wavelengths']['lamps'] = ['Ar_IR_GNIRS'] # par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] diff --git a/pypeit/spectrographs/magellan_fire.py b/pypeit/spectrographs/magellan_fire.py index 7172708e2e..a4286b4558 100644 --- a/pypeit/spectrographs/magellan_fire.py +++ b/pypeit/spectrographs/magellan_fire.py @@ -146,7 +146,6 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with OH lines par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.2 - par['calibrations']['wavelengths']['fwhm'] = 5. par['calibrations']['wavelengths']['sigdetect']=[5,10,10,10,10,20,30,30,30,30,30,10,30,30,60,30,30,10,20,30,10] par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=[3,3,3,2,4,4,4,3,4,4,4,3,4,4,4,4,4,4,6,6,4] @@ -392,7 +391,7 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with arc lines - par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.1 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.05 par['calibrations']['wavelengths']['sigdetect']=3 par['calibrations']['wavelengths']['fwhm'] = 10 par['calibrations']['wavelengths']['n_first']=2 diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index e760108f3c..64590b4201 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -611,17 +611,12 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Save for redo's self.measured_fwhms = measured_fwhms - # determine rms threshold - _med_fwhm = np.median(measured_fwhms[measured_fwhms!=0]) if np.any(measured_fwhms!=0) else None - fwhm_for_thresh = autoid.set_fwhm(self.par, measured_fwhm=_med_fwhm) - self.wave_rms_thresh = round(self.par['rms_thresh_frac_fwhm'] * fwhm_for_thresh,3) # Obtain calibration for all slits if method == 'holy-grail': # Sometimes works, sometimes fails arcfitter = autoid.HolyGrail(arccen, self.lamps, par=self.par, ok_mask=ok_mask_idx, - rms_thresh=self.wave_rms_thresh, measured_fwhms=self.measured_fwhms, nonlinear_counts=self.nonlinear_counts, spectrograph=self.spectrograph.name) @@ -642,7 +637,7 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Now preferred # Slit positions arcfitter = autoid.ArchiveReid( - arccen, self.lamps, self.par, self.wave_rms_thresh, + arccen, self.lamps, self.par, ech_fixed_format=self.spectrograph.ech_fixed_format, ok_mask=ok_mask_idx, measured_fwhms=self.measured_fwhms, @@ -686,7 +681,7 @@ def build_wv_calib(self, arccen, method, skip_QA=False, #ok_mask_idx = ok_mask_idx[:-1] patt_dict, final_fit = autoid.echelle_wvcalib( arccen, order_vec, arcspec_arxiv, wave_soln_arxiv, - self.lamps, self.par, self.wave_rms_thresh, ok_mask=ok_mask_idx, + self.lamps, self.par, ok_mask=ok_mask_idx, measured_fwhms=self.measured_fwhms, nonlinear_counts=self.nonlinear_counts, debug_all=False, @@ -710,7 +705,11 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # Save the new fits (if they meet tolerance) for key in final_fit.keys(): idx = int(key) - if final_fit[key]['rms'] < self.wave_rms_thresh: + # get FWHM for this slit + fwhm = autoid.set_fwhm(self.par, measured_fwhm=self.measured_fwhms[idx], verbose=True) + # get rms threshold for this slit + wave_rms_thresh = round(par['rms_thresh_frac_fwhm'] * fwhm, 3) + if final_fit[key]['rms'] < wave_rms_thresh: self.wv_calib.wv_fits[idx] = final_fit[key] self.wv_calib.wv_fits[idx].spat_id = self.slits.spat_id[idx] self.wv_calib.wv_fits[idx].fwhm = self.measured_fwhms[idx] @@ -810,6 +809,8 @@ def redo_echelle_orders(self, bad_orders:np.ndarray, dets:np.ndarray, order_dets x2=np.ones_like(spec_vec_norm)*order)/order # get FWHM for this order fwhm = autoid.set_fwhm(self.par, measured_fwhm=self.measured_fwhms[iord], verbose=True) + # get rms threshold for this order + wave_rms_thresh = round(self.par['rms_thresh_frac_fwhm'] * fwhm, 3) # Link me tcent, spec_cont_sub, patt_dict_slit, tot_llist = autoid.match_to_arxiv( self.lamps, self.arccen[:,iord], wv_order_mod, @@ -841,7 +842,7 @@ def redo_echelle_orders(self, bad_orders:np.ndarray, dets:np.ndarray, order_dets # Keep? # TODO -- Make this a parameter? increase_rms = 1.5 - if final_fit['rms'] < increase_rms*self.wave_rms_thresh: + if final_fit['rms'] < increase_rms*wave_rms_thresh: # TODO -- This is repeated from build_wv_calib() # Would be nice to consolidate # QA @@ -1084,8 +1085,15 @@ def run(self, skip_QA=False, debug=False, if self.par['echelle']: # Assess the fits rms = np.array([999. if wvfit.rms is None else wvfit.rms for wvfit in self.wv_calib.wv_fits]) - # Test and scale by measured_fwhms - bad_rms = rms > self.wave_rms_thresh + # Test and scale by measured_fwhms + # get used FWHM + if self.measured_fwhms is not None: + fwhm = [autoid.set_fwhm(self.par, measured_fwhm=m_fwhm) for m_fwhm in self.measured_fwhms] + else: + fwhm = autoid.set_fwhm(self.par, measured_fwhm=None) + # get rms threshold for each slit + wave_rms_thresh = np.round(self.par['rms_thresh_frac_fwhm'] * fwhm, 3) + bad_rms = rms > wave_rms_thresh #embed(header='line 975 of wavecalib.py') if np.any(bad_rms): self.wvc_bpm[bad_rms] = True From 222545b3f77b32fc3b83b33901254ac25af60faa Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Tue, 19 Sep 2023 22:30:22 -0700 Subject: [PATCH 08/12] more comments --- pypeit/core/wavecal/autoid.py | 12 ++++++------ pypeit/wavecalib.py | 10 +++------- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index fea9b87f00..e9f9e40bd8 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -1321,7 +1321,7 @@ def echelle_wvcalib(spec, orders, spec_arxiv, wave_arxiv, lamps, par, fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[iord], verbose=True) # get rms threshold for this slit rms_thresh = round(par['rms_thresh_frac_fwhm'] * fwhm, 3) - msgs.info(f"Using rms_threshold = {rms_thresh} ({par['rms_thresh_frac_fwhm']} * median FWHM)") + msgs.info(f"Using RMS threshold = {rms_thresh} (pixels); RMS/FWHM threshold = {par['rms_thresh_frac_fwhm']}") detections[str(iord)], spec_cont_sub[:, iord], all_patt_dict[str(iord)] = reidentify( spec[:, iord], spec_arxiv[:, iord], wave_arxiv[:, iord], tot_line_list, par['nreid_min'], cc_thresh=cc_thresh, match_toler=par['match_toler'], @@ -1637,7 +1637,7 @@ def __init__(self, spec, lamps, par, ech_fixed_format=False, ok_mask=None, fwhm = set_fwhm(self.par, measured_fwhm=measured_fwhms[slit], verbose=True) # get rms threshold for this slit rms_thresh = round(par['rms_thresh_frac_fwhm'] * fwhm, 3) - msgs.info(f"Using rms_threshold = {rms_thresh} ({self.par['rms_thresh_frac_fwhm']} * median FWHM)") + msgs.info(f"Using RMS threshold = {rms_thresh} (pixels); RMS/FWHM threshold = {self.par['rms_thresh_frac_fwhm']}") self.detections[str(slit)], self.spec_cont_sub[:,slit], self.all_patt_dict[str(slit)] = \ reidentify(self.spec[:,slit], self.spec_arxiv[:,ind_sp], self.wave_soln_arxiv[:,ind_sp], @@ -1884,7 +1884,7 @@ def run_brute_loop(self, slit, tcent_ecent, rms_thresh, wavedata=None): idthresh = 0.5 # Criteria for early return (at least this fraction of lines must have # an ID on either side of the spectrum) - msgs.info(f"Using rms_threshold = {rms_thresh} ({self._par['rms_thresh_frac_fwhm']} * median FWHM)") + msgs.info(f"Using RMS threshold = {rms_thresh} (pixels); RMS/FWHM threshold = {self._par['rms_thresh_frac_fwhm']}") best_patt_dict, best_final_fit = None, None # Loop through parameter space for poly in rng_poly: @@ -2375,7 +2375,7 @@ def cross_match(self, good_fit, detections): continue # get FWHM for this slit - fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[bs]) # get rms threshold for this slit rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) @@ -2926,7 +2926,7 @@ def solve_slit(self, slit, psols, msols, tcent_ecent, nstore=1, nselw=3, nseld=3 # tfinal_dict = self.fit_slit(slit, tpatt_dict, use_tcent) # get FWHM for this slit - fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit]) # get rms threshold for this slit rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) @@ -3060,7 +3060,7 @@ def finalize_fit(self, detections): def report_prelim(self, slit, best_patt_dict, best_final_fit): # get FWHM for this slit - fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit], verbose=True) + fwhm = set_fwhm(self._par, measured_fwhm=self._measured_fwhms[slit]) # get rms threshold for this slit rms_thresh = round(self._par['rms_thresh_frac_fwhm'] * fwhm, 3) diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 64590b4201..90a10c3fc4 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -1085,16 +1085,12 @@ def run(self, skip_QA=False, debug=False, if self.par['echelle']: # Assess the fits rms = np.array([999. if wvfit.rms is None else wvfit.rms for wvfit in self.wv_calib.wv_fits]) - # Test and scale by measured_fwhms # get used FWHM - if self.measured_fwhms is not None: - fwhm = [autoid.set_fwhm(self.par, measured_fwhm=m_fwhm) for m_fwhm in self.measured_fwhms] - else: - fwhm = autoid.set_fwhm(self.par, measured_fwhm=None) - # get rms threshold for each slit + fwhm = self.par['fwhm'] if self.measured_fwhms is None or self.par['fwhm_fromlines'] is False \ + else self.measured_fwhms + # get rms threshold for all orders wave_rms_thresh = np.round(self.par['rms_thresh_frac_fwhm'] * fwhm, 3) bad_rms = rms > wave_rms_thresh - #embed(header='line 975 of wavecalib.py') if np.any(bad_rms): self.wvc_bpm[bad_rms] = True msgs.warn("Masking one or more bad orders (RMS)") From 254c7f8686cb5a6d553a0f5d77f14a36f2e302a5 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Tue, 19 Sep 2023 22:56:49 -0700 Subject: [PATCH 09/12] added docstring --- pypeit/core/wavecal/autoid.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/pypeit/core/wavecal/autoid.py b/pypeit/core/wavecal/autoid.py index e9f9e40bd8..334c3223fe 100644 --- a/pypeit/core/wavecal/autoid.py +++ b/pypeit/core/wavecal/autoid.py @@ -1869,11 +1869,19 @@ def run_brute_loop(self, slit, tcent_ecent, rms_thresh, wavedata=None): Args: slit (int): - tcent_ecent: - rms_thresh: - wavedata: + Slit number + tcent_ecent (list): + List of `numpy.ndarray`_ objects, [tcent, ecent], which are the + centroids and errors of the detections to be used. + rms_thresh (float): + RMS threshold for the wavelength solution fit + wavedata (`numpy.ndarray`_, optional): + Line list; see ``linelist`` argument in, e.g., + :func:`~pypeit.core.wavecal.patterns.triangles`. Returns: + tuple: Returns two dictionaries, one containing information about the best pattern, + and the other containing the information about the best final fit. """ # Set the parameter space that gets searched From 1207c19508bcea4c5d57c87b8feacae492aaa92e Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 20 Sep 2023 10:56:00 -1000 Subject: [PATCH 10/12] fix small bug --- pypeit/wavecalib.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 90a10c3fc4..694c50f6da 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -1087,7 +1087,7 @@ def run(self, skip_QA=False, debug=False, rms = np.array([999. if wvfit.rms is None else wvfit.rms for wvfit in self.wv_calib.wv_fits]) # get used FWHM fwhm = self.par['fwhm'] if self.measured_fwhms is None or self.par['fwhm_fromlines'] is False \ - else self.measured_fwhms + else self.measured_fwhms.astype(float) # get rms threshold for all orders wave_rms_thresh = np.round(self.par['rms_thresh_frac_fwhm'] * fwhm, 3) bad_rms = rms > wave_rms_thresh From 92195df9ee19cb132753ce9b80e772e1febdf288 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Sun, 24 Sep 2023 23:41:51 -1000 Subject: [PATCH 11/12] fix some spectrographs --- pypeit/spectrographs/gemini_gnirs.py | 5 ++++- pypeit/spectrographs/gtc_osiris.py | 3 +++ pypeit/spectrographs/keck_kcwi.py | 1 + pypeit/spectrographs/magellan_fire.py | 3 ++- pypeit/spectrographs/p200_dbsp.py | 1 + pypeit/spectrographs/shane_kast.py | 1 + pypeit/spectrographs/vlt_sinfoni.py | 3 ++- pypeit/tests/test_runpypeit.py | 2 +- pypeit/wavecalib.py | 2 +- 9 files changed, 16 insertions(+), 5 deletions(-) diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py index b38e731e50..675e39d81c 100644 --- a/pypeit/spectrographs/gemini_gnirs.py +++ b/pypeit/spectrographs/gemini_gnirs.py @@ -323,7 +323,7 @@ def config_specific_par(self, scifile, inp_par=None): # Wavelengths par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.4 - par['calibrations']['wavelengths']['sigdetect'] = 5.0 + par['calibrations']['wavelengths']['sigdetect'] = 10.0 par['calibrations']['wavelengths']['lamps'] = ['OH_GNIRS'] # par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] par['calibrations']['wavelengths']['n_first'] = 2 @@ -612,6 +612,9 @@ def default_pypeit_par(cls): par['calibrations']['tilts']['spat_order'] = 1 par['calibrations']['tilts']['spec_order'] = 1 + # wavecalib + par['calibrations']['wavelengths']['fwhm_fromlines'] = False + # Make sure that this is reduced as a slit (as opposed to fiber) spectrograph par['reduce']['cube']['slit_spec'] = True par['reduce']['cube']['grating_corr'] = False diff --git a/pypeit/spectrographs/gtc_osiris.py b/pypeit/spectrographs/gtc_osiris.py index 9c8777746d..95c6e8fb2b 100644 --- a/pypeit/spectrographs/gtc_osiris.py +++ b/pypeit/spectrographs/gtc_osiris.py @@ -460,6 +460,9 @@ def default_pypeit_par(cls): par['calibrations']['tilts']['spat_order'] = 1 par['calibrations']['tilts']['spec_order'] = 1 + # wavecalib + par['calibrations']['wavelengths']['fwhm_fromlines'] = False + # Make sure that this is reduced as a slit (as opposed to fiber) spectrograph par['reduce']['cube']['slit_spec'] = True par['reduce']['cube']['combine'] = False # Make separate spec3d files from the input spec2d files diff --git a/pypeit/spectrographs/keck_kcwi.py b/pypeit/spectrographs/keck_kcwi.py index f917f15011..456ba7a700 100644 --- a/pypeit/spectrographs/keck_kcwi.py +++ b/pypeit/spectrographs/keck_kcwi.py @@ -1221,6 +1221,7 @@ def default_pypeit_par(cls): # KCWI has non-uniform spectral resolution across the field-of-view par['calibrations']['wavelengths']['fwhm_spec_order'] = 1 par['calibrations']['wavelengths']['fwhm_spat_order'] = 2 + par['calibrations']['wavelengths']['fwhm_fromlines'] = False # Alter the method used to combine pixel flats par['calibrations']['pixelflatframe']['process']['combine'] = 'median' diff --git a/pypeit/spectrographs/magellan_fire.py b/pypeit/spectrographs/magellan_fire.py index a4286b4558..f9d1ef1c4e 100644 --- a/pypeit/spectrographs/magellan_fire.py +++ b/pypeit/spectrographs/magellan_fire.py @@ -145,7 +145,8 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with OH lines - par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.2 + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.25 + par['calibrations']['wavelengths']['fwhm_fromlines'] = False par['calibrations']['wavelengths']['sigdetect']=[5,10,10,10,10,20,30,30,30,30,30,10,30,30,60,30,30,10,20,30,10] par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=[3,3,3,2,4,4,4,3,4,4,4,3,4,4,4,4,4,4,6,6,4] diff --git a/pypeit/spectrographs/p200_dbsp.py b/pypeit/spectrographs/p200_dbsp.py index caecea72f0..433398d233 100644 --- a/pypeit/spectrographs/p200_dbsp.py +++ b/pypeit/spectrographs/p200_dbsp.py @@ -385,6 +385,7 @@ def config_specific_par(self, scifile, inp_par=None): par['calibrations']['wavelengths']['reid_arxiv'] = 'p200_dbsp_blue_600_4000_d55.fits' elif grating == '600/4000' and dichroic == 'D68': par['calibrations']['wavelengths']['reid_arxiv'] = 'p200_dbsp_blue_600_4000_d68.fits' + par['calibrations']['wavelengths']['fwhm_fromlines'] = False elif grating == '300/3990' and dichroic == 'D55': par['calibrations']['wavelengths']['reid_arxiv'] = 'p200_dbsp_blue_300_3990_d55.fits' else: diff --git a/pypeit/spectrographs/shane_kast.py b/pypeit/spectrographs/shane_kast.py index e58691757f..18bc22d8b5 100644 --- a/pypeit/spectrographs/shane_kast.py +++ b/pypeit/spectrographs/shane_kast.py @@ -457,6 +457,7 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['NeI','HgI','HeI','ArI'] + par['calibrations']['wavelengths']['fwhm_fromlines'] = False # TODO In case someone wants to use the IR algorithm for shane kast this is the telluric file. Note the IR # algorithm is not the default. diff --git a/pypeit/spectrographs/vlt_sinfoni.py b/pypeit/spectrographs/vlt_sinfoni.py index 42ba6cf71c..49d54d4ec6 100644 --- a/pypeit/spectrographs/vlt_sinfoni.py +++ b/pypeit/spectrographs/vlt_sinfoni.py @@ -76,7 +76,8 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution - par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.23 #0.20 # Might be grating dependent.. + par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.1 + par['calibrations']['wavelengths']['fwhm_fromlines'] = False par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 5.0 par['calibrations']['wavelengths']['n_final']= 4 diff --git a/pypeit/tests/test_runpypeit.py b/pypeit/tests/test_runpypeit.py index eb1c157868..6ac009840c 100644 --- a/pypeit/tests/test_runpypeit.py +++ b/pypeit/tests/test_runpypeit.py @@ -90,7 +90,7 @@ def test_run_pypeit(): specObjs = specobjs.SpecObjs.from_fitsfile(spec1d_file) # Check RMS - assert specObjs[0].WAVE_RMS < 0.02 # difference must be less than 0.02 pixels + assert specObjs[0].WAVE_RMS < 0.08 # RMS must be less than 0.08 pixels # Flexure assert abs(-0.03 - specObjs[0].FLEX_SHIFT_TOTAL) < 0.1 # difference must be less than 0.1 pixels diff --git a/pypeit/wavecalib.py b/pypeit/wavecalib.py index 694c50f6da..102dfa0218 100644 --- a/pypeit/wavecalib.py +++ b/pypeit/wavecalib.py @@ -708,7 +708,7 @@ def build_wv_calib(self, arccen, method, skip_QA=False, # get FWHM for this slit fwhm = autoid.set_fwhm(self.par, measured_fwhm=self.measured_fwhms[idx], verbose=True) # get rms threshold for this slit - wave_rms_thresh = round(par['rms_thresh_frac_fwhm'] * fwhm, 3) + wave_rms_thresh = round(self.par['rms_thresh_frac_fwhm'] * fwhm, 3) if final_fit[key]['rms'] < wave_rms_thresh: self.wv_calib.wv_fits[idx] = final_fit[key] self.wv_calib.wv_fits[idx].spat_id = self.slits.spat_id[idx] From a110d410c84c481313afa6c0527a11a6a000da51 Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Thu, 28 Sep 2023 23:58:21 -1000 Subject: [PATCH 12/12] all use fwhm_fromlines --- pypeit/spectrographs/gemini_gmos.py | 1 - pypeit/spectrographs/gemini_gnirs.py | 3 --- pypeit/spectrographs/gtc_osiris.py | 3 --- pypeit/spectrographs/keck_deimos.py | 1 - pypeit/spectrographs/keck_kcwi.py | 1 - pypeit/spectrographs/keck_lris.py | 1 - pypeit/spectrographs/keck_mosfire.py | 2 -- pypeit/spectrographs/keck_nires.py | 1 - pypeit/spectrographs/ldt_deveny.py | 1 - pypeit/spectrographs/magellan_fire.py | 1 - pypeit/spectrographs/p200_dbsp.py | 1 - pypeit/spectrographs/shane_kast.py | 1 - pypeit/spectrographs/vlt_sinfoni.py | 1 - pypeit/spectrographs/vlt_xshooter.py | 1 - 14 files changed, 19 deletions(-) diff --git a/pypeit/spectrographs/gemini_gmos.py b/pypeit/spectrographs/gemini_gmos.py index 3012af71d7..b16d279c50 100644 --- a/pypeit/spectrographs/gemini_gmos.py +++ b/pypeit/spectrographs/gemini_gmos.py @@ -326,7 +326,6 @@ def config_specific_par(self, scifile, inp_par=None): # Allow for various binning binning = parse.parse_binning(self.get_meta_value(headarr, 'binning')) - par['calibrations']['wavelengths']['fwhm_fromlines'] = True return par diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py index 675e39d81c..73ad5fdea8 100644 --- a/pypeit/spectrographs/gemini_gnirs.py +++ b/pypeit/spectrographs/gemini_gnirs.py @@ -612,9 +612,6 @@ def default_pypeit_par(cls): par['calibrations']['tilts']['spat_order'] = 1 par['calibrations']['tilts']['spec_order'] = 1 - # wavecalib - par['calibrations']['wavelengths']['fwhm_fromlines'] = False - # Make sure that this is reduced as a slit (as opposed to fiber) spectrograph par['reduce']['cube']['slit_spec'] = True par['reduce']['cube']['grating_corr'] = False diff --git a/pypeit/spectrographs/gtc_osiris.py b/pypeit/spectrographs/gtc_osiris.py index 95c6e8fb2b..9c8777746d 100644 --- a/pypeit/spectrographs/gtc_osiris.py +++ b/pypeit/spectrographs/gtc_osiris.py @@ -460,9 +460,6 @@ def default_pypeit_par(cls): par['calibrations']['tilts']['spat_order'] = 1 par['calibrations']['tilts']['spec_order'] = 1 - # wavecalib - par['calibrations']['wavelengths']['fwhm_fromlines'] = False - # Make sure that this is reduced as a slit (as opposed to fiber) spectrograph par['reduce']['cube']['slit_spec'] = True par['reduce']['cube']['combine'] = False # Make separate spec3d files from the input spec2d files diff --git a/pypeit/spectrographs/keck_deimos.py b/pypeit/spectrographs/keck_deimos.py index 7dc87a1054..6e7044ef90 100644 --- a/pypeit/spectrographs/keck_deimos.py +++ b/pypeit/spectrographs/keck_deimos.py @@ -395,7 +395,6 @@ def config_specific_par(self, scifile, inp_par=None): # Wavelength FWHM binning = parse.parse_binning(self.get_meta_value(headarr, 'binning')) par['calibrations']['wavelengths']['fwhm'] = 6.0 / binning[1] - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # Objects FWHM # Find objects diff --git a/pypeit/spectrographs/keck_kcwi.py b/pypeit/spectrographs/keck_kcwi.py index 456ba7a700..f917f15011 100644 --- a/pypeit/spectrographs/keck_kcwi.py +++ b/pypeit/spectrographs/keck_kcwi.py @@ -1221,7 +1221,6 @@ def default_pypeit_par(cls): # KCWI has non-uniform spectral resolution across the field-of-view par['calibrations']['wavelengths']['fwhm_spec_order'] = 1 par['calibrations']['wavelengths']['fwhm_spat_order'] = 2 - par['calibrations']['wavelengths']['fwhm_fromlines'] = False # Alter the method used to combine pixel flats par['calibrations']['pixelflatframe']['process']['combine'] = 'median' diff --git a/pypeit/spectrographs/keck_lris.py b/pypeit/spectrographs/keck_lris.py index 96c51874a3..b2e7cb5421 100644 --- a/pypeit/spectrographs/keck_lris.py +++ b/pypeit/spectrographs/keck_lris.py @@ -136,7 +136,6 @@ def config_specific_par(self, scifile, inp_par=None): # Wave FWHM binning = parse.parse_binning(self.get_meta_value(scifile, 'binning')) par['calibrations']['wavelengths']['fwhm'] = 8.0 / binning[0] - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # Arc lamps list from header par['calibrations']['wavelengths']['lamps'] = ['use_header'] diff --git a/pypeit/spectrographs/keck_mosfire.py b/pypeit/spectrographs/keck_mosfire.py index 177a74134e..9059901347 100644 --- a/pypeit/spectrographs/keck_mosfire.py +++ b/pypeit/spectrographs/keck_mosfire.py @@ -210,7 +210,6 @@ def config_specific_par(self, scifile, inp_par=None): # using OH lines if 'long2pos_specphot' not in decker and filter in supported_filters: par['calibrations']['wavelengths']['method'] = 'full_template' - par['calibrations']['wavelengths']['fwhm_fromlines'] = True par['calibrations']['wavelengths']['sigdetect'] = 10. # templates if filter == 'Y': @@ -233,7 +232,6 @@ def config_specific_par(self, scifile, inp_par=None): elif 'long2pos_specphot' in decker and filter in supported_filters: par['calibrations']['wavelengths']['lamps'] = ['Ar_IR_MOSFIRE', 'Ne_IR_MOSFIRE'] par['calibrations']['wavelengths']['method'] = 'full_template' - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # templates if filter == 'Y': par['calibrations']['wavelengths']['reid_arxiv'] = 'keck_mosfire_arcs_Y.fits' diff --git a/pypeit/spectrographs/keck_nires.py b/pypeit/spectrographs/keck_nires.py index af527ccfb1..bd4b3b1138 100644 --- a/pypeit/spectrographs/keck_nires.py +++ b/pypeit/spectrographs/keck_nires.py @@ -83,7 +83,6 @@ def default_pypeit_par(cls): par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.136 par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 2.2 # Measured - par['calibrations']['wavelengths']['fwhm_fromlines'] = True par['calibrations']['wavelengths']['n_final']= [3,4,4,4,4] par['calibrations']['wavelengths']['lamps'] = ['OH_NIRES'] #par['calibrations']['wavelengths']['nonlinear_counts'] = self.detector[0]['nonlinear'] * self.detector[0]['saturation'] diff --git a/pypeit/spectrographs/ldt_deveny.py b/pypeit/spectrographs/ldt_deveny.py index 2d5ce4d884..d18c4f6971 100644 --- a/pypeit/spectrographs/ldt_deveny.py +++ b/pypeit/spectrographs/ldt_deveny.py @@ -289,7 +289,6 @@ def default_pypeit_par(cls): # Set this as default... but use `holy-grail` for DV4, DV8 par['calibrations']['wavelengths']['method'] = 'full_template' # The DeVeny arc line FWHM varies based on slitwidth used - par['calibrations']['wavelengths']['fwhm_fromlines'] = True par['calibrations']['wavelengths']['nsnippet'] = 1 # Default: 2 # Because of the wide wavelength range, solution more non-linear; user higher orders par['calibrations']['wavelengths']['n_first'] = 3 # Default: 2 diff --git a/pypeit/spectrographs/magellan_fire.py b/pypeit/spectrographs/magellan_fire.py index f9d1ef1c4e..a081a8971e 100644 --- a/pypeit/spectrographs/magellan_fire.py +++ b/pypeit/spectrographs/magellan_fire.py @@ -146,7 +146,6 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution with OH lines par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.25 - par['calibrations']['wavelengths']['fwhm_fromlines'] = False par['calibrations']['wavelengths']['sigdetect']=[5,10,10,10,10,20,30,30,30,30,30,10,30,30,60,30,30,10,20,30,10] par['calibrations']['wavelengths']['n_first']=2 par['calibrations']['wavelengths']['n_final']=[3,3,3,2,4,4,4,3,4,4,4,3,4,4,4,4,4,4,6,6,4] diff --git a/pypeit/spectrographs/p200_dbsp.py b/pypeit/spectrographs/p200_dbsp.py index 433398d233..caecea72f0 100644 --- a/pypeit/spectrographs/p200_dbsp.py +++ b/pypeit/spectrographs/p200_dbsp.py @@ -385,7 +385,6 @@ def config_specific_par(self, scifile, inp_par=None): par['calibrations']['wavelengths']['reid_arxiv'] = 'p200_dbsp_blue_600_4000_d55.fits' elif grating == '600/4000' and dichroic == 'D68': par['calibrations']['wavelengths']['reid_arxiv'] = 'p200_dbsp_blue_600_4000_d68.fits' - par['calibrations']['wavelengths']['fwhm_fromlines'] = False elif grating == '300/3990' and dichroic == 'D55': par['calibrations']['wavelengths']['reid_arxiv'] = 'p200_dbsp_blue_300_3990_d55.fits' else: diff --git a/pypeit/spectrographs/shane_kast.py b/pypeit/spectrographs/shane_kast.py index 18bc22d8b5..e58691757f 100644 --- a/pypeit/spectrographs/shane_kast.py +++ b/pypeit/spectrographs/shane_kast.py @@ -457,7 +457,6 @@ def default_pypeit_par(cls): # 1D wavelength solution par['calibrations']['wavelengths']['lamps'] = ['NeI','HgI','HeI','ArI'] - par['calibrations']['wavelengths']['fwhm_fromlines'] = False # TODO In case someone wants to use the IR algorithm for shane kast this is the telluric file. Note the IR # algorithm is not the default. diff --git a/pypeit/spectrographs/vlt_sinfoni.py b/pypeit/spectrographs/vlt_sinfoni.py index 49d54d4ec6..68f40af6b3 100644 --- a/pypeit/spectrographs/vlt_sinfoni.py +++ b/pypeit/spectrographs/vlt_sinfoni.py @@ -77,7 +77,6 @@ def default_pypeit_par(cls): # Wavelengths # 1D wavelength solution par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.1 - par['calibrations']['wavelengths']['fwhm_fromlines'] = False par['calibrations']['wavelengths']['sigdetect']=5.0 par['calibrations']['wavelengths']['fwhm']= 5.0 par['calibrations']['wavelengths']['n_final']= 4 diff --git a/pypeit/spectrographs/vlt_xshooter.py b/pypeit/spectrographs/vlt_xshooter.py index 7ec5db8c21..b2ffa58ceb 100644 --- a/pypeit/spectrographs/vlt_xshooter.py +++ b/pypeit/spectrographs/vlt_xshooter.py @@ -956,7 +956,6 @@ def default_pypeit_par(cls): # This is for 1x1 par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.184 par['calibrations']['wavelengths']['fwhm'] = 3.8 - par['calibrations']['wavelengths']['fwhm_fromlines'] = True # par['calibrations']['wavelengths']['sigdetect'] = 3.0 # Pretty faint lines in places # Reidentification parameters