I am trying to sample a continuous distribution of f(x)=E2exp(-E/T)(1/T3).
Starting from the distribution of f(x)=exp(-E/3) according to the Fluka course presentation, and implement it in the source.f file, please find it in the attachment. The program terminated in the middle.
Supposing incorrections with source.f file, can anyone please help with this? Or providing similar source.f file so I can refer to?
Thanks in advance!
You should clarify the notation (I don’t see any x in your f(x),
parameters appear loosely and normalization is not clear), and
exactly what you want to sample. It looks like you want to sample
from an exponential distribution with parameter T:
f(x) = (1/T) exp(-x/T).
Forget the approach in your source.f (there’s an underlying confusion
between tabulating a probability density function and sampling from
it). Sampling from an exponential distribution is analytical, so there’s
no need to tabulate anything. This distribution is normalized as it
stands:
\int_0^\infty dx (1/T) exp(-x/T) = 1
And moreover you can solve the inverse sampling problem analytically.
Let xi be a uniformly distributed random number in [0,1) i.e. what you
get from FLRNDM(XDUM). You then need to solve:
\int_0^x dx' (1/T) exp(-x'/T) = xi
This yields (double check it)
x = - T * ln(1-xi)
Or, since 1-xi and xi are both distributed homogeneously, you may take
x = -T * ln (xi)
Your x appears to be an energy, so you may want to discard sampled
values outside of a certain domain [E1,E2]. Try along the following
lines (maybe more elegant with a while loop but it does not matter
now):
E1 = <insert your lower energy bound here>
E2 = <insert your upper energy bound here>
666 XI = FLRNDM(XDUM)
XE = - T * LOG(XI)
IF (XE.LT.E1 .OR. XE.GT.E2) GO TO 666
This should work, but double check the math. It is more
expedient/cleaner/efficient to do the little bit of algebra above
restricting x in [E1,E2] from the onset (which affects both the
normalization and the inverse sampling equation) so as not
to waste time rejecting sampled values (depending on what
your E1 and E2 are this can be critical).
Within running, flair prompts an error “FINISHED WITH ERRORS” and the input.log shows:
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error: #0 0x1350ba3b4 #1 0x1350b9ad6 #2 0x7fffa8140b39
Is there anything wrong with my understanding? Could you please explain how the particle energy and probability pass to the program.
Thanks for your reply. After reading the ReadMe, it seems to me that the program is probably suitable for sampling discrete distribution, not for continuous distribution.
After David’s comment the sampling is fine (if the 0 to 100000 GeV domain is really what you want…). The problem remains that as it stands now, the block of code intended for sampling is called only at initialization (it’s called only once when LFIRST is true)…
Note that you have not told us what are your primary particles and what you are scoring.
PS: replace 3 with 3.0D0, 1 with 1.0D0, etc for double-precision…
I corrected the typo in the source routine. Rerun the program within the flair.
The scoring spectrum of primaries is a mono-energtice 0.1GeV, which is exactly the energy setting in the BEAM card, not the continuous distribution as we wished.
It seems to me the source routine didn’t properly pass the energy information to the program. Can you help with this?
Please find the source routine and the input file in the attachment.
As pointed out by @cesc just before, you are doing your sampling only once at initialization, under the condition
| First call initializations:
IF ( LFIRST ) THEN
Instead you should move it down, right after the line
… to this point: don’t change anything
Moreover, your XE variable does not mean anything to FLUKA, so you cannot expect that its value is considered by the code. If it’s meant to carry the kinetic energy in GeV, you should assign it to the relevant FLUKA variable:
TKEFLK (NPFLKA) = XE
and also redefine later the particle momentum accordingly:
PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
& + TWOTWO * AM (IONID) ) )
Note that the assumed particle type is the one set in the BEAM card.
I am simulating an exponential distribution of electron. By scoring the primary electron spectrum to check the source routine.
The source routine is modified according to your and @ceruttif suggestion: change the E2 to an appropriate value(0.1GeV, 100MeV), and put in the right place. Moreover, transfer energy unit from GeV to MeV before passing to TKEFLK.
And now the scoring primary spectrum is more like an exponential distribution between E1 and E2, however not meticulous enough, wondering if there still something I missed or misunderstood? Please find the spectrum and source routine in the attachment.
