*$ CREATE MGDRAW.FOR
*COPY MGDRAW
*                                                                      *
*=== mgdraw ===========================================================*
*                                                                     *
* 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.
*
*
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
      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 '(CTITLE)' 
      INCLUDE '(EMFSTK)'    ! * single electron produced with index NP
*
***      DIMENSION DTQUEN ( MXTRCK, MAXQMG )
      DIMENSION DTQUEN ( MXTRCK, MAXQMG ), AMHELP ( MFSTCK )
*
*
      CHARACTER*20 FILNAM
      LOGICAL LFCOPE
      SAVE LFCOPE
      DATA LFCOPE / .FALSE. /
*
* ====>///////// my part ///////////////
*
      INTEGER IONA, IONZ, IONM            ! my
      DOUBLE PRECISION EDEP, TRKL
      DOUBLE PRECISION  WX, WY, WZ, WMIN, RMIN

      REAL Ekin, Tmks                     ! in eV  and microsecond
      REAL DEL, STEP
      REAL PMAS

      LOGICAL myFIRST, myFIRST2
      DATA myFIRST / .TRUE. /
      DATA myFIRST2 / .TRUE. /
      SAVE myFIRST, myFIRST2

      INTEGER ICH

      CHARACTER*8 reg1, reg2

*       ----------------------------
* save in /HIST/ 5 variables for energy-resolution :
*                E1, E2 and L, E, T for L= 1-100 cm
*              - for L : entry input
*                others 2 : E(L)  and T(L) -> then we know <E> and <T>
*      each histos : H1, H2, ... H14  energy intervals :
*          H1 : 1-4 ev
*          H2 : 4-15 ev
*          H3 : 15-50 ev
*          H4 : 50-120 ev
*          H5 : 120-300 ev
*          H6 : 300-10^3 ev
*          H7 : 10^3-3*10^3 ev
*          H8 : 3*10^3-7*10^3 ev
*          H9 : 7*10^3-15*10^3 ev
*         H10 : 15*10^3-3*10^4 ev
*         H11 : 3*10^4-6*10^4 ev
*         H12 : 6*10^4-10^5 ev
*         H13 : 10^5-2*10^5 ev
*         H14 : 2*10^5-4*10*5 ev
*
**      COMMON /ECUT/ EC(14,2), HH(14,3,100)
cc      COMMON /ECUT/ EC(14,2), HH(14,3,200)     ! 200 bin with 1 cm step
      COMMON /ECUT/ EC(14,2), HH(14,400)     ! 400 bin with 0.5 cm step

*      COMMON /HIST/ H1(100,3), H2(100,3), H3(100,3), H4(100,3),
*     &               H5(100,3), H6(100,3), H7(100,3), H8(100,3),
*     &               H9(100,3), H10(100,3), H11(100,3), H12(100,3),
*     &               H13(100,3), H14(100,3)
*
      SAVE /ECUT/
*      SAVE /HIST/
*
* <=====//////////// my part - end ///////////
*
**      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
***  +-------------------------------------------------------------------*
* --> check  SNGL (ETRACK) procedure only ... return real from double !
*     --------------------------
*

      IF ( .NOT. LFCOPE ) THEN
         LFCOPE = .TRUE.
         IF ( KOMPUT .EQ. 2 ) THEN
            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
         ELSE
            FILNAM = CFDRAW
         END IF

**         WRITE(LUNOUT,*) 'Version Qi.W of Routine MGDRAW_TRACK',
**     &        NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
**     &        SNGL (CMTRCK), SNGL (ATRACK)
**         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' )
*
* ===>///////////// -my-start  ////////////////////
*
*   1 mks = 10^-6 s   (микросекунда)
*   1 мс  = 10^-3 s   (миллисекунда)
*
* 
         OPEN ( UNIT = 90, FILE = '90-my-infor', STATUS = 'NEW', FORM =
     &          'FORMATTED' )

         write (90,*) "...MGDRAW start.. : LUNOUT= , IODRAW= :", 
     &                 LUNOUT,IODRAW
         WRITE(90,'(A32)') RUNTIM
         write (90,*) 

**         write (90,*) "...print info for N at W-foil (after concrete) "
**         write (90,*) "  L-1 : NTRACK, MTRACK, JTRACK, ETRACK, WTRACK,",
**     &                "  CMTRCK(length),  ATRACK(time) "
**         write (90,*) "  L-2 : XTRACK, YTRACK, ZTRACK .. * NTRACK " 
**         write (90,*) 

** comment: LUNOUT = 11 (from output)
**          IODRAW = 49 - it is set at USERDUMP card ...
**
**
**  my-90infor - file is large : 5 000 ev -> 2 Gb ! size ...
**                ...flair is closed ! when open this file ...
**
*
* ...open files for production ....
*    ------------------------
*
         OPEN ( UNIT = 85, FILE = '85-C1', STATUS = 'NEW', 
     &                     FORM ='FORMATTED' )
         OPEN ( UNIT = 86, FILE = '86-C2', STATUS = 'NEW', 
     &                     FORM ='FORMATTED' )
         OPEN ( UNIT = 87, FILE = '87-C3', STATUS = 'NEW', 
     &                     FORM ='FORMATTED' )

         OPEN ( UNIT = 89, FILE = '89-C4', STATUS = 'NEW', 
     &                     FORM ='FORMATTED' )


         write (90,*) 
         write (90,*) "... open files : 85... -  89.... "
