function [newWaves,time,Vmaxs,imaxs,rawevstats,afterpulsestats,tg,baseMean,baseRMS,badWaves,RTFWHMbadwaves,DBLbadwaves,notMaxCh2badwaves,amplifications] = L2process(whichset)

    [pname,fname,~,~,~,~,~,~,~,~,~,~,~,~,thresholdPercent,~,~,~] = L2Params(whichset);

    TTn= [pname fname];
    load(TTn);
    
%     numsets=1;
%     
%     if(numsets==1)
% %      theWaves = theWaves1;
% %      badWaves = badWaves1;   
% %      baseMean=baseMean1;
% %      baseRMS=baseRMS1;
% %      Vmaxs=Vmaxs1;
% %      imaxs=imaxs1;
%     elseif(numsets==2)
%      theWaves=[theWaves1 theWaves2];
%      badWaves=[badWaves1 badWaves2];
%      baseMean=[baseMean1 baseMean2];
%      baseRMS=[baseRMS1 baseRMS2];
%      Vmaxs=[Vmaxs1; Vmaxs2];
%      imaxs=[imaxs1; imaxs2];
%     elseif(numsets==3)
%      theWaves=[theWaves1 theWaves2 theWaves3];
%      badWaves=[badWaves1 badWaves2 badWaves3];
%      baseMean=[baseMean1 baseMean2 baseMean3];
%      baseRMS=[baseRMS1 baseRMS2 baseRMS3];
%      Vmaxs=[Vmaxs1; Vmaxs2; Vmaxs3];
%      imaxs=[imaxs1; imaxs2; imaxs3];
%     elseif(numsets==4)
%      theWaves=[theWaves1 theWaves2 theWaves3 theWaves4];
%      badWaves=[badWaves1 badWaves2 badWaves3 badWaves4];
%      baseMean=[baseMean1 baseMean2 baseMean3 baseMean4];
%      baseRMS=[baseRMS1 baseRMS2 baseRMS3 baseRMS4];
%      Vmaxs=[Vmaxs1; Vmaxs2; Vmaxs3; Vmaxs4];
%      imaxs=[imaxs1; imaxs2; imaxs3; imaxs4];
%     end
    
    %time=time1;
    time=theTimescale;
    tinc=time(5)-time(4);
        
    %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
%**************************************************************************
%**************************************************************************
%********************
% Cut for max on ch2
%     isMaxCh2=(Vmaxs(:,2)>Vmaxs(:,1)).*(Vmaxs(:,2)>Vmaxs(:,3));
%     
%     theWaves1=theWaves(:,logical(isMaxCh2),1);
%     theWaves2=theWaves(:,logical(isMaxCh2),2);
%     theWaves3=theWaves(:,logical(isMaxCh2),3);
%     theWaves4=theWaves(:,logical(isMaxCh2),4);
%     
%     % waves that have double pulses go in DBLbadwaves
%     badwaves1=theWaves(:,~logical(isMaxCh2),1);
%     badwaves2=theWaves(:,~logical(isMaxCh2),2);
%     badwaves3=theWaves(:,~logical(isMaxCh2),3);
%     badwaves4=theWaves(:,~logical(isMaxCh2),4);    
%     
%     clear theWaves;
%     
%     theWaves(:,:,1) = theWaves1;
%     notMaxCh2badwaves(:,:,1) = badwaves1;
%     clear theWaves1;
%     clear badwaves1;
%     theWaves(:,:,2) = theWaves2;
%     notMaxCh2badwaves(:,:,2) = badwaves2;
%     clear theWaves2;
%     clear badwaves2;
%     theWaves(:,:,3) = theWaves3;
%     notMaxCh2badwaves(:,:,3) = badwaves3;
%     clear theWaves3;
%     clear badwaves3;
%     theWaves(:,:,4) = theWaves4;
%     notMaxCh2badwaves(:,:,4) = badwaves4;
%     clear theWaves4;
%     clear badwaves4;
    
