]>
Commit | Line | Data |
---|---|---|
0795afa3 | 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 |