*         write (90,*) "... open files : 85...  and  86.., 87... "
         write (90,*) 


*
*  print head into files ...
*

*
         write (85,*) " C1 : all particles "
*         write (85,*) " Nev - Ekin - Ltot- Time - weight ",
*     &                  "- X - Y - Z  "
**         write (85,*) " Det  -  Nev - Ekin - Ltot- Time - weight ",
**     &                  "- X - Y - Z  "
         write (85,*) " Det - Type -  Nev - Ekin - Ltot- L.Orig -",
     &                " Time - w - X - Y - Z  - cosX - cosY - cosZ "
         write (85,*) " NEWREG, JTRACK, NCASE, Ekin, ",
     &                  " CTRACK, CMTRCK, ATRACK, ",
     &                  "- WTRACK, XSCO, YSCO, ZSCO ",
     &                  "  CXTRCK, CYTRCK, CZTRCK "
         write (85,*) "  "

*
         write (86,*) " C2 : all particles "
**         write (86,*) " Nev - Ekin - Ltot- Time - weight ",
**     &                  "- X - Y - Z  "
***         write (86,*) " Det  -  Nev - Ekin - Ltot- Time - weight ",
***     &                  "- X - Y - Z  "
         write (86,*) " Det - Type -  Nev - Ekin - Ltot- L.Orig -",
     &                " Time - w - X - Y - Z  - cosX - cosY - cosZ "
         write (86,*) " NEWREG, JTRACK, NCASE, Ekin, ",
     &                  " CTRACK, CMTRCK, ATRACK, ",
     &                  "- WTRACK, XSCO, YSCO, ZSCO ",
     &                  "  CXTRCK, CYTRCK, CZTRCK "
         write (86,*) "  "
*

*
         write (87,*) " C3 : "
         write (87,*) " Det - Type -  Nev - Ekin - Ltot- L.Orig -",
     &                " Time - w - X - Y - Z  - cosX - cosY - cosZ "
         write (87,*) " NEWREG, JTRACK, NCASE, Ekin, ",
     &                  " CTRACK, CMTRCK, ATRACK, ",
     &                  "- WTRACK, XSCO, YSCO, ZSCO ",
     &                  "  CXTRCK, CYTRCK, CZTRCK "
         write (87,*) "  "

*
         write (89,*) " C4 : "
         write (89,*) " Det - Type -  Nev - Ekin - Ltot- L.Orig -",
     &                " Time - w - X - Y - Z  - cosX - cosY - cosZ "
         write (89,*) " NEWREG, JTRACK, NCASE, Ekin, ",
     &                  " CTRACK, CMTRCK, ATRACK, ",
     &                  "- WTRACK, XSCO, YSCO, ZSCO ",
     &                  "  CXTRCK, CYTRCK, CZTRCK "
         write (89,*) "  "

*
* ...for histos :
*
         OPEN ( UNIT = 88, FILE = '88-hist', STATUS = 'NEW', 
     &                     FORM ='FORMATTED' )

         write (88,*) " each N= ??? ev : arrays-histos: "
         write (88,*) "      del L - distributions : N, E, T "
         write (88,*) "  "

*
*        .... and fill /ECUT/  for resolution (eV) :
*
         EC(1,1) = 1.0                         ! ev
         EC(1,2) = 4.0
           EC(2,1) = 4.0
           EC(2,2) = 15.0
             EC(3,1) = 15.0
             EC(3,2) = 50.0
               EC(4,1) = 50.0
               EC(4,2) = 120.0
         EC(5,1) = 120.0                         ! ev
         EC(5,2) = 300.0
           EC(6,1) = 300.0
           EC(6,2) = 1E+03
             EC(7,1) = 1E+03
             EC(7,2) = 3E+03
               EC(8,1) = 3E+03
               EC(8,2) = 7E+03
         EC(9,1) = 7E+03                         ! ev
         EC(9,2) = 15E+03
           EC(10,1) = 15E+03
           EC(10,2) = 3E+04
             EC(11,1) = 3E+04
             EC(11,2) = 6E+04
               EC(12,1) = 6E+04
               EC(12,2) = 1E+05
                 EC(13,1) = 1E+05
                 EC(13,2) = 2E+05
                   EC(14,1) = 2E+05
                   EC(14,2) = 4E+05

        DO I=1,14
cc          DO J=1,3
**            DO K=1,100
cc            DO K=1,200                   ! 200 bins
cc               HH(I,J,K) = 0.0
cc            ENDDO
            DO K=1,400                   ! 400 bins
               HH(I,K) = 0.0
            ENDDO
