ISAJET/code/drllyn.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE DRLLYN
   3 C
   4 C          Generate QMW (and QTW) for DRELLYAN or HIGGS event using
   5 C          integrated cross section. Then generate decay -- for HIGGS,
   6 C          the mode must be chosen using the integrated cross sections
   7 C          because of interference with W+W->W+W scattering.
   8 C
   9 C          Note that NOGOOD calls the cross section.
  10 C
  11 C          Ver. 6.40: Add technicolor resonances. Use logs for QDEN,
  12 C          PTDEN, WTFAC, etc. Also scale QMW generation by QMAX.
  13 C
  14 C          Ver. 7.01: Correct QDEN to correspond to correct fit form:
  15 C          SIGMA = ANOMR(K)*(QMAX**2/Q**2)**QPOW
  16 C          See QFUNC.
  17 C
  18 C          Ver. 7.14: Add SUSY Higgs
  19 C          Ver. 7.15: Fix bug with THETAW limits by adding epsilon to
  20 C          allowed range. Check for possible invalid Higgs decays.
  21 C
  22 #if defined(CERNLIB_IMPNONE)
  23       IMPLICIT NONE
  24 #endif
  25 #include "isajet/itapes.inc"
  26 #include "isajet/jetsig.inc"
  27 #include "isajet/totals.inc"
  28 #include "isajet/q1q2.inc"
  29 #include "isajet/partcl.inc"
  30 #include "isajet/pjets.inc"
  31 #include "isajet/pinits.inc"
  32 #include "isajet/wcon.inc"
  33 #include "isajet/primar.inc"
  34 #include "isajet/dylim.inc"
  35 #include "isajet/const.inc"
  36 #include "isajet/jetpar.inc"
  37 #include "isajet/jetlim.inc"
  38 #include "isajet/wgen.inc"
  39 #include "isajet/dypar.inc"
  40 #include "isajet/keys.inc"
  41 #include "isajet/hcon.inc"
  42 #include "isajet/isloop.inc"
  43 #include "isajet/idrun.inc"
  44 #include "isajet/xmssm.inc"
  45 #include "isajet/listss.inc"
  46 C
  47       DIMENSION X(2)
  48       EQUIVALENCE (X(1),X1)
  49       DIMENSION PREST(5),PL(5),EL(3),EML(3),EMSQL(3)
  50       DIMENSION WTFAC(3)
  51       LOGICAL NOGOOD
  52       LOGICAL YGENJ
  53       DIMENSION BRANCH(29),LISTJ(29),LISTW(5)
  54       REAL ACOSH,XXX,ASINH,CHOOSE,RANF,SUM,WTFAC,PTDEN,QDEN,ETA,QPW,
  55      $S12,BRANCH,SUMBR,BRMODE,AMASS,BRINV,TRY,EMSQL,EL,PL12,PREST,
  56      $COSTHL,THL,PHL,PTL,SGN,PL,BP,PLPL,PLMN,AMINI,AMFIN,PINI,PFIN,
  57      $ QPL,QMN,AM1SQ,AM2SQ,ROOT,P1PL,P1MN,P2PL,P2MN,X,EML
  58       INTEGER NTRY,K,IQ1,IQ2,IFL1,IFL2,LISTJ,IQ,NTRY2,IFL,LISTW,I
  59       REAL ZZSTAR
  60       INTEGER IZSTAR,JVIR,N0J
  61 C
  62       DATA LISTJ/
  63      $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
  64      $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,
  65      $10,80,-80,90/
  66       DATA LISTW/10,80,-80,90,92/
  67       ACOSH(XXX)=ALOG(XXX+SQRT(XXX**2-1.))
  68       ASINH(XXX)=ALOG(XXX+SQRT(XXX**2+1.))
  69 C
  70 C          Entry
  71 C
  72       NPTCL=0
  73       NTRY=0
  74 200   CONTINUE
  75       SIGMA=0.
  76       WT=1.
  77     1 CONTINUE
  78       NTRY=NTRY+1
  79       IF(NTRY.GT.NTRIES) GO TO 999
  80       SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
  81       NKINPT=NKINPT+1
  82       SIGMA=0.
  83       WT=1.
  84       DO 2 K=1,3
  85     2 SIGSL(K)=0
  86 C            Choose interval for cross section calculation
  87       CHOOSE=RANF()
  88       SUM=0.
  89       DO 3 K=NKL,NKH
  90         SUM=SUM+QSELWT(K)
  91         IF(CHOOSE.LE.SUM) GO TO 30
  92 3     CONTINUE
  93 30    KSEL=K
  94 C          Generate QTW in selected region
  95       IF(.NOT.FIXQT) THEN
  96         ETA=(PTGN(1,K)+PTGN(2,K)*RANF())**PTGN(3,K)
  97         PTSEL(K)=SQRT(ETA-RNU2(K))
  98         PTDEN=ALOG(ETA)*PTPOW(K)
  99         WTFAC(1)=ALOG(ABS(PTGN(2,K)))+ALOG(ABS(PTGN(3,K)))
 100      1  +ALOG(ABS(PTSEL(K)**2+RNU2(K)))*((PTGN(3,K)-1.)/PTGN(3,K))
 101         PT(3)=PTSEL(K)
 102       ELSE
 103         PTDEN=0.
 104         WTFAC(1)=-1000.
 105       ENDIF
 106 C          Generate QMW
 107       IF(.NOT.FIXQM) THEN
 108         IF(.NOT.K.EQ.2) THEN
 109           QSEL(K)=QMAX**2*(QGEN(1,K)+QGEN(2,K)*RANF())**QGEN(3,K)
 110           QDEN=ALOG(QSEL(K)/QMAX**2)*QPOW(K)
 111           WTFAC(2)=ALOG(ABS(QGEN(2,K)))+ALOG(ABS(QGEN(3,K)))
 112      1    +ALOG(QSEL(K)/QMAX**2)*((QGEN(3,K)-1.)/QGEN(3,K))
 113      2    +ALOG(QMAX**2)
 114           QSEL(K)=SQRT(QSEL(K))
 115           QMW=QSEL(K)
 116         ELSE
 117           ETA=QGEN(3,K)*TAN(QGEN(1,K)+QGEN(2,K)*RANF())
 118           QSEL(K)=SQRT(ETA+EMSQ)
 119           WTFAC(2)=ALOG(QGEN(2,K))+ALOG(QGEN(3,K))
 120      1    +ALOG((ETA/QGEN(3,K))**2+1.)
 121           QMW=QSEL(K)
 122           QDEN=ALOG((QMW**2-EMSQ)**2+EMGAM**2)
 123         ENDIF
 124       ELSE
 125         QDEN=0.
 126         WTFAC(2)=-1000.
 127       ENDIF
 128       SIGSL(K)=EXP(ANORM(K)-PTDEN-QDEN)
 129 C
 130       IF(STDDY) THEN
 131         WT=EXP(WTFAC(2)-ALOG(QSELWT(K)))
 132       ELSE
 133         WT=EXP(WTFAC(1)+WTFAC(2)-ALOG(QSELWT(K)))
 134       ENDIF
 135       QTW=PT(3)
 136       YW=YWMIN+(YWMAX-YWMIN)*RANF()
 137       WT=WT*(YWMAX-YWMIN)
 138       PHIW=PHWMIN+(PHWMAX-PHWMIN)*RANF()
 139       PHI(3)=AMOD(PHIW+PI,2.*PI)
 140       QPW=SQRT(QTW**2+QMW**2)*SINH(YW)
 141       QW=SQRT(QTW**2+QPW**2)
 142       THW=QPW/QW
 143       IF(ABS(THW).GT.1.) THW=SIGN(1.,THW)
 144       THW=ACOS(THW)
 145       IF(THW.LT.THWMIN-1.E-6.OR.THW.GT.THWMAX+1.E-6) GOTO 1
 146       XW=QPW/HALFE
 147       IF(XW.LT.XWMIN.OR.XW.GT.XWMAX) GOTO 1
 148       IF(.NOT.STDDY) THEN
 149         IF(.NOT.YGENJ(3)) GOTO 1
 150         P(3)=PT(3)/STH(3)
 151         XJ(3)=P(3)*CTH(3)/HALFE
 152         IF(XJ(3).LT.XJMIN(3).OR.XJ(3).GT.XJMAX(3)) GOTO 1
 153       ENDIF
 154 C
 155 C          Check integrated cross section
 156 C
 157       IF(NOGOOD(2)) GO TO 1
 158       SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
 159       NWGEN=NWGEN+1
 160       S12=QMW**2
 161 C
 162 C          No decay for KKG:
 163 C          For compatibility reasons, the jet is still the 3rd one.
 164 C          Jets 1 and 2 (W decay products) are voided; no decay step.
 165 C
 166       IF(KEYS(11)) THEN
 167         DO 50 I=1,2
 168           P(I)=0.
 169           PT(I)=0.
 170           CTH(I)=0.
 171           PHI(I)=0.
 172           EMSQL(I)=0.
 173           IDJETS(I)=0
 174 50      CONTINUE
 175         GOTO 350
 176       ENDIF
 177 C
 178 C          Select W decay mode
 179 C          QMW dependence neglected in branching ratios
 180 C          BRANCH is cum. br. with heavy modes subtracted.
 181 C
 182       IF(KEYS(3)) THEN
 183         BRANCH(1)=0.
 184         SUMBR=0.
 185         DO 105 IQ1=2,25
 186           IQ2=MATCH(IQ1,JWTYP)
 187           IF(IQ2.EQ.0) THEN
 188             BRMODE=0.
 189           ELSE
 190             BRMODE=WCBR(IQ1,JWTYP)-WCBR(IQ1-1,JWTYP)
 191             IFL1=LISTJ(IQ1)
 192             IFL2=LISTJ(IQ2)
 193             IF(S12.LE.(AMASS(IFL1)+AMASS(IFL2))**2) BRMODE=0.
 194           ENDIF
 195           BRANCH(IQ1)=BRANCH(IQ1-1)+BRMODE
 196           SUMBR=SUMBR+BRMODE
 197 105     CONTINUE
 198         BRINV=1./SUMBR
 199 C
 200         TRY=RANF()
 201         DO 110 IQ=1,25
 202           IF(TRY.LT.BRANCH(IQ)*BRINV.AND.MATCH(IQ,JWTYP).NE.0) THEN
 203             JETTYP(1)=IQ
 204             JETTYP(2)=MATCH(IQ,JWTYP)
 205             GO TO 120
 206           ENDIF
 207 110     CONTINUE
 208       ENDIF
 209 C
 210 120   IF(GOMSSM) THEN
 211         IFL1=LISTSS(JETTYP(1))
 212         IFL2=LISTSS(JETTYP(2))
 213       ELSE
 214         IFL1=LISTJ(JETTYP(1))
 215         IFL2=LISTJ(JETTYP(2))
 216       ENDIF
 217 C
 218 C          Select masses of decay products. These are just normal masses
 219 C          except for Z Z* decay of Higgs, where one is virtual.
 220 C
 221       EML(1)=AMASS(IFL1)
 222       EML(2)=AMASS(IFL2)
 223       IF(KEYS(7).AND.EML(1)+EML(2).GT.QMW) THEN
 224 C          WW* or ZZ* decay - generate/check W* or Z* mass
 225         IF((IABS(IFL1).EQ.80.AND.IABS(IFL2).EQ.80)
 226      $  .OR.(IFL1.EQ.90.AND.IFL2.EQ.90)) THEN
 227           IZSTAR=3-2*RANF()
 228           IF(GOMSSM) THEN
 229             JVIR=JETTYP(IZSTAR)-76
 230           ELSE
 231             JVIR=JETTYP(IZSTAR)-25
 232           ENDIF
 233           EML(IZSTAR)=ZZSTAR(QMW,JVIR)
 234           IF(EML(IZSTAR).LT.ZSTARS(JVIR,IZSTAR)) GO TO 200
 235 C          Other decay - invalid for this QMW
 236         ELSE
 237           GO TO 200
 238         ENDIF
 239       ENDIF
 240 C
 241 C          Generate W decay in its rest frame and compare with SIGDY2.
 242 C          First set up momenta of decay products:
 243 C
 244       EMSQL(1)=EML(1)**2
 245       EMSQL(2)=EML(2)**2
 246       EL(1)=(S12+EMSQL(1)-EMSQL(2))/(2.*QMW)
 247       EL(2)=(S12+EMSQL(2)-EMSQL(1))/(2.*QMW)
 248       PL12=SQRT((S12-(EML(1)+EML(2))**2)*(S12-(EML(1)-EML(2))**2))
 249      $/(2.*QMW)
 250 C          W momentum
 251       PREST(1)=QTW*COS(PHIW)
 252       PREST(2)=QTW*SIN(PHIW)
 253       PREST(3)=QPW
 254       PREST(4)=SQRT(QW**2+QMW**2)
 255       PREST(5)=QMW
 256       NTRY2=0
 257 C          Generate next W decay
 258 20    CONTINUE
 259       NTRY2=NTRY2+1
 260       IF(NTRY2.GT.NTRIES) GO TO 999
 261       COSTHL=2.*RANF()-1.
 262       THL=ACOS(COSTHL)
 263       PHL=2.*PI*RANF()
 264       PTL=PL12*SIN(THL)
 265 C
 266       DO 300 I=1,2
 267         SGN=3-2*I
 268         PL(1)=SGN*PTL*COS(PHL)
 269         PL(2)=SGN*PTL*SIN(PHL)
 270         PL(3)=SGN*PL12*COSTHL
 271         PL(4)=EL(I)
 272         PL(5)=EML(I)
 273 C          Boost with W momentum
 274         BP=0.
 275         DO 310 K=1,3
 276 310     BP=BP+PL(K)*PREST(K)
 277         BP=BP/PREST(5)
 278         DO 320 K=1,3
 279 320     PL(K)=PL(K)+PREST(K)*PL(4)/PREST(5)
 280      $  +PREST(K)*BP/(PREST(4)+PREST(5))
 281         PL(4)=PL(4)*PREST(4)/PREST(5)+BP
 282 C          Fill common blocks
 283         PT(I)=SQRT(PL(1)**2+PL(2)**2)
 284         P(I)=SQRT(PT(I)**2+PL(3)**2)
 285         IF(PT(I).GT.0.) THEN
 286           PHI(I)=ATAN2(PL(2),PL(1))
 287         ELSE
 288           PHI(I)=(I-1)*PI
 289         ENDIF
 290         IF(PHI(I).LT.0.) PHI(I)=PHI(I)+2.*PI
 291         CTH(I)=PL(3)/P(I)
 292         STH(I)=PT(I)/P(I)
 293         TH(I)=ACOS(CTH(I))
 294         XJ(I)=PL(3)/HALFE
 295         IF(CTH(I).GT.0.) THEN
 296           PLPL=PL(4)+PL(3)
 297           PLMN=(PT(I)**2+EMSQL(I))/PLPL
 298         ELSE
 299           PLMN=PL(4)-PL(3)
 300           PLPL=(PT(I)**2+EMSQL(I))/PLMN
 301         ENDIF
 302         YJ(I)=.5*ALOG(PLPL/PLMN)
 303 300   CONTINUE
 304 C
 305 C          Test cross section
 306 C          Extra kinematics for W+W->W+W
 307 C
 308       IF(KEYS(7).OR.KEYS(9)) THEN
 309         SHAT=S12
 310         IF(GOMSSM) THEN
 311           AMINI=AMASS(LISTSS(INITYP(1)))
 312         ELSE
 313           AMINI=AMASS(LISTJ(INITYP(1)))
 314         ENDIF
 315         AMFIN=EML(1)
 316         PINI=.5*SQRT(S12-4.*AMINI**2)
 317         PFIN=PL12
 318         THAT=AMINI**2+AMFIN**2-.5*S12+2.*PINI*PFIN*COSTHL
 319         UHAT=AMINI**2+AMFIN**2-.5*S12-2.*PINI*PFIN*COSTHL
 320       ENDIF
 321 C
 322 C          Check W decay
 323 C
 324       IF(NOGOOD(3)) GO TO 20
 325 C
 326 C          Check W decay with kinematic limits
 327 C
 328       IF(NOGOOD(4)) GO TO 200
 329 350   NKEEP=NKEEP+1
 330 C
 331 C            Set PBEAM
 332 C
 333       PBEAM(1)=(1.-X1)*HALFE
 334       PBEAM(2)=(1.-X2)*HALFE
 335       IF(NJET.LT.3) GO TO 502
 336       IFL=LISTJ(JETTYP(3))
 337       EMSQL(3)=AMASS(IFL)**2
 338 502   CONTINUE
 339 C
 340 C          Set PJETS
 341 C
 342       IF(KEYS(11)) THEN
 343         N0J=3
 344       ELSE
 345         N0J=1
 346       ENDIF
 347       DO 501 I=N0J,NJET
 348         PJETS(3,I)=P(I)*CTH(I)
 349         PJETS(1,I)=PT(I)*COS(PHI(I))
 350         PJETS(2,I)=PT(I)*SIN(PHI(I))
 351         PJETS(4,I)=SQRT(P(I)**2+EMSQL(I))
 352         PJETS(5,I)=SQRT(EMSQL(I))
 353         IF(KEYS(7).AND.GOMSSM) THEN
 354           IDJETS(I)=LISTSS(JETTYP(I))
 355         ELSE
 356           IDJETS(I)=LISTJ(JETTYP(I))
 357         ENDIF
 358 501   CONTINUE
 359 C          No technicolor IDENT's defined, so...
 360       IF(KEYS(3)) THEN
 361         IDENTW=LISTW(JWTYP)
 362       ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
 363         IDENTW=81
 364       ELSEIF(KEYS(7).AND.GOMSSM) THEN
 365         IDENTW=IHTYPE
 366       ELSEIF(KEYS(11)) THEN
 367         IDENTW=92
 368       ELSE
 369         IDENTW=0
 370       ENDIF
 371 C          W momentum in /PJETS/
 372       IF(KEYS(11)) THEN
 373         QWJET(1)=QTW*COS(PHIW)
 374         QWJET(2)=QTW*SIN(PHIW)
 375         QWJET(3)=QPW
 376         QWJET(4)=SQRT(QW**2+QMW**2)
 377         QWJET(5)=QMW
 378       ELSE
 379         DO 503 K=1,4
 380 503     QWJET(K)=PJETS(K,1)+PJETS(K,2)
 381         QWJET(5)=QMW
 382       ENDIF
 383 C
 384 C          Set PINITS
 385       DO 504 I=1,2
 386         IF(KEYS(7).AND.GOMSSM) THEN
 387           IDINIT(I)=LISTSS(INITYP(I))
 388         ELSE
 389           IDINIT(I)=LISTJ(INITYP(I))
 390         ENDIF
 391         PINITS(5,I)=AMASS(IDINIT(I))
 392         PINITS(1,I)=0.
 393         PINITS(2,I)=0.
 394 504   CONTINUE
 395 C          Calculate total momentum
 396       QPL=QWJET(4)+QWJET(3)
 397       QMN=QWJET(4)-QWJET(3)
 398       IF(NJET.EQ.3) THEN
 399         QPL=QPL+PJETS(4,3)+PJETS(3,3)
 400         QMN=QMN+PJETS(4,3)-PJETS(3,3)
 401       ENDIF
 402 C          and solve initial kinematics
 403       AM1SQ=PINITS(5,1)**2
 404       AM2SQ=PINITS(5,2)**2
 405       ROOT=SQRT((QPL*QMN-AM1SQ-AM2SQ)**2-4.*AM1SQ*AM2SQ)
 406       P1PL=(QPL*QMN+AM1SQ-AM2SQ+ROOT)/(2.*QMN)
 407       P1MN=AM1SQ/P1PL
 408       P2MN=(QPL*QMN+AM2SQ-AM1SQ+ROOT)/(2.*QPL)
 409       P2PL=AM2SQ/P2MN
 410       PINITS(3,1)=.5*(P1PL-P1MN)
 411       PINITS(4,1)=.5*(P1PL+P1MN)
 412       PINITS(3,2)=.5*(P2PL-P2MN)
 413       PINITS(4,2)=.5*(P2PL+P2MN)
 414       RETURN
 415 C
 416 999   CALL PRTEVT(0)
 417       WRITE(ITLIS,9999) NTRIES
 418 9999  FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE AN',
 419      C' EVENT. CHECK LIMITS OR INCREASE NTRIES')
 420       STOP 99
 421       END