Fitting Gaussian components to compact sources

Contents

  1. Fitting Gaussian components to compact sources
    1. Test Data
    2. Manual component-fitting
    3. Automatic component fitting
      1. Setting up the inputs - total intensity data
      2. Setting up the inputs - total intensity data
      3. Building AIPS procedure and loop
      4. SAD Outputs
    4. Unresolved masers
      1. Unresolved OH masers
      2. Fitting to polarized data cubes
    5. Simulations

The physical meaning of maser measurements can be interpreted as described in e.g.

Several AIPS tasks exist which fit 2-dimensional elliptical Gaussian components to patches of emission. This discussion concentrates on their use for maser images. If the emission has an instrisically Gaussian distribution (such as maser emission beamed through an approximately spherical cloud) and the patches of emission are unresolved or of dimensions comparable to the beamsize (say not more than 2-3 times larger) then the fitting can be very accurate. If masers originate from a very asymmetric structure such as a very elongated shock front or (partial) shell, such that emission in a single channel is very elongated or banana-like, then it won't work so well.

Theoretically, the position accuracy of a Gaussian component depends on the (beamsize/2)/(signal-to-noise ratio) (Condon 1997, PASP, 109, 166). This may work for very well-filled apertures (e.g. full track with VLA) but is unrealistic for poorly filled apertures and in any case hits a lower limit when calibration and pixellization errors dominate. The response in wide fields and sparsely filled uv planes was analysed for the NVSS (VLA snapshots) by Condon et al. 1998, AJ, 115, 1693 and for MERLIN by Richards et al. 1999, MNRAS, 306, 954 in Section 2.4.

Whatever the scaling factor for the uncertainties, the fact that the errors depend on

means that using the most uniform weighting and/or smallest possible restoring beam may not give the best position accuracy. If you are trying to separate close but distinct components then you need just sufficient resolution to achieve this, but not so that you increase the noise excessively.

The fitted component size is the size of the emission region convolved with the restoring beam. By deconvolving the beam, the size of bright components much smaller than the beam can be measured. As AIPS does not always do a good job of fitting to very ragged beams, I prefer to print out the undeconvolved size and then subtract a model beam from each component (if the resulting apparent size is smaller than the uncertainty then the component is unresolved). You should experiment with simulations or various methods if necessary. Masers are particularly suitable for this as the brighter and more highly beamed the maser, the more compact it appears, but the signal to noise is better so you still have a chance of resolving it.

As a rule of thumb, for sparse apertures like MERLIN and the EVN, multiplying the SAD errors by a factor of 3, or at least in the range 2-5 (depending on just how bad the dynamic range limitations are) usually gives realistic uncertainties, but you should check

  1. By use of formulae such as in the references above or of your own derivation;
  2. By using simulations or by comparing discrepancies in measurements of different hands of data for unpolarized sources.

The advantage of using SAD errors as the basis for error estimates is that they also include fitting errors. The two examples below outline strategies for resolved and unresolved masers, including fitting to polarized emission. The use of local noise-dependent cutoffs helps to avoid fitting to sidelobes although human checking is always necessary to be sure, and also to identify occasions (such as blended components) where a fit has failed and has to be supported by an input model or better boxing.

Note that this semi-empirical approach to accuracy is reliable if you are measuring many components and can average over the properties of groups of masers. If you are studying one or two individual sources in detail then some of the factors (such as the sensitivity of the fitting algorithm to the size of the box; non-Gaussianity of the source) should be analysed in more detail and other methods such as fitting a slice might be more appropriate. I describe JMFIT and its automated incarnation SAD; there are other JMFIT-like tasks such as IMFIT.

These tasks give you an option of attempting to resolve components or assuming that the beamsize is the upper size limit in all cases. The decision which to use can be based on trials but here we use a priori knowledge;

Test Data

You can try JMFIT and SAD on your own data cubes or single images, which will probably be better than the ones I prepared quickly, but if you don't have any of your own you can use:

