function qproperties(startTime, stopTime, parametersFile, triggersFile, ...
                     clustersFile, livetimeFile, channelName, timeStamp, ...
                     outputDirectory)
% QPROPERTIES display and print plots of trigger properties
%
% QPROPERTIES reads the specified QPipeline unclustered and clustered
% trigger files and produces plots of the followng trigger properties
% for both the clustered and unclustered triggers
%
%   trigger rate as a function of snr threshold
%   histogrammed trigger rate as a function of time
%   scatter plot of trigger center frequency vs. center time
%   scatter plot of trigger snr vs. center time
%
% usage:
%
%   qproperties(startTime, stopTime, parametersFile, triggersFile, ...
%               clustersFile, livetimeFile, channelName, timeStamp, ...
%               outputDirectory);
%
%   startTime        start time for computing and displaying properties
%   stopTime         stop time for computing and displaying properties
%   parametersFile   qpipeline parameter file
%   triggersFile     unclustered triggers file
%   clustersFile     clustered triggers file
%   livetimeFile     livetime segment file
%   channelName      channel name
%   timeStamp        time stamp corresponding to start time
%   outputDirectory  output directory
%
% The livetime file must consist of non-overlapping livetime segments
% sorted by increasing start time.  The livetime file should be a two
% column ASCII text list of start and stop times.
%
% The triggers and clusters files are 5 column ASCII lists of the following
% trigger properties.
%
%   central time [gps seconds]
%   central frequency [Hz]
%   duration [seconds]
%   bandwidth [Hz]
%   normalized energy []
%
% Additional columns are permitted, and will simply be ignored.
%
% The parameter file should be a QPipeline parameter file.

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

% $Id: qproperties.m,v 1.12 2007/06/11 09:39:18 shourov Exp $

% convert string times to numbers
if isstr(startTime),
  startTimeString = startTime;
  startTime = str2num(startTimeString);
else
  startTimeString = num2str(startTime);
end
if isstr(stopTime),
  stopTimeString = stopTime;
  stopTime = str2num(stopTimeString);
else
  stopTimeString = num2str(stopTime);
end

% latex compatible channel name
channelName = strrep(channelName, '_', '\_');

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

% signal to noise ratio threshold
rho0 = 5;

% load livetime segments
disp('loading livetime segments...');
livetime = dataread('file', livetimeFile, '', 'commentstyle', 'shell');

% restrict livetime to requested time period
disp('masking livetime segments...');
if ~isempty(livetime),
  keepIndices = find((livetime(:, 1) < stopTime) & ...
                     (livetime(:, 2) > startTime));
  livetime(:, 1) = max(startTime, livetime(:, 1));
  livetime(:, 2) = min(stopTime, livetime(:, 2));
  livetime = livetime(keepIndices, :);
end

% number of livetime segments to display
numberOfSegments = size(livetime, 1);

% total livetime to display
analysisLivetime = max(1, sum(diff(livetime, 1, 2)));

% load unclustered triggers
disp('loading unclustered triggers...');
[triggersTime, triggersFrequency, triggersDuration, ...
 triggersBandwidth, triggersNormalizedEnergy] = ...
    dataread('file', triggersFile, '%f %f %f %f %f%*[^\n]', ...
             'commentstyle', 'matlab');
triggersQ = 2 * sqrt(pi) * triggersFrequency ./ triggersBandwidth;
triggersRho = sqrt(2 * triggersNormalizedEnergy);
clear triggersDuration;
clear triggersBandwidth;
clear triggersNormalizedEnergy;

% load clustered triggers
disp('loading clustered triggers...');
[clustersTime, clustersFrequency, clustersDuration, ...
 clustersBandwidth, clustersNormalizedEnergy] = ...
    dataread('file', clustersFile, '%f %f %f %f %f%*[^\n]', ...
             'commentstyle', 'matlab');
clustersQ = 2 * sqrt(pi) * clustersFrequency ./ clustersBandwidth;
clustersRho = sqrt(2 * clustersNormalizedEnergy);
clear clustersDuration;
clear clustersBandwidth;
clear clustersNormalizedEnergy;

% restrict unclustered triggers to requested time period
disp('masking unclustered triggers...');
keepIndices = find((triggersTime >= startTime) & ...
                   (triggersTime < stopTime));
