%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\LAMP2\LAMP2_002_190405_154550',]
dataname = ['F:\LAMP37\LAMP37_005_201111_204459',
            'F:\LAMP37\LAMP37_006_201201_192319']
dataname = [ % 'F:\LIMP009\LIMP009_002_210714_091509',
            'F:\LIMP011\LIMP011_010_210705_120340',
            'F:\LIMP013\LIMP013_009_210630_101849',
            'F:\LIMP014\LIMP014_012_211005_112625',
            'F:\LIMP016\LIMP016_004_220101_144210',
            'F:\LIMP019\LIMP019_004_220303_115500']
oxytREM = [];oxytWAKE = [];oxytNREM = [];
for ii = 1:size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.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(1,:)>0.5,1);
%         figure;plot(digitalin);
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(digitalin(1,:)>0.6,1,'last');
        end      
        delete digitalin
    else analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
         lfp_fbtime = find(analogin(1,:)>0.3,1);
         delete analogin
    end
    % sitmulation and signal
    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]);
    
    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)
        if REM(jj,1) > 0 & REM(jj,2) < length(dFoverF)
            oxt_remin = [];
            oxt_remin = mean(dFoverF(REM(jj,1):REM(jj,2)));
            oxt_rem = [oxt_rem;oxt_remin];
        end
    end
    oxt_wake = [];
    for jj = 2:length(WAKE)-1
        if WAKE(jj,1) > 0& WAKE(jj,2) < length(dFoverF)
            oxt_wakein = [];
            oxt_wakein = mean(dFoverF(WAKE(jj,1):WAKE(jj,2)));
            oxt_wake = [oxt_wake;oxt_wakein];
        end
    end

    oxt_nrem = [];
    for jj = 2:length(NREM)-1
        if NREM(jj,1) > 0 & NREM(jj,2) < length(dFoverF)
            oxt_nremin = [];
            oxt_nremin = mean(dFoverF(NREM(jj,1):NREM(jj,2)));
            oxt_nrem = [oxt_nrem;oxt_nremin];
        end
    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'}); 