Complex output after inverse FFT of a real signal
Essentially, your code does not respect the inherent Hermitian symmetry of the output of the FFT. Here, your signal is odd-sized $2K+1$. Hence, this FFT yields a complex vector of coefficients $d$ (real) and $a_k$ (generally complex), arranged as:

$$ \left[d,a_1,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},\overline{a_1} \right]$$

If you want to zero-out some coefficients, you should remove then by conjugate pairs (except for $d$), like in

$$ \left[0,0,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},0 \right].$$

Note: for even sizes, a real coefficient $n_{K+1}$ is inserted between $a_K$ and $\overline{a_K}$.

Your line:

X([1:6, 4780:4801]) = 0


does not satisfy this symmetry. Try instead this one (only valid for odd length):

nACBin = 3;
X([1:nACBin+1, end-nACBin+1:end]) = 0; % 


and you will get a much lower imaginary residue (here due to numerical issues). Here, it somehow works (see below) because of the very simple context (two sines with lucky interval length), but this is rarely recommended in general: 



Especially with short 0 strides, this will convolve your data with a cardinal sine, both broad and wiggling.

As said in comments, you are considering prototyping a filter in the frequency domain using windows. This is a well-covered domain in DSP. However you can avoid those complexities by first looking at filtering in the time domain (Butterworth, Chebyshev, Cauer, Papoulis filters) that already have well-designed implementations. Other filters can be useful for peaks, like polynomial interpolators (Savitzky-Golay, for instance).

Side note: it probably would be more efficient (though a bit more complicated) to embed those frequency constraints into the deconvolution, instead of trying to remove the artifacts afterward.