
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 12);
set(0,'defaultAxesFontSize',12);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')
set(0, 'DefaultAxesFontName', 'Arial');
set(0, 'DefaultLineLineWidth', 1.5);
%%
clear all
close all
clc
fname = 'LIMP014_014_211007_122021';
load([fname,'ACh_popRate_noRip.mat']');
load([fname, '.SleepState.states.mat']);
%
N_dFF = length(dFoverF);

fs = 20e3;
fs_dff = 100;
Tstart = floor(lfp_fbtime(1)/fs);
Nstart = floor(lfp_fbtime(1)/fs*fs_dff);
Tend = ceil(lfp_fbtime(2)/fs);
Nend = ceil(lfp_fbtime(2)/fs*fs_dff);
popRateAligned = double(popRate(Nstart:Nend));

%
sleep_ints = bz_INTtoIDX(SleepState.ints,'statename',{'WAKE','','NREM','','REM'});
sleeps = sleep_ints.states;
% sleeps(sleeps ==3) =5;
% sleeps(sleeps ==2) =3;
% sleepAligned = round(interp(sleeps(Tstart:Tend),fs_dff)*2)/2;
% N_dff = length(dFoverF);
% if N_dff< length(sleepAligned)
%     sleepAligned = sleepAligned(1:N_dff);
% end
%

figure;
plot(sleeps,'or');hold on
tplot = (1:1:N_dFF)/fs_dff;
plot(tplot,zscore(dFoverF),'b')

%%

windowSize = 1;
thresholds =[0.5,1];
time_window = [200,20];
[dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_event_mean] = getFiringRateEvents(popRateAligned,dFoverF,sleeps,thresholds,time_window);
figure
t_plot = -time_window(1):1/fs_dff:time_window(2);
plot(t_plot,zscore(mean(dFF_event_mean,2)),'k')
hold on
plot(t_plot,zscore(mean(dFF_eventNREM,2)),'b')
plot(t_plot,zscore(mean(dFF_eventWAKE,2)),'g')
plot(t_plot,zscore(mean(dFF_eventREM,2)),'r')
legend('all','NREM','WAKE','REM')
%%
%%
clear all
clc
fnames = {'LIMP019_004_220303_115500';
    'LIMP018_003_220127_150338';
    'LIMP014_013_211006_133329';
    'LIMP014_014_211007_122021';};
%
%
fs = 20e3;
fs_dff = 100;
time_window = [15,10];

fig1 = figure();
set(fig1,'Position',[100 100 400 1000])

ths =[0.5,1;1,2;2,3;3,20];
alphas_ = [0.2,0.4,0.7,1.0];

fsAcc = [2,3,2,2];
for kk = 1:size(ths,1)
    dFF_NREM_all = zeros([sum(time_window*fs_dff)+1,length(fnames) ]);
    dFF_REM_all = zeros([sum(time_window*fs_dff)+1,length(fnames) ]);
    dFF_WAKE_all = zeros([sum(time_window*fs_dff)+1,length(fnames) ]);
    for i = 1:length(fnames)
        load([fnames{i},'_Ach_popRate_noRip_renew.mat']');
        load([fnames{i},'.accel.behavior.mat']');
        acc = accel.acceleration.motion;
        acc = downsample(interp(acc,fs_dff),fsAcc(i));
        dFoverF = ACh_popRate_noRip.signal;
        lfp_fbtime = ACh_popRate_noRip.lfp_fbtime;
        N_dFF = length(dFoverF);

        fs = 20e3;
        fs_dff = 100;
        Tstart = floor(lfp_fbtime(1)/fs);
        Nstart = floor(lfp_fbtime(1)/fs*fs_dff);
        Tend = floor(lfp_fbtime(2)/fs);
        Nend = floor(lfp_fbtime(2)/fs*fs_dff);
        popRate = ACh_popRate_noRip.popRate;
        popRateAligned = double(popRate(Nstart:Nend));
        
        %
        accAligned = acc(Nstart:Nend);
        %
        thresh_ = min(ACh_popRate_noRip.popNREM);
        sleeps = 3*ones(size(popRate));
        sleeps(ACh_popRate_noRip.popNREM>thresh_) =3;
        sleeps(ACh_popRate_noRip.popWAKE >thresh_) =1;
        sleeps(ACh_popRate_noRip.popREM >thresh_) =5;
        sleepAligned = sleeps(Nstart:Nend);
        
        if N_dFF< length(sleepAligned)
            sleepAligned = sleepAligned(1:N_dff);
        end

        windowSize = 1;
        thresholds = ths(kk,:);
        [dFF_eventNREM,dFF_eventREM,dFF_eventWAKE,dFF_event_mean] = getFiringRateEventsAch(popRateAligned,(dFoverF),(accAligned),sleepAligned,thresholds,time_window);

        if ~isempty(dFF_eventREM)
            dFF_REM_all (:,i)  =  mean(dFF_eventREM,2);
        end

        if ~isempty(dFF_eventNREM)
            dFF_NREM_all (:,i)  =  mean(dFF_eventNREM,2);
        end
        if ~isempty(dFF_eventWAKE)
            dFF_WAKE_all (:,i)  =  mean(dFF_eventWAKE,2);
        end
        i
    end
    %
    t_plot = -1*time_window(1):1/fs_dff:time_window(2);

    tbaseline = 5;
    subplot(2,1,1)
    dFF_event_all = dFF_NREM_all;
    meanBaseLines = squeeze(mean(dFF_event_all(1:tbaseline*fs_dff,:),1));
    dFF_plot_ACH= dFF_event_all-meanBaseLines;
    plot(t_plot,mean(dFF_plot_ACH,2),'Color',[1.0,0,0, alphas_(kk)]);
    hold on
    title('NREM');
    ave_NREM(:,kk) = mean(dFF_plot_ACH,2);


    box off
    xlim([-10,10])
    subplot(2,1,2)
    dFF_event_all = dFF_WAKE_all;
    meanBaseLines = squeeze(mean(dFF_event_all(1:tbaseline*fs_dff,:),1));
    dFF_plot_ACH= dFF_event_all-meanBaseLines;
    plot(t_plot,mean(dFF_plot_ACH,2),'Color',[1.0,0,0, alphas_(kk)]);
    hold on
    title('WAKE')
    box off
    ave_WAKE(:,kk) = mean(dFF_plot_ACH,2);

    % subplot(3,1,3)
    % dFF_event_all = dFF_REM_all ;
    % meanBaseLines = squeeze(mean(dFF_event_all(1:tbaseline*fs_dff,:),1));
    % dFF_plot_Oxt= dFF_event_all-meanBaseLines;
    % plot(t_plot,mean(dFF_plot_Oxt,2),'Color',[1.0,0, 0, alphas_(kk)]);
    % hold on
    % title('REM')
    % box off
    xlim([-10,10])

end
legend('0.5-1', '1-2', '2-3','3-max',Location='best')
%exportgraphics(gcf,'All_Ach_MUAs.pdf','BackgroundColor','none','ContentType','vector')



%%

figure;plot(sleepAligned)
hold on


yline(mean(popRateAligned)+2*std(popRateAligned))

plot(zscore(dFoverF))
plot((zscore(popRateAligned)))
plot((zscore(accAligned)))
