3 C*********************************************************************
5 SUBROUTINE PYRESD_HIJING
7 C...Allows resonances to decay (including parton showers for hadronic
9 IMPLICIT DOUBLE PRECISION(D)
10 #include "lujets_hijing.inc"
11 #include "ludat1_hijing.inc"
12 #include "ludat2_hijing.inc"
13 #include "ludat3_hijing.inc"
14 #include "pysubs_hijing.inc"
15 #include "pypars_hijing.inc"
16 #include "pyint1_hijing.inc"
17 #include "pyint2_hijing.inc"
18 #include "pyint4_hijing.inc"
19 DIMENSION IREF(10,6),KDCY(2),KFL1(2),KFL2(2),NSD(2),ILIN(6),
20 &COUP(6,4),PK(6,4),PKK(6,6),CTHE(2),PHI(2),WDTP(0:40),
22 COMPLEX FGK,HA(6,6),HC(6,6)
24 C...The F, Xi and Xj functions of Gunion and Kunszt
25 C...(Phys. Rev. D33, 665, plus errata from the authors).
26 FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
27 &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
28 DIGK(DT,DU)=-4.*D34*D56+DT*(3.*DT+4.*DU)+DT**2*(DT*DU/(D34*D56)-
29 &2.*(1./D34+1./D56)*(DT+DU)+2.*(D34/D56+D56/D34))
30 DJGK(DT,DU)=8.*(D34+D56)**2-8.*(D34+D56)*(DT+DU)-6.*DT*DU-
31 &2.*DT*DU*(DT*DU/(D34*D56)-2.*(1./D34+1./D56)*(DT+DU)+
32 &2.*(D34/D56+D56/D34))
34 C...Define initial two objects, initialize loop.
39 IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
40 IREF(1,1)=MINT(84)+2+ISET(ISUB)
42 IREF(1,3)=MINT(83)+6+ISET(ISUB)
44 ELSEIF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
45 IREF(1,1)=MINT(84)+1+ISET(ISUB)
46 IREF(1,2)=MINT(84)+2+ISET(ISUB)
47 IREF(1,3)=MINT(83)+5+ISET(ISUB)
48 IREF(1,4)=MINT(83)+6+ISET(ISUB)
55 C...Loop over one/two resonances; reset decay rates.
57 IF(IP.EQ.1.AND.(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3)) JTMAX=1
66 IF(KFA.LT.23.OR.KFA.GT.40) GOTO 140
67 IF(MDCY(KFA,1).NE.0) THEN
68 IF(ISUB.EQ.1.OR.ISUB.EQ.141) MINT(61)=1
69 CALL PYWIDT_HIJING(KFA,P(ID,5),WDTP,WDTE)
70 IF(KCHG(KFA,3).EQ.0) THEN
73 IPM=(5+ISIGN(1,K(ID,2)))/2
75 IF(JTMAX.EQ.1.OR.IABS(K(IREF(IP,1),2)).NE.IABS(K(IREF(IP,2),2)))
79 IF(JT.EQ.1) I12=INT(4.5+RLU_HIJING(0))
82 RKFL=(WDTE(0,1)+WDTE(0,IPM)+WDTE(0,I12))*RLU_HIJING(0)
83 DO 120 I=1,MDCY(KFA,3)
85 KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
86 KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
87 RKFL=RKFL-(WDTE(I,1)+WDTE(I,IPM)+WDTE(I,I12))
88 IF(RKFL.LE.0.) GOTO 130
93 C...Summarize result on decay channel chosen.
94 IF((KFA.EQ.23.OR.KFA.EQ.24).AND.KFL1(JT).EQ.0) NINH=NINH+1
95 IF(KFL1(JT).EQ.0) GOTO 140
97 IF(IABS(KFL1(JT)).LE.10.OR.KFL1(JT).EQ.21) KDCY(JT)=1
98 IF((IABS(KFL1(JT)).GE.23.AND.IABS(KFL1(JT)).LE.25).OR.
99 &(IABS(KFL1(JT)).EQ.37)) KDCY(JT)=3
102 C...Fill decay products, prepared for parton showers for quarks.
103 IF(KDCY(JT).EQ.1) THEN
104 CALL LU2ENT_HIJING(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
106 CALL LU2ENT_HIJING(N+1,KFL1(JT),KFL2(JT),P(ID,5))
109 CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))
111 IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
114 CTHE(JT)=2.*RLU_HIJING(0)-1.
115 PHI(JT)=PARU(2)*RLU_HIJING(0)
118 IF(MINT(3).EQ.1.AND.IP.EQ.1) THEN
122 IF(JTMAX.EQ.1.AND.KDCY(1).EQ.0) GOTO 530
123 IF(JTMAX.EQ.2.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 530
124 IF(MSTP(45).LE.0.OR.IREF(IP,2).EQ.0.OR.NINH.GE.1) GOTO 500
125 IF(K(IREF(1,1),2).EQ.25.AND.IP.EQ.1) GOTO 500
126 IF(K(IREF(1,1),2).EQ.25.AND.KDCY(1)*KDCY(2).EQ.0) GOTO 500
128 C...Order incoming partons and outgoing resonances.
130 IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
131 IF(K(ILIN(1),2).EQ.21) ILIN(1)=2*MINT(84)+3-ILIN(1)
132 ILIN(2)=2*MINT(84)+3-ILIN(1)
134 IF(IREF(IP,5).EQ.25) IMIN=3
137 IF(K(IREF(IP,1),2).EQ.23) IORD=2
138 IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
139 IF(IABS(K(IREF(IP,IORD),2)).EQ.25) IORD=3-IORD
140 IF(KDCY(IORD).EQ.0) IORD=3-IORD
142 C...Order decay products of resonances.
143 DO 390 JT=IORD,3-IORD,3-2*IORD
144 IF(KDCY(JT).EQ.0) THEN
147 ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
148 ILIN(IMAX+1)=N+2*JT-1
151 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
152 K(N+2*JT,2)=K(NSD(JT)+2,2)
155 ILIN(IMAX+2)=N+2*JT-1
157 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
158 K(N+2*JT,2)=K(NSD(JT)+2,2)
162 C...Find charge, isospin, left- and righthanded couplings.
167 KFA=IABS(K(ILIN(I),2))
168 IF(KFA.GT.20) GOTO 410
169 COUP(I,1)=LUCHGE_HIJING(KFA)/3.
170 COUP(I,2)=(-1)**MOD(KFA,2)
171 COUP(I,4)=-2.*COUP(I,1)*XW
172 COUP(I,3)=COUP(I,2)+COUP(I,4)
175 GZMZ=PMAS(23,1)*PMAS(23,2)
177 GZMW=PMAS(24,1)*PMAS(24,2)
179 GZMZP=PMAS(32,1)*PMAS(32,2)
181 C...Select random angles; construct massless four-vectors.
187 IF(KDCY(JT).EQ.0) GOTO 440
189 P(N+2*JT-1,3)=0.5*P(ID,5)
190 P(N+2*JT-1,4)=0.5*P(ID,5)
191 P(N+2*JT,3)=-0.5*P(ID,5)
192 P(N+2*JT,4)=0.5*P(ID,5)
193 CTHE(JT)=2.*RLU_HIJING(0)-1.
194 PHI(JT)=PARU(2)*RLU_HIJING(0)
195 CALL LUDBRB_HIJING(N+2*JT-1,N+2*JT,ACOS(CTHE(JT)),PHI(JT),
196 &DBLE(P(ID,1)/P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
199 C...Store incoming and outgoing momenta, with random rotation to
200 C...avoid accidental zeroes in HA expressions.
203 P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
205 P(N+4+I,5)=P(ILIN(I),5)
207 450 P(N+4+I,J)=P(ILIN(I),J)
208 THERR=ACOS(2.*RLU_HIJING(0)-1.)
209 PHIRR=PARU(2)*RLU_HIJING(0)
210 CALL LUDBRB_HIJING(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
213 460 PK(I,J)=P(N+4+I,J)
215 C...Calculate internal products.
216 IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25) THEN
217 DO 470 I1=IMIN,IMAX-1
219 HA(I1,I2)=SQRT((PK(I1,4)-PK(I1,3))*(PK(I2,4)+PK(I2,3))/
220 & (1E-20+PK(I1,1)**2+PK(I1,2)**2))*CMPLX(PK(I1,1),PK(I1,2))-
221 & SQRT((PK(I1,4)+PK(I1,3))*(PK(I2,4)-PK(I2,3))/
222 & (1E-20+PK(I2,1)**2+PK(I2,2)**2))*CMPLX(PK(I2,1),PK(I2,2))
223 HC(I1,I2)=CONJG(HA(I1,I2))
224 IF(I1.LE.2) HA(I1,I2)=CMPLX(0.,1.)*HA(I1,I2)
225 IF(I1.LE.2) HC(I1,I2)=CMPLX(0.,1.)*HC(I1,I2)
227 470 HC(I2,I1)=-HC(I1,I2)
232 DO 490 I1=IMIN,IMAX-1
234 PKK(I1,I2)=2.*(PK(I1,4)*PK(I2,4)-PK(I1,1)*PK(I2,1)-
235 &PK(I1,2)*PK(I2,2)-PK(I1,3)*PK(I2,3))
236 490 PKK(I2,I1)=PKK(I1,I2)
238 IF(IREF(IP,5).EQ.25) THEN
239 C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons
240 WT=16.*PKK(3,5)*PKK(4,6)
241 IF(IP.EQ.1) WTMAX=SH**2
242 IF(IP.GE.2) WTMAX=P(IREF(IP,6),5)**4
244 ELSEIF(ISUB.EQ.1) THEN
246 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons
247 EI=KCHG(IABS(MINT(15)),1)/3.
254 GZ=1./(8.*XW*(1.-XW))*SH*(SH-SQMZ)/((SH-SQMZ)**2+GZMZ**2)
255 ZZ=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZ)**2+GZMZ**2)
256 IF(MSTP(43).EQ.1) THEN
257 C...Only gamma* production included
260 ELSEIF(MSTP(43).EQ.2) THEN
261 C...Only Z0 production included
265 ASYM=2.*(EI*AI*GZ*EF*AF+4.*VI*AI*ZZ*VF*AF)/(EI**2*GG*EF**2+
266 & EI*VI*GZ*EF*VF+(VI**2+AI**2)*ZZ*(VF**2+AF**2))
267 WT=1.+ASYM*CTHE(JT)+CTHE(JT)**2
270 C...Angular weight for gamma*/Z0 -> H+ + H-
275 ELSEIF(ISUB.EQ.2) THEN
276 C...Angular weight for W+/- -> 2 quarks/leptons
280 ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
281 C...Angular weight for f + fb -> gluon/gamma + Z0 ->
282 C...-> gluon/gamma + 2 quarks/leptons
283 WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
284 & (PKK(1,3)**2+PKK(2,4)**2)+((COUP(1,3)*COUP(3,4))**2+
285 & (COUP(1,4)*COUP(3,3))**2)*(PKK(1,4)**2+PKK(2,3)**2)
286 WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
287 & ((PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2)
289 ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
290 C...Angular weight for f + fb' -> gluon/gamma + W+/- ->
291 C...-> gluon/gamma + 2 quarks/leptons
292 WT=PKK(1,3)**2+PKK(2,4)**2
293 WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2
295 ELSEIF(ISUB.EQ.22) THEN
296 C...Angular weight for f + fb -> Z0 + Z0 -> 4 quarks/leptons
297 S34=P(IREF(IP,IORD),5)**2
298 S56=P(IREF(IP,3-IORD),5)**2
299 TI=PKK(1,3)+PKK(1,4)+S34
300 UI=PKK(1,5)+PKK(1,6)+S56
301 WT=COUP(1,3)**4*((COUP(3,3)*COUP(5,3)*ABS(FGK(1,2,3,4,5,6)/
302 & TI+FGK(1,2,5,6,3,4)/UI))**2+(COUP(3,4)*COUP(5,3)*ABS(
303 & FGK(1,2,4,3,5,6)/TI+FGK(1,2,5,6,4,3)/UI))**2+(COUP(3,3)*
304 & COUP(5,4)*ABS(FGK(1,2,3,4,6,5)/TI+FGK(1,2,6,5,3,4)/UI))**2+
305 & (COUP(3,4)*COUP(5,4)*ABS(FGK(1,2,4,3,6,5)/TI+FGK(1,2,6,5,4,3)/
306 & UI))**2)+COUP(1,4)**4*((COUP(3,3)*COUP(5,3)*ABS(
307 & FGK(2,1,5,6,3,4)/TI+FGK(2,1,3,4,5,6)/UI))**2+(COUP(3,4)*
308 & COUP(5,3)*ABS(FGK(2,1,6,5,3,4)/TI+FGK(2,1,3,4,6,5)/UI))**2+
309 & (COUP(3,3)*COUP(5,4)*ABS(FGK(2,1,5,6,4,3)/TI+FGK(2,1,4,3,5,6)/
310 & UI))**2+(COUP(3,4)*COUP(5,4)*ABS(FGK(2,1,6,5,4,3)/TI+
311 & FGK(2,1,4,3,6,5)/UI))**2)
312 WTMAX=4.*S34*S56*(COUP(1,3)**4+COUP(1,4)**4)*(COUP(3,3)**2+
313 & COUP(3,4)**2)*(COUP(5,3)**2+COUP(5,4)**2)*4.*(TI/UI+UI/TI+
314 & 2.*SH*(S34+S56)/(TI*UI)-S34*S56*(1./TI**2+1./UI**2))
316 ELSEIF(ISUB.EQ.23) THEN
317 C...Angular weight for f + fb' -> Z0 + W +/- -> 4 quarks/leptons
318 D34=P(IREF(IP,IORD),5)**2
319 D56=P(IREF(IP,3-IORD),5)**2
320 DT=PKK(1,3)+PKK(1,4)+D34
321 DU=PKK(1,5)+PKK(1,6)+D56
322 CAWZ=COUP(2,3)/SNGL(DT)-2.*(1.-XW)*COUP(1,2)/(SH-SQMW)
323 CBWZ=COUP(1,3)/SNGL(DU)+2.*(1.-XW)*COUP(1,2)/(SH-SQMW)
324 WT=COUP(5,3)**2*ABS(CAWZ*FGK(1,2,3,4,5,6)+CBWZ*
325 & FGK(1,2,5,6,3,4))**2+COUP(5,4)**2*ABS(CAWZ*
326 & FGK(1,2,3,4,6,5)+CBWZ*FGK(1,2,6,5,3,4))**2
327 WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*(CAWZ**2*
328 & DIGK(DT,DU)+CBWZ**2*DIGK(DU,DT)+CAWZ*CBWZ*DJGK(DT,DU))
330 ELSEIF(ISUB.EQ.24) THEN
331 C...Angular weight for f + fb -> Z0 + H0 -> 2 quarks/leptons + H0
332 WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
333 & PKK(1,3)*PKK(2,4)+((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*
334 & COUP(3,3))**2)*PKK(1,4)*PKK(2,3)
335 WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
336 & (PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
338 ELSEIF(ISUB.EQ.25) THEN
339 C...Angular weight for f + fb -> W+ + W- -> 4 quarks/leptons
340 D34=P(IREF(IP,IORD),5)**2
341 D56=P(IREF(IP,3-IORD),5)**2
342 DT=PKK(1,3)+PKK(1,4)+D34
343 DU=PKK(1,5)+PKK(1,6)+D56
344 CDWW=(COUP(1,3)*SQMZ/(SH-SQMZ)+COUP(1,2))/SH
345 CAWW=CDWW+0.5*(COUP(1,2)+1.)/SNGL(DT)
346 CBWW=CDWW+0.5*(COUP(1,2)-1.)/SNGL(DU)
347 CCWW=COUP(1,4)*SQMZ/(SH-SQMZ)/SH
348 WT=ABS(CAWW*FGK(1,2,3,4,5,6)-CBWW*FGK(1,2,5,6,3,4))**2+
349 & CCWW**2*ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))**2
350 WTMAX=4.*D34*D56*(CAWW**2*DIGK(DT,DU)+CBWW**2*DIGK(DU,DT)-CAWW*
351 & CBWW*DJGK(DT,DU)+CCWW**2*(DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU)))
353 ELSEIF(ISUB.EQ.26) THEN
354 C...Angular weight for f + fb' -> W+/- + H0 -> 2 quarks/leptons + H0
356 WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
358 ELSEIF(ISUB.EQ.30) THEN
359 C...Angular weight for f + g -> f + Z0 -> f + 2 quarks/leptons
360 IF(K(ILIN(1),2).GT.0) WT=((COUP(1,3)*COUP(3,3))**2+
361 & (COUP(1,4)*COUP(3,4))**2)*(PKK(1,4)**2+PKK(3,5)**2)+
362 & ((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*COUP(3,3))**2)*
363 & (PKK(1,3)**2+PKK(4,5)**2)
364 IF(K(ILIN(1),2).LT.0) WT=((COUP(1,3)*COUP(3,3))**2+
365 & (COUP(1,4)*COUP(3,4))**2)*(PKK(1,3)**2+PKK(4,5)**2)+
366 & ((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*COUP(3,3))**2)*
367 & (PKK(1,4)**2+PKK(3,5)**2)
368 WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
369 & ((PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2)
371 ELSEIF(ISUB.EQ.31) THEN
372 C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons
373 IF(K(ILIN(1),2).GT.0) WT=PKK(1,4)**2+PKK(3,5)**2
374 IF(K(ILIN(1),2).LT.0) WT=PKK(1,3)**2+PKK(4,5)**2
375 WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2
377 ELSEIF(ISUB.EQ.141) THEN
378 C...Angular weight for gamma*/Z0/Z'0 -> 2 quarks/leptons
379 EI=KCHG(IABS(MINT(15)),1)/3.
390 GZ=1./(8.*XW*(1.-XW))*SH*(SH-SQMZ)/((SH-SQMZ)**2+GZMZ**2)
391 GZP=1./(8.*XW*(1.-XW))*SH*(SH-SQMZP)/((SH-SQMZP)**2+GZMZP**2)
392 ZZ=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZ)**2+GZMZ**2)
393 ZZP=2./(16.*XW*(1.-XW))**2*
394 & SH**2*((SH-SQMZ)*(SH-SQMZP)+GZMZ*GZMZP)/
395 & (((SH-SQMZ)**2+GZMZ**2)*((SH-SQMZP)**2+GZMZP**2))
396 ZPZP=1./(16.*XW*(1.-XW))**2*SH**2/((SH-SQMZP)**2+GZMZP**2)
397 IF(MSTP(44).EQ.1) THEN
398 C...Only gamma* production included
404 ELSEIF(MSTP(44).EQ.2) THEN
405 C...Only Z0 production included
411 ELSEIF(MSTP(44).EQ.3) THEN
412 C...Only Z'0 production included
418 ELSEIF(MSTP(44).EQ.4) THEN
419 C...Only gamma*/Z0 production included
423 ELSEIF(MSTP(44).EQ.5) THEN
424 C...Only gamma*/Z'0 production included
428 ELSEIF(MSTP(44).EQ.6) THEN
429 C...Only Z0/Z'0 production included
434 ASYM=2.*(EI*AI*GZ*EF*AF+EI*API*GZP*EF*APF+4.*VI*AI*ZZ*VF*AF+
435 & (VI*API+VPI*AI)*ZZP*(VF*APF+VPF*AF)+4.*VPI*API*ZPZP*VPF*APF)/
436 & (EI**2*GG*EF**2+EI*VI*GZ*EF*VF+EI*VPI*GZP*EF*VPF+
437 & (VI**2+AI**2)*ZZ*(VF**2+AF**2)+(VI*VPI+AI*API)*ZZP*
438 & (VF*VPF+AF*APF)+(VPI**2+API**2)*ZPZP*(VPF**2+APF**2))
439 WT=1.+ASYM*CTHE(JT)+CTHE(JT)**2
446 C...Obtain correct angular distribution by rejection techniques.
447 IF(WT.LT.RLU_HIJING(0)*WTMAX) GOTO 420
449 C...Construct massive four-vectors using angles chosen. Mark decayed
450 C...resonances, add documentation lines. Shower evolution.
451 500 DO 520 JT=1,JTMAX
452 IF(KDCY(JT).EQ.0) GOTO 520
454 CALL LUDBRB_HIJING(NSD(JT)+1,NSD(JT)+2,ACOS(CTHE(JT)),PHI(JT),
455 &DBLE(P(ID,1)/P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
459 IDOC=MINT(83)+MINT(4)
460 DO 510 I=NSD(JT)+1,NSD(JT)+2
466 K(I1,3)=IREF(IP,JT+2)
470 MINT(7)=MINT(83)+6+2*ISET(ISUB)
471 MINT(8)=MINT(83)+7+2*ISET(ISUB)
473 IF(MSTP(71).GE.1.AND.KDCY(JT).EQ.1) CALL LUSHOW_HIJING(NSD(JT)+1,
476 C...Check if new resonances were produced, loop back if needed.
477 IF(KDCY(JT).NE.3) GOTO 520
483 IREF(NP,5)=K(IREF(IP,JT),2)
484 IREF(NP,6)=IREF(IP,JT)
486 530 IF(IP.LT.NP) GOTO 100