2 C*********************************************************************
6 C...Purpose: to handle the decay of unstable particles.
7 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
8 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
9 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
10 COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
11 SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
12 DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),
13 &WTCOR(10),PTAU(4),PCMTAU(4)
14 DOUBLE PRECISION DBETAU(3)
15 DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./
17 C...Functions: momentum in two-particle decays, four-product and
18 C...matrix element times phase space in weak decays.
19 PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)
20 FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3)
21 HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))*
22 &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA)
32 C...Choose lifetime and determine decay vertex.
35 ELSEIF(K(IP,1).NE.4) THEN
36 V(IP,5)=-PMAS(KC,4)*LOG(RLU(0))
39 VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
42 C...Determine whether decay allowed or not.
44 IF(MSTJ(22).EQ.2) THEN
45 IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1
46 ELSEIF(MSTJ(22).EQ.3) THEN
47 IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
48 ELSEIF(MSTJ(22).EQ.4) THEN
49 IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
50 IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
52 IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN
57 C...Interface to external tau decay library (for tau polarization).
58 IF(KFA.EQ.15.AND.MSTJ(28).GE.1) THEN
60 C...Starting values for pointers and momenta.
67 C...Iterate to find position and code of mother of tau.
72 C...If no known origin then impossible to do anything further.
76 ELSEIF(K(IMTAU,2).EQ.K(ITAU,2)) THEN
77 C...If tau -> tau + gamma then add gamma energy and loop.
78 IF(K(K(IMTAU,4),2).EQ.22) THEN
80 PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,4),J)
82 ELSEIF(K(K(IMTAU,5),2).EQ.22) THEN
84 PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,5),J)
89 ELSEIF(IABS(K(IMTAU,2)).GT.100) THEN
90 C...If coming from weak decay of hadron then W is not stored in record,
91 C...but can be reconstructed by adding neutrino momentum.
92 KFORIG=-ISIGN(24,K(ITAU,2))
94 DO 160 II=K(IMTAU,4),K(IMTAU,5)
95 IF(K(II,2)*ISIGN(1,K(ITAU,2)).EQ.-16) THEN
97 PCMTAU(J)=PCMTAU(J)+P(II,J)
103 C...If coming from resonance decay then find latest copy of this
104 C...resonance (may not completely agree).
107 DO 170 II=IMTAU+1,IP-1
108 IF(K(II,2).EQ.KFORIG.AND.K(II,3).EQ.IORIG.AND.
109 & ABS(P(II,5)-P(IORIG,5)).LT.1E-5*P(IORIG,5)) IORIG=II
116 C...Boost tau to rest frame of production process (where known)
117 C...and rotate it to sit along +z axis.
119 DBETAU(J)=PCMTAU(J)/PCMTAU(4)
121 IF(KFORIG.NE.0) CALL LUDBRB(ITAU,ITAU,0.,0.,-DBETAU(1),
122 & -DBETAU(2),-DBETAU(3))
123 PHITAU=ULANGL(P(ITAU,1),P(ITAU,2))
124 CALL LUDBRB(ITAU,ITAU,0.,-PHITAU,0D0,0D0,0D0)
125 THETAU=ULANGL(P(ITAU,3),P(ITAU,1))
126 CALL LUDBRB(ITAU,ITAU,-THETAU,0.,0D0,0D0,0D0)
128 C...Call tau decay routine (if meaningful) and fill extra info.
129 IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
130 CALL LUTAUD(ITAU,IORIG,KFORIG,NDECAY)
131 DO 200 II=NSAV+1,NSAV+NDECAY
140 C...Boost back decay tau and decay products.
144 IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
145 CALL LUDBRB(NSAV+1,N,THETAU,PHITAU,0D0,0D0,0D0)
146 IF(KFORIG.NE.0) CALL LUDBRB(NSAV+1,N,0.,0.,DBETAU(1),
147 & DBETAU(2),DBETAU(3))
149 C...Skip past ordinary tau decay treatment.
157 C...B-B~ mixing: flip sign of meson appropriately.
159 IF((KFA.EQ.511.OR.KFA.EQ.531).AND.MSTJ(26).GE.1) THEN
161 IF(KFA.EQ.531) XBBMIX=PARJ(77)
162 IF(SIN(0.5*XBBMIX*V(IP,5)/PMAS(KC,4))**2.GT.RLU(0)) MMIX=1
163 IF(MMIX.EQ.1) KFS=-KFS
166 C...Check existence of decay channels. Particle/antiparticle rules.
168 IF(MDCY(KC,2).GT.0) THEN
169 MDMDCY=MDME(MDCY(KC,2),2)
170 IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY
172 IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN
173 CALL LUERRM(9,'(LUDECY:) no decay channel defined')
176 IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS
177 IF(KCHG(KC,3).EQ.0) THEN
180 IF(RLU(0).GT.0.5) KFS=-KFS
181 ELSEIF(KFS.GT.0) THEN
189 C...Sum branching ratios of allowed decay channels.
192 DO 230 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1
193 IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
194 &KFSN*MDME(IDL,1).NE.3) GOTO 230
195 IF(MDME(IDL,2).GT.100) GOTO 230
200 CALL LUERRM(2,'(LUDECY:) all decay channels closed by user')
204 C...Select decay channel among allowed ones.
208 IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
209 &KFSN*MDME(IDL,1).NE.3) THEN
210 IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
211 ELSEIF(MDME(IDL,2).GT.100) THEN
212 IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
216 IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 250
219 C...Start readout of decay channel: matrix element, reset counters.
222 IF(NTRY.GT.1000) THEN
223 CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
224 IF(MSTU(21).GE.1) RETURN
230 IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1
233 IF(MBST.EQ.0) PV(1,J)=P(IP,J)
235 IF(MBST.EQ.1) PV(1,4)=P(IP,5)
241 IF(KFA.GT.80) MHADDY=1
243 C...Read out decay products. Convert to standard flavour code.
245 IF(MDME(IDC+1,2).EQ.101) JTMAX=10
247 IF(JT.LE.5) KP=KFDP(IDC,JT)
248 IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)
252 IF(KPA.GT.80) MHADDY=1
253 IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN
255 ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN
257 ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN
258 KFP=-KFS*MOD(KFA/10,10)
259 ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN
260 KFP=KFS*(100*MOD(KFA/10,100)+3)
261 ELSEIF(KPA.EQ.81) THEN
262 KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1)
263 ELSEIF(KP.EQ.82) THEN
264 CALL LUKFDI(-KFS*INT(1.+(2.+PARJ(2))*RLU(0)),0,KFP,KDUMP)
265 IF(KFP.EQ.0) GOTO 260
267 IF(PV(1,5).LT.PARJ(32)+2.*ULMASS(KFP)) GOTO 260
268 ELSEIF(KP.EQ.-82) THEN
270 IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP)
272 IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP(KFP)
274 C...Add decay product to event record or to quark flavour list.
277 IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN
281 PSQ=PSQ+ULMASS(KFLO(NQ))
282 ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.48).AND.NP.EQ.3.AND.
283 &MOD(NQ,2).EQ.1) THEN
288 CALL LUKFDI(KFP,KFI,KFLDMP,K(I,2))
289 IF(K(I,2).EQ.0) GOTO 260
291 P(I,5)=ULMASS(K(I,2))
296 IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1
297 IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1
299 IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2
300 IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1
306 IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32)
311 C...Check masses for resonance decays.
313 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 240
316 C...Choose decay multiplicity in phase space model.
317 290 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN
319 CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1))
320 IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)
322 IF(NTRY.GT.1000) THEN
323 CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
324 IF(MSTU(21).GE.1) RETURN
327 GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU(0))))*
328 & SIN(PARU(2)*RLU(0))
329 ND=0.5+0.5*NP+0.25*NQ+CNDE+GAUSS
330 IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 300
331 IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 300
332 IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 300
333 IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 300
338 C...Form hadrons from flavour content.
342 IF(ND.EQ.NP+NQ/2) GOTO 330
343 DO 320 I=N+NP+1,N+ND-NQ/2
344 JT=1+INT((NQ-1)*RLU(0))
345 CALL LUKFDI(KFL1(JT),0,KFL2,K(I,2))
346 IF(K(I,2).EQ.0) GOTO 300
352 IF(NQ.EQ.4.AND.RLU(0).LT.PARJ(66)) JT=4
353 IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))*
354 & ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3
357 CALL LUKFDI(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))
358 IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 300
359 IF(NQ.EQ.4) CALL LUKFDI(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))
360 IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 300
362 C...Check that sum of decay product masses not too large.
369 P(I,5)=ULMASS(K(I,2))
372 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 300
374 C...Rescale energy to subtract off spectator quark mass.
375 ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45)
378 PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)
380 P(N+NP,J)=PQT*PV(1,J)
381 PV(1,J)=(1.-PQT)*PV(1,J)
383 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
387 C...Phase space factors imposed in W decay.
388 ELSEIF(MMAT.EQ.46) THEN
390 PSMC=ULMASS(K(N+1,2))
392 PSMC=PSMC+ULMASS(K(N+2,2))
393 IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 240
394 HR1=(P(N+1,5)/PV(1,5))**2
395 HR2=(P(N+2,5)/PV(1,5))**2
396 IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2)
397 & .LT.2.*RLU(0)) GOTO 240
400 C...Fully specified final state: check mass broadening effects.
402 IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 260
406 C...Select W mass in decay Q -> W + q, without W propagator.
407 IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN
408 HLQ=(PARJ(32)/PV(1,5))**2
409 HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2
410 HRQ=(P(N+2,5)/PV(1,5))**2
411 360 HW=HLQ+RLU(0)*(HUQ-HLQ)
412 IF(HMEPS(HW).LT.RLU(0)) GOTO 360
413 P(N+1,5)=PV(1,5)*SQRT(HW)
415 C...Ditto, including W propagator. Divide mass range into three regions.
416 ELSEIF(MMAT.EQ.45) THEN
417 HQW=(PV(1,5)/PMAS(24,1))**2
418 HLW=(PARJ(32)/PMAS(24,1))**2
419 HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2
420 HRQ=(P(N+2,5)/PV(1,5))**2
421 HG=PMAS(24,2)/PMAS(24,1)
422 HATL=ATAN((HLW-1.)/HG)
424 HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
426 HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
427 IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN
431 HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2)
432 HM1=1.-SQRT(1./HMV-HG**2)
433 IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN
435 ELSEIF(HMV2.LE.HMV1) THEN
436 HM=MAX(HLW,HM-MIN(0.1,1.-HM))
438 HATM=ATAN((HM-1.)/HG)
440 HWT2=HMV*(MIN(1.,HUW)-HM)
443 HATU=ATAN((HUW-1.)/HG)
448 C...Select mass region and W mass there. Accept according to weight.
449 380 HREG=RLU(0)*(HWT1+HWT2+HWT3)
450 IF(HREG.LE.HWT1) THEN
451 HW=1.+HG*TAN(HATL+RLU(0)*(HATM-HATL))
453 ELSEIF(HREG.LE.HWT1+HWT2) THEN
454 HW=HM+RLU(0)*(MIN(1.,HUW)-HM)
455 HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV
457 HW=1.+HG*TAN(RLU(0)*HATU)
458 HACC=HMEPS(HW/HQW)/HMP1
460 IF(HACC.LT.RLU(0)) GOTO 380
461 P(N+1,5)=PMAS(24,1)*SQRT(HW)
464 C...Determine position of grandmother, number of sisters, Q -> W sign.
468 IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN
470 IF(IM.LT.0.OR.IM.GE.IP) IM=0
471 IF(MMAT.EQ.46.AND.MSTJ(27).EQ.1) THEN
473 ELSEIF(MMAT.EQ.46.AND.MSTJ(27).GE.2.AND.IM.NE.0) THEN
474 IF(K(IM,2).EQ.94) THEN
476 IF(IM.LT.0.OR.IM.GE.IP) IM=0
479 IF(IM.NE.0) KFAM=IABS(K(IM,2))
480 IF(IM.NE.0.AND.MMAT.EQ.3) THEN
481 DO 390 IL=MAX(IP-2,IM+1),MIN(IP+2,N)
482 IF(K(IL,3).EQ.IM) NM=NM+1
483 IF(K(IL,3).EQ.IM.AND.IL.NE.IP) ISIS=IL
485 IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.
486 & MOD(KFAM/1000,10).NE.0) NM=0
489 IF((KFAS.LE.100.OR.MOD(KFAS,10).NE.1.OR.
490 & MOD(KFAS/1000,10).NE.0).AND.KFAS.NE.22) NM=0
492 ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN
493 MSGN=ISIGN(1,K(IM,2)*K(IP,2))
494 IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN=
495 & MSGN*(-1)**MOD(KFAM/100,10)
499 C...Kinematics of one-particle decays.
507 C...Calculate maximum weight ND-particle decay.
511 PMAX=PV(1,5)-PS+P(N+ND,5)
515 PMIN=PMIN+P(N+IL+1,5)
516 WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))
520 C...Find virtual gamma mass in Dalitz decay.
522 ELSEIF(MMAT.EQ.2) THEN
523 PMES=4.*PMAS(11,1)**2
524 PMRHO2=PMAS(131,1)**2
525 PGRHO2=PMAS(131,2)**2
526 430 PMST=PMES*(P(IP,5)**2/PMES)**RLU(0)
527 WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))*
528 & (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/
529 & ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2)
530 IF(WT.LT.RLU(0)) GOTO 430
531 PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST))
533 C...M-generator gives weight. If rejected, try again.
538 DO 450 IL2=IL1-1,1,-1
539 IF(RSAV.LE.RORD(IL2)) GOTO 460
540 RORD(IL2+1)=RORD(IL2)
547 PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
548 WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
550 IF(WT.LT.RLU(0)*WTMAX) GOTO 440
553 C...Perform two-particle decays in respective CM frame.
555 PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
558 UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
559 UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
564 P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)
565 PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
568 C...Lorentz transform decay products to lab frame.
574 BE(J)=PV(IL,J)/PV(IL,4)
578 BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
580 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
582 P(I,4)=GA*(P(I,4)+BEP)
586 C...Check that no infinite loop in matrix element weight.
588 IF(NTRY.GT.800) GOTO 590
590 C...Matrix elements for omega and phi decays.
592 WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2
593 & -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2
594 & +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)
595 IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU(0)) GOTO 420
597 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
598 ELSEIF(MMAT.EQ.2) THEN
601 WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+
602 & PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)
603 IF(WT.LT.RLU(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 490
605 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
606 C...V vector), of form cos**2(theta02) in V1 rest frame, and for
607 C...S0 -> gamma + V1 -> gamma + S2 + S3, of form sin**2(theta02).
608 ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN
615 IF(KFAS.NE.22) HNUM=(FOUR10*FOUR12-PMS1*FOUR02)**2
616 IF(KFAS.EQ.22) HNUM=PMS1*(2.*FOUR10*FOUR12*FOUR02-
617 & PMS1*FOUR02**2-PMS0*FOUR12**2-PMS2*FOUR10**2+PMS1*PMS0*PMS2)
618 HNUM=MAX(1E-6*PMS1**2*PMS0*PMS2,HNUM)
619 HDEN=(FOUR10**2-PMS1*PMS0)*(FOUR12**2-PMS1*PMS2)
620 IF(HNUM.LT.RLU(0)*HDEN) GOTO 490
622 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
623 ELSEIF(MMAT.EQ.4) THEN
624 HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
625 HX2=2.*FOUR(IP,N+2)/P(IP,5)**2
626 HX3=2.*FOUR(IP,N+3)/P(IP,5)**2
627 WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+
628 & ((1.-HX3)/(HX1*HX2))**2
629 IF(WT.LT.2.*RLU(0)) GOTO 420
630 IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2)
633 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
634 ELSEIF(MMAT.EQ.41) THEN
635 HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
636 HXM=MIN(0.75,2.*(1.-PS/P(IP,5)))
637 IF(HX1*(3.-2.*HX1).LT.RLU(0)*HXM*(3.-2.*HXM)) GOTO 420
639 C...Matrix elements for weak decays (only semileptonic for c and b)
640 ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
642 IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3)
643 IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3)
644 IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 420
645 ELSEIF(MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48) THEN
649 P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J)
652 IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)
653 IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)
654 IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 420
656 C...Angular distribution in W decay.
657 ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN
658 IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1)
659 IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1)
660 IF(WT.LT.RLU(0)*P(IM,5)**4/WTCOR(10)) GOTO 490
663 C...Scale back energy and reattach spectator.
664 590 IF(MREM.EQ.1) THEN
666 PV(1,J)=PV(1,J)/(1.-PQT)
672 C...Low invariant mass for system with spectator quark gives particle,
673 C...not two jets. Readjust momenta accordingly.
674 IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN
679 IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE.
680 & (PARJ(32)+PM2+PM3)**2) GOTO 660
683 CALL LUKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))
684 IF(K(N+2,2).EQ.0) GOTO 260
685 P(N+2,5)=ULMASS(K(N+2,2))
691 ELSEIF(MMAT.EQ.44) THEN
696 IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE.
697 & (PARJ(32)+PM3+PM4)**2) GOTO 630
700 CALL LUKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))
701 IF(K(N+3,2).EQ.0) GOTO 260
702 P(N+3,5)=ULMASS(K(N+3,2))
704 P(N+3,J)=P(N+3,J)+P(N+4,J)
706 P(N+3,4)=SQRT(P(N+3,1)**2+P(N+3,2)**2+P(N+3,3)**2+P(N+3,5)**2)
707 HA=P(N+1,4)**2-P(N+2,4)**2
708 HB=HA-(P(N+1,5)**2-P(N+2,5)**2)
709 HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+
710 & (P(N+1,3)-P(N+2,3))**2
711 HD=(PV(1,4)-P(N+3,4))**2
712 HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2
715 HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF)
717 PCOR=HH*(P(N+1,J)-P(N+2,J))
718 P(N+1,J)=P(N+1,J)+PCOR
719 P(N+2,J)=P(N+2,J)-PCOR
721 P(N+1,4)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2+P(N+1,5)**2)
722 P(N+2,4)=SQRT(P(N+2,1)**2+P(N+2,2)**2+P(N+2,3)**2+P(N+2,5)**2)
726 C...Check invariant mass of W jets. May give one particle or start over.
727 630 IF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
728 &.AND.IABS(K(N+1,2)).LT.10) THEN
729 PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2)))
734 IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 640
735 KFLDUM=INT(1.5+RLU(0))
736 CALL LUKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)
737 CALL LUKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)
738 IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 260
739 PSM=ULMASS(KF1)+ULMASS(KF2)
740 IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.PMR.GT.PARJ(64)+PSM) GOTO 640
741 IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 640
742 IF(MMAT.EQ.48) GOTO 420
743 IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 260
746 CALL LUKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))
747 IF(K(N+1,2).EQ.0) GOTO 260
748 P(N+1,5)=ULMASS(K(N+1,2))
752 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
759 C...Phase space decay of partons from W decay.
760 640 IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.IABS(K(N+1,2)).LT.10) THEN
766 PV(1,J)=P(N+1,J)+P(N+2,J)
777 PSQ=PSQ+ULMASS(KFLO(2))
782 C...Boost back for rapidly moving particle.
786 BE(J)=P(IP,J)/P(IP,4)
790 BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
792 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
794 P(I,4)=GA*(P(I,4)+BEP)
798 C...Fill in position of decay vertex.
806 C...Set up for parton shower evolution from jets.
807 IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN
811 K(NSAV+1,4)=MSTU(5)*(NSAV+2)
812 K(NSAV+1,5)=MSTU(5)*(NSAV+3)
813 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
814 K(NSAV+2,5)=MSTU(5)*(NSAV+1)
815 K(NSAV+3,4)=MSTU(5)*(NSAV+1)
816 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
818 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN
821 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
822 K(NSAV+2,5)=MSTU(5)*(NSAV+3)
823 K(NSAV+3,4)=MSTU(5)*(NSAV+2)
824 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
826 ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46)
827 &.AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN
830 K(NSAV+1,4)=MSTU(5)*(NSAV+2)
831 K(NSAV+1,5)=MSTU(5)*(NSAV+2)
832 K(NSAV+2,4)=MSTU(5)*(NSAV+1)
833 K(NSAV+2,5)=MSTU(5)*(NSAV+1)
835 ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46)
836 &.AND.IABS(K(NSAV+1,2)).LE.20.AND.IABS(K(NSAV+2,2)).LE.20) THEN
838 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)
843 KCP=LUCOMP(K(NSAV+1,2))
844 KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))
847 K(NSAV+1,JCON)=MSTU(5)*(NSAV+2)
848 K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)
849 K(NSAV+2,JCON)=MSTU(5)*(NSAV+3)
850 K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)
852 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN
855 K(NSAV+1,4)=MSTU(5)*(NSAV+3)
856 K(NSAV+1,5)=MSTU(5)*(NSAV+3)
857 K(NSAV+3,4)=MSTU(5)*(NSAV+1)
858 K(NSAV+3,5)=MSTU(5)*(NSAV+1)
861 C...Set up for parton shower evolution in t -> W + b.
862 ELSEIF(MSTJ(27).GE.1.AND.MMAT.EQ.45.AND.ND.EQ.3) THEN
865 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
866 K(NSAV+2,5)=MSTU(5)*(NSAV+3)
867 K(NSAV+3,4)=MSTU(5)*(NSAV+2)
868 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
872 C...Mark decayed particle; special option for B-B~ mixing.
873 IF(K(IP,1).EQ.5) K(IP,1)=15
874 IF(K(IP,1).LE.10) K(IP,1)=11
875 IF(MMIX.EQ.1.AND.MSTJ(26).EQ.2.AND.K(IP,1).EQ.11) K(IP,1)=12