Protection in calc. of ZDR1, ZDR2 as in Pythia 6.3
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / lushow_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE LUSHOW_HIJING(IP1,IP2,QMAX)   
6     
7 C...Purpose: to generate timelike parton showers from given partons.    
8       IMPLICIT DOUBLE PRECISION(D)  
9 #include "lujets_hijing.inc"
10 #include "ludat1_hijing.inc"
11 #include "ludat2_hijing.inc"
12       DIMENSION PMTH(5,40),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4),  
13      &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4)   
14     
15 C...Initialization of cutoff masses etc.    
16       IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR.  
17      &QMAX.LE.MIN(PARJ(82),PARJ(83)).OR.MSTJ(41).GE.3) RETURN   
18       PMTH(1,21)=ULMASS_HIJING(21) 
19       PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25*PARJ(82)**2)   
20       PMTH(3,21)=2.*PMTH(2,21)  
21       PMTH(4,21)=PMTH(3,21) 
22       PMTH(5,21)=PMTH(3,21) 
23       PMTH(1,22)=ULMASS_HIJING(22) 
24       PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25*PARJ(83)**2)   
25       PMTH(3,22)=2.*PMTH(2,22)  
26       PMTH(4,22)=PMTH(3,22) 
27       PMTH(5,22)=PMTH(3,22) 
28       PMQTH1=PARJ(82)   
29       IF(MSTJ(41).EQ.2) PMQTH1=MIN(PARJ(82),PARJ(83))   
30       PMQTH2=PMTH(2,21) 
31       IF(MSTJ(41).EQ.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))   
32       DO 100 IF=1,8 
33       PMTH(1,IF)=ULMASS_HIJING(IF) 
34       PMTH(2,IF)=SQRT(PMTH(1,IF)**2+0.25*PMQTH1**2) 
35       PMTH(3,IF)=PMTH(2,IF)+PMQTH2  
36       PMTH(4,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(82)**2)+PMTH(2,21)    
37   100 PMTH(5,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(83)**2)+PMTH(2,22)    
38       PT2MIN=MAX(0.5*PARJ(82),1.1*PARJ(81))**2  
39       ALAMS=PARJ(81)**2 
40       ALFM=LOG(PT2MIN/ALAMS)    
41     
42 C...Store positions of shower initiating partons.   
43       M3JC=0    
44       IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN 
45         NPA=1   
46         IPA(1)=IP1  
47       ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-   
48      &MSTU(32))) THEN   
49         NPA=2   
50         IPA(1)=IP1  
51         IPA(2)=IP2  
52       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0.  
53      &AND.IP2.GE.-3) THEN   
54         NPA=IABS(IP2)   
55         DO 110 I=1,NPA  
56   110   IPA(I)=IP1+I-1  
57       ELSE  
58         CALL LUERRM_HIJING(12, 
59      &  '(LUSHOW_HIJING:) failed to reconstruct showering system') 
60         IF(MSTU(21).GE.1) RETURN    
61       ENDIF 
62     
63 C...Check on phase space available for emission.    
64       IREJ=0    
65       DO 120 J=1,5  
66   120 PS(J)=0.  
67       PM=0. 
68       DO 130 I=1,NPA    
69       KFLA(I)=IABS(K(IPA(I),2)) 
70       PMA(I)=P(IPA(I),5)    
71       IF(KFLA(I).NE.0.AND.(KFLA(I).LE.8.OR.KFLA(I).EQ.21))  
72      &PMA(I)=PMTH(3,KFLA(I))    
73       PM=PM+PMA(I)  
74       IF(KFLA(I).EQ.0.OR.(KFLA(I).GT.8.AND.KFLA(I).NE.21).OR.   
75      &PMA(I).GT.QMAX) IREJ=IREJ+1   
76       DO 130 J=1,4  
77   130 PS(J)=PS(J)+P(IPA(I),J)   
78       IF(IREJ.EQ.NPA) RETURN    
79       PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))   
80       IF(NPA.EQ.1) PS(5)=PS(4)  
81       IF(PS(5).LE.PM+PMQTH1) RETURN 
82       IF(NPA.EQ.2.AND.MSTJ(47).GE.1) THEN   
83         IF(KFLA(1).GE.1.AND.KFLA(1).LE.8.AND.KFLA(2).GE.1.AND.  
84      &  KFLA(2).LE.8) M3JC=1    
85         IF(MSTJ(47).GE.2) M3JC=1    
86       ENDIF 
87     
88 C...Define imagined single initiator of shower for parton system.   
89       NS=N  
90       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN  
91          CALL LUERRM_HIJING(11
92      $        ,'(LUSHOW_HIJING:) no more memory left in LUJETS_HIJING')   
93         IF(MSTU(21).GE.1) RETURN    
94       ENDIF 
95       IF(NPA.GE.2) THEN 
96         K(N+1,1)=11 
97         K(N+1,2)=21 
98         K(N+1,3)=0  
99         K(N+1,4)=0  
100         K(N+1,5)=0  
101         P(N+1,1)=0. 
102         P(N+1,2)=0. 
103         P(N+1,3)=0. 
104         P(N+1,4)=PS(5)  
105         P(N+1,5)=PS(5)  
106         V(N+1,5)=PS(5)**2   
107         N=N+1   
108       ENDIF 
109     
110 C...Loop over partons that may branch.  
111       NEP=NPA   
112       IM=NS 
113       IF(NPA.EQ.1) IM=NS-1  
114   140 IM=IM+1   
115       IF(N.GT.NS) THEN  
116         IF(IM.GT.N) GOTO 380    
117         KFLM=IABS(K(IM,2))  
118         IF(KFLM.EQ.0.OR.(KFLM.GT.8.AND.KFLM.NE.21)) GOTO 140    
119         IF(P(IM,5).LT.PMTH(2,KFLM)) GOTO 140    
120         IGM=K(IM,3) 
121       ELSE  
122         IGM=-1  
123       ENDIF 
124       IF(N+NEP.GT.MSTU(4)-MSTU(32)-5) THEN  
125          CALL LUERRM_HIJING(11
126      $        ,'(LUSHOW_HIJING:) no more memory left in LUJETS_HIJING')   
127         IF(MSTU(21).GE.1) RETURN    
128       ENDIF 
129     
130 C...Position of aunt (sister to branching parton).  
131 C...Origin and flavour of daughters.    
132       IAU=0 
133       IF(IGM.GT.0) THEN 
134         IF(K(IM-1,3).EQ.IGM) IAU=IM-1   
135         IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1 
136       ENDIF 
137       IF(IGM.GE.0) THEN 
138         K(IM,4)=N+1 
139         DO 150 I=1,NEP  
140   150   K(N+I,3)=IM 
141       ELSE  
142         K(N+1,3)=IPA(1) 
143       ENDIF 
144       IF(IGM.LE.0) THEN 
145         DO 160 I=1,NEP  
146   160   K(N+I,2)=K(IPA(I),2)    
147       ELSEIF(KFLM.NE.21) THEN   
148         K(N+1,2)=K(IM,2)    
149         K(N+2,2)=K(IM,5)    
150       ELSEIF(K(IM,5).EQ.21) THEN    
151         K(N+1,2)=21 
152         K(N+2,2)=21 
153       ELSE  
154         K(N+1,2)=K(IM,5)    
155         K(N+2,2)=-K(IM,5)   
156       ENDIF 
157     
158 C...Reset flags on daughers and tries made. 
159       DO 170 IP=1,NEP   
160       K(N+IP,1)=3   
161       K(N+IP,4)=0   
162       K(N+IP,5)=0   
163       KFLD(IP)=IABS(K(N+IP,2))  
164       ITRY(IP)=0    
165       ISL(IP)=0 
166       ISI(IP)=0 
167   170 IF(KFLD(IP).GT.0.AND.(KFLD(IP).LE.8.OR.KFLD(IP).EQ.21)) ISI(IP)=1 
168       ISLM=0    
169     
170 C...Maximum virtuality of daughters.    
171       IF(IGM.LE.0) THEN 
172         DO 180 I=1,NPA  
173         IF(NPA.GE.3) P(N+I,4)=(PS(4)*P(IPA(I),4)-PS(1)*P(IPA(I),1)- 
174      &  PS(2)*P(IPA(I),2)-PS(3)*P(IPA(I),3))/PS(5)  
175         P(N+I,5)=MIN(QMAX,PS(5))    
176         IF(NPA.GE.3) P(N+I,5)=MIN(P(N+I,5),P(N+I,4))    
177   180   IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)    
178       ELSE  
179         IF(MSTJ(43).LE.2) PEM=V(IM,2)   
180         IF(MSTJ(43).GE.3) PEM=P(IM,4)   
181         P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)   
182         P(N+2,5)=MIN(P(IM,5),(1.-V(IM,1))*PEM)  
183         IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)  
184       ENDIF 
185       DO 190 I=1,NEP    
186       PMSD(I)=P(N+I,5)  
187       IF(ISI(I).EQ.1) THEN  
188         IF(P(N+I,5).LE.PMTH(3,KFLD(I))) P(N+I,5)=PMTH(1,KFLD(I))    
189       ENDIF 
190   190 V(N+I,5)=P(N+I,5)**2  
191     
192 C...Choose one of the daughters for evolution.  
193   200 INUM=0    
194       IF(NEP.EQ.1) INUM=1   
195       DO 210 I=1,NEP    
196   210 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I  
197       DO 220 I=1,NEP    
198       IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN   
199         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) INUM=I  
200       ENDIF 
201   220 CONTINUE  
202       IF(INUM.EQ.0) THEN    
203         RMAX=0. 
204         DO 230 I=1,NEP  
205         IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQTH2) THEN  
206           RPM=P(N+I,5)/PMSD(I)  
207           IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,KFLD(I))) THEN  
208             RMAX=RPM    
209             INUM=I  
210           ENDIF 
211         ENDIF   
212   230   CONTINUE    
213       ENDIF 
214     
215 C...Store information on choice of evolving daughter.   
216       INUM=MAX(1,INUM)  
217       IEP(1)=N+INUM 
218       DO 240 I=2,NEP    
219       IEP(I)=IEP(I-1)+1 
220   240 IF(IEP(I).GT.N+NEP) IEP(I)=N+1    
221       DO 250 I=1,NEP    
222   250 KFL(I)=IABS(K(IEP(I),2))  
223       ITRY(INUM)=ITRY(INUM)+1   
224       IF(ITRY(INUM).GT.200) THEN    
225          CALL LUERRM_HIJING(14
226      $        ,'(LUSHOW_HIJING:) caught in infinite loop') 
227         IF(MSTU(21).GE.1) RETURN    
228       ENDIF 
229       Z=0.5 
230       IF(KFL(1).EQ.0.OR.(KFL(1).GT.8.AND.KFL(1).NE.21)) GOTO 300    
231       IF(P(IEP(1),5).LT.PMTH(2,KFL(1))) GOTO 300    
232     
233 C...Calculate allowed z range.  
234       IF(NEP.EQ.1) THEN 
235         PMED=PS(4)  
236       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN    
237         PMED=P(IM,5)    
238       ELSE  
239         IF(INUM.EQ.1) PMED=V(IM,1)*PEM  
240         IF(INUM.EQ.2) PMED=(1.-V(IM,1))*PEM 
241       ENDIF 
242       IF(MOD(MSTJ(43),2).EQ.1) THEN 
243         ZC=PMTH(2,21)/PMED  
244         ZCE=PMTH(2,22)/PMED 
245       ELSE  
246         ZC=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,21)/PMED)**2)))    
247         IF(ZC.LT.1E-4) ZC=(PMTH(2,21)/PMED)**2  
248         ZCE=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,22)/PMED)**2)))   
249         IF(ZCE.LT.1E-4) ZCE=(PMTH(2,22)/PMED)**2    
250       ENDIF 
251       ZC=MIN(ZC,0.491)  
252       ZCE=MIN(ZCE,0.491)    
253       IF((MSTJ(41).EQ.1.AND.ZC.GT.0.49).OR.(MSTJ(41).EQ.2.AND.  
254      &MIN(ZC,ZCE).GT.0.49)) THEN    
255         P(IEP(1),5)=PMTH(1,KFL(1))  
256         V(IEP(1),5)=P(IEP(1),5)**2  
257         GOTO 300    
258       ENDIF 
259     
260 C...Integral of Altarelli-Parisi z kernel for QCD.  
261       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN   
262         FBR=6.*LOG((1.-ZC)/ZC)+MSTJ(45)*(0.5-ZC)    
263       ELSEIF(MSTJ(49).EQ.0) THEN    
264         FBR=(8./3.)*LOG((1.-ZC)/ZC) 
265     
266 C...Integral of Altarelli-Parisi z kernel for scalar gluon. 
267       ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN   
268         FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1.-2.*ZC) 
269       ELSEIF(MSTJ(49).EQ.1) THEN    
270         FBR=(1.-2.*ZC)/3.   
271         IF(IGM.EQ.0.AND.M3JC.EQ.1) FBR=4.*FBR   
272     
273 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon. 
274       ELSEIF(KFL(1).EQ.21) THEN 
275         FBR=6.*MSTJ(45)*(0.5-ZC)    
276       ELSE  
277         FBR=2.*LOG((1.-ZC)/ZC)  
278       ENDIF 
279     
280 C...Integral of Altarelli-Parisi kernel for photon emission.    
281       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) 
282      &FBRE=(KCHG(KFL(1),1)/3.)**2*2.*LOG((1.-ZCE)/ZCE)  
283     
284 C...Inner veto algorithm starts. Find maximum mass for evolution.   
285   260 PMS=V(IEP(1),5)   
286       IF(IGM.GE.0) THEN 
287         PM2=0.  
288         DO 270 I=2,NEP  
289         PM=P(IEP(I),5)  
290         IF(KFL(I).GT.0.AND.(KFL(I).LE.8.OR.KFL(I).EQ.21)) PM=   
291      &  PMTH(2,KFL(I))  
292   270   PM2=PM2+PM  
293         PMS=MIN(PMS,(P(IM,5)-PM2)**2)   
294       ENDIF 
295     
296 C...Select mass for daughter in QCD evolution.  
297       B0=27./6. 
298       DO 280 IF=4,MSTJ(45)  
299   280 IF(PMS.GT.4.*PMTH(2,IF)**2) B0=(33.-2.*IF)/6. 
300       IF(MSTJ(44).LE.0) THEN    
301          PMSQCD=PMS*EXP(MAX(-100.,LOG(RLU_HIJING(0))*PARU(2)/(PARU(111)
302      $        *FBR)))  
303       ELSEIF(MSTJ(44).EQ.1) THEN    
304         PMSQCD=4.*ALAMS*(0.25*PMS/ALAMS)**(RLU_HIJING(0)**(B0/FBR))    
305       ELSE  
306         PMSQCD=PMS*RLU_HIJING(0)**(ALFM*B0/FBR)    
307       ENDIF 
308       IF(ZC.GT.0.49.OR.PMSQCD.LE.PMTH(4,KFL(1))**2) PMSQCD= 
309      &PMTH(2,KFL(1))**2 
310       V(IEP(1),5)=PMSQCD    
311       MCE=1 
312     
313 C...Select mass for daughter in QED evolution.  
314       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) THEN    
315          PMSQED=PMS*EXP(MAX(-100.,LOG(RLU_HIJING(0))*PARU(2)/(PARU(101)
316      $        *FBRE))) 
317         IF(ZCE.GT.0.49.OR.PMSQED.LE.PMTH(5,KFL(1))**2) PMSQED=  
318      &  PMTH(2,KFL(1))**2   
319         IF(PMSQED.GT.PMSQCD) THEN   
320           V(IEP(1),5)=PMSQED    
321           MCE=2 
322         ENDIF   
323       ENDIF 
324     
325 C...Check whether daughter mass below cutoff.   
326       P(IEP(1),5)=SQRT(V(IEP(1),5)) 
327       IF(P(IEP(1),5).LE.PMTH(3,KFL(1))) THEN    
328         P(IEP(1),5)=PMTH(1,KFL(1))  
329         V(IEP(1),5)=P(IEP(1),5)**2  
330         GOTO 300    
331       ENDIF 
332     
333 C...Select z value of branching: q -> qgamma.   
334       IF(MCE.EQ.2) THEN 
335         Z=1.-(1.-ZCE)*(ZCE/(1.-ZCE))**RLU_HIJING(0)    
336         IF(1.+Z**2.LT.2.*RLU_HIJING(0)) GOTO 260   
337         K(IEP(1),5)=22  
338     
339 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.  
340       ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN   
341         Z=1.-(1.-ZC)*(ZC/(1.-ZC))**RLU_HIJING(0)   
342         IF(1.+Z**2.LT.2.*RLU_HIJING(0)) GOTO 260   
343         K(IEP(1),5)=21  
344       ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*(0.5-ZC).LT.RLU_HIJING(0)*FBR)
345      $       THEN    
346         Z=(1.-ZC)*(ZC/(1.-ZC))**RLU_HIJING(0)  
347         IF(RLU_HIJING(0).GT.0.5) Z=1.-Z    
348         IF((1.-Z*(1.-Z))**2.LT.RLU_HIJING(0)) GOTO 260 
349         K(IEP(1),5)=21  
350       ELSEIF(MSTJ(49).NE.1) THEN    
351         Z=ZC+(1.-2.*ZC)*RLU_HIJING(0)  
352         IF(Z**2+(1.-Z)**2.LT.RLU_HIJING(0)) GOTO 260   
353         KFLB=1+INT(MSTJ(45)*RLU_HIJING(0)) 
354         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)  
355         IF(PMQ.GE.1.) GOTO 260  
356         PMQ0=4.*PMTH(2,21)**2/V(IEP(1),5)   
357         IF(MOD(MSTJ(43),2).EQ.0.AND.(1.+0.5*PMQ)*SQRT(1.-PMQ).LT.   
358      &  RLU_HIJING(0)*(1.+0.5*PMQ0)*SQRT(1.-PMQ0)) GOTO 260    
359         K(IEP(1),5)=KFLB    
360     
361 C...Ditto for scalar gluon model.   
362       ELSEIF(KFL(1).NE.21) THEN 
363         Z=1.-SQRT(ZC**2+RLU_HIJING(0)*(1.-2.*ZC))  
364         K(IEP(1),5)=21  
365       ELSEIF(RLU_HIJING(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87))
366      $       THEN  
367         Z=ZC+(1.-2.*ZC)*RLU_HIJING(0)  
368         K(IEP(1),5)=21  
369       ELSE  
370         Z=ZC+(1.-2.*ZC)*RLU_HIJING(0)  
371         KFLB=1+INT(MSTJ(45)*RLU_HIJING(0)) 
372         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)  
373         IF(PMQ.GE.1.) GOTO 260  
374         K(IEP(1),5)=KFLB    
375       ENDIF 
376       IF(MCE.EQ.1.AND.MSTJ(44).GE.2) THEN   
377         IF(Z*(1.-Z)*V(IEP(1),5).LT.PT2MIN) GOTO 260 
378         IF(ALFM/LOG(V(IEP(1),5)*Z*(1.-Z)/ALAMS).LT.RLU_HIJING(0)) GOTO
379      $       260 
380       ENDIF 
381     
382 C...Check if z consistent with chosen m.    
383       IF(KFL(1).EQ.21) THEN 
384         KFLGD1=IABS(K(IEP(1),5))    
385         KFLGD2=KFLGD1   
386       ELSE  
387         KFLGD1=KFL(1)   
388         KFLGD2=IABS(K(IEP(1),5))    
389       ENDIF 
390       IF(NEP.EQ.1) THEN 
391         PED=PS(4)   
392       ELSEIF(NEP.GE.3) THEN 
393         PED=P(IEP(1),4) 
394       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN    
395         PED=0.5*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)    
396       ELSE  
397         IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM   
398         IF(IEP(1).EQ.N+2) PED=(1.-V(IM,1))*PEM  
399       ENDIF 
400       IF(MOD(MSTJ(43),2).EQ.1) THEN 
401         PMQTH3=0.5*PARJ(82) 
402         IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)    
403         PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(IEP(1),5)  
404         PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(IEP(1),5)  
405         ZD=SQRT(MAX(0.,(1.-V(IEP(1),5)/PED**2)*((1.-PMQ1-PMQ2)**2-  
406      &  4.*PMQ1*PMQ2))) 
407         ZH=1.+PMQ1-PMQ2 
408       ELSE  
409         ZD=SQRT(MAX(0.,1.-V(IEP(1),5)/PED**2))  
410         ZH=1.   
411       ENDIF 
412       ZL=0.5*(ZH-ZD)    
413       ZU=0.5*(ZH+ZD)    
414       IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 260   
415       IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL* 
416      &(1.-ZU))) 
417       IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))    
418     
419 C...Three-jet matrix element correction.    
420       IF(IGM.EQ.0.AND.M3JC.EQ.1) THEN   
421         X1=Z*(1.+V(IEP(1),5)/V(NS+1,5)) 
422         X2=1.-V(IEP(1),5)/V(NS+1,5) 
423         X3=(1.-X1)+(1.-X2)  
424         IF(MCE.EQ.2) THEN   
425           KI1=K(IPA(INUM),2)    
426           KI2=K(IPA(3-INUM),2)  
427           QF1=KCHG(IABS(KI1),1)*ISIGN(1,KI1)/3. 
428           QF2=KCHG(IABS(KI2),1)*ISIGN(1,KI2)/3. 
429           WSHOW=QF1**2*(1.-X1)/X3*(1.+(X1/(2.-X2))**2)+ 
430      &    QF2**2*(1.-X2)/X3*(1.+(X2/(2.-X1))**2)    
431           WME=(QF1*(1.-X1)/X3-QF2*(1.-X2)/X3)**2*(X1**2+X2**2)  
432         ELSEIF(MSTJ(49).NE.1) THEN  
433           WSHOW=1.+(1.-X1)/X3*(X1/(2.-X2))**2+  
434      &    (1.-X2)/X3*(X2/(2.-X1))**2    
435           WME=X1**2+X2**2   
436         ELSE    
437           WSHOW=4.*X3*((1.-X1)/(2.-X2)**2+(1.-X2)/(2.-X1)**2)   
438           WME=X3**2 
439         ENDIF   
440         IF(WME.LT.RLU_HIJING(0)*WSHOW) GOTO 260    
441     
442 C...Impose angular ordering by rejection of nonordered emission.    
443       ELSEIF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2) THEN  
444         MAOM=1  
445         ZM=V(IM,1)  
446         IF(IEP(1).EQ.N+2) ZM=1.-V(IM,1) 
447         THE2ID=Z*(1.-Z)*(ZM*P(IM,4))**2/V(IEP(1),5) 
448         IAOM=IM 
449   290   IF(K(IAOM,5).EQ.22) THEN    
450           IAOM=K(IAOM,3)    
451           IF(K(IAOM,3).LE.NS) MAOM=0    
452           IF(MAOM.EQ.1) GOTO 290    
453         ENDIF   
454         IF(MAOM.EQ.1) THEN  
455           THE2IM=V(IAOM,1)*(1.-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)    
456           IF(THE2ID.LT.THE2IM) GOTO 260 
457         ENDIF   
458       ENDIF 
459     
460 C...Impose user-defined maximum angle at first branching.   
461       IF(MSTJ(48).EQ.1) THEN    
462         IF(NEP.EQ.1.AND.IM.EQ.NS) THEN  
463           THE2ID=Z*(1.-Z)*PS(4)**2/V(IEP(1),5)  
464           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260 
465         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN    
466           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)  
467           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260 
468         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN    
469           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)  
470           IF(THE2ID.LT.1./PARJ(86)**2) GOTO 260 
471         ENDIF   
472       ENDIF 
473     
474 C...End of inner veto algorithm. Check if only one leg evolved so far.  
475   300 V(IEP(1),1)=Z 
476       ISL(1)=0  
477       ISL(2)=0  
478       IF(NEP.EQ.1) GOTO 330 
479       IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 200  
480       DO 310 I=1,NEP    
481       IF(ITRY(I).EQ.0.AND.KFLD(I).GT.0.AND.(KFLD(I).LE.8.OR.KFLD(I).EQ. 
482      &21)) THEN 
483         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) GOTO 200    
484       ENDIF 
485   310 CONTINUE  
486     
487 C...Check if chosen multiplet m1,m2,z1,z2 is physical.  
488       IF(NEP.EQ.3) THEN 
489         PA1S=(P(N+1,4)+P(N+1,5))*(P(N+1,4)-P(N+1,5))    
490         PA2S=(P(N+2,4)+P(N+2,5))*(P(N+2,4)-P(N+2,5))    
491         PA3S=(P(N+3,4)+P(N+3,5))*(P(N+3,4)-P(N+3,5))    
492         PTS=0.25*(2.*PA1S*PA2S+2.*PA1S*PA3S+2.*PA2S*PA3S-   
493      &  PA1S**2-PA2S**2-PA3S**2)/PA1S   
494         IF(PTS.LE.0.) GOTO 200  
495       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN    
496         DO 320 I1=N+1,N+2   
497         KFLDA=IABS(K(I1,2)) 
498         IF(KFLDA.EQ.0.OR.(KFLDA.GT.8.AND.KFLDA.NE.21)) GOTO 320 
499         IF(P(I1,5).LT.PMTH(2,KFLDA)) GOTO 320   
500         IF(KFLDA.EQ.21) THEN    
501           KFLGD1=IABS(K(I1,5))  
502           KFLGD2=KFLGD1 
503         ELSE    
504           KFLGD1=KFLDA  
505           KFLGD2=IABS(K(I1,5))  
506         ENDIF   
507         I2=2*N+3-I1 
508         IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN  
509           PED=0.5*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5) 
510         ELSE    
511           IF(I1.EQ.N+1) ZM=V(IM,1)  
512           IF(I1.EQ.N+2) ZM=1.-V(IM,1)   
513           PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-  
514      &    4.*V(N+1,5)*V(N+2,5)) 
515           PED=PEM*(0.5*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/V(IM,5)    
516         ENDIF   
517         IF(MOD(MSTJ(43),2).EQ.1) THEN   
518           PMQTH3=0.5*PARJ(82)   
519           IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)  
520           PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(I1,5)    
521           PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(I1,5)    
522           ZD=SQRT(MAX(0.,(1.-V(I1,5)/PED**2)*((1.-PMQ1-PMQ2)**2-    
523      &    4.*PMQ1*PMQ2)))   
524           ZH=1.+PMQ1-PMQ2   
525         ELSE    
526           ZD=SQRT(MAX(0.,1.-V(I1,5)/PED**2))    
527           ZH=1. 
528         ENDIF   
529         ZL=0.5*(ZH-ZD)  
530         ZU=0.5*(ZH+ZD)  
531         IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(1)=1 
532         IF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(2)=1 
533         IF(KFLDA.EQ.21) V(I1,4)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*(1.-ZU)))   
534         IF(KFLDA.NE.21) V(I1,4)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))   
535   320   CONTINUE    
536         IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN  
537           ISL(3-ISLM)=0 
538           ISLM=3-ISLM   
539         ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN    
540           ZDR1=MAX(0D0,V(N+1,3)/MAX(1D-6,V(N+1,4))-1D0)
541           ZDR2=MAX(0D0,V(N+2,3)/MAX(1D-6,V(N+2,4))-1D0)
542           IF(ZDR2.GT.RLU_HIJING(0)*(ZDR1+ZDR2)) ISL(1)=0   
543           IF(ISL(1).EQ.1) ISL(2)=0  
544           IF(ISL(1).EQ.0) ISLM=1    
545           IF(ISL(2).EQ.0) ISLM=2    
546         ENDIF   
547         IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 200 
548       ENDIF 
549       IF(IGM.GT.0.AND.MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.    
550      &PMTH(2,KFLD(1)).OR.P(N+2,5).GE.PMTH(2,KFLD(2)))) THEN 
551         PMQ1=V(N+1,5)/V(IM,5)   
552         PMQ2=V(N+2,5)/V(IM,5)   
553         ZD=SQRT(MAX(0.,(1.-V(IM,5)/PEM**2)*((1.-PMQ1-PMQ2)**2-  
554      &  4.*PMQ1*PMQ2))) 
555         ZH=1.+PMQ1-PMQ2 
556         ZL=0.5*(ZH-ZD)  
557         ZU=0.5*(ZH+ZD)  
558         IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 200 
559       ENDIF 
560     
561 C...Accepted branch. Construct four-momentum for initial partons.   
562   330 MAZIP=0   
563       MAZIC=0   
564       IF(NEP.EQ.1) THEN 
565         P(N+1,1)=0. 
566         P(N+1,2)=0. 
567         P(N+1,3)=SQRT(MAX(0.,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-   
568      &  P(N+1,5)))) 
569         P(N+1,4)=P(IPA(1),4)    
570         V(N+1,2)=P(N+1,4)   
571       ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN    
572         PED1=0.5*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)    
573         P(N+1,1)=0. 
574         P(N+1,2)=0. 
575         P(N+1,3)=SQRT(MAX(0.,(PED1+P(N+1,5))*(PED1-P(N+1,5))))  
576         P(N+1,4)=PED1   
577         P(N+2,1)=0. 
578         P(N+2,2)=0. 
579         P(N+2,3)=-P(N+1,3)  
580         P(N+2,4)=P(IM,5)-PED1   
581         V(N+1,2)=P(N+1,4)   
582         V(N+2,2)=P(N+2,4)   
583       ELSEIF(NEP.EQ.3) THEN 
584         P(N+1,1)=0. 
585         P(N+1,2)=0. 
586         P(N+1,3)=SQRT(MAX(0.,PA1S)) 
587         P(N+2,1)=SQRT(PTS)  
588         P(N+2,2)=0. 
589         P(N+2,3)=0.5*(PA3S-PA2S-PA1S)/P(N+1,3)  
590         P(N+3,1)=-P(N+2,1)  
591         P(N+3,2)=0. 
592         P(N+3,3)=-(P(N+1,3)+P(N+2,3))   
593         V(N+1,2)=P(N+1,4)   
594         V(N+2,2)=P(N+2,4)   
595         V(N+3,2)=P(N+3,4)   
596     
597 C...Construct transverse momentum for ordinary branching in shower. 
598       ELSE  
599         ZM=V(IM,1)  
600         PZM=SQRT(MAX(0.,(PEM+P(IM,5))*(PEM-P(IM,5))))   
601         PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4.*V(N+1,5)*V(N+2,5)    
602         IF(PZM.LE.0.) THEN  
603           PTS=0.    
604         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN   
605           PTS=(PEM**2*(ZM*(1.-ZM)*V(IM,5)-(1.-ZM)*V(N+1,5)- 
606      &    ZM*V(N+2,5))-0.25*PMLS)/PZM**2    
607         ELSE    
608           PTS=PMLS*(ZM*(1.-ZM)*PEM**2/V(IM,5)-0.25)/PZM**2  
609         ENDIF   
610         PT=SQRT(MAX(0.,PTS))    
611     
612 C...Find coefficient of azimuthal asymmetry due to gluon polarization.  
613         HAZIP=0.    
614         IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21.    
615      &  AND.IAU.NE.0) THEN  
616           IF(K(IGM,3).NE.0) MAZIP=1 
617           ZAU=V(IGM,1)  
618           IF(IAU.EQ.IM+1) ZAU=1.-V(IGM,1)   
619           IF(MAZIP.EQ.0) ZAU=0. 
620           IF(K(IGM,2).NE.21) THEN   
621             HAZIP=2.*ZAU/(1.+ZAU**2)    
622           ELSE  
623             HAZIP=(ZAU/(1.-ZAU*(1.-ZAU)))**2    
624           ENDIF 
625           IF(K(N+1,2).NE.21) THEN   
626             HAZIP=HAZIP*(-2.*ZM*(1.-ZM))/(1.-2.*ZM*(1.-ZM)) 
627           ELSE  
628             HAZIP=HAZIP*(ZM*(1.-ZM)/(1.-ZM*(1.-ZM)))**2 
629           ENDIF 
630         ENDIF   
631     
632 C...Find coefficient of azimuthal asymmetry due to soft gluon   
633 C...interference.   
634         HAZIC=0.    
635         IF(MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.K(N+2,2).EQ.21).    
636      &  AND.IAU.NE.0) THEN  
637           IF(K(IGM,3).NE.0) MAZIC=N+1   
638           IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2    
639           IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.   
640      &    ZM.GT.0.5) MAZIC=N+2  
641           IF(K(IAU,2).EQ.22) MAZIC=0    
642           ZS=ZM 
643           IF(MAZIC.EQ.N+2) ZS=1.-ZM 
644           ZGM=V(IGM,1)  
645           IF(IAU.EQ.IM-1) ZGM=1.-V(IGM,1)   
646           IF(MAZIC.EQ.0) ZGM=1. 
647           HAZIC=(P(IM,5)/P(IGM,5))*SQRT((1.-ZS)*(1.-ZGM)/(ZS*ZGM))  
648           HAZIC=MIN(0.95,HAZIC) 
649         ENDIF   
650       ENDIF 
651     
652 C...Construct kinematics for ordinary branching in shower.  
653   340 IF(NEP.EQ.2.AND.IGM.GT.0) THEN    
654         IF(MOD(MSTJ(43),2).EQ.1) THEN   
655           P(N+1,4)=PEM*V(IM,1)  
656         ELSE    
657           P(N+1,4)=PEM*(0.5*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+ 
658      &    SQRT(PMLS)*ZM)/V(IM,5)    
659         ENDIF   
660         PHI=PARU(2)*RLU_HIJING(0)  
661         P(N+1,1)=PT*COS(PHI)    
662         P(N+1,2)=PT*SIN(PHI)    
663         IF(PZM.GT.0.) THEN  
664           P(N+1,3)=0.5*(V(N+2,5)-V(N+1,5)-V(IM,5)+2.*PEM*P(N+1,4))/PZM  
665         ELSE    
666           P(N+1,3)=0.   
667         ENDIF   
668         P(N+2,1)=-P(N+1,1)  
669         P(N+2,2)=-P(N+1,2)  
670         P(N+2,3)=PZM-P(N+1,3)   
671         P(N+2,4)=PEM-P(N+1,4)   
672         IF(MSTJ(43).LE.2) THEN  
673           V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)  
674           V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)  
675         ENDIF   
676       ENDIF 
677     
678 C...Rotate and boost daughters. 
679       IF(IGM.GT.0) THEN 
680         IF(MSTJ(43).LE.2) THEN  
681           BEX=P(IGM,1)/P(IGM,4) 
682           BEY=P(IGM,2)/P(IGM,4) 
683           BEZ=P(IGM,3)/P(IGM,4) 
684           GA=P(IGM,4)/P(IGM,5)  
685           GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1.+GA)-   
686      &    P(IM,4))  
687         ELSE    
688           BEX=0.    
689           BEY=0.    
690           BEZ=0.    
691           GA=1. 
692           GABEP=0.  
693         ENDIF   
694         THE=ULANGL_HIJING(P(IM,3)+GABEP*BEZ,SQRT((P(IM,1)+GABEP*BEX)**2+   
695      &  (P(IM,2)+GABEP*BEY)**2))    
696         PHI=ULANGL_HIJING(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY) 
697         DO 350 I=N+1,N+2    
698         DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+ 
699      &  SIN(THE)*COS(PHI)*P(I,3)    
700         DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+ 
701      &  SIN(THE)*SIN(PHI)*P(I,3)    
702         DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)  
703         DP(4)=P(I,4)    
704         DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)   
705         DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))    
706         P(I,1)=DP(1)+DGABP*BEX  
707         P(I,2)=DP(2)+DGABP*BEY  
708         P(I,3)=DP(3)+DGABP*BEZ  
709   350   P(I,4)=GA*(DP(4)+DBP)   
710       ENDIF 
711     
712 C...Weight with azimuthal distribution, if required.    
713       IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN 
714         DO 360 J=1,3    
715         DPT(1,J)=P(IM,J)    
716         DPT(2,J)=P(IAU,J)   
717   360   DPT(3,J)=P(N+1,J)   
718         DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)  
719         DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)  
720         DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2    
721         DO 370 J=1,3    
722         DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/DPMM    
723   370   DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/DPMM    
724         DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)  
725         DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)  
726         IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN 
727           CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+ 
728      &    DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))    
729           IF(MAZIP.NE.0) THEN   
730             IF(1.+HAZIP*(2.*CAD**2-1.).LT.RLU_HIJING(0)*(1.+ABS(HAZIP)))   
731      &      GOTO 340    
732           ENDIF 
733           IF(MAZIC.NE.0) THEN   
734             IF(MAZIC.EQ.N+2) CAD=-CAD   
735             IF((1.-HAZIC)*(1.-HAZIC*CAD)/(1.+HAZIC**2-2.*HAZIC*CAD).    
736      &      LT.RLU_HIJING(0)) GOTO 340 
737           ENDIF 
738         ENDIF   
739       ENDIF 
740     
741 C...Continue loop over partons that may branch, until none left.    
742       IF(IGM.GE.0) K(IM,1)=14   
743       N=N+NEP   
744       NEP=2 
745       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN  
746          CALL LUERRM_HIJING(11
747      $        ,'(LUSHOW_HIJING:) no more memory left in LUJETS_HIJING')   
748         IF(MSTU(21).GE.1) N=NS  
749         IF(MSTU(21).GE.1) RETURN    
750       ENDIF 
751       GOTO 140  
752     
753 C...Set information on imagined shower initiator.   
754   380 IF(NPA.GE.2) THEN 
755         K(NS+1,1)=11    
756         K(NS+1,2)=94    
757         K(NS+1,3)=IP1   
758         IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2   
759         K(NS+1,4)=NS+2  
760         K(NS+1,5)=NS+1+NPA  
761         IIM=1   
762       ELSE  
763         IIM=0   
764       ENDIF 
765     
766 C...Reconstruct string drawing information. 
767       DO 390 I=NS+1+IIM,N   
768       IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN    
769         K(I,1)=1    
770       ELSEIF(K(I,1).LE.10) THEN 
771         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5)) 
772         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5)) 
773       ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN 
774         ID1=MOD(K(I,4),MSTU(5)) 
775         IF(K(I,2).GE.1.AND.K(I,2).LE.8) ID1=MOD(K(I,4),MSTU(5))+1   
776         ID2=2*MOD(K(I,4),MSTU(5))+1-ID1 
777         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1 
778         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2 
779         K(ID1,4)=K(ID1,4)+MSTU(5)*I 
780         K(ID1,5)=K(ID1,5)+MSTU(5)*ID2   
781         K(ID2,4)=K(ID2,4)+MSTU(5)*ID1   
782         K(ID2,5)=K(ID2,5)+MSTU(5)*I 
783       ELSE  
784         ID1=MOD(K(I,4),MSTU(5)) 
785         ID2=ID1+1   
786         K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1 
787         K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1 
788         K(ID1,4)=K(ID1,4)+MSTU(5)*I 
789         K(ID1,5)=K(ID1,5)+MSTU(5)*I 
790         K(ID2,4)=0  
791         K(ID2,5)=0  
792       ENDIF 
793   390 CONTINUE  
794     
795 C...Transformation from CM frame.   
796       IF(NPA.GE.2) THEN 
797         BEX=PS(1)/PS(4) 
798         BEY=PS(2)/PS(4) 
799         BEZ=PS(3)/PS(4) 
800         GA=PS(4)/PS(5)  
801         GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))  
802      &  /(1.+GA)-P(IPA(1),4))   
803       ELSE  
804         BEX=0.  
805         BEY=0.  
806         BEZ=0.  
807         GABEP=0.    
808       ENDIF 
809       THE=ULANGL_HIJING(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)    
810      &+GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))   
811       PHI=ULANGL_HIJING(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)   
812       IF(NPA.EQ.3) THEN 
813          CHI=ULANGL_HIJING(COS(THE)*COS(PHI)*(P(IPA(2),1)+GABEP*BEX)
814      $        +COS(THE)*SIN(PHI)*(P(IPA(2),2)+GABEP*BEY)-SIN(THE)
815      $        *(P(IPA(2),3)+GABEP*BEZ),-SIN(PHI)*(P(IPA(2),1)+GABEP*BEX)
816      $        +COS(PHI)*(P(IPA(2),2)+GABEP*BEY)) 
817         CALL LUDBRB_HIJING(NS+1,N,0.,CHI,0D0,0D0,0D0)  
818       ENDIF 
819       DBEX=DBLE(BEX)    
820       DBEY=DBLE(BEY)    
821       DBEZ=DBLE(BEZ)    
822       CALL LUDBRB_HIJING(NS+1,N,THE,PHI,DBEX,DBEY,DBEZ)    
823     
824 C...Decay vertex of shower. 
825       DO 400 I=NS+1,N   
826       DO 400 J=1,5  
827   400 V(I,J)=V(IP1,J)   
828     
829 C...Delete trivial shower, else connect initiators. 
830       IF(N.EQ.NS+NPA+IIM) THEN  
831         N=NS    
832       ELSE  
833         DO 410 IP=1,NPA 
834         K(IPA(IP),1)=14 
835         K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP 
836         K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP 
837         K(NS+IIM+IP,3)=IPA(IP)  
838         IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1  
839         K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)   
840   410   K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)   
841       ENDIF 
842     
843       RETURN    
844       END