Unexpected results of deconvolution with scipy.deconvolve

Question

Below I have plotted the signal (Lifetime decay) I am trying to deconvolve from a known impulse response function (IRF), as well as the IRF itself. I'm using scipy.signal.deconvolve.

Please note for my deconvolution code, I am using only the peak of the IRF, not the entire array sequence as shown on the plot above.

Here is the code I am trying to use:

IRF = IRF * (max(decay)/max(IRF))
# replace zeros to avoid error message 
IRF = np.where(IRF == 0, 0.1, IRF)    
decay = np.where(decay == 0, 0.1, decay)  
# take only the quotient part of the result 
deconv = scipy.signal.deconvolve(decay, IRF)[0]
# "padding" the deconvolved signal so it has the same size as the original signal 
s = int((len(decay)-(len(deconv)))/2)  ## difference on each side 
deconv_res = np.zeros(len(decay))   
end = int(len(decay)-s-1)  # final index
deconv_res[s:end] = deconv
deconv = deconv_res
# convolved normalized to decay height for plotting 
deconv_n = deconv * (max(decay)/max(deconv))

The IRF is an array of float64, the signal is an array of uint16.

I admit I'm not familiar with the maths of deconvolution in detail (I'm just writing a software for microscopy image analysis), so I am blindly changing different aspects of the code, like trying different divider functions, normalization, number of points of the IRF, etc, but nothing is producing anywhere near as expected. The last result I got looks like this (see plot of the original signal and what the signal it tried to deconvolve..)

Could anyone just point me in the right direction as to what could be the issue or what I can read to solve my problem? I see scipy uses inverse filter for the deconvolution, is it not the right method for time domain data? Could you recommend any other way to do this impulse response deconvolution in Python? (I couldn't find any other package for this) ..

Thanks in advance!