function qevent(eventTime, eventDuration, parameterFile, frameCacheFile, ...
                outputDirectory)
% QEVENT Diagnostic version of QPipeline analysis for a single event
%
% QEVENT is a modified version of the Q Pipeline that graphically displays
% intermediate data products around a specified event time.
%
% usage: qevent(eventTime, eventDuration, parameterFile, frameCacheFile, ...
%               outputDirectory);
%
%   eventTime           gps center time of analysis
%   eventDuration       duration of analysis
%   parameterFile       qpipeline parameter file
%   frameCacheFile      readframedata formatted frame cache file
%   outputDirectory     path name of directory to write results in
%
% QEVENT places the resulting image files and an index.html file for web
% based viewing of the results in a subdirectory of the specified output
% directory.  The subdirectory is named by the specified event time.  If no
% output directory is specified, QEVENT uses the current working directory
% by default.
%
% If a parameter file or frame cache file are not specified, QEVENT looks
% for the parameter file parameters.txt and frame cache file framecache.txt
% in the current directory.
%
% See also QPIPELINE, 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: qevent.m,v 1.29 2007/07/03 00:27:26 shourov Exp $

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

% start time of analysis
analysisTimer = clock;

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

% report analysis name and time
fprintf(1, 'QEvent 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(outputDirectory),
  outputDirectory = '.';
end

% ensure numeric event time and duration
if isstr(eventTime),
  eventTimeString = eventTime;
  eventTime = str2num(eventTime);
else
  eventTimeString = sprintf('%#019.9f', eventTime);
end
if isstr(eventDuration),
  eventDuration = str2num(eventDuration);
end

% results subdirectory for this event
eventDirectory = [outputDirectory '/' eventTimeString];

% report command line arguments
fprintf(1, '  eventTime:               %#019.9f\n', eventTime);
fprintf(1, '  eventDuration:           %#019.9f\n', eventDuration);
fprintf(1, '  parameterFile:           %s\n', parameterFile);
fprintf(1, '  frameCacheFile:          %s\n', frameCacheFile);
fprintf(1, '  outputDirectory:         %s\n', outputDirectory);
fprintf(1, '  eventDirectory:          %s\n', eventDirectory);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               event directory                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% report status
fprintf(1, 'creating event directory\n', eventDirectory);

% create event directory
unix(['mkdir -p ' eventDirectory]);

% copy configuration file
unix(['cp ' parameterFile ' ' eventDirectory '/parameters.txt']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               open html report                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% open html report
htmlFile = [eventDirectory '/index.html'];
htmlFID = fopen(htmlFile, 'wt');

% begin html report
fprintf(htmlFID, '<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN"\n');
fprintf(htmlFID, '"http://www.w3.org/TR/html4/strict.dtd">\n');
fprintf(htmlFID, '<html>\n');
fprintf(htmlFID, '<head>\n');
fprintf(htmlFID, '<title>Event at %s</title>\n', eventTimeString);
fprintf(htmlFID, '<link rel="stylesheet" href="qstyle.css" type="text/css">\n');
fprintf(htmlFID, '</head>\n');
fprintf(htmlFID, '<body>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div class="title">\n');
fprintf(htmlFID, '<h1>Event at %s</h1>\n', eventTimeString);
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div class="index">');
fprintf(htmlFID, '<h2>Index</h2>');
fprintf(htmlFID, '<p>');
fprintf(htmlFID, '<a href="parameters.txt">');
fprintf(htmlFID, 'Parameters</a><br />');
fprintf(htmlFID, '<a href="log.txt">');
fprintf(htmlFID, 'Scan Log</a><br />');
fprintf(htmlFID, '<a href="http://www.ligo.caltech.edu/~shourov/q/qevent/">');
fprintf(htmlFID, 'Documentation</a><br />');
fprintf(htmlFID, '</p>');
fprintf(htmlFID, '</div>');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div class="content">');
fprintf(htmlFID, '\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               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 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         segment start and stop times                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

% gps start time of block
blockStartTime = floor(eventTime - blockDuration / 2);

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

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

% 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),
  error('error reading detector data');
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),
    error('error reading injection data');
  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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             produce spectrograms                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% spectrogram parameters
referenceTime = eventTime;
timeRange = [-1 +1] * eventDuration / 2;
normalizedEnergyRange = [];

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

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

% print spectrograms and append to html report
fprintf(htmlFID, '<h2>QSpectrograms of individual detector data</h2>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div>\n');
fprintf(htmlFID, '\n');
for channelNumber = 1 : numberOfChannels,
  fprintf(htmlFID, '<h3>%s</h3><br />', channelNames{channelNumber});
  fprintf(htmlFID, '<div>\n');
  baseName = [eventTimeString '_' channelNames{channelNumber} '_spectrogram_'];
  for plane = 1 : tiling.numberOfPlanes,
    figure((channelNumber - 1) * tiling.numberOfPlanes + plane);
    figureFile = [baseName int2str(plane) '.png'];
    print('-dpng', [eventDirectory '/' figureFile]);
    thumbnailFile = [baseName int2str(plane) '_thumbnail.png'];
    unix(['convert -resize 300x ' eventDirectory '/' figureFile ...
          ' -strip -depth 8 -colors 256 ' eventDirectory '/' thumbnailFile]);
    fprintf(htmlFID, '<a href="%s">\n<img src="%s"></a>\n', ...
            figureFile, thumbnailFile);
  end
  fprintf(htmlFID, '</div>');
end
fprintf(htmlFID, '\n');
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             produce spectrograms                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if analysis mode is collocated or coherent
if any(strcmpi(analysisMode, {'h1h2', 'collocated', 'coherent'})),

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

  % spectrogram parameters
  referenceTime = eventTime;
  timeRange = [-1 +1] * eventDuration / 2;
  normalizedEnergyRange = [];

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

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

  % print spectrograms and append to html report
  fprintf(htmlFID, '<h2>QSpectrograms of coherent sum and null streams</h2>\n');
  fprintf(htmlFID, '\n');
  fprintf(htmlFID, '<div>\n');
  fprintf(htmlFID, '\n');
  for channelNumber = 1 : numberOfChannels,
    fprintf(htmlFID, '<h3>%s</h3><br />', channelNames{channelNumber});
    fprintf(htmlFID, '<div>\n');
    baseName = [eventTimeString '_' channelNames{channelNumber} '_spectrogram_'];
    for plane = 1 : tiling.numberOfPlanes,
      figure((channelNumber - 1) * tiling.numberOfPlanes + plane);
      figureFile = [baseName int2str(plane) '.png'];
      print('-dpng', [eventDirectory '/' figureFile]);
      thumbnailFile = [baseName int2str(plane) '_thumbnail.png'];
      unix(['convert -resize 300x ' eventDirectory '/' figureFile ...
            ' -strip -depth 8 -colors 256 ' eventDirectory '/' thumbnailFile]);
      fprintf(htmlFID, '<a href="%s">\n<img src="%s"></a>\n', ...
              figureFile, thumbnailFile);
    end
    fprintf(htmlFID, '</div>');
  end
  fprintf(htmlFID, '\n');
  fprintf(htmlFID, '</div>\n');
  fprintf(htmlFID, '\n');

% end test for collocated analysis mode
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            producing eventgrams                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% eventgram parameters
referenceTime = eventTime;
timeRange = [-1 +1] * eventDuration / 2;
normalizedEnergyRange = [];

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

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

% print eventgrams and append to html report
fprintf(htmlFID, '<h2>QEventgrams of significant tiles</h2>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div>\n');
fprintf(htmlFID, '\n');
for channelNumber = 1 : numberOfChannels,
  fprintf(htmlFID, '<h3>%s</h3><br />', channelNames{channelNumber});
  fprintf(htmlFID, '<div>\n');
  figure(channelNumber);
  baseName = [eventTimeString '_' channelNames{channelNumber} '_significants'];
  figureFile = [baseName '.png'];
  print('-dpng', [eventDirectory '/' figureFile]);
  thumbnailFile = [baseName '_thumbnail.png'];
  unix(['convert -resize 300x ' eventDirectory '/' figureFile ...
        ' -strip -depth 8 -colors 256 ' eventDirectory '/' thumbnailFile]);
  fprintf(htmlFID, '<a href="%s">\n<img src="%s"></a>\n', ...
          figureFile, thumbnailFile);
  fprintf(htmlFID, '</div>\n');
end
fprintf(htmlFID, '\n');
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            producing eventgrams                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% eventgram parameters
referenceTime = eventTime;
timeRange = [-1 +1] * eventDuration / 2;
normalizedEnergyRange = [];

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

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

% print eventgrams and append to html report
fprintf(htmlFID, '<h2>QEventgrams of non-overlapping significant tiles</h2>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div>\n');
fprintf(htmlFID, '\n');
for channelNumber = 1 : numberOfChannels,
  fprintf(htmlFID, '<h3>%s</h3><br />', channelNames{channelNumber});
  fprintf(htmlFID, '<div>\n');
  figure(channelNumber);
  baseName = [eventTimeString '_' channelNames{channelNumber} '_mosaics'];
  figureFile = [baseName '.png'];
  print('-dpng', [eventDirectory '/' figureFile]);
  thumbnailFile = [baseName '_thumbnail.png'];
  unix(['convert -resize 300x ' eventDirectory '/' figureFile ...
        ' -strip -depth 8 -colors 256 ' eventDirectory '/' thumbnailFile]);
  fprintf(htmlFID, '<a href="%s">\n<img src="%s"></a>\n', ...
          figureFile, thumbnailFile);
  fprintf(htmlFID, '</div>\n');
end
fprintf(htmlFID, '\n');
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');

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

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

% write triggers
fprintf(htmlFID, '<h2>text lists of non-overlapping significant tiles</h2>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div>\n');
fprintf(htmlFID, '\n');
triggerFiles = cell(1, numberOfChannels);
for channelNumber = 1 : numberOfChannels,
  triggerFiles{channelNumber} = ...
      [eventTimeString '_' channelNames{channelNumber} '.' triggerFormat];
  fprintf(htmlFID, '<h3><a href="%s">%s</a></h3><br />', ...
          triggerFiles{channelNumber}, channelNames{channelNumber});
  triggerFiles{channelNumber} = ...
      [eventDirectory '/' triggerFiles{channelNumber}];
end
qwritetriggers(triggers, triggerFiles, triggerFields, triggerFormat);
fprintf(htmlFID, '\n');
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');

% free triggers memory
clear triggers;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              close html report                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% end html report
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '<div class="footer">\n');
fprintf(htmlFID, 'Created by user %s on %s at %s<br />\n', ...
        getenv('USER'), datestr(clock, 29), datestr(clock, 13));
fprintf(htmlFID, '</div>\n');
fprintf(htmlFID, '\n');
fprintf(htmlFID, '</body>\n');
fprintf(htmlFID, '</html>\n');

% close html report
fclose(htmlFID);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 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                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% close all figures and quit
exit