cc          ENDDO
        ENDDO

*
*   <===//////////////// -my-start  - end  /////////////////
*

      END IF

*  +-------------------------------------------------------------------*
*  |  Quenching is activated
      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
*  +-------------------------------------------------------------------*



* =============================
*  ... make dump-file output small...

 
***      IF ( NCASE .LT. 100) THEN       ! * print at the start only ...
***      IF ( NCASE .LT. 10) THEN       ! * print at the start only ...
***      IF ( NCASE .LT. 2) THEN       ! * print at the start only ...
      IF ( NCASE .LT. 2 .AND. NTRACK.LT.1 ) THEN   
                                            ! * print at the start only ...
                                            ! * it is for DUMP channel ..
                               ! standart output for name in the card !

        WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
     &               SNGL (WTRACK), SNGL (CMTRCK), SNGL (ATRACK)

        WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
     &                 SNGL (ZTRACK (I)), I = 0, NTRACK ),
     &               ( SNGL (DTRACK (I)), I = 1, MTRACK ),
     &                 SNGL (CTRACK)
      ENDIF

*
* --------------

* ... it is for information only ...


**      IF ( JTRACK .EQ. 8) THEN            ! --- check neutron
*      IF ( JTRACK .EQ. 8 .AND. MREG.EQ.10) THEN   ! --- check neutron
*                                                  ! after concrete
*
*         WRITE (90,*) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
*     &               SNGL (WTRACK), SNGL (CTRACK), 
*     &               SNGL (CMTRCK), SNGL (ATRACK)
*
**         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 (90,*) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
*     &                  SNGL (ZTRACK (I)) , I = 0, NTRACK )
*         WRITE (90,*) ( SNGL (DTRACK (I)), I = 1, MTRACK )
*
*         
*         WRITE (90, '(3(A,G10.3),1x)') "  EV= ",NCASE, 
*     &              '  J: ', JTRACK, ' L: ', LTRACK
*         WRITE (90, '(7(A,E10.3))') ( ' X: ', SNGL (XTRACK (I)),
*     &                 ' Y: ', SNGL (YTRACK (I)),
*     &                 ' Z: ', SNGL (ZTRACK (I)), I = 0, NTRACK ),
*     &                 ' C: ', SNGL (CTRACK)
**         WRITE (90, '(2(A, G10.3))') 'M: ', MTRACK, ' E: ', ETRACK 
*         write (90,*) 
*      END IF

*--------------------------
*
*   Region n.        1 BLKBODY 
*   Region n.        2 VOID    
*
*   Region n.        3 Rcore
*   Region n.        4 Rmod
*
*   Region n.        5 Rc1   
*   Region n.        6 Rc1   
*   Region n.        7 Rc1   
*
*   Region n.        8 Rbox    
*
*
*  so, we need particles at Regions =  5 - 6 - 7
*                --------
* =================================================

* it is for information only ...

*      IF (MREG .EQ. 10) THEN       ! select W foil... any code...
*      IF (JTRACK.EQ.7 .OR. JTRACK.EQ.8) THEN   !  get photon and neutron
*              
*        IF ( myFIRST2 ) THEN
*           myFIRST2 = .FALSE.
*           open (unit = 91, file = 'dat-all.dat',STATUS='UNKNOWN', FORM=
*     &    'FORMATTED')
**           WRITE(91, * ) '-my- general'  
*           WRITE ( 91,*) " ..MGDRAW 91 ===>  1-st entry:  general"
*           write ( 91,*) "            for W-foil, for N and Photon .." 
*           write ( 91,*) 
*        ENDIF
*     
**        IONA = 0 
**        IONZ = 0
**        IONM = 0
**
**        CALL USRDCI( JTRACK, IONA, IONZ, IONM)   ! - check ions ...
**
**        WRITE(91, '(A5,2x,i3,2x,A5,2x,i3,2x,A5,i3)')
**     &         'ion=',IONA,'Z = ',IONZ,'flag=',IONM          
*
*
*        WRITE(91, '(A,2x,i10,2x,A,2x,i5,2x,I5)')
*     &      " ....Nev= ",NCASE, "JTRACK, Ntrack= ",NTRACK,JTRACK
*        WRITE(91, "(A, E10.3, 2x E10.3 )")
*     &   "                       Ltot=, Time=  ", CMTRCK, ATRACK
*   
*
*        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, '(A,2x,D10.3,2x,D10.3,2x,3(D12.3,1x))')
*     &   '  -M - Ekin - L - len - time: ',
*     &       AM(JTRACK), ETRACK-AM(JTRACK), TRKL, CMTRCK, ATRACK           
*
*      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)

* it is for information only ...

