function MakeFactors_S4V4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MakeFactors_S4V4 - Matlab script to get calibration factors and DQ lists
%                    for all three IFOs for S4 Ligo data. (Version 4)
%                    %% This is not exactly the way how the work was done.
%                       However I put all steps for this work together in
%                       this file. If anyone likes to try to redo, he/she
%                       may need some arrangement and minor changes, but 
%                       without too much difficulty.
%              December 16, 2005               M. Sung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Read the output files from 'ComputeFactors' and put them into a .mat
%%% file all together. - This was done in the grid cluster machine after 
%%% condor jobs were fishished. I just have these lines here to show how
%%% this was done.
% SaveCalibChannel('H1', 'M', 3, 'V4') ;
% SaveCalibChannel('H1', 'M', 3, 'V4') ;
% SaveCalibChannel('H2', 'S', 3, 'V4') ;
% SaveCalibChannel('H2', 'S', 3, 'V4') ;
% SaveCalibChannel('L1', 'M', 3, 'V4') ;
% SaveCalibChannel('L1', 'S', 3, 'V4') ;

%   Calculate the calibration factors for each IFO.
%%   This step is done in local machines after the .mat files resulting
%%   from above were downloaded. - The .mat files from this step will be
%%   avaible in CVS. These files will keep complex values of those
%%   calibration coefficients.
GetFactorsV4('H1','M',3,'V4') ;
GetFactorsV4('H1','S',3,'V4') ;
GetFactorsV4('H2','M',3,'V4') ;
GetFactorsV4('H2','S',3,'V4') ;
GetFactorsV4('L1','M',3,'V4') ;
GetFactorsV4('L1','S',3,'V4') ;

%%% Write out the calibration coefficient files and DQ lists.
write_it_out('H1','M') ;
write_it_out('H1','S') ;
write_it_out('H2','M') ;
write_it_out('H2','S') ;
write_it_out('L1','M') ;
write_it_out('L1','S') ;

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function SaveCalibChannel(ifo, ms, ifcl, runflag)
%% SaveCalibChannel - Read the output text files from 'ComputeFactors' and
%% put them all together in one .mat file.

clear T ASQ DARM EXC fline

pname = 'SaveCalibChannels: ' ;

run = 'S4';
if strcmp(ifo, 'L1')
    fcline = [54.7; 396.7; 1151.5]; 
elseif strcmp(ifo, 'H1')
    fcline = [46.7; 393.1; 1144.3];
elseif strcmp(ifo, 'H2')
    fcline = [54.1; 407.3; 1159.7];
end
cline = [cellstr(int2str(floor(fcline(1)))) cellstr(int2str(floor(fcline(2)))) cellstr(int2str(floor(fcline(3))))];

freq = fcline(ifcl) ; cfreq = char(cline(ifcl)) ;

%% Input files... 
prefix = ['outputs/' run '_' ifo '_' ms '_f'];
files  = [prefix cfreq '_' runflag '-*.dat'];
nfiles = length(dir(files));

if nfiles>0
    T = [ ]; ASQ = [ ]; DARM = [ ]; EXC = [ ];
    for ifile = 0:nfiles-1
        infile = [prefix cfreq '_' runflag '-' num2str(ifile) '.dat'];

        disp([pname 'Reading data from ' infile '...'])
        [Ti, ASQi, DARMi, EXCi] = ReadChannels(infile);
        if Ti ~=-1
            T    = [T; Ti]      ; ASQ  = [ASQ;  ASQi];
            DARM = [DARM; DARMi]; EXC  = [EXC;  EXCi];
        end
    end
else        % for early data without subsegment.
    infile = [prefix char(cline(ifcl)) '_' runflag '.dat'];
    disp([pname 'Reading data from ' infile ' ...'])
    [T, ASQ, DARM, EXC] = ReadChannels(infile);    
end

