Is $\pi$ the best constant in this inequality?

To answer the question in the title, no, $\pi$ is not the best constant in the inequality. We can approximate the optimal $M$ from below simply by choosing finite sets of completely monotonous functions and finding the smallest value of $M$ which works for convex combinations of these. For example, considering $f(x)=(1+x)^{-2}$ shows that $M\ge3$. Choosing $f(x)=(1+x)^{-5/2}$ improves this to $M\ge 3.0792$. As also stated in the comments, choosing convex combinations of functions of the form $(1+x)^{-\alpha}$ gives $M\ge 3.12935$. However, we do have $M < \pi$. In fact, as I'll show, it is possible to transform the dual problem which approximates $M$ from above, and the fact that it is (very slightly) less than $\pi$ will be immediate from this.

The class, $\mathcal{F}$, of totally monotonous functions from $(0,\infty)$ to $\mathbb{R}$ is closed under scaling the domain. That is, if $f\in\mathcal{F}$ then $x\mapsto f(ax)$ is also in $\mathcal{F}$ for each fixed $a > 0$. This means that, for any $M > 0$, the following inequalities over $f\in\mathcal{F}$ are equivalent. \begin{align} &\left(\int_0^\infty f(x)dx\right)^2\le M\sqrt{\int_0^\infty f^2(x)dx\int_0^\infty x^2f^2(x)dx}&&{\rm(1)}\\ &\left(\int_0^\infty f(x)dx\right)^2\le\frac12M\int_0^\infty(1+x^2)f^2(x)dx&&{\rm(2)} \end{align} We can assume that $\int_0^\infty(1+x^2)f^2(x)dx$ is finite, otherwise the inequalities above are both trivial as the right hand sides are infinite. Suppose that (1) holds. Then, the AM-GM inequality $\sqrt{PQ}\le\frac12(P+Q)$ implies that (2) holds. Conversely, suppose that (2) holds. Applying this to $x\mapsto f(x/a)$, changing variables (replace $x$ by $ax$ in the integral), and dividing through by $a^2$ gives $$ \left(\int_0^\infty f(x)dx\right)^2\le\frac12M\int_0^\infty(a^{-1}+ax^2)f^2(x)dx. $$ Taking $a=\sqrt{\int f^2dx/\int x^2f^2dx}$ gives (1). So (1) and (2) are equivalent. In what follows, I'll consider the inequality in the form (2).

Now, the space $\mathcal{H}$ of measurable functions $f\colon(0,\infty)\to\mathbb{R}$ such that $\int(1+x^2)f^2dx < \infty$ is a Hilbert space (after identifying functions which are equal almost everywhere), with inner product defined by $$ \left\langle f,g\right\rangle = \frac\pi2\int_0^\infty(1+x^2)f(x)g(x)dx. $$ Let us now define $\varphi\in\mathcal{H}$ by $\varphi(x)=\frac2\pi(1+x^2)^{-1}$. This is normalised so that $\lVert\varphi\rVert=1$, and $\int f(x)dx=\langle f,\varphi\rangle$ for all $f\in\mathcal{H}$. Inequality (2) is equivalent to $\langle f,\varphi\rangle^2\le\frac M\pi\lVert f\rVert^2$.

Using this notation, we can write $$ \langle f,\varphi\rangle^2 =\lVert f\rVert^2 - \lVert f - \langle f,\varphi\rangle\varphi\rVert^2. $$ Dividing through by $\lVert f\rVert^2$ and optimizing over $f$ gives the optimal value of $M$, \begin{align} \frac M\pi &=\sup_{f\in\mathcal{F}\setminus\{0\}}\lVert f\rVert^{-2}\langle f,\varphi\rangle^2\\ &=\sup_{f\in\mathcal{F},\lVert f\rVert=1}\langle f,\varphi\rangle^2\\ &=1-\inf_{f\in\mathcal{F},\lVert f\rVert=1}\lVert f-\langle f,\varphi\rangle\varphi\rVert^2&&{\rm(3)}\\ &=1-\inf_{f\in\mathcal{F},c\in\mathbb{R},\lVert f\rVert=1}\lVert f-c\varphi\rVert^2&&{\rm(4)} \end{align} The second equality above uses the fact that $\mathcal{F}$ is closed under scaling the range. That is, if $f\in\mathcal{F}$ then so is $x\mapsto af(x)$ for any fixed $a > 0$. Taking $a=1/\lVert f\rVert$ allows us to restrict the optimisation to functions $f\in\mathcal{F}$ normalised so that $\lVert f\rVert=1$. The final term on the right hand side of this sequence of equalities, $\lVert f-c\varphi\rVert^2$, is a quadratic in $c$ minimised at $c=\langle f,\varphi\rangle$, which gives the final inequality.

