ISAJET/code/wpair.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE WPAIR
   3 C
   4 C          Finish generation of wpair 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          Also generate massless decay vectors PZERO for matrix
   9 C          element -- double precision for 32-bit machines.
  10 C
  11 C          Ver 6.26: Check kinematics for W -> ff decay, since Z0 from
  12 C                    Higgs decay can be virtual.
  13 C          Ver. 6.30: Added check in loop 201.
  14 C          Ver. 7.14: Add MSSM Higgs hooks
  15 C
  16 #if defined(CERNLIB_IMPNONE)
  17       IMPLICIT NONE
  18 #endif
  19 #include "isajet/itapes.inc"
  20 #include "isajet/qcdpar.inc"
  21 #include "isajet/jetpar.inc"
  22 #include "isajet/primar.inc"
  23 #include "isajet/q1q2.inc"
  24 #include "isajet/jetsig.inc"
  25 #include "isajet/wwsig.inc"
  26 #include "isajet/wwpar.inc"
  27 #include "isajet/const.inc"
  28 #include "isajet/qsave.inc"
  29 #include "isajet/wcon.inc"
  30 #include "isajet/pjets.inc"
  31 #include "isajet/pinits.inc"
  32 #include "isajet/keys.inc"
  33 #include "isajet/wsig.inc"
  34 #include "isajet/hcon.inc"
  35 #include "isajet/xmssm.inc"
  36 C
  37       DIMENSION X(2),LIST(25),P1WCM(4),P2WCM(4),P1LAB(4),P2LAB(4)
  38      $,P1CM0(4),P2CM0(4),P1LAB0(4),P2LAB0(4)
  39      1,PBOOST(4)
  40       EQUIVALENCE (X(1),X1)
  41       DIMENSION PWW(5,2)
  42       EQUIVALENCE (PWW(1,1),P3WW(1))
  43       DIMENSION JWWTYP(2),THWFF(2),PHIWFF(2)
  44 #if defined(CERNLIB_SINGLE)
  45       REAL P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST
  46 #endif
  47 #if defined(CERNLIB_DOUBLE)
  48       DOUBLE PRECISION P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST
  49 #endif
  50       REAL AMWW1,AMWW2,X,STRUC,STRUCW,RND,RANF,CBRWW,AMASS,AM0,AM1,AM2,
  51      $E1CM,E2CM,P12CM,CTHCM,STHCM,PHICM,CPHICM,SPHICM,P1WCM,P2WCM,
  52      $PBOOST,P1LAB,P2LAB,AFX,SGWWMX,P1LAB0,P2LAB0,THWFF,PHIWFF
  53       INTEGER IH,IQ,JWWTYP,JET,JWT,JQ,IQ1,IQ2,LIST,NREJ,NJ0,K
  54       REAL BRANCH(12),SUMBR
  55       INTEGER IDABS,IDABS1,IDABS2
  56 C
  57       DATA LIST/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
  58      $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16/
  59 C
  60 C          Initialize for given W type.
  61       AMWW1=PJETS(5,1)
  62       AMWW2=PJETS(5,2)
  63       CALL WWKIN(AMWW1,AMWW2)
  64       NPAIR=0
  65 C
  66 C          Calculate and save structure functions.
  67       DO 120 IH=1,2
  68       DO 121 IQ=1,13
  69 121   QSAVE(IQ,IH)=STRUC(X(IH),QSQ,IQ,IDIN(IH))/X(IH)
  70       DO 122 IQ=14,25
  71 122   QSAVE(IQ,IH)=0.
  72       IF(KEYS(7).OR.KEYS(9)) THEN
  73         DO 123 IQ=26,29
  74 123     QSAVE(IQ,IH)=STRUCW(X(IH),IQ-25,IDIN(IH))/X(IH)
  75       ENDIF
  76 120   CONTINUE
  77 C          JWWTYP points to W types 1,2,3,4
  78       IF(KEYS(6)) THEN
  79         JWWTYP(1)=JETTYP(1)
  80         JWWTYP(2)=JETTYP(2)
  81       ELSEIF((KEYS(7).AND..NOT.GOMSSM).OR.KEYS(9)) THEN
  82         JWWTYP(1)=JETTYP(1)-25
  83         JWWTYP(2)=JETTYP(2)-25
  84       ELSEIF(KEYS(7).AND.GOMSSM) THEN
  85         JWWTYP(1)=JETTYP(1)-76
  86         JWWTYP(2)=JETTYP(2)-76
  87       ENDIF
  88 C
  89 C          Select W decay modes and put in /JETSET/. First particle
  90 C          is always the fermion.
  91
  92       DO 200 JET=1,2
  93         IDABS=IABS(IDJETS(JET))
  94         IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 200
  95         RND=RANF()
  96         JWT=JWWTYP(JET)
  97 C          Must only consider allowed decays for this mass
  98         SUMBR=0.
  99         DO 201 JQ=1,12
 100           IQ1=2*JQ
 101           IQ2=MATCH(IQ1,JWT)
 102           IF(IQ2.EQ.0) THEN
 103             BRANCH(JQ)=0.
 104             GO TO 201
 105           ENDIF
 106           AM1=AMASS(LIST(IQ1))
 107           AM2=AMASS(LIST(IQ2))
 108           IF(AM1+AM2.LT.PJETS(5,JET)) THEN
 109             BRANCH(JQ)=RBRWW(JQ,JWT,JET)
 110             SUMBR=SUMBR+BRANCH(JQ)
 111           ELSE
 112             BRANCH(JQ)=0.
 113           ENDIF
 114 201     CONTINUE
 115         IF(SUMBR.LE.0.) GO TO 998
 116         DO 202 JQ=1,12
 117 202     BRANCH(JQ)=BRANCH(JQ)/SUMBR
 118 C
 119         CBRWW=0.
 120         DO 210 JQ=1,12
 121           CBRWW=CBRWW+BRANCH(JQ)
 122           IF(RND.GT.CBRWW) GO TO 210
 123           IQ1=2*JQ
 124           IQ2=MATCH(IQ1,JWT)
 125           IDPAIR(NPAIR+1)=LIST(IQ1)
 126           IDPAIR(NPAIR+2)=LIST(IQ2)
 127           PPAIR(5,NPAIR+1)=AMASS(LIST(IQ1))
 128           PPAIR(5,NPAIR+2)=AMASS(LIST(IQ2))
 129           JPAIR(NPAIR+1)=JET
 130           JPAIR(NPAIR+2)=JET
 131           NPAIR=NPAIR+2
 132           JQWW(JET)=JQ
 133           GO TO 200
 134 210     CONTINUE
 135 200   CONTINUE
 136 C
 137 C          Generate decay uniformly in angle and put in PPAIR.
 138 C          Will check cross section later.
 139 C
 140       NREJ=0
 141 300   NJ0=2
 142       DO 310 JET=1,2
 143         IDABS=IABS(IDJETS(JET))
 144         IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 310
 145 C          Construct W com momenta.
 146         AM0=PJETS(5,JET)
 147         AM1=PPAIR(5,NJ0-1)
 148         AM2=PPAIR(5,NJ0)
 149         E1CM=(AM0**2+AM1**2-AM2**2)/(2.*AM0)
 150         E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)
 151         P12CM=(AM0**2-AM1**2-AM2**2)**2-4.*(AM1*AM2)**2
 152         P12CM=SQRT(P12CM)/(2.*AM0)
 153         CTHCM=2.*RANF()-1.
 154         STHCM=SQRT(1.-CTHCM**2)
 155         PHICM=2.*PI*RANF()
 156         CPHICM=COS(PHICM)
 157         SPHICM=SIN(PHICM)
 158         P1WCM(1)=P12CM*STHCM*CPHICM
 159         P2WCM(1)=-P1WCM(1)
 160         P1WCM(2)=P12CM*STHCM*SPHICM
 161         P2WCM(2)=-P1WCM(2)
 162         P1WCM(3)=P12CM*CTHCM
 163         P2WCM(3)=-P1WCM(3)
 164         P1WCM(4)=E1CM
 165         P2WCM(4)=E2CM
 166 C          Also construct zero mass vectors at same angle
 167 #if defined(CERNLIB_SINGLE)
 168 C          Single precision.
 169         P1CM0(1)=.5*AM0*STHCM*CPHICM
 170         P2CM0(1)=-P1CM0(1)
 171         P1CM0(2)=.5*AM0*STHCM*SPHICM
 172         P2CM0(2)=-P1CM0(2)
 173         P1CM0(3)=.5*AM0*CTHCM
 174         P2CM0(3)=-P1CM0(3)
 175         P1CM0(4)=.5*AM0
 176         P2CM0(4)=P1CM0(4)
 177 #endif
 178 #if defined(CERNLIB_DOUBLE)
 179 C          Double precision.
 180         DAM0=AM0
 181         DCTH=CTHCM
 182         DSTH=DSQRT(1.D0-DCTH**2)
 183         DPHI=PHICM
 184         P1CM0(1)=.5*AM0*DSTH*DCOS(DPHI)
 185         P2CM0(1)=-P1CM0(1)
 186         P1CM0(2)=.5*AM0*DSTH*DSIN(DPHI)
 187         P2CM0(2)=-P1CM0(2)
 188         P1CM0(3)=.5*AM0*DCTH
 189         P2CM0(3)=-P1CM0(3)
 190         P1CM0(4)=.5*AM0
 191         P2CM0(4)=P1CM0(4)
 192 #endif
 193 C          Boost to lab frame.
 194         DO 320 K=1,3
 195 320     PBOOST(K)=-PJETS(K,JET)
 196         PBOOST(4)=PJETS(4,JET)
 197         CALL LBOOST(PBOOST,1,P1WCM,P1LAB)
 198         CALL LBOOST(PBOOST,1,P2WCM,P2LAB)
 199         DO 330 K=1,4
 200           PPAIR(K,NJ0-1)=P1LAB(K)
 201           PPAIR(K,NJ0)=P2LAB(K)
 202 330     CONTINUE
 203 C          Boost zero mass vectors -- double precision for 32 bits.
 204         PZERO(4,NJ0-1)=(P1CM0(4)*PWW(4,JET)+P1CM0(1)*PWW(1,JET)
 205      $  +P1CM0(2)*PWW(2,JET)+P1CM0(3)*PWW(3,JET))/PWW(5,JET)
 206         BOOST=(P1CM0(4)+PZERO(4,NJ0-1))/(PWW(4,JET)+PWW(5,JET))
 207         DO 340 K=1,3
 208 340     PZERO(K,NJ0-1)=P1CM0(K)+BOOST*PWW(K,JET)
 209         PZERO(4,NJ0)=(P2CM0(4)*PWW(4,JET)+P2CM0(1)*PWW(1,JET)
 210      $  +P2CM0(2)*PWW(2,JET)+P2CM0(3)*PWW(3,JET))/PWW(5,JET)
 211         BOOST=(P2CM0(4)+PZERO(4,NJ0))/(PWW(4,JET)+PWW(5,JET))
 212         DO 350 K=1,3
 213 350     PZERO(K,NJ0)=P2CM0(K)+BOOST*PWW(K,JET)
 214         NJ0=NJ0+2
 215 310   CONTINUE
 216 C
 217 C          Calculate cross section SIGWW2 containing TBRWW*RBRWW.
 218 C          Compare with WW cross section containing TBRWW. Ratio
 219 C          must be bounded by 3/(4*PI) for each W.
 220 C
 221       AFX=3./(2.*PI)
 222       IF(KEYS(6)) THEN
 223         CALL SIGWW2
 224         SGWWMX=SIGEVT
 225         IF(IDJETS(1).NE.10) SGWWMX=SGWWMX*RBRWW(JQWW(1),JWWTYP(1),1)*AFX
 226         IF(IDJETS(2).NE.10) SGWWMX=SGWWMX*RBRWW(JQWW(2),JWWTYP(2),2)*AFX
 227       ELSEIF(KEYS(7)) THEN
 228 C          Note that except for WW -> WW processes, SIGH3 just computes
 229 C          the decay angular distribution, so it can be used for both
 230 C          for SM and SUSY HL0/HH0 decays; HA0 -> WW is forbidden.
 231 C          For Z + HL0 decays, we just return, ie use phase space.
 232         IDABS1=IABS(IDJETS(1))
 233         IDABS2=IABS(IDJETS(2))
 234         IF(.NOT.(IDABS1.EQ.10.OR.IDABS1.EQ.80.OR.IDABS1.EQ.90)) RETURN
 235         IF(.NOT.(IDABS2.EQ.10.OR.IDABS2.EQ.80.OR.IDABS2.EQ.90)) RETURN
 236         CALL SIGH3
 237         SGWWMX=SIGLLQ*AFX**2
 238       ELSEIF(KEYS(9)) THEN
 239         CALL SIGTC3
 240         SGWWMX=SIGLLQ*AFX**2
 241       ENDIF
 242       IF(WWSIG.GT.SGWWMX*RANF()) GO TO 400
 243       NREJ=NREJ+1
 244       IF(NREJ.LT.NTRIES) GO TO 300
 245       GO TO 999
 246 C
 247 C          Good event
 248 C
 249 400   CONTINUE
 250       RETURN
 251 C
 252 999   CALL PRTEVT(0)
 253       WRITE(ITLIS,9991) NREJ
 254 9991  FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE ',
 255      1'A GOOD WPAIR EVENT.'/' CHECK LIMITS OR INCREASE NTRIES.')
 256       STOP 99
 257 998   CALL PRTEVT(0)
 258       WRITE(ITLIS,9981) JET
 259 9981  FORMAT(//' ERROR IN WPAIR ... NO DECAY POSSIBLE FOR JET',I3)
 260       STOP 99
 261       END