1 *CMZ : 17/07/98 15.44.33 by Federico Carminati
3 C*********************************************************************
5 SUBROUTINE LUSHOW(IP1,IP2,QMAX)
7 C...Purpose: to generate timelike parton showers from given partons.
8 IMPLICIT DOUBLE PRECISION(D)
10 COMMON /LUJETS/ N,K(200000,5),P(200000,5),V(200000,5)
13 COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
16 COMMON /LUDAT2/ KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
19 DIMENSION PMTH(5,40),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4),
20 &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4)
22 C...Initialization of cutoff masses etc.
23 IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR.
24 &QMAX.LE.MIN(PARJ(82),PARJ(83)).OR.MSTJ(41).GE.3) RETURN
26 PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25*PARJ(82)**2)
27 PMTH(3,21)=2.*PMTH(2,21)
31 PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25*PARJ(83)**2)
32 PMTH(3,22)=2.*PMTH(2,22)
36 IF(MSTJ(41).EQ.2) PMQTH1=MIN(PARJ(82),PARJ(83))
38 IF(MSTJ(41).EQ.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
41 PMTH(2,IF)=SQRT(PMTH(1,IF)**2+0.25*PMQTH1**2)
42 PMTH(3,IF)=PMTH(2,IF)+PMQTH2
43 PMTH(4,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(82)**2)+PMTH(2,21)
44 100 PMTH(5,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(83)**2)+PMTH(2,22)
45 PT2MIN=MAX(0.5*PARJ(82),1.1*PARJ(81))**2
47 ALFM=LOG(PT2MIN/ALAMS)
49 C...Store positions of shower initiating partons.
51 IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
54 ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
59 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0.
66 & '(LUSHOW:) failed to reconstruct showering system')
67 IF(MSTU(21).GE.1) RETURN
70 C...Check on phase space available for emission.
76 KFLA(I)=IABS(K(IPA(I),2))
78 IF(KFLA(I).NE.0.AND.(KFLA(I).LE.8.OR.KFLA(I).EQ.21))
79 &PMA(I)=PMTH(3,KFLA(I))
81 IF(KFLA(I).EQ.0.OR.(KFLA(I).GT.8.AND.KFLA(I).NE.21).OR.
82 &PMA(I).GT.QMAX) IREJ=IREJ+1
84 130 PS(J)=PS(J)+P(IPA(I),J)
85 IF(IREJ.EQ.NPA) RETURN
86 PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
87 IF(NPA.EQ.1) PS(5)=PS(4)
88 IF(PS(5).LE.PM+PMQTH1) RETURN
89 IF(NPA.EQ.2.AND.MSTJ(47).GE.1) THEN
90 IF(KFLA(1).GE.1.AND.KFLA(1).LE.8.AND.KFLA(2).GE.1.AND.
91 & KFLA(2).LE.8) M3JC=1
92 IF(MSTJ(47).GE.2) M3JC=1
95 C...Define imagined single initiator of shower for parton system.
97 IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
98 CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
99 IF(MSTU(21).GE.1) RETURN
116 C...Loop over partons that may branch.
124 IF(KFLM.EQ.0.OR.(KFLM.GT.8.AND.KFLM.NE.21)) GOTO 140
125 IF(P(IM,5).LT.PMTH(2,KFLM)) GOTO 140
130 IF(N+NEP.GT.MSTU(4)-MSTU(32)-5) THEN
131 CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
132 IF(MSTU(21).GE.1) RETURN
135 C...Position of aunt (sister to branching parton).
136 C...Origin and flavour of daughters.
139 IF(K(IM-1,3).EQ.IGM) IAU=IM-1
140 IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
151 160 K(N+I,2)=K(IPA(I),2)
152 ELSEIF(KFLM.NE.21) THEN
155 ELSEIF(K(IM,5).EQ.21) THEN
163 C...Reset flags on daughers and tries made.
168 KFLD(IP)=IABS(K(N+IP,2))
172 170 IF(KFLD(IP).GT.0.AND.(KFLD(IP).LE.8.OR.KFLD(IP).EQ.21)) ISI(IP)=1
175 C...Maximum virtuality of daughters.
178 IF(NPA.GE.3) P(N+I,4)=(PS(4)*P(IPA(I),4)-PS(1)*P(IPA(I),1)-
179 & PS(2)*P(IPA(I),2)-PS(3)*P(IPA(I),3))/PS(5)
180 P(N+I,5)=MIN(QMAX,PS(5))
181 IF(NPA.GE.3) P(N+I,5)=MIN(P(N+I,5),P(N+I,4))
182 180 IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
184 IF(MSTJ(43).LE.2) PEM=V(IM,2)
185 IF(MSTJ(43).GE.3) PEM=P(IM,4)
186 P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
187 P(N+2,5)=MIN(P(IM,5),(1.-V(IM,1))*PEM)
188 IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
193 IF(P(N+I,5).LE.PMTH(3,KFLD(I))) P(N+I,5)=PMTH(1,KFLD(I))
195 190 V(N+I,5)=P(N+I,5)**2
197 C...Choose one of the daughters for evolution.
201 210 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
203 IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
204 IF(P(N+I,5).GE.PMTH(2,KFLD(I))) INUM=I
210 IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQTH2) THEN
212 IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,KFLD(I))) THEN
220 C...Store information on choice of evolving daughter.
225 240 IF(IEP(I).GT.N+NEP) IEP(I)=N+1
227 250 KFL(I)=IABS(K(IEP(I),2))
228 ITRY(INUM)=ITRY(INUM)+1
229 IF(ITRY(INUM).GT.200) THEN
230 CALL LUERRM(14,'(LUSHOW:) caught in infinite loop')
231 IF(MSTU(21).GE.1) RETURN
234 IF(KFL(1).EQ.0.OR.(KFL(1).GT.8.AND.KFL(1).NE.21)) GOTO 300
235 IF(P(IEP(1),5).LT.PMTH(2,KFL(1))) GOTO 300
237 C...Calculate allowed z range.
240 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
243 IF(INUM.EQ.1) PMED=V(IM,1)*PEM
244 IF(INUM.EQ.2) PMED=(1.-V(IM,1))*PEM
246 IF(MOD(MSTJ(43),2).EQ.1) THEN
250 ZC=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,21)/PMED)**2)))
251 IF(ZC.LT.1E-4) ZC=(PMTH(2,21)/PMED)**2
252 ZCE=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,22)/PMED)**2)))
253 IF(ZCE.LT.1E-4) ZCE=(PMTH(2,22)/PMED)**2
257 IF((MSTJ(41).EQ.1.AND.ZC.GT.0.49).OR.(MSTJ(41).EQ.2.AND.
258 &MIN(ZC,ZCE).GT.0.49)) THEN
259 P(IEP(1),5)=PMTH(1,KFL(1))
260 V(IEP(1),5)=P(IEP(1),5)**2
264 C...Integral of Altarelli-Parisi z kernel for QCD.
265 IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
266 FBR=6.*LOG((1.-ZC)/ZC)+MSTJ(45)*(0.5-ZC)
267 ELSEIF(MSTJ(49).EQ.0) THEN
268 FBR=(8./3.)*LOG((1.-ZC)/ZC)
270 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
271 ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
272 FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1.-2.*ZC)
273 ELSEIF(MSTJ(49).EQ.1) THEN
275 IF(IGM.EQ.0.AND.M3JC.EQ.1) FBR=4.*FBR
277 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
278 ELSEIF(KFL(1).EQ.21) THEN
279 FBR=6.*MSTJ(45)*(0.5-ZC)
281 FBR=2.*LOG((1.-ZC)/ZC)
284 C...Integral of Altarelli-Parisi kernel for photon emission.
285 IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8)
286 &FBRE=(KCHG(KFL(1),1)/3.)**2*2.*LOG((1.-ZCE)/ZCE)
288 C...Inner veto algorithm starts. Find maximum mass for evolution.
294 IF(KFL(I).GT.0.AND.(KFL(I).LE.8.OR.KFL(I).EQ.21)) PM=
297 PMS=MIN(PMS,(P(IM,5)-PM2)**2)
300 C...Select mass for daughter in QCD evolution.
303 280 IF(PMS.GT.4.*PMTH(2,IF)**2) B0=(33.-2.*IF)/6.
304 IF(MSTJ(44).LE.0) THEN
305 PMSQCD=PMS*EXP(MAX(-80.,LOG(RLU(0))*PARU(2)/(PARU(111)*FBR)))
306 ELSEIF(MSTJ(44).EQ.1) THEN
307 PMSQCD=4.*ALAMS*(0.25*PMS/ALAMS)**(RLU(0)**(B0/FBR))
309 PMSQCD=PMS*RLU(0)**(ALFM*B0/FBR)
311 IF(ZC.GT.0.49.OR.PMSQCD.LE.PMTH(4,KFL(1))**2) PMSQCD=
316 C...Select mass for daughter in QED evolution.
317 IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) THEN
318 PMSQED=PMS*EXP(MAX(-80.,LOG(RLU(0))*PARU(2)/(PARU(101)*FBRE)))
319 IF(ZCE.GT.0.49.OR.PMSQED.LE.PMTH(5,KFL(1))**2) PMSQED=
321 IF(PMSQED.GT.PMSQCD) THEN
327 C...Check whether daughter mass below cutoff.
328 P(IEP(1),5)=SQRT(V(IEP(1),5))
329 IF(P(IEP(1),5).LE.PMTH(3,KFL(1))) THEN
330 P(IEP(1),5)=PMTH(1,KFL(1))
331 V(IEP(1),5)=P(IEP(1),5)**2
335 C...Select z value of branching: q -> qgamma.
337 Z=1.-(1.-ZCE)*(ZCE/(1.-ZCE))**RLU(0)
338 IF(1.+Z**2.LT.2.*RLU(0)) GOTO 260
341 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
342 ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
343 Z=1.-(1.-ZC)*(ZC/(1.-ZC))**RLU(0)
344 IF(1.+Z**2.LT.2.*RLU(0)) GOTO 260
346 ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*(0.5-ZC).LT.RLU(0)*FBR) THEN
347 Z=(1.-ZC)*(ZC/(1.-ZC))**RLU(0)
348 IF(RLU(0).GT.0.5) Z=1.-Z
349 IF((1.-Z*(1.-Z))**2.LT.RLU(0)) GOTO 260
351 ELSEIF(MSTJ(49).NE.1) THEN
352 Z=ZC+(1.-2.*ZC)*RLU(0)
353 IF(Z**2+(1.-Z)**2.LT.RLU(0)) GOTO 260
354 KFLB=1+INT(MSTJ(45)*RLU(0))
355 PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
356 IF(PMQ.GE.1.) GOTO 260
357 PMQ0=4.*PMTH(2,21)**2/V(IEP(1),5)
358 IF(MOD(MSTJ(43),2).EQ.0.AND.(1.+0.5*PMQ)*SQRT(1.-PMQ).LT.
359 & RLU(0)*(1.+0.5*PMQ0)*SQRT(1.-PMQ0)) GOTO 260
362 C...Ditto for scalar gluon model.
363 ELSEIF(KFL(1).NE.21) THEN
364 Z=1.-SQRT(ZC**2+RLU(0)*(1.-2.*ZC))
366 ELSEIF(RLU(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
367 Z=ZC+(1.-2.*ZC)*RLU(0)
370 Z=ZC+(1.-2.*ZC)*RLU(0)
371 KFLB=1+INT(MSTJ(45)*RLU(0))
372 PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
373 IF(PMQ.GE.1.) GOTO 260
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(0)) GOTO 260
381 C...Check if z consistent with chosen m.
382 IF(KFL(1).EQ.21) THEN
383 KFLGD1=IABS(K(IEP(1),5))
387 KFLGD2=IABS(K(IEP(1),5))
391 ELSEIF(NEP.GE.3) THEN
393 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
394 PED=0.5*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
396 IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
397 IF(IEP(1).EQ.N+2) PED=(1.-V(IM,1))*PEM
399 IF(MOD(MSTJ(43),2).EQ.1) THEN
401 IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
402 PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(IEP(1),5)
403 PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(IEP(1),5)
404 ZD=SQRT(MAX(0.,(1.-V(IEP(1),5)/PED**2)*((1.-PMQ1-PMQ2)**2-
408 ZD=SQRT(MAX(0.,1.-V(IEP(1),5)/PED**2))
413 IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 260
414 IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*
416 IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
418 C...Three-jet matrix element correction.
419 IF(IGM.EQ.0.AND.M3JC.EQ.1) THEN
420 X1=Z*(1.+V(IEP(1),5)/V(NS+1,5))
421 X2=1.-V(IEP(1),5)/V(NS+1,5)
426 QF1=KCHG(IABS(KI1),1)*ISIGN(1,KI1)/3.
427 QF2=KCHG(IABS(KI2),1)*ISIGN(1,KI2)/3.
428 WSHOW=QF1**2*(1.-X1)/X3*(1.+(X1/(2.-X2))**2)+
429 & QF2**2*(1.-X2)/X3*(1.+(X2/(2.-X1))**2)
430 WME=(QF1*(1.-X1)/X3-QF2*(1.-X2)/X3)**2*(X1**2+X2**2)
431 ELSEIF(MSTJ(49).NE.1) THEN
432 WSHOW=1.+(1.-X1)/X3*(X1/(2.-X2))**2+
433 & (1.-X2)/X3*(X2/(2.-X1))**2
436 WSHOW=4.*X3*((1.-X1)/(2.-X2)**2+(1.-X2)/(2.-X1)**2)
439 IF(WME.LT.RLU(0)*WSHOW) GOTO 260
441 C...Impose angular ordering by rejection of nonordered emission.
442 ELSEIF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2) THEN
445 IF(IEP(1).EQ.N+2) ZM=1.-V(IM,1)
446 THE2ID=Z*(1.-Z)*(ZM*P(IM,4))**2/V(IEP(1),5)
448 290 IF(K(IAOM,5).EQ.22) THEN
450 IF(K(IAOM,3).LE.NS) MAOM=0
451 IF(MAOM.EQ.1) GOTO 290
454 THE2IM=V(IAOM,1)*(1.-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
455 IF(THE2ID.LT.THE2IM) GOTO 260
459 C...Impose user-defined maximum angle at first branching.
460 IF(MSTJ(48).EQ.1) THEN
461 IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
462 THE2ID=Z*(1.-Z)*PS(4)**2/V(IEP(1),5)
463 IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260
464 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
465 THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
466 IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260
467 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
468 THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
469 IF(THE2ID.LT.1./PARJ(86)**2) GOTO 260
473 C...End of inner veto algorithm. Check if only one leg evolved so far.
477 IF(NEP.EQ.1) GOTO 330
478 IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 200
480 IF(ITRY(I).EQ.0.AND.KFLD(I).GT.0.AND.(KFLD(I).LE.8.OR.KFLD(I).EQ.
482 IF(P(N+I,5).GE.PMTH(2,KFLD(I))) GOTO 200
486 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
488 PA1S=(P(N+1,4)+P(N+1,5))*(P(N+1,4)-P(N+1,5))
489 PA2S=(P(N+2,4)+P(N+2,5))*(P(N+2,4)-P(N+2,5))
490 PA3S=(P(N+3,4)+P(N+3,5))*(P(N+3,4)-P(N+3,5))
491 PTS=0.25*(2.*PA1S*PA2S+2.*PA1S*PA3S+2.*PA2S*PA3S-
492 & PA1S**2-PA2S**2-PA3S**2)/PA1S
493 IF(PTS.LE.0.) GOTO 200
494 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
497 IF(KFLDA.EQ.0.OR.(KFLDA.GT.8.AND.KFLDA.NE.21)) GOTO 320
498 IF(P(I1,5).LT.PMTH(2,KFLDA)) GOTO 320
507 IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
508 PED=0.5*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
510 IF(I1.EQ.N+1) ZM=V(IM,1)
511 IF(I1.EQ.N+2) ZM=1.-V(IM,1)
512 PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
513 & 4.*V(N+1,5)*V(N+2,5))
514 PED=PEM*(0.5*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/V(IM,5)
516 IF(MOD(MSTJ(43),2).EQ.1) THEN
518 IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
519 PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(I1,5)
520 PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(I1,5)
521 ZD=SQRT(MAX(0.,(1.-V(I1,5)/PED**2)*((1.-PMQ1-PMQ2)**2-
525 ZD=SQRT(MAX(0.,1.-V(I1,5)/PED**2))
530 IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(1)=1
531 IF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(2)=1
532 IF(KFLDA.EQ.21) V(I1,4)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*(1.-ZU)))
533 IF(KFLDA.NE.21) V(I1,4)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
535 IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
538 ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
539 ZDR1=MAX(0.,V(N+1,3)/V(N+1,4)-1.)
540 ZDR2=MAX(0.,V(N+2,3)/V(N+2,4)-1.)
541 IF(ZDR2.GT.RLU(0)*(ZDR1+ZDR2)) ISL(1)=0
542 IF(ISL(1).EQ.1) ISL(2)=0
543 IF(ISL(1).EQ.0) ISLM=1
544 IF(ISL(2).EQ.0) ISLM=2
546 IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 200
548 IF(IGM.GT.0.AND.MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
549 &PMTH(2,KFLD(1)).OR.P(N+2,5).GE.PMTH(2,KFLD(2)))) THEN
550 PMQ1=V(N+1,5)/V(IM,5)
551 PMQ2=V(N+2,5)/V(IM,5)
552 ZD=SQRT(MAX(0.,(1.-V(IM,5)/PEM**2)*((1.-PMQ1-PMQ2)**2-
557 IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 200
560 C...Accepted branch. Construct four-momentum for initial partons.
566 P(N+1,3)=SQRT(MAX(0.,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
570 ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
571 PED1=0.5*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
574 P(N+1,3)=SQRT(MAX(0.,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
579 P(N+2,4)=P(IM,5)-PED1
582 ELSEIF(NEP.EQ.3) THEN
585 P(N+1,3)=SQRT(MAX(0.,PA1S))
588 P(N+2,3)=0.5*(PA3S-PA2S-PA1S)/P(N+1,3)
591 P(N+3,3)=-(P(N+1,3)+P(N+2,3))
596 C...Construct transverse momentum for ordinary branching in shower.
599 PZM=SQRT(MAX(0.,(PEM+P(IM,5))*(PEM-P(IM,5))))
600 PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4.*V(N+1,5)*V(N+2,5)
603 ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
604 PTS=(PEM**2*(ZM*(1.-ZM)*V(IM,5)-(1.-ZM)*V(N+1,5)-
605 & ZM*V(N+2,5))-0.25*PMLS)/PZM**2
607 PTS=PMLS*(ZM*(1.-ZM)*PEM**2/V(IM,5)-0.25)/PZM**2
611 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
613 IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21.
615 IF(K(IGM,3).NE.0) MAZIP=1
617 IF(IAU.EQ.IM+1) ZAU=1.-V(IGM,1)
618 IF(MAZIP.EQ.0) ZAU=0.
619 IF(K(IGM,2).NE.21) THEN
620 HAZIP=2.*ZAU/(1.+ZAU**2)
622 HAZIP=(ZAU/(1.-ZAU*(1.-ZAU)))**2
624 IF(K(N+1,2).NE.21) THEN
625 HAZIP=HAZIP*(-2.*ZM*(1.-ZM))/(1.-2.*ZM*(1.-ZM))
627 HAZIP=HAZIP*(ZM*(1.-ZM)/(1.-ZM*(1.-ZM)))**2
631 C...Find coefficient of azimuthal asymmetry due to soft gluon
634 IF(MSTJ(49).NE.2.AND.MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.
635 & K(N+2,2).EQ.21).AND.IAU.NE.0) THEN
636 IF(K(IGM,3).NE.0) MAZIC=N+1
637 IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
638 IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
639 & ZM.GT.0.5) MAZIC=N+2
640 IF(K(IAU,2).EQ.22) MAZIC=0
642 IF(MAZIC.EQ.N+2) ZS=1.-ZM
644 IF(IAU.EQ.IM-1) ZGM=1.-V(IGM,1)
645 IF(MAZIC.EQ.0) ZGM=1.
646 HAZIC=(P(IM,5)/P(IGM,5))*SQRT((1.-ZS)*(1.-ZGM)/(ZS*ZGM))
647 HAZIC=MIN(0.95,HAZIC)
651 C...Construct kinematics for ordinary branching in shower.
652 340 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
653 IF(MOD(MSTJ(43),2).EQ.1) THEN
656 P(N+1,4)=PEM*(0.5*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
657 & SQRT(PMLS)*ZM)/V(IM,5)
663 P(N+1,3)=0.5*(V(N+2,5)-V(N+1,5)-V(IM,5)+2.*PEM*P(N+1,4))/PZM
669 P(N+2,3)=PZM-P(N+1,3)
670 P(N+2,4)=PEM-P(N+1,4)
671 IF(MSTJ(43).LE.2) THEN
672 V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
673 V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
677 C...Rotate and boost daughters.
679 IF(MSTJ(43).LE.2) THEN
680 BEX=P(IGM,1)/P(IGM,4)
681 BEY=P(IGM,2)/P(IGM,4)
682 BEZ=P(IGM,3)/P(IGM,4)
684 GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1.+GA)-
693 THE=ULANGL(P(IM,3)+GABEP*BEZ,SQRT((P(IM,1)+GABEP*BEX)**2+
694 & (P(IM,2)+GABEP*BEY)**2))
695 PHI=ULANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
697 DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
698 & SIN(THE)*COS(PHI)*P(I,3)
699 DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
700 & SIN(THE)*SIN(PHI)*P(I,3)
701 DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
703 DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
704 DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
705 P(I,1)=DP(1)+DGABP*BEX
706 P(I,2)=DP(2)+DGABP*BEY
707 P(I,3)=DP(3)+DGABP*BEZ
708 350 P(I,4)=GA*(DP(4)+DBP)
711 C...Weight with azimuthal distribution, if required.
712 IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
716 360 DPT(3,J)=P(N+1,J)
717 DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
718 DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
719 DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
721 DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/DPMM
722 370 DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/DPMM
723 DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
724 DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
725 IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN
726 CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
727 & DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
729 IF(1.+HAZIP*(2.*CAD**2-1.).LT.RLU(0)*(1.+ABS(HAZIP)))
733 IF(MAZIC.EQ.N+2) CAD=-CAD
734 IF((1.-HAZIC)*(1.-HAZIC*CAD)/(1.+HAZIC**2-2.*HAZIC*CAD).
735 & LT.RLU(0)) GOTO 340
740 C...Continue loop over partons that may branch, until none left.
741 IF(IGM.GE.0) K(IM,1)=14
744 IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
745 CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
746 IF(MSTU(21).GE.1) N=NS
747 IF(MSTU(21).GE.1) RETURN
751 C...Set information on imagined shower initiator.
752 380 IF(NPA.GE.2) THEN
756 IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
764 C...Reconstruct string drawing information.
766 IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
768 ELSEIF(K(I,1).LE.10) THEN
769 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
770 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
771 ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
772 ID1=MOD(K(I,4),MSTU(5))
773 IF(K(I,2).GE.1.AND.K(I,2).LE.8) ID1=MOD(K(I,4),MSTU(5))+1
774 ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
775 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
776 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
777 K(ID1,4)=K(ID1,4)+MSTU(5)*I
778 K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
779 K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
780 K(ID2,5)=K(ID2,5)+MSTU(5)*I
782 ID1=MOD(K(I,4),MSTU(5))
784 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
785 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
786 K(ID1,4)=K(ID1,4)+MSTU(5)*I
787 K(ID1,5)=K(ID1,5)+MSTU(5)*I
793 C...Transformation from CM frame.
799 GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
800 & /(1.+GA)-P(IPA(1),4))
807 THE=ULANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
808 &+GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
809 PHI=ULANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
811 CHI=ULANGL(COS(THE)*COS(PHI)*(P(IPA(2),1)+GABEP*BEX)+COS(THE)*
812 & SIN(PHI)*(P(IPA(2),2)+GABEP*BEY)-SIN(THE)*(P(IPA(2),3)+GABEP*
813 & BEZ),-SIN(PHI)*(P(IPA(2),1)+GABEP*BEX)+COS(PHI)*(P(IPA(2),2)+
816 CALL LUDBRB(NS+1,N,0.,CHI,0D0,0D0,0D0)
822 CALL LUDBRB(NS+1,N,THE,PHI,DBEX,DBEY,DBEZ)
824 C...Decay vertex of shower.
829 C...Delete trivial shower, else connect initiators.
830 IF(N.EQ.NS+NPA+IIM) THEN
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)