Creating a continuous evaluation of a phase graph

I have a function

$$R(\xi) := \prod_{i=0}^{n} (e^{-i\xi} - z_k)$$ where $$z_{k} \in \mathbb{C}$$. The $$z_{k}$$ are random, but clustered vaguely around the unit circle, if that is relevant. I need to calculate the phase $$\Psi$$ of $$R$$, where we define $$\Psi$$ via $$R(\xi) = \exp(-i\Psi(\xi))|R(\xi)|$$ If we do not care about continuity of $$\Psi$$ in $$\xi$$, then this is easy: We could write (in C++)

auto Psi = [&](Real xi)->Real {
Real phi = 0;
Complex ii = {0,1};
Complex z = exp(-ii*xi);
for (zk in roots)
{
phi -= arg(z-zk);
}
return phi;
};

However, I have an additional constraint that the phase $$\Psi$$ must be continuous on the interval $$[0, 2\pi)$$. The arg function always returns a value in $$[-\pi, \pi]$$, and as such the phase has discontinuities at seemingly random locations, depending on the values of the roots $$\{z_{k}\}$$.

How can I patch up my code/figure out how to make $$\Psi$$ a continuous function of $$\xi$$ on the interval $$[0, 2\pi)$$?