I thank you in advance for your kind support and help.
I am trying to use FLUKA (v. 4-3.1) to perform quick simulations of standard X-ray radiographies. I am only interested in the evaluation of the 2D attenuation pattern of the beam through the object.
Therefore, I want to take advantage of the FLUKA capability to perform ray-tracing through a geometry and save the μ * length product, when RAY pseudoparticles are set as primaries.
I made a simple test, all the files are attached to the present post. I defined a source subroutine in which a series of monenergetic RAYs are sent one-by-one from a point source to each pixel of a detector (which is simply a GADOX screen), scanning all the rows and columns. Then, following this example, I made a Fortran program that saves EXP(-ALGTOT) for each RAY hitting the detector, hence for every pixel. I get the correct attenuation image in about 8 seconds and everything works fine.
The problem is that my real source is an X-ray tube, so the monoenergetic image is not accurate enough. The usual solution to this problem is to obtain monoenergetic images for each energy bin of the spectrum and then calculate the total image as a weighted sum of the various images, were the weights are based on the spectral shape.
However, this way it takes at least 20 min (for 130 energy bins) to get the overall image, a giant fort.10 file is generated, and also it is kinda useless to trace the RAYs through the same geometry. I could as well trace the RAYs only once, get the path lengths in each region and multiply them by the μ(E) vector of each region in a separate program, and then weight them.
I wonder if you could help me to find a better solution. Also, I ask you the following questions:
Is it possible to save only ONE variable in the unformatted Fortran file with RAYs pseudoparticles? For example, to save only the ALGTOT?
Is there any other way to speed up the simulation? Maybe the calculation of the trajectory in the source routine slows down the simulation?
Do you have any tips on opening an unformatted Fortran file (as the one produced by the RAY simulation) with Python or Matlab?
I wonder if the most convenient way to achieve your goal is then to forget about the RAY unformatted file and make use of the src/user/mgdraw_empty.f routine (to be activated by USERDUMP). There, in the ENTRY BXDRAW, you can write on a formatted file the RAY coordinates XSCO, YSCO, ZSCO at the boundary between the region MREG and the region NEWREG (for info, ICODE should be 69). The NCASE variable, giving the event number, allows you to distinguish different rays. This looks like fulfilling your plan. Let me know if I missed something. (If required, a solution not to produce the fort.10 file can be provided).
I can rely for sure on the BXDRAW routine, and I will try to implement it to get what I need. I did not immediately went into BXDRAW because, as I saw from the RAY example, when using RAYs as primaries FLUKA saves all the data anyway. Moreover, it saves exactly the data I need to do the calculation, without a lot of post-processing on my side.
Still, I would like to understand if there are ways to speed up RAY simulations, maybe also allowing to save only a few variables in the fort.10 file or, as suggested by you, even without the generation a fort.10 file.
I would really appreciate any further information on this issue.
In fact, your geometry is pretty simple and the average simulation time per ray is below 0.1 ms, not offering a big margin for improvement. So, the time cost depends on the required resolution, i.e. the number of primaries.
I attach a modified binary version of the concerned core routine, to be renamed kasray.o and then linked into your executable (let me know if you meet compilation compatibility issues), which generates a formatted fort.10 file containing only EXP(-ALGTOT) for each ray entering your region #3.
Hope it helps
I have no problem in the compilation and everything works perfectly. Also, the simulation time is almost half of the previous one (i.e., without the kasray.o file). This solves the issue for sure, so thanks a lot.
In the meanwhile, I wrote an mgdraw.f routine to save the RAY path length in every region (it works with the previous input file if the detector is region is put in #2). I would like to add a small part to this routine, based on the example presented in this post.
As you can see in the attached file, I tried to write a part of the SODRAW in which I scan the MEDFLK array to get the FLUKA codes of the used materials, and then I try to get the arrays of the attenuation coefficients through the SGTTEP function (as in the mentioned example). This would be useful, as I could calculate the attenuation for a polyenergetic beam.
The mgdraw.f is compiled, but when I try to run the simulation i get the following error:
“At line 99 of file emf/sgttep.f
Fortran runtime error: Index ‘0’ of dimension 1 of array ‘ge1’ below lower bound of 1”
I hope you can help me understand the error.
Thanks a lot again, I truly appreciate your help and support
Sorry if I insist, but is there any way to have no fort.10 file output at all?
It would be nice to have to option to only produce the mgdraw.f output, without extra data being generated, to save space on the disk in case of a large dataset of polyenergetic images.