**      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,*) " 92===> 1-st entry: BXRRAW - boundary cross."
**            WRITE(92,*) "        W-foil, for N and Photon ..."
**            WRITE(92,*) "  "
**         ENDIF
**
**         WRITE(92,"(A,3(I10,1x) )")
**     &         " B: Nev, Ntrack, Jtrack = ",NCASE,NTRACK,JTRACK 
**         WRITE(92, "(A, E10.3, 2x E10.3 )")
**     &         "     Ltot=, Time=  ", CMTRCK, ATRACK
**
**         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,'(A,2x,F7.3,1x,E10.3,1x,E10.3)')
**     &       '-PRT-: M, E, Ekin=',AM(JTRACK),ETRACK,ETRACK - AM(JTRACK)
**           WRITE(92,'(A5,1X,3(E15.3,1x))')
**     &             '-POS-',XSCO,YSCO,ZSCO
**           WRITE(92,'(A5,1X,3(E10.3,1x))')
**     &            "-COS-",CXTRCK,CYRCK,CZTRCK
** 
****           WRITE(92,'(A,2x,F7.3,1x,E10.3,1x,E10.3),/,
****     &       A5,1X,3(E15.3,1x),/,
****     &       A5,1X,3(E10.3,1x))')
****     &       '-PRT-: M, E, Ekin=',AM(JTRACK),ETRACK,ETRACK - AM(JTRACK),
****     &       '-POS-',XSCO,YSCO,ZSCO,
****     &       '-COS-',CXTRCK,CYRCK,CZTRCK
**         ENDIF
**
**      ENDIF
**      ENDIF

* -----------------------------------
*
*                particles
*
* PROTON       1    0.9382723   1.000E+18      1      1      0      1      2212
* APROTON      2    0.9382723   1.000E+18     -1     -1      0      1     -2212
* ELECTRON     3    0.0005110   1.000E+18     -1      0      0      1        11
* POSITRON     4    0.0005110   1.000E+18      1      0      0      1       -11
* NEUTRIE      5    0.0000000   1.000E+18      0      0      1      1        12
* ANEUTRIE     6    0.0000000   1.000E+18      0      0      1      1       -12
* PHOTON       7    0.0000000   1.000E+18      0      0      0      1        22
* NEUTRON      8    0.9395656    889.          0      1      0      1      2112
* ANEUTRON     9    0.9395656    889.          0     -1      0      1     -2112
* MUON+       10    0.1056584   2.197E-06      1      0      0      1       -13
* MUON-       11    0.1056584   2.197E-06     -1      0      0      1        13
* KAONLONG    12    0.4976710   5.170E-08      0      0      0      1       130
* PION+       13    0.1395699   2.603E-08      1      0      0      1       211
* PION-       14    0.1395699   2.603E-08     -1      0      0      1      -211
* KAON+       15    0.4936460   1.237E-08      1      0      0      1       321
* KAON-       16    0.4936460   1.237E-08     -1      0      0      1      -321
* LAMBDA      17    1.1156300   2.632E-10      0      1      0      1      3122
* ALAMBDA     18    1.1156300   2.632E-10      0     -1      0      1     -3122
* KAONSHRT    19    0.4976710   8.922E-11      0      0      0      1       310
* SIGMA-      20    1.1974300   1.479E-10     -1      1      0      1      3112
* SIGMA+      21    1.1893700   7.990E-11      1      1      0      1      3222
* SIGMAZER    22    1.1925500   7.400E-20      0      1      0      1      3212
* PIZERO      23    0.1349764   8.400E-17      0      0      0      1       111
* KAONZERO    24    0.4976710    0.00          0      0      0      1       311
*
*
*
* Generalised particle   Number
*
* ALL-PART                 201
* ALL-CHAR                 202
* ALL-NEUT                 203
* ALL-NEGA                 204
* ALL-POSI                 205
* NUCLEONS                 206
* NUC&PI+-                 207
* ENERGY                   208
* PIONS+-                  209
* BEAMPART                 210
* EM-ENRGY                 211
* MUONS                    212
* E+&E-                    213
* AP&AN                    214
* KAONS                    215
* STRANGE                  216
* KAONS+-                  217
* HAD-CHAR                 218
* FISSIONS                 219
* HE-FISS                  220
* LE-FISS                  221
* NEU-BALA                 222
* HAD-NEUT                 223
* KAONS0                   224
*
*
* --------------------------------------------
*  select C1-3     Regions =  5 - 6 - 7
*   -------------------------------------------
*
* ------------------------------------------------
*  error : ind < -6 for particle !!!
*
*        The transported particle is an heavy ion, carrying partile index=-7. 
*      Its properties are NOT stored in the same arrays as for "standard" 
*      particles (JTRACK>=-6)
*      The currently TRANSPORTED ion has JTRACK=-2 or JTRACK <-6
*      (if approx ion transport) . Properties can be retrieved from
*          - PAPROP with index –2 if JTRACK=-2
*       or
*          - FHEAVY with index -JTRACK if JTRACK<-6
*       You can insert a check in your code, something like
*
*            IF (JTRACK .LT. -6) THEN
*              PARTMAS = AMNHEA (-JTRACK)
*            ELSE
*              PARTMAS = AM (JTRACK)
*            END IF
****
*---------------------------------------------------



      STEP = 0.5           !   Leff step for histos


      IF((NCASE.LT.2) .AND. 
     &   (NEWREG.EQ.5 .OR. NEWREG.EQ.6 .OR. NEWREG.EQ.7)) THEN

        CALL GEOR2N ( MREG, reg1, IERR1 )
        CALL GEOR2N ( NEWREG, reg2, IERR2 )

