1 #include "isajet/pilot.h"
4 C Generate minbias event or beam jets for high-pt event using
5 C parameters set in MBSET:
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.
17 C Note that multiple cut Pomerons give approximate KNO scaling.
18 C The only short-range correlations are from resonances.
20 C Ver. 7.09: Add traps on free loops and IMPLICIT NONE.
22 #if defined(CERNLIB_IMPNONE)
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"
34 DIMENSION IFL(3),IFLEXC(2),PXEXC(2),PYEXC(2),SIGN(2)
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/
54 C Select number of cut Pomerons.
60 IF(POMGEN(N).GT.TRY) GO TO 20
64 C Decide if diffractive event
65 IF(RANF().LT.PDIFFR) THEN
66 IDIFF=INT(1.99999*RANF())+1
68 LDIFFR(3-IDIFF)=.FALSE.
74 C Generate leading baryons.
79 C Special treatment for diffractive beam.
83 CALL FLAVOR(IDIN(IB),IFL(1),IFL(2),IFL(3),JSPIN,INDEX)
88 DXBARY(IB)=EPSDIF**RANF()
89 XBARY(IB)=1.-DXBARY(IB)
93 C If not diffractive, construct IDENT of leading baryon
94 CALL FLAVOR(IDIN(IB),IFL(1),IFL(2),IFL(3),JSPIN,INDEX)
96 IFLNEW=ISIGN(INT(RANF()/PUD0)+1,IDIN(IB))
100 IF(IABS(IFL(1)).GT.IABS(IFL(2))) THEN
105 IF(IABS(IFL(2)).GT.IABS(IFL(3))) THEN
110 IF(IABS(IFL(1)).GT.IABS(IFL(2))) THEN
116 IF(IFL(1).EQ.IFL(2).AND.IFL(2).EQ.IFL(3)) THEN
119 JSPIN=INT(RANF()+PJSPN)
121 IF(JSPIN.EQ.0.AND.IFL(1).NE.IFL(2).AND.IFL(2).NE.IFL(3)) THEN
122 IF(RANF().GT.PISPN) THEN
128 IDHAD=1000*IABS(IFL(1))+100*IABS(IFL(2))+10*IABS(IFL(3))+JSPIN
129 IDHAD=ISIGN(IDHAD,IDIN(IB))
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)
138 C Select transverse momentum of baryon
139 115 CALL GETPT(PT1,SIGQT0)
145 CALL GETPT(PT2,SIGQT0)
153 AMT2=PX**2+PY**2+AM**2
155 QPLUS=XBARY(IB)*PPLUS
156 QPLUS=AMAX1(QPLUS,1.E-6)
161 C Add baryon to /PARTCL/ if PZ>0.
162 IF(NPTCL.GE.MXPTCL) GO TO 9999
167 PPTCL(3,NPTCL)=PZ*SIGN(IB)
175 IF(IPASSB.LT.MXPASS) GO TO 110
176 C Just give up if it fails MXPASS times
178 998 FORMAT(//5X,'ERROR IN MBIAS ... COULD NOT MAKE BARYON')
183 C Having accepted baryon, set up XPOM array for cut Pomerons,
184 C rescaling to 1.-XBARY(IB).
193 XPOM(N,IB)=XSUM*XPOM(N,IB)*DXBARY(IB)
197 C Fragment each Pomeron into mesons and baryon pairs in the
198 C Pomeron-Pomeron center of mass.
202 PPOM=SQRT(PBEAM(1)*XPOM(IPOM,1)*PBEAM(2)*XPOM(IPOM,2))
208 ITWIST=INT(1.99999*RANF())+1
211 C Select new quark or diquark. Old diquark implies new quark.
212 C Old quark implies new diquark with probability PBARY0.
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)
219 ID1=INT(RANF()/PUD0)+1
220 ID2=INT(RANF()/PUD0)+1
221 IF(IABS(ID1).GT.IABS(ID2)) THEN
226 IFL2=ISIGN(1000*ID1+100*ID2,IFL1)
229 C Construct meson from quark+antiquark. Else, construct baryon
230 C IDENT from quark+diquark.
233 IF(IMOD1.NE.0.AND.IMOD2.NE.0) THEN
234 JSPIN=INT(RANF()+PJSPN)
237 IF(ID1+ID2.EQ.0) THEN
240 ID1=INT(PMIX01(ID1,JSPIN+1)+RND)
241 $ +INT(PMIX02(ID1,JSPIN+1)+RND)+1
243 ELSEIF(IABS(ID1).GT.IABS(ID2)) THEN
248 IDHAD=ISIGN(100*IABS(ID1)+10*IABS(ID2)+JSPIN,ID1)
259 IF(IABS(ID1).GT.IABS(ID2)) THEN
264 IF(IABS(ID2).GT.IABS(ID3)) THEN
269 IF(IABS(ID1).GT.IABS(ID2)) THEN
274 IF(ID1.EQ.ID2.AND.ID2.EQ.ID3) THEN
277 JSPIN=INT(RANF()+PJSPN)
279 IF(JSPIN.EQ.0.AND.ID1.NE.ID2.AND.ID2.NE.ID3) THEN
280 IF(RANF().LT.PISPN) THEN
286 IDHAD=1000*IABS(ID1)+100*IABS(ID2)+10*IABS(ID3)+JSPIN
287 IDHAD=ISIGN(IDHAD,IFL1)
293 CALL GETPT(PT2,SIGQT0)
301 AMT2=PX**2+PY**2+AM**2
303 C Select x -- same distribution for all particles.
305 IF(RANF().LT.XGEN0(1)) X=1.-X**(1./(XGEN0(2)+1.))
307 QPLUS=AMAX1(QPLUS,1.E-6)
312 C Add particle to /PARTCL/ if PZ>0.
313 IF(NPTCL.GE.MXPTCL) GO TO 9999
318 PPTCL(3,NPTCL)=PZ*SIGN(IB)
326 C Continue if sufficient pplus
328 IF(PPLUS.GT.PEND0.AND.LOOP.LT.MXPTCL) GO TO 200
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)
344 C Rescale hadron momenta for correct four-momentum.
350 CALL RESCAL(NBEGIN,NPTLV1,PSUM,IFAIL)
352 CALL RESCAL(NBEGIN,NPTLV1,PBEAMS,IFAIL)
354 IF(IFAIL.NE.0.AND.IPASS.LT.MXPASS) THEN
368 IF(NP1.LE.NPTCL) GO TO 500
372 WRITE(ITLIS,999) NPTCL
373 999 FORMAT(//5X,'ERROR IN MBIAS...NPTCL >',I5)