dataname = ['F:\oxytocin\ACh_OXT004\ACh_OXT004_007_230306_152436', %% corrected
            'F:\oxytocin\ACh_OXT004\ACh_OXT004_011_230325_130235',   %% corrected
            'F:\oxytocin\ACh_OXT007\ACh_OXT007_015_230507_172920',   %% corrected
            'F:\oxytocin\ACh_OXT004\ACh_OXT004_019_230507_140912']   %% corrected no REM


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);
    digitalin = read_Intan_digitalin_file('digitalin.dat',info.num_board_dig_in_channels)';
    lfp_fbtime(1) = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1);lfp_fbtime(2) = find(digitalin(1,:)>(max(digitalin)+min(digitalin))/2,1,'last');
    delete digitalin
    % sitmulation and signal
    if exist([basename,'_Fluo.mat'])==0
        FluoData = Get_OxtAch_Fluos();
    else load([basename,'_Fluo.mat'])
    end

    figure
    plot(FluoData.time,FluoData.Ach_dFF); hold on;
    plot(FluoData.time,FluoData.Oxt_dFF);xlim([0 inf]);
    %% computing the pop_rate
    load([basename,'.cell_metrics.cellinfo.mat']);
    load([basename,'.spikes.cellinfo.mat']);
    
    SWS =  round (SleepState.ints.NREMstate*100);
    WAKE =  round (SleepState.ints.WAKEstate*100);

    binSize = 0.01;
    spkMat = bz_SpktToSpkmat(spikes,'binSize',binSize,'overlap',1);
    smoothingWindow = 10;
    
    types = convertCharsToStrings(cell_metrics.putativeCellType);
    PyrN = find(strcmp(types,'Pyramidal Cell'));
    WIntN = find(strcmp(types,'Wide Interneuron'));
    NIntN = find(strcmp(types,'Narrow Interneuron'));

    % save memory and add to single vector...
    data = spkMat.data(:,PyrN);
    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);
    h = figure();
    plot(spkMat.timestamps,popRate/5);
    hold on;plot(FluoData.time+lfp_fbtime(1)/2e4, FluoData.Ach_dFF*10);xlim([0 inf]);
    hold on;plot(FluoData.time+lfp_fbtime(1)/2e4, FluoData.Oxt_dFF*10+0.2);xlim([0 inf]);

    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig2\Data\Figure2DEF';
    Fluo_dF.Ach = FluoData.Ach_dFF;
    Fluo_dF.Oxt = FluoData.Oxt_dFF;
    save(fullfile(datpath,[basename,'ACh_OXT_popRate.mat']),'Fluo_dF','popRate','lfp_fbtime')
    close all
end
