function qpipeline(startTime, stopTime, parameterFile, frameCacheFile, ...
                   resultsDirectory)
% QPIPELINE Top level function for the Q transform burst search
%
% QPIPELINE applies the discrete Q transform to search for statistically
% significant transient events in data from interferometric
% gravitational-wave detectors.
%
% usage: qpipeline(startTime, stopTime, parameterFile, frameCacheFile, ...
%                  outputDirectory)
%
%   startTime           gps start time of analysis
%   stopTime            gps stop time of analysis
%   parameterFile       qpipeline parameter file
%   frameCacheFile      readframedata formatted frame cache file
%   resultsDirectory    directory to write results to
%
%
% QPIPELINE produces the following output files in the specified results
% directory.
%
%   livetime-<startTime>-<stopTime>.txt
%   <channelName>-<startTime>-<stopTime>.txt
%
% If not specified, QPIPELINE writes its results to the current directory
% and looks for the parameter file parameters.txt and frame cache file
% framecache.txt in the current directory.
%
% If the specified stop time is less than the specified start time, it is
% instead assumed to be the duration of the analysis.  Non-integer start
% and stop times are truncated to the nearest integer.
%
% See also QTILE, QREADCALIBRATION, QREADDATA, QRESAMPLE, QCONDITION,
% QTRANSFORM, QCOLLOCATE, QTHRESHOLD, QSELECT, QCOINCIDE, and QWRITETRIGGERS.

% Shourov K. Chatterji
% shourov@ligo.caltech.edu

% Leo C. Stein
% lstein@ligo.mit.edu

% $Id: qpipeline.m,v 1.31 2007/06/28 16:35:42 lstein Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             start analysis timer                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% start time of analysis
analysisTimer = clock;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 write header                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report analysis name and time
fprintf(1, 'QPipeline analysis\n');
fprintf(1, 'Run by %s on %s at %s\n', ...
        getenv('USER'), datestr(clock, 29), datestr(clock, 13));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         parse command line arguments                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'parsing command line arguments\n');

% verify correct number of input arguments
error(nargchk(2, 5, nargin));

% apply default arguments
if (nargin < 3) || isempty(parameterFile),
  parameterFile = 'parameters.txt';
elseif (nargin < 4) || isempty(frameCacheFile),
  frameCacheFile = 'framecache.txt';
elseif (nargin < 5) || isempty(resultsDirectory),
  resultsDirectory = '.';
end

% ensure numeric start and stop times
if isstr(startTime),
  startTime = str2num(startTime);
end
if isstr(stopTime),
  stopTime = str2num(stopTime);
end

% truncate to integer start and stop times
startTime = floor(startTime);
stopTime = floor(stopTime);

% handle duration argument in place of stop time
if stopTime < startTime,
  stopTime = startTime + stopTime;
end

% string start and stop times
startTimeString = int2str(startTime);
stopTimeString = int2str(stopTime);
segmentString = [startTimeString '-' stopTimeString];

% report command line arguments
fprintf(1, '  startTime:               %#09d\n', startTime);
fprintf(1, '  stopTime:                %#09d\n', stopTime);
fprintf(1, '  parameterFile:           %s\n', parameterFile);
fprintf(1, '  frameCacheFile:          %s\n', frameCacheFile);
fprintf(1, '  resultsDirectory:        %s\n', resultsDirectory);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               read parameters                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% reoprt status
fprintf(1, 'reading parameter file\n');

% open parameter file for reading
parameterFileID = fopen(parameterFile, 'r');

% initialize parameters
channelNames = [];
frameTypes = [];
timeShifts = [];
injectionNames = [];
injectionTypes = [];
injectionFactors = [];
injectionTimeShifts = [];
calibrationFiles = [];
sampleFrequency = [];
highPassCutoff = [];
lowPassCutoff = [];
whiteningDuration = [];
blockDuration = [];
extraBlockOverlap = [];
qRange = [];
frequencyRange = [];
maximumMismatch = [];
analysisMode = [];
outlierFactor = [];
falseEventRate = [];
falseVetoRate = [];
uncertaintyFactor = [];
applyVeto = [];
maximumSignificants = [];
maximumTriggers = [];
durationInflation = [];
bandwidthInflation = [];
triggerFields = [];
triggerFormat = [];
randomSeed = [];
debugLevel = [];

