Re: Fluka2008 - manual setting of max. particles in source routine

(smtp1.mi.infn.it [192.84.138.69]); Thu, 12 Feb 2009 20:17:47 +0100 (CET)
Sender: owner-fluka-discuss_at_mi.infn.it

On Tue, 10 Feb 2009, Alfredo Ferrari wrote:

> b) statistical errors are ALWAYS WRONG when an external, pregenerated file
> is used. The limiting example is a file with just few particles that
> you repeat again and again. Statistics will rapidly converge to fake
> small errors, while of course real errors will be dominated by the
> paucity of the initial file. There is no simple way around this issue,
> the user has to apply considerable wisdom and expertise in judging
> the statistical robustness of the results of simulations based
> on pregenerated source files

There is a way to get correct errors when doing calculations in two steps
using an intermediate phase space file, but I admit that it is a bit
complicated and generally people don't do it that way.
You can do this:

First step: run 5 independent jobs using independent random number seeds
             (via RANDOMIZE/WHAT(2)) and create 5 separate phase space files
Second step: run 5 jobs each with a different source which reads one of the
             5 phase space files (you can make a single source, directed by
             a parameter in input option SOURCE to read one of the files)
Then you can calculate errors from the 5 last outputs.

I would like to add a consideration on normalization.
The results of the second step R are normalized to unit weight of the
phase space file particles USED in the second step. But we need them to be
normalized to unit weight of the particles used in the first step.
To do that, the results must be multiplied by the total weight P of all the
particles in the phase space file, and divided by the total weight F of the
particles used in the first step: R = R * P / F
Now, F is easily known, because it is output at the end of the .out file
in the first step:
   Total number of primaries run: N for a weight of: F
   !!! Please remember that all results are normalized per unit weight !!!
But to know P, one must read all the phase space file and calculate the
sum of the weights. One can write a little program to do that, or better
include it in the source routine for the second step.
As Alfredo rightly said, the sampling from the phase space file must not
be done sequentially, but randomly: therefore it is necessary in the
initialization phase (IF (LFIRST) ....) to read the content of the phase
space file (X(I), Y(I), Z(I), ... and WEIGHT(I)) into a number of arrays.
It is easy at this stage to keep a sum P of the WEIGHT. It is even possible
then to get the results of the second step alreadt normalized correctly:
Instead of
       WEIPRI = WEIPRI + WTFLK (NPFLKA)
one can put
       WEIPRI = WEIPRI + WTFLK (NPFLKA) * P / F
(where F is copied from the output of the first step, and F calculated
at initialization as explained above).

Alberto

-- 
Alberto Fasso`
SLAC-RP, MS 48, 2575 Sand Hill Road, Menlo Park CA 94025
Phone: (1 650) 926 4762   Fax: (1 650) 926 3569
fasso_at_slac.stanford.edu