C...Program to turn gluinos into gluino-hadrons, assuming
C...the gluino is the stable LSP. Furthermore, optionally, 
C...stop is assumed long-lived, so that intermediate 
C...stop-hadrons are formed.

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=200D0

C...Pick stop and gluino masses.
      PMST=80D0
      PMGL=50D0

C...Decide whether to form stop hadrons (1) or not (0).
      MSTOPH=1

C...Pick process.
      MSEL=0
      MSUB(261)=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...Force gluino mass = lightest neutralino (almost).
      RMSS(1)=PMGL
      RMSS(2)=PMGL
      RMSS(4)=1D4
C...Stop_1 mass.
      IMSS(5)=1
      RMSS(12)=PMST

C...Convenient shorthand.
      KFST=KSUSY1+6
      KCST=PYCOMP(KFST)
      KFGL=KSUSY1+21

C...Only allow decay channel for stop -> neutralino_1 + c. All else off.
      DO 100 IDC=MDCY(KCST,2),MDCY(KCST,2)+MDCY(KCST,3)-1
        IF(KFDP(IDC,1).EQ.KFGL+1.AND.KFDP(IDC,2).EQ.4) THEN
        ELSE
          MDME(IDC,1)=MIN(0,MDME(IDC,1))
        ENDIF
  100 CONTINUE

C...Initialize.
      CALL PYINIT('CMS','e+','e-',ECM)

C...Now overwrite neutralino in decay channel above with gluino.
C...Thereby we artificially create stop -> gluino + c channel.
      DO 110 IDC=MDCY(KCST,2),MDCY(KCST,2)+MDCY(KCST,3)-1
        IF(KFDP(IDC,1).EQ.KFGL+1.AND.KFDP(IDC,2).EQ.4) THEN
          KFDP(IDC,1)=KFGL 
        ELSE
        ENDIF
  110 CONTINUE

C...Put gluino stable.
      MDCY(PYCOMP(KFGL),1)=0
  
C...Show particle data.
c      CALL PYSTAT(2)
C      CALL PYLIST(12)

C...Switch off hadronization in normal PYEVNT call.
      MSTP(111)=0 

C...Histograms.
      CALL PYBOOK(1,'E-P-Q conservation check',100,0D0,0.0001D0)
      CALL PYBOOK(2,'energy in ISR photons',100,0D0,ECM-2D0*PMST)
      CALL PYBOOK(3,'charged hadron multiplicity (excl. R hadrons)',
     &50,-0.5D0,49.5D0)
      CALL PYBOOK(11,'stop squark momentum spectrum',100,0D0,100D0)
      CALL PYBOOK(12,'stop hadron momentum spectrum',100,0D0,100D0)
      CALL PYBOOK(13,'stop hadron mass spectrum, central',
     &100,PMST,PMST+2D0)
      CALL PYBOOK(14,'stop hadron mass spectrum, off central',
     &100,0D0,100D0)
      CALL PYBOOK(21,'gluino parton momentum spectrum',100,0D0,100D0)
      CALL PYBOOK(22,'gluino hadron momentum spectrum',100,0D0,100D0)
      CALL PYBOOK(23,'Gluino-hadron mass spectrum, central',
     &100,PMGL,PMGL+2D0)
      CALL PYBOOK(24,'Gluino-hadron mass spectrum, off central',
     &100,0D0,100D0)
 
C...Event generation loop.
      DO 200 IEV=1,NEV
      if(mod(iev,100).eq.0) write(*,*) 'begin event no', iev

C...Put stop stable.
        IF(MSTOPH.EQ.1) MDCY(KCST,1)=0

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...Restore stop decay capability.
        IF(MSTOPH.EQ.1) MDCY(KCST,1)=1

C...Now perform treatment of stop hadronization and decay.
        IF(MSTOPH.EQ.1) CALL PYSTFR

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

