The convolution you want to do is possible in general. The problem in your code seems to me that you defined the Gaussian for x=pi0_ErrM and the resolution model and convolution for x=pi0_M_prefit.
If you want to a convolution, you got to define the inputs for the same x axis. Maybe you meant to use pi0_ErrM for everything up to the convolution? And then pi0_M_prefit only for the RooJohnson in the other dimension?
Let me know if this goes in the right direction and if you have further questions!