EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lshowr.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lshowr.F
1 
2 C#######################################################################
3 C
4 C The following routines for parton cascades were made together
5 C with M. Bengtsson and T. Sjostrand (Z. Phys. C37 (1988) 465,
6 C Nucl. Phys. B301 (1988) 554). Contain modifications of
7 C routines in PYTHIA 4.8 (Sjostrand, Bengtsson, CPC 46 (1987) 43).
8 C
9 C **********************************************************************
10 
11  SUBROUTINE lshowr(ICALL)
12 
13  IMPLICIT NONE
14 
15  INTEGER nlupdm,nplbuf
16  parameter(nlupdm=4000,nplbuf=5)
17  common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
18  INTEGER n,k
19  REAL p,v
20  SAVE /lujets/
21 
22 *
23 * to avoid variable conflictions, a second keep element is necessary
24 * with the same common block name (see LEPTO2)
25 *
26 
27  COMMON /leptou/ cut(14),lst(40),parl(30),
28  & xlp,ylp,w2lp,q2lp,ulp
29  REAL cut,parl,xlp,ylp,w2lp,q2lp,ulp
30  INTEGER lst
31  SAVE /leptou/
32 
33  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
34  INTEGER mstu,mstj
35  REAL paru,parj
36  SAVE /ludat1/
37 
38  COMMON /lboost/ dbeta(2,3),stheta(2),sphi(2),pb(5),phir
39  DOUBLE PRECISION dbeta
40  REAL stheta,sphi,pb,phir
41  SAVE /lboost/
42 
43  COMMON /pypara/ ipy(80),pypar(80),pyvar(80)
44  REAL pypar,pyvar
45  INTEGER ipy
46  SAVE /pypara/
47 
48 *
49 * to avoid variable conflictions, a second keep element is necessary
50 * with the same common block name (see LYPRO2)
51 *
52  COMMON /lyproc/ isub,kfl(3,2),x(2),sh,th,uh,q2,xsec(0:40)
53  REAL x,sh,th,uh,q2,xsec
54  INTEGER isub,kfl
55  SAVE /lyproc/
56 
57  COMMON /lyint1/ xq(2,-6:6),dsig(-6:6,-6:6,5),fsig(10,10,3)
58  REAL xq,dsig,fsig
59  SAVE /lyint1/
60 
61 
62 
63  INTEGER icall,i,j,ns,ifl,ifl1,ifl2,it,ipu1,ipu2
64  REAL xr,pma1,pma2,ulmass,qmax,ulangl,t1,t2
65 
66  DOUBLE PRECISION dpq2,dpb(3),dpa(3),dcthet,drobo(5)
67 
68  INTEGER ks(9,5)
69  REAL ps(9,5),robo(5),xpq(-6:6)
70  SAVE ks,ps
71 
72  IF(icall.EQ.0) THEN
73 C...Initialize cascade for each event, save event record in overall cms.
74  DO 10 i=1,9
75  DO 10 j=1,5
76  ks(i,j)=0
77  10 ps(i,j)=0.
78  DO 20 j=1,5
79  ks(1,j)=k(1,j)
80  ps(1,j)=p(1,j)
81  ks(2,j)=k(2,j)
82  ps(2,j)=p(2,j)
83  ks(5,j)=k(3,j)
84  ps(5,j)=p(3,j)
85  ks(7,j)=k(4,j)
86  20 ps(7,j)=p(4,j)
87  ks(5,3)=3
88  ks(7,1)=21
89  ks(7,3)=5
90 C CALL GULIST(0,2)
91  RETURN
92  ENDIF
93 
94 C CALL GULIST(1,2)
95 C...Apply parton cascade on QPM event.
96 C...Save incoming and outgoing quark as well as scattered lepton.
97  ks(6,1)=21
98  ks(6,2)=lst(25)
99  ks(6,3)=4
100  ks(8,1)=21
101  ks(8,2)=k(5,2)
102  ks(8,3)=6
103  ks(9,1)=0
104  ks(9,2)=k(4,2)
105  ks(9,3)=5
106  DO 110 j=1,5
107  ps(6,j)=0.
108  ps(8,j)=p(5,j)
109  110 ps(9,j)=p(4,j)
110  xr=xlp
111  dpq2=dble(q2lp)
112  pma1=0.
113  ps(6,5)=pma1
114  pma2=ps(8,5)
115  dpb(1)=0.5d0*(dpq2*(1d0/xr-1d0)+dble(ps(1,5))**2-
116  &ulmass(iabs(ks(7,2)))**2)/(ps(1,4)+ps(2,4))
117  dpb(2)=dsqrt(dpb(1)**2+dpq2)
118  dcthet=(dble(ps(2,4))*dpb(1)-dpq2/(2d0*xr))/(dble(ps(2,3))*
119  &dpb(2))
120  dpa(1)=(dpb(2)*dcthet)**2-dpb(1)**2
121  dpa(2)=dpq2-dble(pma1)**2+dble(pma2)**2
122  ps(6,4)=-(dpa(2)*dpb(1)-dpb(2)*dcthet*dsqrt(dpa(2)**2+4d0*
123  &dble(pma1)**2*dpa(1)))/(2d0*dpa(1))
124  ps(6,3)=-sqrt((ps(6,4)+pma1)*(ps(6,4)-pma1))
125 C...Partons with colour information in hadronic cms frame.
126  DO 120 i=10,26
127  DO 120 j=1,5
128  k(i,j)=0
129  p(i,j)=0.
130  120 v(i,j)=0.
131  ns=20
132  k(ns+1,1)=21
133  k(ns+1,2)=k(3,2)
134  k(ns+1,3)=3
135  k(ns+2,1)=-1
136  k(ns+2,3)=ns+1
137  k(ns+3,2)=ks(6,2)
138  DO 130 j=1,5
139  130 p(ns+1,j)=p(3,j)
140  k(ns+3,1)=13
141  k(ns+3,3)=2
142  p(ns+3,5)=0.
143  k(ns+4,1)=-1
144  k(ns+4,3)=ns+3
145  k(ns+3,4)=ns+5
146  k(ns+3,5)=ns+5
147  p(ns+4,3)=ns+5
148  p(ns+4,4)=ns+5
149  k(ns+5,1)=3
150  k(ns+5,3)=8
151  k(ns+5,2)=ks(8,2)
152  k(ns+6,1)=-1
153  k(ns+6,3)=ns+5
154  DO 140 j=1,4
155  p(ns+5,j)=p(5,j)
156  p(ns+3,j)=p(ns+5,j)-p(ns+1,j)
157  140 CONTINUE
158  p(ns+5,5)=pma2
159  p(ns+6,1)=ns+3
160  p(ns+6,2)=ns+3
161  k(ns+5,4)=(ns+3)*mstu(5)
162  k(ns+5,5)=(ns+3)*mstu(5)
163  n=ns+6
164 C CALL GULIST(2,2)
165 C...Copy saved record in overall cms to line 1 through 9.
166 C...Lines 1,2,5,6,7 in ep cms, 8,9 in hadronic cms
167  DO 150 i=1,9
168  DO 150 j=1,5
169  k(i,j)=ks(i,j)
170  150 p(i,j)=ps(i,j)
171 C CALL GULIST(3,2)
172 C...Scale for bremsstrahlung etc.
173  q2=q2lp
174  ipy(40)=8
175  ipy(47)=n
176 C...Save quantities for later use.
177  x(1)=1.
178  x(2)=xr
179  CALL lystfu(k(2,2),xr,q2lp,xpq)
180  DO 160 ifl=-6,6
181  160 xq(2,ifl)=xpq(ifl)
182  IF(lst(23).EQ.1) THEN
183  isub=39
184  ipy(11)=1
185  ELSEIF(lst(23).EQ.3) THEN
186  isub=39
187  ipy(11)=2
188  ELSEIF(lst(23).EQ.4) THEN
189  isub=39
190  ipy(11)=3
191  ELSEIF(lst(23).EQ.2) THEN
192  isub=40
193  ENDIF
194  kfl(2,1)=k(5,2)
195  ifl1=k(6,2)
196  ifl2=k(8,2)
197  kfl(2,2)=ifl1
198  kfl(1,1)=kfl(2,1)
199  kfl(1,2)=kfl(2,2)
200  IF(isub.EQ.39) kfl(3,1)=k(1,2)
201  IF(isub.EQ.40) kfl(3,1)=k(1,2)+isign(1,k(1,2))
202  kfl(3,2)=ifl2
203  pyvar(2)=(p(1,4)+p(2,4))**2
204  pyvar(1)=sqrt(pyvar(2))
205  pyvar(3)=p(1,5)
206  pyvar(4)=p(2,5)
207  pyvar(5)=p(1,3)
208  ipy(41)=k(1,2)
209  ipy(42)=k(2,2)
210  ipy(48)=0
211 
212 C...Generate timelike parton shower (if required)
213  IF(ipy(13).EQ.1) THEN
214  CALL lscale(1,qmax)
215  qmax=min(qmax,p(25,4))
216  CALL lushow(25,0,qmax)
217  ENDIF
218  it=25
219  IF(n.GE.27) it=27
220  ns=n
221 C CALL GULIST(4,2)
222 
223 C...Generate spacelike parton shower (if required)
224  ipu1=0
225  ipu2=23
226  IF(x(2)*(1.+(p(it,5)**2+pypar(22))/p(21,5)**2).GT.0.999) THEN
227  lst(21)=7
228  RETURN
229  ENDIF
230  IF(ipy(14).GE.1) THEN
231  CALL lysspa(ipu1,ipu2)
232  IF(lst(21).NE.0) RETURN
233  ENDIF
234  IF(n.EQ.ns) THEN
235  DO 220 i=ns+1,ns+4
236  DO 220 j=1,5
237  k(i,j)=0
238  p(i,j)=0.
239  220 v(i,j)=0.
240  k(ns+1,1)=11
241  k(ns+1,2)=kfl(2,1)
242  k(ns+1,3)=21
243  DO 230 j=1,5
244  230 p(ns+1,j)=p(21,j)
245  k(ns+2,1)=-1
246  k(ns+2,3)=ns+1
247  k(ns+3,1)=13
248  k(ns+3,2)=kfl(2,2)
249  k(ns+3,3)=23
250  k(ns+3,4)=23
251  k(ns+3,5)=23
252  p(ns+3,3)=(p(it,5)**2+q2)*(p(21,4)-p(21,3))/(2.*q2)
253  p(ns+3,4)=-p(ns+3,3)
254  k(ns+4,1)=-1
255  k(ns+4,3)=ns+3
256  p(ns+4,3)=23
257  p(ns+4,4)=23
258  p(24,1)=ns+3
259  p(24,2)=ns+3
260  k(23,4)=k(23,4)+(ns+3)*mstu(5)
261  k(23,5)=k(23,5)+(ns+3)*mstu(5)
262  ipu1=0
263  ipu2=ns+3
264  n=n+4
265  ENDIF
266 C CALL GULIST(5,2)
267 
268  IF(abs(p(it,1)).GT.0.1.OR.abs(p(it,2)).GT.0.1) THEN
269 C WRITE(6,*) 'Warning: non-zero pt on final shower initiator'
270 C WRITE(6,*) '0:',IT,K(IT,2),P(IT,1),P(IT,2),P(IT,3),P(IT,4),P(IT,5)
271  lst(21)=8
272  RETURN
273  ENDIF
274  p(it,1)=0.
275  p(it,2)=0.
276 
277 C...Rotate and boost outgoing parton shower
278  IF(n.GT.30) THEN
279  k(n+1,1)=0
280  DO 210 j=1,4
281  210 p(n+1,j)=p(ns+1,j)+p(ns+3,j)
282  IF(p(n+1,4).LE.1.01*p(it,5)) THEN
283  lst(21)=9
284  RETURN
285  ENDIF
286  robo(1)=ulangl(p(it,3),sqrt(p(it,1)**2+p(it,2)**2))
287  robo(2)=ulangl(p(it,1),p(it,2))
288  IF(abs(robo(1)).GT.0.001.OR.abs(robo(2)).GT.0.001) THEN
289  WRITE(6,*) '0:',it,k(it,2),p(it,1),p(it,2),p(it,3),p(it,4),p(it,5)
290  WRITE(6,*) ' ROBO(1-2)=',robo(1),robo(2)
291  ENDIF
292  CALL ludbrb(25,ns,0.,-robo(2),0.d0,0.d0,0.d0)
293  CALL ludbrb(25,ns,-robo(1),0.,0.d0,0.d0,0.d0)
294  drobo(5)=-(p(it,3)*p(it,4)-p(n+1,4)*sqrt(p(n+1,4)**2-
295  & p(it,4)**2+p(it,3)**2))/(p(it,3)**2+p(n+1,4)**2)
296  CALL ludbrb(25,ns,0.,0.,0.d0,0.d0,drobo(5))
297  robo(1)=ulangl(p(n+1,3),sqrt(p(n+1,1)**2+p(n+1,2)**2))
298  robo(2)=ulangl(p(n+1,1),p(n+1,2))
299  CALL ludbrb(25,ns,robo(1),robo(2),0.d0,0.d0,0.d0)
300  ENDIF
301 C CALL GULIST(6,2)
302 
303  q2=q2lp
304 C...Hadron remnant and primordial kt
305  ipy(47)=n
306  CALL lyremn(ipu1,ipu2)
307  IF(ipy(48).EQ.1) THEN
308  lst(21)=10
309  RETURN
310  ENDIF
311 C CALL GULIST(7,2)
312 
313 C...Transform line 1,2 and 5-7 to hadronic cms frame.
314  CALL ludbrb(1,2,0.,0.,-dbeta(2,1),-dbeta(2,2),-dbeta(2,3))
315  CALL ludbrb(1,2,-stheta(2),0.,0.d0,0.d0,0.d0)
316  CALL ludbrb(5,7,0.,0.,-dbeta(2,1),-dbeta(2,2),-dbeta(2,3))
317  CALL ludbrb(5,7,-stheta(2),0.,0.d0,0.d0,0.d0)
318 C CALL GULIST(8,2)
319 
320 C...Rearrange partons along strings
321  mstu(24)=0
322  CALL luprep(0)
323  IF(mstu(24).NE.0) THEN
324 C CALL GULIST(88,2)
325  IF(lst(3).GE.1) WRITE(6,*) ' LUPREP error MSTU(24)= ',mstu(24)
326  lst(21)=11
327  RETURN
328  ENDIF
329 C CALL GULIST(9,2)
330 
331 C...Clean up event record -> order:
332 C...1=inc. lepton; 2=inc. nucleon; 3=exch boson; 4=scat. lepton;
333 C...5=inc. parton before initial shower; 6=inc. quark at boson vertex
334 C...after shower; 7=scat. quark at boson vertex before final shower
335  lst(26)=7
336  DO 510 j=1,5
337  k(n+1,j)=k(4,j)
338  510 p(n+1,j)=p(4,j)
339  DO 520 j=1,5
340  k(3,j)=k(5,j)
341  p(3,j)=p(5,j)
342  k(4,j)=k(9,j)
343  p(4,j)=p(9,j)
344  k(5,j)=k(n+1,j)
345  p(5,j)=p(n+1,j)
346 C K(7,J)=K(8,J)
347 C P(7,J)=P(8,J)
348  k(6,j)=k(ns+3,j)
349  p(6,j)=p(ns+3,j)
350  k(7,j)=k(it,j)
351  p(7,j)=p(it,j)
352  520 CONTINUE
353  k(3,3)=1
354  k(4,3)=1
355  k(6,1)=21
356  k(6,3)=5
357  k(6,4)=0
358  k(6,5)=0
359  k(7,1)=21
360  k(7,3)=6
361  k(7,4)=0
362  k(7,5)=0
363 C...Activate line with scattered lepton.
364  k(4,1)=1
365 C...Deactivate obsolete lines 8, 9 and 21, NS+1 (extra lines with boson)
366  k(8,1)=0
367  k(9,1)=0
368  k(21,1)=0
369  IF(k(ns+1,2).EQ.k(3,2)) k(ns+1,1)=0
370 C...Zero irrelevant lines with K(I,1)<0
371  DO 540 i=1,n
372  IF(k(i,1).LT.0) THEN
373  DO 530 j=1,5
374  k(i,j)=0
375  530 p(i,j)=0.
376  ENDIF
377  540 CONTINUE
378 C CALL GULIST(10,2)
379 C...Delete internal parton lines, i.e. with K(I,1)=13,14
380  IF(mod(lst(4)/10,10).EQ.0) THEN
381  CALL ltimex(t1)
382  CALL luedit(14)
383  CALL ltimex(t2)
384 C CALL GULIST(11,2)
385  ENDIF
386 C...Delete empty lines
387  CALL ltimex(t1)
388  CALL luedit(12)
389  CALL ltimex(t2)
390 C CALL GULIST(12,2)
391 
392  RETURN
393  END