clear all;close all;

dataname = ['F:\LIMP003\LIMP003_006_210413_112545', %% good ach and good CA1  optitrack         
           'F:\LIMP018\LIMP018_003_220127_150338', %% good ach and good CA1  
           'F:\LIMP011\LIMP011_008_210622_092807', %% good ach-signal but CA3 neurons   optitrack
           'F:\LIMP013\LIMP013_011_210630_142011', %% good ach-signal but CA3 neurons   optitrack
           'F:\LIMP016\LIMP016_005_220103_130720', %% need correction of signal but good CA1 neurons
           'F:\LIMP019\LIMP019_007_220307_155644'] %% no ach-signal but has CA1 neurons;

RipCh = [58,26,nan,nan,53,1]
%% inhibtion on ACh signal
for ii = [1,2,5,6]
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.pulses.events.mat']);
    load([basename,'.accel.behavior.mat'])
    
    pulse = pulses.intsPeriods;
    %% get new 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>25) = 0;

    ripples = bz_FindRipples(RipLFPs.data,RipLFPs.timestamps);

    ripple_optoC = [];ripple_optoS = [];ripple_optoP = [];
    ripple_optoCC = [];ripple_optoSS = [];ripple_optoPP = [];
    rippleCC = {};rippleSS = {};ripplePP = {}
    for jj = 1:length(pulse)
        ripple_optoinC = find(ripples.peaks > (pulse(jj,1)-5) & 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)+5))
        ripple_optoP = [ripple_optoP;ripple_optoinP];
    
        ripple_optoCC(jj,1) = length(find(ripples.peaks > (pulse(jj,1)-5) & 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)+5)));
        
        rippleCC{jj,1} = ripples.peaks(find(ripples.peaks > (pulse(jj,1)-5) & 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)+5)))
        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);

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