If $X\sim \Gamma(a,\sigma_x^2)$ and $Y\sim \Gamma(b,\sigma_y^2)$. What will be the probability density function of R? Where $R=\frac{X+C}{X+Y}$, here $C$ is a positive constant, $\Gamma(.,.)$ denotes standard gamma probability density function and '$\sim$' represents 'distributed as'. X and Y are independent random variables.

First of all, since $R=(X+C)/(X+Y)$ (and $X$ and $Y$ are independent gamma variables), then the valid range of $R$ is a priori $(0,\infty)$. The density of $(R,S)$ is given by $f_{R,S} (r,s) = f_X (sr  C)f_Y (s  sr + C)s$, where $f_X$ and $f_Y$ are the densities of $X$ and $Y$ (see the remark below). This leads to $srC > 0$ and $ssr+C > 0$. Thus, $s > C/r$ and $s(r1) < C$. Hence, if $r > 1$, then $C/r < s < C/(r1)$, while if $0 < r < 1$, then $C/r < s < \infty$. So, for $0 < r < 1$ you would use $f_R (r) = \int_{C/r}^\infty {f_{R,S} (r,s)\,{\rm d}s}$, and for $r > 1$, $f_R (r) = \int_{C/r}^{C/(r1)} {f_{R,S} (r,s)\,{\rm d}s}$. In general, a nice expression for $f_R (r)$ is not likely to be found. [Minor point: in view of the title, note that $X+Y$ is, in general, not gamma distributed.] Remark: As Didier observed, a factor of $s$ was missing in the original expression for $f_{R,S}(r,s)$ (now fixed). The density $f_{R,S}(r,s)$ is found as follows. $(X,Y)$ has density $f_{X,Y}(x,y)=f_X (x) f_Y (y)$. Noting that $X=SRC$ and $Y=SSR+C$, it follows by the standard formula for transformation of rv's that $f_{R,S} (r,s)$ is given by $f_{X,Y}(srC,ssr+C)$ times $J(r,s)$, where $J(r,s)$ is given by the determinant $\frac{{\partial (sr  C)}}{{\partial r}}\frac{{\partial (s  sr + C)}}{{\partial s}}  \frac{{\partial (sr  C)}}{{\partial s}}\frac{{\partial (s  sr + C)}}{{\partial r}} = s$. Thus, $f_{R,S} (r,s)$ is $s$ times the original expression. 


The OP asked for the probability density function $f_R$ of $R=(X+C)/(X+Y)$ for a positive $C$. Shai recalled the classical strategy to compute $f_R$ and explained why a nice expression of $f_R$ is unlikely to emerge in the general case. The limit case $C=0$ is nicer. Then $R=X/(X+Y)$, thus $R\in(0,1)$ almost surely and $$ f_R(r)=\frac{cr^{a1}(1r)^{b1}}{\left(\sigma_y^2r+\sigma_y^2(1r)\right)^{a+b}}, $$ for a suitable normalizing constant $c$ depending on $(a,b,\sigma_x^2,\sigma_y^2)$. As expected, if $\sigma_x^2=\sigma_y^2$, $f_R$ is Beta$(a,b)$ density. More generally, $R=T(R_0)$ where $R_0$ is Beta$(a,b)$ and $$ T(r)=\frac{\sigma_x^2r}{\sigma_x^2r+\sigma_y^2(1r)}. $$ 