% begin loop over parameter file
while ~feof(parameterFileID),

  % read one line from parameter file
  parameterLine = fgetl(parameterFileID);

  % remove any comments
  commentIndices = findstr(parameterLine, '#');
  if ~isempty(commentIndices),
    parameterLine = parameterLine(1 : (commentIndices(1) - 1));
  end

  % remove leading and trailing blanks
  parameterLine = fliplr(deblank(fliplr(deblank(parameterLine))));

  % if empty line, skip to the next line
  if isempty(parameterLine),
    continue;
  end

  % locate field separator
  colonIndex = strfind(parameterLine, ':');

  % if field separator not located, report syntax error
  if isempty(colonIndex),
    error('syntax error processing parameter file:\n%s\n', ...
          parameterLine);
  end

  % parse parameter line
  colonIndex = colonIndex(1);
  parameterName = parameterLine(1 : colonIndex - 1);
  parameterValue = parameterLine((colonIndex + 1) : end);
  parameterName = fliplr(deblank(fliplr(deblank(parameterName))));
  parameterValue = fliplr(deblank(fliplr(deblank(parameterValue))));

  % report parameter settings
  fprintf(1, '  %-25s%s\n', [parameterName ':'], parameterValue);

  % assign parameters based on name
  switch parameterName,

    % assign parameters
    case 'channelNames',
      channelNames = eval(parameterValue);
    case 'frameTypes',
      frameTypes = eval(parameterValue);
    case 'timeShifts',
      timeShifts = eval(parameterValue);
    case 'injectionNames',
      injectionNames = eval(parameterValue);
    case 'injectionTypes',
      injectionTypes = eval(parameterValue);
    case 'injectionFactors',
      injectionFactors = eval(parameterValue);
    case 'injectionTimeShifts',
      injectionTimeShifts = eval(parameterValue);
    case 'calibrationFiles',
      calibrationFiles = eval(parameterValue);
    case 'sampleFrequency',
      sampleFrequency = eval(parameterValue);
    case 'highPassCutoff',
      highPassCutoff = eval(parameterValue);
    case 'lowPassCutoff',
      lowPassCutoff = eval(parameterValue);
    case 'whiteningDuration',
      whiteningDuration = eval(parameterValue);
    case 'blockDuration',
      blockDuration = eval(parameterValue);
    case 'extraBlockOverlap',
      extraBlockOverlap = eval(parameterValue);
    case 'qRange',
      qRange = eval(parameterValue);
    case 'frequencyRange',
      frequencyRange = eval(parameterValue);
    case 'maximumMismatch',
      maximumMismatch = eval(parameterValue);
    case 'analysisMode',
      analysisMode = eval(parameterValue);
    case 'outlierFactor',
      outlierFactor = eval(parameterValue);
    case 'falseEventRate',
      falseEventRate = eval(parameterValue);
    case 'falseVetoRate',
      falseVetoRate = eval(parameterValue);
    case 'uncertaintyFactor',
      uncertaintyFactor = eval(parameterValue);
    case 'applyVeto',
      applyVeto = eval(parameterValue);
    case 'maximumSignificants',
      maximumSignificants = eval(parameterValue);
    case 'maximumTriggers',
      maximumTriggers = eval(parameterValue);
    case 'durationInflation',
      durationInflation = eval(parameterValue);
    case 'bandwidthInflation',
      bandwidthInflation = eval(parameterValue);
    case 'triggerFields',
      triggerFields = eval(parameterValue);
    case 'triggerFormat',
      triggerFormat = eval(parameterValue);
    case 'randomSeed',
      randomSeed = eval(parameterValue);
    case 'debugLevel',
      debugLevel = eval(parameterValue);

    % handle unknown parameters
    otherwise,
     error('unknown parameter %s\n', parameterName);

  % end assign parameters based on name
  end

% end loop over parameter file entries
end

% close parameter file
fclose(parameterFileID);

% test for unspecified parameters
if isempty(channelNames),
  error('channelNames not specified');
end
if isempty(frameTypes),
  error('frameTypes not specified');
end
if isempty(sampleFrequency),
  error('sampleFrequency not specified');
end
if isempty(blockDuration),
  error('blockDuration not specified');
end
if isempty(qRange),
  error('qRange not specified');
end
if isempty(frequencyRange),
  error('frequencyRange not specified');
end
if isempty(maximumMismatch),
  error('maximumMismatch not specified');
end
if isempty(analysisMode),
  error('analysisMode not specified');
