ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hijing.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hijing.f
1 c Version 1.383
2 c The variables I_SNG in HIJSFT and JL in ATTRAD were not initialized.
3 c The version initialize them. (as found by Fernando Marroquim)
4 c
5 c
6 c
7 c Version 1.382
8 c Nuclear distribution for deuteron is taken as the Hulthen wave
9 c function as provided by Brian Cole (Columbia)
10 c
11 c
12 c Version 1.381
13 c
14 c The parameters for Wood-Saxon distribution for deuteron are
15 c constrained to give the right rms ratius 2.116 fm
16 c (R=0.0, D=0.5882)
17 c
18 c
19 c Version 1.38
20 c
21 c The following common block is added to record the number of elastic
22 c (NELT, NELP) and inelastic (NINT, NINP) participants
23 c
24 c COMMON/HIJGLBR/NELT,NINT,NELP,NINP
25 c SAVE /HIJGLBR/
26 c
27 c Version 1.37
28 c
29 c A bug in the quenching subroutine is corrected. When calculating the
30 c distance between two wounded nucleons, the displacement of the
31 c impact parameter was not inculded. This bug was discovered by
32 c Dr. V.Uzhinskii JINR, Dubna, Russia
33 c
34 c
35 C Version 1.36
36 c
37 c Modification Oct. 8, 1998. In hijing, log(ran(nseed)) occasionally
38 c causes overfloat. It is modified to log(max(ran(nseed),1.0e-20)).
39 c
40 c
41 C Nothing important has been changed here. A few 'garbage' has been
42 C cleaned up here, like common block HIJJET3 for the sea quark strings
43 C which were originally created to implement the DPM scheme which
44 C later was abadoned in the final version. The lines which operate
45 C on these data are also deleted in the program.
46 C
47 C
48 C Version 1.35
49 C There are some changes in the program: subroutine HARDJET is now
50 C consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
51 C for the triggered hard processes. Now that is done altogether
52 C with other normal hard processes in modified JETINI. In the new
53 C version one calls JETINI every time one calls HIJHRD. In the new
54 C version the effect of the isospin of the nucleon on hard processes,
55 C especially direct photons is correctly considered.
56 C For A+A collisions, one has to initilize pythia
57 C separately for each type of collisions, pp, pn,np and nn,
58 C or hp and hn for hA collisions. In JETINI we use the following
59 C catalogue for different types of collisions:
60 C h+h: h+h (I_TYPE=1)
61 C h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
62 C A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
63 C A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
64 C*****************************************************************
65 c
66 C
67 C Version 1.34
68 C Last modification on January 5, 1998. Two mistakes are corrected in
69 C function G. A Mistake in the subroutine Parton is also corrected.
70 C (These are pointed out by Ysushi Nara).
71 C
72 C
73 C Last modifcation on April 10, 1996. To conduct final
74 C state radiation, PYTHIA reorganize the two scattered
75 C partons and their final momenta will be a little
76 C different. The summed total momenta of the partons
77 C from the final state radiation are stored in HINT1(26-29)
78 C and HINT1(36-39) which are little different from
79 C HINT1(21-24) and HINT1(41-44).
80 C
81 C Version 1.33
82 C
83 C Last modfication on September 11, 1995. When HIJING and
84 C PYTHIA are initialized, the shadowing is evaluated at
85 C b=0 which is the maximum. This will cause overestimate
86 C of shadowing for peripheral interactions. To correct this
87 C problem, shadowing is set to zero when initializing. Then
88 C use these maximum cross section without shadowing as a
89 C normalization of the Monte Carlo. This however increase
90 C the computing time. IHNT2(16) is used to indicate whether
91 C the sturcture function is called for (IHNT2(16)=1) initialization
92 C or for (IHNT2(16)=0)normal collisions simulation
93 C
94 C Last modification on Aagust 28, 1994. Two bugs associate
95 C with the impact parameter dependence of the shadowing is
96 C corrected.
97 C
98 C
99 c Last modification on October 14, 1994. One bug is corrected
100 c in the direct photon production option in subroutine
101 C HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
102 C Another bug associated with keeping the decay history
103 C in the particle information is also corrected.(this problem
104 C was reported by Matt Bloomer)
105 C
106 C
107 C Last modification on July 15, 1994. The option to trig on
108 C heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1)
109 C is added. To do this, set IHPR2(3)=3. For inclusive production,
110 C one should reset HIPR1(10)=0.0. One can also trig larger pt
111 C QQbar production by giving HIPR1(10) a nonvanishing value.
112 C The mass of the heavy quark in the calculation of the cross
113 C section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
114 C default is the charm mass D=1.5). We also include a separate
115 C K-factor for heavy quark and direct photon production by
116 C HIPR1(23)(D=2.0).
117 C
118 C Last modification on May 24, 1994. The option to
119 C retain the information of all particles including those
120 C who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is
121 C added to contain the line number of the parent particle
122 C of the current line which is produced via a decay.
123 C KATT(I,4) is the status number of the particle: 11=particle
124 C which has decayed; 1=finally produced particle.
125 C
126 C
127 C Last modification on May 24, 1994( in HIJSFT when valence quark
128 C is quenched, the following error is corrected. 1.2*IHNT2(1) -->
129 C 1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
130 C
131 C
132 C Last modification on March 16, 1994 (heavy flavor production
133 C processes MSUB(81)=1 MSUB(82)=1 have been switched on,
134 C charm production is the default, B-quark option is
135 C IHPR2(18), when it is switched on, charm quark is
136 C automatically off)
137 C
138 C
139 C Last modification on March 23, 1994 (an error is corrected
140 C in the impact parameter dependence of the jet cross section)
141 C
142 C Last modification Oct. 1993 to comply with non-vax
143 C machines' compiler
144 C
145 C*********************************************
146 C LAST MODIFICATION April 5, 1991
147 CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B
148 C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
149 C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
150 C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18),
151 C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
152 C
153 C gluon jets can form a single string system
154 C
155 C initial state radiation is included
156 C
157 C all QCD subprocesses are included
158 c
159 c direct particles production is included(currently only direct
160 C photon)
161 c
162 C Effect of high P_T trigger bias on multiple jets distribution
163 c
164 C******************************************************************
165 C HIJING.10 *
166 C Heavy Ion Jet INteraction Generator *
167 C by *
168 C X. N. Wang and M. Gyulassy *
169 C Lawrence Berkeley Laboratory *
170 C *
171 C******************************************************************
172 C
173 C******************************************************************
174 C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of *
175 C proj, NFP(K,4) original ID of proj. NFP(K,5)=colli status(0=no,*
176 C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
177 C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it *
178 C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
179 C (0=has not been,1=is going to be, -1=has already been scattered *
180 C NFP(k,11) total number of interactions this proj has suffered *
181 C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant *
182 C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and *
183 C diquark,PP(K,10)=PT of the hard scattering between the valence *
184 C quarks; PP(K,14,15)=the mass of quark,diquark. *
185 C******************************************************************
186 C
187 C****************************************************************
188 C
189 C SUBROUTINE HIJING
190 C
191 C****************************************************************
192  SUBROUTINE hijing(FRAME,BMIN0,BMAX0)
193  CHARACTER frame*8
194  dimension scip(300,300),rnip(300,300),sjip(300,300),jtp(3),
195  & ipcol(90000),itcol(90000)
196  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
197  SAVE /hiparnt/
198 C
199  common/hijcrdn/yp(3,300),yt(3,300)
200  SAVE /hijcrdn/
201  common/hijglbr/nelt,nint,nelp,ninp
202  SAVE /hijglbr/
203 
204 c+++BAC
205 c
206 c Add an error entry (IERRSTAT) to HIMAIN1 common block
207 c
208 c COMMON/HIMAIN1/NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
209 c
210  common/himain1/natt,eatt,jatt,nt,np,n0,n01,n10,n11,ierrstat
211 c---BAC
212  SAVE /himain1/
213 
214 C+++BAC
215 C COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
216 C SAVE /HIMAIN2/
217 C
218  common/himain2/katt(130000,4),patt(130000,4), vatt(130000,4)
219  SAVE /himain2/
220 C---BAC
221 
222  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
223  SAVE /histrng/
224  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
225  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
226  & pjpm(300,500),ntj(300),kftj(300,500),
227  & pjtx(300,500),pjty(300,500),pjtz(300,500),
228  & pjte(300,500),pjtm(300,500)
229  SAVE /hijjet1/
230  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
231  & k2sg(900,100),pxsg(900,100),pysg(900,100),
232  & pzsg(900,100),pesg(900,100),pmsg(900,100)
233  SAVE /hijjet2/
234 C+++BAC
235 C
236 C COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
237 C SAVE /HIJJET4/
238 C
239  common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5), vdr(900,5)
240  SAVE /hijjet4/
241 
242 C---BAC
243 
244  common/ranseed/nseed
245  SAVE /ranseed/
246 C
247  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
248  SAVE /lujets/
249  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
250  SAVE /ludat1/
251 
252 c+++BAC
253 c
254 c Initialize error return code
255 c
256 c---BAC
257  ierrstat = 0
258 
259  bmax=min(bmax0,hipr1(34)+hipr1(35))
260  bmin=min(bmin0,bmax)
261  IF(ihnt2(1).LE.1 .AND. ihnt2(3).LE.1) THEN
262  bmin=0.0
263  bmax=2.5*sqrt(hipr1(31)*0.1/hipr1(40))
264  ENDIF
265 C ********HIPR1(31) is in mb =0.1fm**2
266 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
267 C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
268 C
269  yp(1,1)=0.0
270  yp(2,1)=0.0
271  yp(3,1)=0.0
272  IF(ihnt2(1).LE.1) go to 14
273  DO 10 kp=1,ihnt2(1)
274 5 r=hirnd(1)
275 c
276  if(ihnt2(1).EQ.2) then
277  rnd1=max(atl_ran(nseed),1.0e-20)
278  rnd2=max(atl_ran(nseed),1.0e-20)
279  rnd3=max(atl_ran(nseed),1.0e-20)
280  r=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
281  & +4.38*0.85*log(rnd3)/(4.38+0.85))
282  endif
283 c
284  x=atl_ran(nseed)
285  cx=2.0*x-1.0
286  sx=sqrt(1.0-cx*cx)
287 C ********choose theta from uniform cos(theta) distr
288  phi=atl_ran(nseed)*2.0*hipr1(40)
289 C ********choose phi form uniform phi distr 0 to 2*pi
290  yp(1,kp)=r*sx*cos(phi)
291  yp(2,kp)=r*sx*sin(phi)
292  yp(3,kp)=r*cx
293  IF(hipr1(29).EQ.0.0) go to 10
294  DO 8 kp2=1,kp-1
295  dnbp1=(yp(1,kp)-yp(1,kp2))**2
296  dnbp2=(yp(2,kp)-yp(2,kp2))**2
297  dnbp3=(yp(3,kp)-yp(3,kp2))**2
298  dnbp=dnbp1+dnbp2+dnbp3
299  IF(dnbp.LT.hipr1(29)*hipr1(29)) go to 5
300 C ********two neighbors cannot be closer than
301 C HIPR1(29)
302 8 CONTINUE
303 10 CONTINUE
304 c*******************************
305  if(ihnt2(1).EQ.2) then
306  yp(1,2)=-yp(1,1)
307  yp(2,2)=-yp(2,1)
308  yp(3,2)=-yp(3,1)
309  endif
310 c********************************
311  DO 12 i=1,ihnt2(1)-1
312  DO 12 j=i+1,ihnt2(1)
313  IF(yp(3,i).GT.yp(3,j)) go to 12
314  y1=yp(1,i)
315  y2=yp(2,i)
316  y3=yp(3,i)
317  yp(1,i)=yp(1,j)
318  yp(2,i)=yp(2,j)
319  yp(3,i)=yp(3,j)
320  yp(1,j)=y1
321  yp(2,j)=y2
322  yp(3,j)=y3
323 12 CONTINUE
324 C
325 C******************************
326 14 yt(1,1)=0.0
327  yt(2,1)=0.0
328  yt(3,1)=0.0
329  IF(ihnt2(3).LE.1) go to 24
330  DO 20 kt=1,ihnt2(3)
331 15 r=hirnd(2)
332 c
333  if(ihnt2(3).EQ.2) then
334  rnd1=max(atl_ran(nseed),1.0e-20)
335  rnd2=max(atl_ran(nseed),1.0e-20)
336  rnd3=max(atl_ran(nseed),1.0e-20)
337  r=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
338  & +4.38*0.85*log(rnd3)/(4.38+0.85))
339  endif
340 c
341  x=atl_ran(nseed)
342  cx=2.0*x-1.0
343  sx=sqrt(1.0-cx*cx)
344 C ********choose theta from uniform cos(theta) distr
345  phi=atl_ran(nseed)*2.0*hipr1(40)
346 C ********chose phi form uniform phi distr 0 to 2*pi
347  yt(1,kt)=r*sx*cos(phi)
348  yt(2,kt)=r*sx*sin(phi)
349  yt(3,kt)=r*cx
350  IF(hipr1(29).EQ.0.0) go to 20
351  DO 18 kt2=1,kt-1
352  dnbt1=(yt(1,kt)-yt(1,kt2))**2
353  dnbt2=(yt(2,kt)-yt(2,kt2))**2
354  dnbt3=(yt(3,kt)-yt(3,kt2))**2
355  dnbt=dnbt1+dnbt2+dnbt3
356  IF(dnbt.LT.hipr1(29)*hipr1(29)) go to 15
357 C ********two neighbors cannot be closer than
358 C HIPR1(29)
359 18 CONTINUE
360 20 CONTINUE
361 c**********************************
362  if(ihnt2(3).EQ.2) then
363  yt(1,2)=-yt(1,1)
364  yt(2,2)=-yt(2,1)
365  yt(3,2)=-yt(3,1)
366  endif
367 c*********************************
368  DO 22 i=1,ihnt2(3)-1
369  DO 22 j=i+1,ihnt2(3)
370  IF(yt(3,i).LT.yt(3,j)) go to 22
371  y1=yt(1,i)
372  y2=yt(2,i)
373  y3=yt(3,i)
374  yt(1,i)=yt(1,j)
375  yt(2,i)=yt(2,j)
376  yt(3,i)=yt(3,j)
377  yt(1,j)=y1
378  yt(2,j)=y2
379  yt(3,j)=y3
380 22 CONTINUE
381 C********************
382 24 miss=-1
383 
384 50 miss=miss+1
385  IF(miss.GT.50) THEN
386  WRITE(6,*) 'infinite loop happened in HIJING'
387  stop
388  ENDIF
389 
390  natt=0
391  jatt=0
392  eatt=0.0
393  CALL hijini
394  nlop=0
395 C ********Initialize for a new event
396 60 nt=0
397  np=0
398  n0=0
399  n01=0
400  n10=0
401  n11=0
402  nelt=0
403  nint=0
404  nelp=0
405  ninp=0
406  nsg=0
407  ncolt=0
408 
409 C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
410 C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
411 C BY THE ANGLE PHI FOR EACH COLLISION.******************
412 C
413  bb=sqrt(bmin**2+atl_ran(nseed)*(bmax**2-bmin**2))
414  phi=2.0*hipr1(40)*atl_ran(nseed)
415  bbx=bb*cos(phi)
416  bby=bb*sin(phi)
417  hint1(19)=bb
418  hint1(20)=phi
419 C
420  DO 70 jp=1,ihnt2(1)
421  DO 70 jt=1,ihnt2(3)
422  scip(jp,jt)=-1.0
423  b2=(yp(1,jp)+bbx-yt(1,jt))**2+(yp(2,jp)+bby-yt(2,jt))**2
424  r2=b2*hipr1(40)/hipr1(31)/0.1
425 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
426  rrb1=min((yp(1,jp)**2+yp(2,jp)**2)
427  & /1.2**2/REAL(ihnt2(1))**0.6666667,1.0)
428  rrb2=min((yt(1,jt)**2+yt(2,jt)**2)
429  & /1.2**2/REAL(ihnt2(3))**0.6666667,1.0)
430  aphx1=hipr1(6)*4.0/3.0*(ihnt2(1)**0.3333333-1.0)
431  & *sqrt(1.0-rrb1)
432  aphx2=hipr1(6)*4.0/3.0*(ihnt2(3)**0.3333333-1.0)
433  & *sqrt(1.0-rrb2)
434  hint1(18)=hint1(14)-aphx1*hint1(15)
435  & -aphx2*hint1(16)+aphx1*aphx2*hint1(17)
436  IF(ihpr2(14).EQ.0.OR.
437  & (ihnt2(1).EQ.1.AND.ihnt2(3).EQ.1)) THEN
438  gs=1.0-exp(-(hipr1(30)+hint1(18))*romg(r2)/hipr1(31))
439  rantot=atl_ran(nseed)
440  IF(rantot.GT.gs) go to 70
441  go to 65
442  ENDIF
443  gstot_0=2.0*(1.0-exp(-(hipr1(30)+hint1(18))
444  & /hipr1(31)/2.0*romg(0.0)))
445  r2=r2/gstot_0
446  gs=1.0-exp(-(hipr1(30)+hint1(18))/hipr1(31)*romg(r2))
447  gstot=2.0*(1.0-sqrt(1.0-gs))
448  rantot=atl_ran(nseed)*gstot_0
449  IF(rantot.GT.gstot) go to 70
450  IF(rantot.GT.gs) THEN
451  CALL hijcsc(jp,jt)
452  go to 70
453 C ********perform elastic collisions
454  ENDIF
455  65 scip(jp,jt)=r2
456  rnip(jp,jt)=rantot
457  sjip(jp,jt)=hint1(18)
458  ncolt=ncolt+1
459  ipcol(ncolt)=jp
460  itcol(ncolt)=jt
461 70 CONTINUE
462 C ********total number interactions proj and targ has
463 C suffered
464  IF(ncolt.EQ.0) THEN
465  nlop=nlop+1
466  IF(nlop.LE.20.OR.
467  & (ihnt2(1).EQ.1.AND.ihnt2(3).EQ.1)) go to 60
468  RETURN
469  ENDIF
470 C ********At large impact parameter, there maybe no
471 C interaction at all. For NN collision
472 C repeat the event until interaction happens
473 C
474  IF(ihpr2(3).NE.0) THEN
475  nhard=1+int(atl_ran(nseed)*(ncolt-1)+0.5)
476  nhard=min(nhard,ncolt)
477  jphard=ipcol(nhard)
478  jthard=itcol(nhard)
479  ENDIF
480 C
481  IF(ihpr2(9).EQ.1) THEN
482  nmini=1+int(atl_ran(nseed)*(ncolt-1)+0.5)
483  nmini=min(nmini,ncolt)
484  jpmini=ipcol(nmini)
485  jtmini=itcol(nmini)
486  ENDIF
487 C ********Specifying the location of the hard and
488 C minijet if they are enforced by user
489 C
490  DO 200 jp=1,ihnt2(1)
491  DO 200 jt=1,ihnt2(3)
492  IF(scip(jp,jt).EQ.-1.0) go to 200
493  nfp(jp,11)=nfp(jp,11)+1
494  nft(jt,11)=nft(jt,11)+1
495  IF(nfp(jp,5).LE.1 .AND. nft(jt,5).GT.1) THEN
496  np=np+1
497  n01=n01+1
498  ELSE IF(nfp(jp,5).GT.1 .AND. nft(jt,5).LE.1) THEN
499  nt=nt+1
500  n10=n10+1
501  ELSE IF(nfp(jp,5).LE.1 .AND. nft(jt,5).LE.1) THEN
502  np=np+1
503  nt=nt+1
504  n0=n0+1
505  ELSE IF(nfp(jp,5).GT.1 .AND. nft(jt,5).GT.1) THEN
506  n11=n11+1
507  ENDIF
508 c
509  jout=0
510  nfp(jp,10)=0
511  nft(jt,10)=0
512 C*****************************************************************
513  IF(ihpr2(8).EQ.0 .AND. ihpr2(3).EQ.0) go to 160
514 C ********When IHPR2(8)=0 no jets are produced
515  IF(nfp(jp,6).LT.0 .OR. nft(jt,6).LT.0) go to 160
516 C ********jets can not be produced for (JP,JT)
517 C because not enough energy avaible for
518 C JP or JT
519  r2=scip(jp,jt)
520  hint1(18)=sjip(jp,jt)
521  tt=romg(r2)*hint1(18)/hipr1(31)
522  tts=hipr1(30)*romg(r2)/hipr1(31)
523  njet=0
524  IF(ihpr2(3).NE.0 .AND. jp.EQ.jphard .AND. jt.EQ.jthard) THEN
525  CALL jetini(jp,jt,1)
526  CALL hijhrd(jp,jt,0,jflg,0)
527  hint1(26)=hint1(47)
528  hint1(27)=hint1(48)
529  hint1(28)=hint1(49)
530  hint1(29)=hint1(50)
531  hint1(36)=hint1(67)
532  hint1(37)=hint1(68)
533  hint1(38)=hint1(69)
534  hint1(39)=hint1(70)
535 C
536  IF(abs(hint1(46)).GT.hipr1(11).AND.jflg.EQ.2) nfp(jp,7)=1
537  IF(abs(hint1(56)).GT.hipr1(11).AND.jflg.EQ.2) nft(jt,7)=1
538  IF(max(abs(hint1(46)),abs(hint1(56))).GT.hipr1(11).AND.
539  & jflg.GE.3) iasg(nsg,3)=1
540  ihnt2(9)=ihnt2(14)
541  ihnt2(10)=ihnt2(15)
542  DO 105 i05=1,5
543  hint1(20+i05)=hint1(40+i05)
544  hint1(30+i05)=hint1(50+i05)
545  105 CONTINUE
546  jout=1
547  IF(ihpr2(8).EQ.0) go to 160
548  rrb1=min((yp(1,jp)**2+yp(2,jp)**2)/1.2**2
549  & /REAL(ihnt2(1))**0.6666667,1.0)
550  rrb2=min((yt(1,jt)**2+yt(2,jt)**2)/1.2**2
551  & /REAL(ihnt2(3))**0.6666667,1.0)
552  aphx1=hipr1(6)*4.0/3.0*(ihnt2(1)**0.3333333-1.0)
553  & *sqrt(1.0-rrb1)
554  aphx2=hipr1(6)*4.0/3.0*(ihnt2(3)**0.3333333-1.0)
555  & *sqrt(1.0-rrb2)
556  hint1(65)=hint1(61)-aphx1*hint1(62)
557  & -aphx2*hint1(63)+aphx1*aphx2*hint1(64)
558  ttrig=romg(r2)*hint1(65)/hipr1(31)
559  njet=-1
560 C ********subtract the trigger jet from total number
561 C of jet production to be done since it has
562 C already been produced here
563  xr1=-alog(exp(-ttrig)+atl_ran(nseed)*(1.0-exp(-ttrig)))
564  106 njet=njet+1
565  xr1=xr1-alog(max(atl_ran(nseed),1.0e-20))
566  IF(xr1.LT.ttrig) go to 106
567  xr=0.0
568  107 njet=njet+1
569  xr=xr-alog(max(atl_ran(nseed),1.0e-20))
570  IF(xr.LT.tt-ttrig) go to 107
571  njet=njet-1
572  go to 112
573  ENDIF
574 C ********create a hard interaction with specified P_T
575 c when IHPR2(3)>0
576  IF(ihpr2(9).EQ.1.AND.jp.EQ.jpmini.AND.jt.EQ.jtmini) go to 110
577 C ********create at least one pair of mini jets
578 C when IHPR2(9)=1
579 C
580  IF(ihpr2(8).GT.0 .AND.rnip(jp,jt).LT.exp(-tt)*
581  & (1.0-exp(-tts))) go to 160
582 C ********this is the probability for no jet production
583 110 xr=-alog(exp(-tt)+atl_ran(nseed)*(1.0-exp(-tt)))
584 111 njet=njet+1
585  xr=xr-alog(max(atl_ran(nseed),1.0e-20))
586  IF(xr.LT.tt) go to 111
587 112 njet=min(njet,ihpr2(8))
588  IF(ihpr2(8).LT.0) njet=abs(ihpr2(8))
589 C ******** Determine number of mini jet production
590 C
591  DO 150 i_jet=1,njet
592  CALL jetini(jp,jt,0)
593  CALL hijhrd(jp,jt,jout,jflg,1)
594 C ********JFLG=1 jets valence quarks, JFLG=2 with
595 C gluon jet, JFLG=3 with q-qbar prod for
596 C (JP,JT). If JFLG=0 jets can not be produced
597 C this time. If JFLG=-1, error occured abandon
598 C this event. JOUT is the total hard scat for
599 C (JP,JT) up to now.
600  IF(jflg.EQ.0) go to 160
601  IF(jflg.LT.0) THEN
602  IF(ihpr2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
603  go to 50
604  ENDIF
605  jout=jout+1
606  IF(abs(hint1(46)).GT.hipr1(11).AND.jflg.EQ.2) nfp(jp,7)=1
607  IF(abs(hint1(56)).GT.hipr1(11).AND.jflg.EQ.2) nft(jt,7)=1
608  IF(max(abs(hint1(46)),abs(hint1(56))).GT.hipr1(11).AND.
609  & jflg.GE.3) iasg(nsg,3)=1
610 C ******** jet with PT>HIPR1(11) will be quenched
611  150 CONTINUE
612  160 CONTINUE
613  CALL hijsft(jp,jt,jout,ierror)
614  IF(ierror.NE.0) THEN
615  IF(ihpr2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
616  go to 50
617  ENDIF
618 C
619 C ********conduct soft scattering between JP and JT
620  jatt=jatt+jout
621 
622 200 CONTINUE
623 c
624 c**************************
625 c
626  DO 201 jp=1,ihnt2(1)
627  IF(nfp(jp,5).GT.2) THEN
628  ninp=ninp+1
629  ELSE IF(nfp(jp,5).EQ.2.OR.nfp(jp,5).EQ.1) THEN
630  nelp=nelp+1
631  ENDIF
632  201 continue
633  DO 202 jt=1,ihnt2(3)
634  IF(nft(jt,5).GT.2) THEN
635  nint=nint+1
636  ELSE IF(nft(jt,5).EQ.2.OR.nft(jt,5).EQ.1) THEN
637  nelt=nelt+1
638  ENDIF
639  202 continue
640 c
641 c*******************************
642 
643 
644 C********perform jet quenching for jets with PT>HIPR1(11)**********
645 
646  IF((ihpr2(8).NE.0.OR.ihpr2(3).NE.0).AND.ihpr2(4).GT.0.AND.
647  & ihnt2(1).GT.1.AND.ihnt2(3).GT.1) THEN
648  DO 271 i=1,ihnt2(1)
649  IF(nfp(i,7).EQ.1) CALL quench(i,1)
650 271 CONTINUE
651  DO 272 i=1,ihnt2(3)
652  IF(nft(i,7).EQ.1) CALL quench(i,2)
653 272 CONTINUE
654  DO 273 isg=1,nsg
655  IF(iasg(isg,3).EQ.1) CALL quench(isg,3)
656 273 CONTINUE
657  ENDIF
658 C
659 C**************fragment all the string systems in the following*****
660 C
661 C********N_ST is where particle information starts
662 C********N_STR+1 is the number of strings in fragmentation
663 C********the number of strings before a line is stored in K(I,4)
664 C********IDSTR is id number of the string system (91,92 or 93)
665 C
666  IF(ihpr2(20).NE.0) THEN
667  DO 360 isg=1,nsg
668  CALL hijfrg(isg,3,ierror)
669  IF(mstu(24).NE.0 .OR.ierror.GT.0) THEN
670  mstu(24)=0
671  mstu(28)=0
672  IF(ihpr2(10).NE.0) THEN
673  call lulist(1)
674  WRITE(6,*) 'error occured, repeat the event'
675  ENDIF
676  go to 50
677  ENDIF
678 C ********Check errors
679 C
680 C
681 C DPM: call fastjet and record "jets" in LUJETS
682 C
683 c NSAVE = N
684  CALL hijfst(n,9000,k,p,v)
685 c WRITE(*,*)
686 c WRITE(*,*) N-NSAVE, ' "jets" produced'
687 c WRITE(*,*) ' id px py pz E m '
688 c WRITE(*,*) '==========================================='
689 c DO I = NSAVE+1, N
690 c WRITE(*,2718) K(I,1), (P(I,J), J=1,5)
691 c ENDDO
692 c 2718 FORMAT(I4,5(2X,F6.3))
693 C
694 C
695 C
696 
697  n_st=1
698  idstr=92
699  IF(ihpr2(21).EQ.0) THEN
700  CALL luedit(2)
701 c+++BAC
702 c
703 c Don't perform the search for string if N=1 because it will search beyond the end of the valid particle list
704 c
705 c ELSE
706 c
707  ELSE IF (n .GT. 1) THEN
708 c---BAC
709 
710 351 n_st=n_st+1
711 
712 c+++BAC
713 c
714 c Check for inconsistency -- no string line found
715 c
716 c---BAC
717  if (n_st .GT. n) then
718  ierrstat = 2
719  RETURN
720  ENDIF
721 
722  IF(k(n_st,2).LT.91.OR.k(n_st,2).GT.93) go to 351
723  idstr=k(n_st,2)
724  n_st=n_st+1
725  ENDIF
726 C
727  IF(frame.EQ.'LAB') THEN
728  CALL hiboost
729  ENDIF
730 C ******** boost back to lab frame(if it was in)
731 C
732  n_str=0
733  DO 361 i=n_st,n
734  IF(k(i,2).EQ.idstr) THEN
735 C+++BAC
736  IF (k(i,3) .LT. n_st) THEN
737 C---BAC
738  n_str=n_str+1
739  go to 360
740  ENDIF
741  ENDIF
742 
743  k(i,4)=n_str
744  natt=natt+1
745 c BAC+++
746 c
747 c Add a check on array overflow
748 c
749 c BAC---
750  if (natt .GT. 130000) THEN
751  ierrstat = 1
752  RETURN
753  ENDIF
754 
755  katt(natt,1)=k(i,2)
756  katt(natt,2)=20
757  katt(natt,4)=k(i,1)
758 C+++BAC
759 C IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
760  IF(k(i,3).EQ.0 .OR.
761  & k(i,3).LT.n_st .OR.
762  & (k(k(i,3),2) .EQ. idstr .AND.
763  & k(k(i,3),3) .LT. n_st)) THEN
764 C---BAC
765  katt(natt,3)=0
766  ELSE
767  katt(natt,3)=natt-i+k(i,3)+n_str-k(k(i,3),4)
768  ENDIF
769 C ****** identify the mother particle
770  patt(natt,1)=p(i,1)
771  patt(natt,2)=p(i,2)
772  patt(natt,3)=p(i,3)
773  patt(natt,4)=p(i,4)
774  eatt=eatt+p(i,4)
775 
776  vatt(natt,1)=v(i,1)
777  vatt(natt,2)=v(i,2)
778  vatt(natt,3)=v(i,3)
779 
780  IF ((abs(vatt(natt,3)) .GT. 0.00001) .and.
781  & (katt(natt,3) .eq. 0 )) THEN
782  CALL lulist(3)
783  ENDIF
784 
785 C+++BAC
786  kparent = katt(natt,3)
787  if (kparent .ne. 0) then
788  r = sqrt(vatt(natt,1)**2 + vatt(natt,2)**2 +
789  & vatt(natt,3)**2)
790 
791  rparent = sqrt(vatt(kparent,1)**2 +
792  & vatt(kparent,2)**2 +
793  & vatt(kparent,3)**2)
794  IF (r/rparent .LT. 0.85) THEN
795  CALL lulist(3)
796  ENDIF
797  ENDIF
798 C---BAC
799 
800  vatt(natt,4)=v(i,4)
801  361 CONTINUE
802  360 CONTINUE
803 C ********Fragment the q-qbar jets systems *****
804 C
805  jtp(1)=ihnt2(1)
806  jtp(2)=ihnt2(3)
807  DO 400 ntp=1,2
808  DO 400 j_jtp=1,jtp(ntp)
809  CALL hijfrg(j_jtp,ntp,ierror)
810  IF(mstu(24).NE.0 .OR. ierror.GT.0) THEN
811  mstu(24)=0
812  mstu(28)=0
813  IF(ihpr2(10).NE.0) THEN
814  call lulist(1)
815  WRITE(6,*) 'error occured, repeat the event'
816  ENDIF
817  go to 50
818  ENDIF
819 C ********check errors
820 C
821 C
822 C DPM: call fastjet
823 C
824  CALL hijfst(n,9000,k,p,v)
825 C
826 C
827 C
828 
829  n_st=1
830  idstr=92
831  IF(ihpr2(21).EQ.0) THEN
832  CALL luedit(2)
833 c+++BAC
834 c
835 c Don't perform the search for string if N=1 because it will search beyond the end of the valid particle list
836 c
837 c ELSE
838 c
839  ELSE IF (n .GT. 1) THEN
840 c---BAC
841 381 n_st=n_st+1
842 
843 c+++BAC
844 c
845 c Check for inconsistency -- no string line found
846 c
847 c---BAC
848  if (n_st .GT. n) then
849  print *, 'inconsistency, n = ', n, ', n_st = ', n_st
850  ierrstat = 2
851  RETURN
852  ENDIF
853 
854  IF(k(n_st,2).LT.91.OR.k(n_st,2).GT.93) go to 381
855  idstr=k(n_st,2)
856  n_st=n_st+1
857  ENDIF
858 
859  IF(frame.EQ.'LAB') THEN
860  CALL hiboost
861  ENDIF
862 C ******** boost back to lab frame(if it was in)
863 C
864  nftp=nfp(j_jtp,5)
865  IF(ntp.EQ.2) nftp=10+nft(j_jtp,5)
866  n_str=0
867  DO 390 i=n_st,n
868  IF(k(i,2).EQ.idstr) THEN
869 C+++BAC
870  IF (k(i,3) .LT. n_st) THEN
871 C---BAC
872  n_str=n_str+1
873  go to 390
874  ENDIF
875  ENDIF
876  k(i,4)=n_str
877  natt=natt+1
878 
879 c BAC+++
880 c
881 c Add a check on array overflow
882 c
883 c BAC---
884  if (natt .GT. 130000) THEN
885  ierrstat = 1
886  RETURN
887  ENDIF
888 
889  katt(natt,1)=k(i,2)
890  katt(natt,2)=nftp
891  katt(natt,4)=k(i,1)
892 C+++BAC
893 C IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
894  IF(k(i,3).EQ.0 .OR.
895  & k(i,3).LT.n_st .OR.
896  & (k(k(i,3),2) .EQ. idstr .AND.
897  & k(k(i,3),3) .LT. n_st)) THEN
898 C---BAC
899  katt(natt,3)=0
900  ELSE
901  katt(natt,3)=natt-i+k(i,3)+n_str-k(k(i,3),4)
902  ENDIF
903 C ****** identify the mother particle
904  patt(natt,1)=p(i,1)
905  patt(natt,2)=p(i,2)
906  patt(natt,3)=p(i,3)
907  patt(natt,4)=p(i,4)
908  eatt=eatt+p(i,4)
909 
910  vatt(natt,1)=v(i,1)
911  vatt(natt,2)=v(i,2)
912  vatt(natt,3)=v(i,3)
913 
914  if ((abs(vatt(natt,3)) .gt. 0.00001) .and.
915  & (katt(natt,3) .eq. 0 )) then
916  call lulist(3)
917  endif
918 
919 C+++BAC
920  kparent = katt(natt,3)
921  if (kparent .ne. 0) then
922  r = sqrt(vatt(natt,1)**2 + vatt(natt,2)**2 +
923  & vatt(natt,3)**2)
924 
925  rparent = sqrt(vatt(kparent,1)**2 +
926  & vatt(kparent,2)**2 +
927  & vatt(kparent,3)**2)
928  IF (r/rparent .LT. 0.85) THEN
929  CALL lulist(3)
930  ENDIF
931  ENDIF
932 C---BAC
933 
934  vatt(natt,4)=v(i,4)
935 390 CONTINUE
936 400 CONTINUE
937 C ********Fragment the q-qq related string systems
938  ENDIF
939 
940  DO 450 i=1,ndr
941  natt=natt+1
942 
943 c BAC+++
944 c
945 c Add a check on array overflow
946 c
947 c BAC---
948  if (natt .GT. 130000) THEN
949  ierrstat = 1
950  RETURN
951  ENDIF
952 
953  katt(natt,1)=kfdr(i)
954  katt(natt,2)=40
955  katt(natt,3)=0
956  patt(natt,1)=pdr(i,1)
957  patt(natt,2)=pdr(i,2)
958  patt(natt,3)=pdr(i,3)
959  patt(natt,4)=pdr(i,4)
960  eatt=eatt+pdr(i,4)
961 
962  vatt(natt,1)=vdr(i,1)
963  vatt(natt,2)=vdr(i,2)
964  vatt(natt,3)=vdr(i,3)
965  vatt(natt,4)=vdr(i,4)
966 450 CONTINUE
967 C ********store the direct-produced particles
968 C
969  dengy=eatt/(ihnt2(1)*hint1(6)+ihnt2(3)*hint1(7))-1.0
970  IF(abs(dengy).GT.hipr1(43).AND.ihpr2(20).NE.0
971  & .AND.ihpr2(21).EQ.0) THEN
972  IF(ihpr2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat event'
973 C call lulist(1)
974  go to 50
975  ENDIF
976  RETURN
977  END