8 IMPLICIT DOUBLE PRECISION(d)
9 common/lujets/
n,
k(9000,5),
p(9000,5),
v(9000,5)
11 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
13 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
15 dimension dps(5),psi(4),nfi(3),nfl(3),ifet(3),kflf(3),
16 &kflo(2),pxo(2),pyo(2),wo(2)
26 IF(
i.GT.min(
n,mstu(4)-mstu(32)))
THEN
27 CALL
luerrm(12,
'(LUINDF:) failed to reconstruct jet system')
28 IF(mstu(21).GE.1)
RETURN
30 IF(
k(
i,1).NE.1.AND.
k(
i,1).NE.2) goto 110
33 kq=kchg(kc,2)*isign(1,
k(
i,2))
36 IF(kq.NE.2) kqsum=kqsum+kq
40 120 dps(
j)=dps(
j)+
p(
i,
j)
42 IF(
k(
i,1).EQ.2.OR.(mstj(3).LE.5.AND.
n.GT.
i.AND.
43 &
k(
i+1,1).EQ.2)) goto 110
44 IF(njet.NE.1.AND.kqsum.NE.0)
THEN
45 CALL
luerrm(12,
'(LUINDF:) unphysical flavour combination')
46 IF(mstu(21).GE.1)
RETURN
50 IF(njet.NE.1) CALL ludbrb(nsav+1,nsav+njet,0.,0.,-dps(1)/dps(4),
51 &-dps(2)/dps(4),-dps(3)/dps(4))
55 DO 140
i=nsav+1,nsav+njet
59 nfi(kfa)=nfi(kfa)+isign(1,
k(
i,2))
60 ELSEIF(kfa.GT.1000)
THEN
63 IF(kfla.LE.3) nfi(kfla)=nfi(kfla)+isign(1,
k(
i,2))
64 IF(kflb.LE.3) nfi(kflb)=nfi(kflb)+isign(1,
k(
i,2))
73 CALL
luerrm(14,
'(LUINDF:) caught in infinite loop')
74 IF(mstu(21).GE.1)
RETURN
82 DO 230 ip1=nsav+1,nsav+njet
88 IF(kflh.GT.10) kflh=mod(kflh/1000,10)
90 wf=
p(ip1,4)+sqrt(
p(ip1,1)**2+
p(ip1,2)**2+
p(ip1,3)**2)
93 170
IF(iabs(
k(ip1,2)).NE.21)
THEN
96 CALL
luptdi(0,pxo(1),pyo(1))
100 ELSEIF(mstj(2).LE.2)
THEN
102 IF(mstj(2).EQ.2) mstj(91)=1
103 kflo(1)=int(1.+(2.+parj(2))*
rlu(0))*(-1)**int(
rlu(0)+0.5)
104 CALL
luptdi(0,pxo(1),pyo(1))
111 IF(mstj(2).EQ.4) mstj(91)=1
112 kflo(1)=int(1.+(2.+parj(2))*
rlu(0))*(-1)**int(
rlu(0)+0.5)
114 CALL
luptdi(0,pxo(1),pyo(1))
117 wo(1)=wf*
rlu(0)**(1./3.)
132 IF(
i.GE.mstu(4)-mstu(32)-njet-5)
THEN
133 CALL
luerrm(11,
'(LUINDF:) no more memory left in LUJETS')
134 IF(mstu(21).GE.1)
RETURN
142 IF(
k(
i,2).EQ.0) goto 180
143 IF(mstj(12).GE.3.AND.irank.EQ.1.AND.iabs(kfl1).LE.10.AND.
144 &iabs(kfl2).GT.10)
THEN
145 IF(
rlu(0).GT.parj(19)) goto 200
153 pr=
p(
i,5)**2+
p(
i,1)**2+
p(
i,2)**2
157 IF(mstj(3).GE.1.AND.irank.EQ.1.AND.kflh.GE.4.AND.
158 &
p(
i,3).LE.0.001)
THEN
159 IF(
w.GE.
p(
i,5)+0.5*parj(32)) goto 180
174 IF(mstj(3).GE.0.AND.
p(
i,3).LT.0.)
i=
i-1
175 IF(
w.GT.parj(31)) goto 190
177 IF(mod(mstj(3),5).EQ.4.AND.
n.EQ.nsav1) wf=wf+0.1*parj(32)
178 IF(mod(mstj(3),5).EQ.4.AND.
n.EQ.nsav1) goto 170
181 the=
ulangl(
p(ip1,3),sqrt(
p(ip1,1)**2+
p(ip1,2)**2))
183 CALL ludbrb(nsav1+1,
n,the,
phi,0d0,0d0,0d0)
184 k(
k(ip1,3),4)=nsav1+1
189 IF(njet.EQ.1.OR.mstj(3).LE.0) goto 470
190 IF(mod(mstj(3),5).NE.0.AND.
n-nsav-njet.LT.2) goto 150
193 DO 240
i=nsav+njet+1,
n
195 kfla=mod(kfa/1000,10)
199 IF(kflb.LE.3) nfl(kflb)=nfl(kflb)-isign(1,
k(
i,2))*(-1)**kflb
200 IF(kflc.LE.3) nfl(kflc)=nfl(kflc)+isign(1,
k(
i,2))*(-1)**kflb
202 IF(kfla.LE.3) nfl(kfla)=nfl(kfla)-isign(1,
k(
i,2))
203 IF(kflb.LE.3) nfl(kflb)=nfl(kflb)-isign(1,
k(
i,2))
204 IF(kflc.LE.3) nfl(kflc)=nfl(kflc)-isign(1,
k(
i,2))
207 nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+
208 &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3
209 IF(nreq.EQ.0) goto 320
215 DO 260
i=nsav+njet+1,
n
216 p2=
p(
i,1)**2+
p(
i,2)**2+
p(
i,3)**2
217 IF(
k(
i,1).EQ.1.AND.p2.LT.p2min) irem=
i
218 260
IF(
k(
i,1).EQ.1.AND.p2.LT.p2min) p2min=p2
219 IF(irem.EQ.0) goto 150
222 kfla=mod(kfa/1000,10)
225 IF(kfla.GE.4.OR.kflb.GE.4)
k(irem,1)=8
226 IF(
k(irem,1).EQ.8) goto 250
228 isgn=isign(1,
k(irem,2))*(-1)**kflb
229 IF(kflb.LE.3) nfl(kflb)=nfl(kflb)+isgn
230 IF(kflc.LE.3) nfl(kflc)=nfl(kflc)-isgn
232 IF(kfla.LE.3) nfl(kfla)=nfl(kfla)+isign(1,
k(irem,2))
233 IF(kflb.LE.3) nfl(kflb)=nfl(kflb)+isign(1,
k(irem,2))
234 IF(kflc.LE.3) nfl(kflc)=nfl(kflc)+isign(1,
k(irem,2))
237 nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+
238 &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3
239 IF(nreq.GT.nrem) goto 250
240 DO 270
i=nsav+njet+1,
n
241 270
IF(
k(
i,1).EQ.8)
k(
i,1)=1
245 IF(nfl(1)+nfl(2)+nfl(3).NE.0) nfet=3
246 IF(nreq.LT.nrem) nfet=1
247 IF(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3)).EQ.0) nfet=0
249 ifet(
j)=1+(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3)))*
rlu(0)
250 kflf(
j)=isign(1,nfl(1))
251 IF(ifet(
j).GT.iabs(nfl(1))) kflf(
j)=isign(2,nfl(2))
252 290
IF(ifet(
j).GT.iabs(nfl(1))+iabs(nfl(2))) kflf(
j)=isign(3,nfl(3))
253 IF(nfet.EQ.2.AND.(ifet(1).EQ.ifet(2).OR.kflf(1)*kflf(2).GT.0))
255 IF(nfet.EQ.3.AND.(ifet(1).EQ.ifet(2).OR.ifet(1).EQ.ifet(3).OR.
256 &ifet(2).EQ.ifet(3).OR.kflf(1)*kflf(2).LT.0.OR.kflf(1)*kflf(3).
257 <.0.OR.kflf(1)*(nfl(1)+nfl(2)+nfl(3)).LT.0)) goto 280
258 IF(nfet.EQ.0) kflf(1)=1+int((2.+parj(2))*
rlu(0))
259 IF(nfet.EQ.0) kflf(2)=-kflf(1)
260 IF(nfet.EQ.1) kflf(2)=isign(1+int((2.+parj(2))*
rlu(0)),-kflf(1))
261 IF(nfet.LE.2) kflf(3)=0
262 IF(kflf(3).NE.0)
THEN
263 kflfc=isign(1000*
max(iabs(kflf(1)),iabs(kflf(3)))+
264 & 100*min(iabs(kflf(1)),iabs(kflf(3)))+1,kflf(1))
265 IF(kflf(1).EQ.kflf(3).OR.(1.+3.*parj(4))*
rlu(0).GT.1.)
266 & kflfc=kflfc+isign(2,kflfc)
270 CALL
lukfdi(kflfc,kflf(2),kfldmp,kf)
272 DO 300
j=1,
max(2,nfet)
273 300 nfl(iabs(kflf(
j)))=nfl(iabs(kflf(
j)))-isign(1,kflf(
j))
276 npos=min(1+int(
rlu(0)*nrem),nrem)
277 DO 310
i=nsav+njet+1,
n
278 IF(
k(
i,1).EQ.7) npos=npos-1
279 IF(
k(
i,1).EQ.1.OR.npos.NE.0) goto 310
283 p(
i,4)=sqrt(
p(
i,1)**2+
p(
i,2)**2+
p(
i,3)**2+
p(
i,5)**2)
286 nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+
287 &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3
288 IF(nrem.GT.0) goto 280
291 320
IF(mod(mstj(3),5).NE.0.AND.mod(mstj(3),5).NE.4)
THEN
294 DO 330
i=nsav+njet+1,
n
295 330 psi(
j)=psi(
j)+
p(
i,
j)
296 psi(4)=psi(1)**2+psi(2)**2+psi(3)**2
298 DO 340
i=nsav+njet+1,
n
299 IF(mod(mstj(3),5).EQ.1) pws=pws+
p(
i,4)
300 IF(mod(mstj(3),5).EQ.2) pws=pws+sqrt(
p(
i,5)**2+(psi(1)*
p(
i,1)+
301 & psi(2)*
p(
i,2)+psi(3)*
p(
i,3))**2/psi(4))
302 340
IF(mod(mstj(3),5).EQ.3) pws=pws+1.
303 DO 360
i=nsav+njet+1,
n
304 IF(mod(mstj(3),5).EQ.1) pw=
p(
i,4)
305 IF(mod(mstj(3),5).EQ.2) pw=sqrt(
p(
i,5)**2+(psi(1)*
p(
i,1)+
306 & psi(2)*
p(
i,2)+psi(3)*
p(
i,3))**2/psi(4))
307 IF(mod(mstj(3),5).EQ.3) pw=1.
309 350
p(
i,
j)=
p(
i,
j)-psi(
j)*pw/pws
310 360
p(
i,4)=sqrt(
p(
i,1)**2+
p(
i,2)**2+
p(
i,3)**2+
p(
i,5)**2)
313 ELSEIF(mod(mstj(3),5).EQ.4)
THEN
318 DO 390
i=nsav+njet+1,
n
322 pls=(
p(
i,1)*
p(ir1,1)+
p(
i,2)*
p(ir1,2)+
p(
i,3)*
p(ir1,3))/
323 & (
p(ir1,1)**2+
p(ir1,2)**2+
p(ir1,3)**2)
325 380
p(ir2,
j)=
p(ir2,
j)+
p(
i,
j)-pls*
p(ir1,
j)
326 p(ir2,4)=
p(ir2,4)+
p(
i,4)
327 390
p(ir2,5)=
p(ir2,5)+pls
330 400
IF(
k(
i,1).NE.0) pss=pss+
p(
i,4)/(pecm*(0.8*
p(
i,5)+0.2))
331 DO 420
i=nsav+njet+1,
n
334 pls=(
p(
i,1)*
p(ir1,1)+
p(
i,2)*
p(ir1,2)+
p(
i,3)*
p(ir1,3))/
335 & (
p(ir1,1)**2+
p(ir1,2)**2+
p(ir1,3)**2)
337 410
p(
i,
j)=
p(
i,
j)-
p(ir2,
j)/
k(ir2,1)+(1./(
p(ir2,5)*pss)-1.)*pls*
339 420
p(
i,4)=sqrt(
p(
i,1)**2+
p(
i,2)**2+
p(
i,3)**2+
p(
i,5)**2)
343 IF(mod(mstj(3),5).NE.0)
THEN
347 DO 430
i=nsav+njet+1,
n
350 430 pqs=pqs+
p(
i,5)**2/
p(
i,4)
351 IF(pms.GE.pecm) goto 150
354 pfac=(pecm-pqs)/(pes-pqs)
357 DO 460
i=nsav+njet+1,
n
360 p(
i,4)=sqrt(
p(
i,1)**2+
p(
i,2)**2+
p(
i,3)**2+
p(
i,5)**2)
362 460 pqs=pqs+
p(
i,5)**2/
p(
i,4)
363 IF(neco.LT.10.AND.abs(pecm-pes).GT.2
e-6*pecm) goto 440
367 470
DO 480
i=nsav+njet+1,
n
368 IF(mstu(16).NE.2)
k(
i,3)=nsav+1
369 480
IF(mstu(16).EQ.2)
k(
i,3)=
k(
k(
i,3),3)
370 DO 490
i=nsav+1,nsav+njet
373 IF(mstu(16).NE.2)
THEN
379 IF(
k(
i1,5).LT.
k(
i1,4))
THEN
395 500
v(nsav,
j)=
v(ip,
j)
396 p(nsav,5)=sqrt(
max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))
406 IF(njet.NE.1) CALL ludbrb(nsav+1,
n,0.,0.,dps(1)/dps(4),
407 &dps(2)/dps(4),dps(3)/dps(4))