Evaluation of FEROS pipeline            first draft

H. Hensberge, Royal Observatory of Belgium

Purpose

At moderate and high signal-to-noise ratio, the quality of spectra is limited by systematic errors rather than by random errors. In order to evaluate the limiting factors of the quality of the final spectra, the basic characteristics of the original data and their sensitivity to the pipeline software were checked. The latter was checked partly on artificial data whenever useful.

This report summarises:

1. characteristics of the data which may invalidate specific reduction steps
2. efficiency of subsequent reduction steps of the pipeline (bad column correction, order centring, background subtraction, order extraction, unblazing & flat-fielding, wavelength calibration, rebinning & merging) and easy-to-implement improvements to enhance their efficiency.
3. Our data reduction philosophy in contrast to the FEROS pipeline

The existence of the pipeline is very much appreciated with relation to its official purpose (quick-look data reduction). However, many researchers appear to use the spectra it produces, albeit with applying some parameter updates, as end-products.

Hence, there may be some interest in discussing it in more detail.

The data

The data were obtained in the object+sky mode, with the object in fiber 1 and the sky in fiber 2. For a given spectral order, fiber 1 is on the raw frame left of fiber 2. Calibration images included 6 series of 10 flat-fields with 34 s integration times and 4 series with increasingly shorter integration times of 16 s, 7 s, 3 s and 1 s (allowing to study non-linear responses); 6 series of 10 ThAr spectra, most with integration times of 6s or 7s in order to avoid the saturated lines at redder wavelengths obtained with the officially suggested integration times and regularly individual ThAr spectra during the night. All data were obtained on Feb. 17-19 and Feb. 28 - March 4, 2001. Furthermore, a number of bias frames and long-exposure dark frames were obtained.

1. Science frames

Most of the intrinsic characteristics of the images are easier visualised in the flat-field and wavelength calibration frames and are discussed in those sections. One difference between science and calibration frames to be kept in mind is that in a science frame fiber 1 receives usually much more light than fiber 2 while they receive a comparable amount of light in calibration unit spectra. This difference implies a different sensitivity to systematic errors produced by background and extraction algorithms.

2. ThAr lamp spectra

The light from each fiber is split in two slices to improve the spectral resolution.

Magnification of the frame around strong, unsaturated lines, up to the level where individual pixels are easily discerned, shows that the 2 slices of the same fiber produce spectra that are shifted relative to each other. The slices produce a cross-order profile that is resolved (two maxima are seen in each fiber, see Table 1). In the wavelength direction, the maximum intensity of the sharp ThAr lines in one slice is clearly shifted over several tenths of a pixel with respect to the other slice. The effect is visible in both fibers. It is unimportant as far as resolution is considered, but invalidates algorithms based on the use of the cross-order profile. Indeed, the local cross-order profile is biased strongly by the strong intensity gradient with wavelength. Therefore, the optimal extraction technique will give false results, especially in sharp features and with increasing signal-to-noise ratio. Also, the detection algorithm for cosmic rays (based on the cross-order profile and an error model that only counts read-out and random noise) leads to false detections, especially in the case of large spectral gradients. In particular, false rejections lead finally to artificial asymmetries in strong, sharp absorption lines of late-type stars.

The positions of the ThAr lines are very stable, at least over a few weeks but presumably also over longer time intervals. Changes are slow in time except for a cyclic 300 m/s semi-amplitude variation with a characteristic time of some 35 minutes. Might this be due to a cycle in the temperature control of the coudé room? Anyway, terrestrial absorption lines in the science spectra allow to reduce that effect up to the 100 m/s level and, moreover, the object/calibration mode is available for applications needing higher accuracy.

Finally, faint ghost orders are present (at the level of 100 ADU, from optical reflections?) which cross the aimed spectral orders. They do not have any significant influence on the wavelength calibration itself, but would be disturbing in early-type object spectra taken in object/calibration mode. This is due to the fact that the ghosts intensity is not negligible with respect to the depth of rotationally broadened absorption lines in such stars.

3. Flat-field lamp spectra

Looking at a cross-cut through a raw image at constant Y coordinate (i.e. a row), on an intensity scale optimised for the background, shows that:

1. the spectra of the two fibers overlap slightly (in the bluer orders, the inter-fiber background is systematically slightly higher than the inter-order background).
2. in the red, the overlap between the spectra from fiber 1 in a given order and fiber 2 in the next order to the red is much stronger than the inter-fiber overlap (here, the inter-fiber background is much lower than the inter-order background).

This format implies that the correctness of the extracted spectra depends more on the background model and on the location and length of the extraction slit than if the fibers and orders were projected with sufficient separation.

Figure  Ratio of two extracted flat-field images in different orders. The projection of the orders differed by about half a pixel in cross-order direction. Notice the different intensity ratios in the overlap regions.

