*$ CREATE MGDRAW.FOR
*COPY MGDRAW
*                                                                      *
*=== mgdraw ===========================================================*
*                                                                     *
*-my-
* MGDRAW : to dump the whole history of particle showering, trajectories, 
*         interactions, etc ->  option USERDUMP is requested in input
*         boundary crossings, physical interactions, ...
*         A typical simple task is to draw particles trajectories.
*         to perform a calculation in two steps, where the second step 
*            uses the collision file as a source.
*  USDRAW is called after each particle interaction 
*        (if requested by the user with option USERDUMP, WHAT(4) >= 1.0)
*
*  MGDRAW, activated by option USERDUMP with WHAT(1) >= 100.,
*   usually writes a "collision tape", i.e. a file where all or selected
*   transport events are recorded. The default version (unmodified by the
*   user) offers several possibilities, selected by WHAT(3)
*   in USERDUMP.
*    Additional flexibility is offered by a user entry USDRAW, selected
*    by WHAT(4) in USERDUMP, interfaced with the most important physical
*    events happening during particle transport.
*
*    BXDRAW   - called at boundary crossings, 
*    EEDRAW   - called at event end, 
*    MGDRAW   - for trajectory drawing, 
*    ENDRAW   - for recording of energy depositions and 
*    SODRAW   - for recording of source events
*
*    the five  entries can be activated at the same time by setting 
*    USERDUMP with WHAT(3) = 0.0 and WHAT(4) >= 1.0
*     constitute a complete interface to the whole FLUKA transport.
*

