Dear Gian Luigi,
What I need is a primary beam with a rectangular shape in X and Y and with 11 mrad divergence in X and 0 mrad in Y.
Now, THAT should have been the opening statement in the original question.
I understand you are shooting particles starting along the z axis with a certain polar angle theta from 0 to 0.5*DIVX, and azimuthal angle either 0 or pi, such the particle’s initial direction is in the XZ plane.
What you need to do (adapt/optimize as needed) is sample COS(THETA) homogeneously between 1 and 1-COS(0.5D0*DIVX). The sampling is analytical (double check in case I slipped here) and should give you something along the following ~FORTRAN lines:
XI = FLRNDM(XI) ! Homogeneously sampled number in [0,1)
CTH = 1.0D0 - XI * (1.0D0 - COS(0.5D0*DIVX)) ! Sampling cos(theta), note DIVX must be in radians.
STH = SQRT((1.0D0+CTH) * (1.0D0-CTH)) ! sin(theta)
XI2 = FLRNDM(XI2) ! One more homogeneously sampled number in [0,1)
IF (XI2 .LT. 0.5D0) THEN
TXFLK(NPFLKA) = STH ! 50% of the times along +x
TXFLK(NPFLKA) = -STH ! 50% of the times along -x
TYFLK(NPFLKA) = ZERZER ! Zero in FLUKA
TZFLK(NPFLKA) = CTH
Regarding the rectangular (spatial) shape, you can define XMIN,XMAX,YMIN,YMAX in your source.f and then do
XI3 = FLRNDM(XI3) ! One more homogeneously sampled number in [0,1)
XI4 = FLRNDM(XI4) ! One more homogeneously sampled number in [0,1)
XFLK(NPFLKA) = XMIN + XI3 * (XMAX - XMIN)
YFLK(NPFLKA) = YMIN + XI4 * (YMAX - YMIN)
Adapt/correct as needed, and test profusely.
PS: I multiply all divergences by 0.5 to remain in line with the way they are used everywhere else in FLUKA.