The projection of the orders on the detector changes from night to night only by few tenths of a pixel at most and stays within the same pixel for very long times. However, division of an extracted flat-field by one taken during another night (Fig. 1), with a slightly different projection of the orders on the detector, leads to an interesting correlation with respect to the mismatch of the intensity in the order overlap region. The more the order is displaced, the worse the intensities fit in this 'unblazed' flat-field. If there is a problem unblazing flat-fields, then one may expect problems when unblazing object spectra. Analysis shows that this mismatch is identical for each pair of subsequent orders when expressed in terms of the relative position in each blaze profile (an effect that is transported by cross-dispersion to all orders).Hence, the effect is caused by a temporal change in the blaze function produced possibly by the same (optical?) effect that shifts the spectral orders to a slightly different projected position on the detector. Merging of spectral orders without removal of this kind of effect ,which grows in many orders to roughly 5% for a displacement of half a pixel, produces bad continua for stellar spectra.

Figure  The shape of the 'sinusoidal' disturbance as a function of ln(wavelength). The global intensity distribution of the flat-field spectra, scaled and shifted, is overplotted in red.

Alternatively, a series of consecutive flat-fields may be rectified by division through a specific flat-field from that series. In that case, the order overlap is perfect, but the intensity as a function of wavelength is not as smooth as expected in the wavelength range 420-650 nm. Superposed is a sinusoidal function with a 'periodicity' slightly longer than the free spectral range in the concerned spectral orders. This 'periodicity' changes proportional to the square of the wavelength. (while). The sinusoidal function in the rectified flat-fields dies out gradually to a negligible amplitude for the later members of the series i.e. it stabilises after the lamp has been lit for a sufficient time (something like half an hour!). However, comparison with the science spectra strongly suggests that even then the flat-fields are modulated by a sinusoidal intensity pattern, which produces a wavy continuum in the rectified, unblazed object spectra. It is interesting that the amplitude of this sinusoidal pattern rises steeply from the blue, to reach a maximum at 460-500 nm and then disappears gradually (Fig. 2, 3). It appears that this pattern originates exactly in the wavelength region where the light of the stronger flat-field lamp is attenuated by a combination of three filters. (Note that the flat-field is produced by two identical lamps, one of which is five times stronger with the red light blocked out by the filters). The pattern is similar, but not at all identical in the two fibers (Fig. 3). Understanding this problem correctly would significantly ease the rectification of stellar spectra and, especially, the disentangling of rotationally broadened metal-lines in early-type double-lined spectroscopic binaries.

Figure  Comparison of the 'sinusoidal' disturbance as seen in fiber 1 (black) and fiber 2 (red) in the ratio of two flat-fields from a consecutive series of exposures.

In addition to the correctly identified bad columns, one can identify a bunch of lower sensitivity columns in the red part of the detector. Also, in bias frames one identifies more hot pixels and (parts) of columns, but these are only slightly abnormal in comparison to the listed bad columns. Interesting is, however, that read-out noise and bias level derived by excluding these marginally bad columns are constant to a much higher degree and consistent over the whole frame (see Appendix 1).

The study of flat-fields taken with integration times from 1 s to 34 s (the fainter ones corresponding to exposure levels for bright stars) show that the columns 371 and 377 (dark), 1015, 1349+1350 introduce strong non-linear effects while the columns 697, 951, 1113+1115, 1752 don't. The effects of the latter ones can be flat-fielded away if one takes care of a possible small difference in the position of the orders relative to these columns in the science frame with respect to the flat-field frame. Even a small, sub-pixel offset is sufficient to pass on the effect to a different row as the orders are almost parallel to the detector columns.

More details and figures on the change of the blaze profile and the removal of bad-column blemishes in the case of slightly different projected science and flat-field orders are described in a contribution to ADASS XI (De Cuyper & Hensberge). A postscript version is available by anonymous ftp to midas.oma.be (directory feros).

Evaluation of the FEROS pipeline

Bad-column correction

In order to avoid problems with the detection of the orders, the FEROS pipeline 'corrects' the data in listed bad columns (colmean.prg). The applied interpolation procedure is sufficient for its purpose, but leaves residuals in the extracted spectra and flat-field. These cancel out in the unblazing process, except for non-linear responses, if the orders in the science and flat-field frames were exactly at the same position. If this is not the case, they only cancel out if the unblazing is preceded by a division of star and flat-field spectra in the coordinate space that expresses the position of the centre of the order relative to the bad column. It was checked and proven that such a procedure works perfectly. The procedure is based on:

(a) the identification of the bad columns,

(b) the descriptors fit{num} describing where the order centres are, and 

1. a model for each of the bad-column crossings mentioned above.

Presently, the code was developed simply as a MIDAS procedure.

Suggestion: The stability of the spectral format allows to communicate to the observers before their observing run which spectral regions are heavily affected by non-linear detector response. In some cases, it could be useful to know even before applying that the 414-420 nm region suffers significantly from non-linear detector response. If that region is important for a given programme, flat-fields should also be taken at the level of the stellar exposure in that wavelength region.

