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

From: Alfredo Ferrari <alfredo.ferrari@cern.ch>
Date: Fri Apr 18 2008 - 09:30:59 CEST

Hi all

I believe Noemi is right. There is a mistake in Alberto example, indeed it
should be or like Noemi is proposing, or (I believe this is the intended
original version of Alberto)

XI = FLRNDM(DUMMY)
DO 1 I = 1, NPOINT - 1
IF (CUMUL(I+1) .GT. XI) 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

note the

IF (XI .GE. CUMUL(I)) THEN

going into

IF (CUMUL(I+1) .GT. XI) THEN

where, apart the obvious mistake, GT instead of GE should be used to be
coherent with the FLUKA random number generator which is zero included,
one excluded

Ciao
Alfredo

On Thu, 17 Apr 2008, Alberto Fasso' wrote:

> 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.
>
> 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
>>
>> 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:"
>> 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--
>

--
+----------------------------------------------------------------------+
|  Alfredo Ferrari                ||  Tel.: +41.22.76.76119            |
|  CERN-AB                        ||  Fax.: +41.22.76.69474            |
|  1211 Geneva 23                 ||  e-mail: Alfredo.Ferrari@cern.ch  |
|  Switzerland                    ||                                   |
+----------------------------------------------------------------------+
Received on Fri Apr 18 09:39:43 2008

This archive was generated by hypermail 2.1.8 : Fri Apr 18 2008 - 09:39:43 CEST