clear all;close all;

dataname = ['F:\oxytocin\ACh_OXT004\ACh_OXT004_012_230328_141439',   %% 1,5,10 seconds different amplitude stimulation
            'F:\oxytocin\ACh_OXT005\ACh_OXT005_010_230329_104411',   %% 1,5,10 seconds different amplitude stimulation
            'F:\oxytocin\ACh_OXT006\ACh_OXT006_006_230311_181030',   %% 1,2,5,10 seconds same amplitude stimulation
            'F:\oxytocin\ACh_OXT007\ACh_OXT007_009_230330_104326',   %% 1,5,10 seconds different amplitude stimulation
            'F:\oxytocin\ACh_OXT004\ACh_OXT004_021_230522_143652']   %% 1,5,20 seconds different amplitude stimulation

RipCh = [58,19,1,27,62];

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

    
    % sitmulation and signal
    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);
    delete digitalin
    % sitmulation and signal
    if exist([basename,'_FluoData.mat'])==0
        FluoData = Get_OxtAch_Fluos();
        figure
        plot(FluoData.time,FluoData.Ach_dFF); hold on;
        plot(FluoData.time,FluoData.Oxt_dFF);xlim([0 inf]);
    else load([basename,'_FluoData.mat']);plot(FluoData.Ach_dFF);
    end
    save([basename,'_FluoData.mat'],'FluoData');

%     ACh_Opto = []; Oxt_Opto = [];
%     pulse = pulses.intsPeriods-lfp_fbtime(1)/info.sample_rate;
%     pulses.duration = pulses.intsPeriods(:,2) - pulses.intsPeriods(:,1);
%     pulse = pulse(pulses.duration>4 & pulses.duration<6);
%     for jj = 1:length(pulse)
%         win = round(pulse(jj,1)*100) + [-10 20]*100;
%         if win(1>0) & win(2)<length(FluoData.Ach_dFF)
%             ACh_Opto(jj,:) = FluoData.Ach_dFF(win(1):win(2));
%             ACh_Opto(jj,995:1015) = nan;
%             ACh_Opto(jj,1995:2015) = nan;
%             ACh_Opto(jj,:) = fillmissing(ACh_Opto(jj,:),'linear');
% 
%             Oxt_Opto(jj,:) = FluoData.Oxt_dFF(win(1):win(2));
%             Oxt_Opto(jj,995:1015) = nan;
%             Oxt_Opto(jj,1995:2015) = nan;
%             Oxt_Opto(jj,:) = fillmissing(Oxt_Opto(jj,:),'linear');
%         end
%     end
%     figure;
%     options.color_area = [128 193 219]./255;    % Blue theme
%     options.color_line = [ 52 148 186]./255;
%     options.alpha      = 0.5;
%     options.line_width = 2;
%     options.error      = 'sem';
%     plot_areaerrorbar(ACh_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
%     xlim([1 3000]);
%     figure;
%     options.color_area = [128 193 219]./255;    % Blue theme
%     options.color_line = [ 52 148 186]./255;
%     options.alpha      = 0.5;
%     options.line_width = 2;
%     options.error      = 'sem';
%     plot_areaerrorbar(Oxt_Opto, options);set(gcf,'Position',  [100, 100, 400, 300])
%     xlim([1 3000]);

    %% Get PulseTable
    PulseData =  CreatePulseTable();
%% get the 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>200) = 0;
    load([basename,'.ripples.events.mat']);
%     ripples = bz_FindRipples(RipLFPs.data,RipLFPs.timestamps);
%% save to paper
    datpath = 'G:\My Drive\Paper_ACh_OXT\Fig6';
    ACh_OXT_Act = {};
    ACh_OXT_Act.PulseData =  PulseData;
    ACh_OXT_Act.pulse =  pulses.intsPeriods;
    ACh_OXT_Act.lfp_fbtime = lfp_fbtime;
    ACh_OXT_Act.Ach_dFF = FluoData.Ach_dFF;
    ACh_OXT_Act.Oxt_dFF = FluoData.Oxt_dFF;
    ACh_OXT_Act.LFP = RipLFPs;
    ACh_OXT_Act.ripples = ripples;
    save(fullfile(datpath,[basename,'ACh_OXT_Act.mat']),'ACh_OXT_Act');
end
