3C147 Calibration Tutorial

This is a step-by-step tutorial for calibrating WSRT data using Calico and MeqTrees. This assumes a certain level of familiarity with MeqTrees on the part of the reader -- you should at least be able to run the browser and meqserver, load scripts, look at bookmarks, and generally know what bit goes where. You should also know what the Measurement Equation looks like, and why. You needn't know anything about Calico per se though, as pretty much every nut and bolt (on the user side, at least) will be documented here. Please contact me by e-mail on <smirnov AT astron.nl> should you have any questions.

Note that as our calibration scripts are constantly evolving, some of the screenshots here may be out of date. Specifically, new GUI options may have been added. If you encounter any new options, it is usually safe to leave them at their default value.

For starters, I recommend making yourself a new directory to play in. In this directory, you will need to copy (or symlink) the following three scripts from Cattery/Calico: calico-wsrt.py, calico-view-ms.py and calico-flagger.py. To find the path to your Cattery installation, use the Help|About... option in the meqbrowser.

Now, download and unpack 3C147_spw0.MS.tgz. If you don't have aips++ installed, you will not be able to do automatic flagging on the MS, so you'd better download a pre-flagged MS (3C147_spw0_preflagged.MS.tgz) instead.

Also, download 3C147.MDL -- this is a NEWSTAR sky model for this field that we'll be using for calibration. (Building up your own sky model is a separate can of worms, best left to another tutorial.)

We'll start by looking at the DATA column (observed visibilities) of your MS. Load up the calico-view-ms.py script and set compile-time options as per the attached screenshot (the MS contains garbage XY and YX correlations, so we tell the script to read diagonal terms -- XX and YY -- only). Compile the script, then set run-time options as per the attached screenshot (an explanation of all the imaging options will be given later on.) Note especially the TaQL selection string we supply to the imager -- this causes the imager to use fixed-movable baselines only (WSRT fixed antennas are 0-9, movable are 10-13). This eliminates redundant short baselines, which produces lower sidelobes.

You can now click on "Make a dirty image" to make an image of the DATA column. The image should look something like the attached fits file. All you can really tell here is that we have a bright source, which is 3C147 itself. If you've got the patience to look up and compare coordinates, you will find that in this map the source position is mirrored w.r.t. the phase center. This is due to lwimager (and the aips++/casa imager it is derived from) using a UVW sign convention that is opposite to that of the WSRT. (For the same reason, we will later tell our calibration script to "Invert phases in input data" -- thus producing residuals with inverted phases, which will cause the imager to un-mirror the sky and put all sources at the correct positions.)

Now, open the "Inspect input visibilities" bookmark, and open a bookmark or two from the "Input visibilities by baseline" submenu. Click on "view MS" to run through the MS. Look at the browser displays. You will see a set of per-baseline amplitude tracks (frequency-averaged), and also some 2D time-frequency plots of selected baselines. You can use the right-click menu in the plotters to get at things like phase plots, amplitude-phase displays, etc. Also, the middle mouse button will display cross-sections through the 2D time-frequency plots.

Useful information you can glean from these plots:

* we have data!

* we need to flag -- the amplitude plots have obvious outliers, especially baselines with antenna C.

* the edges of the bandpass are not very well-behaved (you can establish this by looking at cross-sections), and so some channels at the edges of the band should not be used for calibration.

We're going to use channels 3-57 for calibration to avoid the band edges (NB: WSRT channels are numbered in order of DECREASING frequency, so channel 0 is on the RIGHT in our time-frequency plots), and apply Hanning tapering. Tapering will be done later on-the-fly without changing the MS. After tapering only every other channel is going to be used (so only the odd channels -- 3,5,7,...,57 in our case). However, an old bug in the imager prevents it from skipping channels properly, so we're going to have to flag the even channels to force the imager to ignore them.

Load the calico-flagger.py script and set compile-time options as per the screenshot. Compile the script. You will see runtime options as per the screenshot (runtime options 1). The first item is called "Initialize bitflag columns in this MS". MeqTrees can extend the MS with a BITFLAG column, thus allowing for 32 independent "flagsets" instead of the single boolean flagset supported by aips++/casa. Before you do any further flagging, click on the "Initialize" item to insert BITFLAGs into your MS.

Your runtime options will now have new options available. We're going to flag the even channels and store these flags into a flagset called "even channels". Setup your runtime options as per screenshot "runtime options 2", and click on "Run the flagger".

