                Propagating optical photons with FLUKA

1) Introduction
===============
FLUKA can be used to generate and propagate optical photons of Cerenkov,
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. This is true also for the optical properties of
materials. These include the refraction index as a function of wave-length 
(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 physical check on the assumptions
about the light yield and the properties of material.

Optical photons (FLUKA id = -1) are treated according 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 their right energy (sampled according to the 
generation parameters given by the user) and carry along their polarization
information. Cerenkov photons are produced with their expected
polarization, while scintillation photons are assumed to be unpolarized. At
each reflection or refraction polarization is assigned or modified
according to optics laws derived from Maxwell equations.

At a boundary between 2 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 within the framework of Rayleigh scattering (if the
user defines a non zero diffusion coefficient for the material under
consideration.

In order to deal with optical photon problems, the user has a set of
specific data cards to build the input file, and some user routines.

The data cards for optical photon management are of 2 types: 1) OPT-PROP:
to set optical properties of materials.  2) OPT-PROD: to manage light
generation.

The optional user routines concerning optical photons are: i) rfrndx.f: to
specify a refraction index as a function of wavelength, frequency or energy
ii) rflctv.f: to specify the reflectivity of a material) iii) ophbdx.f: to
set optical properties of a boundary surface iv) frghns.f: to set a
possible degree of surface roughness, in order to have both diffusive and
specular reflectivity from a given material surface.  v) queffc.f: to
introduce Quantum efficiency

All running values of optical photon tracking is contained in the TRACKR
common block as for the other ordinary elementary particles.

2) Summary of OPTical PRODuction cards
======================================
*----------------------------------------------------------------------*
*                                                                      *
*     ReaDing material CeRenkov and Transition Radiation properties:   *
*                                                                      *
*     Created on 21 September 1997 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 07-Nov-02     by    Alfredo Ferrari               *
*                                                                      *
*               control card: codewd = opt-prod                        *
*               allows to control/activation Cerenkov, Transition      *
*               Radiation and Scintillation production                 *
*                                                                      *
*     Documentation and usage:                                         *
*                                                                      *
* for sdum    = CERE-OFF (switch off Cerenkov production)              *
* what(1,2,3) = no meaning                                             *
* for sdum    = TRD-OFF  (switch off Transition Radiation production)  *
* what(1,2,3) = no meaning                                             *
* for sdum    = SCIN-OFF (switch off Scintillation photon production)  *
* what(1,2,3) = no meaning                                             *
* for sdum    = CERENKOV (switch on Cerenkov photon production)        *
*     what(1) = minimum Cerenkov photon emission energy                *
*     what(2) = maximum Cerenkov photon emission energy                *
*     what(3) = no meaning                                             *
* for sdum    = CEREN-WV (switch on Cerenkov photon production)        *
*     what(1) = minimum Cerenkov photon emission wavelength            *
*     what(2) = maximum Cerenkov photon emission wavelength            *
*     what(3) = no meaning                                             *
* for sdum    = CEREN-OM (switch on Cerenkov photon production)        *
*     what(1) = minimum Cerenkov photon emission 2 x pi x frequency    *
*     what(2) = maximum Cerenkov photon emission 2 x pi x frequency    *
*     what(3) = no meaning                                             *
* for sdum    = TR-RADIA (switch on Transition Radiation production)   *
*     what(1) = minimum TRD photon emission energy                     *
*     what(2) = maximum TRD photon emission energy                     *
*     what(3) = no meaning                                             *
* for sdum    = SCINTILL (switch on scintillation photon production)   *
*     what(1) = (i_th) scintillation photon emission energy            *
*     what(2) = fraction of deposited energy going into (i_th)         *
*               scintillation photon emission (-99=<what(2)=<0 ignored,*
*               what(2)=<-100 forces to use an user routine)           *
*     what(3) = time constant (s)                                      *
* for sdum    = SCINT-WV (switch on scintillation photon production)   *
*     what(1) = scintillation photon emission wavelength               *
*     what(2) = fraction of deposited energy glint into (i_th)         *
*               scintillation photon emission (-99=<what(2)=<0 ignored,*
*               what(2)=<-100 forces to use an user routine)           *
*     what(3) = time constant (s)                                      *
* for sdum    = SCINT-OM (switch on scintillation photon production)   *
*     what(1) = scintillation photon emission 2 x pi x frequency       *
*     what(2) = fraction of deposited energy going into (i_th)         *
*               scintillation photon emission (-99=<what(2)=<0 ignored,*
*               what(2)=<-100 forces to use an user routine)           *
*     what(3) = time constant (s)                                      *
* for all previous sdum's                                              *
*     what(4) = from material what (4) (def. 3) ............           *
*     what(5) = ....... to material what (5) (def. what(4)) .........  *
*     what(6) = ....... in step of  what (6) (def. 1)                  *
*                                                                      *
*----------------------------------------------------------------------*

