-
Notifications
You must be signed in to change notification settings - Fork 16
VLA Data Reduction
We will reduce VLA data from an 8 hour observation of the field around the 3C147 source at 1.267 MHz using the VLA's C configuration. The total bandwidth of the data is 256MHz, with a channel width of 4MHz and an integration time of 5s.
The data has been pre-flagged for RFI and the absolute flux scale and bandpass have been calibrated. In this reduction, we implement a MeqTrees based self calibration (self-cal) recipe that incoporates a VLA primary beam model (from Walter Brisken's cassbeam tool) as well a full treatment of direction dependent effects during both imaging and calibration.
We start with an initial calibration model which only consists of the 3C147 source; which lies at the phase tracking centre. This model was obtained from an observation of the same field with the Very Long Baseline Array (VLBA) telescope (R. Perley, priv. comm.).
The visibility data (MS format), primary beam model and the initial model can be found [here](place link here). It is import to organise the I/O flow of the pipeline as early as possible. In this case I create an input folder, "input", and an output folder "output" and a directory to store the visibility data (MS), "msdir". I then place the primary beam model (FITS files representing the Jones matrix of the primary beam), in a sub-directory, "fits", in the input folder and the initial model in the input folder, and place the MS in the folder "msdir".
Now, lets write the script. Create and open a Python file (mine is called stimela_vla-c_Lband_3c147.py). The first thing to do in as stimela script is to import stimela
import stimela
Then we organise the I/O flow of the pipeline in the script.
INPUT = "input"
output = "output"
MSDIR = "msdir"
This is the format that the primary beam model has to be specified for both MeqTrees (calibration) and DDFacet (imaging)
BEAM_FILES = "fits/JVLA-L-centred-$(xy)_$(reim).fits"
BEAM_L_AXIS = "X"
BEAM_M_AXIS= "Y"
INIT_MODEL = "3C147-PTVLACC.lsm.html"
Prefix for pipeline output products
PREFIX = "3C147_LBand_256MHz_VLAC"
This is a convenience function that creates filenames for the output images and sky models
def get_names(step, alpha=False):
prefix = "%s-s%d"%(PREFIX, step)
image = "%s.int.restored.fits"%prefix
model_tmp = "%s.lsm.html"%prefix
model = "%s.int.combined.lsm.html"%prefix
if alpha:
alpha = "%s.alpha.fits"%prefix
return prefix, image, model_tmp, model, alpha
else:
return prefix, image, model_tmp, model
Now, use the function above to set the filenames of the output products at various stages of the pipeline.
IMAGE1_PREFIX, IMAGE1, MODEL1_TMP, MODEL1 = get_names(1)
IMAGE2_PREFIX, IMAGE2, MODEL2_TMP, MODEL2, IMAGE2_SPI = get_names(2, alpha=True)
IMAGE3_PREFIX, IMAGE3, MODEL3_TMP, MODEL3 = get_names(3)
IMAGE4_PREFIX, IMAGE4, MODEL4_TMP, MODEL3 = get_names(4)
Solution/Smoothing intervals for the calibration steps
Gjones_INTERVALS = 1,1 # Solution intervals (time, freq) for G Jones solution
DDjones_SMOOTH = 60,24 # Smoothing intervals (time and freq) for dE Jones (differential gains) solution
These are the default options for the calibration steps options. I create this template so that I only have to update a handful of parameters for each calibration step.
calibrate_Template = {
"msname" : MS_List,
"ncpu" : 1,
"column" : "DATA",
"output" : "CORR_RES",
"Gjones" : True,
"Gjones_intervals" : Gjones_INTERVALS,
"Ejones" : True,
"beam_files_pattern" : BEAM_FILES,
"beam_l_axis" : BEAM_L_AXIS,
"beam_m_axis" : BEAM_M_AXIS,
"args" : ["stefcal_gain.delta=1e-6"],
}
A template of default options for the imaging steps.
NPIX = 8192
NFACETS = 21
CELL = 0.55
image_Template = {
"msname" : MS_List,
"PredictColName" : "MODEL_DATA",
"weight" : "briggs",
"robust" : 0,
"cellsize" : CELL,
"npix" : NPIX,
"nfacets" : NFACETS,
"add_beam" : True,
"beam_files_pattern" : BEAM_FILES,
"beam_l_axis" : BEAM_L_AXIS,
"beam_m_axis" : BEAM_M_AXIS,
"NCPU" : 20,
"NFreqBands" : 3,
"DegridBandMHz" : 4,
"SaveOnly" : "dPMRIieS",
"SaveCubes" : "dIiMX",
"PeakFactor" : 0.1,
"ResetDirty" : 1,
"ResetPSF" : 1,
"CycleFactor" : 0,
"MaxMajorIter" : 10,
"FluxThreshold" : 8.5e-5,
"clean_iterations" : 200000,
"SaveIntermediateDirtyImages" : 0,
}
def update(dict1, dict2):
tmp = dict(dict2)
tmp.update(dict1)
return tmp
Source finder options
REL = 0.9
THRESHOLD = 5
Lets start the pipeline
recipe = stimela.Recipe("VLA-C 3C147 reduction", ms_dir=MSDIR)
We first remove any flag sets from the MS, which may be there from previous reductions. We use the flagms tool from MeqTrees owlcat tool.
clear_flags = {
"msname" : MS_List,
"command" : " -r stefcal,stefcal"
}
recipe.add("cab/flagms", "clear_flags", clear_flags,
input=INPUT, output=OUTPUT,
label="CFL:: Clear flags")
cal_INIT_MODEL_G = {
"skymodel" : INIT_MODEL,
"reset" : "all", # Reset all existing gains
"DDjones" : False, # We not doing DD calibration at this stage
"Ejones" : True, # The is the jones term for the primary beam
"label" : "cal1",
}
recipe.add("cab/calibrator",
"subtract_INIT_MODEL_G",
update(cal_INIT_MODEL_G, cal_Template) ,
input=INPUT,
output=OUTPUT,
label="CAL1:: Subtracting initial model (G only)")
image_CAL_INIT_MODEL_G = {
"imageprefix" : IMAGE1_PREFIX,
}
recipe.add("cab/ddfacet", "image_CAL1_RES", update(image_CAL_INIT_MODEL_G, image),
input=INPUT, output=OUTPUT,
shared_memory="256gb",
label="IM1:: Image CAL1 CORR RES")
extract_sources_IM1 = {
"imagename" : IMAGE1,
"psfname" : IMAGE1_PREFIX + ".psf.fits",
"prefix" : IMAGE1_PREFIX,
"thresh_pix" : 8,
"thresh_isl" : 5,
"neg_thresh_pix": 4,
"neg_thresh_isl": 3,
"dd_tagging" : True,
"no_smooth" : True,
}
recipe.add("cab/sourcery", "exctract_sources_CAL1_RES", extract_source_IM1,
input=OUTPUT, output=INPUT,
label="SF:: Extract sources from CAL1 CORR_RES")
addspi = {
"image" : "${OUTPUT}/%s.cube.int.restored.fits"%IMAGE1_PREFIX,
"make_spi" : True,
"tol" : (-4,4),
"spi_image" : "${OUTPUT}/%s-specfit.fits"%IMAGE1_PREFIX,
"add_spi" : True,
"skymodel" : MODEL1,
"freq0" : FREQ0,
"sigma" : 20,
}
recipe.add("cab/specfit", "add_SPIs", addspi,
input=INPUT,
output=OUTPUT,
label="SPI::Add SPIs to MODEL1")
combine_models_1_and_2 = {
"operation" : "tigger-convert",
"command" : "--append ${INPUT}/%s ${INPUT}/%s --rename -f "%(REL, INIT_MODEL, MODEL1)
}
recipe.add("cab/tigger", "combine_models_and_2", combine_models_1_and_2,
input=INPUT, output=INPUT,
label="REL2:: Combine models one and two")
cal_MODEL1_G_dE = {
"skymodel" : MODEL1,
"reset" : "all",
"DDjones" : True,
"DDjones_smoothing" : DE_SMOOTH,
"IFRjones" : True,
"Ejones" : True,
"label" : "cal2",
}
recipe.add("cab/calibrator", "DDjones_subtract_MODEL1_G_dE", update(cal_MODEL1_dE, calibrate),
input=INPUT, output=OUTPUT,
label="CAL2:: Calibrate with MODEL1 (G+dE)")
image_CAL_MODEL1 = {
"imageprefix" : IMAGE2_PREFIX,
}
recipe.add("cab/ddfacet", "image_CAL_MODEL1", update(image_CAL_MODEL1, image),
input=INPUT, output=OUTPUT,
shared_memory=SHM,
label="IM2:: Image CAL2 CORR RES")
cal_MODEL2_G_DD_dE = {
"skymodel" : MODEL2,
"reset" : "all",
"DDjones" : True,
"DDjones_smoothing" : DE_SMOOTH,
"Ejones" : True,
"add_uvmodel" : True, # This add the deconvolution model (stored in the MODEL_DATA column of MS)
"label" : "cal3",
}
recipe.add("cab/calibrator", "DDjones_subtract_LSM_2", update(cal_MODEL2_G_DM_dE, calibrate),
input=INPUT, output=OUTPUT,
label="CAL3:: Final Cal Step (with dE). Subtract MODEL1 and Dec. Model")
image_CAL_MODEL1_DM = {
"prefix" : IMAGE3_PREFIX,
}
recipe.add("cab/ddfacet", "image_CAL_MODEL1_DM", update(image_CAL_MODEL1_DM, image),
input=INPUT, output=OUTPUT,
shared_memory=SHM,
label="IM3:: Image CAL3 CORR RES")