function [theWaves,time,Vmaxs,imaxs,rawevstats,afterpulsestats,tg,baseMean,baseRMS,badWaves] = processWaves2(whichset)
    
    afterpulsestats = 1;
    %Get the various parameters from whatParameters
    [pname,fname,nfiles,sfile,lowCut,Tfromedge,THcountcut,TWcenter,TWsize,maxTfromWcenter,dhumpTWidth,Tcutonafterpulse,reducerange,zeroORmax,thresholdPercent] = whatParameters(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'];
    [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 to determine size*********************
 
    preallocator = 0;
    
    for nf=sfile:(sfile+nfiles-1)
       
        nft = int2str(nf);
        
        %Open files for all 4 scope channels
        Tfname1= [pname fname nft '.wfmCH1.wfm'];
        Tfname2= [pname fname nft '.wfmCH2.wfm'];
        Tfname3= [pname fname nft '.wfmCH3.wfm'];
        Tfname4= [pname fname nft '.wfmCH4.wfm'];
        
        [y1,~] = ReadTekWFMFile(Tfname1);
        [y2,~] = ReadTekWFMFile(Tfname2); 
        [y3,~] = ReadTekWFMFile(Tfname3);       
        [y4,~] = ReadTekWFMFile(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));
 
        %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);
        passesL2CutB = (maxB>lowR).*(maxB<hiR);
  
        
        preallocator = preallocator + sum(logical(passesL2CutB)); 
        
        
        
    end
    disp(preallocator)
    %clear some variables we no longer need
    clear yTL2a
    clear y1L2a
    clear y2L2a
    clear y2L2a
    clear y3L2a
    clear y4L2a
    
    clear y1L2b
    clear y2L2b
    clear y3L2b
    clear y4L2b
    
    clear nft
    clear Tfname1
    clear Tfname2
    clear Tfname3
    clear Tfname4
    
    yT1 = zeros(1000,preallocator);
    yT2 = zeros(1000,preallocator);
    yT3 = zeros(1000,preallocator);
    yT4 = zeros(1000,preallocator);
    
    
    
    
    
    
    
    
    %Loop over total number of files for real*****************************
 
    rawevcount=0;
    
    for nf=sfile:(sfile+nfiles-1)
       
        nft = int2str(nf);
        
        %Open files for all 4 scope channels
        Tfname1= [pname fname nft '.wfmCH1.wfm'];
        Tfname2= [pname fname nft '.wfmCH2.wfm'];
        Tfname3= [pname fname nft '.wfmCH3.wfm'];
        Tfname4= [pname fname nft '.wfmCH4.wfm'];
        
        [y1,~] = ReadTekWFMFile(Tfname1);
        [y2,~] = ReadTekWFMFile(Tfname2); 
        [y3,~] = ReadTekWFMFile(Tfname3);       
        [y4,~] = ReadTekWFMFile(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);
        passesL2CutB = (maxB>lowR).*(maxB<hiR);

        y1L2b = y1L2a(:,logical(passesL2CutB));
        y2L2b = y2L2a(:,logical(passesL2CutB));
        y3L2b = y3L2a(:,logical(passesL2CutB));
        y4L2b = y4L2a(:,logical(passesL2CutB));   
        
        %collect up all the waves that pass both cuts
        yT1 = [yT1 y1L2b];
        yT2 = [yT2 y2L2b];
        yT3 = [yT3 y3L2b];
        yT4 = [yT4 y4L2b];  
        
        
        size(yT2)
        
        [~,nnnw]=size(y2) ;    
        
        rawevcount=rawevcount+nnnw;
        
    end
    
    %clear some variables we no longer need
    clear yTL2a
    clear y1L2a
    clear y2L2a
    clear y2L2a
    clear y3L2a
    clear y4L2a
    
    clear y1L2b
    clear y2L2b
    clear y3L2b
    clear y4L2b
    
    clear nft
    clear Tfname1
    clear Tfname2
    clear Tfname3
    clear Tfname4
    
    
    
    
    
    
    
    
    
    
    %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' )
    %Here, we apply some basic quality cuts to count the total number of
    %events with after pulses, and to cut out those events where the
    %afterpulse is too close to the initial pulse
    
    %loop over all waves, count number of multiple pulses (in terms of the
    %number of times the waveform crosses a threshold), the number of
    %double humps, and the time between the main pusle and the first
    %after-pulse
    
    thresholdcount = zeros(1,as2);
    thresholdgap = zeros(1,as2);
    dhump = zeros(1,as2);
    
    for j=1:as2
        [tc,tg,dh]=findThresholds2((yT1(:,j)+yT2(:,j)+yT3(:,j)+yT4(:,j)),THcountcut,TWcenter,tt,dhumpTWidth);
        thresholdcount(j)=tc;
        thresholdgap(j)=tg*tinc;
        dhump(j)=dh;
    end


%    variables to cut on how far signal max occurs from expected pulse time
    [~,mmi]=max(yT1+yT2+yT3+yT4);
    ITWcenter = int32((TWcenter-tt(1))/tinc);
    NmaxTfromWcenter = int32(maxTfromWcenter/tinc);

    disp( 'Achieved code level 2' )

%     requirement that "good" waves DON'T have either 
%         -more than one threshold (ie afterpulse) within some time cut
%         -a double hump, formed by the pulse and after pulse
%         -a signal max  at some time much later than the expected arrival time for the initial signal
 
    passesGapCut=~( (((thresholdcount>1).*(thresholdgap<Tcutonafterpulse)) + (dhump==1) + (mmi>(ITWcenter+NmaxTfromWcenter)) )) ;
    
    finalsize=sum(passesGapCut);
    
    oldWaves = zeros(nps,finalsize,4);
    oldWaves(:,:,1) = yT1(:,logical(passesGapCut));
    oldWaves(:,:,2) = yT2(:,logical(passesGapCut));
    oldWaves(:,:,3) = yT3(:,logical(passesGapCut));
    oldWaves(:,:,4) = yT4(:,logical(passesGapCut));

    disp( 'Achieved code level 3' )
    
    clear yT1;
    clear yT2;
    clear yT3;
    clear yT4;
    