Manual component-fitting

Get G23.389.ICL001.1 and load the channel at 74.8 KM/S (chan 42 in /home/user/datastore/EVN_LINE/G23.389.ICL001.1 ) into the TV, using functy 'lg' Set a box around the brightest spot using TVWIN and get an estimate of the peak to the nearest pixel using IMST. This is a useful sanity check as Gaussian component fitting can occasionally fail and give spurious results. Set up JMFIT; it will inherit the blc/trc from TVWIN: 

AIPS 1:task 'JMFIT'
AIPS 1:INNAME  'G23.389';INCL 'SUBIM'
AIPS 1:BLC 85 143 42;TRC 102 162 42

Initially leave the other settings default - this will fit a single resolved component. Note that 1 pixel = 1 mas This gave me 

uathac> JMFIT1: Component   1-Gaussian
uathac> JMFIT1:   Peak intensity    = 2.1029E+01 +/-  6.78E-03 JY/BEAM
uathac> JMFIT1:   Integral intensity= 2.7202E+01 +/-  1.40E-02 JANSKYS
uathac> JMFIT1:   X-position        =     92.022 +/-   0.0009 pixels
uathac> JMFIT1:   Y-position        =    152.562 +/-   0.0010 pixels
...
uathac> JMFIT1:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
uathac> JMFIT1:                   Deconvolution of component in pixels
uathac> JMFIT1:                        Nominal     minimum     maximum
uathac> JMFIT1:     Major ax             4.411       4.407       4.416
uathac> JMFIT1:     Minor ax             1.760       1.753       1.768
uathac> JMFIT1:     Pos ang              1.274       1.199       1.348

The semi-analytic formula for sparse arrays gives

Position_error = Noise/Peak * Dirty_Beam/2 * sqrt(Nant/2)

where Nant is the total number of antennas used; in this case 8. The dirty beam had a major axis of 12 mas. Measure the off-source noise in this channel (I got 6.25 mJy/beam in this channel for a 6-mas restoring beam).

This gives Position_error = Noise/Peak * 12

The 'formulaic' error in component size is sqrt(2) * Position_error, or about 0.005 mas, which is similar to half the (max-min) difference in component size reported by SAD. The apparent size of the component is >> the error and the integrated flux density is significantly higher than the peak, suggesting that the componet is genuinely resolved.

Repeat for the faintest patch of emission in this channel and experiment with setting NGAUS 1, 2 or 4. You may have to increase NITER. I fitted two distinct components separated by 6 mas but trying to fit more just gives spurious duplicates.

Experiment with other box sizes, models etc. However, if you try and box the region containing all the emission, JMF fails completely.

Automatic component fitting

Setting up the inputs - total intensity data

Setting up the inputs - total intensity data

If you did the Spectral Line VLBI tutorial you may already have a cube for G23.389. In this example we use a subimage of the 85 channels which contain masers, in the central (275x275) pixels; the exact size does not matter but you should either create a similar subimage or change the blc and trc values below suitably. Alternatively you can load G23.389.ICL001.1 from datastore/EVN_LINE.

Use HELP SAD for more information about how it works. It has its own jargon; an 'island' is a region contained within a given contour (there may be many islands within a single channel). Within an island, the task finds the highest peak, subtracts it and then looks for possible other emission down to user-set limits. These flux/noise based limits can produce less components that the maximum specified by NGAUSS (number of components fitted), unlike in JMFIT. SAD can be used to search large regions but it is useful to restrict this to the area where emission is likely to be found, if known, or at least to the CLEANed region.

  1. Finding the noise
    As spectral line interferometry maps are usually dynamic-range-limited, it is best to measure the noise in each channel to limit the faintest components fitted. Where you set the boxes can make a difference, especially if there might be genuine emission near the edges of an image. The solution used here is to choose the lowest rms in a choice of 4 regions near the edges of each channel (high rms in a box might indicate genuine emission). For these data the regions cover the N and S edges because the worst sidelobes are in these directions.

  2. Dynamic range limitation
    A further safeguard is to take the whole area to be searched in each channel and measure (-1*(min/max)) which gives the minimum percentage of the peak which produces a believable component (this is not a suitable criteria for e.g. V images with possible intrinsically negative emission).

  3. Iterative component fiinding
    It can be more reliable to subtract the brightest components first and then look down to progressively fainter flux cutoffs (CPARMs).

  4. Failed fits
    Bad (i.e. mathematically unsound) fits can be rejected using the DPARMs.

