4

User routines

Unlike some other Monte Carlo particle transport codes, FLUKA gets its input mainly from a simple file. It offers a rich choice of options for scoring most quantities of possible interest and for applying different variance reduction techniques, without requiring the users to write a single line of code. However, although normally there is no need for any “user code”​, there are special cases where this is unavoidable, either because of the complexity of the problem, or because the desired information is too unusual or too problem-specific to be offered as a standard option.

And on the other hand, even when this is not strictly necessary, experienced programmers may like to create customised input/output interfaces. A number of user routines (available on LINUX and UNIX platforms in directory usermvax) allow to define non-standard input​ and output​, and in some cases even to modify to a limited extent the normal particle transport. Most of them are already present in the FLUKA library as dummy or template routines, and require a special command in the standard input file to be activated. Users can modify any one of these routines, and even insert into them further calls to their own private ones, or to external packages​ (at their own risk!). This increased flexibility must be balanced against the advantage of using as far as possible the FLUKA standard facilities, which are known to be reliable and well tested.

Description of available user routines

ABSCFF: user-defined ABSorption CoeFFicient

​ ​

Summary

Argument list (all variables are input only)

WVLNGT :

photon wavelength​ (in cm)

OMGPHO :

angular frequency​ (\(\omega = 2 \pi \nu\)) of the photon (in s\(^{-1}\))

MMAT :

material!index

Function ABSCFF returns a user-defined absorption coefficient for optical photons ​. It is activated by setting WHAT(2) < -99 in command OPT–PROP, with SDUM = blank. See and Generating and propagating optical photons for more information.

COMSCW: weighting deposited energy or stars

​ ​​

Summary

Argument list

IJ :

particle type (1 = proton, 8 = neutron, etc.: see code in Particles transported by FLUKA). Input only, cannot be modified.

XA,YA,ZA  :  current particle position
MREG :

current geometry region

RULL :

amount to be deposited (unweighted)

LLO :

particle generation. Input only, cannot be modified.

ICALL :

internal code calling flag (not for general use)

This function is activated by option USERWEIG with WHAT(6)> 0.0. Energy and star densities obtained via SCORE and USRBIN, energy and stars obtained via EVENTBIN and production of residual nuclei obtained via RESNUCLEi are multiplied by the value returned by this function. The user can implement any desired logic to differentiate the returned value according to any information contained in the argument list (particle type, position, region, amount deposited, particle generation), or information available in COMMON SCOHLP​ (binning number, type of scored quantity). The scored quantity is given by the flag ISCRNG​ (in SCOHLP​):

ISCRNG =

1 \(\longrightarrow\) Energy density binning

ISCRNG =

2 \(\longrightarrow\) Star density binning

ISCRNG =

3 \(\longrightarrow\) Residual nuclei scoring

The binning/detector number is given by JSCRNG​ (in SCOHLP​) and is printed in output between the estimator type and the detector name:

Res. nuclei n. 3  "any-name" , "high" energy products, region n. 4
R-Phi-Z  binning n. 5  "other-name" , generalised particle n.    1
Note that an detector of residual nuclei can have the same JSCRNG number as a binning (use the value of ISCRNG to discriminate). Further information can be obtained including COMMON TRACKR​ (for instance particle’s total energy, direction cosines, age​). TRACKR contains also special user variables (both integer and in double precision) which can be used to save information about particles which have undergone some particular event. If data concerning the current material are needed, it can be accessed as MEDIUM(MREG) if file (FLKMAT)​ is included. Indeed, a common simple application of COMSCW is to score dose according to the local density (especially useful to get the correct average dose in bins straddling a boundary​ between two different media):
         ..................
         INCLUDE '(FLKMAT)'
         INCLUDE '(SCOHLP)'
         ..................
* ========== In order to compute doses ========= *
*        Medium(n) is the material number of region n
*        Rho(m) is the density of material m (in g/cm3)
*        Iscrng = 1 means we are depositing energy (not stars)
         IF ( ISCRNG .EQ. 1 ) THEN
*           to get dose in Gy (elcmks is the electron charge in C)
            COMSCW = ELCMKS * 1.D12 / RHO (MEDIUM(MREG))
         ELSE
*           oneone is defined as 1.D0 in include DBLPRC
            COMSCW = ONEONE
         ENDIF
         ..................

Note that the variables in the argument list, with the exception of IJ, LLO and ICALL, are local copies of those used for particle transport, and therefore can be modified to have an effect on scoring, without affecting transport.

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

Note: setting the variable LSCZER = .TRUE. before RETURN (LSCZER is in COMMON SCOHLP), will cause zero scoring whatever the value returned by COMSCW. This is more efficient than returning a zero value.

DETSCW: weighting/applying a resolution to pulse-height depositions

​

Summary

Argument list:
     EDPDET : energy deposition in the "detector" regions (GeV), the one
              that is going to be scored in the Jscrng_th pulse height
              (DETECT) estimator, if the conditions are fulfilled
     EDPCOI : energy deposition in the "(anti)coincidence regions (GeV),
              of the Jscrng_th pulse height (DETECT) estimator
     ICOIFL : flag for anti-coincidence (-1), no coincidence (0, in this
              case Edpcoi is meaningless and possibly undefined),
              coincidence (+1) (input only)

Subroutine DETSCW allows to modify a pulse-height event-by-event deposition before it is scored, for example by applying a user defined detector resolution. Both the wished for deposition (EDPDET) and the (anti)coincidence one (if any, EDPCOI) can be modfied

DFFCFF: user-defined DiFFusion CoeFFicient

​ ​

Summary

Argument list (all variables are input only)

WVLNGT :

photon wavelength​ (in cm)

OMGPHO :

angular frequency​ (\(\omega = 2 \pi \nu\)) of the photon (in s\(^{-1}\))

MMAT :

material!index

Function DFFCFF returns a user-defined diffusion coefficient for optical photons. It is activated by setting WHAT(3) < -99 in command OPT–PROP, with SDUM = blank. See OPT–PROP and Generating and propagating optical photons for more information.

ELEFLD: definition of an electric field

Summary

Argument list

X, Y, Z  :  current position (input only)
T  :  particle current age (s) (input only)
ETX, ETY, ETZ  :  direction cosines​ of the electric field vector (returned)
E :

electric field intensity​ in MV/m (returned). If the field at X,Y,Z,T is zero, a negative value implies zero field; \(|E|\) (cm) overrides the default safe distance a particle can travel before the field is non-negligible.

NREG :

current region (input only)

IDISC :

if returned = 1, the particle will be discarded

ELEFLD is activated by option ELCFIELD with WHAT(4–6) = 0.0 and is used to return intensity and direction​ of an electric field based on the current position, time, and region. It is called only if the current region has been flagged as having a non-zero electric field by option ASSIGNMAt, with WHAT(5)= 2.0 or = 3.0.

The electric field spatial distribution is often read and interpolated from an external field map​. Note that in any case the direction cosines must be properly normalised​ in double precision (e.g. ETX = SQRT(ONEONE - ETY**2 - ETZ**2)), even if E = 0.0.

Please read carefully the notes on option ELCFIELD.

If name-based input is being used, the name corresponding to NREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

ENDSCP: ENergy density DiStributed — Change of Positions

​ ​​

Summary

Argument list

IJ :

particle type (input only)

NTRUCK :

number of step points (input only)

XTRUCK, YTRUCK, ZTRUCK  : particle step points, can be modified by user
TTRUCK :

step length; should be modified to remain consistent with XTRUCK, YTRUCK, ZTRUCK.

MREG :

region number (input only)

LLO :

particle generation​ (input only)

ICALL :

internal code calling flag (not for general use)

Subroutine ENDSCP allows to shift by a user-defined distance the energy which is being deposited along a step or several step binning portions, by providing new segment endpoints​. A typical application is to simulate an instrument drift​.

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

FLDSCP: FLuence DiStributed — Change of Positions

​ ​​

Summary

Argument list

IJ :

particle type (input only)

PLA :

particle momentum (if > 0), or kinetic energy (if < 0) (input only)

TXX, TYY, TZZ  : particle direction cosines, can be modified by user
NTRUCK :

number of step points (input only)

XTRUCK, YTRUCK, ZTRUCK  : particle step points, can be modified by user
TTRUCK :

step length; should be modified to remain consistent with XTRUCK, YTRUCK, ZTRUCK.

NREG :

new region number (input only)

IOLREG :

old region number (input only)

LLO :

