1 #include "isajet/pilot.h"
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
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.
13 C Note that WT contains various constant factors that were
16 C Jacobean for dQ = (EHMG-SUM)*dRANF
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. :-(
25 #if defined(CERNLIB_IMPNONE)
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"
36 C Local variables; MXJETS defined in /PJETS/
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
43 INTEGER NJET1,I,JJ1,NTRY,J,JSAVE,II,K
48 PCM(A,B,C)=SQRT((A-B-C)*(A+B+C)*(A-B+C)*(A+B-C))/(2.*A)
51 C Generate COM mass and rapidity
55 SHMG=EHMGMN**2+(EHMGMX**2-EHMGMN**2)*RANF()
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
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))
73 PJETS8(3,1)=SYHMG*E1+CYHMG*P12
74 PJETS8(0,1)=CYHMG*E1+SYHMG*P12
77 PJETS8(3,2)=SYHMG*E1-CYHMG*P12
78 PJETS8(0,2)=CYHMG*E1-SYHMG*P12
85 IF(SUM.GE.EHMG) GO TO 999
88 C Generate masses for uniform NJET-body phase space.
93 IF(NTRY.GT.NTRIES) THEN
94 WRITE(ITLIS,9999) NTRY
95 9999 FORMAT(//2X,'ERROR IN MULJET ... NTRY = ',I8)
104 IF(RNEW.LE.RND(J)) GO TO 210
109 C Jacobean for d(shmg)d(yhmg) and overall 1/(2*shmg)
110 WT=(EHMGMX**2-EHMGMN**2)*(YHMGMX-YHMGMN)/(2*SHMG)
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)
119 C Jacobean for final 2-body decay differs by this factor
120 WT=WT*PI/(DELTAQ*EHMG)
122 C Carry out 2-body decays
125 QCM=PCM(AMGEN(I),AMGEN(I+1),AMJET8(I+2))
128 U(1)=SQRT(1-U(3)**2)*COS(PHI)
129 U(2)=SQRT(1-U(3)**2)*SIN(PHI)
131 PJETS8(J,I+2)=QCM*U(J)
132 PGEN(J,I+1)=-PJETS8(J,I+2)
134 PJETS8(0,I+2)=SQRT(QCM**2+AMJET8(I+2)**2)
135 PGEN(0,I+1)=SQRT(QCM**2+AMGEN(I+1)**2)
139 PJETS8(J,NJET+2)=PGEN(J,NJET)
142 C Boost PGEN frames to lab frame.
147 BETA(J)=PGEN(J,I)/PGEN(0,I)
149 GAMMA=PGEN(0,I)/AMGEN(I)
151 BP=BETA(1)*PJETS8(1,K+2)+BETA(2)*PJETS8(2,K+2)+
152 $ BETA(3)*PJETS8(3,K+2)
154 PJETS8(J,K+2)=PJETS8(J,K+2)+GAMMA*BETA(J)*(PJETS8(0,K+2)
155 $ +BP*GAMMA/(GAMMA+1.))
157 PJETS8(0,K+2)=GAMMA*(PJETS8(0,K+2)+BP)
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
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