diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index b3aaac0dbf..8dbf0ec5a8 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -27,8 +27,8 @@ License
 #include "volFields.H"
 #include "regionSplit.H"
 #include "fvcVolumeIntegrate.H"
-#include "Histogram.H"
 #include "mathematicalConstants.H"
+#include "stringListOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -56,6 +56,281 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::regionSizeDistribution::writeGraph
+(
+    const coordSet& coords,
+    const word& valueName,
+    const scalarField& values
+) const
+{
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    const wordList valNames(1, valueName);
+
+    fileName outputPath;
+    if (Pstream::parRun())
+    {
+        outputPath = mesh.time().path()/".."/name_;
+    }
+    else
+    {
+        outputPath = mesh.time().path()/name_;
+    }
+
+    if (mesh.name() != fvMesh::defaultRegion)
+    {
+        outputPath = outputPath/mesh.name();
+    }
+
+    mkDir(outputPath/mesh.time().timeName());
+    OFstream str
+    (
+        outputPath
+      / mesh.time().timeName()
+      / formatterPtr_().getFileName(coords, valNames)
+    );
+    Info<< "Writing distribution of " << valueName << " to " << str.name()
+        << endl;
+
+    List<const scalarField*> valPtrs(1);
+    valPtrs[0] = &values;
+    formatterPtr_().write(coords, valNames, valPtrs, str);
+}
+
+
+void Foam::regionSizeDistribution::writeAlphaFields
+(
+    const regionSplit& regions,
+    const Map<label>& patchRegions,
+    const Map<scalar>& regionVolume,
+    const volScalarField& alpha
+) const
+{
+    const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
+
+    // Split alpha field
+    // ~~~~~~~~~~~~~~~~~
+    // Split into
+    //  - liquidCore            : region connected to inlet patches
+    //  - per region a volume   : for all other regions
+    //  - backgroundAlpha       : remaining alpha
+
+
+    // Construct field
+    volScalarField liquidCore
+    (
+        IOobject
+        (
+            alphaName_ + "_liquidCore",
+            obr_.time().timeName(),
+            obr_,
+            IOobject::NO_READ
+        ),
+        alpha,
+        fvPatchField<scalar>::calculatedType()
+    );
+
+    volScalarField backgroundAlpha
+    (
+        IOobject
+        (
+            alphaName_ + "_background",
+            obr_.time().timeName(),
+            obr_,
+            IOobject::NO_READ
+        ),
+        alpha,
+        fvPatchField<scalar>::calculatedType()
+    );
+
+
+    // Knock out any cell not in patchRegions
+    forAll(liquidCore, cellI)
+    {
+        label regionI = regions[cellI];
+        if (patchRegions.found(regionI))
+        {
+            backgroundAlpha[cellI] = 0;
+        }
+        else
+        {
+            liquidCore[cellI] = 0;
+
+            scalar regionVol = regionVolume[regionI];
+            if (regionVol < maxDropletVol)
+            {
+                backgroundAlpha[cellI] = 0;
+            }
+        }
+    }
+    liquidCore.correctBoundaryConditions();
+    backgroundAlpha.correctBoundaryConditions();
+
+    Info<< "Volume of liquid-core = "
+        << fvc::domainIntegrate(liquidCore).value()
+        << endl;
+    Info<< "Volume of background  = "
+        << fvc::domainIntegrate(backgroundAlpha).value()
+        << endl;
+
+    Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
+    liquidCore.write();
+    Info<< "Writing background field to " << backgroundAlpha.name() << endl;
+    backgroundAlpha.write();
+}
+
+
+Foam::Map<Foam::label> Foam::regionSizeDistribution::findPatchRegions
+(
+    const polyMesh& mesh,
+    const regionSplit& regions
+) const
+{
+    // Mark all regions starting at patches
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Count number of patch faces (just for initial sizing)
+    const labelHashSet patchIDs(mesh.boundaryMesh().patchSet(patchNames_));
+
+    label nPatchFaces = 0;
+    forAllConstIter(labelHashSet, patchIDs, iter)
+    {
+        nPatchFaces += mesh.boundaryMesh()[iter.key()].size();
+    }
+
+
+    Map<label> patchRegions(nPatchFaces);
+    forAllConstIter(labelHashSet, patchIDs, iter)
+    {
+        const polyPatch& pp = mesh.boundaryMesh()[iter.key()];
+
+        // Collect all regions on the patch
+        const labelList& faceCells = pp.faceCells();
+
+        forAll(faceCells, i)
+        {
+            patchRegions.insert
+            (
+                regions[faceCells[i]],
+                Pstream::myProcNo()     // dummy value
+            );
+        }
+    }
+
+
+    // Make sure all the processors have the same set of regions
+    Pstream::mapCombineGather(patchRegions, minEqOp<label>());
+    Pstream::mapCombineScatter(patchRegions);
+
+    return patchRegions;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::regionSizeDistribution::divide
+(
+    const scalarField& num,
+    const scalarField& denom
+)
+{
+    tmp<scalarField> tresult(new scalarField(num.size()));
+    scalarField& result = tresult();
+
+    forAll(denom, i)
+    {
+        if (denom[i] != 0)
+        {
+            result[i] = num[i]/denom[i];
+        }
+        else
+        {
+            result[i] = 0.0;
+        }
+    }
+    return tresult;
+}
+
+
+void Foam::regionSizeDistribution::writeGraphs
+(
+    const word& fieldName,              // name of field
+    const labelList& indices,           // index of bin for each region
+    const scalarField& sortedField,     // per region field data
+    const scalarField& binCount,        // per bin number of regions
+    const coordSet& coords              // graph data for bins
+) const
+{
+    if (Pstream::master())
+    {
+        // Calculate per-bin average
+        scalarField binSum(nBins_, 0.0);
+        forAll(sortedField, i)
+        {
+            binSum[indices[i]] += sortedField[i];
+        }
+
+        scalarField binAvg(divide(binSum, binCount));
+
+        // Per bin deviation
+        scalarField binSqrSum(nBins_, 0.0);
+        forAll(sortedField, i)
+        {
+            binSqrSum[indices[i]] += Foam::sqr(sortedField[i]);
+        }
+        scalarField binDev
+        (
+            sqrt(divide(binSqrSum, binCount) - Foam::sqr(binAvg))
+        );
+
+        // Write average
+        writeGraph(coords, fieldName + "_sum", binSum);
+        // Write average
+        writeGraph(coords, fieldName + "_avg", binAvg);
+        // Write deviation
+        writeGraph(coords, fieldName + "_dev", binDev);
+    }
+}
+
+
+void Foam::regionSizeDistribution::writeGraphs
+(
+    const word& fieldName,              // name of field
+    const scalarField& cellField,       // per cell field data
+    const regionSplit& regions,         // per cell the region(=droplet)
+    const labelList& sortedRegions,     // valid regions in sorted order
+    const scalarField& sortedNormalisation,
+
+    const labelList& indices,           // per region index of bin
+    const scalarField& binCount,        // per bin number of regions
+    const coordSet& coords              // graph data for bins
+) const
+{
+    // Sum on a per-region basis. Parallel reduced.
+    Map<scalar> regionField(regionSum(regions, cellField));
+
+    // Extract in region order
+    scalarField sortedField
+    (
+        sortedNormalisation
+      * extractData
+        (
+            sortedRegions,
+            regionField
+        )
+    );
+
+    writeGraphs
+    (
+        fieldName,      // name of field
+        indices,        // index of bin for each region
+        sortedField,    // per region field data
+        binCount,       // per bin number of regions
+        coords          // graph data for bins
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::regionSizeDistribution::regionSizeDistribution
@@ -103,8 +378,11 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
         dict.lookup("field") >> alphaName_;
         dict.lookup("patches") >> patchNames_;
         dict.lookup("threshold") >> threshold_;
-        dict.lookup("volFraction") >> volFraction_;
+        dict.lookup("maxDiameter") >> maxDiam_;
+        minDiam_ = 0.0;
+        dict.readIfPresent("minDiameter", minDiam_);
         dict.lookup("nBins") >> nBins_;
+        dict.lookup("fields") >> fields_;
 
         word format(dict.lookup("setFormat"));
         formatterPtr_ = writer<scalar>::New(format);
@@ -163,14 +441,17 @@ void Foam::regionSizeDistribution::write()
            : obr_.lookupObject<volScalarField>(alphaName_)
         );
 
-        Info<< "Volume of alpha = "
+        Info<< "Volume of alpha          = "
             << fvc::domainIntegrate(alpha).value()
             << endl;
 
         const scalar meshVol = gSum(mesh.V());
-        Info<< "Mesh volume = " << meshVol << endl;
-        Info<< "Background region volume limit = " << volFraction_*meshVol
-            << endl;
+        const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
+        const scalar delta = (maxDiam_-minDiam_)/nBins_;
+
+        Info<< "Mesh volume              = " << meshVol << endl;
+        Info<< "Maximum droplet diameter = " << maxDiam_ << endl;
+        Info<< "Maximum droplet volume   = " << maxDropletVol << endl;
 
 
         // Determine blocked faces
@@ -260,138 +541,115 @@ void Foam::regionSizeDistribution::write()
         }
 
 
+        // Determine regions connected to supplied patches
+        Map<label> patchRegions(findPatchRegions(mesh, regions));
+
+
+
         // Sum all regions
-        Map<Pair<scalar> > regionVolume(regions.nRegions()/Pstream::nProcs());
-        forAll(alpha, cellI)
-        {
-            scalar cellVol = mesh.V()[cellI];
-            scalar alphaVol = alpha[cellI]*cellVol;
-
-            label regionI = regions[cellI];
-
-            Map<Pair<scalar> >::iterator fnd = regionVolume.find(regionI);
-            if (fnd == regionVolume.end())
-            {
-                regionVolume.insert
-                (
-                    regionI,
-                    Pair<scalar>(cellVol, alphaVol)
-                );
-            }
-            else
-            {
-                fnd().first() += cellVol;
-                fnd().second() += alphaVol;
-            }
-        }
-        Pstream::mapCombineGather(regionVolume, ListPlusEqOp<scalar, 2>());
-        Pstream::mapCombineScatter(regionVolume);
-
+        const scalarField alphaVol = alpha.internalField()*mesh.V();
+        Map<scalar> allRegionVolume(regionSum(regions, mesh.V()));
+        Map<scalar> allRegionAlphaVolume(regionSum(regions, alphaVol));
+        Map<label> allRegionNumCells
+        (
+            regionSum
+            (
+                regions,
+                labelField(mesh.nCells(), 1.0)
+            )
+        );
 
         if (debug)
         {
             Info<< token::TAB << "Region"
                 << token::TAB << "Volume(mesh)"
                 << token::TAB << "Volume(" << alpha.name() << "):"
+                << token::TAB << "nCells"
                 << endl;
             scalar meshSumVol = 0.0;
             scalar alphaSumVol = 0.0;
+            label nCells = 0;
 
-            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            Map<scalar>::const_iterator vIter = allRegionVolume.begin();
+            Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
+            Map<label>::const_iterator numIter = allRegionNumCells.begin();
+            for
+            (
+                ;
+                vIter != allRegionVolume.end()
+             && aIter != allRegionAlphaVolume.end();
+                ++vIter, ++aIter, ++numIter
+            )
             {
-                Info<< token::TAB << iter.key()
-                    << token::TAB << iter().first()
-                    << token::TAB << iter().second() << endl;
+                Info<< token::TAB << vIter.key()
+                    << token::TAB << vIter()
+                    << token::TAB << aIter()
+                    << token::TAB << numIter()
+                    << endl;
 
-                meshSumVol += iter().first();
-                alphaSumVol += iter().second();
+                meshSumVol += vIter();
+                alphaSumVol += aIter();
+                nCells += numIter();
             }
             Info<< token::TAB << "Total:"
                 << token::TAB << meshSumVol
-                << token::TAB << alphaSumVol << endl;
+                << token::TAB << alphaSumVol
+                << token::TAB << nCells
+                << endl;
             Info<< endl;
         }
 
 
 
-        // Mark all regions starting at patches
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        // Count number of patch faces (just for initial sizing)
-        label nPatchFaces = 0;
-        forAll(patchNames_, i)
         {
-            const word& pName = patchNames_[i];
-            label patchI = mesh.boundaryMesh().findPatchID(pName);
-            if (patchI == -1)
+            Info<< "Patch connected regions (liquid core):" << endl;
+            Info<< token::TAB << "Region"
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            forAllConstIter(Map<label>, patchRegions, iter)
             {
-                WarningIn("regionSizeDistribution::write()")
-                    << "Cannot find patch " << pName << ". Valid patches are "
-                    << mesh.boundaryMesh().names()
-                    << endl;
-            }
-            else
-            {
-                nPatchFaces += mesh.boundaryMesh()[patchI].size();
+                label regionI = iter.key();
+                Info<< token::TAB << iter.key()
+                    << token::TAB << allRegionVolume[regionI]
+                    << token::TAB << allRegionAlphaVolume[regionI] << endl;
+
             }
+            Info<< endl;
         }
 
-        Map<label> keepRegions(nPatchFaces);
-        forAll(patchNames_, i)
         {
-            const word& pName = patchNames_[i];
+            Info<< "Background regions:" << endl;
+            Info<< token::TAB << "Region"
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            Map<scalar>::const_iterator vIter = allRegionVolume.begin();
+            Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
 
-            label patchI = mesh.boundaryMesh().findPatchID(pName);
-            if (patchI != -1)
-            {
-                const polyPatch& pp = mesh.boundaryMesh()[patchI];
-
-                // Collect all regions on the patch
-                const labelList& faceCells = pp.faceCells();
-
-                forAll(faceCells, i)
-                {
-                    keepRegions.insert
-                    (
-                        regions[faceCells[i]],
-                        Pstream::myProcNo()
-                    );
-                }
-            }
-        }
-
-
-        // Make sure all the processors have the same set of regions
-        Pstream::mapCombineGather(keepRegions, minEqOp<label>());
-        Pstream::mapCombineScatter(keepRegions);
-
-        Info<< "Patch connected regions (liquid core):" << endl;
-        forAllConstIter(Map<label>, keepRegions, iter)
-        {
-            label regionI = iter.key();
-            Pair<scalar>& vols = regionVolume[regionI];
-            Info<< token::TAB << iter.key()
-                << token::TAB << vols.first()
-                << token::TAB << vols.second() << endl;
-
-        }
-        Info<< endl;
-
-        Info<< "Background regions:" << endl;
-        forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
-        {
-            if
+            for
             (
-               !keepRegions.found(iter.key())
-             && iter().first() >= volFraction_*meshVol
+                ;
+                vIter != allRegionVolume.end()
+             && aIter != allRegionAlphaVolume.end();
+                ++vIter, ++aIter
             )
             {
-                Info<< token::TAB << iter.key()
-                    << token::TAB << iter().first()
-                    << token::TAB << iter().second() << endl;
+                if
+                (
+                   !patchRegions.found(vIter.key())
+                 && vIter() >= maxDropletVol
+                )
+                {
+                    Info<< token::TAB << vIter.key()
+                        << token::TAB << vIter()
+                        << token::TAB << aIter() << endl;
+                }
             }
+            Info<< endl;
         }
-        Info<< endl;
+
 
 
         // Split alpha field
@@ -400,185 +658,197 @@ void Foam::regionSizeDistribution::write()
         //  - liquidCore            : region connected to inlet patches
         //  - per region a volume   : for all other regions
         //  - backgroundAlpha       : remaining alpha
+        writeAlphaFields(regions, patchRegions, allRegionVolume, alpha);
 
 
-        // Construct field
-        volScalarField liquidCore
-        (
-            IOobject
-            (
-                alphaName_ + "_liquidCore",
-                obr_.time().timeName(),
-                obr_,
-                IOobject::NO_READ
-            ),
-            alpha,
-            fvPatchField<scalar>::calculatedType()
-        );
-
-        volScalarField backgroundAlpha
-        (
-            IOobject
-            (
-                alphaName_ + "_background",
-                obr_.time().timeName(),
-                obr_,
-                IOobject::NO_READ
-            ),
-            alpha,
-            fvPatchField<scalar>::calculatedType()
-        );
-
-
-        // Knock out any cell not in keepRegions
-        forAll(liquidCore, cellI)
+        // Extract droplet-only allRegionVolume, i.e. delete liquid core
+        // (patchRegions) and background regions from maps.
+        // Note that we have to use mesh volume (allRegionVolume) and not
+        // allRegionAlphaVolume since background might not have alpha in it.
+        forAllIter(Map<scalar>, allRegionVolume, vIter)
         {
-            label regionI = regions[cellI];
-            if (keepRegions.found(regionI))
+            label regionI = vIter.key();
+            if
+            (
+                patchRegions.found(regionI)
+             || vIter() >= maxDropletVol
+            )
             {
-                backgroundAlpha[cellI] = 0;
-            }
-            else
-            {
-                liquidCore[cellI] = 0;
-
-                scalar regionVol = regionVolume[regionI].first();
-                if (regionVol < volFraction_*meshVol)
-                {
-                    backgroundAlpha[cellI] = 0;
-                }
+                allRegionVolume.erase(vIter);
+                allRegionAlphaVolume.erase(regionI);
+                allRegionNumCells.erase(regionI);
             }
         }
-        liquidCore.correctBoundaryConditions();
-        backgroundAlpha.correctBoundaryConditions();
 
-        Info<< "Volume of liquid-core = "
-            << fvc::domainIntegrate(liquidCore).value()
-            << endl;
-
-        Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
-        liquidCore.write();
-
-        Info<< "Volume of background = "
-            << fvc::domainIntegrate(backgroundAlpha).value()
-            << endl;
-
-        Info<< "Writing background field to " << backgroundAlpha.name() << endl;
-        backgroundAlpha.write();
-
-
-
-        // Collect histogram
-        if (Pstream::master())
+        if (allRegionVolume.size())
         {
-            DynamicList<scalar> diameters(regionVolume.size());
-            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            // Construct mids of bins for plotting
+            pointField xBin(nBins_);
+
+            scalar x = 0.5*delta;
+            forAll(xBin, i)
             {
-                if (!keepRegions.found(iter.key()))
-                {
-                    if (iter().first() < volFraction_*meshVol)
-                    {
-                        scalar v = iter().second();
-                      //scalar diam = Foam::cbrt(v*6/mathematicalConstant::pi);
-                        scalar diam =
-                            Foam::cbrt(v*6/constant::mathematical::pi);
-                        diameters.append(diam);
-                    }
-                }
+                xBin[i] = point(x, 0, 0);
+                x += delta;
             }
 
-            if (diameters.size())
+            const coordSet coords("diameter", "x", xBin, mag(xBin));
+
+
+            // Get in region order the alpha*volume and diameter
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList sortedRegions = allRegionAlphaVolume.sortedToc();
+
+            scalarField sortedVols
+            (
+                extractData
+                (
+                    sortedRegions,
+                    allRegionAlphaVolume
+                )
+            );
+
+            // Calculate the diameters
+            scalarField sortedDiameters(sortedVols.size());
+            forAll(sortedDiameters, i)
             {
-                scalar maxDiam = max(diameters);
-                scalar minDiam = 0.0;
-
-                Info<< "Maximum diameter:" << maxDiam << endl;
-
-                Histogram<List<scalar> > bins
+                sortedDiameters[i] = Foam::cbrt
                 (
-                    minDiam,
-                    maxDiam,
-                    nBins_,
-                    diameters
+                    sortedVols[i]
+                   *6/constant::mathematical::pi
                 );
+            }
 
-                /* 1.7.x
-                scalarField xBin(nBins_);
+            // Determine the bin index for all the diameters
+            labelList indices(sortedDiameters.size());
+            forAll(sortedDiameters, i)
+            {
+                indices[i] = (sortedDiameters[i]-minDiam_)/delta;
+            }
 
-                scalar dx = (maxDiam-minDiam)/nBins_;
-                scalar x = 0.5*dx;
-                forAll(bins.counts(), i)
+            // Calculate the counts per diameter bin
+            scalarField binCount(nBins_, 0.0);
+            forAll(sortedDiameters, i)
+            {
+                binCount[indices[i]] += 1.0;
+            }
+
+            // Write counts
+            if (Pstream::master())
+            {
+                writeGraph(coords, "count", binCount);
+            }
+
+            // Write to screen
+            {
+                Info<< "Bins:" << endl;
+                Info<< token::TAB << "Bin"
+                    << token::TAB << "Min diameter"
+                    << token::TAB << "Count:"
+                    << endl;
+
+                scalar diam = 0.0;
+                forAll(binCount, binI)
                 {
-                    xBin[i] = x;
-                    x += dx;
+                    Info<< token::TAB << binI
+                        << token::TAB << diam
+                        << token::TAB << binCount[binI] << endl;
+                    diam += delta;
                 }
+                Info<< endl;
+            }
 
-                scalarField normalisedCount(bins.counts().size());
-                forAll(bins.counts(), i)
+
+            // Write average and deviation of droplet volume.
+            writeGraphs
+            (
+                "volume",           // name of field
+                indices,            // per region the bin index
+                sortedVols,         // per region field data
+                binCount,           // per bin number of regions
+                coords              // graph data for bins
+            );
+
+            // Collect some more field
+            {
+                wordList scalarNames(obr_.names(volScalarField::typeName));
+                labelList selected = findStrings(fields_, scalarNames);
+
+                forAll(selected, i)
                 {
-                    normalisedCount[i] = 1.0*bins.counts()[i];
+                    const word& fldName = scalarNames[selected[i]];
+                    Info<< "Scalar field " << fldName << endl;
+
+                    const scalarField& fld = obr_.lookupObject
+                    <
+                        volScalarField
+                    >(fldName).internalField();
+
+                    writeGraphs
+                    (
+                        fldName,            // name of field
+                        alphaVol*fld,       // per cell field data
+
+                        regions,            // per cell the region(=droplet)
+                        sortedRegions,      // valid regions in sorted order
+                        1.0/sortedVols,     // per region normalisation
+
+                        indices,            // index of bin for each region
+                        binCount,           // per bin number of regions
+                        coords              // graph data for bins
+                    );
                 }
+            }
+            {
+                wordList vectorNames(obr_.names(volVectorField::typeName));
+                labelList selected = findStrings(fields_, vectorNames);
 
-                const coordSet coords
-                (
-                    "diameter",
-                    "x",
-                    xBin
-                );
-                */
-
-                pointField xBin(nBins_);
-                scalar dx = (maxDiam - minDiam)/nBins_;
-                scalar x = 0.5*dx;
-                forAll(bins.counts(), i)
+                forAll(selected, i)
                 {
-                    xBin[i] = point(x, 0, 0);
-                    x += dx;
+                    const word& fldName = vectorNames[selected[i]];
+                    Info<< "Vector field " << fldName << endl;
+
+                    const vectorField& fld = obr_.lookupObject
+                    <
+                        volVectorField
+                    >(fldName).internalField();
+
+
+                    // Components
+
+                    for (direction cmp = 0; cmp < vector::nComponents; cmp++)
+                    {
+                        writeGraphs
+                        (
+                            fldName + vector::componentNames[cmp],
+                            alphaVol*fld.component(cmp),// per cell field data
+
+                            regions,        // per cell the region(=droplet)
+                            sortedRegions,  // valid regions in sorted order
+                            1.0/sortedVols, // per region normalisation
+
+                            indices,        // index of bin for each region
+                            binCount,       // per bin number of regions
+                            coords          // graph data for bins
+                        );
+                    }
+
+                    // Magnitude
+                    writeGraphs
+                    (
+                        fldName + "mag",    // name of field
+                        alphaVol*mag(fld),  // per cell field data
+
+                        regions,            // per cell the region(=droplet)
+                        sortedRegions,      // valid regions in sorted order
+                        1.0/sortedVols,     // per region normalisation
+
+                        indices,            // index of bin for each region
+                        binCount,           // per bin number of regions
+                        coords              // graph data for bins
+                    );
                 }
-
-                scalarField normalisedCount(bins.counts().size());
-                forAll(bins.counts(), i)
-                {
-                    normalisedCount[i] = 1.0*bins.counts()[i];
-                }
-
-                const coordSet coords
-                (
-                    "diameter",
-                    "x",
-                    xBin,
-                    mag(xBin)
-                );
-                const wordList valNames(1, "count");
-
-
-                fileName outputPath;
-                if (Pstream::parRun())
-                {
-                    outputPath = mesh.time().path()/".."/name_;
-                }
-                else
-                {
-                    outputPath = mesh.time().path()/name_;
-                }
-
-                if (mesh.name() != fvMesh::defaultRegion)
-                {
-                    outputPath = outputPath/mesh.name();
-                }
-
-                mkDir(outputPath/mesh.time().timeName());
-                OFstream str
-                (
-                    outputPath
-                  / mesh.time().timeName()
-                  / formatterPtr_().getFileName(coords, valNames)
-                );
-                Info<< "Writing distribution to " << str.name() << endl;
-
-                List<const scalarField*> valPtrs(1);
-                valPtrs[0] = &normalisedCount;
-                formatterPtr_().write(coords, valNames, valPtrs, str);
             }
         }
     }
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
index 033b5f9aec..7802733249 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
@@ -25,11 +25,67 @@ Class
     Foam::regionSizeDistribution
 
 Description
-    Looks up a field, interpolates it to the faces and determines a connected
-    region from a patch where the field is above a certain value.
-    - Writes a field containing all regions starting at given patch
-      ('liquid core')
-    - All other regions are summed for volume and a histogram is calculated.
+    Droplet size distribution calculation.
+
+    Looks up a void-fraction (alpha) field and splits the mesh into regions
+    based on where the field is below the threshold value. These
+    regions ("droplets") can now be analysed.
+
+    Regions:
+    - (debug) write regions as a volScalarField
+    - (debug) print for all regions the sum of volume and alpha*volume
+    - print the regions connected to a user-defined set of patches.
+      (in spray calculation these form the liquid core)
+    - print the regions with too large volume. These are the 'background'
+      regions.
+
+    Fields:
+    - write volScalarField alpha_liquidCore : alpha with outside liquid core
+                                              set to 0.
+                           alpha_background : alpha with outside background
+                                              set to 0.
+
+    Histogram:
+    - determine histogram of diameter (given minDiameter, maxDiameter, nBins)
+    - write graph of number of droplets per bin
+    - write graph of sum, average and deviation of droplet volume per bin
+    - write graph of sum, average and deviation of user-defined fields. For
+      volVectorFields these are those of the 3 components and the magnitude.
+
+    Sample input:
+
+    functions
+    {
+        regionSizeDistribution
+        {
+            type            regionSizeDistribution;
+
+            outputControl   timeStep;
+            outputInterval  1;
+
+            // Field to determine regions from
+            field           alpha;
+            // Patches that provide the liquid core
+            patches         (inlet);
+            // Delimit alpha regions
+            threshold       0.4;
+
+            // Fields to sample (no need to include alpha)
+            fields          (p U);
+
+            // Number of bins for histogram
+            nBins           100;
+            // Max droplet diameter
+            maxDiameter     0.5e-4;
+            //// Min droplet diameter (default is 0)
+            //minDiameter     0;
+
+            // Writing format
+            setFormat       gnuplot;
+        }
+    }
+
+
 
 SourceFiles
     regionSizeDistribution.C
