103998: Fix in Hijing: to be committed to master, and then to Release
[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / hijing.F
CommitLineData
e74335a4 1* $Id$
2C Version 1.36
3C Nothing important has been changed here. A few 'garbage' has been
4C cleaned up here, like common block HIJJET3 for the sea quark strings
5C which were originally created to implement the DPM scheme which
6C later was abadoned in the final version. The lines which operate
7C on these data are also deleted in the program.
8C
9C
10C Version 1.35
11C There are some changes in the program: subroutine HARDJET is now
12C consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
13C for the triggered hard processes. Now that is done altogether
14C with other normal hard processes in modified JETINI. In the new
15C version one calls JETINI every time one calls HIJHRD. In the new
16C version the effect of the isospin of the nucleon on hard processes,
17C especially direct photons is correctly considered.
18C For A+A collisions, one has to initilize pythia
19C separately for each type of collisions, pp, pn,np and nn,
20C or hp and hn for hA collisions. In JETINI we use the following
21C catalogue for different types of collisions:
22C h+h: h+h (I_TYPE=1)
23C h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
24C A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
25C A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
26C*****************************************************************
27c
28C
29C Version 1.34
30C Last modification on January 5, 1998. Two misstakes are corrected in
31C function G. A Misstake in the subroutine Parton is also corrected.
32C (These are pointed out by Ysushi Nara).
33C
34C
35C Last modifcation on April 10, 1996. To conduct final
36C state radiation, PYTHIA reorganize the two scattered
37C partons and their final momenta will be a little
38C different. The summed total momenta of the partons
39C from the final state radiation are stored in HINT1(26-29)
40C and HINT1(36-39) which are little different from
41C HINT1(21-24) and HINT1(41-44).
42C
43C Version 1.33
44C
45C Last modfication on September 11, 1995. When HIJING and
46C PYTHIA are initialized, the shadowing is evaluated at
47C b=0 which is the maximum. This will cause overestimate
48C of shadowing for peripheral interactions. To correct this
49C problem, shadowing is set to zero when initializing. Then
50C use these maximum cross section without shadowing as a
51C normalization of the Monte Carlo. This however increase
52C the computing time. IHNT2(16) is used to indicate whether
53C the sturcture function is called for (IHNT2(16)=1) initialization
54C or for (IHNT2(16)=0)normal collisions simulation
55C
56C Last modification on Aagust 28, 1994. Two bugs associate
57C with the impact parameter dependence of the shadowing is
58C corrected.
59C
60C
61c Last modification on October 14, 1994. One bug is corrected
62c in the direct photon production option in subroutine
63C HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
64C Another bug associated with keeping the decay history
65C in the particle information is also corrected.(this problem
66C was reported by Matt Bloomer)
67C
68C
69C Last modification on July 15, 1994. The option to trig on
70C heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1)
71C is added. To do this, set IHPR2(3)=3. For inclusive production,
72C one should reset HIPR1(10)=0.0. One can also trig larger pt
73C QQbar production by giving HIPR1(10) a nonvanishing value.
74C The mass of the heavy quark in the calculation of the cross
75C section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
76C default is the charm mass D=1.5). We also include a separate
77C K-factor for heavy quark and direct photon production by
78C HIPR1(23)(D=2.0).
79C
80C Last modification on May 24, 1994. The option to
81C retain the information of all particles including those
82C who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is
83C added to contain the line number of the parent particle
84C of the current line which is produced via a decay.
85C KATT(I,4) is the status number of the particle: 11=particle
86C which has decayed; 1=finally produced particle.
87C
88C
89C Last modification on May 24, 1994( in HIJSFT when valence quark
90C is quenched, the following error is corrected. 1.2*IHNT2(1) -->
91C 1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
92C
93C
94C Last modification on March 16, 1994 (heavy flavor production
95C processes MSUB(81)=1 MSUB(82)=1 have been switched on,
96C charm production is the default, B-quark option is
97C IHPR2(18), when it is switched on, charm quark is
98C automatically off)
99C
100C
101C Last modification on March 23, 1994 (an error is corrected
102C in the impact parameter dependence of the jet cross section)
103C
104C Last modification Oct. 1993 to comply with non-vax
105C machines' compiler
106C
107C*********************************************
108C LAST MODIFICATION April 5, 1991
109CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B
110C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
111C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
112C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18),
113C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
114C
115C gluon jets can form a single string system
116C
117C initial state radiation is included
118C
119C all QCD subprocesses are included
120c
121c direct particles production is included(currently only direct
122C photon)
123c
124C Effect of high P_T trigger bias on multiple jets distribution
125c
126C******************************************************************
127C HIJING.10 *
128C Heavy Ion Jet INteraction Generator *
129C by *
130C X. N. Wang and M. Gyulassy *
131C Lawrence Berkeley Laboratory *
132C *
133C******************************************************************
134C
135C******************************************************************
136C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of *
137C proj, NFP(K,4) original ID of proj. NFP(K,5)=colli status(0=no,*
138C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
139C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it *
140C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
141C (0=has not been,1=is going to be, -1=has already been scattered *
142C NFP(k,11) total number of interactions this proj has suffered *
143C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant *
144C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and *
145C diquark,PP(K,10)=PT of the hard scattering between the valence *
146C quarks; PP(K,14,15)=the mass of quark,diquark. *
147C******************************************************************
148C
149C****************************************************************
150C
151C SUBROUTINE HIJING
152C
153C****************************************************************
154 SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
155 CHARACTER FRAME*8
156 DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
7dd84b41 157 & IPCOL(90000),ITCOL(90000),SCIP2(300,300)
bc676b8e 158#define BLANKET_SAVE
e74335a4 159#include "hiparnt.inc"
160C
161#include "hijcrdn.inc"
66a95620 162#include "hijglbr.inc"
e74335a4 163#include "himain1.inc"
164#include "himain2.inc"
165#include "histrng.inc"
166#include "hijjet1.inc"
167#include "hijjet2.inc"
168#include "hijjet4.inc"
169C
170#include "lujets_hijing.inc"
171#include "ludat1_hijing.inc"
172 SAVE
173
174 BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
175 BMIN=MIN(BMIN0,BMAX)
176 IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
177 BMIN=0.0
178 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
179 ENDIF
180C ********HIPR1(31) is in mb =0.1fm**2
181C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
182C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
183C
184 YP(1,1)=0.0
185 YP(2,1)=0.0
186 YP(3,1)=0.0
187 IF(IHNT2(1).LE.1) GO TO 14
188 DO 10 KP=1,IHNT2(1)
1895 R=HIRND(1)
cb3642a4 190c
191 if(IHNT2(1).EQ.2) then
6082051f 192 rnd1=max(RLU_HIJING(NSEED),1.0e-20)
193 rnd2=max(RLU_HIJING(NSEED),1.0e-20)
194 rnd3=max(RLU_HIJING(NSEED),1.0e-20)
cb3642a4 195 R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
196 & +4.38*0.85*log(rnd3)/(4.38+0.85))
197 endif
198c
e74335a4 199 X=RLU_HIJING(0)
200 CX=2.0*X-1.0
201 SX=SQRT(1.0-CX*CX)
202C ********choose theta from uniform cos(theta) distr
203 PHI=RLU_HIJING(0)*2.0*HIPR1(40)
204C ********choose phi form uniform phi distr 0 to 2*pi
205 YP(1,KP)=R*SX*COS(PHI)
206 YP(2,KP)=R*SX*SIN(PHI)
207 YP(3,KP)=R*CX
208 IF(HIPR1(29).EQ.0.0) GO TO 10
209 DO 8 KP2=1,KP-1
210 DNBP1=(YP(1,KP)-YP(1,KP2))**2
211 DNBP2=(YP(2,KP)-YP(2,KP2))**2
212 DNBP3=(YP(3,KP)-YP(3,KP2))**2
213 DNBP=DNBP1+DNBP2+DNBP3
214 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
215C ********two neighbors cannot be closer than
216C HIPR1(29)
2178 CONTINUE
21810 CONTINUE
cb3642a4 219c*******************************
220 if(IHNT2(1).EQ.2) then
221 YP(1,2)=-YP(1,1)
222 YP(2,2)=-YP(2,1)
223 YP(3,2)=-YP(3,1)
224 endif
225c********************************
e74335a4 226 DO 12 I=1,IHNT2(1)-1
227 DO 12 J=I+1,IHNT2(1)
228 IF(YP(3,I).GT.YP(3,J)) GO TO 12
229 Y1=YP(1,I)
230 Y2=YP(2,I)
231 Y3=YP(3,I)
232 YP(1,I)=YP(1,J)
233 YP(2,I)=YP(2,J)
234 YP(3,I)=YP(3,J)
235 YP(1,J)=Y1
236 YP(2,J)=Y2
237 YP(3,J)=Y3
23812 CONTINUE
239C
240C******************************
24114 YT(1,1)=0.0
242 YT(2,1)=0.0
243 YT(3,1)=0.0
244 IF(IHNT2(3).LE.1) GO TO 24
245 DO 20 KT=1,IHNT2(3)
24615 R=HIRND(2)
cb3642a4 247c
248 if(IHNT2(3).EQ.2) then
6082051f 249 rnd1=max(RLU_HIJING(NSEED),1.0e-20)
250 rnd2=max(RLU_HIJING(NSEED),1.0e-20)
251 rnd3=max(RLU_HIJING(NSEED),1.0e-20)
cb3642a4 252 R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
253 & +4.38*0.85*log(rnd3)/(4.38+0.85))
254 endif
255c
e74335a4 256 X=RLU_HIJING(0)
257 CX=2.0*X-1.0
258 SX=SQRT(1.0-CX*CX)
259C ********choose theta from uniform cos(theta) distr
260 PHI=RLU_HIJING(0)*2.0*HIPR1(40)
261C ********chose phi form uniform phi distr 0 to 2*pi
262 YT(1,KT)=R*SX*COS(PHI)
263 YT(2,KT)=R*SX*SIN(PHI)
264 YT(3,KT)=R*CX
265 IF(HIPR1(29).EQ.0.0) GO TO 20
266 DO 18 KT2=1,KT-1
267 DNBT1=(YT(1,KT)-YT(1,KT2))**2
268 DNBT2=(YT(2,KT)-YT(2,KT2))**2
269 DNBT3=(YT(3,KT)-YT(3,KT2))**2
270 DNBT=DNBT1+DNBT2+DNBT3
271 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
272C ********two neighbors cannot be closer than
273C HIPR1(29)
27418 CONTINUE
27520 CONTINUE
cb3642a4 276c**********************************
277 if(IHNT2(3).EQ.2) then
278 YT(1,2)=-YT(1,1)
279 YT(2,2)=-YT(2,1)
280 YT(3,2)=-YT(3,1)
281 endif
282c*********************************
e74335a4 283 DO 22 I=1,IHNT2(3)-1
284 DO 22 J=I+1,IHNT2(3)
285 IF(YT(3,I).LT.YT(3,J)) GO TO 22
286 Y1=YT(1,I)
287 Y2=YT(2,I)
288 Y3=YT(3,I)
289 YT(1,I)=YT(1,J)
290 YT(2,I)=YT(2,J)
291 YT(3,I)=YT(3,J)
292 YT(1,J)=Y1
293 YT(2,J)=Y2
294 YT(3,J)=Y3
29522 CONTINUE
296
297C********************
29824 MISS=-1
31960fe6
AM
299C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
300C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
301C BY THE ANGLE PHI FOR EACH COLLISION.******************
302C
303 BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
304 PHI=2.0*HIPR1(40)*RLU_HIJING(0)
e74335a4 305
30650 MISS=MISS+1
307 IF(MISS.GT.50) THEN
308 WRITE(6,*) 'infinite loop happened in HIJING'
309 STOP
310 ENDIF
311
312 NATT=0
313 JATT=0
314 EATT=0.0
315 CALL HIJINI
316 NLOP=0
317C ********Initialize for a new event
31860 NT=0
319 NP=0
320 N0=0
321 N01=0
322 N10=0
323 N11=0
8e529af3 324 NELT=0
325 NINT=0
326 NELP=0
327 NINP=0
e74335a4 328 NSG=0
329 NCOLT=0
65a17253 330 NPSPECP=0
331 NNSPECP=0
332 NPSPECT=0
333 NNSPECT=0
e74335a4 334
e74335a4 335 BBX=BB*COS(PHI)
336 BBY=BB*SIN(PHI)
337 HINT1(19)=BB
338 HINT1(20)=PHI
339C
340 DO 70 JP=1,IHNT2(1)
341 DO 70 JT=1,IHNT2(3)
342 SCIP(JP,JT)=-1.0
343 B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
344 R2=B2*HIPR1(40)/HIPR1(31)/0.1
7dd84b41 345 SCIP2(JP,JT)=R2
e74335a4 346C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
347 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
348 & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
349 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
350 & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
351 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
352 & *SQRT(1.0-RRB1)
353 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
354 & *SQRT(1.0-RRB2)
355 HINT1(18)=HINT1(14)-APHX1*HINT1(15)
356 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
357 IF(IHPR2(14).EQ.0.OR.
358 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
359 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
360 RANTOT=RLU_HIJING(0)
361 IF(RANTOT.GT.GS) GO TO 70
362 GO TO 65
363 ENDIF
364 GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
365 & /HIPR1(31)/2.0*ROMG(0.0)))
366 R2=R2/GSTOT_0
367 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
368 GSTOT=2.0*(1.0-SQRT(1.0-GS))
369 RANTOT=RLU_HIJING(0)*GSTOT_0
370 IF(RANTOT.GT.GSTOT) GO TO 70
371 IF(RANTOT.GT.GS) THEN
372 CALL HIJCSC(JP,JT)
373 GO TO 70
374C ********perform elastic collisions
375 ENDIF
376 65 SCIP(JP,JT)=R2
377 RNIP(JP,JT)=RANTOT
378 SJIP(JP,JT)=HINT1(18)
379 NCOLT=NCOLT+1
7dd84b41 380 if (R2.GT.2.D0) THEN
381 write (8,*) R2
382 ENDIF
e74335a4 383 IPCOL(NCOLT)=JP
384 ITCOL(NCOLT)=JT
38570 CONTINUE
db0ed1b7 386
db0ed1b7 387c *** cl glauber ***
db0ed1b7 388 npart=0
389 xmeana=0D0
390 ymeana=0D0
391 xmeanb=0D0
392 ymeanb=0D0
393 xmeanp=0D0
394 ymeanp=0D0
395 xm2=0D0
396 ym2=0D0
397 xym=0D0
8b7ea311 398c
fea3949c 399 IF(NCOLT>0) THEN
400 DO 1110 JP=1,IHNT2(1)
8b7ea311
AM
401 xmeana=xmeana+YP(1,JP)+BBX
402 ymeana=ymeana+YP(2,JP)+BBY
fea3949c 403 DO 1120 JT=1,IHNT2(3)
7dd84b41 404 IF(SCIP2(JP,JT).LT.2.0D0) THEN
fea3949c 405 npart=npart+1
8b7ea311
AM
406 xmeanp=xmeanp+YP(1,JP)+BBX
407 ymeanp=ymeanp+YP(2,JP)+BBY
408 xm2=xm2+(YP(1,JP)+BBX)*(YP(1,JP)+BBX)
409 ym2=ym2+(YP(2,JP)+BBY)*(YP(2,JP)+BBY)
410 xym=xym+(YP(1,JP)+BBX)*(YP(2,JP)+BBY)
fea3949c 411 goto 1110
412 end if
413 1120 continue
414 1110 continue
8b7ea311 415c
fea3949c 416 DO 1130 JT=1,IHNT2(3)
417 xmeanb=xmeanb+YT(1,JT)
418 ymeanb=ymeanb+YT(2,JT)
419 DO 1140 JP=1,IHNT2(1)
7dd84b41 420 IF(SCIP2(JP,JT).LT.2.0D0) THEN
fea3949c 421 npart=npart+1
422 xmeanp=xmeanp+YT(1,JT)
423 ymeanp=ymeanp+YT(2,JT)
424 xm2=xm2+YT(1,JT)*YT(1,JT)
425 ym2=ym2+YT(2,JT)*YT(2,JT)
426 xym=xym+YT(1,JT)*YT(2,JT)
427 goto 1130
428 end if
429 1140 continue
430 1130 continue
db0ed1b7 431
8832289a 432 IF (npart.GT.0) THEN
433 xmeana=xmeana/IHNT2(1)
434 ymeana=ymeana/IHNT2(1)
435 xmeanb=xmeanb/IHNT2(3)
436 ymeanb=ymeanb/IHNT2(3)
437 xmeanp=xmeanp/npart
438 ymeanp=ymeanp/npart
439 xm2=xm2/npart
440 ym2=ym2/npart
441 xym=xym/npart
8b7ea311 442c
8832289a 443 sx2=xm2-xmeanp*xmeanp
444 sy2=ym2-ymeanp*ymeanp
445 sxy=xym-xmeanp*ymeanp
8b7ea311 446c
8832289a 447 delx=xmeanb-xmeana
448 dely=ymeanb-ymeana
449 dtmp=delx**2+dely**2
450 bbtrue=sqrt(dtmp)
451 dnumt=(sy2-sx2)*(delx**2-dely**2)-4D0*sxy*delx*dely
452 ddent=(sy2+sx2)*bbtrue**2
453 eccrp=dnumt/ddent
454 dtmp=(sy2-sx2)*(sy2-sx2)+4D0*sxy*sxy
455 eccpart=sqrt(dtmp)/(sx2+sy2)
456 eccmc=(sy2-sx2)/(sy2+sx2)
457 write(*,*),'HOUT: ',bb,' ',bbtrue,' ',ncolt,' ',npart,
8b7ea311 458 1 ' ',eccrp,' ',eccpart, bbx, bby
8832289a 459 end if
fea3949c 460 end if
8b7ea311 461c
e74335a4 462C ********total number interactions proj and targ has
463C suffered
7dd84b41 464
e74335a4 465 IF(NCOLT.EQ.0) THEN
466 NLOP=NLOP+1
467 IF(NLOP.LE.20.OR.
31960fe6
AM
468 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
469 BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
470 GO TO 60
471 ENDIF
46f33c91 472
e74335a4 473 RETURN
474 ENDIF
475C ********At large impact parameter, there maybe no
476C interaction at all. For NN collision
477C repeat the event until interaction happens
478C
479 IF(IHPR2(3).NE.0) THEN
480 NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
481 NHARD=MIN(NHARD,NCOLT)
482 JPHARD=IPCOL(NHARD)
483 JTHARD=ITCOL(NHARD)
484 ENDIF
485C
486 IF(IHPR2(9).EQ.1) THEN
487 NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
488 NMINI=MIN(NMINI,NCOLT)
489 JPMINI=IPCOL(NMINI)
490 JTMINI=ITCOL(NMINI)
491 ENDIF
492C ********Specifying the location of the hard and
493C minijet if they are enforced by user
494C
495 DO 200 JP=1,IHNT2(1)
496 DO 200 JT=1,IHNT2(3)
497 IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
498 NFP(JP,11)=NFP(JP,11)+1
499 NFT(JT,11)=NFT(JT,11)+1
500 IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
501 NP=NP+1
502 N01=N01+1
503 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
504 NT=NT+1
505 N10=N10+1
506 ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
507 NP=NP+1
508 NT=NT+1
509 N0=N0+1
510 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
511 N11=N11+1
512 ENDIF
513 JOUT=0
514 NFP(JP,10)=0
515 NFT(JT,10)=0
516C*****************************************************************
517 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
518C ********When IHPR2(8)=0 no jets are produced
519 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
520C ********jets can not be produced for (JP,JT)
521C because not enough energy avaible for
522C JP or JT
523 R2=SCIP(JP,JT)
524 HINT1(18)=SJIP(JP,JT)
525 TT=ROMG(R2)*HINT1(18)/HIPR1(31)
526 TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
527 NJET=0
528 IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
529 CALL JETINI(JP,JT,1)
530 CALL HIJHRD(JP,JT,0,JFLG,0)
531 HINT1(26)=HINT1(47)
532 HINT1(27)=HINT1(48)
533 HINT1(28)=HINT1(49)
534 HINT1(29)=HINT1(50)
535 HINT1(36)=HINT1(67)
536 HINT1(37)=HINT1(68)
537 HINT1(38)=HINT1(69)
538 HINT1(39)=HINT1(70)
539C
540 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
541 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
542 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
543 & JFLG.GE.3) IASG(NSG,3)=1
544 IHNT2(9)=IHNT2(14)
545 IHNT2(10)=IHNT2(15)
546 DO 105 I05=1,5
547 HINT1(20+I05)=HINT1(40+I05)
548 HINT1(30+I05)=HINT1(50+I05)
549 105 CONTINUE
550 JOUT=1
551 IF(IHPR2(8).EQ.0) GO TO 160
552 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
553 & /REAL(IHNT2(1))**0.6666667,1.0)
554 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
555 & /REAL(IHNT2(3))**0.6666667,1.0)
556 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
557 & *SQRT(1.0-RRB1)
558 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
559 & *SQRT(1.0-RRB2)
560 HINT1(65)=HINT1(61)-APHX1*HINT1(62)
561 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
562 TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
563 NJET=-1
564C ********subtract the trigger jet from total number
565C of jet production to be done since it has
566C already been produced here
567 XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
568 106 NJET=NJET+1
569 XR1=XR1-ALOG(RLU_HIJING(0))
570 IF(XR1.LT.TTRIG) GO TO 106
571 XR=0.0
572 107 NJET=NJET+1
573 XR=XR-ALOG(RLU_HIJING(0))
574 IF(XR.LT.TT-TTRIG) GO TO 107
575 NJET=NJET-1
576 GO TO 112
577 ENDIF
578C ********create a hard interaction with specified P_T
579c when IHPR2(3)>0
580 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
581C ********create at least one pair of mini jets
582C when IHPR2(9)=1
583C
584 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
585 & (1.0-EXP(-TTS))) GO TO 160
586C ********this is the probability for no jet production
587110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
588111 NJET=NJET+1
589 XR=XR-ALOG(RLU_HIJING(0))
590 IF(XR.LT.TT) GO TO 111
591112 NJET=MIN(NJET,IHPR2(8))
592 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
593C ******** Determine number of mini jet production
594C
595 DO 150 I_JET=1,NJET
596 CALL JETINI(JP,JT,0)
597 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
598C ********JFLG=1 jets valence quarks, JFLG=2 with
599C gluon jet, JFLG=3 with q-qbar prod for
600C (JP,JT). If JFLG=0 jets can not be produced
601C this time. If JFLG=-1, error occured abandon
602C this event. JOUT is the total hard scat for
603C (JP,JT) up to now.
604 IF(JFLG.EQ.0) GO TO 160
605 IF(JFLG.LT.0) THEN
606 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
607 GO TO 50
608 ENDIF
609 JOUT=JOUT+1
610 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
611 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
612 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
613 & JFLG.GE.3) IASG(NSG,3)=1
614C ******** jet with PT>HIPR1(11) will be quenched
615 150 CONTINUE
616 160 CONTINUE
617 CALL HIJSFT(JP,JT,JOUT,IERROR)
618 IF(IERROR.NE.0) THEN
619 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
620 GO TO 50
621 ENDIF
622C
623C ********conduct soft scattering between JP and JT
624 JATT=JATT+JOUT
625
626200 CONTINUE
627
8e529af3 628c
629c**************************
630c
631 DO 201 JP=1,IHNT2(1)
65a17253 632c write(6,*) JP, NFP(JP,3), NFP(JP,4), NFP(JP,5)
8e529af3 633 IF(NFP(JP,5).GT.2) THEN
634 NINP=NINP+1
635 ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
636 NELP=NELP+1
637 ENDIF
65a17253 638
639 IF(NFP(JP,5).LE.2) THEN
640 IF (NFP(JP,3) .EQ. 2212) THEN
641 NPSPECP = NPSPECP + 1
642 ELSE IF (NFP(JP,3) .EQ. 2112) THEN
643 NNSPECP = NNSPECP + 1
644 ENDIF
645 ENDIF
8e529af3 646 201 continue
647 DO 202 JT=1,IHNT2(3)
648 IF(NFT(JT,5).GT.2) THEN
649 NINT=NINT+1
650 ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
651 NELT=NELT+1
652 ENDIF
65a17253 653
654 IF(NFT(JT,5).LE.2) THEN
655 IF (NFT(JT,3) .EQ. 2212) THEN
656 NPSPECT = NPSPECT + 1
657 ELSE IF (NFT(JT,3) .EQ. 2112) THEN
658 NNSPECT = NNSPECT + 1
659 ENDIF
660 ENDIF
8e529af3 661 202 continue
662c
663c*******************************
664
e74335a4 665C********perform jet quenching for jets with PT>HIPR1(11)**********
666
667 IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
668 & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
669 DO 271 I=1,IHNT2(1)
670 IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
671271 CONTINUE
672 DO 272 I=1,IHNT2(3)
673 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
674272 CONTINUE
675 DO 273 ISG=1,NSG
676 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
677273 CONTINUE
678 ENDIF
679C
680C**************fragment all the string systems in the following*****
681C
682C********N_ST is where particle information starts
683C********N_STR+1 is the number of strings in fragmentation
684C********the number of strings before a line is stored in K(I,4)
685C********IDSTR is id number of the string system (91,92 or 93)
686C
cb3642a4 687c
e74335a4 688 IF(IHPR2(20).NE.0) THEN
689 DO 360 ISG=1,NSG
690 CALL HIJFRG(ISG,3,IERROR)
691 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
692 MSTU(24)=0
693 MSTU(28)=0
694 IF(IHPR2(10).NE.0) THEN
695 call LULIST_HIJING(1)
696 WRITE(6,*) 'error occured, repeat the event'
697 ENDIF
698 GO TO 50
699 ENDIF
700C ********Check errors
701C
702 N_ST=1
703 IDSTR=92
704 IF(IHPR2(21).EQ.0) THEN
705 CALL LUEDIT_HIJING(2)
dca9e44f 706 ELSE
db0ed1b7 707351 N_ST=N_ST+1
96200d6d 708 IF((N_ST .LE. 9000)
709 & .AND. (K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93)) GO TO 351
db0ed1b7 710 IDSTR=K(N_ST,2)
711 N_ST=N_ST+1
e74335a4 712 ENDIF
713C
714 IF(FRAME.EQ.'LAB') THEN
715 CALL HIBOOST
716 ENDIF
717C ******** boost back to lab frame(if it was in)
718C
719 N_STR=0
720 DO 360 I=N_ST,N
721 IF(K(I,2).EQ.IDSTR) THEN
722 N_STR=N_STR+1
723 GO TO 360
724 ENDIF
725 K(I,4)=N_STR
726 NATT=NATT+1
727 KATT(NATT,1)=K(I,2)
728 KATT(NATT,2)=20
729 KATT(NATT,4)=K(I,1)
65a17253 730 IF(K(I,3).EQ.0 ) THEN
731 KATT(NATT,3)=0
732 ELSE
733 IF(K(K(I,3),2).EQ.IDSTR) THEN
734 KATT(NATT,3)=0
735 ELSE
736 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
737 ENDIF
738 ENDIF
e74335a4 739C ****** identify the mother particle
740 PATT(NATT,1)=P(I,1)
741 PATT(NATT,2)=P(I,2)
742 PATT(NATT,3)=P(I,3)
743 PATT(NATT,4)=P(I,4)
1a0dd691 744 VATT(NATT,1)=V(I,1)
745 VATT(NATT,2)=V(I,2)
746 VATT(NATT,3)=V(I,3)
747 VATT(NATT,4)=V(I,4)
748
e74335a4 749 EATT=EATT+P(I,4)
750360 CONTINUE
751C ********Fragment the q-qbar jets systems *****
752C
753 JTP(1)=IHNT2(1)
754 JTP(2)=IHNT2(3)
755 DO 400 NTP=1,2
756 DO 400 J_JTP=1,JTP(NTP)
757 CALL HIJFRG(J_JTP,NTP,IERROR)
758 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
759 MSTU(24)=0
760 MSTU(28)=0
761 IF(IHPR2(10).NE.0) THEN
762 call LULIST_HIJING(1)
763 WRITE(6,*) 'error occured, repeat the event'
764 ENDIF
765 GO TO 50
766 ENDIF
767C ********check errors
768C
769 N_ST=1
770 IDSTR=92
46f33c91 771
65a17253 772 NFTP=NFP(J_JTP,5)
773 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
774
e74335a4 775 IF(IHPR2(21).EQ.0) THEN
776 CALL LUEDIT_HIJING(2)
65a17253 777 ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN
e74335a4 778381 N_ST=N_ST+1
779 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
780 IDSTR=K(N_ST,2)
781 N_ST=N_ST+1
46f33c91 782 ENDIF
e74335a4 783 IF(FRAME.EQ.'LAB') THEN
784 CALL HIBOOST
785 ENDIF
786C ******** boost back to lab frame(if it was in)
787C
65a17253 788
e74335a4 789 N_STR=0
790 DO 390 I=N_ST,N
791 IF(K(I,2).EQ.IDSTR) THEN
792 N_STR=N_STR+1
793 GO TO 390
794 ENDIF
795 K(I,4)=N_STR
796 NATT=NATT+1
797 KATT(NATT,1)=K(I,2)
798 KATT(NATT,2)=NFTP
799 KATT(NATT,4)=K(I,1)
cb3642a4 800 IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
e74335a4 801 KATT(NATT,3)=0
cb3642a4 802 ELSE
803 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
e74335a4 804 ENDIF
805C ****** identify the mother particle
806 PATT(NATT,1)=P(I,1)
807 PATT(NATT,2)=P(I,2)
808 PATT(NATT,3)=P(I,3)
809 PATT(NATT,4)=P(I,4)
810 EATT=EATT+P(I,4)
4181e24f 811 VATT(NATT,1)=V(I,1)
812 VATT(NATT,2)=V(I,2)
813 VATT(NATT,3)=V(I,3)
814 VATT(NATT,4)=V(I,4)
e74335a4 815390 CONTINUE
816400 CONTINUE
817C ********Fragment the q-qq related string systems
818 ENDIF
819
820 DO 450 I=1,NDR
821 NATT=NATT+1
822 KATT(NATT,1)=KFDR(I)
823 KATT(NATT,2)=40
824 KATT(NATT,3)=0
825 PATT(NATT,1)=PDR(I,1)
826 PATT(NATT,2)=PDR(I,2)
827 PATT(NATT,3)=PDR(I,3)
828 PATT(NATT,4)=PDR(I,4)
4181e24f 829 VATT(NATT,1)=V(I,1)
830 VATT(NATT,2)=V(I,2)
831 VATT(NATT,3)=V(I,3)
832 VATT(NATT,4)=V(I,4)
e74335a4 833 EATT=EATT+PDR(I,4)
834450 CONTINUE
835C ********store the direct-produced particles
836C
837 DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
838 IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
839 & .AND.IHPR2(21).EQ.0) THEN
840 IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
841 & event'
842C call LULIST_HIJING(1)
843 GO TO 50
844 ENDIF
46f33c91 845
e74335a4 846 RETURN
847 END