outfile = ['mout/' run '_' ifo '_' ms '_' runflag '_cl' num2str(ifcl)];
disp([pname 'Writing calibration channels into ' outfile '.mat ...'])
save(outfile, 'freq', 'T', 'ASQ', 'DARM', 'EXC')

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T, ASQ, DARM, EXC] = ReadChannels(filename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pname = 'ReadChannels: ';

%Read factor files for demodulated lines
dline=textread(filename,'','commentstyle','matlab');

nd = length(dline) ;
if nd == 0
    display([pname 'No data for this segment!!! Return.'])
    T = -1.; ASQ = -1.; DARM = -1.; EXC = -1.; 
else
    T    = dline(:,1) ;
    ASQ  = dline(:,8)  + i*dline(:,9);
    DARM = dline(:,10) + i*dline(:,11);  
    EXC  = dline(:,12) + i*dline(:,13);
end

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function GetFactorsV4(ifo, ms, ifcl, runflag)
%%% GetFactors - Calculate and Save alpha*beta in a .mat file

%% Some constants...
pname = 'GetFactorsV4: ';

run = 'S4';
%% Frequency of the calibration line.
if strcmp(ifo,'L1')
    fcline = [54.7 396.7 1151.5] ;
elseif strcmp(ifo,'H1')
    fcline = [46.7 393.1 1144.3] ;
elseif strcmp(ifo,'H2')
    fcline = [54.1 407.3 1159.7] ;    
end

clfreq = fcline(ifcl) ;
%% Get the open loop gain and digital filter funciton for the frequency.
[G0, D0, bound] = GetG0D0(ifo, clfreq, 0);

disp([pname 'Getting time dependent calibration factors from S4 ' ifo ' data'])

%% Input file... 
infile = [ 'mout/' run '_' ifo '_' ms '_' runflag '_cl' num2str(ifcl) '.mat'];
disp([pname 'Reading data from ' infile '...'])
load(infile)

ntotal = length(T);
display([pname '  Total ' num2str(ntotal) ' entries'])

%% Calculate alpha, betas from demodulated lines using these values
disp([pname 'Calculating alpha and beta...'])
alpha = -(D0/G0)*(ASQ./DARM) ; beta = (1/D0)*((DARM-EXC)./ASQ);

%% Calculate G0 which can vanish the imaginary part of alpha*beta.
xab = -(DARM-EXC)./DARM ; mxab = mean(xab) ;
cxab = real(mxab)/imag(mxab); %cG0 = real(G0)/imag(G0) ;
Gi = imag(mxab) ; Gr = Gi*cxab ; Gest = complex(Gr,Gi) ;
display([pname 'G = ' num2str(Gest) ])

%% find nans, zero them out
inan = find(isnan(alpha) | isnan(beta) | isinf(alpha) | isinf(beta)) ;
alpha(inan) = 0 ;      beta(inan) = 0;

%% alpha*beta
albt = alpha.*beta ;

%Find unphysical values and set them to zero....
nlow = 0                        ; nhigh = 0;
ilow = find(real(albt)<bound.lo); ihigh = find(real(albt)>bound.hi);
if isempty(ilow) ~= 1
    alpha(ilow) = 0 ;   beta(ilow) = 0 ; albt(ilow) = 0 ;
    nlow = length(ilow);
    disp([pname 'Found ' num2str(nlow) ...
    ' alpha*beta values too low for f = ' num2str(fcline(ifcl)) ' hz.'])
end
if isempty(ihigh) ~= 1
    alpha(ihigh) = 0 ;  beta(ihigh) = 0; albt(ihigh) = 0; 
    nhigh = length(ihigh);
    disp([pname 'Found ' num2str(nhigh) ...
    ' alpha*beta values too high for f = ' num2str(fcline(ifcl)) 'hz.'])
end
Bad.nan = T(inan) ; Bad.low = T(ilow) ; Bad.high = T(ihigh);

ibad    = unique([inan ; ilow ; ihigh]); nout    = length(ibad);
Bad.all = T(ibad) ;       
%ngood   = length(T) - nout ;

if strcmp(ms,'16hz')
    clear T albt
end

%%% Save some values in a .mat file
fsave=['dat/' run '_' ifo '_' ms '_' runflag '_cl' num2str(ifcl) '.mat'];
save( fsave,'T','albt', 'alpha', 'beta', 'clfreq', 'Bad')

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [G0, D0, bound] = GetG0D0(ifo, freq, lplot)
%% GetG0D0 - D0, G0 values for S4

pname = 'GetG0D0: ';
disp([pname 'Getting Open Loop Gain and Digital Feedback Filter '])
disp([pname '  at frequency = ' num2str(freq) 'Hz.'])
if strcmp(ifo,'L1')
    fid = fopen('calib/L1/V4try2/L-CAL_REF_OLOOP_GAIN_S4-V4_L1-795727983-8.txt') ;
elseif strcmp(ifo,'H1')
%    fid = fopen('calib/H1/V4try2/H-CAL_REF_OLOOP_GAIN_S4_V4_H1-793099715-8.txt.prov');
    fid = fopen('calib/H1/V4/H-CAL_REF_OLOOP_GAIN_S4_V4_H1-793099715-8.txt');
elseif strcmp(ifo,'H2')
%    fid = fopen('calib/H2/V2/S04-H2-CAL-OLOOP_GAIN-793064176.txt');
    fid = fopen('calib/H2/V4/H-CAL_REF_OLOOP_GAIN_S4_V4_H2-793064176-8.txt');
end
tmp = textscan(fid,'%f%f%f','commentStyle','%');
G   = cell2struct(tmp,{'freq','amp','ang'},2);
clear Gcell
fclose(fid);

%% Get the boundary for the physical region for alpha*bet
bound = GetBound(G, lplot);

if strcmp(ifo,'L1')
    fid = fopen('calib/L1/V4try2/L-CAL_REF_DIGFLT_AS_Q_S4-V4_L1-795727983-8.txt');
    tmp = textscan(fid,'%f%f%f','commentStyle','%');
    D = cell2struct(tmp,{'freq','amp','ang'},2);
    clear Gcell
    fclose(fid);
    
    G0 = interp1(G.freq(2:end),G.amp(2:end).*exp(i*G.ang(2:end)),freq);
    D0 = interp1(D.freq(2:end),D.amp(2:end).*exp(i*D.ang(2:end)),freq);
    
elseif strcmp(ifo,'H1')
%    fid = fopen('calib/H1/S04-H1-CAL-ACTUATION-793099715.txt');
%    tmp = textscan(fid,'%f%f%f');
%    A   = cell2struct(tmp,{'freq','amp','ang'},2);
%    fclose(fid);

%    fid = fopen('calib/H1/S04-H1-ASQCAL-CAV_GAIN-793099715.txt');
%    fid = fopen('calib/H1/V2/S04-H1-ASQCAL-DIGFLT-793099715.txt');
%    fid = fopen('calib/H1/V4try2/H-CAL_REF_DIGFLT_AS_Q_S4_V4_H1-793099715-8.txt.prov');
    fid = fopen('calib/H1/V4/H-CAL_REF_DIGFLT_AS_Q_S4_V4_H1-793099715-8.txt');
    tmp = textscan(fid,'%f%f%f');
    D   = cell2struct(tmp,{'freq','amp','ang'},2);
    fclose(fid);

    G0 = interp1(G.freq(2:end),G.amp(2:end).*exp(i*G.ang(2:end)),freq);
    D0 = interp1(D.freq(2:end),D.amp(2:end).*exp(i*D.ang(2:end)),freq);
%    A0 = (1/4e3)*interp1(A.freq(2:end),A.amp(2:end).*exp(i*A.ang(2:end)),freq);
%    C0 = interp1(C.freq(2:end),C.amp(2:end).*exp(i*C.ang(2:end)),freq);
%    D0 = G0./(A0.*C0);
elseif strcmp(ifo,'H2')
%%    fid = fopen('H2/S04-H2-CAL-ACTUATION-793064176.txt');
%    fid = fopen('CALS4V2/H2/S04-H2-CAL-ACTUATION-793064176.txt');
%    tmp = textscan(fid,'%f%f%f');
%    A   = cell2struct(tmp,{'freq','amp','ang'},2);
%    fclose(fid);

%    fid = fopen('CALS4V2/H2/S04-H2-ASQCAL-CAV_GAIN-793064176.txt');
%    fid = fopen('calib/H2/V2/S04-H2-ASQCAL-DIGFLT-793064176.txt');
    fid = fopen('calib/H2/V4/H-CAL_REF_DIGFLT_AS_Q_S4_V4_H2-793064176-8.txt');
    tmp = textscan(fid,'%f%f%f');
%    C   = cell2struct(tmp,{'freq','amp','ang'},2);
    D   = cell2struct(tmp,{'freq','amp','ang'},2);
    fclose(fid);
    
    G0 = interp1(G.freq(2:end),G.amp(2:end).*exp(i*G.ang(2:end)),freq);
    D0 = interp1(D.freq(2:end),D.amp(2:end).*exp(i*D.ang(2:end)),freq);
%    A0 = (1/4e2)*interp1(A.freq(2:end),A.amp(2:end).*exp(i*A.ang(2:end)),freq);
%    C0 = interp1(C.freq(2:end),C.amp(2:end).*exp(i*C.ang(2:end)),freq);
%    D0 = G0./(A0.*C0);
end

disp(pname)
disp([pname ' G(' num2str(freq) ' Hz) = ' num2str(G0)])
disp([pname ' D(' num2str(freq) ' Hz) = ' num2str(D0)])
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bound = GetBound(G, lplot)
%%% GetBound - Get the bound for (alpha*beta) by using the open loop gain
%%%            plots.

pname = 'GetBound: ';

ifq = 20; ffq = 700;

dang = diff(G.ang);

ifreq1 = find(dang<-5);
ifreq_lo = find(abs(G.freq(ifreq1)-50)<30);
if length(ifreq_lo) ~= 1
    display([pname 'Two solutions for low frequency - confusing! Quit.'])
    return
end

freq_lo = G.freq(ifreq1(ifreq_lo));
bound.lo = 1./G.amp(ifreq1(ifreq_lo));

ifreq2 = find(dang>0.01);
ifreq_hi = find(abs(G.freq(ifreq2)-350)<100);
if length(ifreq_hi) ~= 1
%    G.freq(ifreq2(ifreq_hi))
%    display([pname 'Two solutions for high frequency - Confusing! Quit.'])
    display([pname 'Two solutions for high frequency - Use the smaller freq.'])
%    return
end

freq_hi = G.freq(ifreq2(ifreq_hi(1)));
bound.hi = 1./G.amp(ifreq2(ifreq_hi(1)));

display([pname 'Boundary for alpha*beta: ' num2str(bound.lo) '   ' num2str(bound.hi)])

if lplot == 1;
    clf
    fqbin = [ifq:.125:ffq]*8;
    subplot 211; semilogy(G.freq(fqbin), G.amp(fqbin),'LineWidth',2)
    xlabel('Frequency [Hz]'); ylabel('|G|')
    ax1 = axis; axis([ifq ffq ax1(3) ax1(4)])
    line([ifq ffq], [1 1],'color','g')
    line([freq_lo freq_lo], [ax1(3) ax1(4)],'color', 'r')
    line([ifq 150], [1./bound.lo 1./bound.lo],'color', 'r')
    text(freq_lo+10, 1/bound.lo(1)-.3, ['1/' num2str(bound.lo(1))] )
    text(freq_hi-70, 1/bound.hi(1)-.06, ['1/' num2str(bound.hi(1))] )
    line([freq_hi freq_hi], [ax1(3) ax1(4)],'color', 'r')
    line([500 ffq], [1./bound.hi 1./bound.hi],'color', 'r')

    subplot 212; 
    plot(G.freq(fqbin), G.ang(fqbin),G.freq(fqbin),dang(fqbin),'r','LineWidth',2)
    ax2 = axis; axis([ifq ffq ax2(3) ax2(4)])
    xlabel('Frequency[Hz]'); ylabel('\theta(G)')
    text(freq_lo+10,-pi-1.5, ['f_{lo} = ' num2str(freq_lo) ' Hz'])
    text(freq_hi-100,-pi-1.5, ['f_{hi} = ' num2str(freq_hi) ' Hz'])

    subplot 211 ;

    print('-dpdf','bound.pdf')
    print('-dpng','bound.png')
end

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function endout = write_it_out(ifo, ms)
% write_it_out - Write out alpha, alpha*beta and DQ list in text files

clear
mname = 'write_it_out: ';
if strcmp(ms, 'M')
    ctbin = '060' ;
elseif strcmp(ms,'S')
    ctbin = '001' ;
    dqseg = []    ; dtdq   = [] ; Tdq1s    = [] ;
end

% Read the Segment list file.
fseg = ['S4' ifo 'v00_segs.txt'];
[iseg, tiseg, tfseg, dtseg] = textread(fseg,'%d %d %d %d\n','commentstyle','shell');
nseg = length(iseg); Ttotal = sum(dtseg);

matfile = ['dat/S4_' ifo '_' ms '_conlog_cl3.mat'] ; 
load(matfile,'T','albt','alpha','beta'); 

Tgps = T             ; T0  = 793128400              ; Ts  = T - T0 ;
cti  = num2str(T(1)) ; cdt = num2str(T(end) - T(1)) ;

toutfile = [ifo(1) '-' ifo '_CAL_FAC_S4_V4_' ctbin '-' cti '-' cdt '.txt'];
display([mname 'Writing on ' toutfile '...'])
fid = fopen(toutfile,'W');
fprintf(fid,['%% DeltaT = ' ctbin ' second\n']);
fprintf(fid,'%% $Name  = S4-V4-0-0\n');
fprintf(fid,'%% GPSTime            alpha              gamma\n');
tout = [Tgps'; abs(real(alpha))'; real(albt)'];
fprintf(fid,' %9d %18.9f %18.9f\n',tout);
fclose(fid);

igd = find(albt) ;
figure(1); clf ; 
subplot(211) 
plot(Ts(igd),real(albt(igd)),'.') ; %ylim([.85 1.15])%hline(1,'w') ; 
title(['S4V4: Calibration Factors \gamma - ' ifo ': All conlog segments']) ;
xlabel(['T_{gps} (T_0=' num2str(T0) ')']); ylabel('Re(\gamma)')
subplot(212)
plot(Ts(igd),real(alpha(igd)),'.r', Ts(igd),real(beta(igd)),'.g')
%ylim([-1.2 -.8]) ; %hline(-1,'w') ; 
xlabel(['T_{gps} (T_0=' num2str(T0) ')']); 
legend('Re(\alpha)','Re(\beta)','Location','SouthEast')
title('Re(\alpha), Re(\beta)') ;
fplot = [ifo '_V4_' ctbin '_conlog'];
print('-djpeg',fplot)
    
% Open output DQ list files.
fDQ = ['S4V4_' ifo '_DQ_' ctbin '.txt'] ; fiddq = fopen(fDQ,'W');
% Looping on segments.
Tdq = 0 ; Tgd = [] ; Tout = [];  
albtgd = [] ; alphagd = [] ; betagd = [] ; ab_out = []; al_out = [];

abmax = 1.1 ; abmin = 0.9 ;
if strcmp(ifo,'H2')
    abmax = 1.07 ;
end
abmax1 = 1.5 ; abmin1 = 0.5 ;
for is = 1:nseg
    if strcmp(ms,'S')                 % For 1 second calibration factors
        iT1  = find(Tgps>=tiseg(is) & Tgps<tfseg(is)) ;
        if isempty(iT1)                  % No factor for the segment.
            Tdq1s = [Tdq1s; tiseg(is)] ; dtdq = [dtdq; dtseg(is)] ;
            dqseg = [dqseg; iseg(is)]  ;
        else
            Tseg = Tgps(iT1)           ; Tend = Tseg(end) + 1     ;
            if Tseg(1) > tiseg(is)      % Beginning of the segment.
                Tdq1s = [Tdq1s ; tiseg(is)] ; dqseg = [dqseg; iseg(is)] ;
                dtDQ  = Tseg(1) - tiseg(is) ; dtdq  = [dtdq; dtDQ]      ;
            end

            % Middle of the segment.
            ab1seg  = albt(iT1)                ; 
            al1seg  = alpha(iT1)               ; bt1seg = beta(iT1)                   ;
            izr1    = find(~real(ab1seg))      ; inzr1  = find(real(ab1seg))          ;
            Tgd     = [Tgd; Tseg(inzr1) ]      ; albtgd = [albtgd; ab1seg(inzr1)]     ; 
            alphagd = [alphagd; al1seg(inzr1)] ; betagd = [betagd; bt1seg(inzr1)]     ;

            Tbad    = Tseg(izr1)               ; idq    = find(tfseg(is)-Tbad>15)     ;
            Tdq1s   = [Tdq1s ; Tbad(idq)]      ; dtdq   = [dtdq; ones(length(idq),1)] ;
            dqseg   = [dqseg ; ones(length(izr1),1)*iseg(is)] ;

            % Outliers, especically near the end of science modes.
            clear iout
            iout = find(real(ab1seg)>abmax1 | real(ab1seg)<abmin1) ;
            nout = length(iout) ;
            for io = 1:nout
                if tfseg(is) - Tseg(iout(io)) < 16
                    tout1 = Tseg(iout(io)) ;
                    fprintf(1,'%4d %9d %9d %9d\n', iseg(is),tout1,ab1seg(iout(io)),al1seg(iout(io)));
                    Tout = [Tout; tout1] ; 
                    ab_out = [ab_out; ab1seg(iout(io))]; al_out = [al_out; al1seg(iout(io))];
                    Tdq1s = [Tdq1s ; tout1] ; dqseg = [dqseg; iseg(is)] ;
                    dtdq  = [dtdq; 1 ]     ;
                else
                    fprintf(1,'%4d %10d %10d %s\n', iseg(is),Tseg(iout(io)),ab1seg(iout(io)),'o');
                end
            end

            if Tend < tfseg(is)         % End of the segment.
                Tdq1s = [Tdq1s ; Tend]      ; dqseg = [dqseg; iseg(is)] ;
                dtDQ  = tfseg(is) - Tend    ; dtdq  = [dtdq; dtDQ ]     ;
            end
        end
    elseif strcmp(ms, 'M')             % For 60 second calibration factors.
        iT60 = find(Tgps>tiseg(is)-60 & Tgps<tfseg(is)) ;
        nT = length(iT60) ;
        if isempty(iT60)                 % No factor for the segment.
            fprintf(fiddq,'%4d %10d %10d %7d %s\n', iseg(is),tiseg(is),tfseg(is),dtseg(is),'s');
            Tdq = Tdq + dtseg(is);
        else
            Tseg    = Tgps(iT60)                ; Tend    = Tseg(nT)+60               ;
            ab60seg = albt(iT60)                ; Rab     = real(ab60seg)             ;
            al60seg = alpha(iT60)               ; bt60seg = beta(iT60)                ;
            inzr60  = find(Rab)                 ; 
            Tgd     = [Tgd    ; Tseg(inzr60) ]  ; albtgd  = [albtgd ; ab60seg(inzr60)];
            alphagd = [alphagd; al60seg(inzr60)]; betagd  = [betagd ; bt60seg(inzr60)]; 

            if Tseg(1) > tiseg(is)     % Beginning of the segment.
                dtDQ = Tseg(1) - tiseg(is) ;
                fprintf(fiddq,'%4d %10d %10d %7d %s\n',iseg(is),tiseg(is),Tseg(1),dtDQ,'b');
                Tdq = Tdq + dtDQ60 ;
            end

            % Outliers.
            clear iout
            iout = find(Rab>abmax | Rab<abmin) ; nout = length(iout) ;
            if nout == 1 & iout == nT
                Tend = Tseg(iout) ;
                fprintf(1,'%4d %9d %9d %9d %s\n', iseg(is),Tend,ab60seg(iout),al60seg(iout),'e');
                Tout = [Tout; Tend] ; 
                ab_out = [ab_out; ab60seg(nT)]; al_out = [al_out; al60seg(nT)];
            elseif nout > 1 & iout(end) == nT
                tagout(1:nout-1) = 'm' ; tagout(nout) = 'e' ;
                dout = diff(iout) ; idout = nout - 1 ;
                while idout> 0 & dout(idout) == 1
                    iend  = iout(idout) ; Tend = Tseg(iend) ; tagout(idout) = 'e' ;
                    idout = idout - 1   ;
                end
                for io = 1:nout
                    fprintf(1,'%4d %9d %9d %9d %s\n', iseg(is),Tseg(iout(io)),ab60seg(iout(io)),al60seg(iout(io)),tagout(io));
                    if strcmp(tagout(io),'e')
                        Tout = [Tout; Tseg(iout(io))] ; 
                        ab_out = [ab_out; ab60seg(iout(io))]; al_out = [al_out; al60seg(iout(io))];
                    end
                end
            else
                for io = 1:nout
                    fprintf(1,'%4d %9d %9d %9d %s\n', iseg(is),Tseg(iout(io)),ab60seg(iout(io)),al60seg(iout(io)),'m');
                    if iout(io) == 1
                        TiDQ = tiseg(is)      ; TfDQ = Tseg(2)    ; cflag = 'i' ; 
                    else
                        TiDQ = Tseg(iout(io)) ; TfDQ = TiDQ + 60  ; cflag = 'm' ;
                    end
                    dtDQ = TfDQ - TiDQ        ; Tdq  = Tdq + dtDQ ;
                    fprintf(fiddq,'%4d %10d %10d %7d %s\n',iseg(is),TiDQ,TfDQ,dtDQ,cflag);
                end
            end

            if Tend < tfseg(is)        % End of the segment.
                dtDQ = tfseg(is) - Tend ;
                fprintf(fiddq,'%4d %10d %10d %7d %s\n',iseg(is),Tend,tfseg(is),dtDQ,'e');
                Tdq = Tdq + dtDQ ;
            end
        end
    end
end

figure(2); clf
if strcmp(ms, 'S')
    cdq = FindSequence(Tdq1s, dtdq, dqseg) ;
    for iz1 = 1:length(cdq.times)
        fprintf(fiddq,'%4d %10d %10d %7d\n', cdq.ext(iz1),cdq.times(iz1),...
                                  cdq.times(iz1)+cdq.dts(iz1),cdq.dts(iz1));        
    end
    Tdq = sum(cdq.dts) ; 
    % Making plots.
    ylg = 2 ; yla = 2 ;
    if strcmp(ifo,'L1')
        rabhmn =  0.93   ; rabbsz = 0.0025  ; rabhmx = 1.12   ;
        iabhmn = -0.06   ; iabbsz = 0.0015  ; iabhmx = 0.06   ;
        ralhmn =  0.4    ; ralbsz = 0.01    ; ralhmx = 1.3    ;
        ialhmn = -0.05   ; ialbsz = 0.002   ; ialhmx = 0.05   ;
        rbthmn =  0.8    ; rbtbsz = 0.01    ; rbthmx = 2      ;
        ibthmn = -0.0025 ; ibtbsz = 0.00005 ; ibthmx = 0.0016 ;
    elseif strcmp(ifo,'H1')
        rabhmn =  0.85   ; rabbsz = 0.004   ; rabhmx = 1.12    ;
        iabhmn = -0.05   ; iabbsz = 0.002   ; iabhmx = 0.05    ;
        ralhmn =  0.85   ; ralbsz = 0.005   ; ralhmx = 1.1     ;
        ialhmn = -0.05   ; ialbsz = 0.002   ; ialhmx = 0.05    ;
        rbthmn =  0.9    ; rbtbsz = 0.0025  ; rbthmx = 1.1     ;
        ibthmn = -0.0009 ; ibtbsz = 0.00001 ; ibthmx = -0.0004 ;
    elseif strcmp(ifo,'H2')
        rabhmn =  0.85   ; rabbsz = 0.004   ; rabhmx =  1.12   ;
        iabhmn = -0.1    ; iabbsz = 0.005   ; iabhmx =  0.1    ;
        ralhmn =  0.92   ; ralbsz = 0.005   ; ralhmx =  1.2    ;
        ialhmn = -0.12   ; ialbsz = 0.005   ; ialhmx =  0.09   ;
        rbthmn =  0.83   ; rbtbsz = 0.0025  ; rbthmx =  1.04   ;
        ibthmn = -0.0014 ; ibtbsz = 0.00002 ; ibthmx = -0.0006 ;
    end
    subplot(211);
    plot(Tgd,real(albtgd),'.', Tout,real(ab_out),'.r')
    xlabel('S4 GPS Time') ; ylim([0, ylg]); ylabel('Re(\gamma)')
    title([ifo ': CAL-FAC-S4-V4-' ctbin])
    subplot(212); plot(Tgd,abs(real(alphagd)),'.', Tout,abs(real(al_out)),'.r')
    xlabel('S4 GPS Time') ; ylim([0, yla]); ylabel('|Re(\alpha)|');  shg
elseif strcmp(ms, 'M')
    % Making plots.
    ylg = 1.8 ; yla = 1.3;
    if strcmp(ifo,'L1')
        rabhmn =  0.93   ; rabbsz =.003  ; rabhmx = 1.13  ;
        iabhmn = -0.02   ; iabbsz =.0005 ; iabhmx =  .02  ;
        ralhmn =  0.4    ; ralbsz =.02   ; ralhmx = 1.3   ;
        ialhmn = -0.015  ; ialbsz =.0005 ; ialhmx = .015  ;
        rbthmn =  0.8    ; rbtbsz =.01   ; rbthmx = 2     ;
        ibthmn = -0.0025 ; ibtbsz =.0001 ; ibthmx = .0025 ;
    elseif strcmp(ifo,'H1')
        rabhmn =  0.94   ; rabbsz =.001   ; rabhmx = 1.03  ;
        iabhmn = -0.01   ; iabbsz =.00025 ; iabhmx =  .005 ;
        ralhmn =  0.9    ; ralbsz =.002   ; ralhmx = 1.1   ;
        ialhmn = -0.002  ; ialbsz =.00025 ; ialhmx = .01   ;
        rbthmn =  0.9    ; rbtbsz =.0025  ; rbthmx = 1.1   ;
        ibthmn = -0.0008 ; ibtbsz =.00001 ; ibthmx = -.0003;
    elseif strcmp(ifo,'H2')
        rabhmn =  0.95   ; rabbsz =.001   ; rabhmx = 1.05   ;
        iabhmn = -0.013  ; iabbsz =.0005  ; iabhmx =  .015  ;
        ralhmn =   .98   ; ralbsz =.002   ; ralhmx = 1.15   ;
        ialhmn = -0.02   ; ialbsz =.0005  ; ialhmx =  .015  ;
        rbthmn =  0.8    ; rbtbsz =.001   ; rbthmx = 1.1    ;
        ibthmn = -0.0013 ; ibtbsz =.00001 ; ibthmx = -.0006 ;
    end
    subplot(211);
    plot(Tgd,real(albtgd),'.', Tout,real(ab_out),'.r')
    xlabel('S4 GPS Time') ; ylim([0, ylg]); ylabel('Re(\gamma)')
    title([ifo ': CAL-FAC-S4-V4-' ctbin])
    subplot(212); plot(Tgd,abs(real(alphagd)),'.', Tout,abs(real(al_out)),'.r')
    xlabel('S4 GPS Time') ; ylim([0, yla]); ylabel('|Re(\alpha)|');  shg
end
fplot = [ifo '_V4_' ctbin];
print('-djpeg',fplot)

% Some calculation.
Rdq = Tdq/Ttotal ; 
disp([mname 'Ttotal = ' num2str(Ttotal) 's, dTdq  = ' num2str(Tdq) 's'])
disp([mname 'Rdq   = ' num2str(Rdq)])

figure(3); clf       % Histograms
subplot(321)
rabin = rabhmn:rabbsz:rabhmx ;
hfitg(real(albtgd), rabin, rabhmn, rabhmx) ;
xlim([rabhmn rabhmx]) ; xlabel('Re(\gamma)')
subplot(322)
iabbin = iabhmn:iabbsz:iabhmx ;
hfitg(imag(albtgd), iabbin, iabhmn, iabhmx) ;
xlim([iabhmn iabhmx]) ; xlabel('Im(\gamma)')
subplot(323)
ralbin = ralhmn:ralbsz:ralhmx ;
hfitg(abs(real(alphagd)), ralbin, ralhmn, ralhmx) ;
xlim([ralhmn ralhmx]) ; xlabel('|Re(\alpha)|')
subplot(324)
ialbin = ialhmn:ialbsz:ialhmx ;
hfitg(imag(alphagd), ialbin, ialhmn, ialhmx) ;
xlim([ialhmn ialhmx]) ; xlabel('Im(\alpha)')
subplot(325)
rbtbin = rbthmn:rbtbsz:rbthmx ;
hfitg(abs(real(betagd)), rbtbin, rbthmn, rbthmx) ;
xlim([rbthmn rbthmx]) ; xlabel('|Re(\beta)|')
subplot(326)
ibtbin = ibthmn:ibtbsz:ibthmx ;
hfitg(imag(betagd), ibtbin, ibthmn, ibthmx) ;
xlim([ibthmn ibthmx]) ; xlabel('Im(\beta)')
shg
fhplot = [ifo '_V4_' ctbin '_hist'];
print('-djpeg',fhplot)

% Checkout some bad section with interests.
ibad = [] ;
if strcmp(ifo,'H1') & strcmp(ms,'S')
    ibad = find(Tgps>793521500 & Tgps<793522150) ;   % Near the end of segment 22
elseif strcmp(ifo,'H2') & strcmp(ms,'S')
%    ibad = find(Tgps>793130000 & Tgps<793133500) ;   % segment 1
%    ibad = find(Tgps>793132000 & Tgps<793135300) ;   % segment 2
    ibad = find(Tgps>795643000 & Tgps<795647539) ;   % Near the end of segment 192
end
if ~isempty(ibad)
    figure(4) ; clf
    Tbad = Tgps(ibad) ; albtbd = albt(ibad) ; albad = alpha(ibad) ; btbad = beta(ibad);
    ibdnz = find(albtbd) ;
    Tbadnz  = Tbad(ibdnz)  ; albtbdnz = albtbd(ibdnz) ; 
    albadnz = albad(ibdnz) ; btbadnz  = btbad(ibdnz)  ;
    tmx = (Tbadnz(end) - Tbadnz(1) + 25)/60 ; 
    subplot(321) ; 
    plot((Tbadnz-Tbadnz(1)+1)/60,real(albtbdnz)); 
    xlabel(['Time(min): T0 = ' num2str(Tbadnz(1))]) ; ylabel('Re(\gamma)')
    xlim([0 tmx])
    subplot(322) ; 
    plot((Tbadnz-Tbadnz(1)+1)/60,imag(albtbdnz)); 
    xlabel(['Time(min): T0 = ' num2str(Tbadnz(1))]) ; ylabel('Im(\gamma)')
    xlim([0 tmx])
    subplot(323) ; 
    plot((Tbadnz-Tbadnz(1)+1)/60,real(albadnz)); 
    xlabel(['Time(min): T0 = ' num2str(Tbadnz(1))]) ; ylabel('Re(\alpha)')
    xlim([0 tmx])
    subplot(324) ; 
    plot((Tbadnz-Tbadnz(1)+1)/60,imag(albadnz)); 
    xlabel(['Time(min): T0 = ' num2str(Tbadnz(1))]) ; ylabel('Im(\alpha)')
    xlim([0 tmx])
    subplot(325) ; 
    plot((Tbadnz-Tbadnz(1)+1)/60,real(btbadnz)); 
    xlabel(['Time(min): T0 = ' num2str(Tbadnz(1))]) ; ylabel('Re(\beta)')
    xlim([0 tmx])
    subplot(326) ; 
    plot((Tbadnz-Tbadnz(1)+1)/60,imag(btbadnz)); 
    xlim([0 tmx])
    xlabel(['Time(min): T0 = ' num2str(Tbadnz(1))]) ; ylabel('Im(\beta)')
    print('-djpeg',['check_' ifo '_bad'])
end

fclose(fiddq) ;
endout = 'Nicely Done !!!'

return

function compact = FindSequence(times, dts, extra)

mname = 'FindSequence: ' ;

org.times = times ;
org.dts = dts ;
org.ext = extra ;

ntime = length(times) + 1 ;
while length(org.times)<ntime
    ntime = length(org.times) ;
    org = CombineSequence(org) ;
end

compact = org ;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function compact = CombineSequence(org)

tfs = org.times + org.dts ;

int = 0 ; lskip = 0 ;
for its = 1: length(org.times)
    if lskip == 0
        if its ~= length(org.times) & tfs(its) == org.times(its+1) 
            int = int + 1 ;
            lskip = 1 ;
            newts(int) = org.times(its) ;
            newdt(int) = org.dts(its) + org.dts(its+1) ;
            newex(int) = org.ext(its) ;
        else
            int = int + 1 ;
            newts(int) = org.times(its) ;
            newdt(int) = org.dts(its) ;
            newex(int) = org.ext(its) ;
        end
    elseif lskip == 1
        lskip = 0 ;
    end
end

compact.times = newts ;
compact.dts   = newdt ;
compact.ext   = newex ;

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [parout,chisq]=hfitg(x,n,miny,maxy)
global ndf;

hold off
if (nargin == 2)
  miny=min(min(x));
  maxy=max(max(x));
end
[ny,nx]=hpl(x,n,miny,maxy);
%par=[1,1,1000];
ind=find (x >= miny & x <= maxy);
par(1)=median(x(ind));
par(2)=std(x(ind));
par(3)=max(ny);
[parout,chisq]=chisq_min(par,nx,ny);
hold on
plot(nx,parout(3)*g(nx,parout(1),parout(2)),'r')
hold off
v=axis;

mean_value=['Mean: ',num2str(parout(1),3), '\pm',num2str(parout(2),2)];
%std_value=['\sigma: ',num2str(parout(2),2)];
max_value=['Max: ',num2str(parout(3),3)];
chi2_value=['\chi^2/ndf: ',num2str(chisq/(length(x)-3),2)];

%xper=.01;yper=.9;
xper=.55;yper=.9;
text(xper,yper, mean_value,'FontSize',11,'Units','normalized')
yper=yper-.1;
%text(xper,yper, std_value, 'FontSize',11,'Units','normalized')
%xper=.55;yper=.9;
text(xper,yper, max_value, 'FontSize',11,'Units','normalized')
yper=yper-.1;
text(xper,yper, chi2_value,'FontSize',11,'Units','normalized')

%xline=v(1)*(xper-.1)+v(2)*(.9-xper);
%yline=v(3)*(yper-.05)+v(4)*(.95-yper);
%line([xline,xline],[v(4),yline])
%line([xline,v(2)],[yline,yline])