%% STATISTIC OXYTOCIN BEHAVIRO
dataname = ['F:\LAMP2\LAMP2_002_190405_154550']  % corrected
dataname = ['F:\LAMP3\LAMP3_006_190613_194519']  % corrected channel 53
dataname = [%'G:\data\LAMP37\LAMP37_005_201111_204459',
dataname = {'F:\LAMP37\LAMP37_006_201201_192319', % corrected channel 2
            'F:\LIMP014\LIMP014_012_211005_112625',    %%corrected   optitrack do not use  acclerate 2 times
            'F:\LIMP014\LIMP014_014_211007_122021',   %%corrected   optitrack  acclerate 2 times
            'F:\LIMP016\LIMP016_004_220101_144210',    %%corrected   analogCh  = 2  acclerate 2 times
            'F:\LIMP018\LIMP018_003_220127_150338',    %%corrected  exclude data from 5700 to 11250 acclerate 3 times
            'F:\LIMP019\LIMP019_004_220303_115500'}   %%corrected  acclerate 2 times acclerate 2 times

RipCh = [1,57,57,53, 26,19]
for ii = 1:size(dataname,1)
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.SleepState.states.mat']);
    load([basename,'.accel.behavior.mat'])

    if ii>1 & ii<4
        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,:)>0.6,1,'last');
        end 
    elseif ii>3
        info = read_Intan_RHD2000_file(basepath);
        analogin= read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
        if ii == 3
            DD = analogin(2,:);
        else ii>4
            DD = analogin(1,:);
        end
        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 ii>0 & ii<3
            load([name,ext],[name,'_OpTrTTL']);
            ACh_sti = eval([name,'_OpTrTTL']);
        elseif  ii>2
            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 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)];
    if ii ~=4
        xq = [1:(1/625):length(accel.acceleration.motion)];
    elseif ii == 4
        xq = [1:(3/1250):length(accel.acceleration.motion)];
    end
    Mov_time =  interp1(x,accel.acceleration.motion,xq);

    RipLFPs.data(Mov_time>200) = 0;

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

    
   ACh_ripples = {};
   ACh_ripples.ripples = ripples;
   ACh_ripples.dFoverF = dFoverF;
   ACh_ripples.alignTim = lfp_fbtime;
   ACh_ripples.lfpData = RipLFPs;
   datpath = 'G:\My Drive\Paper_ACh_OXT\Fig4\Figure 4A';
   save(fullfile(datpath,[basename,'_ACh_ripples.mat']),'ACh_ripples');
end
