function [infoTransf_maps_corr,infoTransf_maps_err,infoTransf_mapsSh_corr,infoTransf_mapsSh_err] = ...
    compute_DI_FIT_correrr_trials(neural,behav,sigTime,roiX,roiY,tmproiX,tmproiY,validXlab,validYlab,params,opts,band_label)
% correct and error DI and FIT comput (notice: FIT_S = FIT_C on corr/err trials)

% Notice: band_label is not really useful here, I just left it to simplify moving
% to across-frequencies computations

% Structure initialization
for featIdx = 1:numel(params.selected_features)
    featLab = params.selected_features{featIdx};
    % Correct trials structure
    infoTransf_maps_corr.(featLab).DI = [];
    infoTransf_maps_corr.(featLab).FIT_S = [];
    % Error trials structure
    infoTransf_maps_err.(featLab).DI = [];
    infoTransf_maps_err.(featLab).FIT_S = [];
    % Correct trials shuff structure
    infoTransf_mapsSh_corr.(featLab).DI = [];
    infoTransf_mapsSh_corr.(featLab).FIT_S = [];
    % Error trials shuff structure
    infoTransf_mapsSh_err.(featLab).DI = [];
    infoTransf_mapsSh_err.(featLab).FIT_S = [];
end


for sessIdx = 1:params.nSess
    %% Compute error trials, error first since we want to downsample number of correct trials to match error
    bX = []; bY = []; bhX = []; bhY = []; bXt = []; bYt = [];

    error_idxs = find(behav.v_correct{sessIdx}==0);
    % Behavioral features
    choice = behav.v_choice{sessIdx}(error_idxs);
    stim = choice; % stimulus and choice are equal on correct and error trials
    
    % Binning neural responses
    for t = 1:params.totTimePoints
        edgs = eqpop(neural.data_subj{sessIdx}(error_idxs,t,(2*roiX-1)), opts.n_bins);
        [~,bX.(band_label).left(:,t)] = histc(neural.data_subj{sessIdx}(error_idxs,t,(2*roiX-1)), edgs);
        edgs = eqpop(neural.data_subj{sessIdx}(error_idxs,t,(2*roiY-1)), opts.n_bins);
        [~,bY.(band_label).left(:,t)] = histc(neural.data_subj{sessIdx}(error_idxs,t,(2*roiY-1)), edgs);
        edgs = eqpop(neural.data_subj{sessIdx}(error_idxs,t,(2*roiX)), opts.n_bins);
        [~,bX.(band_label).right(:,t)] = histc(neural.data_subj{sessIdx}(error_idxs,t,(2*roiX)), edgs);
        edgs = eqpop(neural.data_subj{sessIdx}(error_idxs,t,(2*roiY)), opts.n_bins);
        [~,bY.(band_label).right(:,t)] = histc(neural.data_subj{sessIdx}(error_idxs,t,(2*roiY)), edgs);
    end

    % Compute DI and FIT error
    emitFreq = band_label;
    recFreq = band_label;
        
    for featIdx = 1:numel(params.selected_features)
        featLab = params.selected_features{featIdx};

        % Time is w.r.t. receiver Y
        for tidx = 1:numel(sigTime.(recFreq).(validYlab))
            t = sigTime.(recFreq).(validYlab)(tidx);
            % For each time point we select a range of delays so that
            % only the possible infomration-emitting time steps are
            % considered, this increases the computational efficiency
            validDelays = t-sigTime.(emitFreq).(validXlab); 
            validDelays(validDelays<1)=[]; % remove negative delays (emitter in the future)
            validDelays(validDelays>params.maxDelay)=[]; % remove delays longer than params.max_delay
            validDelays = sort(validDelays);

            if ~isempty(validDelays) % check that some valid delays exist
                for d = validDelays
                    % Define hX (h = history --> past), hY, Xt and Yt

                    % Define optim. features
                    % Y present
                    bYt.(recFreq).(featLab) = bY.(recFreq).(featLab)(:,t)';

                    % X past
                    bhX.(emitFreq).(featLab) = bX.(emitFreq).(featLab)(:,t-d)';

                    % Y past
                    bhY.(recFreq).(featLab) = bY.(recFreq).(featLab)(:,t-d)';
    
                    % Compute DI and FIT for error trials (only stimulus
                    % FIT, since stimulus = choice on error trials)
                    [tmpDi, ~, tmpFITs] = compute_FIT_TE(stim, bhX.(emitFreq).(featLab), bYt.(recFreq).(featLab), bhY.(recFreq).(featLab));
                    
                    % DI X to Y
                    infoTransf_maps_err.(featLab).DI(sessIdx,tidx,d) = tmpDi;

                    % FIT_S X to Y 
                    infoTransf_maps_err.(featLab).FIT_S(sessIdx,tidx,d) = tmpFITs;
                    
                end % end delay loop
            end 
        end % end time loop
    end % end features loop
    
    %% Compute correct trials    
    bX = []; bY = []; bhX = []; bhY = []; bXt = []; bYt = [];

    correct_idxs = find(behav.v_correct{sessIdx}==1);
    nCorrTrials = numel(correct_idxs);
    nDiff = nCorrTrials - numel(error_idxs);
    if strcmp(params.corrSubsamp, 'mid')
        correct_idxs(1:ceil(nDiff/2)) = [];
        correct_idxs(end-floor(nDiff/2)+1:end) = [];
    elseif strcmp(params.corrSubsamp, 'rand')
        rng(params.randSeed)
        rIdx = randperm(nCorrTrials);
        rIdx = rIdx(1:numel(error_idxs));
        correct_idxs = correct_idxs(rIdx); 
    end
    % Behavioral features
    choice = behav.v_choice{sessIdx}(correct_idxs);
    stim = choice; % stimulus and choice are equal on correct and error trials

    % Binning neural responses
    for t = 1:params.totTimePoints
        edgs = eqpop(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiX-1)), opts.n_bins);
        [~,bX.(band_label).left(:,t)] = histc(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiX-1)), edgs);
        edgs = eqpop(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiY-1)), opts.n_bins);
        [~,bY.(band_label).left(:,t)] = histc(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiY-1)), edgs);
        edgs = eqpop(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiX)), opts.n_bins);
        [~,bX.(band_label).right(:,t)] = histc(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiX)), edgs);
        edgs = eqpop(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiY)), opts.n_bins);
        [~,bY.(band_label).right(:,t)] = histc(neural.data_subj{sessIdx}(correct_idxs,t,(2*roiY)), edgs);
    end

    % Compute DI and FIT correct
    emitFreq = band_label;
    recFreq = band_label;
    
    for featIdx = 1:numel(params.selected_features)
        featLab = params.selected_features{featIdx};
        % Time is w.r.t. receiver Y
        for tidx = 1:numel(sigTime.(recFreq).(validYlab))
            t = sigTime.(recFreq).(validYlab)(tidx);
            % For each time point we select a range of delays so that
            % only the possible infomration-emitting time steps are
            % considered, this increases the computational efficiency
            validDelays = t-sigTime.(emitFreq).(validXlab); 
            validDelays(validDelays<1)=[]; % remove negative delays (emitter in the future)
            validDelays(validDelays>params.maxDelay)=[]; % remove delays longer than params.max_delay
            validDelays = sort(validDelays);
    
            if ~isempty(validDelays) % check that some valid delays exist
                for d = validDelays
                    % Define hX (h = history --> past), hY, Xt and Yt
    
                    % Y present
                    bYt.(recFreq).(featLab) = bY.(recFreq).(featLab)(:,t)';

                    % X past
                    bhX.(emitFreq).(featLab) = bX.(emitFreq).(featLab)(:,t-d)';

                    % Y past
                    bhY.(recFreq).(featLab) = bY.(recFreq).(featLab)(:,t-d)';
                    
                    % Compute DI and FIT for correct trials (only stimulus
                    % FIT, since stimulus = choice on error trials)
                    [tmpDi, ~, tmpFITs] = compute_FIT_TE(stim, bhX.(emitFreq).(featLab), bYt.(recFreq).(featLab), bhY.(recFreq).(featLab));
                    
                    % DI X to Y
                    infoTransf_maps_corr.(featLab).DI(sessIdx,tidx,d) = tmpDi;

                    % FIT_S X to Y 
                    infoTransf_maps_corr.(featLab).FIT_S(sessIdx,tidx,d) = tmpFITs;
                    
                end % end delay loop
            end
        end % end time loop
    end % end features loop
end % end sess loop

end