*        WRITE (90,*) MREG,reg1
*        WRITE (90,*) NEWREG,reg2,JTRACK, NCASE

        WRITE (90,*) JTRACK, MREG,reg1, " --  ",NEWREG,reg2  
                                               ! ok, it works
        write (90,*) 
      ENDIF
***

        PMAS= 0.0                   ! for heavy part.
        IF (JTRACK .LT. -6) THEN
           PMAS = AMNHEA (-JTRACK)
        ELSE
           PMAS = AM (JTRACK)
        END IF


* ------------------------------------------

      CALL GEOR2N ( NEWREG, reg2, IERR2 )  ! let us use names...

* -------------------------------------------

      IF (reg2.EQ."Rc1" .OR. reg2.EQ."Rc2" .OR. 
     &    reg2.EQ."Rc3" .OR. reg2.EQ."Rc4")     THEN    
                                          ! select C1, C2, C3, C4

        IF ( reg2.EQ."Rc1" )  ICH = 85                !  for C1
        IF ( reg2.EQ."Rc2" )  ICH = 86                !  for C2
        IF ( reg2.EQ."Rc3" )  ICH = 87                !  for C3
        IF ( reg2.EQ."Rc4" )  ICH = 89                !  for C4

*            -------- define particle vector onto the WALL -------
*          ------------------------------------------------
*
*    at the L= 600 cm distance delY = delZ = 50 cm
*               PT= sqrt(50^2 + 50^2) = 70.71 cm
*    then   tan(T) = PT/L = 70.71 / 600 = 0.11785  = max
*           tan(T) =        50/    / 600 = 0.08333 = min = WMIN
*
cc        RMIN = 50.0          ! +- 50 cm 
cc        RMIN = 75.0          ! +- 75 cm  --> 150 cm length /50 cm window 
cc
cc        RR = 600.0/ABS(CXTRCK)              !  scale to get WALL
cc
cc        WY=  ABS(CYTRCK * RR + YSCO)        !  Y at the wall
cc       WZ=  ABS(CZTRCK * RR + ZSCO)        !  Z at the wall
cc
cc        IF ( reg2.EQ."Rc3" ) THEN
cc        IF ( reg2.EQ."Rc3" .OR. reg2.EQ."Rc4" ) THEN
cc*          ... write ALL particles here !
cc*
cc          WRITE(ICH,"( 2(I3,2x),I7,1x,4(E12.5,2x),E10.3,2x,3(E12.4,2x),
cc     &               3(E15.8,2x))")
cc     &             NEWREG, JTRACK, NCASE, Ekin, CTRACK, CMTRCK, Tmks,
cc     &             WTRACK, XSCO, YSCO, ZSCO,
cc     &             CXTRCK, CYTRCK, CZTRCK
cc        END IF

*  .. and for Rc1, Rc2 :
*
cc        IF ( reg2.EQ."Rc1" .OR. reg2.EQ."Rc2") THEN
cc        IF ( (WY .GT. RMIN) .OR. (WZ .GT. RMIN) ) GOTO 100
cc                                            ! kill this particle
cc  ====>  let us write all particles !
cc

*
*  ...check some tracks ...
*  --------------------------------------
*
cc        IF (NCASE.LT.3) THEN
cc           WRITE (90,"(A, 3(1x,E10.3), 3(2x,E12.4), 2(2x,E10.3))") 
cc     &      "track- ", XSCO,YSCO,ZSCO, CXTRCK,CYTRCK,CZTRCK, WY,WZ
cc        ENDIF

cc        PMAS= 0.0                   ! for heavy part.
cc        IF (JTRACK .LT. -6) THEN
cc           PMAS = AMNHEA (-JTRACK)
cc        ELSE
cc           PMAS = AM (JTRACK)
cc        END IF



**        Ekin = ( ETRACK-AM(JTRACK) ) * 1E+09       !  eV 
**        Tmks = ATRACK * 1E+06                      ! mks 
cc        Ekin = ( ETRACK-AM(JTRACK) )               !  GeV 
        Ekin = ( ETRACK- PMAS )               !  GeV 
        IF (Ekin .LT. 1E-14) THEN
           Ekin = 2E-14
            WRITE (90,*) " ============= ev= ",NCASE, " Ekin= ",Ekin
        END IF

        Tmks = ATRACK                              ! sec 

        


