The scattering amplitude can be obtained in QFT from the $n$-point Green's function by the LSZ reduction formula. The $n$-point Green's function are a correlation function of a string of fields :
$$
G(q_1,..,q_n)=\int dx_1...dx_n\,e^{-iq_1x_1}...e^{-iq_nx_n}\langle 0|\mathcal{T}\{A_1(x_1)...A_n(x_n)\}|0\rangle
$$
Now, we can factorize the time-ordering using the theta function, giving terms like:
$$
\int dx_1...dx_n\,e^{-iq_1x_1}...e^{-iq_nx_n}\langle 0|\mathcal{T}\{A_1(x_1)...A_r(x_r)\}\mathcal{T}\{A_{r+1}(x_{r+1})...A_n(x_n)\}|0\rangle \times
$$
$$
\times\theta(\min[x_1^0...x_r^0]-\max[x_{r+1}^0...x_n^0])
$$
The idea is that very old idea in quantum mechanics when we insert the identity operator $I$ in a given representation:
$$
I=\sum_{p,\sigma}|p,\,\sigma\rangle\langle p,\,\sigma|+...
$$
where the first term is one-particle states projector and $...$ are multi-particle states projectors. Keeping just the one particle terms, we have:
$$
\int dx_1...dx_n\,e^{-iq_1x_1}...e^{-iq_nx_n}\int d^3p\sum_\sigma\langle 0|\mathcal{T}\{A_1(x_1)...A_r(x_r)\}|p,\,\sigma\rangle\times
$$
$$
\times\langle p,\,\sigma|\mathcal{T}\{A_{r+1}(x_{r+1})...A_n(x_n)\}|0\rangle \,\theta(\min[x_1^0...x_r^0]-\max[x_{r+1}^0...x_n^0])
$$
we get a factorization of the the correlation function into two pieces, connected just by the theta function. To get ride of this we can use translation symmetry to make:
$$
\langle 0|\mathcal{T}\{A_1(x_1)...A_r(x_r)\}|p,\,\sigma\rangle=e^{ip.x_1}\langle 0|\mathcal{T}\{A_1(0)...A_r(y_r)\}|p,\,\sigma\rangle
$$
$$
\langle p,\,\sigma|\mathcal{T}\{A_{r+1}(x_{r+1})...A_n(x_n)\}|0\rangle=e^{-ip.x_{r+1}}\langle p,\,\sigma|\mathcal{T}\{A_{r+1}(0)...A_n(y_{n})\}|0\rangle
$$
and under this new variables the theta function becomes:
$$
\theta(x_1^0-x_{r+1}^0+\min[0...y_r^0]-\max[0...y_n^0])
$$
using the Fourier representation:
$$
\theta (\tau)=-\frac{1}{2\pi i}\int_{-\infty}^{+\infty}\frac{d\omega e^{-i\omega \tau}}{\omega + i\varepsilon}
$$
now we can perform the integration over $x_1$, $x_{r+1}$ and $p$. Some delta Dirac will shows up enforcing conservation of the momentum between the two blobs and an extra delta enforcing $\omega$ being equal to the energy transferred between the blobs minus the energy $E_p$ of the one-particle state. Then, the pole $(\omega +i\varepsilon)^{-1}$ that comes from the theta function will give rise to a pole $(q^{0}-E_p+i\varepsilon)^{-1}$ where $q^0$ is the energy transferred between the blobs.
Around the pole, we can make $(q^{0}-E_p+i\varepsilon)^{-1}\rightarrow 2E_p (q^2+m^2-i\varepsilon)^{-1}$, with $\vec{q}=\vec{p}$. The term $2E_P$ is absorbed by the integrals to form a relativistic invariant measure. This is how the pole $(q^2+m^2-i\varepsilon)^{-1}$ show up.
Now let us look at the residue. After the LSZ reduction formula, the residue will be precisely the product of two new amplitudes:
$$
\lim_{q^2\rightarrow-m^2}(q^2+m^2-i\varepsilon)A(q_1,...,q_n)=A(q,q_2,...,q_r)\times A(q,q_{r+2},...,q_n)
$$
where $q=q_1+...+q_r=-(q_{r+1}+...+q_{n})$. This means that we have a pole whatever the amplitudes $A(q,q_2,...,q_r,p)$ and $A(q,q_{r+2},...,q_n)$ are non-zero.
For more detailed explanation and calculation see Weinberg, QFT, volume 1, chapter 10.