f8a46fafd5a953fb84a214d577506eb6b7b8e99d
[u/mrichter/AliRoot.git] / TAmpt / AMPT / hijing1.383_ampt.f
1 c.................... hijing1.383_ampt.f
2 c     Version 1.383
3 c     The variables isng in HIJSFT and JL in ATTRAD were not initialized.
4 c     The version initialize them. (as found by Fernando Marroquim)
5 c
6 c
7 c
8 c     Version 1.382
9 c     Nuclear distribution for deuteron is taken as the Hulthen wave
10 c     function as provided by Brian Cole (Columbia)
11 clin     used my own implementation of impact parameter 
12 clin     & proton-neutron distance within a deuteron.
13 c
14 c
15 c     Version 1.381
16 c
17 c     The parameters for Wood-Saxon distribution for deuteron are
18 c     constrained to give the right rms ratius 2.116 fm
19 c     (R=0.0, D=0.5882)
20 c
21 c
22 c     Version 1.38
23 c
24 c     The following common block is added to record the number of elastic
25 c     (NELT, NELP) and inelastic (NINT, NINP) participants
26 c
27 c        COMMON/HJGLBR/NELT,NINT,NELP,NINP
28 c        SAVE /HJGLBR/
29 c
30 c     Version 1.37
31 c
32 c     A bug in the quenching subroutine is corrected. When calculating the
33 c     distance between two wounded nucleons, the displacement of the
34 c     impact parameter was not inculded. This bug was discovered by
35 c     Dr. V.Uzhinskii JINR, Dubna, Russia
36 c
37 c
38 C     Version 1.36
39 c
40 c     Modification Oct. 8, 1998. In hijing, log(ran(nseed)) occasionally
41 c     causes overfloat. It is modified to log(max(ran(nseed),1.0e-20)).
42 c
43 c
44 C     Nothing important has been changed here. A few 'garbage' has been
45 C     cleaned up here, like common block HJJET3 for the sea quark strings
46 C     which were originally created to implement the DPM scheme which
47 C     later was abadoned in the final version. The lines which operate
48 C     on these data are also deleted in the program.
49 C
50 C
51 C     Version 1.35
52 C     There are some changes in the program: subroutine HARDJET is now
53 C     consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
54 C     for the triggered hard processes. Now that is done  altogether
55 C     with other normal hard processes in modified JETINI. In the new
56 C     version one calls JETINI every time one calls HIJHRD. In the new
57 C     version the effect of the isospin of the nucleon on hard processes,
58 C     especially direct photons is correctly considered.
59 C     For A+A collisions, one has to initilize pythia
60 C     separately for each type of collisions, pp, pn,np and nn,
61 C     or hp and hn for hA collisions. In JETINI we use the following
62 C     catalogue for different types of collisions:
63 C     h+h: h+h (itype=1)
64 C     h+A: h+p (itype=1), h+n (itype=2)
65 C     A+h: p+h (itype=1), n+h (itype=2)
66 C     A+A: p+p (itype=1), p+n (itype=2), n+p (itype=3), n+n (itype=4)
67 C*****************************************************************
68 c
69 C
70 C     Version 1.34
71 C     Last modification on January 5, 1998. Two mistakes are corrected in
72 C     function G. A Mistake in the subroutine Parton is also corrected.
73 C     (These are pointed out by Ysushi Nara).
74 C
75 C
76 C       Last modifcation on April 10, 1996. To conduct final
77 C       state radiation, PYTHIA reorganize the two scattered
78 C       partons and their final momenta will be a little
79 C       different. The summed total momenta of the partons
80 C       from the final state radiation are stored in HINT1(26-29)
81 C       and HINT1(36-39) which are little different from 
82 C       HINT1(21-24) and HINT1(41-44).
83 C
84 C       Version 1.33
85 C
86 C       Last modfication  on September 11, 1995. When HIJING and
87 C       PYTHIA are initialized, the shadowing is evaluated at
88 C       b=0 which is the maximum. This will cause overestimate
89 C       of shadowing for peripheral interactions. To correct this
90 C       problem, shadowing is set to zero when initializing. Then
91 C       use these maximum  cross section without shadowing as a
92 C       normalization of the Monte Carlo. This however increase
93 C       the computing time. IHNT2(16) is used to indicate whether
94 C       the sturcture function is called for (IHNT2(16)=1) initialization
95 C       or for (IHNT2(16)=0)normal collisions simulation
96 C
97 C       Last modification on Aagust 28, 1994. Two bugs associate
98 C       with the impact parameter dependence of the shadowing is
99 C       corrected.
100 C
101 C
102 c       Last modification on October 14, 1994. One bug is corrected
103 c       in the direct photon production option in subroutine
104 C       HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
105 C       Another bug associated with keeping the decay history
106 C       in the particle information is also corrected.(this problem
107 C       was reported by Matt Bloomer)
108 C
109 C
110 C       Last modification on July 15, 1994. The option to trig on
111 C       heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1) 
112 C       is added. To do this, set IHPR2(3)=3. For inclusive production,
113 C       one should reset HIPR1(10)=0.0. One can also trig larger pt
114 C       QQbar production by giving HIPR1(10) a nonvanishing value.
115 C       The mass of the heavy quark in the calculation of the cross
116 C       section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
117 C       default is the charm mass D=1.5). We also include a separate
118 C       K-factor for heavy quark and direct photon production by
119 C       HIPR1(23)(D=2.0).
120 C
121 C       Last modification on May 24, 1994.  The option to
122 C       retain the information of all particles including those
123 C       who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is 
124 C       added to contain the line number of the parent particle 
125 C       of the current line which is produced via a decay. 
126 C       KATT(I,4) is the status number of the particle: 11=particle
127 C       which has decayed; 1=finally produced particle.
128 C
129 C
130 C       Last modification on May 24, 1994( in HIJSFT when valence quark
131 C       is quenched, the following error is corrected. 1.2*IHNT2(1) --> 
132 C       1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
133 C
134 C
135 C       Last modification on March 16, 1994 (heavy flavor production
136 C       processes MSUB(81)=1 MSUB(82)=1 have been switched on,
137 C       charm production is the default, B-quark option is
138 C       IHPR2(18), when it is switched on, charm quark is 
139 C       automatically off)
140 C
141 C
142 C       Last modification on March 23, 1994 (an error is corrected
143 C       in the impact parameter dependence of the jet cross section)
144 C
145 C       Last modification Oct. 1993 to comply with non-vax
146 C       machines' compiler 
147 C
148 C*********************************************
149 C       LAST MODIFICATION April 5, 1991
150 CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B 
151 C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
152 C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
153 C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18), 
154 C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
155 C
156 C       gluon jets can form a single string system
157 C
158 C       initial state radiation is included
159 C       
160 C       all QCD subprocesses are included
161 c
162 c       direct particles production is included(currently only direct
163 C               photon)
164 c
165 C       Effect of high P_T trigger bias on multiple jets distribution
166 c
167 C******************************************************************
168 C                               HIJING.10                         *
169 C                 Heavy Ion Jet INteraction Generator             *
170 C                                  by                             *
171 C                  X. N. Wang      and   M. Gyulassy              *
172 C                     Lawrence Berkeley Laboratory                *
173 C                                                                 *
174 C******************************************************************
175 C
176 C******************************************************************
177 C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of  *
178 C proj, NFP(K,4) original ID of proj.  NFP(K,5)=colli status(0=no,*
179 C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
180 C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it   *
181 C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
182 C (0=has not been,1=is going to be, -1=has already been scattered *
183 C NFP(k,11) total number of interactions this proj has suffered   *
184 C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant  *
185 C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and     *
186 C diquark,PP(K,10)=PT of the hard scattering between the valence  *
187 C quarks; PP(K,14,15)=the mass of quark,diquark.                  * 
188 C******************************************************************
189 C
190 C****************************************************************
191 C
192 C       SUBROUTINE HIJING
193 C
194 C****************************************************************
195         SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
196
197 cgsfs   Added following for consistency with AMPT call
198         double precision BMIN0, BMAX0
199
200 cbz1/25/99
201         PARAMETER (MAXPTN=400001)
202 clin-4/20/01        PARAMETER (MAXSTR = 1600)
203         PARAMETER (MAXSTR=150001)
204 cbz1/25/99end
205 clin-4/26/01:
206         PARAMETER (MAXIDL=4001)
207
208 cbz1/31/99
209         DOUBLE PRECISION  GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
210         DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
211         DOUBLE PRECISION  ATAUI, ZT1, ZT2, ZT3
212         DOUBLE PRECISION  xnprod,etprod,xnfrz,etfrz,
213      & dnprod,detpro,dnfrz,detfrz
214
215 cbz1/31/99end
216
217         CHARACTER FRAME*8
218         DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
219      &                        IPCOL(90000),ITCOL(90000)
220         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
221 cc      SAVE /HPARNT/
222 C
223         COMMON/hjcrdn/YP(3,300),YT(3,300)
224 cc      SAVE /hjcrdn/
225 clin-7/16/03 NINT is a intrinsic fortran function, rename it to NINTHJ
226 c        COMMON/HJGLBR/NELT,NINT,NELP,NINP
227         COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
228 cc      SAVE /HJGLBR/
229         COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
230 cc      SAVE /HMAIN1/
231 clin-4/26/01
232 c        COMMON/HMAIN2/KATT(130000,4),PATT(130000,4)
233         COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
234 cc      SAVE /HMAIN2/
235         COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
236 cc      SAVE /HSTRNG/
237         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
238      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
239      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
240      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
241      &                PJTE(300,500),PJTM(300,500)
242 cc      SAVE /HJJET1/
243 clin-4/2008
244 c        COMMON/HJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
245 c     &       K2SG(900,100),PXSG(900,100),PYSG(900,100),
246 c     &       PZSG(900,100),PESG(900,100),PMSG(900,100)
247         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
248      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
249      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
250 cc      SAVE /HJJET2/
251         COMMON/HJJET4/NDR,IADR(MAXSTR,2),KFDR(MAXSTR),PDR(MAXSTR,5)
252 clin-4/2008:
253 c        common/xydr/rtdr(900,2)
254         common/xydr/rtdr(MAXSTR,2)
255 cc      SAVE /HJJET4/
256       COMMON/RNDF77/NSEED
257 cc      SAVE /RNDF77/
258 C
259         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)   
260 cc      SAVE /LUJETSA/
261         COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
262 cc      SAVE /LUDAT1A/
263
264 clin-9/29/03 changed name in order to distinguish from /prec2/
265         COMMON /ARPRC/ ITYPAR(MAXSTR),
266      &       GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
267      &       PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
268      &       XMAR(MAXSTR)
269 ccbz11/11/98
270 c        COMMON /ARPRC/ ITYP(MAXSTR),
271 c     &     GX(MAXSTR), GY(MAXSTR), GZ(MAXSTR), FT(MAXSTR),
272 c     &     PX(MAXSTR), PY(MAXSTR), PZ(MAXSTR), EE(MAXSTR),
273 c     &     XM(MAXSTR)
274 cc      SAVE /ARPRC/
275 ccbz11/11/98end
276
277 cbz1/25/99
278         COMMON /PARA1/ MUL
279 cc      SAVE /PARA1/
280         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
281      &     PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
282      &     XMASS0(MAXPTN), ITYP0(MAXPTN)
283 cc      SAVE /prec1/
284         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
285      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
286      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
287 cc      SAVE /prec2/
288         COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
289 cc      SAVE /ilist7/
290         COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
291 cc      SAVE /ilist8/
292         COMMON /SREC1/ NSP, NST, NSI
293 cc      SAVE /SREC1/
294         COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
295 cc      SAVE /SREC2/
296 cbz1/25/99end
297
298 clin-2/25/00
299         COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
300      & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
301 cc      SAVE /frzout/ 
302 clin-4/11/01 soft:
303       common/anim/nevent,isoft,isflag,izpc
304 cc      SAVE /anim/
305 clin-4/25/01 soft3:
306       DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
307      1     GXSGS,GYSGS,GZSGS,FTSGS
308       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
309      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
310      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
311      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
312 cc      SAVE /SOFT/
313 clin-4/26/01 lepton and photon info:
314         COMMON /NOPREC/ NNOZPC, ITYPN(MAXIDL),
315      &       GXN(MAXIDL), GYN(MAXIDL), GZN(MAXIDL), FTN(MAXIDL),
316      &       PXN(MAXIDL), PYN(MAXIDL), PZN(MAXIDL), EEN(MAXIDL),
317      &       XMN(MAXIDL)
318 cc      SAVE /NOPREC/
319 clin-6/22/01:
320         common /lastt/itimeh,bimp
321 cc      SAVE /lastt/
322         COMMON /AREVT/ IAEVT, IARUN, MISS
323         common/phidcy/iphidcy,pttrig,ntrig,maxmiss
324 cwei        DOUBLE PRECISION PATT
325         SAVE   
326
327 cgsfs      WRITE(*,*) "IN Hijing, FRAME=",FRAME
328 cgsfs      WRITE(*,*) "IN Hijing, BMIN=",BMIN0
329 cgsfs      WRITE(*,*) "IN Hijing, BMAX=",BMAX0
330
331         BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
332         BMIN=MIN(BMIN0,BMAX)
333         IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
334                 BMIN=0.0
335                 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
336         ENDIF
337 C                        ********HIPR1(31) is in mb =0.1fm**2
338 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS 
339 C       BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
340 C
341 cgsfs      WRITE(*,*) "IN Hijing, Modified BMIN=",BMIN
342 cgsfs      WRITE(*,*) "IN Hijing, Modified BMAX=",BMAX
343         YP(1,1)=0.0
344         YP(2,1)=0.0
345         YP(3,1)=0.0
346         IF(IHNT2(1).LE.1) GO TO 14
347         DO 10 KP=1,IHNT2(1)
348 5        R=HIRND(1)
349         X=RANART(NSEED)
350         CX=2.0*X-1.0
351         SX=SQRT(1.0-CX*CX)
352 C                ********choose theta from uniform cos(theta) distr
353         PHI=RANART(NSEED)*2.0*HIPR1(40)
354 C                ********choose phi form uniform phi distr 0 to 2*pi
355         YP(1,KP)=R*SX*COS(PHI)
356         YP(2,KP)=R*SX*SIN(PHI)
357         YP(3,KP)=R*CX
358         IF(HIPR1(29).EQ.0.0) GO TO 10
359         DO 8  KP2=1,KP-1
360                 DNBP1=(YP(1,KP)-YP(1,KP2))**2
361                 DNBP2=(YP(2,KP)-YP(2,KP2))**2
362                 DNBP3=(YP(3,KP)-YP(3,KP2))**2
363                 DNBP=DNBP1+DNBP2+DNBP3
364                 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
365 C                        ********two neighbors cannot be closer than 
366 C                                HIPR1(29)
367 8        CONTINUE
368 10        CONTINUE
369
370 clin-1/27/03 Hulthen wavefn for deuteron borrowed from hijing1.382.f, 
371 c     but modified [divide by 2, & x(p)=-x(n)]: 
372 c     (Note: hijing1.383.f has corrected this bug in hijing1.382.f)
373         if(IHNT2(1).EQ.2) then
374            rnd1=max(RANART(NSEED),1.0e-20)
375            rnd2=max(RANART(NSEED),1.0e-20)
376            rnd3=max(RANART(NSEED),1.0e-20)
377            R=-(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
378      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
379            X=RANART(NSEED)
380            CX=2.0*X-1.0
381            SX=SQRT(1.0-CX*CX)
382            PHI=RANART(NSEED)*2.0*HIPR1(40)
383 c     R above is the relative distance between p & n in a deuteron:
384            R=R/2.
385            YP(1,1)=R*SX*COS(PHI)
386            YP(2,1)=R*SX*SIN(PHI)
387            YP(3,1)=R*CX
388 c     p & n has opposite coordinates in the deuteron frame:
389            YP(1,2)=-YP(1,1)
390            YP(2,2)=-YP(2,1)
391            YP(3,2)=-YP(3,1)
392         endif
393
394         DO 12 I=1,IHNT2(1)-1
395         DO 12 J=I+1,IHNT2(1)
396         IF(YP(3,I).GT.YP(3,J)) GO TO 12
397         Y1=YP(1,I)
398         Y2=YP(2,I)
399         Y3=YP(3,I)
400         YP(1,I)=YP(1,J)
401         YP(2,I)=YP(2,J)
402         YP(3,I)=YP(3,J)
403         YP(1,J)=Y1
404         YP(2,J)=Y2
405         YP(3,J)=Y3
406 12        CONTINUE
407 C
408 C******************************
409 14        YT(1,1)=0.0
410         YT(2,1)=0.0
411         YT(3,1)=0.0
412         IF(IHNT2(3).LE.1) GO TO 24
413         DO 20 KT=1,IHNT2(3)
414 15        R=HIRND(2)
415         X=RANART(NSEED)
416         CX=2.0*X-1.0
417         SX=SQRT(1.0-CX*CX)
418 C                ********choose theta from uniform cos(theta) distr
419         PHI=RANART(NSEED)*2.0*HIPR1(40)
420 C                ********chose phi form uniform phi distr 0 to 2*pi
421         YT(1,KT)=R*SX*COS(PHI)
422         YT(2,KT)=R*SX*SIN(PHI)
423         YT(3,KT)=R*CX
424         IF(HIPR1(29).EQ.0.0) GO TO 20
425         DO 18  KT2=1,KT-1
426                 DNBT1=(YT(1,KT)-YT(1,KT2))**2
427                 DNBT2=(YT(2,KT)-YT(2,KT2))**2
428                 DNBT3=(YT(3,KT)-YT(3,KT2))**2
429                 DNBT=DNBT1+DNBT2+DNBT3
430                 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
431 C                        ********two neighbors cannot be closer than 
432 C                                HIPR1(29)
433 18        CONTINUE
434 20        CONTINUE
435 c
436 clin-1/27/03 Hulthen wavefn for deuteron borrowed from hijing1.382.f, 
437 c     but modified [divide by 2, & x(p)=-x(n)]:
438         if(IHNT2(3).EQ.2) then
439            rnd1=max(RANART(NSEED),1.0e-20)
440            rnd2=max(RANART(NSEED),1.0e-20)
441            rnd3=max(RANART(NSEED),1.0e-20)
442            R=-(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
443      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
444            X=RANART(NSEED)
445            CX=2.0*X-1.0
446            SX=SQRT(1.0-CX*CX)
447            PHI=RANART(NSEED)*2.0*HIPR1(40)
448            R=R/2.
449            YT(1,1)=R*SX*COS(PHI)
450            YT(2,1)=R*SX*SIN(PHI)
451            YT(3,1)=R*CX
452            YT(1,2)=-YT(1,1)
453            YT(2,2)=-YT(2,1)
454            YT(3,2)=-YT(3,1)
455         endif
456 c
457         DO 22 I=1,IHNT2(3)-1
458         DO 22 J=I+1,IHNT2(3)
459         IF(YT(3,I).LT.YT(3,J)) GO TO 22
460         Y1=YT(1,I)
461         Y2=YT(2,I)
462         Y3=YT(3,I)
463         YT(1,I)=YT(1,J)
464         YT(2,I)=YT(2,J)
465         YT(3,I)=YT(3,J)
466         YT(1,J)=Y1
467         YT(2,J)=Y2
468         YT(3,J)=Y3
469 22        CONTINUE
470
471 C********************
472 24        MISS=-1
473 50        MISS=MISS+1
474
475 clin-6/2009 ctest on
476 c        IF(MISS.GT.50) THEN
477         IF(MISS.GT.maxmiss) THEN
478            WRITE(6,*) 'infinite loop happened in  HIJING'
479            STOP
480         ENDIF
481
482 clin-4/30/01:
483         itest=0
484
485         NATT=0
486         JATT=0
487         EATT=0.0
488         CALL HIJINI
489         NLOP=0
490 C                        ********Initialize for a new event
491 60        NT=0
492         NP=0
493         N0=0
494         N01=0
495         N10=0
496         N11=0
497         NELT=0
498         NINTHJ=0
499         NELP=0
500         NINP=0
501         NSG=0
502         NCOLT=0
503
504 C****        BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS 
505 C       RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET 
506 C       BY THE ANGLE PHI  FOR EACH COLLISION.******************
507 C
508         BB=SQRT(BMIN**2+RANART(NSEED)*(BMAX**2-BMIN**2))
509 cbz6/28/99 flow1
510         PHI=2.0*HIPR1(40)*RANART(NSEED)
511 c        PHI=0.
512 cbz6/28/99 flow1 end
513         BBX=BB*COS(PHI)
514         BBY=BB*SIN(PHI)
515         HINT1(19)=BB
516         HINT1(20)=PHI
517 C
518         DO 70 JP=1,IHNT2(1)
519         DO 70 JT=1,IHNT2(3)
520            SCIP(JP,JT)=-1.0
521            B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
522            R2=B2*HIPR1(40)/HIPR1(31)/0.1
523 C                ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
524            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
525      &          /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
526            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
527      &          /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
528            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
529      &           *SQRT(1.0-RRB1)
530            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
531      &           *SQRT(1.0-RRB2)
532            HINT1(18)=HINT1(14)-APHX1*HINT1(15)
533      &                        -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
534            IF(IHPR2(14).EQ.0.OR.
535      &          (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
536               GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
537               RANTOT=RANART(NSEED)
538               IF(RANTOT.GT.GS) GO TO 70
539               GO TO 65
540            ENDIF
541            GSTOT0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
542      &             /HIPR1(31)/2.0*ROMG(0.0)))
543            R2=R2/GSTOT0
544            GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
545            GSTOT=2.0*(1.0-SQRT(1.0-GS))
546            RANTOT=RANART(NSEED)*GSTOT0
547            IF(RANTOT.GT.GSTOT) GO TO 70
548            IF(RANTOT.GT.GS) THEN
549               CALL HIJCSC(JP,JT)
550               GO TO 70
551 C                        ********perform elastic collisions
552            ENDIF
553  65           SCIP(JP,JT)=R2
554            RNIP(JP,JT)=RANTOT
555            SJIP(JP,JT)=HINT1(18)
556            NCOLT=NCOLT+1
557            IPCOL(NCOLT)=JP
558            ITCOL(NCOLT)=JT
559 70        CONTINUE
560 C                ********total number interactions proj and targ has
561 C                                suffered
562
563 clin-5/22/01 write impact parameter:
564         bimp=bb
565 c        write(6,*) '#impact parameter,nlop,ncolt=',bimp,nlop,ncolt
566
567         IF(NCOLT.EQ.0) THEN
568            NLOP=NLOP+1
569            IF(NLOP.LE.20.OR.
570      &           (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
571            RETURN
572         ENDIF
573 C               ********At large impact parameter, there maybe no
574 C                       interaction at all. For NN collision
575 C                       repeat the event until interaction happens
576 C
577         IF(IHPR2(3).NE.0) THEN
578            NHARD=1+INT(RANART(NSEED)*(NCOLT-1)+0.5)
579            NHARD=MIN(NHARD,NCOLT)
580            JPHARD=IPCOL(NHARD)
581            JTHARD=ITCOL(NHARD)
582 clin-6/2009 ctest off:
583 c           write(99,*) IAEVT,NHARD,NCOLT,JPHARD,JTHARD
584         ENDIF
585 C
586         IF(IHPR2(9).EQ.1) THEN
587                 NMINI=1+INT(RANART(NSEED)*(NCOLT-1)+0.5)
588                 NMINI=MIN(NMINI,NCOLT)
589                 JPMINI=IPCOL(NMINI)
590                 JTMINI=ITCOL(NMINI)
591         ENDIF
592 C                ********Specifying the location of the hard and
593 C                        minijet if they are enforced by user
594 C
595         DO 200 JP=1,IHNT2(1)
596         DO 200 JT=1,IHNT2(3)
597         IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
598                 NFP(JP,11)=NFP(JP,11)+1
599                 NFT(JT,11)=NFT(JT,11)+1
600         IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
601                 NP=NP+1
602                 N01=N01+1
603         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
604                 NT=NT+1
605                 N10=N10+1
606         ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
607                 NP=NP+1
608                 NT=NT+1
609                 N0=N0+1
610         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
611                 N11=N11+1
612         ENDIF
613         JOUT=0
614         NFP(JP,10)=0
615         NFT(JT,10)=0
616 C*****************************************************************
617         IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
618 C                ********When IHPR2(8)=0 no jets are produced
619         IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
620 C                ********jets can not be produced for (JP,JT)
621 C                        because not enough energy avaible for 
622 C                                JP or JT 
623         R2=SCIP(JP,JT)
624         HINT1(18)=SJIP(JP,JT)
625         TT=ROMG(R2)*HINT1(18)/HIPR1(31)
626         TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
627         NJET=0
628
629         IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
630            CALL JETINI(JP,JT,1)
631            CALL HIJHRD(JP,JT,0,JFLG,0)
632            HINT1(26)=HINT1(47)
633            HINT1(27)=HINT1(48)
634            HINT1(28)=HINT1(49)
635            HINT1(29)=HINT1(50)
636            HINT1(36)=HINT1(67)
637            HINT1(37)=HINT1(68)
638            HINT1(38)=HINT1(69)
639            HINT1(39)=HINT1(70)
640 C
641            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
642            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
643            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
644      &                                JFLG.GE.3) IASG(NSG,3)=1
645            IHNT2(9)=IHNT2(14)
646            IHNT2(10)=IHNT2(15)
647            DO 105 I05=1,5
648               HINT1(20+I05)=HINT1(40+I05)
649               HINT1(30+I05)=HINT1(50+I05)
650  105           CONTINUE
651 clin-6/2009 ctest off:
652 c           write(99,*) jp,jt,IHPR2(3),HIPR1(10),njet,
653 c     1          ihnt2(9),hint1(21),hint1(22),hint1(23),
654 c     2          ihnt2(10),hint1(31),hint1(32),hint1(33)
655 c           write(99,*) ' '
656            JOUT=1
657            IF(IHPR2(8).EQ.0) GO TO 160
658            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
659      &                /REAL(IHNT2(1))**0.6666667,1.0)
660            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
661      &                /REAL(IHNT2(3))**0.6666667,1.0)
662            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
663      &           *SQRT(1.0-RRB1)
664            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
665      &           *SQRT(1.0-RRB2)
666            HINT1(65)=HINT1(61)-APHX1*HINT1(62)
667      &                        -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
668            TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
669            NJET=-1
670 C                ********subtract the trigger jet from total number
671 C                        of jet production  to be done since it has
672 C                                already been produced here
673            XR1=-ALOG(EXP(-TTRIG)+RANART(NSEED)*(1.0-EXP(-TTRIG)))
674  106           NJET=NJET+1
675            XR1=XR1-ALOG(max(RANART(NSEED),1.0e-20))
676            IF(XR1.LT.TTRIG) GO TO 106
677            XR=0.0
678  107           NJET=NJET+1
679            XR=XR-ALOG(max(RANART(NSEED),1.0e-20))
680            IF(XR.LT.TT-TTRIG) GO TO 107
681            NJET=NJET-1
682            GO TO 112
683         ENDIF
684 C                ********create a hard interaction with specified P_T
685 c                                 when IHPR2(3)>0
686         IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
687 C                ********create at least one pair of mini jets 
688 C                        when IHPR2(9)=1
689 C
690         IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
691      &                (1.0-EXP(-TTS))) GO TO 160
692 C                ********this is the probability for no jet production
693 110        XR=-ALOG(EXP(-TT)+RANART(NSEED)*(1.0-EXP(-TT)))
694 111        NJET=NJET+1
695         XR=XR-ALOG(max(RANART(NSEED),1.0e-20))
696         IF(XR.LT.TT) GO TO 111
697 112        NJET=MIN(NJET,IHPR2(8))
698         IF(IHPR2(8).LT.0)  NJET=ABS(IHPR2(8))
699 C                ******** Determine number of mini jet production
700 C
701         DO 150 ijet=1,NJET
702            CALL JETINI(JP,JT,0)
703            CALL HIJHRD(JP,JT,JOUT,JFLG,1)
704 C                ********JFLG=1 jets valence quarks, JFLG=2 with 
705 C                        gluon jet, JFLG=3 with q-qbar prod for
706 C                        (JP,JT). If JFLG=0 jets can not be produced 
707 C                        this time. If JFLG=-1, error occured abandon
708 C                        this event. JOUT is the total hard scat for
709 C                        (JP,JT) up to now.
710            IF(JFLG.EQ.0) GO TO 160
711            IF(JFLG.LT.0) THEN
712               IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
713               GO TO 50
714            ENDIF
715            JOUT=JOUT+1
716            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
717            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
718            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
719      &                        JFLG.GE.3) IASG(NSG,3)=1
720 C                ******** jet with PT>HIPR1(11) will be quenched
721  150        CONTINUE
722  160        CONTINUE
723
724         CALL HIJSFT(JP,JT,JOUT,IERROR)
725         IF(IERROR.NE.0) THEN
726            IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
727            GO TO 50
728         ENDIF
729 C
730 C                ********conduct soft scattering between JP and JT
731         JATT=JATT+JOUT
732 200        CONTINUE
733 c
734 c**************************
735 c
736 clin-6/2009 write out initial minijet information:
737            call minijet_out(BB)
738            if(pttrig.gt.0.and.ntrig.eq.0) goto 50
739 clin-6/2009 write out initial transverse positions of initial nucleons:
740 c           write(94,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
741         DO 201 JP=1,IHNT2(1)
742 clin-6/2009:
743 c           write(94,203) YP(1,JP)+0.5*BB, YP(2,JP), JP, NFP(JP,5)
744            IF(NFP(JP,5).GT.2) THEN
745               NINP=NINP+1
746            ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
747               NELP=NELP+1
748            ENDIF
749  201    continue
750         DO 202 JT=1,IHNT2(3)
751 clin-6/2009 target nucleon # has a minus sign for distinction from projectile:
752 c           write(94,203) YT(1,JT)-0.5*BB, YT(2,JT), -JT, NFT(JT,5)
753            IF(NFT(JT,5).GT.2) THEN
754               NINTHJ=NINTHJ+1
755            ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
756               NELT=NELT+1
757            ENDIF
758  202    continue
759  203    format(f10.3,1x,f10.3,2(1x,I5))
760 c     
761 c*******************************
762
763
764 C********perform jet quenching for jets with PT>HIPR1(11)**********
765
766         IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
767      &                        IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
768                 DO 271 I=1,IHNT2(1)
769                         IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
770 271                CONTINUE
771                 DO 272 I=1,IHNT2(3)
772                         IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
773 272                CONTINUE
774                 DO 273 ISG=1,NSG
775                         IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
776 273                CONTINUE
777         ENDIF
778
779 clin*****4/09/01-soft1, default way of treating strings:
780         if(isoft.eq.1) then
781 clin-4/16/01 allow fragmentation:
782            isflag=1
783
784 cbz1/25/99
785 c.....transfer data from HIJING to ZPC
786         NSP = IHNT2(1)
787         NST = IHNT2(3)
788         NSI = NSG
789         ISTR = 0
790         NPAR = 0
791         DO 1008 I = 1, IHNT2(1)
792            ISTR = ISTR + 1
793            DO 1007 J = 1, NPJ(I)
794 cbz1/27/99
795 c.....for now only consider gluon cascade
796               IF (KFPJ(I, J) .EQ. 21) THEN
797 cbz1/27/99end
798
799               NPAR = NPAR + 1
800               LSTRG0(NPAR) = ISTR
801               LPART0(NPAR) = J
802               ITYP0(NPAR) = KFPJ(I, J)
803 cbz6/28/99 flow1
804 clin-7/20/01 add dble or sngl to make precisions consistent
805 c              GX0(NPAR) = YP(1, I)
806               GX0(NPAR) = dble(YP(1, I) + 0.5 * BB)
807 cbz6/28/99 flow1 end
808               GY0(NPAR) = dble(YP(2, I))
809               GZ0(NPAR) = 0d0
810               FT0(NPAR) = 0d0
811               PX0(NPAR) = dble(PJPX(I, J))
812               PY0(NPAR) = dble(PJPY(I, J))
813               PZ0(NPAR) = dble(PJPZ(I, J))
814               XMASS0(NPAR) = dble(PJPM(I, J))
815 c              E0(NPAR) = dble(PJPE(I, J))
816               E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
817      1             +PZ0(NPAR)**2+XMASS0(NPAR)**2)
818 clin-7/20/01-end
819
820 cbz1/27/99
821 c.....end gluon selection
822               END IF
823 cbz1/27/99end
824  1007      CONTINUE
825  1008   CONTINUE
826         DO 1010 I = 1, IHNT2(3)
827            ISTR = ISTR + 1
828            DO 1009 J = 1, NTJ(I)
829 cbz1/27/99
830 c.....for now only consider gluon cascade
831               IF (KFTJ(I, J) .EQ. 21) THEN
832 cbz1/27/99end
833               NPAR = NPAR + 1
834               LSTRG0(NPAR) = ISTR
835               LPART0(NPAR) = J
836               ITYP0(NPAR) = KFTJ(I, J)
837 cbz6/28/99 flow1
838 clin-7/20/01 add dble or sngl to make precisions consistent
839 c              GX0(NPAR) = YT(1, I)
840               GX0(NPAR) = dble(YT(1, I) - 0.5 * BB)
841 cbz6/28/99 flow1 end
842               GY0(NPAR) = dble(YT(2, I))
843               GZ0(NPAR) = 0d0
844               FT0(NPAR) = 0d0
845               PX0(NPAR) = dble(PJTX(I, J))
846               PY0(NPAR) = dble(PJTY(I, J))
847               PZ0(NPAR) = dble(PJTZ(I, J))
848               XMASS0(NPAR) = dble(PJTM(I, J))
849 c              E0(NPAR) = dble(PJTE(I, J))
850               E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
851      1             +PZ0(NPAR)**2+XMASS0(NPAR)**2)
852
853 cbz1/27/99
854 c.....end gluon selection
855               END IF
856 cbz1/27/99end
857  1009      CONTINUE
858  1010   CONTINUE
859         DO 1012 I = 1, NSG
860            ISTR = ISTR + 1
861            DO 1011 J = 1, NJSG(I)
862 cbz1/27/99
863 c.....for now only consider gluon cascade
864               IF (K2SG(I, J) .EQ. 21) THEN
865 cbz1/27/99end
866               NPAR = NPAR + 1
867               LSTRG0(NPAR) = ISTR
868               LPART0(NPAR) = J
869               ITYP0(NPAR) = K2SG(I, J)
870 clin-7/20/01 add dble or sngl to make precisions consistent:
871               GX0(NPAR) = 0.5d0 * 
872      1             dble(YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))
873               GY0(NPAR) = 0.5d0 * 
874      2             dble(YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))
875               GZ0(NPAR) = 0d0
876               FT0(NPAR) = 0d0
877               PX0(NPAR) = dble(PXSG(I, J))
878               PY0(NPAR) = dble(PYSG(I, J))
879               PZ0(NPAR) = dble(PZSG(I, J))
880               XMASS0(NPAR) = dble(PMSG(I, J))
881 c              E0(NPAR) = dble(PESG(I, J))
882               E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
883      1             +PZ0(NPAR)**2+XMASS0(NPAR)**2)
884 cbz1/27/99
885 c.....end gluon selection
886               END IF
887 cbz1/27/99end
888  1011      CONTINUE
889  1012   CONTINUE
890         MUL = NPAR
891
892 cbz2/4/99
893         CALL HJANA1
894 cbz2/4/99end
895
896 clin-6/2009:
897         if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
898 c.....call ZPC for parton cascade
899         CALL ZPCMN
900
901 c     write out parton and wounded nucleon information to ana/zpc1.mom:
902 clin-6/2009:
903 c        WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
904 c        WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
905         DO 1013 I = 1, MUL
906 cc           WRITE (14, 411) PX5(I), PY5(I), PZ5(I), ITYP5(I),
907 c     &        XMASS5(I), E5(I)
908            if(dmax1(abs(GX5(I)),abs(GY5(I)),abs(GZ5(I)),abs(FT5(I)))
909      1          .lt.9999) then
910 c              write(14,210) ITYP5(I), PX5(I), PY5(I), PZ5(I), XMASS5(I),
911 c     1             GX5(I), GY5(I), GZ5(I), FT5(I)
912            else
913 c     change format for large numbers:
914 c              write(14,211) ITYP5(I), PX5(I), PY5(I), PZ5(I), XMASS5(I),
915 c     1             GX5(I), GY5(I), GZ5(I), FT5(I)
916            endif
917
918  1013   CONTINUE
919  210    format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
920  211    format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
921  395    format(3I8,f10.4,4I5)
922
923 clin-4/09/01:
924         itest=itest+1
925 c 411    FORMAT(1X, 3F10.3, I6, 2F10.3)
926 cbz3/19/99 end
927
928 clin-5/2009 ctest off:
929 c        call frztm(1,1)
930
931 c.....transfer data back from ZPC to HIJING
932         DO 1014 I = 1, MUL
933            IF (LSTRG1(I) .LE. NSP) THEN
934               NSTRG = LSTRG1(I)
935               NPART = LPART1(I)
936               KFPJ(NSTRG, NPART) = ITYP5(I)
937 clin-7/20/01 add dble or sngl to make precisions consistent
938               PJPX(NSTRG, NPART) = sngl(PX5(I))
939               PJPY(NSTRG, NPART) = sngl(PY5(I))
940               PJPZ(NSTRG, NPART) = sngl(PZ5(I))
941               PJPE(NSTRG, NPART) = sngl(E5(I))
942               PJPM(NSTRG, NPART) = sngl(XMASS5(I))
943            ELSE IF (LSTRG1(I) .LE. NSP + NST) THEN
944               NSTRG = LSTRG1(I) - NSP
945               NPART = LPART1(I)
946               KFTJ(NSTRG, NPART) = ITYP5(I)
947               PJTX(NSTRG, NPART) = sngl(PX5(I))
948               PJTY(NSTRG, NPART) = sngl(PY5(I))
949               PJTZ(NSTRG, NPART) = sngl(PZ5(I))
950               PJTE(NSTRG, NPART) = sngl(E5(I))
951               PJTM(NSTRG, NPART) = sngl(XMASS5(I))
952            ELSE
953               NSTRG = LSTRG1(I) - NSP - NST
954               NPART = LPART1(I)
955               K2SG(NSTRG, NPART) = ITYP5(I)
956               PXSG(NSTRG, NPART) = sngl(PX5(I))
957               PYSG(NSTRG, NPART) = sngl(PY5(I))
958               PZSG(NSTRG, NPART) = sngl(PZ5(I))
959               PESG(NSTRG, NPART) = sngl(E5(I))
960               PMSG(NSTRG, NPART) = sngl(XMASS5(I))
961            END IF
962  1014   CONTINUE
963 cbz1/25/99end
964
965 cbz2/4/99
966         CALL HJANA2
967 cbz2/4/99end
968
969 clin*****4/09/01-soft2, put q+dq+X in strings into ZPC:
970         elseif(isoft.eq.2) then
971         NSP = IHNT2(1)
972         NST = IHNT2(3)
973 clin-4/27/01:
974         NSI = NSG
975         NPAR=0
976         ISTR=0
977 C
978 clin  No fragmentation to hadrons, only on parton level, 
979 c     and transfer minijet and string data from HIJING to ZPC:
980         MSTJ(1)=0
981 clin-4/12/01 forbid soft radiation before ZPC to avoid small-mass strings,
982 c     and forbid jet order reversal before ZPC to avoid unphysical flavors:
983         IHPR2(1)=0
984         isflag=0
985
986         IF(IHPR2(20).NE.0) THEN
987            DO 320 NTP=1,2
988               DO 310 jjtp=1,IHNT2(2*NTP-1)
989                  ISTR = ISTR + 1
990 c change: do gluon kink only once: either here or in fragmentation.
991                  CALL HIJFRG(jjtp,NTP,IERROR)
992 c                 call lulist(1)
993                  if(NTP.eq.1) then
994 c 354                continue
995                     NPJ(jjtp)=MAX0(N-2,0)
996
997 clin-4/12/01:                    NPJ(jjtp)=MAX0(ipartn-2,0)
998                  else
999 c 355                continue
1000                     NTJ(jjtp)=MAX0(N-2,0)
1001 clin-4/12/01:                    NTJ(jjtp)=MAX0(ipartn-2,0)
1002                  endif
1003
1004                  do 300 ii=1,N
1005                  NPAR = NPAR + 1
1006                  LSTRG0(NPAR) = ISTR
1007                  LPART0(NPAR) = II
1008                  ITYP0(NPAR) = K(II,2)
1009                  GZ0(NPAR) = 0d0
1010                  FT0(NPAR) = 0d0
1011 clin-7/20/01 add dble or sngl to make precisions consistent
1012                  PX0(NPAR) = dble(P(II,1))
1013                  PY0(NPAR) = dble(P(II,2))
1014                  PZ0(NPAR) = dble(P(II,3))
1015                  XMASS0(NPAR) = dble(P(II,5))
1016 c                 E0(NPAR) = dble(P(II,4))
1017                  E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
1018      1                +PZ0(NPAR)**2+XMASS0(NPAR)**2)
1019                  IF (NTP .EQ. 1) THEN
1020 clin-7/20/01 add dble or sngl to make precisions consistent
1021                     GX0(NPAR) = dble(YP(1, jjtp)+0.5 * BB)
1022                     GY0(NPAR) = dble(YP(2, jjtp))
1023                     IITYP=ITYP0(NPAR)
1024                     nstrg=LSTRG0(NPAR)
1025                     if(IITYP.eq.2112.or.IITYP.eq.2212) then
1026                     elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1027      1 (II.eq.1.or.II.eq.N)) then
1028                        PP(nstrg,6)=sngl(PX0(NPAR))
1029                        PP(nstrg,7)=sngl(PY0(NPAR))
1030                        PP(nstrg,14)=sngl(XMASS0(NPAR))
1031                     elseif((IITYP.eq.1103.or.IITYP.eq.2101
1032      1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1033      2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1034      3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1035      4 .and.(II.eq.1.or.II.eq.N)) then
1036                        PP(nstrg,8)=sngl(PX0(NPAR))
1037                        PP(nstrg,9)=sngl(PY0(NPAR))
1038                        PP(nstrg,15)=sngl(XMASS0(NPAR))
1039                     else
1040                        NPART = LPART0(NPAR)-1
1041                        KFPJ(NSTRG, NPART) = ITYP0(NPAR)
1042                        PJPX(NSTRG, NPART) = sngl(PX0(NPAR))
1043                        PJPY(NSTRG, NPART) = sngl(PY0(NPAR))
1044                        PJPZ(NSTRG, NPART) = sngl(PZ0(NPAR))
1045                        PJPE(NSTRG, NPART) = sngl(E0(NPAR))
1046                        PJPM(NSTRG, NPART) = sngl(XMASS0(NPAR))
1047                     endif
1048                  ELSE
1049                     GX0(NPAR) = dble(YT(1, jjtp)-0.5 * BB)
1050                     GY0(NPAR) = dble(YT(2, jjtp))
1051                     IITYP=ITYP0(NPAR)
1052                     nstrg=LSTRG0(NPAR)-NSP
1053                     if(IITYP.eq.2112.or.IITYP.eq.2212) then
1054                     elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1055      1 (II.eq.1.or.II.eq.N)) then
1056                        PT(nstrg,6)=sngl(PX0(NPAR))
1057                        PT(nstrg,7)=sngl(PY0(NPAR))
1058                        PT(nstrg,14)=sngl(XMASS0(NPAR))
1059                     elseif((IITYP.eq.1103.or.IITYP.eq.2101
1060      1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1061      2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1062      3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1063      4 .and.(II.eq.1.or.II.eq.N)) then
1064                        PT(nstrg,8)=sngl(PX0(NPAR))
1065                        PT(nstrg,9)=sngl(PY0(NPAR))
1066                        PT(nstrg,15)=sngl(XMASS0(NPAR))
1067                     else
1068                        NPART = LPART0(NPAR)-1
1069                        KFTJ(NSTRG, NPART) = ITYP0(NPAR)
1070                        PJTX(NSTRG, NPART) = sngl(PX0(NPAR))
1071                        PJTY(NSTRG, NPART) = sngl(PY0(NPAR))
1072                        PJTZ(NSTRG, NPART) = sngl(PZ0(NPAR))
1073                        PJTE(NSTRG, NPART) = sngl(E0(NPAR))
1074                        PJTM(NSTRG, NPART) = sngl(XMASS0(NPAR))
1075                     endif
1076                  END IF
1077  300          continue
1078  310          continue
1079  320       continue
1080            DO 330 ISG=1,NSG
1081               ISTR = ISTR + 1
1082               CALL HIJFRG(ISG,3,IERROR)
1083 c              call lulist(2)
1084 c
1085               NJSG(ISG)=N
1086 c
1087               do 1001 ii=1,N
1088                  NPAR = NPAR + 1
1089                  LSTRG0(NPAR) = ISTR
1090                  LPART0(NPAR) = II
1091                  ITYP0(NPAR) = K(II,2)
1092                  GX0(NPAR)=0.5d0*
1093      1                dble(YP(1,IASG(ISG,1))+YT(1,IASG(ISG,2)))
1094                  GY0(NPAR)=0.5d0*
1095      2                dble(YP(2,IASG(ISG,1))+YT(2,IASG(ISG,2)))
1096                  GZ0(NPAR) = 0d0
1097                  FT0(NPAR) = 0d0
1098                  PX0(NPAR) = dble(P(II,1))
1099                  PY0(NPAR) = dble(P(II,2))
1100                  PZ0(NPAR) = dble(P(II,3))
1101                  XMASS0(NPAR) = dble(P(II,5))
1102 c                 E0(NPAR) = dble(P(II,4))
1103                  E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
1104      1                +PZ0(NPAR)**2+XMASS0(NPAR)**2)
1105  1001         continue
1106  330       continue
1107         endif
1108
1109         MUL = NPAR
1110 cbz2/4/99
1111         CALL HJANA1
1112 cbz2/4/99end
1113 clin-6/2009:
1114         if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
1115 c.....call ZPC for parton cascade
1116         CALL ZPCMN
1117 cbz3/19/99
1118 clin-6/2009:
1119 c        WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
1120         WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
1121         itest=itest+1
1122
1123         DO 1015 I = 1, MUL
1124 c           WRITE (14, 311) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1125 c     &        XMASS5(I), E5(I)
1126            WRITE (14, 312) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1127      &        XMASS5(I), E5(I),LSTRG1(I), LPART1(I)
1128  1015   CONTINUE
1129 c 311    FORMAT(1X, 3F10.4, I6, 2F10.4)
1130  312    FORMAT(1X, 3F10.3, I6, 2F10.3,1X,I6,1X,I3)
1131 cbz3/19/99 end
1132
1133 clin-5/2009 ctest off:
1134 c        call frztm(1,1)
1135
1136 clin-4/13/01 initialize four momenta and invariant mass of strings after ZPC:
1137         do 1004 nmom=1,5
1138            do 1002 nstrg=1,nsp
1139               PP(nstrg,nmom)=0.
1140  1002      continue
1141            do 1003 nstrg=1,nst
1142               PT(nstrg,nmom)=0.
1143  1003      continue
1144  1004   continue
1145 clin-4/13/01-end
1146
1147         DO 1005 I = 1, MUL
1148            IITYP=ITYP5(I)
1149            IF (LSTRG1(I) .LE. NSP) THEN
1150               NSTRG = LSTRG1(I)
1151 c     nucleons without interactions:
1152               if(IITYP.eq.2112.or.IITYP.eq.2212) then
1153 clin-7/20/01 add dble or sngl to make precisions consistent
1154                  PP(nstrg,1)=sngl(PX5(I))
1155                  PP(nstrg,2)=sngl(PY5(I))
1156                  PP(nstrg,3)=sngl(PZ5(I))
1157                  PP(nstrg,4)=sngl(E5(I))
1158                  PP(nstrg,5)=sngl(XMASS5(I))
1159 c     valence quark:
1160               elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1161      1 (LPART1(I).eq.1.or.LPART1(I).eq.(NPJ(NSTRG)+2))) then
1162                  PP(nstrg,6)=sngl(PX5(I))
1163                  PP(nstrg,7)=sngl(PY5(I))
1164                  PP(nstrg,14)=sngl(XMASS5(I))
1165                  PP(nstrg,1)=PP(nstrg,1)+sngl(PX5(I))
1166                  PP(nstrg,2)=PP(nstrg,2)+sngl(PY5(I))
1167                  PP(nstrg,3)=PP(nstrg,3)+sngl(PZ5(I))
1168                  PP(nstrg,4)=PP(nstrg,4)+sngl(E5(I))
1169                  PP(nstrg,5)=sqrt(PP(nstrg,4)**2-PP(nstrg,1)**2
1170      1                -PP(nstrg,2)**2-PP(nstrg,3)**2)
1171 c     diquark:
1172               elseif((IITYP.eq.1103.or.IITYP.eq.2101
1173      1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1174      2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1175      3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1176      4 .and.(LPART1(I).eq.1.or.LPART1(I).eq.(NPJ(NSTRG)+2))) then
1177                  PP(nstrg,8)=sngl(PX5(I))
1178                  PP(nstrg,9)=sngl(PY5(I))
1179                  PP(nstrg,15)=sngl(XMASS5(I))
1180                  PP(nstrg,1)=PP(nstrg,1)+sngl(PX5(I))
1181                  PP(nstrg,2)=PP(nstrg,2)+sngl(PY5(I))
1182                  PP(nstrg,3)=PP(nstrg,3)+sngl(PZ5(I))
1183                  PP(nstrg,4)=PP(nstrg,4)+sngl(E5(I))
1184                  PP(nstrg,5)=sqrt(PP(nstrg,4)**2-PP(nstrg,1)**2
1185      1                -PP(nstrg,2)**2-PP(nstrg,3)**2)
1186 c     partons in projectile or target strings:
1187               else
1188                  NPART = LPART1(I)-1
1189                  KFPJ(NSTRG, NPART) = ITYP5(I)
1190                  PJPX(NSTRG, NPART) = sngl(PX5(I))
1191                  PJPY(NSTRG, NPART) = sngl(PY5(I))
1192                  PJPZ(NSTRG, NPART) = sngl(PZ5(I))
1193                  PJPE(NSTRG, NPART) = sngl(E5(I))
1194                  PJPM(NSTRG, NPART) = sngl(XMASS5(I))
1195               endif
1196            ELSE IF (LSTRG1(I) .LE. NSP + NST) THEN
1197               NSTRG = LSTRG1(I) - NSP
1198               if(IITYP.eq.2112.or.IITYP.eq.2212) then
1199                  PT(nstrg,1)=sngl(PX5(I))
1200                  PT(nstrg,2)=sngl(PY5(I))
1201                  PT(nstrg,3)=sngl(PZ5(I))
1202                  PT(nstrg,4)=sngl(E5(I))
1203                  PT(nstrg,5)=sngl(XMASS5(I))
1204               elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1205      1 (LPART1(I).eq.1.or.LPART1(I).eq.(NTJ(NSTRG)+2))) then
1206                  PT(nstrg,6)=sngl(PX5(I))
1207                  PT(nstrg,7)=sngl(PY5(I))
1208                  PT(nstrg,14)=sngl(XMASS5(I))
1209                  PT(nstrg,1)=PT(nstrg,1)+sngl(PX5(I))
1210                  PT(nstrg,2)=PT(nstrg,2)+sngl(PY5(I))
1211                  PT(nstrg,3)=PT(nstrg,3)+sngl(PZ5(I))
1212                  PT(nstrg,4)=PT(nstrg,4)+sngl(E5(I))
1213                  PT(nstrg,5)=sqrt(PT(nstrg,4)**2-PT(nstrg,1)**2
1214      1                -PT(nstrg,2)**2-PT(nstrg,3)**2)
1215               elseif((IITYP.eq.1103.or.IITYP.eq.2101
1216      1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1217      2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1218      3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1219      4 .and.(LPART1(I).eq.1.or.LPART1(I).eq.(NTJ(NSTRG)+2))) then
1220                  PT(nstrg,8)=sngl(PX5(I))
1221                  PT(nstrg,9)=sngl(PY5(I))
1222                  PT(nstrg,15)=sngl(XMASS5(I))
1223                  PT(nstrg,1)=PT(nstrg,1)+sngl(PX5(I))
1224                  PT(nstrg,2)=PT(nstrg,2)+sngl(PY5(I))
1225                  PT(nstrg,3)=PT(nstrg,3)+sngl(PZ5(I))
1226                  PT(nstrg,4)=PT(nstrg,4)+sngl(E5(I))
1227                  PT(nstrg,5)=sqrt(PT(nstrg,4)**2-PT(nstrg,1)**2
1228      1                -PT(nstrg,2)**2-PT(nstrg,3)**2)
1229               else
1230                  NPART = LPART1(I)-1
1231                  KFTJ(NSTRG, NPART) = ITYP5(I)
1232                  PJTX(NSTRG, NPART) = sngl(PX5(I))
1233                  PJTY(NSTRG, NPART) = sngl(PY5(I))
1234                  PJTZ(NSTRG, NPART) = sngl(PZ5(I))
1235                  PJTE(NSTRG, NPART) = sngl(E5(I))
1236                  PJTM(NSTRG, NPART) = sngl(XMASS5(I))
1237               endif
1238            ELSE
1239               NSTRG = LSTRG1(I) - NSP - NST
1240               NPART = LPART1(I)
1241               K2SG(NSTRG, NPART) = ITYP5(I)
1242               PXSG(NSTRG, NPART) = sngl(PX5(I))
1243               PYSG(NSTRG, NPART) = sngl(PY5(I))
1244               PZSG(NSTRG, NPART) = sngl(PZ5(I))
1245               PESG(NSTRG, NPART) = sngl(E5(I))
1246               PMSG(NSTRG, NPART) = sngl(XMASS5(I))
1247            END IF
1248  1005   CONTINUE
1249 cbz1/25/99end
1250
1251 clin-4/09/01  turn on fragmentation with soft radiation 
1252 c     and jet order reversal to form hadrons after ZPC:
1253         MSTJ(1)=1
1254         IHPR2(1)=1
1255         isflag=1
1256 clin-4/13/01 allow small mass strings (D=1.5GeV):
1257         HIPR1(1)=0.94
1258
1259 cbz2/4/99
1260         CALL HJANA2
1261 cbz2/4/99end
1262
1263 clin-4/19/01-soft3, fragment strings, then convert hadrons to partons 
1264 c     and input to ZPC:
1265         elseif(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
1266 clin-4/24/01 normal fragmentation first:
1267         isflag=0
1268
1269         IF(IHPR2(20).NE.0) THEN
1270            DO 560 ISG=1,NSG
1271                 CALL HIJFRG(ISG,3,IERROR)
1272 C
1273                 nsbst=1
1274                 IDSTR=92
1275                 IF(IHPR2(21).EQ.0) THEN
1276                    CALL LUEDIT(2)
1277                 ELSE
1278  551                   nsbst=nsbst+1
1279                    IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO  551
1280                    IDSTR=K(nsbst,2)
1281                    nsbst=nsbst+1
1282                 ENDIF
1283
1284                 IF(FRAME.EQ.'LAB') THEN
1285                         CALL HBOOST
1286                 ENDIF
1287 C                ******** boost back to lab frame(if it was in)
1288 C
1289                 nsbstR=0
1290                 DO 560 I=nsbst,N
1291                    IF(K(I,2).EQ.IDSTR) THEN
1292                       nsbstR=nsbstR+1
1293                       GO TO 560
1294                    ENDIF
1295                    K(I,4)=nsbstR
1296                    NATT=NATT+1
1297                    KATT(NATT,1)=K(I,2)
1298                    KATT(NATT,2)=20
1299                    KATT(NATT,4)=K(I,1)
1300 c     from Yasushi, to avoid violation of array limits:
1301 c                   IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1302 clin-4/2008 to avoid out-of-bound error in K():
1303 c                   IF(K(I,3).EQ.0 .OR. 
1304 c     1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1305 c                      KATT(NATT,3)=0
1306                    IF(K(I,3).EQ.0) THEN
1307                       KATT(NATT,3)=0
1308                    ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1309                       KATT(NATT,3)=0
1310 clin-4/2008-end
1311                    ELSE
1312                       KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1313                    ENDIF
1314
1315 C       ****** identify the mother particle
1316                    PATT(NATT,1)=P(I,1)
1317                    PATT(NATT,2)=P(I,2)
1318                    PATT(NATT,3)=P(I,3)
1319                    PATT(NATT,4)=P(I,4)
1320                    EATT=EATT+P(I,4)
1321                    GXAR(NATT) = 0.5 * (YP(1, IASG(ISG, 1)) +
1322      &                YT(1, IASG(ISG, 2)))
1323                    GYAR(NATT) = 0.5 * (YP(2, IASG(ISG, 1)) +
1324      &                YT(2, IASG(ISG, 2)))
1325                    GZAR(NATT) = 0.
1326                    FTAR(NATT) = 0.
1327                    ITYPAR(NATT) = K(I, 2)
1328                    PXAR(NATT) = P(I, 1)
1329                    PYAR(NATT) = P(I, 2)
1330                    PZAR(NATT) = P(I, 3)
1331                    PEAR(NATT) = P(I, 4)
1332                    XMAR(NATT) = P(I, 5)
1333 cbz11/11/98end
1334
1335  560            CONTINUE
1336 C                ********Fragment the q-qbar jets systems *****
1337 C
1338            JTP(1)=IHNT2(1)
1339            JTP(2)=IHNT2(3)
1340            DO 600 NTP=1,2
1341            DO 600 jjtp=1,JTP(NTP)
1342                 CALL HIJFRG(jjtp,NTP,IERROR)
1343 C
1344                 nsbst=1
1345                 IDSTR=92
1346                 IF(IHPR2(21).EQ.0) THEN
1347                    CALL LUEDIT(2)
1348                 ELSE
1349  581                   nsbst=nsbst+1
1350                    IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO  581
1351                    IDSTR=K(nsbst,2)
1352                    nsbst=nsbst+1
1353                 ENDIF
1354                 IF(FRAME.EQ.'LAB') THEN
1355                         CALL HBOOST
1356                 ENDIF
1357 C                ******** boost back to lab frame(if it was in)
1358 C
1359                 NFTP=NFP(jjtp,5)
1360                 IF(NTP.EQ.2) NFTP=10+NFT(jjtp,5)
1361                 nsbstR=0
1362                 DO 590 I=nsbst,N
1363                    IF(K(I,2).EQ.IDSTR) THEN
1364                       nsbstR=nsbstR+1
1365                       GO TO 590
1366                    ENDIF
1367                    K(I,4)=nsbstR
1368                    NATT=NATT+1
1369                    KATT(NATT,1)=K(I,2)
1370                    KATT(NATT,2)=NFTP
1371                    KATT(NATT,4)=K(I,1)
1372 c                   IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1373 clin-4/2008
1374 c                   IF(K(I,3).EQ.0 .OR.
1375 c     1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1376 c                      KATT(NATT,3)=0
1377                    IF(K(I,3).EQ.0) THEN
1378                       KATT(NATT,3)=0
1379                    ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1380                       KATT(NATT,3)=0
1381 clin-4/2008-end
1382                    ELSE
1383                       KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1384                    ENDIF
1385
1386 C       ****** identify the mother particle
1387                    PATT(NATT,1)=P(I,1)
1388                    PATT(NATT,2)=P(I,2)
1389                    PATT(NATT,3)=P(I,3)
1390                    PATT(NATT,4)=P(I,4)
1391                    EATT=EATT+P(I,4)
1392                    IF (NTP .EQ. 1) THEN
1393                       GXAR(NATT) = YP(1, jjtp)+0.5 * BB
1394                       GYAR(NATT) = YP(2, jjtp)
1395                    ELSE
1396                       GXAR(NATT) = YT(1, jjtp)-0.5 * BB
1397                       GYAR(NATT) = YT(2, jjtp)
1398                    END IF
1399                    GZAR(NATT) = 0.
1400                    FTAR(NATT) = 0.
1401                    ITYPAR(NATT) = K(I, 2)
1402                    PXAR(NATT) = P(I, 1)
1403                    PYAR(NATT) = P(I, 2)
1404                    PZAR(NATT) = P(I, 3)
1405                    PEAR(NATT) = P(I, 4)
1406                    XMAR(NATT) = P(I, 5)
1407 cbz11/11/98end
1408
1409  590                CONTINUE 
1410  600           CONTINUE
1411 C     ********Fragment the q-qq related string systems
1412         ENDIF
1413 clin-4/2008 check for zero NDR value:
1414         if(NDR.ge.1) then
1415 c
1416         DO 650 I=1,NDR
1417                 NATT=NATT+1
1418                 KATT(NATT,1)=KFDR(I)
1419                 KATT(NATT,2)=40
1420                 KATT(NATT,3)=0
1421                 PATT(NATT,1)=PDR(I,1)
1422                 PATT(NATT,2)=PDR(I,2)
1423                 PATT(NATT,3)=PDR(I,3)
1424                 PATT(NATT,4)=PDR(I,4)
1425                 EATT=EATT+PDR(I,4)
1426 clin-11/11/03     set direct photons positions and time at formation:
1427                 GXAR(NATT) = rtdr(I,1)
1428                 GYAR(NATT) = rtdr(I,2)
1429                 GZAR(NATT) = 0.
1430                 FTAR(NATT) = 0.
1431                 ITYPAR(NATT) =KATT(NATT,1) 
1432                 PXAR(NATT) = PATT(NATT,1)
1433                 PYAR(NATT) = PATT(NATT,2)
1434                 PZAR(NATT) = PATT(NATT,3)
1435                 PEAR(NATT) = PATT(NATT,4)
1436                 XMAR(NATT) = PDR(I,5)
1437  650        CONTINUE
1438 clin-4/2008:
1439          endif
1440 clin-6/2009 ctest on:
1441          call embedHighPt
1442 c
1443         CALL HJANA1
1444
1445 clin-4/19/01 convert hadrons to partons for ZPC (with GX0 given):
1446         call htop
1447
1448 clin-7/03/01 move up, used in zpstrg (otherwise not set and incorrect):
1449         nsp=0
1450         nst=0
1451         nsg=natt
1452         NSI=NSG
1453 clin-7/03/01-end
1454
1455 clin-6/2009:
1456         if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
1457 c.....call ZPC for parton cascade
1458         CALL ZPCMN
1459 clin-6/2009:
1460 c        WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
1461         WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
1462         itest=itest+1
1463
1464         DO 1016 I = 1, MUL
1465 c           WRITE (14, 511) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1466 c     &        XMASS5(I), E5(I)
1467            WRITE (14, 512) ITYP5(I), PX5(I), PY5(I), PZ5(I), 
1468      &        XMASS5(I), LSTRG1(I), LPART1(I), FT5(I)
1469
1470  1016   CONTINUE
1471 c 511    FORMAT(1X, 3F10.4, I6, 2F10.4)
1472  512    FORMAT(I6,4(1X,F10.3),1X,I6,1X,I3,1X,F10.3)
1473 c 513    FORMAT(1X, 4F10.4)
1474
1475 clin-5/2009 ctest off:
1476 c        call frztm(1,1)
1477
1478 clin  save data after ZPC for fragmentation purpose:
1479 c.....transfer data back from ZPC to HIJING
1480         DO 1018 I = 1, MAXSTR
1481            DO 1017 J = 1, 3
1482               K1SGS(I, J) = 0
1483               K2SGS(I, J) = 0
1484               PXSGS(I, J) = 0d0
1485               PYSGS(I, J) = 0d0
1486               PZSGS(I, J) = 0d0
1487               PESGS(I, J) = 0d0
1488               PMSGS(I, J) = 0d0
1489               GXSGS(I, J) = 0d0
1490               GYSGS(I, J) = 0d0
1491               GZSGS(I, J) = 0d0
1492               FTSGS(I, J) = 0d0
1493  1017      CONTINUE
1494  1018   CONTINUE
1495         DO 1019 I = 1, MUL
1496            IITYP=ITYP5(I)
1497            NSTRG = LSTRG1(I)
1498            NPART = LPART1(I)
1499            K2SGS(NSTRG, NPART) = ITYP5(I)
1500            PXSGS(NSTRG, NPART) = PX5(I)
1501            PYSGS(NSTRG, NPART) = PY5(I)
1502            PZSGS(NSTRG, NPART) = PZ5(I)
1503            PMSGS(NSTRG, NPART) = XMASS5(I)
1504 clin-7/20/01 E5(I) does no include the finite parton mass XMASS5(I), 
1505 c     so define it anew:
1506 c           PESGS(NSTRG, NPART) = E5(I)
1507 c           if(abs(PZ5(i)/E5(i)).gt.0.9999999d0) 
1508 c     1          write(91,*) 'a',PX5(i),PY5(i),XMASS5(i),PZ5(i),E5(i)
1509            E5(I)=dsqrt(PX5(I)**2+PY5(I)**2+PZ5(I)**2+XMASS5(I)**2)
1510            PESGS(NSTRG, NPART) = E5(I)
1511 c           if(abs(PZ5(i)/E5(i)).gt.0.9999999d0) 
1512 c     1          write(91,*) 'b: new E5(I)=',E5(i)
1513 clin-7/20/01-end
1514            GXSGS(NSTRG, NPART) = GX5(I)
1515            GYSGS(NSTRG, NPART) = GY5(I)
1516            GZSGS(NSTRG, NPART) = GZ5(I)
1517            FTSGS(NSTRG, NPART) = FT5(I)
1518  1019   CONTINUE
1519         CALL HJANA2
1520
1521 clin-4/19/01-end
1522
1523         endif
1524 clin-4/09/01-end
1525
1526 C
1527 C**************fragment all the string systems in the following*****
1528 C
1529 C********nsbst is where particle information starts
1530 C********nsbstR+1 is the number of strings in fragmentation
1531 C********the number of strings before a line is stored in K(I,4)
1532 C********IDSTR is id number of the string system (91,92 or 93)
1533 C
1534 clin-4/30/01 convert partons to hadrons after ZPC:
1535         if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
1536            NATT=0
1537            EATT=0.
1538            call ptoh
1539            do 1006 I=1,nnozpc
1540               NATT=NATT+1
1541               KATT(NATT,1)=ITYPN(I)
1542               PATT(NATT,1)=PXN(I)
1543               PATT(NATT,2)=PYN(I)
1544               PATT(NATT,3)=PZN(I)
1545               PATT(NATT,4)=EEN(I)
1546               EATT=EATT+EEN(I)
1547               GXAR(NATT)=GXN(I)
1548               GYAR(NATT)=GYN(I)
1549               GZAR(NATT)=GZN(I)
1550               FTAR(NATT)=FTN(I)
1551               ITYPAR(NATT)=ITYPN(I)
1552               PXAR(NATT)=PXN(I)
1553               PYAR(NATT)=PYN(I)
1554               PZAR(NATT)=PZN(I)
1555               PEAR(NATT)=EEN(I)
1556               XMAR(NATT)=XMN(I)
1557  1006      continue
1558            goto 565
1559         endif
1560 clin-4/30/01-end        
1561         IF(IHPR2(20).NE.0) THEN
1562            DO 360 ISG=1,NSG
1563                 CALL HIJFRG(ISG,3,IERROR)
1564                 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
1565                    MSTU(24)=0
1566                    MSTU(28)=0
1567                    IF(IHPR2(10).NE.0) THEN
1568 c                      call lulist(2)
1569                       WRITE(6,*) 'error occured ISG, repeat the event'
1570                   write(6,*) ISG
1571
1572                    ENDIF
1573                    GO TO 50
1574                 ENDIF
1575 C                        ********Check errors
1576 C
1577                 nsbst=1
1578                 IDSTR=92
1579                 IF(IHPR2(21).EQ.0) THEN
1580                    CALL LUEDIT(2)
1581                 ELSE
1582 351                   nsbst=nsbst+1
1583                    IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO  351
1584                    IDSTR=K(nsbst,2)
1585                    nsbst=nsbst+1
1586                 ENDIF
1587 C
1588                 IF(FRAME.EQ.'LAB') THEN
1589                         CALL HBOOST
1590                 ENDIF
1591 C                ******** boost back to lab frame(if it was in)
1592 C
1593                 nsbstR=0
1594                 DO 360 I=nsbst,N
1595                    IF(K(I,2).EQ.IDSTR) THEN
1596                       nsbstR=nsbstR+1
1597                       GO TO 360
1598                    ENDIF
1599                    K(I,4)=nsbstR
1600                    NATT=NATT+1
1601                    KATT(NATT,1)=K(I,2)
1602                    KATT(NATT,2)=20
1603                    KATT(NATT,4)=K(I,1)
1604 c                   IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1605 clin-4/2008:
1606 c                   IF(K(I,3).EQ.0 .OR. 
1607 c     1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1608 c                      KATT(NATT,3)=0
1609                    IF(K(I,3).EQ.0) THEN
1610                       KATT(NATT,3)=0
1611                    ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1612                       KATT(NATT,3)=0
1613 clin-4/2008-end
1614                    ELSE
1615                       KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1616                    ENDIF
1617
1618 C       ****** identify the mother particle
1619                    PATT(NATT,1)=P(I,1)
1620                    PATT(NATT,2)=P(I,2)
1621                    PATT(NATT,3)=P(I,3)
1622                    PATT(NATT,4)=P(I,4)
1623                    EATT=EATT+P(I,4)
1624
1625 cbz11/11/98
1626 cbz1/25/99
1627 c                   GXAR(NATT) = 0.5 * (YP(1, IASG(ISG, 1)) +
1628 c     &                YT(1, IASG(ISG, 2)))
1629 c                   GYAR(NATT) = 0.5 * (YP(2, IASG(ISG, 1)) +
1630 c     &                YT(2, IASG(ISG, 2)))
1631                    LSG = NSP + NST + ISG
1632                    GXAR(NATT) = sngl(ZT1(LSG))
1633                    GYAR(NATT) = sngl(ZT2(LSG))
1634                    GZAR(NATT) = sngl(ZT3(LSG))
1635                    FTAR(NATT) = sngl(ATAUI(LSG))
1636 cbz1/25/99end
1637                    ITYPAR(NATT) = K(I, 2)
1638                    PXAR(NATT) = P(I, 1)
1639                    PYAR(NATT) = P(I, 2)
1640                    PZAR(NATT) = P(I, 3)
1641                    PEAR(NATT) = P(I, 4)
1642                    XMAR(NATT) = P(I, 5)
1643 cbz11/11/98end
1644
1645 360           CONTINUE
1646 C                ********Fragment the q-qbar jets systems *****
1647 C
1648            JTP(1)=IHNT2(1)
1649            JTP(2)=IHNT2(3)
1650            DO 400 NTP=1,2
1651            DO 400 jjtp=1,JTP(NTP)
1652                 CALL HIJFRG(jjtp,NTP,IERROR)
1653                 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
1654                    MSTU(24)=0
1655                    MSTU(28)=0
1656                    IF(IHPR2(10).NE.0) THEN
1657 c                  call lulist(2)
1658                   WRITE(6,*) 'error occured P&T, repeat the event'
1659                   WRITE(6,*) NTP,jjtp
1660 clin-6/2009 when this happens, the event will be repeated, 
1661 c     and another record for the same event number will be written into
1662 c     zpc.dat, zpc.res, minijet-initial-beforePropagation.dat,
1663 c     parton-initial-afterPropagation.dat, parton-after-coalescence.dat, 
1664 c     and parton-collisionsHistory.dat. 
1665                    ENDIF
1666                    GO TO 50
1667                 ENDIF
1668 C                        ********check errors
1669 C
1670                 nsbst=1
1671                 IDSTR=92
1672                 IF(IHPR2(21).EQ.0) THEN
1673                    CALL LUEDIT(2)
1674                 ELSE
1675 381                   nsbst=nsbst+1
1676                    IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO  381
1677                    IDSTR=K(nsbst,2)
1678                    nsbst=nsbst+1
1679                 ENDIF
1680                 IF(FRAME.EQ.'LAB') THEN
1681                         CALL HBOOST
1682                 ENDIF
1683 C                ******** boost back to lab frame(if it was in)
1684 C
1685                 NFTP=NFP(jjtp,5)
1686                 IF(NTP.EQ.2) NFTP=10+NFT(jjtp,5)
1687                 nsbstR=0
1688                 DO 390 I=nsbst,N
1689                    IF(K(I,2).EQ.IDSTR) THEN
1690                       nsbstR=nsbstR+1
1691                       GO TO 390
1692                    ENDIF
1693                    K(I,4)=nsbstR
1694                    NATT=NATT+1
1695                    KATT(NATT,1)=K(I,2)
1696                    KATT(NATT,2)=NFTP
1697                    KATT(NATT,4)=K(I,1)
1698 c                   IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1699 clin-4/2008:
1700 c                   IF(K(I,3).EQ.0 .OR. 
1701 c     1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1702 c                      KATT(NATT,3)=0
1703                    IF(K(I,3).EQ.0) THEN
1704                       KATT(NATT,3)=0
1705                    ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1706                       KATT(NATT,3)=0
1707 clin-4/2008-end
1708                    ELSE
1709                       KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1710                    ENDIF
1711 C       ****** identify the mother particle
1712                    PATT(NATT,1)=P(I,1)
1713                    PATT(NATT,2)=P(I,2)
1714                    PATT(NATT,3)=P(I,3)
1715                    PATT(NATT,4)=P(I,4)
1716                    EATT=EATT+P(I,4)
1717 cbz11/11/98
1718 cbz1/25/99
1719 c                   IF (NTP .EQ. 1) THEN
1720 c                      GXAR(NATT) = YP(1, jjtp)
1721 c                   ELSE
1722 c                      GXAR(NATT) = YT(1, jjtp)
1723 c                   END IF
1724 c                   IF (NTP .EQ. 1) THEN
1725 c                      GYAR(NATT) = YP(2, jjtp)
1726 c                   ELSE
1727 c                      GYAR(NATT) = YT(2, jjtp)
1728 c                   END IF
1729                    IF (NTP .EQ. 1) THEN
1730                       LSG = jjtp
1731                    ELSE
1732                       LSG = jjtp + NSP
1733                    END IF
1734                    GXAR(NATT) = sngl(ZT1(LSG))
1735                    GYAR(NATT) = sngl(ZT2(LSG))
1736                    GZAR(NATT) = sngl(ZT3(LSG))
1737                    FTAR(NATT) = sngl(ATAUI(LSG))
1738 cbz1/25/99end
1739                    ITYPAR(NATT) = K(I, 2)
1740                    PXAR(NATT) = P(I, 1)
1741                    PYAR(NATT) = P(I, 2)
1742                    PZAR(NATT) = P(I, 3)
1743                    PEAR(NATT) = P(I, 4)
1744                    XMAR(NATT) = P(I, 5)
1745 cbz11/11/98end
1746
1747 390                CONTINUE 
1748 400           CONTINUE
1749 C     ********Fragment the q-qq related string systems
1750         ENDIF
1751
1752         DO 450 I=1,NDR
1753            NATT=NATT+1
1754            KATT(NATT,1)=KFDR(I)
1755            KATT(NATT,2)=40
1756            KATT(NATT,3)=0
1757            PATT(NATT,1)=PDR(I,1)
1758            PATT(NATT,2)=PDR(I,2)
1759            PATT(NATT,3)=PDR(I,3)
1760            PATT(NATT,4)=PDR(I,4)
1761            EATT=EATT+PDR(I,4)
1762 clin-11/11/03     set direct photons positions and time at formation:
1763            GXAR(NATT) = rtdr(I,1)
1764            GYAR(NATT) = rtdr(I,2)
1765            GZAR(NATT) = 0.
1766            FTAR(NATT) = 0.
1767            ITYPAR(NATT) =KATT(NATT,1) 
1768            PXAR(NATT) = PATT(NATT,1)
1769            PYAR(NATT) = PATT(NATT,2)
1770            PZAR(NATT) = PATT(NATT,3)
1771            PEAR(NATT) = PATT(NATT,4)
1772            XMAR(NATT) = PDR(I,5)
1773  450    CONTINUE
1774
1775 C                        ********store the direct-produced particles
1776 C
1777
1778 clin-4/19/01 soft3:
1779  565    continue
1780
1781         DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
1782         IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
1783      &     .AND.IHPR2(21).EQ.0) THEN
1784          IF(IHPR2(10).NE.0) 
1785      &        WRITE(6,*) 'Energy not conserved, repeat the event'
1786 c                call lulist(1)
1787          write(6,*) 'violated:EATT,NATT,B=',EATT,NATT,bimp
1788          GO TO 50
1789         ENDIF
1790 c        write(6,*) 'satisfied:EATT,NATT,B=',EATT,NATT,bimp
1791 c        write(6,*) ' '
1792
1793         RETURN
1794         END
1795 C
1796 C
1797 C
1798         SUBROUTINE HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
1799         CHARACTER FRAME*4,PROJ*4,TARG*4,EFRAME*4
1800         DOUBLE PRECISION  DD1,DD2,DD3,DD4
1801         COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
1802 cc      SAVE /HSTRNG/
1803         COMMON/hjcrdn/YP(3,300),YT(3,300)
1804 cc      SAVE /hjcrdn/
1805         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1806 cc      SAVE /HPARNT/
1807         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
1808 cc      SAVE /HIJDAT/
1809         COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1810 cc      SAVE /LUDAT1A/
1811         EXTERNAL FNKICK,FNKC2,FNSTRU,FNSTRM,FNSTRS
1812         SAVE   
1813
1814         CALL TITLE
1815         IHNT2(1)=IAP
1816         IHNT2(2)=IZP
1817         IHNT2(3)=IAT
1818         IHNT2(4)=IZT
1819         IHNT2(5)=0
1820         IHNT2(6)=0
1821 C
1822         HINT1(8)=MAX(ULMASS(2112),ULMASS(2212))
1823         HINT1(9)=HINT1(8)
1824 C
1825         IF(PROJ.NE.'A') THEN
1826                 IF(PROJ.EQ.'P') THEN
1827                     IHNT2(5)=2212
1828                 ELSE IF(PROJ.EQ.'PBAR') THEN 
1829                     IHNT2(5)=-2212
1830                 ELSE IF(PROJ.EQ.'PI+') THEN
1831                     IHNT2(5)=211
1832                 ELSE IF(PROJ.EQ.'PI-') THEN
1833                     IHNT2(5)=-211
1834                 ELSE IF(PROJ.EQ.'K+') THEN
1835                     IHNT2(5)=321
1836                 ELSE IF(PROJ.EQ.'K-') THEN
1837                     IHNT2(5)=-321
1838                 ELSE IF(PROJ.EQ.'N') THEN
1839                     IHNT2(5)=2112
1840                 ELSE IF(PROJ.EQ.'NBAR') THEN
1841                     IHNT2(5)=-2112
1842                 ELSE
1843                     WRITE(6,*) PROJ, 'wrong or unavailable proj name'
1844                     STOP
1845                 ENDIF
1846                 HINT1(8)=ULMASS(IHNT2(5))
1847         ENDIF
1848         IF(TARG.NE.'A') THEN
1849                 IF(TARG.EQ.'P') THEN
1850                     IHNT2(6)=2212
1851                 ELSE IF(TARG.EQ.'PBAR') THEN 
1852                     IHNT2(6)=-2212
1853                 ELSE IF(TARG.EQ.'PI+') THEN
1854                     IHNT2(6)=211
1855                 ELSE IF(TARG.EQ.'PI-') THEN
1856                     IHNT2(6)=-211
1857                 ELSE IF(TARG.EQ.'K+') THEN
1858                     IHNT2(6)=321
1859                 ELSE IF(TARG.EQ.'K-') THEN
1860                     IHNT2(6)=-321
1861                 ELSE IF(TARG.EQ.'N') THEN
1862                     IHNT2(6)=2112
1863                 ELSE IF(TARG.EQ.'NBAR') THEN
1864                     IHNT2(6)=-2112
1865                 ELSE
1866                     WRITE(6,*) TARG,'wrong or unavailable targ name'
1867                     STOP
1868                 ENDIF
1869                 HINT1(9)=ULMASS(IHNT2(6))
1870         ENDIF
1871
1872 C...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
1873         IF(IHPR2(12).GT.0) THEN
1874         CALL LUGIVE('MDCY(C221,1)=0')
1875 clin-11/07/00 no K* decays:
1876         CALL LUGIVE('MDCY(C313,1)=0')
1877         CALL LUGIVE('MDCY(C-313,1)=0')
1878         CALL LUGIVE('MDCY(C323,1)=0')
1879         CALL LUGIVE('MDCY(C-323,1)=0')
1880 clin-1/04/01 no K0 and K0bar decays so K0L and K0S do not appear,
1881 c     this way the K/Kbar difference is accounted for exactly:
1882         CALL LUGIVE('MDCY(C311,1)=0')
1883         CALL LUGIVE('MDCY(C-311,1)=0')
1884 clin-11/08/00 no Delta decays:
1885         CALL LUGIVE('MDCY(C1114,1)=0')
1886         CALL LUGIVE('MDCY(C2114,1)=0')
1887         CALL LUGIVE('MDCY(C2214,1)=0')
1888         CALL LUGIVE('MDCY(C2224,1)=0')
1889         CALL LUGIVE('MDCY(C-1114,1)=0')
1890         CALL LUGIVE('MDCY(C-2114,1)=0')
1891         CALL LUGIVE('MDCY(C-2214,1)=0')
1892         CALL LUGIVE('MDCY(C-2224,1)=0')
1893 clin-11/07/00-end
1894 cbz12/4/98
1895         CALL LUGIVE('MDCY(C213,1)=0')
1896         CALL LUGIVE('MDCY(C-213,1)=0')
1897         CALL LUGIVE('MDCY(C113,1)=0')
1898         CALL LUGIVE('MDCY(C223,1)=0')
1899         CALL LUGIVE('MDCY(C333,1)=0')
1900 cbz12/4/98end
1901         CALL LUGIVE('MDCY(C111,1)=0')
1902         CALL LUGIVE('MDCY(C310,1)=0')
1903         CALL LUGIVE('MDCY(C411,1)=0;MDCY(C-411,1)=0')
1904         CALL LUGIVE('MDCY(C421,1)=0;MDCY(C-421,1)=0')
1905         CALL LUGIVE('MDCY(C431,1)=0;MDCY(C-431,1)=0')
1906         CALL LUGIVE('MDCY(C511,1)=0;MDCY(C-511,1)=0')
1907         CALL LUGIVE('MDCY(C521,1)=0;MDCY(C-521,1)=0')
1908         CALL LUGIVE('MDCY(C531,1)=0;MDCY(C-531,1)=0')
1909         CALL LUGIVE('MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
1910         CALL LUGIVE('MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
1911         CALL LUGIVE('MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
1912         CALL LUGIVE('MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
1913         CALL LUGIVE('MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
1914         CALL LUGIVE('MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
1915         CALL LUGIVE('MDCY(C3334,1)=0;MDCY(C-3334,1)=0')
1916         ENDIF
1917         MSTU(12)=0
1918         MSTU(21)=1
1919         IF(IHPR2(10).EQ.0) THEN
1920                 MSTU(22)=0
1921                 MSTU(25)=0
1922                 MSTU(26)=0
1923         ENDIF
1924
1925 clin    parj(41) and (42) are a, b parameters in Lund, read from input.ampt:
1926 c        PARJ(41)=HIPR1(3)
1927 c        PARJ(42)=HIPR1(4)
1928 c        PARJ(41)=2.2
1929 c        PARJ(42)=0.5
1930
1931 clin  2 popcorn parameters read from input.ampt:
1932 c        IHPR2(11) = 3
1933 c        PARJ(5) = 0.5
1934         MSTJ(12)=IHPR2(11)
1935
1936 clin  parj(21) gives the mean gaussian width for hadron Pt:
1937         PARJ(21)=HIPR1(2)
1938 clin  parj(2) is gamma_s=P(s)/P(u), kappa propto 1/b/(2+a) assumed.
1939         rkp=HIPR1(4)*(2+HIPR1(3))/PARJ(42)/(2+PARJ(41))
1940         PARJ(2)=PARJ(2)**(1./rkp)
1941         PARJ(21)=PARJ(21)*sqrt(rkp)
1942 clin-10/31/00 update when string tension is changed:
1943         HIPR1(2)=PARJ(21)
1944
1945 C                        ******** set up for jetset
1946         IF(FRAME.EQ.'LAB') THEN
1947            DD1=dble(EFRM)
1948            DD2=dble(HINT1(8))
1949            DD3=dble(HINT1(9))
1950            HINT1(1)=SQRT(HINT1(8)**2+2.0*HINT1(9)*EFRM+HINT1(9)**2)
1951            DD4=DSQRT(DD1**2-DD2**2)/(DD1+DD3)
1952            HINT1(2)=sngl(DD4)
1953            HINT1(3)=0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1954            DD4=DSQRT(DD1**2-DD2**2)/DD1
1955            HINT1(4)=0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1956            HINT1(5)=0.0
1957            HINT1(6)=EFRM
1958            HINT1(7)=HINT1(9)
1959         ELSE IF(FRAME.EQ.'CMS') THEN
1960            HINT1(1)=EFRM
1961            HINT1(2)=0.0
1962            HINT1(3)=0.0
1963            DD1=dble(HINT1(1))
1964            DD2=dble(HINT1(8))
1965            DD3=dble(HINT1(9))
1966            DD4=DSQRT(1.D0-4.D0*DD2**2/DD1**2)
1967            HINT1(4)=0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1968            DD4=DSQRT(1.D0-4.D0*DD3**2/DD1**2)
1969            HINT1(5)=-0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1970            HINT1(6)=HINT1(1)/2.0
1971            HINT1(7)=HINT1(1)/2.0
1972         ENDIF
1973 C                ********define Lorentz transform to lab frame
1974 c
1975 C                ********calculate the cross sections involved with
1976 C                        nucleon collisions.
1977         IF(IHNT2(1).GT.1) THEN
1978                 CALL HIJWDS(IHNT2(1),1,RMAX)
1979                 HIPR1(34)=RMAX
1980 C                        ********set up Wood-Sax distr for proj.
1981         ENDIF
1982         IF(IHNT2(3).GT.1) THEN
1983                 CALL HIJWDS(IHNT2(3),2,RMAX)
1984                 HIPR1(35)=RMAX
1985 C                        ********set up Wood-Sax distr for  targ.
1986         ENDIF
1987 C
1988 C
1989         I=0
1990 20        I=I+1
1991         IF(I.EQ.10) GO TO 30
1992         IF(HIDAT0(10,I).LE.HINT1(1)) GO TO 20
1993 30        IF(I.EQ.1) I=2
1994         DO 40 J=1,9
1995            HIDAT(J)=HIDAT0(J,I-1)+(HIDAT0(J,I)-HIDAT0(J,I-1))
1996      &          *(HINT1(1)-HIDAT0(10,I-1))/(HIDAT0(10,I)-HIDAT0(10,I-1))
1997 40        CONTINUE
1998         HIPR1(31)=HIDAT(5)
1999         HIPR1(30)=2.0*HIDAT(5)
2000 C
2001 C
2002         CALL HIJCRS
2003 C
2004         IF(IHPR2(5).NE.0) THEN
2005                 CALL HIFUN(3,0.0,36.0,FNKICK)
2006 C                ********booking for generating pt**2 for pt kick
2007         ENDIF
2008         CALL HIFUN(7,0.0,6.0,FNKC2)
2009         CALL HIFUN(4,0.0,1.0,FNSTRU)
2010         CALL HIFUN(5,0.0,1.0,FNSTRM)
2011         CALL HIFUN(6,0.0,1.0,FNSTRS)
2012 C                ********booking for x distribution of valence quarks
2013         EFRAME='Ecm'
2014         IF(FRAME.EQ.'LAB') EFRAME='Elab'
2015         WRITE(6,100) EFRAME,EFRM,PROJ,IHNT2(1),IHNT2(2),
2016      &               TARG,IHNT2(3),IHNT2(4) 
2017 100        FORMAT(
2018      &        10X,'**************************************************'/
2019      &        10X,'*',48X,'*'/
2020      &        10X,'*         HIJING has been initialized at         *'/
2021      &        10X,'*',13X,A4,'= ',F10.2,' GeV/n',13X,'*'/
2022      &        10X,'*',48X,'*'/
2023      &        10X,'*',8X,'for ',
2024      &        A4,'(',I3,',',I3,')',' + ',A4,'(',I3,',',I3,')',7X,'*'/
2025      &        10X,'**************************************************')
2026         RETURN
2027         END
2028 C
2029 C
2030 C
2031         FUNCTION FNKICK(X)
2032         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2033 cc      SAVE /HPARNT/
2034         SAVE   
2035         FNKICK=1.0/(X+HIPR1(19)**2)/(X+HIPR1(20)**2)
2036      &                /(1+EXP((SQRT(X)-HIPR1(20))/0.4))
2037         RETURN
2038         END
2039 C
2040 C
2041         FUNCTION FNKC2(X)
2042         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2043 cc      SAVE /HPARNT/
2044         SAVE   
2045         FNKC2=X*EXP(-2.0*X/HIPR1(42))
2046         RETURN
2047         END
2048 C
2049 C
2050 C
2051         FUNCTION FNSTRU(X)
2052         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2053 cc      SAVE /HPARNT/
2054         SAVE   
2055         FNSTRU=(1.0-X)**HIPR1(44)/
2056      &                (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
2057         RETURN
2058         END
2059 C
2060 C
2061 C
2062         FUNCTION FNSTRM(X)
2063         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2064 cc      SAVE /HPARNT/
2065         SAVE   
2066         FNSTRM=1.0/((1.0-X)**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
2067      &          /(X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
2068         RETURN
2069         END
2070 C
2071 C
2072         FUNCTION FNSTRS(X)
2073         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2074 cc      SAVE /HPARNT/
2075         SAVE   
2076         FNSTRS=(1.0-X)**HIPR1(47)/
2077      &                (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(48)
2078         RETURN
2079         END
2080 C
2081 C
2082 C
2083 C
2084         SUBROUTINE HBOOST
2085               IMPLICIT DOUBLE PRECISION(D)  
2086               COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5) 
2087 cc      SAVE /LUJETSA/ 
2088               COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
2089 cc      SAVE /LUDAT1A/ 
2090         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2091 cc      SAVE /HPARNT/
2092         SAVE   
2093         DO 100 I=1,N
2094            DBETA=dble(P(I,3)/P(I,4))
2095            IF(ABS(DBETA).GE.1.D0) THEN
2096               DB=dble(HINT1(2))
2097               IF(DB.GT.0.99999999D0) THEN 
2098 C                ********Rescale boost vector if too close to unity. 
2099                  WRITE(6,*) '(HIBOOT:) boost vector too large' 
2100                  DB=0.99999999D0
2101               ENDIF 
2102               DGA=1D0/SQRT(1D0-DB**2)
2103               DP3=dble(P(I,3))
2104               DP4=dble(P(I,4))
2105               P(I,3)=sngl((DP3+DB*DP4)*DGA)
2106               P(I,4)=sngl((DP4+DB*DP3)*DGA)
2107               GO TO 100
2108            ENDIF
2109            Y=0.5*sngl(DLOG((1.D0+DBETA)/(1.D0-DBETA)))
2110            AMT=SQRT(P(I,1)**2+P(I,2)**2+P(I,5)**2)
2111            P(I,3)=AMT*SINH(Y+HINT1(3))
2112            P(I,4)=AMT*COSH(Y+HINT1(3))
2113 100        CONTINUE
2114         RETURN
2115         END
2116 C
2117 C
2118 C
2119 C
2120         SUBROUTINE QUENCH(JPJT,NTP)
2121         PARAMETER (MAXSTR=150001)
2122         DIMENSION RDP(300),LQP(300),RDT(300),LQT(300)
2123         COMMON/hjcrdn/YP(3,300),YT(3,300)
2124 cc      SAVE /hjcrdn/
2125         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2126 cc      SAVE /HPARNT/
2127 C
2128         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2129      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2130      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
2131      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2132      &                PJTE(300,500),PJTM(300,500)
2133 cc      SAVE /HJJET1/
2134         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2135      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2136      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2137 cc      SAVE /HJJET2/
2138         COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2139 cc      SAVE /HSTRNG/
2140       COMMON/RNDF77/NSEED
2141 cc      SAVE /RNDF77/
2142         SAVE   
2143 C
2144 c     Uzhi:
2145         BB=HINT1(19)
2146         PHI=HINT1(20)
2147         BBX=BB*COS(PHI)
2148         BBY=BB*SIN(PHI)
2149 c
2150         IF(NTP.EQ.2) GO TO 400
2151         IF(NTP.EQ.3) GO TO 2000 
2152 C*******************************************************
2153 C Jet interaction for proj jet in the direction PHIP
2154 C******************************************************
2155 C
2156         IF(NFP(JPJT,7).NE.1) RETURN
2157
2158         JP=JPJT
2159         DO 290 I=1,NPJ(JP)
2160            PTJET0=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2)
2161            IF(PTJET0.LE.HIPR1(11)) GO TO 290
2162            PTOT=SQRT(PTJET0*PTJET0+PJPZ(JP,I)**2)
2163            IF(PTOT.LT.HIPR1(8)) GO TO 290
2164            PHIP=ULANGL(PJPX(JP,I),PJPY(JP,I))
2165 C******* find the wounded proj which can interact with jet***
2166            KP=0
2167            DO 100 I2=1,IHNT2(1)
2168               IF(NFP(I2,5).NE.3 .OR. I2.EQ.JP) GO TO 100
2169               DX=YP(1,I2)-YP(1,JP)
2170               DY=YP(2,I2)-YP(2,JP)
2171               PHI=ULANGL(DX,DY)
2172               DPHI=ABS(PHI-PHIP)
2173 c     Uzhi:
2174               IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2175               IF(DPHI.GE.HIPR1(40)/2.0) GO TO 100
2176               RD0=SQRT(DX*DX+DY*DY)
2177               IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 100
2178               KP=KP+1
2179               LQP(KP)=I2
2180               RDP(KP)=COS(DPHI)*RD0
2181  100           CONTINUE
2182 C*******        rearrange according decending rd************
2183            DO 110 I2=1,KP-1
2184               DO 110 J2=I2+1,KP
2185                  IF(RDP(I2).LT.RDP(J2)) GO TO 110
2186                  RD=RDP(I2)
2187                  LQ=LQP(I2)
2188                  RDP(I2)=RDP(J2)
2189                  LQP(I2)=LQP(J2)
2190                  RDP(J2)=RD
2191                  LQP(J2)=LQ
2192  110              CONTINUE
2193 C****** find wounded targ which can interact with jet********
2194               KT=0
2195               DO 120 I2=1,IHNT2(3)
2196                  IF(NFT(I2,5).NE.3) GO TO 120
2197                  DX=YT(1,I2)-YP(1,JP)-BBX
2198                  DY=YT(2,I2)-YP(2,JP)-BBY
2199                  PHI=ULANGL(DX,DY)
2200                  DPHI=ABS(PHI-PHIP)
2201 c     Uzhi:
2202                  IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2203                  IF(DPHI.GT.HIPR1(40)/2.0) GO TO 120
2204                  RD0=SQRT(DX*DX+DY*DY)
2205                  IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 120
2206                  KT=KT+1
2207                  LQT(KT)=I2
2208                  RDT(KT)=COS(DPHI)*RD0
2209  120              CONTINUE
2210 C*******        rearrange according decending rd************
2211               DO 130 I2=1,KT-1
2212                  DO 130 J2=I2+1,KT
2213                     IF(RDT(I2).LT.RDT(J2)) GO TO 130
2214                     RD=RDT(I2)
2215                     LQ=LQT(I2)
2216                     RDT(I2)=RDT(J2)
2217                     LQT(I2)=LQT(J2)
2218                     RDT(J2)=RD
2219                     LQT(J2)=LQ
2220  130                 CONTINUE
2221                 
2222                  MP=0
2223                  MT=0
2224                  R0=0.0
2225                  NQ=0
2226                  DP=0.0
2227                  PTOT=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2+PJPZ(JP,I)**2)
2228                  V1=PJPX(JP,I)/PTOT
2229                  V2=PJPY(JP,I)/PTOT
2230                  V3=PJPZ(JP,I)/PTOT
2231
2232  200                 RN=RANART(NSEED)
2233  210                 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 290
2234                  IF(MT.GE.KT) GO TO 220
2235                  IF(MP.GE.KP) GO TO 240
2236                  IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 240
2237  220                 MP=MP+1
2238                  DRR=RDP(MP)-R0
2239                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
2240                  DP=DRR*HIPR1(14)
2241                  IF(KFPJ(JP,I).NE.21) DP=0.5*DP
2242 C        ********string tension of quark jet is 0.5 of gluon's 
2243                  IF(DP.LE.0.2) GO TO 210
2244                  IF(PTOT.LE.0.4) GO TO 290
2245                  IF(PTOT.LE.DP) DP=PTOT-0.2
2246                  DE=DP
2247
2248                  IF(KFPJ(JP,I).NE.21) THEN
2249                     PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
2250      &                   +PP(LQP(MP),3)**2
2251                     DE=SQRT(PJPM(JP,I)**2+PTOT**2)
2252      &                        -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
2253                     ERSHU=(PP(LQP(MP),4)+DE-DP)**2
2254                     AMSHU=ERSHU-PRSHU
2255                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
2256                     PP(LQP(MP),4)=SQRT(ERSHU)
2257                     PP(LQP(MP),5)=SQRT(AMSHU)
2258                  ENDIF
2259 C                ********reshuffle the energy when jet has mass
2260                  R0=RDP(MP)
2261                  DP1=DP*V1
2262                  DP2=DP*V2
2263                  DP3=DP*V3
2264 C                ********momentum and energy transfer from jet
2265                  
2266                  NPJ(LQP(MP))=NPJ(LQP(MP))+1
2267                  KFPJ(LQP(MP),NPJ(LQP(MP)))=21
2268                  PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
2269                  PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
2270                  PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
2271                  PJPE(LQP(MP),NPJ(LQP(MP)))=DP
2272                  PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
2273                  GO TO 260
2274
2275  240                 MT=MT+1
2276                  DRR=RDT(MT)-R0
2277                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
2278                  DP=DRR*HIPR1(14)
2279                  IF(DP.LE.0.2) GO TO 210
2280                  IF(PTOT.LE.0.4) GO TO 290
2281                  IF(PTOT.LE.DP) DP=PTOT-0.2
2282                  DE=DP
2283
2284                  IF(KFPJ(JP,I).NE.21) THEN
2285                     PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
2286      &                   +PT(LQT(MT),3)**2
2287                     DE=SQRT(PJPM(JP,I)**2+PTOT**2)
2288      &                        -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
2289                     ERSHU=(PT(LQT(MT),4)+DE-DP)**2
2290                     AMSHU=ERSHU-PRSHU
2291                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
2292                     PT(LQT(MT),4)=SQRT(ERSHU)
2293                     PT(LQT(MT),5)=SQRT(AMSHU)
2294                  ENDIF
2295 C                ********reshuffle the energy when jet has mass
2296
2297                  R0=RDT(MT)
2298                  DP1=DP*V1
2299                  DP2=DP*V2
2300                  DP3=DP*V3
2301 C                ********momentum and energy transfer from jet
2302                  NTJ(LQT(MT))=NTJ(LQT(MT))+1
2303                  KFTJ(LQT(MT),NTJ(LQT(MT)))=21
2304                  PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
2305                  PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
2306                  PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
2307                  PJTE(LQT(MT),NTJ(LQT(MT)))=DP
2308                  PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
2309
2310  260                 PJPX(JP,I)=(PTOT-DP)*V1
2311                  PJPY(JP,I)=(PTOT-DP)*V2
2312                  PJPZ(JP,I)=(PTOT-DP)*V3
2313                  PJPE(JP,I)=PJPE(JP,I)-DE
2314
2315                  PTOT=PTOT-DP
2316                  NQ=NQ+1
2317                  GO TO 200
2318  290              CONTINUE
2319
2320               RETURN
2321
2322 C*******************************************************
2323 C Jet interaction for target jet in the direction PHIT
2324 C******************************************************
2325 C
2326 C******* find the wounded proj which can interact with jet***
2327
2328  400              IF(NFT(JPJT,7).NE.1) RETURN
2329               JT=JPJT
2330               DO 690 I=1,NTJ(JT)
2331                  PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
2332                  IF(PTJET0.LE.HIPR1(11)) GO TO 690
2333                  PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
2334                  IF(PTOT.LT.HIPR1(8)) GO TO 690
2335                  PHIT=ULANGL(PJTX(JT,I),PJTY(JT,I))
2336                  KP=0
2337                  DO 500 I2=1,IHNT2(1)
2338                     IF(NFP(I2,5).NE.3) GO TO 500
2339                     DX=YP(1,I2)+BBX-YT(1,JT)
2340                     DY=YP(2,I2)+BBY-YT(2,JT)
2341                     PHI=ULANGL(DX,DY)
2342                     DPHI=ABS(PHI-PHIT)
2343 c     Uzhi:
2344                     IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2345                     IF(DPHI.GT.HIPR1(40)/2.0) GO TO 500
2346                     RD0=SQRT(DX*DX+DY*DY)
2347                     IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 500
2348                     KP=KP+1
2349                     LQP(KP)=I2
2350                     RDP(KP)=COS(DPHI)*RD0
2351  500                 CONTINUE
2352 C*******        rearrange according to decending rd************
2353                  DO 510 I2=1,KP-1
2354                     DO 510 J2=I2+1,KP
2355                        IF(RDP(I2).LT.RDP(J2)) GO TO 510
2356                        RD=RDP(I2)
2357                        LQ=LQP(I2)
2358                        RDP(I2)=RDP(J2)
2359                        LQP(I2)=LQP(J2)
2360                        RDP(J2)=RD
2361                        LQP(J2)=LQ
2362  510                    CONTINUE
2363 C****** find wounded targ which can interact with jet********
2364                     KT=0
2365                     DO 520 I2=1,IHNT2(3)
2366                        IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 520
2367                        DX=YT(1,I2)-YT(1,JT)
2368                        DY=YT(2,I2)-YT(2,JT)
2369                        PHI=ULANGL(DX,DY)
2370                        DPHI=ABS(PHI-PHIT)
2371 c     Uzhi:
2372                        IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2373                        IF(DPHI.GT.HIPR1(40)/2.0) GO TO 520
2374                        RD0=SQRT(DX*DX+DY*DY)
2375                        IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 520
2376                        KT=KT+1
2377                        LQT(KT)=I2
2378                        RDT(KT)=COS(DPHI)*RD0
2379  520                    CONTINUE
2380 C*******        rearrange according to decending rd************
2381                     DO 530 I2=1,KT-1
2382                        DO 530 J2=I2+1,KT
2383                           IF(RDT(I2).LT.RDT(J2)) GO TO 530
2384                           RD=RDT(I2)
2385                           LQ=LQT(I2)
2386                           RDT(I2)=RDT(J2)
2387                           LQT(I2)=LQT(J2)
2388                           RDT(J2)=RD
2389                           LQT(J2)=LQ
2390  530                       CONTINUE
2391                        
2392                        MP=0
2393                        MT=0
2394                        NQ=0
2395                        DP=0.0
2396                        R0=0.0
2397                 PTOT=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2+PJTZ(JT,I)**2)
2398                 V1=PJTX(JT,I)/PTOT
2399                 V2=PJTY(JT,I)/PTOT
2400                 V3=PJTZ(JT,I)/PTOT
2401
2402  600                RN=RANART(NSEED)
2403  610                IF(MT.GE.KT .AND. MP.GE.KP) GO TO 690
2404                 IF(MT.GE.KT) GO TO 620
2405                 IF(MP.GE.KP) GO TO 640
2406                 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 640
2407 620                MP=MP+1
2408                 DRR=RDP(MP)-R0
2409                 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
2410                 DP=DRR*HIPR1(14)
2411                 IF(KFTJ(JT,I).NE.21) DP=0.5*DP
2412 C        ********string tension of quark jet is 0.5 of gluon's 
2413                 IF(DP.LE.0.2) GO TO 610
2414                 IF(PTOT.LE.0.4) GO TO 690
2415                 IF(PTOT.LE.DP) DP=PTOT-0.2
2416                 DE=DP
2417 C
2418                 IF(KFTJ(JT,I).NE.21) THEN
2419                    PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
2420      &                   +PP(LQP(MP),3)**2
2421                    DE=SQRT(PJTM(JT,I)**2+PTOT**2)
2422      &                     -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
2423                    ERSHU=(PP(LQP(MP),4)+DE-DP)**2
2424                    AMSHU=ERSHU-PRSHU
2425                    IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
2426                    PP(LQP(MP),4)=SQRT(ERSHU)
2427                    PP(LQP(MP),5)=SQRT(AMSHU)
2428                 ENDIF
2429 C                ********reshuffle the energy when jet has mass
2430 C
2431                 R0=RDP(MP)
2432                 DP1=DP*V1
2433                 DP2=DP*V2
2434                 DP3=DP*V3
2435 C                ********momentum and energy transfer from jet
2436                 NPJ(LQP(MP))=NPJ(LQP(MP))+1
2437                 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
2438                 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
2439                 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
2440                 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
2441                 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
2442                 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
2443
2444                 GO TO 660
2445
2446 640                MT=MT+1
2447                 DRR=RDT(MT)-R0
2448                 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
2449                 DP=DRR*HIPR1(14)
2450                 IF(DP.LE.0.2) GO TO 610
2451                 IF(PTOT.LE.0.4) GO TO 690
2452                 IF(PTOT.LE.DP) DP=PTOT-0.2
2453                 DE=DP
2454
2455                 IF(KFTJ(JT,I).NE.21) THEN
2456                    PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
2457      &                   +PT(LQT(MT),3)**2
2458                    DE=SQRT(PJTM(JT,I)**2+PTOT**2)
2459      &                     -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
2460                    ERSHU=(PT(LQT(MT),4)+DE-DP)**2
2461                    AMSHU=ERSHU-PRSHU
2462                    IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
2463                    PT(LQT(MT),4)=SQRT(ERSHU)
2464                    PT(LQT(MT),5)=SQRT(AMSHU)
2465                 ENDIF
2466 C                ********reshuffle the energy when jet has mass
2467
2468                 R0=RDT(MT)
2469                 DP1=DP*V1
2470                 DP2=DP*V2
2471                 DP3=DP*V3
2472 C                ********momentum and energy transfer from jet
2473                 NTJ(LQT(MT))=NTJ(LQT(MT))+1
2474                 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
2475                 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
2476                 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
2477                 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
2478                 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
2479                 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
2480
2481 660                PJTX(JT,I)=(PTOT-DP)*V1
2482                 PJTY(JT,I)=(PTOT-DP)*V2
2483                 PJTZ(JT,I)=(PTOT-DP)*V3
2484                 PJTE(JT,I)=PJTE(JT,I)-DE
2485
2486                 PTOT=PTOT-DP
2487                 NQ=NQ+1
2488                 GO TO 600
2489 690        CONTINUE
2490         RETURN
2491 C********************************************************
2492 C        Q-QBAR jet interaction
2493 C********************************************************
2494 2000        ISG=JPJT
2495         IF(IASG(ISG,3).NE.1) RETURN
2496 C
2497         JP=IASG(ISG,1)
2498         JT=IASG(ISG,2)
2499         XJ=(YP(1,JP)+BBX+YT(1,JT))/2.0
2500         YJ=(YP(2,JP)+BBY+YT(2,JT))/2.0
2501         DO 2690 I=1,NJSG(ISG)
2502            PTJET0=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2)
2503            IF(PTJET0.LE.HIPR1(11).OR.PESG(ISG,I).LT.HIPR1(1))
2504      &            GO TO 2690
2505            PTOT=SQRT(PTJET0*PTJET0+PZSG(ISG,I)**2)
2506            IF(PTOT.LT.MAX(HIPR1(1),HIPR1(8))) GO TO 2690
2507            PHIQ=ULANGL(PXSG(ISG,I),PYSG(ISG,I))
2508            KP=0
2509            DO 2500 I2=1,IHNT2(1)
2510               IF(NFP(I2,5).NE.3.OR.I2.EQ.JP) GO TO 2500
2511               DX=YP(1,I2)+BBX-XJ
2512               DY=YP(2,I2)+BBY-YJ
2513               PHI=ULANGL(DX,DY)
2514               DPHI=ABS(PHI-PHIQ)
2515 c     Uzhi:
2516               IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2517               IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2500
2518               RD0=SQRT(DX*DX+DY*DY)
2519               IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2500
2520               KP=KP+1
2521               LQP(KP)=I2
2522               RDP(KP)=COS(DPHI)*RD0
2523  2500           CONTINUE
2524 C*******        rearrange according to decending rd************
2525            DO 2510 I2=1,KP-1
2526               DO 2510 J2=I2+1,KP
2527                  IF(RDP(I2).LT.RDP(J2)) GO TO 2510
2528                  RD=RDP(I2)
2529                  LQ=LQP(I2)
2530                  RDP(I2)=RDP(J2)
2531                  LQP(I2)=LQP(J2)
2532                  RDP(J2)=RD
2533                  LQP(J2)=LQ
2534  2510              CONTINUE
2535 C****** find wounded targ which can interact with jet********
2536               KT=0
2537               DO 2520 I2=1,IHNT2(3)
2538                  IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 2520
2539                  DX=YT(1,I2)-XJ
2540                  DY=YT(2,I2)-YJ
2541                  PHI=ULANGL(DX,DY)
2542                  DPHI=ABS(PHI-PHIQ)
2543 c     Uzhi:
2544                  IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2545                  IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2520
2546                  RD0=SQRT(DX*DX+DY*DY)
2547                  IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2520
2548                  KT=KT+1
2549                  LQT(KT)=I2
2550                  RDT(KT)=COS(DPHI)*RD0
2551  2520              CONTINUE
2552 C*******        rearrange according to decending rd************
2553               DO 2530 I2=1,KT-1
2554                  DO 2530 J2=I2+1,KT
2555                     IF(RDT(I2).LT.RDT(J2)) GO TO 2530
2556                     RD=RDT(I2)
2557                     LQ=LQT(I2)
2558                     RDT(I2)=RDT(J2)
2559                     LQT(I2)=LQT(J2)
2560                     RDT(J2)=RD
2561                     LQT(J2)=LQ
2562  2530                 CONTINUE
2563                 
2564                  MP=0
2565                  MT=0
2566                  NQ=0
2567                  DP=0.0
2568                  R0=0.0
2569                  PTOT=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2
2570      &                +PZSG(ISG,I)**2)
2571                  V1=PXSG(ISG,I)/PTOT
2572                  V2=PYSG(ISG,I)/PTOT
2573                  V3=PZSG(ISG,I)/PTOT
2574
2575  2600                 RN=RANART(NSEED)
2576  2610                 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 2690
2577                  IF(MT.GE.KT) GO TO 2620
2578                  IF(MP.GE.KP) GO TO 2640
2579                  IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 2640
2580  2620                 MP=MP+1
2581                  DRR=RDP(MP)-R0
2582                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
2583                  DP=DRR*HIPR1(14)/2.0
2584                  IF(DP.LE.0.2) GO TO 2610
2585                  IF(PTOT.LE.0.4) GO TO 2690
2586                  IF(PTOT.LE.DP) DP=PTOT-0.2
2587                  DE=DP
2588 C
2589                  IF(K2SG(ISG,I).NE.21) THEN
2590                     IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
2591                     PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
2592      &                    +PP(LQP(MP),3)**2
2593                     DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
2594      &                       -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
2595                     ERSHU=(PP(LQP(MP),4)+DE-DP)**2
2596                     AMSHU=ERSHU-PRSHU
2597                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
2598                     PP(LQP(MP),4)=SQRT(ERSHU)
2599                     PP(LQP(MP),5)=SQRT(AMSHU)
2600                  ENDIF
2601 C                ********reshuffle the energy when jet has mass
2602 C
2603                  R0=RDP(MP)
2604                  DP1=DP*V1
2605                  DP2=DP*V2
2606                  DP3=DP*V3
2607 C                ********momentum and energy transfer from jet
2608                  NPJ(LQP(MP))=NPJ(LQP(MP))+1
2609                  KFPJ(LQP(MP),NPJ(LQP(MP)))=21
2610                  PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
2611                  PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
2612                  PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
2613                  PJPE(LQP(MP),NPJ(LQP(MP)))=DP
2614                  PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
2615
2616                  GO TO 2660
2617
2618  2640                 MT=MT+1
2619                  DRR=RDT(MT)-R0
2620                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
2621                  DP=DRR*HIPR1(14)
2622                  IF(DP.LE.0.2) GO TO 2610
2623                  IF(PTOT.LE.0.4) GO TO 2690
2624                  IF(PTOT.LE.DP) DP=PTOT-0.2
2625                  DE=DP
2626
2627                  IF(K2SG(ISG,I).NE.21) THEN
2628                     IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
2629                     PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
2630      &                    +PT(LQT(MT),3)**2
2631                     DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
2632      &                       -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
2633                     ERSHU=(PT(LQT(MT),4)+DE-DP)**2
2634                     AMSHU=ERSHU-PRSHU
2635                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
2636                     PT(LQT(MT),4)=SQRT(ERSHU)
2637                     PT(LQT(MT),5)=SQRT(AMSHU)
2638                  ENDIF
2639 C               ********reshuffle the energy when jet has mass
2640
2641                  R0=RDT(MT)
2642                  DP1=DP*V1
2643                  DP2=DP*V2
2644                  DP3=DP*V3
2645 C                ********momentum and energy transfer from jet
2646                  NTJ(LQT(MT))=NTJ(LQT(MT))+1
2647                  KFTJ(LQT(MT),NTJ(LQT(MT)))=21
2648                  PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
2649                  PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
2650                  PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
2651                  PJTE(LQT(MT),NTJ(LQT(MT)))=DP
2652                  PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
2653
2654  2660                 PXSG(ISG,I)=(PTOT-DP)*V1
2655                  PYSG(ISG,I)=(PTOT-DP)*V2
2656                  PZSG(ISG,I)=(PTOT-DP)*V3
2657                  PESG(ISG,I)=PESG(ISG,I)-DE
2658
2659                  PTOT=PTOT-DP
2660                  NQ=NQ+1
2661                  GO TO 2600
2662  2690        CONTINUE
2663         RETURN
2664         END
2665
2666 C
2667 C
2668 C
2669 C
2670         SUBROUTINE HIJFRG(JTP,NTP,IERROR)
2671 C        NTP=1, fragment proj string, NTP=2, targ string, 
2672 C       NTP=3, independent 
2673 C        strings from jets.  JTP is the line number of the string
2674 C*******Fragment all leadng strings of proj and targ**************
2675 C        IHNT2(1)=atomic #, IHNT2(2)=proton #(=-1 if anti-proton)  *
2676 C******************************************************************
2677         PARAMETER (MAXSTR=150001)
2678         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2679 cc      SAVE /HPARNT/
2680         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
2681 cc      SAVE /HIJDAT/
2682         COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2683 cc      SAVE /HSTRNG/
2684         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2685      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2686      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
2687      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2688      &                PJTE(300,500),PJTM(300,500)
2689 cc      SAVE /HJJET1/
2690         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2691      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2692      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2693 cc      SAVE /HJJET2/
2694 C
2695         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
2696 cc      SAVE /LUJETSA/
2697         COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2698 cc      SAVE /LUDAT1A/
2699       COMMON/RNDF77/NSEED
2700 cc      SAVE /RNDF77/
2701 clin-4/11/01 soft:
2702       common/anim/nevent,isoft,isflag,izpc
2703 cc      SAVE /anim/
2704         SAVE   
2705         
2706 cbz3/12/99
2707 c.....set up fragmentation function according to the number of collisions
2708 c.....a wounded nucleon has suffered
2709 c        IF (NTP .EQ. 1) THEN
2710 c           NCOLL = NFP(JTP, 11)
2711 c        ELSE IF (NTP .EQ. 2) THEN
2712 c           NCOLL = NFT(JTP, 11)
2713 c        ELSE IF (NTP .EQ. 3) THEN
2714 c           NCOLL = (NFP(IASG(JTP,1), 11) + NFT(IASG(JTP,2), 11)) / 2
2715 c        END IF
2716 c        IF (NCOLL .LE. 1) THEN
2717 c           PARJ(5) = 0.5
2718 c        ELSE IF (NCOLL .EQ. 2) THEN
2719 c           PARJ(5) = 0.75
2720 c        ELSE IF (NCOLL .EQ. 3) THEN
2721 c           PARJ(5) = 1.17
2722 c        ELSE IF (NCOLL .EQ. 4) THEN
2723 c           PARJ(5) = 2.0
2724 c        ELSE IF (NCOLL .EQ. 5) THEN
2725 c           PARJ(5) = 4.5
2726 c        ELSE IF (NCOLL .GE. 6) THEN
2727 c           PARJ(5) = 49.5
2728 c        END IF
2729 c        PARJ(5) = 0.5
2730 cbz3/12/99 end
2731
2732         IERROR=0
2733         CALL LUEDIT(0)
2734         N=0
2735 C                        ********initialize the document lines
2736         IF(NTP.EQ.3) THEN
2737                 ISG=JTP
2738                 N=NJSG(ISG)
2739                 DO 100 I=1,NJSG(ISG)
2740                    K(I,1)=K1SG(ISG,I)
2741                    K(I,2)=K2SG(ISG,I)
2742                    P(I,1)=PXSG(ISG,I)
2743                    P(I,2)=PYSG(ISG,I)
2744                    P(I,3)=PZSG(ISG,I)
2745                    P(I,4)=PESG(ISG,I)
2746                    P(I,5)=PMSG(ISG,I)
2747  100            CONTINUE
2748
2749 C                IF(IHPR2(1).GT.0) CALL ATTRAD(IERROR)
2750 c                IF(IERROR.NE.0) RETURN
2751 C                CALL LULIST(1)
2752                 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2753              RETURN
2754         ENDIF
2755 C
2756         IF(NTP.EQ.2) GO TO 200
2757         IF(JTP.GT.IHNT2(1))   RETURN
2758         IF(NFP(JTP,5).NE.3.AND.NFP(JTP,3).NE.0
2759      &            .AND.NPJ(JTP).EQ.0.AND.NFP(JTP,10).EQ.0) GO TO 1000
2760         IF(NFP(JTP,15).EQ.-1) THEN
2761                 KF1=NFP(JTP,2)
2762                 KF2=NFP(JTP,1)
2763                 PQ21=PP(JTP,6)
2764                 PQ22=PP(JTP,7)
2765                 PQ11=PP(JTP,8)
2766                 PQ12=PP(JTP,9)
2767                 AM1=PP(JTP,15)
2768                 AM2=PP(JTP,14)
2769         ELSE
2770                 KF1=NFP(JTP,1)
2771                 KF2=NFP(JTP,2)
2772                 PQ21=PP(JTP,8)
2773                 PQ22=PP(JTP,9)
2774                 PQ11=PP(JTP,6)
2775                 PQ12=PP(JTP,7)
2776                 AM1=PP(JTP,14)
2777                 AM2=PP(JTP,15)        
2778         ENDIF
2779
2780 C        ********for NFP(JTP,15)=-1 NFP(JTP,1) IS IN -Z DIRECTION
2781         PB1=PQ11+PQ21
2782         PB2=PQ12+PQ22
2783         PB3=PP(JTP,3)
2784         PECM=PP(JTP,5)
2785         BTZ=PB3/PP(JTP,4)
2786         IF((ABS(PB1-PP(JTP,1)).GT.0.01.OR.
2787      &     ABS(PB2-PP(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
2788      &     WRITE(6,*) '  Pt of Q and QQ do not sum to the total',jtp
2789      &     ,ntp,pq11,pq21,pb1,'*',pq12,pq22,pb2,'*',pp(JTP,1),pp(JTP,2)
2790         GO TO 300
2791
2792 200        IF(JTP.GT.IHNT2(3))  RETURN
2793         IF(NFT(JTP,5).NE.3.AND.NFT(JTP,3).NE.0
2794      &           .AND.NTJ(JTP).EQ.0.AND.NFT(JTP,10).EQ.0) GO TO 1200
2795         IF(NFT(JTP,15).EQ.1) THEN
2796                 KF1=NFT(JTP,1)
2797                 KF2=NFT(JTP,2)
2798                 PQ11=PT(JTP,6)
2799                 PQ12=PT(JTP,7)
2800                 PQ21=PT(JTP,8)
2801                 PQ22=PT(JTP,9)
2802                 AM1=PT(JTP,14)
2803                 AM2=PT(JTP,15)
2804         ELSE
2805                 KF1=NFT(JTP,2)
2806                 KF2=NFT(JTP,1)
2807                 PQ11=PT(JTP,8)
2808                 PQ12=PT(JTP,9)
2809                 PQ21=PT(JTP,6)
2810                 PQ22=PT(JTP,7)
2811                 AM1=PT(JTP,15)
2812                 AM2=PT(JTP,14)
2813         ENDIF        
2814 C        ********for NFT(JTP,15)=1 NFT(JTP,1) IS IN +Z DIRECTION
2815         PB1=PQ11+PQ21
2816         PB2=PQ12+PQ22
2817         PB3=PT(JTP,3)
2818         PECM=PT(JTP,5)
2819         BTZ=PB3/PT(JTP,4)
2820
2821         IF((ABS(PB1-PT(JTP,1)).GT.0.01.OR.
2822      &     ABS(PB2-PT(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
2823      &     WRITE(6,*) '  Pt of Q and QQ do not sum to the total',jtp
2824      &     ,ntp,pq11,pq21,pb1,'*',pq12,pq22,pb2,'*',pt(JTP,1),pt(JTP,2)
2825 300        IF(PECM.LT.HIPR1(1)) THEN
2826            IERROR=1
2827            IF(IHPR2(10).EQ.0) RETURN
2828            WRITE(6,*) ' ECM=',PECM,' energy of the string is too small'
2829 clin:
2830            write (6,*) 'JTP,NTP,pq=',JTP,NTP,pq11,pq12,pq21,pq22
2831            RETURN
2832         ENDIF
2833         AMT=PECM**2+PB1**2+PB2**2
2834         AMT1=AM1**2+PQ11**2+PQ12**2
2835         AMT2=AM2**2+PQ21**2+PQ22**2
2836         PZCM=SQRT(ABS(AMT**2+AMT1**2+AMT2**2-2.0*AMT*AMT1
2837      &       -2.0*AMT*AMT2-2.0*AMT1*AMT2))/2.0/SQRT(AMT)
2838 C                *******PZ of end-partons in c.m. frame of the string
2839         K(1,1)=2
2840         K(1,2)=KF1
2841         P(1,1)=PQ11
2842         P(1,2)=PQ12
2843         P(1,3)=PZCM
2844         P(1,4)=SQRT(AMT1+PZCM**2)
2845         P(1,5)=AM1
2846         K(2,1)=1
2847         K(2,2)=KF2
2848         P(2,1)=PQ21
2849         P(2,2)=PQ22
2850         P(2,3)=-PZCM
2851         P(2,4)=SQRT(AMT2+PZCM**2)
2852         P(2,5)=AM2
2853         N=2
2854 C*****
2855         CALL HIROBO(0.0,0.0,0.0,0.0,BTZ)
2856         JETOT=0
2857         IF((PQ21**2+PQ22**2).GT.(PQ11**2+PQ12**2)) THEN
2858                 PMAX1=P(2,1)
2859                 PMAX2=P(2,2)
2860                 PMAX3=P(2,3)
2861         ELSE
2862                 PMAX1=P(1,1)
2863                 PMAX2=P(1,2)
2864                 PMAX3=P(1,3)
2865         ENDIF
2866         IF(NTP.EQ.1) THEN
2867                 PP(JTP,10)=PMAX1
2868                 PP(JTP,11)=PMAX2
2869                 PP(JTP,12)=PMAX3
2870         ELSE IF(NTP.EQ.2) THEN
2871                 PT(JTP,10)=PMAX1
2872                 PT(JTP,11)=PMAX2
2873                 PT(JTP,12)=PMAX3
2874         ENDIF
2875 C*******************attach produced jets to the leadng partons****
2876         IF(NTP.EQ.1.AND.NPJ(JTP).NE.0) THEN
2877                 JETOT=NPJ(JTP)
2878 C                IF(NPJ(JTP).GE.2) CALL HIJSRT(JTP,1)
2879 C                        ********sort jets in order of y
2880                 IEX=0
2881                 IF((ABS(KF1).GT.1000.AND.KF1.LT.0)
2882      &                        .OR.(ABS(KF1).LT.1000.AND.KF1.GT.0)) IEX=1
2883                 DO 520 I=N,2,-1
2884                 DO 520 J=1,5
2885                         II=NPJ(JTP)+I
2886                         K(II,J)=K(I,J)
2887                         P(II,J)=P(I,J)
2888                         V(II,J)=V(I,J)
2889 520                CONTINUE
2890
2891                 DO 540 I=1,NPJ(JTP)
2892                         DO 542 J=1,5
2893                                 K(I+1,J)=0
2894                                 V(I+1,J)=0
2895 542                        CONTINUE                                
2896                         I0=I
2897 clin-4/12/01:                        IF(IEX.EQ.1) I0=NPJ(JTP)-I+1
2898                         IF(IEX.EQ.1.and.(isoft.ne.2.or.isflag.ne.0))
2899      1 I0=NPJ(JTP)-I+1
2900 C                                ********reverse the order of jets
2901                         KK1=KFPJ(JTP,I0)
2902                         K(I+1,1)=2
2903                         K(I+1,2)=KK1
2904                         IF(KK1.NE.21 .AND. KK1.NE.0)  K(I+1,1)=
2905      &                          1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
2906                         P(I+1,1)=PJPX(JTP,I0)
2907                         P(I+1,2)=PJPY(JTP,I0)
2908                         P(I+1,3)=PJPZ(JTP,I0)
2909                         P(I+1,4)=PJPE(JTP,I0)
2910                         P(I+1,5)=PJPM(JTP,I0)
2911 540                CONTINUE
2912                 N=N+NPJ(JTP)
2913         ELSE IF(NTP.EQ.2.AND.NTJ(JTP).NE.0) THEN
2914                 JETOT=NTJ(JTP)
2915 c                IF(NTJ(JTP).GE.2)  CALL HIJSRT(JTP,2)
2916 C                        ********sort jets in order of y
2917                 IEX=1
2918                 IF((ABS(KF2).GT.1000.AND.KF2.LT.0)
2919      &                        .OR.(ABS(KF2).LT.1000.AND.KF2.GT.0)) IEX=0
2920                 DO 560 I=N,2,-1
2921                 DO 560 J=1,5
2922                         II=NTJ(JTP)+I
2923                         K(II,J)=K(I,J)
2924                         P(II,J)=P(I,J)
2925                         V(II,J)=V(I,J)
2926 560                CONTINUE
2927                 DO 580 I=1,NTJ(JTP)
2928                         DO 582 J=1,5
2929                                 K(I+1,J)=0
2930                                 V(I+1,J)=0
2931 582                        CONTINUE                                
2932                         I0=I
2933 clin-4/12/01:                        IF(IEX.EQ.1) I0=NTJ(JTP)-I+1
2934                         IF(IEX.EQ.1.and.(isoft.ne.2.or.isflag.ne.0))
2935      1 I0=NTJ(JTP)-I+1
2936 C                                ********reverse the order of jets
2937                         KK1=KFTJ(JTP,I0)
2938                         K(I+1,1)=2
2939                         K(I+1,2)=KK1
2940                         IF(KK1.NE.21 .AND. KK1.NE.0) K(I+1,1)=
2941      &                           1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
2942                         P(I+1,1)=PJTX(JTP,I0)
2943                         P(I+1,2)=PJTY(JTP,I0)
2944                         P(I+1,3)=PJTZ(JTP,I0)
2945                         P(I+1,4)=PJTE(JTP,I0)
2946                         P(I+1,5)=PJTM(JTP,I0)
2947 580                CONTINUE
2948                 N=N+NTJ(JTP)
2949         ENDIF
2950         IF(IHPR2(1).GT.0.AND.RANART(NSEED).LE.HIDAT(3)) THEN
2951              HDAT20=HIDAT(2)
2952              HPR150=HIPR1(5)
2953              IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
2954      &                        HIDAT(2)=2.0
2955              IF(HINT1(1).GE.1000.0.AND.JETOT.EQ.0)THEN
2956                 HIDAT(2)=3.0
2957                 HIPR1(5)=5.0
2958              ENDIF
2959              CALL ATTRAD(IERROR)
2960              HIDAT(2)=HDAT20
2961              HIPR1(5)=HPR150
2962         ELSE IF(JETOT.EQ.0.AND.IHPR2(1).GT.0.AND.
2963      &                       HINT1(1).GE.1000.0.AND.
2964      &                RANART(NSEED).LE.0.8) THEN
2965                 HDAT20=HIDAT(2)
2966                 HPR150=HIPR1(5)
2967                 HIDAT(2)=3.0
2968                 HIPR1(5)=5.0
2969              IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
2970      &                        HIDAT(2)=2.0
2971                 CALL ATTRAD(IERROR)
2972                 HIDAT(2)=HDAT20
2973                 HIPR1(5)=HPR150
2974         ENDIF
2975         IF(IERROR.NE.0) RETURN
2976 C                ******** conduct soft radiations
2977 C****************************
2978 C
2979 C
2980 clin-4/11/01 soft:
2981 c        CALL LUEXEC
2982         if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2983
2984         RETURN
2985
2986 1000        N=1
2987         K(1,1)=1
2988                K(1,2)=NFP(JTP,3)
2989         DO 1100 JJ=1,5
2990                        P(1,JJ)=PP(JTP,JJ)
2991 1100                CONTINUE
2992 C                        ********proj remain as a nucleon or delta
2993 clin-4/11/01 soft:
2994 c        CALL LUEXEC
2995         if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2996
2997 C        call lulist(1)
2998         RETURN
2999 C
3000 1200        N=1
3001         K(1,1)=1
3002         K(1,2)=NFT(JTP,3)
3003         DO 1300 JJ=1,5
3004                 P(1,JJ)=PT(JTP,JJ)
3005 1300        CONTINUE
3006 C                        ********targ remain as a nucleon or delta
3007 clin-4/11/01 soft:
3008 c        CALL LUEXEC
3009         if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
3010
3011 C        call lulist(1)
3012         RETURN
3013         END
3014 C
3015 C
3016 C
3017 C
3018 C****************************************************************
3019 C        conduct soft radiation according to dipole approxiamtion
3020 C****************************************************************
3021         SUBROUTINE ATTRAD(IERROR)
3022 C
3023         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3024 cc      SAVE /HPARNT/
3025         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3026 cc      SAVE /HIJDAT/
3027         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3028 cc      SAVE /LUJETSA/
3029       COMMON/RNDF77/NSEED
3030 cc      SAVE /RNDF77/
3031         SAVE   
3032         IERROR=0
3033
3034 C.....S INVARIANT MASS-SQUARED BETWEEN PARTONS I AND I+1......
3035 C.....SM IS THE LARGEST MASS-SQUARED....
3036
3037 40        SM=0.
3038         JL=1
3039         DO 30 I=1,N-1
3040            S=2.*(P(I,4)*P(I+1,4)-P(I,1)*P(I+1,1)-P(I,2)*P(I+1,2)
3041      &                -P(I,3)*P(I+1,3))+P(I,5)**2+P(I+1,5)**2
3042            IF(S.LT.0.) S=0.
3043            WP=SQRT(S)-1.5*(P(I,5)+P(I+1,5))
3044            IF(WP.GT.SM) THEN
3045               PBT1=P(I,1)+P(I+1,1)
3046               PBT2=P(I,2)+P(I+1,2)
3047               PBT3=P(I,3)+P(I+1,3)
3048               PBT4=P(I,4)+P(I+1,4)
3049               BTT=(PBT1**2+PBT2**2+PBT3**2)/PBT4**2
3050               IF(BTT.GE.1.0-1.0E-10) GO TO 30
3051               IF((I.NE.1.OR.I.NE.N-1).AND.
3052      &             (K(I,2).NE.21.AND.K(I+1,2).NE.21)) GO TO 30
3053               JL=I
3054               SM=WP
3055            ENDIF
3056 30        CONTINUE
3057         S=(SM+1.5*(P(JL,5)+P(JL+1,5)))**2
3058               IF(SM.LT.HIPR1(5)) GOTO 2
3059      
3060 C.....MAKE PLACE FOR ONE GLUON.....
3061               IF(JL+1.EQ.N) GOTO 190
3062               DO 160 J=N,JL+2,-1
3063                       K(J+1,1)=K(J,1)
3064                 K(J+1,2)=K(J,2)
3065                       DO 150 M=1,5
3066 150                           P(J+1,M)=P(J,M)
3067 160                   CONTINUE
3068 190           N=N+1
3069      
3070 C.....BOOST TO REST SYSTEM FOR PARTICLES JL AND JL+1.....
3071               P1=P(JL,1)+P(JL+1,1)
3072               P2=P(JL,2)+P(JL+1,2)
3073               P3=P(JL,3)+P(JL+1,3)
3074               P4=P(JL,4)+P(JL+1,4)
3075               BEX=-P1/P4
3076               BEY=-P2/P4
3077               BEZ=-P3/P4
3078         IMIN=JL
3079         IMAX=JL+1
3080               CALL ATROBO(0.,0.,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
3081         IF(IERROR.NE.0) RETURN
3082 C.....ROTATE TO Z-AXIS....
3083               CTH=P(JL,3)/SQRT(P(JL,4)**2-P(JL,5)**2)
3084               IF(ABS(CTH).GT.1.0)  CTH=MAX(-1.,MIN(1.,CTH))
3085               THETA=ACOS(CTH)
3086               PHI=ULANGL(P(JL,1),P(JL,2))
3087               CALL ATROBO(0.,-PHI,0.,0.,0.,IMIN,IMAX,IERROR)
3088               CALL ATROBO(-THETA,0.,0.,0.,0.,IMIN,IMAX,IERROR)
3089      
3090 C.....CREATE ONE GLUON AND ORIENTATE.....
3091      
3092 1        CALL AR3JET(S,X1,X3,JL)
3093               CALL ARORIE(S,X1,X3,JL)                
3094         IF(HIDAT(2).GT.0.0) THEN
3095                  PTG1=SQRT(P(JL,1)**2+P(JL,2)**2)
3096                  PTG2=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
3097                  PTG3=SQRT(P(JL+2,1)**2+P(JL+2,2)**2)
3098            PTG=MAX(PTG1,PTG2,PTG3)
3099            IF(PTG.GT.HIDAT(2)) THEN
3100               FMFACT=EXP(-(PTG**2-HIDAT(2)**2)/HIPR1(2)**2)
3101               IF(RANART(NSEED).GT.FMFACT) GO TO 1
3102            ENDIF
3103         ENDIF
3104 C.....ROTATE AND BOOST BACK.....
3105         IMIN=JL
3106         IMAX=JL+2
3107               CALL ATROBO(THETA,PHI,-BEX,-BEY,-BEZ,IMIN,IMAX,IERROR)
3108         IF(IERROR.NE.0) RETURN
3109 C.....ENUMERATE THE GLUONS.....
3110               K(JL+2,1)=K(JL+1,1)
3111         K(JL+2,2)=K(JL+1,2)
3112         K(JL+2,3)=K(JL+1,3)
3113         K(JL+2,4)=K(JL+1,4)
3114         K(JL+2,5)=K(JL+1,5)
3115               P(JL+2,5)=P(JL+1,5)
3116               K(JL+1,1)=2
3117         K(JL+1,2)=21
3118         K(JL+1,3)=0
3119         K(JL+1,4)=0
3120         K(JL+1,5)=0
3121               P(JL+1,5)=0.
3122 C----THETA FUNCTION DAMPING OF THE EMITTED GLUONS. FOR HADRON-HADRON.
3123 C----R0=VFR(2)
3124 C              IF(VFR(2).GT.0.) THEN
3125 C              PTG=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
3126 C              PTGMAX=WSTRI/2.
3127 C              DOPT=SQRT((4.*PAR(71)*VFR(2))/WSTRI)
3128 C              PTOPT=(DOPT*WSTRI)/(2.*VFR(2))
3129 C              IF(PTG.GT.PTOPT) IORDER=IORDER-1
3130 C              IF(PTG.GT.PTOPT) GOTO 1
3131 C              ENDIF
3132 C-----
3133              IF(SM.GE.HIPR1(5)) GOTO 40
3134
3135 2              K(1,1)=2
3136         K(1,3)=0
3137         K(1,4)=0
3138         K(1,5)=0
3139               K(N,1)=1
3140         K(N,3)=0
3141         K(N,4)=0
3142         K(N,5)=0
3143
3144               RETURN
3145               END
3146
3147
3148         SUBROUTINE AR3JET(S,X1,X3,JL)
3149 C     
3150         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3151 cc      SAVE /HPARNT/
3152         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3153 cc      SAVE /LUJETSA/
3154       COMMON/RNDF77/NSEED
3155 cc      SAVE /RNDF77/
3156         SAVE   
3157 C     
3158         C=1./3.
3159               IF(K(JL,2).NE.21 .AND. K(JL+1,2).NE.21) C=8./27.
3160               EXP1=3
3161               EXP3=3
3162               IF(K(JL,2).NE.21) EXP1=2
3163               IF(K(JL+1,2).NE.21) EXP3=2
3164               A=0.24**2/S
3165               YMA=ALOG(.5/SQRT(A)+SQRT(.25/A-1))
3166               D=4.*C*YMA
3167               SM1=P(JL,5)**2/S
3168               SM3=P(JL+1,5)**2/S
3169               XT2M=(1.-2.*SQRT(SM1)+SM1-SM3)*(1.-2.*SQRT(SM3)-SM1+SM3)
3170               XT2M=MIN(.25,XT2M)
3171               NTRY=0
3172 1             IF(NTRY.EQ.5000) THEN
3173                 X1=.5*(2.*SQRT(SM1)+1.+SM1-SM3)
3174                 X3=.5*(2.*SQRT(SM3)+1.-SM1+SM3)
3175                 RETURN
3176               ENDIF
3177               NTRY=NTRY+1
3178      
3179               XT2=A*(XT2M/A)**(RANART(NSEED)**(1./D))
3180      
3181               YMAX=ALOG(.5/SQRT(XT2)+SQRT(.25/XT2-1.))
3182               Y=(2.*RANART(NSEED)-1.)*YMAX
3183               X1=1.-SQRT(XT2)*EXP(Y)
3184               X3=1.-SQRT(XT2)*EXP(-Y)
3185               X2=2.-X1-X3
3186               NEG=0
3187               IF(K(JL,2).NE.21 .OR. K(JL+1,2).NE.21) THEN
3188         IF((1.-X1)*(1.-X2)*(1.-X3)-X2*SM1*(1.-X1)-X2*SM3*(1.-X3).
3189      &  LE.0..OR.X1.LE.2.*SQRT(SM1)-SM1+SM3.OR.X3.LE.2.*SQRT(SM3)
3190      &  -SM3+SM1) NEG=1
3191         X1=X1+SM1-SM3
3192         X3=X3-SM1+SM3
3193              ENDIF
3194               IF(NEG.EQ.1) GOTO 1
3195      
3196               FG=2.*YMAX*C*(X1**EXP1+X3**EXP3)/D
3197               XT2M=XT2
3198               IF(FG.LT.RANART(NSEED)) GOTO 1
3199      
3200               RETURN
3201               END
3202 C*************************************************************
3203
3204
3205         SUBROUTINE ARORIE(S,X1,X3,JL)
3206 C     
3207         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3208 cc      SAVE /HPARNT/
3209         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3210 cc      SAVE /LUJETSA/
3211       COMMON/RNDF77/NSEED
3212 cc      SAVE /RNDF77/
3213         SAVE   
3214 C     
3215              W=SQRT(S)
3216              X2=2.-X1-X3
3217              E1=.5*X1*W
3218              E3=.5*X3*W
3219              P1=SQRT(E1**2-P(JL,5)**2)
3220         P3=SQRT(E3**2-P(JL+1,5)**2)
3221         CBET=1.
3222         IF(P1.GT.0..AND.P3.GT.0.) CBET=(P(JL,5)**2
3223      &           +P(JL+1,5)**2+2.*E1*E3-S*(1.-X2))/(2.*P1*P3)
3224               IF(ABS(CBET).GT.1.0) CBET=MAX(-1.,MIN(1.,CBET))
3225               BET=ACOS(CBET)
3226      
3227 C.....MINIMIZE PT1-SQUARED PLUS PT3-SQUARED.....
3228               IF(P1.GE.P3) THEN
3229            PSI=.5*ULANGL(P1**2+P3**2*COS(2.*BET),-P3**2*SIN(2.*BET))
3230            PT1=P1*SIN(PSI)
3231            PZ1=P1*COS(PSI)
3232            PT3=P3*SIN(PSI+BET)
3233            PZ3=P3*COS(PSI+BET)
3234               ELSE IF(P3.GT.P1) THEN
3235            PSI=.5*ULANGL(P3**2+P1**2*COS(2.*BET),-P1**2*SIN(2.*BET))
3236            PT1=P1*SIN(BET+PSI)
3237            PZ1=-P1*COS(BET+PSI)
3238            PT3=P3*SIN(PSI)
3239            PZ3=-P3*COS(PSI)
3240               ENDIF
3241      
3242               DEL=2.0*HIPR1(40)*RANART(NSEED)
3243               P(JL,4)=E1
3244               P(JL,1)=PT1*SIN(DEL)
3245               P(JL,2)=-PT1*COS(DEL)
3246               P(JL,3)=PZ1
3247               P(JL+2,4)=E3
3248               P(JL+2,1)=PT3*SIN(DEL)
3249               P(JL+2,2)=-PT3*COS(DEL)
3250               P(JL+2,3)=PZ3
3251               P(JL+1,4)=W-E1-E3
3252               P(JL+1,1)=-P(JL,1)-P(JL+2,1)
3253               P(JL+1,2)=-P(JL,2)-P(JL+2,2)
3254               P(JL+1,3)=-P(JL,3)-P(JL+2,3)
3255               RETURN
3256               END
3257
3258
3259 C
3260 C*******************************************************************
3261 C        make  boost and rotation to entries from IMIN to IMAX
3262 C*******************************************************************
3263         SUBROUTINE ATROBO(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
3264         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3265 cc      SAVE /LUJETSA/
3266         DIMENSION ROT(3,3),PV(3)
3267         DOUBLE PRECISION DP(4),DBEX,DBEY,DBEZ,DGA,DGA2,DBEP,DGABEP
3268         SAVE   
3269         IERROR=0
3270      
3271               IF(IMIN.LE.0 .OR. IMAX.GT.N .OR. IMIN.GT.IMAX) RETURN
3272
3273               IF(THE**2+PHI**2.GT.1E-20) THEN
3274 C...ROTATE (TYPICALLY FROM Z AXIS TO DIRECTION THETA,PHI)
3275            ROT(1,1)=COS(THE)*COS(PHI)
3276            ROT(1,2)=-SIN(PHI)
3277            ROT(1,3)=SIN(THE)*COS(PHI)
3278            ROT(2,1)=COS(THE)*SIN(PHI)
3279            ROT(2,2)=COS(PHI)
3280            ROT(2,3)=SIN(THE)*SIN(PHI)
3281            ROT(3,1)=-SIN(THE)
3282            ROT(3,2)=0.
3283            ROT(3,3)=COS(THE)
3284            DO 120 I=IMIN,IMAX
3285 C**************           IF(MOD(K(I,1)/10000,10).GE.6) GOTO 120
3286               DO 100 J=1,3
3287  100                 PV(J)=P(I,J)
3288                  DO 110 J=1,3
3289  110                    P(I,J)=ROT(J,1)*PV(1)+ROT(J,2)*PV(2)
3290      &                     +ROT(J,3)*PV(3)
3291  120                 CONTINUE
3292         ENDIF
3293      
3294               IF(BEX**2+BEY**2+BEZ**2.GT.1E-20) THEN
3295 C...LORENTZ BOOST (TYPICALLY FROM REST TO MOMENTUM/ENERGY=BETA)
3296                 DBEX=dble(BEX)
3297                 DBEY=dble(BEY)
3298                 DBEZ=dble(BEZ)
3299                 DGA2=1D0-DBEX**2-DBEY**2-DBEZ**2
3300                 IF(DGA2.LE.0D0) THEN
3301                         IERROR=1
3302                         RETURN
3303                 ENDIF
3304                 DGA=1D0/DSQRT(DGA2)
3305                 DO 140 I=IMIN,IMAX
3306 C*************           IF(MOD(K(I,1)/10000,10).GE.6) GOTO 140
3307                    DO 130 J=1,4
3308  130                  DP(J)=dble(P(I,J))
3309                    DBEP=DBEX*DP(1)+DBEY*DP(2)+DBEZ*DP(3)
3310                    DGABEP=DGA*(DGA*DBEP/(1D0+DGA)+DP(4))
3311                    P(I,1)=sngl(DP(1)+DGABEP*DBEX)
3312                    P(I,2)=sngl(DP(2)+DGABEP*DBEY)
3313                    P(I,3)=sngl(DP(3)+DGABEP*DBEZ)
3314                    P(I,4)=sngl(DGA*(DP(4)+DBEP))
3315 140                   CONTINUE
3316               ENDIF
3317      
3318               RETURN
3319               END
3320 C
3321 C
3322 C
3323         SUBROUTINE HIJHRD(JP,JT,JOUT,JFLG,IOPJT0)
3324 C
3325 C        IOPTJET=1, ALL JET WILL FORM SINGLE STRING SYSTEM
3326 C                0, ONLY Q-QBAR JET FORM SINGLE STRING SYSTEM
3327 C*******Perform jets production and fragmentation when JP JT *******
3328 C     scatter. JOUT-> number of hard scatterings precede this one  *
3329 C     for the the same pair(JP,JT). JFLG->a flag to show whether   *
3330 C     jets can be produced (with valence quark=1,gluon=2, q-qbar=3)*
3331 C     or not(0). Information of jets are in  COMMON/ATTJET and     *
3332 C     /MINJET. ABS(NFP(JP,6)) is the total number jets produced by *
3333 C    JP. If NFP(JP,6)<0 JP can not produce jet anymore.                   *
3334 C*******************************************************************
3335         PARAMETER (MAXSTR=150001)
3336         DIMENSION IP(100,2),IPQ(50),IPB(50),IT(100,2),ITQ(50),ITB(50)
3337         COMMON/hjcrdn/YP(3,300),YT(3,300)
3338 cc      SAVE /hjcrdn/
3339         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3340 cc      SAVE /HPARNT/
3341         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3342 cc      SAVE /HIJDAT/
3343         COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3344 cc      SAVE /HSTRNG/
3345         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3346      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3347      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
3348      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3349      &                PJTE(300,500),PJTM(300,500)
3350 cc      SAVE /HJJET1/
3351         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
3352      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
3353      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
3354 cc      SAVE /HJJET2/
3355 c        COMMON/HJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
3356         COMMON/HJJET4/NDR,IADR(MAXSTR,2),KFDR(MAXSTR),PDR(MAXSTR,5)
3357         common/xydr/rtdr(MAXSTR,2)
3358 cc      SAVE /HJJET4/
3359       COMMON/RNDF77/NSEED
3360 cc      SAVE /RNDF77/
3361 C************************************ HIJING common block
3362         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3363 cc      SAVE /LUJETSA/
3364         COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3365 cc      SAVE /LUDAT1A/
3366         COMMON/PYSUBSA/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
3367 cc      SAVE /PYSUBSA/
3368         COMMON/PYPARSA/MSTP(200),PARP(200),MSTI(200),PARI(200)
3369 cc      SAVE /PYPARSA/
3370         COMMON/PYINT1A/MINT(400),VINT(400)
3371 cc      SAVE /PYINT1A/
3372         COMMON/PYINT2A/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
3373 cc      SAVE /PYINT2A/
3374         COMMON/PYINT5A/NGEN(0:200,3),XSEC(0:200,3)
3375 cc      SAVE /PYINT5A/
3376         COMMON/HPINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
3377 cc      SAVE /HPINT/
3378         SAVE   
3379 C*********************************** LU common block
3380         MXJT=500
3381 C                SIZE OF COMMON BLOCK FOR # OF PARTON PER STRING
3382         MXSG=900
3383 C                SIZE OF COMMON BLOCK FOR # OF SINGLE STRINGS
3384         MXSJ=100
3385 C                SIZE OF COMMON BLOCK FOR # OF PARTON PER SINGLE
3386 C                STRING
3387         JFLG=0
3388         IHNT2(11)=JP
3389         IHNT2(12)=JT
3390 C
3391         IOPJET=IOPJT0
3392         IF(IOPJET.EQ.1.AND.(NFP(JP,6).NE.0.OR.NFT(JT,6).NE.0))
3393      &                   IOPJET=0
3394         IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
3395         IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) RETURN
3396 C                ******** JP or JT can not produce jet anymore
3397 C
3398         IF(JOUT.EQ.0) THEN
3399                 EPP=PP(JP,4)+PP(JP,3)
3400                 EPM=PP(JP,4)-PP(JP,3)
3401                 ETP=PT(JT,4)+PT(JT,3)
3402                 ETM=PT(JT,4)-PT(JT,3)
3403                 IF(EPP.LT.0.0) GO TO 1000
3404                 IF(EPM.LT.0.0) GO TO 1000
3405                 IF(ETP.LT.0.0) GO TO 1000
3406                 IF(ETM.LT.0.0) GO TO 1000
3407                 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
3408         ENDIF
3409 C                ********for the first hard scattering of (JP,JT)
3410 C                        have collision only when Ycm(JP)>Ycm(JT)
3411
3412         ECUT1=HIPR1(1)+HIPR1(8)+PP(JP,14)+PP(JP,15)
3413         ECUT2=HIPR1(1)+HIPR1(8)+PT(JT,14)+PT(JT,15)
3414         IF(PP(JP,4).LE.ECUT1) THEN
3415                 NFP(JP,6)=-ABS(NFP(JP,6))
3416                 RETURN
3417         ENDIF
3418         IF(PT(JT,4).LE.ECUT2) THEN
3419                 NFT(JT,6)=-ABS(NFT(JT,6))
3420                 RETURN
3421         ENDIF
3422 C                *********must have enough energy to produce jets
3423
3424         MISS=0
3425         MISP=0
3426         MIST=0
3427 C
3428         IF(NFP(JP,10).EQ.0 .AND. NFT(JT,10).EQ.0) THEN
3429                 MINT(44)=MINT4
3430                 MINT(45)=MINT5
3431                 XSEC(0,1)=ATXS(0)
3432                 XSEC(11,1)=ATXS(11)
3433                 XSEC(12,1)=ATXS(12)
3434                 XSEC(28,1)=ATXS(28)
3435                 DO 120 I=1,20
3436                 COEF(11,I)=ATCO(11,I)
3437                 COEF(12,I)=ATCO(12,I)
3438                 COEF(28,I)=ATCO(28,I)
3439 120                CONTINUE
3440         ELSE
3441                 ISUB11=0
3442                 ISUB12=0
3443                 ISUB28=0
3444                 IF(XSEC(11,1).NE.0) ISUB11=1
3445                 IF(XSEC(12,1).NE.0) ISUB12=1
3446                 IF(XSEC(28,1).NE.0) ISUB28=1                
3447                 MINT(44)=MINT4-ISUB11-ISUB12-ISUB28
3448                 MINT(45)=MINT5-ISUB11-ISUB12-ISUB28
3449                 XSEC(0,1)=ATXS(0)-ATXS(11)-ATXS(12)-ATXS(28)
3450                 XSEC(11,1)=0.0
3451                 XSEC(12,1)=0.0
3452                 XSEC(28,1)=0.0        
3453                 DO 110 I=1,20
3454                 COEF(11,I)=0.0
3455                 COEF(12,I)=0.0
3456                 COEF(28,I)=0.0
3457 110                CONTINUE
3458         ENDIF                
3459 C        ********Scatter the valence quarks only once per NN 
3460 C       collision,
3461 C                afterwards only gluon can have hard scattering.
3462  155        CALL PYTHIAA
3463         JJ=MINT(31)
3464         IF(JJ.NE.1) GO TO 155
3465 C                *********one hard collision at a time
3466         IF(K(7,2).EQ.-K(8,2)) THEN
3467                 QMASS2=(P(7,4)+P(8,4))**2-(P(7,1)+P(8,1))**2
3468      &                        -(P(7,2)+P(8,2))**2-(P(7,3)+P(8,3))**2
3469                 QM=ULMASS(K(7,2))
3470                 IF(QMASS2.LT.(2.0*QM+HIPR1(1))**2) GO TO 155
3471         ENDIF
3472 C                ********q-qbar jets must has minimum mass HIPR1(1)
3473         PXP=PP(JP,1)-P(3,1)
3474         PYP=PP(JP,2)-P(3,2)
3475         PZP=PP(JP,3)-P(3,3)
3476         PEP=PP(JP,4)-P(3,4)
3477         PXT=PT(JT,1)-P(4,1)
3478         PYT=PT(JT,2)-P(4,2)
3479         PZT=PT(JT,3)-P(4,3)
3480         PET=PT(JT,4)-P(4,4)
3481
3482         IF(PEP.LE.ECUT1) THEN
3483                 MISP=MISP+1
3484                 IF(MISP.LT.50) GO TO 155
3485                 NFP(JP,6)=-ABS(NFP(JP,6))
3486                 RETURN
3487         ENDIF
3488         IF(PET.LE.ECUT2) THEN
3489                 MIST=MIST+1
3490                 IF(MIST.LT.50) GO TO 155
3491                 NFT(JT,6)=-ABS(NFT(JT,6))
3492                 RETURN
3493         ENDIF
3494 C                ******** if the remain energy<ECUT the proj or targ
3495 C                         can not produce jet anymore
3496
3497         WP=PEP+PZP+PET+PZT
3498         WM=PEP-PZP+PET-PZT
3499         IF(WP.LT.0.0 .OR. WM.LT.0.0) THEN
3500                 MISS=MISS+1
3501 clin-6/2009 Let user set the limit when selecting high-Pt events 
3502 c     because more attempts may be needed:
3503 c                IF(MISS.LT.50) GO TO 155
3504                 if(pttrig.gt.0) then
3505                    if(MISS.LT.maxmiss) then
3506                 write(6,*) 'Failed to generate minijet Pt>',pttrig,'GeV'
3507                       GO TO 155
3508                    endif
3509                 else
3510                    IF(MISS.LT.50) GO TO 155
3511                 endif
3512
3513                 RETURN
3514         ENDIF
3515 C                ********the total W+, W- must be positive
3516         SW=WP*WM
3517         AMPX=SQRT((ECUT1-HIPR1(8))**2+PXP**2+PYP**2+0.01)
3518         AMTX=SQRT((ECUT2-HIPR1(8))**2+PXT**2+PYT**2+0.01)
3519         SXX=(AMPX+AMTX)**2
3520         IF(SW.LT.SXX.OR.VINT(43).LT.HIPR1(1)) THEN
3521                 MISS=MISS+1
3522 clin-6/2009 ctest on
3523 c                IF(MISS.LT.50) GO TO 155
3524                 IF(MISS.GT.maxmiss) GO TO 155
3525                 RETURN
3526         ENDIF  
3527 C                ********the proj and targ remnants must have at least
3528 C                        a CM energy that can produce two strings
3529 C                        with minimum mass HIPR1(1)(see HIJSFT HIJFRG)
3530 C
3531         HINT1(41)=P(7,1)
3532         HINT1(42)=P(7,2)
3533         HINT1(43)=P(7,3)
3534         HINT1(44)=P(7,4)
3535         HINT1(45)=P(7,5)
3536         HINT1(46)=SQRT(P(7,1)**2+P(7,2)**2)
3537         HINT1(51)=P(8,1)
3538         HINT1(52)=P(8,2)
3539         HINT1(53)=P(8,3)
3540         HINT1(54)=P(8,4)
3541         HINT1(55)=P(8,5)
3542         HINT1(56)=SQRT(P(8,1)**2+P(8,2)**2) 
3543         IHNT2(14)=K(7,2)
3544         IHNT2(15)=K(8,2)
3545 C
3546         PINIRD=(1.0-EXP(-2.0*(VINT(47)-HIDAT(1))))
3547      &                /(1.0+EXP(-2.0*(VINT(47)-HIDAT(1))))
3548         IINIRD=0
3549         IF(RANART(NSEED).LE.PINIRD) IINIRD=1
3550         IF(K(7,2).EQ.-K(8,2)) GO TO 190
3551         IF(K(7,2).EQ.21.AND.K(8,2).EQ.21.AND.IOPJET.EQ.1) GO TO 190
3552 C*******************************************************************
3553 C        gluon  jets are going to be connectd with
3554 C        the final leadng string of quark-aintquark
3555 C*******************************************************************
3556         JFLG=2
3557         JPP=0
3558         LPQ=0
3559         LPB=0
3560         JTT=0
3561         LTQ=0
3562         LTB=0
3563         IS7=0
3564         IS8=0
3565         HINT1(47)=0.0
3566         HINT1(48)=0.0
3567         HINT1(49)=0.0
3568         HINT1(50)=0.0
3569         HINT1(67)=0.0
3570         HINT1(68)=0.0
3571         HINT1(69)=0.0
3572         HINT1(70)=0.0
3573         DO 180 I=9,N
3574            IF(K(I,3).EQ.1 .OR. K(I,3).EQ.2.OR.
3575      &                   ABS(K(I,2)).GT.30) GO TO 180
3576 C************************************************************
3577            IF(K(I,3).EQ.7) THEN
3578               HINT1(47)=HINT1(47)+P(I,1)
3579               HINT1(48)=HINT1(48)+P(I,2)
3580               HINT1(49)=HINT1(49)+P(I,3)
3581               HINT1(50)=HINT1(50)+P(I,4)
3582            ENDIF
3583            IF(K(I,3).EQ.8) THEN
3584               HINT1(67)=HINT1(67)+P(I,1)
3585               HINT1(68)=HINT1(68)+P(I,2)
3586               HINT1(69)=HINT1(69)+P(I,3)
3587               HINT1(70)=HINT1(70)+P(I,4)
3588            ENDIF
3589 C************************modifcation made on Apr 10. 1996*****
3590            IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
3591               NDR=NDR+1
3592               IADR(NDR,1)=JP
3593               IADR(NDR,2)=JT
3594               KFDR(NDR)=K(I,2)
3595               PDR(NDR,1)=P(I,1)
3596               PDR(NDR,2)=P(I,2)
3597               PDR(NDR,3)=P(I,3)
3598               PDR(NDR,4)=P(I,4)
3599               PDR(NDR,5)=P(I,5)
3600               rtdr(NDR,1)=0.5*(YP(1,JP)+YT(1,JT))
3601               rtdr(NDR,2)=0.5*(YP(2,JP)+YT(2,JT))
3602 C************************************************************
3603               GO TO 180
3604 C************************correction made on Oct. 14,1994*****
3605            ENDIF
3606            IF(K(I,3).EQ.7.OR.K(I,3).EQ.3) THEN
3607               IF(K(I,3).EQ.7.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(7,2)
3608      &                     .AND.IS7.EQ.0) THEN
3609                  PP(JP,10)=P(I,1)
3610                  PP(JP,11)=P(I,2)
3611                  PP(JP,12)=P(I,3)
3612                  PZP=PZP+P(I,3)
3613                  PEP=PEP+P(I,4)
3614                  NFP(JP,10)=1
3615                  IS7=1
3616                  GO TO 180
3617               ENDIF
3618               IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
3619      &                               IINIRD.EQ.0)) THEN
3620                  PXP=PXP+P(I,1)
3621                  PYP=PYP+P(I,2)
3622                  PZP=PZP+P(I,3)
3623                  PEP=PEP+P(I,4)
3624                  GO TO 180 
3625               ENDIF
3626               JPP=JPP+1
3627               IP(JPP,1)=I
3628               IP(JPP,2)=0
3629               IF(K(I,2).NE.21) THEN
3630                  IF(K(I,2).GT.0) THEN
3631                     LPQ=LPQ+1
3632                     IPQ(LPQ)=JPP
3633                     IP(JPP,2)=LPQ
3634                  ELSE IF(K(I,2).LT.0) THEN
3635                     LPB=LPB+1
3636                     IPB(LPB)=JPP
3637                     IP(JPP,2)=-LPB
3638                  ENDIF
3639               ENDIF
3640            ELSE IF(K(I,3).EQ.8.OR.K(I,3).EQ.4) THEN
3641               IF(K(I,3).EQ.8.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(8,2)
3642      &                                .AND.IS8.EQ.0) THEN
3643                  PT(JT,10)=P(I,1)
3644                  PT(JT,11)=P(I,2)