PYTHIA6/QPYTHIA/q-pyshow.1.0.F

   1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   2 C
   3 C     AUXILIARY ROUTINES FOR Q-PYTHIA version 1.0.
   4 C
   5 C     DATE: 26.09.2008.
   6 C
   7 C     AUTHORS: N. Armesto, L. Cunqueiro and C. A. Salgado
   8 C              Departamento de Fisica de Particulas and IGFAE
   9 C              Universidade de Santiago de Compostela
  10 C              15706 Santiago de Compostela, Spain
  11 C
  12 C     EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es,
  13 C             Carlos.Salgado@cern.ch
  14 C
  15 C     CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18.
  16 C              NOT to be modified by user.
  17 C
  18 C     WHEN USING Q-PYTHIA, PLEASE QUOTE:
  19 C
  20 C     1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado,
  21 C        in preparation.
  22 C     2) T. Sjostrand, S. Mrenna and P. Skands,
  23 C        ``PYTHIA 6.4 physics and manual,''
  24 C        JHEP 0605 (2006) 026 [arXiv:hep-ph/0603175].
  25 C
  26 C     DISCLAIMER: this program comes without any guarantees. Beware of
  27 C                 errors and use common sense when interpreting results.
  28 C                 Any modifications are done under exclusive
  29 C                 maker's resposibility.
  30 C
  31 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  32 C
  33 c     VERSION WITH THE LARGE X BEHAVIOR OF THE MEDIUM PART INTRODUCED
  34 C     BY MULTIPLYING BY THE NUMERATOR OF THE COLLINEAR PART OF THE
  35 C     VACUUM SPLITTING FUNCTION
  36 c
  37       function splitq1(w)
  38 c     to integrate, adding vaccum plus medium q -> qg
  39       implicit double precision (a-h,o-z)
  40       z=w
  41       auxz=z*(1.d0-z)
  42       auxq=splitgq(z)+splitmedq1(z)
  43       if (auxq .gt. 0.d0) then
  44          splitq1=auxq
  45       else
  46          splitq1=0.d0
  47       endif
  48       return
  49       end
  50 c
  51       function splitg1(w)
  52 c     to integrate, adding vaccum plus medium g -> gg, and g -> qqbar
  53       implicit double precision (a-h,o-z)
  54       z=w
  55       auxz=z*(1.d0-z)
  56 c     argument of running coupling is taken as kt of emission
  57       auxg=splitgg(z)+splitmedg1(z)
  58       if (auxg .gt. 0.d0) then
  59          splitg1=(auxg+splitqqbar(z))
  60       else
  61          splitg1=splitqqbar(z)
  62       endif
  63       return
  64       end
  65 c
  66       function splitq2(w)
  67 c     to integrate, adding vaccum plus medium q -> qg
  68       implicit double precision (a-h,o-z)
  69       z=w
  70       auxq=splitgq(z)+splitmedq2(z)
  71       if (auxq .gt. 0.d0) then
  72          splitq2=auxq
  73       else
  74          splitq2=0.d0
  75       endif
  76       if(splitmedq2(z).eq.0.) then
  77       endif
  78       return
  79       end
  80 c
  81       function splitg2(z)
  82 c     to integrate, adding vaccum plus medium g -> gg, and g -> qqbar
  83       implicit double precision (a-h,o-z)
  84       auxg=splitgg(z)+splitmedg2(z)
  85       if(auxg.gt.0.d0) then
  86       splitg2=auxg
  87       else
  88       splitg2=0.d0
  89       endif
  90       if(splitmedg2(z).eq.0.) then
  91       endif
  92       return
  93       end
  94 c
  95       function splitgq(z)
  96 c     q -> qg splitting kernel at 1 loop for the vacuum
  97       implicit double precision (a-h,o-z)
  98       xnc=3.d0
  99       splitgq=(0.5d0*(xnc-1.d0/xnc))*(1.d0+z*z)/(1.d0-z)
 100       return
 101       end
 102 c
 103       function splitgg(z)
 104 c     g -> gg splitting kernel at 1 loop for the vacuum
 105       implicit double precision (a-h,o-z)
 106       xnc=3.d0
 107       auxz=z*(1.d0-z)
 108       auxz2=1.d0-auxz
 109       splitgg=xnc*auxz2*auxz2/auxz
 110       return
 111       end
 112 c
 113       function splitqqbar(z)
 114 c     g -> qqbar splitting kernel at 1 loop
 115       implicit double precision (a-h,o-z)
 116       xnf=5.d0
 117       auxz=1.d0-z
 118       splitqqbar=0.5d0*xnf*(z*z+auxz*auxz)
 119       return
 120       end
 121 c
 122       function splitmedg1(z)
 123 c     g -> gg splitting kernel at 1 loop for the medium
 124       implicit double precision (a-h,o-z)
 125       common/qpc1/eee,qhatl,omegac
 126       common/qpvir1/pmed
 127       xnc=3.d0
 128       pi=dacos(-1.d0)
 129       if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
 130          splitmedg1=0.d0
 131       else
 132 c     symmetrized by hand with respect to 1/2
 133          if (z .ge. 0.5d0) then
 134             zz=z
 135          else
 136             zz=1.d0-z
 137          endif
 138          t=pmed*pmed
 139          auxz=1.d0-zz
 140          auxz2=zz*auxz
 141          ome=eee*auxz/omegac
 142          xkappa2=auxz2*t/qhatl
 143          fff=genspec(ome,xkappa2)
 144 cc     1/2 to avoid double counting
 145 c         splitmedg=0.5d0*xnc*2.d0*pi*zz*t*fff/qhatl
 146 c     we multiply by max(z,1-z) to introduce the large z behavior from the
 147 c     numerator in the vacuum
 148          flx=max(zz,auxz)
 149 c     1/2 to avoid double counting
 150          splitmedg1=0.5*flx*xnc*2.d0*pi*zz*t*fff/qhatl
 151       endif
 152       return
 153       end
 154 c
 155       function splitmedq1(z)
 156 c     q -> qg splitting kernel at 1 loop for the medium
 157       implicit double precision (a-h,o-z)
 158       common/qpc1/eee,qhatl,omegac
 159       common/qpvir1/pmed
 160       xnc=3.d0
 161       pi=dacos(-1.d0)
 162       if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
 163          splitmedq1=0.d0
 164       else
 165          t=pmed*pmed
 166          auxz=1.d0-z
 167          auxz2=z*auxz
 168          ome=eee*auxz/omegac
 169          xkappa2=auxz2*t/qhatl
 170          fff=genspec(ome,xkappa2)
 171 c         splitmedq=(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
 172 c     we multiply by 1+z**2 to introduce the large z behavior from the
 173 c     numerator in the vacuum
 174          flx=0.5d0*(1.d0+z*z)
 175          splitmedq1=flx*(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
 176       endif
 177       return
 178       end
 179 c
 180       function splitmedg2(z)
 181 c     g -> gg splitting kernel at 1 loop for the medium
 182       implicit double precision (a-h,o-z)
 183       common/qpc1/eee,qhatl,omegac
 184       common/qpvir2/virt
 185       xnc=3.d0
 186       pi=dacos(-1.d0)
 187       if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
 188          splitmedg2=0.d0
 189       else
 190 c     symmetrized by hand with respect to 1/2
 191          if (z .ge. 0.5d0) then
 192             zz=z
 193          else
 194             zz=1.d0-z
 195          endif
 196          t=virt
 197          auxz=1.d0-zz
 198          auxz2=zz*auxz
 199          ome=eee*auxz/omegac
 200          xkappa2=auxz2*t/qhatl
 201          fff=genspec(ome,xkappa2)
 202 cc     1/2 to avoid double counting
 203 c         splitmedg=0.5d0*xnc*2.d0*pi*zz*t*fff/qhatl
 204 c     we multiply by max(z,1-z) to introduce the large z behavior from the
 205 c     numerator in the vacuum
 206          flx=max(zz,auxz)
 207 c     1/2 to avoid double counting
 208          splitmedg2=0.5*flx*xnc*2.d0*pi*zz*t*fff/qhatl
 209       endif
 210       return
 211       end
 212 c
 213       function splitmedq2(z)
 214 c     q -> qg splitting kernel at 1 loop for the medium
 215       implicit double precision (a-h,o-z)
 216       common/qpc1/eee,qhatl,omegac
 217       common/qpvir2/virt
 218       xnc=3.d0
 219       pi=dacos(-1.d0)
 220       if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
 221          splitmedq2=0.d0
 222       else
 223          t=virt
 224          auxz=1.d0-z
 225          auxz2=z*auxz
 226          ome=eee*auxz/omegac
 227          xkappa2=auxz2*t/qhatl
 228          fff=genspec(ome,xkappa2)
 229 c         splitmedq=(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
 230 c     we multiply by 1+z**2 to introduce the large z behavior from the
 231 c     numerator in the vacuum
 232          flx=0.5d0*(1.d0+z*z)
 233          splitmedq2=flx*(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
 234       endif
 235       return
 236       end
 237 c
 238       function genspec(ome,xk2)
 239 C     THIS FUNCTION GENERATES (omega/omegac) dI/d(omega/omegac) dkappa2,
 240 C     omegac=qhat L/2, kappa2=kt2/(qhat L), in the mss approximation for m=0,
 241 c     using interpolation and extrapolation. It reads file grid-qp.dat.
 242 c     ome=omega/omegac, xk2=kappa2.
 243 C     MAXIMUM GRID 101 TIMES 101, MODIFY ARRAY DIMENSIONS IF EXCEEDED.
 244 c     alphas=1, cr=1.
 245       implicit double precision (a-h,o-z)
 246       dimension xkap2(101), xlkap2(101), xome(101), xlome(101)
 247       dimension xspec(101,101)
 248       dimension aux1(101), aux2(101)
 249       save xkap2, xlkap2, xome, xlome, xspec, npkap, npome
 250       DATA IFLAG/0/
 251 c     WE READ THE GRID ONLY THE FIRST TIME.
 252       IF (IFLAG .EQ. 0) THEN
 253 c         print*, 'reading qgrid'
 254          open(11,file=
 255      +'/ed22/dfs/work/ALICE/offline/AliRoot/trunk/PYTHIA6/qgrid',
 256      +status='old')
 257          read(11,*) npkap
 258          read(11,*) npome
 259          npkap=npkap+1
 260          npome=npome+1
 261          do 10 i=1, npkap, 1
 262             read(11,*) xkap2(i), xlkap2(i)
 263 10       continue
 264          do 20 i=1, npome, 1
 265             read(11,*) xome(i), xlome(i)
 266 20       continue
 267          do 30 j=1, npome, 1
 268             do 40 i=1, npkap, 1
 269                read(11,*) xspec(i,j)
 270 40          continue
 271 30       continue
 272          close(11)
 273          iflag=1
 274       ENDIF
 275 c     cases
 276 c     for ome>largest value set to 0,
 277 c     for xk2< smallest value frozen,
 278 c     for xk2> largest value 1/kappa4 extrapolation.
 279       if (ome .gt. xome(npome)) then
 280          genspec=0.d0
 281       elseif (ome .lt. xome(1)) then
 282          scal=.05648d0*dexp(1.674d0*ome)*dlog(.136d0/ome)/(ome**.5397d0)
 283          scal=0.25d0*9.d0*scal/xspec(1,1)
 284          if (xk2 .le. xkap2(1)) then
 285             genspec=scal*xspec(1,1)
 286          elseif (xk2 .eq. xkap2(npkap)) then
 287             genspec=scal*xspec(npkap,1)
 288          elseif (xk2 .gt. xkap2(npkap)) then
 289             genspec=scal*xspec(npkap,1)*
 290      >              xkap2(npkap)*xkap2(npkap)/(xk2*xk2)
 291          else
 292             do 50 i=1, npkap, 1
 293                aux1(i)=xspec(i,1)
 294 50          continue
 295             genspec=scal*ddivdif(aux1,xlkap2,npkap,dlog(xk2),4)
 296          endif
 297       else
 298          iexact=-1
 299          if (ome .eq. xome(1)) then
 300             iexact=1
 301             goto 70
 302          else
 303             do 60 i=1, npome-1, 1
 304                if (ome .eq. xome(i+1)) then
 305                   iexact=i+1
 306                   goto 70
 307                elseif (ome .lt. xome(i+1)) then
 308                   iprev=i
 309                   ipost=i+1
 310                   goto 70
 311                endif
 312 60          continue
 313 70          continue
 314          endif
 315          if (iexact .gt. 0) then
 316             if (xk2 .le. xkap2(1)) then
 317                genspec=xspec(1,iexact)
 318             elseif (xk2 .eq. xkap2(npkap)) then
 319                genspec=xspec(npkap,iexact)
 320             elseif (xk2 .gt. xkap2(npkap)) then
 321                genspec=xspec(npkap,iexact)*
 322      >                 xkap2(npkap)*xkap2(npkap)/(xk2*xk2)
 323             else
 324                do 80 i=1, npkap, 1
 325                   aux1(i)=xspec(i,iexact)
 326 80             continue
 327                genspec=ddivdif(aux1,xlkap2,npkap,dlog(xk2),4)
 328             endif
 329          else
 330             if (xk2 .le. xkap2(1)) then
 331                genprev=xspec(1,iprev)
 332                genpost=xspec(1,ipost)
 333             elseif (xk2 .eq. xkap2(npkap)) then
 334                genprev=xspec(npkap,iprev)
 335                genpost=xspec(npkap,ipost)
 336             elseif (xk2 .gt. xkap2(npkap)) then
 337                genprev=xspec(npkap,iprev)*
 338      >                 xkap2(npkap)*xkap2(npkap)/(xk2*xk2)
 339                genpost=xspec(npkap,ipost)*
 340      >                 xkap2(npkap)*xkap2(npkap)/(xk2*xk2)
 341             else
 342                do 90 i=1, npkap, 1
 343                   aux1(i)=xspec(i,iprev)
 344                   aux2(i)=xspec(i,ipost)
 345 90             continue
 346                genprev=ddivdif(aux1,xlkap2,npkap,dlog(xk2),4)
 347                genpost=ddivdif(aux2,xlkap2,npkap,dlog(xk2),4)
 348             endif
 349             g12=genprev-genpost
 350             xl12=xlome(iprev)-xlome(ipost)
 351             c1=g12/xl12
 352             c2=genprev-c1*xlome(iprev)
 353             genspec=c1*dlog(ome)+c2
 354          endif
 355       endif
 356 c
 357       RETURN
 358       END
 359 C
 360 *
 361 * $Id: divdif.F,v 1.1.1.1 1996/02/15 17:48:36 mclareni Exp $
 362 *
 363 * $Log: divdif.F,v $
 364 * Revision 1.1.1.1  1996/02/15 17:48:36  mclareni
 365 * Kernlib
 366 *
 367 *
 368       FUNCTION DDIVDIF(F,A,NN,X,MM)
 369 c     copy of cernlib divdif in double precision.
 370       implicit double precision (a-h,o-z)
 371       DIMENSION A(NN),F(NN),T(20),D(20)
 372       LOGICAL EXTRA
 373       LOGICAL MFLAG,RFLAG
 374       DATA MMAX/10/
 375 C
 376 C  TABULAR INTERPOLATION USING SYMMETRICALLY PLACED ARGUMENT POINTS.
 377 C
 378 C  START.  FIND SUBSCRIPT IX OF X IN ARRAY A.
 379       IF( (NN.LT.2) .OR. (MM.LT.1) ) GO TO 601
 380       N=NN
 381       M=MIN0(MM,MMAX,N-1)
 382       MPLUS=M+1
 383       IX=0
 384       IY=N+1
 385       IF(A(1).GT.A(N)) GO TO 4
 386 C     (SEARCH INCREASING ARGUMENTS.)
 387     1    MID=(IX+IY)/2
 388          IF(X.GE.A(MID)) GO TO 2
 389             IY=MID
 390             GO TO 3
 391 C        (IF TRUE.)
 392     2       IX=MID
 393     3    IF(IY-IX.GT.1) GO TO 1
 394          GO TO 7
 395 C     (SEARCH DECREASING ARGUMENTS.)
 396     4    MID=(IX+IY)/2
 397          IF(X.LE.A(MID)) GO TO 5
 398             IY=MID
 399             GO TO 6
 400 C        (IF TRUE.)
 401     5       IX=MID
 402     6    IF(IY-IX.GT.1) GO TO 4
 403 C
 404 C  COPY REORDERED INTERPOLATION POINTS INTO (T(I),D(I)), SETTING
 405 C  *EXTRA* TO TRUE IF M+2 POINTS TO BE USED.
 406     7 NPTS=M+2-MOD(M,2)
 407       IP=0
 408       L=0
 409       GO TO 9
 410     8    L=-L
 411          IF(L.GE.0) L=L+1
 412     9    ISUB=IX+L
 413          IF((1.LE.ISUB).AND.(ISUB.LE.N)) GO TO 501
 414 C        (SKIP POINT.)
 415             NPTS=MPLUS
 416             GO TO 11
 417 C        (INSERT POINT.)
 418  501        IP=IP+1
 419             T(IP)=A(ISUB)
 420             D(IP)=F(ISUB)
 421    11    IF(IP.LT.NPTS) GO TO 8
 422       EXTRA=NPTS.NE.MPLUS
 423 C
 424 C  REPLACE D BY THE LEADING DIAGONAL OF A DIVIDED-DIFFERENCE TABLE, SUP-
 425 C  PLEMENTED BY AN EXTRA LINE IF *EXTRA* IS TRUE.
 426       DO 14 L=1,M
 427          IF(.NOT.EXTRA) GO TO 12
 428             ISUB=MPLUS-L
 429             D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB))
 430    12    I=MPLUS
 431          DO 13 J=L,M
 432             ISUB=I-L
 433             D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB))
 434             I=I-1
 435    13    CONTINUE
 436    14 CONTINUE
 437 C
 438 C  EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES
 439 C  OF LAST DIFFERENCE IF *EXTRA* IS TRUE.
 440       SUM=D(MPLUS)
 441       IF(EXTRA) SUM=0.5*(SUM+D(M+2))
 442       J=M
 443       DO 15 L=1,M
 444          SUM=D(J)+(X-T(J))*SUM
 445          J=J-1
 446    15 CONTINUE
 447       DDIVDIF=SUM
 448       RETURN
 449 C
 450  601  CALL KERMTR('E105.1',LGFILE,MFLAG,RFLAG)
 451       DDIVDIF=0
 452       IF(MFLAG) THEN
 453          IF(LGFILE.EQ.0) THEN
 454             IF(MM.LT.1) WRITE(*,101) MM
 455             IF(NN.LT.2) WRITE(*,102) NN
 456          ELSE
 457             IF(MM.LT.1) WRITE(LGFILE,101) MM
 458             IF(NN.LT.2) WRITE(LGFILE,102) NN
 459          ENDIF
 460       ENDIF
 461       IF(.NOT.RFLAG) CALL ABEND
 462       RETURN
 463   101 FORMAT( 7X, 'FUNCTION DDIVDIF ... M =',I6,' IS LESS THAN 1')
 464   102 FORMAT( 7X, 'FUNCTION DDIVDIF ... N =',I6,' IS LESS THAN 2')
 465       END
 466 c
 467 C     COPY OF CERN DGAUSS
 468 C
 469       FUNCTION DGAUSS1(F,A,B,EPS)
 470       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 471       DIMENSION W(12),X(12)
 472       PARAMETER (Z1 = 1.D0, HF = Z1/2.D0, CST = 5.D0*Z1/1000.D0)
 473       DATA X
 474      1        /0.96028 98564 97536 23168 35608 68569 47D0,
 475      2         0.79666 64774 13626 73959 15539 36475 83D0,
 476      3         0.52553 24099 16328 98581 77390 49189 25D0,
 477      4         0.18343 46424 95649 80493 94761 42360 18D0,
 478      5         0.98940 09349 91649 93259 61541 73450 33D0,
 479      6         0.94457 50230 73232 57607 79884 15534 61D0,
 480      7         0.86563 12023 87831 74388 04678 97712 39D0,
 481      8         0.75540 44083 55003 03389 51011 94847 44D0,
 482      9         0.61787 62444 02643 74844 66717 64048 79D0,
 483      A         0.45801 67776 57227 38634 24194 42983 58D0,
 484      B         0.28160 35507 79258 91323 04605 01460 50D0,
 485      C         0.95012 50983 76374 40185 31933 54249 58D-1/
 486
 487       DATA W
 488      1        /0.10122 85362 90376 25915 25313 54309 96D0,
 489      2         0.22238 10344 53374 47054 43559 94426 24D0,
 490      3         0.31370 66458 77887 28733 79622 01986 60D0,
 491      4         0.36268 37833 78361 98296 51504 49277 20D0,
 492      5         0.27152 45941 17540 94851 78057 24560 18D-1,
 493      6         0.62253 52393 86478 92862 84383 69943 78D-1,
 494      7         0.95158 51168 24927 84809 92510 76022 46D-1,
 495      8         0.12462 89712 55533 87205 24762 82192 02D0,
 496      9         0.14959 59888 16576 73208 15017 30547 48D0,
 497      A         0.16915 65193 95002 53818 93120 79030 36D0,
 498      B         0.18260 34150 44923 58886 67636 67969 22D0,
 499      C         0.18945 06104 55068 49628 53967 23208 28D0/
 500       EXTERNAL F
 501       H=0.D0
 502       IF(B .EQ. A) GO TO 99
 503       CONST=CST/ABS(B-A)
 504       BB=A
 505     1 AA=BB
 506       BB=B
 507     2 C1=HF*(BB+AA)
 508       C2=HF*(BB-AA)
 509       S8=0.D0
 510       DO 3 I = 1,4
 511       U=C2*X(I)
 512     3 S8=S8+W(I)*(F(C1+U)+F(C1-U))
 513       S16=0.D0
 514       DO 4 I = 5,12
 515       U=C2*X(I)
 516     4 S16=S16+W(I)*(F(C1+U)+F(C1-U))
 517       S16=C2*S16
 518       IF(ABS(S16-C2*S8) .LE. EPS*(1.D0+ABS(S16))) THEN
 519        H=H+S16
 520        IF(BB .NE. B) GO TO 1
 521       ELSE
 522        BB=C1
 523        IF(1.D0+CONST*ABS(C2) .NE. 1.D0) GO TO 2
 524        H=0.D0
 525        WRITE(6,*) 'DGAUSS1: TOO HIGH ACCURACY REQUIRED'
 526        GO TO 99
 527       END IF
 528    99 DGAUSS1=H
 529       RETURN
 530       END
 531 c
 532       FUNCTION SIMDIS(Numb,zmin,nzur,RI)
 533 C     IT SIMULATES A RANDOM NUMBER ACCORDING TO A DISCRETE DISTRIBUTION GIVEN
 534 C     BY ARRAY YA AT POINTS XA. THOUGHT FOR PYTHIA (PYR(0)).
 535 C     N: NUMBER OF POINTS IN THE ARRAYS.
 536 C     XA: ARRAY OF X-VALUES.
 537 C     YA: ARRAY OF Y-VALUES.
 538 c     RI: VALUE OF THE INTEGRAL.
 539       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 540       DIMENSION XA(500), YA(500)
 541       common/qpc1/eee,qhatl,omegac
 542       dlz=(1.d0-2.d0*zmin)/500.d0
 543       do 1000 no=1,500
 544       xa(no)=zmin+no*dlz
 545       if(nzur.eq.1) ya(no)=splitq2(xa(no))
 546       if(nzur.eq.21) ya(no)=splitg2(xa(no))
 547       if(nzur.eq.3) ya(no)=splitqqbar(xa(no))
 548  1000 continue
 549       RAL=PYR(0)*RI
 550       XAUX=0.D0
 551       xauxold=0.d0
 552       DO 1001 I=2, Numb, 1
 553       XAUX=XAUX+(XA(I)-XA(I-1))*0.5D0*
 554      + (YA(I)+YA(I-1))
 555
 556            IF (XAUX .GE. RAL) GOTO 2011
 557
 558          IF (I .EQ. Numb) THEN
 559             SIMDIS=XA(I)
 560
 561             RETURN
 562          ENDIF
 563          XAUXOLD=XAUX
 564  1001  CONTINUE
 565  2011  SIMDIS=(XA(I)-XA(I-1))*(RAL-XAUXOLD)/(XAUX-XAUXOLD)+
 566      + XA(I-1)
 567
 568       RETURN
 569       END
 570
 571
 572       SUBROUTINE QPYROBO(XI,YI,ZI,TI,THE,PHI,BEX,BEY,BEZ,XP,YP,ZP,TP)
 573 C     N. Armesto, 16.04.2009
 574 C     performs a boost and rotation of (t,x,y,z) to (tp,xp,yp,zp):
 575 C     cut version of PYROBO, angles and boost parameters identical.
 576       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 577 C...Local arrays.
 578       DIMENSION ROT(3,3),VR(3),DV(4)
 579 C
 580       X=XI
 581       Y=YI
 582       Z=ZI
 583       T=TI
 584 C...Rotate, typically from z axis to direction (theta,phi).
 585       IF(THE**2+PHI**2.GT.1D-20) THEN
 586         ROT(1,1)=COS(THE)*COS(PHI)
 587         ROT(1,2)=-SIN(PHI)
 588         ROT(1,3)=SIN(THE)*COS(PHI)
 589         ROT(2,1)=COS(THE)*SIN(PHI)
 590         ROT(2,2)=COS(PHI)
 591         ROT(2,3)=SIN(THE)*SIN(PHI)
 592         ROT(3,1)=-SIN(THE)
 593         ROT(3,2)=0D0
 594         ROT(3,3)=COS(THE)
 595 C   Instead of loop 120 in PYROBO.
 596         VR(1)=X
 597         VR(2)=Y
 598         VR(3)=Z
 599 C   Instead of loop 130 in PYROBO.
 600         J=1
 601         X=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
 602         J=2
 603         Y=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
 604         J=3
 605         Z=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
 606       ENDIF
 607 C   If nothing happens...
 608       XP=X
 609       YP=Y
 610       ZP=Z
 611       TP=T
 612 C...Boost, typically from rest to momentum/energy=beta.
 613       IF(BEX**2+BEY**2+BEZ**2.GT.1D-20) THEN
 614         DBX=BEX
 615         DBY=BEY
 616         DBZ=BEZ
 617         DB=SQRT(DBX**2+DBY**2+DBZ**2)
 618         EPS1=1D0-1D-12
 619         IF(DB.GT.EPS1) THEN
 620 C...Rescale boost vector if too close to unity.
 621           CALL PYERRM(3,'(PYROBO:) boost vector too large')
 622           DBX=DBX*(EPS1/DB)
 623           DBY=DBY*(EPS1/DB)
 624           DBZ=DBZ*(EPS1/DB)
 625           DB=EPS1
 626         ENDIF
 627         DGA=1D0/SQRT(1D0-DB**2)
 628 C    Instead of loop 150 in PYROBO.
 629         DV(1)=X
 630         DV(2)=Y
 631         DV(3)=Z
 632         DV(4)=T
 633         DBV=DBX*DV(1)+DBY*DV(2)+DBZ*DV(3)
 634         DGABV=DGA*(DGA*DBV/(1D0+DGA)+DV(4))
 635         XP=DV(1)+DGABV*DBX
 636         YP=DV(2)+DGABV*DBY
 637         ZP=DV(3)+DGABV*DBZ
 638         TP=DGA*(DV(4)+DBV)
 639       ENDIF
 640 C
 641       RETURN
 642       END
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 653 C
 654 C     PYSHOW ROUTINE FOR Q-PYTHIA version 1.0.
 655 C
 656 C     DATE: 26.09.2008.
 657 C
 658 C     AUTHORS: N. Armesto, L. Cunqueiro and C. A. Salgado
 659 C              Departamento de Fisica de Particulas and IGFAE
 660 C              Universidade de Santiago de Compostela
 661 C              15706 Santiago de Compostela, Spain
 662 C
 663 C     EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es,
 664 C             Carlos.Salgado@cern.ch
 665 C
 666 C     CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18.
 667 C
 668 C     WHEN USING Q-PYTHIA, PLEASE QUOTE:
 669 C
 670 C     1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado,
 671 C        in preparation.
 672 C     2) T. Sjostrand, S. Mrenna and P. Skands,
 673 C        ``PYTHIA 6.4 physics and manual,''
 674 C        JHEP 0605 (2006) 026 [arXiv:hep-ph/0603175].
 675 C
 676 C     INSTRUCTIONS: initial parton position is initialized by a call
 677 C                   to user-defined routine qpygin(x0,y0,z0,t0),
 678 C                   where these are the initial coordinates in the
 679 C                   center-of-mass frame of the hard collision
 680 C                   (if applicable for the type of process you study).
 681 C                   The values of qhatL and omegac have to be computed
 682 C                   by the user, using his preferred medium model, in
 683 C                   routine qpygeo, which takes as input the position
 684 C                   x,y,z,t of the parton to branch, the trajectory
 685 C                   defined by the three-vector betax,betay,betaz,
 686 C                   (all values in the center-of-mass frame of the
 687 C                   hard collision), and  returns the value of qhatL
 688 C                   (in GeV**2) and omegac (in GeV).
 689 C                   Both routines are to be found at the end of this file.
 690 C
 691 C     DISCLAIMER: this program comes without any guarantees. Beware of
 692 C                 errors and use common sense when interpreting results.
 693 C                 Any modifications are done under exclusive
 694 C                 maker's resposibility.
 695 C
 696 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 697 C*********************************************************************
 698
 699 C...PYSHOW
 700 C...Generates timelike parton showers from given partons.
 701
 702       SUBROUTINE PYSHOWQ(IP1,IP2,QMAX)
 703
 704 C...Double precision and integer declarations.
 705       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 706       IMPLICIT INTEGER(I-N)
 707       INTEGER PYK,PYCHGE,PYCOMP
 708 C...Parameter statement to help give large particle numbers.
 709       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
 710      &KEXCIT=4000000,KDIMEN=5000000)
 711       PARAMETER (MAXNUR=500)
 712 Cacs+
 713       PARAMETER (NNPOS=4000)
 714       DIMENSION PPOS(NNPOS,4)
 715 Cacs-
 716 C...Commonblocks.
 717       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
 718       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
 719       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
 720       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
 721       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
 722       COMMON/PYINT1/MINT(400),VINT(400)
 723       SAVE /PYPART/,/PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/
 724 Cacs+
 725       common/qpc1/eee,qhatl,omegac
 726       common/qpvir1/pmed
 727       common/qpvir2/virt
 728       COMMON/QPLT/QPLTA1,QPLTA2,QPLTBX,QPLTBY,QPLTBZ
 729       external splitg1
 730       external splitq1
 731       external splitg2
 732       external splitq2
 733       external splitqqbar
 734       data iflag/0/
 735 Cacs-
 736 C...Local arrays.
 737       DIMENSION PMTH(5,140),PS(5),PMA(100),PMSD(100),IEP(100),IPA(100),
 738      &KFLA(100),KFLD(100),KFL(100),ITRY(100),ISI(100),ISL(100),DP(100),
 739      &DPT(5,4),KSH(0:140),KCII(2),NIIS(2),IIIS(2,2),THEIIS(2,2),
 740      &PHIIIS(2,2),ISII(2),ISSET(2),ISCOL(0:140),ISCHG(0:140),
 741      &IREF(1000)
 742 Cacs+
 743       IF (IFLAG .EQ. 0) THEN
 744          WRITE(MSTU(11),*)
 745          WRITE(MSTU(11),*) '*******************************************'
 746          WRITE(MSTU(11),*)
 747          WRITE(MSTU(11),*) '            Q-PYTHIA version 1.0'
 748          WRITE(MSTU(11),*)
 749          WRITE(MSTU(11),*) 'DATE: 26.09.2008'
 750          WRITE(MSTU(11),*)
 751          WRITE(MSTU(11),*) 'AUTHORS: N. Armesto, L. Cunqueiro and'
 752          WRITE(MSTU(11),*) '         C. A. Salgado'
 753          WRITE(MSTU(11),*) ' Departamento de Fisica de Particulas'
 754          WRITE(MSTU(11),*) ' and IGFAE'
 755          WRITE(MSTU(11),*) ' Universidade de Santiago de Compostela'
 756          WRITE(MSTU(11),*) ' 15706 Santiago de Compostela, Spain'
 757          WRITE(MSTU(11),*)
 758          WRITE(MSTU(11),*) 'EMAILS: nestor@fpaxp1.usc.es,'
 759          WRITE(MSTU(11),*) '        leticia@fpaxp1.usc.es,'
 760          WRITE(MSTU(11),*) '        Carlos.Salgado@cern.ch'
 761          WRITE(MSTU(11),*)
 762          WRITE(MSTU(11),*) 'NOTE: fixed to PYTHIA-6.4.18'
 763          WRITE(MSTU(11),*)
 764          WRITE(MSTU(11),*) 'WHEN USING Q-PYTHIA, PLEASE QUOTE:'
 765          WRITE(MSTU(11),*) '1) N. Armesto, G. Corcella, L. Cunqueiro'
 766          WRITE(MSTU(11),*) '   and C. A. Salgado, in preparation.'
 767          WRITE(MSTU(11),*) '2) T. Sjostrand, S. Mrenna and P. Skands,'
 768          WRITE(MSTU(11),*) '   PYTHIA 6.4 physics and manual,'
 769          WRITE(MSTU(11),*) '   JHEP 0605 (2006) 026'
 770          WRITE(MSTU(11),*) '   [arXiv:hep-ph/0603175].'
 771          WRITE(MSTU(11),*)
 772          WRITE(MSTU(11),*) 'INSTRUCTIONS: look at the web page and'
 773          WRITE(MSTU(11),*) ' header of modfied routine PYSHOW at the'
 774          WRITE(MSTU(11),*) ' end of Q-PYTHIA file.'
 775          WRITE(MSTU(11),*)
 776          WRITE(MSTU(11),*) 'DISCLAIMER: this program comes without any'
 777          WRITE(MSTU(11),*) ' guarantees. Beware of errors and use'
 778          WRITE(MSTU(11),*) ' common sense when interpreting results.'
 779          WRITE(MSTU(11),*) ' Any modifications are done under exclusive'
 780          WRITE(MSTU(11),*) ' makers resposibility.'
 781          WRITE(MSTU(11),*)
 782          WRITE(MSTU(11),*) '*******************************************'
 783          WRITE(MSTU(11),*)
 784          IFLAG=1
 785       ENDIF
 786 Cacs-
 787
 788 C...Check that QMAX not too low.
 789       IF(MSTJ(41).LE.0) THEN
 790         RETURN
 791       ELSEIF(MSTJ(41).EQ.1.OR.MSTJ(41).EQ.11) THEN
 792         IF(QMAX.LE.PARJ(82).AND.IP2.GE.-80) RETURN
 793       ELSE
 794         IF(QMAX.LE.MIN(PARJ(82),PARJ(83),PARJ(90)).AND.IP2.GE.-80)
 795      &  RETURN
 796       ENDIF
 797
 798 C...Store positions of shower initiating partons.
 799       MPSPD=0
 800       IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
 801         NPA=1
 802         IPA(1)=IP1
 803       ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
 804      &  MSTU(32))) THEN
 805         NPA=2
 806         IPA(1)=IP1
 807         IPA(2)=IP2
 808       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0
 809      &  .AND.IP2.GE.-80) THEN
 810         NPA=IABS(IP2)
 811         DO 100 I=1,NPA
 812           IPA(I)=IP1+I-1
 813   100   CONTINUE
 814       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.
 815      &IP2.EQ.-100) THEN
 816         MPSPD=1
 817         NPA=2
 818         IPA(1)=IP1+6
 819         IPA(2)=IP1+7
 820       ELSE
 821         CALL PYERRM(12,
 822      &  '(PYSHOW:) failed to reconstruct showering system')
 823         IF(MSTU(21).GE.1) RETURN
 824       ENDIF
 825
 826 C...Send off to PYPTFS for pT-ordered evolution if requested,
 827 C...if at least 2 partons, and without predefined shower branchings.
 828       IF((MSTJ(41).EQ.11.OR.MSTJ(41).EQ.12).AND.NPA.GE.2.AND.
 829      &MPSPD.EQ.0) THEN
 830         NPART=NPA
 831         DO 110 II=1,NPART
 832           IPART(II)=IPA(II)
 833           PTPART(II)=0.5D0*QMAX
 834   110   CONTINUE
 835         CALL PYPTFS(2,0.5D0*QMAX,0D0,PTGEN)
 836         RETURN
 837       ENDIF
 838
 839 C...Initialization of cutoff masses etc.
 840       DO 120 IFL=0,40
 841         ISCOL(IFL)=0
 842         ISCHG(IFL)=0
 843         KSH(IFL)=0
 844   120 CONTINUE
 845       ISCOL(21)=1
 846       KSH(21)=1
 847       PMTH(1,21)=PYMASS(21)
 848       PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25D0*PARJ(82)**2)
 849       PMTH(3,21)=2D0*PMTH(2,21)
 850       PMTH(4,21)=PMTH(3,21)
 851       PMTH(5,21)=PMTH(3,21)
 852       PMTH(1,22)=PYMASS(22)
 853       PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25D0*PARJ(83)**2)
 854       PMTH(3,22)=2D0*PMTH(2,22)
 855       PMTH(4,22)=PMTH(3,22)
 856       PMTH(5,22)=PMTH(3,22)
 857       PMQTH1=PARJ(82)
 858       IF(MSTJ(41).GE.2) PMQTH1=MIN(PARJ(82),PARJ(83))
 859       PMQT1E=MIN(PMQTH1,PARJ(90))
 860       PMQTH2=PMTH(2,21)
 861       IF(MSTJ(41).GE.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
 862       PMQT2E=MIN(PMQTH2,0.5D0*PARJ(90))
 863       DO 130 IFL=1,5
 864         ISCOL(IFL)=1
 865         IF(MSTJ(41).GE.2) ISCHG(IFL)=1
 866         KSH(IFL)=1
 867         PMTH(1,IFL)=PYMASS(IFL)
 868         PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PMQTH1**2)
 869         PMTH(3,IFL)=PMTH(2,IFL)+PMQTH2
 870         PMTH(4,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(82)**2)+PMTH(2,21)
 871         PMTH(5,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(83)**2)+PMTH(2,22)
 872   130 CONTINUE
 873       DO 140 IFL=11,15,2
 874         IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) ISCHG(IFL)=1
 875         IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) KSH(IFL)=1
 876         PMTH(1,IFL)=PYMASS(IFL)
 877         PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(90)**2)
 878         PMTH(3,IFL)=PMTH(2,IFL)+0.5D0*PARJ(90)
 879         PMTH(4,IFL)=PMTH(3,IFL)
 880         PMTH(5,IFL)=PMTH(3,IFL)
 881   140 CONTINUE
 882       PT2MIN=MAX(0.5D0*PARJ(82),1.1D0*PARJ(81))**2
 883       ALAMS=PARJ(81)**2
 884       ALFM=LOG(PT2MIN/ALAMS)
 885
 886 C...Check on phase space available for emission.
 887       IREJ=0
 888       DO 150 J=1,5
 889         PS(J)=0D0
 890   150 CONTINUE
 891       PM=0D0
 892       KFLA(2)=0
 893       DO 170 I=1,NPA
 894         KFLA(I)=IABS(K(IPA(I),2))
 895         PMA(I)=P(IPA(I),5)
 896 C...Special cutoff masses for initial partons (may be a heavy quark,
 897 C...squark, ..., and need not be on the mass shell).
 898         IR=30+I
 899         IF(NPA.LE.1) IREF(I)=IR
 900         IF(NPA.GE.2) IREF(I+1)=IR
 901         ISCOL(IR)=0
 902         ISCHG(IR)=0
 903         KSH(IR)=0
 904         IF(KFLA(I).LE.8) THEN
 905           ISCOL(IR)=1
 906           IF(MSTJ(41).GE.2) ISCHG(IR)=1
 907         ELSEIF(KFLA(I).EQ.11.OR.KFLA(I).EQ.13.OR.KFLA(I).EQ.15.OR.
 908      &  KFLA(I).EQ.17) THEN
 909           IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) ISCHG(IR)=1
 910         ELSEIF(KFLA(I).EQ.21) THEN
 911           ISCOL(IR)=1
 912         ELSEIF((KFLA(I).GE.KSUSY1+1.AND.KFLA(I).LE.KSUSY1+8).OR.
 913      &  (KFLA(I).GE.KSUSY2+1.AND.KFLA(I).LE.KSUSY2+8)) THEN
 914           ISCOL(IR)=1
 915         ELSEIF(KFLA(I).EQ.KSUSY1+21) THEN
 916           ISCOL(IR)=1
 917 C...QUARKONIA+++
 918 C...same for QQ~[3S18]
 919         ELSEIF(MSTP(148).GE.1.AND.(KFLA(I).EQ.9900443.OR.
 920      &  KFLA(I).EQ.9900553)) THEN
 921           ISCOL(IR)=1
 922 C...QUARKONIA---
 923         ENDIF
 924         IF(ISCOL(IR).EQ.1.OR.ISCHG(IR).EQ.1) KSH(IR)=1
 925         PMTH(1,IR)=PMA(I)
 926         IF(ISCOL(IR).EQ.1.AND.ISCHG(IR).EQ.1) THEN
 927           PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PMQTH1**2)
 928           PMTH(3,IR)=PMTH(2,IR)+PMQTH2
 929           PMTH(4,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)+PMTH(2,21)
 930           PMTH(5,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(83)**2)+PMTH(2,22)
 931         ELSEIF(ISCOL(IR).EQ.1) THEN
 932           PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)
 933           PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(82)
 934           PMTH(4,IR)=PMTH(3,IR)
 935           PMTH(5,IR)=PMTH(3,IR)
 936         ELSEIF(ISCHG(IR).EQ.1) THEN
 937           PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(90)**2)
 938           PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(90)
 939           PMTH(4,IR)=PMTH(3,IR)
 940           PMTH(5,IR)=PMTH(3,IR)
 941         ENDIF
 942         IF(KSH(IR).EQ.1) PMA(I)=PMTH(3,IR)
 943         PM=PM+PMA(I)
 944         IF(KSH(IR).EQ.0.OR.PMA(I).GT.10D0*QMAX) IREJ=IREJ+1
 945         DO 160 J=1,4
 946           PS(J)=PS(J)+P(IPA(I),J)
 947   160   CONTINUE
 948   170 CONTINUE
 949       IF(IREJ.EQ.NPA.AND.IP2.GE.-7) RETURN
 950       PS(5)=SQRT(MAX(0D0,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
 951       IF(NPA.EQ.1) PS(5)=PS(4)
 952       IF(PS(5).LE.PM+PMQT1E) RETURN
 953
 954 C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
 955       KFSRCE=0
 956       IF(IP2.LE.0) THEN
 957       ELSEIF(K(IP1,3).EQ.K(IP2,3).AND.K(IP1,3).GT.0) THEN
 958         KFSRCE=IABS(K(K(IP1,3),2))
 959       ELSE
 960         IPAR1=MAX(1,K(IP1,3))
 961         IPAR2=MAX(1,K(IP2,3))
 962         IF(K(IPAR1,3).EQ.K(IPAR2,3).AND.K(IPAR1,3).GT.0)
 963      &       KFSRCE=IABS(K(K(IPAR1,3),2))
 964       ENDIF
 965       ITYPES=0
 966       IF(KFSRCE.GE.1.AND.KFSRCE.LE.8) ITYPES=1
 967       IF(KFSRCE.GE.KSUSY1+1.AND.KFSRCE.LE.KSUSY1+8) ITYPES=2
 968       IF(KFSRCE.GE.KSUSY2+1.AND.KFSRCE.LE.KSUSY2+8) ITYPES=2
 969       IF(KFSRCE.GE.21.AND.KFSRCE.LE.24) ITYPES=3
 970       IF(KFSRCE.GE.32.AND.KFSRCE.LE.34) ITYPES=3
 971       IF(KFSRCE.EQ.25.OR.(KFSRCE.GE.35.AND.KFSRCE.LE.37)) ITYPES=4
 972       IF(KFSRCE.GE.KSUSY1+22.AND.KFSRCE.LE.KSUSY1+37) ITYPES=5
 973       IF(KFSRCE.EQ.KSUSY1+21) ITYPES=6
 974
 975 C...Identify two primary showerers.
 976       ITYPE1=0
 977       IF(KFLA(1).GE.1.AND.KFLA(1).LE.8) ITYPE1=1
 978       IF(KFLA(1).GE.KSUSY1+1.AND.KFLA(1).LE.KSUSY1+8) ITYPE1=2
 979       IF(KFLA(1).GE.KSUSY2+1.AND.KFLA(1).LE.KSUSY2+8) ITYPE1=2
 980       IF(KFLA(1).GE.21.AND.KFLA(1).LE.24) ITYPE1=3
 981       IF(KFLA(1).GE.32.AND.KFLA(1).LE.34) ITYPE1=3
 982       IF(KFLA(1).EQ.25.OR.(KFLA(1).GE.35.AND.KFLA(1).LE.37)) ITYPE1=4
 983       IF(KFLA(1).GE.KSUSY1+22.AND.KFLA(1).LE.KSUSY1+37) ITYPE1=5
 984       IF(KFLA(1).EQ.KSUSY1+21) ITYPE1=6
 985       ITYPE2=0
 986       IF(KFLA(2).GE.1.AND.KFLA(2).LE.8) ITYPE2=1
 987       IF(KFLA(2).GE.KSUSY1+1.AND.KFLA(2).LE.KSUSY1+8) ITYPE2=2
 988       IF(KFLA(2).GE.KSUSY2+1.AND.KFLA(2).LE.KSUSY2+8) ITYPE2=2
 989       IF(KFLA(2).GE.21.AND.KFLA(2).LE.24) ITYPE2=3
 990       IF(KFLA(2).GE.32.AND.KFLA(2).LE.34) ITYPE2=3
 991       IF(KFLA(2).EQ.25.OR.(KFLA(2).GE.35.AND.KFLA(2).LE.37)) ITYPE2=4
 992       IF(KFLA(2).GE.KSUSY1+22.AND.KFLA(2).LE.KSUSY1+37) ITYPE2=5
 993       IF(KFLA(2).EQ.KSUSY1+21) ITYPE2=6
 994
 995 C...Order of showerers. Presence of gluino.
 996       ITYPMN=MIN(ITYPE1,ITYPE2)
 997       ITYPMX=MAX(ITYPE1,ITYPE2)
 998       IORD=1
 999       IF(ITYPE1.GT.ITYPE2) IORD=2
1000       IGLUI=0
1001       IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1
1002
1003 C...Check if 3-jet matrix elements to be used.
1004       M3JC=0
1005       ALPHA=0.5D0
1006       IF(NPA.EQ.2.AND.MSTJ(47).GE.1.AND.MPSPD.EQ.0) THEN
1007         IF(MSTJ(38).NE.0) THEN
1008           M3JC=MSTJ(38)
1009           ALPHA=PARJ(80)
1010           MSTJ(38)=0
1011         ELSEIF(MSTJ(47).GE.6) THEN
1012           M3JC=MSTJ(47)
1013         ELSE
1014           ICLASS=1
1015           ICOMBI=4
1016
1017 C...Vector/axial vector -> q + qbar; q -> q + V.
1018           IF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.(ITYPES.EQ.0.OR.
1019      &    ITYPES.EQ.3)) THEN
1020             ICLASS=2
1021             IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN
1022               ICOMBI=1
1023             ELSEIF(KFSRCE.EQ.23.OR.(KFSRCE.EQ.0.AND.
1024      &      K(IPA(1),2)+K(IPA(2),2).EQ.0)) THEN
1025 C...gamma*/Z0: assume e+e- initial state if unknown.
1026               EI=-1D0
1027               IF(KFSRCE.EQ.23) THEN
1028                 IANNFL=K(K(IP1,3),3)
1029                 IF(IANNFL.NE.0) THEN
1030                   KANNFL=IABS(K(IANNFL,2))
1031                   IF(KANNFL.GE.1.AND.KANNFL.LE.18) EI=KCHG(KANNFL,1)/3D0
1032                 ENDIF
1033               ENDIF
1034               AI=SIGN(1D0,EI+0.1D0)
1035               VI=AI-4D0*EI*PARU(102)
1036               EF=KCHG(KFLA(1),1)/3D0
1037               AF=SIGN(1D0,EF+0.1D0)
1038               VF=AF-4D0*EF*PARU(102)
1039               XWC=1D0/(16D0*PARU(102)*(1D0-PARU(102)))
1040               SH=PS(5)**2
1041               SQMZ=PMAS(23,1)**2
1042               SQWZ=PS(5)*PMAS(23,2)
1043               SBWZ=1D0/((SH-SQMZ)**2+SQWZ**2)
1044               VECT=EI**2*EF**2+2D0*EI*VI*EF*VF*XWC*SH*(SH-SQMZ)*SBWZ+
1045      &        (VI**2+AI**2)*VF**2*XWC**2*SH**2*SBWZ
1046               AXIV=(VI**2+AI**2)*AF**2*XWC**2*SH**2*SBWZ
1047               ICOMBI=3
1048               ALPHA=VECT/(VECT+AXIV)
1049             ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN
1050               ICOMBI=4
1051             ENDIF
1052 C...For chi -> chi q qbar, use V/A -> q qbar as first approximation.
1053           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.5) THEN
1054             ICLASS=2
1055           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
1056      &    ITYPES.EQ.1)) THEN
1057             ICLASS=3
1058
1059 C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
1060           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN
1061             ICLASS=4
1062             IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN
1063               ICOMBI=1
1064             ELSEIF(KFSRCE.EQ.36) THEN
1065               ICOMBI=2
1066             ENDIF
1067           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
1068      &    ITYPES.EQ.1)) THEN
1069             ICLASS=5
1070
1071 C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
1072           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
1073      &    ITYPES.EQ.3)) THEN
1074             ICLASS=6
1075           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
1076      &    ITYPES.EQ.2)) THEN
1077             ICLASS=7
1078           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN
1079             ICLASS=8
1080           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
1081      &    ITYPES.EQ.2)) THEN
1082             ICLASS=9
1083
1084 C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
1085           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
1086      &    ITYPES.EQ.5)) THEN
1087             ICLASS=10
1088           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
1089      &    ITYPES.EQ.2)) THEN
1090             ICLASS=11
1091           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
1092      &    ITYPES.EQ.1)) THEN
1093             ICLASS=12
1094
1095 C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
1096           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.6) THEN
1097             ICLASS=13
1098           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
1099      &    ITYPES.EQ.2)) THEN
1100             ICLASS=14
1101           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
1102      &    ITYPES.EQ.1)) THEN
1103             ICLASS=15
1104
1105 C...g -> ~g + ~g (eikonal approximation).
1106           ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN
1107             ICLASS=16
1108           ENDIF
1109           M3JC=5*ICLASS+ICOMBI
1110         ENDIF
1111       ENDIF
1112
1113 C...Find if interference with initial state partons.
1114       MIIS=0
1115       IF(MSTJ(50).GE.1.AND.MSTJ(50).LE.3.AND.NPA.EQ.2.AND.KFSRCE.EQ.0
1116      &.AND.MPSPD.EQ.0) MIIS=MSTJ(50)
1117       IF(MSTJ(50).GE.4.AND.MSTJ(50).LE.6.AND.NPA.EQ.2.AND.MPSPD.EQ.0)
1118      &MIIS=MSTJ(50)-3
1119       IF(MIIS.NE.0) THEN
1120         DO 190 I=1,2
1121           KCII(I)=0
1122           KCA=PYCOMP(KFLA(I))
1123           IF(KCA.NE.0) KCII(I)=KCHG(KCA,2)*ISIGN(1,K(IPA(I),2))
1124           NIIS(I)=0
1125           IF(KCII(I).NE.0) THEN
1126             DO 180 J=1,2
1127               ICSI=MOD(K(IPA(I),3+J)/MSTU(5),MSTU(5))
1128               IF(ICSI.GT.0.AND.ICSI.NE.IPA(1).AND.ICSI.NE.IPA(2).AND.
1129      &        (KCII(I).EQ.(-1)**(J+1).OR.KCII(I).EQ.2)) THEN
1130                 NIIS(I)=NIIS(I)+1
1131                 IIIS(I,NIIS(I))=ICSI
1132               ENDIF
1133   180       CONTINUE
1134           ENDIF
1135   190   CONTINUE
1136         IF(NIIS(1)+NIIS(2).EQ.0) MIIS=0
1137       ENDIF
1138
1139 C...Boost interfering initial partons to rest frame
1140 C...and reconstruct their polar and azimuthal angles.
1141 Cacs+
1142         qplta1=0.d0
1143         qplta2=0.d0
1144         qpltbx=0.d0
1145         qpltby=0.d0
1146         qpltbz=0.d0
1147 Cacs-
1148       IF(MIIS.NE.0) THEN
1149         DO 210 I=1,2
1150           DO 200 J=1,5
1151             K(N+I,J)=K(IPA(I),J)
1152             P(N+I,J)=P(IPA(I),J)
1153             V(N+I,J)=0D0
1154   200     CONTINUE
1155   210   CONTINUE
1156         DO 230 I=3,2+NIIS(1)
1157           DO 220 J=1,5
1158             K(N+I,J)=K(IIIS(1,I-2),J)
1159             P(N+I,J)=P(IIIS(1,I-2),J)
1160             V(N+I,J)=0D0
1161   220     CONTINUE
1162   230   CONTINUE
1163         DO 250 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
1164           DO 240 J=1,5
1165             K(N+I,J)=K(IIIS(2,I-2-NIIS(1)),J)
1166             P(N+I,J)=P(IIIS(2,I-2-NIIS(1)),J)
1167             V(N+I,J)=0D0
1168   240     CONTINUE
1169   250   CONTINUE
1170         CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,0D0,-PS(1)/PS(4),
1171      &  -PS(2)/PS(4),-PS(3)/PS(4))
1172         PHI=PYANGL(P(N+1,1),P(N+1,2))
1173         CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,-PHI,0D0,0D0,0D0)
1174         THE=PYANGL(P(N+1,3),P(N+1,1))
1175         CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),-THE,0D0,0D0,0D0,0D0)
1176 Cacs+
1177         qplta1=-the
1178         qplta2=-phi
1179         qpltbx=-PS(1)/PS(4)
1180         qpltby=-PS(2)/PS(4)
1181         qpltbz=-PS(3)/PS(4)
1182 Cacs-
1183         DO 260 I=3,2+NIIS(1)
1184           THEIIS(1,I-2)=PYANGL(P(N+I,3),SQRT(P(N+I,1)**2+P(N+I,2)**2))
1185           PHIIIS(1,I-2)=PYANGL(P(N+I,1),P(N+I,2))
1186   260   CONTINUE
1187         DO 270 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
1188           THEIIS(2,I-2-NIIS(1))=PARU(1)-PYANGL(P(N+I,3),
1189      &    SQRT(P(N+I,1)**2+P(N+I,2)**2))
1190           PHIIIS(2,I-2-NIIS(1))=PYANGL(P(N+I,1),P(N+I,2))
1191   270   CONTINUE
1192       ENDIF
1193
1194 C...Boost 3 or more partons to their rest frame.
1195 Cacs+
1196 c      IF(NPA.GE.3) CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
1197 c     &-PS(2)/PS(4),-PS(3)/PS(4))
1198       IF(NPA.GE.3) THEN
1199         CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
1200      &-PS(2)/PS(4),-PS(3)/PS(4))
1201         qplta1=0.d0
1202         qplta2=0.d0
1203         qpltbx=-PS(1)/PS(4)
1204         qpltby=-PS(2)/PS(4)
1205         qpltbz=-PS(3)/PS(4)
1206
1207       ENDIF
1208 Cacs-
1209
1210 C...Define imagined single initiator of shower for parton system.
1211       NS=N
1212       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
1213         CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
1214         IF(MSTU(21).GE.1) RETURN
1215       ENDIF
1216   280 N=NS
1217       IF(NPA.GE.2) THEN
1218         K(N+1,1)=11
1219         K(N+1,2)=21
1220         K(N+1,3)=0
1221         K(N+1,4)=0
1222         K(N+1,5)=0
1223         P(N+1,1)=0D0
1224         P(N+1,2)=0D0
1225         P(N+1,3)=0D0
1226         P(N+1,4)=PS(5)
1227         P(N+1,5)=PS(5)
1228         V(N+1,5)=PS(5)**2
1229         N=N+1
1230         IREF(1)=21
1231       ENDIF
1232
1233
1234
1235
1236 Cacs+
1237       call qpygin(pposx0,pposy0,pposz0,ppost0) ! in fm
1238       do 10101 iijj=1, nnpos, 1
1239          ppos(iijj,1)=pposx0
1240          ppos(iijj,2)=pposy0
1241          ppos(iijj,3)=pposz0
1242          ppos(iijj,4)=ppost0
1243 10101 continue
1244
1245 Cacs-
1246
1247
1248
1249
1250 C...Loop over partons that may branch.
1251       NEP=NPA
1252       IM=NS
1253       IF(NPA.EQ.1) IM=NS-1
1254   290 IM=IM+1
1255       IF(N.GT.NS) THEN
1256         IF(IM.GT.N) GOTO 600
1257         KFLM=IABS(K(IM,2))
1258         IR=IREF(IM-NS)
1259         IF(KSH(IR).EQ.0) GOTO 290
1260         IF(P(IM,5).LT.PMTH(2,IR)) GOTO 290
1261         IGM=K(IM,3)
1262       ELSE
1263         IGM=-1
1264       ENDIF
1265       IF(N+NEP.GT.MSTU(4)-MSTU(32)-10) THEN
1266         CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
1267         IF(MSTU(21).GE.1) RETURN
1268       ENDIF
1269
1270 C...Position of aunt (sister to branching parton).
1271 C...Origin and flavour of daughters.
1272       IAU=0
1273       IF(IGM.GT.0) THEN
1274         IF(K(IM-1,3).EQ.IGM) IAU=IM-1
1275         IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
1276       ENDIF
1277       IF(IGM.GE.0) THEN
1278         K(IM,4)=N+1
1279         DO 300 I=1,NEP
1280           K(N+I,3)=IM
1281   300   CONTINUE
1282       ELSE
1283         K(N+1,3)=IPA(1)
1284       ENDIF
1285       IF(IGM.LE.0) THEN
1286         DO 310 I=1,NEP
1287           K(N+I,2)=K(IPA(I),2)
1288   310   CONTINUE
1289       ELSEIF(KFLM.NE.21) THEN
1290         K(N+1,2)=K(IM,2)
1291         K(N+2,2)=K(IM,5)
1292         IREF(N+1-NS)=IREF(IM-NS)
1293         IREF(N+2-NS)=IABS(K(N+2,2))
1294       ELSEIF(K(IM,5).EQ.21) THEN
1295         K(N+1,2)=21
1296         K(N+2,2)=21
1297         IREF(N+1-NS)=21
1298         IREF(N+2-NS)=21
1299       ELSE
1300         K(N+1,2)=K(IM,5)
1301         K(N+2,2)=-K(IM,5)
1302         IREF(N+1-NS)=IABS(K(N+1,2))
1303         IREF(N+2-NS)=IABS(K(N+2,2))
1304       ENDIF
1305
1306 C...Reset flags on daughters and tries made.
1307       DO 320 IP=1,NEP
1308         K(N+IP,1)=3
1309         K(N+IP,4)=0
1310         K(N+IP,5)=0
1311         KFLD(IP)=IABS(K(N+IP,2))
1312         IF(KCHG(PYCOMP(KFLD(IP)),2).EQ.0) K(N+IP,1)=1
1313         ITRY(IP)=0
1314         ISL(IP)=0
1315         ISI(IP)=0
1316         IF(KSH(IREF(N+IP-NS)).EQ.1) ISI(IP)=1
1317   320 CONTINUE
1318       ISLM=0
1319
1320 C...Maximum virtuality of daughters.
1321       IF(IGM.LE.0) THEN
1322         DO 330 I=1,NPA
1323           IF(NPA.GE.3) P(N+I,4)=P(IPA(I),4)
1324           P(N+I,5)=MIN(QMAX,PS(5))
1325           IR=IREF(N+I-NS)
1326           IF(IP2.LE.-8) P(N+I,5)=MAX(P(N+I,5),2D0*PMTH(3,IR))
1327           IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
1328   330   CONTINUE
1329       ELSE
1330         IF(MSTJ(43).LE.2) PEM=V(IM,2)
1331         IF(MSTJ(43).GE.3) PEM=P(IM,4)
1332         P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
1333         P(N+2,5)=MIN(P(IM,5),(1D0-V(IM,1))*PEM)
1334         IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
1335       ENDIF
1336       DO 340 I=1,NEP
1337         PMSD(I)=P(N+I,5)
1338         IF(ISI(I).EQ.1) THEN
1339           IR=IREF(N+I-NS)
1340           IF(P(N+I,5).LE.PMTH(3,IR)) P(N+I,5)=PMTH(1,IR)
1341         ENDIF
1342         V(N+I,5)=P(N+I,5)**2
1343   340 CONTINUE
1344
1345 C...Choose one of the daughters for evolution.
1346   350 INUM=0
1347       IF(NEP.EQ.1) INUM=1
1348       DO 360 I=1,NEP
1349         IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
1350   360 CONTINUE
1351       DO 370 I=1,NEP
1352         IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
1353           IR=IREF(N+I-NS)
1354           IF(P(N+I,5).GE.PMTH(2,IR)) INUM=I
1355         ENDIF
1356   370 CONTINUE
1357       IF(INUM.EQ.0) THEN
1358         RMAX=0D0
1359         DO 380 I=1,NEP
1360           IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQT2E) THEN
1361             RPM=P(N+I,5)/PMSD(I)
1362             IR=IREF(N+I-NS)
1363             IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,IR)) THEN
1364               RMAX=RPM
1365               INUM=I
1366             ENDIF
1367           ENDIF
1368   380   CONTINUE
1369       ENDIF
1370
1371 C...Cancel choice of predetermined daughter already treated.
1372       INUM=MAX(1,INUM)
1373       INUMT=INUM
1374       IF(MPSPD.EQ.1.AND.IGM.EQ.0.AND.ITRY(INUMT).GE.1) THEN
1375         IF(K(IP1-1+INUM,4).GT.0) INUM=3-INUM
1376       ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2.AND.ITRY(INUMT).GE.1) THEN
1377         IF(KFLD(INUMT).NE.21.AND.K(IP1+2,4).GT.0) INUM=3-INUM
1378         IF(KFLD(INUMT).EQ.21.AND.K(IP1+3,4).GT.0) INUM=3-INUM
1379       ENDIF
1380
1381 C...Store information on choice of evolving daughter.
1382       IEP(1)=N+INUM
1383 Cacs+
1384       idf=k(iep(1),3)
1385       zz1=v(idf,1)
1386       zzz=zz1
1387       zz2=1.d0-zz1
1388       if (nep .gt. 1 .and. inum .eq. 2) then
1389          zzz=zz2
1390       endif
1391       ttt=v(idf,5)
1392       if(zz1.gt.0.d0) then
1393             eee=zzz*p(idf,4)
1394       else
1395             eee=p(idf,4)
1396       endif
1397       xkt=zz1*zz2*ttt
1398       if (xkt .gt. 0.d0) then
1399          xlcoh=(2.d0*eee/(zz1*zz2*ttt))*0.1973d0
1400       else
1401          xlcoh=0.d0
1402       endif
1403       if (idf .eq. 0) then ! for the initial parton if it has no father
1404          xbx=p(iep(1),1)/p(iep(1),4)
1405          xby=p(iep(1),2)/p(iep(1),4)
1406          xbz=p(iep(1),3)/p(iep(1),4)
1407          call qpygeo(pposx0,pposy0,pposz0,ppost0,
1408      >               xbx,xby,xbz,qhatl,omegac)
1409       else
1410          xbx=p(idf,1)/p(idf,4)
1411          xby=p(idf,2)/p(idf,4)
1412          xbz=p(idf,3)/p(idf,4)
1413
1414
1415
1416          ppos(iep(1),1)=ppos(idf,1)+xbx*xlcoh
1417          ppos(iep(1),2)=ppos(idf,2)+xby*xlcoh
1418          ppos(iep(1),3)=ppos(idf,3)+xbz*xlcoh
1419          ppos(iep(1),4)=ppos(idf,4)+xlcoh
1420          call qpygeo(ppos(iep(1),1),ppos(iep(1),2),ppos(iep(1),3),
1421      >               ppos(iep(1),4),xbx,xby,xbz,qhatl,omegac)
1422       endif
1423 Cacs-
1424
1425       DO 390 I=2,NEP
1426         IEP(I)=IEP(I-1)+1
1427         IF(IEP(I).GT.N+NEP) IEP(I)=N+1
1428   390 CONTINUE
1429       DO 400 I=1,NEP
1430         KFL(I)=IABS(K(IEP(I),2))
1431   400 CONTINUE
1432       ITRY(INUM)=ITRY(INUM)+1
1433       IF(ITRY(INUM).GT.200) THEN
1434         CALL PYERRM(14,'(PYSHOW:) caught in infinite loop')
1435         IF(MSTU(21).GE.1) RETURN
1436       ENDIF
1437       Z=0.5D0
1438       IR=IREF(IEP(1)-NS)
1439       IF(KSH(IR).EQ.0) GOTO 450
1440       IF(P(IEP(1),5).LT.PMTH(2,IR)) GOTO 450
1441
1442 C...Check if evolution already predetermined for daughter.
1443       IPSPD=0
1444       IF(MPSPD.EQ.1.AND.IGM.EQ.0) THEN
1445         IF(K(IP1-1+INUM,4).GT.0) IPSPD=IP1-1+INUM
1446       ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2) THEN
1447         IF(KFL(1).NE.21.AND.K(IP1+2,4).GT.0) IPSPD=IP1+2
1448         IF(KFL(1).EQ.21.AND.K(IP1+3,4).GT.0) IPSPD=IP1+3
1449       ENDIF
1450       IF(INUM.EQ.1.OR.INUM.EQ.2) THEN
1451         ISSET(INUM)=0
1452         IF(IPSPD.NE.0) ISSET(INUM)=1
1453       ENDIF
1454
1455 C...Select side for interference with initial state partons.
1456       IF(MIIS.GE.1.AND.IEP(1).LE.NS+3) THEN
1457         III=IEP(1)-NS-1
1458         ISII(III)=0
1459         IF(IABS(KCII(III)).EQ.1.AND.NIIS(III).EQ.1) THEN
1460           ISII(III)=1
1461         ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.1) THEN
1462           IF(PYR(0).GT.0.5D0) ISII(III)=1
1463         ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.2) THEN
1464           ISII(III)=1
1465           IF(PYR(0).GT.0.5D0) ISII(III)=2
1466         ENDIF
1467       ENDIF
1468
1469 C...Calculate allowed z range.
1470       IF(NEP.EQ.1) THEN
1471         PMED=PS(4)
1472       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1473         PMED=P(IM,5)
1474       ELSE
1475         IF(INUM.EQ.1) PMED=V(IM,1)*PEM
1476         IF(INUM.EQ.2) PMED=(1D0-V(IM,1))*PEM
1477       ENDIF
1478       IF(MOD(MSTJ(43),2).EQ.1) THEN
1479         ZC=PMTH(2,21)/PMED
1480         ZCE=PMTH(2,22)/PMED
1481         IF(ISCOL(IR).EQ.0) ZCE=0.5D0*PARJ(90)/PMED
1482       ELSE
1483         ZC=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTH(2,21)/PMED)**2)))
1484         IF(ZC.LT.1D-6) ZC=(PMTH(2,21)/PMED)**2
1485         PMTMPE=PMTH(2,22)
1486         IF(ISCOL(IR).EQ.0) PMTMPE=0.5D0*PARJ(90)
1487         ZCE=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTMPE/PMED)**2)))
1488         IF(ZCE.LT.1D-6) ZCE=(PMTMPE/PMED)**2
1489       ENDIF
1490       ZC=MIN(ZC,0.491D0)
1491       ZCE=MIN(ZCE,0.49991D0)
1492       IF(((MSTJ(41).EQ.1.AND.ZC.GT.0.49D0).OR.(MSTJ(41).GE.2.AND.
1493      &MIN(ZC,ZCE).GT.0.4999D0)).AND.IPSPD.EQ.0) THEN
1494         P(IEP(1),5)=PMTH(1,IR)
1495         V(IEP(1),5)=P(IEP(1),5)**2
1496         GOTO 450
1497       ENDIF
1498
1499 C...Integral of Altarelli-Parisi z kernel for QCD.
1500 C...(Includes squark and gluino; with factor N_C/C_F extra for latter).
1501       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
1502 Cacs+
1503 C      FBR= 6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
1504       FBR=dgauss1(splitg1,zc,1.d0-zc,1.d-3)
1505 Cacs-
1506 C...QUARKONIA+++
1507 C...Evolution of QQ~[3S18] state if MSTP(148)=1.
1508       ELSEIF(MSTJ(49).EQ.0.AND.MSTP(149).GE.0.AND.
1509      &       (KFL(1).EQ.9900443.OR.KFL(1).EQ.9900553)) THEN
1510         FBR=6D0*LOG((1D0-ZC)/ZC)
1511 C...QUARKONIA---
1512       ELSEIF(MSTJ(49).EQ.0) THEN
1513 Cacs+
1514 C      FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC)
1515       FBR=dgauss1(splitq1,zc,1.d0-zc,1.d-3)
1516
1517 Cacs-
1518         IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0)
1519
1520 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
1521       ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
1522         FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1D0-2D0*ZC)
1523       ELSEIF(MSTJ(49).EQ.1) THEN
1524         FBR=(1D0-2D0*ZC)/3D0
1525         IF(IGM.EQ.0.AND.M3JC.GE.1) FBR=4D0*FBR
1526
1527 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
1528       ELSEIF(KFL(1).EQ.21) THEN
1529         FBR=6D0*MSTJ(45)*(0.5D0-ZC)
1530       ELSE
1531         FBR=2D0*LOG((1D0-ZC)/ZC)
1532       ENDIF
1533
1534 C...Reset QCD probability for colourless.
1535       IF(ISCOL(IR).EQ.0) FBR=0D0
1536
1537 C...Integral of Altarelli-Parisi kernel for photon emission.
1538       FBRE=0D0
1539       IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1) THEN
1540         IF(KFL(1).LE.18) THEN
1541           FBRE=(KCHG(KFL(1),1)/3D0)**2*2D0*LOG((1D0-ZCE)/ZCE)
1542         ENDIF
1543         IF(MSTJ(41).EQ.10) FBRE=PARJ(84)*FBRE
1544       ENDIF
1545
1546 C...Inner veto algorithm starts. Find maximum mass for evolution.
1547   410 PMS=V(IEP(1),5)
1548       IF(IGM.GE.0) THEN
1549         PM2=0D0
1550         DO 420 I=2,NEP
1551           PM=P(IEP(I),5)
1552           IRI=IREF(IEP(I)-NS)
1553           IF(KSH(IRI).EQ.1) PM=PMTH(2,IRI)
1554           PM2=PM2+PM
1555   420   CONTINUE
1556         PMS=MIN(PMS,(P(IM,5)-PM2)**2)
1557       ENDIF
1558
1559 C...Select mass for daughter in QCD evolution.
1560       B0=27D0/6D0
1561       DO 430 IFF=4,MSTJ(45)
1562         IF(PMS.GT.4D0*PMTH(2,IFF)**2) B0=(33D0-2D0*IFF)/6D0
1563   430 CONTINUE
1564 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
1565       PMSC=MAX(0.5D0*PARJ(82),PMS-PMTH(1,IR)**2)
1566 C...Already predetermined choice.
1567       IF(IPSPD.NE.0) THEN
1568         PMSQCD=P(IPSPD,5)**2
1569       ELSEIF(FBR.LT.1D-3) THEN
1570         PMSQCD=0D0
1571       ELSEIF(MSTJ(44).LE.0) THEN
1572         PMSQCD=PMSC*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/(PARU(111)*FBR)))
1573       ELSEIF(MSTJ(44).EQ.1) THEN
1574         PMSQCD=4D0*ALAMS*(0.25D0*PMSC/ALAMS)**(PYR(0)**(B0/FBR))
1575       ELSE
1576         PMSQCD=PMSC*EXP(MAX(-50D0,ALFM*B0*LOG(PYR(0))/FBR))
1577       ENDIF
1578 C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
1579       IF(IPSPD.EQ.0) PMSQCD=PMSQCD+PMTH(1,IR)**2
1580       IF(ZC.GT.0.49D0.OR.PMSQCD.LE.PMTH(4,IR)**2) PMSQCD=PMTH(2,IR)**2
1581       V(IEP(1),5)=PMSQCD
1582       MCE=1
1583
1584 C...Select mass for daughter in QED evolution.
1585       IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1.AND.IPSPD.EQ.0) THEN
1586 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
1587         PMSE=MAX(0.5D0*PARJ(83),PMS-PMTH(1,IR)**2)
1588         IF(FBRE.LT.1D-3) THEN
1589           PMSQED=0D0
1590         ELSE
1591           PMSQED=PMSE*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/
1592      &    (PARU(101)*FBRE)))
1593         ENDIF
1594 C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
1595         PMSQED=PMSQED+PMTH(1,IR)**2
1596         IF(ZCE.GT.0.4999D0.OR.PMSQED.LE.PMTH(5,IR)**2) PMSQED=
1597      &  PMTH(2,IR)**2
1598         IF(PMSQED.GT.PMSQCD) THEN
1599           V(IEP(1),5)=PMSQED
1600           MCE=2
1601         ENDIF
1602       ENDIF
1603
1604 C...Check whether daughter mass below cutoff.
1605       P(IEP(1),5)=SQRT(V(IEP(1),5))
1606       IF(P(IEP(1),5).LE.PMTH(3,IR)) THEN
1607         P(IEP(1),5)=PMTH(1,IR)
1608         V(IEP(1),5)=P(IEP(1),5)**2
1609         GOTO 450
1610       ENDIF
1611 Cacs+
1612        virt=V(IEP(1),5)
1613 Cacs-
1614
1615 C...Already predetermined choice of z, and flavour in g -> qqbar.
1616       IF(IPSPD.NE.0) THEN
1617         IPSGD1=K(IPSPD,4)
1618         IPSGD2=K(IPSPD,5)
1619         PMSGD1=P(IPSGD1,5)**2
1620         PMSGD2=P(IPSGD2,5)**2
1621         ALAMPS=SQRT(MAX(1D-10,(PMSQCD-PMSGD1-PMSGD2)**2-
1622      &  4D0*PMSGD1*PMSGD2))
1623         Z=0.5D0*(PMSQCD*(2D0*P(IPSGD1,4)/P(IPSPD,4)-1D0)+ALAMPS-
1624      &  PMSGD1+PMSGD2)/ALAMPS
1625         Z=MAX(0.00001D0,MIN(0.99999D0,Z))
1626         IF(KFL(1).NE.21) THEN
1627           K(IEP(1),5)=21
1628         ELSE
1629           K(IEP(1),5)=IABS(K(IPSGD1,2))
1630         ENDIF
1631
1632 C...Select z value of branching: q -> qgamma.
1633       ELSEIF(MCE.EQ.2) THEN
1634         Z=1D0-(1D0-ZCE)*(ZCE/(1D0-ZCE))**PYR(0)
1635         IF(1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
1636         K(IEP(1),5)=22
1637
1638 C...QUARKONIA+++
1639 C...Select z value of branching: QQ~[3S18] -> QQ~[3S18]g.
1640       ELSEIF(MSTJ(49).EQ.0.AND.
1641      &       (KFL(1).EQ.9900443.OR.KFL(1).EQ.9900553)) THEN
1642         Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1643 C...Select always the harder 'gluon' if the switch MSTP(149)<=0.
1644         IF(MSTP(149).LE.0.OR.PYR(0).GT.0.5D0) Z=1D0-Z
1645         IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
1646         K(IEP(1),5)=21
1647 C...QUARKONIA---
1648
1649 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
1650       ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
1651 Cacs+
1652 C        Z=1D0-(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1653 C...Only do z weighting when no ME correction afterwards.
1654 C        IF(M3JC.EQ.0.AND.1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
1655 C
1656         anfbr=dgauss1(splitq2,zc,1.d0-zc,1.d-3)
1657         z=simdis(500,zc,1,anfbr)
1658 Cacs-
1659         K(IEP(1),5)=21
1660       ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN
1661 Cacs+
1662 c        Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1663         anfbr=dgauss1(splitg2,zc,1.d0-zc,1.d-3)
1664
1665         z=simdis(500,zc,21,anfbr)
1666
1667         IF(PYR(0).GT.0.5D0) Z=1D0-Z
1668 c        IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
1669 Cacs-
1670         K(IEP(1),5)=21
1671       ELSEIF(MSTJ(49).NE.1) THEN
1672 Cacs+
1673 c        Z=PYR(0)
1674 c        IF(Z**2+(1D0-Z)**2.LT.PYR(0)) GOTO 410
1675         anfbr=dgauss1(splitqqbar,zc,1.d0-zc,1.d-3)
1676         z=simdis(500,zc,3,anfbr)
1677
1678 Cacs-
1679         KFLB=1+INT(MSTJ(45)*PYR(0))
1680         PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5)
1681         IF(PMQ.GE.1D0) GOTO 410
1682         IF(MSTJ(44).LE.2.OR.MSTJ(44).EQ.4) THEN
1683           IF(Z.LT.ZC.OR.Z.GT.1D0-ZC) GOTO 410
1684           PMQ0=4D0*PMTH(2,21)**2/V(IEP(1),5)
1685           IF(MOD(MSTJ(43),2).EQ.0.AND.(1D0+0.5D0*PMQ)*SQRT(1D0-PMQ)
1686      &    .LT.PYR(0)*(1D0+0.5D0*PMQ0)*SQRT(1D0-PMQ0)) GOTO 410
1687         ELSE
1688           IF((1D0+0.5D0*PMQ)*SQRT(1D0-PMQ).LT.PYR(0)) GOTO 410
1689         ENDIF
1690         K(IEP(1),5)=KFLB
1691
1692 C...Ditto for scalar gluon model.
1693       ELSEIF(KFL(1).NE.21) THEN
1694         Z=1D0-SQRT(ZC**2+PYR(0)*(1D0-2D0*ZC))
1695         K(IEP(1),5)=21
1696       ELSEIF(PYR(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
1697         Z=ZC+(1D0-2D0*ZC)*PYR(0)
1698         K(IEP(1),5)=21
1699       ELSE
1700         Z=ZC+(1D0-2D0*ZC)*PYR(0)
1701         KFLB=1+INT(MSTJ(45)*PYR(0))
1702         PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5)
1703         IF(PMQ.GE.1D0) GOTO 410
1704         K(IEP(1),5)=KFLB
1705       ENDIF
1706
1707 C...Correct to alpha_s(pT^2) (optionally m^2/4 for g -> q qbar).
1708       IF(MCE.EQ.1.AND.MSTJ(44).GE.2.AND.IPSPD.EQ.0) THEN
1709         IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND.
1710      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1711           IF(ALFM/LOG(V(IEP(1),5)*0.25D0/ALAMS).LT.PYR(0)) GOTO 410
1712         ELSE
1713           PT2APP=Z*(1D0-Z)*V(IEP(1),5)
1714           IF(MSTJ(44).GE.4) PT2APP=PT2APP*
1715      &    (1D0-PMTH(1,IR)**2/V(IEP(1),5))**2
1716           IF(PT2APP.LT.PT2MIN) GOTO 410
1717           IF(ALFM/LOG(PT2APP/ALAMS).LT.PYR(0)) GOTO 410
1718         ENDIF
1719       ENDIF
1720
1721 C...Check if z consistent with chosen m.
1722       IF(KFL(1).EQ.21) THEN
1723         IRGD1=IABS(K(IEP(1),5))
1724         IRGD2=IRGD1
1725       ELSE
1726         IRGD1=IR
1727         IRGD2=IABS(K(IEP(1),5))
1728       ENDIF
1729       IF(NEP.EQ.1) THEN
1730         PED=PS(4)
1731       ELSEIF(NEP.GE.3) THEN
1732         PED=P(IEP(1),4)
1733       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1734         PED=0.5D0*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
1735       ELSE
1736         IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
1737         IF(IEP(1).EQ.N+2) PED=(1D0-V(IM,1))*PEM
1738       ENDIF
1739       IF(MOD(MSTJ(43),2).EQ.1) THEN
1740         PMQTH3=0.5D0*PARJ(82)
1741         IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83)
1742         IF(IRGD2.EQ.22.AND.ISCOL(IR).EQ.0) PMQTH3=0.5D0*PARJ(90)
1743         PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(IEP(1),5)
1744         PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(IEP(1),5)
1745         ZD=SQRT(MAX(0D0,(1D0-V(IEP(1),5)/PED**2)*((1D0-PMQ1-PMQ2)**2-
1746      &  4D0*PMQ1*PMQ2)))
1747         ZH=1D0+PMQ1-PMQ2
1748       ELSE
1749         ZD=SQRT(MAX(0D0,1D0-V(IEP(1),5)/PED**2))
1750         ZH=1D0
1751       ENDIF
1752       IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND.
1753      &(MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1754       ELSEIF(IPSPD.NE.0) THEN
1755       ELSE
1756         ZL=0.5D0*(ZH-ZD)
1757         ZU=0.5D0*(ZH+ZD)
1758         IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 410
1759       ENDIF
1760       IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1D0-ZL)/MAX(1D-20,ZL*
1761      &(1D0-ZU)))
1762       IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
1763
1764 C...Width suppression for q -> q + g.
1765       IF(MSTJ(40).NE.0.AND.KFL(1).NE.21.AND.IPSPD.EQ.0) THEN
1766         IF(IGM.EQ.0) THEN
1767           EGLU=0.5D0*PS(5)*(1D0-Z)*(1D0+V(IEP(1),5)/V(NS+1,5))
1768         ELSE
1769           EGLU=PMED*(1D0-Z)
1770         ENDIF
1771         CHI=PARJ(89)**2/(PARJ(89)**2+EGLU**2)
1772         IF(MSTJ(40).EQ.1) THEN
1773           IF(CHI.LT.PYR(0)) GOTO 410
1774         ELSEIF(MSTJ(40).EQ.2) THEN
1775           IF(1D0-CHI.LT.PYR(0)) GOTO 410
1776         ENDIF
1777       ENDIF
1778
1779 C...Three-jet matrix element correction.
1780       IF(M3JC.GE.1) THEN
1781         WME=1D0
1782         WSHOW=1D0
1783
1784 C...QED matrix elements: only for massless case so far.
1785         IF(MCE.EQ.2.AND.IGM.EQ.0) THEN
1786           X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5))
1787           X2=1D0-V(IEP(1),5)/V(NS+1,5)
1788           X3=(1D0-X1)+(1D0-X2)
1789           KI1=K(IPA(INUM),2)
1790           KI2=K(IPA(3-INUM),2)
1791           QF1=KCHG(PYCOMP(KI1),1)*ISIGN(1,KI1)/3D0
1792           QF2=KCHG(PYCOMP(KI2),1)*ISIGN(1,KI2)/3D0
1793           WSHOW=QF1**2*(1D0-X1)/X3*(1D0+(X1/(2D0-X2))**2)+
1794      &    QF2**2*(1D0-X2)/X3*(1D0+(X2/(2D0-X1))**2)
1795           WME=(QF1*(1D0-X1)/X3-QF2*(1D0-X2)/X3)**2*(X1**2+X2**2)
1796         ELSEIF(MCE.EQ.2) THEN
1797
1798 C...QCD matrix elements, including mass effects.
1799         ELSEIF(MSTJ(49).NE.1.AND.K(IEP(1),2).NE.21) THEN
1800           PS1ME=V(IEP(1),5)
1801           PM1ME=PMTH(1,IR)
1802           M3JCC=M3JC
1803           IF(IR.GE.31.AND.IGM.EQ.0) THEN
1804 C...QCD ME: original parton, first branching.
1805             PM2ME=PMTH(1,63-IR)
1806             ECMME=PS(5)
1807           ELSEIF(IR.GE.31) THEN
1808 C...QCD ME: original parton, subsequent branchings.
1809             PM2ME=PMTH(1,63-IR)
1810             PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5))
1811             ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1812           ELSEIF(K(IM,2).EQ.21) THEN
1813 C...QCD ME: secondary partons, first branching.
1814             PM2ME=PM1ME
1815             ZMME=V(IM,1)
1816             IF(IEP(1).GT.IEP(2)) ZMME=1D0-ZMME
1817             PMLME=SQRT(MAX(0D0,(V(IM,5)-PS1ME-PM2ME**2)**2-
1818      &      4D0*PS1ME*PM2ME**2))
1819             PEDME=PEM*(0.5D0*(V(IM,5)-PMLME+PS1ME-PM2ME**2)+PMLME*ZMME)/
1820      &      V(IM,5)
1821             ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1822             M3JCC=66
1823           ELSE
1824 C...QCD ME: secondary partons, subsequent branchings.
1825             PM2ME=PM1ME
1826             PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5))
1827             ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1828             M3JCC=66
1829           ENDIF
1830 C...Construct ME variables.
1831           R1ME=PM1ME/ECMME
1832           R2ME=PM2ME/ECMME
1833           X1=(1D0+PS1ME/ECMME**2-R2ME**2)*(Z+(1D0-Z)*PM1ME**2/PS1ME)
1834           X2=1D0+R2ME**2-PS1ME/ECMME**2
1835 C...Call ME, with right order important for two inequivalent showerers.
1836           IF(IR.EQ.IORD+30) THEN
1837             WME=PYMAEL(M3JCC,X1,X2,R1ME,R2ME,ALPHA)
1838           ELSE
1839             WME=PYMAEL(M3JCC,X2,X1,R2ME,R1ME,ALPHA)
1840           ENDIF
1841 C...Split up total ME when two radiating partons.
1842           ISPRAD=1
1843           IF((M3JCC.GE.16.AND.M3JCC.LE.19).OR.
1844      &    (M3JCC.GE.26.AND.M3JCC.LE.29).OR.
1845      &    (M3JCC.GE.36.AND.M3JCC.LE.39).OR.
1846      &    (M3JCC.GE.46.AND.M3JCC.LE.49).OR.
1847      &    (M3JCC.GE.56.AND.M3JCC.LE.64)) ISPRAD=0
1848           IF(ISPRAD.EQ.1) WME=WME*MAX(1D-10,1D0+R1ME**2-R2ME**2-X1)/
1849      &    MAX(1D-10,2D0-X1-X2)
1850 C...Evaluate shower rate to be compared with.
1851           WSHOW=2D0/(MAX(1D-10,2D0-X1-X2)*
1852      &    MAX(1D-10,1D0+R2ME**2-R1ME**2-X2))
1853           IF(IGLUI.EQ.1.AND.IR.GE.31) WSHOW=(9D0/4D0)*WSHOW
1854         ELSEIF(MSTJ(49).NE.1) THEN
1855
1856 C...Toy model scalar theory matrix elements; no mass effects.
1857         ELSE
1858           X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5))
1859           X2=1D0-V(IEP(1),5)/V(NS+1,5)
1860           X3=(1D0-X1)+(1D0-X2)
1861           WSHOW=4D0*X3*((1D0-X1)/(2D0-X2)**2+(1D0-X2)/(2D0-X1)**2)
1862           WME=X3**2
1863           IF(MSTJ(102).GE.2) WME=X3**2-2D0*(1D0+X3)*(1D0-X1)*(1D0-X2)*
1864      &    PARJ(171)
1865         ENDIF
1866
1867         IF(WME.LT.PYR(0)*WSHOW) GOTO 410
1868       ENDIF
1869
1870 C...Impose angular ordering by rejection of nonordered emission.
1871       IF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2.AND.IPSPD.EQ.0) THEN
1872         PEMAO=V(IM,1)*P(IM,4)
1873         IF(IEP(1).EQ.N+2) PEMAO=(1D0-V(IM,1))*P(IM,4)
1874         IF(IR.GE.31.AND.MSTJ(42).GE.5) THEN
1875           MAOD=0
1876         ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.4
1877      &  .OR.MSTJ(42).EQ.7)) THEN
1878           MAOD=0
1879         ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.3
1880      &  .OR.MSTJ(42).EQ.6)) THEN
1881           MAOD=1
1882           PMDAO=PMTH(2,K(IEP(1),5))
1883           THE2ID=Z*(1D0-Z)*PEMAO**2/(V(IEP(1),5)-4D0*PMDAO**2)
1884         ELSE
1885           MAOD=1
1886           THE2ID=Z*(1D0-Z)*PEMAO**2/V(IEP(1),5)
1887           IF(MSTJ(42).GE.3.AND.MSTJ(42).NE.5) THE2ID=THE2ID*
1888      &    (1D0+PMTH(1,IR)**2*(1D0-Z)/(V(IEP(1),5)*Z))**2
1889         ENDIF
1890         MAOM=1
1891         IAOM=IM
1892   440   IF(K(IAOM,5).EQ.22) THEN
1893           IAOM=K(IAOM,3)
1894           IF(K(IAOM,3).LE.NS) MAOM=0
1895           IF(MAOM.EQ.1) GOTO 440
1896         ENDIF
1897         IF(MAOM.EQ.1.AND.MAOD.EQ.1) THEN
1898           THE2IM=V(IAOM,1)*(1D0-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
1899           IF(THE2ID.LT.THE2IM) GOTO 410
1900         ENDIF
1901       ENDIF
1902
1903 C...Impose user-defined maximum angle at first branching.
1904       IF(MSTJ(48).EQ.1.AND.IPSPD.EQ.0) THEN
1905         IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
1906           THE2ID=Z*(1D0-Z)*PS(4)**2/V(IEP(1),5)
1907           IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 410
1908         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
1909           THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5)
1910           IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 410
1911         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
1912           THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5)
1913           IF(PARJ(86)**2*THE2ID.LT.1D0) GOTO 410
1914         ENDIF
1915       ENDIF
1916
1917 C...Impose angular constraint in first branching from interference
1918 C...with initial state partons.
1919       IF(MIIS.GE.2.AND.IEP(1).LE.NS+3) THEN
1920         THE2D=MAX((1D0-Z)/Z,Z/(1D0-Z))*V(IEP(1),5)/(0.5D0*P(IM,4))**2
1921         IF(IEP(1).EQ.NS+2.AND.ISII(1).GE.1) THEN
1922           IF(THE2D.GT.THEIIS(1,ISII(1))**2) GOTO 410
1923         ELSEIF(IEP(1).EQ.NS+3.AND.ISII(2).GE.1) THEN
1924           IF(THE2D.GT.THEIIS(2,ISII(2))**2) GOTO 410
1925         ENDIF
1926       ENDIF
1927
1928 C...End of inner veto algorithm. Check if only one leg evolved so far.
1929   450 V(IEP(1),1)=Z
1930       ISL(1)=0
1931       ISL(2)=0
1932       IF(NEP.EQ.1) GOTO 490
1933       IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 350
1934       DO 460 I=1,NEP
1935         IR=IREF(N+I-NS)
1936         IF(ITRY(I).EQ.0.AND.KSH(IR).EQ.1) THEN
1937           IF(P(N+I,5).GE.PMTH(2,IR)) GOTO 350
1938         ENDIF
1939   460 CONTINUE
1940
1941 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
1942       IF(NEP.GE.3) THEN
1943         PMSUM=0D0
1944         DO 470 I=1,NEP
1945           PMSUM=PMSUM+P(N+I,5)
1946   470   CONTINUE
1947         IF(PMSUM.GE.PS(5)) GOTO 350
1948       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
1949         DO 480 I1=N+1,N+2
1950           IRDA=IREF(I1-NS)
1951           IF(KSH(IRDA).EQ.0) GOTO 480
1952           IF(P(I1,5).LT.PMTH(2,IRDA)) GOTO 480
1953           IF(IRDA.EQ.21) THEN
1954             IRGD1=IABS(K(I1,5))
1955             IRGD2=IRGD1
1956           ELSE
1957             IRGD1=IRDA
1958             IRGD2=IABS(K(I1,5))
1959           ENDIF
1960           I2=2*N+3-I1
1961           IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1962             PED=0.5D0*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
1963           ELSE
1964             IF(I1.EQ.N+1) ZM=V(IM,1)
1965             IF(I1.EQ.N+2) ZM=1D0-V(IM,1)
1966             PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
1967      &      4D0*V(N+1,5)*V(N+2,5))
1968             PED=PEM*(0.5D0*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/
1969      &      V(IM,5)
1970           ENDIF
1971           IF(MOD(MSTJ(43),2).EQ.1) THEN
1972             PMQTH3=0.5D0*PARJ(82)
1973             IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83)
1974             IF(IRGD2.EQ.22.AND.ISCOL(IRDA).EQ.0) PMQTH3=0.5D0*PARJ(90)
1975             PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(I1,5)
1976             PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(I1,5)
1977             ZD=SQRT(MAX(0D0,(1D0-V(I1,5)/PED**2)*((1D0-PMQ1-PMQ2)**2-
1978      &      4D0*PMQ1*PMQ2)))
1979             ZH=1D0+PMQ1-PMQ2
1980           ELSE
1981             ZD=SQRT(MAX(0D0,1D0-V(I1,5)/PED**2))
1982             ZH=1D0
1983           ENDIF
1984           IF(IRDA.EQ.21.AND.IRGD1.LT.10.AND.
1985      &    (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1986           ELSE
1987             ZL=0.5D0*(ZH-ZD)
1988             ZU=0.5D0*(ZH+ZD)
1989             IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND.
1990      &      ISSET(1).EQ.0) THEN
1991               ISL(1)=1
1992             ELSEIF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND.
1993      &      ISSET(2).EQ.0) THEN
1994               ISL(2)=1
1995             ENDIF
1996           ENDIF
1997           IF(IRDA.EQ.21) V(I1,4)=LOG(ZU*(1D0-ZL)/MAX(1D-20,
1998      &    ZL*(1D0-ZU)))
1999           IF(IRDA.NE.21) V(I1,4)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
2000   480   CONTINUE
2001         IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
2002           ISL(3-ISLM)=0
2003           ISLM=3-ISLM
2004         ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
2005           ZDR1=MAX(0D0,V(N+1,3)/MAX(1D-6,V(N+1,4))-1D0)
2006           ZDR2=MAX(0D0,V(N+2,3)/MAX(1D-6,V(N+2,4))-1D0)
2007           IF(ZDR2.GT.PYR(0)*(ZDR1+ZDR2)) ISL(1)=0
2008           IF(ISL(1).EQ.1) ISL(2)=0
2009           IF(ISL(1).EQ.0) ISLM=1
2010           IF(ISL(2).EQ.0) ISLM=2
2011         ENDIF
2012         IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 350
2013       ENDIF
2014       IRD1=IREF(N+1-NS)
2015       IRD2=IREF(N+2-NS)
2016       IF(IGM.GT.0) THEN
2017         IF(MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
2018      &  PMTH(2,IRD1).OR.P(N+2,5).GE.PMTH(2,IRD2))) THEN
2019           PMQ1=V(N+1,5)/V(IM,5)
2020           PMQ2=V(N+2,5)/V(IM,5)
2021           ZD=SQRT(MAX(0D0,(1D0-V(IM,5)/PEM**2)*((1D0-PMQ1-PMQ2)**2-
2022      &    4D0*PMQ1*PMQ2)))
2023           ZH=1D0+PMQ1-PMQ2
2024           ZL=0.5D0*(ZH-ZD)
2025           ZU=0.5D0*(ZH+ZD)
2026           IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 350
2027         ENDIF
2028       ENDIF
2029
2030 C...Accepted branch. Construct four-momentum for initial partons.
2031   490 MAZIP=0
2032       MAZIC=0
2033       IF(NEP.EQ.1) THEN
2034         P(N+1,1)=0D0
2035         P(N+1,2)=0D0
2036         P(N+1,3)=SQRT(MAX(0D0,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
2037      &  P(N+1,5))))
2038         P(N+1,4)=P(IPA(1),4)
2039         V(N+1,2)=P(N+1,4)
2040       ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
2041         PED1=0.5D0*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
2042         P(N+1,1)=0D0
2043         P(N+1,2)=0D0
2044         P(N+1,3)=SQRT(MAX(0D0,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
2045         P(N+1,4)=PED1
2046         P(N+2,1)=0D0
2047         P(N+2,2)=0D0
2048         P(N+2,3)=-P(N+1,3)
2049         P(N+2,4)=P(IM,5)-PED1
2050         V(N+1,2)=P(N+1,4)
2051         V(N+2,2)=P(N+2,4)
2052       ELSEIF(NEP.GE.3) THEN
2053 C...Rescale all momenta for energy conservation.
2054         LOOP=0
2055         PES=0D0
2056         PQS=0D0
2057         DO 510 I=1,NEP
2058           DO 500 J=1,4
2059             P(N+I,J)=P(IPA(I),J)
2060   500     CONTINUE
2061           PES=PES+P(N+I,4)
2062           PQS=PQS+P(N+I,5)**2/P(N+I,4)
2063   510   CONTINUE
2064   520   LOOP=LOOP+1
2065         FAC=(PS(5)-PQS)/(PES-PQS)
2066         PES=0D0
2067         PQS=0D0
2068         DO 540 I=1,NEP
2069           DO 530 J=1,3
2070             P(N+I,J)=FAC*P(N+I,J)
2071   530     CONTINUE
2072           P(N+I,4)=SQRT(P(N+I,5)**2+P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2)
2073           V(N+I,2)=P(N+I,4)
2074           PES=PES+P(N+I,4)
2075           PQS=PQS+P(N+I,5)**2/P(N+I,4)
2076   540   CONTINUE
2077         IF(LOOP.LT.10.AND.ABS(PES-PS(5)).GT.1D-12*PS(5)) GOTO 520
2078
2079 C...Construct transverse momentum for ordinary branching in shower.
2080       ELSE
2081         ZM=V(IM,1)
2082         LOOPPT=0
2083   550   LOOPPT=LOOPPT+1
2084         PZM=SQRT(MAX(0D0,(PEM+P(IM,5))*(PEM-P(IM,5))))
2085         PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4D0*V(N+1,5)*V(N+2,5)
2086         IF(PZM.LE.0D0) THEN
2087           PTS=0D0
2088         ELSEIF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
2089      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
2090           PTS=PMLS*ZM*(1D0-ZM)/V(IM,5)
2091         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
2092           PTS=(PEM**2*(ZM*(1D0-ZM)*V(IM,5)-(1D0-ZM)*V(N+1,5)-
2093      &    ZM*V(N+2,5))-0.25D0*PMLS)/PZM**2
2094         ELSE
2095           PTS=PMLS*(ZM*(1D0-ZM)*PEM**2/V(IM,5)-0.25D0)/PZM**2
2096         ENDIF
2097         IF(PTS.LT.0D0.AND.LOOPPT.LT.10) THEN
2098           ZM=0.05D0+0.9D0*ZM
2099           GOTO 550
2100         ELSEIF(PTS.LT.0D0) THEN
2101           GOTO 280
2102         ENDIF
2103         PT=SQRT(MAX(0D0,PTS))
2104
2105 C...Global statistics.
2106         MINT(353)=MINT(353)+1
2107         VINT(353)=VINT(353)+PT
2108         IF (MINT(353).EQ.1) VINT(358)=PT
2109
2110 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
2111         HAZIP=0D0
2112         IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21
2113      &  .AND.IAU.NE.0) THEN
2114           IF(K(IGM,3).NE.0) MAZIP=1
2115           ZAU=V(IGM,1)
2116           IF(IAU.EQ.IM+1) ZAU=1D0-V(IGM,1)
2117           IF(MAZIP.EQ.0) ZAU=0D0
2118           IF(K(IGM,2).NE.21) THEN
2119             HAZIP=2D0*ZAU/(1D0+ZAU**2)
2120           ELSE
2121             HAZIP=(ZAU/(1D0-ZAU*(1D0-ZAU)))**2
2122           ENDIF
2123           IF(K(N+1,2).NE.21) THEN
2124             HAZIP=HAZIP*(-2D0*ZM*(1D0-ZM))/(1D0-2D0*ZM*(1D0-ZM))
2125           ELSE
2126             HAZIP=HAZIP*(ZM*(1D0-ZM)/(1D0-ZM*(1D0-ZM)))**2
2127           ENDIF
2128         ENDIF
2129
2130 C...Find coefficient of azimuthal asymmetry due to soft gluon
2131 C...interference.
2132         HAZIC=0D0
2133         IF(MSTJ(49).NE.2.AND.MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.
2134      &  K(N+2,2).EQ.21).AND.IAU.NE.0) THEN
2135           IF(K(IGM,3).NE.0) MAZIC=N+1
2136           IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
2137           IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
2138      &    ZM.GT.0.5D0) MAZIC=N+2
2139           IF(K(IAU,2).EQ.22) MAZIC=0
2140           ZS=ZM
2141           IF(MAZIC.EQ.N+2) ZS=1D0-ZM
2142           ZGM=V(IGM,1)
2143           IF(IAU.EQ.IM-1) ZGM=1D0-V(IGM,1)
2144           IF(MAZIC.EQ.0) ZGM=1D0
2145           IF(MAZIC.NE.0) HAZIC=(P(IM,5)/P(IGM,5))*
2146      &    SQRT((1D0-ZS)*(1D0-ZGM)/(ZS*ZGM))
2147           HAZIC=MIN(0.95D0,HAZIC)
2148         ENDIF
2149       ENDIF
2150
2151 C...Construct energies for ordinary branching in shower.
2152   560 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
2153         IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
2154      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
2155           P(N+1,4)=0.5D0*(PEM*(V(IM,5)+V(N+1,5)-V(N+2,5))+
2156      &    PZM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5)
2157         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
2158           P(N+1,4)=PEM*V(IM,1)
2159         ELSE
2160           P(N+1,4)=PEM*(0.5D0*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
2161      &    SQRT(PMLS)*ZM)/V(IM,5)
2162         ENDIF
2163
2164 C...Already predetermined choice of phi angle or not
2165
2166         PHI=PARU(2)*PYR(0)
2167         IF(MPSPD.EQ.1.AND.IGM.EQ.NS+1) THEN
2168           IPSPD=IP1+IM-NS-2
2169           IF(K(IPSPD,4).GT.0) THEN
2170             IPSGD1=K(IPSPD,4)
2171             IF(IM.EQ.NS+2) THEN
2172               PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
2173             ELSE
2174               PHI=PYANGL(-P(IPSGD1,1),P(IPSGD1,2))
2175             ENDIF
2176           ENDIF
2177         ELSEIF(MPSPD.EQ.1.AND.IGM.EQ.NS+2) THEN
2178           IPSPD=IP1+IM-NS-2
2179           IF(K(IPSPD,4).GT.0) THEN
2180             IPSGD1=K(IPSPD,4)
2181             PHIPSM=PYANGL(P(IPSPD,1),P(IPSPD,2))
2182             THEPSM=PYANGL(P(IPSPD,3),SQRT(P(IPSPD,1)**2+P(IPSPD,2)**2))
2183             CALL PYROBO(IPSGD1,IPSGD1,0D0,-PHIPSM,0D0,0D0,0D0)
2184             CALL PYROBO(IPSGD1,IPSGD1,-THEPSM,0D0,0D0,0D0,0D0)
2185             PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
2186             CALL PYROBO(IPSGD1,IPSGD1,THEPSM,PHIPSM,0D0,0D0,0D0)
2187           ENDIF
2188         ENDIF
2189
2190 C...Construct momenta for ordinary branching in shower.
2191         P(N+1,1)=PT*COS(PHI)
2192         P(N+1,2)=PT*SIN(PHI)
2193         IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
2194      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
2195           P(N+1,3)=0.5D0*(PZM*(V(IM,5)+V(N+1,5)-V(N+2,5))+
2196      &    PEM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5)
2197         ELSEIF(PZM.GT.0D0) THEN
2198           P(N+1,3)=0.5D0*(V(N+2,5)-V(N+1,5)-V(IM,5)+
2199      &    2D0*PEM*P(N+1,4))/PZM
2200         ELSE
2201           P(N+1,3)=0D0
2202         ENDIF
2203         P(N+2,1)=-P(N+1,1)
2204         P(N+2,2)=-P(N+1,2)
2205         P(N+2,3)=PZM-P(N+1,3)
2206         P(N+2,4)=PEM-P(N+1,4)
2207         IF(MSTJ(43).LE.2) THEN
2208           V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
2209           V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
2210         ENDIF
2211       ENDIF
2212
2213 C...Rotate and boost daughters.
2214       IF(IGM.GT.0) THEN
2215         IF(MSTJ(43).LE.2) THEN
2216           BEX=P(IGM,1)/P(IGM,4)
2217           BEY=P(IGM,2)/P(IGM,4)
2218           BEZ=P(IGM,3)/P(IGM,4)
2219           GA=P(IGM,4)/P(IGM,5)
2220           GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1D0+GA)-
2221      &    P(IM,4))
2222         ELSE
2223           BEX=0D0
2224           BEY=0D0
2225           BEZ=0D0
2226           GA=1D0
2227           GABEP=0D0
2228         ENDIF
2229         PTIMB=SQRT((P(IM,1)+GABEP*BEX)**2+(P(IM,2)+GABEP*BEY)**2)
2230         THE=PYANGL(P(IM,3)+GABEP*BEZ,PTIMB)
2231         IF(PTIMB.GT.1D-4) THEN
2232           PHI=PYANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
2233         ELSE
2234           PHI=0D0
2235         ENDIF
2236         DO 570 I=N+1,N+2
2237           DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
2238      &    SIN(THE)*COS(PHI)*P(I,3)
2239           DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
2240      &    SIN(THE)*SIN(PHI)*P(I,3)
2241           DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
2242           DP(4)=P(I,4)
2243           DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
2244           DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
2245           P(I,1)=DP(1)+DGABP*BEX
2246           P(I,2)=DP(2)+DGABP*BEY
2247           P(I,3)=DP(3)+DGABP*BEZ
2248           P(I,4)=GA*(DP(4)+DBP)
2249   570   CONTINUE
2250       ENDIF
2251
2252 C...Weight with azimuthal distribution, if required.
2253       IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
2254         DO 580 J=1,3
2255           DPT(1,J)=P(IM,J)
2256           DPT(2,J)=P(IAU,J)
2257           DPT(3,J)=P(N+1,J)
2258   580   CONTINUE
2259         DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
2260         DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
2261         DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
2262         DO 590 J=1,3
2263           DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/MAX(1D-10,DPMM)
2264           DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/MAX(1D-10,DPMM)
2265   590   CONTINUE
2266         DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
2267         DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
2268         IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1D0*PARJ(82)) THEN
2269           CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
2270      &    DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
2271           IF(MAZIP.NE.0) THEN
2272             IF(1D0+HAZIP*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(HAZIP)))
2273      &      GOTO 560
2274           ENDIF
2275           IF(MAZIC.NE.0) THEN
2276             IF(MAZIC.EQ.N+2) CAD=-CAD
2277             IF((1D0-HAZIC)*(1D0-HAZIC*CAD)/(1D0+HAZIC**2-2D0*HAZIC*CAD)
2278      &      .LT.PYR(0)) GOTO 560
2279           ENDIF
2280         ENDIF
2281       ENDIF
2282
2283 C...Azimuthal anisotropy due to interference with initial state partons.
2284       IF(MOD(MIIS,2).EQ.1.AND.IGM.EQ.NS+1.AND.(K(N+1,2).EQ.21.OR.
2285      &K(N+2,2).EQ.21)) THEN
2286         III=IM-NS-1
2287         IF(ISII(III).GE.1) THEN
2288           IAZIID=N+1
2289           IF(K(N+1,2).NE.21) IAZIID=N+2
2290           IF(K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
2291      &    P(N+1,4).GT.P(N+2,4)) IAZIID=N+2
2292           THEIID=PYANGL(P(IAZIID,3),SQRT(P(IAZIID,1)**2+P(IAZIID,2)**2))
2293           IF(III.EQ.2) THEIID=PARU(1)-THEIID
2294           PHIIID=PYANGL(P(IAZIID,1),P(IAZIID,2))
2295           HAZII=MIN(0.95D0,THEIID/THEIIS(III,ISII(III)))
2296           CAD=COS(PHIIID-PHIIIS(III,ISII(III)))
2297           PHIREL=ABS(PHIIID-PHIIIS(III,ISII(III)))
2298           IF(PHIREL.GT.PARU(1)) PHIREL=PARU(2)-PHIREL
2299           IF((1D0-HAZII)*(1D0-HAZII*CAD)/(1D0+HAZII**2-2D0*HAZII*CAD)
2300      &    .LT.PYR(0)) GOTO 560
2301         ENDIF
2302       ENDIF
2303
2304 C...Continue loop over partons that may branch, until none left.
2305       IF(IGM.GE.0) K(IM,1)=14
2306       N=N+NEP
2307       NEP=2
2308       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
2309         CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
2310         IF(MSTU(21).GE.1) N=NS
2311         IF(MSTU(21).GE.1) RETURN
2312       ENDIF
2313       GOTO 290
2314
2315 C...Set information on imagined shower initiator.
2316   600 IF(NPA.GE.2) THEN
2317         K(NS+1,1)=11
2318         K(NS+1,2)=94
2319         K(NS+1,3)=IP1
2320         IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
2321         K(NS+1,4)=NS+2
2322         K(NS+1,5)=NS+1+NPA
2323         IIM=1
2324       ELSE
2325         IIM=0
2326       ENDIF
2327
2328 C...Reconstruct string drawing information.
2329       DO 610 I=NS+1+IIM,N
2330         KQ=KCHG(PYCOMP(K(I,2)),2)
2331         IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
2332           K(I,1)=1
2333         ELSEIF(K(I,1).LE.10.AND.IABS(K(I,2)).GE.11.AND.
2334      &    IABS(K(I,2)).LE.18) THEN
2335           K(I,1)=1
2336         ELSEIF(K(I,1).LE.10) THEN
2337           K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
2338           K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
2339         ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
2340           ID1=MOD(K(I,4),MSTU(5))
2341           IF(KQ.EQ.1.AND.K(I,2).GT.0) ID1=MOD(K(I,4),MSTU(5))+1
2342           IF(KQ.EQ.2.AND.(K(ID1,2).EQ.21.OR.K(ID1+1,2).EQ.21).AND.
2343      &    PYR(0).GT.0.5D0) ID1=MOD(K(I,4),MSTU(5))+1
2344           ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
2345           K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
2346           K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
2347           K(ID1,4)=K(ID1,4)+MSTU(5)*I
2348           K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
2349           K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
2350           K(ID2,5)=K(ID2,5)+MSTU(5)*I
2351         ELSE
2352           ID1=MOD(K(I,4),MSTU(5))
2353           ID2=ID1+1
2354           K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
2355           K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
2356           IF(KQ.EQ.1.OR.K(ID1,1).GE.11) THEN
2357             K(ID1,4)=K(ID1,4)+MSTU(5)*I
2358             K(ID1,5)=K(ID1,5)+MSTU(5)*I
2359           ELSE
2360             K(ID1,4)=0
2361             K(ID1,5)=0
2362           ENDIF
2363           K(ID2,4)=0
2364           K(ID2,5)=0
2365         ENDIF
2366   610 CONTINUE
2367
2368 C...Transformation from CM frame.
2369       IF(NPA.EQ.1) THEN
2370         THE=PYANGL(P(IPA(1),3),SQRT(P(IPA(1),1)**2+P(IPA(1),2)**2))
2371         PHI=PYANGL(P(IPA(1),1),P(IPA(1),2))
2372         MSTU(33)=1
2373         CALL PYROBO(NS+1,N,THE,PHI,0D0,0D0,0D0)
2374       ELSEIF(NPA.EQ.2) THEN
2375         BEX=PS(1)/PS(4)
2376         BEY=PS(2)/PS(4)
2377         BEZ=PS(3)/PS(4)
2378         GA=PS(4)/PS(5)
2379         GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
2380      &  /(1D0+GA)-P(IPA(1),4))
2381         THE=PYANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
2382      &  +GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
2383         PHI=PYANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
2384         MSTU(33)=1
2385         CALL PYROBO(NS+1,N,THE,PHI,BEX,BEY,BEZ)
2386       ELSE
2387         CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),
2388      &  PS(3)/PS(4))
2389         MSTU(33)=1
2390         CALL PYROBO(NS+1,N,0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),PS(3)/PS(4))
2391       ENDIF
2392
2393 C...Decay vertex of shower.
2394       DO 630 I=NS+1,N
2395         DO 620 J=1,5
2396           V(I,J)=V(IP1,J)
2397   620   CONTINUE
2398   630 CONTINUE
2399
2400 C...Delete trivial shower, else connect initiators.
2401       IF(N.LE.NS+NPA+IIM) THEN
2402         N=NS
2403       ELSE
2404         DO 640 IP=1,NPA
2405           K(IPA(IP),1)=14
2406           K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP
2407           K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP
2408           K(NS+IIM+IP,3)=IPA(IP)
2409           IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1
2410           IF(K(NS+IIM+IP,1).NE.1) THEN
2411             K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)
2412             K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)
2413           ENDIF
2414   640   CONTINUE
2415       ENDIF
2416
2417       RETURN
2418       END
2419 C
2420       SUBROUTINE QPYGIN(X0,Y0,Z0,T0)
2421 C     USER-DEFINED ROUTINE: IT SETS THE INITIAL POSITION AND TIME OF THE
2422 C     PARENT BRANCHING PARTON (X, Y, Z, T, IN FM) IN THE CENTER-OF-MASS
2423 C     FRAME OF THE HARD COLLISION (IF APPLICABLE FOR THE TYPE OF EVENTS
2424 C     YOU ARE SIMULATING). INFORMATION ABOUT THE BOOST AND ROTATION IS
2425 C     CONTAINED IN THE IN COMMON QPLT BELOW.
2426       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2427 C     NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST
2428 C     PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE
2429 C     CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD
2430 C     SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
2431       COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ
2432 cforalice+
2433 c     Here the transverse coordinates of the hard scattering are set by
2434 c     glauber geometry.
2435       call GetRandomXY(xrang,yrang)
2436       xin=xrang ! fm
2437       yin=yrang ! fm
2438 cforalice-
2439       zin=0.d0 ! fm
2440       tin=0.d0 ! fm
2441       call qpyrobo(xin,yin,zin,tin,0.d0,0.d0,bbx,bby,bbz,
2442      + x1,y1,z1,t1)
2443       call qpyrobo(x1,y1,z1,t1,0.d0,aa2,0.d0,0.d0,0.d0,
2444      + x2,y2,z2,t2)
2445       call qpyrobo(x2,y2,z2,t2,aa1,0.d0,0.d0,0.d0,0.d0,
2446      + xout,yout,zout,tout)
2447       x0=xout
2448       y0=yout
2449       z0=zout
2450       t0=tout
2451       RETURN
2452       END
2453 C
2454       SUBROUTINE QPYGEO(x,y,z,t,bx,by,bz,qhl,oc)
2455 C     USER-DEFINED ROUTINE:
2456 C     The values of qhatL and omegac have to be computed
2457 C     by the user, using his preferred medium model, in
2458 C     this routine, which takes as input the position
2459 C     x,y,z,t (in fm) of the parton to branch, the trajectory
2460 C     defined by the three-vector bx,by,bz (in units of c),
2461 C     (all values in the center-of-mass frame of the
2462 C     hard collision), and returns the value of qhatL
2463 C     (in GeV**2) and omegac (in GeV).
2464       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2465 C     NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST
2466 C     PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE
2467 C     CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD
2468 C     SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
2469       COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ
2470       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2471       qhat=parj(198)
2472       xl=parj(199)
2473       bimp=parj(197)
2474 c     Here we give five options for the geometry of the medium:
2475
2476 c$$$cforalice+
2477 c$$$     (0) we call routine CalculateI0I1 in AliFastGlauber. Given the position
2478 c$$$     of the parton in the reaction plane (x,y), the direction in the
2479 c$$$     reaction plane phi=atan(py/px) and the impact parameter of the
2480 c$$$     collision bimp, it gives back the transverse path length to the
2481 c$$$     "end" of the medium and the integrated qhat along that path length.
2482 c$$$     See the seter for this option in Configqpythia.C(one can set the
2483 c$$$     value of xkscale by doing SetQhat(xkscale), with xkscale in fm).
2484 c$$$     The set value is passed here through the pythia free parameter parj(198).
2485
2486
2487       xkscale=parj(198)
2488       ellcut=20.d0
2489       xlcero=0.d0
2490       xlone=0.d0
2491
2492       call qpyrobo(x,y,z,t,-aa1,-aa2,-bbx,-bby,-bbz,xout,
2493      + yout,zout,tout)
2494
2495       call qpyrobo(bx,by,bz,1.d0,-aa1,-aa2,-bbx,-bby,-bbz,
2496      + bx1,by1,bz1,bt1)
2497
2498
2499       phi=datan2(by1,bx1)
2500       phia=phi
2501       bimpa=bimp
2502       ellcuta=ellcut
2503       xa=xout
2504       ya=yout
2505       call CalculateI0I1(xlcero,xlone,bimpa,xa,ya,phia,ellcuta)
2506       if(xlcero.eq.0.d0) then
2507            xlp=0.d0
2508            qhl=0.d0
2509       else
2510       xlp=2.d0*xlone/xlcero
2511       qhl=0.1973d0*0.1973d0*xlcero*xkscale
2512       endif
2513
2514
2515 c$$$cforalice-
2516 c     To use any of these folowing 1,2,3 or 4 options the user should specify
2517 c     a constant value for the transport coefficient and an initial in medium length.
2518 c     This can be done in the user Config file by setting: SetQhat(qhat), with qhat in
2519 c     GeV2/fm and SetLength(xl), with xl in fm. Those values are passed here through
2520 c     pythia free parameters parj(198) and parj(199).
2521
2522 c     (1) to fix the length to the initial value, uncomment the next three lines
2523 c     and comment the other definitions of xlp and qhl above and below.
2524 c      xlp=xl
2525 c      if (xlp .lt. 0.d0) xlp=0.d0
2526 c      qhl=xlp*qhat ! GeV**2
2527 c        print*, xlp,qhl
2528 c     (2) simplest ansatz: for an initial parton along the z-axis (approximate)
2529 c      starting in the center of a medium (-xl,+xl) along the z-axis
2530 c       if (bz .gt. 0.d0) then
2531 c         xlp=xl-z
2532 c       else
2533 c         xlp=xl+z
2534 c       endif
2535 c      if (xlp .gt. (2.d0*xl)) xlp=2.d0*xl
2536 c      if (xlp .lt. 0.d0) xlp=0.d0
2537 c      qhl=xlp*qhat ! GeV**2
2538
2539 c     (3) for a parton at midrapidity inside a cylinder (approximate)
2540 c      xlp=xl-dsqrt(x*x+y*y)
2541 c      if (xlp .lt. 0.d0) xlp=0.d0
2542 c      qhl=xlp*qhat ! GeV**2
2543
2544 c     (4) for a brick defined by planes (x,y,0) and (x,y,xl), comment
2545 c     the previous lines and uncomment lines between the comment 'brick'.
2546 c     brick+
2547 c       if (z .ge. 0.d0 .and. z .le. xl)
2548 c     >    then
2549 c            if (bz .gt. 0.d0) then
2550 c               ttpp=(xl-z)/bz
2551 c               xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2552 c     >             (xl-z)**2.d0)
2553 c            else
2554 c               ttpp=z/dabs(bz)
2555 c               xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2556 c     >             (z)**2.d0)
2557 c            endif
2558 c         elseif (z .lt. 0.d0) then
2559 c           if (bz .lt. 0.d0) then
2560 c              xlp=0.d0
2561 c           else
2562 c              ttpp1=-z/bz
2563 c              ttpp2=(xl-z)/bz
2564 c              xxpp1=x+bx*ttpp1
2565 c              xxpp2=x+bx*ttpp2
2566 c              yypp1=y+by*ttpp1
2567 c              yypp2=y+by*ttpp2
2568 c              xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2569 c     >                  xl**2.d0)
2570 c           endif
2571 c         elseif (z .gt. xl) then
2572 c           if (bz .gt. 0.d0) then
2573 c              xlp=0.d0
2574 c           else
2575 c              ttpp1=z/dabs(bz)
2576 c              ttpp2=(-xl+z)/dabs(bz)
2577 c              xxpp1=x+bx*ttpp1
2578 c              xxpp2=x+bx*ttpp2
2579 c              yypp1=y+by*ttpp1
2580 c              yypp2=y+by*ttpp2
2581 c              xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2582 c     >                  xl**2.d0)
2583 c           endif
2584 c         endif
2585 c      if (xlp .lt. 0.d0) xlp=0.d0
2586 c      qhl=xlp*qhat ! GeV**2
2587 c     brick-
2588
2589       oc=0.5d0*qhl*xlp/0.1973d0 ! GeV
2590       RETURN
2591       END