Order centring

The pipeline 

1. looks for the orders (2 fibers) at Y=0 by cross-correlating with a pre-defined template template.bdf;
2. follows the order to lower resp. higher Y and cross-correlates with the observed profile at Y=0 to determine the offset in X relative to X at Y=0.
3. fits each order with a polynomial of degree n-1 (keyword FIT_DEG = n)

The centring algorithm used is the gravity method.

Table   template.bdf (33 elements from top row left to bottom row right)

0.01035 

0.07539

0.05509

1.47449

2.38302

2.59481



2.07112

1.53611

1.76976

2.54428

2.81448

2.24283



1.15652

0.24075

0.02000

0.00555

0.00395

0.00503



0.01377

0.08718

0.51263

1.34174

2.25746

2.38183



1.66776

1.28949

1.72414

2.61115

2.67828

1.90487



0.84115

0.12834

0.01568









The success of this procedure depends 

1. on the quality of template.bdf (any centring error made in this step survives the next, differential step);
2. on the quality of the reference cross-cuts at Y=0 (sometimes containing a 'corrected' bad column); 
3. on the random noise in each order separately (commonly the blue orders are very noisy)
4. on the systematic discretisation errors related to the gravity method (these are 0.15 pixel at most for the cross-order profiles seen in FEROS)

The template.bdf image available at La Silla does neither reproduce the separation between the two fibers, nor the separation between the two overlapping slices in the same fiber (Fig. 4). To reduce centring mismatch errors, I applied instead the template specified in Table 1.

Figure  The template constructed from the flay-field exposures (full line) compared to the FEROS template (dashed line). Relative positions are in mm.

The current version of the pipeline contains, unfortunately, some software bugs producing the following errors:

1. a sign error in center_one_order.prg produces wrong sub-pixel fractions for half of the cases
2. a wrong start value in cuts.bdf relates the selected cross-cuts to wrong rows in the original image (a miscount without serious consequences)
3. quite exceptionally, after correction of these bugs, a position measurement is wrong by just one pixel. However, I did not yet locate the origin of it, but it is easy to eliminate such high-residual points from the final fit.

Discretisation related errors due to the gravity method can easily be removed in the case that the cross-order profile is stable. One computes for FEROS the following correction (easy to enter in center_all_order.prg):

0.5 d  (1.195 + 2.17 abs(d))  (1 - 2 abs(d))  /  (1 - abs(d))

where d is the pixel fraction of the order centre,   max[abs(d)] = 0.5

Discretisation errors are identified by plotting the residuals of the fit against the pixel fraction on which the order centre is positioned (i.e. X - int(X) where  X  is the position of the order centre in pixel units)

Figure  Evaluation of centring errors as a consequence of choice of template and the gravity method. See text for detailed explanation.

A quantitative evaluation of the involved errors, on artificial data, is shown in Fig. 5. The artificial data were produced as 100 discretisations of the high-resolution cross-order template derived from the flat-field observations, shifting the bins with steps of 0.01 pixel relative to the template (the template of Table 1 is just one of these discretisations). One set of data had flux in the two fibers (black curves), and another set (red curves) had only flux in the first fiber (as often is the case in object/sky mode science frames). Centring was performed with the FEROS template (with and without bug) and with our template. The errors in the centring are shown (in pixel units). Our template gives identical results for the two data set (red short-dashed line coincides with black one) as the template for each fiber is 'perfect'. The deviations from 0 are purely due to differences in discretisation between the real data and the template. They are well represented by the formula given earlier in this subsection. The long-dashed curves show similar discretisation effects, but also an offset. Due to a different template mismatch (with the one-fiber or the two-fiber data), the offset is different in both cases. The full lines show the effect of the sign-bug in the code. One should compare this centering errors with the centring accuracy of 0.01 pixel (rms) obtained finally in our data reduction.

Background removal

With the default settings of the FEROS keywords, the subtracted background is overestimated. This is easily seen at the end of the pipeline reduction since the spectra returned for the sky fiber display negative values at most wavelengths. Another obvious way of checking this is to plot a cross-cut of the original frame at a scale that shows the 'feet' of the orders, and overplot the derived background model.

Most of the damage can be avoided using BG_MEDIANX = BGMEDIANY = 1. However, more is needed to treat the fiber and order overlap correctly.

Order extraction

The mentioned lack of alignment in wavelength between the two slices of the same fiber lead us to use only a simple linear extraction (method 's'). Since we have multiple exposures, we prefer to perform the cosmic ray detection using this multiplicity rather than relying on the cross-order profile. Hence, the discussion of the order extraction is restricted to the simplest extraction method.

