ISAJET/code/sigww2.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE SIGWW2
   3 C
   4 C          Calculate WPAIR decay distribution
   5 C          D(SIGMA)/D(PT**2)D(Y1)D(Y2)D(OMEGA1)D(OMEGA2)
   6 C          for modes selected in WPAIR.
   7 C
   8 C          Also fix the initial parton types to those selected.
   9 C
  10 C          Cross sections from SCHOONSCHIP (1980) neglecting W width
  11 C          and quark masses. Hence use zero-mass vectors PZERO from
  12 C          WPAIR to define kinematics.
  13 C          QK(P1) + QB(P2) --> W1(P3) + W2(P4)
  14 C                   W1(P3) --> QK(Q1) + QB(Q2)
  15 C                   W2(P4) --> QK(Q3) + QB(Q4)
  16 C          S=(P3+P4)**2,  T=(P3-P1)**2,  U=(P3-P2)**2
  17 C          S1=(Q1+P4)**2, T1=(Q1-P1)**2, U1=(Q1-P2)**2
  18 C          S3=(Q3+P3)**2, T3=(Q3-P2)**2, U3=(Q3-P1)**2
  19 C          S13=(Q1+Q3)**2
  20 C          Note that the W+- final couplings have been set equal to 1.
  21 C          in the SCHOONSCHIP formulas and must be restored.
  22 C
  23 C          Need double precision for 32-bit machines.
  24 C
  25 C          Ver. 5.35 - correct symmetrization for DN DB -> W+ W-.
  26 C          Ver. 6.22 - use W + GM decay distributions from
  27 C                      Cortes, Hagiwara, and Herzog, NP B278, 26 (1986)
  28 C
  29 #if defined(CERNLIB_IMPNONE)
  30       IMPLICIT NONE
  31 #endif
  32 #include "isajet/itapes.inc"
  33 #include "isajet/qcdpar.inc"
  34 #include "isajet/jetpar.inc"
  35 #include "isajet/primar.inc"
  36 #include "isajet/q1q2.inc"
  37 #include "isajet/const.inc"
  38 #include "isajet/qsave.inc"
  39 #include "isajet/wcon.inc"
  40 #include "isajet/pjets.inc"
  41 #include "isajet/pinits.inc"
  42 #include "isajet/wwsig.inc"
  43 #include "isajet/wwpar.inc"
  44 C
  45       DIMENSION P1(5),P2(5),QSGN(6),PP1(4),PP2(4)
  46       EQUIVALENCE (S,SWW),(T,TWW),(U,UWW)
  47       EQUIVALENCE (P1(1),P1WW(1)),(P2(1),P2WW(1))
  48 C          Double precision kinematics for 32-bit.
  49 #if defined(CERNLIB_SINGLE)
  50       REAL S,T,U,T1,U1,T3,U3,P1,P2
  51      1,TX,UX,TT,UU,TT1,UU1,TT3,UU3,PP1,PP2
  52       REAL TERM,WWSS,WWST,WWTT,ZZALL,WZSS,WZST,WZSU,WZTU
  53      1,WGSS,WGST,WGSU,WGTU
  54 #endif
  55 #if defined(CERNLIB_DOUBLE)
  56       DOUBLE PRECISION S,T,U,T1,U1,T3,U3,P1,P2
  57      1,TX,UX,TT,UU,TT1,UU1,TT3,UU3,PP1,PP2
  58       DOUBLE PRECISION TERM,WWSS,WWST,WWTT,ZZALL,WZSS,WZST,WZSU,WZTU
  59      1,WGSS,WGST,WGSU,WGTU
  60 #endif
  61       REAL P3IS3,P3IS4,FJAC,AMW1,AMW2,GAM1,GAM2,SGN,QSGN,AMASS3
  62       REAL P1DQ2,P2DQ1
  63       REAL A1,B1,A2,B2,ES,SMS,SMSZG,EQ3(12)
  64       REAL Q(5),QB(5),KK(5),E(5),EB(5)
  65       INTEGER K,JQ1,JQ3,JW1,JW2,IW1,IW2,IQ1,IQ2,IQ,ISWAPQ,JW,JZ,ISGN
  66       INTEGER IFLI,IFLJ,JG,IL,IW
  67       LOGICAL LQK1
  68 C
  69       DATA QSGN/1.,-1.,-1.,1.,-1.,1./
  70       DATA EQ3/2.,-1.,-1.,2.,-1.,2.,0.,-3.,0.,-3.,0.,-3./
  71 C
  72 C          Entry
  73 C
  74       ES=4*PI*ALFA
  75       WWSIG=0.
  76       IF(IDJETS(1).EQ.10.OR.IDJETS(2).EQ.10) GO TO 2
  77 C          Normal case
  78       IF((IDJETS(1).EQ.80.AND.IDJETS(2).EQ.-80).OR.
  79      $(IDJETS(1).EQ.90.AND.IDJETS(2).EQ.90).OR.
  80      $(IABS(IDJETS(1)).EQ.80.AND.IDJETS(2).EQ.90)) THEN
  81         DO 10 K=1,4
  82           P3(K)=P3WW(K)
  83           Q1(K)=PZERO(K,1)
  84           Q3(K)=PZERO(K,3)
  85 10      CONTINUE
  86         P3IS3=1.
  87         P3IS4=0.
  88         JQ1=1
  89         JQ3=3
  90         JW1=1
  91         JW2=2
  92         TX=T
  93         UX=U
  94 C          Crossed case
  95       ELSE
  96         DO 20 K=1,4
  97           P3(K)=P4WW(K)
  98           Q1(K)=PZERO(K,3)
  99           Q3(K)=PZERO(K,1)
 100 20      CONTINUE
 101         P3IS3=0.
 102         P3IS4=1.
 103         JQ1=3
 104         JQ3=1
 105         JW1=2
 106         JW2=1
 107         TX=U
 108         UX=T
 109       ENDIF
 110 C          Variables
 111       T1=-2.*(Q1(4)*P1(4)-Q1(1)*P1(1)-Q1(2)*P1(2)-Q1(3)*P1(3))
 112       U1=-2.*(Q1(4)*P2(4)-Q1(1)*P2(1)-Q1(2)*P2(2)-Q1(3)*P2(3))
 113       T3=-2.*(Q3(4)*P2(4)-Q3(1)*P2(1)-Q3(2)*P2(2)-Q3(3)*P2(3))
 114       U3=-2.*(Q3(4)*P1(4)-Q3(1)*P1(1)-Q3(2)*P1(2)-Q3(3)*P1(3))
 115       S13=2.*(Q1(4)*Q3(4)-Q1(1)*Q3(1)-Q1(2)*Q3(2)-Q1(3)*Q3(3))
 116 C          Jacobean for 4-body cross section in terms of squared
 117 C          matrix exement in narrow resonance approximation--
 118 C          1/((P**2-M**2)**2+M**2*GAM**2)=1/(2*M*GAM)*DELTA(P**2-M**2)
 119       FJAC=S/SCM*UNITS
 120       FJAC=FJAC*ALFA**4/(256.*PI*3.*S**2)
 121       AMW1=PJETS(5,1)
 122       AMW2=PJETS(5,2)
 123       GAM1=WGAM(JETTYP(1))
 124       GAM2=WGAM(JETTYP(2))
 125       FJAC=FJAC/(AMW1*GAM1*AMW2*GAM2)
 126       FJAC=FJAC*P(1)*P(2)/SQRT((P(1)**2+AMW1**2)*(P(2)**2+AMW2**2))
 127 C          Color factor
 128       IF(IABS(IDPAIR(1)).LT.10) FJAC=3.*FJAC
 129       IF(IABS(IDPAIR(3)).LT.10) FJAC=3.*FJAC
 130 C
 131 C          W+ W- pair decays
 132 C          Standard order is UP + UB --> W+ + W-
 133 C
 134       IF(.NOT.((JETTYP(1).EQ.2.AND.JETTYP(2).EQ.3).OR.(JETTYP(1).EQ.3
 135      1.AND.JETTYP(2).EQ.2))) GO TO 200
 136       FJAC=.5*FJAC*AQ(2,2)**4
 137 C
 138 C          Select W+ W- OR W- W+, swapping T and U for latter.
 139       IW1=JETTYP(1)
 140       IW2=JETTYP(2)
 141 C
 142 C          Select quarks
 143       IQ1=INITYP(1)
 144       IQ2=INITYP(2)
 145       IQ=IQ1/2
 146       CQ=AQDP(IQ,2)**2
 147       CV=AQDP(IQ,1)/S+EZDP*AQDP(IQ,4)/(S-ZM2)
 148       CA=EZDP*BQDP(IQ,4)/(S-ZM2)
 149       SGN=QSGN(IQ)
 150       ISWAPQ=1
 151       IF(SGN.LT.0.) ISWAPQ=-1
 152       IF(ISWAPQ.GT.0) THEN
 153         TT=TX
 154         UU=UX
 155         TT1=T1
 156         UU1=U1
 157         TT3=T3
 158         UU3=U3
 159         DO 122 K=1,4
 160           PP1(K)=P1(K)
 161           PP2(K)=P2(K)
 162           P3(K)=P3IS3*P3WW(K)+P3IS4*P4WW(K)
 163           Q1(K)=PZERO(K,JQ1)
 164           Q3(K)=PZERO(K,JQ3)
 165 122     CONTINUE
 166       ELSE
 167         TT=UX
 168         UU=TX
 169         TT1=U3
 170         UU1=T3
 171         TT3=U1
 172         UU3=T1
 173         DO 123 K=1,4
 174           PP1(K)=P1(K)
 175           PP2(K)=P2(K)
 176           P3(K)=P3IS4*P3WW(K)+P3IS3*P4WW(K)
 177           Q1(K)=PZERO(K,JQ3)
 178           Q3(K)=PZERO(K,JQ1)
 179 123     CONTINUE
 180       ENDIF
 181 C
 182       IF(IQ1.EQ.2*IQ) THEN
 183         TERM=WWTT(TT,UU,TT1,UU1,TT3,UU3)
 184         TERM=TERM-SGN*WWST(TT,UU,TT1,UU1,TT3,UU3,PP1,PP2)
 185         TERM=TERM+WWSS(TT,UU,TT1,UU1,TT3,UU3)
 186         WWSIG=TERM*QSAVE(2*IQ,1)*QSAVE(2*IQ+1,2)*FJAC
 187       ELSE
 188         TERM=WWTT(UU,TT,UU1,TT1,UU3,TT3)
 189         TERM=TERM-SGN*WWST(UU,TT,UU1,TT1,UU3,TT3,PP2,PP1)
 190         TERM=TERM+WWSS(UU,TT,UU1,TT1,UU3,TT3)
 191         WWSIG=TERM*QSAVE(2*IQ+1,1)*QSAVE(2*IQ,2)*FJAC
 192       ENDIF
 193 C
 194       RETURN
 195 C
 196 C          Z0 Z0 pair decays
 197 C          Standard order is UP + UB --> Z0 + Z0
 198 C
 199 200   IF(.NOT.(JETTYP(1).EQ.4.AND.JETTYP(2).EQ.4)) GO TO 300
 200       FJAC=.5*FJAC
 201 C
 202 C          Select quarks
 203       IQ1=INITYP(1)
 204       IQ2=INITYP(2)
 205       IQ=IQ1/2
 206       CV=AQDP(IQ,4)**2+BQDP(IQ,4)**2
 207       CA=2.*AQDP(IQ,4)*BQDP(IQ,4)
 208       CV1=AQDP(JQWW(1),4)**2+BQDP(JQWW(1),4)**2
 209       CA1=2.*AQDP(JQWW(1),4)*BQDP(JQWW(1),4)
 210       CV3=AQDP(JQWW(2),4)**2+BQDP(JQWW(2),4)**2
 211       CA3=2.*AQDP(JQWW(2),4)*BQDP(JQWW(2),4)
 212 C
 213       TERM=ZZALL(TX,UX,T1,U1,T3,U3,P1,P2)
 214       IF(INITYP(1).EQ.2*IQ) THEN
 215         WWSIG=TERM*QSAVE(2*IQ,1)*QSAVE(2*IQ+1,2)*FJAC
 216       ELSE
 217         WWSIG=TERM*QSAVE(2*IQ+1,1)*QSAVE(2*IQ,2)*FJAC
 218       ENDIF
 219 C
 220       RETURN
 221 C
 222 C          W+- Z0 pair decays
 223 C          Standard order is DN + UB --> W- + Z0
 224 C
 225 300   JW=JW1
 226       JZ=JW2
 227       ISGN=-ISIGN(1,IDJETS(JW))
 228       SGN=ISGN
 229       CV3=AQDP(JQWW(JZ),4)**2+BQDP(JQWW(JZ),4)**2
 230       CA3=2.*AQDP(JQWW(JZ),4)*BQDP(JQWW(JZ),4)
 231       FJAC=.5*FJAC*AQ(1,2)**2
 232 C
 233 C          Select quarks. Formulas are for DN UB --> W- Z0.
 234 C          Use symmetry for other cases.
 235       IQ1=INITYP(1)
 236       IQ2=INITYP(2)
 237       IQ=IQ1/2
 238 C          Find whether IQ1 should be fermion or antifermion.
 239       IF(IQ1.EQ.2*(IQ1/2)) THEN
 240         ISWAPQ=+1
 241         IFLI=IQ1/2
 242         IFLJ=IQ2/2
 243       ELSE
 244         ISWAPQ=-1
 245         IFLI=IQ2/2
 246         IFLJ=IQ1/2
 247       ENDIF
 248 C
 249       CS=AQDP(IQ,JETTYP(JW))*EZDP/(S-WM2)
 250       CT=AQDP(IQ,JETTYP(JW))*(AQDP(IFLJ,4)+BQDP(IFLJ,4))
 251       CU=AQDP(IQ,JETTYP(JW))*(AQDP(IFLI,4)+BQDP(IFLI,4))
 252 C
 253 C          SWAP T AND U AS NEEDED
 254       IF(ISWAPQ*ISGN.GT.0) THEN
 255         TT=TX
 256         UU=UX
 257         TT1=T1
 258         UU1=U1
 259         TT3=T3
 260         UU3=U3
 261         DO 321 K=1,4
 262           PP1(K)=P1(K)
 263           PP2(K)=P2(K)
 264 321     CONTINUE
 265       ELSE
 266         TT=UX
 267         UU=TX
 268         TT1=U1
 269         UU1=T1
 270         TT3=U3
 271         UU3=T3
 272         DO 323 K=1,4
 273           PP1(K)=P2(K)
 274           PP2(K)=P1(K)
 275 323     CONTINUE
 276       ENDIF
 277 C
 278       TERM=WZSS(TT,UU,TT1,UU1,TT3,UU3,PP1,PP2)
 279       TERM=TERM-SGN*WZST(TT,UU,TT1,UU1,TT3,UU3,PP1,PP2)
 280       TERM=TERM-SGN*WZSU(TT,UU,TT1,UU1,TT3,UU3,PP1,PP2)
 281       TERM=TERM+WZTU(TT,UU,TT1,UU1,TT3,UU3,PP1,PP2)
 282       WWSIG=TERM*QSAVE(IQ1,1)*QSAVE(IQ2,2)*FJAC
 283 C
 284       RETURN
 285 C
 286 C     Do Z+gamma or W+gamma 3-body subprocesses
 287 C
 288 2     CONTINUE
 289 C
 290 C          Z+gamma
 291 C          Standard order is UP + UB --> Z0 + gamma
 292 C
 293       IF(.NOT.(JETTYP(1).EQ.4.AND.JETTYP(2).EQ.1)) GO TO 505
 294       FJAC=S/SCM*P(1)/SQRT(P(1)**2+WMASS(4)**2)*UNITS
 295 C
 296 C          Select quarks
 297       IQ1=INITYP(1)
 298       IQ2=INITYP(2)
 299       IQ=IQ1/2
 300       A1=-AQ(IQ,4)
 301       B1=BQ(IQ,4)
 302       A2=-AQ(JQWW(1),4)
 303       B2=BQ(JQWW(1),4)
 304       DO K=1,5
 305         Q(K)=SNGL(P1WW(K))
 306         QB(K)=SNGL(P2WW(K))
 307         KK(K)=SNGL(P4WW(K))
 308         E(K)=SNGL(PZERO(K,1))
 309         EB(K)=SNGL(PZERO(K,2))
 310       END DO
 311 C
 312       IF(INITYP(1).EQ.2*IQ) THEN
 313         SMS=SMSZG(Q,QB,KK,E,EB,A1,B1,A2,B2)
 314         TERM=ES**3*(EQ3(IQ)/3.)**2*SMS/192./PI**4/WMASS(4)/WGAM(4)/S**2
 315         WWSIG=TERM*QSAVE(2*IQ,1)*QSAVE(2*IQ+1,2)*FJAC/2.
 316       ELSE
 317         SMS=SMSZG(QB,Q,KK,E,EB,A1,B1,A2,B2)
 318         TERM=ES**3*(EQ3(IQ)/3.)**2*SMS/192./PI**4/WMASS(4)/WGAM(4)/S**2
 319         WWSIG=TERM*QSAVE(2*IQ+1,1)*QSAVE(2*IQ,2)*FJAC/2.
 320       ENDIF
 321 505   IF(.NOT.(JETTYP(1).EQ.1.AND.JETTYP(2).EQ.4)) GO TO 509
 322       FJAC=S/SCM*P(2)/SQRT(P(2)**2+WMASS(4)**2)*UNITS
 323 C
 324 C          Select quarks
 325       IQ1=INITYP(1)
 326       IQ2=INITYP(2)
 327       IQ=IQ1/2
 328       A1=-AQ(IQ,4)
 329       B1=BQ(IQ,4)
 330       A2=-AQ(JQWW(2),4)
 331       B2=BQ(JQWW(2),4)
 332       DO K=1,5
 333         Q(K)=SNGL(P1WW(K))
 334         QB(K)=SNGL(P2WW(K))
 335         KK(K)=SNGL(P3WW(K))
 336         E(K)=SNGL(PZERO(K,1))
 337         EB(K)=SNGL(PZERO(K,2))
 338       END DO
 339 C
 340       IF(INITYP(1).EQ.2*IQ) THEN
 341         SMS=SMSZG(Q,QB,KK,E,EB,A1,B1,A2,B2)
 342         TERM=ES**3*(EQ3(IQ)/3.)**2*SMS/192./PI**4/WMASS(4)/WGAM(4)/S**2
 343         WWSIG=TERM*QSAVE(2*IQ,1)*QSAVE(2*IQ+1,2)*FJAC/2.
 344       ELSE
 345         SMS=SMSZG(QB,Q,KK,E,EB,A1,B1,A2,B2)
 346         TERM=ES**3*(EQ3(IQ)/3.)**2*SMS/192./PI**4/WMASS(4)/WGAM(4)/S**2
 347         WWSIG=TERM*QSAVE(2*IQ+1,1)*QSAVE(2*IQ,2)*FJAC/2.
 348       ENDIF
 349
 350 C          W+- GM pair decays
 351 C          Standard order is DN + UB --> W- + GM
 352 C
 353 C          Swap if W is jet 2
 354 509   IF (ABS(IDJETS(1)).EQ.80.OR.ABS(IDJETS(2)).EQ.80) THEN
 355       IF(IDJETS(2).EQ.10) THEN
 356         DO 510 K=1,4
 357           P3(K)=P3WW(K)
 358           Q1(K)=PZERO(K,1)
 359 510     CONTINUE
 360         AMASS3=PJETS(5,1)
 361         JW=1
 362         JG=2
 363         TX=T
 364         UX=U
 365       ELSE
 366         DO 520 K=1,4
 367           P3(K)=P4WW(K)
 368           Q1(K)=PZERO(K,1)
 369 520     CONTINUE
 370         AMASS3=PJETS(5,2)
 371         JW=2
 372         JG=1
 373         TX=U
 374         UX=T
 375       ENDIF
 376       IF(IDJETS(JW).EQ.80) THEN
 377         IW=2
 378       ELSE
 379         IW=3
 380       ENDIF
 381 C
 382       T1=-2.*(Q1(4)*P1(4)-Q1(1)*P1(1)-Q1(2)*P1(2)-Q1(3)*P1(3))
 383       U1=-2.*(Q1(4)*P2(4)-Q1(1)*P2(1)-Q1(2)*P2(2)-Q1(3)*P2(3))
 384 C          Jacobean
 385       FJAC=S/SCM*UNITS
 386       FJAC=FJAC*P(JW)/SQRT(P(JW)**2+WM2)
 387 C          Sum over quarks. Formulas are for DN UB --> W- GM.
 388 C          Use symmetry for other cases.
 389       IQ1=INITYP(1)
 390       IQ2=INITYP(2)
 391       IQ=IQ1/2
 392       IF(2*IQ.EQ.IQ1) THEN
 393         LQK1=.TRUE.
 394       ELSE
 395         LQK1=.FALSE.
 396       ENDIF
 397 C          Swap t and u as necessary
 398       IF((LQK1.AND.IW.EQ.3).OR.(.NOT.LQK1.AND.IW.EQ.2)) THEN
 399         TT=TX
 400         UU=UX
 401         TT1=T1
 402         UU1=U1
 403       ELSE
 404         TT=UX
 405         UU=TX
 406         TT1=U1
 407         UU1=T1
 408       ENDIF
 409 C          Lepton or quark pointer
 410       IL=IABS(IDPAIR(1))
 411       IF(IL.GT.6) IL=IL-4
 412 C
 413 C          Matrix element - properly crossed variables.
 414 C          Remember PZERO(K,1) is always the fermion.
 415       IF(LQK1) THEN
 416         P1DQ2=P1(4)*PZERO(4,2)-P1(1)*PZERO(1,2)-P1(2)*PZERO(2,2)
 417      $  -P1(3)*PZERO(3,2)
 418         P2DQ1=P2(4)*PZERO(4,1)-P2(1)*PZERO(1,1)-P2(2)*PZERO(2,1)
 419      $  -P2(3)*PZERO(3,1)
 420       ELSE
 421         P1DQ2=P2(4)*PZERO(4,2)-P2(1)*PZERO(1,2)-P2(2)*PZERO(2,2)
 422      $  -P2(3)*PZERO(3,2)
 423         P2DQ1=P1(4)*PZERO(4,1)-P1(1)*PZERO(1,1)-P1(2)*PZERO(2,1)
 424      $  -P1(3)*PZERO(3,1)
 425       ENDIF
 426       TERM=ALFA**2/(8.*SIN2W*S**2)*TBRWW(IW,JW)*RBRWW(IL,IW,JW)
 427      $*(-1./3.+UU/(TT+UU))**2/(TT*UU)*(4.*P2DQ1**2+4.*P1DQ2**2)
 428       WWSIG=TERM*QSAVE(IQ1,1)*QSAVE(IQ2,2)*FJAC
 429       END IF
 430 C
 431       RETURN
 432       END