C...Program to turn gluinos into gluino-hadrons, assuming
C...the gluino is the (pseudo)stable LSP. Also possibility 
C...to let this gluino-hadron decay later on. Can be used
C...e.g. in split SUSY scenarios, where the gluino-hadron
C...might decay inside the detector.

C******************************************************************

C...All real arithmetic in double precision.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
C...Three Pythia functions return integers, so need declaring.
      INTEGER PYK,PYCHGE,PYCOMP
C...Parameter statement to help give large particle numbers
C...(left- and righthanded SUSY, technicolor, excited fermions,
C...extra dimensions).
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
     &KEXCIT=4000000,KDIMEN=5000000)

C...EXTERNAL statement links PYDATA on most machines.
      EXTERNAL PYDATA

C...Commonblocks.
C...The event record.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
C...Parameters.
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
C...Particle properties + some flavour parameters.
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
C...Decay information.
C...Note that dimensions below grew from 4000 to 8000 in Pythia 6.2!
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
C...Selection of hard scattering subprocesses.
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
C...Parameters. 
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
C...Process information.
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
C...Process names.
      COMMON/PYINT6/PROC(0:500)
      CHARACTER PROC*28
C...Supersymmetry parameters.
      COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
C...Local arrays.
      DIMENSION EPQSUM(6)

C...Number of events, CM energy.
      NEV=10000
      ECM=14000D0

C...Pick gluino mass.
      KFGL=KSUSY1+21
      PMGL=200D0

C...Pick processes: here gluino pair production.
      MSEL=0
      MSUB(243)=1
      MSUB(244)=1

C...Set necessary SUSY parameters (brute force).
      IMSS(1)=1
C...Gluino mass - of "constituent" mass kind.
      IMSS(3)=1
      RMSS(3)=PMGL
C...Pick other ino masses somewhat smaller.
      RMSS(1)=0.5D0*PMGL
      RMSS(2)=0.5D0*PMGL
      RMSS(4)=1D4
C...Stop_1 mass.
      IMSS(5)=1
      RMSS(12)=2D0*PMGL

C...Initialize.
      CALL PYINIT('CMS','p','p',ECM)
  
C...Show particle data.
C      CALL PYSTAT(2)
C      CALL PYLIST(12)

C...Switch off hadronization in normal PYEVNT call. Do not change!
      MSTP(111)=0 

C...Put gluino stable. Do not change!
      MDCY(PYCOMP(KFGL),1)=0

C...Histograms.
      CALL PYBOOK(1,'E-P-Q conservation check',100,0D0,0.01D0)
      CALL PYBOOK(3,'charged hadron multiplicity (excl. R hadrons)',
     &100,-1D0,399D0)
      CALL PYBOOK(5,'extra multiplicity from gluino decays',
     &100,-0.5D0,99.5D0)
      CALL PYBOOK(21,'gluino parton pT spectrum',100,0D0,500D0)
      CALL PYBOOK(22,'gluino hadron pT spectrum',100,0D0,500D0)
      CALL PYBOOK(23,'Gluino-hadron mass spectrum',100,PMGL,PMGL+2D0)
 
C...Event generation loop.
      DO 200 IEV=1,NEV
        if(mod(iev,100).eq.0) write(*,*) 'begin event no', iev

C...Generate event, but without worrying about small string systems.
        MSTJ(14)=-1
        CALL PYEVNT
        MSTJ(14)=1
        DO 140 J=1,6
          EPQSUM(J)=PYP(0,J)
 140    CONTINUE

C...Now perform treatment of gluino hadronization.
        CALL PYGLFR

C...If treatment failed, skip this event.
        IF(MSTU(24).NE.0) GOTO 200

C...List first few events.
c        IF(IEV.LE.2) CALL PYLIST(2)

