Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhisspa.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhisspa.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhisspa(IPU1,IPU2)
5 
6 C...Generates spacelike parton showers.
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  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
15  SAVE /pyhisubs/
16  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
17  SAVE /pyhipars/
18  common/pyhiint1/mint(400),vint(400)
19  SAVE /pyhiint1/
20  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
21  SAVE /pyhiint2/
22  common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
23  SAVE /pyhiint3/
24  dimension kfls(4),is(2),xs(2),zs(2),q2s(2),tevs(2),robo(5),
25  &xfs(2,-6:6),xfa(-6:6),xfb(-6:6),xfn(-6:6),wtap(-6:6),wtsf(-6:6),
26  &the2(2),alam(2),dq2(3),dpc(3),dpd(4),dpb(4)
27 
28 C...Calculate maximum virtuality and check that evolution possible.
29  ipus1=ipu1
30  ipus2=ipu2
31  isub=mint(1)
32  q2e=vint(52)
33  IF(iset(isub).EQ.1) THEN
34  q2e=q2e/parp(67)
35  ELSEIF(iset(isub).EQ.3.OR.iset(isub).EQ.4) THEN
36  q2e=pmas(23,1)**2
37  IF(isub.EQ.8.OR.isub.EQ.76.OR.isub.EQ.77) q2e=pmas(24,1)**2
38  ENDIF
39  tmax=log(parp(67)*parp(63)*q2e/parp(61)**2)
40  IF(parp(67)*q2e.LT.max(parp(62)**2,2.*parp(61)**2).OR.
41  &tmax.LT.0.2) RETURN
42 
43 C...Common constants and initial values. Save normal Lambda value.
44  xe0=2.*parp(65)/vint(1)
45  alams=paru(111)
46  paru(111)=parp(61)
47  ns=n
48  100 n=ns
49  DO 110 jt=1,2
50  kfls(jt)=mint(14+jt)
51  kfls(jt+2)=kfls(jt)
52  xs(jt)=vint(40+jt)
53  zs(jt)=1.
54  q2s(jt)=parp(67)*q2e
55  tevs(jt)=tmax
56  alam(jt)=parp(61)
57  the2(jt)=100.
58  DO 110 kfl=-6,6
59  110 xfs(jt,kfl)=xsfx(jt,kfl)
60  dsh=vint(44)
61  IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) dsh=vint(26)*vint(2)
62 
63 C...Pick up leg with highest virtuality.
64  120 n=n+1
65  jt=1
66  IF(n.GT.ns+1.AND.q2s(2).GT.q2s(1)) jt=2
67  kflb=kfls(jt)
68  xb=xs(jt)
69  DO 130 kfl=-6,6
70  130 xfb(kfl)=xfs(jt,kfl)
71  dshr=2d0*sqrt(dsh)
72  dshz=dsh/dble(zs(jt))
73  xe=max(xe0,xb*(1./(1.-parp(66))-1.))
74  IF(xb+xe.GE.0.999) THEN
75  q2b=0.
76  goto 220
77  ENDIF
78 
79 C...Maximum Q2 without or with Q2 ordering. Effective Lambda and n_f.
80  IF(mstp(62).LE.1) THEN
81  q2b=0.5*(1./zs(jt)+1.)*q2s(jt)+0.5*(1./zs(jt)-1.)*(q2s(3-jt)-
82  & sngl(dsh)+sqrt((sngl(dsh)+q2s(1)+q2s(2))**2+8.*q2s(1)*q2s(2)*
83  & zs(jt)/(1.-zs(jt))))
84  tevb=log(parp(63)*q2b/alam(jt)**2)
85  ELSE
86  q2b=q2s(jt)
87  tevb=tevs(jt)
88  ENDIF
89  alsdum=ulalps(parp(63)*q2b)
90  tevb=tevb+2.*log(alam(jt)/paru(117))
91  tevbsv=tevb
92  alam(jt)=paru(117)
93  b0=(33.-2.*mstu(118))/6.
94 
95 C...Calculate Altarelli-Parisi and structure function weights.
96  DO 140 kfl=-6,6
97  wtap(kfl)=0.
98  140 wtsf(kfl)=0.
99  IF(kflb.EQ.21) THEN
100  wtapq=16.*(1.-sqrt(xb+xe))/(3.*sqrt(xb))
101  DO 150 kfl=-mstp(54),mstp(54)
102  IF(kfl.EQ.0) wtap(kfl)=6.*log((1.-xb)/xe)
103  150 IF(kfl.NE.0) wtap(kfl)=wtapq
104  ELSE
105  wtap(0)=0.5*xb*(1./(xb+xe)-1.)
106  wtap(kflb)=8.*log((1.-xb)*(xb+xe)/xe)/3.
107  ENDIF
108  160 wtsum=0.
109  IF(kflb.NE.21) xfbo=xfb(kflb)
110  IF(kflb.EQ.21) xfbo=xfb(0)
111 C***************************************************************
112 C**********ERROR HAS OCCURED HERE
113  IF(xfbo.EQ.0.0) THEN
114  WRITE(mstu(11),1000)
115  WRITE(mstu(11),1001) kflb,xfb(kflb)
116  xfbo=0.00001
117  ENDIF
118 C****************************************************************
119  DO 170 kfl=-mstp(54),mstp(54)
120  wtsf(kfl)=xfb(kfl)/xfbo
121  170 wtsum=wtsum+wtap(kfl)*wtsf(kfl)
122  wtsum=max(0.0001,wtsum)
123 
124 C...Choose new t: fix alpha_s, alpha_s(Q2), alpha_s(k_T2).
125  180 IF(mstp(64).LE.0) THEN
126  tevb=tevb+log(rlu(0))*paru(2)/(paru(111)*wtsum)
127  ELSEIF(mstp(64).EQ.1) THEN
128  tevb=tevb*exp(max(-100.,log(rlu(0))*b0/wtsum))
129  ELSE
130  tevb=tevb*exp(max(-100.,log(rlu(0))*b0/(5.*wtsum)))
131  ENDIF
132  190 q2ref=alam(jt)**2*exp(tevb)
133  q2b=q2ref/parp(63)
134 
135 C...Evolution ended or select flavour for branching parton.
136  IF(q2b.LT.parp(62)**2) THEN
137  q2b=0.
138  ELSE
139  wtran=rlu(0)*wtsum
140  kfla=-mstp(54)-1
141  200 kfla=kfla+1
142  wtran=wtran-wtap(kfla)*wtsf(kfla)
143  IF(kfla.LT.mstp(54).AND.wtran.GT.0.) goto 200
144  IF(kfla.EQ.0) kfla=21
145 
146 C...Choose z value and corrective weight.
147  IF(kflb.EQ.21.AND.kfla.EQ.21) THEN
148  z=1./(1.+((1.-xb)/xb)*(xe/(1.-xb))**rlu(0))
149  wtz=(1.-z*(1.-z))**2
150  ELSEIF(kflb.EQ.21) THEN
151  z=xb/(1.-rlu(0)*(1.-sqrt(xb+xe)))**2
152  wtz=0.5*(1.+(1.-z)**2)*sqrt(z)
153  ELSEIF(kfla.EQ.21) THEN
154  z=xb*(1.+rlu(0)*(1./(xb+xe)-1.))
155  wtz=1.-2.*z*(1.-z)
156  ELSE
157  z=1.-(1.-xb)*(xe/((xb+xe)*(1.-xb)))**rlu(0)
158  wtz=0.5*(1.+z**2)
159  ENDIF
160 
161 C...Option with resummation of soft gluon emission as effective z shift.
162  IF(mstp(65).GE.1) THEN
163  rsoft=6.
164  IF(kflb.NE.21) rsoft=8./3.
165  z=z*(tevb/tevs(jt))**(rsoft*xe/((xb+xe)*b0))
166  IF(z.LE.xb) goto 180
167  ENDIF
168 
169 C...Option with alpha_s(k_T2)Q2): demand k_T2 > cutoff, reweight.
170  IF(mstp(64).GE.2) THEN
171  IF((1.-z)*q2b.LT.parp(62)**2) goto 180
172  alprat=tevb/(tevb+log(1.-z))
173  IF(alprat.LT.5.*rlu(0)) goto 180
174  IF(alprat.GT.5.) wtz=wtz*alprat/5.
175  ENDIF
176 
177 C...Option with angular ordering requirement.
178  IF(mstp(62).GE.3) THEN
179  the2t=(4.*z**2*q2b)/(vint(2)*(1.-z)*xb**2)
180  IF(the2t.GT.the2(jt)) goto 180
181  ENDIF
182 
183 C...Weighting with new structure functions.
184  CALL pyhistfu(mint(10+jt),xb,q2ref,xfn,jt)
185  IF(kflb.NE.21) xfbn=xfn(kflb)
186  IF(kflb.EQ.21) xfbn=xfn(0)
187  IF(xfbn.LT.1e-20) THEN
188  IF(kfla.EQ.kflb) THEN
189  tevb=tevbsv
190  wtap(kflb)=0.
191  goto 160
192  ELSEIF(tevbsv-tevb.GT.0.2) THEN
193  tevb=0.5*(tevbsv+tevb)
194  goto 190
195  ELSE
196  xfbn=1e-10
197  ENDIF
198  ENDIF
199  DO 210 kfl=-mstp(54),mstp(54)
200  210 xfb(kfl)=xfn(kfl)
201  xa=xb/z
202  CALL pyhistfu(mint(10+jt),xa,q2ref,xfa,jt)
203  IF(kfla.NE.21) xfan=xfa(kfla)
204  IF(kfla.EQ.21) xfan=xfa(0)
205  IF(xfan.LT.1e-20) goto 160
206  IF(kfla.NE.21) wtsfa=wtsf(kfla)
207  IF(kfla.EQ.21) wtsfa=wtsf(0)
208  IF(wtz*xfan/xfbn.LT.rlu(0)*wtsfa) goto 160
209  ENDIF
210 
211 C...Define two hard scatterers in their CM-frame.
212  220 IF(n.EQ.ns+2) THEN
213  dq2(jt)=q2b
214  dplcm=sqrt((dsh+dq2(1)+dq2(2))**2-4d0*dq2(1)*dq2(2))/dshr
215  DO 240 jr=1,2
216  i=ns+jr
217  IF(jr.EQ.1) ipo=ipus1
218  IF(jr.EQ.2) ipo=ipus2
219  DO 230 j=1,5
220  k(i,j)=0
221  p(i,j)=0.
222  230 v(i,j)=0.
223  k(i,1)=14
224  k(i,2)=kfls(jr+2)
225  k(i,4)=ipo
226  k(i,5)=ipo
227  p(i,3)=dplcm*(-1)**(jr+1)
228  p(i,4)=(dsh+dq2(3-jr)-dq2(jr))/dshr
229  p(i,5)=-sqrt(sngl(dq2(jr)))
230  k(ipo,1)=14
231  k(ipo,3)=i
232  k(ipo,4)=mod(k(ipo,4),mstu(5))+mstu(5)*i
233  240 k(ipo,5)=mod(k(ipo,5),mstu(5))+mstu(5)*i
234 
235 C...Find maximum allowed mass of timelike parton.
236  ELSEIF(n.GT.ns+2) THEN
237  jr=3-jt
238  dq2(3)=q2b
239  dpc(1)=p(is(1),4)
240  dpc(2)=p(is(2),4)
241  dpc(3)=0.5*(abs(p(is(1),3))+abs(p(is(2),3)))
242  dpd(1)=dsh+dq2(jr)+dq2(jt)
243  dpd(2)=dshz+dq2(jr)+dq2(3)
244  dpd(3)=sqrt(dpd(1)**2-4d0*dq2(jr)*dq2(jt))
245  dpd(4)=sqrt(dpd(2)**2-4d0*dq2(jr)*dq2(3))
246  ikin=0
247  IF(q2s(jr).GE.(0.5*parp(62))**2.AND.dpd(1)-dpd(3).GE.
248  & 1d-10*dpd(1)) ikin=1
249  IF(ikin.EQ.0) dmsma=(dq2(jt)/dble(zs(jt))-dq2(3))*(dsh/
250  & (dsh+dq2(jt))-dsh/(dshz+dq2(3)))
251  IF(ikin.EQ.1) dmsma=(dpd(1)*dpd(2)-dpd(3)*dpd(4))/(2.*
252  & dq2(jr))-dq2(jt)-dq2(3)
253 
254 C...Generate timelike parton shower (if required).
255  it=n
256  DO 250 j=1,5
257  k(it,j)=0
258  p(it,j)=0.
259  250 v(it,j)=0.
260  k(it,1)=3
261  k(it,2)=21
262  IF(kflb.EQ.21.AND.kfls(jt+2).NE.21) k(it,2)=-kfls(jt+2)
263  IF(kflb.NE.21.AND.kfls(jt+2).EQ.21) k(it,2)=kflb
264  p(it,5)=ulmass(k(it,2))
265  IF(sngl(dmsma).LE.p(it,5)**2) goto 100
266  IF(mstp(63).GE.1) THEN
267  p(it,4)=(dshz-dsh-p(it,5)**2)/dshr
268  p(it,3)=sqrt(p(it,4)**2-p(it,5)**2)
269  IF(mstp(63).EQ.1) THEN
270  q2tim=dmsma
271  ELSEIF(mstp(63).EQ.2) THEN
272  q2tim=min(sngl(dmsma),parp(71)*q2s(jt))
273  ELSE
274 C'''Here remains to introduce angular ordering in first branching.
275  q2tim=dmsma
276  ENDIF
277  CALL lushow(it,0,sqrt(q2tim))
278  IF(n.GE.it+1) p(it,5)=p(it+1,5)
279  ENDIF
280 
281 C...Reconstruct kinematics of branching: timelike parton shower.
282  dms=p(it,5)**2
283  IF(ikin.EQ.0) dpt2=(dmsma-dms)*(dshz+dq2(3))/(dsh+dq2(jt))
284  IF(ikin.EQ.1) dpt2=(dmsma-dms)*(0.5*dpd(1)*dpd(2)+0.5*dpd(3)*
285  & dpd(4)-dq2(jr)*(dq2(jt)+dq2(3)+dms))/(4.*dsh*dpc(3)**2)
286  IF(dpt2.LT.0.) goto 100
287  dpb(1)=(0.5*dpd(2)-dpc(jr)*(dshz+dq2(jr)-dq2(jt)-dms)/
288  & dshr)/dpc(3)-dpc(3)
289  p(it,1)=sqrt(sngl(dpt2))
290  p(it,3)=dpb(1)*(-1)**(jt+1)
291  p(it,4)=(dshz-dsh-dms)/dshr
292  IF(n.GE.it+1) THEN
293  dpb(1)=sqrt(dpb(1)**2+dpt2)
294  dpb(2)=sqrt(dpb(1)**2+dms)
295  dpb(3)=p(it+1,3)
296  dpb(4)=sqrt(dpb(3)**2+dms)
297  dbez=(dpb(4)*dpb(1)-dpb(3)*dpb(2))/(dpb(4)*dpb(2)-dpb(3)*
298  & dpb(1))
299  CALL ludbrb(it+1,n,0.,0.,0d0,0d0,dbez)
300  the=ulangl(p(it,3),p(it,1))
301  CALL ludbrb(it+1,n,the,0.,0d0,0d0,0d0)
302  ENDIF
303 
304 C...Reconstruct kinematics of branching: spacelike parton.
305  DO 260 j=1,5
306  k(n+1,j)=0
307  p(n+1,j)=0.
308  260 v(n+1,j)=0.
309  k(n+1,1)=14
310  k(n+1,2)=kflb
311  p(n+1,1)=p(it,1)
312  p(n+1,3)=p(it,3)+p(is(jt),3)
313  p(n+1,4)=p(it,4)+p(is(jt),4)
314  p(n+1,5)=-sqrt(sngl(dq2(3)))
315 
316 C...Define colour flow of branching.
317  k(is(jt),3)=n+1
318  k(it,3)=n+1
319  id1=it
320  IF((k(n+1,2).GT.0.AND.k(n+1,2).NE.21.AND.k(id1,2).GT.0.AND.
321  & k(id1,2).NE.21).OR.(k(n+1,2).LT.0.AND.k(id1,2).EQ.21).OR.
322  & (k(n+1,2).EQ.21.AND.k(id1,2).EQ.21.AND.rlu(0).GT.0.5).OR.
323  & (k(n+1,2).EQ.21.AND.k(id1,2).LT.0)) id1=is(jt)
324  id2=it+is(jt)-id1
325  k(n+1,4)=k(n+1,4)+id1
326  k(n+1,5)=k(n+1,5)+id2
327  k(id1,4)=k(id1,4)+mstu(5)*(n+1)
328  k(id1,5)=k(id1,5)+mstu(5)*id2
329  k(id2,4)=k(id2,4)+mstu(5)*id1
330  k(id2,5)=k(id2,5)+mstu(5)*(n+1)
331  n=n+1
332 
333 C...Boost to new CM-frame.
334  CALL ludbrb(ns+1,n,0.,0.,-dble((p(n,1)+p(is(jr),1))/(p(n,4)+
335  & p(is(jr),4))),0d0,-dble((p(n,3)+p(is(jr),3))/(p(n,4)+
336  & p(is(jr),4))))
337  ir=n+(jt-1)*(is(1)-n)
338  CALL ludbrb(ns+1,n,-ulangl(p(ir,3),p(ir,1)),paru(2)*rlu(0),
339  & 0d0,0d0,0d0)
340  ENDIF
341 
342 C...Save quantities, loop back.
343  is(jt)=n
344  q2s(jt)=q2b
345  dq2(jt)=q2b
346  IF(mstp(62).GE.3) the2(jt)=the2t
347  dsh=dshz
348  IF(q2b.GE.(0.5*parp(62))**2) THEN
349  kfls(jt+2)=kfls(jt)
350  kfls(jt)=kfla
351  xs(jt)=xa
352  zs(jt)=z
353  DO 270 kfl=-6,6
354  270 xfs(jt,kfl)=xfa(kfl)
355  tevs(jt)=tevb
356  ELSE
357  IF(jt.EQ.1) ipu1=n
358  IF(jt.EQ.2) ipu2=n
359  ENDIF
360  IF(n.GT.mstu(4)-mstu(32)-10) THEN
361  CALL luerrm(11,'(PYHISSPA:) no more memory left in LUJETS')
362  IF(mstu(21).GE.1) n=ns
363  IF(mstu(21).GE.1) RETURN
364  ENDIF
365  IF(max(q2s(1),q2s(2)).GE.(0.5*parp(62))**2.OR.n.LE.ns+1) goto 120
366 
367 C...Boost hard scattering partons to frame of shower initiators.
368  DO 280 j=1,3
369  280 robo(j+2)=(p(ns+1,j)+p(ns+2,j))/(p(ns+1,4)+p(ns+2,4))
370  DO 290 j=1,5
371  290 p(n+2,j)=p(ns+1,j)
372  robot=robo(3)**2+robo(4)**2+robo(5)**2
373  IF(robot.GE.0.999999) THEN
374  robot=1.00001*sqrt(robot)
375  robo(3)=robo(3)/robot
376  robo(4)=robo(4)/robot
377  robo(5)=robo(5)/robot
378  ENDIF
379  CALL ludbrb(n+2,n+2,0.,0.,-dble(robo(3)),-dble(robo(4)),
380  &-dble(robo(5)))
381  robo(2)=ulangl(p(n+2,1),p(n+2,2))
382  robo(1)=ulangl(p(n+2,3),sqrt(p(n+2,1)**2+p(n+2,2)**2))
383  CALL ludbrb(mint(83)+5,ns,robo(1),robo(2),dble(robo(3)),
384  &dble(robo(4)),dble(robo(5)))
385 
386 C...Store user information. Reset Lambda value.
387  k(ipu1,3)=mint(83)+3
388  k(ipu2,3)=mint(83)+4
389  DO 300 jt=1,2
390  mint(12+jt)=kfls(jt)
391  300 vint(140+jt)=xs(jt)
392  paru(111)=alams
393  1000 FORMAT(5x,'structure function has a zero point here')
394  1001 FORMAT(5x,'xf(x,i=',i5,')=',f10.5)
395 
396  RETURN
397  END