I have a numerical problem which requires integrating product of two (real, non-negative) functions
$$ \int_0^t f(x)g(x) \,dx $$
where I can compute $f(x)$ easily, but instead of $g(x)$ I only know its integrals over some number of bins,
$$ G_n = \int_{t_n}^{t_{n+1}} g(x)\, dx. $$
The function $g(x)$ is zero everywhere outside the binned interval $[t_0, t_N]$, and the bins are consecutive and inside the $[0,t]$ interval, so that
$$ \sum_{n=0}^{N-1} G_n = \int_{t_0}^{t_N}g(x)\,dx = \int_0^tg(x)\,dx. $$
The bins are all equally wide, too, so $t_{n+1}-t_n$ is constant.
Is there a better way to approximate the true (but unknown) integral from knowing $f(x)$ and all the $G_n$ than just taking $g(x)$ to be a piecewise-constant function?
The $G_n$ are the areas under $g$ from $t_n$ to $t_{n+1}.$ If we pretend, for a moment, that $g$ is piecewise constant over $[t_n, t_{n+1}],$ then from the rectangle area formula, we have that $G_n \approx g(t_n)(t_{n+1}-t_n)=g(t_n)\,\Delta t, $ since your time intervals are constant. From this, we get that $$g(t_n)\approx\frac{G_n}{\Delta t}. $$ In other words, we can re-construct $g$ at the bin boundaries, approximately, from the values of $G_n$ and the time interval width.
Once you have the $g(t_n),$ if you want a smoother approximating function, why not use cubic splines? That would at least give you a smoother approximation to $g$.
Note: the scheme I've set up here has a "gotcha": drift. What you probably ought to do is set up a linear system of simultaneous equations to solve for the $g(t_n),$ so that you can enforce $g(t)=0$ outside $[0,t]$. That will solve the drift problem for you. Alternatively, you could simply try it first, see what kind of error you got, and then spread the error out over the different approximating values for $g(t_n)$.