The extraction slit used by the pipeline consists of an integer number of pixels, centred as closely as possible. No pixel fractions are considered; hence the centring cannot be perfect and the position of the slit changes discontinuously with a whole pixel after some cross-cuts. No slit contains all flux of a given order/fiber without including unwanted flux of the other fiber or the next order, since the wings of these cross-order distributions overlap. Hence, it is unavoidable that a varying amount of flux is missed in different cross-cuts. The discretisation error introduced thereby can be visualised by changing the descriptors FIT{num} by half a pixel to the left and the right and then inspecting the ratio of such extractions. The discontinuous slit steps then occur at different cross-cuts. By visualising the ratio of the extractions, one sees a pattern that varies cyclically each time the centre of the order in a cross-cut shifts over exactly one pixel. Hence, the produced pattern has a high frequency where the order deviates significantly from vertical and has a low frequency where the order is vertical.

The consequences of the above on the final extraction depend on the slit length and the separation of the orders/fibers. Ideally, the slit should be as long as possible but exclude unwanted flux. My analysis shows that a slit length of 17 pixels has to be preferred for the majority of the orders (physical spectral order 63 to 33 i.e. the 31 bluest orders).However, when the order separation is smaller than the fiber separation a shorter extraction slit is better (15 pixels for spectral orders 32 to 30, 13 pixels for spectral orders 29 to 25). With a fixed slit for all orders, there are wavelengths where the amplitude of the discretisation pattern will exceed the 1% level. The standard setting PROFILE_W=15 is an acceptable compromise when using a fixed slit length for all orders.

It is of some interest to think a bit deeper about centring, especially in the case of asymmetrical cross-order profiles. The 'best' extraction is rather the one that includes as much of the wanted flux (without disturbance of unwanted flux). Or, alternatively, one that includes the same fraction of the flux and at the same time minimises the random noise (the latter is what optimal extraction aims at, but falls beyond our discussion for reasons mentioned earlier). That is not necessarily what is obtained in the case of mathematical centring. For example, when the spacing between the orders is much larger than between the fibers (as it is in the blue), it would not harm to extract flux from further away at the side of the cross-order profile that does not face the other fiber, but one should not go to near to the competing fiber on the other side. Finally, the length of the slit that minimises the discretisation patterns discussed above is larger than the slit that would optimise the random signal-to-noise and is less suited for extraction of low S/N data. The extraction with this large slit is also susceptible to background errors. When the background model is bad, systematic errors of one are another kind (due to background or due to discretisation) are unavoidable.

Technically, the extraction is prepared by producing a set of straightened images, with information on the background-subtracted flux and the pixel fraction of the order centre. The latter information is wrong for orders that run out of the frame i.e. do not reach Y=1. Since that happens only in the bluest spectral orders (62, 63), which do not contain useful flux in my case, I did not search for the reason.

Unblazing & flat-fielding

The extracted science spectra are divided by extracted flat-field spectra. The purpose is multiple:

1. eliminate the blaze profile modulation of each spectral order
2. eliminate fringe patterns
3. eliminate multiplicative, local sensitivity variations on the detector

The latter two require that the spectral orders in both exposures are projected to identical positions. The fringe patterns are sufficiently broad to cancel well even if that condition is only approximately fulfilled. However, the situation is different in the case of bad pixels/columns (i.e. in the case of sharp disturbances). These only cancel in a coordinate system that expresses where the order is relative to the disturbance.

The main concern in echelle spectroscopy is, however, (a). As mentioned in the discussion of the characteristics of the flat-fields, the apparent blaze function changes with time, together with the projection of the spectral orders on the detector. Moreover, it is modulated by a 'periodical' sinusoidal pattern in the wavelength range where light of the stronger flat-field lamp is partly attenuated. Order merging cannot be performed correctly before these effects are removed. At present, I developed methods to treat these defects during the data reduction. The best approach would be to identify their origin and either eliminate or monitor them during the observing run.

Wavelength calibration 

FEROS provides two options for wavelength calibration:

1. object/calibration observing mode for the most precise velocity work

(b) object/sky mode and ThAr calibration lamp images in the two fibers

The discussion here is limited to (b). When describing the calibration frames, the high temporal frequency changes in the positions of the spectral lines have been described already.

The pipeline measures the position of the detected calibration lines and identifies a subset of unblended lines to perform a 2-D fit to their positions. A differential analysis on the positions of individual lines in several ThAr frames reveals a few lines that are systematically deviating, presumably due to unrecognised blemishes. In particular, the complete spectral order 34 should not be considered for the calibration of fiber 2 because of the effects of a bad column. Removal of these lines leads to final fit residuals with a rms below 2.5E-3 A (6E-2 pixel at 5400 A) Which lines to remove will to some extent depend on the positions of the orders when detector blemishes are involved. Hence, it might be more efficient to inform the user about the problem and let him treat it than to hard-code it in the pipeline - unless experience shows that the stability of the spectrograph over long time intervals is causing always the same lines to suffer. The latter could be well the case but should be checked on the whole available database.

