function [widthCell, varargout] = plotParamsThresholdedDistribution(StimTableCell, threshold, qualityMeasure, figDir, varargin) 
% PLOTPARAMSTHRESHOLDEDDISTRIBUTION  Plots Gaussian fit params distributions for
% a given threshold of a quality measure.
%   C = PLOTSCATTERPARAMSVSQUALITYINDEX(STIMTABLE, THRESHOLDARRAY, QUALITYMEASURE, FIGDIR) 
%   plots the STIMTABLE params for the cutoffs of THRESHOLDARRAY for the 
%   selected QUALITYMEASURE, and saves to FIGDIR.
%
%   C = ADDME(A,B) adds A and B together.
%
%   See also SUM, PLUS.

% Maybe still align the yaxis labels and include the correlation coefficent
% if you want to make some claims.
switch nargin
    case 4        
        rSquareThreshold = -20;
        alignToOff = false;
        alignToAll = false;
    case 5
        rSquareThreshold = varargin{1};        
        alignToOff = false;
        alignToAll = false;
    case 6
        rSquareThreshold = varargin{1};        
        alignToOff = varargin{2};
        alignToAll = false;
    case 7
        rSquareThreshold = varargin{1};        
        alignToOff = varargin{2};
        alignToAll = varargin{3};
    otherwise
        warning('Too many arguments, cannot process them meaningfully!');
end

[~, ...
 ~, ...
 fitParams, ...
 flyInds, ...
 isAboveThreshold, ...
 tcExt, ...
 tcArgExt] = getAlignedTuningCurveArrayFromTable(StimTableCell, ...
                                                         threshold, ...
                                                         rSquareThreshold, ...
                                                         alignToOff, ...
                                                         alignToAll);

nStims = numel(StimTableCell);
for iStim = 1: nStims
    StimTable = StimTableCell{iStim};
    [~, ~, ~, responseIndArray, ~] = getParamsArrayFromStimTable(StimTable, qualityMeasure);
    if threshold >= max(responseIndArray)
        warning(['Went over threshold respQuality ' num2str(threshold) ...
                 '. Range is [' num2str(min(responseIndArray)) ...
                 ', ' num2str(max(responseIndArray)) ']']);
%         threshold = min(responseIndArray);
    end
    
    [~, gof] = getFitParamsStructArrayFromTable(StimTable);
    rSquare = [gof.rsquare];
    if rSquareThreshold >= max(rSquare)
        warning(['Went over threshold rSquare ' num2str(rSquareThreshold) ...
                 '. Range is [' num2str(min(rSquare)) ', ' num2str(max(rSquare)) ']']);
%         rSquareThreshold = min(rSquare);
    end
    
    stimDegreeFactor = 2;
    amplitudeCell{iStim} = fitParams.a1{iStim};
    positionCell{iStim} = stimDegreeFactor * (fitParams.b1{iStim} - 1);
    widthCell{iStim} = stimDegreeFactor  * 2 * sqrt(log(2)) * fitParams.c1{iStim}; % full-width half maximum conversion    amplitudeCell{iStim} = a1;
    nRois(iStim) = sum(isAboveThreshold{iStim});
    nTotalRois(iStim) = numel(isAboveThreshold{iStim}); 
    % Check how many flies have good cells. This neeeds to include rSquare.
    nContributingFlies(iStim) = numel(unique(flyInds{iStim}));
    tcArgExt{iStim} = (tcArgExt{iStim} - 1) * stimDegreeFactor;
end

metadata.nRois = nRois;
metadata.nTotalRois = nTotalRois;
metadata.nContributingFlies = nContributingFlies;

varargout{1} = metadata;
legendStr = arrayfun(@(x, y, z) ['N = ' num2str(x) ...
                                ' (n = ' num2str(y) '/' num2str(z) ')'], ...
                     nContributingFlies, nRois, nTotalRois, 'uni', 0);

fileStimSuffix = {'ONOFF', 'XOFF', 'YOFF', 'XON', 'YON'};
%%
dataCell = [amplitudeCell; positionCell; widthCell; tcExt; tcArgExt];
yLabelCell = {'Amplitude \DeltaF/F_0', 'Center position (\circ)', ...
              'FWHM (\circ)', 'Extreme \DeltaF/F_0', ...
              'Extreme position (\circ)'};
fileNameSuffixCell = {'Amp', 'Pos', 'Width', 'NonParamAmp', 'NonParamPos'};
supportCell = {'unbounded', 'positive', 'positive', 'unbounded', 'unbounded'};
supportCell = {'unbounded', [0-eps 60+eps], [0-eps 100], 'unbounded', [0-eps 60+eps]};
if rSquareThreshold < 0
    supportCell{2} = 'unbounded';
end
nParams = size(dataCell, 1);
for iParam = 1: nParams    
    data = dataCell(iParam, :);
    hFig = createPrintFig(10 * [1 3 / 4]);
    hParamsSubAx(1) = axes('Parent', hFig);
% Set to pi in case the array is empty.
isEmptyArray = cellfun(@isempty, data);
if any(isEmptyArray)
    data(isEmptyArray) = {pi};
end
hBeanAx = beanPlot(data, fileStimSuffix(2: end), getONOFFXYColors, 0.5, ...
                   hParamsSubAx(1), supportCell{iParam}, 'vertical', 0, 0, 1);

% axis(hParamsSubAx, 'square');
hYLabel(1) = ylabel(hParamsSubAx(1), yLabelCell{iParam});

hLegend = legend(findall(hBeanAx, 'Type', 'patch'), legendStr, ...
                 'Location', 'northoutside');
setFontForThesis([hParamsSubAx, hBeanAx, hLegend], hFig)
arrayfun(@prettifyAxes, hParamsSubAx); 
arrayfun(@offsetAxes, hParamsSubAx); 
[~, ~, genotype] = parseFlyPath(fileparts(StimTable.timeSeriesPath{1}));

figFileName = ['ReceptiveField' fileNameSuffixCell{iParam} '-' genotype '-' ...
               qualityMeasure '-' num2str(threshold) ...
               '-rSquare-' num2str(rSquareThreshold)...
               '-All' num2str(alignToAll) ...
               '-Off' num2str(alignToOff)];
% set(gcf,'renderer','painters')
% set(gcf, 'PaperPositionMode', 'auto');
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
print(hFig, [figDir figFileName '.png'], '-dpng', '-r300')
% print(hFig, [figDir figFileName '.eps'], '-depsc2', '-r300')
end

end