More volume overlaps corrected
[u/mrichter/AliRoot.git] / ISAJET / code / muljet.F
1 #include "isajet/pilot.h"
2       SUBROUTINE MULJET(WT)
3 C
4 C          Using masses from /MGKIN/, generate NJET<=MXJETS body phase 
5 C          space point satisfying cuts:
6 C          (1) Generate kinematic point using successive 2-body decays
7 C              with Jacobean
8 C              dPhi_N(p1...pN) = dQ1 p1/(4*pi) dPhi_(N-1)(q1...pN)
9 C          (2) Apply individual jet cuts from /JETLIM/ and dijet
10 C              cuts from /MGLIMS/ to ensure IR-safe cross section.
11 C          (3) Return weight WT or 0 if outside limits. 
12 C
13 C          Note that WT contains various constant factors that were 
14 C          dropped in DECAY:
15 C            1/(2*SHMG) Jacobean
16 C            Jacobean for dQ = (EHMG-SUM)*dRANF
17 C            Factors of 4pi
18 C
19 C          MadGraph/Helas notation: 
20 C          PJETS8(0:3,1:2) = initial momenta
21 C          PJETS8(0:3,3:NJET+2) = final momenta
22 C          Note: ANSI extensions, e.g. REAL*8 P(0:3) are required for 
23 C          compatibility with Helas and MadGraph. :-(
24 C
25 #if defined(CERNLIB_IMPNONE)
26       IMPLICIT NONE
27 #endif
28 C
29 #include "isajet/itapes.inc"
30 #include "isajet/jetlim.inc"
31 #include "isajet/mglims.inc"
32 #include "isajet/pjets.inc"
33 #include "isajet/mgkin.inc"
34 #include "isajet/primar.inc"
35 C
36 C          Local variables; MXJETS defined in /PJETS/
37 C
38       REAL*8 PGEN(0:3,MXJETS),AMGEN(MXJETS),RND(MXJETS)
39       REAL*8 SHMG,EHMG,YHMG,SUM,SUM1,RNEW,WT,QCM,PI,
40      $U(3),PHI,BETA(3),GAMMA,BP,PTI,PPI,YI,XJI,PHII,AMIJ
41       REAL*8 CYHMG,SYHMG,E1,E2,P12,DELTAQ
42       REAL*8 PCM,A,B,C
43       INTEGER NJET1,I,JJ1,NTRY,J,JSAVE,II,K
44       REAL RANF
45 C
46 C          Function definition
47 C
48       PCM(A,B,C)=SQRT((A-B-C)*(A+B+C)*(A-B+C)*(A+B-C))/(2.*A)
49 C
50 C
51 C          Generate COM mass and rapidity
52 C
53       PI=4.D0*DATAN(1.D0)
54 100   CONTINUE
55       SHMG=EHMGMN**2+(EHMGMX**2-EHMGMN**2)*RANF()
56       EHMG=SQRT(SHMG)
57       YHMG=YHMGMN+(YHMGMX-YHMGMN)*RANF()
58       IF(EHMG*EXP(ABS(YHMG)).GT.ECM) GO TO 999
59       IF(EHMG.LT.AMJET8(1)+AMJET8(2)) GO TO 999
60       CYHMG=DCOSH(YHMG)
61       SYHMG=SINH(YHMG)
62       AMGEN(1)=EHMG
63       PGEN(1,1)=0
64       PGEN(2,1)=0
65       PGEN(3,1)=EHMG*SYHMG
66       PGEN(0,1)=EHMG*CYHMG
67       E1=(EHMG**2+AMJET8(1)**2-AMJET8(2)**2)/(2*EHMG)
68       E2=(EHMG**2-AMJET8(1)**2+AMJET8(2)**2)/(2*EHMG)
69       P12=PCM(EHMG,AMJET8(1),AMJET8(2))
70 C          Initial momenta
71       PJETS8(1,1)=0
72       PJETS8(2,1)=0
73       PJETS8(3,1)=SYHMG*E1+CYHMG*P12
74       PJETS8(0,1)=CYHMG*E1+SYHMG*P12
75       PJETS8(1,2)=0
76       PJETS8(2,2)=0
77       PJETS8(3,2)=SYHMG*E1-CYHMG*P12
78       PJETS8(0,2)=CYHMG*E1-SYHMG*P12
79 C
80       NJET1=NJET-1
81       SUM=0
82       DO 110 I=1,NJET
83         SUM=SUM+AMJET8(I+2)
84 110   CONTINUE
85       IF(SUM.GE.EHMG) GO TO 999
86       DELTAQ=EHMG-SUM
87 C
88 C          Generate masses for uniform NJET-body phase space.
89 C
90       NTRY=0
91 200   CONTINUE
92       NTRY=NTRY+1
93       IF(NTRY.GT.NTRIES) THEN
94         WRITE(ITLIS,9999) NTRY
95 9999    FORMAT(//2X,'ERROR IN MULJET ... NTRY = ',I8)
96         STOP99
97       ENDIF
98       RND(1)=1
99       DO 210 I=2,NJET1
100         RNEW=RANF()
101         DO 220 JJ1=1,I-1
102           J=I-JJ1
103           JSAVE=J+1
104           IF(RNEW.LE.RND(J)) GO TO 210
105           RND(JSAVE)=RND(J)
106 220     CONTINUE
107 210   RND(JSAVE)=RNEW
108       RND(NJET)=0
109 C          Jacobean for d(shmg)d(yhmg) and overall 1/(2*shmg)
110       WT=(EHMGMX**2-EHMGMN**2)*(YHMGMX-YHMGMN)/(2*SHMG)
111       SUM1=SUM
112       DO 230 I=2,NJET
113         SUM1=SUM1-AMJET8(I-1+2)
114         AMGEN(I)=SUM1+RND(I)*(AMGEN(1)-SUM)
115         IF(AMGEN(I-1).LE.AMGEN(I)+AMJET8(I-1+2)) GO TO 200
116 C          Jacobean for sigma_n -> sigma_n-1
117         WT=WT*PCM(AMGEN(I-1),AMGEN(I),AMJET8(I-1+2))*DELTAQ/(4*PI**2)
118 230   CONTINUE
119 C          Jacobean for final 2-body decay differs by this factor
120       WT=WT*PI/(DELTAQ*EHMG)
121 C
122 C          Carry out 2-body decays
123 C
124       DO 310 I=1,NJET1
125         QCM=PCM(AMGEN(I),AMGEN(I+1),AMJET8(I+2))
126         U(3)=2.*RANF()-1
127         PHI=2*PI*RANF()
128         U(1)=SQRT(1-U(3)**2)*COS(PHI)
129         U(2)=SQRT(1-U(3)**2)*SIN(PHI)
130         DO 320 J=1,3
131           PJETS8(J,I+2)=QCM*U(J)
132           PGEN(J,I+1)=-PJETS8(J,I+2)
133 320     CONTINUE
134         PJETS8(0,I+2)=SQRT(QCM**2+AMJET8(I+2)**2)
135         PGEN(0,I+1)=SQRT(QCM**2+AMGEN(I+1)**2)
136 310   CONTINUE
137 C
138       DO 330 J=0,3
139         PJETS8(J,NJET+2)=PGEN(J,NJET)
140 330   CONTINUE
141 C
142 C          Boost PGEN frames to lab frame.
143 C
144       DO 400 II=1,NJET1
145         I=NJET-II
146         DO 410 J=1,3
147           BETA(J)=PGEN(J,I)/PGEN(0,I)
148 410     CONTINUE
149         GAMMA=PGEN(0,I)/AMGEN(I)
150         DO 420 K=I,NJET
151           BP=BETA(1)*PJETS8(1,K+2)+BETA(2)*PJETS8(2,K+2)+
152      $    BETA(3)*PJETS8(3,K+2)
153           DO 430 J=1,3
154             PJETS8(J,K+2)=PJETS8(J,K+2)+GAMMA*BETA(J)*(PJETS8(0,K+2)
155      $      +BP*GAMMA/(GAMMA+1.))
156 430       CONTINUE
157           PJETS8(0,K+2)=GAMMA*(PJETS8(0,K+2)+BP)
158 420     CONTINUE
159 400   CONTINUE
160 C
161 C          Check limits
162 C
163       DO 500 I=1,NJET
164         PTI=SQRT(PJETS8(1,I+2)**2+PJETS8(2,I+2)**2)
165         IF(PTI.LE.PTMIN(I).OR.PTI.GE.PTMAX(I)) GO TO 999
166         PPI=SQRT(PTI**2+PJETS8(3,I+2)**2)
167         IF(PPI.LE.PMIN(I).OR.PPI.GE.PMAX(I)) GO TO 999
168         XJI=PJETS8(3,I+2)/PPI
169         IF(XJI.LE.XJMIN(I).OR.XJI.GE.XJMAX(I)) GO TO 999
170         PHII=ATAN2(PJETS8(2,I+2),PJETS8(1,I+2))
171         IF(PHII.LT.0) PHII=PHII+2*PI
172         IF(PHII.LE.PHIMIN(I).OR.PHII.GE.PHIMAX(I)) GO TO 999
173         YI=-LOG(TAN(ACOS(XJI)/2))
174         IF(YI.LE.YJMIN(I).OR.YI.GE.YJMAX(I)) GO TO 999
175 500   CONTINUE
176 C
177       DO 510 I=1,NJET
178         DO 520 J=I+1,NJET
179           AMIJ=(PJETS8(0,I+2)+PJETS8(0,J+2))**2
180      $        -(PJETS8(1,I+2)+PJETS8(1,J+2))**2
181      $        -(PJETS8(2,I+2)+PJETS8(2,J+2))**2
182      $        -(PJETS8(3,I+2)+PJETS8(3,J+2))**2
183           AMIJ=SIGN(SQRT(ABS(AMIJ)),AMIJ)
184           IF(AMIJ.LE.AMIJMN(I,J).OR.AMIJ.GE.AMIJMX(I,J)) GO TO 999
185 520     CONTINUE
186 510   CONTINUE
187 C
188       RETURN
189 C
190 999   WT=0
191       RETURN
192       END