]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA6/QPYTHIA/q-pyshow.1.2.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Partial Pythia to sync with master for missing files
[u/mrichter/AliRoot.git] / PYTHIA6 / QPYTHIA / q-pyshow.1.2.F
CommitLineData
6c43eccb 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)
249 save xkap2, xlkap2, xome, xlome, xspec, npkap, npome
250 character*1000 filnam
251 character*1000 chroot
252 DATA IFLAG/0/
253c 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)
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
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
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/
728
729Cacs+
730 common/qpc1/eee,qhatl,omegac
731 common/qpvir1/pmed
732 common/qpvir2/virt
733 COMMON/QPLT/QPLTA1,QPLTA2,QPLTBX,QPLTBY,QPLTBZ
734 external splitg1
735 external splitq1
736 external splitg2
737 external splitq2
738 external splitqqbar
739 data iflag/0/
740Cacs-
741C...Local arrays.
742 DIMENSION PMTH(5,140),PS(5),PMA(100),PMSD(100),IEP(100),IPA(100),
743 &KFLA(100),KFLD(100),KFL(100),ITRY(100),ISI(100),ISL(100),DP(100),
744 &DPT(5,4),KSH(0:140),KCII(2),NIIS(2),IIIS(2,2),THEIIS(2,2),
745 &PHIIIS(2,2),ISII(2),ISSET(2),ISCOL(0:140),ISCHG(0:140),
746 &IREF(1000)
747Cacs+
748
749 IF (IFLAG .EQ. 0) THEN
750 WRITE(MSTU(11),*)
751 WRITE(MSTU(11),*) '*******************************************'
752 WRITE(MSTU(11),*)
753 WRITE(MSTU(11),*) ' Q-PYTHIA version 1.0'
754 WRITE(MSTU(11),*)
755 WRITE(MSTU(11),*) 'DATE: 26.09.2008'
756 WRITE(MSTU(11),*)
757 WRITE(MSTU(11),*) 'AUTHORS: N. Armesto, L. Cunqueiro and'
758 WRITE(MSTU(11),*) ' C. A. Salgado'
759 WRITE(MSTU(11),*) ' Departamento de Fisica de Particulas'
760 WRITE(MSTU(11),*) ' and IGFAE'
761 WRITE(MSTU(11),*) ' Universidade de Santiago de Compostela'
762 WRITE(MSTU(11),*) ' 15706 Santiago de Compostela, Spain'
763 WRITE(MSTU(11),*)
764 WRITE(MSTU(11),*) 'EMAILS: nestor@fpaxp1.usc.es,'
765 WRITE(MSTU(11),*) ' leticia@fpaxp1.usc.es,'
766 WRITE(MSTU(11),*) ' Carlos.Salgado@cern.ch'
767 WRITE(MSTU(11),*)
768 WRITE(MSTU(11),*) 'NOTE: fixed to PYTHIA-6.4.18'
769 WRITE(MSTU(11),*)
770 WRITE(MSTU(11),*) 'WHEN USING Q-PYTHIA, PLEASE QUOTE:'
771 WRITE(MSTU(11),*) '1) N. Armesto, G. Corcella, L. Cunqueiro'
772 WRITE(MSTU(11),*) ' and C. A. Salgado, in preparation.'
773 WRITE(MSTU(11),*) '2) T. Sjostrand, S. Mrenna and P. Skands,'
774 WRITE(MSTU(11),*) ' PYTHIA 6.4 physics and manual,'
775 WRITE(MSTU(11),*) ' JHEP 0605 (2006) 026'
776 WRITE(MSTU(11),*) ' [arXiv:hep-ph/0603175].'
777 WRITE(MSTU(11),*)
778 WRITE(MSTU(11),*) 'INSTRUCTIONS: look at the web page and'
779 WRITE(MSTU(11),*) ' header of modfied routine PYSHOW at the'
780 WRITE(MSTU(11),*) ' end of Q-PYTHIA file.'
781 WRITE(MSTU(11),*)
782 WRITE(MSTU(11),*) 'DISCLAIMER: this program comes without any'
783 WRITE(MSTU(11),*) ' guarantees. Beware of errors and use'
784 WRITE(MSTU(11),*) ' common sense when interpreting results.'
785 WRITE(MSTU(11),*) ' Any modifications are done under exclusive'
786 WRITE(MSTU(11),*) ' makers resposibility.'
787 WRITE(MSTU(11),*)
788 WRITE(MSTU(11),*) '*******************************************'
789 WRITE(MSTU(11),*)
790 IFLAG=1
791 ENDIF
792Cacs-
793
794C...Check that QMAX not too low.
795 IF(MSTJ(41).LE.0) THEN
796 RETURN
797 ELSEIF(MSTJ(41).EQ.1.OR.MSTJ(41).EQ.11) THEN
798 IF(QMAX.LE.PARJ(82).AND.IP2.GE.-80) RETURN
799 ELSE
800 IF(QMAX.LE.MIN(PARJ(82),PARJ(83),PARJ(90)).AND.IP2.GE.-80)
801 & RETURN
802 ENDIF
803
804C...Store positions of shower initiating partons.
805 MPSPD=0
806 IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
807 NPA=1
808 IPA(1)=IP1
809 ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
810 & MSTU(32))) THEN
811 NPA=2
812 IPA(1)=IP1
813 IPA(2)=IP2
814 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0
815 & .AND.IP2.GE.-80) THEN
816 NPA=IABS(IP2)
817 DO 100 I=1,NPA
818 IPA(I)=IP1+I-1
819 100 CONTINUE
820 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.
821 &IP2.EQ.-100) THEN
822 MPSPD=1
823 NPA=2
824 IPA(1)=IP1+6
825 IPA(2)=IP1+7
826 ELSE
827 CALL PYERRM(12,
828 & '(PYSHOW:) failed to reconstruct showering system')
829 IF(MSTU(21).GE.1) RETURN
830 ENDIF
831
832C...Send off to PYPTFS for pT-ordered evolution if requested,
833C...if at least 2 partons, and without predefined shower branchings.
834 IF((MSTJ(41).EQ.11.OR.MSTJ(41).EQ.12).AND.NPA.GE.2.AND.
835 &MPSPD.EQ.0) THEN
836 NPART=NPA
837 DO 110 II=1,NPART
838 IPART(II)=IPA(II)
839 PTPART(II)=0.5D0*QMAX
840 110 CONTINUE
841 CALL PYPTFS(2,0.5D0*QMAX,0D0,PTGEN)
842 RETURN
843 ENDIF
844
845C...Initialization of cutoff masses etc.
846 DO 120 IFL=0,40
847 ISCOL(IFL)=0
848 ISCHG(IFL)=0
849 KSH(IFL)=0
850 120 CONTINUE
851 ISCOL(21)=1
852 KSH(21)=1
853 PMTH(1,21)=PYMASS(21)
854 PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25D0*PARJ(82)**2)
855 PMTH(3,21)=2D0*PMTH(2,21)
856 PMTH(4,21)=PMTH(3,21)
857 PMTH(5,21)=PMTH(3,21)
858 PMTH(1,22)=PYMASS(22)
859 PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25D0*PARJ(83)**2)
860 PMTH(3,22)=2D0*PMTH(2,22)
861 PMTH(4,22)=PMTH(3,22)
862 PMTH(5,22)=PMTH(3,22)
863 PMQTH1=PARJ(82)
864 IF(MSTJ(41).GE.2) PMQTH1=MIN(PARJ(82),PARJ(83))
865 PMQT1E=MIN(PMQTH1,PARJ(90))
866 PMQTH2=PMTH(2,21)
867 IF(MSTJ(41).GE.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
868 PMQT2E=MIN(PMQTH2,0.5D0*PARJ(90))
869 DO 130 IFL=1,5
870 ISCOL(IFL)=1
871 IF(MSTJ(41).GE.2) ISCHG(IFL)=1
872 KSH(IFL)=1
873 PMTH(1,IFL)=PYMASS(IFL)
874 PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PMQTH1**2)
875 PMTH(3,IFL)=PMTH(2,IFL)+PMQTH2
876 PMTH(4,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(82)**2)+PMTH(2,21)
877 PMTH(5,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(83)**2)+PMTH(2,22)
878 130 CONTINUE
879 DO 140 IFL=11,15,2
880 IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) ISCHG(IFL)=1
881 IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) KSH(IFL)=1
882 PMTH(1,IFL)=PYMASS(IFL)
883 PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(90)**2)
884 PMTH(3,IFL)=PMTH(2,IFL)+0.5D0*PARJ(90)
885 PMTH(4,IFL)=PMTH(3,IFL)
886 PMTH(5,IFL)=PMTH(3,IFL)
887 140 CONTINUE
888 PT2MIN=MAX(0.5D0*PARJ(82),1.1D0*PARJ(81))**2
889 ALAMS=PARJ(81)**2
890 ALFM=LOG(PT2MIN/ALAMS)
891
892C...Check on phase space available for emission.
893 IREJ=0
894 DO 150 J=1,5
895 PS(J)=0D0
896 150 CONTINUE
897 PM=0D0
898 KFLA(2)=0
899 DO 170 I=1,NPA
900 KFLA(I)=IABS(K(IPA(I),2))
901 PMA(I)=P(IPA(I),5)
902C...Special cutoff masses for initial partons (may be a heavy quark,
903C...squark, ..., and need not be on the mass shell).
904 IR=30+I
905 IF(NPA.LE.1) IREF(I)=IR
906 IF(NPA.GE.2) IREF(I+1)=IR
907 ISCOL(IR)=0
908 ISCHG(IR)=0
909 KSH(IR)=0
910 IF(KFLA(I).LE.8) THEN
911 ISCOL(IR)=1
912 IF(MSTJ(41).GE.2) ISCHG(IR)=1
913 ELSEIF(KFLA(I).EQ.11.OR.KFLA(I).EQ.13.OR.KFLA(I).EQ.15.OR.
914 & KFLA(I).EQ.17) THEN
915 IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) ISCHG(IR)=1
916 ELSEIF(KFLA(I).EQ.21) THEN
917 ISCOL(IR)=1
918 ELSEIF((KFLA(I).GE.KSUSY1+1.AND.KFLA(I).LE.KSUSY1+8).OR.
919 & (KFLA(I).GE.KSUSY2+1.AND.KFLA(I).LE.KSUSY2+8)) THEN
920 ISCOL(IR)=1
921 ELSEIF(KFLA(I).EQ.KSUSY1+21) THEN
922 ISCOL(IR)=1
923C...QUARKONIA+++
924C...same for QQ~[3S18]
925 ELSEIF(MSTP(148).GE.1.AND.(KFLA(I).EQ.9900443.OR.
926 & KFLA(I).EQ.9900553)) THEN
927 ISCOL(IR)=1
928C...QUARKONIA---
929 ENDIF
930 IF(ISCOL(IR).EQ.1.OR.ISCHG(IR).EQ.1) KSH(IR)=1
931 PMTH(1,IR)=PMA(I)
932 IF(ISCOL(IR).EQ.1.AND.ISCHG(IR).EQ.1) THEN
933 PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PMQTH1**2)
934 PMTH(3,IR)=PMTH(2,IR)+PMQTH2
935 PMTH(4,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)+PMTH(2,21)
936 PMTH(5,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(83)**2)+PMTH(2,22)
937 ELSEIF(ISCOL(IR).EQ.1) THEN
938 PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)
939 PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(82)
940 PMTH(4,IR)=PMTH(3,IR)
941 PMTH(5,IR)=PMTH(3,IR)
942 ELSEIF(ISCHG(IR).EQ.1) THEN
943 PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(90)**2)
944 PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(90)
945 PMTH(4,IR)=PMTH(3,IR)
946 PMTH(5,IR)=PMTH(3,IR)
947 ENDIF
948 IF(KSH(IR).EQ.1) PMA(I)=PMTH(3,IR)
949 PM=PM+PMA(I)
950 IF(KSH(IR).EQ.0.OR.PMA(I).GT.10D0*QMAX) IREJ=IREJ+1
951 DO 160 J=1,4
952 PS(J)=PS(J)+P(IPA(I),J)
953 160 CONTINUE
954 170 CONTINUE
955 IF(IREJ.EQ.NPA.AND.IP2.GE.-7) RETURN
956 PS(5)=SQRT(MAX(0D0,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
957 IF(NPA.EQ.1) PS(5)=PS(4)
958 IF(PS(5).LE.PM+PMQT1E) RETURN
959
960C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
961 KFSRCE=0
962 IF(IP2.LE.0) THEN
963 ELSEIF(K(IP1,3).EQ.K(IP2,3).AND.K(IP1,3).GT.0) THEN
964 KFSRCE=IABS(K(K(IP1,3),2))
965 ELSE
966 IPAR1=MAX(1,K(IP1,3))
967 IPAR2=MAX(1,K(IP2,3))
968 IF(K(IPAR1,3).EQ.K(IPAR2,3).AND.K(IPAR1,3).GT.0)
969 & KFSRCE=IABS(K(K(IPAR1,3),2))
970 ENDIF
971 ITYPES=0
972 IF(KFSRCE.GE.1.AND.KFSRCE.LE.8) ITYPES=1
973 IF(KFSRCE.GE.KSUSY1+1.AND.KFSRCE.LE.KSUSY1+8) ITYPES=2
974 IF(KFSRCE.GE.KSUSY2+1.AND.KFSRCE.LE.KSUSY2+8) ITYPES=2
975 IF(KFSRCE.GE.21.AND.KFSRCE.LE.24) ITYPES=3
976 IF(KFSRCE.GE.32.AND.KFSRCE.LE.34) ITYPES=3
977 IF(KFSRCE.EQ.25.OR.(KFSRCE.GE.35.AND.KFSRCE.LE.37)) ITYPES=4
978 IF(KFSRCE.GE.KSUSY1+22.AND.KFSRCE.LE.KSUSY1+37) ITYPES=5
979 IF(KFSRCE.EQ.KSUSY1+21) ITYPES=6
980
981C...Identify two primary showerers.
982 ITYPE1=0
983 IF(KFLA(1).GE.1.AND.KFLA(1).LE.8) ITYPE1=1
984 IF(KFLA(1).GE.KSUSY1+1.AND.KFLA(1).LE.KSUSY1+8) ITYPE1=2
985 IF(KFLA(1).GE.KSUSY2+1.AND.KFLA(1).LE.KSUSY2+8) ITYPE1=2
986 IF(KFLA(1).GE.21.AND.KFLA(1).LE.24) ITYPE1=3
987 IF(KFLA(1).GE.32.AND.KFLA(1).LE.34) ITYPE1=3
988 IF(KFLA(1).EQ.25.OR.(KFLA(1).GE.35.AND.KFLA(1).LE.37)) ITYPE1=4
989 IF(KFLA(1).GE.KSUSY1+22.AND.KFLA(1).LE.KSUSY1+37) ITYPE1=5
990 IF(KFLA(1).EQ.KSUSY1+21) ITYPE1=6
991 ITYPE2=0
992 IF(KFLA(2).GE.1.AND.KFLA(2).LE.8) ITYPE2=1
993 IF(KFLA(2).GE.KSUSY1+1.AND.KFLA(2).LE.KSUSY1+8) ITYPE2=2
994 IF(KFLA(2).GE.KSUSY2+1.AND.KFLA(2).LE.KSUSY2+8) ITYPE2=2
995 IF(KFLA(2).GE.21.AND.KFLA(2).LE.24) ITYPE2=3
996 IF(KFLA(2).GE.32.AND.KFLA(2).LE.34) ITYPE2=3
997 IF(KFLA(2).EQ.25.OR.(KFLA(2).GE.35.AND.KFLA(2).LE.37)) ITYPE2=4
998 IF(KFLA(2).GE.KSUSY1+22.AND.KFLA(2).LE.KSUSY1+37) ITYPE2=5
999 IF(KFLA(2).EQ.KSUSY1+21) ITYPE2=6
1000
1001C...Order of showerers. Presence of gluino.
1002 ITYPMN=MIN(ITYPE1,ITYPE2)
1003 ITYPMX=MAX(ITYPE1,ITYPE2)
1004 IORD=1
1005 IF(ITYPE1.GT.ITYPE2) IORD=2
1006 IGLUI=0
1007 IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1
1008
1009C...Check if 3-jet matrix elements to be used.
1010 M3JC=0
1011 ALPHA=0.5D0
1012 IF(NPA.EQ.2.AND.MSTJ(47).GE.1.AND.MPSPD.EQ.0) THEN
1013 IF(MSTJ(38).NE.0) THEN
1014 M3JC=MSTJ(38)
1015 ALPHA=PARJ(80)
1016 MSTJ(38)=0
1017 ELSEIF(MSTJ(47).GE.6) THEN
1018 M3JC=MSTJ(47)
1019 ELSE
1020 ICLASS=1
1021 ICOMBI=4
1022
1023C...Vector/axial vector -> q + qbar; q -> q + V.
1024 IF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.(ITYPES.EQ.0.OR.
1025 & ITYPES.EQ.3)) THEN
1026 ICLASS=2
1027 IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN
1028 ICOMBI=1
1029 ELSEIF(KFSRCE.EQ.23.OR.(KFSRCE.EQ.0.AND.
1030 & K(IPA(1),2)+K(IPA(2),2).EQ.0)) THEN
1031C...gamma*/Z0: assume e+e- initial state if unknown.
1032 EI=-1D0
1033 IF(KFSRCE.EQ.23) THEN
1034 IANNFL=K(K(IP1,3),3)
1035 IF(IANNFL.NE.0) THEN
1036 KANNFL=IABS(K(IANNFL,2))
1037 IF(KANNFL.GE.1.AND.KANNFL.LE.18) EI=KCHG(KANNFL,1)/3D0
1038 ENDIF
1039 ENDIF
1040 AI=SIGN(1D0,EI+0.1D0)
1041 VI=AI-4D0*EI*PARU(102)
1042 EF=KCHG(KFLA(1),1)/3D0
1043 AF=SIGN(1D0,EF+0.1D0)
1044 VF=AF-4D0*EF*PARU(102)
1045 XWC=1D0/(16D0*PARU(102)*(1D0-PARU(102)))
1046 SH=PS(5)**2
1047 SQMZ=PMAS(23,1)**2
1048 SQWZ=PS(5)*PMAS(23,2)
1049 SBWZ=1D0/((SH-SQMZ)**2+SQWZ**2)
1050 VECT=EI**2*EF**2+2D0*EI*VI*EF*VF*XWC*SH*(SH-SQMZ)*SBWZ+
1051 & (VI**2+AI**2)*VF**2*XWC**2*SH**2*SBWZ
1052 AXIV=(VI**2+AI**2)*AF**2*XWC**2*SH**2*SBWZ
1053 ICOMBI=3
1054 ALPHA=VECT/(VECT+AXIV)
1055 ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN
1056 ICOMBI=4
1057 ENDIF
1058C...For chi -> chi q qbar, use V/A -> q qbar as first approximation.
1059 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.5) THEN
1060 ICLASS=2
1061 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
1062 & ITYPES.EQ.1)) THEN
1063 ICLASS=3
1064
1065C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
1066 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN
1067 ICLASS=4
1068 IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN
1069 ICOMBI=1
1070 ELSEIF(KFSRCE.EQ.36) THEN
1071 ICOMBI=2
1072 ENDIF
1073 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
1074 & ITYPES.EQ.1)) THEN
1075 ICLASS=5
1076
1077C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
1078 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
1079 & ITYPES.EQ.3)) THEN
1080 ICLASS=6
1081 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
1082 & ITYPES.EQ.2)) THEN
1083 ICLASS=7
1084 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN
1085 ICLASS=8
1086 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
1087 & ITYPES.EQ.2)) THEN
1088 ICLASS=9
1089
1090C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
1091 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
1092 & ITYPES.EQ.5)) THEN
1093 ICLASS=10
1094 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
1095 & ITYPES.EQ.2)) THEN
1096 ICLASS=11
1097 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
1098 & ITYPES.EQ.1)) THEN
1099 ICLASS=12
1100
1101C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
1102 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.6) THEN
1103 ICLASS=13
1104 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
1105 & ITYPES.EQ.2)) THEN
1106 ICLASS=14
1107 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
1108 & ITYPES.EQ.1)) THEN
1109 ICLASS=15
1110
1111C...g -> ~g + ~g (eikonal approximation).
1112 ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN
1113 ICLASS=16
1114 ENDIF
1115 M3JC=5*ICLASS+ICOMBI
1116 ENDIF
1117 ENDIF
1118
1119C...Find if interference with initial state partons.
1120 MIIS=0
1121 IF(MSTJ(50).GE.1.AND.MSTJ(50).LE.3.AND.NPA.EQ.2.AND.KFSRCE.EQ.0
1122 &.AND.MPSPD.EQ.0) MIIS=MSTJ(50)
1123 IF(MSTJ(50).GE.4.AND.MSTJ(50).LE.6.AND.NPA.EQ.2.AND.MPSPD.EQ.0)
1124 &MIIS=MSTJ(50)-3
1125 IF(MIIS.NE.0) THEN
1126 DO 190 I=1,2
1127 KCII(I)=0
1128 KCA=PYCOMP(KFLA(I))
1129 IF(KCA.NE.0) KCII(I)=KCHG(KCA,2)*ISIGN(1,K(IPA(I),2))
1130 NIIS(I)=0
1131 IF(KCII(I).NE.0) THEN
1132 DO 180 J=1,2
1133 ICSI=MOD(K(IPA(I),3+J)/MSTU(5),MSTU(5))
1134 IF(ICSI.GT.0.AND.ICSI.NE.IPA(1).AND.ICSI.NE.IPA(2).AND.
1135 & (KCII(I).EQ.(-1)**(J+1).OR.KCII(I).EQ.2)) THEN
1136 NIIS(I)=NIIS(I)+1
1137 IIIS(I,NIIS(I))=ICSI
1138 ENDIF
1139 180 CONTINUE
1140 ENDIF
1141 190 CONTINUE
1142 IF(NIIS(1)+NIIS(2).EQ.0) MIIS=0
1143 ENDIF
1144
1145C...Boost interfering initial partons to rest frame
1146C...and reconstruct their polar and azimuthal angles.
1147Cacs+
1148 qplta1=0.d0
1149 qplta2=0.d0
1150 qpltbx=0.d0
1151 qpltby=0.d0
1152 qpltbz=0.d0
1153Cacs-
1154 IF(MIIS.NE.0) THEN
1155 DO 210 I=1,2
1156 DO 200 J=1,5
1157 K(N+I,J)=K(IPA(I),J)
1158 P(N+I,J)=P(IPA(I),J)
1159 V(N+I,J)=0D0
1160 200 CONTINUE
1161 210 CONTINUE
1162 DO 230 I=3,2+NIIS(1)
1163 DO 220 J=1,5
1164 K(N+I,J)=K(IIIS(1,I-2),J)
1165 P(N+I,J)=P(IIIS(1,I-2),J)
1166 V(N+I,J)=0D0
1167 220 CONTINUE
1168 230 CONTINUE
1169 DO 250 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
1170 DO 240 J=1,5
1171 K(N+I,J)=K(IIIS(2,I-2-NIIS(1)),J)
1172 P(N+I,J)=P(IIIS(2,I-2-NIIS(1)),J)
1173 V(N+I,J)=0D0
1174 240 CONTINUE
1175 250 CONTINUE
1176 CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,0D0,-PS(1)/PS(4),
1177 & -PS(2)/PS(4),-PS(3)/PS(4))
1178 PHI=PYANGL(P(N+1,1),P(N+1,2))
1179 CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,-PHI,0D0,0D0,0D0)
1180 THE=PYANGL(P(N+1,3),P(N+1,1))
1181 CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),-THE,0D0,0D0,0D0,0D0)
1182Cacs+
1183 qplta1=-the
1184 qplta2=-phi
1185 qpltbx=-PS(1)/PS(4)
1186 qpltby=-PS(2)/PS(4)
1187 qpltbz=-PS(3)/PS(4)
1188Cacs-
1189 DO 260 I=3,2+NIIS(1)
1190 THEIIS(1,I-2)=PYANGL(P(N+I,3),SQRT(P(N+I,1)**2+P(N+I,2)**2))
1191 PHIIIS(1,I-2)=PYANGL(P(N+I,1),P(N+I,2))
1192 260 CONTINUE
1193 DO 270 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
1194 THEIIS(2,I-2-NIIS(1))=PARU(1)-PYANGL(P(N+I,3),
1195 & SQRT(P(N+I,1)**2+P(N+I,2)**2))
1196 PHIIIS(2,I-2-NIIS(1))=PYANGL(P(N+I,1),P(N+I,2))
1197 270 CONTINUE
1198 ENDIF
1199
1200C...Boost 3 or more partons to their rest frame.
1201Cacs+
1202c IF(NPA.GE.3) CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
1203c &-PS(2)/PS(4),-PS(3)/PS(4))
1204 IF(NPA.GE.3) THEN
1205 CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
1206 &-PS(2)/PS(4),-PS(3)/PS(4))
1207 qplta1=0.d0
1208 qplta2=0.d0
1209 qpltbx=-PS(1)/PS(4)
1210 qpltby=-PS(2)/PS(4)
1211 qpltbz=-PS(3)/PS(4)
1212
1213 ENDIF
1214Cacs-
1215
1216C...Define imagined single initiator of shower for parton system.
1217 NS=N
1218 IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
1219 CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
1220 IF(MSTU(21).GE.1) RETURN
1221 ENDIF
1222 280 N=NS
1223 IF(NPA.GE.2) THEN
1224 K(N+1,1)=11
1225 K(N+1,2)=21
1226 K(N+1,3)=0
1227 K(N+1,4)=0
1228 K(N+1,5)=0
1229 P(N+1,1)=0D0
1230 P(N+1,2)=0D0
1231 P(N+1,3)=0D0
1232 P(N+1,4)=PS(5)
1233 P(N+1,5)=PS(5)
1234 V(N+1,5)=PS(5)**2
1235 N=N+1
1236 IREF(1)=21
1237 ENDIF
1238
1239
1240
1241
1242Cacs+
1243 call qpygin(pposx0,pposy0,pposz0,ppost0) ! in fm
1244 do 10101 iijj=1, nnpos, 1
1245 ppos(iijj,1)=pposx0
1246 ppos(iijj,2)=pposy0
1247 ppos(iijj,3)=pposz0
1248 ppos(iijj,4)=ppost0
124910101 continue
1250
1251Cacs-
1252
1253
1254
1255
1256C...Loop over partons that may branch.
1257 NEP=NPA
1258 IM=NS
1259 IF(NPA.EQ.1) IM=NS-1
1260 290 IM=IM+1
1261 IF(N.GT.NS) THEN
1262 IF(IM.GT.N) GOTO 600
1263 KFLM=IABS(K(IM,2))
1264 IR=IREF(IM-NS)
1265 IF(KSH(IR).EQ.0) GOTO 290
1266 IF(P(IM,5).LT.PMTH(2,IR)) GOTO 290
1267 IGM=K(IM,3)
1268 ELSE
1269 IGM=-1
1270 ENDIF
1271 IF(N+NEP.GT.MSTU(4)-MSTU(32)-10) THEN
1272 CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
1273 IF(MSTU(21).GE.1) RETURN
1274 ENDIF
1275
1276C...Position of aunt (sister to branching parton).
1277C...Origin and flavour of daughters.
1278 IAU=0
1279 IF(IGM.GT.0) THEN
1280 IF(K(IM-1,3).EQ.IGM) IAU=IM-1
1281 IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
1282 ENDIF
1283 IF(IGM.GE.0) THEN
1284 K(IM,4)=N+1
1285 DO 300 I=1,NEP
1286 K(N+I,3)=IM
1287 300 CONTINUE
1288 ELSE
1289 K(N+1,3)=IPA(1)
1290 ENDIF
1291 IF(IGM.LE.0) THEN
1292 DO 310 I=1,NEP
1293 K(N+I,2)=K(IPA(I),2)
1294 310 CONTINUE
1295 ELSEIF(KFLM.NE.21) THEN
1296 K(N+1,2)=K(IM,2)
1297 K(N+2,2)=K(IM,5)
1298 IREF(N+1-NS)=IREF(IM-NS)
1299 IREF(N+2-NS)=IABS(K(N+2,2))
1300 ELSEIF(K(IM,5).EQ.21) THEN
1301 K(N+1,2)=21
1302 K(N+2,2)=21
1303 IREF(N+1-NS)=21
1304 IREF(N+2-NS)=21
1305 ELSE
1306 K(N+1,2)=K(IM,5)
1307 K(N+2,2)=-K(IM,5)
1308 IREF(N+1-NS)=IABS(K(N+1,2))
1309 IREF(N+2-NS)=IABS(K(N+2,2))
1310 ENDIF
1311
1312C...Reset flags on daughters and tries made.
1313 DO 320 IP=1,NEP
1314 K(N+IP,1)=3
1315 K(N+IP,4)=0
1316 K(N+IP,5)=0
1317 KFLD(IP)=IABS(K(N+IP,2))
1318 IF(KCHG(PYCOMP(KFLD(IP)),2).EQ.0) K(N+IP,1)=1
1319 ITRY(IP)=0
1320 ISL(IP)=0
1321 ISI(IP)=0
1322 IF(KSH(IREF(N+IP-NS)).EQ.1) ISI(IP)=1
1323 320 CONTINUE
1324 ISLM=0
1325
1326C...Maximum virtuality of daughters.
1327 IF(IGM.LE.0) THEN
1328 DO 330 I=1,NPA
1329 IF(NPA.GE.3) P(N+I,4)=P(IPA(I),4)
1330 P(N+I,5)=MIN(QMAX,PS(5))
1331 IR=IREF(N+I-NS)
1332 IF(IP2.LE.-8) P(N+I,5)=MAX(P(N+I,5),2D0*PMTH(3,IR))
1333 IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
1334 330 CONTINUE
1335 ELSE
1336 IF(MSTJ(43).LE.2) PEM=V(IM,2)
1337 IF(MSTJ(43).GE.3) PEM=P(IM,4)
1338 P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
1339 P(N+2,5)=MIN(P(IM,5),(1D0-V(IM,1))*PEM)
1340 IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
1341 ENDIF
1342 DO 340 I=1,NEP
1343 PMSD(I)=P(N+I,5)
1344 IF(ISI(I).EQ.1) THEN
1345 IR=IREF(N+I-NS)
1346 IF(P(N+I,5).LE.PMTH(3,IR)) P(N+I,5)=PMTH(1,IR)
1347 ENDIF
1348 V(N+I,5)=P(N+I,5)**2
1349 340 CONTINUE
1350
1351C...Choose one of the daughters for evolution.
1352 350 INUM=0
1353 IF(NEP.EQ.1) INUM=1
1354 DO 360 I=1,NEP
1355 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
1356 360 CONTINUE
1357 DO 370 I=1,NEP
1358 IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
1359 IR=IREF(N+I-NS)
1360 IF(P(N+I,5).GE.PMTH(2,IR)) INUM=I
1361 ENDIF
1362 370 CONTINUE
1363 IF(INUM.EQ.0) THEN
1364 RMAX=0D0
1365 DO 380 I=1,NEP
1366 IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQT2E) THEN
1367 RPM=P(N+I,5)/PMSD(I)
1368 IR=IREF(N+I-NS)
1369 IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,IR)) THEN
1370 RMAX=RPM
1371 INUM=I
1372 ENDIF
1373 ENDIF
1374 380 CONTINUE
1375 ENDIF
1376
1377C...Cancel choice of predetermined daughter already treated.
1378 INUM=MAX(1,INUM)
1379 INUMT=INUM
1380 IF(MPSPD.EQ.1.AND.IGM.EQ.0.AND.ITRY(INUMT).GE.1) THEN
1381 IF(K(IP1-1+INUM,4).GT.0) INUM=3-INUM
1382 ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2.AND.ITRY(INUMT).GE.1) THEN
1383 IF(KFLD(INUMT).NE.21.AND.K(IP1+2,4).GT.0) INUM=3-INUM
1384 IF(KFLD(INUMT).EQ.21.AND.K(IP1+3,4).GT.0) INUM=3-INUM
1385 ENDIF
1386
1387C...Store information on choice of evolving daughter.
1388 IEP(1)=N+INUM
1389Cacs+
1390 idf=k(iep(1),3)
1391 zz1=v(idf,1)
1392 zzz=zz1
1393 zz2=1.d0-zz1
1394 if (nep .gt. 1 .and. inum .eq. 2) then
1395 zzz=zz2
1396 endif
1397 ttt=v(idf,5)
1398 if(zz1.gt.0.d0) then
1399 eee=zzz*p(idf,4)
1400 else
1401 eee=p(idf,4)
1402 endif
1403 xkt=zz1*zz2*ttt
1404 if (xkt .gt. 0.d0) then
1405 xlcoh=(2.d0*eee/(zz1*zz2*ttt))*0.1973d0
1406 else
1407 xlcoh=0.d0
1408 endif
1409 if (idf .eq. 0) then ! for the initial parton if it has no father
1410 xbx=p(iep(1),1)/p(iep(1),4)
1411 xby=p(iep(1),2)/p(iep(1),4)
1412 xbz=p(iep(1),3)/p(iep(1),4)
1413 call qpygeo(pposx0,pposy0,pposz0,ppost0,
1414 > xbx,xby,xbz,qhatl,omegac)
1415 else
1416 xbx=p(idf,1)/p(idf,4)
1417 xby=p(idf,2)/p(idf,4)
1418 xbz=p(idf,3)/p(idf,4)
1419
1420
1421
1422 ppos(iep(1),1)=ppos(idf,1)+xbx*xlcoh
1423 ppos(iep(1),2)=ppos(idf,2)+xby*xlcoh
1424 ppos(iep(1),3)=ppos(idf,3)+xbz*xlcoh
1425 ppos(iep(1),4)=ppos(idf,4)+xlcoh
1426 call qpygeo(ppos(iep(1),1),ppos(iep(1),2),ppos(iep(1),3),
1427 > ppos(iep(1),4),xbx,xby,xbz,qhatl,omegac)
1428 endif
1429Cacs-
1430
1431 DO 390 I=2,NEP
1432 IEP(I)=IEP(I-1)+1
1433 IF(IEP(I).GT.N+NEP) IEP(I)=N+1
1434 390 CONTINUE
1435 DO 400 I=1,NEP
1436 KFL(I)=IABS(K(IEP(I),2))
1437 400 CONTINUE
1438 ITRY(INUM)=ITRY(INUM)+1
1439 IF(ITRY(INUM).GT.200) THEN
1440 CALL PYERRM(14,'(PYSHOW:) caught in infinite loop')
1441 IF(MSTU(21).GE.1) RETURN
1442 ENDIF
1443 Z=0.5D0
1444 IR=IREF(IEP(1)-NS)
1445 IF(KSH(IR).EQ.0) GOTO 450
1446 IF(P(IEP(1),5).LT.PMTH(2,IR)) GOTO 450
1447
1448C...Check if evolution already predetermined for daughter.
1449 IPSPD=0
1450 IF(MPSPD.EQ.1.AND.IGM.EQ.0) THEN
1451 IF(K(IP1-1+INUM,4).GT.0) IPSPD=IP1-1+INUM
1452 ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2) THEN
1453 IF(KFL(1).NE.21.AND.K(IP1+2,4).GT.0) IPSPD=IP1+2
1454 IF(KFL(1).EQ.21.AND.K(IP1+3,4).GT.0) IPSPD=IP1+3
1455 ENDIF
1456 IF(INUM.EQ.1.OR.INUM.EQ.2) THEN
1457 ISSET(INUM)=0
1458 IF(IPSPD.NE.0) ISSET(INUM)=1
1459 ENDIF
1460
1461C...Select side for interference with initial state partons.
1462 IF(MIIS.GE.1.AND.IEP(1).LE.NS+3) THEN
1463 III=IEP(1)-NS-1
1464 ISII(III)=0
1465 IF(IABS(KCII(III)).EQ.1.AND.NIIS(III).EQ.1) THEN
1466 ISII(III)=1
1467 ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.1) THEN
1468 IF(PYR(0).GT.0.5D0) ISII(III)=1
1469 ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.2) THEN
1470 ISII(III)=1
1471 IF(PYR(0).GT.0.5D0) ISII(III)=2
1472 ENDIF
1473 ENDIF
1474
1475C...Calculate allowed z range.
1476 IF(NEP.EQ.1) THEN
1477 PMED=PS(4)
1478 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1479 PMED=P(IM,5)
1480 ELSE
1481 IF(INUM.EQ.1) PMED=V(IM,1)*PEM
1482 IF(INUM.EQ.2) PMED=(1D0-V(IM,1))*PEM
1483 ENDIF
1484 IF(MOD(MSTJ(43),2).EQ.1) THEN
1485 ZC=PMTH(2,21)/PMED
1486 ZCE=PMTH(2,22)/PMED
1487 IF(ISCOL(IR).EQ.0) ZCE=0.5D0*PARJ(90)/PMED
1488 ELSE
1489 ZC=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTH(2,21)/PMED)**2)))
1490 IF(ZC.LT.1D-6) ZC=(PMTH(2,21)/PMED)**2
1491 PMTMPE=PMTH(2,22)
1492 IF(ISCOL(IR).EQ.0) PMTMPE=0.5D0*PARJ(90)
1493 ZCE=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTMPE/PMED)**2)))
1494 IF(ZCE.LT.1D-6) ZCE=(PMTMPE/PMED)**2
1495 ENDIF
1496 ZC=MIN(ZC,0.491D0)
1497 ZCE=MIN(ZCE,0.49991D0)
1498 IF(((MSTJ(41).EQ.1.AND.ZC.GT.0.49D0).OR.(MSTJ(41).GE.2.AND.
1499 &MIN(ZC,ZCE).GT.0.4999D0)).AND.IPSPD.EQ.0) THEN
1500 P(IEP(1),5)=PMTH(1,IR)
1501 V(IEP(1),5)=P(IEP(1),5)**2
1502 GOTO 450
1503 ENDIF
1504
1505C...Integral of Altarelli-Parisi z kernel for QCD.
1506C...(Includes squark and gluino; with factor N_C/C_F extra for latter).
1507 IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
1508Cacs+
1509C FBR= 6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
1510 FBR=dgauss1(splitg1,zc,1.d0-zc,1.d-3)
1511Cacs-
1512C...QUARKONIA+++
1513C...Evolution of QQ~[3S18] state if MSTP(148)=1.
1514 ELSEIF(MSTJ(49).EQ.0.AND.MSTP(149).GE.0.AND.
1515 & (KFL(1).EQ.9900443.OR.KFL(1).EQ.9900553)) THEN
1516 FBR=6D0*LOG((1D0-ZC)/ZC)
1517C...QUARKONIA---
1518 ELSEIF(MSTJ(49).EQ.0) THEN
1519Cacs+
1520C FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC)
1521 FBR=dgauss1(splitq1,zc,1.d0-zc,1.d-3)
1522
1523Cacs-
1524 IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0)
1525
1526C...Integral of Altarelli-Parisi z kernel for scalar gluon.
1527 ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
1528 FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1D0-2D0*ZC)
1529 ELSEIF(MSTJ(49).EQ.1) THEN
1530 FBR=(1D0-2D0*ZC)/3D0
1531 IF(IGM.EQ.0.AND.M3JC.GE.1) FBR=4D0*FBR
1532
1533C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
1534 ELSEIF(KFL(1).EQ.21) THEN
1535 FBR=6D0*MSTJ(45)*(0.5D0-ZC)
1536 ELSE
1537 FBR=2D0*LOG((1D0-ZC)/ZC)
1538 ENDIF
1539
1540C...Reset QCD probability for colourless.
1541 IF(ISCOL(IR).EQ.0) FBR=0D0
1542
1543C...Integral of Altarelli-Parisi kernel for photon emission.
1544 FBRE=0D0
1545 IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1) THEN
1546 IF(KFL(1).LE.18) THEN
1547 FBRE=(KCHG(KFL(1),1)/3D0)**2*2D0*LOG((1D0-ZCE)/ZCE)
1548 ENDIF
1549 IF(MSTJ(41).EQ.10) FBRE=PARJ(84)*FBRE
1550 ENDIF
1551
1552C...Inner veto algorithm starts. Find maximum mass for evolution.
1553 410 PMS=V(IEP(1),5)
1554 IF(IGM.GE.0) THEN
1555 PM2=0D0
1556 DO 420 I=2,NEP
1557 PM=P(IEP(I),5)
1558 IRI=IREF(IEP(I)-NS)
1559 IF(KSH(IRI).EQ.1) PM=PMTH(2,IRI)
1560 PM2=PM2+PM
1561 420 CONTINUE
1562 PMS=MIN(PMS,(P(IM,5)-PM2)**2)
1563 ENDIF
1564
1565C...Select mass for daughter in QCD evolution.
1566 B0=27D0/6D0
1567 DO 430 IFF=4,MSTJ(45)
1568 IF(PMS.GT.4D0*PMTH(2,IFF)**2) B0=(33D0-2D0*IFF)/6D0
1569 430 CONTINUE
1570C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
1571 PMSC=MAX(0.5D0*PARJ(82),PMS-PMTH(1,IR)**2)
1572C...Already predetermined choice.
1573 IF(IPSPD.NE.0) THEN
1574 PMSQCD=P(IPSPD,5)**2
1575 ELSEIF(FBR.LT.1D-3) THEN
1576 PMSQCD=0D0
1577 ELSEIF(MSTJ(44).LE.0) THEN
1578 PMSQCD=PMSC*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/(PARU(111)*FBR)))
1579 ELSEIF(MSTJ(44).EQ.1) THEN
1580 PMSQCD=4D0*ALAMS*(0.25D0*PMSC/ALAMS)**(PYR(0)**(B0/FBR))
1581 ELSE
1582 PMSQCD=PMSC*EXP(MAX(-50D0,ALFM*B0*LOG(PYR(0))/FBR))
1583 ENDIF
1584C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
1585 IF(IPSPD.EQ.0) PMSQCD=PMSQCD+PMTH(1,IR)**2
1586 IF(ZC.GT.0.49D0.OR.PMSQCD.LE.PMTH(4,IR)**2) PMSQCD=PMTH(2,IR)**2
1587 V(IEP(1),5)=PMSQCD
1588 MCE=1
1589
1590C...Select mass for daughter in QED evolution.
1591 IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1.AND.IPSPD.EQ.0) THEN
1592C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
1593 PMSE=MAX(0.5D0*PARJ(83),PMS-PMTH(1,IR)**2)
1594 IF(FBRE.LT.1D-3) THEN
1595 PMSQED=0D0
1596 ELSE
1597 PMSQED=PMSE*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/
1598 & (PARU(101)*FBRE)))
1599 ENDIF
1600C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
1601 PMSQED=PMSQED+PMTH(1,IR)**2
1602 IF(ZCE.GT.0.4999D0.OR.PMSQED.LE.PMTH(5,IR)**2) PMSQED=
1603 & PMTH(2,IR)**2
1604 IF(PMSQED.GT.PMSQCD) THEN
1605 V(IEP(1),5)=PMSQED
1606 MCE=2
1607 ENDIF
1608 ENDIF
1609
1610C...Check whether daughter mass below cutoff.
1611 P(IEP(1),5)=SQRT(V(IEP(1),5))
1612 IF(P(IEP(1),5).LE.PMTH(3,IR)) THEN
1613 P(IEP(1),5)=PMTH(1,IR)
1614 V(IEP(1),5)=P(IEP(1),5)**2
1615 GOTO 450
1616 ENDIF
1617Cacs+
1618 virt=V(IEP(1),5)
1619Cacs-
1620
1621C...Already predetermined choice of z, and flavour in g -> qqbar.
1622 IF(IPSPD.NE.0) THEN
1623 IPSGD1=K(IPSPD,4)
1624 IPSGD2=K(IPSPD,5)
1625 PMSGD1=P(IPSGD1,5)**2
1626 PMSGD2=P(IPSGD2,5)**2
1627 ALAMPS=SQRT(MAX(1D-10,(PMSQCD-PMSGD1-PMSGD2)**2-
1628 & 4D0*PMSGD1*PMSGD2))
1629 Z=0.5D0*(PMSQCD*(2D0*P(IPSGD1,4)/P(IPSPD,4)-1D0)+ALAMPS-
1630 & PMSGD1+PMSGD2)/ALAMPS
1631 Z=MAX(0.00001D0,MIN(0.99999D0,Z))
1632 IF(KFL(1).NE.21) THEN
1633 K(IEP(1),5)=21
1634 ELSE
1635 K(IEP(1),5)=IABS(K(IPSGD1,2))
1636 ENDIF
1637
1638C...Select z value of branching: q -> qgamma.
1639 ELSEIF(MCE.EQ.2) THEN
1640 Z=1D0-(1D0-ZCE)*(ZCE/(1D0-ZCE))**PYR(0)
1641 IF(1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
1642 K(IEP(1),5)=22
1643
1644C...QUARKONIA+++
1645C...Select z value of branching: QQ~[3S18] -> QQ~[3S18]g.
1646 ELSEIF(MSTJ(49).EQ.0.AND.
1647 & (KFL(1).EQ.9900443.OR.KFL(1).EQ.9900553)) THEN
1648 Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1649C...Select always the harder 'gluon' if the switch MSTP(149)<=0.
1650 IF(MSTP(149).LE.0.OR.PYR(0).GT.0.5D0) Z=1D0-Z
1651 IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
1652 K(IEP(1),5)=21
1653C...QUARKONIA---
1654
1655C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
1656 ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
1657Cacs+
1658C Z=1D0-(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1659C...Only do z weighting when no ME correction afterwards.
1660C IF(M3JC.EQ.0.AND.1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
1661C
1662 anfbr=dgauss1(splitq2,zc,1.d0-zc,1.d-3)
1663 z=simdis(500,zc,1,anfbr)
1664Cacs-
1665 K(IEP(1),5)=21
1666 ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN
1667Cacs+
1668c Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1669 anfbr=dgauss1(splitg2,zc,1.d0-zc,1.d-3)
1670
1671 z=simdis(500,zc,21,anfbr)
1672
1673 IF(PYR(0).GT.0.5D0) Z=1D0-Z
1674c IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
1675Cacs-
1676 K(IEP(1),5)=21
1677 ELSEIF(MSTJ(49).NE.1) THEN
1678Cacs+
1679c Z=PYR(0)
1680c IF(Z**2+(1D0-Z)**2.LT.PYR(0)) GOTO 410
1681 anfbr=dgauss1(splitqqbar,zc,1.d0-zc,1.d-3)
1682 z=simdis(500,zc,3,anfbr)
1683
1684Cacs-
1685 KFLB=1+INT(MSTJ(45)*PYR(0))
1686 PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5)
1687 IF(PMQ.GE.1D0) GOTO 410
1688 IF(MSTJ(44).LE.2.OR.MSTJ(44).EQ.4) THEN
1689 IF(Z.LT.ZC.OR.Z.GT.1D0-ZC) GOTO 410
1690 PMQ0=4D0*PMTH(2,21)**2/V(IEP(1),5)
1691 IF(MOD(MSTJ(43),2).EQ.0.AND.(1D0+0.5D0*PMQ)*SQRT(1D0-PMQ)
1692 & .LT.PYR(0)*(1D0+0.5D0*PMQ0)*SQRT(1D0-PMQ0)) GOTO 410
1693 ELSE
1694 IF((1D0+0.5D0*PMQ)*SQRT(1D0-PMQ).LT.PYR(0)) GOTO 410
1695 ENDIF
1696 K(IEP(1),5)=KFLB
1697
1698C...Ditto for scalar gluon model.
1699 ELSEIF(KFL(1).NE.21) THEN
1700 Z=1D0-SQRT(ZC**2+PYR(0)*(1D0-2D0*ZC))
1701 K(IEP(1),5)=21
1702 ELSEIF(PYR(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
1703 Z=ZC+(1D0-2D0*ZC)*PYR(0)
1704 K(IEP(1),5)=21
1705 ELSE
1706 Z=ZC+(1D0-2D0*ZC)*PYR(0)
1707 KFLB=1+INT(MSTJ(45)*PYR(0))
1708 PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5)
1709 IF(PMQ.GE.1D0) GOTO 410
1710 K(IEP(1),5)=KFLB
1711 ENDIF
1712
1713C...Correct to alpha_s(pT^2) (optionally m^2/4 for g -> q qbar).
1714 IF(MCE.EQ.1.AND.MSTJ(44).GE.2.AND.IPSPD.EQ.0) THEN
1715 IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND.
1716 & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1717 IF(ALFM/LOG(V(IEP(1),5)*0.25D0/ALAMS).LT.PYR(0)) GOTO 410
1718 ELSE
1719 PT2APP=Z*(1D0-Z)*V(IEP(1),5)
1720 IF(MSTJ(44).GE.4) PT2APP=PT2APP*
1721 & (1D0-PMTH(1,IR)**2/V(IEP(1),5))**2
1722 IF(PT2APP.LT.PT2MIN) GOTO 410
1723 IF(ALFM/LOG(PT2APP/ALAMS).LT.PYR(0)) GOTO 410
1724 ENDIF
1725 ENDIF
1726
1727C...Check if z consistent with chosen m.
1728 IF(KFL(1).EQ.21) THEN
1729 IRGD1=IABS(K(IEP(1),5))
1730 IRGD2=IRGD1
1731 ELSE
1732 IRGD1=IR
1733 IRGD2=IABS(K(IEP(1),5))
1734 ENDIF
1735 IF(NEP.EQ.1) THEN
1736 PED=PS(4)
1737 ELSEIF(NEP.GE.3) THEN
1738 PED=P(IEP(1),4)
1739 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1740 PED=0.5D0*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
1741 ELSE
1742 IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
1743 IF(IEP(1).EQ.N+2) PED=(1D0-V(IM,1))*PEM
1744 ENDIF
1745 IF(MOD(MSTJ(43),2).EQ.1) THEN
1746 PMQTH3=0.5D0*PARJ(82)
1747 IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83)
1748 IF(IRGD2.EQ.22.AND.ISCOL(IR).EQ.0) PMQTH3=0.5D0*PARJ(90)
1749 PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(IEP(1),5)
1750 PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(IEP(1),5)
1751 ZD=SQRT(MAX(0D0,(1D0-V(IEP(1),5)/PED**2)*((1D0-PMQ1-PMQ2)**2-
1752 & 4D0*PMQ1*PMQ2)))
1753 ZH=1D0+PMQ1-PMQ2
1754 ELSE
1755 ZD=SQRT(MAX(0D0,1D0-V(IEP(1),5)/PED**2))
1756 ZH=1D0
1757 ENDIF
1758 IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND.
1759 &(MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1760 ELSEIF(IPSPD.NE.0) THEN
1761 ELSE
1762 ZL=0.5D0*(ZH-ZD)
1763 ZU=0.5D0*(ZH+ZD)
1764 IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 410
1765 ENDIF
1766 IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1D0-ZL)/MAX(1D-20,ZL*
1767 &(1D0-ZU)))
1768 IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
1769
1770C...Width suppression for q -> q + g.
1771 IF(MSTJ(40).NE.0.AND.KFL(1).NE.21.AND.IPSPD.EQ.0) THEN
1772 IF(IGM.EQ.0) THEN
1773 EGLU=0.5D0*PS(5)*(1D0-Z)*(1D0+V(IEP(1),5)/V(NS+1,5))
1774 ELSE
1775 EGLU=PMED*(1D0-Z)
1776 ENDIF
1777 CHI=PARJ(89)**2/(PARJ(89)**2+EGLU**2)
1778 IF(MSTJ(40).EQ.1) THEN
1779 IF(CHI.LT.PYR(0)) GOTO 410
1780 ELSEIF(MSTJ(40).EQ.2) THEN
1781 IF(1D0-CHI.LT.PYR(0)) GOTO 410
1782 ENDIF
1783 ENDIF
1784
1785C...Three-jet matrix element correction.
1786 IF(M3JC.GE.1) THEN
1787 WME=1D0
1788 WSHOW=1D0
1789
1790C...QED matrix elements: only for massless case so far.
1791 IF(MCE.EQ.2.AND.IGM.EQ.0) THEN
1792 X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5))
1793 X2=1D0-V(IEP(1),5)/V(NS+1,5)
1794 X3=(1D0-X1)+(1D0-X2)
1795 KI1=K(IPA(INUM),2)
1796 KI2=K(IPA(3-INUM),2)
1797 QF1=KCHG(PYCOMP(KI1),1)*ISIGN(1,KI1)/3D0
1798 QF2=KCHG(PYCOMP(KI2),1)*ISIGN(1,KI2)/3D0
1799 WSHOW=QF1**2*(1D0-X1)/X3*(1D0+(X1/(2D0-X2))**2)+
1800 & QF2**2*(1D0-X2)/X3*(1D0+(X2/(2D0-X1))**2)
1801 WME=(QF1*(1D0-X1)/X3-QF2*(1D0-X2)/X3)**2*(X1**2+X2**2)
1802 ELSEIF(MCE.EQ.2) THEN
1803
1804C...QCD matrix elements, including mass effects.
1805 ELSEIF(MSTJ(49).NE.1.AND.K(IEP(1),2).NE.21) THEN
1806 PS1ME=V(IEP(1),5)
1807 PM1ME=PMTH(1,IR)
1808 M3JCC=M3JC
1809 IF(IR.GE.31.AND.IGM.EQ.0) THEN
1810C...QCD ME: original parton, first branching.
1811 PM2ME=PMTH(1,63-IR)
1812 ECMME=PS(5)
1813 ELSEIF(IR.GE.31) THEN
1814C...QCD ME: original parton, subsequent branchings.
1815 PM2ME=PMTH(1,63-IR)
1816 PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5))
1817 ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1818 ELSEIF(K(IM,2).EQ.21) THEN
1819C...QCD ME: secondary partons, first branching.
1820 PM2ME=PM1ME
1821 ZMME=V(IM,1)
1822 IF(IEP(1).GT.IEP(2)) ZMME=1D0-ZMME
1823 PMLME=SQRT(MAX(0D0,(V(IM,5)-PS1ME-PM2ME**2)**2-
1824 & 4D0*PS1ME*PM2ME**2))
1825 PEDME=PEM*(0.5D0*(V(IM,5)-PMLME+PS1ME-PM2ME**2)+PMLME*ZMME)/
1826 & V(IM,5)
1827 ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1828 M3JCC=66
1829 ELSE
1830C...QCD ME: secondary partons, subsequent branchings.
1831 PM2ME=PM1ME
1832 PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5))
1833 ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1834 M3JCC=66
1835 ENDIF
1836C...Construct ME variables.
1837 R1ME=PM1ME/ECMME
1838 R2ME=PM2ME/ECMME
1839 X1=(1D0+PS1ME/ECMME**2-R2ME**2)*(Z+(1D0-Z)*PM1ME**2/PS1ME)
1840 X2=1D0+R2ME**2-PS1ME/ECMME**2
1841C...Call ME, with right order important for two inequivalent showerers.
1842 IF(IR.EQ.IORD+30) THEN
1843 WME=PYMAEL(M3JCC,X1,X2,R1ME,R2ME,ALPHA)
1844 ELSE
1845 WME=PYMAEL(M3JCC,X2,X1,R2ME,R1ME,ALPHA)
1846 ENDIF
1847C...Split up total ME when two radiating partons.
1848 ISPRAD=1
1849 IF((M3JCC.GE.16.AND.M3JCC.LE.19).OR.
1850 & (M3JCC.GE.26.AND.M3JCC.LE.29).OR.
1851 & (M3JCC.GE.36.AND.M3JCC.LE.39).OR.
1852 & (M3JCC.GE.46.AND.M3JCC.LE.49).OR.
1853 & (M3JCC.GE.56.AND.M3JCC.LE.64)) ISPRAD=0
1854 IF(ISPRAD.EQ.1) WME=WME*MAX(1D-10,1D0+R1ME**2-R2ME**2-X1)/
1855 & MAX(1D-10,2D0-X1-X2)
1856C...Evaluate shower rate to be compared with.
1857 WSHOW=2D0/(MAX(1D-10,2D0-X1-X2)*
1858 & MAX(1D-10,1D0+R2ME**2-R1ME**2-X2))
1859 IF(IGLUI.EQ.1.AND.IR.GE.31) WSHOW=(9D0/4D0)*WSHOW
1860 ELSEIF(MSTJ(49).NE.1) THEN
1861
1862C...Toy model scalar theory matrix elements; no mass effects.
1863 ELSE
1864 X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5))
1865 X2=1D0-V(IEP(1),5)/V(NS+1,5)
1866 X3=(1D0-X1)+(1D0-X2)
1867 WSHOW=4D0*X3*((1D0-X1)/(2D0-X2)**2+(1D0-X2)/(2D0-X1)**2)
1868 WME=X3**2
1869 IF(MSTJ(102).GE.2) WME=X3**2-2D0*(1D0+X3)*(1D0-X1)*(1D0-X2)*
1870 & PARJ(171)
1871 ENDIF
1872
1873 IF(WME.LT.PYR(0)*WSHOW) GOTO 410
1874 ENDIF
1875
1876C...Impose angular ordering by rejection of nonordered emission.
1877 IF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2.AND.IPSPD.EQ.0) THEN
1878 PEMAO=V(IM,1)*P(IM,4)
1879 IF(IEP(1).EQ.N+2) PEMAO=(1D0-V(IM,1))*P(IM,4)
1880 IF(IR.GE.31.AND.MSTJ(42).GE.5) THEN
1881 MAOD=0
1882 ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.4
1883 & .OR.MSTJ(42).EQ.7)) THEN
1884 MAOD=0
1885 ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.3
1886 & .OR.MSTJ(42).EQ.6)) THEN
1887 MAOD=1
1888 PMDAO=PMTH(2,K(IEP(1),5))
1889 THE2ID=Z*(1D0-Z)*PEMAO**2/(V(IEP(1),5)-4D0*PMDAO**2)
1890 ELSE
1891 MAOD=1
1892 THE2ID=Z*(1D0-Z)*PEMAO**2/V(IEP(1),5)
1893 IF(MSTJ(42).GE.3.AND.MSTJ(42).NE.5) THE2ID=THE2ID*
1894 & (1D0+PMTH(1,IR)**2*(1D0-Z)/(V(IEP(1),5)*Z))**2
1895 ENDIF
1896 MAOM=1
1897 IAOM=IM
1898 440 IF(K(IAOM,5).EQ.22) THEN
1899 IAOM=K(IAOM,3)
1900 IF(K(IAOM,3).LE.NS) MAOM=0
1901 IF(MAOM.EQ.1) GOTO 440
1902 ENDIF
1903 IF(MAOM.EQ.1.AND.MAOD.EQ.1) THEN
1904 THE2IM=V(IAOM,1)*(1D0-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
1905 IF(THE2ID.LT.THE2IM) GOTO 410
1906 ENDIF
1907 ENDIF
1908
1909C...Impose user-defined maximum angle at first branching.
1910 IF(MSTJ(48).EQ.1.AND.IPSPD.EQ.0) THEN
1911 IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
1912 THE2ID=Z*(1D0-Z)*PS(4)**2/V(IEP(1),5)
1913 IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 410
1914 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
1915 THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5)
1916 IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 410
1917 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
1918 THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5)
1919 IF(PARJ(86)**2*THE2ID.LT.1D0) GOTO 410
1920 ENDIF
1921 ENDIF
1922
1923C...Impose angular constraint in first branching from interference
1924C...with initial state partons.
1925 IF(MIIS.GE.2.AND.IEP(1).LE.NS+3) THEN
1926 THE2D=MAX((1D0-Z)/Z,Z/(1D0-Z))*V(IEP(1),5)/(0.5D0*P(IM,4))**2
1927 IF(IEP(1).EQ.NS+2.AND.ISII(1).GE.1) THEN
1928 IF(THE2D.GT.THEIIS(1,ISII(1))**2) GOTO 410
1929 ELSEIF(IEP(1).EQ.NS+3.AND.ISII(2).GE.1) THEN
1930 IF(THE2D.GT.THEIIS(2,ISII(2))**2) GOTO 410
1931 ENDIF
1932 ENDIF
1933
1934C...End of inner veto algorithm. Check if only one leg evolved so far.
1935 450 V(IEP(1),1)=Z
1936 ISL(1)=0
1937 ISL(2)=0
1938 IF(NEP.EQ.1) GOTO 490
1939 IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 350
1940 DO 460 I=1,NEP
1941 IR=IREF(N+I-NS)
1942 IF(ITRY(I).EQ.0.AND.KSH(IR).EQ.1) THEN
1943 IF(P(N+I,5).GE.PMTH(2,IR)) GOTO 350
1944 ENDIF
1945 460 CONTINUE
1946
1947C...Check if chosen multiplet m1,m2,z1,z2 is physical.
1948 IF(NEP.GE.3) THEN
1949 PMSUM=0D0
1950 DO 470 I=1,NEP
1951 PMSUM=PMSUM+P(N+I,5)
1952 470 CONTINUE
1953 IF(PMSUM.GE.PS(5)) GOTO 350
1954 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
1955 DO 480 I1=N+1,N+2
1956 IRDA=IREF(I1-NS)
1957 IF(KSH(IRDA).EQ.0) GOTO 480
1958 IF(P(I1,5).LT.PMTH(2,IRDA)) GOTO 480
1959 IF(IRDA.EQ.21) THEN
1960 IRGD1=IABS(K(I1,5))
1961 IRGD2=IRGD1
1962 ELSE
1963 IRGD1=IRDA
1964 IRGD2=IABS(K(I1,5))
1965 ENDIF
1966 I2=2*N+3-I1
1967 IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1968 PED=0.5D0*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
1969 ELSE
1970 IF(I1.EQ.N+1) ZM=V(IM,1)
1971 IF(I1.EQ.N+2) ZM=1D0-V(IM,1)
1972 PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
1973 & 4D0*V(N+1,5)*V(N+2,5))
1974 PED=PEM*(0.5D0*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/
1975 & V(IM,5)
1976 ENDIF
1977 IF(MOD(MSTJ(43),2).EQ.1) THEN
1978 PMQTH3=0.5D0*PARJ(82)
1979 IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83)
1980 IF(IRGD2.EQ.22.AND.ISCOL(IRDA).EQ.0) PMQTH3=0.5D0*PARJ(90)
1981 PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(I1,5)
1982 PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(I1,5)
1983 ZD=SQRT(MAX(0D0,(1D0-V(I1,5)/PED**2)*((1D0-PMQ1-PMQ2)**2-
1984 & 4D0*PMQ1*PMQ2)))
1985 ZH=1D0+PMQ1-PMQ2
1986 ELSE
1987 ZD=SQRT(MAX(0D0,1D0-V(I1,5)/PED**2))
1988 ZH=1D0
1989 ENDIF
1990 IF(IRDA.EQ.21.AND.IRGD1.LT.10.AND.
1991 & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1992 ELSE
1993 ZL=0.5D0*(ZH-ZD)
1994 ZU=0.5D0*(ZH+ZD)
1995 IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND.
1996 & ISSET(1).EQ.0) THEN
1997 ISL(1)=1
1998 ELSEIF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND.
1999 & ISSET(2).EQ.0) THEN
2000 ISL(2)=1
2001 ENDIF
2002 ENDIF
2003 IF(IRDA.EQ.21) V(I1,4)=LOG(ZU*(1D0-ZL)/MAX(1D-20,
2004 & ZL*(1D0-ZU)))
2005 IF(IRDA.NE.21) V(I1,4)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
2006 480 CONTINUE
2007 IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
2008 ISL(3-ISLM)=0
2009 ISLM=3-ISLM
2010 ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
2011 ZDR1=MAX(0D0,V(N+1,3)/MAX(1D-6,V(N+1,4))-1D0)
2012 ZDR2=MAX(0D0,V(N+2,3)/MAX(1D-6,V(N+2,4))-1D0)
2013 IF(ZDR2.GT.PYR(0)*(ZDR1+ZDR2)) ISL(1)=0
2014 IF(ISL(1).EQ.1) ISL(2)=0
2015 IF(ISL(1).EQ.0) ISLM=1
2016 IF(ISL(2).EQ.0) ISLM=2
2017 ENDIF
2018 IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 350
2019 ENDIF
2020 IRD1=IREF(N+1-NS)
2021 IRD2=IREF(N+2-NS)
2022 IF(IGM.GT.0) THEN
2023 IF(MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
2024 & PMTH(2,IRD1).OR.P(N+2,5).GE.PMTH(2,IRD2))) THEN
2025 PMQ1=V(N+1,5)/V(IM,5)
2026 PMQ2=V(N+2,5)/V(IM,5)
2027 ZD=SQRT(MAX(0D0,(1D0-V(IM,5)/PEM**2)*((1D0-PMQ1-PMQ2)**2-
2028 & 4D0*PMQ1*PMQ2)))
2029 ZH=1D0+PMQ1-PMQ2
2030 ZL=0.5D0*(ZH-ZD)
2031 ZU=0.5D0*(ZH+ZD)
2032 IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 350
2033 ENDIF
2034 ENDIF
2035
2036C...Accepted branch. Construct four-momentum for initial partons.
2037 490 MAZIP=0
2038 MAZIC=0
2039 IF(NEP.EQ.1) THEN
2040 P(N+1,1)=0D0
2041 P(N+1,2)=0D0
2042 P(N+1,3)=SQRT(MAX(0D0,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
2043 & P(N+1,5))))
2044 P(N+1,4)=P(IPA(1),4)
2045 V(N+1,2)=P(N+1,4)
2046 ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
2047 PED1=0.5D0*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
2048 P(N+1,1)=0D0
2049 P(N+1,2)=0D0
2050 P(N+1,3)=SQRT(MAX(0D0,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
2051 P(N+1,4)=PED1
2052 P(N+2,1)=0D0
2053 P(N+2,2)=0D0
2054 P(N+2,3)=-P(N+1,3)
2055 P(N+2,4)=P(IM,5)-PED1
2056 V(N+1,2)=P(N+1,4)
2057 V(N+2,2)=P(N+2,4)
2058 ELSEIF(NEP.GE.3) THEN
2059C...Rescale all momenta for energy conservation.
2060 LOOP=0
2061 PES=0D0
2062 PQS=0D0
2063 DO 510 I=1,NEP
2064 DO 500 J=1,4
2065 P(N+I,J)=P(IPA(I),J)
2066 500 CONTINUE
2067 PES=PES+P(N+I,4)
2068 PQS=PQS+P(N+I,5)**2/P(N+I,4)
2069 510 CONTINUE
2070 520 LOOP=LOOP+1
2071 FAC=(PS(5)-PQS)/(PES-PQS)
2072 PES=0D0
2073 PQS=0D0
2074 DO 540 I=1,NEP
2075 DO 530 J=1,3
2076 P(N+I,J)=FAC*P(N+I,J)
2077 530 CONTINUE
2078 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)
2079 V(N+I,2)=P(N+I,4)
2080 PES=PES+P(N+I,4)
2081 PQS=PQS+P(N+I,5)**2/P(N+I,4)
2082 540 CONTINUE
2083 IF(LOOP.LT.10.AND.ABS(PES-PS(5)).GT.1D-12*PS(5)) GOTO 520
2084
2085C...Construct transverse momentum for ordinary branching in shower.
2086 ELSE
2087 ZM=V(IM,1)
2088 LOOPPT=0
2089 550 LOOPPT=LOOPPT+1
2090 PZM=SQRT(MAX(0D0,(PEM+P(IM,5))*(PEM-P(IM,5))))
2091 PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4D0*V(N+1,5)*V(N+2,5)
2092 IF(PZM.LE.0D0) THEN
2093 PTS=0D0
2094 ELSEIF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
2095 & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
2096 PTS=PMLS*ZM*(1D0-ZM)/V(IM,5)
2097 ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
2098 PTS=(PEM**2*(ZM*(1D0-ZM)*V(IM,5)-(1D0-ZM)*V(N+1,5)-
2099 & ZM*V(N+2,5))-0.25D0*PMLS)/PZM**2
2100 ELSE
2101 PTS=PMLS*(ZM*(1D0-ZM)*PEM**2/V(IM,5)-0.25D0)/PZM**2
2102 ENDIF
2103 IF(PTS.LT.0D0.AND.LOOPPT.LT.10) THEN
2104 ZM=0.05D0+0.9D0*ZM
2105 GOTO 550
2106 ELSEIF(PTS.LT.0D0) THEN
2107 GOTO 280
2108 ENDIF
2109 PT=SQRT(MAX(0D0,PTS))
2110
2111C...Global statistics.
2112 MINT(353)=MINT(353)+1
2113 VINT(353)=VINT(353)+PT
2114 IF (MINT(353).EQ.1) VINT(358)=PT
2115
2116C...Find coefficient of azimuthal asymmetry due to gluon polarization.
2117 HAZIP=0D0
2118 IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21
2119 & .AND.IAU.NE.0) THEN
2120 IF(K(IGM,3).NE.0) MAZIP=1
2121 ZAU=V(IGM,1)
2122 IF(IAU.EQ.IM+1) ZAU=1D0-V(IGM,1)
2123 IF(MAZIP.EQ.0) ZAU=0D0
2124 IF(K(IGM,2).NE.21) THEN
2125 HAZIP=2D0*ZAU/(1D0+ZAU**2)
2126 ELSE
2127 HAZIP=(ZAU/(1D0-ZAU*(1D0-ZAU)))**2
2128 ENDIF
2129 IF(K(N+1,2).NE.21) THEN
2130 HAZIP=HAZIP*(-2D0*ZM*(1D0-ZM))/(1D0-2D0*ZM*(1D0-ZM))
2131 ELSE
2132 HAZIP=HAZIP*(ZM*(1D0-ZM)/(1D0-ZM*(1D0-ZM)))**2
2133 ENDIF
2134 ENDIF
2135
2136C...Find coefficient of azimuthal asymmetry due to soft gluon
2137C...interference.
2138 HAZIC=0D0
2139 IF(MSTJ(49).NE.2.AND.MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.
2140 & K(N+2,2).EQ.21).AND.IAU.NE.0) THEN
2141 IF(K(IGM,3).NE.0) MAZIC=N+1
2142 IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
2143 IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
2144 & ZM.GT.0.5D0) MAZIC=N+2
2145 IF(K(IAU,2).EQ.22) MAZIC=0
2146 ZS=ZM
2147 IF(MAZIC.EQ.N+2) ZS=1D0-ZM
2148 ZGM=V(IGM,1)
2149 IF(IAU.EQ.IM-1) ZGM=1D0-V(IGM,1)
2150 IF(MAZIC.EQ.0) ZGM=1D0
2151 IF(MAZIC.NE.0) HAZIC=(P(IM,5)/P(IGM,5))*
2152 & SQRT((1D0-ZS)*(1D0-ZGM)/(ZS*ZGM))
2153 HAZIC=MIN(0.95D0,HAZIC)
2154 ENDIF
2155 ENDIF
2156
2157C...Construct energies for ordinary branching in shower.
2158 560 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
2159 IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
2160 & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
2161 P(N+1,4)=0.5D0*(PEM*(V(IM,5)+V(N+1,5)-V(N+2,5))+
2162 & PZM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5)
2163 ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
2164 P(N+1,4)=PEM*V(IM,1)
2165 ELSE
2166 P(N+1,4)=PEM*(0.5D0*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
2167 & SQRT(PMLS)*ZM)/V(IM,5)
2168 ENDIF
2169
2170C...Already predetermined choice of phi angle or not
2171
2172 PHI=PARU(2)*PYR(0)
2173 IF(MPSPD.EQ.1.AND.IGM.EQ.NS+1) THEN
2174 IPSPD=IP1+IM-NS-2
2175 IF(K(IPSPD,4).GT.0) THEN
2176 IPSGD1=K(IPSPD,4)
2177 IF(IM.EQ.NS+2) THEN
2178 PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
2179 ELSE
2180 PHI=PYANGL(-P(IPSGD1,1),P(IPSGD1,2))
2181 ENDIF
2182 ENDIF
2183 ELSEIF(MPSPD.EQ.1.AND.IGM.EQ.NS+2) THEN
2184 IPSPD=IP1+IM-NS-2
2185 IF(K(IPSPD,4).GT.0) THEN
2186 IPSGD1=K(IPSPD,4)
2187 PHIPSM=PYANGL(P(IPSPD,1),P(IPSPD,2))
2188 THEPSM=PYANGL(P(IPSPD,3),SQRT(P(IPSPD,1)**2+P(IPSPD,2)**2))
2189 CALL PYROBO(IPSGD1,IPSGD1,0D0,-PHIPSM,0D0,0D0,0D0)
2190 CALL PYROBO(IPSGD1,IPSGD1,-THEPSM,0D0,0D0,0D0,0D0)
2191 PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
2192 CALL PYROBO(IPSGD1,IPSGD1,THEPSM,PHIPSM,0D0,0D0,0D0)
2193 ENDIF
2194 ENDIF
2195
2196C...Construct momenta for ordinary branching in shower.
2197 P(N+1,1)=PT*COS(PHI)
2198 P(N+1,2)=PT*SIN(PHI)
2199 IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
2200 & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
2201 P(N+1,3)=0.5D0*(PZM*(V(IM,5)+V(N+1,5)-V(N+2,5))+
2202 & PEM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5)
2203 ELSEIF(PZM.GT.0D0) THEN
2204 P(N+1,3)=0.5D0*(V(N+2,5)-V(N+1,5)-V(IM,5)+
2205 & 2D0*PEM*P(N+1,4))/PZM
2206 ELSE
2207 P(N+1,3)=0D0
2208 ENDIF
2209 P(N+2,1)=-P(N+1,1)
2210 P(N+2,2)=-P(N+1,2)
2211 P(N+2,3)=PZM-P(N+1,3)
2212 P(N+2,4)=PEM-P(N+1,4)
2213 IF(MSTJ(43).LE.2) THEN
2214 V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
2215 V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
2216 ENDIF
2217 ENDIF
2218
2219C...Rotate and boost daughters.
2220 IF(IGM.GT.0) THEN
2221 IF(MSTJ(43).LE.2) THEN
2222 BEX=P(IGM,1)/P(IGM,4)
2223 BEY=P(IGM,2)/P(IGM,4)
2224 BEZ=P(IGM,3)/P(IGM,4)
2225 GA=P(IGM,4)/P(IGM,5)
2226 GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1D0+GA)-
2227 & P(IM,4))
2228 ELSE
2229 BEX=0D0
2230 BEY=0D0
2231 BEZ=0D0
2232 GA=1D0
2233 GABEP=0D0
2234 ENDIF
2235 PTIMB=SQRT((P(IM,1)+GABEP*BEX)**2+(P(IM,2)+GABEP*BEY)**2)
2236 THE=PYANGL(P(IM,3)+GABEP*BEZ,PTIMB)
2237 IF(PTIMB.GT.1D-4) THEN
2238 PHI=PYANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
2239 ELSE
2240 PHI=0D0
2241 ENDIF
2242 DO 570 I=N+1,N+2
2243 DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
2244 & SIN(THE)*COS(PHI)*P(I,3)
2245 DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
2246 & SIN(THE)*SIN(PHI)*P(I,3)
2247 DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
2248 DP(4)=P(I,4)
2249 DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
2250 DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
2251 P(I,1)=DP(1)+DGABP*BEX
2252 P(I,2)=DP(2)+DGABP*BEY
2253 P(I,3)=DP(3)+DGABP*BEZ
2254 P(I,4)=GA*(DP(4)+DBP)
2255 570 CONTINUE
2256 ENDIF
2257
2258C...Weight with azimuthal distribution, if required.
2259 IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
2260 DO 580 J=1,3
2261 DPT(1,J)=P(IM,J)
2262 DPT(2,J)=P(IAU,J)
2263 DPT(3,J)=P(N+1,J)
2264 580 CONTINUE
2265 DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
2266 DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
2267 DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
2268 DO 590 J=1,3
2269 DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/MAX(1D-10,DPMM)
2270 DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/MAX(1D-10,DPMM)
2271 590 CONTINUE
2272 DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
2273 DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
2274 IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1D0*PARJ(82)) THEN
2275 CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
2276 & DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
2277 IF(MAZIP.NE.0) THEN
2278 IF(1D0+HAZIP*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(HAZIP)))
2279 & GOTO 560
2280 ENDIF
2281 IF(MAZIC.NE.0) THEN
2282 IF(MAZIC.EQ.N+2) CAD=-CAD
2283 IF((1D0-HAZIC)*(1D0-HAZIC*CAD)/(1D0+HAZIC**2-2D0*HAZIC*CAD)
2284 & .LT.PYR(0)) GOTO 560
2285 ENDIF
2286 ENDIF
2287 ENDIF
2288
2289C...Azimuthal anisotropy due to interference with initial state partons.
2290 IF(MOD(MIIS,2).EQ.1.AND.IGM.EQ.NS+1.AND.(K(N+1,2).EQ.21.OR.
2291 &K(N+2,2).EQ.21)) THEN
2292 III=IM-NS-1
2293 IF(ISII(III).GE.1) THEN
2294 IAZIID=N+1
2295 IF(K(N+1,2).NE.21) IAZIID=N+2
2296 IF(K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
2297 & P(N+1,4).GT.P(N+2,4)) IAZIID=N+2
2298 THEIID=PYANGL(P(IAZIID,3),SQRT(P(IAZIID,1)**2+P(IAZIID,2)**2))
2299 IF(III.EQ.2) THEIID=PARU(1)-THEIID
2300 PHIIID=PYANGL(P(IAZIID,1),P(IAZIID,2))
2301 HAZII=MIN(0.95D0,THEIID/THEIIS(III,ISII(III)))
2302 CAD=COS(PHIIID-PHIIIS(III,ISII(III)))
2303 PHIREL=ABS(PHIIID-PHIIIS(III,ISII(III)))
2304 IF(PHIREL.GT.PARU(1)) PHIREL=PARU(2)-PHIREL
2305 IF((1D0-HAZII)*(1D0-HAZII*CAD)/(1D0+HAZII**2-2D0*HAZII*CAD)
2306 & .LT.PYR(0)) GOTO 560
2307 ENDIF
2308 ENDIF
2309
2310C...Continue loop over partons that may branch, until none left.
2311 IF(IGM.GE.0) K(IM,1)=14
2312 N=N+NEP
2313 NEP=2
2314 IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
2315 CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
2316 IF(MSTU(21).GE.1) N=NS
2317 IF(MSTU(21).GE.1) RETURN
2318 ENDIF
2319 GOTO 290
2320
2321C...Set information on imagined shower initiator.
2322 600 IF(NPA.GE.2) THEN
2323 K(NS+1,1)=11
2324 K(NS+1,2)=94
2325 K(NS+1,3)=IP1
2326 IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
2327 K(NS+1,4)=NS+2
2328 K(NS+1,5)=NS+1+NPA
2329 IIM=1
2330 ELSE
2331 IIM=0
2332 ENDIF
2333
2334C...Reconstruct string drawing information.
2335 DO 610 I=NS+1+IIM,N
2336 KQ=KCHG(PYCOMP(K(I,2)),2)
2337 IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
2338 K(I,1)=1
2339 ELSEIF(K(I,1).LE.10.AND.IABS(K(I,2)).GE.11.AND.
2340 & IABS(K(I,2)).LE.18) THEN
2341 K(I,1)=1
2342 ELSEIF(K(I,1).LE.10) THEN
2343 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
2344 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
2345 ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
2346 ID1=MOD(K(I,4),MSTU(5))
2347 IF(KQ.EQ.1.AND.K(I,2).GT.0) ID1=MOD(K(I,4),MSTU(5))+1
2348 IF(KQ.EQ.2.AND.(K(ID1,2).EQ.21.OR.K(ID1+1,2).EQ.21).AND.
2349 & PYR(0).GT.0.5D0) ID1=MOD(K(I,4),MSTU(5))+1
2350 ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
2351 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
2352 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
2353 K(ID1,4)=K(ID1,4)+MSTU(5)*I
2354 K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
2355 K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
2356 K(ID2,5)=K(ID2,5)+MSTU(5)*I
2357 ELSE
2358 ID1=MOD(K(I,4),MSTU(5))
2359 ID2=ID1+1
2360 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
2361 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
2362 IF(KQ.EQ.1.OR.K(ID1,1).GE.11) THEN
2363 K(ID1,4)=K(ID1,4)+MSTU(5)*I
2364 K(ID1,5)=K(ID1,5)+MSTU(5)*I
2365 ELSE
2366 K(ID1,4)=0
2367 K(ID1,5)=0
2368 ENDIF
2369 K(ID2,4)=0
2370 K(ID2,5)=0
2371 ENDIF
2372 610 CONTINUE
2373
2374C...Transformation from CM frame.
2375 IF(NPA.EQ.1) THEN
2376 THE=PYANGL(P(IPA(1),3),SQRT(P(IPA(1),1)**2+P(IPA(1),2)**2))
2377 PHI=PYANGL(P(IPA(1),1),P(IPA(1),2))
2378 MSTU(33)=1
2379 CALL PYROBO(NS+1,N,THE,PHI,0D0,0D0,0D0)
2380 ELSEIF(NPA.EQ.2) THEN
2381 BEX=PS(1)/PS(4)
2382 BEY=PS(2)/PS(4)
2383 BEZ=PS(3)/PS(4)
2384 GA=PS(4)/PS(5)
2385 GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
2386 & /(1D0+GA)-P(IPA(1),4))
2387 THE=PYANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
2388 & +GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
2389 PHI=PYANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
2390 MSTU(33)=1
2391 CALL PYROBO(NS+1,N,THE,PHI,BEX,BEY,BEZ)
2392 ELSE
2393 CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),
2394 & PS(3)/PS(4))
2395 MSTU(33)=1
2396 CALL PYROBO(NS+1,N,0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),PS(3)/PS(4))
2397 ENDIF
2398
2399C...Decay vertex of shower.
2400 DO 630 I=NS+1,N
2401 DO 620 J=1,5
2402 V(I,J)=V(IP1,J)
2403 620 CONTINUE
2404 630 CONTINUE
2405
2406C...Delete trivial shower, else connect initiators.
2407 IF(N.LE.NS+NPA+IIM) THEN
2408 N=NS
2409 ELSE
2410 DO 640 IP=1,NPA
2411 K(IPA(IP),1)=14
2412 K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP
2413 K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP
2414 K(NS+IIM+IP,3)=IPA(IP)
2415 IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1
2416 IF(K(NS+IIM+IP,1).NE.1) THEN
2417 K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)
2418 K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)
2419 ENDIF
2420 640 CONTINUE
2421 ENDIF
2422
2423 RETURN
2424 END
2425
2426
2427 SUBROUTINE QPYGIN(X0,Y0,Z0,T0)
2428C NOT TO BE TOUCHED BY THE USER:
2429C IT SETS THE INITIAL POSITION AND TIME OF THE
2430C PARENT BRANCHING PARTON (X, Y, Z, T, IN FM) IN THE CENTER-OF-MASS
2431C FRAME OF THE HARD COLLISION (IF APPLICABLE FOR THE TYPE OF EVENTS
2432C YOU ARE SIMULATING). INFORMATION ABOUT THE BOOST AND ROTATION IS
2433C CONTAINED IN THE IN COMMON QPLT BELOW.
2434 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2435C NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST
2436C PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE
2437C CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD
2438C SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
2439C XX, YY, ZZ, TT COMING FROM QPYGIN0 IN THE CMS OF THE COLLISION, GENERATED
2440C ONCE PER NUCLEON-NUCLEON COLLISION.
2441 COMMON/TOQPYGIN/XX,YY,ZZ,TT
2442 COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ
2443C Boost and rotation.
2444 call qpyrobo(xx,yy,zz,tt,0.d0,0.d0,bbx,bby,bbz,x1,y1,z1,t1)
2445 call qpyrobo(x1,y1,z1,t1,0.d0,aa2,0.d0,0.d0,0.d0,x2,y2,z2,t2)
2446 call qpyrobo(x2,y2,z2,t2,aa1,0.d0,0.d0,0.d0,0.d0,x0,y0,z0,t0)
2447C
2448 RETURN
2449 END
2450
2451
2452
2453 SUBROUTINE QPYGIN0
2454C USER-DEFINED ROUTINE: IT SETS THE INITIAL POSITION AND TIME OF THE
2455C PARENT BRANCHING PARTON (X, Y, Z, T, IN FM) IN THE CENTER-OF-MASS
2456C FRAME OF THE COLLISION.
2457C IT SHOULD BE CALLED IN THE MAIN PROGRAM, ONCE PER NUCLEON-NUCLEON
2458C COLLISION (PER PP COLLISION IF YOU ARE SIMPLY RUNNING PP).
2459 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2460 COMMON/TOQPYGIN/XX,YY,ZZ,TT
2461C Example: this is to be set by the user e.g. as coming from the overlap
2462C of nuclear profile functions (number of binary nucleon-nucleon
2463C collisions).
2464 call GetRandomXY(xrang,yrang)
2465 xin=xrang ! fm
2466 yin=yrang ! fm
2467 call savexy(xin, yin)
2468 zin=0.d0
2469 tin=0.d0
2470C Passing to the common block.
2471 xx=xin
2472 yy=yin
2473 zz=zin
2474 tt=tin
2475c
2476c print*,"initial coordinates", xx,yy,zz,tt
2477 RETURN
2478 END
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489C
2490 SUBROUTINE QPYGEO(x,y,z,t,bx,by,bz,qhl,oc)
2491C USER-DEFINED ROUTINE:
2492C The values of qhatL and omegac have to be computed
2493C by the user, using his preferred medium model, in
2494C this routine, which takes as input the position
2495C x,y,z,t (in fm) of the parton to branch, the trajectory
2496C defined by the three-vector bx,by,bz (in units of c),
2497C (all values in the center-of-mass frame of the
2498C hard collision), and returns the value of qhatL
2499C (in GeV**2) and omegac (in GeV).
2500 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2501C NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST
2502C PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE
2503C CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD
2504C SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
2505 COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ
2506 COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2507 qhat=parj(198)
2508 xl=parj(199)
2509 bimp=parj(197)
2510c Here we give five options for the geometry of the medium:
2511
2512c$$$cforalice+
2513c$$$ (0) we call routine CalculateI0I1 in AliFastGlauber. Given the position
2514c$$$ of the parton in the reaction plane (x,y), the direction in the
2515c$$$ reaction plane phi=atan(py/px) and the impact parameter of the
2516c$$$ collision bimp, it gives back the transverse path length to the
2517c$$$ "end" of the medium and the integrated qhat along that path length.
2518c$$$ See the seter for this option in Configqpythia.C(one can set the
2519c$$$ value of xkscale by doing SetQhat(xkscale), with xkscale in fm).
2520c$$$ The set value is passed here through the pythia free parameter parj(198).
2521
2522
2523 xkscale=parj(198)
2524 ellcut=20.d0
2525 xlcero=0.d0
2526 xlone=0.d0
2527
2528 call qpyrobo(x,y,z,t,-aa1,-aa2,-bbx,-bby,-bbz,xout,
2529 + yout,zout,tout)
2530
2531 call qpyrobo(bx,by,bz,1.d0,-aa1,-aa2,-bbx,-bby,-bbz,
2532 + bx1,by1,bz1,bt1)
2533
2534
2535 phi=datan2(by1,bx1)
2536 phia=phi
2537 bimpa=bimp
2538 ellcuta=ellcut
2539 xa=xout
2540 ya=yout
2541 call CalculateI0I1(xlcero,xlone,bimpa,xa,ya,phia,ellcuta)
2542 if(xlcero.eq.0.d0) then
2543 xlp=0.d0
2544 qhl=0.d0
2545 else
2546 xlp=2.d0*xlone/xlcero
2547 qhl=0.1973d0*0.1973d0*xlcero*xkscale
2548 call savei0i1(xlcero, xlone)
2549 endif
2550
2551
2552c$$$cforalice-
2553c To use any of these folowing 1,2,3 or 4 options the user should specify
2554c a constant value for the transport coefficient and an initial in medium length.
2555c This can be done in the user Config file by setting: SetQhat(qhat), with qhat in
2556c GeV2/fm and SetLength(xl), with xl in fm. Those values are passed here through
2557c pythia free parameters parj(198) and parj(199).
2558
2559c (1) to fix the length to the initial value, uncomment the next three lines
2560c and comment the other definitions of xlp and qhl above and below.
2561c xlp=xl
2562c if (xlp .lt. 0.d0) xlp=0.d0
2563c qhl=xlp*qhat ! GeV**2
2564c print*, xlp,qhl
2565c (2) simplest ansatz: for an initial parton along the z-axis (approximate)
2566c starting in the center of a medium (-xl,+xl) along the z-axis
2567c if (bz .gt. 0.d0) then
2568c xlp=xl-z
2569c else
2570c xlp=xl+z
2571c endif
2572c if (xlp .gt. (2.d0*xl)) xlp=2.d0*xl
2573c if (xlp .lt. 0.d0) xlp=0.d0
2574c qhl=xlp*qhat ! GeV**2
2575
2576c (3) for a parton at midrapidity inside a cylinder (approximate)
2577c xlp=xl-dsqrt(x*x+y*y)
2578c if (xlp .lt. 0.d0) xlp=0.d0
2579c qhl=xlp*qhat ! GeV**2
2580
2581c (4) for a brick defined by planes (x,y,0) and (x,y,xl), comment
2582c the previous lines and uncomment lines between the comment 'brick'.
2583c brick+
2584c if (z .ge. 0.d0 .and. z .le. xl)
2585c > then
2586c if (bz .gt. 0.d0) then
2587c ttpp=(xl-z)/bz
2588c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2589c > (xl-z)**2.d0)
2590c else
2591c ttpp=z/dabs(bz)
2592c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2593c > (z)**2.d0)
2594c endif
2595c elseif (z .lt. 0.d0) then
2596c if (bz .lt. 0.d0) then
2597c xlp=0.d0
2598c else
2599c ttpp1=-z/bz
2600c ttpp2=(xl-z)/bz
2601c xxpp1=x+bx*ttpp1
2602c xxpp2=x+bx*ttpp2
2603c yypp1=y+by*ttpp1
2604c yypp2=y+by*ttpp2
2605c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2606c > xl**2.d0)
2607c endif
2608c elseif (z .gt. xl) then
2609c if (bz .gt. 0.d0) then
2610c xlp=0.d0
2611c else
2612c ttpp1=z/dabs(bz)
2613c ttpp2=(-xl+z)/dabs(bz)
2614c xxpp1=x+bx*ttpp1
2615c xxpp2=x+bx*ttpp2
2616c yypp1=y+by*ttpp1
2617c yypp2=y+by*ttpp2
2618c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2619c > xl**2.d0)
2620c endif
2621c endif
2622c if (xlp .lt. 0.d0) xlp=0.d0
2623c qhl=xlp*qhat ! GeV**2
2624c brick-
2625
2626 oc=0.5d0*qhl*xlp/0.1973d0 ! GeV
2627 RETURN
2628 END