C...List first few events.
        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))
        IF(EPQ.GT.1D-4) CALL PYFILL(1,EPQ,1D0)
        IF(EPQ.GT.1D-4) NMERR=NMERR+1
        EGAMM=0D0
        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-8) THEN
             WRITE(*,*) 'problem',I,DM2
            IF(DM2.GT.1D-4) NMERR=NMERR+1
          ENDIF
          IF(K(I,1).EQ.1) THEN
            IF(K(I,2).EQ.22.AND.K(I,3).LE.4) EGAMM=EGAMM+P(I,4)
            IF(PYCHGE(K(I,2)).NE.0) NCH=NCH+1
          ENDIF      
          PABS=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)  
C...Statistics: stop.         
          IF((K(I,1).EQ.11.OR.K(I,1).EQ.12).AND.IABS(K(I,2)).EQ.KFST) 
     &    THEN
            CALL PYFILL(11,PABS,1D0)
          ENDIF
          IF(K(I,1).EQ.13.AND.IABS(K(I,2)).EQ.KFST) THEN
            CALL PYFILL(12,PABS,1D0)
            IF(K(I,2).GT.0) THEN
              IOTHER=MOD(K(I,4)/MSTU(5),MSTU(5))
            ELSE
              IOTHER=MOD(K(I,5)/MSTU(5),MSTU(5))
            ENDIF
            PMSTHD=P(I,5)+P(IOTHER,5)
            IF(PMSTHD.LT.PMST.OR.PMSTHD.GT.PMST+2D0) THEN
              CALL PYFILL(14,PMSTHD,1D0)
            ELSE
              CALL PYFILL(13,PMSTHD,1D0)
            ENDIF
          ENDIF
C...Statistics: gluino.
          IF(K(I,1).EQ.12.AND.K(I,2).EQ.KFGL) THEN
            CALL PYFILL(21,PABS,1D0)
          ENDIF
          IF((K(I,1).EQ.6.OR.K(I,1).EQ.7).AND.K(I,2).EQ.KFGL) THEN
            CALL PYFILL(22,PABS,1D0)
            IF(P(I,5).LT.PMGL.OR.P(I,5).GT.PMGL+2D0) THEN
              CALL PYFILL(24,P(I,5),1D0)
            ELSE
              CALL PYFILL(23,P(I,5),1D0)
            ENDIF
          ENDIF
 160    CONTINUE
        CALL PYFILL(2,EGAMM,1D0)
        CALL PYFILL(3,DBLE(NCH),1D0)
        IF(NMERR.GT.0) 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...PYSTFR
C...Fragments the string near to a stop, to form a stop-hadron, 
C...by producing a new q-qbar pair.
 
      SUBROUTINE PYSTFR
 
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),PSAV(5),IJOIN(2),IPOSST(10) 
 
C...Free parameter: max kinetic energy in gluino-hadron.
      PMKIN=0.5D0
 
C...Free parameter: part of stop mass that does not participate
C...in weak decay.
      PMINAC=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.
      KFST=KSUSY1+6
      KCST=PYCOMP(KFST)
      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).
      NOLD=N
      NSTOP=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(IABS(K(I,2)).EQ.KFST) THEN
            NSTOP=NSTOP+1
            IPOSST(NSTOP)=N
          ENDIF   
        ENDIF
  120 CONTINUE
      NTMP=N

C...Loop over (up to) two stops per event.
C...Identify position of stop (randomize order of treatment).
      IRNST=INT(1.5D0+PYR(0))
      DO 300 ISTOP=1,NSTOP
        IST=IPOSST(1)
        IF(NSTOP.EQ.2.AND.ISTOP.NE.IRNST) IST=IPOSST(2)

C...Identify range of partons on string the stop belongs to. 
        IMIN=IST+1
  140   IMIN=IMIN-1
        IF(K(IMIN-1,1).EQ.2) GOTO 140
        IMAX=IST-1
  150   IMAX=IMAX+1
        IF(K(IMAX,1).EQ.2) GOTO 150
        IOTHER=IMAX
        IF(IST.EQ.IMAX) IOTHER=IMIN  
 
