1 #include "isajet/pilot.h"
3 C----------------------------------------------------------------------
5 C- Purpose and Methods :
6 C- decay W's and Z's as done in ISAJET
8 C- Created 6-MAY-1991 Serban D. Protopopescu
10 C----------------------------------------------------------------------
11 #if defined(CERNLIB_IMPNONE)
14 #include "isajet/const.inc"
15 #include "isajet/frame.inc"
16 #include "isajet/jetpar.inc"
17 #include "isajet/jetset.inc"
18 #include "isajet/jwork.inc"
19 #include "isajet/pjets.inc"
20 #include "isajet/partcl.inc"
21 #include "isajet/primar.inc"
22 #include "isajet/wcon.inc"
25 REAL PREST(5),PL(5),EL(3),EML(3),EMSQL(3)
28 INTEGER LISTJ(29),LISTW(4)
29 REAL RANF,SUM,PTDEN,QDEN,ETA,
30 $S12,SUMBR,BRMODE,AMASS,BRINV,TRY,PL12,
31 $COSTHL,THL,PHL,PTL,SGN,BP,PLPL,PLMN,AMINI,AMFIN,PINI,PFIN,
32 $ QPL,QMN,AM1SQ,AM2SQ,ROOT,P1PL,P1MN,P2PL,P2MN
33 INTEGER NADD,K,IQ1,IQ2,IFL1,IFL2,IQ,IFL,I
38 $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
39 $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,
41 DATA LISTW/10,80,-80,90/
42 C----------------------------------------------------------------------
48 C Kinematics. Note that YW is the true rapidity and QW is
49 C the true 3-momentum. See DRLLYN.
52 QTW=SQRT(QWJET(1)**2+QWJET(2)**2)
53 QW=SQRT(QWJET(1)**2+QWJET(2)**2+QWJET(3)**2)
55 PHIW=ATAN2(QWJET(2),QWJET(1))
56 IF(PHIW.LT.0) PHIW=PHIW+2*PI
60 QWPL=QWJET(4)+QWJET(3)
61 QWMN=QWJET(4)-QWJET(3)
62 IF(QWPL.GT.0..AND.QWMN.GT.0.) THEN
63 YW=0.5*ALOG(QWPL/QWMN)
65 YW=999.*SIGN(1.,QWJET(3))
74 C QMW dependence neglected in branching ratios
75 C BRANCH is cum. br. with heavy modes subtracted.
85 BRMODE=WCBR(IQ1,JWTYP)-WCBR(IQ1-1,JWTYP)
88 IF(S12.LE.(AMASS(IFL1)+AMASS(IFL2))**2) BRMODE=0.
90 BRANCH(IQ1)=BRANCH(IQ1-1)+BRMODE
97 IF(TRY.LT.BRANCH(IQ)*BRINV.AND.MATCH(IQ,JWTYP).NE.0) THEN
99 JETTYP(2)=MATCH(IQ,JWTYP)
104 120 IFL1=LISTJ(JETTYP(1))
105 IFL2=LISTJ(JETTYP(2))
107 C Select masses of decay products.
112 C Generate W decay in its rest frame
113 C First set up momenta of decay products:
117 EL(1)=(S12+EMSQL(1)-EMSQL(2))/(2.*QMW)
118 EL(2)=(S12+EMSQL(2)-EMSQL(1))/(2.*QMW)
119 PL12=SQRT((S12-(EML(1)+EML(2))**2)*(S12-(EML(1)-EML(2))**2))
123 140 PREST(K)=QWJET(K)
124 C Generate next W decay
133 PL(1)=SGN*PTL*COS(PHL)
134 PL(2)=SGN*PTL*SIN(PHL)
135 PL(3)=SGN*PL12*COSTHL
138 C Boost with W momentum
141 310 BP=BP+PL(K)*PREST(K)
144 320 PL(K)=PL(K)+PREST(K)*PL(4)/PREST(5)
145 $ +PREST(K)*BP/(PREST(4)+PREST(5))
146 PL(4)=PL(4)*PREST(4)/PREST(5)+BP
148 PT(I)=SQRT(PL(1)**2+PL(2)**2)
149 P(I)=SQRT(PT(I)**2+PL(3)**2)
151 PHI(I)=ATAN2(PL(2),PL(1))
155 IF(PHI(I).LT.0.) PHI(I)=PHI(I)+2.*PI
160 IF(CTH(I).GT.0.) THEN
162 PLMN=(PT(I)**2+EMSQL(I))/PLPL
165 PLPL=(PT(I)**2+EMSQL(I))/PLMN
167 YJ(I)=.5*ALOG(PLPL/PLMN)
173 PJETS(3,I)=P(I)*CTH(I)
174 PJETS(1,I)=PT(I)*COS(PHI(I))
175 PJETS(2,I)=PT(I)*SIN(PHI(I))
176 PJETS(4,I)=SQRT(P(I)**2+EMSQL(I))
177 PJETS(5,I)=SQRT(EMSQL(I))
178 IDJETS(I)=LISTJ(JETTYP(I))