Uncertainty and accumulation during target rotation

Scoring
1

Versions

FLUKA: 4.5.0

Description

Dear FLUKA Experts,

I would like to ask a question regarding the correct statistical uncertainty propagation when accumulating voxel-wise dose (USRBIN) across multiple independent angle runs. I run FLUKA separately for n different target orientations (Flair loop as described here: Accumulated Dosage in X-Ray Tomography - Scoring - FLUKA User Forum).

My understanding is that FLUKA reports, for each voxel/bin, a relative statistical error as:

\epsilon = \frac{\sigma}{\bar{X}}

where \bar{X} is the estimated mean score in the bin and \sigma is the standard deviation of the estimator.

My questions are:

  1. Is it valid to assume that the voxel-wise dose estimators from different orientations are statistically independent if each orientation is simulated as a separate FLUKA run with different random seeds, while keeping the beam parameters identical?

  2. If so, is the correct way to compute the statistical uncertainty of the accumulated dose (per voxel) the following:

First, estimate the accumulated dose for k orientations:

\bar{D}_{\mathrm{acc}} = \frac{1}{n} \sum_{k=1}^{n} D_k

Then convert the relative error for each run into an absolute 1\sigma uncertainty:

\sigma_{k} = \epsilon_{k} D_{k}

Then propagate the uncertainty across orientations by adding variances and applying the same normalisation as for the mean dose:

\sigma_{\mathrm{acc}} = \frac{1}{n} \sqrt{\sum_{k=1}^{n} \sigma_k^2} \quad \Rightarrow \quad \epsilon_{\mathrm{acc}} = \frac{\sigma_{\mathrm{acc}}}{\bar{D}_{\mathrm{acc}}}
  1. If the uncertainties are not independent, would that require an additional covariance term?

Best regards,
Marcin

3

Dear Marcin Sikorski,

Indeed if the only thing that is changing in your simulation is the scoring, then the physics remains exactly the same. You can generally expect large correlations between the bin contents (and hence errors) of two simulations if the only difference between them is a transformation of the USRBIN mesh.

If you think these correlations could be a problem - and it’s not clear to me whether in this case they are - you can change the seeds of the simulations by changing the cycle number in flair.

image

You would then have fully independent simulations.

Your development of error calculation when combining cycles is about right but bear in mind that primary weight oftentimes have to be considered. Take a look to the Monte Carlo lecture of our beginner’s course for more details:

Finally, while it is crucial that you understand how the error is calculated, they are always calculated by FLUKA itself when using built-in scorers such at USRBIN. FLUKA will combine the results from multiple cycles and provide the mean and error of the combined simulations. Is there a reason why you wish to do that? What is your workflow when running FLUKA simulations?

4

Dear Benoit,

Thank you very much for your reply.
My main goal is to quantify the statistical uncertainty of the dose accumulated during target rotation in tomography.

For each angle k, I run a separate input file (5 cycles per run) and obtain a voxel dose estimate D_k(\mathbf{r}) with FLUKA-reported statistical uncertainty ϵ_k​(\mathbf{r})=σ_k​(\mathbf{r})/D_k(\mathbf{r}) where \mathbf{r} denotes the voxel (bin) position in the scoring mesh.

Assuming that:

  • each input uses a different random seed (independent RNG state),
  • the total photon number during the experiment is N_{\mathrm{tot}},
  • photons are uniformly distributed over N_\theta projection angles,

If the dose accumulated is (after conversion to Gy):

D_{\mathrm{acc}}(\mathbf{r}) = \frac{N_{\mathrm{tot}}}{N_\theta} \sum_{k=1}^{N_\theta} D_k(\mathbf{r}).

Would the uncertainty then be correctly evaluated as:

\sigma_{\mathrm{acc}}^2(\mathbf{r}) = \left(\frac{N_{\mathrm{tot}}}{N_\theta}\right)^2 \sum_{k=1}^{N_\theta} \sigma_k^2(\mathbf{r}),

and therefore

\epsilon_{\mathrm{acc}}(\mathbf{r}) = \frac{\sigma_{\mathrm{acc}}(\mathbf{r})} {D_{\mathrm{acc}}(\mathbf{r})} \, ?

Thank you for your help.
Best regards,
Marcin

source_newgen.f (20.8 KB)
spectrum.txt (64.3 KB)
X-Ray-Dose.flair (7.4 KB)
X-Ray-Dose.inp (3.3 KB)

5

Dear Marcin,

Thank you for the details, I understand better what you want to achieve.

I will assume that, for your experimental setup, each angle is equally probable. If it is not the case, you will need to add some weights to account for the fact that some angles are more probable than others. Perhaps the nominal angle is the most probable and other angles should be weighted according to a normal distribution where (\theta-\mu_\theta) is the deviation from the nominal orientation.

I think your approach to calculate the average dose is correct. I would indeed average the dose from all simulations then normalize to the number of expected photons.