particle generation​ (input only)

ICALL :

internal code calling flag (not for general use)

Subroutine FLDSCP allows to shift by a user-defined distance the track whose length is being scored as fluence along a step or several step binning portions, by providing new segment endpoints​. A typical application is to simulate an instrument drift​.

If name-based input is being used, the names corresponding to MREG and IOLREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

FLUSCW: weighting fluence, current and yield

​ ​​

Summary

Argument list:

IJ :

particle type (input only, cannot be modified)

PLA :

particle momentum (if > 0.0)

or   -PLA = kinetic energy (if < 0.0)

TXX, TYY, TZZ  :  particle current direction cosines
WEE :

particle weight

XX, YY, ZZ  :  particle position
NRGFLK :

current region (after boundary crossing)

IOLREG :

previous region (before boundary crossing). Useful only with boundary crossing estimators (for other estimators it has no meaning)

LLO :

particle generation (input only, cannot be modified)​

NSURF :

internal code calling flag (not for general use)

Function FLUSCW is activated by option USERWEIG, with WHAT(3)> 0.0. Yields obtained via USRYIELD, fluences calculated with USRBDX, USRTRACK, USRCOLL, USRBIN (respectively,,,), and currents calculated with USRBDX are multiplied by the value returned by this function. The user can implement any desired logic to differentiate the returned value according to any information contained in the argument list (particle type, energy, direction, weight, position, region, boundary, particle generation), or information available in COMMON SCOHLP​ (binning or detector number, estimator type). The estimator type is given by the flag ISCRNG​ (in COMMON SCOHLP):

ISCRNG =

1 \(\longrightarrow\) Boundary crossing estimator

ISCRNG =

2 \(\longrightarrow\) Track-length binning

ISCRNG =

3 \(\longrightarrow\) Track-length estimator

ISCRNG =

4 \(\longrightarrow\) Collision density estimator

ISCRNG =

5 \(\longrightarrow\) Yield estimator

The binning/detector number is given by JSCRNG (in COMMON SCOHLP) and is printed in output:

Bdrx n. 2  "bdxname" , generalised particle n. 8, from region n. 22 to region n. 78
Track n. 6  "trkname" , generalised particle n. 14, region n.   9
Note that a track-length detector​ can have the same JSCRNG​ number as a boundary crossing one or a binning​ etc. (use the value of ISCRNG​ to discriminate the different estimators). Further information can be obtained including COMMON TRACKR​ (for instance particle age​). TRACKR contains also special user variables (both integer and in double precision) which can be used to save information about particles which have undergone some particular event.

Function FLUSCW has many applications. A common one is conditional scoring​ (score only if within a certain distance from a point, etc.): for instance it is possible to implement a sort of 2-dimensional fluence binning on a plane boundary.

Other interesting applications are based on the fact that FLUSCW is called at every boundary crossing, provided that at least one USRBDX detector has been requested. Although the function has been designed mainly to weight scored quantities, it can be “cheated” to do all sorts of side things, even not directly connected with scoring. Note that the variables in the argument list, with the exception of IJ, LLO and NSURF, are local copies of those used for particle transport, and therefore can be modified to have an effect on scoring, without affecting transport.

If name-based input is being used, the name corresponding to NREG and IOLREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

Note: setting the variable LSCZER = .TRUE. before RETURN (LSCZER is in COMMON SCOHLP), will cause zero scoring whatever the value returned by COMSCW or FLUSCW. This is more efficient than returning a zero value.

For heavy-ion fluence, current or yield scoring, all ions carry the same particle identifier IJ = -2. To distinguish them in FLUSCW, use

CALL USRDCI ( IJ, IONA, IONZ, IONM )

The returned integers are the mass number IONA, atomic number IONZ, and isomeric-state flag IONM. They can be used to accept or reject the score.

FORMFU: nuclear FORM Factor User-defined

​ ​

Summary

Argument list (all variables are input only)

IJ :

particle code, except that it is set to 3 for both \(e^+\) and \(e^-\)

QU2 :

squared momentum transfer​ (GeV/c)\(^2\)

ZMEDIU :

atomic number​ of target nucleus

AMEDIU :

atomic mass​ of target nucleus

Function FORMFU can be used to override the standard value of the nuclear charge form factor​. It must return the squared value of the nuclear charge form factor for particle IJ.

The default version computes the form factor in Born approximation​ for a medium of given composition, using the simple expression given by Tsai [Tsa74]​, and accounts also for the contribution of incoherent scattering​.
The function is called by the multiple​ and single scattering​ routines if option MULSOPT ​ has been issued with WHAT(3) <0.0 for electrons and positrons, or WHAT(2) <0.0 for hadrons and muons. See Note 2 to option MULSOPT.

FRGHNS: material roughness (for optical photons)

​ ​

Summary

Argument list (all variables are input only)

TXX, TYY, TZZ  :  particle direction cosines​
UXSRFC, UYSRFC, UZSRFC  :  direction of the normal​ to the surface
MREG :

region from which the particle is coming

NEWREG :

region to which the particle is going

MMAT :

material of the region from which the particle is coming

MMATNW :

material of the region to which the particle is going

Function FRGHNS can be used to return a non-zero value for the roughness of a boundary between two materials, relevant for optical photon transport (default roughness is zero for all boundaries). Meaningful only if options OPT–PROP or OPT–PROD,) have been requested. See OPT–PROP and Generating and propagating optical photons for more information.

If name-based input is being used, the names corresponding to MREG and NEWREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

MUSRBR, LUSRBL, FUSRBV: user-defined quantities for special binning

These three functions are used to define 3-dimensional fluence distributions to be calculated by special user-defined binnings​ (see option USRBIN with WHAT(1)= 8.0 in the first card).
MUSRBR defines a discrete (integer) variable (by default: region number).

Argument list (all variables are input only)

IJ :

particle type

PCONTR :

particle momentum

XA, YA, ZA  :  particle position
MREG :

current region

LATCLL :

current lattice cell

ICALL :

internal code calling flag (not for general use)

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

LUSRBL defines another discrete (integer) variable (by default: lattice number)

Summary

Argument list (all variables are input only)

IJ :

particle type

PCONTR :

particle momentum

XFLK, YFLK, ZFLK  :  particle position
MREG :

current region

LATCLL :

current lattice cell

ICALL :

internal code calling flag (not for general use)

FUSRBV defines a continuous (double precision) variable (no default)

Summary

Argument list (all variables are input only)

IJ :

particle type

PCONTR :

particle momentum

XFLK, YFLK, ZFLK  :  particle position
MREG :

current region

ICALL :

internal code calling flag (not for general use)

The 3 functions are called at track-length​ events. What is scored is the particle track-length multiplied by the particle’s weight, possibly modified by a user-written FLUSCW​ user routine (see FLUSCW: weighting fluence, current and yield).

LATTIC: symmetry transformation for lattice geometry

Subroutine LATTIC is activated by one or more LATTICE cards in the geometry input (see LATTICE card). It is expected to transform coordinates and direction cosines from any lattice cell​ (defined by card LATTICE) to the reference system in which the basic structure has been described.

The user is expected to provide a transformation of coordinates and vector direction cosines from each lattice cell to the corresponding basic structure (in ENTRY LATTIC) and of direction cosines from the basic structure to each corresponding lattice cell (in ENTRY LATNOR).

Entries:

LATTIC (position and direction symmetry transformation from lattice cell to prototype​ structure)

Argument list

XB(1),

XB(2), XB(3)  :  actual physical position coordinates in IRLTGG lattice cell

WB(1),

WB(2), WB(3)  :  actual physical direction cosines in IRLTGG lattice cell

DIST :

current step length

SB(1), SB(2), SB(3)  :  transformed coordinates in prototype cell
UB(1), UB(2), UB(3)  :  transformed cosines in prototype cell
IR :

region number in prototype cell

IRLTGG :

lattice cell number

IRLT :

array containing region indices corresponding to lattice cells

IFLAG :

reserved variable

LATTIC returns the tracking point coordinates (SB) and direction cosines (UB) in the reference prototype geometrical structure, corresponding to real position/direction XB, WB in the actual cell IRLTGG (defined as input region IR by a LATTICE card, see LATTICE card).

When the lattice option is activated, the tracking proceeds in two different systems: the “real” one, and that of the basic symmetry unit. Particle positions and directions are swapped from their real values to their symmetric ones in the basic cell, to perform the physical transport in the regions and materials that form the prototype geometrical structure and back again to the real world. The correspondence between “real” and “basic” position/direction depends on the symmetry transformation​ and on the lattice cell number.