*  the secondary particles produced is available in COMMON GENSTK,
*  that concerning delta rays produced by heavy ions 
*  (in which case the properties of the single electron produced are 
*  available in COMMON EMFSTK, with index NP).
*
*  heavy evaporation fragments (deuterons, 3H, 3He, , with 
*  JTRACK ID equal respectively to -3, -4, -5, -6) 
*  and ssion/fragmentation products generated in an inelastic 
*  interaction (with JTRACK = -7 to -12), which are all stored 
*  in COMMON FHEAVY, with index NPHEAV. 
*
*  To get the kinetic energy of particles with JTRACK < -6), 
*  one must subtract from their total energy (ETRACK in COMMON TRACKR) 
*  their fully stripped nuclear mass (AMNHEA in COMMON FHEAVY).
*
*  JPART -> see p.43 of manual
*
*      PROTON   1   p      Proton 2212
*      APROTON  2  anti-p  Antiproton -2212
*      ELECTRON 3  e-      Electron 11
*      POSITRON 4  e+      Positron -11
*      PHOTON   7  gamma   Photon 22
*      NEUTRON  8  n       Neutron 2112
*      ANEUTRON 9
*
*
*         ---------
*   /comput/ 
*         komput       = system (1=OpenVms,  2=ibm-vm,   3=ibm-mvs,    *
*                                4=cray,     5=unix-aix, 6=unix-hp,    *
*                                7=unix-sun, 8=DEC-unix, 9=Linux...)   *
*   COMMON / CHMGDD / CFDRAW
*                     ...name from card ?
*   COMMON / MGDDCM / IODRAW, LDDRAW, LEDRAW, LSDRAW, LUDRAW
*     iodraw     = logical unit for the output of the trajectory/energy*
*                  deposition dump                                     *
*       
*        *   (optical photons have always id = -1 ... check for JTRACK = -1
*      
*    SBGL - ?
*                     -------------
*
*    default version of MGDRAW :
*         Case 1 (First variable > 0): continuous energy deposition
*         Case 2 (First variable = 0): point energy deposition
*         Case 3 (First variable < 0): source particles
*
*           In Case 1, the following variables are written
*               ------
* First record:
*      NTRACK, MTRACK, JTRACK, ETRACK, WTRACK 
*                            (three integers and two real variables)
* Next record:
*      (XTRACK(I), YTRACK(I), ZTRACK(I), I = 0, NTRACK),
*      (DTRACK(J), J = 1, MTRACK), CTRACK
*             (that is, NTRACK + MTRACK + 1 real variables)
*
*   NTRACK : number of track segments
*   MTRACK : number of energy deposition events along the track
*   JTRACK : type of particle
*   ETRACK : total energy of the particle
*   WTRACK : weight of the particle
*    XTRACK(I), YTRACK(I), ZTRACK(I) = coordinates defining the upstream 
*                                      end of the (I+1)th segment;
*                    for I = NTRACK, the end of the trajectory
*    DTRACK(J) = energy deposition in the Jth deposition event along 
*                the trajectory
*   CTRACK : total length of the curved path
*         ---------
*
*               In Case 2, the following variables are written:
*                 ---------
* First record:
*       0, ICODE, JTRACK, ETRACK, WTRACK 
*                (three integers and  two real variables)
* Next record:
*      XSCO, YSCO, ZSCO, RULL (4 real variables)
*
*   JTRACK, ETRACK, WTRACK have the meaning explained above
*   XSCO, YSCO, ZSCO = coordinates of the energy deposition point
*   RULL = amount of energy deposited
*   ICODE indicates the type of point event giving rise to energy 
*         deposition, as explained below:  .... see Icode definitions ...
*
*
*             In Case 3, the following variables are written:
*                ------
* First record:
*       -NCASE, NPFLKA, NSTMAX, TKESUM, WEIPRI, 
*                      (three integers and  two real variables)
* Next record:
*      (ILOFLK(I), ETOT(I), WTFLK(I), XFLK(I), YFLK(I), ZFLK(I), TXFLK(I), 
*            TYFLK(I),TZFLK(I), I = 1, NPFLKA) 
*            (NPFLKA times (that is, one integer +  8 real variables))
*
* NCASE     = number of primaries treated so far (including the current one)
* NPFLKA    = number of particles in the stack
* NSTMAX    = maximum number of particles in stack so far
* TKESUM    = total kinetic energy of the primaries of a user written SOURCE
* WEIPRI    = total weight of the primaries handled so far
* ILOFLK(I) = type of the Ith stack particle (see 5.1)
* ETOT(I)   = total energy of Ith stack particle
* XFLK(I), YFLK(I), ZFLK(I) = source coordinates for the Ith stack particle
* TXFLK(I), TYFLK(I), TZFLK(I) = direction cosines of the Ith stack particle
*
*
*     /Genstk/ is the storage for secondaries created in hadronic      *
*              events                                                  *
*        Np        = total number of secondaries                       *
*       Kpart (ip) = (Paprop) id of the ip_th secondary                *
*      Infext (ip) = possible extra infos for the ip_th secondary      *
*         Cxr (ip) = x-axis direction cosine of the ip_th secondary    *
*         Cyr (ip) = y-axis direction cosine of the ip_th secondary    *
*         Czr (ip) = z-axis direction cosine of the ip_th secondary    *
*      Cxrpol (ip) = x-axis direction cosine of the ip_th secondary    *
*                    polarization vector (rest frame when applicable)  *
*      Cyrpol (ip) = y-axis direction cosine of the ip_th secondary    *
*                    polarization vector (rest frame when applicable)  *
*      Czrpol (ip) = z-axis direction cosine of the ip_th secondary    *
*                    polarization vector (rest frame when applicable)  *
*         Tki (ip) = laboratory kinetic energy of ip_th secondary (GeV)*
*         Plr (ip) = laboratory momentum of the ip_th secondary (GeV/c)*
*         Wei (ip) = statistical weight of the ip_th secondary         *
*      Agesec (ip) = "age" of the ip_th secondary with respect to the  *
*                    interaction time                                  *
*        Tv        = excitation energy (GeV)                           *
*        Tvcms     = actual excitation energy of the residual nucleus  *
*        Tvrecl    = recoil kinetic energy of the residual nucleus     *
*        Tvheav    = recoil kinetic energies of heavy (2-H, 3-H, 3-He, *
*                    4-He) fragments after evaporation                 *
*        Tvbind    = approximate energy wasted in nuclear binding      *
*                    effects (not yet operational)                     *
*
*     Common emfstk (EMF STacK) for EMF                                *
*      COMMON /EMFSTK/ ETEMF  (MESTCK), PMEMF  (MESTCK), XEMF   (MESTCK),
*     &                YEMF   (MESTCK), ZEMF   (MESTCK), UEMF   (MESTCK),
*     &                VEMF   (MESTCK), WEMF   (MESTCK), DNEAR  (MESTCK),
*     &                UPOL   (MESTCK), VPOL   (MESTCK), WPOL   (MESTCK),
*     &                USNRML (MESTCK), VSNRML (MESTCK), WSNRML (MESTCK),
*     &                WTEMF  (MESTCK), AGEMF  (MESTCK), CMPEMF (MESTCK),
*     &                RDLYEM (MESTCK), DCHEMF (MESTCK),
*     &                ESPARK (MKBMX1,MESTCK),
*     &                IESPAK (MKBMX2,MESTCK),           ICHEMF (MESTCK),
*     &                IREMF  (MESTCK), IRLATT (MESTCK), NHPEMF (MESTCK),
*     &                LLOEMF (MESTCK), LOUEMF (MESTCK), LRDEMF (MESTCK),
*     &                LCHEMF (MESTCK), NPEMF , NPSTRT
*
*     TRACKs Recording  
*          J0trck = see above                                          *
*          Ifltrk = flag used for internal debugging (trying to solve  *
*                   possible residual issues with Mgdraw driven        *
*                   quenching)                                         *
*          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)                                          *
*          Wcinou = charge  algebraic balance of interactions (for     *
*                   all interactions)                                  *
*          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) *
*          Tkniel = fraction of energy deposition going into (unres-   *
*                   tricted) NIEL                                      *
*          Tkrnie = fraction of energy deposition going into (res-     *
*                   tricted) NIEL                                      *
*          Tkedpa = fraction of energy deposition going into DPAs      *
*          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)   *
*          Iprodc = flag for prompt(=1)/radioactive products(=2)       *
*          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      *
*                                                                      *
*                  -----------------
*
*   > So XTRACK(0), YTRACK(0) and ZTRACK(0) are coordinates of the birth
*   > location of photons generated and I don't need to use COMMON GENSTK in
*   > USDRAW entry?
*
*  that's right. By the way, in (GENSTK) you cannot find the birth location,
*  that is rather XSCO, YSCO, ZSCO in USDRAW. The latter, when needed, allows
*  to identify the generation process (ICODE).
*  A "signal" means energy deposited in one or more trigger regions 
*  by the same primary particle and descendants (i.e., between the same calls).
*  The output of DETECT is a distribution of energy deposited per event 
*  in the region(s) making up the detector in (anti-)coincidence with 
*  a signal larger than a given cuto or threshold in the trigger region(s).
*
*
*
*                 -----------------
*                 -----------------
*                 -----------------
*
*      USERDUMP:  quenching ->  For SDUM = UDQUENCH - Birks parameter
*
*                 -------------------
*       http://www.fluka.org/fluka.php?id=man_onl&sub=99
*
*      CALL USRDCI(IJ,IONA,IONZ,IONM)
*
*       The three integer values returned are the following ion properties:
*
*            IONA = mass number of the ion
*            IONZ = atomic number
*            IONM = flag for isomeric state
*      Based on their values, the user can decide or not to accept the
*             scoring. 
*
*                     ----------------
*
*       BXDRAW - is called at each boundary crossing
*                (if requested by the user  with USERDUMP, WHAT(3) < 7.0)
*                .... my output here ...     
*
*       ICODE : physical compartment originating the call, 
*               as in the MGDRAW entry
*        MREG : region from which the particle is exiting
*      NEWREG : region the particle is entering
*      XSCO, YSCO, ZSCO : coordinates of crossing point
*
*                         -------------
*  If name-based input is being used, the names corresponding 
*    to MREG and NEWREG can be obtained via a call to routine GEOR2N:
*
*            CALL GEOR2N (NUMREG, NAMREG, IERR)
*
*  where NUMREG (input variable) is a region number, and NAMREG 
*  (returned variable) is the corresponding region name 
*  (to be declared as CHARACTER*8). IERR is a returned error code: 
*   if = 0, the conversion is successful. Example:
*
*            .......................................
*            CHARACTER*8 MRGNAM, NRGNAM
*            .......................................
*            ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
*            CALL GEOR2N ( MREG,   MRGNAM, IERR1 )
*            CALL GEOR2N ( NEWREG, NRGNAM, IERR2 )
*            IF(IERR1 .NE. 0 .OR. IERR2 .NE. 0) STOP "Error in name conversion"
*            .......................................
*            IF(MRGNAM .EQ. "MyUpsREG" .AND. NRGNAM .EQ. "MyDwnREG") THEN
*            .......................................
*
*             -------------
*
*      As a summary:
*  1. The energy deposition along the track is scored in the main part of
*    mgdraw with:
*       CALL GEOR2N ( MREG, NRGNAM, IERR2 )
*       IF (NRGNAM(1:6).EQ.'SCINT0') THEN
*        DO J = 1,MTRACK
*         ERULL0 = ERULL0 + DTRACK(J)
*        ENDDO
*       ENDIF
*
*  2. The point-like energy deposition is scored in ENDRAW with:
*       CALL GEOR2N ( MREG, NRGNAM, IERR2 )
*       IF (NRGNAM(1:6).EQ.'SCINT0') THEN
*        ERULL0 = ERULL0 + RULL
*       ENDIF
*
* I want to add just a thing, for the future: take care if you want to use
* this method with other kind of primaries. If you have highly ionizing
* primaries -like protons- you must consider quenching effects in the
* scintillator. In this case, in order to correctly compute the light yield,
* the value you have to consider for the energy deposited in the scintillator-
* both in MGDRAW and in the ENDRAW entry - is the "quenched" one, obtained
* according with the Birks Law.
* You can consider this effect:
* (a) by switching on the quenching with the card USERDUMP with sdum=UDQUENCH
* (see the manual for details about the Birks coefficients, and look at Note 3
* pag. 239)
* (b) by using, instead of DTRACK(j) and RULL , the quenched values:
*   - DTQUEN(j,jbk) in the main MGDRAW   
*    jbk identifies the set of Birks parameters you input via USERDUMP. 
* jbk=1  if you use only one set, with two parameters given via what(1) and
* what(2)  (it is generalyy the case)
*  - RULLL  in the ENDRAW entry (please check this)
*   
*
*
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
      SUBROUTINE MGDRAW ( ICODE, MREG )
