clear all;close all;

dataname = ['F:\oxytocin\ACh_OXT004\ACh_OXT004_012_230328_141439', %% good ach and good CA1  optitrack         
           'F:\oxytocin\ACh_OXT005\ACh_OXT005_006_230309_151131', %% good ach need pick 5s stimulation optitrack   
           'F:\oxytocin\ACh_OXT007\ACh_OXT007_009_230330_104326',  %% good ach need pick 5s stimulation optitrack   
           'F:\oxytocin\ACh_OXT006\ACh_OXT006_006_230311_181030',  %% good ach need pick 5s stimulation optitrack   
           'F:\oxytocin\ACh_OXT004\ACh_OXT004_020_230516_130710',   %% NO ACH LFP need pick 5s stimulation optitrack
           'F:\oxytocin\ACh_OXT007\ACh_OXT007_016_230517_114109']  %% NO ACH LFP need pick 5s stimulation optitrack
RipCh = [58,27,27,8,59,23];
%% inhibtion on ACh signal
for ii = 1:size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.pulses.events.mat']);
    load([basename,'.accel.behavior.mat']);
    
    info = read_Intan_RHD2000_file(basepath);
    RipLFPs = bz_GetLFP(RipCh(ii),'basepath',basepath,'noPrompts',true);
    RipLFPs.data = double(RipLFPs.data);
    RipLFPs.data(abs(RipLFPs.data)>5000) = nan;
    RipLFPs.data = fillmissing(RipLFPs.data,'linear');

    pulse = pulses.intsPeriods;
    pulses.duration = pulses.intsPeriods(:,2) - pulses.intsPeriods(:,1);
    idx = find(pulses.duration>4 & pulses.duration<6);
    pulse = pulse(idx,:)
    %% exclude the move-pulse
  
    pul_NREM = [];
    for jj = 1:length(pulse)
        win = round(pulse(jj,:))+[-5 10];
        pul_in = [];
        if mean(accel.acceleration.motion(win(1):win(2)))<50
            pul_in = pulse(jj,:);
            pul_NREM = [pul_NREM;pul_in];
        end 
    end
   
    %% get LFP and ripples
    delete([basename,'.ripples.events.mat']);
    delete([basename,'.EMGFromLFP.LFP.mat']);
    
    RipLFPs = bz_GetLFP(RipCh(ii),'basepath',basepath,'noPrompts',true);
    RipLFPs.data = double(RipLFPs.data);
    RipLFPs.data(abs(RipLFPs.data)>5000) = nan;
    RipLFPs.data = fillmissing(RipLFPs.data,'linear');

    % get accel into 1250hz
    x = [1:length(accel.acceleration.motion)];
    xq = [1:1/1250:length(accel.acceleration.motion)];
    Mov_time =  interp1(x,accel.acceleration.motion,xq);

    RipLFPs.data(Mov_time>50) = 0;

    ripples = bz_FindRipples(RipLFPs.data,RipLFPs.timestamps);
    
    %% get peri_ripples stas

    if ii == 1;
        pulse = pul_NREM-10;
    else
        pulse = pul_NREM;
    end
    ripple_optoC = [];ripple_optoS = [];ripple_optoP = [];
    ripple_optoCC = [];ripple_optoSS = [];ripple_optoPP = [];
    rippleCC = {};rippleSS = {};ripplePP = {}
    pulse = pul_NREM;
    for jj = 1:length(pulse)
        ripple_optoinC = find(ripples.peaks > (pulse(jj,1)-10) & ripples.peaks < pulse(jj,1))
        ripple_optoC = [ripple_optoC;ripple_optoinC];
        ripple_optoinS = find(ripples.peaks > pulse(jj,1) & ripples.peaks < pulse(jj,2))
        ripple_optoS = [ripple_optoS;ripple_optoinS];
        ripple_optoinP = find(ripples.peaks > pulse(jj,2) & ripples.peaks < (pulse(jj,2)+10))
        ripple_optoP = [ripple_optoP;ripple_optoinP];
    
        ripple_optoCC(jj,1) = length(find(ripples.peaks > (pulse(jj,1)-10) & ripples.peaks < pulse(jj,1)));
        ripple_optoSS(jj,1) = length(find(ripples.peaks > pulse(jj,1) & ripples.peaks < pulse(jj,2)));
        ripple_optoPP(jj,1) = length(find(ripples.peaks > pulse(jj,2) & ripples.peaks < (pulse(jj,2)+10)));
        
        rippleCC{jj,1} = ripples.peaks(find(ripples.peaks > (pulse(jj,1)-10) & ripples.peaks < pulse(jj,1)))
        IRI_C{jj,1} = diff(rippleCC{jj,1});
        rippleSS{jj,1} = ripples.peaks(find(ripples.peaks > pulse(jj,1) & ripples.peaks < pulse(jj,2)))
        IRI_S{jj,1} = diff(rippleSS{jj,1});
        ripplePP{jj,1} = ripples.peaks(find(ripples.peaks > pulse(jj,2) & ripples.peaks < (pulse(jj,2)+10)))
        IRI_P{jj,1} = diff(ripplePP{jj,1});
    end
    
    
    rippleDurC = ripples.timestamps(ripple_optoC,2)-ripples.timestamps(ripple_optoC,1);
    rippleDurS = ripples.timestamps(ripple_optoS,2)-ripples.timestamps(ripple_optoS,1);
    rippleDurP = ripples.timestamps(ripple_optoP,2)-ripples.timestamps(ripple_optoP,1);
    
    ripplePeakPoweC = ripples.peakNormedPower(ripple_optoC);
    ripplePeakPoweS = ripples.peakNormedPower(ripple_optoS);
    ripplePeakPoweP = ripples.peakNormedPower(ripple_optoP);
    
    aa = [];
    aa(:,1) = ripple_optoCC;
    aa(:,2) = ripple_optoSS;
    aa(:,3) = ripple_optoPP;
    figure; boxchart(aa)

    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig5\Figure 5C\ACh_activation\ACh_activation_Ripples';
    save(fullfile(datpath,[basename,'_ripple_opto.mat']),'ripple_optoCC','ripple_optoSS','ripple_optoPP');
   
end
