defines a detector to score a double-differential particle yield
around an extended or a point target
(see also USRBDX)
The full definition of the detector may require two successive cards
(the second card, identified by the character '&' in any column from
71 to 78, must be given unless the corresponding defaults are
acceptable to the user)
First card:
For SDUM = anything but BEAMDEF:
WHAT(1) = ie + ia * 100, where ie and ia indicate the two physical
quantities with respect to which the double-differential
yield is calculated.
If ie > 0, the yield will be analysed in linear intervals
with respect to the first quantity; if < 0, the yield
distribution will be binned logarithmically.
(Note that for rapidity, pseudorapidity and Feynman-x
logarithmic intervals are not available and will be forced
to linear if requested).
For the second quantity, ia, only one interval will be
considered.
|ie| or |ia| = 1 : kinetic energy in GeV
= 2 : total momentum in GeV/c
= 3 : rapidity in the lab frame
(only linear scoring available)
= 4 : rapidity in the c.m.s. frame
(only linear scoring available)
= 5 : pseudorapidity in the lab frame
(only linear scoring available)
= 6 : pseudorapidity in the c.m.s. frame
(only linear scoring available)
= 7 : Feynman-x in the lab frame (E/Ebeam)
(only linear scoring available)
= 8 : Feynman-x in the c.m.s. frame
(only linear scoring available)
= 9 : transverse momentum in GeV/c
= 10 : transverse mass in GeV
= 11 : longitudinal momentum in the lab frame (GeV/c)
= 12 : longitudinal momentum in the c.m.s. frame (GeV/c)
= 13 : total energy in GeV
= 14 : polar angle in the lab frame (#)
= 15 : polar angle in the c.m.s. frame (#)
= 16 : square transverse momentum in (GeV/c)**2
= 17 : weighted angle in the lab frame (#)
= 18 : weighted transverse momentum in GeV/c (#)
= 19 : ratio laboratory momentum/beam momentum
= 20 : transverse kinetic energy
= 21 : excitation energy
= 22 : particle charge
= 23 : particle LET
( = 24 : like 14, but with input data given
in degrees rather than in radians ) (#)
( = 25 : like 15, but with input data given
in degrees rather than in radians ) (#)
= 26 : laboratory kinetic energy/nucleon
= 27 : laboratory momentum/nucleon
# see the notes for what concerns the differential
quantity (steradian instead of angle)
WHAT(2) > 0.0: (generalised) particle type to be scored.
< -80.0 and WHAT(4) = -1.0 and WHAT(5) = -2: the (generalised)
particles of type IJ ENTERING an inelastic hadronic
interaction are scored by setting WHAT(2) = - 100 - IJ
Default = 201.0 (all particles)
WHAT(3) = logical output unit:
> 0.0 : formatted data are written on WHAT(3) unit
< 0.0 : unformatted data are written on |WHAT(3)| unit
Values of |WHAT(3)| < 21.0 should be avoided (with the
exception of +11.0).
Default = 11.0 (standard output unit)
WHAT(4) > 0.0: first region defining the boundary (upstream region)
= -1.0 and WHAT(5) = -2.0: the yield of particles
EMERGING from inelastic hadronic interactions is scored
Default = 1.0
WHAT(5) > 0.0: second region defining the boundary (downstream region)
= -2.0 and WHAT(4) = -1.0: the yield of particles
EMERGING from inelastic hadronic interactions is scored
Default = 2.0
WHAT(6) = normalisation factor (the results will be divided by WHAT(6))
SDUM = detector name (max. 10 characters)
Continuation card:
WHAT(1) = Upper limit of the scoring interval for the first quantity
Default: beam energy value
WHAT(2) = Lower limit of the scoring interval for the first quantity
Default: 0.0 if linear binning, 1.0 otherwise. Note that these
values might not be meaningful for all available quantities.
WHAT(3) = number of scoring intervals for the first quantity
Default: 50.
WHAT(4) = Upper scoring limit for the second quantity
Default: no default!
WHAT(5) = Lower scoring limit for the second quantity
Default: 0.0
WHAT(6) = ixa + 100 * ixm, where ixa indicates the kind of yield or
cross section desired and ixm the target material (if needed
in order to calculate a cross section, otherwise ixm = 0)
ixa = 1 : plain double-differential cross section
d2 sigma / d x1 d x2 where x1, x2 are the first
and second quantity
ixa = 2 : invariant cross section E d3 sigma / dp3
ixa = 3 : plain double differential yield
d2 N / d x1 d x2 where x1, x2 are the first
and second quantity
ixa = 4 : double differential yield
d2 (x2 N) / d x1 d x2 where x1, x2 are the
first and second quantity
ixa = 5 : double differential yield
d2 (x1 N) / d x1 d x2 where x1, x2 are the
first and second quantity
ixa = 6 : double differential fluence yield
1/cos(theta) d2 N / d x1 d x2 where x1, x2 are
the first and second quantity, and theta is
the angle between the particle and the normal
to the surface
ixm = : material number of the target
Default: 1.0 (plain double-differential cross section)
Note that calculating a cross section has little meaning
in case of a thick target.
For SDUM = BEAMDEF:
WHAT(1) = projectile particle index
Default = IJBEAM (beam particle)
WHAT(2) = target particle index (used by the code to define
the c.m.s. frame)
Default: 1.0 (proton)
WHAT(3) = projectile momentum
Default = PBEAM (beam momentum)
WHAT(4,5,6) = projectile direction cosines
Default = VBEAM, VBEAM, WBEAM (beam direction cosines)
Default (option USRYIELD not given): no yield estimator is defined
Notes: While option USRBDX calculates angular distributions WITH
RESPECT TO THE NORMAL to the boundary at the point of crossing,
USRYIELD's distributions are calculated WITH RESPECT TO
A FIXED DIRECTION (the beam direction, or a different direction
specified by the user with SDUM = BEAMDEF).
When scoring thick-target yields, the angle considered is that
between the direction of the particle at the point where it crosses
the target surface and the beam direction (or a different direction
specified by the user, see previous Note). The target surface is
defined as the boundary between two regions (positive values of
WHAT(4) and WHAT(5) of the first USRYIELD card.
Point-target yields, i.e. yields of particles emerging from
inelastic hadronic interactions with single nuclei (including
hadronic interactions by ions and real or virtual photons), are
scored by setting WHAT(4) = -1.0 and WHAT(5) = -2.0 in the first
USRYIELD card. As an alternative, the corresponding cross sections
can be calculated, depending on the value of WHAT(6).
In addition, if WHAT(2) in the same card is < -80.0, the
distributions of particles ENTERING the inelastic hadronic
interactions can be scored.
Differential yields (or cross sections) are scored over any
desired number of intervals for what concerns the first
quantity, but over only one interval for the second quantity.
However, the results are always expressed as second derivatives
(or third derivatives in the case of invariant cross sections),
and NOT as interval-integrated yields.
In order to obtain more intervals for the second quantity, the
user must define further USRYIELD detectors.
In the case of polar angle quantities (|ie| or |ia| = 14,15,
17,18,24,25) the differential yield is always referred to
solid angle in steradian, although input is specified in
radian or degrees.
A USRYIELD card with SDUM = BEAMDEF, if given, does not refer
to a particular detector, but modifies the reference projectile
or target parameters for all USRYIELD detectors of the current
run. No continuation card has to be given after one with
SDUM = BEAMDEF.
The logical output unit for the estimator results (WHAT(3) of
the first USRYIELD card) can be any one of the following:
- the standard output unit 11: estimator results will be
written on the same file as the standard FLUKA output
- a pre-connected unit (via a symbolic link on most UNIX systems,
ASSIGN under VMS, or equivalent commands on other systems)
- a file opened with the FLUKA command OPEN
- a file opened with a Fortran OPEN statement in a user-written
initialisation routine such as USRINI or SOURCE (see 13})
- a dynamically opened file, with a default name assigned by the
Fortran compiler (typically fort.xx or ftn.xx, with xx equal
to the chosen logical output unit number).
The results of several USRYIELD detectors in a same
FLUKA run can be written on the same file, but of course only if
they are all in the same mode (all formatted, or all unformatted).
It is also possible in principle to write on the same file the
results of different kinds of estimators (USRBDX, USRBIN, etc.)
but this is not recommended, especially in the case of an
unformatted file, because it would make very difficult any reading
and analysis.
Not all 17x17 combinations of quantities are accepted by the
code, nor are they all meaningful (for instance one could run
successfully by setting WHAT(1) with ia = ie, but the result
would have no physical meaning). A list of possible quantity
combinations is given below:
An example on how to read USRYIELD unformatted output
is given below. An explanation of the meaning of the different
variables is given in the comments at the beginning of the program.
The program lists the bin boundaries of the first variable
("energy", in increasing order), and the corresponding
double-differential and cumulative yield (integrated over both
variables ("energy" and "angle").
A more complex program USYSUW, which allows to compute also
standard deviations over several runs, is available with the normal
FLUKA code distribution in directory FLUPRO/flutil.
PROGRAM RDYLD
*----------------------------------------------------------------------*
* Up to MXUSYL user-defined yield detectors are allowed *
* *
* AYLHGH = maximum angle (steradian) (WHAT(4) of 2nd card) *
* AYLLOW = minimum angle (steradian) (WHAT(5) of 2nd card) *
* DAYLBN = angle bin width
* DEYLBN = energy bin width (GeV) if linear in energy or referring *
* to a low-energy neutron energy group, otherwise *
* ratio between upper and lower edge of energy intervals *
* EYLHGH = maximum energy (GeV) (WHAT(1) of 2nd card) *
* EYLLOW = minimum energy (GeV) (WHAT(2) of 2nd card: may be *
* re-defined if low-energy neutrons are scored *
* ENGMAX = upper energies (GeV) of the neutron groups *
* IDUSYL = (generalised) particle scored (WHAT(2) of first card) *
* IGMUYL = maximum low-energy neutron group to be scored *
* IJUSYL = projectile identity *
* JTUSYL = target identity *
* ITUSYL = type of binning = WHAT(1) of first card = ie + ia*100 *
* |ie| = 1 : kinetic energy binning *
* |ie| = 2 : total momentum binning *
* |ie| = 3 : lab. rapidity binning *
* |ie| = 4 : cms rapidity binning *
* |ie| = 5 : lab. pseudorap. binning *
* |ie| = 6 : cms pseudorap. binning *
* |ie| = 7 : lab. x binning *
* |ie| = 8 : cms Feynmann x binning *
* |ie| = 9 : transverse mom. binning *
* |ie| =10 : transverse mass binning *
* |ie| =11 : lab. long. mom. binning *
* |ie| =12 : cms long. mom. binning *
* |ie| =13 : total energy binning *
* |ie| =14 : lab. angle binning *
* |ie| =15 : cms angle binning *
* |ie| =16 : p_t squared binning *
* |ie| =17 : lab. angle binning *
* with 1/2pi sin(theta) weight *
* |ie| =18 : p_t binning *
* with 1/(2pi p_t) weight *
* |ie| =19 : frac. lab mom. binning *
* |ie| =20 : trans. kin. en. binning *
* |ie| =21 : excitation en. binning *
* ie > 0 --> linear, ie < 0 --> logarithmic *
* ia has the same meaning but for the 2nd variable *
* IXUSYL = cross section kind, ixa + ixm * 100 *
* LLNUYL = no low energy neutron scoring if false, yes if true *
* MY = id-number of the yield detector *
* NCASE = number of beam particles followed *
* NEYLBN = number of energy intervals (re-defined by the program if *
* low-energy neutrons are scored) *
* NR1UYL = first region (WHAT(4) of first card) *
* NR2UYL = second region (WHAT(5) of first card) *
* PUSRYL = momentum of projectile to be used to define (possible) *
* Lorentz transformations, Feynman X etc. *
* RUNTIM = date and time of the run (as printed on standard output) *
* RUNTIT = title of the run (as given by card TITLE) *
* SCORED = result array *
* SGUSYL = adopted cross section (if any) *
* SQSUYL = cms energy for Lorentz transformation *
* TITUYL = detector name (SDUM in first USBDRX card) *
* USNRYL = user normalisation factor *
* UUSRYL,VUSRYL,WUSRYL = laboratory projectile direction *
* WEIPRI = total weight of primaries *
* IAUSYL,IEUSYL = ausiliary arrays where itusyl is decoded *
* NHIGH = number of energy bins above low-energy neutron limit *
* CUMUL = energy-angle cumulative yield *
* *
*----------------------------------------------------------------------*
PARAMETER ( MXUSYL = 100 ) ! max. number of usryield detectors
PARAMETER ( MXSCOR = 100000 ) ! storage for results
PARAMETER ( NMXGRP = 100 ) ! max # of low-energy neutron groups
PARAMETER ( PIPIPI = 3.141592653589793D+00 )
PARAMETER ( HLFHLF = 0.5D+00 )
PARAMETER ( ONEONE = 1.D+00 )
PARAMETER ( TWOTWO = 2.D+00 )
PARAMETER ( TWOPIP = 2.D+00 * PIPIPI )
LOGICAL LLNUYL
CHARACTER RUNTIT*80,RUNTIM*32,TITUYL*10,FILNAM*80
DOUBLE PRECISION CUMUL, EN1, EN2
DIMENSION EYLLOW(MXUSYL), EYLHGH(MXUSYL), AYLLOW(MXUSYL),
& AYLHGH(MXUSYL), DEYLBN(MXUSYL), NEYLBN(MXUSYL),
& NR1UYL(MXUSYL), NR2UYL(MXUSYL), USNRYL(MXUSYL),
& SGUSYL(MXUSYL), ITUSYL(MXUSYL), IAUSYL(MXUSYL),
& IDUSYL(MXUSYL), IEUSYL(MXUSYL), DAYLBN(MXUSYL),
& IGMUYL(MXUSYL), IXUSYL(MXUSYL), LLNUYL(MXUSYL),
& TITUYL(MXUSYL), ENGMAX (NMXGRP+1), SCORED(MXSCOR),
& MY(MXUSYL), NHIGH(MXUSYL)
WRITE(*,*)' Type the name of the input file:'
READ (*,'(A)') FILNAM
LQ = INDEX(FILNAM,' ') - 1
OPEN (UNIT=1, FILE=FILNAM, STATUS='OLD', FORM='UNFORMATTED')
OPEN (UNIT=2, FILE=FILNAM(1:LQ)//'.txt', STATUS='NEW')
*----------- read and write 1st record ---------------------------------
READ (1) RUNTIT,RUNTIM,WEIPRI,NCASE
WRITE(2,100) RUNTIT, RUNTIM, NCASE, WEIPRI
READ (1) IJUSYL, JTUSYL, PUSRYL, SQSUYL, UUSRYL, VUSRYL, WUSRYL
*----------- loop on bdryx detector data in the present file -----------
DO 1 IX = 1, MXUSYL
NY = IX
* ---------------- read and write 2nd record --------------------
READ (1, END=1000) MY(NY), TITUYL(NY), ITUSYL(NY), IXUSYL(NY),
& IDUSYL(NY), NR1UYL(NY), NR2UYL(NY), USNRYL(NY), SGUSYL(NY),
& LLNUYL(NY), EYLLOW(NY), EYLHGH(NY), NEYLBN(NY), DEYLBN(NY),
& AYLLOW(NY), AYLHGH(NY)
WRITE(2,101) MY(NY), TITUYL(NY), IDUSYL(NY), NR1UYL(NY),
& NR2UYL(NY), USNRYL(NY), SGUSYL(NY)
AAUSYL = 0.01D+00 * DBLE (ITUSYL(NY))
IAUSYL (NY) = NINT (AAUSYL)
IEUSYL (NY) = ITUSYL (NY) - IAUSYL (NY) * 100
* ---------------------- low-energy neutrons --------------------
* ------------ if low-en. neutrons, read group energies ---------
IF ( LLNUYL (NY) ) THEN
READ (1) IGMUYL(NY), (ENGMAX(IG), IG = 1, IGMUYL(NY)+1)
WRITE (2,102) IGMUYL(NY)
ELSE
IGMUYL(NY) = 0
END IF
* ----- if second variable is angle, differential is in sr ------
IF ( IAUSYL (NY) .EQ. 14 .OR. IAUSYL (NY) .EQ. 15 ) THEN
IF ( AYLHGH (NY) .LT. 0.3D+00 * PIPIPI ) THEN
* Small angle (exact!):
TTLOW = TAN (HLFHLF*AYLLOW(NY))
TTHGH = TAN (HLFHLF*AYLHGH(NY))
DAYLBN (NY) = TWOPIP * TWOTWO * ( TTHGH + TTLOW )
& * ( TTHGH - TTLOW ) / ( ONEONE + TTLOW**2 )
& / ( ONEONE + TTHGH**2 )
ELSE
* Large angle:
DAYLBN (NY) = TWOPIP * ( COS (DBLE(AYLLOW(NY)))
& - COS (DBLE(AYLHGH(NY))) )
END IF
ELSE
DAYLBN (NY) = AYLHGH (NY) - AYLLOW (NY)
END IF
* ------------------- linear or log in "energy" -----------------
IF ( IEUSYL (NY) .GT. 0 ) THEN
* Linear 1st variable binning
WRITE (2,103) EYLLOW(NY), EYLHGH(NY), NEYLBN(NY), DEYLBN(NY)
ELSE
* Log 1st variable binning
WRITE (2,104) EYLLOW(NY), EYLHGH(NY), NEYLBN(NY), DEYLBN(NY)
END IF
* -------------- there is 1 angle bin only ----------------------
WRITE (2,105) AYLLOW(NY), AYLHGH(NY), DAYLBN(NY)
* interv = total number of energy intervals
* (intervals above the limit for n-groups + intervals below)
INTERV = NEYLBN(NY) + IGMUYL(NY)
*------------ read the scoring results as a 1-dimensional array --------
READ (1) (SCORED(J), J = 1, INTERV)
WRITE (2,'(" ",
& " Differential Integral" )')
WRITE (2,'(" ",
& " Yield Yield " )')
WRITE (2,'("Lower energy Upper energy",
& " sr**-1 GeV**-1",3X,"number of particles"/)')
CUMUL = 0.D0
IF ( LLNUYL (NY) ) THEN
* low-energy neutron data, if present, are stored backwards
IG = IGMUYL(NY)
EN1 = ENGMAX(IG+1)
* ---------------- loop on low-energy groups -----------------
DO 2 JG = NEYLBN(NY) + IGMUYL(NY), NEYLBN(NY) + 1,-1
EN2 = ENGMAX(IG)
CUMUL = CUMUL + (EN2 - EN1) * SCORED(JG) * DAYLBN(NY)
WRITE(2,106) EN1, EN2, SCORED(JG), CUMUL
IG = IG - 1
EN1 = EN2
2 CONTINUE ! end loop on low-energy groups
* ------------------------------------------------------------
* adjust lower limit of high energies to make upper limit
* coincide with the requested one. Count the high energy bins
NHIGH(NY) = 0
EYLLOW(NY) = EN1
EN1 = EYLHGH(NY)
DO 3 IE = 1, NEYLBN(NY)
IF(IEUSYL(NY) .GT. 0) THEN
EN2 = EN1 - DEYLBN(NY)
ELSE
EN2 = EN1 / DEYLBN(NY)
END IF
EN1 = EN2
NHIGH(NY) = NHIGH(NY) + 1
IF(EN1 .LE. EYLLOW(NY)) GO TO 4
3 CONTINUE
4 CONTINUE
ELSE
EN1 = EYLLOW(NY)
NHIGH(NY) = NEYLBN(NY)
END IF
* ---------------------------------------------------------------
* first bin above the n-group limit (may have different width)
IF(IEUSYL(NY) .GT. 0) THEN
EN2 = EN1 + DEYLBN(NY)
ELSE
EN2 = EN1 * DEYLBN(NY)
END IF
* CUMUL = CUMUL + (EN2 - ETCLOW(NX)) *
CUMUL = CUMUL +
& (EN2 - EN1) * SCORED(NEYLBN(NY) - NHIGH(NY) +1) * DAYLBN(NY)
WRITE(2,106) EYLLOW(NY), EN2,
& SCORED(NEYLBN(NY) - NHIGH(NY) + 1), CUMUL
EN1 = EN2
* ----------- loop on energies above the n-group limit ----------
DO 5 IE = 2, NHIGH(NY)
IF(IEUSYL(NY) .GT. 0) THEN
EN2 = EN1 + DEYLBN(NY)
ELSE
EN2 = EN1 * DEYLBN(NY)
END IF
CUMUL = CUMUL +
& (EN2 - EN1)*SCORED(NEYLBN(NY)-NHIGH(NY)+IE) * DAYLBN(NY)
WRITE(2,106) EN1, EN2,SCORED(NEYLBN(NY)-NHIGH(NY)+IE), CUMUL
EN1 = EN2
5 CONTINUE ! end loop on energies above limit
* ---------------------------------------------------------------
1 CONTINUE ! end loop on detectors
* ------------------------------------------------------------------
1000 CONTINUE
100 FORMAT(/,1X,'*****',2X,A80,2X,'*****',/,/,10X,A32,/,/,
& 10X,'Total number of particles followed ',I9,', for a ',
& 'total weight of ',1P,E15.8,/)
101 FORMAT(/,3X,'Yield n. ',I3,' "',A10,
& '" , generalised particle n. ',I4,', from region n. ',I6,
& ' to region n. ',I6,
& /,6X,'user normalisation: ',1P,E15.8,', adopted cross ',
& 'section (if any): ',1P,E15.8,' mb')
102 FORMAT(6X,'low energy neutrons scored from group 1 to group ', I5)
103 FORMAT(6X,'linear energy binning from ',1P,E11.4,' to ',
& E11.4,' GeV, ',0P,I5,' bins (',1P,E11.4,' GeV wide)')
104 FORMAT(6X,'logar. energy binning from ',1P,E11.4,' to ',
& E11.4,' GeV, ',0P,I5,' bins (ratio :',1P,E11.4,')')
105 FORMAT(6X,'One angular bin from ',1P,E11.4,' to ',
& E11.4,' rad , (',1P,E11.4,' sr wide )'/)
106 FORMAT(1P,2(E11.4,5X),2(E14.7,5X))
END
Example:
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7...+...8
USRYIELD 1399.0 13. 21.0 3.0 2.0 1.0TotPi+(E)
USRYIELD 50.0 0.001 100.03.14159265 0.0 3.0 &
* Score double differential yield of positive pions going from region 3 to
* region 2 with a polar angle between 0 and pi with respect to the beam
* direction. Energy distribution is in 100 logarithmic intervals between 1 MeV
* and 50 GeV. Normalisation factor = 1. Results are written formatted on
* unit 21.
