Re: Sampling from a generic distribution (SOURCE)+MGDRAW

From: Alberto Fasso' <fasso@SLAC.Stanford.EDU>
Date: Thu Apr 17 2008 - 17:57:14 CEST

Dear Noemi,

there was no mistake in my algorithm, although your modified version is equ=
ally
correct. The DO loop scans the cumulative function, which is monotonically=
=20
increasing, one interval at the time. We know that XI (which is a number be=
tween
0 and 1) is located in one of the intervals of CUMUL (which is also between
0 and 1). When the loop finds that XI .GE. CUMUL(I), we KNOW that it is
in that interval, so we are automatically sure that XI is also .LT. CUMUL(I=
+1),
i.e. we know that the XI value is necessarily between the two extremes of t=
he=20
interval.
If we continued to scan, of course XI .LT. CUMUL(I+1) would not be valid
anymore: but we stop here, of course (GO TO 2!)
Therefore, the program finds that the first part of the IF is satisfied, an=
d
never checks the second part. So your addition is useless, but not wrong.
It just does nothing.

Concerning OUTLEVEL: this is a very old command which has been removed
from the program a long time ago. Also this does no harm: it is just ignore=
d
by the program. But I suggest that you remove it.

Concerning your problems with USERDUMP: your input card
USERDUMP 101. 2.0 SEAMU
is incorrect.
I don't know if you are using free format or not. If you do, you cannot lea=
ve
any WHAT empty. If SEAMU is the SDUM, then in free format you must write:
USERDUMP 101. 2.0 0.0 0.0 0.0 0.0 SEAMU
or
USERDUMP 101. 2.0, , , , ,SEAMU
If instead you are using non-free format, you can leave some WHATs blank
but they and SDUM must be properly aligned:
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+=
=2E...8
USERDUMP 101. 2.0 SEAMU

Kind regards,

Alberto

On Thu, 17 Apr 2008, Finetti Noemi wrote:

