]>
Commit | Line | Data |
---|---|---|
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 |