ISAJET/code/frgjet.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE FRGJET(JET)
   3 C
   4 C          Hadronize all partons in /JETSET/ corresponding to jet JET.
   5 C
   6 #if defined(CERNLIB_IMPNONE)
   7       IMPLICIT NONE
   8 #endif
   9 #include "isajet/itapes.inc"
  10 #include "isajet/primar.inc"
  11 #include "isajet/jetpar.inc"
  12 #include "isajet/pjets.inc"
  13 #include "isajet/pinits.inc"
  14 #include "isajet/partcl.inc"
  15 #include "isajet/const.inc"
  16 #include "isajet/jetset.inc"
  17 #include "isajet/jwork.inc"
  18 #include "isajet/keys.inc"
  19 #include "isajet/q1q2.inc"
  20 #include "isajet/frame.inc"
  21 C
  22       REAL ROT(3,3),POLD(5),PNEW(5),PSUM(5)
  23       REAL CPHI,SPHI,AMSUM,ESUM,PJ,CTHJ,STHJ,PTJ
  24       INTEGER K,K1,K2,IP,NPLV1,IFAIL,NBEGIN,JET,NFRAG,NFRGMX,JETJ,
  25      $JTABS,NFIRST,J
  26 C
  27       DATA PSUM/5*0./
  28 C
  29 C          NFRAG counter protects against possible infinite loop.
  30 C
  31       NFRAG=0
  32       NFRGMX=10*MXJSET
  33 201   NBEGIN=NPTCL+1
  34       NFRAG=NFRAG+1
  35 C
  36 C          Loop over partons
  37 C
  38       ESUM=0.
  39       DO 220 J=1,NJSET
  40         IF(JDCAY(J).NE.0) GO TO 220
  41         JETJ=JORIG(J)/JPACK
  42         IF(JETJ.NE.JET) GO TO 220
  43         ESUM=ESUM+PJSET(4,J)
  44 C
  45 C          Generate Field-Feynman jet for each quark or gluon, or...
  46 C
  47         JTABS = IABS(JTYPE(J))
  48         IF(JTABS.LT.10) THEN
  49           NFIRST=NPTCL+1
  50           CALL JETGEN(J)
  51           IF(NPTCL.LT.NFIRST) GO TO 220
  52 C
  53 C          Rotate hadrons to parton direction
  54 C
  55           PTJ=PJSET(1,J)**2+PJSET(2,J)**2
  56           PJ=SQRT(PTJ+PJSET(3,J)**2)
  57           PTJ=SQRT(PTJ)
  58 C          Following is to fix occasional bug on 32-bit machines
  59           IF(PJ.GT.0.) THEN
  60             CTHJ=PJSET(3,J)/PJ
  61             STHJ=PTJ/PJ
  62           ELSE
  63             CTHJ=1.
  64             STHJ=0.
  65           ENDIF
  66           IF(PTJ.GT.0.) THEN
  67             CPHI=PJSET(1,J)/PTJ
  68             SPHI=PJSET(2,J)/PTJ
  69           ELSE
  70             CPHI=SIGN(1.,PJSET(3,J))
  71             SPHI=0.
  72           ENDIF
  73           ROT(1,1)=CPHI*CTHJ
  74           ROT(2,1)=SPHI*CTHJ
  75           ROT(3,1)=-STHJ
  76           ROT(1,2)=-SPHI
  77           ROT(2,2)=CPHI
  78           ROT(3,2)=0.
  79           ROT(1,3)=CPHI*STHJ
  80           ROT(2,3)=SPHI*STHJ
  81           ROT(3,3)=CTHJ
  82           DO 230 IP=NFIRST,NPTCL
  83             DO 235 K=1,3
  84               POLD(K)=PPTCL(K,IP)
  85               PPTCL(K,IP)=0
  86 235         CONTINUE
  87             DO 240 K1=1,3
  88             DO 240 K2=1,3
  89 240         PPTCL(K1,IP)=PPTCL(K1,IP)+ROT(K1,K2)*POLD(K2)
  90 230       CONTINUE
  91 C
  92 C          ... hadronize all other partons with delta function.
  93 C
  94         ELSE
  95           IF((IABS(JTYPE(J)).EQ.80.OR.IABS(JTYPE(J)).EQ.90).AND.
  96      $    .NOT.KEYS(2).AND..NOT.KEYS(12)) GO TO 210
  97           IF(NPTCL.GE.MXPTCL) GO TO 9999
  98           NPTCL=NPTCL+1
  99           DO 255 K=1,5
 100             PPTCL(K,NPTCL)=PJSET(K,J)
 101 255       CONTINUE
 102           IORIG(NPTCL)=-J
 103           IDENT(NPTCL)=JTYPE(J)
 104           IDCAY(NPTCL)=0
 105         ENDIF
 106 220   CONTINUE
 107 C
 108 C          Sum masses and insert jet label
 109 C
 110       AMSUM=0.
 111       DO 260 IP=NBEGIN,NPTCL
 112         AMSUM=AMSUM+PPTCL(5,IP)
 113         IORIG(IP)=ISIGN(IABS(IORIG(IP))+IPACK*JET,IORIG(IP))
 114 260   CONTINUE
 115 C
 116 C          Require sum of masses less than jet energy.
 117 C
 118       IF(AMSUM.GT.ESUM.AND.NBEGIN.NE.NPTCL.AND.NFRAG.LT.NFRGMX) THEN
 119         NPTCL=NBEGIN-1
 120         GO TO 201
 121       ENDIF
 122 C
 123 C          For WPAIR events rescale jet to W mass.
 124 C
 125       IF((KEYS(6).OR.KEYS(7).OR.KEYS(9).OR.KEYS(10)).AND.JET.LT.10)
 126      $ THEN
 127         IF(IABS(JTYPE(JET+N0JETS-1)).LT.80) RETURN
 128         IF(AMSUM.GE.PJSET(5,JET+N0JETS-1)) THEN
 129           IF(NFRAG.GT.NFRGMX) RETURN
 130           NPTCL=NBEGIN-1
 131           GO TO 201
 132         ENDIF
 133         PSUM(4)=PJSET(5,JET+N0JETS-1)
 134         PSUM(5)=PSUM(4)
 135         NPLV1=NPTCL
 136         CALL RESCAL(NBEGIN,NPLV1,PSUM,IFAIL)
 137       ENDIF
 138 C
 139 210   RETURN
 140 C
 141 C          Error
 142 C
 143 9999  CALL PRTEVT(0)
 144       WRITE(ITLIS,9998) NPTCL
 145 9998  FORMAT(//' ERROR IN FRGJET ... NPTCL > ',I6)
 146       RETURN
 147       END