Now let's check the flags we've put in. Open up the "View MS data & flags" menu, and set it up as per screenshot "runtime options 3". Here we're asking it to ignore the standard FLAG column, and only show the flags we've put into flagset "even channels". Open up one of the bookmarks for "Inspect input visibilities by baseline", then click on "View MS". Your time-frequency plots should now look like the attached -- the black vertical bands correspond to the flagged channels.

NB: if you don't see the vertical bands, this may be due to a bug in the browser (if you do, then it has already been fixed). In this case restart the browser and meqserver, reload the calico-flagger.py script, and do "View MS" again.

We will now try to automatically flag bad data using the aips++ autoflag tool. calico-flagger.py provides a graphical interface to this tool.

Open up the runtime options menu, and setup the "Run autoflagger" submenu as per screenshot 1. This will apply a spectral rejection flagger (thus flagging "bad" rows which deviate from the median bandpass shape) a time-median flagger, and a uv-amplitude outlier flagfger. For details on these methods, see the "autoflag" section in the AIPS++ User Reference Manual, http://www.astron.nl/aips++/docs/user/SynthesisRef/node3.html. Note also that we've checked the "Reset existing flags" option. This will clear the current content of the FLAG column first, so as to start with a clean slate. (This doesn't affect the "even channels" flags we set in the previous step, since those have been stored into a flagset -- in the BITFLAG column -- and not in the FLAG column.)

Now click on "Run the autoflagger" and wait for things to finish.

Now open up the "Inspect input visibilities" bookmark, and a per-baseline bookmark (one that contains antenna C, preferrably). Open up "View MS data & flags" and setup options as per screenshot 2 (we now want to read the standard a.k.a. "legacy" FLAG column, since this is where autoflag puts its flags.) Run "View MS" again -- you should see visibility tracks something like the attached image. Note how this has changed w.r.t. the original plot -- a lot of the bad data has been flagged, and no longer shows up. Look also at the per-baseline plots to see individually flagged visibilities. You can also look at the flagreport.ps file generated by the autoflagger for some interesting flag-related plots.

Now we need to transfer these new flags into a flagset. The MeqTree approach is to treat the FLAG column as a "scratch" flag column for autoflag. Once you're happy with the flags in the FLAG column, you can save them as a named flagset. This is done via the "Transfer FLAG/FLAG_ROW column into a flagset" submenu. Set it up as per screenshot 3, then click on "Transfer FLAG/FLAG_ROW".

Finally, we need to take the flagsets we have generated ("even flags" and "autoflag 1"), combine them, and fill the standard FLAG column with the results, so that the imager can see both sets of flags. This is done via the "Change the FLAG/FLAG_ROW columns" submenu. Set it up as per screenshot 4, then click on "Fill FLAG/FLAG_ROW".

If you're confused as to why all this transferring of flags back and forth is necessary, consider the following:

* Having just one FLAG column is extremely inflexible. Imagine what happens if you run the autoflag tool and it flags too much -- there's no way to back out the changes apart from clearing all flags and starting over.

* Bitflags (=flagsets) eliminate this weakness. A flagset can be removed or changed without affecting the other flagsets.

* Unfortunately, the imager (and other aips++/casa tools) are only aware of the FLAG column. To make them "see" our flagsets, we have to copy flagsets over to the FLAG column.

We're now ready to start calibrating. Restart the browser (if you've been playing around a lot, it will be consuming a lot of memory -- unfortunately, the only way around it is periodic restarts.) Now load up the calico-wsrt.py script and observe the multitude of compile-time options.

Calibration is a delicate business at the best of times, all the more so if you're trying to reach 1,000,000:1 dynamic range. You really need to get all the settings right. The screenshot will provide a reference for the initial settings, but it's also important to understand what all these options mean. The rest of this entry is dedicated to a walk-through of the settings.

* "MS" is, naturally, your MS.

* "Antenna subset" allows you to specify a subset of antennas. If your MS contains missing or bad data for one or more antennas, you can specify a subset here. For example, "0:4,6:" means "antennas 0 through 4, and from 6 on", thus skipping antenna 5. (Of course, it is also possible to use calico-flagger.py to flag an antenna wholesale, but specifying a subset here instead will produce a somewhat smaller and faster tree.) "None" means use all antennas.

* "Correlations" should be set to "2x2, diagonal terms only". This MS has garbage XY/YX correlations, this option will cause the calibration tree to ignore off-diagonal (i.e. cross-correlation) terms. Other possibilities are "2x2" (for calibrating all four correlations), "2" (for MSs that contain only XX and YY), and "1" (for MSs that contain only a single correlation.)

