how to get the phasor(frequency domain) representation of a digitalized signal?
This is simply about representing a sinusoidal signal $x(t)$ in the following way:

$$x(t) \ = \ \cos(\omega_0t+\phi) \ = \ a \cos(\omega_0 t) \ - \ b \sin(\omega_0 t)$$

The 'method' described in your question is about identifying the constants $a$ and $b$. The phasor of $x(t)$ is then $X=a+jb$, because

$$x(t)=\text{Re}\{Xe^{j\omega_0t}\}$$

Using basic trigonometric identities you get

$$a=\cos\phi\quad\text{and}\quad b=\sin\phi$$

But you can also compute these values by correlation (as shown in your question), which is equivalent to computing the Fourier coefficients of $x(t)$ (apart from a sign difference for the coefficient $b$). Of course, the only non-zero coefficients are the coefficients corresponding to the fundamental:

$$a=\frac{2}{T}\int_0^Tx(t)\cos(\omega_0 t) \ dt,\quad  T=\frac{2\pi}{\omega_0}\\
b=-\frac{2}{T}\int_0^Tx(t)\sin(\omega_0 t) \ dt$$

These last two equations correspond to the first two steps described in your question.