clear all;close all;
%% data no optitrack
dataname = ['F:\oxytocin\OXYT010\OXYT010_004_220126_112539',
            'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
            'F:\oxytocin\OXYT013\OXYT013_002_220412_124441',
            'F:\oxytocin\OXYT014\OXYT014_001_220429_155806']


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,:)>(max(analogin(1,:))+min(analogin(1,:)))/2,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

    if exist([basename,'.ripples.events.mat'])
        load([basename,'.ripples.events.mat']);
    else
        ripples = bz_FindRipples(basepath,RipCh);
    end
    %% get CA1 popRate
    load([basename,'.cell_metrics.cellinfo.mat']);
    load([basename,'.spikes.cellinfo.mat']);
    

    binSize = 0.01;
    spkMat = bz_SpktToSpkmat(spikes,'binSize',binSize,'overlap',1);
    smoothingWindow = 10;

    spkMat = bz_SpktToSpkmat(spikes,'binSize',binSize,'overlap',1);
    types = convertCharsToStrings(cell_metrics.putativeCellType);
    spkMat_CA1 = spkMat.data(:,cell_metrics.tags.CA1);


    data = spkMat_CA1;
    popRate = zeros(length(data),1);
    Nber = size(data,2);
    for spk = 1:Nber
        popRate = popRate + zscore(Smooth(data(:,spk),smoothingWindow));
    end
    popRate = popRate./spk;
    popRate = popRate + abs(min(popRate));
    % popRate = smoothdata(popRate,'gaussian',250);
    IDX = round(ripples.timestamps*100)
    for kk = 1:length(ripples.peaks)
        popRate(IDX(kk,1):IDX(kk,2)) = min(popRate);
    end

    h = figure();
    plot(spkMat.timestamps,popRate);
    tim = [1:length(dFoverF)]/100;
    hold on;plot(tim+lfp_fbtime(1)/2e4, dFoverF*10);xlim([0 inf]);

    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig4\Figure 4C\Oxt_popRate_noRip';
    save(fullfile(datpath,[basename,'Oxt_popRate_noRip.mat']),'dFoverF','popRate','lfp_fbtime')
    close all
    
end