HIJING/hipyset1_35/pysspa_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE PYSSPA_HIJING(IPU1,IPU2)
   6
   7 C...Generates spacelike parton showers.
   8       IMPLICIT DOUBLE PRECISION(D)
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12 #include "pysubs_hijing.inc"
  13 #include "pypars_hijing.inc"
  14 #include "pyint1_hijing.inc"
  15 #include "pyint2_hijing.inc"
  16 #include "pyint3_hijing.inc"
  17       DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVS(2),ROBO(5),
  18      &XFS(2,-6:6),XFA(-6:6),XFB(-6:6),XFN(-6:6),WTAP(-6:6),WTSF(-6:6),
  19      &THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),DPB(4)
  20
  21 C...Calculate maximum virtuality and check that evolution possible.
  22       IPUS1=IPU1
  23       IPUS2=IPU2
  24       ISUB=MINT(1)
  25       Q2E=VINT(52)
  26       IF(ISET(ISUB).EQ.1) THEN
  27         Q2E=Q2E/PARP(67)
  28       ELSEIF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN
  29         Q2E=PMAS(23,1)**2
  30         IF(ISUB.EQ.8.OR.ISUB.EQ.76.OR.ISUB.EQ.77) Q2E=PMAS(24,1)**2
  31       ENDIF
  32       TMAX=LOG(PARP(67)*PARP(63)*Q2E/PARP(61)**2)
  33       IF(PARP(67)*Q2E.LT.MAX(PARP(62)**2,2.*PARP(61)**2).OR.
  34      &TMAX.LT.0.2) RETURN
  35
  36 C...Common constants and initial values. Save normal Lambda value.
  37       XE0=2.*PARP(65)/VINT(1)
  38       ALAMS=PARU(111)
  39       PARU(111)=PARP(61)
  40       NS=N
  41   100 N=NS
  42       DO 110 JT=1,2
  43       KFLS(JT)=MINT(14+JT)
  44       KFLS(JT+2)=KFLS(JT)
  45       XS(JT)=VINT(40+JT)
  46       ZS(JT)=1.
  47       Q2S(JT)=PARP(67)*Q2E
  48       TEVS(JT)=TMAX
  49       ALAM(JT)=PARP(61)
  50       THE2(JT)=100.
  51       DO 110 KFL=-6,6
  52   110 XFS(JT,KFL)=XSFX(JT,KFL)
  53       DSH=VINT(44)
  54       IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) DSH=VINT(26)*VINT(2)
  55
  56 C...Pick up leg with highest virtuality.
  57   120 N=N+1
  58       JT=1
  59       IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
  60       KFLB=KFLS(JT)
  61       XB=XS(JT)
  62       DO 130 KFL=-6,6
  63   130 XFB(KFL)=XFS(JT,KFL)
  64       DSHR=2D0*SQRT(DSH)
  65       DSHZ=DSH/DBLE(ZS(JT))
  66       XE=MAX(XE0,XB*(1./(1.-PARP(66))-1.))
  67       IF(XB+XE.GE.0.999) THEN
  68         Q2B=0.
  69         GOTO 220
  70       ENDIF
  71
  72 C...Maximum Q2 without or with Q2 ordering. Effective Lambda and n_f.
  73       IF(MSTP(62).LE.1) THEN
  74         Q2B=0.5*(1./ZS(JT)+1.)*Q2S(JT)+0.5*(1./ZS(JT)-1.)*(Q2S(3-JT)-
  75      &  SNGL(DSH)+SQRT((SNGL(DSH)+Q2S(1)+Q2S(2))**2+8.*Q2S(1)*Q2S(2)*
  76      &  ZS(JT)/(1.-ZS(JT))))
  77         TEVB=LOG(PARP(63)*Q2B/ALAM(JT)**2)
  78       ELSE
  79         Q2B=Q2S(JT)
  80         TEVB=TEVS(JT)
  81       ENDIF
  82       ALSDUM=ULALPS_HIJING(PARP(63)*Q2B)
  83       TEVB=TEVB+2.*LOG(ALAM(JT)/PARU(117))
  84       TEVBSV=TEVB
  85       ALAM(JT)=PARU(117)
  86       B0=(33.-2.*MSTU(118))/6.
  87
  88 C...Calculate Altarelli-Parisi and structure function weights.
  89       DO 140 KFL=-6,6
  90       WTAP(KFL)=0.
  91   140 WTSF(KFL)=0.
  92       IF(KFLB.EQ.21) THEN
  93         WTAPQ=16.*(1.-SQRT(XB+XE))/(3.*SQRT(XB))
  94         DO 150 KFL=-MSTP(54),MSTP(54)
  95         IF(KFL.EQ.0) WTAP(KFL)=6.*LOG((1.-XB)/XE)
  96   150   IF(KFL.NE.0) WTAP(KFL)=WTAPQ
  97       ELSE
  98         WTAP(0)=0.5*XB*(1./(XB+XE)-1.)
  99         WTAP(KFLB)=8.*LOG((1.-XB)*(XB+XE)/XE)/3.
 100       ENDIF
 101   160 WTSUM=0.
 102       IF(KFLB.NE.21) XFBO=XFB(KFLB)
 103       IF(KFLB.EQ.21) XFBO=XFB(0)
 104 C***************************************************************
 105 C**********ERROR HAS OCCURED HERE
 106       IF(XFBO.EQ.0.0) THEN
 107                 WRITE(MSTU(11),1000)
 108                 WRITE(MSTU(11),1001) KFLB,XFB(KFLB)
 109                 XFBO=0.00001
 110       ENDIF
 111 C****************************************************************
 112       DO 170 KFL=-MSTP(54),MSTP(54)
 113       WTSF(KFL)=XFB(KFL)/XFBO
 114   170 WTSUM=WTSUM+WTAP(KFL)*WTSF(KFL)
 115       WTSUM=MAX(0.0001,WTSUM)
 116
 117 C...Choose new t: fix alpha_s, alpha_s(Q2), alpha_s(k_T2).
 118   180 IF(MSTP(64).LE.0) THEN
 119         TEVB=TEVB+LOG(RLU_HIJING(0))*PARU(2)/(PARU(111)*WTSUM)
 120       ELSEIF(MSTP(64).EQ.1) THEN
 121         TEVB=TEVB*EXP(MAX(-100.,LOG(RLU_HIJING(0))*B0/WTSUM))
 122       ELSE
 123         TEVB=TEVB*EXP(MAX(-100.,LOG(RLU_HIJING(0))*B0/(5.*WTSUM)))
 124       ENDIF
 125   190 Q2REF=ALAM(JT)**2*EXP(TEVB)
 126       Q2B=Q2REF/PARP(63)
 127
 128 C...Evolution ended or select flavour for branching parton.
 129       IF(Q2B.LT.PARP(62)**2) THEN
 130         Q2B=0.
 131       ELSE
 132         WTRAN=RLU_HIJING(0)*WTSUM
 133         KFLA=-MSTP(54)-1
 134   200   KFLA=KFLA+1
 135         WTRAN=WTRAN-WTAP(KFLA)*WTSF(KFLA)
 136         IF(KFLA.LT.MSTP(54).AND.WTRAN.GT.0.) GOTO 200
 137         IF(KFLA.EQ.0) KFLA=21
 138
 139 C...Choose z value and corrective weight.
 140         IF(KFLB.EQ.21.AND.KFLA.EQ.21) THEN
 141           Z=1./(1.+((1.-XB)/XB)*(XE/(1.-XB))**RLU_HIJING(0))
 142           WTZ=(1.-Z*(1.-Z))**2
 143         ELSEIF(KFLB.EQ.21) THEN
 144           Z=XB/(1.-RLU_HIJING(0)*(1.-SQRT(XB+XE)))**2
 145           WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
 146         ELSEIF(KFLA.EQ.21) THEN
 147           Z=XB*(1.+RLU_HIJING(0)*(1./(XB+XE)-1.))
 148           WTZ=1.-2.*Z*(1.-Z)
 149         ELSE
 150           Z=1.-(1.-XB)*(XE/((XB+XE)*(1.-XB)))**RLU_HIJING(0)
 151           WTZ=0.5*(1.+Z**2)
 152         ENDIF
 153
 154 C...Option with resummation of soft gluon emission as effective z shift.
 155         IF(MSTP(65).GE.1) THEN
 156           RSOFT=6.
 157           IF(KFLB.NE.21) RSOFT=8./3.
 158           Z=Z*(TEVB/TEVS(JT))**(RSOFT*XE/((XB+XE)*B0))
 159           IF(Z.LE.XB) GOTO 180
 160         ENDIF
 161
 162 C...Option with alpha_s(k_T2)Q2): demand k_T2 > cutoff, reweight.
 163         IF(MSTP(64).GE.2) THEN
 164           IF((1.-Z)*Q2B.LT.PARP(62)**2) GOTO 180
 165           ALPRAT=TEVB/(TEVB+LOG(1.-Z))
 166           IF(ALPRAT.LT.5.*RLU_HIJING(0)) GOTO 180
 167           IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
 168         ENDIF
 169
 170 C...Option with angular ordering requirement.
 171         IF(MSTP(62).GE.3) THEN
 172           THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
 173           IF(THE2T.GT.THE2(JT)) GOTO 180
 174         ENDIF
 175
 176 C...Weighting with new structure functions.
 177         CALL PYSTFU_HIJING(MINT(10+JT),XB,Q2REF,XFN,JT)
 178         IF(KFLB.NE.21) XFBN=XFN(KFLB)
 179         IF(KFLB.EQ.21) XFBN=XFN(0)
 180         IF(XFBN.LT.1E-20) THEN
 181           IF(KFLA.EQ.KFLB) THEN
 182             TEVB=TEVBSV
 183             WTAP(KFLB)=0.
 184             GOTO 160
 185           ELSEIF(TEVBSV-TEVB.GT.0.2) THEN
 186             TEVB=0.5*(TEVBSV+TEVB)
 187             GOTO 190
 188           ELSE
 189             XFBN=1E-10
 190           ENDIF
 191         ENDIF
 192         DO 210 KFL=-MSTP(54),MSTP(54)
 193   210   XFB(KFL)=XFN(KFL)
 194         XA=XB/Z
 195         CALL PYSTFU_HIJING(MINT(10+JT),XA,Q2REF,XFA,JT)
 196         IF(KFLA.NE.21) XFAN=XFA(KFLA)
 197         IF(KFLA.EQ.21) XFAN=XFA(0)
 198         IF(XFAN.LT.1E-20) GOTO 160
 199         IF(KFLA.NE.21) WTSFA=WTSF(KFLA)
 200         IF(KFLA.EQ.21) WTSFA=WTSF(0)
 201         IF(WTZ*XFAN/XFBN.LT.RLU_HIJING(0)*WTSFA) GOTO 160
 202       ENDIF
 203
 204 C...Define two hard scatterers in their CM-frame.
 205   220 IF(N.EQ.NS+2) THEN
 206         DQ2(JT)=Q2B
 207         DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
 208         DO 240 JR=1,2
 209         I=NS+JR
 210         IF(JR.EQ.1) IPO=IPUS1
 211         IF(JR.EQ.2) IPO=IPUS2
 212         DO 230 J=1,5
 213         K(I,J)=0
 214         P(I,J)=0.
 215   230   V(I,J)=0.
 216         K(I,1)=14
 217         K(I,2)=KFLS(JR+2)
 218         K(I,4)=IPO
 219         K(I,5)=IPO
 220         P(I,3)=DPLCM*(-1)**(JR+1)
 221         P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR
 222         P(I,5)=-SQRT(SNGL(DQ2(JR)))
 223         K(IPO,1)=14
 224         K(IPO,3)=I
 225         K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
 226   240   K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
 227
 228 C...Find maximum allowed mass of timelike parton.
 229       ELSEIF(N.GT.NS+2) THEN
 230         JR=3-JT
 231         DQ2(3)=Q2B
 232         DPC(1)=P(IS(1),4)
 233         DPC(2)=P(IS(2),4)
 234         DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
 235         DPD(1)=DSH+DQ2(JR)+DQ2(JT)
 236         DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
 237         DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
 238         DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
 239         IKIN=0
 240         IF(Q2S(JR).GE.(0.5*PARP(62))**2.AND.DPD(1)-DPD(3).GE.
 241      &  1D-10*DPD(1)) IKIN=1
 242         IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/DBLE(ZS(JT))-DQ2(3))*(DSH/
 243      &  (DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3)))
 244         IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/(2.*
 245      &  DQ2(JR))-DQ2(JT)-DQ2(3)
 246
 247 C...Generate timelike parton shower (if required).
 248         IT=N
 249         DO 250 J=1,5
 250         K(IT,J)=0
 251         P(IT,J)=0.
 252   250   V(IT,J)=0.
 253         K(IT,1)=3
 254         K(IT,2)=21
 255         IF(KFLB.EQ.21.AND.KFLS(JT+2).NE.21) K(IT,2)=-KFLS(JT+2)
 256         IF(KFLB.NE.21.AND.KFLS(JT+2).EQ.21) K(IT,2)=KFLB
 257         P(IT,5)=ULMASS_HIJING(K(IT,2))
 258         IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
 259         IF(MSTP(63).GE.1) THEN
 260           P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
 261           P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
 262           IF(MSTP(63).EQ.1) THEN
 263             Q2TIM=DMSMA
 264           ELSEIF(MSTP(63).EQ.2) THEN
 265             Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
 266           ELSE
 267 C'''Here remains to introduce angular ordering in first branching.
 268             Q2TIM=DMSMA
 269           ENDIF
 270           CALL LUSHOW_HIJING(IT,0,SQRT(Q2TIM))
 271           IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
 272         ENDIF
 273
 274 C...Reconstruct kinematics of branching: timelike parton shower.
 275         DMS=P(IT,5)**2
 276         IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
 277         IF(IKIN.EQ.1) DPT2=(DMSMA-DMS)*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
 278      &  DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/(4.*DSH*DPC(3)**2)
 279         IF(DPT2.LT.0.) GOTO 100
 280         DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
 281      &  DSHR)/DPC(3)-DPC(3)
 282         P(IT,1)=SQRT(SNGL(DPT2))
 283         P(IT,3)=DPB(1)*(-1)**(JT+1)
 284         P(IT,4)=(DSHZ-DSH-DMS)/DSHR
 285         IF(N.GE.IT+1) THEN
 286           DPB(1)=SQRT(DPB(1)**2+DPT2)
 287           DPB(2)=SQRT(DPB(1)**2+DMS)
 288           DPB(3)=P(IT+1,3)
 289           DPB(4)=SQRT(DPB(3)**2+DMS)
 290           DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
 291      &    DPB(1))
 292           CALL LUDBRB_HIJING(IT+1,N,0.,0.,0D0,0D0,DBEZ)
 293           THE=ULANGL_HIJING(P(IT,3),P(IT,1))
 294           CALL LUDBRB_HIJING(IT+1,N,THE,0.,0D0,0D0,0D0)
 295         ENDIF
 296
 297 C...Reconstruct kinematics of branching: spacelike parton.
 298         DO 260 J=1,5
 299         K(N+1,J)=0
 300         P(N+1,J)=0.
 301   260   V(N+1,J)=0.
 302         K(N+1,1)=14
 303         K(N+1,2)=KFLB
 304         P(N+1,1)=P(IT,1)
 305         P(N+1,3)=P(IT,3)+P(IS(JT),3)
 306         P(N+1,4)=P(IT,4)+P(IS(JT),4)
 307         P(N+1,5)=-SQRT(SNGL(DQ2(3)))
 308
 309 C...Define colour flow of branching.
 310         K(IS(JT),3)=N+1
 311         K(IT,3)=N+1
 312         ID1=IT
 313         IF((K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(ID1,2).GT.0.AND.
 314      &  K(ID1,2).NE.21).OR.(K(N+1,2).LT.0.AND.K(ID1,2).EQ.21).OR.
 315      &  (K(N+1,2).EQ.21.AND.K(ID1,2).EQ.21.AND.RLU_HIJING(0).GT.0.5).OR.
 316      &  (K(N+1,2).EQ.21.AND.K(ID1,2).LT.0)) ID1=IS(JT)
 317         ID2=IT+IS(JT)-ID1
 318         K(N+1,4)=K(N+1,4)+ID1
 319         K(N+1,5)=K(N+1,5)+ID2
 320         K(ID1,4)=K(ID1,4)+MSTU(5)*(N+1)
 321         K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
 322         K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
 323         K(ID2,5)=K(ID2,5)+MSTU(5)*(N+1)
 324         N=N+1
 325
 326 C...Boost to new CM-frame.
 327         CALL LUDBRB_HIJING(NS+1,N,0.,0.,-DBLE((P(N,1)+P(IS(JR),1))/(P(N
 328      $       ,4)+P(IS(JR),4))),0D0,-DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)
 329      $       +P(IS(JR),4))))
 330         IR=N+(JT-1)*(IS(1)-N)
 331         CALL LUDBRB_HIJING(NS+1,N,-ULANGL_HIJING(P(IR,3),P(IR,1)),PARU(2
 332      $       )*RLU_HIJING(0),0D0,0D0,0D0)
 333       ENDIF
 334
 335 C...Save quantities, loop back.
 336       IS(JT)=N
 337       Q2S(JT)=Q2B
 338       DQ2(JT)=Q2B
 339       IF(MSTP(62).GE.3) THE2(JT)=THE2T
 340       DSH=DSHZ
 341       IF(Q2B.GE.(0.5*PARP(62))**2) THEN
 342         KFLS(JT+2)=KFLS(JT)
 343         KFLS(JT)=KFLA
 344         XS(JT)=XA
 345         ZS(JT)=Z
 346         DO 270 KFL=-6,6
 347   270   XFS(JT,KFL)=XFA(KFL)
 348         TEVS(JT)=TEVB
 349       ELSE
 350         IF(JT.EQ.1) IPU1=N
 351         IF(JT.EQ.2) IPU2=N
 352       ENDIF
 353       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
 354          CALL LUERRM_HIJING(11
 355      $        ,'(PYSSPA_HIJING:) no more memory left in LUJETS_HIJING')
 356         IF(MSTU(21).GE.1) N=NS
 357         IF(MSTU(21).GE.1) RETURN
 358       ENDIF
 359       IF(MAX(Q2S(1),Q2S(2)).GE.(0.5*PARP(62))**2.OR.N.LE.NS+1) GOTO 120
 360
 361 C...Boost hard scattering partons to frame of shower initiators.
 362       DO 280 J=1,3
 363   280 ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
 364       DO 290 J=1,5
 365   290 P(N+2,J)=P(NS+1,J)
 366       ROBOT=ROBO(3)**2+ROBO(4)**2+ROBO(5)**2
 367       IF(ROBOT.GE.0.999999) THEN
 368         ROBOT=1.00001*SQRT(ROBOT)
 369         ROBO(3)=ROBO(3)/ROBOT
 370         ROBO(4)=ROBO(4)/ROBOT
 371         ROBO(5)=ROBO(5)/ROBOT
 372       ENDIF
 373       CALL LUDBRB_HIJING(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
 374      &-DBLE(ROBO(5)))
 375       ROBO(2)=ULANGL_HIJING(P(N+2,1),P(N+2,2))
 376       ROBO(1)=ULANGL_HIJING(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
 377       CALL LUDBRB_HIJING(MINT(83)+5,NS,ROBO(1),ROBO(2),DBLE(ROBO(3)),
 378      &DBLE(ROBO(4)),DBLE(ROBO(5)))
 379
 380 C...Store user information. Reset Lambda value.
 381       K(IPU1,3)=MINT(83)+3
 382       K(IPU2,3)=MINT(83)+4
 383       DO 300 JT=1,2
 384       MINT(12+JT)=KFLS(JT)
 385   300 VINT(140+JT)=XS(JT)
 386       PARU(111)=ALAMS
 387  1000 FORMAT(5X,'structure function has a zero point here')
 388  1001 FORMAT(5X,'xf(x,i=',I5,')=',F10.5)
 389
 390       RETURN
 391       END