Probability Random Variable Greater than Maximum Order Statistic of Correlated RV's

Let $\epsilon_{0}$, $\epsilon_{1}$ and $\epsilon_{2}$ be independent standard normal random variables.

I would like to compute $\mathbb{P}[b(p \epsilon_{0} + (1-p)\epsilon_{1}) > \max\{\epsilon_{0}p + (1-p)\epsilon_{2}, C + a(p \epsilon_{0} + (1-p)\epsilon_{1}))]$,

For $C,b>0$, $p \in (0,1)$, $a \in (0,1)$.

How do I write this problem down in integral form so that I can plug it into a numerical simulator? I'm having trouble because of the $\epsilon_{0}$...

-

You want to integrate the joint density $f(x,y,z)$ of $\epsilon_0, \epsilon_1, \epsilon_2$ over the intersection of the regions $R_1$ where $b (p x + (1-p) y) > p x + (1-p) z$ and $R_2$ where $b (p x + (1-p) y) > C + a (p x + (1-p) y)$. You can rewrite $R_1$ as $z < \dfrac{(b-1) p x}{1-p} + b y$, and $R_2$ as $(b-a)(p x + (1-p)y) > C$; if $b>a$ this says $x > \dfrac{C}{p(b-a)} - \dfrac{1-p}{p} y$ (change $>$ to $<$ if $b<a$). So if $b>a$ the integration can be set up as $$\int_{-\infty}^\infty dy \int_{{C}/(p(b-a)) - (1-p)y/{p}}^\infty dx \int_{-\infty}^{(b-1)px/(1-p)+by} dz\ f(x,y,z)$$