end
if isempty(outlierFactor),
  error('outlierFactor not specified');
end
if isempty(falseEventRate),
  error('falseEventRate not specified');
end
if isempty(maximumSignificants),
  error('maximumSignificants not specified');
end
if isempty(maximumTriggers),
  error('maximumTriggers not specified');
end
if isempty(durationInflation),
  error('durationInflation not specified');
end
if isempty(bandwidthInflation),
  error('bandwidthInflation not specified');
end
if any(strcmpi(analysisMode, {'h1h2', 'h1h2veto', 'collocated', 'coherent'})),
  if isempty(falseVetoRate),
    error('falseVetoRate not specified');
  end
  if isempty(uncertaintyFactor),
    error('uncertaintyFactor not specified');
  end
end

% force cell arrays and vectors
if ~iscell(channelNames),
  channelNames = mat2cell(channelNames, size(channelNames, 1), ...
                          size(channelNames, 2));
end
if ~iscell(frameTypes),
  frameTypes = mat2cell(frameTypes, size(frameTypes, 1), ...
                        size(frameTypes, 2));
end
if ~iscell(injectionNames) && ~isempty(injectionNames),
  injectionNames = mat2cell(injectionNames, size(injectionNames, 1), ...
                            size(injectionNames, 2));
end
if ~iscell(injectionTypes) && ~isempty(injectionTypes),
  injectionTypes = mat2cell(injectionTypes, size(injectionTypes, 1), ...
                            size(injectionTypes, 2));
end
if ~iscell(calibrationFiles) && ~isempty(calibrationFiles),
  calibrationFiles = mat2cell(calibrationFiles, size(calibrationFiles, 1), ...
                              size(calibrationFiles, 2));
end
if iscell(timeShifts) && ~isempty(timeShifts),
  timeShifts = [timeShifts{:}];
end
if iscell(injectionFactors) && ~isempty(injectionFactors),
  injectionFactors = [injectionFactors{:}];
end
if iscell(injectionTimeShifts) && ~isempty(injectionTimeShifts),
  injectionTimeShifts = [injectionTimeShifts{:}];
end
if ~iscell(triggerFields) && ~isempty(triggerFields),
  triggerFields = mat2cell(triggerFields, size(triggerFields, 1), ...
                           size(triggerFields, 2));
end

% number of channels
numberOfChannels = length(channelNames);

% apply default values
if isempty(timeShifts),
  timeShifts = zeros(1, numberOfChannels);
end
if isempty(calibrationFiles),
  calibrationFiles = cell(1, numberOfChannels);
  [calibrationFiles{:}] = deal('NONE');
end
if isempty(injectionNames),
  injectionNames = cell(1, numberOfChannels);
  [injectionNames{:}] = deal('NONE');
end
if isempty(injectionTypes),
  injectionTypes = cell(1, numberOfChannels);
  [injectionTypes{:}] = deal('NONE');
end
if isempty(injectionFactors),
  injectionFactors = zeros(1, numberOfChannels);
end
if isempty(injectionTimeShifts),
  injectionTimeShifts = zeros(1, numberOfChannels);
end
if isempty(falseVetoRate),
  falseVetoRate = 0;
end
if isempty(uncertaintyFactor),
  uncertaintyFactor = 0;
end
if isempty(applyVeto),
  applyVeto = true;
end
if isempty(extraBlockOverlap),
  extraBlockOverlap = 0;
end
if isempty(triggerFields),
  triggerFields = {'time', 'frequency', 'duration', 'bandwidth', 'normalizedEnergy'};
  if any(strcmpi(analysisMode, {'h1h2', 'h1h2veto', 'collocated', 'coherent'})),
    triggerFields{end + 1} = 'incoherentEnergy';
  end
end
if isempty(triggerFormat),
  triggerFormat = 'txt';
end
if isempty(randomSeed),
  randomSeed = sum(1e6 * clock);
end
if isempty(debugLevel),
  debugLevel = 1;
end

% set random number generator seeds
rand('state', randomSeed);
randn('state', randomSeed);

% channels with requested injections
injectionChannels = find(~strcmp(upper(injectionNames), 'NONE') & ...
                         ~strcmp(upper(injectionTypes), 'NONE') & ...
                         (injectionFactors ~= 0));

