clear all, close all
MedialSeptum_Recordings
% MS10: 61, 62, 63, 64 % gamma id: 62
% MS12: 78 (Peter_MS12_170714_122034_concat), 79 (Peter_MS12_170715_111545_concat), 80 (Peter_MS12_170716_172307_concat), 81 (Peter_MS12_170717_111614_concat), 83 (Peter_MS12_170719_095305_concat)
% MS13: 92 (Peter_MS13_171129_105507_concat), 93 (Peter_MS13_171130_121758_concat), 88 (Peter_MS13_171110_163224_concat)
% MS21: 126 (Peter_MS21_180629_110332_concat), 140 (Peter_MS21_180627_143449_concat), 143 (control: Peter_MS21_180719_155941_concat)
% MS22: 139 (Peter_MS22_180628_120341_concat), 127 (Peter_MS22_180629_110319_concat), 144 (control: Peter_MS22_180719_122813_concat)
id = 126; % 62
datasets = [61,62,63,64,78,79,80,81,83,92,93,88,126,140,143,139,127,144];
for k = 9:length(datasets)
    close all
    recording = recordings(datasets(k));
    disp(['Processing ', num2str(k),'/',num2str(length(datasets)),' recordingID = ' num2str(datasets(k)),', basename = ' recording.name])
    cd(fullfile(datapath, recording.animal_id, recording.name))
    load([recording.name,'.animal.behavior.mat'])
%     load([recording.name,'.cooling.manipulation.mat'])
    load([recording.name,'.trials.behavior.mat'])
    spikes = loadSpikes('clusteringpath',recording.SpikeSorting.path,'clusteringformat',recording.SpikeSorting.method,'basename',recording.name);
    [session, basename, basepath, clusteringpath] = db_set_path('session',recording.name);
    
    prebehaviortime = 0;
if recording.concat_behavior_nb > 0
    prebehaviortime = 0;
    nChannels = session.extracellular.nChannels;
    sr = session.extracellular.sr;
    if all(recording.concat_behavior_nb > 1)
        for i = 1:recording.concat_behavior_nb-1
            fullpath = fullfile([datapath,recording.animal_id], recording.concat_recordings{i}, [recording.concat_recordings{i}, '.dat']);
            temp2_ = dir(fullpath);
            prebehaviortime = prebehaviortime + temp2_.bytes/nChannels/2/sr;
        end
    end
    behaviortime = 0;
    for i = 1:length(recording.concat_behavior_nb)
        i1 = recording.concat_behavior_nb(i);
        fullpath = fullfile([datapath, recording.animal_id], recording.concat_recordings{i1}, [recording.concat_recordings{i1}, '.dat']);
        temp2_ = dir(fullpath);
        behaviortime = behaviortime+temp2_.bytes/nChannels/2/sr;
    end
    
else
    temp_ = dir(fname);
    behaviortime = temp_.bytes/nChannels/2/sr;
end
    
    
    % Firing rate linearized
    behavior = animal;
    behavior.state = trials.state;
    behavior.pos = animal.pos_linearized;
    behavior.rim = zeros(size(animal.arm));
    behavior.rim(find(animal.arm | animal.rim)) = 1;
    behavior.pos_limits = animal.pos_linearized_limits;
    behavior.maze.boundaries = [diff(animal.pos_y_limits),diff(animal.pos_y_limits)+abs(animal.maze.polar_theta_limits(1))-5];
    behavior.speed_th = 10;
    
    spikes2 = [];
    
    for i = 1:size(spikes.ts,2)
        spikes2.ts{i} = spikes.ts{i};
        spikes2.times{i} = spikes.times{i};
        spikes2.cluID(i) = spikes.cluID(i);
        spikes2.total(i) = spikes.total(i);
        spikes2.pos{i} = interp1(behavior.time,behavior.pos,spikes2.times{i});
        %     spikes2.pos_linearized{i} = interp1(behavior.time,behavior.pos_linearized,spikes2.times{i});
        spikes2.pos{i} = interp1(behavior.time,behavior.pos,spikes2.times{i});
        spikes2.speed{i} = interp1(behavior.time,behavior.speed,spikes2.times{i});
        spikes2.rim{i} = interp1(behavior.time,behavior.rim,spikes2.times{i},'nearest');
        spikes2.trials{i} = spikes.trials{i};
        spikes2.theta_freq{i} = spikes.theta_freq{i};
        spikes2.speed{i} = spikes.speed{i};
    end
    
    
    firingRateMap = plot_FiringRateMapAverage('animal',behavior,'spikes',spikes2);
    for i= 1:size(firingRateMap.map,2)
        firingRateMap.map{i}(isnan(firingRateMap.map{i}))= 0;
    end
    firingRateMap.total = spikes2.total./behaviortime;
    firingRateMap.boundaries = behavior.maze.boundaries;
    
    sucess = saveStruct(firingRateMap,'session',session,'datatype','firingRateMap');
    % save([recording.name, '.firingRateMap.firingRateMap.mat'],'firingRateMap')
    
    CoolingStates = plot_FiringRateMap2('animal',behavior,'spikes',spikes2,'trials',trials);
    CoolingStates.boundaries = behavior.maze.boundaries;
    CoolingStates.labels = behavior.state_labels;
    for i=1:length(behavior.state_labels)
        CoolingStates.trial_count(i) = length(unique(trials.trials2(find(trials.state==i))));
    end
    sucess = saveStruct(CoolingStates,'session',session,'datatype','firingRateMap');
    % save([recording.name, '.CoolingStates.firingRateMap.mat'],'CoolingStates')
    
    % Linearized
    % PhasePrecessionSlope_linearized = plot_FiringRateMap(behavior,spikes2,trials,theta,sr,'linearized_cooling_states');
    
    % Linearized and colored and grouped by left/right trials
    trials2 = trials;
    trials2.state(find(~isnan(trials2.state))) = trials.stat(trials.trials2(find(~isnan(trials2.state))));
    behavior2 = behavior;
    behavior2.state(find(~isnan(trials2.state))) = trials.stat(trials.trials2(find(~isnan(trials2.state))));
    behavior2.state_labels = trials.labels;
    % PhasePrecessionSlope_linearized_left_right = plot_FiringRateMap(behavior2,spikes2,trials2,theta,sr,'linearized_left_right');
    
    LeftRight = plot_FiringRateMap2('animal',behavior2,'spikes',spikes2,'trials',trials);
    LeftRight.boundaries = behavior.maze.boundaries;
    LeftRight.labels = {'Left','Right'};
    sucess = saveStruct(LeftRight,'session',session,'datatype','firingRateMap');
    % save([recording.name, '.LeftRight.firingRateMap.mat'],'LeftRight')
end