C...Statistics: generic.
        NMERR=0
        EPQ=ABS(PYP(0,4)-EPQSUM(4))+ABS(PYP(0,1)-EPQSUM(1))+
     &  ABS(PYP(0,2)-EPQSUM(2))+ABS(PYP(0,3)-EPQSUM(3))+
     &  ABS(PYP(0,6)-EPQSUM(6))
        CALL PYFILL(1,EPQ,1D0)
        IF(EPQ.GT.1D-2) NMERR=NMERR+1
        NCH=0
        DO 160 I=1,N
          DM2=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
          IF(DM2.GT.1D-4) THEN
             WRITE(*,*) 'problem particle',I,K(I,2),DM2
            IF(DM2.GT.1D-2) NMERR=NMERR+1
          ENDIF
          IF(K(I,1).EQ.1.AND.PYCHGE(K(I,2)).NE.0) NCH=NCH+1
          PT=SQRT(P(I,1)**2+P(I,2)**2) 
 
C...Statistics: gluino.
          IF(K(I,1).EQ.12.AND.K(I,2).EQ.KFGL) THEN
            CALL PYFILL(21,PT,1D0)
          ENDIF
          IF((K(I,1).EQ.6.OR.K(I,1).EQ.7).AND.K(I,2).EQ.KFGL) THEN
            CALL PYFILL(22,PT,1D0)
            CALL PYFILL(23,P(I,5),1D0)
          ENDIF
 160    CONTINUE
        CALL PYFILL(3,DBLE(NCH),1D0)
        IF(NMERR.GT.0) CALL PYLIST(2)

C...Now force gluinos to decay. Look in PYGLDC for options to make 
C...only one of them decay inside the detector.
        NSAV=N
        CALL PYGLDC

C...Charged multiplicity in gluino decays.
        NCH=0
        DO 180 I=NSAV+1,N
          IF(K(I,1).EQ.1.AND.PYCHGE(K(I,2)).NE.0) NCH=NCH+1
 180    CONTINUE
        CALL PYFILL(5,DBLE(NCH),1D0)

C...List first few events.
        IF(IEV.LE.2) CALL PYLIST(2)

C...End of event generation loop.
 200  CONTINUE

C...Cross section - not relevant in this case. Histograms.
      CALL PYSTAT(1)      
      CALL PYHIST

      END
 
C*********************************************************************
 
C...PYGLFR
C...Fragments the string near to a gluino, to form a gluino-hadron, 
C...either by producing a new g-g pair or two new q-qbar ones.
 
      SUBROUTINE PYGLFR
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      INTEGER PYK,PYCHGE,PYCOMP
C...Parameter statement to help give large particle numbers
C...(left- and righthanded SUSY, excited fermions).
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
C...Commonblocks.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
C...Note that dimensions below grew from 4000 to 8000 in Pythia 6.2!
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
      SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYINT1/,
     &/PYINT2/
C...Local array.
      DIMENSION PSUM(5),KFSAV(2),PMSAV(2),PSAV(2,5) 

C...Free parameter: relative probability for gluino-gluon-ball.
C...(But occasional low-mass string will never become it anyway.)
      PROBGG=0.1D0
 
C...Free parameter: gluon constituent mass.
      PMGLU=0.7D0
 
C...Free parameter: max kinetic energy in gluino-hadron.
      PMKIN=0.5D0

C...Switch off popcorn baryon production. (Not imperative, but more
C...failed events when popcorn is allowed.)
      MSTJ12=MSTJ(12)
      MSTJ(12)=1

C...Convenient shorthand.
      KFGL=KSUSY1+21

C...Loopback point for serious problems, with new try.
      LOOP=0
      CALL PYEDIT(21)
      CHGSAV=PYP(0,6)
   90 LOOP=LOOP+1
      IF(LOOP.GT.1) CALL PYEDIT(22)

C...Take copy of string system(s), leaving extra free slot after gluino.
C...(Eventually to be overwritten by one q and one qbar string break.)
      NOLD=N
      NGLUI=0
      DO 120 I=1,NOLD
        ICOPY=0
        IF(K(I,1).EQ.2) ICOPY=1
        IF(K(I,1).EQ.1.AND.I.GE.2) THEN
          IF(K(I-1,1).EQ.12) ICOPY=1
        ENDIF
        IF(ICOPY.EQ.1) THEN  
          N=N+1
          DO 100 J=1,5
            K(N,J)=K(I,J)
            P(N,J)=P(I,J)
            V(N,J)=V(I,J)
  100     CONTINUE
          K(I,1)=K(I,1)+10
          K(I,4)=N
          K(I,5)=N
          K(N,3)=I
          IF(K(I,2).EQ.KFGL) THEN
            NGLUI=NGLUI+1  
            N=N+1
            DO 110 J=1,5
              K(N,J)=K(N-1,J)
              P(N,J)=0D0
              V(N,J)=V(I,J)
  110       CONTINUE
            K(I,5)=N
            K(N,2)=21
          ENDIF
        ENDIF
  120 CONTINUE