* "Start Purr", if checked, will launch the Purr tool (which is what I used to create this tutorial). Purr is a great way to keep a detailed log of all your calibration activities, so you really want this option on.

* "What do we want to do" determines the kind of tree that will be built. These options can be enabled in any combination. For this calibration, we want all three enabled (later on we'll see another combo):

- "Calibrate" will include calibration options. We select "Calibrate on: visibilities" (the other calibration options being somewhat experimental.) Furthermore, we want to omit short baselines from this calibration -- there's no extended sources in our sky model, so the shorter baselines will not contribute anything useful.

- "Subtract sky model and generate residuals" will take our sky model, apply instrumental corruptions determined during the calibration step, and subtract this from the data, thus producing "corrupted residuals", which can be written to a different column of the MS.

- "Correct the data or residuals" will apply corrections to the residuals (or to the original data, if the "Subtract" option is not checked), thus producing a set of "corrected" visibilities, which can be written to a different column of the MS.

* "Measurement Equation options" contains some global ME-related options. The important one enables time & bandwidth smearing on predicted model sources. Since this is relatively expensive, and makes little difference for the fainter sources, we can tell Calico to restrict smearing to only the 50 brightest sources in the sky model.

* "Sky model" specifies what sky model to use. The one we want is the "MeowLSM" module. Within that module, a number of options are provided:

- "LSM file" is simply the model file. Select the 3C147.MDL, which you should've downloaded along with the MS.

- For "File format", select NEWSTAR.

- "Primary beam expression" specifies an approximate model for the primary beam. Note that this is NOT the actual PB model applied during calibration (the latter is provided by the "E Jones" module below). Because our sky model contains INTRINSIC source flux, the LSM module needs some way of mapping that to approximate APPARENT flux, for the purposes of sorting sources from brightest to faintest. The "PB expression" option gives it a Python expression for this conversion. We'll specify the "NEWSTAR" expression here, since our sky model comes from NEWSTAR.

- "Use subset of LSM sources" can be set to limit the number of sources included in the sky model. Building a model with too many sources in it will make your trees rapidly eat up CPU and RAM, while the fainter sources contribute practically nothing to calibration solutions. For this reason, we tell the LSM to use the 50 brightest (apparently brightest -- this is where the PB expression above comes in) sources only, by specifying a subset of "0:49".

- "Make solvable source parameters" should be enabled if you plan to solve for source parameters, which we don't.

- "Save LSM in native format" will convert our NEWSTAR-format file to a native LSM format file. We will not need this option.

- "Save LSM in text format" will convert our NEWSTAR-format file to a plain ASCII table. This allows you to modify your sky model with a simple text editor. We will not need this option.

- "Show LSM GUI" will cause a little GUI to pop up, which you can use to look at your sky model. We will not be needing this option.

* "Export sky model as kvis annotations" will cause the current sky model to be exported into a Karma-compatible annotations file. The kvis viewer can then display these on top of an image. Annotations will be fantastically useful when we start looking at residual maps.

- The "Filename" option specifies a file to write the annotations to.

- The minimum annotation is a cross at every source position. The "Label format" option can be used to include text labels next to the cross (check options help for details -- the label string selected here will include the source name and its apparent flux.) Don't use labels for large (100+) sky models, as it produces overcrowded annotations.

- "Export each N sources as a separate file" can be handy for very large sky models, as it will split the annotations into chunks (in order of decreasing apparent brightness). This lets you load them into kvis chunk-by-chunk. With only 50 sources in our model, we will not be needing this option.

Subsequent options deal with various Jones terms that you can plug into your measurement equation. Every Calico script can define its own set of available Jones terms. The WSRT script we're dealing with defines the following:

* "Use E Jones (primary beam)." This is the PB model used for calibration. If your sky model is defined in terms of apparent fluxes (which would be the case if you were building it up from residual maps), you do not need to include a PB. Our sky model is defined in terms of intrinsic fluxes, so we need to enable this option.

- "Beam model" specifies the kind of beam model. The only option (currently) is the WSRT cos^3 beam. (NB: this is the model for the voltage beam; the corresponding power beam model is cos^6.)

- The two "Apply pointing error" options only come into play if you're trying to solve for antenna-based pointing offsets. At time of writing (Dec 2008), this has never worked on real data for anyone, but the options are there for experimentation. We leave them off.

- "Beam scale factor" is a magical empirically-established number that determines the size of the beam. At 21cm it is known to be 65 or thereabouts. At higher or lower frequencies it can be on the order of 62~63. Ask your nearest friendly WSRT guru for details; we use 65 here.

- "Solving for beam scale" allows one to solve for the above-mentioned magic number. This is experimental, we're not going to go there.

- "Clip beam" is important because the cos^3 PB model gets inaccurate at the edge of the beam. A setting of .1 will clip the beam gain at .1 (thus clipping the power beam at .01). Normally, a gain of .1 would then be applied to all sources beyond the gain=.1 circle, but if "NEWSTAR-compatible beam clipping" is enabled, then clipping is only done where gain<|.1|. As the value of |cos^3| increases further out, this produces PB "sidelobes" which are at the same level as the main lobe. NEWSTAR does this on purpose to have at least a rough approximation of the first sidelobe. Since we're using a NEWSTAR-derived sky model (and so all sources have their intrinsic brightness determined via NEWSTAR's naughty clipping method), we need to have this option enabled.

