]>
Commit | Line | Data |
---|---|---|
6c43eccb | 1 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
2 | C | |
3 | C AUXILIARY ROUTINES FOR Q-PYTHIA version 1.0. | |
4 | C | |
5 | C DATE: 26.09.2008. | |
6 | C | |
7 | C AUTHORS: N. Armesto, L. Cunqueiro and C. A. Salgado | |
8 | C Departamento de Fisica de Particulas and IGFAE | |
9 | C Universidade de Santiago de Compostela | |
10 | C 15706 Santiago de Compostela, Spain | |
11 | C | |
12 | C EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es, | |
13 | C Carlos.Salgado@cern.ch | |
14 | C | |
15 | C CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18. | |
16 | C NOT to be modified by user. | |
17 | C | |
18 | C WHEN USING Q-PYTHIA, PLEASE QUOTE: | |
19 | C | |
20 | C 1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado, | |
21 | C in preparation. | |
22 | C 2) T. Sjostrand, S. Mrenna and P. Skands, | |
23 | C ``PYTHIA 6.4 physics and manual,'' | |
24 | C JHEP 0605 (2006) 026 [arXiv:hep-ph/0603175]. | |
25 | C | |
26 | C DISCLAIMER: this program comes without any guarantees. Beware of | |
27 | C errors and use common sense when interpreting results. | |
28 | C Any modifications are done under exclusive | |
29 | C maker's resposibility. | |
30 | C | |
31 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC | |
32 | C | |
33 | c VERSION WITH THE LARGE X BEHAVIOR OF THE MEDIUM PART INTRODUCED | |
34 | C BY MULTIPLYING BY THE NUMERATOR OF THE COLLINEAR PART OF THE | |
35 | C VACUUM SPLITTING FUNCTION | |
36 | c | |
37 | function splitq1(w) | |
38 | c to integrate, adding vaccum plus medium q -> qg | |
39 | implicit double precision (a-h,o-z) | |
40 | z=w | |
41 | auxz=z*(1.d0-z) | |
42 | auxq=splitgq(z)+splitmedq1(z) | |
43 | if (auxq .gt. 0.d0) then | |
44 | splitq1=auxq | |
45 | else | |
46 | splitq1=0.d0 | |
47 | endif | |
48 | return | |
49 | end | |
50 | c | |
51 | function splitg1(w) | |
52 | c to integrate, adding vaccum plus medium g -> gg, and g -> qqbar | |
53 | implicit double precision (a-h,o-z) | |
54 | z=w | |
55 | auxz=z*(1.d0-z) | |
56 | c argument of running coupling is taken as kt of emission | |
57 | auxg=splitgg(z)+splitmedg1(z) | |
58 | if (auxg .gt. 0.d0) then | |
59 | splitg1=(auxg+splitqqbar(z)) | |
60 | else | |
61 | splitg1=splitqqbar(z) | |
62 | endif | |
63 | return | |
64 | end | |
65 | c | |
66 | function splitq2(w) | |
67 | c to integrate, adding vaccum plus medium q -> qg | |
68 | implicit double precision (a-h,o-z) | |
69 | z=w | |
70 | auxq=splitgq(z)+splitmedq2(z) | |
71 | if (auxq .gt. 0.d0) then | |
72 | splitq2=auxq | |
73 | else | |
74 | splitq2=0.d0 | |
75 | endif | |
76 | if(splitmedq2(z).eq.0.) then | |
77 | endif | |
78 | return | |
79 | end | |
80 | c | |
81 | function splitg2(z) | |
82 | c to integrate, adding vaccum plus medium g -> gg, and g -> qqbar | |
83 | implicit double precision (a-h,o-z) | |
84 | auxg=splitgg(z)+splitmedg2(z) | |
85 | if(auxg.gt.0.d0) then | |
86 | splitg2=auxg | |
87 | else | |
88 | splitg2=0.d0 | |
89 | endif | |
90 | if(splitmedg2(z).eq.0.) then | |
91 | endif | |
92 | return | |
93 | end | |
94 | c | |
95 | function splitgq(z) | |
96 | c q -> qg splitting kernel at 1 loop for the vacuum | |
97 | implicit double precision (a-h,o-z) | |
98 | xnc=3.d0 | |
99 | splitgq=(0.5d0*(xnc-1.d0/xnc))*(1.d0+z*z)/(1.d0-z) | |
100 | return | |
101 | end | |
102 | c | |
103 | function splitgg(z) | |
104 | c g -> gg splitting kernel at 1 loop for the vacuum | |
105 | implicit double precision (a-h,o-z) | |
106 | xnc=3.d0 | |
107 | auxz=z*(1.d0-z) | |
108 | auxz2=1.d0-auxz | |
109 | splitgg=xnc*auxz2*auxz2/auxz | |
110 | return | |
111 | end | |
112 | c | |
113 | function splitqqbar(z) | |
114 | c g -> qqbar splitting kernel at 1 loop | |
115 | implicit double precision (a-h,o-z) | |
116 | xnf=5.d0 | |
117 | auxz=1.d0-z | |
118 | splitqqbar=0.5d0*xnf*(z*z+auxz*auxz) | |
119 | return | |
120 | end | |
121 | c | |
122 | function splitmedg1(z) | |
123 | c g -> gg splitting kernel at 1 loop for the medium | |
124 | implicit double precision (a-h,o-z) | |
125 | common/qpc1/eee,qhatl,omegac | |
126 | common/qpvir1/pmed | |
127 | xnc=3.d0 | |
128 | pi=dacos(-1.d0) | |
129 | if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then | |
130 | splitmedg1=0.d0 | |
131 | else | |
132 | c symmetrized by hand with respect to 1/2 | |
133 | if (z .ge. 0.5d0) then | |
134 | zz=z | |
135 | else | |
136 | zz=1.d0-z | |
137 | endif | |
138 | t=pmed*pmed | |
139 | auxz=1.d0-zz | |
140 | auxz2=zz*auxz | |
141 | ome=eee*auxz/omegac | |
142 | xkappa2=auxz2*t/qhatl | |
143 | fff=genspec(ome,xkappa2) | |
144 | cc 1/2 to avoid double counting | |
145 | c splitmedg=0.5d0*xnc*2.d0*pi*zz*t*fff/qhatl | |
146 | c we multiply by max(z,1-z) to introduce the large z behavior from the | |
147 | c numerator in the vacuum | |
148 | flx=max(zz,auxz) | |
149 | c 1/2 to avoid double counting | |
150 | splitmedg1=0.5*flx*xnc*2.d0*pi*zz*t*fff/qhatl | |
151 | endif | |
152 | return | |
153 | end | |
154 | c | |
155 | function splitmedq1(z) | |
156 | c q -> qg splitting kernel at 1 loop for the medium | |
157 | implicit double precision (a-h,o-z) | |
158 | common/qpc1/eee,qhatl,omegac | |
159 | common/qpvir1/pmed | |
160 | xnc=3.d0 | |
161 | pi=dacos(-1.d0) | |
162 | if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then | |
163 | splitmedq1=0.d0 | |
164 | else | |
165 | t=pmed*pmed | |
166 | auxz=1.d0-z | |
167 | auxz2=z*auxz | |
168 | ome=eee*auxz/omegac | |
169 | xkappa2=auxz2*t/qhatl | |
170 | fff=genspec(ome,xkappa2) | |
171 | c splitmedq=(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl | |
172 | c we multiply by 1+z**2 to introduce the large z behavior from the | |
173 | c numerator in the vacuum | |
174 | flx=0.5d0*(1.d0+z*z) | |
175 | splitmedq1=flx*(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl | |
176 | endif | |
177 | return | |
178 | end | |
179 | c | |
180 | function splitmedg2(z) | |
181 | c g -> gg splitting kernel at 1 loop for the medium | |
182 | implicit double precision (a-h,o-z) | |
183 | common/qpc1/eee,qhatl,omegac | |
184 | common/qpvir2/virt | |
185 | xnc=3.d0 | |
186 | pi=dacos(-1.d0) | |
187 | if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then | |
188 | splitmedg2=0.d0 | |
189 | else | |
190 | c symmetrized by hand with respect to 1/2 | |
191 | if (z .ge. 0.5d0) then | |
192 | zz=z | |
193 | else | |
194 | zz=1.d0-z | |
195 | endif | |
196 | t=virt | |
197 | auxz=1.d0-zz | |
198 | auxz2=zz*auxz | |
199 | ome=eee*auxz/omegac | |
200 | xkappa2=auxz2*t/qhatl | |
201 | fff=genspec(ome,xkappa2) | |
202 | cc 1/2 to avoid double counting | |
203 | c splitmedg=0.5d0*xnc*2.d0*pi*zz*t*fff/qhatl | |
204 | c we multiply by max(z,1-z) to introduce the large z behavior from the | |
205 | c numerator in the vacuum | |
206 | flx=max(zz,auxz) | |
207 | c 1/2 to avoid double counting | |
208 | splitmedg2=0.5*flx*xnc*2.d0*pi*zz*t*fff/qhatl | |
209 | endif | |
210 | return | |
211 | end | |
212 | c | |
213 | function splitmedq2(z) | |
214 | c q -> qg splitting kernel at 1 loop for the medium | |
215 | implicit double precision (a-h,o-z) | |
216 | common/qpc1/eee,qhatl,omegac | |
217 | common/qpvir2/virt | |
218 | xnc=3.d0 | |
219 | pi=dacos(-1.d0) | |
220 | if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then | |
221 | splitmedq2=0.d0 | |
222 | else | |
223 | t=virt | |
224 | auxz=1.d0-z | |
225 | auxz2=z*auxz | |
226 | ome=eee*auxz/omegac | |
227 | xkappa2=auxz2*t/qhatl | |
228 | fff=genspec(ome,xkappa2) | |
229 | c splitmedq=(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl | |
230 | c we multiply by 1+z**2 to introduce the large z behavior from the | |
231 | c numerator in the vacuum | |
232 | flx=0.5d0*(1.d0+z*z) | |
233 | splitmedq2=flx*(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl | |
234 | endif | |
235 | return | |
236 | end | |
237 | c | |
238 | function genspec(ome,xk2) | |
239 | C THIS FUNCTION GENERATES (omega/omegac) dI/d(omega/omegac) dkappa2, | |
240 | C omegac=qhat L/2, kappa2=kt2/(qhat L), in the mss approximation for m=0, | |
241 | c using interpolation and extrapolation. It reads file grid-qp.dat. | |
242 | c ome=omega/omegac, xk2=kappa2. | |
243 | C MAXIMUM GRID 101 TIMES 101, MODIFY ARRAY DIMENSIONS IF EXCEEDED. | |
244 | c alphas=1, cr=1. | |
245 | implicit double precision (a-h,o-z) | |
246 | dimension xkap2(101), xlkap2(101), xome(101), xlome(101) | |
247 | dimension xspec(101,101) | |
248 | dimension aux1(101), aux2(101) | |
249 | save xkap2, xlkap2, xome, xlome, xspec, npkap, npome | |
250 | character*1000 filnam | |
251 | character*1000 chroot | |
252 | DATA IFLAG/0/ | |
253 | c WE READ THE GRID ONLY THE FIRST TIME. | |
254 | IF (IFLAG .EQ. 0) THEN | |
255 | chroot=' ' | |
256 | call getenvf('ALICE_ROOT',chroot) | |
257 | lnroot= lnblnk(chroot) | |
258 | filnam=chroot(1:lnroot)//'/PYTHIA6/QPYTHIA/qgrid' | |
259 | write(6,*) "Opening file ", filnam | |
260 | open(11,file=filnam, status='old') | |
261 | read(11,*) npkap | |
262 | read(11,*) npome | |
263 | npkap=npkap+1 | |
264 | npome=npome+1 | |
265 | do 10 i=1, npkap, 1 | |
266 | read(11,*) xkap2(i), xlkap2(i) | |
267 | 10 continue | |
268 | do 20 i=1, npome, 1 | |
269 | read(11,*) xome(i), xlome(i) | |
270 | 20 continue | |
271 | do 30 j=1, npome, 1 | |
272 | do 40 i=1, npkap, 1 | |
273 | read(11,*) xspec(i,j) | |
274 | 40 continue | |
275 | 30 continue | |
276 | close(11) | |
277 | iflag=1 | |
278 | ENDIF | |
279 | c cases | |
280 | c for ome>largest value set to 0, | |
281 | c for xk2< smallest value frozen, | |
282 | c for xk2> largest value 1/kappa4 extrapolation. | |
283 | if (ome .gt. xome(npome)) then | |
284 | genspec=0.d0 | |
285 | elseif (ome .lt. xome(1)) then | |
286 | scal=.05648d0*dexp(1.674d0*ome)*dlog(.136d0/ome)/(ome**.5397d0) | |
287 | scal=0.25d0*9.d0*scal/xspec(1,1) | |
288 | if (xk2 .le. xkap2(1)) then | |
289 | genspec=scal*xspec(1,1) | |
290 | elseif (xk2 .eq. xkap2(npkap)) then | |
291 | genspec=scal*xspec(npkap,1) | |
292 | elseif (xk2 .gt. xkap2(npkap)) then | |
293 | genspec=scal*xspec(npkap,1)* | |
294 | > xkap2(npkap)*xkap2(npkap)/(xk2*xk2) | |
295 | else | |
296 | do 50 i=1, npkap, 1 | |
297 | aux1(i)=xspec(i,1) | |
298 | 50 continue | |
299 | genspec=scal*ddivdif(aux1,xlkap2,npkap,dlog(xk2),4) | |
300 | endif | |
301 | else | |
302 | iexact=-1 | |
303 | if (ome .eq. xome(1)) then | |
304 | iexact=1 | |
305 | goto 70 | |
306 | else | |
307 | do 60 i=1, npome-1, 1 | |
308 | if (ome .eq. xome(i+1)) then | |
309 | iexact=i+1 | |
310 | goto 70 | |
311 | elseif (ome .lt. xome(i+1)) then | |
312 | iprev=i | |
313 | ipost=i+1 | |
314 | goto 70 | |
315 | endif | |
316 | 60 continue | |
317 | 70 continue | |
318 | endif | |
319 | if (iexact .gt. 0) then | |
320 | if (xk2 .le. xkap2(1)) then | |
321 | genspec=xspec(1,iexact) | |
322 | elseif (xk2 .eq. xkap2(npkap)) then | |
323 | genspec=xspec(npkap,iexact) | |
324 | elseif (xk2 .gt. xkap2(npkap)) then | |
325 | genspec=xspec(npkap,iexact)* | |
326 | > xkap2(npkap)*xkap2(npkap)/(xk2*xk2) | |
327 | else | |
328 | do 80 i=1, npkap, 1 | |
329 | aux1(i)=xspec(i,iexact) | |
330 | 80 continue | |
331 | genspec=ddivdif(aux1,xlkap2,npkap,dlog(xk2),4) | |
332 | endif | |
333 | else | |
334 | if (xk2 .le. xkap2(1)) then | |
335 | genprev=xspec(1,iprev) | |
336 | genpost=xspec(1,ipost) | |
337 | elseif (xk2 .eq. xkap2(npkap)) then | |
338 | genprev=xspec(npkap,iprev) | |
339 | genpost=xspec(npkap,ipost) | |
340 | elseif (xk2 .gt. xkap2(npkap)) then | |
341 | genprev=xspec(npkap,iprev)* | |
342 | > xkap2(npkap)*xkap2(npkap)/(xk2*xk2) | |
343 | genpost=xspec(npkap,ipost)* | |
344 | > xkap2(npkap)*xkap2(npkap)/(xk2*xk2) | |
345 | else | |
346 | do 90 i=1, npkap, 1 | |
347 | aux1(i)=xspec(i,iprev) | |
348 | aux2(i)=xspec(i,ipost) | |
349 | 90 continue | |
350 | genprev=ddivdif(aux1,xlkap2,npkap,dlog(xk2),4) | |
351 | genpost=ddivdif(aux2,xlkap2,npkap,dlog(xk2),4) | |
352 | endif | |
353 | g12=genprev-genpost | |
354 | xl12=xlome(iprev)-xlome(ipost) | |
355 | c1=g12/xl12 | |
356 | c2=genprev-c1*xlome(iprev) | |
357 | genspec=c1*dlog(ome)+c2 | |
358 | endif | |
359 | endif | |
360 | c | |
361 | RETURN | |
362 | END | |
363 | C | |
364 | * | |
365 | * $Id: divdif.F,v 1.1.1.1 1996/02/15 17:48:36 mclareni Exp $ | |
366 | * | |
367 | * $Log: divdif.F,v $ | |
368 | * Revision 1.1.1.1 1996/02/15 17:48:36 mclareni | |
369 | * Kernlib | |
370 | * | |
371 | * | |
372 | FUNCTION DDIVDIF(F,A,NN,X,MM) | |
373 | c copy of cernlib divdif in double precision. | |
374 | implicit double precision (a-h,o-z) | |
375 | DIMENSION A(NN),F(NN),T(20),D(20) | |
376 | LOGICAL EXTRA | |
377 | LOGICAL MFLAG,RFLAG | |
378 | DATA MMAX/10/ | |
379 | C | |
380 | C TABULAR INTERPOLATION USING SYMMETRICALLY PLACED ARGUMENT POINTS. | |
381 | C | |
382 | C START. FIND SUBSCRIPT IX OF X IN ARRAY A. | |
383 | IF( (NN.LT.2) .OR. (MM.LT.1) ) GO TO 601 | |
384 | N=NN | |
385 | M=MIN0(MM,MMAX,N-1) | |
386 | MPLUS=M+1 | |
387 | IX=0 | |
388 | IY=N+1 | |
389 | IF(A(1).GT.A(N)) GO TO 4 | |
390 | C (SEARCH INCREASING ARGUMENTS.) | |
391 | 1 MID=(IX+IY)/2 | |
392 | IF(X.GE.A(MID)) GO TO 2 | |
393 | IY=MID | |
394 | GO TO 3 | |
395 | C (IF TRUE.) | |
396 | 2 IX=MID | |
397 | 3 IF(IY-IX.GT.1) GO TO 1 | |
398 | GO TO 7 | |
399 | C (SEARCH DECREASING ARGUMENTS.) | |
400 | 4 MID=(IX+IY)/2 | |
401 | IF(X.LE.A(MID)) GO TO 5 | |
402 | IY=MID | |
403 | GO TO 6 | |
404 | C (IF TRUE.) | |
405 | 5 IX=MID | |
406 | 6 IF(IY-IX.GT.1) GO TO 4 | |
407 | C | |
408 | C COPY REORDERED INTERPOLATION POINTS INTO (T(I),D(I)), SETTING | |
409 | C *EXTRA* TO TRUE IF M+2 POINTS TO BE USED. | |
410 | 7 NPTS=M+2-MOD(M,2) | |
411 | IP=0 | |
412 | L=0 | |
413 | GO TO 9 | |
414 | 8 L=-L | |
415 | IF(L.GE.0) L=L+1 | |
416 | 9 ISUB=IX+L | |
417 | IF((1.LE.ISUB).AND.(ISUB.LE.N)) GO TO 501 | |
418 | C (SKIP POINT.) | |
419 | NPTS=MPLUS | |
420 | GO TO 11 | |
421 | C (INSERT POINT.) | |
422 | 501 IP=IP+1 | |
423 | T(IP)=A(ISUB) | |
424 | D(IP)=F(ISUB) | |
425 | 11 IF(IP.LT.NPTS) GO TO 8 | |
426 | EXTRA=NPTS.NE.MPLUS | |
427 | C | |
428 | C REPLACE D BY THE LEADING DIAGONAL OF A DIVIDED-DIFFERENCE TABLE, SUP- | |
429 | C PLEMENTED BY AN EXTRA LINE IF *EXTRA* IS TRUE. | |
430 | DO 14 L=1,M | |
431 | IF(.NOT.EXTRA) GO TO 12 | |
432 | ISUB=MPLUS-L | |
433 | D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB)) | |
434 | 12 I=MPLUS | |
435 | DO 13 J=L,M | |
436 | ISUB=I-L | |
437 | D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB)) | |
438 | I=I-1 | |
439 | 13 CONTINUE | |
440 | 14 CONTINUE | |
441 | C | |
442 | C EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES | |
443 | C OF LAST DIFFERENCE IF *EXTRA* IS TRUE. | |
444 | SUM=D(MPLUS) | |
445 | IF(EXTRA) SUM=0.5*(SUM+D(M+2)) | |
446 | J=M | |
447 | DO 15 L=1,M | |
448 | SUM=D(J)+(X-T(J))*SUM | |
449 | J=J-1 | |
450 | 15 CONTINUE | |
451 | DDIVDIF=SUM | |
452 | RETURN | |
453 | C | |
454 | 601 CALL KERMTR('E105.1',LGFILE,MFLAG,RFLAG) | |
455 | DDIVDIF=0 | |
456 | IF(MFLAG) THEN | |
457 | IF(LGFILE.EQ.0) THEN | |
458 | IF(MM.LT.1) WRITE(*,101) MM | |
459 | IF(NN.LT.2) WRITE(*,102) NN | |
460 | ELSE | |
461 | IF(MM.LT.1) WRITE(LGFILE,101) MM | |
462 | IF(NN.LT.2) WRITE(LGFILE,102) NN | |
463 | ENDIF | |
464 | ENDIF | |
465 | IF(.NOT.RFLAG) CALL ABEND | |
466 | RETURN | |
467 | 101 FORMAT( 7X, 'FUNCTION DDIVDIF ... M =',I6,' IS LESS THAN 1') | |
468 | 102 FORMAT( 7X, 'FUNCTION DDIVDIF ... N =',I6,' IS LESS THAN 2') | |
469 | END | |
470 | c | |
471 | C COPY OF CERN DGAUSS | |
472 | C | |
473 | FUNCTION DGAUSS1(F,A,B,EPS) | |
474 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
475 | DIMENSION W(12),X(12) | |
476 | PARAMETER (Z1 = 1.D0, HF = Z1/2.D0, CST = 5.D0*Z1/1000.D0) | |
477 | DATA X | |
478 | 1 /0.96028 98564 97536 23168 35608 68569 47D0, | |
479 | 2 0.79666 64774 13626 73959 15539 36475 83D0, | |
480 | 3 0.52553 24099 16328 98581 77390 49189 25D0, | |
481 | 4 0.18343 46424 95649 80493 94761 42360 18D0, | |
482 | 5 0.98940 09349 91649 93259 61541 73450 33D0, | |
483 | 6 0.94457 50230 73232 57607 79884 15534 61D0, | |
484 | 7 0.86563 12023 87831 74388 04678 97712 39D0, | |
485 | 8 0.75540 44083 55003 03389 51011 94847 44D0, | |
486 | 9 0.61787 62444 02643 74844 66717 64048 79D0, | |
487 | A 0.45801 67776 57227 38634 24194 42983 58D0, | |
488 | B 0.28160 35507 79258 91323 04605 01460 50D0, | |
489 | C 0.95012 50983 76374 40185 31933 54249 58D-1/ | |
490 | ||
491 | DATA W | |
492 | 1 /0.10122 85362 90376 25915 25313 54309 96D0, | |
493 | 2 0.22238 10344 53374 47054 43559 94426 24D0, | |
494 | 3 0.31370 66458 77887 28733 79622 01986 60D0, | |
495 | 4 0.36268 37833 78361 98296 51504 49277 20D0, | |
496 | 5 0.27152 45941 17540 94851 78057 24560 18D-1, | |
497 | 6 0.62253 52393 86478 92862 84383 69943 78D-1, | |
498 | 7 0.95158 51168 24927 84809 92510 76022 46D-1, | |
499 | 8 0.12462 89712 55533 87205 24762 82192 02D0, | |
500 | 9 0.14959 59888 16576 73208 15017 30547 48D0, | |
501 | A 0.16915 65193 95002 53818 93120 79030 36D0, | |
502 | B 0.18260 34150 44923 58886 67636 67969 22D0, | |
503 | C 0.18945 06104 55068 49628 53967 23208 28D0/ | |
504 | EXTERNAL F | |
505 | H=0.D0 | |
506 | IF(B .EQ. A) GO TO 99 | |
507 | CONST=CST/ABS(B-A) | |
508 | BB=A | |
509 | 1 AA=BB | |
510 | BB=B | |
511 | 2 C1=HF*(BB+AA) | |
512 | C2=HF*(BB-AA) | |
513 | S8=0.D0 | |
514 | DO 3 I = 1,4 | |
515 | U=C2*X(I) | |
516 | 3 S8=S8+W(I)*(F(C1+U)+F(C1-U)) | |
517 | S16=0.D0 | |
518 | DO 4 I = 5,12 | |
519 | U=C2*X(I) | |
520 | 4 S16=S16+W(I)*(F(C1+U)+F(C1-U)) | |
521 | S16=C2*S16 | |
522 | IF(ABS(S16-C2*S8) .LE. EPS*(1.D0+ABS(S16))) THEN | |
523 | H=H+S16 | |
524 | IF(BB .NE. B) GO TO 1 | |
525 | ELSE | |
526 | BB=C1 | |
527 | IF(1.D0+CONST*ABS(C2) .NE. 1.D0) GO TO 2 | |
528 | H=0.D0 | |
529 | WRITE(6,*) 'DGAUSS1: TOO HIGH ACCURACY REQUIRED' | |
530 | GO TO 99 | |
531 | END IF | |
532 | 99 DGAUSS1=H | |
533 | RETURN | |
534 | END | |
535 | c | |
536 | FUNCTION SIMDIS(Numb,zmin,nzur,RI) | |
537 | C IT SIMULATES A RANDOM NUMBER ACCORDING TO A DISCRETE DISTRIBUTION GIVEN | |
538 | C BY ARRAY YA AT POINTS XA. THOUGHT FOR PYTHIA (PYR(0)). | |
539 | C N: NUMBER OF POINTS IN THE ARRAYS. | |
540 | C XA: ARRAY OF X-VALUES. | |
541 | C YA: ARRAY OF Y-VALUES. | |
542 | c RI: VALUE OF THE INTEGRAL. | |
543 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
544 | DIMENSION XA(500), YA(500) | |
545 | common/qpc1/eee,qhatl,omegac | |
546 | dlz=(1.d0-2.d0*zmin)/500.d0 | |
547 | do 1000 no=1,500 | |
548 | xa(no)=zmin+no*dlz | |
549 | if(nzur.eq.1) ya(no)=splitq2(xa(no)) | |
550 | if(nzur.eq.21) ya(no)=splitg2(xa(no)) | |
551 | if(nzur.eq.3) ya(no)=splitqqbar(xa(no)) | |
552 | 1000 continue | |
553 | RAL=PYR(0)*RI | |
554 | XAUX=0.D0 | |
555 | xauxold=0.d0 | |
556 | DO 1001 I=2, Numb, 1 | |
557 | XAUX=XAUX+(XA(I)-XA(I-1))*0.5D0* | |
558 | + (YA(I)+YA(I-1)) | |
559 | ||
560 | IF (XAUX .GE. RAL) GOTO 2011 | |
561 | ||
562 | IF (I .EQ. Numb) THEN | |
563 | SIMDIS=XA(I) | |
564 | ||
565 | RETURN | |
566 | ENDIF | |
567 | XAUXOLD=XAUX | |
568 | 1001 CONTINUE | |
569 | 2011 SIMDIS=(XA(I)-XA(I-1))*(RAL-XAUXOLD)/(XAUX-XAUXOLD)+ | |
570 | + XA(I-1) | |
571 | ||
572 | RETURN | |
573 | END | |
574 | ||
575 | ||
576 | SUBROUTINE QPYROBO(XI,YI,ZI,TI,THE,PHI,BEX,BEY,BEZ,XP,YP,ZP,TP) | |
577 | C N. Armesto, 16.04.2009 | |
578 | C performs a boost and rotation of (t,x,y,z) to (tp,xp,yp,zp): | |
579 | C cut version of PYROBO, angles and boost parameters identical. | |
580 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
581 | C...Local arrays. | |
582 | DIMENSION ROT(3,3),VR(3),DV(4) | |
583 | C | |
584 | X=XI | |
585 | Y=YI | |
586 | Z=ZI | |
587 | T=TI | |
588 | C...Rotate, typically from z axis to direction (theta,phi). | |
589 | IF(THE**2+PHI**2.GT.1D-20) THEN | |
590 | ROT(1,1)=COS(THE)*COS(PHI) | |
591 | ROT(1,2)=-SIN(PHI) | |
592 | ROT(1,3)=SIN(THE)*COS(PHI) | |
593 | ROT(2,1)=COS(THE)*SIN(PHI) | |
594 | ROT(2,2)=COS(PHI) | |
595 | ROT(2,3)=SIN(THE)*SIN(PHI) | |
596 | ROT(3,1)=-SIN(THE) | |
597 | ROT(3,2)=0D0 | |
598 | ROT(3,3)=COS(THE) | |
599 | C Instead of loop 120 in PYROBO. | |
600 | VR(1)=X | |
601 | VR(2)=Y | |
602 | VR(3)=Z | |
603 | C Instead of loop 130 in PYROBO. | |
604 | J=1 | |
605 | X=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3) | |
606 | J=2 | |
607 | Y=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3) | |
608 | J=3 | |
609 | Z=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3) | |
610 | ENDIF | |
611 | C If nothing happens... | |
612 | XP=X | |
613 | YP=Y | |
614 | ZP=Z | |
615 | TP=T | |
616 | C...Boost, typically from rest to momentum/energy=beta. | |
617 | IF(BEX**2+BEY**2+BEZ**2.GT.1D-20) THEN | |
618 | DBX=BEX | |
619 | DBY=BEY | |
620 | DBZ=BEZ | |
621 | DB=SQRT(DBX**2+DBY**2+DBZ**2) | |
622 | EPS1=1D0-1D-12 | |
623 | IF(DB.GT.EPS1) THEN | |
624 | C...Rescale boost vector if too close to unity. | |
625 | CALL PYERRM(3,'(PYROBO:) boost vector too large') | |
626 | DBX=DBX*(EPS1/DB) | |
627 | DBY=DBY*(EPS1/DB) | |
628 | DBZ=DBZ*(EPS1/DB) | |
629 | DB=EPS1 | |
630 | ENDIF | |
631 | DGA=1D0/SQRT(1D0-DB**2) | |
632 | C Instead of loop 150 in PYROBO. | |
633 | DV(1)=X | |
634 | DV(2)=Y | |
635 | DV(3)=Z | |
636 | DV(4)=T | |
637 | DBV=DBX*DV(1)+DBY*DV(2)+DBZ*DV(3) | |
638 | DGABV=DGA*(DGA*DBV/(1D0+DGA)+DV(4)) | |
639 | XP=DV(1)+DGABV*DBX | |
640 | YP=DV(2)+DGABV*DBY | |
641 | ZP=DV(3)+DGABV*DBZ | |
642 | TP=DGA*(DV(4)+DBV) | |
643 | ENDIF | |
644 | C | |
645 | RETURN | |
646 | END | |
647 | ||
648 | ||
649 | ||
650 | ||
651 | ||
652 | ||
653 | ||
654 | ||
655 | ||
656 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC | |
657 | C | |
658 | C PYSHOW ROUTINE FOR Q-PYTHIA version 1.0. | |
659 | C | |
660 | C DATE: 26.09.2008. | |
661 | C | |
662 | C AUTHORS: N. Armesto, L. Cunqueiro and C. A. Salgado | |
663 | C Departamento de Fisica de Particulas and IGFAE | |
664 | C Universidade de Santiago de Compostela | |
665 | C 15706 Santiago de Compostela, Spain | |
666 | C | |
667 | C EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es, | |
668 | C Carlos.Salgado@cern.ch | |
669 | C | |
670 | C CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18. | |
671 | C | |
672 | C WHEN USING Q-PYTHIA, PLEASE QUOTE: | |
673 | C | |
674 | C 1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado, | |
675 | C in preparation. | |
676 | C 2) T. Sjostrand, S. Mrenna and P. Skands, | |
677 | C ``PYTHIA 6.4 physics and manual,'' | |
678 | C JHEP 0605 (2006) 026 [arXiv:hep-ph/0603175]. | |
679 | C | |
680 | C INSTRUCTIONS: initial parton position is initialized by a call | |
681 | C to user-defined routine qpygin(x0,y0,z0,t0), | |
682 | C where these are the initial coordinates in the | |
683 | C center-of-mass frame of the hard collision | |
684 | C (if applicable for the type of process you study). | |
685 | C The values of qhatL and omegac have to be computed | |
686 | C by the user, using his preferred medium model, in | |
687 | C routine qpygeo, which takes as input the position | |
688 | C x,y,z,t of the parton to branch, the trajectory | |
689 | C defined by the three-vector betax,betay,betaz, | |
690 | C (all values in the center-of-mass frame of the | |
691 | C hard collision), and returns the value of qhatL | |
692 | C (in GeV**2) and omegac (in GeV). | |
693 | C Both routines are to be found at the end of this file. | |
694 | C | |
695 | C DISCLAIMER: this program comes without any guarantees. Beware of | |
696 | C errors and use common sense when interpreting results. | |
697 | C Any modifications are done under exclusive | |
698 | C maker's resposibility. | |
699 | C | |
700 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC | |
701 | C********************************************************************* | |
702 | ||
703 | C...PYSHOW | |
704 | C...Generates timelike parton showers from given partons. | |
705 | ||
706 | SUBROUTINE PYSHOWQ(IP1,IP2,QMAX) | |
707 | ||
708 | C...Double precision and integer declarations. | |
709 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
710 | IMPLICIT INTEGER(I-N) | |
711 | INTEGER PYK,PYCHGE,PYCOMP | |
712 | C...Parameter statement to help give large particle numbers. | |
713 | PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000, | |
714 | &KEXCIT=4000000,KDIMEN=5000000) | |
715 | PARAMETER (MAXNUR=500) | |
716 | Cacs+ | |
717 | PARAMETER (NNPOS=4000) | |
718 | DIMENSION PPOS(NNPOS,4) | |
719 | Cacs- | |
720 | C...Commonblocks. | |
721 | COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR) | |
722 | COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) | |
723 | COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) | |
724 | COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) | |
725 | COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) | |
726 | COMMON/PYINT1/MINT(400),VINT(400) | |
727 | SAVE /PYPART/,/PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/ | |
728 | ||
729 | Cacs+ | |
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/ | |
740 | Cacs- | |
741 | C...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) | |
747 | Cacs+ | |
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 | |
792 | Cacs- | |
793 | ||
794 | C...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 | ||
804 | C...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 | ||
832 | C...Send off to PYPTFS for pT-ordered evolution if requested, | |
833 | C...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 | ||
845 | C...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 | ||
892 | C...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) | |
902 | C...Special cutoff masses for initial partons (may be a heavy quark, | |
903 | C...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 | |
923 | C...QUARKONIA+++ | |
924 | C...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 | |
928 | C...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 | ||
960 | C...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 | ||
981 | C...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 | ||
1001 | C...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 | ||
1009 | C...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 | ||
1023 | C...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 | |
1031 | C...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 | |
1058 | C...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 | ||
1065 | C...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 | ||
1077 | C...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 | ||
1090 | C...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 | ||
1101 | C...~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 | ||
1111 | C...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 | ||
1119 | C...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 | ||
1145 | C...Boost interfering initial partons to rest frame | |
1146 | C...and reconstruct their polar and azimuthal angles. | |
1147 | Cacs+ | |
1148 | qplta1=0.d0 | |
1149 | qplta2=0.d0 | |
1150 | qpltbx=0.d0 | |
1151 | qpltby=0.d0 | |
1152 | qpltbz=0.d0 | |
1153 | Cacs- | |
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) | |
1182 | Cacs+ | |
1183 | qplta1=-the | |
1184 | qplta2=-phi | |
1185 | qpltbx=-PS(1)/PS(4) | |
1186 | qpltby=-PS(2)/PS(4) | |
1187 | qpltbz=-PS(3)/PS(4) | |
1188 | Cacs- | |
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 | ||
1200 | C...Boost 3 or more partons to their rest frame. | |
1201 | Cacs+ | |
1202 | c IF(NPA.GE.3) CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4), | |
1203 | c &-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 | |
1214 | Cacs- | |
1215 | ||
1216 | C...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 | ||
1242 | Cacs+ | |
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 | |
1249 | 10101 continue | |
1250 | ||
1251 | Cacs- | |
1252 | ||
1253 | ||
1254 | ||
1255 | ||
1256 | C...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 | ||
1276 | C...Position of aunt (sister to branching parton). | |
1277 | C...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 | ||
1312 | C...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 | ||
1326 | C...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 | ||
1351 | C...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 | ||
1377 | C...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 | ||
1387 | C...Store information on choice of evolving daughter. | |
1388 | IEP(1)=N+INUM | |
1389 | Cacs+ | |
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 | |
1429 | Cacs- | |
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 | ||
1448 | C...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 | ||
1461 | C...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 | ||
1475 | C...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 | ||
1505 | C...Integral of Altarelli-Parisi z kernel for QCD. | |
1506 | C...(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 | |
1508 | Cacs+ | |
1509 | C FBR= 6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0 | |
1510 | FBR=dgauss1(splitg1,zc,1.d0-zc,1.d-3) | |
1511 | Cacs- | |
1512 | C...QUARKONIA+++ | |
1513 | C...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) | |
1517 | C...QUARKONIA--- | |
1518 | ELSEIF(MSTJ(49).EQ.0) THEN | |
1519 | Cacs+ | |
1520 | C FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC) | |
1521 | FBR=dgauss1(splitq1,zc,1.d0-zc,1.d-3) | |
1522 | ||
1523 | Cacs- | |
1524 | IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0) | |
1525 | ||
1526 | C...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 | ||
1533 | C...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 | ||
1540 | C...Reset QCD probability for colourless. | |
1541 | IF(ISCOL(IR).EQ.0) FBR=0D0 | |
1542 | ||
1543 | C...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 | ||
1552 | C...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 | ||
1565 | C...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 | |
1570 | C...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) | |
1572 | C...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 | |
1584 | C...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 | ||
1590 | C...Select mass for daughter in QED evolution. | |
1591 | IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1.AND.IPSPD.EQ.0) THEN | |
1592 | C...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 | |
1600 | C...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 | ||
1610 | C...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 | |
1617 | Cacs+ | |
1618 | virt=V(IEP(1),5) | |
1619 | Cacs- | |
1620 | ||
1621 | C...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 | ||
1638 | C...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 | ||
1644 | C...QUARKONIA+++ | |
1645 | C...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) | |
1649 | C...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 | |
1653 | C...QUARKONIA--- | |
1654 | ||
1655 | C...Select z value of branching: q -> qg, g -> gg, g -> qqbar. | |
1656 | ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN | |
1657 | Cacs+ | |
1658 | C Z=1D0-(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0) | |
1659 | C...Only do z weighting when no ME correction afterwards. | |
1660 | C IF(M3JC.EQ.0.AND.1D0+Z**2.LT.2D0*PYR(0)) GOTO 410 | |
1661 | C | |
1662 | anfbr=dgauss1(splitq2,zc,1.d0-zc,1.d-3) | |
1663 | z=simdis(500,zc,1,anfbr) | |
1664 | Cacs- | |
1665 | K(IEP(1),5)=21 | |
1666 | ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN | |
1667 | Cacs+ | |
1668 | c 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 | |
1674 | c IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410 | |
1675 | Cacs- | |
1676 | K(IEP(1),5)=21 | |
1677 | ELSEIF(MSTJ(49).NE.1) THEN | |
1678 | Cacs+ | |
1679 | c Z=PYR(0) | |
1680 | c 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 | ||
1684 | Cacs- | |
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 | ||
1698 | C...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 | ||
1713 | C...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 | ||
1727 | C...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 | ||
1770 | C...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 | ||
1785 | C...Three-jet matrix element correction. | |
1786 | IF(M3JC.GE.1) THEN | |
1787 | WME=1D0 | |
1788 | WSHOW=1D0 | |
1789 | ||
1790 | C...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 | ||
1804 | C...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 | |
1810 | C...QCD ME: original parton, first branching. | |
1811 | PM2ME=PMTH(1,63-IR) | |
1812 | ECMME=PS(5) | |
1813 | ELSEIF(IR.GE.31) THEN | |
1814 | C...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 | |
1819 | C...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 | |
1830 | C...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 | |
1836 | C...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 | |
1841 | C...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 | |
1847 | C...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) | |
1856 | C...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 | ||
1862 | C...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 | ||
1876 | C...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 | ||
1909 | C...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 | ||
1923 | C...Impose angular constraint in first branching from interference | |
1924 | C...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 | ||
1934 | C...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 | ||
1947 | C...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 | ||
2036 | C...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 | |
2059 | C...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 | ||
2085 | C...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 | ||
2111 | C...Global statistics. | |
2112 | MINT(353)=MINT(353)+1 | |
2113 | VINT(353)=VINT(353)+PT | |
2114 | IF (MINT(353).EQ.1) VINT(358)=PT | |
2115 | ||
2116 | C...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 | ||
2136 | C...Find coefficient of azimuthal asymmetry due to soft gluon | |
2137 | C...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 | ||
2157 | C...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 | ||
2170 | C...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 | ||
2196 | C...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 | ||
2219 | C...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 | ||
2258 | C...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 | ||
2289 | C...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 | ||
2310 | C...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 | ||
2321 | C...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 | ||
2334 | C...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 | ||
2374 | C...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 | ||
2399 | C...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 | ||
2406 | C...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) | |
2428 | C NOT TO BE TOUCHED BY THE USER: | |
2429 | C IT SETS THE INITIAL POSITION AND TIME OF THE | |
2430 | C PARENT BRANCHING PARTON (X, Y, Z, T, IN FM) IN THE CENTER-OF-MASS | |
2431 | C FRAME OF THE HARD COLLISION (IF APPLICABLE FOR THE TYPE OF EVENTS | |
2432 | C YOU ARE SIMULATING). INFORMATION ABOUT THE BOOST AND ROTATION IS | |
2433 | C CONTAINED IN THE IN COMMON QPLT BELOW. | |
2434 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
2435 | C NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST | |
2436 | C PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE | |
2437 | C CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD | |
2438 | C SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO. | |
2439 | C XX, YY, ZZ, TT COMING FROM QPYGIN0 IN THE CMS OF THE COLLISION, GENERATED | |
2440 | C ONCE PER NUCLEON-NUCLEON COLLISION. | |
2441 | COMMON/TOQPYGIN/XX,YY,ZZ,TT | |
2442 | COMMON/QPLT/AA1,AA2,BBX,BBY,BBZ | |
2443 | C 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) | |
2447 | C | |
2448 | RETURN | |
2449 | END | |
2450 | ||
2451 | ||
2452 | ||
2453 | SUBROUTINE QPYGIN0 | |
2454 | C USER-DEFINED ROUTINE: IT SETS THE INITIAL POSITION AND TIME OF THE | |
2455 | C PARENT BRANCHING PARTON (X, Y, Z, T, IN FM) IN THE CENTER-OF-MASS | |
2456 | C FRAME OF THE COLLISION. | |
2457 | C IT SHOULD BE CALLED IN THE MAIN PROGRAM, ONCE PER NUCLEON-NUCLEON | |
2458 | C 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 | |
2461 | C Example: this is to be set by the user e.g. as coming from the overlap | |
2462 | C of nuclear profile functions (number of binary nucleon-nucleon | |
2463 | C 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 | |
2470 | C Passing to the common block. | |
2471 | xx=xin | |
2472 | yy=yin | |
2473 | zz=zin | |
2474 | tt=tin | |
2475 | c | |
2476 | c print*,"initial coordinates", xx,yy,zz,tt | |
2477 | RETURN | |
2478 | END | |
2479 | ||
2480 | ||
2481 | ||
2482 | ||
2483 | ||
2484 | ||
2485 | ||
2486 | ||
2487 | ||
2488 | ||
2489 | C | |
2490 | SUBROUTINE QPYGEO(x,y,z,t,bx,by,bz,qhl,oc) | |
2491 | C USER-DEFINED ROUTINE: | |
2492 | C The values of qhatL and omegac have to be computed | |
2493 | C by the user, using his preferred medium model, in | |
2494 | C this routine, which takes as input the position | |
2495 | C x,y,z,t (in fm) of the parton to branch, the trajectory | |
2496 | C defined by the three-vector bx,by,bz (in units of c), | |
2497 | C (all values in the center-of-mass frame of the | |
2498 | C hard collision), and returns the value of qhatL | |
2499 | C (in GeV**2) and omegac (in GeV). | |
2500 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
2501 | C NOW THE COMMON CONTAINING THE VALUES OF THE TWO ANGLES AND THREE BOSST | |
2502 | C PARAMETERS USED, IN PYSHOW, TO CHANGE THROUGH PYROBO FROM THE | |
2503 | C CENTER-OF-MASS OF THE COLLISION TO THE CENTER-OF-MASS OF THE HARD | |
2504 | C 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) | |
2510 | c Here we give five options for the geometry of the medium: | |
2511 | ||
2512 | c$$$cforalice+ | |
2513 | c$$$ (0) we call routine CalculateI0I1 in AliFastGlauber. Given the position | |
2514 | c$$$ of the parton in the reaction plane (x,y), the direction in the | |
2515 | c$$$ reaction plane phi=atan(py/px) and the impact parameter of the | |
2516 | c$$$ collision bimp, it gives back the transverse path length to the | |
2517 | c$$$ "end" of the medium and the integrated qhat along that path length. | |
2518 | c$$$ See the seter for this option in Configqpythia.C(one can set the | |
2519 | c$$$ value of xkscale by doing SetQhat(xkscale), with xkscale in fm). | |
2520 | c$$$ 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 | ||
2552 | c$$$cforalice- | |
2553 | c To use any of these folowing 1,2,3 or 4 options the user should specify | |
2554 | c a constant value for the transport coefficient and an initial in medium length. | |
2555 | c This can be done in the user Config file by setting: SetQhat(qhat), with qhat in | |
2556 | c GeV2/fm and SetLength(xl), with xl in fm. Those values are passed here through | |
2557 | c pythia free parameters parj(198) and parj(199). | |
2558 | ||
2559 | c (1) to fix the length to the initial value, uncomment the next three lines | |
2560 | c and comment the other definitions of xlp and qhl above and below. | |
2561 | c xlp=xl | |
2562 | c if (xlp .lt. 0.d0) xlp=0.d0 | |
2563 | c qhl=xlp*qhat ! GeV**2 | |
2564 | c print*, xlp,qhl | |
2565 | c (2) simplest ansatz: for an initial parton along the z-axis (approximate) | |
2566 | c starting in the center of a medium (-xl,+xl) along the z-axis | |
2567 | c if (bz .gt. 0.d0) then | |
2568 | c xlp=xl-z | |
2569 | c else | |
2570 | c xlp=xl+z | |
2571 | c endif | |
2572 | c if (xlp .gt. (2.d0*xl)) xlp=2.d0*xl | |
2573 | c if (xlp .lt. 0.d0) xlp=0.d0 | |
2574 | c qhl=xlp*qhat ! GeV**2 | |
2575 | ||
2576 | c (3) for a parton at midrapidity inside a cylinder (approximate) | |
2577 | c xlp=xl-dsqrt(x*x+y*y) | |
2578 | c if (xlp .lt. 0.d0) xlp=0.d0 | |
2579 | c qhl=xlp*qhat ! GeV**2 | |
2580 | ||
2581 | c (4) for a brick defined by planes (x,y,0) and (x,y,xl), comment | |
2582 | c the previous lines and uncomment lines between the comment 'brick'. | |
2583 | c brick+ | |
2584 | c if (z .ge. 0.d0 .and. z .le. xl) | |
2585 | c > then | |
2586 | c if (bz .gt. 0.d0) then | |
2587 | c ttpp=(xl-z)/bz | |
2588 | c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+ | |
2589 | c > (xl-z)**2.d0) | |
2590 | c else | |
2591 | c ttpp=z/dabs(bz) | |
2592 | c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+ | |
2593 | c > (z)**2.d0) | |
2594 | c endif | |
2595 | c elseif (z .lt. 0.d0) then | |
2596 | c if (bz .lt. 0.d0) then | |
2597 | c xlp=0.d0 | |
2598 | c else | |
2599 | c ttpp1=-z/bz | |
2600 | c ttpp2=(xl-z)/bz | |
2601 | c xxpp1=x+bx*ttpp1 | |
2602 | c xxpp2=x+bx*ttpp2 | |
2603 | c yypp1=y+by*ttpp1 | |
2604 | c yypp2=y+by*ttpp2 | |
2605 | c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+ | |
2606 | c > xl**2.d0) | |
2607 | c endif | |
2608 | c elseif (z .gt. xl) then | |
2609 | c if (bz .gt. 0.d0) then | |
2610 | c xlp=0.d0 | |
2611 | c else | |
2612 | c ttpp1=z/dabs(bz) | |
2613 | c ttpp2=(-xl+z)/dabs(bz) | |
2614 | c xxpp1=x+bx*ttpp1 | |
2615 | c xxpp2=x+bx*ttpp2 | |
2616 | c yypp1=y+by*ttpp1 | |
2617 | c yypp2=y+by*ttpp2 | |
2618 | c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+ | |
2619 | c > xl**2.d0) | |
2620 | c endif | |
2621 | c endif | |
2622 | c if (xlp .lt. 0.d0) xlp=0.d0 | |
2623 | c qhl=xlp*qhat ! GeV**2 | |
2624 | c brick- | |
2625 | ||
2626 | oc=0.5d0*qhl*xlp/0.1973d0 ! GeV | |
2627 | RETURN | |
2628 | END |