% This script gathers together the noise info of each channel over all the
% the test measurements.

%%% add paths for SRF analysis tools
addpath /home/hannon/SRF_results_UMBC/SRC/Programs
addpath /home/hannon/SRF_results_UMBC/SRC/Functions

temp=input('Enter temperature name (161, 155, 149, or all) : ','s');
%%%%%%%%%%%
switch temp
   case '161'
      testnums=[ 21403  21406  21408  21409  21414  21415  21416  21418 ];
      gain    =[ 'A'    'B'    'O'    'A'    'B'    'O'    'A'    'O'];
   case '155'
      testnums=[ 21479  21480  21481  21484 ];
      gain    =[ 'B'    'O'    'A'    'B'];
   case '149'
      testnums=[ 21616  21619  21620  21621  21622  21625 ];
      gain    =[ 'A'    'A'    'B'    'O'    'A'    'B'];
   case 'all'
      testnums=[ 21403  21406  21408  21409  21414  21415  21416  21418 ...
                 21479  21480  21481  21484 ...
                 21616  21619  21620  21621  21622  21625 ];
      gain    =[ 'A'    'B'    'O'    'A'    'B'    'O'    'A'    'O' ...
                 'B'    'O'    'A'    'B' ...
                 'A'    'A'    'B'    'O'    'A'    'B'];
   otherwise
   error('invalid temperature name')
end
%%%%%%%%%%%


modulex=input('Enter module name (without the M) : ','s');
%%%%%%%%%%%%%
switch modulex
   case '12'
      module='12';
      nchan=130;
      sdir='srf';
      noiselev=2.5E-2;
      indbad=[118  119  82  83  9  10];
   case '11'
      module='11';
      nchan=144;
      sdir='srf';
      noiselev=1.5E-2;
      indbad=[37  141:143];
   case '10'
      module='10';
      nchan=167;
      sdir='srf';
      noiselev=2.2E-2;
      indsum=117:142;
%      indbad=[29  30  37  38];
      indbad=[];
   case '9'
      module='9';
      nchan=167;
      sdir='srf';
      noiselev=2E-2;
      indbad=[115];
   case '8'
      module='8';
      nchan=161;
      sdir='srf';
      noiselev=1.2E-2;
      indbad=[];
   case '7'
      module='7';
      nchan=167;
      sdir='srf';
      noiselev=6E-3;
      indbad=[];
   case '6'
      module='6';
      nchan=167;
      sdir='srf';
      noiselev=7E-3;
      indbad=[];
   case '5'
      module='5';
      nchan=159;
      sdir='srf';
      noiselev=1E-2;
      indbad=[];
   case '4d'
      module='4d';
      nchan=106;
      sdir='srf';
      noiselev=3E-3;
      indbad=[14  15  70:82  86:88];
   case '4d_nozf'
      module='4d';
      nchan=106;
      sdir='srf_nozf';
      noiselev=3E-3;
      indbad=[14  15  70:82  86:88];
   case '4c'
      module='4c';
      nchan=94;
      sdir='srf';
      noiselev=4E-3;
      indbad=[];
   case '4c_nozf'
      module='4c';
      nchan=94;
      sdir='srf_nozf';
      noiselev=4E-3;
      indbad=[];
   case '3'
      module='3';
      nchan=192;
      sdir='srf';
      noiselev=5E-3;
      indbad=[162:172];
   case '3_nozf'
      module='3';
      nchan=192;
      sdir='srf_nozf';
      noiselev=2E-3;
      indbad=[162:172];
   case '4b'
      module='4b';
      nchan=106;
      sdir='srf';
      noiselev=5E-3;
      indbad=[27:35  49];
   case '4b_nozf'
      module='4b';
      nchan=106;
      sdir='srf_nozf';
      noiselev=3E-3;
      indbad=[27:35  49];
   case '4a'
      module='4a';
      nchan=104;
      sdir='srf';
      noiselev=3E-3;
      indbad=[78];
   case '4a_nozf'
      module='4a';
      nchan=104;
      sdir='srf_nozf';
      noiselev=3E-3;
      indbad=[78];
   case '2b'
      module='2b';
      nchan=150;
      sdir='srf_nozf';
      noiselev=2E-3;
      indbad=[1:12  116  117];
   case '1b'
      module='1b';
      nchan=130;
      sdir='srf_nozf';
      noiselev=1E-3;
      indbad=[120:130];
   case '2a'
      module='2a';
      nchan=116;
      sdir='srf_nozf';
      noiselev=1E-3;
      indbad=[1:7  18  23  107  110];
   case '1a'
      module='1a';
      nchan=118;
      sdir='srf_nozf';
      noiselev=1E-3;
      indbad=[20  107:118];
   otherwise
      error('invalid module name');
end
%%%%%%%%%%%%%
stdlev=13;

ntests=length(testnums);

%%% specify complete dir for data
sdir=['/home/hannon/SRF_results_UMBC/' sdir];
%%%

% Make up arrays for output & initialize with bad data flags
fmidall=-1*ones(nchan,ntests);
noise1all=-1*ones(nchan,ntests);
noise2all=-1*ones(nchan,ntests);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over the test numbers and read in data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:ntests
   itest=testnums(ii);
   disp(['loading data for test ' int2str(itest)])
   [srf, freq, noisearray, indok, id, s]=loadsrf(sdir, itest, module);
   %
   noise1all(indok,ii)=noisearray(1,indok)';
   noise2all(indok,ii)=noisearray(2,indok)';
   %
   % Pull out the not bad channels
   [indgood]=get_good(srf, noisearray, indok, 2*noiselev, stdlev);
   %
   % Measure the centroids of the good channels
   [fmid, fpeak, width]=measure_width(freq, srf, indgood, 0.5);
   %
   fmidall(indgood,ii)=fmid(indgood)';
   %
