clear all;close all;

dataname = ['F:\oxytocin\OXYT014\OXYT014_005_220515_133849', %% good signal and good CA1/CA2/CA3  optitrack         
%             'F:\oxytocin\OXYT025\OXYT025_005_230110_110227', %% good signal and good CA1/CA2/CA3  optitrack   OXT025 NO response for OXT signa\l      
%             'F:\oxytocin\OXYT025\OXYT025_006_230116_123923', %% good signal and good CA1/CA2/CA3  optitrack         
%             'F:\oxytocin\OXYT025\OXYT025_010_230123_170927', %% good signal and good CA1/CA2/CA3  optitrack
%             'F:\oxytocin\OXYT025\OXYT025_013_230207_105757', %% good signal and good CA1/CA2/CA3  optitrack  
%             'F:\oxytocin\OXYT025\OXYT025_016_230215_163917', %% good signal and good CA3  optitrack OxtSig_H
            'F:\oxytocin\OXYT026\OXYT026_002_230103_172212', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_006_230114_141420', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_007_230117_135633', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_010_230208_180921', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_013_230220_120816', %%  good signal and good CA3  optitrack OxtSig_H
            'F:\oxytocin\OXYT027\OXYT027_002_230104_110308'] %%  good signal and good CA1/CA2/CA3  optitrack

%% inhibtion on ACh signal
for ii = 1:size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.pulses.events.mat']);
    if ii>1
        load([basename,'.accel.behavior.mat']);
    end
    
    info = read_Intan_RHD2000_file(basepath);
    if ii == 1
        analogin= read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        DD = analogin(1,:);
        lfp_fbtime = find(DD>(max(DD)+min(DD))/2,1);
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(DD>(max(DD)+min(DD))/2,1,'last');            
        end 
    else 
        digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
        lfp_fbtime = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1);
        if isempty(lfp_fbtime) == 1
            lfp_fbtime = [0]
        else lfp_fbtime(1,2) = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1,'last');
        end 
    end

    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat';
        if ii == 1;
            load([name,ext],[name,'_AChSti']);
            ACh_sti = eval([name,'_AChSti']); 
        else 
            load([name,ext],[name,'_OptoTTL']);
            ACh_sti = eval([name,'_OptoTTL']); 
        end
        otr = 0.3;
        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);
        if ii == 1
            load([name,ext],[name,'_AChSig1']);
            ACh_signal = eval([name,'_AChSig1']);
        elseif ii == 6 | ii == 11
            load([name,ext],[name,'_OxtSig_H']);
            ACh_signal = eval([name,'_OxtSig_H']);
        else 
            load([name,ext],[name,'_AChSig_H']);
            ACh_signal = eval([name,'_AChSig_H']);
        end
        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 peri_sti signal based on behavior 
    Oxt_Opto = [];
    pulse = pulses.intsPeriods-lfp_fbtime(1)/info.sample_rate;
    pul_NREM = [];
    
    if ii == 1
        SleepState.ints.NREMstate = SleepState.ints.NREMstate - lfp_fbtime(1)/info.sample_rate;
        NREM.opto = [];
        for jj = 1:length(SleepState.ints.NREMstate)    
            NREM.optoin = find(pulse(:,1)>SleepState.ints.NREMstate(jj,1)& pulse(:,2)< SleepState.ints.NREMstate(jj,2))
            NREM.opto = [NREM.opto;NREM.optoin];
        end 
        pul_NREM = pulse(NREM.opto,:);
    else    
        for jj = 1:length(pulse)
            win = round(pulse(jj,:)+lfp_fbtime(1)/info.sample_rate)+[-10 20];
            pul_in = [];
            if mean(accel.acceleration.motion(win(1):win(2)))<100
                pul_in = pulse(jj,:);
                pul_NREM = [pul_NREM;pul_in];
            end 
        end
    end

    for jj = 1:size(pul_NREM,1)
        win = round(pul_NREM(jj,1)*100) + [-10 40]*100;
        if win(1>0) & win(2)<length(dFoverF)
            Oxt_Opto(jj,:) = dFoverF(win(1):win(2));
            Oxt_Opto(jj,995:1012) = nan;
            Oxt_Opto(jj,1995:2012) = nan;
            Oxt_Opto(jj,:) = fillmissing(Oxt_Opto(jj,:),'linear');
        end
    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(Oxt_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
    xlim([1 5000]);title([basename])
%% divide by behavior
    Oxt_NREMopto  = Oxt_Opto;
    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig5\Figure 5D\Oxt_activation\Oxt_activation_signal';
    save(fullfile(datpath,[basename,'Oxt_peri_Act.mat']),'Oxt_NREMopto');
%     close all;
end
