More volume overlaps corrected
[u/mrichter/AliRoot.git] / ISAJET / code / drllyn.F
1 #include "isajet/pilot.h"
2       SUBROUTINE DRLLYN
3 C
4 C          Generate QMW (and QTW) for DRELLYAN or HIGGS event using
5 C          integrated cross section. Then generate decay -- for HIGGS,
6 C          the mode must be chosen using the integrated cross sections
7 C          because of interference with W+W->W+W scattering.
8 C
9 C          Note that NOGOOD calls the cross section.
10 C
11 C          Ver. 6.40: Add technicolor resonances. Use logs for QDEN,
12 C          PTDEN, WTFAC, etc. Also scale QMW generation by QMAX.
13 C
14 C          Ver. 7.01: Correct QDEN to correspond to correct fit form:
15 C          SIGMA = ANOMR(K)*(QMAX**2/Q**2)**QPOW
16 C          See QFUNC.
17 C
18 C          Ver. 7.14: Add SUSY Higgs
19 C          Ver. 7.15: Fix bug with THETAW limits by adding epsilon to
20 C          allowed range. Check for possible invalid Higgs decays.
21 C
22 #if defined(CERNLIB_IMPNONE)
23       IMPLICIT NONE
24 #endif
25 #include "isajet/itapes.inc"
26 #include "isajet/jetsig.inc"
27 #include "isajet/totals.inc"
28 #include "isajet/q1q2.inc"
29 #include "isajet/partcl.inc"
30 #include "isajet/pjets.inc"
31 #include "isajet/pinits.inc"
32 #include "isajet/wcon.inc"
33 #include "isajet/primar.inc"
34 #include "isajet/dylim.inc"
35 #include "isajet/const.inc"
36 #include "isajet/jetpar.inc"
37 #include "isajet/jetlim.inc"
38 #include "isajet/wgen.inc"
39 #include "isajet/dypar.inc"
40 #include "isajet/keys.inc"
41 #include "isajet/hcon.inc"
42 #include "isajet/isloop.inc"
43 #include "isajet/idrun.inc"
44 #include "isajet/xmssm.inc"
45 #include "isajet/listss.inc"
46 C
47       DIMENSION X(2)
48       EQUIVALENCE (X(1),X1)
49       DIMENSION PREST(5),PL(5),EL(3),EML(3),EMSQL(3)
50       DIMENSION WTFAC(3)
51       LOGICAL NOGOOD
52       LOGICAL YGENJ
53       DIMENSION BRANCH(29),LISTJ(29),LISTW(5)
54       REAL ACOSH,XXX,ASINH,CHOOSE,RANF,SUM,WTFAC,PTDEN,QDEN,ETA,QPW,
55      $S12,BRANCH,SUMBR,BRMODE,AMASS,BRINV,TRY,EMSQL,EL,PL12,PREST,
56      $COSTHL,THL,PHL,PTL,SGN,PL,BP,PLPL,PLMN,AMINI,AMFIN,PINI,PFIN,
57      $ QPL,QMN,AM1SQ,AM2SQ,ROOT,P1PL,P1MN,P2PL,P2MN,X,EML
58       INTEGER NTRY,K,IQ1,IQ2,IFL1,IFL2,LISTJ,IQ,NTRY2,IFL,LISTW,I
59       REAL ZZSTAR
60       INTEGER IZSTAR,JVIR,N0J
61 C
62       DATA LISTJ/
63      $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
64      $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,
65      $10,80,-80,90/
66       DATA LISTW/10,80,-80,90,92/
67       ACOSH(XXX)=ALOG(XXX+SQRT(XXX**2-1.))
68       ASINH(XXX)=ALOG(XXX+SQRT(XXX**2+1.))
69 C
70 C          Entry
71 C
72       NPTCL=0
73       NTRY=0
74 200   CONTINUE
75       SIGMA=0.
76       WT=1.
77     1 CONTINUE
78       NTRY=NTRY+1
79       IF(NTRY.GT.NTRIES) GO TO 999
80       SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
81       NKINPT=NKINPT+1
82       SIGMA=0.
83       WT=1.
84       DO 2 K=1,3
85     2 SIGSL(K)=0
86 C            Choose interval for cross section calculation
87       CHOOSE=RANF()
88       SUM=0.
89       DO 3 K=NKL,NKH
90         SUM=SUM+QSELWT(K)
91         IF(CHOOSE.LE.SUM) GO TO 30
92 3     CONTINUE
93 30    KSEL=K
94 C          Generate QTW in selected region
95       IF(.NOT.FIXQT) THEN
96         ETA=(PTGN(1,K)+PTGN(2,K)*RANF())**PTGN(3,K)
97         PTSEL(K)=SQRT(ETA-RNU2(K))
98         PTDEN=ALOG(ETA)*PTPOW(K)
99         WTFAC(1)=ALOG(ABS(PTGN(2,K)))+ALOG(ABS(PTGN(3,K)))
100      1  +ALOG(ABS(PTSEL(K)**2+RNU2(K)))*((PTGN(3,K)-1.)/PTGN(3,K))
101         PT(3)=PTSEL(K)
102       ELSE
103         PTDEN=0.
104         WTFAC(1)=-1000.
105       ENDIF
106 C          Generate QMW
107       IF(.NOT.FIXQM) THEN
108         IF(.NOT.K.EQ.2) THEN
109           QSEL(K)=QMAX**2*(QGEN(1,K)+QGEN(2,K)*RANF())**QGEN(3,K)
110           QDEN=ALOG(QSEL(K)/QMAX**2)*QPOW(K)
111           WTFAC(2)=ALOG(ABS(QGEN(2,K)))+ALOG(ABS(QGEN(3,K)))
112      1    +ALOG(QSEL(K)/QMAX**2)*((QGEN(3,K)-1.)/QGEN(3,K))
113      2    +ALOG(QMAX**2)
114           QSEL(K)=SQRT(QSEL(K))
115           QMW=QSEL(K)
116         ELSE
117           ETA=QGEN(3,K)*TAN(QGEN(1,K)+QGEN(2,K)*RANF())
118           QSEL(K)=SQRT(ETA+EMSQ)
119           WTFAC(2)=ALOG(QGEN(2,K))+ALOG(QGEN(3,K))
120      1    +ALOG((ETA/QGEN(3,K))**2+1.)
121           QMW=QSEL(K)
122           QDEN=ALOG((QMW**2-EMSQ)**2+EMGAM**2)
123         ENDIF
124       ELSE
125         QDEN=0.
126         WTFAC(2)=-1000.
127       ENDIF
128       SIGSL(K)=EXP(ANORM(K)-PTDEN-QDEN)
129 C
130       IF(STDDY) THEN
131         WT=EXP(WTFAC(2)-ALOG(QSELWT(K)))
132       ELSE
133         WT=EXP(WTFAC(1)+WTFAC(2)-ALOG(QSELWT(K)))
134       ENDIF      
135       QTW=PT(3)
136       YW=YWMIN+(YWMAX-YWMIN)*RANF()
137       WT=WT*(YWMAX-YWMIN)
138       PHIW=PHWMIN+(PHWMAX-PHWMIN)*RANF()
139       PHI(3)=AMOD(PHIW+PI,2.*PI)
140       QPW=SQRT(QTW**2+QMW**2)*SINH(YW)
141       QW=SQRT(QTW**2+QPW**2)
142       THW=QPW/QW
143       IF(ABS(THW).GT.1.) THW=SIGN(1.,THW)
144       THW=ACOS(THW)
145       IF(THW.LT.THWMIN-1.E-6.OR.THW.GT.THWMAX+1.E-6) GOTO 1
146       XW=QPW/HALFE
147       IF(XW.LT.XWMIN.OR.XW.GT.XWMAX) GOTO 1
148       IF(.NOT.STDDY) THEN
149         IF(.NOT.YGENJ(3)) GOTO 1
150         P(3)=PT(3)/STH(3)
151         XJ(3)=P(3)*CTH(3)/HALFE
152         IF(XJ(3).LT.XJMIN(3).OR.XJ(3).GT.XJMAX(3)) GOTO 1
153       ENDIF
154 C
155 C          Check integrated cross section
156 C
157       IF(NOGOOD(2)) GO TO 1
158       SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
159       NWGEN=NWGEN+1
160       S12=QMW**2
161 C
162 C          No decay for KKG:
163 C          For compatibility reasons, the jet is still the 3rd one.
164 C          Jets 1 and 2 (W decay products) are voided; no decay step.
165 C
166       IF(KEYS(11)) THEN
167         DO 50 I=1,2
168           P(I)=0.
169           PT(I)=0.
170           CTH(I)=0.
171           PHI(I)=0.
172           EMSQL(I)=0.
173           IDJETS(I)=0
174 50      CONTINUE
175         GOTO 350
176       ENDIF
177 C
178 C          Select W decay mode
179 C          QMW dependence neglected in branching ratios
180 C          BRANCH is cum. br. with heavy modes subtracted.
181 C
182       IF(KEYS(3)) THEN
183         BRANCH(1)=0.
184         SUMBR=0.
185         DO 105 IQ1=2,25
186           IQ2=MATCH(IQ1,JWTYP)
187           IF(IQ2.EQ.0) THEN
188             BRMODE=0.
189           ELSE
190             BRMODE=WCBR(IQ1,JWTYP)-WCBR(IQ1-1,JWTYP)
191             IFL1=LISTJ(IQ1)
192             IFL2=LISTJ(IQ2)
193             IF(S12.LE.(AMASS(IFL1)+AMASS(IFL2))**2) BRMODE=0.
194           ENDIF
195           BRANCH(IQ1)=BRANCH(IQ1-1)+BRMODE
196           SUMBR=SUMBR+BRMODE
197 105     CONTINUE
198         BRINV=1./SUMBR
199 C
200         TRY=RANF()
201         DO 110 IQ=1,25
202           IF(TRY.LT.BRANCH(IQ)*BRINV.AND.MATCH(IQ,JWTYP).NE.0) THEN
203             JETTYP(1)=IQ
204             JETTYP(2)=MATCH(IQ,JWTYP)
205             GO TO 120
206           ENDIF
207 110     CONTINUE
208       ENDIF
209 C
210 120   IF(GOMSSM) THEN
211         IFL1=LISTSS(JETTYP(1))
212         IFL2=LISTSS(JETTYP(2))
213       ELSE
214         IFL1=LISTJ(JETTYP(1))
215         IFL2=LISTJ(JETTYP(2))
216       ENDIF
217 C
218 C          Select masses of decay products. These are just normal masses
219 C          except for Z Z* decay of Higgs, where one is virtual.
220 C
221       EML(1)=AMASS(IFL1)
222       EML(2)=AMASS(IFL2)
223       IF(KEYS(7).AND.EML(1)+EML(2).GT.QMW) THEN
224 C          WW* or ZZ* decay - generate/check W* or Z* mass
225         IF((IABS(IFL1).EQ.80.AND.IABS(IFL2).EQ.80)
226      $  .OR.(IFL1.EQ.90.AND.IFL2.EQ.90)) THEN
227           IZSTAR=3-2*RANF()
228           IF(GOMSSM) THEN
229             JVIR=JETTYP(IZSTAR)-76
230           ELSE
231             JVIR=JETTYP(IZSTAR)-25
232           ENDIF
233           EML(IZSTAR)=ZZSTAR(QMW,JVIR)
234           IF(EML(IZSTAR).LT.ZSTARS(JVIR,IZSTAR)) GO TO 200
235 C          Other decay - invalid for this QMW
236         ELSE
237           GO TO 200
238         ENDIF
239       ENDIF
240 C
241 C          Generate W decay in its rest frame and compare with SIGDY2.
242 C          First set up momenta of decay products:
243 C
244       EMSQL(1)=EML(1)**2
245       EMSQL(2)=EML(2)**2
246       EL(1)=(S12+EMSQL(1)-EMSQL(2))/(2.*QMW)
247       EL(2)=(S12+EMSQL(2)-EMSQL(1))/(2.*QMW)
248       PL12=SQRT((S12-(EML(1)+EML(2))**2)*(S12-(EML(1)-EML(2))**2))
249      $/(2.*QMW)
250 C          W momentum
251       PREST(1)=QTW*COS(PHIW)
252       PREST(2)=QTW*SIN(PHIW)
253       PREST(3)=QPW
254       PREST(4)=SQRT(QW**2+QMW**2)
255       PREST(5)=QMW
256       NTRY2=0
257 C          Generate next W decay
258 20    CONTINUE
259       NTRY2=NTRY2+1
260       IF(NTRY2.GT.NTRIES) GO TO 999
261       COSTHL=2.*RANF()-1.
262       THL=ACOS(COSTHL)
263       PHL=2.*PI*RANF()
264       PTL=PL12*SIN(THL)
265 C
266       DO 300 I=1,2
267         SGN=3-2*I
268         PL(1)=SGN*PTL*COS(PHL)
269         PL(2)=SGN*PTL*SIN(PHL)
270         PL(3)=SGN*PL12*COSTHL
271         PL(4)=EL(I)
272         PL(5)=EML(I)
273 C          Boost with W momentum
274         BP=0.
275         DO 310 K=1,3
276 310     BP=BP+PL(K)*PREST(K)
277         BP=BP/PREST(5)
278         DO 320 K=1,3
279 320     PL(K)=PL(K)+PREST(K)*PL(4)/PREST(5)
280      $  +PREST(K)*BP/(PREST(4)+PREST(5))
281         PL(4)=PL(4)*PREST(4)/PREST(5)+BP
282 C          Fill common blocks
283         PT(I)=SQRT(PL(1)**2+PL(2)**2)
284         P(I)=SQRT(PT(I)**2+PL(3)**2)
285         IF(PT(I).GT.0.) THEN
286           PHI(I)=ATAN2(PL(2),PL(1))
287         ELSE
288           PHI(I)=(I-1)*PI
289         ENDIF
290         IF(PHI(I).LT.0.) PHI(I)=PHI(I)+2.*PI
291         CTH(I)=PL(3)/P(I)
292         STH(I)=PT(I)/P(I)
293         TH(I)=ACOS(CTH(I))
294         XJ(I)=PL(3)/HALFE
295         IF(CTH(I).GT.0.) THEN
296           PLPL=PL(4)+PL(3)
297           PLMN=(PT(I)**2+EMSQL(I))/PLPL
298         ELSE
299           PLMN=PL(4)-PL(3)
300           PLPL=(PT(I)**2+EMSQL(I))/PLMN
301         ENDIF
302         YJ(I)=.5*ALOG(PLPL/PLMN)
303 300   CONTINUE
304 C
305 C          Test cross section
306 C          Extra kinematics for W+W->W+W
307 C
308       IF(KEYS(7).OR.KEYS(9)) THEN
309         SHAT=S12
310         IF(GOMSSM) THEN
311           AMINI=AMASS(LISTSS(INITYP(1)))
312         ELSE
313           AMINI=AMASS(LISTJ(INITYP(1)))
314         ENDIF
315         AMFIN=EML(1)
316         PINI=.5*SQRT(S12-4.*AMINI**2)
317         PFIN=PL12
318         THAT=AMINI**2+AMFIN**2-.5*S12+2.*PINI*PFIN*COSTHL
319         UHAT=AMINI**2+AMFIN**2-.5*S12-2.*PINI*PFIN*COSTHL
320       ENDIF
321 C
322 C          Check W decay
323 C
324       IF(NOGOOD(3)) GO TO 20
325 C
326 C          Check W decay with kinematic limits
327 C
328       IF(NOGOOD(4)) GO TO 200
329 350   NKEEP=NKEEP+1
330 C
331 C            Set PBEAM
332 C
333       PBEAM(1)=(1.-X1)*HALFE
334       PBEAM(2)=(1.-X2)*HALFE
335       IF(NJET.LT.3) GO TO 502
336       IFL=LISTJ(JETTYP(3))
337       EMSQL(3)=AMASS(IFL)**2
338 502   CONTINUE
339 C
340 C          Set PJETS
341 C
342       IF(KEYS(11)) THEN
343         N0J=3
344       ELSE
345         N0J=1
346       ENDIF
347       DO 501 I=N0J,NJET
348         PJETS(3,I)=P(I)*CTH(I)
349         PJETS(1,I)=PT(I)*COS(PHI(I))
350         PJETS(2,I)=PT(I)*SIN(PHI(I))
351         PJETS(4,I)=SQRT(P(I)**2+EMSQL(I))
352         PJETS(5,I)=SQRT(EMSQL(I))
353         IF(KEYS(7).AND.GOMSSM) THEN
354           IDJETS(I)=LISTSS(JETTYP(I))
355         ELSE
356           IDJETS(I)=LISTJ(JETTYP(I))
357         ENDIF
358 501   CONTINUE
359 C          No technicolor IDENT's defined, so...
360       IF(KEYS(3)) THEN
361         IDENTW=LISTW(JWTYP)
362       ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
363         IDENTW=81
364       ELSEIF(KEYS(7).AND.GOMSSM) THEN
365         IDENTW=IHTYPE
366       ELSEIF(KEYS(11)) THEN
367         IDENTW=92
368       ELSE
369         IDENTW=0
370       ENDIF
371 C          W momentum in /PJETS/
372       IF(KEYS(11)) THEN
373         QWJET(1)=QTW*COS(PHIW)
374         QWJET(2)=QTW*SIN(PHIW)
375         QWJET(3)=QPW
376         QWJET(4)=SQRT(QW**2+QMW**2)
377         QWJET(5)=QMW
378       ELSE
379         DO 503 K=1,4
380 503     QWJET(K)=PJETS(K,1)+PJETS(K,2)
381         QWJET(5)=QMW
382       ENDIF
383 C
384 C          Set PINITS
385       DO 504 I=1,2
386         IF(KEYS(7).AND.GOMSSM) THEN
387           IDINIT(I)=LISTSS(INITYP(I))
388         ELSE
389           IDINIT(I)=LISTJ(INITYP(I))
390         ENDIF
391         PINITS(5,I)=AMASS(IDINIT(I))
392         PINITS(1,I)=0.
393         PINITS(2,I)=0.
394 504   CONTINUE
395 C          Calculate total momentum
396       QPL=QWJET(4)+QWJET(3)
397       QMN=QWJET(4)-QWJET(3)
398       IF(NJET.EQ.3) THEN
399         QPL=QPL+PJETS(4,3)+PJETS(3,3)
400         QMN=QMN+PJETS(4,3)-PJETS(3,3)
401       ENDIF
402 C          and solve initial kinematics
403       AM1SQ=PINITS(5,1)**2
404       AM2SQ=PINITS(5,2)**2
405       ROOT=SQRT((QPL*QMN-AM1SQ-AM2SQ)**2-4.*AM1SQ*AM2SQ)
406       P1PL=(QPL*QMN+AM1SQ-AM2SQ+ROOT)/(2.*QMN)
407       P1MN=AM1SQ/P1PL
408       P2MN=(QPL*QMN+AM2SQ-AM1SQ+ROOT)/(2.*QPL)
409       P2PL=AM2SQ/P2MN
410       PINITS(3,1)=.5*(P1PL-P1MN)
411       PINITS(4,1)=.5*(P1PL+P1MN)
412       PINITS(3,2)=.5*(P2PL-P2MN)
413       PINITS(4,2)=.5*(P2PL+P2MN)
414       RETURN
415 C
416 999   CALL PRTEVT(0)
417       WRITE(ITLIS,9999) NTRIES
418 9999  FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE AN',
419      C' EVENT. CHECK LIMITS OR INCREASE NTRIES')
420       STOP 99
421       END