*$ CREATE MGDRAW.FOR
*COPY MGDRAW
*                                                                      *
*=== mgdraw ===========================================================*
*                                                                      *
      SUBROUTINE MGDRAW ( ICODE, MREG )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2006      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  05-may-06       by        Alfredo Ferrari           *
*                                              INFN - Milan            *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(CASLIM)'
      INCLUDE '(COMPUT)'
      INCLUDE '(SOURCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(GENSTK)'
      INCLUDE '(MGDDCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(QUEMGD)'
      INCLUDE '(SUMCOU)'
      INCLUDE '(TRACKR)'
      INCLUDE '(SOUEVT)'
*
      DIMENSION DTQUEN ( MXTRCK, MAXQMG )
*
      CHARACTER*20 FILNAM
      LOGICAL LFCOPE
      SAVE LFCOPE
      DATA LFCOPE / .FALSE. /
      DIMENSION XTUP(10)
      REAL EDEPCON(6,2000)
      SAVE EDEPCON
      DATA EDEPCON /12000*0.0D0/
      REAL EDEPPOI(6,2000)
      SAVE EDEPPOI
      DATA EDEPPOI/12000*0.0D0/
      REAL ENERIN(2000)
      SAVE ENERIN
      DATA ENERIN/2000*0.0D0/
      REAL ENEROUT(2000)
      SAVE ENEROUT
      DATA ENEROUT/2000*0.0D0/
      REAL MOMINAR(12,2000)
      SAVE MOMINAR
      DATA MOMINAR/24000*0.0D0/
      REAL MOMOUTAR(12,2000)
      SAVE MOMOUTAR
      DATA MOMOUTAR/24000*0.0D0/
      INTEGER MOMINCO
      SAVE MOMINCO
      DATA MOMINCO/0/
      INTEGER MOMOUTCO
      SAVE MOMOUTCO
      DATA MOMOUTCO/0/
      REAL MOMTARAR(12,2000)
      SAVE MOMTARAR
      DATA MOMTARAR/24000*0.0D0/
      INTEGER MOMTARCO
      SAVE MOMTARCO
      DATA MOMTARCO/0/
      INTEGER TRE
      
      INTEGER IONA, IONZ, IONM
      
      INTEGER BESEEN
      INTEGER UNIQMAX
      SAVE UNIQMAX
      DATA UNIQMAX/0/
      INTEGER UNIQ(2000)
      SAVE UNIQ
      DATA UNIQ/2000*0.0D0/
*
*----------------------------------------------------------------------*
*                                                                      *
*     Icode = 1: call from Kaskad                                      *
*     Icode = 2: call from Emfsco                                      *
*     Icode = 3: call from Kasneu                                      *
*     Icode = 4: call from Kashea                                      *
*     Icode = 5: call from Kasoph                                      *
*                                                                      *
*----------------------------------------------------------------------*
*                                                                      *
      IF ( .NOT. LFCOPE ) THEN
         LFCOPE = .TRUE.
         IF ( KOMPUT .EQ. 2 ) THEN
            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
         ELSE
            FILNAM = CFDRAW
         END IF
         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
     &          'UNFORMATTED' )
      END IF
      
      IF ((JTRACK.LE.-2).AND.(MREG.EQ.4).AND.(MTRACK.GT.0)) THEN
          BESEEN=-1
	  DO I=1,UNIQMAX
	     IF(UNIQ(I).EQ.ISPUSR(MKBMX2)) THEN
	        BESEEN=I
	     END IF
	  END DO

	  IF (BESEEN.EQ.-1) THEN
	     UNIQMAX=UNIQMAX+1
	     UNIQ(UNIQMAX)=ISPUSR(MKBMX2)
	     BESEEN=UNIQMAX
	  END IF

	  IF (JTRACK.GE.-6) THEN
	  CALL USRDCI(JTRACK,IONA,IONZ,IONM)
	  ELSE
	  IONA=IBHEAV(-JTRACK)
	  IONZ=ICHEAV(-JTRACK)
	  IONM=IMHEAV(-JTRACK)
	  END IF

          DO I=1,MTRACK
             EDEPCON(4,BESEEN)=EDEPCON(4,BESEEN)+DTRACK(I)
          END DO
            EDEPCON(1,BESEEN)=IONZ
            EDEPCON(2,BESEEN)=IONA
            EDEPCON(3,BESEEN)=IONM
	    EDEPCON(5,BESEEN)=ISPUSR(MKBMX2)
	    EDEPCON(6,BESEEN)=ISPUSR(1)
      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
      IF ( LQEMGD ) THEN
         RULLL = ZERZER
         CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