C...Loop over (up to) two gluinos per event.
      DO 300 IGLUI=1,NGLUI

C...Identify position of gluino (randomize order of treatment).
        IGL=0
        NGL=0
        DO 130 I=1,N
          IF(K(I,1).EQ.2.AND.K(I,2).EQ.KFGL) THEN
            NGL=NGL+1
            IF(IGLUI.EQ.NGLUI) THEN
              IGL=I
            ELSEIF(NGL.EQ.1) THEN
              IF(PYR(0).LT.0.5D0) IGL=I
            ELSEIF(IGL.EQ.0) THEN
              IGL=I
            ENDIF
          ENDIF
  130   CONTINUE

C...Identify range of partons on string the gluino belongs to. 
        IMIN=IGL
  140   IMIN=IMIN-1
        IF(K(IMIN-1,1).EQ.2) GOTO 140
        IMAX=IGL
  150   IMAX=IMAX+1
        IF(K(IMAX,1).EQ.2) GOTO 150
 
C...Find mass of this gluino-string. 
        DO 170 J=1,5
          PSUM(J)=0D0
          DO 160 I=IMIN,IMAX
            PSUM(J)=PSUM(J)+P(I,J)
  160     CONTINUE
  170   CONTINUE
        PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-
     &  PSUM(3)**2))
 
C...If low-mass, then consider gluino-hadron already formed.
        IF(PSUM(5).LE.P(IGL,5)+P(IMIN,5)+P(IMAX,5)+PMKIN) THEN
          DO 180 I=IMIN,IMAX
            K(I,1)=15+IGLUI
  180     CONTINUE
          GOTO 300
        ENDIF    

C...Else break string by production of new gg or two new qqbar pairs.
C...(Also diquarks allowed, but not popcorn, and not two adjacent.)
        IF(PYR(0).LT.PROBGG) THEN
C...Let a gluon occupy two slots, to make administration work the same
C...way as for the qqbar case.
          KFSAV(1)=21
          KFSAV(2)=21
          PMSAV(1)=0.5D0*PMGLU  
          PMSAV(2)=0.5D0*PMGLU  
        ELSE
  185     CALL PYDCYK(K(IMIN,2),0,KFSAV(1),KFDUM)
          CALL PYDCYK(K(IMAX,2),0,KFSAV(2),KFDUM)
          IF(IABS(KFSAV(1)).GT.10.AND.IABS(KFSAV(2)).GT.10) GOTO 185
          IF(IABS(KFSAV(1)).GT.10.AND.IABS(K(IGL-1,2)).GT.10) GOTO 185
          IF(IABS(KFSAV(2)).GT.10.AND.IABS(K(IGL+2,2)).GT.10) GOTO 185
          KFSAV(1)=ISIGN(MOD(IABS(KFSAV(1)),10000),KFSAV(1))
          KFSAV(2)=ISIGN(MOD(IABS(KFSAV(2)),10000),KFSAV(2))
          MSTJ(93)=1 
          PMSAV(1)=PYMASS(KFSAV(1))         
          MSTJ(93)=1 
          PMSAV(2)=PYMASS(KFSAV(2))
        ENDIF         

C...Mass of gluino-hadron.
        PMGSAV=P(IGL,5)
        PMGB=P(IGL,5)+PMSAV(1)+PMSAV(2)

C...Pick at random order in which both sides of gluino string break.
        ISIDE=INT(1.5D0+PYR(0))
        DO 220 ISDE=1,2
          IF(ISDE.EQ.1) IS=ISIDE
          IF(ISDE.EQ.2) IS=3-ISIDE

