function [dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_eventAll] = getFiringRateEvents(popRate,dFoverF,sleepState,thresholds,time_window)
%windowSize = 0.5; % Size of the sliding window
%threshold = 0; % Frequency threshold
% Calculate event frequencies using sliding window approach
fs_dff = 100;
N_dFF = length(dFoverF);

%%
durationPeakDiff = 10;
[pka, locs_a, ~] = findpeaks(popRate,'minpeakheight',mean(popRate)+thresholds(1)*std(popRate),...
    'WidthReference','halfprom','MinPeakDistance',durationPeakDiff*2);
[pkb, locs_b, ~] = findpeaks(popRate,'minpeakheight',mean(popRate)+thresholds(2)*std(popRate),...
    'WidthReference','halfprom','MinPeakDistance',durationPeakDiff*2);
[locs,pks] = setdiff(locs_a,locs_b);
% [locs12,pks12] = setdiff(locs_a,locs_b);
% figure
% t_plot = 1:length(popRate);
% plot(t_plot,popRate,"Color", [0,0,0.5,0.3]);
% colrs=colorGradient([0.75 0.75 0.75],[0,0,0],4);
% 
% hold on
% plot(locs12,pka(pks12),'kx')
% yline(mean(popRate)+1*std(popRate),"Color",colrs(2,:),"LineWidth",2)
% 
% [pk3, locs_3, ~] = findpeaks(popRate,'minpeakheight',mean(popRate)+3*std(popRate),...
%     'WidthReference','halfprom','MinPeakDistance',durationPeakDiff*2);
% [locs3m,pks3m] = intersect(locs_3,locs_b);
% plot(locs3m,pk3(pks3m),'kx')
% yline(mean(popRate)+3*std(popRate),"Color",colrs(4,:),"LineWidth",2)
% 
% [locs23,pks23] = setdiff(locs_b,locs_3);
% plot(locs23,pkb(pks23),'kx')
% yline(mean(popRate)+2*std(popRate),"Color",colrs(3,:),"LineWidth",2)
% 
% [pk4, locs_4, ~] = findpeaks(popRate,'minpeakheight',mean(popRate)+0.5*std(popRate),...
%     'WidthReference','halfprom','MinPeakDistance',durationPeakDiff*2);
% [locs44,pks44] = setdiff(locs_4,locs_a);
% plot(locs44,pk4(pks44),'kx')
% yline(mean(popRate)+0.5*std(popRate),"Color",colrs(1,:),"LineWidth",2)
% xlim([3e4 4e4])
% box off
%%

eventTimes = locs;

dFF_eventNREM = [];
dFF_eventWAKE = [];
dFF_eventREM = [];
dFF_eventAll = [];
for i = 1:length(eventTimes)
    idx = eventTimes(i);
    if (idx-sum(time_window)*fs_dff >0) && (idx+sum(time_window)*fs_dff<N_dFF)
        isState = sleepState(idx-time_window(1)*fs_dff:idx+time_window(2)*fs_dff);
        a = dFoverF(idx-time_window(1)*fs_dff:idx+time_window(2)*fs_dff);
        if ~any(isState ~=3,"all")
            dFF_eventNREM = [dFF_eventNREM a];
        elseif ~any(isState ~=1,"all")
            dFF_eventWAKE = [dFF_eventWAKE a];
        elseif ~any(isState ~=5,"all")
            dFF_eventREM = [dFF_eventREM a];
        end
        dFF_eventAll = [dFF_eventAll a];

    end
end

miu = mean(dFF_eventAll(1:1000,:));
sigma_ = std(dFF_eventAll(1:1000,:));
dFF_eventAll = (dFF_eventAll - miu)./sigma_;

if ~isempty(dFF_eventWAKE)
    miu = mean(dFF_eventWAKE(1:1000,:));
    sigma_ = std(dFF_eventWAKE(1:1000,:));
    dFF_eventWAKE = (dFF_eventWAKE - miu)./sigma_;
end
if ~isempty(dFF_eventNREM)
    miu = mean(dFF_eventNREM(1:1000,:));
    sigma_ = std(dFF_eventNREM(1:1000,:));
    dFF_eventNREM = (dFF_eventNREM - miu)./sigma_;
end
if ~isempty(dFF_eventREM)
    miu = mean(dFF_eventREM(1:1000,:));
    sigma_ = std(dFF_eventREM(1:1000,:));
    dFF_eventREM = (dFF_eventREM - miu)./sigma_;
end