LATNOR (LATtice cell NORmal transformation from prototype structure to lattice cell)

Argument list

UN(1),

UN(2), UN(3)  :  direction cosines​ of the vector normal to the surface, in the prototype cell​ (entry values) and in the lattice cell (returned values)

IRLTNO :

present lattice cell number

ENTRY LATNOR transforms the direction cosines stored in the vector UN(3) from the system of the basic prototype unit to that of the real world in lattice cell number IRLTNO. Therefore, this cosine transformation must be the inverse of that performed on the cosines by the LATTIC entry: but while LATTIC maps vector UB to a different vector WB, LATNOR maps the UN vector to itself.

Note that if the transformation implies a rotation​, it is necessary to save first the incoming UN cosines to local variables, to avoid overwriting the vector before all transformation statements are executed.


Notes

​

  1. Different symmetry transformations can of course be implemented in the same LATTIC routine (each being activated by a different cell number or range of cell numbers).

  2. The advantage of the lattice geometry is to avoid describing in detail the geometry of repetitive multi-modular structures​. It must be realised, however, that a penalty is generally paid in computer efficiency​.

  3. Also, a region contained in the prototype cell​ and all those “mapped” to it inside lattice cells​ are treated by the program as if they were connected with “non-overlapping ORs” (see Meaning of the + - OR operators, Meaning of the + - | operators) into a single region. Therefore, any region-based scoring (options SCORE, USRTRACK, etc.) can only provide quantities averaged over the whole structure. More detailed information must be obtained by region-independent options such as USRBIN or by user-written routines (MGDRAW ​, see MGDRAW: general event interface). The USRBIN and EVENTBIN options,), with WHAT(1) = 8, can also be used to request a special binning type which activates the MUSRBR, LUSRBL, FUSRBV​ ​, ​ user routines to recover lattice information (see MUSRBR, LUSRBL, FUSRBV: user-defined quantities for special binning).

  4. A transformation between a lattice cell and a prototype region can alternatively be defined without resorting to the LATTIC user routine. In this case, the transformation is defined via a ROT-DEFIni card and the correspondence is established by giving the transformation index in the SDUM of the LATTICE card (see LATTICE card).

MAGFLD: definition of a magnetic field

Summary

Argument list

X, Y, Z  :  current position (input only)
T  :  particle current age (s) (input only)
BTX, BTY, BTZ  :  direction cosines​ of the magnetic field vector (returned)
B :

magnetic field intensity​ in tesla (returned). If the field at X,Y,Z,T is zero, a negative value implies zero field; \(|B|\) (cm) overrides the default safe distance a particle can travel before the field is non-negligible.

NREG :

current region (input only)

IDISC :

if returned = 1, the particle will be discarded

MAGFLD is activated by option MGNFIELD with WHAT(4–6) = 0.0 and is used to return intensity and direction​ of a magnetic field based on the current position, time, and region. It is called only if the current region has been flagged as having a non-zero magnetic field by option ASSIGNMAt, with WHAT(5)= 1.0 or = 3.0.

The magnetic field spatial distribution is often read and interpolated from an external field map​. Note that in any case the direction cosines must be properly normalised​ in double precision (e.g. BTX = SQRT(ONEONE - BTY**2 - BTZ**2)), even if B = 0.0.

Please read carefully the notes on option MGNFIELD.

If name-based input is being used, the name corresponding to NREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

MDSTCK: management of the stack of secondaries

Summary

Argument list

IFLAG :
type of nuclear interaction which has produced secondaries:

1:

inelastic

2:

elastic

3:

low-energy neutron

NUMSEC :

number of secondary particles produced in the interaction

MDSTCK is called after a nuclear interaction​ in which at least one secondary particle​ has been produced, before any biasing is applied, to decide which secondary will be loaded in the main stack for further transport. The properties of the secondaries are stored in the secondary stack (COMMON GENSTK)​. With MDSTCK, users can analyse those secondaries, write them to a file, or even modify the content of GENSTK (for instance applying their own biasing). In the latter case, however, it is their responsibility to make sure that energy is conserved, the various physical quantities are still consistent, etc.

MGDRAW: general event interface

Subroutine MGDRAW, activated by option USERDUMP with WHAT(1)≥ 100.0, usually writes a “collision tape”, i.e. a file where all or selected transport events are recorded. The default version (unmodified by the user) offers several possibilities, selected by WHAT(3) in USERDUMP. Details are given in Collision tape.

Additional flexibility is offered by a user entry USDRAW, interfaced with the most important physical events happening during particle transport. The user can modify of course also any other entry of this subroutine (BXDRAW called at boundary crossings, EEDRAW called at event end, MGDRAW for trajectory drawing, ENDRAW for recording of energy depositions and SODRAW for recording of source events): for instance the format of the output file can be changed, and different combinations of events can be written to file.

No information is written by default at EEDRAW and BXDRAW calls, but the entries are called for any value of WHAT(3) < 7.0 in USERDUMP (EEDRAW also for WHAT(4) ≥ 1).

But the most interesting aspect of the routine is that the six entries (all of which, if desired, can be activated at the same time by setting USERDUMP with WHAT(3)= 0.0 and WHAT(4)≥ 1.0) constitute a complete interface to the whole FLUKA transport. Therefore, MGDRAW can be used not only to write a collision tape, but to do any kind of complex analysis (for instance studying correlations​ between events).

Entries:

MGDRAW (trajectory dumping for drawing)​
​​

Summary

Argument list (all variables are input only)

ICODE :
FLUKA physical compartment originating the call

= 1:

call from subroutine KASKAD (hadrons and muons)

= 2:

call from subroutine EMFSCO (\(e^-\), \(e^+\) and photons)

= 3:

call from subroutine KASNEU (low-energy neutrons)​

= 4:

call from subroutine KASHEA (heavy ions)​

= 5:

call from subroutine KASOPH (optical photons)​

MREG :

current region

MGDRAW writes by default, for each trajectory, the following variables (contained in COMMON TRACKR):
NTRACK :

number of track segments​

MTRACK :

number of continuous energy deposition events​ along the track. Local energy deposition events, i.e. energy deposition at a point, such as that from heavy recoils, particles below threshold and low energy neutron kerma, are written instead by the entry ENDRAW (see below)

JTRACK :

type of particle

ETRACK :

total energy of the particle

WTRACK :

weight of the particle

NTRACK values of XTRACK, YTRACK, ZTRACK: end of each track segment​
MTRACK values of DTRACK: energy deposited at each deposition event ​
CTRACK :

total length of the curved path​

