- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: Francesc Salvat-Pujol <francesc.salvat.pujol_at_cern.ch>

Date: Fri, 9 Nov 2018 21:48:30 +0100

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
*

Date: Fri, 9 Nov 2018 21:48:30 +0100

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:

-- 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_infoReceived on Fri Nov 09 2018 - 23:30:49 CET

*
This archive was generated by hypermail 2.3.0
: Fri Nov 09 2018 - 23:30:51 CET
*