end

inda=find(gain == 'A'); na=length(inda);
indb=find(gain == 'B'); nb=length(indb);
indo=find(gain == 'O'); no=length(indo);

% Remove temporary variables
clear s srf freq noisearray ii itest fmid fpeak width indgood indok sdir stdlev

id=id';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calc some summary statistics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
amean=-1*ones(nchan,1);
bmean=-1*ones(nchan,1);
omean=-1*ones(nchan,1);
astd=-1*ones(nchan,1);
bstd=-1*ones(nchan,1);
ostd=-1*ones(nchan,1);
alr=zeros(nchan,1);
blr=zeros(nchan,1);
olr=zeros(nchan,1);
an=zeros(nchan,1);
bn=zeros(nchan,1);
on=zeros(nchan,1);
apop=zeros(nchan,1);
bpop=zeros(nchan,1);
opop=zeros(nchan,1);
% Loop over the channels
chanind=(1:nchan)';
for ichan=1:nchan
   %
   % Do A gain
   % Pull out data for this channel
   junk1=noise1all(ichan,inda)/noiselev;
   junk2=noise2all(ichan,inda)/noiselev;
   % Look for non-measurements
   ibad=find(junk1 < 0);
   nbad=length(ibad);
   igood=1:na;
   if (nbad > 0)
      igood=setdiff(igood,ibad);
   end
   n=length(igood);
   junk1=junk1(igood);
   junk2=junk2(igood);
   junk12=0.5*(junk1 + junk2);
   % Look for atypically large noise "pop"
   if (n > 2)
      ipop=find(junk12 > 3*min(junk12));
      apop(ichan)=length(ipop);
      if (apop(ichan) > 0)
         iok=setdiff(1:n, ipop);
      else
         iok=1:n;
      end
   else
      iok=1:n;
   end
   an(ichan)=length(iok);
   if (an(ichan) > 0)
      amean(ichan)=mean( junk12(iok) );
      alr(ichan)=mean(junk1(iok))/mean(junk2(iok));
      if (an(ichan) > 2)
         astd(ichan)=std( junk12(iok) );
      end
   end

   % Do B gain
   % Pull out data for this channel
   junk1=noise1all(ichan,indb)/noiselev;
   junk2=noise2all(ichan,indb)/noiselev;
   % Look for non-measurements
   ibad=find(junk1 < 0);
   nbad=length(ibad);
   igood=1:nb;
   if (nbad > 0)
      igood=setdiff(igood,ibad);
   end
   n=length(igood);
   junk1=junk1(igood);
   junk2=junk2(igood);
   junk12=0.5*(junk1 + junk2);
   % Look for atypically large noise "pop"
   if (n > 2)
      ipop=find(junk12 > 3*min(junk12));
      bpop(ichan)=length(ipop);
      if (bpop(ichan) > 0)
         iok=setdiff(1:n, ipop);
      else
         iok=1:n;
      end
   else
      iok=1:n;
   end
   bn(ichan)=length(iok);
   if (bn(ichan) > 0)
      bmean(ichan)=mean( junk12(iok) );
      blr(ichan)=mean(junk1(iok))/mean(junk2(iok));
      if (bn(ichan) > 2)
         bstd(ichan)=std( junk12(iok) );
      end
   end

   % Do Opt gain
   % Pull out data for this channel
   junk1=noise1all(ichan,indo)/noiselev;
   junk2=noise2all(ichan,indo)/noiselev;
   % Look for non-measurements
   ibad=find(junk1 < 0);
   nbad=length(ibad);
   igood=1:no;
   if (nbad > 0)
      igood=setdiff(igood,ibad);
   end
   n=length(igood);
   junk1=junk1(igood);
   junk2=junk2(igood);
   junk12=0.5*(junk1 + junk2);
   % Look for atypically large noise "pop"
   if (n > 2)
      ipop=find(junk12 > 3*min(junk12));
      opop(ichan)=length(ipop);
      if (opop(ichan) > 0)
         iok=setdiff(1:n, ipop);
      else
         iok=1:n;
      end
   else
      iok=1:n;
   end
   on(ichan)=length(iok);
   if (on(ichan) > 0)
      omean(ichan)=mean( junk12(iok) );
      olr(ichan)=mean(junk1(iok))/mean(junk2(iok));
      if (on(ichan) > 2)
         ostd(ichan)=std( junk12(iok) );
      end
   end

end

fmid=max(fmidall')';

clear ii iok n junk1 junk2 junk12 nbad ibad igood

%%% convert relative noise to absolute noise
amean=amean*noiselev;
bmean=bmean*noiselev;
omean=omean*noiselev;
astd=astd*noiselev;
bstd=bstd*noiselev;
ostd=ostd*noiselev;
%%%

clf
subplot(211)
plot(fmid,amean,'b.',fmid,bmean,'g.',fmid,omean,'r.'),grid
title(['M' modulex])
subplot(212)
plot(fmid,alr,'b.',fmid,blr,'g.',fmid,olr,'r.'),grid
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Write out summary report
%%%%%%%%%%%%%%%%%%%%%%%%%%
report=[chanind fmid amean astd an apop ...
                     bmean bstd bn bpop ...
                     omean ostd on opop];

fname=input('Enter name of noise report file to create : ','s');

eval(['save ' fname ' report -ascii']);

%%% end of function %%%