***        IF ( JTRACK.EQ.201 ) THEN      !  get  ALL particles
*        IF ( JTRACK.EQ.7 .OR. JTRACK.EQ.8 ) THEN      !  get  G  and N
**          WRITE(85, " (I7,1x,7(E12.4,1x)) ")
***          WRITE(85, " (I5,2x,I7,1x,7(E12.4,1x)) ")
***     &             NEWREG, NCASE, Ekin, CMTRCK, Tmks,
***     &             WTRACK, XSCO, YSCO, ZSCO 
**          WRITE(85, " (I5,2x,I7,1x,10(E12.4,1x)) ")
          WRITE(ICH,"( 2(I3,2x),I7,1x,4(E12.5,2x),E10.3,2x,3(E12.4,2x),
     &               3(E15.8,2x))")
     &             NEWREG, JTRACK, NCASE, Ekin, CTRACK, CMTRCK, Tmks,
     &             WTRACK, XSCO, YSCO, ZSCO,
     &             CXTRCK, CYTRCK, CZTRCK
***        END IF                                 ! G and N
*  
cc      END IF           ! detector-region
      END IF         
* --------------------------------------------

cc 100  CONTINUE                       ! particle is out of the window

*
* --------------------------------------------
*  select C2
*   -------------------------------------------
*
        IF (reg2.EQ."Rc2" .AND. JTRACK.EQ.8) THEN
                                                    ! sel. N in C2
                                                    ! to write L,t...

**        Ekin = ( ETRACK-AM(JTRACK) ) * 1E+09       !  eV 
**        Tmks = ATRACK * 1E+06                      ! mks 
cc        Ekin = ( ETRACK-AM(JTRACK) )               !  GeV 
        Ekin = ( ETRACK- PMAS )               !  GeV 

        IF (Ekin .LT. 1E-14) THEN
           Ekin = 1E-14
cc            WRITE (90,*) " ************** ev= ",NCASE, " Ekin= ",Ekin
        END IF

        Tmks = ATRACK                              ! sec 

**        IF ( JTRACK.EQ.8 ) THEN      !  get  N
**          WRITE(ICH,"( 2(I3,2x),I7,1x,4(E12.5,2x),E10.3,2x,3(E12.4,2x),
**     &               3(E15.8,2x))")
**     &             NEWREG, JTRACK, NCASE, Ekin, CTRACK, CMTRCK, Tmks,
**     &             WTRACK, XSCO, YSCO, ZSCO,
**     &             CXTRCK, CYTRCK, CZTRCK
**        END IF                                 ! G and N


*
*  ..... now fill histos for delayed N in moderator : 
*        ============================================
*
*           del L = V * t 
*        
*          take V from  formula :   E= m*V2/2 ==> V= sqrt( 2*E / m) 
*
*          vv= sqrt(2.0 * emev / amn) * vnorm         // cm/mks
*              vnorm = 1E+05        // for V-speed: m/ns -> cm/mks
*

*   E-bins:
*
*N, Ei= 1-4 ev
*N, Ei= 4-15 ev      n vs delL(cm) ,  100, 0.0,100.0)
*N, Ei= 15-50 ev,
*N, Ei= 50-120 ev,
*N, Ei= 120-300 ev,
*N, Ei= 300-1000 ev,
*N, Ei= 1000-3000 ev, 
*N, Ei= 3000-7000 ev,
* N, Ei= 7000-15 000 ev,
*N, Ei= 15 10^3- 3 10^4 ev
* N, Ei= 3 10^4 - 6*10^4,
*N, Ei= 6 10^4 -  10^5 ev,
*N, Ei= 10^5 - 2 10^5 ev,
*N, Ei= 2 10^5 - 4 10^5 ev,
*
**
***  En= 0.025 eV --> 2200 m/s -> from books
***
***  vv = sqrt(2.0*0.025*1E-06 /  939.565)  ... i.e MeV / MeV
***   vv= 7.2949377e-06    then  == 2200 m/s = 22E+05/1E+06 cm/mks = 0.22 cm/mks
***   so k = 0.22 / 7.2949377e-06 =  30157.9
***  
***  c = 299792458 м/с = 2.998E+08 m/s
***  then   sqrt(2*0.0253E-06 / 939.565) *c = 
***          sqrt(2*0.0253E-06 / 939.565) * 2.998E+08 = 2200 m/s  !
***
***  V = 1.38429E-03 * sqrt(E) , cm/ns , E in eV
***
***  v = 1.3842933 sqrt(E) cm/mks, E in eV 
***         1.3842933 * sqrt(0.0253) = 2202 m/s
***



      ev =   Ekin * 1E+09            ! E , eV
      ekev = Ekin * 1E+06            ! E , KeV
      emev = Ekin * 1E+03            ! E , MeV

      ttt = ATRACK * 1E+06           ! T, mks

      amn= 939.565                   ! MeV
cc      vnorm = 1E+05                  ! for V-speed: m/ns -> cm/mks
      vnorm = 0.2998E+05                  ! for V-speed: m/ns -> cm/mks


cc      WRITE (IODRAW,*) ".. ev= ",NCASE,JTRACK," Ekin, ETRACK= ",Ekin,
cc     &      SNGL (ETRACK),emev,amn
cc      WRITE (90,*) ".. ev= ",NCASE,JTRACK," Ekin, ETRACK= ",Ekin,
cc     &      SNGL (ETRACK),emev,amn

