% INPUT parameters
% ----------------------------------------------------------------------
pname1 = '/home/preston/matlabs3/bin/workspaces/L2_processed/MCP_72_88_scan/';
%pname1 = '/home/preston/matlabs3/bin/workspaces/L2_processed/5-27_MCPscriptbasedMeasurementsy72-78/';

fname1 = '200_1500_2800_2mvsc1p32mvtr/data.mat'; 

pname2 = '/home/preston/matlabs3/bin/workspaces/L2_processed/MCP_72_88_scan/';
%pname2 = '/home/preston/matlabs3/bin/workspaces/L2_processed/5-27_MCPscriptbasedMeasurementsy72-78/';

fname2 = '200_1500_2800_7mvsc4p76mvtr/data.mat';

cut1 = 1e6;
cut2 = 2e6;
% ---------------------------------------------------------------------
disp([pname1 fname1])
disp([pname2 fname2])

TTn= [pname1 fname1];
load(TTn);

ampsLow = amplifications(:,5);
VmaxsLow = Vmaxs;
wavesLow = newWaves;

% Need to remove some low amplification events before we do rt fwhm and arrival time
% data with low amplification events removed becomes timing*
ampMin = 1.2e6;
[timingWavesLow,timingAmplificationsLow,~] = amplificationControl(newWaves,time,ampMin,2e7,rawevstatsL2);
[timingVmaxsLow,~] = getMaxs(timingWavesLow); 
TTn= [pname2 fname2];
load(TTn);

ampsHigh = amplifications(:,5);
VmaxsHigh = Vmaxs;
wavesHigh = newWaves;

% Need to remove some high amplification events before we do rt fwhm and arrival time

[timingWavesHigh,timingAmplificationsHigh,~] = amplificationControl(newWaves,time,ampMin,2e7,rawevstatsL2);
[timingVmaxsHigh,~] = getMaxs(timingWavesHigh); 

clear TTn pname1 pname2 fname1 fname2
clear DBLbadwaves DHwaves RTFWHMbadwaves rawevstatsL2 oldWaves
clear flattops scopescale triglevel newWaves Vmaxs imaxs

% At this point, we have 2 files opened and extracted the amplifications
% and voltage maxs for the high and low scale.  Cleared everything else.

ww = ampsLow < cut2;
xx = ampsLow > cut1;
yy = ampsHigh > cut1;
zz = ampsHigh < cut2;

clear cut1 cut2

n1 = sum(ww.*xx);
n2 = sum(yy.*zz);

if ( n1 ~= 0 )
    normfac = n2/n1;
end

clear ww xx yy zz

% rise time and fwhm info gathering


[~,gg,~] = size(timingWavesLow);
[~,hh,~] = size(timingWavesHigh);
tinc = time(5)-time(4);

risetimeLow = zeros(1,gg);
risetimeHigh = zeros(1,hh);
fwhmsLow = zeros(1,gg);
fwhmsHigh = zeros(1,hh);
for j = 1:gg
    [~,~,twidthhm,~]=findThresholds(timingWavesLow(:,j,2),.5*timingVmaxsLow(j,2),time);        

    fwhmsLow(j)=twidthhm*tinc; 

    [~,~,rt]=findCFThreshold(timingWavesLow(:,j,2),time,0.2);
    risetimeLow(j)=rt;
end

for j = 1:hh
    [~,~,twidthhm,~]=findThresholds(timingWavesHigh(:,j,2),.5*timingVmaxsHigh(j,2),time);        

    fwhmsHigh(j)=twidthhm*tinc; 

    [~,~,rt]=findCFThreshold(timingWavesHigh(:,j,2),time,0.2);
    risetimeHigh(j)=rt;
end

% Don't want rise times less than .2 ns or fwhms less than 1 ns, they are noisy events contaminating my plots
xx = risetimeLow < .2e-9;
ww = fwhmsLow < .5e-9;
xx = xx+ww;
yy = risetimeHigh < .2e-9;
ww = fwhmsHigh < .5e-9;
yy = yy+ww;


risetimeLow = risetimeLow(~xx);
risetimeHigh = risetimeHigh(~yy);
fwhmsLow = fwhmsLow(~xx);
fwhmsHigh = fwhmsHigh(~yy);
timingWavesLow = timingWavesLow(:,~xx,:);
timingWavesHigh = timingWavesHigh(:,~yy,:);
[timingVmaxsLow,~] = getMaxs(timingWavesLow);
[timingVmaxsHigh,~] = getMaxs(timingWavesHigh);


clear xx yy ww
evstats(1) = ( sum(risetimeLow) + sum(risetimeHigh) )/( gg + hh );
evstats(2) = ( sum(fwhmsLow) + sum(fwhmsHigh) )/( gg + hh );

% Plotting
performL3plots

clear twidthhm rt j hh gg

% error on a bin is sqrt of units within bin