Other variables are available in TRACKR​ (but not written by MGDRAW unless the latter is modified by the user: particle momentum, direction cosines, cosines of the polarisation vector, age​, generation​, etc. (see a full list in the comment in the INCLUDE file).

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW below.

BXDRAW (boundary crossing dumping)​
​

Summary

Argument list (all variables are input only)

ICODE :

physical compartment originating the call, as in the MGDRAW entry

MREG :

region from which the particle is exiting

NEWREG :

region the particle is entering

XSCO, YSCO, ZSCO  :  coordinates of crossing point

BXDRAW is called at each boundary crossing (if requested by the user with USERDUMP, WHAT(3) < 7.0). There is no default output: any output must be supplied by the user.

If name-based input is being used, the names corresponding to MREG and NEWREG can be obtained via a call to routine GEOR2N​:

CALL GEOR2N (NUMREG, NAMREG, IERR)
where NUMREG (input variable) is a region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful.
Example:
.......................................
CHARACTER*8 MRGNAM, NRGNAM
.......................................
ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
CALL GEOR2N ( MREG,   MRGNAM, IERR1 )
CALL GEOR2N ( NEWREG, NRGNAM, IERR2 )
IF(IERR1 .NE. 0 .OR. IERR2 .NE. 0) STOP "Error in name conversion"
.......................................
IF(MRGNAM .EQ. "MyUpsREG" .AND. NRGNAM .EQ. "MyDwnREG") THEN
.......................................
EEDRAW (event end dumping)
​

Summary

Argument list (all variables are input only)

ICODE =

-1: event not completed

= 0: normal event termination

= 4: stack overflow

EEDRAW is called at the end of each event, or primary history, (if requested by the user with USERDUMP, WHAT(3) ≤ 0.0). There is no default output: any output must be supplied by the user.

ENDRAW (dumping of energy deposition at a point) ​
​

Summary

Argument list (all variables are input only)

ICODE :

type of event originating energy deposition

ICODE

= 1\(x\): call from subroutine KASKAD​ (hadrons and muons);

= 10: elastic interaction recoil​

= 11: inelastic interaction recoil​

= 12: stopping particle​

= 14: particle escaping​ (energy deposited in blackhole)

= 15: time kill​

ICODE

= 2\(x\): call from subroutine EMFSCO​ (electrons, positrons and photons)

= 20: local energy deposition​ (i.e. photoelectric)

= 21 or 22: particle below threshold​

= 23: particle escaping​ (energy deposited in blackhole)​

= 24: time kill​

ICODE

= 3\(x\): call from subroutine KASNEU​ (low-energy neutrons)​

= 30: target recoil​

= 31: neutron below threshold​

= 32: neutron escaping (energy deposited in blackhole)

= 33: time kill​

ICODE

= 4\(x\): call from subroutine KASHEA​ (heavy ions)​

= 40: ion escaping​ (energy deposited in blackhole)​

= 41: time kill​

= 42: delta ray stack overflow​

ICODE

= 5\(x\): call from subroutine KASOPH​ (optical photons)​

= 50: optical photon absorption​

= 51: optical photon escaping​ (energy deposited in blackhole)

= 52: time kill​

MREG :

current region

RULL :

energy amount deposited

XSCO, YSCO, ZSCO  :  point where energy is deposited

ENDRAW writes by default, for each energy deposition point:
0 :

flag identifying ENDRAW output from that of other entries

ICODE :

see argument list

JTRACK, ETRACK, WTRACK  :  see MGDRAW above. Note that for recoils, electrons, positrons and photons below threshold and kerma deposition, JTRACK can be outside the allowed particle ID range, assuming values like:

208: heavy recoil

211: electron, positron or photon below threshold

308: low energy neutron kerma

In those cases the ID of the particle originating the interaction is saved in the TRACKR variable J0TRK (which otherwise has value zero)
XSCO, YSCO, ZSCO, RULL  :  see argument list.

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

SODRAW (source particle dumping)​

Summary

Argument list


No arguments

​

SODRAW writes by default, for each source​ or beam particle:
-NCASE

(in COMMON CASLIM , with a minus sign to identify SODRAW output): number of primaries followed so far​

NPFLKA

(in COMMON FLKSTK) : stack pointer​​

NSTMAX

(in COMMON FLKSTK): highest value of the stack pointer encountered so far​

TKESUM

(in COMMON SOURCM) : total kinetic energy of the primaries of a user written source (see user subroutine SOURCE in SOURCE: user-written source), if applicable. Otherwise = 0.0​

WEIPRI

(in COMMON SUMCOU) : total weight of the primaries handled so far​

NPFLKA times:

ILOFLK:

type of source particle

(all variables in

TKEFLK + AM:

total particle energy (kinetic+mass)

COMMON FLKSTK)

WTFLK:

source particle weight

XFLK, YFLK, ZFLK:

source particle position

TXFLK, TYFLK, TZFLK:

source particle direction cosines

USDRAW (user-defined dumping)​

Argument list (all variables are input only)


ICODE :

type of event

ICODE

= 10\(x\): call from subroutine KASKAD​ (hadron and muon interactions);

= 100: elastic interaction​ secondaries

= 101: inelastic interaction​ secondaries

= 102: particle decay​ secondaries

= 103: delta ray​ generation secondaries

= 104: pair production​ secondaries

= 105: bremsstrahlung​ secondaries

= 110: radioactive decay products​

ICODE

= 20\(x\): call from subroutine EMFSCO​ (electron, positron and photon interactions)

= 208: bremsstrahlung secondaries

= 210: M​ secondaries

= 212: Bhabha​ secondaries

= 214: in-flight annihilation​ secondaries

= 215: annihilation at rest​ secondaries

= 217: pair production secondaries

= 219: Compton scattering​ secondaries

= 221: photoelectric​ secondaries

= 225: Rayleigh scattering​ secondaries

ICODE

= 30\(x\): call from subroutine KASNEU​ (low-energy neutron interactions)​

= 300: neutron interaction​ secondaries

ICODE

= 40\(x\): call from subroutine KASHEA​ (heavy ion interactions)

= 400: delta ray generation secondaries


MREG :

current region

XSCO, YSCO, ZSCO  :  interaction point

USDRAW is called after each particle interaction (if requested by the user with option USERDUMP, WHAT(4)≥ 1.0). There is no default output: any output must be supplied by the user.

Information about the secondary particles produced is available in COMMON GENSTK​, except that concerning delta rays produced by heavy ions (in which case the properties of the single electron produced are available in COMMON EMFSTK​, with index NP). Another exception is that about heavy evaporation fragments (deuterons, \(^3\)H, \(^3\)He, \(\alpha\), with JTRACK ID equal respectively to -3, -4, -5, -6) and fission/fragmentation products generated in an inelastic interaction (with JTRACK = -7 to -12), which are all stored in COMMON FHEAVY​, with index NPHEAV. To get the kinetic energy of particles with JTRACK < -6), one must subtract from their total energy (ETRACK in COMMON TRACKR) their fully stripped nuclear mass (AMNHEA in COMMON FHEAVY).

Information about the interacting particle and its trajectory can be found in COMMON TRACKR​ (see description under the MGDRAW entry above). In TRACKR there are also some spare variables at the user’s disposal: LLOUSE​ (integer), ISPUSR​ (integer array) and SPAUSR​ (double precision array). Like many other TRACKR variables, each of them has a correspondent in the particle stacks, i.e. the COMMONs from which the particles are unloaded at the beginning of their transport: FLKSTK, EMFSTK​ ​ and OPPHST​ (respectively, the stack​ of hadrons/muons, electrons/photons, and optical photons). The correspondence with TRACKR is shown below under STUPRF/STUPRE​ ​ (STUPRE, STUPRF, STUPRR: SeT User PRoperties for Emf and Fluka particles). When a particle is generated, its properties (weight, momentum, energy, coordinates etc., as well as the values of the user flags​) are loaded into one of the stacks. The user can write a STUPRF or STUPRE subroutine (see description below in STUPRE, STUPRF, STUPRR: SeT User PRoperties for Emf and Fluka particles) to change anyone of such flags just before it is saved in stack.

When a particle starts to be transported, its stack variables are copied to the corresponding TRACKR ones. Unlike the other TRACKR variables, which in general become modified during transport due to energy loss, scattering etc., the user flags keep their original value copied from stack until they are changed by the user himself (generally under the USDRAW entry).

One common application is the following: after an interaction which has produced secondaries, let USDRAW copy some properties of the interacting particle into the TRACKR​ user variables. When STUPRF is called next to load the secondaries into stack, by default it copies the TRACKR user variables to the stack ones. In this way, information about the parent can be still carried by its daughters (and possibly by further descendants). This technique is sometimes referred to as “latching”​.

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

OPHBDX: user-defined Optical PHoton BounDary-(X)crossing properties

​ ​​ ​

Summary

Argument list (all variables are input only)

OMGPHO :

angular frequency​ (\(\omega = 2 \pi \nu\)) of the photon (in s\(^{-1}\))

WVLNGT :

photon wavelength​ (in cm)

MREG :

old region number

NEWREG :

new region number

SIGANW :

absorption coefficient​ in the new region (cm\(^{-1}\))

SIGDNW :

diffusion coefficient​ in the new region (cm\(^{-1}\))

RFNDPR :

refractive index​ in the new region

VGRPNW :

group velocity​ in the new region (cm s\(^{-1}\))

LPHKLL :

if .TRUE., the photon will be absorbed on the boundary​

Subroutine OPHBDX sets the optical properties of a boundary surface. The call is activated by command OPT–PROP, with SDUM = SPEC–BDX. See OPT–PROP and Generating and propagating optical photons for more information.

If name-based input is being used, the name corresponding to MREG and NEWREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

PSHCKP: user-call when PuSHing CerenKov Photons to the stack

​

Summary

Argument list (all variables are input only):
   PX/Y/ZCR : momentum components of the photon (GeV/c)
     EPHSMP : energy of the photon (GeV)
 X/Y/ZTRKCR : coordinates of the production point
     ATRCKR : photon age at production (s)
  POX/Y/ZCR : polarization components of the photon
     WTRACK : statistical weight of the photon
     LLOUSE : user defined variable (see common (TRACKR))