C...Find mass of this stop-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 stop-hadron already formed.
        IF(PSUM(5).LE.P(IST,5)+P(IOTHER,5)+PMKIN) THEN
          DO 180 I=IMIN,IMAX
            K(I,1)=2
            IF(I.EQ.IMAX) K(I,1)=1
            IF(I.NE.IST) THEN
              DO 175 J=1,5
                P(IST,J)=P(IST,J)+P(I,J)
                P(I,J)=0D0
  175         CONTINUE
            ENDIF
  180     CONTINUE
          P(IST,5)=SQRT(MAX(0D0,P(IST,4)**2-P(IST,1)**2-P(IST,2)**2-
     &    P(IST,3)**2))
          GOTO 300
        ENDIF    

C...Else break string by production of new qqbar pair.
C...(Also diquarks allowed, but not popcorn.)
        INFLAV=ISIGN(4,K(IST,2))
        CALL PYDCYK(INFLAV,0,KFSAV,KFDUM)
        KFSAV=ISIGN(MOD(IABS(KFSAV),10000),KFSAV)
        MSTJ(93)=1 
        PMSAV=PYMASS(KFSAV)         

C...Mass of stop-hadron.
        PMSSAV=P(IST,5)
        PMSHAD=P(IST,5)+PMSAV

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

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

C...Boost to rest frame of system, and align stop along +z axis.
          CALL PYROBO(IST,IST,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(IST,1),P(IST,2))
          CALL PYROBO(IST,IST,0D0,-PHI,0D0,0D0,0D0)
          CALL PYROBO(IREC,IREC,0D0,-PHI,0D0,0D0,0D0)
          THETA=PYANGL(P(IST,3),P(IST,1)) 
          CALL PYROBO(IST,IST,-THETA,0D0,0D0,0D0,0D0)
          CALL PYROBO(IREC,IREC,-THETA,0D0,0D0,0D0,0D0)

C...Calculate new kinematics in this frame, for desired stop hadron mass.
          ETOT=P(IST,4)+P(IREC,4)
          PMREC=P(IREC,5)
          IF(K(IREC,2).NE.21.AND.IABS(K(IREC,2)).NE.KFST) THEN
            MSTJ(93)=1 
            PMREC=PYMASS(K(IREC,2))         
          ENDIF 
          IF(ETOT.GT.PMSHAD+PMREC) THEN
            IFAIL=0
            PZNEW=0.5D0*SQRT(MAX(0D0,(ETOT**2-PMSHAD**2-PMREC**2)**2-
     &      4D0*PMSHAD**2*PMREC**2))/ETOT
            P(IST,3)=PZNEW
            P(IST,4)=SQRT(PZNEW**2+PMSHAD**2)
            P(IST,5)=PMSHAD
            P(IREC,3)=-PZNEW
            P(IREC,4)=SQRT(PZNEW**2+PMREC**2)
            P(IREC,5)=PMREC

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

C...Bost back to lab frame.
          CALL PYROBO(IST,IST,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.)
          IF(IFAIL.EQ.1) THEN
            IF(IST.EQ.IMIN.AND.IREC.LT.IMAX) THEN
              IREC=IREC+1
              GOTO 200
            ELSEIF(IST.EQ.IMAX.AND.IREC.GT.IMIN) THEN
              IREC=IREC-1
              GOTO 200
            ENDIF
          ENDIF

C...New slots at end of record for squark + new antiquark.
        DO 230 J=1,5
          K(N+1,J)=0
          P(N+1,J)=P(IST,J)
          V(N+1,J)=V(IST,J)
          K(N+2,J)=0
          P(N+2,J)=0D0
          V(N+2,J)=V(IST,J)
  230   CONTINUE
 
C...Status and code of these slots.
        K(N+1,1)=2
        K(N+2,1)=1
        K(N+1,2)=K(IST,2)
        K(N+2,2)=KFSAV
        K(N+1,3)=K(IST,3)
        K(N+2,3)=K(IST,3)
        N=N+2
        
C...Code and momentum of new string endpoint.
        K(IST,2)=-KFSAV
        DO 240 J=1,5
          P(IST,J)=PSAV(J)
  240   CONTINUE
 
C...End of loop over two stops.
  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
      NNOW=N

