]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/lushow_hijing.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Using EINCLUDE instead of EXPORT to access HIJING/AliHijingRndm.h from THijing
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / lushow_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE LUSHOW_HIJING(IP1,IP2,QMAX)
6
7C...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
15C...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
42C...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
63C...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
88C...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
110C...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
130C...Position of aunt (sister to branching parton).
131C...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
158C...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
170C...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
192C...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
215C...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
233C...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
260C...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
266C...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
273C...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
280C...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
284C...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
296C...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
313C...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
325C...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
333C...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
339C...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
361C...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
382C...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
419C...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
442C...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
460C...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
474C...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
487C...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(0.,V(N+1,3)/V(N+1,4)-1.)
541 ZDR2=MAX(0.,V(N+2,3)/V(N+2,4)-1.)
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
561C...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
597C...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
612C...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
632C...Find coefficient of azimuthal asymmetry due to soft gluon
633C...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
652C...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
678C...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
712C...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
741C...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
753C...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
766C...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
795C...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
824C...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
829C...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