Correct ran function.
[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),
157 & IPCOL(90000),ITCOL(90000)
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
299
30050 MISS=MISS+1
301 IF(MISS.GT.50) THEN
302 WRITE(6,*) 'infinite loop happened in HIJING'
303 STOP
304 ENDIF
305
306 NATT=0
307 JATT=0
308 EATT=0.0
309 CALL HIJINI
310 NLOP=0
311C ********Initialize for a new event
31260 NT=0
313 NP=0
314 N0=0
315 N01=0
316 N10=0
317 N11=0
8e529af3 318 NELT=0
319 NINT=0
320 NELP=0
321 NINP=0
e74335a4 322 NSG=0
323 NCOLT=0
324
325C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
326C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
327C BY THE ANGLE PHI FOR EACH COLLISION.******************
328C
329 BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
330 PHI=2.0*HIPR1(40)*RLU_HIJING(0)
331 BBX=BB*COS(PHI)
332 BBY=BB*SIN(PHI)
333 HINT1(19)=BB
334 HINT1(20)=PHI
335C
336 DO 70 JP=1,IHNT2(1)
337 DO 70 JT=1,IHNT2(3)
338 SCIP(JP,JT)=-1.0
339 B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
340 R2=B2*HIPR1(40)/HIPR1(31)/0.1
341C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
342 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
343 & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
344 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
345 & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
346 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
347 & *SQRT(1.0-RRB1)
348 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
349 & *SQRT(1.0-RRB2)
350 HINT1(18)=HINT1(14)-APHX1*HINT1(15)
351 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
352 IF(IHPR2(14).EQ.0.OR.
353 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
354 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
355 RANTOT=RLU_HIJING(0)
356 IF(RANTOT.GT.GS) GO TO 70
357 GO TO 65
358 ENDIF
359 GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
360 & /HIPR1(31)/2.0*ROMG(0.0)))
361 R2=R2/GSTOT_0
362 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
363 GSTOT=2.0*(1.0-SQRT(1.0-GS))
364 RANTOT=RLU_HIJING(0)*GSTOT_0
365 IF(RANTOT.GT.GSTOT) GO TO 70
366 IF(RANTOT.GT.GS) THEN
367 CALL HIJCSC(JP,JT)
368 GO TO 70
369C ********perform elastic collisions
370 ENDIF
371 65 SCIP(JP,JT)=R2
372 RNIP(JP,JT)=RANTOT
373 SJIP(JP,JT)=HINT1(18)
374 NCOLT=NCOLT+1
375 IPCOL(NCOLT)=JP
376 ITCOL(NCOLT)=JT
37770 CONTINUE
378C ********total number interactions proj and targ has
379C suffered
380 IF(NCOLT.EQ.0) THEN
381 NLOP=NLOP+1
382 IF(NLOP.LE.20.OR.
383 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
46f33c91 384
385
e74335a4 386 RETURN
387 ENDIF
388C ********At large impact parameter, there maybe no
389C interaction at all. For NN collision
390C repeat the event until interaction happens
391C
392 IF(IHPR2(3).NE.0) THEN
393 NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
394 NHARD=MIN(NHARD,NCOLT)
395 JPHARD=IPCOL(NHARD)
396 JTHARD=ITCOL(NHARD)
397 ENDIF
398C
399 IF(IHPR2(9).EQ.1) THEN
400 NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
401 NMINI=MIN(NMINI,NCOLT)
402 JPMINI=IPCOL(NMINI)
403 JTMINI=ITCOL(NMINI)
404 ENDIF
405C ********Specifying the location of the hard and
406C minijet if they are enforced by user
407C
408 DO 200 JP=1,IHNT2(1)
409 DO 200 JT=1,IHNT2(3)
410 IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
411 NFP(JP,11)=NFP(JP,11)+1
412 NFT(JT,11)=NFT(JT,11)+1
413 IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
414 NP=NP+1
415 N01=N01+1
416 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
417 NT=NT+1
418 N10=N10+1
419 ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
420 NP=NP+1
421 NT=NT+1
422 N0=N0+1
423 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
424 N11=N11+1
425 ENDIF
426 JOUT=0
427 NFP(JP,10)=0
428 NFT(JT,10)=0
429C*****************************************************************
430 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
431C ********When IHPR2(8)=0 no jets are produced
432 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
433C ********jets can not be produced for (JP,JT)
434C because not enough energy avaible for
435C JP or JT
436 R2=SCIP(JP,JT)
437 HINT1(18)=SJIP(JP,JT)
438 TT=ROMG(R2)*HINT1(18)/HIPR1(31)
439 TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
440 NJET=0
441 IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
442 CALL JETINI(JP,JT,1)
443 CALL HIJHRD(JP,JT,0,JFLG,0)
444 HINT1(26)=HINT1(47)
445 HINT1(27)=HINT1(48)
446 HINT1(28)=HINT1(49)
447 HINT1(29)=HINT1(50)
448 HINT1(36)=HINT1(67)
449 HINT1(37)=HINT1(68)
450 HINT1(38)=HINT1(69)
451 HINT1(39)=HINT1(70)
452C
453 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
454 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
455 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
456 & JFLG.GE.3) IASG(NSG,3)=1
457 IHNT2(9)=IHNT2(14)
458 IHNT2(10)=IHNT2(15)
459 DO 105 I05=1,5
460 HINT1(20+I05)=HINT1(40+I05)
461 HINT1(30+I05)=HINT1(50+I05)
462 105 CONTINUE
463 JOUT=1
464 IF(IHPR2(8).EQ.0) GO TO 160
465 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
466 & /REAL(IHNT2(1))**0.6666667,1.0)
467 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
468 & /REAL(IHNT2(3))**0.6666667,1.0)
469 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
470 & *SQRT(1.0-RRB1)
471 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
472 & *SQRT(1.0-RRB2)
473 HINT1(65)=HINT1(61)-APHX1*HINT1(62)
474 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
475 TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
476 NJET=-1
477C ********subtract the trigger jet from total number
478C of jet production to be done since it has
479C already been produced here
480 XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
481 106 NJET=NJET+1
482 XR1=XR1-ALOG(RLU_HIJING(0))
483 IF(XR1.LT.TTRIG) GO TO 106
484 XR=0.0
485 107 NJET=NJET+1
486 XR=XR-ALOG(RLU_HIJING(0))
487 IF(XR.LT.TT-TTRIG) GO TO 107
488 NJET=NJET-1
489 GO TO 112
490 ENDIF
491C ********create a hard interaction with specified P_T
492c when IHPR2(3)>0
493 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
494C ********create at least one pair of mini jets
495C when IHPR2(9)=1
496C
497 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
498 & (1.0-EXP(-TTS))) GO TO 160
499C ********this is the probability for no jet production
500110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
501111 NJET=NJET+1
502 XR=XR-ALOG(RLU_HIJING(0))
503 IF(XR.LT.TT) GO TO 111
504112 NJET=MIN(NJET,IHPR2(8))
505 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
506C ******** Determine number of mini jet production
507C
508 DO 150 I_JET=1,NJET
509 CALL JETINI(JP,JT,0)
510 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
511C ********JFLG=1 jets valence quarks, JFLG=2 with
512C gluon jet, JFLG=3 with q-qbar prod for
513C (JP,JT). If JFLG=0 jets can not be produced
514C this time. If JFLG=-1, error occured abandon
515C this event. JOUT is the total hard scat for
516C (JP,JT) up to now.
517 IF(JFLG.EQ.0) GO TO 160
518 IF(JFLG.LT.0) THEN
519 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
520 GO TO 50
521 ENDIF
522 JOUT=JOUT+1
523 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
524 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
525 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
526 & JFLG.GE.3) IASG(NSG,3)=1
527C ******** jet with PT>HIPR1(11) will be quenched
528 150 CONTINUE
529 160 CONTINUE
530 CALL HIJSFT(JP,JT,JOUT,IERROR)
531 IF(IERROR.NE.0) THEN
532 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
533 GO TO 50
534 ENDIF
535C
536C ********conduct soft scattering between JP and JT
537 JATT=JATT+JOUT
538
539200 CONTINUE
540
8e529af3 541c
542c**************************
543c
544 DO 201 JP=1,IHNT2(1)
545 IF(NFP(JP,5).GT.2) THEN
546 NINP=NINP+1
547 ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
548 NELP=NELP+1
549 ENDIF
550 201 continue
551 DO 202 JT=1,IHNT2(3)
552 IF(NFT(JT,5).GT.2) THEN
553 NINT=NINT+1
554 ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
555 NELT=NELT+1
556 ENDIF
557 202 continue
558c
559c*******************************
560
e74335a4 561C********perform jet quenching for jets with PT>HIPR1(11)**********
562
563 IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
564 & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
565 DO 271 I=1,IHNT2(1)
566 IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
567271 CONTINUE
568 DO 272 I=1,IHNT2(3)
569 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
570272 CONTINUE
571 DO 273 ISG=1,NSG
572 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
573273 CONTINUE
574 ENDIF
575C
576C**************fragment all the string systems in the following*****
577C
578C********N_ST is where particle information starts
579C********N_STR+1 is the number of strings in fragmentation
580C********the number of strings before a line is stored in K(I,4)
581C********IDSTR is id number of the string system (91,92 or 93)
582C
cb3642a4 583 GAMMA0=HINT1(1)/2.0/0.93827
584 BETA0=sqrt(GAMMA0**2-1.0)/GAMMA0
585c
e74335a4 586 IF(IHPR2(20).NE.0) THEN
587 DO 360 ISG=1,NSG
588 CALL HIJFRG(ISG,3,IERROR)
589 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
590 MSTU(24)=0
591 MSTU(28)=0
592 IF(IHPR2(10).NE.0) THEN
593 call LULIST_HIJING(1)
594 WRITE(6,*) 'error occured, repeat the event'
595 ENDIF
596 GO TO 50
597 ENDIF
598C ********Check errors
599C
600 N_ST=1
601 IDSTR=92
602 IF(IHPR2(21).EQ.0) THEN
603 CALL LUEDIT_HIJING(2)
604 ELSE
605351 N_ST=N_ST+1
606 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351
607 IDSTR=K(N_ST,2)
608 N_ST=N_ST+1
609 ENDIF
610C
611 IF(FRAME.EQ.'LAB') THEN
612 CALL HIBOOST
613 ENDIF
614C ******** boost back to lab frame(if it was in)
615C
616 N_STR=0
617 DO 360 I=N_ST,N
618 IF(K(I,2).EQ.IDSTR) THEN
619 N_STR=N_STR+1
620 GO TO 360
621 ENDIF
622 K(I,4)=N_STR
623 NATT=NATT+1
624 KATT(NATT,1)=K(I,2)
625 KATT(NATT,2)=20
626 KATT(NATT,4)=K(I,1)
627 IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
628 KATT(NATT,3)=0
629 ELSE
630 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
631 ENDIF
632C ****** identify the mother particle
633 PATT(NATT,1)=P(I,1)
634 PATT(NATT,2)=P(I,2)
635 PATT(NATT,3)=P(I,3)
636 PATT(NATT,4)=P(I,4)
637 EATT=EATT+P(I,4)
cb3642a4 638 VATT01=0.5*(YP(1,IASG(ISG,1))+YT(1,IASG(ISG,2)))
639 VATT02=0.5*(YP(2,IASG(ISG,1))+YT(2,IASG(ISG,2)))
640 VATT03=0.5*(YP(3,IASG(ISG,1))+YT(3,IASG(ISG,2)))
641 & /GAMMA0
642 VATT04=-0.5*(YP(3,IASG(ISG,1))
643 & -YT(3,IASG(ISG,2)))/GAMMA0/BETA0
644 RARB=1.12*(IHNT2(1)**0.33333+IHNT2(3)**0.33333)
645 V3MIN1=RARB/GAMMA0
646 V3MIN2=1.0/MAX(1.0,5.08*ABS(PATT(I,3)))
647 VATT_MIN=MAX(V3MIN1,V3MIN2)
6082051f 648 VATT03=VATT03+(0.5-RLU_HIJING(0))*VATT_MIN
cb3642a4 649 amt2=P(I,1)**2+P(I,2)**2+P(I,5)**2
650 IF(amt2.GT.0.0) THEN
651 tauf=0.2*2.0*P(I,3)/amt2
652 VATT(NATT,1)=VATT01
653 VATT(NATT,2)=VATT02
654 VATT(NATT,3)=VATT03+tauf
655 VATT(NATT,4)=abs(VATT(NATT,3))
656 ELSE
657 VATT(NATT,4)=abs(VATT03)
658 VATT(NATT,1)=VATT01
659 VATT(NATT,2)=VATT02
660 VATT(NATT,3)=VATT03
661 ENDIF
e74335a4 662360 CONTINUE
663C ********Fragment the q-qbar jets systems *****
664C
665 JTP(1)=IHNT2(1)
666 JTP(2)=IHNT2(3)
667 DO 400 NTP=1,2
668 DO 400 J_JTP=1,JTP(NTP)
669 CALL HIJFRG(J_JTP,NTP,IERROR)
670 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
671 MSTU(24)=0
672 MSTU(28)=0
673 IF(IHPR2(10).NE.0) THEN
674 call LULIST_HIJING(1)
675 WRITE(6,*) 'error occured, repeat the event'
676 ENDIF
677 GO TO 50
678 ENDIF
679C ********check errors
680C
681 N_ST=1
682 IDSTR=92
46f33c91 683
e74335a4 684 IF(IHPR2(21).EQ.0) THEN
685 CALL LUEDIT_HIJING(2)
cb3642a4 686 ELSE
e74335a4 687381 N_ST=N_ST+1
688 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
689 IDSTR=K(N_ST,2)
690 N_ST=N_ST+1
46f33c91 691 ENDIF
e74335a4 692 IF(FRAME.EQ.'LAB') THEN
693 CALL HIBOOST
694 ENDIF
695C ******** boost back to lab frame(if it was in)
696C
cb3642a4 697 NFTP=NFP(J_JTP,5)
698 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
e74335a4 699 N_STR=0
700 DO 390 I=N_ST,N
701 IF(K(I,2).EQ.IDSTR) THEN
702 N_STR=N_STR+1
703 GO TO 390
704 ENDIF
705 K(I,4)=N_STR
706 NATT=NATT+1
707 KATT(NATT,1)=K(I,2)
708 KATT(NATT,2)=NFTP
709 KATT(NATT,4)=K(I,1)
cb3642a4 710 IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
e74335a4 711 KATT(NATT,3)=0
cb3642a4 712 ELSE
713 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
e74335a4 714 ENDIF
715C ****** identify the mother particle
716 PATT(NATT,1)=P(I,1)
717 PATT(NATT,2)=P(I,2)
718 PATT(NATT,3)=P(I,3)
719 PATT(NATT,4)=P(I,4)
720 EATT=EATT+P(I,4)
cb3642a4 721C ******* the following calculate the formation time and position
722C of the produced hadrons
723 IF(NTP.EQ.1) THEN
724 VATT01=YP(1,J_JTP)
725 VATT02=YP(2,J_JTP)
726 VATT03=YP(3,J_JTP)/GAMMA0
727 ELSE IF(NTP.EQ.2) THEN
728 VATT01=YT(1,J_JTP)
729 VATT02=YT(2,J_JTP)
730 VATT03=YT(3,J_JTP)/GAMMA0
731 ENDIF
732 RARB=1.12*(IHNT2(1)**0.33333+IHNT2(3)**0.33333)
733 V3MIN1=RARB/GAMMA0
734 V3MIN2=1.0/MAX(1.0,5.08*ABS(PATT(I,3)))
735 VATT_MIN=MAX(V3MIN1,V3MIN2)
6082051f 736 VATT03=VATT03+(0.5-RLU_HIJING(0))*VATT_MIN
cb3642a4 737 amt2=P(I,1)**2+P(I,2)**2+P(I,5)**2
738 IF(amt2.GT.0.0) THEN
739 tauf=0.2*2.0*P(I,3)/amt2
740 VATT(NATT,1)=VATT01
741 VATT(NATT,2)=VATT02
742 VATT(NATT,3)=VATT03+tauf
743 VATT(NATT,4)=abs(VATT(NATT,3))
744 ELSE
745 VATT(NATT,4)=abs(VATT03)
746 VATT(NATT,1)=VATT01
747 VATT(NATT,2)=VATT02
748 VATT(NATT,3)=VATT03
749 ENDIF
e74335a4 750390 CONTINUE
751400 CONTINUE
752C ********Fragment the q-qq related string systems
753 ENDIF
754
755 DO 450 I=1,NDR
756 NATT=NATT+1
757 KATT(NATT,1)=KFDR(I)
758 KATT(NATT,2)=40
759 KATT(NATT,3)=0
760 PATT(NATT,1)=PDR(I,1)
761 PATT(NATT,2)=PDR(I,2)
762 PATT(NATT,3)=PDR(I,3)
763 PATT(NATT,4)=PDR(I,4)
764 EATT=EATT+PDR(I,4)
cb3642a4 765 amt2=P(I,1)**2+P(I,2)**2+P(I,5)**2
766 IF(amt2.GT.0.0) THEN
767 tauf=0.2*2.0*P(I,3)/amt2
768 VATT(NATT,1)=0.0
769 VATT(NATT,2)=0.0
770 VATT(NATT,3)=tauf
771 VATT(NATT,4)=abs(VATT(NATT,3))
772 ELSE
773 VATT(NATT,4)=0.0
774 VATT(NATT,1)=0.0
775 VATT(NATT,2)=0.0
776 VATT(NATT,3)=0.0
777 ENDIF
e74335a4 778450 CONTINUE
779C ********store the direct-produced particles
780C
781 DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
782 IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
783 & .AND.IHPR2(21).EQ.0) THEN
784 IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
785 & event'
786C call LULIST_HIJING(1)
787 GO TO 50
788 ENDIF
46f33c91 789
e74335a4 790 RETURN
791 END