label this version hijing ampt in printout so that we can see that standard HIJING...
[u/mrichter/AliRoot.git] / TAmpt / AMPT / hijing1.383_ampt.f
CommitLineData
0119ef9a 1c.................... hijing1.383_ampt.f
2c Version 1.383
3c The variables isng in HIJSFT and JL in ATTRAD were not initialized.
4c The version initialize them. (as found by Fernando Marroquim)
5c
6c
7c
8c Version 1.382
9c Nuclear distribution for deuteron is taken as the Hulthen wave
10c function as provided by Brian Cole (Columbia)
11clin used my own implementation of impact parameter
12clin & proton-neutron distance within a deuteron.
13c
14c
15c Version 1.381
16c
17c The parameters for Wood-Saxon distribution for deuteron are
18c constrained to give the right rms ratius 2.116 fm
19c (R=0.0, D=0.5882)
20c
21c
22c Version 1.38
23c
24c The following common block is added to record the number of elastic
25c (NELT, NELP) and inelastic (NINT, NINP) participants
26c
27c COMMON/HJGLBR/NELT,NINT,NELP,NINP
28c SAVE /HJGLBR/
29c
30c Version 1.37
31c
32c A bug in the quenching subroutine is corrected. When calculating the
33c distance between two wounded nucleons, the displacement of the
34c impact parameter was not inculded. This bug was discovered by
35c Dr. V.Uzhinskii JINR, Dubna, Russia
36c
37c
38C Version 1.36
39c
40c Modification Oct. 8, 1998. In hijing, log(ran(nseed)) occasionally
41c causes overfloat. It is modified to log(max(ran(nseed),1.0e-20)).
42c
43c
44C Nothing important has been changed here. A few 'garbage' has been
45C cleaned up here, like common block HJJET3 for the sea quark strings
46C which were originally created to implement the DPM scheme which
47C later was abadoned in the final version. The lines which operate
48C on these data are also deleted in the program.
49C
50C
51C Version 1.35
52C There are some changes in the program: subroutine HARDJET is now
53C consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
54C for the triggered hard processes. Now that is done altogether
55C with other normal hard processes in modified JETINI. In the new
56C version one calls JETINI every time one calls HIJHRD. In the new
57C version the effect of the isospin of the nucleon on hard processes,
58C especially direct photons is correctly considered.
59C For A+A collisions, one has to initilize pythia
60C separately for each type of collisions, pp, pn,np and nn,
61C or hp and hn for hA collisions. In JETINI we use the following
62C catalogue for different types of collisions:
63C h+h: h+h (itype=1)
64C h+A: h+p (itype=1), h+n (itype=2)
65C A+h: p+h (itype=1), n+h (itype=2)
66C A+A: p+p (itype=1), p+n (itype=2), n+p (itype=3), n+n (itype=4)
67C*****************************************************************
68c
69C
70C Version 1.34
71C Last modification on January 5, 1998. Two mistakes are corrected in
72C function G. A Mistake in the subroutine Parton is also corrected.
73C (These are pointed out by Ysushi Nara).
74C
75C
76C Last modifcation on April 10, 1996. To conduct final
77C state radiation, PYTHIA reorganize the two scattered
78C partons and their final momenta will be a little
79C different. The summed total momenta of the partons
80C from the final state radiation are stored in HINT1(26-29)
81C and HINT1(36-39) which are little different from
82C HINT1(21-24) and HINT1(41-44).
83C
84C Version 1.33
85C
86C Last modfication on September 11, 1995. When HIJING and
87C PYTHIA are initialized, the shadowing is evaluated at
88C b=0 which is the maximum. This will cause overestimate
89C of shadowing for peripheral interactions. To correct this
90C problem, shadowing is set to zero when initializing. Then
91C use these maximum cross section without shadowing as a
92C normalization of the Monte Carlo. This however increase
93C the computing time. IHNT2(16) is used to indicate whether
94C the sturcture function is called for (IHNT2(16)=1) initialization
95C or for (IHNT2(16)=0)normal collisions simulation
96C
97C Last modification on Aagust 28, 1994. Two bugs associate
98C with the impact parameter dependence of the shadowing is
99C corrected.
100C
101C
102c Last modification on October 14, 1994. One bug is corrected
103c in the direct photon production option in subroutine
104C HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
105C Another bug associated with keeping the decay history
106C in the particle information is also corrected.(this problem
107C was reported by Matt Bloomer)
108C
109C
110C Last modification on July 15, 1994. The option to trig on
111C heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1)
112C is added. To do this, set IHPR2(3)=3. For inclusive production,
113C one should reset HIPR1(10)=0.0. One can also trig larger pt
114C QQbar production by giving HIPR1(10) a nonvanishing value.
115C The mass of the heavy quark in the calculation of the cross
116C section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
117C default is the charm mass D=1.5). We also include a separate
118C K-factor for heavy quark and direct photon production by
119C HIPR1(23)(D=2.0).
120C
121C Last modification on May 24, 1994. The option to
122C retain the information of all particles including those
123C who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is
124C added to contain the line number of the parent particle
125C of the current line which is produced via a decay.
126C KATT(I,4) is the status number of the particle: 11=particle
127C which has decayed; 1=finally produced particle.
128C
129C
130C Last modification on May 24, 1994( in HIJSFT when valence quark
131C is quenched, the following error is corrected. 1.2*IHNT2(1) -->
132C 1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
133C
134C
135C Last modification on March 16, 1994 (heavy flavor production
136C processes MSUB(81)=1 MSUB(82)=1 have been switched on,
137C charm production is the default, B-quark option is
138C IHPR2(18), when it is switched on, charm quark is
139C automatically off)
140C
141C
142C Last modification on March 23, 1994 (an error is corrected
143C in the impact parameter dependence of the jet cross section)
144C
145C Last modification Oct. 1993 to comply with non-vax
146C machines' compiler
147C
148C*********************************************
149C LAST MODIFICATION April 5, 1991
150CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B
151C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
152C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
153C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18),
154C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
155C
156C gluon jets can form a single string system
157C
158C initial state radiation is included
159C
160C all QCD subprocesses are included
161c
162c direct particles production is included(currently only direct
163C photon)
164c
165C Effect of high P_T trigger bias on multiple jets distribution
166c
167C******************************************************************
168C HIJING.10 *
169C Heavy Ion Jet INteraction Generator *
170C by *
171C X. N. Wang and M. Gyulassy *
172C Lawrence Berkeley Laboratory *
173C *
174C******************************************************************
175C
176C******************************************************************
177C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of *
178C proj, NFP(K,4) original ID of proj. NFP(K,5)=colli status(0=no,*
179C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
180C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it *
181C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
182C (0=has not been,1=is going to be, -1=has already been scattered *
183C NFP(k,11) total number of interactions this proj has suffered *
184C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant *
185C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and *
186C diquark,PP(K,10)=PT of the hard scattering between the valence *
187C quarks; PP(K,14,15)=the mass of quark,diquark. *
188C******************************************************************
189C
190C****************************************************************
191C
192C SUBROUTINE HIJING
193C
194C****************************************************************
195 SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
196
197cgsfs Added following for consistency with AMPT call
198 double precision BMIN0, BMAX0
199
200cbz1/25/99
201 PARAMETER (MAXPTN=400001)
202clin-4/20/01 PARAMETER (MAXSTR = 1600)
203 PARAMETER (MAXSTR=150001)
204cbz1/25/99end
205clin-4/26/01:
206 PARAMETER (MAXIDL=4001)
207
208cbz1/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
215cbz1/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)
221cc SAVE /HPARNT/
222C
223 COMMON/hjcrdn/YP(3,300),YT(3,300)
224cc SAVE /hjcrdn/
225clin-7/16/03 NINT is a intrinsic fortran function, rename it to NINTHJ
226c COMMON/HJGLBR/NELT,NINT,NELP,NINP
227 COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
228cc SAVE /HJGLBR/
229 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
230cc SAVE /HMAIN1/
231clin-4/26/01
232c COMMON/HMAIN2/KATT(130000,4),PATT(130000,4)
233 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
234cc SAVE /HMAIN2/
235 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
236cc 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)
242cc SAVE /HJJET1/
243clin-4/2008
244c COMMON/HJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
245c & K2SG(900,100),PXSG(900,100),PYSG(900,100),
246c & 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)
250cc SAVE /HJJET2/
251 COMMON/HJJET4/NDR,IADR(MAXSTR,2),KFDR(MAXSTR),PDR(MAXSTR,5)
252clin-4/2008:
253c common/xydr/rtdr(900,2)
254 common/xydr/rtdr(MAXSTR,2)
255cc SAVE /HJJET4/
256 COMMON/RNDF77/NSEED
257cc SAVE /RNDF77/
258C
259 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
260cc SAVE /LUJETSA/
261 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
262cc SAVE /LUDAT1A/
263
264clin-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)
269ccbz11/11/98
270c COMMON /ARPRC/ ITYP(MAXSTR),
271c & GX(MAXSTR), GY(MAXSTR), GZ(MAXSTR), FT(MAXSTR),
272c & PX(MAXSTR), PY(MAXSTR), PZ(MAXSTR), EE(MAXSTR),
273c & XM(MAXSTR)
274cc SAVE /ARPRC/
275ccbz11/11/98end
276
277cbz1/25/99
278 COMMON /PARA1/ MUL
279cc 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)
283cc 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)
287cc SAVE /prec2/
288 COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
289cc SAVE /ilist7/
290 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
291cc SAVE /ilist8/
292 COMMON /SREC1/ NSP, NST, NSI
293cc SAVE /SREC1/
294 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
295cc SAVE /SREC2/
296cbz1/25/99end
297
298clin-2/25/00
299 COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
300 & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
301cc SAVE /frzout/
302clin-4/11/01 soft:
303 common/anim/nevent,isoft,isflag,izpc
304cc SAVE /anim/
305clin-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)
312cc SAVE /SOFT/
313clin-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)
318cc SAVE /NOPREC/
319clin-6/22/01:
320 common /lastt/itimeh,bimp
321cc SAVE /lastt/
322 COMMON /AREVT/ IAEVT, IARUN, MISS
323 common/phidcy/iphidcy,pttrig,ntrig,maxmiss
324cwei DOUBLE PRECISION PATT
325 SAVE
326
327cgsfs WRITE(*,*) "IN Hijing, FRAME=",FRAME
328cgsfs WRITE(*,*) "IN Hijing, BMIN=",BMIN0
329cgsfs 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
337C ********HIPR1(31) is in mb =0.1fm**2
338C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
339C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
340C
341cgsfs WRITE(*,*) "IN Hijing, Modified BMIN=",BMIN
342cgsfs 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)
3485 R=HIRND(1)
349 X=RANART(NSEED)
350 CX=2.0*X-1.0
351 SX=SQRT(1.0-CX*CX)
352C ********choose theta from uniform cos(theta) distr
353 PHI=RANART(NSEED)*2.0*HIPR1(40)
354C ********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
365C ********two neighbors cannot be closer than
366C HIPR1(29)
3678 CONTINUE
36810 CONTINUE
369
370clin-1/27/03 Hulthen wavefn for deuteron borrowed from hijing1.382.f,
371c but modified [divide by 2, & x(p)=-x(n)]:
372c (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)
383c 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
388c 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
40612 CONTINUE
407C
408C******************************
40914 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)
41415 R=HIRND(2)
415 X=RANART(NSEED)
416 CX=2.0*X-1.0
417 SX=SQRT(1.0-CX*CX)
418C ********choose theta from uniform cos(theta) distr
419 PHI=RANART(NSEED)*2.0*HIPR1(40)
420C ********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
431C ********two neighbors cannot be closer than
432C HIPR1(29)
43318 CONTINUE
43420 CONTINUE
435c
436clin-1/27/03 Hulthen wavefn for deuteron borrowed from hijing1.382.f,
437c 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
456c
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
46922 CONTINUE
470
471C********************
47224 MISS=-1
47350 MISS=MISS+1
474
475clin-6/2009 ctest on
476c IF(MISS.GT.50) THEN
477 IF(MISS.GT.maxmiss) THEN
478 WRITE(6,*) 'infinite loop happened in HIJING'
479 STOP
480 ENDIF
481
482clin-4/30/01:
483 itest=0
484
485 NATT=0
486 JATT=0
487 EATT=0.0
488 CALL HIJINI
489 NLOP=0
490C ********Initialize for a new event
49160 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
504C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
505C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
506C BY THE ANGLE PHI FOR EACH COLLISION.******************
507C
508 BB=SQRT(BMIN**2+RANART(NSEED)*(BMAX**2-BMIN**2))
509cbz6/28/99 flow1
d9dcac94 510 PHI=2.0*HIPR1(40)*RANART(NSEED)
511c PHI=0.
0119ef9a 512cbz6/28/99 flow1 end
513 BBX=BB*COS(PHI)
514 BBY=BB*SIN(PHI)
515 HINT1(19)=BB
516 HINT1(20)=PHI
517C
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
523C ********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
551C ********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
55970 CONTINUE
560C ********total number interactions proj and targ has
561C suffered
562
563clin-5/22/01 write impact parameter:
564 bimp=bb
565c 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
573C ********At large impact parameter, there maybe no
574C interaction at all. For NN collision
575C repeat the event until interaction happens
576C
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)
582clin-6/2009 ctest off:
583c write(99,*) IAEVT,NHARD,NCOLT,JPHARD,JTHARD
584 ENDIF
585C
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
592C ********Specifying the location of the hard and
593C minijet if they are enforced by user
594C
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
616C*****************************************************************
617 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
618C ********When IHPR2(8)=0 no jets are produced
619 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
620C ********jets can not be produced for (JP,JT)
621C because not enough energy avaible for
622C 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)
640C
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
651clin-6/2009 ctest off:
652c write(99,*) jp,jt,IHPR2(3),HIPR1(10),njet,
653c 1 ihnt2(9),hint1(21),hint1(22),hint1(23),
654c 2 ihnt2(10),hint1(31),hint1(32),hint1(33)
655c 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
670C ********subtract the trigger jet from total number
671C of jet production to be done since it has
672C 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
684C ********create a hard interaction with specified P_T
685c when IHPR2(3)>0
686 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
687C ********create at least one pair of mini jets
688C when IHPR2(9)=1
689C
690 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
691 & (1.0-EXP(-TTS))) GO TO 160
692C ********this is the probability for no jet production
693110 XR=-ALOG(EXP(-TT)+RANART(NSEED)*(1.0-EXP(-TT)))
694111 NJET=NJET+1
695 XR=XR-ALOG(max(RANART(NSEED),1.0e-20))
696 IF(XR.LT.TT) GO TO 111
697112 NJET=MIN(NJET,IHPR2(8))
698 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
699C ******** Determine number of mini jet production
700C
701 DO 150 ijet=1,NJET
702 CALL JETINI(JP,JT,0)
703 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
704C ********JFLG=1 jets valence quarks, JFLG=2 with
705C gluon jet, JFLG=3 with q-qbar prod for
706C (JP,JT). If JFLG=0 jets can not be produced
707C this time. If JFLG=-1, error occured abandon
708C this event. JOUT is the total hard scat for
709C (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
720C ******** 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
729C
730C ********conduct soft scattering between JP and JT
731 JATT=JATT+JOUT
732200 CONTINUE
733c
734c**************************
735c
736clin-6/2009 write out initial minijet information:
737 call minijet_out(BB)
738 if(pttrig.gt.0.and.ntrig.eq.0) goto 50
739clin-6/2009 write out initial transverse positions of initial nucleons:
740c write(94,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
741 DO 201 JP=1,IHNT2(1)
742clin-6/2009:
743c 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)
751clin-6/2009 target nucleon # has a minus sign for distinction from projectile:
752c 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))
760c
761c*******************************
762
763
764C********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)
770271 CONTINUE
771 DO 272 I=1,IHNT2(3)
772 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
773272 CONTINUE
774 DO 273 ISG=1,NSG
775 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
776273 CONTINUE
777 ENDIF
778
779clin*****4/09/01-soft1, default way of treating strings:
780 if(isoft.eq.1) then
781clin-4/16/01 allow fragmentation:
782 isflag=1
783
784cbz1/25/99
785c.....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)
794cbz1/27/99
795c.....for now only consider gluon cascade
796 IF (KFPJ(I, J) .EQ. 21) THEN
797cbz1/27/99end
798
799 NPAR = NPAR + 1
800 LSTRG0(NPAR) = ISTR
801 LPART0(NPAR) = J
802 ITYP0(NPAR) = KFPJ(I, J)
803cbz6/28/99 flow1
804clin-7/20/01 add dble or sngl to make precisions consistent
805c GX0(NPAR) = YP(1, I)
806 GX0(NPAR) = dble(YP(1, I) + 0.5 * BB)
807cbz6/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))
815c E0(NPAR) = dble(PJPE(I, J))
816 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
817 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
818clin-7/20/01-end
819
820cbz1/27/99
821c.....end gluon selection
822 END IF
823cbz1/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)
829cbz1/27/99
830c.....for now only consider gluon cascade
831 IF (KFTJ(I, J) .EQ. 21) THEN
832cbz1/27/99end
833 NPAR = NPAR + 1
834 LSTRG0(NPAR) = ISTR
835 LPART0(NPAR) = J
836 ITYP0(NPAR) = KFTJ(I, J)
837cbz6/28/99 flow1
838clin-7/20/01 add dble or sngl to make precisions consistent
839c GX0(NPAR) = YT(1, I)
840 GX0(NPAR) = dble(YT(1, I) - 0.5 * BB)
841cbz6/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))
849c 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
853cbz1/27/99
854c.....end gluon selection
855 END IF
856cbz1/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)
862cbz1/27/99
863c.....for now only consider gluon cascade
864 IF (K2SG(I, J) .EQ. 21) THEN
865cbz1/27/99end
866 NPAR = NPAR + 1
867 LSTRG0(NPAR) = ISTR
868 LPART0(NPAR) = J
869 ITYP0(NPAR) = K2SG(I, J)
870clin-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))
881c E0(NPAR) = dble(PESG(I, J))
882 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
883 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
884cbz1/27/99
885c.....end gluon selection
886 END IF
887cbz1/27/99end
888 1011 CONTINUE
889 1012 CONTINUE
890 MUL = NPAR
891
892cbz2/4/99
893 CALL HJANA1
894cbz2/4/99end
895
896clin-6/2009:
897 if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
898c.....call ZPC for parton cascade
899 CALL ZPCMN
900
901c write out parton and wounded nucleon information to ana/zpc1.mom:
902clin-6/2009:
903c WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
904c WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
905 DO 1013 I = 1, MUL
906cc WRITE (14, 411) PX5(I), PY5(I), PZ5(I), ITYP5(I),
907c & XMASS5(I), E5(I)
908 if(dmax1(abs(GX5(I)),abs(GY5(I)),abs(GZ5(I)),abs(FT5(I)))
909 1 .lt.9999) then
910c write(14,210) ITYP5(I), PX5(I), PY5(I), PZ5(I), XMASS5(I),
911c 1 GX5(I), GY5(I), GZ5(I), FT5(I)
912 else
913c change format for large numbers:
914c write(14,211) ITYP5(I), PX5(I), PY5(I), PZ5(I), XMASS5(I),
915c 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
923clin-4/09/01:
924 itest=itest+1
925c 411 FORMAT(1X, 3F10.3, I6, 2F10.3)
926cbz3/19/99 end
927
928clin-5/2009 ctest off:
929c call frztm(1,1)
930
931c.....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)
937clin-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
963cbz1/25/99end
964
965cbz2/4/99
966 CALL HJANA2
967cbz2/4/99end
968
969clin*****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)
973clin-4/27/01:
974 NSI = NSG
975 NPAR=0
976 ISTR=0
977C
978clin No fragmentation to hadrons, only on parton level,
979c and transfer minijet and string data from HIJING to ZPC:
980 MSTJ(1)=0
981clin-4/12/01 forbid soft radiation before ZPC to avoid small-mass strings,
982c 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
990c change: do gluon kink only once: either here or in fragmentation.
991 CALL HIJFRG(jjtp,NTP,IERROR)
992c call lulist(1)
993 if(NTP.eq.1) then
994c 354 continue
995 NPJ(jjtp)=MAX0(N-2,0)
996
997clin-4/12/01: NPJ(jjtp)=MAX0(ipartn-2,0)
998 else
999c 355 continue
1000 NTJ(jjtp)=MAX0(N-2,0)
1001clin-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
1011clin-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))
1016c 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
1020clin-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)
1083c call lulist(2)
1084c
1085 NJSG(ISG)=N
1086c
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))
1102c 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
1110cbz2/4/99
1111 CALL HJANA1
1112cbz2/4/99end
1113clin-6/2009:
1114 if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
1115c.....call ZPC for parton cascade
1116 CALL ZPCMN
1117cbz3/19/99
1118clin-6/2009:
1119c 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
1124c WRITE (14, 311) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1125c & 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
1129c 311 FORMAT(1X, 3F10.4, I6, 2F10.4)
1130 312 FORMAT(1X, 3F10.3, I6, 2F10.3,1X,I6,1X,I3)
1131cbz3/19/99 end
1132
1133clin-5/2009 ctest off:
1134c call frztm(1,1)
1135
1136clin-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
1145clin-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)
1151c nucleons without interactions:
1152 if(IITYP.eq.2112.or.IITYP.eq.2212) then
1153clin-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))
1159c 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)
1171c 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)
1186c 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
1249cbz1/25/99end
1250
1251clin-4/09/01 turn on fragmentation with soft radiation
1252c and jet order reversal to form hadrons after ZPC:
1253 MSTJ(1)=1
1254 IHPR2(1)=1
1255 isflag=1
1256clin-4/13/01 allow small mass strings (D=1.5GeV):
1257 HIPR1(1)=0.94
1258
1259cbz2/4/99
1260 CALL HJANA2
1261cbz2/4/99end
1262
1263clin-4/19/01-soft3, fragment strings, then convert hadrons to partons
1264c and input to ZPC:
1265 elseif(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
1266clin-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)
1272C
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
1287C ******** boost back to lab frame(if it was in)
1288C
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)
1300c from Yasushi, to avoid violation of array limits:
1301c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1302clin-4/2008 to avoid out-of-bound error in K():
1303c IF(K(I,3).EQ.0 .OR.
1304c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1305c 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
1310clin-4/2008-end
1311 ELSE
1312 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1313 ENDIF
1314
1315C ****** 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)
1333cbz11/11/98end
1334
1335 560 CONTINUE
1336C ********Fragment the q-qbar jets systems *****
1337C
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)
1343C
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
1357C ******** boost back to lab frame(if it was in)
1358C
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)
1372c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1373clin-4/2008
1374c IF(K(I,3).EQ.0 .OR.
1375c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1376c 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
1381clin-4/2008-end
1382 ELSE
1383 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1384 ENDIF
1385
1386C ****** 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)
1407cbz11/11/98end
1408
1409 590 CONTINUE
1410 600 CONTINUE
1411C ********Fragment the q-qq related string systems
1412 ENDIF
1413clin-4/2008 check for zero NDR value:
1414 if(NDR.ge.1) then
1415c
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)
1426clin-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
1438clin-4/2008:
1439 endif
1440clin-6/2009 ctest on:
1441 call embedHighPt
1442c
1443 CALL HJANA1
1444
1445clin-4/19/01 convert hadrons to partons for ZPC (with GX0 given):
1446 call htop
1447
1448clin-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
1453clin-7/03/01-end
1454
1455clin-6/2009:
1456 if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
1457c.....call ZPC for parton cascade
1458 CALL ZPCMN
1459clin-6/2009:
1460c 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
1465c WRITE (14, 511) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1466c & 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
1471c 511 FORMAT(1X, 3F10.4, I6, 2F10.4)
1472 512 FORMAT(I6,4(1X,F10.3),1X,I6,1X,I3,1X,F10.3)
1473c 513 FORMAT(1X, 4F10.4)
1474
1475clin-5/2009 ctest off:
1476c call frztm(1,1)
1477
1478clin save data after ZPC for fragmentation purpose:
1479c.....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)
1504clin-7/20/01 E5(I) does no include the finite parton mass XMASS5(I),
1505c so define it anew:
1506c PESGS(NSTRG, NPART) = E5(I)
1507c if(abs(PZ5(i)/E5(i)).gt.0.9999999d0)
1508c 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)
1511c if(abs(PZ5(i)/E5(i)).gt.0.9999999d0)
1512c 1 write(91,*) 'b: new E5(I)=',E5(i)
1513clin-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
1521clin-4/19/01-end
1522
1523 endif
1524clin-4/09/01-end
1525
1526C
1527C**************fragment all the string systems in the following*****
1528C
1529C********nsbst is where particle information starts
1530C********nsbstR+1 is the number of strings in fragmentation
1531C********the number of strings before a line is stored in K(I,4)
1532C********IDSTR is id number of the string system (91,92 or 93)
1533C
1534clin-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
1560clin-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
1568c 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
1575C ********Check errors
1576C
1577 nsbst=1
1578 IDSTR=92
1579 IF(IHPR2(21).EQ.0) THEN
1580 CALL LUEDIT(2)
1581 ELSE
1582351 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
1587C
1588 IF(FRAME.EQ.'LAB') THEN
1589 CALL HBOOST
1590 ENDIF
1591C ******** boost back to lab frame(if it was in)
1592C
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)
1604c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1605clin-4/2008:
1606c IF(K(I,3).EQ.0 .OR.
1607c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1608c 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
1613clin-4/2008-end
1614 ELSE
1615 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1616 ENDIF
1617
1618C ****** 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
1625cbz11/11/98
1626cbz1/25/99
1627c GXAR(NATT) = 0.5 * (YP(1, IASG(ISG, 1)) +
1628c & YT(1, IASG(ISG, 2)))
1629c GYAR(NATT) = 0.5 * (YP(2, IASG(ISG, 1)) +
1630c & 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))
1636cbz1/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)
1643cbz11/11/98end
1644
1645360 CONTINUE
1646C ********Fragment the q-qbar jets systems *****
1647C
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
1657c call lulist(2)
1658 WRITE(6,*) 'error occured P&T, repeat the event'
1659 WRITE(6,*) NTP,jjtp
1660clin-6/2009 when this happens, the event will be repeated,
1661c and another record for the same event number will be written into
1662c zpc.dat, zpc.res, minijet-initial-beforePropagation.dat,
1663c parton-initial-afterPropagation.dat, parton-after-coalescence.dat,
1664c and parton-collisionsHistory.dat.
1665 ENDIF
1666 GO TO 50
1667 ENDIF
1668C ********check errors
1669C
1670 nsbst=1
1671 IDSTR=92
1672 IF(IHPR2(21).EQ.0) THEN
1673 CALL LUEDIT(2)
1674 ELSE
1675381 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
1683C ******** boost back to lab frame(if it was in)
1684C
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)
1698c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1699clin-4/2008:
1700c IF(K(I,3).EQ.0 .OR.
1701c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1702c 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
1707clin-4/2008-end
1708 ELSE
1709 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1710 ENDIF
1711C ****** 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)
1717cbz11/11/98
1718cbz1/25/99
1719c IF (NTP .EQ. 1) THEN
1720c GXAR(NATT) = YP(1, jjtp)
1721c ELSE
1722c GXAR(NATT) = YT(1, jjtp)
1723c END IF
1724c IF (NTP .EQ. 1) THEN
1725c GYAR(NATT) = YP(2, jjtp)
1726c ELSE
1727c GYAR(NATT) = YT(2, jjtp)
1728c 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))
1738cbz1/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)
1745cbz11/11/98end
1746
1747390 CONTINUE
1748400 CONTINUE
1749C ********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)
1762clin-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
1775C ********store the direct-produced particles
1776C
1777
1778clin-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'
1786c call lulist(1)
1787 write(6,*) 'violated:EATT,NATT,B=',EATT,NATT,bimp
1788 GO TO 50
1789 ENDIF
1790c write(6,*) 'satisfied:EATT,NATT,B=',EATT,NATT,bimp
1791c write(6,*) ' '
1792
1793 RETURN
1794 END
1795C
1796C
1797C
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)
1802cc SAVE /HSTRNG/
1803 COMMON/hjcrdn/YP(3,300),YT(3,300)
1804cc SAVE /hjcrdn/
1805 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1806cc SAVE /HPARNT/
1807 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
1808cc SAVE /HIJDAT/
1809 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1810cc 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
1821C
1822 HINT1(8)=MAX(ULMASS(2112),ULMASS(2212))
1823 HINT1(9)=HINT1(8)
1824C
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
1872C...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
1873 IF(IHPR2(12).GT.0) THEN
1874 CALL LUGIVE('MDCY(C221,1)=0')
1875clin-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')
1880clin-1/04/01 no K0 and K0bar decays so K0L and K0S do not appear,
1881c 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')
1884clin-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')
1893clin-11/07/00-end
1894cbz12/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')
1900cbz12/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
1925clin parj(41) and (42) are a, b parameters in Lund, read from input.ampt:
1926c PARJ(41)=HIPR1(3)
1927c PARJ(42)=HIPR1(4)
1928c PARJ(41)=2.2
1929c PARJ(42)=0.5
1930
1931clin 2 popcorn parameters read from input.ampt:
1932c IHPR2(11) = 3
1933c PARJ(5) = 0.5
1934 MSTJ(12)=IHPR2(11)
1935
1936clin parj(21) gives the mean gaussian width for hadron Pt:
1937 PARJ(21)=HIPR1(2)
1938clin 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)
1942clin-10/31/00 update when string tension is changed:
1943 HIPR1(2)=PARJ(21)
1944
1945C ******** 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
1973C ********define Lorentz transform to lab frame
1974c
1975C ********calculate the cross sections involved with
1976C nucleon collisions.
1977 IF(IHNT2(1).GT.1) THEN
1978 CALL HIJWDS(IHNT2(1),1,RMAX)
1979 HIPR1(34)=RMAX
1980C ********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
1985C ********set up Wood-Sax distr for targ.
1986 ENDIF
1987C
1988C
1989 I=0
199020 I=I+1
1991 IF(I.EQ.10) GO TO 30
1992 IF(HIDAT0(10,I).LE.HINT1(1)) GO TO 20
199330 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))
199740 CONTINUE
1998 HIPR1(31)=HIDAT(5)
1999 HIPR1(30)=2.0*HIDAT(5)
2000C
2001C
2002 CALL HIJCRS
2003C
2004 IF(IHPR2(5).NE.0) THEN
2005 CALL HIFUN(3,0.0,36.0,FNKICK)
2006C ********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)
2012C ********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)
2017100 FORMAT(
2018 & 10X,'**************************************************'/
2019 & 10X,'*',48X,'*'/
7a129c8c 2020 & 10X,'* HIJING for AMPT initialized at *'/
0119ef9a 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
2028C
2029C
2030C
2031 FUNCTION FNKICK(X)
2032 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2033cc 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
2039C
2040C
2041 FUNCTION FNKC2(X)
2042 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2043cc SAVE /HPARNT/
2044 SAVE
2045 FNKC2=X*EXP(-2.0*X/HIPR1(42))
2046 RETURN
2047 END
2048C
2049C
2050C
2051 FUNCTION FNSTRU(X)
2052 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2053cc 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
2059C
2060C
2061C
2062 FUNCTION FNSTRM(X)
2063 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2064cc 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
2070C
2071C
2072 FUNCTION FNSTRS(X)
2073 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2074cc 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
2080C
2081C
2082C
2083C
2084 SUBROUTINE HBOOST
2085 IMPLICIT DOUBLE PRECISION(D)
2086 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
2087cc SAVE /LUJETSA/
2088 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2089cc SAVE /LUDAT1A/
2090 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2091cc 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
2098C ********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))
2113100 CONTINUE
2114 RETURN
2115 END
2116C
2117C
2118C
2119C
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)
2124cc SAVE /hjcrdn/
2125 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2126cc SAVE /HPARNT/
2127C
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)
2133cc 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)
2137cc SAVE /HJJET2/
2138 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2139cc SAVE /HSTRNG/
2140 COMMON/RNDF77/NSEED
2141cc SAVE /RNDF77/
2142 SAVE
2143C
2144c Uzhi:
2145 BB=HINT1(19)
2146 PHI=HINT1(20)
2147 BBX=BB*COS(PHI)
2148 BBY=BB*SIN(PHI)
2149c
2150 IF(NTP.EQ.2) GO TO 400
2151 IF(NTP.EQ.3) GO TO 2000
2152C*******************************************************
2153C Jet interaction for proj jet in the direction PHIP
2154C******************************************************
2155C
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))
2165C******* 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)
2173c 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
2182C******* 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
2193C****** 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)
2201c 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
2210C******* 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
2242C ********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
2259C ********reshuffle the energy when jet has mass
2260 R0=RDP(MP)
2261 DP1=DP*V1
2262 DP2=DP*V2
2263 DP3=DP*V3
2264C ********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
2295C ********reshuffle the energy when jet has mass
2296
2297 R0=RDT(MT)
2298 DP1=DP*V1
2299 DP2=DP*V2
2300 DP3=DP*V3
2301C ********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
2322C*******************************************************
2323C Jet interaction for target jet in the direction PHIT
2324C******************************************************
2325C
2326C******* 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)
2343c 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
2352C******* 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
2363C****** 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)
2371c 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
2380C******* 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
2407620 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
2412C ********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
2417C
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
2429C ********reshuffle the energy when jet has mass
2430C
2431 R0=RDP(MP)
2432 DP1=DP*V1
2433 DP2=DP*V2
2434 DP3=DP*V3
2435C ********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
2446640 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
2466C ********reshuffle the energy when jet has mass
2467
2468 R0=RDT(MT)
2469 DP1=DP*V1
2470 DP2=DP*V2
2471 DP3=DP*V3
2472C ********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
2481660 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
2489690 CONTINUE
2490 RETURN
2491C********************************************************
2492C Q-QBAR jet interaction
2493C********************************************************
24942000 ISG=JPJT
2495 IF(IASG(ISG,3).NE.1) RETURN
2496C
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)
2515c 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
2524C******* 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
2535C****** 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)
2543c 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
2552C******* 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
2588C
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
2601C ********reshuffle the energy when jet has mass
2602C
2603 R0=RDP(MP)
2604 DP1=DP*V1
2605 DP2=DP*V2
2606 DP3=DP*V3
2607C ********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
2639C ********reshuffle the energy when jet has mass
2640
2641 R0=RDT(MT)
2642 DP1=DP*V1
2643 DP2=DP*V2
2644 DP3=DP*V3
2645C ********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
2666C
2667C
2668C
2669C
2670 SUBROUTINE HIJFRG(JTP,NTP,IERROR)
2671C NTP=1, fragment proj string, NTP=2, targ string,
2672C NTP=3, independent
2673C strings from jets. JTP is the line number of the string
2674C*******Fragment all leadng strings of proj and targ**************
2675C IHNT2(1)=atomic #, IHNT2(2)=proton #(=-1 if anti-proton) *
2676C******************************************************************
2677 PARAMETER (MAXSTR=150001)
2678 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2679cc SAVE /HPARNT/
2680 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
2681cc SAVE /HIJDAT/
2682 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2683cc 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)
2689cc 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)
2693cc SAVE /HJJET2/
2694C
2695 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
2696cc SAVE /LUJETSA/
2697 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2698cc SAVE /LUDAT1A/
2699 COMMON/RNDF77/NSEED
2700cc SAVE /RNDF77/
2701clin-4/11/01 soft:
2702 common/anim/nevent,isoft,isflag,izpc
2703cc SAVE /anim/
2704 SAVE
2705
2706cbz3/12/99
2707c.....set up fragmentation function according to the number of collisions
2708c.....a wounded nucleon has suffered
2709c IF (NTP .EQ. 1) THEN
2710c NCOLL = NFP(JTP, 11)
2711c ELSE IF (NTP .EQ. 2) THEN
2712c NCOLL = NFT(JTP, 11)
2713c ELSE IF (NTP .EQ. 3) THEN
2714c NCOLL = (NFP(IASG(JTP,1), 11) + NFT(IASG(JTP,2), 11)) / 2
2715c END IF
2716c IF (NCOLL .LE. 1) THEN
2717c PARJ(5) = 0.5
2718c ELSE IF (NCOLL .EQ. 2) THEN
2719c PARJ(5) = 0.75
2720c ELSE IF (NCOLL .EQ. 3) THEN
2721c PARJ(5) = 1.17
2722c ELSE IF (NCOLL .EQ. 4) THEN
2723c PARJ(5) = 2.0
2724c ELSE IF (NCOLL .EQ. 5) THEN
2725c PARJ(5) = 4.5
2726c ELSE IF (NCOLL .GE. 6) THEN
2727c PARJ(5) = 49.5
2728c END IF
2729c PARJ(5) = 0.5
2730cbz3/12/99 end
2731
2732 IERROR=0
2733 CALL LUEDIT(0)
2734 N=0
2735C ********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
2749C IF(IHPR2(1).GT.0) CALL ATTRAD(IERROR)
2750c IF(IERROR.NE.0) RETURN
2751C CALL LULIST(1)
2752 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2753 RETURN
2754 ENDIF
2755C
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
2780C ********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
2792200 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
2814C ********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)
2825300 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'
2829clin:
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)
2838C *******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
2854C*****
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
2875C*******************attach produced jets to the leadng partons****
2876 IF(NTP.EQ.1.AND.NPJ(JTP).NE.0) THEN
2877 JETOT=NPJ(JTP)
2878C IF(NPJ(JTP).GE.2) CALL HIJSRT(JTP,1)
2879C ********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)
2889520 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
2895542 CONTINUE
2896 I0=I
2897clin-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
2900C ********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)
2911540 CONTINUE
2912 N=N+NPJ(JTP)
2913 ELSE IF(NTP.EQ.2.AND.NTJ(JTP).NE.0) THEN
2914 JETOT=NTJ(JTP)
2915c IF(NTJ(JTP).GE.2) CALL HIJSRT(JTP,2)
2916C ********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)
2926560 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
2931582 CONTINUE
2932 I0=I
2933clin-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
2936C ********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)
2947580 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
2976C ******** conduct soft radiations
2977C****************************
2978C
2979C
2980clin-4/11/01 soft:
2981c CALL LUEXEC
2982 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2983
2984 RETURN
2985
29861000 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)
29911100 CONTINUE
2992C ********proj remain as a nucleon or delta
2993clin-4/11/01 soft:
2994c CALL LUEXEC
2995 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2996
2997C call lulist(1)
2998 RETURN
2999C
30001200 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)
30051300 CONTINUE
3006C ********targ remain as a nucleon or delta
3007clin-4/11/01 soft:
3008c CALL LUEXEC
3009 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
3010
3011C call lulist(1)
3012 RETURN
3013 END
3014C
3015C
3016C
3017C
3018C****************************************************************
3019C conduct soft radiation according to dipole approxiamtion
3020C****************************************************************
3021 SUBROUTINE ATTRAD(IERROR)
3022C
3023 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3024cc SAVE /HPARNT/
3025 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3026cc SAVE /HIJDAT/
3027 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3028cc SAVE /LUJETSA/
3029 COMMON/RNDF77/NSEED
3030cc SAVE /RNDF77/
3031 SAVE
3032 IERROR=0
3033
3034C.....S INVARIANT MASS-SQUARED BETWEEN PARTONS I AND I+1......
3035C.....SM IS THE LARGEST MASS-SQUARED....
3036
303740 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
305630 CONTINUE
3057 S=(SM+1.5*(P(JL,5)+P(JL+1,5)))**2
3058 IF(SM.LT.HIPR1(5)) GOTO 2
3059
3060C.....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
3066150 P(J+1,M)=P(J,M)
3067160 CONTINUE
3068190 N=N+1
3069
3070C.....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
3082C.....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
3090C.....CREATE ONE GLUON AND ORIENTATE.....
3091
30921 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
3104C.....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
3109C.....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.
3122C----THETA FUNCTION DAMPING OF THE EMITTED GLUONS. FOR HADRON-HADRON.
3123C----R0=VFR(2)
3124C IF(VFR(2).GT.0.) THEN
3125C PTG=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
3126C PTGMAX=WSTRI/2.
3127C DOPT=SQRT((4.*PAR(71)*VFR(2))/WSTRI)
3128C PTOPT=(DOPT*WSTRI)/(2.*VFR(2))
3129C IF(PTG.GT.PTOPT) IORDER=IORDER-1
3130C IF(PTG.GT.PTOPT) GOTO 1
3131C ENDIF
3132C-----
3133 IF(SM.GE.HIPR1(5)) GOTO 40
3134
31352 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)
3149C
3150 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3151cc SAVE /HPARNT/
3152 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3153cc SAVE /LUJETSA/
3154 COMMON/RNDF77/NSEED
3155cc SAVE /RNDF77/
3156 SAVE
3157C
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
31721 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
3202C*************************************************************
3203
3204
3205 SUBROUTINE ARORIE(S,X1,X3,JL)
3206C
3207 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3208cc SAVE /HPARNT/
3209 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3210cc SAVE /LUJETSA/
3211 COMMON/RNDF77/NSEED
3212cc SAVE /RNDF77/
3213 SAVE
3214C
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
3227C.....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
3259C
3260C*******************************************************************
3261C make boost and rotation to entries from IMIN to IMAX
3262C*******************************************************************
3263 SUBROUTINE ATROBO(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
3264 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3265cc 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
3274C...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
3285C************** 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
3295C...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
3306C************* 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))
3315140 CONTINUE
3316 ENDIF
3317
3318 RETURN
3319 END
3320C
3321C
3322C
3323 SUBROUTINE HIJHRD(JP,JT,JOUT,JFLG,IOPJT0)
3324C
3325C IOPTJET=1, ALL JET WILL FORM SINGLE STRING SYSTEM
3326C 0, ONLY Q-QBAR JET FORM SINGLE STRING SYSTEM
3327C*******Perform jets production and fragmentation when JP JT *******
3328C scatter. JOUT-> number of hard scatterings precede this one *
3329C for the the same pair(JP,JT). JFLG->a flag to show whether *
3330C jets can be produced (with valence quark=1,gluon=2, q-qbar=3)*
3331C or not(0). Information of jets are in COMMON/ATTJET and *
3332C /MINJET. ABS(NFP(JP,6)) is the total number jets produced by *
3333C JP. If NFP(JP,6)<0 JP can not produce jet anymore. *
3334C*******************************************************************
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)
3338cc SAVE /hjcrdn/
3339 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3340cc SAVE /HPARNT/
3341 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3342cc SAVE /HIJDAT/
3343 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3344cc 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)
3350cc 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)
3354cc SAVE /HJJET2/
3355c 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)
3358cc SAVE /HJJET4/
3359 COMMON/RNDF77/NSEED
3360cc SAVE /RNDF77/
3361C************************************ HIJING common block
3362 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3363cc SAVE /LUJETSA/
3364 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3365cc SAVE /LUDAT1A/
3366 COMMON/PYSUBSA/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
3367cc SAVE /PYSUBSA/
3368 COMMON/PYPARSA/MSTP(200),PARP(200),MSTI(200),PARI(200)
3369cc SAVE /PYPARSA/
3370 COMMON/PYINT1A/MINT(400),VINT(400)
3371cc SAVE /PYINT1A/
3372 COMMON/PYINT2A/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
3373cc SAVE /PYINT2A/
3374 COMMON/PYINT5A/NGEN(0:200,3),XSEC(0:200,3)
3375cc SAVE /PYINT5A/
3376 COMMON/HPINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
3377cc SAVE /HPINT/
3378 SAVE
3379C*********************************** LU common block
3380 MXJT=500
3381C SIZE OF COMMON BLOCK FOR # OF PARTON PER STRING
3382 MXSG=900
3383C SIZE OF COMMON BLOCK FOR # OF SINGLE STRINGS
3384 MXSJ=100
3385C SIZE OF COMMON BLOCK FOR # OF PARTON PER SINGLE
3386C STRING
3387 JFLG=0
3388 IHNT2(11)=JP
3389 IHNT2(12)=JT
3390C
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
3396C ******** JP or JT can not produce jet anymore
3397C
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
3409C ********for the first hard scattering of (JP,JT)
3410C 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
3422C *********must have enough energy to produce jets
3423
3424 MISS=0
3425 MISP=0
3426 MIST=0
3427C
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)
3439120 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
3457110 CONTINUE
3458 ENDIF
3459C ********Scatter the valence quarks only once per NN
3460C collision,
3461C afterwards only gluon can have hard scattering.
440e3d40 3462 155 CALL PYTHIAA
0119ef9a 3463 JJ=MINT(31)
3464 IF(JJ.NE.1) GO TO 155
3465C *********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
3472C ********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
3494C ******** if the remain energy<ECUT the proj or targ
3495C 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
3501clin-6/2009 Let user set the limit when selecting high-Pt events
3502c because more attempts may be needed:
3503c 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
3515C ********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
3522clin-6/2009 ctest on
3523c IF(MISS.LT.50) GO TO 155
3524 IF(MISS.GT.maxmiss) GO TO 155
3525 RETURN
3526 ENDIF
3527C ********the proj and targ remnants must have at least
3528C a CM energy that can produce two strings
3529C with minimum mass HIPR1(1)(see HIJSFT HIJFRG)
3530C
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)
3545C
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
3552C*******************************************************************
3553C gluon jets are going to be connectd with
3554C the final leadng string of quark-aintquark
3555C*******************************************************************
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
3576C************************************************************
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
3589C************************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))
3602C************************************************************
3603 GO TO 180
3604C************************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)
3645 PT(JT,12)=P(I,3)
3646 PZT=PZT+P(I,3)
3647 PET=PET+P(I,4)
3648 NFT(JT,10)=1
3649 IS8=1
3650 GO TO 180
3651 ENDIF
3652 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
3653 & IINIRD.EQ.0)) THEN
3654 PXT=PXT+P(I,1)
3655 PYT=PYT+P(I,2)
3656 PZT=PZT+P(I,3)
3657 PET=PET+P(I,4)
3658 GO TO 180
3659 ENDIF
3660 JTT=JTT+1
3661 IT(JTT,1)=I
3662 IT(JTT,2)=0
3663 IF(K(I,2).NE.21) THEN
3664 IF(K(I,2).GT.0) THEN
3665 LTQ=LTQ+1
3666 ITQ(LTQ)=JTT
3667 IT(JTT,2)=LTQ
3668 ELSE IF(K(I,2).LT.0) THEN
3669 LTB=LTB+1
3670 ITB(LTB)=JTT
3671 IT(JTT,2)=-LTB
3672 ENDIF
3673 ENDIF
3674 ENDIF
3675 180 CONTINUE
3676c
3677c
3678 IF(LPQ.NE.LPB .OR. LTQ.NE.LTB) THEN
3679 MISS=MISS+1
3680clin-6/2009 ctest on
3681c IF(MISS.LE.50) GO TO 155
3682 IF(MISS.LE.maxmiss) GO TO 155
3683 WRITE(6,*) ' Q -QBAR NOT MATCHED IN HIJHRD'
3684 JFLG=0
3685 RETURN
3686 ENDIF
3687C****The following will rearrange the partons so that a quark is***
3688C****allways followed by an anti-quark ****************************
3689
3690 J=0
3691181 J=J+1
3692 IF(J.GT.JPP) GO TO 182
3693 IF(IP(J,2).EQ.0) THEN
3694 GO TO 181
3695 ELSE IF(IP(J,2).NE.0) THEN
3696 LP=ABS(IP(J,2))
3697 IP1=IP(J,1)
3698 IP2=IP(J,2)
3699 IP(J,1)=IP(IPQ(LP),1)
3700 IP(J,2)=IP(IPQ(LP),2)
3701 IP(IPQ(LP),1)=IP1
3702 IP(IPQ(LP),2)=IP2
3703 IF(IP2.GT.0) THEN
3704 IPQ(IP2)=IPQ(LP)
3705 ELSE IF(IP2.LT.0) THEN