*                MGDRAW (trajectory dumping for drawing)

      INCLUDE '(DBLPRC)' ! * mathematical and   physical constants 
      INCLUDE '(DIMPAR)' ! * DIMension PARameters:  
      INCLUDE '(IOUNIT)' ! * files parameters
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2013      by        Alfredo Ferrari           *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     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   12-Nov-13      by        Alfredo Ferrari           *
*                                              INFN - Milan            *
*                                                                      *
*----------------------------------------------------------------------*
* http://www.fluka.org/content/manuals/online/13.2.13.html
*
      INCLUDE '(CASLIM)'    ! * to decide when to stop the run        
      INCLUDE '(COMPUT)'    ! * information about the computer used 
      INCLUDE '(SOURCM)'    ! * include some function
      INCLUDE '(FHEAVY)'    ! * storage for heavy secondaries 
      INCLUDE '(FLKSTK)'    ! * stack for the primaries
      INCLUDE '(GENSTK)'    ! * hadron GENerator STacK storage for secondaries 
      INCLUDE '(MGDDCM)'    ! * MGDraw Dump CoMmon
      INCLUDE '(PAPROP)'    ! * particle parameters
      INCLUDE '(QUEMGD)'    ! * BRKMG1 (MAXQMG), BRKMG2 (MAXQMG), NQEMGD, LQEMGD
      INCLUDE '(SUMCOU)'    ! * SUMmary COUnters common - inter. parameters
      INCLUDE '(TRACKR)'    ! * TRACKs Recording - track parameters
