ISAJET/code/whiggs.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE WHIGGS
   3 C
   4 C          Finish generation of whiggs events started bY TWOJET.
   5 C          Select W decay modes as allowed by WMODE1, WMODE2.
   6 C          Generate W decay angles and put vectors in PPAIR.
   7 C
   8 C
   9 #if defined(CERNLIB_IMPNONE)
  10       IMPLICIT NONE
  11 #endif
  12 #include "isajet/itapes.inc"
  13 #include "isajet/qcdpar.inc"
  14 #include "isajet/jetpar.inc"
  15 #include "isajet/primar.inc"
  16 #include "isajet/q1q2.inc"
  17 #include "isajet/jetsig.inc"
  18 #include "isajet/const.inc"
  19 #include "isajet/qsave.inc"
  20 #include "isajet/wcon.inc"
  21 #include "isajet/pjets.inc"
  22 #include "isajet/pinits.inc"
  23 #include "isajet/keys.inc"
  24 #include "isajet/hcon.inc"
  25 #include "isajet/wwpar.inc"
  26 #include "isajet/xmssm.inc"
  27 C
  28       DIMENSION X(2),LIST(25),P1WCM(4),P2WCM(4),P1LAB(4),P2LAB(4)
  29      1,PBOOST(4)
  30       EQUIVALENCE (X(1),X1)
  31       DIMENSION JWWTYP(2)
  32       REAL GVQ(2),GAQ(2),GVL(2),GAL(2)
  33       REAL X,RND,RANF,CBRWW,AMASS,AM0,AM1,AM2,
  34      $E1CM,E2CM,P12CM,CTHCM,STHCM,PHICM,CPHICM,SPHICM,P1WCM,P2WCM,
  35      $PBOOST,P1LAB,P2LAB,ZHSIG,ZHMAX
  36       INTEGER JWWTYP,JET,JWT,JQ,IQ1,IQ2,LIST,NREJ,NJ0,K
  37       REAL BRANCH(12),SUMBR,BETAWH,GAMWH,PZWHCM,CTHD,WHSIG
  38       INTEGER IDABS,IDABS1,IDABSJ,IDIABS
  39 C
  40       DATA LIST/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
  41      $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16/
  42 C
  43       GVQ(1)=.25-2*SIN2W/3.
  44       GVQ(2)=-.25+SIN2W/3.
  45       GAQ(1)=-.25
  46       GAQ(2)=.25
  47       GVL(1)=.25
  48       GVL(2)=-.25+SIN2W
  49       GAL(1)=-.25
  50       GAL(2)=.25
  51       NPAIR=0
  52       IF(KEYS(10).AND..NOT.GOMSSM) THEN
  53         JWWTYP(1)=JETTYP(1)-25
  54         JWWTYP(2)=JETTYP(2)-25
  55       ELSEIF(KEYS(10).AND.GOMSSM) THEN
  56         JWWTYP(1)=JETTYP(1)-76
  57         JWWTYP(2)=JETTYP(2)-76
  58       ENDIF
  59 C
  60 C          Select W decay modes and put in /JETSET/. First particle
  61 C          is always the fermion.
  62
  63       DO 200 JET=1,2
  64         IDABS=IABS(IDJETS(JET))
  65         IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 200
  66         RND=RANF()
  67         JWT=JWWTYP(JET)
  68 C          Must only consider allowed decays for this mass
  69         SUMBR=0.
  70         DO 201 JQ=1,12
  71           IQ1=2*JQ
  72           IQ2=MATCH(IQ1,JWT)
  73           IF(IQ2.EQ.0) THEN
  74             BRANCH(JQ)=0.
  75             GO TO 201
  76           ENDIF
  77           AM1=AMASS(LIST(IQ1))
  78           AM2=AMASS(LIST(IQ2))
  79           IF(AM1+AM2.LT.PJETS(5,JET)) THEN
  80             BRANCH(JQ)=RBRWW(JQ,JWT,JET)
  81             SUMBR=SUMBR+BRANCH(JQ)
  82           ELSE
  83             BRANCH(JQ)=0.
  84           ENDIF
  85 201     CONTINUE
  86         IF(SUMBR.LE.0.) GO TO 998
  87         DO 202 JQ=1,12
  88 202     BRANCH(JQ)=BRANCH(JQ)/SUMBR
  89 C
  90         CBRWW=0.
  91         DO 210 JQ=1,12
  92           CBRWW=CBRWW+BRANCH(JQ)
  93           IF(RND.GT.CBRWW) GO TO 210
  94           IQ1=2*JQ
  95           IQ2=MATCH(IQ1,JWT)
  96           IDPAIR(NPAIR+1)=LIST(IQ1)
  97           IDPAIR(NPAIR+2)=LIST(IQ2)
  98           PPAIR(5,NPAIR+1)=AMASS(LIST(IQ1))
  99           PPAIR(5,NPAIR+2)=AMASS(LIST(IQ2))
 100           JPAIR(NPAIR+1)=JET
 101           JPAIR(NPAIR+2)=JET
 102           NPAIR=NPAIR+2
 103           JQWW(JET)=JQ
 104           GO TO 200
 105 210     CONTINUE
 106 200   CONTINUE
 107 C
 108 C          Generate decay uniformly in angle and put in PPAIR.
 109 C          Will check cross section later.
 110 C
 111       NREJ=0
 112 300   NJ0=2
 113       DO 310 JET=1,2
 114         IDABS=IABS(IDJETS(JET))
 115         IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 310
 116 C          Construct W com momenta.
 117         IDABSJ=IDABS
 118         AM0=PJETS(5,JET)
 119         AM1=PPAIR(5,NJ0-1)
 120         AM2=PPAIR(5,NJ0)
 121         E1CM=(AM0**2+AM1**2-AM2**2)/(2.*AM0)
 122         E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)
 123         P12CM=(AM0**2-AM1**2-AM2**2)**2-4.*(AM1*AM2)**2
 124         P12CM=SQRT(P12CM)/(2.*AM0)
 125         CTHCM=2.*RANF()-1.
 126         STHCM=SQRT(1.-CTHCM**2)
 127         PHICM=2.*PI*RANF()
 128         CPHICM=COS(PHICM)
 129         SPHICM=SIN(PHICM)
 130         P1WCM(1)=P12CM*STHCM*CPHICM
 131         P2WCM(1)=-P1WCM(1)
 132         P1WCM(2)=P12CM*STHCM*SPHICM
 133         P2WCM(2)=-P1WCM(2)
 134         P1WCM(3)=P12CM*CTHCM
 135         P2WCM(3)=-P1WCM(3)
 136         P1WCM(4)=E1CM
 137         P2WCM(4)=E2CM
 138 C          Boost to lab frame.
 139         DO 320 K=1,3
 140 320     PBOOST(K)=-PJETS(K,JET)
 141         PBOOST(4)=PJETS(4,JET)
 142         CALL LBOOST(PBOOST,1,P1WCM,P1LAB)
 143         CALL LBOOST(PBOOST,1,P2WCM,P2LAB)
 144         DO 330 K=1,4
 145           PPAIR(K,NJ0-1)=P1LAB(K)
 146           PPAIR(K,NJ0)=P2LAB(K)
 147 330     CONTINUE
 148       NJ0=NJ0+2
 149 310   CONTINUE
 150 C
 151 C          Impose simple (1+-cos(theta))**2 decay distribution for WH
 152 C          Must use P1 in WH CoM frame
 153       IF (IDABSJ.NE.80.AND.IDABSJ.NE.90) GO TO 400
 154       BETAWH=(PJETS(3,1)+PJETS(3,2))/(PJETS(4,1)+PJETS(4,2))
 155       GAMWH=1./SQRT(1.-BETAWH**2)
 156       PZWHCM=GAMWH*(P1LAB(3)-BETAWH*P1LAB(4))
 157       CTHD=PZWHCM/SQRT(P1LAB(1)**2+P1LAB(2)**2+PZWHCM**2)
 158       IF (IDINIT(1).LT.0) CTHD=-CTHD
 159       IDIABS=IABS(IDINIT(1))
 160       IDABS1=IABS(IDPAIR(1))
 161       IF (IDABSJ.EQ.80) THEN
 162         WHSIG=(1.+CTHD)**2
 163         IF(WHSIG.GT.4*RANF()) GO TO 400
 164       END IF
 165       IF (IDABSJ.EQ.90) THEN
 166         IF (IDIABS.EQ.1.OR.IDIABS.EQ.4) THEN
 167           IF (IDABS1.EQ.1.OR.IDABS1.EQ.4) THEN
 168             ZHSIG=(GVQ(1)**2+GAQ(1)**2)**2*(1.+CTHD**2)
 169      $             +8*GVQ(1)*GAQ(1)*GVQ(1)*GAQ(1)*CTHD
 170             ZHMAX=2*(GVQ(1)**2+GAQ(1)**2)**2
 171      $             +8*GVQ(1)*GAQ(1)*GVQ(1)*GAQ(1)
 172           ELSEIF (IDABS1.EQ.2.OR.IDABS1.EQ.3.OR.IDABS1.EQ.5) THEN
 173             ZHSIG=(GVQ(1)**2+GAQ(1))*(GVQ(2)**2+GAQ(2)**2)*(1.+CTHD**2)
 174      $             +8*GVQ(1)*GAQ(1)*GVQ(2)*GAQ(2)*CTHD
 175             ZHMAX=(GVQ(1)**2+GAQ(1))*(GVQ(2)**2+GAQ(2)**2)*2
 176      $             +8*GVQ(1)*GAQ(1)*GVQ(2)*GAQ(2)
 177           ELSEIF (IDABS1.EQ.11.OR.IDABS1.EQ.13.OR.IDABS1.EQ.15) THEN
 178             ZHSIG=(GVQ(1)**2+GAQ(1))*(GVL(1)**2+GAL(1)**2)*(1.+CTHD**2)
 179      $             +8*GVQ(1)*GAQ(1)*GVL(1)*GAL(1)*CTHD
 180             ZHMAX=(GVQ(1)**2+GAQ(1))*(GVL(1)**2+GAL(1)**2)*2
 181      $             +8*GVQ(1)*GAQ(1)*GVL(1)*GAL(1)
 182           ELSEIF (IDABS1.EQ.12.OR.IDABS1.EQ.14.OR.IDABS1.EQ.16) THEN
 183             ZHSIG=(GVQ(1)**2+GAQ(1))*(GVL(2)**2+GAL(2)**2)*(1.+CTHD**2)
 184      $             +8*GVQ(1)*GAQ(1)*GVL(2)*GAL(2)*CTHD
 185             ZHMAX=(GVQ(1)**2+GAQ(1))*(GVL(2)**2+GAL(2)**2)*2
 186      $             +8*GVQ(1)*GAQ(1)*GVL(2)*GAL(2)
 187           END IF
 188         ELSE IF (IDIABS.EQ.2.OR.IDIABS.EQ.3.OR.IDIABS.EQ.5) THEN
 189           IF (IDABS1.EQ.1.OR.IDABS1.EQ.4) THEN
 190             ZHSIG=(GVQ(2)**2+GAQ(2)**2)**2*(1.+CTHD**2)
 191      $             +8*GVQ(2)*GAQ(2)*GVQ(1)*GAQ(1)*CTHD
 192             ZHMAX=(GVQ(2)**2+GAQ(2)**2)**2*2
 193      $             +8*GVQ(2)*GAQ(2)*GVQ(1)*GAQ(1)
 194           ELSEIF (IDABS1.EQ.2.OR.IDABS1.EQ.3.OR.IDABS1.EQ.5) THEN
 195             ZHSIG=(GVQ(2)**2+GAQ(2))*(GVQ(2)**2+GAQ(2)**2)*(1.+CTHD**2)
 196      $             +8*GVQ(2)*GAQ(2)*GVQ(2)*GAQ(2)*CTHD
 197             ZHMAX=(GVQ(2)**2+GAQ(2))*(GVQ(2)**2+GAQ(2)**2)*2
 198      $             +8*GVQ(2)*GAQ(2)*GVQ(2)*GAQ(2)
 199           ELSEIF (IDABS1.EQ.11.OR.IDABS1.EQ.13.OR.IDABS1.EQ.15) THEN
 200             ZHSIG=(GVQ(2)**2+GAQ(2))*(GVL(1)**2+GAL(1)**2)*(1.+CTHD**2)
 201      $             +8*GVQ(2)*GAQ(2)*GVL(1)*GAL(1)*CTHD
 202             ZHMAX=(GVQ(2)**2+GAQ(2))*(GVL(1)**2+GAL(1)**2)*2
 203      $             +8*GVQ(2)*GAQ(2)*GVL(1)*GAL(1)
 204           ELSEIF (IDABS1.EQ.12.OR.IDABS1.EQ.14.OR.IDABS1.EQ.16) THEN
 205             ZHSIG=(GVQ(2)**2+GAQ(2))*(GVL(2)**2+GAL(2)**2)*(1.+CTHD**2)
 206      $             +8*GVQ(2)*GAQ(2)*GVL(2)*GAL(2)*CTHD
 207             ZHMAX=(GVQ(2)**2+GAQ(2))*(GVL(2)**2+GAL(2)**2)*2
 208      $             +8*GVQ(2)*GAQ(2)*GVL(2)*GAL(2)
 209           END IF
 210         END IF
 211         IF(ZHSIG.GT.RANF()*ZHMAX) GO TO 400
 212       END IF
 213       NREJ=NREJ+1
 214       IF(NREJ.LT.NTRIES) GO TO 300
 215       GO TO 999
 216 C
 217 C          Good event
 218 C
 219 400   CONTINUE
 220       RETURN
 221 C
 222 999   CALL PRTEVT(0)
 223       WRITE(ITLIS,9991) NREJ
 224 9991  FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE ',
 225      1'A GOOD WHIGGS EVENT.'/' CHECK LIMITS OR INCREASE NTRIES.')
 226       STOP 99
 227 998   CALL PRTEVT(0)
 228       WRITE(ITLIS,9981) JET
 229 9981  FORMAT(//' ERROR IN WHIGGS ... NO DECAY POSSIBLE FOR JET',I3)
 230       STOP 99
 231       END