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