**From:** Alberto Fasso' (*fasso@SLAC.Stanford.EDU*)

**Date:** Mon Jul 02 2007 - 16:39:48 CEST

**Previous message:**Alfredo Ferrari: "(Further) Fluka2006.3b respin"**In reply to:**Chris Theis: "RE: A weight Problem in uniform sampling"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

On Wed, 27 Jun 2007, Chris Theis wrote:

*> Hello Laurent,
*

*>
*

*>>
*

*>> I try to sample a synchrotron source for Fluka calculations and i use
*

*>> for that source.f . I perform an uniform sampling over my energy range
*

*>> [1:500] keV and this sampling is very efficient (I manage to sample
*

*> all
*

*>> the energies).
*

*>
*

*> Actually your sampling is efficient only at the first glance because you
*

*> have the same chance to sample parts of your spectrum which are not very
*

*> important (=low weight), as you have sampling the important ones (=high
*

*> weight). Depending on your spectrum and the respective weights this
*

*> approach might work, but I'd recommend a different way to elude the
*

*> problems you're running into.
*

Chris, your answer to Laurent may be correct from a general standpoint,

but in practice it is misleading. The "flat" sampling described by Laurent

is a widely used technique to sample from a synchrotron radiation spectrum.

The fact is that the important part of the spectrum _is_ the one with

low weights, while photons with a large weight have very small energies

and are easily absorbed. In most practical problems, it is the most

energetic photons which matter. Their low weight is compensated by

their penetrating power.

*>> My problem comes from the weight calculation. I use the formula w =
*

*>> N(E)/A*(Emax-Emin) where A is the integral of the spectrum over the
*

*>> energy interval [Emin:Emax] and N(E) the number of photons per second
*

*>> per keV as function of energy.
*

The correct formula is w = N(E)*(Emax-Emin). Since FLUKA normalizes

all results to unit weight, the presence of a constant factor 1/A should

not matter (but be careful that it does not lead to underflow!).

See the FLUKA Manual, Chap. 13 (User routines) for a derivation of

the formula, in the description of routine SOURCE.

*>> 2) If an uniform sampling is not desirable for this type of spectra
*

*>> (because of the huge weight interval), do you recommends a better
*

*>> function that 1/N(E) to biase the spectra ?
*

*>
*

*> Instead of sampling the energy interval uniformly and tampering with the
*

*> weights of the particles you could do it the other way round. Namely,
*

*> you sample the important parts of your spectrum (those with the high
*

*> weights) with a higher probability than those that have a lower
*

*> contribution to the result. This way you do not need to change the
*

*> weight and introduce large weight fluctuations.
*

As I said, the important parts of the spectrum ARE NOT those with the

high weights: they contribute very little to the result.

I have used that formula many times and I have never had large weight

fluctuations. These are most likely due to a wrong weight assignment.

Anyway, weight fluctuations can be damped using a weight window: but

this is certainly not necessary here.

*> There are a number of ways that you could to this:
*

*>
*

*> *) Inversion method: You calculate the the integral of your spectrum
*

*> iteratively over all bins. In the end you'll end up with a table of N
*

*> bins, like your original spectrum, and each bin X contains the integral
*

*> from the start-bin up to bin index X. After normalization this will give
*

*> you the tabulated cumulative distribution function of your spectrum. Now
*

*> you can sample this function uniformly and mapping the result back to
*

*> the original spectrum you can determine the energy of the particle that
*

*> should be started.
*

This suggestion is correct, but only provided the spectrum is very strongly

biased towards the high energies. The reason is the same expressed above:

low-energy photons don't matter, and the sampling should favor the high

energy ones. I have used this technique for many years to sample from a

synchrotron radiation spectrum, biasing the spectrum by a function of

energy (generally I multiplied the spectrum by E**8, other people prefer

an exponential exp(E*a), with a = constant). However, I now prefer the

other method, which is much easier to implement. Multiplying by E**8

may create overflow/underflow problems depending on the units used, and it is

often necessary to "adjust" the spectrum by re-normalizing it by a factor

such that values (and their inverses) stay within overflow limits.

Using powers smaller than E**8 becomes very inefficient: I have tried

E**4 and E**6, they give the same results but with much longer computer times.

And values larger than 8 cannot be managed because of the overflows.

*> *) Rejection sampling for your spectrum either with a uniform function
*

*> as a delimiter or a biased function. However, these things are a little
*

*> trickier to implement and the efficiency strongly depends on the
*

*> delimiting function.
*

*>
*

*> *) Transformation method, if your spectrum can be described by an
*

*> analytical function that is bijective and integrable.
*

These are standard sampling techniques, but they are not suited to

this particular kind of problems. There are hundreds of synchrotron radiation

machines in the world, and for shielding problems of this kind the only

sampling techniques used in practices are the two outlined above (or better, in

95% of cases people use the "flat" sampling, and only in a couple of cases

I have seen using the biased-cumulative spectrum).

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@slac.stanford.edu

**Previous message:**Alfredo Ferrari: "(Further) Fluka2006.3b respin"**In reply to:**Chris Theis: "RE: A weight Problem in uniform sampling"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

*
This archive was generated by hypermail 2.1.6
: Mon Jul 02 2007 - 18:07:11 CEST
*