First, as Dave Tweed has answered, you generally use a lock-in to recover a small signal buried in noise.
That said, your script is not properly implementing a lock-in amplifier, as evidenced by your second trace. Your problem is that the DC component of your original signal needs to be suppressed (the signal should be AC-coupled). If your reference sine wave has a DC component of zero (which it should) then for a signal with zero degrees of phase shift and an average of zero, the output will be a sine-squared wave (plus noise). Note that this will be rectified, with no signal component negative.This will allow a low-pass filter to recover the amplitude of the desired frequency, but not its shape.
What you seem to be trying to do is simple noise rejection, and there are two possibilities. Either your noise is broadband, with significant noise energy both above and below your frequency of interest, or the noise is only significant above your fundamental.
Assuming the latter, you can process your signal using only a high-pass filter, made arbitrarily sharp and close to your fundamental. If the former, you need a bandpass filter.
In either case, looking at the crossover distortion shown in the last figure is going to be very, very difficult. That's a low-energy, high-frequency artifact, and may not be easily recovered from the noise. If you really want to try, the first thing you need to do is simulate your signal, then perform an FFT on it to establish the frequency response your filter needs in order not to exclude the signal of interest. Then compare this to the noise spectrum and you'll probably see that they overlap.
Other than an extremely large averaging filter (many, many waveforms averaged), I don't see any good way to recover your feature of interest.
EDIT - Having stated that a signal in noise needs a bandpass filter to recover it, I should explain that the multiplier used acts as what is called "mixer" in the RF world, and its effect is to frequency shift the signal by the reference frequency. This is useful in the case of the lock-in amplifier because it shifts the signal frequency to DC. In this case, a bandpass filter on the original signal becomes a low-pass filter on the processed signal, and the trick of the lock-in is that it's MUCH easier to make a very sharp, narrow lowpass filter than it is a very sharp, narrow bandpass filter. To begin with, the lowpass filter response is intrinsically referenced to DC, or zero Hz. This means that there is no central frequency of the filter to drift with time and temperature, which is a major problem with bandpass filters.
On the other hand, since the desired signal is now DC, you cannot recover the signal shape. Every deviation from the fundamental frequency (sine wave) shows up as a frequency deviation in the processed signal. If the artifact of interest is part of the signal at the base frequency, the frequency deviations show up as harmonics, and the closest to the fundamental is at twice the fundamental. This means that any close filtering will eliminate the part of the signal which corresponds to the glitch.
In the frequency domain, AM modulation ("multiplication") is the convolution of the modulated signal with the carrier. Thus, the Fourier Transform energy is at the carrier frequency, and has the bandwidth of the original signal.
So, as long as the width of your lock-in filter is as wide as the bandwidth of the original signal, you don't lose signal.
Best Answer
Trig theory to extract the phase & magnitude.
If you only care about phase...
Say you have a signal \$Acos(\omega t + \phi)\$ and you want to extract \$\phi\$. You can use an oscillator of the same frequency to extract this info BUT the issue is the phase.
\$V_{sig} = Acos(\omega t + \phi)\$
\$V_{osc} = cos(\omega t)\$
\$V_{sig}V_{osc} = Acos(\omega t + \phi) Cos(\omega t)\$
By the double angle identity:
\$ = \frac{1}{2}Acos(\phi) + \frac{1}{2}ACos(2 \omega t + \phi)\$
a DC term relating to the phase can be realised as well as a component at twice teh freqency of the carrier. The phase can then be extracted by a moving average filter at the carrier frequency of a simple low pass filter.
\$V_{sig}V_{osc}Filt = \frac{1}{2}Acos(\phi)\$
If you care about phase and magnitude
To clearly extract the phase and magnitude then two oscillators, in quadrature, are required.
\$V_{sig} = Asin(\omega t +\phi)\$
\$V_{osc0} = Xsin(\omega t)\$
\$V_{osc90} = Xcos(\omega t)\$
\$V_0 = Xsin(\omega t)Asin(\omega t +\phi) = \frac{XA}{2}(cos(\phi) - cos(2\omega t + \phi))\$ \$V_{90} = Xcos(\omega t)Asin(\omega t +\phi) = \frac{XA}{2}(sin(\phi) + sin(2\omega t + \phi))\$
Filter these signals to remove the twice carrier component
\$V_{0f} = \frac{XA}{2}(cos(\phi) ) \$
\$V_{90f} = \frac{XA}{2}(sin(\phi) ) \$
via trig:
\$\phi = atan( \frac{V_{90f}}{V_{0f}} )\$ \$A = \frac{2}{X}\sqrt{V_{0f}^2 + V_{90f}^2 } \$