* "Use dE Jones" includes a differential gain term in the ME. We will enable this option later on, but leave it off for now.

* "Use B Jones" includes a bandpass Jones term. We'll be doing per-channel selfcal, so this option stays off.

* "Use G Jones" includes a gain/phase Jones term. Enable this, and check the "FullRealImag" option. With this option, we're going to simultaneously solve for the real and imaginary parts of a 2x2 complex gain matrix. The other option ("DiagAmplPhase") is to solve for amplitudes and phases separately, but (in my experience) solving for the real and imaginary parts converges faster.

* Via "Use interferometer gains" and "biases", we can add multiplicative and additive interferometer-based errors (a.k.a. "closure errors") into our ME. We will deal with these later, for now leave both options off, since it makes the tree smaller.

Now press "compile" to build your calibration tree.

You should now see the Runtime Options dialog. Set it up as per attached screenshot.

The first option submenu determines how data will be read and written:

* "Input MS column" should be DATA, since that's where the observed visibilities reside.

* "Output MS column" is where the corrected residuals are going to go. CORRECTED_DATA is a natural choice. (Note that column names have no particular significance to MeqTrees, we could just as easily choose to write to the MODEL_DATA column, in fact later we will.)

* "Apply Hanning taper", if checked, will aplly Hanning tapering on-the-fly. We want this on.

* "Invert phases in input" should be on. WSRT MSs use an opposite UVW sign convention w.r.t. the imager, which produces mirrored images unless we flip the phases.

* "Data description ID" usually corresponds to the spectral window, in this case 0. If you run the tree and get messages about "MS selection returns no rows", it's ususally due to an incorrect DDID setting.

* "Channel selection" can be used to select a subset of the channels. As mentioned previously, we want to use channels 3-57 with a stepping of 2 (since Hanning tapering makes half the channels redundant).