% force row arrays and vectors
channelNames = channelNames(:);
frameTypes = frameTypes(:);
calibrationFiles = calibrationFiles(:);
timeShifts = timeShifts(:);
injectionNames = injectionNames(:);
injectionTypes = injectionTypes(:);
injectionFactors = injectionFactors(:);
injectionTimeShifts = injectionTimeShifts(:);

% store requested channel names
requestedChannelNames = channelNames;

% list of sites
sites = unique(regexprep(channelNames, '.:.*$', ''));
numberOfSites = length(sites);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             validate parameters                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% validate number of frame types
if length(frameTypes) ~= numberOfChannels,
  error('number of frameTypes and channelNames are inconsistent');
end

% validate number of time shifts
if (length(timeShifts) ~= numberOfChannels) && ~isempty(timeShifts),
  error('number of timeShifts and channelNames are inconsistent');
end

% validate number of injection names
if (length(injectionNames) ~= numberOfChannels) && ~isempty(injectionNames),
  error('number of injectionNames and channelNames are inconsistent');
end

% validate number of injection types
if (length(injectionTypes) ~= numberOfChannels) && ~isempty(injectionTypes),
  error('number of injectionTypes and channelNames are inconsistent');
end

% validate number of injection factors
if (length(injectionFactors) ~= numberOfChannels) && ~isempty(injectionFactors),
  error('number of injectionFactors and channelNames are inconsistent');
end

% validate number of injection time shifts
if (length(injectionTimeShifts) ~= numberOfChannels) && ...
      ~isempty(injectionTimeShifts),
  error('number of injectionTimeShifts and channelNames are inconsistent');
end

% validate number of calibration names
if (length(calibrationFiles) ~= numberOfChannels) && ~isempty(calibrationFiles),
  error('number of calibrationFiles and channelNames are inconsistent');
end

% validate analysis mode
switch lower(analysisMode),
  case 'single',
  case 'coincident',
    if (numberOfChannels < 2),
      error([analysisMode ' analysis mode requires at least two channels']);
    end
  case {'h1h2', 'h1h2veto'},
    if ((numberOfChannels ~= 2) || ~all(strcmp(sites, 'H'))),
      error([analysisMode ' analysis mode requires two hanford detectors']);
    end
  case 'collocated',
    if (numberOfSites > 1),
      error([analysisMode ' analysis mode permits only one site']);
    end
    if (numberOfChannels < 2),
      error([analysisMode ' analysis mode requires at least two channels']);
    end
  case 'coherent',
    if (numberOfSites < 3),
      error([analysisMode ' analysis mode requires at least three sites']);
    end
  otherwise,
    error('unknown analysis mode %s\n', analysisMode);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              results directory                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'creating results directory\n', resultsDirectory);

% create output directory
unix(['mkdir -p ' resultsDirectory]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              tile signal space                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'tiling search space\n');

% generate Q transform tiling
transientFactor = 4;
tiling = qtile(blockDuration, qRange, frequencyRange, sampleFrequency, ...
               maximumMismatch, highPassCutoff, lowPassCutoff, ...
               whiteningDuration, transientFactor);

% determine minimum block overlap
minimumBlockOverlap = 2 * tiling.transientDuration + extraBlockOverlap;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            read frame file cache                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'reading frame cache\n');

% read frame file cache
frameCache = loadframecache(frameCacheFile);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               read calibration                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if calibration files are specified
if ~all(strcmp(upper(calibrationFiles), 'NONE')),

  % report status
  fprintf(1, 'reading calibration\n');

  % read calibration information
  responseFunctions = qreadcalibration(calibrationFiles, tiling, debugLevel);

% otherwise
else

  % set empty response functions
  responseFunctions = cell(numberOfChannels, 1);

% end test for calibration files
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              partition segment                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'partioning segment\n');

% segment duration
segmentDuration = stopTime - startTime;

% force time shift alignment with sample interval
timeShifts = round(timeShifts * sampleFrequency) / sampleFrequency;

% maximum and minimum time shifts
maximumTimeShift = max(timeShifts);
minimumTimeShift = min(timeShifts);

% total segment livetime loss due to time shifts
timeShiftLoss = maximumTimeShift - minimumTimeShift;

% if segment duration is less than block duration
if (segmentDuration - timeShiftLoss) < blockDuration,

  % write error to log file
  error('segment to short to analyze');

% otherwise, continue
end

% number of blocks in segments
numberOfBlocks = ceil((segmentDuration - timeShiftLoss - ...
                       minimumBlockOverlap) / ...
                      (blockDuration - minimumBlockOverlap));

