# Sampling from continuous distribution

Dear fluka users and experts:

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?

Best regards,
Honghu

source.f (9.3 KB)

Dear Honghu,

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).

With kind regards,

Cesc

Dear Cesc,

Sorry for the ambiguities with my expression, it is an exponential distribution with parameter T:
f(x) = (1/T) exp(-x/T)*

I checked the math, and replace the following lines with yours.

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
#2 0x7fffa8140b39

Is there anything wrong with my understanding? Could you please explain how the particle energy and probability pass to the program.

Best regards,
Honghu

Hello,

you might be interested to take a look at Sampling from energy spectra

Cheers
Chris

Dear Chris,

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.

Best regards,
Honghu

I recheck the source.f file, the flair can now correctly compile and run the executable program.

However, the output of spectrum scoring is all Zero. Assuming no energy information of primaries pass to the program.

Dear Honghu,

it seems there is a typo in your source routine.
Your random number is X1, but when you calculate XE, you use XI.

Cheers,
David

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…

Dear David,

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.

phoyield.inp (2.5 KB)

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.

Dear cesc:

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.

input-0: Plot #5.eps (83.8 KB)
source.f (9.4 KB)

Hello,

It looks all right. You can delete the “PARAMETER (NMAX=1000)” line, as it is no longer needed.

For good measure, replace 0 by 0.0D0 and 1000 by 1000.0D0.

Cheers,

Cesc

Thanks for your help!

Honghu