Building AIPS procedure and loop

This procedure is used to find the off-source noise and then the dynamic range limit for G23.389.ICL001.1. This is done separately for each of the 85 velocity planes, (275x275) pixels per plane at 1 mas/pixel, restoring beam 6-mas circular. As it is a subimage there is no need to avoid the edges when measuring the noise.

AIPS has i and j as built-in variables. In this procedure:

For more information see HELP POPSY and HELP PROC in AIPS or look them up in the AIPS Cookbook. 

proc DETNOISE (i)
blc 1 1 i; trc 137 50 i;imst;j=pixstd
blc 138 1 i; trc 275 50 i;imst
j = min (j, pixstd)
blc 1 226 i; trc 137 275 i;imst
j = min (j, pixstd)
blc 138 226 i; trc 275 275 i;imst
j = min (j, pixstd)
j = 5*j;
blc 0 0 i;trc 0 0 i;imst
j = max (j, -1*pix2val)
return;finish

The peak in cube (from the image header) is 21.05 Jy/beam.1 The first and second pass cut-offs (cparm(1), (2)) are chosen empirically, the same for all channels. The lowest cutoff (cparm(3) and dparm(1), (2)) are set to the noise/dynamic range limit per channel j.

Note that the last IMST in the procedure uses the same box size as is used in SAD to search for components. If you have an image which is much larger than the total extent of the maser region blc(1),(2) and trc(1),(2) should be set to just enclose the maser region (found using eg. SQASH or MOMNT or KNTR).

DOALL refers to the number of peaks likely to be found within a single island. Here, the elongation of some components in the brightest channel suggests that this is multiple, and so we use DOALL 1; however this can lead to artificial fragmentation of bright components and is one of the things which has to be decided from experience with any particular dataset.

For masers likely to be resolved, DPARM(4),(5) and possibly (6) should be slightly larger than the maximum feasible size. All other parameters are straightforward; in some cases extreme values are set to disable selection e.g. we don't use the DPARMs to reject components with high residuals as we have DOALL 1. 

AIPS 1:task 'SAD'
AIPS 1:blc 0 0 i; trc 0 0 i
AIPS 1:inver i
AIPS 1:ngauss 100
AIPS 1:cparm 5 0.5 j
AIPS 1:docrt 0
AIPS 1:sort 'S'
AIPS 1:outver -1;doresid -1
AIPS 1:doall 1;dowid 1
AIPS 1:dparm j,j,1000,20,20,20,1000

To test on a few channels e.g. around the peak, 

AIPS 1:dowait 1
AIPS 1:for i=41:45;detnoise(i);inver i;cparm(3)=j;dparm(1)=j;dparm(2)=j;go sad;end

AIPS 1:for i=15:19;detnoise(i);inver i;cparm(3)=j;dparm(1)=j;dparm(2)=j;go sad;end

SAD Outputs

The fitted parameters are written to an MF table. inver i produces a separate components table for each plane of the cube. This is useful for overplotting components or if you also want to measure polarised flux (see below) but slightly complicates writing components out; use inver 1 to write all to the same MF table.

Note: If you are re-running SAD it will not re-find existing components, so delete old MF tables first.

