]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/pyresd_hijing.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master_patch'
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pyresd_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE PYRESD_HIJING
6
7C...Allows resonances to decay (including parton showers for hadronic
8C...channels).
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),
21 &WDTE(0:40,0:5)
22 COMPLEX FGK,HA(6,6),HC(6,6)
23
24C...The F, Xi and Xj functions of Gunion and Kunszt
25C...(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))
33
34C...Define initial two objects, initialize loop.
35 ISUB=MINT(1)
36 SH=VINT(44)
37 IREF(1,5)=0
38 IREF(1,6)=0
39 IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
40 IREF(1,1)=MINT(84)+2+ISET(ISUB)
41 IREF(1,2)=0
42 IREF(1,3)=MINT(83)+6+ISET(ISUB)
43 IREF(1,4)=0
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)
49 ENDIF
50 NP=1
51 IP=0
52 100 IP=IP+1
53 NINH=0
54
55C...Loop over one/two resonances; reset decay rates.
56 JTMAX=2
57 IF(IP.EQ.1.AND.(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3)) JTMAX=1
58 DO 140 JT=1,JTMAX
59 KDCY(JT)=0
60 KFL1(JT)=0
61 KFL2(JT)=0
62 NSD(JT)=IREF(IP,JT)
63 ID=IREF(IP,JT)
64 IF(ID.EQ.0) GOTO 140
65 KFA=IABS(K(ID,2))
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
71 IPM=2
72 ELSE
73 IPM=(5+ISIGN(1,K(ID,2)))/2
74 ENDIF
75 IF(JTMAX.EQ.1.OR.IABS(K(IREF(IP,1),2)).NE.IABS(K(IREF(IP,2),2)))
76 & THEN
77 I12=4
78 ELSE
79 IF(JT.EQ.1) I12=INT(4.5+RLU_HIJING(0))
80 I12=9-I12
81 ENDIF
82 RKFL=(WDTE(0,1)+WDTE(0,IPM)+WDTE(0,I12))*RLU_HIJING(0)
83 DO 120 I=1,MDCY(KFA,3)
84 IDC=I+MDCY(KFA,2)-1
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
89 120 CONTINUE
90 130 CONTINUE
91 ENDIF
92
93C...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
96 KDCY(JT)=2
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
100 NSD(JT)=N
101
102C...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))
105 ELSE
106 CALL LU2ENT_HIJING(N+1,KFL1(JT),KFL2(JT),P(ID,5))
107 ENDIF
108 IF(JTMAX.EQ.1) THEN
109 CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))
110 $ *RLU_HIJING(0)
111 IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
112 PHI(JT)=VINT(24)
113 ELSE
114 CTHE(JT)=2.*RLU_HIJING(0)-1.
115 PHI(JT)=PARU(2)*RLU_HIJING(0)
116 ENDIF
117 140 CONTINUE
118 IF(MINT(3).EQ.1.AND.IP.EQ.1) THEN
119 MINT(25)=KFL1(1)
120 MINT(26)=KFL2(1)
121 ENDIF
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
127
128C...Order incoming partons and outgoing resonances.
129 ILIN(1)=MINT(84)+1
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)
133 IMIN=1
134 IF(IREF(IP,5).EQ.25) IMIN=3
135 IMAX=2
136 IORD=1
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
141
142C...Order decay products of resonances.
143 DO 390 JT=IORD,3-IORD,3-2*IORD
144 IF(KDCY(JT).EQ.0) THEN
145 ILIN(IMAX+1)=NSD(JT)
146 IMAX=IMAX+1
147 ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
148 ILIN(IMAX+1)=N+2*JT-1
149 ILIN(IMAX+2)=N+2*JT
150 IMAX=IMAX+2
151 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
152 K(N+2*JT,2)=K(NSD(JT)+2,2)
153 ELSE
154 ILIN(IMAX+1)=N+2*JT
155 ILIN(IMAX+2)=N+2*JT-1
156 IMAX=IMAX+2
157 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
158 K(N+2*JT,2)=K(NSD(JT)+2,2)
159 ENDIF
160 390 CONTINUE
161
162C...Find charge, isospin, left- and righthanded couplings.
163 XW=PARU(102)
164 DO 410 I=IMIN,IMAX
165 DO 400 J=1,4
166 400 COUP(I,J)=0.
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)
173 410 CONTINUE
174 SQMZ=PMAS(23,1)**2
175 GZMZ=PMAS(23,1)*PMAS(23,2)
176 SQMW=PMAS(24,1)**2
177 GZMW=PMAS(24,1)*PMAS(24,2)
178 SQMZP=PMAS(32,1)**2
179 GZMZP=PMAS(32,1)*PMAS(32,2)
180
181C...Select random angles; construct massless four-vectors.
182 420 DO 430 I=N+1,N+4
183 K(I,1)=1
184 DO 430 J=1,5
185 430 P(I,J)=0.
186 DO 440 JT=1,JTMAX
187 IF(KDCY(JT).EQ.0) GOTO 440
188 ID=IREF(IP,JT)
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)))
197 440 CONTINUE
198
199C...Store incoming and outgoing momenta, with random rotation to
200C...avoid accidental zeroes in HA expressions.
201 DO 450 I=1,IMAX
202 K(N+4+I,1)=1
203 P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
204 &P(ILIN(I),5)**2)
205 P(N+4+I,5)=P(ILIN(I),5)
206 DO 450 J=1,3
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)
211 DO 460 I=1,IMAX
212 DO 460 J=1,4
213 460 PK(I,J)=P(N+4+I,J)
214
215C...Calculate internal products.
216 IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25) THEN
217 DO 470 I1=IMIN,IMAX-1
218 DO 470 I2=I1+1,IMAX
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)
226 HA(I2,I1)=-HA(I1,I2)
227 470 HC(I2,I1)=-HC(I1,I2)
228 ENDIF
229 DO 480 I=1,2
230 DO 480 J=1,4
231 480 PK(I,J)=-PK(I,J)
232 DO 490 I1=IMIN,IMAX-1
233 DO 490 I2=I1+1,IMAX
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)
237
238 IF(IREF(IP,5).EQ.25) THEN
239C...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
243
244 ELSEIF(ISUB.EQ.1) THEN
245 IF(KFA.NE.37) THEN
246C...Angular weight for gamma*/Z0 -> 2 quarks/leptons
247 EI=KCHG(IABS(MINT(15)),1)/3.
248 AI=SIGN(1.,EI+0.1)
249 VI=AI-4.*EI*XW
250 EF=KCHG(KFA,1)/3.
251 AF=SIGN(1.,EF+0.1)
252 VF=AF-4.*EF*XW
253 GG=1.
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
257C...Only gamma* production included
258 GZ=0.
259 ZZ=0.
260 ELSEIF(MSTP(43).EQ.2) THEN
261C...Only Z0 production included
262 GG=0.
263 GZ=0.
264 ENDIF
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
268 WTMAX=2.+ABS(ASYM)
269 ELSE
270C...Angular weight for gamma*/Z0 -> H+ + H-
271 WT=1.-CTHE(JT)**2
272 WTMAX=1.
273 ENDIF
274
275 ELSEIF(ISUB.EQ.2) THEN
276C...Angular weight for W+/- -> 2 quarks/leptons
277 WT=(1.+CTHE(JT))**2
278 WTMAX=4.
279
280 ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
281C...Angular weight for f + fb -> gluon/gamma + Z0 ->
282C...-> 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)
288
289 ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
290C...Angular weight for f + fb' -> gluon/gamma + W+/- ->
291C...-> 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
294
295 ELSEIF(ISUB.EQ.22) THEN
296C...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))
315
316 ELSEIF(ISUB.EQ.23) THEN
317C...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))
329
330 ELSEIF(ISUB.EQ.24) THEN
331C...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))
337
338 ELSEIF(ISUB.EQ.25) THEN
339C...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)))
352
353 ELSEIF(ISUB.EQ.26) THEN
354C...Angular weight for f + fb' -> W+/- + H0 -> 2 quarks/leptons + H0
355 WT=PKK(1,3)*PKK(2,4)
356 WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
357
358 ELSEIF(ISUB.EQ.30) THEN
359C...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)
370
371 ELSEIF(ISUB.EQ.31) THEN
372C...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
376
377 ELSEIF(ISUB.EQ.141) THEN
378C...Angular weight for gamma*/Z0/Z'0 -> 2 quarks/leptons
379 EI=KCHG(IABS(MINT(15)),1)/3.
380 AI=SIGN(1.,EI+0.1)
381 VI=AI-4.*EI*XW
382 API=SIGN(1.,EI+0.1)
383 VPI=API-4.*EI*XW
384 EF=KCHG(KFA,1)/3.
385 AF=SIGN(1.,EF+0.1)
386 VF=AF-4.*EF*XW
387 APF=SIGN(1.,EF+0.1)
388 VPF=APF-4.*EF*XW
389 GG=1.
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
398C...Only gamma* production included
399 GZ=0.
400 GZP=0.
401 ZZ=0.
402 ZZP=0.
403 ZPZP=0.
404 ELSEIF(MSTP(44).EQ.2) THEN
405C...Only Z0 production included
406 GG=0.
407 GZ=0.
408 GZP=0.
409 ZZP=0.
410 ZPZP=0.
411 ELSEIF(MSTP(44).EQ.3) THEN
412C...Only Z'0 production included
413 GG=0.
414 GZ=0.
415 GZP=0.
416 ZZ=0.
417 ZZP=0.
418 ELSEIF(MSTP(44).EQ.4) THEN
419C...Only gamma*/Z0 production included
420 GZP=0.
421 ZZP=0.
422 ZPZP=0.
423 ELSEIF(MSTP(44).EQ.5) THEN
424C...Only gamma*/Z'0 production included
425 GZ=0.
426 ZZ=0.
427 ZZP=0.
428 ELSEIF(MSTP(44).EQ.6) THEN
429C...Only Z0/Z'0 production included
430 GG=0.
431 GZ=0.
432 GZP=0.
433 ENDIF
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
440 WTMAX=2.+ABS(ASYM)
441
442 ELSE
443 WT=1.
444 WTMAX=1.
445 ENDIF
446C...Obtain correct angular distribution by rejection techniques.
447 IF(WT.LT.RLU_HIJING(0)*WTMAX) GOTO 420
448
449C...Construct massive four-vectors using angles chosen. Mark decayed
450C...resonances, add documentation lines. Shower evolution.
451 500 DO 520 JT=1,JTMAX
452 IF(KDCY(JT).EQ.0) GOTO 520
453 ID=IREF(IP,JT)
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)))
456 K(ID,1)=K(ID,1)+10
457 K(ID,4)=NSD(JT)+1
458 K(ID,5)=NSD(JT)+2
459 IDOC=MINT(83)+MINT(4)
460 DO 510 I=NSD(JT)+1,NSD(JT)+2
461 MINT(4)=MINT(4)+1
462 I1=MINT(83)+MINT(4)
463 K(I,3)=I1
464 K(I1,1)=21
465 K(I1,2)=K(I,2)
466 K(I1,3)=IREF(IP,JT+2)
467 DO 510 J=1,5
468 510 P(I1,J)=P(I,J)
469 IF(JTMAX.EQ.1) THEN
470 MINT(7)=MINT(83)+6+2*ISET(ISUB)
471 MINT(8)=MINT(83)+7+2*ISET(ISUB)
472 ENDIF
473 IF(MSTP(71).GE.1.AND.KDCY(JT).EQ.1) CALL LUSHOW_HIJING(NSD(JT)+1,
474 &NSD(JT)+2,P(ID,5))
475
476C...Check if new resonances were produced, loop back if needed.
477 IF(KDCY(JT).NE.3) GOTO 520
478 NP=NP+1
479 IREF(NP,1)=NSD(JT)+1
480 IREF(NP,2)=NSD(JT)+2
481 IREF(NP,3)=IDOC+1
482 IREF(NP,4)=IDOC+2
483 IREF(NP,5)=K(IREF(IP,JT),2)
484 IREF(NP,6)=IREF(IP,JT)
485 520 CONTINUE
486 530 IF(IP.LT.NP) GOTO 100
487
488 RETURN
489 END