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