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