Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lushow.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lushow.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lushow(IP1,IP2,QMAX)
5 
6 C...Purpose: to generate timelike parton showers from given partons.
7  IMPLICIT DOUBLE PRECISION(d)
8  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14  dimension pmth(5,40),ps(5),pma(4),pmsd(4),iep(4),ipa(4),
15  &kfla(4),kfld(4),kfl(4),itry(4),isi(4),isl(4),dp(4),dpt(5,4)
16 
17 C...Initialization of cutoff masses etc.
18  IF(mstj(41).LE.0.OR.(mstj(41).EQ.1.AND.qmax.LE.parj(82)).OR.
19  &qmax.LE.min(parj(82),parj(83)).OR.mstj(41).GE.3) RETURN
20  pmth(1,21)=ulmass(21)
21  pmth(2,21)=sqrt(pmth(1,21)**2+0.25*parj(82)**2)
22  pmth(3,21)=2.*pmth(2,21)
23  pmth(4,21)=pmth(3,21)
24  pmth(5,21)=pmth(3,21)
25  pmth(1,22)=ulmass(22)
26  pmth(2,22)=sqrt(pmth(1,22)**2+0.25*parj(83)**2)
27  pmth(3,22)=2.*pmth(2,22)
28  pmth(4,22)=pmth(3,22)
29  pmth(5,22)=pmth(3,22)
30  pmqth1=parj(82)
31  IF(mstj(41).EQ.2) pmqth1=min(parj(82),parj(83))
32  pmqth2=pmth(2,21)
33  IF(mstj(41).EQ.2) pmqth2=min(pmth(2,21),pmth(2,22))
34  DO 100 if=1,8
35  pmth(1,if)=ulmass(if)
36  pmth(2,if)=sqrt(pmth(1,if)**2+0.25*pmqth1**2)
37  pmth(3,if)=pmth(2,if)+pmqth2
38  pmth(4,if)=sqrt(pmth(1,if)**2+0.25*parj(82)**2)+pmth(2,21)
39  100 pmth(5,if)=sqrt(pmth(1,if)**2+0.25*parj(83)**2)+pmth(2,22)
40  pt2min=max(0.5*parj(82),1.1*parj(81))**2
41  alams=parj(81)**2
42  alfm=log(pt2min/alams)
43 
44 C...Store positions of shower initiating partons.
45  m3jc=0
46  IF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.ip2.EQ.0) THEN
47  npa=1
48  ipa(1)=ip1
49  ELSEIF(min(ip1,ip2).GT.0.AND.max(ip1,ip2).LE.min(n,mstu(4)-
50  &mstu(32))) THEN
51  npa=2
52  ipa(1)=ip1
53  ipa(2)=ip2
54  ELSEIF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.ip2.LT.0.
55  &and.ip2.GE.-3) THEN
56  npa=iabs(ip2)
57  DO 110 i=1,npa
58  110 ipa(i)=ip1+i-1
59  ELSE
60  CALL luerrm(12,
61  & '(LUSHOW:) failed to reconstruct showering system')
62  IF(mstu(21).GE.1) RETURN
63  ENDIF
64 
65 C...Check on phase space available for emission.
66  irej=0
67  DO 120 j=1,5
68  120 ps(j)=0.
69  pm=0.
70  DO 130 i=1,npa
71  kfla(i)=iabs(k(ipa(i),2))
72  pma(i)=p(ipa(i),5)
73  IF(kfla(i).NE.0.AND.(kfla(i).LE.8.OR.kfla(i).EQ.21))
74  &pma(i)=pmth(3,kfla(i))
75  pm=pm+pma(i)
76  IF(kfla(i).EQ.0.OR.(kfla(i).GT.8.AND.kfla(i).NE.21).OR.
77  &pma(i).GT.qmax) irej=irej+1
78  DO 130 j=1,4
79  130 ps(j)=ps(j)+p(ipa(i),j)
80  IF(irej.EQ.npa) RETURN
81  ps(5)=sqrt(max(0.,ps(4)**2-ps(1)**2-ps(2)**2-ps(3)**2))
82  IF(npa.EQ.1) ps(5)=ps(4)
83  IF(ps(5).LE.pm+pmqth1) RETURN
84  IF(npa.EQ.2.AND.mstj(47).GE.1) THEN
85  IF(kfla(1).GE.1.AND.kfla(1).LE.8.AND.kfla(2).GE.1.AND.
86  & kfla(2).LE.8) m3jc=1
87  IF(mstj(47).GE.2) m3jc=1
88  ENDIF
89 
90 C...Define imagined single initiator of shower for parton system.
91  ns=n
92  IF(n.GT.mstu(4)-mstu(32)-5) THEN
93  CALL luerrm(11,'(LUSHOW:) no more memory left in LUJETS')
94  IF(mstu(21).GE.1) RETURN
95  ENDIF
96  IF(npa.GE.2) THEN
97  k(n+1,1)=11
98  k(n+1,2)=21
99  k(n+1,3)=0
100  k(n+1,4)=0
101  k(n+1,5)=0
102  p(n+1,1)=0.
103  p(n+1,2)=0.
104  p(n+1,3)=0.
105  p(n+1,4)=ps(5)
106  p(n+1,5)=ps(5)
107  v(n+1,5)=ps(5)**2
108  n=n+1
109  ENDIF
110 
111 C...Loop over partons that may branch.
112  nep=npa
113  im=ns
114  IF(npa.EQ.1) im=ns-1
115  140 im=im+1
116  IF(n.GT.ns) THEN
117  IF(im.GT.n) goto 380
118  kflm=iabs(k(im,2))
119  IF(kflm.EQ.0.OR.(kflm.GT.8.AND.kflm.NE.21)) goto 140
120  IF(p(im,5).LT.pmth(2,kflm)) goto 140
121  igm=k(im,3)
122  ELSE
123  igm=-1
124  ENDIF
125  IF(n+nep.GT.mstu(4)-mstu(32)-5) THEN
126  CALL luerrm(11,'(LUSHOW:) no more memory left in LUJETS')
127  IF(mstu(21).GE.1) RETURN
128  ENDIF
129 
130 C...Position of aunt (sister to branching parton).
131 C...Origin and flavour of daughters.
132  iau=0
133  IF(igm.GT.0) THEN
134  IF(k(im-1,3).EQ.igm) iau=im-1
135  IF(n.GE.im+1.AND.k(im+1,3).EQ.igm) iau=im+1
136  ENDIF
137  IF(igm.GE.0) THEN
138  k(im,4)=n+1
139  DO 150 i=1,nep
140  150 k(n+i,3)=im
141  ELSE
142  k(n+1,3)=ipa(1)
143  ENDIF
144  IF(igm.LE.0) THEN
145  DO 160 i=1,nep
146  160 k(n+i,2)=k(ipa(i),2)
147  ELSEIF(kflm.NE.21) THEN
148  k(n+1,2)=k(im,2)
149  k(n+2,2)=k(im,5)
150  ELSEIF(k(im,5).EQ.21) THEN
151  k(n+1,2)=21
152  k(n+2,2)=21
153  ELSE
154  k(n+1,2)=k(im,5)
155  k(n+2,2)=-k(im,5)
156  ENDIF
157 
158 C...Reset flags on daughers and tries made.
159  DO 170 ip=1,nep
160  k(n+ip,1)=3
161  k(n+ip,4)=0
162  k(n+ip,5)=0
163  kfld(ip)=iabs(k(n+ip,2))
164  itry(ip)=0
165  isl(ip)=0
166  isi(ip)=0
167  170 IF(kfld(ip).GT.0.AND.(kfld(ip).LE.8.OR.kfld(ip).EQ.21)) isi(ip)=1
168  islm=0
169 
170 C...Maximum virtuality of daughters.
171  IF(igm.LE.0) THEN
172  DO 180 i=1,npa
173  IF(npa.GE.3) p(n+i,4)=(ps(4)*p(ipa(i),4)-ps(1)*p(ipa(i),1)-
174  & ps(2)*p(ipa(i),2)-ps(3)*p(ipa(i),3))/ps(5)
175  p(n+i,5)=min(qmax,ps(5))
176  IF(npa.GE.3) p(n+i,5)=min(p(n+i,5),p(n+i,4))
177  180 IF(isi(i).EQ.0) p(n+i,5)=p(ipa(i),5)
178  ELSE
179  IF(mstj(43).LE.2) pem=v(im,2)
180  IF(mstj(43).GE.3) pem=p(im,4)
181  p(n+1,5)=min(p(im,5),v(im,1)*pem)
182  p(n+2,5)=min(p(im,5),(1.-v(im,1))*pem)
183  IF(k(n+2,2).EQ.22) p(n+2,5)=pmth(1,22)
184  ENDIF
185  DO 190 i=1,nep
186  pmsd(i)=p(n+i,5)
187  IF(isi(i).EQ.1) THEN
188  IF(p(n+i,5).LE.pmth(3,kfld(i))) p(n+i,5)=pmth(1,kfld(i))
189  ENDIF
190  190 v(n+i,5)=p(n+i,5)**2
191 
192 C...Choose one of the daughters for evolution.
193  200 inum=0
194  IF(nep.EQ.1) inum=1
195  DO 210 i=1,nep
196  210 IF(inum.EQ.0.AND.isl(i).EQ.1) inum=i
197  DO 220 i=1,nep
198  IF(inum.EQ.0.AND.itry(i).EQ.0.AND.isi(i).EQ.1) THEN
199  IF(p(n+i,5).GE.pmth(2,kfld(i))) inum=i
200  ENDIF
201  220 CONTINUE
202  IF(inum.EQ.0) THEN
203  rmax=0.
204  DO 230 i=1,nep
205  IF(isi(i).EQ.1.AND.pmsd(i).GE.pmqth2) THEN
206  rpm=p(n+i,5)/pmsd(i)
207  IF(rpm.GT.rmax.AND.p(n+i,5).GE.pmth(2,kfld(i))) THEN
208  rmax=rpm
209  inum=i
210  ENDIF
211  ENDIF
212  230 CONTINUE
213  ENDIF
214 
215 C...Store information on choice of evolving daughter.
216  inum=max(1,inum)
217  iep(1)=n+inum
218  DO 240 i=2,nep
219  iep(i)=iep(i-1)+1
220  240 IF(iep(i).GT.n+nep) iep(i)=n+1
221  DO 250 i=1,nep
222  250 kfl(i)=iabs(k(iep(i),2))
223  itry(inum)=itry(inum)+1
224  IF(itry(inum).GT.200) THEN
225  CALL luerrm(14,'(LUSHOW:) caught in infinite loop')
226  IF(mstu(21).GE.1) RETURN
227  ENDIF
228  z=0.5
229  IF(kfl(1).EQ.0.OR.(kfl(1).GT.8.AND.kfl(1).NE.21)) goto 300
230  IF(p(iep(1),5).LT.pmth(2,kfl(1))) goto 300
231 
232 C...Calculate allowed z range.
233  IF(nep.EQ.1) THEN
234  pmed=ps(4)
235  ELSEIF(igm.EQ.0.OR.mstj(43).LE.2) THEN
236  pmed=p(im,5)
237  ELSE
238  IF(inum.EQ.1) pmed=v(im,1)*pem
239  IF(inum.EQ.2) pmed=(1.-v(im,1))*pem
240  ENDIF
241  IF(mod(mstj(43),2).EQ.1) THEN
242  zc=pmth(2,21)/pmed
243  zce=pmth(2,22)/pmed
244  ELSE
245  zc=0.5*(1.-sqrt(max(0.,1.-(2.*pmth(2,21)/pmed)**2)))
246  IF(zc.LT.1e-4) zc=(pmth(2,21)/pmed)**2
247  zce=0.5*(1.-sqrt(max(0.,1.-(2.*pmth(2,22)/pmed)**2)))
248  IF(zce.LT.1e-4) zce=(pmth(2,22)/pmed)**2
249  ENDIF
250  zc=min(zc,0.491)
251  zce=min(zce,0.491)
252  IF((mstj(41).EQ.1.AND.zc.GT.0.49).OR.(mstj(41).EQ.2.AND.
253  &min(zc,zce).GT.0.49)) THEN
254  p(iep(1),5)=pmth(1,kfl(1))
255  v(iep(1),5)=p(iep(1),5)**2
256  goto 300
257  ENDIF
258 
259 C...Integral of Altarelli-Parisi z kernel for QCD.
260  IF(mstj(49).EQ.0.AND.kfl(1).EQ.21) THEN
261  fbr=6.*log((1.-zc)/zc)+mstj(45)*(0.5-zc)
262  ELSEIF(mstj(49).EQ.0) THEN
263  fbr=(8./3.)*log((1.-zc)/zc)
264 
265 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
266  ELSEIF(mstj(49).EQ.1.AND.kfl(1).EQ.21) THEN
267  fbr=(parj(87)+mstj(45)*parj(88))*(1.-2.*zc)
268  ELSEIF(mstj(49).EQ.1) THEN
269  fbr=(1.-2.*zc)/3.
270  IF(igm.EQ.0.AND.m3jc.EQ.1) fbr=4.*fbr
271 
272 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
273  ELSEIF(kfl(1).EQ.21) THEN
274  fbr=6.*mstj(45)*(0.5-zc)
275  ELSE
276  fbr=2.*log((1.-zc)/zc)
277  ENDIF
278 
279 C...Integral of Altarelli-Parisi kernel for photon emission.
280  IF(mstj(41).EQ.2.AND.kfl(1).GE.1.AND.kfl(1).LE.8)
281  &fbre=(kchg(kfl(1),1)/3.)**2*2.*log((1.-zce)/zce)
282 
283 C...Inner veto algorithm starts. Find maximum mass for evolution.
284  260 pms=v(iep(1),5)
285  IF(igm.GE.0) THEN
286  pm2=0.
287  DO 270 i=2,nep
288  pm=p(iep(i),5)
289  IF(kfl(i).GT.0.AND.(kfl(i).LE.8.OR.kfl(i).EQ.21)) pm=
290  & pmth(2,kfl(i))
291  270 pm2=pm2+pm
292  pms=min(pms,(p(im,5)-pm2)**2)
293  ENDIF
294 
295 C...Select mass for daughter in QCD evolution.
296  b0=27./6.
297  DO 280 if=4,mstj(45)
298  280 IF(pms.GT.4.*pmth(2,if)**2) b0=(33.-2.*if)/6.
299  IF(mstj(44).LE.0) THEN
300  pmsqcd=pms*exp(max(-100.,log(rlu(0))*paru(2)/(paru(111)*fbr)))
301  ELSEIF(mstj(44).EQ.1) THEN
302  pmsqcd=4.*alams*(0.25*pms/alams)**(rlu(0)**(b0/fbr))
303  ELSE
304  pmsqcd=pms*rlu(0)**(alfm*b0/fbr)
305  ENDIF
306  IF(zc.GT.0.49.OR.pmsqcd.LE.pmth(4,kfl(1))**2) pmsqcd=
307  &pmth(2,kfl(1))**2
308  v(iep(1),5)=pmsqcd
309  mce=1
310 
311 C...Select mass for daughter in QED evolution.
312  IF(mstj(41).EQ.2.AND.kfl(1).GE.1.AND.kfl(1).LE.8) THEN
313  pmsqed=pms*exp(max(-100.,log(rlu(0))*paru(2)/(paru(101)*fbre)))
314  IF(zce.GT.0.49.OR.pmsqed.LE.pmth(5,kfl(1))**2) pmsqed=
315  & pmth(2,kfl(1))**2
316  IF(pmsqed.GT.pmsqcd) THEN
317  v(iep(1),5)=pmsqed
318  mce=2
319  ENDIF
320  ENDIF
321 
322 C...Check whether daughter mass below cutoff.
323  p(iep(1),5)=sqrt(v(iep(1),5))
324  IF(p(iep(1),5).LE.pmth(3,kfl(1))) THEN
325  p(iep(1),5)=pmth(1,kfl(1))
326  v(iep(1),5)=p(iep(1),5)**2
327  goto 300
328  ENDIF
329 
330 C...Select z value of branching: q -> qgamma.
331  IF(mce.EQ.2) THEN
332  z=1.-(1.-zce)*(zce/(1.-zce))**rlu(0)
333  IF(1.+z**2.LT.2.*rlu(0)) goto 260
334  k(iep(1),5)=22
335 
336 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
337  ELSEIF(mstj(49).NE.1.AND.kfl(1).NE.21) THEN
338  z=1.-(1.-zc)*(zc/(1.-zc))**rlu(0)
339  IF(1.+z**2.LT.2.*rlu(0)) goto 260
340  k(iep(1),5)=21
341  ELSEIF(mstj(49).EQ.0.AND.mstj(45)*(0.5-zc).LT.rlu(0)*fbr) THEN
342  z=(1.-zc)*(zc/(1.-zc))**rlu(0)
343  IF(rlu(0).GT.0.5) z=1.-z
344  IF((1.-z*(1.-z))**2.LT.rlu(0)) goto 260
345  k(iep(1),5)=21
346  ELSEIF(mstj(49).NE.1) THEN
347  z=zc+(1.-2.*zc)*rlu(0)
348  IF(z**2+(1.-z)**2.LT.rlu(0)) goto 260
349  kflb=1+int(mstj(45)*rlu(0))
350  pmq=4.*pmth(2,kflb)**2/v(iep(1),5)
351  IF(pmq.GE.1.) goto 260
352  pmq0=4.*pmth(2,21)**2/v(iep(1),5)
353  IF(mod(mstj(43),2).EQ.0.AND.(1.+0.5*pmq)*sqrt(1.-pmq).LT.
354  & rlu(0)*(1.+0.5*pmq0)*sqrt(1.-pmq0)) goto 260
355  k(iep(1),5)=kflb
356 
357 C...Ditto for scalar gluon model.
358  ELSEIF(kfl(1).NE.21) THEN
359  z=1.-sqrt(zc**2+rlu(0)*(1.-2.*zc))
360  k(iep(1),5)=21
361  ELSEIF(rlu(0)*(parj(87)+mstj(45)*parj(88)).LE.parj(87)) THEN
362  z=zc+(1.-2.*zc)*rlu(0)
363  k(iep(1),5)=21
364  ELSE
365  z=zc+(1.-2.*zc)*rlu(0)
366  kflb=1+int(mstj(45)*rlu(0))
367  pmq=4.*pmth(2,kflb)**2/v(iep(1),5)
368  IF(pmq.GE.1.) goto 260
369  k(iep(1),5)=kflb
370  ENDIF
371  IF(mce.EQ.1.AND.mstj(44).GE.2) THEN
372  IF(z*(1.-z)*v(iep(1),5).LT.pt2min) goto 260
373  IF(alfm/log(v(iep(1),5)*z*(1.-z)/alams).LT.rlu(0)) goto 260
374  ENDIF
375 
376 C...Check if z consistent with chosen m.
377  IF(kfl(1).EQ.21) THEN
378  kflgd1=iabs(k(iep(1),5))
379  kflgd2=kflgd1
380  ELSE
381  kflgd1=kfl(1)
382  kflgd2=iabs(k(iep(1),5))
383  ENDIF
384  IF(nep.EQ.1) THEN
385  ped=ps(4)
386  ELSEIF(nep.GE.3) THEN
387  ped=p(iep(1),4)
388  ELSEIF(igm.EQ.0.OR.mstj(43).LE.2) THEN
389  ped=0.5*(v(im,5)+v(iep(1),5)-pm2**2)/p(im,5)
390  ELSE
391  IF(iep(1).EQ.n+1) ped=v(im,1)*pem
392  IF(iep(1).EQ.n+2) ped=(1.-v(im,1))*pem
393  ENDIF
394  IF(mod(mstj(43),2).EQ.1) THEN
395  pmqth3=0.5*parj(82)
396  IF(kflgd2.EQ.22) pmqth3=0.5*parj(83)
397  pmq1=(pmth(1,kflgd1)**2+pmqth3**2)/v(iep(1),5)
398  pmq2=(pmth(1,kflgd2)**2+pmqth3**2)/v(iep(1),5)
399  zd=sqrt(max(0.,(1.-v(iep(1),5)/ped**2)*((1.-pmq1-pmq2)**2-
400  & 4.*pmq1*pmq2)))
401  zh=1.+pmq1-pmq2
402  ELSE
403  zd=sqrt(max(0.,1.-v(iep(1),5)/ped**2))
404  zh=1.
405  ENDIF
406  zl=0.5*(zh-zd)
407  zu=0.5*(zh+zd)
408  IF(z.LT.zl.OR.z.GT.zu) goto 260
409  IF(kfl(1).EQ.21) v(iep(1),3)=log(zu*(1.-zl)/max(1e-20,zl*
410  &(1.-zu)))
411  IF(kfl(1).NE.21) v(iep(1),3)=log((1.-zl)/max(1e-10,1.-zu))
412 
413 C...Three-jet matrix element correction.
414  IF(igm.EQ.0.AND.m3jc.EQ.1) THEN
415  x1=z*(1.+v(iep(1),5)/v(ns+1,5))
416  x2=1.-v(iep(1),5)/v(ns+1,5)
417  x3=(1.-x1)+(1.-x2)
418  IF(mce.EQ.2) THEN
419  ki1=k(ipa(inum),2)
420  ki2=k(ipa(3-inum),2)
421  qf1=kchg(iabs(ki1),1)*isign(1,ki1)/3.
422  qf2=kchg(iabs(ki2),1)*isign(1,ki2)/3.
423  wshow=qf1**2*(1.-x1)/x3*(1.+(x1/(2.-x2))**2)+
424  & qf2**2*(1.-x2)/x3*(1.+(x2/(2.-x1))**2)
425  wme=(qf1*(1.-x1)/x3-qf2*(1.-x2)/x3)**2*(x1**2+x2**2)
426  ELSEIF(mstj(49).NE.1) THEN
427  wshow=1.+(1.-x1)/x3*(x1/(2.-x2))**2+
428  & (1.-x2)/x3*(x2/(2.-x1))**2
429  wme=x1**2+x2**2
430  ELSE
431  wshow=4.*x3*((1.-x1)/(2.-x2)**2+(1.-x2)/(2.-x1)**2)
432  wme=x3**2
433  ENDIF
434  IF(wme.LT.rlu(0)*wshow) goto 260
435 
436 C...Impose angular ordering by rejection of nonordered emission.
437  ELSEIF(mce.EQ.1.AND.igm.GT.0.AND.mstj(42).GE.2) THEN
438  maom=1
439  zm=v(im,1)
440  IF(iep(1).EQ.n+2) zm=1.-v(im,1)
441  the2id=z*(1.-z)*(zm*p(im,4))**2/v(iep(1),5)
442  iaom=im
443  290 IF(k(iaom,5).EQ.22) THEN
444  iaom=k(iaom,3)
445  IF(k(iaom,3).LE.ns) maom=0
446  IF(maom.EQ.1) goto 290
447  ENDIF
448  IF(maom.EQ.1) THEN
449  the2im=v(iaom,1)*(1.-v(iaom,1))*p(iaom,4)**2/v(iaom,5)
450  IF(the2id.LT.the2im) goto 260
451  ENDIF
452  ENDIF
453 
454 C...Impose user-defined maximum angle at first branching.
455  IF(mstj(48).EQ.1) THEN
456  IF(nep.EQ.1.AND.im.EQ.ns) THEN
457  the2id=z*(1.-z)*ps(4)**2/v(iep(1),5)
458  IF(the2id.LT.1./parj(85)**2) goto 260
459  ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+2) THEN
460  the2id=z*(1.-z)*(0.5*p(im,4))**2/v(iep(1),5)
461  IF(the2id.LT.1./parj(85)**2) goto 260
462  ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+3) THEN
463  the2id=z*(1.-z)*(0.5*p(im,4))**2/v(iep(1),5)
464  IF(the2id.LT.1./parj(86)**2) goto 260
465  ENDIF
466  ENDIF
467 
468 C...End of inner veto algorithm. Check if only one leg evolved so far.
469  300 v(iep(1),1)=z
470  isl(1)=0
471  isl(2)=0
472  IF(nep.EQ.1) goto 330
473  IF(nep.EQ.2.AND.p(iep(1),5)+p(iep(2),5).GE.p(im,5)) goto 200
474  DO 310 i=1,nep
475  IF(itry(i).EQ.0.AND.kfld(i).GT.0.AND.(kfld(i).LE.8.OR.kfld(i).EQ.
476  &21)) THEN
477  IF(p(n+i,5).GE.pmth(2,kfld(i))) goto 200
478  ENDIF
479  310 CONTINUE
480 
481 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
482  IF(nep.EQ.3) THEN
483  pa1s=(p(n+1,4)+p(n+1,5))*(p(n+1,4)-p(n+1,5))
484  pa2s=(p(n+2,4)+p(n+2,5))*(p(n+2,4)-p(n+2,5))
485  pa3s=(p(n+3,4)+p(n+3,5))*(p(n+3,4)-p(n+3,5))
486  pts=0.25*(2.*pa1s*pa2s+2.*pa1s*pa3s+2.*pa2s*pa3s-
487  & pa1s**2-pa2s**2-pa3s**2)/pa1s
488  IF(pts.LE.0.) goto 200
489  ELSEIF(igm.EQ.0.OR.mstj(43).LE.2.OR.mod(mstj(43),2).EQ.0) THEN
490  DO 320 i1=n+1,n+2
491  kflda=iabs(k(i1,2))
492  IF(kflda.EQ.0.OR.(kflda.GT.8.AND.kflda.NE.21)) goto 320
493  IF(p(i1,5).LT.pmth(2,kflda)) goto 320
494  IF(kflda.EQ.21) THEN
495  kflgd1=iabs(k(i1,5))
496  kflgd2=kflgd1
497  ELSE
498  kflgd1=kflda
499  kflgd2=iabs(k(i1,5))
500  ENDIF
501  i2=2*n+3-i1
502  IF(igm.EQ.0.OR.mstj(43).LE.2) THEN
503  ped=0.5*(v(im,5)+v(i1,5)-v(i2,5))/p(im,5)
504  ELSE
505  IF(i1.EQ.n+1) zm=v(im,1)
506  IF(i1.EQ.n+2) zm=1.-v(im,1)
507  pml=sqrt((v(im,5)-v(n+1,5)-v(n+2,5))**2-
508  & 4.*v(n+1,5)*v(n+2,5))
509  ped=pem*(0.5*(v(im,5)-pml+v(i1,5)-v(i2,5))+pml*zm)/v(im,5)
510  ENDIF
511  IF(mod(mstj(43),2).EQ.1) THEN
512  pmqth3=0.5*parj(82)
513  IF(kflgd2.EQ.22) pmqth3=0.5*parj(83)
514  pmq1=(pmth(1,kflgd1)**2+pmqth3**2)/v(i1,5)
515  pmq2=(pmth(1,kflgd2)**2+pmqth3**2)/v(i1,5)
516  zd=sqrt(max(0.,(1.-v(i1,5)/ped**2)*((1.-pmq1-pmq2)**2-
517  & 4.*pmq1*pmq2)))
518  zh=1.+pmq1-pmq2
519  ELSE
520  zd=sqrt(max(0.,1.-v(i1,5)/ped**2))
521  zh=1.
522  ENDIF
523  zl=0.5*(zh-zd)
524  zu=0.5*(zh+zd)
525  IF(i1.EQ.n+1.AND.(v(i1,1).LT.zl.OR.v(i1,1).GT.zu)) isl(1)=1
526  IF(i1.EQ.n+2.AND.(v(i1,1).LT.zl.OR.v(i1,1).GT.zu)) isl(2)=1
527  IF(kflda.EQ.21) v(i1,4)=log(zu*(1.-zl)/max(1e-20,zl*(1.-zu)))
528  IF(kflda.NE.21) v(i1,4)=log((1.-zl)/max(1e-10,1.-zu))
529  320 CONTINUE
530  IF(isl(1).EQ.1.AND.isl(2).EQ.1.AND.islm.NE.0) THEN
531  isl(3-islm)=0
532  islm=3-islm
533  ELSEIF(isl(1).EQ.1.AND.isl(2).EQ.1) THEN
534  zdr1=max(0.,v(n+1,3)/v(n+1,4)-1.)
535  zdr2=max(0.,v(n+2,3)/v(n+2,4)-1.)
536  IF(zdr2.GT.rlu(0)*(zdr1+zdr2)) isl(1)=0
537  IF(isl(1).EQ.1) isl(2)=0
538  IF(isl(1).EQ.0) islm=1
539  IF(isl(2).EQ.0) islm=2
540  ENDIF
541  IF(isl(1).EQ.1.OR.isl(2).EQ.1) goto 200
542  ENDIF
543  IF(igm.GT.0.AND.mod(mstj(43),2).EQ.1.AND.(p(n+1,5).GE.
544  &pmth(2,kfld(1)).OR.p(n+2,5).GE.pmth(2,kfld(2)))) THEN
545  pmq1=v(n+1,5)/v(im,5)
546  pmq2=v(n+2,5)/v(im,5)
547  zd=sqrt(max(0.,(1.-v(im,5)/pem**2)*((1.-pmq1-pmq2)**2-
548  & 4.*pmq1*pmq2)))
549  zh=1.+pmq1-pmq2
550  zl=0.5*(zh-zd)
551  zu=0.5*(zh+zd)
552  IF(v(im,1).LT.zl.OR.v(im,1).GT.zu) goto 200
553  ENDIF
554 
555 C...Accepted branch. Construct four-momentum for initial partons.
556  330 mazip=0
557  mazic=0
558  IF(nep.EQ.1) THEN
559  p(n+1,1)=0.
560  p(n+1,2)=0.
561  p(n+1,3)=sqrt(max(0.,(p(ipa(1),4)+p(n+1,5))*(p(ipa(1),4)-
562  & p(n+1,5))))
563  p(n+1,4)=p(ipa(1),4)
564  v(n+1,2)=p(n+1,4)
565  ELSEIF(igm.EQ.0.AND.nep.EQ.2) THEN
566  ped1=0.5*(v(im,5)+v(n+1,5)-v(n+2,5))/p(im,5)
567  p(n+1,1)=0.
568  p(n+1,2)=0.
569  p(n+1,3)=sqrt(max(0.,(ped1+p(n+1,5))*(ped1-p(n+1,5))))
570  p(n+1,4)=ped1
571  p(n+2,1)=0.
572  p(n+2,2)=0.
573  p(n+2,3)=-p(n+1,3)
574  p(n+2,4)=p(im,5)-ped1
575  v(n+1,2)=p(n+1,4)
576  v(n+2,2)=p(n+2,4)
577  ELSEIF(nep.EQ.3) THEN
578  p(n+1,1)=0.
579  p(n+1,2)=0.
580  p(n+1,3)=sqrt(max(0.,pa1s))
581  p(n+2,1)=sqrt(pts)
582  p(n+2,2)=0.
583  p(n+2,3)=0.5*(pa3s-pa2s-pa1s)/p(n+1,3)
584  p(n+3,1)=-p(n+2,1)
585  p(n+3,2)=0.
586  p(n+3,3)=-(p(n+1,3)+p(n+2,3))
587  v(n+1,2)=p(n+1,4)
588  v(n+2,2)=p(n+2,4)
589  v(n+3,2)=p(n+3,4)
590 
591 C...Construct transverse momentum for ordinary branching in shower.
592  ELSE
593  zm=v(im,1)
594  pzm=sqrt(max(0.,(pem+p(im,5))*(pem-p(im,5))))
595  pmls=(v(im,5)-v(n+1,5)-v(n+2,5))**2-4.*v(n+1,5)*v(n+2,5)
596  IF(pzm.LE.0.) THEN
597  pts=0.
598  ELSEIF(mod(mstj(43),2).EQ.1) THEN
599  pts=(pem**2*(zm*(1.-zm)*v(im,5)-(1.-zm)*v(n+1,5)-
600  & zm*v(n+2,5))-0.25*pmls)/pzm**2
601  ELSE
602  pts=pmls*(zm*(1.-zm)*pem**2/v(im,5)-0.25)/pzm**2
603  ENDIF
604  pt=sqrt(max(0.,pts))
605 
606 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
607  hazip=0.
608  IF(mstj(49).NE.1.AND.mod(mstj(46),2).EQ.1.AND.k(im,2).EQ.21.
609  & and.iau.NE.0) THEN
610  IF(k(igm,3).NE.0) mazip=1
611  zau=v(igm,1)
612  IF(iau.EQ.im+1) zau=1.-v(igm,1)
613  IF(mazip.EQ.0) zau=0.
614  IF(k(igm,2).NE.21) THEN
615  hazip=2.*zau/(1.+zau**2)
616  ELSE
617  hazip=(zau/(1.-zau*(1.-zau)))**2
618  ENDIF
619  IF(k(n+1,2).NE.21) THEN
620  hazip=hazip*(-2.*zm*(1.-zm))/(1.-2.*zm*(1.-zm))
621  ELSE
622  hazip=hazip*(zm*(1.-zm)/(1.-zm*(1.-zm)))**2
623  ENDIF
624  ENDIF
625 
626 C...Find coefficient of azimuthal asymmetry due to soft gluon
627 C...interference.
628  hazic=0.
629  IF(mstj(46).GE.2.AND.(k(n+1,2).EQ.21.OR.k(n+2,2).EQ.21).
630  & and.iau.NE.0) THEN
631  IF(k(igm,3).NE.0) mazic=n+1
632  IF(k(igm,3).NE.0.AND.k(n+1,2).NE.21) mazic=n+2
633  IF(k(igm,3).NE.0.AND.k(n+1,2).EQ.21.AND.k(n+2,2).EQ.21.AND.
634  & zm.GT.0.5) mazic=n+2
635  IF(k(iau,2).EQ.22) mazic=0
636  zs=zm
637  IF(mazic.EQ.n+2) zs=1.-zm
638  zgm=v(igm,1)
639  IF(iau.EQ.im-1) zgm=1.-v(igm,1)
640  IF(mazic.EQ.0) zgm=1.
641  hazic=(p(im,5)/p(igm,5))*sqrt((1.-zs)*(1.-zgm)/(zs*zgm))
642  hazic=min(0.95,hazic)
643  ENDIF
644  ENDIF
645 
646 C...Construct kinematics for ordinary branching in shower.
647  340 IF(nep.EQ.2.AND.igm.GT.0) THEN
648  IF(mod(mstj(43),2).EQ.1) THEN
649  p(n+1,4)=pem*v(im,1)
650  ELSE
651  p(n+1,4)=pem*(0.5*(v(im,5)-sqrt(pmls)+v(n+1,5)-v(n+2,5))+
652  & sqrt(pmls)*zm)/v(im,5)
653  ENDIF
654  phi=paru(2)*rlu(0)
655  p(n+1,1)=pt*cos(phi)
656  p(n+1,2)=pt*sin(phi)
657  IF(pzm.GT.0.) THEN
658  p(n+1,3)=0.5*(v(n+2,5)-v(n+1,5)-v(im,5)+2.*pem*p(n+1,4))/pzm
659  ELSE
660  p(n+1,3)=0.
661  ENDIF
662  p(n+2,1)=-p(n+1,1)
663  p(n+2,2)=-p(n+1,2)
664  p(n+2,3)=pzm-p(n+1,3)
665  p(n+2,4)=pem-p(n+1,4)
666  IF(mstj(43).LE.2) THEN
667  v(n+1,2)=(pem*p(n+1,4)-pzm*p(n+1,3))/p(im,5)
668  v(n+2,2)=(pem*p(n+2,4)-pzm*p(n+2,3))/p(im,5)
669  ENDIF
670  ENDIF
671 
672 C...Rotate and boost daughters.
673  IF(igm.GT.0) THEN
674  IF(mstj(43).LE.2) THEN
675  bex=p(igm,1)/p(igm,4)
676  bey=p(igm,2)/p(igm,4)
677  bez=p(igm,3)/p(igm,4)
678  ga=p(igm,4)/p(igm,5)
679  gabep=ga*(ga*(bex*p(im,1)+bey*p(im,2)+bez*p(im,3))/(1.+ga)-
680  & p(im,4))
681  ELSE
682  bex=0.
683  bey=0.
684  bez=0.
685  ga=1.
686  gabep=0.
687  ENDIF
688  the=ulangl(p(im,3)+gabep*bez,sqrt((p(im,1)+gabep*bex)**2+
689  & (p(im,2)+gabep*bey)**2))
690  phi=ulangl(p(im,1)+gabep*bex,p(im,2)+gabep*bey)
691  DO 350 i=n+1,n+2
692  dp(1)=cos(the)*cos(phi)*p(i,1)-sin(phi)*p(i,2)+
693  & sin(the)*cos(phi)*p(i,3)
694  dp(2)=cos(the)*sin(phi)*p(i,1)+cos(phi)*p(i,2)+
695  & sin(the)*sin(phi)*p(i,3)
696  dp(3)=-sin(the)*p(i,1)+cos(the)*p(i,3)
697  dp(4)=p(i,4)
698  dbp=bex*dp(1)+bey*dp(2)+bez*dp(3)
699  dgabp=ga*(ga*dbp/(1d0+ga)+dp(4))
700  p(i,1)=dp(1)+dgabp*bex
701  p(i,2)=dp(2)+dgabp*bey
702  p(i,3)=dp(3)+dgabp*bez
703  350 p(i,4)=ga*(dp(4)+dbp)
704  ENDIF
705 
706 C...Weight with azimuthal distribution, if required.
707  IF(mazip.NE.0.OR.mazic.NE.0) THEN
708  DO 360 j=1,3
709  dpt(1,j)=p(im,j)
710  dpt(2,j)=p(iau,j)
711  360 dpt(3,j)=p(n+1,j)
712  dpma=dpt(1,1)*dpt(2,1)+dpt(1,2)*dpt(2,2)+dpt(1,3)*dpt(2,3)
713  dpmd=dpt(1,1)*dpt(3,1)+dpt(1,2)*dpt(3,2)+dpt(1,3)*dpt(3,3)
714  dpmm=dpt(1,1)**2+dpt(1,2)**2+dpt(1,3)**2
715  DO 370 j=1,3
716  dpt(4,j)=dpt(2,j)-dpma*dpt(1,j)/dpmm
717  370 dpt(5,j)=dpt(3,j)-dpmd*dpt(1,j)/dpmm
718  dpt(4,4)=sqrt(dpt(4,1)**2+dpt(4,2)**2+dpt(4,3)**2)
719  dpt(5,4)=sqrt(dpt(5,1)**2+dpt(5,2)**2+dpt(5,3)**2)
720  IF(min(dpt(4,4),dpt(5,4)).GT.0.1*parj(82)) THEN
721  cad=(dpt(4,1)*dpt(5,1)+dpt(4,2)*dpt(5,2)+
722  & dpt(4,3)*dpt(5,3))/(dpt(4,4)*dpt(5,4))
723  IF(mazip.NE.0) THEN
724  IF(1.+hazip*(2.*cad**2-1.).LT.rlu(0)*(1.+abs(hazip)))
725  & goto 340
726  ENDIF
727  IF(mazic.NE.0) THEN
728  IF(mazic.EQ.n+2) cad=-cad
729  IF((1.-hazic)*(1.-hazic*cad)/(1.+hazic**2-2.*hazic*cad).
730  & lt.rlu(0)) goto 340
731  ENDIF
732  ENDIF
733  ENDIF
734 
735 C...Continue loop over partons that may branch, until none left.
736  IF(igm.GE.0) k(im,1)=14
737  n=n+nep
738  nep=2
739  IF(n.GT.mstu(4)-mstu(32)-5) THEN
740  CALL luerrm(11,'(LUSHOW:) no more memory left in LUJETS')
741  IF(mstu(21).GE.1) n=ns
742  IF(mstu(21).GE.1) RETURN
743  ENDIF
744  goto 140
745 
746 C...Set information on imagined shower initiator.
747  380 IF(npa.GE.2) THEN
748  k(ns+1,1)=11
749  k(ns+1,2)=94
750  k(ns+1,3)=ip1
751  IF(ip2.GT.0.AND.ip2.LT.ip1) k(ns+1,3)=ip2
752  k(ns+1,4)=ns+2
753  k(ns+1,5)=ns+1+npa
754  iim=1
755  ELSE
756  iim=0
757  ENDIF
758 
759 C...Reconstruct string drawing information.
760  DO 390 i=ns+1+iim,n
761  IF(k(i,1).LE.10.AND.k(i,2).EQ.22) THEN
762  k(i,1)=1
763  ELSEIF(k(i,1).LE.10) THEN
764  k(i,4)=mstu(5)*(k(i,4)/mstu(5))
765  k(i,5)=mstu(5)*(k(i,5)/mstu(5))
766  ELSEIF(k(mod(k(i,4),mstu(5))+1,2).NE.22) THEN
767  id1=mod(k(i,4),mstu(5))
768  IF(k(i,2).GE.1.AND.k(i,2).LE.8) id1=mod(k(i,4),mstu(5))+1
769  id2=2*mod(k(i,4),mstu(5))+1-id1
770  k(i,4)=mstu(5)*(k(i,4)/mstu(5))+id1
771  k(i,5)=mstu(5)*(k(i,5)/mstu(5))+id2
772  k(id1,4)=k(id1,4)+mstu(5)*i
773  k(id1,5)=k(id1,5)+mstu(5)*id2
774  k(id2,4)=k(id2,4)+mstu(5)*id1
775  k(id2,5)=k(id2,5)+mstu(5)*i
776  ELSE
777  id1=mod(k(i,4),mstu(5))
778  id2=id1+1
779  k(i,4)=mstu(5)*(k(i,4)/mstu(5))+id1
780  k(i,5)=mstu(5)*(k(i,5)/mstu(5))+id1
781  k(id1,4)=k(id1,4)+mstu(5)*i
782  k(id1,5)=k(id1,5)+mstu(5)*i
783  k(id2,4)=0
784  k(id2,5)=0
785  ENDIF
786  390 CONTINUE
787 
788 C...Transformation from CM frame.
789  IF(npa.GE.2) THEN
790  bex=ps(1)/ps(4)
791  bey=ps(2)/ps(4)
792  bez=ps(3)/ps(4)
793  ga=ps(4)/ps(5)
794  gabep=ga*(ga*(bex*p(ipa(1),1)+bey*p(ipa(1),2)+bez*p(ipa(1),3))
795  & /(1.+ga)-p(ipa(1),4))
796  ELSE
797  bex=0.
798  bey=0.
799  bez=0.
800  gabep=0.
801  ENDIF
802  the=ulangl(p(ipa(1),3)+gabep*bez,sqrt((p(ipa(1),1)
803  &+gabep*bex)**2+(p(ipa(1),2)+gabep*bey)**2))
804  phi=ulangl(p(ipa(1),1)+gabep*bex,p(ipa(1),2)+gabep*bey)
805  IF(npa.EQ.3) THEN
806  chi=ulangl(cos(the)*cos(phi)*(p(ipa(2),1)+gabep*bex)+cos(the)*
807  & sin(phi)*(p(ipa(2),2)+gabep*bey)-sin(the)*(p(ipa(2),3)+gabep*
808  & bez),-sin(phi)*(p(ipa(2),1)+gabep*bex)+cos(phi)*(p(ipa(2),2)+
809  & gabep*bey))
810  CALL ludbrb(ns+1,n,0.,chi,0d0,0d0,0d0)
811  ENDIF
812  dbex=dble(bex)
813  dbey=dble(bey)
814  dbez=dble(bez)
815  CALL ludbrb(ns+1,n,the,phi,dbex,dbey,dbez)
816 
817 C...Decay vertex of shower.
818  DO 400 i=ns+1,n
819  DO 400 j=1,5
820  400 v(i,j)=v(ip1,j)
821 
822 C...Delete trivial shower, else connect initiators.
823  IF(n.EQ.ns+npa+iim) THEN
824  n=ns
825  ELSE
826  DO 410 ip=1,npa
827  k(ipa(ip),1)=14
828  k(ipa(ip),4)=k(ipa(ip),4)+ns+iim+ip
829  k(ipa(ip),5)=k(ipa(ip),5)+ns+iim+ip
830  k(ns+iim+ip,3)=ipa(ip)
831  IF(iim.EQ.1.AND.mstu(16).NE.2) k(ns+iim+ip,3)=ns+1
832  k(ns+iim+ip,4)=mstu(5)*ipa(ip)+k(ns+iim+ip,4)
833  410 k(ns+iim+ip,5)=mstu(5)*ipa(ip)+k(ns+iim+ip,5)
834  ENDIF
835 
836  RETURN
837  END