HERWIG/src/hwccut.f

   1
   2 CDECK  ID>, HWCCUT.
   3
   4 *CMZ :-        -26/04/91  14.29.39  by  Federico Carminati
   5
   6 *-- Author :    Bryan Webber
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWCCUT(JHEP,KHEP,PCL,BTCLUS,SPLIT)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     Cuts into 2 the cluster, momentum PCL, made of partons JHEP & KHEP
  15
  16 C-----------------------------------------------------------------------
  17
  18       INCLUDE 'HERWIG61.INC'
  19
  20       DOUBLE PRECISION HWREXQ,HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,EMX,EMY,
  21
  22      & QM3,PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5),SKAPPA,DELTM,
  23
  24      & VSCA,VTMP(4),RKAPPA,VCLUS
  25
  26       INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,ID1,ID2,ID3,NTRY,NTRYMX,J,IB
  27
  28       LOGICAL BTCLUS,SPLIT
  29
  30       EXTERNAL HWREXQ,HWUPCM,HWR,HWVDOT,HWRINT
  31
  32       COMMON/HWCFRM/VCLUS(4,NMXHEP)
  33
  34       PARAMETER (SKAPPA=1.,NTRYMX=100)
  35
  36       IF (IERROR.NE.0) RETURN
  37
  38       EMC=PCL(5)
  39
  40       ID1=IDHW(JHEP)
  41
  42       ID2=IDHW(KHEP)
  43
  44       QM1=RMASS(ID1)
  45
  46       QM2=RMASS(ID2)
  47
  48       SPLIT=.FALSE.
  49
  50       NTRY=0
  51
  52 C Decide if cluster contains a b-(anti)quark
  53
  54       IF (ID1.EQ.5.OR.ID1.EQ.11.OR.ID2.EQ.5.OR.ID2.EQ.11) THEN
  55
  56         IB=2
  57
  58       ELSE
  59
  60         IB=1
  61
  62       ENDIF
  63
  64       IF (BTCLUS) THEN
  65
  66 C Split beam and target clusters as soft clusters
  67
  68 C Both (remnant) children treated like soft clusters if IOPREM=0(1)
  69
  70   10    ID3=HWRINT(1,2)
  71
  72         QM3=RMASS(ID3)
  73
  74         IF (EMC.LE.QM1+QM2+2.*QM3) THEN
  75
  76           ID3=3-ID3
  77
  78           QM3=RMASS(ID3)
  79
  80           IF (EMC.LE.QM1+QM2+2.*QM3) RETURN
  81
  82         ENDIF
  83
  84         PXY=EMC-QM1-QM2-TWO*QM3
  85
  86         IF (ISTHEP(JHEP).EQ.153.OR.ISTHEP(JHEP).EQ.154.OR.
  87
  88      &      IOPREM.EQ.0) THEN
  89
  90           EMX=QM1+QM3+HWREXQ(BTCLM,PXY)
  91
  92         ELSE
  93
  94           EMX=QM1+QM3+PXY*HWR()**PSPLT(IB)
  95
  96         ENDIF
  97
  98         IF (ISTHEP(KHEP).EQ.153.OR.ISTHEP(KHEP).EQ.154.OR.
  99
 100      &      IOPREM.EQ.0) THEN
 101
 102           EMY=QM2+QM3+HWREXQ(BTCLM,PXY)
 103
 104         ELSE
 105
 106           EMY=QM2+QM3+PXY*HWR()**PSPLT(IB)
 107
 108         ENDIF
 109
 110         IF (EMX+EMY.GE.EMC) THEN
 111
 112           NTRY=NTRY+1
 113
 114           IF (NTRY.GT.NTRYMX) RETURN
 115
 116           GOTO 10
 117
 118         ENDIF
 119
 120         PCX=HWUPCM(EMX,QM1,QM3)
 121
 122         PCY=HWUPCM(EMY,QM2,QM3)
 123
 124       ELSE
 125
 126 C Choose fragment masses for ordinary cluster
 127
 128         PXY=EMC-QM1-QM2
 129
 130   20    NTRY=NTRY+1
 131
 132         IF (NTRY.GT.NTRYMX) RETURN
 133
 134   30    EMX=QM1+PXY*HWR()**PSPLT(IB)
 135
 136         EMY=QM2+PXY*HWR()**PSPLT(IB)
 137
 138         IF (EMX+EMY.GE.EMC) GOTO 30
 139
 140 C u,d,s pair production with weights QWT
 141
 142   40    ID3=HWRINT(1,3)
 143
 144         IF (QWT(ID3).LT.HWR()) GOTO 40
 145
 146         QM3=RMASS(ID3)
 147
 148         PCX=HWUPCM(EMX,QM1,QM3)
 149
 150         IF (PCX.LT.ZERO) GOTO 20
 151
 152         PCY=HWUPCM(EMY,QM2,QM3)
 153
 154         IF (PCY.LT.ZERO) GOTO 20
 155
 156         SPLIT=.TRUE.
 157
 158       ENDIF
 159
 160 C Boost antiquark to CoM frame to find axis
 161
 162       CALL HWULOF(PCL,PHEP(1,KHEP),AX)
 163
 164       RCM=1./SQRT(HWVDOT(3,AX,AX))
 165
 166       CALL HWVSCA(3,RCM,AX,AX)
 167
 168 C Construct new CoM momenta (collinear)
 169
 170       PXY=HWUPCM(EMC,EMX,EMY)
 171
 172       CALL HWVSCA(3,PXY,AX,PC)
 173
 174       PC(4)=SQRT(PXY**2+EMY**2)
 175
 176       PC(5)=EMY
 177
 178       CALL HWVSCA(3,PCY,AX,PA)
 179
 180       PA(4)=SQRT(PCY**2+QM2**2)
 181
 182       PA(5)=QM2
 183
 184       CALL HWULOB(PC,PA,PB)
 185
 186       CALL HWVDIF(4,PC,PB,PA)
 187
 188       PA(5)=QM3
 189
 190       LHEP=NHEP+1
 191
 192       MHEP=NHEP+2
 193
 194       CALL HWULOB(PCL,PB,PHEP(1,KHEP))
 195
 196       CALL HWULOB(PCL,PA,PHEP(1,MHEP))
 197
 198       CALL HWVSCA(3,-ONE,PC,PC)
 199
 200       PC(4)=EMC-PC(4)
 201
 202       PC(5)=EMX
 203
 204       CALL HWVSCA(3,PCX,AX,PA)
 205
 206       PA(4)=SQRT(PCX**2+QM3**2)
 207
 208       CALL HWULOB(PC,PA,PB)
 209
 210       CALL HWULOB(PCL,PB,PHEP(1,LHEP))
 211
 212       DO 50 J=1,4
 213
 214   50  PHEP(J,JHEP)=PCL(J)-PHEP(J,KHEP)-PHEP(J,LHEP)-PHEP(J,MHEP)
 215
 216       PHEP(5,JHEP)=QM1
 217
 218       CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))
 219
 220 C Construct new vertex positions
 221
 222       RKAPPA=GEV2MM/SKAPPA
 223
 224       CALL HWVSCA(3,RKAPPA,AX,AX)
 225
 226       DELTM=(EMX-EMY)*(EMX+EMY)/(TWO*EMC)
 227
 228       CALL HWVSCA(3,DELTM,AX,VTMP)
 229
 230       VTMP(4)=(HALF*EMC-PXY)*RKAPPA
 231
 232       CALL HWULB4(PCL,VTMP,VTMP)
 233
 234       CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VHEP(1,LHEP))
 235
 236       CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))
 237
 238       VSCA=0.25*EMC+HALF*(PXY+DELTM)
 239
 240       CALL HWVSCA(3,VSCA,AX,VTMP)
 241
 242       VTMP(4)=(EMC-VSCA)*RKAPPA
 243
 244       CALL HWULB4(PCL,VTMP,VTMP)
 245
 246       CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,MHEP))
 247
 248       VSCA=-0.25*EMC+HALF*(DELTM-PXY)
 249
 250       CALL HWVSCA(3,VSCA,AX,VTMP)
 251
 252       VTMP(4)=(EMC+VSCA)*RKAPPA
 253
 254       CALL HWULB4(PCL,VTMP,VTMP)
 255
 256       CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,JHEP))
 257
 258 C (Re-)label quarks
 259
 260       IDHW(LHEP)=ID3+6
 261
 262       IDHW(MHEP)=ID3
 263
 264       IDHEP(MHEP)= IDPDG(ID3)
 265
 266       IDHEP(LHEP)=-IDPDG(ID3)
 267
 268       ISTHEP(LHEP)=151
 269
 270       ISTHEP(MHEP)=151
 271
 272       JMOHEP(2,JHEP)=LHEP
 273
 274       JDAHEP(2,KHEP)=MHEP
 275
 276       JMOHEP(1,LHEP)=JMOHEP(1,KHEP)
 277
 278       JMOHEP(2,LHEP)=MHEP
 279
 280       JDAHEP(1,LHEP)=0
 281
 282       JDAHEP(2,LHEP)=JHEP
 283
 284       JMOHEP(1,MHEP)=JMOHEP(1,JHEP)
 285
 286       JMOHEP(2,MHEP)=KHEP
 287
 288       JDAHEP(1,MHEP)=0
 289
 290       JDAHEP(2,MHEP)=LHEP
 291
 292       NHEP=NHEP+2
 293
 294   999 END
 295
 296 CDECK  ID>, HWCDEC.
 297
 298 *CMZ :-        -26/04/91  10.18.56  by  Bryan Webber
 299
 300 *-- Author :    Bryan Webber
 301
 302 C-----------------------------------------------------------------------
 303
 304       SUBROUTINE HWCDEC
 305
 306 C-----------------------------------------------------------------------
 307
 308 C     DECAYS CLUSTERS INTO PRIMARY HADRONS
 309
 310 C-----------------------------------------------------------------------
 311
 312       INCLUDE 'HERWIG61.INC'
 313
 314       INTEGER JCL,KCL,IP,JP,KP,IST,ID1,ID2,ID3
 315
 316       IF (IERROR.NE.0) RETURN
 317
 318       IF (IPROC/1000.EQ.9.OR.IPROC/1000.EQ.5) THEN
 319
 320 C---RELABEL CLUSTER CONNECTED TO REMNANT IN DIS
 321
 322         DO 10 JCL=2,NHEP
 323
 324         IF (ISTHEP(JCL).EQ.164) GOTO 20
 325
 326         IF (ISTHEP(JCL).EQ.165) THEN
 327
 328           IP=JMOHEP(1,JCL)
 329
 330           JP=JMOHEP(2,JCL)
 331
 332           KP=IP
 333
 334           IF (ISTHEP(IP).EQ.162) THEN
 335
 336             KP=JP
 337
 338             JP=IP
 339
 340           ENDIF
 341
 342           IF (JMOHEP(2,KP).NE.JP) THEN
 343
 344             IP=JMOHEP(2,KP)
 345
 346           ELSE
 347
 348             IP=JDAHEP(2,KP)
 349
 350           ENDIF
 351
 352           KCL=JDAHEP(1,IP)
 353
 354           IF (ISTHEP(KCL)/10.NE.16) CALL HWWARN('HWCDEC',100,*999)
 355
 356           ISTHEP(KCL)=164
 357
 358           GOTO 20
 359
 360         ENDIF
 361
 362    10   CONTINUE
 363
 364       ENDIF
 365
 366    20 CONTINUE
 367
 368       DO 30 JCL=1,NHEP
 369
 370       IST=ISTHEP(JCL)
 371
 372       IF (IST.GT.162.AND.IST.LT.166) THEN
 373
 374 C---DON'T HADRONIZE BEAM/TARGET CLUSTERS
 375
 376         IF (IST.EQ.163.OR..NOT.GENSOF) THEN
 377
 378 C---SET UP FLAVOURS FOR CLUSTER DECAY
 379
 380           CALL HWCFLA(IDHW(JMOHEP(1,JCL)),IDHW(JMOHEP(2,JCL)),ID1,ID3)
 381
 382           CALL HWCHAD(JCL,ID1,ID3,ID2)
 383
 384         ENDIF
 385
 386       ENDIF
 387
 388    30 CONTINUE
 389
 390       ISTAT=50
 391
 392   999 END
 393
 394 CDECK  ID>, HWCFLA.
 395
 396 *CMZ :-        -26/04/91  10.18.56  by  Bryan Webber
 397
 398 *-- Author :    Bryan Webber
 399
 400 C-----------------------------------------------------------------------
 401
 402       SUBROUTINE HWCFLA(JD1,JD2,ID1,ID2)
 403
 404 C-----------------------------------------------------------------------
 405
 406 C     SETS UP FLAVOURS FOR CLUSTER DECAY
 407
 408 C-----------------------------------------------------------------------
 409
 410       INTEGER JD1,JD2,ID1,ID2,JD,JDEC(12)
 411
 412       DATA JDEC/1,2,3,10,11,12,4,5,6,7,8,9/
 413
 414       JD=JD1
 415
 416       IF (JD.GT.12) JD=JD-108
 417
 418       ID1=JDEC(JD)
 419
 420       JD=JD2
 421
 422       IF (JD.GT.12) JD=JD-96
 423
 424       ID2=JDEC(JD-6)
 425
 426       END