]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ISAJET/code/wpair.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- The part of JETAN dealing with ESD data has been separated from the one using MC...
[u/mrichter/AliRoot.git] / ISAJET / code / wpair.F
CommitLineData
0795afa3 1#include "isajet/pilot.h"
2 SUBROUTINE WPAIR
3C
4C Finish generation of wpair events started bY TWOJET.
5C Select W decay modes as allowed by WMODE1, WMODE2.
6C Generate W decay angles and put vectors in PPAIR.
7C
8C Also generate massless decay vectors PZERO for matrix
9C element -- double precision for 32-bit machines.
10C
11C Ver 6.26: Check kinematics for W -> ff decay, since Z0 from
12C Higgs decay can be virtual.
13C Ver. 6.30: Added check in loop 201.
14C Ver. 7.14: Add MSSM Higgs hooks
15C
16#if defined(CERNLIB_IMPNONE)
17 IMPLICIT NONE
18#endif
19#include "isajet/itapes.inc"
20#include "isajet/qcdpar.inc"
21#include "isajet/jetpar.inc"
22#include "isajet/primar.inc"
23#include "isajet/q1q2.inc"
24#include "isajet/jetsig.inc"
25#include "isajet/wwsig.inc"
26#include "isajet/wwpar.inc"
27#include "isajet/const.inc"
28#include "isajet/qsave.inc"
29#include "isajet/wcon.inc"
30#include "isajet/pjets.inc"
31#include "isajet/pinits.inc"
32#include "isajet/keys.inc"
33#include "isajet/wsig.inc"
34#include "isajet/hcon.inc"
35#include "isajet/xmssm.inc"
36C
37 DIMENSION X(2),LIST(25),P1WCM(4),P2WCM(4),P1LAB(4),P2LAB(4)
38 $,P1CM0(4),P2CM0(4),P1LAB0(4),P2LAB0(4)
39 1,PBOOST(4)
40 EQUIVALENCE (X(1),X1)
41 DIMENSION PWW(5,2)
42 EQUIVALENCE (PWW(1,1),P3WW(1))
43 DIMENSION JWWTYP(2),THWFF(2),PHIWFF(2)
44#if defined(CERNLIB_SINGLE)
45 REAL P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST
46#endif
47#if defined(CERNLIB_DOUBLE)
48 DOUBLE PRECISION P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST
49#endif
50 REAL AMWW1,AMWW2,X,STRUC,STRUCW,RND,RANF,CBRWW,AMASS,AM0,AM1,AM2,
51 $E1CM,E2CM,P12CM,CTHCM,STHCM,PHICM,CPHICM,SPHICM,P1WCM,P2WCM,
52 $PBOOST,P1LAB,P2LAB,AFX,SGWWMX,P1LAB0,P2LAB0,THWFF,PHIWFF
53 INTEGER IH,IQ,JWWTYP,JET,JWT,JQ,IQ1,IQ2,LIST,NREJ,NJ0,K
54 REAL BRANCH(12),SUMBR
55 INTEGER IDABS,IDABS1,IDABS2
56C
57 DATA LIST/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
58 $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16/
59C
60C Initialize for given W type.
61 AMWW1=PJETS(5,1)
62 AMWW2=PJETS(5,2)
63 CALL WWKIN(AMWW1,AMWW2)
64 NPAIR=0
65C
66C Calculate and save structure functions.
67 DO 120 IH=1,2
68 DO 121 IQ=1,13
69121 QSAVE(IQ,IH)=STRUC(X(IH),QSQ,IQ,IDIN(IH))/X(IH)
70 DO 122 IQ=14,25
71122 QSAVE(IQ,IH)=0.
72 IF(KEYS(7).OR.KEYS(9)) THEN
73 DO 123 IQ=26,29
74123 QSAVE(IQ,IH)=STRUCW(X(IH),IQ-25,IDIN(IH))/X(IH)
75 ENDIF
76120 CONTINUE
77C JWWTYP points to W types 1,2,3,4
78 IF(KEYS(6)) THEN
79 JWWTYP(1)=JETTYP(1)
80 JWWTYP(2)=JETTYP(2)
81 ELSEIF((KEYS(7).AND..NOT.GOMSSM).OR.KEYS(9)) THEN
82 JWWTYP(1)=JETTYP(1)-25
83 JWWTYP(2)=JETTYP(2)-25
84 ELSEIF(KEYS(7).AND.GOMSSM) THEN
85 JWWTYP(1)=JETTYP(1)-76
86 JWWTYP(2)=JETTYP(2)-76
87 ENDIF
88C
89C Select W decay modes and put in /JETSET/. First particle
90C is always the fermion.
91
92 DO 200 JET=1,2
93 IDABS=IABS(IDJETS(JET))
94 IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 200
95 RND=RANF()
96 JWT=JWWTYP(JET)
97C Must only consider allowed decays for this mass
98 SUMBR=0.
99 DO 201 JQ=1,12
100 IQ1=2*JQ
101 IQ2=MATCH(IQ1,JWT)
102 IF(IQ2.EQ.0) THEN
103 BRANCH(JQ)=0.
104 GO TO 201
105 ENDIF
106 AM1=AMASS(LIST(IQ1))
107 AM2=AMASS(LIST(IQ2))
108 IF(AM1+AM2.LT.PJETS(5,JET)) THEN
109 BRANCH(JQ)=RBRWW(JQ,JWT,JET)
110 SUMBR=SUMBR+BRANCH(JQ)
111 ELSE
112 BRANCH(JQ)=0.
113 ENDIF
114201 CONTINUE
115 IF(SUMBR.LE.0.) GO TO 998
116 DO 202 JQ=1,12
117202 BRANCH(JQ)=BRANCH(JQ)/SUMBR
118C
119 CBRWW=0.
120 DO 210 JQ=1,12
121 CBRWW=CBRWW+BRANCH(JQ)
122 IF(RND.GT.CBRWW) GO TO 210
123 IQ1=2*JQ
124 IQ2=MATCH(IQ1,JWT)
125 IDPAIR(NPAIR+1)=LIST(IQ1)
126 IDPAIR(NPAIR+2)=LIST(IQ2)
127 PPAIR(5,NPAIR+1)=AMASS(LIST(IQ1))
128 PPAIR(5,NPAIR+2)=AMASS(LIST(IQ2))
129 JPAIR(NPAIR+1)=JET
130 JPAIR(NPAIR+2)=JET
131 NPAIR=NPAIR+2
132 JQWW(JET)=JQ
133 GO TO 200
134210 CONTINUE
135200 CONTINUE
136C
137C Generate decay uniformly in angle and put in PPAIR.
138C Will check cross section later.
139C
140 NREJ=0
141300 NJ0=2
142 DO 310 JET=1,2
143 IDABS=IABS(IDJETS(JET))
144 IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 310
145C Construct W com momenta.
146 AM0=PJETS(5,JET)
147 AM1=PPAIR(5,NJ0-1)
148 AM2=PPAIR(5,NJ0)
149 E1CM=(AM0**2+AM1**2-AM2**2)/(2.*AM0)
150 E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)
151 P12CM=(AM0**2-AM1**2-AM2**2)**2-4.*(AM1*AM2)**2
152 P12CM=SQRT(P12CM)/(2.*AM0)
153 CTHCM=2.*RANF()-1.
154 STHCM=SQRT(1.-CTHCM**2)
155 PHICM=2.*PI*RANF()
156 CPHICM=COS(PHICM)
157 SPHICM=SIN(PHICM)
158 P1WCM(1)=P12CM*STHCM*CPHICM
159 P2WCM(1)=-P1WCM(1)
160 P1WCM(2)=P12CM*STHCM*SPHICM
161 P2WCM(2)=-P1WCM(2)
162 P1WCM(3)=P12CM*CTHCM
163 P2WCM(3)=-P1WCM(3)
164 P1WCM(4)=E1CM
165 P2WCM(4)=E2CM
166C Also construct zero mass vectors at same angle
167#if defined(CERNLIB_SINGLE)
168C Single precision.
169 P1CM0(1)=.5*AM0*STHCM*CPHICM
170 P2CM0(1)=-P1CM0(1)
171 P1CM0(2)=.5*AM0*STHCM*SPHICM
172 P2CM0(2)=-P1CM0(2)
173 P1CM0(3)=.5*AM0*CTHCM
174 P2CM0(3)=-P1CM0(3)
175 P1CM0(4)=.5*AM0
176 P2CM0(4)=P1CM0(4)
177#endif
178#if defined(CERNLIB_DOUBLE)
179C Double precision.
180 DAM0=AM0
181 DCTH=CTHCM
182 DSTH=DSQRT(1.D0-DCTH**2)
183 DPHI=PHICM
184 P1CM0(1)=.5*AM0*DSTH*DCOS(DPHI)
185 P2CM0(1)=-P1CM0(1)
186 P1CM0(2)=.5*AM0*DSTH*DSIN(DPHI)
187 P2CM0(2)=-P1CM0(2)
188 P1CM0(3)=.5*AM0*DCTH
189 P2CM0(3)=-P1CM0(3)
190 P1CM0(4)=.5*AM0
191 P2CM0(4)=P1CM0(4)
192#endif
193C Boost to lab frame.
194 DO 320 K=1,3
195320 PBOOST(K)=-PJETS(K,JET)
196 PBOOST(4)=PJETS(4,JET)
197 CALL LBOOST(PBOOST,1,P1WCM,P1LAB)
198 CALL LBOOST(PBOOST,1,P2WCM,P2LAB)
199 DO 330 K=1,4
200 PPAIR(K,NJ0-1)=P1LAB(K)
201 PPAIR(K,NJ0)=P2LAB(K)
202330 CONTINUE
203C Boost zero mass vectors -- double precision for 32 bits.
204 PZERO(4,NJ0-1)=(P1CM0(4)*PWW(4,JET)+P1CM0(1)*PWW(1,JET)
205 $ +P1CM0(2)*PWW(2,JET)+P1CM0(3)*PWW(3,JET))/PWW(5,JET)
206 BOOST=(P1CM0(4)+PZERO(4,NJ0-1))/(PWW(4,JET)+PWW(5,JET))
207 DO 340 K=1,3
208340 PZERO(K,NJ0-1)=P1CM0(K)+BOOST*PWW(K,JET)
209 PZERO(4,NJ0)=(P2CM0(4)*PWW(4,JET)+P2CM0(1)*PWW(1,JET)
210 $ +P2CM0(2)*PWW(2,JET)+P2CM0(3)*PWW(3,JET))/PWW(5,JET)
211 BOOST=(P2CM0(4)+PZERO(4,NJ0))/(PWW(4,JET)+PWW(5,JET))
212 DO 350 K=1,3
213350 PZERO(K,NJ0)=P2CM0(K)+BOOST*PWW(K,JET)
214 NJ0=NJ0+2
215310 CONTINUE
216C
217C Calculate cross section SIGWW2 containing TBRWW*RBRWW.
218C Compare with WW cross section containing TBRWW. Ratio
219C must be bounded by 3/(4*PI) for each W.
220C
221 AFX=3./(2.*PI)
222 IF(KEYS(6)) THEN
223 CALL SIGWW2
224 SGWWMX=SIGEVT
225 IF(IDJETS(1).NE.10) SGWWMX=SGWWMX*RBRWW(JQWW(1),JWWTYP(1),1)*AFX
226 IF(IDJETS(2).NE.10) SGWWMX=SGWWMX*RBRWW(JQWW(2),JWWTYP(2),2)*AFX
227 ELSEIF(KEYS(7)) THEN
228C Note that except for WW -> WW processes, SIGH3 just computes
229C the decay angular distribution, so it can be used for both
230C for SM and SUSY HL0/HH0 decays; HA0 -> WW is forbidden.
231C For Z + HL0 decays, we just return, ie use phase space.
232 IDABS1=IABS(IDJETS(1))
233 IDABS2=IABS(IDJETS(2))
234 IF(.NOT.(IDABS1.EQ.10.OR.IDABS1.EQ.80.OR.IDABS1.EQ.90)) RETURN
235 IF(.NOT.(IDABS2.EQ.10.OR.IDABS2.EQ.80.OR.IDABS2.EQ.90)) RETURN
236 CALL SIGH3
237 SGWWMX=SIGLLQ*AFX**2
238 ELSEIF(KEYS(9)) THEN
239 CALL SIGTC3
240 SGWWMX=SIGLLQ*AFX**2
241 ENDIF
242 IF(WWSIG.GT.SGWWMX*RANF()) GO TO 400
243 NREJ=NREJ+1
244 IF(NREJ.LT.NTRIES) GO TO 300
245 GO TO 999
246C
247C Good event
248C
249400 CONTINUE
250 RETURN
251C
252999 CALL PRTEVT(0)
253 WRITE(ITLIS,9991) NREJ
2549991 FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE ',
255 1'A GOOD WPAIR EVENT.'/' CHECK LIMITS OR INCREASE NTRIES.')
256 STOP 99
257998 CALL PRTEVT(0)
258 WRITE(ITLIS,9981) JET
2599981 FORMAT(//' ERROR IN WPAIR ... NO DECAY POSSIBLE FOR JET',I3)
260 STOP 99
261 END