function [dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_eventAll] = getFiringRateEvents(popRate,dFoverF,acc,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);

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-10*fs_dff:idx+10*fs_dff);
        isAcc = acc(idx-10*fs_dff:idx+10*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