It would be preferable to minimise the residuals to the fit in log wavelength rather than in wavelength, as pixel size is proportional to log wavelength and not at all to wavelength. Minimisation in wavelength space is too tolerant for the blue relative to the red.

The high stability of the spectrograph keeps line shifts in the sub-pixel range during a typical observing run. Then, differential work relative to a template increases the robustness of the wavelength calibration. Differential line positions are successfully fitted by a second degree polynomial in the position along the order, independent of the order number. That is probably just a bit oversimplifying. However, with an rms from less than 5E-3 pixel for the most similar frames to 2E-2 pixel for the most shifted frames, the gain in consistency with respect to an independent absolute calibration (rms of the order of 5E-2 pixel) is obvious.

Figure  Global intensity distribution of the flat-field lamp spectra with ln(wavelength). The blaze profile in each order is removed.Vertical lines indicate start of the useful part of a spectral order (full line) and end of it (dashed line).

Rebinning and merging

Using the results of the wavelength calibration, the spectra are rebinned either in steps constant in wavelength or in the natural logarithm of wavelength. The latter choice has two advantages:

1. pixel size on the detector is in a very good approximation proportional to the logarithm of the wavelength. A step of 0.909E-5 in log wavelength preserves the resolution on the detector;
2. velocity shifts are proportional to the logarithm of the wavelength

Merging is done by combining data that refer to the same wavelength in overlapping orders. High-noise edges of orders are cut off (reject overlapping outside 1.5 times the free spectral range). If anything, I would advise to reject even more. Indeed, little is gained in terms of signal-to-noise when the intensity in one order is significantly lower than in the other and systematic errors in the less intense order, due either to unblazing or to background subtraction, may also be significantly higher. Fluxes of different orders are combined using linearly decreasing weights towards the edges of the orders.

An unnecessarily weak point in this approach (besides the lack of overlap consistency mentioned in the section about unblazing) is that one sometimes combines data in one order with data affected by detector blemishes (bad columns etc.) in another. Since bad columns and order positions are known, this can be avoided.

Presentation of the final spectra

The merged spectra are presented in relative intensity units i.e. the intensity of the science object relative to the flat-field spectra. It is left to the user to produce the normalised (continuum level = 1) object spectra used in astrophysical analyses.

Spectra in these relatively intensity units are far from 'flat'. Partly because of the energy distribution of the object (modified by the earth atmosphere), but dominantly by the bimodal intensity distribution of the flat-field spectra.  I derived an approximate intensity distribution in wavelength space for the flat-field spectra (Fig. 6). Dividing that away (i.e. equivalent to dividing the object spectra by normalised flat-field spectra) guarantees a consistent treatment of the general flat-field light distribution, with remaining temporal variations in that distribution being much smaller than changes implied by external observing conditions. This leads to object spectra that are much 'flatter' than the presently delivered ones. This eases the direct display of the spectra and helps the user in determining a good final normalisation for the object spectra.

Our reduction philosophy & concrete reduction scheme

The reduction chain I use is based on two main considerations:

1. work as much as possible differentially
2. use the available FEROS software as much as possible

Subsequent steps are as follows:

(1) Pre-reduction (bias subtraction, bad column treatment, removal of overscan pixels) as in the FEROS pipeline, average the consecutively taken flat-fields and ThAr frames.

(2) Construct template.bdf from the cross-order profile in several orders of the flat-field frames (around the location where the orders are vertical, allowing to average the local cross-order profile without loss of resolution). As different orders are discretised in a different way, the derived cross-order profiles can be combined to a higher resolution profile (as with dithered images).

(3) Apply the order centring procedure of the pipeline (debugged) on a master flat-field. However, all cross-cuts are used and not just one out of ten (note that the variables MAXLINES in echelle.h should be changed to 4102 to allow that; presently it is, arbitrarily, equal to 1000). Gravity discretisation errors are removed. A selection of cross-cuts is made in the table centers.tbl in order to avoid the use of low accuracy data in the right lower part of the image. The 2-D fit made is a bi-variate polynomial of degree 5 along the order and degree 7 in the order number (note that one has to increase the dimension of the output keys in the MIDAS software to fit such a polynomial). The residuals relative to this model have an rms lower than 0.01 pixel in all well-exposed orders! Some systematics in the residuals is detected locally because of two reasons: (a) bad columns and other blemishes, and (b) imperfections in CCD production (systematics at the level of 0.01 pixel in every subpart of 512 pixels in vertical direction show how the detector was build up).

The same procedure is applied to one reference spectrum of each star for a given observing run.