@@ -41,6 +97,9 @@ SourceFiles
 
 #include "pointFieldFwd.H"
 #include "writer.H"
+#include "Map.H"
+#include "volFieldsFwd.H"
+#include "wordReList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +110,8 @@ namespace Foam
 class objectRegistry;
 class dictionary;
 class mapPolyMesh;
+class regionSplit;
+class polyMesh;
 
 /*---------------------------------------------------------------------------*\
                            Class regionSizeDistribution Declaration
@@ -72,23 +133,85 @@ class regionSizeDistribution
         word alphaName_;
 
         //- Patches to walk from
-        wordList patchNames_;
+        wordReList patchNames_;
 
         //- Clip value
         scalar threshold_;
 
-        //- Background region volFraction
-        scalar volFraction_;
+        //- Maximum droplet diameter
+        scalar maxDiam_;
+
+        //- Minimum droplet diameter
+        scalar minDiam_;
 
         //- Mumber of bins
         label nBins_;
 
+        //- Names of fields to sample on regions
+        wordReList fields_;
+
         //- Output formatter to write
         autoPtr<writer<scalar> > formatterPtr_;
 
 
     // Private Member Functions
 
+        template<class Type>
+        Map<Type> regionSum(const regionSplit&, const Field<Type>&) const;
+
+        //- Get data in order
+        template<class Type>
+        List<Type> extractData(const UList<label>& keys, const Map<Type>&)
+        const;
+
+        void writeGraph
+        (
+            const coordSet& coords,
+            const word& valueName,
+            const scalarField& values
+        ) const;
+
+        //- Write volfields with the parts of alpha which are not
+        //  droplets (liquidCore, backGround)
+        void writeAlphaFields
+        (
+            const regionSplit& regions,
+            const Map<label>& keepRegions,
+            const Map<scalar>& regionVolume,
+            const volScalarField& alpha
+        ) const;
+
+        //- Mark all regions starting at patches
+        Map<label> findPatchRegions(const polyMesh&, const regionSplit&) const;
+
+        //- Helper: divide if denom != 0
+        static tmp<scalarField> divide(const scalarField&, const scalarField&);
+
+        //- Given per-region data calculate per-bin average/deviation and graph
+        void writeGraphs
+        (
+            const word& fieldName,              // name of field
+            const labelList& indices,           // index of bin for each region
+            const scalarField& sortedField,     // per region field data
+            const scalarField& binCount,        // per bin number of regions
+            const coordSet& coords              // graph data for bins
+        ) const;
+
+        //- Given per-cell data calculate per-bin average/deviation and graph
+        void writeGraphs
+        (
+            const word& fieldName,              // name of field
+            const scalarField& cellField,       // per cell field data
+
+            const regionSplit& regions,         // per cell the region(=droplet)
+            const labelList& sortedRegions,     // valid regions in sorted order
+            const scalarField& sortedNormalisation,
+
+            const labelList& indices,           // index of bin for each region
+            const scalarField& binCount,        // per bin number of regions
+            const coordSet& coords              // graph data for bins
+        ) const;
+
         //- Disallow default bitwise copy construct
         regionSizeDistribution(const regionSizeDistribution&);
 
@@ -156,6 +279,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "regionSizeDistributionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C
new file mode 100644
index 0000000000..2f55434315
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSizeDistribution.H"
+#include "regionSplit.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Map<Type> Foam::regionSizeDistribution::regionSum
+(
+    const regionSplit& regions,
+    const Field<Type>& fld
+) const
+{
+    // Per region the sum of fld
+    Map<Type> regionToSum(regions.nRegions()/Pstream::nProcs());
+
+    forAll(fld, cellI)
+    {
+        label regionI = regions[cellI];
+
+        typename Map<Type>::iterator fnd = regionToSum.find(regionI);
+        if (fnd == regionToSum.end())
+        {
+            regionToSum.insert(regionI, fld[cellI]);
+        }
+        else
+        {
+            fnd() += fld[cellI];
+        }
+    }
+    Pstream::mapCombineGather(regionToSum, plusEqOp<Type>());
+    Pstream::mapCombineScatter(regionToSum);
+
+    return regionToSum;
+}
+
+
+// Get data in sortedToc order
+template<class Type>
+Foam::List<Type> Foam::regionSizeDistribution::extractData
+(
+    const UList<label>& keys,
+    const Map<Type>& regionData
+) const
+{
+    List<Type> sortedData(keys.size());
+
+    forAll(keys, i)
+    {
+        sortedData[i] = regionData[keys[i]];
+    }
+    return sortedData;
+}
+
+
+// ************************************************************************* //