% if only one block,
if numberOfBlocks == 1,

  % ignore block overlap
  blockOverlap = 0;

% otherwise,
else

  % overlap in seconds between blocks
  blockOverlap = (segmentDuration - timeShiftLoss - ...
                  blockDuration * numberOfBlocks) ./ ...
                 (1 - numberOfBlocks);
% continue
end

% report segment information
fprintf(1, '  segment duration:    %9.2f seconds\n', segmentDuration);
fprintf(1, '  time shift loss:     %9.2f seconds\n', timeShiftLoss);
fprintf(1, '  block duration:      %9.2f seconds\n', blockDuration);
fprintf(1, '  block overlap:       %9.2f seconds\n', blockOverlap);
fprintf(1, '  number of blocks:    %9u blocks\n', numberOfBlocks);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               analyzed blocks                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'identifying analyzed blocks\n');

% name of livetime file
livetimeFile = [resultsDirectory '/livetime_' segmentString '.txt'];

% read list of previously analyzed blocks
if exist(livetimeFile) == 2,
  analyzedBlocks = dataread('file', livetimeFile, '%u%*[^\n]');
else
  analyzedBlocks = [];
end

% list of blocks to analyze
blocksToAnalyze = setdiff(1 : numberOfBlocks, analyzedBlocks);

% report number of previously analyzed blocks
fprintf(1, '  analyzed blocks:     %9u blocks\n', length(analyzedBlocks));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              begin block loop                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialize error flag
errorFlag = 0;

% begin loop over blocks
for blockNumber = blocksToAnalyze,

  % start time of block analysis
  blockTimer = clock;

  % gps start time of block
  blockStartTime = startTime + maximumTimeShift + ...
                   (blockNumber - 1) * (blockDuration - blockOverlap);

  % force start time alignment with data samples
  blockStartTime = floor(blockStartTime * sampleFrequency) / ...
                   sampleFrequency;

  % gps stop time of block
  blockStopTime = blockStartTime + blockDuration;

  % report status
  fprintf(1, 'analyzing block %u of %u at %#019.9f\n', ...
          blockNumber, numberOfBlocks, blockStartTime);

  % reset channel names and number of channels
  channelNames = requestedChannelNames;
  numberOfChannels = length(channelNames);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                            read detector data                              %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  reading detector data\n');

  % read detector data
  [data, sampleFrequencies] = ...
      qreaddata(frameCache, channelNames, frameTypes, blockStartTime, ...
                blockStopTime, timeShifts, debugLevel);

  % test for missing detector data
  if any(sampleFrequencies == 0),
    fprintf(1, 'ERROR: error reading detector data\n');
    errorFlag = 1;
    continue;
  end

  % standardize channel names
  channelNames = strrep(channelNames, ';', ':');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                           read injection data                              %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % if any injections are requested
  if any(injectionChannels),

    % report status
    fprintf(1, '  reading injection data\n');

    % read injection data
    [injectionData, injectionSampleFrequencies] = ...
        qreaddata(frameCache, injectionNames, injectionTypes, blockStartTime, ...
                  blockStopTime, injectionTimeShifts, debugLevel);

    % test for missing injection data
    if any(injectionSampleFrequencies(injectionChannels) == 0),
      fprintf(1, 'ERROR: error reading injection data\n');
      errorFlag = 1;
      continue;
    end

    % test for inconsistent injection data
    if any(injectionSampleFrequencies(injectionChannels) ~= ...
           sampleFrequencies(injectionChannels)),
      error('inconsistent detector and injection data');
    end

    % add injection data to detector data
    for channelNumber = injectionChannels,
      data{channelNumber} = data{channelNumber} + ...
          injectionFactors(channelNumber) * injectionData{channelNumber};
    end

    % free injection memory
    clear injectionData;
    clear injectionSampleFrequencies;

  % end test for requested injections
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                        apply h1h2 analysis modes                           %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % if either h1h2 analysis mode
  if strncmpi(analysisMode, 'h1h2', 4),

    % report status
    fprintf(1, '  summing data\n');

    % compute sum and difference time series
    data{3} = data{1} - data{2};
    data{4} = data{1} + data{2};

    % define sum and difference channel names
    channelNames{3} = [channelNames{1}(1) ':H1H2-NULL'];
    channelNames{4} = [channelNames{1}(1) ':H1H2-SUM'];

    % sample frequencies of sum and difference time series
    sampleFrequencies(3) = sampleFrequencies(1);
    sampleFrequencies(4) = sampleFrequencies(2);

    % define placeholder response functions
    responseFunctions{3} = [];
    responseFunctions{4} = [];

    % increase the number of channels
    numberOfChannels = 4;

  % end test for either h1h2 analysis mode
  end

  % if h1h2veto analysis mode
  if strcmpi(analysisMode, 'h1h2veto'),

    % discard single detector data
    data = data(3:4);
    channelNames = channelNames(3:4);
    sampleFrequencies = sampleFrequencies(3:4);
    responseFunctions = responseFunctions(3:4);
    numberOfChannels = 2;

  % end test for h1h2veto analysis mode
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              resample data                                 %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  resampling data\n');

  % resample data
  data = qresample(data, sampleFrequencies, sampleFrequency);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              condition data                                %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  conditioning data\n');

  % if either h1h2 analysis mode
  if strncmpi(analysisMode, 'h1h2', 4),

    % initialize calibration structure
    calibration = cell(size(data));

    % condition coherent sum channels
    [data(1:2), calibration(1:2)] = ...
        qcondition(data(1:2), tiling, responseFunctions(1:2));

    % temporarily disable whitening for null stream channels
    whiteningDuration = tiling.whiteningDuration;
    tiling.whiteningDuration = 0;

    % condition null stream channels
    [data(3:4), calibration(3:4)] = ...
        qcondition(data(3:4), tiling, responseFunctions(3:4));

    % re-enable whitening
    tiling.whiteningDuration = whiteningDuration;

  % otherwise
  else

    % condition data
    [data, calibration] = qcondition(data, tiling, responseFunctions);
    
  % end test for either h1h2 analysis mode
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              transform data                                %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  transforming data\n');

  % discrete q transform of data
  transforms = qtransform(data, tiling, calibration, outlierFactor);

  % free data memory
  clear data;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                           display spectrograms                             %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   % report status
