| 0795afa3 |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame