std::piecewise_constant_distribution gives wrong values? [BUG?]

I have the following strange phenomenon which puzzles me!:

I have a piecewise constant probability density given as

 using RandomGenType = std::mt19937_64;
RandomGenType gen(1230123);

std::array<double,5> intervals {0.59e-3, 0.7e-3, 0.85e-3, 1e-3, 1.18e-3};
std::array<double,4> weights {1.36815, 1.99139, 0.29116, 0.039562};
std::piecewise_constant_distribution<double>
distribution (intervals.begin(),intervals.end(),weights.begin());

When I draw n random numbers (radius of sphere) from this distribution and compute the mass of the sphere and sum them up like:

unsigned int n = 1000000;
double density = 2400;
double mass = 0;

for(int i=0;i<n;i++){
    auto d = 2* distribution(gen);
    mass += d*d*d/3.0*M_PI_2*density;
}

I get mass = 4.3283 kg (see LIVE here)

Doing the EXACT identical thing in Mathematica like:

Gives the assumably correct value of 4.5287 kg.

Which is not the same, also with different seeds I never get the same? Question : What the hack is wrong with the sampling in C++?

Simple Mathematica Code:

df[r_] := Piecewise[{{0, r <= 0.59}, {1.36814, 
      Inequality[0.59, Less, r, 
               LessEqual, 0.7]}, {0, r > 0.7}}, Indeterminate] + 
       Piecewise[{{0, r <= 0.7}, {1.99139, 
      Inequality[0.7, Less, r, LessEqual, 
               0.85]}, {0, r > 0.85}}, Indeterminate] + 
       Piecewise[{{0., r <= 0.85}, {0.29116, 
      Inequality[0.85, Less, r, LessEqual, 
               1]}, {0., r > 1}}, Indeterminate] + 
       Piecewise[{{0, r <= 1}, {0.039562, 
      Inequality[1, Less, r, LessEqual, 1.18]}, 
           {0, r > 1.18}}, Indeterminate]; 

 PDFr = ProbabilityDistribution[pdf[r], {r, 0.59, 1.18}];
 dataR = RandomVariate[PDFr, 1000000];
 Fold[#1 + (2*#2*10^-3)^3 (\[Rho] \[Pi]/2)/3 &, 0, dataR] /. \[Rho] -> 
 2400