File indexing completed on 2024-04-06 12:13:45
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192 SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
0193 CHARACTER FRAME*8
0194 DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
0195 & IPCOL(90000),ITCOL(90000)
0196 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0197
0198 COMMON/HIJCRDN/YP(3,300),YT(3,300)
0199 COMMON/HIJGLBR/NELT,NINT,NELP,NINP
0200
0201 COMMON/HIMAIN1/NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
0202 COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4),VATT(130000,4)
0203 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
0204 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
0205 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
0206 & PJPM(300,500),NTJ(300),KFTJ(300,500),
0207 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
0208 & PJTE(300,500),PJTM(300,500)
0209 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
0210 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
0211 & PZSG(900,100),PESG(900,100),PMSG(900,100)
0212 COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5),VDR(900,4)
0213 COMMON/RANSEED/NSEED
0214
0215 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0216 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0217 SAVE
0218
0219
0220 IERRSTAT = 0
0221 BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
0222 BMIN=MIN(BMIN0,BMAX)
0223 IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
0224 BMIN=0.0
0225 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
0226 ENDIF
0227
0228
0229
0230
0231 YP(1,1)=0.0
0232 YP(2,1)=0.0
0233 YP(3,1)=0.0
0234 IF(IHNT2(1).LE.1) GO TO 14
0235 DO 10 KP=1,IHNT2(1)
0236 5 R=HIRND(1)
0237
0238 if(IHNT2(1).EQ.2) then
0239 rnd1=max(RLU(NSEED),1.0e-20)
0240 rnd2=max(RLU(NSEED),1.0e-20)
0241 rnd3=max(RLU(NSEED),1.0e-20)
0242 R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
0243 & +4.38*0.85*log(rnd3)/(4.38+0.85))
0244 endif
0245
0246 X=RLU(NSEED)
0247 CX=2.0*X-1.0
0248 SX=SQRT(1.0-CX*CX)
0249
0250 PHI=RLU(NSEED)*2.0*HIPR1(40)
0251
0252 YP(1,KP)=R*SX*COS(PHI)
0253 YP(2,KP)=R*SX*SIN(PHI)
0254 YP(3,KP)=R*CX
0255 IF(HIPR1(29).EQ.0.0) GO TO 10
0256 DO 8 KP2=1,KP-1
0257 DNBP1=(YP(1,KP)-YP(1,KP2))**2
0258 DNBP2=(YP(2,KP)-YP(2,KP2))**2
0259 DNBP3=(YP(3,KP)-YP(3,KP2))**2
0260 DNBP=DNBP1+DNBP2+DNBP3
0261 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
0262
0263
0264 8 CONTINUE
0265 10 CONTINUE
0266
0267 if(IHNT2(1).EQ.2) then
0268 YP(1,2)=-YP(1,1)
0269 YP(2,2)=-YP(2,1)
0270 YP(3,2)=-YP(3,1)
0271 endif
0272
0273 DO 12 I=1,IHNT2(1)-1
0274 DO 12 J=I+1,IHNT2(1)
0275 IF(YP(3,I).GT.YP(3,J)) GO TO 12
0276 Y1=YP(1,I)
0277 Y2=YP(2,I)
0278 Y3=YP(3,I)
0279 YP(1,I)=YP(1,J)
0280 YP(2,I)=YP(2,J)
0281 YP(3,I)=YP(3,J)
0282 YP(1,J)=Y1
0283 YP(2,J)=Y2
0284 YP(3,J)=Y3
0285 12 CONTINUE
0286
0287
0288 14 YT(1,1)=0.0
0289 YT(2,1)=0.0
0290 YT(3,1)=0.0
0291 IF(IHNT2(3).LE.1) GO TO 24
0292 DO 20 KT=1,IHNT2(3)
0293 15 R=HIRND(2)
0294
0295 if(IHNT2(3).EQ.2) then
0296 rnd1=max(RLU(NSEED),1.0e-20)
0297 rnd2=max(RLU(NSEED),1.0e-20)
0298 rnd3=max(RLU(NSEED),1.0e-20)
0299 R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
0300 & +4.38*0.85*log(rnd3)/(4.38+0.85))
0301 endif
0302
0303 X=RLU(NSEED)
0304 CX=2.0*X-1.0
0305 SX=SQRT(1.0-CX*CX)
0306
0307 PHI=RLU(NSEED)*2.0*HIPR1(40)
0308
0309 YT(1,KT)=R*SX*COS(PHI)
0310 YT(2,KT)=R*SX*SIN(PHI)
0311 YT(3,KT)=R*CX
0312 IF(HIPR1(29).EQ.0.0) GO TO 20
0313 DO 18 KT2=1,KT-1
0314 DNBT1=(YT(1,KT)-YT(1,KT2))**2
0315 DNBT2=(YT(2,KT)-YT(2,KT2))**2
0316 DNBT3=(YT(3,KT)-YT(3,KT2))**2
0317 DNBT=DNBT1+DNBT2+DNBT3
0318 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
0319
0320
0321 18 CONTINUE
0322 20 CONTINUE
0323
0324 if(IHNT2(3).EQ.2) then
0325 YT(1,2)=-YT(1,1)
0326 YT(2,2)=-YT(2,1)
0327 YT(3,2)=-YT(3,1)
0328 endif
0329
0330 DO 22 I=1,IHNT2(3)-1
0331 DO 22 J=I+1,IHNT2(3)
0332 IF(YT(3,I).LT.YT(3,J)) GO TO 22
0333 Y1=YT(1,I)
0334 Y2=YT(2,I)
0335 Y3=YT(3,I)
0336 YT(1,I)=YT(1,J)
0337 YT(2,I)=YT(2,J)
0338 YT(3,I)=YT(3,J)
0339 YT(1,J)=Y1
0340 YT(2,J)=Y2
0341 YT(3,J)=Y3
0342 22 CONTINUE
0343
0344 24 MISS=-1
0345
0346 50 MISS=MISS+1
0347 IF(MISS.GT.50) THEN
0348 WRITE(6,*) 'infinite loop happened in HIJING'
0349 STOP
0350 ENDIF
0351
0352 NATT=0
0353 JATT=0
0354 EATT=0.0
0355 CALL HIJINI
0356 NLOP=0
0357
0358 60 NT=0
0359 NP=0
0360 N0=0
0361 N01=0
0362 N10=0
0363 N11=0
0364 NELT=0
0365 NINT=0
0366 NELP=0
0367 NINP=0
0368 NSG=0
0369 NCOLT=0
0370
0371
0372
0373
0374
0375 BB=SQRT(BMIN**2+RLU(NSEED)*(BMAX**2-BMIN**2))
0376 PHI=2.0*HIPR1(40)*RLU(NSEED)
0377 BBX=BB*COS(PHI)
0378 BBY=BB*SIN(PHI)
0379 HINT1(19)=BB
0380 HINT1(20)=PHI
0381
0382 DO 70 JP=1,IHNT2(1)
0383 DO 70 JT=1,IHNT2(3)
0384 SCIP(JP,JT)=-1.0
0385 B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
0386 R2=B2*HIPR1(40)/HIPR1(31)/0.1
0387
0388 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
0389 & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
0390 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
0391 & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
0392 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
0393 & *SQRT(1.0-RRB1)
0394 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
0395 & *SQRT(1.0-RRB2)
0396 HINT1(18)=HINT1(14)-APHX1*HINT1(15)
0397 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
0398 IF(IHPR2(14).EQ.0.OR.
0399 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
0400 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
0401 RANTOT=RLU(NSEED)
0402 IF(RANTOT.GT.GS) GO TO 70
0403 GO TO 65
0404 ENDIF
0405 GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
0406 & /HIPR1(31)/2.0*ROMG(0.0)))
0407 R2=R2/GSTOT_0
0408 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
0409 GSTOT=2.0*(1.0-SQRT(1.0-GS))
0410 RANTOT=RLU(NSEED)*GSTOT_0
0411 IF(RANTOT.GT.GSTOT) GO TO 70
0412 IF(RANTOT.GT.GS) THEN
0413 CALL HIJCSC(JP,JT)
0414 GO TO 70
0415
0416 ENDIF
0417 65 SCIP(JP,JT)=R2
0418 RNIP(JP,JT)=RANTOT
0419 SJIP(JP,JT)=HINT1(18)
0420 NCOLT=NCOLT+1
0421 IPCOL(NCOLT)=JP
0422 ITCOL(NCOLT)=JT
0423 70 CONTINUE
0424
0425
0426 IF(NCOLT.EQ.0) THEN
0427 NLOP=NLOP+1
0428 IF(NLOP.LE.20.OR.
0429 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
0430 RETURN
0431 ENDIF
0432
0433
0434
0435
0436 IF(IHPR2(3).NE.0) THEN
0437 NHARD=1+INT(RLU(NSEED)*(NCOLT-1)+0.5)
0438 NHARD=MIN(NHARD,NCOLT)
0439 JPHARD=IPCOL(NHARD)
0440 JTHARD=ITCOL(NHARD)
0441 ENDIF
0442
0443 IF(IHPR2(9).EQ.1) THEN
0444 NMINI=1+INT(RLU(NSEED)*(NCOLT-1)+0.5)
0445 NMINI=MIN(NMINI,NCOLT)
0446 JPMINI=IPCOL(NMINI)
0447 JTMINI=ITCOL(NMINI)
0448 ENDIF
0449
0450
0451
0452 DO 200 JP=1,IHNT2(1)
0453 DO 200 JT=1,IHNT2(3)
0454 IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
0455 NFP(JP,11)=NFP(JP,11)+1
0456 NFT(JT,11)=NFT(JT,11)+1
0457 IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
0458 NP=NP+1
0459 N01=N01+1
0460 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
0461 NT=NT+1
0462 N10=N10+1
0463 ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
0464 NP=NP+1
0465 NT=NT+1
0466 N0=N0+1
0467 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
0468 N11=N11+1
0469 ENDIF
0470
0471 JOUT=0
0472 NFP(JP,10)=0
0473 NFT(JT,10)=0
0474
0475 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
0476
0477 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
0478
0479
0480
0481 R2=SCIP(JP,JT)
0482 HINT1(18)=SJIP(JP,JT)
0483 TT=ROMG(R2)*HINT1(18)/HIPR1(31)
0484 TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
0485 NJET=0
0486 IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
0487 CALL JETINI(JP,JT,1)
0488 CALL HIJHRD(JP,JT,0,JFLG,0)
0489 HINT1(26)=HINT1(47)
0490 HINT1(27)=HINT1(48)
0491 HINT1(28)=HINT1(49)
0492 HINT1(29)=HINT1(50)
0493 HINT1(36)=HINT1(67)
0494 HINT1(37)=HINT1(68)
0495 HINT1(38)=HINT1(69)
0496 HINT1(39)=HINT1(70)
0497
0498 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
0499 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
0500 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
0501 & JFLG.GE.3) IASG(NSG,3)=1
0502 IHNT2(9)=IHNT2(14)
0503 IHNT2(10)=IHNT2(15)
0504 DO 105 I05=1,5
0505 HINT1(20+I05)=HINT1(40+I05)
0506 HINT1(30+I05)=HINT1(50+I05)
0507 105 CONTINUE
0508 JOUT=1
0509 IF(IHPR2(8).EQ.0) GO TO 160
0510 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
0511 & /REAL(IHNT2(1))**0.6666667,1.0)
0512 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
0513 & /REAL(IHNT2(3))**0.6666667,1.0)
0514 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
0515 & *SQRT(1.0-RRB1)
0516 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
0517 & *SQRT(1.0-RRB2)
0518 HINT1(65)=HINT1(61)-APHX1*HINT1(62)
0519 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
0520 TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
0521 NJET=-1
0522
0523
0524
0525 XR1=-ALOG(EXP(-TTRIG)+RLU(NSEED)*(1.0-EXP(-TTRIG)))
0526 106 NJET=NJET+1
0527 XR1=XR1-ALOG(max(RLU(NSEED),1.0e-20))
0528 IF(XR1.LT.TTRIG) GO TO 106
0529 XR=0.0
0530 107 NJET=NJET+1
0531 XR=XR-ALOG(max(RLU(NSEED),1.0e-20))
0532 IF(XR.LT.TT-TTRIG) GO TO 107
0533 NJET=NJET-1
0534 GO TO 112
0535 ENDIF
0536
0537
0538 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
0539
0540
0541
0542 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
0543 & (1.0-EXP(-TTS))) GO TO 160
0544
0545 110 XTMP=EXP(-TT)
0546 XR=-ALOG(XTMP+HIJRAN(NSEED)*(1.0-XTMP))
0547 111 NJET=NJET+1
0548 XR=XR-ALOG(max(RLU(NSEED),1.0e-20))
0549 IF(XR.LT.TT) GO TO 111
0550 112 NJET=MIN(NJET,IHPR2(8))
0551 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
0552
0553
0554 DO 150 I_JET=1,NJET
0555 CALL JETINI(JP,JT,0)
0556 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
0557
0558
0559
0560
0561
0562
0563 IF(JFLG.EQ.0) GO TO 160
0564 IF(JFLG.LT.0) THEN
0565 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
0566 GO TO 50
0567 ENDIF
0568 JOUT=JOUT+1
0569 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
0570 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
0571 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
0572 & JFLG.GE.3) IASG(NSG,3)=1
0573
0574 150 CONTINUE
0575 160 CONTINUE
0576 CALL HIJSFT(JP,JT,JOUT,IERROR)
0577 IF(IERROR.NE.0) THEN
0578 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
0579 GO TO 50
0580 ENDIF
0581
0582
0583 JATT=JATT+JOUT
0584
0585 200 CONTINUE
0586
0587
0588
0589 DO 201 JP=1,IHNT2(1)
0590 IF(NFP(JP,5).GT.2) THEN
0591 NINP=NINP+1
0592 ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
0593 NELP=NELP+1
0594 ENDIF
0595 201 continue
0596 DO 202 JT=1,IHNT2(3)
0597 IF(NFT(JT,5).GT.2) THEN
0598 NINT=NINT+1
0599 ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
0600 NELT=NELT+1
0601 ENDIF
0602 202 continue
0603
0604
0605
0606
0607
0608
0609 IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
0610 & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
0611 DO 271 I=1,IHNT2(1)
0612 IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
0613 271 CONTINUE
0614 DO 272 I=1,IHNT2(3)
0615 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
0616 272 CONTINUE
0617 DO 273 ISG=1,NSG
0618 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
0619 273 CONTINUE
0620 ENDIF
0621
0622
0623
0624
0625
0626
0627
0628
0629 IF(IHPR2(20).NE.0) THEN
0630 DO 360 ISG=1,NSG
0631 CALL HIJFRG(ISG,3,IERROR)
0632 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
0633 MSTU(24)=0
0634 MSTU(28)=0
0635 IF(IHPR2(10).NE.0) THEN
0636 call lulist(1)
0637 WRITE(6,*) 'error occured, repeat the event'
0638 ENDIF
0639 GO TO 50
0640 ENDIF
0641
0642
0643 N_ST=1
0644 IDSTR=92
0645 IF(IHPR2(21).EQ.0) THEN
0646 CALL LUEDIT(2)
0647
0648
0649 ELSE IF(N.GT.1) THEN
0650 351 N_ST=N_ST+1
0651
0652 IF(N_ST.GT.N) THEN
0653 IERRSTAT=2
0654 RETURN
0655 ENDIF
0656
0657 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351
0658 IDSTR=K(N_ST,2)
0659 N_ST=N_ST+1
0660 ENDIF
0661
0662 IF(FRAME.EQ.'LAB') THEN
0663 CALL HIBOOST
0664 ENDIF
0665
0666
0667 N_STR=0
0668 DO 360 I=N_ST,N
0669 IF(K(I,2).EQ.IDSTR) THEN
0670 IF(K(I,3).LT.N_ST) THEN
0671 N_STR=N_STR+1
0672 GO TO 360
0673 ENDIF
0674 ENDIF
0675 K(I,4)=N_STR
0676 NATT=NATT+1
0677
0678 IF(NATT.GT.130000) THEN
0679 IERRSTAT=1
0680 RETURN
0681 ENDIF
0682 KATT(NATT,1)=K(I,2)
0683 KATT(NATT,2)=20
0684 KATT(NATT,4)=K(I,1)
0685 IF(K(I,3).EQ.0 .OR. K(I,3).LT.N_ST .OR.
0686 & (K(K(I,3),2).EQ.IDSTR .AND.
0687 & K(K(I,3),3).LT.N_ST)) THEN
0688 KATT(NATT,3)=0
0689 ELSE
0690 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
0691 ENDIF
0692
0693 PATT(NATT,1)=P(I,1)
0694 PATT(NATT,2)=P(I,2)
0695 PATT(NATT,3)=P(I,3)
0696 PATT(NATT,4)=P(I,4)
0697 EATT=EATT+P(I,4)
0698 VATT(NATT,1)=V(I,1)
0699 VATT(NATT,2)=V(I,2)
0700 VATT(NATT,3)=V(I,3)
0701 VATT(NATT,4)=V(I,4)
0702 360 CONTINUE
0703
0704
0705 JTP(1)=IHNT2(1)
0706 JTP(2)=IHNT2(3)
0707 DO 400 NTP=1,2
0708 DO 400 J_JTP=1,JTP(NTP)
0709 CALL HIJFRG(J_JTP,NTP,IERROR)
0710 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
0711 MSTU(24)=0
0712 MSTU(28)=0
0713 IF(IHPR2(10).NE.0) THEN
0714 call lulist(1)
0715 WRITE(6,*) 'error occured, repeat the event'
0716 ENDIF
0717 GO TO 50
0718 ENDIF
0719
0720
0721 N_ST=1
0722 IDSTR=92
0723 IF(IHPR2(21).EQ.0) THEN
0724 CALL LUEDIT(2)
0725
0726
0727 ELSE IF(N.GT.1) THEN
0728 381 N_ST=N_ST+1
0729
0730 IF(N_ST.GT.N) THEN
0731 IERRSTAT=2
0732 RETURN
0733 ENDIF
0734 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
0735 IDSTR=K(N_ST,2)
0736 N_ST=N_ST+1
0737 ENDIF
0738 IF(FRAME.EQ.'LAB') THEN
0739 CALL HIBOOST
0740 ENDIF
0741
0742
0743 NFTP=NFP(J_JTP,5)
0744 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
0745 N_STR=0
0746 DO 390 I=N_ST,N
0747 IF(K(I,2).EQ.IDSTR) THEN
0748 IF(K(I,3).LT.N_ST) THEN
0749 N_STR=N_STR+1
0750 GO TO 390
0751 ENDIF
0752 ENDIF
0753 K(I,4)=N_STR
0754 NATT=NATT+1
0755
0756 IF(NATT.GT.130000) THEN
0757 IERRSTAT=1
0758 RETURN
0759 ENDIF
0760 KATT(NATT,1)=K(I,2)
0761 KATT(NATT,2)=NFTP
0762 KATT(NATT,4)=K(I,1)
0763 IF(K(I,3).EQ.0 .OR. K(I,3).LT.N_ST .OR.
0764 & (K(K(I,3),2).EQ.IDSTR .AND.
0765 & K(K(I,3),3).LT.N_ST)) THEN
0766 KATT(NATT,3)=0
0767 ELSE
0768 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
0769 ENDIF
0770
0771 PATT(NATT,1)=P(I,1)
0772 PATT(NATT,2)=P(I,2)
0773 PATT(NATT,3)=P(I,3)
0774 PATT(NATT,4)=P(I,4)
0775 EATT=EATT+P(I,4)
0776 VATT(NATT,1)=V(I,1)
0777 VATT(NATT,2)=V(I,2)
0778 VATT(NATT,3)=V(I,3)
0779 VATT(NATT,4)=V(I,4)
0780 390 CONTINUE
0781 400 CONTINUE
0782
0783 ENDIF
0784
0785 DO 450 I=1,NDR
0786 NATT=NATT+1
0787
0788 IF(NATT.GT.130000) THEN
0789 IERRSTAT=1
0790 RETURN
0791 ENDIF
0792 KATT(NATT,1)=KFDR(I)
0793 KATT(NATT,2)=40
0794 KATT(NATT,3)=0
0795 PATT(NATT,1)=PDR(I,1)
0796 PATT(NATT,2)=PDR(I,2)
0797 PATT(NATT,3)=PDR(I,3)
0798 PATT(NATT,4)=PDR(I,4)
0799 EATT=EATT+PDR(I,4)
0800 VATT(NATT,1)=VDR(I,1)
0801 VATT(NATT,2)=VDR(I,2)
0802 VATT(NATT,3)=VDR(I,3)
0803 VATT(NATT,4)=VDR(I,4)
0804 450 CONTINUE
0805
0806
0807 DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
0808 IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
0809 & .AND.IHPR2(21).EQ.0) THEN
0810 IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, '//
0811 & 'repeat the event'
0812
0813 GO TO 50
0814 ENDIF
0815 RETURN
0816 END
0817
0818
0819
0820 SUBROUTINE HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
0821 CHARACTER FRAME*4,PROJ*4,TARG*4,EFRAME*4
0822 DOUBLE PRECISION DD1,DD2,DD3,DD4
0823 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
0824 COMMON/HIJCRDN/YP(3,300),YT(3,300)
0825 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0826 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
0827 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0828 EXTERNAL FNKICK,FNKICK2,FNSTRU,FNSTRUM,FNSTRUS
0829 SAVE
0830
0831 CALL TITLE
0832 IHNT2(1)=IAP
0833 IHNT2(2)=IZP
0834 IHNT2(3)=IAT
0835 IHNT2(4)=IZT
0836 IHNT2(5)=0
0837 IHNT2(6)=0
0838
0839 HINT1(8)=MAX(ULMASS(2112),ULMASS(2212))
0840 HINT1(9)=HINT1(8)
0841
0842 IF(PROJ.NE.'A') THEN
0843 IF(PROJ.EQ.'P') THEN
0844 IHNT2(5)=2212
0845 ELSE IF(PROJ.EQ.'PBAR') THEN
0846 IHNT2(5)=-2212
0847 ELSE IF(PROJ.EQ.'PI+') THEN
0848 IHNT2(5)=211
0849 ELSE IF(PROJ.EQ.'PI-') THEN
0850 IHNT2(5)=-211
0851 ELSE IF(PROJ.EQ.'K+') THEN
0852 IHNT2(5)=321
0853 ELSE IF(PROJ.EQ.'K-') THEN
0854 IHNT2(5)=-321
0855 ELSE IF(PROJ.EQ.'N') THEN
0856 IHNT2(5)=2112
0857 ELSE IF(PROJ.EQ.'NBAR') THEN
0858 IHNT2(5)=-2112
0859 ELSE
0860 WRITE(6,*) PROJ, 'wrong or unavailable proj name'
0861 STOP
0862 ENDIF
0863 HINT1(8)=ULMASS(IHNT2(5))
0864 ENDIF
0865 IF(TARG.NE.'A') THEN
0866 IF(TARG.EQ.'P') THEN
0867 IHNT2(6)=2212
0868 ELSE IF(TARG.EQ.'PBAR') THEN
0869 IHNT2(6)=-2212
0870 ELSE IF(TARG.EQ.'PI+') THEN
0871 IHNT2(6)=211
0872 ELSE IF(TARG.EQ.'PI-') THEN
0873 IHNT2(6)=-211
0874 ELSE IF(TARG.EQ.'K+') THEN
0875 IHNT2(6)=321
0876 ELSE IF(TARG.EQ.'K-') THEN
0877 IHNT2(6)=-321
0878 ELSE IF(TARG.EQ.'N') THEN
0879 IHNT2(6)=2112
0880 ELSE IF(TARG.EQ.'NBAR') THEN
0881 IHNT2(6)=-2112
0882 ELSE
0883 WRITE(6,*) TARG,'wrong or unavailable targ name'
0884 STOP
0885 ENDIF
0886 HINT1(9)=ULMASS(IHNT2(6))
0887 ENDIF
0888
0889
0890
0891 IF(IHPR2(12).GT.0) THEN
0892 CALL LUGIVE('MDCY(C111,1)=0')
0893 CALL LUGIVE('MDCY(C310,1)=0')
0894 CALL LUGIVE('MDCY(C411,1)=0;MDCY(C-411,1)=0')
0895 CALL LUGIVE('MDCY(C421,1)=0;MDCY(C-421,1)=0')
0896 CALL LUGIVE('MDCY(C431,1)=0;MDCY(C-431,1)=0')
0897 CALL LUGIVE('MDCY(C511,1)=0;MDCY(C-511,1)=0')
0898 CALL LUGIVE('MDCY(C521,1)=0;MDCY(C-521,1)=0')
0899 CALL LUGIVE('MDCY(C531,1)=0;MDCY(C-531,1)=0')
0900 CALL LUGIVE('MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
0901 CALL LUGIVE('MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
0902 CALL LUGIVE('MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
0903 CALL LUGIVE('MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
0904 CALL LUGIVE('MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
0905 CALL LUGIVE('MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
0906 CALL LUGIVE('MDCY(C3334,1)=0;MDCY(C-3334,1)=0')
0907 ENDIF
0908 MSTU(12)=0
0909 MSTU(21)=1
0910 IF(IHPR2(10).EQ.0) THEN
0911 MSTU(22)=0
0912 MSTU(25)=0
0913 MSTU(26)=0
0914 ENDIF
0915 MSTJ(12)=IHPR2(11)
0916 PARJ(21)=HIPR1(2)
0917 PARJ(41)=HIPR1(3)
0918 PARJ(42)=HIPR1(4)
0919
0920 IF(FRAME.EQ.'LAB') THEN
0921 DD1=EFRM
0922 DD2=HINT1(8)
0923 DD3=HINT1(9)
0924 HINT1(1)=SQRT(HINT1(8)**2+2.0*HINT1(9)*EFRM+HINT1(9)**2)
0925 DD4=DSQRT(DD1**2-DD2**2)/(DD1+DD3)
0926 HINT1(2)=DD4
0927 HINT1(3)=0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
0928 DD4=DSQRT(DD1**2-DD2**2)/DD1
0929 HINT1(4)=0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
0930 HINT1(5)=0.0
0931 HINT1(6)=EFRM
0932 HINT1(7)=HINT1(9)
0933 ELSE IF(FRAME.EQ.'CMS') THEN
0934 HINT1(1)=EFRM
0935 HINT1(2)=0.0
0936 HINT1(3)=0.0
0937 DD1=HINT1(1)
0938 DD2=HINT1(8)
0939 DD3=HINT1(9)
0940 DD4=DSQRT(1.D0-4.D0*DD2**2/DD1**2)
0941 HINT1(4)=0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
0942 DD4=DSQRT(1.D0-4.D0*DD3**2/DD1**2)
0943 HINT1(5)=-0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
0944 HINT1(6)=HINT1(1)/2.0
0945 HINT1(7)=HINT1(1)/2.0
0946 ENDIF
0947
0948
0949
0950
0951 IF(IHNT2(1).GT.1) THEN
0952 CALL HIJWDS(IHNT2(1),1,RMAX)
0953 HIPR1(34)=RMAX
0954
0955 ENDIF
0956 IF(IHNT2(3).GT.1) THEN
0957 CALL HIJWDS(IHNT2(3),2,RMAX)
0958 HIPR1(35)=RMAX
0959
0960 ENDIF
0961
0962
0963 I=0
0964 20 I=I+1
0965 IF(I.EQ.10) GO TO 30
0966 IF(HIDAT0(10,I).LE.HINT1(1)) GO TO 20
0967 30 IF(I.EQ.1) I=2
0968 DO 40 J=1,9
0969 HIDAT(J)=HIDAT0(J,I-1)+(HIDAT0(J,I)-HIDAT0(J,I-1))
0970 & *(HINT1(1)-HIDAT0(10,I-1))/(HIDAT0(10,I)-HIDAT0(10,I-1))
0971 40 CONTINUE
0972 HIPR1(31)=HIDAT(5)
0973 HIPR1(30)=2.0*HIDAT(5)
0974
0975
0976 CALL HIJCRS
0977
0978 IF(IHPR2(5).NE.0) THEN
0979 CALL HIFUN(3,0.0,36.0,FNKICK)
0980
0981 ENDIF
0982 CALL HIFUN(7,0.0,6.0,FNKICK2)
0983 CALL HIFUN(4,0.0,1.0,FNSTRU)
0984 CALL HIFUN(5,0.0,1.0,FNSTRUM)
0985 CALL HIFUN(6,0.0,1.0,FNSTRUS)
0986
0987 EFRAME='Ecm'
0988 IF(FRAME.EQ.'LAB') EFRAME='Elab'
0989 WRITE(6,100) EFRAME,EFRM,PROJ,IHNT2(1),IHNT2(2),
0990 & TARG,IHNT2(3),IHNT2(4)
0991 100 FORMAT(//10X,'****************************************
0992 & **********'/
0993 & 10X,'*',48X,'*'/
0994 & 10X,'* HIJING has been initialized at *'/
0995 & 10X,'*',13X,A4,'= ',F10.2,' GeV/n',13X,'*'/
0996 & 10X,'*',48X,'*'/
0997 & 10X,'*',8X,'for ',
0998 & A4,'(',I3,',',I3,')',' + ',A4,'(',I3,',',I3,')',7X,'*'/
0999 & 10X,'**************************************************')
1000 RETURN
1001 END
1002
1003
1004
1005 FUNCTION FNKICK(X)
1006 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1007 SAVE
1008 FNKICK=1.0/(X+HIPR1(19)**2)/(X+HIPR1(20)**2)
1009 & /(1+EXP((SQRT(X)-HIPR1(20))/0.4))
1010 RETURN
1011 END
1012
1013
1014 FUNCTION FNKICK2(X)
1015 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1016 SAVE
1017 FNKICK2=X*EXP(-2.0*X/HIPR1(42))
1018 RETURN
1019 END
1020
1021
1022
1023 FUNCTION FNSTRU(X)
1024 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1025 SAVE
1026 FNSTRU=(1.0-X)**HIPR1(44)/
1027 & (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
1028 RETURN
1029 END
1030
1031
1032
1033 FUNCTION FNSTRUM(X)
1034 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1035 SAVE
1036 FNSTRUM=1.0/((1.0-X)**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
1037 & /(X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
1038 RETURN
1039 END
1040
1041
1042 FUNCTION FNSTRUS(X)
1043 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1044 SAVE
1045 FNSTRUS=(1.0-X)**HIPR1(47)/
1046 & (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(48)
1047 RETURN
1048 END
1049
1050
1051
1052
1053 SUBROUTINE HIBOOST
1054 IMPLICIT DOUBLE PRECISION(D)
1055 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
1056 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1057 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1058 SAVE
1059 DO 100 I=1,N
1060 DBETA=P(I,3)/P(I,4)
1061 IF(ABS(DBETA).GE.1.D0) THEN
1062 DB=HINT1(2)
1063 IF(DB.GT.0.99999999D0) THEN
1064
1065 WRITE(6,*) '(HIBOOT:) boost vector too large'
1066 DB=0.99999999D0
1067 ENDIF
1068 DGA=1D0/SQRT(1D0-DB**2)
1069 DP3=P(I,3)
1070 DP4=P(I,4)
1071 P(I,3)=(DP3+DB*DP4)*DGA
1072 P(I,4)=(DP4+DB*DP3)*DGA
1073 GO TO 100
1074 ENDIF
1075 Y=0.5*DLOG((1.D0+DBETA)/(1.D0-DBETA))
1076 AMT=SQRT(P(I,1)**2+P(I,2)**2+P(I,5)**2)
1077 P(I,3)=AMT*SINH(Y+HINT1(3))
1078 P(I,4)=AMT*COSH(Y+HINT1(3))
1079 100 CONTINUE
1080 RETURN
1081 END
1082
1083
1084
1085
1086 SUBROUTINE QUENCH(JPJT,NTP)
1087 DIMENSION RDP(300),LQP(300),RDT(300),LQT(300)
1088 COMMON/HIJCRDN/YP(3,300),YT(3,300)
1089 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1090
1091 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
1092 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
1093 & PJPM(300,500),NTJ(300),KFTJ(300,500),
1094 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
1095 & PJTE(300,500),PJTM(300,500)
1096 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
1097 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
1098 & PZSG(900,100),PESG(900,100),PMSG(900,100)
1099 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
1100 COMMON/RANSEED/NSEED
1101 SAVE
1102
1103 BB=HINT1(19) ! Uzhi
1104 PHI=HINT1(20) ! Uzhi
1105 BBX=BB*COS(PHI) ! Uzhi
1106 BBY=BB*SIN(PHI) ! Uzhi
1107
1108 IF(NTP.EQ.2) GO TO 400
1109 IF(NTP.EQ.3) GO TO 2000
1110
1111
1112
1113
1114 IF(NFP(JPJT,7).NE.1) RETURN
1115
1116 JP=JPJT
1117 DO 290 I=1,NPJ(JP)
1118 PTJET0=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2)
1119 IF(PTJET0.LE.HIPR1(11)) GO TO 290
1120 PTOT=SQRT(PTJET0*PTJET0+PJPZ(JP,I)**2)
1121 IF(PTOT.LT.HIPR1(8)) GO TO 290
1122 PHIP=ULANGL(PJPX(JP,I),PJPY(JP,I))
1123
1124 KP=0
1125 DO 100 I2=1,IHNT2(1)
1126 IF(NFP(I2,5).NE.3 .OR. I2.EQ.JP) GO TO 100
1127 DX=YP(1,I2)-YP(1,JP)
1128 DY=YP(2,I2)-YP(2,JP)
1129 PHI=ULANGL(DX,DY)
1130 DPHI=ABS(PHI-PHIP)
1131 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
1132 IF(DPHI.GE.HIPR1(40)/2.0) GO TO 100
1133 RD0=SQRT(DX*DX+DY*DY)
1134 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 100
1135 KP=KP+1
1136 LQP(KP)=I2
1137 RDP(KP)=COS(DPHI)*RD0
1138 100 CONTINUE
1139
1140 DO 110 I2=1,KP-1
1141 DO 110 J2=I2+1,KP
1142 IF(RDP(I2).LT.RDP(J2)) GO TO 110
1143 RD=RDP(I2)
1144 LQ=LQP(I2)
1145 RDP(I2)=RDP(J2)
1146 LQP(I2)=LQP(J2)
1147 RDP(J2)=RD
1148 LQP(J2)=LQ
1149 110 CONTINUE
1150
1151 KT=0
1152 DO 120 I2=1,IHNT2(3)
1153 IF(NFT(I2,5).NE.3) GO TO 120
1154 DX=YT(1,I2)-YP(1,JP)-BBX
1155 DY=YT(2,I2)-YP(2,JP)-BBY
1156 PHI=ULANGL(DX,DY)
1157 DPHI=ABS(PHI-PHIP)
1158 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
1159 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 120
1160 RD0=SQRT(DX*DX+DY*DY)
1161 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 120
1162 KT=KT+1
1163 LQT(KT)=I2
1164 RDT(KT)=COS(DPHI)*RD0
1165 120 CONTINUE
1166
1167 DO 130 I2=1,KT-1
1168 DO 130 J2=I2+1,KT
1169 IF(RDT(I2).LT.RDT(J2)) GO TO 130
1170 RD=RDT(I2)
1171 LQ=LQT(I2)
1172 RDT(I2)=RDT(J2)
1173 LQT(I2)=LQT(J2)
1174 RDT(J2)=RD
1175 LQT(J2)=LQ
1176 130 CONTINUE
1177
1178 MP=0
1179 MT=0
1180 R0=0.0
1181 NQ=0
1182 DP=0.0
1183 PTOT=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2+PJPZ(JP,I)**2)
1184 V1=PJPX(JP,I)/PTOT
1185 V2=PJPY(JP,I)/PTOT
1186 V3=PJPZ(JP,I)/PTOT
1187
1188 200 RN=RLU(NSEED)
1189 210 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 290
1190 IF(MT.GE.KT) GO TO 220
1191 IF(MP.GE.KP) GO TO 240
1192 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 240
1193 220 MP=MP+1
1194 DRR=RDP(MP)-R0
1195 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
1196 DP=DRR*HIPR1(14)
1197 IF(KFPJ(JP,I).NE.21) DP=0.5*DP
1198
1199 IF(DP.LE.0.2) GO TO 210
1200 IF(PTOT.LE.0.4) GO TO 290
1201 IF(PTOT.LE.DP) DP=PTOT-0.2
1202 DE=DP
1203
1204 IF(KFPJ(JP,I).NE.21) THEN
1205 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
1206 & +PP(LQP(MP),3)**2
1207 DE=SQRT(PJPM(JP,I)**2+PTOT**2)
1208 & -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
1209 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
1210 AMSHU=ERSHU-PRSHU
1211 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
1212 PP(LQP(MP),4)=SQRT(ERSHU)
1213 PP(LQP(MP),5)=SQRT(AMSHU)
1214 ENDIF
1215
1216 R0=RDP(MP)
1217 DP1=DP*V1
1218 DP2=DP*V2
1219 DP3=DP*V3
1220
1221
1222 NPJ(LQP(MP))=NPJ(LQP(MP))+1
1223 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
1224 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
1225 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
1226 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
1227 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
1228 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
1229 GO TO 260
1230
1231 240 MT=MT+1
1232 DRR=RDT(MT)-R0
1233 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
1234 DP=DRR*HIPR1(14)
1235 IF(DP.LE.0.2) GO TO 210
1236 IF(PTOT.LE.0.4) GO TO 290
1237 IF(PTOT.LE.DP) DP=PTOT-0.2
1238 DE=DP
1239
1240 IF(KFPJ(JP,I).NE.21) THEN
1241 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
1242 & +PT(LQT(MT),3)**2
1243 DE=SQRT(PJPM(JP,I)**2+PTOT**2)
1244 & -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
1245 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
1246 AMSHU=ERSHU-PRSHU
1247 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
1248 PT(LQT(MT),4)=SQRT(ERSHU)
1249 PT(LQT(MT),5)=SQRT(AMSHU)
1250 ENDIF
1251
1252
1253 R0=RDT(MT)
1254 DP1=DP*V1
1255 DP2=DP*V2
1256 DP3=DP*V3
1257
1258 NTJ(LQT(MT))=NTJ(LQT(MT))+1
1259 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
1260 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
1261 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
1262 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
1263 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
1264 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
1265
1266 260 PJPX(JP,I)=(PTOT-DP)*V1
1267 PJPY(JP,I)=(PTOT-DP)*V2
1268 PJPZ(JP,I)=(PTOT-DP)*V3
1269 PJPE(JP,I)=PJPE(JP,I)-DE
1270
1271 PTOT=PTOT-DP
1272 NQ=NQ+1
1273 GO TO 200
1274 290 CONTINUE
1275
1276 RETURN
1277
1278
1279
1280
1281
1282
1283
1284 400 IF(NFT(JPJT,7).NE.1) RETURN
1285 JT=JPJT
1286 DO 690 I=1,NTJ(JT)
1287 PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
1288 IF(PTJET0.LE.HIPR1(11)) GO TO 690
1289 PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
1290 IF(PTOT.LT.HIPR1(8)) GO TO 690
1291 PHIT=ULANGL(PJTX(JT,I),PJTY(JT,I))
1292 KP=0
1293 DO 500 I2=1,IHNT2(1)
1294 IF(NFP(I2,5).NE.3) GO TO 500
1295 DX=YP(1,I2)+BBX-YT(1,JT)
1296 DY=YP(2,I2)+BBY-YT(2,JT)
1297 PHI=ULANGL(DX,DY)
1298 DPHI=ABS(PHI-PHIT)
1299 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
1300 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 500
1301 RD0=SQRT(DX*DX+DY*DY)
1302 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 500
1303 KP=KP+1
1304 LQP(KP)=I2
1305 RDP(KP)=COS(DPHI)*RD0
1306 500 CONTINUE
1307
1308 DO 510 I2=1,KP-1
1309 DO 510 J2=I2+1,KP
1310 IF(RDP(I2).LT.RDP(J2)) GO TO 510
1311 RD=RDP(I2)
1312 LQ=LQP(I2)
1313 RDP(I2)=RDP(J2)
1314 LQP(I2)=LQP(J2)
1315 RDP(J2)=RD
1316 LQP(J2)=LQ
1317 510 CONTINUE
1318
1319 KT=0
1320 DO 520 I2=1,IHNT2(3)
1321 IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 520
1322 DX=YT(1,I2)-YT(1,JT)
1323 DY=YT(2,I2)-YT(2,JT)
1324 PHI=ULANGL(DX,DY)
1325 DPHI=ABS(PHI-PHIT)
1326 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
1327 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 520
1328 RD0=SQRT(DX*DX+DY*DY)
1329 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 520
1330 KT=KT+1
1331 LQT(KT)=I2
1332 RDT(KT)=COS(DPHI)*RD0
1333 520 CONTINUE
1334
1335 DO 530 I2=1,KT-1
1336 DO 530 J2=I2+1,KT
1337 IF(RDT(I2).LT.RDT(J2)) GO TO 530
1338 RD=RDT(I2)
1339 LQ=LQT(I2)
1340 RDT(I2)=RDT(J2)
1341 LQT(I2)=LQT(J2)
1342 RDT(J2)=RD
1343 LQT(J2)=LQ
1344 530 CONTINUE
1345
1346 MP=0
1347 MT=0
1348 NQ=0
1349 DP=0.0
1350 R0=0.0
1351 PTOT=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2+PJTZ(JT,I)**2)
1352 V1=PJTX(JT,I)/PTOT
1353 V2=PJTY(JT,I)/PTOT
1354 V3=PJTZ(JT,I)/PTOT
1355
1356 600 RN=RLU(NSEED)
1357 610 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 690
1358 IF(MT.GE.KT) GO TO 620
1359 IF(MP.GE.KP) GO TO 640
1360 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 640
1361 620 MP=MP+1
1362 DRR=RDP(MP)-R0
1363 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
1364 DP=DRR*HIPR1(14)
1365 IF(KFTJ(JT,I).NE.21) DP=0.5*DP
1366
1367 IF(DP.LE.0.2) GO TO 610
1368 IF(PTOT.LE.0.4) GO TO 690
1369 IF(PTOT.LE.DP) DP=PTOT-0.2
1370 DE=DP
1371
1372 IF(KFTJ(JT,I).NE.21) THEN
1373 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
1374 & +PP(LQP(MP),3)**2
1375 DE=SQRT(PJTM(JT,I)**2+PTOT**2)
1376 & -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
1377 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
1378 AMSHU=ERSHU-PRSHU
1379 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
1380 PP(LQP(MP),4)=SQRT(ERSHU)
1381 PP(LQP(MP),5)=SQRT(AMSHU)
1382 ENDIF
1383
1384
1385 R0=RDP(MP)
1386 DP1=DP*V1
1387 DP2=DP*V2
1388 DP3=DP*V3
1389
1390 NPJ(LQP(MP))=NPJ(LQP(MP))+1
1391 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
1392 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
1393 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
1394 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
1395 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
1396 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
1397
1398 GO TO 660
1399
1400 640 MT=MT+1
1401 DRR=RDT(MT)-R0
1402 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
1403 DP=DRR*HIPR1(14)
1404 IF(DP.LE.0.2) GO TO 610
1405 IF(PTOT.LE.0.4) GO TO 690
1406 IF(PTOT.LE.DP) DP=PTOT-0.2
1407 DE=DP
1408
1409 IF(KFTJ(JT,I).NE.21) THEN
1410 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
1411 & +PT(LQT(MT),3)**2
1412 DE=SQRT(PJTM(JT,I)**2+PTOT**2)
1413 & -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
1414 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
1415 AMSHU=ERSHU-PRSHU
1416 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
1417 PT(LQT(MT),4)=SQRT(ERSHU)
1418 PT(LQT(MT),5)=SQRT(AMSHU)
1419 ENDIF
1420
1421
1422 R0=RDT(MT)
1423 DP1=DP*V1
1424 DP2=DP*V2
1425 DP3=DP*V3
1426
1427 NTJ(LQT(MT))=NTJ(LQT(MT))+1
1428 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
1429 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
1430 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
1431 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
1432 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
1433 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
1434
1435 660 PJTX(JT,I)=(PTOT-DP)*V1
1436 PJTY(JT,I)=(PTOT-DP)*V2
1437 PJTZ(JT,I)=(PTOT-DP)*V3
1438 PJTE(JT,I)=PJTE(JT,I)-DE
1439
1440 PTOT=PTOT-DP
1441 NQ=NQ+1
1442 GO TO 600
1443 690 CONTINUE
1444 RETURN
1445
1446
1447
1448 2000 ISG=JPJT
1449 IF(IASG(ISG,3).NE.1) RETURN
1450
1451 JP=IASG(ISG,1)
1452 JT=IASG(ISG,2)
1453 XJ=(YP(1,JP)+BBX+YT(1,JT))/2.0
1454 YJ=(YP(2,JP)+BBY+YT(2,JT))/2.0
1455 DO 2690 I=1,NJSG(ISG)
1456 PTJET0=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2)
1457 IF(PTJET0.LE.HIPR1(11).OR.PESG(ISG,I).LT.HIPR1(1))
1458 & GO TO 2690
1459 PTOT=SQRT(PTJET0*PTJET0+PZSG(ISG,I)**2)
1460 IF(PTOT.LT.MAX(HIPR1(1),HIPR1(8))) GO TO 2690
1461 PHIQ=ULANGL(PXSG(ISG,I),PYSG(ISG,I))
1462 KP=0
1463 DO 2500 I2=1,IHNT2(1)
1464 IF(NFP(I2,5).NE.3.OR.I2.EQ.JP) GO TO 2500
1465 DX=YP(1,I2)+BBX-XJ
1466 DY=YP(2,I2)+BBY-YJ
1467 PHI=ULANGL(DX,DY)
1468 DPHI=ABS(PHI-PHIQ)
1469 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
1470 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2500
1471 RD0=SQRT(DX*DX+DY*DY)
1472 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2500
1473 KP=KP+1
1474 LQP(KP)=I2
1475 RDP(KP)=COS(DPHI)*RD0
1476 2500 CONTINUE
1477
1478 DO 2510 I2=1,KP-1
1479 DO 2510 J2=I2+1,KP
1480 IF(RDP(I2).LT.RDP(J2)) GO TO 2510
1481 RD=RDP(I2)
1482 LQ=LQP(I2)
1483 RDP(I2)=RDP(J2)
1484 LQP(I2)=LQP(J2)
1485 RDP(J2)=RD
1486 LQP(J2)=LQ
1487 2510 CONTINUE
1488
1489 KT=0
1490 DO 2520 I2=1,IHNT2(3)
1491 IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 2520
1492 DX=YT(1,I2)-XJ
1493 DY=YT(2,I2)-YJ
1494 PHI=ULANGL(DX,DY)
1495 DPHI=ABS(PHI-PHIQ)
1496 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
1497 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2520
1498 RD0=SQRT(DX*DX+DY*DY)
1499 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2520
1500 KT=KT+1
1501 LQT(KT)=I2
1502 RDT(KT)=COS(DPHI)*RD0
1503 2520 CONTINUE
1504
1505 DO 2530 I2=1,KT-1
1506 DO 2530 J2=I2+1,KT
1507 IF(RDT(I2).LT.RDT(J2)) GO TO 2530
1508 RD=RDT(I2)
1509 LQ=LQT(I2)
1510 RDT(I2)=RDT(J2)
1511 LQT(I2)=LQT(J2)
1512 RDT(J2)=RD
1513 LQT(J2)=LQ
1514 2530 CONTINUE
1515
1516 MP=0
1517 MT=0
1518 NQ=0
1519 DP=0.0
1520 R0=0.0
1521 PTOT=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2
1522 & +PZSG(ISG,I)**2)
1523 V1=PXSG(ISG,I)/PTOT
1524 V2=PYSG(ISG,I)/PTOT
1525 V3=PZSG(ISG,I)/PTOT
1526
1527 2600 RN=RLU(NSEED)
1528 2610 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 2690
1529 IF(MT.GE.KT) GO TO 2620
1530 IF(MP.GE.KP) GO TO 2640
1531 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 2640
1532 2620 MP=MP+1
1533 DRR=RDP(MP)-R0
1534 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
1535 DP=DRR*HIPR1(14)/2.0
1536 IF(DP.LE.0.2) GO TO 2610
1537 IF(PTOT.LE.0.4) GO TO 2690
1538 IF(PTOT.LE.DP) DP=PTOT-0.2
1539 DE=DP
1540
1541 IF(K2SG(ISG,I).NE.21) THEN
1542 IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
1543 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
1544 & +PP(LQP(MP),3)**2
1545 DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
1546 & -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
1547 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
1548 AMSHU=ERSHU-PRSHU
1549 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
1550 PP(LQP(MP),4)=SQRT(ERSHU)
1551 PP(LQP(MP),5)=SQRT(AMSHU)
1552 ENDIF
1553
1554
1555 R0=RDP(MP)
1556 DP1=DP*V1
1557 DP2=DP*V2
1558 DP3=DP*V3
1559
1560 NPJ(LQP(MP))=NPJ(LQP(MP))+1
1561 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
1562 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
1563 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
1564 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
1565 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
1566 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
1567
1568 GO TO 2660
1569
1570 2640 MT=MT+1
1571 DRR=RDT(MT)-R0
1572 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
1573 DP=DRR*HIPR1(14)
1574 IF(DP.LE.0.2) GO TO 2610
1575 IF(PTOT.LE.0.4) GO TO 2690
1576 IF(PTOT.LE.DP) DP=PTOT-0.2
1577 DE=DP
1578
1579 IF(K2SG(ISG,I).NE.21) THEN
1580 IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
1581 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
1582 & +PT(LQT(MT),3)**2
1583 DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
1584 & -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
1585 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
1586 AMSHU=ERSHU-PRSHU
1587 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
1588 PT(LQT(MT),4)=SQRT(ERSHU)
1589 PT(LQT(MT),5)=SQRT(AMSHU)
1590 ENDIF
1591
1592
1593 R0=RDT(MT)
1594 DP1=DP*V1
1595 DP2=DP*V2
1596 DP3=DP*V3
1597
1598 NTJ(LQT(MT))=NTJ(LQT(MT))+1
1599 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
1600 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
1601 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
1602 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
1603 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
1604 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
1605
1606 2660 PXSG(ISG,I)=(PTOT-DP)*V1
1607 PYSG(ISG,I)=(PTOT-DP)*V2
1608 PZSG(ISG,I)=(PTOT-DP)*V3
1609 PESG(ISG,I)=PESG(ISG,I)-DE
1610
1611 PTOT=PTOT-DP
1612 NQ=NQ+1
1613 GO TO 2600
1614 2690 CONTINUE
1615 RETURN
1616 END
1617
1618
1619
1620
1621
1622 SUBROUTINE HIJFRG(JTP,NTP,IERROR)
1623
1624
1625
1626
1627
1628
1629 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1630 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
1631 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
1632 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
1633 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
1634 & PJPM(300,500),NTJ(300),KFTJ(300,500),
1635 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
1636 & PJTE(300,500),PJTM(300,500)
1637 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
1638 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
1639 & PZSG(900,100),PESG(900,100),PMSG(900,100)
1640
1641 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
1642 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1643 COMMON/RANSEED/NSEED
1644
1645 IERROR=0
1646 CALL LUEDIT(0)
1647 N=0
1648
1649 IF(NTP.EQ.3) THEN
1650 ISG=JTP
1651 N=NJSG(ISG)
1652 DO 100 I=1,NJSG(ISG)
1653 K(I,1)=K1SG(ISG,I)
1654 K(I,2)=K2SG(ISG,I)
1655 P(I,1)=PXSG(ISG,I)
1656 P(I,2)=PYSG(ISG,I)
1657 P(I,3)=PZSG(ISG,I)
1658 P(I,4)=PESG(ISG,I)
1659 P(I,5)=PMSG(ISG,I)
1660
1661 V(I,1)=0
1662 V(I,2)=0
1663 V(I,3)=0
1664 V(I,4)=0
1665 V(I,5)=0
1666 100 CONTINUE
1667
1668
1669
1670 CALL LUEXEC
1671 RETURN
1672 ENDIF
1673
1674 IF(NTP.EQ.2) GO TO 200
1675 IF(JTP.GT.IHNT2(1)) RETURN
1676 IF(NFP(JTP,5).NE.3.AND.NFP(JTP,3).NE.0
1677 & .AND.NPJ(JTP).EQ.0.AND.NFP(JTP,10).EQ.0) GO TO 1000
1678 IF(NFP(JTP,15).EQ.-1) THEN
1679 KF1=NFP(JTP,2)
1680 KF2=NFP(JTP,1)
1681 PQ21=PP(JTP,6)
1682 PQ22=PP(JTP,7)
1683 PQ11=PP(JTP,8)
1684 PQ12=PP(JTP,9)
1685 AM1=PP(JTP,15)
1686 AM2=PP(JTP,14)
1687 ELSE
1688 KF1=NFP(JTP,1)
1689 KF2=NFP(JTP,2)
1690 PQ21=PP(JTP,8)
1691 PQ22=PP(JTP,9)
1692 PQ11=PP(JTP,6)
1693 PQ12=PP(JTP,7)
1694 AM1=PP(JTP,14)
1695 AM2=PP(JTP,15)
1696 ENDIF
1697
1698 PB1=PQ11+PQ21
1699 PB2=PQ12+PQ22
1700 PB3=PP(JTP,3)
1701 PECM=PP(JTP,5)
1702 BTZ=PB3/PP(JTP,4)
1703 IF((ABS(PB1-PP(JTP,1)).GT.0.01.OR.
1704 & ABS(PB2-PP(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
1705 & WRITE(6,*) ' Pt of Q and QQ do not sum to the total'
1706
1707 GO TO 300
1708
1709 200 IF(JTP.GT.IHNT2(3)) RETURN
1710 IF(NFT(JTP,5).NE.3.AND.NFT(JTP,3).NE.0
1711 & .AND.NTJ(JTP).EQ.0.AND.NFT(JTP,10).EQ.0) GO TO 1200
1712 IF(NFT(JTP,15).EQ.1) THEN
1713 KF1=NFT(JTP,1)
1714 KF2=NFT(JTP,2)
1715 PQ11=PT(JTP,6)
1716 PQ12=PT(JTP,7)
1717 PQ21=PT(JTP,8)
1718 PQ22=PT(JTP,9)
1719 AM1=PT(JTP,14)
1720 AM2=PT(JTP,15)
1721 ELSE
1722 KF1=NFT(JTP,2)
1723 KF2=NFT(JTP,1)
1724 PQ11=PT(JTP,8)
1725 PQ12=PT(JTP,9)
1726 PQ21=PT(JTP,6)
1727 PQ22=PT(JTP,7)
1728 AM1=PT(JTP,15)
1729 AM2=PT(JTP,14)
1730 ENDIF
1731
1732 PB1=PQ11+PQ21
1733 PB2=PQ12+PQ22
1734 PB3=PT(JTP,3)
1735 PECM=PT(JTP,5)
1736 BTZ=PB3/PT(JTP,4)
1737
1738 IF((ABS(PB1-PT(JTP,1)).GT.0.01.OR.
1739 & ABS(PB2-PT(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
1740 & WRITE(6,*) ' Pt of Q and QQ do not sum to the total'
1741
1742 300 IF(PECM.LT.HIPR1(1)) THEN
1743 IERROR=1
1744 IF(IHPR2(10).EQ.0) RETURN
1745 WRITE(6,*) ' ECM=',PECM,' energy of the string is too small'
1746 RETURN
1747 ENDIF
1748 AMT=PECM**2+PB1**2+PB2**2
1749 AMT1=AM1**2+PQ11**2+PQ12**2
1750 AMT2=AM2**2+PQ21**2+PQ22**2
1751 PZCM=SQRT(ABS(AMT**2+AMT1**2+AMT2**2-2.0*AMT*AMT1
1752 & -2.0*AMT*AMT2-2.0*AMT1*AMT2))/2.0/SQRT(AMT)
1753
1754 K(1,1)=2
1755 K(1,2)=KF1
1756 P(1,1)=PQ11
1757 P(1,2)=PQ12
1758 P(1,3)=PZCM
1759 P(1,4)=SQRT(AMT1+PZCM**2)
1760 P(1,5)=AM1
1761 K(2,1)=1
1762 K(2,2)=KF2
1763 P(2,1)=PQ21
1764 P(2,2)=PQ22
1765 P(2,3)=-PZCM
1766 P(2,4)=SQRT(AMT2+PZCM**2)
1767 P(2,5)=AM2
1768
1769 V(1,1)=0
1770 V(1,2)=0
1771 V(1,3)=0
1772 V(1,4)=0
1773 V(1,5)=0
1774 V(2,1)=0
1775 V(2,2)=0
1776 V(2,3)=0
1777 V(2,4)=0
1778 V(2,5)=0
1779 N=2
1780
1781 CALL HIROBO(0.0,0.0,0.0,0.0,BTZ)
1782 JETOT=0
1783 IF((PQ21**2+PQ22**2).GT.(PQ11**2+PQ12**2)) THEN
1784 PMAX1=P(2,1)
1785 PMAX2=P(2,2)
1786 PMAX3=P(2,3)
1787 ELSE
1788 PMAX1=P(1,1)
1789 PMAX2=P(1,2)
1790 PMAX3=P(1,3)
1791 ENDIF
1792 IF(NTP.EQ.1) THEN
1793 PP(JTP,10)=PMAX1
1794 PP(JTP,11)=PMAX2
1795 PP(JTP,12)=PMAX3
1796 ELSE IF(NTP.EQ.2) THEN
1797 PT(JTP,10)=PMAX1
1798 PT(JTP,11)=PMAX2
1799 PT(JTP,12)=PMAX3
1800 ENDIF
1801
1802 IF(NTP.EQ.1.AND.NPJ(JTP).NE.0) THEN
1803 JETOT=NPJ(JTP)
1804
1805
1806 IEX=0
1807 IF((ABS(KF1).GT.1000.AND.KF1.LT.0)
1808 & .OR.(ABS(KF1).LT.1000.AND.KF1.GT.0)) IEX=1
1809 DO 520 I=N,2,-1
1810 DO 520 J=1,5
1811 II=NPJ(JTP)+I
1812 K(II,J)=K(I,J)
1813 P(II,J)=P(I,J)
1814 V(II,J)=V(I,J)
1815 520 CONTINUE
1816 DO 540 I=1,NPJ(JTP)
1817 DO 542 J=1,5
1818 K(I+1,J)=0
1819 V(I+1,J)=0
1820 542 CONTINUE
1821 I0=I
1822 IF(IEX.EQ.1) I0=NPJ(JTP)-I+1
1823
1824 KK1=KFPJ(JTP,I0)
1825 K(I+1,1)=2
1826 K(I+1,2)=KK1
1827 IF(KK1.NE.21 .AND. KK1.NE.0) K(I+1,1)=
1828 & 1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
1829 P(I+1,1)=PJPX(JTP,I0)
1830 P(I+1,2)=PJPY(JTP,I0)
1831 P(I+1,3)=PJPZ(JTP,I0)
1832 P(I+1,4)=PJPE(JTP,I0)
1833 P(I+1,5)=PJPM(JTP,I0)
1834 540 CONTINUE
1835 N=N+NPJ(JTP)
1836 ELSE IF(NTP.EQ.2.AND.NTJ(JTP).NE.0) THEN
1837 JETOT=NTJ(JTP)
1838
1839
1840 IEX=1
1841 IF((ABS(KF2).GT.1000.AND.KF2.LT.0)
1842 & .OR.(ABS(KF2).LT.1000.AND.KF2.GT.0)) IEX=0
1843 DO 560 I=N,2,-1
1844 DO 560 J=1,5
1845 II=NTJ(JTP)+I
1846 K(II,J)=K(I,J)
1847 P(II,J)=P(I,J)
1848 V(II,J)=V(I,J)
1849 560 CONTINUE
1850 DO 580 I=1,NTJ(JTP)
1851 DO 582 J=1,5
1852 K(I+1,J)=0
1853 V(I+1,J)=0
1854 582 CONTINUE
1855 I0=I
1856 IF(IEX.EQ.1) I0=NTJ(JTP)-I+1
1857
1858 KK1=KFTJ(JTP,I0)
1859 K(I+1,1)=2
1860 K(I+1,2)=KK1
1861 IF(KK1.NE.21 .AND. KK1.NE.0) K(I+1,1)=
1862 & 1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
1863 P(I+1,1)=PJTX(JTP,I0)
1864 P(I+1,2)=PJTY(JTP,I0)
1865 P(I+1,3)=PJTZ(JTP,I0)
1866 P(I+1,4)=PJTE(JTP,I0)
1867 P(I+1,5)=PJTM(JTP,I0)
1868 580 CONTINUE
1869 N=N+NTJ(JTP)
1870 ENDIF
1871 IF(IHPR2(1).GT.0.AND.RLU(NSEED).LE.HIDAT(3)) THEN
1872 HIDAT20=HIDAT(2)
1873 HIPR150=HIPR1(5)
1874 IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
1875 & HIDAT(2)=2.0
1876 IF(HINT1(1).GE.1000.0.AND.JETOT.EQ.0)THEN
1877 HIDAT(2)=3.0
1878 HIPR1(5)=5.0
1879 ENDIF
1880 CALL ATTRAD(IERROR)
1881 HIDAT(2)=HIDAT20
1882 HIPR1(5)=HIPR150
1883 ELSE IF(JETOT.EQ.0.AND.IHPR2(1).GT.0.AND.
1884 & HINT1(1).GE.1000.0.AND.
1885 & RLU(NSEED).LE.0.8) THEN
1886 HIDAT20=HIDAT(2)
1887 HIPR150=HIPR1(5)
1888 HIDAT(2)=3.0
1889 HIPR1(5)=5.0
1890 IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
1891 & HIDAT(2)=2.0
1892 CALL ATTRAD(IERROR)
1893 HIDAT(2)=HIDAT20
1894 HIPR1(5)=HIPR150
1895 ENDIF
1896 IF(IERROR.NE.0) RETURN
1897
1898
1899
1900
1901
1902 CALL LUEXEC
1903 RETURN
1904
1905 1000 N=1
1906 K(1,1)=1
1907 K(1,2)=NFP(JTP,3)
1908 DO 1100 JJ=1,5
1909 P(1,JJ)=PP(JTP,JJ)
1910
1911 V(1,JJ)=0
1912 1100 CONTINUE
1913
1914 CALL LUEXEC
1915
1916 RETURN
1917
1918 1200 N=1
1919 K(1,1)=1
1920 K(1,2)=NFT(JTP,3)
1921 DO 1300 JJ=1,5
1922 P(1,JJ)=PT(JTP,JJ)
1923
1924 V(1,JJ)=0
1925 1300 CONTINUE
1926
1927 CALL LUEXEC
1928
1929 RETURN
1930 END
1931
1932
1933
1934
1935
1936
1937
1938 SUBROUTINE HIJSRT(JPJT,NPT)
1939 DIMENSION KF(100),PX(100),PY(100),PZ(100),PE(100),PM(100)
1940 DIMENSION Y(100),IP(100,2)
1941 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
1942 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
1943 & PJPM(300,500),NTJ(300),KFTJ(300,500),
1944 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
1945 & PJTE(300,500),PJTM(300,500)
1946 SAVE
1947 IF(NPT.EQ.2) GO TO 500
1948 JP=JPJT
1949 IQ=0
1950 I=1
1951 100 KF(I)=KFPJ(JP,I)
1952 PX(I)=PJPX(JP,I)
1953 PY(I)=PJPY(JP,I)
1954 PZ(I)=PJPZ(JP,I)
1955 PE(I)=PJPE(JP,I)
1956 PM(I)=PJPM(JP,I)
1957 Y(I-IQ)=0.5*ALOG((ABS(PE(I)+PZ(I))+1.E-5)
1958 & /(ABS(PE(I)-PZ(I))+1.E-5))
1959 IP(I-IQ,1)=I
1960 IP(I-IQ,2)=0
1961 IF(KF(I).NE.21) THEN
1962 IP(I-IQ,2)=1
1963 IQ=IQ+1
1964 I=I+1
1965 KF(I)=KFPJ(JP,I)
1966 PX(I)=PJPX(JP,I)
1967 PY(I)=PJPY(JP,I)
1968 PZ(I)=PJPZ(JP,I)
1969 PE(I)=PJPE(JP,I)
1970 PM(I)=PJPM(JP,I)
1971 ENDIF
1972 I=I+1
1973 IF(I.LE.NPJ(JP)) GO TO 100
1974
1975 DO 200 I=1,NPJ(JP)-IQ
1976 DO 200 J=I+1,NPJ(JP)-IQ
1977 IF(Y(I).GT.Y(J)) GO TO 200
1978 IP1=IP(I,1)
1979 IP2=IP(I,2)
1980 IP(I,1)=IP(J,1)
1981 IP(I,2)=IP(J,2)
1982 IP(J,1)=IP1
1983 IP(J,2)=IP2
1984 200 CONTINUE
1985
1986 IQQ=0
1987 I=1
1988 300 KFPJ(JP,I)=KF(IP(I-IQQ,1))
1989 PJPX(JP,I)=PX(IP(I-IQQ,1))
1990 PJPY(JP,I)=PY(IP(I-IQQ,1))
1991 PJPZ(JP,I)=PZ(IP(I-IQQ,1))
1992 PJPE(JP,I)=PE(IP(I-IQQ,1))
1993 PJPM(JP,I)=PM(IP(I-IQQ,1))
1994 IF(IP(I-IQQ,2).EQ.1) THEN
1995 KFPJ(JP,I+1)=KF(IP(I-IQQ,1)+1)
1996 PJPX(JP,I+1)=PX(IP(I-IQQ,1)+1)
1997 PJPY(JP,I+1)=PY(IP(I-IQQ,1)+1)
1998 PJPZ(JP,I+1)=PZ(IP(I-IQQ,1)+1)
1999 PJPE(JP,I+1)=PE(IP(I-IQQ,1)+1)
2000 PJPM(JP,I+1)=PM(IP(I-IQQ,1)+1)
2001 I=I+1
2002 IQQ=IQQ+1
2003 ENDIF
2004 I=I+1
2005 IF(I.LE.NPJ(JP)) GO TO 300
2006
2007 RETURN
2008
2009 500 JT=JPJT
2010 IQ=0
2011 I=1
2012 600 KF(I)=KFTJ(JT,I)
2013 PX(I)=PJTX(JT,I)
2014 PY(I)=PJTY(JT,I)
2015 PZ(I)=PJTZ(JT,I)
2016 PE(I)=PJTE(JT,I)
2017 PM(I)=PJTM(JT,I)
2018 Y(I-IQ)=0.5*ALOG((ABS(PE(I)+PZ(I))+1.E-5)
2019 & /(ABS(PE(I)-PZ(I))+1.E-5))
2020 IP(I-IQ,1)=I
2021 IP(I-IQ,2)=0
2022 IF(KF(I).NE.21) THEN
2023 IP(I-IQ,2)=1
2024 IQ=IQ+1
2025 I=I+1
2026 KF(I)=KFTJ(JT,I)
2027 PX(I)=PJTX(JT,I)
2028 PY(I)=PJTY(JT,I)
2029 PZ(I)=PJTZ(JT,I)
2030 PE(I)=PJTE(JT,I)
2031 PM(I)=PJTM(JT,I)
2032 ENDIF
2033 I=I+1
2034 IF(I.LE.NTJ(JT)) GO TO 600
2035
2036 DO 700 I=1,NTJ(JT)-IQ
2037 DO 700 J=I+1,NTJ(JT)-IQ
2038 IF(Y(I).LT.Y(J)) GO TO 700
2039 IP1=IP(I,1)
2040 IP2=IP(I,2)
2041 IP(I,1)=IP(J,1)
2042 IP(I,2)=IP(J,2)
2043 IP(J,1)=IP1
2044 IP(J,2)=IP2
2045 700 CONTINUE
2046
2047 IQQ=0
2048 I=1
2049 800 KFTJ(JT,I)=KF(IP(I-IQQ,1))
2050 PJTX(JT,I)=PX(IP(I-IQQ,1))
2051 PJTY(JT,I)=PY(IP(I-IQQ,1))
2052 PJTZ(JT,I)=PZ(IP(I-IQQ,1))
2053 PJTE(JT,I)=PE(IP(I-IQQ,1))
2054 PJTM(JT,I)=PM(IP(I-IQQ,1))
2055 IF(IP(I-IQQ,2).EQ.1) THEN
2056 KFTJ(JT,I+1)=KF(IP(I-IQQ,1)+1)
2057 PJTX(JT,I+1)=PX(IP(I-IQQ,1)+1)
2058 PJTY(JT,I+1)=PY(IP(I-IQQ,1)+1)
2059 PJTZ(JT,I+1)=PZ(IP(I-IQQ,1)+1)
2060 PJTE(JT,I+1)=PE(IP(I-IQQ,1)+1)
2061 PJTM(JT,I+1)=PM(IP(I-IQQ,1)+1)
2062 I=I+1
2063 IQQ=IQQ+1
2064 ENDIF
2065 I=I+1
2066 IF(I.LE.NTJ(JT)) GO TO 800
2067 RETURN
2068 END
2069
2070
2071
2072
2073
2074
2075
2076 SUBROUTINE ATTRAD(IERROR)
2077
2078 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2079 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
2080 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
2081 COMMON/RANSEED/NSEED
2082 SAVE
2083 IERROR=0
2084
2085
2086
2087
2088 40 SM=0.
2089 JL=1
2090 DO 30 I=1,N-1
2091 S=2.*(P(I,4)*P(I+1,4)-P(I,1)*P(I+1,1)-P(I,2)*P(I+1,2)
2092 & -P(I,3)*P(I+1,3))+P(I,5)**2+P(I+1,5)**2
2093 IF(S.LT.0.) S=0.
2094 WP=SQRT(S)-1.5*(P(I,5)+P(I+1,5))
2095 IF(WP.GT.SM) THEN
2096 PBT1=P(I,1)+P(I+1,1)
2097 PBT2=P(I,2)+P(I+1,2)
2098 PBT3=P(I,3)+P(I+1,3)
2099 PBT4=P(I,4)+P(I+1,4)
2100 BTT=(PBT1**2+PBT2**2+PBT3**2)/PBT4**2
2101 IF(BTT.GE.1.0-1.0E-10) GO TO 30
2102 IF((I.NE.1.OR.I.NE.N-1).AND.
2103 & (K(I,2).NE.21.AND.K(I+1,2).NE.21)) GO TO 30
2104 JL=I
2105 SM=WP
2106 ENDIF
2107 30 CONTINUE
2108 S=(SM+1.5*(P(JL,5)+P(JL+1,5)))**2
2109 IF(SM.LT.HIPR1(5)) GOTO 2
2110
2111
2112 IF(JL+1.EQ.N) GOTO 190
2113 DO 160 J=N,JL+2,-1
2114 K(J+1,1)=K(J,1)
2115 K(J+1,2)=K(J,2)
2116 DO 150 M=1,5
2117 150 P(J+1,M)=P(J,M)
2118 160 CONTINUE
2119 190 N=N+1
2120
2121
2122 P1=P(JL,1)+P(JL+1,1)
2123 P2=P(JL,2)+P(JL+1,2)
2124 P3=P(JL,3)+P(JL+1,3)
2125 P4=P(JL,4)+P(JL+1,4)
2126 BEX=-P1/P4
2127 BEY=-P2/P4
2128 BEZ=-P3/P4
2129 IMIN=JL
2130 IMAX=JL+1
2131 CALL ATROBO(0.,0.,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
2132 IF(IERROR.NE.0) RETURN
2133
2134 CTH=P(JL,3)/SQRT(P(JL,4)**2-P(JL,5)**2)
2135 IF(ABS(CTH).GT.1.0) CTH=MAX(-1.,MIN(1.,CTH))
2136 THETA=ACOS(CTH)
2137 PHI=ULANGL(P(JL,1),P(JL,2))
2138 CALL ATROBO(0.,-PHI,0.,0.,0.,IMIN,IMAX,IERROR)
2139 CALL ATROBO(-THETA,0.,0.,0.,0.,IMIN,IMAX,IERROR)
2140
2141
2142
2143 1 CALL AR3JET(S,X1,X3,JL)
2144 CALL ARORIE(S,X1,X3,JL)
2145 IF(HIDAT(2).GT.0.0) THEN
2146 PTG1=SQRT(P(JL,1)**2+P(JL,2)**2)
2147 PTG2=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
2148 PTG3=SQRT(P(JL+2,1)**2+P(JL+2,2)**2)
2149 PTG=MAX(PTG1,PTG2,PTG3)
2150 IF(PTG.GT.HIDAT(2)) THEN
2151 FMFACT=EXP(-(PTG**2-HIDAT(2)**2)/HIPR1(2)**2)
2152 IF(RLU(NSEED).GT.FMFACT) GO TO 1
2153 ENDIF
2154 ENDIF
2155
2156 IMIN=JL
2157 IMAX=JL+2
2158 CALL ATROBO(THETA,PHI,-BEX,-BEY,-BEZ,IMIN,IMAX,IERROR)
2159 IF(IERROR.NE.0) RETURN
2160
2161 K(JL+2,1)=K(JL+1,1)
2162 K(JL+2,2)=K(JL+1,2)
2163 K(JL+2,3)=K(JL+1,3)
2164 K(JL+2,4)=K(JL+1,4)
2165 K(JL+2,5)=K(JL+1,5)
2166 P(JL+2,5)=P(JL+1,5)
2167 K(JL+1,1)=2
2168 K(JL+1,2)=21
2169 K(JL+1,3)=0
2170 K(JL+1,4)=0
2171 K(JL+1,5)=0
2172 P(JL+1,5)=0.
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184 IF(SM.GE.HIPR1(5)) GOTO 40
2185
2186 2 K(1,1)=2
2187 K(1,3)=0
2188 K(1,4)=0
2189 K(1,5)=0
2190 K(N,1)=1
2191 K(N,3)=0
2192 K(N,4)=0
2193 K(N,5)=0
2194
2195 RETURN
2196 END
2197
2198
2199 SUBROUTINE AR3JET(S,X1,X3,JL)
2200
2201 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2202 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
2203 COMMON/RANSEED/NSEED
2204 SAVE
2205
2206 C=1./3.
2207 IF(K(JL,2).NE.21 .AND. K(JL+1,2).NE.21) C=8./27.
2208 EXP1=3
2209 EXP3=3
2210 IF(K(JL,2).NE.21) EXP1=2
2211 IF(K(JL+1,2).NE.21) EXP3=2
2212 A=0.24**2/S
2213 YMA=ALOG(.5/SQRT(A)+SQRT(.25/A-1))
2214 D=4.*C*YMA
2215 SM1=P(JL,5)**2/S
2216 SM3=P(JL+1,5)**2/S
2217 XT2M=(1.-2.*SQRT(SM1)+SM1-SM3)*(1.-2.*SQRT(SM3)-SM1+SM3)
2218 XT2M=MIN(.25,XT2M)
2219 NTRY=0
2220 1 IF(NTRY.EQ.5000) THEN
2221 X1=.5*(2.*SQRT(SM1)+1.+SM1-SM3)
2222 X3=.5*(2.*SQRT(SM3)+1.-SM1+SM3)
2223 RETURN
2224 ENDIF
2225 NTRY=NTRY+1
2226
2227 XT2=A*(XT2M/A)**(RLU(NSEED)**(1./D))
2228
2229 YMAX=ALOG(.5/SQRT(XT2)+SQRT(.25/XT2-1.))
2230 Y=(2.*RLU(NSEED)-1.)*YMAX
2231 X1=1.-SQRT(XT2)*EXP(Y)
2232 X3=1.-SQRT(XT2)*EXP(-Y)
2233 X2=2.-X1-X3
2234 NEG=0
2235 IF(K(JL,2).NE.21 .OR. K(JL+1,2).NE.21) THEN
2236 IF((1.-X1)*(1.-X2)*(1.-X3)-X2*SM1*(1.-X1)-X2*SM3*(1.-X3).
2237 & LE.0..OR.X1.LE.2.*SQRT(SM1)-SM1+SM3.OR.X3.LE.2.*SQRT(SM3)
2238 & -SM3+SM1) NEG=1
2239 X1=X1+SM1-SM3
2240 X3=X3-SM1+SM3
2241 ENDIF
2242 IF(NEG.EQ.1) GOTO 1
2243
2244 FG=2.*YMAX*C*(X1**EXP1+X3**EXP3)/D
2245 XT2M=XT2
2246 IF(FG.LT.RLU(NSEED)) GOTO 1
2247
2248 RETURN
2249 END
2250
2251
2252
2253 SUBROUTINE ARORIE(S,X1,X3,JL)
2254
2255 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2256 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
2257 COMMON/RANSEED/NSEED
2258 SAVE
2259
2260 W=SQRT(S)
2261 X2=2.-X1-X3
2262 E1=.5*X1*W
2263 E3=.5*X3*W
2264 P1=SQRT(E1**2-P(JL,5)**2)
2265 P3=SQRT(E3**2-P(JL+1,5)**2)
2266 CBET=1.
2267 IF(P1.GT.0..AND.P3.GT.0.) CBET=(P(JL,5)**2
2268 & +P(JL+1,5)**2+2.*E1*E3-S*(1.-X2))/(2.*P1*P3)
2269 IF(ABS(CBET).GT.1.0) CBET=MAX(-1.,MIN(1.,CBET))
2270 BET=ACOS(CBET)
2271
2272
2273 IF(P1.GE.P3) THEN
2274 PSI=.5*ULANGL(P1**2+P3**2*COS(2.*BET),-P3**2*SIN(2.*BET))
2275 PT1=P1*SIN(PSI)
2276 PZ1=P1*COS(PSI)
2277 PT3=P3*SIN(PSI+BET)
2278 PZ3=P3*COS(PSI+BET)
2279 ELSE IF(P3.GT.P1) THEN
2280 PSI=.5*ULANGL(P3**2+P1**2*COS(2.*BET),-P1**2*SIN(2.*BET))
2281 PT1=P1*SIN(BET+PSI)
2282 PZ1=-P1*COS(BET+PSI)
2283 PT3=P3*SIN(PSI)
2284 PZ3=-P3*COS(PSI)
2285 ENDIF
2286
2287 DEL=2.0*HIPR1(40)*RLU(NSEED)
2288 P(JL,4)=E1
2289 P(JL,1)=PT1*SIN(DEL)
2290 P(JL,2)=-PT1*COS(DEL)
2291 P(JL,3)=PZ1
2292 P(JL+2,4)=E3
2293 P(JL+2,1)=PT3*SIN(DEL)
2294 P(JL+2,2)=-PT3*COS(DEL)
2295 P(JL+2,3)=PZ3
2296 P(JL+1,4)=W-E1-E3
2297 P(JL+1,1)=-P(JL,1)-P(JL+2,1)
2298 P(JL+1,2)=-P(JL,2)-P(JL+2,2)
2299 P(JL+1,3)=-P(JL,3)-P(JL+2,3)
2300 RETURN
2301 END
2302
2303
2304
2305
2306
2307
2308 SUBROUTINE ATROBO(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
2309 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
2310 DIMENSION ROT(3,3),PV(3)
2311 DOUBLE PRECISION DP(4),DBEX,DBEY,DBEZ,DGA,DGA2,DBEP,DGABEP
2312 SAVE
2313 IERROR=0
2314
2315 IF(IMIN.LE.0 .OR. IMAX.GT.N .OR. IMIN.GT.IMAX) RETURN
2316
2317 IF(THE**2+PHI**2.GT.1E-20) THEN
2318
2319 ROT(1,1)=COS(THE)*COS(PHI)
2320 ROT(1,2)=-SIN(PHI)
2321 ROT(1,3)=SIN(THE)*COS(PHI)
2322 ROT(2,1)=COS(THE)*SIN(PHI)
2323 ROT(2,2)=COS(PHI)
2324 ROT(2,3)=SIN(THE)*SIN(PHI)
2325 ROT(3,1)=-SIN(THE)
2326 ROT(3,2)=0.
2327 ROT(3,3)=COS(THE)
2328 DO 120 I=IMIN,IMAX
2329
2330 DO 100 J=1,3
2331 100 PV(J)=P(I,J)
2332 DO 110 J=1,3
2333 110 P(I,J)=ROT(J,1)*PV(1)+ROT(J,2)*PV(2)
2334 & +ROT(J,3)*PV(3)
2335 120 CONTINUE
2336 ENDIF
2337
2338 IF(BEX**2+BEY**2+BEZ**2.GT.1E-20) THEN
2339
2340 DBEX=BEX
2341 DBEY=BEY
2342 DBEZ=BEZ
2343 DGA2=1D0-DBEX**2-DBEY**2-DBEZ**2
2344 IF(DGA2.LE.0D0) THEN
2345 IERROR=1
2346 RETURN
2347 ENDIF
2348 DGA=1D0/DSQRT(DGA2)
2349 DO 140 I=IMIN,IMAX
2350
2351 DO 130 J=1,4
2352 130 DP(J)=P(I,J)
2353 DBEP=DBEX*DP(1)+DBEY*DP(2)+DBEZ*DP(3)
2354 DGABEP=DGA*(DGA*DBEP/(1D0+DGA)+DP(4))
2355 P(I,1)=DP(1)+DGABEP*DBEX
2356 P(I,2)=DP(2)+DGABEP*DBEY
2357 P(I,3)=DP(3)+DGABEP*DBEZ
2358 P(I,4)=DGA*(DP(4)+DBEP)
2359 140 CONTINUE
2360 ENDIF
2361
2362 RETURN
2363 END
2364
2365
2366
2367 SUBROUTINE HIJHRD(JP,JT,JOUT,JFLG,IOPJET0)
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379 DIMENSION IP(100,2),IPQ(50),IPB(50),IT(100,2),ITQ(50),ITB(50)
2380 COMMON/HIJCRDN/YP(3,300),YT(3,300)
2381 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2382 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
2383 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2384 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2385 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2386 & PJPM(300,500),NTJ(300),KFTJ(300,500),
2387 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2388 & PJTE(300,500),PJTM(300,500)
2389 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
2390 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
2391 & PZSG(900,100),PESG(900,100),PMSG(900,100)
2392 COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5),VDR(900,4)
2393 COMMON/RANSEED/NSEED
2394
2395 COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
2396 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2397 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
2398 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2399 COMMON/PYINT1/MINT(400),VINT(400)
2400 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
2401 COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
2402 COMMON/HIPYINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
2403 SAVE
2404
2405 MXJT=500
2406
2407 MXSG=900
2408
2409 MXSJ=100
2410
2411
2412 JFLG=0
2413 IHNT2(11)=JP
2414 IHNT2(12)=JT
2415
2416 IOPJET=IOPJET0
2417 IF(IOPJET.EQ.1.AND.(NFP(JP,6).NE.0.OR.NFT(JT,6).NE.0))
2418 & IOPJET=0
2419 IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
2420 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) RETURN
2421
2422
2423 IF(JOUT.EQ.0) THEN
2424 EPP=PP(JP,4)+PP(JP,3)
2425 EPM=PP(JP,4)-PP(JP,3)
2426 ETP=PT(JT,4)+PT(JT,3)
2427 ETM=PT(JT,4)-PT(JT,3)
2428 IF(EPP.LT.0.0) GO TO 1000
2429 IF(EPM.LT.0.0) GO TO 1000
2430 IF(ETP.LT.0.0) GO TO 1000
2431 IF(ETM.LT.0.0) GO TO 1000
2432 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
2433 ENDIF
2434
2435
2436
2437 ECUT1=HIPR1(1)+HIPR1(8)+PP(JP,14)+PP(JP,15)
2438 ECUT2=HIPR1(1)+HIPR1(8)+PT(JT,14)+PT(JT,15)
2439 IF(PP(JP,4).LE.ECUT1) THEN
2440 NFP(JP,6)=-ABS(NFP(JP,6))
2441 RETURN
2442 ENDIF
2443 IF(PT(JT,4).LE.ECUT2) THEN
2444 NFT(JT,6)=-ABS(NFT(JT,6))
2445 RETURN
2446 ENDIF
2447
2448
2449 MISS=0
2450 MISP=0
2451 MIST=0
2452
2453 IF(NFP(JP,10).EQ.0 .AND. NFT(JT,10).EQ.0) THEN
2454 MINT(44)=MINT4
2455 MINT(45)=MINT5
2456 XSEC(0,1)=ATXS(0)
2457 XSEC(11,1)=ATXS(11)
2458 XSEC(12,1)=ATXS(12)
2459 XSEC(28,1)=ATXS(28)
2460 DO 120 I=1,20
2461 COEF(11,I)=ATCO(11,I)
2462 COEF(12,I)=ATCO(12,I)
2463 COEF(28,I)=ATCO(28,I)
2464 120 CONTINUE
2465 ELSE
2466 ISUB11=0
2467 ISUB12=0
2468 ISUB28=0
2469 IF(XSEC(11,1).NE.0) ISUB11=1
2470 IF(XSEC(12,1).NE.0) ISUB12=1
2471 IF(XSEC(28,1).NE.0) ISUB28=1
2472 MINT(44)=MINT4-ISUB11-ISUB12-ISUB28
2473 MINT(45)=MINT5-ISUB11-ISUB12-ISUB28
2474 XSEC(0,1)=ATXS(0)-ATXS(11)-ATXS(12)-ATXS(28)
2475 XSEC(11,1)=0.0
2476 XSEC(12,1)=0.0
2477 XSEC(28,1)=0.0
2478 DO 110 I=1,20
2479 COEF(11,I)=0.0
2480 COEF(12,I)=0.0
2481 COEF(28,I)=0.0
2482 110 CONTINUE
2483 ENDIF
2484
2485
2486
2487 155 CALL PYTHIA
2488 JJ=MINT(31)
2489 IF(JJ.NE.1) GO TO 155
2490
2491 IF(K(7,2).EQ.-K(8,2)) THEN
2492 QMASS2=(P(7,4)+P(8,4))**2-(P(7,1)+P(8,1))**2
2493 & -(P(7,2)+P(8,2))**2-(P(7,3)+P(8,3))**2
2494 QM=ULMASS(K(7,2))
2495 IF(QMASS2.LT.(2.0*QM+HIPR1(1))**2) GO TO 155
2496 ENDIF
2497
2498 PXP=PP(JP,1)-P(3,1)
2499 PYP=PP(JP,2)-P(3,2)
2500 PZP=PP(JP,3)-P(3,3)
2501 PEP=PP(JP,4)-P(3,4)
2502 PXT=PT(JT,1)-P(4,1)
2503 PYT=PT(JT,2)-P(4,2)
2504 PZT=PT(JT,3)-P(4,3)
2505 PET=PT(JT,4)-P(4,4)
2506
2507 IF(PEP.LE.ECUT1) THEN
2508 MISP=MISP+1
2509 IF(MISP.LT.50) GO TO 155
2510 NFP(JP,6)=-ABS(NFP(JP,6))
2511 RETURN
2512 ENDIF
2513 IF(PET.LE.ECUT2) THEN
2514 MIST=MIST+1
2515 IF(MIST.LT.50) GO TO 155
2516 NFT(JT,6)=-ABS(NFT(JT,6))
2517 RETURN
2518 ENDIF
2519
2520
2521
2522 WP=PEP+PZP+PET+PZT
2523 WM=PEP-PZP+PET-PZT
2524 IF(WP.LT.0.0 .OR. WM.LT.0.0) THEN
2525 MISS=MISS+1
2526 IF(MISS.LT.50) GO TO 155
2527 RETURN
2528 ENDIF
2529
2530 SW=WP*WM
2531 AMPX=SQRT((ECUT1-HIPR1(8))**2+PXP**2+PYP**2+0.01)
2532 AMTX=SQRT((ECUT2-HIPR1(8))**2+PXT**2+PYT**2+0.01)
2533 SXX=(AMPX+AMTX)**2
2534 IF(SW.LT.SXX.OR.VINT(43).LT.HIPR1(1)) THEN
2535 MISS=MISS+1
2536 IF(MISS.LT.50) GO TO 155
2537 RETURN
2538 ENDIF
2539
2540
2541
2542
2543 HINT1(41)=P(7,1)
2544 HINT1(42)=P(7,2)
2545 HINT1(43)=P(7,3)
2546 HINT1(44)=P(7,4)
2547 HINT1(45)=P(7,5)
2548 HINT1(46)=SQRT(P(7,1)**2+P(7,2)**2)
2549 HINT1(51)=P(8,1)
2550 HINT1(52)=P(8,2)
2551 HINT1(53)=P(8,3)
2552 HINT1(54)=P(8,4)
2553 HINT1(55)=P(8,5)
2554 HINT1(56)=SQRT(P(8,1)**2+P(8,2)**2)
2555 IHNT2(14)=K(7,2)
2556 IHNT2(15)=K(8,2)
2557
2558 PINIRAD=(1.0-EXP(-2.0*(VINT(47)-HIDAT(1))))
2559 & /(1.0+EXP(-2.0*(VINT(47)-HIDAT(1))))
2560 I_INIRAD=0
2561 IF(RLU(NSEED).LE.PINIRAD) I_INIRAD=1
2562 IF(K(7,2).EQ.-K(8,2)) GO TO 190
2563 IF(K(7,2).EQ.21.AND.K(8,2).EQ.21.AND.IOPJET.EQ.1) GO TO 190
2564
2565
2566
2567
2568 JFLG=2
2569 JPP=0
2570 LPQ=0
2571 LPB=0
2572 JTT=0
2573 LTQ=0
2574 LTB=0
2575 IS7=0
2576 IS8=0
2577 HINT1(47)=0.0
2578 HINT1(48)=0.0
2579 HINT1(49)=0.0
2580 HINT1(50)=0.0
2581 HINT1(67)=0.0
2582 HINT1(68)=0.0
2583 HINT1(69)=0.0
2584 HINT1(70)=0.0
2585 DO 180 I=9,N
2586 IF(K(I,3).EQ.1 .OR. K(I,3).EQ.2.OR.
2587 & ABS(K(I,2)).GT.30) GO TO 180
2588
2589 IF(K(I,3).EQ.7) THEN
2590 HINT1(47)=HINT1(47)+P(I,1)
2591 HINT1(48)=HINT1(48)+P(I,2)
2592 HINT1(49)=HINT1(49)+P(I,3)
2593 HINT1(50)=HINT1(50)+P(I,4)
2594 ENDIF
2595 IF(K(I,3).EQ.8) THEN
2596 HINT1(67)=HINT1(67)+P(I,1)
2597 HINT1(68)=HINT1(68)+P(I,2)
2598 HINT1(69)=HINT1(69)+P(I,3)
2599 HINT1(70)=HINT1(70)+P(I,4)
2600 ENDIF
2601
2602 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
2603 NDR=NDR+1
2604 IADR(NDR,1)=JP
2605 IADR(NDR,2)=JT
2606 KFDR(NDR)=K(I,2)
2607 PDR(NDR,1)=P(I,1)
2608 PDR(NDR,2)=P(I,2)
2609 PDR(NDR,3)=P(I,3)
2610 PDR(NDR,4)=P(I,4)
2611 PDR(NDR,5)=P(I,5)
2612 VDR(NDR,1)=V(I,1)
2613 VDR(NDR,2)=V(I,2)
2614 VDR(NDR,3)=V(I,3)
2615 VDR(NDR,4)=V(I,4)
2616
2617 GO TO 180
2618
2619 ENDIF
2620 IF(K(I,3).EQ.7.OR.K(I,3).EQ.3) THEN
2621 IF(K(I,3).EQ.7.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(7,2)
2622 & .AND.IS7.EQ.0) THEN
2623 PP(JP,10)=P(I,1)
2624 PP(JP,11)=P(I,2)
2625 PP(JP,12)=P(I,3)
2626 PZP=PZP+P(I,3)
2627 PEP=PEP+P(I,4)
2628 NFP(JP,10)=1
2629 IS7=1
2630 GO TO 180
2631 ENDIF
2632 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
2633 & I_INIRAD.EQ.0)) THEN
2634 PXP=PXP+P(I,1)
2635 PYP=PYP+P(I,2)
2636 PZP=PZP+P(I,3)
2637 PEP=PEP+P(I,4)
2638 GO TO 180
2639 ENDIF
2640 JPP=JPP+1
2641 IP(JPP,1)=I
2642 IP(JPP,2)=0
2643 IF(K(I,2).NE.21) THEN
2644 IF(K(I,2).GT.0) THEN
2645 LPQ=LPQ+1
2646 IPQ(LPQ)=JPP
2647 IP(JPP,2)=LPQ
2648 ELSE IF(K(I,2).LT.0) THEN
2649 LPB=LPB+1
2650 IPB(LPB)=JPP
2651 IP(JPP,2)=-LPB
2652 ENDIF
2653 ENDIF
2654 ELSE IF(K(I,3).EQ.8.OR.K(I,3).EQ.4) THEN
2655 IF(K(I,3).EQ.8.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(8,2)
2656 & .AND.IS8.EQ.0) THEN
2657 PT(JT,10)=P(I,1)
2658 PT(JT,11)=P(I,2)
2659 PT(JT,12)=P(I,3)
2660 PZT=PZT+P(I,3)
2661 PET=PET+P(I,4)
2662 NFT(JT,10)=1
2663 IS8=1
2664 GO TO 180
2665 ENDIF
2666 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
2667 & I_INIRAD.EQ.0)) THEN
2668 PXT=PXT+P(I,1)
2669 PYT=PYT+P(I,2)
2670 PZT=PZT+P(I,3)
2671 PET=PET+P(I,4)
2672 GO TO 180
2673 ENDIF
2674 JTT=JTT+1
2675 IT(JTT,1)=I
2676 IT(JTT,2)=0
2677 IF(K(I,2).NE.21) THEN
2678 IF(K(I,2).GT.0) THEN
2679 LTQ=LTQ+1
2680 ITQ(LTQ)=JTT
2681 IT(JTT,2)=LTQ
2682 ELSE IF(K(I,2).LT.0) THEN
2683 LTB=LTB+1
2684 ITB(LTB)=JTT
2685 IT(JTT,2)=-LTB
2686 ENDIF
2687 ENDIF
2688 ENDIF
2689 180 CONTINUE
2690
2691
2692 IF(LPQ.NE.LPB .OR. LTQ.NE.LTB) THEN
2693 MISS=MISS+1
2694 IF(MISS.LE.50) GO TO 155
2695 WRITE(6,*) ' Q -QBAR NOT MATCHED IN HIJHRD'
2696 JFLG=0
2697 RETURN
2698 ENDIF
2699
2700
2701
2702 J=0
2703 181 J=J+1
2704 IF(J.GT.JPP) GO TO 182
2705 IF(IP(J,2).EQ.0) THEN
2706 GO TO 181
2707 ELSE IF(IP(J,2).NE.0) THEN
2708 LP=ABS(IP(J,2))
2709 IP1=IP(J,1)
2710 IP2=IP(J,2)
2711 IP(J,1)=IP(IPQ(LP),1)
2712 IP(J,2)=IP(IPQ(LP),2)
2713 IP(IPQ(LP),1)=IP1
2714 IP(IPQ(LP),2)=IP2
2715 IF(IP2.GT.0) THEN
2716 IPQ(IP2)=IPQ(LP)
2717 ELSE IF(IP2.LT.0) THEN
2718 IPB(-IP2)=IPQ(LP)
2719 ENDIF
2720
2721 IP1=IP(J+1,1)
2722 IP2=IP(J+1,2)
2723 IP(J+1,1)=IP(IPB(LP),1)
2724 IP(J+1,2)=IP(IPB(LP),2)
2725 IP(IPB(LP),1)=IP1
2726 IP(IPB(LP),2)=IP2
2727 IF(IP2.GT.0) THEN
2728 IPQ(IP2)=IPB(LP)
2729 ELSE IF(IP2.LT.0) THEN
2730 IPB(-IP2)=IPB(LP)
2731 ENDIF
2732
2733 J=J+1
2734 GO TO 181
2735 ENDIF
2736
2737 182 J=0
2738 183 J=J+1
2739 IF(J.GT.JTT) GO TO 184
2740 IF(IT(J,2).EQ.0) THEN
2741 GO TO 183
2742 ELSE IF(IT(J,2).NE.0) THEN
2743 LT=ABS(IT(J,2))
2744 IT1=IT(J,1)
2745 IT2=IT(J,2)
2746 IT(J,1)=IT(ITQ(LT),1)
2747 IT(J,2)=IT(ITQ(LT),2)
2748 IT(ITQ(LT),1)=IT1
2749 IT(ITQ(LT),2)=IT2
2750 IF(IT2.GT.0) THEN
2751 ITQ(IT2)=ITQ(LT)
2752 ELSE IF(IT2.LT.0) THEN
2753 ITB(-IT2)=ITQ(LT)
2754 ENDIF
2755
2756 IT1=IT(J+1,1)
2757 IT2=IT(J+1,2)
2758 IT(J+1,1)=IT(ITB(LT),1)
2759 IT(J+1,2)=IT(ITB(LT),2)
2760 IT(ITB(LT),1)=IT1
2761 IT(ITB(LT),2)=IT2
2762 IF(IT2.GT.0) THEN
2763 ITQ(IT2)=ITB(LT)
2764 ELSE IF(IT2.LT.0) THEN
2765 ITB(-IT2)=ITB(LT)
2766 ENDIF
2767
2768 J=J+1
2769 GO TO 183
2770
2771 ENDIF
2772
2773 184 CONTINUE
2774 IF(NPJ(JP)+JPP.GT.MXJT.OR.NTJ(JT)+JTT.GT.MXJT) THEN
2775 JFLG=0
2776 WRITE(6,*) 'number of partons per string exceeds'
2777 WRITE(6,*) 'the common block size'
2778 RETURN
2779 ENDIF
2780
2781 DO 186 J=1,JPP
2782 KFPJ(JP,NPJ(JP)+J)=K(IP(J,1),2)
2783 PJPX(JP,NPJ(JP)+J)=P(IP(J,1),1)
2784 PJPY(JP,NPJ(JP)+J)=P(IP(J,1),2)
2785 PJPZ(JP,NPJ(JP)+J)=P(IP(J,1),3)
2786 PJPE(JP,NPJ(JP)+J)=P(IP(J,1),4)
2787 PJPM(JP,NPJ(JP)+J)=P(IP(J,1),5)
2788 186 CONTINUE
2789 NPJ(JP)=NPJ(JP)+JPP
2790 DO 188 J=1,JTT
2791 KFTJ(JT,NTJ(JT)+J)=K(IT(J,1),2)
2792 PJTX(JT,NTJ(JT)+J)=P(IT(J,1),1)
2793 PJTY(JT,NTJ(JT)+J)=P(IT(J,1),2)
2794 PJTZ(JT,NTJ(JT)+J)=P(IT(J,1),3)
2795 PJTE(JT,NTJ(JT)+J)=P(IT(J,1),4)
2796 PJTM(JT,NTJ(JT)+J)=P(IT(J,1),5)
2797 188 CONTINUE
2798 NTJ(JT)=NTJ(JT)+JTT
2799 GO TO 900
2800
2801
2802
2803 190 JFLG=3
2804 IF(K(7,2).NE.21.AND.K(8,2).NE.21.AND.
2805 & K(7,2)*K(8,2).GT.0) GO TO 155
2806 JPP=0
2807 LPQ=0
2808 LPB=0
2809 DO 200 I=9,N
2810 IF(K(I,3).EQ.1.OR.K(I,3).EQ.2.OR.
2811 & ABS(K(I,2)).GT.30) GO TO 200
2812 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
2813 NDR=NDR+1
2814 IADR(NDR,1)=JP
2815 IADR(NDR,2)=JT
2816 KFDR(NDR)=K(I,2)
2817 PDR(NDR,1)=P(I,1)
2818 PDR(NDR,2)=P(I,2)
2819 PDR(NDR,3)=P(I,3)
2820 PDR(NDR,4)=P(I,4)
2821 PDR(NDR,5)=P(I,5)
2822 VDR(NDR,1)=V(I,1)
2823 VDR(NDR,2)=V(I,2)
2824 VDR(NDR,3)=V(I,3)
2825 VDR(NDR,4)=V(I,4)
2826
2827 GO TO 200
2828
2829 ENDIF
2830 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
2831 & I_INIRAD.EQ.0)) THEN
2832 PXP=PXP+P(I,1)
2833 PYP=PYP+P(I,2)
2834 PZP=PZP+P(I,3)
2835 PEP=PEP+P(I,4)
2836 GO TO 200
2837 ENDIF
2838 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
2839 & I_INIRAD.EQ.0)) THEN
2840 PXT=PXT+P(I,1)
2841 PYT=PYT+P(I,2)
2842 PZT=PZT+P(I,3)
2843 PET=PET+P(I,4)
2844 GO TO 200
2845 ENDIF
2846 JPP=JPP+1
2847 IP(JPP,1)=I
2848 IP(JPP,2)=0
2849 IF(K(I,2).NE.21) THEN
2850 IF(K(I,2).GT.0) THEN
2851 LPQ=LPQ+1
2852 IPQ(LPQ)=JPP
2853 IP(JPP,2)=LPQ
2854 ELSE IF(K(I,2).LT.0) THEN
2855 LPB=LPB+1
2856 IPB(LPB)=JPP
2857 IP(JPP,2)=-LPB
2858 ENDIF
2859 ENDIF
2860 200 CONTINUE
2861 IF(LPQ.NE.LPB) THEN
2862 MISS=MISS+1
2863 IF(MISS.LE.50) GO TO 155
2864 WRITE(6,*) LPQ,LPB, 'Q-QBAR NOT CONSERVED OR NOT MATCHED'
2865 JFLG=0
2866 RETURN
2867 ENDIF
2868
2869
2870
2871 J=0
2872 220 J=J+1
2873 IF(J.GT.JPP) GO TO 222
2874 IF(IP(J,2).EQ.0) GO TO 220
2875 LP=ABS(IP(J,2))
2876 IP1=IP(J,1)
2877 IP2=IP(J,2)
2878 IP(J,1)=IP(IPQ(LP),1)
2879 IP(J,2)=IP(IPQ(LP),2)
2880 IP(IPQ(LP),1)=IP1
2881 IP(IPQ(LP),2)=IP2
2882 IF(IP2.GT.0) THEN
2883 IPQ(IP2)=IPQ(LP)
2884 ELSE IF(IP2.LT.0) THEN
2885 IPB(-IP2)=IPQ(LP)
2886 ENDIF
2887 IPQ(LP)=J
2888
2889 IP1=IP(J+1,1)
2890 IP2=IP(J+1,2)
2891 IP(J+1,1)=IP(IPB(LP),1)
2892 IP(J+1,2)=IP(IPB(LP),2)
2893 IP(IPB(LP),1)=IP1
2894 IP(IPB(LP),2)=IP2
2895 IF(IP2.GT.0) THEN
2896 IPQ(IP2)=IPB(LP)
2897 ELSE IF(IP2.LT.0) THEN
2898 IPB(-IP2)=IPB(LP)
2899 ENDIF
2900
2901 IPB(LP)=J+1
2902 J=J+1
2903 GO TO 220
2904
2905 222 CONTINUE
2906 IF(LPQ.GE.1) THEN
2907 DO 240 L0=2,LPQ
2908 IP1=IP(2*L0-3,1)
2909 IP2=IP(2*L0-3,2)
2910 IP(2*L0-3,1)=IP(IPQ(L0),1)
2911 IP(2*L0-3,2)=IP(IPQ(L0),2)
2912 IP(IPQ(L0),1)=IP1
2913 IP(IPQ(L0),2)=IP2
2914 IF(IP2.GT.0) THEN
2915 IPQ(IP2)=IPQ(L0)
2916 ELSE IF(IP2.LT.0) THEN
2917 IPB(-IP2)=IPQ(L0)
2918 ENDIF
2919 IPQ(L0)=2*L0-3
2920
2921 IP1=IP(2*L0-2,1)
2922 IP2=IP(2*L0-2,2)
2923 IP(2*L0-2,1)=IP(IPB(L0),1)
2924 IP(2*L0-2,2)=IP(IPB(L0),2)
2925 IP(IPB(L0),1)=IP1
2926 IP(IPB(L0),2)=IP2
2927 IF(IP2.GT.0) THEN
2928 IPQ(IP2)=IPB(L0)
2929 ELSE IF(IP2.LT.0) THEN
2930 IPB(-IP2)=IPB(L0)
2931 ENDIF
2932 IPB(L0)=2*L0-2
2933 240 CONTINUE
2934
2935
2936 IP1=IP(2*LPQ-1,1)
2937 IP2=IP(2*LPQ-1,2)
2938 IP(2*LPQ-1,1)=IP(IPQ(1),1)
2939 IP(2*LPQ-1,2)=IP(IPQ(1),2)
2940 IP(IPQ(1),1)=IP1
2941 IP(IPQ(1),2)=IP2
2942 IF(IP2.GT.0) THEN
2943 IPQ(IP2)=IPQ(1)
2944 ELSE IF(IP2.LT.0) THEN
2945 IPB(-IP2)=IPQ(1)
2946 ENDIF
2947 IPQ(1)=2*LPQ-1
2948
2949
2950 IP1=IP(JPP,1)
2951 IP2=IP(JPP,2)
2952 IP(JPP,1)=IP(IPB(1),1)
2953 IP(JPP,2)=IP(IPB(1),2)
2954 IP(IPB(1),1)=IP1
2955 IP(IPB(1),2)=IP2
2956 IF(IP2.GT.0) THEN
2957 IPQ(IP2)=IPB(1)
2958 ELSE IF(IP2.LT.0) THEN
2959 IPB(-IP2)=IPB(1)
2960 ENDIF
2961 IPB(1)=JPP
2962
2963
2964 ENDIF
2965 IF(NSG.GE.MXSG) THEN
2966 JFLG=0
2967 WRITE(6,*) 'number of jets forming single strings exceeds'
2968 WRITE(6,*) 'the common block size'
2969 RETURN
2970 ENDIF
2971 IF(JPP.GT.MXSJ) THEN
2972 JFLG=0
2973 WRITE(6,*) 'number of partons per single jet system'
2974 WRITE(6,*) 'exceeds the common block size'
2975 RETURN
2976 ENDIF
2977
2978 NSG=NSG+1
2979 NJSG(NSG)=JPP
2980 IASG(NSG,1)=JP
2981 IASG(NSG,2)=JT
2982 IASG(NSG,3)=0
2983 DO 300 I=1,JPP
2984 K1SG(NSG,I)=2
2985 K2SG(NSG,I)=K(IP(I,1),2)
2986 IF(K2SG(NSG,I).LT.0) K1SG(NSG,I)=1
2987 PXSG(NSG,I)=P(IP(I,1),1)
2988 PYSG(NSG,I)=P(IP(I,1),2)
2989 PZSG(NSG,I)=P(IP(I,1),3)
2990 PESG(NSG,I)=P(IP(I,1),4)
2991 PMSG(NSG,I)=P(IP(I,1),5)
2992 300 CONTINUE
2993 K1SG(NSG,1)=2
2994 K1SG(NSG,JPP)=1
2995
2996 900 PP(JP,1)=PXP
2997 PP(JP,2)=PYP
2998 PP(JP,3)=PZP
2999 PP(JP,4)=PEP
3000 PP(JP,5)=0.0
3001 PT(JT,1)=PXT
3002 PT(JT,2)=PYT
3003 PT(JT,3)=PZT
3004 PT(JT,4)=PET
3005 PT(JT,5)=0.0
3006
3007 NFP(JP,6)=NFP(JP,6)+1
3008 NFT(JT,6)=NFT(JT,6)+1
3009 RETURN
3010
3011 1000 JFLG=-1
3012 IF(IHPR2(10).EQ.0) RETURN
3013 WRITE(6,*) 'Fatal HIJHRD error'
3014 WRITE(6,*) JP, ' proj E+,E-',EPP,EPM,' status',NFP(JP,5)
3015 WRITE(6,*) JT, ' targ E+,E_',ETP,ETM,' status',NFT(JT,5)
3016 RETURN
3017 END
3018
3019
3020
3021
3022
3023 SUBROUTINE JETINI(JP,JT,I_TRIG)
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038 CHARACTER BEAM*16,TARG*16
3039 DIMENSION XSEC0(8,0:200),COEF0(8,200,20),INI(8),
3040 & MINT44(8),MINT45(8)
3041 COMMON/HIJCRDN/YP(3,300),YT(3,300)
3042 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3043 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3044 COMMON/HIPYINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
3045
3046 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3047 COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
3048 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
3049 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
3050 COMMON/PYINT1/MINT(400),VINT(400)
3051 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
3052 COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
3053 SAVE
3054 DATA INI/8*0/I_LAST/-1/
3055
3056 IHNT2(11)=JP
3057 IHNT2(12)=JT
3058 IF(IHNT2(5).NE.0 .AND. IHNT2(6).NE.0) THEN
3059 I_TYPE=1
3060 ELSE IF(IHNT2(5).NE.0 .AND. IHNT2(6).EQ.0) THEN
3061 I_TYPE=1
3062 IF(NFT(JT,4).EQ.2112) I_TYPE=2
3063 ELSE IF(IHNT2(5).EQ.0 .AND. IHNT2(6).NE.0) THEN
3064 I_TYPE=1
3065 IF(NFP(JP,4).EQ.2112) I_TYPE=2
3066 ELSE
3067 IF(NFP(JP,4).EQ.2212 .AND. NFT(JT,4).EQ.2212) THEN
3068 I_TYPE=1
3069 ELSE IF(NFP(JP,4).EQ.2212 .AND. NFT(JT,4).EQ.2112) THEN
3070 I_TYPE=2
3071 ELSE IF(NFP(JP,4).EQ.2112 .AND. NFT(JT,4).EQ.2212) THEN
3072 I_TYPE=3
3073 ELSE
3074 I_TYPE=4
3075 ENDIF
3076 ENDIF
3077
3078 IF(I_TRIG.NE.0) GO TO 160
3079 IF(I_TRIG.EQ.I_LAST) GO TO 150
3080 MSTP(2)=2
3081
3082 MSTP(33)=1
3083 PARP(31)=HIPR1(17)
3084
3085 MSTP(51)=3
3086
3087 MSTP(61)=1
3088
3089 MSTP(71)=1
3090
3091 IF(IHPR2(2).EQ.0.OR.IHPR2(2).EQ.2) MSTP(61)=0
3092 IF(IHPR2(2).EQ.0.OR.IHPR2(2).EQ.1) MSTP(71)=0
3093
3094 MSTP(81)=0
3095
3096 MSTP(82)=1
3097
3098 MSTP(111)=0
3099
3100 IF(IHPR2(10).EQ.0) MSTP(122)=0
3101
3102 PARP(81)=HIPR1(8)
3103 CKIN(5)=HIPR1(8)
3104 CKIN(3)=HIPR1(8)
3105 CKIN(4)=HIPR1(9)
3106 IF(HIPR1(9).LE.HIPR1(8)) CKIN(4)=-1.0
3107 CKIN(9)=-10.0
3108 CKIN(10)=10.0
3109 MSEL=0
3110 DO 100 ISUB=1,200
3111 MSUB(ISUB)=0
3112 100 CONTINUE
3113 MSUB(11)=1
3114 MSUB(12)=1
3115 MSUB(13)=1
3116 MSUB(28)=1
3117 MSUB(53)=1
3118 MSUB(68)=1
3119 MSUB(81)=1
3120 MSUB(82)=1
3121 DO 110 J=1,MIN(8,MDCY(21,3))
3122 110 MDME(MDCY(21,2)+J-1,1)=0
3123 ISEL=4
3124 IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
3125 MDME(MDCY(21,2)+ISEL-1,1)=1
3126
3127 MSUB(14)=1
3128 MSUB(18)=1
3129 MSUB(29)=1
3130
3131 150 IF(INI(I_TYPE).NE.0) GO TO 800
3132 GO TO 400
3133
3134
3135
3136 160 I_TYPE=4+I_TYPE
3137 IF(I_TRIG.EQ.I_LAST) GO TO 260
3138 PARP(81)=ABS(HIPR1(10))-0.25
3139 CKIN(5)=ABS(HIPR1(10))-0.25
3140 CKIN(3)=ABS(HIPR1(10))-0.25
3141 CKIN(4)=ABS(HIPR1(10))+0.25
3142 IF(HIPR1(10).LT.HIPR1(8)) CKIN(4)=-1.0
3143
3144 MSEL=0
3145 DO 101 ISUB=1,200
3146 MSUB(ISUB)=0
3147 101 CONTINUE
3148 IF(IHPR2(3).EQ.1) THEN
3149 MSUB(11)=1
3150 MSUB(12)=1
3151 MSUB(13)=1
3152 MSUB(28)=1
3153 MSUB(53)=1
3154 MSUB(68)=1
3155 MSUB(81)=1
3156 MSUB(82)=1
3157 MSUB(14)=1
3158 MSUB(18)=1
3159 MSUB(29)=1
3160 DO 102 J=1,MIN(8,MDCY(21,3))
3161 102 MDME(MDCY(21,2)+J-1,1)=0
3162 ISEL=4
3163 IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
3164 MDME(MDCY(21,2)+ISEL-1,1)=1
3165
3166 ELSE IF(IHPR2(3).EQ.2) THEN
3167 MSUB(14)=1
3168 MSUB(18)=1
3169 MSUB(29)=1
3170
3171
3172 ELSE IF(IHPR2(3).EQ.3) THEN
3173 CKIN(3)=MAX(0.0,HIPR1(10))
3174 CKIN(5)=HIPR1(8)
3175 PARP(81)=HIPR1(8)
3176 MSUB(81)=1
3177 MSUB(82)=1
3178 DO 105 J=1,MIN(8,MDCY(21,3))
3179 105 MDME(MDCY(21,2)+J-1,1)=0
3180 ISEL=4
3181 IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
3182 MDME(MDCY(21,2)+ISEL-1,1)=1
3183
3184 ENDIF
3185 260 IF(INI(I_TYPE).NE.0) GO TO 800
3186
3187
3188 400 INI(I_TYPE)=1
3189 IF(IHPR2(10).EQ.0) MSTP(122)=0
3190 IF(NFP(JP,4).EQ.2212) THEN
3191 BEAM='P'
3192 ELSE IF(NFP(JP,4).EQ.-2212) THEN
3193 BEAM='P~'
3194 ELSE IF(NFP(JP,4).EQ.2112) THEN
3195 BEAM='N'
3196 ELSE IF(NFP(JP,4).EQ.-2112) THEN
3197 BEAM='N~'
3198 ELSE IF(NFP(JP,4).EQ.211) THEN
3199 BEAM='PI+'
3200 ELSE IF(NFP(JP,4).EQ.-211) THEN
3201 BEAM='PI-'
3202 ELSE IF(NFP(JP,4).EQ.321) THEN
3203 BEAM='PI+'
3204 ELSE IF(NFP(JP,4).EQ.-321) THEN
3205 BEAM='PI-'
3206 ELSE
3207 WRITE(6,*) 'unavailable beam type', NFP(JP,4)
3208 ENDIF
3209 IF(NFT(JT,4).EQ.2212) THEN
3210 TARG='P'
3211 ELSE IF(NFT(JT,4).EQ.-2212) THEN
3212 TARG='P~'
3213 ELSE IF(NFT(JT,4).EQ.2112) THEN
3214 TARG='N'
3215 ELSE IF(NFT(JT,4).EQ.-2112) THEN
3216 TARG='N~'
3217 ELSE IF(NFT(JT,4).EQ.211) THEN
3218 TARG='PI+'
3219 ELSE IF(NFT(JT,4).EQ.-211) THEN
3220 TARG='PI-'
3221 ELSE IF(NFT(JT,4).EQ.321) THEN
3222 TARG='PI+'
3223 ELSE IF(NFT(JT,4).EQ.-321) THEN
3224 TARG='PI-'
3225 ELSE
3226 WRITE(6,*) 'unavailable target type', NFT(JT,4)
3227 ENDIF
3228
3229 IHNT2(16)=1
3230
3231
3232
3233 CALL PYINIT('CMS',BEAM,TARG,HINT1(1))
3234 MINT4=MINT(44)
3235 MINT5=MINT(45)
3236 MINT44(I_TYPE)=MINT(44)
3237 MINT45(I_TYPE)=MINT(45)
3238 ATXS(0)=XSEC(0,1)
3239 XSEC0(I_TYPE,0)=XSEC(0,1)
3240 DO 500 I=1,200
3241 ATXS(I)=XSEC(I,1)
3242 XSEC0(I_TYPE,I)=XSEC(I,1)
3243 DO 500 J=1,20
3244 ATCO(I,J)=COEF(I,J)
3245 COEF0(I_TYPE,I,J)=COEF(I,J)
3246 500 CONTINUE
3247
3248 I_LAST=I_TRIG
3249 IHNT2(16)=0
3250
3251 RETURN
3252
3253
3254
3255
3256 800 MINT(44)=MINT44(I_TYPE)
3257 MINT(45)=MINT45(I_TYPE)
3258 MINT4=MINT(44)
3259 MINT5=MINT(45)
3260 XSEC(0,1)=XSEC0(I_TYPE,0)
3261 ATXS(0)=XSEC(0,1)
3262 DO 900 I=1,200
3263 XSEC(I,1)=XSEC0(I_TYPE,I)
3264 ATXS(I)=XSEC(I,1)
3265 DO 900 J=1,20
3266 COEF(I,J)=COEF0(I_TYPE,I,J)
3267 ATCO(I,J)=COEF(I,J)
3268 900 CONTINUE
3269 I_LAST=I_TRIG
3270 MINT(11)=NFP(JP,4)
3271 MINT(12)=NFT(JT,4)
3272 RETURN
3273 END
3274
3275
3276
3277 SUBROUTINE HIJINI
3278 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3279 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3280 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3281 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3282 & PJPM(300,500),NTJ(300),KFTJ(300,500),
3283 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3284 & PJTE(300,500),PJTM(300,500)
3285 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
3286 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
3287 & PZSG(900,100),PESG(900,100),PMSG(900,100)
3288 COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5),VDR(900,4)
3289 COMMON/RANSEED/NSEED
3290 SAVE
3291
3292
3293
3294 NSG=0
3295 NDR=0
3296 IPP=2212
3297 IPT=2212
3298 IF(IHNT2(5).NE.0) IPP=IHNT2(5)
3299 IF(IHNT2(6).NE.0) IPT=IHNT2(6)
3300
3301
3302 DO 100 I=1,IHNT2(1)
3303 PP(I,1)=0.0
3304 PP(I,2)=0.0
3305 PP(I,3)=SQRT(HINT1(1)**2/4.0-HINT1(8)**2)
3306 PP(I,4)=HINT1(1)/2
3307 PP(I,5)=HINT1(8)
3308 PP(I,6)=0.0
3309 PP(I,7)=0.0
3310 PP(I,8)=0.0
3311 PP(I,9)=0.0
3312 PP(I,10)=0.0
3313 NFP(I,3)=IPP
3314 NFP(I,4)=IPP
3315 NFP(I,5)=0
3316 NFP(I,6)=0
3317 NFP(I,7)=0
3318 NFP(I,8)=0
3319 NFP(I,9)=0
3320 NFP(I,10)=0
3321 NFP(I,11)=0
3322 NPJ(I)=0
3323 IF(I.GT.ABS(IHNT2(2))) NFP(I,3)=2112
3324 CALL ATTFLV(NFP(I,3),IDQ,IDQQ)
3325 NFP(I,1)=IDQ
3326 NFP(I,2)=IDQQ
3327 NFP(I,15)=-1
3328 IF(ABS(IDQ).GT.1000.OR.(ABS(IDQ*IDQQ).LT.100.AND.
3329 & RLU(NSEED).LT.0.5)) NFP(I,15)=1
3330 PP(I,14)=ULMASS(IDQ)
3331 PP(I,15)=ULMASS(IDQQ)
3332 100 CONTINUE
3333
3334 DO 200 I=1,IHNT2(3)
3335 PT(I,1)=0.0
3336 PT(I,2)=0.0
3337 PT(I,3)=-SQRT(HINT1(1)**2/4.0-HINT1(9)**2)
3338 PT(I,4)=HINT1(1)/2.0
3339 PT(I,5)=HINT1(9)
3340 PT(I,6)=0.0
3341 PT(I,7)=0.0
3342 PT(I,8)=0.0
3343 PT(I,9)=0.0
3344 PT(I,10)=0.0
3345 NFT(I,3)=IPT
3346 NFT(I,4)=IPT
3347 NFT(I,5)=0
3348 NFT(I,6)=0
3349 NFT(I,7)=0
3350 NFT(I,8)=0
3351 NFT(I,9)=0
3352 NFT(I,10)=0
3353 NFT(I,11)=0
3354 NTJ(I)=0
3355 IF(I.GT.ABS(IHNT2(4))) NFT(I,3)=2112
3356 CALL ATTFLV(NFT(I,3),IDQ,IDQQ)
3357 NFT(I,1)=IDQ
3358 NFT(I,2)=IDQQ
3359 NFT(I,15)=1
3360 IF(ABS(IDQ).GT.1000.OR.(ABS(IDQ*IDQQ).LT.100.AND.
3361 & RLU(NSEED).LT.0.5)) NFT(I,15)=-1
3362 PT(I,14)=ULMASS(IDQ)
3363 PT(I,15)=ULMASS(IDQQ)
3364 200 CONTINUE
3365 RETURN
3366 END
3367
3368
3369
3370 SUBROUTINE ATTFLV(ID,IDQ,IDQQ)
3371 COMMON/RANSEED/NSEED
3372 SAVE
3373
3374 IF(ABS(ID).LT.100) THEN
3375 NSIGN=1
3376 IDQ=ID/100
3377 IDQQ=-ID/10+IDQ*10
3378 IF(ABS(IDQ).EQ.3) NSIGN=-1
3379 IDQ=NSIGN*IDQ
3380 IDQQ=NSIGN*IDQQ
3381 IF(IDQ.LT.0) THEN
3382 ID0=IDQ
3383 IDQ=IDQQ
3384 IDQQ=ID0
3385 ENDIF
3386 RETURN
3387 ENDIF
3388
3389
3390
3391
3392
3393
3394
3395
3396
3397
3398 IDQ=2
3399 IF(ABS(ID).EQ.2112) IDQ=1
3400 IDQQ=2101
3401 X=RLU(NSEED)
3402 IF(X.LE.0.5) GO TO 30
3403 IF(X.GT.0.666667) GO TO 10
3404 IDQQ=2103
3405 GO TO 30
3406 10 IDQ=1
3407 IDQQ=2203
3408 IF(ABS(ID).EQ.2112) THEN
3409 IDQ=2
3410 IDQQ=1103
3411 ENDIF
3412 30 IF(ID.LT.0) THEN
3413 ID00=IDQQ
3414 IDQQ=-IDQ
3415 IDQ=-ID00
3416 ENDIF
3417 RETURN
3418 END
3419
3420
3421
3422
3423
3424 SUBROUTINE HIJCSC(JP,JT)
3425 DIMENSION PSC1(5),PSC2(5)
3426 COMMON/HIJCRDN/YP(3,300),YT(3,300)
3427 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3428 COMMON/RANSEED/NSEED
3429 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3430 SAVE
3431 IF(JP.EQ.0 .OR. JT.EQ.0) GO TO 25
3432 DO 10 I=1,5
3433 PSC1(I)=PP(JP,I)
3434 PSC2(I)=PT(JT,I)
3435 10 CONTINUE
3436 CALL HIJELS(PSC1,PSC2)
3437 DPP1=PSC1(1)-PP(JP,1)
3438 DPP2=PSC1(2)-PP(JP,2)
3439 DPT1=PSC2(1)-PT(JT,1)
3440 DPT2=PSC2(2)-PT(JT,2)
3441 PP(JP,6)=PP(JP,6)+DPP1/2.0
3442 PP(JP,7)=PP(JP,7)+DPP2/2.0
3443 PP(JP,8)=PP(JP,8)+DPP1/2.0
3444 PP(JP,9)=PP(JP,9)+DPP2/2.0
3445 PT(JT,6)=PT(JT,6)+DPT1/2.0
3446 PT(JT,7)=PT(JT,7)+DPT2/2.0
3447 PT(JT,8)=PT(JT,8)+DPT1/2.0
3448 PT(JT,9)=PT(JT,9)+DPT2/2.0
3449 DO 20 I=1,4
3450 PP(JP,I)=PSC1(I)
3451 PT(JT,I)=PSC2(I)
3452 20 CONTINUE
3453 NFP(JP,5)=MAX(1,NFP(JP,5))
3454 NFT(JT,5)=MAX(1,NFT(JT,5))
3455
3456 RETURN
3457
3458
3459 25 IF(JP.EQ.0) GO TO 45
3460 PABS=SQRT(PP(JP,1)**2+PP(JP,2)**2+PP(JP,3)**2)
3461 BX=PP(JP,1)/PABS
3462 BY=PP(JP,2)/PABS
3463 BZ=PP(JP,3)/PABS
3464 DO 40 I=1,IHNT2(1)
3465 IF(I.EQ.JP) GO TO 40
3466 DX=YP(1,I)-YP(1,JP)
3467 DY=YP(2,I)-YP(2,JP)
3468 DZ=YP(3,I)-YP(3,JP)
3469 DIS=DX*BX+DY*BY+DZ*BZ
3470 IF(DIS.LE.0) GO TO 40
3471 BB=DX**2+DY**2+DZ**2-DIS**2
3472 R2=BB*HIPR1(40)/HIPR1(31)/0.1
3473
3474 GS=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
3475 & *ROMG(R2))**2
3476 GS0=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
3477 & *ROMG(0.0))**2
3478 IF(RLU(NSEED).GT.GS/GS0) GO TO 40
3479 DO 30 K=1,5
3480 PSC1(K)=PP(JP,K)
3481 PSC2(K)=PP(I,K)
3482 30 CONTINUE
3483 CALL HIJELS(PSC1,PSC2)
3484 DPP1=PSC1(1)-PP(JP,1)
3485 DPP2=PSC1(2)-PP(JP,2)
3486 DPT1=PSC2(1)-PP(I,1)
3487 DPT2=PSC2(2)-PP(I,2)
3488 PP(JP,6)=PP(JP,6)+DPP1/2.0
3489 PP(JP,7)=PP(JP,7)+DPP2/2.0
3490 PP(JP,8)=PP(JP,8)+DPP1/2.0
3491 PP(JP,9)=PP(JP,9)+DPP2/2.0
3492 PP(I,6)=PP(I,6)+DPT1/2.0
3493 PP(I,7)=PP(I,7)+DPT2/2.0
3494 PP(I,8)=PP(I,8)+DPT1/2.0
3495 PP(I,9)=PP(I,9)+DPT2/2.0
3496 DO 35 K=1,5
3497 PP(JP,K)=PSC1(K)
3498 PP(I,K)=PSC2(K)
3499 35 CONTINUE
3500 NFP(I,5)=MAX(1,NFP(I,5))
3501 GO TO 45
3502 40 CONTINUE
3503 45 IF(JT.EQ.0) GO TO 80
3504 PABS=SQRT(PT(JT,1)**2+PT(JT,2)**2+PT(JT,3)**2)
3505 BX=PT(JT,1)/PABS
3506 BY=PT(JT,2)/PABS
3507 BZ=PT(JT,3)/PABS
3508 DO 70 I=1,IHNT2(3)
3509 IF(I.EQ.JT) GO TO 70
3510 DX=YT(1,I)-YT(1,JT)
3511 DY=YT(2,I)-YT(2,JT)
3512 DZ=YT(3,I)-YT(3,JT)
3513 DIS=DX*BX+DY*BY+DZ*BZ
3514 IF(DIS.LE.0) GO TO 70
3515 BB=DX**2+DY**2+DZ**2-DIS**2
3516 R2=BB*HIPR1(40)/HIPR1(31)/0.1
3517
3518 GS=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
3519 & *ROMG(R2)))**2
3520 GS0=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
3521 & *ROMG(0.0)))**2
3522 IF(RLU(NSEED).GT.GS/GS0) GO TO 70
3523 DO 60 K=1,5
3524 PSC1(K)=PT(JT,K)
3525 PSC2(K)=PT(I,K)
3526 60 CONTINUE
3527 CALL HIJELS(PSC1,PSC2)
3528 DPP1=PSC1(1)-PT(JT,1)
3529 DPP2=PSC1(2)-PT(JT,2)
3530 DPT1=PSC2(1)-PT(I,1)
3531 DPT2=PSC2(2)-PT(I,2)
3532 PT(JT,6)=PT(JT,6)+DPP1/2.0
3533 PT(JT,7)=PT(JT,7)+DPP2/2.0
3534 PT(JT,8)=PT(JT,8)+DPP1/2.0
3535 PT(JT,9)=PT(JT,9)+DPP2/2.0
3536 PT(I,6)=PT(I,6)+DPT1/2.0
3537 PT(I,7)=PT(I,7)+DPT2/2.0
3538 PT(I,8)=PT(I,8)+DPT1/2.0
3539 PT(I,9)=PT(I,9)+DPT2/2.0
3540 DO 65 K=1,5
3541 PT(JT,K)=PSC1(K)
3542 PT(I,K)=PSC2(K)
3543 65 CONTINUE
3544 NFT(I,5)=MAX(1,NFT(I,5))
3545 GO TO 80
3546 70 CONTINUE
3547 80 RETURN
3548 END
3549
3550
3551
3552
3553
3554
3555 SUBROUTINE HIJELS(PSC1,PSC2)
3556 IMPLICIT DOUBLE PRECISION(D)
3557 DIMENSION PSC1(5),PSC2(5)
3558 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3559 COMMON/RANSEED/NSEED
3560 SAVE
3561
3562 CC=1.0-HINT1(12)/HINT1(13)
3563 RR=(1.0-CC)*HINT1(13)/HINT1(12)/(1.0-HIPR1(33))-1.0
3564 BB=0.5*(3.0+RR+SQRT(9.0+10.0*RR+RR**2))
3565 EP=SQRT((PSC1(1)-PSC2(1))**2+(PSC1(2)-PSC2(2))**2
3566 & +(PSC1(3)-PSC2(3))**2)
3567 IF(EP.LE.0.1) RETURN
3568 ELS0=98.0/EP+52.0*(1.0+RR)**2
3569 PCM1=PSC1(1)+PSC2(1)
3570 PCM2=PSC1(2)+PSC2(2)
3571 PCM3=PSC1(3)+PSC2(3)
3572 ECM=PSC1(4)+PSC2(4)
3573 AM1=PSC1(5)**2
3574 AM2=PSC2(5)**2
3575 AMM=ECM**2-PCM1**2-PCM2**2-PCM3**2
3576 IF(AMM.LE.PSC1(5)+PSC2(5)) RETURN
3577
3578
3579 PMAX=(AMM**2+AM1**2+AM2**2-2.0*AMM*AM1-2.0*AMM*AM2
3580 & -2.0*AM1*AM2)/4.0/AMM
3581 PMAX=ABS(PMAX)
3582 20 TT=RLU(NSEED)*MIN(PMAX,1.5)
3583 ELS=98.0*EXP(-2.8*TT)/EP
3584 & +52.0*EXP(-9.2*TT)*(1.0+RR*EXP(-4.6*(BB-1.0)*TT))**2
3585 IF(RLU(NSEED).GT.ELS/ELS0) GO TO 20
3586 PHI=2.0*HIPR1(40)*RLU(NSEED)
3587
3588 DBX=PCM1/ECM
3589 DBY=PCM2/ECM
3590 DBZ=PCM3/ECM
3591 DB=SQRT(DBX**2+DBY**2+DBZ**2)
3592 IF(DB.GT.0.99999999D0) THEN
3593 DBX=DBX*(0.99999999D0/DB)
3594 DBY=DBY*(0.99999999D0/DB)
3595 DBZ=DBZ*(0.99999999D0/DB)
3596 DB=0.99999999D0
3597 WRITE(6,*) ' (HIJELS) boost vector too large'
3598
3599 ENDIF
3600 DGA=1D0/SQRT(1D0-DB**2)
3601
3602 DP1=SQRT(TT)*SIN(PHI)
3603 DP2=SQRT(TT)*COS(PHI)
3604 DP3=SQRT(PMAX-TT)
3605 DP4=SQRT(PMAX+AM1)
3606 DBP=DBX*DP1+DBY*DP2+DBZ*DP3
3607 DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4)
3608 PSC1(1)=DP1+DGABP*DBX
3609 PSC1(2)=DP2+DGABP*DBY
3610 PSC1(3)=DP3+DGABP*DBZ
3611 PSC1(4)=DGA*(DP4+DBP)
3612
3613 DP1=-SQRT(TT)*SIN(PHI)
3614 DP2=-SQRT(TT)*COS(PHI)
3615 DP3=-SQRT(PMAX-TT)
3616 DP4=SQRT(PMAX+AM2)
3617 DBP=DBX*DP1+DBY*DP2+DBZ*DP3
3618 DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4)
3619 PSC2(1)=DP1+DGABP*DBX
3620 PSC2(2)=DP2+DGABP*DBY
3621 PSC2(3)=DP3+DGABP*DBZ
3622 PSC2(4)=DGA*(DP4+DBP)
3623 RETURN
3624 END
3625
3626
3627
3628
3629
3630
3631
3632
3633 SUBROUTINE HIJSFT(JP,JT,JOUT,IERROR)
3634 COMMON/HIJCRDN/YP(3,300),YT(3,300)
3635 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3636 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3637 COMMON/RANSEED/NSEED
3638 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3639 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3640 & PJPM(300,500),NTJ(300),KFTJ(300,500),
3641 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3642 & PJTE(300,500),PJTM(300,500)
3643 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
3644 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
3645 & PZSG(900,100),PESG(900,100),PMSG(900,100)
3646 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3647 COMMON/DPMCOM1/JJP,JJT,AMP,AMT,APX0,ATX0,AMPN,AMTN,AMP0,AMT0,
3648 & NFDP,NFDT,WP,WM,SW,XREMP,XREMT,DPKC1,DPKC2,PP11,PP12,
3649 & PT11,PT12,PTP2,PTT2
3650 COMMON/DPMCOM2/NDPM,KDPM(20,2),PDPM1(20,5),PDPM2(20,5)
3651 SAVE
3652
3653
3654
3655
3656
3657
3658 IERROR=0
3659 JJP=JP
3660 JJT=JT
3661 NDPM=0
3662 IOPMAIN=0
3663 IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
3664 I_SNG=0
3665
3666 EPP=PP(JP,4)+PP(JP,3)
3667 EPM=PP(JP,4)-PP(JP,3)
3668 ETP=PT(JT,4)+PT(JT,3)
3669 ETM=PT(JT,4)-PT(JT,3)
3670
3671 WP=EPP+ETP
3672 WM=EPM+ETM
3673 SW=WP*WM
3674
3675
3676 IF(WP.LT.0.0 .OR. WM.LT.0.0) GO TO 1000
3677
3678 IF(JOUT.EQ.0) THEN
3679 IF(EPP.LT.0.0) GO TO 1000
3680 IF(EPM.LT.0.0) GO TO 1000
3681 IF(ETP.LT.0.0) GO TO 1000
3682 IF(ETM.LT.0.0) GO TO 1000
3683 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
3684 ENDIF
3685
3686
3687
3688
3689
3690 IHNT2(11)=JP
3691 IHNT2(12)=JT
3692
3693
3694
3695 MISS=0
3696 PKC1=0.0
3697 PKC2=0.0
3698 PKC11=0.0
3699 PKC12=0.0
3700 PKC21=0.0
3701 PKC22=0.0
3702 DPKC11=0.0
3703 DPKC12=0.0
3704 DPKC21=0.0
3705 DPKC22=0.0
3706 IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
3707 IF(NFP(JP,10).EQ.1) THEN
3708 PHI1=ULANGL(PP(JP,10),PP(JP,11))
3709 PPJET=SQRT(PP(JP,10)**2+PP(JP,11)**2)
3710 PKC1=PPJET
3711 PKC11=PP(JP,10)
3712 PKC12=PP(JP,11)
3713 ENDIF
3714 IF(NFT(JT,10).EQ.1) THEN
3715 PHI2=ULANGL(PT(JT,10),PT(JT,11))
3716 PTJET=SQRT(PT(JT,10)**2+PT(JT,11)**2)
3717 PKC2=PTJET
3718 PKC21=PT(JT,10)
3719 PKC22=PT(JT,11)
3720 ENDIF
3721 IF(IHPR2(4).GT.0.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
3722 IF(NFP(JP,10).EQ.0) THEN
3723 PHI=-PHI2
3724 ELSE IF(NFT(JT,10).EQ.0) THEN
3725 PHI=PHI1
3726 ELSE
3727 PHI=(PHI1+PHI2-HIPR1(40))/2.0
3728 ENDIF
3729 BX=HINT1(19)*COS(HINT1(20))
3730 BY=HINT1(19)*SIN(HINT1(20))
3731 XP0=YP(1,JP)
3732 YP0=YP(2,JP)
3733 XT0=YT(1,JT)+BX
3734 YT0=YT(2,JT)+BY
3735 R1=MAX(1.2*IHNT2(1)**0.3333333,
3736 & SQRT(XP0**2+YP0**2))
3737 R2=MAX(1.2*IHNT2(3)**0.3333333,
3738 & SQRT((XT0-BX)**2+(YT0-BY)**2))
3739 IF(ABS(COS(PHI)).LT.1.0E-5) THEN
3740 DD1=R1
3741 DD2=R1
3742 DD3=ABS(BY+SQRT(R2**2-(XP0-BX)**2)-YP0)
3743 DD4=ABS(BY-SQRT(R2**2-(XP0-BX)**2)-YP0)
3744 GO TO 5
3745 ENDIF
3746 BB=2.0*SIN(PHI)*(COS(PHI)*YP0-SIN(PHI)*XP0)
3747 CC=(YP0**2-R1**2)*COS(PHI)**2+XP0*SIN(PHI)*(
3748 & XP0*SIN(PHI)-2.0*YP0*COS(PHI))
3749 DD=BB**2-4.0*CC
3750 IF(DD.LT.0.0) GO TO 10
3751 XX1=(-BB+SQRT(DD))/2.0
3752 XX2=(-BB-SQRT(DD))/2.0
3753 DD1=ABS((XX1-XP0)/COS(PHI))
3754 DD2=ABS((XX2-XP0)/COS(PHI))
3755
3756 BB=2.0*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)-2.0*BX
3757 CC=(BX**2+(YT0-BY)**2-R2**2)*COS(PHI)**2+XT0*SIN(PHI)
3758 & *(XT0*SIN(PHI)-2.0*COS(PHI)*(YT0-BY))
3759 & -2.0*BX*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)
3760 DD=BB**2-4.0*CC
3761 IF(DD.LT.0.0) GO TO 10
3762 XX1=(-BB+SQRT(DD))/2.0
3763 XX2=(-BB-SQRT(DD))/2.0
3764 DD3=ABS((XX1-XT0)/COS(PHI))
3765 DD4=ABS((XX2-XT0)/COS(PHI))
3766
3767 5 DD1=MIN(DD1,DD3)
3768 DD2=MIN(DD2,DD4)
3769 IF(DD1.LT.HIPR1(13)) DD1=0.0
3770 IF(DD2.LT.HIPR1(13)) DD2=0.0
3771 IF(NFP(JP,10).EQ.1.AND.PPJET.GT.HIPR1(11)) THEN
3772 DP1=DD1*HIPR1(14)/2.0
3773 DP1=MIN(DP1,PPJET-HIPR1(11))
3774 PKC1=PPJET-DP1
3775 DPX1=COS(PHI1)*DP1
3776 DPY1=SIN(PHI1)*DP1
3777 PKC11=PP(JP,10)-DPX1
3778 PKC12=PP(JP,11)-DPY1
3779 IF(DP1.GT.0.0) THEN
3780 CTHEP=PP(JP,12)/SQRT(PP(JP,12)**2+PPJET**2)
3781 DPZ1=DP1*CTHEP/SQRT(1.0-CTHEP**2)
3782 DPE1=SQRT(DPX1**2+DPY1**2+DPZ1**2)
3783 EPPPRM=PP(JP,4)+PP(JP,3)-DPE1-DPZ1
3784 EPMPRM=PP(JP,4)-PP(JP,3)-DPE1+DPZ1
3785 IF(EPPPRM.LE.0.0.OR.EPMPRM.LE.0.0) GO TO 15
3786 EPP=EPPPRM
3787 EPM=EPMPRM
3788 PP(JP,10)=PKC11
3789 PP(JP,11)=PKC12
3790 NPJ(JP)=NPJ(JP)+1
3791 KFPJ(JP,NPJ(JP))=21
3792 PJPX(JP,NPJ(JP))=DPX1
3793 PJPY(JP,NPJ(JP))=DPY1
3794 PJPZ(JP,NPJ(JP))=DPZ1
3795 PJPE(JP,NPJ(JP))=DPE1
3796 PJPM(JP,NPJ(JP))=0.0
3797 PP(JP,3)=PP(JP,3)-DPZ1
3798 PP(JP,4)=PP(JP,4)-DPE1
3799 ENDIF
3800 ENDIF
3801 15 IF(NFT(JT,10).EQ.1.AND.PTJET.GT.HIPR1(11)) THEN
3802 DP2=DD2*HIPR1(14)/2.0
3803 DP2=MIN(DP2,PTJET-HIPR1(11))
3804 PKC2=PTJET-DP2
3805 DPX2=COS(PHI2)*DP2
3806 DPY2=SIN(PHI2)*DP2
3807 PKC21=PT(JT,10)-DPX2
3808 PKC22=PT(JT,11)-DPY2
3809 IF(DP2.GT.0.0) THEN
3810 CTHET=PT(JT,12)/SQRT(PT(JT,12)**2+PTJET**2)
3811 DPZ2=DP2*CTHET/SQRT(1.0-CTHET**2)
3812 DPE2=SQRT(DPX2**2+DPY2**2+DPZ2**2)
3813 ETPPRM=PT(JT,4)+PT(JT,3)-DPE2-DPZ2
3814 ETMPRM=PT(JT,4)-PT(JT,3)-DPE2+DPZ2
3815 IF(ETPPRM.LE.0.0.OR.ETMPRM.LE.0.0) GO TO 16
3816 ETP=ETPPRM
3817 ETM=ETMPRM
3818 PT(JT,10)=PKC21
3819 PT(JT,11)=PKC22
3820 NTJ(JT)=NTJ(JT)+1
3821 KFTJ(JT,NTJ(JT))=21
3822 PJTX(JT,NTJ(JT))=DPX2
3823 PJTY(JT,NTJ(JT))=DPY2
3824 PJTZ(JT,NTJ(JT))=DPZ2
3825 PJTE(JT,NTJ(JT))=DPE2
3826 PJTM(JT,NTJ(JT))=0.0
3827 PT(JT,3)=PT(JT,3)-DPZ2
3828 PT(JT,4)=PT(JT,4)-DPE2
3829 ENDIF
3830 ENDIF
3831 16 DPKC11=-(PP(JP,10)-PKC11)/2.0
3832 DPKC12=-(PP(JP,11)-PKC12)/2.0
3833 DPKC21=-(PT(JT,10)-PKC21)/2.0
3834 DPKC22=-(PT(JT,11)-PKC22)/2.0
3835 WP=EPP+ETP
3836 WM=EPM+ETM
3837 SW=WP*WM
3838 ENDIF
3839 ENDIF
3840
3841
3842
3843
3844 10 PTP02=PP(JP,1)**2+PP(JP,2)**2
3845 PTT02=PT(JT,1)**2+PT(JT,2)**2
3846
3847 AMQ=MAX(PP(JP,14)+PP(JP,15),PT(JT,14)+PT(JT,15))
3848 AMX=HIPR1(1)+AMQ
3849
3850
3851 AMP0=AMX
3852 DPM0=AMX
3853 NFDP=0
3854 IF(NFP(JP,5).LE.2.AND.NFP(JP,3).NE.0) THEN
3855 AMP0=ULMASS(NFP(JP,3))
3856 NFDP=NFP(JP,3)+2*NFP(JP,3)/ABS(NFP(JP,3))
3857 DPM0=ULMASS(NFDP)
3858 IF(DPM0.LE.0.0) THEN
3859 NFDP=NFDP-2*NFDP/ABS(NFDP)
3860 DPM0=ULMASS(NFDP)
3861 ENDIF
3862 ENDIF
3863 AMT0=AMX
3864 DTM0=AMX
3865 NFDT=0
3866 IF(NFT(JT,5).LE.2.AND.NFT(JT,3).NE.0) THEN
3867 AMT0=ULMASS(NFT(JT,3))
3868 NFDT=NFT(JT,3)+2*NFT(JT,3)/ABS(NFT(JT,3))
3869 DTM0=ULMASS(NFDT)
3870 IF(DTM0.LE.0.0) THEN
3871 NFDT=NFDT-2*NFDT/ABS(NFDT)
3872 DTM0=ULMASS(NFDT)
3873 ENDIF
3874 ENDIF
3875
3876 AMPN=SQRT(AMP0**2+PTP02)
3877 AMTN=SQRT(AMT0**2+PTT02)
3878 SNN=(AMPN+AMTN)**2+0.001
3879
3880 IF(SW.LT.SNN+0.001) GO TO 4000
3881
3882
3883 SWPTN=4.0*(MAX(AMP0,AMT0)**2+MAX(PTP02,PTT02))
3884 SWPTD=4.0*(MAX(DPM0,DTM0)**2+MAX(PTP02,PTT02))
3885 SWPTX=4.0*(AMX**2+MAX(PTP02,PTT02))
3886 IF(SW.LE.SWPTN) THEN
3887 PKCMX=0.0
3888 ELSE IF(SW.GT.SWPTN .AND. SW.LE.SWPTD
3889 & .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
3890 PKCMX=SQRT(SW/4.0-MAX(AMP0,AMT0)**2)
3891 & -SQRT(MAX(PTP02,PTT02))
3892 ELSE IF(SW.GT.SWPTD .AND. SW.LE.SWPTX
3893 & .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
3894 PKCMX=SQRT(SW/4.0-MAX(DPM0,DTM0)**2)
3895 & -SQRT(MAX(PTP02,PTT02))
3896 ELSE IF(SW.GT.SWPTX) THEN
3897 PKCMX=SQRT(SW/4.0-AMX**2)-SQRT(MAX(PTP02,PTT02))
3898 ENDIF
3899
3900
3901
3902 IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
3903 IF(PKC1.GT.PKCMX) THEN
3904 PKC1=PKCMX
3905 PKC11=PKC1*COS(PHI1)
3906 PKC12=PKC1*SIN(PHI1)
3907 DPKC11=-(PP(JP,10)-PKC11)/2.0
3908 DPKC12=-(PP(JP,11)-PKC12)/2.0
3909 ENDIF
3910 IF(PKC2.GT.PKCMX) THEN
3911 PKC2=PKCMX
3912 PKC21=PKC2*COS(PHI2)
3913 PKC22=PKC2*SIN(PHI2)
3914 DPKC21=-(PT(JT,10)-PKC21)/2.0
3915 DPKC22=-(PT(JT,11)-PKC22)/2.0
3916 ENDIF
3917 DPKC1=DPKC11+DPKC21
3918 DPKC2=DPKC12+DPKC22
3919 NFP(JP,10)=-NFP(JP,10)
3920 NFT(JT,10)=-NFT(JT,10)
3921 GO TO 40
3922 ENDIF
3923
3924
3925 I_SNG=0
3926 IF(IHPR2(13).NE.0 .AND. RLU(NSEED).LE.HIDAT(4)) I_SNG=1
3927 IF((NFP(JP,5).EQ.3 .OR.NFT(JT,5).EQ.3).OR.
3928 & (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
3929 & (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) I_SNG=0
3930
3931
3932 IF(IHPR2(5).EQ.0) THEN
3933 PKC=HIPR1(2)*SQRT(-ALOG(1.0-RLU(NSEED)
3934 & *(1.0-EXP(-PKCMX**2/HIPR1(2)**2))))
3935 GO TO 30
3936 ENDIF
3937 PKC=HIRND2(3,0.0,PKCMX**2)
3938 PKC=SQRT(PKC)
3939 IF(PKC.GT.HIPR1(20))
3940 & PKC=HIPR1(2)*SQRT(-ALOG(EXP(-HIPR1(20)**2/HIPR1(2)**2)
3941 & -RLU(NSEED)*(EXP(-HIPR1(20)**2/HIPR1(2)**2)-
3942 & EXP(-PKCMX**2/HIPR1(2)**2))))
3943
3944 IF(I_SNG.EQ.1) PKC=0.65*SQRT(
3945 & -ALOG(1.0-RLU(NSEED)*(1.0-EXP(-PKCMX**2/0.65**2))))
3946
3947 30 PHI0=2.0*HIPR1(40)*RLU(NSEED)
3948 PKC11=PKC*SIN(PHI0)
3949 PKC12=PKC*COS(PHI0)
3950 PKC21=-PKC11
3951 PKC22=-PKC12
3952 DPKC1=0.0
3953 DPKC2=0.0
3954 40 PP11=PP(JP,1)+PKC11-DPKC1
3955 PP12=PP(JP,2)+PKC12-DPKC2
3956 PT11=PT(JT,1)+PKC21-DPKC1
3957 PT12=PT(JT,2)+PKC22-DPKC2
3958 PTP2=PP11**2+PP12**2
3959 PTT2=PT11**2+PT12**2
3960
3961 AMPN=SQRT(AMP0**2+PTP2)
3962 AMTN=SQRT(AMT0**2+PTT2)
3963 SNN=(AMPN+AMTN)**2+0.001
3964
3965 WP=EPP+ETP
3966 WM=EPM+ETM
3967 SW=WP*WM
3968
3969 IF(SW.LT.SNN) THEN
3970 MISS=MISS+1
3971 IF(MISS.LE.100) then
3972 PKC=0.0
3973 GO TO 30
3974 ENDIF
3975 IF(IHPR2(10).NE.0)
3976 & WRITE(6,*) 'Error occured in Pt kick section of HIJSFT'
3977 GO TO 4000
3978 ENDIF
3979
3980 AMPD=SQRT(DPM0**2+PTP2)
3981 AMTD=SQRT(DTM0**2+PTT2)
3982
3983 AMPX=SQRT(AMX**2+PTP2)
3984 AMTX=SQRT(AMX**2+PTT2)
3985
3986 DPN=AMPN**2/SW
3987 DTN=AMTN**2/SW
3988 DPD=AMPD**2/SW
3989 DTD=AMTD**2/SW
3990 DPX=AMPX**2/SW
3991 DTX=AMTX**2/SW
3992
3993 SPNTD=(AMPN+AMTD)**2
3994 SPNTX=(AMPN+AMTX)**2
3995
3996 SPDTN=(AMPD+AMTN)**2
3997 SPXTN=(AMPX+AMTN)**2
3998
3999 SPDTX=(AMPD+AMTX)**2
4000 SPXTD=(AMPX+AMTD)**2
4001 SDD=(AMPD+AMTD)**2
4002 SXX=(AMPX+AMTX)**2
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012 CONTINUE
4013 IF(SW.GT.SXX+0.001) THEN
4014 IF(I_SNG.EQ.0) THEN
4015 D1=DPX
4016 D2=DTX
4017 NFP3=0
4018 NFT3=0
4019 GO TO 400
4020 ELSE
4021
4022
4023 IF((NFP(JP,5).EQ.3 .AND.NFT(JT,5).EQ.3).OR.
4024 & (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
4025 & (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
4026 D1=DPX
4027 D2=DTX
4028 NFP3=0
4029 NFT3=0
4030 GO TO 400
4031 ENDIF
4032
4033
4034 IF(RLU(NSEED).GT.0.5.OR.(NFT(JT,5).GT.2.OR.
4035 & NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
4036 D1=DPN
4037 D2=DTX
4038 NFP3=NFP(JP,3)
4039 NFT3=0
4040 GO TO 220
4041 ELSE
4042 D1=DPX
4043 D2=DTN
4044 NFP3=0
4045 NFT3=NFT(JT,3)
4046 GO TO 240
4047 ENDIF
4048
4049 ENDIF
4050 ELSE IF(SW.GT.MAX(SPDTX,SPXTD)+0.001 .AND.
4051 & SW.LE.SXX+0.001) THEN
4052 IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0.AND.
4053 & RLU(NSEED).GT.0.5).OR.(NPJ(JP).EQ.0
4054 & .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
4055 D1=DPD
4056 D2=DTX
4057 NFP3=NFDP
4058 NFT3=0
4059 GO TO 220
4060 ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
4061 D1=DPX
4062 D2=DTD
4063 NFP3=0
4064 NFT3=NFDT
4065 GO TO 240
4066 ENDIF
4067 GO TO 4000
4068 ELSE IF(SW.GT.MIN(SPDTX,SPXTD)+0.001.AND.
4069 & SW.LE.MAX(SPDTX,SPXTD)+0.001) THEN
4070 IF(SPDTX.LE.SPXTD.AND.NPJ(JP).EQ.0
4071 & .AND.NFP(JP,5).LE.2) THEN
4072 D1=DPD
4073 D2=DTX
4074 NFP3=NFDP
4075 NFT3=0
4076 GO TO 220
4077 ELSE IF(SPDTX.GT.SPXTD.AND.NTJ(JT).EQ.0
4078 & .AND.NFT(JT,5).LE.2) THEN
4079 D1=DPX
4080 D2=DTD
4081 NFP3=0
4082 NFT3=NFDT
4083 GO TO 240
4084 ENDIF
4085
4086
4087 IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
4088 & .AND.RLU(NSEED).GT.0.5).OR.(NPJ(JP).EQ.0
4089 & .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
4090 D1=DPN
4091 D2=DTX
4092 NFP3=NFP(JP,3)
4093 NFT3=0
4094 GO TO 220
4095 ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
4096 D1=DPX
4097 D2=DTN
4098 NFP3=0
4099 NFT3=NFT(JT,3)
4100 GO TO 240
4101 ENDIF
4102 GO TO 4000
4103 ELSE IF(SW.GT.MAX(SPNTX,SPXTN)+0.001 .AND.
4104 & SW.LE.MIN(SPDTX,SPXTD)+0.001) THEN
4105 IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
4106 & .AND.RLU(NSEED).GT.0.5).OR.(NPJ(JP).EQ.0
4107 & .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
4108 D1=DPN
4109 D2=DTX
4110 NFP3=NFP(JP,3)
4111 NFT3=0
4112 GO TO 220
4113 ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
4114 D1=DPX
4115 D2=DTN
4116 NFP3=0
4117 NFT3=NFT(JT,3)
4118 GO TO 240
4119 ENDIF
4120 GO TO 4000
4121 ELSE IF(SW.GT.MIN(SPNTX,SPXTN)+0.001 .AND.
4122 & SW.LE.MAX(SPNTX,SPXTN)+0.001) THEN
4123 IF(SPNTX.LE.SPXTN.AND.NPJ(JP).EQ.0
4124 & .AND.NFP(JP,5).LE.2) THEN
4125 D1=DPN
4126 D2=DTX
4127 NFP3=NFP(JP,3)
4128 NFT3=0
4129 GO TO 220
4130 ELSEIF(SPNTX.GT.SPXTN.AND.NTJ(JT).EQ.0
4131 & .AND.NFT(JT,5).LE.2) THEN
4132 D1=DPX
4133 D2=DTN
4134 NFP3=0
4135 NFT3=NFT(JT,3)
4136 GO TO 240
4137 ENDIF
4138 GO TO 4000
4139 ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
4140 & (NPJ(JP).NE.0 .OR.NTJ(JT).NE.0)) THEN
4141 GO TO 4000
4142 ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
4143 & NFP(JP,5).GT.2.AND.NFT(JT,5).GT.2) THEN
4144 GO TO 4000
4145 ELSE IF(SW.GT.SDD+0.001.AND.SW.LE.
4146 & MIN(SPNTX,SPXTN)+0.001) THEN
4147 D1=DPD
4148 D2=DTD
4149 NFP3=NFDP
4150 NFT3=NFDT
4151 GO TO 100
4152 ELSE IF(SW.GT.MAX(SPNTD,SPDTN)+0.001
4153 & .AND. SW.LE.SDD+0.001) THEN
4154 IF(RLU(NSEED).GT.0.5) THEN
4155 D1=DPD
4156 D2=DTN
4157 NFP3=NFDP
4158 NFT3=NFT(JT,3)
4159 GO TO 100
4160 ELSE
4161 D1=DPN
4162 D2=DTD
4163 NFP3=NFP(JP,3)
4164 NFT3=NFDT
4165 GO TO 100
4166 ENDIF
4167 ELSE IF(SW.GT.MIN(SPNTD,SPDTN)+0.001
4168 & .AND. SW.LE.MAX(SPNTD,SPDTN)+0.001) THEN
4169 IF(SPNTD.GT.SPDTN) THEN
4170 D1=DPD
4171 D2=DTN
4172 NFP3=NFDP
4173 NFT3=NFT(JT,3)
4174 GO TO 100
4175 ELSE
4176 D1=DPN
4177 D2=DTD
4178 NFP3=NFP(JP,3)
4179 NFT3=NFDT
4180 GO TO 100
4181 ENDIF
4182 ELSE IF(SW.LE.MIN(SPNTD,SPDTN)+0.001) THEN
4183 D1=DPN
4184 D2=DTN
4185 NFP3=NFP(JP,3)
4186 NFT3=NFT(JT,3)
4187 GO TO 100
4188 ENDIF
4189 WRITE(6,*) ' Error in HIJSFT: There is no path to here'
4190 RETURN
4191
4192
4193
4194
4195
4196 100 NFP5=MAX(2,NFP(JP,5))
4197 NFT5=MAX(2,NFT(JT,5))
4198 BB1=1.0+D1-D2
4199 BB2=1.0+D2-D1
4200 IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
4201 MISS=MISS+1
4202 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
4203 PKC=PKC*0.5
4204 GO TO 30
4205 ENDIF
4206 IF(RLU(NSEED).LT.0.5) THEN
4207 X1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
4208 X2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
4209 ELSE
4210 X1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
4211 X2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
4212 ENDIF
4213 IHNT2(13)=2
4214 GO TO 600
4215
4216
4217
4218
4219 220 NFP5=MAX(2,NFP(JP,5))
4220 NFT5=3
4221 IF(NFP3.EQ.0) NFP5=3
4222 BB2=1.0+D2-D1
4223 IF(BB2**2.LT.4.0*D2) THEN
4224 MISS=MISS+1
4225 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
4226 PKC=PKC*0.5
4227 GO TO 30
4228 ENDIF
4229 XMIN=(BB2-SQRT(BB2**2-4.0*D2))/2.0
4230 XMAX=(BB2+SQRT(BB2**2-4.0*D2))/2.0
4231 MISS4=0
4232 222 X2=HIRND2(6,XMIN,XMAX)
4233 X1=D1/(1.0-X2)
4234 IF(X2*(1.0-X1).LT.(D2+1.E-4/SW)) THEN
4235 MISS4=MISS4+1
4236 IF(MISS4.LE.1000) GO TO 222
4237 GO TO 5000
4238 ENDIF
4239 IHNT2(13)=2
4240 GO TO 600
4241
4242 240 NFP5=3
4243 NFT5=MAX(2,NFT(JT,5))
4244 IF(NFT3.EQ.0) NFT5=3
4245 BB1=1.0+D1-D2
4246 IF(BB1**2.LT.4.0*D1) THEN
4247 MISS=MISS+1
4248 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
4249 PKC=PKC*0.5
4250 GO TO 30
4251 ENDIF
4252 XMIN=(BB1-SQRT(BB1**2-4.0*D1))/2.0
4253 XMAX=(BB1+SQRT(BB1**2-4.0*D1))/2.0
4254 MISS4=0
4255 242 X1=HIRND2(6,XMIN,XMAX)
4256 X2=D2/(1.0-X1)
4257 IF(X1*(1.0-X2).LT.(D1+1.E-4/SW)) THEN
4258 MISS4=MISS4+1
4259 IF(MISS4.LE.1000) GO TO 242
4260 GO TO 5000
4261 ENDIF
4262 IHNT2(13)=2
4263 GO TO 600
4264
4265
4266
4267
4268
4269
4270 400 NFP5=3
4271 NFT5=3
4272 BB1=1.0+D1-D2
4273 BB2=1.0+D2-D1
4274 IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
4275 MISS=MISS+1
4276 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
4277 PKC=PKC*0.5
4278 GO TO 30
4279 ENDIF
4280 XMIN1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
4281 XMAX1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
4282 XMIN2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
4283 XMAX2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
4284 MISS4=0
4285 410 X1=HIRND2(4,XMIN1,XMAX1)
4286 X2=HIRND2(4,XMIN2,XMAX2)
4287 IF(NFP(JP,5).EQ.3.OR.NFT(JT,5).EQ.3) THEN
4288 X1=HIRND2(6,XMIN1,XMAX1)
4289 X2=HIRND2(6,XMIN2,XMAX2)
4290 ENDIF
4291
4292 IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
4293 & ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
4294 X1=HIRND2(5,XMIN1,XMAX1)
4295 ENDIF
4296 IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
4297 & ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
4298 X2=HIRND2(5,XMIN2,XMAX2)
4299 ENDIF
4300
4301
4302
4303
4304
4305
4306 IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000) X1=1.0-X1
4307 XXP=X1*(1.0-X2)
4308 XXT=X2*(1.0-X1)
4309 IF(XXP.LT.(D1+1.E-4/SW) .OR. XXT.LT.(D2+1.E-4/SW)) THEN
4310 MISS4=MISS4+1
4311 IF(MISS4.LE.1000) GO TO 410
4312 GO TO 5000
4313 ENDIF
4314 IHNT2(13)=3
4315
4316
4317 600 CONTINUE
4318 IF(X1*(1.0-X2).LT.(AMPN**2-1.E-4)/SW.OR.
4319 & X2*(1.0-X1).LT.(AMTN**2-1.E-4)/SW) THEN
4320 MISS=MISS+1
4321 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 2000
4322 PKC=0.0
4323 GO TO 30
4324 ENDIF
4325
4326 EPP=(1.0-X2)*WP
4327 EPM=X1*WM
4328 ETP=X2*WP
4329 ETM=(1.0-X1)*WM
4330 PP(JP,3)=(EPP-EPM)/2.0
4331 PP(JP,4)=(EPP+EPM)/2.0
4332 IF(EPP*EPM-PTP2.LT.0.0) GO TO 6000
4333 PP(JP,5)=SQRT(EPP*EPM-PTP2)
4334 NFP(JP,3)=NFP3
4335 NFP(JP,5)=NFP5
4336
4337 PT(JT,3)=(ETP-ETM)/2.0
4338 PT(JT,4)=(ETP+ETM)/2.0
4339 IF(ETP*ETM-PTT2.LT.0.0) GO TO 6000
4340 PT(JT,5)=SQRT(ETP*ETM-PTT2)
4341 NFT(JT,3)=NFT3
4342 NFT(JT,5)=NFT5
4343
4344
4345 PP(JP,1)=PP11-PKC11
4346 PP(JP,2)=PP12-PKC12
4347
4348 KICKDIP=1
4349 KICKDIT=1
4350 IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
4351 & ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
4352 KICKDIP=0
4353 ENDIF
4354 IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
4355 & ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
4356 KICKDIT=0
4357 ENDIF
4358 IF((KICKDIP.EQ.0.AND.RLU(NSEED).LT.0.5)
4359 & .OR.(KICKDIP.NE.0.AND.RLU(NSEED)
4360 & .LT.0.5/(1.0+(PKC11**2+PKC12**2)/HIPR1(22)**2))) THEN
4361 PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,6)
4362 PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,7)
4363 PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
4364 & +PP(JP,8)+PKC11
4365 PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
4366 & +PP(JP,9)+PKC12
4367 ELSE
4368 PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,8)
4369 PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,9)
4370 PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
4371 & +PP(JP,6)+PKC11
4372 PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
4373 & +PP(JP,7)+PKC12
4374 ENDIF
4375 PP(JP,1)=PP(JP,6)+PP(JP,8)
4376 PP(JP,2)=PP(JP,7)+PP(JP,9)
4377
4378 PT(JT,1)=PT11-PKC21
4379 PT(JT,2)=PT12-PKC22
4380 IF((KICKDIT.EQ.0.AND.RLU(NSEED).LT.0.5)
4381 & .OR.(KICKDIT.NE.0.AND.RLU(NSEED)
4382 & .LT.0.5/(1.0+(PKC21**2+PKC22**2)/HIPR1(22)**2))) THEN
4383 PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,6)
4384 PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,7)
4385 PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
4386 & +PT(JT,8)+PKC21
4387 PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
4388 & +PT(JT,9)+PKC22
4389 ELSE
4390 PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,8)
4391 PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,9)
4392 PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
4393 & +PT(JT,6)+PKC21
4394 PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
4395 & +PT(JT,7)+PKC22
4396 ENDIF
4397 PT(JT,1)=PT(JT,6)+PT(JT,8)
4398 PT(JT,2)=PT(JT,7)+PT(JT,9)
4399
4400
4401 IF(NPJ(JP).NE.0) NFP(JP,5)=3
4402 IF(NTJ(JT).NE.0) NFT(JT,5)=3
4403
4404 IF(EPP/(EPM+0.0001).LT.ETP/(ETM+0.0001).AND.
4405 & ABS(NFP(JP,1)*NFP(JP,2)).LT.1000000)THEN
4406 DO 620 JSB=1,15
4407 PSB=PP(JP,JSB)
4408 PP(JP,JSB)=PT(JT,JSB)
4409 PT(JT,JSB)=PSB
4410 NSB=NFP(JP,JSB)
4411 NFP(JP,JSB)=NFT(JT,JSB)
4412 NFT(JT,JSB)=NSB
4413 620 CONTINUE
4414
4415
4416 ENDIF
4417
4418 RETURN
4419
4420
4421 1000 IERROR=1
4422 IF(IHPR2(10).EQ.0) RETURN
4423 WRITE(6,*) ' Fatal HIJSFT start error,abandon this event'
4424 WRITE(6,*) ' PROJ E+,E-,W+',EPP,EPM,WP
4425 WRITE(6,*) ' TARG E+,E-,W-',ETP,ETM,WM
4426 WRITE(6,*) ' W+*W-, (APN+ATN)^2',SW,SNN
4427 RETURN
4428 2000 IERROR=0
4429 IF(IHPR2(10).EQ.0) RETURN
4430 WRITE(6,*) ' (2)energy partition fail,'
4431 WRITE(6,*) ' HIJSFT not performed, but continue'
4432 WRITE(6,*) ' MP1,MPN',X1*(1.0-X2)*SW,AMPN**2
4433 WRITE(6,*) ' MT2,MTN',X2*(1.0-X1)*SW,AMTN**2
4434 RETURN
4435 3000 IERROR=0
4436 IF(IHPR2(10).EQ.0) RETURN
4437 WRITE(6,*) ' (3)something is wrong with the pt kick, '
4438 WRITE(6,*) ' HIJSFT not performed, but continue'
4439 WRITE(6,*) ' D1=',D1,' D2=',D2,' SW=',SW
4440 WRITE(6,*) ' HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
4441 WRITE(6,*) ' THIS COLLISON NFP5=',NFP5, ' NFT5=',NFT5
4442 WRITE(6,*) ' # OF JET IN PROJ',NPJ(JP),' IN TARG',NTJ(JT)
4443 RETURN
4444 4000 IERROR=0
4445 IF(IHPR2(10).EQ.0) RETURN
4446 WRITE(6,*) ' (4)unable to choose process, but not harmful'
4447 WRITE(6,*) ' HIJSFT not performed, but continue'
4448 WRITE(6,*) ' PTP=',SQRT(PTP2),' PTT=',SQRT(PTT2),' SW=',SW
4449 WRITE(6,*) ' AMCUT=',AMX,' JP=',JP,' JT=',JT
4450 WRITE(6,*) ' HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
4451 RETURN
4452 5000 IERROR=0
4453 IF(IHPR2(10).EQ.0) RETURN
4454 WRITE(6,*) ' energy partition failed(5),for limited try'
4455 WRITE(6,*) ' HIJSFT not performed, but continue'
4456 WRITE(6,*) ' NFP5=',NFP5,' NFT5=',NFT5
4457 WRITE(6,*) ' D1',D1,' X1(1-X2)',X1*(1.0-X2)
4458 WRITE(6,*) ' D2',D2,' X2(1-X1)',X2*(1.0-X1)
4459 RETURN
4460 6000 PKC=0.0
4461 MISS=MISS+1
4462 IF(MISS.LT.100) GO TO 30
4463 IERROR=1
4464 IF(IHPR2(10).EQ.0) RETURN
4465 WRITE(6,*) ' ERROR OCCURED, HIJSFT NOT PERFORMED'
4466 WRITE(6,*) ' Abort this event'
4467 WRITE(6,*) 'MTP,PTP2',EPP*EPM,PTP2,' MTT,PTT2',ETP*ETM,PTT2
4468 RETURN
4469 END
4470
4471
4472
4473 SUBROUTINE HIJFLV(ID)
4474 COMMON/RANSEED/NSEED
4475 SAVE
4476 ID=1
4477 RNID=RLU(NSEED)
4478 IF(RNID.GT.0.43478) THEN
4479 ID=2
4480 IF(RNID.GT.0.86956) ID=3
4481 ENDIF
4482 RETURN
4483 END
4484
4485
4486
4487 SUBROUTINE HIPTDI(PT,PTMAX,IOPT)
4488 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4489 COMMON/RANSEED/NSEED
4490 SAVE
4491 IF(IOPT.EQ.2) THEN
4492 PT=HIRND2(7,0.0,PTMAX)
4493 IF(PT.GT.HIPR1(8))
4494 & PT=HIPR1(2)*SQRT(-ALOG(EXP(-HIPR1(8)**2/HIPR1(2)**2)
4495 & -RLU(NSEED)*(EXP(-HIPR1(8)**2/HIPR1(2)**2)-
4496 & EXP(-PTMAX**2/HIPR1(2)**2))))
4497
4498 ELSE
4499 PT=HIPR1(2)*SQRT(-ALOG(1.0-RLU(NSEED)*
4500 & (1.0-EXP(-PTMAX**2/HIPR1(2)**2))))
4501 ENDIF
4502 PTMAX0=MAX(PTMAX,0.01)
4503 PT=MIN(PTMAX0-0.01,PT)
4504 RETURN
4505 END
4506
4507
4508
4509
4510
4511
4512
4513 SUBROUTINE HIJWDS(IA,IDH,XHIGH)
4514
4515
4516 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4517 COMMON/WOOD/R,D,FNORM,W
4518 DIMENSION IAA(20),RR(20),DD(20),WW(20),RMS(20)
4519 EXTERNAL RWDSAX,WDSAX
4520 SAVE
4521
4522
4523
4524 DATA IAA/2,4,12,16,27,32,40,56,63,93,184,197,208,7*0./
4525 DATA RR/0.01,.964,2.355,2.608,2.84,3.458,3.766,3.971,4.214,
4526 1 4.87,6.51,6.38,6.624,7*0./
4527 DATA DD/0.5882,.322,.522,.513,.569,.61,.586,.5935,.586,.573,
4528 1 .535,.535,.549,7*0./
4529 DATA WW/0.0,.517,-0.149,-0.051,0.,-0.208,-0.161,13*0./
4530 DATA RMS/2.11,1.71,2.46,2.73,3.05,3.247,3.482,3.737,3.925,4.31,
4531 1 5.42,5.33,5.521,7*0./
4532
4533 A=IA
4534
4535
4536 D=0.54
4537
4538 R=1.19*A**(1./3.) - 1.61*A**(-1./3.)
4539
4540 W=0.
4541
4542
4543
4544 DO 10 I=1,13
4545 IF (IA.EQ.IAA(I)) THEN
4546 R=RR(I)
4547 D=DD(I)
4548 W=WW(I)
4549 RS=RMS(I)
4550 END IF
4551 10 CONTINUE
4552
4553 FNORM=1.0
4554 XLOW=0.
4555 XHIGH=R+ 12.*D
4556 IF (W.LT.-0.01) THEN
4557 IF (XHIGH.GT.R/SQRT(ABS(W))) XHIGH=R/SQRT(ABS(W))
4558 END IF
4559 FGAUS=GAUSS1(RWDSAX,XLOW,XHIGH,0.001)
4560 FNORM=1./FGAUS
4561
4562 IF (IDH.EQ.1) THEN
4563 HINT1(72)=R
4564 HINT1(73)=D
4565 HINT1(74)=W
4566 HINT1(75)=FNORM/4.0/HIPR1(40)
4567 ELSE IF (IDH.EQ.2) THEN
4568 HINT1(76)=R
4569 HINT1(77)=D
4570 HINT1(78)=W
4571 HINT1(79)=FNORM/4.0/HIPR1(40)
4572 ENDIF
4573
4574
4575
4576 CALL HIFUN(IDH,XLOW,XHIGH,RWDSAX)
4577 RETURN
4578 END
4579
4580
4581 FUNCTION WDSAX(X)
4582
4583 COMMON/WOOD/R,D,FNORM,W
4584 SAVE
4585 WDSAX=FNORM*(1.+W*(X/R)**2)/(1+EXP((X-R)/D))
4586 IF (W.LT.0.) THEN
4587 IF (X.GE.R/SQRT(ABS(W))) WDSAX=0.
4588 ENDIF
4589 RETURN
4590 END
4591
4592
4593 FUNCTION RWDSAX(X)
4594 RWDSAX=X*X*WDSAX(X)
4595 RETURN
4596 END
4597
4598
4599
4600
4601 FUNCTION WDSAX1(X)
4602
4603
4604
4605
4606 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4607 SAVE
4608 WDSAX1=HINT1(75)*(1.+HINT1(74)*(X/HINT1(72))**2)/
4609 & (1+EXP((X-HINT1(72))/HINT1(73)))
4610 IF (HINT1(74).LT.0.) THEN
4611 IF (X.GE.HINT1(72)/SQRT(ABS(HINT1(74)))) WDSAX1=0.
4612 ENDIF
4613 RETURN
4614 END
4615
4616
4617 FUNCTION WDSAX2(X)
4618
4619
4620
4621
4622 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4623 SAVE
4624 WDSAX2=HINT1(79)*(1.+HINT1(78)*(X/HINT1(76))**2)/
4625 & (1+EXP((X-HINT1(76))/HINT1(77)))
4626 IF (HINT1(78).LT.0.) THEN
4627 IF (X.GE.HINT1(76)/SQRT(ABS(HINT1(78)))) WDSAX2=0.
4628 ENDIF
4629 RETURN
4630 END
4631
4632
4633
4634
4635
4636 FUNCTION PROFILE(XB)
4637 COMMON/PACT/BB,B1,PHI,Z1
4638 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4639 EXTERNAL FLAP, FLAP1, FLAP2
4640 SAVE
4641
4642 BB=XB
4643 PROFILE=1.0
4644 IF(IHNT2(1).GT.1 .AND. IHNT2(3).GT.1) THEN
4645 PROFILE=float(IHNT2(1))*float(IHNT2(3))*0.1*
4646 & GAUSS1(FLAP,0.0,HIPR1(34),0.01)
4647 ELSE IF(IHNT2(1).EQ.1 .AND. IHNT2(3).GT.1) THEN
4648 PROFILE=0.2*float(IHNT2(3))*
4649 & GAUSS1(FLAP2,0.0,HIPR1(35),0.001)
4650 ELSE IF(IHNT2(1).GT.1 .AND. IHNT2(3).EQ.1) THEN
4651 PROFILE=0.2*float(IHNT2(1))*
4652 & GAUSS1(FLAP1,0.0,HIPR1(34),0.001)
4653 ENDIF
4654 RETURN
4655 END
4656
4657
4658 FUNCTION FLAP(X)
4659 COMMON/PACT/BB,B1,PHI,Z1
4660 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4661 EXTERNAL FGP1
4662 SAVE
4663 B1=X
4664 FLAP=GAUSS2(FGP1,0.0,2.0*HIPR1(40),0.01)
4665 RETURN
4666 END
4667
4668 FUNCTION FGP1(X)
4669 COMMON/PACT/BB,B1,PHI,Z1
4670 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4671 EXTERNAL FGP2
4672 SAVE
4673 PHI=X
4674 FGP1=2.0*GAUSS3(FGP2,0.0,HIPR1(34),0.01)
4675 RETURN
4676 END
4677
4678 FUNCTION FGP2(X)
4679 COMMON/PACT/BB,B1,PHI,Z1
4680 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4681 EXTERNAL FGP3
4682 SAVE
4683 Z1=X
4684 FGP2=2.0*GAUSS4(FGP3,0.0,HIPR1(35),0.01)
4685 RETURN
4686 END
4687
4688 FUNCTION FGP3(X)
4689 COMMON/PACT/BB,B1,PHI,Z1
4690 SAVE
4691 R1=SQRT(B1**2+Z1**2)
4692 R2=SQRT(BB**2+B1**2-2.0*B1*BB*COS(PHI)+X**2)
4693 FGP3=B1*WDSAX1(R1)*WDSAX2(R2)
4694 RETURN
4695 END
4696
4697
4698 FUNCTION FLAP1(X)
4699 COMMON/PACT/BB,B1,PHI,Z1
4700 SAVE
4701 R=SQRT(BB**2+X**2)
4702 FLAP1=WDSAX1(R)
4703 RETURN
4704 END
4705
4706
4707 FUNCTION FLAP2(X)
4708 COMMON/PACT/BB,B1,PHI,Z1
4709 SAVE
4710 R=SQRT(BB**2+X**2)
4711 FLAP2=WDSAX2(R)
4712 RETURN
4713 END
4714
4715
4716
4717
4718
4719
4720
4721
4722 SUBROUTINE HIFUN(I,XMIN,XMAX,FHB)
4723 COMMON/HIJHB/RR(10,4097),XX(10,4097)
4724 EXTERNAL FHB
4725 SAVE
4726 FNORM=GAUSS1(FHB,XMIN,XMAX,0.001)
4727 DO 100 J=1,4097
4728 XX(I,J)=XMIN+(XMAX-XMIN)*(J-1)/4096.0
4729 XDD=XX(I,J)
4730 RR(I,J)=GAUSS1(FHB,XMIN,XDD,0.001)/FNORM
4731 100 CONTINUE
4732 RETURN
4733 END
4734
4735
4736
4737 FUNCTION HIRND(I)
4738 COMMON/HIJHB/RR(10,4097),XX(10,4097)
4739 COMMON/RANSEED/NSEED
4740 SAVE
4741 RX=RLU(NSEED)
4742 JL=0
4743 JU=4098
4744 10 IF(JU-JL.GT.1) THEN
4745 JM=(JU+JL)/2
4746 IF((RR(I,4097).GT.RR(I,1)).EQV.(RX.GT.RR(I,JM))) THEN
4747 JL=JM
4748 ELSE
4749 JU=JM
4750 ENDIF
4751 GO TO 10
4752 ENDIF
4753 J=JL
4754 IF(J.LT.1) J=1
4755 IF(J.GE.4097) J=4096
4756
4757
4758 FRAC=(RX-RR(I,J))/(RR(I,J+1)-RR(I,J))
4759 HIRND=XX(I,J)+FRAC*(XX(I,J+1)-XX(I,J))
4760 RETURN
4761 END
4762
4763
4764
4765
4766
4767 FUNCTION HIRND2(I,XMIN,XMAX)
4768 COMMON/HIJHB/RR(10,4097),XX(10,4097)
4769 COMMON/RANSEED/NSEED
4770 SAVE
4771 IF(XMIN.LT.XX(I,1)) XMIN=XX(I,1)
4772 IF(XMAX.GT.XX(I,4097)) XMAX=XX(I,4097)
4773 JMIN=1+4096*(XMIN-XX(I,1))/(XX(I,4097)-XX(I,1))
4774 JMAX=1+4096*(XMAX-XX(I,1))/(XX(I,4097)-XX(I,1))
4775 RX=RR(I,JMIN)+(RR(I,JMAX)-RR(I,JMIN))*RLU(NSEED)
4776 JL=0
4777 JU=4098
4778 10 IF(JU-JL.GT.1) THEN
4779 JM=(JU+JL)/2
4780 IF((RR(I,4097).GT.RR(I,1)).EQV.(RX.GT.RR(I,JM))) THEN
4781 JL=JM
4782 ELSE
4783 JU=JM
4784 ENDIF
4785 GO TO 10
4786 ENDIF
4787 J=JL
4788 IF(J.LT.1) J=1
4789 IF(J.GE.4097) J=4096
4790
4791
4792 FRAC=(RX-RR(I,J))/(RR(I,J+1)-RR(I,J))
4793 HIRND2=XX(I,J)+FRAC*(XX(I,J+1)-XX(I,J))
4794 RETURN
4795 END
4796
4797
4798
4799
4800 SUBROUTINE HIJCRS
4801
4802
4803 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4804 COMMON/NJET/N,IP_CRS
4805 EXTERNAL FHIN,FTOT,FNJET,FTOTJET,FTOTRIG
4806 SAVE
4807 IF(HINT1(1).GE.10.0) CALL CRSJET
4808
4809
4810 APHX1=HIPR1(6)*(IHNT2(1)**0.3333333-1.0)
4811 APHX2=HIPR1(6)*(IHNT2(3)**0.3333333-1.0)
4812 HINT1(11)=HINT1(14)-APHX1*HINT1(15)
4813 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
4814 HINT1(10)=GAUSS1(FTOTJET,0.0,20.0,0.01)
4815 HINT1(12)=GAUSS1(FHIN,0.0,20.0,0.01)
4816 HINT1(13)=GAUSS1(FTOT,0.0,20.0,0.01)
4817 HINT1(60)=HINT1(61)-APHX1*HINT1(62)
4818 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
4819 HINT1(59)=GAUSS1(FTOTRIG,0.0,20.0,0.01)
4820 IF(HINT1(59).EQ.0.0) HINT1(59)=HINT1(60)
4821 IF(HINT1(1).GE.10.0) Then
4822 DO 20 I=0,20
4823 N=I
4824 HINT1(80+I)=GAUSS1(FNJET,0.0,20.0,0.01)/HINT1(12)
4825 20 CONTINUE
4826 ENDIF
4827 HINT1(10)=HINT1(10)*HIPR1(31)
4828 HINT1(12)=HINT1(12)*HIPR1(31)
4829 HINT1(13)=HINT1(13)*HIPR1(31)
4830 HINT1(59)=HINT1(59)*HIPR1(31)
4831
4832
4833 IF(IHPR2(13).NE.0) THEN
4834 HIPR1(33)=1.36*(1.0+36.0/HINT1(1)**2)
4835 & *ALOG(0.6+0.1*HINT1(1)**2)
4836 HIPR1(33)=HIPR1(33)/HINT1(12)
4837 ENDIF
4838
4839
4840 RETURN
4841 END
4842
4843
4844
4845
4846 FUNCTION FTOT(X)
4847 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4848 SAVE
4849 OMG=OMG0(X)*(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
4850 FTOT=2.0*(1.0-EXP(-OMG))
4851 RETURN
4852 END
4853
4854
4855
4856 FUNCTION FHIN(X)
4857 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4858 SAVE
4859 OMG=OMG0(X)*(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
4860 FHIN=1.0-EXP(-2.0*OMG)
4861 RETURN
4862 END
4863
4864
4865
4866 FUNCTION FTOTJET(X)
4867 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4868 SAVE
4869 OMG=OMG0(X)*HINT1(11)/HIPR1(31)/2.0
4870 FTOTJET=1.0-EXP(-2.0*OMG)
4871 RETURN
4872 END
4873
4874
4875
4876 FUNCTION FTOTRIG(X)
4877 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4878 SAVE
4879 OMG=OMG0(X)*HINT1(60)/HIPR1(31)/2.0
4880 FTOTRIG=1.0-EXP(-2.0*OMG)
4881 RETURN
4882 END
4883
4884
4885
4886
4887 FUNCTION FNJET(X)
4888 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4889 COMMON/NJET/N,IP_CRS
4890 SAVE
4891 OMG1=OMG0(X)*HINT1(11)/HIPR1(31)
4892 C0=EXP(N*ALOG(OMG1)-SGMIN(N+1))
4893 IF(N.EQ.0) C0=1.0-EXP(-2.0*OMG0(X)*HIPR1(30)/HIPR1(31)/2.0)
4894 FNJET=C0*EXP(-OMG1)
4895 RETURN
4896 END
4897
4898
4899
4900
4901
4902 FUNCTION SGMIN(N)
4903 GA=0.
4904 IF(N.LE.2) GO TO 20
4905 DO 10 I=1,N-1
4906 Z=I
4907 GA=GA+ALOG(Z)
4908 10 CONTINUE
4909 20 SGMIN=GA
4910 RETURN
4911 END
4912
4913
4914
4915 FUNCTION OMG0(X)
4916 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4917 COMMON /BESEL/X4
4918 EXTERNAL BK
4919 SAVE
4920 X4=HIPR1(32)*SQRT(X)
4921 OMG0=HIPR1(32)**2*GAUSS2(BK,X4,X4+20.0,0.01)/96.0
4922 RETURN
4923 END
4924
4925
4926
4927 FUNCTION ROMG(X)
4928
4929
4930 DIMENSION FR(0:1000)
4931 SAVE
4932 DATA I0/0/
4933 IF(I0.NE.0) GO TO 100
4934 DO 50 I=1,1001
4935 XR=(I-1)*0.01
4936 FR(I-1)=OMG0(XR)
4937 50 CONTINUE
4938 100 I0=1
4939 IF(X.GE.10.0) THEN
4940 ROMG=0.0
4941 RETURN
4942 ENDIF
4943 IX=INT(X*100)
4944 ROMG=(FR(IX)*((IX+1)*0.01-X)+FR(IX+1)*(X-IX*0.01))/0.01
4945 RETURN
4946 END
4947
4948
4949
4950 FUNCTION BK(X)
4951 COMMON /BESEL/X4
4952 SAVE
4953 BK=EXP(-X)*(X**2-X4**2)**2.50/15.0
4954 RETURN
4955 END
4956
4957
4958
4959
4960
4961 SUBROUTINE CRSJET
4962 IMPLICIT REAL*8(A-H,O-Z)
4963 REAL HIPR1(100),HINT1(100)
4964 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
4965 COMMON/NJET/N,IP_CRS
4966 COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
4967 COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
4968 COMMON/BVEG3/F,TI,TSI
4969 COMMON/SEEDVAX/NUM1
4970 EXTERNAL FJET,FJETRIG
4971 SAVE
4972
4973
4974
4975
4976
4977
4978 NDIM=3
4979 IP_CRS=0
4980 CALL VEGAS(FJET,AVGI,SD,CHI2A)
4981 HINT1(14)=AVGI/2.5682
4982 IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
4983 IP_CRS=1
4984 CALL VEGAS(FJET,AVGI,SD,CHI2A)
4985 HINT1(15)=AVGI/2.5682
4986 ENDIF
4987 IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
4988 IP_CRS=2
4989 CALL VEGAS(FJET,AVGI,SD,CHI2A)
4990 HINT1(16)=AVGI/2.5682
4991 ENDIF
4992 IF(IHPR2(6).EQ.1.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
4993 IP_CRS=3
4994 CALL VEGAS(FJET,AVGI,SD,CHI2A)
4995 HINT1(17)=AVGI/2.5682
4996 ENDIF
4997
4998
4999 IF(IHPR2(3).NE.0) THEN
5000 IP_CRS=0
5001 CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
5002 HINT1(61)=AVGI/2.5682
5003 IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
5004 IP_CRS=1
5005 CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
5006 HINT1(62)=AVGI/2.5682
5007 ENDIF
5008 IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
5009 IP_CRS=2
5010 CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
5011 HINT1(63)=AVGI/2.5682
5012 ENDIF
5013 IF(IHPR2(6).EQ.1.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
5014 IP_CRS=3
5015 CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
5016 HINT1(64)=AVGI/2.5682
5017 ENDIF
5018 ENDIF
5019
5020
5021 RETURN
5022 END
5023
5024
5025
5026 FUNCTION FJET(X,WGT)
5027 IMPLICIT REAL*8(A-H,O-Z)
5028 REAL HIPR1(100),HINT1(100)
5029 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5030 DIMENSION X(10)
5031 SAVE
5032 WGT=WGT
5033 PT2=(HINT1(1)**2/4.0-HIPR1(8)**2)*X(1)+HIPR1(8)**2
5034 XT=2.0*DSQRT(PT2)/HINT1(1)
5035 YMX1=DLOG(1.0/XT+DSQRT(1.0/XT**2-1.0))
5036 Y1=2.0*YMX1*X(2)-YMX1
5037 YMX2=DLOG(2.0/XT-DEXP(Y1))
5038 YMN2=DLOG(2.0/XT-DEXP(-Y1))
5039 Y2=(YMX2+YMN2)*X(3)-YMN2
5040 FJET=2.0*YMX1*(YMX2+YMN2)*(HINT1(1)**2/4.0-HIPR1(8)**2)
5041 & *G(Y1,Y2,PT2)/2.0
5042 RETURN
5043 END
5044
5045
5046
5047 FUNCTION FJETRIG(X,WGT)
5048 IMPLICIT REAL*8(A-H,O-Z)
5049 REAL HIPR1(100),HINT1(100),PTMAX,PTMIN
5050 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5051 DIMENSION X(10)
5052 SAVE
5053 PTMIN=ABS(HIPR1(10))-0.25
5054 PTMIN=MAX(PTMIN,HIPR1(8))
5055 AM2=0.D0
5056 IF(IHPR2(3).EQ.3) THEN
5057 AM2=HIPR1(7)**2
5058 PTMIN=MAX(0.0,HIPR1(10))
5059 ENDIF
5060 PTMAX=ABS(HIPR1(10))+0.25
5061 IF(HIPR1(10).LE.0.0) PTMAX=HINT1(1)/2.0-AM2
5062 IF(PTMAX.LE.PTMIN) PTMAX=PTMIN+0.25
5063 PT2=(PTMAX**2-PTMIN**2)*X(1)+PTMIN**2
5064 AMT2=PT2+AM2
5065 XT=2.0*DSQRT(AMT2)/HINT1(1)
5066 YMX1=DLOG(1.0/XT+DSQRT(1.0/XT**2-1.0))
5067 Y1=2.0*YMX1*X(2)-YMX1
5068 YMX2=DLOG(2.0/XT-DEXP(Y1))
5069 YMN2=DLOG(2.0/XT-DEXP(-Y1))
5070 Y2=(YMX2+YMN2)*X(3)-YMN2
5071 IF(IHPR2(3).EQ.3) THEN
5072 GTRIG=2.0*GHVQ(Y1,Y2,AMT2)
5073 ELSE IF(IHPR2(3).EQ.2) THEN
5074 GTRIG=2.0*GPHOTON(Y1,Y2,PT2)
5075 ELSE
5076 GTRIG=G(Y1,Y2,PT2)
5077 ENDIF
5078 FJETRIG=2.0*YMX1*(YMX2+YMN2)*(PTMAX**2-PTMIN**2)
5079 & *GTRIG/2.0
5080 RETURN
5081 END
5082
5083
5084
5085 FUNCTION GHVQ(Y1,Y2,AMT2)
5086 IMPLICIT REAL*8 (A-H,O-Z)
5087 REAL HIPR1(100),HINT1(100)
5088 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5089 DIMENSION F(2,7)
5090 SAVE
5091 XT=2.0*DSQRT(AMT2)/HINT1(1)
5092 X1=0.50*XT*(DEXP(Y1)+DEXP(Y2))
5093 X2=0.50*XT*(DEXP(-Y1)+DEXP(-Y2))
5094 SS=X1*X2*HINT1(1)**2
5095 AF=4.0
5096 IF(IHPR2(18).NE.0) AF=5.0
5097 DLAM=HIPR1(15)
5098 APH=12.0*3.1415926/(33.0-2.0*AF)/DLOG(AMT2/DLAM**2)
5099
5100 CALL PARTON(F,X1,X2,AMT2)
5101
5102 Gqq=4.0*(COSH(Y1-Y2)+HIPR1(7)**2/AMT2)/(1.D0+COSH(Y1-Y2))/9.0
5103 & *(F(1,1)*F(2,2)+F(1,2)*F(2,1)+F(1,3)*F(2,4)
5104 & +F(1,4)*F(2,3)+F(1,5)*F(2,6)+F(1,6)*F(2,5))
5105 Ggg=(8.D0*COSH(Y1-Y2)-1.D0)*(COSH(Y1-Y2)+2.0*HIPR1(7)**2/AMT2
5106 & -2.0*HIPR1(7)**4/AMT2**2)/(1.0+COSH(Y1-Y2))/24.D0
5107 & *F(1,7)*F(2,7)
5108
5109 GHVQ=(Gqq+Ggg)*HIPR1(23)*3.14159*APH**2/SS**2
5110 RETURN
5111 END
5112
5113
5114
5115 FUNCTION GPHOTON(Y1,Y2,PT2)
5116 IMPLICIT REAL*8 (A-H,O-Z)
5117 REAL HIPR1(100),HINT1(100)
5118 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5119 DIMENSION F(2,7)
5120 SAVE
5121 XT=2.0*DSQRT(PT2)/HINT1(1)
5122 X1=0.50*XT*(DEXP(Y1)+DEXP(Y2))
5123 X2=0.50*XT*(DEXP(-Y1)+DEXP(-Y2))
5124 Z=DSQRT(1.D0-XT**2/X1/X2)
5125 SS=X1*X2*HINT1(1)**2
5126 T=-(1.0-Z)/2.0
5127 U=-(1.0+Z)/2.0
5128 AF=3.0
5129 DLAM=HIPR1(15)
5130 APH=12.0*3.1415926/(33.0-2.0*AF)/DLOG(PT2/DLAM**2)
5131 APHEM=1.0/137.0
5132
5133 CALL PARTON(F,X1,X2,PT2)
5134
5135 G11=-(U**2+1.0)/U/3.0*F(1,7)*(4.0*F(2,1)+4.0*F(2,2)
5136 & +F(2,3)+F(2,4)+F(2,5)+F(2,6))/9.0
5137 G12=-(T**2+1.0)/T/3.0*F(2,7)*(4.0*F(1,1)+4.0*F(1,2)
5138 & +F(1,3)+F(1,4)+F(1,5)+F(1,6))/9.0
5139 G2=8.0*(U**2+T**2)/U/T/9.0*(4.0*F(1,1)*F(2,2)
5140 & +4.0*F(1,2)*F(2,1)+F(1,3)*F(2,4)+F(1,4)*F(2,3)
5141 & +F(1,5)*F(2,6)+F(1,6)*F(2,5))/9.0
5142
5143 GPHOTON=(G11+G12+G2)*HIPR1(23)*3.14159*APH*APHEM/SS**2
5144 RETURN
5145 END
5146
5147
5148
5149
5150 FUNCTION G(Y1,Y2,PT2)
5151 IMPLICIT REAL*8 (A-H,O-Z)
5152 REAL HIPR1(100),HINT1(100)
5153 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5154 DIMENSION F(2,7)
5155 SAVE
5156 XT=2.0*DSQRT(PT2)/HINT1(1)
5157 X1=0.50*XT*(DEXP(Y1)+DEXP(Y2))
5158 X2=0.50*XT*(DEXP(-Y1)+DEXP(-Y2))
5159 Z=DSQRT(1.D0-XT**2/X1/X2)
5160 SS=X1*X2*HINT1(1)**2
5161 T=-(1.0-Z)/2.0
5162 U=-(1.0+Z)/2.0
5163 AF=3.0
5164 DLAM=HIPR1(15)
5165 APH=12.0*3.1415926/(33.0-2.0*AF)/DLOG(PT2/DLAM**2)
5166
5167 CALL PARTON(F,X1,X2,PT2)
5168
5169 G11=( (F(1,1)+F(1,2))*(F(2,3)+F(2,4)+F(2,5)+F(2,6))
5170 & +(F(1,3)+F(1,4))*(F(2,5)+F(2,6)) )*SUBCRS1(T,U)
5171
5172 G12=( (F(2,1)+F(2,2))*(F(1,3)+F(1,4)+F(1,5)+F(1,6))
5173 & +(F(2,3)+F(2,4))*(F(1,5)+F(1,6)) )*SUBCRS1(U,T)
5174
5175 G13=(F(1,1)*F(2,1)+F(1,2)*F(2,2)+F(1,3)*F(2,3)+F(1,4)*F(2,4)
5176 & +F(1,5)*F(2,5)+F(1,6)*F(2,6))*(SUBCRS1(U,T)
5177 & +SUBCRS1(T,U)-8.D0/T/U/27.D0)
5178
5179 G2=(AF-1)*(F(1,1)*F(2,2)+F(2,1)*F(1,2)+F(1,3)*F(2,4)
5180 & +F(2,3)*F(1,4)+F(1,5)*F(2,6)+F(2,5)*F(1,6))*SUBCRS2(T,U)
5181
5182 G31=(F(1,1)*F(2,2)+F(1,3)*F(2,4)+F(1,5)*F(2,6))*SUBCRS3(T,U)
5183 G32=(F(2,1)*F(1,2)+F(2,3)*F(1,4)+F(2,5)*F(1,6))*SUBCRS3(U,T)
5184
5185 G4=(F(1,1)*F(2,2)+F(2,1)*F(1,2)+F(1,3)*F(2,4)+F(2,3)*F(1,4)+
5186 1 F(1,5)*F(2,6)+F(2,5)*F(1,6))*SUBCRS4(T,U)
5187
5188 G5=AF*F(1,7)*F(2,7)*SUBCRS5(T,U)
5189
5190 G61=F(1,7)*(F(2,1)+F(2,2)+F(2,3)+F(2,4)+F(2,5)
5191 & +F(2,6))*SUBCRS6(T,U)
5192 G62=F(2,7)*(F(1,1)+F(1,2)+F(1,3)+F(1,4)+F(1,5)
5193 & +F(1,6))*SUBCRS6(U,T)
5194
5195 G7=F(1,7)*F(2,7)*SUBCRS7(T,U)
5196
5197 G=(G11+G12+G13+G2+G31+G32+G4+G5+G61+G62+G7)*HIPR1(17)*
5198 1 3.14159D0*APH**2/SS**2
5199 RETURN
5200 END
5201
5202
5203
5204 FUNCTION SUBCRS1(T,U)
5205 IMPLICIT REAL*8(A-H,O-Z)
5206 SUBCRS1=4.D0/9.D0*(1.D0+U**2)/T**2
5207 RETURN
5208 END
5209
5210
5211 FUNCTION SUBCRS2(T,U)
5212 IMPLICIT REAL*8(A-H,O-Z)
5213 SUBCRS2=4.D0/9.D0*(T**2+U**2)
5214 RETURN
5215 END
5216
5217
5218 FUNCTION SUBCRS3(T,U)
5219 IMPLICIT REAL*8(A-H,O-Z)
5220 SUBCRS3=4.D0/9.D0*(T**2+U**2+(1.D0+U**2)/T**2
5221 1 -2.D0*U**2/3.D0/T)
5222 RETURN
5223 END
5224
5225
5226 FUNCTION SUBCRS4(T,U)
5227 IMPLICIT REAL*8(A-H,O-Z)
5228 SUBCRS4=8.D0/3.D0*(T**2+U**2)*(4.D0/9.D0/T/U-1.D0)
5229 RETURN
5230 END
5231
5232
5233
5234 FUNCTION SUBCRS5(T,U)
5235 IMPLICIT REAL*8(A-H,O-Z)
5236 SUBCRS5=3.D0/8.D0*(T**2+U**2)*(4.D0/9.D0/T/U-1.D0)
5237 RETURN
5238 END
5239
5240
5241 FUNCTION SUBCRS6(T,U)
5242 IMPLICIT REAL*8(A-H,O-Z)
5243 SUBCRS6=(1.D0+U**2)*(1.D0/T**2-4.D0/U/9.D0)
5244 RETURN
5245 END
5246
5247
5248 FUNCTION SUBCRS7(T,U)
5249 IMPLICIT REAL*8(A-H,O-Z)
5250 SUBCRS7=9.D0/2.D0*(3.D0-T*U-U/T**2-T/U**2)
5251 RETURN
5252 END
5253
5254
5255
5256 SUBROUTINE PARTON(F,X1,X2,QQ)
5257 IMPLICIT REAL*8(A-H,O-Z)
5258 REAL HIPR1(100),HINT1(100)
5259 COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5260 COMMON/NJET/N,IP_CRS
5261 DIMENSION F(2,7)
5262 SAVE
5263 DLAM=HIPR1(15)
5264 Q0=HIPR1(16)
5265 S=DLOG(DLOG(QQ/DLAM**2)/DLOG(Q0**2/DLAM**2))
5266 IF(IHPR2(7).EQ.2) GO TO 200
5267
5268 AT1=0.419+0.004*S-0.007*S**2
5269 AT2=3.460+0.724*S-0.066*S**2
5270 GMUD=4.40-4.86*S+1.33*S**2
5271 AT3=0.763-0.237*S+0.026*S**2
5272 AT4=4.00+0.627*S-0.019*S**2
5273 GMD=-0.421*S+0.033*S**2
5274
5275 CAS=1.265-1.132*S+0.293*S**2
5276 AS=-0.372*S-0.029*S**2
5277 BS=8.05+1.59*S-0.153*S**2
5278 APHS=6.31*S-0.273*S**2
5279 BTAS=-10.5*S-3.17*S**2
5280 GMS=14.7*S+9.80*S**2
5281
5282
5283
5284
5285
5286
5287
5288
5289 CAG=1.56-1.71*S+0.638*S**2
5290 AG=-0.949*S+0.325*S**2
5291 BG=6.0+1.44*S-1.05*S**2
5292 APHG=9.0-7.19*S+0.255*S**2
5293 BTAG=-16.5*S+10.9*S**2
5294 GMG=15.3*S-10.1*S**2
5295 GO TO 300
5296
5297 200 AT1=0.374+0.014*S
5298 AT2=3.33+0.753*S-0.076*S**2
5299 GMUD=6.03-6.22*S+1.56*S**2
5300 AT3=0.761-0.232*S+0.023*S**2
5301 AT4=3.83+0.627*S-0.019*S**2
5302 GMD=-0.418*S+0.036*S**2
5303
5304 CAS=1.67-1.92*S+0.582*S**2
5305 AS=-0.273*S-0.164*S**2
5306 BS=9.15+0.530*S-0.763*S**2
5307 APHS=15.7*S-2.83*S**2
5308 BTAS=-101.0*S+44.7*S**2
5309 GMS=223.0*S-117.0*S**2
5310
5311
5312
5313
5314
5315
5316
5317
5318 CAG=0.879-0.971*S+0.434*S**2
5319 AG=-1.16*S+0.476*S**2
5320 BG=4.0+1.23*S-0.254*S**2
5321 APHG=9.0-5.64*S-0.817*S**2
5322 BTAG=-7.54*S+5.50*S**2
5323 GMG=-0.596*S+1.26*S**2
5324
5325 300 B12=DEXP(GMRE(AT1)+GMRE(AT2+1.D0)-GMRE(AT1+AT2+1.D0))
5326 B34=DEXP(GMRE(AT3)+GMRE(AT4+1.D0)-GMRE(AT3+AT4+1.D0))
5327 CNUD=3.D0/B12/(1.D0+GMUD*AT1/(AT1+AT2+1.D0))
5328 CND=1.D0/B34/(1.D0+GMD*AT3/(AT3+AT4+1.D0))
5329
5330
5331
5332
5333 FUD1=CNUD*X1**AT1*(1.D0-X1)**AT2*(1.D0+GMUD*X1)
5334 FS1=CAS*X1**AS*(1.D0-X1)**BS*(1.D0+APHS*X1
5335 & +BTAS*X1**2+GMS*X1**3)
5336 F(1,3)=CND*X1**AT3*(1.D0-X1)**AT4*(1.D0+GMD*X1)+FS1/6.D0
5337 F(1,1)=FUD1-F(1,3)+FS1/3.D0
5338 F(1,2)=FS1/6.D0
5339 F(1,4)=FS1/6.D0
5340 F(1,5)=FS1/6.D0
5341 F(1,6)=FS1/6.D0
5342 F(1,7)=CAG*X1**AG*(1.D0-X1)**BG*(1.D0+APHG*X1
5343 & +BTAG*X1**2+GMG*X1**3)
5344
5345 FUD2=CNUD*X2**AT1*(1.D0-X2)**AT2*(1.D0+GMUD*X2)
5346 FS2=CAS*X2**AS*(1.D0-X2)**BS*(1.D0+APHS*X2
5347 & +BTAS*X2**2+GMS*X2**3)
5348 F(2,3)=CND*X2**AT3*(1.D0-X2)**AT4*(1.D0+GMD*X2)+FS2/6.D0
5349 F(2,1)=FUD2-F(2,3)+FS2/3.D0
5350 F(2,2)=FS2/6.D0
5351 F(2,4)=FS2/6.D0
5352 F(2,5)=FS2/6.D0
5353 F(2,6)=FS2/6.D0
5354 F(2,7)=CAG*X2**AG*(1.D0-X2)**BG*(1.D0+APHG*X2
5355 & +BTAG*X2**2+GMG*X2**3)
5356
5357
5358 IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
5359 AAX=1.193*ALOG(FLOAT(IHNT2(1)))**0.16666666
5360 RRX=AAX*(X1**3-1.2*X1**2+0.21*X1)+1.0
5361 & +1.079*(FLOAT(IHNT2(1))**0.33333333-1.0)
5362 & /DLOG(IHNT2(1)+1.0D0)*DSQRT(X1)*DEXP(-X1**2/0.01)
5363 IF(IP_CRS.EQ.1 .OR.IP_CRS.EQ.3) RRX=DEXP(-X1**2/0.01)
5364 DO 400 I=1,7
5365 F(1,I)=RRX*F(1,I)
5366 400 CONTINUE
5367 ENDIF
5368 IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
5369 AAX=1.193*ALOG(FLOAT(IHNT2(3)))**0.16666666
5370 RRX=AAX*(X2**3-1.2*X2**2+0.21*X2)+1.0
5371 & +1.079*(FLOAT(IHNT2(3))**0.33333-1.0)
5372 & /DLOG(IHNT2(3)+1.0D0)*DSQRT(X2)*DEXP(-X2**2/0.01)
5373 IF(IP_CRS.EQ.2 .OR. IP_CRS.EQ.3) RRX=DEXP(-X2**2/0.01)
5374 DO 500 I=1,7
5375 F(2,I)=RRX*F(2,I)
5376 500 CONTINUE
5377 ENDIF
5378
5379 RETURN
5380 END
5381
5382
5383
5384 FUNCTION GMRE(X)
5385 IMPLICIT REAL*8(A-H,O-Z)
5386 Z=X
5387 IF(X.GT.3.0D0) GO TO 10
5388 Z=X+3.D0
5389 10 GMRE=0.5D0*DLOG(2.D0*3.14159265D0/Z)+Z*DLOG(Z)-Z+DLOG(1.D0
5390 1 +1.D0/12.D0/Z+1.D0/288.D0/Z**2-139.D0/51840.D0/Z**3
5391 1 -571.D0/2488320.D0/Z**4)
5392 IF(Z.EQ.X) GO TO 20
5393 GMRE=GMRE-DLOG(Z-1.D0)-DLOG(Z-2.D0)-DLOG(Z-3.D0)
5394 20 CONTINUE
5395 RETURN
5396 END
5397
5398
5399
5400 FUNCTION GMIN(N)
5401 IMPLICIT REAL*8(A-H,O-Z)
5402 GA=0.
5403 IF(N.LE.2) GO TO 20
5404 DO 10 I=1,N-1
5405 Z=I
5406 GA=GA+DLOG(Z)
5407 10 CONTINUE
5408 20 GMIN=GA
5409 RETURN
5410 END
5411
5412
5413
5414
5415 BLOCK DATA HIDATA
5416 REAL*8 XL(10),XU(10),ACC
5417 COMMON/BVEG1/XL,XU,ACC,NDIM,NCALL,ITMX,NPRN
5418 COMMON/SEEDVAX/NUM1
5419 COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5420 COMMON/RANSEED/NSEED
5421 COMMON/HIMAIN1/ NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
5422 COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4),VATT(130000,4)
5423 COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
5424 COMMON/HIJCRDN/YP(3,300),YT(3,300)
5425 COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
5426 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
5427 & PJPM(300,500),NTJ(300),KFTJ(300,500),
5428 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
5429 & PJTE(300,500),PJTM(300,500)
5430 COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
5431 & K2SG(900,100),PXSG(900,100),PYSG(900,100),
5432 & PZSG(900,100),PESG(900,100),PMSG(900,100)
5433 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
5434 COMMON/HIPYINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
5435 SAVE
5436 DATA NUM1/30123984/,XL/10*0.D0/,XU/10*1.D0/
5437 DATA NCALL/1000/ITMX/100/ACC/0.01/NPRN/0/
5438
5439 DATA NSEED/74769375/
5440 DATA HIPR1/
5441 & 1.5, 0.35, 0.5, 0.9, 2.0, 0.1, 1.5, 2.0, -1.0, -2.25,
5442 & 2.0, 0.5, 1.0, 2.0, 0.2, 2.0, 2.5, 0.3, 0.1, 1.4,
5443 & 1.6, 1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 57.0,
5444 & 28.5, 3.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
5445 & 3.141592654,
5446 & 0.0, 0.4, 0.1, 1.5, 0.1, 0.25, 0.0, 0.5, 0.0, 0.0,
5447 & 50*0.0/
5448
5449 DATA IHPR2/
5450 & 1, 3, 0, 1, 1, 1, 1, 10, 0, 0,
5451 & 1, 0, 1, 1, 0, 0, 1, 0, 0, 1,
5452 & 1, 29*0/
5453
5454 DATA HINT1/100*0/
5455 DATA IHNT2/50*0/
5456
5457
5458 DATA NATT/0/EATT/0.0/JATT/0/NT/0/NP/0/N0/0/N01/0/N10/0/N11/0/
5459 DATA KATT/520000*0/PATT/520000*0.0/
5460
5461 DATA NFP/4500*0/PP/4500*0.0/NFT/4500*0/PT/4500*0.0/
5462
5463 DATA YP/900*0.0/YT/900*0.0/
5464
5465 DATA NPJ/300*0/KFPJ/150000*0/PJPX/150000*0.0/PJPY/150000*0.0/
5466 & PJPZ/150000*0.0/PJPE/150000*0.0/PJPM/150000*0.0/
5467 DATA NTJ/300*0/KFTJ/150000*0/PJTX/150000*0.0/PJTY/150000*0.0/
5468 & PJTZ/150000*0.0/PJTE/150000*0.0/PJTM/150000*0.0/
5469
5470 DATA NSG/0/NJSG/900*0/IASG/2700*0/K1SG/90000*0/K2SG/90000*0/
5471 & PXSG/90000*0.0/PYSG/90000*0.0/PZSG/90000*0.0/PESG/90000*0.0/
5472 & PMSG/90000*0.0/
5473 DATA MINT4/0/MINT5/0/ATCO/4000*0.0/ATXS/201*0.0/
5474 DATA (HIDAT0(1,I),I=1,10)/0.0,0.0,0.0,0.0,0.0,0.0,2.25,
5475 & 2.5,4.0,4.1/
5476 DATA (HIDAT0(2,I),I=1,10)/2.0,3.0,5.0,6.0,7.0,8.0,8.0,10.0,
5477 & 10.0,10.0/
5478 DATA (HIDAT0(3,I),I=1,10)/1.0,0.8,0.8,0.7,0.45,0.215,
5479 & 0.21,0.19,0.19,0.19/
5480 DATA (HIDAT0(4,I),I=1,10)/0.35,0.35,0.3,0.3,0.3,0.3,
5481 & 0.5,0.6,0.6,0.6/
5482 DATA (HIDAT0(5,I),I=1,10)/23.8,24.0,26.0,26.2,27.0,28.5,28.5,
5483 & 28.5,28.5,28.5/
5484 DATA ((HIDAT0(J,I),I=1,10),J=6,9)/40*0.0/
5485 DATA (HIDAT0(10,I),I=1,10)/5.0,20.0,53.0,62.0,100.0,200.0,
5486 & 546.0,900.0,1800.0,4000.0/
5487 DATA HIDAT/10*0.0/
5488 END
5489
5490
5491
5492
5493
5494
5495
5496
5497
5498
5499 SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A)
5500 IMPLICIT REAL*8(A-H,O-Z)
5501 COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
5502 COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
5503 COMMON/BVEG3/F,TI,TSI
5504 EXTERNAL FXN
5505 DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10)
5506 1 ,KG(10),IA(10)
5507 REAL*4 QRAN(10)
5508 SAVE
5509 DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/
5510
5511 NDO=1
5512 DO 1 J=1,NDIM
5513 1 XI(1,J)=ONE
5514
5515 ENTRY VEGAS1(FXN,AVGI,SD,CHI2A)
5516
5517 IT=0
5518 SI=0.
5519 SI2=SI
5520 SWGT=SI
5521 SCHI=SI
5522
5523 ENTRY VEGAS2(FXN,AVGI,SD,CHI2A)
5524
5525 ND=NDMX
5526 NG=1
5527 IF(MDS.EQ.0) GO TO 2
5528 NG=(NCALL/2.)**(1./NDIM)
5529 MDS=1
5530 IF((2*NG-NDMX).LT.0) GO TO 2
5531 MDS=-1
5532 NPG=NG/NDMX+1
5533 ND=NG/NPG
5534 NG=NPG*ND
5535 2 K=NG**NDIM
5536 NPG=NCALL/K
5537 IF(NPG.LT.2) NPG=2
5538 CALLS=NPG*K
5539 DXG=ONE/NG
5540 DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
5541 XND=ND
5542 NDM=ND-1
5543 DXG=DXG*XND
5544 XJAC=ONE/CALLS
5545 DO 3 J=1,NDIM
5546
5547 DX(J)=XU(J)-XL(J)
5548 3 XJAC=XJAC*DX(J)
5549
5550
5551
5552 IF(ND.EQ.NDO) GO TO 8
5553 RC=NDO/XND
5554 DO 7 J=1,NDIM
5555 K=0
5556 XN=0.
5557 DR=XN
5558 I=K
5559 4 K=K+1
5560 DR=DR+ONE
5561 XO=XN
5562 XN=XI(K,J)
5563 5 IF(RC.GT.DR) GO TO 4
5564 I=I+1
5565 DR=DR-RC
5566 XIN(I)=XN-(XN-XO)*DR
5567 IF(I.LT.NDM) GO TO 5
5568 DO 6 I=1,NDM
5569 6 XI(I,J)=XIN(I)
5570 7 XI(ND,J)=ONE
5571 NDO=ND
5572
5573 8 CONTINUE
5574 IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND
5575 1 ,(XL(J),XU(J),J=1,NDIM)
5576
5577 ENTRY VEGAS3(FXN,AVGI,SD,CHI2A)
5578
5579 9 IT=IT+1
5580 TI=0.
5581 TSI=TI
5582 DO 10 J=1,NDIM
5583 KG(J)=1
5584 DO 10 I=1,ND
5585 D(I,J)=TI
5586 10 DI(I,J)=TI
5587
5588 11 FB=0.
5589 F2B=FB
5590 K=0
5591 12 K=K+1
5592 CALL ARAN9(QRAN,NDIM)
5593 WGT=XJAC
5594 DO 15 J=1,NDIM
5595 XN=(KG(J)-QRAN(J))*DXG+ONE
5596
5597 IA(J)=XN
5598 IF(IA(J).GT.1) GO TO 13
5599 XO=XI(IA(J),J)
5600 RC=(XN-IA(J))*XO
5601 GO TO 14
5602 13 XO=XI(IA(J),J)-XI(IA(J)-1,J)
5603 RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
5604 14 X(J)=XL(J)+RC*DX(J)
5605 WGT=WGT*XO*XND
5606 15 CONTINUE
5607
5608 F=WGT
5609 F=F*FXN(X,WGT)
5610 F2=F*F
5611 FB=FB+F
5612 F2B=F2B+F2
5613 DO 16 J=1,NDIM
5614 DI(IA(J),J)=DI(IA(J),J)+F
5615 16 IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
5616 IF(K.LT.NPG) GO TO 12
5617
5618 F2B=DSQRT(F2B*NPG)
5619 F2B=(F2B-FB)*(F2B+FB)
5620 TI=TI+FB
5621 TSI=TSI+F2B
5622 IF(MDS.GE.0) GO TO 18
5623 DO 17 J=1,NDIM
5624 17 D(IA(J),J)=D(IA(J),J)+F2B
5625 18 K=NDIM
5626 19 KG(K)=MOD(KG(K),NG)+1
5627 IF(KG(K).NE.1) GO TO 11
5628 K=K-1
5629 IF(K.GT.0) GO TO 19
5630
5631
5632
5633 TSI=TSI*DV2G
5634 TI2=TI*TI
5635 WGT=TI2/(TSI+1.0d-37)
5636 SI=SI+TI*WGT
5637 SI2=SI2+TI2
5638 SWGT=SWGT+WGT
5639 SWGT=SWGT+1.0D-37
5640 SI2=SI2+1.0D-37
5641 SCHI=SCHI+TI2*WGT
5642 AVGI=SI/(SWGT)
5643 SD=SWGT*IT/(SI2)
5644 CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-.999)
5645 SD=DSQRT(ONE/SD)
5646
5647 IF(NPRN.EQ.0) GO TO 21
5648 TSI=DSQRT(TSI)
5649 WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A
5650 IF(NPRN.GE.0) GO TO 21
5651 DO 20 J=1,NDIM
5652 20 WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND)
5653
5654
5655
5656 21 DO 23 J=1,NDIM
5657 XO=D(1,J)
5658 XN=D(2,J)
5659 D(1,J)=(XO+XN)/2.
5660 DT(J)=D(1,J)
5661 DO 22 I=2,NDM
5662 D(I,J)=XO+XN
5663 XO=XN
5664 XN=D(I+1,J)
5665 D(I,J)=(D(I,J)+XN)/3.
5666 22 DT(J)=DT(J)+D(I,J)
5667 D(ND,J)=(XN+XO)/2.
5668 23 DT(J)=DT(J)+D(ND,J)
5669
5670 DO 28 J=1,NDIM
5671 RC=0.
5672 DO 24 I=1,ND
5673 R(I)=0.
5674 IF (DT(J).GE.1.0D18) THEN
5675 WRITE(6,*) '************** A SINGULARITY >1.0D18'
5676
5677
5678
5679
5680
5681
5682
5683
5684 END IF
5685 IF(D(I,J).LE.1.0D-18) GO TO 24
5686 XO=DT(J)/D(I,J)
5687 R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH
5688 24 RC=RC+R(I)
5689 RC=RC/XND
5690 K=0
5691 XN=0.
5692 DR=XN
5693 I=K
5694 25 K=K+1
5695 DR=DR+R(K)
5696 XO=XN
5697
5698 XN=XI(K,J)
5699 26 IF(RC.GT.DR) GO TO 25
5700 I=I+1
5701 DR=DR-RC
5702 XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30)
5703 IF(I.LT.NDM) GO TO 26
5704 DO 27 I=1,NDM
5705 27 XI(I,J)=XIN(I)
5706 28 XI(ND,J)=ONE
5707
5708 IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9
5709 200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',I3,' NCALL=',F8.0
5710 1 /28X,' IT=',I5,' ITMX=',I5/28X,' ACC=',G9.3
5711 2 /28X,' MDS=',I3,' ND=',I4/28X,' (XL,XU)=',
5712 3 (T40,'( ',G12.6,' , ',G12.6,' )'))
5713 201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3,
5714 1 ': INTEGRAL =',G14.8/21X,'STD DEV =',G10.4 /
5715 2 ' ACCUMULATED RESULTS: INTEGRAL =',G14.8 /
5716 3 24X,'STD DEV =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4)
5717 202 FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,' DELT I ',
5718 1 2X,' CONV''CE ',11X,'X',7X,' DELT I ',2X,' CONV''CE '
5719 2 ,11X,'X',7X,' DELT I ',2X,' CONV''CE ' /
5720 2 (' ',3G12.4,5X,3G12.4,5X,3G12.4))
5721 RETURN
5722 END
5723
5724
5725 SUBROUTINE ARAN9(QRAN,NDIM)
5726 DIMENSION QRAN(10)
5727 COMMON/SEEDVAX/NUM1
5728 DO 1 I=1,NDIM
5729 1 QRAN(I)=RLU(NUM1)
5730 RETURN
5731 END
5732
5733
5734
5735
5736
5737 FUNCTION GAUSS1(F,A,B,EPS)
5738 EXTERNAL F
5739 DIMENSION W(12),X(12)
5740 DATA CONST/1.0E-12/
5741 DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
5742 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5743 & .1826034,.1894506/
5744 DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
5745 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5746 & .2816036,.0950125/
5747 DELTA=CONST*ABS(A-B)
5748 GAUSS1=0.0
5749 AA=A
5750 5 Y=B-AA
5751 IF(ABS(Y).LE.DELTA) RETURN
5752 2 BB=AA+Y
5753 C1=0.5*(AA+BB)
5754 C2=C1-AA
5755 S8=0.0
5756 S16=0.0
5757 DO 1 I=1,4
5758 U=X(I)*C2
5759 1 S8=S8+W(I)*(F(C1+U)+F(C1-U))
5760 DO 3 I=5,12
5761 U=X(I)*C2
5762 3 S16=S16+W(I)*(F(C1+U)+F(C1-U))
5763 S8=S8*C2
5764 S16=S16*C2
5765 IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
5766 GAUSS1=GAUSS1+S16
5767 AA=BB
5768 GOTO 5
5769 4 Y=0.5*Y
5770 IF(ABS(Y).GT.DELTA) GOTO 2
5771 WRITE(6,7)
5772 GAUSS1=0.0
5773 RETURN
5774 7 FORMAT(1X,'GAUSS1....TOO HIGH ACURACY REQUIRED')
5775 END
5776
5777
5778
5779 FUNCTION GAUSS2(F,A,B,EPS)
5780 EXTERNAL F
5781 DIMENSION W(12),X(12)
5782 DATA CONST/1.0E-12/
5783 DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
5784 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5785 & .1826034,.1894506/
5786 DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
5787 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5788 & .2816036,.0950125/
5789 DELTA=CONST*ABS(A-B)
5790 GAUSS2=0.0
5791 AA=A
5792 5 Y=B-AA
5793 IF(ABS(Y).LE.DELTA) RETURN
5794 2 BB=AA+Y
5795 C1=0.5*(AA+BB)
5796 C2=C1-AA
5797 S8=0.0
5798 S16=0.0
5799 DO 1 I=1,4
5800 U=X(I)*C2
5801 1 S8=S8+W(I)*(F(C1+U)+F(C1-U))
5802 DO 3 I=5,12
5803 U=X(I)*C2
5804 3 S16=S16+W(I)*(F(C1+U)+F(C1-U))
5805 S8=S8*C2
5806 S16=S16*C2
5807 IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
5808 GAUSS2=GAUSS2+S16
5809 AA=BB
5810 GOTO 5
5811 4 Y=0.5*Y
5812 IF(ABS(Y).GT.DELTA) GOTO 2
5813 WRITE(6,7)
5814 GAUSS2=0.0
5815 RETURN
5816 7 FORMAT(1X,'GAUSS2....TOO HIGH ACURACY REQUIRED')
5817 END
5818
5819
5820
5821 FUNCTION GAUSS3(F,A,B,EPS)
5822 EXTERNAL F
5823 DIMENSION W(12),X(12)
5824 DATA CONST/1.0E-12/
5825 DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
5826 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5827 & .1826034,.1894506/
5828 DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
5829 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5830 & .2816036,.0950125/
5831 DELTA=CONST*ABS(A-B)
5832 GAUSS3=0.0
5833 AA=A
5834 5 Y=B-AA
5835 IF(ABS(Y).LE.DELTA) RETURN
5836 2 BB=AA+Y
5837 C1=0.5*(AA+BB)
5838 C2=C1-AA
5839 S8=0.0
5840 S16=0.0
5841 DO 1 I=1,4
5842 U=X(I)*C2
5843 1 S8=S8+W(I)*(F(C1+U)+F(C1-U))
5844 DO 3 I=5,12
5845 U=X(I)*C2
5846 3 S16=S16+W(I)*(F(C1+U)+F(C1-U))
5847 S8=S8*C2
5848 S16=S16*C2
5849 IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
5850 GAUSS3=GAUSS3+S16
5851 AA=BB
5852 GOTO 5
5853 4 Y=0.5*Y
5854 IF(ABS(Y).GT.DELTA) GOTO 2
5855 WRITE(6,7)
5856 GAUSS3=0.0
5857 RETURN
5858 7 FORMAT(1X,'GAUSS3....TOO HIGH ACURACY REQUIRED')
5859 END
5860
5861
5862
5863
5864 FUNCTION GAUSS4(F,A,B,EPS)
5865 EXTERNAL F
5866 DIMENSION W(12),X(12)
5867 DATA CONST/1.0E-12/
5868 DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
5869 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5870 & .1826034,.1894506/
5871 DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
5872 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5873 & .2816036,.0950125/
5874 DELTA=CONST*ABS(A-B)
5875 GAUSS4=0.0
5876 AA=A
5877 5 Y=B-AA
5878 IF(ABS(Y).LE.DELTA) RETURN
5879 2 BB=AA+Y
5880 C1=0.5*(AA+BB)
5881 C2=C1-AA
5882 S8=0.0
5883 S16=0.0
5884 DO 1 I=1,4
5885 U=X(I)*C2
5886 1 S8=S8+W(I)*(F(C1+U)+F(C1-U))
5887 DO 3 I=5,12
5888 U=X(I)*C2
5889 3 S16=S16+W(I)*(F(C1+U)+F(C1-U))
5890 S8=S8*C2
5891 S16=S16*C2
5892 IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
5893 GAUSS4=GAUSS4+S16
5894 AA=BB
5895 GOTO 5
5896 4 Y=0.5*Y
5897 IF(ABS(Y).GT.DELTA) GOTO 2
5898 WRITE(6,7)
5899 GAUSS4=0.0
5900 RETURN
5901 7 FORMAT(1X,'GAUSS4....TOO HIGH ACURACY REQUIRED')
5902 END
5903
5904
5905
5906
5907
5908 SUBROUTINE TITLE
5909 WRITE(6,200)
5910 200 FORMAT(//10X,
5911 & '**************************************************'/10X,
5912 & '* | \\ _______ / ------/ *'/10X,
5913 & '* ----- ------ |_____| /_/ / *'/10X,
5914 & '* ||\\ / |_____| / / \\ *'/10X,
5915 & '* /| \\ /_/ /_______ /_ / \\_ *'/10X,
5916 & '* / | / / / / / | ------- *'/10X,
5917 & '* | / /\\ / / | / | *'/10X,
5918 & '* | / / \\ / / \\_| / ------- *'/10X,
5919 & '* *'/10X,
5920 & '**************************************************'/10X,
5921 & ' HIJING '/10X,
5922 & ' Heavy Ion Jet INteraction Generator '/10X,
5923 & ' by '/10X,
5924 & ' X. N. Wang and M. Gyulassy '/10X,
5925 & ' Lawrence Berkeley Laboratory '//)
5926 RETURN
5927 END