*
      INCLUDE '(EMFSTK)'    ! * single electron produced with index NP
*
      DIMENSION DTQUEN ( MXTRCK, MAXQMG )
*
      CHARACTER*20 FILNAM
      LOGICAL LFCOPE
      SAVE LFCOPE
      DATA LFCOPE / .FALSE. /
*
      INTEGER IONA, IONZ, IONM            ! my
      DOUBLE PRECISION EDEP, TRKL

      LOGICAL myFIRST
      DATA myFIRST / .TRUE. /
      SAVE myFIRST

**      SAVE T0

*
*----------------------------------------------------------------------*
*      Case 1 (First variable > 0): continuous energy deposition       *
*
*     Icode = 1: call from Kaskad   (hadrons and muons)                *
*     Icode = 2: call from Emfsco   (e-, e+ and photons)               *
*     Icode = 3: call from Kasneu   (low-energy neutrons)              *
*     Icode = 4: call from Kashea   (heavy ions)                       *
*     Icode = 5: call from Kasoph   (optical photons)                  *
*
*  MREG : current region                                               *
*
*----------------------------------------------------------------------*
*                                                                      *
**      IF ( .NOT. LFCOPE ) THEN
**         LFCOPE = .TRUE.
**         IF ( KOMPUT .EQ. 2 ) THEN
**            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
**         ELSE
**            FILNAM = CFDRAW
**         END IF
C IODRAW is the logical unit for the output of the trajectory/energy
c deposition dump, maybe set by USERDUMP
**         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
**     &          'UNFORMATTED' )
**      END IF
**      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)
***  +-------------------------------------------------------------------*
***  |  Quenching is activated         (Birks law).
**      IF ( LQEMGD ) THEN
**         IF ( MTRACK .GT. 0 ) THEN
**            RULLL  = ZERZER
**            CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
**            WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)), I = 1, MTRACK ),
**     &                         JBK = 1, NQEMGD )
**         END IF
**      END IF
***  |  End of quenching
***  +-------------------------------------------------------------------*
*
* -my-start
*
* --> check  SNGL (ETRACK) procedure only ...
*     --------------------------
      IF ( .NOT. LFCOPE ) THEN
         WRITE(LUNOUT,*) 'Version Qi.W of Routine MGDRAW_TRACK',
     &                   NTRACK, MTRACK, JTRACK, SNGL (ETRACK)
         LFCOPE = .TRUE.
         IF ( KOMPUT .EQ. 2 ) THEN