c      PRINT *, ".. ev= ",NCASE,JTRACK," Ekin, ETRACK= ",Ekin,
c     &      SNGL (ETRACK),emev,amn

      vv= sqrt(2.0 * emev / amn) * vnorm      ! cm/mks


      DO I = 1, 14
**        IF ((emev .GE. EC(I,1)) .AND. (emev .LT. EC(I,2))) THEN
        IF ((ev .GE. EC(I,1)) .AND. (ev .LT. EC(I,2))) THEN
           DEL= vv * ttt                ! del L definition for target           
*           AIND = DEL / 100.
*           IND2 = AIND * 100
*           IND = INT( DEL / 100.)            ! number of bins
*          
*           write(90, "(A,3(2x,I8))") "I,IND,IND2= ", I, IND, IND2
*           write(90, "(A,3(2x,E10.3))") "DEL, V,T = ",DEL,vv,ttt 

cc           IND= INT(DEL)
**           IND= INT(DEL/2.)            ! L bin = 2 cm
           IND=  DEL / STEP


*           IF(IND .LT. 1) CONTINUE
*           IF(IND .GT. 100) CONTINUE

*           IF ( (IND .GE.1)  .AND. (IND.LE.100)) THEN
cc           IF ( (IND .GE.1)  .AND. (IND.LE.200)) THEN
cc             HH(I,1,IND) = HH(I,1,IND) + 1
cc             HH(I,2,IND) = HH(I,2,IND) + emev        ! MeV
cc             HH(I,3,IND) = HH(I,3,IND) + ttt         ! mks
cc           END IF
           IF ( (IND .GT.0)  .AND. (IND.LE.400)) THEN
             HH(I,IND) = HH(I,IND) + 1
           END IF


        END IF
      ENDDO

* ------------------------------------------------
*  
      END IF           ! detector-region
* --------------------------------------------


*
* <==== ///////////// -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,A5,1X,D10.3,/,A5,1X,D10.3)') 
**     &         'END = ',JTRACK,  'EDP= ', EDEP, 
**     &                           'TRK= ', TRKL
***        END IF
**      END IF

*      EDEP=0               ! update values
*      TRKL=0

*
* print event number here
*
**      IF (NCASE .LT. 100) THEN
      IF (NCASE .LT. 3) THEN
         WRITE(90,"(A,2x,I10,2x,I10)")
     &        "*** EEDRAW - END of EVENT : Ncase=", NCASE

      ELSE IF ( MOD (NCASE, 1000) .EQ. 0)  THEN
              WRITE(90,"(A,2x,I10,2x,I10)")
     &        "*** EEDRAW - END of EVENT : Nevet= ", NCASE
              WRITE(90,'(A32)') RUNTIM
      END IF  

*
*  type my histos :
*
*      IF ( MOD (NCASE, 5000) .EQ. 0) THEN
*      IF ( MOD (NCASE, 20000) .EQ. 0) THEN
*      IF ( MOD (NCASE, 10000) .EQ. 0) THEN
cc      IF ( MOD (NCASE, 1000) .EQ. 0) THEN
      IF ( MOD (NCASE, NCASES) .EQ. 0) THEN

          write (88,*) "  "
          write (88,*) "    NCASE =  ",NCASE
          write (88,*) "    ================  "


**        DO I= 1,14
**
**          write (88,*) "  "
***          write (88,*) " E1 - E2 : ", EC(I,1), EC(I,2)
**          write (88,"(A,2(2x,E8.3))") " E1 - E2 : ", EC(I,1), EC(I,2)
**          write (88,*) "  "
**         
***          DO J= 1,3
**          DO J= 1,1
**
**            IF(J.EQ.1) THEN
**              write (88,*) "......... hist with ENTRY ...."
**              DO K=1,100 
***                 write (88,*) K, HH(I,J,K)
**                 write (88,"(I6, 1(2x,E8.3))") K, HH(I,J,K)
**              ENDDO
**            END IF
**
**            IF(J.EQ.2) THEN
**              write (88,*) "......... hist with <E> ...."
**              DO K=1,100 
**                 srn= 0.0
**                 if(HH(I,1,K) .GE. 1) srn= HH(I,J,K)/HH(I,1,K)
**                 write (88,*) K, srn
**              ENDDO
**            END IF
**
**            IF(J.EQ.3) THEN
**               write (88,*) "......... hist with <T> ...."
**               DO K=1,100 
**                 srn= 0.0
**                 if(HH(I,1,K) .GE. 1) srn= HH(I,J,K)/HH(I,1,K)
**                 write (88,*) K, srn
**               ENDDO
**             END IF
**
**          ENDDO              ! J = 3
**
**        ENDDO                 ! I = 14

