11 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
17 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
18 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
22 common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
23 & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
24 & xmi(2,240),pt2mi(240),imisep(0:240)
27 parameter(nersiz=4000)
28 COMMON /pycbls/mco(nersiz,2),ncc,jcco(nersiz,2),jccn(nersiz,2)
30 COMMON /pyctag/nct,mct(nersiz,2)
31 SAVE /pycbls/,/pyctag/
32 dimension jst(2,3),iv(2,3),idq(3),nvsum(2),nbrtot(2),ng(2)
33 & ,itjunc(2),mout(2),insr(1000,3),istr(6),ymi(240)
42 DO 110
i=
mint(84)+1,nersiz
49 k(
i,4)=mod(
k(
i,4),mstu(5)**2)
50 k(
i,5)=mod(
k(
i,5),mstu(5)**2)
56 DO 130 mg=
mint(84)+1,nersiz
71 IF (
mstp(89).EQ.1)
THEN
73 IF (im.LE.
mint(31))
THEN
74 ymi(im)=
log(xmi(1,im)/xmi(2,im))
88 kfs=isign(1,
mint(10+js))
91 IF(kfival(js,1).EQ.0)
THEN
92 IF(
mint(10+js).EQ.111)
THEN
93 kfival(js,1)=int(1.5d0+
pyr(0))
94 kfival(js,2)=-kfival(js,1)
95 ELSEIF(
mint(10+js).EQ.22)
THEN
98 IF(pyrkf.GT.0.1d0) kfival(js,1)=2
99 IF(pyrkf.GT.0.5d0) kfival(js,1)=3
100 IF(pyrkf.GT.0.6d0) kfival(js,1)=4
101 kfival(js,2)=-kfival(js,1)
102 ELSEIF(
mint(10+js).EQ.130.OR.
mint(10+js).EQ.310)
THEN
103 IF(
pyr(0).GT.0.5d0)
THEN
122 IF(kfival(js,
j).EQ.jfa) nvalq=nvalq+1
123 IF(kfival(js,
j).EQ.-jfa) nvalqb=nvalqb+1
125 nvsum(js)=nvsum(js)+nvalq+nvalqb
128 ifl =
k(imi(js,im,1),2)
131 IF (ifl.EQ.jfa.AND.imi(js,im,2).EQ.0)
THEN
134 jv=nvsum(js)-nvalq-nvalqb
135 iv(js,jv)=imi(js,im,1)
136 ELSEIF (ifl.EQ.-jfa.AND.imi(js,im,2).EQ.0)
THEN
139 jv=nvsum(js)-nvalq-nvalqb
140 iv(js,jv)=imi(js,im,1)
141 ELSEIF (ifa.EQ.jfa)
THEN
143 IF (imi(js,im,2).LT.0) nsea=nsea-ifs
147 nflsum=iabs(nsea)+nvalq+nvalqb
148 nbrtot(js)=nbrtot(js)+nflsum
149 IF (
n+nflsum+1.GT.mstu(4))
THEN
150 CALL
pyerrm(11,
'(PYMIHK:) no more memory left in PYJETS')
155 IF (nflsum.GT.0)
THEN
165 k(
n,2)=isign(jfa,nsea)
166 IF (ia.LE.nvalq)
k(
n,2)=jfa
167 IF (ia.GT.nvalq.AND.ia.LE.nvalq+nvalqb)
k(
n,2)=-jfa
173 IF (ia.LE.nvalq+nvalqb)
THEN
176 iv(js,jv)=imi(js,nmi(js),1)
184 IF (im.LE.nmi(js))
THEN
185 IF (
k(imi(js,im,1),2).EQ.21)
THEN
188 ELSEIF (imi(js,im,2).NE.0.AND.
k(imi(js,im,1),2).GT.0)
THEN
190 IF (imi(js,im,2).LT.0)
THEN
192 230 imc=mod(imc,nmi(js))+1
193 IF (
k(imi(js,imc,1),2).NE.-
k(imi(js,im,1),2)) goto 230
194 IF (imi(js,imc,2).GE.0) goto 230
195 imi(js, im,2) = imi(js,imc,1)
196 imi(js,imc,2) = imi(js, im,1)
212 k(imi(js, im,2),5)=
k(imi(js, im,2),5)+mstu(5)*
n
213 k(imi(js, im,1),4)=
k(imi(js, im,1),4)+mstu(5)*
n
225 IF (iabs(
mint(10+js)).GT.1000)
THEN
227 itjunc(js) = (3-kfs)/2
240 k(
n,4)=itjunc(js)*mstu(5)
247 k(iv(js,jv),4+manti)=mod(
k(iv(js,jv),4+manti),mstu(5))
255 IF (
k(iv(js,1),2).GT.0)
THEN
266 k(iq,4)=mod(
k(iq,4),mstu(5))+mstu(5)*iqbar
267 k(iqbar,5)=mod(
k(iqbar,5),mstu(5))+mstu(5)*iq
269 IF (nbrtot(js).EQ.0)
THEN
293 IF (jst(js,jv).LE.
mint(53).AND.jst(js,jv).GT.0)
294 & mout(js)=mout(js)+1
310 IF (
k(
i1,2).NE.21.AND.(9-2*jcs).NE.isign(1,
k(
i1,2)))
312 IF (
k(
i1,jcs)/mstu(5)**2.NE.0) goto 300
316 IF(
mint(51).NE.0)
RETURN
324 DO 330 im=
mint(31)+1,nmi(js)
326 IF (
k(ip,2).NE.21)
THEN
327 jc=(3-isign(1,
k(ip,2)))/2
328 IF (mct(ip,jc).EQ.0)
THEN
359 IF (
k(
i,1).NE.3) goto 350
363 IF(kq.EQ.0.OR.(mqgst.EQ.1.AND.kq.EQ.2)) goto 350
367 IF(kq*isign(1,
k(
i,2)).LT.0) kcs=5
368 IF(mct(
i,kcs-3).NE.0) goto 350
371 IF(
mint(51).NE.0)
RETURN
394 beta=(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
395 CALL
pyrobo(imisep(im-1)+1,imisep(im),0d0,0d0,0d0,0d0,
beta)
402 IF (
k(
i,2).NE.88)
THEN
405 IF (js.EQ.2)
p(
i,3)=-
p(
i,3)
415 k(
i,4)=mod(
k(
i,4),mstu(5)**2)
416 k(
i,5)=mod(
k(
i,5),mstu(5)**2)
422 IF (itjunc(js).EQ.0)
THEN
426 jcco(ncc,1)=
max(jc1,jc2)
427 jcco(ncc,2)=min(jc1,jc2)
430 IF (mct(
i,1).EQ.jcco(ncc,1)) mct(
i,1)=jcco(ncc,2)
431 IF (mct(
i,2).EQ.jcco(ncc,1)) mct(
i,2)=jcco(ncc,2)
437 IF (
pyr(0).GT.0.5d0.OR.ng(1).EQ.0) js=2
438 IF (ng(js).GT.0)
THEN
443 imgl=
pyr(0)*nmi(js)+1
444 450 imgl=mod(imgl,nmi(js))+1
447 IF (nmgl.LE.nmi(js).AND.nopt.LE.3)
THEN
450 IF (
k(igl,2).NE.21.OR.
k(igl,4)/mstu(5).NE.0
451 & .OR.
k(igl,5)/mstu(5).NE.0) goto 450
454 imp1=
pyr(0)*nmi(js)+1
455 460 imp1=mod(imp1,nmi(js))+1
457 IF (imp1.EQ.imgl) goto 460
459 IF (nmp1.LE.nmi(js).AND.nopt.LE.3)
THEN
465 470 manti=mod(manti+1,2)
468 ip2 =mod(
k(ip1,4+manti)/mstu(5),mstu(5))
471 IF (ip2.LE.0) goto 470
472 IF (
k(ip2,2).EQ.21.AND.ip2.GT.
mint(53)) goto 470
474 IF (
k(ip1,4+manti)/mstu(5)**2.EQ.2) goto 470
475 IF (
k(ip2,5-manti)/mstu(5)**2.EQ.2) goto 470
478 k(ip1,4+manti)=
k(ip1,4+manti)+2*mstu(5)**2
479 IF (
k(ip2,2).NE.88)
THEN
480 k(ip2,5-manti)=
k(ip2,5-manti)+2*mstu(5)**2
487 jcg1=mco(igl,2-manti)
488 jcg2=mco(igl,1+manti)
489 jcp1=mco(ip1,1+manti)
490 jcp2=mco(ip2,2-manti)
492 CALL
pymihg(jcp1,jcg1,jcp2,jcg2)
494 IF (maccpt.EQ.0) goto 470
497 jcg1=mct(igl,2-manti)
498 jcg2=mct(igl,1+manti)
499 jcp1=mct(ip1,1+manti)
500 jcp2=mct(ip2,2-manti)
503 IF (
mstp(89).EQ.0)
THEN
508 IF (ip1.GT.
mint(53).AND.ip2.GT.
mint(53)
509 & .AND.mout(js).NE.0.AND.
pyr(0).GT.parp80)
THEN
514 ELSEIF (
mstp(89).EQ.1)
THEN
519 IF (ygl.EQ.100d0)
THEN
521 ida1=mod(
k(igl,4),mstu(5))
522 ida2=mod(
k(igl,5),mstu(5))
525 IF (imi(js,imt,1).EQ.ida1.OR.imi(js,imt,1).EQ.ida2)
527 IF (abs(ygl).GT.abs(ymi(imt))) ygl=ymi(imt)
534 IF (yp1.EQ.100d0)
THEN
536 ida1=mod(
k(ip1,4),mstu(5))
537 ida2=mod(
k(ip1,5),mstu(5))
540 IF (imi(js,imt,1).EQ.ida1.OR.imi(js,imt,1).EQ.ida2)
542 IF (abs(yp1).GT.abs(ymi(imt))) yp1=ymi(imt)
547 IF (
k(ip2,2).NE.88)
THEN
549 IF (imi(js,imt,1).EQ.ip2) yp2=ymi(imt)
552 IF (yp2.EQ.100d0)
THEN
554 ida1=mod(
k(ip2,4),mstu(5))
555 ida2=mod(
k(ip2,5),mstu(5))
558 IF (imi(js,imt,1).EQ.ida1.OR.imi(js,imt,1).EQ.ida2
560 IF (abs(yp2).GT.abs(ymi(imt))) yp2=ymi(imt)
568 rl=abs(ygl-yp1)+abs(ygl-yp2)
569 ELSEIF (
mstp(89).EQ.2)
THEN
578 itju=mod(
k(ip2,4)/mstu(5),mstu(5))
583 IF (
k(
i,1).LT.10)
THEN
585 IF (mct(
i,1).EQ.jcg1) istr(1)=
i
586 IF (mct(
i,2).EQ.jcg1) istr(2)=
i
587 IF (mct(
i,1).EQ.jcg2) istr(3)=
i
588 IF (mct(
i,2).EQ.jcg2) istr(4)=
i
593 icmo=mod(
k(
i,4)/mstu(5),mstu(5))
594 iamo=mod(
k(
i,5)/mstu(5),mstu(5))
596 IF (
k(icmo,1).EQ.42.AND.mct(
i,1).EQ.jcg1.AND.istr(2)
597 & .EQ.0) istr(2) = icmo
598 IF (
k(iamo,1).EQ.42.AND.mct(
i,2).EQ.jcg1.AND.istr(1)
599 & .EQ.0) istr(1) = iamo
600 IF (
k(icmo,1).EQ.42.AND.mct(
i,1).EQ.jcg2.AND.istr(4)
601 & .EQ.0) istr(4) = icmo
602 IF (
k(iamo,1).EQ.42.AND.mct(
i,2).EQ.jcg2.AND.istr(3)
603 & .EQ.0) istr(3) = iamo
606 istr(5)=istr(1+2*manti)
607 istr(6)=istr(4-2*manti)
608 rl=
max(1d0,four(istr(1),istr(2)))*
max(1d0,four(istr(3)
609 & ,istr(4)))/
max(1d0,four(istr(5),istr(6)))
613 IF (abs(1d0-rl/rlopt).LT.0.05d0)
THEN
615 ELSEIF (rl.GT.rlopt)
THEN
624 IF (nopt.GT.1000) goto 470
625 insr(nopt,1+2*manti)=ip2
627 insr(nopt,3-2*manti)=ip1
628 IF (
mstp(89).GT.0.OR.nopt.EQ.0) goto 470
630 IF (
mstp(89).GT.0.OR.nopt.EQ.0) goto 460
634 k(
i,4)=mod(
k(
i,4),mstu(5)**2)
635 k(
i,5)=mod(
k(
i,5),mstu(5)**2)
637 IF (
mstp(89).GT.0.OR.nopt.EQ.0) goto 450
647 IF (mretry.LE.10.AND.(itjunc(1).NE.0.OR.jst(1,3).EQ.0).and
648 & .(itjunc(2).NE.0.OR.jst(2,3).EQ.0))
THEN
651 IF (itjunc(js).NE.0)
THEN
657 IF (
k(
i,2).EQ.88.AND.
k(
i,3).EQ.js) iju=
i
662 k(iv(js,jv),4+manti)=mod(
k(iv(js,jv),4+manti),mstu(5))
669 k(iq,4)=mod(
k(iq,4),mstu(5))+mstu(5)*iqbar
670 k(iqbar,5)=mod(
k(iqbar,5),mstu(5))+mstu(5)*iq
676 IF (
k(
i,2).EQ.21)
THEN
677 k(
i,4)=mod(
k(
i,4),mstu(5))
678 k(
i,5)=mod(
k(
i,5),mstu(5))
693 CALL
pyerrm(19,
'(PYMIHK:) No physical colour flow found!')
694 WRITE(mstu(11),*)
'NG:', ng,
' MOUT:', mout(js)
703 610 iin=mod(iin,nopt)+1
704 IF (insr(iin,1).GT.
mint(53).AND.insr(iin,3).GT.
mint(53)
705 & .AND.mout(js).NE.0.AND.
pyr(0).GT.parp80) goto 610
724 CALL
pymihg(jacm,jacg,jcm,jcg)
725 IF (maccpt.EQ.0)
THEN
729 CALL
pyerrm(11,
'(PYMIHK:) Unphysical colour flow!')
730 WRITE(mstu(11),*)
'attaching', igl,
' between', icm, iacm
739 jcco(icc,1)=jccn(icc,1)
740 jcco(icc,2)=jccn(icc,2)
747 k(igl,4)=mod(
k(igl,4),mstu(5))+mstu(5)*icm
749 IF (
k(icm,2).NE.88)
THEN
750 k(icm,5)=mod(
k(icm,5),mstu(5))+mstu(5)*igl
754 IF (jst(js,msj).EQ.iacm) jst(js,msj)=igl
758 k(igl,5)=mod(
k(igl,5),mstu(5))+mstu(5)*iacm
760 IF (
k(iacm,2).NE.88)
THEN
761 k(iacm,4)=mod(
k(iacm,4),mstu(5))+mstu(5)*igl
765 IF (jst(js,msj).EQ.icm) jst(js,msj)=igl
771 IF (ng(1).GT.0.OR.ng(2).GT.0) goto 440
775 DO 670 igl=
mint(53)+1,
n
776 IF (
k(igl,2).EQ.21.AND.
k(igl,3).EQ.
mint(83)+js.AND.
777 &
k(igl,1).EQ.14)
THEN
780 icd=mod(
k(igl,4),mstu(5))
781 iad=mod(
k(igl,5),mstu(5))
783 k(iad,5)=mod(
k(iad,5),mstu(5))+mstu(5)*iam
784 k(icd,4)=mod(
k(icd,4),mstu(5))+mstu(5)*icm
786 IF (
k(icm,2).NE.88)
THEN
787 k(icm,5)=mod(
k(icm,5),mstu(5))+mstu(5)*icd
791 IF (jst(js,msj).EQ.igl) jst(js,msj)=icd
794 IF (
k(iam,2).NE.88)
THEN
795 k(iam,4)=mod(
k(iam,4),mstu(5))+mstu(5)*iad
798 IF (jst(js,msj).EQ.igl) jst(js,msj)=iad
807 IF (im.GT.
mint(31).AND.
k(imi(js,im,1),2).NE.21) goto 680
808 IF (im.GT.
mint(31))
THEN
810 DO 690 imr=im,nmi(js)
811 imi(js,imr,1)=imi(js,imr+1,1)
812 imi(js,imr,2)=imi(js,imr+1,2)
818 IF (itjunc(js).NE.0)
THEN
820 IF (
k(
i,2).EQ.88.AND.
k(
i,3).EQ.
mint(83)+js) iju=
i
828 IF (jst(js,msj).GT.
mint(53).AND.iabs(
k(jst(js,msj),2)).LE.5)
832 idq(nbrjq)=-jst(js,msj)
834 IF (iv(js,jv).EQ.jst(js,msj))
THEN
835 idq(nbrjq)=jst(js,msj)
843 k(iju,i45)=
k(iju,i45)+(mstu(5)**i12)*jst(js,msj)
847 IF ((
mstp(88).GE.0.AND.nbrvq.GE.2).OR.(nbrjq.GE.2.AND.
mstp(88)
853 730 jflip=nbrjq*
pyr(0)+1
854 IF (idq(jflip).LT.0)
THEN
855 idq(jflip)=-idq(jflip)
858 IF (ndiq.LE.1) goto 730
862 IF (idq(jdq).LE.0)
THEN
866 IF (iabs(
k(idq(1),2)).LT.iabs(
k(idq(2),2)))
THEN
882 IF (idq(jdq).EQ.iv(js,jv)) jko=jko-jv
892 kfdq=1000*
k(idq(1),2)+100*
k(idq(2),2)
894 IF (
k(idq(1),2).NE.
k(idq(2),2).AND.
895 & (1d0+3d0*parj(4))*
pyr(0).LT.1d0) is=1
896 kfdq=kfdq+isign(is,kfdq)
902 770
IF (idq(3).NE.0.AND.
mstp(88).GE.2)
THEN
905 CALL
pykfdi(kfdq,
k(iabs(idq(3)),2),kfdum,kfbar)
906 IF (kfbar.EQ.0.AND.
ntry.LE.100)
THEN
908 ELSEIF(
ntry.GT.100)
THEN
937 IF (ip.NE.idq(1).AND.ip.NE.idq(2))
THEN
939 k(iju,5-manti)=ip*mstu(5)
940 k(ip,4+manti)=mod(
k(ip,4+manti),mstu(5))+
942 mct(iju,2-manti)=mct(ip,1+manti)
953 IF ((
k(
i,3).EQ.
mint(83)+js.OR.
k(
i,3).EQ.
mint(83)+2+js).and
956 imo=
k(
i,isid)/mstu(5)
957 ida=mod(
k(
i,isid),mstu(5))
959 IF (
k(imo,1).EQ.-1) imo=iju
962 IF (
k(ida,1).EQ.-1) ida=iju
964 k(
i,isid)=ida+mstu(5)*imo
974 IF (nbrtot(js).EQ.0)
THEN
985 IF (itjunc(js).NE.0)
THEN
987 jleg=
pyr(0)*nvsum(js)+1
990 jct=mct(
i1,itjunc(js))
991 mct(igl,3-itjunc(js))=jct
993 mct(igl,itjunc(js))=nct
997 CALL
pyerrm(19,
'(PYMIHK:) Empty meson beam remnant')
999 WRITE(mstu(11),*)
'This should not have been possible!'
1006 i2=mod(
k(
i1,4+manti)/mstu(5),mstu(5))
1007 k(
i1,4+manti)=mod(
k(
i1,4+manti),mstu(5))+mstu(5)*igl
1008 k(igl,5-manti)=mod(
k(igl,5-manti),mstu(5))+mstu(5)*
i1
1009 k(igl,4+manti)=mod(
k(igl,4+manti),mstu(5))+mstu(5)*
i2
1010 IF (
k(
i2,2).NE.88)
THEN
1011 k(
i2,5-manti)=mod(
k(
i2,5-manti),mstu(5))+mstu(5)*igl
1013 IF (mod(
k(
i2,4),mstu(5)).EQ.
i1)
THEN
1014 k(
i2,4)=(
k(
i2,4)/mstu(5))*mstu(5)+igl
1015 ELSEIF(mod(
k(
i2,5)/mstu(5),mstu(5)).EQ.
i1)
THEN
1016 k(
i2,5)=mod(
k(
i2,5),mstu(5))+mstu(5)*igl
1018 k(
i2,5)=(
k(
i2,5)/mstu(5))*mstu(5)+igl
1031 IF (
k(
i,1).LE.0) goto 850
1033 IF ((
k(
i,2).NE.21.OR.
k(
i,1).NE.14).AND.
k(
i,2).NE.88)
THEN
1043 IF (
k(
i,2).EQ.21.AND.
k(
i,1).EQ.14)
THEN
1045 jcd=mod(
k(
i,4),mstu(5))
1046 jad=mod(
k(
i,5),mstu(5))
1048 IF (imi(js,im,1).EQ.jcd) imc=im
1049 IF (imi(js,im,1).EQ.jad) ima=im
1051 imi(js,imc,2)=imi(js,ima,1)
1052 imi(js,ima,2)=imi(js,imc,1)
1064 IF (mfound.EQ.0) jcd=jcd+1
1065 ELSEIF (mct(
i,1).EQ.jct.AND.
k(
i,1).GE.1)
THEN
1068 ELSEIF (mct(
i,2).EQ.jct.AND.
k(
i,1).GE.1)
THEN
1072 IF (
i.LE.
n) goto 890
1073 IF (jct.LT.nct) goto 880
1077 IF (iboost.EQ.1)
THEN
1078 DO 900 im=1,
mint(31)
1079 beta=-(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
1080 CALL
pyrobo(imisep(im-1)+1,imisep(im),0d0,0d0,0d0,0d0,
beta)
1089 & ,
'(PYMIHK:) Inconsistent kinematics. Too many boosts.')