Now if you goal is the find the systematic error on the dose measurement related to rotating the tomography target (uncertainty on the positioning), you should really use the standard deviation on the dose of the measurements meaning: \sigma^2_\text{acc} (\textbf{r}) = \frac{N_\text{tot}^2}{N_\theta} \sum^{N_\theta}_{k=1} (D_{k}(\textbf{r}) - \left\langle D_{k}(\textbf{r}) \right\rangle )^2 (note that your formula should have N^1 not N^2 on the denominator). The statistical error of the FLUKA error does not factor into this (unless you would like to calculate the error on this estimation of the systematic error).

To convince yourself, keep in mind that the systematic error on the rotation angle has a “true” value which depends on your physical setup, independent of your simulation. It makes no sense that it would vary as a function of the number of simulated primaries (as in your formula).

Or perhaps you want to rotate the target “on purpose” during the tomography or want to understand random movements? In this case, supposing that the current orientation is independent of the history of previous orientations, or that all orientation are possible over time independently of the orientation at t_0, then, over a long enough period of time, your error on the accumulated dose after an infinite time is zero. You irradiate a target one time, then another, there is no reason the dose would be different between the two “experiments” if they have been exposed for a long enough time. In any case, the FLUKA statistical error is not relevant at your error, as I understand it, is a statistical effect. It is more related to the conceptualization of your problem.

As I said before, the FLUKA statistical error may, however, factor in when calculating the error on your estimation of the systematic error.

6

Dear Benoit,

Thank you very much for your reply.

1.Just to confirm: the assumption that each angle is equally probable is valid in the case of fly-scan tomography when the angular increment (\Delta \theta ) and exposure time (\Delta t) per frame are constant, so that each projection corresponds to the same photon number incoming from the beam?
This remains true even if the sample itself is irradiated differently at different angles due to geometric overlap effects, provided that the total number of incident photons per projection is the same and the USRBIN mesh is identical for all angles?

2.To clarify, the rotation I refer to is the standard sample rotation during tomography (not motion during mechanical loading). The uncertainty I am interested in concerns the uncertainty of Monte Carlo simulations of the total accumulated dose over the full scan. I understand that there are two distinct sources of “uncertainty”:
(i) Monte Carlo statistical uncertainty of the scored dose (decreasing with the number of simulated primaries), and
(ii) positioning/setup uncertainty (e.g., angle mis-positioning), which is a physical effect and does not depend on the number of simulated primaries.
And correct me if I’m wrong, but depending on which we are interested, the denominator would be N^1 or N^2?
My reasoning for estimating the statistical uncertainty was that if the runs are independent then:
Given that:

D_{\mathrm{acc}}(\mathbf r)=\frac{N_{\mathrm{tot}}}{N_\theta} \sum_{k=1}^{N_\theta} D_k(\mathbf r)

And:

\text{Var}(D_{k}(\mathbf{r}))=[\epsilon_k(\mathbf{r})D_{k}(\mathbf{r})]^2

Following:

\mathrm{Var}(aX)=a^2\mathrm{Var}(X)

Then:

\mathrm{Var}\!\left(D_{\mathrm{acc}}(\mathbf r)\right) = \left(\frac{N_{\mathrm{tot}}}{N_\theta}\right)^2 \mathrm{Var}\!\left( \sum_{k=1}^{N_\theta} D_k(\mathbf r) \right) \rightarrow \epsilon_{\mathrm{acc}}(\mathbf r) = \frac{ \sqrt{ \mathrm{Var}\!\left(D_{\mathrm{acc}}(\mathbf r)\right) } }{ D_{\mathrm{acc}}(\mathbf r) }

The reason I’m interested in the uncertainty of the estimate is to determine if the estimate of D_{acc} is valid.

Thank you again and sorry for any trouble.

Regards,
Marcin

7

Dear Marcin Sikorski,

Thank you for the explanation, I did not know it was for fly-scan tomography.

In that case, you should indeed carry out multiple simulations at different angles spread-out evenly. Then, your expression for D_\text{acc} is correct provided that N_\text{tot} is the number number of source particules during the scan (in other words, the integrated beam current). In this context, your error on the average on the simulation error is correct because it is enough to first order to assume that the relative simulation error is the average.

Your approach is, however, only correct if the simulation errors are all equal. It would be better to use a weighted average with weights w_i = \sigma^{-2}. See Statistics and the Treatment of Experimental Data

One problem with your approach is it assumes that the dose varies smoothly and evenly as a function of the angle (in other words, it only works if the second derivative of the dose variations is small compared to the first derivative). To improve your methodology, I would plot the simulated dose values as a function of the angle D=D(\theta). I would fit a sum of sin and cos functions to these points using the simulation errors as error on y. The errors on the amplitudes of the trigonometric will be related to the simulation errors used for the fit. Then, you proceed with an analytical integration of your fit function over the range of angles used for the scan, thereby obtaining D_\text{acc}. With error propagation you obtain the error on the integral, thereby obtaining the error on D_\text{acc}.