How do I normalize a predesigned IIR filter? Obviously  my question here is very similar to  this previously asked question: How to normalize predesigned IIR filter coefficients to have unity gain?

Unfortunately, in the responses available to that question, it isn't super clear how the answers are obtained.
I have a pre-designed IIR filter that is applied in 2 steps:

Step 1:

    Yn =  Xn + a1*Yn-1 - b1*Yn-2


Step 2:

    Zn =  Yn + a2*Zn-1 - b2*Zn-2


with Xn being my initial input and Zn my final output.
Given these 2 pairs of coefficients (a1,b1) and (a2,b2), how do I go about finding the gain of this IIR filter to ensure it has unity gain.

Thanks for the help
As mentioned in the comments, you have to choose a frequency $\omega_0$ at which you want to normalize the gain. This could be DC (i.e., $\omega_0=0$) or any other frequency, depending on the filter's characteristic.

You have two filters in series, each with transfer function

$$H_i(z)=\frac{1}{1-a_iz^{-1}+b_iz^{-2}},\qquad i=1,2\tag{1}$$

The total transfer function is just the product of the two individual transfer functions:

$$H(z)=H_1(z)\cdot H_2(z)\tag{2}$$

The gain at frequency $\omega_0$ equals the magnitude of the total transfer function $(2)$ evaluated at $z=e^{j\omega_0}$. So if you scale one of the two transfer functions $H_i(z)$ by the inverse of that gain, your total system will have unity gain at $\omega_0$.

Scaling in your time domain equations simply means scaling the input variable.

Note that $\omega_0$ is a normalized frequency in radians:

$$\omega_0=2\pi\frac{f_0}{f_s}\tag{3}$$

where $f_s$ is the sampling frequency.