The easiest way to think about this is as a chain of two Normal random variables.

You've said your sensor has a known error variance (or at least estimated from a previous experiment) and known no bias (implied it's symmetric about zero, I'll show you how to incorporate non-zero bias either way though). So you have basically defined the distribution of errors as $$e_i \sim N(b, \Delta^2)$$ where $b$ is the known bias of the sensor and $\Delta$ is the known standard deviation of the errors.

If you think of your measurements using this sensor as Normal, then, assuming you know the error for each observation $e_i$, your measurements are simply $$x_i | e_i \sim N(\mu + e_i, \sigma^2).$$ If we remember that the joint distribution is just the product of the conditional distribution and prior then $P(x_i, e_i, \mu) = P(x_i | e_i) P(e_i)$. We can get the marginal distribution for the $x_i$ variables by integrating out $e_i$. Using a common result about Normal distributions (called conjugacy if you want to look it up), we can see $$x_i \sim N(\mu + b, \sigma^2 + \Delta^2)$$.

The maximum likelihood estimate of $\mu$ and $\sigma^2$ are derived from this distribution using standard techniques. In short, $$\hat{\mu} = \bar{x}-b$$ and $$\hat{\sigma^2} = \frac{1}{N}\sum_i (x_i - \bar{x})^2 - \Delta^2.$$

So, if you're interested in the interval for $\mu$ (or $\tau$ in your derivation, I believe), then you just apply these MLEs and the results about averages of normal random variables to get $$ \bar{x} \sim N(\mu, \frac{\sigma^2 + \Delta^2}{N}) $$ and apply this to the typical Normal interval calculation: $$\mu = \bar{x} - b \pm z_\alpha \sqrt{\frac{\hat{\sigma^2} + \Delta^2}{N}}$$ which is simply $$\mu = \bar{x} - b \pm z_\alpha \sqrt{\frac{\frac{1}{N}\sum_i (x_i - \bar{x})^2 - \Delta^2 + \Delta^2}{N}}$$ which simplifies easily enough to $$\mu = \bar{x} - b \pm z_\alpha \sqrt{\frac{\frac{1}{N}\sum_i (x_i - \bar{x})^2}{N}}$$

In your case, you've stated $b = 0$, so that's easy enough to replace!