C...Check that no low-mass system of diquark-antidiquark kind,
C...or very low-mass of any kind.
      KFBEG=0
      DO 332 J=1,5
        PSUM(J)=0D0
  332 CONTINUE
      DO 338 I=NOLD+1,NNOW
        DO 334 J=1,4
          PSUM(J)=PSUM(J)+P(I,J)
  334   CONTINUE
        IF(KFBEG.EQ.0) THEN
          KFBEG=IABS(K(I,2))
          MSTJ(93)=1 
          PSUM(5)=PSUM(5)+PYMASS(K(I,2))         
        ELSEIF(K(I,1).EQ.1) THEN
          KFEND=IABS(K(I,2))
          MSTJ(93)=1 
          PSUM(5)=PSUM(5)+PYMASS(K(I,2))         
          DELTA=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-
     &    PSUM(3)**2))-PSUM(5)
          IF(KFBEG.GT.10.AND.KFBEG.LT.10000.AND.KFEND.GT.10.AND.
     &    KFEND.LT.10000.AND.DELTA.LT.PARJ(32)) GOTO 90
          IF(DELTA.LT.0D0) GOTO 90
          KFBEG=0
          DO 336 J=1,5
            PSUM(J)=0D0
  336     CONTINUE
        ENDIF
  338 CONTINUE

C...Finished with stop hadronization. Restore baryon production model.
      MSTJ(12)=MSTJ12

C...Now hadronize small string piece. Some cheating to allow sensible
C...momentum shuffling.
      IF(N-NOLD.GT.4) THEN
        MDCY(8,1)=0
        DO 340 I=NOLD+1,N
          IF(IABS(K(I,2)).EQ.KFST) THEN
            K(I,1)=1
            K(I,2)=ISIGN(8,K(I,2))
          ELSEIF(IABS(K(I-1,2)).EQ.8) THEN
            K(I,1)=11
          ENDIF
  340   CONTINUE
        MSTJ16=MSTJ(16)
        MSTJ(16)=0
        CALL PYEXEC
        MSTJ(16)=MSTJ16
        DO 350 I=NOLD+1,NNOW
          IF(IABS(K(I,2)).EQ.8) THEN
            K(I,1)=2
            K(I,2)=ISIGN(KFST,K(I,2))
          ELSEIF(IABS(K(I-1,2)).EQ.KFST.AND.K(I-1,1).EQ.2) THEN
            K(I,1)=1
          ENDIF
  350   CONTINUE
      ENDIF

C...Split momentum inside stop-hadron.
      DO 370 I=NOLD+1,NNOW
        IF(IABS(K(I,2)).EQ.KFST) THEN
          ISPEC=I+1
          IF(K(I,1).EQ.1) ISPEC=I-1
          MSTJ(93)=1 
          PMSPEC=PYMASS(K(ISPEC,2))         
          FAC=(PMSPEC+PMINAC)/P(I,5)
          DO 360 J=1,5
            P(ISPEC,J)=FAC*P(I,J)
            P(I,J)=(1D0-FAC)*P(I,J)
  360     CONTINUE
    
C...Hook up colours and decay stop to gluino.
          IJOIN(1)=I
          IJOIN(2)=ISPEC
          CALL PYJOIN(2,IJOIN)
          CALL PYRESD(I)
        ENDIF      
  370 CONTINUE

C...Rearrange partons along string.
      MSTJ(14)=-1
      CALL PYPREP(0)
      MSTJ(14)=1
          
      RETURN
      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
          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.)
      DO 340 I=1,N
        IF(K(I,1).EQ.16.OR.K(I,1).EQ.17) K(I,1)=K(I,1)-10
  340 CONTINUE

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
        WRITE(*,*) '     ...give it new try...'
        MSTU(23)=MSTU(23)-1
        GOTO 90
      ELSEIF(MSTU(24).NE.0) THEN
        WRITE(*,*) '     ...but still fail after repeated attempts!'
      ELSEIF(MSTU(24).EQ.0.AND.LOOP.GT.1) THEN
        WRITE(*,*) '     ...and now it worked!'
      ENDIF

C...Finished. Restore baryon production model.
      MSTJ(12)=MSTJ12
      
      RETURN
      END