Hi Phillippe,
If you want to get some sense of the model fitting across all the subjects, you could convert the full F to zscore (normal distribution) by running 3dmerge with option -lzscore, and then do group analysis with 3dttest on the zscore values by testing whether the zscore is signicificant at each voxel. I think this could work for your purpose.
As I said before, converting full F to t is problematic because it is hard to decide what sign you would assign to t with t = +/- F^(1/2) at each voxel. And also 3dttest runs on data presumably from normal distribution, but not from t.
Your question No. 1: The assumption of normality is about the noise at each voxel, but not the real signal across all voxels in the brain. Theoretically if you model the real signal at each voxel with all the necessary components (drifting baseline and all stimuli), the leftover part would be white noise, and such noise is independent across all voxels. Probably you could test whether the residuals at each voxel are N(0, 1)-distributed or not.
3dhistog would give you the histogram of a subbrick in a dataset, thus it is the distribution of that quantity (i.e., F, t, zscore, etc.) across all the voxels. But this is not what you want to test. Instead, it is the residual at each voxel you would want to check its normality.
Question No. 2: I am not so sure whether there is such a plugin available or not. Hope somebody else could chime in about this.
Gang