clear all;close all;

dataname = ['F:\oxytocin\OXYT014\OXYT014_005_220515_133849', %% good signal and good CA1/CA2/CA3  optitrack         
%             'F:\oxytocin\OXYT025\OXYT025_005_230110_110227', %% good signal and good CA1/CA2/CA3  optitrack   OXT025 NO response for OXT signa\l      
%             'F:\oxytocin\OXYT025\OXYT025_006_230116_123923', %% good signal and good CA1/CA2/CA3  optitrack         
%             'F:\oxytocin\OXYT025\OXYT025_010_230123_170927', %% good signal and good CA1/CA2/CA3  optitrack
%             'F:\oxytocin\OXYT025\OXYT025_013_230207_105757', %% good signal and good CA1/CA2/CA3  optitrack  
%             'F:\oxytocin\OXYT025\OXYT025_016_230215_163917', %% good signal and good CA3  optitrack OxtSig_H
            'F:\oxytocin\OXYT026\OXYT026_002_230103_172212', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_006_230114_141420', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_007_230117_135633', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_010_230208_180921', %%  good signal and good CA1/CA2/CA3  optitrack
            'F:\oxytocin\OXYT026\OXYT026_013_230220_120816', %%  good signal and good CA3  optitrack OxtSig_H
            'F:\oxytocin\OXYT027\OXYT027_002_230104_110308'] %%  good signal and good CA1/CA2/CA3  optitrack

RipCh = [93,8,8,27,46,42,11];
%% 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,'.SleepState.states.mat']);
    if ii>1
        load([basename,'.accel.behavior.mat']);
    end
    
    
    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
  
    pulse = pulses.intsPeriods;
    pul_NREM = [];
    
    if ii == 1
        SleepState.ints.NREMstate = SleepState.ints.NREMstate;
        NREM.opto = [];
        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 
        pul_NREM = pulse(NREM.opto,:);
    else    
        for jj = 1:length(pulse)
            win = round(pulse(jj,:))+[-10 20];
            pul_in = [];
            if mean(accel.acceleration.motion(win(1):win(2)))<100
                pul_in = pulse(jj,:);
                pul_NREM = [pul_NREM;pul_in];
            end 
        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 or NREM into 1250hz
    if ii == 1
        x = [1:length(SleepState.idx.states)];
        xq = [1:1/1250:length(SleepState.idx.states)];
        Mov_time =  interp1(x,SleepState.idx.states,xq);
        RipLFPs.data(Mov_time~=3) = 0;
    else 
        x = [1:length(accel.acceleration.motion)];
        xq = [1:1/625:length(accel.acceleration.motion)];
        Mov_time =  interp1(x,accel.acceleration.motion,xq);
        RipLFPs.data(Mov_time>200) = 0;
    end

    ripples = bz_FindRipples(RipLFPs.data,RipLFPs.timestamps);
%% get opto-ripple
    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)+10 & ripples.peaks < (pulse(jj,2)+20))
        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)+10 & ripples.peaks < (pulse(jj,2)+20)));
        
        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)+10 & ripples.peaks < (pulse(jj,2)+20)))
        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 5D\Oxt_activation\Oxt_activation_ripples';
    save(fullfile(datpath,[basename,'_ripple_opto.mat']),'ripple_optoCC','ripple_optoSS','ripple_optoPP');
end
