Bottom Line of the Qualitative Explanation: The Fourier Transform is a correlation of our arbitrary signal $x(t)$ with all frequencies, with each frequency given as a complex exponential: $e^{j\omega t}$ (this is a rotating phasor of magnitude $1$ and angle $\omega t$). This holds for $x(t)$ as being real as in the OP's case or complex. The general form for a correlation computation (that we see in so many places in signal processing: matched filter receiver, Viterbi decoding, and yes the Fourier Transform!) is the integration of a complex conjugate product (or sum of products in discrete signals), and such a computation will beautifully have minimum noise in the presence of white stationary noise. This is explained in more detail below with my hope of providing a more intuitive understanding. Such a correlation will maximize the estimate of frequency and phase for each $\omega$ in the presence of white noise (but may not necessarily provide an optimum estimate in the presence of other noise types).
Details:
I actually find the Fourier Transform described with sines and cosines much more cumbersome and the expression with the exponentials quite intuitive once the mystery of the form $e^{j\omega t}$ is resolved.
Here are a few main take-aways:
First and foremost, it may not be apparent that the expression $Ke^{j\phi}$ is just a complex phasor with real magnitude $K$ and real angle $\phi$. So $Ke^{j\phi} = K\angle\phi$.
Next, and also important, consider a single Fourier frequency tone as a single rotating phasor, not a sinusoid, rotating at a constant rate with a constant magnitude. Understanding this makes signal processing a lot more intuitive in my opinion. Similar to a bicycle wheel rotating where the frequency is it's rate of rotation in cycles per second; we will also get signed frequency to denote the direction of rotation: If the phasor is rotating counter-clockwise; this is a positive frequency. If the phasor is rotating clockwise, this is a negative frequency. This is consistent with the definition of instantaneous frequency as the time derivative of phase: a change in phase versus a change in time.
Below is a plot demonstrating this where we see a phasor rotating counter-clockwise on a complex plane ("I" is the in-phase or real axis and "Q" is the quadrature or imaginary axis).
Now Euler's relationship between $\cos(\omega t)$, $\sin(\omega t)$ and $e^{j \omega t}$ should make a little more sense; for example we can now see how:
$$2\cos(\omega t) = e^{j\omega t} + e^{-j \omega t}$$
These are just two phasors rotating with equal rates and opposite phase (a positive and negative frequency) that are added together. When you add two phasors, the angles cancel, or graphically place on phasor on the end of the other as shown in green below. In either case the result is the cosine which is always real (oscillates back and forth along the real axis).
The general formula for a correlation measure of two functions is the integration of the complex conjugate product:
$$corr = \int f_1(t)f_2^*(t)dt$$
Conjugation means invert the phase and is denoted by "*" as in $x^*(t)$ is the conjugate of $x(t)$)For example, the two phasors as shown above for the cosine are in a conjugate relationship. One phasor is the conjugate of the other.
The correlation operation shown above has the wonderful property that "like" signals will grow in the integration to relatively large quantities while noise will reduce significantly through averaging. To see how "like" signals grow: consider a signal that is rotating in phase (a constant frequency); if we multiply it by its complex conjugate we will de-rotate it back to the positive real axis, and thus keep growing along the real axis in the integration (since $e^{j\omega t}e^{-j \omega t} =1$. As another example consider a sequence: 1 -1 1 1 -1, if we multiply it with itself, all samples will become 1, and the sum (integration) will grow to 5. Similarly if each sample was complex and also had a phase angle, if we multiplied it by it's complex conjugate, we would get the same result of each product being 1 and the sum growing to 5. Any other sequences with the same magnitude on each sample but other arbitrary angles will not grow as large.
We do the same with the Fourier Transform:
$$F(\omega) = \int x(t)e^{-j \omega t}dt$$
For each single frequency $\omega$, resulting in a spinning phasor given by $e^{j \omega t}$, we multiply our arbitrary waveform x(t) by the complex conjugate $e^{-j \omega t}$, and integrate: we are correlating our waveform to each each single frequency given by $e^{j\omega t}$ and through that determining it's magnitude and phase within $x(t)$! As we repeat that sweeping through each $\omega$ of interest, we get the complete spectrum of all frequencies. This is the Fourier Transform.
The OP then went on to ask, and again in the comments why choosing the same $\omega$ will maximize the result for the given phase angle versus any other $\omega$. Consider the full relationship given by Euler as detailed in the plot below as a continuation to my attempts to make $e^{j\omega t}$ more intuitive:
So when we correlate to $e^{j \omega t}$, we are correlating to both the real and imaginary terms (cosine and sine), resulting in a complex result, whose angle will be the angle of that tone in $x(t)$!
$$F(\omega) = \int x(t)e^{-j \omega t}dt$$
$$= \int (x_i(t)+jx_q(t))(\cos (\omega t) + j\sin(\omega t)dt $$
If we multiply out the real and imaginary terms, we see how when a real tone is shifted by any phase angle ($cos(\omega_o t + \phi)$, correlating it with $e^{j\omega t}$ will have a maximum magnitude when $\omega = \omega_o$ for any $\phi$, and that the resulting complex output of the correlation (the real and imaginary terms above) will be at that same angle!
For the case of the OP's question where $x(t)$ is real the above reduces to:
$$F(\omega) = \int (\cos(\omega t - \phi))(\cos (\omega t) + j\sin(\omega t)dt $$
$$ = \int \cos(\omega t - \phi)\cos (\omega t) dt + j \int \cos(\omega t - \phi)\sin (\omega t) dt$$
I prefer to not resort to the sine and cosine representation even in my own head to make sense of what is going on, but to just consider $e^{j\omega_o t}$ on its own as a single spinning phasor. If our signal $x(t)$ was also just a single spinning phasor also at $e^{j\omega_o t}$, then when we do the complex conjugate product: $x(t)e^{-j\omega_o t} = e^{j\omega_o t}e^{-j\omega_o t} = e^{j0} = 1$!! I show this graphically in the link at the bottom, but with no added phase to $x(t)$ we rotated to the real axis $(1)$ at all samples, which we integrate over time period $T$ leading to the maximum result of $T$ with no imaginary term, so it is $T\angle 0$. If we then add a phase rotation to $x(t)$ to get $e^{j\omega_o t + \phi}$ The product becomes $x(t)e^{-j\omega_o t} = e^{j\omega_o t + \phi}e^{-j\omega_o t} = e^{j\phi} = 1\angle \phi$. This would also grow to the same magnitude $T$ but with the angle $\phi$ Thus our output is complex with real and imaginary terms, that we could also get to the long way (in our head) with sines and cosines. You can carry this same thing out to real signals ($x(t)= cos(\omega t+\phi)$ by doing what I did above with both positive and negative frequency tones ($e^{j\omega_o t + \phi}$ amd $e^{-(j\omega_o t + \phi)}$). I show this all graphically and intuitively in the link given below.
To then see further details in how the correlation process itself will maximize the SNR of the $e^{j\omega t}$ signal (this is true under the condition of white noise), please see my response in this post, which includes showing how the correlation falls off as we deviate slightly from $\omega_o$:
Derivation of the Optimal Matched Filter - Convolution vs. Correlation
