PHOS/shaker/lux4jt.f

   1 *CMZ :          17/07/98  15.44.35  by  Federico Carminati
   2 *-- Author :
   3 C*********************************************************************
   4
   5       SUBROUTINE LUX4JT(NJET,CUT,KFL,ECM,KFLN,X1,X2,X4,X12,X14)
   6
   7 C...Purpose: to select the kinematical variables of four-jet events.
   8 *KEEP,LUDAT1.
   9       COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  10       SAVE /LUDAT1/
  11 *KEND.
  12       DIMENSION WTA(4),WTB(4),WTC(4),WTD(4),WTE(4)
  13
  14 C...Common constants. Colour factors for QCD and Abelian gluon theory.
  15       PMQ=ULMASS(KFL)
  16       QME=(2.*PMQ/ECM)**2
  17       CT=LOG(1./CUT-5.)
  18       IF(MSTJ(109).EQ.0) THEN
  19         CF=4./3.
  20         CN=3.
  21         TR=2.5
  22       ELSE
  23         CF=1.
  24         CN=0.
  25         TR=15.
  26       ENDIF
  27
  28 C...Choice of process (qqbargg or qqbarqqbar).
  29   100 NJET=4
  30       IT=1
  31       IF(PARJ(155).GT.RLU(0)) IT=2
  32       IF(MSTJ(101).LE.-3) IT=-MSTJ(101)-2
  33       IF(IT.EQ.1) WTMX=0.7/CUT**2
  34       IF(IT.EQ.1.AND.MSTJ(109).EQ.2) WTMX=0.6/CUT**2
  35       IF(IT.EQ.2) WTMX=0.1125*CF*TR/CUT**2
  36       ID=1
  37
  38 C...Sample the five kinematical variables (for qqgg preweighted in y34).
  39   110 Y134=3.*CUT+(1.-6.*CUT)*RLU(0)
  40       Y234=3.*CUT+(1.-6.*CUT)*RLU(0)
  41       IF(IT.EQ.1) Y34=(1.-5.*CUT)*EXP(-CT*RLU(0))
  42       IF(IT.EQ.2) Y34=CUT+(1.-6.*CUT)*RLU(0)
  43       IF(Y34.LE.Y134+Y234-1..OR.Y34.GE.Y134*Y234) GOTO 110
  44       VT=RLU(0)
  45       CP=COS(PARU(1)*RLU(0))
  46       Y14=(Y134-Y34)*VT
  47       Y13=Y134-Y14-Y34
  48       VB=Y34*(1.-Y134-Y234+Y34)/((Y134-Y34)*(Y234-Y34))
  49       Y24=0.5*(Y234-Y34)*(1.-4.*SQRT(MAX(0.,VT*(1.-VT)*VB*(1.-VB)))*
  50      &CP-(1.-2.*VT)*(1.-2.*VB))
  51       Y23=Y234-Y34-Y24
  52       Y12=1.-Y134-Y23-Y24
  53       IF(MIN(Y12,Y13,Y14,Y23,Y24).LE.CUT) GOTO 110
  54       Y123=Y12+Y13+Y23
  55       Y124=Y12+Y14+Y24
  56
  57 C...Calculate matrix elements for qqgg or qqqq process.
  58       IC=0
  59       WTTOT=0.
  60   120 IC=IC+1
  61       IF(IT.EQ.1) THEN
  62         WTA(IC)=(Y12*Y34**2-Y13*Y24*Y34+Y14*Y23*Y34+3.*Y12*Y23*Y34+
  63      &  3.*Y12*Y14*Y34+4.*Y12**2*Y34-Y13*Y23*Y24+2.*Y12*Y23*Y24-
  64      &  Y13*Y14*Y24-2.*Y12*Y13*Y24+2.*Y12**2*Y24+Y14*Y23**2+2.*Y12*
  65      &  Y23**2+Y14**2*Y23+4.*Y12*Y14*Y23+4.*Y12**2*Y23+2.*Y12*Y14**2+
  66      &  2.*Y12*Y13*Y14+4.*Y12**2*Y14+2.*Y12**2*Y13+2.*Y12**3)/(2.*Y13*
  67      &  Y134*Y234*Y24)+(Y24*Y34+Y12*Y34+Y13*Y24-Y14*Y23+Y12*Y13)/(Y13*
  68      &  Y134**2)+2.*Y23*(1.-Y13)/(Y13*Y134*Y24)+Y34/(2.*Y13*Y24)
  69         WTB(IC)=(Y12*Y24*Y34+Y12*Y14*Y34-Y13*Y24**2+Y13*Y14*Y24+2.*Y12*
  70      &  Y14*Y24)/(Y13*Y134*Y23*Y14)+Y12*(1.+Y34)*Y124/(Y134*Y234*Y14*
  71      &  Y24)-(2.*Y13*Y24+Y14**2+Y13*Y23+2.*Y12*Y13)/(Y13*Y134*Y14)+
  72      &  Y12*Y123*Y124/(2.*Y13*Y14*Y23*Y24)
  73         WTC(IC)=-(5.*Y12*Y34**2+2.*Y12*Y24*Y34+2.*Y12*Y23*Y34+2.*Y12*
  74      &  Y14*Y34+2.*Y12*Y13*Y34+4.*Y12**2*Y34-Y13*Y24**2+Y14*Y23*Y24+
  75      &  Y13*Y23*Y24+Y13*Y14*Y24-Y12*Y14*Y24-Y13**2*Y24-3.*Y12*Y13*Y24-
  76      &  Y14*Y23**2-Y14**2*Y23+Y13*Y14*Y23-3.*Y12*Y14*Y23-Y12*Y13*Y23)/
  77      &  (4.*Y134*Y234*Y34**2)+(3.*Y12*Y34**2-3.*Y13*Y24*Y34+3.*Y12*Y24*
  78      &  Y34+3.*Y14*Y23*Y34-Y13*Y24**2-Y12*Y23*Y34+6.*Y12*Y14*Y34+2.*Y12*
  79      &  Y13*Y34-2.*Y12**2*Y34+Y14*Y23*Y24-3.*Y13*Y23*Y24-2.*Y13*Y14*
  80      &  Y24+4.*Y12*Y14*Y24+2.*Y12*Y13*Y24+3.*Y14*Y23**2+2.*Y14**2*Y23+
  81      &  2.*Y14**2*Y12+2.*Y12**2*Y14+6.*Y12*Y14*Y23-2.*Y12*Y13**2-
  82      &  2.*Y12**2*Y13)/(4.*Y13*Y134*Y234*Y34)
  83         WTC(IC)=WTC(IC)+(2.*Y12*Y34**2-2.*Y13*Y24*Y34+Y12*Y24*Y34+
  84      &  4.*Y13*Y23*Y34+4.*Y12*Y14*Y34+2.*Y12*Y13*Y34+2.*Y12**2*Y34-
  85      &  Y13*Y24**2+3.*Y14*Y23*Y24+4.*Y13*Y23*Y24-2.*Y13*Y14*Y24+
  86      &  4.*Y12*Y14*Y24+2.*Y12*Y13*Y24+2.*Y14*Y23**2+4.*Y13*Y23**2+
  87      &  2.*Y13*Y14*Y23+2.*Y12*Y14*Y23+4.*Y12*Y13*Y23+2.*Y12*Y14**2+4.*
  88      &  Y12**2*Y13+4.*Y12*Y13*Y14+2.*Y12**2*Y14)/(4.*Y13*Y134*Y24*Y34)-
  89      &  (Y12*Y34**2-2.*Y14*Y24*Y34-2.*Y13*Y24*Y34-Y14*Y23*Y34+Y13*Y23*
  90      &  Y34+Y12*Y14*Y34+2.*Y12*Y13*Y34-2.*Y14**2*Y24-4.*Y13*Y14*Y24-
  91      &  4.*Y13**2*Y24-Y14**2*Y23-Y13**2*Y23+Y12*Y13*Y14-Y12*Y13**2)/
  92      &  (2.*Y13*Y34*Y134**2)+(Y12*Y34**2-4.*Y14*Y24*Y34-2.*Y13*Y24*Y34-
  93      &  2.*Y14*Y23*Y34-4.*Y13*Y23*Y34-4.*Y12*Y14*Y34-4.*Y12*Y13*Y34-
  94      &  2.*Y13*Y14*Y24+2.*Y13**2*Y24+2.*Y14**2*Y23-2.*Y13*Y14*Y23-
  95      &  Y12*Y14**2-6.*Y12*Y13*Y14-Y12*Y13**2)/(4.*Y34**2*Y134**2)
  96         WTTOT=WTTOT+Y34*CF*(CF*WTA(IC)+(CF-0.5*CN)*WTB(IC)+CN*WTC(IC))/
  97      &  8.
  98       ELSE
  99         WTD(IC)=(Y13*Y23*Y34+Y12*Y23*Y34-Y12**2*Y34+Y13*Y23*Y24+2.*Y12*
 100      &  Y23*Y24-Y14*Y23**2+Y12*Y13*Y24+Y12*Y14*Y23+Y12*Y13*Y14)/(Y13**2*
 101      &  Y123**2)-(Y12*Y34**2-Y13*Y24*Y34+Y12*Y24*Y34-Y14*Y23*Y34-Y12*
 102      &  Y23*Y34-Y13*Y24**2+Y14*Y23*Y24-Y13*Y23*Y24-Y13**2*Y24+Y14*
 103      &  Y23**2)/(Y13**2*Y123*Y134)+(Y13*Y14*Y12+Y34*Y14*Y12-Y34**2*Y12+
 104      &  Y13*Y14*Y24+2.*Y34*Y14*Y24-Y23*Y14**2+Y34*Y13*Y24+Y34*Y23*Y14+
 105      &  Y34*Y13*Y23)/(Y13**2*Y134**2)-(Y34*Y12**2-Y13*Y24*Y12+Y34*Y24*
 106      &  Y12-Y23*Y14*Y12-Y34*Y14*Y12-Y13*Y24**2+Y23*Y14*Y24-Y13*Y14*Y24-
 107      &  Y13**2*Y24+Y23*Y14**2)/(Y13**2*Y134*Y123)
 108         WTE(IC)=(Y12*Y34*(Y23-Y24+Y14+Y13)+Y13*Y24**2-Y14*Y23*Y24+Y13*
 109      &  Y23*Y24+Y13*Y14*Y24+Y13**2*Y24-Y14*Y23*(Y14+Y23+Y13))/(Y13*Y23*
 110      &  Y123*Y134)-Y12*(Y12*Y34-Y23*Y24-Y13*Y24-Y14*Y23-Y14*Y13)/(Y13*
 111      &  Y23*Y123**2)-(Y14+Y13)*(Y24+Y23)*Y34/(Y13*Y23*Y134*Y234)+
 112      &  (Y12*Y34*(Y14-Y24+Y23+Y13)+Y13*Y24**2-Y23*Y14*Y24+Y13*Y14*Y24+
 113      &  Y13*Y23*Y24+Y13**2*Y24-Y23*Y14*(Y14+Y23+Y13))/(Y13*Y14*Y134*
 114      &  Y123)-Y34*(Y34*Y12-Y14*Y24-Y13*Y24-Y23*Y14-Y23*Y13)/(Y13*Y14*
 115      &  Y134**2)-(Y23+Y13)*(Y24+Y14)*Y12/(Y13*Y14*Y123*Y124)
 116         WTTOT=WTTOT+CF*(TR*WTD(IC)+(CF-0.5*CN)*WTE(IC))/16.
 117       ENDIF
 118
 119 C...Permutations of momenta in matrix element. Weighting.
 120   130 IF(IC.EQ.1.OR.IC.EQ.3.OR.ID.EQ.2.OR.ID.EQ.3) THEN
 121         YSAV=Y13
 122         Y13=Y14
 123         Y14=YSAV
 124         YSAV=Y23
 125         Y23=Y24
 126         Y24=YSAV
 127         YSAV=Y123
 128         Y123=Y124
 129         Y124=YSAV
 130       ENDIF
 131       IF(IC.EQ.2.OR.IC.EQ.4.OR.ID.EQ.3.OR.ID.EQ.4) THEN
 132         YSAV=Y13
 133         Y13=Y23
 134         Y23=YSAV
 135         YSAV=Y14
 136         Y14=Y24
 137         Y24=YSAV
 138         YSAV=Y134
 139         Y134=Y234
 140         Y234=YSAV
 141       ENDIF
 142       IF(IC.LE.3) GOTO 120
 143       IF(ID.EQ.1.AND.WTTOT.LT.RLU(0)*WTMX) GOTO 110
 144       IC=5
 145
 146 C...qqgg events: string configuration and event type.
 147       IF(IT.EQ.1) THEN
 148         IF(MSTJ(109).EQ.0.AND.ID.EQ.1) THEN
 149           PARJ(156)=Y34*(2.*(WTA(1)+WTA(2)+WTA(3)+WTA(4))+4.*(WTC(1)+
 150      &    WTC(2)+WTC(3)+WTC(4)))/(9.*WTTOT)
 151           IF(WTA(2)+WTA(4)+2.*(WTC(2)+WTC(4)).GT.RLU(0)*(WTA(1)+WTA(2)+
 152      &    WTA(3)+WTA(4)+2.*(WTC(1)+WTC(2)+WTC(3)+WTC(4)))) ID=2
 153           IF(ID.EQ.2) GOTO 130
 154         ELSEIF(MSTJ(109).EQ.2.AND.ID.EQ.1) THEN
 155           PARJ(156)=Y34*(WTA(1)+WTA(2)+WTA(3)+WTA(4))/(8.*WTTOT)
 156           IF(WTA(2)+WTA(4).GT.RLU(0)*(WTA(1)+WTA(2)+WTA(3)+WTA(4))) ID=2
 157           IF(ID.EQ.2) GOTO 130
 158         ENDIF
 159         MSTJ(120)=3
 160         IF(MSTJ(109).EQ.0.AND.0.5*Y34*(WTC(1)+WTC(2)+WTC(3)+WTC(4)).GT.
 161      &  RLU(0)*WTTOT) MSTJ(120)=4
 162         KFLN=21
 163
 164 C...Mass cuts. Kinematical variables out.
 165         IF(Y12.LE.CUT+QME) NJET=2
 166         IF(NJET.EQ.2) GOTO 150
 167         Q12=0.5*(1.-SQRT(1.-QME/Y12))
 168         X1=1.-(1.-Q12)*Y234-Q12*Y134
 169         X4=1.-(1.-Q12)*Y134-Q12*Y234
 170         X2=1.-Y124
 171         X12=(1.-Q12)*Y13+Q12*Y23
 172         X14=Y12-0.5*QME
 173         IF(Y134*Y234/((1.-X1)*(1.-X4)).LE.RLU(0)) NJET=2
 174
 175 C...qqbarqqbar events: string configuration, choose new flavour.
 176       ELSE
 177         IF(ID.EQ.1) THEN
 178           WTR=RLU(0)*(WTD(1)+WTD(2)+WTD(3)+WTD(4))
 179           IF(WTR.LT.WTD(2)+WTD(3)+WTD(4)) ID=2
 180           IF(WTR.LT.WTD(3)+WTD(4)) ID=3
 181           IF(WTR.LT.WTD(4)) ID=4
 182           IF(ID.GE.2) GOTO 130
 183         ENDIF
 184         MSTJ(120)=5
 185         PARJ(156)=CF*TR*(WTD(1)+WTD(2)+WTD(3)+WTD(4))/(16.*WTTOT)
 186   140   KFLN=1+INT(5.*RLU(0))
 187         IF(KFLN.NE.KFL.AND.0.2*PARJ(156).LE.RLU(0)) GOTO 140
 188         IF(KFLN.EQ.KFL.AND.1.-0.8*PARJ(156).LE.RLU(0)) GOTO 140
 189         IF(KFLN.GT.MSTJ(104)) NJET=2
 190         PMQN=ULMASS(KFLN)
 191         QMEN=(2.*PMQN/ECM)**2
 192
 193 C...Mass cuts. Kinematical variables out.
 194         IF(Y24.LE.CUT+QME.OR.Y13.LE.1.1*QMEN) NJET=2
 195         IF(NJET.EQ.2) GOTO 150
 196         Q24=0.5*(1.-SQRT(1.-QME/Y24))
 197         Q13=0.5*(1.-SQRT(1.-QMEN/Y13))
 198         X1=1.-(1.-Q24)*Y123-Q24*Y134
 199         X4=1.-(1.-Q24)*Y134-Q24*Y123
 200         X2=1.-(1.-Q13)*Y234-Q13*Y124
 201         X12=(1.-Q24)*((1.-Q13)*Y14+Q13*Y34)+Q24*((1.-Q13)*Y12+Q13*Y23)
 202         X14=Y24-0.5*QME
 203         X34=(1.-Q24)*((1.-Q13)*Y23+Q13*Y12)+Q24*((1.-Q13)*Y34+Q13*Y14)
 204         IF(PMQ**2+PMQN**2+MIN(X12,X34)*ECM**2.LE.
 205      &  (PARJ(127)+PMQ+PMQN)**2) NJET=2
 206         IF(Y123*Y134/((1.-X1)*(1.-X4)).LE.RLU(0)) NJET=2
 207       ENDIF
 208   150 IF(MSTJ(101).LE.-2.AND.NJET.EQ.2) GOTO 100
 209
 210       RETURN
 211       END