pythiaeRHIC/pythia/py6frm.f

0001
0002 C*********************************************************************
0003
0004 C...PY6FRM
0005 C...An interface from a six-fermion generator to include
0006 C...parton showers and hadronization.
0007
0008       SUBROUTINE PY6FRM(P12,P13,P21,P23,P31,P32,PTOP,IRAD,ITAU,ICOM)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       SAVE /PYJETS/,/PYDAT1/
0018 C...Local arrays.
0019       DIMENSION IJOIN(2),INTAU(6),BETA(3),BETAO(3),BETAN(3)
0020
0021 C...Call PYHEPC to convert input from HEPEVT to PYJETS common.
0022       IF(ICOM.EQ.0) THEN
0023         MSTU(28)=0
0024         CALL PYHEPC(2)
0025       ENDIF
0026
0027 C...Loop through entries and pick up all final fermions/antifermions.
0028       I1=0
0029       I2=0
0030       I3=0
0031       I4=0
0032       I5=0
0033       I6=0
0034       DO 100 I=1,N
0035       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 100
0036       KFA=IABS(K(I,2))
0037       IF((KFA.GE.1.AND.KFA.LE.6).OR.(KFA.GE.11.AND.KFA.LE.16)) THEN
0038         IF(K(I,2).GT.0) THEN
0039           IF(I1.EQ.0) THEN
0040             I1=I
0041           ELSEIF(I3.EQ.0) THEN
0042             I3=I
0043           ELSEIF(I5.EQ.0) THEN
0044             I5=I
0045           ELSE
0046             CALL PYERRM(16,'(PY6FRM:) more than three fermions')
0047           ENDIF
0048         ELSE
0049           IF(I2.EQ.0) THEN
0050             I2=I
0051           ELSEIF(I4.EQ.0) THEN
0052             I4=I
0053           ELSEIF(I6.EQ.0) THEN
0054             I6=I
0055           ELSE
0056             CALL PYERRM(16,'(PY6FRM:) more than three antifermions')
0057           ENDIF
0058         ENDIF
0059       ENDIF
0060   100 CONTINUE
0061
0062 C...Check that event is arranged according to conventions.
0063       IF(I5.EQ.0.OR.I6.EQ.0) THEN
0064         CALL PYERRM(16,'(PY6FRM:) event contains too few fermions')
0065       ENDIF
0066       IF(I2.LT.I1.OR.I3.LT.I2.OR.I4.LT.I3.OR.I5.LT.I4.OR.I6.LT.I5) THEN
0067         CALL PYERRM(6,'(PY6FRM:) fermions arranged in wrong order')
0068       ENDIF
0069
0070 C...Check which fermion pairs are quarks and which leptons.
0071       IF(IABS(K(I1,2)).LT.10.AND.IABS(K(I2,2)).LT.10) THEN
0072         IQL12=1
0073       ELSEIF(IABS(K(I1,2)).GT.10.AND.IABS(K(I2,2)).GT.10) THEN
0074         IQL12=2
0075       ELSE
0076         CALL PYERRM(16,'(PY6FRM:) first fermion pair inconsistent')
0077       ENDIF
0078       IF(IABS(K(I3,2)).LT.10.AND.IABS(K(I4,2)).LT.10) THEN
0079         IQL34=1
0080       ELSEIF(IABS(K(I3,2)).GT.10.AND.IABS(K(I4,2)).GT.10) THEN
0081         IQL34=2
0082       ELSE
0083         CALL PYERRM(16,'(PY6FRM:) second fermion pair inconsistent')
0084       ENDIF
0085       IF(IABS(K(I5,2)).LT.10.AND.IABS(K(I6,2)).LT.10) THEN
0086         IQL56=1
0087       ELSEIF(IABS(K(I5,2)).GT.10.AND.IABS(K(I6,2)).GT.10) THEN
0088         IQL56=2
0089       ELSE
0090         CALL PYERRM(16,'(PY6FRM:) third fermion pair inconsistent')
0091       ENDIF
0092
0093 C...Decide whether to allow or not photon radiation in showers.
0094       MSTJ(41)=2
0095       IF(IRAD.EQ.0) MSTJ(41)=1
0096
0097 C...Allow dipole pairings only among leptons and quarks separately.
0098       P12D=P12
0099       P13D=0D0
0100       IF(IQL34.EQ.IQL56) P13D=P13
0101       P21D=0D0
0102       IF(IQL12.EQ.IQL34) P21D=P21
0103       P23D=0D0
0104       IF(IQL12.EQ.IQL34.AND.IQL12.EQ.IQL56) P23D=P23
0105       P31D=0D0
0106       IF(IQL12.EQ.IQL34.AND.IQL12.EQ.IQL56) P31D=P31
0107       P32D=0D0
0108       IF(IQL12.EQ.IQL56) P32D=P32
0109
0110 C...Decide whether t+tbar.
0111       ITOP=0
0112       IF(PYR(0).LT.PTOP) THEN
0113         ITOP=1
0114
0115 C...If t+tbar: reconstruct t's.
0116         IT=N+1
0117         ITB=N+2
0118         DO 110 J=1,5
0119           K(IT,J)=0
0120           K(ITB,J)=0
0121           P(IT,J)=P(I1,J)+P(I3,J)+P(I4,J)
0122           P(ITB,J)=P(I2,J)+P(I5,J)+P(I6,J)
0123           V(IT,J)=0D0
0124           V(ITB,J)=0D0
0125   110   CONTINUE
0126         K(IT,1)=1
0127         K(ITB,1)=1
0128         K(IT,2)=6
0129         K(ITB,2)=-6
0130         P(IT,5)=SQRT(MAX(0D0,P(IT,4)**2-P(IT,1)**2-P(IT,2)**2-
0131      &  P(IT,3)**2))
0132         P(ITB,5)=SQRT(MAX(0D0,P(ITB,4)**2-P(ITB,1)**2-P(ITB,2)**2-
0133      &  P(ITB,3)**2))
0134         N=N+2
0135
0136 C...If t+tbar: colour join t's and let them shower.
0137         IJOIN(1)=IT
0138         IJOIN(2)=ITB
0139         CALL PYJOIN(2,IJOIN)
0140         PMTTS=(P(IT,4)+P(ITB,4))**2-(P(IT,1)+P(ITB,1))**2-
0141      &  (P(IT,2)+P(ITB,2))**2-(P(IT,3)+P(ITB,3))**2
0142         CALL PYSHOW(IT,ITB,SQRT(MAX(0D0,PMTTS)))
0143
0144 C...If t+tbar: pick up the t's after shower.
0145         ITNEW=IT
0146         ITBNEW=ITB
0147         DO 120 I=ITB+1,N
0148           IF(K(I,2).EQ.6) ITNEW=I
0149           IF(K(I,2).EQ.-6) ITBNEW=I
0150   120   CONTINUE
0151
0152 C...If t+tbar: loop over two top systems.
0153         DO 200 IT1=1,2
0154           IF(IT1.EQ.1) THEN
0155             ITO=IT
0156             ITN=ITNEW
0157             IBO=I1
0158             IW1=I3
0159             IW2=I4
0160           ELSE
0161             ITO=ITB
0162             ITN=ITBNEW
0163             IBO=I2
0164             IW1=I5
0165             IW2=I6
0166           ENDIF
0167           IF(IABS(K(IBO,2)).NE.5) CALL PYERRM(6,
0168      &    '(PY6FRM:) not b in t decay')
0169
0170 C...If t+tbar: find boost from original to new top frame.
0171           DO 130 J=1,3
0172             BETAO(J)=P(ITO,J)/P(ITO,4)
0173             BETAN(J)=P(ITN,J)/P(ITN,4)
0174   130     CONTINUE
0175
0176 C...If t+tbar: boost copy of b by t shower and connect it in colour.
0177           N=N+1
0178           IB=N
0179           K(IB,1)=3
0180           K(IB,2)=K(IBO,2)
0181           K(IB,3)=ITN
0182           DO 140 J=1,5
0183             P(IB,J)=P(IBO,J)
0184             V(IB,J)=0D0
0185   140     CONTINUE
0186           CALL PYROBO(IB,IB,0D0,0D0,-BETAO(1),-BETAO(2),-BETAO(3))
0187           CALL PYROBO(IB,IB,0D0,0D0,BETAN(1),BETAN(2),BETAN(3))
0188           K(IB,4)=MSTU(5)*ITN
0189           K(IB,5)=MSTU(5)*ITN
0190           K(ITN,4)=K(ITN,4)+IB
0191           K(ITN,5)=K(ITN,5)+IB
0192           K(ITN,1)=K(ITN,1)+10
0193           K(IBO,1)=K(IBO,1)+10
0194
0195 C...If t+tbar: construct W recoiling against b.
0196           N=N+1
0197           IW=N
0198           DO 150 J=1,5
0199             K(IW,J)=0
0200             V(IW,J)=0D0
0201   150     CONTINUE
0202           K(IW,1)=1
0203           KCHW=PYCHGE(K(IW1,2))+PYCHGE(K(IW2,2))
0204           IF(IABS(KCHW).EQ.3) THEN
0205             K(IW,2)=ISIGN(24,KCHW)
0206           ELSE
0207             CALL PYERRM(16,'(PY6FRM:) fermion pair inconsistent with W')
0208           ENDIF
0209           K(IW,3)=IW1
0210
0211 C...If t+tbar: construct W momentum, including boost by t shower.
0212           DO 160 J=1,4
0213             P(IW,J)=P(IW1,J)+P(IW2,J)
0214   160     CONTINUE
0215           P(IW,5)=SQRT(MAX(0D0,P(IW,4)**2-P(IW,1)**2-P(IW,2)**2-
0216      &    P(IW,3)**2))
0217           CALL PYROBO(IW,IW,0D0,0D0,-BETAO(1),-BETAO(2),-BETAO(3))
0218           CALL PYROBO(IW,IW,0D0,0D0,BETAN(1),BETAN(2),BETAN(3))
0219
0220 C...If t+tbar: boost b and W to top rest frame.
0221           DO 170 J=1,3
0222             BETA(J)=(P(IB,J)+P(IW,J))/(P(IB,4)+P(IW,4))
0223   170     CONTINUE
0224           CALL PYROBO(IB,IB,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0225           CALL PYROBO(IW,IW,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0226
0227 C...If t+tbar: let b shower and pick up modified W.
0228           PMTS=(P(IB,4)+P(IW,4))**2-(P(IB,1)+P(IW,1))**2-
0229      &    (P(IB,2)+P(IW,2))**2-(P(IB,3)+P(IW,3))**2
0230           CALL PYSHOW(IB,IW,SQRT(MAX(0D0,PMTS)))
0231           DO 180 I=IW,N
0232             IF(IABS(K(I,2)).EQ.24) IWM=I
0233   180     CONTINUE
0234
0235 C...If t+tbar: take copy of W decay products.
0236           DO 190 J=1,5
0237             K(N+1,J)=K(IW1,J)
0238             P(N+1,J)=P(IW1,J)
0239             V(N+1,J)=V(IW1,J)
0240             K(N+2,J)=K(IW2,J)
0241             P(N+2,J)=P(IW2,J)
0242             V(N+2,J)=V(IW2,J)
0243   190     CONTINUE
0244           K(IW1,1)=K(IW1,1)+10
0245           K(IW2,1)=K(IW2,1)+10
0246           K(IWM,1)=K(IWM,1)+10
0247           K(IWM,4)=N+1
0248           K(IWM,5)=N+2
0249           K(N+1,3)=IWM
0250           K(N+2,3)=IWM
0251           IF(IT1.EQ.1) THEN
0252             I3=N+1
0253             I4=N+2
0254           ELSE
0255             I5=N+1
0256             I6=N+2
0257           ENDIF
0258           N=N+2
0259
0260 C...If t+tbar: boost W decay products, first by effects of t shower,
0261 C...then by those of b shower. b and its shower simple boost back.
0262           CALL PYROBO(N-1,N,0D0,0D0,-BETAO(1),-BETAO(2),-BETAO(3))
0263           CALL PYROBO(N-1,N,0D0,0D0,BETAN(1),BETAN(2),BETAN(3))
0264           CALL PYROBO(N-1,N,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0265           CALL PYROBO(N-1,N,0D0,0D0,-P(IW,1)/P(IW,4),
0266      &    -P(IW,2)/P(IW,4),-P(IW,3)/P(IW,4))
0267           CALL PYROBO(N-1,N,0D0,0D0,P(IWM,1)/P(IWM,4),
0268      &    P(IWM,2)/P(IWM,4),P(IWM,3)/P(IWM,4))
0269           CALL PYROBO(IB,IB,0D0,0D0,BETA(1),BETA(2),BETA(3))
0270           CALL PYROBO(IW,N,0D0,0D0,BETA(1),BETA(2),BETA(3))
0271   200   CONTINUE
0272       ENDIF
0273
0274 C...Decide on dipole pairing.
0275       IP1=I1
0276       IP3=I3
0277       IP5=I5
0278       PRN=PYR(0)*(P12D+P13D+P21D+P23D+P31D+P32D)
0279       IF(ITOP.EQ.1.OR.PRN.LT.P12D) THEN
0280         IP2=I2
0281         IP4=I4
0282         IP6=I6
0283       ELSEIF(PRN.LT.P12D+P13D) THEN
0284         IP2=I2
0285         IP4=I6
0286         IP6=I4
0287       ELSEIF(PRN.LT.P12D+P13D+P21D) THEN
0288         IP2=I4
0289         IP4=I2
0290         IP6=I6
0291       ELSEIF(PRN.LT.P12D+P13D+P21D+P23D) THEN
0292         IP2=I4
0293         IP4=I6
0294         IP6=I2
0295       ELSEIF(PRN.LT.P12D+P13D+P21D+P23D+P31D) THEN
0296         IP2=I6
0297         IP4=I2
0298         IP6=I4
0299       ELSE
0300         IP2=I6
0301         IP4=I4
0302         IP6=I2
0303       ENDIF
0304
0305 C...Do colour joinings and parton showers
0306 C...(except ones already made for t+tbar).
0307       IF(ITOP.EQ.0) THEN
0308         IF(IQL12.EQ.1) THEN
0309           IJOIN(1)=IP1
0310           IJOIN(2)=IP2
0311           CALL PYJOIN(2,IJOIN)
0312         ENDIF
0313         IF(IQL12.EQ.1.OR.IRAD.EQ.1) THEN
0314           PM12S=(P(IP1,4)+P(IP2,4))**2-(P(IP1,1)+P(IP2,1))**2-
0315      &    (P(IP1,2)+P(IP2,2))**2-(P(IP1,3)+P(IP2,3))**2
0316           CALL PYSHOW(IP1,IP2,SQRT(MAX(0D0,PM12S)))
0317         ENDIF
0318       ENDIF
0319       IF(IQL34.EQ.1) THEN
0320         IJOIN(1)=IP3
0321         IJOIN(2)=IP4
0322         CALL PYJOIN(2,IJOIN)
0323       ENDIF
0324       IF(IQL34.EQ.1.OR.IRAD.EQ.1) THEN
0325         PM34S=(P(IP3,4)+P(IP4,4))**2-(P(IP3,1)+P(IP4,1))**2-
0326      &  (P(IP3,2)+P(IP4,2))**2-(P(IP3,3)+P(IP4,3))**2
0327         CALL PYSHOW(IP3,IP4,SQRT(MAX(0D0,PM34S)))
0328       ENDIF
0329       IF(IQL56.EQ.1) THEN
0330         IJOIN(1)=IP5
0331         IJOIN(2)=IP6
0332         CALL PYJOIN(2,IJOIN)
0333       ENDIF
0334       IF(IQL56.EQ.1.OR.IRAD.EQ.1) THEN
0335         PM56S=(P(IP5,4)+P(IP6,4))**2-(P(IP5,1)+P(IP6,1))**2-
0336      &  (P(IP5,2)+P(IP6,2))**2-(P(IP5,3)+P(IP6,3))**2
0337         CALL PYSHOW(IP5,IP6,SQRT(MAX(0D0,PM56S)))
0338       ENDIF
0339
0340 C...Do fragmentation and decays. Possibly except tau decay.
0341       IF(ITAU.EQ.0) THEN
0342         NTAU=0
0343         DO 210 I=1,N
0344         IF(IABS(K(I,2)).EQ.15.AND.K(I,1).EQ.1) THEN
0345           NTAU=NTAU+1
0346           INTAU(NTAU)=I
0347           K(I,1)=11
0348         ENDIF
0349   210   CONTINUE
0350       ENDIF
0351       CALL PYEXEC
0352       IF(ITAU.EQ.0) THEN
0353         DO 220 I=1,NTAU
0354         K(INTAU(I),1)=1
0355   220   CONTINUE
0356       ENDIF
0357
0358 C...Call PYHEPC to convert output from PYJETS to HEPEVT common.
0359       IF(ICOM.EQ.0) THEN
0360         MSTU(28)=0
0361         CALL PYHEPC(1)
0362       ENDIF
0363
0364       END