function [dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_eventAll] = getFiringRateEventsAch(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-15*fs_dff:idx+15*fs_dff);
        a = dFoverF(idx-time_window(1)*fs_dff:idx+time_window(2)*fs_dff);
        if ~any(isState ~=3,"all")
            if ~any(isAcc >20)
                dFF_eventNREM = [dFF_eventNREM a];
            end
        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

wbaseline = fs_dff*10;
miu = mean(dFF_eventAll(1:wbaseline,:));
sigma_ = std(dFF_eventAll(1:wbaseline,:));
dFF_eventAll = (dFF_eventAll - miu)./sigma_;

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



