pythiaeRHIC/pythia/pytbdy.f

0001
0002 C*********************************************************************
0003
0004 C...PYTBDY
0005 C...Generates 3-body decays of gauginos.
0006
0007       SUBROUTINE PYTBDY(IDIN)
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...Parameter statement to help give large particle numbers.
0014       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0015      &KEXCIT=4000000,KDIMEN=5000000)
0016 C...Commonblocks.
0017       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0018       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0019       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0020 C     COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0021 C     COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0022       COMMON/PYSSMT/ZMIX(4,4),UMIX(2,2),VMIX(2,2),SMZ(4),SMW(2),
0023      &SFMIX(16,4),ZMIXI(4,4),UMIXI(2,2),VMIXI(2,2)
0024 C     SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYSSMT/
0025       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSSMT/
0026
0027 C...Local variables.
0028       DOUBLE PRECISION XM(5)
0029       COMPLEX*16 OLPP,ORPP,QLL,QLR,QRR,QRL,GLIJ,GRIJ,PROPZ
0030       COMPLEX*16 QLLS,QRRS,QLRS,QRLS,QLLU,QRRU,QLRT,QRLT
0031       COMPLEX*16 ZMIXC(4,4),UMIXC(2,2),VMIXC(2,2)
0032       DOUBLE PRECISION S12MIN,S12MAX,YJACO1,S23AVE,S23DF1,S23DF2
0033       DOUBLE PRECISION D1,D2,D3,P1,P2,P3,CTHE1,STHE1,CTHE3,STHE3
0034       DOUBLE PRECISION CPHI1,SPHI1
0035       DOUBLE PRECISION S23DEL,EPS
0036       DOUBLE PRECISION GOLDEN,AX,BX,CX,TOL,XMIN,R,C
0037       PARAMETER (R=0.61803399D0,C=1D0-R,TOL=1D-3)
0038       DOUBLE PRECISION F1,F2,X0,X1,X2,X3
0039       INTEGER INOID(4)
0040       DATA INOID/22,23,25,35/
0041       DATA EPS/1D-6/
0042
0043       ID=IDIN
0044       ISKIP=1
0045       XM(1)=P(N+1,5)
0046       XM(2)=P(N+2,5)
0047       XM(3)=P(N+3,5)
0048       XM(5)=P(ID,5)
0049
0050 C...GENERATE S12
0051       S12MIN=(XM(1)+XM(2))**2
0052       S12MAX=(XM(5)-XM(3))**2
0053       YJACO1=S12MAX-S12MIN
0054
0055 C...Initialize some parameters
0056       XW=PARU(102)
0057       XW1=1D0-XW
0058       TANW=SQRT(XW/XW1)
0059       IZID1=0
0060       IWID1=0
0061       IZID2=0
0062       IWID2=0
0063
0064       IA=K(N+2,2)
0065       JA=K(N+3,2)
0066
0067 C...Mrenna: check that we are indeed decaying a SUSY particle
0068       IF(IABS(K(ID,2)).LT.KSUSY1.OR.IABS(K(ID,2)).GE.3000000) THEN
0069
0070       ELSE
0071         DO 100 I1=1,4
0072           IF(MOD(K(N+1,2),KSUSY1).EQ.INOID(I1)) IZID1=I1
0073           IF(MOD(K(ID,2),KSUSY1).EQ.INOID(I1)) IZID2=I1
0074  100    CONTINUE
0075         IF(MOD(K(N+1,2),KSUSY1).EQ.24) IWID1=1
0076         IF(MOD(K(N+1,2),KSUSY1).EQ.37) IWID1=2
0077         IF(MOD(K(ID,2),KSUSY1).EQ.24) IWID2=1
0078         IF(MOD(K(ID,2),KSUSY1).EQ.37) IWID2=2
0079         ZM12=XM(5)**2
0080         ZM22=XM(1)**2
0081         EI=KCHG(PYCOMP(IABS(IA)),1)/3D0
0082         T3I=SIGN(1D0,EI+1D-6)/2D0
0083       ENDIF
0084
0085       IF(MAX(ABS(IA),ABS(JA)).EQ.6) THEN
0086         ISKIP=0
0087       ELSEIF(IZID1*IZID2.NE.0) THEN
0088         SQMZ=PMAS(23,1)**2
0089         GMMZ=PMAS(23,1)*PMAS(23,2)
0090         DO 110 I=1,4
0091           ZMIXC(IZID1,I)=DCMPLX(ZMIX(IZID1,I),ZMIXI(IZID1,I))
0092           ZMIXC(IZID2,I)=DCMPLX(ZMIX(IZID2,I),ZMIXI(IZID2,I))
0093   110   CONTINUE
0094         OLPP=(ZMIXC(IZID1,3)*DCONJG(ZMIXC(IZID2,3))-
0095      &  ZMIXC(IZID1,4)*DCONJG(ZMIXC(IZID2,4)))/2D0
0096         ORPP=DCONJG(OLPP)
0097         XLL2=PMAS(PYCOMP(KSUSY1+IABS(IA)),1)**2
0098         XLR2=XLL2
0099         XRR2=PMAS(PYCOMP(KSUSY2+IABS(IA)),1)**2
0100         XRL2=XRR2
0101         GLIJ=(T3I*ZMIXC(IZID1,2)-TANW*(T3I-EI)*ZMIXC(IZID1,1))*
0102      &  DCONJG(T3I*ZMIXC(IZID2,2)-TANW*(T3I-EI)*ZMIXC(IZID2,1))
0103         GRIJ=ZMIXC(IZID1,1)*DCONJG(ZMIXC(IZID2,1))*(EI*TANW)**2
0104         XM1M2=SMZ(IZID1)*SMZ(IZID2)
0105         QLLS=DCMPLX((T3I-EI*XW)/XW1)*OLPP
0106         QLLU=-GLIJ
0107         QLRS=-DCMPLX((T3I-EI*XW)/XW1)*ORPP
0108         QLRT=DCONJG(GLIJ)
0109         QRLS=-DCMPLX((EI*XW)/XW1)*OLPP
0110         QRLT=GRIJ
0111         QRRS=DCMPLX((EI*XW)/XW1)*ORPP
0112         QRRU=-DCONJG(GRIJ)
0113       ELSEIF(IZID1*IWID2.NE.0.OR.IZID2*IWID1.NE.0) THEN
0114         IF(IZID1.NE.0) THEN
0115           XM1M2=SMZ(IZID1)*SMW(IWID2)
0116           IZID1=IWID2
0117           IZID2=IZID1
0118         ELSE
0119           XM1M2=SMZ(IZID2)*SMW(IWID1)
0120           IZID1=IWID1
0121         ENDIF
0122         RT2I = 1D0/SQRT(2D0)
0123         SQMZ=PMAS(24,1)**2
0124         GMMZ=PMAS(24,1)*PMAS(24,2)
0125         DO 120 I=1,2
0126           VMIXC(IZID1,I)=DCMPLX(VMIX(IZID1,I),VMIXI(IZID1,I))
0127           UMIXC(IZID1,I)=DCMPLX(UMIX(IZID1,I),UMIXI(IZID1,I))
0128   120   CONTINUE
0129         DO 130 I=1,4
0130           ZMIXC(IZID2,I)=DCMPLX(ZMIX(IZID2,I),ZMIXI(IZID2,I))
0131   130   CONTINUE
0132         QLLS=(DCONJG(ZMIXC(IZID2,2))*VMIXC(IZID1,1)-
0133      &  DCONJG(ZMIXC(IZID2,4))*VMIXC(IZID1,2)*RT2I)
0134         QLRS=(ZMIXC(IZID2,2)*DCONJG(UMIXC(IZID1,1))+
0135      &  ZMIXC(IZID2,3)*DCONJG(UMIXC(IZID1,2))*RT2I)
0136         EJ=KCHG(IABS(JA),1)/3D0
0137         T3J=SIGN(1D0,EJ+1D-6)/2D0
0138         QRLS=DCMPLX(0D0,0D0)
0139         QRLT=QRLS
0140         QRRS=QRLS
0141         QRRU=QRLS
0142         XRR2=1D6**2
0143         XRL2=XRR2
0144         XLR2  = PMAS(PYCOMP(KSUSY1+IABS(JA)),1)**2
0145         XLL2  = PMAS(PYCOMP(KSUSY1+IABS(IA)),1)**2
0146         IF(MOD(IA,2).EQ.0) THEN
0147           QLLU=VMIXC(IZID1,1)*DCONJG(ZMIXC(IZID2,1)*(EI-T3I)*
0148      &    TANW+ZMIXC(IZID2,2)*T3I)
0149           QLRT=-DCONJG(UMIXC(IZID1,1))*(
0150      &    ZMIXC(IZID2,1)*(EJ-T3J)*TANW+ZMIXC(IZID2,2)*T3J)
0151         ELSE
0152           QLLU=VMIXC(IZID1,1)*DCONJG(ZMIXC(IZID2,1)*(EJ-T3J)*
0153      &    TANW+ZMIXC(IZID2,2)*T3J)
0154           QLRT=-DCONJG(UMIXC(IZID1,1))*(
0155      &    ZMIXC(IZID2,1)*(EI-T3I)*TANW+ZMIXC(IZID2,2)*T3I)
0156         ENDIF
0157       ELSEIF(IWID1*IWID2.NE.0) THEN
0158         IZID1=IWID1
0159         IZID2=IWID2
0160         XM1M2=SMW(IWID1)*SMW(IWID2)
0161         SQMZ=PMAS(23,1)**2
0162         GMMZ=PMAS(23,1)*PMAS(23,2)
0163         DO 140 I=1,2
0164           VMIXC(IZID1,I)=DCMPLX(VMIX(IZID1,I),VMIXI(IZID1,I))
0165           UMIXC(IZID1,I)=DCMPLX(UMIX(IZID1,I),UMIXI(IZID1,I))
0166           VMIXC(IZID2,I)=DCMPLX(VMIX(IZID2,I),VMIXI(IZID2,I))
0167           UMIXC(IZID2,I)=DCMPLX(UMIX(IZID2,I),UMIXI(IZID2,I))
0168   140   CONTINUE
0169         OLPP=-VMIXC(IZID2,1)*DCONJG(VMIXC(IZID1,1))-
0170      &  VMIXC(IZID2,2)*DCONJG(VMIXC(IZID1,2))/2D0
0171         ORPP=-UMIXC(IZID1,1)*DCONJG(UMIXC(IZID2,1))-
0172      &  UMIXC(IZID1,2)*DCONJG(UMIXC(IZID2,2))/2D0
0173         QRLS=-DCMPLX(EI/XW1)*ORPP
0174         QLLS=DCMPLX((T3I-XW*EI)/XW/XW1)*ORPP
0175         QRRS=-DCMPLX(EI/XW1)*OLPP
0176         QLRS=DCMPLX((T3I-XW*EI)/XW/XW1)*OLPP
0177         IF(MOD(IA,2).EQ.0) THEN
0178           XLR2=PMAS(PYCOMP(KSUSY1+IABS(IA)-1),1)**2
0179           QLRT=-UMIXC(IZID2,1)*DCONJG(UMIXC(IZID1,1))*DCMPLX(T3I/XW)
0180         ELSE
0181           XLR2=PMAS(PYCOMP(KSUSY1+IABS(IA)+1),1)**2
0182           QLRT=-VMIXC(IZID2,1)*DCONJG(VMIXC(IZID1,1))*DCMPLX(T3I/XW)
0183         ENDIF
0184       ELSEIF(MOD(K(N+1,2),KSUSY1).EQ.21.OR.MOD(K(ID,2),KSUSY1).EQ.21)
0185      &THEN
0186         ISKIP=0
0187       ELSE
0188         ISKIP=0
0189       ENDIF
0190
0191       IF(ISKIP.NE.0) THEN
0192         WTMAX=0D0
0193         DO 160 KT=1,100
0194           S12=S12MIN+YJACO1*(KT-1)/99
0195           S23AVE=XM(2)**2+XM(3)**2-(S12+XM(2)**2-XM(1)**2)
0196      &    *(S12+XM(3)**2-XM(5)**2)/(2D0*S12)
0197           S23DF1=(S12-XM(2)**2-XM(1)**2)**2
0198      &    -(2D0*XM(1)*XM(2))**2
0199           S23DF2=(S12-XM(3)**2-XM(5)**2)**2
0200      &    -(2D0*XM(3)*XM(5))**2
0201           S23DF1=S23DF1*EPS
0202           S23DF2=S23DF2*EPS
0203           S23DEL=SQRT(MAX(0D0,S23DF1*S23DF2))/(2D0*S12)
0204           S23DEL=S23DEL/EPS
0205           S23MIN=S23AVE-S23DEL
0206           S23MAX=S23AVE+S23DEL
0207           YJACO2=S23MAX-S23MIN
0208           TH=S12
0209           DO 150 KS=1,100
0210             S23=S23MIN+YJACO2*(KS-1)/99
0211             SH=S23
0212             UH=ZM12+ZM22-SH-TH
0213             WU2 = (UH-ZM12)*(UH-ZM22)
0214             WT2 = (TH-ZM12)*(TH-ZM22)
0215             WS2 = XM1M2*SH
0216             PROPZ2 = (SH-SQMZ)**2 + GMMZ**2
0217             PROPZ=DCMPLX(SH-SQMZ,-GMMZ)/DCMPLX(PROPZ2)
0218             QLL=QLLS*PROPZ+QLLU/DCMPLX(UH-XLL2)
0219             QLR=QLRS*PROPZ+QLRT/DCMPLX(TH-XLR2)
0220             QRL=QRLS*PROPZ+QRLT/DCMPLX(TH-XRL2)
0221             QRR=QRRS*PROPZ+QRRU/DCMPLX(UH-XRR2)
0222             WT0=-((ABS(QLL)**2+ABS(QRR)**2)*WU2+
0223      &      (ABS(QRL)**2+ABS(QLR)**2)*WT2+
0224      &      2D0*DBLE(QLR*DCONJG(QLL)+QRL*DCONJG(QRR))*WS2)
0225             IF(WT0.GT.WTMAX) WTMAX=WT0
0226   150     CONTINUE
0227   160   CONTINUE
0228
0229         WTMAX=WTMAX*1.05D0
0230       ENDIF
0231
0232 C...FIND S12*
0233       AX=S12MIN
0234       CX=S12MAX
0235       BX=S12MIN+0.5D0*YJACO1
0236       X0=AX
0237       X3=CX
0238       IF(ABS(CX-BX).GT.ABS(BX-AX))THEN
0239         X1=BX
0240         X2=BX+C*(CX-BX)
0241       ELSE
0242         X2=BX
0243         X1=BX-C*(BX-AX)
0244       ENDIF
0245
0246 C...SOLVE FOR F1 AND F2
0247       S23DF1=(X1-XM(2)**2-XM(1)**2)**2
0248      &-(2D0*XM(1)*XM(2))**2
0249       S23DF2=(X1-XM(3)**2-XM(5)**2)**2
0250      &-(2D0*XM(3)*XM(5))**2
0251       S23DF1=S23DF1*EPS
0252       S23DF2=S23DF2*EPS
0253       S23DEL=SQRT(MAX(0D0,S23DF1*S23DF2))/(2D0*X1)
0254       F1=-2D0*S23DEL/EPS
0255       S23DF1=(X2-XM(2)**2-XM(1)**2)**2
0256      &-(2D0*XM(1)*XM(2))**2
0257       S23DF2=(X2-XM(3)**2-XM(5)**2)**2
0258      &-(2D0*XM(3)*XM(5))**2
0259       S23DF1=S23DF1*EPS
0260       S23DF2=S23DF2*EPS
0261       S23DEL=SQRT(MAX(0D0,S23DF1*S23DF2))/(2D0*X2)
0262       F2=-2D0*S23DEL/EPS
0263
0264   170 IF(ABS(X3-X0).GT.TOL*(ABS(X1)+ABS(X2)))THEN
0265 C...Possibility of infinite loop with .LT.; changed to .LE. (SKANDS)
0266         IF(F2.LE.F1)THEN
0267           X0=X1
0268           X1=X2
0269           X2=R*X1+C*X3
0270           F1=F2
0271           S23DF1=(X2-XM(2)**2-XM(1)**2)**2
0272      &    -(2D0*XM(1)*XM(2))**2
0273           S23DF2=(X2-XM(3)**2-XM(5)**2)**2
0274      &    -(2D0*XM(3)*XM(5))**2
0275           S23DF1=S23DF1*EPS
0276           S23DF2=S23DF2*EPS
0277           S23DEL=SQRT(MAX(0D0,S23DF1*S23DF2))/(2D0*X2)
0278           F2=-2D0*S23DEL/EPS
0279         ELSE
0280           X3=X2
0281           X2=X1
0282           X1=R*X2+C*X0
0283           F2=F1
0284           S23DF1=(X1-XM(2)**2-XM(1)**2)**2
0285      &    -(2D0*XM(1)*XM(2))**2
0286           S23DF2=(X1-XM(3)**2-XM(5)**2)**2
0287      &    -(2D0*XM(3)*XM(5))**2
0288           S23DF1=S23DF1*EPS
0289           S23DF2=S23DF2*EPS
0290           S23DEL=SQRT(MAX(0D0,S23DF1*S23DF2))/(2D0*X1)
0291           F1=-2D0*S23DEL/EPS
0292         ENDIF
0293         GOTO 170
0294       ENDIF
0295 C...WE WANT THE MAXIMUM, NOT THE MINIMUM
0296       IF(F1.LT.F2)THEN
0297         GOLDEN=-F1
0298         XMIN=X1
0299       ELSE
0300         GOLDEN=-F2
0301         XMIN=X2
0302       ENDIF
0303
0304       IKNT=0
0305   180 S12=S12MIN+PYR(0)*YJACO1
0306       IKNT=IKNT+1
0307 C...GENERATE S23
0308       S23AVE=XM(2)**2+XM(3)**2-(S12+XM(2)**2-XM(1)**2)
0309      &*(S12+XM(3)**2-XM(5)**2)/(2D0*S12)
0310       S23DF1=(S12-XM(2)**2-XM(1)**2)**2
0311      &-(2D0*XM(1)*XM(2))**2
0312       S23DF2=(S12-XM(3)**2-XM(5)**2)**2
0313      &-(2D0*XM(3)*XM(5))**2
0314       S23DF1=S23DF1*EPS
0315       S23DF2=S23DF2*EPS
0316       S23DEL=SQRT(MAX(0D0,S23DF1*S23DF2))/(2D0*S12)
0317       S23DEL=S23DEL/EPS
0318       S23MIN=S23AVE-S23DEL
0319       S23MAX=S23AVE+S23DEL
0320       YJACO2=S23MAX-S23MIN
0321       S23=S23MIN+PYR(0)*YJACO2
0322
0323 C...CHECK THE SAMPLING
0324       IF(IKNT.GT.100) THEN
0325         WRITE(MSTU(11),*) ' IKNT > 100 IN PYTBDY '
0326         GOTO 190
0327       ENDIF
0328       IF(YJACO2.LT.PYR(0)*GOLDEN) GOTO 180
0329
0330       IF(ISKIP.EQ.0) GOTO 190
0331
0332       SH=S23
0333       TH=S12
0334       UH=ZM12+ZM22-SH-TH
0335
0336       WU2 = (UH-ZM12)*(UH-ZM22)
0337       WT2 = (TH-ZM12)*(TH-ZM22)
0338       WS2 = XM1M2*SH
0339       PROPZ2 = (SH-SQMZ)**2 + GMMZ**2
0340       PROPZ=DCMPLX(SH-SQMZ,-GMMZ)/DCMPLX(PROPZ2)
0341
0342       QLL=QLLS*PROPZ+QLLU/DCMPLX(UH-XLL2)
0343       QLR=QLRS*PROPZ+QLRT/DCMPLX(TH-XLR2)
0344       QRL=QRLS*PROPZ+QRLT/DCMPLX(TH-XRL2)
0345       QRR=QRRS*PROPZ+QRRU/DCMPLX(UH-XRR2)
0346 c      QLL=DCMPLX((T3I-EI*XW)/XW1)*OLPP*PROPZ-GLIJ/DCMPLX(UH-XML2)
0347 c      QLR=-DCMPLX((T3I-EI*XW)/XW1)*ORPP*PROPZ+DCONJG(GLIJ)
0348 c     &/DCMPLX(TH-XML2)
0349 c      QRL=-DCMPLX((EI*XW)/XW1)*OLPP*PROPZ+GRIJ/DCMPLX(TH-XMR2)
0350 c      QRR=DCMPLX((EI*XW)/XW1)*ORPP*PROPZ
0351 c     &-DCONJG(GRIJ)/DCMPLX(UH-XMR2)
0352       WT=-((ABS(QLL)**2+ABS(QRR)**2)*WU2+
0353      &(ABS(QRL)**2+ABS(QLR)**2)*WT2+
0354      &2D0*DBLE(QLR*DCONJG(QLL)+QRL*DCONJG(QRR))*WS2)
0355
0356       IF(WT.LT.PYR(0)*WTMAX) GOTO 180
0357       IF(WT.GT.WTMAX) PRINT*,' WT > WTMAX ',WT,WTMAX
0358
0359   190 D3=(XM(5)**2+XM(3)**2-S12)/(2D0*XM(5))
0360       D1=(XM(5)**2+XM(1)**2-S23)/(2D0*XM(5))
0361       D2=XM(5)-D1-D3
0362       P1=SQRT(D1*D1-XM(1)**2)
0363       P2=SQRT(D2*D2-XM(2)**2)
0364       P3=SQRT(D3*D3-XM(3)**2)
0365       CTHE1=2D0*PYR(0)-1D0
0366       ANG1=2D0*PYR(0)*PARU(1)
0367       CPHI1=COS(ANG1)
0368       SPHI1=SIN(ANG1)
0369       ARG=1D0-CTHE1**2
0370       IF(ARG.LT.0D0.AND.ARG.GT.-1D-3) ARG=0D0
0371       STHE1=SQRT(ARG)
0372       P(N+1,1)=P1*STHE1*CPHI1
0373       P(N+1,2)=P1*STHE1*SPHI1
0374       P(N+1,3)=P1*CTHE1
0375       P(N+1,4)=D1
0376
0377 C...GET CPHI3
0378       ANG3=2D0*PYR(0)*PARU(1)
0379       CPHI3=COS(ANG3)
0380       SPHI3=SIN(ANG3)
0381       CTHE3=(P2**2-P1**2-P3**2)/2D0/P1/P3
0382       ARG=1D0-CTHE3**2
0383       IF(ARG.LT.0D0.AND.ARG.GT.-1D-3) ARG=0D0
0384       STHE3=SQRT(ARG)
0385       P(N+3,1)=-P3*STHE3*CPHI3*CTHE1*CPHI1
0386      &+P3*STHE3*SPHI3*SPHI1
0387      &+P3*CTHE3*STHE1*CPHI1
0388       P(N+3,2)=-P3*STHE3*CPHI3*CTHE1*SPHI1
0389      &-P3*STHE3*SPHI3*CPHI1
0390      &+P3*CTHE3*STHE1*SPHI1
0391       P(N+3,3)=P3*STHE3*CPHI3*STHE1
0392      &+P3*CTHE3*CTHE1
0393       P(N+3,4)=D3
0394
0395       DO 200 I=1,3
0396         P(N+2,I)=-P(N+1,I)-P(N+3,I)
0397   200 CONTINUE
0398       P(N+2,4)=D2
0399
0400       RETURN
0401       END