Question related to final tabulation of lines/calibrating the FIES spectrograph

[jsinkbaek]

Dear Mirek,

I noticed that the number of tabulated lines is higher than the number of lines used for the refined wavelength solution. I assume this is due to sigma-clipping. Is there information stored so that I can select only those retained for the final solution?

Relevant lines in my feed: 2022-09-01 20:14:56,257 - IdentifyArcEmissionLines - do_stage_fiber - INFO - 3529 emission lines identified from 91 unique diffraction orders

2022-09-01 20:15:08,054 - SolutionRefineFinal - do_stage_fiber - INFO - 2372 lines within 4.5 median absolute deviations and 3523 lines within 4.5 standard deviations

Oddly enough, the amount stored is not exactly equal to 3529, but instead 3523. Is there another reason for this?

All the best, Jeppe

[gmbrandt]

Dear Jeppe,

It is great to see you using the code, and I am glad to see that it (appears) to be working. For my own curiosity, does it appear to be giving a sensible solution for your instrument?

For your first question, it would help me to know what the output that you currently have is. Do you have the wavelengths (wavelengths output by the xwavecal model) of those 3523 stored lines?

Related, what is your 'final_mad_clip' set to in your config file?

And I am unsure as to the answer to the second question about why 6 detected lines are missing. I just surveyed the code, and I did not see any place where emission lines are forever removed from the list. We could zoom about this issue (email me at gmbrandt@ucsb.edu), but I suspect that may take some time to track down. Or it could be some small feature (like lines with 'nan' or negative flux are removed), that I forgot is in the code. I suggest we resolve the above issue first, then save the missing 6 lines issue for later.

Best,

Mirek

[gmbrandt]

Just as a quick followup, I just now re-ran the e2e test case for NRES. At least for that case, the number of lines detected (and printed out to console) matches the number of lines saved to the final table that is written to disc. So I think focusing on the first issue is probably best.

[jsinkbaek]

Hi Mirek,

Thank you for the quick reply!

At the moment, I don't think I'm getting something fully sensible out. However it's difficult to verify, since at a surface glance I am getting wavelength vs. reference wavelength of the identified lines which seems semi-sensible.

I'm using data from the FIES spectrograph. In its high-res fiber setting, it has a resolution of 67000 and a wavelength coverage of approximately 3600-8700Å spread over 91 orders.

Since the red orders are very over-exposed, I'm using only 81 orders (orders 10-90). I have kept the original ref-ids and ids such that the principal order number still makes sense. With this, xwavecal finds lines from 3200-6600Å, which doesn't seem right.

I attach a sample file and config in the link: https://drive.google.com/drive/folders/1TMDO0yy_-5Czi9ji2cyRI5ah23UKid0E?usp=sharing

The stdout from the run I include here: log.txt

[jsinkbaek]

For clarity:

The input spectrum is FIFg300264_interm.fits

The output spectrum is LaPalma_NOT_20220731_FIFg300264.fits_wavecal_100.fits

All the best, Jeppe

[gmbrandt]

Thanks! I requested access.

[jsinkbaek]

I think I allowed you access immediately upon your request no? Anyway, I just opened it for all now, since that was my original intention.

[gmbrandt]

Yep, I downloaded the files. Thanks. I'll take a look.

[gmbrandt]

OK I did some significant digging. At first, the overlap fitting was failing. There were many false matches.

The only major problem I saw, was actually with the data itself:

the estimate of the standard error appeared way off (many factors of 10, as far as I can tell). Accurate errors (i.e., correct to within a factor of 2-3) are crucial to the performance of the algorithm (in particular, the overlap fitting stage). So I overrode the error estimate you included with sqrt(counts) + a light read noise. The most important part is getting the Poisson noise right, because High S/N lines matter the most for us.

Minor changes that went along with that change:

min_num_matches = 10 # so that for an overlap to be verified, 10 or more peaks need to line up. I looked at the spectrum. We have so many beautiful matches in this spectrum, that we might as well be picky.

increased flux tol to flux_tol = 0.4

