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