triggersTime = triggersTime(keepIndices);
triggersFrequency = triggersFrequency(keepIndices);
triggersQ = triggersQ(keepIndices);
triggersRho = triggersRho(keepIndices);

% restrict clustered triggers to requested time period
disp('masking clustered triggers...');
keepIndices = find((clustersTime >= startTime) & ...
                   (clustersTime < stopTime));
clustersTime = clustersTime(keepIndices);
clustersFrequency = clustersFrequency(keepIndices);
clustersQ = clustersQ(keepIndices);
clustersRho = clustersRho(keepIndices);

% sort unclustered triggers
disp('sorting unclustered triggers...');
[ignore, sortedIndices] = sort(triggersRho);
triggersTime = triggersTime(sortedIndices);
triggersFrequency = triggersFrequency(sortedIndices);
triggersQ = triggersQ(sortedIndices);
triggersRho = triggersRho(sortedIndices);

% sort clustered triggers
disp('sorting clustered triggers...');
[ignore, sortedIndices] = sort(clustersRho);
clustersTime = clustersTime(sortedIndices);
clustersFrequency = clustersFrequency(sortedIndices);
clustersQ = clustersQ(sortedIndices);
clustersRho = clustersRho(sortedIndices);

% threshold unclustered triggers
disp('thresholding unclustered triggers...');
keepIndices = find(triggersRho > rho0);
numberOfTriggers = length(keepIndices);
triggersTime = triggersTime(keepIndices);
triggersFrequency = triggersFrequency(keepIndices);
triggersQ = triggersQ(keepIndices);
triggersRho = triggersRho(keepIndices);

% threshold clustered triggers
disp('thresholding clustered triggers...');
keepIndices = find(clustersRho > rho0);
numberOfClusters = length(keepIndices);
clustersTime = clustersTime(keepIndices);
clustersFrequency = clustersFrequency(keepIndices);
clustersQ = clustersQ(keepIndices);
clustersRho = clustersRho(keepIndices);

% identify trigger subsets by snr threshold
triggers10Indices = find(triggersRho > 10);
triggers20Indices = find(triggersRho > 20);
[ignore, triggersLoudestIndex] = max(triggersRho);

% identify cluster subsets by snr threshold
clusters10Indices = find(clustersRho > 10);
clusters20Indices = find(clustersRho > 20);
[ignore, clustersLoudestIndex] = max(clustersRho);

% renormalize time scale
disp('renormalizing time scale...');
analysisDuration = stopTime - startTime;
if analysisDuration >= 3 * 24 * 60 * 60,
  timeScale = 24 * 60 * 60;
  timeLabel = 'days';
elseif analysisDuration >= 3 * 60 * 60,
  timeScale = 60 * 60;
  timeLabel = 'hours';
elseif analysisDuration >= 3 * 60,
  timeScale = 60;
  timeLabel = 'minutes';
else
  timeScale = 1;
  timeLabel = 'seconds';
end
triggersTime = (triggersTime - startTime) / timeScale;
clustersTime = (clustersTime - startTime) / timeScale;
livetime = (livetime - startTime) / timeScale;

% plot unclustered trigger rate vs. snr threshold
disp('plotting unclustered trigger rate vs. snr threshold...');
triggersRate = (numberOfTriggers : -1 : 1) / analysisLivetime;
clf;
set(gca, 'FontSize', 18);
handles = loglog(triggersRho, triggersRate, 'bo-');
set(gca, 'XScale', 'log');
set(gca, 'YScale', 'log');
set(handles, 'MarkerSize', 4);
set(handles, 'LineWidth', 2);
axis([rho0 1e3 10^floor(log10(1 / analysisDuration)) 1e1]);
grid on;
xlabel('Signal to noise ratio threshold');
ylabel('Trigger rate exceeding threshold [Hz]');
title([channelName ' unclustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/unclustered_snr_rate_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/unclustered_snr_rate_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/unclustered_snr_rate_' num2str(rho0) '_thumbnail.png']);