*         WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)), I = 1, MTRACK ),
*     &                      JBK = 1, NQEMGD )
      END IF
*  |  End of quenching
*  +-------------------------------------------------------------------*
      RETURN
*
*======================================================================*
*                                                                      *
*     Boundary-(X)crossing DRAWing:                                    *
*                                                                      *
*     Icode = 1x: call from Kaskad                                     *
*             19: boundary crossing                                    *
*     Icode = 2x: call from Emfsco                                     *
*             29: boundary crossing                                    *
*     Icode = 3x: call from Kasneu                                     *
*             39: boundary crossing                                    *
*     Icode = 4x: call from Kashea                                     *
*             49: boundary crossing                                    *
*     Icode = 5x: call from Kasoph                                     *
*             59: boundary crossing                                    *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
*      PRINT *,"FROM THE BXDRAW ROUTINE"
*      PRINT *,MREG,NEWREG,JTRACK,ETRACK,AM(JTRACK)

      IF(MREG.EQ.3.AND.NEWREG.EQ.2.AND.JTRACK.LE.-2) THEN
         IF( JTRACK.EQ.1) THEN
            IF(ETRACK.GT.AM(JTRACK)) THEN
               XTUP(1) = JTRACK
               XTUP(2) = ETRACK-AM(JTRACK)
               XTUP(3) = XSCO
               XTUP(4) = YSCO
               XTUP(5) = ZSCO
               XTUP(6) = CXTRCK
               XTUP(7) = CYTRCK
               XTUP(8) = CZTRCK
*               PRINT *,"BXDRAW: Boundary crossing at event: ",NCASE
*               CALL treefill(JTRACK,
*     +                       (ETRACK-AM(JTRACK)),XSCO,YSCO,ZSCO,
*     +                       CXTRCK,CYTRCK,CZTRCK)
*               DO I=1,Npsoev
*      PRINT *,"Xpos=",Xsoevt(I),"  Ypos=",Ysoevt(I),"  Zpos=",Zsoevt(I)
*      PRINT *,"Xdir=",TXsoev(I),"  Ydir=",TYsoev(I),"  Zdir=",TZsoev(I)
*               END DO

            ENDIF
         ENDIF
	 IF (JTRACK.GE.-6) THEN
	 CALL USRDCI(JTRACK,IONA,IONZ,IONM)
	 ELSE 
         IONA=IBHEAV(-JTRACK)
	 IONZ=ICHEAV(-JTRACK)
	 IONM=IMHEAV(-JTRACK)
	 END IF
	 MOMTARCO=MOMTARCO+1
	 MOMTARAR(1,MOMTARCO)=IONZ
	 MOMTARAR(2,MOMTARCO)=IONA
	 MOMTARAR(3,MOMTARCO)=IONM
	 MOMTARAR(4,MOMTARCO)=Ptrack/IONA
	 MOMTARAR(5,MOMTARCO)=Cxtrck
	 MOMTARAR(6,MOMTARCO)=Cytrck
	 MOMTARAR(7,MOMTARCO)=Cztrck
	 MOMTARAR(8,MOMTARCO)=Xtrack(Mtrack)
	 MOMTARAR(9,MOMTARCO)=Ytrack(Mtrack)
	 MOMTARAR(10,MOMTARCO)=Ztrack(Mtrack)
	 MOMTARAR(11,MOMTARCO)=ISPUSR(MKBMX2)
	 MOMTARAR(12,MOMTARCO)=ISPUSR(1)
      ENDIF
      IF (MREG.EQ.2.AND.NEWREG.EQ.4.AND.JTRACK.LE.-2) THEN
	 IF (JTRACK.GE.-6) THEN
	 CALL USRDCI(JTRACK,IONA,IONZ,IONM)
	 ELSE 
         IONA=IBHEAV(-JTRACK)
	 IONZ=ICHEAV(-JTRACK)
	 IONM=IMHEAV(-JTRACK)
	 END IF
*	 ENERIN(IONZ)=ENERIN(IONZ)+ETRACK-AM(JTRACK)
	 MOMINCO=MOMINCO+1
	 MOMINAR(1,MOMINCO)=IONZ
	 MOMINAR(2,MOMINCO)=IONA
	 MOMINAR(3,MOMINCO)=IONM
	 MOMINAR(4,MOMINCO)=Ptrack/IONA
	 MOMINAR(5,MOMINCO)=Cxtrck
	 MOMINAR(6,MOMINCO)=Cytrck
	 MOMINAR(7,MOMINCO)=Cztrck
	 MOMINAR(8,MOMINCO)=Xtrack(Mtrack)
	 MOMINAR(9,MOMINCO)=Ytrack(Mtrack)
	 MOMINAR(10,MOMINCO)=Ztrack(Mtrack)
	 MOMINAR(11,MOMINCO)=ISPUSR(MKBMX2)
	 MOMINAR(12,MOMINCO)=ISPUSR(1)
