clear all;clc
dataname = ['F:\LIMP003\LIMP003_007_210420_103320'  %%corrected   optitrack shank1-4 CA1
            'F:\LIMP014\LIMP014_012_211005_112625',    %%corrected   optitrack shank1-4 CA1
            'F:\LIMP014\LIMP014_014_211007_122021',   %%corrected   optitrack shank1-4 CA1 
            'F:\LIMP014\LIMP014_013_211006_133329',   %%corrected   optitrack shank1-4 CA1 
            'F:\LIMP018\LIMP018_003_220127_150338',    %%corrected  all shanks CA1 
            'F:\LIMP019\LIMP019_004_220303_115500']    %%corrected  shank1-4 CA1
RipCh = [58,57,57,57, 26,1]


for kk = 5:6 %1:size(dataname,1)
    cd(dataname(kk,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.accel.behavior.mat'])

    %%%%%%%%%%%%%%%%%% delete short periods less than windows
    idx = diff(SleepState.ints.NREMstate,1,2)<30;
    SleepState.ints.NREMstate(idx,:) = [];
    idx = diff(SleepState.ints.REMstate,1,2)<30;
    SleepState.ints.REMstate(idx,:) = [];
    idx = diff(SleepState.ints.WAKEstate,1,2)<30;
    SleepState.ints.WAKEstate(idx,:) = [];
    sleepstates = bz_INTtoIDX(SleepState.ints);
   
    %%%%%%%%%%%%%%%%% get alignment time
    if kk>0 & kk<5
        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 kk>4
        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 kk>0 & kk<5
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);          
        elseif  kk>4
            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 ripples

    if exist([basename,'.ripples.events.mat'])
        load([basename,'.ripples.events.mat']);
    else
        ripples = bz_FindRipples(basepath,RipCh);
    end
   
    %% get CA1 popRate
    load([basename,'.spikes.cellinfo.mat']);
%     load([basename,'.spikes.cellinfo.mat']);
    spikes = spikes.times;
    
    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;
    

    
    Spikes_CA1 = spkMat.data;
    

    data = Spikes_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',500);

%     IDX = round(ripples.timestamps*100)
%     for ii = 1:length(ripples.peaks)
%         popRate(IDX(ii,1):IDX(ii,2)) = min(popRate);
%     end



    h = figure();
    plot([1:length(popRate)]/100,popRate);
    tim = [1:length(dFoverF)]/100;
    hold on;plot(tim+lfp_fbtime(1)/2e4, dFoverF*3);xlim([0 inf]);
    hold on;plot(sleepstates.timestamps, sleepstates.states);ylim([-1 3.5]);
    % hold on;plot(accel.timestamps/2,accel.acceleration.motion/1000-1)


    %% find peaks 
    threshold =[1,2,3];
    alpha_val = linspace(0.2,1.0, length(threshold));
    for ii = 1:length(threshold)
        durations = 10;
        [pks locs width] = findpeaks(popRate,'minpeakheight',mean(popRate)+threshold(ii)*std(popRate),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        figure;findpeaks(popRate,'minpeakheight',mean(popRate)+threshold(ii)*std(popRate),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        hold on; plot([1:length(dFoverF)]+lfp_fbtime(1)/2e2, dFoverF);
        hold on; plot(sleepstates.timestamps*100, sleepstates.states);ylim([-1 3.5])
        %get peri-popBurst 
        events = round(locs+[-10 15]*100-lfp_fbtime(1)/2e2);
        events = events(events(:,1)>0,:);events = events(events(:,2)<length(dFoverF),:); 
    
        dF_pop = [];
        for jj = 1:length(events)
           dF_pop(:,jj) = dFoverF(events(jj,1):events(jj,2));
        end
        miu = mean(dF_pop(1:1000,:));
        sigma = std(dF_pop(1:1000,:));
        dF_pop_norm = (dF_pop - miu)./sigma;
        p = figure;%hold on;
        p(ii) = plot(mean(dF_pop_norm,2), 'k');
        p(ii).Color = [p(ii).Color alpha_val(ii)];
        %plot_areaerrorbar(dF_pop_norm', options);set(gcf,'Position',  [100, 100, 400, 300])
        title('oxt');xlim([1 inf]);
    end
    legendCell = strcat('std=',string(num2cell(threshold)));
    legend(legendCell)
    yline(0,'k--')
    xline(1000,'k--')
    %% find peaks based on behavior states
    % select popRate besed on diffirent behavior, have to double check the
    % idx
   
    FR = popRate;
    statesIDX = interp1([1:length(sleepstates.states)],sleepstates.states,[1:0.01:length(sleepstates.states)]);
    FR(statesIDX ~= 1) = min(popRate);popWAKE = FR;
    FR = popRate;
    FR(statesIDX ~= 2) = min(popRate);popNREM = FR;
    FR = popRate;
    FR(statesIDX ~= 3) = min(popRate);popREM = FR;
    %%%%%%%%%%%%%%%%%%%%%%%

    threshold =[1,2,3];
    alpha_val = linspace(0.2,1.0, length(threshold));
    %%%%% for NREM
    for ii = 1:length(threshold)
        durations = 10;
        [pks locs width] = findpeaks(popNREM,'minpeakheight',mean(popNREM)+threshold(ii)*std(popNREM),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        figure;findpeaks(popNREM,'minpeakheight',mean(popNREM)+threshold(ii)*std(popNREM),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        hold on; plot([1:length(dFoverF)]+lfp_fbtime(1)/2e2, dFoverF);
        hold on; plot(sleepstates.timestamps*100, sleepstates.states);ylim([-1 3.5])
        %get peri-popBurst 
        events = round(locs+[-10 15]*100-lfp_fbtime(1)/2e2);
        events = events(events(:,1)>0,:);events = events(events(:,2)<length(dFoverF),:); 
    
        dF_pop = [];
        for jj = 1:length(events)
           dF_pop(:,jj) = dFoverF(events(jj,1):events(jj,2));
        end
        miu = mean(dF_pop(1:1000,:));
        sigma = std(dF_pop(1:1000,:));
        dF_pop_norm = (dF_pop - miu)./sigma;
        p = figure; %hold on;
        p(ii) = plot(mean(dF_pop_norm,2), 'k');
        p(ii).Color = [p(ii).Color alpha_val(ii)];
        %plot_areaerrorbar(dF_pop_norm', options);set(gcf,'Position',  [100, 100, 400, 300])
        title('oxt');xlim([1 inf]);
    end
    legendCell = strcat('std=',string(num2cell(threshold)));
    legend(legendCell)
    yline(0,'k--')
    xline(1000,'k--')

    %%%%% for WAKE
    for ii = 1:length(threshold)
        durations = 10;
        [pks locs width] = findpeaks(popWAKE,'minpeakheight',mean(popWAKE)+threshold(ii)*std(popWAKE),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        figure;findpeaks(popWAKE,'minpeakheight',mean(popWAKE)+threshold(ii)*std(popWAKE),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        hold on; plot([1:length(dFoverF)]+lfp_fbtime(1)/2e2, dFoverF);
        hold on; plot(sleepstates.timestamps*100, sleepstates.states);ylim([-1 3.5])
        %get peri-popBurst 
        events = round(locs+[-10 15]*100-lfp_fbtime(1)/2e2);
        events = events(events(:,1)>0,:);events = events(events(:,2)<length(dFoverF),:); 
    
        dF_pop = [];
        for jj = 1:length(events)
           dF_pop(:,jj) = dFoverF(events(jj,1):events(jj,2));
        end
        miu = mean(dF_pop(1:1000,:));
        sigma = std(dF_pop(1:1000,:));
        dF_pop_norm = (dF_pop - miu)./sigma;
        p = figure; %hold on;
        p(ii) = plot(mean(dF_pop_norm,2), 'k');
        p(ii).Color = [p(ii).Color alpha_val(ii)];
        %plot_areaerrorbar(dF_pop_norm', options);set(gcf,'Position',  [100, 100, 400, 300])
        title('oxt');xlim([1 inf]);
    end
    legendCell = strcat('std=',string(num2cell(threshold)));
    legend(legendCell)
    yline(0,'k--')
    xline(1000,'k--')

    %%%%% for REM
    for ii = 1:length(threshold)
        durations = 10;
        [pks locs width] = findpeaks(popREM,'minpeakheight',mean(popREM)+threshold(ii)*std(popREM),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        figure;findpeaks(popREM,'minpeakheight',mean(popREM)+threshold(ii)*std(popREM),...
            'WidthReference','halfprom','MinPeakDistance',durations(1)*2);
        hold on; plot([1:length(dFoverF)]+lfp_fbtime(1)/2e2, dFoverF);
        hold on; plot(sleepstates.timestamps*100, sleepstates.states);ylim([-1 3.5])
        %get peri-popBurst 
        events = round(locs+[-10 15]*100-lfp_fbtime(1)/2e2);
        events = events(events(:,1)>0,:);events = events(events(:,2)<length(dFoverF),:); 
    
        dF_pop = [];
        for jj = 1:length(events)
           dF_pop(:,jj) = dFoverF(events(jj,1):events(jj,2));
        end
        miu = mean(dF_pop(1:1000,:));
        sigma = std(dF_pop(1:1000,:));
        dF_pop_norm = (dF_pop - miu)./sigma;
        p = figure; %hold on;
        p(ii) = plot(mean(dF_pop_norm,2), 'k');
        p(ii).Color = [p(ii).Color alpha_val(ii)];
        %plot_areaerrorbar(dF_pop_norm', options);set(gcf,'Position',  [100, 100, 400, 300])
        title('oxt');xlim([1 inf]);
    end
    legendCell = strcat('std=',string(num2cell(threshold)));
    legend(legendCell)
    yline(0,'k--')
    xline(1000,'k--')
%% save data
    datpath = 'G:\My Drive\Papers_Buzsaki\Paper_ACh_OXT\Fig4\Figure 4C\ACh_popRate_noRip';
    ACh_popRate_noRip = {};
    ACh_popRate_noRip.signal = dFoverF;
    ACh_popRate_noRip.lfp_fbtime =  lfp_fbtime;
    ACh_popRate_noRip.popRate = popRate;
    ACh_popRate_noRip.popNREM = popNREM;
    ACh_popRate_noRip.popREM = popREM;
    ACh_popRate_noRip.popWAKE = popWAKE;
    save(fullfile(datpath,[basename,'_ACh_popRate_noRip.mat']),'ACh_popRate_noRip')
    close all
    
end