c KOMPUT is the computer system type.
            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
         ELSE
            FILNAM = CFDRAW
         END IF
C IODRAW is the logical unit for the output of the trajectory/energy
c deposition dump, maybe set by USERDUMP
         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
     &          'UNFORMATTED' )
 
         OPEN ( UNIT = 90, FILE = 'my-inform', STATUS = 'NEW', FORM =
     &          'FORMATTED' )
      END IF
* -------
      IF ( JTRACK .EQ. 8) THEN            ! --- check neutron
         WRITE (IODRAW, '(3(A,G10.3))') 'J: ', JTRACK, ' L: ', LTRACK
         WRITE (IODRAW, '(7(A,G10.3))') ( ' X: ', SNGL (XTRACK (I)),
     &                 ' Y: ', SNGL (YTRACK (I)),
     &                 ' Z: ', SNGL (ZTRACK (I)), I = 0, NTRACK ),
     &                 ' C: ', SNGL (CTRACK)
          WRITE (IODRAW, '(2(A, G10.3))') 'M: ', MTRACK, ' E: ', ETRACK 
      END IF
*----------


      WRITE (90,*) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
     &               SNGL (WTRACK)
      WRITE (90,*) NTRACK, MTRACK, JTRACK, ETRACK,
     &               WTRACK
      WRITE (90,*) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
     &                 SNGL (ZTRACK (I)), I = 0, NTRACK ),
     &               ( SNGL (DTRACK (I)), I = 1, MTRACK ),
     &                 SNGL (CTRACK)

      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)
*--------------------------
*
*       get neutrons/photons at W (tungsten foil) and print it
* regions:
*
*     3  Rbox   
*     4  Rtarg 
*     5  Rcore 
*     6  Rwall 
*     7  Rhole 
*     8  Rscin 
*     9  Rcount1 
*    10  RW            - W foil after concrete
*                --------
*
      IF (MREG .EQ. 10) THEN       ! select W foil... any code...
      IF (JTRACK.EQ.7 .OR. JTRACK.EQ.8) THEN   !  get photon and neutron
              
        IF ( myFIRST ) THEN
           myFIRST = .FALSE.
*           open (unit = 91, file = 'dat-all.dat',STATUS='UNKNOWN', FORM=
*     &    'FORMATTED')
           WRITE(91,*) " ===> 1-st entry: USDRAW - general"
        ENDIF
     

        IONA = 0 
        IONZ = 0
        IONM = 0

        CALL USRDCI( JTRACK, IONA, IONZ, IONM)   ! - check ions ...
        WRITE(91, * ) '-my- general'  
        WRITE(91, '(A5,2x,i3,2x,A5,2x,i3,2x,A5,i3)')
     &         'ion=',IONA,'Z = ',IONZ,'flag=',IONM          

        TRKL= 0.0

        do 50 I=1, NTRACK                       ! total L for track
          TRKL=TRKL+
     &    SQRT( (XTRACK(I)-XTRACK(I-1))**2 +
     &    (YTRACK(I)-YTRACK(I-1))**2 + (ZTRACK(I)-ZTRACK(I-1))**2)
50      continue                            ! total edep
        do 40 I=1, MTRACK
           EDEP=EDEP+DTRACK(I)
40      continue

        WRITE(91, '(A10,2x,i3,2x,D10.3,2x,D10.3,2x,D10.3)')
     &            'track,M,E,L=',JTRACK,AM(JTRACK),ETRACK,TRKL           

      ENDIF                      ! my particles
      END IF                     ! end of W foil

*
* -my- end
*
      RETURN
