From: Finetti Noemi <noemi.finetti@aquila.infn.it>

Date: Thu Apr 17 2008 - 10:27:22 CEST

Date: Thu Apr 17 2008 - 10:27:22 CEST

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 MONTALEReceived on Thu Apr 17 14:45:06 2008

*
This archive was generated by hypermail 2.1.8
: Thu Apr 17 2008 - 14:45:06 CEST
*