%% STATISTIC OXYTOCIN BEHAVIRO
clear all; close all;

dataname = ['F:\oxytocin\OXYT010\OXYT010_004_220126_112539',
            'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
            'F:\oxytocin\OXYT014\OXYT014_001_220429_155806',
            'F:\oxytocin\OXYT012\OXYT012_001_220323_130702',
            ]
RipCh = [23,70,19,44];Ripple_wake_ACh =[];
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);
    analogin= read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
    lfp_fbtime = find(analogin(1,:)>(max(analogin(1,:))+min(analogin(1,:)))/2,1)
    if isempty(lfp_fbtime) == 1
        lfp_fbtime = [0]
    else lfp_fbtime(1,2) = find(analogin(1,:)>0.6,1,'last');
    end 
    
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
        load([name,ext],[name,'_AChSti']);
        ACh_sti = eval([name,'_AChSti']);
        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
    
    %% get RIPPLES
    %load([basename,'.ripples.events.mat']);
    ripples = bz_FindRipples(basepath, RipCh(ii));

    Wake = SleepState.ints.WAKEstate;
    Ripple_wake =[];
    for k = 1:length(Wake)
        Rip_wake = [];
        Rip_wake = ripples.peaks(ripples.peaks>Wake(k,1)&ripples.peaks<Wake(k,2));
        Ripple_wake = [Ripple_wake;Rip_wake];
    end
    
    events = Ripple_wake- lfp_fbtime(1)/20000;
    events = events(find(events-20 > 1 & events+20 < length(dFoverF)/100));
    win = round((events-20)*100); win(:,2) = win(:,1)+4001;
    peri_dF = [];
    for kk = 1:length(win)
        peri_dF(kk,:) = dFoverF(win(kk,1):win(kk,2));
    end
    h = figure;
    options.color_area = [128 193 219]./255;    % Blue theme
    options.color_line = [ 52 148 186]./255;
    options.alpha      = 0.5;
    options.line_width = 2;
    options.error      = 'sem';
    plot_areaerrorbar(peri_dF, options);set(gcf,'Position',  [100, 100, 400, 300])
    xlim([1 4001]);
  
    Ripple_wake_ACh = [Ripple_wake_ACh;peri_dF(:,1:4001)];
end

figure;
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
plot_areaerrorbar(Ripple_wake_ACh, options);set(gcf,'Position',  [100, 100, 400, 300])
xlim([1 4001]);
