Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luprep.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luprep.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luprep(IP)
5 
6 C...Purpose: to rearrange partons along strings, to allow small systems
7 C...to collapse into one or two particles and to check flavours.
8  IMPLICIT DOUBLE PRECISION(d)
9  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
10  SAVE /lujets/
11  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
12  SAVE /ludat1/
13  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
14  SAVE /ludat2/
15  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
16  SAVE /ludat3/
17  dimension dps(5),dpc(5),ue(3)
18 
19 C...Rearrange parton shower product listing along strings: begin loop.
20  i1=n
21  DO 130 mqgst=1,2
22  DO 120 i=max(1,ip),n
23  IF(k(i,1).NE.3) goto 120
24  kc=lucomp(k(i,2))
25  IF(kc.EQ.0) goto 120
26  kq=kchg(kc,2)
27  IF(kq.EQ.0.OR.(mqgst.EQ.1.AND.kq.EQ.2)) goto 120
28 
29 C...Pick up loose string end.
30  kcs=4
31  IF(kq*isign(1,k(i,2)).LT.0) kcs=5
32  ia=i
33  nstp=0
34  100 nstp=nstp+1
35  IF(nstp.GT.4*n) THEN
36  CALL luerrm(14,'(LUPREP:) caught in infinite loop')
37  RETURN
38  ENDIF
39 
40 C...Copy undecayed parton.
41  IF(k(ia,1).EQ.3) THEN
42  IF(i1.GE.mstu(4)-mstu(32)-5) THEN
43  CALL luerrm(11,'(LUPREP:) no more memory left in LUJETS')
44  RETURN
45  ENDIF
46  i1=i1+1
47  k(i1,1)=2
48  IF(nstp.GE.2.AND.iabs(k(ia,2)).NE.21) k(i1,1)=1
49  k(i1,2)=k(ia,2)
50  k(i1,3)=ia
51  k(i1,4)=0
52  k(i1,5)=0
53  DO 110 j=1,5
54  p(i1,j)=p(ia,j)
55  110 v(i1,j)=v(ia,j)
56  k(ia,1)=k(ia,1)+10
57  IF(k(i1,1).EQ.1) goto 120
58  ENDIF
59 
60 C...Go to next parton in colour space.
61  ib=ia
62  IF(mod(k(ib,kcs)/mstu(5)**2,2).EQ.0.AND.mod(k(ib,kcs),mstu(5)).
63  &ne.0) THEN
64  ia=mod(k(ib,kcs),mstu(5))
65  k(ib,kcs)=k(ib,kcs)+mstu(5)**2
66  mrev=0
67  ELSE
68  IF(k(ib,kcs).GE.2*mstu(5)**2.OR.mod(k(ib,kcs)/mstu(5),mstu(5)).
69  & eq.0) kcs=9-kcs
70  ia=mod(k(ib,kcs)/mstu(5),mstu(5))
71  k(ib,kcs)=k(ib,kcs)+2*mstu(5)**2
72  mrev=1
73  ENDIF
74  IF(ia.LE.0.OR.ia.GT.n) THEN
75  CALL luerrm(12,'(LUPREP:) colour rearrangement failed')
76  RETURN
77  ENDIF
78  IF(mod(k(ia,4)/mstu(5),mstu(5)).EQ.ib.OR.mod(k(ia,5)/mstu(5),
79  &mstu(5)).EQ.ib) THEN
80  IF(mrev.EQ.1) kcs=9-kcs
81  IF(mod(k(ia,kcs)/mstu(5),mstu(5)).NE.ib) kcs=9-kcs
82  k(ia,kcs)=k(ia,kcs)+2*mstu(5)**2
83  ELSE
84  IF(mrev.EQ.0) kcs=9-kcs
85  IF(mod(k(ia,kcs),mstu(5)).NE.ib) kcs=9-kcs
86  k(ia,kcs)=k(ia,kcs)+mstu(5)**2
87  ENDIF
88  IF(ia.NE.i) goto 100
89  k(i1,1)=1
90  120 CONTINUE
91  130 CONTINUE
92  n=i1
93 
94 C...Find lowest-mass colour singlet jet system, OK if above threshold.
95  IF(mstj(14).LE.0) goto 320
96  ns=n
97  140 nsin=n-ns
98  pdm=1.+parj(32)
99  ic=0
100  DO 190 i=max(1,ip),ns
101  IF(k(i,1).NE.1.AND.k(i,1).NE.2) THEN
102  ELSEIF(k(i,1).EQ.2.AND.ic.EQ.0) THEN
103  nsin=nsin+1
104  ic=i
105  DO 150 j=1,4
106  150 dps(j)=p(i,j)
107  mstj(93)=1
108  dps(5)=ulmass(k(i,2))
109  ELSEIF(k(i,1).EQ.2) THEN
110  DO 160 j=1,4
111  160 dps(j)=dps(j)+p(i,j)
112  ELSEIF(ic.NE.0.AND.kchg(lucomp(k(i,2)),2).NE.0) THEN
113  DO 170 j=1,4
114  170 dps(j)=dps(j)+p(i,j)
115  mstj(93)=1
116  dps(5)=dps(5)+ulmass(k(i,2))
117  pd=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))-dps(5)
118  IF(pd.LT.pdm) THEN
119  pdm=pd
120  DO 180 j=1,5
121  180 dpc(j)=dps(j)
122  ic1=ic
123  ic2=i
124  ENDIF
125  ic=0
126  ELSE
127  nsin=nsin+1
128  ENDIF
129  190 CONTINUE
130  IF(pdm.GE.parj(32)) goto 320
131 
132 C...Fill small-mass system as cluster.
133  nsav=n
134  pecm=sqrt(max(0d0,dpc(4)**2-dpc(1)**2-dpc(2)**2-dpc(3)**2))
135  k(n+1,1)=11
136  k(n+1,2)=91
137  k(n+1,3)=ic1
138  k(n+1,4)=n+2
139  k(n+1,5)=n+3
140  p(n+1,1)=dpc(1)
141  p(n+1,2)=dpc(2)
142  p(n+1,3)=dpc(3)
143  p(n+1,4)=dpc(4)
144  p(n+1,5)=pecm
145 
146 C...Form two particles from flavours of lowest-mass system, if feasible.
147  k(n+2,1)=1
148  k(n+3,1)=1
149  IF(mstu(16).NE.2) THEN
150  k(n+2,3)=n+1
151  k(n+3,3)=n+1
152  ELSE
153  k(n+2,3)=ic1
154  k(n+3,3)=ic2
155  ENDIF
156  k(n+2,4)=0
157  k(n+3,4)=0
158  k(n+2,5)=0
159  k(n+3,5)=0
160  IF(iabs(k(ic1,2)).NE.21) THEN
161  kc1=lucomp(k(ic1,2))
162  kc2=lucomp(k(ic2,2))
163  IF(kc1.EQ.0.OR.kc2.EQ.0) goto 320
164  kq1=kchg(kc1,2)*isign(1,k(ic1,2))
165  kq2=kchg(kc2,2)*isign(1,k(ic2,2))
166  IF(kq1+kq2.NE.0) goto 320
167  200 CALL lukfdi(k(ic1,2),0,kfln,k(n+2,2))
168  CALL lukfdi(k(ic2,2),-kfln,kfldmp,k(n+3,2))
169  IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 200
170  ELSE
171  IF(iabs(k(ic2,2)).NE.21) goto 320
172  210 CALL lukfdi(1+int((2.+parj(2))*rlu(0)),0,kfln,kfdmp)
173  CALL lukfdi(kfln,0,kflm,k(n+2,2))
174  CALL lukfdi(-kfln,-kflm,kfldmp,k(n+3,2))
175  IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 210
176  ENDIF
177  p(n+2,5)=ulmass(k(n+2,2))
178  p(n+3,5)=ulmass(k(n+3,2))
179  IF(p(n+2,5)+p(n+3,5)+parj(64).GE.pecm.AND.nsin.EQ.1) goto 320
180  IF(p(n+2,5)+p(n+3,5)+parj(64).GE.pecm) goto 260
181 
182 C...Perform two-particle decay of jet system, if possible.
183  IF(pecm.GE.0.02*dpc(4)) THEN
184  pa=sqrt((pecm**2-(p(n+2,5)+p(n+3,5))**2)*(pecm**2-
185  & (p(n+2,5)-p(n+3,5))**2))/(2.*pecm)
186  ue(3)=2.*rlu(0)-1.
187  phi=paru(2)*rlu(0)
188  ue(1)=sqrt(1.-ue(3)**2)*cos(phi)
189  ue(2)=sqrt(1.-ue(3)**2)*sin(phi)
190  DO 220 j=1,3
191  p(n+2,j)=pa*ue(j)
192  220 p(n+3,j)=-pa*ue(j)
193  p(n+2,4)=sqrt(pa**2+p(n+2,5)**2)
194  p(n+3,4)=sqrt(pa**2+p(n+3,5)**2)
195  CALL ludbrb(n+2,n+3,0.,0.,dpc(1)/dpc(4),dpc(2)/dpc(4),
196  & dpc(3)/dpc(4))
197  ELSE
198  np=0
199  DO 230 i=ic1,ic2
200  230 IF(k(i,1).EQ.1.OR.k(i,1).EQ.2) np=np+1
201  ha=p(ic1,4)*p(ic2,4)-p(ic1,1)*p(ic2,1)-p(ic1,2)*p(ic2,2)-
202  & p(ic1,3)*p(ic2,3)
203  IF(np.GE.3.OR.ha.LE.1.25*p(ic1,5)*p(ic2,5)) goto 260
204  hd1=0.5*(p(n+2,5)**2-p(ic1,5)**2)
205  hd2=0.5*(p(n+3,5)**2-p(ic2,5)**2)
206  hr=sqrt(max(0.,((ha-hd1-hd2)**2-(p(n+2,5)*p(n+3,5))**2)/
207  & (ha**2-(p(ic1,5)*p(ic2,5))**2)))-1.
208  hc=p(ic1,5)**2+2.*ha+p(ic2,5)**2
209  hk1=((p(ic2,5)**2+ha)*hr+hd1-hd2)/hc
210  hk2=((p(ic1,5)**2+ha)*hr+hd2-hd1)/hc
211  DO 240 j=1,4
212  p(n+2,j)=(1.+hk1)*p(ic1,j)-hk2*p(ic2,j)
213  240 p(n+3,j)=(1.+hk2)*p(ic2,j)-hk1*p(ic1,j)
214  ENDIF
215  DO 250 j=1,4
216  v(n+1,j)=v(ic1,j)
217  v(n+2,j)=v(ic1,j)
218  250 v(n+3,j)=v(ic2,j)
219  v(n+1,5)=0.
220  v(n+2,5)=0.
221  v(n+3,5)=0.
222  n=n+3
223  goto 300
224 
225 C...Else form one particle from the flavours available, if possible.
226  260 k(n+1,5)=n+2
227  IF(iabs(k(ic1,2)).GT.100.AND.iabs(k(ic2,2)).GT.100) THEN
228  goto 320
229  ELSEIF(iabs(k(ic1,2)).NE.21) THEN
230  CALL lukfdi(k(ic1,2),k(ic2,2),kfldmp,k(n+2,2))
231  ELSE
232  kfln=1+int((2.+parj(2))*rlu(0))
233  CALL lukfdi(kfln,-kfln,kfldmp,k(n+2,2))
234  ENDIF
235  IF(k(n+2,2).EQ.0) goto 260
236  p(n+2,5)=ulmass(k(n+2,2))
237 
238 C...Find parton/particle which combines to largest extra mass.
239  ir=0
240  ha=0.
241  DO 280 mcomb=1,3
242  IF(ir.NE.0) goto 280
243  DO 270 i=max(1,ip),n
244  IF(k(i,1).LE.0.OR.k(i,1).GT.10.OR.(i.GE.ic1.AND.i.LE.ic2.
245  &and.k(i,1).GE.1.AND.k(i,1).LE.2)) goto 270
246  IF(mcomb.EQ.1) kci=lucomp(k(i,2))
247  IF(mcomb.EQ.1.AND.kci.EQ.0) goto 270
248  IF(mcomb.EQ.1.AND.kchg(kci,2).EQ.0.AND.i.LE.ns) goto 270
249  IF(mcomb.EQ.2.AND.iabs(k(i,2)).GT.10.AND.iabs(k(i,2)).LE.100)
250  &goto 270
251  hcr=dpc(4)*p(i,4)-dpc(1)*p(i,1)-dpc(2)*p(i,2)-dpc(3)*p(i,3)
252  IF(hcr.GT.ha) THEN
253  ir=i
254  ha=hcr
255  ENDIF
256  270 CONTINUE
257  280 CONTINUE
258 
259 C...Shuffle energy and momentum to put new particle on mass shell.
260  hb=pecm**2+ha
261  hc=p(n+2,5)**2+ha
262  hd=p(ir,5)**2+ha
263 C******************CHANGES BY HIJING************
264  hk2=0.0
265  IF(ha**2-(pecm*p(ir,5))**2.EQ.0.0.OR.hb+hd.EQ.0.0) go to 285
266 C******************
267  hk2=0.5*(hb*sqrt(((hb+hc)**2-4.*(hb+hd)*p(n+2,5)**2)/
268  &(ha**2-(pecm*p(ir,5))**2))-(hb+hc))/(hb+hd)
269  285 hk1=(0.5*(p(n+2,5)**2-pecm**2)+hd*hk2)/hb
270  DO 290 j=1,4
271  p(n+2,j)=(1.+hk1)*dpc(j)-hk2*p(ir,j)
272  p(ir,j)=(1.+hk2)*p(ir,j)-hk1*dpc(j)
273  v(n+1,j)=v(ic1,j)
274  290 v(n+2,j)=v(ic1,j)
275  v(n+1,5)=0.
276  v(n+2,5)=0.
277  n=n+2
278 
279 C...Mark collapsed system and store daughter pointers. Iterate.
280  300 DO 310 i=ic1,ic2
281  IF((k(i,1).EQ.1.OR.k(i,1).EQ.2).AND.kchg(lucomp(k(i,2)),2).NE.0)
282  &THEN
283  k(i,1)=k(i,1)+10
284  IF(mstu(16).NE.2) THEN
285  k(i,4)=nsav+1
286  k(i,5)=nsav+1
287  ELSE
288  k(i,4)=nsav+2
289  k(i,5)=n
290  ENDIF
291  ENDIF
292  310 CONTINUE
293  IF(n.LT.mstu(4)-mstu(32)-5) goto 140
294 
295 C...Check flavours and invariant masses in parton systems.
296  320 np=0
297  kfn=0
298  kqs=0
299  DO 330 j=1,5
300  330 dps(j)=0.
301  DO 360 i=max(1,ip),n
302  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 360
303  kc=lucomp(k(i,2))
304  IF(kc.EQ.0) goto 360
305  kq=kchg(kc,2)*isign(1,k(i,2))
306  IF(kq.EQ.0) goto 360
307  np=np+1
308  IF(kq.NE.2) THEN
309  kfn=kfn+1
310  kqs=kqs+kq
311  mstj(93)=1
312  dps(5)=dps(5)+ulmass(k(i,2))
313  ENDIF
314  DO 340 j=1,4
315  340 dps(j)=dps(j)+p(i,j)
316  IF(k(i,1).EQ.1) THEN
317  IF(np.NE.1.AND.(kfn.EQ.1.OR.kfn.GE.3.OR.kqs.NE.0)) CALL
318  & luerrm(2,'(LUPREP:) unphysical flavour combination')
319  IF(np.NE.1.AND.dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2.LT.
320  & (0.9*parj(32)+dps(5))**2) CALL luerrm(3,
321  & '(LUPREP:) too small mass in jet system')
322  np=0
323  kfn=0
324  kqs=0
325  DO 350 j=1,5
326  350 dps(j)=0.
327  ENDIF
328  360 CONTINUE
329 
330  RETURN
331  END