C----- at the beginning of GDRAW
      INCLUDE '(GENSTK)' ! daughter particles common block
      INCLUDE '(TRACKR)' ! mother particle common block
C-------------------------------

      INTEGER PIDKPIU,PIDKMENO,PIDE,PIDPOS
      PARAMETER(PIDKPIU = 15 )
      PARAMETER(PIDKMENO = 16 )
      PARAMETER(PIDE=3)
      PARAMETER(PIDPOS=4)

      REAL BRKE2
      PARAMETER(BRKE2=1.55D-5)


      ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
     

C check if the routine is called by a decay
      IF(ICODE.eq.102) THEN
C check if the decaying particle is a K+
        IF(JTRACK.eq.PIDKPIU.or.JTRACK.eq.PIDKMENO) THEN

          IF(FLRNDM(XXX).gt.BRKE2) GOTO 100

        MMOTH = AM(JTRACK) ! set the mass of the mother particle
        IF(PTRACK.ge.0.d0) THEN ! set the momentum of the mother particle
          PMMOTH=PTRACK
        ELSE
          PMMOTH = ETRACK**2-MMOTH**2 ! the same if PTRACK is not defined
          PMMOTH = sqrt(PMMOTH)
        ENDIF

        PPMOTH(1)=cxtrck*PMMOTH  ! set the x,y,z component of the momenta of the mother particle
        PPMOTH(2)=cytrck*PMMOTH
        PPMOTH(3)=cztrck*PMMOTH

        IF(JTRACK.eq.PIDKMENO) THEN
           MPART(1) = AM(PIDE)        ! set the type of the daughter particle
           Kpart(1) = PIDE            ! in GENSTK e+ for K+ decay e- for K- decay (overwriting current assignment)
        ELSE
           MPART(1) = AM(PIDPOS)      ! set also the mass of the first particle to use in the decay simulator
           Kpart(1) = PIDPOS
        ENDIF

        MPART(2) = 0D0 ! set the mass of the second particle (neutrino)

        CALL NBODYDEC(MMOTH,PPMOTH,2,MPART,PPART) ! Two body decay simulator
        NP=1                                      ! overwrite the number of decaying particles in GENSTK
        Cxr(1) = PPART(1,1)/PPART(5,1)            
        CXr(1) = PPART(2,1)/PPART(5,1)            ! overwrite the cosine director of the decaying particles
        Cxr(1) = PPART(3,1)/PPART(5,1)

        Tki(1) = PPART(4,1)-MPART(1)              ! kinetic energy
        Plr(1) = PPART(5,1)                       ! momentum
        Wei(1) = 1d0
        Agesec(1) = 0d0                           ! age respect to the decay time
 100    CONTINUE

      ENDIF
      ENDIF

