------------------------------------- Argument list: NOMORE : if set = 1, no more calls will occur (the run will be terminated after exhausting the primary particles loaded into FLKSTK 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, polarisation or mixture of particles) too complicated to be described with the BEAM, BEAMPOS and POLARIZAti cards alone. For each phase-space variable, a value must be loaded into 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 fileReading from a file is needed, for instance, when the particle data are taken from a collision file, written by FLUKA or by another program. The user must open the file with a unit number > 20 (unit numbers lower than 20 are reserved), in one of the following ways: - Using option OPEN, withSDUM= OLD - In a user subroutine USRINI or USRGLO (see below), with a Fortran OPEN statement. Option USRICALL (resp. USRGCALL) is needed to activate the call to the routine. - 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 assignmentDirect 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 asWHAT(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 14}): 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 distributionTo 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 distributionOne way to sample a value XX from a generic distribution f(x) is the following. First integrate the distribution function, analytically or numerically, and normalise to 1 to obtain the normalised cumulative distribution: /x /xmax F(x) = | f(x)dx / | f(x)dx /xmin /xmin Then, sample a uniform pseudo-random number t using FLRNDM and get the desired result by finding the inverse value: -1 XX = F (t) (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 distributionThe 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: /x /xmax /x G(x) = | g(x)dx / | g(x)dx = | f(x)*h(x)dx / B /xmin /xmin /xmin and we need also the integral of f(x) over the whole interval: /xmax A = | f(x)dx /xmin 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 t as before, and we get the sampled value XX by inverting G(x): -1 XX = G (t) The particle is assigned a weight B/(A * h(XX)) A special case of importance sampling is when the biasing function chosen is the inverse of the unbiased distribution function: h(x) = 1/f(x) g(x) = f(x) * h(x) = 1 G(x) = (x - xmin) / (xmax - xmin) In this case we sample a uniform pseudo-random number t using FLRNDM as shown above. The sampled value XX is simply given by: XX = xmin + (xmax - xmin)*t and the particle is assigned a weight: /xmax B/(A h(X)) = f(XX)*(xmax - xmin) / | f(x)dx /xmin 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 XX is sampled with the same probability over all possible values of x, independently of the value f(XX) 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) and USREIN (called before unloading from stack the first source particle of an event, i.e., just after the call to SOURCE) 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 into 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 variables concerning the particle being loaded: WTFLK (NPFLKA) = ONEONE sets the weight of the particle = 1. 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 Section 13.2.19}). 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 not to do it here DO 200 ISPR = 1, MKBMX2 ISPARK (ISPR,NPFLKA) = 0 More user variables (integer) 200 CONTINUE NPARMA = NPARMA + 1 Updating the maximum particle number NUMPAR (NPFLKA) = NPARMA Setting the current particle number NEVENT (NPFLKA) = 0 Resetting the current event number DFNEAR (NPFLKA) = +ZERZER Resetting the distance to the nearest boundary * ... to this point: don't change anything 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 FLKSTK, 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 above) to save the stack for possible further use. IMPORTANT! ---------- Note 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. Note 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 at the BEAMAXES command).