HERWIG/src/hwbtim.f

   1
   2 CDECK  ID>, HWBTIM.
   3
   4 *CMZ :-        -26/04/91  14.27.17  by  Federico Carminati
   5
   6 *-- Author :    Ian Knowles
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWBTIM(INITBR,INTERF)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     Constructs full 4-momentum & production vertices in time-like jet
  15
  16 C     initiated by INITBR, interference partner INTERF and spin density
  17
  18 C     RHOPAR(INITBR). DECPAR(INITBR) returns jet's spin density matrix.
  19
  20 C     Includes azimuthal angular correlations between branching planes
  21
  22 C     due to spin (if AZSPIN) using the algorithm of Knowles & Collins.
  23
  24 C     Ses Nucl. Phys. B304 (1988) 794 & Comp. Phys. Comm. 58 (1990) 271.
  25
  26 C-----------------------------------------------------------------------
  27
  28       INCLUDE 'HERWIG61.INC'
  29
  30       DOUBLE PRECISION HWR,DMIN,PT,EIKON,EINUM,EIDEN1,EIDEN2,EISCR,
  31
  32      & WT,SPIN,Z1,Z2,PRMAX,CAZ,CX,SX,ROHEP(3),RMAT(3,3),ZERO2(2)
  33
  34       INTEGER INITBR,INTERF,IPAR,JPAR,KPAR,LPAR,MPAR,NTRY,JOLD
  35
  36       LOGICAL EICOR,SWAP
  37
  38       EXTERNAL HWR
  39
  40       DATA ZERO2,DMIN/ZERO,ZERO,1.D-15/
  41
  42       IF (IERROR.NE.0) RETURN
  43
  44       JPAR=INITBR
  45
  46       KPAR=INTERF
  47
  48       IF ((JDAPAR(1,JPAR).NE.0).OR.(IDPAR(JPAR).EQ.13)) GOTO 30
  49
  50 C No branching, assign decay matrix
  51
  52       CALL HWVZRO(2,DECPAR(1,JPAR))
  53
  54       RETURN
  55
  56 C Advance up the leader
  57
  58 C     Find the parent and partner of J
  59
  60   10  IPAR=JMOPAR(1,JPAR)
  61
  62       KPAR=JPAR+1
  63
  64 C Generate new Rho
  65
  66       IF (JMOPAR(1,KPAR).EQ.IPAR) THEN
  67
  68 C        Generate Rho'
  69
  70          CALL HWBAZF(IPAR,JPAR,PHIPAR(1,IPAR),RHOPAR(1,IPAR),
  71
  72      &                                   ZERO2,RHOPAR(1,JPAR))
  73
  74       ELSE
  75
  76          KPAR=JPAR-1
  77
  78          IF (JMOPAR(1,KPAR).NE.IPAR)
  79
  80      &   CALL HWWARN('HWBTIM',100,*999)
  81
  82 C        Generate Rho''
  83
  84          CALL HWBAZF(IPAR,KPAR,RHOPAR(1,IPAR),PHIPAR(1,IPAR),
  85
  86      &                         DECPAR(1,KPAR),RHOPAR(1,JPAR))
  87
  88       ENDIF
  89
  90 C Generate azimuthal angle of J's branching
  91
  92   30  IF (JDAPAR(1,JPAR).EQ.0) THEN
  93
  94 C        Final state gluon
  95
  96          CALL HWVZRO(2,DECPAR(1,JPAR))
  97
  98          IF (JPAR.EQ.INITBR) RETURN
  99
 100          GOTO 70
 101
 102       ELSE
 103
 104 C Assign an angle to a branching using an M-function
 105
 106 C        Find the daughters of J
 107
 108          LPAR=JDAPAR(1,JPAR)
 109
 110          MPAR=JDAPAR(2,JPAR)
 111
 112 C Soft correlations
 113
 114          CALL HWUROT(PPAR(1,JPAR), ONE,ZERO,RMAT)
 115
 116          CALL HWUROF(RMAT,PPAR(1,KPAR),ROHEP)
 117
 118          PT=MAX(SQRT(ROHEP(1)*ROHEP(1)+ROHEP(2)*ROHEP(2)),DMIN)
 119
 120          EIKON=1.
 121
 122          SWAP=.FALSE.
 123
 124          EICOR=AZSOFT.AND.((IDPAR(LPAR).EQ.13).OR.(IDPAR(MPAR).EQ.13))
 125
 126          IF (EICOR) THEN
 127
 128 C           Rearrange s.t. LPAR is the (softest) gluon
 129
 130             IF (IDPAR(MPAR).EQ.13) THEN
 131
 132                IF (IDPAR(LPAR).NE.13.OR.
 133
 134      &             PPAR(4,MPAR).LT.PPAR(4,LPAR)) THEN
 135
 136                   SWAP=.TRUE.
 137
 138                   LPAR=MPAR
 139
 140                   MPAR=LPAR-1
 141
 142                ENDIF
 143
 144             ENDIF
 145
 146             EINUM=(PPAR(4,KPAR)*PPAR(4,LPAR))
 147
 148      &        *ABS(PPAR(2,LPAR)-PPAR(2,MPAR))
 149
 150             EIDEN1=(PPAR(4,KPAR)*PPAR(4,LPAR))-ROHEP(3)*PPAR(3,LPAR)
 151
 152             EIDEN2=PT*ABS(PPAR(1,LPAR))
 153
 154             EISCR=1.-(PPAR(5,MPAR)/PPAR(4,MPAR))**2
 155
 156      &           /MIN(PPAR(2,LPAR),PPAR(2,MPAR))
 157
 158             EIKON=EISCR+EINUM/MAX(EIDEN1-EIDEN2,DMIN)
 159
 160          ENDIF
 161
 162 C Spin correlations
 163
 164          WT=0.
 165
 166          SPIN=1.
 167
 168          IF (AZSPIN) THEN
 169
 170             Z1=PPAR(4,LPAR)/PPAR(4,JPAR)
 171
 172             Z2=1.-Z1
 173
 174             IF (IDPAR(JPAR).EQ.13.AND.IDPAR(LPAR).EQ.13) THEN
 175
 176                WT=Z1*Z2/(Z1/Z2+Z2/Z1+Z1*Z2)
 177
 178             ELSEIF (IDPAR(JPAR).EQ.13.AND.IDPAR(LPAR).LT.13) THEN
 179
 180                WT=-2.*Z1*Z2/(Z1*Z1+Z2*Z2)
 181
 182             ENDIF
 183
 184          ENDIF
 185
 186 C Assign the azimuthal angle
 187
 188          PRMAX=(1.+ABS(WT))*EIKON
 189
 190          NTRY=0
 191
 192    50    NTRY=NTRY+1
 193
 194          IF (NTRY.GT.NBTRY) CALL HWWARN('HWBTIM',101,*999)
 195
 196          CALL HWRAZM( ONE,CX,SX)
 197
 198          CALL HWUROT(PPAR(1,JPAR),CX,SX,RMAT)
 199
 200 C Determine the angle between the branching planes
 201
 202          CALL HWUROF(RMAT,PPAR(1,KPAR),ROHEP)
 203
 204          CAZ=ROHEP(1)/PT
 205
 206          PHIPAR(1,JPAR)=2.*CAZ*CAZ-1.
 207
 208          PHIPAR(2,JPAR)=2.*CAZ*ROHEP(2)/PT
 209
 210          IF (EICOR) EIKON=EISCR+EINUM/MAX(EIDEN1-EIDEN2*CAZ,DMIN)
 211
 212          IF (AZSPIN) SPIN=1.+WT*(RHOPAR(1,JPAR)*PHIPAR(1,JPAR)
 213
 214      &                          +RHOPAR(2,JPAR)*PHIPAR(2,JPAR))
 215
 216          IF (SPIN*EIKON.LT.HWR()*PRMAX) GOTO 50
 217
 218 C Construct full 4-momentum of L and M
 219
 220          JOLD=JPAR
 221
 222          IF (SWAP) THEN
 223
 224            PPAR(1,LPAR)=-PPAR(1,LPAR)
 225
 226            PPAR(1,MPAR)=-PPAR(1,MPAR)
 227
 228            JPAR=MPAR
 229
 230          ELSE
 231
 232            JPAR=LPAR
 233
 234          ENDIF
 235
 236          PPAR(2,LPAR)=0.
 237
 238          CALL HWUROB(RMAT,PPAR(1,LPAR),PPAR(1,LPAR))
 239
 240          PPAR(2,MPAR)=0.
 241
 242          CALL HWUROB(RMAT,PPAR(1,MPAR),PPAR(1,MPAR))
 243
 244 C Assign production vertex to L and M
 245
 246          CALL HWUDKL(IDPAR(JOLD),PPAR(1,JOLD),VPAR(1,LPAR))
 247
 248          CALL HWVSUM(4,VPAR(1,JOLD),VPAR(1,LPAR),VPAR(1,LPAR))
 249
 250          CALL HWVEQU(4,VPAR(1,LPAR),VPAR(1,MPAR))
 251
 252       ENDIF
 253
 254   60  IF (JDAPAR(1,JPAR).NE.0) GOTO 10
 255
 256 C Assign decay matrix
 257
 258       CALL HWVZRO(2,DECPAR(1,JPAR))
 259
 260 C Backtrack down the leader
 261
 262   70  IPAR=JMOPAR(1,JPAR)
 263
 264       KPAR=JDAPAR(1,IPAR)
 265
 266       IF (KPAR.EQ.JPAR) THEN
 267
 268 C        Develop the side branch
 269
 270          JPAR=JDAPAR(2,IPAR)
 271
 272          GOTO 60
 273
 274       ELSE
 275
 276 C        Construct decay matrix
 277
 278          CALL HWBAZF(IPAR,KPAR,DECPAR(1,JPAR),DECPAR(1,KPAR),
 279
 280      &                         PHIPAR(1,IPAR),DECPAR(1,IPAR))
 281
 282       ENDIF
 283
 284       IF (IPAR.EQ.INITBR) RETURN
 285
 286       JPAR=IPAR
 287
 288       GOTO 70
 289
 290   999 END
 291
 292 CDECK  ID>, HWBTOP.
 293
 294 *CMZ :-        -14/10/99  18.04.56  by  Mike Seymour
 295
 296 *-- Author :    Gennaro Corcella
 297
 298 C-----------------------------------------------------------------------
 299
 300       SUBROUTINE HWBTOP
 301
 302 C-----------------------------------------------------------------------
 303
 304       INCLUDE 'HERWIG61.INC'
 305
 306       DOUBLE PRECISION HWBVMC,HWR,HWUALF,HWUSQR,X(3),W,
 307
 308      & X3MIN,X3MAX,X1MIN,X1MAX,QSCALE,GLUFAC,R(3,3),M(3),
 309
 310      & E(3),AW,PTSQ,EM,EPS,MASDEP,A,B,C,GAMDEP,LAMBDA,
 311
 312      & PW(5),PT(5),PW1(5),CS,SN,EPG,QQ,RR,CC
 313
 314       INTEGER ID,ID3,IHEP,KHEP,WHEP,ICMF,K
 315
 316       EXTERNAL HWBVMC,HWUALF,HWUSQR,HWR
 317
 318       LAMBDA(A,B,C)=(A**2+B**2+C**2-2*A*B-2*B*C-2*C*A)/(4*A)
 319
 320 C---FIND AN UNTREATED CMF
 321
 322       ICMF=0
 323
 324       DO 10 IHEP=1,NHEP
 325
 326 C----FIND A DECAYING TOP QUARK
 327
 328  10     IF (ISTHEP(IHEP).EQ.155.AND.ISTHEP(JDAHEP(1,IHEP)).EQ.113
 329
 330      &       .AND.(IDHW(IHEP).EQ.6.OR.IDHW(IHEP).EQ.12))
 331
 332      &       ICMF=IHEP
 333
 334       IF (ICMF.EQ.0) RETURN
 335
 336       EM=PHEP(5,ICMF)
 337
 338       X3MIN=2*GCUTME/EM
 339
 340 C---GENERATE X(1),X(3) ACCORDING TO 1/((1-X(1))*X(3)**2)
 341
 342  100  CONTINUE
 343
 344 C-----AW=(MW/MT)**2
 345
 346       AW=(PHEP(5,JDAHEP(1,ICMF))/EM)**2
 347
 348 C---CHOOSE X3
 349
 350       X3MAX=1-AW
 351
 352       X(3)=X3MIN*X3MAX/(X3MIN+(X3MAX-X3MIN)*HWR())
 353
 354 C--CC, QQ AND RR ARE THE VARIABLE DEFINED IN OUR PAPER
 355
 356 C--IN ORDER TO SOLVE THE CUBIC EQUATION
 357
 358       CC=(1-AW)**2/4
 359
 360       QQ=(AW**2-4*(1-X(3))*(2-CC-X(3))-2*AW*(3+2*X(3)))/3
 361
 362      &     -((3+2*AW-4*X(3))**2)/9
 363
 364       RR=((3+2*AW-4*X(3))*(AW**2-4*(1-X(3))*(2-CC-X(3))
 365
 366      &     -2*AW*(3+2*X(3)))-3*(AW*(4-AW)*(2-CC)+(1-CC)
 367
 368      &     *(2*(1-X(3))-AW)**2))/6-(ONE/27)*(3+2*AW-4*X(3))**3
 369
 370 C---CHOOSE X1
 371
 372       X1MAX=2*(-QQ**3)**(ONE/6)*COS(ACOS(RR/SQRT(-QQ**3))/3)
 373
 374      &     -(3+2*AW-4*X(3))/3
 375
 376       X1MIN=1-X(3)+(AW*X(3))/(1-X(3))
 377
 378       IF (X1MAX.GE.1.OR.X1MIN.GE.1.OR.X1MAX.LE.X1MIN) GOTO 100
 379
 380       X(1)=1-(1-X1MAX)*((1-X1MIN)/(1-X1MAX))**HWR()
 381
 382 C---CALCULATE WEIGHT
 383
 384       W=((1+1/AW-2*AW)*((1-AW)*X(3)-(1-X(1))*(1-X(3))-X(3)**2)
 385
 386      &     +(1+1/(2*AW))*X(3)*(X(1)+X(3)-1)**2+2*X(3)**2*(1-X(1)))
 387
 388      &     *(1/X3MIN-1/X3MAX)*LOG((1-X1MIN)/(1-X1MAX))
 389
 390 C---QSCALE=DURHAM-LIKE TRANSVERSE MOMENTUM OF THE GLUON
 391
 392       QSCALE=EM*HWUSQR(X(3)*(1-X(1))/(2-X(1)-X(3)-AW))
 393
 394 C---FACTOR FOR GLUON EMISSION
 395
 396       ID=IDHW(JDAHEP(2,ICMF))
 397
 398       GLUFAC=0
 399
 400       IF (QSCALE.GT.HWBVMC(13)) GLUFAC=CFFAC*HWUALF(1,QSCALE)
 401
 402      &     /(PIFAC*(1-AW)*(1-2*AW+1/AW))
 403
 404 C---IN FRACTION GLUFAC*W OF EVENTS ADD A GLUON
 405
 406       IF (GLUFAC*W.GT.HWR()) THEN
 407
 408         ID3=13
 409
 410       ELSE
 411
 412         RETURN
 413
 414       ENDIF
 415
 416 C---CHECK INFRA-RED CUT-OFF FOR GLUON
 417
 418       M(1)=PHEP(5,JDAHEP(1,ICMF))
 419
 420       M(2)=HWBVMC(ID)
 421
 422       M(3)=HWBVMC(ID3)
 423
 424       E(1)=HALF*EM*(X(1)+AW+(-M(2)**2-M(3)**2)/EM**2)
 425
 426       E(3)=HALF*EM*X(3)
 427
 428       E(2)=EM-E(1)-E(3)
 429
 430       PTSQ=-LAMBDA(E(1)**2-M(1)**2,E(3)**2-M(3)**2,
 431
 432      &     E(2)**2-M(2)**2)
 433
 434       IF (PTSQ.LE.0.OR.E(1).LE.M(1).OR.E(2).LE.M(2).OR.E(3).LE.M(3))
 435
 436      $     RETURN
 437
 438 C---CALCULATE MASS-DEPENDENT SUPPRESSION
 439
 440       EPS=(RMASS(ID)/EM)**2
 441
 442       EPG=(RMASS(ID3)/EM)**2
 443
 444       GAMDEP=(1-AW)*(1+1/AW-2*AW)/(SQRT(1+AW**2+EPS**2
 445
 446      &     -2*AW-2*EPS-2*AW*EPS)*(1+EPS+(1-EPS)**2/AW-2*AW))
 447
 448       MASDEP=GAMDEP/(1-X(1))*((1+EPS+(1-EPS)**2/AW-2*AW)
 449
 450      &     *((1-AW+EPS)*X(3)*(1-X(1))-(1-X(1))**2*(1-X(3))
 451
 452      &     -X(3)**2*(1-X(1)+EPS))+(1+(1+EPS)/(2*AW))*X(3)
 453
 454      &     *(1-X(1))*(X(1)+X(3)-1)**2+2*X(3)**2*(1-X(1))**2)
 455
 456       IF (MASDEP.LT.HWR()*((1+1/AW-2*AW)*((1-AW)*X(3)
 457
 458      &     -(1-X(1))*(1-X(3))-X(3)**2)+(1+1/(2*AW))*X(3)
 459
 460      &     *(X(1)+X(3)-1)**2+2*X(3)**2*(1-X(1)))) RETURN
 461
 462 C---STORE OLD MOMENTA
 463
 464 c---PT = TOP MOMENTUM, PW= W MOMENTUM
 465
 466       CALL HWVEQU(5,PHEP(1,ICMF),PT)
 467
 468       CALL HWVEQU(5,PHEP(1,JDAHEP(1,ICMF)),PW)
 469
 470 C--------GET THE NON-EMITTING PARTON CMF DIRECTION
 471
 472       CALL HWULOF(PHEP(1,ICMF),PW,PW)
 473
 474       CALL HWRAZM(ONE,CS,SN)
 475
 476       CALL HWUROT(PW,CS,SN,R)
 477
 478       CALL HWUROF(R,PW,PW)
 479
 480       CALL HWUMAS(PW)
 481
 482 C---REORDER ENTRIES: IHEP=EMITTER,  KHEP=EMITTED
 483
 484       NHEP=NHEP+1
 485
 486       IHEP=JDAHEP(2,ICMF)
 487
 488       WHEP=JDAHEP(1,ICMF)
 489
 490       KHEP=NHEP
 491
 492 C---SET UP MOMENTA IN TOP REST FRAME
 493
 494       PHEP(1,ICMF)=0
 495
 496       PHEP(2,ICMF)=0
 497
 498       PHEP(3,ICMF)=0
 499
 500       PHEP(4,ICMF)=EM
 501
 502       PHEP(5,ICMF)=EM
 503
 504       PHEP(4,IHEP)=HALF*EM*(2-X(1)-X(3)+EPS-AW+EPG)
 505
 506       PHEP(4,KHEP)=HALF*EM*X(3)
 507
 508       PHEP(5,IHEP)=RMASS(ID)
 509
 510       PHEP(5,KHEP)=RMASS(ID3)
 511
 512       PHEP(3,KHEP)=HALF*EM*((X(1)+AW-EPS-EPG)*X(3)-2*(1+EPS-AW
 513
 514      $     -EPG-(2+EPS+EPG-AW-X(1)-X(3))))/HWUSQR((X(1)+AW
 515
 516      $     -EPS-EPG)**2-4*AW)
 517
 518       PHEP(3,IHEP)=-PHEP(3,KHEP)-HALF*EM
 519
 520      $     *HWUSQR((X(1)+AW-EPS-EPG)**2-4*AW)
 521
 522       PHEP(2,IHEP)=0
 523
 524       PHEP(1,KHEP)=HWUSQR(PHEP(4,KHEP)**2-PHEP(5,KHEP)**2
 525
 526      $     -PHEP(3,KHEP)**2)
 527
 528       PHEP(1,IHEP)=-PHEP(1,KHEP)
 529
 530       PHEP(2,KHEP)=0
 531
 532       CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,KHEP),PW1)
 533
 534       CALL HWVDIF(4,PHEP(1,ICMF),PW1,PW1)
 535
 536       CALL HWUMAS(PW1)
 537
 538       DO K=1,5
 539
 540         PHEP(K,WHEP)=PW1(K)
 541
 542       ENDDO
 543
 544 C---ORIENT IN CMF, THEN BOOST TO LAB
 545
 546       CALL HWUROB(R,PHEP(1,ICMF),PHEP(1,ICMF))
 547
 548       CALL HWUROB(R,PHEP(1,IHEP),PHEP(1,IHEP))
 549
 550       CALL HWUROB(R,PHEP(1,WHEP),PHEP(1,WHEP))
 551
 552       CALL HWUROB(R,PHEP(1,KHEP),PHEP(1,KHEP))
 553
 554       CALL HWULOB(PT,PHEP(1,IHEP),PHEP(1,IHEP))
 555
 556       CALL HWULOB(PT,PHEP(1,KHEP),PHEP(1,KHEP))
 557
 558       CALL HWULOB(PT,PHEP(1,ICMF),PHEP(1,ICMF))
 559
 560       CALL HWULOB(PT,PHEP(1,WHEP),PHEP(1,WHEP))
 561
 562 C---STATUS AND COLOUR CONNECTION
 563
 564       ISTHEP(KHEP)=114
 565
 566       IDHW(KHEP)=ID3
 567
 568       IDHEP(KHEP)=IDPDG(ID3)
 569
 570       JDAHEP(2,ICMF)=KHEP
 571
 572       JMOHEP(1,KHEP)=ICMF
 573
 574       JMOHEP(1,IHEP)=ICMF
 575
 576       JDAHEP(1,KHEP)=0
 577
 578       JMOHEP(2,IHEP)=ICMF
 579
 580       JDAHEP(2,IHEP)=KHEP
 581
 582       JMOHEP(2,KHEP)=IHEP
 583
 584       JDAHEP(2,KHEP)=ICMF
 585
 586  999  END
 587
 588 CDECK  ID>, HWBVMC.
 589
 590 *CMZ :-        -26/04/91  11.11.54  by  Bryan Webber
 591
 592 *-- Author :    Bryan Webber
 593
 594 C-----------------------------------------------------------------------
 595
 596       FUNCTION HWBVMC(ID)
 597
 598 C-----------------------------------------------------------------------
 599
 600 C     VIRTUAL MASS CUTOFF FOR PARTON TYPE ID
 601
 602 C-----------------------------------------------------------------------
 603
 604       INCLUDE 'HERWIG61.INC'
 605
 606       DOUBLE PRECISION HWBVMC
 607
 608       INTEGER ID
 609
 610       IF (ID.EQ.13) THEN
 611
 612         HWBVMC=RMASS(ID)+VGCUT
 613
 614       ELSEIF (ID.LT.13) THEN
 615
 616         HWBVMC=RMASS(ID)+VQCUT
 617
 618       ELSEIF (ID.EQ.59) THEN
 619
 620         HWBVMC=RMASS(ID)+VPCUT
 621
 622       ELSE
 623
 624         HWBVMC=RMASS(ID)
 625
 626       ENDIF
 627
 628       END