3) Summary of OPTical PROPperty cards is the following:
=======================================================
*----------------------------------------------------------------------*
*                                                                      *
*     ReaDing material OPtical PRoperties:                             *
*                                                                      *
*     Created on 21 September 1997 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 10-Nov-02     by    Alfredo Ferrari               *
*                                                                      *
*               control card: codewd = opt-prop                        *
*               allows to control material optical properties,         *
*               as well as detection sensitivities                     *
*                                                                      *
*     Documentation and usage:                                         *
*                                                                      *
* for sdum    = WV-LIMIT                                               *
*     what(1) = minimum wavelength for opt. phot. transport            *
*     what(2) = central wavelength for opt. phot. transport            *
*     what(3) = maximum wavelength for opt. phot. transport            *
*          (all min,cen,max reset to default if < 0, untouched if = 0) *
* for sdum    = OM-LIMIT                                               *
*     what(1) = minimum 2pi x freq. (s^-1) for opt. phot. transport    *
*     what(2) = central 2pi x freq. (s^-1) for opt. phot. transport    *
*     what(3) = maximum 2pi x freq. (s^-1) for opt. phot. transport    *
*          (all min,cen,max reset to default if < 0, untouched if = 0) *
* for sdum    = RESET --> all optical properties zeroed                *
*     what(1,2,3) = no meaning                                         *
* for sdum    = METAL --> flag the material as a metal                 *
*     what(1) = 1st optical property (not used for the moment)         *
*     what(2) = 2nd optical property (not used for the moment)         *
*     what(3) = 3rd optical property (1-reflectivity index at          *
*               Wv/Omcnop)                                             *
* for sdum different from WV-LIMIT, OM-LIMIT, METAL, RESET and sdum not*
*     containing &                                                     *
*     what(1) = 1st optical property (refr. index at Wv/Omcnop)        *
*     what(2) = 2nd optical property (abs. coeff. (cm^-1) at Wv/Omcnop)*
*     what(3) = 3rd optical property (dif. coeff. (cm^-1) at Wv/Omcnop)*
* for sdum containing &1 (&2...) etc                                   *
*     what(1) = 4(7...)th optical property of the material             *
*     what(2) = 5(8...)th optical property of the material             *
*     what(3) = 6(9...)th optical property of the material             *
*     The 4th,5th,6th material properties are:                         *
*         d n_refr/ dx, d mu_abs/dx, d mu_dif/dx                       *
*     the 7th,8th,9th material properties are:                         *
*         d2 n_refr/ dx^2, d2 mu_abs/dx^2, d2 mu_dif/dx^2              *
*     the 10th,11th,12th material properties are:                      *
*         d3 n_refr/ dx^3, d3 mu_abs/dx^3, d3 mu_dif/dx^3              *
*     (the same applies to the reflectivity for metals)                *
*     where x is:                                                      *
*            x = ( lamda - lamda_cen ) / lamda_cen    or               *
*            x = ( omega - omega_cen ) / omega_cen                     *
* for all previous sdum's                                              *
*     what(4) = from material what (4) (def. 3) ............           *
*     what(5) = ....... to material what (5) (def. what(4)) .........  *
*     what(6) = ....... in step of  what (6) (def. 1)                  *
*                                                                      *
* for sdum    = WV-SENSI                                               *
*     what(1) = minimum wavelength for opt. phot. sensitivity          *
*     what(2) = central wavelength for opt. phot. sensitivity          *
*     what(3) = maximum wavelength for opt. phot. sensitivity          *
*          (all min,cen,max reset to default if < 0, untouched if = 0) *
*     what(4) = not used                                               *
*     what(5) = not used                                               *
*     what(6) = not used                                               *
* for sdum    = OM-SENSI                                               *
*     what(1) = minimum 2pi x freq. (s^-1) for opt. phot. sensitivity  *
*     what(2) = central 2pi x freq. (s^-1) for opt. phot. sensitivity  *
*     what(3) = maximum 2pi x freq. (s^-1) for opt. phot. sensitivity  *
*          (all min,cen,max reset to default if < 0, untouched if = 0) *
*     what(4) = not used                                               *
*     what(5) = not used                                               *
*     what(6) = not used                                               *
* for sdum    = SENSITIV (setup the optical photon detection sensiti-  *
*                         vity parameters)                             *
*     what(1) = 0th optical photon sensitivity parameter               *
*     what(2) = 1st optical photon sensitivity parameter               *
*     what(3) = 2nd optical photon sensitivity parameter               *
*     what(4) = 3rd optical photon sensitivity parameter               *
*     what(5) = maximum optical photon sensitivity over the allowed    *
*               range (must be consistent with the previous values,    *
*               can be overestimated, def.=1)                          *
*     what(6) = not used                                               *
*     The optical photon sensitivity parameters are:                   *
*         sens(0), d sens/dx, d2 sens/dx^2, d3 sens/dx^3               *
*     where x is:                                                      *
*            x = ( lamda - lamda_sens_cen ) / lamda_sens_cen    or     *
*            x = ( omega - omega_sens_cen ) / omega_sens_cen           *
*                                                                      *
*     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   *
*     !!! Setting the central value for a given optical property !!!   *
*     !!! to <-99 forces the use of user supplied routines       !!!   *
*     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   *
*                                                                      *
* for sdum    = SPEC-BDX (flag special boundary crossings where special*
*                         user defined properties must be used by      *
*                         calling the user routine Ophbdx)             *
*     what(1) = >= 1: special boundary treatment activated             *
*                    (n_th+1 boundary)                                 *
*               =  0: ignored                                          *
*               =<-1: special boundary treatment deactivated           *
*                    (n_th+1 boundary)                                 *
*     what(2) = From/to region Nint(what(2)) (n_th+1 boundary)         *
*     what(3) = To/from region Nint(what(3)) (n_th+1 boundary)         *
*     what(4) = >= 1: special boundary treatment activated             *
*                    (n_th+2 boundary)                                 *
*               =  0: ignored                                          *
*               =<-1: special boundary treatment deactivated           *
*                    (n_th+2 boundary)                                 *
*     what(5) = From/to region Nint(what(5)) (n_th+2 boundary)         *
*     what(6) = To/from region Nint(what(6)) (n_th+2 boundary)         *
*                                                                      *
*----------------------------------------------------------------------*