overlap_min_peak_snr = 30 # from visual inspection of the spectra, after fixing the standard error estimates. Basically, the "noise floor" on most orders was around a S/N 10. This could be due to a constant offset above 0 that isn't corrected for. I don't really want to mess with trying to fix that. I just tweaked this number a little bit until I got about 5000 lines in the "xwavecal.wavelength.IdentifyArcEmissionLinesLowSN" stage.

After these changes, the overlap fitting is working very well.

Here is an example of two orders, where one has been stretched according to the mapping xwavecal found.

[gmbrandt]

But still, the wavelength solution was failing. Which is very bizarre. For a standard ThAr spectrum, fitting the overlaps should be the hardest part. I checked, and the solution is failing on the FindGlobalScale stage.

This should only happen in cases where the line list is a very bad match to the actual lines present on the detector.

So @jsinkbaek , is there any reason to suspect that the observed lines on the detector are substantially different from that normally present in ThAr? Is the ThAr lamp very old? Are there other trace elements in there?

[gmbrandt]

P.S. here is the updated config file:

[data]
# The data class to use for the data. The class must implement all of the methods of xwavecal.images.DataProduct
data_class = xwavecal.images.Image
# the extension number or name in the .fits file where the raw data for the exposure lies. Typically 0
# for non-compressed files, and 1 for .fits.fz files.
primary_data_extension = 0
# elements that must be contained in the filenames of the files to reduce. E.g. ['.fits', '20190405']
files_contain = ['.fits', 'interm_correcterrs']
# the fits header keys which gives various information about the data frame.
# instrument2 is just another specifier for the instrument. If not needed, set it to the same as 'instrument'
header_keys = {'type': 'OBJECT',
               'gain': 'GAIN',
               'read_noise': 'RDNOISE',
               'data_section': 'TRIMSEC',
               'overscan_section': 'BIASSEC',
               'fiber_state': 'FIFMSKNM',
               'observation_date': 'DATE-OBS',
               'instrument': 'TELESCOP',
               'instrument2': 'INSTRUME',
               'site_name': 'OBSERVAT',
              # unique_id should give some kind of running tag for each raw frame, per reduction batch.
              # E.g. 1 or 14 or 52 etc..
               'unique_id': 'FILENAME'}
# the corresponding possible entries of 'type' and which ones correspond to wavecal and lampflat.
# format as a python dict.
type_keys = {'LAMPFLAT': 'lampflat',
             'ThAr': 'wavecal',
             'F4 HiRes': 'thar&none&none'
            }

[stages]
# Reduction stages for a wavelength calibration frame, in order. Comment out MakeFiberTemplate on successive runs.
wavecal = [
          'xwavecal.fibers.MakeFiberTemplate',
          'xwavecal.fibers.IdentifyFibers',
          'xwavecal.wavelength.Initialize',
          'xwavecal.wavelength.LoadReferenceLineList',
          'xwavecal.wavelength.IdentifyArcEmissionLinesLowSN',
          'xwavecal.wavelength.BlazeCorrectArcEmissionLines',
          'xwavecal.wavelength.FitOverlaps',
          'xwavecal.wavelength.IdentifyArcEmissionLines',
          #'xwavecal.wavelength.IdentifyPrincipleOrderNumber',
          'xwavecal.wavelength.SolveFromOverlaps',
          'xwavecal.wavelength.FindGlobalScale',
          'xwavecal.wavelength.SolutionRefineInitial',
          'xwavecal.wavelength.SolutionRefineFinal',
          'xwavecal.wavelength.IdentifyArcEmissionLinesLowSN',
          'xwavecal.wavelength.ApplyToSpectrum',
          'xwavecal.wavelength.TabulateArcEmissionLines']

[reduction]
# path to the reference line list
line_list_path = "/home/gmbrandt/Documents/Repositories/xwavecal/xwavecal/data/ThAr_atlas_ESO.txt"
# database which records processed images, calibration images etc.
database_path = "/home/gmbrandt/Downloads/FIES_spectrograph_jsinkbaekspectra/pipeline.db"
# time format for 'observation_date' in the fits headers. Must comply with datetime.datetime.strftime etc..
# For NRES below, a timestamp in 'DATE-OBS' looks like '2019-04-10T12:56:44.466'. Use %% in place of %.
time_format = "%%Y-%%m-%%dT%%H:%%M:%%S.%%f"

# settings for each wavelength reduction stage

