The formula $$r(z) \sim \int_0^{N_S(z)} N_T(z-N_S^{-1}(x)) dx$$ can be re-written in a more appealing way. With the change of variable $u=N_S^{-1}(x)$ it becomes $$r(z) \sim \int_0^{z} N_T(z-u)N'_S(u) du,$$ where $N'_S(u)$ is the derivative of $N_S(u)$ with respect to $u$. With an additional change of variable $u=zv$ it becomes $$r(z) \sim z\int_0^{1} N_T(z(1-v))N'_S(zv) dv.$$ Likewise $$t(z) \sim r'(z) = \frac{dr(z)}{dz} =z\int_0^{1} N'_T(z(1-v))N'_S(zv) dv .$$
An interesting case is when $S=T$ and $$N_S(u) \sim \frac{a u^b}{(\log u)^c}, \mbox{ with } 0<a, 0<b<1, \mbox{ and } c>0.$$ This covers sums of two primes ($a=1, b=1, c=1$) and sums of two squares ($a=1, b=\frac{1}{2}, c=0$). We have: $$r(z) \sim \frac{a^2b z^{2b}}{(\log z)^{2c}}\cdot \int_0^1 (1-v)^b v^{b-1}dv = \frac{a^2b z^{2b}}{(\log z)^{2c}}\cdot \frac{\Gamma(b)\Gamma(b+1)}{\Gamma(2b+1)}$$
$$r'(z) \sim \frac{2 a^2 b^2 z^{2b-1}}{(\log z)^{2c}}\cdot \int_0^1 (1-v)^b v^{b-1}dv = \frac{2a^2 b^2 z^{2b-1}}{(\log z)^{2c}}\cdot \frac{\Gamma(b)\Gamma(b+1)}{\Gamma(2b+1)}$$
Notes
Solutions such as $z=x+y$ and $z=y+x$ count as two solutions: $(x,y)$ and $(y, x)$.
The asymptotic formula for $t(z) \sim r'(z)$, representing the number of solutions to $z=x+y$ with $x\in S, y\in T$ is true only on average, as $z$ becomes larger and larger. There may still be infinitely many integer $z$'s for which $t(z)=0$ even if $r'(z) \rightarrow\infty$ as $z\rightarrow\infty$.
We assume that the sets $S$ and $T$ are "well balanced", both for small and large values. For instance, if you remove the first $10^{5000}$ elements of $S$, the asymptotic formula for $N_S(u)$ remains unchanged, but this is likely to cause many formulas to fail.
On some tests, I noticed that there are more solutions (on average) to $z=x+y$ with $x\in S, y\in T$ (here $x, y, z$ are integers), if $z$ is even.
If $S=T$ is the set of primes, some adjustments must be made because the primes are not "well balanced", they are less random than they seem (for instance the sum of two odd primes can not be an odd number, but there are also more subtle issues). This is best described in the Wikipedia entry about Goldbach's conjecture (see section about heuristics).
To generate a set like $S$, one way is as follows. Use a random number generator function $U$ returning independent uniform deviates on $[0, 1]$. If $U(k) < N'_S(k)$ then add the integer $k$ to the set $S$, otherwise discard it. Do that for all integers.
For sums involving three terms, say $R+S+T$, you can proceed as follows: first work on $S'=R+S$ and derive all the asymptotics for $S'$ using the methodology proposed here. Then work on $S'+T$.
If there are singularities in the functions $N_S$ or $N_S'$, they must be handled properly in the integral formulas, unless the integrals are improper but converging.
Generalization of the formula
It also works if $S\neq T$. Say
$$N_S(u) \sim \frac{a_1 u^{b_1}}{(\log u)^{c_1}}, N_T(u) \sim \frac{a_2 u^{b_2}}{(\log u)^{c_2}}$$ with $0<a_1,a_2, 0<b_1, b_2 <1$, and $c_1, c_2>0$. I am working on getting the formula right, I will publish here when double-checked.
In particular, it applies to sums of a square and a prime, see here.