% plot clustered trigger rate vs. snr threshold
disp('plotting cumulative snr distribution...');
clustersRate = (numberOfClusters : -1 : 1) / analysisLivetime;
clf;
set(gca, 'FontSize', 18);
handles = loglog(clustersRho, clustersRate, 'bo-');
set(gca, 'XScale', 'log');
set(gca, 'YScale', 'log');
set(handles, 'MarkerSize', 4);
set(handles, 'LineWidth', 2);
axis([rho0 1e3 10^floor(log10(1 / analysisDuration)) 1e1]);
grid on;
xlabel('Signal to noise ratio threshold');
ylabel('Trigger rate exceeding threshold [Hz]');
title([channelName ' clustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/clustered_snr_rate_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/clustered_snr_rate_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/clustered_snr_rate_' num2str(rho0) '_thumbnail.png']);

% plot unclustered trigger rate vs. time
disp('plotting unclustered trigger rate vs. time...');
clf;
set(gca, 'FontSize', 18);
binWidth = 60 / timeScale;
binEdges = 0 : binWidth : analysisDuration / timeScale;
binCenters = binWidth / 2 : binWidth : analysisDuration / timeScale - binWidth / 2;
numberOfBins = length(binCenters);
if isempty(livetime),
  binLivetime = zeros(size(binCenters));
else
  binLivetime = ...
      sum(max(0, min(ones(numberOfSegments, 1) * binEdges(2 : end), ...
                     livetime(:, 2) * ones(1, numberOfBins)) - ...
                 max(ones(numberOfSegments, 1) * binEdges(1 : end - 1), ...
                     livetime(:, 1) * ones(1, numberOfBins))), 1);
end
binLivetime(find(binLivetime == 0)) = 1;
if isempty(triggersTime),
  triggersRate = zeros(length(binCenters), 1);
  triggersRateMinimum = zeros(length(binCenters), 1);
  triggersRateMaximum = zeros(length(binCenters), 1);
else
  triggersCount = histc(triggersTime, binEdges);
  triggersRate = triggersCount(1 : end - 1) ./ (binLivetime * timeScale).';
  triggersRateMinimum = triggersRate - ...
                        sqrt(triggersCount(1 : end - 1)) ./ ...
                        (binLivetime * timeScale).';
  triggersRateMaximum = triggersRate + ...
                        sqrt(triggersCount(1 : end - 1)) ./ ...
                        (binLivetime * timeScale).';
  triggersRateMinimum = max(1e-2, triggersRateMinimum);
  triggersRateMaximum = min(1e+2, triggersRateMaximum);
end
handles = semilogy([binCenters; binCenters;], ...
                   [triggersRateMinimum.'; triggersRateMaximum.';], ...
                   'k-', ...
                   binCenters, triggersRate.', 'bo');
set(handles(1 : end - 1), 'LineWidth', 1);
set(handles(end), 'MarkerSize', 5, ...
                  'MarkerEdgeColor', 'b', ...
                  'MarkerFaceColor', 'b');
grid on;
axis([0 analysisDuration / timeScale 1e-2 1e2])
xlabel(['Time [' timeLabel ']']);
ylabel('Trigger rate [Hz]');
title([channelName ' unclustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/unclustered_time_rate_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/unclustered_time_rate_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/unclustered_time_rate_' num2str(rho0) '_thumbnail.png']);

% plot clustered trigger rate vs. time
disp('plotting clustered trigger rate vs. time...');
clf;
set(gca, 'FontSize', 18);
binWidth = 60 / timeScale;
binEdges = 0 : binWidth : analysisDuration / timeScale;
binCenters = binWidth / 2 : binWidth : analysisDuration / timeScale - binWidth / 2;
numberOfBins = length(binCenters);
if isempty(livetime),
  binLivetime = zeros(size(binCenters));
else
  binLivetime = ...
      sum(max(0, min(ones(numberOfSegments, 1) * binEdges(2 : end), ...
                     livetime(:, 2) * ones(1, numberOfBins)) - ...
                 max(ones(numberOfSegments, 1) * binEdges(1 : end - 1), ...
                     livetime(:, 1) * ones(1, numberOfBins))), 1);
end
binLivetime(find(binLivetime == 0)) = 1;
if isempty(clustersTime),
  clustersRate = zeros(length(binCenters), 1);
  clustersRateMinimum = zeros(length(binCenters), 1);
  clustersRateMaximum = zeros(length(binCenters), 1);
else
  clustersCount = histc(clustersTime, binEdges);
  clustersRate = clustersCount(1 : end - 1) ./ (binLivetime * timeScale).';
  clustersRateMinimum = clustersRate - ...
                        sqrt(clustersCount(1 : end - 1)) ./ ...
                        (binLivetime * timeScale).';
  clustersRateMaximum = clustersRate + ...
                        sqrt(clustersCount(1 : end - 1)) ./ ...
                        (binLivetime * timeScale).';
  clustersRateMinimum = max(1e-2, clustersRateMinimum);
  clustersRateMaximum = min(1e+2, clustersRateMaximum);
end
handles = semilogy([binCenters; binCenters;], ...
                   [clustersRateMinimum.'; clustersRateMaximum.';], ...
                   'k-', ...
                   binCenters, clustersRate.', 'bo');
set(handles(1 : end - 1), 'LineWidth', 1);
set(handles(end), 'MarkerSize', 5, ...
                  'MarkerEdgeColor', 'b', ...
                  'MarkerFaceColor', 'b');
grid on;
axis([0 analysisDuration / timeScale 1e-2 1e2])
xlabel(['Time [' timeLabel ']']);
ylabel('Trigger rate [Hz]');
title([channelName ' clustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/clustered_time_rate_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/clustered_time_rate_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/clustered_time_rate_' num2str(rho0) '_thumbnail.png']);

% plot unclustered trigger snr vs. time
disp('plotting unclustered trigger snr vs. time...');
clf;
hold on;
set(gca, 'FontSize', 18);
handles1 = semilogy(triggersTime, ...
                    triggersRho, ...
                    'bo');
handles2 = semilogy(triggersTime(triggers10Indices), ...
                    triggersRho(triggers10Indices), ...
                    'go');
handles3 = semilogy(triggersTime(triggers20Indices), ...
                    triggersRho(triggers20Indices), ...
                    'ro');
handles4 = semilogy(triggersTime(triggersLoudestIndex), ...
                    triggersRho(triggersLoudestIndex), ...
                    'ko');
set(gca, 'XScale', 'linear');
set(gca, 'YScale', 'log');
set(handles1, 'MarkerSize', 3, ...
              'MarkerEdgeColor', 'b', ...
              'MarkerFaceColor', 'b');
set(handles2, 'MarkerSize', 6, ...
              'MarkerEdgeColor', 'g', ...
              'MarkerFaceColor', 'g');
set(handles3, 'MarkerSize', 9, ...
              'MarkerEdgeColor', 'r', ...
              'MarkerFaceColor', 'r');
set(handles4, 'MarkerSize', 20, 'LineWidth', 4);
axis([0 analysisDuration / timeScale rho0 1e3]);
xlabel(['Time [' timeLabel ']']);
ylabel('Signal to noise ratio');
title([channelName ' unclustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/unclustered_time_snr_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/unclustered_time_snr_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/unclustered_time_snr_' num2str(rho0) '_thumbnail.png']);

% plot clustered trigger snr vs. time
disp('plotting clustered trigger snr vs. time...');
clf;
hold on;
set(gca, 'FontSize', 18);
handles1 = semilogy(clustersTime, ...
                    clustersRho, ...
                    'bo');
handles2 = semilogy(clustersTime(clusters10Indices), ...
                    clustersRho(clusters10Indices), ...
                    'go');
handles3 = semilogy(clustersTime(clusters20Indices), ...
                    clustersRho(clusters20Indices), ...
                    'ro');
handles4 = semilogy(clustersTime(clustersLoudestIndex), ...
                    clustersRho(clustersLoudestIndex), ...
                    'ko');
set(gca, 'XScale', 'linear');
set(gca, 'YScale', 'log');
set(handles1, 'MarkerSize', 3, ...
              'MarkerEdgeColor', 'b', ...
              'MarkerFaceColor', 'b');
set(handles2, 'MarkerSize', 6, ...
              'MarkerEdgeColor', 'g', ...
              'MarkerFaceColor', 'g');
set(handles3, 'MarkerSize', 9, ...
              'MarkerEdgeColor', 'r', ...
              'MarkerFaceColor', 'r');
set(handles4, 'MarkerSize', 20, 'LineWidth', 4);
axis([0 analysisDuration / timeScale rho0 1e3]);
xlabel(['Time [' timeLabel ']']);
ylabel('Signal to noise ratio');
title([channelName ' clustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/clustered_time_snr_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/clustered_time_snr_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/clustered_time_snr_' num2str(rho0) '_thumbnail.png']);

% plot unclustered trigger frequency vs. time
disp('plotting unclustered trigger frequency vs. time...');
clf;
hold on;
set(gca, 'FontSize', 18);
handles1 = plot(triggersTime, ...
                log2(triggersFrequency), ...
                'bo');
handles2 = plot(triggersTime(triggers10Indices), ...
                log2(triggersFrequency(triggers10Indices)), ...
                'go');
handles3 = plot(triggersTime(triggers20Indices), ...
                log2(triggersFrequency(triggers20Indices)), ...
                'ro');
handles4 = plot(triggersTime(triggersLoudestIndex), ...
                log2(triggersFrequency(triggersLoudestIndex)), ...
                'ko');
set(gca, 'XScale', 'linear');
set(gca, 'YScale', 'linear');
set(handles1, 'MarkerSize', 3, ...
              'MarkerEdgeColor', 'b', ...
              'MarkerFaceColor', 'b');
set(handles2, 'MarkerSize', 6, ...
              'MarkerEdgeColor', 'g', ...
              'MarkerFaceColor', 'g');
set(handles3, 'MarkerSize', 9, ...
              'MarkerEdgeColor', 'r', ...
              'MarkerFaceColor', 'r');
set(handles4, 'MarkerSize', 20, 'LineWidth', 4);
axis([0 analysisDuration / timeScale 5 11]);
set(gca, 'YTick', 5 : 1 : 11);
set(gca, 'YTickLabel', 2.^(5 : 1 : 11));
xlabel(['Time [' timeLabel ']']);
ylabel('Frequency [Hz]');
title([channelName ' unclustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/unclustered_time_frequency_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/unclustered_time_frequency_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/unclustered_time_frequency_' num2str(rho0) '_thumbnail.png']);

% plot clustered trigger frequency vs. time
disp('plotting clustered trigger frequency vs. time...');
clf;
hold on;
set(gca, 'FontSize', 18);
handles1 = plot(clustersTime, ...
                log2(clustersFrequency), ...
                'bo');
handles2 = plot(clustersTime(clusters10Indices), ...
                log2(clustersFrequency(clusters10Indices)), ...
                'go');
handles3 = plot(clustersTime(clusters20Indices), ...
                log2(clustersFrequency(clusters20Indices)), ...
                'ro');
handles4 = plot(clustersTime(clustersLoudestIndex), ...
                log2(clustersFrequency(clustersLoudestIndex)), ...
                'ko');
set(gca, 'XScale', 'linear');
set(gca, 'YScale', 'linear');
set(handles1, 'MarkerSize', 3, ...
              'MarkerEdgeColor', 'b', ...
              'MarkerFaceColor', 'b');
set(handles2, 'MarkerSize', 6, ...
              'MarkerEdgeColor', 'g', ...
              'MarkerFaceColor', 'g');
set(handles3, 'MarkerSize', 9, ...
              'MarkerEdgeColor', 'r', ...
              'MarkerFaceColor', 'r');
set(handles4, 'MarkerSize', 20, 'LineWidth', 4);
axis([0 analysisDuration / timeScale 5 11]);
set(gca, 'YTick', 5 : 1 : 11);
set(gca, 'YTickLabel', 2.^(5 : 1 : 11));
xlabel(['Time [' timeLabel ']']);
ylabel('Frequency [Hz]');
title([channelName ' clustered triggers with SNR > ' num2str(rho0)]);
handle = axes('Position', get(gca, 'Position'));
set(handle, 'FontSize', 18, 'YTick', [], 'XTick', [], ...
            'Color', 'none', 'YAxisLocation', 'right');
ylabel(timeStamp);
print('-dpng', '-r100', [outputDirectory '/clustered_time_frequency_' num2str(rho0) '.png']);
unix(['convert -resize 350x ' ...
      outputDirectory '/clustered_time_frequency_' num2str(rho0) '.png ' ...
      '-strip -depth 8 -colors 256 ' ...
      outputDirectory '/clustered_time_frequency_' num2str(rho0) '_thumbnail.png']);

% close all figures and quit
disp('done.');
exit