clear all;clc   %% check the alighment 
%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\LIMP013\LIMP013_009_210630_101849',
            'F:\LIMP014\LIMP014_013_211006_133329',
            'F:\LIMP014\LIMP014_014_211007_122021',
            'F:\LIMP018\LIMP018_003_220127_150338', % align with analogin
            'F:\LIMP019\LIMP019_004_220303_115500'] % align with analogin
% Ch = [57;17;17;90;112]; %% theta channel
Ch = [57,50,27,26,19];  %% Rip channel
ACh_trans_all = [];
for ii = 1: size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.SleepScoreLFP.LFP.mat']);
    load([basename,'.accel.behavior.mat']);

        
    info = read_Intan_RHD2000_file(basepath);
    if ii<4
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels);
        lfp_fbtime = find(digitalin> (min(digitalin)+max(digitalin))/2,1);lfp_fbtime(1,2) = find(digitalin> (min(digitalin)+max(digitalin))/2,1,'last');
        delete digitalin
        if ii>1 & ii< 4
            lfp_fbtime = lfp_fbtime*3;
        end
    else 
        analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        lfp_fbtime = find(analogin(1,:)> (min(analogin(1,:))+max(analogin(1,:)))/2,1);lfp_fbtime(1,2) = find(analogin(1,:)> (min(analogin(1,:))+max(analogin(1,:)))/2,1,'last');
        delete analogin   
    end
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
         if ii<4
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);
            fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>1,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>1,1,'last'));
        else 
            load([name,ext],[name,'_AChSti']);
            ACh_sti = eval([name,'_AChSti']);
            fbtime = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1));fbtime(1,2) = ACh_sti.times(find(ACh_sti.values(:,1)>0.5,1,'last'));
        end
        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,2);
        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')
    end
    
%     GammaCh = Ch(ii);
%     THLFPs = bz_GetLFP(GammaCh,'basepath',basepath,'noPrompts',true);
%     thetarange = linspace(0.5,250,1000);
%     lfp = double(THLFPs.data);
%     lfp(abs(lfp)>8000) = nan;
%     lfp = fillmissing(lfp,'linear');
%     [ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
%     [spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',2);
%     bands = SpectrogramBands(spectrogrammt,fm);
%     figure;
%     colormap('jet')
%     imagesc(tt,thetarange,abs(ss));clim([0,8e4]);
%     yticks([0.5,20,50,100,200]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
%     set(gca,'YDir','normal');
%     hold on; plot(tm,(smoothdata(bands.spindles,10))/1000+110,'color','r');xlim([0 inf]);
%     hold on; plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, dFoverF*100+100,'g');xlim([0 inf])
%     hold on; plot(accel.timestamps/3, accel.acceleration.motion/50+150,'w');xlim([0 inf]);
%     figure;plot(SleepState.idx.states);xlim([0 inf]);

    %% get transit from NREM to Micro
    MAfile = 'G:\My Drive\Papers_Buzsaki\Paper_ACh_OXT\Fig2\nonNREMpackets_data\ACh_NREM'
    load(fullfile(MAfile,[basename,'_nNREM_Data.mat']),'nNREM_Data');
    MA_AChidx = round(nNREM_Data.MA.timestamps *100 - lfp_fbtime(1)/200);
    MA_AChidx = MA_AChidx(MA_AChidx(:,1)>0 & MA_AChidx(:,2)<length(dFoverF) ,:);
    ACh_trans = [];
    for i = 1:size(MA_AChidx,1)
        win = MA_AChidx(i,1) + [-10 10]*100;
        ACh_trans(i,:) = dFoverF(win(1):win(2))';
    end

    ACh_trans_all = [ACh_trans_all;ACh_trans];
    
end
%% sort and plot
MA_duration = nNREM_Data.MA.timestamps(:,2)-nNREM_Data.MA.timestamps(:,1);
[a_sort, a_order] = sort(MA_duration);
ACh_trans_sorted = ACh_trans_all(a_order,:);
colormap('jet');
ts = 0;te = (a_sort+100)*100;


Ntrial = size(ACh_trans_sorted,1);
ave_oxt = mean(ACh_trans_sorted(:,1:500),2);
miu = std(ACh_trans_sorted(:,1:500),[],2);
ACh_trans_sorted_norm = (ACh_trans_sorted - ave_oxt)./miu;

figure;
subplot(1,2,1)
imagesc(ACh_trans_sorted_norm)
zlim = [-8 3];
caxis(zlim)
colormap(NegativeEnhancingColormap(256, zlim, ...
        [0 0 1], [1 0 0],[1 1 1],1.25));
ylabel('Cell Number');colorbar
hold on
plot(te,1:Ntrial,'-.k','LineWidth',1);
hold on;
xline(1000,'-.k')

subplot(1,2,2)
plot(1:size(ACh_trans_sorted,2),mean(smoothdata(ACh_trans_sorted_norm,2)));
xlim([1 inf])

exportgraphics(gcf,'MA_transit_ACh.pdf','BackgroundColor','none','ContentType','vector')