4) Examples of use
==================
In order to help the user understand how to deal with optical photons,
in the following we include 2 input files and, describing respectively the
production of Cerenkov and Scintillation light in Liquid Argon. A specific
user routine, giving the refraction index of Liquid Argon as a function of
wave-length is also shown.

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 assign OPTICAL PHOTONS
to the BEAM card, therefore light must be generated starting from ordinary
particles. The examples here enclosed consider 0.5 GeV muons in a box
of 4 x 4 x 4 m**3. In order to avoid unnecessary complications in this
example, secondary particle production by muons is switched off. This is of
course not required in real problems.

As far as the output is concerned, the following example propose a standard
binning (USRBDX) applied to optical photons, together with a user-specific
output built starting from the mgdraw.f user routine, where a dump of
optical photon tracking is inserted. At the end of this paragraph we shall 
propose the relevant code lines to be inserted in mgdraw.f (activated by the
USEMPDUMP card), together with an example of readout.

4.1) Routine to assign a continuous Refraction index of Liq. Argon as a
     function of Wave-Length
 
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 '(MAPA)'
*
* 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


4.2) Example input no. 1: Only Cerenkov light is generated

Cerenkov 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 diffusion coefficient for Rayleigh scattering is the result obtained
from measurements performed in the framework of the Icarus collaboration.

