*$ 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. /
*
      integer lablpdg(-9:40)
      data lablpdg /    0,     0,     0,    47,    49    
     &             ,   46,    45,     0,     1,    48
     &             ,   14,    15,     3,     2,     4 
     &             ,    4,     1,    13,    25,     5
     &             ,    6,    10,     8,     9,    11
     &             ,   12,    18,    26,    16,    21
     &             ,   19,    20,     7,    10,    10
     &             ,    0,     4,     4,     0,     0
     &             ,   27,    28,    29,    22,    30
     &             ,   23,    31,    24,    32,     0/
*variables for Edep check, Npart_CV check, and seeds recording for Event
      real ecv(36),ecc00(2),ecc02(8),emb(32),efb(16),ecc03(6),
     &     theta,
     &     xseed(200),yseed(200),zseed(200),pxseed(200),
     &     pyseed(200),pzseed(200),pidseed(200),tseed(200)
      integer Nsecondaries,Ncsi,cvid,cc00id,cc02id,mbid,fbid,cc03id
      
*
*----------------------------------------------------------------------*
*                                                                      *
*     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
c      WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
c     &               SNGL (WTRACK)
c      WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
c     &                 SNGL (ZTRACK (I)), I = 0, NTRACK ),
c     &               ( SNGL (DTRACK (I)), I = 1, MTRACK ),
c     &                 SNGL (CTRACK)
*  +-------------------------------------------------------------------*
*  |  Quenching is activated
      IF ( LQEMGD ) THEN
         RULLL = ZERZER
         CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
c         WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)), I = 1, MTRACK ),
c     &                      JBK = 1, NQEMGD )
      END IF
*  |  End of quenching
*  +-------------------------------------------------------------------*
*add Dtrack
      if (MTRACK.gt.0) then
         do I=1,MTRACK
            theta=atan2(SNGL(YTRACk(I)),SNGL(XTRACK(I)))+3.14159265359
            if (MREG.ge.681.and.MREG.le.749) then
               if ((abs(SNGL(YTRACk(I))).gt.6.0).or.
     +              (abs(SNGL(XTRACk(I))).gt.6.0)) then
                  cvid=int(theta/(2.0*3.14159265359/32.0))+1
               else
                  cvid=int(theta/(2.0*3.14159265359/4.0))+33
               endif
               ecv(cvid)=ecv(cvid)+DTRACK(I)
            else if (MREG.ge.198.and.MREG.le.208) then
               cc00id=int(theta/(2.0*3.14159265359/2.0))+1
               ecc00(cc00id)=ecc00(cc00id)+DTRACK(I)
            else if (MREG.ge.247.and.MREG.le.583) then
               if (mod(MREG-247,8).eq.0) then
                  cc02id=int(theta/(2.0*3.14159265359/8.0))+1
                  ecc02(cc02id)=ecc02(cc02id)+DTRACK(I)
               endif
            else if (MREG.ge.3.and.MREG.le.92) then
               if ((mod(MREG-3,2).eq.0).or.(MREG.eq.92)) then
                  mbid=int(theta/(2.0*3.14159265359/32.0))+1
                  emb(mbid)=emb(mbid)+DTRACK(I)
               endif
            else if (MREG.ge.94.and.MREG.le.194) then
               if (mod(MREG-94,2).eq.0) then
                  fbid=int(theta/(2.0*3.14159265359/16.0))+1
                  efb(fbid)=efb(fbid)+DTRACK(I)
               endif
            else if (MREG.ge.593.and.MREG.le.634) then
                  cc03id=mod((MREG-593),6)+1
                  ecc03(cc03id)=ecc03(cc03id)+DTRACK(I)
            endif
         enddo
      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 )
      if ( (mreg.eq.93).and.((newreg.eq.589)
     +     .or.(newreg.eq.590)) ) then
c      if ( (mreg.eq.196).and.(newreg.eq.93) ) then
*      if ( (mreg.eq.2).and.(newreg.eq.2) ) then
*        if (lablpdg(jtrack).eq.13) go to 100
*        write(49,*) 'BXDRAW:',ncase,icode,mreg,newreg,xsco,ysco,zsco,atrack
*        write(49,*) npsoev,(nevent(i),i=1,5)
*        write(49,*) ntrack,mtrack,jtrack,etrack,cxtrck,cytrck,cztrck,wtrack
*        if (cztrck.lt.0.0) go to 100
        asupx=ptrack*cxtrck
        asupy=ptrack*cytrck
        asupz=ptrack*cztrck
        Nsecondaries=Nsecondaries+1
        pidseed(Nsecondaries)=jtrack
        tseed(Nsecondaries)=atrack
        xseed(Nsecondaries)=xsco
        yseed(Nsecondaries)=ysco
        zseed(Nsecondaries)=zsco
        pxseed(Nsecondaries)=asupx
        pyseed(Nsecondaries)=asupy
        pzseed(Nsecondaries)=asupz
