clear all;clc   %% check the alighment 
%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\oxytocin\OXYT010\OXYT010_004_220126_112539',
            'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
            %'F:\oxytocin\OXYT012\OXYT012_001_220323_130702',
            'F:\oxytocin\OXYT013\OXYT013_002_220412_124441',
            'F:\oxytocin\OXYT014\OXYT014_001_220429_155806']
ACh_trans_all  = [];
for ii = 1: size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
        
    info = read_Intan_RHD2000_file(basepath);
    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    
        
    % sitmulation and signal
   if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat';
        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'));
        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
    
%    
    %% get transit from NREM to Micro
    MAfile = 'G:\My Drive\Papers_Buzsaki\Paper_ACh_OXT\Fig2\nonNREMpackets_data\OXT_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;
        if win(1)>0 & win(2)<length(dFoverF)
           ACh_trans(i,:) = dFoverF(win(1):win(2))';
        end
    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');
hold on
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_OXT.pdf','BackgroundColor','none','ContentType','vector')
