1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3 C AUXILIARY ROUTINES FOR Q-PYTHIA version 1.0.
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
12 C EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es,
13 C Carlos.Salgado@cern.ch
15 C CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18.
16 C NOT to be modified by user.
18 C WHEN USING Q-PYTHIA, PLEASE QUOTE:
20 C 1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado,
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].
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.
31 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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
38 c to integrate, adding vaccum plus medium q -> qg
39 implicit double precision (a-h,o-z)
42 auxq=splitgq(z)+splitmedq1(z)
43 if (auxq .gt. 0.d0) then
52 c to integrate, adding vaccum plus medium g -> gg, and g -> qqbar
53 implicit double precision (a-h,o-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))
67 c to integrate, adding vaccum plus medium q -> qg
68 implicit double precision (a-h,o-z)
70 auxq=splitgq(z)+splitmedq2(z)
71 if (auxq .gt. 0.d0) then
76 if(splitmedq2(z).eq.0.) then
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)
90 if(splitmedg2(z).eq.0.) then
96 c q -> qg splitting kernel at 1 loop for the vacuum
97 implicit double precision (a-h,o-z)
99 splitgq=(0.5d0*(xnc-1.d0/xnc))*(1.d0+z*z)/(1.d0-z)
104 c g -> gg splitting kernel at 1 loop for the vacuum
105 implicit double precision (a-h,o-z)
109 splitgg=xnc*auxz2*auxz2/auxz
113 function splitqqbar(z)
114 c g -> qqbar splitting kernel at 1 loop
115 implicit double precision (a-h,o-z)
118 splitqqbar=0.5d0*xnf*(z*z+auxz*auxz)
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
129 if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
132 c symmetrized by hand with respect to 1/2
133 if (z .ge. 0.5d0) then
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
149 c 1/2 to avoid double counting
150 splitmedg1=0.5*flx*xnc*2.d0*pi*zz*t*fff/qhatl
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
162 if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
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
175 splitmedq1=flx*(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
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
187 if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
190 c symmetrized by hand with respect to 1/2
191 if (z .ge. 0.5d0) then
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
207 c 1/2 to avoid double counting
208 splitmedg2=0.5*flx*xnc*2.d0*pi*zz*t*fff/qhatl
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
220 if (qhatl .le. 0.d0 .or. omegac .le. 0.d0) then
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
233 splitmedq2=flx*(0.5d0*(xnc-1.d0/xnc))*2.d0*pi*z*t*fff/qhatl
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.
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
253 c WE READ THE GRID ONLY THE FIRST TIME.
254 IF (IFLAG .EQ. 0) THEN
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')
266 read(11,*) xkap2(i), xlkap2(i)
269 read(11,*) xome(i), xlome(i)
273 read(11,*) xspec(i,j)
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
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)
299 genspec=scal*ddivdif(aux1,xlkap2,npkap,dlog(xk2),4)
303 if (ome .eq. xome(1)) then
307 do 60 i=1, npome-1, 1
308 if (ome .eq. xome(i+1)) then
311 elseif (ome .lt. xome(i+1)) then
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)
329 aux1(i)=xspec(i,iexact)
331 genspec=ddivdif(aux1,xlkap2,npkap,dlog(xk2),4)
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)
347 aux1(i)=xspec(i,iprev)
348 aux2(i)=xspec(i,ipost)
350 genprev=ddivdif(aux1,xlkap2,npkap,dlog(xk2),4)
351 genpost=ddivdif(aux2,xlkap2,npkap,dlog(xk2),4)
354 xl12=xlome(iprev)-xlome(ipost)
356 c2=genprev-c1*xlome(iprev)
357 genspec=c1*dlog(ome)+c2
365 * $Id: divdif.F,v 1.1.1.1 1996/02/15 17:48:36 mclareni Exp $
368 * Revision 1.1.1.1 1996/02/15 17:48:36 mclareni
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)
380 C TABULAR INTERPOLATION USING SYMMETRICALLY PLACED ARGUMENT POINTS.
382 C START. FIND SUBSCRIPT IX OF X IN ARRAY A.
383 IF( (NN.LT.2) .OR. (MM.LT.1) ) GO TO 601
389 IF(A(1).GT.A(N)) GO TO 4
390 C (SEARCH INCREASING ARGUMENTS.)
392 IF(X.GE.A(MID)) GO TO 2
397 3 IF(IY-IX.GT.1) GO TO 1
399 C (SEARCH DECREASING ARGUMENTS.)
401 IF(X.LE.A(MID)) GO TO 5
406 6 IF(IY-IX.GT.1) GO TO 4
408 C COPY REORDERED INTERPOLATION POINTS INTO (T(I),D(I)), SETTING
409 C *EXTRA* TO TRUE IF M+2 POINTS TO BE USED.
417 IF((1.LE.ISUB).AND.(ISUB.LE.N)) GO TO 501
425 11 IF(IP.LT.NPTS) GO TO 8
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.
431 IF(.NOT.EXTRA) GO TO 12
433 D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB))
437 D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB))
442 C EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES
443 C OF LAST DIFFERENCE IF *EXTRA* IS TRUE.
445 IF(EXTRA) SUM=0.5*(SUM+D(M+2))
448 SUM=D(J)+(X-T(J))*SUM
454 601 CALL KERMTR('E105.1',LGFILE,MFLAG,RFLAG)
458 IF(MM.LT.1) WRITE(*,101) MM
459 IF(NN.LT.2) WRITE(*,102) NN
461 IF(MM.LT.1) WRITE(LGFILE,101) MM
462 IF(NN.LT.2) WRITE(LGFILE,102) NN
465 IF(.NOT.RFLAG) CALL ABEND
467 101 FORMAT( 7X, 'FUNCTION DDIVDIF ... M =',I6,' IS LESS THAN 1')
468 102 FORMAT( 7X, 'FUNCTION DDIVDIF ... N =',I6,' IS LESS THAN 2')
471 C COPY OF CERN DGAUSS
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)
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/
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/
506 IF(B .EQ. A) GO TO 99
516 3 S8=S8+W(I)*(F(C1+U)+F(C1-U))
520 4 S16=S16+W(I)*(F(C1+U)+F(C1-U))
522 IF(ABS(S16-C2*S8) .LE. EPS*(1.D0+ABS(S16))) THEN
524 IF(BB .NE. B) GO TO 1
527 IF(1.D0+CONST*ABS(C2) .NE. 1.D0) GO TO 2
529 WRITE(6,*) 'DGAUSS1: TOO HIGH ACCURACY REQUIRED'
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
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))
557 XAUX=XAUX+(XA(I)-XA(I-1))*0.5D0*
560 IF (XAUX .GE. RAL) GOTO 2011
562 IF (I .EQ. Numb) THEN
569 2011 SIMDIS=(XA(I)-XA(I-1))*(RAL-XAUXOLD)/(XAUX-XAUXOLD)+
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)
582 DIMENSION ROT(3,3),VR(3),DV(4)
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)
592 ROT(1,3)=SIN(THE)*COS(PHI)
593 ROT(2,1)=COS(THE)*SIN(PHI)
595 ROT(2,3)=SIN(THE)*SIN(PHI)
599 C Instead of loop 120 in PYROBO.
603 C Instead of loop 130 in PYROBO.
605 X=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
607 Y=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
609 Z=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
611 C If nothing happens...
616 C...Boost, typically from rest to momentum/energy=beta.
617 IF(BEX**2+BEY**2+BEZ**2.GT.1D-20) THEN
621 DB=SQRT(DBX**2+DBY**2+DBZ**2)
624 C...Rescale boost vector if too close to unity.
625 CALL PYERRM(3,'(PYROBO:) boost vector too large')
631 DGA=1D0/SQRT(1D0-DB**2)
632 C Instead of loop 150 in PYROBO.
637 DBV=DBX*DV(1)+DBY*DV(2)+DBZ*DV(3)
638 DGABV=DGA*(DGA*DBV/(1D0+DGA)+DV(4))
656 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
658 C PYSHOW ROUTINE FOR Q-PYTHIA version 1.0.
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
667 C EMAILS: nestor@fpaxp1.usc.es, leticia@fpaxp1.usc.es,
668 C Carlos.Salgado@cern.ch
670 C CONTENT: auxiliary files for modified PYSHOW, fixed to PYTHIA-6.4.18.
672 C WHEN USING Q-PYTHIA, PLEASE QUOTE:
674 C 1) N. Armesto, G. Corcella, L. Cunqueiro and C. A. Salgado,
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].
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.
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.
700 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
701 C*********************************************************************
704 C...Generates timelike parton showers from given partons.
706 SUBROUTINE PYSHOWQ(IP1,IP2,QMAX)
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)
717 PARAMETER (NNPOS=4000)
718 DIMENSION PPOS(NNPOS,4)
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/
730 common/qpc1/eee,qhatl,omegac
733 COMMON/QPLT/QPLTA1,QPLTA2,QPLTBX,QPLTBY,QPLTBZ
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),
749 IF (IFLAG .EQ. 0) THEN
751 WRITE(MSTU(11),*) '*******************************************'
753 WRITE(MSTU(11),*) ' Q-PYTHIA version 1.0'
755 WRITE(MSTU(11),*) 'DATE: 26.09.2008'
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'
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'
768 WRITE(MSTU(11),*) 'NOTE: fixed to PYTHIA-6.4.18'
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].'
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.'
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.'
788 WRITE(MSTU(11),*) '*******************************************'
794 C...Check that QMAX not too low.
795 IF(MSTJ(41).LE.0) THEN
797 ELSEIF(MSTJ(41).EQ.1.OR.MSTJ(41).EQ.11) THEN
798 IF(QMAX.LE.PARJ(82).AND.IP2.GE.-80) RETURN
800 IF(QMAX.LE.MIN(PARJ(82),PARJ(83),PARJ(90)).AND.IP2.GE.-80)
804 C...Store positions of shower initiating partons.
806 IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
809 ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
814 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0
815 & .AND.IP2.GE.-80) THEN
820 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.
828 & '(PYSHOW:) failed to reconstruct showering system')
829 IF(MSTU(21).GE.1) RETURN
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.
839 PTPART(II)=0.5D0*QMAX
841 CALL PYPTFS(2,0.5D0*QMAX,0D0,PTGEN)
845 C...Initialization of cutoff masses etc.
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)
864 IF(MSTJ(41).GE.2) PMQTH1=MIN(PARJ(82),PARJ(83))
865 PMQT1E=MIN(PMQTH1,PARJ(90))
867 IF(MSTJ(41).GE.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
868 PMQT2E=MIN(PMQTH2,0.5D0*PARJ(90))
871 IF(MSTJ(41).GE.2) ISCHG(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)
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)
888 PT2MIN=MAX(0.5D0*PARJ(82),1.1D0*PARJ(81))**2
890 ALFM=LOG(PT2MIN/ALAMS)
892 C...Check on phase space available for emission.
900 KFLA(I)=IABS(K(IPA(I),2))
902 C...Special cutoff masses for initial partons (may be a heavy quark,
903 C...squark, ..., and need not be on the mass shell).
905 IF(NPA.LE.1) IREF(I)=IR
906 IF(NPA.GE.2) IREF(I+1)=IR
910 IF(KFLA(I).LE.8) THEN
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
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
921 ELSEIF(KFLA(I).EQ.KSUSY1+21) THEN
924 C...same for QQ~[3S18]
925 ELSEIF(MSTP(148).GE.1.AND.(KFLA(I).EQ.9900443.OR.
926 & KFLA(I).EQ.9900553)) THEN
930 IF(ISCOL(IR).EQ.1.OR.ISCHG(IR).EQ.1) KSH(IR)=1
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)
948 IF(KSH(IR).EQ.1) PMA(I)=PMTH(3,IR)
950 IF(KSH(IR).EQ.0.OR.PMA(I).GT.10D0*QMAX) IREJ=IREJ+1
952 PS(J)=PS(J)+P(IPA(I),J)
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
960 C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
963 ELSEIF(K(IP1,3).EQ.K(IP2,3).AND.K(IP1,3).GT.0) THEN
964 KFSRCE=IABS(K(K(IP1,3),2))
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))
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
981 C...Identify two primary showerers.
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
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
1001 C...Order of showerers. Presence of gluino.
1002 ITYPMN=MIN(ITYPE1,ITYPE2)
1003 ITYPMX=MAX(ITYPE1,ITYPE2)
1005 IF(ITYPE1.GT.ITYPE2) IORD=2
1007 IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1
1009 C...Check if 3-jet matrix elements to be used.
1012 IF(NPA.EQ.2.AND.MSTJ(47).GE.1.AND.MPSPD.EQ.0) THEN
1013 IF(MSTJ(38).NE.0) THEN
1017 ELSEIF(MSTJ(47).GE.6) THEN
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
1027 IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN
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.
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
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)))
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
1054 ALPHA=VECT/(VECT+AXIV)
1055 ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN
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
1061 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
1062 & ITYPES.EQ.1)) THEN
1065 C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
1066 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN
1068 IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN
1070 ELSEIF(KFSRCE.EQ.36) THEN
1073 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
1074 & ITYPES.EQ.1)) THEN
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
1081 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
1082 & ITYPES.EQ.2)) THEN
1084 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN
1086 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
1087 & ITYPES.EQ.2)) THEN
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
1094 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
1095 & ITYPES.EQ.2)) THEN
1097 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
1098 & ITYPES.EQ.1)) THEN
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
1104 ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
1105 & ITYPES.EQ.2)) THEN
1107 ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
1108 & ITYPES.EQ.1)) THEN
1111 C...g -> ~g + ~g (eikonal approximation).
1112 ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN
1115 M3JC=5*ICLASS+ICOMBI
1119 C...Find if interference with initial state partons.
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)
1129 IF(KCA.NE.0) KCII(I)=KCHG(KCA,2)*ISIGN(1,K(IPA(I),2))
1131 IF(KCII(I).NE.0) THEN
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
1137 IIIS(I,NIIS(I))=ICSI
1142 IF(NIIS(1)+NIIS(2).EQ.0) MIIS=0
1145 C...Boost interfering initial partons to rest frame
1146 C...and reconstruct their polar and azimuthal angles.
1157 K(N+I,J)=K(IPA(I),J)
1158 P(N+I,J)=P(IPA(I),J)
1162 DO 230 I=3,2+NIIS(1)
1164 K(N+I,J)=K(IIIS(1,I-2),J)
1165 P(N+I,J)=P(IIIS(1,I-2),J)
1169 DO 250 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
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)
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)
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))
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))
1200 C...Boost 3 or more partons to their rest frame.
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))
1205 CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
1206 &-PS(2)/PS(4),-PS(3)/PS(4))
1216 C...Define imagined single initiator of shower for parton system.
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
1243 call qpygin(pposx0,pposy0,pposz0,ppost0) ! in fm
1244 do 10101 iijj=1, nnpos, 1
1256 C...Loop over partons that may branch.
1259 IF(NPA.EQ.1) IM=NS-1
1262 IF(IM.GT.N) GOTO 600
1265 IF(KSH(IR).EQ.0) GOTO 290
1266 IF(P(IM,5).LT.PMTH(2,IR)) GOTO 290
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
1276 C...Position of aunt (sister to branching parton).
1277 C...Origin and flavour of daughters.
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
1293 K(N+I,2)=K(IPA(I),2)
1295 ELSEIF(KFLM.NE.21) THEN
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
1308 IREF(N+1-NS)=IABS(K(N+1,2))
1309 IREF(N+2-NS)=IABS(K(N+2,2))
1312 C...Reset flags on daughters and tries made.
1317 KFLD(IP)=IABS(K(N+IP,2))
1318 IF(KCHG(PYCOMP(KFLD(IP)),2).EQ.0) K(N+IP,1)=1
1322 IF(KSH(IREF(N+IP-NS)).EQ.1) ISI(IP)=1
1326 C...Maximum virtuality of daughters.
1329 IF(NPA.GE.3) P(N+I,4)=P(IPA(I),4)
1330 P(N+I,5)=MIN(QMAX,PS(5))
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)
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)
1344 IF(ISI(I).EQ.1) THEN
1346 IF(P(N+I,5).LE.PMTH(3,IR)) P(N+I,5)=PMTH(1,IR)
1348 V(N+I,5)=P(N+I,5)**2
1351 C...Choose one of the daughters for evolution.
1355 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
1358 IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
1360 IF(P(N+I,5).GE.PMTH(2,IR)) INUM=I
1366 IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQT2E) THEN
1367 RPM=P(N+I,5)/PMSD(I)
1369 IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,IR)) THEN
1377 C...Cancel choice of predetermined daughter already treated.
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
1387 C...Store information on choice of evolving daughter.
1394 if (nep .gt. 1 .and. inum .eq. 2) then
1398 if(zz1.gt.0.d0) then
1404 if (xkt .gt. 0.d0) then
1405 xlcoh=(2.d0*eee/(zz1*zz2*ttt))*0.1973d0
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)
1416 xbx=p(idf,1)/p(idf,4)
1417 xby=p(idf,2)/p(idf,4)
1418 xbz=p(idf,3)/p(idf,4)
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)
1433 IF(IEP(I).GT.N+NEP) IEP(I)=N+1
1436 KFL(I)=IABS(K(IEP(I),2))
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
1445 IF(KSH(IR).EQ.0) GOTO 450
1446 IF(P(IEP(1),5).LT.PMTH(2,IR)) GOTO 450
1448 C...Check if evolution already predetermined for daughter.
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
1456 IF(INUM.EQ.1.OR.INUM.EQ.2) THEN
1458 IF(IPSPD.NE.0) ISSET(INUM)=1
1461 C...Select side for interference with initial state partons.
1462 IF(MIIS.GE.1.AND.IEP(1).LE.NS+3) THEN
1465 IF(IABS(KCII(III)).EQ.1.AND.NIIS(III).EQ.1) THEN
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
1471 IF(PYR(0).GT.0.5D0) ISII(III)=2
1475 C...Calculate allowed z range.
1478 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1481 IF(INUM.EQ.1) PMED=V(IM,1)*PEM
1482 IF(INUM.EQ.2) PMED=(1D0-V(IM,1))*PEM
1484 IF(MOD(MSTJ(43),2).EQ.1) THEN
1487 IF(ISCOL(IR).EQ.0) ZCE=0.5D0*PARJ(90)/PMED
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
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
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
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
1509 C FBR= 6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
1510 FBR=dgauss1(splitg1,zc,1.d0-zc,1.d-3)
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)
1518 ELSEIF(MSTJ(49).EQ.0) THEN
1520 C FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC)
1521 FBR=dgauss1(splitq1,zc,1.d0-zc,1.d-3)
1524 IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0)
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
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)
1537 FBR=2D0*LOG((1D0-ZC)/ZC)
1540 C...Reset QCD probability for colourless.
1541 IF(ISCOL(IR).EQ.0) FBR=0D0
1543 C...Integral of Altarelli-Parisi kernel for photon emission.
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)
1549 IF(MSTJ(41).EQ.10) FBRE=PARJ(84)*FBRE
1552 C...Inner veto algorithm starts. Find maximum mass for evolution.
1559 IF(KSH(IRI).EQ.1) PM=PMTH(2,IRI)
1562 PMS=MIN(PMS,(P(IM,5)-PM2)**2)
1565 C...Select mass for daughter in QCD evolution.
1567 DO 430 IFF=4,MSTJ(45)
1568 IF(PMS.GT.4D0*PMTH(2,IFF)**2) B0=(33D0-2D0*IFF)/6D0
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.
1574 PMSQCD=P(IPSPD,5)**2
1575 ELSEIF(FBR.LT.1D-3) THEN
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))
1582 PMSQCD=PMSC*EXP(MAX(-50D0,ALFM*B0*LOG(PYR(0))/FBR))
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
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
1597 PMSQED=PMSE*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/
1598 & (PARU(101)*FBRE)))
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=
1604 IF(PMSQED.GT.PMSQCD) THEN
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
1621 C...Already predetermined choice of z, and flavour in g -> qqbar.
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
1635 K(IEP(1),5)=IABS(K(IPSGD1,2))
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
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
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
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
1662 anfbr=dgauss1(splitq2,zc,1.d0-zc,1.d-3)
1663 z=simdis(500,zc,1,anfbr)
1666 ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN
1668 c Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
1669 anfbr=dgauss1(splitg2,zc,1.d0-zc,1.d-3)
1671 z=simdis(500,zc,21,anfbr)
1673 IF(PYR(0).GT.0.5D0) Z=1D0-Z
1674 c IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
1677 ELSEIF(MSTJ(49).NE.1) THEN
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)
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
1694 IF((1D0+0.5D0*PMQ)*SQRT(1D0-PMQ).LT.PYR(0)) GOTO 410
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))
1702 ELSEIF(PYR(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
1703 Z=ZC+(1D0-2D0*ZC)*PYR(0)
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
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
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
1727 C...Check if z consistent with chosen m.
1728 IF(KFL(1).EQ.21) THEN
1729 IRGD1=IABS(K(IEP(1),5))
1733 IRGD2=IABS(K(IEP(1),5))
1737 ELSEIF(NEP.GE.3) THEN
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)
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
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-
1755 ZD=SQRT(MAX(0D0,1D0-V(IEP(1),5)/PED**2))
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
1764 IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 410
1766 IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1D0-ZL)/MAX(1D-20,ZL*
1768 IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
1770 C...Width suppression for q -> q + g.
1771 IF(MSTJ(40).NE.0.AND.KFL(1).NE.21.AND.IPSPD.EQ.0) THEN
1773 EGLU=0.5D0*PS(5)*(1D0-Z)*(1D0+V(IEP(1),5)/V(NS+1,5))
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
1785 C...Three-jet matrix element correction.
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)
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
1804 C...QCD matrix elements, including mass effects.
1805 ELSEIF(MSTJ(49).NE.1.AND.K(IEP(1),2).NE.21) THEN
1809 IF(IR.GE.31.AND.IGM.EQ.0) THEN
1810 C...QCD ME: original parton, first branching.
1813 ELSEIF(IR.GE.31) THEN
1814 C...QCD ME: original parton, subsequent branchings.
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.
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)/
1827 ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
1830 C...QCD ME: secondary partons, subsequent branchings.
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))
1836 C...Construct ME variables.
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)
1845 WME=PYMAEL(M3JCC,X2,X1,R2ME,R1ME,ALPHA)
1847 C...Split up total ME when two radiating partons.
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
1862 C...Toy model scalar theory matrix elements; no mass effects.
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)
1869 IF(MSTJ(102).GE.2) WME=X3**2-2D0*(1D0+X3)*(1D0-X1)*(1D0-X2)*
1873 IF(WME.LT.PYR(0)*WSHOW) GOTO 410
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
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
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
1888 PMDAO=PMTH(2,K(IEP(1),5))
1889 THE2ID=Z*(1D0-Z)*PEMAO**2/(V(IEP(1),5)-4D0*PMDAO**2)
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
1898 440 IF(K(IAOM,5).EQ.22) THEN
1900 IF(K(IAOM,3).LE.NS) MAOM=0
1901 IF(MAOM.EQ.1) GOTO 440
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
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
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
1934 C...End of inner veto algorithm. Check if only one leg evolved so far.
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
1942 IF(ITRY(I).EQ.0.AND.KSH(IR).EQ.1) THEN
1943 IF(P(N+I,5).GE.PMTH(2,IR)) GOTO 350
1947 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
1951 PMSUM=PMSUM+P(N+I,5)
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
1957 IF(KSH(IRDA).EQ.0) GOTO 480
1958 IF(P(I1,5).LT.PMTH(2,IRDA)) GOTO 480
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)
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)/
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-
1987 ZD=SQRT(MAX(0D0,1D0-V(I1,5)/PED**2))
1990 IF(IRDA.EQ.21.AND.IRGD1.LT.10.AND.
1991 & (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
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
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
2003 IF(IRDA.EQ.21) V(I1,4)=LOG(ZU*(1D0-ZL)/MAX(1D-20,
2005 IF(IRDA.NE.21) V(I1,4)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
2007 IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
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
2018 IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 350
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-
2032 IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 350
2036 C...Accepted branch. Construct four-momentum for initial partons.
2042 P(N+1,3)=SQRT(MAX(0D0,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
2044 P(N+1,4)=P(IPA(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)
2050 P(N+1,3)=SQRT(MAX(0D0,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
2055 P(N+2,4)=P(IM,5)-PED1
2058 ELSEIF(NEP.GE.3) THEN
2059 C...Rescale all momenta for energy conservation.
2065 P(N+I,J)=P(IPA(I),J)
2068 PQS=PQS+P(N+I,5)**2/P(N+I,4)
2071 FAC=(PS(5)-PQS)/(PES-PQS)
2076 P(N+I,J)=FAC*P(N+I,J)
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)
2081 PQS=PQS+P(N+I,5)**2/P(N+I,4)
2083 IF(LOOP.LT.10.AND.ABS(PES-PS(5)).GT.1D-12*PS(5)) GOTO 520
2085 C...Construct transverse momentum for ordinary branching in shower.
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)
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
2101 PTS=PMLS*(ZM*(1D0-ZM)*PEM**2/V(IM,5)-0.25D0)/PZM**2
2103 IF(PTS.LT.0D0.AND.LOOPPT.LT.10) THEN
2106 ELSEIF(PTS.LT.0D0) THEN
2109 PT=SQRT(MAX(0D0,PTS))
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
2116 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
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
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)
2127 HAZIP=(ZAU/(1D0-ZAU*(1D0-ZAU)))**2
2129 IF(K(N+1,2).NE.21) THEN
2130 HAZIP=HAZIP*(-2D0*ZM*(1D0-ZM))/(1D0-2D0*ZM*(1D0-ZM))
2132 HAZIP=HAZIP*(ZM*(1D0-ZM)/(1D0-ZM*(1D0-ZM)))**2
2136 C...Find coefficient of azimuthal asymmetry due to soft gluon
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
2147 IF(MAZIC.EQ.N+2) ZS=1D0-ZM
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)
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)
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)
2170 C...Already predetermined choice of phi angle or not
2173 IF(MPSPD.EQ.1.AND.IGM.EQ.NS+1) THEN
2175 IF(K(IPSPD,4).GT.0) THEN
2178 PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
2180 PHI=PYANGL(-P(IPSGD1,1),P(IPSGD1,2))
2183 ELSEIF(MPSPD.EQ.1.AND.IGM.EQ.NS+2) THEN
2185 IF(K(IPSPD,4).GT.0) THEN
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)
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
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)
2219 C...Rotate and boost daughters.
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)-
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)
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)
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)
2258 C...Weight with azimuthal distribution, if required.
2259 IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
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
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)
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))
2278 IF(1D0+HAZIP*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(HAZIP)))
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
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
2293 IF(ISII(III).GE.1) THEN
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
2310 C...Continue loop over partons that may branch, until none left.
2311 IF(IGM.GE.0) K(IM,1)=14
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
2321 C...Set information on imagined shower initiator.
2322 600 IF(NPA.GE.2) THEN
2326 IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
2334 C...Reconstruct string drawing information.
2336 KQ=KCHG(PYCOMP(K(I,2)),2)
2337 IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
2339 ELSEIF(K(I,1).LE.10.AND.IABS(K(I,2)).GE.11.AND.
2340 & IABS(K(I,2)).LE.18) THEN
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
2358 ID1=MOD(K(I,4),MSTU(5))
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
2374 C...Transformation from CM frame.
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))
2379 CALL PYROBO(NS+1,N,THE,PHI,0D0,0D0,0D0)
2380 ELSEIF(NPA.EQ.2) THEN
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)
2391 CALL PYROBO(NS+1,N,THE,PHI,BEX,BEY,BEZ)
2393 CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),
2396 CALL PYROBO(NS+1,N,0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),PS(3)/PS(4))
2399 C...Decay vertex of shower.
2406 C...Delete trivial shower, else connect initiators.
2407 IF(N.LE.NS+NPA+IIM) THEN
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)
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)
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
2464 call GetRandomXY(xrang,yrang)
2467 call savexy(xin, yin)
2470 C Passing to the common block.
2476 c print*,"initial coordinates", xx,yy,zz,tt
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)
2510 c Here we give five options for the geometry of the medium:
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).
2528 call qpyrobo(x,y,z,t,-aa1,-aa2,-bbx,-bby,-bbz,xout,
2531 call qpyrobo(bx,by,bz,1.d0,-aa1,-aa2,-bbx,-bby,-bbz,
2541 call CalculateI0I1(xlcero,xlone,bimpa,xa,ya,phia,ellcuta)
2542 if(xlcero.eq.0.d0) then
2546 xlp=2.d0*xlone/xlcero
2547 qhl=0.1973d0*0.1973d0*xlcero*xkscale
2548 call savei0i1(xlcero, xlone)
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).
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.
2562 c if (xlp .lt. 0.d0) xlp=0.d0
2563 c qhl=xlp*qhat ! GeV**2
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
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
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
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'.
2584 c if (z .ge. 0.d0 .and. z .le. xl)
2586 c if (bz .gt. 0.d0) then
2588 c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2592 c xlp=dsqrt((bx*ttpp)**2.d0+(by*ttpp)**2.d0+
2595 c elseif (z .lt. 0.d0) then
2596 c if (bz .lt. 0.d0) then
2605 c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2608 c elseif (z .gt. xl) then
2609 c if (bz .gt. 0.d0) then
2613 c ttpp2=(-xl+z)/dabs(bz)
2618 c xlp=dsqrt((xxpp1-xxpp2)**2.d0+(yypp1-yypp2)**2.d0+
2622 c if (xlp .lt. 0.d0) xlp=0.d0
2623 c qhl=xlp*qhat ! GeV**2
2626 oc=0.5d0*qhl*xlp/0.1973d0 ! GeV