*
*======================================================================*
*                                                                      *
*     Boundary-(X)crossing DRAWing:                                    *
*                                                                      *
*     Icode = 1x: call from Kaskad       (hadrons and muons)           *
*             19: boundary crossing                                    *
*     Icode = 2x: call from Emfsco       (e-, e+ and photons)          *
*             29: boundary crossing                                    *
*     Icode = 3x: call from Kasneu       (low-energy neutrons)         *
*             39: boundary crossing                                    *
*     Icode = 4x: call from Kashea                                     *
*             49: boundary crossing                                    *
*     Icode = 5x: call from Kasoph                                     *
*             59: boundary crossing                                    *
*                                                                      *
*======================================================================*
* http://www.fluka.org/content/manuals/online/13.2.13.html
*                                                                      *
      ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
*           BXDRAW (boundary crossing dumping)
*
* -my-start
*
*      CHARACTER*8 reg1, reg2
*
*      CALL GEOR2N ( MREG, reg1, IERR1 )
*      CALL GEOR2N ( NEWREG, reg2, IERR2 )
*      IF(IERR1 .NE. 0 .OR. IERR2 .NE. 0) STOP "Error in name conversion"
*
*      IF((JTRACK .EQ. -6) .OR. (JTRACK .EQ. -2)) THEN
*        IF (reg1 .EQ. "Rbox" .AND. reg2 .EQ. "RW") THEN     
*           WRITE(81,7) JTRACK, XSCO, YSCO, ZSCO, CXTRCK, CYTRCK, 
*     &          CZTRCK,ETRACK-AM(JTRACK),WTRACK
*7               FORMAT(I8, 8E25.5)
*        ENDIF
*      ENDIF
* --------- 

      CALL USRDCI( JTRACK, IONA, IONZ, IONM)

      IF (NEWREG.EQ.10) THEN                   ! select W foil only
      IF (JTRACK.EQ.7 .OR. JTRACK.EQ.8) THEN   !  get photon and neutron

         IF ( myFIRST ) THEN
           myFIRST = .FALSE.
*           open (unit = 92, file = 'dat-b.dat',STATUS='UNKNOWN', FORM=
*     &    'FORMATTED')
            WRITE(92,*) " ===> 1-st entry: USDRAW - boundary crossing"
         ENDIF

         IF(ETRACK.GT.AM(JTRACK)) THEN 
           WRITE(92, * ) '-my-'  
*           WRITE(92,'(
*     &       A5,2X,I4,1X,I2,13X,I2,13X,D10.3,D10.3,/,
*     &       A5,1X,D10.3,1X,D10.3,1X,D10.3,/,
*     &       A5,1X,D10.3,1X,D10.3,1X,D10.3)')
           WRITE(92,'(A5,2X,3(I3,1x),3(F10.3,1x),/,
     &       A5,1X,3(F10.3,1x),/,
     &       A5,1X,3(F10.3,1x))')
     &       '-PRT-',JTRACK,IONA,IONZ,AM(JTRACK),ETRACK,
     &       ETRACK - AM(JTRACK),
     &       '-POS-',XSCO,YSCO,ZSCO,
     &       '-COS-',CXTRCK,CYRCK,CZTRCK
         ENDIF

      ENDIF
      ENDIF

*
* -my-end
*
      RETURN
*
*======================================================================*
*                                                                      *
*     Event End DRAWing:                                               *
*                                                                      *
*======================================================================*
* http://www.fluka.org/content/manuals/online/13.2.13.html
*                                                                      *
      ENTRY EEDRAW ( ICODE )
*           EEDRAW (event end)
*
* -my-start
*
*           open (unit = 93, file = 'dat-e.dat',STATUS='UNKNOWN', FORM=
*     &    'FORMATTED')

      IF((ETRACK.GT.AM(JTRACK))) THEN      ! check end of event...
*        IF((EDEP.GT.1.0E-9)) THEN
           WRITE(93,'(A5,I3,A3,1X,D10.0,/,A3,1X,D10.3)') 
     &         'END = ',JTRACK,  'EDP', EDEP, 
     &                           'TRK', TRKL
*        END IF
      END IF

      EDEP=0               ! update values
      TRKL=0
  
*
* -my-end
*
      RETURN
*
*======================================================================*
*           Case 2 (First variable = 0): point energy deposition       *
*
*     ENergy deposition DRAWing:                                       *
*                                                                      *
*     Icode = 1x: call from Kaskad            (hadrons and muons)      *
*             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            (e-, e+ and photons)     *
*             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            (low-energy neutrons)    *
*             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                                            *
*                                                                      *
*======================================================================*
* http://www.fluka.org/content/manuals/online/13.2.13.html
*                                                                      *
      ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO )
