clear all;close all;

dataname = ['F:\LIMP003\LIMP003_006_210413_112545', %% good ach and good CA1  optitrack         
           'F:\LIMP018\LIMP018_003_220127_150338', %% good ach and good CA1  
           'F:\LIMP011\LIMP011_008_210622_092807', %% good ach-signal but CA3 neurons   optitrack
           'F:\LIMP013\LIMP013_011_210630_142011', %% good ach-signal but CA3 neurons   optitrack
           'F:\LIMP016\LIMP016_005_220103_130720', %% need correction of signal but good CA1 neurons
           'F:\LIMP019\LIMP019_007_220307_155644'] %% no ach-signal but has CA1 neurons;

%% inhibtion on ACh signal
for ii = 1:4
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.pulses.events.mat'])
    
    if ii~=2
        info = read_Intan_RHD2000_file(basepath);
        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 
    elseif ii == 2
        info = read_Intan_RHD2000_file(basepath);
        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 
    end
    % sitmulation and signal
    if exist('dFoverF.mat')==0
        name = basename(1:11);
        ext = '.mat'
        if ii~=2
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);          
        elseif  ii == 2
            load([name,ext],[name,'_AChSti']);
            ACh_sti = eval([name,'_AChSti']);
        end
        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 peri_sti signal
    ACh_Opto = [];
    if ii~=4;
        aa= 0;
    elseif ii ==4
        aa = 4;
    end
    pulse = pulses.intsPeriods-lfp_fbtime(1)/info.sample_rate+aa;
    for jj = 1:length(pulse)
        win = round(pulse(jj,1)*100) + [-5 10]*100;
        if win(1>0) & win(2)<length(dFoverF)
            ACh_Opto(jj,:) = dFoverF(win(1):win(2));
            ACh_Opto(jj,495:545) = nan;
            ACh_Opto(jj,995:1012) = nan;
            ACh_Opto(jj,:) = fillmissing(ACh_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(ACh_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
    xlim([1 1500]);
%% divide by behavior
    SleepState.ints.REMstate = SleepState.ints.REMstate - lfp_fbtime(1)/info.sample_rate+aa;
    SleepState.ints.NREMstate = SleepState.ints.NREMstate - lfp_fbtime(1)/info.sample_rate+aa;
    SleepState.ints.WAKEstate = SleepState.ints.WAKEstate - lfp_fbtime(1)/info.sample_rate+aa;
    REM.opto = [];NREM.opto = [];WAKE.opto = [];
    for jj = 1:length(SleepState.ints.REMstate)
        REM.optoin = find(pulse(:,1)>SleepState.ints.REMstate(jj,1)& pulse(:,2)< SleepState.ints.REMstate(jj,2))
        REM.opto = [REM.opto;REM.optoin];
    end
    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
    for jj = 1:length(SleepState.ints.WAKEstate)
        WAKE.optoin = find(pulse(:,1)>SleepState.ints.WAKEstate(jj,1)& pulse(:,2)< SleepState.ints.WAKEstate(jj,2))
        WAKE.opto = [WAKE.opto;WAKE.optoin]
    end
    ACh_REMopto = ACh_Opto(REM.opto,:);
    ACh_NREMopto = ACh_Opto(NREM.opto,:);
    ACh_WAKEopto = ACh_Opto(WAKE.opto,:);

    g = figure(2);
    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(ACh_REMopto, options);set(gcf,'Position',  [100, 100, 400, 300])
    xlim([0 1500]);title('REM')
    
    o = figure(3);
    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(ACh_NREMopto, options);set(gcf,'Position',  [100, 100, 400, 300])
    xlim([0 1500]);title('NREM')
    
    k = figure(4);
    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(ACh_WAKEopto, options);set(gcf,'Position',  [100, 100, 400, 300])
    xlim([0 1500]);title('WAKE')

    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig5\Figure 5B';
    save(fullfile(datpath,[basename,'ACh_peri_Inh.mat']),'ACh_Opto','ACh_REMopto','ACh_NREMopto','ACh_WAKEopto');
    close all;
end
