RE: A weight Problem in uniform sampling

From: Alberto Fasso' (fasso@SLAC.Stanford.EDU)
Date: Mon Jul 02 2007 - 16:39:48 CEST

  • Next message: Yung-Shun Yeh: "Re: about muon coincidence in undreground lab"

    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 Fasso`
    SLAC-RP, MS 48, 2575 Sand Hill Road, Menlo Park CA 94025
    Phone: (1 650) 926 4762   Fax: (1 650) 926 3569

  • Next message: Yung-Shun Yeh: "Re: about muon coincidence in undreground lab"

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