*           ENDRAW (dumping of energy deposition at a point)
*
*      IF ( .NOT. LFCOPE ) THEN
*         LFCOPE = .TRUE.
*         IF ( KOMPUT .EQ. 2 ) THEN
*            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
*         ELSE
*            FILNAM = CFDRAW
*         END IF
*         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
*     &          'UNFORMATTED' )
*      END IF
*      WRITE (IODRAW)  0, ICODE, JTRACK, SNGL (ETRACK), SNGL (WTRACK)
*      WRITE (IODRAW)  SNGL (XSCO), SNGL (YSCO), SNGL (ZSCO), SNGL (RULL)
**  +-------------------------------------------------------------------*
**  |  Quenching is activated : calculate quenching factor
**  |  and store quenched energy in DTQUEN(1, jbk)
*      IF ( LQEMGD ) THEN
*         RULLL = RULL
*         CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
*         WRITE (IODRAW) ( SNGL (DTQUEN(1, JBK)), JBK = 1, NQEMGD )
*      END IF
**  |  end quenching
**  +-------------------------------------------------------------------*

*
* -my-start
*
      IF (MREG.EQ.10) THEN             ! W-foil
      if (JTRACK.EQ.7 .OR. JTRACK.EQ.8) THEN   !  get photon and neutron

          EDEP=EDEP+RULL

      ENDIF
      END IF

*
* -my-end
*
      RETURN
*
*======================================================================*
*                                                                      *
*     SOurce particle DRAWing:                                         *
*               Case 3 (First variable < 0): source particles          *
*======================================================================*
* http://www.fluka.org/content/manuals/online/13.2.13.html
*
      ENTRY SODRAW
*           SODRAW (source particle dumping)
*
*      IF ( .NOT. LFCOPE ) THEN
*         LFCOPE = .TRUE.
*         IF ( KOMPUT .EQ. 2 ) THEN
*            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
*         ELSE
*            FILNAM = CFDRAW
*         END IF
*         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
*     &          'UNFORMATTED' )
*      END IF
*      WRITE (IODRAW) -NCASE, NPFLKA, NSTMAX, SNGL (TKESUM),
*     &                SNGL (WEIPRI)
**  +-------------------------------------------------------------------*
**  |  (Radioactive) isotope: it works only for 1 source particle on
**  |  the stack for the time being
*      IF ( ILOFLK (NPFLKA) .GE. 100000 .AND. LRADDC (NPFLKA) ) THEN
*         IARES  = MOD ( ILOFLK (NPFLKA), 100000  )  / 100
*         IZRES  = MOD ( ILOFLK (NPFLKA), 10000000 ) / 100000
*         IISRES = ILOFLK (NPFLKA) / 10000000
*         IONID  = ILOFLK (NPFLKA)
*         WRITE (IODRAW) ( IONID,SNGL(-TKEFLK(I)),
*     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
*     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
*     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
*     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
**  |
**  +-------------------------------------------------------------------*
**  |  Patch for heavy ions: it works only for 1 source particle on
**  |  the stack for the time being
*      ELSE IF ( ABS (ILOFLK (NPFLKA)) .GE. 10000 ) THEN
*         IONID = ILOFLK (NPFLKA)
*         CALL DCDION ( IONID )
*         WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-IONID)),
*     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
*     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
*     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
*     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
**  |
**  +-------------------------------------------------------------------*
**  |  Patch for heavy ions: ???
*      ELSE IF ( ILOFLK (NPFLKA) .LT. -6 ) THEN
*         WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-ILOFLK(NPFLKA))),
*     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
*     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
*     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
*     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
**  |
**  +-------------------------------------------------------------------*
**  |
*      ELSE
*         WRITE (IODRAW) ( ILOFLK(I), SNGL (TKEFLK(I)+AM(ILOFLK(I))),
*     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
*     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
*     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
*     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
*      END IF
**  |
**  +-------------------------------------------------------------------*

      RETURN
*
*======================================================================*
*                                                                      *
*     USer dependent DRAWing:                                          *
*                                                                      *
*     Icode = 10x: call from Kaskad           (hadrons and muons)      *
*             100: elastic   interaction secondaries                   *
*             101: inelastic interaction secondaries                   *
*             102: particle decay  secondaries                         *
*             103: delta ray  generation secondaries                   *
*             104: pair production secondaries                         *
*             105: bremsstrahlung  secondaries                         *
*             110: radioactive decay products                          *
*     Icode = 20x: call from Emfsco            (e-, e+ and photons)    *
*             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                  *
*             237: mu pair     production secondaries                  *
*     Icode = 30x: call from Kasneu           (low-energy neutrons)    *
*             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          *
*                                                                      *
*======================================================================*
* http://www.fluka.org/content/manuals/online/13.2.13.html
*
      ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )


      IF (MREG.EQ.10) THEN             ! select W foil only
      IF (JTRACK.EQ.7 .OR. JTRACK.EQ.8) THEN   !  get photon and neutron

        IF ( myFIRST ) THEN
           myFIRST = .FALSE.
