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