%% run this fro final analysis
clearvars -except allSess
close all

%% loading the necessary files and variable  names
targetFile = dir('*.MergePoints.events.mat'); load(targetFile(1).name);  
targetFile = dir('*.hippocampalLayers.channelinfo.mat'); load(targetFile(1).name); 
targetFile = dir('*sessionInfo.mat'); load(targetFile(1).name);
figure_title = sessionInfo.FileName;
ripplech = hippocampalLayers.bestShankLayers.pyramidal;
cutwind = 50; % number of datapoints taken before and after the event

%% creating two sets of detected events with different filters

[Epi2]= findRipples_a(ripplech,'thresholds',[3 10],'passband',[150 300],...
    'EMGThresh',1,'durations',[50 250], 'saveMat',false,'restrict',[MergePoints.timestamps(1,1) MergePoints.timestamps(1,2) ],'frequency',1250);


[Epi1]= findRipples_a(ripplech,'thresholds',[3 10],'passband',[50 80],...
    'EMGThresh',1,'durations',[50 250], 'saveMat',false,'restrict',[MergePoints.timestamps(1,1) MergePoints.timestamps(1,2) ],'frequency',1250);

%% selecting the common ones

distance = abs(Epi1.peaks - Epi2.peaks');
[epi_ind,epi2_ind] =find(distance <= 0.1);


Epi.timestamps = Epi1.timestamps(epi_ind,:);
Epi.peaks = Epi1.peaks(epi_ind);
Epi.peakNormedPower = Epi1.peakNormedPower(epi_ind,:);
Epi.lfp = Epi1.lfp;
Epi.signal =Epi1.signal;
%% adapted from IEDplot

validEvents = find(Epi.peakNormedPower>=5);
Epi_Events.timestamps = Epi.timestamps(validEvents,:);
try 
   Epi_Events.peaks = Epi.peaks(validEvents,:);
   Epi_Events.peakNormedPower = Epi.peakNormedPower(validEvents,:);
end

%%
initial_index = zeros(1,length(Epi_Events.peaks));
single_events = zeros(2*cutwind+1,length(Epi_Events.peaks));
single_lowpass = zeros(2*cutwind+1,length(Epi_Events.peaks));
for yy = 1:length(Epi_Events.peaks)
initial_index(yy)=find(Epi.lfp.timestamps == Epi_Events.peaks(yy));
single_events(:,yy) = Epi.lfp.data((initial_index(yy)-cutwind):(initial_index(yy)+cutwind));
single_lowpass(:,yy) = Epi.signal((initial_index(yy)-cutwind):(initial_index(yy)+cutwind)); 
end

[M,I] = max(single_events); %maximum values and indices of single events
index = (initial_index - cutwind)+I-1; %creates an index of the peaks

%%
% remove big distances between initially registered and real
h = figure;
hist(I,100);
title('Select threshold with the mouse and press right click...');
[cutind1,thr1] = ginput(1);
hold on
line(cutind1*[1 1],ylim,'Color','r');
title('Select threshold with the mouse and press right click...');
[cutind2,thr2] = ginput(1);
hold on
line(cutind2*[1 1],ylim,'Color','c');
pause(1);
close(h);
A = find(I>=cutind1 & I<=cutind2);
index = index(A);
initial_index = initial_index(A);nts.peakNormedPower = Epi.peakNormedPower(A,:);


Peaks_timestamps = Epi.lfp.timestamps(index);
Epi_Events.peaks = Epi_Events.peaks(A);
Epi_Events.timestamps = Epi_Events.timestamps(A,:);
Epi_Events.peakNormedPower = Epi_Events.peakNormedPower(A);
%% plot and save

mkdir('Analysis');

%% figure for testing
for ii = 1:9:9*(fix(size(A,2)/9))
h = figure;
multiplot;
prompt = 'Do you want more? Y/N [Y][N]: ';
str = input(prompt, 's');
if isempty(str)
    str = 'Y';
else str = 'N';
end
    if str == 'Y';
      
        continue
    else
break
    end
end


%%fig
%%
h = figure;
plot (Epi.lfp.timestamps, Epi.lfp.data);
hold on
plot (Epi.lfp.timestamps, Epi.signal);
plot(Peaks_timestamps,Epi.lfp.data(index),'o');
title("session " + figure_title(end-1:end) + "")
% scatter(Epi_Events.timestamps(:,1),Epi_Events.peakNormedPower);
% scatter(Epi_Events.timestamps(:,2),Epi_Events.peakNormedPower);
xlabel ('s');
axis tight
saveas(gcf,'Analysis\IED_sess.fig');  
saveas(gcf,'Analysis\IED_sess.jpeg');