pythiaeRHIC/pythia/pyremn.f

0001
0002 C*********************************************************************
0003
0004 C...PYREMN
0005 C...Adds on target remnants (one or two from each side) and
0006 C...includes primordial kT for hadron beams.
0007
0008       SUBROUTINE PYREMN(IPU1,IPU2)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0019       COMMON/PYINT1/MINT(400),VINT(400)
0020       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/
0021 C...Local arrays.
0022       DIMENSION KFLCH(2),KFLSP(2),CHI(2),PMS(0:6),IS(2),ISN(2),ROBO(5),
0023      &PSYS(0:2,5),PMIN(0:2),QOLD(4),QNEW(4),DBE(3),PSUM(4)
0024
0025 C...Find event type and remaining energy.
0026       ISUB=MINT(1)
0027       NS=N
0028       IF(MINT(50).EQ.0.OR.MOD(MSTP(81),10).LE.0) THEN
0029         VINT(143)=1D0-VINT(141)
0030         VINT(144)=1D0-VINT(142)
0031       ENDIF
0032
0033 C...Define initial partons.
0034       NTRY=0
0035   100 NTRY=NTRY+1
0036       DO 130 JT=1,2
0037         I=MINT(83)+JT+2
0038         IF(JT.EQ.1) IPU=IPU1
0039         IF(JT.EQ.2) IPU=IPU2
0040         K(I,1)=21
0041         K(I,2)=K(IPU,2)
0042         K(I,3)=I-2
0043         PMS(JT)=0D0
0044         VINT(156+JT)=0D0
0045         VINT(158+JT)=0D0
0046         IF(MINT(47).EQ.1) THEN
0047           DO 110 J=1,5
0048             P(I,J)=P(I-2,J)
0049   110     CONTINUE
0050         ELSEIF(ISUB.EQ.95) THEN
0051           K(I,2)=21
0052         ELSE
0053           P(I,5)=P(IPU,5)
0054
0055 C...No primordial kT, or chosen according to truncated Gaussian or
0056 C...exponential, or (for photon) predetermined or power law.
0057   120     IF(MINT(40+JT).EQ.2.AND.MINT(10+JT).NE.22) THEN
0058             IF(MSTP(91).LE.0) THEN
0059               PT=0D0
0060             ELSEIF(MSTP(91).EQ.1) THEN
0061               PT=PARP(91)*SQRT(-LOG(PYR(0)))
0062             ELSE
0063               RPT1=PYR(0)
0064               RPT2=PYR(0)
0065               PT=-PARP(92)*LOG(RPT1*RPT2)
0066             ENDIF
0067             IF(PT.GT.PARP(93)) GOTO 120
0068           ELSEIF(MINT(106+JT).EQ.3) THEN
0069             PTA=SQRT(VINT(282+JT))
0070             PTB=0D0
0071             IF(MSTP(66).EQ.5.AND.MSTP(93).EQ.1) THEN
0072               PTB=PARP(99)*SQRT(-LOG(PYR(0)))
0073             ELSEIF(MSTP(66).EQ.5.AND.MSTP(93).EQ.2) THEN
0074               RPT1=PYR(0)
0075               RPT2=PYR(0)
0076               PTB=-PARP(99)*LOG(RPT1*RPT2)
0077             ENDIF
0078             IF(PTB.GT.PARP(100)) GOTO 120
0079             PT=SQRT(PTA**2+PTB**2+2D0*PTA*PTB*COS(PARU(2)*PYR(0)))
0080             PT=PT*0.8D0**MINT(57)
0081             IF(NTRY.GT.10) PT=PT*0.8D0**(NTRY-10)
0082           ELSEIF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) THEN
0083             IF(MSTP(93).LE.0) THEN
0084               PT=0D0
0085             ELSEIF(MSTP(93).EQ.1) THEN
0086               PT=PARP(99)*SQRT(-LOG(PYR(0)))
0087             ELSEIF(MSTP(93).EQ.2) THEN
0088               RPT1=PYR(0)
0089               RPT2=PYR(0)
0090               PT=-PARP(99)*LOG(RPT1*RPT2)
0091             ELSEIF(MSTP(93).EQ.3) THEN
0092               HA=PARP(99)**2
0093               HB=PARP(100)**2
0094               PT=SQRT(MAX(0D0,HA*(HA+HB)/(HA+HB-PYR(0)*HB)-HA))
0095             ELSE
0096               HA=PARP(99)**2
0097               HB=PARP(100)**2
0098               IF(MSTP(93).EQ.5) HB=MIN(VINT(48),PARP(100)**2)
0099               PT=SQRT(MAX(0D0,HA*((HA+HB)/HA)**PYR(0)-HA))
0100             ENDIF
0101             IF(PT.GT.PARP(100)) GOTO 120
0102           ELSE
0103             PT=0D0
0104           ENDIF
0105           VINT(156+JT)=PT
0106           PHI=PARU(2)*PYR(0)
0107           P(I,1)=PT*COS(PHI)
0108           P(I,2)=PT*SIN(PHI)
0109           PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0110         ENDIF
0111   130 CONTINUE
0112       IF(MINT(47).EQ.1) RETURN
0113
0114 C...Kinematics construction for initial partons.
0115       I1=MINT(83)+3
0116       I2=MINT(83)+4
0117       IF(ISUB.EQ.95) THEN
0118         SHS=0D0
0119         SHR=0D0
0120       ELSE
0121         SHS=VINT(141)*VINT(142)*VINT(2)+(P(I1,1)+P(I2,1))**2+
0122      &  (P(I1,2)+P(I2,2))**2
0123         SHR=SQRT(MAX(0D0,SHS))
0124         IF((SHS-PMS(1)-PMS(2))**2-4D0*PMS(1)*PMS(2).LE.0D0) GOTO 100
0125         P(I1,4)=0.5D0*(SHR+(PMS(1)-PMS(2))/SHR)
0126         P(I1,3)=SQRT(MAX(0D0,P(I1,4)**2-PMS(1)))
0127         P(I2,4)=SHR-P(I1,4)
0128         P(I2,3)=-P(I1,3)
0129
0130 C...Transform partons to overall CM-frame.
0131         ROBO(3)=(P(I1,1)+P(I2,1))/SHR
0132         ROBO(4)=(P(I1,2)+P(I2,2))/SHR
0133         CALL PYROBO(I1,I2,0D0,0D0,-ROBO(3),-ROBO(4),0D0)
0134         ROBO(2)=PYANGL(P(I1,1),P(I1,2))
0135         CALL PYROBO(I1,I2,0D0,-ROBO(2),0D0,0D0,0D0)
0136         ROBO(1)=PYANGL(P(I1,3),P(I1,1))
0137         CALL PYROBO(I1,I2,-ROBO(1),0D0,0D0,0D0,0D0)
0138         CALL PYROBO(I2+1,MINT(52),0D0,-ROBO(2),0D0,0D0,0D0)
0139         CALL PYROBO(I1,MINT(52),ROBO(1),ROBO(2),ROBO(3),ROBO(4),0D0)
0140         ROBO(5)=(VINT(141)-VINT(142))/(VINT(141)+VINT(142))
0141         CALL PYROBO(I1,MINT(52),0D0,0D0,0D0,0D0,ROBO(5))
0142       ENDIF
0143
0144 C...Optionally fix up x and Q2 definitions for leptoproduction.
0145       IDISXQ=0
0146       IF((MINT(43).EQ.2.OR.MINT(43).EQ.3).AND.((ISUB.EQ.10.AND.
0147      &MSTP(23).GE.1).OR.(ISUB.EQ.83.AND.MSTP(23).GE.2))) IDISXQ=1
0148       IF(IDISXQ.EQ.1) THEN
0149
0150 C...Find where incoming and outgoing leptons/partons are sitting.
0151         LESD=1
0152         IF(MINT(42).EQ.1) LESD=2
0153         LPIN=MINT(83)+3-LESD
0154         LEIN=MINT(84)+LESD
0155         LQIN=MINT(84)+3-LESD
0156         LEOUT=MINT(84)+2+LESD
0157         LQOUT=MINT(84)+5-LESD
0158         IF(K(LEIN,3).GT.LEIN) LEIN=K(LEIN,3)
0159         IF(K(LQIN,3).GT.LQIN) LQIN=K(LQIN,3)
0160         LSCMS=0
0161         DO 140 I=MINT(84)+5,N
0162           IF(K(I,2).EQ.94) THEN
0163             LSCMS=I
0164             LEOUT=I+LESD
0165             LQOUT=I+3-LESD
0166           ENDIF
0167   140   CONTINUE
0168         LQBG=IPU1
0169         IF(LESD.EQ.1) LQBG=IPU2
0170
0171 C...Calculate actual and wanted momentum transfer.
0172         XNOM=VINT(43-LESD)
0173         Q2NOM=-VINT(45)
0174         HPK=2D0*(P(LPIN,4)*P(LEIN,4)-P(LPIN,1)*P(LEIN,1)-
0175      &  P(LPIN,2)*P(LEIN,2)-P(LPIN,3)*P(LEIN,3))*
0176      &  (P(MINT(83)+LESD,4)*VINT(40+LESD)/P(LEIN,4))
0177         HPT2=MAX(0D0,Q2NOM*(1D0-Q2NOM/(XNOM*HPK)))
0178         FAC=SQRT(HPT2/(P(LEOUT,1)**2+P(LEOUT,2)**2))
0179         P(N+1,1)=FAC*P(LEOUT,1)
0180         P(N+1,2)=FAC*P(LEOUT,2)
0181         P(N+1,3)=0.25D0*((HPK-Q2NOM/XNOM)/P(LPIN,4)-
0182      &  Q2NOM/(P(MINT(83)+LESD,4)*VINT(40+LESD)))*(-1)**(LESD+1)
0183         P(N+1,4)=SQRT(P(LEOUT,5)**2+P(N+1,1)**2+P(N+1,2)**2+
0184      &  P(N+1,3)**2)
0185         DO 150 J=1,4
0186           QOLD(J)=P(LEIN,J)-P(LEOUT,J)
0187           QNEW(J)=P(LEIN,J)-P(N+1,J)
0188   150   CONTINUE
0189
0190 C...Boost outgoing electron and daughters.
0191         IF(LSCMS.EQ.0) THEN
0192           DO 160 J=1,4
0193             P(LEOUT,J)=P(N+1,J)
0194   160     CONTINUE
0195         ELSE
0196           DO 170 J=1,3
0197             P(N+2,J)=(P(N+1,J)-P(LEOUT,J))/(P(N+1,4)+P(LEOUT,4))
0198   170     CONTINUE
0199           PINV=2D0/(1D0+P(N+2,1)**2+P(N+2,2)**2+P(N+2,3)**2)
0200           DO 180 J=1,3
0201             DBE(J)=PINV*P(N+2,J)
0202   180     CONTINUE
0203           DO 200 I=LSCMS+1,N
0204             IORIG=I
0205   190       IORIG=K(IORIG,3)
0206             IF(IORIG.GT.LEOUT) GOTO 190
0207             IF(I.EQ.LEOUT.OR.IORIG.EQ.LEOUT)
0208      &      CALL PYROBO(I,I,0D0,0D0,DBE(1),DBE(2),DBE(3))
0209   200     CONTINUE
0210         ENDIF
0211
0212 C...Copy shower initiator and all outgoing partons.
0213         NCOP=N+1
0214         K(NCOP,3)=LQBG
0215         DO 210 J=1,5
0216           P(NCOP,J)=P(LQBG,J)
0217   210   CONTINUE
0218         DO 240 I=MINT(84)+1,N
0219           ICOP=0
0220           IF(K(I,1).GT.10) GOTO 240
0221           IF(I.EQ.LQBG.OR.I.EQ.LQOUT) THEN
0222             ICOP=I
0223           ELSE
0224             IORIG=I
0225   220       IORIG=K(IORIG,3)
0226             IF(IORIG.EQ.LQBG.OR.IORIG.EQ.LQOUT) THEN
0227               ICOP=IORIG
0228             ELSEIF(IORIG.GT.MINT(84).AND.IORIG.LE.N) THEN
0229               GOTO 220
0230             ENDIF
0231           ENDIF
0232           IF(ICOP.NE.0) THEN
0233             NCOP=NCOP+1
0234             K(NCOP,3)=I
0235             DO 230 J=1,5
0236               P(NCOP,J)=P(I,J)
0237   230       CONTINUE
0238           ENDIF
0239   240   CONTINUE
0240
0241 C...Calculate relative rescaling factors.
0242         SLC=3-2*LESD
0243         PLCSUM=0D0
0244         DO 250 I=N+2,NCOP
0245           PLCSUM=PLCSUM+(P(I,4)+SLC*P(I,3))
0246   250   CONTINUE
0247         DO 260 I=N+2,NCOP
0248           V(I,1)=(P(I,4)+SLC*P(I,3))/PLCSUM
0249   260   CONTINUE
0250
0251 C...Transfer extra three-momentum of current.
0252         DO 280 I=N+2,NCOP
0253           DO 270 J=1,3
0254             P(I,J)=P(I,J)+V(I,1)*(QNEW(J)-QOLD(J))
0255   270     CONTINUE
0256           P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0257   280   CONTINUE
0258
0259 C...Iterate change of initiator momentum to get energy right.
0260         ITER=0
0261   290   ITER=ITER+1
0262         PEEX=-P(N+1,4)-QNEW(4)
0263         PEMV=-P(N+1,3)/P(N+1,4)
0264         DO 300 I=N+2,NCOP
0265           PEEX=PEEX+P(I,4)
0266           PEMV=PEMV+V(I,1)*P(I,3)/P(I,4)
0267   300   CONTINUE
0268         IF(ABS(PEMV).LT.1D-10) THEN
0269           MINT(51)=1
0270           MINT(57)=MINT(57)+1
0271           RETURN
0272         ENDIF
0273         PZCH=-PEEX/PEMV
0274         P(N+1,3)=P(N+1,3)+PZCH
0275         P(N+1,4)=SQRT(P(N+1,5)**2+P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)
0276         DO 310 I=N+2,NCOP
0277           P(I,3)=P(I,3)+V(I,1)*PZCH
0278           P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0279   310   CONTINUE
0280         IF(ITER.LT.10.AND.ABS(PEEX).GT.1D-6*P(N+1,4)) GOTO 290
0281
0282 C...Modify momenta in event record.
0283         HBE=2D0*(P(N+1,4)+P(LQBG,4))*(P(N+1,3)-P(LQBG,3))/
0284      &  ((P(N+1,4)+P(LQBG,4))**2+(P(N+1,3)-P(LQBG,3))**2)
0285         IF(ABS(HBE).GE.1D0) THEN
0286           MINT(51)=1
0287           MINT(57)=MINT(57)+1
0288           RETURN
0289         ENDIF
0290         I=MINT(83)+5-LESD
0291         CALL PYROBO(I,I,0D0,0D0,0D0,0D0,HBE)
0292         DO 330 I=N+1,NCOP
0293           ICOP=K(I,3)
0294           DO 320 J=1,4
0295             P(ICOP,J)=P(I,J)
0296   320     CONTINUE
0297   330   CONTINUE
0298       ENDIF
0299
0300 C...Check minimum invariant mass of remnant system(s).
0301       PSYS(0,4)=P(I1,4)+P(I2,4)+0.5D0*VINT(1)*(VINT(151)+VINT(152))
0302       PSYS(0,3)=P(I1,3)+P(I2,3)+0.5D0*VINT(1)*(VINT(151)-VINT(152))
0303       PMS(0)=MAX(0D0,PSYS(0,4)**2-PSYS(0,3)**2)
0304       PMIN(0)=SQRT(PMS(0))
0305       DO 340 JT=1,2
0306         PSYS(JT,4)=0.5D0*VINT(1)*VINT(142+JT)
0307         PSYS(JT,3)=PSYS(JT,4)*(-1)**(JT-1)
0308         PMIN(JT)=0D0
0309         IF(MINT(44+JT).EQ.1) GOTO 340
0310         MINT(105)=MINT(102+JT)
0311         MINT(109)=MINT(106+JT)
0312         CALL PYSPLI(MINT(10+JT),MINT(12+JT),KFLCH(JT),KFLSP(JT))
0313         IF(MINT(51).NE.0) THEN
0314           MINT(57)=MINT(57)+1
0315           RETURN
0316         ENDIF
0317         IF(KFLCH(JT).NE.0) PMIN(JT)=PMIN(JT)+PYMASS(KFLCH(JT))
0318         IF(KFLSP(JT).NE.0) PMIN(JT)=PMIN(JT)+PYMASS(KFLSP(JT))
0319         IF(KFLCH(JT)*KFLSP(JT).NE.0) PMIN(JT)=PMIN(JT)+0.5D0*PARP(111)
0320         PMIN(JT)=SQRT(PMIN(JT)**2+P(MINT(83)+JT+2,1)**2+
0321      &  P(MINT(83)+JT+2,2)**2)
0322   340 CONTINUE
0323       IF(PMIN(0)+PMIN(1)+PMIN(2).GT.VINT(1).OR.(MINT(45).GE.2.AND.
0324      &PMIN(1).GT.PSYS(1,4)).OR.(MINT(46).GE.2.AND.PMIN(2).GT.
0325      &PSYS(2,4))) THEN
0326         MINT(51)=1
0327         MINT(57)=MINT(57)+1
0328         RETURN
0329       ENDIF
0330
0331 C...Loop over two remnants; skip if none there.
0332       I=NS
0333       DO 410 JT=1,2
0334         ISN(JT)=0
0335         IF(MINT(44+JT).EQ.1) GOTO 410
0336         IF(JT.EQ.1) IPU=IPU1
0337         IF(JT.EQ.2) IPU=IPU2
0338
0339 C...Store first remnant parton.
0340         I=I+1
0341         IS(JT)=I
0342         ISN(JT)=1
0343         DO 350 J=1,5
0344           K(I,J)=0
0345           P(I,J)=0D0
0346           V(I,J)=0D0
0347   350   CONTINUE
0348         K(I,1)=1
0349         K(I,2)=KFLSP(JT)
0350         K(I,3)=MINT(83)+JT
0351         P(I,5)=PYMASS(K(I,2))
0352
0353 C...First parton colour connections and kinematics.
0354         KCOL=KCHG(PYCOMP(KFLSP(JT)),2)
0355         IF(KCOL.EQ.2) THEN
0356           K(I,1)=3
0357           K(I,4)=MSTU(5)*IPU+IPU
0358           K(I,5)=MSTU(5)*IPU+IPU
0359           K(IPU,4)=MOD(K(IPU,4),MSTU(5))+MSTU(5)*I
0360           K(IPU,5)=MOD(K(IPU,5),MSTU(5))+MSTU(5)*I
0361         ELSEIF(KCOL.NE.0) THEN
0362           K(I,1)=3
0363           KFLS=(3-KCOL*ISIGN(1,KFLSP(JT)))/2
0364           K(I,KFLS+3)=IPU
0365           K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
0366         ENDIF
0367         IF(KFLCH(JT).EQ.0) THEN
0368           P(I,1)=-P(MINT(83)+JT+2,1)
0369           P(I,2)=-P(MINT(83)+JT+2,2)
0370           PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0371           PSYS(JT,3)=SQRT(MAX(0D0,PSYS(JT,4)**2-PMS(JT)))*(-1)**(JT-1)
0372           P(I,3)=PSYS(JT,3)
0373           P(I,4)=PSYS(JT,4)
0374
0375 C...When extra remnant parton or hadron: store extra remnant.
0376         ELSE
0377           I=I+1
0378           ISN(JT)=2
0379           DO 360 J=1,5
0380             K(I,J)=0
0381             P(I,J)=0D0
0382             V(I,J)=0D0
0383   360     CONTINUE
0384           K(I,1)=1
0385           K(I,2)=KFLCH(JT)
0386           K(I,3)=MINT(83)+JT
0387           P(I,5)=PYMASS(K(I,2))
0388
0389 C...Find parton colour connections of extra remnant.
0390           KCOL=KCHG(PYCOMP(KFLCH(JT)),2)
0391           IF(KCOL.EQ.2) THEN
0392             K(I,1)=3
0393             K(I,4)=MSTU(5)*IPU+IPU
0394             K(I,5)=MSTU(5)*IPU+IPU
0395             K(IPU,4)=MOD(K(IPU,4),MSTU(5))+MSTU(5)*I
0396             K(IPU,5)=MOD(K(IPU,5),MSTU(5))+MSTU(5)*I
0397           ELSEIF(KCOL.NE.0) THEN
0398             K(I,1)=3
0399             KFLS=(3-KCOL*ISIGN(1,KFLCH(JT)))/2
0400             K(I,KFLS+3)=IPU
0401             K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
0402           ENDIF
0403
0404 C...Relative transverse momentum when two remnants.
0405           LOOP=0
0406   370     LOOP=LOOP+1
0407           CALL PYPTDI(1,P(I-1,1),P(I-1,2))
0408           IF(IABS(MINT(10+JT)).LT.20) THEN
0409             P(I-1,1)=0D0
0410             P(I-1,2)=0D0
0411           ELSE
0412             P(I-1,1)=P(I-1,1)-0.5D0*P(MINT(83)+JT+2,1)
0413             P(I-1,2)=P(I-1,2)-0.5D0*P(MINT(83)+JT+2,2)
0414           ENDIF
0415           PMS(JT+2)=P(I-1,5)**2+P(I-1,1)**2+P(I-1,2)**2
0416           P(I,1)=-P(MINT(83)+JT+2,1)-P(I-1,1)
0417           P(I,2)=-P(MINT(83)+JT+2,2)-P(I-1,2)
0418           PMS(JT+4)=P(I,5)**2+P(I,1)**2+P(I,2)**2
0419
0420 C...Meson or baryon; photon as meson. For splitup below.
0421           IMB=1
0422           IF(MOD(MINT(10+JT)/1000,10).NE.0) IMB=2
0423
0424 C***Relative distribution for electron into two electrons. Temporary!
0425           IF(IABS(MINT(10+JT)).LT.20.AND.MINT(14+JT).EQ.-MINT(10+JT))
0426      &    THEN
0427             CHI(JT)=PYR(0)
0428
0429 C...Relative distribution of electron energy into electron plus parton.
0430           ELSEIF(IABS(MINT(10+JT)).LT.20) THEN
0431             XHRD=VINT(140+JT)
0432             XE=VINT(154+JT)
0433             CHI(JT)=(XE-XHRD)/(1D0-XHRD)
0434
0435 C...Relative distribution of energy for particle into two jets.
0436           ELSEIF(IABS(KFLCH(JT)).LE.10.OR.KFLCH(JT).EQ.21) THEN
0437             CHIK=PARP(92+2*IMB)
0438             IF(MSTP(92).LE.1) THEN
0439               IF(IMB.EQ.1) CHI(JT)=PYR(0)
0440               IF(IMB.EQ.2) CHI(JT)=1D0-SQRT(PYR(0))
0441             ELSEIF(MSTP(92).EQ.2) THEN
0442               CHI(JT)=1D0-PYR(0)**(1D0/(1D0+CHIK))
0443             ELSEIF(MSTP(92).EQ.3) THEN
0444               CUT=2D0*0.3D0/VINT(1)
0445   380         CHI(JT)=PYR(0)**2
0446               IF((CHI(JT)**2/(CHI(JT)**2+CUT**2))**0.25D0*
0447      &        (1D0-CHI(JT))**CHIK.LT.PYR(0)) GOTO 380
0448             ELSEIF(MSTP(92).EQ.4) THEN
0449               CUT=2D0*0.3D0/VINT(1)
0450               CUTR=(1D0+SQRT(1D0+CUT**2))/CUT
0451   390         CHIR=CUT*CUTR**PYR(0)
0452               CHI(JT)=(CHIR**2-CUT**2)/(2D0*CHIR)
0453               IF((1D0-CHI(JT))**CHIK.LT.PYR(0)) GOTO 390
0454             ELSE
0455               CUT=2D0*0.3D0/VINT(1)
0456               CUTA=CUT**(1D0-PARP(98))
0457               CUTB=(1D0+CUT)**(1D0-PARP(98))
0458   400         CHI(JT)=(CUTA+PYR(0)*(CUTB-CUTA))**(1D0/(1D0-PARP(98)))
0459               IF(((CHI(JT)+CUT)**2/(2D0*(CHI(JT)**2+CUT**2)))**
0460      &        (0.5D0*PARP(98))*(1D0-CHI(JT))**CHIK.LT.PYR(0)) GOTO 400
0461             ENDIF
0462
0463 C...Relative distribution of energy for particle into jet plus particle.
0464           ELSE
0465             IF(MSTP(94).LE.1) THEN
0466               IF(IMB.EQ.1) CHI(JT)=PYR(0)
0467               IF(IMB.EQ.2) CHI(JT)=1D0-SQRT(PYR(0))
0468               IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1D0-CHI(JT)
0469             ELSEIF(MSTP(94).EQ.2) THEN
0470               CHI(JT)=1D0-PYR(0)**(1D0/(1D0+PARP(93+2*IMB)))
0471               IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1D0-CHI(JT)
0472             ELSEIF(MSTP(94).EQ.3) THEN
0473               CALL PYZDIS(1,0,PMS(JT+4),ZZ)
0474               CHI(JT)=ZZ
0475             ELSE
0476               CALL PYZDIS(1000,0,PMS(JT+4),ZZ)
0477               CHI(JT)=ZZ
0478             ENDIF
0479           ENDIF
0480
0481 C...Construct total transverse mass; reject if too large.
0482           CHI(JT)=MAX(1D-8,MIN(1D0-1D-8,CHI(JT)))
0483           PMS(JT)=PMS(JT+4)/CHI(JT)+PMS(JT+2)/(1D0-CHI(JT))
0484           IF(PMS(JT).GT.PSYS(JT,4)**2) THEN
0485             IF(LOOP.LT.100) THEN
0486               GOTO 370
0487             ELSE
0488               MINT(51)=1
0489               MINT(57)=MINT(57)+1
0490               RETURN
0491             ENDIF
0492           ENDIF
0493           PSYS(JT,3)=SQRT(MAX(0D0,PSYS(JT,4)**2-PMS(JT)))*(-1)**(JT-1)
0494           VINT(158+JT)=CHI(JT)
0495
0496 C...Subdivide longitudinal momentum according to value selected above.
0497           PW1=CHI(JT)*(PSYS(JT,4)+ABS(PSYS(JT,3)))
0498           P(IS(JT)+1,4)=0.5D0*(PW1+PMS(JT+4)/PW1)
0499           P(IS(JT)+1,3)=0.5D0*(PW1-PMS(JT+4)/PW1)*(-1)**(JT-1)
0500           P(IS(JT),4)=PSYS(JT,4)-P(IS(JT)+1,4)
0501           P(IS(JT),3)=PSYS(JT,3)-P(IS(JT)+1,3)
0502         ENDIF
0503   410 CONTINUE
0504       N=I
0505
0506 C...Check if longitudinal boosts needed - if so pick two systems.
0507       PDEV=ABS(PSYS(0,4)+PSYS(1,4)+PSYS(2,4)-VINT(1))+
0508      &ABS(PSYS(0,3)+PSYS(1,3)+PSYS(2,3))
0509       IF(PDEV.LE.1D-6*VINT(1)) RETURN
0510       IF(ISN(1).EQ.0) THEN
0511         IR=0
0512         IL=2
0513       ELSEIF(ISN(2).EQ.0) THEN
0514         IR=1
0515         IL=0
0516       ELSEIF(VINT(143).GT.0.2D0.AND.VINT(144).GT.0.2D0) THEN
0517         IR=1
0518         IL=2
0519       ELSEIF(VINT(143).GT.0.2D0) THEN
0520         IR=1
0521         IL=0
0522       ELSEIF(VINT(144).GT.0.2D0) THEN
0523         IR=0
0524         IL=2
0525       ELSEIF(PMS(1)/PSYS(1,4)**2.GT.PMS(2)/PSYS(2,4)**2) THEN
0526         IR=1
0527         IL=0
0528       ELSE
0529         IR=0
0530         IL=2
0531       ENDIF
0532       IG=3-IR-IL
0533
0534 C...E+-pL wanted for system to be modified.
0535       IF((IG.EQ.1.AND.ISN(1).EQ.0).OR.(IG.EQ.2.AND.ISN(2).EQ.0)) THEN
0536         PPB=VINT(1)
0537         PNB=VINT(1)
0538       ELSE
0539         PPB=VINT(1)-(PSYS(IG,4)+PSYS(IG,3))
0540         PNB=VINT(1)-(PSYS(IG,4)-PSYS(IG,3))
0541       ENDIF
0542
0543 C...To keep x and Q2 in leptoproduction: do not count scattered lepton.
0544       IF(IDISXQ.EQ.1.AND.IG.NE.0) THEN
0545         PPB=PPB-(PSYS(0,4)+PSYS(0,3))
0546         PNB=PNB-(PSYS(0,4)-PSYS(0,3))
0547         DO 420 J=1,4
0548           PSYS(0,J)=0D0
0549   420   CONTINUE
0550         DO 450 I=MINT(84)+1,NS
0551           IF(K(I,1).GT.10) GOTO 450
0552           INCL=0
0553           IORIG=I
0554   430     IF(IORIG.EQ.LQOUT.OR.IORIG.EQ.LPIN+2) INCL=1
0555           IORIG=K(IORIG,3)
0556           IF(IORIG.GT.LPIN) GOTO 430
0557           IF(INCL.EQ.0) GOTO 450
0558           DO 440 J=1,4
0559             PSYS(0,J)=PSYS(0,J)+P(I,J)
0560   440     CONTINUE
0561   450   CONTINUE
0562         PMS(0)=MAX(0D0,PSYS(0,4)**2-PSYS(0,3)**2)
0563         PPB=PPB+(PSYS(0,4)+PSYS(0,3))
0564         PNB=PNB+(PSYS(0,4)-PSYS(0,3))
0565       ENDIF
0566
0567 C...Construct longitudinal boosts.
0568       DPMTB=PPB*PNB
0569       DPMTR=PMS(IR)
0570       DPMTL=PMS(IL)
0571       DSQLAM=SQRT(MAX(0D0,(DPMTB-DPMTR-DPMTL)**2-4D0*DPMTR*DPMTL))
0572       IF(DSQLAM.LE.1D-6*DPMTB) THEN
0573         MINT(51)=1
0574         MINT(57)=MINT(57)+1
0575         RETURN
0576       ENDIF
0577       DSQSGN=SIGN(1D0,PSYS(IR,3)*PSYS(IL,4)-PSYS(IL,3)*PSYS(IR,4))
0578       DRKR=(DPMTB+DPMTR-DPMTL+DSQLAM*DSQSGN)/
0579      &(2D0*(PSYS(IR,4)+PSYS(IR,3))*PNB)
0580       DRKL=(DPMTB+DPMTL-DPMTR+DSQLAM*DSQSGN)/
0581      &(2D0*(PSYS(IL,4)-PSYS(IL,3))*PPB)
0582       DBER=(DRKR**2-1D0)/(DRKR**2+1D0)
0583       DBEL=-(DRKL**2-1D0)/(DRKL**2+1D0)
0584
0585 C...Perform longitudinal boosts.
0586       IF(IR.EQ.1.AND.ISN(1).EQ.1.AND.DBER.LE.-0.99999999D0) THEN
0587         P(IS(1),3)=0D0
0588         P(IS(1),4)=SQRT(P(IS(1),5)**2+P(IS(1),1)**2+P(IS(1),2)**2)
0589       ELSEIF(IR.EQ.1) THEN
0590         CALL PYROBO(IS(1),IS(1)+ISN(1)-1,0D0,0D0,0D0,0D0,DBER)
0591       ELSEIF(IDISXQ.EQ.1) THEN
0592         DO 470 I=I1,NS
0593           INCL=0
0594           IORIG=I
0595   460     IF(IORIG.EQ.LQOUT.OR.IORIG.EQ.LPIN+2) INCL=1
0596           IORIG=K(IORIG,3)
0597           IF(IORIG.GT.LPIN) GOTO 460
0598           IF(INCL.EQ.1) CALL PYROBO(I,I,0D0,0D0,0D0,0D0,DBER)
0599   470   CONTINUE
0600       ELSE
0601         CALL PYROBO(I1,NS,0D0,0D0,0D0,0D0,DBER)
0602       ENDIF
0603       IF(IL.EQ.2.AND.ISN(2).EQ.1.AND.DBEL.GE.0.99999999D0) THEN
0604         P(IS(2),3)=0D0
0605         P(IS(2),4)=SQRT(P(IS(2),5)**2+P(IS(2),1)**2+P(IS(2),2)**2)
0606       ELSEIF(IL.EQ.2) THEN
0607         CALL PYROBO(IS(2),IS(2)+ISN(2)-1,0D0,0D0,0D0,0D0,DBEL)
0608       ELSEIF(IDISXQ.EQ.1) THEN
0609         DO 490 I=I1,NS
0610           INCL=0
0611           IORIG=I
0612   480     IF(IORIG.EQ.LQOUT.OR.IORIG.EQ.LPIN+2) INCL=1
0613           IORIG=K(IORIG,3)
0614           IF(IORIG.GT.LPIN) GOTO 480
0615           IF(INCL.EQ.1) CALL PYROBO(I,I,0D0,0D0,0D0,0D0,DBEL)
0616   490   CONTINUE
0617       ELSE
0618         CALL PYROBO(I1,NS,0D0,0D0,0D0,0D0,DBEL)
0619       ENDIF
0620
0621 C...Final check that energy-momentum conservation worked.
0622       PESUM=0D0
0623       PZSUM=0D0
0624       DO 500 I=MINT(84)+1,N
0625         IF(K(I,1).GT.10) GOTO 500
0626         PESUM=PESUM+P(I,4)
0627         PZSUM=PZSUM+P(I,3)
0628   500 CONTINUE
0629       PDEV=ABS(PESUM-VINT(1))+ABS(PZSUM)
0630       IF(PDEV.GT.1D-4*VINT(1)) THEN
0631         MINT(51)=1
0632         MINT(57)=MINT(57)+1
0633         RETURN
0634       ENDIF
0635
0636 C...Calculate rotation and boost from overall CM frame to
0637 C...hadronic CM frame in leptoproduction.
0638       MINT(91)=0
0639       IF(MINT(82).EQ.1.AND.(MINT(43).EQ.2.OR.MINT(43).EQ.3)) THEN
0640         MINT(91)=1
0641         LESD=1
0642         IF(MINT(42).EQ.1) LESD=2
0643         LPIN=MINT(83)+3-LESD
0644
0645 C...Sum upp momenta of everything not lepton or photon to define boost.
0646         DO 510 J=1,4
0647           PSUM(J)=0D0
0648   510   CONTINUE
0649         DO 530 I=1,N
0650           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 530
0651           IF(IABS(K(I,2)).GE.11.AND.IABS(K(I,2)).LE.20) GOTO 530
0652           IF(K(I,2).EQ.22) GOTO 530
0653           DO 520 J=1,4
0654             PSUM(J)=PSUM(J)+P(I,J)
0655   520     CONTINUE
0656   530   CONTINUE
0657         VINT(223)=-PSUM(1)/PSUM(4)
0658         VINT(224)=-PSUM(2)/PSUM(4)
0659         VINT(225)=-PSUM(3)/PSUM(4)
0660
0661 C...Boost incoming hadron to hadronic CM frame to determine rotations.
0662         K(N+1,1)=1
0663         DO 540 J=1,5
0664           P(N+1,J)=P(LPIN,J)
0665           V(N+1,J)=V(LPIN,J)
0666   540   CONTINUE
0667         CALL PYROBO(N+1,N+1,0D0,0D0,VINT(223),VINT(224),VINT(225))
0668         VINT(222)=-PYANGL(P(N+1,1),P(N+1,2))
0669         CALL PYROBO(N+1,N+1,0D0,VINT(222),0D0,0D0,0D0)
0670         IF(LESD.EQ.2) THEN
0671           VINT(221)=-PYANGL(P(N+1,3),P(N+1,1))
0672         ELSE
0673           VINT(221)=PYANGL(-P(N+1,3),P(N+1,1))
0674         ENDIF
0675       ENDIF
0676
0677       RETURN
0678       END