%   fprintf(1, '  displaying spectrograms\n');

%   % spectrogram parameters
%   referenceTime = blockStartTime + tiling.transientDuration;
%   timeRange = [0 blockDuration - 2 * tiling.transientDuration];
%   normalizedEnergyRange = [0 25.5];

%   % display spectrograms
%   qspectrogram(transforms, tiling, blockStartTime, referenceTime, ...
%                timeRange, frequencyRange, qRange, normalizedEnergyRange, ...
%                channelNames);

%   % flush display
%   drawnow;

%   % wait for user input
%   pause;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                           apply analysis mode                              %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % apply analysis mode
  switch lower(analysisMode),

    % handle case of single detector analysis
    case 'single',

    % handle case of coincident analysis
    case 'coincident',

    % handle case of h1h2veto analysis
    case 'h1h2veto',

    % handle case of h1h2 analysis
    case 'h1h2',

      % report status
      fprintf(1, '  collocating data\n');

      % produce coherent sum and null stream transforms for collocated data
      [transforms(1:2), channelNames(1:2)] = ...
          qcollocate(transforms(1:2), tiling, channelNames(1:2));

      % discard post-transform collocated null stream
      transforms = transforms([1 3 4]);
      channelNames = channelNames([1 3 4]);
      numberOfChannels = length(channelNames);

    % handle case of collocated detectors
    case 'collocated',

      % report status
      fprintf(1, '  collocating data\n');

      % produce coherent sum and null stream transforms for collocated data
      [transforms, channelNames] = qcollocate(transforms, tiling, channelNames);

    % handle case of non-aligned detectors
    case 'coherent',

      % report error
      error('coherent analysis not supported');

    % unknown analysis mode
    otherwise,

      % report error
      error('unknown analysis mode %s\n', analysisMode);

  % end test for analysis mode
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                           display spectrograms                             %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   % report status
%   fprintf(1, '  displaying spectrograms\n');

%   % spectrogram parameters
%   referenceTime = blockStartTime + tiling.transientDuration;
%   timeRange = [0 blockDuration - 2 * tiling.transientDuration];
%   normalizedEnergyRange = [];

%   % display spectrograms
%   qspectrogram(transforms, tiling, blockStartTime, referenceTime, ...
%                timeRange, frequencyRange, qRange, normalizedEnergyRange, ...
%                channelNames);

