pythiaeRHIC/pythia/py4frm.f

0001
0002 C*********************************************************************
0003
0004 C...PY4FRM
0005 C...An interface from a four-fermion generator to include
0006 C...parton showers and hadronization.
0007
0008       SUBROUTINE PY4FRM(ATOTSQ,A1SQ,A2SQ,ISTRAT,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       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0018       COMMON/PYINT1/MINT(400),VINT(400)
0019       SAVE /PYJETS/,/PYDAT1/,/PYPARS/,/PYINT1/
0020 C...Local arrays.
0021       DIMENSION IJOIN(2),INTAU(4)
0022
0023 C...Call PYHEPC to convert input from HEPEVT to PYJETS common.
0024       IF(ICOM.EQ.0) THEN
0025         MSTU(28)=0
0026         CALL PYHEPC(2)
0027       ENDIF
0028
0029 C...Loop through entries and pick up all final fermions/antifermions.
0030       I1=0
0031       I2=0
0032       I3=0
0033       I4=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           ELSE
0044             CALL PYERRM(16,'(PY4FRM:) more than two fermions')
0045           ENDIF
0046         ELSE
0047           IF(I2.EQ.0) THEN
0048             I2=I
0049           ELSEIF(I4.EQ.0) THEN
0050             I4=I
0051           ELSE
0052             CALL PYERRM(16,'(PY4FRM:) more than two antifermions')
0053           ENDIF
0054         ENDIF
0055       ENDIF
0056   100 CONTINUE
0057
0058 C...Check that event is arranged according to conventions.
0059       IF(I3.EQ.0.OR.I4.EQ.0) THEN
0060         CALL PYERRM(16,'(PY4FRM:) event contains too few fermions')
0061       ENDIF
0062       IF(I2.LT.I1.OR.I3.LT.I2.OR.I4.LT.I3) THEN
0063         CALL PYERRM(6,'(PY4FRM:) fermions arranged in wrong order')
0064       ENDIF
0065
0066 C...Check which fermion pairs are quarks and which leptons.
0067       IF(IABS(K(I1,2)).LT.10.AND.IABS(K(I2,2)).LT.10) THEN
0068         IQL12=1
0069       ELSEIF(IABS(K(I1,2)).GT.10.AND.IABS(K(I2,2)).GT.10) THEN
0070         IQL12=2
0071       ELSE
0072         CALL PYERRM(16,'(PY4FRM:) first fermion pair inconsistent')
0073       ENDIF
0074       IF(IABS(K(I3,2)).LT.10.AND.IABS(K(I4,2)).LT.10) THEN
0075         IQL34=1
0076       ELSEIF(IABS(K(I3,2)).GT.10.AND.IABS(K(I4,2)).GT.10) THEN
0077         IQL34=2
0078       ELSE
0079         CALL PYERRM(16,'(PY4FRM:) second fermion pair inconsistent')
0080       ENDIF
0081
0082 C...Decide whether to allow or not photon radiation in showers.
0083       MSTJ(41)=2
0084       IF(IRAD.EQ.0) MSTJ(41)=1
0085
0086 C...Decide on dipole pairing.
0087       IP1=I1
0088       IP2=I2
0089       IP3=I3
0090       IP4=I4
0091       IF(IQL12.EQ.IQL34) THEN
0092         R1SQ=A1SQ
0093         R2SQ=A2SQ
0094         DELTA=ATOTSQ-A1SQ-A2SQ
0095         IF(ISTRAT.EQ.1) THEN
0096           IF(DELTA.GT.0D0) R1SQ=R1SQ+DELTA
0097           IF(DELTA.LT.0D0) R2SQ=MAX(0D0,R2SQ+DELTA)
0098         ELSEIF(ISTRAT.EQ.2) THEN
0099           IF(DELTA.GT.0D0) R2SQ=R2SQ+DELTA
0100           IF(DELTA.LT.0D0) R1SQ=MAX(0D0,R1SQ+DELTA)
0101         ENDIF
0102         IF(R2SQ.GT.PYR(0)*(R1SQ+R2SQ)) THEN
0103           IP2=I4
0104           IP4=I2
0105         ENDIF
0106       ENDIF
0107
0108 C...If colour reconnection then bookkeep W+W- or Z0Z0
0109 C...and copy q qbar q qbar consecutively.
0110       IF(MSTP(115).GE.1.AND.IQL12.EQ.1.AND.IQL34.EQ.1) THEN
0111         K(N+1,1)=11
0112         K(N+1,3)=IP1
0113         K(N+1,4)=N+3
0114         K(N+1,5)=N+4
0115         K(N+2,1)=11
0116         K(N+2,3)=IP3
0117         K(N+2,4)=N+5
0118         K(N+2,5)=N+6
0119         IF(K(IP1,2)+K(IP2,2).EQ.0) THEN
0120           K(N+1,2)=23
0121           K(N+2,2)=23
0122           MINT(1)=22
0123         ELSEIF(PYCHGE(K(IP1,2)).GT.0) THEN
0124           K(N+1,2)=24
0125           K(N+2,2)=-24
0126           MINT(1)=25
0127         ELSE
0128           K(N+1,2)=-24
0129           K(N+2,2)=24
0130           MINT(1)=25
0131         ENDIF
0132         DO 110 J=1,5
0133           K(N+3,J)=K(IP1,J)
0134           K(N+4,J)=K(IP2,J)
0135           K(N+5,J)=K(IP3,J)
0136           K(N+6,J)=K(IP4,J)
0137           P(N+1,J)=P(IP1,J)+P(IP2,J)
0138           P(N+2,J)=P(IP3,J)+P(IP4,J)
0139           P(N+3,J)=P(IP1,J)
0140           P(N+4,J)=P(IP2,J)
0141           P(N+5,J)=P(IP3,J)
0142           P(N+6,J)=P(IP4,J)
0143           V(N+1,J)=V(IP1,J)
0144           V(N+2,J)=V(IP3,J)
0145           V(N+3,J)=V(IP1,J)
0146           V(N+4,J)=V(IP2,J)
0147           V(N+5,J)=V(IP3,J)
0148           V(N+6,J)=V(IP4,J)
0149   110   CONTINUE
0150         P(N+1,5)=SQRT(MAX(0D0,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2-
0151      &  P(N+1,3)**2))
0152         P(N+2,5)=SQRT(MAX(0D0,P(N+2,4)**2-P(N+2,1)**2-P(N+2,2)**2-
0153      &  P(N+2,3)**2))
0154         K(N+3,3)=N+1
0155         K(N+4,3)=N+1
0156         K(N+5,3)=N+2
0157         K(N+6,3)=N+2
0158 C...Remove original q qbar q qbar and update counters.
0159         K(IP1,1)=K(IP1,1)+10
0160         K(IP2,1)=K(IP2,1)+10
0161         K(IP3,1)=K(IP3,1)+10
0162         K(IP4,1)=K(IP4,1)+10
0163         IW1=N+1
0164         IW2=N+2
0165         NSD1=N+2
0166         IP1=N+3
0167         IP2=N+4
0168         IP3=N+5
0169         IP4=N+6
0170         N=N+6
0171       ENDIF
0172
0173 C...Do colour joinings and parton showers.
0174       IF(IQL12.EQ.1) THEN
0175         IJOIN(1)=IP1
0176         IJOIN(2)=IP2
0177         CALL PYJOIN(2,IJOIN)
0178       ENDIF
0179       IF(IQL12.EQ.1.OR.IRAD.EQ.1) THEN
0180         PM12S=(P(IP1,4)+P(IP2,4))**2-(P(IP1,1)+P(IP2,1))**2-
0181      &  (P(IP1,2)+P(IP2,2))**2-(P(IP1,3)+P(IP2,3))**2
0182         CALL PYSHOW(IP1,IP2,SQRT(MAX(0D0,PM12S)))
0183       ENDIF
0184       NAFT1=N
0185       IF(IQL34.EQ.1) THEN
0186         IJOIN(1)=IP3
0187         IJOIN(2)=IP4
0188         CALL PYJOIN(2,IJOIN)
0189       ENDIF
0190       IF(IQL34.EQ.1.OR.IRAD.EQ.1) THEN
0191         PM34S=(P(IP3,4)+P(IP4,4))**2-(P(IP3,1)+P(IP4,1))**2-
0192      &  (P(IP3,2)+P(IP4,2))**2-(P(IP3,3)+P(IP4,3))**2
0193         CALL PYSHOW(IP3,IP4,SQRT(MAX(0D0,PM34S)))
0194       ENDIF
0195
0196 C...Optionally do colour reconnection.
0197       MINT(32)=0
0198       MSTI(32)=0
0199       IF(MSTP(115).GE.1.AND.IQL12.EQ.1.AND.IQL34.EQ.1) THEN
0200         CALL PYRECO(IW1,IW2,NSD1,NAFT1)
0201         MSTI(32)=MINT(32)
0202       ENDIF
0203
0204 C...Do fragmentation and decays. Possibly except tau decay.
0205       IF(ITAU.EQ.0) THEN
0206         NTAU=0
0207         DO 120 I=1,N
0208         IF(IABS(K(I,2)).EQ.15.AND.K(I,1).EQ.1) THEN
0209           NTAU=NTAU+1
0210           INTAU(NTAU)=I
0211           K(I,1)=11
0212         ENDIF
0213   120   CONTINUE
0214       ENDIF
0215       CALL PYEXEC
0216       IF(ITAU.EQ.0) THEN
0217         DO 130 I=1,NTAU
0218         K(INTAU(I),1)=1
0219   130   CONTINUE
0220       ENDIF
0221
0222 C...Call PYHEPC to convert output from PYJETS to HEPEVT common.
0223       IF(ICOM.EQ.0) THEN
0224         MSTU(28)=0
0225         CALL PYHEPC(1)
0226       ENDIF
0227
0228       END