MF tables can be inspected using the task MFPRT; however this does not give all parameters. PRTAB shows you all the outputs and you can use the BOX input to select the columns you want. You can write to disc as an ascii file; you will then need your own programs to strip out the line continuation marks etc. and convert to more useful units (including converting chan No to velocity). This shows slightly edited output for channel 17 using BOX  1 2 12 3 13 4 14 5 15 23 24 26 27 29 30. DELTAX/Y are the X and Y offsets and D0/-/+ are the deconvolved component mean, upper and lower axes. 

 1   2         12      3           13       4               14           5          15                 23            24             26           27             29            30
      PEAK   ERR   I FLUX   ERR    DELTAX     ERR        DELTAY     ERR         D0 MAJ    D0 MIN      D- MAJ   D-MIN       D+MAJ     D+MIN
          JY/BEAM         JY              DEGREES                   DEGREES              DEG         DEG          DEG        DEG          DEG         DEG
 17 3.345   0.0041 4.308   0.0084 -4.023E-05 1.065E-09 -5.542E-05 9.107E-10 1.216E-06  4.788E-07 1.212E-06 4.710E-07 1.220E-06 4.864E-07
 17 0.6124 0.0041 0.6295 0.0073  -1.522E-05 4.927E-09 -3.561E-05 4.823E-09 3.811E-07 1.085E-07 3.251E-07 0.000E+00 4.303E-07 2.232E-07
 17 0.8274 0.0041 0.1023 0.0082  -4.382E-05 4.303E-08 -5.726E-05 3.574E-08 1.210E-06 5.045E-08 1.027E-06 0.000E+00 1.378E-06 5.346E-07
 17 0.6966 0.0040 0.1021 0.0091  -3.765E-05 4.888E-08 -5.716E-05 5.604E-08 1.948E-06 0.000E+00 1.747E-06 0.000E+00 2.140E-06 5.620E-07

To get an idea of the sizes, 1E-06 deg is 3.6 mas.

DORESID 1 writes a residual file for each plane which is useful to check how well components have been subtracted, for example in chan 17 the weakest component fitted was 70 mJy/beam, above a cut-off of 24 mJy/beam; the residual image has a range of -23 - 25 mJy/beam so in this case the fitting has worked well in identifying peaks. You can also make a residual cube (task 'MCUBE') and subtract it from the original (task 'COMB'), for example in order to make a SQASH which is not affected by noise in bright channels.

In channel 17 SAD recovered a total of 5.1 Jy; IMST gives a total flux of 5 Jy in the measured area of chan 17, probably slightly less because of negative sidelobes.

You can plot MF components like CC using CCNTR; the components for 5 channels are shown over a single 0.025 mJy/beam contour. Yellow, green, pink, blue and red correspond to channels 41-45. The more northerly lobe is two distinct series of components. The more southerly shows a systematic displacement of position with velocity for the peak around 14.3269 s, 57".540-57".545 but some of the other components may be sidelobes.

The size of individual components is the beamed size and tends to be smaller for brighter components (see e.g. Elitzur 1992). The brightest component in chan 17 (above) has an integrated flux density of 4.308(0.008) Jy and a (geometric) mean FWHM of 2.75(0.03) mas. One way to get an honest estimate of the brighness temperature is to take the upper limit of the component size and the lower limit of the integrated flux density S which gives a lower limit on the brightness temperature Tb:

Gaussian Component Area A = pi * [(FWHM+FWHMerr)_radians]2 /(4 ln2)

Tb = [(S-Serr) * 10-26 * wavelength2]/[2 * KBoltzmann * A]

The total extent of a series of components with a systematic position-velocity gradient gives a measure of the unbeamed extent of an emision region, which may be bounded by density, by the velocity gradient, or by other factors. The bright southerly component in channels around 43 has an extent (above our flux cut-off) of 3.5(0.1) mas. These measurements can thus be used to investigate both the AU-scale structure of star-forming regions, and the maser physics such as whether amplification is saturated or unsaturated.