Subroutine PSHCKP is for informational purposes, it is called every time a Cerenkov photon is pushed to the stack

PSHSCP: user-call when PuSHing SCintillation Photons to the stack

​

Summary

Argument list (all variables are input only):
   PX/Y/ZCR : momentum components of the photon (GeV/c)
     EPHSMP : energy of the photon (GeV)
 X/Y/ZTRKCR : coordinates of the production point
     ATRCKR : photon age at production (s)
  POX/Y/ZCR : polarization components of the photon
     WTRACK : statistical weight of the photon
     LLOUSE : user defined variable (see common (TRACKR))

Subroutine PSHSCP is for informational purposes, it is called every time a scintillation photon is pushed to the stack

QUEFFC: user-defined QUantum EFFiCiency

​ ​ ​

Summary

Argument list (all variables are input only)

WVLNGT :

photon wavelength​ (in cm)

OMGPHO :

angular frequency​ (\(\omega = 2 \pi \nu\)) of the photon (in s\(^{-1}\))

Function QUEFFC returns a user-defined quantum efficiency for an optical photon of the given wavelength or frequency.
It is activated with option OPT–PROP with SDUM = SENSITIV, by setting the 0\(^{th}\) photon sensitivity parameter to a value < -99. See OPT–PROP and Generating and propagating optical photons for more information.

RFLCTV: user-defined ReFLeCTiVity

​ ​ ​

Summary

Argument list (all variables are input only)

WVLNGT :

photon wavelength​ (in cm)

OMGPHO :

angular frequency​ (\(\omega = 2 \pi \nu\)) of the photon (in s\(^{-1}\))

MMAT :

material!index

Function RFLCTV returns a user-defined value equal to \(1-r\), where \(r\) is the reflectivity of the current material for an optical photon of the given wavelength or frequency..
It is activated by command OPT–PROP with SDUM = METAL and WHAT(3) < -99. See OPT–PROP and Generating and propagating optical photons for more information.

RFRNDX: user-defined ReFRaction iNDeX

​ ​ ​

Summary

Argument list (all variables are input only)

WVLNGT :

photon wavelength​ (in cm)

OMGPHO :

angular frequency​ (\(\omega = 2 \pi \nu\)) of the photon (in s\(^{-1}\))

MMAT :

material!index

Function RFRNDX returns a user-defined refraction index of the current material for an optical photon of the given wavelength or frequency.
It is activated by command OPT–PROP with SDUM = blank and WHAT(1) < -99. See OPT–PROP and Generating and propagating optical photons for more information.

SOEVSV: SOurce EVent SaVing

​​

Summary

Argument list

No arguments

Subroutine SOEVSV is always called after a beam particle is loaded onto stack​, but a call to SOEVSV can be inserted by the user anywhere in a user routine.

SOEVSV copies the whole COMMON FLKSTK​ to another COMMON, SOUEVT​, which can be included in other user routines. In other words, this routine is used to “take a snapshot”​ of the particle bank at a particular time for further use (interfacing to independent generators, etc.)

SOURCE: user-written source

​​

Summary

Argument list

NOMORE :

if set = 1, no more calls will occur (the run will be terminated after exhausting the primary particles​ loaded onto stack​ in the present call). The history number limit set with option START will be overridden

Subroutine SOURCE is probably the most frequently used user routine. It is activated by option SOURCE and is used to sample primary particle properties from distributions (in space, energy, time, direction or mixture of particles) too complicated to be described with the BEAM, BEAMPOS and BEAMAXES cards alone. For each phase-space variable, a value must be loaded onto COMMON FLKSTK​ (particle bank​) before returning control. These values can be read from a file​, generated by some sampling algorithm​, or just assigned​.

Reading from a file

Reading from a file is needed, for instance, when the particle data are taken from a collision file​, written by FLUKA or by another program (see Collision tape). The user must open the file with a unit number > 20.0 (unit numbers lower than 20 are reserved​), in one of the following ways:

  1. Using option OPEN, with SDUM = OLD

  2. In a user subroutine USRINI​ or USRGLO​ (see USRINI: USeR INItialisation or USRGLO: USeR GLObal settings below), with a Fortran OPEN statement. Option USRICALL (resp. USRGCALL) is needed to activate the call to the routine.

  3. With an OPEN statement in the initialisation part of subroutine SOURCE itself.

Then, a READ statement in SOURCE can be used to get the data to load in stack, for instance:

READ(21,*) IPART, X, Y, Z, COSX, COSY, COSZ, ENERGY, WEIGHT
ILOFLK (NPFLKA) = IPART
XFLK   (NPFLKA) = X
YFLK   (NPFLKA) = Y
ZFLK   (NPFLKA) = Z
TXFLK  (NPFLKA) = COSX
…etc…
(NPFLKA is the current stack index​​).

Direct assignment

Direct assignment can be done explicitly, for instance:​

PMOFLK (NPFLKA) = 305.2D0

or implicitly, leaving unmodified values input with BEAM or BEAMPOS:

PMOFLK (NPFLKA) = PBEAM

(PBEAM is the momentum value input as WHAT(1) in option BEAM). A set of direct assignments, one for each of several different stack entries, can be useful, for example, to define a series of RAYs​ through the geometry (see Use of RAY pseudo-particles):

    DO 10 I = 1, 20
       NPFLKA = NPFLKA + 1
       ILOFLK (NPFLKA) = 0               ! (0 is the RAY particle id number)
       XFLK  (NPFLKA)  = 500.D0 + DBLE(I) * 40.D0
       YFLK  (NPFLKA)  = 200.D0
       ...etc...
10  CONTINUE

Sampling from a uniform distribution

To sample from a uniform distribution, the user must use the function FLRNDM(DUMMY)​, which returns a double precision pseudo-random number​​ uniformly distributed between 0 (included) and 1 (not included). Actually, DUMMY can be any variable name. A simple example of sampling from a uniform distribution is that of a linear source​ along the Z axis, between Z = 10 and Z = 80:

Z1 = 10.D0
Z2 = 80.D0
ZFLK (NPFLKA) = 10.D0 + (Z2 - Z1) * FLRNDM(XXX)

Sampling from a generic distribution

One way to sample a value \(X\) from a generic distribution \(f(x)\) is the following.
First integrate the distribution function, analytically or numerically, and normalise​ to 1 the obtained cumulative distribution​:
\[F(x)~=~\frac{\int_{x_{min}}^x\,f(x)\mathrm{d} x}{\int_{x_{min}}^{x_{max}}\,f(x)\mathrm{d} x}\]

Then, sample a uniform pseudo-random number \(\xi\) using FLRNDM and get the desired result by finding the inverse value \(X\,=\,F^{-1}(\xi)\) (analytically or most often by interpolation).

A FLUKA subroutine is available to sample directly from a Gaussian distribution​​:

CALL FLNRRN (RGAUSS)

or, if two independent Gaussian distributed numbers are needed:

CALL FLNRR2 (RGAUS1, RGAUS2)

(faster than calling FLNRRN twice).

Sampling from a biased distribution

The technique for sampling from a generic distribution described above can be extended to modify the probability of sampling in different parts of the interval (importance sampling)​. We replace \(f(x)\) by a weighted function \(g(x)\,=\,f(x)\,h(x)\), where \(h(x)\) is any appropriate function of \(x\) we like to choose. We normalise​ \(g(x)\) in the same way as \(f(x)\) before:

\[G(x)~=~\frac{\int_{x_{min}}^x\,g(x)\mathrm{d} x}{\int_{x_{min}}^{x_{max}}\,g(x)\mathrm{d} x} ~=~\frac{\int_{x_{min}}^x\,f(x)\mathrm{d} x\,h(x)\mathrm{d} x}{B}\]

and we need also the integral of \(f(x)\) over the whole interval:

\[A~=~\int_{x_{min}}^{x_{max}}\,f(x)\mathrm{d} x\]

All the sampling is done using the biased cumulative normalised ​ function \(G\) instead of the original unbiased \(F\): we sample a uniform pseudo-random number \(\xi\) as before, and we get the sampled value \(X\) by inverting \(G(x)\):

\[X~=~G^{-1}(\xi)\]

The particle is assigned a weight \({ \frac{B}{A\,h(X)}}\).

