I am simulating X-Ray tomography for a sample made of bone, immersed in water, and enclosed in a plastic container. The sample rotates from 0° to 180° during the simulation, and I have two main questions:
Estimating Accumulated Dosage:
To calculate the accumulated dosage over the entire rotation, should I simulate each orientation step separately (e.g., by running a loop with different angle steps) and sum the results? Or is there a more direct approach in FLUKA for handling rotational simulations and calculating accumulated dosage?
USRBIN Resolution vs. Beam Size:
I am using the USRBIN scorer to map the dose distribution. Is there a recommended guideline for setting the number of bins in a USRBIN card relative to the beam size or other simulation parameters?
thank you for your question. I can give this a closer look in the next days, but here are some quick replies that might already solve the issue.
While it is not possible (or it might be, but would be very hard to implement, and I do not recommend doing so) to implement a time-dependent geometry, or to sample rotated geometries in a single run, there might be a simple way that exploits the inherent properties of your geometry. In particular, if I view your simulations correctly, there is no distinguishable difference between rotating your sample and rotating your source about the same axis, but in opposite directions (i.e. clockwise sample rotation of x degrees = counterclockwise source rotation of x degrees).
If that is the case, you can define a source routine which picks up a different angle of the source at each iteration. The result of the simulation will already be summed over all the sampled angles, so all you need to do is make sure the sampling is done as you wish (uniformly between 0 and 180 degrees, I assume)
Otherwise, in a more general scenario where cylindrical symmetry is broken, looping in FLAIR through the variable ANGLE and summing the results is a robust solution.
Some guidelines to keep in mind when choosing the mesh:
a) if your bins do not align with the object you want to score, you might get edge effects. In some cases these are unavoidable. In your case, the mesh can line up well with the sample.
b) once the shape of your scoring region is fixed, the total and the average of the scored variables will not change, as any subdivisions are integrated over
c) while there are many ways of choosing the number of bins for visualizing a dose map, or to score the maximum dose, I would recommend starting with 1/2 sigma of the beam (remember that FWHM is 2.35*sigma, and for a uniform distribution sigma ~ interval/3.5), which should give you a pretty good start for estimating the maximum dose. Along the longitudinal axis, you can pick something in the order of a mm/cm similar to the range of the secondaries you would generate. (of course, in rotating the source we don’t really have a distinction between longitudinal and transversal, so just use the smaller one for both). From there, you can increase/decrease the mesh depending on how smooth your distribution looks, in general you do not want the scored bins to have large relative differences (subjective, but say >30% for illustration) between neighboring bins, unless the underlying geometry is different as well.
I hope this solves the issue. If not, let me know and I can address any concerns in more details in the next few days.
What if I want to place a filter between the beam and the sample? If I rotate the source, I also need to rotate the filter.
With looping the variable ANGLE, is it possible to use some user routine where that loop would be defined, along with the sum of the dose? If so, is there any user routine available that can be modified in order to achieve that?
Also, just to be sure, if I want to estimate the average accumulated dosage in X-ray tomography, would it be the average of the sum of the results (i.e., sum all results first, then compute the average based on a number of voxels)?
If you break cylindrical symmetry, the solution involving rotating the source will not work. Of course, you could still place a cylindrical filter around the setup (the difference in terms of scoring would be quite small, possibly trivial.
It is not a FLUKA routine, but rather a functionality of FLAIR. To do so, make sure you define the variable ANGLE in the appropriate place (at the top of the input usually, under the title). Then after you create a new run in FLAIR, you will find an option “loop” which allows you to create multiple runs each with a different value of angle. You would need to add the score of each run manually at the end.
To answer your question, I am going to assume you proceed through the loop route (if you implement a source routine, you would get the results per primary averaged over all orientations of the sample). If you assume that all angles are equally likely, to achieve the results you want, average dose per primary photon (to be scaled after by the strength of your source), you can do the operations in whichever order you prefer, since the number of source particles incident at any angle is assumed uniform. Just remember to scale the sum by 1/n, where n is the number of angles in your loop, so that your final score is still per source particle. Thus, you can either first find the average across all voxels, then sum together across angles, and divide by number of angles, or first sum together across angles the mesh, then divide by the number of angles, and finally average across the all voxels. The two will give you identical results.
Thank you for your reply. I just have some final questions:
Regarding the assumption that all angles are equally likely in achieving the results—since the number of source particles incident at any angle is assumed uniform—would this still hold if the beam has angular divergence in X and Y and is polychromatic?
When using the LOOP option, do we need to specify the values for primaries, time, and random numbers manually, or are these handled automatically according to the INPUT tab? Also, would I need to run the first “default” input file indicated by < >, or only the input files generated by the loop or simply select all and run?
they are handled automatically. However, if you want to change these or other properties of many files at the same time, you can select all of them, right click on the property you want to select to unlock it, change it, and righ click again to lock it again. You may skip the default run.
I have tried to obtain the dose values (Gy) for a simple scenario with no rotation, but I realised I’m confused.
How can I obtain the dose distribution that can be plotted in 3D?
Would I need to multiply the values by voxel volume and then by the number of photons and conversion factor 1.602176462E-07?
And for the total deposited dose in the entire volume it would be then simply a sum?
I have tried to obtain the dose values for a simple scenario with no rotation, but I’m not sure if the procedure shown in the attached Python code is correct. I would appreciate any help.
If you select the DOSE scoring, it will give you the dose in GeV/g/primary. To get the dose in Gy, you need to do the following:
multiply by 1.602e-7 (Gy/(Gev/g) ) the value of each voxel, such that now each voxel represents Gy/primary in that region
multiply by the photon current, i.e. the number of photons per second in your actual beam. This will give you a map where each voxel is Gy/s
you can plot this as a 3d map now. You do not need to multiply by the voxel volume.
If you want to compute the total dose, in theory you would have to be careful because dose (Gy or GeV/g) are densities, so they don’t add like normal scalars would. So you would have to first multiply each voxel by the voxel mass to get a deposited energy. Then add together all the resulting scaled voxel energies, and divide by the total mass. However, because you are using a cartesian mesh on a uniform sample all the voxel volumes and masses are the same, so you can simplify the operations a lot and simply sum together all the doses as they are provided, and dividing by the number of voxels (i.e. take the average of the voxel doses), and this will give you the dose over a set of voxels.
The python code is mostly correct:
change the line total_dose = np.sum(value*conversion_factor)
to total_dose = np.mean(value*conversion_factor)
which will divide by the number of voxels at the end
in the comments you mention a few times “dose per second”, but your conversion factor includes a time interval (1 s) so while numerically it is identical, the operations would indicate you are trying to compute the integrated dose.
Thank you for your reply, however, I’m still unsure about the meaning of the obtained values.
When I multiply each voxel by 1.602e-7 and the number of photons in 1 second: converted_values = value * number_of_photons * conversion_factor
Then, the value corresponding to one voxel is mostly higher than the total dose computed with: total_dose_mean = np.mean(value * number_of_photons * conversion_factor)
However, my reasoning was that if there is a dose deposited in each voxel and if I sum those values, then I would have the total dose in the entire volume, and I don’t understand why the values in each voxel would be higher than the total dose.
Also, I realised that if I use the following formulas, then I would have the total dose based on the sum of the values: converted_values_with_bins = value * number_of_photons * conversion_factor * (1/number_of_bins) sum_converted_values = np.sum(converted_values_with_bins)
Would that be the way to represent the dose deposition?
Additionally, when I divide the output from the SRBIN Region by the volume of the target and multiply by 1.602e-7 and the number of photons in 1 second, then I get an average value slightly higher than the total dose obtained over the set of voxels:
Does it mean that if my voxel size is decreasing (dxdydz->0), then the values obtained from the USRBIN region would be equivalent to the total dose over the voxel-based volume?
Sorry for so many questions, but as a FLUKA beginner, I’m feeling a bit lost.
Dose is an intensive quantity, like density and opposite to energy deposition. The latter [GeV] will always increase over a larger volume, as you expect. This is not the case for density [g/cm3], dose [GeV/g] or energy deposition density [GeV/cm3]: over a larger volume they can easily decrease, since their value expresses by definition the average over the volume. In this regard, total dose intended as local dose sum over a larger volume is a flawed concept, as earlier discussed.
Thank you for the replies. To summarize, correct me if I’m wrong:
If I want to plot the dose values in 3D and obtain the total dose in the entire volume for a simple case with 1 orientation then I would use the following: converted_voxel_values = voxel_values * number_of_photons * 1.602176462E-07 total_dose = np.mean(converted_voxel_values)
For n orientations (n input files) then I would need to use the following formulas: accumulated_values = summed_voxel_values / n total_dose = np.mean(summed_voxel_values / n)
Is the division by n compensating for the addition of volumes/masses for each orientation?
yes, to convert the value of each voxel from GeV/g to Gy the procedure you implemented is correct
yes, in the case of a uniform (cartesian) mesh, and a uniform density, the “total dose” - meaning the dose computed for the entire sample volume - is the average of all the voxels (you should convince yourself that is the same result as finding the total energy deposited in the sample divided by the total mass of the sample). You can check the result of the computation by adding an extra USRBIN with the same limits in x, y, z but with only a single division in each dimension, i.e. one voxel for the entire sample.
Assuming summed_voxel_values is converted_voxel_values summed across all different runs, your computation of the accumulated_value is correct, but the explanation you gave is not. The division by n is not compensating for volumes or masses, but rather for the number of photons.
If number_of_photons is the total number of photons emitted by the source in 1s, and if during that time interval the source has rotated uniformly across n different orientations, then the number of photons for each orientation is number_of_photons/n. This is the factor of 1/n you are seeing, simply to scale each of the runs by the number of photons emitted from that angle. If instead you are irradiating for 1 s at each angle, then you don’t need the factor of 1/n. Just make sure you understand the source before scaling the simulation result.
the total dose is again correct, see point 2. And to tie it back to point 3, it is actually the fact that you are taking the mean, rather than the sum, that is compensating for the addition of volumes/masses.
I hope these help out, you are mostly on the right track. I recommend you clarify for yourself how many photons are incident on the sample at each orientation, and to convince yourself that the formulas for total dose are simply computing total energy deposited in the sample / total mass of the sample .
