Include VATT to store vertex off particles.
[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
346 RETURN
347 ENDIF
348C ********At large impact parameter, there maybe no
349C interaction at all. For NN collision
350C repeat the event until interaction happens
351C
352 IF(IHPR2(3).NE.0) THEN
353 NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
354 NHARD=MIN(NHARD,NCOLT)
355 JPHARD=IPCOL(NHARD)
356 JTHARD=ITCOL(NHARD)
357 ENDIF
358C
359 IF(IHPR2(9).EQ.1) THEN
360 NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
361 NMINI=MIN(NMINI,NCOLT)
362 JPMINI=IPCOL(NMINI)
363 JTMINI=ITCOL(NMINI)
364 ENDIF
365C ********Specifying the location of the hard and
366C minijet if they are enforced by user
367C
368 DO 200 JP=1,IHNT2(1)
369 DO 200 JT=1,IHNT2(3)
370 IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
371 NFP(JP,11)=NFP(JP,11)+1
372 NFT(JT,11)=NFT(JT,11)+1
373 IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
374 NP=NP+1
375 N01=N01+1
376 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
377 NT=NT+1
378 N10=N10+1
379 ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
380 NP=NP+1
381 NT=NT+1
382 N0=N0+1
383 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
384 N11=N11+1
385 ENDIF
386 JOUT=0
387 NFP(JP,10)=0
388 NFT(JT,10)=0
389C*****************************************************************
390 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
391C ********When IHPR2(8)=0 no jets are produced
392 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
393C ********jets can not be produced for (JP,JT)
394C because not enough energy avaible for
395C JP or JT
396 R2=SCIP(JP,JT)
397 HINT1(18)=SJIP(JP,JT)
398 TT=ROMG(R2)*HINT1(18)/HIPR1(31)
399 TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
400 NJET=0
401 IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
402 CALL JETINI(JP,JT,1)
403 CALL HIJHRD(JP,JT,0,JFLG,0)
404 HINT1(26)=HINT1(47)
405 HINT1(27)=HINT1(48)
406 HINT1(28)=HINT1(49)
407 HINT1(29)=HINT1(50)
408 HINT1(36)=HINT1(67)
409 HINT1(37)=HINT1(68)
410 HINT1(38)=HINT1(69)
411 HINT1(39)=HINT1(70)
412C
413 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
414 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
415 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
416 & JFLG.GE.3) IASG(NSG,3)=1
417 IHNT2(9)=IHNT2(14)
418 IHNT2(10)=IHNT2(15)
419 DO 105 I05=1,5
420 HINT1(20+I05)=HINT1(40+I05)
421 HINT1(30+I05)=HINT1(50+I05)
422 105 CONTINUE
423 JOUT=1
424 IF(IHPR2(8).EQ.0) GO TO 160
425 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
426 & /REAL(IHNT2(1))**0.6666667,1.0)
427 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
428 & /REAL(IHNT2(3))**0.6666667,1.0)
429 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
430 & *SQRT(1.0-RRB1)
431 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
432 & *SQRT(1.0-RRB2)
433 HINT1(65)=HINT1(61)-APHX1*HINT1(62)
434 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
435 TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
436 NJET=-1
437C ********subtract the trigger jet from total number
438C of jet production to be done since it has
439C already been produced here
440 XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
441 106 NJET=NJET+1
442 XR1=XR1-ALOG(RLU_HIJING(0))
443 IF(XR1.LT.TTRIG) GO TO 106
444 XR=0.0
445 107 NJET=NJET+1
446 XR=XR-ALOG(RLU_HIJING(0))
447 IF(XR.LT.TT-TTRIG) GO TO 107
448 NJET=NJET-1
449 GO TO 112
450 ENDIF
451C ********create a hard interaction with specified P_T
452c when IHPR2(3)>0
453 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
454C ********create at least one pair of mini jets
455C when IHPR2(9)=1
456C
457 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
458 & (1.0-EXP(-TTS))) GO TO 160
459C ********this is the probability for no jet production
460110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
461111 NJET=NJET+1
462 XR=XR-ALOG(RLU_HIJING(0))
463 IF(XR.LT.TT) GO TO 111
464112 NJET=MIN(NJET,IHPR2(8))
465 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
466C ******** Determine number of mini jet production
467C
468 DO 150 I_JET=1,NJET
469 CALL JETINI(JP,JT,0)
470 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
471C ********JFLG=1 jets valence quarks, JFLG=2 with
472C gluon jet, JFLG=3 with q-qbar prod for
473C (JP,JT). If JFLG=0 jets can not be produced
474C this time. If JFLG=-1, error occured abandon
475C this event. JOUT is the total hard scat for
476C (JP,JT) up to now.
477 IF(JFLG.EQ.0) GO TO 160
478 IF(JFLG.LT.0) THEN
479 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
480 GO TO 50
481 ENDIF
482 JOUT=JOUT+1
483 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
484 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
485 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
486 & JFLG.GE.3) IASG(NSG,3)=1
487C ******** jet with PT>HIPR1(11) will be quenched
488 150 CONTINUE
489 160 CONTINUE
490 CALL HIJSFT(JP,JT,JOUT,IERROR)
491 IF(IERROR.NE.0) THEN
492 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
493 GO TO 50
494 ENDIF
495C
496C ********conduct soft scattering between JP and JT
497 JATT=JATT+JOUT
498
499200 CONTINUE
500
501C********perform jet quenching for jets with PT>HIPR1(11)**********
502
503 IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
504 & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
505 DO 271 I=1,IHNT2(1)
506 IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
507271 CONTINUE
508 DO 272 I=1,IHNT2(3)
509 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
510272 CONTINUE
511 DO 273 ISG=1,NSG
512 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
513273 CONTINUE
514 ENDIF
515C
516C**************fragment all the string systems in the following*****
517C
518C********N_ST is where particle information starts
519C********N_STR+1 is the number of strings in fragmentation
520C********the number of strings before a line is stored in K(I,4)
521C********IDSTR is id number of the string system (91,92 or 93)
522C
523 IF(IHPR2(20).NE.0) THEN
524 DO 360 ISG=1,NSG
525 CALL HIJFRG(ISG,3,IERROR)
526 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
527 MSTU(24)=0
528 MSTU(28)=0
529 IF(IHPR2(10).NE.0) THEN
530 call LULIST_HIJING(1)
531 WRITE(6,*) 'error occured, repeat the event'
532 ENDIF
533 GO TO 50
534 ENDIF
535C ********Check errors
536C
537 N_ST=1
538 IDSTR=92
539 IF(IHPR2(21).EQ.0) THEN
540 CALL LUEDIT_HIJING(2)
541 ELSE
542351 N_ST=N_ST+1
543 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351
544 IDSTR=K(N_ST,2)
545 N_ST=N_ST+1
546 ENDIF
547C
548 IF(FRAME.EQ.'LAB') THEN
549 CALL HIBOOST
550 ENDIF
551C ******** boost back to lab frame(if it was in)
552C
553 N_STR=0
554 DO 360 I=N_ST,N
555 IF(K(I,2).EQ.IDSTR) THEN
556 N_STR=N_STR+1
557 GO TO 360
558 ENDIF
559 K(I,4)=N_STR
560 NATT=NATT+1
561 KATT(NATT,1)=K(I,2)
562 KATT(NATT,2)=20
563 KATT(NATT,4)=K(I,1)
564 IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
565 KATT(NATT,3)=0
566 ELSE
567 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
568 ENDIF
569C ****** identify the mother particle
570 PATT(NATT,1)=P(I,1)
571 PATT(NATT,2)=P(I,2)
572 PATT(NATT,3)=P(I,3)
573 PATT(NATT,4)=P(I,4)
574 EATT=EATT+P(I,4)
575360 CONTINUE
576C ********Fragment the q-qbar jets systems *****
577C
578 JTP(1)=IHNT2(1)
579 JTP(2)=IHNT2(3)
580 DO 400 NTP=1,2
581 DO 400 J_JTP=1,JTP(NTP)
582 CALL HIJFRG(J_JTP,NTP,IERROR)
583 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
584 MSTU(24)=0
585 MSTU(28)=0
586 IF(IHPR2(10).NE.0) THEN
587 call LULIST_HIJING(1)
588 WRITE(6,*) 'error occured, repeat the event'
589 ENDIF
590 GO TO 50
591 ENDIF
592C ********check errors
593C
594 N_ST=1
595 IDSTR=92
596 IF(IHPR2(21).EQ.0) THEN
597 CALL LUEDIT_HIJING(2)
598 ELSE
599381 N_ST=N_ST+1
600 IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
601 IDSTR=K(N_ST,2)
602 N_ST=N_ST+1
603 ENDIF
604 IF(FRAME.EQ.'LAB') THEN
605 CALL HIBOOST
606 ENDIF
607C ******** boost back to lab frame(if it was in)
608C
609 NFTP=NFP(J_JTP,5)
610 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
611 N_STR=0
612 DO 390 I=N_ST,N
613 IF(K(I,2).EQ.IDSTR) THEN
614 N_STR=N_STR+1
615 GO TO 390
616 ENDIF
617 K(I,4)=N_STR
618 NATT=NATT+1
619 KATT(NATT,1)=K(I,2)
620 KATT(NATT,2)=NFTP
621 KATT(NATT,4)=K(I,1)
622 IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
623 KATT(NATT,3)=0
624 ELSE
625 KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
626 ENDIF
627C ****** identify the mother particle
628 PATT(NATT,1)=P(I,1)
629 PATT(NATT,2)=P(I,2)
630 PATT(NATT,3)=P(I,3)
631 PATT(NATT,4)=P(I,4)
632 EATT=EATT+P(I,4)
633390 CONTINUE
634400 CONTINUE
635C ********Fragment the q-qq related string systems
636 ENDIF
637
638 DO 450 I=1,NDR
639 NATT=NATT+1
640 KATT(NATT,1)=KFDR(I)
641 KATT(NATT,2)=40
642 KATT(NATT,3)=0
643 PATT(NATT,1)=PDR(I,1)
644 PATT(NATT,2)=PDR(I,2)
645 PATT(NATT,3)=PDR(I,3)
646 PATT(NATT,4)=PDR(I,4)
647 EATT=EATT+PDR(I,4)
648450 CONTINUE
649C ********store the direct-produced particles
650C
651 DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
652 IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
653 & .AND.IHPR2(21).EQ.0) THEN
654 IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
655 & event'
656C call LULIST_HIJING(1)
657 GO TO 50
658 ENDIF
659 RETURN
660 END