mzechmeister / serval

calculate radial velocities from stellar spectra
MIT License
36 stars 9 forks source link

new DRS support #49

Closed saigrain closed 2 years ago

saigrain commented 2 years ago

Hello,

I'm wondering if you have any plans to introduce support for the new (ESPRESSO) DRS? As the full HARPS-N dataset has now been reprocessed with this pipeline, and the HARPS South data is a few months down the line, this would be quite useful. I had a half-hearted go at doing it myself but realised the format changes between the old and new DRS data products are more significant than I had initially realised.

Many thanks in advance

Suzanne

mzechmeister commented 2 years ago

Hi,

I can try. Can you share (publicly or privately) some test files with me? Say like 10 spectra?

Mathias

saigrain commented 2 years ago

Awesome. I will try to identify a suitable set of files that are in the public domain, so that there is no possible problem with me sharing them with you, and get back to you ASAP.

saigrain commented 2 years ago

Hi again,

I asked Xavier Dumusque about whether there are any data processed by the new DRS that are public and that you could use to get SERVAL running on those, and here is his reply:

"SERVAL could indeed be applied to the new DRS data. Note however that new DRS data (except the Sun) are only available to the HARPS-N consortium for the time being. However, we are currently discussing with INAF to make does data public at one points. Mathias could use the published solar data to test SERVAL, as those are derived with the new pipeline. He can download a few spectra here https://dace.unige.ch/observationSearch/?observationType=[%22solarSpectroscopy%22]"

I also asked him whether the format of the output files of the new DRS is documented somewhere, as there are some significant differences with the old one. He suggested to use the ESPRESSO pipeline manual (available here https://ftp.eso.org/pub/dfs/pipelines/instruments/espresso/espdr-pipeline-manual-2.3.3.pdf) - 200 pages, but hopefully you can find in there the information you will need. He also said he's very happy for you to ask him any questions you may have.

Thanks for agreeing to have a look into this, let me know how you get on and whether I can help in any way.

Suzanne

saigrain commented 2 years ago

Xavier provided the following additional information which may be helpful:

Here is a brief description of the s2d file S2D The main difference with the old pipeline e2ds products is that the wavelength solution is know given as an output (and not in the header), and this wavelength solution is already corrected from instrumental drift (it was not the case before).

mzechmeister commented 2 years ago

Ok, serval is now running on the new DRS for HARPN. (Don't be confused that the modified file is inst_HARPS.py; inst_HARPN.py just uses that script, too.)

I tested it with eight random spectra.

> cat sun.lis
/data/HARPN/DRS/sun/2015-08-17/r.HARPN.2015-08-17T11-03-36.560_S2D_A.fits
/data/HARPN/DRS/sun/2016-07-30/r.HARPN.2016-07-30T10-55-07.796_S2D_A.fits
/data/HARPN/DRS/sun/2016-09-30/r.HARPN.2016-09-30T13-44-33.172_S2D_A.fits
/data/HARPN/DRS/sun/2016-10-02/r.HARPN.2016-10-02T12-22-53.147_S2D_A.fits
/data/HARPN/DRS/sun/2017-03-09/r.HARPN.2017-03-09T12-53-43.004_S2D_A.fits
/data/HARPN/DRS/sun/2017-04-16/r.HARPN.2017-04-16T15-09-10.830_S2D_A.fits
/data/HARPN/DRS/sun/2017-06-04/r.HARPN.2017-06-04T15-34-42.027_S2D_A.fits
/data/HARPN/DRS/sun/2018-04-16/r.HARPN.2018-04-16T15-08-57.357_S2D_A.fits
>
> serval sun sun.lis -inst HARPN -snmax 460 -look

The solar HARPN spectra looks good. Formal uncertainties of ~1 m/s per order:

image

And final RV uncertainties are like 30 cm/s, nice. Yet, I was first surprised by the large scatter for the RVs (rms ~ 6.8 m/s). I thought it were problems with the barycentric correction. Thereby, I realised that the wavemap is already corrected for BERV and that the DRS RVs looks very similar.

srv sun -drs

image

Note, I undo the DRS BERV correction from the wavemap in https://github.com/mzechmeister/serval/blob/57c4bf51e368c011212a1f52193f8aacf131327a/src/inst_HARPS.py#L162 If -targ is not specified, the DRS BERV will be used again later. DRS BERV values are good and quick. But the past has shown, that it is safer to recompute BERV consistently (example unsafe DRS BERV: users specified no or different proper motions). However, the -targ option however does not accept the Sun, which is not resolved by simbad and anyway somewhat special. At least I computed quickly some BERVs for the Sun

>>> from astropy.coordinates import get_sun, EarthLocation
>>> from astropy.time import Time
>>> 
>>> bjd = [2457251.9575, 2457599.95161, 2457662.06932, 2457664.01261, 2457822.03411, 2457860.12814, 2457909.14581, 2458225.12798]
>>> 
>>> t = Time(bjd,  format='jd', scale='tdb')
>>> sun = get_sun(t)
>>> tng = EarthLocation.of_site("lapalma")
>>> sun.radial_velocity_correction( location=tng)
<Quantity [ 564.72575819,  447.86926282,  434.50564285,  582.52057852,
           -389.71130622, -676.7929822 , -433.19018654, -671.43068784] m / s>

If that is useful, it could also fiddled into serval. And here is a comparison with DRS BERV stored in sun.brv.dat image Some 3 m/s differences on top of 14 offset. I don't know the reason.

If your are happy with the changes, feel free to close this issue.

saigrain commented 2 years ago

Thanks for the quick turnaround! I'm testing it now.

One thing I noticed is that serval.py doesn't search for S2D files if you pass it a directory. Using a list gets around that, but it might be something to add before closing the issue?

Right now I'm also getting an error in read_spec.py at line 707 because I'm running python3 and str.decode() no longer exists. Did you test the changes under python2 only?

saigrain commented 2 years ago

PS: The scatter in the solar RVs is not completely crazy. There's a contribution from the solar system planets, mainly Jupiter with amplitude ~11m/s and period ~6 months (synodic period of Earth and Jupiter) and a contribution from activity with peak to peak amplitude of several m/s.

I do not know the origin of the difference between the BERV you compute and those from the DRS. Could be worth checking with the Geneva people.

mzechmeister commented 2 years ago

Yeah, I still use mostly python 2.7.

str.decode() is now fixed.

Note, instead of a file list, serval accepts also fifos, e.g.

serval sun <(ls /data/HARPN/DRS/sun/*/*S2D_A.fits) -inst HARPN

Here one could put custom wildcards or make other parsing on the fly with head, tail, awk, etc.. Furthermore, -nset and -n_excl are serval options for filtering the file input.

Anyway, for convenience S2D_A.fits is now included as a search pattern. https://github.com/mzechmeister/serval/blob/fe67f54b275d50a1dda44215c691e09e043d9dc4/src/inst_HARPS.py#L11

Since S2D spectra come in nite subfolders, wildcards for pathname are likely desired as well. Those needs to be quoted:

serval sun "/data/HARPN/DRS/sun/*" -inst HARPN
saigrain commented 2 years ago

It works! Thank you! Closing this issue.