ISAJET/code/twojet.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE TWOJET
   3 C
   4 C          Driving routine to generate initial parameters for jets,
   5 C          assuming zero initial transverse momentum, ie PT(1)=PT(2).
   6 C
   7 C          Parameters are PT,YJ,PHI with P,YJ,XJ as dependent variables,
   8 C          where YJ=RAPIDITY, XJ=Feynman X.
   9 C          All parameters are stored in COMMON/JETPAR/.
  10 C          Cross section is called from NOGOOD.
  11 C
  12 #if defined(CERNLIB_IMPNONE)
  13       IMPLICIT NONE
  14 #endif
  15 #include "isajet/idrun.inc"
  16 #include "isajet/itapes.inc"
  17 #include "isajet/keys.inc"
  18 #include "isajet/mbpar.inc"
  19 #include "isajet/pjets.inc"
  20 #include "isajet/pinits.inc"
  21 #include "isajet/jetlim.inc"
  22 #include "isajet/ptpar.inc"
  23 #include "isajet/jetpar.inc"
  24 #include "isajet/primar.inc"
  25 #include "isajet/partcl.inc"
  26 #include "isajet/const.inc"
  27 #include "isajet/jetsig.inc"
  28 #include "isajet/totals.inc"
  29 #include "isajet/isloop.inc"
  30 #include "isajet/sstype.inc"
  31 #include "isajet/xmssm.inc"
  32 C
  33       REAL ACOSH,XXX,WTFCN,PPP,RANF,SIGN,SGN,AMQ1,AMASS,AMQ2
  34       REAL PPLUS,PMINUS,PSUM3,PSUM4,PPL,PMN,SQ1,SQ2,ROOT,P1PL,P1MN
  35       REAL P2PL,P2MN,AMI1,AMI2
  36       INTEGER NREJ,I,II,IS,IFL1,IFL2
  37       REAL X(2)
  38       EQUIVALENCE (X(1),X1)
  39       LOGICAL NOGOOD
  40       LOGICAL YGENJ
  41       INTEGER LISTJ(17),LISTW(4),LISTSS(85),LISTSM(30)
  42 C
  43 C          SUSY IDENT codes from /SSTYPE/. (Fortran 77 allows - signs
  44 C          in parameter statements but not data statements.)
  45       INTEGER MSUPL,MSDNL,MSSTL,MSCHL,MSBT1,MSTP1,
  46      $MSUPR,MSDNR,MSSTR,MSCHR,MSBT2,MSTP2,MSW1,MSW2,
  47      $MSNEL,MSEL,MSNML,MSMUL,MSNTL,MSTAU1,MSER,MSMUR,MSTAU2
  48       PARAMETER (MSUPL=-ISUPL)
  49       PARAMETER (MSDNL=-ISDNL)
  50       PARAMETER (MSSTL=-ISSTL)
  51       PARAMETER (MSCHL=-ISCHL)
  52       PARAMETER (MSBT1=-ISBT1)
  53       PARAMETER (MSTP1=-ISTP1)
  54       PARAMETER (MSUPR=-ISUPR)
  55       PARAMETER (MSDNR=-ISDNR)
  56       PARAMETER (MSSTR=-ISSTR)
  57       PARAMETER (MSCHR=-ISCHR)
  58       PARAMETER (MSBT2=-ISBT2)
  59       PARAMETER (MSTP2=-ISTP2)
  60       PARAMETER (MSW1=-ISW1)
  61       PARAMETER (MSW2=-ISW2)
  62       PARAMETER (MSNEL=-ISNEL)
  63       PARAMETER (MSEL=-ISEL)
  64       PARAMETER (MSNML=-ISNML)
  65       PARAMETER (MSMUL=-ISMUL)
  66       PARAMETER (MSNTL=-ISNTL)
  67       PARAMETER (MSTAU1=-ISTAU1)
  68       PARAMETER (MSER=-ISER)
  69       PARAMETER (MSMUR=-ISMUR)
  70       PARAMETER (MSTAU2=-ISTAU2)
  71 C
  72       DATA LISTSS/ISGL,
  73      $ISUPL,MSUPL,ISDNL,MSDNL,ISSTL,MSSTL,ISCHL,MSCHL,ISBT1,MSBT1,
  74      $ISTP1,MSTP1,
  75      $ISUPR,MSUPR,ISDNR,MSDNR,ISSTR,MSSTR,ISCHR,MSCHR,ISBT2,MSBT2,
  76      $ISTP2,MSTP2,
  77      $ISW1,MSW1,ISW2,MSW2,ISZ1,ISZ2,ISZ3,ISZ4,
  78      $ISNEL,MSNEL,ISEL,MSEL,ISNML,MSNML,ISMUL,MSMUL,ISNTL,MSNTL,
  79      $ISTAU1,MSTAU1,ISER,MSER,ISMUR,MSMUR,ISTAU2,MSTAU2,
  80      $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,11,-11,12,-12,13,-13,
  81      $14,-14,15,-15,16,-16,10,80,-80,90,82,83,84,86,-86/
  82       DATA LISTSM/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,11,-11,12,-12,13,-13,
  83      $14,-14,15,-15,16,-16,10,80,-80,90,81/
  84       DATA LISTJ/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7,8,-8/
  85       DATA LISTW/10,80,-80,90/
  86 C          Inverse hyperbolic cosine function
  87       ACOSH(XXX)=ALOG(XXX+SQRT(XXX**2-1.))
  88       WTFCN(PPP)=2.*PPP*PTGEN2*PTGEN3*PPP**((PTGEN3-1.)/PTGEN3)
  89 C
  90 C          Initialize
  91 C
  92       NPTCL=0
  93       PHI(1)=PHIMIN(1)+(PHIMAX(1)-PHIMIN(1))*RANF()
  94       PHI(2)=AMOD(PHI(1)+PI,2.*PI)
  95       NREJ=-1
  96       SIGMA=0.
  97       WT=1.
  98       IF(.NOT.FIXPT(2)) GOTO 101
  99       FIXPT(1)=.TRUE.
 100       PT(1)=PT(2)
 101   101 CONTINUE
 102       IF(FIXPT(1)) GOTO 400
 103       DO 110 I=1,2
 104       IF(FIXP(I)) GOTO 200
 105       IF(FIXXJ(I)) GOTO 300
 106   110 CONTINUE
 107 C
 108 C          Genetate PT and YJ with no variables fixed
 109 C
 110   111 NREJ=NREJ+1
 111       IF(NREJ.GT.NTRIES) GO TO 910
 112       SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
 113       NKINPT=NKINPT+1
 114       SIGMA=0.
 115       WT=1.
 116 C            Generate PT with a power law distribution
 117       PT(1)=(PTGEN1+PTGEN2*RANF())**PTGEN3
 118       PT(2)=PT(1)
 119       SIGMAX=PTFUN1*PT(1)**PTFUN2
 120 C          GENERATE FLAT IN YJ, CALCULATE CORRESPONDING TH
 121       DO 115 I=1,2
 122       IF(FIXYJ(I)) GOTO 115
 123       IF(.NOT.YGENJ(I)) GOTO 111
 124   115 CONTINUE
 125       DO 116 I=1,2
 126       P(I)=PT(I)/STH(I)
 127       IF(P(I).LT.PMIN(I).OR.P(I).GT.PMAX(I)) GOTO 111
 128       XJ(I)=P(I)*CTH(I)/HALFE
 129       IF(XJ(I).LT.XJMIN(I).OR.XJ(I).GT.XJMAX(I)) GOTO 111
 130   116 CONTINUE
 131       WT=WT*WTFCN(PT(1))
 132       IF(NOGOOD(1)) GOTO 111
 133       SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
 134       NKEEP=NKEEP+1
 135       GO TO 500
 136 C
 137 C          Generate PT and YJ fixing P
 138 C
 139   200 CONTINUE
 140       II=3-I
 141   211 NREJ=NREJ+1
 142       IF(NREJ.GT.NTRIES) GO TO 910
 143       NKINPT=NKINPT+1
 144       WT=0.
 145       IF(FIXYJ(I)) GOTO 212
 146 C          Generate PT with a power law distribution
 147       PT(1)=(PTGEN1+PTGEN2*RANF())**PTGEN3
 148       SIGMAX=PTFUN1*PT(1)**PTFUN2
 149       PT(2)=PT(1)
 150 C          Given PT, TH is fixed except for a sign
 151       STH(I)=PT(I)/P(I)
 152       SIGN=1.0
 153       IF(RANF().GT.0.5) SIGN=-1.0
 154       CTH(I)=SIGN*SQRT(1.-STH(I)**2)
 155       TH(I)=ATAN2(STH(I),CTH(I))
 156       YJ(I)=-ALOG(TAN(TH(I)/2.))
 157       IF(YJ(I).LT.YJMIN(I).OR.YJ(I).GT.YJMAX(I)) GOTO 211
 158       GOTO 213
 159   212 PT(1)=P(I)*STH(I)
 160   213 CONTINUE
 161       XJ(I)=P(I)*CTH(I)/HALFE
 162       IF(XJ(I).LT.XJMIN(I).OR.XJ(I).GT.XJMAX(I)) GOTO 211
 163       IF(FIXP(II)) GOTO 220
 164       IF(FIXXJ(II)) GOTO 230
 165       IF(FIXYJ(II)) GOTO 215
 166       IF(.NOT.YGENJ(II)) GOTO 211
 167   215 CONTINUE
 168       P(II)=PT(II)/STH(II)
 169       IF(P(II).LT.PMIN(II).OR.P(II).GT.PMAX(II)) GOTO 211
 170       XJ(II)=P(II)*CTH(II)/HALFE
 171       IF(XJ(II).LT.XJMIN(II).OR.XJ(II).GT.XJMAX(II)) GOTO 211
 172       GOTO 250
 173 220   STH(II)=PT(II)/P(II)
 174       SGN=1.0
 175       IF(RANF().GT.0.5) SGN=-1.0
 176       CTH(II)=SGN*SQRT(1.-STH(II)**2)
 177       TH(II)=ATAN2(STH(II),CTH(II))
 178       YJ(II)=-ALOG(TAN(TH(II)/2.))
 179       IF(YJ(II).LT.YJMIN(II).OR.YJ(II).GT.YJMAX(II)) GOTO 211
 180       XJ(II)=P(II)*CTH(II)/HALFE
 181       IF(XJ(II).LT.XJMIN(II).OR.XJ(II).GT.XJMAX(II)) GOTO 211
 182       GOTO 250
 183   230 TH(II)=ATAN2(PT(II),XJ(II)*HALFE)
 184       YJ(II)=-ALOG(TAN(TH(II)/2.))
 185       IF(YJ(II).LT.YJMIN(II).OR.YJ(II).GT.YJMAX(II)) GOTO 211
 186       CTH(II)=COS(TH(II))
 187       STH(II)=SIN(TH(II))
 188   250 CONTINUE
 189       IF(NOGOOD(1)) GOTO 211
 190       NKEEP=NKEEP+1
 191       GO TO 500
 192 C
 193 C          Generate PT and YJ at fixed XJ
 194 C
 195   300 CONTINUE
 196       II=3-I
 197   311 NREJ=NREJ+1
 198       IF(NREJ.GT.NTRIES) GO TO 910
 199       NKINPT=NKINPT+1
 200       WT=0.
 201 C          Generate PT with a power law distribution
 202       PT(1)=(PTGEN1+PTGEN2*RANF())**PTGEN3
 203       SIGMAX=PTFUN1*PT(1)**PTFUN2
 204       PT(2)=PT(1)
 205       TH(I)=ATAN2(PT(I),XJ(I)*HALFE)
 206       YJ(I)=-ALOG(TAN(TH(I)/2.))
 207       IF(YJ(I).LT.YJMIN(I).OR.YJ(I).GT.YJMAX(I)) GOTO 311
 208       CTH(I)=COS(TH(I))
 209       STH(I)=SIN(TH(I))
 210       P(I)=PT(I)/STH(I)
 211       IF(FIXYJ(II)) GOTO 315
 212       IF(FIXP(II)) GOTO 314
 213       YJ(II)=YJMIN(II)+(YJMAX(II)-YJMIN(II))*RANF()
 214       TH(II)=2.*ATAN(EXP(-YJ(II)))
 215       CTH(II)=COS(TH(II))
 216       STH(II)=SIN(TH(II))
 217       GOTO 315
 218   314 CONTINUE
 219       STH(II)=PT(II)/P(II)
 220       CTH(II)=SQRT(1.-STH(II)**2)
 221       IF(RANF().GT.0.5) CTH(II)=-CTH(II)
 222       TH(II)=ATAN2(STH(II),CTH(II))
 223       YJ(II)=-ALOG(TAN(TH(II)/2.))
 224   315 CONTINUE
 225       P(II)=PT(II)/STH(II)
 226       XJ(II)=P(II)*CTH(II)/HALFE
 227       IF(XJ(II).LT.XJMIN(II).OR.XJ(II).GT.XJMAX(II)) GOTO 311
 228       IF(NOGOOD(1)) GOTO 311
 229       NKEEP=NKEEP+1
 230       GO TO 500
 231 C
 232 C          Generate YJ at fixed PT
 233 C
 234   400 CONTINUE
 235       PT(2)=PT(1)
 236   411 NREJ=NREJ+1
 237       IF(NREJ.GT.NTRIES) GO TO 910
 238       NKINPT=NKINPT+1
 239       WT=0.
 240       DO 415 I=1,2
 241       IF(FIXYJ(I)) GOTO 415
 242       IF(FIXP(I)) GOTO 413
 243       IF(.NOT.YGENJ(I)) GO TO 411
 244       GOTO 414
 245   413 CONTINUE
 246       IS=1
 247       IF(RANF().GT.0.5) IS=2
 248       CTH(I)=CTHS(IS,I)
 249       TH(I)=THS(IS,I)
 250       YJ(I)=YJS(IS,I)
 251   414 CONTINUE
 252       P(I)=PT(I)/STH(I)
 253       XJ(I)=P(I)*CTH(I)/HALFE
 254   415 CONTINUE
 255       IF(NOGOOD(1)) GOTO 411
 256       NKEEP=NKEEP+1
 257 C
 258 C          Reset /JETPAR/
 259 C
 260   500 CONTINUE
 261       IF(KEYS(1)) THEN
 262         IFL1=LISTJ(JETTYP(1))
 263         IFL2=LISTJ(JETTYP(2))
 264         AMQ1=AMASS(IFL1)
 265         AMQ2=AMASS(IFL2)
 266         AMI1=AMASS(LISTJ(INITYP(1)))
 267         AMI2=AMASS(LISTJ(INITYP(2)))
 268         CALL TWOKIN(AMI1,AMI2,AMQ1,AMQ2)
 269       ELSEIF(KEYS(5).OR.(KEYS(10).AND.GOMSSM)) THEN
 270         IFL1=LISTSS(JETTYP(1))
 271         IFL2=LISTSS(JETTYP(2))
 272         AMQ1=AMASS(IFL1)
 273         AMQ2=AMASS(IFL2)
 274         CALL TWOKIN(0.,0.,AMQ1,AMQ2)
 275       ELSEIF(KEYS(6)) THEN
 276         IFL1=LISTW(JETTYP(1))
 277         IFL2=LISTW(JETTYP(2))
 278         AMQ1=AMASS(IFL1)
 279         AMQ2=AMASS(IFL2)
 280         CALL TWOKIN(0.,0.,AMQ1,AMQ2)
 281       ELSEIF(KEYS(8)) THEN
 282         IF(JETTYP(1).LE.13) THEN
 283           IFL1=LISTJ(JETTYP(1))
 284         ELSE
 285           IFL1=10
 286         ENDIF
 287         IF(JETTYP(2).LE.13) THEN
 288           IFL2=LISTJ(JETTYP(2))
 289         ELSE
 290           IFL2=10
 291         ENDIF
 292         AMQ1=AMASS(IFL1)
 293         AMQ2=AMASS(IFL2)
 294         CALL TWOKIN(0.,0.,AMQ1,AMQ2)
 295       ELSEIF(KEYS(10).AND.(.NOT.GOMSSM)) THEN
 296         IFL1=LISTSM(JETTYP(1))
 297         IFL2=LISTSM(JETTYP(2))
 298         AMQ1=AMASS(IFL1)
 299         AMQ2=AMASS(IFL2)
 300         CALL TWOKIN(0.,0.,AMQ1,AMQ2)
 301       ENDIF
 302 C
 303 C            Set PBEAM and PJETS
 304 C
 305       PBEAM(1)=(1.-X1)*HALFE
 306       PBEAM(2)=(1.-X2)*HALFE
 307       DO 501 I=1,2
 308         PJETS(3,I)=P(I)*CTH(I)
 309         PJETS(1,I)=PT(I)*COS(PHI(I))
 310         PJETS(2,I)=PT(I)*SIN(PHI(I))
 311         IF(KEYS(1)) THEN
 312           IDJETS(I)=LISTJ(JETTYP(I))
 313         ELSEIF(KEYS(5).OR.(KEYS(10).AND.GOMSSM)) THEN
 314           IDJETS(I)=LISTSS(JETTYP(I))
 315         ELSEIF(KEYS(6)) THEN
 316           IDJETS(I)=LISTW(JETTYP(I))
 317         ELSEIF(KEYS(8)) THEN
 318           IDJETS(1)=IFL1
 319           IDJETS(2)=IFL2
 320         ELSEIF(KEYS(10)) THEN
 321           IDJETS(I)=LISTSM(JETTYP(I))
 322         ENDIF
 323         PJETS(5,I)=AMASS(IDJETS(I))
 324         PJETS(4,I)=SQRT(P(I)**2+PJETS(5,I)**2)
 325   501 CONTINUE
 326 C
 327 C          Set PINITS
 328 C
 329       DO 600 I=1,2
 330       IDINIT(I)=LISTJ(INITYP(I))
 331       PINITS(5,I)=AMASS(IDINIT(I))
 332       PPLUS=X(I)*ECM
 333       PMINUS=PINITS(5,I)**2/PPLUS
 334       PINITS(4,I)=.5*(PPLUS+PMINUS)
 335       PINITS(3,I)=.5*(PPLUS-PMINUS)*(3-2*I)
 336       PINITS(2,I)=0.
 337       PINITS(1,I)=0.
 338 600   CONTINUE
 339 C          Calculate PINITS exactly.
 340       PSUM3=PJETS(3,1)+PJETS(3,2)
 341       PSUM4=PJETS(4,1)+PJETS(4,2)
 342       IF(PSUM3.GT.0.) THEN
 343         PPL=PSUM4+PSUM3
 344         PMN=SHAT/PPL
 345       ELSE
 346         PMN=PSUM4-PSUM3
 347         PPL=SHAT/PMN
 348       ENDIF
 349       SQ1=PINITS(5,1)**2
 350       SQ2=PINITS(5,2)**2
 351       ROOT=SQRT((PPL*PMN-SQ1-SQ2)**2-4.*SQ1*SQ2)
 352       P1PL=(PPL*PMN+SQ1-SQ2+ROOT)/(2.*PMN)
 353       P1MN=SQ1/P1PL
 354       P2MN=(PPL*PMN+SQ2-SQ1+ROOT)/(2.*PPL)
 355       P2PL=SQ2/P2MN
 356       PINITS(4,1)=.5*(P1PL+P1MN)
 357       PINITS(3,1)=.5*(P1PL-P1MN)
 358       PINITS(4,2)=.5*(P2PL+P2MN)
 359       PINITS(3,2)=.5*(P2PL-P2MN)
 360       RETURN
 361 C
 362 C          Error
 363 C
 364 910   CALL PRTEVT(0)
 365       WRITE(ITLIS,1000) NREJ
 366  1000 FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE AN',
 367      $' EVENT. CHECK LIMITS OR INCREASE NTRIES.')
 368       STOP 99
 369       END