pythiaeRHIC/pythia/pyonia.f

0001
0002 C*********************************************************************
0003
0004 C...PYONIA
0005 C...Generates Upsilon and toponium decays into three gluons
0006 C...or two gluons and a photon.
0007
0008       SUBROUTINE PYONIA(KFL,ECM)
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/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0019
0020 C...Printout. Check input parameters.
0021       IF(MSTU(12).NE.12345) CALL PYLIST(0)
0022       IF(KFL.LT.0.OR.KFL.GT.8) THEN
0023         CALL PYERRM(16,'(PYONIA:) called with unknown flavour code')
0024         IF(MSTU(21).GE.1) RETURN
0025       ENDIF
0026       IF(ECM.LT.PARJ(127)+2.02D0*PARF(101)) THEN
0027         CALL PYERRM(16,'(PYONIA:) called with too small CM energy')
0028         IF(MSTU(21).GE.1) RETURN
0029       ENDIF
0030
0031 C...Initial e+e- and onium state (optional).
0032       NC=0
0033       IF(MSTJ(115).GE.2) THEN
0034         NC=NC+2
0035         CALL PY1ENT(NC-1,11,0.5D0*ECM,0D0,0D0)
0036         K(NC-1,1)=21
0037         CALL PY1ENT(NC,-11,0.5D0*ECM,PARU(1),0D0)
0038         K(NC,1)=21
0039       ENDIF
0040       KFLC=IABS(KFL)
0041       IF(MSTJ(115).GE.3.AND.KFLC.GE.5) THEN
0042         NC=NC+1
0043         KF=110*KFLC+3
0044         MSTU10=MSTU(10)
0045         MSTU(10)=1
0046         P(NC,5)=ECM
0047         CALL PY1ENT(NC,KF,ECM,0D0,0D0)
0048         K(NC,1)=21
0049         K(NC,3)=1
0050         MSTU(10)=MSTU10
0051       ENDIF
0052
0053 C...Choose x1 and x2 according to matrix element.
0054       NTRY=0
0055   100 X1=PYR(0)
0056       X2=PYR(0)
0057       X3=2D0-X1-X2
0058       IF(X3.GE.1D0.OR.((1D0-X1)/(X2*X3))**2+((1D0-X2)/(X1*X3))**2+
0059      &((1D0-X3)/(X1*X2))**2.LE.2D0*PYR(0)) GOTO 100
0060       NTRY=NTRY+1
0061       NJET=3
0062       IF(MSTJ(101).LE.4) CALL PY3ENT(NC+1,21,21,21,ECM,X1,X3)
0063       IF(MSTJ(101).GE.5) CALL PY3ENT(-(NC+1),21,21,21,ECM,X1,X3)
0064
0065 C...Photon-gluon-gluon events. Small system modifications. Jet origin.
0066       MSTU(111)=MSTJ(108)
0067       IF(MSTJ(108).EQ.2.AND.(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.1))
0068      &MSTU(111)=1
0069       PARU(112)=PARJ(121)
0070       IF(MSTU(111).EQ.2) PARU(112)=PARJ(122)
0071       QF=0D0
0072       IF(KFLC.NE.0) QF=KCHG(KFLC,1)/3D0
0073       RGAM=7.2D0*QF**2*PARU(101)/PYALPS(ECM**2)
0074       MK=0
0075       ECMC=ECM
0076       IF(PYR(0).GT.RGAM/(1D0+RGAM)) THEN
0077         IF(1D0-MAX(X1,X2,X3).LE.MAX((PARJ(126)/ECM)**2,PARJ(125)))
0078      &  NJET=2
0079         IF(NJET.EQ.2.AND.MSTJ(101).LE.4) CALL PY2ENT(NC+1,21,21,ECM)
0080         IF(NJET.EQ.2.AND.MSTJ(101).GE.5) CALL PY2ENT(-(NC+1),21,21,ECM)
0081       ELSE
0082         MK=1
0083         ECMC=SQRT(1D0-X1)*ECM
0084         IF(ECMC.LT.2D0*PARJ(127)) GOTO 100
0085         K(NC+1,1)=1
0086         K(NC+1,2)=22
0087         K(NC+1,4)=0
0088         K(NC+1,5)=0
0089         IF(MSTJ(101).GE.5) K(NC+2,4)=MSTU(5)*(NC+3)
0090         IF(MSTJ(101).GE.5) K(NC+2,5)=MSTU(5)*(NC+3)
0091         IF(MSTJ(101).GE.5) K(NC+3,4)=MSTU(5)*(NC+2)
0092         IF(MSTJ(101).GE.5) K(NC+3,5)=MSTU(5)*(NC+2)
0093         NJET=2
0094         IF(ECMC.LT.4D0*PARJ(127)) THEN
0095           MSTU10=MSTU(10)
0096           MSTU(10)=1
0097           P(NC+2,5)=ECMC
0098           CALL PY1ENT(NC+2,83,0.5D0*(X2+X3)*ECM,PARU(1),0D0)
0099           MSTU(10)=MSTU10
0100           NJET=0
0101         ENDIF
0102       ENDIF
0103       DO 110 IP=NC+1,N
0104         K(IP,3)=K(IP,3)+(MSTJ(115)/2)+(KFLC/5)*(MSTJ(115)/3)*(NC-1)
0105   110 CONTINUE
0106
0107 C...Differential cross-sections. Upper limit for cross-section.
0108       IF(MSTJ(106).EQ.1) THEN
0109         SQ2=SQRT(2D0)
0110         HF1=1D0-PARJ(131)*PARJ(132)
0111         HF3=PARJ(133)**2
0112         CT13=(X1*X3-2D0*X1-2D0*X3+2D0)/(X1*X3)
0113         ST13=SQRT(1D0-CT13**2)
0114         SIGL=0.5D0*X3**2*((1D0-X2)**2+(1D0-X3)**2)*ST13**2
0115         SIGU=(X1*(1D0-X1))**2+(X2*(1D0-X2))**2+(X3*(1D0-X3))**2-SIGL
0116         SIGT=0.5D0*SIGL
0117         SIGI=(SIGL*CT13/ST13+0.5D0*X1*X3*(1D0-X2)**2*ST13)/SQ2
0118         SIGMAX=(2D0*HF1+HF3)*ABS(SIGU)+2D0*(HF1+HF3)*ABS(SIGL)+2D0*(HF1+
0119      &  2D0*HF3)*ABS(SIGT)+2D0*SQ2*(HF1+2D0*HF3)*ABS(SIGI)
0120
0121 C...Angular orientation of event.
0122   120   CHI=PARU(2)*PYR(0)
0123         CTHE=2D0*PYR(0)-1D0
0124         PHI=PARU(2)*PYR(0)
0125         CCHI=COS(CHI)
0126         SCHI=SIN(CHI)
0127         C2CHI=COS(2D0*CHI)
0128         S2CHI=SIN(2D0*CHI)
0129         THE=ACOS(CTHE)
0130         STHE=SIN(THE)
0131         C2PHI=COS(2D0*(PHI-PARJ(134)))
0132         S2PHI=SIN(2D0*(PHI-PARJ(134)))
0133         SIG=((1D0+CTHE**2)*HF1+STHE**2*C2PHI*HF3)*SIGU+2D0*(STHE**2*HF1-
0134      &  STHE**2*C2PHI*HF3)*SIGL+2D0*(STHE**2*C2CHI*HF1+((1D0+CTHE**2)*
0135      &  C2CHI*C2PHI-2D0*CTHE*S2CHI*S2PHI)*HF3)*SIGT-
0136      &  2D0*SQ2*(2D0*STHE*CTHE*CCHI*HF1-2D0*STHE*
0137      &  (CTHE*CCHI*C2PHI-SCHI*S2PHI)*HF3)*SIGI
0138         IF(SIG.LT.SIGMAX*PYR(0)) GOTO 120
0139         CALL PYROBO(NC+1,N,0D0,CHI,0D0,0D0,0D0)
0140         CALL PYROBO(NC+1,N,THE,PHI,0D0,0D0,0D0)
0141       ENDIF
0142
0143 C...Generate parton shower. Rearrange along strings and check.
0144       IF(MSTJ(101).GE.5.AND.NJET.GE.2) THEN
0145         CALL PYSHOW(NC+MK+1,-NJET,ECMC)
0146         MSTJ14=MSTJ(14)
0147         IF(MSTJ(105).EQ.-1) MSTJ(14)=-1
0148         IF(MSTJ(105).GE.0) MSTU(28)=0
0149         CALL PYPREP(0)
0150         MSTJ(14)=MSTJ14
0151         IF(MSTJ(105).GE.0.AND.MSTU(28).NE.0) GOTO 100
0152       ENDIF
0153
0154 C...Generate fragmentation. Information for PYTABU:
0155       IF(MSTJ(105).EQ.1) CALL PYEXEC
0156       MSTU(161)=110*KFLC+3
0157       MSTU(162)=0
0158
0159       RETURN
0160       END