# the name of the blaze corrected spectrum.
blaze_corrected_spectrum_name = "BLZCORR"

# Extraction settings #
box_spectrum_name = "SPECBOX"
ivar_spectrum_name = "SPECIVR"

# order identification settings #
# reference id of the central diffraction order in the database arc template.
ref_id = 45
# trace id of the diffraction order to use to create an arc template (only necessary if using MakeFiberTemplate)
template_trace_id = 45
# Set low_fiber_first = False if you find that the fiber designation in extracted spectra is the exact
# opposite of what it should be. low_fiber_first = True means the lowest lying fiber (lowest y coordinate)
# will be identified with the first lit fiber. See xwavecal.fibers.IdentifyFibers.build_fiber_column
# for a more thorough explanation with examples.
low_fiber_first = False

# Wavelength calibration settings #
# spectrum to wavelength calibrate. E.g. SPECBOX.
main_spectrum_name = "SPECBOX"

overlap_table_name = "OVERLAP"
# the min signal to noise for an emission peak to be used to fit the overlaps.
overlap_min_peak_snr = 30
# float between 0 and 1. Two peaks must have blaze corrected fluxes which agree within
# this tolerance to be counted as a match.
flux_tol = 0.4

emission_lines_table_name = "LINES"
# default below: 1000
max_red_overlap = 1000
# default below: 2000
max_blue_overlap = 1200
overlap_linear_scale_range = (0.5, 2)
min_num_matches = 10

# Default below: (0.8, 1.2)
global_scale_range = (0.8, 1.2)
# Default for R=50000: 10. Default for R=500000: 1.
global_scale_spacing = 5
approx_detector_range_angstroms = 5000
approx_num_orders = 81
# principle order number (m0) settings for the IdentifyPrincipleOrderNumber stage
# start (inclusive), stop (exclusive)
m0_range = (60, 70)
principle_order_number = 64
# the minimum number of overlaps which must be fit for the wavelength solution to proceed.
min_num_overlaps = 5
# the min signal to noise for an emission peak to be considered in the wavelength solution. Default: 20
min_peak_snr = 30

# initial model that is constrained via the overlaps and is used to find the global scale:
# format is: {xpower: [ipower, ipower,...],..}
# e.g. {1: [0, 1, 2]} represents (1/(m0+i)) * (a*x * (b + c*i + d*i^2)). The 1/(m0+i) prefactor is always included.
initial_wavelength_model = {1: [0, 1, 2],
                            2: [0, 1, 2]}
# wavelength model for the initial refine stage:
intermediate_wavelength_model = {0: [0, 1, 2],
                                 1: [0, 1, 2],
                                 2: [0, 1, 2]}
# wavelength model that the final refine stage will end at (stage begins with the intermediate model):
final_wavelength_model = {0: [0, 1, 2, 3, 4, 5],
                          1: [0, 1, 2, 3, 4, 5],
                          2: [0, 1, 2, 3, 4, 5],
                          3: [0, 1, 2, 3, 4, 5],
                          4: [0]}
# outliers will be ignored when solving if their residuals are k m.a.d.s away from 0. k set by initial_mad_clip in SolutionRefineInitial
initial_mad_clip = 6
# k set by initial_mad_clip in SolutionRefineFinal
final_mad_clip = 4

[jsinkbaek]

Dear Mirek,

Thank you for the very detailed report. I will look it over and attempt again sometime next week.

The error makes a lot of sense! I lacked an error estimate, so I manually made one by measuring the standard deviation of the nearest 9 points and dividing with sqrt(9-1). The last part might be the problem, as that means errors decrease if I used more than 9 points, and increase if I use less, which is not realistic in this case.

I originally estimated the flux error from the blaze function. Problem there is that it only measures the overall efficiency, and not the noise of an individual line.

Noise-floor also is a reasonable culprit. I have been told by the instrument astronomer that the extraction routine does not background subtract ThAr spectra.

I know that the lamp has very recently been changed after these spectra were taken as it was failing rapidly after. However, these same spectra have been wavelength calibrated with in-house software without issue, but using a curated list of around 600-750 lines. It might be that they are less sensitive to lamp age. It is not related to over-exposed lines right? Because there are plenty of those in the red.

