Accurate description of the spectrum of sinusoid waveform added with WGN?
Notice that in your first image, there are peaks where they would be expected to be, so that's fine. I guess you think that "the spectrum obtained through MATLAB seems not like the theroetical one" because of the noisy spectrum you get, instead of a straight line at $\sigma^2$ on the vertical axis.

The random signal added to $s(t)$ is not random actually, it's just a realization you get when calling awgn(). So the FFT will return just a noisy spectrum for that signal. If you generate lots of realizations of $r(t)$ and plot the average of all their FFTs, you'll obtain what you were expecting for.

For your second question on how to calculate the threshold SNR, some criteria is needed. If you average lots of FFTs the variance of the noisy part of the spectrum will drastically decrease, so the SNR needed to see the peaks will be much smaller. If you use just one realization, then the SNR will be higher. There is not a fixed number for the needed SNR, it depends on what you refer to "sinusoid covered by noise". If the spectrum observed was just like the theoretical one (i.e. you averaged a lot of FFTs), a sinusoid would be "completely covered by noise" when its amplitude is such that the corresponding peak in the frequency domain is below $\sigma^2$.