function [theWaves,theTimescale,Vmaxs,imaxs,rawevstats,baseMean,baseRMS,badWaves,lightsource,scopescale,triglevel] = processWaves(whichset)
    

    %Get the various parameters from procWavesParams
    [pname,fname,nfiles,sfile,lowCut,Tfromedge,THcountcut,TWcenter,TWsize,~,~,~,reducerange,zeroORmax,thresholdPercent,reducedRangeWidth,lightsource,scopescale,triglevel] = procWavesParams(whichset);

    %Open one test file, to determine the number of size of the arrays
    %storing the pulse information are stored
    %TTn= [pname fname '1.wfmCH1.wfm'];
%    TTn= [pname fname '1_CH1.wfm'];
    TTn= [pname fname '1.wfmCH1.wfm'];
    [yt,tt] = ReadTekWFMFile(TTn);   
    [nps,~] = size(yt);

    %Initiate the variables to store the pulse data
    yT1 = zeros(nps,1);
    yT2 = zeros(nps,1);
    yT3 = zeros(nps,1);
    yT4 = zeros(nps,1);
    
    %the sampling increment
    tinc = tt(5,1)-tt(4,1);
    
    %When we used pulsed LED data, the peak times vary over the whole
    %time interval being recorded by the scope. To make sure that we 
    %have the entire pulse contained within a given time window, we 
    %require the signal max to be a certain amount of time away from 
    %the edges of the scopre's time range...This is done in terms of 
    %the indices of the pulse array
    NstepsFromEdge=int32(Tfromedge/tinc);
    lowR=NstepsFromEdge;
    hiR=nps-NstepsFromEdge;
    
    %Loop over total number of files**************************************
 
    rawevcount=0;
    
    for nf=sfile:(sfile+nfiles-1)
       
        nft = int2str(nf);
        
        %Open files for all 4 scope channels
%         Tfname1= [pname fname nft '_CH1.wfm'];
%         Tfname2= [pname fname nft '_CH2.wfm'];
%         Tfname3= [pname fname nft '_CH3.wfm'];
%         Tfname4= [pname fname nft '_CH4.wfm'];
        
        Tfname1= [pname fname nft '.wfmCH1.wfm'];
        Tfname2= [pname fname nft '.wfmCH2.wfm'];
        Tfname3= [pname fname nft '.wfmCH3.wfm'];
        Tfname4= [pname fname nft '.wfmCH4.wfm'];
        clear nft
        
        [y1,~] = ReadTekWFMFile(Tfname1);
        clear Tfname1
        [y2,~] = ReadTekWFMFile(Tfname2);
        clear Tfname2
        [y3,~] = ReadTekWFMFile(Tfname3);
        clear Tfname3
        [y4,~] = ReadTekWFMFile(Tfname4);  
        clear Tfname4
        
        %sum the signal on all 4 channels to apply offline quality cut on
        %total signal
        ytotall = y4+y3+y2+y1;
        passesL2Cut = max(ytotall(:,:) > lowCut);
        
        yTL2a = ytotall(:,logical(passesL2Cut));
        y1L2a = y1(:,logical(passesL2Cut));
        y2L2a = y2(:,logical(passesL2Cut));
        y3L2a = y3(:,logical(passesL2Cut));
        y4L2a = y4(:,logical(passesL2Cut));
 
        %require that the index of the pulse maximum be a certain number of
        %steps away from the time interval measured by the sope
        [~,maxB] = max(yTL2a);
        clear yTL2a
        passesL2CutB = (maxB>lowR).*(maxB<hiR);

        y1L2b = y1L2a(:,logical(passesL2CutB));
        clear y1L2a
        y2L2b = y2L2a(:,logical(passesL2CutB));
        clear y2L2a
        y3L2b = y3L2a(:,logical(passesL2CutB));
        clear y3L2a
        y4L2b = y4L2a(:,logical(passesL2CutB));   
        clear y4L2a
        
        %collect up all the waves that pass both cuts
        yT1 = [yT1 y1L2b];
        clear y1L2b
        yT2 = [yT2 y2L2b];
        clear y2L2b
        yT3 = [yT3 y3L2b];
        clear y3L2b
        yT4 = [yT4 y4L2b];  
        clear y4L2b        
        
        size(yT2)
        
        [~,nnnw]=size(y2) ;    
        
        rawevcount=rawevcount+nnnw;
        
    end
   
    
    %determine the number of waves stored after all the initial quality
    %cuts
    rawevstats(1)=rawevcount;
    [~, as2] = size(yT2);
    
    % First wave is ignored
    yT1=yT1(:,2:as2);
    yT2=yT2(:,2:as2);
    yT3=yT3(:,2:as2);
    yT4=yT4(:,2:as2);
    
    [~, as2] = size(yT2);%necesary repeat
    rawevstats(2) = as2;
    %**********************************************************************
    disp( 'Achieved code level 1' )

    oldWaves = zeros(nps,as2,4);
    oldWaves(:,:,1) = yT1;
    clear yT1
    oldWaves(:,:,2) = yT2;
    clear yT2
    oldWaves(:,:,3) = yT3;
    clear yT3
    oldWaves(:,:,4) = yT4;
    clear yT4
    
    
    thresholdcount = zeros(1,as2);
    thresholdgap = zeros(1,as2);
    fwhms = zeros(1,as2);
    risetime = zeros(1,as2);
    threshwidths = zeros(1,as2);
    maxTs = zeros(1,as2);
    edgeTs = zeros(1,as2);
    maxTdiffs = zeros(1,as2);
    
    
    for j=1:as2

         [tc,tg,twidth,edgeT]=findThresholds((oldWaves(:,j,1)+oldWaves(:,j,2)+oldWaves(:,j,3)+oldWaves(:,j,4)),THcountcut,tt);
         thresholdcount(j)=tc;
         thresholdgap(j)=tg*tinc;