*
*     ====  type info as a table here ==========
*

        write (88,*) "  "
        DO I= 1,14
           write (88,"(I3,2x,A,2(2x,E8.3))") 
     &                I," E1 - E2 : ", EC(I,1), EC(I,2)
        ENDDO
        write (88,*) "  "

        write (88,*) "  "
        write (88,*) "  ===Entries into N vs L, (bin= 2 cm)"
        write (88,*) "  "
        write (88,*) "   BIN  - L,cm  - delE=1  -   delE=2  - ",
     &              "  delE=3 ",
     &      " -  delE=4   -   delE=5  -   delE=6  -   delE=7   - ",
     &      "  delE=8  -  delE=9  -  delE=10 - ",
     &      "   delE=11  -  delE=12  -  delE=13  -  delE=14"
        write (88,*) "  "


cc        DO I = 1, 200                   ! 200 bins
        DO I = 1, 400                   ! 200 bins
**        DO I = 1, 100
*          DO J = 1,7
*             write (88,"(I6, 1(2x,E8.3))") I, I*2, HH(J,1,I)
*          ENDDO

cc           write (88,"(2I6,2x, 14(2x,E10.3))") 
cc     &           I, I, HH(1,1,I), HH(2,1,I),
cc     &           HH(3,1,I), HH(4,1,I), HH(5,1,I), HH(6,1,I), HH(7,1,I),
cc     &         HH(8,1,I), HH(9,1,I), HH(10,1,I), HH(11,1,I), HH(12,1,I),
cc     &         HH(13,1,I), HH(14,1,I)

           write (88,"(I6,2x, f9.2, 14(2x,E10.3))") 
cc     &         I, HL(I,1), HH(1,I), HH(2,I), HH(3,I), HH(4,I), HH(5,I),
     &         I, 0.0 + I*STEP, 
     &                 HH(1,I), HH(2,I), HH(3,I), HH(4,I), HH(5,I),
     &                 HH(6,I), HH(7,I), HH(8,I), HH(9,I), HH(10,I),
     &                 HH(11,I), HH(12,I), HH(13,I), HH(14,I)

*           write (88,"(2I6,2x, 14(2x,E10.3))") 
*     &           I, I*2, HH(1,1,I), HH(2,1,I),
*     &           HH(3,1,I), HH(4,1,I), HH(5,1,I), HH(6,1,I), HH(7,1,I),
*     &         HH(8,1,I), HH(9,1,I), HH(10,1,I), HH(11,1,I), HH(12,1,I),
*     &         HH(13,1,I), HH(14,1,I)

*           write (88,"(2I6, 7(2x,E8.3))") I, I*2, HH(8,1,I), HH(9,1,I),
*     &        HH(10,1,I), HH(11,1,I), HH(12,1,I), HH(13,1,I), HH(14,1,I)
        ENDDO


      END IF          !  end of N-cycle

*
*  <=====/////////////  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
**  |
**  +-------------------------------------------------------------------*

*
* print event number here ... COMMON / FLKSTK /
*
**      IF (NCASE .LT. 100) THEN
*      IF (NCASE .LT. 10) THEN
*         WRITE(95,'(A,2x,I10,2x,I10,/)') 
*     &          " -- SODRAW - start beam part -- , NCASE= ",NCASE
*
*      ELSE IF ( MOD (NCASE, 1000) .EQ. 0)  THEN
*         WRITE(95,'(A,2x,I10,2x,I10,/)') 
*     &          " -- SODRAW - start beam part -- , NCASE= ",NCASE
*   
*      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 ( .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


*
*            print out interaction coordinates:
*    Icode = 101 :  inelastic interaction secondaries  
*

*        IF (ICODE .EQ. 101) THEN
*
***        IF( NCASE .LT. 100) THEN
*        IF( NCASE .LT. 10 .AND. JTRACK.EQ.1) THEN   ! proton only
*
*        IF ( myFIRST ) THEN
*           myFIRST = .FALSE.
*           open (unit = 94, file = 'Inter.dat',STATUS='UNKNOWN', FORM=
*     &    'FORMATTED')
*           WRITE(94,*) " ===> 1-st entry: USDRAW - user"
*        ENDIF
*
***        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
*        WRITE (94, "(A, 2(E10.3, 2x), E10.3 )")
*     &   "    Ltot=, Time=  ", CTRACK, CMTRCK, ATRACK
*
*        ENDIF
*
*        ENDIF        ! icode
*


*
*  ...for information only
*

**      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
*8**              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
**        WRITE (94, "(A, E10.3, 2x E10.3 )")
**     &   "    Ltot=, Time=  ", CMTRCK, ATRACK

*
* 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))") 
**           WRITE ( 94,"(A,2x,5(I5,1x),6(F10.3,1x))")
**     &           'Secon-- ',NCASE,JTRACK,ip,ICODE,ZPART,
**     &          APART,TKI(ip),CXR(ip),CYR(ip),CZR(ip),WEI(ip)
**           WRITE( 94, "(A, E10.3, 2x E10.3 )")
**     &             "    Ltot=, Time=  ", CMTRCK, ATRACK
**
**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


