This improper integral comes from a problem of periodic orbit. The integral evaluates one half of the period.
In a special case, the integral is $$I=\int_{r_1}^{r_2}\frac{dr}{r\sqrt{\Phi^2(r,r_1)-1}}$$ where $$\Phi(u,v)=\frac{u\exp{(-u)}}{v\exp{(-v)}}$$
The interval follows $\Phi(r_1,r_2)=1$, $r_1<r_2$.
I have found a solution to a special case (by applying perturbation method to the original ODE), which is $$\lim_{r_1\rightarrow r_2} I =\pi$$ When $r_1 \rightarrow r_2$, we have $r_1, r_2 \rightarrow r_0$, where $r_0$ is the peak position of $g(r)=r\exp{(-r)}$.
The numerical verification is shown below:

$\uparrow$ The interval of the integral and the integrand
$\uparrow$ The integral as a function of $r_2$
My problem is to derive a closed form for $I(r_1)$, or even just a Taylor expansion about $r_0$. I appreciate any hint.
Thanks!
If you are interested, here is the general form of the integral: $$I=\int_{r_1}^{r_2}\frac{dr}{r\sqrt{\Phi^2(r,r_1)-1}}$$ where $$\Phi(u,v)=\frac{u\exp{(k(u))}}{v\exp{(k(v))}}$$ and $k$ is a decreasing function. The interval follows $\Phi(r_1,r_2)=1$, $r_1<r_2$.
By solving the original ODE using perturbation method, the solution to a special case is $$\lim_{r_1\rightarrow r_2} I =\frac{\pi}{\sqrt{1+r_0 k''(r_0)/k'(r_0)}}$$
When $k(r)=-r$, it reduces to $\pi$.
In fact, $\lim_{r_1 \rightarrow r_2} I (k(r)=-C\cdot r^n) = \pi/\sqrt{n}$.
Thanks to Fabian, the second derivative at $r=1$ matches $I=\pi-\frac{\pi}{12}\epsilon^2+O(\epsilon^3)$:

$\uparrow$ The above is the numerical second order derivative of figure 2.

I do not know how to obtain an explicit solution to the problem. However, it is possible to have a Taylor series of the integral $I(r_1)$ around $r_1=1$.
Let us first perform the substitution $$ r= r_1 (1-x) + r_2 x$$ such that the boundaries of the integral do not depend on $r_1$. In particular, we obtain the expression $$ I(r_1) =\int_0^1 \frac{r_2 -r_1}{(r_1 (1-x) + r_2 x) [\Phi(r_1 (1-x) + r_2 x,r_1)^2 -1 ]^{1/2}}\,dx\,.$$
Next, we need a relation between $r_2$ and $r_1$. If you look at the function $r\exp(-r)$ you see that it is monotonous on the interval $r\in[0,1]$ and $r\in[1,\infty]$. The inverse of this function is commonly called the Lambert W function. In particular, the inverse of the respective branches are denoted by $$ r=- W(-x) \in [0,1], \qquad r=-W_{-1}(-x) \in [1,\infty]\,.$$ With this notation, we have $$ r_2 =-W_{-1}(-r_1e^{-r_1}), \quad r_1 =-W(-r_2e^{-r_2})\,.$$
For $r_1$ close to $1$, we need the expansion of $W$ close to the branch point (see (4.26) of this paper). We obtain $$ r_2 = 1+ \epsilon + \frac{2}{3} \epsilon^2 + \frac{4}{9} \epsilon^3 + \frac{44}{135}\epsilon^4 + O(\epsilon^5) \tag{1}$$ with $\epsilon = 1-r_1$.
Investigating first the point $r_1=r_2=1$. We set $r_1 = 1-\epsilon$, $r_2 = 1+\epsilon$ (we know from (1) that $r_1$ and $r_2$ approach 1 from below and above at equal rate). To zeroth order in $\epsilon$, we obtain $$ I(1) =\int_0^1 \frac{1}{\sqrt{x(1-x)}}\,dx = \pi \,$$ as you have already observed.
In a next step, we look at $I(1) -I(r_1)$ for $r_1$ close to $1$. Using (1), we expand to third order in $\epsilon$. We obtain $$I(1)- I(r_1) = \int_0^1\left[\frac{\left(-30 x^2+34 x-5\right) \epsilon ^2}{9 \sqrt{(1-x) x}}+\frac{2 \left(472 x^3-858 x^2+422 x-33\right) \epsilon ^3}{135 \sqrt{(1-x) x}}+\frac{2 (2 x-1) \epsilon }{3 \sqrt{(1-x) x}}\right]dx= \frac{\pi}{12} \epsilon^2 + \frac{\pi}{18} \epsilon^3+O(\epsilon^4)\,.$$
To obtain a higher order approximation, we need more terms in (1). In particular, we have $$ r_2 = 1+\epsilon +\frac{2 \epsilon ^2}{3}+\frac{4 \epsilon ^3}{9}+\frac{44 \epsilon ^4}{135}+\frac{104 \epsilon ^5}{405}+\frac{40 \epsilon ^6}{189}+\frac{7648 \epsilon ^7}{42525}+\frac{2848 \epsilon ^8}{18225}+\frac{31712 \epsilon ^9}{229635}+\frac{23429344 \epsilon ^{10}}{189448875} +O(\epsilon^{11})\,.$$
Now, the expansion of the integral in $\epsilon$ yields $$ I(r_1) = \pi -\frac{\pi \epsilon ^2}{12}-\frac{\pi \epsilon ^3}{18}-\frac{23 \pi \epsilon ^4}{576}-\frac{67 \pi \epsilon ^5}{2160}-\frac{7613 \pi \epsilon ^6}{311040}-\frac{21419 \pi \epsilon ^7}{1088640}-\frac{320153 \pi \epsilon ^8}{19906560}-\frac{31342051 \pi \epsilon ^9}{2351462400} + O(\epsilon^{10})\,.$$