HIJING/hipyset1_35/pyrand_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE PYRAND_HIJING
   6
   7 C...Generates quantities characterizing the high-pT scattering at the
   8 C...parton level according to the matrix elements. Chooses incoming,
   9 C...reacting partons, their momentum fractions and one of the possible
  10 C...subprocesses.
  11 #include "ludat1_hijing.inc"
  12 #include "ludat2_hijing.inc"
  13 #include "pysubs_hijing.inc"
  14 #include "pypars_hijing.inc"
  15 #include "pyint1_hijing.inc"
  16 #include "pyint2_hijing.inc"
  17 #include "pyint3_hijing.inc"
  18 #include "pyint4_hijing.inc"
  19 #include "pyint5_hijing.inc"
  20
  21 C...Initial values, specifically for (first) semihard interaction.
  22       MINT(17)=0
  23       MINT(18)=0
  24       VINT(143)=1.
  25       VINT(144)=1.
  26       IF(MSUB(95).EQ.1.OR.MINT(82).GE.2) CALL PYMULT_HIJING(2)
  27       ISUB=0
  28   100 MINT(51)=0
  29
  30 C...Choice of process type - first event of overlay.
  31       IF(MINT(82).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) THEN
  32         RSUB=XSEC(0,1)*RLU_HIJING(0)
  33         DO 110 I=1,200
  34         IF(MSUB(I).NE.1) GOTO 110
  35         ISUB=I
  36         RSUB=RSUB-XSEC(I,1)
  37         IF(RSUB.LE.0.) GOTO 120
  38   110   CONTINUE
  39   120   IF(ISUB.EQ.95) ISUB=96
  40
  41 C...Choice of inclusive process type - overlayed events.
  42       ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
  43         RSUB=VINT(131)*RLU_HIJING(0)
  44         ISUB=96
  45         IF(RSUB.GT.VINT(106)) ISUB=93
  46         IF(RSUB.GT.VINT(106)+VINT(104)) ISUB=92
  47         IF(RSUB.GT.VINT(106)+VINT(104)+VINT(103)) ISUB=91
  48       ENDIF
  49       IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+1
  50       IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+1
  51       MINT(1)=ISUB
  52
  53 C...Find resonances (explicit or implicit in cross-section).
  54       MINT(72)=0
  55       KFR1=0
  56       IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
  57         KFR1=KFPR(ISUB,1)
  58       ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
  59         KFR1=25
  60       ENDIF
  61       IF(KFR1.NE.0) THEN
  62         TAUR1=PMAS(KFR1,1)**2/VINT(2)
  63         GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
  64         MINT(72)=1
  65         MINT(73)=KFR1
  66         VINT(73)=TAUR1
  67         VINT(74)=GAMR1
  68       ENDIF
  69       IF(ISUB.EQ.141) THEN
  70         KFR2=23
  71         TAUR2=PMAS(KFR2,1)**2/VINT(2)
  72         GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
  73         MINT(72)=2
  74         MINT(74)=KFR2
  75         VINT(75)=TAUR2
  76         VINT(76)=GAMR2
  77       ENDIF
  78
  79 C...Find product masses and minimum pT of process,
  80 C...optionally with broadening according to a truncated Breit-Wigner.
  81       VINT(63)=0.
  82       VINT(64)=0.
  83       MINT(71)=0
  84       VINT(71)=CKIN(3)
  85       IF(MINT(82).GE.2) VINT(71)=0.
  86       IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
  87         DO 130 I=1,2
  88         IF(KFPR(ISUB,I).EQ.0) THEN
  89         ELSEIF(MSTP(42).LE.0) THEN
  90           VINT(62+I)=PMAS(KFPR(ISUB,I),1)**2
  91         ELSE
  92           VINT(62+I)=ULMASS_HIJING(KFPR(ISUB,I))**2
  93         ENDIF
  94   130   CONTINUE
  95         IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
  96         IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
  97       ENDIF
  98
  99       IF(ISET(ISUB).EQ.0) THEN
 100 C...Double or single diffractive, or elastic scattering:
 101 C...choose m^2 according to 1/m^2 (diffractive), constant (elastic)
 102         IS=INT(1.5+RLU_HIJING(0))
 103         VINT(63)=VINT(3)**2
 104         VINT(64)=VINT(4)**2
 105         IF(ISUB.EQ.92.OR.ISUB.EQ.93) VINT(62+IS)=PARP(111)**2
 106         IF(ISUB.EQ.93) VINT(65-IS)=PARP(111)**2
 107         SH=VINT(2)
 108         SQM1=VINT(3)**2
 109         SQM2=VINT(4)**2
 110         SQM3=VINT(63)
 111         SQM4=VINT(64)
 112         SQLA12=(SH-SQM1-SQM2)**2-4.*SQM1*SQM2
 113         SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
 114         THTER1=SQM1+SQM2+SQM3+SQM4-(SQM1-SQM2)*(SQM3-SQM4)/SH-SH
 115         THTER2=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
 116         THL=0.5*(THTER1-THTER2)
 117         THU=0.5*(THTER1+THTER2)
 118         THM=MIN(MAX(THL,PARP(101)),THU)
 119         JTMAX=0
 120         IF(ISUB.EQ.92.OR.ISUB.EQ.93) JTMAX=ISUB-91
 121         DO 140 JT=1,JTMAX
 122         MINT(13+3*JT-IS*(2*JT-3))=1
 123         SQMMIN=VINT(59+3*JT-IS*(2*JT-3))
 124         SQMI=VINT(8-3*JT+IS*(2*JT-3))**2
 125         SQMJ=VINT(3*JT-1-IS*(2*JT-3))**2
 126         SQMF=VINT(68-3*JT+IS*(2*JT-3))
 127         SQUA=0.5*SH/SQMI*((1.+(SQMI-SQMJ)/SH)*THM+SQMI-SQMF-
 128      &  SQMJ**2/SH+(SQMI+SQMJ)*SQMF/SH+(SQMI-SQMJ)**2/SH**2*SQMF)
 129         QUAR=SH/SQMI*(THM*(THM+SH-SQMI-SQMJ-SQMF*(1.-(SQMI-SQMJ)/SH))+
 130      &  SQMI*SQMJ-SQMJ*SQMF*(1.+(SQMI-SQMJ-SQMF)/SH))
 131         SQMMAX=SQUA+SQRT(MAX(0.,SQUA**2-QUAR))
 132         IF(ABS(QUAR/SQUA**2).LT.1.E-06) SQMMAX=0.5*QUAR/SQUA
 133         SQMMAX=MIN(SQMMAX,(VINT(1)-SQRT(SQMF))**2)
 134         VINT(59+3*JT-IS*(2*JT-3))=SQMMIN*(SQMMAX/SQMMIN)**RLU_HIJING(0)
 135   140   CONTINUE
 136 C...Choose t-hat according to exp(B*t-hat+C*t-hat^2).
 137         SQM3=VINT(63)
 138         SQM4=VINT(64)
 139         SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
 140         THTER1=SQM1+SQM2+SQM3+SQM4-(SQM1-SQM2)*(SQM3-SQM4)/SH-SH
 141         THTER2=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
 142         THL=0.5*(THTER1-THTER2)
 143         THU=0.5*(THTER1+THTER2)
 144         B=VINT(121)
 145         C=VINT(122)
 146         IF(ISUB.EQ.92.OR.ISUB.EQ.93) THEN
 147           B=0.5*B
 148           C=0.5*C
 149         ENDIF
 150         THM=MIN(MAX(THL,PARP(101)),THU)
 151         EXPTH=0.
 152         THARG=B*(THM-THU)
 153         IF(THARG.GT.-20.) EXPTH=EXP(THARG)
 154   150   TH=THU+LOG(EXPTH+(1.-EXPTH)*RLU_HIJING(0))/B
 155         TH=MAX(THM,MIN(THU,TH))
 156         RATLOG=MIN((B+C*(TH+THM))*(TH-THM),(B+C*(TH+THU))*(TH-THU))
 157         IF(RATLOG.LT.LOG(RLU_HIJING(0))) GOTO 150
 158         VINT(21)=1.
 159         VINT(22)=0.
 160         VINT(23)=MIN(1.,MAX(-1.,(2.*TH-THTER1)/THTER2))
 161
 162 C...Note: in the following, by In is meant the integral over the
 163 C...quantity multiplying coefficient cn.
 164 C...Choose tau according to h1(tau)/tau, where
 165 C...h1(tau) = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau+tau_R) +
 166 C...I0/I3*c3*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
 167 C...I0/I4*c4*1/(tau+tau_R') +
 168 C...I0/I5*c5*tau/((s*tau-m'^2)^2+(m'*Gamma')^2), and
 169 C...c0 + c1 + c2 + c3 + c4 + c5 = 1
 170       ELSEIF(ISET(ISUB).GE.1.AND.ISET(ISUB).LE.4) THEN
 171         CALL PYKLIM_HIJING(1)
 172         IF(MINT(51).NE.0) GOTO 100
 173         RTAU=RLU_HIJING(0)
 174         MTAU=1
 175         IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
 176         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
 177         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
 178         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
 179      &  MTAU=5
 180         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
 181      &  COEF(ISUB,5)) MTAU=6
 182         CALL PYKMAP_HIJING(1,MTAU,RLU_HIJING(0))
 183
 184 C...2 -> 3, 4 processes:
 185 C...Choose tau' according to h4(tau,tau')/tau', where
 186 C...h4(tau,tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau', and
 187 C...c0 + c1 = 1.
 188         IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN
 189           CALL PYKLIM_HIJING(4)
 190           IF(MINT(51).NE.0) GOTO 100
 191           RTAUP=RLU_HIJING(0)
 192           MTAUP=1
 193           IF(RTAUP.GT.COEF(ISUB,15)) MTAUP=2
 194           CALL PYKMAP_HIJING(4,MTAUP,RLU_HIJING(0))
 195         ENDIF
 196
 197 C...Choose y* according to h2(y*), where
 198 C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
 199 C...I0/I3*c3*1/cosh(y*), I0 = y*max-y*min, and c1 + c2 + c3 = 1.
 200         CALL PYKLIM_HIJING(2)
 201         IF(MINT(51).NE.0) GOTO 100
 202         RYST=RLU_HIJING(0)
 203         MYST=1
 204         IF(RYST.GT.COEF(ISUB,7)) MYST=2
 205         IF(RYST.GT.COEF(ISUB,7)+COEF(ISUB,8)) MYST=3
 206         CALL PYKMAP_HIJING(2,MYST,RLU_HIJING(0))
 207
 208 C...2 -> 2 processes:
 209 C...Choose cos(theta-hat) (cth) according to h3(cth), where
 210 C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
 211 C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
 212 C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
 213 C...and c0 + c1 + c2 + c3 + c4 = 1.
 214         CALL PYKLIM_HIJING(3)
 215         IF(MINT(51).NE.0) GOTO 100
 216         IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
 217           RCTH=RLU_HIJING(0)
 218           MCTH=1
 219           IF(RCTH.GT.COEF(ISUB,10)) MCTH=2
 220           IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)) MCTH=3
 221           IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)+COEF(ISUB,12)) MCTH=4
 222           IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)+COEF(ISUB,12)+
 223      &    COEF(ISUB,13)) MCTH=5
 224           CALL PYKMAP_HIJING(3,MCTH,RLU_HIJING(0))
 225         ENDIF
 226
 227 C...Low-pT or multiple interactions (first semihard interaction).
 228       ELSEIF(ISET(ISUB).EQ.5) THEN
 229         CALL PYMULT_HIJING(3)
 230         ISUB=MINT(1)
 231       ENDIF
 232
 233 C...Choose azimuthal angle.
 234       VINT(24)=PARU(2)*RLU_HIJING(0)
 235
 236 C...Check against user cuts on kinematics at parton level.
 237       MINT(51)=0
 238       IF(ISUB.LE.90.OR.ISUB.GT.100) CALL PYKLIM_HIJING(0)
 239       IF(MINT(51).NE.0) GOTO 100
 240       IF(MINT(82).EQ.1.AND.MSTP(141).GE.1) THEN
 241         MCUT=0
 242         IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
 243      &  CALL PYKCUT_HIJING(MCUT)
 244         IF(MCUT.NE.0) GOTO 100
 245       ENDIF
 246
 247 C...Calculate differential cross-section for different subprocesses.
 248       CALL PYSIGH_HIJING(NCHN,SIGS)
 249
 250 C...Calculations for Monte Carlo estimate of all cross-sections.
 251       IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
 252         XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
 253       ELSEIF(MINT(82).EQ.1) THEN
 254         XSEC(ISUB,2)=XSEC(ISUB,2)+XSEC(ISUB,1)
 255       ENDIF
 256
 257 C...Multiple interactions: store results of cross-section calculation.
 258       IF(MINT(43).EQ.4.AND.MSTP(82).GE.3) THEN
 259         VINT(153)=SIGS
 260         CALL PYMULT_HIJING(4)
 261       ENDIF
 262
 263 C...Weighting using estimate of maximum of differential cross-section.
 264       VIOL=SIGS/XSEC(ISUB,1)
 265       IF(VIOL.LT.RLU_HIJING(0)) GOTO 100
 266
 267 C...Check for possible violation of estimated maximum of differential
 268 C...cross-section used in weighting.
 269       IF(MSTP(123).LE.0) THEN
 270         IF(VIOL.GT.1.) THEN
 271           WRITE(MSTU(11),1000) VIOL,NGEN(0,3)+1
 272           WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
 273           STOP
 274         ENDIF
 275       ELSEIF(MSTP(123).EQ.1) THEN
 276         IF(VIOL.GT.VINT(108)) THEN
 277           VINT(108)=VIOL
 278 C          IF(VIOL.GT.1.) THEN
 279 C            WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
 280 C            WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),
 281 C     &      VINT(26)
 282 C          ENDIF
 283         ENDIF
 284       ELSEIF(VIOL.GT.VINT(108)) THEN
 285         VINT(108)=VIOL
 286         IF(VIOL.GT.1.) THEN
 287           XDIF=XSEC(ISUB,1)*(VIOL-1.)
 288           XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
 289           IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
 290      &    XSEC(0,1)=XSEC(0,1)+XDIF
 291 C          WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
 292 C          WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
 293 C          IF(ISUB.LE.9) THEN
 294 C            WRITE(MSTU(11),1300) ISUB,XSEC(ISUB,1)
 295 C          ELSEIF(ISUB.LE.99) THEN
 296 C            WRITE(MSTU(11),1400) ISUB,XSEC(ISUB,1)
 297 C          ELSE
 298 C            WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1)
 299 C          ENDIF
 300           VINT(108)=1.
 301         ENDIF
 302       ENDIF
 303
 304 C...Multiple interactions: choose impact parameter.
 305       VINT(148)=1.
 306       IF(MINT(43).EQ.4.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3)
 307      &THEN
 308         CALL PYMULT_HIJING(5)
 309         IF(VINT(150).LT.RLU_HIJING(0)) GOTO 100
 310       ENDIF
 311       IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
 312         IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+1
 313         IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
 314       ENDIF
 315       IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
 316
 317 C...Choose flavour of reacting partons (and subprocess).
 318       RSIGS=SIGS*RLU_HIJING(0)
 319       QT2=VINT(48)
 320       RQQBAR=PARP(87)*(1.-(QT2/(QT2+(PARP(88)*PARP(82))**2))**2)
 321       IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
 322      &RLU_HIJING(0).GT.RQQBAR)) THEN
 323         DO 190 ICHN=1,NCHN
 324         KFL1=ISIG(ICHN,1)
 325         KFL2=ISIG(ICHN,2)
 326         MINT(2)=ISIG(ICHN,3)
 327         RSIGS=RSIGS-SIGH(ICHN)
 328         IF(RSIGS.LE.0.) GOTO 210
 329   190   CONTINUE
 330
 331 C...Multiple interactions: choose qqbar preferentially at small pT.
 332       ELSEIF(ISUB.EQ.96) THEN
 333         CALL PYSPLI_HIJING(MINT(11),21,KFL1,KFLDUM)
 334         CALL PYSPLI_HIJING(MINT(12),21,KFL2,KFLDUM)
 335         MINT(1)=11
 336         MINT(2)=1
 337         IF(KFL1.EQ.KFL2.AND.RLU_HIJING(0).LT.0.5) MINT(2)=2
 338
 339 C...Low-pT: choose string drawing configuration.
 340       ELSE
 341         KFL1=21
 342         KFL2=21
 343         RSIGS=6.*RLU_HIJING(0)
 344         MINT(2)=1
 345         IF(RSIGS.GT.1.) MINT(2)=2
 346         IF(RSIGS.GT.2.) MINT(2)=3
 347       ENDIF
 348
 349 C...Reassign QCD process. Partons before initial state radiation.
 350   210 IF(MINT(2).GT.10) THEN
 351         MINT(1)=MINT(2)/10
 352         MINT(2)=MOD(MINT(2),10)
 353       ENDIF
 354       MINT(15)=KFL1
 355       MINT(16)=KFL2
 356       MINT(13)=MINT(15)
 357       MINT(14)=MINT(16)
 358       VINT(141)=VINT(41)
 359       VINT(142)=VINT(42)
 360
 361 C...Format statements for differential cross-section maximum violations.
 362  1000 FORMAT(1X,'Error: maximum violated by',1P,E11.3,1X,
 363      &'in event',1X,I7,'.'/1X,'Execution stopped!')
 364  1100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
 365      &E11.3,', y* =',E11.3,', cthe = ',0P,F11.7,', tau'' =',1P,E11.3)
 366  1200 FORMAT(1X,'Warning: maximum violated by',1P,E11.3,1X,
 367      &'in event',1X,I7)
 368  1300 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,E11.3)
 369  1400 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,E11.3)
 370  1500 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,E11.3)
 371
 372       RETURN
 373       END