*        write(49,10) ncase,lablpdg(jtrack),atrack,xsco,ysco,zsco
*     &    ,asupx,asupy,asupz,wtrack
      endif
*   10 format(i8,i4,1pe13.5,3(0pf12.6),3(1pe14.5),0pf5.1)
  100 continue
      RETURN
*
*======================================================================*
*                                                                      *
*     Event End DRAWing:                                               *
*                                                                      *
*======================================================================*
*                                                                      *
      ENTRY EEDRAW ( ICODE )
      if (Nsecondaries.gt.1) then
c      if (ecc02.gt.0.000001) then
c         do I=1,Nsecondaries
c            write(49,10) ncase,lablpdg(pidseed(I)),tseed(I),xseed(I),
c     &          yseed(I),zseed(I),pxseed(I),pyseed(I),pzseed(I),wtrack
c         enddo
         write(49,10) ecc02(1),ecc00(1),ecv(1)
      endif
c   10 format(i8,i4,1pe13.5,3(0pf12.6),3(1pe14.5),0pf5.1)
   10 format(3(1pe14.5))
      do I=1,36
         ecv(I)=0.0
         if (I.le.32) then
            emb(I)=0.0
         endif
         if (I.le.16) then
            efb(I)=0.0
         endif
         if (I.le.8) then
            ecc02(I)=0.0
         endif
         if (I.le.6) then
            ecc03(I)=0.0
         endif
         if (I.le.2) then
            ecc00(I)=0.0
         endif
      enddo
      Nsecondaries=0
      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 ( .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
c      WRITE (IODRAW)  0, ICODE, JTRACK, SNGL (ETRACK), SNGL (WTRACK)
c      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 )
c         WRITE (IODRAW) ( SNGL (DTQUEN(1, JBK)), JBK = 1, NQEMGD )
      END IF
*  |  end quenching
*  +-------------------------------------------------------------------*
*add RULL
      theta=atan2(SNGL(YSCO),SNGL(XSCO))+3.14159265359
      if (MREG.ge.681.and.MREG.le.749) then
         if ((abs(SNGL(YSCO)).gt.6.0).or.
     +        (abs(SNGL(XSCO)).gt.6.0)) then
            cvid=int(theta/(2.0*3.14159265359/32.0))+1
         else
            cvid=int(theta/(2.0*3.14159265359/4.0))+33
         endif
         ecv(cvid)=ecv(cvid)+RULL
      else if (MREG.ge.198.and.MREG.le.208) then
         cc00id=int(theta/(2.0*3.14159265359/2.0))+1
         ecc00(cc00id)=ecc00(cc00id)+RULL
      else if (MREG.ge.247.and.MREG.le.583) then
         if (mod(MREG-247,8).eq.0) then
            cc02id=int(theta/(2.0*3.14159265359/8.0))+1
            ecc02(cc02id)=ecc02(cc02id)+RULL
         endif
      else if (MREG.ge.3.and.MREG.le.92) then
         if ((mod(MREG-3,2).eq.0).or.(MREG.eq.92)) then
            mbid=int(theta/(2.0*3.14159265359/32.0))+1
            emb(mbid)=emb(mbid)+RULL
         endif
      else if (MREG.ge.94.and.MREG.le.194) then
         if (mod(MREG-94,2).eq.0) then
            fbid=int(theta/(2.0*3.14159265359/16.0))+1
            efb(fbid)=efb(fbid)+RULL
         endif
      else if (MREG.ge.593.and.MREG.le.634) then
         cc03id=mod((MREG-593),6)+1
         ecc03(cc03id)=ecc03(cc03id)+RULL
      endif
      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)
c         WRITE (IODRAW) ( IONID,SNGL(-TKEFLK(I)),
c     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
c     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
c     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
c     &                    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 )
c         WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-IONID)),
c     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
c     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
c     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
c     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
*  |
*  +-------------------------------------------------------------------*
*  |  Patch for heavy ions: ???
      ELSE IF ( ILOFLK (NPFLKA) .LT. -6 ) THEN
c         WRITE (IODRAW) ( IONID,SNGL(TKEFLK(I)+AMNHEA(-ILOFLK(NPFLKA))),
c     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
c     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
c     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
c     &                    SNGL (TZFLK(I)), I = 1, NPFLKA )
*  |
*  +-------------------------------------------------------------------*
*  |
      ELSE
c         WRITE (IODRAW) ( ILOFLK(I), SNGL (TKEFLK(I)+AM(ILOFLK(I))),
c     &                    SNGL (WTFLK(I)), SNGL (XFLK (I)),
c     &                    SNGL (YFLK (I)), SNGL (ZFLK (I)),
c     &                    SNGL (TXFLK(I)), SNGL (TYFLK(I)),
c     &                    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