You can plot the spectral profile of individual series of components in AIPS, such as for the bright southerly component 

AIPS 1:task 'ISPEC'
AIPS 1:INNAME 'G23.389';incl 'ICL001'
AIPS 1:blc 88 151 30;trc 95 158 
AIPS 1:opty 'FLUX'

but it can be more accurate to plot the fitted components as a function of channel, especially if their position moves noticably, since there is no simple way of choosing the optimum box size in ISPEC.

Finally, it is important to examine and sort the output of SAD carefully, e.g. after printing out selections from the MF table with PRTAB, and compare the output with KNTR maps, to ensure that

Most recidivist maser astronomers or research groups have their favourite software to do this, if anyone is interested this can be discussed at ERIS.

Unresolved masers

Some of the terminology to be used was introduced in the preceeding section on resolved masers.

This example uses MERLIN OH cubes, either as you have made or loaded from /home/user/datastore/MERLIN_CUBES  SPER_99_67.ICL001.1 etc.

Unresolved OH masers

The blc's etc. here are appropriate for the ready-made SPER_99_67.ICL001 etc. made from the original channels 31-480, with 512 x 40 mas pixels. If your data have different dimensions use appropriate limits.

The plot of the cube, use of TVMOVIE etc. shows that the masers appear to be confined to within the central 10 pixels or so in any direction. It might be worth checking carefully for more extended emission but here we assume that there are no more far-flung masers. The cube has been CLEANed out to within a few pixels of each edge. As it is a high-Dec source the sidelobes are fairly symmetric.

We write a noise- and dynamic range-finding procedure within appropriate regions: 

proc DETNOISOH (i)
blc 31 31 i; trc 200 200 i;imst;j=pixstd
blc 313 31 i; trc 472 200 i;imst
j = min (j, pixstd)
blc 31 313 i; trc 200 472 i;imst
j = min (j, pixstd)
blc 313 313 i; trc 472 472 i;imst
j = min (j, pixstd)
j = 4*j;
blc 246 246 i;trc 267 267 i;imst
j = max (j, -1*pix2val)
return;finish

As the MERLIN aperture is better filled at 1.6 GHz and S Per is high-Dec, we try a 4 sigma cutoff.

The image header shows that the peak is ~5.5 Jy/beam for a 150-mas beam. At this resolution the masers are probably unresolved; these are typical SAD settings: 

AIPS 1:task 'SAD'
AIPS 1:blc 246 246 i; trc 267 267 i
AIPS 1:inver i
AIPS 1:ngauss 100
AIPS 1:cparm 1 0.2 j
AIPS 1:docrt 0
AIPS 1:sort 'S'
AIPS 1:outver -1;doresid -1
AIPS 1:doall 1;dowid -1
AIPS 1:dparm j,j,1000,20,20,20,1000

Test on a few channels e.g. around the blue-shifted peak 

AIPS 1:dowait 1
AIPS 1:for i=291:310;detnoisoh(i);inver i;cparm(3)=j;dparm(1)=j;dparm(2)=j;go sad;end

You will see that it produces a series of components close in position with offsets <100 mas and also the odd faint, isolated component which is probably spurious.

This shows MFPRT output for separate series of red- and blue-shifted channels (note for any real analysis you should use PRTAB to get uncertainties also!). I have edited the output for easy reading and also rejected any component without a counterpart in an adjacent channel and pixel. 

         Peak int   X-offset   Y-offset   
 Channel   JY/BEAM        mas        mas  
red-shifted 
     127     0.062     35.236      2.804  
     128     0.043     10.298    -42.318  
     129     0.070     43.168    -48.443  
     130     0.072     49.351    -55.783  
     131     0.064     55.092    -71.439  

     133     0.071     88.784    -51.531  
     134     0.067     29.371    -48.264  
     135     0.074     78.143    -70.435  
     136     0.125    116.733    -16.603  
     137     0.079     65.007    -78.373  

     139     0.074     69.537    -24.348  
     139     0.043    275.133    212.889  
     140     0.156     33.994    -43.957  
     140     0.043    268.116    200.928  