A special case of importance sampling is when the biasing function chosen is the inverse of the unbiased distribution function​:
\[h(x)~=~\frac{1}{f(x)}\]
\[g(x)~=~f(x)\,h(x)~=~1\]
\[B~=~\int_{x_{min}}^{x_{max}}\,g(x)\mathrm{d} x~=~\int_{x_{min}}^{x_{max}}\,\mathrm{d} x~=~x_{max}\,-\,x_{min}\]
\[G(x)~=~\frac{x\,-\,x_{min}}{x_{max}\,-\,x_{min}}\]

In this case we sample a uniform pseudo-random number t using FLRNDM as shown above. The sampled value \(X\) is simply given by:

\[X~=~x_{min}\,+\,(x_{max}\,-\,x_{min})\,t\]

and the particle is assigned a weight

\[\frac{B}{A\,h(X)}~=~f(X)\frac{x_{max}\,-\,x_{min}}{\int_{x_{min}}^{x_{max}}\,f(x)\mathrm{d} x}\]

But since FLUKA normalizes all results per unit primary weight, any constant factor is eliminated in the normalization. Therefore it is sufficient to assign each particle a weight \(f(x)\).

Because \(X\) is sampled with the same probability over all possible values of \(x\), independently of the value \(f(X)\) of the function, this technique is used to ensure that sampling is done uniformly over the whole interval, even though \(f(x)\) might have very small values somewhere. For instance it may be important to avoid undersampling​ in the high-energy tail of a spectrum​, steeply falling with energy but more penetrating, such as that of cosmic rays​ or synchrotron radiation​.

Option SOURCE allows the user to input up to 12 numerical values (WHASOU(1),(2)…(12)) and one 8-character string (SDUSOU) which can be accessed by the subroutine by including the following line:​

INCLUDE '(SOURCM)'

These values can be used as parameters or switches for a multi-source​ routine capable to handle several cases, or to identify an external file to be read, etc., without having to compile and link again the routine.

In the SOURCE​ routine there are a number of mandatory statements, (clearly marked as such in accompanying comments) which must not be removed or modified. The following IF block initialises​ the total kinetic energy of the primary particles and sets two flags: the first to skip the IF block in all next calls, and the second to remind the program, when writing the final output, that a user source has been used:​​

*  +-------------------------------------------------------------------*
*  |  First call initialisations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialisation ***
      END IF
*  |
*  +-------------------------------------------------------------------*

The user can insert into the above IF block any other initialisation needed, for instance the preparation of a cumulative spectrum array from which to sample the energy of the source particles. Note that user initialisation can take place also in routines USRINI​ and USRGLO​ (activated at input time by input options USRICALL and USRGCALL, see USRINI: USeR INItialisation and USRGLO: USeR GLObal settings, and USREIN​ (called before unloading from stack​ the first source particle of an event​, i.e. just after the call to SOURCE: see USREIN: USeR Event INitialisation (called before each event))

At the time SOURCE is called, the particle bank FLKSTK​ is always empty and the stack pointer​ NPFLKA​ has value 0.

The user can load onto the FLKSTK stack one or more source particles at each call: for each particle loaded the pointer must be increased by 1. The template version of SOURCE​ loads only one particle: if several are loaded the following sequence, until the statement CALL SOEVSV​ not included, must be repeated once for each particle, possibly inside a DO loop:

NPFLKA = NPFLKA + 1           !  increases the pointer

The following statements assign a value to each of the FLKSTK stack variables concerning the particle being loaded.

WTFLK (NPFLKA) = ONEONE
sets the weight of the particle​ = 1.0
This must be changed if the sampling of one or more of the particle properties are biased​. In that case, generally the weight must be set after the sampling, and its value depends on the sampling outcome.
WEIPRI = WEIPRI + WTFLK (NPFLKA)

updates the total weight of the primaries​ (don’t change)

ILOFLK (NPFLKA) = IJBEAM

by default sets the type of particle equal to the one defined by the BEAM card. If no BEAM card is given in input, IJBEAM is = 1 (proton), but it is strongly recommended to always provide a BEAM command in input (see Note 1) at the end of SOURCE: user-written source).

The above statement is followed by several others that must not be changed or removed. In the template routine, they are encompassed by the comment lines:

From this point .... / ... to this point: don't change anything

These statements are:

 * From this point ....
       LOFLK (NPFLKA) = 1             !  Generation is 1 for source particles
       LOUSE (NPFLKA) = 0             !  User variables: the user can set
       DO 100 ISPR = 1, MKBMX1        !  different values in the STUPRF or
          SPAREK (1,NPFLKA) = ZERZER  !  STUPRE routine, but it is better
100    CONTINUE
       DO 200 ISPR = 1, MKBMX2
          ISPARK (ISPR,NPFLKA) = 0
200    CONTINUE
       NPARMA = NPARMA + 1
       NUMPAR (NPFLKA) = NPARMA
       NEVENT (NPFLKA) = 0
       DFNEAR (NPFLKA) = +ZERZER

The following statements can be overridden or rewritten by the user, assigning new values or sampling them from problem-dependent distributions.

First three statements which are rarely modified:

AGESTK (NPFLKA) = +ZERZER       !  Particle age is zero by default
AKNSHR (NPFLKA) = -TWOTWO       !  Resets the Kshort component of
                                !  K0/K0bar. Usually not to be changed.
IGROUP (NPFLKA) = 0             !  Group number for low-energy
                                !  neutrons: if set to 0, the program
                                !  derives it from the kinetic energy

​​​

Then the most frequently changed lines: both energy and momentum of the particle must be loaded onto the FLKSTK​ stack​, but the two cannot be defined independently. Appropriate kinematical (relativistic) relations must be applied to derive one from the other.

In the template routine, the momentum is assumed to be assigned by BEAM option (its value, PBEAM, is taken from COMMON BEAMCM ​, which contains all values defined by options BEAM and BEAMPOS).

PMOFLK (NPFLKA) = PBEAM

Therefore, the kinetic energy (in GeV) must be derived:

TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IJBEAM)**2 ) - AM (IJBEAM)

(where AM is the rest mass, in COMMON PAPROP​, and IJBEAM is the particle type, in COMMON BEAMCM​).

If instead the energy had been sampled first from some spectrum, and ENSAMP would be the sampled value, the two statements above would become:

TKEFLK (NPFLKA) = ENSAMP
PMOFLK (NPFLKA) = SQRT(ENSAMP * (ENSAMP + TWOTWO * AM(IJBEAM)))

The direction cosines​ are loaded next:

TXFLK (NPFLKA) = UBEAM             !  Assumed here to be the same as
TYFLK (NPFLKA) = VBEAM             !  defined by option BEAMPOS. UBEAM,
TZFLK (NPFLKA) = WBEAM             !  VBEAM, WBEAM are some among the beam
                                   !  properties in COMMON BEAMCM

(If BEAMPOS is not given, by default UBEAM = VBEAM = 0.0, WBEAM = 1.0) ​

Remember to make sure that the cosines are normalised​! One could replace the last statement by:

TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK(NPFLKA)**2 - TYFLK(NPFLKA)**2 )

The polarisation cosines​ are not set by default:

TXPOL (NPFLKA) = -TWOTWO        ! -2 is a flag for "no polarisation"
TYPOL (NPFLKA) = +ZERZER
TZPOL (NPFLKA) = +ZERZER

but appropriate values need to be given in some cases, for instance in synchrotron radiation​ shielding problems.

Finally the particle coordinates​, set again by default equal to those input with BEAMPOS:

XFLK (NPFLKA) = XBEAM             !  Assumed here to be the same as
YFLK (NPFLKA) = YBEAM             !  defined by option BEAMPOS. XBEAM,
ZFLK (NPFLKA) = ZBEAM             !  YBEAM, ZBEAM are also in COMMON BEAMCM

(If BEAMPOS is not given, by default XBEAM = YBEAM = ZBEAM = 0.0).

If for example our problem required instead a linear source uniformly​ distributed along Z between Z1 and Z2, we could replace the last statement by:

ZFLK (NPFLKA) = Z1 + FLRNDM(UGH) * (Z2 - Z1)
The following lines in the template SOURCE​ routine should never be changed. They calculate the total energy of the primary particles​, define the remaining properties of the particles (starting region and lattice cell) and do some geometry initialisation​.
The last line calls the SOEVSV​ user routine (see description in SOEVSV: SOurce EVent SaVing above) to save the stack for possible further use.

