pythiaeRHIC/pythia/pymsin.f

0001
0002 C*********************************************************************
0003
0004 C...PYMSIN
0005 C...Initializes supersymmetry: finds sparticle masses and
0006 C...branching ratios and stores this information.
0007 C...AUTHOR: STEPHEN MRENNA
0008 C...Author: P. Skands (SLHA + RPV + ISASUSY Interface, NMSSM)
0009
0010       SUBROUTINE PYMSIN
0011
0012 C...Double precision and integer declarations.
0013       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0014       IMPLICIT INTEGER(I-N)
0015       INTEGER PYK,PYCHGE,PYCOMP
0016 C...Parameter statement to help give large particle numbers.
0017       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0018      &KEXCIT=4000000,KDIMEN=5000000)
0019 C...Commonblocks.
0020       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0021       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0022       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0023       COMMON/PYDAT4/CHAF(500,2)
0024       CHARACTER CHAF*16
0025       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0026       COMMON/PYINT4/MWID(500),WIDS(500,5)
0027       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
0028       COMMON/PYMSRV/RVLAM(3,3,3), RVLAMP(3,3,3), RVLAMB(3,3,3)
0029       COMMON/PYSSMT/ZMIX(4,4),UMIX(2,2),VMIX(2,2),SMZ(4),SMW(2),
0030      &SFMIX(16,4),ZMIXI(4,4),UMIXI(2,2),VMIXI(2,2)
0031       COMMON/PYHTRI/HHH(7)
0032       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYPARS/,/PYINT4/,
0033      &/PYMSSM/,/PYMSRV/,/PYSSMT/
0034
0035 C...Local variables.
0036       DOUBLE PRECISION ALFA,BETA
0037       DOUBLE PRECISION TANB,AL,BE,COSA,COSB,SINA,SINB,XW
0038       INTEGER I,J,J1,I1,K1
0039       INTEGER KC,LKNT,IDLAM(400,3)
0040       DOUBLE PRECISION XLAM(0:400)
0041       DOUBLE PRECISION WDTP(0:400),WDTE(0:400,0:5)
0042       DOUBLE PRECISION XARG,COS2B,XMW2,XMZ2
0043       DOUBLE PRECISION DELM,XMDIF
0044       DOUBLE PRECISION DX,DY,DS,DMU2,DMA2,DQ2,DU2,DD2,DL2,DE2,DHU2,DHD2
0045       DOUBLE PRECISION ARG,SGNMU,R
0046       INTEGER IMSSM
0047       INTEGER IRPRTY
0048       INTEGER KFSUSY(50),MWIDSU(36),MDCYSU(36)
0049       SAVE MWIDSU,MDCYSU
0050       DATA KFSUSY/
0051      &1000001,2000001,1000002,2000002,1000003,2000003,
0052      &1000004,2000004,1000005,2000005,1000006,2000006,
0053      &1000011,2000011,1000012,2000012,1000013,2000013,
0054      &1000014,2000014,1000015,2000015,1000016,2000016,
0055      &1000021,1000022,1000023,1000025,1000035,1000024,
0056      &1000037,1000039,     25,     35,     36,     37,
0057      &      6,     24,     45,     46,1000045, 9*0/
0058       DATA INIT/0/
0059
0060 C...Do nothing if SUSY not requested.
0061       IMSSM=IMSS(1)
0062       IF(IMSSM.EQ.0) RETURN
0063
0064 C...Save copy of MWID(KC) and MDCY(KC,1) values before
0065 C...they are set to zero for the LSP.
0066       IF(INIT.EQ.0) THEN
0067         INIT=1
0068         DO 100 I=1,36
0069           KF=KFSUSY(I)
0070           KC=PYCOMP(KF)
0071           MWIDSU(I)=MWID(KC)
0072           MDCYSU(I)=MDCY(KC,1)
0073   100   CONTINUE
0074       ENDIF
0075
0076 C...Restore MWID(KC) and MDCY(KC,1) values previously zeroed for LSP.
0077       DO 110 I=1,36
0078         KF=KFSUSY(I)
0079         KC=PYCOMP(KF)
0080         IF(MDCY(KC,1).EQ.0.AND.MDCYSU(I).NE.0) THEN
0081           MWID(KC)=MWIDSU(I)
0082           MDCY(KC,1)=MDCYSU(I)
0083         ENDIF
0084   110 CONTINUE
0085
0086 C...First part of routine: set masses and couplings.
0087
0088 C...Reset mixing values in sfermion sector to pure left/right.
0089       DO 120 I=1,16
0090         SFMIX(I,1)=1D0
0091         SFMIX(I,4)=1D0
0092         SFMIX(I,2)=0D0
0093         SFMIX(I,3)=0D0
0094   120 CONTINUE
0095
0096 C...Add NMSSM states if NMSSM switched on, and change old names.
0097       IF (IMSS(13).NE.0) THEN
0098 C...  Switch on NMSSM
0099         WRITE(MSTU(11),*) '(PYMSIN:) switching on NMSSM'
0100
0101         KFN=25
0102         KCN=KFN
0103         CHAF(KCN,1)='H_10'
0104         CHAF(KCN,2)=' '
0105
0106         KFN=35
0107         KCN=KFN
0108         CHAF(KCN,1)='H_20'
0109         CHAF(KCN,2)=' '
0110
0111         KFN=45
0112         KCN=KFN
0113         CHAF(KCN,1)='H_30'
0114         CHAF(KCN,2)=' '
0115
0116         KFN=36
0117         KCN=KFN
0118         CHAF(KCN,1)='A_10'
0119         CHAF(KCN,2)=' '
0120
0121         KFN=46
0122         KCN=KFN
0123         CHAF(KCN,1)='A_20'
0124         CHAF(KCN,2)=' '
0125
0126         KFN=1000045
0127         KCN=PYCOMP(KFN)
0128         IF (KCN.EQ.0) THEN
0129           DO 123 KCT=100,MSTU(6)
0130             IF(KCHG(KCT,4).GT.100) KCN=KCT
0131  123      CONTINUE
0132           KCN=KCN+1
0133           KCHG(KCN,4)=KFN
0134           MSTU(20)=0
0135         ENDIF
0136 C...  Set stable for now
0137         PMAS(KCN,2)=1D-6
0138         MWID(KCN)=0
0139         MDCY(KCN,1)=0
0140         MDCY(KCN,2)=0
0141         MDCY(KCN,3)=0
0142         CHAF(KCN,1)='~chi_50'
0143         CHAF(KCN,2)=' '
0144       ENDIF
0145
0146 C...Read spectrum from SLHA file.
0147       IF (IMSSM.EQ.11.AND.IMSS(21).NE.0) THEN
0148 C...First check for new states
0149         CALL PYSLHA(0,0,IFAIL)
0150 C...Then read spectrum
0151         CALL PYSLHA(1,0,IFAIL)
0152       ELSEIF (IMSS(21).NE.0) THEN
0153 C...Check for new states but don't read spectrum
0154         CALL PYSLHA(0,0,IFAIL)
0155       ENDIF
0156
0157 C...Common couplings.
0158       TANB=RMSS(5)
0159       BETA=ATAN(TANB)
0160       COSB=COS(BETA)
0161       SINB=TANB*COSB
0162       COS2B=COS(2D0*BETA)
0163       ALFA=RMSS(18)
0164       XMW2=PMAS(24,1)**2
0165       XMZ2=PMAS(23,1)**2
0166       XW=PARU(102)
0167
0168 C...Define sparticle masses for a general MSSM simulation.
0169       IF(IMSSM.EQ.1) THEN
0170         IF(IMSS(9).EQ.0) RMSS(22)=RMSS(9)
0171         DO 130 I=1,5,2
0172           KC=PYCOMP(KSUSY1+I)
0173           PMAS(KC,1)=SQRT(RMSS(8)**2-(2D0*XMW2+XMZ2)*COS2B/6D0)
0174           KC=PYCOMP(KSUSY2+I)
0175           PMAS(KC,1)=SQRT(RMSS(9)**2+(XMW2-XMZ2)*COS2B/3D0)
0176           KC=PYCOMP(KSUSY1+I+1)
0177           PMAS(KC,1)=SQRT(RMSS(8)**2+(4D0*XMW2-XMZ2)*COS2B/6D0)
0178           KC=PYCOMP(KSUSY2+I+1)
0179           PMAS(KC,1)=SQRT(RMSS(22)**2-(XMW2-XMZ2)*COS2B*2D0/3D0)
0180   130   CONTINUE
0181         XARG=RMSS(6)**2-PMAS(24,1)**2*ABS(COS(2D0*BETA))
0182         IF(XARG.LT.0D0) THEN
0183           WRITE(MSTU(11),*) ' SNEUTRINO MASS IS NEGATIVE'//
0184      &    ' FROM THE SUM RULE. '
0185           WRITE(MSTU(11),*) '  TRY A SMALLER VALUE OF TAN(BETA). '
0186           RETURN
0187         ELSE
0188           XARG=SQRT(XARG)
0189         ENDIF
0190         DO 140 I=11,15,2
0191           PMAS(PYCOMP(KSUSY1+I),1)=RMSS(6)
0192           PMAS(PYCOMP(KSUSY2+I),1)=RMSS(7)
0193           PMAS(PYCOMP(KSUSY1+I+1),1)=XARG
0194           PMAS(PYCOMP(KSUSY2+I+1),1)=9999D0
0195   140   CONTINUE
0196         IF(IMSS(8).EQ.1) THEN
0197           RMSS(13)=RMSS(6)
0198           RMSS(14)=RMSS(7)
0199         ENDIF
0200
0201 C...Alternatively derive masses from SUGRA relations.
0202       ELSEIF(IMSSM.EQ.2) THEN
0203         RMSS(36)=RMSS(16)
0204         CALL PYAPPS
0205 C...Or use ISASUSY
0206       ELSEIF(IMSSM.EQ.12.OR.IMSSM.EQ.13) THEN
0207         RMSS(36)=RMSS(16)
0208         CALL PYSUGI
0209         ALFA=RMSS(18)
0210         GOTO 170
0211       ELSE
0212         GOTO 170
0213       ENDIF
0214
0215 C...Add in extra D-term contributions.
0216       IF(IMSS(7).EQ.1) THEN
0217         R=0.43D0
0218         DX=RMSS(23)
0219         DY=RMSS(24)
0220         DS=RMSS(25)
0221         WRITE(MSTU(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
0222         WRITE(MSTU(11),*) 'C  NEW DTERMS ADDED TO SCALAR MASSES   '
0223         WRITE(MSTU(11),*) 'C   IN A U(B-L) THEORY                 '
0224         WRITE(MSTU(11),*) 'C   DX = ',DX
0225         WRITE(MSTU(11),*) 'C   DY = ',DY
0226         WRITE(MSTU(11),*) 'C   DS = ',DS
0227         WRITE(MSTU(11),*) 'C                                      '
0228         DY=R*DY-4D0/33D0*(1D0-R)*DX+(1D0-R)/33D0*DS
0229         WRITE(MSTU(11),*) 'C   DY AT THE WEAK SCALE = ',DY
0230         WRITE(MSTU(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
0231         DQ2=DY/6D0-DX/3D0-DS/3D0
0232         DU2=-2D0*DY/3D0-DX/3D0-DS/3D0
0233         DD2=DY/3D0+DX-2D0*DS/3D0
0234         DL2=-DY/2D0+DX-2D0*DS/3D0
0235         DE2=DY-DX/3D0-DS/3D0
0236         DHU2=DY/2D0+2D0*DX/3D0+2D0*DS/3D0
0237         DHD2=-DY/2D0-2D0*DX/3D0+DS
0238         DMU2=(-DY/2D0-2D0/3D0*DX+(COSB**2-2D0*SINB**2/3D0)*DS)
0239      &  /ABS(COS2B)
0240         DMA2 = 2D0*DMU2+DHU2+DHD2
0241         DO 150 I=1,5,2
0242           KC=PYCOMP(KSUSY1+I)
0243           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DQ2)
0244           KC=PYCOMP(KSUSY2+I)
0245           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DD2)
0246           KC=PYCOMP(KSUSY1+I+1)
0247           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DQ2)
0248           KC=PYCOMP(KSUSY2+I+1)
0249           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DU2)
0250   150   CONTINUE
0251         DO 160 I=11,15,2
0252           KC=PYCOMP(KSUSY1+I)
0253           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DL2)
0254           KC=PYCOMP(KSUSY2+I)
0255           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DE2)
0256           KC=PYCOMP(KSUSY1+I+1)
0257           PMAS(KC,1)=SQRT(PMAS(KC,1)**2+DL2)
0258   160   CONTINUE
0259         IF(RMSS(4)**2+DMU2.LT.0D0) THEN
0260           WRITE(MSTU(11),*) ' MU2 DRIVEN NEGATIVE '
0261           CALL PYSTOP(104)
0262         ENDIF
0263         SGNMU=SIGN(1D0,RMSS(4))
0264         RMSS(4)=SGNMU*SQRT(RMSS(4)**2+DMU2)
0265         ARG=RMSS(10)**2*SIGN(1D0,RMSS(10))+DQ2
0266         RMSS(10)=SIGN(SQRT(ABS(ARG)),ARG)
0267         ARG=RMSS(11)**2*SIGN(1D0,RMSS(11))+DD2
0268         RMSS(11)=SIGN(SQRT(ABS(ARG)),ARG)
0269         ARG=RMSS(12)**2*SIGN(1D0,RMSS(12))+DU2
0270         RMSS(12)=SIGN(SQRT(ABS(ARG)),ARG)
0271         ARG=RMSS(13)**2*SIGN(1D0,RMSS(13))+DL2
0272         RMSS(13)=SIGN(SQRT(ABS(ARG)),ARG)
0273         ARG=RMSS(14)**2*SIGN(1D0,RMSS(14))+DE2
0274         RMSS(14)=SIGN(SQRT(ABS(ARG)),ARG)
0275         IF( RMSS(19)**2 + DMA2 .LE. 50D0 ) THEN
0276           WRITE(MSTU(11),*) ' MA DRIVEN TOO LOW '
0277           CALL PYSTOP(104)
0278         ENDIF
0279         RMSS(19)=SQRT(RMSS(19)**2+DMA2)
0280         RMSS(6)=SQRT(RMSS(6)**2+DL2)
0281         RMSS(7)=SQRT(RMSS(7)**2+DE2)
0282         WRITE(MSTU(11),*) ' MTL = ',RMSS(10)
0283         WRITE(MSTU(11),*) ' MBR = ',RMSS(11)
0284         WRITE(MSTU(11),*) ' MTR = ',RMSS(12)
0285         WRITE(MSTU(11),*) ' SEL = ',RMSS(6),RMSS(13)
0286         WRITE(MSTU(11),*) ' SER = ',RMSS(7),RMSS(14)
0287       ENDIF
0288
0289 C...Fix the third generation sfermions.
0290       CALL PYTHRG
0291
0292 C...Fix the neutralino--chargino--gluino sector.
0293       CALL PYINOM
0294
0295 C...Fix the Higgs sector.
0296       CALL PYHGGM(ALFA)
0297
0298 C...Choose the Gunion-Haber convention.
0299       ALFA=-ALFA
0300       RMSS(18)=ALFA
0301
0302 C...Print information on mass parameters.
0303       IF(IMSSM.EQ.2.AND.MSTP(122).GT.0) THEN
0304         WRITE(MSTU(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
0305         WRITE(MSTU(11),*) ' USING APPROXIMATE SUGRA RELATIONS '
0306         WRITE(MSTU(11),*) ' M0 = ',RMSS(8)
0307         WRITE(MSTU(11),*) ' M1/2=',RMSS(1)
0308         WRITE(MSTU(11),*) ' TANB=',RMSS(5)
0309         WRITE(MSTU(11),*) ' MU = ',RMSS(4)
0310         WRITE(MSTU(11),*) ' AT = ',RMSS(16)
0311         WRITE(MSTU(11),*) ' MA = ',RMSS(19)
0312         WRITE(MSTU(11),*) ' MTOP=',PMAS(6,1)
0313         WRITE(MSTU(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
0314       ENDIF
0315       IF(IMSS(20).EQ.1) THEN
0316         WRITE(MSTU(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
0317         WRITE(MSTU(11),*) ' DEBUG MODE '
0318         WRITE(MSTU(11),*) ' UMIX = ',UMIX(1,1),UMIX(1,2),
0319      &  UMIX(2,1),UMIX(2,2)
0320         WRITE(MSTU(11),*) ' UMIXI = ',UMIXI(1,1),UMIXI(1,2),
0321      &  UMIXI(2,1),UMIXI(2,2)
0322         WRITE(MSTU(11),*) ' VMIX = ',VMIX(1,1),VMIX(1,2),
0323      &  VMIX(2,1),VMIX(2,2)
0324         WRITE(MSTU(11),*) ' VMIXI = ',VMIXI(1,1),VMIXI(1,2),
0325      &  VMIXI(2,1),VMIXI(2,2)
0326         WRITE(MSTU(11),*) ' ZMIX = ',(ZMIX(1,I),I=1,4)
0327         WRITE(MSTU(11),*) ' ZMIXI = ',(ZMIXI(1,I),I=1,4)
0328         WRITE(MSTU(11),*) ' ZMIX = ',(ZMIX(2,I),I=1,4)
0329         WRITE(MSTU(11),*) ' ZMIXI = ',(ZMIXI(2,I),I=1,4)
0330         WRITE(MSTU(11),*) ' ZMIX = ',(ZMIX(3,I),I=1,4)
0331         WRITE(MSTU(11),*) ' ZMIXI = ',(ZMIXI(3,I),I=1,4)
0332         WRITE(MSTU(11),*) ' ZMIX = ',(ZMIX(4,I),I=1,4)
0333         WRITE(MSTU(11),*) ' ZMIXI = ',(ZMIXI(4,I),I=1,4)
0334         WRITE(MSTU(11),*) ' ALFA = ',ALFA
0335         WRITE(MSTU(11),*) ' BETA = ',BETA
0336         WRITE(MSTU(11),*) ' STOP = ',(SFMIX(6,I),I=1,4)
0337         WRITE(MSTU(11),*) ' SBOT = ',(SFMIX(5,I),I=1,4)
0338         WRITE(MSTU(11),*) 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
0339       ENDIF
0340
0341 C...Set up the Higgs couplings - needed here since initialization
0342 C...in PYINRE did not yet occur when PYWIDT is called below.
0343   170 AL=ALFA
0344       BE=BETA
0345       SINA=SIN(AL)
0346       COSA=COS(AL)
0347       COSB=COS(BE)
0348       SINB=TANB*COSB
0349       SBMA=SIN(BE-AL)
0350       SAPB=SIN(AL+BE)
0351       CAPB=COS(AL+BE)
0352       CBMA=COS(BE-AL)
0353       C2A=COS(2D0*AL)
0354       C2B=COSB**2-SINB**2
0355 C...tanb (used for H+)
0356       PARU(141)=TANB
0357
0358 C...Firstly: h
0359 C...Coupling to d-type quarks
0360       PARU(161)=SINA/COSB
0361 C...Coupling to u-type quarks
0362       PARU(162)=-COSA/SINB
0363 C...Coupling to leptons
0364       PARU(163)=PARU(161)
0365 C...Coupling to Z
0366       PARU(164)=SBMA
0367 C...Coupling to W
0368       PARU(165)=PARU(164)
0369
0370 C...Secondly: H
0371 C...Coupling to d-type quarks
0372       PARU(171)=-COSA/COSB
0373 C...Coupling to u-type quarks
0374       PARU(172)=-SINA/SINB
0375 C...Coupling to leptons
0376       PARU(173)=PARU(171)
0377 C...Coupling to Z
0378       PARU(174)=CBMA
0379 C...Coupling to W
0380       PARU(175)=PARU(174)
0381 C...Coupling to h
0382       IF(IMSS(4).GE.2) THEN
0383         PARU(176)=COS(2D0*AL)*COS(BE+AL)-2D0*SIN(2D0*AL)*SIN(BE+AL)
0384       ELSE
0385         HHH(3)=HHH(3)+HHH(4)+HHH(5)
0386         PARU(176)=-3D0/HHH(1)*(HHH(1)*SINA**2*COSB*COSA+
0387      1  HHH(2)*COSA**2*SINB*SINA+HHH(3)*(SINA**3*SINB+COSA**3*COSB-
0388      2  2D0/3D0*CBMA)-HHH(6)*SINA*(COSB*C2A+COSA*CAPB)+
0389      3  HHH(7)*COSA*(SINB*C2A+SINA*CAPB))
0390       ENDIF
0391 C...Coupling to H+
0392 C...Define later
0393       IF(IMSS(4).GE.2) THEN
0394         PARU(168)=-SBMA-COS(2D0*BE)*SAPB/2D0/(1D0-XW)
0395       ELSE
0396         PARU(168)=1D0/HHH(1)*(HHH(1)*SINB**2*COSB*SINA-
0397      1 HHH(2)*COSB**2*SINB*COSA-HHH(3)*(SINB**3*COSA-COSB**3*SINA)+
0398      2 2D0*HHH(5)*SBMA-HHH(6)*SINB*(COSB*SAPB+SINA*C2B)-
0399      3 HHH(7)*COSB*(COSA*C2B-SINB*SAPB)-(HHH(5)-HHH(4))*SBMA)
0400       ENDIF
0401 C...Coupling to A
0402       IF(IMSS(4).GE.2) THEN
0403         PARU(177)=COS(2D0*BE)*COS(BE+AL)
0404       ELSE
0405         PARU(177)=-1D0/HHH(1)*(HHH(1)*SINB**2*COSB*COSA+
0406      1 HHH(2)*COSB**2*SINB*SINA+HHH(3)*(SINB**3*SINA+COSB**3*COSA)-
0407      2 2D0*HHH(5)*CBMA-HHH(6)*SINB*(COSB*CAPB+COSA*C2B)+
0408      3 HHH(7)*COSB*(SINB*CAPB+SINA*C2B))
0409       ENDIF
0410 C...Coupling to H+
0411       IF(IMSS(4).GE.2) THEN
0412         PARU(178)=PARU(177)
0413       ELSE
0414         PARU(178)=PARU(177)-(HHH(5)-HHH(4))/HHH(1)*CBMA
0415       ENDIF
0416 C...Thirdly, A
0417 C...Coupling to d-type quarks
0418       PARU(181)=TANB
0419 C...Coupling to u-type quarks
0420       PARU(182)=1D0/PARU(181)
0421 C...Coupling to leptons
0422       PARU(183)=PARU(181)
0423       PARU(184)=0D0
0424       PARU(185)=0D0
0425 C...Coupling to Z h
0426       PARU(186)=COS(BE-AL)
0427 C...Coupling to Z H
0428       PARU(187)=SIN(BE-AL)
0429       PARU(188)=0D0
0430       PARU(189)=0D0
0431       PARU(190)=0D0
0432
0433 C...Finally: H+
0434 C...Coupling to W h
0435       PARU(195)=COS(BE-AL)
0436
0437 C...Tell that all Higgs couplings have been set.
0438       MSTP(4)=1
0439
0440 C...Set R-Violating couplings.
0441 C...Set lambda couplings to common value or "natural values".
0442       IF ((IMSS(51).NE.3).AND.(IMSS(51).NE.0)) THEN
0443         VIR3=1D0/(126D0)**3
0444         DO 200 IRK=1,3
0445           DO 190 IRI=1,3
0446             DO 180 IRJ=1,3
0447               IF (IRI.NE.IRJ) THEN
0448                 IF (IRI.LT.IRJ) THEN
0449                   RVLAM(IRI,IRJ,IRK)=RMSS(51)
0450                   IF (IMSS(51).EQ.2) RVLAM(IRI,IRJ,IRK)=RMSS(51)*
0451      &              SQRT(PMAS(9+2*IRI,1)*PMAS(9+2*IRJ,1)*
0452      &              PMAS(9+2*IRK,1)*VIR3)
0453                 ELSE
0454                   RVLAM(IRI,IRJ,IRK)=-RVLAM(IRJ,IRI,IRK)
0455                 ENDIF
0456               ELSE
0457                 RVLAM(IRI,IRJ,IRK)=0D0
0458               ENDIF
0459   180       CONTINUE
0460   190     CONTINUE
0461   200   CONTINUE
0462       ENDIF
0463 C...Set lambda' couplings to common value or "natural values".
0464       IF ((IMSS(52).NE.3).AND.(IMSS(52).NE.0)) THEN
0465         VIR3=1D0/(126D0)**3
0466         DO 230 IRI=1,3
0467           DO 220 IRJ=1,3
0468             DO 210 IRK=1,3
0469               RVLAMP(IRI,IRJ,IRK)=RMSS(52)
0470               IF (IMSS(52).EQ.2) RVLAMP(IRI,IRJ,IRK)=RMSS(52)*
0471      &          SQRT(PMAS(9+2*IRI,1)*0.5D0*(PMAS(2*IRJ,1)+
0472      &          PMAS(2*IRJ-1,1))*PMAS(2*IRK-1,1)*VIR3)
0473   210       CONTINUE
0474   220     CONTINUE
0475   230   CONTINUE
0476       ENDIF
0477 C...Set lambda'' couplings to common value or "natural values".
0478       IF ((IMSS(53).NE.3).AND.(IMSS(53).NE.0)) THEN
0479         VIR3=1D0/(126D0)**3
0480         DO 260 IRI=1,3
0481           DO 250 IRJ=1,3
0482             DO 240 IRK=1,3
0483               IF (IRJ.NE.IRK) THEN
0484                 IF (IRJ.LT.IRK) THEN
0485                   RVLAMB(IRI,IRJ,IRK)=RMSS(53)
0486                   IF (IMSS(53).EQ.2) RVLAMB(IRI,IRJ,IRK)=
0487      &              RMSS(53)*SQRT(PMAS(2*IRI,1)*PMAS(2*IRJ-1,1)*
0488      &              PMAS(2*IRK-1,1)*VIR3)
0489                 ELSE
0490                   RVLAMB(IRI,IRJ,IRK)=-RVLAMB(IRI,IRK,IRJ)
0491                 ENDIF
0492               ELSE
0493                 RVLAMB(IRI,IRJ,IRK) = 0D0
0494               ENDIF
0495   240       CONTINUE
0496   250     CONTINUE
0497   260   CONTINUE
0498       ENDIF
0499
0500 C...Antisymmetrize couplings set by user
0501       IF (IMSS(51).EQ.3.OR.IMSS(53).EQ.3) THEN
0502         DO 290 IRI=1,3
0503           DO 280 IRJ=1,3
0504             DO 270 IRK=1,3
0505               IF (RVLAM(IRI,IRJ,IRK).NE.-RVLAM(IRJ,IRI,IRK)) THEN
0506                 RVLAM(IRJ,IRI,IRK)=-RVLAM(IRI,IRJ,IRK)
0507                 IF (IRI.EQ.IRJ) RVLAM(IRI,IRJ,IRK)=0D0
0508               ENDIF
0509               IF (RVLAMB(IRI,IRJ,IRK).NE.-RVLAMB(IRI,IRK,IRJ)) THEN
0510                 RVLAMB(IRI,IRK,IRJ)=-RVLAMB(IRI,IRJ,IRK)
0511                 IF (IRJ.EQ.IRK) RVLAMB(IRI,IRJ,IRK)=0D0
0512               ENDIF
0513   270       CONTINUE
0514   280     CONTINUE
0515   290   CONTINUE
0516       ENDIF
0517
0518 C...Write spectrum to SLHA file
0519       IF (IMSS(23).NE.0) THEN
0520         IFAIL=0
0521         CALL PYSLHA(3,0,IFAIL)
0522       ENDIF
0523
0524 C...Second part of routine: set decay modes and branching ratios.
0525
0526 C...Allow chi10 -> gravitino + gamma or not.
0527       KC=PYCOMP(KSUSY1+39)
0528       IF( IMSS(11) .NE. 0 ) THEN
0529         PMAS(KC,1)=RMSS(21)/1D9
0530         PMAS(KC,2)=0D0
0531         IRPRTY=0
0532         WRITE(MSTU(11),*) ' ALLOWING DECAYS TO GRAVITINOS '
0533       ELSE IF (IMSS(51).GE.1.OR.IMSS(52).GE.1.OR.IMSS(53).GE.1) THEN
0534         IRPRTY=0
0535         IF (IMSS(51).GE.1) WRITE(MSTU(11),*)
0536      &       ' ALLOWING SUSY LLE DECAYS'
0537         IF (IMSS(52).GE.1) WRITE(MSTU(11),*)
0538      &       ' ALLOWING SUSY LQD DECAYS'
0539         IF (IMSS(53).GE.1) WRITE(MSTU(11),*)
0540      &       ' ALLOWING SUSY UDD DECAYS'
0541         IF (IMSS(53).GE.1.AND.IMSS(52).GE.1) WRITE(MSTU(11),*)
0542      &   ' --- Warning: R-Violating couplings possibly',
0543      &       ' incompatible with proton decay'
0544       ELSE
0545         PMAS(KC,1)=9999D0
0546         IRPRTY=1
0547       ENDIF
0548
0549 C...Loop over sparticle and Higgs species.
0550       PMCHI1=PMAS(PYCOMP(KSUSY1+22),1)
0551 C...Find the LSP or NLSP for a gravitino LSP
0552       ILSP=0
0553       PMLSP=1D20
0554       DO 300 I=1,36
0555         KF=KFSUSY(I)
0556         IF(KF.EQ.1000039) GOTO 300
0557         KC=PYCOMP(KF)
0558         IF(PMAS(KC,1).LT.PMLSP) THEN
0559           ILSP=I
0560           PMLSP=PMAS(KC,1)
0561         ENDIF
0562   300 CONTINUE
0563       DO 370 I=1,50
0564         IF (I.GT.39.AND.IMSS(13).NE.1) GOTO 370
0565         KF=KFSUSY(I)
0566         IF (KF.EQ.0) GOTO 370
0567         KC=PYCOMP(KF)
0568         LKNT=0
0569
0570 C...Check if there are any decays listed for this sparticle
0571 C...in a file
0572         IF (IMSS(22).NE.0) THEN
0573           IFAIL=0
0574 C...First look for MASS entry if not already done
0575           IF (IMSS(1).NE.11.AND.IMSS(21).NE.0) CALL PYSLHA(5,KF,IFAIL)
0576 C...Then look for decay info
0577           IFAIL=0
0578           CALL PYSLHA(2,KF,IFAIL)
0579           IF (IFAIL.EQ.0.OR.KF.EQ.6.OR.KF.EQ.24) GOTO 370
0580         ELSEIF (I.GE.37) THEN
0581           GOTO 370
0582         ENDIF
0583
0584 C...Sfermion decays.
0585         IF(I.LE.24) THEN
0586 C...First check to see if sneutrino is lighter than chi10.
0587           IF((I.EQ.15.OR.I.EQ.19.OR.I.EQ.23).AND.
0588      &    PMAS(KC,1).LT.PMCHI1) THEN
0589           ELSE
0590             CALL PYSFDC(KF,XLAM,IDLAM,LKNT)
0591           ENDIF
0592
0593 C...Gluino decays.
0594         ELSEIF(I.EQ.25) THEN
0595           CALL PYGLUI(KF,XLAM,IDLAM,LKNT)
0596           IF(I.EQ.ILSP.AND.IRPRTY.EQ.1) LKNT=0
0597
0598 C...Neutralino decays.
0599         ELSEIF(I.GE.26.AND.I.LE.29) THEN
0600           CALL PYNJDC(KF,XLAM,IDLAM,LKNT)
0601 C...chi10 stable or chi10 -> gravitino + gamma.
0602           IF(I.EQ.26.AND.IRPRTY.EQ.1) THEN
0603             PMAS(KC,2)=1D-6
0604             MDCY(KC,1)=0
0605             MWID(KC)=0
0606           ENDIF
0607
0608 C...Chargino decays.
0609         ELSEIF(I.GE.30.AND.I.LE.31) THEN
0610           CALL PYCJDC(KF,XLAM,IDLAM,LKNT)
0611
0612 C...Gravitino is stable.
0613         ELSEIF(I.EQ.32) THEN
0614           MDCY(KC,1)=0
0615           MWID(KC)=0
0616
0617 C...Higgs decays.
0618         ELSEIF(I.GE.33.AND.I.LE.36) THEN
0619 C...Calculate decays to non-SUSY particles.
0620           CALL PYWIDT(KF,PMAS(KC,1)**2,WDTP,WDTE)
0621           LKNT=0
0622           DO 310 I1=0,100
0623             XLAM(I1)=0D0
0624   310     CONTINUE
0625           DO 330 I1=1,MDCY(KC,3)
0626             K1=MDCY(KC,2)+I1-1
0627             IF(IABS(KFDP(K1,1)).GT.KSUSY1.OR.
0628      &      IABS(KFDP(K1,2)).GT.KSUSY1) GOTO 330
0629             XLAM(I1)=WDTP(I1)
0630             XLAM(0)=XLAM(0)+XLAM(I1)
0631             DO 320 J1=1,3
0632               IDLAM(I1,J1)=KFDP(K1,J1)
0633   320       CONTINUE
0634             LKNT=LKNT+1
0635   330     CONTINUE
0636 C...Add the decays to SUSY particles.
0637           CALL PYHEXT(KF,XLAM,IDLAM,LKNT)
0638         ENDIF
0639 C...Zero the branching ratios for use in loop mode
0640 C...thanks to K. Matchev (FNAL)
0641         DO 340 IDC=MDCY(KC,2),MDCY(KC,2)+MDCY(KC,3)-1
0642           BRAT(IDC)=0D0
0643   340   CONTINUE
0644
0645 C...Set stable particles.
0646         IF(LKNT.EQ.0) THEN
0647           MDCY(KC,1)=0
0648           MWID(KC)=0
0649           PMAS(KC,2)=1D-6
0650           PMAS(KC,3)=1D-5
0651           PMAS(KC,4)=0D0
0652
0653 C...Store branching ratios in the standard tables.
0654         ELSE
0655           IDC=MDCY(KC,2)+MDCY(KC,3)-1
0656           DELM=1D6
0657           DO 360 IL=1,LKNT
0658             IDCSV=IDC
0659   350       IDC=IDC+1
0660             BRAT(IDC)=0D0
0661             IF(IDC.EQ.MDCY(KC,2)+MDCY(KC,3)) IDC=MDCY(KC,2)
0662             IF(IDLAM(IL,1).EQ.KFDP(IDC,1).AND.IDLAM(IL,2).EQ.
0663      &      KFDP(IDC,2).AND.IDLAM(IL,3).EQ.KFDP(IDC,3)) THEN
0664               BRAT(IDC)=XLAM(IL)/XLAM(0)
0665               XMDIF=PMAS(KC,1)
0666               IF(MDME(IDC,1).GE.1) THEN
0667                 XMDIF=XMDIF-PMAS(PYCOMP(KFDP(IDC,1)),1)-
0668      &          PMAS(PYCOMP(KFDP(IDC,2)),1)
0669                 IF(KFDP(IDC,3).NE.0) XMDIF=XMDIF-
0670      &          PMAS(PYCOMP(KFDP(IDC,3)),1)
0671               ENDIF
0672               IF(I.LE.32) THEN
0673                 IF(XMDIF.GE.0D0) THEN
0674                   DELM=MIN(DELM,XMDIF)
0675                 ELSE
0676                   WRITE(MSTU(11),*) ' ERROR WITH DELM ',DELM,XMDIF
0677                   WRITE(MSTU(11),*) ' KF = ',KF
0678                   WRITE(MSTU(11),*) ' KF(decay) = ',(KFDP(IDC,J),J=1,3)
0679                 ENDIF
0680               ENDIF
0681               GOTO 360
0682             ELSEIF(IDC.EQ.IDCSV) THEN
0683               WRITE(MSTU(11),*) ' Error in PYMSIN: SUSY decay ',
0684      &        'channel not recognized:'
0685               WRITE(MSTU(11),*) KF,' -> ',(IDLAM(IL,J),J=1,3)
0686               GOTO 360
0687             ELSE
0688               GOTO 350
0689             ENDIF
0690   360     CONTINUE
0691
0692 C...Store width, cutoff and lifetime.
0693           PMAS(KC,2)=XLAM(0)
0694           IF(PMAS(KC,2).LT.0.1D0*DELM) THEN
0695             PMAS(KC,3)=PMAS(KC,2)*10D0
0696           ELSE
0697             PMAS(KC,3)=0.95D0*DELM
0698           ENDIF
0699           IF(PMAS(KC,2).NE.0D0) THEN
0700             PMAS(KC,4)=PARU(3)/PMAS(KC,2)*1D-12
0701           ENDIF
0702 C...Write decays to SLHA file
0703           IF (IMSS(24).NE.0) THEN
0704             IFAIL=0
0705             CALL PYSLHA(4,KF,IFAIL)
0706           ENDIF
0707
0708         ENDIF
0709   370 CONTINUE
0710
0711       RETURN
0712       END