blue-shifted
     296     0.538     41.365      -3.260 
     297     0.979     46.165      -2.706 
     298     1.653     43.662      -2.694 
     299     2.706     38.789      -0.459 
     300     4.343     37.238      -0.044 
     301     6.529     38.702       0.358 
     302     8.095     40.023      -0.011 
     303     6.811     40.015      -0.800 
     304     4.177     38.796       0.114 
     305     2.847     35.831      -0.018 
     305     0.050     39.448     202.498 
     306     2.130     36.956       1.117 
     306     0.085     45.212     183.251 
     307     2.360     36.301      -1.072 
     307     0.062     38.218     159.653 
     308     3.555     36.271      -2.189 
     309     5.466     36.624      -2.758 
     310     6.125     36.109      -3.485

The brighter components are much closer together which might be due to a number of causes e.g.:

However, the majority of the red-shifted components are clustered around (50,-50) and the blue-shifted ones around (40,0); a proper analysis using errors shows that the difference is significant and thus component fitting can resolve maser shells which have a total size similar to the beam size or less.

Fitting to polarized data cubes

LL and RR cubes can be treated the same as total intensity - however, the noise will be (root 2) higher - as can any cube in which you expect all positive emission, e.g. total linear polarization made from Q and U in COMB. However, Stokes Q U and V cubes, and also quantities such as polarization angle, can legitimately be positive or negative, but the cutoffs in SAD are all positive. JMFIT can be used to fit manually using a negative cutoff but even that does not allow both positive and negative components to be found at once.

Phil Diamond has written an AIPS task 'MFQUV' which can be used to measure the flux density at the position of components in an existing MF table. This is valid because there cannot be polarized flux in the absense of total intensity.

Simulations

You can add artificial components to data in either the image or the uv plane. Task 'IMMOD' is used for the former but the latter, using task 'UVMOD', is a more stringent test as you can investigate the effects of varying the weighting in IMAGR as well as the accuracy of fitting.

In either case, start with a single maser-free plane of the image cube or calibrated uv data (well clear of any edge effects). For example, if G23.389.SPLIT was made from channels 721-830 of the original EVN methanol data, then channel 1 can be extracted from G23.389.ICL001 using SUBIM or from G23.389.SPLIT using SPLIT again. In the latter case, let's call this G23.389MOD.SPLIT.

Use the parameters of fitted components as guidelines for components to add, experimenting with spacings etc. This is a typical experiment: 

AIPS 1:TASK 'UVMOD'
AIPS 1:INN 'G23.389MOD'; INCL 'SPLIT'
AIPS 1:OUTN 'G23.389M_1';OUTCL 'UVMOD'
AIPS 1:NGAUSS 4; CTYPE 1
AIPS 1:GMAX  4 0.1 0.5, 0.2 
AIPS 1:GPOS -0.05,-0.05, -0.06,-0.05,  0.02,0.02, 0.02,0.023
AIPS 1:GWIDTH 0.002 0.002 0, 0.005 0.005 0, 0.003 0.003 0, 0.003 0.003 0
AIPS 1:ZEROSP 0;flux 0; factor 1

This adds 4 circular Gaussian components typical for G23.389, one pair separated by 10 mas with a 40:1 peak contrast, and the other pair closer in flux density and only 3 mas apart.

IMAGR the output G23.389M_1.UVMOD, have a look at the map in the TV and use IMST to check that it is roughly what you expect and run SAD on the map.

Depending the nature of the source you want to compare with, you could:

You can use this method for various tests including:

  • 1 If the end channels or noisy edges are retained then the header peak may not be useful, in which case check the part of the cube with real signal.

GaussianComponentFitting (last edited 2007-09-27 09:45:48 by AntonSmit)