%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\LAMP2\LAMP2_002_190405_154550']  % corrected sleepscore corrected
dataname = ['F:\LAMP3\LAMP3_006_190613_194519']  % corrected
dataname = [%'G:\data\LAMP37\LAMP37_005_201111_204459',
            'F:\LAMP37\LAMP37_006_201201_192319'] % corrected
dataname = [%'F:\LIMP009\LIMP009_002_210714_091509',
            'F:\LIMP011\LIMP011_014_210723_121540',    %%correct sleepSCORE, REM of ACh not good optitrack 
            'F:\LIMP013\LIMP013_009_210630_101849',    %%corrected sleepSCORE  optitrack
            'F:\LIMP014\LIMP014_012_211005_112625',    %%corrected sleepSCORE   optitrack
            'F:\LIMP014\LIMP014_014_211007_122021',    %%corrected sleepSCORE   optitrack
            'F:\LIMP016\LIMP016_004_220101_144210',    %%corrected sleepSCORE   analogCh  = 2
            'F:\LIMP018\LIMP018_003_220127_150338',    %%corrected sleepSCORE  exclude data from 5700 to 11250
            'F:\LIMP019\LIMP019_004_220303_115500']    %%corrected
dF_NR_tranAll = [];dF_NW_tranAll = [];dF_RW_tranAll = [];
for ii = 1:size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    
    if ii>0 & ii<5
        info = read_Intan_RHD2000_file(basepath);
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
        lfp_fbtime = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1);
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(digitalin(1,:)>0.6,1,'last');
        end 
    elseif ii>4
        info = read_Intan_RHD2000_file(basepath);
        analogin= read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        if ii == 5
            DD = analogin(2,:);
        else ii>5
            DD = analogin(1,:);
        end
        lfp_fbtime = find(DD>(max(DD)+min(DD))/2,1);
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(DD>(max(DD)+min(DD))/2,1,'last');            
        end 
    end
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
        if ii>0 & ii<5
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);
        elseif  ii>4
            load([name,ext],[name,'_AChSti']);
            ACh_sti = eval([name,'_AChSti']);
        end
        otr = 0.5;
        fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>otr,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>otr,1,'last'));
        fbIDX = round(fbtime*20000);
        load([name,ext],[name,'_AChSig1']);
        ACh_signal = eval([name,'_AChSig1']);
        ACh_signal.data = downsample(ACh_signal.values((fbIDX(1,1)):(fbIDX(1,2))),200);
        ACh_signal.time = downsample(ACh_signal.times((fbIDX(1,1)):(fbIDX(1,2))),200)- fbtime(1,1);
        p = [];
        p = polyfit(ACh_signal.time,ACh_signal.data,5);
        Ybaseline = polyval(p,ACh_signal.time);
        figure
        plot(ACh_signal.time, ACh_signal.data); hold on;
        plot(ACh_signal.time,Ybaseline);
        dFoverF = smoothdata((ACh_signal.data - Ybaseline)./Ybaseline,'movmean',10);
        hold on;
        plot(ACh_signal.time, dFoverF);xlim([0 inf]);
    else load('dFoverF.mat'); figure;plot(dFoverF);
    end
    
    %% find transit index over 20s 
    sleep_ints = bz_INTtoIDX(SleepState.ints,statenames ={'WAKE','','NREM','','REM'});
    idx = find(ischange(sleep_ints.states)==1) ;
    aa = diff(idx);idx_ff =[];
    for k = 1:length(aa)-1
        if aa(k)>20 & aa(k+1)>20
            ff = idx(k+1);
            idx_ff = [idx_ff;ff]
        end
    end
    idx_ff = idx_ff(idx_ff < length(dFoverF)/100 + lfp_fbtime(1)/20000 & idx_ff > lfp_fbtime(1)/20000);
    %% seperate into NREM-REM, NREM-WAKE, REM-WAKE
    NR_tim = []; NW_tim = []; RW_tim = [];
    for j = 1:length(idx_ff)
        if sleep_ints.states(idx_ff(j)-1) == 3  & sleep_ints.states(idx_ff(j)) == 5
            NR_tim = [NR_tim;idx_ff(j)];
        elseif sleep_ints.states(idx_ff(j)-1) == 3  & sleep_ints.states(idx_ff(j)) == 1
            NW_tim = [NW_tim;idx_ff(j)];
        else sleep_ints.states(idx_ff(j)-1) == 5  & sleep_ints.states(idx_ff(j)) == 1
            RW_tim = [RW_tim;idx_ff(j)];
        end
    end
    %% get NREM-REM DATA

    NR_tim_dff = NR_tim -lfp_fbtime(1)/20000;
    dF_NR_transit = [];
    for jj = 1:length(NR_tim_dff)
        idx = round((NR_tim_dff(jj) + [-20 20])*100);
        dF_NR_transit(:,jj) = dFoverF(idx(1):idx(2));
    end
    dF_NR_transit = dF_NR_transit';
    dF_NR_tranAll = [dF_NR_tranAll;dF_NR_transit];
    
   %% get NREM-WAKE
    NW_tim_dff= NW_tim -lfp_fbtime(1)/20000;
    dF_NW_transit = [];
    for jj = 1:length(NW_tim_dff)
        idx = round((NW_tim_dff(jj) + [-20 20])*100);
        dF_NW_transit(:,jj) = dFoverF(idx(1):idx(2));
    end
    dF_NW_transit = dF_NW_transit';
    dF_NW_tranAll = [dF_NW_tranAll;dF_NW_transit];
  %% get REM-WAKE
    RW_tim_dff= RW_tim -lfp_fbtime(1)/20000;
    dF_RW_transit = [];
    for jj = 1:length(RW_tim_dff)
        idx = round((RW_tim_dff(jj) + [-20 20])*100);
        dF_RW_transit(:,jj) = dFoverF(idx(1):idx(2));
    end
    dF_RW_transit = dF_RW_transit';
    dF_RW_tranAll = [dF_RW_tranAll;dF_RW_transit];

end
save('G:\My Drive\Paper_ACh_OXT\Fig2\Data\Figure2B\ACh_behaviorTrans.mat','dF_NR_tranAll','dF_NW_tranAll','dF_RW_tranAll')


figure;
outsr = 100;
imagesc([1:size(dF_NR_tranAll,2)./outsr]'-20,[1:size(dF_NR_tranAll,1)]',...
    zscore(dF_NR_tranAll,[],2));
xlabel('Times from transit (s)'); ylabel('TRANSIT rumber #'); title('NREM-REM')

figure;
outsr = 100;
imagesc([1:size(dF_NW_tranAll,2)./outsr]'-20,[1:size(dF_NW_tranAll,1)]',...
    zscore(dF_NW_tranAll,[],2));
xlabel('Times from transit (s)'); ylabel('TRANSIT rumber #');title('NREM-WAKE')

figure;
outsr = 100;
imagesc([1:size(dF_RW_tranAll,2)./outsr]'-20,[1:size(dF_RW_tranAll,1)]',...
    zscore(dF_RW_tranAll,[],2));
xlabel('Times from transit (s)'); ylabel('TRANSIT rumber #');title('REM-WAKE')