PYTHIA/pythia/pyremn.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYREMN(IPU1,IPU2)
   5
   6 C...Adds on target remnants (one or two from each side) and
   7 C...includes primordial kT for hadron beams.
   8       IMPLICIT DOUBLE PRECISION(D)
   9       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  13       COMMON/PYINT1/MINT(400),VINT(400)
  14       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
  15       SAVE /PYPARS/,/PYINT1/
  16       DIMENSION KFLCH(2),KFLSP(2),CHI(2),PMS(0:6),IS(2),ISN(2),ROBO(5),
  17      &PSYS(0:2,5),PMIN(0:2),QOLD(4),QNEW(4),DBE(3),PSUM(4)
  18
  19 C...Find event type and remaining energy.
  20       ISUB=MINT(1)
  21       NS=N
  22       IF(MINT(50).EQ.0.OR.MSTP(81).LE.0) THEN
  23         VINT(143)=1.-VINT(141)
  24         VINT(144)=1.-VINT(142)
  25       ENDIF
  26
  27 C...Define initial partons.
  28       NTRY=0
  29   100 NTRY=NTRY+1
  30       DO 130 JT=1,2
  31       I=MINT(83)+JT+2
  32       IF(JT.EQ.1) IPU=IPU1
  33       IF(JT.EQ.2) IPU=IPU2
  34       K(I,1)=21
  35       K(I,2)=K(IPU,2)
  36       K(I,3)=I-2
  37       PMS(JT)=0.
  38       VINT(156+JT)=0.
  39       VINT(158+JT)=0.
  40       IF(MINT(47).EQ.1) THEN
  41         DO 110 J=1,5
  42         P(I,J)=P(I-2,J)
  43   110   CONTINUE
  44       ELSEIF(ISUB.EQ.95) THEN
  45         K(I,2)=21
  46       ELSE
  47         P(I,5)=P(IPU,5)
  48
  49 C...No primordial kT, or chosen according to truncated Gaussian or
  50 C...exponential, or (for photon) predetermined or power law.
  51   120   IF(MINT(40+JT).EQ.2.AND.MINT(10+JT).NE.22) THEN
  52           IF(MSTP(91).LE.0) THEN
  53             PT=0.
  54           ELSEIF(MSTP(91).EQ.1) THEN
  55             PT=PARP(91)*SQRT(-LOG(RLU(0)))
  56           ELSE
  57             RPT1=RLU(0)
  58             RPT2=RLU(0)
  59             PT=-PARP(92)*LOG(RPT1*RPT2)
  60           ENDIF
  61           IF(PT.GT.PARP(93)) GOTO 120
  62         ELSEIF(MINT(106+JT).EQ.3) THEN
  63           PT=SQRT(VINT(282+JT))
  64           PT=PT*0.8**MINT(57)
  65           IF(NTRY.GT.10) PT=PT*0.8**(NTRY-10)
  66         ELSEIF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) THEN
  67           IF(MSTP(93).LE.0) THEN
  68             PT=0.
  69           ELSEIF(MSTP(93).EQ.1) THEN
  70             PT=PARP(99)*SQRT(-LOG(RLU(0)))
  71           ELSEIF(MSTP(93).EQ.2) THEN
  72             RPT1=RLU(0)
  73             RPT2=RLU(0)
  74             PT=-PARP(99)*LOG(RPT1*RPT2)
  75           ELSEIF(MSTP(93).EQ.3) THEN
  76             HA=PARP(99)**2
  77             HB=PARP(100)**2
  78             PT=SQRT(MAX(0.,HA*(HA+HB)/(HA+HB-RLU(0)*HB)-HA))
  79           ELSE
  80             HA=PARP(99)**2
  81             HB=PARP(100)**2
  82             IF(MSTP(93).EQ.5) HB=MIN(VINT(48),PARP(100)**2)
  83             PT=SQRT(MAX(0.,HA*((HA+HB)/HA)**RLU(0)-HA))
  84           ENDIF
  85           IF(PT.GT.PARP(100)) GOTO 120
  86         ELSE
  87           PT=0.
  88         ENDIF
  89         VINT(156+JT)=PT
  90         PHI=PARU(2)*RLU(0)
  91         P(I,1)=PT*COS(PHI)
  92         P(I,2)=PT*SIN(PHI)
  93         PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
  94       ENDIF
  95   130 CONTINUE
  96       IF(MINT(47).EQ.1) RETURN
  97
  98 C...Kinematics construction for initial partons.
  99       I1=MINT(83)+3
 100       I2=MINT(83)+4
 101       IF(ISUB.EQ.95) THEN
 102         SHS=0.
 103         SHR=0.
 104       ELSE
 105         SHS=VINT(141)*VINT(142)*VINT(2)+(P(I1,1)+P(I2,1))**2+
 106      &  (P(I1,2)+P(I2,2))**2
 107         SHR=SQRT(MAX(0.,SHS))
 108         IF((SHS-PMS(1)-PMS(2))**2-4.*PMS(1)*PMS(2).LE.0.) GOTO 100
 109         P(I1,4)=0.5*(SHR+(PMS(1)-PMS(2))/SHR)
 110         P(I1,3)=SQRT(MAX(0.,P(I1,4)**2-PMS(1)))
 111         P(I2,4)=SHR-P(I1,4)
 112         P(I2,3)=-P(I1,3)
 113
 114 C...Transform partons to overall CM-frame.
 115         ROBO(3)=(P(I1,1)+P(I2,1))/SHR
 116         ROBO(4)=(P(I1,2)+P(I2,2))/SHR
 117         CALL LUDBRB(I1,I2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),0D0)
 118         ROBO(2)=ULANGL(P(I1,1),P(I1,2))
 119         CALL LUDBRB(I1,I2,0.,-ROBO(2),0D0,0D0,0D0)
 120         ROBO(1)=ULANGL(P(I1,3),P(I1,1))
 121         CALL LUDBRB(I1,I2,-ROBO(1),0.,0D0,0D0,0D0)
 122         CALL LUDBRB(I1,MINT(52),ROBO(1),ROBO(2),DBLE(ROBO(3)),
 123      &  DBLE(ROBO(4)),0D0)
 124         ROBO(5)=MAX(-0.999999,MIN(0.999999,(VINT(141)-VINT(142))/
 125      &  (VINT(141)+VINT(142))))
 126         CALL LUDBRB(I1,MINT(52),0.,0.,0D0,0D0,DBLE(ROBO(5)))
 127       ENDIF
 128
 129 C...Optionally fix up x and Q2 definitions for leptoproduction.
 130       IDISXQ=0
 131       IF((MINT(43).EQ.2.OR.MINT(43).EQ.3).AND.((ISUB.EQ.10.AND.
 132      &MSTP(23).GE.1).OR.(ISUB.EQ.83.AND.MSTP(23).GE.2))) IDISXQ=1
 133       IF(IDISXQ.EQ.1) THEN
 134
 135 C...Find where incoming and outgoing leptons/partons are sitting.
 136         LESD=1
 137         IF(MINT(42).EQ.1) LESD=2
 138         LPIN=MINT(83)+3-LESD
 139         LEIN=MINT(84)+LESD
 140         LQIN=MINT(84)+3-LESD
 141         LEOUT=MINT(84)+2+LESD
 142         LQOUT=MINT(84)+5-LESD
 143         IF(K(LEIN,3).GT.LEIN) LEIN=K(LEIN,3)
 144         IF(K(LQIN,3).GT.LQIN) LQIN=K(LQIN,3)
 145         LSCMS=0
 146         DO 140 I=MINT(84)+5,N
 147         IF(K(I,2).EQ.94) THEN
 148           LSCMS=I
 149           LEOUT=I+LESD
 150           LQOUT=I+3-LESD
 151         ENDIF
 152   140   CONTINUE
 153         LQBG=IPU1
 154         IF(LESD.EQ.1) LQBG=IPU2
 155
 156 C...Calculate actual and wanted momentum transfer.
 157         XNOM=VINT(43-LESD)
 158         Q2NOM=-VINT(45)
 159         HPK=2.*(P(LPIN,4)*P(LEIN,4)-P(LPIN,1)*P(LEIN,1)-
 160      &  P(LPIN,2)*P(LEIN,2)-P(LPIN,3)*P(LEIN,3))*
 161      &  (P(MINT(83)+LESD,4)*VINT(40+LESD)/P(LEIN,4))
 162         HPT2=MAX(0.,Q2NOM*(1.-Q2NOM/(XNOM*HPK)))
 163         FAC=SQRT(HPT2/(P(LEOUT,1)**2+P(LEOUT,2)**2))
 164         P(N+1,1)=FAC*P(LEOUT,1)
 165         P(N+1,2)=FAC*P(LEOUT,2)
 166         P(N+1,3)=0.25*((HPK-Q2NOM/XNOM)/P(LPIN,4)-
 167      &  Q2NOM/(P(MINT(83)+LESD,4)*VINT(40+LESD)))*(-1)**(LESD+1)
 168         P(N+1,4)=SQRT(P(LEOUT,5)**2+P(N+1,1)**2+P(N+1,2)**2+
 169      &  P(N+1,3)**2)
 170         DO 150 J=1,4
 171         QOLD(J)=P(LEIN,J)-P(LEOUT,J)
 172         QNEW(J)=P(LEIN,J)-P(N+1,J)
 173   150   CONTINUE
 174
 175 C...Boost outgoing electron and daughters.
 176         IF(LSCMS.EQ.0) THEN
 177           DO 160 J=1,4
 178           P(LEOUT,J)=P(N+1,J)
 179   160     CONTINUE
 180         ELSE
 181           DO 170 J=1,3
 182           P(N+2,J)=(P(N+1,J)-P(LEOUT,J))/(P(N+1,4)+P(LEOUT,4))
 183   170     CONTINUE
 184           PINV=2./(1.+P(N+2,1)**2+P(N+2,2)**2+P(N+2,3)**2)
 185           DO 180 J=1,3
 186           DBE(J)=PINV*P(N+2,J)
 187   180     CONTINUE
 188           DO 200 I=LSCMS+1,N
 189           IORIG=I
 190   190     IORIG=K(IORIG,3)
 191           IF(IORIG.GT.LEOUT) GOTO 190
 192           IF(I.EQ.LEOUT.OR.IORIG.EQ.LEOUT)
 193      &    CALL LUDBRB(I,I,0.,0.,DBE(1),DBE(2),DBE(3))
 194   200     CONTINUE
 195         ENDIF
 196
 197 C...Copy shower initiator and all outgoing partons.
 198         NCOP=N+1
 199         K(NCOP,3)=LQBG
 200         DO 210 J=1,5
 201         P(NCOP,J)=P(LQBG,J)
 202   210   CONTINUE
 203         DO 240 I=MINT(84)+1,N
 204         ICOP=0
 205         IF(K(I,1).GT.10) GOTO 240
 206         IF(I.EQ.LQBG.OR.I.EQ.LQOUT) THEN
 207           ICOP=I
 208         ELSE
 209           IORIG=I
 210   220     IORIG=K(IORIG,3)
 211           IF(IORIG.EQ.LQBG.OR.IORIG.EQ.LQOUT) THEN
 212             ICOP=IORIG
 213           ELSEIF(IORIG.GT.MINT(84).AND.IORIG.LE.N) THEN
 214             GOTO 220
 215           ENDIF
 216         ENDIF
 217         IF(ICOP.NE.0) THEN
 218           NCOP=NCOP+1
 219           K(NCOP,3)=I
 220           DO 230 J=1,5
 221           P(NCOP,J)=P(I,J)
 222   230     CONTINUE
 223         ENDIF
 224   240   CONTINUE
 225
 226 C...Calculate relative rescaling factors.
 227         SLC=3-2*LESD
 228         PLCSUM=0.
 229         DO 250 I=N+2,NCOP
 230         PLCSUM=PLCSUM+(P(I,4)+SLC*P(I,3))
 231   250   CONTINUE
 232         DO 260 I=N+2,NCOP
 233         V(I,1)=(P(I,4)+SLC*P(I,3))/PLCSUM
 234   260   CONTINUE
 235
 236 C...Transfer extra three-momentum of current.
 237         DO 280 I=N+2,NCOP
 238         DO 270 J=1,3
 239         P(I,J)=P(I,J)+V(I,1)*(QNEW(J)-QOLD(J))
 240   270   CONTINUE
 241         P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
 242   280   CONTINUE
 243
 244 C...Iterate change of initiator momentum to get energy right.
 245         ITER=0
 246   290   ITER=ITER+1
 247         PEEX=-P(N+1,4)-QNEW(4)
 248         PEMV=-P(N+1,3)/P(N+1,4)
 249         DO 300 I=N+2,NCOP
 250         PEEX=PEEX+P(I,4)
 251         PEMV=PEMV+V(I,1)*P(I,3)/P(I,4)
 252   300   CONTINUE
 253         IF(ABS(PEMV).LT.1E-10) THEN
 254           MINT(51)=1
 255           MINT(57)=MINT(57)+1
 256           RETURN
 257         ENDIF
 258         PZCH=-PEEX/PEMV
 259         P(N+1,3)=P(N+1,3)+PZCH
 260         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)
 261         DO 310 I=N+2,NCOP
 262         P(I,3)=P(I,3)+V(I,1)*PZCH
 263         P(I,4)=SQRT(P(I,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
 264   310   CONTINUE
 265         IF(ITER.LT.10.AND.ABS(PEEX).GT.1E-6*P(N+1,4)) GOTO 290
 266
 267 C...Modify momenta in event record.
 268         HBE=2.*(P(N+1,4)+P(LQBG,4))*(P(N+1,3)-P(LQBG,3))/
 269      &  ((P(N+1,4)+P(LQBG,4))**2+(P(N+1,3)-P(LQBG,3))**2)
 270         IF(ABS(HBE).GT.0.999999) THEN
 271           MINT(51)=1
 272           MINT(57)=MINT(57)+1
 273           RETURN
 274         ENDIF
 275         I=MINT(83)+5-LESD
 276         CALL LUDBRB(I,I,0.,0.,0D0,0D0,DBLE(HBE))
 277         DO 330 I=N+1,NCOP
 278         ICOP=K(I,3)
 279         DO 320 J=1,4
 280         P(ICOP,J)=P(I,J)
 281   320   CONTINUE
 282   330   CONTINUE
 283       ENDIF
 284
 285 C...Check minimum invariant mass of remnant system(s).
 286       PSYS(0,4)=P(I1,4)+P(I2,4)+0.5*VINT(1)*(VINT(151)+VINT(152))
 287       PSYS(0,3)=P(I1,3)+P(I2,3)+0.5*VINT(1)*(VINT(151)-VINT(152))
 288       PMS(0)=MAX(0.,PSYS(0,4)**2-PSYS(0,3)**2)
 289       PMIN(0)=SQRT(PMS(0))
 290       DO 340 JT=1,2
 291       PSYS(JT,4)=0.5*VINT(1)*VINT(142+JT)
 292       PSYS(JT,3)=PSYS(JT,4)*(-1)**(JT-1)
 293       PMIN(JT)=0.
 294       IF(MINT(44+JT).EQ.1) GOTO 340
 295       MINT(105)=MINT(102+JT)
 296       MINT(109)=MINT(106+JT)
 297       CALL PYSPLI(MINT(10+JT),MINT(12+JT),KFLCH(JT),KFLSP(JT))
 298       IF(KFLCH(JT).NE.0) PMIN(JT)=PMIN(JT)+ULMASS(KFLCH(JT))
 299       IF(KFLSP(JT).NE.0) PMIN(JT)=PMIN(JT)+ULMASS(KFLSP(JT))
 300       IF(KFLCH(JT)*KFLSP(JT).NE.0) PMIN(JT)=PMIN(JT)+0.5*PARP(111)
 301       PMIN(JT)=SQRT(PMIN(JT)**2+P(MINT(83)+JT+2,1)**2+
 302      &P(MINT(83)+JT+2,2)**2)
 303   340 CONTINUE
 304       IF(PMIN(0)+PMIN(1)+PMIN(2).GT.VINT(1).OR.(MINT(45).GE.2.AND.
 305      &PMIN(1).GT.PSYS(1,4)).OR.(MINT(46).GE.2.AND.PMIN(2).GT.
 306      &PSYS(2,4))) THEN
 307         MINT(51)=1
 308         MINT(57)=MINT(57)+1
 309         RETURN
 310       ENDIF
 311
 312 C...Loop over two remnants; skip if none there.
 313       I=NS
 314       DO 410 JT=1,2
 315       ISN(JT)=0
 316       IF(MINT(44+JT).EQ.1) GOTO 410
 317       IF(JT.EQ.1) IPU=IPU1
 318       IF(JT.EQ.2) IPU=IPU2
 319
 320 C...Store first remnant parton.
 321       I=I+1
 322       IS(JT)=I
 323       ISN(JT)=1
 324       DO 350 J=1,5
 325       K(I,J)=0
 326       P(I,J)=0.
 327       V(I,J)=0.
 328   350 CONTINUE
 329       K(I,1)=1
 330       K(I,2)=KFLSP(JT)
 331       K(I,3)=MINT(83)+JT
 332       P(I,5)=ULMASS(K(I,2))
 333
 334 C...First parton colour connections and kinematics.
 335       KCOL=KCHG(LUCOMP(KFLSP(JT)),2)
 336       IF(KCOL.EQ.2) THEN
 337         K(I,1)=3
 338         K(I,4)=MSTU(5)*IPU+IPU
 339         K(I,5)=MSTU(5)*IPU+IPU
 340         K(IPU,4)=MOD(K(IPU,4),MSTU(5))+MSTU(5)*I
 341         K(IPU,5)=MOD(K(IPU,5),MSTU(5))+MSTU(5)*I
 342       ELSEIF(KCOL.NE.0) THEN
 343         K(I,1)=3
 344         KFLS=(3-KCOL*ISIGN(1,KFLSP(JT)))/2
 345         K(I,KFLS+3)=IPU
 346         K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
 347       ENDIF
 348       IF(KFLCH(JT).EQ.0) THEN
 349         P(I,1)=-P(MINT(83)+JT+2,1)
 350         P(I,2)=-P(MINT(83)+JT+2,2)
 351         PMS(JT)=P(I,5)**2+P(I,1)**2+P(I,2)**2
 352         PSYS(JT,3)=SQRT(MAX(0.,PSYS(JT,4)**2-PMS(JT)))*(-1)**(JT-1)
 353         P(I,3)=PSYS(JT,3)
 354         P(I,4)=PSYS(JT,4)
 355
 356 C...When extra remnant parton or hadron: store extra remnant.
 357       ELSE
 358         I=I+1
 359         ISN(JT)=2
 360         DO 360 J=1,5
 361         K(I,J)=0
 362         P(I,J)=0.
 363         V(I,J)=0.
 364   360   CONTINUE
 365         K(I,1)=1
 366         K(I,2)=KFLCH(JT)
 367         K(I,3)=MINT(83)+JT
 368         P(I,5)=ULMASS(K(I,2))
 369
 370 C...Find parton colour connections of extra remnant.
 371         KCOL=KCHG(LUCOMP(KFLCH(JT)),2)
 372         IF(KCOL.EQ.2) THEN
 373           K(I,1)=3
 374           K(I,4)=MSTU(5)*IPU+IPU
 375           K(I,5)=MSTU(5)*IPU+IPU
 376           K(IPU,4)=MOD(K(IPU,4),MSTU(5))+MSTU(5)*I
 377           K(IPU,5)=MOD(K(IPU,5),MSTU(5))+MSTU(5)*I
 378         ELSEIF(KCOL.NE.0) THEN
 379           K(I,1)=3
 380           KFLS=(3-KCOL*ISIGN(1,KFLCH(JT)))/2
 381           K(I,KFLS+3)=IPU
 382           K(IPU,6-KFLS)=MOD(K(IPU,6-KFLS),MSTU(5))+MSTU(5)*I
 383         ENDIF
 384
 385 C...Relative transverse momentum when two remnants.
 386         LOOP=0
 387   370   LOOP=LOOP+1
 388         CALL LUPTDI(1,P(I-1,1),P(I-1,2))
 389         IF(IABS(MINT(10+JT)).LT.20) THEN
 390           P(I-1,1)=0.
 391           P(I-1,2)=0.
 392         ENDIF
 393         PMS(JT+2)=P(I-1,5)**2+P(I-1,1)**2+P(I-1,2)**2
 394         P(I,1)=-P(MINT(83)+JT+2,1)-P(I-1,1)
 395         P(I,2)=-P(MINT(83)+JT+2,2)-P(I-1,2)
 396         PMS(JT+4)=P(I,5)**2+P(I,1)**2+P(I,2)**2
 397
 398 C...Meson or baryon; photon as meson. For splitup below.
 399         IMB=1
 400         IF(MOD(MINT(10+JT)/1000,10).NE.0) IMB=2
 401
 402 C***Relative distribution for electron into two electrons. Temporary!
 403         IF(IABS(MINT(10+JT)).LT.20.AND.MINT(14+JT).EQ.-MINT(10+JT))
 404      &  THEN
 405           CHI(JT)=RLU(0)
 406
 407 C...Relative distribution of electron energy into electron plus parton.
 408         ELSEIF(IABS(MINT(10+JT)).LT.20) THEN
 409           XHRD=VINT(140+JT)
 410           XE=VINT(154+JT)
 411           CHI(JT)=(XE-XHRD)/(1.-XHRD)
 412
 413 C...Relative distribution of energy for particle into two jets.
 414         ELSEIF(IABS(KFLCH(JT)).LE.10.OR.KFLCH(JT).EQ.21) THEN
 415           CHIK=PARP(92+2*IMB)
 416           IF(MSTP(92).LE.1) THEN
 417             IF(IMB.EQ.1) CHI(JT)=RLU(0)
 418             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU(0))
 419           ELSEIF(MSTP(92).EQ.2) THEN
 420             CHI(JT)=1.-RLU(0)**(1./(1.+CHIK))
 421           ELSEIF(MSTP(92).EQ.3) THEN
 422             CUT=2.*0.3/VINT(1)
 423   380       CHI(JT)=RLU(0)**2
 424             IF((CHI(JT)**2/(CHI(JT)**2+CUT**2))**0.25*(1.-CHI(JT))**CHIK
 425      &      .LT.RLU(0)) GOTO 380
 426           ELSEIF(MSTP(92).EQ.4) THEN
 427             CUT=2.*0.3/VINT(1)
 428             CUTR=(1.+SQRT(1.+CUT**2))/CUT
 429   390       CHIR=CUT*CUTR**RLU(0)
 430             CHI(JT)=(CHIR**2-CUT**2)/(2.*CHIR)
 431             IF((1.-CHI(JT))**CHIK.LT.RLU(0)) GOTO 390
 432           ELSE
 433             CUT=2.*0.3/VINT(1)
 434             CUTA=CUT**(1.-PARP(98))
 435             CUTB=(1.+CUT)**(1.-PARP(98))
 436   400       CHI(JT)=(CUTA+RLU(0)*(CUTB-CUTA))**(1./(1.-PARP(98)))
 437             IF(((CHI(JT)+CUT)**2/(2.*(CHI(JT)**2+CUT**2)))**
 438      &      (0.5*PARP(98))*(1.-CHI(JT))**CHIK.LT.RLU(0)) GOTO 400
 439           ENDIF
 440
 441 C...Relative distribution of energy for particle into jet plus particle.
 442         ELSE
 443           IF(MSTP(94).LE.1) THEN
 444             IF(IMB.EQ.1) CHI(JT)=RLU(0)
 445             IF(IMB.EQ.2) CHI(JT)=1.-SQRT(RLU(0))
 446             IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1.-CHI(JT)
 447           ELSEIF(MSTP(94).EQ.2) THEN
 448             CHI(JT)=1.-RLU(0)**(1./(1.+PARP(93+2*IMB)))
 449             IF(MOD(KFLCH(JT)/1000,10).NE.0) CHI(JT)=1.-CHI(JT)
 450           ELSEIF(MSTP(94).EQ.3) THEN
 451             CALL LUZDIS(1,0,PMS(JT+4),ZZ)
 452             CHI(JT)=ZZ
 453           ELSE
 454             CALL LUZDIS(1000,0,PMS(JT+4),ZZ)
 455             CHI(JT)=ZZ
 456           ENDIF
 457         ENDIF
 458
 459 C...Construct total transverse mass; reject if too large.
 460         PMS(JT)=PMS(JT+4)/CHI(JT)+PMS(JT+2)/(1.-CHI(JT))
 461         IF(PMS(JT).GT.PSYS(JT,4)**2) THEN
 462           IF(LOOP.LT.10) THEN
 463             GOTO 370
 464           ELSE
 465             MINT(51)=1
 466             MINT(57)=MINT(57)+1
 467             RETURN
 468           ENDIF
 469         ENDIF
 470         PSYS(JT,3)=SQRT(MAX(0.,PSYS(JT,4)**2-PMS(JT)))*(-1)**(JT-1)
 471         VINT(158+JT)=CHI(JT)
 472
 473 C...Subdivide longitudinal momentum according to value selected above.
 474         PW1=CHI(JT)*(PSYS(JT,4)+ABS(PSYS(JT,3)))
 475         P(IS(JT)+1,4)=0.5*(PW1+PMS(JT+4)/PW1)
 476         P(IS(JT)+1,3)=0.5*(PW1-PMS(JT+4)/PW1)*(-1)**(JT-1)
 477         P(IS(JT),4)=PSYS(JT,4)-P(IS(JT)+1,4)
 478         P(IS(JT),3)=PSYS(JT,3)-P(IS(JT)+1,3)
 479       ENDIF
 480   410 CONTINUE
 481       N=I
 482
 483 C...Check if longitudinal boosts needed - if so pick two systems.
 484       PDEV=ABS(PSYS(0,4)+PSYS(1,4)+PSYS(2,4)-VINT(1))+
 485      &ABS(PSYS(0,3)+PSYS(1,3)+PSYS(2,3))
 486       IF(PDEV.LE.1E-6*VINT(1)) RETURN
 487       IF(ISN(1).EQ.0) THEN
 488         IR=0
 489         IL=2
 490       ELSEIF(ISN(2).EQ.0) THEN
 491         IR=1
 492         IL=0
 493       ELSEIF(VINT(143).GT.0.2.AND.VINT(144).GT.0.2) THEN
 494         IR=1
 495         IL=2
 496       ELSEIF(VINT(143).GT.0.2) THEN
 497         IR=1
 498         IL=0
 499       ELSEIF(VINT(144).GT.0.2) THEN
 500         IR=0
 501         IL=2
 502       ELSEIF(PMS(1)/PSYS(1,4)**2.GT.PMS(2)/PSYS(2,4)**2) THEN
 503         IR=1
 504         IL=0
 505       ELSE
 506         IR=0
 507         IL=2
 508       ENDIF
 509       IG=3-IR-IL
 510
 511 C...E+-pL wanted for system to be modified.
 512       IF((IG.EQ.1.AND.ISN(1).EQ.0).OR.(IG.EQ.2.AND.ISN(2).EQ.0)) THEN
 513         PPB=VINT(1)
 514         PNB=VINT(1)
 515       ELSE
 516         PPB=VINT(1)-(PSYS(IG,4)+PSYS(IG,3))
 517         PNB=VINT(1)-(PSYS(IG,4)-PSYS(IG,3))
 518       ENDIF
 519
 520 C...To keep x and Q2 in leptoproduction: do not count scattered lepton.
 521       IF(IDISXQ.EQ.1.AND.IG.NE.0) THEN
 522         PMTB=PPB*PNB
 523         PMTR=PMS(IR)
 524         PMTL=PMS(IL)
 525         SQLAM=SQRT(MAX(0.,(PMTB-PMTR-PMTL)**2-4.*PMTR*PMTL))
 526         SQSGN=SIGN(1.,PSYS(IR,3)*PSYS(IL,4)-PSYS(IL,3)*PSYS(IR,4))
 527         RKR=(PMTB+PMTR-PMTL+SQLAM*SQSGN)/(2.*(PSYS(IR,4)+PSYS(IR,3))
 528      &  *PNB)
 529         RKL=(PMTB+PMTL-PMTR+SQLAM*SQSGN)/(2.*(PSYS(IL,4)-PSYS(IL,3))
 530      &  *PPB)
 531         BER=(RKR**2-1.)/(RKR**2+1.)
 532         BEL=-(RKL**2-1.)/(RKL**2+1.)
 533         PPB=PPB-(PSYS(0,4)+PSYS(0,3))
 534         PNB=PNB-(PSYS(0,4)-PSYS(0,3))
 535         DO 420 J=1,4
 536         PSYS(0,J)=0.
 537   420   CONTINUE
 538         DO 450 I=MINT(84)+1,NS
 539         IF(K(I,1).GT.10) GOTO 450
 540         INCL=0
 541         IORIG=I
 542   430   IF(IORIG.EQ.LQOUT.OR.IORIG.EQ.LPIN+2) INCL=1
 543         IORIG=K(IORIG,3)
 544         IF(IORIG.GT.LPIN) GOTO 430
 545         IF(INCL.EQ.0) GOTO 450
 546         DO 440 J=1,4
 547         PSYS(0,J)=PSYS(0,J)+P(I,J)
 548   440   CONTINUE
 549   450   CONTINUE
 550         PMS(0)=MAX(0.,PSYS(0,4)**2-PSYS(0,3)**2)
 551         PPB=PPB+(PSYS(0,4)+PSYS(0,3))
 552         PNB=PNB+(PSYS(0,4)-PSYS(0,3))
 553       ENDIF
 554
 555 C...Construct longitudinal boosts.
 556       DPMTB=PPB*PNB
 557       DPMTR=PMS(IR)
 558       DPMTL=PMS(IL)
 559       DSQLAM=SQRT(MAX(0D0,(DPMTB-DPMTR-DPMTL)**2-4D0*DPMTR*DPMTL))
 560       IF(DSQLAM.LE.1D-6*DPMTB) THEN
 561         MINT(51)=1
 562         MINT(57)=MINT(57)+1
 563         RETURN
 564       ENDIF
 565       DSQSGN=SIGN(1D0,DBLE(PSYS(IR,3)*PSYS(IL,4)-PSYS(IL,3)*PSYS(IR,4)))
 566       DRKR=(DPMTB+DPMTR-DPMTL+DSQLAM*DSQSGN)/
 567      &(2.*(PSYS(IR,4)+PSYS(IR,3))*PNB)
 568       DRKL=(DPMTB+DPMTL-DPMTR+DSQLAM*DSQSGN)/
 569      &(2.*(PSYS(IL,4)-PSYS(IL,3))*PPB)
 570       DBER=(DRKR**2-1D0)/(DRKR**2+1D0)
 571       DBEL=-(DRKL**2-1D0)/(DRKL**2+1D0)
 572
 573 C...Perform longitudinal boosts.
 574       IF(IR.EQ.1.AND.ISN(1).EQ.1.AND.DBER.LE.-0.99999999D0) THEN
 575         P(IS(1),3)=0.
 576         P(IS(1),4)=SQRT(P(IS(1),5)**2+P(IS(1),1)**2+P(IS(1),2)**2)
 577       ELSEIF(IR.EQ.1) THEN
 578         CALL LUDBRB(IS(1),IS(1)+ISN(1)-1,0.,0.,0D0,0D0,DBER)
 579       ELSEIF(IDISXQ.EQ.1) THEN
 580         DO 470 I=I1,NS
 581         INCL=0
 582         IORIG=I
 583   460   IF(IORIG.EQ.LQOUT.OR.IORIG.EQ.LPIN+2) INCL=1
 584         IORIG=K(IORIG,3)
 585         IF(IORIG.GT.LPIN) GOTO 460
 586         IF(INCL.EQ.1) CALL LUDBRB(I,I,0.,0.,0D0,0D0,DBER)
 587   470   CONTINUE
 588       ELSE
 589         CALL LUDBRB(I1,NS,0.,0.,0D0,0D0,DBER)
 590       ENDIF
 591       IF(IL.EQ.2.AND.ISN(2).EQ.1.AND.DBEL.GE.0.99999999D0) THEN
 592         P(IS(2),3)=0.
 593         P(IS(2),4)=SQRT(P(IS(2),5)**2+P(IS(2),1)**2+P(IS(2),2)**2)
 594       ELSEIF(IL.EQ.2) THEN
 595         CALL LUDBRB(IS(2),IS(2)+ISN(2)-1,0.,0.,0D0,0D0,DBEL)
 596       ELSEIF(IDISXQ.EQ.1) THEN
 597         DO 490 I=I1,NS
 598         INCL=0
 599         IORIG=I
 600   480   IF(IORIG.EQ.LQOUT.OR.IORIG.EQ.LPIN+2) INCL=1
 601         IORIG=K(IORIG,3)
 602         IF(IORIG.GT.LPIN) GOTO 480
 603         IF(INCL.EQ.1) CALL LUDBRB(I,I,0.,0.,0D0,0D0,DBEL)
 604   490   CONTINUE
 605       ELSE
 606         CALL LUDBRB(I1,NS,0.,0.,0D0,0D0,DBEL)
 607       ENDIF
 608
 609 C...Final check that energy-momentum conservation worked.
 610       PESUM=0.
 611       PZSUM=0.
 612       DO 500 I=MINT(84)+1,N
 613       IF(K(I,1).GT.10) GOTO 500
 614       PESUM=PESUM+P(I,4)
 615       PZSUM=PZSUM+P(I,3)
 616   500 CONTINUE
 617       PDEV=ABS(PESUM-VINT(1))+ABS(PZSUM)
 618       IF(PDEV.GT.1E-4*VINT(1)) THEN
 619         MINT(51)=1
 620         MINT(57)=MINT(57)+1
 621         RETURN
 622       ENDIF
 623
 624 C...Calculate rotation and boost from overall CM frame to
 625 C...hadronic CM frame in leptoproduction.
 626       MINT(91)=0
 627       IF(MINT(82).EQ.1.AND.(MINT(43).EQ.2.OR.MINT(43).EQ.3)) THEN
 628         MINT(91)=1
 629         LESD=1
 630         IF(MINT(42).EQ.1) LESD=2
 631         LPIN=MINT(83)+3-LESD
 632
 633 C...Sum upp momenta of everything not lepton or photon to define boost.
 634         DO 510 J=1,4
 635         PSUM(J)=0.
 636   510   CONTINUE
 637         DO 530 I=1,N
 638         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 530
 639         IF(IABS(K(I,2)).GE.11.AND.IABS(K(I,2)).LE.20) GOTO 530
 640         IF(K(I,2).EQ.22) GOTO 530
 641         DO 520 J=1,4
 642         PSUM(J)=PSUM(J)+P(I,J)
 643   520   CONTINUE
 644   530   CONTINUE
 645         VINT(223)=-PSUM(1)/PSUM(4)
 646         VINT(224)=-PSUM(2)/PSUM(4)
 647         VINT(225)=-PSUM(3)/PSUM(4)
 648
 649 C...Boost incoming hadron to hadronic CM frame to determine rotations.
 650         K(N+1,1)=1
 651         DO 540 J=1,5
 652         P(N+1,J)=P(LPIN,J)
 653         V(N+1,J)=V(LPIN,J)
 654   540   CONTINUE
 655         CALL LUDBRB(N+1,N+1,0.,0.,DBLE(VINT(223)),DBLE(VINT(224)),
 656      &  DBLE(VINT(225)))
 657         VINT(222)=-ULANGL(P(N+1,1),P(N+1,2))
 658         CALL LUDBRB(N+1,N+1,0.,VINT(222),0D0,0D0,0D0)
 659         IF(LESD.EQ.2) THEN
 660           VINT(221)=-ULANGL(P(N+1,3),P(N+1,1))
 661         ELSE
 662           VINT(221)=ULANGL(-P(N+1,3),P(N+1,1))
 663         ENDIF
 664       ENDIF
 665
 666       RETURN
 667       END