%                  thresholdgap(j)=tg;
         threshwidths(j)=twidth*tinc;
         edgeTs(j)=edgeT;
         
%         disp( 'half max now' )
         
        [mval,im]=max( (oldWaves(:,j,1)+oldWaves(:,j,2)+oldWaves(:,j,3)+oldWaves(:,j,4)) );       
        hmax = 0.5*mval;
               
        [~,~,twidthhm,~]=findThresholds((oldWaves(:,j,1)+oldWaves(:,j,2)+oldWaves(:,j,3)+oldWaves(:,j,4)),hmax,tt);        

        fwhms(j)=twidthhm*tinc; 
        
        [~,~,rt]=findCFThreshold(oldWaves(:,j,2),tt,0.3);
        risetime(j)=rt;
        
        maxTs(j)=tt(im);
         
        maxTdiffs(j)=maxTs(j)-edgeT;

     end
    
    
    disp( 'Achieved code level 2' )

%***********************************************Begin reduce range***
    
 %   dude=size(oldWaves)

    if(reducerange==1)
        % Here we cut the waves down to the range centered around the
        % typical pulse peak or centered around the pulse max (depending on
        % the value of zeroORmax)...
        [theWaveso,time,~,baseMeano,baseRMSo,badWaves,reducedTime] = ReduceRange(oldWaves,tt,TWsize,TWcenter,zeroORmax,thresholdPercent,reducedRangeWidth);
    else
        %...we don't cut down the length of the waves
        theWaveso=oldWaves;
        time=tt;
        [baseMeano,baseRMSo]=baselineAna(oldWaves,100);
        badWaves=-555;    
    end
    disp( 'Achieved code level 3' )

    %make arrays of the max signal, and index of the max signal of each of
    %the waves
    
    [vmaxs1o,imaxs1o] = max(theWaveso(:,:,1));
    [vmaxs2o,imaxs2o] = max(theWaveso(:,:,2));
    [vmaxs3o,imaxs3o] = max(theWaveso(:,:,3));
    [vmaxs4o,imaxs4o] = max(theWaveso(:,:,4));
    
    maxmax(1)=max(vmaxs1o);
    maxmax(2)=max(vmaxs2o);
    maxmax(3)=max(vmaxs3o);
    maxmax(4)=max(vmaxs4o);
    measuredscale = (max(maxmax))/8.;
    
    minmax(1)=min(vmaxs1o);
    minmax(2)=min(vmaxs2o);   
    minmax(3)=min(vmaxs3o);    
    minmax(4)=min(vmaxs4o);
    measuredtrigger = minmax(2);
    
    if(scopescale>0) 
        scalediscrepancy = abs(scopescale-measuredscale)/scopescale;
    else
        disp('Whoa! Scope scale set to a nonsensical number!!!');
    end
    
    if(scalediscrepancy>0.25)
        disp(' ');
        disp('************************************************************************************');
        disp('************************************************************************************');
        disp('ALERT!!!! Just to let you know: the effective scope scale measured in the data of');
        disp(measuredscale)
        disp('does not agree to your hard coded value of');
        disp(scopescale)
        disp('to within 0.25');                
        disp('************************************************************************************');                        
        disp('************************************************************************************');
        disp(' ');        
    end
    
    if(triglevel>0)
        
        disp(triglevel)
        disp(measuredtrigger)
        disp(((triglevel-measuredtrigger)/triglevel))
        
        
    
        if( ((triglevel-measuredtrigger)/triglevel)<-0.1 )
            
                disp(' ');
                disp('************************************************************************************');
                disp('************************************************************************************');
                disp('ALERT!!!! Just to let you know: the effective trigger level measured in the data of');
                disp(measuredtrigger)
                disp('is higher than your hard coded value of');
                disp(triglevel)
                disp('by more than 0.1');
                disp('************************************************************************************');
                disp('************************************************************************************');
                disp(' ');
        end
        
    end
    
    
    [nmaxso,~]=size(vmaxs1o');
    
    rawevstats(3) = nmaxso;
    
    vmaxs1o;
    nmaxso;
    
    Vmaxso=zeros(squeeze(nmaxso),4);
    imaxso=zeros(squeeze(nmaxso),4);
    
    Vmaxso(:,1)=vmaxs1o;
    Vmaxso(:,2)=vmaxs2o;
    Vmaxso(:,3)=vmaxs3o;    
    Vmaxso(:,4)=vmaxs4o;
    
    imaxso(:,1)=imaxs1o;
    imaxso(:,2)=imaxs2o;
    imaxso(:,3)=imaxs3o;    
    imaxso(:,4)=imaxs4o;
        
    ChannelMaxOccurences = zeros(4);

    MaxOnOneCriteria = ( Vmaxso(:,1) > Vmaxso(:,2) ) .* ( Vmaxso(:,1) > Vmaxso(:,3) ) .* ( Vmaxso(:,1) > Vmaxso(:,4) );

    MaxOnTwoCriteria = ( Vmaxso(:,2) > Vmaxso(:,1) ) .* ( Vmaxso(:,2) > Vmaxso(:,3) ) .* ( Vmaxso(:,2) > Vmaxso(:,4) );
%    MaxOnTwoCriteria = ( Vmaxso(:,2) > Vmaxso(:,1) ) .* ( Vmaxso(:,2) >
%    Vmaxso(:,3) ) ;
    MaxOnThreeCriteria = ( Vmaxso(:,3) > Vmaxso(:,1) ) .* ( Vmaxso(:,3) > Vmaxso(:,2) ) .* ( Vmaxso(:,3) > Vmaxso(:,4) );

    MaxOnFourCriteria = ( Vmaxso(:,4) >= Vmaxso(:,1) ) .* ( Vmaxso(:,4) >= Vmaxso(:,2) ) .* ( Vmaxso(:,4) >= Vmaxso(:,3) );

    
    ChannelMaxOccurences(1) = sum( MaxOnOneCriteria );
    ChannelMaxOccurences(2) = sum( MaxOnTwoCriteria ) + sum( MaxOnFourCriteria );
    ChannelMaxOccurences(3) = sum( MaxOnThreeCriteria );
    ChannelMaxOccurences(4) = sum( MaxOnTwoCriteria );
    ChannelMaxOccurences(5) = sum( MaxOnFourCriteria );
        
    MaxOnTwoCriteria = MaxOnTwoCriteria + MaxOnFourCriteria;
    
    theWaves=theWaveso(:,logical(MaxOnTwoCriteria),:);
    
    [vmaxs1,imaxs1] = max(theWaves(:,:,1));
    [vmaxs2,imaxs2] = max(theWaves(:,:,2));
    [vmaxs3,imaxs3] = max(theWaves(:,:,3));
    [vmaxs4,imaxs4] = max(theWaves(:,:,4));
    
    [nmaxs,~]=size(vmaxs1');
    
    rawevstats(4) = nmaxs;
    
    vmaxs1;
    nmaxs;
    
    Vmaxs=zeros(squeeze(nmaxs),4);
    imaxs=zeros(squeeze(nmaxs),4);
    
    Vmaxs(:,1)=vmaxs1;
    Vmaxs(:,2)=vmaxs2;
    Vmaxs(:,3)=vmaxs3;    
    Vmaxs(:,4)=vmaxs4;
    
    imaxs(:,1)=imaxs1;
    imaxs(:,2)=imaxs2;
    imaxs(:,3)=imaxs3;    
    imaxs(:,4)=imaxs4;

    baseMean=baseMeano(:,logical(MaxOnTwoCriteria));
    baseRMS=baseRMSo(:,logical(MaxOnTwoCriteria));  

%    baseMean=baseMeano;
%    baseRMS=baseRMSo;      
   
    size(MaxOnTwoCriteria)
    size(baseMeano)
    
    
    disp( 'Achieved code level 4' )
    

    if (reducerange == 1)
        theTimescale = reducedTime';
    else
        theTimescale = time;
    end

    [timeSize,~] = size(theTimescale);
    
    [~,MaxWaves,dm2] = size(theWaves);

    figure(1)
    subplot(2,2,1)
    plot(theTimescale,theWaves(:,:,1));
    title('Channel 1')

    ylabel('Voltage (V)')
    xlabel('time (seconds)')
    subplot(2,2,2)
    plot(theTimescale,theWaves(:,:,2))
    title('Channel 2')
    ylabel('Voltage (V)')
    xlabel('time (seconds)')
    subplot(2,2,3)
    plot(theTimescale,theWaves(:,:,3))
    title('Channel 3')
    ylabel('Voltage (V)')
    xlabel('time (seconds)')
    subplot(2,2,4)
    plot(theTimescale,theWaves(:,:,4))
    title('Channel 4')
    ylabel('Voltage (V)')
    xlabel('time (seconds)')

    %: Which channel saw Vmax most frequently?

    figure(2)

    subplot(2,2,1)
    bar( ChannelMaxOccurences(1:3))
    title('Precision of Light Being Centered about Channel Two')
    ylabel('Number of Incidents Where Max Signal Occurs at a Given Channel')
    xlabel('Stripline Number')
    
    subplot(2,2,2)
    centroid= (Vmaxs(:,1)+2*Vmaxs(:,2)+3*Vmaxs(:,3))./((Vmaxs(:,1)+Vmaxs(:,2)+Vmaxs(:,3)));
    centbins=1:0.1:3;
    hist(centroid,centbins);
    
    h3 = findobj(gca,'Type','patch');
    set(h3,'FaceColor','r','EdgeColor','w');
    title('Centroids of Pulses');
    xlabel('strip coordinate');
    
    
    % At what time did Max tend to occur?
    % Have not confirmed correctness of time for reduceRange = 1
    % scenario!!!!


    maxonCH1 = ( (Vmaxs(:,1) > Vmaxs(:,2)).*(Vmaxs(:,1) > Vmaxs(:,3)).*(Vmaxs(:,1) > Vmaxs(:,4)) );
    maxonCH2 = ( (Vmaxs(:,2) > Vmaxs(:,1)).*(Vmaxs(:,2) > Vmaxs(:,3)).*(Vmaxs(:,2) > Vmaxs(:,4)) );
    maxonCH3 = ( (Vmaxs(:,3) > Vmaxs(:,1)).*(Vmaxs(:,3) > Vmaxs(:,2)).*(Vmaxs(:,3) > Vmaxs(:,4)) );

    subplot(2,2,3:4)

    CH1imaxs = squeeze(imaxs(logical(maxonCH1),1));
    CH3imaxs = squeeze(imaxs(logical(maxonCH3),1));
    CH2imaxs = squeeze(imaxs(logical(maxonCH2),2));
    CH4imaxs = squeeze(imaxs(logical(maxonCH2),4));
    imaxbins = 0.5:1:(timeSize-0.5);

    ChannelTwoEventCounts = hist(CH2imaxs,imaxbins);
    ChannelOneEventCounts = hist(CH1imaxs,imaxbins);
    ChannelThreeEventCounts = hist(CH3imaxs,imaxbins);
    ChannelFourEventCounts = hist(CH4imaxs,imaxbins);

    taxis=(theTimescale(1:timeSize));
    stairs(taxis,ChannelTwoEventCounts,'k');
    hold on
    stairs(taxis,ChannelOneEventCounts,'r');
    stairs(taxis,ChannelThreeEventCounts,'b');
    stairs(taxis,ChannelFourEventCounts,'m');

    title('Time Corresponding to Point of Max Signal for Each Channel')
    ylabel('Number of Signals that Max for a Given Time')
    xlabel('Time (seconds)')

    h = legend('CH2','CH1','CH3','CH4');
    set(h,'Interpreter','none');

    figure(3)

    baseMbins=-2e-3:5e-6:8e-3;
    subplot(2,1,1)
    hist(baseMean(2,:),baseMbins)
    title('Histogram of Mean Pedistals')
    xlabel('Mean Pedistal (Volts)')
    
    subplot(2,1,2)
    hist(baseRMS(2,:),baseMbins);
    title('Histogram of Pedistal RMSs (noise)')
    xlabel('Pedistal RMS')


    figure(4)
%     scrsz = get(0,'ScreenSize');
%     figure('Position',[1 scrsz(4) scrsz(3) scrsz(4)])
    subplot(2,2,1);
    nedgebins=0:1:10;
    hist(thresholdcount,nedgebins)
    title('Histogram of Threshold Counts')
    xlabel('# of Thresholds')

    subplot(2,2,2);
    isbiggerthanzero=(thresholdgap(:)>0.);
    %      isbiggerthanzero=((thresholdgap(:)>0.).*(thresholdgap(:)<30.));
    newthreshgap=thresholdgap(logical(isbiggerthanzero));

    tbinsyo=-1.e-8:1e-9:9e-8;
    %   tbinsyo=0:1:200;

    hist(newthreshgap,tbinsyo);
    h2 = findobj(gca,'Type','patch');
    set(h2,'FaceColor','g','EdgeColor','b');
    title('Histogram of Time Delay Till Second Pulse')
    xlabel('time delay (seconds)')
    

    subplot(2,2,3);
    tbinsyo2=-2e-9:2e-10:8e-9;
%    hist(fwhms,tbinsyo2,'y');
    hist(fwhms(logical(MaxOnTwoCriteria)),tbinsyo2,'y');
    title('FWHM Distribution')
    xlabel('FWHM (seconds)')
    
    subplot(2,2,4)
    
    tbinsyo2=0:4e-10:70e-9;
%    hist(risetime,tbinsyo2);
    hist(risetime(logical(MaxOnTwoCriteria)),tbinsyo2);
    title('Rise Time Distribution')
    xlabel('rise time (seconds)')
    

    
    % *****************END OF CONTRIBUTION
    disp(fname);
end

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    