As per my understanding, for 1d bins in z direction for a specified region, the bin volume should be dz*(xhi-xlo)*(yhi-ylo), where lo and hi are from the region specifications.

However, after a quick peek into the source code, I guess the the bin volume, for 1-d bins in z dim, is computed based on the xy surface area of the box, and it does not consider any region condition.

So to get desired densities in the specified region, I should scale the densities as rho(z) * xy surface are of the box / xy surface area of the region. With this, the density values that I get are in the expected value range.

I would appreciate if someone could comment on this.

There is no logic in the code to calculate the intersection
of a region with a 1d slab, 2d pencil, 3d brick to get what
you are thinking of as the “exact” volume. That would be
very hard to do for arbitrary regions, even a sphere is not
simple. If the region is a union/intersection, then much
harder.

We should put a comment on the doc page about this …

There is no logic in the code to calculate the intersection
of a region with a 1d slab, 2d pencil, 3d brick to get what
you are thinking of as the “exact” volume. That would be
very hard to do for arbitrary regions, even a sphere is not
simple. If the region is a union/intersection, then much
harder.

While doing and exact determination is hard, getting a decent approximation should be straightforward via monte carlo integration. Since binned data has large statistical error bars anyway, it would not be as bad even with a small number of samples.

I guess, since determining volume of a generic region may not be straightforward, in that case fix ave/chunk … density/* … should only write average particle number or mass in the bins and let user scale it with appropriate volume.