ISAJET/code/mbias.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE MBIAS
   3 C
   4 C          Generate minbias event or beam jets for high-pt event using
   5 C          parameters set in MBSET:
   6 C
   7 C          (1) Select number NPOM of cut pomerons -- cf cut Reggeon
   8 C              field theory of Abramovskii, Kanchelli, and Gribov.
   9 C          (2) Generate xf for leading baryons including 1/(1-xf)
  10 C              diffractive term and guessed NPOM dependence,
  11 C              F(XF)=(1-XF)**(A+B/NPOM)
  12 C          (3) Select xf for each half of each Pomeron. then fragment
  13 C              each half Pomeron into mesons and baryons independently
  14 C              in the Pomeron-Pomeron center of mass. This avoids
  15 C              making xf=0 a singular point.
  16 C
  17 C          Note that multiple cut Pomerons give approximate KNO scaling.
  18 C          The only short-range correlations are from resonances.
  19 C
  20 C          Ver. 7.09: Add traps on free loops and IMPLICIT NONE.
  21 C
  22 #if defined(CERNLIB_IMPNONE)
  23       IMPLICIT NONE
  24 #endif
  25 #include "isajet/itapes.inc"
  26 #include "isajet/keys.inc"
  27 #include "isajet/mbgen.inc"
  28 #include "isajet/primar.inc"
  29 #include "isajet/jetpar.inc"
  30 #include "isajet/const.inc"
  31 #include "isajet/partcl.inc"
  32 #include "isajet/mbpar.inc"
  33 C
  34       DIMENSION IFL(3),IFLEXC(2),PXEXC(2),PYEXC(2),SIGN(2)
  35       DIMENSION PSUM(5)
  36       DIMENSION LDIFFR(2)
  37       LOGICAL LDIFFR
  38       REAL RANF,AMASS
  39       REAL RND,XX,XSUM,P0,PPOM,PXEXC,PSUM,SIGN,PYEXC,GAM,BETA,X,
  40      $AM,PPLUS,EPSDIF,PEND0,TRY,PX,PY,PX2,PY2,QMINUS,PZ,QPLUS,PT1,
  41      $PHI1,XBGEN,PT2,PHI2,PX1,PY1,AMT2
  42       INTEGER ID1,ID2,IFL1,IFL2,IMOD1,IMOD2,ITWIST,IPOM,LOOP,NFIRST,
  43      $ID3,IFLEXC,IFL,I,NP2,IP,NP1,IFAIL,NPTLV1,IDHAD,IB,NEW,JSPIN,
  44      $INDEX,NBEGIN,IPASS,MXPASS,N,IDIFF,IPASSB,IFLNEW,ISWAP
  45       DATA SIGN/1.,-1./,PEND0/.14/
  46       DATA PSUM/5*0./
  47       DATA MXPASS/200/
  48 C
  49 C          Start
  50       NBEGIN=NPTCL+1
  51       IPASS=1
  52       IPASSB=1
  53 C
  54 C          Select number of cut Pomerons.
  55 C
  56 1     CONTINUE
  57       TRY=RANF()
  58       DO 10 N=MNPOM,MXPOM
  59         NPOM=N
  60         IF(POMGEN(N).GT.TRY) GO TO 20
  61 10    CONTINUE
  62 20    CONTINUE
  63 C
  64 C          Decide if diffractive event
  65       IF(RANF().LT.PDIFFR) THEN
  66         IDIFF=INT(1.99999*RANF())+1
  67         LDIFFR(IDIFF)=.TRUE.
  68         LDIFFR(3-IDIFF)=.FALSE.
  69       ELSE
  70         LDIFFR(1)=.FALSE.
  71         LDIFFR(2)=.FALSE.
  72       ENDIF
  73 C
  74 C          Generate leading baryons.
  75 C
  76       DO 100 IB=1,2
  77         PPLUS=2.*PBEAM(IB)
  78 C
  79 C          Special treatment for diffractive beam.
  80         IF(LDIFFR(IB)) THEN
  81           IDHAD=IDIN(IB)
  82           AM=AMASS(IDHAD)
  83           CALL FLAVOR(IDIN(IB),IFL(1),IFL(2),IFL(3),JSPIN,INDEX)
  84           NEW=INT(3.*RANF())+1
  85           IFLEXC(1)=+IFL(NEW)
  86           IFLEXC(2)=-IFL(NEW)
  87           EPSDIF=2./SCM
  88           DXBARY(IB)=EPSDIF**RANF()
  89           XBARY(IB)=1.-DXBARY(IB)
  90           GO TO 115
  91         ENDIF
  92 C
  93 C          If not diffractive, construct IDENT of leading baryon
  94         CALL FLAVOR(IDIN(IB),IFL(1),IFL(2),IFL(3),JSPIN,INDEX)
  95         NEW=INT(3.*RANF())+1
  96         IFLNEW=ISIGN(INT(RANF()/PUD0)+1,IDIN(IB))
  97         IFLEXC(1)=IFL(NEW)
  98         IFLEXC(2)=-IFLNEW
  99         IFL(NEW)=IFLNEW
 100         IF(IABS(IFL(1)).GT.IABS(IFL(2))) THEN
 101           ISWAP=IFL(1)
 102           IFL(1)=IFL(2)
 103           IFL(2)=ISWAP
 104         ENDIF
 105         IF(IABS(IFL(2)).GT.IABS(IFL(3))) THEN
 106           ISWAP=IFL(2)
 107           IFL(2)=IFL(3)
 108           IFL(3)=ISWAP
 109         ENDIF
 110         IF(IABS(IFL(1)).GT.IABS(IFL(2))) THEN
 111           ISWAP=IFL(1)
 112           IFL(1)=IFL(2)
 113           IFL(2)=ISWAP
 114         ENDIF
 115         JSPIN=1
 116         IF(IFL(1).EQ.IFL(2).AND.IFL(2).EQ.IFL(3)) THEN
 117           JSPIN=1
 118         ELSE
 119           JSPIN=INT(RANF()+PJSPN)
 120         ENDIF
 121         IF(JSPIN.EQ.0.AND.IFL(1).NE.IFL(2).AND.IFL(2).NE.IFL(3)) THEN
 122           IF(RANF().GT.PISPN) THEN
 123             ISWAP=IFL(1)
 124             IFL(1)=IFL(2)
 125             IFL(2)=ISWAP
 126           ENDIF
 127         ENDIF
 128         IDHAD=1000*IABS(IFL(1))+100*IABS(IFL(2))+10*IABS(IFL(3))+JSPIN
 129         IDHAD=ISIGN(IDHAD,IDIN(IB))
 130         AM=AMASS(IDHAD)
 131 C
 132 C          Select xf for nondiffractive baryon, flat for NPOM=1 and
 133 C          like mesons for NPOM=infinity.
 134 110     XBGEN=XGEN0(2)*(1.-1./NPOM)
 135         DXBARY(IB)=RANF()**(1./(XBGEN+1.))
 136         XBARY(IB)=1.-DXBARY(IB)
 137 C
 138 C          Select transverse momentum of baryon
 139 115     CALL GETPT(PT1,SIGQT0)
 140         PHI1=2.*PI*RANF()
 141         PX1=PT1*COS(PHI1)
 142         PY1=PT1*SIN(PHI1)
 143         PXEXC(1)=PX1
 144         PYEXC(1)=PY1
 145         CALL GETPT(PT2,SIGQT0)
 146         PHI2=2.*PI*RANF()
 147         PX2=PT2*COS(PHI2)
 148         PY2=PT2*SIN(PHI2)
 149         PXEXC(2)=PX2
 150         PYEXC(2)=PY2
 151         PX=-PX1-PX2
 152         PY=-PY1-PY2
 153         AMT2=PX**2+PY**2+AM**2
 154 C
 155         QPLUS=XBARY(IB)*PPLUS
 156         QPLUS=AMAX1(QPLUS,1.E-6)
 157         QMINUS=AMT2/QPLUS
 158         PZ=.5*(QPLUS-QMINUS)
 159         P0=.5*(QPLUS+QMINUS)
 160 C
 161 C          Add baryon to /PARTCL/ if PZ>0.
 162         IF(NPTCL.GE.MXPTCL) GO TO 9999
 163         IF(PZ.GE.0.) THEN
 164           NPTCL=NPTCL+1
 165           PPTCL(1,NPTCL)=PX
 166           PPTCL(2,NPTCL)=PY
 167           PPTCL(3,NPTCL)=PZ*SIGN(IB)
 168           PPTCL(4,NPTCL)=P0
 169           PPTCL(5,NPTCL)=AM
 170           IORIG(NPTCL)=0
 171           IDCAY(NPTCL)=0
 172           IDENT(NPTCL)=IDHAD
 173         ELSE
 174           IPASSB=IPASSB+1
 175           IF(IPASSB.LT.MXPASS) GO TO 110
 176 C          Just give up if it fails MXPASS times
 177           WRITE(ITLIS,998)
 178 998       FORMAT(//5X,'ERROR IN MBIAS ... COULD NOT MAKE BARYON')
 179           XBARY(IB)=0.
 180           DXBARY(IB)=1.
 181         ENDIF
 182 C
 183 C          Having accepted baryon, set up XPOM array for cut Pomerons,
 184 C          rescaling to 1.-XBARY(IB).
 185         XSUM=0.
 186         DO 120 N=1,NPOM
 187           XX=RANF()
 188           XPOM(N,IB)=XX
 189           XSUM=XSUM+XX
 190 120     CONTINUE
 191         XSUM=1./XSUM
 192         DO 130 N=1,NPOM
 193           XPOM(N,IB)=XSUM*XPOM(N,IB)*DXBARY(IB)
 194 130     CONTINUE
 195 100   CONTINUE
 196 C
 197 C          Fragment each Pomeron into mesons and baryon pairs in the
 198 C          Pomeron-Pomeron center of mass.
 199 C
 200       DO 1000 IB=1,2
 201         DO 2000 IPOM=1,NPOM
 202           PPOM=SQRT(PBEAM(1)*XPOM(IPOM,1)*PBEAM(2)*XPOM(IPOM,2))
 203           PPLUS=2.*PPOM
 204           NFIRST=NPTCL+1
 205           LOOP=0
 206 C
 207 200       CONTINUE
 208           ITWIST=INT(1.99999*RANF())+1
 209           LOOP=LOOP+1
 210 C
 211 C          Select new quark or diquark. Old diquark implies new quark.
 212 C          Old quark implies new diquark with probability PBARY0.
 213           IFL1=IFLEXC(ITWIST)
 214           IF(MOD(IFL1,100).EQ.0) THEN
 215             IFL2=ISIGN(INT(RANF()/PUD0)+1,+IFL1)
 216           ELSEIF(RANF().GT.PBARY0) THEN
 217             IFL2=ISIGN(INT(RANF()/PUD0)+1,-IFL1)
 218           ELSE
 219             ID1=INT(RANF()/PUD0)+1
 220             ID2=INT(RANF()/PUD0)+1
 221             IF(IABS(ID1).GT.IABS(ID2)) THEN
 222               ISWAP=ID1
 223               ID1=ID2
 224               ID2=ISWAP
 225             ENDIF
 226             IFL2=ISIGN(1000*ID1+100*ID2,IFL1)
 227           ENDIF
 228           IFLEXC(ITWIST)=-IFL2
 229 C          Construct meson from quark+antiquark. Else, construct baryon
 230 C          IDENT from quark+diquark.
 231           IMOD1=MOD(IFL1,100)
 232           IMOD2=MOD(IFL2,100)
 233           IF(IMOD1.NE.0.AND.IMOD2.NE.0) THEN
 234             JSPIN=INT(RANF()+PJSPN)
 235             ID1=IFL1
 236             ID2=IFL2
 237             IF(ID1+ID2.EQ.0) THEN
 238               RND=RANF()
 239               ID1=IABS(ID1)
 240               ID1=INT(PMIX01(ID1,JSPIN+1)+RND)
 241      $        +INT(PMIX02(ID1,JSPIN+1)+RND)+1
 242               ID2=-ID1
 243             ELSEIF(IABS(ID1).GT.IABS(ID2)) THEN
 244               ISWAP=ID1
 245               ID1=ID2
 246               ID2=ISWAP
 247             ENDIF
 248             IDHAD=ISIGN(100*IABS(ID1)+10*IABS(ID2)+JSPIN,ID1)
 249           ELSE
 250             IF(IMOD1.EQ.0) THEN
 251               ID3=MOD(IFL1/100,10)
 252               ID2=IFL1/1000
 253               ID1=IFL2
 254             ELSE
 255               ID3=MOD(IFL2/100,10)
 256               ID2=IFL2/1000
 257               ID1=IFL1
 258             ENDIF
 259             IF(IABS(ID1).GT.IABS(ID2)) THEN
 260               ISWAP=ID1
 261               ID1=ID2
 262               ID2=ISWAP
 263             ENDIF
 264             IF(IABS(ID2).GT.IABS(ID3)) THEN
 265               ISWAP=ID2
 266               ID2=ID3
 267               ID3=ISWAP
 268             ENDIF
 269             IF(IABS(ID1).GT.IABS(ID2)) THEN
 270               ISWAP=ID1
 271               ID1=ID2
 272               ID2=ISWAP
 273             ENDIF
 274             IF(ID1.EQ.ID2.AND.ID2.EQ.ID3) THEN
 275               JSPIN=1
 276             ELSE
 277               JSPIN=INT(RANF()+PJSPN)
 278             ENDIF
 279             IF(JSPIN.EQ.0.AND.ID1.NE.ID2.AND.ID2.NE.ID3) THEN
 280               IF(RANF().LT.PISPN) THEN
 281                 ISWAP=ID1
 282                 ID1=ID2
 283                 ID2=ISWAP
 284               ENDIF
 285             ENDIF
 286             IDHAD=1000*IABS(ID1)+100*IABS(ID2)+10*IABS(ID3)+JSPIN
 287             IDHAD=ISIGN(IDHAD,IFL1)
 288           ENDIF
 289 C
 290           AM=AMASS(IDHAD)
 291           PX1=PXEXC(ITWIST)
 292           PY1=PYEXC(ITWIST)
 293           CALL GETPT(PT2,SIGQT0)
 294           PHI2=2.*PI*RANF()
 295           PX2=PT2*COS(PHI2)
 296           PY2=PT2*SIN(PHI2)
 297           PXEXC(ITWIST)=PX2
 298           PYEXC(ITWIST)=PY2
 299           PX=PX1-PX2
 300           PY=PY1-PY2
 301           AMT2=PX**2+PY**2+AM**2
 302 C
 303 C          Select x -- same distribution for all particles.
 304           X=RANF()
 305           IF(RANF().LT.XGEN0(1)) X=1.-X**(1./(XGEN0(2)+1.))
 306           QPLUS=X*PPLUS
 307           QPLUS=AMAX1(QPLUS,1.E-6)
 308           QMINUS=AMT2/QPLUS
 309           P0=.5*(QPLUS+QMINUS)
 310           PZ=.5*(QPLUS-QMINUS)
 311 C
 312 C          Add particle to /PARTCL/ if PZ>0.
 313           IF(NPTCL.GE.MXPTCL) GO TO 9999
 314           IF(PZ.GE.0.) THEN
 315             NPTCL=NPTCL+1
 316             PPTCL(1,NPTCL)=PX
 317             PPTCL(2,NPTCL)=PY
 318             PPTCL(3,NPTCL)=PZ*SIGN(IB)
 319             PPTCL(4,NPTCL)=P0
 320             PPTCL(5,NPTCL)=AM
 321             IORIG(NPTCL)=0
 322             IDCAY(NPTCL)=0
 323             IDENT(NPTCL)=IDHAD
 324           ENDIF
 325 C
 326 C          Continue if sufficient pplus
 327           PPLUS=(1.-X)*PPLUS
 328           IF(PPLUS.GT.PEND0.AND.LOOP.LT.MXPTCL) GO TO 200
 329 C
 330 C          Boost hadrons to lab frame.
 331           IF(NPTCL.LT.NFIRST) GO TO 2000
 332           BETA=(XPOM(IPOM,1)*PBEAM(1)-XPOM(IPOM,2)*PBEAM(2))/(2.*PPOM)
 333           GAM=(XPOM(IPOM,1)*PBEAM(1)+XPOM(IPOM,2)*PBEAM(2))/(2.*PPOM)
 334           DO 400 IP=NFIRST,NPTCL
 335             P0=GAM*PPTCL(4,IP)+BETA*PPTCL(3,IP)
 336             PZ=BETA*PPTCL(4,IP)+GAM*PPTCL(3,IP)
 337             PPTCL(3,IP)=PZ
 338             PPTCL(4,IP)=P0
 339 400       CONTINUE
 340 C
 341 2000    CONTINUE
 342 1000  CONTINUE
 343 C
 344 C          Rescale hadron momenta for correct four-momentum.
 345 C
 346       NPTLV1=NPTCL
 347       IF(KEYS(4)) THEN
 348         PSUM(4)=ECM
 349         PSUM(5)=ECM
 350         CALL RESCAL(NBEGIN,NPTLV1,PSUM,IFAIL)
 351       ELSE
 352         CALL RESCAL(NBEGIN,NPTLV1,PBEAMS,IFAIL)
 353       ENDIF
 354       IF(IFAIL.NE.0.AND.IPASS.LT.MXPASS) THEN
 355         IPASS=IPASS+1
 356         NPTCL=NBEGIN-1
 357         GO TO 1
 358       ENDIF
 359 C
 360 C          Decay hadrons
 361 C
 362       NP1=NBEGIN
 363 500   NP2=NPTCL
 364       DO 510 I=NP1,NP2
 365       CALL DECAY(I)
 366 510   CONTINUE
 367       NP1=NP2+1
 368       IF(NP1.LE.NPTCL) GO TO 500
 369       RETURN
 370 C
 371 9999  CALL PRTEVT(0)
 372       WRITE(ITLIS,999) NPTCL
 373 999   FORMAT(//5X,'ERROR IN MBIAS...NPTCL >',I5)
 374       RETURN
 375       END