ISAJET/code/sigdy.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE SIGDY
   3 C
   4 C          Compute the Drell-Yan and Drell-Yan plus jet cross sections
   5 C          d(sigma)/d(qmw**2)d(qtw**2)d(yw)d(yj)
   6 C
   7 C          SIGMA    = cross section summed over quark types allowed by
   8 C                     JETTYPE3 and WTYPE cards.
   9 C          SIGS(I)  = partial cross section for I1 + I2 --> I3 + I4.
  10 C          INOUT(I) = IOPAK**3*I4 + IOPAK**2*I3 + IOPAK*I2 + I1
  11 C                     using JETTYPE code.
  12 C
  13 C          QT cutoff for W+JET taken from Parisi and Petronzio,
  14 C          Nucl Phys B154, 427
  15 C          qk + gl --> qk + w suppressed at low QTW by extra factor
  16 C          of qtw**2/(qtw**2+qt2cut(qmw))
  17 C
  18 C          Ver 7.17: include top mass for gb --> Wt and gt --> Zt
  19 C          with no extra qt suppression factor. Note we do NOT include
  20 C          gt --> Wb; while this process makes sense for qt >> m_t,
  21 C          it has a pole in the physical region at low qt from the
  22 C          on-shell decay t --> Wb. We let Q**2 --> Q**2 + m_t**2
  23 C          in the scale for the parton distributions.
  24 C
  25 C          Ver 7.32: Rewrite AJLWT for gb --> Wt, etc., in terms of
  26 C          scaled variables, and restore SWT**5 later to avoid
  27 C          floating errors on VMS.
  28 C
  29 C          Ver 7.41: Recalculate COUT for each mass(!).
  30 C
  31 #if defined(CERNLIB_IMPNONE)
  32       IMPLICIT NONE
  33 #endif
  34 #include "isajet/itapes.inc"
  35 #include "isajet/qcdpar.inc"
  36 #include "isajet/jetpar.inc"
  37 #include "isajet/primar.inc"
  38 #include "isajet/q1q2.inc"
  39 #include "isajet/jetsig.inc"
  40 #include "isajet/qsave.inc"
  41 #include "isajet/wcon.inc"
  42 #include "isajet/const.inc"
  43 #include "isajet/nodcay.inc"
  44 C
  45       REAL X(2)
  46       REAL Z,S,T,U,QMW2,QZW,EHAT,Q2SAVE,YHAT,EY,P3Z,P1,P2,AMASS,ANEFF,
  47      $SIG0,DENOM,QT2CUT,SIGT,SIGU,FAC,PROP,FACTOR,SIG,AMT,AMT2,SWT,
  48      $P1WT,P2WT,X1WT,X2WT,TWT,UWT,Q2,QFCN,STRUC,XX,ACOSH,ATANH,P2M,P1M
  49       REAL AMI2,AMF2,EFWT
  50       REAL AJLWT,AJLZT1,AJLZT2,A2,A2B2,QQ,TM2
  51       INTEGER I,IQ,IH,IQ1,IFL,IQ2,IW
  52       INTEGER NZERO(4)
  53       REAL AMFAC(13)
  54       INTEGER NUTYP(25)
  55       INTEGER IFL1,IFL2
  56       REAL TERM
  57       EQUIVALENCE (S,SHAT),(T,THAT),(U,UHAT),(X(1),X1)
  58 C
  59       DATA NZERO/11,9,9,11/
  60       DATA AMFAC/11*0.,2*1./
  61       DATA NUTYP/13*0,1,1,0,0,1,1,0,0,1,1,0,0/
  62 C
  63 C          Functions
  64       ACOSH(Z)=ALOG(Z+SQRT(Z**2-1.))
  65       ATANH(Z)=.5*ALOG((1.+Z)/(1.-Z))
  66       PROP(I)=1./((QMW2-WMASS(I)**2)**2+(WMASS(I)*WGAM(I))**2)
  67 C          Qt cutoff function
  68       QT2CUT(QMW)=CUTOFF*QMW**CUTPOW
  69 C          Parton distributions
  70       QFCN(XX,IQ,IH)=STRUC(XX,QSQ+AMT2,IQ,IDIN(IH))/XX
  71 C          Integrated matrix elements JLint from FORM
  72       AJLWT(S,T,QQ,TM2)=
  73      $ - 32*QQ**3*S*T + 32*QQ**3*S*TM2 + 32*QQ**2*S**2*T
  74      $ + 32*QQ**2*S*T**2 - 16*QQ**2*S*T*TM2 - 16*QQ**2*S*TM2**2
  75      $ - 16*QQ*S**3*T + 16*QQ*S**3*TM2 - 16*QQ*S**2*T*TM2
  76      $ - 16*QQ*S*T**3 + 32*QQ*S*T**2*TM2 - 16*QQ*S*T*TM2**2
  77      $ - 8*S**3*T*TM2 + 8*S**3*TM2**2 - 16*S**2*T**2*TM2
  78      $ + 16*S**2*T*TM2**2 - 16*S**2*TM2**3 - 8*S*T**3*TM2
  79      $ + 8*S*T**2*TM2**2 - 8*S*T*TM2**3 + 8*S*TM2**4
  80 C
  81       AJLZT1(S,T,QQ,TM2)=
  82      $ + A2 * ( - 96*QQ**2*S*T*TM2 + 96*QQ**2*S*TM2**2
  83      $ + 96*QQ**2*T*TM2**2 - 96*QQ**2*TM2**3 + 96*QQ*S**2*T*TM2
  84      $ + 96*QQ*S*T**2*TM2 - 192*QQ*S*T*TM2**2 - 96*QQ*S*TM2**3
  85      $ - 96*QQ*T*TM2**3 + 192*QQ*TM2**4 + 16*S**3*T*TM2
  86      $ - 16*S**3*TM2**2 + 32*S**2*T**2*TM2 - 112*S**2*T*TM2**2
  87      $ + 80*S**2*TM2**3 + 16*S*T**3*TM2 - 112*S*T**2*TM2**2
  88      $ + 224*S*T*TM2**3 - 128*S*TM2**4 - 16*T**3*TM2**2
  89      $ + 80*T**2*TM2**3 - 128*T*TM2**4 + 64*TM2**5 )
  90       AJLZT2(S,T,QQ,TM2)=
  91      $ + A2B2 * ( - 16*QQ**3*S*T + 16*QQ**3*S*TM2 + 16*QQ**3*T*TM2
  92      $ - 16*QQ**3*TM2**2 + 16*QQ**2*S**2*T + 16*QQ**2*S*T**2
  93      $ + 32*QQ**2*S*T*TM2 - 80*QQ**2*S*TM2**2 - 80*QQ**2*T*TM2**2
  94      $ + 96*QQ**2*TM2**3 - 8*QQ*S**3*T + 8*QQ*S**3*TM2
  95      $ - 40*QQ*S**2*T*TM2 - 24*QQ*S**2*TM2**2 - 8*QQ*S*T**3
  96      $ - 40*QQ*S*T**2*TM2 + 80*QQ*S*T*TM2**2 + 96*QQ*S*TM2**3
  97      $ + 8*QQ*T**3*TM2 - 24*QQ*T**2*TM2**2 + 96*QQ*T*TM2**3
  98      $ - 144*QQ*TM2**4 - 16*S**3*T*TM2 + 16*S**3*TM2**2
  99      $ - 32*S**2*T**2*TM2 + 112*S**2*T*TM2**2 - 80*S**2*TM2**3
 100      $ - 16*S*T**3*TM2 + 112*S*T**2*TM2**2 - 224*S*T*TM2**3
 101      $ + 128*S*TM2**4 + 16*T**3*TM2**2 - 80*T**2*TM2**3
 102      $ + 128*T*TM2**4 - 64*TM2**5 )
 103 C
 104 C          Kinematics
 105 C
 106       QMW2=QMW**2
 107       QTMW=SQRT(QMW2+QTW**2)
 108       Q0W=QTMW*COSH(YW)
 109       QZW=QTMW*SINH(YW)
 110       QW=SQRT(QZW**2+QTW**2)
 111 C          Protect against errors
 112       IF(QW.NE.0.) THEN
 113         CTHW=QZW/QW
 114         STHW=QTW/QW
 115         IF(ABS(CTHW).LT.1.) THEN
 116           THW=ACOS(CTHW)
 117         ELSE
 118           CTHW=0.
 119           STHW=1.
 120           THW=.5*PI
 121         ENDIF
 122       ELSE
 123         CTHW=0.
 124         STHW=1.
 125         THW=.5*PI
 126       ENDIF
 127 C
 128       IF(STDDY) THEN
 129 C          Kinematics for standard Drell-Yan
 130         EHAT=QMW
 131         SHAT=QMW**2
 132         QSQ=SHAT
 133         Q2SAVE=QSQ
 134         YHAT=YW
 135         EY=EXP(YHAT)
 136         X1=EHAT/ECM*EY
 137         X2=EHAT/(ECM*EY)
 138       ELSE
 139 C          Kinematics for Drell-Yan plus jet
 140         P3Z=P(3)*CTH(3)
 141         SHAT=QMW2+2.*Q0W*P(3)-2.*QZW*P3Z+2.*PT(3)**2
 142         P1=.5*(P(3)+P3Z+Q0W+QZW)
 143         P2=.5*(P(3)-P3Z+Q0W-QZW)
 144         X1=P1/HALFE
 145         X2=P2/HALFE
 146         THAT=-2.*P1*(P(3)-P3Z)
 147         UHAT=-2.*P2*(P(3)+P3Z)
 148         QSQ=QTW**2
 149         QSQ=AMAX1(QSQ,4.)
 150         ANEFF=4.+QSQ/(QSQ+AMASS(5)**2)+QSQ/(QSQ+AMASS(6)**2)
 151         ALFQSQ=12.*PI/((33.-2.*ANEFF)*ALOG(QSQ/ALAM2))
 152         Q2SAVE=QSQ
 153         QSQ=SHAT
 154       ENDIF
 155 C
 156 C          Initialize
 157 C
 158       SIGMA=0.
 159       NSIGS=0
 160       DO 100 I=1,MXSIGS
 161         SIGS(I)=0.
 162 100   CONTINUE
 163       IF(X1.GE.1..OR.X2.GE.1.) RETURN
 164 C
 165 C          Compute structure functions
 166 C
 167       DO 110 IH=1,2
 168         DO 120 IQ=1,11
 169           QSAVE(IQ,IH)=STRUC(X(IH),QSQ,IQ,IDIN(IH))/X(IH)
 170 120     CONTINUE
 171         QSAVE(12,IH)=0
 172         QSAVE(13,IH)=0
 173 110   CONTINUE
 174       QSQ=Q2SAVE
 175 C
 176 C          Recompute COUT for this mass
 177 C
 178       DO 130 IW=1,4
 179         COUT(IW)=0.
 180         IF(.NOT.GODY(IW)) GO TO 130
 181         DO 140 IQ1=2,25
 182           IQ2=MATCH(IQ1,IW)
 183           IF(IQ2.EQ.0) GO TO 140
 184           IF(.NOT.(GOQ(IQ1,1).AND.GOQ(IQ2,2))) GO TO 140
 185           IF(NUTYP(IQ1)*NUTYP(IQ2).EQ.1.AND.NONUNU) GO TO 140
 186           IFL1=IQ1/2
 187           IFL2=IQ2/2
 188           IF(AMASS(IFL1)+AMASS(IFL2).GE.QMW) GO TO 140
 189           TERM=.5*(AQ(IFL1,IW)**2+BQ(IFL1,IW)**2)
 190           IF(IQ1.LE.13) TERM=3.*TERM
 191           COUT(IW)=COUT(IW)+TERM
 192 140     CONTINUE
 193 130   CONTINUE
 194 C
 195       IF(STDDY) GO TO 400
 196 C
 197 C          Compute cross section for types allowed by WTYPE and
 198 C          JETTYPE cards.
 199 C
 200 C          qk + gl --> qk + W
 201 C
 202       SIG0=ALFA**2*ALFQSQ*QMW2/(9.*SCM*S)*UNITS
 203       DENOM=S**2*EXP(.5*ALOG(QTW**4+QT2CUT(QMW)**2))
 204       SIGT=SIG0*(S**2+U**2+2.*T*QMW2)*(-T)/DENOM
 205       SIGU=SIG0*(S**2+T**2+2.*U*QMW2)*(-U)/DENOM
 206       DO 200 IW=1,4
 207         IF(.NOT.GODY(IW)) GO TO 200
 208         FAC=COUT(IW)*PROP(IW)
 209         DO 210 IQ=2,NZERO(IW)
 210           IF(.NOT.GOQ(IQ,3)) GO TO 210
 211           IQ1=MATCH(IQ,4)
 212           IQ1=MATCH(IQ1,IW)
 213           IF(IQ1.EQ.0) GO TO 210
 214           IFL=IQ/2
 215           FACTOR=FAC*(AQ(IFL,IW)**2+BQ(IFL,IW)**2)
 216      $    *QTW**2/(QTW**2+QT2CUT(QMW))
 217           SIG=FACTOR*SIGT*QSAVE(IQ1,1)*QSAVE(1,2)
 218           CALL SIGFIL(SIG,IQ1,1,IW,IQ)
 219           SIG=FACTOR*SIGU*QSAVE(IQ1,2)*QSAVE(1,1)
 220           CALL SIGFIL(SIG,1,IQ1,IW,IQ)
 221 210     CONTINUE
 222 200   CONTINUE
 223 C
 224 C          bt,tp + gl -> bt,tp + W,Z
 225 C
 226       AMT=AMASS(6)
 227       AMT2=AMT**2
 228       Q2=QMW2
 229       DO 220 IW=2,4
 230         IF(.NOT.GODY(IW)) GO TO 220
 231         DO 230 IQ=NZERO(IW)+1,13
 232           IF(.NOT.GOQ(IQ,3)) GO TO 230
 233           IQ1=MATCH(IQ,4)
 234           IQ1=MATCH(IQ1,IW)
 235           IF(IQ1.EQ.0) GO TO 230
 236           IF(IQ1.GE.12.AND.IW.NE.4) GO TO 230
 237 C          Assign zero or top masses for initial/final quarks
 238           AMF2=AMT2*AMFAC(IQ)
 239           AMI2=AMT2*AMFAC(IQ1)
 240           EFWT=SQRT(P(3)**2+AMF2)
 241           SWT=QMW2+AMF2+2.*Q0W*EFWT-2.*QZW*P3Z+2.*PT(3)**2
 242 C
 243 C          qk + gl initial state
 244 C          Do kinematics using p(small) = 0 for gluon
 245 C
 246           P1WT=EFWT+P3Z+Q0W+QZW
 247           P1M=AMI2/P1WT
 248           P2WT=EFWT-P3Z+Q0W-QZW-P1M
 249           X1WT=.5*P1WT/HALFE
 250           X2WT=.5*P2WT/HALFE
 251           TWT=-P1WT*(EFWT-P3Z)-P1M*(EFWT+P3Z)+AMI2+AMF2
 252           UWT=-P2WT*(EFWT+P3Z)+AMF2
 253           IF(X1WT.LT.0.OR.X1WT.GT.1.OR.X2WT.LT.0.OR.X2WT.GT.1)
 254      $    GO TO 240
 255 C          Cross sections
 256           IF(IW.EQ.2.OR.IW.EQ.3) THEN
 257             SIG0=ALFA**2*ALFQSQ/(144*SCM*SWT)*UNITS
 258             SIG0=SIG0*(AQ(5,IW)**2+BQ(5,IW)**2)*COUT(IW)*PROP(IW)
 259             SIGU=SIG0*AJLWT(SWT/SWT,UWT/SWT,Q2/SWT,AMT2/SWT)*SWT*
 260      $      (SWT/(SWT-AMI2))**2*(SWT/(UWT-AMF2))**2
 261             SIG=SIGU*QFCN(X1WT,IQ1,1)*QFCN(X2WT,1,2)
 262             CALL SIGFIL(SIG,IQ1,1,IW,IQ)
 263           ELSEIF(IW.EQ.4) THEN
 264             SIG0=ALFA**2*ALFQSQ/(144*SCM*SWT)*UNITS
 265             SIG0=SIG0*COUT(IW)*PROP(IW)
 266             A2=AQ(6,IW)**2
 267             A2B2=AQ(6,IW)**2+BQ(6,IW)**2
 268             SIGU=SIG0*(AJLZT1(SWT/SWT,UWT/SWT,Q2/SWT,AMT2/SWT)+
 269      $      AJLZT2(SWT/SWT,UWT/SWT,Q2/SWT,AMT2/SWT))*SWT*
 270      $      (SWT/(SWT-AMI2))**2*(SWT/(UWT-AMF2))**2
 271             SIG=SIGU*QFCN(X1WT,IQ1,1)*QFCN(X2WT,1,2)
 272             CALL SIGFIL(SIG,IQ1,1,IW,IQ)
 273           ENDIF
 274 C
 275 C          gl + qk initial state
 276 C          Do kinematics  using p(small) = 0 for gluon
 277 C
 278 240       P2WT=EFWT-P3Z+Q0W-QZW
 279           P2M=AMI2/P2WT
 280           P1WT=EFWT+P3Z+Q0W+QZW-P2M
 281           X1WT=.5*P1WT/HALFE
 282           X2WT=.5*P2WT/HALFE
 283           TWT=-P1WT*(EFWT-P3Z)+AMF2
 284           UWT=-P2WT*(EFWT+P3Z)-P2M*(EFWT-P3Z)+AMI2+AMF2
 285           IF(X1WT.LT.0.OR.X1WT.GT.1.OR.X2WT.LT.0.OR.X2WT.GT.1)
 286      $    GO TO 230
 287 C          Cross sections
 288           IF(IW.EQ.2.OR.IW.EQ.3) THEN
 289             SIG0=ALFA**2*ALFQSQ/(144*SCM*SWT)*UNITS
 290             SIG0=SIG0*(AQ(5,IW)**2+BQ(5,IW)**2)*COUT(IW)*PROP(IW)
 291             SIGT=SIG0*AJLWT(SWT/SWT,TWT/SWT,Q2/SWT,AMT2/SWT)*SWT*
 292      $      (SWT/(SWT-AMI2))**2*(SWT/(TWT-AMF2)**2)
 293             SIG=SIGT*QFCN(X1WT,1,1)*QFCN(X2WT,IQ1,2)
 294             CALL SIGFIL(SIG,1,IQ1,IW,IQ)
 295           ELSEIF(IW.EQ.4) THEN
 296             SIG0=ALFA**2*ALFQSQ/(144*SCM*SWT)*UNITS
 297             SIG0=SIG0*COUT(IW)*PROP(IW)
 298             A2=AQ(6,IW)**2
 299             A2B2=AQ(6,IW)**2+BQ(6,IW)**2
 300             SIGU=SIG0*(AJLZT1(SWT/SWT,TWT/SWT,Q2/SWT,AMT2/SWT)+
 301      $      AJLZT2(SWT/SWT,TWT/SWT,Q2/SWT,AMT2/SWT))*SWT*
 302      $      (SWT/(SWT-AMI2))**2*(SWT/(TWT-AMF2))**2
 303             SIG=SIGU*QFCN(X1WT,1,1)*QFCN(X2WT,IQ1,2)
 304             CALL SIGFIL(SIG,1,IQ1,IW,IQ)
 305           ENDIF
 306 230     CONTINUE
 307 220   CONTINUE
 308 C
 309 C          qk + qb --> gl + W
 310 C
 311       IF(.NOT.GOQ(1,3)) RETURN
 312       SIG0=8.*ALFA**2*ALFQSQ*QMW2/(27.*SCM*S)*UNITS
 313       DENOM=S*EXP(.5*ALOG(QTW**4+QT2CUT(QMW)**2))
 314       SIG0=SIG0*(T**2+U**2+2.*S*QMW2)/DENOM
 315       DO 300 IW=1,4
 316         IF(.NOT.GODY(IW)) GO TO 300
 317         FAC=COUT(IW)*PROP(IW)
 318         DO 310 IQ1=2,11
 319           IQ2=MATCH(IQ1,IW)
 320           IF(IQ2.EQ.0) GO TO 310
 321           IFL=IQ1/2
 322           SIG=FAC*SIG0*(AQ(IFL,IW)**2+BQ(IFL,IW)**2)
 323      $    *QSAVE(IQ1,1)*QSAVE(IQ2,2)
 324           CALL SIGFIL(SIG,IQ1,IQ2,IW,1)
 325 310     CONTINUE
 326 300   CONTINUE
 327       RETURN
 328 C
 329 C          Standard Drell-Yan for QT=0.
 330 C
 331 400   CONTINUE
 332       SIG0=4.*PI*ALFA**2*QMW2/(9.*SCM)*UNITS
 333       DO 410 IW=1,4
 334         IF(.NOT.GODY(IW)) GO TO 410
 335         FAC=COUT(IW)*PROP(IW)
 336         DO 420 IQ1=2,13
 337           IQ2=MATCH(IQ1,IW)
 338           IF(IQ2.EQ.0) GO TO 420
 339           IFL=IQ1/2
 340           SIG=FAC*SIG0*(AQ(IFL,IW)**2+BQ(IFL,IW)**2)
 341      $    *QSAVE(IQ1,1)*QSAVE(IQ2,2)
 342           CALL SIGFIL(SIG,IQ1,IQ2,IW,0)
 343 420     CONTINUE
 344 410   CONTINUE
 345 C
 346       RETURN
 347       END