(4) Determine for all remaining flat-field frames the positions of the spectral orders relative to the chosen master by cross-correlation with that master in each cross-cut. The underlying reason for this choice is that the master should be a perfect template. The position of the cross-correlation peak is determined by taking into account the corrected for discretisation errors proposed by by David & Verschueren (1995, A&AS 111, 183). Few outliers are removed. The differential positions are fitted with a bi-variate polynomial of second degree both along the order and in order number ), but the dominating term is the constant. The contribution of all other terms remains below 0.03 pixel. The differential fit is used to construct the absolute positions and to save them in the FIT{num} descriptors needed further-on by the pipeline.

(5) Apply the same differential work for all stellar spectra relative to the stellar reference frame (the perfect template)  and for the stellar reference frame relative to the master flat-field (possibly slightly less perfect). The execution of (4) and (5) is presently time-consuming because it was, for testing purposes, programmed in MIDAS command language.

(6) Determine and subtract the background. A filtering/correction  procedure provides a frame with our background estimate at the positions where the pipeline will fit the background. In particular, the procedure applies a first-order correction to inter-fiber and inter-order intensities which were affected by nearby spectral orders. This frame is then subjected to the pipeline back/feros command for an improved smoothing and further compatibility with the pipeline software, using

BG_MEDIANX = BG_MEDIANY =1 and BG_STEPX = BG_STEPY = 2.

(7) Produce the straightened images with the pipeline software

(8) Extract the orders choosing parameters in accordance with the spacing of orders and fibers and with the reduction of discretisation errors. Use slits with a length of 13, 15 and 17 pixels (17 pixels used in spectral orders 63 to 33, 15 pixels in 32 to 30, 13 pixels in 29 to 25) and use order positions as determined earlier, and also shifted by half a pixel to the left and to the right, in order to be able to diminish the earlier discussed discretisation effects which vary with improper centring (by interpolating between the resulting extraction using the back_str_D images as weights).

(9) Wavelength calibration. Apply the pipeline procedure to the calibration frames

(some may be the average of series) . Compare the results of all frames differentially in order to detect and eliminate lines with inconsistent positions in many frames. Eliminate such lines from the measurement tables (presently done interactively). Fit the standard 2-D relation to a chosen reference calibration frame For all other calibration frames, fit the differences in position relative to the reference calibration by a second degree polynomial in the position along the order (the same polynomial in each order). Update the line positions in the measurements tables to their values in the reference frame plus the differential correction and fit the absolute relation to these updated positions in order to provide the needed dcoef descriptor

(10) Rebinning to log wavelength with a step of 0.909E-5 using the pipeline software.

The next steps are performed quite independent of the pipeline. These steps have been tested but some of them were not yet applied routinely to many spectra. In different steps, one goes forth and back to the spectra in individual order in pixel space or log wavelength space, as required.

(11) Prepare normalised flat-fields for unblazing. Two steps are necessary:

1. construction of  an image that describes (approximately) how the flat-field flux changes with wavelength or log-wavelength. It is done once and interactively. In principle, sensitivities to any factor that is present in flat-field spectra and object spectra should not enter in the constructed flux distribution, but factors that apply only to the flat-field spectra should.
2. removal of bad-column blemishes. If the position of the flat-field orders differs somewhat from the ones in the science frame, the blemishes in the flux distribution need to be modelled as a function of the distance of order centre to such a  bad column (using all available flats, which differ somewhat in the positions of the orders)
3. application of the models developed in (b) to all flat-field and object spectra.

(12) Unblazing by division of the object spectrum through a normalised flat-field. Check the order overlap consistency i.e. the end of the preceding order and the start of the next order must have consistent intensities. If not, compute the ratio of these fluxes in all overlap regions (over that part of the overlap that will be used effectively in the merging process). That ratio turns out to vary little and at most slowly with wavelength such that it can be represented by a linear function of spectral order number. Using this linear fit, a correction function is constructed that 'tilts' the intensities in the spectral orders to an unbiased overlap.

(13) Only in the case of multiple exposures (in my case often 2 or 3 consecutive exposures) which should be combined to a single higher signal-to-noise spectrum:

1. model the smooth low-frequency ratio of the spectra relative to the chosen reference exposure for the multiple exposure series
2. by inter-comparison of the exposures, taking (a) into consideration, detect cosmic rays in each of the spectra
3. combine the multiple exposures using proper weights. (a) is taken into account to preserve the global characteristics of one of the original exposures; cosmic rays are masked by giving the concerned pixels weight 0.

(14) Merge the spectral orders:

1. using the information about the rebinning available in the wstart and nptot descriptors of the rebinned image; 
2. specifying for each order the part of the overlap region to be considered, and avoiding, whenever possible, to deteriorate correct data by data from blemished parts of the detector.

The choice of the useful parts of the overlap region has been made interactively (once for all frames), taking into account:

1. relative intensities in the different orders ;
2. relative importance of sources of systematic errors;
3. blemishes visible in the extracted data (in many frames).

(15) Rectification of the spectra to a normalised continuum level

