% PROCESS2 Prepare segment lists for S4
%
% PROCESS2 prepares single detector and coincident segment lists for
% LIGO's S4 run.

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

% $Id: process2.m,v 1.1 2007/05/15 15:48:31 shourov Exp $

% include segments toolbox in matlab path
path('../matlab', path);

% load concatenated science and injection mode segments
G1_segments = loadsegments('G1_segments.txt');
H1_segments = loadsegments('H1_segments.txt');
H2_segments = loadsegments('H2_segments.txt');
L1_segments = loadsegments('L1_segments.txt');

% apply epoch cut
analysis_segment = [793130413 795679213];
G1_segments = andsegments(G1_segments, analysis_segment);
H1_segments = andsegments(H1_segments, analysis_segment);
H2_segments = andsegments(H2_segments, analysis_segment);
L1_segments = andsegments(L1_segments, analysis_segment);

% produce coincident segment lists
H1H2_segments = andsegments(H1_segments, H2_segments);
H1L1_segments = andsegments(H1_segments, L1_segments);
H2L1_segments = andsegments(H2_segments, L1_segments);
H1H2L1_segments = andsegments(H1_segments, H2_segments, L1_segments);

% cull short segments
minimumDuration = 63;
G1_segments = cullsegments(G1_segments, minimumDuration);
H1_segments = cullsegments(H1_segments, minimumDuration);
H2_segments = cullsegments(H2_segments, minimumDuration);
L1_segments = cullsegments(L1_segments, minimumDuration);
H1H2_segments = cullsegments(H1H2_segments, minimumDuration);
H1L1_segments = cullsegments(H1L1_segments, minimumDuration);
H2L1_segments = cullsegments(H2L1_segments, minimumDuration);
H1H2L1_segments = cullsegments(H1H2L1_segments, minimumDuration);

% split long segments
maximumDuration = 4 * 60 * 60;
minimumOverlap = 16;
G1_segments = splitsegments(G1_segments, maximumDuration, minimumOverlap);
H1_segments = splitsegments(H1_segments, maximumDuration, minimumOverlap);
H2_segments = splitsegments(H2_segments, maximumDuration, minimumOverlap);
L1_segments = splitsegments(L1_segments, maximumDuration, minimumOverlap);
H1H2_segments = splitsegments(H1H2_segments, maximumDuration, minimumOverlap);
H1L1_segments = splitsegments(H1L1_segments, maximumDuration, minimumOverlap);
H2L1_segments = splitsegments(H2L1_segments, maximumDuration, minimumOverlap);
H1H2L1_segments = splitsegments(H1H2L1_segments, maximumDuration, minimumOverlap);

% round to nearest integer seconds
G1_segments = floor(G1_segments);
H1_segments = floor(H1_segments);
H2_segments = floor(H2_segments);
L1_segments = floor(L1_segments);
H1H2_segments = floor(H1H2_segments);
H1L1_segments = floor(H1L1_segments);
H2L1_segments = floor(H2L1_segments);
H1H2L1_segments = floor(H1H2L1_segments);

% remove raw segments
!rm -f G1_segments.txt
!rm -f H1_segments.txt
!rm -f H2_segments.txt
!rm -f L1_segments.txt

% write segments
writesegments(G1_segments, 'G1.txt', 'two');
writesegments(H1_segments, 'H1.txt', 'two');
writesegments(H2_segments, 'H2.txt', 'two');
writesegments(L1_segments, 'L1.txt', 'two');
writesegments(H1H2_segments, 'H1H2.txt', 'two');
writesegments(H1L1_segments, 'H1L1.txt', 'two');
writesegments(H2L1_segments, 'H2L1.txt', 'two');
writesegments(H1H2L1_segments, 'H1H2L1.txt', 'two');

% produce exclusively single and double coincident segment lists
H1h2l1_segments = andsegments(H1_segments, ...
                              notsegments(H2_segments), ...
                              notsegments(L1_segments));
h1H2l1_segments = andsegments(notsegments(H1_segments), ...
                              H2_segments, ...
                              notsegments(L1_segments));
h1h2L1_segments = andsegments(notsegments(H1_segments), ...
                              notsegments(H2_segments), ...
                              L1_segments);
H1H2l1_segments = andsegments(H1_segments, ...
                              H2_segments, ...
                              notsegments(L1_segments));
H1h2L1_segments = andsegments(H1_segments, ...
                              notsegments(H2_segments), ...
                              L1_segments);
h1H2L1_segments = andsegments(notsegments(H1_segments), ...
                              H2_segments, ...
                              L1_segments);

