clear all;close all;
%%%%%%%%%taking time of 1-3hr recording after CNO
dataname = ['F:\oxytocin\OXYT025\OXYT025_003_230101_131022',
            'F:\oxytocin\OXYT025\OXYT025_004_230105_135258',%%% 10mg/kg optitrack signal wrong no alignment
            'F:\oxytocin\OXYT025\OXYT025_012_230127_163348', %%% 10mg/kg optitrack signal wrong no alignment
            'F:\oxytocin\OXYT026\OXYT026_003_230105_094900',
            'F:\oxytocin\OXYT026\OXYT026_005_230113_165251',  %%% 10mg/kg optitrack signal name: AChSig_H 
            'F:\oxytocin\OXYT026\OXYT026_008_230117_172559',  %%% 10mg/kg optitrack signal name: AChSig_H 
            'F:\oxytocin\OXYT027\OXYT027_003_230104_141857']  %%% 5mg/kg optitrack signal name: AChSig_H 

RipCh = [23,23,58,27,8,27,11];
TheCh = [29,27,27,27,11,11,11];

for ii = 1:size(dataname,1);
    cd(dataname(ii,:))
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    load([basename,'.accel.behavior.mat']);
    

    %%%%%%%%%%%%%%%%% synchronize the fiber photometry signal
    % syn
    info = read_Intan_RHD2000_file(basepath);
     if ii ==2
        analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels)';
        lfp_fbtime = find(analogin(:,1)>(max(analogin(:,1))+min(analogin(:,1)))/2,1);
        delete analogin

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

    if ii ~=2
        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
    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\CNO';
    OXT_CNO = {};
    
    OXT_CNO.lfp_fbtime = lfp_fbtime;
    
    OXT_CNO.Oxt_dFF = dFoverF;
    OXT_CNO.RipLFP = RipLFPs;
    OXT_CNO.THLFP = THLFPs;
    OXT_CNO.ripples = ripples;
    OXT_CNO.spikes = spkMat;
    OXT_CNO.acel = accel.acceleration.motion;
    save(fullfile(datpath,[basename,'OXT_CNO.mat']),'OXT_CNO');
end                    
