*$ 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. /
      COMMON/NUMBERIN/NUMSIZE
*
*----------------------------------------------------------------------*
*                                                                      *
*     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                                      *
*                                                                      *
*----------------------------------------------------------------------*

*  |  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 )

      IF (.NOT. LFCOPE) THEN
        LFCOPE = .TRUE.
        OPEN(UNIT=60, FILE = 'Ang.60',FORM = 'FORMATTED',
     &        STATUS = 'UNKNOWN')

      END IF



      IF( JTRACK.EQ.10) THEN ! Select muon+
        IF(ETRACK.GT.AM(JTRACK)) THEN ! Muon has survived
          TPPTRACK = ETRACK - AM(JTRACK)

*           IF((TPPTRACK .GE. 3.9E-3) .AND. (TPPTRACK .GE. 4.1E-3)) THEN
* Calculate the x' and y' in phase diagram.
* NOTE: Now the beam direction is x axis for 90 degree detecter
*       y and z is the horizotal and vertical directions respectively
*

************************ 0 degree ***************************
      IF(MREG.EQ.6 .AND. NEWREG.EQ.7) THEN ! 0 degree
* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      cosn = CZTRCK
         IF(cosn .GE. coszTH) THEN

      IANG = 0
      xn = XSCO*10.
      yn = YSCO*10.

      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK

      XPP = 1000. * PPx / PPz
      YPP = 1000. * PPy / PPz

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK
         END IF

************************ 15 degree ***************************
      ELSE IF(MREG.EQ.8 .AND. NEWREG.EQ.9) THEN ! 15 degree
      IANG = 15

* Rotate an angle (AngR) with repect to y axis
      AngR = 15.*3.1415926/180.


* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK
      CALL ROTSYS(PPx,PPy,PPz,ANGR,Pxn,Pyn,Pzn)
* direction cosine in new coordinate system
      cosn = Pzn / PTRACK

         IF(cosn .GE. coszTH) THEN

      CALL ROTSYS(XSCO,YSCO,ZSCO,ANGR,xn,yn,zn)

      xn = xn*10.
      yn = yn*10.


      XPP = 1000. * Pxn / Pzn
      YPP = 1000. * pyn / Pzn

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK

          END IF

************************ 30 degree ***************************
      ELSE IF(MREG.EQ.10 .AND. NEWREG.EQ.11) THEN ! 30 degree
      IANG = 30

* Rotate an angle (AngR) with repect to y axis
      AngR = 30.*3.1415926/180.

* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK
      CALL ROTSYS(PPx,PPy,PPz,ANGR,Pxn,Pyn,Pzn)
* direction cosine in new coordinate system
      cosn = Pzn / PTRACK

         IF(cosn .GE. coszTH) THEN

      CALL ROTSYS(XSCO,YSCO,ZSCO,ANGR,xn,yn,zn)
      xn = xn*10.
      yn = yn*10.


      XPP = 1000. * Pxn / Pzn
      YPP = 1000. * pyn / Pzn

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK

          END IF


************************ 45 degree ***************************
      ELSE IF(MREG.EQ.12 .AND. NEWREG.EQ.13) THEN ! 45 degree
      IANG = 45

* Rotate an angle (AngR) with repect to y axis
      AngR = 45.*3.1415926/180.

* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK
      CALL ROTSYS(PPx,PPy,PPz,ANGR,Pxn,Pyn,Pzn)
* direction cosine in new coordinate system
      cosn = Pzn / PTRACK

         IF(cosn .GE. coszTH) THEN

      CALL ROTSYS(XSCO,YSCO,ZSCO,ANGR,xn,yn,zn)
      xn = xn*10.
      yn = yn*10.


      XPP = 1000. * Pxn / Pzn
      YPP = 1000. * pyn / Pzn

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK

         END IF

************************ 60 degree ***************************
      ELSE IF(MREG.EQ.14 .AND. NEWREG.EQ.15) THEN ! 60 degree
      IANG = 60

* Rotate an angle (AngR) with repect to y axis
      AngR = 60.*3.1415926/180.

* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK
      CALL ROTSYS(PPx,PPy,PPz,ANGR,Pxn,Pyn,Pzn)
* direction cosine in new coordinate system
      cosn = Pzn / PTRACK

         IF(cosn .GE. coszTH) THEN

      CALL ROTSYS(XSCO,YSCO,ZSCO,ANGR,xn,yn,zn)
      xn = xn*10.
      yn = yn*10.


      XPP = 1000. * Pxn / Pzn
      YPP = 1000. * pyn / Pzn

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK

          END IF

************************ 75 degree ***************************
      ELSE IF(MREG.EQ.16 .AND. NEWREG.EQ.17) THEN ! 75 degree
      IANG = 75

* Rotate an angle (AngR) with repect to y axis
      AngR = 75.*3.1415926/180.

* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK
      CALL ROTSYS(PPx,PPy,PPz,ANGR,Pxn,Pyn,Pzn)
* direction cosine in new coordinate system
      cosn = Pzn / PTRACK

         IF(cosn .GE. coszTH) THEN
      CALL ROTSYS(XSCO,YSCO,ZSCO,ANGR,xn,yn,zn)
      xn = xn*10.
      yn = yn*10.


      XPP = 1000. * Pxn / Pzn
      YPP = 1000. * pyn / Pzn

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK

          END IF

************************ 90 degree ***************************
      ELSE IF(MREG.EQ.18 .AND. NEWREG.EQ.19) THEN ! 90 degree
      IANG = 90

* Rotate an angle (AngR) with repect to y axis
      AngR = 90.*3.1415926/180.

* Solid angle (unit: sr) and x direction cosine
      srad = 0.8D0
      coszTH = 1.0D0 - srad/(2. * 3.1415926535897932384626433832795)
      PPx = CXTRCK * PTRACK
      PPy = CYTRCK * PTRACK
      PPz = CZTRCK * PTRACK
      CALL ROTSYS(PPx,PPy,PPz,ANGR,Pxn,Pyn,Pzn)
* direction cosine in new coordinate system
      cosn = Pzn / PTRACK

         IF(cosn .GE. coszTH) THEN

      CALL ROTSYS(XSCO,YSCO,ZSCO,ANGR,xn,yn,zn)
      xn = xn*10.
      yn = yn*10.


      XPP = 1000. * Pxn / Pzn
      YPP = 1000. * pyn / Pzn

      WRITE(60,200)IANG, xn, XPP, yn, YPP, TPPTRACK, ATRACK
           END IF

      END IF




        ENDIF
      ENDIF

  200 format(I4, 6(1x, E12.5))



      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 )


*  |  end quenching
*  +-------------------------------------------------------------------*
      RETURN
*
*======================================================================*
*                                                                      *
*     SOurce particle DRAWing:                                         *
*                                                                      *
*======================================================================*
*
      ENTRY SODRAW

*  |
*  +-------------------------------------------------------------------*
      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 )

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





      SUBROUTINE ROTSYS(XO,YO,ZO,ANG,XN,YN,ZN)
*Rotating an angle around the y axis
	DOUBLE PRECISION XO,YO,ZO,ANG,XN,YN,ZN

      XN = XO * COS(Ang) - ZO * SIN(Ang)
	YN = YO
	ZN = ZO * COS(Ang) + XO * SIN(Ang)

	RETURN
	END


