--- ******************************************************************* Dr. Jilberto Zamora Saa Dzhelepov Laboratory of Nuclear Problems Joint Institute for Nuclear Research 141980 Dubna, Moscow region, Russia phone: +(496)2162589 (Office) phone: +7(915)4671294 (Mobile) email: jzamorasaa_at_jinr.ru ******************************************************************** El 09-11-2018 5:48 pm, Francesc Salvat-Pujol escribió: > 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
This archive was generated by hypermail 2.3.0 : Thu Dec 20 2018 - 19:36:54 CET