Understanding truncation and rounding error in IEEE floating point system?
As listed your reasoning is indeed circular: you start by assuming $h \approx \sqrt{\epsilon}$, then use the results in point 3 to infer that we should take $h \approx \sqrt{\epsilon}$.  So that's not convincing.  So let's redo the analysis more carefully, then we'll look at the meaning of the constant.

A more careful analysis

As in your analysis, we'll assume $f'(x) \approx 1$ and $f''(x) \approx 1$.

Your analysis of truncation error is fine; the truncation error is indeed $\frac{h}{2} f''(x) + O(h^2)$, which given your assumptions is $\frac{h}{2} + O(h^2)$.  In our analysis $h$ will be small enough that $O(h^2)$is negligible, so we'll call the truncation error $\frac{h}{2}$ for simplicity.

What about rounding error?  The rounding error of representing $f(h)$ in floating point is $\epsilon$ in relative error, where $\epsilon$ is the machine epsilon; in absolute error, that is $\epsilon f(h)$.  Since you assumed $f(h) \approx 1$, that is $\approx \epsilon$.  Similarly, the rounding error of representing $f(x+h)$ is $\epsilon f(x+h)$ in absolute error; since $f(x+h) \approx f(x) \approx 1$, that is $\approx \epsilon$ as well.  Now the absolute error of $f(x+h)-f(x)$ is at most about $2\epsilon + \epsilon (f(x+h)-f(x))$; the $2\epsilon$ comes from the error in each of the two subterms (and the triangle inequality), and the $\epsilon (f(x+h)-f(x))$ comes from roundoff error in representing the result.  Since $f(x+h)-f(x) \approx h f'(x) \approx h$, that means the absolute error of $f(x+h)-f(x)$ is at most about $2\epsilon + \epsilon h$.

Then we divide by $h$.  The absolute error becomes $2\epsilon/h + \epsilon$, plus a term due to representing the result in floating point; that extra term is $\epsilon ((f(x+h)-f(x)/h)$, and since $(f(x+h)-f(x)/h \approx f'(x) \approx 1$, we see that the extra term is $\approx \epsilon$.  Therefore the absolute error in computing $(f(x+h)-f(x)/h$ is about
$2\epsilon/h + \epsilon + \epsilon$.  So, the total roundoff error is about $2\epsilon/h + 2\epsilon$, in absolute error.

The total error is at most the sum of the truncation error and the roundoff error, or about
$$\frac{h}{2} + {2\epsilon \over h} + 2\epsilon.$$
We want to know what value of $h$ minimizes this expression.  We can figure that out by letting $g(h) = \frac{h}{2} + 2\epsilon/h + 2\epsilon$ and minimizing $g(h)$; the minima must be at a point where the derivative is zero.  So, differentiate:
$$g'(h) = \frac{1}{2} - {2\epsilon \over h^2}.$$
Setting $g'(h)=0$ and solving for $h$ yields $h^2 =4 \epsilon$, i.e., the function $g(h)$ is minimized at $h=2\sqrt{\epsilon}$; at that point, we have a total (absolute) error of at most $g(2\sqrt{\epsilon}) \approx 2 \sqrt{\epsilon} + 2 \epsilon \approx 2\sqrt{\epsilon}$.

The summary: assuming $f(x)\approx 1$ and $f'(x) \approx 1$, the total error is minimized when you choose $h$ to be about $h\approx 2\sqrt{\epsilon}$.  If you choose $h$ too much larger than that, the truncation error increases and starts dominating.  If you choose $h$ too much smaller than that, the roundoff error increases and starts dominating.

So, if you know $f(x)\approx 1$ and $f'(x) \approx 1$, then it's reasonable to choose $h\approx 2\sqrt{\epsilon}$.

The constant

Now, you might be wondering a bunch of things.  Why the constant 2?  Where did that come from?  How critical is it for the constant to be 2?  What if I use $h \approx \sqrt{\epsilon}$ instead of $h \approx \sqrt{2\epsilon}$?

Well, if you redo the calculation with a different value of $h$, you'll see that the total error is not too sensitive to the constant 2.  For instance, if you choose $h=\epsilon$, then the total error will be at most about $g(\epsilon) \approx  2.5 \sqrt{\epsilon} + 2 \epsilon \approx 2.5 \sqrt{\epsilon}$, which is not that different.  So don't worry too much about the exact value of the constant.

Where did the constant 2 came from?  Well, partly it is coming from the fact that we are doing multiple operations, and you can get a little bit of roundoff error from each operation, so the error can accumulate.  And, partly it is coming from the assumptions that $f'(x) \approx 1$ and $f''(x) \approx 1$.  Of course, those assumptions are pretty sketchy.  In practice, we don't know what $f'(x)$ is before doing the computation (that's the whole point of this computation; we are trying to estimate/compute $f'(x)$), and we don't know what $f''(x)$ is.  So, that assumption also looks circular.

But here's why in this case it might be OK.  You could redo the calculation with different assumptions, e.g., $f'(x) \approx 2$ and $f''(x) \approx 3$.  If you do that, you'll find a slightly different optimal value of $h$; instead of $2 \sqrt{\epsilon}$, you'll have a slightly different constant replacing the 2.  So we can see that the 2 comes from the assumptions about the magnitude of $f'(x)$ and $f''(x)$.  But, crucially, if $f'(x)$ and $f''(x)$ aren't too large or too small, then the resulting constant in the optimal choice of $h$ won't be too different.  It might be $5 \sqrt{\epsilon}$ or $0.5 \sqrt{\epsilon}$, but that's not so bad.

And, in practice, we're often dealing with smooth functions where the derivatives don't get too large, so while we don't know exactly what $f'(x)$ or $f''(x)$ are, we might be able to assume they won't be too large.  As a result, taking $h \approx c \sqrt{\epsilon}$ for some small value of $c$ will be optimal.  We don't know what the optimal value of $c$ is exactly, but we can just heuristically pick a small value for $c$, like $c=1$ or $c=2$, and even if that's not exactly optimal, we'll probably do reasonably well -- it won't be too much worse.

So, the final answer is: no, the 2 is not essential.  It's perfectly reasonable to use $h = \sqrt{\epsilon}$ instead of $h = 2 \sqrt{\epsilon}$ -- but you probably want to choose $h$ somewhere in that neighborhood.