*$ 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)'
*
      DIMENSION DTQUEN ( MXTRCK, MAXQMG )
*
      CHARACTER*20 FILNAM
      LOGICAL LFCOPE
      SAVE LFCOPE
      DATA LFCOPE / .FALSE. /
*
*----------------------------------------------------------------------*
*                                                                      *
*     Icode = 1: call from Kaskad                                      *
*     Icode = 2: call from Emfsco                                      *
*     Icode = 3: call from Kasneu                                      *
*     Icode = 4: call from Kashea                                      *
*     Icode = 5: call from Kasoph                                      *
*                                                                      *
*----------------------------------------------------------------------*
*                                                                      *
      IF ( .NOT. LFCOPE ) THEN
         LFCOPE = .TRUE.
         IF ( KOMPUT .EQ. 2 ) THEN
            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
         ELSE
            FILNAM = CFDRAW
         END IF
         WRITE(*,*) 'TRAJECTORY OPEN!'
         WRITE(*,'(A)') 'FILNAM = ',FILNAM
         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
     &          'UNFORMATTED' )
      END IF
C
C     Write trajectories of optical photons
C
      IF(JTRACK .EQ. 11) THEN
         WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
     &        SNGL (WTRACK)
         WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
     &        SNGL (ZTRACK (I)), I = 0, NTRACK ),
     &        ( SNGL (DTRACK (I)), I = 1,MTRACK ),
     &        SNGL (CTRACK)
         WRITE(IODRAW) SNGL(CXTRCK),SNGL(CYTRCK),SNGL(CZTRCK)
      ENDIF
      RETURN      

*
*======================================================================*
*                                                                      *
*     Boundary-(X)crossing DRAWing:                                    *
*                                                                      *
*     Icode = 1x: call from Kaskad                                     *
*             19: boundary crossing                                    *
*     Icode = 2x: call from Emfsco                                     *
*             29: boundary crossing                                    *
*     Icode = 3x: call from Kasneu                                     *
*             39: boundary crossing                                    *
*     Icode = 4x: call from Kashea                                     *
*             49: boundary crossing                                    *
*     Icode = 5x: call from Kasoph                                     *
*             59: boundary crossing                                    *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
      
cc      DO I=7,3,-1
cc         IF(JTRACK.EQ.8) THEN
cc              IF(MREG.EQ.I.and.newreg.eq.I-1) THEN 
cc                  IF(NPFLKA.EQ.0.0) THEN 
cc                      WRITE(50.-I,*) ETRACK-AM(JTRACK)
cc                  ENDIF
cc              ENDIF
cc	 ENDIF     
cc      ENDDO	   

      IF(MREG.EQ.7.and.newreg.eq.8) THEN 
         IF(JTRACK.EQ.11) THEN 
	      IF(ETRACK.GT.AM(JTRACK)) THEN
	          WRITE(43.0,*) mreg, ETRACK-AM(JTRACK),
     +		                XSCO,YSCO,ZSCO, NCASE, NPFLKA
c                  WRITE(41.0,*) ETRACK-AM(JTRACK)
              ENDIF
	  ENDIF
      ENDIF
*//     
*//      IF(MREG.EQ.3.AND.NEWREG.EQ.2) THEN 
*//         IF(JTRACK.EQ.11) THEN 
*//	      IF(ETRACK.GT.AM(JTRACK)) THEN
*//	          WRITE(44.0,*) mreg, ETRACK-AM(JTRACK),
*//     +		                XSCO,YSCO,ZSCO
*//c                  WRITE(41.0,*) ETRACK-AM(JTRACK)
*//              ENDIF
*//	  ENDIF
*//      ENDIF
      
      RETURN
*
*======================================================================*
*                                                                      *
*     Event End DRAWing:                                               *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY EEDRAW ( ICODE )
      RETURN
