function printDistributionStatistics(pbeDefinition)

    disp('--- DISTRIBUTION STATISTICS ---')
    % precalculations
    positionMatrix = pbeDefinition.seeding.disD.positionMatrix;
    numberVector = 1/pbeDefinition.seeding.nSample * ...
        pbeDefinition.seeding.fhandle(positionMatrix) ./ ...
        pbeDefinition.seeding.disD.probabilityVector;
    % have we obtained the right distribution
    populationCenter = mean(positionMatrix, 2)'
    populationStd = std(positionMatrix, [], 2)'
    maxDistanceFromCenter = max(sqrt(sum(...
          bsxfun(@minus, positionMatrix, pbeDefinition.seeding.mean).^2, 1)))
    
    disp('- fractions in different regions')
    marginFraction = ...
        sum(pbeDefinition.seeding.marginFilter(positionMatrix)) / ...
        length(numberVector)
    validFraction = ...
        sum(isValidDec(positionMatrix, pbeDefinition.dec)) / ...
        length(numberVector)
    totalFraction = sum(isValidDec(positionMatrix, pbeDefinition.dec) & ...
        pbeDefinition.seeding.marginFilter(positionMatrix)) / ...
        length(numberVector)

    cutoff10Fraction = sum(...
        sqrt(sum(...
          bsxfun(@minus, positionMatrix, pbeDefinition.seeding.mean).^2, 1)...
        ) <= 1.0*pbeDefinition.seeding.std(1)) / length(numberVector)
    cutoff20Fraction = sum(...
        sqrt(sum(...
          bsxfun(@minus, positionMatrix, pbeDefinition.seeding.mean).^2, 1)...
        ) <= 2.0*pbeDefinition.seeding.std(1)) / length(numberVector)
    cutoff30Fraction = sum(...
        sqrt(sum(...
        bsxfun(@minus, positionMatrix, pbeDefinition.seeding.mean).^2, 1)...
        ) <= 3.0*pbeDefinition.seeding.std(1)) / length(numberVector)
        
    disp('- integral properties')
    totalNumber = sum(numberVector)
%     totalNumberWithMargin = ...
%         pbeDefinition.seeding.validityFilter(positionMatrix) * numberVector'
%     numberFraction = sum(computeMeasureDec(pbeDefinition.dec, 'hC', ...
%         positionMatrix, 'volume', [], [], -1) > 0) / length(numberVector)
    totalVolume = (pbeDefinition.seeding.marginFilter(positionMatrix) .* ...
        computeMeasureDec(pbeDefinition.dec, 'hC', ...
          positionMatrix, 'volume', [], [], -1)) * numberVector'
    totalSurface = computeMeasureDec(pbeDefinition.dec, 'hC', ...
        positionMatrix, 'surface', [], [], -1) * numberVector'
    totalMeanWidth = computeMeasureDec(pbeDefinition.dec, 'hC', ...
        positionMatrix, 'meanWidth', [], [], -1) * numberVector'
%     totalMeanWidth2 = 0;
%     for iSample = 1:pbeDefinition.seeding.nSample
%         totalMeanWidth2 = totalMeanWidth2 + ...
%             1/pbeDefinition.seeding.nSample * ...
%         pbeDefinition.seeding.meanWidthHandle(positionMatrix(:, iSample)) ./ ...
%         pbeDefinition.seeding.disD.probabilityVector(iSample);
%     end
%     totalMeanWidth2
    
    totalMoment3 = sum(positionMatrix.^3, 1) * numberVector'

    disp('- plotting 1D distributions')
%     % particle number over mean width
%     binLower = 0; binUpper = 1e-17; binSpacing = 1e-18;
% %     binLower = 0; binUpper = 4e-6; binSpacing = 1e-7;
%     % bounds in [lower bounds; upper bounds]
%     boundDefinitionMatrix = [binLower:binSpacing:(binUpper - binSpacing);...
%         (binLower + binSpacing):binSpacing:binUpper]
%     % pivot position
%     pivotPositionVector = mean(boundDefinitionMatrix, 1)
%     nPivot = length(pivotPositionVector);
%     
%     % handles for integration
%     disp('..measure y')
% 	yVector = numberVector;
%     size(yVector)
%     disp('..measure x')
%     xVector = computeMeasureDec(pbeDefinition.dec, 'hC', ...
%         positionMatrix, 'volume', [], [], -1);
%     size(xVector)
%     disp('..measure finished')
%     
%     %   integrand looks like:
%     %     yProperty*(xProperty > lowerBound)*(xProperty <= upperBound)
%     %   Now consider that yProperty is the vector as given by
%     %   yAxisHandle. While yProperty is a row vector, the boundary
%     %   filters can be written as matrices that are multiplied from the
%     %   right:
%     %     yProperty*boundaryFilterMatrix
%     %   so that the result of the integral will become a column vector.
%     % Hence, we construct the boundary filter, first.
%     boundaryFilter = ...
%             bsxfun(@gt, xVector, boundDefinitionMatrix(1,:)') .* ...
%             bsxfun(@le, xVector, boundDefinitionMatrix(2,:)');
%     size(boundaryFilter)
%     % integrand
%     integrandVector = boundaryFilter .* ...
%                      (ones(size(boundaryFilter, 1), 1)*yVector) * ...
%                      numberVector';
%     size(integrandVector)
%     figure, plot(pivotPositionVector, integrandVector, 'bx-')    
% %     title('particle number distr. density over mean width')
%     ylabel('number distr. density in #/m^3')
%     xlabel('mean width in m^3')
%     % % set(gca,'xscale', 'log')

end