C...Pick momentum sharing according to fragmentation function as if bottom.
          PMBSAV=PARF(105)
          PARF(105)=PMGSAV
          CALL PYZDIS(5,0,PMGB**2,ZGL)
          PARF(105)=PMBSAV 
          ZGL=MAX(0.9D0,MIN(0.9999D0,ZGL)) 
          DO 190 J=1,5
            PSAV(IS,J)=(1D0-ZGL)*P(IGL,J)
            P(IGL,J)=ZGL*P(IGL,J)
  190    CONTINUE

C...Target gluino-hadron mass for this stage of momentum reshuffling.
          PMOLD=P(IGL,5)
          IF(ISDE.EQ.1) PMNEW=PMGSAV+PMSAV(IS)
          IF(ISDE.EQ.2) PMNEW=PMGB 

C...Recoiling parton from which to shuffle momentum. System momentum.
          IF(IS.EQ.1) IREC=IGL-1
          IF(IS.EQ.2) IREC=IGL+2
  200     DO 210 J=1,4
            PSUM(J)=P(IGL,J)+P(IREC,J)
  210     CONTINUE           

C...Boost to rest frame of system, and align gluino along +z axis.
          CALL PYROBO(IGL,IGL,0D0,0D0,-PSUM(1)/PSUM(4),
     &    -PSUM(2)/PSUM(4),-PSUM(3)/PSUM(4))
          CALL PYROBO(IREC,IREC,0D0,0D0,-PSUM(1)/PSUM(4),
     &    -PSUM(2)/PSUM(4),-PSUM(3)/PSUM(4))
          PHI=PYANGL(P(IGL,1),P(IGL,2))
          CALL PYROBO(IGL,IGL,0D0,-PHI,0D0,0D0,0D0)
          CALL PYROBO(IREC,IREC,0D0,-PHI,0D0,0D0,0D0)
          THETA=PYANGL(P(IGL,3),P(IGL,1)) 
          CALL PYROBO(IGL,IGL,-THETA,0D0,0D0,0D0,0D0)
          CALL PYROBO(IREC,IREC,-THETA,0D0,0D0,0D0,0D0)

C...Calculate new kinematics in this frame, for desired gluino mass.
          ETOT=P(IGL,4)+P(IREC,4)
          IF(ETOT.GT.PMNEW+P(IREC,5)) THEN
            IFAIL=0
            PZNEW=0.5D0*SQRT(MAX(0D0,(ETOT**2-PMNEW**2-P(IREC,5)**2)**2-
     &      4D0*PMNEW**2*P(IREC,5)**2))/ETOT
            P(IGL,3)=PZNEW
            P(IGL,4)=SQRT(PZNEW**2+PMNEW**2)
            P(IGL,5)=PMNEW
            P(IREC,3)=-PZNEW
            P(IREC,4)=SQRT(PZNEW**2+P(IREC,5)**2)

C...If not enough momentum, take what can be taken.
          ELSE
            IFAIL=1
            PMOLD=ETOT-P(IREC,5)
            P(IGL,3)=0D0
            P(IGL,4)=PMOLD
            P(IGL,5)=PMOLD
            P(IREC,3)=0D0
            P(IREC,4)=P(IREC,5)
          ENDIF

C...Bost back to lab frame.
          CALL PYROBO(IGL,IGL,THETA,PHI,PSUM(1)/PSUM(4),
     &    PSUM(2)/PSUM(4),PSUM(3)/PSUM(4))
          CALL PYROBO(IREC,IREC,THETA,PHI,PSUM(1)/PSUM(4),
     &    PSUM(2)/PSUM(4),PSUM(3)/PSUM(4))

C...Loop back when not enough momentum could be shuffled.
C...(As long as there is something left on either side.)
          IF(IFAIL.EQ.1) THEN
  215       IF(IS.EQ.1.AND.IREC.GT.IMIN) THEN
              IREC=IREC-1
              GOTO 200
            ELSEIF(IS.EQ.2.AND.IREC.LT.IMAX) THEN
              IREC=IREC+1
              GOTO 200
            ELSEIF(ISDE.EQ.2.AND.IS.EQ.3-ISIDE) THEN
              IS=ISIDE
              IREC=IRECSV
              GOTO 215
            ENDIF
          ENDIF