Now, it can be seen that the class $\mathcal{F}$ of totally monotonous functions is a closed convex subset of $\mathcal{H}$. This means that the right hand sides of (3) and (4) attain their bounds for some $f\in\mathcal{F}$ and $c\in\mathbb{R}$. At this bound, we cannot have $f=c\varphi$ as neither $\varphi$ nor $-\varphi$ are totally monotonous. So, from (4), we have $M < \pi$ as claimed.

I'll now show how we can be more precise about approximating $M$ from above, and a dual problem can be used find arbitrarily close upper bounds for $M$. First, $f$ is totally monotonous (almost everywhere) iff $\int \theta^{(r)}fdx\ge0$ for all $r\in\mathbb{Z}_{\ge0}$ and smooth non-negative functions $\theta$ with compact support in $(0,\infty)$. Taking $\chi=\varphi\theta^{(r)}$, this is equivalent to $\langle\chi,f\rangle\ge0$. So, $\mathcal{F}$ is the intersection of sets of the form $\{f\colon\langle\chi,f\rangle\ge0\}$, so is closed and convex in $\mathcal{H}$. Next, the second derivative of $\varphi(x)$ is $\frac4\pi(3x^2-1)(1+x^2)^{-3}$, so it is concave on $(0,1/\sqrt{3})$. Therefore, if $\theta$ is a smooth nonnegative function with compact support in $(0,1/\sqrt{3})$ then $\int\theta^{(2)}\varphi dx < 0$. Setting $\chi=\varphi\theta^{(2)}$ we have $\langle\chi,\varphi\rangle < 0$ and $\langle\chi,f\rangle\ge0$ for all $f\in\mathcal{F}$. As we'll see, this implies that $M$ is bounded above by $\pi(1-\langle\xi,f\rangle^2/\lVert\xi\rVert^2) < \pi$.

Moving on to the dual problem, let $\mathcal{F}^*$ be the closure in $\mathcal{H}$ of positive linear combinations of functions of the form $\varphi\theta^{(r)}$ where $r\in\mathbb{Z}_{\ge0}$ and $\theta$ is a smooth nonnegative function with compact support in $(0,\infty)$. We have a kind of duality between $\mathcal{F}$ and $\mathcal{F}^*$. \begin{align} \mathcal{F}&=\left\{f\in\mathcal{H}\colon\langle\chi,f\rangle\ge0\ \forall\chi\in\mathcal{F}^*\right\},\\ \mathcal{F}^*&=\left\{\chi\in\mathcal{H}\colon\langle\chi,f\rangle\ge0\ \forall f\in\mathcal{F}\right\}. \end{align} The first of these equalities follows from from the definition of totally monotonous functions, as explained above. The second follows from the first, by the Hahn-Banach theorem. I'll base the dual problem on the following lemma.

Lemma: Suppose that $f\in\mathcal{F}$ and $\chi\in\mathcal{F}^*$ satisfy $\lVert f\rVert=\lVert\chi\rVert=1$ and $\langle f,\varphi\rangle\ge0\ge \langle \chi,\varphi\rangle$. Then, \begin{align} \langle f,\varphi\rangle^2+\langle \chi,\varphi\rangle^2\le1.&&{\rm(5)} \end{align} Furthermore, this is an equality if and only if, subject to the conditions above, $f$ maximises $\langle f,\varphi\rangle$ and $\chi$ minimises $\langle \chi,\varphi\rangle$.