The way to perform this step was studied on a spectrum of a sharp-lined early B-type star that was observed earlier also with CASPEC (albeit over a shorter wavelength range, but with higher signal-to-noise ratio). This made it easier to define wavelengths where the continuum is seen. A spline function allowing some smoothing (in accord with the noise in the spectra) is fitted and the result inspected. Iteration on the continuum points is allowed (primarily in the wavelength region without CASPEC data).

Such a procedure confirms the superposition of a sinusoidal function on the unblazed stellar spectrum. This function is introduced by the flat-field spectrum and needs to be removed. Once such a function is removed (see the interactive procedure below) from a single spectrum, it will be straightforward to remove it from other spectra if the assumption or suspicion is valid that it is due to the flat-field spectrum and not at all to the stellar one. Indeed, the appropriate function for spectra unblazed with the same flat-field will be identical to the existing one, while for spectra unblazed with other flat fields the appropriate function is derived straightforwardly from the ratio of the concerned flat-field images.

It is easy to measure the sinusoidal function present in the ratio of two flat-field spectra since, as was discussed in the section on the flat-field data, at least the amplitude of the function changes in time. Due to the fact that the superimposed sinusoidal function is a small disturbance on a strong signal, the function seen in such a ratio is in first approximation also the function seen in an individual flat-field albeit with a different global amplitude. This is not anymore true if the phase changes (the phase seems to shift slightly in the start of each series, but the first frame can be skipped) or the functional form itself. At present, the shape of the sinusoidal function in an unblazed flat-field is used as an approximation to the sinusoidal function present in the unblazed stellar spectrum. The global amplitude of that function and the phase are used as free parameters and fixed by trial-and-error. The unblazed stellar spectrum is divided by such functions, while varying amplitude and phase. Wavelengths at which the intensity is supposed to reach the continuum (CASPEC spectra, and literature data at longer wavelengths) are visually inspected. Promising trials are normalised to a continuum equal to one, permitting further iteration on the definition of the data subjected to the spline fit for the continuum. The final choice for amplitude and phase is made by visual inspection of these test spectra.

Figure  The 'sinusoidal' disturbance in the continuum of NGC224-201 compared to the predicted correction function (red). Notice how hydrogen line profiles are affected.

An extract of the constructed function is presented in Fig. 7 (Fig. 2 shows the global shape), together with the sharp-lined early-B stellar spectrum, normalised without correcting for the 'sinusoidal' disturbance. The broader, deep absorption is Hß with a diffuse interstellar band in the red wing. It is obvious that the red, sinusoidal line represents the continuum better than the presently adopted one. It was not yet checked whether that functions, adapted as explained above, would remove the pattern superposed on other stellar spectra. In fact, with a forthcoming FEROS observing run, it is hoped to solve the problem in a more straightforward way as suggested in footnote 12.

Appendix 1   Summary of hot columns/pixels on the detector 

NORMAL AVERAGE ON CCD = 203.40, ON OVERSCAN = 203.33, rms=1.12

ALL DATA REFER TO THE MEDIAN OF 5 CONSECUTIVE BIAS FRAMES

---------------------------------------------------------------------

col    row        mean      rms        row     mean        rms      notes   

256    <1689   203.39    1.12      >1689   204.18    1.14    1689=220 

270    <1679   203.45    1.13      >1678   204.66    1.21  

271    all          203.50    1.14

272    all          203.42    1.10

273    all          203.55    1.12

274    all          203.56    1.13

369    <1501   203.40    1.07      >1999   204.03    1.14    irrelevant difference

385    all          203.51    1.12

386    <1617   203.73    1.15      >1619   232.13    2.36    FEROS      1618=324

387    <1617   203.44    1.10      >1630   203.41    1.15    1617=219, 1618-30=205.62

                                                                                                                          tail

388    <1611   203.37    1.16      >1639   205.32    1.25    1631=210 insignificant?

389    <1623   203.40    1.19      >1630   203.98    1.15    irrelevant difference

393    all         203.53    1.13                       

394    <1623   203.44    1.16      >1630   205.05    1.30    1623=217, 1624=215,

                                                                                            1625=215, 1626=214,

                                                                                            1627=212 1628=209

                                                                                            1629=210 1630=211

696    all         203.58    1.14

697    < 870   204.43    1.19       >869     259.67    2.59      FEROS

894    <1514   203.42    1.13      >1521   204.30    1.18     1514=217, 1515-20=208:

                                                                                             marginal difference

934    <1494   203.47    1.11      >1493   204.30    1.17

951    <1476   203.55    1.11      >1487   209.90    1.49      FEROS smooth connection

952    <1463   203.42    1.10      >1463   204.12    1.19      1463=217 marginal diff.

953    all          203.39    1.13

954    all          203.52    1.12 

967    <1451   203.41    1.13      >1468   204.45    1.20       marginal difference

1011    <1424   203.51    1.12      >1438   203.91    1.14      irrelevant difference

