pythiaeRHIC/pythia/pyinre.f

0001
0002 C*********************************************************************
0003
0004 C...PYINRE
0005 C...Calculates full and effective widths of gauge bosons, stores
0006 C...masses and widths, rescales coefficients to be used for
0007 C...resonance production generation.
0008
0009       SUBROUTINE PYINRE
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Parameter statement to help give large particle numbers.
0016       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0017      &KEXCIT=4000000,KDIMEN=5000000)
0018 C...Commonblocks.
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0022       COMMON/PYDAT4/CHAF(500,2)
0023       CHARACTER CHAF*16
0024       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0025       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0026       COMMON/PYINT1/MINT(400),VINT(400)
0027       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0028       COMMON/PYINT4/MWID(500),WIDS(500,5)
0029       COMMON/PYINT6/PROC(0:500)
0030       CHARACTER PROC*28
0031       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
0032       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYSUBS/,/PYPARS/,
0033      &/PYINT1/,/PYINT2/,/PYINT4/,/PYINT6/,/PYMSSM/
0034 C...Local arrays and data.
0035       DIMENSION WDTP(0:400),WDTE(0:400,0:5),WDTPM(0:400),
0036      &WDTEM(0:400,0:5),KCORD(500),PMORD(500)
0037
0038 C...Born level couplings in MSSM Higgs doublet sector.
0039       XW=PARU(102)
0040       XWV=XW
0041       IF(MSTP(8).GE.2) XW=1D0-(PMAS(24,1)/PMAS(23,1))**2
0042       XW1=1D0-XW
0043       IF(MSTP(4).EQ.2) THEN
0044         TANBE=PARU(141)
0045         RATBE=((1D0-TANBE**2)/(1D0+TANBE**2))**2
0046         SQMZ=PMAS(23,1)**2
0047         SQMW=PMAS(24,1)**2
0048         SQMH=PMAS(25,1)**2
0049         SQMA=SQMH*(SQMZ-SQMH)/(SQMZ*RATBE-SQMH)
0050         SQMHP=0.5D0*(SQMA+SQMZ+SQRT((SQMA+SQMZ)**2-4D0*SQMA*SQMZ*RATBE))
0051         SQMHC=SQMA+SQMW
0052         IF(SQMH.GE.SQMZ.OR.MIN(SQMA,SQMHP,SQMHC).LE.0D0) THEN
0053           WRITE(MSTU(11),5000)
0054           CALL PYSTOP(101)
0055         ENDIF
0056         PMAS(35,1)=SQRT(SQMHP)
0057         PMAS(36,1)=SQRT(SQMA)
0058         PMAS(37,1)=SQRT(SQMHC)
0059         ALSU=0.5D0*ATAN(2D0*TANBE*(SQMA+SQMZ)/((1D0-TANBE**2)*
0060      &  (SQMA-SQMZ)))
0061         BESU=ATAN(TANBE)
0062         PARU(142)=1D0
0063         PARU(143)=1D0
0064         PARU(161)=-SIN(ALSU)/COS(BESU)
0065         PARU(162)=COS(ALSU)/SIN(BESU)
0066         PARU(163)=PARU(161)
0067         PARU(164)=SIN(BESU-ALSU)
0068         PARU(165)=PARU(164)
0069         PARU(168)=SIN(BESU-ALSU)+0.5D0*COS(2D0*BESU)*SIN(BESU+ALSU)/XW
0070         PARU(171)=COS(ALSU)/COS(BESU)
0071         PARU(172)=SIN(ALSU)/SIN(BESU)
0072         PARU(173)=PARU(171)
0073         PARU(174)=COS(BESU-ALSU)
0074         PARU(175)=PARU(174)
0075         PARU(176)=COS(2D0*ALSU)*COS(BESU+ALSU)-2D0*SIN(2D0*ALSU)*
0076      &  SIN(BESU+ALSU)
0077         PARU(177)=COS(2D0*BESU)*COS(BESU+ALSU)
0078         PARU(178)=COS(BESU-ALSU)-0.5D0*COS(2D0*BESU)*COS(BESU+ALSU)/XW
0079         PARU(181)=TANBE
0080         PARU(182)=1D0/TANBE
0081         PARU(183)=PARU(181)
0082         PARU(184)=0D0
0083         PARU(185)=PARU(184)
0084         PARU(186)=COS(BESU-ALSU)
0085         PARU(187)=SIN(BESU-ALSU)
0086         PARU(188)=PARU(186)
0087         PARU(189)=PARU(187)
0088         PARU(190)=0D0
0089         PARU(195)=COS(BESU-ALSU)
0090       ENDIF
0091
0092 C...Reset effective widths of gauge bosons.
0093       DO 110 I=1,500
0094         DO 100 J=1,5
0095           WIDS(I,J)=1D0
0096   100   CONTINUE
0097   110 CONTINUE
0098
0099 C...Order resonances by increasing mass (except Z0 and W+/-).
0100       NRES=0
0101       DO 140 KC=1,500
0102         KF=KCHG(KC,4)
0103         IF(KF.EQ.0) GOTO 140
0104         IF(MWID(KC).EQ.0) GOTO 140
0105         IF(KC.EQ.7.OR.KC.EQ.8.OR.KC.EQ.17.OR.KC.EQ.18) THEN
0106           IF(MSTP(1).LE.3) GOTO 140
0107         ENDIF
0108         IF(KF/KSUSY1.EQ.1.OR.KF/KSUSY1.EQ.2) THEN
0109           IF(IMSS(1).LE.0) GOTO 140
0110         ENDIF
0111         NRES=NRES+1
0112         PMRES=PMAS(KC,1)
0113         IF(KC.EQ.23.OR.KC.EQ.24) PMRES=0D0
0114         DO 120 I1=NRES-1,1,-1
0115           IF(PMRES.GE.PMORD(I1)) GOTO 130
0116           KCORD(I1+1)=KCORD(I1)
0117           PMORD(I1+1)=PMORD(I1)
0118   120   CONTINUE
0119   130   KCORD(I1+1)=KC
0120         PMORD(I1+1)=PMRES
0121   140 CONTINUE
0122
0123 C...Loop over possible resonances.
0124       DO 180 I=1,NRES
0125         KC=KCORD(I)
0126         KF=KCHG(KC,4)
0127
0128 C...Check that no fourth generation channels on by mistake.
0129         IF(MSTP(1).LE.3) THEN
0130           DO 150 J=1,MDCY(KC,3)
0131             IDC=J+MDCY(KC,2)-1
0132             KFA1=IABS(KFDP(IDC,1))
0133             KFA2=IABS(KFDP(IDC,2))
0134             IF(KFA1.EQ.7.OR.KFA1.EQ.8.OR.KFA1.EQ.17.OR.KFA1.EQ.18.OR.
0135      &      KFA2.EQ.7.OR.KFA2.EQ.8.OR.KFA2.EQ.17.OR.KFA2.EQ.18)
0136      &      MDME(IDC,1)=-1
0137   150     CONTINUE
0138         ENDIF
0139
0140 C...Check that no supersymmetric channels on by mistake.
0141         IF(IMSS(1).LE.0) THEN
0142           DO 160 J=1,MDCY(KC,3)
0143             IDC=J+MDCY(KC,2)-1
0144             KFA1S=IABS(KFDP(IDC,1))/KSUSY1
0145             KFA2S=IABS(KFDP(IDC,2))/KSUSY1
0146             IF(KFA1S.EQ.1.OR.KFA1S.EQ.2.OR.KFA2S.EQ.1.OR.KFA2S.EQ.2)
0147      &      MDME(IDC,1)=-1
0148   160     CONTINUE
0149         ENDIF
0150
0151 C...Find mass and evaluate width.
0152         PMR=PMAS(KC,1)
0153         IF(KF.EQ.25.OR.KF.EQ.35.OR.KF.EQ.36) MINT(62)=1
0154         IF(MWID(KC).EQ.3) MINT(63)=1
0155         CALL PYWIDT(KF,PMR**2,WDTP,WDTE)
0156         MINT(51)=0
0157
0158 C...Evaluate suppression factors due to non-simulated channels.
0159         IF(KCHG(KC,3).EQ.0) THEN
0160           WDTP0I=0D0
0161           IF(WDTP(0).GT.0D0) WDTP0I=1D0/WDTP(0)
0162           WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))**2+
0163      &    2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
0164      &    2D0*WDTE(0,4)*WDTE(0,5))*WDTP0I**2
0165           WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))*WDTP0I
0166           WIDS(KC,3)=0D0
0167           WIDS(KC,4)=0D0
0168           WIDS(KC,5)=0D0
0169         ELSE
0170           IF(MWID(KC).EQ.3) MINT(63)=1
0171           CALL PYWIDT(-KF,PMR**2,WDTPM,WDTEM)
0172           MINT(51)=0
0173           WDTP0I=0D0
0174           IF(WDTP(0).GT.0D0) WDTP0I=1D0/WDTP(0)
0175           WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))*(WDTEM(0,1)+WDTEM(0,3))+
0176      &    (WDTE(0,1)+WDTE(0,2))*(WDTEM(0,4)+WDTEM(0,5))+
0177      &    (WDTE(0,4)+WDTE(0,5))*(WDTEM(0,1)+WDTEM(0,3))+
0178      &    WDTE(0,4)*WDTEM(0,5)+WDTE(0,5)*WDTEM(0,4))*WDTP0I**2
0179           WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))*WDTP0I
0180           WIDS(KC,3)=(WDTEM(0,1)+WDTEM(0,3)+WDTEM(0,4))*WDTP0I
0181           WIDS(KC,4)=((WDTE(0,1)+WDTE(0,2))**2+
0182      &    2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
0183      &    2D0*WDTE(0,4)*WDTE(0,5))*WDTP0I**2
0184           WIDS(KC,5)=((WDTEM(0,1)+WDTEM(0,3))**2+
0185      &    2D0*(WDTEM(0,1)+WDTEM(0,3))*(WDTEM(0,4)+WDTEM(0,5))+
0186      &    2D0*WDTEM(0,4)*WDTEM(0,5))*WDTP0I**2
0187         ENDIF
0188
0189 C...Set resonance widths and branching ratios;
0190 C...also on/off switch for decays.
0191         IF(MWID(KC).EQ.1.OR.MWID(KC).EQ.3) THEN
0192           PMAS(KC,2)=WDTP(0)
0193           PMAS(KC,3)=MIN(0.9D0*PMAS(KC,1),10D0*PMAS(KC,2))
0194           IF(MSTP(41).EQ.0.OR.MSTP(41).EQ.1) MDCY(KC,1)=MSTP(41)
0195           DO 170 J=1,MDCY(KC,3)
0196             IDC=J+MDCY(KC,2)-1
0197             BRAT(IDC)=0D0
0198             IF(WDTP(0).GT.0D0) BRAT(IDC)=WDTP(J)/WDTP(0)
0199   170     CONTINUE
0200         ENDIF
0201   180 CONTINUE
0202
0203 C...Flavours of leptoquark: redefine charge and name.
0204       KFLQQ=KFDP(MDCY(42,2),1)
0205       KFLQL=KFDP(MDCY(42,2),2)
0206       KCHG(42,1)=KCHG(PYCOMP(KFLQQ),1)*ISIGN(1,KFLQQ)+
0207      &KCHG(PYCOMP(KFLQL),1)*ISIGN(1,KFLQL)
0208       LL=1
0209       IF(IABS(KFLQL).EQ.13) LL=2
0210       IF(IABS(KFLQL).EQ.15) LL=3
0211       CHAF(42,1)='LQ_'//CHAF(IABS(KFLQQ),1)(1:1)//
0212      &CHAF(IABS(KFLQL),1)(1:LL)//' '
0213       CHAF(42,2)=CHAF(42,2)(1:4+LL)//'bar '
0214
0215 C...Special cases in treatment of gamma*/Z0: redefine process name.
0216       IF(MSTP(43).EQ.1) THEN
0217         PROC(1)='f + fbar -> gamma*'
0218         PROC(15)='f + fbar -> g + gamma*'
0219         PROC(19)='f + fbar -> gamma + gamma*'
0220         PROC(30)='f + g -> f + gamma*'
0221         PROC(35)='f + gamma -> f + gamma*'
0222       ELSEIF(MSTP(43).EQ.2) THEN
0223         PROC(1)='f + fbar -> Z0'
0224         PROC(15)='f + fbar -> g + Z0'
0225         PROC(19)='f + fbar -> gamma + Z0'
0226         PROC(30)='f + g -> f + Z0'
0227         PROC(35)='f + gamma -> f + Z0'
0228       ELSEIF(MSTP(43).EQ.3) THEN
0229         PROC(1)='f + fbar -> gamma*/Z0'
0230         PROC(15)='f + fbar -> g + gamma*/Z0'
0231         PROC(19)='f+ fbar -> gamma + gamma*/Z0'
0232         PROC(30)='f + g -> f + gamma*/Z0'
0233         PROC(35)='f + gamma -> f + gamma*/Z0'
0234       ENDIF
0235
0236 C...Special cases in treatment of gamma*/Z0/Z'0: redefine process name.
0237       IF(MSTP(44).EQ.1) THEN
0238         PROC(141)='f + fbar -> gamma*'
0239       ELSEIF(MSTP(44).EQ.2) THEN
0240         PROC(141)='f + fbar -> Z0'
0241       ELSEIF(MSTP(44).EQ.3) THEN
0242         PROC(141)='f + fbar -> Z''0'
0243       ELSEIF(MSTP(44).EQ.4) THEN
0244         PROC(141)='f + fbar -> gamma*/Z0'
0245       ELSEIF(MSTP(44).EQ.5) THEN
0246         PROC(141)='f + fbar -> gamma*/Z''0'
0247       ELSEIF(MSTP(44).EQ.6) THEN
0248         PROC(141)='f + fbar -> Z0/Z''0'
0249       ELSEIF(MSTP(44).EQ.7) THEN
0250         PROC(141)='f + fbar -> gamma*/Z0/Z''0'
0251       ENDIF
0252
0253 C...Special cases in treatment of WW -> WW: redefine process name.
0254       IF(MSTP(45).EQ.1) THEN
0255         PROC(77)='W+ + W+ -> W+ + W+'
0256       ELSEIF(MSTP(45).EQ.2) THEN
0257         PROC(77)='W+ + W- -> W+ + W-'
0258       ELSEIF(MSTP(45).EQ.3) THEN
0259         PROC(77)='W+/- + W+/- -> W+/- + W+/-'
0260       ENDIF
0261
0262 C...Format for error information.
0263  5000 FORMAT(1X,'Error: unphysical input tan^2(beta) and m_H ',
0264      &'combination'/1X,'Execution stopped!')
0265
0266       RETURN
0267       END