pythiaeRHIC/pythia/pyxdif.f

0001
0002 C*********************************************************************
0003
0004 C...PYXDIF
0005 C...Gives the angular orientation of events.
0006
0007       SUBROUTINE PYXDIF(NC,NJET,KFL,ECM,CHI,THE,PHI)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0018
0019 C...Charge. Factors depending on polarization for QED case.
0020       QF=KCHG(KFL,1)/3D0
0021       POLL=1D0-PARJ(131)*PARJ(132)
0022       POLD=PARJ(132)-PARJ(131)
0023       IF(MSTJ(102).LE.1.OR.MSTJ(109).EQ.1) THEN
0024         HF1=POLL
0025         HF2=0D0
0026         HF3=PARJ(133)**2
0027         HF4=0D0
0028
0029 C...Factors depending on flavour, energy and polarization for QFD case.
0030       ELSE
0031         SFF=1D0/(16D0*PARU(102)*(1D0-PARU(102)))
0032         SFW=ECM**4/((ECM**2-PARJ(123)**2)**2+(PARJ(123)*PARJ(124))**2)
0033         SFI=SFW*(1D0-(PARJ(123)/ECM)**2)
0034         AE=-1D0
0035         VE=4D0*PARU(102)-1D0
0036         AF=SIGN(1D0,QF)
0037         VF=AF-4D0*QF*PARU(102)
0038         HF1=QF**2*POLL-2D0*QF*VF*SFI*SFF*(VE*POLL-AE*POLD)+
0039      &  (VF**2+AF**2)*SFW*SFF**2*((VE**2+AE**2)*POLL-2D0*VE*AE*POLD)
0040         HF2=-2D0*QF*AF*SFI*SFF*(AE*POLL-VE*POLD)+2D0*VF*AF*SFW*SFF**2*
0041      &  (2D0*VE*AE*POLL-(VE**2+AE**2)*POLD)
0042         HF3=PARJ(133)**2*(QF**2-2D0*QF*VF*SFI*SFF*VE+(VF**2+AF**2)*
0043      &  SFW*SFF**2*(VE**2-AE**2))
0044         HF4=-PARJ(133)**2*2D0*QF*VF*SFW*(PARJ(123)*PARJ(124)/ECM**2)*
0045      &  SFF*AE
0046       ENDIF
0047
0048 C...Mass factor. Differential cross-sections for two-jet events.
0049       SQ2=SQRT(2D0)
0050       QME=0D0
0051       IF(MSTJ(103).GE.4.AND.IABS(MSTJ(101)).LE.1.AND.MSTJ(102).LE.1.AND.
0052      &MSTJ(109).NE.1) QME=(2D0*PYMASS(KFL)/ECM)**2
0053       IF(NJET.EQ.2) THEN
0054         SIGU=4D0*SQRT(1D0-QME)
0055         SIGL=2D0*QME*SQRT(1D0-QME)
0056         SIGT=0D0
0057         SIGI=0D0
0058         SIGA=0D0
0059         SIGP=4D0
0060
0061 C...Kinematical variables. Reduce four-jet event to three-jet one.
0062       ELSE
0063         IF(NJET.EQ.3) THEN
0064           X1=2D0*P(NC+1,4)/ECM
0065           X2=2D0*P(NC+3,4)/ECM
0066         ELSE
0067           ECMR=P(NC+1,4)+P(NC+4,4)+SQRT((P(NC+2,1)+P(NC+3,1))**2+
0068      &    (P(NC+2,2)+P(NC+3,2))**2+(P(NC+2,3)+P(NC+3,3))**2)
0069           X1=2D0*P(NC+1,4)/ECMR
0070           X2=2D0*P(NC+4,4)/ECMR
0071         ENDIF
0072
0073 C...Differential cross-sections for three-jet (or reduced four-jet).
0074         XQ=(1D0-X1)/(1D0-X2)
0075         CT12=(X1*X2-2D0*X1-2D0*X2+2D0+QME)/SQRT((X1**2-QME)*(X2**2-QME))
0076         ST12=SQRT(1D0-CT12**2)
0077         IF(MSTJ(109).NE.1) THEN
0078           SIGU=2D0*X1**2+X2**2*(1D0+CT12**2)-QME*(3D0+CT12**2-X1-X2)-
0079      &    QME*X1/XQ+0.5D0*QME*((X2**2-QME)*ST12**2-2D0*X2)*XQ
0080           SIGL=(X2*ST12)**2-QME*(3D0-CT12**2-2.5D0*(X1+X2)+X1*X2+QME)+
0081      &    0.5D0*QME*(X1**2-X1-QME)/XQ+0.5D0*QME*((X2**2-QME)*CT12**2-
0082      &    X2)*XQ
0083           SIGT=0.5D0*(X2**2-QME-0.5D0*QME*(X2**2-QME)/XQ)*ST12**2
0084           SIGI=((1D0-0.5D0*QME*XQ)*(X2**2-QME)*ST12*CT12+
0085      &    QME*(1D0-X1-X2+0.5D0*X1*X2+0.5D0*QME)*ST12/CT12)/SQ2
0086           SIGA=X2**2*ST12/SQ2
0087           SIGP=2D0*(X1**2-X2**2*CT12)
0088
0089 C...Differential cross-sect for scalar gluons (no mass effects).
0090         ELSE
0091           X3=2D0-X1-X2
0092           XT=X2*ST12
0093           CT13=SQRT(MAX(0D0,1D0-(XT/X3)**2))
0094           SIGU=(1D0-PARJ(171))*(X3**2-0.5D0*XT**2)+
0095      &    PARJ(171)*(X3**2-0.5D0*XT**2-4D0*(1D0-X1)*(1D0-X2)**2/X1)
0096           SIGL=(1D0-PARJ(171))*0.5D0*XT**2+
0097      &    PARJ(171)*0.5D0*(1D0-X1)**2*XT**2
0098           SIGT=(1D0-PARJ(171))*0.25D0*XT**2+
0099      &    PARJ(171)*0.25D0*XT**2*(1D0-2D0*X1)
0100           SIGI=-(0.5D0/SQ2)*((1D0-PARJ(171))*XT*X3*CT13+
0101      &    PARJ(171)*XT*((1D0-2D0*X1)*X3*CT13-X1*(X1-X2)))
0102           SIGA=(0.25D0/SQ2)*XT*(2D0*(1D0-X1)-X1*X3)
0103           SIGP=X3**2-2D0*(1D0-X1)*(1D0-X2)/X1
0104         ENDIF
0105       ENDIF
0106
0107 C...Upper bounds for differential cross-section.
0108       HF1A=ABS(HF1)
0109       HF2A=ABS(HF2)
0110       HF3A=ABS(HF3)
0111       HF4A=ABS(HF4)
0112       SIGMAX=(2D0*HF1A+HF3A+HF4A)*ABS(SIGU)+2D0*(HF1A+HF3A+HF4A)*
0113      &ABS(SIGL)+2D0*(HF1A+2D0*HF3A+2D0*HF4A)*ABS(SIGT)+2D0*SQ2*
0114      &(HF1A+2D0*HF3A+2D0*HF4A)*ABS(SIGI)+4D0*SQ2*HF2A*ABS(SIGA)+
0115      &2D0*HF2A*ABS(SIGP)
0116
0117 C...Generate angular orientation according to differential cross-sect.
0118   100 CHI=PARU(2)*PYR(0)
0119       CTHE=2D0*PYR(0)-1D0
0120       PHI=PARU(2)*PYR(0)
0121       CCHI=COS(CHI)
0122       SCHI=SIN(CHI)
0123       C2CHI=COS(2D0*CHI)
0124       S2CHI=SIN(2D0*CHI)
0125       THE=ACOS(CTHE)
0126       STHE=SIN(THE)
0127       C2PHI=COS(2D0*(PHI-PARJ(134)))
0128       S2PHI=SIN(2D0*(PHI-PARJ(134)))
0129       SIG=((1D0+CTHE**2)*HF1+STHE**2*(C2PHI*HF3-S2PHI*HF4))*SIGU+
0130      &2D0*(STHE**2*HF1-STHE**2*(C2PHI*HF3-S2PHI*HF4))*SIGL+
0131      &2D0*(STHE**2*C2CHI*HF1+((1D0+CTHE**2)*C2CHI*C2PHI-2D0*CTHE*S2CHI*
0132      &S2PHI)*HF3-((1D0+CTHE**2)*C2CHI*S2PHI+2D0*CTHE*S2CHI*C2PHI)*HF4)*
0133      &SIGT-2D0*SQ2*(2D0*STHE*CTHE*CCHI*HF1-2D0*STHE*(CTHE*CCHI*C2PHI-
0134      &SCHI*S2PHI)*HF3+2D0*STHE*(CTHE*CCHI*S2PHI+SCHI*C2PHI)*HF4)*SIGI+
0135      &4D0*SQ2*STHE*CCHI*HF2*SIGA+2D0*CTHE*HF2*SIGP
0136       IF(SIG.LT.SIGMAX*PYR(0)) GOTO 100
0137
0138       RETURN
0139       END