1012    all         203.49    1.12

1013    <1424   203.43    1.11      >1438   203.55    1.14      irrelevant difference

1014    <1424   203.43    1.10      >1438   205.14    1.26

1015    <1428   203.47    1.15      >1430   211.02    1.32      FEROS   1429=226

1016    all          203.41    1.12

1017    all          203.47    1.12

1018    <1424   203.51    1.10      >1438   204.50    1.19      marginal difference

1019    all          203.43    1.13

1020    <1636   203.45    1.11      >1635   204.51    1.23      marginal difference

1021    all          203.59    1.12

1022    all          203.39    1.12

1023    <1636   203.57    1.17      >1635   204.15    1.15      marginal difference

1112    <1401   203.38    1.14      >1400   204.34    1.18      marginal difference

1113    <1404   203.49    1.11      >1405   216.18    1.57      FEROS 1404=299,

                                                                                                             1405=222

1114    all          203.59    1.16

1115    <1405   203.51    1.13      >1407   210.64    1.45      FEROS 1406=214, 

                                                                                                             1407=255

1116    all          203.54    1.12

1120    all          203.64    1.15

1346    all          203.65    1.14

1347    all          203.81    1.12

1348    < 610    203.51    1.19       > 609   204.70    1.17      marginal difference

1349    < 609    203.72    1.05       > 611   265.33    1.39      FEROS 609=249,

                                                                                              610=332, 611=329 

                                                       not929-940                     929-940=259.33_rms1.83

                                                                                               410=226

1350    < 609    293.50    2.77       > 611   33570     54.7      609=24438, 610=65535,

                                                                                               611=65535 FEROS

                                                       not926-943                      927-940=29076_rms52.1

                                                                                               926=31253, 941=32941,

                                                                                               942=32728, 943=33233

                                                                                               410=12261

1752    <1530   203.45   1.10      >1531   207.74    1.37       FEROS 1530=272

-----------------------------------------------------------------------------

columns < about 60 have a slightly lower bias (203.33 rms 1.12)

columns > 1960 have a bias marginally different from the previous columns

--------------------------------------------------

compare with regions of undisturbed bias columns:

275- 368 : 203.39  rms 1.12

395- 695 : 203.39  rms 1.12

698- 893 : 203.40  rms 1.12

1024-1111 : 203.40  rms 1.12

1121-1345 : 203.41  rms 1.12

1351-1751 : 203.41  rms 1.12

1753-1960 : 203.41  rms 1.12

-------------------------------------------------   ******

compare with overscan region:   203.33  rms 1.12    **OK**

-------------------------------------------------   ******

The physical order number changes inversely proportional to the wavelength itself. Hence, this sinusoidal pattern is in different orders not in the same relative position with respect to the wavelength of maximum blaze intensity.

such a test reveals also possible centring problems i.e. the 'centred' slit should produce the highest flux

or to unrecognised blending; a definitive answer on the origin of these deviating positions could be given from comparing frames in which the wavelengths are significantly displaced

note that they are not removed due to a high residual in one frame, but because of lack of consistency in differential positions in many frames; this criterion is much stronger than a simple removal of apparent outliers in a single fit

they deteriorate the fit in regions where better data are available; anyway, these regions with extremely low flux are not used in the final spectra

) our fit has only 48 degrees of freedom instead of the 195 used by the pipeline; coupling the orders makes the fit towards the blue orders much more robust

The preparatory background procedure should be seen as a quick fix, not as a final solution, and is more sensitive to dark blemishes than one would wish. However, any remaining systematic errors are lower than when applying the pipeline straightforwardly.

All steps up to the extraction are done exactly as for flat-field and science frames rather than as done in the pipeline, but using the order positions derived from the flat-fields.

As the use of filters apparently introduces features on the scale of a few free spectral ranges, the construction of such an image requires iterations. But, once available, it should be useful as long as the present configuration of the instrument and calibration unit survives.

On unblazed flat-field spectra, such a correction eliminates to a good approximation the observed variation of the apparent blaze profile (between the considered flat-field and the one used for unblazing), giving thus some support to this method.

in the case of a single exposure, the detection and removal of cosmic rays is handled in the standard way by the pipeline, despite the remarks made earlier, but the resulting spectrum is compared in detail with the linear extraction using all pixels.

However, it is obvious that it would be preferable to measure that function empirically in a direct way, as stated earlier. Flat-fields taken without using the filter-affected lamp and stacked in a sufficient number of them in order to have the necessary flux in the blue could offer a way out. Note that the interest is only in relatively low frequency behaviour. Hence, the signal can be improved by filtering out high frequencies.

Strictly speaking, the correction must be made before the final normalisation that depends on the choice of continuum points. The comparison shown here is not yet optimal because the chosen continuum points should be re-determined on the non-normalised spectra corrected for the 'sinusoidal' disturbance.