]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/QPYTHIA/q-pyshow.1.0.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added static method to retrieve a QA data object from a given list
[u/mrichter/AliRoot.git] / PYTHIA6 / QPYTHIA / q-pyshow.1.0.F
CommitLineData
b527e4b2 1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2C
3C AUXILIARY ROUTINES FOR Q-PYTHIA version 1.0.
4C
5C DATE: 26.09.2008.
6C
7C AUTHORS: N. Armesto, L. Cunqueiro and C. A. Salgado
8C Departamento de Fisica de Particulas and IGFAE
9C Universidade de Santiago de Compostela
10C 15706 Santiago de Compostela, Spain
11C
12C EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es,
13C Carlos.Salgado@cern.ch
14C
15C CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18.
16C NOT to be modified by user.
17C
18C WHEN USING Q-PYTHIA, PLEASE QUOTE:
19C
20C 1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado,
21C in preparation.
22C 2) T. Sjostrand, S. Mrenna and P. Skands,
23C ``PYTHIA 6.4 physics and manual,''
24C JHEP 0605 (2006) 026 [arXiv:hep-ph/0603175].
25C
26C DISCLAIMER: this program comes without any guarantees. Beware of
27C errors and use common sense when interpreting results.
28C Any modifications are done under exclusive
29C maker's resposibility.
30C
31CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
32C
33c VERSION WITH THE LARGE X BEHAVIOR OF THE MEDIUM PART INTRODUCED
34C BY MULTIPLYING BY THE NUMERATOR OF THE COLLINEAR PART OF THE
35C VACUUM SPLITTING FUNCTION
36c
37 function splitq1(w)
38c 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
50c
51 function splitg1(w)
52c 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)
56c 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
65c
66 function splitq2(w)
67c 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
80c
81 function splitg2(z)
82c 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
94c
95 function splitgq(z)
96c 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
102c
103 function splitgg(z)
104c 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
112c
113 function splitqqbar(z)
114c 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
121c
122 function splitmedg1(z)
123c 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
132c 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)
144cc 1/2 to avoid double counting
145c splitmedg=0.5d0*xnc*2.d0*pi*zz*t*fff/qhatl
146c we multiply by max(z,1-z) to introduce the large z behavior from the
147c numerator in the vacuum
148 flx=max(zz,auxz)
149c 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
154c
155 function splitmedq1(z)
156c 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)
171c splitmedq=(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
172c we multiply by 1+z**2 to introduce the large z behavior from the
173c 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
179c
180 function splitmedg2(z)
181c 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
190c 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)
202cc 1/2 to avoid double counting
203c splitmedg=0.5d0*xnc*2.d0*pi*zz*t*fff/qhatl
204c we multiply by max(z,1-z) to introduce the large z behavior from the
205c numerator in the vacuum
206 flx=max(zz,auxz)
207c 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
212c
213 function splitmedq2(z)
214c 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)
229c splitmedq=(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
230c we multiply by 1+z**2 to introduce the large z behavior from the
231c 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
237c
238 function genspec(ome,xk2)
239C THIS FUNCTION GENERATES (omega/omegac) dI/d(omega/omegac) dkappa2,
240C omegac=qhat L/2, kappa2=kt2/(qhat L), in the mss approximation for m=0,
241c using interpolation and extrapolation. It reads file grid-qp.dat.
242c ome=omega/omegac, xk2=kappa2.
243C MAXIMUM GRID 101 TIMES 101, MODIFY ARRAY DIMENSIONS IF EXCEEDED.
244c 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)
b527e4b2 249 save xkap2, xlkap2, xome, xlome, xspec, npkap, npome
6baf84ad 250 character*1000 filnam
251 character*1000 chroot
b527e4b2 252 DATA IFLAG/0/
253c WE READ THE GRID ONLY THE FIRST TIME.
254 IF (IFLAG .EQ. 0) THEN
6baf84ad 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')
b527e4b2 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)
26710 continue
268 do 20 i=1, npome, 1
269 read(11,*) xome(i), xlome(i)
27020 continue
271 do 30 j=1, npome, 1
272 do 40 i=1, npkap, 1
273 read(11,*) xspec(i,j)
27440 continue
27530 continue
276 close(11)
277 iflag=1
278 ENDIF
279c cases
280c for ome>largest value set to 0,
281c for xk2< smallest value frozen,
282c 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)
29850 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
31660 continue
31770 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)
33080 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)
34990 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
360c
361 RETURN
362 END
363C
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)
373c 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/
379C
380C TABULAR INTERPOLATION USING SYMMETRICALLY PLACED ARGUMENT POINTS.
381C
382C 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
390C (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
395C (IF TRUE.)
396 2 IX=MID
397 3 IF(IY-IX.GT.1) GO TO 1
398 GO TO 7
399C (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
404C (IF TRUE.)
405 5 IX=MID
406 6 IF(IY-IX.GT.1) GO TO 4
407C
408C COPY REORDERED INTERPOLATION POINTS INTO (T(I),D(I)), SETTING
409C *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
418C (SKIP POINT.)
419 NPTS=MPLUS
420 GO TO 11
421C (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
427C
428C REPLACE D BY THE LEADING DIAGONAL OF A DIVIDED-DIFFERENCE TABLE, SUP-
429C 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
441C
442C EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES
443C 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
453C
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
470c
471C COPY OF CERN DGAUSS
472C
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
535c
536 FUNCTION SIMDIS(Numb,zmin,nzur,RI)
537C IT SIMULATES A RANDOM NUMBER ACCORDING TO A DISCRETE DISTRIBUTION GIVEN
538C BY ARRAY YA AT POINTS XA. THOUGHT FOR PYTHIA (PYR(0)).
539C N: NUMBER OF POINTS IN THE ARRAYS.
540C XA: ARRAY OF X-VALUES.
541C YA: ARRAY OF Y-VALUES.
542c 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
cac2b066 574
575
576 SUBROUTINE QPYROBO(XI,YI,ZI,TI,THE,PHI,BEX,BEY,BEZ,XP,YP,ZP,TP)
577C N. Armesto, 16.04.2009
578C performs a boost and rotation of (t,x,y,z) to (tp,xp,yp,zp):
579C cut version of PYROBO, angles and boost parameters identical.
580 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
581C...Local arrays.
582 DIMENSION ROT(3,3),VR(3),DV(4)
583C
584 X=XI
585 Y=YI
586 Z=ZI
587 T=TI
588C...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)
599C Instead of loop 120 in PYROBO.
600 VR(1)=X
601 VR(2)=Y
602 VR(3)=Z
603C 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
611C If nothing happens...
612 XP=X
613 YP=Y
614 ZP=Z
615 TP=T
616C...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
624C...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)
632C 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
644C
645 RETURN
646 END
647
648
649
650
651
652
653
654
655
b527e4b2 656CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
657C
658C PYSHOW ROUTINE FOR Q-PYTHIA version 1.0.
659C
660C DATE: 26.09.2008.
661C
662C AUTHORS: N. Armesto, L. Cunqueiro and C. A. Salgado
663C Departamento de Fisica de Particulas and IGFAE
664C Universidade de Santiago de Compostela
665C 15706 Santiago de Compostela, Spain
666C
667C EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es,
668C Carlos.Salgado@cern.ch
669C
670C CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18.
671C
672C WHEN USING Q-PYTHIA, PLEASE QUOTE:
673C
674C 1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado,
675C in preparation.
676C 2) T. Sjostrand, S. Mrenna and P. Skands,
677C ``PYTHIA 6.4 physics and manual,''
678C JHEP 0605 (2006) 026 [arXiv:hep-ph/0603175].
679C
680C INSTRUCTIONS: initial parton position is initialized by a call
681C to user-defined routine qpygin(x0,y0,z0,t0),
682C where these are the initial coordinates in the
683C center-of-mass frame of the hard collision
684C (if applicable for the type of process you study).
685C The values of qhatL and omegac have to be computed
686C by the user, using his preferred medium model, in
687C routine qpygeo, which takes as input the position
688C x,y,z,t of the parton to branch, the trajectory
689C defined by the three-vector betax,betay,betaz,
690C (all values in the center-of-mass frame of the
691C hard collision), and returns the value of qhatL
692C (in GeV**2) and omegac (in GeV).
693C Both routines are to be found at the end of this file.
694C
695C DISCLAIMER: this program comes without any guarantees. Beware of
696C errors and use common sense when interpreting results.
697C Any modifications are done under exclusive
698C maker's resposibility.
699C
700CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
701C*********************************************************************
702
703C...PYSHOW
704C...Generates timelike parton showers from given partons.
705
706 SUBROUTINE PYSHOWQ(IP1,IP2,QMAX)
707
708C...Double precision and integer declarations.
709 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
710 IMPLICIT INTEGER(I-N)
711 INTEGER PYK,PYCHGE,PYCOMP
712C...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)
716Cacs+
717 PARAMETER (NNPOS=4000)
718 DIMENSION PPOS(NNPOS,4)
719Cacs-
720C...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/
728Cacs+
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/
739Cacs-
740C...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)
746Cacs+
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
790Cacs-
791
792C...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
802C...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
830C...Send off to PYPTFS for pT-ordered evolution if requested,
831C...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
843C...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
890C...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)
900C...Special cutoff masses for initial partons (may be a heavy quark,
901C...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
921C...QUARKONIA+++
922C...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
926C...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
958C...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
979C...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
999C...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
1007C...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
1021C...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
1029C...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
1056C...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
1063C...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
1075C...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
1088C...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
1099C...~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
1109C...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
1117C...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
1143C...Boost interfering initial partons to rest frame
1144C...and reconstruct their polar and azimuthal angles.
1145Cacs+
1146 qplta1=0.d0
1147 qplta2=0.d0
1148 qpltbx=0.d0
1149 qpltby=0.d0
1150 qpltbz=0.d0
1151Cacs-
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)
1180Cacs+
1181 qplta1=-the
1182 qplta2=-phi
1183 qpltbx=-PS(1)/PS(4)
1184 qpltby=-PS(2)/PS(4)
1185 qpltbz=-PS(3)/PS(4)
1186Cacs-
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
1198C...Boost 3 or more partons to their rest frame.
1199Cacs+
1200c IF(NPA.GE.3) CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
1201c &-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)
cac2b066 1210
b527e4b2 1211 ENDIF
1212Cacs-
1213
1214C...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
1240Cacs+
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
124710101 continue
1248
1249Cacs-
1250
1251
1252
1253
1254C...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
1274C...Position of aunt (sister to branching parton).
1275C...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
1310C...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
1324C...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
1349C...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
1375C...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
1385C...Store information on choice of evolving daughter.
1386 IEP(1)=N+INUM
1387Cacs+
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)
cac2b066 1417
1418
1419
b527e4b2 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
1427Cacs-
cac2b066 1428
b527e4b2 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
1446C...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
1459C...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
1473C...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
1503C...Integral of Altarelli-Parisi z kernel for QCD.
1504C...(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
1506Cacs+
1507C FBR= 6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
1508 FBR=dgauss1(splitg1,zc,1.d0-zc,1.d-3)
1509Cacs-
1510C...QUARKONIA+++
1511C...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)
1515C...QUARKONIA---
1516 ELSEIF(MSTJ(49).EQ.0) THEN
1517Cacs+
1518C FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC)
1519 FBR=dgauss1(splitq1,zc,1.d0-zc,1.d-3)
1520
1521Cacs-
1522 IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0)
1523
1524C...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
1531C...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
1538C...Reset QCD probability for colourless.
1539 IF(ISCOL(IR).EQ.0) FBR=0D0
1540
1541C...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
1550C...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
1563C...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
1568C...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)
1570C...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
1582C...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
1588C...Select mass for daughter in QED evolution.
1589 IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1.AND.IPSPD.EQ.0) THEN
1590C...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
1598C...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
1608C...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
1615Cacs+
1616 virt=V(IEP(1),5)
1617Cacs-
1618
1619C...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
1636C...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
1642C...QUARKONIA+++
1643C...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)
1647C...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
1651C...QUARKONIA---
1652
1653C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
1654 ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
1655Cacs+
1656C Z=1D0-(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1657C...Only do z weighting when no ME correction afterwards.
1658C IF(M3JC.EQ.0.AND.1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
1659C
1660 anfbr=dgauss1(splitq2,zc,1.d0-zc,1.d-3)
1661 z=simdis(500,zc,1,anfbr)
1662Cacs-
1663 K(IEP(1),5)=21
1664 ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN
1665Cacs+
1666c 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
1672c IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
1673Cacs-
1674 K(IEP(1),5)=21
1675 ELSEIF(MSTJ(49).NE.1) THEN
1676Cacs+
1677c Z=PYR(0)
1678c 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
1682Cacs-
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
1696C...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
1711C...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
1725C...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
1768C...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
1783C...Three-jet matrix element correction.
1784 IF(M3JC.GE.1) THEN
1785 WME=1D0
1786 WSHOW=1D0
1787
1788C...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
1802C...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
1808C...QCD ME: original parton, first branching.
1809 PM2ME=PMTH(1,63-IR)
1810 ECMME=PS(5)
1811 ELSEIF(IR.GE.31) THEN
1812C...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
1817C...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
1828C...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
1834C...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
1839C...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
1845C...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)
1854C...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
1860C...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
1874C...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
1907C...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
1921C...Impose angular constraint in first branching from interference
1922C...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
1932C...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
1945C...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
2034C...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
2057C...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
2083C...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
2109C...Global statistics.
2110 MINT(353)=MINT(353)+1
2111 VINT(353)=VINT(353)+PT
2112 IF (MINT(353).EQ.1) VINT(358)=PT
2113
2114C...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
2134C...Find coefficient of azimuthal asymmetry due to soft gluon
2135C...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
2155C...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
2168C...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
2194C...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
2217C...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
2256C...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
2287C...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
2308C...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
2319C...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
2332C...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
2372C...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
2397C...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
2404C...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
2423C
2424 SUBROUTINE QPYGIN(X0,Y0,Z0,T0)
2425C USER-DEFINED ROUTINE: IT SETS THE INITIAL POSITION AND TIME OF THE
2426C PARENT BRANCHING PARTON (X, Y, Z, T, IN FM) IN THE CENTER-OF-MASS
2427C FRAME OF THE HARD COLLISION (IF APPLICABLE FOR THE TYPE OF EVENTS
2428C YOU ARE SIMULATING). INFORMATION ABOUT THE BOOST AND ROTATION IS
2429C CONTAINED IN THE IN COMMON QPLT BELOW.
2430 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2431C NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST
2432C PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE
2433C CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD
2434C SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
2435 COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ
2436cforalice+
1bab4b79 2437 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
99994a98 2438 SAVE /PYPARS/
b527e4b2 2439c Here the transverse coordinates of the hard scattering are set by
2440c glauber geometry.
2441 call GetRandomXY(xrang,yrang)
cac2b066 2442 xin=xrang ! fm
2443 yin=yrang ! fm
e9719084 2444 call savexy(xin, yin)
b527e4b2 2445cforalice-
cac2b066 2446 zin=0.d0 ! fm
2447 tin=0.d0 ! fm
2448 call qpyrobo(xin,yin,zin,tin,0.d0,0.d0,bbx,bby,bbz,
2449 + x1,y1,z1,t1)
2450 call qpyrobo(x1,y1,z1,t1,0.d0,aa2,0.d0,0.d0,0.d0,
2451 + x2,y2,z2,t2)
2452 call qpyrobo(x2,y2,z2,t2,aa1,0.d0,0.d0,0.d0,0.d0,
2453 + xout,yout,zout,tout)
2454 x0=xout
2455 y0=yout
2456 z0=zout
2457 t0=tout
b527e4b2 2458 RETURN
2459 END
2460C
2461 SUBROUTINE QPYGEO(x,y,z,t,bx,by,bz,qhl,oc)
2462C USER-DEFINED ROUTINE:
2463C The values of qhatL and omegac have to be computed
2464C by the user, using his preferred medium model, in
2465C this routine, which takes as input the position
2466C x,y,z,t (in fm) of the parton to branch, the trajectory
2467C defined by the three-vector bx,by,bz (in units of c),
2468C (all values in the center-of-mass frame of the
2469C hard collision), and returns the value of qhatL
2470C (in GeV**2) and omegac (in GeV).
2471 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2472C NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST
2473C PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE
2474C CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD
2475C SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
2476 COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ
2477 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2478 qhat=parj(198)
2479 xl=parj(199)
2480 bimp=parj(197)
2481c Here we give five options for the geometry of the medium:
2482
cac2b066 2483c$$$cforalice+
2484c$$$ (0) we call routine CalculateI0I1 in AliFastGlauber. Given the position
2485c$$$ of the parton in the reaction plane (x,y), the direction in the
2486c$$$ reaction plane phi=atan(py/px) and the impact parameter of the
2487c$$$ collision bimp, it gives back the transverse path length to the
2488c$$$ "end" of the medium and the integrated qhat along that path length.
2489c$$$ See the seter for this option in Configqpythia.C(one can set the
2490c$$$ value of xkscale by doing SetQhat(xkscale), with xkscale in fm).
2491c$$$ The set value is passed here through the pythia free parameter parj(198).
2492
b527e4b2 2493
2494 xkscale=parj(198)
2495 ellcut=20.d0
2496 xlcero=0.d0
2497 xlone=0.d0
cac2b066 2498
2499 call qpyrobo(x,y,z,t,-aa1,-aa2,-bbx,-bby,-bbz,xout,
2500 + yout,zout,tout)
2501
2502 call qpyrobo(bx,by,bz,1.d0,-aa1,-aa2,-bbx,-bby,-bbz,
2503 + bx1,by1,bz1,bt1)
2504
b527e4b2 2505
b527e4b2 2506 phi=datan2(by1,bx1)
b527e4b2 2507 phia=phi
b527e4b2 2508 bimpa=bimp
2509 ellcuta=ellcut
cac2b066 2510 xa=xout
2511 ya=yout
b527e4b2 2512 call CalculateI0I1(xlcero,xlone,bimpa,xa,ya,phia,ellcuta)
2513 if(xlcero.eq.0.d0) then
2514 xlp=0.d0
2515 qhl=0.d0
2516 else
2517 xlp=2.d0*xlone/xlcero
cac2b066 2518 qhl=0.1973d0*0.1973d0*xlcero*xkscale
e9719084 2519 call savei0i1(xlcero, xlone)
cac2b066 2520 endif
2521
2522
2523c$$$cforalice-
b527e4b2 2524c To use any of these folowing 1,2,3 or 4 options the user should specify
2525c a constant value for the transport coefficient and an initial in medium length.
2526c This can be done in the user Config file by setting: SetQhat(qhat), with qhat in
2527c GeV2/fm and SetLength(xl), with xl in fm. Those values are passed here through
2528c pythia free parameters parj(198) and parj(199).
2529
2530c (1) to fix the length to the initial value, uncomment the next three lines
2531c and comment the other definitions of xlp and qhl above and below.
2532c xlp=xl
2533c if (xlp .lt. 0.d0) xlp=0.d0
2534c qhl=xlp*qhat ! GeV**2
cac2b066 2535c print*, xlp,qhl
b527e4b2 2536c (2) simplest ansatz: for an initial parton along the z-axis (approximate)
2537c starting in the center of a medium (-xl,+xl) along the z-axis
2538c if (bz .gt. 0.d0) then
2539c xlp=xl-z
2540c else
2541c xlp=xl+z
2542c endif
2543c if (xlp .gt. (2.d0*xl)) xlp=2.d0*xl
2544c if (xlp .lt. 0.d0) xlp=0.d0
2545c qhl=xlp*qhat ! GeV**2
2546
2547c (3) for a parton at midrapidity inside a cylinder (approximate)
2548c xlp=xl-dsqrt(x*x+y*y)
2549c if (xlp .lt. 0.d0) xlp=0.d0
2550c qhl=xlp*qhat ! GeV**2
2551
2552c (4) for a brick defined by planes (x,y,0) and (x,y,xl), comment
2553c the previous lines and uncomment lines between the comment 'brick'.
2554c brick+
2555c if (z .ge. 0.d0 .and. z .le. xl)
2556c > then
2557c if (bz .gt. 0.d0) then
2558c ttpp=(xl-z)/bz
2559c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2560c > (xl-z)**2.d0)
2561c else
2562c ttpp=z/dabs(bz)
2563c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2564c > (z)**2.d0)
2565c endif
2566c elseif (z .lt. 0.d0) then
2567c if (bz .lt. 0.d0) then
2568c xlp=0.d0
2569c else
2570c ttpp1=-z/bz
2571c ttpp2=(xl-z)/bz
2572c xxpp1=x+bx*ttpp1
2573c xxpp2=x+bx*ttpp2
2574c yypp1=y+by*ttpp1
2575c yypp2=y+by*ttpp2
2576c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2577c > xl**2.d0)
2578c endif
2579c elseif (z .gt. xl) then
2580c if (bz .gt. 0.d0) then
2581c xlp=0.d0
2582c else
2583c ttpp1=z/dabs(bz)
2584c ttpp2=(-xl+z)/dabs(bz)
2585c xxpp1=x+bx*ttpp1
2586c xxpp2=x+bx*ttpp2
2587c yypp1=y+by*ttpp1
2588c yypp2=y+by*ttpp2
2589c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2590c > xl**2.d0)
2591c endif
2592c endif
2593c if (xlp .lt. 0.d0) xlp=0.d0
2594c qhl=xlp*qhat ! GeV**2
2595c brick-
2596
2597 oc=0.5d0*qhl*xlp/0.1973d0 ! GeV
2598 RETURN
2599 END