> Dear Giuseppe, Alberto and Francesco,
> thank you very much for your advice! In my SOURCE routine, now I am calli=
ng=20
> the subroutine SFLOOD in order to describe an uniform isotropic flux insi=
de a=20
> sphere of radius 7.1 cm and I implemented the Alberto algorithm (reported=
 in=20
> the following) in order to sampling from a generic distribution and it wo=
rks!=20
> By the way, there was a mistake in this algorithm:
> -----------------------------------------------------------------
> XI =3D FLRNDM(DUMMY)
> DO 1 I =3D 1, NPOINT - 1
> IF (XI .GE. CUMUL(I)) THEN
> * come here when cumul(i) =3D< xi < cumul(i+1). Interpolate
> SAMPLE =3D ENERGY(I) + (ENERGY(I+1) - ENERGY(I))
> & * (XI - CUMUL(I)) / (CUMUL(I+1) - CUMUL(I))
> GO TO 2
> END IF
> 1 CONTINUE
> 2 CONTINUE
> --------------------------------------------------------------------
> must be changed in:
> --------------------------------------------------------------------
> XI =3D FLRNDM(DUMMY)
> DO 1 I =3D 0, NPOINT - 1
> IF ( (XI .GE. CUMUL(I)).and.(XI.LT.CUMUL(I+1)) ) THEN
> * come here when cumul(i) =3D< xi < cumul(i+1). Interpolate
> SAMPLE =3D ENERGY(I) + (ENERGY(I+1) - ENERGY(I))
> & * (XI - CUMUL(I)) / (CUMUL(I+1) - CUMUL(I))
> GO TO 2
> END IF
> 1 CONTINUE
> 2 CONTINUE
> -------------------------------------------------------
> Am I right?
>
> Moreover, by following the example reported in=20
> http://lxmi.mi.infn.it/~battist/fluka_hbk/ I was able to provide the ntup=
la=20
> of HBOOK in order to verify the distribution of the kinetic energy, impac=
t=20
> point and cosines of "beam'' particles (thanks Francesco for your advice=
=20
> about how to prepare the executable module linking to CERNLIB!).
>
> By following the example reported in=20
> http://lxmi.mi.infn.it/~battist/fluka_hbk/node2.html I am not able to=20
> activate the routine MGDRAW (what(1) =3D or > 100.) with what(3)=3D2 in o=
rder to=20
> activate the entry BXDRAW in MGDRAW because my input data card does not w=
ork=20
> properly for two reasons:
> 1) the card: OUTLEVEL 1.0 7.0
> seams to be an unknown control code and in the out file I read:
> "*** Unknown control code: OUTLEVEL, ignored ***"
> 2) the card: USERDUMP 101. 2.0 SEAMU
> aborted the run (core dumped) and in the log file I read:
> "sfe: formatted io not allowed
> apparent state: unit 15 named fort.15
> lately writing sequential formatted external IO
> nditions/requests:"
> Could you please help me to solve these problems?
> Kind regards,
> noemi
>
> Alberto Fasso' wrote:
>
>> On Fri, 28 Mar 2008, Finetti Noemi wrote:
>>=20
>>> One more question ... from Fluka on-line manual I read something about=
=20
>>> "Sampling from a generic distribution" but it is not clear to me how I =
can=20
>>> find with Fluka the inverse value XX=3DF-1(t) of the cumulative distrib=
ution=20
>>> of a given function ("fu" in my source subroutine) by interpolation.
>>> noemi
>>=20
>>=20
>> Dear Noemi,
>>=20
>> 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 p=
art
>> of the source routine, you calculate the normalized cumulative function.=
=20
>> (Here
>> I have done the integral by simple trapeziums, but you can use something=
=20
>> 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 inter=
val
>> it is located. Interpolate linearly between the two energies correspondi=
ng
>> 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=20
>> faster search of the cumulative function without the need of a DO loop.
>>=20
>> Kind regards,
>>=20
>> Alberto
>>=20
>>=20
>> PARAMETER(NPOINT =3D 100)
>>=20
>> DIMENSION ENERGY(0:NPOINT)
>> DIMENSION SPECTR(0:NPOINT)
>> DIMENSION CUMUL(0:NPOINT)
>>=20
>> * initialize npoint+1 values of energy and corresponding spectrum
>> DATA (ENERGY(I, I =3D 0, NPOINT) /123.D0, 456.D0, 789.D0, .... ,
>> & .... , ...., ...., ...., ...., ...., ...., ...., ....,
>> & .... /
>> DATA (SPECTR(I, I =3D 0, NPOINT) /987.D0, 654.D0, 321.D0, .... ,
>> & .... , ...., ...., ...., ...., ...., ...., ...., ....,
>> & .... /
>>=20
>> IF(LFIRST) THEN
>> ........ * calculate the cumulative function (here simple trapezium=20
>> integration)
>> CUMUL(0) =3D ZERZER
>> DO I =3D 1, NPOINT
>> CUMUL(I) =3D CUMUL(I-1) + HLFHLF * (SPECTR(I) + SPECTR(I-1))
>> & * (ENERGY(I) - ENERGY(I-1))
>> END DO
>> * normalize the cumulative function to 1
>> SUM =3D CUMUL(NPOINT)
>> DO I =3D 1, NPOINT
>> CUMUL(I) =3D CUMUL(I)/SUM
>> END DO
>> ........
>> END IF
>> ........................................................................
>> ........................................................................
>> * sample an energy from the spectrum
>> * get a random number between 0 and 1
>> XI =3D FLRNDM(DUMMY)
>> DO 1 I =3D 1, NPOINT - 1
>> IF (XI .GE. CUMUL(I)) THEN
>> * come here when cumul(i) =3D< xi < cumul(i+1). Interpolate
>> SAMPLE =3D 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) =3D SAMPLE
>> .......................................................................
>>=20
>
>

--=20
Alberto Fass=F2
SLAC-RP, MS 48, 2575 Sand Hill Road, Menlo Park CA 94025
Phone: (1 650) 926 4762 Fax: (1 650) 926 3569
fasso@slac.stanford.edu
--1334196310-982552566-1208447834=3D:30002--
Received on Thu Apr 17 22:40:19 2008

This archive was generated by hypermail 2.1.8 : Thu Apr 17 2008 - 22:40:19 CEST