Plugging in the definition, $X_1$ following $\operatorname{Gamma}(a+b,1)$ means its density is
$$f_{X_1}(x_1) = \frac1{ \Gamma(a+b)}\, x_1^{a+b-1} e^{-x_1} \qquad \text{for}~~ 0 < x_1 < \infty $$
The density of $X_2$ is
$$f_{X_2}(x_2) = \frac{ \Gamma(a+b) }{ \Gamma(a) \Gamma(b) }\, x_2^{a-1} (1 - x_2)^{b-1} \qquad \text{for}~~ 0<x_2<1$$
The fact that $X_1 \perp X_2$ means their joint density is just the direct product
$$f_{X_1X_2}(x_1,\, x_2) = \frac{ x_1^{a+b-1} e^{-x_1} \, x_2^{a-1} (1 - x_2)^{b-1} }{ \Gamma(a) \Gamma(b) } \qquad \text{for}~~ \begin{cases}
0<x_1<\infty \\
0<x_2<1 \end{cases} \tag{1}\label{joint density}$$
The 2-dim transformation is
$$\begin{cases}
Y_1 = X_1 X_2 \\[1.5ex]
Y_2 = X_1 (1 - X_2)
\end{cases} \Longleftrightarrow \begin{cases}
X_1 = Y_1 + Y_2 \\[2ex]
X_2 = \dfrac{ Y_1 }{ Y_1 + Y_2}
\end{cases} \qquad \text{where}~~ \begin{cases}
0<y_1<\infty \\
0<y_2<\infty \end{cases}$$
with the Jacobian (of the inverse mapping) as
$$J = \left| \begin{matrix} \dfrac{ \partial x_1}{ \partial y_1} & \dfrac{ \partial x_1}{ \partial y_2} \\
\dfrac{ \partial x_2}{ \partial y_1} & \dfrac{ \partial x_2}{ \partial y_2}\end{matrix} \right| = \left| \begin{matrix} 1 & 1 \\
\dfrac{ y_2 }{ (y_1 +y_2)^2 } & \dfrac{ -y_1 }{ (y_1 +y_2)^2 } \end{matrix} \right| = \frac{-1}{ y_1 + y_2 }$$
The transformed joint density for $Y_1$ and $Y_2$ is
\begin{align}
f_{Y_1Y_2}( y_1 ,~y_2 ) &= |J| \cdot f_{X_1X_2}( x_1,\, x_2)\Bigg|_{\substack{x_1 = y_1+y_2 \\ x_2 = \frac{y_1}{y_1 + y_2}}} \qquad \text{, plug in Eq.(\ref{joint density})}\\
&= \frac1{ y_1 + y_2} \cdot \frac{ (y_1 + y_2)^{a+b-1} e^{-(y_1 + y_2)} }{ \Gamma(a) \Gamma(b) }\, \left(\frac{y_1}{ y_1 + y_2}\right)^{a-1} \left(\frac{y_2}{ y_1 + y_2} \right)^{b-1} \\
&= \frac1{ \Gamma(a) \Gamma(b) }\, y_1^{a-1} y_2^{b-1} e^{-(y_1 + y_2)} \qquad \text{for}~~0<y_1<\infty ,~0<y_2<\infty
\end{align}
The marginal density of $Y_1$ can be obtained from the joint as
\begin{align}
f_{Y_1}(y_1) &= \int_{y_2 = 0}^{\infty} f_{Y_1Y_2}( y_1 ,~y_2 ) \,\mathrm{d}y_2 \\
&= \frac1{\Gamma(a)} y_1^{a-1} e^{-y_1} \int_{y_2 = 0}^{\infty} \frac1{\Gamma(b)} y_2^{b-1} e^{-y_2} \,\mathrm{d}y_2 \qquad \small\begin{aligned}[c]
&\text{integral is just the} \\
&\text{kernel of Gamma distribution}\end{aligned} \\
&= \frac1{\Gamma(a)} y_1^{a-1} e^{-y_1}
\end{align}
Thus one identifies the distribution of $Y_1$ as $\operatorname{Gamma}(a,1)$.
Similarly, or noting the symmetry in the joint $f_{Y_1Y_2}( y_1 ,~y_2 )$, we have $Y_2$ follows $\operatorname{Gamma}(b,1)$.
First I would recommend you to move the indicator function into the integral by splitting the integral into seperate cases:
$$\int_0^\infty \lambda e^{-\lambda y} \cdot \underbrace{\mathbf{1}_{[0,1]} (z-y)}_{= \mathbf{1}_{[z, z-1]}(y)} \, dy = \mathbf{1}_{[0, 1]}(z) \cdot \int_0^z \lambda e^{-\lambda y} \, dy \; + \mathbf{1}_{(1, \infty)}(z) \int_{z-1}^z \lambda e^{-\lambda y} \, dy$$
Now continue by integrating the terms.
Edit (Explanation)
Basically we are splitting up the integral in three different cases as our indicatorfunction $\mathbf{1}_{[z-1, z]}(y)$ behaves different for different $z$.
So technically
$$\int_0^\infty \lambda e^{-\lambda y} \cdot \mathbf{1}_{[z, z-1]}(y) \, dy = \underbrace{\mathbf{1}_{(-\infty, 0]}(z) \cdot \int_0^\infty \lambda e^{-\lambda y} \cdot \mathbf{1}_{[z, z-1]}(y) \, dy}_{= 0} \\
+ \mathbf{1}_{(0, 1)}(z) \cdot \int_0^\infty \lambda e^{-\lambda y} \cdot \mathbf{1}_{[z, z-1]}(y) \, dy +\mathbf{1}_{[1, \infty)}(z) \int_0^\infty \lambda e^{-\lambda y} \cdot \mathbf{1}_{[z, z-1]}(y) \, dy$$
The first integral is $0$ because $y$ has to be $\ge 0$ and $\in [z-1, z]$ but $z < 0$.
The choice of the splitting borders is useful because for $z \in [0, 1]$ the lower bound of the integral has to respect both the $0$ from $\mathbf{1}_{[0, \infty)}(y)$ but also the $z-1$ from $\mathbf{1}_{[z-1, z]}(y)$ so we have to take the maximum of both, which is exactly $0$ if $z \in [0,1]$. For $z \in [1, \infty)$ follows $z-1 > 0$ so the maximum of the lower bound is $z-1$ not $0$.
Best Answer
Use the following result:
Assuming $Y$ and $Z$ are independent, the PDF of $X = YZ$ is given by:
$$f_X(x) = \int_{-\infty}^{\infty} \frac{1}{|u|} f_{Y}(u) f_Z\left(\frac{x}{u}\right) du$$