PYTHIA/pythia/pyinre.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYINRE
   5
   6 C...Calculates full and effective widths of gauge bosons, stores
   7 C...masses and widths, rescales coefficients to be used for
   8 C...resonance production generation.
   9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  10       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  11       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  12       COMMON/LUDAT4/CHAF(500)
  13       CHARACTER CHAF*8
  14       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  15       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
  16       COMMON/PYINT1/MINT(400),VINT(400)
  17       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
  18       COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
  19       COMMON/PYINT6/PROC(0:200)
  20       CHARACTER PROC*28
  21       SAVE /LUDAT1/,/LUDAT2/,/LUDAT3/,/LUDAT4/
  22       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT4/,/PYINT6/
  23       DIMENSION WDTP(0:40),WDTE(0:40,0:5),WDTPM(0:40),WDTEM(0:40,0:5)
  24       DIMENSION KCINP(16),KCORD(16),PMORD(16)
  25       DATA KCINP/23,24,25,6,7,8,17,18,32,34,35,36,37,38,39,40/
  26
  27 C...Born level couplings in MSSM Higgs doublet sector.
  28       XW=PARU(102)
  29       XWV=XW
  30       IF(MSTP(8).GE.2) XW=1.-(PMAS(24,1)/PMAS(23,1))**2
  31       XW1=1.-XW
  32       IF(MSTP(4).EQ.2) THEN
  33         TANBE=PARU(141)
  34         RATBE=((1.-TANBE**2)/(1.+TANBE**2))**2
  35         SQMZ=PMAS(23,1)**2
  36         SQMW=PMAS(24,1)**2
  37         SQMH=PMAS(25,1)**2
  38         SQMA=SQMH*(SQMZ-SQMH)/(SQMZ*RATBE-SQMH)
  39         SQMHP=0.5*(SQMA+SQMZ+SQRT((SQMA+SQMZ)**2-4.*SQMA*SQMZ*RATBE))
  40         SQMHC=SQMA+SQMW
  41         IF(SQMH.GE.SQMZ.OR.MIN(SQMA,SQMHP,SQMHC).LE.0.) THEN
  42           WRITE(MSTU(11),5000)
  43           STOP
  44         ENDIF
  45         PMAS(35,1)=SQRT(SQMHP)
  46         PMAS(36,1)=SQRT(SQMA)
  47         PMAS(37,1)=SQRT(SQMHC)
  48         ALSU=0.5*ATAN(2.*TANBE*(SQMA+SQMZ)/((1.-TANBE**2)*
  49      &  (SQMA-SQMZ)))
  50         BESU=ATAN(TANBE)
  51         PARU(142)=1.
  52         PARU(143)=1.
  53         PARU(161)=-SIN(ALSU)/COS(BESU)
  54         PARU(162)=COS(ALSU)/SIN(BESU)
  55         PARU(163)=PARU(161)
  56         PARU(164)=SIN(BESU-ALSU)
  57         PARU(165)=PARU(164)
  58         PARU(168)=SIN(BESU-ALSU)+0.5*COS(2.*BESU)*SIN(BESU+ALSU)/XW
  59         PARU(171)=COS(ALSU)/COS(BESU)
  60         PARU(172)=SIN(ALSU)/SIN(BESU)
  61         PARU(173)=PARU(171)
  62         PARU(174)=COS(BESU-ALSU)
  63         PARU(175)=PARU(174)
  64         PARU(176)=COS(2.*ALSU)*COS(BESU+ALSU)-2.*SIN(2.*ALSU)*
  65      &  SIN(BESU+ALSU)
  66         PARU(177)=COS(2.*BESU)*COS(BESU+ALSU)
  67         PARU(178)=COS(BESU-ALSU)-0.5*COS(2.*BESU)*COS(BESU+ALSU)/XW
  68         PARU(181)=TANBE
  69         PARU(182)=1./TANBE
  70         PARU(183)=PARU(181)
  71         PARU(184)=0.
  72         PARU(185)=PARU(184)
  73         PARU(186)=COS(BESU-ALSU)
  74         PARU(187)=SIN(BESU-ALSU)
  75         PARU(188)=PARU(186)
  76         PARU(189)=PARU(187)
  77         PARU(190)=0.
  78         PARU(195)=COS(BESU-ALSU)
  79       ENDIF
  80
  81 C...Change matrix element codes when top and 4th generation
  82 C...decay before fragmentation.
  83       IF(MSTP(48).GE.1) THEN
  84         IOFF=MDCY(6,2)-1
  85         DO 100 I=4,7
  86         MDME(IOFF+I,2)=0
  87   100   CONTINUE
  88         MDME(IOFF+9,2)=0
  89       ENDIF
  90       IF(MSTP(6).EQ.1) THEN
  91         IOFF=MDCY(7,2)-1
  92         DO 110 I=1,4
  93         MDME(IOFF+I,2)=0
  94   110   CONTINUE
  95         IOFF=MDCY(8,2)-1
  96         DO 120 I=1,4
  97         MDME(IOFF+I,2)=0
  98   120   CONTINUE
  99         IOFF=MDCY(17,2)-1
 100         MDME(IOFF+2,2)=0
 101         MDME(IOFF+3,2)=0
 102         MDME(IOFF+4,2)=0
 103         IOFF=MDCY(18,2)-1
 104         MDME(IOFF+1,2)=0
 105         MDME(IOFF+2,2)=0
 106       ELSEIF(MSTP(49).GE.1) THEN
 107         IOFF=MDCY(7,2)-1
 108         DO 130 I=4,7
 109         MDME(IOFF+I,2)=0
 110   130   CONTINUE
 111         MDME(IOFF+9,2)=0
 112         MDME(IOFF+10,2)=0
 113         IOFF=MDCY(8,2)-1
 114         DO 140 I=4,7
 115         MDME(IOFF+I,2)=0
 116   140   CONTINUE
 117         MDME(IOFF+9,2)=0
 118         MDME(IOFF+10,2)=0
 119         IOFF=MDCY(17,2)-1
 120         MDME(IOFF+4,2)=0
 121         MDME(IOFF+6,2)=0
 122         IOFF=MDCY(18,2)-1
 123         MDME(IOFF+2,2)=0
 124         MDME(IOFF+3,2)=0
 125       ENDIF
 126
 127 C...Reset full and effective widths of gauge bosons.
 128       DO 160 I=21,40
 129       DO 150 J=0,40
 130       WIDP(I,J)=0.
 131       WIDE(I,J)=0.
 132   150 CONTINUE
 133       WIDS(I,1)=1.
 134       WIDS(I,2)=1.
 135       WIDS(I,3)=1.
 136   160 CONTINUE
 137
 138 C...Order resonances by increasing mass (except Z0 and W+/-).
 139       DO 170 I=1,3
 140       KCORD(I)=KCINP(I)
 141       PMORD(I)=PMAS(KCORD(I),1)
 142   170 CONTINUE
 143       DO 200 I=4,16
 144       KCIN=KCINP(I)
 145       PMIN=PMAS(KCIN,1)
 146       DO 180 I1=I-1,3,-1
 147       IF(PMIN.GE.PMORD(I1)) GOTO 190
 148       KCORD(I1+1)=KCORD(I1)
 149       PMORD(I1+1)=PMORD(I1)
 150   180 CONTINUE
 151   190 KCORD(I1+1)=KCIN
 152       PMORD(I1+1)=PMIN
 153   200 CONTINUE
 154
 155 C...Loop over possible resonances.
 156       DO 250 I=1,16
 157       KC=KCORD(I)
 158       IF(KC.EQ.6.AND.MSTP(48).LE.0) GOTO 250
 159       IF(KC.EQ.7.OR.KC.EQ.8.OR.KC.EQ.17.OR.KC.EQ.18) THEN
 160         IF(MSTP(6).NE.1.AND.(MSTP(49).LE.0.OR.MSTP(1).LE.3)) GOTO 250
 161         IF(KC.EQ.18.AND.PMORD(I).LT.1.) GOTO 250
 162       ENDIF
 163       KCL=KC
 164       IF(KC.GE.6.AND.KC.LE.8) KCL=KC+20
 165       IF(KC.EQ.17.OR.KC.EQ.18) KCL=KC+12
 166
 167 C...Change decay modes for q* and l*.
 168       IF(MSTP(6).EQ.1.AND.(KC.EQ.7.OR.KC.EQ.8.OR.KC.EQ.17.OR.
 169      &KC.EQ.18)) THEN
 170         DO 210 J=1,MDCY(KC,3)
 171         IDC=J+MDCY(KC,2)-1
 172         KF2=KFDP(IDC,2)
 173         IF(KF2.EQ.7.OR.KF2.EQ.8.OR.KF2.EQ.17.OR.KF2.EQ.18)
 174      &  KFDP(IDC,2)=KF2-6
 175   210   CONTINUE
 176       ENDIF
 177
 178 C...Check that no fourth generation channels on by mistake.
 179       IF(MSTP(1).LE.3) THEN
 180         DO 220 J=1,MDCY(KC,3)
 181         IDC=J+MDCY(KC,2)-1
 182         KFA1=IABS(KFDP(IDC,1))
 183         KFA2=IABS(KFDP(IDC,2))
 184         IF(KFA1.EQ.7.OR.KFA1.EQ.8.OR.KFA1.EQ.17.OR.KFA1.EQ.18.OR.KFA2
 185      &  .EQ.7.OR.KFA2.EQ.8.OR.KFA2.EQ.17.OR.KFA2.EQ.18) MDME(IDC,1)=-1
 186   220   CONTINUE
 187       ENDIF
 188
 189 C...Find mass and evaluate width.
 190       PMR=PMAS(KC,1)
 191       IF(KC.EQ.25.OR.KC.EQ.35.OR.KC.EQ.36) MINT(62)=1
 192       CALL PYWIDT(KC,PMR**2,WDTP,WDTE)
 193       IF(KC.EQ.6.OR.KC.EQ.7.OR.KC.EQ.8.OR.KC.EQ.17.OR.KC.EQ.18)
 194      &CALL PYWIDT(-KC,PMR**2,WDTPM,WDTEM)
 195       MINT(51)=0
 196
 197 C...Evaluate suppression factors due to non-simulated channels.
 198       IF(KCHG(KC,3).EQ.0) THEN
 199         WIDS(KCL,1)=((WDTE(0,1)+WDTE(0,2))**2+
 200      &  2.*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
 201      &  2.*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
 202         WIDS(KCL,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
 203         WIDS(KCL,3)=0.
 204       ELSEIF(KC.EQ.6.OR.KC.EQ.7.OR.KC.EQ.8.OR.KC.EQ.17.OR.KC.EQ.18) THEN
 205         WIDS(KCL,1)=((WDTE(0,1)+WDTE(0,2))*(WDTEM(0,1)+WDTEM(0,3))+
 206      &  (WDTE(0,1)+WDTE(0,2))*(WDTEM(0,4)+WDTEM(0,5))+
 207      &  (WDTE(0,4)+WDTE(0,5))*(WDTEM(0,1)+WDTEM(0,3))+
 208      &  WDTE(0,4)*WDTEM(0,5)+WDTE(0,5)*WDTEM(0,4))/WDTP(0)**2
 209         WIDS(KCL,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
 210         WIDS(KCL,3)=(WDTEM(0,1)+WDTEM(0,3)+WDTEM(0,4))/WDTP(0)
 211       ELSE
 212         WIDS(KCL,1)=((WDTE(0,1)+WDTE(0,2))*(WDTE(0,1)+WDTE(0,3))+
 213      &  (WDTE(0,1)+WDTE(0,2)+WDTE(0,1)+WDTE(0,3))*(WDTE(0,4)+WDTE(0,5))+
 214      &  2.*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
 215         WIDS(KCL,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
 216         WIDS(KCL,3)=(WDTE(0,1)+WDTE(0,3)+WDTE(0,4))/WDTP(0)
 217         IF(KC.EQ.24) THEN
 218           VINT(91)=((WDTE(0,1)+WDTE(0,2))**2+2.*(WDTE(0,1)+WDTE(0,2))*
 219      &    (WDTE(0,4)+WDTE(0,5))+2.*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
 220           VINT(92)=((WDTE(0,1)+WDTE(0,3))**2+2.*(WDTE(0,1)+WDTE(0,3))*
 221      &    (WDTE(0,4)+WDTE(0,5))+2.*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
 222         ENDIF
 223       ENDIF
 224
 225 C...Find factors to give widths in GeV.
 226       AEM=ULALEM(PMR**2)
 227       IF(MSTP(8).GE.1) AEM=SQRT(2.)*PARU(105)*PMAS(24,1)**2*XW/PARU(1)
 228       IF(KC.LE.20) THEN
 229         FAC=PMR
 230       ELSEIF(KC.EQ.23.OR.KC.EQ.32) THEN
 231         FAC=AEM/(48.*XW*XW1)*PMR
 232       ELSEIF(KC.EQ.24.OR.KC.EQ.34) THEN
 233         FAC=AEM/(24.*XW)*PMR
 234       ELSEIF(KC.EQ.25.OR.KC.EQ.35.OR.KC.EQ.36.OR.KC.EQ.37) THEN
 235         FAC=AEM/(8.*XW)*(PMR/PMAS(24,1))**2*PMR
 236       ELSEIF(KC.EQ.38) THEN
 237         FAC=PMR
 238       ELSEIF(KC.EQ.39) THEN
 239         FAC=AEM/4.*PMR
 240       ELSEIF(KC.EQ.40) THEN
 241         FAC=AEM/(12.*XW)*PMR
 242       ENDIF
 243
 244 C...Translate widths into GeV and save them.
 245       DO 230 J=0,40
 246       WIDP(KCL,J)=FAC*WDTP(J)
 247       WIDE(KCL,J)=FAC*WDTE(J,0)
 248   230 CONTINUE
 249
 250 C...Set resonance widths and branching ratios in JETSET;
 251 C...also on/off switch for decays in PYTHIA/JETSET.
 252       PMAS(KC,2)=WIDP(KCL,0)
 253       PMAS(KC,3)=MIN(0.9*PMAS(KC,1),10.*PMAS(KC,2))
 254       MDCY(KC,1)=MSTP(41)
 255       DO 240 J=1,MDCY(KC,3)
 256       IDC=J+MDCY(KC,2)-1
 257       BRAT(IDC)=0.
 258       IF(WIDE(KCL,0).GT.0.) BRAT(IDC)=WIDE(KCL,J)/WIDE(KCL,0)
 259   240 CONTINUE
 260   250 CONTINUE
 261
 262 C...Flavours of leptoquark: redefine charge and name.
 263       KFLQQ=KFDP(MDCY(39,2),1)
 264       KFLQL=KFDP(MDCY(39,2),2)
 265       KCHG(39,1)=KCHG(IABS(KFLQQ),1)*ISIGN(1,KFLQQ)+
 266      &KCHG(IABS(KFLQL),1)*ISIGN(1,KFLQL)
 267       CHAF(39)(4:4)=CHAF(IABS(KFLQQ))(1:1)
 268       CHAF(39)(5:7)=CHAF(IABS(KFLQL))(1:3)
 269
 270 C...Scenario with q* and l*: redefine names.
 271       IF(MSTP(6).EQ.1) THEN
 272         CHAF(7)='d*'
 273         CHAF(8)='u*'
 274         CHAF(17)='e*'
 275         CHAF(18)='nu*_e'
 276       ENDIF
 277
 278 C...Special cases in treatment of gamma*/Z0: redefine process name.
 279       IF(MSTP(43).EQ.1) THEN
 280         PROC(1)='f + f~ -> gamma*'
 281         PROC(15)='f + f~ -> g + gamma*'
 282         PROC(19)='f + f~ -> gamma + gamma*'
 283         PROC(30)='f + g -> f + gamma*'
 284         PROC(35)='f + gamma -> f + gamma*'
 285       ELSEIF(MSTP(43).EQ.2) THEN
 286         PROC(1)='f + f~ -> Z0'
 287         PROC(15)='f + f~ -> g + Z0'
 288         PROC(19)='f + f~ -> gamma + Z0'
 289         PROC(30)='f + g -> f + Z0'
 290         PROC(35)='f + gamma -> f + Z0'
 291       ELSEIF(MSTP(43).EQ.3) THEN
 292         PROC(1)='f + f~ -> gamma*/Z0'
 293         PROC(15)='f + f~ -> g + gamma*/Z0'
 294         PROC(19)='f + f~ -> gamma + gamma*/Z0'
 295         PROC(30)='f + g -> f + gamma*/Z0'
 296         PROC(35)='f + gamma -> f + gamma*/Z0'
 297       ENDIF
 298
 299 C...Special cases in treatment of gamma*/Z0/Z'0: redefine process name.
 300       IF(MSTP(44).EQ.1) THEN
 301         PROC(141)='f + f~ -> gamma*'
 302       ELSEIF(MSTP(44).EQ.2) THEN
 303         PROC(141)='f + f~ -> Z0'
 304       ELSEIF(MSTP(44).EQ.3) THEN
 305         PROC(141)='f + f~ -> Z''0'
 306       ELSEIF(MSTP(44).EQ.4) THEN
 307         PROC(141)='f + f~ -> gamma*/Z0'
 308       ELSEIF(MSTP(44).EQ.5) THEN
 309         PROC(141)='f + f~ -> gamma*/Z''0'
 310       ELSEIF(MSTP(44).EQ.6) THEN
 311         PROC(141)='f + f~ -> Z0/Z''0'
 312       ELSEIF(MSTP(44).EQ.7) THEN
 313         PROC(141)='f + f~ -> gamma*/Z0/Z''0'
 314       ENDIF
 315
 316 C...Special cases in treatment of WW -> WW: redefine process name.
 317       IF(MSTP(45).EQ.1) THEN
 318         PROC(77)='W+ + W+ -> W+ + W+'
 319       ELSEIF(MSTP(45).EQ.2) THEN
 320         PROC(77)='W+ + W- -> W+ + W-'
 321       ELSEIF(MSTP(45).EQ.3) THEN
 322         PROC(77)='W+/- + W+/- -> W+/- + W+/-'
 323       ENDIF
 324
 325 C...Format for error information.
 326  5000 FORMAT(1X,'Error: unphysical input tan^2(beta) and m_H ',
 327      &'combination'/1X,'Execution stopped!')
 328
 329       RETURN
 330       END