Numerical computation of a double integral

The function F decays so rapidly that it should be okay to truncate the domain of integration. For example, we have

NIntegrate[Q[v1], {v1, -30, 30}, PrecisionGoal -> 14]/ NIntegrate[Q[v1], {v1, -100, 100}, PrecisionGoal -> 14] - 1

2.66454*10^-15

So the following should have quite a few digits of precision:

F[v_] := -v + 1/(1 + v/(-1 + E^v)) - 1/(1 + (E^v v)/(-1 + E^v))
Q[v1_?NumericQ] := Exp[NIntegrate[F[v], {v, 0, v1}]];
const = NIntegrate[Q[v1], {v1, -30, 30}, PrecisionGoal -> 14]
P[v1_] := Q[v1]/const;

There's no need for finesse or truncation here:

F[v_] = -v + 1/(1 + v/(E^v - 1)) - 1/(1 + (v E^v)/(E^v - 1));
G[v1_?NumericQ] := Exp[NIntegrate[F[v], {v, 0, v1}]]
Gnorm = NIntegrate[G[v1], {v1, -\[Infinity], \[Infinity]}]
(*    2.87582    *)
P[v0_?NumericQ] := G[v0]/Gnorm