\documentclass[a4paper,12pt]{article}

%nicer layout
\usepackage[activate={true,nocompatibility},final,tracking=true,kerning=true,spacing=true,factor=1100,stretch=10,shrink=10]{microtype}

\newcommand{\pipe}{\mbox{\texttt{rPICARD}}}

\usepackage{listings}

\usepackage{setspace}
\onehalfspacing

\usepackage{amsmath}

\usepackage{longtable}

%better margins
\usepackage{microtype}

%for \uline - underline with line breaks
\usepackage{ulem}

%outer enumeration can start with letters
\usepackage{enumitem}

\usepackage{listings}
\lstset{breaklines=true, frame=single, basicstyle=\scriptsize}

\usepackage[top=55pt,bottom=55pt,left=64pt,right=62pt]{geometry}
\emergencystretch 6em%

\usepackage{xcolor}
\usepackage{hyperref}
\hypersetup{
    colorlinks,
    linkcolor={blue!90!black},
    urlcolor={green!60!black}
}

\usepackage{rotating} %for vertical tables
\usepackage{caption,tabularx,booktabs}
\newcolumntype{P}[1]{>{\centering\arraybackslash}p{#1}}

%skip items in enumerations
\makeatletter
\newcommand{\skipitems}[1]{%
  \addtocounter{\@enumctr}{#1}%
}
\makeatother

\title{Documentation for the Radboud PIpeline for the Calibration of high Angular Resolution Data
\href{https://bitbucket.org/M_Janssen/picard}{(\texttt{\pipe{}})}}
\author{Michael Janssen}
\date{\today}

\begin{document}

\maketitle

\tableofcontents

\newpage
\section{Preface}
\label{intro}

This documentation is a manual for the pipeline.
It is not meant be a guide that explains the methods and intricacies of radio astronomy.
If the reader needs more background information, I recommend two references:
`Interferometry and Synthesis in Radio Astronomy' by A. Richard Thompson, James M. Moran, and George W. Swenson Jr.
(the \href{http://www.springer.com/gp/book/9783319444291}{third edition} is available online for free)
and `Synthesis Imaging in Radio Astronomy II', edited by G. B. Taylor, C. L. Carilli, and R. A. Perley
(volume 180 of Astronomical Society of the Pacific Conference Series. 1999).
Additionally, the \href{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}
paper can be consulted for a formal documentation of the pipeline, explaining calibration choice made
in more detail and providing more information about the basic principles of VLBI data reduction.

This pipeline makes use of the \href{https://casa.nrao.edu/}{CASA} framework for radio astronomy data
processing. The purpose of the pipeline is to calibrate very long baseline interferometry (VLBI) data
from many different arrays. This became possible within CASA, thanks to a project initiated by 
\href{https://blackholecam.org/}{BlackHoleCam}. This project led to the addition
of new data reduction tasks to the CASA software package, which are crucial for the processing 
of VLBI data.\footnote{The most important new task is \texttt{fringefit()}, which is available since CASA 5.3.}
These tasks were developed by Ilse van Bemmel, Mark Kettenis, and Des Small at
\href{http://www.jive.eu/}{JIVE}.

\newpage
\section{The Basic Code Philosophy}
\label{codephilo}

\begin{enumerate}
\item No parameter is hard-coded. The user can fine-tune every knob to achieve an optimal calibration.
\item All parameters have sensible or self-tuning default values. It will suffice to set only a few
basic parameters (path to the data and the names of the observed science target and calibrators) to
achieve a good calibration completely hands-free.
\item Every step of the pipeline will produce diagnostic plots by default, which show if the
calibration worked as intended. The plots plots are saved in special directories -- it will be evident
when the pipeline was run and which set of input parameters were used.
\item Every step of the pipeline can be re-run separately and without any effort for optimization
(e.g., after trying different input parameters or after manually editing bad data).
\end{enumerate}

Following these points, experienced users can easily customize the pipeline to their needs, while
inexperienced users will still get good results when running the pipeline `blindly' with the default
settings.\footnote{If there are no severe issues with the dataset that need to be dealt with in a special way,
the pipeline should produce perfectly calibrated data when using the default settings -- thanks to self-tuning
parameters and the automatic flagging routines.}
Thanks to the many diagnostics of the pipeline, users can easily understand how and if the
calibration steps worked. If the diagnostics show that a certain step did not work properly,
then it can be re-run indefinite times without effort, allowing the user to optimize that step
and fine-tune the calibration procedure.

\pipe{} is a highly modular code, with a standardized wrapper for each calibration task.
The wrapper automatically takes care of potential smoothing and plotting of calibration solutions
as set by the user via a set of standardized input parameters,
on-the-fly calibration, and the passing of any required
complementary metadata from the Measurement Set (e.g., frequency, antenna, and scan information), in single objects.
This makes it straightforward to adjust any calibration strategy or to add new functionalities, see \autoref{strategy}.

\newpage
\section{Code Overview}
\label{codeoverview}

The code is organized into the main script \texttt{picard/main\_picard.py}, which handles the sequential calls to the
different CASA calibration tasks, and several modules in the \texttt{picard/pipe\_modules} folder:
For phase calibration, amplitude calibration, polarization calibration, diagnostics,
opacity corrections, automatic flagging, auxiliary functionalities, and
generic functionalities around the CASA calibration framework.
Additional scripts from JIVE for the amplitude calibration can be found in
\texttt{picard/pipe\_modules/JIVE\_scripts}.

\texttt{picard/picard} is a wrapper around \texttt{picard/main\_picard.py}, handling a proper call to \textit{mpicasa -c}.
\texttt{setup.py} must be used to link the pipeline to your CASA installation (saving the PATH to your CASA binary in
a \texttt{your\_casapath.txt} file). The setup script will also establish the mpi environment for your machine and
use the \texttt{input\_template} folder (\autoref{inpparams}) to write a first set of default input parameters
to \texttt{picard/input}. For the docker images, \texttt{setup.py} was already executed automatically.

By default, \texttt{setup.py} will search your whole file system for a CASA installation. 
Passing \textit{-d $<$path$>$} as command line argument will narrow down the search to the directory
\texttt{$<$path$>$} (useful if CASA is installed on a network disk or if the local file system is very big).

\texttt{picard/interactive\_utils.py} is a module that can be loaded into interactive CASA sessions
independently from the rest of the pipeline (if added to your PYTHONPATH). See \autoref{interactive}.

\texttt{picard/scripts/calibration} contains example scripts where \pipe{} is run to calibrate data and \texttt{picard/scripts/imaging} contains scripts where \pipe{} is run to image calibrated data.

\newpage
\section{MPI Scalability}
\label{mpi}

The pipeline makes use of the built-in \href{https://casa.nrao.edu/casadocs/@@search?SearchableText=mpi}{MPI}
scalability of CASA. Some tasks (like \texttt{applycal()} or \texttt{flagdata()} for example) are internally parallelized and will
run faster if more CPU power is available.

Fringe-fitting is the most time-consuming task. It is done on a per-scan basis; for speed-up each scan
will be fringe-fitted separately as a separate MPI job.

The \texttt{setup.py} script will prepare a default \texttt{picard/input/mpi\_host\_file}, filling it
with the maximum number of cores available on the host machine.
By default, the call to \textit{mpicasa} in \texttt{picard/picard} reads from this mpi host file via
\mbox{\textit{-hostfile \$\{pipedir\}input/mpi\_host\_file}}.
However, if a custom \texttt{mpi\_host\_file} exists in a local input folder (\autoref{inpparams}) that the code is told
to use via the \textit{-{}-input} or \textit{-p} command line arguments (\autoref{cmdline}), the computing resources in that host file
are used instead.

A \texttt{mpi\_host\_file} simply contains the number of cores that are to be used per machine. I.e., for a single machine,
a single line is given with the name of the machine and the number of cores (`slots') that CASA is allowed to use.
For an MPI-able computing cluster, multiple computing nodes/machines can be used by specifying the name of each
node/machine together with the number of available cores on a separate line.
Another option is to use the \textit{-n} command line option (\autoref{cmdline}) to specify the number of
cores that are to be used.

\subsection{Memory Safeguard}
\label{memsafe}

For datasets with many visibilities per scan (many baselines, long scans, large bandwidth), it can happen that 
the parallelized \texttt{applycal()} and \texttt{fringefit()} steps (\autoref{overview}) cause memory starvation if
the number of used CPU cores times the amount of memory needed to calibrate a single scan\footnote{Different
scans can contain different numbers of visibilities.} is larger than the system’s total amount of available memory.
If \pipe{} is run on a system with many cores but not a lot of RAM, the \texttt{mpi\_memory\_safety}
parameter in \texttt{picard/input/constants.inp} can be set to limit the number of scans that are
simultaneously calibrated by \texttt{applycal()} and to wait until enough memory is available before submitting
a \texttt{fringefit()} job to an MPI server.
The overhead of this memory monitoring and resource scheduling will decrease the performance a little bit.

\newpage
\section{Quick Start Guide}
\label{quickstart}

This section provides instructions on how to run the pipeline on the provided example EVN fits-idi file
in the \texttt{testing} folder. The purpose of doing this test is twofold:
\begin{itemize}
\item Test if the pipeline works as intended on your machine.
\item Follow the typical actions that need to be taken when reducing a dataset with the pipeline:
\end{itemize}

\begin{enumerate}
\item Pick a location on your file system where you want to work with your data (\texttt{testing/}).
\item Prepare a set of input parameters (\autoref{inpparams}), starting from the \texttt{picard/input} file
in the pipeline folder (assuming \textit{\$ ./picard/setup.py} was done already). 
It is recommended to have the input folder always together with the data itself.
So, copy the \texttt{picard/input/*.inp} files to a new input folder in your workdir (\texttt{testing/input/}). There, edit the
\texttt{observation.inp} file, setting the absolute path for the working directory.
For example, \texttt{workdir = $\sim$/Software/picard/testing} would be a typical input. Or use 
\texttt{workdir = \$pwd} to set the working directory to \mbox{\texttt{$<$/path/to/input/folder$>$/../}}
(which is specified with the \textit{-{}- input} command line argument, see below).
\item Your uv-data must be present in the working directory either as fits-idi files (or links to these files) or as measurement set. 
In \texttt{testing/}, a single example EVN fits-idi file is present.
\item Also, all useful metadata should be put in the workdir (or links to these files), see \autoref{filenames}.
An ANTAB table (\autoref{antab}) must be present for the amplitude calibration.
In \texttt{testing/}, there is an ANTAB table \texttt{example.antab}, a uvflag file \texttt{example.flag}, a source model file 
\texttt{3C84.smodel} (\autoref{smodel}), and a file with receiver temperature information \texttt{example.trx} (\autoref{opac}).\footnote{
No opacity correction is done for EVN data. The  Trx file is just a dummy example here.}
\item The fits-idi file in \texttt{testing/} is an unknown dataset, so we first want to know what is in there before setting any other input 
parameters. So, we run the pipeline in `inspection mode' (just loading the data and writing a listobs file, 
see \autoref{cmdline}): \textit{\$ picard -{}-input testing/input -l e}. For the test case you should see something similar to
listing 1 below.
\item Inspect the diagnostics folder (\autoref{diagnostics}) from the previous action, which was created in the workdir and read the
\texttt{list.obs} file that was created to get an overview of the observation: 
Pick the best reference antenna and note down science targets and calibrators.
\item Edit \texttt{testing/input/observation.inp}, defining science targets and calibrators and \texttt{testing/input/array.inp}, 
defining \texttt{array\_type} and \texttt{refant}. For the test case, set \texttt{array\_type = EVN} ;
\texttt{refant = ON} ; \texttt{science\_target = None} and all calibrators \texttt{= 3C84}.
\item Make a dry run (\autoref{quickmode}), while also making an initial flag version backup (\autoref{flags})
and without creating a new diagnostics folder this time (\autoref{cmdline}):
\textit{\$ picard -{}-input testing/input -q -r}. For the test case you should see something similar to the
listing 2 below.
\item For a real dataset, you could now perform the the complete calibration:
\textit{\$ picard -{}-input $<$path-to-input-files$>$}. If you know your data
beforehand and the pipeline has been verified to work on your system, you can skip actions 5,6, and 8.
\end{enumerate}

\lstinputlisting[caption=Typical terminal output for action 5 when using the example data in \texttt{testing/}]{example_terminal1.prt}
\lstinputlisting[caption=Typical terminal output for action 8 when using the example data in \texttt{testing/}]{example_terminal2.prt}

\subsection{Usage Examples}
\label{usage_examples}

This subsection gives some examples of things you can do with the pipeline.
Here, I assume that you have followed the above steps and that you are in a working directory, where the data, metadata, and an input folder
are present.
For an explanation of the command line arguments, see \autoref{cmdline}.

\begin{itemize}
\item Only load the data and create an initial flag version: \newline \textit{\$ picard -prq}
\item Make plots of uncalibrated and unflagged data, without creating a new diagnostics folder: \newline \textit{\$ picard -pdrq h,i,k}
\item Print overview of flagged data: \newline \textit{\$ picard -pq j}
\item Run flagging steps to see how much data they flag: \newline \textit{\$ picard -prq b,c,d,e,f,g,i}
\item Undo all flags, calibrate step 0, and make plots: \newline \textit{\$ picard -prq 0,h,i,k}
\item Full run of pipeline: \newline \textit{\$ picard -p}
\item Reload metadata (necessary after changing the source selection in 
the \texttt{input/observation.inp} input file for example) without creating a new diagnostics folder:
\newline \textit{\$ picard -pmdq}
\item Determine fringe solints again (assuming it is step 4 here, this can be different for different
array types set in the \texttt{input/array.inp} input file): \newline \textit{\$ picard -pfq 4}
\item Full run of pipeline excluding step 5 (if the calibration table still exists from an earlier run it will be applied nonetheless;
in that case it may be useful to first delete the calibration table): \newline \textit{\$ picard -prq x,0$\sim$4,6$\sim$99}
\item Interactive imaging of a VLBI.ms.avg dataset created by \pipe{} in the current working directory.
The source is 3C279: \newline
\textit{\$ mpicasa -n 3 $<$path-to/bin/casa$>$} \newline
\textit{import interactive\_utils} \newline
\textit{interactive\_utils.imager(`3C279')}
\item Interactive imaging of a EVN.uvfits file in the current working directory.
Perform an additional correction for the parallactic angle and allow the user to turn off interactive imaging
once a final set of CLEAN boxes are in place.
The source is 3C273: \newline
\textit{\$ mpicasa -n 3 $<$path-to/bin/casa$>$} \newline
\textit{import interactive\_utils} \newline
\textit{interactive\_utils.imager(`3C273', `EVN.uvfits', parang=True, goautomated=True)}
\end{itemize}

\newpage
\section{Overview of the Default Pipeline Steps}
\label{overview}

Below, a brief overview of the basic pipeline procedures is given. For a more detailed documentation, see \autoref{papers}. Actions which will potentially flag data are underlined.

\begin{itemize}

\item \textbf{Preparation procedures}. These are executed every time the pipeline is run (except if the \textit{-h} command line argument
(\autoref{cmdline}) is used.

\begin{itemize}

\item Start the MPI client.

\item Read all input parameters from \texttt{picard/input/*.inp} (unless overwritten)
      and store them in the \texttt{inp\_params} object, see \autoref{inpparams}.

\item Change directory to the workdir specified in \texttt{picard/input/observation.inp}.

\item Make a diagnostics folder with the current timestamp, see \autoref{diagnostics}.
Can be turned off with the \textit{-d} command line argument, see \autoref{cmdline}.

\item Write out the used input parameters in the diagnostics folder.

\item Grab an ANTAB table in the workdir folder
(searched for by file extensions, see \autoref{filenames}).
This step can be skipped for phase-only calibration, see the example in \autoref{quickmode}.

\item Make a VLA-type gain curve table from the gain block of the ANTAB table.
Skipped if no ANTAB table is present.

\item Get all fits-idi files that are to be loaded from all 
files with certain extensions (\autoref{filenames}) in the workdir folder. Use the tsys block in the ANTAB
table to attach a \texttt{SYSTEM\_TEMPERATURE} table to the first fits-idi file and \uline{load in all files 
to get a single measurement set}\footnote{Will flag visibilities which have a weight of zero.} if there is enough free space. Exits otherwise.
This step is skipped if a measurement set is already present.

\item Make an initial backup flagversion. Can fall back to that version if the pipeline is run again
(\autoref{flags}).

\item Save the output of \texttt{listobs()} to a file in the diagnostics folder. Exit here if
\textit{-l e} was passed as command line argument.

\item Determine a set of MS metadata
(either from scratch or by reading a previously determined one from disk).
This set contains info about the correlations, antennas, scans, and frequency setup.
Also, a few selected scans are determined for every source. These are as uniformly distributed
over the observations while scans with the most antennas present are picked. These scans are useful
for diagnostics.

\item Print a quick overview of all steps and exit here in the case of a dry run 
(\textit{-q} command line option without a list specified, see \autoref{cmdline}).

\end{itemize}

\item \textbf{Run steps part 1.} These are executed before the actual calibration tasks.
Quickmode (\textit{-q}) can be used to perform only a subset of these steps.
More information on the underlined flagging steps can be found in \autoref{flaggingalgo}.

\begin{enumerate}[label=(\alph*)]

\item Look for model files of the observed sources (searched for by file extensions, see \autoref{filenames})
and load them into the \texttt{MODEL\_DATA} column of the
measurement set as described in \autoref{smodel}.

\item \uline{Apply flag tables attached to the fits-idi files.}

\item \uline {Look for files with certain flag extensions (section \ref{filenames}) in the
workdir folder and apply the flags stored in these files 
(if necessary convert uvflag to casa flag format).}

\item \uline{Automatic flagging algorithm based on outliers in the autocorrelation vs time data}.

\item \uline{Automatic flagging algorithm based on outliers in the autocorrelation vs frequency data.}

\item \uline{Flag edge channels.}

\item \uline{Flag start and end segments of scans (quacking).}

\end{enumerate}

\item \textbf{Calibration tasks.} Perform incremental on-the-fly calibration 
      (for every task, the calibration solutions from all previous steps are applied on the fly where applicable),
      potentially smooth the solutions, and save plots of the calibration
      results to the diagnostics folder for every step.
      Quickmode (\textit{-q}) can be used to perform only a subset of these tasks.
      The steps below are examples for a polarization continuum VLBA
      dataset. For other observations the steps may be different (see \texttt{picard/main\_picard.py} and \texttt{picard/input/observation.inp}).
 
 \begin{enumerate}[label=(\arabic*)]
 
 \setcounter{enumi}{-1}

  \item Generate an \texttt{accor()} calibration table. Corrects amplitude errors introduced by faulty sampler 
        thresholds. The corrections are estimated by scaling the auto-correlations to unity. Not required for EVN data which
        has this correction already applied.
        
  \item Generate a scalar \texttt{bandpass()} calibration table. It uses the autocorrelations to calibrate the amplitude frequency 
        response of every station.
        
%  \item Generate an \texttt{acscl()} calibration table to scale the autocorrelations to unity on a per-integration basis, correcting for
%        potential errors introduced by the scalar bandpass. To do so, it only uses the inner $75 \%$ of the channels per spectral window.
%        This methods is based on (but not identical to) the new recommended VLBA amplitude calibration strategy devised by Craig Walker
%        (see VLBA scientific memo number 37).

  \item Generate a tsys calibration table.
    If necessary, correct for the opacity-induced source attenuation,
    using all available trec files for better fits (\autoref{filenames} and \autoref{opac}).

  \item Generate a gc calibration table. The system temperature and gain curve information are
        used to multiply the correlation coefficients in the raw data with system equivalent flux densities (SEFDs)
        to have the visibilities scaled to units of Jansky. See \autoref{antab}.

  \item Generate a single-band \texttt{fringefit()} calibration table to correct for instrumental effects.
        The instrumental phase and delay calibration is done by fringe-fitting every scan of the sources set via \texttt{calibrators\_instrphase}
        in \texttt{picard/input/observation.inp}
        (unless an explicit list of scans is given in \texttt{picard/input/array\_finetune.inp})
        and then interpolating across scans between all solutions which have a high enough SNR
        (the threshold is defined in \texttt{picard/input/array\_finetune.inp}).
        By default, scans will be fringe-fitter over their entire duration in this step.
        The basic principle of fringe-fitting
        is outlined in \autoref{fring}.

  \item Determine optimal \texttt{fringefit()} solution intervals for every scan on a calibrator source
    (either from scratch or by reading in previously determined solints from disk).
        This is done by looking at the distributions of SNR vs solution interval on every baseline
        to the reference station and taking the shortest solint within the coherence time that still yields to 
        detections on all baselines.

  \item Generate a multi-band \texttt{fringefit()} calibration table for the calibrator sources.
        Write out the same SNR-weighted average of RPC+LCP rates for both polarizations
        and smooth the multi-band delays in time (per scan and per antenna).
        Optimal solution intervals obtained from the previous step are used here.

  \item Generate a complex \texttt{bandpass()} calibration table. Corrects for the phase response
        over the bandpass for every antenna. (An amplitude correction is also possible, but not recommended, 
        because this is already done by the scalar bandpass. See \texttt{picard/input/array\_finetune.inp} for details.)

  \item Generate a RL delay calibration table using \texttt{gaincal()}. Calibrates the
        delay differences between the RCP and LCP receivers. Necessary for polarization experiments.

  \item Generate a RL phase calibration table using \texttt{polcal()}.
        Calibrates the RL phase for polarization experiments. The absolute RL phase
        establishes the EVPA of the source.

  \item Generate a D-terms calibration table using \texttt{polcal()}. Corrects for
        the instrumental leakage between the RCP and LCP feeds of every station.
        
  \item Generate a multi-band \texttt{fringefit()} calibration table for the science targets
        on scan-based solution interval with open search windows to take out the bulk delay and rate
        offsets of each scan.
        
  \item Determine optimal \texttt{fringefit()} solution intervals for all scans on science targets,
        with the same method as outlined for the calibrator sources above.
  
  \item Generate another multi-band \texttt{fringefit()} calibration table for the science targets.
        As opposed to the previous fringe-fit step, where the bulk delay and rate offsets are taken out,
        this step solves for residual intra-scan atmospheric phase distortions on the short optimal
        solution intervals determined by the previous step. Typically, narrow search windows are used here.
        Same post-processing as for the calibrator sources multi-band \texttt{fringefit()} step (see above).
  
  \end{enumerate}
  
\item \textbf{Run steps part 2.} These are executed after the calibration tasks.
Quickmode (\textit{-q}) can be used to perform only a subset of these steps.

\begin{enumerate}[label=(\alph*)]

\skipitems{7}

\item Clear all existing calibrated data from previous \texttt{applycal()} runs.

\item \uline{Run \texttt{applycal()}, using the tables generated by all calibration steps.
Flags all uncalibrated data (e.g., if no fringes were found in a scan).}

\item Write the amount of flagged data to a file in the diagnostics folder.

\item Save plots for a set of calibrated visibilities in the diagnostics folder.

\item Take the calibrated data and average it into a new MS.
Then, create uvfits files per source from the calibrated and averaged data.

\end{enumerate}

\end{itemize}

\newpage
\section{Input Parameters}
\label{inpparams}

By default, input parameters are read from all \texttt{*.inp} files in the \texttt{picard/input} folder. They are stored in a
single object \texttt{inp\_params}, which is passed around in the code. A small description is given for each
input parameter in the input files. One should definitely set the parameters in \texttt{observation.in}p each time the pipeline
is run for a new VLBI experiment. For different arrays, different frequencies, or different
calibration strategies, \texttt{array.inp} should be edited. Usually \texttt{constants.inp} can be left as the
default. For a-typical calibration strategies or severe issues in the dataset which require non-standard
calibration methods, the \texttt{array\_finetune.inp} parameters can be adjusted.
New input values can easily be added as
\texttt{name\_of\_parameter = value\_of\_parameter}.

Multiple values can be separated by a semicolon, so that \texttt{name\_of\_parameter} becomes an array within
the code. If not given as string, numbers are recognized as floats if they contain a decimal point
and as integers otherwise. Booleans and strings are trivially recognized. \texttt{value\_of\_parameter} can
then easily be accessed as \texttt{inp\_params.name\_of\_parameter} in the code.
For readability, lines can be continued with a backslash (\textbackslash) character before a new line.

With the \textit{-{}-input} or \textit{-p} command line arguments, a set of input parameters can be read from any directory, instead of the default
\texttt{picard/input} folder (\autoref{cmdline}).

The idea behind the input files is that users can keep and modify their own set of inputs.
Therefore, \texttt{picard/input} has been added to \texttt{.gitignore}. However, an up to date
developer's copy of the input files are kept as template in the repository's \texttt{input\_template} folder.
If these files are updated, then users will have to adjust their own input files accordingly,
as described in \autoref{updates}.

Additionally, \pipe{} determines a set of ancillary or `internal' metadata derived from
the MS data itself. This data is stored in a single object for quick access from the different
calibration modules. The stored information is about the polarization, frequency setup,
data integration time, stations, scans, and observed sources.

\subsection{Parameter Examples}
\label{paperinput}

Here, more detailed information is given about some of the input parameters
mentioned in \href{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}.

\subsubsection{Parameters set by setup.py}

The code uses the following prioritized reference station lists for the phase
calibration by default:
\begin{enumerate}
\item \textbf{EHT}: ALMA, LMT, APEX, SMT, IRAM30m, SPT
\item \textbf{GMVA}: ALMA, Effelsberg, IRAM30m, Los Alamos, Fort Davis, Pie Town, Kitt Peak
\item \textbf{VLBA}: Los Alamos, Fort Davis, Pie Town, Kitt Peak
\item \textbf{EVN}: Effelsberg, Yebes, Medicina, Noto
\end{enumerate}
The user could alter these lists based on the actual performance of the stations
in an observations. The most important effects that need to be taken into account
are weather and technical issues. These will become evident from the diagnostic output
of the pipeline; most importantly from plots of system temperatures.

The following default search ranges are used when looking for fringes to
calibrate atmospherically induced in-scan phase fluctuations:
\begin{enumerate}
\item \textbf{EHT}: 2\,s -- 30\,s
\item \textbf{GMVA}: 2\,s -- 180\,s
\item \textbf{VLBA}: 30\,s -- 400\,s
\item \textbf{EVN}: 60\,s -- 400\,s
\end{enumerate}
Depending on the overall weather conditions, exact observing frequency, and 
source brightness, the user could adjust these search ranges. The best indicator
is the number of fringe-non detections within the search windows, which are reported
in the diagnostic plots for the fringe-fit solution interval searches.
It is possible to specify different search ranges for calibrator sources and science targets.
In the mm regime, \pipe{} will try a few longer solution intervals to get detections
for stations in scans outside of the ranges given above. For the EHT, 60\,s and 120\,s
are tried. For the GMVA, 240\,s and 300\,s are tried.

For both the list of reference stations and phase calibration solution
intervals, new lists of solution intervals for other arrays can easily be added.

\subsubsection{Imaging Parameters}

\autoref{imager_params} gives a brief overview of the important parameters for \pipe{} imager.

\begin{longtable}{c c P{0.525\linewidth}}
\caption{Description of the most important \pipe{} imager parameters and how they are updated in each imaging plus self-calibration iteration.
For more information about the \textit{tclean} parameter, see \url{
https://casa.nrao.edu/casadocs-devel/stable/global-task-list/task_tclean/about}.}             
\label{imager_params}
\endfirsthead
\endhead
\hline\hline                           
Parameter Name & Default & Description\\ 
\hline  
   \textit{cellsize} & None & Size of an image pixel. If not set, CASA will determine it based on the longest baseline.\\
   \textit{imsize} & None & Size of the image in pixels. If not set, CASA will determine it based on the field of view.\\
   \textit{niter0} & 500 & Initial value for the maximum number of CLEAN iterations.\\
   \textit{cleaniterations} & `shallow' & Parameter to update the \textit{tclean niter} parameter, starting from the \textit{niter0} parameter.
   If set to `constant', niter will not be updated and for `shallow', \textit{niter} will be increased by \textit{niter0} after each iteration.\\
   \textit{gain} & 0.1 & The CLEAN gain (fraction of source flux subtracted frm residual image).\\
   \textit{multiscale} & [0,2,6] & Pixel sizes for multi-scale deconvolution.\\
   \textit{nterms} & 1 & Number of Taylor coefficients in the spectral model for multi-frequency deconvolution.\\
   \textit{robust} & 0.5 & Briggs weighting robustness parameter (Briggs, D. S. 1995, 27, 1444).\\
   \textit{threshold} & `auto' & Stop cleaning when the residual peak flux reaches this value.
   The default (`auto') uses the point source sensitivity determined from the data by the CASA \textit{imager} tool.\\
   \textit{nsigma} & 3.0 & CLEAN stopping criterion based on the median absolute deviation.\\ 
   \textit{startmod\_sc} & 0 & Do an initial phase self-calibration to a point source to align phases.
   The timescale can be set in seconds. The default (0) is to use the data integration time $t_\mathrm{int}$.\\
   \textit{timeavg} & False & Can average the data in time after \textit{startmod\_sc}.\\
   \textit{phase\_only\_selfcal} & [300s,10s] & A list of phase-only self-calibration steps to be done before the amplitude self-calibration
   and after \textit{startmod\_sc} and \textit{timeavg}.\\               
   \textit{startsolint} & 10 hours & The starting solution interval for the amplitude calibration.\\  
   \textit{solint\_denominator} & 2 & The factor by which the amplitude self-calibration timescale $T_\mathrm{amp}$ is lowered after each iteration.\\
   \textit{N\_sciter} & 0 & Maximum number of self-calibration iterations. The default (0) is to stop when $T_\mathrm{amp}$ reaches $t_\mathrm{int}$.
   If set to $-1$, only the phases will be self-calibrated. If set to $-2$, a single image without self-calibration will be made.\\
   \textit{minsnr\_sc\_phase} & 3.0 & The SNR cut for phase self-calibration solutions.\\
   \textit{minsnr\_sc\_amp} & 5.0 & The SNR cut for amplitude self-calibration solutions.\\
   \textit{flag\_last\_sc} & True & Flag data for which no self-calibration solutions are obtained in the last iteration.\\
   \textit{amp\_selfcal\_ants} & False & Can define a subset of station for which amplitude self-calibration is to be done.
   The default is to obtain solutions for all stations.\\
   \textit{station\_constraints} & None & Can constrain the range for amplitude self-calibration gains.\\
   \textit{station\_weights} & None & Can modify the data weighting of specific antennas.\\
   \textit{uvrange\_sc} & None & Can constrain the u-v range for which amplitude self-calibration solutions are obtained.\\
   \textit{uvzero\_mod} & None & Can anchor the large scale flux in the model visibilities to a specified zero-spacing source flux
                          density up to a maximum u-v range.\\
   \textit{mask} & None & Can supply files with CLEAN box masks for the source that are to be used for different iterations.\\
   \textit{interactive} & True & Build CLEAN boxes interactively in \textit{tclean}.\\
   \textit{goautomated} & False & Option to set \textit{interactive} to False at runtime for a final set of CLEAN boxes.\\
   \textit{usemask} & `auto-multithresh' & Enable the \textit{tclean} CLEAN auto-boxing capabilities (see text).\\
\hline
\end{longtable}

\newpage
\section{Command Line Arguments}
\label{cmdline}

The available command line arguments for \texttt{picard} are\\
\textit{[-p,-{}-input $<$new-input-folder$>$] [-f] [-m] [-l (e)] [-q ($<$list$>$)] [-r (a)] [-d] [-s] [-n $<$\#cores$>$] [-i] [-h]}.
They are all optional.
\begin{itemize}
\item With \textit{-{}-input}, the pipeline will look for *.inp files and a \texttt{mpi\_host\_file} in
\texttt{new-input-folder} instead of the default \texttt{picard/input/} folder (\autoref{inpparams}).
If the shorter \textit{-p} is used instead, the code will look for an \texttt{input/} or \texttt{input\_template/}
folder in the current working directory to use as \texttt{new-input-folder}.
With this, every dataset can have a custom set of input parameters (as it should be), while using the appropriate computational resources.
\item \textit{-f} overwrites the \texttt{inp\_params.fringe\_params\_load} parameter, forcing it 
to False: If there is an already stored set of optimal fringe-fitting parameters, then it will not be loaded
and instead be determined again from scratch. The fringe-fitting parameters contain the information
for each scan, which refant to pick and the optimal fringe-fit solution interval.
\item \textit{-m} overwrites the \texttt{inp\_params.ms\_metadata\_load} parameter, forcing it
to False: If there is an already stored set of metadata from the measurement set, then it will not be loaded
and instead be determined again from scratch.
Typically, this should be used when the science targets or calibrators in
\texttt{picard/input/observation.inp} or \texttt{num\_selected\_scans} in \texttt{picard/input/constants.inp}
are altered.
\item \textit{-l (e)} overwrites the default behaviour, where the \texttt{listobs()} task is not run again
if a listobs file already exists. With \textit{-l}, the pipeline will enforce to run \texttt{listobs()}
again instead of quickly copying over an old file. Normally, running \texttt{listobs()} takes some time 
and should not be done multiple times unless the MS changed internally (after running \texttt{split()} for example).
If not only \textit{-l} is passed as command line argument, but \textit{-l e}, the pipeline will exit after writing 
the listobs file (`inspection only'). 
This is useful if you do not know the names of the sources in the fits-idi files (\autoref{quickstart}).
\item \textit{-q}, the `quickmode', is explained in \autoref{quickmode}.
\item \textit{-r}, an extra option for the handing of flag versions for the quickmode, is explained in \autoref{flags}.
\item \textit{-d} overwrites the default behaviour, where all diagnostics (\autoref{diagnostics}) are stored in a new folder
with the current datetime attached to its name for this run. With \textit{-d}, the diagnostics folder that was last modified 
(typically, this should be the one from the previous run) will be used.
This is useful when you are running the pipeline step by step using \textit{-q}: combined with \textit{-d}, all diagnostics will still go 
in the same directory then.
If the same step is executed again in this mode, the diagnostics from the previous run will be overwritten.
\item If \textit{-s} is given, all calibration tables that were generated from previous runs and that would have been applied
for the \texttt{array\_type}\footnote{The \texttt{array\_type} is specified in \texttt{picard/input/array.inp} (\autoref{inpparams})
and defines which calibration steps are executed according to the \texttt{array\_specific\_steps} defined in 
\texttt{picard/main\_picard.py}.}
of this run are deleted before any new calibration tasks are executed (\autoref{overview}).
\item With \textit{-n}, the pipeline will not use the \texttt{mpi\_host\_file} file to establish the
MPI environment (\autoref{mpi}). Instead, \textit{$<$\#cores$>$} (must be some integer number $>$ 1) CPU cores will be used for this run.
\item \textit{-i} activates the interactive/supervised mode. The pipeline will operate at the user's pace,
waiting for a keypress for every step -- either to execute the step or to exit at this stage. 
This is useful if the user wants to examine the results from every step before advancing.
The module presented in \autoref{interactive} offers a suite of functions that can be used to
post-process calibration solutions from several pipeline steps 
in an interactive CASA shell before continuing to the next step.
\item \textit{-h} (or \textit{-{}-h}, \textit{-help}, \textit{-{}-help}): Print help message and exit.
\end{itemize}

The order of the command line arguments is irrelevant, for example:\\
\textit{\$ picard -m -q 1,3$\sim$5,b,e -r a -{}-input /data/VLBI/VLBA/M87/2008/input/ -f}\\
and\\
\textit{\$ picard -{}-input /data/VLBI/VLBA/M87/2008/input/ -f -m -r a -q e,1,b,3$\sim$5}\\
are equivalent. The only constraint is that the list after \textit{-q} must consist of comma separated characters without spaces
(\autoref{quickmode}).

Single-hyphen arguments can be combined as usual: \textit{-fmdq ...} will be expanded into \textit{-f -m -d -q ...} for example.

\newpage
\section{Quickmode}
\label{quickmode}

Quickmode (\textit{-q}) can be used to re-run certain steps of the pipeline.

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

Pass \textit{-q $<$list$>$} to \textit{picard} (see also \autoref{cmdline}).
For example:\\
\textit{\$ picard -q 1$\sim$4,6,b,c,e,d}

No spaces are allowed in \textit{$<$list$>$} and it must consist of comma separated values.
A $\sim$ can be used to indicate a range of numbers;
the example above is interpreted as 1,2,3,4,6,b,c,e,d.
If \textit{x} is given, it is interpreted as a,b,c,d,e,f,g,h,i,j,k,l,m,n... (for the lazy).

Information about which steps correspond to which number and letter are printed
for every run of the pipeline.
If \textit{-q} is used without a list specified: \textit{\$ picard -q}, \\
the pipeline will exit after the information about steps is printed
and not perform any calibration (dry run). Only the preparation procedures (\autoref{overview}) are executed.
A dry run is meant for testing and to quickly see which steps are executed in which order.

Rationale:\\
For multiple runs of the pipeline on the same dataset,
one can quickly redo certain steps after adjusting calibration parameters in \texttt{picard/input/array.inp}
or \texttt{picard/input/array\_finetunes.inp} (\autoref{inpparams} and \autoref{strategy}).
Or, a subset of steps can be executed, then the user can manually edit the calibration tables
(flag delay outliers in the single-band \texttt{fringefit()} table for example, see \autoref{interactive})
or flag visibilities and then run the remaining steps (see also the \textit{-i} command line option, \autoref{cmdline}).
Note that, if you are re-running single steps, all other previous existing calibration tables
will still be applied on-the-fly and they will still all be applied with \texttt{applycal()} as
described in \autoref{overview}. If this is not what you want, then you can either use the \textit{-s}
command line option (\autoref{cmdline}) or to move the old calibration tables to a different place.

Another pathway would be to run the pipeline in the standard mode,
then image the exported data and save source models to disk (with CASA \texttt{tclean()} for example, see \autoref{interactive}).
These can then be read in by a subsequent run of the pipeline (\autoref{smodel}).
Then, one can for example redo the fringe-fitting step with the model using \\
\textit{\$ picard -q $<$fringe-fit step number$>$ -r a}\footnote{See \autoref{flags} for the meaning of \textit{-r a}.}

Or, quickmode can be used to exclude the amplitude calibrations steps, when no ANTAB table (\autoref{antab})
is present\footnote{In that case it is necessary to set \texttt{pass\_missing\_antab = True} in \texttt{picard/input/constants.inp}.},
or to first do the phase calibration steps and the amplitude calibration steps at the end
(see also \autoref{strategy} on how to do that in a more systematic way).

Also, quickmode can be used to plot uncalibrated visibility data. In the default pipeline run (\autoref{overview}), one of the last steps
(after the calibration tables have been applied) is to plot visibility data. If you want to plot uncalibrated data, quickmode can be used
to run the \texttt{clearcal()} and plotting step without \texttt{applycal()}:
\newline \mbox{\textit{\$ picard -q $<$clearcal step number$>$,$<$plotting step number$>$ -r}}

\newpage
\section{Interactivity Capabilities}
\label{interactive}

\texttt{picard/interactive\_utils.py} contains functions that are independent from the rest of the pipeline.
The functions can be imported in an interactive CASA shell if the module is added to your PYTHONPATH.


The following functionalities are currently implemented:
\begin{itemize}
\item An interactive flagging GUI based on the CASA \texttt{plotcal()} task for a quick and easy 
post-processing of calibration solutions for:
\begin{itemize}
\item Single-band \texttt{fringefit()} calibration tables.
\item Multi-band \texttt{fringefit()} calibration tables.
\item \texttt{bandpass()} calibration tables.
\end{itemize}
\item An imager that uses CASA's multi-scale, multi-frequency \texttt{tclean()} image reconstruction algorithm
      together with self-calibration loops to make images of the observed sources, creating FITS files.
      It requires an interactive mpicasa session, which can be started with
      \newline
      \$ mpicasa -n $<$number\_of\_cores$>$ path\_to\_casa/casa [casa\_options]
\item A function to concatenate UVFITS files (useful when comparing CASA data with products from the EHT HOPS pipeline).
\end{itemize}

\newpage
\section{Fringe-fitting}
\label{fring}

Correlators employ sophisticated geometric models and make use of clock searches when they let the signals
from a pair of two antennas interfere -- i.e. cross multiplying their signals to form baseline-based
visibilities. However, the correlator models are never perfect -- in almost every case residual delays
(phase slopes versus frequency) and rates (phase slopes versus time)
will still be present in the data. These have to be taken out in post-processing, using a technique that
is called `fringe-fitting' (otherwise there will be large coherence losses when averaging in time and frequency).

The CASA \texttt{fringefit()} task will FFT the visibilities versus frequency and versus time and find 
the baseline-based residual delays and rates from the FFT peaks. These are globalized (going
from baseline quantities to station solutions) with a least squares approach and referenced to a common
reference antenna.

Typically, the first step is to align the phases across the separate frequency windows (IFs or spectral
windows) by fringe-fitting each window individually, referenced to the same frequency.\footnote{
At high frequency observations it may be necessary to take out the fast phase rotations beforehand. 
In CASA it is possible to do that with \texttt{fringefit()} using \texttt{combine=`spw'} as the
FFT will be done over the whole band, which is very sensitive.} The signal must
be strong in each window, so a bright source should be used.
Rates are zeroed and the phase plus delay solutions will be applied to all scans.
This procedure is called `instrumental phase and delay calibration'.
In the code, the method for the instrumental phase and delay calibration can be specified in
\texttt{picard/input/array\_finetune.inp} (\autoref{inpparams}), as explained in the comments for the different
parameters. One can either select specific scans or fringe-fit all scans on the instrumental phase 
calibrator specified in \texttt{picard/input/observation.inp}. Next, one can either fringe-fit over scan durations,
give a specific solution interval, or use the estimated optimal solution interval for each scan (see 
\autoref{ffopti} below). For the two latter cases, one should make use of the smoothing options in
\texttt{picard/input/array\_finetune.inp} to obtain single solutions per scan to correct for this instrumental
effect. It is recommended to integrate over scan durations.
Lastly, one can set which solutions are to be used. First an SNR cut is made (can be set
by the user, the default is 10). From the remaining solutions
one can either select only a single solution (with the highest SNR) per station and spectral window
that is to be applied to the whole dataset. Or interpolate across scans between all
solutions that made the initial SNR cut. This could be used if there are drifts in the electronics during
very long observing sessions or when recording equipment was restarted during a run for example.

Once the phases are aligned (coherent) across the whole frequency band, \texttt{fringefit()} can run over
the whole band for a higher signal to noise ratio -- so the weaker sources can be fringe-fitted as well. For high frequencies,
a source is not detected if you cannot find fringes on it. At lower frequencies it is possible to transfer
phases from a bright nearby calibrator, so that weaker sources can be detected (at lower frequencies
the troposphere does not distort the phases that much, so that solutions derived from another part of
the sky can be used). Phase-referencing with this pipeline is explained in \autoref{phref}.

\subsection{Picking a Reference Station}
\label{refantpick}

The \texttt{refant} and \texttt{refant\_minvaliddata} parameters in \texttt{picard/input/array.inp}
and \texttt{picard/input/array\_finetune.inp} respectively
determine which reference station will be used for which scan.
The \texttt{refant} parameter gives a prioritized list of antennas as options for the reference stations.
For each scan, the first antenna in the \texttt{refant} list which has a fraction of valid data
(number of unflagged visibilities over the number of flagged visibilities) larger than \texttt{refant\_minvaliddata}, unless less than $15\,\%$ of additional valid data was gained with respect to a station higher up in the \texttt{refant} priority list.
If none of the \texttt{refant} antennas makes this cut-off, the one antenna with the most valid data from that list is
picked as refant.

For each scan, fringes are referenced to the refant of that scan.
In the end, all fringe solutions will be re-referenced to a common refant with the CASA \textit{rerefant} task.
For polarization experiments, low values for \texttt{refant\_minvaliddata} should be used to avoid frequent re-referencing.
Also, if there is a single central and sensitive station is present in the array, this station should be first in the \texttt{refant} list
and \texttt{refant\_minvaliddata} could be set to a low number to ensure that this best station is picked most of the time.

Within scans, \pipe{} has the option to perform an exhaustive fringe search.
The correlated signal power is usually a strong function of baseline length for resolved sources --
the shorter the baseline length, the higher the SNR. For global 3\,mm observations for example,
picking Effelsberg as reference station may yield detections to all other European stations
and to the GBT across the Atlantic, but not to each individual VLBA station. However, fringes may
be present between the GBT and the VLBA. So both the European and the North American stations
can be calibrated by performing two fringe-fits -- one with Effelsberg as reference and one with
GBT as reference, where all fringe solutions to the GBT are re-referenced to Effelsberg using the
phase, delay, and rate relation between the two stations from the Effelsberg-GBT fringe detection in
the first fringe-fit run. A generalized version of this process is implemented in \pipe{}.
The algorithm will go through each scan $s$ in the observation and do the following:
\begin{enumerate}
\item Fringe-fit the data using the primary refant $P_s$. This refant is determined based on the
amount of valid data in each scan, see above. Store all fringes with insufficient SNR in a list $N$ = \{non-detections\}
and make a static copy $N_P$ of that list to keep track of the non-detections to the primary reference station.
\item Compute the list of secondary reference stations $S_s$:
If we denote all possible reference
stations, set by the \texttt{refant} parameter in \texttt{picard/input/array.inp}, as $R$, the secondary
refants for a particular scan are $S_s = \{R\} - \{P_s\}$.
\item If no non-detections are found in step 1, or when no detections are found between $P_s$ and any
station in $S_s$, the algorithm will go to step 6.
\item For each station $s_s$ in $S_s$:
\begin{enumerate}
\item Fringe-fit the data using $s_s$ as reference station and gather all stations, which can be connected
to $s_s$ via fringe detections, in a list called $f_s$.
\item Compute the connections to other reference stations as $c_1 = R \cap f_s$ and to stations not yet
connected as $c_2 = N \cap f_s$.
\item If
\begin{equation}
\left|c_1\right| \geq2
\end{equation}
or
\begin{equation}
\left|c_1\right|\geq1 \wedge \left|c_2\right| \geq1 \;,
\end{equation}
$s_s$ is added to a list $U$ of useful secondary reference stations as $U(s_s) \equiv [U_1(s_s), U_2(s_s)] = [c_1, c_2]$. The entries of
$U$ are referred to as `sub-clusters' in the code.
\item Update the list of non-detections: $N \rightarrow \{N\} - \{c_2\}$. Go to step 5 when $N$ is empty.
\end{enumerate}
\item Connect all sub-clusters:
\begin{enumerate}
\item The $c_1$ entries in each $U$ determine which $s_s$ can be connected.
All entries are discarded for which no path -- via any combination of multiple $c_1$ connections --
can be established to $P_s$. If that leaves $U$ empty, the algorithm will go to step 6.
\item For each station $n_p$ in $N_P$ -- if it is in any $c_2$ of $U$ -- the shortest path through $c_1$ connections, along the prioritized
list of reference stations, is stored as
\begin{equation}
E_s = \{ [n_p, P_s\leftarrow x_1 \leftarrow x_2 \leftarrow ... \leftarrow x_M] \; | \; n_p \in N_P \cap U_2 \;\, \mathrm{and} \;\, x_i \in U_1 \} \;.
\end{equation}
These $E_s$ paths determine the exhaustive fringe search for each scan $s$ employed by \pipe{}. In the notation above,
fringes for $n_p$ are obtained from a $x_M$ reference station, which can be connected to $P_s$ via a chain of other secondary
reference stations $x_{j \neq M}$.
\end{enumerate}
\item The FFT delay and rate solutions and their SNR of all fringe-fit steps in this exhaustive fringe search process
are stored as diagnostic output in csv ASCII format by the pipeline (\autoref{diagnostics} and \autoref{filenames}).
This information is usful to get a quick overview of all detections and non-detections in the dataset and to identify
gross residual post-correlation delay and rate errors.
\end{enumerate}
In the above process,
each fringe-fit is performed over full scan durations and only FFTs are computed to determined detections, the
least-squares is skipped.
SNR cuts for non-detections are based on $1.2\times$\texttt{fringe\_minSNR\_mb\_long\_sci} for science targets
and \texttt{fringe\_minSNR\_sb\_instrumental} for calibrator sources from \texttt{picard/input/array\_finetune.inp}
(\autoref{ffsteps}). The SNR simply set limits for scan-based SNR thresholds below which an exhaustive fringe
search needs to be employed -- no data will be flagged in this step. \pipe{} does the reference station search
and exhaustive fringe search before any other phase calibration steps. In all subsequent fringe-fit runs, all fringe
solutions are re-referenced according to the $E_s$ paths. That is, for every fringe-fit of a scan $s$, the phase, delay,
and rate solutions are re-referenced for every calibration table entry of $n_p$ in $E_s$.
Along the re-referencing chain, all flags are accumulated -- if a single fringe along the chain fails,
a flag will be written for the final re-referenced fringe as well.
It should be noted that fringe solutions have uncertainties $S$ due to the thermal noise in the system.
The total uncertainty along a re-referencing chain grows as
\begin{equation}
S_\mathrm{total} = \sqrt{\sum_{i=1}^M S_i^2} \;,
\end{equation}
where each $S_i$ represents the fringe solutions uncertainty of a specific station $x_i$ along the chain.
Considering current VLBI arrays, maybe one or two secondary reference stations $x_M$
are typically chosen for an exhaustive fringe search and they will have direct connections to $P_s$ in most cases.
The \texttt{fringe\_exhaustive\_refant\_search\_depth} parameter in \texttt{picard/input/array\_finetune.inp}
can be set to control the maximum depth of the exhaustive fringe search or to disable it completely.

In the future, a FFT baseline stacking functionality will be added to the \texttt{fringefit} task, which
will achieve similar results as the exhaustive fringe search.

\subsection{Finding the Optimal Solution Intervals}
\label{ffopti}

Optimal fringe-fit solution intervals should be long enough to have enough SNR for a detection and short enough to capture
atmospheric phase fluctuations. This is done first in the \texttt{task\_fringefit\_solint\_cal} step for the calibrator sources and
for the science targets in a later \texttt{task\_fringefit\_solint\_sci} step, where solutions derived from calibrator scans 
for instrumental effects like a phase bandpass and delay plus phase offsets between spectral windows can be applied on-the-fly.

The method used is to just run the FFT
(not the globalization) for many different solints for every scan and all baselines. First, a minimum detection SNR threshold
is set in \texttt{picard/input/array\_finetune.inp} (the default is 5.5). For all baselines that make it above this threshold, the smallest
solint that yields an SNR above the minimum is stored. From that collection the longest solint is picked -- this is the minimum solint
to have all possible detections on the different baselines. If \texttt{fft\_solint\_estimation=sqrt} is set in \texttt{picard/input/array\_finetune.inp},
the smallest solint that yields detections is taken as starting from which the expected coherent increase in SNR (by the square root
of the solution interval) is compared to the actual increase in SNR from the FFT for each baseline. Once the actual SNR drops below
the expected parabolic increase, the solint corresponding to that breaking point is taken as the coherence time on that baseline for that scan.
Finally, the smallest number from the collection of all coherence times on each baselines per scan is taken as the optimal solution interval.
Currently, a simpler and more robust method is implemented as default (\texttt{fft\_solint\_estimation=1.01}), where baselines are ignored if 
their SNRs do not increase with solint after a few iterations and the stopping criterion is set by SNR $>$ threshold.
The determined optimal scan-based solints will be used for all subsequent global fringe-fitting steps. There is the additional option
to fringe-fit certain scans again with different (typically longer) solints, if there are non-detections on certain baselines
-- see the \texttt{fringe\_solint\_mb\_reiterate} parameter in \texttt{picard/input/array.inp}.

\begin{figure}[t!]
   \centering
   \includegraphics[width=0.48\textwidth]{images/solint_estimation/detections_KP--BR.png}
   \includegraphics[width=0.48\textwidth]{images/solint_estimation/fringe_snr_vs_solint_KP--BR.png}
      \caption{Example of diagnostic plots from a fringe-fit solution interval optimization
               procedure. \textit{Left:} The number of detections within a specific scan for
               the KP--BR baseline for different solution intervals. It can be seen how the total
               number of detections goes down, as less solutions can be found within a fixed
               scan duration for longer segmentation times.
               \textit{Right:} The change in FFT SNR for different solution intervals. The minimum,
               median, and maximum SNR from all solution segments for each specific integration
               time are plotted for the polarization with the smaller overall SNR.
               The horizontal red line marks a detection threshold of $\mathrm{SNR} = 5.5$ and
               the vertical red line marks the minimum solution interval to have detections
               on all possible baselines, which was driven by a different baseline with a lower SNR.
              }
         \label{fig:solintsearch}
\end{figure}

All related data and plots from this optimization procedure are saved in the diagnostics (\autoref{diagnostics}) folder.
An example is shown in \autoref{fig:solintsearch}.

Dictionaries of the determined solution intervals and reference stations (\autoref{refantpick}) are stored on disk.
For re-runs of fringe-fitting steps, these parameters are loaded and used again, unless the \textit{-f} command line option is set (\autoref{cmdline}).

\subsection{The Individual Fringe-fitting steps}
\label{ffsteps}

Table \ref{tab:ff_params} gives an overview of the different
fringe-fitting steps implemented in \pipe{} together with
their main input parameter options, which are in
\texttt{picard/input/array.inp} and \texttt{picard/input/array\_finetune.inp}.

Science targets are first fringe-fitted on very long timescales to maximize the chance of getting a detection
with the \texttt{phase\_calibration.task\_fringefit\_multi\_sci\_long} step.
This is not done for the brighter calibrators to save computational time.
Each solution interval search following the method introduced in \autoref{ffopti} is done with the appropriate
FFT windows. For example, narrow fringe search windows are used for the calibration of intra-scan
atmospheric residual effects after the \texttt{phase\_calibration.task\_fringefit\_multi\_sci\_long} step for the science targets.
If the \texttt{fringe\_SNR\_cutoff\_float\_Perr} parameter in \texttt{picard/input/constants.inp} is set to True,
the min SNR fringe detection threshold is lowered based on a reduction of the probability of false detection
if narrow search windows are used.

\begin{sidewaystable}
\caption{Overview fringe-fitting steps and main input parameters in \texttt{picard/input/array.inp} and \texttt{picard/input/array\_finetune.inp}.}             
\label{tab:ff_params}      
\centering          
\begin{tabular}{c c c c P{0.45\linewidth}}
\hline\hline       
                      
Function\textsuperscript{\textit{a}} & search windows\textsuperscript{\textit{b}} & solint\textsuperscript{\textit{c}} & min SNR\textsuperscript{\textit{d}} & step description\\ 
\hline
   \texttt{solint\_cal} & \texttt{initial} & NA & \texttt{mb\_shortFFT} & Solint estimation following the \autoref{ffopti} method for calibrators.\\ \\
   \texttt{multi\_cal\_coher}\textsuperscript{\textit{e}} & \texttt{initial} & \texttt{mb\_coher} & \texttt{mb\_coher} & `Coherence' calibration
for the calibrators on \autoref{ffopti} timescales, where intra-scan phases and rates are solved before doing
the instrumental phase and delay calibration. This is done only at high frequencies.\\ \\
   \texttt{single} & \texttt{initial} & \texttt{sb\_instrumental} & \texttt{sb\_instrumental} & Instrumental phase and delay calibration.\\ \\
   \texttt{multi\_cal\_short}\textsuperscript{\textit{e}} & \texttt{initial} & \texttt{mb\_short} & \texttt{mb\_short\_cal} & Multi-band
fringe-fit on \autoref{ffopti} timescales, solving for phases, rates, and delays for calibrator sources.
No long integration like \texttt{multi\_solint\_sci} is done for the calibrators to save computational time --
the calibrator sources should be bright enough that this is not necessary.\\ \\
   \texttt{multi\_sci\_long} & \texttt{initial} & \texttt{mb\_long} & \texttt{mb\_long\_sci} & Multi-band
fringe-fit on long timescales for the science targets. This will determine if a source can be detected or not.\\ \\
   \texttt{multi\_solint\_sci} & \texttt{mb\_sci\_short} & NA & \texttt{mb\_shortFFT} & Solint estimation following the \autoref{ffopti} method for science targets.\\ \\
   \texttt{multi\_sci\_short}\textsuperscript{\textit{e}} & \texttt{mb\_sci\_short} & \texttt{mb\_short} & \texttt{mb\_short\_sci} & Residual multi-band fringe-fit
on \autoref{ffopti} timescales, solving for residual (after \texttt{multi\_sci\_long}) intra-scan phase, rate, and delay variations for science targets.
Typically, narrow search windows are used for the residuals.\\
\hline
\end{tabular}
\newline
\textsuperscript{\textit{a}} All function names have a \texttt{phase\_calibration.task\_fringefit\_} prefix.\\
\textsuperscript{\textit{b}} All FFT search window input parameters have a \texttt{fringe\_*\_window\_} prefix.\\
\textsuperscript{\textit{c}} All fringe-fit solution interval parameters have a \texttt{fringe\_solint\_} prefix.\\
\textsuperscript{\textit{d}} All fringe-fit SNR threshold parameters have a \texttt{fringe\_minSNR\_} prefix.\\
\textsuperscript{\textit{e}} If for some stations no fringes are found within the solint ranges for these steps,
                             these stations are fringe-fitted again on longer timescales, given by the
                             \texttt{fringe\_solint\_mb\_reiterate} parameter in \texttt{picard/input/array.inp}.
\end{sidewaystable}


\newpage
\section{ANTAB Tables: A-Priori Information for the Amplitude Calibration}
\label{antab}

You should get ANTAB tables from the observatories. They contain information about the DPFU, gain curve, and Tsys measurements
from every station during the run. In the list below, specific information about ANTAB tables for different arrays is provided.
You should put the ANTAB table (or a link to it) in your working directory (\autoref{filenames}).

\begin{enumerate}
\item \textbf{EHT}: You should have received an ANTAB table with you data. If not, then please contact \href{mailto:M.Janssen@astro.ru.nl}{me}.
If you want to calibrate the sidebands tgether, you may want to edit the INDEX rows of the ANTAB table.
\item \textbf{GMVA}: You should have received an ANTAB tables stations with you data. If not, then you
could contact \href{mailto:tkrichbaum@mpifr-bonn.mpg.de}{Thomas Krichbaum}. Otherwise, you may find your metadata
\href{http://www.vlba.nrao.edu/astro/VOBS/astronomy}{here}.
\item \textbf{VLBA}: \texttt{SYSTEM\_TEMPERATURE} tables should already be attached to your fits-idi files.
Also a \texttt{GAIN\_CURVE} extension should be attached as well.
If not, you will still need an ANTAB table with DPFU and gain curve information per station.
In that case, you should compile your own ANTAB table using entries from the \texttt{vlba\_gains.key} file.
Instructions on where to find that file on the public NRAO ftp server are outlined  
\href{https://science.nrao.edu/facilities/vlba/calibration-and-tools/caliblogs}{here}. 
You will need to select DPFU and POLY entries from every station for the right frequency range and observing time.
Moreover, you will have to adjust the format a little bit (switch $<$stationcode$>$ with GAIN, to get the format outlined below).
\item \textbf{EVN}: You should have received an ANTAB table with your data. If not, you will find one
on the \href{http://archive.jive.nl/scripts/listarch.php}{EVN archive} under your experiment code.
\end{enumerate}

A detailed description of the ANTAB format can be found \href{http://www.aips.nrao.edu/cgi-bin/ZXHLP2.PL?ANTAB}{here}.
In the simplest form, an ANTAB table looks like this:

1) First, the \textit{gain group} -- a single line for each station giving the DPFU and gain curve:
\begin{lstlisting}
GAIN <stationcode> <mount-type> DPFU = <rcp>, <lcp> POLY = <a0>, <a1>, <a2>, <a3>, ... /
\end{lstlisting}
Here, $<$mount-type$>$ is typically ALTAZ, $<$rcp$>$ and $<$lcp$>$ are the DPFU values for the RCP and LCP receiver,
and $<$$a_i$$>$ are the polynomial coefficients of the gain curve. \newline

2.1) Second, the \textit{Tsys group} -- blocks of system temperature with columns per IF (spw) and rows for the time spacing.
For each station a block starts with
\begin{lstlisting}
TSYS <stationcode> INDEX = `i1', `i2', `i3', ... /
\end{lstlisting}
Where, the `$i_j$' assign columns to RCP, LCP and IFs (spws). \newline

2.2) Next, the Tsys values for the station block are given; one line for every timestamp:
\begin{lstlisting}
<UT-day> <UT-time> <T1> <T2> <T3> ...
\end{lstlisting}
With $<$UT-day$>$ in the ddd format, $<$UT-time$>$ in a hh:mm.mm format, and $<$$T_j$$>$ 
the $T_\mathrm{sys}$ values in Kelvin assigned by $i_j$ to a receiver and IF (spw).\newline

2.3) Lastly, a station block is terminated with a / and the next one can start:
\begin{lstlisting}
/
TSYS <next stationcode> INDEX = ... /
...
\end{lstlisting}
\vspace{1cm}
The ANTAB table formatting requirements are more strict for CASA than for AIPS.
Every entry must be on a single line with keys and keyvalues written in the order shown above.
Example:
\begin{lstlisting}
GAIN AA ELEV DPFU = 1.0 POLY = 1.0 /
GAIN AP ELEV DPFU = 0.02478, 0.02521 POLY = 0.95515, 0.0022795, -2.8953E-05 /
GAIN AZ ALTAZ DPFU =  0.01683, 0.01681 POLY = 0.727089, 0.947364E-02, -0.822152E-04 /
GAIN PV ELEV DPFU = 0.127 POLY = 0.36029, 0.026008, -0.00026431 /
TSYS PV INDEX = 'R2', 'L2', 'R1', 'L1' /
141 03:00.317 131.89 142.49 132.70 142.70
141 03:01.867 132.07 142.76 132.77 142.76
!... and so on (exclamation mark = comment char)
/
TSYS AZ INDEX = 'L1', 'L2', 'R1', 'R2' /
141 03:01.456 122.24 291.24 117.59 332.75
141 03:02.956 122.18 292.14 117.66 331.82
!... and so on
/
!... idem for the other station blocks
\end{lstlisting}

The amplitude calibration scales visibilities formed from voltages in the correlator
to a physical scale in Jansky (Jy). 
The system equivalent flux density of a single telescope is given as 
\begin{equation}
\label{eq_sefd}
\mathrm{SEFD} = \frac{\mathrm{Tsys} \cdot \exp(\tau)}{\mathrm{DPFU} \cdot \mathrm{gc}} \; .
\end{equation}

Where Tsys is the system noise temperature in Kelvin (function of time), $\tau$ is the opacity (\autoref{opac}), DPFU
is the telescope gain in Kelvin per Jansky (the DPFU depends on the effective geometric area of the dish; it is mostly constant -- differences may occur
between day and night time), and gc is the normalized gain curve, which describes the variation of the antenna gain
with elevation.

If a calibrated correlation coefficient (visibility amplitude) on a baseline between stations 1
and 2 is given by $r_{1,2}$, then the correlated flux density $S_{1,2}$ in Jy will be
\begin{equation}
S_{1,2} = r_{1,2} \sqrt{\mathrm{SEFD}_1 \cdot \mathrm{SEFD}_2} \; ,
\end{equation}
where $\mathrm{SEFD_i}$ is the system equivalent flux density of station $i$ given by \autoref{eq_sefd}.

\newpage
\section{Diagnostics}
\label{diagnostics}

By default, the code produces logs and many diagnostic plots from the different calibration steps,
unless these options are turned off in \texttt{picard/input/array\_finetune.inp} (\autoref{inpparams}). The diagnostics are stored in a folder in
\texttt{inp\_params.workdir} that also has the current datetime (UTC) attached to its name.\footnote{Unless \textit{-d}
is used, see \autoref{cmdline}.}
This and the fact
that also the input parameters are logged, make it easy to run the pipeline several times on the
same dataset with different input parameters and to compare the results.

For every run of the pipeline, \pipe{} writes the usual \texttt{casa.log} file with the datetime attached.
Moreover, (because of the talkative nature of mpicasa), stderr is redirected to a 
\texttt{mpi\_and\_err.out} file with the current datetime attached as well.
If the pipeline finishes properly, both \texttt{casa.log} and \texttt{mpi\_and\_err.out} are moved over to the 
diagnostics folder of the run. If for some reason the pipeline crashes
(not terminating with a \mbox{\textit{- FINISHED -}} print to the terminal), both files are kept in the directory
where \pipe{} was called from.
They can be examined to find the reason for the crash.

Moreover, the command line call with the used command line arguments (\autoref{cmdline}) are saved to a file
in the diagnostics folder (as specified in \texttt{picard/input/constants.inp}, see \autoref{inpparams}) for every run. If
the \textit{-d} command line argument is used (writing incrementally to the same diagnostics folder), all
incremental calls are appended to the same file.

\subsection{Jplotter}
\label{jplotter}

By default, the code will produce plots of (calibrated) visibilities as one of the last steps (\autoref{overview}).
The recommendation is to use Harro Verkouter's jplotter program (See \texttt{README.md}) as it is faster and produces plots
with a better layout than the alternative standard CASA \texttt{plotms()} tool. Moreover, jplotter does not require
X-forwarding when saving plots directly to disk (it is unclear why this seems to be the case for \texttt{plotms()}...)
and it has a simple and intuitive syntax exemplified in
the \texttt{make\_jplotter\_plot()} function in \texttt{picard/pipe\_modules/diagnostics.py}.

The jplotter figures made by \pipe{} will show amplitudes and phases as a function of time (averaged over frequency)
and as a function of frequency (averaged over scan durations) for each baseline.
Different correlations (RR, LL and if specified also RL, LR) will be depicted with different colors. Unflagged
datapoints are shown as dots. Fully flagged datapoints are plotted as crosses.

\newpage
\section{Special Filenames}
\label{filenames}

For the default input parameters, the pipeline will recursively search the workdir specified in
\texttt{picard/input/observation.inp} (\autoref{inpparams}) for files with certain extensions.\footnote{Links will also be followed.}
These extensions are defined in \texttt{picard/input/constants.inp}:

\begin{enumerate}

\item \texttt{flagfile\_extensions}: All files with these extensions are used as metadata flag files.

\item \texttt{antab\_extensions}: A file which has one of these extensions is used as ANTAB metadata, see \autoref{antab}.

\item \texttt{fitsidi\_extensions}: All files with these extensions are used as fits-idi input files.

\item \texttt{modelfile\_extensions}: All files with these extensions are used as models for the observed sources, see \autoref{smodel}.

\item \texttt{trecfile\_extensions}: All files with these extensions are used to look for $T_\mathrm{rx}$
(receiver temperature) info per station, see \autoref{opac}.

\item \texttt{weatherfile\_extensions}: A file which has one of these extensions is used as weather metadata, needed for \autoref{opac}.
Such an external file should only be necessary for data from the GMVA.

\item \texttt{antenna\_mount\_corrections\_file}: If this file exists, it will be used to overwrite antenna's mount types in the 
MS. This is only necessary if erroneous mount types are written in the fits-idi input files.

\end{enumerate}

The other special filenames used by the pipeline in the working directory are:

\begin{enumerate}

\item \texttt{ms\_name} in \texttt{picard/input/observation.inp}: Name of the measurement set.

\item \texttt{gc\_dpfu\_fromidi\_file} in \texttt{picard/input/constants.inp}: Name of the file that will be created to convert DPFUs and gain curves from a GAIN\_CURVE FITS-IDI extension into an ANTAB ASCII file.

\item \texttt{gc\_conversion\_file} in \texttt{picard/input/constants.inp}: Name of the file that will be created to make a VLA-type gain curve from the ANTAB data.

\item \texttt{store\_optimal\_fringe\_params} in \texttt{picard/input/constants.inp}: Files (with additional .sci and .cal extensions) that will be created to store fringe-fit parameters.

\item \texttt{store\_scan\_refants} in \texttt{picard/input/constants.inp}: Files (with additional .sci and .cal extensions) that will be created to store dictionaries with the chosen reference station for each scan.

\item \texttt{diag\_fringe\_overview} in \texttt{picard/input/constants.inp}: Files (with additional .sci and .cal extensions) that will show an overview of FFT-fringes for all scans across all baselines to all reference stations used (including secondary ones, \autoref{refantpick}).

\item \texttt{store\_ms\_metadata} in \texttt{picard/input/observation.inp}: Name of the file that will be created to store internal (ancillary) metadata.

\item \texttt{diagdir} in \texttt{picard/input/observation.inp}: Name of the diagnostics folder that will be created with the UTC time attached to its name.

\item \texttt{diag\_*} in \texttt{picard/input/observation.inp}: Name of several files that will be created in the diagnostics folder.

\item \texttt{calib\_*} in \texttt{picard/input/array\_finetune.inp}: These parameters contain names for the calibration tables that will be created. For phase-referencing, tables with a .phaseref extension will also be created.

\end{enumerate}

\newpage
\section{Phase-Referencing}
\label{phref}

Phase-referencing is enabled if \texttt{calibrators\_phaseref} in \texttt{picard/input/observation.inp} (\autoref{inpparams}) is set.
For every \texttt{science\_target},
there must be a corresponding \texttt{calirbators\_phaseref} source specified.
An extra (.phaseref) table from the phase-referencing calibrators is created where all flagged(failed) solutions are removed and
all solutions are smoothed my a median filter on a per-scan, per-antenna, per-spw basis.
The fringe solutions from this phaseref table will be transferred to the science targets.
If \texttt{phaseref\_ff\_science = True} in \texttt{picard/input/observation.inp}, a fringe-fit is also done on the science targets 
after the solutions from the phase-referencing sources have been applied. This can be done if the science targets are strong enough to be
fringe-fitted for residual phases/rates/delays. Else, set \texttt{phaseref\_ff\_science = False}.

\newpage
\section{Opacity Correction}
\label{opac}

At high observing frequencies ($ > 15$ GHz), it is necessary to correct for the attenuation of the
source caused by the atmospheric opacity $\tau$.
This correction can be applied by raising the system temperature: $T_\mathrm{sys} \rightarrow T_\mathrm{sys} \exp(\tau)$.
If the system temperature is measured with a chopper,
then this correction is already applied and no further post-processing is necessary.
This is the case for data from the EHT array for example.
For other arrays, like the VLBA and GMVA, it is necessary to perform this correction.

The pipeline can solve for the opacity by solving the (approximate) $T_\mathrm{sys}$ equation for $\tau$:
\begin{equation}
\label{taufromtsys}
T_\mathrm{sys} = T_\mathrm{rx} + (1-e^{-\tau})T_\mathrm{atm} \Leftrightarrow \tau = -\log\left(1- \frac{T_\mathrm{sys}-T_\mathrm{rx}}{T_\mathrm{atm}}\right).
\end{equation}

Here, $T_\mathrm{atm}$ is the mean atmospheric temperature and $T_\mathrm{rx}$ is the receiver temperature.
The equation should be evaluated separately for RCP and LCP receivers.

For \autoref{taufromtsys}, $T_\mathrm{atm}$ must be estimated from the available weather (WX) data in the
\texttt{WEATHER} subtable of the measurement set. Users can easily add their favourite weather models in the 
\texttt{opacity\_correction.py} module of the pipeline.
Currently available are the accurate and recommended 
\href{http://cab.inta-csic.es/users/jrpardo/class_atm.html}{Pardo et al. model} 
and the simple weather model from 
Altshuler, Falcone, and Wulfsberg (1968): `Atmospheric effects on propagation at millimeter wavelengths'\footnote{
Here, $T_\mathrm{atm}$ is estimated from the ambient temperature $T_\mathrm{amb}$, which is
measured by the local telescope weather station, as: $T_\mathrm{atm} = 1.12 \, T_\mathrm{amb} - 50$ K.}.

The other unknown is the receiver temperature $T_\mathrm{rx}$.
By default, the code will estimate $T_\mathrm{rx}$ from an extrapolation of a $T_\mathrm{sys}$ vs airmass fit to zero
airmass.
The \texttt{opac\_corr\_airm\_max} parameter in \texttt{picard/input/array\_finetune.inp} can be set to exclude
low elevation (high airmass) data from the fit.
The fit may not be reliable if the atmospheric conditions change significantly over the
course of the observations.
By default, plots of the fits are made which can be inspected.
Also, the pipeline will print warnings if there seems to be problems with the $T_\mathrm{rx}$ fits.

There are two possible failure modes, both of which are typically caused by faulty estimates of $T_\mathrm{rx}$:

1) $T_\mathrm{rx} < \min{(T_\mathrm{sys})}$.\footnote{$T_\mathrm{sys}$ below 2.73 K are regarded as invalid and will be ignored.}
This can occur when the weather was bad during high elevation observations.
If this happens, the code will try to guess a new value for $T_\mathrm{rx}$,
based on the $\min{(T_\mathrm{sys})}$ measurement while assuming that $\tau = 0.05$ and $T_\mathrm{amb} = 273.15$ K locally.

2) $\left( T_\mathrm{sys}-T_\mathrm{rx} \right)/T_\mathrm{atm} >1$. This can occur when the fit underestimated 
$T_\mathrm{rx}$ or when $T_\mathrm{atm}$ is erroneously small. This typically takes place when $T_\mathrm{sys}$ is large --
i.e. at low elevations with a high opacity. If this happens, the code will assume $\tau = 1$ instead of using \autoref{taufromtsys}.

Additionally, whenever $\tau$ exceeds 1.0 at a certain time for a certain source,
the code will use the previously determined value for the opacity.

The $T_\mathrm{rx}$ values determined from fits as described above can be overwritten by values manually entered
in a \texttt{trec} file (\autoref{filenames}). Obviously, this should be done when observatories provided accurate
$T_\mathrm{rx}$ values. Otherwise, the pipeline will print warnings if problems with the opacity 
correction are encountered. In that case, the user should definitely inspect the diagnostic plots for the 
$T_\mathrm{rx}$ fits. For problematic fits, a better receiver temperature estimate can generally be obtained by looking at the
fits and guessing what $T_\mathrm{rx}$ should be.

Once, a \texttt{trec} file has been written with improved estimates, quickmode (\autoref{quickmode}) can be used to
re-run \texttt{task\_tsys\_add\_exptau()}.

The \texttt{trec} files should be ASCII and contain lines with the following format:\\
$<$station code$>$ $<$receiver$>$ $<${}$T_\mathrm{rx}$ value in K$>$

For example:

\begin{lstlisting}
BR RCP 113.5
BR LCP 110.5
GB RCP 65.5
KP LCP 111.5
KP RCP 97.8
\end{lstlisting}

\newpage
\section{Source Models}
\label{smodel}

CASA will store source models in a \texttt{MODEL\_DATA} data column of the MS.
For many calibration solvers, the data is divided by the model.
Having a reliable source model therefore helps with the calibration.

By default, the pipeline will look for files with modelfile\_extensions in the workdir (\autoref{filenames}).
For every $<$source$>$ in the data,
the pipeline will search for $<$source$>$.$<$modelfile\_extension$>$ files
and use them as source model for that source.

An optimal calibration strategy would be to run the pipeline,
make first images from the products (using CASA's \texttt{tclean()} task for example, see \autoref{interactive}),
save the resulting source models (modelfiles are automatically generated by \texttt{tclean()}),
and run the pipeline again.
This second run will then automatically make use of the determined source models
for an optimized calibration (e.g., \texttt{fringefit()} or \texttt{polcal()}).
Then, new images can be made with the improved calibration.

\newpage
\section{Flagging Algorithms}
\label{flaggingalgo}

By default the code will first look through any available metadata from flagging information.
Then, the edges of the bandpasses will be flagged as specified in \texttt{picard/input/flagging.inp} (\autoref{inpparams}).
See also \autoref{overview} and \autoref{flags}.

Additionally, two experimental flagging algorithms are implemented, which worked satisfactorily for the EHT array.
However, they have not been tested thoroughly for any other array. It is noteworthy that the options in \texttt{picard/input/flagging.inp}
allow for flagging dry runs, where no written flags are actually applied. Every flagging step will write the flags to files specified in
\texttt{picard/input/flagging.inp}.

\texttt{fg\_autocorr\_vs\_freq()} looks for outliers in frequency space (channel outliers in station's bandpasses). It will flag outliers based on the difference
in amplitude values between individual (or groups of) channels compared to the median derivative in the autocorrelation frequency spectrum.

\texttt{fg\_autocorr\_vs\_time()} looks at the station based frequency-averaged autocorrelation spectrum as a function of
time and flags around integrations whose amplitude differs by too much from the median amplitude per scan.

Both algorithms are described in more detail in the function's docstrings in \\ \texttt{picard/pipe\_modules/flagging\_algorithm.py}.

\newpage
\section{Flag Versions}
\label{flags}

By default, the pipeline will save an initial blank flag version when run for the first time.
Following \autoref{overview}, there are multiple steps (marked by being underlined) which will write flags.
Moreover, when the calibration is applied with \texttt{applycal()} at the end of the pipeline,
all uncalibrated data will be flagged, while a backup is saved.

In quickmode (\textit{-q}) (\autoref{quickmode}), the current flag status is used by default.
If the pipeline is run not in quickmode, the flags are saved/reverted to the initial blank version by default.

These default approaches are overwritten with the \textit{-r} command line option. Using just \textit{-r}, flags are saved/reverted
to the initial version even in quickmode. Using \textit{-r a}, flags are reverted to the latest flag version from \texttt{applycal()},
i.e. the flags before \texttt{applycal()} flagged any data. Ergo, this version should contain flags from the explicit flagging steps of the pipeline exclusively.

Whenever -\textit{q} is used to redo certain calibration tasks\footnote{Every time the pipeline is run, all steps are printed to the terminal.
The steps which are enumerated by numbers (and not letters) are the tasks. This convention is also adopted for the list in \autoref{overview}.
Moreover, all tasks follow the naming convention \texttt{task\_$^*$} when printed to the terminal.}
of the pipeline, it is recommended to use \textit{-r a}, to have the flags from the flagging algorithms applied, while ignoring flags from previous \texttt{applycal()} runs.

%\newpage
%\section{Akin-Incremental Calibration}
%\label{akinincremental}

%Sometimes it is useful to do incremental calibration for the same calibration tasks.
%For the EHT for example, a possible calibration strategy involves four fringe-fitting stages:
%single band $\rightarrow$ multi-band $\rightarrow$ single band again $\rightarrow$ multi-band again.
%This strategy is employed because of the fast phase rotations at 230 GHz; see \ref{papers} for details.

%For this type of incremental calibration, it may be tempting to copy-paste calibration functions
%and just give them a new name. But this would lead to duplicate code and there is an easier option:

%Say you want to run \texttt{task\_foo()} two times with different values for the parameter \texttt{bar}:
%\begin{enumerate}
%\item Add \texttt{bar2 = something else} next to \texttt{bar = something} in one of the \\ \texttt{picard/input/*.inp} files (see \ref{inpparams}).
%\item Define a function \texttt{task\_foo2()} next to \texttt{task\_foo()} and add a \texttt{calib\_foo2} entry
%      next to the \texttt{calib\_foo} entry in \texttt{picard/input/array.inp}, with unique names for the calibration
%      tables and diagnostics folders (see \ref{strategy}).
%\item Add \texttt{overwrite\_bar = $'\,'$} as default parameter to the inputs for \texttt{task\_foo()}
%      (c.f. \texttt{task\_fringefit\_multi\_general()} in \texttt{phase\_calibration.py} as a practical example).
%      In the body of code for \texttt{task\_foo()}, use \texttt{\_inp\_params.bar} if \texttt{overwrite\_bar} is not set and
%      else, \texttt{overwrite\_bar} for the \texttt{bar} parameter.
%\item In \texttt{task\_foo2()}, just do \texttt{return task\_foo(...,overwrite\_bar=\_inp\_params.bar2)}.
%\item Add \texttt{task\_foo2} in the desired position to \texttt{all\_calibration\_steps} in \texttt{main\_picard.py}
%      after \texttt{task\_foo} to have the correct incremental calibration (see \ref{strategy}).
%\end{enumerate}

%Following this strategy, an arbitrary number of akin-incremental calibration steps for the same task can be trivially implemented.

\newpage
\section{Change of Calibration Strategy}
\label{strategy}

A new calibration strategy (for a different telescope for example) can easily be implemented by
adding a new block for \texttt{inp\_params.array\_type} in the main file of the pipeline, ideally with a
customized \texttt{picard/input/array\_finetune.inp} file (\autoref{inpparams}). Of course, it is also possible to change the order of calibration steps in \texttt{picard/main\_picard.py}.
For example, the amplitude calibration steps could be done after fringe-fitting. Generally -- for any calibration task -- all calibration
tables (if written successfully) from the previous steps will be applied on-the-fly where applicable.\footnote{
For some tasks like \texttt{accor()}, on-the-fly calibration is ignored. Also, the multi-band fringe-fitting is split between 
science targets and calibrators.}
A new calibration task can be defined in any module and if it obeys to the convention
\texttt{task\_taskname()}$\leftrightarrow$\texttt{calib\_taskname} outlined in \texttt{picard/input/array\_finetune.inp}, it can be wrapped in
\texttt{calibration.go\_calibrate()} for on-the-fly calibration, potential smoothing, and generation of
diagnostic plots.

\sloppy

Example:\\
Let's say you have developed a method to refine the bandpass calibration, which you want to apply after the standard bandpass calibration.
Then you should do the following:
\begin{enumerate}
\item Add a custom module, let's say \texttt{my\_functions.py} in the \texttt{pipe\_modules} folder. Add any required imports to that module.
\item Put your code in a function in that module, obeying the naming convention \texttt{task\_taskname()} -- let's say we call it \texttt{task\_super\_bandpass()}.
      Make sure that the function returns \texttt{True} after successful execution.
      The task must take \texttt{\_inp\_params, \_ms\_metadata, caltable, \_fly\_calib\_tables, \_fly\_calib\_interp,
      \_fly\_calib\_gainfd, \_mpi\_client} as non-default arguments
      and any additional number of default arguments to overwrite parameters if necessary.
      Any user-defined input parameter must be defined in a \texttt{picard/input/*.inp} file -- let's say you add a \texttt{normalize\_super\_bandpass = True}
      parameter to a \texttt{picard/input/my\_inputs.inp} file, which you can then retrieve as \texttt{\_inp\_params.normalize\_super\_bandpass} 
      in \texttt{task\_super\_bandpass()}. \texttt{\_ms\_metadata} can be used to retrieve any kind of required meta data (see the \texttt{diagnostics.py} module).
      \texttt{\_fly\_calib\_tables} and \texttt{\_fly\_calib\_interp} are the calibration tables from all previous calibration steps (see 3. below),
      that should be applied on-the-fly (which is trivial for CASA tasks). \texttt{caltable} is the name of the calibration table that will be written 
      by \texttt{task\_super\_bandpass()} (see 4. below).
\item Import \texttt{my\_functions.py} in \texttt{main\_picard.py} and add \texttt{my\_functions.task\_super\_bandpass} to \texttt{all\_calibration\_steps}
      for the arrays for which this calibration is to be done. Remember that the order within \texttt{all\_calibration\_steps} matters -- the steps are 
      executed in sequence while applying the calibration tables from previous steps on-the-fly.
\item Add the mandatory input parameters to the ``Block of calibration parameters for each calibration step" in \texttt{picard/input/array\_finetune.inp}.
      You have to add it as \texttt{calib\_super\_bandpass = ...} to define the name of the output calibration table, together with
      smoothing, diagnostics, and interpolation options for the calibration solutions.
\end{enumerate}

\newpage
\section{Code Updates}
\label{updates}

To always have the  latest version of the code, running \textit{git pull} should suffice.

As described under `Versioning' in the \texttt{README.md} file, no additional action is required when pulling changes that go along with
minor updates or patches.

For major updates, either \texttt{README.md} changed since a new recommended CASA version became available 
(run \textit{./setup.py} again in this case), or the files in \texttt{input\_template} (\autoref{inpparams}) are updated
(then \textit{git diff} should be examined -- you will have to adapt your local input files to the changes).


\newpage
\section{Frequently Asked Questions}
\label{faq}

\begin{enumerate}
\item Why are you doing a scalar bandpass? \\
The scalar bandpass is used to accurately correct for the amplitude bandpass. Corrections from the complex bandpass are limited by the SNR of the cross-correlations and therefore, data from all scans of the specified calibrator sources are aggregated for the calibration task. If a scalar bandpass correction has been performed, the complex bandpass is used to correct only the phases. The scalar bandpass calibration should be skipped in the presence of strong RFI that has not been flagged.
Additional configuration options are available with the \texttt{solvemode\_scalar\_bandpass} parameter in \texttt{picard/input/array\_finetune.inp}, e.g., one can choose to not obtain solutions scan or to skip the scalar bandpass entirely.

\item What is the difference between the \texttt{VLBAhi} and \texttt{VLBAlo} options for the \texttt{array\_type} input parameter?\\
When \texttt{VLBAhi} is selected, an additional opacity correction is performed (\autoref{opac}). For \texttt{VLBAlo}, the natively measured system temperatures are applied directly. I recommend to use \texttt{VLBAhi} for observing frequencies above 15 GHz, and \texttt{VLBAlo} below.

\item What is the difference between the \texttt{input} and \texttt{input\_template} folders?\\
This is described in \autoref{codeoverview} and \autoref{inpparams}.
\end{enumerate}

\newpage
\section{Known Issues}
\label{bugs}

\begin{enumerate}
\item Issue: CASA aborts with the message that a \texttt{$\sim$/.matplotlib} path does not exist or that a \texttt{$\sim$/.matplotlib/tex.cache} file does already exist.\\
      Solution: Start CASA or \pipe{} again.
\item Issue: `UnboundLocalError' when initializing the metadata collection. \\
      Possible solution: Verify that all sources set in the \texttt{observation.inp} input file are actually present in the
      data, especially when working with a subset of the data created with the CASA \texttt{split()} task.
\item Issue: `IOError: Too many open files'. \\
      Context: In the CASA MPI implementation, the MS is split into several sub-MSs. The number of sub-MSs created are limited
      by the number $N$ of files that can be opened at the same time (c.f. $N$ = \$ ulimit -n) when the data is loaded.
      If $N$ changes (data copied to different system or someone changed the system settings), the above IOError can occur. \\
      Solution: Raise $N$ to the number used when the data was loaded.
\item Issue: `fatal: Not a git repository: ...'. \\
      Solution: Use git pull instead of downloading a tarball of the repository.
\item Issue: `no array in row x of column MODEL\_DATA'. \\
      Context: This can happen when a source model is used.
      The CASA MPI implementation seems to handle source models incorrectly in some cases. \\
      Solution: Load the data again after setting \texttt{MS\_partitioning = None} in \texttt{picard/input/constants.inp}.
\item Issue: \pipe{} aborts with a message about `ORTE was unable to reliably start one or more daemons...'\\
      Solution: Verify that your \texttt{input/mpi\_host\_file} is correct and that all servers can be reached.
\item Issue: `Exception: MPI is not enabled'. \\
      Solution: Run \pipe{} with \textit{-n N\_cores}, where N\_cores $\geq$ 2.
\item The pipeline starts in a home folder of a different machine and may complain that no input files can be found. \\
      Solution: 1) Verify your \texttt{mpi\_host\_file} (\autoref{mpi}). Is the current machine that you want to use in there? \\
      Or 2): Run \pipe{} with \textit{-n N\_cores}, where N\_cores $\geq$ 2.
\item Issue when using docker or singularity: `All nodes which are allocated for this job are already filled.' \\
      Solution: Run \pipe{} with \textit{-n N\_cores}, where N\_cores $\geq$ 2.
\item Issue: Docker complains about permissions or other obscure errors. \\
      Solution: Use root-less docker without root/sudo. It is available starting with docker version 19.03.
\item Issue when using docker or singularity: `KeyError: 'getpwuid(): uid not found' \\
      Context: This can happen when there is a mismatch between user information on the local system and docker. \\
      Solution \#1: Use singularity instead of docker. \\
      Solution \#2: Also pass \textit{--env HOME=/data --user \$(id -u) -v /etc/passwd:/etc/passwd} to \textit{docker run} (see \texttt{README}). \\
      If that does not work (can happen for macOS): \\
      Solution \#3: It could be that the user information is not actually stored in the \mbox{\texttt{/etc/passwd}} file on your local system.
      For macOS, user information can be retrieved with the \textit{dscl} or \textit{nidump} commands and a custom \texttt{passwd} file can be created with the information about the user.
      That file can then be shared with docker by using \mbox{\textit{-v passwd:/etc/passwd}} instead of \textit{-v /etc/passwd:/etc/passwd} for \textit{docker run}.
\item Issue when using singularity for \pipe{} or for the generation of diagnostic plots with \href{https://github.com/haavee/jiveplot}{jplotter}: `Unknown image format/type'.\\
      Solution: Run the plotting command again. If that does not help, upgrade your singularity installation to version 3.x.x.
\item Issue: When processing data with docker or singularity containers, software from the host system, as specified in the PATH, are used instead of the software from inside the containers. \\
      Solution: Set the environment variable \textit{PYTHONNOUSERSITE=1}.
\item Issue: Trying to plot the solutions in a calibration table does not work and this message is printed: `Note: Either your CalTable pre-dates name-based selection, or
      does not (yet) support selection, or the MS associated
      with this cal table does not exist.  All antennas,
      fields, spws are being selected for plotting.' \\
      Context: Calibration tables are directly associated with the MS that they were obtained from. \\
      Solution: You will need to retain the original data structure on your file system.
      By default this is a working directory where the MS is present and the calibration table resides inside a \texttt{calibration\_tables/} subfolder.
\item Issue: Antenna names are `unknown' in an exported UVFITS file. \\
      Context: Sometimes the CASA exportuvfits() code messes up station names. \\
      Solution: Run the export data step again.
\item The flux density calibration gains and/or system temperatures from an ANTAB table are not applied [to some antennas]. \\
      Context: The formatting restrictions for ANTAB tables are more strict for CASA compared to AIPS. \\
      Solution: Make sure that the ANTAB format follows exactly the file format exemplified in \autoref{antab}.
\item Issue: `TypeError: nan\_to\_num() takes exactly 1 argument (2 given)', `Invalid rgb arg', or a raised AttributeError about a colormap when \pipe{} is generating diagnostic plots.\\
      Solution: Upgrade to the CASA version given in the README file.
\item Issue: Python OS Errors occur and no new terminals can be opened. \\
      Context: When \textit{mpicasa} is killed, zombie /tmp/CASA\_MPIServer processes can be left over. \\
      Solution: Verify that no instance of \pipe{} is running and kill all zombies with \mbox{\$ \textit{pkill -f /tmp/CASA\_MPIServer}}.
\item Issue: The imager of \texttt{picard/interactive\_utils.py} (\autoref{interactive}) crashes with errors related to the \texttt{tclean} task (e.g., `TypeError: Exception from task\_tclean'). \\
      Context: There may be compatibility issues with newer versions of CASA. \\
      Solution: Try again with an older CASA version, for example CASA 5.6.
\end{enumerate}

\subsection{Subarrays}

Subarrays are not yet fully supported in CASA. If scans overlap in time, the calibration solutions from one scan
can influence the other. In this case, the recommendation is to split out the antennas from the different subarrays
into new measurement sets and process these individually. If necessary, a join phase calibration can be done beforehand
on overlapping scans for calibration purposes. The phae calibration table should be applied to the data and the corrected
data should then be split out.

\newpage
\section{Publications -- for Detailed Information}
\label{papers}

\begin{itemize}
\item Paper describing the pipeline in detail:\newline\href{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}.\newline
The DOI is \href{https://doi.org/10.1051/0004-6361/201935181}{https://doi.org/10.1051/0004-6361/201935181}.
This paper contains plots which illustrate the secondary reference antenna selection for the exhaustive fringe search (\autoref{refantpick}), the fringe-fit solution interval optimization search (\autoref{ffopti}), and the opacity correction (\autoref{opac}) for example.
Please cite this paper when you are using the pipeline. The software is also on \textsc{ascl}: \href{https://ascl.net/1905.015}{https://ascl.net/1905.015}.
\end{itemize}

\newpage
\section{Comparisons with AIPS-based Calibration}
\label{casaaipscomp}

Appendix A of \href{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}{https://ui.adsabs.harvard.edu/abs/2019A\%26A...626A..75J}
shows some comparisons between an \pipe{}-based calibration and a calibration using standard AIPS calibration methods.
Here, some additional information regarding this comparison is given.

   \autoref{vlbaununcalib} shows examples of uncalibrated data from the VLBA BW0106 M87 experiment.
   Uncorrected post-correlation delays are clearly visible. \autoref{vlbacalib} shows the same data after
   calibration by \pipe{}. Edge channels are flagged, delays are taken out, bandpasses are flattened,
   and the phase scatter has been reduced as atmospheric phase variations have been calibrated.
   There is still some phase scatter present in the Owens Valley -- St. Croix baseline 
   due to the SNR-limited timescales on which intra-scan phase variations can be corrected for.
   The low SNR is also reflected in the amplitudes which are affected by decoherence from the
   scan-average. 
   The data corresponding to \autoref{vlbacalib} after 
   the standard AIPS calibration
   is shown in Figures \autoref{vlbacaliba1r} and \autoref{vlbacaliba2l}, separately for the RR and LL.
   The AIPS calibrated data is in good agreement with the \pipe{} calibrate data.
   
                        \begin{figure}[h]
   \centering\offinterlineskip
   \includegraphics[width=0.55\textwidth,angle=270,origin=c]{images/vlba_calib/KPMKuncal-cropped.pdf}
      \includegraphics[width=0.55\textwidth,angle=270,origin=c]{images/vlba_calib/OVSCuncal-cropped.pdf}
      \vspace{-1cm}
      \caption{Unprocessed raw data of M87. Scan-averaged phases and amplitudes are shown for the two spectral windows in the dataset as 
      a function of frequency for the parallel-hand correlations (color-coded) of two baselines (Kitt Peak -- 
      Mauna Kea in the top panel and Owens Valley -- St. Croix in the bottom panel).
      The plots were made with \href{https://github.com/haavee/jiveplot}{jplotter}
      as part of the diagnostics automatically generated by \pipe{}.
              }
         \label{vlbaununcalib}
   \end{figure}
   

                     \begin{figure}[h]
   \centering\offinterlineskip
   \includegraphics[width=0.55\textwidth,angle=270,origin=c]{images/vlba_calib/KPMK-cropped.pdf}
      \includegraphics[width=0.55\textwidth,angle=270,origin=c]{images/vlba_calib/OVSC-cropped.pdf}
            \vspace{-1cm}
      \caption{\pipe{}'s pipeline processed data of M87. Scan-averaged phases and amplitudes are shown for the two spectral windows in the dataset as 
      a function of frequency for the parallel-hand correlations (color-coded) of two baselines (Kitt Peak -- 
      Mauna Kea in the top panel and Owens Valley -- St. Croix in the bottom panel). Crosses indicate flagged data.
      The plots were made with \href{https://github.com/haavee/jiveplot}{jplotter}
      as part of the diagnostics automatically generated by \pipe{}.
              }
         \label{vlbacalib}
   \end{figure}

   \newpage
                     \begin{figure}[h]
   \centering\offinterlineskip
   \includegraphics[width=0.95\textwidth]{images/vlba_calib/KPMKaipsrcl.pdf}
   \includegraphics[width=0.95\textwidth]{images/vlba_calib/KPMKaipslcl.pdf}
      \caption{AIPS processed data of M87 corresponding to the data shown in the top panel
               (KP-MK baseline) of Fig.~\ref{vlbacalib}.
               The AIPS data is shown separately for the RR and LL correlations
               in the top and bottom panels respectively
               as a function of frequency channel here.
              }
         \label{vlbacaliba1r}
   \end{figure}
  
   \newpage
                     \begin{figure}[h]
   \centering\offinterlineskip
   \includegraphics[width=0.95\textwidth]{images/vlba_calib/SCOVaipsrcl.pdf}
   \includegraphics[width=0.95\textwidth]{images/vlba_calib/OVSCaipslcl.pdf}
      \caption{AIPS processed data of M87 corresponding to the data shown in the bottom panel
               (OV-SC baseline) of Fig.~\ref{vlbacalib}.
               The AIPS data is shown separately for the RR and LL correlations
               in the top and bottom panels respectively
               as a function of frequency channel here.
              }
         \label{vlbacaliba2l}
   \end{figure}


\end{document}
