%% detect ID Spike
% clear all
clearvars -except allSess
close all

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

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


%figure; subplot(1,2,1); hist(Epi.peakNormedPower,100);
%  h = figure;
%         hist(Epi.peakNormedPower,100);
%         title('Select threshold with the mouse and press right click...');
%         [cutpoint,thr] = ginput(1);
%         hold on
%         line(cutpoint*[1 1],ylim,'Color','r');
%         pause(1);
%         close(h);
% validEvents = find(Epi.peakNormedPower>=cutpoint);
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);
Epi_Events.peaks=[];
Epi_Events.peaks=Peaks_timestamps;

mkdir('Analysis')
%% plot a 3x3 figure with individual events
 
if size(A,2) >= 36 
    z = 3;
else z = fix(size(A,2)/9);
end


for ii = 1:9:9*z
h = figure;
multiplot;
title("session " + figure_title(end-1:end) + "")
%saveas(h,sprintf('Analysis\multiplot%d.png',ii))
saveas(h,'Analysis\multiplot.fig');  
saveas(h,'Analysis\multiplot.jpeg'); 
end


% for ii = 1:9:9*(fix(size(A,2)/9))
% figure;
% gcf;
% 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

%% plotting the whole trace;

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');  

%% wavelet;

lfp = Epi.lfp;

if size(Peaks_timestamps) >= 50;
[wavAvg, lfpAvg ] = bz_eventWavelet (lfp, Peaks_timestamps (1:50), 'twin',[.3 .3], 'frange',[6 90], 'whitening',false);
caxis([1000 30000])
else
[wavAvg, lfpAvg ] = bz_eventWavelet (lfp, Peaks_timestamps, 'twin',[.3 .3], 'frange',[10 300], 'whitening',false);
caxis([1000 20000])
end
title (" session " + figure_title(end-1:end) + " f_m_a_x " +wavAvg.maxFreq+ " pow_m_a_x " +wavAvg.maxPow+ "")
saveas(gcf,'Analysis\IED.fig');  
saveas(gcf,'Analysis\IED.jpeg');  


%% isolate events outside behviour

%ms removal of behavior
a = InIntervals(Epi.peaks,[MergePoints.timestamps(1,1),MergePoints.timestamps(1,end);MergePoints.timestamps(end,1), MergePoints.timestamps(end, end)]);
fieldsR = fields(Epi);
            for j = 1:size(fieldsR,1)
                if ~isstruct(Epi.(fieldsR{j})) && size(Epi.(fieldsR{j}),1) > 3
                    Epi.(fieldsR{j})(find(a==0),:) = [];
                end
            end
% get Epi Stats             
 Epi = computeRippleStats('Ripples',Epi);
 
 
% check for  noise
figure;
    plot(Epi.rippleStats.maps.timestamps, Epi.rippleStats.maps.ripples_raw/1E3,'color', [.1 .1 .1 .1]); hold on;
    plot(Epi.rippleStats.maps.timestamps, median(Epi.rippleStats.maps.ripples_raw)/1E3,'color', [0 0 0],'LineWidth',1.5);
    axis tight
    xlabel('s'); ylabel('Amp');title('Not Clean');
    
%% saving

filename = split(pwd,filesep); filename = filename{end};
save([filename '.IED.mat'],'Epi_Events');



