How to simulate spread -out bragg peak

Dear experts,
I want to simulate spread -out bragg peak of proton therapy, a source.f should be defined (FLUKA manual P258), but I do not know the fortran code and how to define the source?
Best Regards

Dear Mike,

Indeed you will have to customize your own source.f.

I understand you have a well defined distribution of primary proton energies. In what form do you have it? Is it a spectrum tabulated in a file? Do you want to consider discrete energies with certain weights? Do you have an analytical expression for it?

Please clarify so as to better gauge the reply.


Dear cesc,
It was proved that the source file source.f should be written by the user of FLUKA, where Idefine the weightof the proton energy (a file named SOBP.dat).How to read a spectrum from an existing file SOBP.dat? In others words, how to read and properly normalized at inizialization in source.f?

Dear Mike,

source.f (10.2 KB)

Attached an example source.f (make a diff with user/source.f to see
the basic changes):

  • Line 45: set number of lines in SOBP.DAT (there are more automatic
    approaches but nevermind now).

  • Lines 63-79: load the spectrum from SOBP.DAT (I think this is
    case-sensitive). It is assumed that energies are sorted in ascending
    order, are given in GeV, and that weights are non-negative. The
    cumulative distribution is evaluated using a trapezoidal integration
    (refine as/if needed). The output file *_fort.99 contains the read in
    SOBP.DAT and the evaluated cumulative distribution.

  • Line 82: define particle species ID = 1 (proton).

  • Line 152-164: sample kinetic energy by sampling a random value
    homogeneously in [0,1), looking up the interval where the cumulative
    distribution bounds the sampled number and linearly interpolating so
    as to yield a kinetic energy. If you uncomment line 163 the sampled
    values will be printed out in *_fort.63 (useful to debug and to cross
    check you get what you want).

You’ll need to compile the source.f and produce your custom binary
(check out the Compile tab in Flair).

If you look in source.f you’ll see that the initial position
(XFLK,YFLK,ZFLK) and direction (TXFLK,TYFLK,TZFLK) are taken from
the BEAMPOS card, without applying any spread, regardless of what you
requested in the BEAM card. You can of course override them. Refine
as needed.

See if this flies for you.


PS: a “purist” comment… The procedure above solves the
inverse sampling problem by linearly interpolating the cumulative
distribution. This implies that the sampled probability distribution
(the derivative of the cumulative) will be step-wise constant. I.e., if
you make a histogram (with bins much thinner than the stepsize in
SOBP.DAT) of many sampled kinetic energies you’ll see steps.
In practice, this does not matter, provided your SOBP.DAT is
tabulated on a reasonably dense grid.

PS2: it is assumed above that you have a continuous proton-energy
distribution tabulated at some energy grid. If what you have are
discrete proton energies, the scheme above does not apply.