%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\oxytocin\OXYT010\OXYT010_004_220126_112539',
            'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
            'F:\oxytocin\OXYT013\OXYT013_001_220411_204913',
            'F:\oxytocin\OXYT012\OXYT012_001_220323_130702',
            'F:\oxytocin\OXYT014\OXYT014_001_220429_155806']
oxytREM = [];oxytWAKE = [];oxytNREM = [];
for ii = 1:5 %size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    
    info = read_Intan_RHD2000_file(basepath);
    analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
    lfp_fbtime = find(analogin(1,:)>0.6,1);lfp_fbtime(1,2) = find(analogin(1,:)>0.6,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,10);
        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
    REM = SleepState.ints.REMstate - lfp_fbtime(1)/20000; WAKE = SleepState.ints.WAKEstate - lfp_fbtime(1)/20000;NREM = SleepState.ints.NREMstate - lfp_fbtime(1)/20000;
    REM = round(REM*100);WAKE = round(WAKE*100);NREM = round(NREM*100);
    oxt_rem = []
    for jj = 1:length(REM)
        oxt_remin = [];
        oxt_remin = mean(dFoverF(REM(jj,1):REM(jj,2)));
        oxt_rem = [oxt_rem;oxt_remin];
    end
    oxt_wake = [];
    for jj = 2:length(WAKE)-1
        oxt_wakein = [];
        oxt_wakein = mean(dFoverF(WAKE(jj,1):WAKE(jj,2)));
        oxt_wake = [oxt_wake;oxt_wakein];
    end

    oxt_nrem = [];
    for jj = 2:length(NREM)-1
        oxt_nremin = [];
        oxt_nremin = mean(dFoverF(NREM(jj,1):NREM(jj,2)));
        oxt_nrem = [oxt_nrem;oxt_nremin];
    end
    oxytREM = [oxytREM;oxt_rem];
    oxytWAKE = [oxytWAKE;oxt_wake];
    oxytNREM = [oxytNREM;oxt_nrem];
end
G = [ones(size(oxytREM));2*ones(size(oxytNREM));3*ones(size(oxytWAKE))];
X = [oxytREM;oxytNREM;oxytWAKE];
figure;boxplot(X,G,'notch','on','colors',[0 0 0],'symbol','','labels',{'data1','data2','data3'}); 