pythiaeRHIC/pythia/pyboei.f

0001
0002 C*********************************************************************
0003
0004 C...PYBOEI
0005 C...Modifies an event so as to approximately take into account
0006 C...Bose-Einstein effects according to a simple phenomenological
0007 C...parametrization.
0008
0009       SUBROUTINE PYBOEI(NSAV)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Parameter statement to help give large particle numbers.
0016       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0017      &KEXCIT=4000000,KDIMEN=5000000)
0018 C...Commonblocks.
0019       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0020       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0021       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0022       COMMON/PYINT1/MINT(400),VINT(400)
0023       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYINT1/
0024 C...Local arrays and data.
0025       DIMENSION DPS(4),KFBE(9),NBE(0:10),BEI(100),BEI3(100),
0026      &BEIW(100),BEI3W(100)
0027       DATA KFBE/211,-211,111,321,-321,130,310,221,331/
0028 C...Statement function: squared invariant mass.
0029       SDIP(I,J)=((P(I,4)+P(J,4))**2-(P(I,3)+P(J,3))**2-
0030      &(P(I,2)+P(J,2))**2-(P(I,1)+P(J,1))**2)
0031
0032 C...Boost event to overall CM frame. Calculate CM energy.
0033       IF((MSTJ(51).NE.1.AND.MSTJ(51).NE.2).OR.N-NSAV.LE.1) RETURN
0034       DO 100 J=1,4
0035         DPS(J)=0D0
0036   100 CONTINUE
0037       DO 120 I=1,N
0038         KFA=IABS(K(I,2))
0039         IF(K(I,1).LE.10.AND.((KFA.GT.10.AND.KFA.LE.20).OR.KFA.EQ.22)
0040      &  .AND.K(I,3).GT.0) THEN
0041           KFMA=IABS(K(K(I,3),2))
0042           IF(KFMA.GT.10.AND.KFMA.LE.80) K(I,1)=-K(I,1)
0043         ENDIF
0044         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
0045         DO 110 J=1,4
0046           DPS(J)=DPS(J)+P(I,J)
0047   110   CONTINUE
0048   120 CONTINUE
0049       CALL PYROBO(0,0,0D0,0D0,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
0050      &-DPS(3)/DPS(4))
0051       PECM=0D0
0052       DO 130 I=1,N
0053         IF(K(I,1).GE.1.AND.K(I,1).LE.10) PECM=PECM+P(I,4)
0054   130 CONTINUE
0055
0056 C...Check if we have separated strings
0057
0058 C...Reserve copy of particles by species at end of record.
0059       IWP=0
0060       IWN=0
0061       NBE(0)=N+MSTU(3)
0062       NMAX=NBE(0)
0063       SMMIN=PECM
0064       DO 190 IBE=1,MIN(10,MSTJ(52)+1)
0065         NBE(IBE)=NBE(IBE-1)
0066         DO 180 I=NSAV+1,N
0067           IF(IBE.EQ.MIN(10,MSTJ(52)+1)) THEN
0068             DO 140 IIBE=1,IBE-1
0069               IF(K(I,2).EQ.KFBE(IIBE)) GOTO 180
0070   140       CONTINUE
0071           ELSE
0072             IF(K(I,2).NE.KFBE(IBE)) GOTO 180
0073           ENDIF
0074           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 180
0075           IF(NBE(IBE).GE.MSTU(4)-MSTU(32)-5) THEN
0076             CALL PYERRM(11,'(PYBOEI:) no more memory left in PYJETS')
0077             RETURN
0078           ENDIF
0079           NBE(IBE)=NBE(IBE)+1
0080           NMAX=NBE(IBE)
0081           K(NBE(IBE),1)=I
0082           K(NBE(IBE),2)=0
0083           K(NBE(IBE),3)=0
0084           K(NBE(IBE),4)=0
0085           K(NBE(IBE),5)=0
0086           P(NBE(IBE),1)=0.0D0
0087           P(NBE(IBE),2)=0.0D0
0088           P(NBE(IBE),3)=0.0D0
0089           P(NBE(IBE),4)=0.0D0
0090           P(NBE(IBE),5)=0.0D0
0091           SMMIN=MIN(SMMIN,P(I,5))
0092 C...Check if particles comes from different W's or Z's
0093           IF((MSTJ(53).NE.0.OR.MSTJ(56).GT.0).AND.MINT(32).EQ.0) THEN
0094             IM=I
0095   150       IF(K(IM,3).GT.0) THEN
0096               IM=K(IM,3)
0097               IF(ABS(K(IM,2)).NE.24.AND.K(IM,2).NE.23) GOTO 150
0098               K(NBE(IBE),5)=IM
0099               IF(IWP.EQ.0.AND.K(IM,2).EQ.24) IWP=IM
0100               IF(IWN.EQ.0.AND.K(IM,2).EQ.-24) IWN=IM
0101               IF(IWP.EQ.0.AND.K(IM,2).EQ.23) IWP=IM
0102               IF(IWN.EQ.0.AND.K(IM,2).EQ.23.AND.IM.NE.IWP) IWN=IM
0103             ENDIF
0104           ENDIF
0105 C...Check if particles comes from different strings.
0106           IF(PARJ(94).GT.0.0D0) THEN
0107             IM=I
0108   160       IF(K(IM,3).GT.0) THEN
0109               IM=K(IM,3)
0110               IF(K(IM,2).NE.92.AND.K(IM,2).NE.91) GOTO 160
0111               K(NBE(IBE),5)=IM
0112             ENDIF
0113           ENDIF
0114           DO 170 J=1,3
0115             P(NBE(IBE),J)=0D0
0116             V(NBE(IBE),J)=0D0
0117   170     CONTINUE
0118           P(NBE(IBE),5)=-1.0D0
0119   180   CONTINUE
0120   190 CONTINUE
0121       IF(NBE(MIN(9,MSTJ(52)))-NBE(0).LE.1) GOTO 510
0122
0123 C...Calculate separation between W+ and W- or between two Z0's.
0124 C...No separation if there has been re-connections.
0125       SIGW=PARJ(93)
0126       IF(IWP.GT.0.AND.IWN.GT.0.AND.MSTJ(56).GT.0.AND.MINT(32).EQ.0) THEN
0127         IF(K(IWP,2).EQ.23) THEN
0128           DMW=PMAS(23,1)
0129           DGW=PMAS(23,2)
0130         ELSE
0131           DMW=PMAS(24,1)
0132           DGW=PMAS(24,2)
0133         ENDIF
0134         DMP=P(IWP,5)
0135         DMN=P(IWN,5)
0136         TAUPD=DMP/SQRT((DMP**2-DMW**2)**2+(DGW*(DMP**2)/DMW)**2)
0137         TAUND=DMN/SQRT((DMN**2-DMW**2)**2+(DGW*(DMN**2)/DMW)**2)
0138         TAUP=-TAUPD*LOG(PYR(IDUM))
0139         TAUN=-TAUND*LOG(PYR(IDUM))
0140         DXP=TAUP*PYP(IWP,8)/DMP
0141         DXN=TAUN*PYP(IWN,8)/DMN
0142         DX=DXP+DXN
0143         SIGW=1.0D0/(1.0D0/PARJ(93)+REAL(MSTJ(56))*DX)
0144         IF(PARJ(94).LT.0.0D0) SIGW=1.0D0/(1.0D0/SIGW-1.0D0/PARJ(94))
0145       ENDIF
0146
0147 C...Add separation between strings.
0148       IF(PARJ(94).GT.0.0D0) THEN
0149         SIGW=1.0D0/(1.0D0/SIGW+1.0D0/PARJ(94))
0150         IWP=-1
0151         IWN=-1
0152       ENDIF
0153
0154       IF(MSTJ(57).EQ.1.AND.MSTJ(54).LT.0) THEN
0155         DO 220 IBE=1,MIN(9,MSTJ(52))
0156           DO 210 I1M=NBE(IBE-1)+1,NBE(IBE)
0157             Q2MIN=PECM**2
0158             I1=K(I1M,1)
0159             DO 200 I2M=NBE(IBE-1)+1,NBE(IBE)
0160               IF(I2M.EQ.I1M) GOTO 200
0161               I2=K(I2M,1)
0162               Q2=(P(I1,4)+P(I2,4))**2-(P(I1,1)+P(I2,1))**2-
0163      &        (P(I1,2)+P(I2,2))**2-(P(I1,3)+P(I2,3))**2-
0164      &        (P(I1,5)+P(I2,5))**2
0165               IF(Q2.GT.0.0D0.AND.Q2.LT.Q2MIN) THEN
0166                 Q2MIN=Q2
0167               ENDIF
0168   200       CONTINUE
0169             P(I1M,5)=Q2MIN
0170   210     CONTINUE
0171   220   CONTINUE
0172       ENDIF
0173
0174 C...Tabulate integral for subsequent momentum shift.
0175       DO 400 IBE=1,MIN(9,MSTJ(52))
0176         IF(IBE.NE.1.AND.IBE.NE.4.AND.IBE.LE.7) GOTO 270
0177         IF(IBE.EQ.1.AND.MAX(NBE(1)-NBE(0),NBE(2)-NBE(1),NBE(3)-NBE(2))
0178      &  .LE.1) GOTO 270
0179         IF(IBE.EQ.4.AND.MAX(NBE(4)-NBE(3),NBE(5)-NBE(4),NBE(6)-NBE(5),
0180      &  NBE(7)-NBE(6)).LE.1) GOTO 270
0181         IF(IBE.GE.8.AND.NBE(IBE)-NBE(IBE-1).LE.1) GOTO 270
0182         IF(IBE.EQ.1) PMHQ=2D0*PYMASS(211)
0183         IF(IBE.EQ.4) PMHQ=2D0*PYMASS(321)
0184         IF(IBE.EQ.8) PMHQ=2D0*PYMASS(221)
0185         IF(IBE.EQ.9) PMHQ=2D0*PYMASS(331)
0186         QDEL=0.1D0*MIN(PMHQ,PARJ(93))
0187         QDEL3=0.1D0*MIN(PMHQ,PARJ(93)*3.0D0)
0188         QDELW=0.1D0*MIN(PMHQ,SIGW)
0189         QDEL3W=0.1D0*MIN(PMHQ,SIGW*3.0D0)
0190         IF(MSTJ(51).EQ.1) THEN
0191           NBIN=MIN(100,NINT(9D0*PARJ(93)/QDEL))
0192           NBIN3=MIN(100,NINT(27D0*PARJ(93)/QDEL3))
0193           NBINW=MIN(100,NINT(9D0*SIGW/QDELW))
0194           NBIN3W=MIN(100,NINT(27D0*SIGW/QDEL3W))
0195           BEEX=EXP(0.5D0*QDEL/PARJ(93))
0196           BEEX3=EXP(0.5D0*QDEL3/(3.0D0*PARJ(93)))
0197           BEEXW=EXP(0.5D0*QDELW/SIGW)
0198           BEEX3W=EXP(0.5D0*QDEL3W/(3.0D0*SIGW))
0199           BERT=EXP(-QDEL/PARJ(93))
0200           BERT3=EXP(-QDEL3/(3.0D0*PARJ(93)))
0201           BERTW=EXP(-QDELW/SIGW)
0202           BERT3W=EXP(-QDEL3W/(3.0D0*SIGW))
0203         ELSE
0204           NBIN=MIN(100,NINT(3D0*PARJ(93)/QDEL))
0205           NBIN3=MIN(100,NINT(9D0*PARJ(93)/QDEL3))
0206           NBINW=MIN(100,NINT(3D0*SIGW/QDELW))
0207           NBIN3W=MIN(100,NINT(9D0*SIGW/QDEL3W))
0208         ENDIF
0209         DO 230 IBIN=1,NBIN
0210           QBIN=QDEL*(IBIN-0.5D0)
0211           BEI(IBIN)=QDEL*(QBIN**2+QDEL**2/12D0)/SQRT(QBIN**2+PMHQ**2)
0212           IF(MSTJ(51).EQ.1) THEN
0213             BEEX=BEEX*BERT
0214             BEI(IBIN)=BEI(IBIN)*BEEX
0215           ELSE
0216             BEI(IBIN)=BEI(IBIN)*EXP(-(QBIN/PARJ(93))**2)
0217           ENDIF
0218           IF(IBIN.GE.2) BEI(IBIN)=BEI(IBIN)+BEI(IBIN-1)
0219   230   CONTINUE
0220         DO 240 IBIN=1,NBIN3
0221           QBIN=QDEL3*(IBIN-0.5D0)
0222           BEI3(IBIN)=QDEL3*(QBIN**2+QDEL3**2/12D0)/SQRT(QBIN**2+PMHQ**2)
0223           IF(MSTJ(51).EQ.1) THEN
0224             BEEX3=BEEX3*BERT3
0225             BEI3(IBIN)=BEI3(IBIN)*BEEX3
0226           ELSE
0227             BEI3(IBIN)=BEI3(IBIN)*EXP(-(QBIN/(3.0D0*PARJ(93)))**2)
0228           ENDIF
0229           IF(IBIN.GE.2) BEI3(IBIN)=BEI3(IBIN)+BEI3(IBIN-1)
0230   240   CONTINUE
0231         DO 250 IBIN=1,NBINW
0232           QBIN=QDELW*(IBIN-0.5D0)
0233           BEIW(IBIN)=QDELW*(QBIN**2+QDELW**2/12D0)/SQRT(QBIN**2+PMHQ**2)
0234           IF(MSTJ(51).EQ.1) THEN
0235             BEEXW=BEEXW*BERTW
0236             BEIW(IBIN)=BEIW(IBIN)*BEEXW
0237           ELSE
0238             BEIW(IBIN)=BEIW(IBIN)*EXP(-(QBIN/SIGW)**2)
0239           ENDIF
0240           IF(IBIN.GE.2) BEIW(IBIN)=BEIW(IBIN)+BEIW(IBIN-1)
0241   250   CONTINUE
0242         DO 260 IBIN=1,NBIN3W
0243           QBIN=QDEL3W*(IBIN-0.5D0)
0244           BEI3W(IBIN)=QDEL3W*(QBIN**2+QDEL3W**2/12D0)/
0245      &    SQRT(QBIN**2+PMHQ**2)
0246           IF(MSTJ(51).EQ.1) THEN
0247             BEEX3W=BEEX3W*BERT3W
0248             BEI3W(IBIN)=BEI3W(IBIN)*BEEX3W
0249           ELSE
0250             BEI3W(IBIN)=BEI3W(IBIN)*EXP(-(QBIN/(3.0D0*SIGW))**2)
0251           ENDIF
0252           IF(IBIN.GE.2) BEI3W(IBIN)=BEI3W(IBIN)+BEI3W(IBIN-1)
0253   260   CONTINUE
0254
0255 C...Loop through particle pairs and find old relative momentum.
0256   270   DO 390 I1M=NBE(IBE-1)+1,NBE(IBE)-1
0257           I1=K(I1M,1)
0258           DO 380 I2M=I1M+1,NBE(IBE)
0259             IF(MSTJ(53).EQ.1.AND.K(I1M,5).NE.K(I2M,5)) GOTO 380
0260             IF(MSTJ(53).EQ.2.AND.K(I1M,5).EQ.K(I2M,5)) GOTO 380
0261             I2=K(I2M,1)
0262             Q2OLD=(P(I1,4)+P(I2,4))**2-(P(I1,1)+P(I2,1))**2-(P(I1,2)+
0263      &      P(I2,2))**2-(P(I1,3)+P(I2,3))**2-(P(I1,5)+P(I2,5))**2
0264             IF(Q2OLD.LE.0.0D0) GOTO 380
0265             QOLD=SQRT(Q2OLD)
0266
0267 C...Calculate new relative momentum.
0268             QMOV=0.0D0
0269             QMOV3=0.0D0
0270             QMOVW=0.0D0
0271             QMOV3W=0.0D0
0272             IF(QOLD.LT.1D-3*QDEL) THEN
0273               GOTO 280
0274             ELSEIF(QOLD.LE.QDEL) THEN
0275               QMOV=QOLD/3D0
0276             ELSEIF(QOLD.LT.(NBIN-0.1D0)*QDEL) THEN
0277               RBIN=QOLD/QDEL
0278               IBIN=RBIN
0279               RINP=(RBIN**3-IBIN**3)/(3*IBIN*(IBIN+1)+1)
0280               QMOV=(BEI(IBIN)+RINP*(BEI(IBIN+1)-BEI(IBIN)))*
0281      &        SQRT(Q2OLD+PMHQ**2)/Q2OLD
0282             ELSE
0283               QMOV=BEI(NBIN)*SQRT(Q2OLD+PMHQ**2)/Q2OLD
0284             ENDIF
0285   280       Q2NEW=Q2OLD*(QOLD/(QOLD+3D0*PARJ(92)*QMOV))**(2D0/3D0)
0286             IF(QOLD.LT.1D-3*QDEL3) THEN
0287               GOTO 290
0288             ELSEIF(QOLD.LE.QDEL3) THEN
0289               QMOV3=QOLD/3D0
0290             ELSEIF(QOLD.LT.(NBIN3-0.1D0)*QDEL3) THEN
0291               RBIN3=QOLD/QDEL3
0292               IBIN3=RBIN3
0293               RINP3=(RBIN3**3-IBIN3**3)/(3*IBIN3*(IBIN3+1)+1)
0294               QMOV3=(BEI3(IBIN3)+RINP3*(BEI3(IBIN3+1)-BEI3(IBIN3)))*
0295      &        SQRT(Q2OLD+PMHQ**2)/Q2OLD
0296             ELSE
0297               QMOV3=BEI3(NBIN3)*SQRT(Q2OLD+PMHQ**2)/Q2OLD
0298             ENDIF
0299   290       Q2NEW3=Q2OLD*(QOLD/(QOLD+3D0*PARJ(92)*QMOV3))**(2D0/3D0)
0300             RSCALE=1.0D0
0301             IF(MSTJ(54).EQ.2)
0302      &      RSCALE=1.0D0-EXP(-(QOLD/(2D0*PARJ(93)))**2)
0303             IF((IWP.NE.-1.AND.MSTJ(56).LE.0).OR.IWP.EQ.0.OR.IWN.EQ.0.OR.
0304      &      K(I1M,5).EQ.K(I2M,5)) GOTO 320
0305
0306             IF(QOLD.LT.1D-3*QDELW) THEN
0307               GOTO 300
0308             ELSEIF(QOLD.LE.QDELW) THEN
0309               QMOVW=QOLD/3D0
0310             ELSEIF(QOLD.LT.(NBINW-0.1D0)*QDELW) THEN
0311               RBINW=QOLD/QDELW
0312               IBINW=RBINW
0313               RINPW=(RBINW**3-IBINW**3)/(3*IBINW*(IBINW+1)+1)
0314               QMOVW=(BEIW(IBINW)+RINPW*(BEIW(IBINW+1)-BEIW(IBINW)))*
0315      &        SQRT(Q2OLD+PMHQ**2)/Q2OLD
0316             ELSE
0317               QMOVW=BEIW(NBINW)*SQRT(Q2OLD+PMHQ**2)/Q2OLD
0318             ENDIF
0319   300       Q2NEW=Q2OLD*(QOLD/(QOLD+3D0*PARJ(92)*QMOVW))**(2D0/3D0)
0320             IF(QOLD.LT.1D-3*QDEL3W) THEN
0321               GOTO 310
0322             ELSEIF(QOLD.LE.QDEL3W) THEN
0323               QMOV3W=QOLD/3D0
0324             ELSEIF(QOLD.LT.(NBIN3W-0.1D0)*QDEL3W) THEN
0325               RBIN3W=QOLD/QDEL3W
0326               IBIN3W=RBIN3W
0327               RINP3W=(RBIN3W**3-IBIN3W**3)/(3*IBIN3W*(IBIN3W+1)+1)
0328               QMOV3W=(BEI3W(IBIN3W)+RINP3W*(BEI3W(IBIN3W+1)-
0329      &        BEI3W(IBIN3W)))*SQRT(Q2OLD+PMHQ**2)/Q2OLD
0330             ELSE
0331               QMOV3W=BEI3W(NBIN3W)*SQRT(Q2OLD+PMHQ**2)/Q2OLD
0332             ENDIF
0333   310       Q2NEW3=Q2OLD*(QOLD/(QOLD+3D0*PARJ(92)*QMOV3W))**(2D0/3D0)
0334             IF(MSTJ(54).EQ.2)
0335      &      RSCALE=1.0D0-EXP(-(QOLD/(2D0*SIGW))**2)
0336
0337   320       CALL PYBESQ(I1,I2,NMAX,Q2OLD,Q2NEW)
0338             DO 330 J=1,3
0339               P(I1M,J)=P(I1M,J)+P(NMAX+1,J)
0340               P(I2M,J)=P(I2M,J)+P(NMAX+2,J)
0341   330       CONTINUE
0342             IF(MSTJ(54).GE.1) THEN
0343               CALL PYBESQ(I1,I2,NMAX,Q2OLD,Q2NEW3)
0344               DO 340 J=1,3
0345                 V(I1M,J)=V(I1M,J)+P(NMAX+1,J)*RSCALE
0346                 V(I2M,J)=V(I2M,J)+P(NMAX+2,J)*RSCALE
0347   340         CONTINUE
0348             ELSEIF(MSTJ(54).LE.-1) THEN
0349               EDEL=P(I1,4)+P(I2,4)-
0350      &        SQRT(MAX(Q2NEW-Q2OLD+(P(I1,4)+P(I2,4))**2,0.0D0))
0351               A2=(P(I1,1)-P(I2,1))**2+(P(I1,2)-P(I2,2))**2+
0352      &        (P(I1,3)-P(I2,3))**2
0353               WMAX=-1.0D20
0354               MI3=0
0355               MI4=0
0356               S12=SDIP(I1,I2)
0357               SM1=(P(I1,5)+SMMIN)**2
0358               DO 360 I3M=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1))
0359                 IF(I3M.EQ.I1M.OR.I3M.EQ.I2M) GOTO 360
0360                 IF(MSTJ(53).EQ.1.AND.K(I3M,5).NE.K(I1M,5)) GOTO 360
0361                 IF(MSTJ(53).EQ.-2.AND.K(I1M,5).EQ.K(I2M,5).AND.
0362      &          K(I3M,5).NE.K(I1M,5)) GOTO 360
0363                 I3=K(I3M,1)
0364                 IF(K(I3,2).EQ.K(I1,2)) GOTO 360
0365                 S13=SDIP(I1,I3)
0366                 S23=SDIP(I2,I3)
0367                 SM3=(P(I3,5)+SMMIN)**2
0368                 IF(MSTJ(54).EQ.-2) THEN
0369                   WI=(MIN(S12*SM3,S13*MIN(SM1,SM3),
0370      &            S23*MIN(SM1,SM3))*SM1)
0371                 ELSE
0372                   WI=((P(I1,4)+P(I2,4)+P(I3,4))**2-
0373      &            (P(I1,3)+P(I2,3)+P(I3,3))**2-
0374      &            (P(I1,2)+P(I2,2)+P(I3,2))**2-
0375      &            (P(I1,1)+P(I2,1)+P(I3,1))**2)
0376                 ENDIF
0377                 IF(MSTJ(57).EQ.1.AND.P(I3M,5).GT.0) THEN
0378                   IF (WMAX*WI.GE.(1.0D0-EXP(-P(I3M,5)/(PARJ(93)**2))))
0379      &                 GOTO 360
0380                 ELSE
0381                   IF(WMAX*WI.GE.1.0) GOTO 360
0382                 ENDIF
0383                 DO 350 I4M=I3M+1,NBE(MIN(10,MSTJ(52)+1))
0384                   IF(I4M.EQ.I1M.OR.I4M.EQ.I2M) GOTO 350
0385                   IF(MSTJ(53).EQ.1.AND.K(I4M,5).NE.K(I1M,5)) GOTO 350
0386                   IF(MSTJ(53).EQ.-2.AND.K(I1M,5).EQ.K(I2M,5).AND.
0387      &            K(I4M,5).NE.K(I1M,5)) GOTO 350
0388                   I4=K(I4M,1)
0389                   IF(K(I3,2).EQ.K(I4,2).OR.K(I4,2).EQ.K(I1,2))
0390      &            GOTO 350
0391                   IF((P(I3,4)+P(I4,4)+EDEL)**2.LT.
0392      &            (P(I3,1)+P(I4,1))**2+(P(I3,2)+P(I4,2))**2+
0393      &            (P(I3,3)+P(I4,3))**2+(P(I3,5)+P(I4,5))**2)
0394      &            GOTO 350
0395                   IF(MSTJ(54).EQ.-2) THEN
0396                     S14=SDIP(I1,I4)
0397                     S24=SDIP(I2,I4)
0398                     S34=SDIP(I3,I4)
0399                     W=S12*MIN(MIN(S23,S24),MIN(S13,S14))*S34
0400                     W=MIN(W,S13*MIN(MIN(S23,S34),S12)*S24)
0401                     W=MIN(W,S14*MIN(MIN(S24,S34),S12)*S23)
0402                     W=MIN(W,MIN(S23,S24)*S13*S14)
0403                     W=1.0D0/W
0404                   ELSE
0405 C...weight=1-cos(theta)/mtot2
0406                     S1234=(P(I1,4)+P(I2,4)+P(I3,4)+P(I4,4))**2-
0407      &              (P(I1,3)+P(I2,3)+P(I3,3)+P(I4,3))**2-
0408      &              (P(I1,2)+P(I2,2)+P(I3,2)+P(I4,2))**2-
0409      &              (P(I1,1)+P(I2,1)+P(I3,1)+P(I4,1))**2
0410                     W=1.0D0/S1234
0411                     IF(W.LE.WMAX) GOTO 350
0412                   ENDIF
0413                   IF(MSTJ(57).EQ.1.AND.P(I3M,5).GT.0)
0414      &            W=W*(1.0D0-EXP(-P(I3M,5)/(PARJ(93)**2)))
0415                   IF(MSTJ(57).EQ.1.AND.P(I4M,5).GT.0)
0416      &            W=W*(1.0D0-EXP(-P(I4M,5)/(PARJ(93)**2)))
0417                   IF(W.LE.WMAX) GOTO 350
0418                   MI3=I3M
0419                   MI4=I4M
0420                   WMAX=W
0421   350           CONTINUE
0422   360         CONTINUE
0423               IF(MI4.EQ.0) GOTO 380
0424               I3=K(MI3,1)
0425               I4=K(MI4,1)
0426               EOLD=P(I3,4)+P(I4,4)
0427               ENEW=EOLD+EDEL
0428               P2=(P(I3,1)+P(I4,1))**2+(P(I3,2)+P(I4,2))**2+
0429      &        (P(I3,3)+P(I4,3))**2
0430               Q2NEWP=MAX(0.0D0,ENEW**2-P2-(P(I3,5)+P(I4,5))**2)
0431               Q2OLDP=MAX(0.0D0,EOLD**2-P2-(P(I3,5)+P(I4,5))**2)
0432               CALL PYBESQ(I3,I4,NMAX,Q2OLDP,Q2NEWP)
0433               DO 370 J=1,3
0434                 V(MI3,J)=V(MI3,J)+P(NMAX+1,J)
0435                 V(MI4,J)=V(MI4,J)+P(NMAX+2,J)
0436   370         CONTINUE
0437             ENDIF
0438   380     CONTINUE
0439   390   CONTINUE
0440   400 CONTINUE
0441
0442 C...Shift momenta and recalculate energies.
0443       ESUMP=0.0D0
0444       ESUM=0.0D0
0445       PROD=0.0D0
0446       DO 430 IM=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1))
0447         I=K(IM,1)
0448         ESUMP=ESUMP+P(I,4)
0449         DO 410 J=1,3
0450           P(I,J)=P(I,J)+P(IM,J)
0451   410   CONTINUE
0452         P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0453         ESUM=ESUM+P(I,4)
0454         DO 420 J=1,3
0455           PROD=PROD+V(IM,J)*P(I,J)/P(I,4)
0456   420   CONTINUE
0457   430 CONTINUE
0458
0459       PARJ(96)=0.0D0
0460       IF(MSTJ(54).NE.0.AND.PROD.NE.0.0D0) THEN
0461   440   ALPHA=(ESUMP-ESUM)/PROD
0462         PARJ(96)=PARJ(96)+ALPHA
0463         PROD=0.0D0
0464         ESUM=0.0D0
0465         DO 470 IM=NBE(0)+1,NBE(MIN(10,MSTJ(52)+1))
0466           I=K(IM,1)
0467           DO 450 J=1,3
0468             P(I,J)=P(I,J)+ALPHA*V(IM,J)
0469   450     CONTINUE
0470           P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0471           ESUM=ESUM+P(I,4)
0472           DO 460 J=1,3
0473             PROD=PROD+V(IM,J)*P(I,J)/P(I,4)
0474   460     CONTINUE
0475   470   CONTINUE
0476         IF(PROD.NE.0.0D0.AND.ABS(ESUMP-ESUM)/PECM.GT.0.00001D0)
0477      &  GOTO 440
0478       ENDIF
0479
0480 C...Rescale all momenta for energy conservation.
0481       PES=0D0
0482       PQS=0D0
0483       DO 480 I=1,N
0484         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 480
0485         PES=PES+P(I,4)
0486         PQS=PQS+P(I,5)**2/P(I,4)
0487   480 CONTINUE
0488       PARJ(95)=PES-PECM
0489       FAC=(PECM-PQS)/(PES-PQS)
0490       DO 500 I=1,N
0491         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 500
0492         DO 490 J=1,3
0493           P(I,J)=FAC*P(I,J)
0494   490   CONTINUE
0495         P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0496   500 CONTINUE
0497
0498 C...Boost back to correct reference frame.
0499   510 CALL PYROBO(0,0,0D0,0D0,DPS(1)/DPS(4),DPS(2)/DPS(4),DPS(3)/DPS(4))
0500       DO 520 I=1,N
0501         IF(K(I,1).LT.0) K(I,1)=-K(I,1)
0502   520 CONTINUE
0503
0504       RETURN
0505       END