196 dimension scip(300,300),rnip(300,300),sjip(300,300),jtp(3),
197 & ipcol(90000),itcol(90000)
198 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
205 common/himain1/natt,eatt,jatt,nt,np,n0,n01,n10,n11
207 common/himain2/katt(130000,4),patt(130000,4)
209 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
211 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
212 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
213 & pjpm(300,500),ntj(300),kftj(300,500),
214 & pjtx(300,500),pjty(300,500),pjtz(300,500),
215 & pjte(300,500),pjtm(300,500)
217 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
218 & k2sg(900,100),pxsg(900,100),pysg(900,100),
219 & pzsg(900,100),pesg(900,100),pmsg(900,100)
221 common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5)
226 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
228 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
231 bmax=min(bmax0,hipr1(34)+hipr1(35))
233 IF(ihnt2(1).LE.1 .AND. ihnt2(3).LE.1)
THEN
235 bmax=2.5*
sqrt(hipr1(31)*0.1/hipr1(40))
247 IF(ihnt2(1).LE.1) go to 14
251 if(ihnt2(1).EQ.2)
then
252 rnd1=max(
rlu(0),1.0
e-20)
253 rnd2=max(
rlu(0),1.0
e-20)
254 rnd3=max(
rlu(0),1.0
e-20)
255 r=-0.5*(
log(rnd1)*4.38/2.0+
log(rnd2)*0.85/2.0
256 & +4.38*0.85*
log(rnd3)/(4.38+0.85))
268 IF(hipr1(29).EQ.0.0) go to 10
270 dnbp1=(yp(1,kp)-yp(1,kp2))**2
271 dnbp2=(yp(2,kp)-yp(2,kp2))**2
272 dnbp3=(yp(3,kp)-yp(3,kp2))**2
273 dnbp=dnbp1+dnbp2+dnbp3
274 IF(dnbp.LT.hipr1(29)*hipr1(29)) go to 5
280 if(ihnt2(1).EQ.2)
then
288 IF(yp(3,i).GT.yp(3,j)) go to 12
304 IF(ihnt2(3).LE.1) go to 24
308 if(ihnt2(3).EQ.2)
then
309 rnd1=max(
rlu(0),1.0
e-20)
310 rnd2=max(
rlu(0),1.0
e-20)
311 rnd3=max(
rlu(0),1.0
e-20)
312 r=-0.5*(
log(rnd1)*4.38/2.0+
log(rnd2)*0.85/2.0
313 & +4.38*0.85*
log(rnd3)/(4.38+0.85))
325 IF(hipr1(29).EQ.0.0) go to 20
327 dnbt1=(
yt(1,kt)-
yt(1,kt2))**2
328 dnbt2=(
yt(2,kt)-
yt(2,kt2))**2
329 dnbt3=(
yt(3,kt)-
yt(3,kt2))**2
330 dnbt=dnbt1+dnbt2+dnbt3
331 IF(dnbt.LT.hipr1(29)*hipr1(29)) go to 15
337 if(ihnt2(3).EQ.2)
then
345 IF(
yt(3,i).LT.
yt(3,j)) go to 22
361 WRITE(6,*)
'infinite loop happened in HIJING'
389 bb=
sqrt(bmin**2+
rlu(0)*(bmax**2-bmin**2))
400 b2=(yp(1,jp)+bbx-
yt(1,jt))**2+(yp(2,jp)+bby-
yt(2,jt))**2
401 r2=b2*hipr1(40)/hipr1(31)/0.1
403 rrb1=min((yp(1,jp)**2+yp(2,jp)**2)
404 & /1.2**2/
REAL(ihnt2(1))**0.6666667,1.0)
405 rrb2=min((
yt(1,jt)**2+
yt(2,jt)**2)
406 & /1.2**2/
REAL(ihnt2(3))**0.6666667,1.0)
407 aphx1=hipr1(6)*4.0/3.0*(ihnt2(1)**0.3333333-1.0)
409 aphx2=hipr1(6)*4.0/3.0*(ihnt2(3)**0.3333333-1.0)
411 hint1(18)=hint1(14)-aphx1*hint1(15)
412 & -aphx2*hint1(16)+aphx1*aphx2*hint1(17)
414 IF(ihpr2(14).EQ.0.OR.
415 & (ihnt2(1).EQ.1.AND.ihnt2(3).EQ.1))
THEN
417 gs=1.0-
exp(-(hipr1(30)+hint1(18))*
romg(r2)/hipr1(31))
420 IF(rantot.GT.gs) go to 70
424 gstot_0=2.0*(1.0-
exp(-(hipr1(30)+hint1(18))
425 & /hipr1(31)/2.0*
romg(0.0)))
428 gs=1.0-
exp(-(hipr1(30)+hint1(18))/hipr1(31)*
romg(r2))
429 gstot=2.0*(1.0-
sqrt(1.0-gs))
430 rantot=
rlu(0)*gstot_0
432 IF(rantot.GT.gstot) go to 70
433 IF(rantot.GT.gs)
THEN
442 sjip(jp,jt)=hint1(18)
452 & (ihnt2(1).EQ.1.AND.ihnt2(3).EQ.1)) go to 60
459 IF(ihpr2(3).NE.0)
THEN
460 nhard=1+
int(
rlu(0)*(ncolt-1)+0.5)
461 nhard=min(nhard,ncolt)
466 IF(ihpr2(9).EQ.1)
THEN
467 nmini=1+
int(
rlu(0)*(ncolt-1)+0.5)
468 nmini=min(nmini,ncolt)
477 IF(scip(jp,jt).EQ.-1.0) go to 200
478 nfp(jp,11)=nfp(jp,11)+1
479 nft(jt,11)=nft(jt,11)+1
480 IF(nfp(jp,5).LE.1 .AND. nft(jt,5).GT.1)
THEN
483 ELSE IF(nfp(jp,5).GT.1 .AND. nft(jt,5).LE.1)
THEN
486 ELSE IF(nfp(jp,5).LE.1 .AND. nft(jt,5).LE.1)
THEN
490 ELSE IF(nfp(jp,5).GT.1 .AND. nft(jt,5).GT.1)
THEN
498 IF(ihpr2(8).EQ.0 .AND. ihpr2(3).EQ.0) go to 160
500 IF(nfp(jp,6).LT.0 .OR. nft(jt,6).LT.0) go to 160
505 hint1(18)=sjip(jp,jt)
506 tt=
romg(r2)*hint1(18)/hipr1(31)
507 tts=hipr1(30)*
romg(r2)/hipr1(31)
509 IF(ihpr2(3).NE.0 .AND. jp.EQ.jphard .AND. jt.EQ.jthard)
THEN
511 CALL
hijhrd(jp,jt,0,jflg,0)
521 IF(abs(hint1(46)).GT.hipr1(11).AND.jflg.EQ.2) nfp(jp,7)=1
522 IF(abs(hint1(56)).GT.hipr1(11).AND.jflg.EQ.2) nft(jt,7)=1
523 IF(max(abs(hint1(46)),abs(hint1(56))).GT.hipr1(11).AND.
524 & jflg.GE.3) iasg(nsg,3)=1
528 hint1(20+i05)=hint1(40+i05)
529 hint1(30+i05)=hint1(50+i05)
532 IF(ihpr2(8).EQ.0) go to 160
533 rrb1=min((yp(1,jp)**2+yp(2,jp)**2)/1.2**2
534 & /
REAL(ihnt2(1))**0.6666667,1.0)
535 rrb2=min((
yt(1,jt)**2+
yt(2,jt)**2)/1.2**2
536 & /
REAL(ihnt2(3))**0.6666667,1.0)
537 aphx1=hipr1(6)*4.0/3.0*(ihnt2(1)**0.3333333-1.0)
539 aphx2=hipr1(6)*4.0/3.0*(ihnt2(3)**0.3333333-1.0)
541 hint1(65)=hint1(61)-aphx1*hint1(62)
542 & -aphx2*hint1(63)+aphx1*aphx2*hint1(64)
543 ttrig=
romg(r2)*hint1(65)/hipr1(31)
548 xr1=-alog(
exp(-ttrig)+
rlu(0)*(1.0-
exp(-ttrig)))
550 xr1=xr1-alog(max(
rlu(0),1.0
e-20))
551 IF(xr1.LT.ttrig) go to 106
554 xr=xr-alog(max(
rlu(0),1.0
e-20))
555 IF(xr.LT.
tt-ttrig) go to 107
561 IF(ihpr2(9).EQ.1.AND.jp.EQ.jpmini.AND.jt.EQ.jtmini) go to 110
565 IF(ihpr2(8).GT.0 .AND.rnip(jp,jt).LT.
exp(-
tt)*
566 & (1.0-
exp(-tts))) go to 160
570 xr=xr-alog(max(
rlu(0),1.0
e-20))
571 IF(xr.LT.
tt) go to 111
572 112 njet=min(njet,ihpr2(8))
573 IF(ihpr2(8).LT.0) njet=abs(ihpr2(8))
578 CALL
hijhrd(jp,jt,jout,jflg,1)
585 IF(jflg.EQ.0) go to 160
587 IF(ihpr2(10).NE.0)
WRITE(6,*)
'error occured in HIJHRD'
591 IF(abs(hint1(46)).GT.hipr1(11).AND.jflg.EQ.2) nfp(jp,7)=1
592 IF(abs(hint1(56)).GT.hipr1(11).AND.jflg.EQ.2) nft(jt,7)=1
593 IF(max(abs(hint1(46)),abs(hint1(56))).GT.hipr1(11).AND.
594 & jflg.GE.3) iasg(nsg,3)=1
598 CALL
hijsft(jp,jt,jout,ierror)
600 IF(ihpr2(10).NE.0)
WRITE(6,*)
'error occured in HIJSFT'
612 IF(nfp(jp,5).GT.2)
THEN
614 ELSE IF(nfp(jp,5).EQ.2.OR.nfp(jp,5).EQ.1)
THEN
619 IF(nft(jt,5).GT.2)
THEN
621 ELSE IF(nft(jt,5).EQ.2.OR.nft(jt,5).EQ.1)
THEN
631 IF((ihpr2(8).NE.0.OR.ihpr2(3).NE.0).AND.ihpr2(4).GT.0.AND.
632 & ihnt2(1).GT.1.AND.ihnt2(3).GT.1)
THEN
634 IF(nfp(i,7).EQ.1) CALL
quench(i,1)
637 IF(nft(i,7).EQ.1) CALL
quench(i,2)
640 IF(iasg(isg,3).EQ.1) CALL
quench(isg,3)
651 IF(ihpr2(20).NE.0)
THEN
654 IF(mstu(24).NE.0 .OR.ierror.GT.0)
THEN
657 IF(ihpr2(10).NE.0)
THEN
659 WRITE(6,*)
'error occured, repeat the event'
667 IF(ihpr2(21).EQ.0)
THEN
671 IF(k(n_st,2).LT.91.OR.k(n_st,2).GT.93) go to 351
676 IF(frame.EQ.
'LAB')
THEN
683 IF(k(i,2).EQ.idstr)
THEN
692 IF(k(i,3).EQ.0 .OR. k(k(i,3),2).EQ.idstr)
THEN
695 katt(natt,3)=natt-i+k(i,3)+n_str-k(k(i,3),4)
709 DO 400 j_jtp=1,jtp(ntp)
710 CALL
hijfrg(j_jtp,ntp,ierror)
711 IF(mstu(24).NE.0 .OR. ierror.GT.0)
THEN
714 IF(ihpr2(10).NE.0)
THEN
716 WRITE(6,*)
'error occured, repeat the event'
724 IF(ihpr2(21).EQ.0)
THEN
728 IF(k(n_st,2).LT.91.OR.k(n_st,2).GT.93) go to 381
732 IF(frame.EQ.
'LAB')
THEN
738 IF(ntp.EQ.2) nftp=10+nft(j_jtp,5)
741 IF(k(i,2).EQ.idstr)
THEN
750 IF(k(i,3).EQ.0 .OR. k(k(i,3),2).EQ.idstr)
THEN
753 katt(natt,3)=natt-i+k(i,3)+n_str-k(k(i,3),4)
771 patt(natt,1)=pdr(i,1)
772 patt(natt,2)=pdr(i,2)
773 patt(natt,3)=pdr(i,3)
774 patt(natt,4)=pdr(i,4)
779 dengy=eatt/(ihnt2(1)*hint1(6)+ihnt2(3)*hint1(7))-1.0
780 IF(abs(dengy).GT.hipr1(43).AND.ihpr2(20).NE.0
781 & .AND.ihpr2(21).EQ.0)
THEN
782 IF(ihpr2(10).NE.0)
THEN
783 WRITE(6,*)
'Energy not conserved, repeat the event'
798 CHARACTER frame*4,proj*4,targ*4,eframe*4
799 DOUBLE PRECISION dd1,dd2,dd3,dd4
800 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
804 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
806 common/hijdat/hidat0(10,10),hidat(10)
808 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
886 IF(ihpr2(12).GT.0)
THEN
887 CALL
lugive(
'MDCY(C111,1)=0')
888 CALL
lugive(
'MDCY(C310,1)=0')
889 CALL
lugive(
'MDCY(C411,1)=0;MDCY(C-411,1)=0')
890 CALL
lugive(
'MDCY(C421,1)=0;MDCY(C-421,1)=0')
891 CALL
lugive(
'MDCY(C431,1)=0;MDCY(C-431,1)=0')
892 CALL
lugive(
'MDCY(C511,1)=0;MDCY(C-511,1)=0')
893 CALL
lugive(
'MDCY(C521,1)=0;MDCY(C-521,1)=0')
894 CALL
lugive(
'MDCY(C531,1)=0;MDCY(C-531,1)=0')
895 CALL
lugive(
'MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
896 CALL
lugive(
'MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
897 CALL
lugive(
'MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
898 CALL
lugive(
'MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
899 CALL
lugive(
'MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
900 CALL
lugive(
'MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
901 CALL
lugive(
'MDCY(C3334,1)=0;MDCY(C-3334,1)=0')
905 IF(ihpr2(10).EQ.0)
THEN
915 IF(frame.EQ.
'LAB')
THEN
919 hint1(1)=
sqrt(hint1(8)**2+2.0*hint1(9)*efrm+hint1(9)**2)
920 dd4=dsqrt(dd1**2-dd2**2)/(dd1+dd3)
922 hint1(3)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
923 dd4=dsqrt(dd1**2-dd2**2)/dd1
924 hint1(4)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
928 ELSE IF(frame.EQ.
'CMS')
THEN
935 dd4=dsqrt(1.d0-4.d0*dd2**2/dd1**2)
936 hint1(4)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
937 dd4=dsqrt(1.d0-4.d0*dd3**2/dd1**2)
938 hint1(5)=-0.5*dlog((1.d0+dd4)/(1.d0-dd4))
939 hint1(6)=hint1(1)/2.0
940 hint1(7)=hint1(1)/2.0
946 IF(ihnt2(1).GT.1)
THEN
947 CALL
hijwds(ihnt2(1),1,rmax)
951 IF(ihnt2(3).GT.1)
THEN
952 CALL
hijwds(ihnt2(3),2,rmax)
963 IF(hidat0(10,i).LE.hint1(1)) go to 20
966 hidat(j)=hidat0(j,i-1)+(hidat0(j,i)-hidat0(j,i-1))
967 & *(hint1(1)-hidat0(10,i-1))/(hidat0(10,i)-hidat0(10,i-1))
970 hipr1(30)=2.0*hidat(5)
978 IF(ihpr2(5).NE.0)
THEN
988 IF(frame.EQ.
'LAB') eframe=
'Elab'
989 WRITE(6,100) eframe,efrm,proj,ihnt2(1),ihnt2(2),
990 & targ,ihnt2(3),ihnt2(4)
991 100
FORMAT(//10x,
'****************************************
994 & 10x,
'* HIJING has been initialized at *'/
995 & 10x,
'*',13x,a4,
'= ',f10.2,
' GeV/n',13x,
'*'/
998 & a4,
'(',i3,
',',i3,
')',
' + ',a4,
'(',i3,
',',i3,
')',7x,
'*'/
999 & 10x,
'**************************************************')
1006 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1008 fnkick=1.0/(x+hipr1(19)**2)/(x+hipr1(20)**2)
1009 & /(1+
exp((
sqrt(x)-hipr1(20))/0.4))
1015 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1024 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1026 fnstru=(1.0-x)**hipr1(44)/
1027 & (x**2+hipr1(45)**2/hint1(1)**2)**hipr1(46)
1034 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
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)
1043 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1046 & (x**2+hipr1(45)**2/hint1(1)**2)**hipr1(48)
1054 IMPLICIT DOUBLE PRECISION(
d)
1055 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
1057 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
1059 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1063 IF(abs(dbeta).GE.1.d0)
THEN
1065 IF(db.GT.0.99999999d0)
THEN
1067 WRITE(6,*)
'(HIBOOT:) boost vector too large'
1070 dga=1d0/
sqrt(1d0-db**2)
1073 p(i,3)=(dp3+db*dp4)*dga
1074 p(i,4)=(dp4+db*dp3)*dga
1077 y=0.5*dlog((1.d0+dbeta)/(1.d0-dbeta))
1078 amt=
sqrt(
p(i,1)**2+
p(i,2)**2+
p(i,5)**2)
1079 p(i,3)=amt*sinh(
y+hint1(3))
1080 p(i,4)=amt*cosh(
y+hint1(3))
1089 dimension rdp(300),lqp(300),rdt(300),lqt(300)
1092 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1095 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
1096 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
1097 & pjpm(300,500),ntj(300),kftj(300,500),
1098 & pjtx(300,500),pjty(300,500),pjtz(300,500),
1099 & pjte(300,500),pjtm(300,500)
1101 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
1102 & k2sg(900,100),pxsg(900,100),pysg(900,100),
1103 & pzsg(900,100),pesg(900,100),pmsg(900,100)
1105 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
1115 IF(ntp.EQ.2) go to 400
1116 IF(ntp.EQ.3) go to 2000
1121 IF(nfp(jpjt,7).NE.1)
RETURN
1125 ptjet0=
sqrt(pjpx(jp,i)**2+pjpy(jp,i)**2)
1126 IF(ptjet0.LE.hipr1(11)) go to 290
1127 ptot=
sqrt(ptjet0*ptjet0+pjpz(jp,i)**2)
1128 IF(ptot.LT.hipr1(8)) go to 290
1129 phip=
ulangl(pjpx(jp,i),pjpy(jp,i))
1132 DO 100 i2=1,ihnt2(1)
1133 IF(nfp(i2,5).NE.3 .OR. i2.EQ.jp) go to 100
1134 dx=yp(1,i2)-yp(1,jp)
1135 dy=yp(2,i2)-yp(2,jp)
1138 IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi
1139 IF(dphi.GE.hipr1(40)/2.0) go to 100
1141 IF(rd0*
sin(dphi).GT.hipr1(12)) go to 100
1144 rdp(kp)=
cos(dphi)*rd0
1149 IF(rdp(i2).LT.rdp(j2)) go to 110
1159 DO 120 i2=1,ihnt2(3)
1160 IF(nft(i2,5).NE.3) go to 120
1161 dx=
yt(1,i2)-yp(1,jp)-bbx
1162 dy=
yt(2,i2)-yp(2,jp)-bby
1165 IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi
1166 IF(dphi.GT.hipr1(40)/2.0) go to 120
1168 IF(rd0*
sin(dphi).GT.hipr1(12)) go to 120
1171 rdt(kt)=
cos(dphi)*rd0
1176 IF(rdt(i2).LT.rdt(j2)) go to 130
1190 ptot=
sqrt(pjpx(jp,i)**2+pjpy(jp,i)**2+pjpz(jp,i)**2)
1196 210
IF(
mt.GE.kt .AND. mp.GE.kp) go to 290
1197 IF(
mt.GE.kt) go to 220
1198 IF(mp.GE.kp) go to 240
1199 IF(rdp(mp+1).GT.rdt(
mt+1)) go to 240
1202 IF(rn.GE.1.0-
exp(-drr/hipr1(13))) go to 210
1204 IF(kfpj(jp,i).NE.21) dp=0.5*dp
1206 IF(dp.LE.0.2) go to 210
1207 IF(ptot.LE.0.4) go to 290
1208 IF(ptot.LE.dp) dp=ptot-0.2
1211 IF(kfpj(jp,i).NE.21)
THEN
1212 prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
1214 de=
sqrt(pjpm(jp,i)**2+ptot**2)
1215 & -
sqrt(pjpm(jp,i)**2+(ptot-dp)**2)
1216 ershu=(pp(lqp(mp),4)+de-dp)**2
1218 IF(amshu.LT.hipr1(1)*hipr1(1)) go to 210
1219 pp(lqp(mp),4)=
sqrt(ershu)
1220 pp(lqp(mp),5)=
sqrt(amshu)
1229 npj(lqp(mp))=npj(lqp(mp))+1
1230 kfpj(lqp(mp),npj(lqp(mp)))=21
1231 pjpx(lqp(mp),npj(lqp(mp)))=dp1
1232 pjpy(lqp(mp),npj(lqp(mp)))=dp2
1233 pjpz(lqp(mp),npj(lqp(mp)))=dp3
1234 pjpe(lqp(mp),npj(lqp(mp)))=dp
1235 pjpm(lqp(mp),npj(lqp(mp)))=0.0
1240 IF(rn.GE.1.0-
exp(-drr/hipr1(13))) go to 210
1242 IF(dp.LE.0.2) go to 210
1243 IF(ptot.LE.0.4) go to 290
1244 IF(ptot.LE.dp) dp=ptot-0.2
1247 IF(kfpj(jp,i).NE.21)
THEN
1248 prshu=pt(lqt(
mt),1)**2+pt(lqt(
mt),2)**2
1250 de=
sqrt(pjpm(jp,i)**2+ptot**2)
1251 & -
sqrt(pjpm(jp,i)**2+(ptot-dp)**2)
1252 ershu=(pt(lqt(
mt),4)+de-dp)**2
1254 IF(amshu.LT.hipr1(1)*hipr1(1)) go to 210
1255 pt(lqt(
mt),4)=
sqrt(ershu)
1256 pt(lqt(
mt),5)=
sqrt(amshu)
1265 ntj(lqt(
mt))=ntj(lqt(
mt))+1
1266 kftj(lqt(
mt),ntj(lqt(
mt)))=21
1267 pjtx(lqt(
mt),ntj(lqt(
mt)))=dp1
1268 pjty(lqt(
mt),ntj(lqt(
mt)))=dp2
1269 pjtz(lqt(
mt),ntj(lqt(
mt)))=dp3
1270 pjte(lqt(
mt),ntj(lqt(
mt)))=dp
1271 pjtm(lqt(
mt),ntj(lqt(
mt)))=0.0
1273 260 pjpx(jp,i)=(ptot-dp)*v1
1274 pjpy(jp,i)=(ptot-dp)*v2
1275 pjpz(jp,i)=(ptot-dp)*v3
1276 pjpe(jp,i)=pjpe(jp,i)-de
1291 400
IF(nft(jpjt,7).NE.1)
RETURN
1294 ptjet0=
sqrt(pjtx(jt,i)**2+pjty(jt,i)**2)
1295 IF(ptjet0.LE.hipr1(11)) go to 690
1296 ptot=
sqrt(ptjet0*ptjet0+pjtz(jt,i)**2)
1297 IF(ptot.LT.hipr1(8)) go to 690
1298 phit=
ulangl(pjtx(jt,i),pjty(jt,i))
1300 DO 500 i2=1,ihnt2(1)
1301 IF(nfp(i2,5).NE.3) go to 500
1302 dx=yp(1,i2)+bbx-
yt(1,jt)
1303 dy=yp(2,i2)+bby-
yt(2,jt)
1306 IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi
1307 IF(dphi.GT.hipr1(40)/2.0) go to 500
1309 IF(rd0*
sin(dphi).GT.hipr1(12)) go to 500
1312 rdp(kp)=
cos(dphi)*rd0
1317 IF(rdp(i2).LT.rdp(j2)) go to 510
1327 DO 520 i2=1,ihnt2(3)
1328 IF(nft(i2,5).NE.3 .OR. i2.EQ.jt) go to 520
1333 IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi
1334 IF(dphi.GT.hipr1(40)/2.0) go to 520
1336 IF(rd0*
sin(dphi).GT.hipr1(12)) go to 520
1339 rdt(kt)=
cos(dphi)*rd0
1344 IF(rdt(i2).LT.rdt(j2)) go to 530
1358 ptot=
sqrt(pjtx(jt,i)**2+pjty(jt,i)**2+pjtz(jt,i)**2)
1364 610
IF(
mt.GE.kt .AND. mp.GE.kp) go to 690
1365 IF(
mt.GE.kt) go to 620
1366 IF(mp.GE.kp) go to 640
1367 IF(rdp(mp+1).GT.rdt(
mt+1)) go to 640
1370 IF(rn.GE.1.0-
exp(-drr/hipr1(13))) go to 610
1372 IF(kftj(jt,i).NE.21) dp=0.5*dp
1374 IF(dp.LE.0.2) go to 610
1375 IF(ptot.LE.0.4) go to 690
1376 IF(ptot.LE.dp) dp=ptot-0.2
1379 IF(kftj(jt,i).NE.21)
THEN
1380 prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
1382 de=
sqrt(pjtm(jt,i)**2+ptot**2)
1383 & -
sqrt(pjtm(jt,i)**2+(ptot-dp)**2)
1384 ershu=(pp(lqp(mp),4)+de-dp)**2
1386 IF(amshu.LT.hipr1(1)*hipr1(1)) go to 610
1387 pp(lqp(mp),4)=
sqrt(ershu)
1388 pp(lqp(mp),5)=
sqrt(amshu)
1397 npj(lqp(mp))=npj(lqp(mp))+1
1398 kfpj(lqp(mp),npj(lqp(mp)))=21
1399 pjpx(lqp(mp),npj(lqp(mp)))=dp1
1400 pjpy(lqp(mp),npj(lqp(mp)))=dp2
1401 pjpz(lqp(mp),npj(lqp(mp)))=dp3
1402 pjpe(lqp(mp),npj(lqp(mp)))=dp
1403 pjpm(lqp(mp),npj(lqp(mp)))=0.0
1409 IF(rn.GE.1.0-
exp(-drr/hipr1(13))) go to 610
1411 IF(dp.LE.0.2) go to 610
1412 IF(ptot.LE.0.4) go to 690
1413 IF(ptot.LE.dp) dp=ptot-0.2
1416 IF(kftj(jt,i).NE.21)
THEN
1417 prshu=pt(lqt(
mt),1)**2+pt(lqt(
mt),2)**2
1419 de=
sqrt(pjtm(jt,i)**2+ptot**2)
1420 & -
sqrt(pjtm(jt,i)**2+(ptot-dp)**2)
1421 ershu=(pt(lqt(
mt),4)+de-dp)**2
1423 IF(amshu.LT.hipr1(1)*hipr1(1)) go to 610
1424 pt(lqt(
mt),4)=
sqrt(ershu)
1425 pt(lqt(
mt),5)=
sqrt(amshu)
1434 ntj(lqt(
mt))=ntj(lqt(
mt))+1
1435 kftj(lqt(
mt),ntj(lqt(
mt)))=21
1436 pjtx(lqt(
mt),ntj(lqt(
mt)))=dp1
1437 pjty(lqt(
mt),ntj(lqt(
mt)))=dp2
1438 pjtz(lqt(
mt),ntj(lqt(
mt)))=dp3
1439 pjte(lqt(
mt),ntj(lqt(
mt)))=dp
1440 pjtm(lqt(
mt),ntj(lqt(
mt)))=0.0
1442 660 pjtx(jt,i)=(ptot-dp)*v1
1443 pjty(jt,i)=(ptot-dp)*v2
1444 pjtz(jt,i)=(ptot-dp)*v3
1445 pjte(jt,i)=pjte(jt,i)-de
1456 IF(iasg(isg,3).NE.1)
RETURN
1460 xj=(yp(1,jp)+bbx+
yt(1,jt))/2.0
1461 yj=(yp(2,jp)+bby+
yt(2,jt))/2.0
1462 DO 2690 i=1,njsg(isg)
1463 ptjet0=
sqrt(pxsg(isg,i)**2+pysg(isg,i)**2)
1464 IF(ptjet0.LE.hipr1(11).OR.pesg(isg,i).LT.hipr1(1))
1466 ptot=
sqrt(ptjet0*ptjet0+pzsg(isg,i)**2)
1467 IF(ptot.LT.max(hipr1(1),hipr1(8))) go to 2690
1468 phiq=
ulangl(pxsg(isg,i),pysg(isg,i))
1470 DO 2500 i2=1,ihnt2(1)
1471 IF(nfp(i2,5).NE.3.OR.i2.EQ.jp) go to 2500
1476 IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi
1477 IF(dphi.GT.hipr1(40)/2.0) go to 2500
1479 IF(rd0*
sin(dphi).GT.hipr1(12)) go to 2500
1482 rdp(kp)=
cos(dphi)*rd0
1487 IF(rdp(i2).LT.rdp(j2)) go to 2510
1497 DO 2520 i2=1,ihnt2(3)
1498 IF(nft(i2,5).NE.3 .OR. i2.EQ.jt) go to 2520
1503 IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi
1504 IF(dphi.GT.hipr1(40)/2.0) go to 2520
1506 IF(rd0*
sin(dphi).GT.hipr1(12)) go to 2520
1509 rdt(kt)=
cos(dphi)*rd0
1514 IF(rdt(i2).LT.rdt(j2)) go to 2530
1528 ptot=
sqrt(pxsg(isg,i)**2+pysg(isg,i)**2
1535 2610
IF(
mt.GE.kt .AND. mp.GE.kp) go to 2690
1536 IF(
mt.GE.kt) go to 2620
1537 IF(mp.GE.kp) go to 2640
1538 IF(rdp(mp+1).GT.rdt(
mt+1)) go to 2640
1541 IF(rn.GE.1.0-
exp(-drr/hipr1(13))) go to 2610
1542 dp=drr*hipr1(14)/2.0
1543 IF(dp.LE.0.2) go to 2610
1544 IF(ptot.LE.0.4) go to 2690
1545 IF(ptot.LE.dp) dp=ptot-0.2
1548 IF(k2sg(isg,i).NE.21)
THEN
1549 IF(ptot.LT.dp+hipr1(1)) go to 2690
1550 prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
1552 de=
sqrt(pmsg(isg,i)**2+ptot**2)
1553 & -
sqrt(pmsg(isg,i)**2+(ptot-dp)**2)
1554 ershu=(pp(lqp(mp),4)+de-dp)**2
1556 IF(amshu.LT.hipr1(1)*hipr1(1)) go to 2610
1557 pp(lqp(mp),4)=
sqrt(ershu)
1558 pp(lqp(mp),5)=
sqrt(amshu)
1567 npj(lqp(mp))=npj(lqp(mp))+1
1568 kfpj(lqp(mp),npj(lqp(mp)))=21
1569 pjpx(lqp(mp),npj(lqp(mp)))=dp1
1570 pjpy(lqp(mp),npj(lqp(mp)))=dp2
1571 pjpz(lqp(mp),npj(lqp(mp)))=dp3
1572 pjpe(lqp(mp),npj(lqp(mp)))=dp
1573 pjpm(lqp(mp),npj(lqp(mp)))=0.0
1579 IF(rn.GE.1.0-
exp(-drr/hipr1(13))) go to 2610
1581 IF(dp.LE.0.2) go to 2610
1582 IF(ptot.LE.0.4) go to 2690
1583 IF(ptot.LE.dp) dp=ptot-0.2
1586 IF(k2sg(isg,i).NE.21)
THEN
1587 IF(ptot.LT.dp+hipr1(1)) go to 2690
1588 prshu=pt(lqt(
mt),1)**2+pt(lqt(
mt),2)**2
1590 de=
sqrt(pmsg(isg,i)**2+ptot**2)
1591 & -
sqrt(pmsg(isg,i)**2+(ptot-dp)**2)
1592 ershu=(pt(lqt(
mt),4)+de-dp)**2
1594 IF(amshu.LT.hipr1(1)*hipr1(1)) go to 2610
1595 pt(lqt(
mt),4)=
sqrt(ershu)
1596 pt(lqt(
mt),5)=
sqrt(amshu)
1605 ntj(lqt(
mt))=ntj(lqt(
mt))+1
1606 kftj(lqt(
mt),ntj(lqt(
mt)))=21
1607 pjtx(lqt(
mt),ntj(lqt(
mt)))=dp1
1608 pjty(lqt(
mt),ntj(lqt(
mt)))=dp2
1609 pjtz(lqt(
mt),ntj(lqt(
mt)))=dp3
1610 pjte(lqt(
mt),ntj(lqt(
mt)))=dp
1611 pjtm(lqt(
mt),ntj(lqt(
mt)))=0.0
1613 2660 pxsg(isg,i)=(ptot-dp)*v1
1614 pysg(isg,i)=(ptot-dp)*v2
1615 pzsg(isg,i)=(ptot-dp)*v3
1616 pesg(isg,i)=pesg(isg,i)-de
1636 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1638 common/hijdat/hidat0(10,10),hidat(10)
1640 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
1642 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
1643 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
1644 & pjpm(300,500),ntj(300),kftj(300,500),
1645 & pjtx(300,500),pjty(300,500),pjtz(300,500),
1646 & pjte(300,500),pjtm(300,500)
1648 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
1649 & k2sg(900,100),pxsg(900,100),pysg(900,100),
1650 & pzsg(900,100),pesg(900,100),pmsg(900,100)
1653 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
1655 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
1667 DO 100 i=1,njsg(isg)
1683 IF(ntp.EQ.2) go to 200
1684 IF(jtp.GT.ihnt2(1))
RETURN
1685 IF(nfp(jtp,5).NE.3.AND.nfp(jtp,3).NE.0
1686 & .AND.npj(jtp).EQ.0.AND.nfp(jtp,10).EQ.0) go to 1000
1687 IF(nfp(jtp,15).EQ.-1)
THEN
1712 IF((abs(pb1-pp(jtp,1)).GT.0.01.OR.
1713 & abs(pb2-pp(jtp,2)).GT.0.01).AND.ihpr2(10).NE.0)
1714 &
WRITE(6,*)
' Pt of Q and QQ do not sum to the total'
1718 200
IF(jtp.GT.ihnt2(3))
RETURN
1719 IF(nft(jtp,5).NE.3.AND.nft(jtp,3).NE.0
1720 & .AND.ntj(jtp).EQ.0.AND.nft(jtp,10).EQ.0) go to 1200
1721 IF(nft(jtp,15).EQ.1)
THEN
1747 IF((abs(pb1-pt(jtp,1)).GT.0.01.OR.
1748 & abs(pb2-pt(jtp,2)).GT.0.01).AND.ihpr2(10).NE.0)
1749 &
WRITE(6,*)
' Pt of Q and QQ do not sum to the total'
1751 300
IF(pecm.LT.hipr1(1))
THEN
1753 IF(ihpr2(10).EQ.0)
RETURN
1754 WRITE(6,*)
' ECM=',pecm,
' energy of the string is too small'
1757 amt=pecm**2+pb1**2+pb2**2
1758 amt1=am1**2+pq11**2+pq12**2
1759 amt2=am2**2+pq21**2+pq22**2
1760 pzcm=
sqrt(abs(amt**2+amt1**2+amt2**2-2.0*amt*amt1
1761 & -2.0*amt*amt2-2.0*amt1*amt2))/2.0/
sqrt(amt)
1768 p(1,4)=
sqrt(amt1+pzcm**2)
1775 p(2,4)=
sqrt(amt2+pzcm**2)
1779 CALL
hirobo(0.0,0.0,0.0,0.0,btz)
1781 IF((pq21**2+pq22**2).GT.(pq11**2+pq12**2))
THEN
1794 ELSE IF(ntp.EQ.2)
THEN
1800 IF(ntp.EQ.1.AND.npj(jtp).NE.0)
THEN
1805 IF((abs(kf1).GT.1000.AND.kf1.LT.0)
1806 & .OR.(abs(kf1).LT.1000.AND.kf1.GT.0)) iex=1
1820 IF(iex.EQ.1) i0=npj(jtp)-i+1
1825 IF(kk1.NE.21 .AND. kk1.NE.0) k(i+1,1)=
1826 & 1+(abs(kk1)+(2*iex-1)*kk1)/2/abs(kk1)
1827 p(i+1,1)=pjpx(jtp,i0)
1828 p(i+1,2)=pjpy(jtp,i0)
1829 p(i+1,3)=pjpz(jtp,i0)
1830 p(i+1,4)=pjpe(jtp,i0)
1831 p(i+1,5)=pjpm(jtp,i0)
1834 ELSE IF(ntp.EQ.2.AND.ntj(jtp).NE.0)
THEN
1839 IF((abs(kf2).GT.1000.AND.kf2.LT.0)
1840 & .OR.(abs(kf2).LT.1000.AND.kf2.GT.0)) iex=0
1854 IF(iex.EQ.1) i0=ntj(jtp)-i+1
1859 IF(kk1.NE.21 .AND. kk1.NE.0) k(i+1,1)=
1860 & 1+(abs(kk1)+(2*iex-1)*kk1)/2/abs(kk1)
1861 p(i+1,1)=pjtx(jtp,i0)
1862 p(i+1,2)=pjty(jtp,i0)
1863 p(i+1,3)=pjtz(jtp,i0)
1864 p(i+1,4)=pjte(jtp,i0)
1865 p(i+1,5)=pjtm(jtp,i0)
1869 IF(ihpr2(1).GT.0.AND.
rlu(0).LE.hidat(3))
THEN
1872 IF(ihpr2(8).EQ.0.AND.ihpr2(3).EQ.0.AND.ihpr2(9).EQ.0)
1874 IF(hint1(1).GE.1000.0.AND.jetot.EQ.0)
THEN
1881 ELSE IF(jetot.EQ.0.AND.ihpr2(1).GT.0.AND.
1882 & hint1(1).GE.1000.0.AND.
1883 &
rlu(0).LE.0.8)
THEN
1888 IF(ihpr2(8).EQ.0.AND.ihpr2(3).EQ.0.AND.ihpr2(9).EQ.0)
1894 IF(ierror.NE.0)
RETURN
1933 dimension kf(100),
px(100),
py(100),
pz(100),pe(100),pm(100)
1934 dimension
y(100),ip(100,2)
1935 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
1936 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
1937 & pjpm(300,500),ntj(300),kftj(300,500),
1938 & pjtx(300,500),pjty(300,500),pjtz(300,500),
1939 & pjte(300,500),pjtm(300,500)
1941 IF(npt.EQ.2) go to 500
1945 100 kf(i)=kfpj(jp,i)
1951 y(i-iq)=0.5*alog((abs(pe(i)+
pz(i))+1.
e-5)
1952 & /(abs(pe(i)-
pz(i))+1.
e-5))
1955 IF(kf(i).NE.21)
THEN
1967 IF(i.LE.npj(jp)) go to 100
1969 DO 200 i=1,npj(jp)-iq
1970 DO 200 j=i+1,npj(jp)-iq
1971 IF(
y(i).GT.
y(j)) go to 200
1982 300 kfpj(jp,i)=kf(ip(i-iqq,1))
1983 pjpx(jp,i)=
px(ip(i-iqq,1))
1984 pjpy(jp,i)=
py(ip(i-iqq,1))
1985 pjpz(jp,i)=
pz(ip(i-iqq,1))
1986 pjpe(jp,i)=pe(ip(i-iqq,1))
1987 pjpm(jp,i)=pm(ip(i-iqq,1))
1988 IF(ip(i-iqq,2).EQ.1)
THEN
1989 kfpj(jp,i+1)=kf(ip(i-iqq,1)+1)
1990 pjpx(jp,i+1)=
px(ip(i-iqq,1)+1)
1991 pjpy(jp,i+1)=
py(ip(i-iqq,1)+1)
1992 pjpz(jp,i+1)=
pz(ip(i-iqq,1)+1)
1993 pjpe(jp,i+1)=pe(ip(i-iqq,1)+1)
1994 pjpm(jp,i+1)=pm(ip(i-iqq,1)+1)
1999 IF(i.LE.npj(jp)) go to 300
2006 600 kf(i)=kftj(jt,i)
2012 y(i-iq)=0.5*alog((abs(pe(i)+
pz(i))+1.
e-5)
2013 & /(abs(pe(i)-
pz(i))+1.
e-5))
2016 IF(kf(i).NE.21)
THEN
2028 IF(i.LE.ntj(jt)) go to 600
2030 DO 700 i=1,ntj(jt)-iq
2031 DO 700 j=i+1,ntj(jt)-iq
2032 IF(
y(i).LT.
y(j)) go to 700
2043 800 kftj(jt,i)=kf(ip(i-iqq,1))
2044 pjtx(jt,i)=
px(ip(i-iqq,1))
2045 pjty(jt,i)=
py(ip(i-iqq,1))
2046 pjtz(jt,i)=
pz(ip(i-iqq,1))
2047 pjte(jt,i)=pe(ip(i-iqq,1))
2048 pjtm(jt,i)=pm(ip(i-iqq,1))
2049 IF(ip(i-iqq,2).EQ.1)
THEN
2050 kftj(jt,i+1)=kf(ip(i-iqq,1)+1)
2051 pjtx(jt,i+1)=
px(ip(i-iqq,1)+1)
2052 pjty(jt,i+1)=
py(ip(i-iqq,1)+1)
2053 pjtz(jt,i+1)=
pz(ip(i-iqq,1)+1)
2054 pjte(jt,i+1)=pe(ip(i-iqq,1)+1)
2055 pjtm(jt,i+1)=pm(ip(i-iqq,1)+1)
2060 IF(i.LE.ntj(jt)) go to 800
2072 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2074 common/hijdat/hidat0(10,10),hidat(10)
2076 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
2086 s=2.*(
p(i,4)*
p(i+1,4)-
p(i,1)*
p(i+1,1)-
p(i,2)*
p(i+1,2)
2087 & -
p(i,3)*
p(i+1,3))+
p(i,5)**2+
p(i+1,5)**2
2091 pbt1=
p(i,1)+
p(i+1,1)
2092 pbt2=
p(i,2)+
p(i+1,2)
2093 pbt3=
p(i,3)+
p(i+1,3)
2094 pbt4=
p(i,4)+
p(i+1,4)
2095 btt=(pbt1**2+pbt2**2+pbt3**2)/pbt4**2
2096 IF(btt.GE.1.0-1.0
e-10) go to 30
2097 IF((i.NE.1.OR.i.NE.
n-1).AND.
2098 & (k(i,2).NE.21.AND.k(i+1,2).NE.21)) go to 30
2103 s=(sm+1.5*(
p(jl,5)+
p(jl+1,5)))**2
2104 IF(sm.LT.hipr1(5)) goto 2
2107 IF(jl+1.EQ.
n) goto 190
2117 p1=
p(jl,1)+
p(jl+1,1)
2118 p2=
p(jl,2)+
p(jl+1,2)
2119 p3=
p(jl,3)+
p(jl+1,3)
2120 p4=
p(jl,4)+
p(jl+1,4)
2126 CALL
atrobo(0.,0.,bex,bey,bez,imin,
imax,ierror)
2127 IF(ierror.NE.0)
RETURN
2129 cth=
p(jl,3)/
sqrt(
p(jl,4)**2-
p(jl,5)**2)
2130 IF(abs(cth).GT.1.0) cth=max(-1.,min(1.,cth))
2140 IF(hidat(2).GT.0.0)
THEN
2141 ptg1=
sqrt(
p(jl,1)**2+
p(jl,2)**2)
2142 ptg2=
sqrt(
p(jl+1,1)**2+
p(jl+1,2)**2)
2143 ptg3=
sqrt(
p(jl+2,1)**2+
p(jl+2,2)**2)
2144 ptg=max(ptg1,ptg2,ptg3)
2145 IF(ptg.GT.hidat(2))
THEN
2146 fmfact=
exp(-(ptg**2-hidat(2)**2)/hipr1(2)**2)
2147 IF(
rlu(0).GT.fmfact) go to 1
2154 IF(ierror.NE.0)
RETURN
2179 IF(sm.GE.hipr1(5)) goto 40
2196 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2198 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
2204 IF(k(jl,2).NE.21 .AND. k(jl+1,2).NE.21) c=8./27.
2207 IF(k(jl,2).NE.21) exp1=2
2208 IF(k(jl+1,2).NE.21) exp3=2
2214 xt2m=(1.-2.*
sqrt(sm1)+sm1-sm3)*(1.-2.*
sqrt(sm3)-sm1+sm3)
2217 1
IF(ntry.EQ.5000)
THEN
2218 x1=.5*(2.*
sqrt(sm1)+1.+sm1-sm3)
2219 x3=.5*(2.*
sqrt(sm3)+1.-sm1+sm3)
2224 xt2=
a*(xt2m/
a)**(
rlu(0)**(1./
d))
2226 ymax=alog(.5/
sqrt(xt2)+
sqrt(.25/xt2-1.))
2227 y=(2.*
rlu(0)-1.)*ymax
2232 IF(k(jl,2).NE.21 .OR. k(jl+1,2).NE.21)
THEN
2233 IF((1.-x1)*(1.-x2)*(1.-x3)-x2*sm1*(1.-x1)-x2*sm3*(1.-x3).
2234 & le.0..OR.x1.LE.2.*
sqrt(sm1)-sm1+sm3.OR.x3.LE.2.*
sqrt(sm3)
2241 fg=2.*ymax*c*(x1**exp1+x3**exp3)/
d
2243 IF(fg.LT.
rlu(0)) goto 1
2252 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2254 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
2266 IF(
p1.GT.0..AND.
p3.GT.0.) cbet=(
p(jl,5)**2
2267 & +
p(jl+1,5)**2+2.*e1*e3-
s*(1.-x2))/(2.*
p1*
p3)
2268 IF(abs(cbet).GT.1.0) cbet=max(-1.,min(1.,cbet))
2278 ELSE IF(
p3.GT.
p1)
THEN
2286 del=2.0*hipr1(40)*
rlu(0)
2288 p(jl,1)=pt1*
sin(del)
2289 p(jl,2)=-pt1*
cos(del)
2292 p(jl+2,1)=pt3*
sin(del)
2293 p(jl+2,2)=-pt3*
cos(del)
2296 p(jl+1,1)=-
p(jl,1)-
p(jl+2,1)
2297 p(jl+1,2)=-
p(jl,2)-
p(jl+2,2)
2298 p(jl+1,3)=-
p(jl,3)-
p(jl+2,3)
2307 SUBROUTINE atrobo(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
2308 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
2314 IF(imin.LE.0 .OR.
imax.GT.
n .OR. imin.GT.
imax)
RETURN
2316 IF(the**2+
phi**2.GT.1
e-20)
THEN
2332 110
p(i,j)=rot(j,1)*pv(1)+rot(j,2)*pv(2)
2337 IF(bex**2+bey**2+bez**2.GT.1
e-20)
THEN
2342 dga2=1d0-dbex**2-dbey**2-dbez**2
2343 IF(dga2.LE.0d0)
THEN
2352 dbep=dbex*dp(1)+dbey*dp(2)+dbez*dp(3)
2353 dgabep=dga*(dga*dbep/(1d0+dga)+dp(4))
2354 p(i,1)=dp(1)+dgabep*dbex
2355 p(i,2)=dp(2)+dgabep*dbey
2356 p(i,3)=dp(3)+dgabep*dbez
2357 p(i,4)=dga*(dp(4)+dbep)
2378 dimension ip(100,2),ipq(50),ipb(50),it(100,2),itq(50),itb(50)
2381 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2383 common/hijdat/hidat0(10,10),hidat(10)
2385 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
2387 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
2388 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
2389 & pjpm(300,500),ntj(300),kftj(300,500),
2390 & pjtx(300,500),pjty(300,500),pjtz(300,500),
2391 & pjte(300,500),pjtm(300,500)
2393 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
2394 & k2sg(900,100),pxsg(900,100),pysg(900,100),
2395 & pzsg(900,100),pesg(900,100),pmsg(900,100)
2397 common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5)
2402 common/lujets/
n,k(9000,5),
p(9000,5),v(9000,5)
2404 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
2406 common/pysubs/msel,msub(200),kfin(2,-40:40),ckin(200)
2408 common/pypars/mstp(200),parp(200),msti(200),pari(200)
2410 common/pyint1/mint(400),vint(400)
2412 common/pyint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
2414 common/pyint5/ngen(0:200,3),xsec(0:200,3)
2416 common/hipyint/mint4,mint5,atco(200,20),atxs(0:200)
2431 IF(iopjet.EQ.1.AND.(nfp(jp,6).NE.0.OR.nft(jt,6).NE.0))
2433 IF(jp.GT.ihnt2(1) .OR. jt.GT.ihnt2(3))
RETURN
2434 IF(nfp(jp,6).LT.0 .OR. nft(jt,6).LT.0)
RETURN
2438 epp=pp(jp,4)+pp(jp,3)
2439 epm=pp(jp,4)-pp(jp,3)
2440 etp=pt(jt,4)+pt(jt,3)
2441 etm=pt(jt,4)-pt(jt,3)
2442 IF(epp.LT.0.0) go to 1000
2443 IF(epm.LT.0.0) go to 1000
2444 IF(etp.LT.0.0) go to 1000
2445 IF(etm.LT.0.0) go to 1000
2446 IF(epp/(epm+0.01).LE.etp/(etm+0.01))
RETURN
2451 ecut1=hipr1(1)+hipr1(8)+pp(jp,14)+pp(jp,15)
2452 ecut2=hipr1(1)+hipr1(8)+pt(jt,14)+pt(jt,15)
2453 IF(pp(jp,4).LE.ecut1)
THEN
2454 nfp(jp,6)=-abs(nfp(jp,6))
2457 IF(pt(jt,4).LE.ecut2)
THEN
2458 nft(jt,6)=-abs(nft(jt,6))
2467 IF(nfp(jp,10).EQ.0 .AND. nft(jt,10).EQ.0)
THEN
2475 coef(11,i)=atco(11,i)
2476 coef(12,i)=atco(12,i)
2477 coef(28,i)=atco(28,i)
2483 IF(xsec(11,1).NE.0) isub11=1
2484 IF(xsec(12,1).NE.0) isub12=1
2485 IF(xsec(28,1).NE.0) isub28=1
2486 mint(44)=mint4-isub11-isub12-isub28
2487 mint(45)=mint5-isub11-isub12-isub28
2488 xsec(0,1)=atxs(0)-atxs(11)-atxs(12)-atxs(28)
2503 IF(jj.NE.1) go to 155
2505 IF(k(7,2).EQ.-k(8,2))
THEN
2506 qmass2=(
p(7,4)+
p(8,4))**2-(
p(7,1)+
p(8,1))**2
2507 & -(
p(7,2)+
p(8,2))**2-(
p(7,3)+
p(8,3))**2
2509 IF(qmass2.LT.(2.0*qm+hipr1(1))**2) go to 155
2521 IF(pep.LE.ecut1)
THEN
2523 IF(misp.LT.50) go to 155
2524 nfp(jp,6)=-abs(nfp(jp,6))
2527 IF(pet.LE.ecut2)
THEN
2529 IF(mist.LT.50) go to 155
2530 nft(jt,6)=-abs(nft(jt,6))
2538 IF(
wp.LT.0.0 .OR. wm.LT.0.0)
THEN
2540 IF(miss.LT.50) go to 155
2545 ampx=
sqrt((ecut1-hipr1(8))**2+pxp**2+pyp**2+0.01)
2546 amtx=
sqrt((ecut2-hipr1(8))**2+pxt**2+pyt**2+0.01)
2548 IF(sw.LT.sxx.OR.vint(43).LT.hipr1(1))
THEN
2550 IF(miss.LT.50) go to 155
2562 hint1(46)=
sqrt(
p(7,1)**2+
p(7,2)**2)
2568 hint1(56)=
sqrt(
p(8,1)**2+
p(8,2)**2)
2572 pinirad=(1.0-
exp(-2.0*(vint(47)-hidat(1))))
2573 & /(1.0+
exp(-2.0*(vint(47)-hidat(1))))
2575 IF(
rlu(0).LE.pinirad) i_inirad=1
2576 IF(k(7,2).EQ.-k(8,2)) go to 190
2577 IF(k(7,2).EQ.21.AND.k(8,2).EQ.21.AND.iopjet.EQ.1) go to 190
2600 IF(k(i,3).EQ.1 .OR. k(i,3).EQ.2.OR.
2601 & abs(k(i,2)).GT.30) go to 180
2603 IF(k(i,3).EQ.7)
THEN
2604 hint1(47)=hint1(47)+
p(i,1)
2605 hint1(48)=hint1(48)+
p(i,2)
2606 hint1(49)=hint1(49)+
p(i,3)
2607 hint1(50)=hint1(50)+
p(i,4)
2609 IF(k(i,3).EQ.8)
THEN
2610 hint1(67)=hint1(67)+
p(i,1)
2611 hint1(68)=hint1(68)+
p(i,2)
2612 hint1(69)=hint1(69)+
p(i,3)
2613 hint1(70)=hint1(70)+
p(i,4)
2616 IF(k(i,2).GT.21.AND.k(i,2).LE.30)
THEN
2630 IF(k(i,3).EQ.7.OR.k(i,3).EQ.3)
THEN
2631 IF(k(i,3).EQ.7.AND.k(i,2).NE.21.AND.k(i,2).EQ.k(7,2)
2632 & .AND.is7.EQ.0)
THEN
2642 IF(k(i,3).EQ.3.AND.(k(i,2).NE.21.OR.
2643 & i_inirad.EQ.0))
THEN
2653 IF(k(i,2).NE.21)
THEN
2654 IF(k(i,2).GT.0)
THEN
2658 ELSE IF(k(i,2).LT.0)
THEN
2664 ELSE IF(k(i,3).EQ.8.OR.k(i,3).EQ.4)
THEN
2665 IF(k(i,3).EQ.8.AND.k(i,2).NE.21.AND.k(i,2).EQ.k(8,2)
2666 & .AND.is8.EQ.0)
THEN
2676 IF(k(i,3).EQ.4.AND.(k(i,2).NE.21.OR.
2677 & i_inirad.EQ.0))
THEN
2687 IF(k(i,2).NE.21)
THEN
2688 IF(k(i,2).GT.0)
THEN
2692 ELSE IF(k(i,2).LT.0)
THEN
2702 IF(lpq.NE.lpb .OR. ltq.NE.ltb)
THEN
2704 IF(miss.LE.50) go to 155
2705 WRITE(6,*)
' Q -QBAR NOT MATCHED IN HIJHRD'
2714 IF(j.GT.jpp) go to 182
2715 IF(ip(j,2).EQ.0)
THEN
2717 ELSE IF(ip(j,2).NE.0)
THEN
2721 ip(j,1)=ip(ipq(lp),1)
2722 ip(j,2)=ip(ipq(lp),2)
2727 ELSE IF(ip2.LT.0)
THEN
2733 ip(j+1,1)=ip(ipb(lp),1)
2734 ip(j+1,2)=ip(ipb(lp),2)
2739 ELSE IF(ip2.LT.0)
THEN
2749 IF(j.GT.jtt) go to 184
2750 IF(it(j,2).EQ.0)
THEN
2752 ELSE IF(it(j,2).NE.0)
THEN
2756 it(j,1)=it(itq(lt),1)
2757 it(j,2)=it(itq(lt),2)
2762 ELSE IF(it2.LT.0)
THEN
2768 it(j+1,1)=it(itb(lt),1)
2769 it(j+1,2)=it(itb(lt),2)
2774 ELSE IF(it2.LT.0)
THEN
2784 IF(npj(jp)+jpp.GT.mxjt.OR.ntj(jt)+jtt.GT.mxjt)
THEN
2786 WRITE(6,*)
'number of partons per string exceeds'
2787 WRITE(6,*)
'the common block size'
2792 kfpj(jp,npj(jp)+j)=k(ip(j,1),2)
2793 pjpx(jp,npj(jp)+j)=
p(ip(j,1),1)
2794 pjpy(jp,npj(jp)+j)=
p(ip(j,1),2)
2795 pjpz(jp,npj(jp)+j)=
p(ip(j,1),3)
2796 pjpe(jp,npj(jp)+j)=
p(ip(j,1),4)
2797 pjpm(jp,npj(jp)+j)=
p(ip(j,1),5)
2801 kftj(jt,ntj(jt)+j)=k(it(j,1),2)
2802 pjtx(jt,ntj(jt)+j)=
p(it(j,1),1)
2803 pjty(jt,ntj(jt)+j)=
p(it(j,1),2)
2804 pjtz(jt,ntj(jt)+j)=
p(it(j,1),3)
2805 pjte(jt,ntj(jt)+j)=
p(it(j,1),4)
2806 pjtm(jt,ntj(jt)+j)=
p(it(j,1),5)
2814 IF(k(7,2).NE.21.AND.k(8,2).NE.21.AND.
2815 & k(7,2)*k(8,2).GT.0) go to 155
2820 IF(k(i,3).EQ.1.OR.k(i,3).EQ.2.OR.
2821 & abs(k(i,2)).GT.30) go to 200
2822 IF(k(i,2).GT.21.AND.k(i,2).LE.30)
THEN
2836 IF(k(i,3).EQ.3.AND.(k(i,2).NE.21.OR.
2837 & i_inirad.EQ.0))
THEN
2844 IF(k(i,3).EQ.4.AND.(k(i,2).NE.21.OR.
2845 & i_inirad.EQ.0))
THEN
2855 IF(k(i,2).NE.21)
THEN
2856 IF(k(i,2).GT.0)
THEN
2860 ELSE IF(k(i,2).LT.0)
THEN
2869 IF(miss.LE.50) go to 155
2870 WRITE(6,*) lpq,lpb,
'Q-QBAR NOT CONSERVED OR NOT MATCHED'
2879 IF(j.GT.jpp) go to 222
2880 IF(ip(j,2).EQ.0) go to 220
2884 ip(j,1)=ip(ipq(lp),1)
2885 ip(j,2)=ip(ipq(lp),2)
2890 ELSE IF(ip2.LT.0)
THEN
2897 ip(j+1,1)=ip(ipb(lp),1)
2898 ip(j+1,2)=ip(ipb(lp),2)
2903 ELSE IF(ip2.LT.0)
THEN
2916 ip(2*l0-3,1)=ip(ipq(l0),1)
2917 ip(2*l0-3,2)=ip(ipq(l0),2)
2922 ELSE IF(ip2.LT.0)
THEN
2929 ip(2*l0-2,1)=ip(ipb(l0),1)
2930 ip(2*l0-2,2)=ip(ipb(l0),2)
2935 ELSE IF(ip2.LT.0)
THEN
2944 ip(2*lpq-1,1)=ip(ipq(1),1)
2945 ip(2*lpq-1,2)=ip(ipq(1),2)
2950 ELSE IF(ip2.LT.0)
THEN
2958 ip(jpp,1)=ip(ipb(1),1)
2959 ip(jpp,2)=ip(ipb(1),2)
2964 ELSE IF(ip2.LT.0)
THEN
2971 IF(nsg.GE.mxsg)
THEN
2973 WRITE(6,*)
'number of jets forming single strings exceeds'
2974 WRITE(6,*)
'the common block size'
2977 IF(jpp.GT.mxsj)
THEN
2979 WRITE(6,*)
'number of partons per single jet system'
2980 WRITE(6,*)
'exceeds the common block size'
2991 k2sg(nsg,i)=k(ip(i,1),2)
2992 IF(k2sg(nsg,i).LT.0) k1sg(nsg,i)=1
2993 pxsg(nsg,i)=
p(ip(i,1),1)
2994 pysg(nsg,i)=
p(ip(i,1),2)
2995 pzsg(nsg,i)=
p(ip(i,1),3)
2996 pesg(nsg,i)=
p(ip(i,1),4)
2997 pmsg(nsg,i)=
p(ip(i,1),5)
3013 nfp(jp,6)=nfp(jp,6)+1
3014 nft(jt,6)=nft(jt,6)+1
3018 IF(ihpr2(10).EQ.0)
RETURN
3019 WRITE(6,*)
'Fatal HIJHRD error'
3020 WRITE(6,*) jp,
' proj E+,E-',epp,epm,
' status',nfp(jp,5)
3021 WRITE(6,*) jt,
' targ E+,E_',etp,etm,
' status',nft(jt,5)
3044 CHARACTER beam*16,targ*16
3045 dimension xsec0(8,0:200),coef0(8,200,20),ini(8),
3046 & mint44(8),mint45(8)
3049 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3051 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3053 common/hipyint/mint4,mint5,atco(200,20),atxs(0:200)
3056 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
3058 common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
3060 common/pysubs/msel,msub(200),kfin(2,-40:40),ckin(200)
3062 common/pypars/mstp(200),parp(200),msti(200),pari(200)
3064 common/pyint1/mint(400),vint(400)
3066 common/pyint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
3068 common/pyint5/ngen(0:200,3),xsec(0:200,3)
3070 DATA ini/8*0/i_last/-1/
3076 IF(ihnt2(5).NE.0 .AND. ihnt2(6).NE.0)
THEN
3078 ELSE IF(ihnt2(5).NE.0 .AND. ihnt2(6).EQ.0)
THEN
3080 IF(nft(jt,4).EQ.2112) i_type=2
3081 ELSE IF(ihnt2(5).EQ.0 .AND. ihnt2(6).NE.0)
THEN
3083 IF(nfp(jp,4).EQ.2112) i_type=2
3085 IF(nfp(jp,4).EQ.2212 .AND. nft(jt,4).EQ.2212)
THEN
3087 ELSE IF(nfp(jp,4).EQ.2212 .AND. nft(jt,4).EQ.2112)
THEN
3089 ELSE IF(nfp(jp,4).EQ.2112 .AND. nft(jt,4).EQ.2212)
THEN
3096 IF(i_trig.NE.0) go to 160
3097 IF(i_trig.EQ.i_last) go to 150
3109 IF(ihpr2(2).EQ.0.OR.ihpr2(2).EQ.2) mstp(61)=0
3110 IF(ihpr2(2).EQ.0.OR.ihpr2(2).EQ.1) mstp(71)=0
3118 IF(ihpr2(10).EQ.0) mstp(122)=0
3124 IF(hipr1(9).LE.hipr1(8)) ckin(4)=-1.0
3139 DO 110 j=1,min(8,mdcy(21,3))
3140 110 mdme(mdcy(21,2)+j-1,1)=0
3142 IF(hint1(1).GE.20.0 .and. ihpr2(18).EQ.1) isel=5
3143 mdme(mdcy(21,2)+isel-1,1)=1
3149 150
IF(ini(i_type).NE.0) go to 800
3155 IF(i_trig.EQ.i_last) go to 260
3156 parp(81)=abs(hipr1(10))-0.25
3157 ckin(5)=abs(hipr1(10))-0.25
3158 ckin(3)=abs(hipr1(10))-0.25
3159 ckin(4)=abs(hipr1(10))+0.25
3160 IF(hipr1(10).LT.hipr1(8)) ckin(4)=-1.0
3166 IF(ihpr2(3).EQ.1)
THEN
3178 DO 102 j=1,min(8,mdcy(21,3))
3179 102 mdme(mdcy(21,2)+j-1,1)=0
3181 IF(hint1(1).GE.20.0 .and. ihpr2(18).EQ.1) isel=5
3182 mdme(mdcy(21,2)+isel-1,1)=1
3184 ELSE IF(ihpr2(3).EQ.2)
THEN
3190 ELSE IF(ihpr2(3).EQ.3)
THEN
3191 ckin(3)=max(0.0,hipr1(10))
3196 DO 105 j=1,min(8,mdcy(21,3))
3197 105 mdme(mdcy(21,2)+j-1,1)=0
3199 IF(hint1(1).GE.20.0 .and. ihpr2(18).EQ.1) isel=5
3200 mdme(mdcy(21,2)+isel-1,1)=1
3203 260
IF(ini(i_type).NE.0) go to 800
3207 IF(ihpr2(10).EQ.0) mstp(122)=0
3208 IF(nfp(jp,4).EQ.2212)
THEN
3210 ELSE IF(nfp(jp,4).EQ.-2212)
THEN
3212 ELSE IF(nfp(jp,4).EQ.2112)
THEN
3214 ELSE IF(nfp(jp,4).EQ.-2112)
THEN
3216 ELSE IF(nfp(jp,4).EQ.211)
THEN
3218 ELSE IF(nfp(jp,4).EQ.-211)
THEN
3220 ELSE IF(nfp(jp,4).EQ.321)
THEN
3222 ELSE IF(nfp(jp,4).EQ.-321)
THEN
3225 WRITE(6,*)
'unavailable beam type', nfp(jp,4)
3227 IF(nft(jt,4).EQ.2212)
THEN
3229 ELSE IF(nft(jt,4).EQ.-2212)
THEN
3231 ELSE IF(nft(jt,4).EQ.2112)
THEN
3233 ELSE IF(nft(jt,4).EQ.-2112)
THEN
3235 ELSE IF(nft(jt,4).EQ.211)
THEN
3237 ELSE IF(nft(jt,4).EQ.-211)
THEN
3239 ELSE IF(nft(jt,4).EQ.321)
THEN
3241 ELSE IF(nft(jt,4).EQ.-321)
THEN
3244 WRITE(6,*)
'unavailable target type', nft(jt,4)
3251 CALL
pyinit(
'CMS',beam,targ,hint1(1))
3254 mint44(i_type)=mint(44)
3255 mint45(i_type)=mint(45)
3257 xsec0(i_type,0)=xsec(0,1)
3260 xsec0(i_type,i)=xsec(i,1)
3263 coef0(i_type,i,j)=coef(i,j)
3273 800 mint(44)=mint44(i_type)
3274 mint(45)=mint45(i_type)
3277 xsec(0,1)=xsec0(i_type,0)
3280 xsec(i,1)=xsec0(i_type,i)
3283 coef(i,j)=coef0(i_type,i,j)
3295 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3297 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3299 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
3300 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
3301 & pjpm(300,500),ntj(300),kftj(300,500),
3302 & pjtx(300,500),pjty(300,500),pjtz(300,500),
3303 & pjte(300,500),pjtm(300,500)
3305 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
3306 & k2sg(900,100),pxsg(900,100),pysg(900,100),
3307 & pzsg(900,100),pesg(900,100),pmsg(900,100)
3309 common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5)
3320 IF(ihnt2(5).NE.0) ipp=ihnt2(5)
3321 IF(ihnt2(6).NE.0) ipt=ihnt2(6)
3327 pp(i,3)=
sqrt(hint1(1)**2/4.0-hint1(8)**2)
3345 IF(i.GT.abs(ihnt2(2))) nfp(i,3)=2112
3346 CALL
attflv(nfp(i,3),idq,idqq)
3350 IF(abs(idq).GT.1000.OR.(abs(idq*idqq).LT.100.AND.
3351 &
rlu(0).LT.0.5)) nfp(i,15)=1
3359 pt(i,3)=-
sqrt(hint1(1)**2/4.0-hint1(9)**2)
3360 pt(i,4)=hint1(1)/2.0
3377 IF(i.GT.abs(ihnt2(4))) nft(i,3)=2112
3378 CALL
attflv(nft(i,3),idq,idqq)
3382 IF(abs(idq).GT.1000.OR.(abs(idq*idqq).LT.100.AND.
3383 &
rlu(0).LT.0.5)) nft(i,15)=-1
3396 IF(abs(id).LT.100)
THEN
3400 IF(abs(idq).EQ.3) nsign=-1
3421 IF(abs(id).EQ.2112) idq=1
3424 IF(x.LE.0.5) go to 30
3425 IF(x.GT.0.666667) go to 10
3430 IF(abs(id).EQ.2112)
THEN
3447 dimension psc1(5),psc2(5)
3450 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3454 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3456 IF(jp.EQ.0 .OR. jt.EQ.0) go to 25
3462 dpp1=psc1(1)-pp(jp,1)
3463 dpp2=psc1(2)-pp(jp,2)
3464 dpt1=psc2(1)-pt(jt,1)
3465 dpt2=psc2(2)-pt(jt,2)
3466 pp(jp,6)=pp(jp,6)+dpp1/2.0
3467 pp(jp,7)=pp(jp,7)+dpp2/2.0
3468 pp(jp,8)=pp(jp,8)+dpp1/2.0
3469 pp(jp,9)=pp(jp,9)+dpp2/2.0
3470 pt(jt,6)=pt(jt,6)+dpt1/2.0
3471 pt(jt,7)=pt(jt,7)+dpt2/2.0
3472 pt(jt,8)=pt(jt,8)+dpt1/2.0
3473 pt(jt,9)=pt(jt,9)+dpt2/2.0
3478 nfp(jp,5)=max(1,nfp(jp,5))
3479 nft(jt,5)=max(1,nft(jt,5))
3484 25
IF(jp.EQ.0) go to 45
3485 pabs=
sqrt(pp(jp,1)**2+pp(jp,2)**2+pp(jp,3)**2)
3490 IF(i.EQ.jp) go to 40
3495 IF(dis.LE.0) go to 40
3497 r2=
bb*hipr1(40)/hipr1(31)/0.1
3499 gs=1.0-
exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3501 gs0=1.0-
exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3503 IF(
rlu(0).GT.gs/gs0) go to 40
3509 dpp1=psc1(1)-pp(jp,1)
3510 dpp2=psc1(2)-pp(jp,2)
3511 dpt1=psc2(1)-pp(i,1)
3512 dpt2=psc2(2)-pp(i,2)
3513 pp(jp,6)=pp(jp,6)+dpp1/2.0
3514 pp(jp,7)=pp(jp,7)+dpp2/2.0
3515 pp(jp,8)=pp(jp,8)+dpp1/2.0
3516 pp(jp,9)=pp(jp,9)+dpp2/2.0
3517 pp(i,6)=pp(i,6)+dpt1/2.0
3518 pp(i,7)=pp(i,7)+dpt2/2.0
3519 pp(i,8)=pp(i,8)+dpt1/2.0
3520 pp(i,9)=pp(i,9)+dpt2/2.0
3525 nfp(i,5)=max(1,nfp(i,5))
3528 45
IF(jt.EQ.0) go to 80
3529 50 pabs=
sqrt(pt(jt,1)**2+pt(jt,2)**2+pt(jt,3)**2)
3534 IF(i.EQ.jt) go to 70
3539 IF(dis.LE.0) go to 70
3541 r2=
bb*hipr1(40)/hipr1(31)/0.1
3543 gs=(1.0-
exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3545 gs0=(1.0-
exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3547 IF(
rlu(0).GT.gs/gs0) go to 70
3553 dpp1=psc1(1)-pt(jt,1)
3554 dpp2=psc1(2)-pt(jt,2)
3555 dpt1=psc2(1)-pt(i,1)
3556 dpt2=psc2(2)-pt(i,2)
3557 pt(jt,6)=pt(jt,6)+dpp1/2.0
3558 pt(jt,7)=pt(jt,7)+dpp2/2.0
3559 pt(jt,8)=pt(jt,8)+dpp1/2.0
3560 pt(jt,9)=pt(jt,9)+dpp2/2.0
3561 pt(i,6)=pt(i,6)+dpt1/2.0
3562 pt(i,7)=pt(i,7)+dpt2/2.0
3563 pt(i,8)=pt(i,8)+dpt1/2.0
3564 pt(i,9)=pt(i,9)+dpt2/2.0
3569 nft(i,5)=max(1,nft(i,5))
3581 IMPLICIT DOUBLE PRECISION(
d)
3582 dimension psc1(5),psc2(5)
3583 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3588 cc=1.0-hint1(12)/hint1(13)
3589 rr=(1.0-cc)*hint1(13)/hint1(12)/(1.0-hipr1(33))-1.0
3590 bb=0.5*(3.0+rr+
sqrt(9.0+10.0*rr+rr**2))
3591 ep=
sqrt((psc1(1)-psc2(1))**2+(psc1(2)-psc2(2))**2
3592 & +(psc1(3)-psc2(3))**2)
3593 IF(ep.LE.0.1)
RETURN
3594 els0=98.0/ep+52.0*(1.0+rr)**2
3595 pcm1=psc1(1)+psc2(1)
3596 pcm2=psc1(2)+psc2(2)
3597 pcm3=psc1(3)+psc2(3)
3601 amm=ecm**2-pcm1**2-pcm2**2-pcm3**2
3602 IF(amm.LE.psc1(5)+psc2(5))
RETURN
3605 pmax=(amm**2+am1**2+am2**2-2.0*amm*am1-2.0*amm*am2
3606 & -2.0*am1*am2)/4.0/amm
3608 20
tt=
rlu(0)*min(pmax,1.5)
3609 els=98.0*
exp(-2.8*
tt)/ep
3611 IF(
rlu(0).GT.els/els0) go to 20
3617 db=
sqrt(dbx**2+dby**2+dbz**2)
3618 IF(db.GT.0.99999999d0)
THEN
3619 dbx=dbx*(0.99999999d0/db)
3620 dby=dby*(0.99999999d0/db)
3621 dbz=dbz*(0.99999999d0/db)
3623 WRITE(6,*)
' (HIJELS) boost vector too large'
3626 dga=1d0/
sqrt(1d0-db**2)
3632 dbp=dbx*dp1+dby*dp2+dbz*dp3
3633 dgabp=dga*(dga*dbp/(1d0+dga)+dp4)
3634 psc1(1)=dp1+dgabp*dbx
3635 psc1(2)=dp2+dgabp*dby
3636 psc1(3)=dp3+dgabp*dbz
3637 psc1(4)=dga*(dp4+dbp)
3643 dbp=dbx*dp1+dby*dp2+dbz*dp3
3644 dgabp=dga*(dga*dbp/(1d0+dga)+dp4)
3645 psc2(1)=dp1+dgabp*dbx
3646 psc2(2)=dp2+dgabp*dby
3647 psc2(3)=dp3+dgabp*dbz
3648 psc2(4)=dga*(dp4+dbp)
3662 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3664 common/hijdat/hidat0(10,10),hidat(10)
3668 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
3669 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
3670 & pjpm(300,500),ntj(300),kftj(300,500),
3671 & pjtx(300,500),pjty(300,500),pjtz(300,500),
3672 & pjte(300,500),pjtm(300,500)
3674 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
3675 & k2sg(900,100),pxsg(900,100),pysg(900,100),
3676 & pzsg(900,100),pesg(900,100),pmsg(900,100)
3678 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3680 common/dpmcom1/jjp,jjt,amp,amt,apx0,atx0,ampn,amtn,amp0,amt0,
3681 & nfdp,nfdt,
wp,wm,sw,xremp,xremt,dpkc1,dpkc2,pp11,pp12,
3682 & pt11,pt12,ptp2,ptt2
3684 common/dpmcom2/ndpm,kdpm(20,2),pdpm1(20,5),pdpm2(20,5)
3697 IF(jp.GT.ihnt2(1) .OR. jt.GT.ihnt2(3))
RETURN
3699 epp=pp(jp,4)+pp(jp,3)
3700 epm=pp(jp,4)-pp(jp,3)
3701 etp=pt(jt,4)+pt(jt,3)
3702 etm=pt(jt,4)-pt(jt,3)
3709 IF(
wp.LT.0.0 .OR. wm.LT.0.0) go to 1000
3712 IF(epp.LT.0.0) go to 1000
3713 IF(epm.LT.0.0) go to 1000
3714 IF(etp.LT.0.0) go to 1000
3715 IF(etm.LT.0.0) go to 1000
3716 IF(epp/(epm+0.01).LE.etp/(etm+0.01))
RETURN
3739 IF(nfp(jp,10).EQ.1.OR.nft(jt,10).EQ.1)
THEN
3740 IF(nfp(jp,10).EQ.1)
THEN
3741 phi1=
ulangl(pp(jp,10),pp(jp,11))
3742 ppjet=
sqrt(pp(jp,10)**2+pp(jp,11)**2)
3747 IF(nft(jt,10).EQ.1)
THEN
3748 phi2=
ulangl(pt(jt,10),pt(jt,11))
3749 ptjet=
sqrt(pt(jt,10)**2+pt(jt,11)**2)
3754 IF(ihpr2(4).GT.0.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1)
THEN
3755 IF(nfp(jp,10).EQ.0)
THEN
3757 ELSE IF(nft(jt,10).EQ.0)
THEN
3760 phi=(phi1+phi2-hipr1(40))/2.0
3762 bx=hint1(19)*
cos(hint1(20))
3763 by=hint1(19)*
sin(hint1(20))
3768 r1=max(1.2*ihnt2(1)**0.3333333,
3769 &
sqrt(xp0**2+yp0**2))
3770 r2=max(1.2*ihnt2(3)**0.3333333,
3771 &
sqrt((xt0-bx)**2+(yt0-by)**2))
3772 IF(abs(
cos(
phi)).LT.1.0
e-5)
THEN
3775 dd3=abs(by+
sqrt(r2**2-(xp0-bx)**2)-yp0)
3776 dd4=abs(by-
sqrt(r2**2-(xp0-bx)**2)-yp0)
3783 IF(dd.LT.0.0) go to 10
3786 dd1=abs((xx1-xp0)/
cos(
phi))
3787 dd2=abs((xx2-xp0)/
cos(
phi))
3794 IF(dd.LT.0.0) go to 10
3797 dd3=abs((xx1-xt0)/
cos(
phi))
3798 dd4=abs((xx2-xt0)/
cos(
phi))
3802 IF(dd1.LT.hipr1(13)) dd1=0.0
3803 IF(dd2.LT.hipr1(13)) dd2=0.0
3804 IF(nfp(jp,10).EQ.1.AND.ppjet.GT.hipr1(11))
THEN
3805 dp1=dd1*hipr1(14)/2.0
3806 dp1=min(dp1,ppjet-hipr1(11))
3810 pkc11=pp(jp,10)-dpx1
3811 pkc12=pp(jp,11)-dpy1
3813 cthep=pp(jp,12)/
sqrt(pp(jp,12)**2+ppjet**2)
3814 dpz1=dp1*cthep/
sqrt(1.0-cthep**2)
3815 dpe1=
sqrt(dpx1**2+dpy1**2+dpz1**2)
3816 eppprm=pp(jp,4)+pp(jp,3)-dpe1-dpz1
3817 epmprm=pp(jp,4)-pp(jp,3)-dpe1+dpz1
3818 IF(eppprm.LE.0.0.OR.epmprm.LE.0.0) go to 15
3825 pjpx(jp,npj(jp))=dpx1
3826 pjpy(jp,npj(jp))=dpy1
3827 pjpz(jp,npj(jp))=dpz1
3828 pjpe(jp,npj(jp))=dpe1
3829 pjpm(jp,npj(jp))=0.0
3830 pp(jp,3)=pp(jp,3)-dpz1
3831 pp(jp,4)=pp(jp,4)-dpe1
3834 15
IF(nft(jt,10).EQ.1.AND.ptjet.GT.hipr1(11))
THEN
3835 dp2=dd2*hipr1(14)/2.0
3836 dp2=min(dp2,ptjet-hipr1(11))
3840 pkc21=pt(jt,10)-dpx2
3841 pkc22=pt(jt,11)-dpy2
3843 cthet=pt(jt,12)/
sqrt(pt(jt,12)**2+ptjet**2)
3844 dpz2=dp2*cthet/
sqrt(1.0-cthet**2)
3845 dpe2=
sqrt(dpx2**2+dpy2**2+dpz2**2)
3846 etpprm=pt(jt,4)+pt(jt,3)-dpe2-dpz2
3847 etmprm=pt(jt,4)-pt(jt,3)-dpe2+dpz2
3848 IF(etpprm.LE.0.0.OR.etmprm.LE.0.0) go to 16
3855 pjtx(jt,ntj(jt))=dpx2
3856 pjty(jt,ntj(jt))=dpy2
3857 pjtz(jt,ntj(jt))=dpz2
3858 pjte(jt,ntj(jt))=dpe2
3859 pjtm(jt,ntj(jt))=0.0
3860 pt(jt,3)=pt(jt,3)-dpz2
3861 pt(jt,4)=pt(jt,4)-dpe2
3864 16 dpkc11=-(pp(jp,10)-pkc11)/2.0
3865 dpkc12=-(pp(jp,11)-pkc12)/2.0
3866 dpkc21=-(pt(jt,10)-pkc21)/2.0
3867 dpkc22=-(pt(jt,11)-pkc22)/2.0
3877 10 ptp02=pp(jp,1)**2+pp(jp,2)**2
3878 ptt02=pt(jt,1)**2+pt(jt,2)**2
3880 amq=max(pp(jp,14)+pp(jp,15),pt(jt,14)+pt(jt,15))
3887 IF(nfp(jp,5).LE.2.AND.nfp(jp,3).NE.0)
THEN
3889 nfdp=nfp(jp,3)+2*nfp(jp,3)/abs(nfp(jp,3))
3891 IF(dpm0.LE.0.0)
THEN
3892 nfdp=nfdp-2*nfdp/abs(nfdp)
3899 IF(nft(jt,5).LE.2.AND.nft(jt,3).NE.0)
THEN
3901 nfdt=nft(jt,3)+2*nft(jt,3)/abs(nft(jt,3))
3903 IF(dtm0.LE.0.0)
THEN
3904 nfdt=nfdt-2*nfdt/abs(nfdt)
3909 ampn=
sqrt(amp0**2+ptp02)
3910 amtn=
sqrt(amt0**2+ptt02)
3911 snn=(ampn+amtn)**2+0.001
3913 IF(sw.LT.snn+0.001) go to 4000
3916 20 swptn=4.0*(max(amp0,amt0)**2+max(ptp02,ptt02))
3917 swptd=4.0*(max(dpm0,dtm0)**2+max(ptp02,ptt02))
3918 swptx=4.0*(amx**2+max(ptp02,ptt02))
3919 IF(sw.LE.swptn)
THEN
3921 ELSE IF(sw.GT.swptn .AND. sw.LE.swptd
3922 & .AND.npj(jp).EQ.0.AND.ntj(jt).EQ.0)
THEN
3923 pkcmx=
sqrt(sw/4.0-max(amp0,amt0)**2)
3924 & -
sqrt(max(ptp02,ptt02))
3925 ELSE IF(sw.GT.swptd .AND. sw.LE.swptx
3926 & .AND.npj(jp).EQ.0.AND.ntj(jt).EQ.0)
THEN
3927 pkcmx=
sqrt(sw/4.0-max(dpm0,dtm0)**2)
3928 & -
sqrt(max(ptp02,ptt02))
3929 ELSE IF(sw.GT.swptx)
THEN
3930 pkcmx=
sqrt(sw/4.0-amx**2)-
sqrt(max(ptp02,ptt02))
3935 IF(nfp(jp,10).EQ.1.OR.nft(jt,10).EQ.1)
THEN
3936 IF(pkc1.GT.pkcmx)
THEN
3938 pkc11=pkc1*
cos(phi1)
3939 pkc12=pkc1*
sin(phi1)
3940 dpkc11=-(pp(jp,10)-pkc11)/2.0
3941 dpkc12=-(pp(jp,11)-pkc12)/2.0
3943 IF(pkc2.GT.pkcmx)
THEN
3945 pkc21=pkc2*
cos(phi2)
3946 pkc22=pkc2*
sin(phi2)
3947 dpkc21=-(pt(jt,10)-pkc21)/2.0
3948 dpkc22=-(pt(jt,11)-pkc22)/2.0
3952 nfp(jp,10)=-nfp(jp,10)
3953 nft(jt,10)=-nft(jt,10)
3959 IF(ihpr2(13).NE.0 .AND.
rlu(0).LE.hidat(4)) i_sng=1
3960 IF((nfp(jp,5).EQ.3 .OR.nft(jt,5).EQ.3).OR.
3961 & (npj(jp).NE.0.OR.nfp(jp,10).NE.0).OR.
3962 & (ntj(jt).NE.0.OR.nft(jt,10).NE.0)) i_sng=0
3965 IF(ihpr2(5).EQ.0)
THEN
3966 pkc=hipr1(2)*
sqrt(-alog(1.0-
rlu(0)
3967 & *(1.0-
exp(-pkcmx**2/hipr1(2)**2))))
3970 pkc=
hirnd2(3,0.0,pkcmx**2)
3972 IF(pkc.GT.hipr1(20))
3973 & pkc=hipr1(2)*
sqrt(-alog(
exp(-hipr1(20)**2/hipr1(2)**2)
3974 & -
rlu(0)*(
exp(-hipr1(20)**2/hipr1(2)**2)-
3975 &
exp(-pkcmx**2/hipr1(2)**2))))
3977 IF(i_sng.EQ.1) pkc=0.65*
sqrt(
3978 & -alog(1.0-
rlu(0)*(1.0-
exp(-pkcmx**2/0.65**2))))
3980 30 phi0=2.0*hipr1(40)*
rlu(0)
3987 40 pp11=pp(jp,1)+pkc11-dpkc1
3988 pp12=pp(jp,2)+pkc12-dpkc2
3989 pt11=pt(jt,1)+pkc21-dpkc1
3990 pt12=pt(jt,2)+pkc22-dpkc2
3991 ptp2=pp11**2+pp12**2
3992 ptt2=pt11**2+pt12**2
3994 ampn=
sqrt(amp0**2+ptp2)
3995 amtn=
sqrt(amt0**2+ptt2)
3996 snn=(ampn+amtn)**2+0.001
4004 IF(miss.LE.100)
then
4009 &
WRITE(6,*)
'Error occured in Pt kick section of HIJSFT'
4013 ampd=
sqrt(dpm0**2+ptp2)
4014 amtd=
sqrt(dtm0**2+ptt2)
4016 ampx=
sqrt(amx**2+ptp2)
4017 amtx=
sqrt(amx**2+ptt2)
4026 spntd=(ampn+amtd)**2
4027 spntx=(ampn+amtx)**2
4029 spdtn=(ampd+amtn)**2
4030 spxtn=(ampx+amtn)**2
4032 spdtx=(ampd+amtx)**2
4033 spxtd=(ampx+amtd)**2
4046 IF(sw.GT.sxx+0.001)
THEN
4056 IF((nfp(jp,5).EQ.3 .AND.nft(jt,5).EQ.3).OR.
4057 & (npj(jp).NE.0.OR.nfp(jp,10).NE.0).OR.
4058 & (ntj(jt).NE.0.OR.nft(jt,10).NE.0))
THEN
4067 IF(
rlu(0).GT.0.5.OR.(nft(jt,5).GT.2.OR.
4068 & ntj(jt).NE.0.OR.nft(jt,10).NE.0))
THEN
4083 ELSE IF(sw.GT.max(spdtx,spxtd)+0.001 .AND.
4084 & sw.LE.sxx+0.001)
THEN
4085 IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0.AND.
4086 &
rlu(0).GT.0.5).OR.(npj(jp).EQ.0
4087 & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2)
THEN
4093 ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2)
THEN
4101 ELSE IF(sw.GT.min(spdtx,spxtd)+0.001.AND.
4102 & sw.LE.max(spdtx,spxtd)+0.001)
THEN
4103 IF(spdtx.LE.spxtd.AND.npj(jp).EQ.0
4104 & .AND.nfp(jp,5).LE.2)
THEN
4110 ELSE IF(spdtx.GT.spxtd.AND.ntj(jt).EQ.0
4111 & .AND.nft(jt,5).LE.2)
THEN
4120 IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0
4121 & .AND.
rlu(0).GT.0.5).OR.(npj(jp).EQ.0
4122 & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2)
THEN
4128 ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2)
THEN
4136 ELSE IF(sw.GT.max(spntx,spxtn)+0.001 .AND.
4137 & sw.LE.min(spdtx,spxtd)+0.001)
THEN
4138 IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0
4139 & .AND.
rlu(0).GT.0.5).OR.(npj(jp).EQ.0
4140 & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2)
THEN
4146 ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2)
THEN
4154 ELSE IF(sw.GT.min(spntx,spxtn)+0.001 .AND.
4155 & sw.LE.max(spntx,spxtn)+0.001)
THEN
4156 IF(spntx.LE.spxtn.AND.npj(jp).EQ.0
4157 & .AND.nfp(jp,5).LE.2)
THEN
4163 ELSEIF(spntx.GT.spxtn.AND.ntj(jt).EQ.0
4164 & .AND.nft(jt,5).LE.2)
THEN
4172 ELSE IF(sw.LE.min(spntx,spxtn)+0.001 .AND.
4173 & (npj(jp).NE.0 .OR.ntj(jt).NE.0))
THEN
4175 ELSE IF(sw.LE.min(spntx,spxtn)+0.001 .AND.
4176 & nfp(jp,5).GT.2.AND.nft(jt,5).GT.2)
THEN
4178 ELSE IF(sw.GT.sdd+0.001.AND.sw.LE.
4179 & min(spntx,spxtn)+0.001)
THEN
4185 ELSE IF(sw.GT.max(spntd,spdtn)+0.001
4186 & .AND. sw.LE.sdd+0.001)
THEN
4187 IF(
rlu(0).GT.0.5)
THEN
4200 ELSE IF(sw.GT.min(spntd,spdtn)+0.001
4201 & .AND. sw.LE.max(spntd,spdtn)+0.001)
THEN
4202 IF(spntd.GT.spdtn)
THEN
4215 ELSE IF(sw.LE.min(spntd,spdtn)+0.001)
THEN
4222 WRITE(6,*)
' Error in HIJSFT: There is no path to here'
4229 100 nfp5=max(2,nfp(jp,5))
4230 nft5=max(2,nft(jt,5))
4233 IF(bb1**2.LT.4.0*
d1 .OR. bb2**2.LT.4.0*
d2)
THEN
4235 IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4239 IF(
rlu(0).LT.0.5)
THEN
4240 x1=(bb1-
sqrt(bb1**2-4.0*
d1))/2.0
4241 x2=(bb2-
sqrt(bb2**2-4.0*
d2))/2.0
4243 x1=(bb1+
sqrt(bb1**2-4.0*
d1))/2.0
4244 x2=(bb2+
sqrt(bb2**2-4.0*
d2))/2.0
4252 220 nfp5=max(2,nfp(jp,5))
4254 IF(nfp3.EQ.0) nfp5=3
4256 IF(bb2**2.LT.4.0*
d2)
THEN
4258 IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4267 IF(x2*(1.0-x1).LT.(
d2+1.
e-4/sw))
THEN
4269 IF(miss4.LE.1000) go to 222
4276 nft5=max(2,nft(jt,5))
4277 IF(nft3.EQ.0) nft5=3
4279 IF(bb1**2.LT.4.0*
d1)
THEN
4281 IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4290 IF(x1*(1.0-x2).LT.(
d1+1.
e-4/sw))
THEN
4292 IF(miss4.LE.1000) go to 242
4307 IF(bb1**2.LT.4.0*
d1 .OR. bb2**2.LT.4.0*
d2)
THEN
4309 IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4313 xmin1=(bb1-
sqrt(bb1**2-4.0*
d1))/2.0
4314 xmax1=(bb1+
sqrt(bb1**2-4.0*
d1))/2.0
4315 xmin2=(bb2-
sqrt(bb2**2-4.0*
d2))/2.0
4316 xmax2=(bb2+
sqrt(bb2**2-4.0*
d2))/2.0
4318 410 x1=
hirnd2(4,xmin1,xmax1)
4320 IF(nfp(jp,5).EQ.3.OR.nft(jt,5).EQ.3)
THEN
4325 IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000.OR.
4326 & abs(nfp(jp,1)*nfp(jp,2)).LT.100)
THEN
4329 IF(abs(nft(jt,1)*nft(jt,2)).GT.1000000.OR.
4330 & abs(nft(jt,1)*nft(jt,2)).LT.100)
THEN
4339 IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000) x1=1.0-x1
4342 IF(xxp.LT.(
d1+1.
e-4/sw) .OR. xxt.LT.(
d2+1.
e-4/sw))
THEN
4344 IF(miss4.LE.1000) go to 410
4351 IF(x1*(1.0-x2).LT.(ampn**2-1.
e-4)/sw.OR.
4352 & x2*(1.0-x1).LT.(amtn**2-1.
e-4)/sw)
THEN
4354 IF(miss.GT.100.OR.pkc.EQ.0.0) go to 2000
4363 pp(jp,3)=(epp-epm)/2.0
4364 pp(jp,4)=(epp+epm)/2.0
4365 IF(epp*epm-ptp2.LT.0.0) go to 6000
4366 pp(jp,5)=
sqrt(epp*epm-ptp2)
4370 pt(jt,3)=(etp-etm)/2.0
4371 pt(jt,4)=(etp+etm)/2.0
4372 IF(etp*etm-ptt2.LT.0.0) go to 6000
4373 pt(jt,5)=
sqrt(etp*etm-ptt2)
4383 IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000.OR.
4384 & abs(nfp(jp,1)*nfp(jp,2)).LT.100)
THEN
4387 IF(abs(nft(jt,1)*nft(jt,2)).GT.1000000.OR.
4388 & abs(nft(jt,1)*nft(jt,2)).LT.100)
THEN
4391 IF((kickdip.EQ.0.AND.
rlu(0).LT.0.5)
4392 & .OR.(kickdip.NE.0.AND.
rlu(0)
4393 & .LT.0.5/(1.0+(pkc11**2+pkc12**2)/hipr1(22)**2)))
THEN
4394 pp(jp,6)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0+pp(jp,6)
4395 pp(jp,7)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0+pp(jp,7)
4396 pp(jp,8)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0
4398 pp(jp,9)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0
4401 pp(jp,8)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0+pp(jp,8)
4402 pp(jp,9)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0+pp(jp,9)
4403 pp(jp,6)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0
4405 pp(jp,7)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0
4408 pp(jp,1)=pp(jp,6)+pp(jp,8)
4409 pp(jp,2)=pp(jp,7)+pp(jp,9)
4413 IF((kickdit.EQ.0.AND.
rlu(0).LT.0.5)
4414 & .OR.(kickdit.NE.0.AND.
rlu(0)
4415 & .LT.0.5/(1.0+(pkc21**2+pkc22**2)/hipr1(22)**2)))
THEN
4416 pt(jt,6)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0+pt(jt,6)
4417 pt(jt,7)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0+pt(jt,7)
4418 pt(jt,8)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0
4420 pt(jt,9)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0
4423 pt(jt,8)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0+pt(jt,8)
4424 pt(jt,9)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0+pt(jt,9)
4425 pt(jt,6)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0
4427 pt(jt,7)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0
4430 pt(jt,1)=pt(jt,6)+pt(jt,8)
4431 pt(jt,2)=pt(jt,7)+pt(jt,9)
4434 IF(npj(jp).NE.0) nfp(jp,5)=3
4435 IF(ntj(jt).NE.0) nft(jt,5)=3
4437 IF(epp/(epm+0.0001).LT.etp/(etm+0.0001).AND.
4438 & abs(nfp(jp,1)*nfp(jp,2)).LT.1000000)
THEN
4441 pp(jp,jsb)=pt(jt,jsb)
4444 nfp(jp,jsb)=nft(jt,jsb)
4455 IF(ihpr2(10).EQ.0)
RETURN
4456 WRITE(6,*)
' Fatal HIJSFT start error,abandon this event'
4457 WRITE(6,*)
' PROJ E+,E-,W+',epp,epm,
wp
4458 WRITE(6,*)
' TARG E+,E-,W-',etp,etm,wm
4459 WRITE(6,*)
' W+*W-, (APN+ATN)^2',sw,snn
4462 IF(ihpr2(10).EQ.0)
RETURN
4463 WRITE(6,*)
' (2)energy partition fail,'
4464 WRITE(6,*)
' HIJSFT not performed, but continue'
4465 WRITE(6,*)
' MP1,MPN',x1*(1.0-x2)*sw,ampn**2
4466 WRITE(6,*)
' MT2,MTN',x2*(1.0-x1)*sw,amtn**2
4469 IF(ihpr2(10).EQ.0)
RETURN
4470 WRITE(6,*)
' (3)something is wrong with the pt kick, '
4471 WRITE(6,*)
' HIJSFT not performed, but continue'
4472 WRITE(6,*)
' D1=',
d1,
' D2=',
d2,
' SW=',sw
4473 WRITE(6,*)
' HISTORY NFP5=',nfp(jp,5),
' NFT5=',nft(jt,5)
4474 WRITE(6,*)
' THIS COLLISON NFP5=',nfp5,
' NFT5=',nft5
4475 WRITE(6,*)
' # OF JET IN PROJ',npj(jp),
' IN TARG',ntj(jt)
4478 IF(ihpr2(10).EQ.0)
RETURN
4479 WRITE(6,*)
' (4)unable to choose process, but not harmful'
4480 WRITE(6,*)
' HIJSFT not performed, but continue'
4481 WRITE(6,*)
' PTP=',
sqrt(ptp2),
' PTT=',
sqrt(ptt2),
' SW=',sw
4482 WRITE(6,*)
' AMCUT=',amx,
' JP=',jp,
' JT=',jt
4483 WRITE(6,*)
' HISTORY NFP5=',nfp(jp,5),
' NFT5=',nft(jt,5)
4486 IF(ihpr2(10).EQ.0)
RETURN
4487 WRITE(6,*)
' energy partition failed(5),for limited try'
4488 WRITE(6,*)
' HIJSFT not performed, but continue'
4489 WRITE(6,*)
' NFP5=',nfp5,
' NFT5=',nft5
4490 WRITE(6,*)
' D1',
d1,
' X1(1-X2)',x1*(1.0-x2)
4491 WRITE(6,*)
' D2',
d2,
' X2(1-X1)',x2*(1.0-x1)
4495 IF(miss.LT.100) go to 30
4497 IF(ihpr2(10).EQ.0)
RETURN
4498 WRITE(6,*)
' ERROR OCCURED, HIJSFT NOT PERFORMED'
4499 WRITE(6,*)
' Abort this event'
4500 WRITE(6,*)
'MTP,PTP2',epp*epm,ptp2,
' MTT,PTT2',etp*etm,ptt2
4511 IF(rnid.GT.0.43478)
THEN
4513 IF(rnid.GT.0.86956) id=3
4521 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4528 & pt=hipr1(2)*
sqrt(-alog(
exp(-hipr1(8)**2/hipr1(2)**2)
4529 & -
rlu(0)*(
exp(-hipr1(8)**2/hipr1(2)**2)-
4530 &
exp(-ptmax**2/hipr1(2)**2))))
4533 pt=hipr1(2)*
sqrt(-alog(1.0-
rlu(0)*
4534 & (1.0-
exp(-ptmax**2/hipr1(2)**2))))
4536 ptmax0=max(ptmax,0.01)
4537 pt=min(ptmax0-0.01,pt)
4550 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4554 dimension iaa(20),rr(20),dd(20),ww(20),
rms(20)
4559 DATA iaa/2,4,12,16,27,32,40,56,63,93,184,197,208,7*0./
4560 DATA rr/0.01,.964,2.355,2.608,2.84,3.458,3.766,3.971,4.214,
4561 1 4.87,6.51,6.38,6.624,7*0./
4562 DATA dd/0.5882,.322,.522,.513,.569,.61,.586,.5935,.586,.573,
4563 1 .535,.535,.549,7*0./
4564 DATA ww/0.0,.517,-0.149,-0.051,0.,-0.208,-0.161,13*0./
4565 DATA rms/2.11,1.71,2.46,2.73,3.05,3.247,3.482,3.737,3.925,4.31,
4566 1 5.42,5.33,5.521,7*0./
4568 SAVE iaa, rr, dd, ww,
rms
4575 r=1.19*
a**(1./3.) - 1.61*
a**(-1./3.)
4582 IF (ia.EQ.iaa(i))
THEN
4593 IF (w.LT.-0.01)
THEN
4594 IF (xhigh.GT.
r/
sqrt(abs(w))) xhigh=
r/
sqrt(abs(w))
4603 hint1(75)=fnorm/4.0/hipr1(40)
4604 ELSE IF (idh.EQ.2)
THEN
4608 hint1(79)=fnorm/4.0/hipr1(40)
4643 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4645 wdsax1=hint1(75)*(1.+hint1(74)*(x/hint1(72))**2)/
4646 & (1+
exp((x-hint1(72))/hint1(73)))
4647 IF (hint1(74).LT.0.)
THEN
4648 IF (x.GE.hint1(72)/
sqrt(abs(hint1(74))))
wdsax1=0.
4659 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4661 wdsax2=hint1(79)*(1.+hint1(78)*(x/hint1(76))**2)/
4662 & (1+
exp((x-hint1(76))/hint1(77)))
4663 IF (hint1(78).LT.0.)
THEN
4664 IF (x.GE.hint1(76)/
sqrt(abs(hint1(78))))
wdsax2=0.
4676 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4682 IF(ihnt2(1).GT.1 .AND. ihnt2(3).GT.1)
THEN
4683 profile=float(ihnt2(1))*float(ihnt2(3))*0.1*
4685 ELSE IF(ihnt2(1).EQ.1 .AND. ihnt2(3).GT.1)
THEN
4688 ELSE IF(ihnt2(1).GT.1 .AND. ihnt2(3).EQ.1)
THEN
4699 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4710 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4721 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4732 r1=
sqrt(b1**2+z1**2)
4786 10
IF(ju-jl.GT.1)
THEN
4788 IF((rr(i,201).GT.rr(i,1)).EQV.(rx.GT.rr(i,jm)))
THEN
4813 jmin=1+200*(
xmin-
xx(i,1))/(
xx(i,201)-
xx(i,1))
4815 rx=rr(i,jmin)+(rr(i,
jmax)-rr(i,jmin))*
rlu(0)
4818 10
IF(ju-jl.GT.1)
THEN
4820 IF((rr(i,201).GT.rr(i,1)).EQV.(rx.GT.rr(i,jm)))
THEN
4840 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4845 IF(hint1(1).GE.10.0) CALL
crsjet
4848 aphx1=hipr1(6)*(ihnt2(1)**0.3333333-1.0)
4849 aphx2=hipr1(6)*(ihnt2(3)**0.3333333-1.0)
4850 hint1(11)=hint1(14)-aphx1*hint1(15)
4851 & -aphx2*hint1(16)+aphx1*aphx2*hint1(17)
4855 hint1(60)=hint1(61)-aphx1*hint1(62)
4856 & -aphx2*hint1(63)+aphx1*aphx2*hint1(64)
4858 IF(hint1(59).EQ.0.0) hint1(59)=hint1(60)
4859 IF(hint1(1).GE.10.0)
Then
4865 hint1(10)=hint1(10)*hipr1(31)
4866 hint1(12)=hint1(12)*hipr1(31)
4867 hint1(13)=hint1(13)*hipr1(31)
4868 hint1(59)=hint1(59)*hipr1(31)
4871 IF(ihpr2(13).NE.0)
THEN
4872 hipr1(33)=1.36*(1.0+36.0/hint1(1)**2)
4873 & *alog(0.6+0.1*hint1(1)**2)
4874 hipr1(33)=hipr1(33)/hint1(12)
4885 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4887 omg=
omg0(x)*(hipr1(30)+hint1(11))/hipr1(31)/2.0
4895 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4897 omg=
omg0(x)*(hipr1(30)+hint1(11))/hipr1(31)/2.0
4905 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4907 omg=
omg0(x)*hint1(11)/hipr1(31)/2.0
4915 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4917 omg=
omg0(x)*hint1(60)/hipr1(31)/2.0
4926 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4930 omg1=
omg0(x)*hint1(11)/hipr1(31)
4932 IF(
n.EQ.0) c0=1.0-
exp(-2.0*
omg0(x)*hipr1(30)/hipr1(31)/2.0)
4955 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4960 x4=hipr1(32)*
sqrt(x)
4971 dimension fr(0:1000)
4979 IF(i0.NE.0) go to 100
4993 romg=(fr(ix)*((ix+1)*0.01-x)+fr(ix+1)*(x-ix*0.01))/0.01
5004 bk=
exp(-x)*(x**2-x4**2)**2.50/15.0
5013 IMPLICIT REAL*8(
a-h,o-z)
5014 REAL hipr1(100),hint1(100)
5015 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5019 common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn
5021 common/bveg2/xi(50,10),si,si2,swgt,schi,ndo,it
5037 hint1(14)=avgi/2.5682
5038 IF(ihpr2(6).EQ.1 .AND. ihnt2(1).GT.1)
THEN
5041 hint1(15)=avgi/2.5682
5043 IF(ihpr2(6).EQ.1 .AND. ihnt2(3).GT.1)
THEN
5046 hint1(16)=avgi/2.5682
5048 IF(ihpr2(6).EQ.1.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1)
THEN
5051 hint1(17)=avgi/2.5682
5055 IF(ihpr2(3).NE.0)
THEN
5058 hint1(61)=avgi/2.5682
5059 IF(ihpr2(6).EQ.1 .AND. ihnt2(1).GT.1)
THEN
5062 hint1(62)=avgi/2.5682
5064 IF(ihpr2(6).EQ.1 .AND. ihnt2(3).GT.1)
THEN
5067 hint1(63)=avgi/2.5682
5069 IF(ihpr2(6).EQ.1.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1)
THEN
5072 hint1(64)=avgi/2.5682
5083 IMPLICIT REAL*8(
a-h,o-z)
5084 REAL hipr1(100),hint1(100)
5085 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5088 pt2=(hint1(1)**2/4.0-hipr1(8)**2)*x(1)+hipr1(8)**2
5089 xt=2.0*dsqrt(pt2)/hint1(1)
5090 ymx1=dlog(1.0/
xt+dsqrt(1.0/
xt**2-1.0))
5091 y1=2.0*ymx1*x(2)-ymx1
5092 ymx2=dlog(2.0/
xt-dexp(y1))
5093 ymn2=dlog(2.0/
xt-dexp(-y1))
5094 y2=(ymx2+ymn2)*x(3)-ymn2
5095 fjet=2.0*ymx1*(ymx2+ymn2)*(hint1(1)**2/4.0-hipr1(8)**2)
5103 IMPLICIT REAL*8(
a-h,o-z)
5104 REAL hipr1(100),hint1(100),ptmax,ptmin
5105 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5108 ptmin=abs(hipr1(10))-0.25
5109 ptmin=max(ptmin,hipr1(8))
5111 IF(ihpr2(3).EQ.3)
THEN
5113 ptmin=max(0.0,hipr1(10))
5115 ptmax=abs(hipr1(10))+0.25
5116 IF(hipr1(10).LE.0.0) ptmax=hint1(1)/2.0-am2
5117 IF(ptmax.LE.ptmin) ptmax=ptmin+0.25
5118 pt2=(ptmax**2-ptmin**2)*x(1)+ptmin**2
5120 xt=2.0*dsqrt(amt2)/hint1(1)
5121 ymx1=dlog(1.0/
xt+dsqrt(1.0/
xt**2-1.0))
5122 y1=2.0*ymx1*x(2)-ymx1
5123 ymx2=dlog(2.0/
xt-dexp(y1))
5124 ymn2=dlog(2.0/
xt-dexp(-y1))
5125 y2=(ymx2+ymn2)*x(3)-ymn2
5126 IF(ihpr2(3).EQ.3)
THEN
5127 gtrig=2.0*
ghvq(y1,y2,amt2)
5128 ELSE IF(ihpr2(3).EQ.2)
THEN
5133 fjetrig=2.0*ymx1*(ymx2+ymn2)*(ptmax**2-ptmin**2)
5141 IMPLICIT REAL*8 (
a-h,o-z)
5142 REAL hipr1(100),hint1(100)
5143 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5146 xt=2.0*dsqrt(amt2)/hint1(1)
5147 x1=0.50*
xt*(dexp(y1)+dexp(y2))
5148 x2=0.50*
xt*(dexp(-y1)+dexp(-y2))
5149 ss=x1*x2*hint1(1)**2
5151 IF(ihpr2(18).NE.0) af=5.0
5153 aph=12.0*3.1415926/(33.0-2.0*af)/dlog(amt2/dlam**2)
5157 gqq=4.0*(cosh(y1-y2)+hipr1(7)**2/amt2)/(1.d0+cosh(y1-y2))/9.0
5158 & *(
f(1,1)*
f(2,2)+
f(1,2)*
f(2,1)+
f(1,3)*
f(2,4)
5159 & +
f(1,4)*
f(2,3)+
f(1,5)*
f(2,6)+
f(1,6)*
f(2,5))
5160 ggg=(8.d0*cosh(y1-y2)-1.d0)*(cosh(y1-y2)+2.0*hipr1(7)**2/amt2
5161 & -2.0*hipr1(7)**4/amt2**2)/(1.0+cosh(y1-y2))/24.d0
5164 ghvq=(gqq+ggg)*hipr1(23)*3.14159*aph**2/ss**2
5171 IMPLICIT REAL*8 (
a-h,o-z)
5172 REAL hipr1(100),hint1(100)
5173 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5176 xt=2.0*dsqrt(pt2)/hint1(1)
5177 x1=0.50*
xt*(dexp(y1)+dexp(y2))
5178 x2=0.50*
xt*(dexp(-y1)+dexp(-y2))
5179 z=dsqrt(1.d0-
xt**2/x1/x2)
5180 ss=x1*x2*hint1(1)**2
5185 aph=12.0*3.1415926/(33.0-2.0*af)/dlog(pt2/dlam**2)
5190 g11=-(u**2+1.0)/u/3.0*
f(1,7)*(4.0*
f(2,1)+4.0*
f(2,2)
5191 & +
f(2,3)+
f(2,4)+
f(2,5)+
f(2,6))/9.0
5192 g12=-(
t**2+1.0)/
t/3.0*
f(2,7)*(4.0*
f(1,1)+4.0*
f(1,2)
5193 & +
f(1,3)+
f(1,4)+
f(1,5)+
f(1,6))/9.0
5194 g2=8.0*(u**2+
t**2)/u/
t/9.0*(4.0*
f(1,1)*
f(2,2)
5195 & +4.0*
f(1,2)*
f(2,1)+
f(1,3)*
f(2,4)+
f(1,4)*
f(2,3)
5196 & +
f(1,5)*
f(2,6)+
f(1,6)*
f(2,5))/9.0
5198 gphoton=(g11+g12+g2)*hipr1(23)*3.14159*aph*aphem/ss**2
5205 FUNCTION g(Y1,Y2,PT2)
5206 IMPLICIT REAL*8 (
a-h,o-z)
5207 REAL hipr1(100),hint1(100)
5208 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5211 xt=2.0*dsqrt(pt2)/hint1(1)
5212 x1=0.50*
xt*(dexp(y1)+dexp(y2))
5213 x2=0.50*
xt*(dexp(-y1)+dexp(-y2))
5214 z=dsqrt(1.d0-
xt**2/x1/x2)
5215 ss=x1*x2*hint1(1)**2
5220 aph=12.0*3.1415926/(33.0-2.0*af)/dlog(pt2/dlam**2)
5224 g11=( (
f(1,1)+
f(1,2))*(
f(2,3)+
f(2,4)+
f(2,5)+
f(2,6))
5227 g12=( (
f(2,1)+
f(2,2))*(
f(1,3)+
f(1,4)+
f(1,5)+
f(1,6))
5230 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)
5234 g2=(af-1)*(
f(1,1)*
f(2,2)+
f(2,1)*
f(1,2)+
f(1,3)*
f(2,4)
5240 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)+
5245 g61=
f(1,7)*(
f(2,1)+
f(2,2)+
f(2,3)+
f(2,4)+
f(2,5)
5247 g62=
f(2,7)*(
f(1,1)+
f(1,2)+
f(1,3)+
f(1,4)+
f(1,5)
5252 g=(g11+g12+g13+g2+g31+g32+g4+g5+g61+g62+g7)*hipr1(17)*
5253 1 3.14159d0*aph**2/ss**2
5260 IMPLICIT REAL*8(
a-h,o-z)
5267 IMPLICIT REAL*8(
a-h,o-z)
5274 IMPLICIT REAL*8(
a-h,o-z)
5275 subcrs3=4.d0/9.d0*(
t**2+u**2+(1.d0+u**2)/
t**2
5276 1 -2.d0*u**2/3.d0/
t)
5282 IMPLICIT REAL*8(
a-h,o-z)
5283 subcrs4=8.d0/3.d0*(
t**2+u**2)*(4.d0/9.d0/
t/u-1.d0)
5290 IMPLICIT REAL*8(
a-h,o-z)
5291 subcrs5=3.d0/8.d0*(
t**2+u**2)*(4.d0/9.d0/
t/u-1.d0)
5297 IMPLICIT REAL*8(
a-h,o-z)
5298 subcrs6=(1.d0+u**2)*(1.d0/
t**2-4.d0/u/9.d0)
5304 IMPLICIT REAL*8(
a-h,o-z)
5312 IMPLICIT REAL*8(
a-h,o-z)
5313 REAL hipr1(100),hint1(100)
5314 common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5321 s=dlog(dlog(qq/dlam**2)/dlog(q0**2/dlam**2))
5322 IF(ihpr2(7).EQ.2) go to 200
5324 at1=0.419+0.004*
s-0.007*
s**2
5325 at2=3.460+0.724*
s-0.066*
s**2
5326 gmud=4.40-4.86*
s+1.33*
s**2
5327 at3=0.763-0.237*
s+0.026*
s**2
5328 at4=4.00+0.627*
s-0.019*
s**2
5329 gmd=-0.421*
s+0.033*
s**2
5331 cas=1.265-1.132*
s+0.293*
s**2
5332 as=-0.372*
s-0.029*
s**2
5333 bs=8.05+1.59*
s-0.153*
s**2
5334 aphs=6.31*
s-0.273*
s**2
5335 btas=-10.5*
s-3.17*
s**2
5336 gms=14.7*
s+9.80*
s**2
5345 cag=1.56-1.71*
s+0.638*
s**2
5346 ag=-0.949*
s+0.325*
s**2
5347 bg=6.0+1.44*
s-1.05*
s**2
5348 aphg=9.0-7.19*
s+0.255*
s**2
5349 btag=-16.5*
s+10.9*
s**2
5350 gmg=15.3*
s-10.1*
s**2
5353 200 at1=0.374+0.014*
s
5354 at2=3.33+0.753*
s-0.076*
s**2
5355 gmud=6.03-6.22*
s+1.56*
s**2
5356 at3=0.761-0.232*
s+0.023*
s**2
5357 at4=3.83+0.627*
s-0.019*
s**2
5358 gmd=-0.418*
s+0.036*
s**2
5360 cas=1.67-1.92*
s+0.582*
s**2
5361 as=-0.273*
s-0.164*
s**2
5362 bs=9.15+0.530*
s-0.763*
s**2
5363 aphs=15.7*
s-2.83*
s**2
5364 btas=-101.0*
s+44.7*
s**2
5365 gms=223.0*
s-117.0*
s**2
5374 cag=0.879-0.971*
s+0.434*
s**2
5375 ag=-1.16*
s+0.476*
s**2
5376 bg=4.0+1.23*
s-0.254*
s**2
5377 aphg=9.0-5.64*
s-0.817*
s**2
5378 btag=-7.54*
s+5.50*
s**2
5379 gmg=-0.596*
s+1.26*
s**2
5381 300 b12=dexp(
gmre(at1)+
gmre(at2+1.d0)-
gmre(at1+at2+1.d0))
5383 cnud=3.d0/b12/(1.d0+gmud*at1/(at1+at2+1.d0))
5384 cnd=1.d0/b34/(1.d0+gmd*at3/(at3+at4+1.d0))
5389 fud1=cnud*x1**at1*(1.d0-x1)**at2*(1.d0+gmud*x1)
5390 fs1=cas*x1**as*(1.d0-x1)**bs*(1.d0+aphs*x1
5391 & +btas*x1**2+gms*x1**3)
5392 f(1,3)=cnd*x1**at3*(1.d0-x1)**at4*(1.d0+gmd*x1)+fs1/6.d0
5393 f(1,1)=fud1-
f(1,3)+fs1/3.d0
5398 f(1,7)=cag*x1**ag*(1.d0-x1)**bg*(1.d0+aphg*x1
5399 & +btag*x1**2+gmg*x1**3)
5401 fud2=cnud*x2**at1*(1.d0-x2)**at2*(1.d0+gmud*x2)
5402 fs2=cas*x2**as*(1.d0-x2)**bs*(1.d0+aphs*x2
5403 & +btas*x2**2+gms*x2**3)
5404 f(2,3)=cnd*x2**at3*(1.d0-x2)**at4*(1.d0+gmd*x2)+fs2/6.d0
5405 f(2,1)=fud2-
f(2,3)+fs2/3.d0
5410 f(2,7)=cag*x2**ag*(1.d0-x2)**bg*(1.d0+aphg*x2
5411 & +btag*x2**2+gmg*x2**3)
5414 IF(ihpr2(6).EQ.1 .AND. ihnt2(1).GT.1)
THEN
5415 aax=1.193*alog(float(ihnt2(1)))**0.16666666
5416 rrx=aax*(x1**3-1.2*x1**2+0.21*x1)+1.0
5417 & +1.079*(float(ihnt2(1))**0.33333333-1.0)
5418 & /dlog(ihnt2(1)+1.0d0)*dsqrt(x1)*dexp(-x1**2/0.01)
5419 IF(ip_crs.EQ.1 .OR.ip_crs.EQ.3) rrx=dexp(-x1**2/0.01)
5424 IF(ihpr2(6).EQ.1 .AND. ihnt2(3).GT.1)
THEN
5425 aax=1.193*alog(float(ihnt2(3)))**0.16666666
5426 rrx=aax*(x2**3-1.2*x2**2+0.21*x2)+1.0
5427 & +1.079*(float(ihnt2(3))**0.33333-1.0)
5428 & /dlog(ihnt2(3)+1.0d0)*dsqrt(x2)*dexp(-x2**2/0.01)
5429 IF(ip_crs.EQ.2 .OR. ip_crs.EQ.3) rrx=dexp(-x2**2/0.01)
5441 IMPLICIT REAL*8(
a-h,o-z)
5443 IF(x.GT.3.0d0) go to 10
5445 10
gmre=0.5d0*dlog(2.d0*3.14159265d0/z)+z*dlog(z)-z+dlog(1.d0
5446 1 +1.d0/12.d0/z+1.d0/288.d0/z**2-139.d0/51840.d0/z**3
5447 1 -571.d0/2488320.d0/z**4)
5449 gmre=
gmre-dlog(z-1.d0)-dlog(z-2.d0)-dlog(z-3.d0)
5457 IMPLICIT REAL*8(
a-h,o-z)
5472 REAL*8 xl(10),xu(10),acc
5473 common/bveg1/xl,xu,acc,ndim,ncall,itmx,nprn
5477 common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
5481 common/himain1/ natt,eatt,jatt,nt,np,n0,n01,n10,n11
5483 common/himain2/katt(130000,4),patt(130000,4)
5485 common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
5489 common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
5490 & pjpy(300,500),pjpz(300,500),pjpe(300,500),
5491 & pjpm(300,500),ntj(300),kftj(300,500),
5492 & pjtx(300,500),pjty(300,500),pjtz(300,500),
5493 & pjte(300,500),pjtm(300,500)
5495 common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
5496 & k2sg(900,100),pxsg(900,100),pysg(900,100),
5497 & pzsg(900,100),pesg(900,100),pmsg(900,100)
5499 common/hijdat/hidat0(10,10),hidat(10)
5501 common/hipyint/mint4,mint5,atco(200,20),atxs(0:200)
5503 DATA num1/30123984/,xl/10*0.d0/,xu/10*1.d0/
5504 DATA ncall/1000/itmx/100/acc/0.01/nprn/0/
5507 DATA nseed/74769375/
5509 & 1.5, 0.35, 0.5, 0.9, 2.0, 0.1, 1.5, 2.0, -1.0, -2.25,
5510 & 2.0, 0.5, 1.0, 2.0, 0.2, 2.0, 2.5, 0.3, 0.1, 1.4,
5511 & 1.6, 1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 57.0,
5512 & 28.5, 3.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
5514 & 0.0, 0.4, 0.1, 1.5, 0.1, 0.25, 0.0, 0.5, 0.0, 0.0,
5518 & 1, 3, 0, 1, 1, 1, 1, 10, 0, 0,
5519 & 1, 1, 1, 1, 0, 0, 1, 0, 0, 1,
5526 DATA natt/0/eatt/0.0/jatt/0/nt/0/np/0/n0/0/n01/0/n10/0/n11/0/
5527 DATA katt/520000*0/patt/520000*0.0/
5529 DATA nfp/4500*0/pp/4500*0.0/nft/4500*0/pt/4500*0.0/
5531 DATA yp/900*0.0/
yt/900*0.0/
5533 DATA npj/300*0/kfpj/150000*0/pjpx/150000*0.0/pjpy/150000*0.0/
5534 & pjpz/150000*0.0/pjpe/150000*0.0/pjpm/150000*0.0/
5535 DATA ntj/300*0/kftj/150000*0/pjtx/150000*0.0/pjty/150000*0.0/
5536 & pjtz/150000*0.0/pjte/150000*0.0/pjtm/150000*0.0/
5538 DATA nsg/0/njsg/900*0/iasg/2700*0/k1sg/90000*0/k2sg/90000*0/
5539 & pxsg/90000*0.0/pysg/90000*0.0/pzsg/90000*0.0/pesg/90000*0.0/
5541 DATA mint4/0/mint5/0/atco/4000*0.0/atxs/201*0.0/
5542 DATA (hidat0(1,i),i=1,10)/0.0,0.0,0.0,0.0,0.0,0.0,2.25,
5544 DATA (hidat0(2,i),i=1,10)/2.0,3.0,5.0,6.0,7.0,8.0,8.0,10.0,
5546 DATA (hidat0(3,i),i=1,10)/1.0,0.8,0.8,0.7,0.45,0.215,
5547 & 0.21,0.19,0.19,0.19/
5548 DATA (hidat0(4,i),i=1,10)/0.35,0.35,0.3,0.3,0.3,0.3,
5550 DATA (hidat0(5,i),i=1,10)/23.8,24.0,26.0,26.2,27.0,28.5,28.5,
5552 DATA ((hidat0(j,i),i=1,10),j=6,9)/40*0.0/
5553 DATA (hidat0(10,i),i=1,10)/5.0,20.0,53.0,62.0,100.0,200.0,
5554 & 546.0,900.0,1800.0,4000.0/
5568 IMPLICIT REAL*8(
a-h,o-z)
5569 common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn
5571 common/bveg2/xi(50,10),si,si2,swgt,schi,ndo,it
5576 dimension
d(50,10),di(50,10),xin(50),
r(50),
dx(10),dt(10),x(10)
5579 DATA ndmx/50/,alph/1.5d0/,
one/1.d0/,mds/-1/
5581 save ndmx, alph,
one, mds
5587 entry vegas1(fxn,avgi,sd,chi2a)
5595 entry vegas2(fxn,avgi,sd,chi2a)
5599 IF(mds.EQ.0) go to 2
5600 ng=(ncall/2.)**(1./ndim)
5602 IF((2*ng-ndmx).LT.0) go to 2
5612 dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-
one)
5624 IF(nd.EQ.ndo) go to 8
5635 5
IF(rc.GT.dr) go to 4
5638 xin(i)=xn-(xn-xo)*dr
5639 IF(i.LT.ndm) go to 5
5646 IF(nprn.NE.0)
WRITE(16,200) ndim,calls,it,itmx,acc,mds,nd
5647 1 ,(xl(j),xu(j),j=1,ndim)
5649 entry vegas3(fxn,avgi,sd,chi2a)
5664 CALL
aran9(qran,ndim)
5667 xn=(
kg(j)-qran(j))*dxg+
one
5670 IF(ia(j).GT.1) go to 13
5674 13 xo=xi(ia(j),j)-xi(ia(j)-1,j)
5675 rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
5676 14 x(j)=xl(j)+rc*
dx(j)
5686 di(ia(j),j)=di(ia(j),j)+
f
5687 16
IF(mds.GE.0)
d(ia(j),j)=
d(ia(j),j)+f2
5688 IF(k.LT.npg) go to 12
5691 f2b=(f2b-fb)*(f2b+fb)
5694 IF(mds.GE.0) go to 18
5696 17
d(ia(j),j)=
d(ia(j),j)+f2b
5699 IF(
kg(k).NE.1) go to 11
5707 wgt=ti2/(tsi+1.0
d-37)
5716 chi2a=sd*(schi/swgt-avgi*avgi)/(it-.999)
5719 IF(nprn.EQ.0) go to 21
5721 WRITE(16,201) it,ti,tsi,avgi,sd,chi2a
5722 IF(nprn.GE.0) go to 21
5724 20
WRITE(16,202) j,(xi(i,j),di(i,j),
d(i,j),i=1,nd)
5737 d(i,j)=(
d(i,j)+xn)/3.
5738 22 dt(j)=dt(j)+
d(i,j)
5740 23 dt(j)=dt(j)+
d(nd,j)
5746 IF (dt(j).GE.1.0d18)
THEN
5747 WRITE(6,*)
'************** A SINGULARITY >1.0D18'
5757 IF(
d(i,j).LE.1.0
d-18) go to 24
5759 r(i)=((xo-
one)/xo/dlog(xo))**alph
5771 26
IF(rc.GT.dr) go to 25
5774 xin(i)=xn-(xn-xo)*dr/(
r(k)+1.0
d-30)
5775 IF(i.LT.ndm) go to 26
5780 IF(it.LT.itmx.AND.acc*dabs(avgi).LT.sd) go to 9
5781 200
FORMAT(
'0INPUT PARAMETERS FOR VEGAS: NDIM=',i3,
' NCALL=',f8.0
5782 1 /28x,
' IT=',i5,
' ITMX=',i5/28x,
' ACC=',g9.3
5783 2 /28x,
' MDS=',i3,
' ND=',i4/28x,
' (XL,XU)=',
5784 3 (t40,
'( ',g12.6,
' , ',g12.6,
' )'))
5785 201
FORMAT(///
' INTEGRATION BY VEGAS' /
'0ITERATION NO.',i3,
5786 1
': INTEGRAL =',g14.8/21x,
'STD DEV =',g10.4 /
5787 2
' ACCUMULATED RESULTS: INTEGRAL =',g14.8 /
5788 3 24x,
'STD DEV =',g10.4 / 24x,
'CHI**2 PER IT''N =',g10.4)
5789 202
FORMAT(
'0DATA FOR AXIS',i2 /
' ',6x,
'X',7x,
' DELT I ',
5790 1 2x,
' CONV''CE ',11x,
'X',7x,
' DELT I ',2x,
' CONV''CE '
5791 2 ,11x,
'X',7x,
' DELT I ',2x,
' CONV''CE ' /
5792 2 (
' ',3g12.4,5x,3g12.4,5x,3g12.4))
5811 dimension w(12),x(12)
5813 DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5814 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5815 & .1826034,.1894506/
5816 DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5817 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5818 & .2816036,.0950125/
5823 IF(abs(
y).LE.
delta)
RETURN
5831 1 s8=s8+w(i)*(
f(c1+u)+
f(c1-u))
5834 3 s16=s16+w(i)*(
f(c1+u)+
f(c1-u))
5837 IF(abs(s16-s8).GT.
eps*(1.+abs(s16))) goto 4
5842 IF(abs(
y).GT.
delta) goto 2
5846 7
FORMAT(1x,
'GAUSS1....TOO HIGH ACURACY REQUIRED')
5853 dimension w(12),x(12)
5855 DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5856 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5857 & .1826034,.1894506/
5858 DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5859 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5860 & .2816036,.0950125/
5865 IF(abs(
y).LE.
delta)
RETURN
5873 1 s8=s8+w(i)*(
f(c1+u)+
f(c1-u))
5876 3 s16=s16+w(i)*(
f(c1+u)+
f(c1-u))
5879 IF(abs(s16-s8).GT.
eps*(1.+abs(s16))) goto 4
5884 IF(abs(
y).GT.
delta) goto 2
5888 7
FORMAT(1x,
'GAUSS2....TOO HIGH ACURACY REQUIRED')
5895 dimension w(12),x(12)
5897 DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5898 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5899 & .1826034,.1894506/
5900 DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5901 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5902 & .2816036,.0950125/
5907 IF(abs(
y).LE.
delta)
RETURN
5915 1 s8=s8+w(i)*(
f(c1+u)+
f(c1-u))
5918 3 s16=s16+w(i)*(
f(c1+u)+
f(c1-u))
5921 IF(abs(s16-s8).GT.
eps*(1.+abs(s16))) goto 4
5926 IF(abs(
y).GT.
delta) goto 2
5930 7
FORMAT(1x,
'GAUSS3....TOO HIGH ACURACY REQUIRED')
5938 dimension w(12),x(12)
5940 DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5941 & .0622535,0.0951585,.1246290,.1495960,.1691565,
5942 & .1826034,.1894506/
5943 DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5944 & .9445750,0.8656312,.7554044,.6178762,.4580168,
5945 & .2816036,.0950125/
5950 IF(abs(
y).LE.
delta)
RETURN
5958 1 s8=s8+w(i)*(
f(c1+u)+
f(c1-u))
5961 3 s16=s16+w(i)*(
f(c1+u)+
f(c1-u))
5964 IF(abs(s16-s8).GT.
eps*(1.+abs(s16))) goto 4
5969 IF(abs(
y).GT.
delta) goto 2
5973 7
FORMAT(1x,
'GAUSS4....TOO HIGH ACURACY REQUIRED')
5983 &
'**************************************************'/10x,
5984 &
'* | \\ _______ / ------/ *'/10x,
5985 &
'* ----- ------ |_____| /_/ / *'/10x,
5986 &
'* ||\\ / |_____| / / \\ *'/10x,
5987 &
'* /| \\ /_/ /_______ /_ / \\_ *'/10x,
5988 &
'* / | / / / / / | ------- *'/10x,
5989 &
'* | / /\\ / / | / | *'/10x,
5990 &
'* | / / \\ / / \\_| / ------- *'/10x,
5992 &
'**************************************************'/10x,
5994 &
' Heavy Ion Jet INteraction Generator '/10x,
5996 &
' X. N. Wang and M. Gyulassy '/10x,
5997 &
' Lawrence Berkeley Laboratory '//)
static c2_factory< G4double > c2
G4int nint(G4double number)
subroutine ar3jet(S, X1, X3, JL)
subroutine hifun(I, XMIN, XMAX, FHB)
subroutine pyinit(FRAME, BEAM, TARGET, WIN)
std::vector< ExP01TrackerHit * > a
subroutine hijfrg(JTP, NTP, IERROR)
subroutine hijwds(IA, IDH, XHIGH)
function gauss1(F, A, B, EPS)
function gauss3(F, A, B, EPS)
function gphoton(Y1, Y2, PT2)
subroutine aran9(QRAN, NDIM)
static const G4double eps
static int FASTCALL common(PROLOG_STATE *state, int tok)
subroutine hijcsc(JP, JT)
static constexpr double g
G4int mod(G4int a, G4int b)
subroutine vegas(FXN, AVGI, SD, CHI2A)
subroutine hijing(BMIN0, BMAX0)
static constexpr double m
function gauss2(F, A, B, EPS)
subroutine atrobo(THE, PHI, BEX, BEY, BEZ, IMIN, IMAX, IERROR)
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
function hirnd2(I, XMIN, XMAX)
subroutine hijels(PSC1, PSC2)
static constexpr double kg
subroutine attrad(IERROR)
subroutine hiptdi(PT, PTMAX, IOPT)
subroutine quench(JPJT, NTP)
subroutine hijhrd(JP, JT, JOUT, JFLG, IOPJET0)
function gauss4(F, A, B, EPS)
subroutine hijsrt(JPJT, NPT)
static c2_log_p< float_type > & log()
make a *new object
subroutine arorie(S, X1, X3, JL)
subroutine attflv(ID, IDQ, IDQQ)
subroutine hirobo(THE, PHI, BEX, BEY, BEZ)
static c2_sqrt_p< float_type > & sqrt()
make a *new object
static c2_cos_p< float_type > & cos()
make a *new object
function ghvq(Y1, Y2, AMT2)
subroutine parton(F, X1, X2, QQ)
subroutine jetini(JP, JT, I_TRIG)
subroutine hijsft(JP, JT, JOUT, IERROR)
static c2_sin_p< float_type > & sin()
make a *new object
static c2_exp_p< float_type > & exp()
make a *new object