% determine livetimes
analysis_livetime = sum(lensegments(mergesegments(analysis_segment)));
G1_livetime = sum(lensegments(mergesegments(G1_segments)));
H1_livetime = sum(lensegments(mergesegments(H1_segments)));
H2_livetime = sum(lensegments(mergesegments(H2_segments)));
L1_livetime = sum(lensegments(mergesegments(L1_segments)));
H1H2_livetime = sum(lensegments(mergesegments(H1H2_segments)));
H1L1_livetime = sum(lensegments(mergesegments(H1L1_segments)));
H2L1_livetime = sum(lensegments(mergesegments(H2L1_segments)));
H1H2L1_livetime = sum(lensegments(mergesegments(H1H2L1_segments)));
H1h2l1_livetime = sum(lensegments(mergesegments(H1h2l1_segments)));
h1H2l1_livetime = sum(lensegments(mergesegments(h1H2l1_segments)));
h1h2L1_livetime = sum(lensegments(mergesegments(h1h2L1_segments)));
H1H2l1_livetime = sum(lensegments(mergesegments(H1H2l1_segments)));
H1h2L1_livetime = sum(lensegments(mergesegments(H1h2L1_segments)));
h1H2L1_livetime = sum(lensegments(mergesegments(h1H2L1_segments)));

% determine appropriate time unit
if analysis_livetime < 5 * 60,
  timeUnitValue = 1;
  timeUnitName = 's';
elseif analysis_livetime < 5 * 60 * 60,
  timeUnitValue = 60;
  timeUnitName = 'm';
elseif analysis_livetime < 3 * 24 * 60 * 60,
  timeUnitValue= 60 * 60;
  timeUnitName = 'h';
elseif analysis_livetime < 2 * 365.25 * 24 * 60 * 60,
  timeUnitValue = 24 * 60 * 60;
  timeUnitName = 'd';
else,
  timeUnitValue = 365.25 * 24 * 60 * 60;
  timeUnitName = 'y';
end

% convert livetimes to days
analysis_livetime = analysis_livetime / timeUnitValue;
G1_livetime = G1_livetime / timeUnitValue;
H1_livetime = H1_livetime / timeUnitValue;
H2_livetime = H2_livetime / timeUnitValue;
L1_livetime = L1_livetime / timeUnitValue;
H1H2_livetime = H1H2_livetime / timeUnitValue;
H1L1_livetime = H1L1_livetime / timeUnitValue;
H2L1_livetime = H2L1_livetime / timeUnitValue;
H1H2L1_livetime = H1H2L1_livetime / timeUnitValue;
H1h2l1_livetime = H1h2l1_livetime / timeUnitValue;
h1H2l1_livetime = h1H2l1_livetime / timeUnitValue;
h1h2L1_livetime = h1h2L1_livetime / timeUnitValue;
H1H2l1_livetime = H1H2l1_livetime / timeUnitValue;
H1h2L1_livetime = H1h2L1_livetime / timeUnitValue;
h1H2L1_livetime = h1H2L1_livetime / timeUnitValue;

% determine duty cycles
G1_dutycycle = 100 * G1_livetime / analysis_livetime;
H1_dutycycle = 100 * H1_livetime / analysis_livetime;
H2_dutycycle = 100 * H2_livetime / analysis_livetime;
L1_dutycycle = 100 * L1_livetime / analysis_livetime;
H1H2_dutycycle = 100 * H1H2_livetime / analysis_livetime;
H1L1_dutycycle = 100 * H1L1_livetime / analysis_livetime;
H2L1_dutycycle = 100 * H2L1_livetime / analysis_livetime;
H1H2L1_dutycycle = 100 * H1H2L1_livetime / analysis_livetime;
H1h2l1_dutycycle = 100 * H1h2l1_livetime / analysis_livetime;
h1H2l1_dutycycle = 100 * h1H2l1_livetime / analysis_livetime;
h1h2L1_dutycycle = 100 * h1h2L1_livetime / analysis_livetime;
H1H2l1_dutycycle = 100 * H1H2l1_livetime / analysis_livetime;
H1h2L1_dutycycle = 100 * H1h2L1_livetime / analysis_livetime;
h1H2L1_dutycycle = 100 * h1H2L1_livetime / analysis_livetime;

% write livetimes and duty cycles
fid = fopen('livetime.txt', 'w');
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'G1:', G1_livetime, timeUnitName, G1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1:', H1_livetime, timeUnitName, H1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H2:', H2_livetime, timeUnitName, H2_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'L1:', L1_livetime, timeUnitName, L1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1H2:', H1H2_livetime, timeUnitName, H1H2_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1L1:', H1L1_livetime, timeUnitName, H1L1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H2L1:', H2L1_livetime, timeUnitName, H2L1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1H2L1:', H1H2L1_livetime, timeUnitName, H1H2L1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1h2l1:', H1h2l1_livetime, timeUnitName, H1h2l1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'h1H2l1:', h1H2l1_livetime, timeUnitName, h1H2l1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'h1h2L1:', h1h2L1_livetime, timeUnitName, h1h2L1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1H2l1:', H1H2l1_livetime, timeUnitName, H1H2l1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'H1h2L1:', H1h2L1_livetime, timeUnitName, H1h2L1_dutycycle);
fprintf(fid, '%-8s %5.1f%s %5.1f%%\n', ...
        'h1H2L1:', h1H2L1_livetime, timeUnitName, h1H2L1_dutycycle);
fclose(fid);