In continuation to the previous discussion Biasing in FLUKA - #2 by horvathd,
I tried to modify the source.f, transported photons were registered for some histories, then it got stopped due to some error. Here I am attaching the .out and .err file. Can you please tell me what is the error showing here ?
To test whether the sampling rules are written properly or not, I took a print of XFLKA, YFLKA, ZFLKA for 1000 histories for a single cycle. The plot is according to my intention. (please note a cylindrical volume has been rejected from sampling, since a phantom will be placed there in second step). The sampling is less in outer area and sampling is more as the region is closer to the cylindrical region. To compensate that, weights have been assigned accordingly.
There was one mistake in radius sampling, I have modified it and hence the revised plot and revised log file have been edited in the previous post. I will mail you the source file to you personally.
In regards to the previous post Source sampling with variable weight and probability, I have sent the source file. Can you please confirm whether you have received it ? I tried to sent it to @amario personally. From my end, there was some internet issue while uploading it.
FLUKA stores the total primary weight in the WEIPRI variable.
However, looking at your source routine I have some pointers:
It’s a bit confusing, that you are calculating the area of a circle, using it to set probabilities and weights, but you sample location is a sphere.
You are using the FLRNDM(XDUMMY) function without the XDUMMY variable. (Any double precision variable is acceptable as argument)
I couldn’t figure out why did you set the probabilities / weights as you did, but I think the probabilities needs to be independent from the system. Setting all probabilities to 1.0 is ok.
You used an old source.f template. You can find the current one in /usr/local/fluka/scr/user directory.
You can use a USRBIN scoring to plot a weighted source distribution. For this trick, you need to:
Assign to all regions a non VACUUM material.
Set a photon transport threshold higher that the source energy. This will force all photon primaries stop immediately after they are created and deposit their energy on the spot.
Score ENERGY with USRBIN
If the weighted sampling work correctly, each bin should have the same value (depending on the statistics).
You coded the direction and spherical sampling by hand. There are build in FLUKA function to take care of this.
Based on these points, I reworked your source routine with the new template. This should work correctly: cloud_source_v4.f (11.3 KB)
As far as I remember, I have converted the source.f using the python script and there was no error while building and executing it.
Also, my intention was to sample points as performed in ICRP 144, as was discussed in Biasing in FLUKA - #2 by horvathd, so I tried to divide the sphere in 7 shells and adjusted the sampling (forcing maximum sampling near the cylindrical target) , and since by this the uniformity is lost, I tried to adjust their weights (1 being the outer most shell and then gradually decreasing as it comes closer to the cylinder) . Since this cannot be done using built in cards, I wrote direction sampling and position sampling.
the converted old source routines should work in most cases without any problem, but it is nice to keep them up to date, and more importantly if you start writing a new one from scratch.
I was aware the goal of the source routine, but I didn’t understand how the calculations worked in yours.
Thank you @horvathd, I have just gone through your modified routine. It is written so well and compact, thank you again. Mine was soo childish. Also, I thought the way/rate the surface is increasing, that rate I should incorporate as the probability/ weight. But that was not at all necessary. So, thank you again for the help.
SAMPLING THETA AND PHI
CALL SFECFE(sin_phi, cos_phi)
cos_theta = -ONEONE + TWOTWO * FLRNDM(XDUMMY)
sin_theta = SQRT( ONEONE - cos_theta ** TWOTWO )
In this case, isn’t it needed to take minus sign of the SQRT also (sin_theta) ?
The same notations C1… C7 have been used to define two separate parameters (cumulative and normalizing cumulative probabilities). Is it required to use different notations or will it be automatically replaced each time ?
select_layer = FLRNDM(XDUMMY)
IF (select_layer .LT. C1) THEN
ROUT = R1
RIN = ZERZER
weight = W1
ELSE IF (select_layer .LT. C2) THEN
ROUT = R2
RIN = R1
weight = W2
Isn’t it required to use both GT and LT condition otherwise both IF and ELSEIF conditions are overlapping?
What is the role of XDUMMY in the argument ? What is the difference between the outputs of FLRNDM if this argument is not written ?
Since \theta is sampled from the [0, \pi] range (i.e. -1 \le cos(\theta) \le 1), sin(\theta) is always positive.
This is just a programming trick, to save some memory by reusing the same variables. The values will be correctly updated, but using new variables is also fine.
It is not required to check if the variable is greater that the previous check, because it is tested in sequence.
If the variable satisfies the test in the IF section, then the ELSE IF part is skipped.
If the variable doesn’t satisfies the IF section (i.e. it is already greater that value), the ELSE IF only need to check if it is lower than the next one.
This is also valid for the consecutive ELSE IFs.
From the FLUKA manual:
Note: This function (FLRNDM) requires a dummy double precision variable as argument. The used variable doesn’t influence the returned random number: it is intended to avoid compiler optimizations which could inhibit the actual sampling.
Plus, here are a couple other optional points which can be implemented to further optimize the source routine:
The source routine can be made independent from the number of layers, by using arrays to store the values, and for cycles to do the calculations
Most of the calculation (Volumes, probabilities, weights) will be always the same. These can be moved to the initialization section. The values can be saved between the different calls to the source routine by using the SAVE statement during variable declaration. For example see LFIRST variable.
