ISAJET/code/frgmnt.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE FRGMNT
   3 C
   4 C          Control jet fragmentation.  Boost to frames defined in
   5 C          EVOLVE and call JETGEN.
   6 C
   7 C          EVOLVE initializes /PJSET/ as follows--
   8 C          1      - 2            = PINITS (except for E+E-)
   9 C          N0W    - N0W          = QWJET  (for DRELLYAN, NJET=3)
  10 C          N0JETS - N0JETS+NJET  = PJETS
  11 C          N0PAIR - N0PAIR+NPAIR = PPAIR  (for WPAIR)
  12 C
  13 #if defined(CERNLIB_IMPNONE)
  14       IMPLICIT NONE
  15 #endif
  16 #include "isajet/itapes.inc"
  17 #include "isajet/primar.inc"
  18 #include "isajet/jetpar.inc"
  19 #include "isajet/pjets.inc"
  20 #include "isajet/pinits.inc"
  21 #include "isajet/partcl.inc"
  22 #include "isajet/const.inc"
  23 #include "isajet/jetset.inc"
  24 #include "isajet/jwork.inc"
  25 #include "isajet/keys.inc"
  26 #include "isajet/q1q2.inc"
  27 #include "isajet/frame.inc"
  28 C
  29       REAL PSUM(5),PALLJ(5),P12(5),PIN(5,2),PWREST(5),PADD(5)
  30       REAL POLD(5),PNEW(5)
  31       REAL PINPL,PINMN,BP,PT2AVE,PTADD,RANF,PHIADD,PALLPL,PALLMN
  32       REAL PALLX,PALLY
  33       INTEGER K,J,JJET,NZERO,IB,NPTCL1,NPTCL2,IFAIL,JET,NPJET1,NPLV1
  34       INTEGER NPJET3,IP,NP1,NP2,NFIRST,IP1,IFR,NLJ
  35       DOUBLE PRECISION DSUM(5),DPASS(5)
  36 C
  37 C          Initialize
  38       DO 100 K=1,5
  39 100   DSUM(K)=0.
  40       NLJ=NJET
  41       IF(KEYS(3)) NLJ=NJET+1
  42       DO 101 J=1,NLJ
  43         JJET=N0JETS+J-1
  44         IF(JJET.EQ.N0W) GOTO 101
  45         CALL DBLVEC(PJSET(1,JJET),DPASS)
  46         DO 102 K=1,4
  47 102     DSUM(K)=DSUM(K)+DPASS(K)
  48 101   CONTINUE
  49       DSUM(5)=DSQRT(DSUM(4)**2-DSUM(1)**2-DSUM(2)**2-DSUM(3)**2)
  50       DO 103 K=1,5
  51 103   PALLJ(K)=DSUM(K)
  52 C
  53       NZERO=NPTCL+1
  54 C
  55 C          Fragment partons from initial state shower
  56 C
  57       IF(.NOT.KEYS(2)) THEN
  58         DO 110 J=1,NJSET
  59           IF(JDCAY(J).EQ.JPACK*J+J) THEN
  60             IB=JORIG(J)/JPACK-10
  61             DO 120 K=1,5
  62 120         PIN(K,IB)=PJSET(K,J)
  63           ENDIF
  64 110     CONTINUE
  65 C
  66         CALL FRGJET(11)
  67         CALL FRGJET(12)
  68 C
  69         NPTCL1=NPTCL+1
  70         NPTCL2=NPTCL1+1
  71         IF(NPTCL1.GT.MXPTCL) GO TO 9999
  72         PINPL=.5*(PIN(4,1)+PIN(3,1)+PIN(4,2)+PIN(3,2))
  73         PINMN=.5*(PIN(4,1)-PIN(3,1)+PIN(4,2)-PIN(3,2))
  74         PPTCL(1,NPTCL1)=0.
  75         PPTCL(2,NPTCL1)=0.
  76         PPTCL(3,NPTCL1)=HALFE-PINPL
  77         PPTCL(4,NPTCL1)=HALFE-PINPL
  78         PPTCL(5,NPTCL1)=0.
  79         PPTCL(1,NPTCL2)=0.
  80         PPTCL(2,NPTCL2)=0.
  81         PPTCL(3,NPTCL2)=-(HALFE-PINMN)
  82         PPTCL(4,NPTCL2)=HALFE-PINMN
  83         PPTCL(5,NPTCL2)=0.
  84         DO 130 K=1,4
  85 130     PSUM(K)=-PALLJ(K)
  86         PSUM(4)=PSUM(4)+ECM
  87         PSUM(5)=PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2
  88         IF(PSUM(5).GE.0.) THEN
  89           PSUM(5)=SQRT(PSUM(5))
  90           CALL RESCAL(NZERO,NPTCL2,PSUM,IFAIL)
  91         ENDIF
  92 C
  93         DO 140 K=1,4
  94 140     PBEAMS(K)=PPTCL(K,NPTCL1)+PPTCL(K,NPTCL2)
  95         PBEAMS(5)=SQRT(PBEAMS(4)**2-PBEAMS(1)**2-PBEAMS(2)**2
  96      $  -PBEAMS(3)**2)
  97       ENDIF
  98 C
  99 C          Boost partons from final jets with -FRAME
 100 C
 101 200   DO 210 J=1,NJSET
 102         JET=JORIG(J)/JPACK
 103         IF ( JET.EQ.0 ) THEN
 104           IFR=1
 105         ELSE
 106           IF(JET.GT.10) GO TO 210
 107           IF(KEYS(6)) THEN
 108             IF(IDJETS(JET).EQ.10) GO TO 210
 109           ENDIF
 110           IFR=IFRAME(JET)
 111         ENDIF
 112 C
 113 C          Do this boost in double precision for 32-bit machines
 114         CALL DBOOST(-1,FRAME(1,IFR),PJSET(1,J))
 115 210   CONTINUE
 116 C
 117 C          Fragment partons from final jets
 118 C
 119       NPJET1=NPTCL+1
 120       DO 220 K=1,4
 121 220   PSUM(K)=0
 122 C
 123 C          Conserve mass of 1+2 for DRELLYAN (automatic for WPAIR)
 124 C
 125       IF(KEYS(3)) THEN
 126         CALL FRGJET(1)
 127         CALL FRGJET(2)
 128         IF(STDDY) THEN
 129           DO 230 K=1,4
 130             PSUM(K)=PJSET(K,3)+PJSET(K,4)
 131   230     CONTINUE
 132         ELSE
 133           DO 240 K=1,4
 134             PSUM(K)=PJSET(K,N0W+1)+PJSET(K,N0W+2)
 135   240     CONTINUE
 136         ENDIF
 137         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
 138         NPLV1=NPTCL
 139         CALL RESCAL(NPJET1,NPLV1,PSUM,IFAIL)
 140 C          EXTRADIM has only jet3 + graviton
 141       ELSEIF(KEYS(11)) THEN
 142         CALL FRGJET(3)
 143         CALL FRGJET(0)
 144         NPLV1=NPTCL
 145         DO 241 K=1,4
 146           PSUM(K)=PJSET(K,3)+PJSET(K,4)
 147 241     CONTINUE
 148         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
 149         CALL RESCAL (NPJET1,NPLV1,PSUM,IFAIL)
 150       ELSE
 151 C          All other processes
 152         DO 242 J=1,NJET
 153           JJET=N0JETS+J-1
 154           CALL FRGJET(J)
 155           DO 243 K=1,4
 156   243     PSUM(K)=PSUM(K)+PJSET(K,JJET)
 157   242   CONTINUE
 158         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
 159         NPLV1=NPTCL
 160         CALL RESCAL(NPJET1,NPLV1,PSUM,IFAIL)
 161       ENDIF
 162 C
 163 C           Add extra jets for DRELLYAN
 164       IF(KEYS(3).AND..NOT.STDDY) THEN
 165         NPJET3=NPTCL+1
 166         DO 245 J=3,NJET
 167 245     CALL FRGJET(J)
 168         NPTCL1=NPTCL+1
 169         IF(NPTCL1.GT.MXPTCL) GO TO 9999
 170         DO 250 K=1,4
 171           PPTCL(K,NPTCL1)=PJSET(K,N0W)
 172 250     PSUM(K)=PJSET(K,N0W)
 173         PPTCL(5,NPTCL1)=PJSET(5,N0W)
 174         DO 246 J=3,NJET
 175           JJET=N0JETS+J-3
 176           DO 246 K=1,4
 177             PSUM(K)=PSUM(K)+PJSET(K,JJET)
 178   246   CONTINUE
 179         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
 180         CALL RESCAL(NPJET3,NPTCL1,PSUM,IFAIL)
 181         DO 260 K=1,5
 182 260     PWREST(K)=PPTCL(K,NPTCL1)
 183       ENDIF
 184 C
 185 C          Boost partons back to PP COM
 186 C
 187       DO 300 J=1,NJSET
 188         JET=JORIG(J)/JPACK
 189         IF ( JET.EQ.0 ) THEN
 190           IFR=1
 191         ELSE
 192           IF(JET.GT.10) GO TO 300
 193           IF(KEYS(6)) THEN
 194             IF(IDJETS(JET).EQ.10) GO TO 300
 195           ENDIF
 196           IFR=IFRAME(JET)
 197         ENDIF
 198         BP=0.
 199         DO 305 K=1,3
 200 305     BP=BP+FRAME(K,IFR)*PJSET(K,J)
 201         BP=BP/FRAME(5,IFR)
 202         DO 310 K=1,3
 203 310     PJSET(K,J)=PJSET(K,J)+FRAME(K,IFR)*PJSET(4,J)/FRAME(5,IFR)
 204      $  +FRAME(K,IFR)*BP/(FRAME(4,IFR)+FRAME(5,IFR))
 205         PJSET(4,J)=FRAME(4,IFR)*PJSET(4,J)/FRAME(5,IFR)+BP
 206 300   CONTINUE
 207 C
 208 C          Reset FRAME to boost hadrons to PP COM
 209 C
 210       IF(KEYS(1).OR.KEYS(2).OR.(KEYS(3).AND.NJET.EQ.2).OR.KEYS(5)
 211      $.OR.(KEYS(7).AND.NPAIR.EQ.0).OR.KEYS(8)) THEN
 212         DO 410 K=1,5
 213           FRAME(K,1)=PALLJ(K)
 214 410     CONTINUE
 215       ELSEIF(KEYS(3).AND.NJET.GT.2) THEN
 216         DO 420 K=1,5
 217 420     FRAME(K,1)=PALLJ(K)
 218         BP=0.
 219         DO 430 K=1,3
 220 430     BP=BP+FRAME(K,1)*PWREST(K)
 221         BP=BP/FRAME(5,1)
 222         DO 440 K=1,3
 223           FRAME(K,2)=PWREST(K)+FRAME(K,1)*PWREST(4)/FRAME(5,1)
 224      $    +FRAME(K,1)*BP/(FRAME(4,1)+FRAME(5,1))
 225 440     CONTINUE
 226         FRAME(4,2)=FRAME(4,1)*PWREST(4)/FRAME(5,1)+BP
 227       ENDIF
 228 C
 229 C          Boost hadrons back to PP COM
 230 C
 231       DO 500 IP=NZERO,NPTCL
 232         JET=IABS(IORIG(IP))/IPACK
 233         IF(JET.GT.10) GO TO 500
 234         IF(KEYS(6)) THEN
 235           IF(IDJETS(JET).EQ.10) GO TO 500
 236         ENDIF
 237         IF(JET.EQ.0) THEN
 238           IFR=1
 239         ELSE
 240           IFR=IFRAME(JET)
 241         ENDIF
 242         BP=0.
 243         DO 510 K=1,3
 244 510     BP=BP+FRAME(K,IFR)*PPTCL(K,IP)
 245         BP=BP/FRAME(5,IFR)
 246         DO 520 K=1,3
 247 520     PPTCL(K,IP)=PPTCL(K,IP)+FRAME(K,IFR)*PPTCL(4,IP)/FRAME(5,IFR)
 248      $  +FRAME(K,IFR)*BP/(FRAME(4,IFR)+FRAME(5,IFR))
 249         PPTCL(4,IP)=FRAME(4,IFR)*PPTCL(4,IP)/FRAME(5,IFR)+BP
 250 500   CONTINUE
 251 C
 252 C          Add intrinsic PT
 253 C
 254       IF(.NOT.KEYS(2)) THEN
 255         PT2AVE=.1*SQRT(QSQ)
 256         PTADD=SQRT(-PT2AVE*ALOG(RANF()))
 257         PHIADD=2.*PI*RANF()
 258         PADD(1)=2.*PTADD*COS(PHIADD)
 259         PADD(2)=2.*PTADD*SIN(PHIADD)
 260 C          Must use large and small components carefully to calculate
 261 C          mass on 32-bit machines.
 262         PALLPL=0.
 263         PALLMN=0.
 264         PALLX=0.
 265         PALLY=0.
 266         DO 525 IP=NZERO,NPTCL
 267           PALLX=PALLX+PPTCL(1,IP)
 268           PALLY=PALLY+PPTCL(2,IP)
 269           IF(PPTCL(3,IP).GT.0.) THEN
 270             PALLPL=PALLPL+(PPTCL(4,IP)+PPTCL(3,IP))
 271             PALLMN=PALLMN+(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(5,IP)**2)
 272      $      /(PPTCL(4,IP)+PPTCL(3,IP))
 273           ELSE
 274             PALLMN=PALLMN+(PPTCL(4,IP)-PPTCL(3,IP))
 275             PALLPL=PALLPL+(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(5,IP)**2)
 276      $      /(PPTCL(4,IP)-PPTCL(3,IP))
 277           ENDIF
 278 525     CONTINUE
 279         POLD(1)=PALLX
 280         POLD(2)=PALLY
 281         POLD(3)=.5*(PALLPL-PALLMN)
 282         POLD(4)=.5*(PALLPL+PALLMN)
 283         POLD(5)=SQRT(PALLPL*PALLMN-PALLX**2-PALLY**2)
 284         PNEW(1)=PADD(1)+POLD(1)
 285         PNEW(2)=PADD(2)+POLD(2)
 286         PNEW(3)=POLD(3)
 287         PNEW(4)=SQRT(PNEW(1)**2+PNEW(2)**2+PNEW(3)**2+POLD(5)**2)
 288         PNEW(5)=POLD(5)
 289 C
 290         DO 530 IP=NZERO,NPTCL
 291           BP=0.
 292           DO 531 K=1,3
 293 531       BP=BP+POLD(K)*PPTCL(K,IP)
 294           BP=BP/POLD(5)
 295           DO 532 K=1,3
 296 532       PPTCL(K,IP)=PPTCL(K,IP)-POLD(K)*PPTCL(4,IP)/POLD(5)
 297      $    +POLD(K)*BP/(POLD(4)+POLD(5))
 298           PPTCL(4,IP)=PPTCL(4,IP)*POLD(4)/POLD(5)-BP
 299 C
 300           BP=0.
 301           DO 533 K=1,3
 302 533       BP=BP+PNEW(K)*PPTCL(K,IP)
 303           BP=BP/PNEW(5)
 304           DO 534 K=1,3
 305 534       PPTCL(K,IP)=PPTCL(K,IP)+PNEW(K)*PPTCL(4,IP)/PNEW(5)
 306      $    +PNEW(K)*BP/(PNEW(4)+PNEW(5))
 307           PPTCL(4,IP)=PPTCL(4,IP)*PNEW(4)/PNEW(5)+BP
 308 530     CONTINUE
 309 C
 310 C            Add opposite PT to beam jets
 311         DO 541 K=1,4
 312 541     PBEAMS(K)=-PNEW(K)
 313         PBEAMS(4)=PBEAMS(4)+ECM
 314         PBEAMS(5)=PBEAMS(4)**2-PBEAMS(1)**2-PBEAMS(2)**2 -PBEAMS(3)**2
 315         IF ( PBEAMS(5).GT.0 ) THEN
 316           PBEAMS(5)=SQRT(PBEAMS(5))
 317         ELSE
 318           PBEAMS(4)=SQRT(PBEAMS(4)**2-PBEAMS(5)+4.)
 319           PBEAMS(5)=2.
 320         ENDIF
 321       ENDIF
 322 C
 323 C          Decay hadrons
 324 C
 325       NP1=NZERO
 326 600   NP2=NPTCL
 327       DO 610 IP=NP1,NP2
 328         NFIRST=NPTCL+1
 329         JET=IABS(IORIG(IP))/IPACK
 330         CALL DECAY(IP)
 331         DO 620 IP1=NFIRST,NPTCL
 332 620     IORIG(IP1)=ISIGN(IABS(IORIG(IP1))+IPACK*JET,IORIG(IP1))
 333 610   CONTINUE
 334       NP1=NP2+1
 335       IF(NP1.LE.NPTCL) GO TO 600
 336       RETURN
 337 C
 338 C          Error
 339 C
 340 9999  CALL PRTEVT(0)
 341       WRITE(ITLIS,9998) NPTCL
 342 9998  FORMAT(//' ERROR IN FRGMNT ... NPTCL > ',I6)
 343       RETURN
 344       END