% Cuts on double pulses
%     [~, as2, ~] = size(theWaves);
%     doubleWaves=zeros(as2,1);
%     for j=1:as2
%         [thcount,~,~,~] = findThresholds(theWaves(:,j,2),thresholdPercent*max(theWaves(:,j,2)),time);
%         if ( thcount == 1 )
%             doubleWaves(j) = 1;
%         else
%             doubleWaves(j) = 0;
%         end
%     end
%     
% 
%     
%     theWaves1=theWaves(:,logical(doubleWaves),1);
%     theWaves2=theWaves(:,logical(doubleWaves),2);
%     theWaves3=theWaves(:,logical(doubleWaves),3);
%     theWaves4=theWaves(:,logical(doubleWaves),4);
%     
%     % waves that have double pulses go in DBLbadwaves
%     badwaves1=theWaves(:,~logical(doubleWaves),1);
%     badwaves2=theWaves(:,~logical(doubleWaves),2);
%     badwaves3=theWaves(:,~logical(doubleWaves),3);
%     badwaves4=theWaves(:,~logical(doubleWaves),4);    
%     
%     clear theWaves;
%     
%     theWaves(:,:,1) = theWaves1;
%     DBLbadwaves(:,:,1) = badwaves1;
%     clear theWaves1;
%     clear badwaves1;
%     theWaves(:,:,2) = theWaves2;
%     DBLbadwaves(:,:,2) = badwaves2;
%     clear theWaves2;
%     clear badwaves2;
%     theWaves(:,:,3) = theWaves3;
%     DBLbadwaves(:,:,3) = badwaves3;
%     clear theWaves3;
%     clear badwaves3;
%     theWaves(:,:,4) = theWaves4;
%     DBLbadwaves(:,:,4) = badwaves4;
%     clear theWaves4;
%     clear badwaves4;

    % Cuts on rise time & FWHM ********************************************
%     [~, as2, ~] = size(theWaves);
%     riseTimes=zeros(as2,1);
%     FWHMs=zeros(as2,1);
%     passes=zeros(as2,1);
%     for j=1:as2
%         riseTimes(j) = getRiseTime(theWaves(:,j,2),thresholdPercent*max(theWaves(:,j,2)),tinc);
%         [~,~,twidth,~] = findThresholds(theWaves(:,j,2),.5*max(theWaves(:,j,2)),time);
%         FWHMs(j) = twidth*tinc;
%         
%         if ( (riseTimes(j) < riseTimeMax) && (FWHMs(j) > fwhmMin) )
%             passes(j) = 1;
%         end
%     end
%     
%     [vmaxs1,~] = max(theWaves(:,:,1));
%     [vmaxs2,~] = max(theWaves(:,:,2));
%     [vmaxs3,~] = max(theWaves(:,:,3));
%     [vmaxs4,~] = max(theWaves(:,:,4));
%     
%     
%     newVmaxs(:,1)=vmaxs1;
%     newVmaxs(:,2)=vmaxs2;
%     newVmaxs(:,3)=vmaxs3;    
%     newVmaxs(:,4)=vmaxs4;
%     
%     [centroids] = getCentroids( newVmaxs );
%     goodCentroids = (centroids<2.4).*(centroids>1.6);
%     
%     passes = passes.*goodCentroids;
%     
%     
%     theWaves1=theWaves(:,logical(passes),1);
%     theWaves2=theWaves(:,logical(passes),2);
%     theWaves3=theWaves(:,logical(passes),3);
%     theWaves4=theWaves(:,logical(passes),4);
%     
%     % waves that have too long a rise time or FWHM go into RTFWHMbadwaves
%     badwaves1=theWaves(:,~logical(passes),1);
%     badwaves2=theWaves(:,~logical(passes),2);
%     badwaves3=theWaves(:,~logical(passes),3);
%     badwaves4=theWaves(:,~logical(passes),4);    
%     
%     clear theWaves;
%     
%     theWaves(:,:,1) = theWaves1;
%     RTFWHMbadwaves(:,:,1) = badwaves1;
%     clear theWaves1;
%     clear badwaves1;
%     theWaves(:,:,2) = theWaves2;
%     RTFWHMbadwaves(:,:,2) = badwaves2;
%     clear theWaves2;
%     clear badwaves2;
%     theWaves(:,:,3) = theWaves3;
%     RTFWHMbadwaves(:,:,3) = badwaves3;
%     clear theWaves3;
%     clear badwaves3;
%     theWaves(:,:,4) = theWaves4;
%     RTFWHMbadwaves(:,:,4) = badwaves4;
%     clear theWaves4;
%     clear badwaves4;
%     
%     
    
    
    % *********************************************************************
    
    [~, as2, ~] = size(theWaves);
    thresholdcount = zeros(1,as2);
    thresholdgap = zeros(1,as2);
    dhump = zeros(1,as2);
    
    for j=1:as2
        [tc,tg,dh]=findThresholds((theWaves(:,j,1)+theWaves(:,j,2)+theWaves(:,j,2)+theWaves(:,j,2)),.2,theTimescale);
        thresholdcount(j)=tc;
        thresholdgap(j)=tg*tinc;
        dhump(j)=dh;
    end




