| fe4da5cc |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame