Generating and propagating optical photons¶
FLUKA can be used to generate and propagate optical photons of Cherenkov, scintillation and transition radiation light. Light generation is switched off by default and is activated and totally controlled by the user by means of data cards and user routines. This is true also for the optical properties of materials. These include the refraction index as a function of wavelength (or frequency or energy), the reflection coefficient of a given material, etc.
In this respect, the user has the responsibility of issuing the right input directives: the code does not perform any physics check on the assumptions about the light yield and the properties of material.
Optical photons (FLUKA id = -1) are treated according to the laws of geometrical optics and therefore can be reflected and refracted at boundaries between different materials. From the physics point of view, optical photons have a certain energy (sampled according to the generation parameters given by the user) and carry along their polarisation information. Cherenkov photons are produced with their expected polarisation, while scintillation photons are assumed to be unpolarised. At each reflection or refraction, polarisation is assigned or modified according to optics laws derived from Maxwell equations.
At a boundary between two materials with different refraction index, an optical photon is propagated (refracted) or reflected with a relative probability calculated according to the laws of optics.
Furthermore, optical photons can be absorbed in flight (if the user defines a non zero absorption coefficient for the material under consideration) or elastically scattered (Rayleigh scattering) if the user defines a non zero diffusion coefficient for the material under consideration).
In order to deal with optical photon problems, two specific input commands are available to the user:
OPT–PROP: to set optical properties of materials.OPT–PROD: to manage light generation.
RFRNDX: to specify a refraction index as a function of wavelength, frequency or energyRFLCTV: to specify the reflectivity of a material. This can be activated by cardOPT–PROPwithSDUM = METALandWHAT(3)< -99.OPHBDX: to set optical properties of a boundary surface. The call is activated by cardOPT–PROPwithSDUM = SPEC–BDX.FRGHNS: to set a possible degree of surface roughness , in order to have both diffusive and specular reflectivity from a given material surface.QUEFFC: to request a detailed quantum efficiency treatment. This is activated by cardOPT–PROPwithSDUM = SENSITIV, setting the 0\(^{th}\) optical photon sensitivity parameter to a value lesser than-99(WHAT(1)< -99).
All running values of optical photon tracking are contained in COMMON TRACKR , just as for the other ordinary elementary particles (see INCLUDE files, MGDRAW: general event interface).
Cherenkov transport and quantum efficiency¶
In order to use quantum efficiency (via of the QUEFFC routine) the user must input a sensitivity < -100 by means of the OPT–PROP option with SDUM = SENSITIV).
That option sets the quantum efficiency as a function of photon energy overall through the problem and it is not material/region dependent. The reason is that it is applied “a priori” at photon generation time (for obvious time saving reasons). Here below is an explanation taken directly from the code.
current particle is at a given position, in a given material
Cherenkov (or scintillation) photons are going to be produced
the photon generation probability is immediately reduced over the whole energy spectrum according to the maximum quantum efficiency of the problem. The latter, set as
WHAT(5)in cardOPT–PROD, is meaningful even if routine ``QUEFFC`` is used, see belowthe detailed efficiency, set again by the
OPT–PRODcommand withSDUM = SENSITIVor by routineQUEFFC, is used for a further reduction when the actual energy of each individual photon is selected (rejecting it against Q.E.\((\omega )\)/Q.E.\(_{max}\), where \(\omega =\) photon angular frequency)
Summarising, the yes/no detection check is done at production and not at detection: this in order to substantially cut down CPU time. If one wants all photons to be produced the sensitivity must be set = 1. Then it is still possible to apply a quantum efficiency curve at detection, by means of the user weighting routine FLUSCW (see FLUSCW: weighting fluence, current and yield) or by a user-written off-line code.
Since the quantum efficiency curve provided by OPT–PROD with SDUM = SENSITIV is applied at production and not at detection, it is not known which material the photon will eventually end up in.
Furthermore, WHAT(5) must be set anyway equal to the maximum quantum efficiency over the photon energy range under consideration. One cannot use the QUEFFC routine as a way to provide an initial screening on the produced photons, i.e. to use a “safe” initial guess for the quantum efficiency (say, for instance 20%) and then, at detection, refine it through more sophisticated curves, i.e. rejecting against the actual quantum efficiency/0.2 (this again can be done in routine FLUSCW). This makes sense of course if the user has different quantum efficiency curves for different detectors (one should use in QUEFFC the curve that maximises all of them and then refine it by rejection case by case), or if the quantum efficiency is position/angle dependent upon arrival on the photomultiplier (again one should use inside QUEFFC the quantum efficiency for the most efficient position/angle and refine by rejection at detection time.
Optical photons are absorbed in those materials where the user selected properties dictate absorption, i.e. metals or materials with a non zero absorption cross section. These absorption events can be detected in different ways. For instance:
through energy deposition by particle
-1(optical photons have always id= -1), photons usually deposit all energy in one step (since only absorption and coherent scattering are implemented). So one can check forJTRACK = -1and energy deposition (RULL) in a given region (e.g. the photo-cathode of the PMT). One can also apply an extra quantum efficiency selection, e.g. using theCOMSCW user routine.through boundary crossing of particles
-1into the given region, however this is correct if and only if absorption is set such that the photon will not survive crossing the region. Again further selections can be performed, e.g. using theFLUSCW user routine.
Several practical examples of handling optical photons are presented in Handling optical photons.
Handling optical photons¶
In order to help the user to understand how to deal with optical photons, in the following we describe two input files respectively concerning the production in Liquid Argon of Cherenkov (Input Example no. 1: Only Cherenkov light is generated) and Scintillation light (Input Example no. 2: Only Scintillation light is concerned). A specific user routine, giving the refraction index of Liquid Argon as a function of wavelength is also shown (Routine assigning a continuous Refraction Index as a function of Wavelength).
It is a very simple case, in which muons are generated inside a box filled with Liquid Argon. Notice that at present it is not yet possible to request optical photons as primary particles via the BEAM card. Therefore light must be generated starting from ordinary particles, or by a special user-written SOURCE routine, where optical photons are loaded into their dedicated stack (OPPHST) instead of that of ordinary particles (FLKSTK). An example of such SOURCE is shown in Example of SOURCE routine for optical photons.
The examples presented here consider 0.5 GeV muons in a box of \(4 \times 4 \times 4\) m\(^3\). In order to avoid unnecessary complications in the example, secondary particle production by muons is switched off. Of course this is not required in real problems.
As far as the output is concerned, the following example proposes a standard energy spectrum scoring at a boundary (option USRBDX) applied to optical photons, together with a user-specific output built via the MGDRAW user routine (see MGDRAW: general event interface), where a dump of optical photon tracking is inserted. At the end of this section (in Example of SOURCE routine for optical photons) we will propose the relevant code lines to be inserted in MGDRAW (activated by the USERDUMP card), together with an example of readout (Readout of the sample user output).
Example of SOURCE routine for optical photons¶
*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
SUBROUTINE SOURCE ( NOMORE )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1990-2009 by Alfredo Ferrari & Paola Sala *
* All Rights Reserved. *
* *
* *
* New source for FLUKA9x-FLUKA20xy: *
* *
* Created on 07 january 1990 by Alfredo Ferrari & Paola Sala *
* Infn - Milan *
* *
* Last change on 08-feb-09 by Alfredo Ferrari *
* *
* This is just an example of a possible user written source routine. *
* note that the beam card still has some meaning - in the scoring the *
* maximum momentum used in deciding the binning is taken from the *
* beam momentum. Other beam card parameters are obsolete. *
* *
* Output variables: *
* *
* Nomore = if > 0 the run will be terminated *
* *
*----------------------------------------------------------------------*
*
INCLUDE '(BEAMCM)'
INCLUDE '(FHEAVY)'
INCLUDE '(FLKSTK)'
INCLUDE '(IOIOCM)'
INCLUDE '(LTCLCM)'
INCLUDE '(PAPROP)'
INCLUDE '(SOURCM)'
INCLUDE '(SUMCOU)'
INCLUDE '(OPPHST)'
INCLUDE '(TRACKR)'
*
LOGICAL LFIRST
*
SAVE LFIRST
DATA LFIRST / .TRUE. /
*======================================================================*
* *
* BASIC VERSION *
* *
*======================================================================*
NOMORE = 0
* +-------------------------------------------------------------------*
* | First call initializations:
IF ( LFIRST ) THEN
* | *** The following 3 cards are mandatory ***
TKESUM = ZERZER
LFIRST = .FALSE.
LUSSRC = .TRUE.
* | *** User initialization ***
END IF
* |
* +-------------------------------------------------------------------*
* Push one source particle to the stack. Note that you could as well
* push many but this way we reserve a maximum amount of space in the
* stack for the secondaries to be generated
* LSTOPP is the stack counter: of course any time source is called it
* must be =0
IJBEAM = -1
LSTOPP = LSTOPP + 1
* Weight of optical photon
WTOPPH (LSTOPP) = ONEONE
WEIPRI = WEIPRI + WTOPPH (LSTOPP)
NUMOPH = NUMOPH + 1
IF ( NUMOPH .GT. 1000000000 ) THEN
MUMOPH = MUMOPH + 1
NUMOPH = NUMOPH - 1000000000
END IF
WOPTPH = WOPTPH + ONEONE
*
* Insert in POPTPH (LSTOPP) the proper energy for optical photon
*
POPTPH (LSTOPP) = 4.D-09
DONEAR (LSTOPP) = ZERZER
* Injection coordinates of optical photon
XOPTPH (LSTOPP) = XBEAM
YOPTPH (LSTOPP) = YBEAM
ZOPTPH (LSTOPP) = ZBEAM
* Initial direction cosines of optical photon
TXOPPH (LSTOPP) = UBEAM
TYOPPH (LSTOPP) = VBEAM
TZOPPH (LSTOPP) = WBEAM
* Set-up the polarization vector
TXPOPP (LSTOPP) = -TWOTWO
TYPOPP (LSTOPP) = ZERZER
TZPOPP (LSTOPP) = ZERZER
* age
AGOPPH (LSTOPP) = ZERZER
* total path
CMPOPP (LSTOPP) = ZERZER
* Particle generation
LOOPPH (LSTOPP) = 1
LOUOPP (LSTOPP) = LLOUSE
DO 2100 ISPR = 1, MKBMX1
SPAROK (ISPR,LSTOPP) = ZERZER
2100 CONTINUE
DO 2200 ISPR = 1, MKBMX2
ISPORK (ISPR,LSTOPP) = 0
2200 CONTINUE
TKESUM = TKESUM + POPTPH (LSTOPP) * WTOPPH (LSTOPP)
*
CALL GEOCRS ( TXOPPH (LSTOPP), TYOPPH (LSTOPP), TZOPPH (LSTOPP)
$ )
CALL GEOREG ( XOPTPH (LSTOPP), YOPTPH (LSTOPP), ZOPTPH (LSTOPP)
$ ,NREGOP (LSTOPP), IDISC )
* Do not change these cards:
CALL GEOHSM ( IHSPNT, 1, -11, MLATTC )
NLATOP (LSTOPP) = MLATTC
CALL SOEVSV
RETURN
*=== End of subroutine Source =========================================*
END
Routine assigning a continuous Refraction Index as a function of Wavelength¶
Notice that in this example a check is performed on the material number. In the following problems, the light will be generated on material no. 18. In order to avoid problems a FLUKA abort is generated if the routine is called by mistake for a different material.
*
*=== Rfrndx ===========================================================*
*
DOUBLE PRECISION FUNCTION RFRNDX ( WVLNGT, OMGPHO, MMAT )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* user-defined ReFRaction iNDeX: *
* *
* Created on 19 September 1998 by Alfredo Ferrari & Paola Sala *
* Infn - Milan *
* *
* Last change on 25-Oct-02 by Alfredo Ferrari *
* *
* *
*----------------------------------------------------------------------*
*
INCLUDE '(FLKMAT)'
*
* Check on the material number
*
IF ( MMAT .NE. 18 ) THEN
CALL FLABRT ( 'RFRNDX', 'MMAT IS NOT SCINTILLATOR!' )
END IF
*
WL = WVLNGT * 1.D+07
RFRNDX = ONEONE
& + 9.39373D+00*(4.15D-08/(0.000087892D+00 - WL**(-2))
& + 2.075D-07 / (0.000091012D+00 - WL**(-2))
& + 4.333D-06 / (0.00021402 D+00 - WL**(-2)))
RETURN
*=== End of function Rfrndx ===========================================*
END
Input Example no. 1: Only Cherenkov light is generated¶
Cherenkov light generation depends on the refraction index. Among the different possibilities, here we have chosen to give the refraction index by means of the user routine shown above.
The relevant data cards are commented.
The value inserted for light absorption in this example is arbitrary, while the mean free path for Rayleigh scattering is the result obtained from measurements performed in the framework of the ICARUS collaboration.
TITLE
Test of Cherenkov light production in Liquid Argon
DEFAULTS PRECISIO
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
BEAM -10.000 MUON+
BEAMPOS 0.0 0.0 190.0 NEGATIVE
DELTARAY -1.0 18.0 18.0
PAIRBREM -3.0 18.0 18.0
MUPHOTON -1.0 18.0 18.0
PHOTONUC -1.0 3.0 100.0
IONTRANS -6.0
DISCARD 27.0 28.0 43.0 44.0 5.0 6.0
GEOBEGIN COMBINAT
Test
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
* A large box for the blackhole
RPP 1 -9999999. +9999999. -9999999. +9999999. -9999999. +9999999.
* A smaller box for for liquid argon
RPP 2 -200.0 +200.0 -200.0 +200.0 -200.0 +200.0
END
*== Region Definitions =================================================
* 1) Blackhole
BL1 +1 -2
* 2) Liquid Argon
LG3 +2
END
GEOEND
* Switch off electron and photon transport
EMF EMF-OFF
*
MATERIAL 18.0 0.0 1.400 18.0 ARGON
* Select neutron cross sections at liquid argon temperature
LOW-MAT 18.0 18.0 -2.0 87.0 ARGON
*
ASSIGNMAT 1.0 1.0 500. 1.0 0.0
ASSIGNMAT 18.0 2.0 2.0
*
* Set Light production/transport properties: from 100 to 600 nm in all materials
OPT-PROP 1.000E-05 3.500E-05 6.000E-05 3.0 100.0 WV-LIMIT
* Set all materials to "metal" with 0 reflectivity:
OPT-PROP 1.0 3.0 100.0 METAL
* resets all previous properties for material n. 18 (Liquid Argon)
OPT-PROP 18.0 RESET
* switches off scintillation light production in material n. 18 (Liq. Argon)
OPT-PROD 18.0 SCIN-OFF
* defines Cherenkov production for material n. 18 (Liq. Argon)
OPT-PROD 1.100E-05 5.500E-05 18.0 CEREN-WV
* The following card restores the wave-length limits for material n. 18
OPT-PROP 1.000E-05 3.500E-05 6.000E-05 18.0 WV-LIMIT
* The following card, for material n. 18:
* a) calls the RFRNDX user routine (to define the refraction index
* vs wave-length (WHAT(1)< -99)
* b) sets to 1000 cm the mean free path for absorption.
* c) sets to 90 cm the mean free path for Rayleigh scattering.
OPT-PROP -100.0 0.001 0.01111 18.0
* The following card defines the "Sensitivity" in order to introduce the
* maximum Quantum Efficiency at generation level.
* Here 1/10 of photons is actually generated.
* Fluctuations are properly sampled
OPT-PROP 0.1 0.1 SENSITIV
SCORE 208.0 211.0 201.0 210.0
RANDOMIZ 1.0
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
USRBDX 1.0 -1.0 -55.0 2.0 1.0 Opt.Phot
USRBDX 12.0E-09 0.0 120.0 &
USERDUMP 111. 2. MGDRAW
START 10000.0
STOP
Input Example no. 2: Only Scintillation light is concerned¶
Here it is necessary to point out that, at present, FLUKA can generate scintillation lines only for monochromatic lines. A maximum number of 3 different lines is possible. The value inserted here (128 nm) is the correct one for Liquid Argon. The fraction of deposited energy going into scintillation light depends on the degree of recombination after ionisation.
Again, the value used here is a parameter justified in the framework of the ICARUS collaboration, where about 20000 photons/MeV of deposited energy have been measured for the electric field of 500 V/cm (the field used in ICARUS). A different electric field intensity will change the degree of recombination and therefore the light yield.
TITLE
Test of scintillation light production in Liquid Argon
DEFAULTS PRECISIO
BEAM -0.5000 MUON+
BEAMPOS 0.0 0.0 199.0 NEGATIVE
*
DELTARAY -1.0 18.0 18.0
PAIRBREM -3.0 18.0 18.0
MUPHOTON -1.0 18.0 18.0
PHOTONUC -1.0 3.0 100.0
IONTRANS -6.0
DISCARD 27.0 28.0 43.0 44.0 5.0 6.0
GEOBEGIN COMBINAT
Test
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
* A large box for the blackhole
RPP 1 -9999999. +9999999. -9999999. +9999999. -9999999. +9999999.
* A SMALLER BOX FOR FOR LIQUID ARGON
RPP 2 -200.0 +200.0 -200.0 +200.0 -200.0 +200.0
END
*== Region Definitions =================================================
* 1) Blackhole
BL1 +1 -2
* 2) Liquid Argon
LG3 +2
END
GEOEND
*
EMF EMF-OFF
*
MATERIAL 18.0 0.0 1.400 18.0 ARGON
LOW-MAT 18.0 18.0 -2.0 87.0 ARGON
ASSIGNMAT 1.0 1.0 500. 1.0 0.0
ASSIGNMAT 18.0 2.0 2.0
*
* Set Light production/transport properties: from 100 to 600 nm in all materials
OPT-PROP 1.000E-05 1.280E-05 6.000E-05 3.0 100.0 WV-LIMIT
* Set all materials to "metal" with 0 reflectivity:
OPT-PROP 1.0 3.0 100.0 METAL
* resets all previous properties for material n. 18 (Liquid Argon)
OPT-PROP 18.0 RESET
* switches off Cherenkov light production in material n. 18 (Liquid Argon)
OPT-PROD 18.0 CERE-OFF
* defines Scint. light production for material n. 18 (Liq. Argon). Parameters:
* a) wavelength (cm) of first scintillation line.
* b) fraction of deposited energy going into scint. light
* (in Liquid Argon ~ 2 10**4 photons/MeV)
OPT-PROD 1.280E-05 1.937E-01 18.0 SCINT-WV
* The following card restores the wave-length limits for material n. 18
OPT-PROP 1.000E-05 1.280E-05 6.000E-05 18.0 WV-LIMIT
* The following card, for material n. 18:
* a) calls the RFRNDX user routine (to define the refraction index
* vs wave-length (WHAT(1)< -99)
* b) sets to 1000 cm the mean free path for absorption.
* c) sets to 90 cm the mean free path for Rayleigh scattering
OPT-PROP -100.0 0.001 0.01111 18.0
* The following card defines the "Sensitivity" in order to introduce the
* maximum Quantum Efficiency at generation level. Here 1/10 of photons are
* actually generated.
* Fluctuations are properly sampled
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
OPT-PROP 0.1 0.1 SENSITIV
SCORE 208.0 211.0 201.0 210.0
RANDOMIZ 1.0
*
USRBDX 1.0 -1.0 -55.0 2.0 1.0 Opt.Phot
USRBDX 12.0E-09 0.0 120.0 &
USERDUMP 111. 2. MGDRAW
START 10000.
STOP
User output for optical photons from USERDUMP¶
The user can request any kind of standard FLUKA output for optical photons and also a user specific output, starting for instance from the MGDRAW user routine. Here an example follows, where a few variables are simply recorded in the output “collision tape” (dump file) at each step in the tracking only for particle id = -1 (optical photons).
It can be useful, in order to exploit the flags available in this routine, and explained in the different comments, to know that the FLUKA routine which drives the transport of optical photons is KASOPH.
The content of COMMON TRACKR can be used to take full advantage of the possibilities offered by the MGDRAW routine (see MGDRAW: general event interface).
Warning: in the present version of FLUKA, there is not yet the possibility of using the User Particle Properties for optical photons (variables SPAUSR, ISPUSR and the STUPRF user routine)
* *
*=== mgdraw ===========================================================*
* *
SUBROUTINE MGDRAW ( ICODE, MREG )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* MaGnetic field trajectory DRAWing: actually this entry manages *
* all trajectory dumping for *
* drawing *
* *
* Created on 01 march 1990 by Alfredo Ferrari *
* INFN - Milan *
* last change 05-may-06 by Alfredo Ferrari *
* INFN - Milan *
* *
*----------------------------------------------------------------------*
*
INCLUDE '(CASLIM)'
INCLUDE '(COMPUT)'
INCLUDE '(SOURCM)'
INCLUDE '(FHEAVY)'
INCLUDE '(FLKSTK)'
INCLUDE '(GENSTK)'
INCLUDE '(MGDDCM)'
INCLUDE '(PAPROP)'
INCLUDE '(QUEMGD)'
INCLUDE '(SUMCOU)'
INCLUDE '(TRACKR)'
*
DIMENSION DTQUEN ( MXTRCK, MAXQMG )
*
CHARACTER*20 FILNAM
LOGICAL LFCOPE
SAVE LFCOPE
DATA LFCOPE / .FALSE. /
*
*----------------------------------------------------------------------*
* *
* Icode = 1: call from Kaskad *
* Icode = 2: call from Emfsco *
* Icode = 3: call from Kasneu *
* Icode = 4: call from Kashea *
* Icode = 5: call from Kasoph *
* *
*----------------------------------------------------------------------*
* *
IF ( .NOT. LFCOPE ) THEN
LFCOPE = .TRUE.
IF ( KOMPUT .EQ. 2 ) THEN
FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
ELSE
FILNAM = CFDRAW
END IF
WRITE(*,*) 'TRAJECTORY OPEN!'
WRITE(*,'(A)') 'FILNAM = ',FILNAM
OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
& 'UNFORMATTED' )
END IF
C
C Write trajectories of optical photons
C
IF(JTRACK .EQ. -1) THEN
WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
& SNGL (WTRACK)
WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
& SNGL (ZTRACK (I)), I = 0, NTRACK ),
& ( SNGL (DTRACK (I)), I = 1,MTRACK ),
& SNGL (CTRACK)
WRITE(IODRAW) SNGL(CXTRCK),SNGL(CYTRCK),SNGL(CZTRCK)
ENDIF
RETURN
*
*======================================================================*
* *
* Boundary-(X)crossing DRAWing: *
* *
* Icode = 1x: call from Kaskad *
* 19: boundary crossing *
* Icode = 2x: call from Emfsco *
* 29: boundary crossing *
* Icode = 3x: call from Kasneu *
* 39: boundary crossing *
* Icode = 4x: call from Kashea *
* 49: boundary crossing *
* Icode = 5x: call from Kasoph *
* 59: boundary crossing *
* *
*======================================================================*
* *
ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
RETURN
*
*======================================================================*
* *
* Event End DRAWing: *
* *
*======================================================================*
* *
ENTRY EEDRAW ( ICODE )
RETURN
*
*======================================================================*
* *
* ENergy deposition DRAWing: *
* *
* Icode = 1x: call from Kaskad *
* 10: elastic interaction recoil *
* 11: inelastic interaction recoil *
* 12: stopping particle *
* 13: pseudo-neutron deposition *
* 14: escape *
* 15: time kill *
* Icode = 2x: call from Emfsco *
* 20: local energy deposition (i.e. photoelectric) *
* 21: below threshold, iarg=1 *
* 22: below threshold, iarg=2 *
* 23: escape *
* 24: time kill *
* Icode = 3x: call from Kasneu *
* 30: target recoil *
* 31: below threshold *
* 32: escape *
* 33: time kill *
* Icode = 4x: call from Kashea *
* 40: escape *
* 41: time kill *
* 42: delta ray stack overflow *
* Icode = 5x: call from Kasoph *
* 50: optical photon absorption *
* 51: escape *
* 52: time kill *
* *
*======================================================================*
* *
ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO )
RETURN
*
*======================================================================*
* *
* SOurce particle DRAWing: *
* *
*======================================================================*
*
ENTRY SODRAW
* |
* +-------------------------------------------------------------------*
RETURN
*
*======================================================================*
* *
* USer dependent DRAWing: *
* *
* Icode = 10x: call from Kaskad *
* 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: decay products *
* Icode = 20x: call from Emfsco *
* 208: bremsstrahlung secondaries *
* 210: Moller 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 = 30x: call from Kasneu *
* 300: interaction secondaries *
* Icode = 40x: call from Kashea *
* 400: delta ray generation secondaries *
* For all interactions secondaries are put on GENSTK common (kp=1,np) *
* but for KASHEA delta ray generation where only the secondary elec- *
* tron is present and stacked on FLKSTK common for kp=npflka *
* *
*======================================================================*
*
ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
* No output by default:
RETURN
*=== End of subroutine Mgdraw ==========================================*
END
Readout of the sample user output¶
A sample program to readout the output obtained from the previously shown MGDRAW routine is presented here below. In this example the key routine in the one called VXREAD, where some trivial output is sent to logical units 66 and 67. Of course the user must adapt such a readout program to his own needs.
PROGRAM MGREAD
CHARACTER FILE*80
*
WRITE (*,*)' Name of the binary file?'
READ (*,'(A)') FILE
OPEN ( UNIT = 33, FILE = FILE, STATUS ='OLD',
& FORM = 'UNFORMATTED' )
1000 CONTINUE
WRITE (*,*)' Event number?'
READ (*,*) NCASE
IF ( NCASE .LE. 0 ) STOP
CALL VXREAD (NCASE)
GO TO 1000
END
SUBROUTINE VXREAD (NCASE)
PARAMETER ( MXH = 2000 )
PARAMETER ( MXPR = 300 )
DIMENSION XH (MXH), YH (MXH), ZH (MXH), DH (MXH),
& EPR (MXPR), WPR (MXPR), XPR (MXPR), YPR (MXPR),
& ZPR (MXPR), TXP (MXPR), TYP (MXPR), TZP (MXPR),
& IPR (MXPR)
*
LUNSCR = 33
REWIND (LUNSCR)
*
* +-------------------------------------------------------------------*
* |
NEVT=0
DO 4000 I=1,2000000000
READ (LUNSCR,END=4100) NDUM,MDUM,JDUM,EDUM,WDUM
IF(I.EQ.1) WRITE(*,*) 'NDUM,MDUM,JDUM,EDUM,WDUM',NDUM,MDUM,JDUM
& ,EDUM,WDUM
NEVT = NEVT + 1
* | +----------------------------------------------------------------*
* | | Real tracking data:
* | +----------------------------------------------------------------*
IF ( NDUM .GT. 0 ) THEN
NTRACK=NDUM
MTRACK=MDUM
JTRACK=JDUM
ETRACK=EDUM
WTRACK=WDUM
IF(NTRACK.GT.1) WRITE(67,*) 'NTRACK = ',NTRACK
READ (LUNSCR)(XH(J),YH(J),ZH(J),J=1,NTRACK+1),
& (DH(J),J=1,MTRACK), CTRACK
READ (LUNSCR) CXX,CYY,CZZ
IF(I.EQ.1) THEN
WRITE(67,*) (XH(J),YH(J),ZH(J),J=1,NTRACK+1),
& (DH(J),J=1,MTRACK), CTRACK
WRITE(67,*) CXX,CYY,CZZ
ENDIF
DO J=1,NTRACK+1
WRITE(67,*) XH(J),YH(J),ZH(J),
& CXX,CYY,CZZ
END DO
IF ( NEVT.EQ.NCASE ) THEN
WRITE(66,*)' New step:'
WRITE(66,*)' Part.id.:',JTRACK,' Kin.En.:',ETRACK,
& ' N.of substep:', NTRACK
WRITE(66,*)' X, Y, Z, i=0, # substep'
WRITE(*,*)' New step:'
WRITE(*,*)' Part.id.:',JTRACK,' Kin.En.:',ETRACK,
& ' N.of substep:', NTRACK
WRITE(*,*)' X, Y, Z, i=0, # substep'
END IF
* | |
* | +----------------------------------------------------------------*
* | | Energy deposition data:
ELSE IF ( NDUM .EQ. 0 ) THEN
ICODE1=MDUM/10
ICODE2=MDUM-ICODE1*10
IJDEPO=JDUM
ENPART=EDUM
WDEPOS=WDUM
READ (LUNSCR) XSCO, YSCO, ZSCO, ENDEPO
IF ( NEVT.EQ.NCASE ) THEN
WRITE(66,*) ' En. dep. code n.:',MDUM
WRITE(66,*) IJDEPO,' Tot. en. proj.:', ENPART,
& ' Weight:',WDEPOS
WRITE(66,*) ' Position:',XSCO,YSCO,ZSCO,
& ' En. Dep.:',ENDEPO
END IF
* | |
* | +----------------------------------------------------------------*
* | | Source particle:
ELSE
NEVT =-NDUM
LPRIMA = MDUM
NSTMAX = JDUM
TKESUM = EDUM
WEIPRI = WDUM
READ (LUNSCR) ( IPR(J),EPR(J),WPR(J),XPR(J),YPR(J),
& ZPR(J),TXP(J),TYP(J),TZP(J),J=1,LPRIMA )
DO J = 1, LPRIMA
IF ( ABS(IPR(J)) .LT. 10000 ) THEN
LPTRUE=J
END IF
END DO
LPROJ = LPRIMA - LPTRUE
LPRIMA = LPTRUE
IF (NEVT .EQ. NCASE) THEN
WRITE(66,*)' Event #',NEVT
IF ( LPROJ .GT. 0) THEN
WRITE(66,*)
& ' Original projectile(s),n. of:',LPROJ
DO IJ = 1, LPROJ
J=LPRIMA+IJ
IPR(J) = IPR(J)/10000
WRITE(66,*) ' Part.id.:',IPR(J),' Kin.en.:',
& EPR(J),' Weight:',WPR(J)
WRITE(66,*) IPR(J),EPR(J),WPR(J)
WRITE(66,*) ' Position :', XPR(J),YPR(J),ZPR(J)
WRITE(66,*) ' Direction:', TXP(J),TYP(J),TZP(J)
END DO
END IF
WRITE(66,*)' Source particle(s), n. of:',LPRIMA
DO J = 1, LPRIMA
WRITE(66,*) ' Part.id.:',IPR(J),' Kin.en.:',
& EPR(J),' Weight:',WPR(J)
WRITE(66,*) ' Position :', XPR(J),YPR(J),ZPR(J)
WRITE(66,*) ' Direction:', TXP(J),TYP(J),TZP(J)
C WRITE(67,*) XPR(J)/1.E+05,YPR(J)/1.E+05,ZPR(J)/1.E+05
END DO
END IF
IF (NEVT.GT.NCASE) GO TO 4100
END IF
* | |
* | +----------------------------------------------------------------*
4000 CONTINUE
* |
* +-------------------------------------------------------------------*
4100 CONTINUE
RETURN
END