*	 ELSE
*	    PRINT *,"HUGE PROBLEM!!!not an ion!!",JTRACK
      END IF
      IF (MREG.EQ.4.AND.NEWREG.EQ.2.AND.JTRACK.LE.-2) THEN
	 IF (JTRACK.GE.-6) THEN
	 CALL USRDCI(JTRACK,IONA,IONZ,IONM)
	 ELSE 
         IONA=IBHEAV(-JTRACK)
	 IONZ=ICHEAV(-JTRACK)
	 IONM=IMHEAV(-JTRACK)
	 END IF
*	 ENEROUT(IONZ)=ENEROUT(IONZ)+ETRACK-AM(JTRACK)
	 MOMOUTCO=MOMOUTCO+1
	 MOMOUTAR(1,MOMOUTCO)=IONZ
	 MOMOUTAR(2,MOMOUTCO)=IONA
	 MOMOUTAR(3,MOMOUTCO)=IONM
	 MOMOUTAR(4,MOMOUTCO)=Ptrack/IONA
	 MOMOUTAR(5,MOMOUTCO)=Cxtrck
	 MOMOUTAR(6,MOMOUTCO)=Cytrck
	 MOMOUTAR(7,MOMOUTCO)=Cztrck
	 MOMOUTAR(8,MOMOUTCO)=Xtrack(Mtrack)
	 MOMOUTAR(9,MOMOUTCO)=Ytrack(Mtrack)
	 MOMOUTAR(10,MOMOUTCO)=Ztrack(Mtrack)
	 MOMOUTAR(11,MOMOUTCO)=ISPUSR(MKBMX2)
	 MOMOUTAR(12,MOMOUTCO)=ISPUSR(1)
      END IF
      
*      IF (MREG.EQ.3.AND.NEWREG.EQ.2) THEN   
*         IF (JTRACK.GT.-2) THEN
*            WTRACK=ZERZER
*	  ELSE 
*	  CALL USRDCI(JTRACK,IONA,IONZ,IONM)
*	  PRINT *,"IONA=",IONA," IONZ=",IONZ," IONM=",IONM
*         END IF
*      END IF 
      RETURN
*
*======================================================================*
*                                                                      *
*     Event End DRAWing:                                               *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY EEDRAW ( ICODE )
*      PRINT *,"EEDRAW:",Ncase," number of ions",UNIQMAX
*      PRINT *,"EEDRAW:",Ncase," unique array:"
*      DO I=1,UNIQMAX
*      PRINT *,"UNIQ[",I,"]=",UNIQ(I)
*      END DO
*      PRINT *,"EEDRAW:",Ncase," continious deposition:",EDEPCON
*      DO J=1,UNIQMAX
*      DO I=1,5
*      PRINT *,"EDEPCON[",I,"][",J,"]=",EDEPCON(I,J)
*      END DO
*      END DO
*      PRINT *,"EEDRAW:",Ncase," point deposition:",EDEPPOI
*      PRINT *,"EEDRAW:",Ncase," input energy:",ENERIN
*      PRINT *,"EEDRAW:",Ncase," output energy:",ENEROUT
      PRINT *,"EEDRAW:",Ncase," Number of particles tout:",MOMTARCO
      PRINT *,"EEDRAW:",Ncase," array of particles tout:",MOMTARAR
      PRINT *,"EEDRAW:",Ncase," Number of particles din:",MOMINCO
      PRINT *,"EEDRAW:",Ncase," Array of particles din:",MOMINAR
      PRINT *,"EEDRAW:",Ncase," Number of particles dout:",MOMOUTCO
      PRINT *,"EEDRAW:",Ncase," Array of particles dout:",MOMOUTAR
*      CALL targetpt(MOMTARCO,MOMTARAR)
*      CALL degptin(MOMINCO,MOMINAR)
*      CALL degptout(MOMOUTCO,MOMOUTAR)
*      CALL degencon(UNIQMAX,EDEPCON)
*      CALL degenpoi(UNIQMAX,EDEPPOI)
*      CALL treefill()