I did orient the spectrum correctly, right?

[gmbrandt]

"It is not related to over-exposed lines right? Because there are plenty of those in the red." This very well could be. This is a good idea. I'm going to try the following today: rerun the reduction with just order id's 30-91, instead of 10-91.

Yep you did orient the spectrum perfectly!

"I know that the lamp has very recently been changed after these spectra were taken as it was failing rapidly after." Oh interesting. This is unrelated, but I would be highly interested in trying to calibrate a new spectrum if you have one on hand.

[gmbrandt]

Oh WOW, ok truncating the spectrum to orders 35 and onward worked immediately. You get what appears to be an excellent wavelength solution, with zero effort. Sweet.

Heres a gdrive link to the data and input files etc. https://drive.google.com/drive/folders/18LZMf48ULXv9iVedJeSij4zyVouhiyRy?usp=sharing

I included the scripts to "fix" the errors. I ran first "correct_error_estimates.py" then "create_truncated_spectrum.py" then I reduced xwavecal using the config config_220731_truncated_spectrum.ini The output files are in output/

Could you inspect that output wavelength calibrated spectrum and let me know if it is indeed great quality? From my inspection, it looks excellent across the entire range of orders 35-90 (wavelength range 3670 to 5810 angstroms. The average velocity dispersion I get is roughly 11 m/s (converting delta lambda/lambda to delta v, then dividing by the sqrt of the number of inlier lines; there are about 1500 inlier lines).

Cheers,

Mirek

[jsinkbaek]

HI Mirek,

Thank you for the hard work! That sounds correct, see:

http://www.not.iac.es/instruments/fies/orderlist

Since we are starting on order 64 (I checked the IRAF solution for the original wavelength calibration), 35 orders up is order 99, which is 5731-5801 Å. Since it is an old order-list and the ccd has been moved since, that range could very well be a couple of 10s of ångstrom off.

I will try and inspect your output and see!

[jsinkbaek]

Hi Mirek,

I see in your truncate_spectrum.py you have from order 25: onward?

[gmbrandt]

(edited) Yeah but the spectrum you sent me only had orders 10 onward (in terms of relative id). I clipped the first 25, so now we just have 35 and onward. (edit): the reference ID is not updated according to the principle order number (64).

[jsinkbaek]

Ahh ok that makes sense!

I looked at my spectra, and it looks like linearity already starts to happen around the 20th order. I have plotted order 15-22:

[gmbrandt]

Ah ok, great. I would go ahead and give it a try then with orders 20 onward.

If it doesn't work immediately, here are a couple options:

1. restricting your RV analysis to <5800 Angstroms I think will be ok (so using 35 onward). The red will be littered with tellurics anyway. Although it would be a shame if you missed H alpha.
2. If 1., is not an option, and you need 1-35, then mask all the saturated lines in the spectrum in the red. You'll want to be sure that the masked values = the background, so that there are no divots missing in the spectrum.

P.S. just for posterity, I ran the Identify Order Number stage on this spectrum. There are actually two principle order numbers which work, 63 and 64. I don't think there is any significance to this, just that the model is sufficiently flexible in order number (apparently) to make up for the off-by-one error induced by using 63. But I think it's valuable to point this out. This is the first time I have come across a spectrum where 2 principle order numbers both work equally as well. I would probably just stick with m0=64, if that is what the instrument scientists have said.

P.P.S. if you want to make a jupyter notebook tutorial of calibrating FIES using xwavecal, that would be a welcome contribution. Just go ahead and make a branch and let me know when you are ready for a pull request.

[jsinkbaek]

Hi Mirek,

Thank you so much for your help. It has been an amazing experience to get this much feedback and assistance.

Yes of course, for stellar work it should be okay. However, my specific interest is in the individual ThAr lines, so the more lines I can identify the more data I have.

I am surprised that it works for 2 different starting orders. Maybe it's because the orders are very compact on the CCD? It is only around 2k X 2k pixels.

I have made some standardized ThAr .fits file preparation tools, which should be relatively easy to adapt to a jupyter notebook. I will append it to my list of stuff to do.

I suspect that we should be able to get better internal precision than 11 m/s with this method. Using a global RV fitting method to an empirical template (co-added ThArs), I was able to get an average 1.8 m/s internal precision for each spectrum. The error was attained by applying the method to each order independently. So either some of the lines could be misidentified, or there is a wealth of information hidden in hard-to-identify weak lines below the S/N limit.

What do you think?

All the best, Jeppe

[jsinkbaek]

Or perhaps the lines just need to be properly weighted to reduce the error?

[gmbrandt]

It is my pleasure! It is very nice and rewarding that you want to use xwavecal, and I am really excited that it is working.

I am surprised that it works for 2 different starting orders. Maybe it's because the orders are very compact on the CCD? It is only around 2k X 2k pixels.

I am also, but I think there is little significance here to that. Certainly no physical significance. I believe the reason "it works" is purely computational. The xwavecal model is a large Legendre polynomial series in both order i and pixel x, so I think it is saying that that model has enough flexibility to overcome the error of using m0=63. My guess is that if you are able to add in orders 1-35, that this will NOT be the case anymore, and that m0=64 will be the only winner. This is because i/(m0+1) and i/(m0) are not very different for high order numbers, but for low order numbers they are. Because we have so far excluded low order numbers, we were able to get away with two m0's that worked.

I suspect that we should be able to get better internal precision than 11 m/s with this method. Using a global RV fitting method to an empirical template (co-added ThArs), I was able to get an average 1.8 m/s internal precision for each spectrum. The error was attained by applying the method to each order independently. So either some of the lines could be misidentified, or there is a wealth of information hidden in hard-to-identify weak lines below the S/N limit.

I think I have an answer to this, but I might be misunderstanding. How exactly is "internal precision" defined here?

[jsinkbaek]

I perform a least-squares fit of an empirical template to each order separately, by measuring a chi^2 for several trial RVs and then fitting the minimum with a polynomial. Each order RV then has a relative error defined by the shape of the polynomial.

I sigma-clip the order RVs with 4 or 5 sigma.

My final global RV measurement is then the weighted average of each order, and the error is a weighted standard deviation / sqrt(n_orders -1)

[gmbrandt]

Ahhh, OK. I am following. Yes so the 11 m/s value I reported to you is the absolute RV precision. I probably should stop quoting the absolute RV precision because it is not that useful in the end. What you have measured is the more relevant quantity, which is the relative RV precision.

You could repeat the same measurement on the xwavecal calibrated spectrum (and after calibrating the stacked, empirical ThAr template). I suspect you will approach this 2 m/s number.

Best,

Mirek

[jsinkbaek]

Ahh so you mean what average RV offset we are measuring from the line list?

Yes so my plan is to use 1 wavelength calibration (to measure RVs relative to that), and then apply that to all the measured pixel values of the lines of all the other spectra.

All the best, Jeppe

[gmbrandt]

Ahh so you mean what average RV offset we are measuring from the line list?

Sort of. The standard deviation of this, as opposed to the offset. The offset is roughly 0 (as it should be). When I quoted the absolute precision, I calculated: (wavelength of line j, under the xwavecal best-fit model) - (wavelength of line j, from the reference line list) For all the lines, j, that were used in the solution (so 2000 or whatever). Then I calculate the standard deviation of that list of residuals, and call that delta_lambda. Then I do delta_v = delta_lamba/lambda * c and then divide by sqrt(N).

Yes so my plan is to use 1 wavelength calibration (to measure RVs relative to that), and then apply that to all the measured pixel values of the lines of all the other spectra.

Sounds reasonable. I have done plenty of RV work where I do a single xwavecal solution and apply it to all the spectra in the time series. This works pretty well when all the spectra are from ~1 month (over which NRES is pretty stable).

[jsinkbaek]

Would it not be simpler to convert all of the delta_wavelengths to RV and then do the standard deviation?

I get 23 m/s, but I'm also using 4400 lines across order 20-90. How do you select only the best-fit ones? Or are the best-fit lines already the ones that are output to hdul[4] in the fits file?

[gmbrandt]

Would it not be simpler to convert all of the delta_wavelengths to RV and then do the standard deviation?

That should be fine too. I think this ends up being equivalent, as std(RV) = c* std(delta_lambda/lambda) -- but I am not sure. You could try both ways.

How do you select only the best-fit ones?

To select the best fit lines, you have to do this manually. But is not hard. You would do:

hdul = fits.open('output/LaPalma_NOT_20220731_FIFg300264.fits_wavecal_100.fits')
lines = Table(hdul['lines'].data)
resids = lines['wavelength'] - lines['reference_wavelength']
plt.figure()
plt.hist(resids, bins=np.linspace(-0.03, 0.03, 100))
plt.show()

which would yield this plot:

The distribution is roughly concentrated between -0.015 Ang and +0.015 Ang. So you could restrict your analysis to those lines, by doing:

where_best_residuals = np.abs(resids) < 0.015
best_fit_lines = lines[where_best_residuals]

[jsinkbaek]

Essentially a manual sigma rejection, sounds simple enough. Thank you!

[gmbrandt]

Yep, sorry I have a habit of giving long winded answers to every question.

[jsinkbaek]

No problem. By combining 170 spectra into one mean master (and setting the S/N limits to around 200), I manage to get an absolute rv error of 5 m/s by following your prescription. However, in this I see an absolute mean drift of ~-5 m/s for the ESO line-list, and ~50 m/s when using a line-list from Redman et al 2014.

[jsinkbaek]

I was also wondering if it was possible to use the flux levels for the line identification?

Right now my solution is to manually discard lines that are either below 0.3 or above 2 times the expected line intensity. But I am wondering if that information could be used as a criterium along with the pixel values?

[gmbrandt]

I was also wondering if it was possible to use the flux levels for the line identification? Right now my solution is to manually discard lines that are either below 0.3 or above 2 times the expected line intensity. But I am wondering if that information could be used as a criterium along with the pixel values?

Ah, are you referring to tailoring the line list to the instrument? You could do such a thing with line intensities. irrc, the estimated intensities are usually pretty good (like within 50%, sort of like the range you suggest), but they are almost never spot on. But yeah, I mean it's probably worth a shot. Such a capability is not built in to xwavecal though.

You should also checkout the line list from Murphy et al. 2007. This is an improvement over the original ESO ThAr atlas. We don't use it for NRES, but we should. I included in that google drive one version that is ready to go (Murphy2007_ll.txt ; note these wavelengths have been converted to vacuum). The calibration fails with this line list, I think because of the larger density of lines, but maybe one could get it to work if you curated it as you suggest.

[jsinkbaek]

As mentioned, I had success with using a (VERY) high S/N requirement while using a combined master frame and an extensive line list by Redman 2014. Have you used that list before?

Then I rejected any line detections where the median normalized flux within each order was less than 0.3 or larger than 2 times the median normalized line intensity. After, I sigma-rejected by 1 sigma, which corresponded to a wavelength deviation of 0.17Å.

Then, I went through the remaining lines, and each time I had a reference line detected 2 or more times in 1 order, I kept only the value with the lowest pixel error (is this best practice?). If I then had a reference line detected multiple times across multiple orders after this, I took the average weighted by 1/pix_err**2.

From this, I obtained around 2000 final lines in a reference line list with instrument specific intensities, that I will then try and use in a similar fashion for the spectra individually.

[jsinkbaek]

Okay, I found an improved solution:

I start with 3284 detected lines in my master, all with S/N > 300.

• First I reject using flux/intensities like before. Here I decide on the limits manually such that it minimizes my standard deviation.

• Then I reject or merge lines in case of multiples, like before.

• Finally I do 5 sigma clipping (instead of normal rejection), which brings my std of wl-wl_ref down from an insane 417Å to 0.006Å.

With this I am left with 1182 unique lines for my master, which all have S/N > 300.

[jsinkbaek]

Interestingly enough, this also translates to an absolute RV error of 11 m/s.

[gmbrandt]

I have not used Redman 2014 before. I really have not spent much time at all in the game of tailoring line lists to a specific instrument. It is somewhat a can of worms.

You are the expert here for FIES/this line list investigation -- So I don't have any more to offer on that front. I will say that what you have outlined above seems reasonable to me.

Also, I am closing this issue because we are stepping out of the domain of things directly related to xwavecal. Feel free to email me, although I don't think I can be of much help in this particular line of investigation.

Thanks so much for using xwavecal!