Important Notes

  1. Even though a user-written source is used, it is recommended to issue in input a card BEAM with an energy larger than the maximum energy that can be sampled by the user source. This value is used to set up tables such as stopping powers, cross sections etc., and if not provided, crashes can occur.

  2. The values of beam characteristics defined by commands BEAM and POLARIZAti are available in COMMON BEAMCM​: the angular divergence (variable DIVBM), beam width (XSPOT and YSPOT), and the polarisation vector (UBMPOL, VBMPOL, WBMPOL) can help to set up a scheme to sample the corresponding quantities from user-defined distributions. But sampling from the distributions pre-defined by BEAM and POLARIZAti is not simply inherited by subroutine SOURCE: it is the responsibility of the user to write such a scheme!
    For this task, it may be useful to define a “beam reference frame”​ by means of option BEAMAXES (see more details on).

SPHSPC: user-defined Sampling optical PHoton SPectrum

​

Summary

Argument list (all variables are input only):
     IP     : ip-th  scintillation emission for material Mmat
     MMAT   : material index

Function SPHSPC returns the energy (GeV) of a scintillation optical photon sampled from an user-defined spectrum, corresponding to the ip-th scintillation emission for material Mmat. It is activated by OPT-PROD, with SDUM = SCINTILL/SCINT-WVb/SCINT-OM, and WHAT(2) < -99. See option OPT-PROD and Generating and propagating optical photons and usfsci.f for more information.

STUPRE, STUPRF, STUPRR: SeT User PRoperties for Emf and Fluka particles

These three functions are used to assign a value to one or more stack​ user variables when the corresponding particle is loaded onto one of the stacks (FLKSTK​ for hadrons/muons, EMFSTK​ for electrons/photons, OPPHST for optical photons, and RDPSTK for residual nuclei).

In each of these stacks the user has access to one integer variable, one integer array and one double precision array. Each of them is copied to a correspondent variable or array in COMMON TRACKR​ at the beginning of transport:

Correspondence

FLKSTK

EMFSTK

TRACKR

integer variable:

LOUSE

LOUEMF

\(\longrightarrow\)

LLOUSE

integer array:

ISPARK

IESPAK

\(\longrightarrow\)

ISPUSR

double precision array:

SPAREK

ESPARK

\(\longrightarrow\)

SPAUSR

The user can access and modify the TRACKR variables via subroutine MGDRAW​ and its entries ENDRAW​, SODRAW​ and especially USDRAW​ (see description above). STUPRF, STUPRE and STUPRR can be used to do the reverse, namely to copy TRACKR user variables to those of the relevant stack (see USDRAW above).

Note that a stack OPPHST​ exists also for optical photons​, containing similar user variables and arrays LOUOPP​, ISPORK​ and SPAROK​. These optical-photon user properties are handled by STUPRF.

STUPRE is called before loading into stack electrons, positrons and photons.

Summary

Argument list

No arguments

