Should removing DC component (substracting average) be after applying window?

I noticed that in the current code (https://github.com/endolith/waveform_analysis/blob/master/waveform_analysis/thd.py#L61) the DC component is removed before applying window. The calculated DC offset would be different after the window has been applied, causing the resulting FFT to have a rather large DC component in some cases.

For example, I have a 1kHz -3dBFS sine tone 32-bit wave file at 48kHz sample rate. I feed only the first 16K samples to the code (obviously the last sine cycle won't be complete, and that's why window functions are used). The result I got is a THD+N of -77.8dB while I expect it to be much better. If I increase the FFT points to 1M point, the THD+N gets improved to <-100dB. By inspecting the noise spectrum in the 16K point case, I noticed that the DC bin has a -77.8dB magnitude, explaining the rather poor THD+N. If I move the subtraction to after the windowing (windowed -= mean(windowed)), the THD+N immediately gets improved to -189.7dB, and the performance is consistent across different FFT sizes between 16K and 2M. (while in the previous cases, performances improve with larger FFT size).

These make sense to me however I am not sure if my modifications are correct or not. Do you have any idea about this?