*
*======================================================================*
*                                                                      *
*     ENergy deposition DRAWing:                                       *
*                                                                      *
*     Icode = 1x: call from Kaskad                                     *
*             10: elastic interaction recoil                           *
*             11: inelastic interaction recoil                         *
*             12: stopping particle                                    *
*             13: pseudo-neutron deposition                            *
*             14: escape                                               *
*             15: time kill                                            *
*     Icode = 2x: call from Emfsco                                     *
*             20: local energy deposition (i.e. photoelectric)         *
*             21: below threshold, iarg=1                              *
*             22: below threshold, iarg=2                              *
*             23: escape                                               *
*             24: time kill                                            *
*     Icode = 3x: call from Kasneu                                     *
*             30: target recoil                                        *
*             31: below threshold                                      *
*             32: escape                                               *
*             33: time kill                                            *
*     Icode = 4x: call from Kashea                                     *
*             40: escape                                               *
*             41: time kill                                            *
*             42: delta ray stack overflow                             *
*     Icode = 5x: call from Kasoph                                     *
*             50: optical photon absorption                            *
*             51: escape                                               *
*             52: time kill                                            *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO )
ccc      IF ( .NOT. LFCOPE ) THEN
ccc         LFCOPE = .TRUE.
ccc         IF ( KOMPUT .EQ. 2 ) THEN
ccc            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
ccc         ELSE
ccc            FILNAM = CFDRAW
ccc         END IF
ccc         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
ccc     &          'UNFORMATTED' )
ccc      END IF
ccc      WRITE (IODRAW)  0, ICODE, JTRACK, SNGL (ETRACK), SNGL (WTRACK)
ccc      WRITE (IODRAW)  SNGL (XSCO), SNGL (YSCO), SNGL (ZSCO), SNGL (RULL)
ccc*  +-------------------------------------------------------------------*
ccc*  |  Quenching is activated : calculate quenching factor
ccc*  |  and store quenched energy in DTQUEN(1, jbk)
ccc      IF ( LQEMGD ) THEN
ccc         RULLL = RULL
ccc         CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
ccc         WRITE (IODRAW) ( SNGL (DTQUEN(1, JBK)), JBK = 1, NQEMGD )
ccc      END IF
ccc*  |  end quenching
ccc*  +-------------------------------------------------------------------*
      RETURN
*
*======================================================================*
*                                                                      *
*     SOurce particle DRAWing:                                         *
*                                                                      *
*======================================================================*
*
      ENTRY SODRAW