%***********************************************Begin reduce range***
    
    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)...
        [theWaves,time,~,baseMean,baseRMS,badWaves] = ReduceRange2(oldWaves,tt,TWsize,TWcenter,zeroORmax,thresholdPercent);
    else
        %...we don't cut down the length of the waves
        theWaves=oldWaves;
        time=tt;
        baseMean=-555;
        baseRMS=-555;
    
    end
    disp( 'Achieved code level 4' )

    %make arrays of the max signal, and index of the max signal of each of
    %the waves
    
    [vmaxs1,imaxs1] = max(theWaves(:,:,1));
    [vmaxs2,imaxs2] = max(theWaves(:,:,2));
    [vmaxs3,imaxs3] = max(theWaves(:,:,3));
    [vmaxs4,imaxs4] = max(theWaves(:,:,4));
    
    [nmaxs,~]=size(vmaxs1');
    
    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;
    
    %************** ADDED BY PRESTON: The following plots all four channels into fig. 1
    [~,MaxWaves,~] = size(theWaves);
    figure(1)
    
    subplot(2,2,1)
    plot(theWaves(:,:,1))
    title('Channel 1')
    ylabel('Voltage (V)')
    xlabel('time (ps)')
    xlabel('Time (seconds)')
    subplot(2,2,2)
    plot(theWaves(:,:,2))
    title('Channel 2')
    ylabel('Voltage (V)')
    xlabel('Time (seconds)')
    subplot(2,2,3)
    plot(theWaves(:,:,3))
    title('Channel 3')
    ylabel('Voltage (V)')
    xlabel('Time (seconds)')
    subplot(2,2,4)
    plot(theWaves(:,:,4))
    title('Channel 4')
    ylabel('Voltage (V)')
    xlabel('Time (seconds)')
    
    
    % : Which channel saw Vmax most frequently?

    figure(2)
    ChannelMaxOccurences = zeros(4);

    MaxOnOneCriteria = ( Vmaxs(:,1) > Vmaxs(:,2) ) .* ( Vmaxs(:,1) > Vmaxs(:,3) ) .* ( Vmaxs(:,1) > Vmaxs(:,4) ); 
    ChannelMaxOccurences(1) = sum( MaxOnOneCriteria );
    MaxOnTwoCriteria = ( Vmaxs(:,2) > Vmaxs(:,1) ) .* ( Vmaxs(:,2) > Vmaxs(:,3) ) .* ( Vmaxs(:,2) > Vmaxs(:,4) ); 
    ChannelMaxOccurences(2) = sum( MaxOnTwoCriteria );
    MaxOnThreeCriteria = ( Vmaxs(:,3) > Vmaxs(:,1) ) .* ( Vmaxs(:,3) > Vmaxs(:,2) ) .* ( Vmaxs(:,3) > Vmaxs(:,4) ); 
    ChannelMaxOccurences(3) = sum( MaxOnThreeCriteria );
    MaxOnFourCriteria = ( Vmaxs(:,4) >= Vmaxs(:,1) ) .* ( Vmaxs(:,4) >= Vmaxs(:,2) ) .* ( Vmaxs(:,4) >= Vmaxs(:,3) ); 
    ChannelMaxOccurences(4) = sum( MaxOnFourCriteria );
    

   
    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('Channel Number')
    
    % At what time did Max tend to occur?
   
    
    ChannelOneEventCounts = zeros(80);
    ChannelTwoEventCounts = zeros(80);
    ChannelThreeEventCounts = zeros(80);

    
    for i = 1 : MaxWaves
        if ( Vmaxs(i,1) > Vmaxs(i,2) && Vmaxs(i,1) > Vmaxs(i,3) && Vmaxs(i,1) > Vmaxs(i,4) )
            ChannelOneEventCounts( imaxs(i,1) ) = ChannelOneEventCounts( imaxs(i,1) ) + 1;
        elseif ( Vmaxs(i,2) > Vmaxs(i,1) && Vmaxs(i,2) > Vmaxs(i,3) && Vmaxs(i,2) > Vmaxs(i,4) )
            ChannelTwoEventCounts( imaxs(i,2) ) = ChannelTwoEventCounts( imaxs(i,2) ) + 1;
        elseif ( Vmaxs(i,3) > Vmaxs(i,1) && Vmaxs(i,3) > Vmaxs(i,2) && Vmaxs(i,3) > Vmaxs(i,4) )
            ChannelThreeEventCounts( imaxs(i,3) ) = ChannelThreeEventCounts( imaxs(i,3) ) + 1;
        end           
    end
%     figure(3)
%     bar(time,ChannelOneEventCounts);
%     title('Time Corresponding to Point of Max Signal for Channel One')
%     ylabel('Number of Signals that Max for a Given Time')
%     xlabel('Time (seconds)')
%     
%     figure(4)
%     bar(time,ChannelTwoEventCounts);
%     title('Time Corresponding to Point of Max Signal for Channel Two')
%     ylabel('Number of Signals that Max for a Given Time')
%     xlabel('Time (seconds)')
%     
%     figure(5)
%     bar(time,ChannelThreeEventCounts)
%     title('Time Corresponding to Point of Max Signal for Channel Three')
%     ylabel('Number of Signals that Max for a Given Time')
%     xlabel('Time (seconds)')
    

    % *****************END OF CONTRIBUTION

end

    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    