The default version does nothing (the user variables of the parent particle are already set equal to the original projectile by the various electromagnetic interaction routines. Also the region/position etc. are already set inside the stack arrays.

STUPRF is called before loading into stack hadrons, muons, neutrinos, low-energy neutrons, heavy ions and optical photons

Argument list


IJ :

type of the parent particle

MREG :

current region

XX, YY, ZZ  :  particle position
NPSECN :

index in COMMON GENSTK of the secondary being loaded onto stack​

NPPRMR :

if > 0, the secondary being loaded is actually still the interacting particle (it can happen in some biasing situations)

All heavy ions carry the same identifier IJ = -2. The characteristics of primary ions are characterised by option HI-PROPE. To obtain the those of a secondary ion, call the routine USRDCI as follows: CALL USRDCI(IJ,IONA,IONZ,IONM) The three integer values returned are the following ion properties: IONA = mass number of the ion IONZ = atomic number IONM = flag for isomeric state

The default version copies to stack the user flags of the parent.

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

UBSSET: User BiaSing SETting

Summary

Argument list


IR :

region number

RRHADR :

multiplicity biasing​ factor to be applied to the secondaries from hadronic interactions in region IR (WHAT(2) of card BIASING)

HMPHAD :

Importance of region​ IR for hadrons and muons (WHAT(3) of card BIASING, with WHAT(1)= 0.0 or 1.0). Actually the routine argument is an integer, IMPHAD, equal to importance multiplied by 10000, but the user should consider only the double precision version HMPHAD (a conversion from and to the integer version is provided at the beginning and at the end of the routine, and should not be changed)

HMPLOW :

Importance of region IR for low-energy neutrons (WHAT(3) of card BIASING, with WHAT(1)= 0.0 or 3.0). Actually the routine argument is an integer, IMPLOW, equal to importance multiplied by 10000, but the user should consider only the double precision version HMPLOW (a conversion from and to the integer version is provided at the beginning and at the end of the routine, and should not be changed)

HMPEMF :

Importance of region IR for electrons and photons (WHAT(3) of card BIASING, with WHAT(1)= 0.0 or 2.0). Actually the routine argument is an integer, IMPEMF, equal to importance multiplied by 10000, but the user should consider only the double precision version HMPEMF (a conversion from and to the integer version is provided at the beginning and at the end of the routine, and should not be changed)

IGCUTO :

Cutoff group​ index for low-energy neutrons in region IR (WHAT(1) in card LOW–BIAS)

IGNONA :

Non-analogue absorption​ group limit for low-energy neutrons in region IR (WHAT(2) in card LOW–BIAS)

PNONAN :

Non-analogue survival probability​ for low-energy neutrons in region IR (WHAT(3) in card LOW–BIAS)

IGDWSC :

Group limit for biased downscattering​ for low-energy neutrons in region IR (WHAT(1) in card LOW–DOWN)

FDOWSC :

Biased downscattering factor​ for low-energy neutrons in region IR (WHAT(2) in card LOW–DOWN)

JWSHPP :

Weight Window/importance profile​ index for low-energy neutrons in region IR (SDUM in WW–FACTOr)

WWLOW :

Weight Window lower level​ in region IR (WHAT(1) in card WW–FACTOr, possibly modified by WHAT(4) in WW–THRESh or WHAT(2) in WW–PROFIle)

WWHIG :

Weight-Window upper level in region IR (WHAT(2) in card WW–FACTOr, possibly modified by WHAT(4) in WW–THRESh or WHAT(2) in WW–PROFIle)

WWMUL :

Weight-Window multiplicative factor applied to the two energy thresholds defined with WW–THRESh, for region IR (WHAT(3) in card WW–FACTOr)

EXPTR :

Exponential transform parameter for region IR (WHAT(2) in card EXPTRANS) (not implemented yet!!!!!!!!)

ELECUT :

\(e^+\), \(e^-\) cutoff in region IR (WHAT(1) in card EMFCUT)

GAMCUT :

Photon cutoff in region IR (WHAT(2) in card EMFCUT)

LPEMF :

Leading Particle Biasing flag in region IR (SDUM = LPBEMF in card EMF–BIAS, or WHAT(3) in card EMFCUT)

ELPEMF :

Maximum \(e^+\)/\(e^-\) energy for applying Leading Particle Biasing (WHAT(2) in card EMF–BIAS with SDUM = LPBEMF)

PLPEMF :

Maximum photon energy for applying leading particle biasing (WHAT(3) in card EMF–BIAS with SDUM = LPBEMF)

Subroutine UBSSET does not require a special command to be activated: is always called several times for each region: (once for every biasing option or suboption) after the end of input reading and before starting the calculations. The default version is a dummy and does nothing. The user can replace it to override any biasing parameters specified in input.

The UBSSET subroutine is used especially in cases with a large number of regions, because it allows to derive the biasing parameters from simple algorithms instead of entering each input value by hand. Choosing an appropriate numbering scheme for the geometry regions can often facilitate the task.

For instance, assuming a simple slab geometry with an expected exponential hadron attenuation from region 3 to region 20, each region being one half-value-layer thick, one could write the following in order to set importances that would keep the hadron number about constant in all regions:

IF(IR .GE. 3 .AND. IR .LE. 20) HMPHAD = ONEONE * TWOTWO**(IR-3)
It is important, however, not to do recursive assignments of the type:
GAMCUT(5) = GAMCUT(5) * HLFHLF
(HLFHLF is a FLUKA constant = 0.5D0) because that would halve the value of the photon cutoff for Region 5 at every call, and the number of calls is not known to the user.

If name-based input is being used, the name corresponding to IR can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

UDCDRL: User defined DeCay DiRection biasing and Lambda (for \(\nu\) only)

Summary

Argument list

input:
IJ :

type of decaying particle​

KPB :

outgoing neutrino, the one direction biasing is asked for

NDCY :

number of decay products

output:
UDCDRB, VDCDRB, WDCDRB  :  cosines of the preferential outgoing direction for the neutrino
UDCDRL width of the distribution around the preferential direction ​

Function UDCDRL is used to bias the direction of a neutrino emitted by a decaying particle of type IJ event by event. The preferential direction axis is returned in the UDCDRB, VDCDRB, WDCDRB variables, and the value returned by UDCDRL is the \(\lambda\) for direction biasing ​: the zenith angle​ around the selected axis is sampled according to \(\exp[(1-\cos(\theta))/\lambda\)] .

USFSCI: USer-defined Fraction for SCIntillation PHotons

​

Summary

Argument list (all variables are input only):
     IP     : ip-th  scintillation emission for material Mmat
     MMAT   : material index

Function USFSCI returns the number of scintillation photons per deposited GeV corresponding to the ip-th scintillation emission for material Mmat. It is activated by OPT-PROD, with SDUM = SCINTILL/SCINT-WVb/SCINT-OM, and WHAT(2) < -99. See option OPT-PROD and Generating and propagating optical photons, and sphspc.f for more information.

USIMBS: USer defined IMportance BiaSing

​ ​

Summary

Argument list

input:
MREG :

region at the beginning of the step​

NEWREG :

region at the end of the step

output:
FIMP :

returns the user-defined importance ratio between the position at the end and at the beginning of the step

Subroutine USIMBS is activated by card BIASING with SDUM = USER. The routine is called at every particle step. It can be used to implement any importance biasing scheme based on region number and on phase space coordinates and other information provided by COMMON TRACKR ​.
Warning: The user must balance the very effective biasing power offered by the routine with the important demand on CPU time due to the large number of calls.

If name-based input is being used, the names corresponding to MREG and NEWREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

USREIN: USeR Event INitialisation (called before each event)

​ ​

Summary

Argument list

No arguments

Subroutine USREIN is called just before the first source particle​ of an event​ is unloaded from stack​ and begins to be transported. An event​ is the full history of a group of related particles and their descendants. If primaries are loaded into stack by the input option BEAM, there is only one source particle per event; but there can be more if the user routine SOURCE​ is used to load particles into stack. USREIN does not need any special command to be activated, but the default version of USREIN does nothing: the user can write here any kind of initialisation​.

USREOU: USeR Event OUtput (called at the end of each event)

​​

Summary

Argument list

No arguments

Subroutine USREOU is called at the end of each event​, namely after all event primary particles​ and their descendants have been transported. (See USREIN above for a definition of an event).
USREOU does not need any special command to be activated, but the default version of USREOU does nothing: the user can write here any kind of event analysis, output, etc.

USRGLO: USeR GLObal settings

​ ​

Summary

Argument list

WHAT(1),(2),(3),(4),(5),(6)  :  user-provided numerical parameters
SDUM :

user-provided character string (8 characters)

Subroutine USRGLO is called before any other initialisation is done by the program, provided a command USRGCALL is present anywhere in the input. It can be used to do any kind of initialisation​: reading and manipulating data from one or more files, calling other private routines, etc.
The calling parameters can carry any kind of useful information or can be used as flags to choose between different possible actions to be performed before any particle transport takes place.

USRHSC: USer provided RHo (density) SCaling factors

​

Summary

Argument list:
     input:
       MREG   : current region number
       REGNAM : current region name
       MMAT   : current master material number
       MMANAM : current master material name
       SFRHIO : input dE/dx density scaling factor
       SFRHOT : input other processes density scaling factor
     output:
       SFRHIO : new dE/dx density scaling factor
       SFRHOT : new other processes density scaling factor

Subroutine USRHSC allows to modify run-time the scaling factors defined in input, and or in the Voxel phantom file .

USRINI: USeR INItialisation

​ ​

Summary

Argument list

WHAT(1),(2),(3),(4),(5),(6)  :  user-provided numerical parameters
SDUM :

user-provided character string (8 characters)

Subroutine USRINI is called every time a USRICALL card is read in input. It can be used to do any kind of initialisation​: reading and manipulating data from one or more files, calling other private routines, etc.
The calling parameters can carry any kind of useful information or can be used as flags to choose between different possible actions to be performed before any particle transport takes place.

USRMED: USeR MEDium dependent directives

​ ​

Summary

Argument list

IJ :

particle type

EKSCO :

particle kinetic energy (GeV)

PLA :

particle momentum (GeV/c)

WEE :

particle weight

MREG :

previous region number

NEWREG :

current region number

XX, YY, ZZ  :  particle position
TXX, TYY, TZZ  :  particle direction

Subroutine USRMED is activated by option MAT–PROP with SDUM = USERDIRE, for one or more materials​ indicated by the user. It is called every time a particle is going to be transported in one of the user-flagged materials​.

Two cases are possible:

  1. MREG = NEWREG: the particle is going to move from a point inside the medium. The user is allowed to change only the particle weight​. Typical application: simulating attenuation of optical photons​ in an absorbing medium​ by reducing the photon weight.

  2. MREG ≠ NEWREG: the particle is going to move from a point on a boundary between two different regions. The user may change any of the following: particle weight, current region number, direction cosines.

Typical applications:

  • simulating refraction​, by changing the direction cosines so that the particle is still inside the same region. To do this, one generally needs the direction cosines of the normal to the surface: TXNOR(NPFLKA), TYNOR(NPFLKA), TZNOR(NPFLKA) (COMMON FLKSTK must be included).

  • simulating reflection​ (albedo​) at a boundary. The direction cosines​ must be modified according to some reflection law or albedo angular distribution​, and NEWREG must be set = MREG.

In both cases the weight can also be reduced to account for surface reflectivity​ or similar (if the particle is an optical photon, the FRGHNS​ user function (FRGHNS: material roughness (for optical photons)) can be called to establish a surface roughness.​) Also, setting the weight WEE to zero is a way to kill the particle.

If name-based input is being used, the names corresponding to MREG and NEWREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

USROUT: USeR OUTput

​​

Summary

Argument list

WHAT(1),(2),(3),(4),(5),(6)  :  user-given numerical parameters
SDUM :

user-given character string (8 characters)

Subroutine USROUT is called every time a USROCALL card is read in input. It is used to print special user-written output in addition to the standard one provided by default.
The calling parameters can carry any kind of useful information or can be used as flags to choose between different possible actions to be performed after all particle transport has taken place.

USRRNC: USeR Residual NuClei

​ ​

Summary

Argument list

IZ :

Atomic number​ of the residual nucleus

IA :

Mass number​ of the residual nucleus

IS :

Isomeric state​ of the residual nucleus

X, Y, Z :  particle position
MREG :

current region

WEE :

particle weight

ICALL :

internal code calling flag (not for general use)

Subroutine USRRNC is called every time a residual nucleus​ is produced, if option USERWEIG has been requested with WHAT(5) > 0.

If name-based input is being used, the name corresponding to MREG can be obtained via a call to routine GEOR2N:​

CALL GEOR2N (NUMREG, NAMREG, IERR)

where NUMREG (input variable) is the region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0 the conversion is successful. See example in the description of BXDRAW above.

USTCKV: User STepping CerenKoV

​

Summary

Argument list (input only):
     NPROD      : Total number of Cerenkov photons produced (their features
                  can be retrieved from the last NPROD entries in common
                  OPPHST)
     MREG       : Current region
     X/Y/ZTRKCR : Final production position

USTSCN: User STepping SCiNtillation

​

Summary

Argument list (input only):
     NPROD      : Total number of scintillation photons produced (their
                  features can be retrieved from the last NPROD entries in
                  common OPPHST)
     MREG       : Current region
     X/Y/ZTRKCR : Final production position

WVLNSH: WaVe LeNgth SHifted photon production

​

Summary

Argument list:
     input:
       EKPHOT     : absorbed photon energy (GeV)
       WVLNGT     : absorbed photon wave-length (cm)
       MREG       : number of the current region
       MMAT       : number of the current rmaterial
       MXWSPH     : maximum number of secondary photons
     output:
       NWVSHP     : number of produced secondary photons
       WVSHPH (I) : > 0 -> energy    (GeV) of the i_th produced photon
                    < 0 -> wavelength (cm) of the i_th produced photon
       DWVSHP (I) : production delay  (s)  of the i_th produced photon

Subroutine WVLNSH allows to produce secondary, wave-shifted, photons every time an optical photon is absorbed