C...End loop over fragmentation of two sides around gluino.
         IRECSV=IREC
  220   CONTINUE

C...New slots at end of record for gluino + new quarks/gluons.
        DO 230 J=1,5
          K(N+1,J)=0
          P(N+1,J)=P(IGL,J)
          V(N+1,J)=V(IGL,J)
          K(N+2,J)=0
          P(N+2,J)=0D0
          V(N+2,J)=V(IGL,J)
          K(N+3,J)=0
          P(N+3,J)=0D0
          V(N+3,J)=V(IGL,J)
  230   CONTINUE
 
C...Status and code of these slots.
        K(N+1,1)=15+IGLUI
        K(N+2,1)=15+IGLUI
        K(N+3,1)=15+IGLUI
        K(N+1,2)=K(IGL,2)
        IF(KFSAV(1).EQ.21) THEN
          K(N+2,2)=KFSAV(1)
          K(N+3,2)=KFSAV(2)
        ELSE
          K(N+2,2)=-KFSAV(1)
          K(N+3,2)=-KFSAV(2)
        ENDIF
        K(N+1,3)=K(IGL,3)
        K(N+2,3)=K(IGL,3)
        K(N+3,3)=K(IGL,3)
        N=N+3
        
C...Code and momentum of two new string endpoints.
        K(IGL,2)=KFSAV(1)
        K(IGL+1,2)=KFSAV(2)
        IF(KFSAV(1).NE.21) K(IGL,1)=1
        DO 240 J=1,5
          P(IGL,J)=PSAV(1,J)
          P(IGL+1,J)=PSAV(2,J)
  240   CONTINUE
 
C...End of loop over two gluinos.
  300 CONTINUE

C...Cleanup: remove zero-energy gluons.
      NNOW=N
      N=NOLD
      DO 330 I=NOLD+1,NNOW
        IF(K(I,2).EQ.21.AND.P(I,4).LT.1D-10) THEN
        ELSEIF(I.EQ.N+1) THEN
          N=N+1
        ELSE
          N=N+1
          DO 320 J=1,5
            K(N,J)=K(I,J)
            P(N,J)=P(I,J)
            V(N,J)=V(I,J)
  320     CONTINUE
        ENDIF
  330 CONTINUE

C...Finish off with standard hadronization.
      CALL PYEXEC

C...Restore code of gluino-hadrons to non-fragmented numbers.
C...(Status code K(I,1) = 6 and 7 can now be used to identify 
C...constituents of the two gluino-hadrons.)
C...Check that exactly two extrapartons accompanying each gluino.
      N6=0
      N7=0
      DO 340 I=1,N
        IF(K(I,1).EQ.16.OR.K(I,1).EQ.17) K(I,1)=K(I,1)-10
        IF(K(I,1).EQ.6) N6=N6+1
        IF(K(I,1).EQ.7) N7=N7+1
  340 CONTINUE
      IF(N6.NE.3.OR.N7.NE.3) MSTU(24)=1

C...Extracheck charge.
      CHGNEW=PYP(0,6)
      IF(ABS(CHGNEW-CHGSAV).GT.0.1D0) MSTU(24)=1

C...In case of trouble, make up to five attempts.
      IF(MSTU(24).NE.0.AND.LOOP.LT.5) THEN
        MSTU(23)=MSTU(23)-1
        GOTO 90
      ELSEIF(MSTU(24).NE.0) THEN
        WRITE(*,*) 'Failed after repeated attempts!'
      ELSEIF(MSTU(24).EQ.0.AND.LOOP.GT.1) THEN
      ENDIF

C...Finished. Restore baryon production model.
      MSTJ(12)=MSTJ12
      
      RETURN
      END
 
C*********************************************************************
 
C...PYGLDC
C...Decays the gluino inside a gluino-hadron.
C...Optionally you may pick to let only one or none decay, see below.
 
      SUBROUTINE PYGLDC

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      INTEGER PYK,PYCHGE,PYCOMP
C...Parameter statement to help give large particle numbers
C...(left- and righthanded SUSY, excited fermions).
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
C...Commonblocks.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
C...Note that dimensions below grew from 4000 to 8000 in Pythia 6.2!
      COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
      SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYINT1/,
     &/PYINT2/