TITLE
Test of Cerenkov light production in Liquid Argon
DEFAULTS                                                              PRECISIO
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
EVENTYPE                       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...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      39.948    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 and 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 Cerenkov 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 of wave-length (WHAT(1)< -99)
* b) puts to 200 cm the absorption coefficient
* c) puts to 90 cm the diffusion (Rayleigh) coefficient
OPT-PROP   -100.0    0.001      0.01111    18.0                    
* The following card defines the "Sensitivity" in order to introduce the maximum Quantum Eff.
* 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
OUTLEVEL  1.0       4.0
*
RANDOMIZ  1.0
*
USRBDX    -1.0      213.0     -48.0     2.0        1.0       1.0      Electrons
USRBDX     0.100    0.01       30.0                           1.         &
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      1.                 10000.0   1.0        0.0
STOP


4.3) Example input 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 ionization.
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, as it is the one used in
Icarus. A different electric field intensity will change the degree of
recombination and therefore the light yield.


TITLE
Test of scint. 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
EVENTYPE                       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...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      39.948    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 and 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 Cerenkov light production in material n. 18 (Liq. 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 Liq. Argon ~ 2 10**4 ph/MeV)
OPT-PROD   1.280E-05 9.686E-02           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 of wave-length (WHAT(1)< -99)
* b) puts to 200 cm the absorption coefficient
* c) puts to 90 cm the diffusion (Rayleigh) coefficient
OPT-PROP   -100.0    0.001     0.01111    18.0                    
* The following card defines the "Sensitivity" in order to introduce the maximum Quantum Eff.
* 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
OUTLEVEL  1.0       4.0
*
RANDOMIZ  1.0
*
USRBDX    -1.0      213.0     -48.0     2.0        1.0       1.0      Electrons
USRBDX     0.100    0.01       30.0                           1.         &
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      1.                 10000.0   1.0        0.0
STOP

4.4) User output from USERDUMP

The user can build any kind of standard FLUKA output for optical photons
and also a user specific output, starting for instance from the mgdraw.f
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, that the FLUKA routine which drives the transport of optical
photons is KASOPH.

*                                                                      *
*=== 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 by                     Alfredo Ferrari                   *
*                                    INFN - Milan                      *
*     last change  10-Nov-02  by     Alfredo Ferrari                   *
*                                    INFN - Milan                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(CASLIM)'
      INCLUDE '(COMPUT)'
      INCLUDE '(EPISOR)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FINUC)'
      INCLUDE '(MGDDCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(STACK)'
      INCLUDE '(STARS)'
      INCLUDE '(TRACKR)'
*
      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     Selects 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                                            *
*     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                         *
*     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 FINUC common (kp=1,np)  *
*  but for KASHEA delta ray generation where only the secondary elec-  *
*  tron is present and stacked on STACK common for kp=lstack           *
*                                                                      *
*======================================================================*
*
      ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
* No output by default:
      RETURN
*=== End of subroutine Mgdraw ==========================================*
      END

4.5) Readout of the sample user output

A sample program to readout the output obtained from the previously
shown mgdraw.f, is here presented. 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
*
      CALL CMSPPR
      WRITE (*,*)' Binary file?'
      READ  (*,'(A)') FILE
* Start_UNIX_seq
      OPEN ( UNIT = 33, FILE = FILE, STATUS ='OLD',
     &       FORM = 'UNFORMATTED' )
* End_UNIX_seq
 1000 CONTINUE
         WRITE (*,*)' Event #:'
         READ  (*,*) NCASE
         IF ( NCASE .LE. 0 ) GO TO 9999
         CALL VXREAD (NCASE)
      GO TO 1000
 9999 CONTINUE
      STOP
      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
            LSTMAX = 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

4.6) The  TRACKR common

The content of the TRACKR common is here recalled to help users to exploit
as much as possible the possibilities offered from the use of mgdraw routine.
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 stupfr.f user routine)