* "Additional TaQL selection" can be used to further restrict the subset of the MS being operated on, by supplying a Table Query Language string (see http://www.astron.nl/aips++/docs/notes/199/199.html for details.)

* "Read flags from MS" determines which flags will be in effect during calibration. We definitely DO want to include the "autoflag 1" flagset. We DON'T want to include the "even channels" flagset, since Hanning tapering will flag a channel if any of its immediate neighbours are flagged (the taper having a width of 3), so including "even channels" and then doing tapering will, in effect, cause all channels to appear flagged! Likewise, we don't want to read the standard FLAG column, since that currently includes even channels.

The next submenu is called "Calibrate G diagonal terms". This is where we set up everything related to our selfcal gains solution:

* "Tile size" sets the size of a processing chunk. The MS is read timeslot by timeslot, and processed in chunks (tiles) of the given size. The effect of tile size is fourfold:

(a) tile size determines the maximum size of a solution interval (in time). In fact, if subtiling (see below) is not used, a single tile IS the solution interval.

(b) larger tiles are faster (to a point). It is not efficient to process the data in very small chunks. Once a tile is big enough though (~few hundred time/freq points), there's little to be gained in making it bigger. In our case, a tile size of 20 (28 channels, so 20x28=560 time/freq points) is already efficient enough.

(c) larger tiles consume more memory. If you make a tile too large (where the exact value of "too large" depends on the complexity of your tree), you WILL run out of memory. In our case, anything above 200 is unreasonable (on a 2Gb machine).

(d) when using subtiling (more on this below), larger tiles mean more discontinuity between solutions, thus making the solver iterate more.

Taking all of the above into account, 20 is a reasonable tile size for our G solutions.

* "Solve for G_diag" should be on: we're solving for the diagonal terms of the G-Jones term, which correspond to gains.

* The "Select solvability" options allow you to restrict your solution to specific subsets of parameters (by antenna, etc.) We will not be needing this.

* "Override initial value" allows you to specify a non-default initial value for parameters. We will not need this option.

* "Polynomial degree" allows us to fit polynomials (in frequency and time) over the solution domain. Since our solution domain is going to be 1x1 time/freq points, we leave the degrees at 0.

* "Solution subinterval" works together with tile size (above) to determine the solution domain. By setting 1 for both time and frequency, we're telling MeqTrees to do an independent gain solution for each time/freq point (a.k.a. per-channel selfcal). By judicious setting of tile sizes and subintervals (and polynomial degrees), you can give your parameters any time-frequency behaviour you desire. For example, 1 for freq and None for time will give an independent solution in each frequency channel, but constant in time (per each tile). This would be appropriate if we were solving for bandpasses. 1 for time and None for freq would do an independent solution per each timeslot, but common across all frequencies (i.e. if we were solving for freq-independent antenna gains). We will see more examples of this later on in our calibration.

* "MEP table name" specifies a table where the solutions are going to go. Tables reside within the MS. The default name is almost always appropriate.

* "Initialize with solution from MEP table" will cause each solvable parameter to start at whatever previous value was stored in the table (if any). This is useful if you're, e.g., refining previously-obtained solutions. You almost always want this option on, unless your solution table contains junk (in which case it's probably easier to remove the table in the first place.)

* "...even if time domains don't match" should be on if we're transferring calibrator solutions. Normally, solutions in the table are only loaded if their time/freq domains match. Since calibrator solutions will almost always have a different time domain, this option is required to use them.

* "Start from solution of previous time interval" should be on when solving for parameters that have temporal continuity (as almost all physical parameters do.)

* "Save solutions to MEP table even if not converged" will cause solutions to be saved even if the solver exceeds the maximum number of iterations (see below) without converging. Non-convergence is usually the sign of a poor model, but can often happen during difficult calibrations. Still, a poorly converged solution is usually better than none at all, so we leave this option on. We shouldn't run into non-convergence in the course of this exercise anyway.

* "Parameter constraints" are highly experimental, don't go there for now.

* "Clear out all previous solutions", if clicked, will remove solution tables (for this particular parameter -- G diagonal terms -- only) from disk. This is useful if you want to solve for something from scratch.

* "Simultaneously solve for other parameters" allows you to throw other parameters into the solution. We won't be needing this.

* "Solver options" contain a number of tweaks that affect solver behaviour. The defaults are usually good enough, but we do change three options here:

- "Convergence threshold" determines when the solution is considered to have converged. 1e-5 is a good value, but in low-SNR situations you may want to set this to 1e-6 or even lower.

- "Max iterations" tells the solver when to give up if things are not converging.

- "Subtiling convergence quota" tells the solver how many subdomains (within a tile) need to converge before it considers the whole solution to have converged. Our tile is 20 timeslots and 28 channels, our solution subdomain is 1 timeslot and 1 channel (thus giving 560 subdomains per tile). By setting the quota to 1.0, we tell the solver to keep iterating until all 560 subdomains have converged. With difficult calibrations, you may want to set a lower quota to speed things up, but 1 will do just fine in our case.

Finally, the "Calibrate G diagonal terms" button will set the ball rolling. Before you press it however, double-check your options, and open at least two bookmarks:

* "Inspect corrected data/residuals" is something you should ALWAYS have open, during any calibration. A quick look at the residuals will tell you practically everything you need to know about how your calibration went.

* "Inspector:G" will show your gain solutions.

Now, click on "Calibrate G diagonal terms", and go get a coffee or something, as this will take a while. As the calibration progresses, you can look at the status display at the bottom of the browser (see screenshot). This tells you:

* Approximately how long the calibration will go on

* Which tile it's currently on (in terms of tile number, and in terms of time relative to start of observation)

* Current solver iteration and chi-square value. Dropping chi-square = good, bouncing chi-square means the solver is having trouble converging. Iteration number turns greem if things converge, or red if things don't converge but the maximum iteration number has been reached.

