Distribution of sample sum given partial noisy observation Let $X$ be a random variable of unknown expectation $\mu$ and standard-deviation $\sigma$.
Let $m$, $n$ $\in \mathbb{N}$ and $p = m + n$.
Let $(X_i) = (X_1, \dots , X_m)$ and $(X'_i) = (X'_1, \dots , X'_n)$ be two samples of $X$, that is the $X_i$ and $X'_i$ are iid as $X$.
We observed a realisation of the first sample $(X_i)$, corrupted by an independent additive noise $B$.
Denoting the observations $\tilde{X} = X  + B$, the available information are the $(\tilde{x}_i)$.
We are interested into predicting
$$ Z = \frac{1}{p} \left( \sum_{i=1}^{m} X_i + \sum_{i=1}^{n} X'_i \right).$$
More precisely, we search bounding quantiles of $Z|(\tilde{x}_i)$.
We denote $\omega = (\mu, \sigma)$ and $\omega_B = (\mu_B, \sigma_B)$ and assume we are able to simulate both associated random variables, $\Omega$ and $\Omega_B$.
I proposed an answer below but would welcome criticism, or insight on how to achieve the actual computation efficiently.
 A: Your notation is non-standard for a sampling theory problem, and it is also a bit overcomplicated; this makes your question more difficult than it needs to be.  Essentially you have a finite population of values, you take a random sample that is corrupted by white noise, and you want to make an inference about the unknown population mean.  This is a variation on standard sampling theory results that lead you to a confidence interval for the mean of a finite population, but this is a tricky case, since it raises some identifiability problems in your analysis.
The first thing to note in this kind of problem is that you have no way from to decompose the mean of the uncorrupted sample values and the mean of the noise; your data does not give information to do this.  So you are only going to be able to do this kind of analysis at all if you are willing to assume that you have genuine noise, in the sense that it has zero mean.  If your "noise" has non-zero mean then the population mean you are trying to estimate is non-identifiable.  So, I am going to proceed under the assumption that $\mu_B = 0$, noting that if this is not the case then there is no solution.
The next thing to note is that, even with this assumption, you encounter another problem because you have no information that allows you to decompose the variance of the uncorrupted sample values from the variance of the noise.  Because you only observe the additive corrupted sample values, you can only identify the sum of these two variances, not their individual contributions.

Obtaining an interval estimate for the population mean: For simplicity, I'm going to eschew your notation, and answer this question using standard notation from sampling theory, so I'm going to take $N$ to be the population size and $n$ the sample size.$^\dagger$  I'm also going to use standard bar notation for sample means.  If we pollute the sample with white noise $B_1, ..., B_n \sim \text{IID N}(0, \sigma_B^2)$ then we have:
$$\begin{equation} \begin{aligned}
\bar{\tilde{X}}_n - \bar{X}_N 
&= \frac{1}{n} \sum_{i=1}^n (X_i + B_i) - \frac{1}{N} \sum_{i=1}^N X_i \\[6pt]
&= \frac{1}{n} \sum_{i=1}^n B_i + \Bigg( \frac{1}{n} - \frac{1}{N} \Bigg) \sum_{i=1}^n X_i - \frac{1}{N} \sum_{i=n+1}^N X_i \\[6pt]
&= \frac{1}{n} \sum_{i=1}^n B_i + \frac{N-n}{Nn} \sum_{i=1}^n X_i - \frac{1}{N} \sum_{i=n+1}^N X_i. \\[6pt]
\end{aligned} \end{equation}$$
Since $\mathbb{E}(B_i) = 0$ we have $\mathbb{E}(\bar{\tilde{X}}_n - \bar{X}_N) = 0$ so that the corrupted sample mean is an unbiased estimator of the uncorrupted population mean.  The variance of the difference is:
$$\begin{equation} \begin{aligned}
\mathbb{V}(\bar{\tilde{X}}_n - \bar{X}_N) 
&= \frac{\sigma_B^2}{n} + \Big( \frac{N-n}{Nn} \Big)^2 n \sigma^2 + \Big( \frac{1}{N} \Big)^2 (N-n) \sigma^2 \\[6pt]
&= \frac{\sigma_B^2}{n} + \Bigg[ \frac{N-n}{Nn} +\frac{1}{N} \Bigg] \frac{N-n}{N} \sigma^2 \\[6pt]
&= \frac{\sigma_B^2}{n} + \frac{N-n}{N} \cdot \frac{\sigma^2}{n}. \\[6pt]
\end{aligned} \end{equation}$$
If we define the parameter $\phi = \sigma_B^2 / \sigma^2$ to measure the noise-to-signal variance ratio, we have:
$$\mathbb{V}(\bar{\tilde{X}}_n - \bar{X}_N) = \Bigg( \phi + \frac{N-n}{N} \Bigg) \frac{\sigma^2}{n}.$$
Assuming an underlying distribution with finite variance and kurtosis, we can apply the CLT to obtain the asymptotic distributions:
$$\frac{\bar{\tilde{X}}_n - \bar{X}_N}{\sigma / \sqrt{n}} \overset{\text{Approx}}{\sim} \text{ N} \Bigg(0,  \phi + \frac{N-n}{N} \Bigg) \quad \quad \quad \frac{\tilde{S}_n^2}{\sigma^2} \overset{\text{approx}}{\sim} (1+\phi) \cdot \text{Chi-Sq}(n-1).$$
This gives the asymptotic distributional result:
$$\frac{\bar{\tilde{X}}_n - \bar{X}_N}{\tilde{S}_n / \sqrt{n}} \overset{\text{Approx}}{\sim} \text{ } \sqrt{\phi + \frac{N-n}{N}} \cdot \frac{\text{St}(n-1)}{1+\phi}.$$
Now, in the special case where you have no noise, you have $\sigma_B = 0$ so that $\phi = 0$, and this gives you the standard pivotal quantity for inference for a finite population mean in a sampling problem.  However, in the more general case of present concern we have a big problem, since the parameter $\phi$ is non-identifiable.  If we are willing to assume a value for this parameter (and thus assume we know the noise-to-signal ratio) then we can use the above result to obtain the confidence interval for the unknown population mean:
$$\text{CI}_n(1-\alpha) = \bar{\tilde{x}}_n \pm \sqrt{\phi + \frac{N-n}{N}} \frac{t_{n-1, \alpha / 2}}{(1+\phi)\sqrt{n}} \cdot \tilde{s}_n.$$
In the case where you are willing to assume a value for the (non-identifiable) noise-to-signal variance ratio, this should give you a reasonable confidence interval for the unknown population mean.  The interval accounts for the population size through the "finite population correction" term, but this is augmented by a correction for the corruption of our sample by the noise term.  We can see that as $\phi \rightarrow \infty$ the interval expands to the whole real number line, as expected by intuition (since this represents the case where we are essentially just observing noise).

$^\dagger$ So my $N$ is your $p$ and my $n$ is your $m$ (and my $N-n$ is your $n$).
A: In the following, $f(y)$ denotes the probability density of any random variable $Y$.
The density of interest is
$$ f(z, \omega, \omega_B | \tilde{\Sigma}) =
   f(z, | \tilde{\Sigma} , \omega, \omega_B)   \;  f( \omega, \omega_B | \tilde{\Sigma} ),
\quad Eq. 1
$$
with $\tilde{\Sigma} = \sum_{i=1}^{m} \tilde{x}_i$.
For sufficiently large $n$ and $m$, the central limit theorem states that
$$ f(z, | \tilde{\Sigma} , \omega, \omega_B) \approx \phi(z, \mu_Z, \sigma_Z) $$
where $z \to \phi(z, \mu_Z, \sigma_Z)$ is the density of a gaussian random variable, with expectation
$$ \mu_Z = \frac{1}{p} \left( \sum_{i=1}^{m} \tilde{x}_i  - m \, \mu_B + n \, \mu \right) , $$
and standard deviation
$$ \sigma_Z = \frac{1}{p} \sqrt{m \, \sigma^2_B  +  n \, \sigma^2} .$$
Applying Bayes formula, we get
$$ f( \omega, \omega_B | \tilde{\Sigma}  )= 
\frac{ f(\tilde{\Sigma} | \omega, \omega_B ) \; f( \omega, \omega_B )}
{f(\tilde{\Sigma})}. 
\quad Eq. 2
$$
As $\Omega$ and $\Omega_B$ are independent, $f( \omega, \omega_B ) = f( \omega) \; f( \omega_B )$.
For large $m$, the central limit theorem implies
$$ f(\tilde{\Sigma} | \omega, \omega_B ) \approx
\phi(\tilde{\Sigma}, m \, (\mu + \mu_B),\sqrt{m \, (\sigma^2 + \sigma_B^2)}).$$
A Monte Carlo estimate of the looked for quantiles can be devised by combining equation 1 and 2.
Indeed, 
$$ Pr(Z|\tilde{\Sigma} < a) = \frac{1}{M} \int_0^a\!\!\!\!\int\!\!\!\!\int 
f(z|\tilde{\Sigma},\omega, \omega_B) \; f(\tilde{\Sigma} | \omega, \omega_B) \; f(\omega) \, f(\omega_B)
\; d\omega d\omega_B dz,
\quad Eq. 3
$$
with 
$$M = \int\!\!\!\!\int f(\tilde{\Sigma} | \omega, \omega_B) \; f(\omega) \, f(\omega_B)
\; d\omega d\omega_B .
\quad Eq. 4
$$
We can for instance


*

*sample $\Omega$ and $\Omega_b$;

*compute $M$ with the empirical equivalent of equation 4;

*compute the outermost integrand of equation 3 for increasing values of $z$ until the sum of those values weighted by the increments in $z$ reaches $a$.