*=== trackr ===========================================================*
*                                                                      *
*----------------------------------------------------------------------*
*                                                                      *
*     TRACKs Recording       by  Alfredo Ferrari, INFN - Milan         *
*                                                                      *
*     last change    31 January 2001    by   Alfredo Ferrari           *
*                                                                      *
*            included in :                                             *
*                          electr                                      *
*                          emfsco                                      *
*                          kaskad (new version)                        *
*                          kashea                                      *
*                          kasneu                                      *
*                          geoden (new version)                        *
*                          mageas                                      *
*                          magmov                                      *
*                          magnew                                      *
*                          move                                        *
*                          photon                                      *
*                          usrsco                                      *
*                                                                      *
*          Ntrack = number of track segments                           *
*          Mtrack = number of energy deposition events along the track *
*   0 < i < Ntrack                                                     *
*          Xtrack = end x-point of the ith track segment               *
*          Ytrack = end y-point of the ith track segment               *
*          Ztrack = end z-point of the ith track segment               *
*   1 < i < Ntrack                                                     *
*          Ttrack = length of the ith track segment                    *
*   1 < j < Mtrack                                                     *
*          Dtrack = energy deposition of the jth deposition event      *
*          Dptrck = momentum loss of the jth deposition event          *
*                                                                      *
*          Jtrack = identity number of the particle                    *
*          Etrack = total energy of the particle                       *
*          Ptrack = momentum of the particle (not always defined, if   *
*                 < 0 must be obtained from Etrack)                    *
*      Cx,y,ztrck = direction cosines of the current particle          *
*      Cx,y,ztrpl = polarization cosines of the current particle       *
*          Wtrack = weight of the particle                             *
*          Wscrng = scoring weight: it can differ from Wtrack if some  *
*                   biasing techniques are used (for example inelastic *
*                   interaction length biasing)                        *
*          Ctrack = total curved path                                  *
*          Cmtrck = cumulative curved path since particle birth        *
*          Zfftrk = <Z_eff> of the particle                            *
*          Zfrttk = actual Z_eff of the particle                       *
*          Atrack = age of the particle                                *
*          Akshrt = Kshrt amplitude for K0/K0bar                       *
*          Aklong = Klong amplitude for K0/K0bar                       *
*          Wninou = neutron algebraic balance of interactions (both    *
*                   for "high" energy particles and "low" energy       *
*                   neutrons)                                          *
*          Spausr = user defined spare variables for the current       *
*                   particle                                           *
*          Sttrck = macroscopic total cross section for low energy     *
*                   neutron collisions                                 *
*          Satrck = macroscopic absorption cross section for low energy*
*                   neutron collisions (it can be negative for Pnab>1) *
*          Ktrack = if > 0 neutron group of the particle (neutron)     *
*                                                                      *
*          Ntrack > 0, Mtrack > 0 : energy loss distributed along the  *
*                                   track                              *
*          Ntrack > 0, Mtrack = 0 : no energy loss along the track     *
*          Ntrack = 0, Mtrack = 0 : local energy deposition (the       *
*                                   value and the point are not re-    *
*                                   corded in Trackr)                  *
*          Mmtrck = flag recording the material index for low energy   *
*                   neutron collisions                                 *
*          Lt1trk = initial lattice cell of the current track          *
*                  (or lattice cell for a point energy deposition)     *
*          Lt2trk = final   lattice cell of the current track          *
*          Ihspnt = current geometry history pointer (not set if -1)   *
*          Ltrack = flag recording the generation number               *
*          Llouse = user defined flag for the current particle         *
*          Ispusr = user defined spare flags for the current particle  *
*          Lfsssc = logical flag for inelastic interactions ending with*
*                   fission (used also for low energy neutrons)        *
*          Lpkill = logical (user) flag for sudden particle death      *
*                                                                      *
*----------------------------------------------------------------------*
*                                                                      *
      PARAMETER ( MXTRCK = 2500 )
      LOGICAL LFSSSC, LPKILL
      COMMON / TRACKR /  XTRACK ( 0:MXTRCK ), YTRACK ( 0:MXTRCK ),
     &                   ZTRACK ( 0:MXTRCK ), TTRACK   ( MXTRCK ),
     &                   DTRACK   ( MXTRCK ), DPTRCK ( 3,MXTRCK ),
     &                   ETRACK, PTRACK, CXTRCK, CYTRCK, CZTRCK, WTRACK, 
     &                   CXTRPL, CYTRPL, CZTRPL, ZFFTRK, ZFRTTK, ATRACK,
     &                   CTRACK, CMTRCK, AKSHRT, AKLONG, WSCRNG, WNINOU,
     &                   SPAUSR(MKBMX1), STTRCK, SATRCK, NTRACK, MTRACK,
     &                   JTRACK, KTRACK, MMTRCK, LT1TRK, LT2TRK, IHSPNT,
     &                   LTRACK, LLOUSE, ISPUSR(MKBMX2), LFSSSC, LPKILL
      EQUIVALENCE ( SPAUSE, SPAUSR (1) )
      EQUIVALENCE ( ISPUSE, ISPUSR (1) )

