
clear all;

datafile = {'F:\oxytocin\ACh_OXT004\ACh_OXT004_015_230411_150537',
            'F:\oxytocin\ACh_OXT004\ACh_OXT004_016_230412_145104',
            'F:\oxytocin\ACh_OXT004\ACh_OXT004_017_230412_151435'}
Ch = [53,53,53]

for ii = 1:length(datafile)
    cd(datafile{ii,:})
    
    basepath = pwd;
    basename = bz_BasenameFromBasepath(basepath);
    

    % get alighment
    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);
    % get name 
    load([basename,'_Fluo.mat']);
    %% get LFPs
    THLFPs = bz_GetLFP(Ch(ii),'basepath',basepath,'noPrompts',true);
    thetarange = linspace(0,250,1000);
    lfp = double(THLFPs.data);
    lfp(abs(lfp)>5000) = nan;
    lfp = fillmissing(lfp,'linear');
    [ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
    [spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',1);
    bands = SpectrogramBands(spectrogrammt,fm);



    figure;
    colormap('jet')
    imagesc(tt,thetarange,abs(ss));clim([0,2e5]);ylim([0 50])
    yticks([0.5,20,50,100,250]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
    set(gca,'YDir','normal');
    hold on;plot(FluoData.time+lfp_fbtime(1)/20000, (FluoData.Ach_dFF)*200+30,'w-');
    hold on;plot(tm, bands.theta/8000+20,'g-');
    %%

%     datpath = 'G:\My Drive\Paper_ACh_OXT\Fig1\Data\Figure1F\OXT_theta'
%     save(fullfile(datpath,[basename,'_LFPs.mat']),'THLFPs','lfp_fbtime');
end
