Is there any neat way to calculate this contour integral $PV\int_{\mathbb{R}^n}\frac{e^{i(b_1z_1+\cdots+b_nz_n)}}{\prod_{j$$\mathrm {PV}\int_{\mathbb{R}^n} \frac{e^{i(b_1z_1+\cdots+b_nz_n)}}{\prod_{j<k}(z_k-z_j)}dz=?$$
Other than integrate this term by term? It is in fact the Fourier transform of the inverse Vandermonde determinant.
I also found this question:
Fourier transform of the inverse of Vandermonde determinant
which is asking for the $n=3$ case.

So finally you meant the Fourier transform of the distribution $$F(x) = PV(\frac{1}{\prod_{j < k} (x_k-x_j)})$$
that is $$\hat{F}(\omega) = \int_{\mathbb{R}^n} F(x) e^{-i \langle \omega,x \rangle}d^n x$$
Then 
$$1=F(x)\prod_{j < k}  (x_k-x_j)$$
implies $$(2\pi)^n i^{n(n-1)/2} \delta(\omega)=\prod_{j < k}(\partial_j-\partial_k) \widehat{F}(\omega) $$
so that, with $\ast$ the convolution and $T_{j,k}$ the Dirac delta distribution on the subspace orthogonal to $x_j-x_k$ $$\widehat{F}(\omega)  = (2\pi)^n i^{n(n-1)/2} \quad T_{1,2} 1_{x_1-x_2 > 0} \ast \ldots \ast T_{j,k}  1_{x_j-x_k > 0} \ast \ldots \ast  T_{n-1,n} 1_{x_{n-1}-x_n > 0}  + P(\omega)$$
for some polynomial $P$ such that $(\prod_{j < k}(\partial_j-\partial_k))P = 0$, and we may deduce the coefficients of $P$ from symmetry considerations.
