clear all;close all;
%%%%%%% Lateral septum mediated effect comparison
dataname = [ 'F:\oxytocin\OXYT025\OXYT025_001_221231_123130', %% optitrack,  shank 2,3,4 is in hippocampus
             'F:\oxytocin\OXYT026\OXYT026_001_230103_151806', %% optitrack, shank 3,4 is in hippocampus
             'F:\oxytocin\OXYT027\OXYT027_001_230103_201733', %% optitrack, only shank 4 is in hippocampus
             'F:\oxytocin\OXYT027\OXYT027_002_230104_110308',
             'F:\oxytocin\OXYT025\OXYT025_005_230110_110227', %% only 4000s
             'F:\oxytocin\OXYT025\OXYT025_006_230116_123923'] %% only 1-6000s
RipCh = [27,11,11,23,23];
TheCh = [11,11,11,11,11];

for ii = 1:size(dataname,1);
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.accel.behavior.mat']);
    if exist([basename,'.pulses.events.mat'])>0
        load([basename,'.pulses.events.mat']);
    end 
    

    %%%%%%%%%%%%%%%%% synchronize the fiber photometry signal
    % syn
    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);

    name = basename(1:11);
    if exist('dFoverF.mat')==0
        ext = '.mat'
        load([name,ext],[name,'_OptoTTL']);
        Opti = eval([name,'_OptoTTL']);
        Opti.data = downsample(Opti.values,10);
        Opti.time = downsample(Opti.times,10);
        time_st = round(Opti.time(find(Opti.data>1,1))*2e4);time_nd = round(Opti.time(find(Opti.data>1,1,'last'))*2e4);
        % signal
        load([name,ext],[name,'_AChSig_H']);
        ACh_signal = eval([name,'_AChSig_H']);
        ACh_signal.data = downsample(ACh_signal.values(time_st:time_nd),200);
        ACh_signal.time = downsample(ACh_signal.times(time_st:time_nd),200)- Opti.time(find(Opti.data>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]);
        title([name]);
    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');

    THLFPs = bz_GetLFP(TheCh(ii),'basepath',basepath,'noPrompts',true);
    THLFPs.data = double(THLFPs.data);
    THLFPs.data(abs(THLFPs.data)>5000) = nan;
    THLFPs.data = fillmissing(THLFPs.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 spikes
    load([basename,'.spikes.cellinfo.mat']);
    binSize = 0.01;
    spkMat = bz_SpktToSpkmat(spikes,'binSize',binSize,'overlap',1);
    %% arrange structure
    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig8\Control';
    OXT_CNO_Ctrl = {};
    if exist([basename,'.pulses.events.mat'])>0
        OXT_CNO_Ctrl.pulse = pulses.intsPeriods;
    end 
    
    OXT_CNO_Ctrl.lfp_fbtime = lfp_fbtime;
    OXT_CNO_Ctrl.Oxt_dFF = dFoverF;
    OXT_CNO_Ctrl.RipLFP = RipLFPs;
    OXT_CNO_Ctrl.THLFP = THLFPs;
    OXT_CNO_Ctrl.ripples = ripples;
    OXT_CNO_Ctrl.spikes = spkMat;
    OXT_CNO_Ctrl.acel = accel.acceleration.motion;
    save(fullfile(datpath,[basename,'OXT_CNO_Ctrl.mat']),'OXT_CNO_Ctrl');
end