Re: use of CERNLIB with Fluka (source)

On Fri, 28 Mar 2008, Finetti Noemi wrote:

> One more question ...
> from Fluka on-line manual I read something about "Sampling from a generic
> distribution" but it is not clear to me how I can find with Fluka the inverse
> value XX=F-1(t) of the cumulative distribution of a given function ("fu" in
> my source subroutine) by interpolation.
> noemi

Dear Noemi,

I send you here a simple example. You need to prepare two tables: one
of energies and one of the corresponding spectral values (that you can
calculate apart using your "fu" function). Then, in the initialization part
of the source routine, you calculate the normalized cumulative function. (Here
I have done the integral by simple trapeziums, but you can use something more
sophisticated if you want).
At the time to sample an energy, you take a random number and compare it
with the values of the cumulative function until you find in which interval
it is located. Interpolate linearly between the two energies corresponding
to the extremes of that interval.
Of course you can make variants of this scheme: for instance instead of
using DATA you can initialize the energy and spectrum tables inside the
IF(LFIRST) block. And if the energies are equidistant, you can make a faster
search of the cumulative function without the need of a DO loop.

Kind regards,

Alberto

       PARAMETER(NPOINT = 100)

       DIMENSION ENERGY(0:NPOINT)
       DIMENSION SPECTR(0:NPOINT)
       DIMENSION CUMUL(0:NPOINT)

* initialize npoint+1 values of energy and corresponding spectrum
       DATA (ENERGY(I, I = 0, NPOINT) /123.D0, 456.D0, 789.D0, .... ,
      & .... , ...., ...., ...., ...., ...., ...., ...., ....,
      & .... /
       DATA (SPECTR(I, I = 0, NPOINT) /987.D0, 654.D0, 321.D0, .... ,
      & .... , ...., ...., ...., ...., ...., ...., ...., ....,
      & .... /

       IF(LFIRST) THEN
       ........
* calculate the cumulative function (here simple trapezium integration)
           CUMUL(0) = ZERZER
           DO I = 1, NPOINT
              CUMUL(I) = CUMUL(I-1) + HLFHLF * (SPECTR(I) + SPECTR(I-1))
      & * (ENERGY(I) - ENERGY(I-1))
           END DO
* normalize the cumulative function to 1
           SUM = CUMUL(NPOINT)
           DO I = 1, NPOINT
              CUMUL(I) = CUMUL(I)/SUM
           END DO
       ........
       END IF
........................................................................
........................................................................
* sample an energy from the spectrum
* get a random number between 0 and 1
       XI = FLRNDM(DUMMY)
       DO 1 I = 1, NPOINT - 1
          IF (XI .GE. CUMUL(I)) THEN
* come here when cumul(i) =< xi < cumul(i+1). Interpolate
             SAMPLE = ENERGY(I) + (ENERGY(I+1) - ENERGY(I))
      & * (XI - CUMUL(I)) / (CUMUL(I+1) - CUMUL(I))
             GO TO 2
          END IF
  1 CONTINUE
  2 CONTINUE
.......................................................................
       TKEFLK (NPFLKA) = SAMPLE
.......................................................................
Received on Sat Mar 29 12:10:27 2008