Proof: For any pair $u_1,u_2$ of orthonormal elements of $\mathcal{H}$, \begin{align} 1=\lVert \varphi\rVert^2&=\langle u_1,\varphi\rangle^2+\langle u_2,\varphi\rangle^2+\lVert\varphi-\langle u_1,\varphi\rangle u_1-\langle u_2,\varphi\rangle u_2\rVert^2\\ &\ge\langle u_1,\varphi\rangle^2+\langle u_2,\varphi\rangle^2. \end{align} Apply this with $u_1=(f-\langle\chi,f\rangle\chi)/\lVert f-\langle\chi,f\rangle\chi\rVert$ and $u_2=\chi$, and multiply through by $\lVert f-\langle\chi,f\rangle\chi\rVert^2$, $$ \langle f-\langle\chi,f\rangle\chi,\varphi\rangle^2+\langle\chi,\varphi\rangle^2\lVert f-\langle\chi,f\rangle\chi\rVert^2\le\lVert f-\langle\chi,f\rangle\chi\rVert^2 $$ Expanding out the inner products, using $\lVert f\rVert=\lVert\chi\rVert=1$, and rearringing the terms gives, $$ \langle f,\varphi\rangle^2+\langle\chi,\varphi\rangle^2\le1-\langle\chi,f\rangle^2+2\langle\chi,f\rangle\langle f,\varphi\rangle\langle\chi,\varphi\rangle. $$ Applying the conditions that $\langle\chi,f\rangle$ and $\langle f,\varphi\rangle$ are positive and $\langle\chi,\varphi\rangle$ is negative shows that the right hand side is bounded by $1$, giving inequality (5).

Equality is only possible in (5) when both terms on the left hand side are maximised or, equivalently, when $\langle f,\varphi\rangle$ is maximised and $\langle\chi,\varphi\rangle$ is minimised. So, it only remains to be shown that equality can be achieved.

Choose $f\in\mathcal{F}$ with $\lVert f\rVert=1$ to maximise $\langle f,\varphi\rangle$. Then, for any $g\in\mathcal{F}^*$ and $\epsilon > 0$ consider $h\in\mathcal{F}$ defined by, \begin{align} h&\equiv\left(f+\epsilon g\right)/\lVert f+\epsilon g\rVert= \left(f+\epsilon g\right)\left( 1+2\epsilon \langle f,g\rangle+\epsilon^2\right)^{-1/2}\\ &=f+\left(g-\langle f,g\rangle f\right)\epsilon+(\epsilon^2). \end{align} The optimality of $f$ gives $\langle f-h,\varphi\rangle\ge0$ and, in order for this to hold for arbitrarily small $\epsilon$, $$ \langle\langle f,\varphi\rangle f-\varphi,g\rangle =\langle -g+\langle f,g\rangle f,\varphi\rangle\ge0. $$ It follows that $\langle f,\varphi\rangle f-\varphi$ is in $\mathcal{F}^*$, and its inner product with $\varphi$ is $\langle f,\varphi\rangle^2-1$, which is negative. So, normalising this, $\chi=(\langle f,\varphi\rangle f-\varphi)/\lVert\langle f,\varphi\rangle f-\varphi\rVert$ satisfies the conditions of the lemma, and it can be checked that it gives equality in (5). QED

Now, it follows immediately from this lemma that if we set \begin{align} \alpha &= \sup_{f\in\mathcal{F},\lVert f\rVert=1}\langle f,\varphi\rangle > 0\\ \beta &= \inf_{\chi\in\mathcal{F}^*,\lVert\chi\rVert=1}\langle\chi,\varphi\rangle < 0 \end{align} then $\alpha^2+\beta^2=1$ and, $$ M = \pi\alpha^2 = \pi(1-\beta^2) < \pi. $$ So, we have the dual problem of trying to find $\beta$ by searching over $\chi\in\mathcal{F}^*$, which approximates $M$ from above.