* Solution rank and number of unknowns. 31360 unknowns corresponds to 14*2*2 unknowns (two reals per complex gain per 2 receptors per 14 antennas) per 560 subdomains. The rank tells you how many unknowns are actually constrained by the data. A slight drop in rank is normal for this calibration, as, barring closure errors, there's always 1 redundant phase term per 14 antennas. (In more difficult cases, dropping rank can be a symptom of a poorly constrained solution.)

* Number of subdomains per this tile (560), and how many of them have already converged to a solution.

* "Enough!" is actually a secret button. Clicking this button will tell the solver to stop iterating and go on to the next tile. The solution will be considered unconverged, as if the maximum number of iterations had been reached. When things go right, you don't need this button.

Now wait for the calibration to finish and study your residual plots.

Once things finish, you residual plot should look something like the attached. Apart from some crud on short baselines, the residuals look entirely reasonable. Since these shortest baselines are not used for our calibration or imaging, we don't need to worry about the crud (but we will shortly flag it anyway, just to keep it from cluttering up further residual plots.) Note also the G solution plots -- the phase of antenna C is drifting all over the place for some reason.

But first, let's make an image of the residuals. Open up the "Make an image for this MS" runtime submenu, and refer to the attached screenshot to set up options. Options are as follows:

* "Imager to use" is either "lwimager" or "aips++", depending on what yopu have installed. Either will work.

* "Image type or column" should be set to CORRECTED_DATA, since that's where we wrote our residuals.

* "Frequency channels in image" can be used to ask for per-channel images, or average of all channels, or multi-frequency synthesis. "Average all" will do in this case, since the band is rather narrow.

* "Imaging weights" should be set to "radial" for the WSRT.

* "Apply Gaussian taper" applies an additional taper to visibility weights. Our setting of 12"x12" (which corresponds to the synthesized beam size of the WSRT at 21cm) will emphasize point sources at the expense of extended structure. See http://www.astron.nl/aips++/docs/user/SynthesisRef/node142.html for details.

* "Stokes parameters" can be used to make polarization images. We only want "I" for now.

* "Image size in pixels": larger images take longer to make, but 2048x2048 is reasonably quick.

* "Image size in arcmin" sets the field of view. We'll go for 60'x60' for now, as the WSRT primary beam drops off beyond that size. You can also try 120'x120' -- there are a few sufficiently bright sources poking through that far out.

* "Image padding factor for FFTs" can eliminate aliasing in images. 1.2 is a good magic number (but try 2 or something if you see aliased sources.)

* "Enable w-projection" is not needed for WSRT, as it is a coplanar array.

* "Phase center" can be used to center the image on a specific location. The default is the phase center.

* "Custom MS selection" can be enabled to image a subset of the data (by default, the subset specified in the "Data selection" submenu is used.) We need to enable it so as to specify a channel stepping of 1 (remember the discussion of imager bugs, and why we have flagged all the even channels), and to supply a TaQL string: "ANTENNA1 in 0:10 and ANTENNA2 in 10:14". The latter is a standard WSRT trick that restricts imaging to fixed-movable baselines only, which cuts down on the grating rings caused by redundant baselines.

Now press "Make a dirty image". After a short while, your residual image should pop up, and it should look like the attached. If you're using Karma, try "Overlay:Load annotations" to load our LSM annotation file.

Note the following features of this map:

* Some of the central source (3C147 itself) is left behind, although well below 20,000:1 level already. This can be eliminated by solving for closure errors, which we will do shortly.

* The map is dominated by grating rings from off-center sources. Some of these exhibit irregular (i.e. not radially symmetric) patterns, which suggests that they're due to pointing errors and other image-plane effects, rather than errors in the sky model.

We're now going to do some quick flagging on the residuals. Load up calico-flagger.py again, compile it, then set up autoflag options as per attached screenshot. We're going to do spectral rejection again (with a 7-sigma threshold, as we only want to catch obvious outliers), and uv-binning. Note also that MS column is set to CORRECTED_DATA, and "Reset exisiting flags" is OFF -- we do want to retain previous flags (and we want the even channel flags to be taken into account, since there are no valid residuals on the even channels.)

Click on "Run the autoflagger" and wait for things to finish. Now let's look at the resulting flags. Open up "View MS data & flags", and set it up as per attached screenshot. We want to look at CORRECTED_DATA, with proper channel selection and stepping, and we only want to read the FLAG column. Open up the "Inspect input visibilities" bookmark, then click on "View MS data & flags".

