Hola Jilberto,
For inspiration you may first want to examine a slightly simpler 1D case
(sampling from a tabulated energy spectrum), as shown e.g. in
http://www.fluka.org/web_archive/earchive/new-fluka-discuss/6798.html
which contains an example spectrum and a custom source.f routine that
samples from it. Note also the modification of the input file mentioned
in the link. If you make a diff (or e.g. vim -d) with the distributed
source.f (under usermvax/ in your FLUKA directory) you will see the
modifications. In a nutshell, there is:
- a block that loads a source spectrum from the provided file at
initialization and prepares a cumulative function for sampling
later,
- a block which samples the kinetic energy from the histogram and
- a block which carefully sets the sampled kinetic energy and
calculates the modulus of the linear momentum (of course with the
relativistic expression).
You then compile your custom source.f
$FLUPRO/flutil/fff source.f
If there are no errors, this produces an object file which you then
bundle into your custom binary.
$FLUPRO/flutil/lfluka -m fluka -o tuejecutable source.o
Depending on what you do, link with ldpmqmd (see the manual). Run as
usual, but passing "-e tuejecutable" additionally.
If I follow you correctly, you need instead to sample from a tabulated
2D distribution h(E,theta), given in terms of the energy E and a
(polar?) angle theta. So it's a slight generalization of the foregoing.
Stay alert for the integrals and the sampling of the polar angle (see
below).
I assume that the theta dependency is not trivial, i.e., that you have
sufficiently different behaviors in theta at different energies
(otherwise the stuff below can be greatly simplified...).
Presumably you have h(E,theta) tabulated on a reasonably dense grid of
energies E_i and angles theta_j. On paper, at each tabular energy you
would do an integration on the sphere (azimuthal angle aside), à la
g(E_i) = \int_0^pi d theta sin(theta) h(E_i,theta) ,
so that g(E) gives you an energy distribution regardless of angle
(treated below). In practice you evaluate it numerically in as
reasonable a way as you can, but in any event not forgetting about the
solid-angle element. In the same spirit as in the 1D case above, you may
now use the tabulated g(E_i) to generate a cumulative distribution for
the sampling of the kinetic energy.
Once you have a sampled kinetic energy Esampled, you look for the
"active" tabular interval verifying E_i <= Esampled < E_{i+1}. Use a
homogeneously distributed random number (search in the manual for
FLRNDM) to take E_i with probability (E_{i+1}-Esampled)/(E_{i+1}-E_i),
and E_{i+1} otherwise (I do not write down exactly how because every
time I get it wrong, but the idea is along these lines...).
At each tabular energy E_i, you can keep the cumulative angular
distribution in memory (do not forget the solid-angle element when
integrating...) and use it to sample the polar angle in exactly the same
spirit as above. It then remains to sample the azimuthal angle. Unless
you have something else in mind, you can sample e.g. homogeneously in
[0,2pi). With this you have all you need to set the initial direction in
your source.f in terms of the director cosines using the variables
TXFLK, TYFLK, TZFLK.
===============================================
*************** WARNING... ********************
===============================================
Things can go wrong when setting up schemes like the above by hand...
Debug and test intensively before running any simulation. Make e.g. a 2D
histogram of the sampled energies and angles and convince yourself that
you are indeed sampling what was intended.
As a final note, I do not know what degrees of freedom you have in
preparing the input 2D source histogram, but any time there's a polar
angle theta involved, a way to sleep well is to tabulate, integrate, and
sample in terms of mu=cos(theta). Doing the change of variables you
immediately see e.g.
g(E_i) = \int_{-1}^1 d mu h(E_i,mu) ,
that is, the term "sin(theta) d theta" in the solid-angle element
becomes trivially "d mu" and you can mostly forget about faux pas :)
Hope this is reasonably helpful/accurate!
Un saludo,
Cesc
PS: scoring-wise, there should be nothing special.
On Fri, Nov 09 2018, at 12:17 -0300, jilberto Zamora Saa wrote:
>
>Dear FLUKA experts,
>
>How I should do in case I want to define my own beam of particles, let say, I would like to provide a spectrum of muons which depend on energy and Zenit angle and then use it to see the fluence in a detector.
>
>any help is welcome
>
>regards,
>Jilberto
--
Francesc Salvat Pujol
CERN-EN/STI
CH-1211 Geneva 23
Switzerland
Tel: +41 22 76 64011
Fax: +41 22 76 69474
__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info
Received on Fri Nov 09 2018 - 23:30:49 CET