function [dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_eventAll] = getFiringRateEvents(ripRate,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(ripRate,'minpeakheight',mean(ripRate)+thresholds(1)*std(ripRate),...
    'WidthReference','halfprom','MinPeakDistance',durationPeakDiff*2);
[pkb, locs_b, ~] = findpeaks(ripRate,'minpeakheight',mean(ripRate)+thresholds(2)*std(ripRate),...
    '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-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




