updates in MeanTh scan analysis and new method added to display hitmaps (A. Mastroserio)
[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)
190 X=RLU_HIJING(0)
191 CX=2.0*X-1.0
192 SX=SQRT(1.0-CX*CX)
193C ********choose theta from uniform cos(theta) distr
194 PHI=RLU_HIJING(0)*2.0*HIPR1(40)
195C ********choose phi form uniform phi distr 0 to 2*pi
196 YP(1,KP)=R*SX*COS(PHI)
197 YP(2,KP)=R*SX*SIN(PHI)
198 YP(3,KP)=R*CX
199 IF(HIPR1(29).EQ.0.0) GO TO 10
200 DO 8 KP2=1,KP-1
201 DNBP1=(YP(1,KP)-YP(1,KP2))**2
202 DNBP2=(YP(2,KP)-YP(2,KP2))**2
203 DNBP3=(YP(3,KP)-YP(3,KP2))**2
204 DNBP=DNBP1+DNBP2+DNBP3
205 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
206C ********two neighbors cannot be closer than
207C HIPR1(29)
2088 CONTINUE
20910 CONTINUE
210 DO 12 I=1,IHNT2(1)-1
211 DO 12 J=I+1,IHNT2(1)
212 IF(YP(3,I).GT.YP(3,J)) GO TO 12
213 Y1=YP(1,I)
214 Y2=YP(2,I)
215 Y3=YP(3,I)
216 YP(1,I)=YP(1,J)
217 YP(2,I)=YP(2,J)
218 YP(3,I)=YP(3,J)
219 YP(1,J)=Y1
220 YP(2,J)=Y2
221 YP(3,J)=Y3
22212 CONTINUE
223C
224C******************************
22514 YT(1,1)=0.0
226 YT(2,1)=0.0
227 YT(3,1)=0.0
228 IF(IHNT2(3).LE.1) GO TO 24
229 DO 20 KT=1,IHNT2(3)
23015 R=HIRND(2)
231 X=RLU_HIJING(0)
232 CX=2.0*X-1.0
233 SX=SQRT(1.0-CX*CX)
234C ********choose theta from uniform cos(theta) distr
235 PHI=RLU_HIJING(0)*2.0*HIPR1(40)
236C ********chose phi form uniform phi distr 0 to 2*pi
237 YT(1,KT)=R*SX*COS(PHI)
238 YT(2,KT)=R*SX*SIN(PHI)
239 YT(3,KT)=R*CX
240 IF(HIPR1(29).EQ.0.0) GO TO 20
241 DO 18 KT2=1,KT-1
242 DNBT1=(YT(1,KT)-YT(1,KT2))**2
243 DNBT2=(YT(2,KT)-YT(2,KT2))**2
244 DNBT3=(YT(3,KT)-YT(3,KT2))**2
245 DNBT=DNBT1+DNBT2+DNBT3
246 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
247C ********two neighbors cannot be closer than
248C HIPR1(29)
24918 CONTINUE
25020 CONTINUE
251 DO 22 I=1,IHNT2(3)-1
252 DO 22 J=I+1,IHNT2(3)
253 IF(YT(3,I).LT.YT(3,J)) GO TO 22
254 Y1=YT(1,I)
255 Y2=YT(2,I)
256 Y3=YT(3,I)
257 YT(1,I)=YT(1,J)
258 YT(2,I)=YT(2,J)
259 YT(3,I)=YT(3,J)
260 YT(1,J)=Y1
261 YT(2,J)=Y2
262 YT(3,J)=Y3
26322 CONTINUE
264
265C********************
26624 MISS=-1
267
26850 MISS=MISS+1
269 IF(MISS.GT.50) THEN
270 WRITE(6,*) 'infinite loop happened in HIJING'
271 STOP
272 ENDIF
273
274 NATT=0
275 JATT=0
276 EATT=0.0
277 CALL HIJINI
278 NLOP=0
279C ********Initialize for a new event
28060 NT=0
281 NP=0
282 N0=0
283 N01=0
284 N10=0
285 N11=0
8e529af3 286 NELT=0
287 NINT=0
288 NELP=0
289 NINP=0
e74335a4 290 NSG=0
291 NCOLT=0
292
293C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
294C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
295C BY THE ANGLE PHI FOR EACH COLLISION.******************
296C
297 BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
298 PHI=2.0*HIPR1(40)*RLU_HIJING(0)
299 BBX=BB*COS(PHI)
300 BBY=BB*SIN(PHI)
301 HINT1(19)=BB
302 HINT1(20)=PHI
303C
304 DO 70 JP=1,IHNT2(1)
305 DO 70 JT=1,IHNT2(3)
306 SCIP(JP,JT)=-1.0
307 B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
308 R2=B2*HIPR1(40)/HIPR1(31)/0.1
309C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
310 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
311 & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
312 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
313 & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
314 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
315 & *SQRT(1.0-RRB1)
316 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
317 & *SQRT(1.0-RRB2)
318 HINT1(18)=HINT1(14)-APHX1*HINT1(15)
319 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
320 IF(IHPR2(14).EQ.0.OR.
321 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
322 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
323 RANTOT=RLU_HIJING(0)
324 IF(RANTOT.GT.GS) GO TO 70
325 GO TO 65
326 ENDIF
327 GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
328 & /HIPR1(31)/2.0*ROMG(0.0)))
329 R2=R2/GSTOT_0
330 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
331 GSTOT=2.0*(1.0-SQRT(1.0-GS))
332 RANTOT=RLU_HIJING(0)*GSTOT_0
333 IF(RANTOT.GT.GSTOT) GO TO 70
334 IF(RANTOT.GT.GS) THEN
335 CALL HIJCSC(JP,JT)
336 GO TO 70
337C ********perform elastic collisions
338 ENDIF
339 65 SCIP(JP,JT)=R2
340 RNIP(JP,JT)=RANTOT
341 SJIP(JP,JT)=HINT1(18)
342 NCOLT=NCOLT+1
343 IPCOL(NCOLT)=JP
344 ITCOL(NCOLT)=JT
34570 CONTINUE
346C ********total number interactions proj and targ has
347C suffered
348 IF(NCOLT.EQ.0) THEN
349 NLOP=NLOP+1
350 IF(NLOP.LE.20.OR.
351 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
46f33c91 352
353
e74335a4 354 RETURN
355 ENDIF
356C ********At large impact parameter, there maybe no
357C interaction at all. For NN collision
358C repeat the event until interaction happens
359C
360 IF(IHPR2(3).NE.0) THEN
361 NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
362 NHARD=MIN(NHARD,NCOLT)
363 JPHARD=IPCOL(NHARD)
364 JTHARD=ITCOL(NHARD)
365 ENDIF
366C
367 IF(IHPR2(9).EQ.1) THEN
368 NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
369 NMINI=MIN(NMINI,NCOLT)
370 JPMINI=IPCOL(NMINI)
371 JTMINI=ITCOL(NMINI)
372 ENDIF
373C ********Specifying the location of the hard and
374C minijet if they are enforced by user
375C
376 DO 200 JP=1,IHNT2(1)
377 DO 200 JT=1,IHNT2(3)
378 IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
379 NFP(JP,11)=NFP(JP,11)+1
380 NFT(JT,11)=NFT(JT,11)+1
381 IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
382 NP=NP+1
383 N01=N01+1
384 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
385 NT=NT+1
386 N10=N10+1
387 ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
388 NP=NP+1
389 NT=NT+1
390 N0=N0+1
391 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
392 N11=N11+1
393 ENDIF
394 JOUT=0
395 NFP(JP,10)=0
396 NFT(JT,10)=0
397C*****************************************************************
398 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
399C ********When IHPR2(8)=0 no jets are produced
400 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
401C ********jets can not be produced for (JP,JT)
402C because not enough energy avaible for
403C JP or JT
404 R2=SCIP(JP,JT)
405 HINT1(18)=SJIP(JP,JT)
406 TT=ROMG(R2)*HINT1(18)/HIPR1(31)
407 TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
408 NJET=0
409 IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
410 CALL JETINI(JP,JT,1)
411 CALL HIJHRD(JP,JT,0,JFLG,0)
412 HINT1(26)=HINT1(47)
413 HINT1(27)=HINT1(48)
414 HINT1(28)=HINT1(49)
415 HINT1(29)=HINT1(50)
416 HINT1(36)=HINT1(67)
417 HINT1(37)=HINT1(68)
418 HINT1(38)=HINT1(69)
419 HINT1(39)=HINT1(70)
420C
421 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
422 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
423 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
424 & JFLG.GE.3) IASG(NSG,3)=1
425 IHNT2(9)=IHNT2(14)
426 IHNT2(10)=IHNT2(15)
427 DO 105 I05=1,5
428 HINT1(20+I05)=HINT1(40+I05)
429 HINT1(30+I05)=HINT1(50+I05)
430 105 CONTINUE
431 JOUT=1
432 IF(IHPR2(8).EQ.0) GO TO 160
433 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
434 & /REAL(IHNT2(1))**0.6666667,1.0)
435 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
436 & /REAL(IHNT2(3))**0.6666667,1.0)
437 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
438 & *SQRT(1.0-RRB1)
439 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
440 & *SQRT(1.0-RRB2)
441 HINT1(65)=HINT1(61)-APHX1*HINT1(62)
442 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
443 TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
444 NJET=-1
445C ********subtract the trigger jet from total number
446C of jet production to be done since it has
447C already been produced here
448 XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
449 106 NJET=NJET+1
450 XR1=XR1-ALOG(RLU_HIJING(0))
451 IF(XR1.LT.TTRIG) GO TO 106
452 XR=0.0
453 107 NJET=NJET+1
454 XR=XR-ALOG(RLU_HIJING(0))
455 IF(XR.LT.TT-TTRIG) GO TO 107
456 NJET=NJET-1
457 GO TO 112
458 ENDIF
459C ********create a hard interaction with specified P_T
460c when IHPR2(3)>0
461 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
462C ********create at least one pair of mini jets
463C when IHPR2(9)=1
464C
465 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
466 & (1.0-EXP(-TTS))) GO TO 160
467C ********this is the probability for no jet production
468110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
469111 NJET=NJET+1
470 XR=XR-ALOG(RLU_HIJING(0))
471 IF(XR.LT.TT) GO TO 111
472112 NJET=MIN(NJET,IHPR2(8))
473 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
474C ******** Determine number of mini jet production
475C
476 DO 150 I_JET=1,NJET
477 CALL JETINI(JP,JT,0)
478 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
479C ********JFLG=1 jets valence quarks, JFLG=2 with
480C gluon jet, JFLG=3 with q-qbar prod for
481C (JP,JT). If JFLG=0 jets can not be produced
482C this time. If JFLG=-1, error occured abandon
483C this event. JOUT is the total hard scat for
484C (JP,JT) up to now.
485 IF(JFLG.EQ.0) GO TO 160
486 IF(JFLG.LT.0) THEN
487 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
488 GO TO 50
489 ENDIF
490 JOUT=JOUT+1
491 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
492 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
493 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
494 & JFLG.GE.3) IASG(NSG,3)=1
495C ******** jet with PT>HIPR1(11) will be quenched
496 150 CONTINUE
497 160 CONTINUE
498 CALL HIJSFT(JP,JT,JOUT,IERROR)
499 IF(IERROR.NE.0) THEN
500 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
501 GO TO 50
502 ENDIF
503C
504C ********conduct soft scattering between JP and JT
505 JATT=JATT+JOUT
506
507200 CONTINUE
508
8e529af3 509c
510c**************************
511c
512 DO 201 JP=1,IHNT2(1)
513 IF(NFP(JP,5).GT.2) THEN
514 NINP=NINP+1
515 ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
516 NELP=NELP+1
517 ENDIF
518 201 continue
519 DO 202 JT=1,IHNT2(3)
520 IF(NFT(JT,5).GT.2) THEN
521 NINT=NINT+1
522 ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
523 NELT=NELT+1
524 ENDIF
525 202 continue
526c
527c*******************************
528
e74335a4 529C********perform jet quenching for jets with PT>HIPR1(11)**********
530
531 IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
532 & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
533 DO 271 I=1,IHNT2(1)
534 IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
535271 CONTINUE
536 DO 272 I=1,IHNT2(3)
537 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
538272 CONTINUE
539 DO 273 ISG=1,NSG
540 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
541273 CONTINUE
542 ENDIF
543C
544C**************fragment all the string systems in the following*****
545C
546C********N_ST is where particle information starts
547C********N_STR+1 is the number of strings in fragmentation
548C********the number of strings before a line is stored in K(I,4)
549C********IDSTR is id number of the string system (91,92 or 93)
550C
551 IF(IHPR2(20).NE.0) THEN
552 DO 360 ISG=1,NSG
553 CALL HIJFRG(ISG,3,IERROR)
554 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
555 MSTU(24)=0
556 MSTU(28)=0
557 IF(IHPR2(10).NE.0) THEN
558 call LULIST_HIJING(1)
559 WRITE(6,*) 'error occured, repeat the event'
560 ENDIF
561 GO TO 50
562 ENDIF
563C ********Check errors
564C
565 N_ST=1
566 IDSTR=92
567 IF(IHPR2(21).EQ.0) THEN
568 CALL LUEDIT_HIJING(2)
569 ELSE
570351 N_ST=N_ST+1
571 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351
572 IDSTR=K(N_ST,2)
573 N_ST=N_ST+1
574 ENDIF
575C
576 IF(FRAME.EQ.'LAB') THEN
577 CALL HIBOOST
578 ENDIF
579C ******** boost back to lab frame(if it was in)
580C
581 N_STR=0
582 DO 360 I=N_ST,N
583 IF(K(I,2).EQ.IDSTR) THEN
584 N_STR=N_STR+1
585 GO TO 360
586 ENDIF
587 K(I,4)=N_STR
588 NATT=NATT+1
589 KATT(NATT,1)=K(I,2)
590 KATT(NATT,2)=20
591 KATT(NATT,4)=K(I,1)
592 IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
593 KATT(NATT,3)=0
594 ELSE
595 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
596 ENDIF
597C ****** identify the mother particle
598 PATT(NATT,1)=P(I,1)
599 PATT(NATT,2)=P(I,2)
600 PATT(NATT,3)=P(I,3)
601 PATT(NATT,4)=P(I,4)
2818cab6 602 VATT(NATT,1)=V(I,1)
603 VATT(NATT,2)=V(I,2)
604 VATT(NATT,3)=V(I,3)
605 VATT(NATT,4)=V(I,4)
606
e74335a4 607 EATT=EATT+P(I,4)
608360 CONTINUE
609C ********Fragment the q-qbar jets systems *****
610C
611 JTP(1)=IHNT2(1)
612 JTP(2)=IHNT2(3)
613 DO 400 NTP=1,2
614 DO 400 J_JTP=1,JTP(NTP)
615 CALL HIJFRG(J_JTP,NTP,IERROR)
616 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
617 MSTU(24)=0
618 MSTU(28)=0
619 IF(IHPR2(10).NE.0) THEN
620 call LULIST_HIJING(1)
621 WRITE(6,*) 'error occured, repeat the event'
622 ENDIF
623 GO TO 50
624 ENDIF
625C ********check errors
626C
627 N_ST=1
628 IDSTR=92
46f33c91 629
630 NFTP=NFP(J_JTP,5)
631 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
632
e74335a4 633 IF(IHPR2(21).EQ.0) THEN
634 CALL LUEDIT_HIJING(2)
46f33c91 635 ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN
e74335a4 636381 N_ST=N_ST+1
637 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
638 IDSTR=K(N_ST,2)
639 N_ST=N_ST+1
46f33c91 640 ENDIF
e74335a4 641 IF(FRAME.EQ.'LAB') THEN
642 CALL HIBOOST
643 ENDIF
644C ******** boost back to lab frame(if it was in)
645C
e74335a4 646 N_STR=0
647 DO 390 I=N_ST,N
648 IF(K(I,2).EQ.IDSTR) THEN
649 N_STR=N_STR+1
650 GO TO 390
651 ENDIF
652 K(I,4)=N_STR
653 NATT=NATT+1
654 KATT(NATT,1)=K(I,2)
655 KATT(NATT,2)=NFTP
656 KATT(NATT,4)=K(I,1)
3bc18c05 657 IF(K(I,3).EQ.0 ) THEN
e74335a4 658 KATT(NATT,3)=0
3bc18c05 659 ELSE
660 IF(K(K(I,3),2).EQ.IDSTR) THEN
661 KATT(NATT,3)=0
662 ELSE
663 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
664 ENDIF
e74335a4 665 ENDIF
666C ****** identify the mother particle
667 PATT(NATT,1)=P(I,1)
668 PATT(NATT,2)=P(I,2)
669 PATT(NATT,3)=P(I,3)
670 PATT(NATT,4)=P(I,4)
2818cab6 671 VATT(NATT,1)=V(I,1)
672 VATT(NATT,2)=V(I,2)
673 VATT(NATT,3)=V(I,3)
674 VATT(NATT,4)=V(I,4)
675
e74335a4 676 EATT=EATT+P(I,4)
677390 CONTINUE
678400 CONTINUE
679C ********Fragment the q-qq related string systems
680 ENDIF
681
682 DO 450 I=1,NDR
683 NATT=NATT+1
684 KATT(NATT,1)=KFDR(I)
685 KATT(NATT,2)=40
686 KATT(NATT,3)=0
687 PATT(NATT,1)=PDR(I,1)
688 PATT(NATT,2)=PDR(I,2)
689 PATT(NATT,3)=PDR(I,3)
690 PATT(NATT,4)=PDR(I,4)
2818cab6 691 VATT(NATT,1)=V(I,1)
692 VATT(NATT,2)=V(I,2)
693 VATT(NATT,3)=V(I,3)
694 VATT(NATT,4)=V(I,4)
695
e74335a4 696 EATT=EATT+PDR(I,4)
697450 CONTINUE
698C ********store the direct-produced particles
699C
700 DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
701 IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
702 & .AND.IHPR2(21).EQ.0) THEN
703 IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
704 & event'
705C call LULIST_HIJING(1)
706 GO TO 50
707 ENDIF
46f33c91 708
e74335a4 709 RETURN
710 END