Once the script finishes, you should see visibility tracks similar to the attached. Note how most outliers have been eliminated. We now want to again save this set of flags into a new flagset, which we'll call "autoflag residuals". Open up the "Transfer FLAG/FLAG_ROW column into a flagset" submenu, set it up as per screenshot, and click on the "Transfer" button.

Now is a good time to restart the browser, to clear any accumulated crud. Next, load up calico-wsrt.py again, and modify compile-time options as per screenshot:

* Enable "Use interferometer gains" (a.k.a. closure errors)

* Calibrate on all baselines now, without excluding the short ones (we could probably skip them, as they don't take part in imaging, but may as well solve for them for completeness.)

Now, compile the script and go into runtime options. Open up the "Calibrate IG (interferometer-based gains)" submenu and set the tile size to 120 (see attached screenshot). This means we're going to be solving for one gain per interferometer, per 120 timeslots (i.e. 1 hour), with no frequency variation. Interferometer gains are known to vary very slowly; also, any small frequency variation will be averaged out in our image, so it makes sense to keep the parameter count as low as possible.

Now open up the usual "Inspect corrected data/residuals" bookmark. If you're curious, also open a bookmark from the "IG (interferometer-based gains)" bookmark submenu, as this will show you the actual gain solutions on a number of baselines. Now click on the "Calibrate IG" button to start the solution. Get a short coffee, as this solution will be considerably faster.

Your residual plot should look something like the attached. Note how residual level has dropped w.r.t. to the previous plot. If you had bookmarks open for individual interferometer gain solutions, you'd see gains within ~1e-4 of unity (as you would expect from a decent interferometer).

Finally, use the runtime menu to make another residual image (see attached). Note that the central source has now completely disappeared, and the dominant remaining feature are grating rings from off-center sources. Congratulations, you have reached the NEWSTAR level of calibration. It is now time to go further and do something about those off-center sources.

Residuals from off-center sources can be eliminated by solving for differential gains in the direction of a handful of troublesome sources. But before we do this, let's subtract out all the really faint sources first. The rationale for doing this is as follows:

* Predicting many sources is computationally expensive (our model has 298...) It is for this reason that we've been calibrating using the 50 brightest sources only. We got away with this because sources beyond the first 50 are all fainter than 1 mJy, so selfcal is completely dominated by 3C147 itself (22 Jy).

* The residual errors we're looking at now are produced by a handful of off-center sources in the 5-30 mJy range. Once we start solving for such subtle effects, the contribution of the sub-mJy sources may become significant.

* On the other hand, differential gain (a.k.a. dE) solutions do not affect sources that do not have a dE term associated with them. So we may as well subtract all the sub-mJy sources in the uv-plane once and for all (using the current selfcal solution), and then continue calibrating what's left (which, in theory, will only contain contributions from the first 50 sources.)

This procedure will also nicely demonstrate the flexibilty of Calico scripts. Now then, restart the browser (for luck), load calico-wsrt.py again, and setup compile-time options as follows (refer to screenshot):

* "What do we want to do" should be set to subtract only. We want to take the observed visibilities, subtract the prediction for faint sources corrupted by our estimated instrumental errors, and write out the residuals. We do NOT want to correct the residuals at this stage, since we want them to be the "corrupted" observed visibilities minus the "corrupted" faint source model.

* "Time and bandwidth smearing" can be switched off, it's expensive and unnecessary for sources as faint as these.

* Under "Sky model", we set "subset of LSM sources" to "50:", meaning from source 50 on.

* Under "Export sky model as kvis annotations", we set a different filename (to have a different set of annotations for the faint sources), and set label format to None, so as to not overcrowd the annotations with labels.

Now, compile the script. Setup runtime options as follows (refer to screenshot):

* "Output column" should be set to MODEL_DATA. This isn't really "model data" we're producing, but remember, to us the column names are just labels (they have some meaning to aips++ itself, but -- fortunately -- we're not using aips++.) What we need is a column to store "observed data minus faint sources", and MODEL_DATA will do the job just fine.

* Under "Generate residuals", set the tile size to 20. No solving is going on here, so the tile size simply determines the size of a processing chunk. Since our tree is quite a bit more hefty than previous calibration trees (due to the large number of sources being predicted), anything over 20 will probably consume too much memory. 20 is fine on a 2Gb 32-bit machine, but if you're running on a 64-bit machine with RAM to spare, feel free to increase this size accordingly.

Now click on the "Generate residuals" button, and go get a medium-sized coffee.

By studying the last residual image, with LSM annotations overlaid, you can determine a number of "offensive" sources that fail to subtract cleanly. For starters, these are NEWS1073 NEWS1065 NEWS1000 NEWS1037 NEWS7 NEWS1003_1 NEWS1017. Let's now solve for differential gains in the direction of these sources.

Restart the browser for luck, and load up calico-wsrt.py. Setup compile-time options as follows (refer to screenshot):

* "What do we want to do" should be set to calibrate, subtract and correct once again. Baselines shorter than 144m should be ignored.

* "Time & bandwidth smearing" should be re-enabled.

* Under "Sky model", set the source subset to 0:49 again.

* "Export sky model as kvis annotations" should be on, as we will be modifying our annotations.

* We now enable "Use dE Jones (differential gains)". We also want to apply dE to only a subset of sources (applying dE to all sources is expensive and unrobust, mathematically speaking). Enter the list of sources given above, separated by spaces.

* We also want to annotate the selected sources using a different colour for the marker and label.

Now compile the script, then setup runtime options carefully (refer to screenshot "runtime options 1"):

* Input column is now MODEL_DATA, since that's where our "observed minus faint" visibilities reside.

* Output column should be CORRECTED_DATA as usual.

* Hanning tapering and phase inversion is now OFF. The previous subtract step has already applied Hanning tapering and inverted phases, so the visibilites in MODEL_DATA do NOT need these operations.

* Read flags: include flagsets "autoflag 1" and "autoflag residuals".

* Under "Calibrate dE diagonal terms", set the tile size to 60. Do not set any subintervals -- we'll be solving for just one differential gain per 60 timeslots and all frequencies (and per source, per antenna, of course.) This will be sufficient to drive off-center errors below the noise level, since our off-center sources are relatively faint (only x1000 brighter than noise.) Should you ever need to calibrate for differential gain towards a brighter source, you may need to give the gain term more time/freq variability.

* Under Solver options, set convergence threshold to 1e-6. Differential gains are a low-SNR parameter (again, because of the relative faintness of the sources), so they appreciate tighter convergence.

Now, load up bookmarks for "Inspect corrected data/residuals" and "inspector:dE". You can also open a "Solver" bookmark, to fully appreciate the minuteness of the effect being solved for here.

Now click on the "Calibrate dE diagonal terms". Coffee time, make it a large one.

At the end of the solution, your residuals and dE plots should look something like the attached images. Now go and make a residual map (one is attached here for reference.)

It helps to load two images in to Karma simulatenously -- the map we made now, and our previous-best residual map (the one made after calibrating for closure errors). Overlay your newly-made LSM annotations (these should clearly indicate what sources have differential gain solutions), and use the Histogram control of Karma to set intensity clipping to +/-0.15mJy (you will need to do this for both images, or else use "Freeze displayed intensity range" in the "View" menu). Now flip between the two images and marvel at the difference. If everything went well, it should be obvious that all sources with dE solutions have been comprehensively flyswatted!

A number of other "offensive" sources have now become visible, namely NEWS1009 NEWS3159 NEWS1001_2 NEWS1043. Let's solve for their differential gains as well. Restart the browser for luck. Reload calico-wsrt.py. In compile-time options, make just one change: under "Use dE Jones", go to the source subset option, and append "NEWS1009 NEWS3159 NEWS1001_2 NEWS1043" to the list of sources given there. Be careful to not erase the other sources, as this list needs to contain ALL sources that have an associated dE term.

Compile the script, then refer to the screenshot for runtime options. Go into "Calibrate dE diagonal terms", and enable "Select solvability by subgroup". What we want to do here is tell the system to only solve for the four new dE terms that we have added, without touching the previously solved for dE's. This we specify via the "solvability by subgroup" submenu.

Load up the usual bookmarks, and click on the "Calibrate dE" button. Time for a last coffee.

Your dE solutions ("data 2" tab) should look something like the attached. The residual plot should not be discernibly different from that of the previous step -- the effects we're solving here are truly tiny -- but is still worth looking at to make sure nothing's gone haywire. Of course, the real proof is in the residual image, so make it, and compare it to the previous two. Our four new "candidates" should now be eliminated.

Congratulations, you have beaten NEWSTAR! You can, in principle, keep adding more dE terms and eliminating more sources until you get bored, but this tutorial is now over.

Thanks for making it all the way to the end!