5) Questions on Cerenkov transport and quantum efficiency
=========================================================
 
In order to use queffc, you have to input a sensitivity < -100 => the routine gets 
called (with the OPT-PROP with SDUM SENSITIV). 


>(1) When I use SENSITIV for one material, queffc gets actually called for 
>all materials for which optical photon transport has been switched on.
>Why ? 
>Here the input card (my best guess):

>OPT-PROP      -100.0    -100.0    -100.0      34.0      34.0       
>0.0SENSITIV 

SENSITIV is setting the quantum efficiency as a function of photon energy
*OVERALL* thru 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). I appended below an explaination taken
directly from the code of these options, hoping it can be understood.
So the use is:
a) you are at position whatever, material whatever
b) Cerenkov (or scintillation) photons are going to be produced
c) the maximum quantum efficiency of the problem (set in the same card
   as WHAT(5), IT IS MEANINGFUL EVEN IF queffc is used, see below) 
   is used immediately to reduce over the whole energy spectrum the 
   photon generation probability
d) the detailed efficiency set again by the OPT-PROD card with SENSITIV
   or by queffc is used for a further reduction when the actual energy
   of each individual photon is selected (rejecting it against
   quantum_efficiency_omega/max_quantum_efficiency)

Summarizing, the yes/no detection check is done AT PRODUCTION not at
detection, this in order to substantially cut down CPU times. If you want
all photons produced set the sensitivity to 1. If then you want to apply
a quantum efficiency curve AT DETECTION, you can do it with whichever
of the various user scoring options (ie fluscw/comscw) or off-line
if you wish.


>(2) queffc is called without a #region nor #material in the argument list. 
>I can get the current material from the geometry manager in it makes sense.
>Is this the right way to use it ?

I think this is already answered in the previous point. It is wrong, the
quantum efficiency curve provided by OPT-PROD with SENSITIV is applied 
AT PRODUCTION, not at detection, so it cannot know which material the
photon will eventually end up in. Further WHAT(5) must be anyway set
to the maximum quantum efficiency over the photon energy range under
consideration. If you prefer, you can use think queffc as a way
to provide a starting screening on the produced photons, ie using a "safe"
initial guess for the quantum efficiency (ie 20%) and then, at detection,
refining it thru more sophisticated curves, ie rejecting against
actual_quantum_efficiency/0.2 (this again can be done in any of the scoring
user routines). This makes sense of course if you have different quantum
efficiency curves for different detectors (use in queffc one that 
maximizes all them and then refine by rejection case by case), or if your
quantum efficiency is position/angle dependent upon arrival on the PMT
(again use inside queffc the one for the most efficient position/angle and 
refine by rejection at detection time, see below)

>(3) In Geant3 the detection of a photon is signaled by an energy 
>deposition in the first step inside the sensitive volume. Inefficiency is 
>simulated by having the energy deposition 0. and the photon is stopped.

>How is the quantum efficiency used in FLUKA ? 

I believe the previous point explains how the quantum efficiency is used in
FLUKA. Actually how to "detect" a photon detection is very different from
Geant3. FLUKA doesn't have the concept of "active" volumes, photons are
absorbed in those material where the user selected properties dictate
absorption, ie metals or materials with a non zero absorption cross section.
You can detect this absorption events in many different ways,
a) 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 for JTRACK=-1 and energy deposition (or ETRACK) in a given region
   (ie the photo cathode of the PMT). One can also apply whichever (extra)
   quantum efficiency selection, ie using COMSCW
b) through boundary crossing of particles -1 into 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, ie using FLUSCW
  