*           open (unit = 94, file = 'dat-us.dat',STATUS='UNKNOWN', FORM=
*     &    'FORMATTED')
           WRITE(94,*) " ===> 1-st entry: USDRAW - user"
        ENDIF

**        IF ( .NOT. LFCOPE ) THEN
**           LFCOPE = .TRUE.
**           IF ( KOMPUT .EQ. 2 ) THEN
**              FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
**           ELSE
**              FILNAM = CFDRAW
**           END IF
**           OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
**     &          'UNFORMATTED' )
**        END IF

*        WRITE (IODRAW,*) 'Prim-- ',NCASE,JTRACK,NPFLKA,ICODE,
*        WRITE (94,*) 'Prim-- ',NCASE,JTRACK,NPFLKA,ICODE,
*     &        NP,NPHEAV,ETRACK,AM(JTRACK),XSCO,YSCO,ZSCO,WTRACK
        WRITE (94,"(A7,2x,6(I3,1x),6(f10.3,1x))") 
     &        'Prim-- ',NCASE,JTRACK,NPFLKA,ICODE,
     &        NP,NPHEAV,ETRACK,AM(JTRACK),XSCO,YSCO,ZSCO,WTRACK

*
* secondaries ...
*
        do 10 ip = 1, NP
*          WRITE (IODRAW,*) 'Seond-- ',NCASE,JTRACK,
*          WRITE (94,*) 'Seond-- ',NCASE,JTRACK,
*     &          ip,ICODE,ZPART,APART,TKI(ip),CXR(ip),CYR(ip),
*     &          CZR(ip),WEI(ip)
           WRITE (94,"(A9,2x,5(I3,1x),6(F10.3),1x))") 
     &           'Secon-- ',NCASE,JTRACK,ip,ICODE,ZPART,
     &          APART,TKI(ip),CXR(ip),CYR(ip),CZR(ip),WEI(ip)
10      continue
*
* electrons-photons
*
        do 20 ip=1,NP
*          WRITE (IODRAW,*)'Elec-phot-- ', ICHEMF(ip),
*          WRITE (94,*)'Elec-phot-- ', ICHEMF(ip),
*     &         ICODE, ETEMF(ip), XEMF(ip), YEMF(ip),ZEMF(ip)
          WRITE (94,"(A12,2x,2(I3,1x),4(F10.3,1x))")
     &          'Elec-phot-- ', ICHEMF(ip),ICODE,
     &           ETEMF(ip), XEMF(ip), YEMF(ip),ZEMF(ip)
20      continue
*
* ions
*
*        do 30 ip=1,NPHEAV
*           WRITE (IODRAW,*) 'Heavy--', NCASE, JTRACK,
*           WRITE (94,*) 'Heavy--', NCASE, JTRACK,
*     &         ip,ICODE,KHEAVY(ip),ICHEAV(KHEAVY(ip)),IBHEAV(KHEAVY(ip)),
*     &         TKHEAV(ip)
*30      continue
*
      ENDIF                         ! my particles
      ENDIF                         ! end of foil

****************************************************
**    'MGDRAW: dump neutrons for photoelectric int.'
**
**   select photonuclear interaction
**   IF ( ICODE.EQ.101 .AND. JTRACK.EQ.7 ) THEN
**       N = 0
**        DO 10 IP = 1, NP
**    *     check that at least a neutron is among the daughter
**          IF ( KPART(IP) .EQ. 8 ) N = N + 1
** 10      CONTINUE
**
**        IF (N .GT. 0) THEN
**           WRITE (94,*)  "NCASE, ATRACK, N_neutrons", NCASE, ATRACK, N
**  *     store time of the photonuclear interaction
**         T0 = ATRACK
**       END IF
**     END IF
******************************

* No output by default:
      RETURN


*=== End of subrutine Mgdraw ==========================================*
      END


