General Idea
The built-in function mathutils.geometry.box_fit_2d
could help you with that.

from left to right: original, without correction, with correction
box_fit_2d
is a function that calculates the angle by which two-dimensional coordinates have to be rotated to best fit an axis-aligned rectangle.
For every object you take a list of all vertices and project them to a plane perpendicular to the axis around which you want to rotate. Then you convert these projected vertices' coordinates, which are still three-dimensional, into two-dimensional coordinates so that box_fit_2d
accepts them. If you rotate around one of the main axes then this is simply a matter of dropping one of the coordinates. E.g. vertex [x, y, z] becomes [y, z] when you rotate around the x axis.
This two-dimensional vector is passed to the function which returns the angle that you can use to rotate the object.
As far as I'm aware, box_fit_2d
doesn't do checks about which side is longer. So, if you want the bars to point in the same direction, you have to do an additional check, e.g. whether the range of y values of the rotated object's vertices is higher than the range of z values, like it's the case for the box in front in the image in the middle. And if that's the case increase the rotation angle by 90 degrees.
Implementation
I generally use the add-on Animation Nodes for this kind of calculation, but this method could be done with a Python script in Blender without any additional tools.
The naïve Animation Nodes node tree that was used to generate the middle image is pictured below. Note that you have to would have to list the module mathutils
in the advanced node settings for the expression node, if you want to use the module in an expression.


The correction that would increase the rotational angle by another 90 degrees when the object is horizontally aligned was wired like below:
