Sampling from a generic distribution (SOURCE)+MGDRAW

Dear Giuseppe, Alberto and Francesco,
thank you very much for your advice! In my SOURCE routine, now I am
calling the subroutine SFLOOD in order to describe an uniform isotropic
flux inside a sphere of radius 7.1 cm and I implemented the Alberto
algorithm (reported in the following) in order to sampling from a
generic distribution and it works! By the way, there was a mistake in
this algorithm:
-----------------------------------------------------------------
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
--------------------------------------------------------------------
must be changed in:
--------------------------------------------------------------------
XI = FLRNDM(DUMMY)
DO 1 I = 0, NPOINT - 1
IF ( (XI .GE. CUMUL(I)).and.(XI.LT.CUMUL(I+1)) ) 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
-------------------------------------------------------
Am I right?

Moreover, by following the example reported in
http://lxmi.mi.infn.it/~battist/fluka_hbk/ I was able to provide the
ntupla of HBOOK in order to verify the distribution of the kinetic
energy, impact point and cosines of "beam'' particles (thanks Francesco
for your advice about how to prepare the executable module linking to
CERNLIB!).

By following the example reported in
http://lxmi.mi.infn.it/~battist/fluka_hbk/node2.html I am not able to
activate the routine MGDRAW (what(1) = or > 100.) with what(3)=2 in
order to activate the entry BXDRAW in MGDRAW because my input data card
does not work 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:
>
>> 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
> .......................................................................
>

-- 
Address: Dott.ssa Noemi Finetti
c/o Dipartimento di Fisica dell'Universita' di L'Aquila
Via Vetoio - 67010 Coppito - L'Aquila - Italy
Phone(office):+39-0862-433051; Fax(department):+39-0862-433033
"NON SI NASCONDE FUORI DEL MONDO CHI LO SALVA E NON LO SA. 
E' UNO COME NOI, NON DEI MIGLIORI." EUGENIO MONTALE