C...Set light quark mass, which defines kinematics.
      PMLQ = 0.5D0

C...Loop through to find gluinos in systems (at most two).
      KFGL=KSUSY1+21
      IGL1=0
      IGL2=0
      DO 100 I=1,N
        IF(K(I,1).EQ.6.AND.K(I,2).EQ.KFGL) IGL1=I
        IF(K(I,1).EQ.7.AND.K(I,2).EQ.KFGL) IGL2=I
 100  CONTINUE  

C...The two gluino(-hadrons) are now stored in lines IGL1 AND IGL2
C...(with either zero if only one exists). You may here use the
C...momentum information to decide whether either or both should
C...decay inside the detector, and set that information.
C...Default behaviour is that both decay, =1, if not set =0.
C...This would also be an excellent place if you want to store
C...information on a displaced event vertex in the V array.
      MGL1DC=1
      MGL2DC=1

C...Loop over the two gluino-hadrons. Skip for intended stable.
      DO 150 IGL=1,2
        IF(IGL.EQ.1.AND.MGL1DC.LE.0) GOTO 150
        IF(IGL.EQ.2.AND.MGL2DC.LE.0) GOTO 150
        I=IGL1
        IF(IGL.EQ.2) I=IGL2

C...Check that system is of expected kind, else warn.
        IF(K(I+1,1).NE.K(I,1).OR.K(I+2,1).NE.K(I,1)) THEN
          WRITE(*,*) 'PROBLEM WITH GLUINO IN LINE =',I
          CALL PYLIST(1)
        ENDIF

C...Share gluino"hadron" momentum with two light quarks.
        FAC=PMLQ/P(I,5)
        DO 120 J=1,5
          P(I+1,J)=FAC*P(I,J)
          P(I+2,J)=FAC*P(I,J)
          P(I,J)=(1D0-2D0*FAC)*P(I,J)
 120    CONTINUE

C...Set up colour flow.
        K(I,1)=3
        K(I+1,1)=3
        K(I+2,1)=3
        IF((K(I+1,2).GE.1.AND.K(I+1,2).LT.10).OR.K(I+1,2).LT.-10) THEN
          K(I,4)=MSTU(5)*(I+2)
          K(I,5)=MSTU(5)*(I+1)
          K(I+1,4)=MSTU(5)*I
          K(I+2,5)=MSTU(5)*I
        ELSE
          K(I,5)=MSTU(5)*(I+2)
          K(I,4)=MSTU(5)*(I+1)
          K(I+1,5)=MSTU(5)*I
          K(I+2,4)=MSTU(5)*I
        ENDIF

C..Put gluino unstable.
        MDCY(PYCOMP(KFGL),1)=1

C...Let gluino decay now.
        CALL PYRESD(I)

C...Put gluino stable.
        MDCY(PYCOMP(KFGL),1)=0

 150  CONTINUE

C...Temporarily put undecayed gluinos as decayed, 
C...to avoid PYEXEC from trying (and failing) to decay them.
      IF(MGL1DC.LE.0) THEN 
        K(IGL1,1)=K(IGL1,1)+10
        K(IGL1+1,1)=K(IGL1+1,1)+10
        K(IGL1+2,1)=K(IGL1+2,1)+10
      ENDIF
      IF(MGL2DC.LE.0) THEN 
        K(IGL2,1)=K(IGL2,1)+10
        K(IGL2+1,1)=K(IGL2+1,1)+10
        K(IGL2+2,1)=K(IGL2+2,1)+10
      ENDIF

C...Let the gluino decay products fragment, etc.
      CALL PYEXEC

C...Undo temporary change of status code above.
      IF(MGL1DC.LE.0) THEN 
        K(IGL1,1)=K(IGL1,1)-10
        K(IGL1+1,1)=K(IGL1+1,1)-10
        K(IGL1+2,1)=K(IGL1+2,1)-10
      ENDIF
      IF(MGL2DC.LE.0) THEN 
        K(IGL2,1)=K(IGL2,1)-10
        K(IGL2+1,1)=K(IGL2+1,1)-10
        K(IGL2+2,1)=K(IGL2+2,1)-10
      ENDIF

      RETURN
      END