ccc      IF ( .NOT. LFCOPE ) THEN
ccc         LFCOPE = .TRUE.
ccc         IF ( KOMPUT .EQ. 2 ) THEN
ccc            FILNAM = '/'//CFDRAW(1:8)//' DUMP A'
ccc         ELSE
ccc            FILNAM = CFDRAW
ccc         END IF
ccc         OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM =
ccc     &          'UNFORMATTED' )
ccc      END IF
ccc      WRITE (IODRAW) -NCASE, NPFLKA, NSTMAX, SNGL (TKESUM),
ccc     &                SNGL (WEIPRI)
ccc*  +-------------------------------------------------------------------*
ccc*  |  (Radioactive) isotope: it works only for 1 source particle on
ccc*  |  the stack for the time being
ccc      IF ( ILOFLK (NPFLKA) .GE. 100000 .AND. LRADDC (NPFLKA) ) THEN
ccc         IARES  = MOD ( ILOFLK (NPFLKA), 100000  )  / 100
ccc         IZRES  = MOD ( ILOFLK (NPFLKA), 10000000 ) / 100000
ccc         IISRES = ILOFLK (NPFLKA) / 10000000
ccc         IONID  = ILOFLK (NPFLKA)
ccc         WRITE (IODRAW) ( IONID,SNGL(-TKEFLK(I)),
ccc     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
ccc     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
ccc     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
ccc     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
ccc*  |
ccc*  +-------------------------------------------------------------------*
ccc*  |  Patch for heavy ions: it works only for 1 source particle on
ccc*  |  the stack for the time being
ccc      ELSE IF ( ABS (ILOFLK (NPFLKA)) .GE. 10000 ) THEN
ccc         IONID = ILOFLK (NPFLKA)
ccc         CALL DCDION ( IONID )
ccc         WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-IONID)),
ccc     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
ccc     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
ccc     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
ccc     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
ccc*  |
ccc*  +-------------------------------------------------------------------*
ccc*  |  Patch for heavy ions: ???
ccc      ELSE IF ( ILOFLK (NPFLKA) .LT. -6 ) THEN
ccc         WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-ILOFLK(NPFLKA))),
ccc     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
ccc     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
ccc     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
ccc     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
ccc*  |
ccc*  +-------------------------------------------------------------------*
ccc*  |
ccc      ELSE
ccc         WRITE (IODRAW) ( ILOFLK(I), SNGL (TKEFLK(I)+AM(ILOFLK(I))),
ccc     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
ccc     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
ccc     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
ccc     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
ccc      END IF
ccc*  |
ccc*  +-------------------------------------------------------------------*
ccc      
cccc         WRITE (42.0,*) ( ILOFLK(I), SNGL (TKEFLK(I)+AM(ILOFLK(I))),
cccc     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
cccc     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
cccc     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
cccc     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
ccc      write(42.0,*) (sngl(txflk(i)), sngl(tyflk(i)), sngl(tzflk(i)),
ccc     +               sngl(xflk(i)), sngl(yflk(i)), sngl(zflk(i)),
ccc     +               i=1,npflka)
      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 )

      nn=0
      hikin=0
      nocount=99999
      if (mreg.eq.8.0) then
          do 10 ip = 1, NP
c               WRITE( 41.0, *) "  ", ip,KPART(ip),TKI(ip),WEI(ip)
c               write( 41.0, *) "  ", cxr(ip), cyr(ip), czr(ip)
                if(kpart(ip).eq.8) nn=nn+1
 10       continue

          if (nn.gt.0 .and. jtrack.eq.8) then
              do ip=1,np
                  if(kpart(ip).eq.8 .and. tki(ip).gt.hikin) then
                      nocount=ip
                      hikin=tki(ip)
                  endif
              enddo
              do ip=1,np
                  if(ip.ne.nocount .and. kpart(ip).eq.8) then
*                     write(41.0,*) jtrack, xtrack(ntrack), ytrack(ntrack),
*         +                          ztrack(nctrak),tki(ip),cxr(ip),cyr(ip),
*         +                          czr(ip)
                      write(41.0,*)  jtrack, xsco, ysco, zsco,tki(ip),
     +                               NCASE, NPFLKA
c     +                              tki(ip),cxr(ip),cyr(ip),
c     +                              czr(ip)
                  endif
              enddo
          else if (nn.gt.0.and.jtrack.ne.8) then
               do ip=1,np
                   if (kpart(ip).eq.8) then
*                       write(41.0,*) jtrack,xtrack(ntrack),ytrack(ntrack),
*         +                          ztrack(nctrak),tki(ip),cxr(ip),cyr(ip),
*         +                          czr(ip)
                      write(41.0,*)  jtrack, xsco, ysco, zsco, tki(ip),
     +                               NCASE, NPFLKA
c     +                              tki(ip),cxr(ip),cyr(ip),
c     +                              czr(ip)
                   endif
                enddo

          endif
      endif


      if(mreg.eq.9.0.and.jtrack.eq.11.0) then
          write(44.0,*) mreg, ncase, etrack-am(jtrack), icode
      else if(mreg.eq.10.0.and.jtrack.eq.11.0) then
          write(44.0,*) mreg, ncase, etrack-am(jtrack), icode
      else if(mreg.eq.11.0.and.jtrack.eq.11.0) then
          write(44.0,*) mreg, ncase, etrack-am(jtrack), icode
      else if(mreg.eq.8.0.and.jtrack.eq.11.0) then
          write(44.0,*) mreg, ncase, etrack-am(jtrack), icode
      endif
             
          
cc* Here is where I write my output.
cc      WRITE( 41.0, * ) 'Routine USDRAW called'
cc      WRITE( 41.0, * )  NCASE, NPFLKA, JTRACK, ICODE, NP
cc      WRITE( 41.0, * )  ETRACK, PTRACK, WTRACK, MMTRCK, KTRACK
cc* Loop over Secondaries.
cc      do 10 ip = 1, NP
cc            WRITE( 41.0, *) ip,KPART(ip),TKI(ip),WEI(ip)
cc 10   continue
cc* Loop over Heavy Secondaries.
cc      WRITE( 41.0, * )  NPHEAV, ICRES, IBRES, ERES, EKRES
cc      do 20 ip = 1, NPHEAV
cc            WRITE( 41.0, *) ip,KHEAVY(ip),TKHEAV(ip)
cc            WRITE( 41.0, *) ip,ICHEAV(KHEAVY(ip)),IBHEAV(KHEAVY(ip))
cc 20   continue

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