%   % flush display
%   drawnow;

%   % wait for user input
%   pause;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              threshold data                                %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  thresholding data\n');

  % identify significant tiles
  thresholdReferenceTime = [];
  thresholdTimeRange = [];
  thresholdFrequencyRange = [];
  thresholdQRange = [];
  [triggers, channelNames] = ...
      qthreshold(transforms, tiling, blockStartTime, falseEventRate, ...
                 thresholdReferenceTime, thresholdTimeRange, ...
                 thresholdFrequencyRange, thresholdQRange, ...
                 maximumSignificants, analysisMode, falseVetoRate, ...
                 uncertaintyFactor, applyVeto, channelNames, debugLevel);

  % number of trigger channels
  numberOfChannels = length(channelNames);

  % free transform memory
  clear transforms;

  % report significant tile numbers
  for channelNumber = 1 : numberOfChannels,
    fprintf(1, '    %-21s %6u tiles\n', [channelNames{channelNumber} ':'], ...
            length(triggers{channelNumber}.time));
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                            display eventgrams                              %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   % report status
%   fprintf(1, '  displaying eventgrams\n');

%   % eventgram parameters
%   referenceTime = blockStartTime + tiling.transientDuration;
%   timeRange = [0 blockDuration - 2 * tiling.transientDuration];
%   normalizedEnergyRange = [];

%   % display eventgrams
%   qeventgram(triggers, tiling, blockStartTime, referenceTime, ...
%              timeRange, frequencyRange, ...
%              durationInflation, bandwidthInflation, ...
%              normalizedEnergyRange, channelNames);

%   % flush display
%   drawnow;

%   % wait for user input
%   pause;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                           downselect triggers                              %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  downselecting triggers\n');

  % excluding overlapping significant tiles
  triggers = qselect(triggers, durationInflation, bandwidthInflation, ...
                     maximumTriggers, channelNames, debugLevel);

  % report downselected tile numbers
  for channelNumber = 1 : numberOfChannels,
    fprintf(1, '    %-21s %6u tiles\n', [channelNames{channelNumber} ':'], ...
            length(triggers{channelNumber}.time));
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                            display eventgrams                              %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   % report status
%   fprintf(1, '  displaying eventgrams\n');

%   % eventgram parameters
%   referenceTime = blockStartTime + tiling.transientDuration;
%   timeRange = [0 blockDuration - 2 * tiling.transientDuration];
%   normalizedEnergyRange = [];

%   % display eventgrams
%   qeventgram(triggers, tiling, blockStartTime, referenceTime, ...
%              timeRange, frequencyRange, ...
%              durationInflation, bandwidthInflation, ...
%              normalizedEnergyRange, channelNames);

%   % flush display
%   drawnow;

%   % wait for user input
%   pause;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              write triggers                                %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % write status to log file
  fprintf(1, '  writing triggers\n');

  % determine file names
  triggerFiles = cell(1, numberOfChannels);
  for channelNumber = 1 : numberOfChannels,
    triggerFiles{channelNumber} = ...
        [resultsDirectory '/' channelNames{channelNumber} '_' ...
         segmentString '.' triggerFormat];
  end

  % write triggers to file
  qwritetriggers(triggers, triggerFiles, triggerFields, triggerFormat);

  % free triggers memory
  clear triggers;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              write livetime                                %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % write status to log file
  fprintf(1, '  writing livetime\n');

  % open livetime file for appending
  livetimeFileFID = fopen(livetimeFile, 'at');

  % test for error
  if livetimeFileFID == -1,
    error('cannot open livetime file for writing');
  end

  % write block livetime to livetime file
  fprintf(livetimeFileFID, '%#04u %#019.9f %#019.9f\n', blockNumber, ...
          blockStartTime + tiling.transientDuration, ...
          blockStopTime - tiling.transientDuration);

  % close livetime file
  fclose(livetimeFileFID);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                              end block loop                                %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % report status
  fprintf(1, '  block complete\n');

  % report total elapsed time
  fprintf(1, '    elapsed time:          %5.0f seconds\n', ...
          etime(clock, blockTimer));

% end loop over blocks
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 write footer                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'segment complete\n');

% report total elapsed time
fprintf(1, '  elapsed time:            %5.0f seconds\n', ...
        etime(clock, analysisTimer));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    return                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% exit with error status
if errorFlag,
  error('analysis had errors.');
else
  return;
end