In EWS all reduction steps, and the normalization by the photometry, can be executed in Pipelines which give you a few displays, but save all the intermediate steps for later diagnoses of reliability. We will now show the overall calls, list the steps actually carried out, and the intermediate files. See the UserGuide for more details. There are two pipelines, midiVisPipe if no photometric data is being used, and midiPipe when photometric data is to be used. Please look at these links for details about the calling sequences. Below I give a short summary. For the non-photometric case there is also a higher level pipeline midiProcess2Vis that applies midiVisPipe to both target data and calibrator data. It them combines the two results to give calibrated output spectra.
midiVisPipe, tag, filelist [,mask=mask] [,smooth=smooth] [,gsmooth=gsmooth] [,msmooth=msmooth][/twopass][,ave=ave][,maxopd=maxopd] midiPipe, tag ,filelist [,mask=mask] [...same parameters as midiVisPipe][,dSky=dSky]
midiPipe checks your input list a little, and then calls midiVisPipe to deal with the interferometric data, midiPhotoPipe to deal with the photometry, and a C-program, oirRedCal to combine the two into visibilities. midiVisPipe calls the following C-programs to process your raw interferometric data:
oir1dCompressData
oirFormFringes
oirRotateInsOpd
oirFGroupDelay
(optionally oirPowerDelay and another iteration of oirFGroupDelay)
oirFormFringes (2nd iteration)
oirRotateGroupDelay
oirAutoFlag
midiVisPipe calls the following C-programs to process your raw interferometric data: <\p>
oirChopPhotoImages oirMakePhotoSpectraIn the minimum specification of these scripts, there are only two required input parameters: a character string, tag, which gives a prefix that will be attached to all output files produced by the pipe. This will typically also include a directory. The second input parameter is a file list. So a typical midipipe routine call would be :
inputFiles = midiGui(dir=rawdatadir) tag = '/cool2/reduceddata/HD12345/firstTry' midipipe, tag, inputFilesmidiGui/s are routines to selected data files from a directory full of files.mask can be used to specify a mask file as described above. If you don't specify it, the system looks up some default masks based on the date of your observations. It is often better to create a new mask based on your calibrator observations using midiMakeMask This is done automatically in the higher level script: midiProcess2Vis I will describe twopass and the various smooth parameters later.
Now a description of the individual steps:
Compression
The first reduction step is compression, exactly as described earlier. The actual C-function is called oir1dCompressData, and the output is a filetag+'.compressed.fits'which is in the same format as regular data (and can be accessed with oirGetData()) except that the DATA1, DATA2.. arrays are only 1 pixel high. The actual shell call is:oir1dCompressData "file1 file2 file3" maskfilename tag.compressed.fitsClick here to see an example of the output.
Formation of fringes and high-pass filter oirFormFringes
In the second step the two interferometric channels are subtracted (they are 180 degrees out of phase), which reduces the background by about 90%. I assume the wavelength channels are sufficiently aligned. At this point I also apply a high-pass filter in time to reduce sky and instrumental backgrounds. For each pixel in the compressed data, I run a boxcar filter of given width in frame numbers over the data, and then subtract this smoothed version from the original version. The default boxcar width is 50 frames, but this can be overwritten by the smooth parameter in the call. If you use a very small smooth parameter (e.g. 2 or 3) some of the real interferometric signal will be removed by the filter, so your sensitivity goes down. But if your source is very weak and you discover that the group delay finding step (below) finds low frequency GARBAGE left over from the sky, then you should use a smaller value of smooth to suppress this. Note: be sure to use the same value of smooth for your target and for your calibrator! A second step to supress backgrounds is to fit and remove a polynomial in wavenumber from each spectrum. If the actual OPD is at least several wavelengths away from zero, the coherent fringe will appear as a sine wave in wavenumber, while the background will be smooth and nearly constant with wavenumber. So fitting a low order polynomial and subtracting it removes the background and little of the signal, unless you are close to zero OPD. In this case you remove signal too. Not a good idea. The parameter -averageOrder specifies which order of polynomial to remove. For non-zero OPD tracking, order=0 or 1 are pretty effective for weak sources. For zero OPD tracking, it is better to omit this parameter, so no subtraction is done. The output of this step is tag.fringes.fits, which is also a simple imagedata structure, with now only a single DATA1 data array.oirFormFringes tag.compressed.fits tag.fringes.fits smoothClick here to see an example of the output.
Remove Instrumental, known, OPD
The MIDI and VLTI hardware insert known OPD changes into the optical path, partly to keep the total OPD within reasonable bounds, partly to allow background removal. The effects of these changes needs to be removed from the data before further processing. These instrumental OPDs can be read out of the data arrays in human units (e.g. microns) withopd = midigetopd(inputtag)This works with almost all the data files during the reduction, since the OPD is copied from one to the next, but it is faster with data that is already compressed. The known instrumental OPD is removed by multiplying the data by exp(-ikD) with the program oirRotateInsOpd. There are no optional parameters.oirRotateInsOpd tag.fringes.fits tag.insopd.fitsThe output is again an imagedata structure, but the actual data are now complex, but represented by two float numbers. To convert the IDL structure to complex values, and transpose the array so that time is in the X-direction and channel in the Y-direction, which is usually nicer for display, use midiGetComplexrotdata = midiGetComplex(tag,'insopd') tvscl,abs(rotdata[indgen(1000),*])
Atmospheric group delay analysis oirFGroupDelay et al.
We next find the atmospheric group OPD: Take the FFT of the output of the previous step. Average several frames together to suppress the image peak and increase the S/N. Then find the OPD of the peak in the absolute value and write it down for later use. This is done by: oirFGroupDelay infile outfile -smooth gsmooth -medsmooth medsmooth ... . There are more parameters, but you should look at the users manual for these. There is also a multipass version of oirFGroupDelay described below. gsmooth defines the time interval (in seconds) over which data is averaged before looking for the OPD. The best choice of gsmooth depends on the source and the weather. Strong sources in rapidly varying OPD weather should use a small value of gsmooth, but probably not smaller than aboud 0.06 sec. Alternatively weak sources in good weather should use a larger value say 0.5 s. For weak sources in bad weather, you're out of luck. medsmooth defines the time interval (seconds) of a median smoothing applied a postiori to the OPD solutions before storing them. This processing is based on the assumption that the atmopheric OPD cannot suddenly jump 20 microns, so that if you see such a jump in the solutions, it must be a glitch or an error where the signal is too weak. Good values for medsmooth are probably about 0.4 s for strong sources and 2-3 seconds for weak sources. The result, tag.groupdelay.fits contains two interesting tables. The main data table is a pseudocomplex data set, with delay function (complex) as a function of delay (the actual value of delay as a function of pixel number is too difficult to explain here).
There is also a DELAY table in this file, which you can access with:
delay = midiGetDelay(tag)or more primatively: delay=oirgetdelay(tag+'.groupdelay.fits')This table has two columns: the time in the usual units, and the estimated Group Delay in meters, but if you use midiGetDelay you get the answer in microns. It is very useful to plot Group Delay against tracking OPD, because if they don't follow each other, the on-line tracking mechanism failed, and your data are probably worthless. So you can look at the data that the delay finding algorithm tried to use by
df = midigetcomplex(tag+'groupdelay') tvscl,abs(df[1000+indgen(1000),*]) tvscl,rebin(float(df[1000+indgen(5000),*]), 1000, 512)Note that the data returned by midiGetComplex has not been smoothed in the time direction, so the S/N is lower and the time resolution higher than in the search algorithm. In particularly you will also see the image peak zooming up and down in delay as the piezos scan. If you average these in time they will get much weaker. You can try:dfsmooth = csmooth2(df,s) tvscl,rebin(float(dfsmooth[1000+indgen(5000),*]), 1000, 512)(s is the sigma of a gaussian smoothing kernel). It is a very good idea to look at these delay functions. If you do not see a smooth, well behaved peak, your data are probably worthless!
MultiPass Group Delay Processing
If we have to search a large OPD region for a weak source, we may easily pick up a noise peak. To avoid this problem we can try to first localize the peak crudely by averaging the delay function (the FFT of the spectrum) over times longer than the atmospheric coherence time. We average the absolute value of the df (squared); otherwise the varying phase would decimate the signal. We record the peak of this heavily smoothed peak of the df as a function of time. Then we run a new search for the df peak, averaging coherently over times of order the coherence time, but we restrict the search to be near the just found peak of the incoherently smoothed function. We actually run three programs: 1. oirFGroupDelay. This only performs FFTs of the input spectra to form the delay function. 2. oirPowerDelay. This first averages the df coherently over a time gsmooth and then incoherently over a time ampsmooth and finally runs a median smooth with interval medsmooth. For weak sources typical values are gsmooth~0.5sec, ampsmooth~8 sec, medsmooth~16 sec. 3. oirFGroupDelay again. Here gsmooth and medsmooth should have values as specified in the single pass mode, but additionally the parameter maxopd should be specified with a value of ~2 wavelengths, to limit the OPD search to regions near the peak of step 2.
Remove atmospheric OPDs/align frames:oirRotateGroupDelay
Now we remove the Group Delay and the instrumental delay simultaneously from the original fringe data. At the same time, we determine and remove phase fluctuations caused by water vapor, which are nearly independent of wavelength, and hence not found by oirFGroupDelay.oirRotateGroupDelay tag.fringes.fits tag.groupdelay.fits tag.ungroupdelay.fits[-smooth smooth]The only optional input is smooth, which specifies the time interval over which to smooth data before estimating the water vapor phase. You might choose the same value as gsmooth. The output is once again pseudocomplex.
Editing: oirAutoFlag
Flag individual frames that seem fishy. This is not done on the basis of the estimated fringe amplitude (since dropping low amplitude frames biases the result), but on two other criteria: The distance between the tracking OPD and the Group Delay OPD, and the existence of OPD jumps in the Group Delay estimate. A large difference between the tracking OPD and the Group Delay OPD indicates that the on-line MIDI system had not (yet) found the fringe. The maximum allowable difference depends on the spectral resolution of MIDI being used: spectral resolution divided by the OPD difference, measured in compatable units, should be less than one radian. In practice this means that ΔOPD should be less than 100 μm for the PRISM and about 500μm for the GRISM. In particular, there are 50 or so frames at the beginning of each tracking run where we deliberately set the OPD way off, and these should be rejected from the estimation. Secondly, if there is a big OPD jump from one frame to the next, caused e.g. by the UT auto focusing mechanism, then the data near that jump is probably garbage, because the OPD was varying substantially during single integrations. Thus the call goes:oirAutoFlag tag.ungroupdelay.fits tag.delayfile.fits tag.flag.fits deltaOpd jumpOpd sideDropwhere deltaOpd is the allowable maximum ΔOPD (microns), jumpOpd is the maximum allowable jump (microns) and sideDrop is the number of points to drop on each side of a jump. The current defaults are:deltaOpd= 150 microns (PRISM) deltaOpd= 800 microns (GRISM) jumpOpd = 10 microns (from one frame to the next) sideDrop= 1The output file tag+'.flag.fits' contains a FLAG table as described in the original OIR FITS document, and is essentially a list of times that contain data flagged as nogood, and the reason why.
Coherent average oirAverageVis
Finally all "good" frames are averaged together to form the estimated complex visibility.oirAverageVis tag.ungroupdelay.fits tag.flag.fits tag.corr.fitsThe output: tag+'.corr.fits' is in OI_VISIBILITY format, at discussed above.
Photometry
Average Photometric Images oirChopPhotoImages
Now we have to compute the photometry in order to normalize the visibility. This takes place:oirChopPhotoImages -in "infile1 infile2" -out outFile [...etc] For chopped photometric files with either the A or B shutter closed, the sky frames and target frames are averaged separately, and the average sky is subtracted from the average target frames. If requested, an edge-of-slit estimate of the residual sky background is made and subtracted.Extract Photometric Spectra oirMakePhotoSpectra
oirMakePhotoSpectra -A Afile -B Bfile -mask maskFile -out outFile The data from oirChopPhotoImages for both telescopes are compressed in the Y-direction n various ways. This yields a 1-dimension imagingdata file with 6-rows, each for 1 variation, plus 6-rows of error estimates, not particularly reliable, in a file called tag+'.photometry.fits':Note that the maskfile should be identical to the one used in computing visibilities!.
- UnMasked Flux (ADU/s/channel) for AOPEN shutter position! (Average of two channels I1 + I2)
- UnMasked Flux (ADU/s/channel) for BOPEN shutter position
- UnMasked Flux GeometricMean of A and B (i.e. of rows 1 and 2). A and B are multiplied pixel by pixel, sqrted, then added in y-dir
- Masked Flux (ADU/s/channel) for AOPEN shutter position
- Masked Flux (ADU/s/channel) for BOPEN shutter position
- Masked Flux GeometricMean of A and B (i.e. of rows 4 and 5).
Photometric Normalization oirRedCal
Finally! If photometric data is available, the correlated flux computed from oirAverageVis is divided by the masked, geometric mean photometry (row 6 of tag.photometry.fits.oirRedCal tagwhere the tag should be enough to point to the two files needed. The output is automatically named tag.redcal.fits. At this point midiPipe is finished and will have put two plots on the IDL screen:
midi(Vis)Pipe,caltag,calfiles...which produces caltag.redcal.fits and all the other files. We then type:
midiCalibrate, tag, calTag [,calDatabase=cdb][,calFlux10=cf10] [,diam=d] [,/print] [/nophot]This calls the C-program oirCalibrateVis which divides a lot of things by a lot of other things and produces several output files: If photometry data is being used, these FITS files are created: tag.calvis.fits (a VISIBILITY file with calibrated visibilities)
.
Resolution:If you know that your calibrator is somewhat resolved, you can specify
its diameter as diam=diameter, in millarcsec, and a uniform disk model
will be used during the calibration. If the calibrator is in the
EWS/v.Boekel database, its diameter will be taken from there.
The program produces a large postscript files tag+'calPlots.ps' which plots just about all the calibrated
and uncalibrated quantities against wavelength, with somewhat
unreliable error estimates. Look at these carefully!
The /print option causes the program to produces about 12
ASCII files named with tag----.dat, which contain the same
information as the plots in case you don't want to deal with the
binary FITS files for further processing of your data.