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.