%     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)) )) ;
%     passesGapCut=~( (((thresholdcount>1).*(thresholdgap<Tcutonafterpulse)) )) ;
% %    passesGapCut=~( (((thresholdcount>1).*(thresholdgap<Tcutonafterpulse)) + (dhump==1) )) ;
% 
%     finalsize=sum(passesGapCut);
%     
%     
%    
%     newWaves = zeros(nps,finalsize,4);
%     newWaves(:,:,1) = theWaves(:,logical(passesGapCut),1);
%     clear yT1
%     newWaves(:,:,2) = theWaves(:,logical(passesGapCut),2);
%     clear yT2
%     newWaves(:,:,3) = theWaves(:,logical(passesGapCut),3);
%     clear yT3
%     newWaves(:,:,4) = theWaves(:,logical(passesGapCut),4);
%     clear yT4

    
    newWaves=theWaves;
    disp( 'Achieved code level 2' )
    [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;
    
% ****************** Amplification work ****************************

[asize bsize]=size(theWaves(:,:,2));

AmpFact1=zeros(bsize,1);
AmpFact2=zeros(bsize,1);
AmpFact3=zeros(bsize,1);
AmpFact4=zeros(bsize,1);
AmpFactTot=zeros(bsize,1);


for i=1:bsize
    totWaves = theWaves(:,i,1) + theWaves(:,i,2) + theWaves(:,i,3) + theWaves(:,i,4);
    amp1=IntegrateCharge(theWaves(:,i,1),time,1,asize,12);
    amp2=IntegrateCharge(theWaves(:,i,2),time,1,asize,12);
    amp3=IntegrateCharge(theWaves(:,i,3),time,1,asize,12);
    amp4=IntegrateCharge(theWaves(:,i,4),time,1,asize,12);
    ampTot=IntegrateCharge(totWaves,time,1,asize,12);
    AmpFact1(i,1)=(amp1)/(1.60217646e-19);
    AmpFact2(i,1)=(amp2)/(1.60217646e-19);
    AmpFact3(i,1)=(amp3)/(1.60217646e-19);
    AmpFact4(i,1)=(amp4)/(1.60217646e-19);
    AmpFactTot(i,1)=(ampTot)/(1.60217646e-19);

end

amplifications(:,1) = AmpFact1;
clear AmpFact1;
amplifications(:,2) = AmpFact2;
clear AmpFact2;
amplifications(:,3) = AmpFact3;
clear AmpFact3;
amplifications(:,4) = AmpFact4;
clear AmpFact4;
amplifications(:,5) = AmpFactTot;
clear AmpFactTot;

figure(1);
scatter(Vmaxs(:,1),amplifications(:,1),1)
title('Amplification vs Signal Max Voltage on Channel 1');
ylabel('Amplification');
xlabel('Max Signal Voltage');


figure(2);
scatter(Vmaxs(:,2),amplifications(:,2),1)
title('Amplification vs Signal Max Voltage on Channel 2');
ylabel('Amplification');
xlabel('Max Signal Voltage');


figure(3);
scatter(Vmaxs(:,3),amplifications(:,3),1)
title('Amplification vs Signal Max Voltage on Channel 3');
ylabel('Amplification');
xlabel('Max Signal Voltage');


figure(4);
scatter(Vmaxs(:,4),amplifications(:,4),1)
title('Amplification vs Signal Max Voltage on Channel 4');
ylabel('Amplification');
xlabel('Max Signal Voltage');


figure(5);
scatter(Vmaxs(:,1)+Vmaxs(:,2)+Vmaxs(:,3)+Vmaxs(:,4),amplifications(:,5),1)
title('Amplification vs Signal Max Voltage on the Sum of All Channels');
ylabel('Amplification');
xlabel('Max Signal Voltage');


%**************************************************************************
%**************************************************************************
%******************double hump counting end

rawevstats=1;
afterpulsestats=1;

[~,as2,~] = size(theWaves);
priseTimes=zeros(as2,1);
pFWHMs=zeros(as2,1);

    for j=1:as2
        priseTimes(j) = getRiseTime(theWaves(:,j,2),thresholdPercent*max(theWaves(:,j,2)),tinc);
        [~,~,twidth,~] = findThresholds(theWaves(:,j,2),.5*max(theWaves(:,j,2)),time);
        pFWHMs(j) = twidth*tinc;
        

    end
    figure(6)
    scatter(amplifications(:,2),priseTimes);
    title('risetime')
    
    figure(7)
    rtbins=0:4e-10:70e-9;
    hist(priseTimes,rtbins)
    title('risetime')

    figure(8)
    scatter(amplifications(:,2),pFWHMs);
    title('fwhm')
    
    figure(9)
    rttbins=-2e-9:2e-10:8e-9;
    hist(pFWHMs,rttbins)
    title('fwhm')
end