*      DO I=1,MOMTARCO
*      PRINT *,MOMTARCO,MOMTARAR
*      END DO

      
      UNIQMAX=0
      DO I=1,200
      UNIQ(I)=0
      END DO  
  
      DO I=1,2000
      ENERIN(I)=0.0D0
      ENEROUT(I)=0.0D0
      END DO
      
      DO J=1,2000
      DO I=1,6
      EDEPCON(I,J)=0.0D0
      EDEPPOI(I,J)=0.0D0
      END DO
      END DO
      
      MOMINCO=0
      MOMOUTCO=0
      DO J=1,2000
         DO I=1,12
            MOMINAR(I,J)=0.0D0
            MOMOUTAR(I,J)=0.0D0
	 END DO
      END DO
      MOMTARCO=0
      DO J=1,2000
         DO I=1,12
            MOMTARAR(I,J)=0.0D0
	 END DO
      END DO
      PRINT *,"N=",Ncase
      TRE=Ncase-(INT(Ncase/1000)*1000)
      IF (TRE.EQ.0) THEN
*      CALL saveauto()
      END IF
      RETURN
*
*======================================================================*
*                                                                      *
*     ENergy deposition DRAWing:                                       *
*                                                                      *
*     Icode = 1x: call from Kaskad                                     *
*             10: elastic interaction recoil                           *
*             11: inelastic interaction recoil                         *
*             12: stopping particle                                    *
*             13: pseudo-neutron deposition                            *
*             14: escape                                               *
*             15: time kill                                            *
*     Icode = 2x: call from Emfsco                                     *
*             20: local energy deposition (i.e. photoelectric)         *
*             21: below threshold, iarg=1                              *
*             22: below threshold, iarg=2                              *
*             23: escape                                               *
*             24: time kill                                            *
*     Icode = 3x: call from Kasneu                                     *
*             30: target recoil                                        *
*             31: below threshold                                      *
*             32: escape                                               *
*             33: time kill                                            *
*     Icode = 4x: call from Kashea                                     *
*             40: escape                                               *
*             41: time kill                                            *
*             42: delta ray stack overflow                             *
*     Icode = 5x: call from Kasoph                                     *
*             50: optical photon absorption                            *
*             51: escape                                               *
*             52: time kill                                            *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO )
      IF ((JTRACK.EQ.-2).AND.(MREG.EQ.4)) THEN
          BESEEN=-1
	  DO I=1,UNIQMAX
	     IF(UNIQ(I).EQ.ISPUSR(MKBMX2)) THEN
	        BESEEN=I
	     END IF
	  END DO
	  IF (BESEEN.EQ.-1) THEN
	     UNIQMAX=UNIQMAX+1
	     UNIQ(UNIQMAX)=ISPUSR(MKBMX2)
	     BESEEN=UNIQMAX
	  END IF
	  CALL USRDCI(JTRACK,IONA,IONZ,IONM)
          EDEPPOI(4,BESEEN)=EDEPPOI(4,BESEEN)+RULL
          EDEPCON(1,BESEEN)=IONZ
          EDEPCON(2,BESEEN)=IONA
          EDEPCON(3,BESEEN)=IONM
	  EDEPCON(5,BESEEN)=ISPUSR(MKBMX2)
	  EDEPCON(6,BESEEN)=ISPUSR(1)
      END IF
      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
*  +-------------------------------------------------------------------*
      RETURN
*
*======================================================================*
*                                                                      *
*     SOurce particle DRAWing:                                         *
*                                                                      *
*======================================================================*
*
      ENTRY SODRAW
      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 )
*  |
*  +-------------------------------------------------------------------*
*  |
*      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                                    *
*             100: elastic   interaction secondaries                   *
*             101: inelastic interaction secondaries                   *
*             102: particle decay  secondaries                         *
*             103: delta ray  generation secondaries                   *
*             104: pair production secondaries                         *
*             105: bremsstrahlung  secondaries                         *
*             110: decay products                                      *
*     Icode = 20x: call from Emfsco                                    *
*             208: bremsstrahlung secondaries                          *
*             210: Moller secondaries                                  *
*             212: Bhabha secondaries                                  *
*             214: in-flight annihilation secondaries                  *
*             215: annihilation at rest   secondaries                  *
*             217: pair production        secondaries                  *
*             219: Compton scattering     secondaries                  *
*             221: photoelectric          secondaries                  *
*             225: Rayleigh scattering    secondaries                  *
*     Icode = 30x: call from Kasneu                                    *
*             300: interaction secondaries                             *
*     Icode = 40x: call from Kashea                                    *
*             400: delta ray  generation secondaries                   *
*  For all interactions secondaries are put on GENSTK common (kp=1,np) *
*  but for KASHEA delta ray generation where only the secondary elec-  *
*  tron is present and stacked on FLKSTK common for kp=npflka          *
*                                                                      *
*======================================================================*
*
      ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
      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
* No output by default:
      RETURN
*=== End of subrutine Mgdraw ==========================================*
      END

