ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hijing1.383.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hijing1.383.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 c SUBROUTINE HIJING(FRAME,BMIN0,BMAX0) !khaled
193  SUBROUTINE hijing(BMIN0,BMAX0)
194 
195  CHARACTER frame*8
196  dimension scip(300,300),rnip(300,300),sjip(300,300),jtp(3),
197  & ipcol(90000),itcol(90000)
198  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
199  SAVE /hiparnt/
200 C
201  common/hijcrdn/yp(3,300),yt(3,300)
202  SAVE /hijcrdn/
203  common/hijglbr/nelt,nint,nelp,ninp
204  SAVE /hijglbr/
205  common/himain1/natt,eatt,jatt,nt,np,n0,n01,n10,n11
206  SAVE /himain1/
207  common/himain2/katt(130000,4),patt(130000,4)
208  SAVE /himain2/
209  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
210  SAVE /histrng/
211  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
212  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
213  & pjpm(300,500),ntj(300),kftj(300,500),
214  & pjtx(300,500),pjty(300,500),pjtz(300,500),
215  & pjte(300,500),pjtm(300,500)
216  SAVE /hijjet1/
217  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
218  & k2sg(900,100),pxsg(900,100),pysg(900,100),
219  & pzsg(900,100),pesg(900,100),pmsg(900,100)
220  SAVE /hijjet2/
221  common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5)
222  SAVE /hijjet4/
223  common/ranseed/nseed
224  SAVE /ranseed/
225 C
226  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
227  SAVE /lujets/
228  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
229  SAVE /ludat1/
230 
231  bmax=min(bmax0,hipr1(34)+hipr1(35))
232  bmin=min(bmin0,bmax)
233  IF(ihnt2(1).LE.1 .AND. ihnt2(3).LE.1) THEN
234  bmin=0.0
235  bmax=2.5*sqrt(hipr1(31)*0.1/hipr1(40))
236 c write(*,*)"HIJING bmin bmax",BMIN,"=====", BMAX
237  ENDIF
238 c ----------------------------
239  frame='CMS' !khaled //CMS
240 C ********HIPR1(31) is in mb =0.1fm**2
241 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
242 C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
243 C
244  yp(1,1)=0.0
245  yp(2,1)=0.0
246  yp(3,1)=0.0
247  IF(ihnt2(1).LE.1) go to 14
248  DO 10 kp=1,ihnt2(1)
249 5 r=hirnd(1)
250 c
251  if(ihnt2(1).EQ.2) then
252  rnd1=max(rlu(0),1.0e-20)
253  rnd2=max(rlu(0),1.0e-20)
254  rnd3=max(rlu(0),1.0e-20)
255  r=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
256  & +4.38*0.85*log(rnd3)/(4.38+0.85))
257  endif
258 c
259  x=rlu(0)
260  cx=2.0*x-1.0
261  sx=sqrt(1.0-cx*cx)
262 C ********choose theta from uniform cos(theta) distr
263  phi=rlu(0)*2.0*hipr1(40)
264 C ********choose phi form uniform phi distr 0 to 2*pi
265  yp(1,kp)=r*sx*cos(phi)
266  yp(2,kp)=r*sx*sin(phi)
267  yp(3,kp)=r*cx
268  IF(hipr1(29).EQ.0.0) go to 10
269  DO 8 kp2=1,kp-1
270  dnbp1=(yp(1,kp)-yp(1,kp2))**2
271  dnbp2=(yp(2,kp)-yp(2,kp2))**2
272  dnbp3=(yp(3,kp)-yp(3,kp2))**2
273  dnbp=dnbp1+dnbp2+dnbp3
274  IF(dnbp.LT.hipr1(29)*hipr1(29)) go to 5
275 C ********two neighbors cannot be closer than
276 C HIPR1(29)
277 8 CONTINUE
278 10 CONTINUE
279 c*******************************
280  if(ihnt2(1).EQ.2) then
281  yp(1,2)=-yp(1,1)
282  yp(2,2)=-yp(2,1)
283  yp(3,2)=-yp(3,1)
284  endif
285 c********************************
286  DO 12 i=1,ihnt2(1)-1
287  DO 12 j=i+1,ihnt2(1)
288  IF(yp(3,i).GT.yp(3,j)) go to 12
289  y1=yp(1,i)
290  y2=yp(2,i)
291  y3=yp(3,i)
292  yp(1,i)=yp(1,j)
293  yp(2,i)=yp(2,j)
294  yp(3,i)=yp(3,j)
295  yp(1,j)=y1
296  yp(2,j)=y2
297  yp(3,j)=y3
298 12 CONTINUE
299 C
300 C******************************
301 14 yt(1,1)=0.0
302  yt(2,1)=0.0
303  yt(3,1)=0.0
304  IF(ihnt2(3).LE.1) go to 24
305  DO 20 kt=1,ihnt2(3)
306 15 r=hirnd(2)
307 c
308  if(ihnt2(3).EQ.2) then
309  rnd1=max(rlu(0),1.0e-20)
310  rnd2=max(rlu(0),1.0e-20)
311  rnd3=max(rlu(0),1.0e-20)
312  r=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
313  & +4.38*0.85*log(rnd3)/(4.38+0.85))
314  endif
315 c
316  x=rlu(0)
317  cx=2.0*x-1.0
318  sx=sqrt(1.0-cx*cx)
319 C ********choose theta from uniform cos(theta) distr
320  phi=rlu(0)*2.0*hipr1(40)
321 C ********chose phi form uniform phi distr 0 to 2*pi
322  yt(1,kt)=r*sx*cos(phi)
323  yt(2,kt)=r*sx*sin(phi)
324  yt(3,kt)=r*cx
325  IF(hipr1(29).EQ.0.0) go to 20
326  DO 18 kt2=1,kt-1
327  dnbt1=(yt(1,kt)-yt(1,kt2))**2
328  dnbt2=(yt(2,kt)-yt(2,kt2))**2
329  dnbt3=(yt(3,kt)-yt(3,kt2))**2
330  dnbt=dnbt1+dnbt2+dnbt3
331  IF(dnbt.LT.hipr1(29)*hipr1(29)) go to 15
332 C ********two neighbors cannot be closer than
333 C HIPR1(29)
334 18 CONTINUE
335 20 CONTINUE
336 c**********************************
337  if(ihnt2(3).EQ.2) then
338  yt(1,2)=-yt(1,1)
339  yt(2,2)=-yt(2,1)
340  yt(3,2)=-yt(3,1)
341  endif
342 c*********************************
343  DO 22 i=1,ihnt2(3)-1
344  DO 22 j=i+1,ihnt2(3)
345  IF(yt(3,i).LT.yt(3,j)) go to 22
346  y1=yt(1,i)
347  y2=yt(2,i)
348  y3=yt(3,i)
349  yt(1,i)=yt(1,j)
350  yt(2,i)=yt(2,j)
351  yt(3,i)=yt(3,j)
352  yt(1,j)=y1
353  yt(2,j)=y2
354  yt(3,j)=y3
355 22 CONTINUE
356 C********************
357 24 miss=-1
358 
359 50 miss=miss+1
360  IF(miss.GT.50) THEN
361  WRITE(6,*) 'infinite loop happened in HIJING'
362  stop
363  ENDIF
364 
365  natt=0
366  jatt=0
367  eatt=0.0
368  CALL hijini
369  nlop=0
370 C ********Initialize for a new event
371 60 nt=0
372  np=0
373  n0=0
374  n01=0
375  n10=0
376  n11=0
377  nelt=0
378  nint=0
379  nelp=0
380  ninp=0
381  nsg=0
382  ncolt=0
383 
384 C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
385 C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
386 C BY THE ANGLE PHI FOR EACH COLLISION.******************
387 C
388 
389  bb=sqrt(bmin**2+rlu(0)*(bmax**2-bmin**2))
390 c write(*,*)"HIJING bmin bmax bb", bmin, bmax,bb
391  phi=2.0*hipr1(40)*rlu(0)
392  bbx=bb*cos(phi)
393  bby=bb*sin(phi)
394  hint1(19)=bb
395  hint1(20)=phi
396 C
397  DO 70 jp=1,ihnt2(1)
398  DO 70 jt=1,ihnt2(3)
399  scip(jp,jt)=-1.0
400  b2=(yp(1,jp)+bbx-yt(1,jt))**2+(yp(2,jp)+bby-yt(2,jt))**2
401  r2=b2*hipr1(40)/hipr1(31)/0.1
402 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
403  rrb1=min((yp(1,jp)**2+yp(2,jp)**2)
404  & /1.2**2/REAL(ihnt2(1))**0.6666667,1.0)
405  rrb2=min((yt(1,jt)**2+yt(2,jt)**2)
406  & /1.2**2/REAL(ihnt2(3))**0.6666667,1.0)
407  aphx1=hipr1(6)*4.0/3.0*(ihnt2(1)**0.3333333-1.0)
408  & *sqrt(1.0-rrb1)
409  aphx2=hipr1(6)*4.0/3.0*(ihnt2(3)**0.3333333-1.0)
410  & *sqrt(1.0-rrb2)
411  hint1(18)=hint1(14)-aphx1*hint1(15)
412  & -aphx2*hint1(16)+aphx1*aphx2*hint1(17)
413 
414  IF(ihpr2(14).EQ.0.OR.
415  & (ihnt2(1).EQ.1.AND.ihnt2(3).EQ.1)) THEN
416 c write(*,*)'R2 ROMG(R2)',R2, ROMG(R2), OMG0(R2) KHal
417  gs=1.0-exp(-(hipr1(30)+hint1(18))*romg(r2)/hipr1(31))
418  rantot=rlu(0)
419 
420  IF(rantot.GT.gs) go to 70
421  go to 65
422  ENDIF
423 c write(*,*)'ROMG(0.0)' !kha
424  gstot_0=2.0*(1.0-exp(-(hipr1(30)+hint1(18))
425  & /hipr1(31)/2.0*romg(0.0)))
426 c write(*,*)GSTOT_0, ROMG(0.0)
427  r2=r2/gstot_0
428  gs=1.0-exp(-(hipr1(30)+hint1(18))/hipr1(31)*romg(r2))
429  gstot=2.0*(1.0-sqrt(1.0-gs))
430  rantot=rlu(0)*gstot_0
431 c write(*,*)'rantot, gstot',rantot, gstot
432  IF(rantot.GT.gstot) go to 70
433  IF(rantot.GT.gs) THEN
434 c write(*,*)"----------------------------------start HIJCSC"
435  CALL hijcsc(jp,jt)
436 c write(*,*)"----------------------------------end HIJCSC"
437  go to 70
438 C ********perform elastic collisions
439  ENDIF
440  65 scip(jp,jt)=r2
441  rnip(jp,jt)=rantot
442  sjip(jp,jt)=hint1(18)
443  ncolt=ncolt+1
444  ipcol(ncolt)=jp
445  itcol(ncolt)=jt
446 70 CONTINUE
447 C ********total number interactions proj and targ has
448 C suffered
449  IF(ncolt.EQ.0) THEN
450  nlop=nlop+1
451  IF(nlop.LE.20.OR.
452  & (ihnt2(1).EQ.1.AND.ihnt2(3).EQ.1)) go to 60
453  RETURN
454  ENDIF
455 C ********At large impact parameter, there maybe no
456 C interaction at all. For NN collision
457 C repeat the event until interaction happens
458 C
459  IF(ihpr2(3).NE.0) THEN
460  nhard=1+int(rlu(0)*(ncolt-1)+0.5)
461  nhard=min(nhard,ncolt)
462  jphard=ipcol(nhard)
463  jthard=itcol(nhard)
464  ENDIF
465 C
466  IF(ihpr2(9).EQ.1) THEN
467  nmini=1+int(rlu(0)*(ncolt-1)+0.5)
468  nmini=min(nmini,ncolt)
469  jpmini=ipcol(nmini)
470  jtmini=itcol(nmini)
471  ENDIF
472 C ********Specifying the location of the hard and
473 C minijet if they are enforced by user
474 C
475  DO 200 jp=1,ihnt2(1)
476  DO 200 jt=1,ihnt2(3)
477  IF(scip(jp,jt).EQ.-1.0) go to 200
478  nfp(jp,11)=nfp(jp,11)+1
479  nft(jt,11)=nft(jt,11)+1
480  IF(nfp(jp,5).LE.1 .AND. nft(jt,5).GT.1) THEN
481  np=np+1
482  n01=n01+1
483  ELSE IF(nfp(jp,5).GT.1 .AND. nft(jt,5).LE.1) THEN
484  nt=nt+1
485  n10=n10+1
486  ELSE IF(nfp(jp,5).LE.1 .AND. nft(jt,5).LE.1) THEN
487  np=np+1
488  nt=nt+1
489  n0=n0+1
490  ELSE IF(nfp(jp,5).GT.1 .AND. nft(jt,5).GT.1) THEN
491  n11=n11+1
492  ENDIF
493 c
494  jout=0
495  nfp(jp,10)=0
496  nft(jt,10)=0
497 C*****************************************************************
498  IF(ihpr2(8).EQ.0 .AND. ihpr2(3).EQ.0) go to 160
499 C ********When IHPR2(8)=0 no jets are produced
500  IF(nfp(jp,6).LT.0 .OR. nft(jt,6).LT.0) go to 160
501 C ********jets can not be produced for (JP,JT)
502 C because not enough energy avaible for
503 C JP or JT
504  r2=scip(jp,jt)
505  hint1(18)=sjip(jp,jt)
506  tt=romg(r2)*hint1(18)/hipr1(31)
507  tts=hipr1(30)*romg(r2)/hipr1(31)
508  njet=0
509  IF(ihpr2(3).NE.0 .AND. jp.EQ.jphard .AND. jt.EQ.jthard) THEN
510  CALL jetini(jp,jt,1)
511  CALL hijhrd(jp,jt,0,jflg,0)
512  hint1(26)=hint1(47)
513  hint1(27)=hint1(48)
514  hint1(28)=hint1(49)
515  hint1(29)=hint1(50)
516  hint1(36)=hint1(67)
517  hint1(37)=hint1(68)
518  hint1(38)=hint1(69)
519  hint1(39)=hint1(70)
520 C
521  IF(abs(hint1(46)).GT.hipr1(11).AND.jflg.EQ.2) nfp(jp,7)=1
522  IF(abs(hint1(56)).GT.hipr1(11).AND.jflg.EQ.2) nft(jt,7)=1
523  IF(max(abs(hint1(46)),abs(hint1(56))).GT.hipr1(11).AND.
524  & jflg.GE.3) iasg(nsg,3)=1
525  ihnt2(9)=ihnt2(14)
526  ihnt2(10)=ihnt2(15)
527  DO 105 i05=1,5
528  hint1(20+i05)=hint1(40+i05)
529  hint1(30+i05)=hint1(50+i05)
530  105 CONTINUE
531  jout=1
532  IF(ihpr2(8).EQ.0) go to 160
533  rrb1=min((yp(1,jp)**2+yp(2,jp)**2)/1.2**2
534  & /REAL(ihnt2(1))**0.6666667,1.0)
535  rrb2=min((yt(1,jt)**2+yt(2,jt)**2)/1.2**2
536  & /REAL(ihnt2(3))**0.6666667,1.0)
537  aphx1=hipr1(6)*4.0/3.0*(ihnt2(1)**0.3333333-1.0)
538  & *sqrt(1.0-rrb1)
539  aphx2=hipr1(6)*4.0/3.0*(ihnt2(3)**0.3333333-1.0)
540  & *sqrt(1.0-rrb2)
541  hint1(65)=hint1(61)-aphx1*hint1(62)
542  & -aphx2*hint1(63)+aphx1*aphx2*hint1(64)
543  ttrig=romg(r2)*hint1(65)/hipr1(31)
544  njet=-1
545 C ********subtract the trigger jet from total number
546 C of jet production to be done since it has
547 C already been produced here
548  xr1=-alog(exp(-ttrig)+rlu(0)*(1.0-exp(-ttrig)))
549  106 njet=njet+1
550  xr1=xr1-alog(max(rlu(0),1.0e-20))
551  IF(xr1.LT.ttrig) go to 106
552  xr=0.0
553  107 njet=njet+1
554  xr=xr-alog(max(rlu(0),1.0e-20))
555  IF(xr.LT.tt-ttrig) go to 107
556  njet=njet-1
557  go to 112
558  ENDIF
559 C ********create a hard interaction with specified P_T
560 c when IHPR2(3)>0
561  IF(ihpr2(9).EQ.1.AND.jp.EQ.jpmini.AND.jt.EQ.jtmini) go to 110
562 C ********create at least one pair of mini jets
563 C when IHPR2(9)=1
564 C
565  IF(ihpr2(8).GT.0 .AND.rnip(jp,jt).LT.exp(-tt)*
566  & (1.0-exp(-tts))) go to 160
567 C ********this is the probability for no jet production
568 110 xr=-alog(exp(-tt)+rlu(0)*(1.0-exp(-tt)))
569 111 njet=njet+1
570  xr=xr-alog(max(rlu(0),1.0e-20))
571  IF(xr.LT.tt) go to 111
572 112 njet=min(njet,ihpr2(8))
573  IF(ihpr2(8).LT.0) njet=abs(ihpr2(8))
574 C ******** Determine number of mini jet production
575 C
576  DO 150 i_jet=1,njet
577  CALL jetini(jp,jt,0)
578  CALL hijhrd(jp,jt,jout,jflg,1)
579 C ********JFLG=1 jets valence quarks, JFLG=2 with
580 C gluon jet, JFLG=3 with q-qbar prod for
581 C (JP,JT). If JFLG=0 jets can not be produced
582 C this time. If JFLG=-1, error occured abandon
583 C this event. JOUT is the total hard scat for
584 C (JP,JT) up to now.
585  IF(jflg.EQ.0) go to 160
586  IF(jflg.LT.0) THEN
587  IF(ihpr2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
588  go to 50
589  ENDIF
590  jout=jout+1
591  IF(abs(hint1(46)).GT.hipr1(11).AND.jflg.EQ.2) nfp(jp,7)=1
592  IF(abs(hint1(56)).GT.hipr1(11).AND.jflg.EQ.2) nft(jt,7)=1
593  IF(max(abs(hint1(46)),abs(hint1(56))).GT.hipr1(11).AND.
594  & jflg.GE.3) iasg(nsg,3)=1
595 C ******** jet with PT>HIPR1(11) will be quenched
596  150 CONTINUE
597  160 CONTINUE
598  CALL hijsft(jp,jt,jout,ierror)
599  IF(ierror.NE.0) THEN
600  IF(ihpr2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
601  go to 50
602  ENDIF
603 C
604 C ********conduct soft scattering between JP and JT
605  jatt=jatt+jout
606 
607 200 CONTINUE
608 c
609 c**************************
610 c
611  DO 201 jp=1,ihnt2(1)
612  IF(nfp(jp,5).GT.2) THEN
613  ninp=ninp+1
614  ELSE IF(nfp(jp,5).EQ.2.OR.nfp(jp,5).EQ.1) THEN
615  nelp=nelp+1
616  ENDIF
617  201 continue
618  DO 202 jt=1,ihnt2(3)
619  IF(nft(jt,5).GT.2) THEN
620  nint=nint+1
621  ELSE IF(nft(jt,5).EQ.2.OR.nft(jt,5).EQ.1) THEN
622  nelt=nelt+1
623  ENDIF
624  202 continue
625 c
626 c*******************************
627 
628 
629 C********perform jet quenching for jets with PT>HIPR1(11)**********
630 
631  IF((ihpr2(8).NE.0.OR.ihpr2(3).NE.0).AND.ihpr2(4).GT.0.AND.
632  & ihnt2(1).GT.1.AND.ihnt2(3).GT.1) THEN
633  DO 271 i=1,ihnt2(1)
634  IF(nfp(i,7).EQ.1) CALL quench(i,1)
635 271 CONTINUE
636  DO 272 i=1,ihnt2(3)
637  IF(nft(i,7).EQ.1) CALL quench(i,2)
638 272 CONTINUE
639  DO 273 isg=1,nsg
640  IF(iasg(isg,3).EQ.1) CALL quench(isg,3)
641 273 CONTINUE
642  ENDIF
643 C
644 C**************fragment all the string systems in the following*****
645 C
646 C********N_ST is where particle information starts
647 C********N_STR+1 is the number of strings in fragmentation
648 C********the number of strings before a line is stored in K(I,4)
649 C********IDSTR is id number of the string system (91,92 or 93)
650 C
651  IF(ihpr2(20).NE.0) THEN
652  DO 360 isg=1,nsg
653  CALL hijfrg(isg,3,ierror)
654  IF(mstu(24).NE.0 .OR.ierror.GT.0) THEN
655  mstu(24)=0
656  mstu(28)=0
657  IF(ihpr2(10).NE.0) THEN
658  call lulist(1)
659  WRITE(6,*) 'error occured, repeat the event'
660  ENDIF
661  go to 50
662  ENDIF
663 C ********Check errors
664 C
665  n_st=1
666  idstr=92
667  IF(ihpr2(21).EQ.0) THEN
668  CALL luedit(2)
669  ELSE
670 351 n_st=n_st+1
671  IF(k(n_st,2).LT.91.OR.k(n_st,2).GT.93) go to 351
672  idstr=k(n_st,2)
673  n_st=n_st+1
674  ENDIF
675 C
676  IF(frame.EQ.'LAB') THEN
677  CALL hiboost
678  ENDIF
679 C ******** boost back to lab frame(if it was in)
680 C
681  n_str=0
682  DO 360 i=n_st,n
683  IF(k(i,2).EQ.idstr) THEN
684  n_str=n_str+1
685  go to 360
686  ENDIF
687  k(i,4)=n_str
688  natt=natt+1
689  katt(natt,1)=k(i,2)
690  katt(natt,2)=20
691  katt(natt,4)=k(i,1)
692  IF(k(i,3).EQ.0 .OR. k(k(i,3),2).EQ.idstr) THEN
693  katt(natt,3)=0
694  ELSE
695  katt(natt,3)=natt-i+k(i,3)+n_str-k(k(i,3),4)
696  ENDIF
697 C ****** identify the mother particle
698  patt(natt,1)=p(i,1)
699  patt(natt,2)=p(i,2)
700  patt(natt,3)=p(i,3)
701  patt(natt,4)=p(i,4)
702  eatt=eatt+p(i,4)
703 360 CONTINUE
704 C ********Fragment the q-qbar jets systems *****
705 C
706  jtp(1)=ihnt2(1)
707  jtp(2)=ihnt2(3)
708  DO 400 ntp=1,2
709  DO 400 j_jtp=1,jtp(ntp)
710  CALL hijfrg(j_jtp,ntp,ierror)
711  IF(mstu(24).NE.0 .OR. ierror.GT.0) THEN
712  mstu(24)=0
713  mstu(28)=0
714  IF(ihpr2(10).NE.0) THEN
715  call lulist(1)
716  WRITE(6,*) 'error occured, repeat the event'
717  ENDIF
718  go to 50
719  ENDIF
720 C ********check errors
721 C
722  n_st=1
723  idstr=92
724  IF(ihpr2(21).EQ.0) THEN
725  CALL luedit(2)
726  ELSE
727 381 n_st=n_st+1
728  IF(k(n_st,2).LT.91.OR.k(n_st,2).GT.93) go to 381
729  idstr=k(n_st,2)
730  n_st=n_st+1
731  ENDIF
732  IF(frame.EQ.'LAB') THEN
733  CALL hiboost
734  ENDIF
735 C ******** boost back to lab frame(if it was in)
736 C
737  nftp=nfp(j_jtp,5)
738  IF(ntp.EQ.2) nftp=10+nft(j_jtp,5)
739  n_str=0
740  DO 390 i=n_st,n
741  IF(k(i,2).EQ.idstr) THEN
742  n_str=n_str+1
743  go to 390
744  ENDIF
745  k(i,4)=n_str
746  natt=natt+1
747  katt(natt,1)=k(i,2)
748  katt(natt,2)=nftp
749  katt(natt,4)=k(i,1)
750  IF(k(i,3).EQ.0 .OR. k(k(i,3),2).EQ.idstr) THEN
751  katt(natt,3)=0
752  ELSE
753  katt(natt,3)=natt-i+k(i,3)+n_str-k(k(i,3),4)
754  ENDIF
755 C ****** identify the mother particle
756  patt(natt,1)=p(i,1)
757  patt(natt,2)=p(i,2)
758  patt(natt,3)=p(i,3)
759  patt(natt,4)=p(i,4)
760  eatt=eatt+p(i,4)
761 390 CONTINUE
762 400 CONTINUE
763 C ********Fragment the q-qq related string systems
764  ENDIF
765 
766  DO 450 i=1,ndr
767  natt=natt+1
768  katt(natt,1)=kfdr(i)
769  katt(natt,2)=40
770  katt(natt,3)=0
771  patt(natt,1)=pdr(i,1)
772  patt(natt,2)=pdr(i,2)
773  patt(natt,3)=pdr(i,3)
774  patt(natt,4)=pdr(i,4)
775  eatt=eatt+pdr(i,4)
776 450 CONTINUE
777 C ********store the direct-produced particles
778 C
779  dengy=eatt/(ihnt2(1)*hint1(6)+ihnt2(3)*hint1(7))-1.0
780  IF(abs(dengy).GT.hipr1(43).AND.ihpr2(20).NE.0
781  & .AND.ihpr2(21).EQ.0) THEN
782  IF(ihpr2(10).NE.0) THEN
783  WRITE(6,*) 'Energy not conserved, repeat the event'
784  ENDIF
785 C call lulist(1)
786  go to 50
787  ENDIF
788  RETURN
789  END
790 C
791 C
792 C
793 c SUBROUTINE HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT) !khaled
794 c SUBROUTINE HIJSET(EFRM,IAP,IZP,IAT,IZT) !khaled
795  SUBROUTINE hijset(efrm)
796 
797 c real*4 efrm
798  CHARACTER frame*4,proj*4,targ*4,eframe*4
799  DOUBLE PRECISION dd1,dd2,dd3,dd4
800  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
801  SAVE /histrng/
802  common/hijcrdn/yp(3,300),yt(3,300)
803  SAVE /hijcrdn/
804  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
805  SAVE /hiparnt/
806  common/hijdat/hidat0(10,10),hidat(10)
807  SAVE /hijdat/
808  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
809  SAVE /ludat1/
810  EXTERNAL fnkick,fnkick2,fnstru,fnstrum,fnstrus
811  CALL title
812 c ---------->khaled
813  frame='CMS'
814 c FRAME='LAB'
815  proj="A"
816  targ="B"
817 c write(*,*)"Fortran khaled"
818 c write(*,*)EFRM, IHNT2(1),IHNT2(2),IHNT2(3),IHNT2(4)
819 c write(*,*)IHNT2(5), IHNT2(6), HINT1(8), HINT1(9)
820 c write(*,*)iap, izp, iat,izt
821 c write(*,*)hipr1(1), hipr1(2),hipr1(3),hipr1(4), ihpr2(12)
822 c ---------> closed by khaled-------------->
823 c ------------------------------------
824 c IHNT2(1)=IAP
825 c IHNT2(2)=IZP
826 c IHNT2(3)=IAT
827 c IHNT2(4)=IZT
828 c IHNT2(5)=0
829 c IHNT2(6)=0
830 C
831 
832 
833 c HINT1(8)=MAX(ULMASS(2112),ULMASS(2212))
834 c HINT1(9)=HINT1(8)
835 C
836 c ---------> closed by khaled
837 c IF(PROJ.NE.'A') THEN
838 c IF(PROJ.EQ.'P') THEN
839 c IHNT2(5)=2212
840 c ELSE IF(PROJ.EQ.'PBAR') THEN
841 c IHNT2(5)=-2212
842 c ELSE IF(PROJ.EQ.'PI+') THEN
843 c IHNT2(5)=211
844 c ELSE IF(PROJ.EQ.'PI-') THEN
845 c IHNT2(5)=-211
846 c ELSE IF(PROJ.EQ.'K+') THEN
847 c IHNT2(5)=321
848 c ELSE IF(PROJ.EQ.'K-') THEN
849 c IHNT2(5)=-321
850 c ELSE IF(PROJ.EQ.'N') THEN
851 c IHNT2(5)=2112
852 c ELSE IF(PROJ.EQ.'NBAR') THEN
853 c IHNT2(5)=-2112
854 c ELSE
855 c WRITE(6,*) PROJ, 'wrong or unavailable proj name'
856 c STOP
857 c ENDIF
858 c HINT1(8)=ULMASS(IHNT2(5))
859 c ENDIF
860 c IF(TARG.NE.'A') THEN
861 c IF(TARG.EQ.'P') THEN
862 c IHNT2(6)=2212
863 c ELSE IF(TARG.EQ.'PBAR') THEN
864 c IHNT2(6)=-2212
865 c ELSE IF(TARG.EQ.'PI+') THEN
866 c IHNT2(6)=211
867 c ELSE IF(TARG.EQ.'PI-') THEN
868 c IHNT2(6)=-211
869 c ELSE IF(TARG.EQ.'K+') THEN
870 c IHNT2(6)=321
871 c ELSE IF(TARG.EQ.'K-') THEN
872 c IHNT2(6)=-321
873 c ELSE IF(TARG.EQ.'N') THEN
874 c IHNT2(6)=2112
875 c ELSE IF(TARG.EQ.'NBAR') THEN
876 c IHNT2(6)=-2112
877 c ELSE
878 c WRITE(6,*) TARG,'wrong or unavailable targ name'
879 c STOP
880 c ENDIF
881 c HINT1(9)=ULMASS(IHNT2(6))
882 c ENDIF
883 c ------------------------>end close by khaled
884 
885 C...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
886  IF(ihpr2(12).GT.0) THEN
887  CALL lugive('MDCY(C111,1)=0')
888  CALL lugive('MDCY(C310,1)=0')
889  CALL lugive('MDCY(C411,1)=0;MDCY(C-411,1)=0')
890  CALL lugive('MDCY(C421,1)=0;MDCY(C-421,1)=0')
891  CALL lugive('MDCY(C431,1)=0;MDCY(C-431,1)=0')
892  CALL lugive('MDCY(C511,1)=0;MDCY(C-511,1)=0')
893  CALL lugive('MDCY(C521,1)=0;MDCY(C-521,1)=0')
894  CALL lugive('MDCY(C531,1)=0;MDCY(C-531,1)=0')
895  CALL lugive('MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
896  CALL lugive('MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
897  CALL lugive('MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
898  CALL lugive('MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
899  CALL lugive('MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
900  CALL lugive('MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
901  CALL lugive('MDCY(C3334,1)=0;MDCY(C-3334,1)=0')
902  ENDIF
903  mstu(12)=0
904  mstu(21)=1
905  IF(ihpr2(10).EQ.0) THEN
906  mstu(22)=0
907  mstu(25)=0
908  mstu(26)=0
909  ENDIF
910  mstj(12)=ihpr2(11)
911  parj(21)=hipr1(2)
912  parj(41)=hipr1(3)
913  parj(42)=hipr1(4)
914 C ******** set up for jetset
915  IF(frame.EQ.'LAB') THEN
916  dd1=efrm
917  dd2=hint1(8)
918  dd3=hint1(9)
919  hint1(1)=sqrt(hint1(8)**2+2.0*hint1(9)*efrm+hint1(9)**2)
920  dd4=dsqrt(dd1**2-dd2**2)/(dd1+dd3)
921  hint1(2)=dd4
922  hint1(3)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
923  dd4=dsqrt(dd1**2-dd2**2)/dd1
924  hint1(4)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
925  hint1(5)=0.0
926  hint1(6)=efrm
927  hint1(7)=hint1(9)
928  ELSE IF(frame.EQ.'CMS') THEN
929  hint1(1)=efrm
930  hint1(2)=0.0
931  hint1(3)=0.0
932  dd1=hint1(1)
933  dd2=hint1(8)
934  dd3=hint1(9)
935  dd4=dsqrt(1.d0-4.d0*dd2**2/dd1**2)
936  hint1(4)=0.5*dlog((1.d0+dd4)/(1.d0-dd4))
937  dd4=dsqrt(1.d0-4.d0*dd3**2/dd1**2)
938  hint1(5)=-0.5*dlog((1.d0+dd4)/(1.d0-dd4))
939  hint1(6)=hint1(1)/2.0
940  hint1(7)=hint1(1)/2.0
941  ENDIF
942 C ********define Lorentz transform to lab frame
943 c
944 C ********calculate the cross sections involved with
945 C nucleon collisions.
946  IF(ihnt2(1).GT.1) THEN
947  CALL hijwds(ihnt2(1),1,rmax)
948  hipr1(34)=rmax
949 C ********set up Wood-Sax distr for proj.
950  ENDIF
951  IF(ihnt2(3).GT.1) THEN
952  CALL hijwds(ihnt2(3),2,rmax)
953  hipr1(35)=rmax
954 C ********set up Wood-Sax distr for targ.
955  ENDIF
956 c write(*,*)"radii hipr1(34--35)", hipr1(34),"<------>", hipr1(35) !khaled
957 c write(*,*)"hint1(6)--Hint1(7)", hint1(6),"<------>", hint1(7) !khaled
958 C
959 C
960  i=0
961 20 i=i+1
962  IF(i.EQ.10) go to 30
963  IF(hidat0(10,i).LE.hint1(1)) go to 20
964 30 IF(i.EQ.1) i=2
965  DO 40 j=1,9
966  hidat(j)=hidat0(j,i-1)+(hidat0(j,i)-hidat0(j,i-1))
967  & *(hint1(1)-hidat0(10,i-1))/(hidat0(10,i)-hidat0(10,i-1))
968 40 CONTINUE
969  hipr1(31)=hidat(5)
970  hipr1(30)=2.0*hidat(5)
971 c write(*,*)"hipr1(31--30)", hipr1(31),"<------>", hipr1(30) !khaled
972 C
973 C
974 
975  CALL hijcrs
976 
977 C
978  IF(ihpr2(5).NE.0) THEN
979  CALL hifun(3,0.0,36.0,fnkick)
980 C ********booking for generating pt**2 for pt kick
981  ENDIF
982  CALL hifun(7,0.0,6.0,fnkick2)
983  CALL hifun(4,0.0,1.0,fnstru)
984  CALL hifun(5,0.0,1.0,fnstrum)
985  CALL hifun(6,0.0,1.0,fnstrus)
986 C ********booking for x distribution of valence quarks
987  eframe='Ecm'
988  IF(frame.EQ.'LAB') eframe='Elab'
989  WRITE(6,100) eframe,efrm,proj,ihnt2(1),ihnt2(2),
990  & targ,ihnt2(3),ihnt2(4)
991 100 FORMAT(//10x,'****************************************
992  & **********'/
993  & 10x,'*',48x,'*'/
994  & 10x,'* HIJING has been initialized at *'/
995  & 10x,'*',13x,a4,'= ',f10.2,' GeV/n',13x,'*'/
996  & 10x,'*',48x,'*'/
997  & 10x,'*',8x,'for ',
998  & a4,'(',i3,',',i3,')',' + ',a4,'(',i3,',',i3,')',7x,'*'/
999  & 10x,'**************************************************')
1000  RETURN
1001  END
1002 C
1003 C
1004 C
1005  FUNCTION fnkick(X)
1006  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1007  SAVE /hiparnt/
1008  fnkick=1.0/(x+hipr1(19)**2)/(x+hipr1(20)**2)
1009  & /(1+exp((sqrt(x)-hipr1(20))/0.4))
1010  RETURN
1011  END
1012 C
1013 C
1014  FUNCTION fnkick2(X)
1015  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1016  SAVE /hiparnt/
1017  fnkick2=x*exp(-2.0*x/hipr1(42))
1018  RETURN
1019  END
1020 C
1021 C
1022 C
1023  FUNCTION fnstru(X)
1024  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1025  SAVE /hiparnt/
1026  fnstru=(1.0-x)**hipr1(44)/
1027  & (x**2+hipr1(45)**2/hint1(1)**2)**hipr1(46)
1028  RETURN
1029  END
1030 C
1031 C
1032 C
1033  FUNCTION fnstrum(X)
1034  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1035  SAVE /hiparnt/
1036  fnstrum=1.0/((1.0-x)**2+hipr1(45)**2/hint1(1)**2)**hipr1(46)
1037  & /(x**2+hipr1(45)**2/hint1(1)**2)**hipr1(46)
1038  RETURN
1039  END
1040 C
1041 C
1042  FUNCTION fnstrus(X)
1043  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1044  SAVE /hiparnt/
1045  fnstrus=(1.0-x)**hipr1(47)/
1046  & (x**2+hipr1(45)**2/hint1(1)**2)**hipr1(48)
1047  RETURN
1048  END
1049 C
1050 C
1051 C
1052 C
1053  SUBROUTINE hiboost
1054  IMPLICIT DOUBLE PRECISION(d)
1055  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
1056  SAVE /lujets/
1057  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
1058  SAVE /ludat1/
1059  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1060  SAVE /hiparnt/
1061  DO 100 i=1,n
1062  dbeta=p(i,3)/p(i,4)
1063  IF(abs(dbeta).GE.1.d0) THEN
1064  db=hint1(2)
1065  IF(db.GT.0.99999999d0) THEN
1066 C ********Rescale boost vector if too close to unity.
1067  WRITE(6,*) '(HIBOOT:) boost vector too large'
1068  db=0.99999999d0
1069  ENDIF
1070  dga=1d0/sqrt(1d0-db**2)
1071  dp3=p(i,3)
1072  dp4=p(i,4)
1073  p(i,3)=(dp3+db*dp4)*dga
1074  p(i,4)=(dp4+db*dp3)*dga
1075  go to 100
1076  ENDIF
1077  y=0.5*dlog((1.d0+dbeta)/(1.d0-dbeta))
1078  amt=sqrt(p(i,1)**2+p(i,2)**2+p(i,5)**2)
1079  p(i,3)=amt*sinh(y+hint1(3))
1080  p(i,4)=amt*cosh(y+hint1(3))
1081 100 CONTINUE
1082  RETURN
1083  END
1084 C
1085 C
1086 C
1087 C
1088  SUBROUTINE quench(JPJT,NTP)
1089  dimension rdp(300),lqp(300),rdt(300),lqt(300)
1090  common/hijcrdn/yp(3,300),yt(3,300)
1091  SAVE /hijcrdn/
1092  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1093  SAVE /hiparnt/
1094 C
1095  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
1096  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
1097  & pjpm(300,500),ntj(300),kftj(300,500),
1098  & pjtx(300,500),pjty(300,500),pjtz(300,500),
1099  & pjte(300,500),pjtm(300,500)
1100  SAVE /hijjet1/
1101  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
1102  & k2sg(900,100),pxsg(900,100),pysg(900,100),
1103  & pzsg(900,100),pesg(900,100),pmsg(900,100)
1104  SAVE /hijjet2/
1105  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
1106  SAVE /histrng/
1107  common/ranseed/nseed
1108  SAVE /ranseed/
1109 C
1110  bb=hint1(19) ! Uzhi
1111  phi=hint1(20) ! Uzhi
1112  bbx=bb*cos(phi) ! Uzhi
1113  bby=bb*sin(phi) ! Uzhi
1114 c
1115  IF(ntp.EQ.2) go to 400
1116  IF(ntp.EQ.3) go to 2000
1117 C*******************************************************
1118 C Jet interaction for proj jet in the direction PHIP
1119 C******************************************************
1120 C
1121  IF(nfp(jpjt,7).NE.1) RETURN
1122 
1123  jp=jpjt
1124  DO 290 i=1,npj(jp)
1125  ptjet0=sqrt(pjpx(jp,i)**2+pjpy(jp,i)**2)
1126  IF(ptjet0.LE.hipr1(11)) go to 290
1127  ptot=sqrt(ptjet0*ptjet0+pjpz(jp,i)**2)
1128  IF(ptot.LT.hipr1(8)) go to 290
1129  phip=ulangl(pjpx(jp,i),pjpy(jp,i))
1130 C******* find the wounded proj which can interact with jet***
1131  kp=0
1132  DO 100 i2=1,ihnt2(1)
1133  IF(nfp(i2,5).NE.3 .OR. i2.EQ.jp) go to 100
1134  dx=yp(1,i2)-yp(1,jp)
1135  dy=yp(2,i2)-yp(2,jp)
1136  phi=ulangl(dx,dy)
1137  dphi=abs(phi-phip)
1138  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
1139  IF(dphi.GE.hipr1(40)/2.0) go to 100
1140  rd0=sqrt(dx*dx+dy*dy)
1141  IF(rd0*sin(dphi).GT.hipr1(12)) go to 100
1142  kp=kp+1
1143  lqp(kp)=i2
1144  rdp(kp)=cos(dphi)*rd0
1145  100 CONTINUE
1146 C******* rearrange according decending rd************
1147  DO 110 i2=1,kp-1
1148  DO 110 j2=i2+1,kp
1149  IF(rdp(i2).LT.rdp(j2)) go to 110
1150  rd=rdp(i2)
1151  lq=lqp(i2)
1152  rdp(i2)=rdp(j2)
1153  lqp(i2)=lqp(j2)
1154  rdp(j2)=rd
1155  lqp(j2)=lq
1156  110 CONTINUE
1157 C****** find wounded targ which can interact with jet********
1158  kt=0
1159  DO 120 i2=1,ihnt2(3)
1160  IF(nft(i2,5).NE.3) go to 120
1161  dx=yt(1,i2)-yp(1,jp)-bbx
1162  dy=yt(2,i2)-yp(2,jp)-bby
1163  phi=ulangl(dx,dy)
1164  dphi=abs(phi-phip)
1165  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
1166  IF(dphi.GT.hipr1(40)/2.0) go to 120
1167  rd0=sqrt(dx*dx+dy*dy)
1168  IF(rd0*sin(dphi).GT.hipr1(12)) go to 120
1169  kt=kt+1
1170  lqt(kt)=i2
1171  rdt(kt)=cos(dphi)*rd0
1172  120 CONTINUE
1173 C******* rearrange according decending rd************
1174  DO 130 i2=1,kt-1
1175  DO 130 j2=i2+1,kt
1176  IF(rdt(i2).LT.rdt(j2)) go to 130
1177  rd=rdt(i2)
1178  lq=lqt(i2)
1179  rdt(i2)=rdt(j2)
1180  lqt(i2)=lqt(j2)
1181  rdt(j2)=rd
1182  lqt(j2)=lq
1183  130 CONTINUE
1184 
1185  mp=0
1186  mt=0
1187  r0=0.0
1188  nq=0
1189  dp=0.0
1190  ptot=sqrt(pjpx(jp,i)**2+pjpy(jp,i)**2+pjpz(jp,i)**2)
1191  v1=pjpx(jp,i)/ptot
1192  v2=pjpy(jp,i)/ptot
1193  v3=pjpz(jp,i)/ptot
1194 
1195  200 rn=rlu(0)
1196  210 IF(mt.GE.kt .AND. mp.GE.kp) go to 290
1197  IF(mt.GE.kt) go to 220
1198  IF(mp.GE.kp) go to 240
1199  IF(rdp(mp+1).GT.rdt(mt+1)) go to 240
1200  220 mp=mp+1
1201  drr=rdp(mp)-r0
1202  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 210
1203  dp=drr*hipr1(14)
1204  IF(kfpj(jp,i).NE.21) dp=0.5*dp
1205 C ********string tension of quark jet is 0.5 of gluon's
1206  IF(dp.LE.0.2) go to 210
1207  IF(ptot.LE.0.4) go to 290
1208  IF(ptot.LE.dp) dp=ptot-0.2
1209  de=dp
1210 
1211  IF(kfpj(jp,i).NE.21) THEN
1212  prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
1213  & +pp(lqp(mp),3)**2
1214  de=sqrt(pjpm(jp,i)**2+ptot**2)
1215  & -sqrt(pjpm(jp,i)**2+(ptot-dp)**2)
1216  ershu=(pp(lqp(mp),4)+de-dp)**2
1217  amshu=ershu-prshu
1218  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 210
1219  pp(lqp(mp),4)=sqrt(ershu)
1220  pp(lqp(mp),5)=sqrt(amshu)
1221  ENDIF
1222 C ********reshuffle the energy when jet has mass
1223  r0=rdp(mp)
1224  dp1=dp*v1
1225  dp2=dp*v2
1226  dp3=dp*v3
1227 C ********momentum and energy transfer from jet
1228 
1229  npj(lqp(mp))=npj(lqp(mp))+1
1230  kfpj(lqp(mp),npj(lqp(mp)))=21
1231  pjpx(lqp(mp),npj(lqp(mp)))=dp1
1232  pjpy(lqp(mp),npj(lqp(mp)))=dp2
1233  pjpz(lqp(mp),npj(lqp(mp)))=dp3
1234  pjpe(lqp(mp),npj(lqp(mp)))=dp
1235  pjpm(lqp(mp),npj(lqp(mp)))=0.0
1236  go to 260
1237 
1238  240 mt=mt+1
1239  drr=rdt(mt)-r0
1240  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 210
1241  dp=drr*hipr1(14)
1242  IF(dp.LE.0.2) go to 210
1243  IF(ptot.LE.0.4) go to 290
1244  IF(ptot.LE.dp) dp=ptot-0.2
1245  de=dp
1246 
1247  IF(kfpj(jp,i).NE.21) THEN
1248  prshu=pt(lqt(mt),1)**2+pt(lqt(mt),2)**2
1249  & +pt(lqt(mt),3)**2
1250  de=sqrt(pjpm(jp,i)**2+ptot**2)
1251  & -sqrt(pjpm(jp,i)**2+(ptot-dp)**2)
1252  ershu=(pt(lqt(mt),4)+de-dp)**2
1253  amshu=ershu-prshu
1254  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 210
1255  pt(lqt(mt),4)=sqrt(ershu)
1256  pt(lqt(mt),5)=sqrt(amshu)
1257  ENDIF
1258 C ********reshuffle the energy when jet has mass
1259 
1260  r0=rdt(mt)
1261  dp1=dp*v1
1262  dp2=dp*v2
1263  dp3=dp*v3
1264 C ********momentum and energy transfer from jet
1265  ntj(lqt(mt))=ntj(lqt(mt))+1
1266  kftj(lqt(mt),ntj(lqt(mt)))=21
1267  pjtx(lqt(mt),ntj(lqt(mt)))=dp1
1268  pjty(lqt(mt),ntj(lqt(mt)))=dp2
1269  pjtz(lqt(mt),ntj(lqt(mt)))=dp3
1270  pjte(lqt(mt),ntj(lqt(mt)))=dp
1271  pjtm(lqt(mt),ntj(lqt(mt)))=0.0
1272 
1273  260 pjpx(jp,i)=(ptot-dp)*v1
1274  pjpy(jp,i)=(ptot-dp)*v2
1275  pjpz(jp,i)=(ptot-dp)*v3
1276  pjpe(jp,i)=pjpe(jp,i)-de
1277 
1278  ptot=ptot-dp
1279  nq=nq+1
1280  go to 200
1281  290 CONTINUE
1282 
1283  RETURN
1284 
1285 C*******************************************************
1286 C Jet interaction for target jet in the direction PHIT
1287 C******************************************************
1288 C
1289 C******* find the wounded proj which can interact with jet***
1290 
1291  400 IF(nft(jpjt,7).NE.1) RETURN
1292  jt=jpjt
1293  DO 690 i=1,ntj(jt)
1294  ptjet0=sqrt(pjtx(jt,i)**2+pjty(jt,i)**2)
1295  IF(ptjet0.LE.hipr1(11)) go to 690
1296  ptot=sqrt(ptjet0*ptjet0+pjtz(jt,i)**2)
1297  IF(ptot.LT.hipr1(8)) go to 690
1298  phit=ulangl(pjtx(jt,i),pjty(jt,i))
1299  kp=0
1300  DO 500 i2=1,ihnt2(1)
1301  IF(nfp(i2,5).NE.3) go to 500
1302  dx=yp(1,i2)+bbx-yt(1,jt)
1303  dy=yp(2,i2)+bby-yt(2,jt)
1304  phi=ulangl(dx,dy)
1305  dphi=abs(phi-phit)
1306  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
1307  IF(dphi.GT.hipr1(40)/2.0) go to 500
1308  rd0=sqrt(dx*dx+dy*dy)
1309  IF(rd0*sin(dphi).GT.hipr1(12)) go to 500
1310  kp=kp+1
1311  lqp(kp)=i2
1312  rdp(kp)=cos(dphi)*rd0
1313  500 CONTINUE
1314 C******* rearrange according to decending rd************
1315  DO 510 i2=1,kp-1
1316  DO 510 j2=i2+1,kp
1317  IF(rdp(i2).LT.rdp(j2)) go to 510
1318  rd=rdp(i2)
1319  lq=lqp(i2)
1320  rdp(i2)=rdp(j2)
1321  lqp(i2)=lqp(j2)
1322  rdp(j2)=rd
1323  lqp(j2)=lq
1324  510 CONTINUE
1325 C****** find wounded targ which can interact with jet********
1326  kt=0
1327  DO 520 i2=1,ihnt2(3)
1328  IF(nft(i2,5).NE.3 .OR. i2.EQ.jt) go to 520
1329  dx=yt(1,i2)-yt(1,jt)
1330  dy=yt(2,i2)-yt(2,jt)
1331  phi=ulangl(dx,dy)
1332  dphi=abs(phi-phit)
1333  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
1334  IF(dphi.GT.hipr1(40)/2.0) go to 520
1335  rd0=sqrt(dx*dx+dy*dy)
1336  IF(rd0*sin(dphi).GT.hipr1(12)) go to 520
1337  kt=kt+1
1338  lqt(kt)=i2
1339  rdt(kt)=cos(dphi)*rd0
1340  520 CONTINUE
1341 C******* rearrange according to decending rd************
1342  DO 530 i2=1,kt-1
1343  DO 530 j2=i2+1,kt
1344  IF(rdt(i2).LT.rdt(j2)) go to 530
1345  rd=rdt(i2)
1346  lq=lqt(i2)
1347  rdt(i2)=rdt(j2)
1348  lqt(i2)=lqt(j2)
1349  rdt(j2)=rd
1350  lqt(j2)=lq
1351  530 CONTINUE
1352 
1353  mp=0
1354  mt=0
1355  nq=0
1356  dp=0.0
1357  r0=0.0
1358  ptot=sqrt(pjtx(jt,i)**2+pjty(jt,i)**2+pjtz(jt,i)**2)
1359  v1=pjtx(jt,i)/ptot
1360  v2=pjty(jt,i)/ptot
1361  v3=pjtz(jt,i)/ptot
1362 
1363  600 rn=rlu(0)
1364  610 IF(mt.GE.kt .AND. mp.GE.kp) go to 690
1365  IF(mt.GE.kt) go to 620
1366  IF(mp.GE.kp) go to 640
1367  IF(rdp(mp+1).GT.rdt(mt+1)) go to 640
1368 620 mp=mp+1
1369  drr=rdp(mp)-r0
1370  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 610
1371  dp=drr*hipr1(14)
1372  IF(kftj(jt,i).NE.21) dp=0.5*dp
1373 C ********string tension of quark jet is 0.5 of gluon's
1374  IF(dp.LE.0.2) go to 610
1375  IF(ptot.LE.0.4) go to 690
1376  IF(ptot.LE.dp) dp=ptot-0.2
1377  de=dp
1378 C
1379  IF(kftj(jt,i).NE.21) THEN
1380  prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
1381  & +pp(lqp(mp),3)**2
1382  de=sqrt(pjtm(jt,i)**2+ptot**2)
1383  & -sqrt(pjtm(jt,i)**2+(ptot-dp)**2)
1384  ershu=(pp(lqp(mp),4)+de-dp)**2
1385  amshu=ershu-prshu
1386  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 610
1387  pp(lqp(mp),4)=sqrt(ershu)
1388  pp(lqp(mp),5)=sqrt(amshu)
1389  ENDIF
1390 C ********reshuffle the energy when jet has mass
1391 C
1392  r0=rdp(mp)
1393  dp1=dp*v1
1394  dp2=dp*v2
1395  dp3=dp*v3
1396 C ********momentum and energy transfer from jet
1397  npj(lqp(mp))=npj(lqp(mp))+1
1398  kfpj(lqp(mp),npj(lqp(mp)))=21
1399  pjpx(lqp(mp),npj(lqp(mp)))=dp1
1400  pjpy(lqp(mp),npj(lqp(mp)))=dp2
1401  pjpz(lqp(mp),npj(lqp(mp)))=dp3
1402  pjpe(lqp(mp),npj(lqp(mp)))=dp
1403  pjpm(lqp(mp),npj(lqp(mp)))=0.0
1404 
1405  go to 660
1406 
1407 640 mt=mt+1
1408  drr=rdt(mt)-r0
1409  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 610
1410  dp=drr*hipr1(14)
1411  IF(dp.LE.0.2) go to 610
1412  IF(ptot.LE.0.4) go to 690
1413  IF(ptot.LE.dp) dp=ptot-0.2
1414  de=dp
1415 
1416  IF(kftj(jt,i).NE.21) THEN
1417  prshu=pt(lqt(mt),1)**2+pt(lqt(mt),2)**2
1418  & +pt(lqt(mt),3)**2
1419  de=sqrt(pjtm(jt,i)**2+ptot**2)
1420  & -sqrt(pjtm(jt,i)**2+(ptot-dp)**2)
1421  ershu=(pt(lqt(mt),4)+de-dp)**2
1422  amshu=ershu-prshu
1423  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 610
1424  pt(lqt(mt),4)=sqrt(ershu)
1425  pt(lqt(mt),5)=sqrt(amshu)
1426  ENDIF
1427 C ********reshuffle the energy when jet has mass
1428 
1429  r0=rdt(mt)
1430  dp1=dp*v1
1431  dp2=dp*v2
1432  dp3=dp*v3
1433 C ********momentum and energy transfer from jet
1434  ntj(lqt(mt))=ntj(lqt(mt))+1
1435  kftj(lqt(mt),ntj(lqt(mt)))=21
1436  pjtx(lqt(mt),ntj(lqt(mt)))=dp1
1437  pjty(lqt(mt),ntj(lqt(mt)))=dp2
1438  pjtz(lqt(mt),ntj(lqt(mt)))=dp3
1439  pjte(lqt(mt),ntj(lqt(mt)))=dp
1440  pjtm(lqt(mt),ntj(lqt(mt)))=0.0
1441 
1442 660 pjtx(jt,i)=(ptot-dp)*v1
1443  pjty(jt,i)=(ptot-dp)*v2
1444  pjtz(jt,i)=(ptot-dp)*v3
1445  pjte(jt,i)=pjte(jt,i)-de
1446 
1447  ptot=ptot-dp
1448  nq=nq+1
1449  go to 600
1450 690 CONTINUE
1451  RETURN
1452 C********************************************************
1453 C Q-QBAR jet interaction
1454 C********************************************************
1455 2000 isg=jpjt
1456  IF(iasg(isg,3).NE.1) RETURN
1457 C
1458  jp=iasg(isg,1)
1459  jt=iasg(isg,2)
1460  xj=(yp(1,jp)+bbx+yt(1,jt))/2.0
1461  yj=(yp(2,jp)+bby+yt(2,jt))/2.0
1462  DO 2690 i=1,njsg(isg)
1463  ptjet0=sqrt(pxsg(isg,i)**2+pysg(isg,i)**2)
1464  IF(ptjet0.LE.hipr1(11).OR.pesg(isg,i).LT.hipr1(1))
1465  & go to 2690
1466  ptot=sqrt(ptjet0*ptjet0+pzsg(isg,i)**2)
1467  IF(ptot.LT.max(hipr1(1),hipr1(8))) go to 2690
1468  phiq=ulangl(pxsg(isg,i),pysg(isg,i))
1469  kp=0
1470  DO 2500 i2=1,ihnt2(1)
1471  IF(nfp(i2,5).NE.3.OR.i2.EQ.jp) go to 2500
1472  dx=yp(1,i2)+bbx-xj
1473  dy=yp(2,i2)+bby-yj
1474  phi=ulangl(dx,dy)
1475  dphi=abs(phi-phiq)
1476  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
1477  IF(dphi.GT.hipr1(40)/2.0) go to 2500
1478  rd0=sqrt(dx*dx+dy*dy)
1479  IF(rd0*sin(dphi).GT.hipr1(12)) go to 2500
1480  kp=kp+1
1481  lqp(kp)=i2
1482  rdp(kp)=cos(dphi)*rd0
1483  2500 CONTINUE
1484 C******* rearrange according to decending rd************
1485  DO 2510 i2=1,kp-1
1486  DO 2510 j2=i2+1,kp
1487  IF(rdp(i2).LT.rdp(j2)) go to 2510
1488  rd=rdp(i2)
1489  lq=lqp(i2)
1490  rdp(i2)=rdp(j2)
1491  lqp(i2)=lqp(j2)
1492  rdp(j2)=rd
1493  lqp(j2)=lq
1494  2510 CONTINUE
1495 C****** find wounded targ which can interact with jet********
1496  kt=0
1497  DO 2520 i2=1,ihnt2(3)
1498  IF(nft(i2,5).NE.3 .OR. i2.EQ.jt) go to 2520
1499  dx=yt(1,i2)-xj
1500  dy=yt(2,i2)-yj
1501  phi=ulangl(dx,dy)
1502  dphi=abs(phi-phiq)
1503  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
1504  IF(dphi.GT.hipr1(40)/2.0) go to 2520
1505  rd0=sqrt(dx*dx+dy*dy)
1506  IF(rd0*sin(dphi).GT.hipr1(12)) go to 2520
1507  kt=kt+1
1508  lqt(kt)=i2
1509  rdt(kt)=cos(dphi)*rd0
1510  2520 CONTINUE
1511 C******* rearrange according to decending rd************
1512  DO 2530 i2=1,kt-1
1513  DO 2530 j2=i2+1,kt
1514  IF(rdt(i2).LT.rdt(j2)) go to 2530
1515  rd=rdt(i2)
1516  lq=lqt(i2)
1517  rdt(i2)=rdt(j2)
1518  lqt(i2)=lqt(j2)
1519  rdt(j2)=rd
1520  lqt(j2)=lq
1521  2530 CONTINUE
1522 
1523  mp=0
1524  mt=0
1525  nq=0
1526  dp=0.0
1527  r0=0.0
1528  ptot=sqrt(pxsg(isg,i)**2+pysg(isg,i)**2
1529  & +pzsg(isg,i)**2)
1530  v1=pxsg(isg,i)/ptot
1531  v2=pysg(isg,i)/ptot
1532  v3=pzsg(isg,i)/ptot
1533 
1534  2600 rn=rlu(0)
1535  2610 IF(mt.GE.kt .AND. mp.GE.kp) go to 2690
1536  IF(mt.GE.kt) go to 2620
1537  IF(mp.GE.kp) go to 2640
1538  IF(rdp(mp+1).GT.rdt(mt+1)) go to 2640
1539  2620 mp=mp+1
1540  drr=rdp(mp)-r0
1541  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 2610
1542  dp=drr*hipr1(14)/2.0
1543  IF(dp.LE.0.2) go to 2610
1544  IF(ptot.LE.0.4) go to 2690
1545  IF(ptot.LE.dp) dp=ptot-0.2
1546  de=dp
1547 C
1548  IF(k2sg(isg,i).NE.21) THEN
1549  IF(ptot.LT.dp+hipr1(1)) go to 2690
1550  prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
1551  & +pp(lqp(mp),3)**2
1552  de=sqrt(pmsg(isg,i)**2+ptot**2)
1553  & -sqrt(pmsg(isg,i)**2+(ptot-dp)**2)
1554  ershu=(pp(lqp(mp),4)+de-dp)**2
1555  amshu=ershu-prshu
1556  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 2610
1557  pp(lqp(mp),4)=sqrt(ershu)
1558  pp(lqp(mp),5)=sqrt(amshu)
1559  ENDIF
1560 C ********reshuffle the energy when jet has mass
1561 C
1562  r0=rdp(mp)
1563  dp1=dp*v1
1564  dp2=dp*v2
1565  dp3=dp*v3
1566 C ********momentum and energy transfer from jet
1567  npj(lqp(mp))=npj(lqp(mp))+1
1568  kfpj(lqp(mp),npj(lqp(mp)))=21
1569  pjpx(lqp(mp),npj(lqp(mp)))=dp1
1570  pjpy(lqp(mp),npj(lqp(mp)))=dp2
1571  pjpz(lqp(mp),npj(lqp(mp)))=dp3
1572  pjpe(lqp(mp),npj(lqp(mp)))=dp
1573  pjpm(lqp(mp),npj(lqp(mp)))=0.0
1574 
1575  go to 2660
1576 
1577  2640 mt=mt+1
1578  drr=rdt(mt)-r0
1579  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 2610
1580  dp=drr*hipr1(14)
1581  IF(dp.LE.0.2) go to 2610
1582  IF(ptot.LE.0.4) go to 2690
1583  IF(ptot.LE.dp) dp=ptot-0.2
1584  de=dp
1585 
1586  IF(k2sg(isg,i).NE.21) THEN
1587  IF(ptot.LT.dp+hipr1(1)) go to 2690
1588  prshu=pt(lqt(mt),1)**2+pt(lqt(mt),2)**2
1589  & +pt(lqt(mt),3)**2
1590  de=sqrt(pmsg(isg,i)**2+ptot**2)
1591  & -sqrt(pmsg(isg,i)**2+(ptot-dp)**2)
1592  ershu=(pt(lqt(mt),4)+de-dp)**2
1593  amshu=ershu-prshu
1594  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 2610
1595  pt(lqt(mt),4)=sqrt(ershu)
1596  pt(lqt(mt),5)=sqrt(amshu)
1597  ENDIF
1598 C ********reshuffle the energy when jet has mass
1599 
1600  r0=rdt(mt)
1601  dp1=dp*v1
1602  dp2=dp*v2
1603  dp3=dp*v3
1604 C ********momentum and energy transfer from jet
1605  ntj(lqt(mt))=ntj(lqt(mt))+1
1606  kftj(lqt(mt),ntj(lqt(mt)))=21
1607  pjtx(lqt(mt),ntj(lqt(mt)))=dp1
1608  pjty(lqt(mt),ntj(lqt(mt)))=dp2
1609  pjtz(lqt(mt),ntj(lqt(mt)))=dp3
1610  pjte(lqt(mt),ntj(lqt(mt)))=dp
1611  pjtm(lqt(mt),ntj(lqt(mt)))=0.0
1612 
1613  2660 pxsg(isg,i)=(ptot-dp)*v1
1614  pysg(isg,i)=(ptot-dp)*v2
1615  pzsg(isg,i)=(ptot-dp)*v3
1616  pesg(isg,i)=pesg(isg,i)-de
1617 
1618  ptot=ptot-dp
1619  nq=nq+1
1620  go to 2600
1621  2690 CONTINUE
1622  RETURN
1623  END
1624 
1625 C
1626 C
1627 C
1628 C
1629  SUBROUTINE hijfrg(JTP,NTP,IERROR)
1630 C NTP=1, fragment proj string, NTP=2, targ string,
1631 C NTP=3, independent
1632 C strings from jets. JTP is the line number of the string
1633 C*******Fragment all leading strings of proj and targ**************
1634 C IHNT2(1)=atomic #, IHNT2(2)=proton #(=-1 if anti-proton) *
1635 C******************************************************************
1636  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
1637  SAVE /hiparnt/
1638  common/hijdat/hidat0(10,10),hidat(10)
1639  SAVE /hijdat/
1640  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
1641  SAVE /histrng/
1642  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
1643  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
1644  & pjpm(300,500),ntj(300),kftj(300,500),
1645  & pjtx(300,500),pjty(300,500),pjtz(300,500),
1646  & pjte(300,500),pjtm(300,500)
1647  SAVE /hijjet1/
1648  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
1649  & k2sg(900,100),pxsg(900,100),pysg(900,100),
1650  & pzsg(900,100),pesg(900,100),pmsg(900,100)
1651  SAVE /hijjet2/
1652 C
1653  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
1654  SAVE /lujets/
1655  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
1656  SAVE /ludat1/
1657  common/ranseed/nseed
1658  SAVE /ranseed/
1659 
1660  ierror=0
1661  CALL luedit(0)
1662  n=0
1663 C ********initialize the document lines
1664  IF(ntp.EQ.3) THEN
1665  isg=jtp
1666  n=njsg(isg)
1667  DO 100 i=1,njsg(isg)
1668  k(i,1)=k1sg(isg,i)
1669  k(i,2)=k2sg(isg,i)
1670  p(i,1)=pxsg(isg,i)
1671  p(i,2)=pysg(isg,i)
1672  p(i,3)=pzsg(isg,i)
1673  p(i,4)=pesg(isg,i)
1674  p(i,5)=pmsg(isg,i)
1675 100 CONTINUE
1676 C IF(IHPR2(1).GT.0) CALL ATTRAD(IERROR)
1677 c IF(IERROR.NE.0) RETURN
1678 C CALL LULIST(1)
1679  CALL luexec
1680  RETURN
1681  ENDIF
1682 C
1683  IF(ntp.EQ.2) go to 200
1684  IF(jtp.GT.ihnt2(1)) RETURN
1685  IF(nfp(jtp,5).NE.3.AND.nfp(jtp,3).NE.0
1686  & .AND.npj(jtp).EQ.0.AND.nfp(jtp,10).EQ.0) go to 1000
1687  IF(nfp(jtp,15).EQ.-1) THEN
1688  kf1=nfp(jtp,2)
1689  kf2=nfp(jtp,1)
1690  pq21=pp(jtp,6)
1691  pq22=pp(jtp,7)
1692  pq11=pp(jtp,8)
1693  pq12=pp(jtp,9)
1694  am1=pp(jtp,15)
1695  am2=pp(jtp,14)
1696  ELSE
1697  kf1=nfp(jtp,1)
1698  kf2=nfp(jtp,2)
1699  pq21=pp(jtp,8)
1700  pq22=pp(jtp,9)
1701  pq11=pp(jtp,6)
1702  pq12=pp(jtp,7)
1703  am1=pp(jtp,14)
1704  am2=pp(jtp,15)
1705  ENDIF
1706 C ********for NFP(JTP,15)=-1 NFP(JTP,1) IS IN -Z DIRECTION
1707  pb1=pq11+pq21
1708  pb2=pq12+pq22
1709  pb3=pp(jtp,3)
1710  pecm=pp(jtp,5)
1711  btz=pb3/pp(jtp,4)
1712  IF((abs(pb1-pp(jtp,1)).GT.0.01.OR.
1713  & abs(pb2-pp(jtp,2)).GT.0.01).AND.ihpr2(10).NE.0)
1714  & WRITE(6,*) ' Pt of Q and QQ do not sum to the total'
1715 
1716  go to 300
1717 
1718 200 IF(jtp.GT.ihnt2(3)) RETURN
1719  IF(nft(jtp,5).NE.3.AND.nft(jtp,3).NE.0
1720  & .AND.ntj(jtp).EQ.0.AND.nft(jtp,10).EQ.0) go to 1200
1721  IF(nft(jtp,15).EQ.1) THEN
1722  kf1=nft(jtp,1)
1723  kf2=nft(jtp,2)
1724  pq11=pt(jtp,6)
1725  pq12=pt(jtp,7)
1726  pq21=pt(jtp,8)
1727  pq22=pt(jtp,9)
1728  am1=pt(jtp,14)
1729  am2=pt(jtp,15)
1730  ELSE
1731  kf1=nft(jtp,2)
1732  kf2=nft(jtp,1)
1733  pq11=pt(jtp,8)
1734  pq12=pt(jtp,9)
1735  pq21=pt(jtp,6)
1736  pq22=pt(jtp,7)
1737  am1=pt(jtp,15)
1738  am2=pt(jtp,14)
1739  ENDIF
1740 C ********for NFT(JTP,15)=1 NFT(JTP,1) IS IN +Z DIRECTION
1741  pb1=pq11+pq21
1742  pb2=pq12+pq22
1743  pb3=pt(jtp,3)
1744  pecm=pt(jtp,5)
1745  btz=pb3/pt(jtp,4)
1746 
1747  IF((abs(pb1-pt(jtp,1)).GT.0.01.OR.
1748  & abs(pb2-pt(jtp,2)).GT.0.01).AND.ihpr2(10).NE.0)
1749  & WRITE(6,*) ' Pt of Q and QQ do not sum to the total'
1750 
1751 300 IF(pecm.LT.hipr1(1)) THEN
1752  ierror=1
1753  IF(ihpr2(10).EQ.0) RETURN
1754  WRITE(6,*) ' ECM=',pecm,' energy of the string is too small'
1755  RETURN
1756  ENDIF
1757  amt=pecm**2+pb1**2+pb2**2
1758  amt1=am1**2+pq11**2+pq12**2
1759  amt2=am2**2+pq21**2+pq22**2
1760  pzcm=sqrt(abs(amt**2+amt1**2+amt2**2-2.0*amt*amt1
1761  & -2.0*amt*amt2-2.0*amt1*amt2))/2.0/sqrt(amt)
1762 C *******PZ of end-partons in c.m. frame of the string
1763  k(1,1)=2
1764  k(1,2)=kf1
1765  p(1,1)=pq11
1766  p(1,2)=pq12
1767  p(1,3)=pzcm
1768  p(1,4)=sqrt(amt1+pzcm**2)
1769  p(1,5)=am1
1770  k(2,1)=1
1771  k(2,2)=kf2
1772  p(2,1)=pq21
1773  p(2,2)=pq22
1774  p(2,3)=-pzcm
1775  p(2,4)=sqrt(amt2+pzcm**2)
1776  p(2,5)=am2
1777  n=2
1778 C*****
1779  CALL hirobo(0.0,0.0,0.0,0.0,btz)
1780  jetot=0
1781  IF((pq21**2+pq22**2).GT.(pq11**2+pq12**2)) THEN
1782  pmax1=p(2,1)
1783  pmax2=p(2,2)
1784  pmax3=p(2,3)
1785  ELSE
1786  pmax1=p(1,1)
1787  pmax2=p(1,2)
1788  pmax3=p(1,3)
1789  ENDIF
1790  IF(ntp.EQ.1) THEN
1791  pp(jtp,10)=pmax1
1792  pp(jtp,11)=pmax2
1793  pp(jtp,12)=pmax3
1794  ELSE IF(ntp.EQ.2) THEN
1795  pt(jtp,10)=pmax1
1796  pt(jtp,11)=pmax2
1797  pt(jtp,12)=pmax3
1798  ENDIF
1799 C*******************attach produced jets to the leading partons****
1800  IF(ntp.EQ.1.AND.npj(jtp).NE.0) THEN
1801  jetot=npj(jtp)
1802 C IF(NPJ(JTP).GE.2) CALL HIJSRT(JTP,1)
1803 C ********sort jets in order of y
1804  iex=0
1805  IF((abs(kf1).GT.1000.AND.kf1.LT.0)
1806  & .OR.(abs(kf1).LT.1000.AND.kf1.GT.0)) iex=1
1807  DO 520 i=n,2,-1
1808  DO 520 j=1,5
1809  ii=npj(jtp)+i
1810  k(ii,j)=k(i,j)
1811  p(ii,j)=p(i,j)
1812  v(ii,j)=v(i,j)
1813 520 CONTINUE
1814  DO 540 i=1,npj(jtp)
1815  DO 542 j=1,5
1816  k(i+1,j)=0
1817  v(i+1,j)=0
1818 542 CONTINUE
1819  i0=i
1820  IF(iex.EQ.1) i0=npj(jtp)-i+1
1821 C ********reverse the order of jets
1822  kk1=kfpj(jtp,i0)
1823  k(i+1,1)=2
1824  k(i+1,2)=kk1
1825  IF(kk1.NE.21 .AND. kk1.NE.0) k(i+1,1)=
1826  & 1+(abs(kk1)+(2*iex-1)*kk1)/2/abs(kk1)
1827  p(i+1,1)=pjpx(jtp,i0)
1828  p(i+1,2)=pjpy(jtp,i0)
1829  p(i+1,3)=pjpz(jtp,i0)
1830  p(i+1,4)=pjpe(jtp,i0)
1831  p(i+1,5)=pjpm(jtp,i0)
1832 540 CONTINUE
1833  n=n+npj(jtp)
1834  ELSE IF(ntp.EQ.2.AND.ntj(jtp).NE.0) THEN
1835  jetot=ntj(jtp)
1836 c IF(NTJ(JTP).GE.2) CALL HIJSRT(JTP,2)
1837 C ********sort jets in order of y
1838  iex=1
1839  IF((abs(kf2).GT.1000.AND.kf2.LT.0)
1840  & .OR.(abs(kf2).LT.1000.AND.kf2.GT.0)) iex=0
1841  DO 560 i=n,2,-1
1842  DO 560 j=1,5
1843  ii=ntj(jtp)+i
1844  k(ii,j)=k(i,j)
1845  p(ii,j)=p(i,j)
1846  v(ii,j)=v(i,j)
1847 560 CONTINUE
1848  DO 580 i=1,ntj(jtp)
1849  DO 582 j=1,5
1850  k(i+1,j)=0
1851  v(i+1,j)=0
1852 582 CONTINUE
1853  i0=i
1854  IF(iex.EQ.1) i0=ntj(jtp)-i+1
1855 C ********reverse the order of jets
1856  kk1=kftj(jtp,i0)
1857  k(i+1,1)=2
1858  k(i+1,2)=kk1
1859  IF(kk1.NE.21 .AND. kk1.NE.0) k(i+1,1)=
1860  & 1+(abs(kk1)+(2*iex-1)*kk1)/2/abs(kk1)
1861  p(i+1,1)=pjtx(jtp,i0)
1862  p(i+1,2)=pjty(jtp,i0)
1863  p(i+1,3)=pjtz(jtp,i0)
1864  p(i+1,4)=pjte(jtp,i0)
1865  p(i+1,5)=pjtm(jtp,i0)
1866 580 CONTINUE
1867  n=n+ntj(jtp)
1868  ENDIF
1869  IF(ihpr2(1).GT.0.AND.rlu(0).LE.hidat(3)) THEN
1870  hidat20=hidat(2)
1871  hipr150=hipr1(5)
1872  IF(ihpr2(8).EQ.0.AND.ihpr2(3).EQ.0.AND.ihpr2(9).EQ.0)
1873  & hidat(2)=2.0
1874  IF(hint1(1).GE.1000.0.AND.jetot.EQ.0)THEN
1875  hidat(2)=3.0
1876  hipr1(5)=5.0
1877  ENDIF
1878  CALL attrad(ierror)
1879  hidat(2)=hidat20
1880  hipr1(5)=hipr150
1881  ELSE IF(jetot.EQ.0.AND.ihpr2(1).GT.0.AND.
1882  & hint1(1).GE.1000.0.AND.
1883  & rlu(0).LE.0.8) THEN
1884  hidat20=hidat(2)
1885  hipr150=hipr1(5)
1886  hidat(2)=3.0
1887  hipr1(5)=5.0
1888  IF(ihpr2(8).EQ.0.AND.ihpr2(3).EQ.0.AND.ihpr2(9).EQ.0)
1889  & hidat(2)=2.0
1890  CALL attrad(ierror)
1891  hidat(2)=hidat20
1892  hipr1(5)=hipr150
1893  ENDIF
1894  IF(ierror.NE.0) RETURN
1895 C ******** conduct soft radiations
1896 C****************************
1897 C
1898 C
1899 C CALL LULIST(1)
1900  CALL luexec
1901  RETURN
1902 
1903 1000 n=1
1904  k(1,1)=1
1905  k(1,2)=nfp(jtp,3)
1906  DO 1100 jj=1,5
1907  p(1,jj)=pp(jtp,jj)
1908 1100 CONTINUE
1909 C ********proj remain as a nucleon or delta
1910  CALL luexec
1911 C call lulist(1)
1912  RETURN
1913 C
1914 1200 n=1
1915  k(1,1)=1
1916  k(1,2)=nft(jtp,3)
1917  DO 1300 jj=1,5
1918  p(1,jj)=pt(jtp,jj)
1919 1300 CONTINUE
1920 C ********targ remain as a nucleon or delta
1921  CALL luexec
1922 C call lulist(1)
1923  RETURN
1924  END
1925 C
1926 C
1927 C
1928 C********************************************************************
1929 C Sort the jets associated with a nucleon in order of their
1930 C rapdities
1931 C********************************************************************
1932  SUBROUTINE hijsrt(JPJT,NPT)
1933  dimension kf(100),px(100),py(100),pz(100),pe(100),pm(100)
1934  dimension y(100),ip(100,2)
1935  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
1936  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
1937  & pjpm(300,500),ntj(300),kftj(300,500),
1938  & pjtx(300,500),pjty(300,500),pjtz(300,500),
1939  & pjte(300,500),pjtm(300,500)
1940  SAVE /hijjet1/
1941  IF(npt.EQ.2) go to 500
1942  jp=jpjt
1943  iq=0
1944  i=1
1945 100 kf(i)=kfpj(jp,i)
1946  px(i)=pjpx(jp,i)
1947  py(i)=pjpy(jp,i)
1948  pz(i)=pjpz(jp,i)
1949  pe(i)=pjpe(jp,i)
1950  pm(i)=pjpm(jp,i)
1951  y(i-iq)=0.5*alog((abs(pe(i)+pz(i))+1.e-5)
1952  & /(abs(pe(i)-pz(i))+1.e-5))
1953  ip(i-iq,1)=i
1954  ip(i-iq,2)=0
1955  IF(kf(i).NE.21) THEN
1956  ip(i-iq,2)=1
1957  iq=iq+1
1958  i=i+1
1959  kf(i)=kfpj(jp,i)
1960  px(i)=pjpx(jp,i)
1961  py(i)=pjpy(jp,i)
1962  pz(i)=pjpz(jp,i)
1963  pe(i)=pjpe(jp,i)
1964  pm(i)=pjpm(jp,i)
1965  ENDIF
1966  i=i+1
1967  IF(i.LE.npj(jp)) go to 100
1968 
1969  DO 200 i=1,npj(jp)-iq
1970  DO 200 j=i+1,npj(jp)-iq
1971  IF(y(i).GT.y(j)) go to 200
1972  ip1=ip(i,1)
1973  ip2=ip(i,2)
1974  ip(i,1)=ip(j,1)
1975  ip(i,2)=ip(j,2)
1976  ip(j,1)=ip1
1977  ip(j,2)=ip2
1978 200 CONTINUE
1979 C ********sort in decending y
1980  iqq=0
1981  i=1
1982 300 kfpj(jp,i)=kf(ip(i-iqq,1))
1983  pjpx(jp,i)=px(ip(i-iqq,1))
1984  pjpy(jp,i)=py(ip(i-iqq,1))
1985  pjpz(jp,i)=pz(ip(i-iqq,1))
1986  pjpe(jp,i)=pe(ip(i-iqq,1))
1987  pjpm(jp,i)=pm(ip(i-iqq,1))
1988  IF(ip(i-iqq,2).EQ.1) THEN
1989  kfpj(jp,i+1)=kf(ip(i-iqq,1)+1)
1990  pjpx(jp,i+1)=px(ip(i-iqq,1)+1)
1991  pjpy(jp,i+1)=py(ip(i-iqq,1)+1)
1992  pjpz(jp,i+1)=pz(ip(i-iqq,1)+1)
1993  pjpe(jp,i+1)=pe(ip(i-iqq,1)+1)
1994  pjpm(jp,i+1)=pm(ip(i-iqq,1)+1)
1995  i=i+1
1996  iqq=iqq+1
1997  ENDIF
1998  i=i+1
1999  IF(i.LE.npj(jp)) go to 300
2000 
2001  RETURN
2002 
2003 500 jt=jpjt
2004  iq=0
2005  i=1
2006 600 kf(i)=kftj(jt,i)
2007  px(i)=pjtx(jt,i)
2008  py(i)=pjty(jt,i)
2009  pz(i)=pjtz(jt,i)
2010  pe(i)=pjte(jt,i)
2011  pm(i)=pjtm(jt,i)
2012  y(i-iq)=0.5*alog((abs(pe(i)+pz(i))+1.e-5)
2013  & /(abs(pe(i)-pz(i))+1.e-5))
2014  ip(i-iq,1)=i
2015  ip(i-iq,2)=0
2016  IF(kf(i).NE.21) THEN
2017  ip(i-iq,2)=1
2018  iq=iq+1
2019  i=i+1
2020  kf(i)=kftj(jt,i)
2021  px(i)=pjtx(jt,i)
2022  py(i)=pjty(jt,i)
2023  pz(i)=pjtz(jt,i)
2024  pe(i)=pjte(jt,i)
2025  pm(i)=pjtm(jt,i)
2026  ENDIF
2027  i=i+1
2028  IF(i.LE.ntj(jt)) go to 600
2029 
2030  DO 700 i=1,ntj(jt)-iq
2031  DO 700 j=i+1,ntj(jt)-iq
2032  IF(y(i).LT.y(j)) go to 700
2033  ip1=ip(i,1)
2034  ip2=ip(i,2)
2035  ip(i,1)=ip(j,1)
2036  ip(i,2)=ip(j,2)
2037  ip(j,1)=ip1
2038  ip(j,2)=ip2
2039 700 CONTINUE
2040 C ********sort in acending y
2041  iqq=0
2042  i=1
2043 800 kftj(jt,i)=kf(ip(i-iqq,1))
2044  pjtx(jt,i)=px(ip(i-iqq,1))
2045  pjty(jt,i)=py(ip(i-iqq,1))
2046  pjtz(jt,i)=pz(ip(i-iqq,1))
2047  pjte(jt,i)=pe(ip(i-iqq,1))
2048  pjtm(jt,i)=pm(ip(i-iqq,1))
2049  IF(ip(i-iqq,2).EQ.1) THEN
2050  kftj(jt,i+1)=kf(ip(i-iqq,1)+1)
2051  pjtx(jt,i+1)=px(ip(i-iqq,1)+1)
2052  pjty(jt,i+1)=py(ip(i-iqq,1)+1)
2053  pjtz(jt,i+1)=pz(ip(i-iqq,1)+1)
2054  pjte(jt,i+1)=pe(ip(i-iqq,1)+1)
2055  pjtm(jt,i+1)=pm(ip(i-iqq,1)+1)
2056  i=i+1
2057  iqq=iqq+1
2058  ENDIF
2059  i=i+1
2060  IF(i.LE.ntj(jt)) go to 800
2061  RETURN
2062  END
2063 
2064 C
2065 C
2066 C
2067 C****************************************************************
2068 C conduct soft radiation according to dipole approxiamtion
2069 C****************************************************************
2070  SUBROUTINE attrad(IERROR)
2071 C
2072  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2073  SAVE /hiparnt/
2074  common/hijdat/hidat0(10,10),hidat(10)
2075  SAVE /hijdat/
2076  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
2077  SAVE /lujets/
2078  ierror=0
2079 
2080 C.....S INVARIANT MASS-SQUARED BETWEEN PARTONS I AND I+1......
2081 C.....SM IS THE LARGEST MASS-SQUARED....
2082 
2083 40 sm=0.
2084  jl=1
2085  DO 30 i=1,n-1
2086  s=2.*(p(i,4)*p(i+1,4)-p(i,1)*p(i+1,1)-p(i,2)*p(i+1,2)
2087  & -p(i,3)*p(i+1,3))+p(i,5)**2+p(i+1,5)**2
2088  IF(s.LT.0.) s=0.
2089  wp=sqrt(s)-1.5*(p(i,5)+p(i+1,5))
2090  IF(wp.GT.sm) THEN
2091  pbt1=p(i,1)+p(i+1,1)
2092  pbt2=p(i,2)+p(i+1,2)
2093  pbt3=p(i,3)+p(i+1,3)
2094  pbt4=p(i,4)+p(i+1,4)
2095  btt=(pbt1**2+pbt2**2+pbt3**2)/pbt4**2
2096  IF(btt.GE.1.0-1.0e-10) go to 30
2097  IF((i.NE.1.OR.i.NE.n-1).AND.
2098  & (k(i,2).NE.21.AND.k(i+1,2).NE.21)) go to 30
2099  jl=i
2100  sm=wp
2101  ENDIF
2102 30 CONTINUE
2103  s=(sm+1.5*(p(jl,5)+p(jl+1,5)))**2
2104  IF(sm.LT.hipr1(5)) goto 2
2105 
2106 C.....MAKE PLACE FOR ONE GLUON.....
2107  IF(jl+1.EQ.n) goto 190
2108  DO 160 j=n,jl+2,-1
2109  k(j+1,1)=k(j,1)
2110  k(j+1,2)=k(j,2)
2111  DO 150 m=1,5
2112 150 p(j+1,m)=p(j,m)
2113 160 CONTINUE
2114 190 n=n+1
2115 
2116 C.....BOOST TO REST SYSTEM FOR PARTICLES JL AND JL+1.....
2117  p1=p(jl,1)+p(jl+1,1)
2118  p2=p(jl,2)+p(jl+1,2)
2119  p3=p(jl,3)+p(jl+1,3)
2120  p4=p(jl,4)+p(jl+1,4)
2121  bex=-p1/p4
2122  bey=-p2/p4
2123  bez=-p3/p4
2124  imin=jl
2125  imax=jl+1
2126  CALL atrobo(0.,0.,bex,bey,bez,imin,imax,ierror)
2127  IF(ierror.NE.0) RETURN
2128 C.....ROTATE TO Z-AXIS....
2129  cth=p(jl,3)/sqrt(p(jl,4)**2-p(jl,5)**2)
2130  IF(abs(cth).GT.1.0) cth=max(-1.,min(1.,cth))
2131  theta=acos(cth)
2132  phi=ulangl(p(jl,1),p(jl,2))
2133  CALL atrobo(0.,-phi,0.,0.,0.,imin,imax,ierror)
2134  CALL atrobo(-theta,0.,0.,0.,0.,imin,imax,ierror)
2135 
2136 C.....CREATE ONE GLUON AND ORIENTATE.....
2137 
2138 1 CALL ar3jet(s,x1,x3,jl)
2139  CALL arorie(s,x1,x3,jl)
2140  IF(hidat(2).GT.0.0) THEN
2141  ptg1=sqrt(p(jl,1)**2+p(jl,2)**2)
2142  ptg2=sqrt(p(jl+1,1)**2+p(jl+1,2)**2)
2143  ptg3=sqrt(p(jl+2,1)**2+p(jl+2,2)**2)
2144  ptg=max(ptg1,ptg2,ptg3)
2145  IF(ptg.GT.hidat(2)) THEN
2146  fmfact=exp(-(ptg**2-hidat(2)**2)/hipr1(2)**2)
2147  IF(rlu(0).GT.fmfact) go to 1
2148  ENDIF
2149  ENDIF
2150 C.....ROTATE AND BOOST BACK.....
2151  imin=jl
2152  imax=jl+2
2153  CALL atrobo(theta,phi,-bex,-bey,-bez,imin,imax,ierror)
2154  IF(ierror.NE.0) RETURN
2155 C.....ENUMERATE THE GLUONS.....
2156  k(jl+2,1)=k(jl+1,1)
2157  k(jl+2,2)=k(jl+1,2)
2158  k(jl+2,3)=k(jl+1,3)
2159  k(jl+2,4)=k(jl+1,4)
2160  k(jl+2,5)=k(jl+1,5)
2161  p(jl+2,5)=p(jl+1,5)
2162  k(jl+1,1)=2
2163  k(jl+1,2)=21
2164  k(jl+1,3)=0
2165  k(jl+1,4)=0
2166  k(jl+1,5)=0
2167  p(jl+1,5)=0.
2168 C----THETA FUNCTION DAMPING OF THE EMITTED GLUONS. FOR HADRON-HADRON.
2169 C----R0=VFR(2)
2170 C IF(VFR(2).GT.0.) THEN
2171 C PTG=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
2172 C PTGMAX=WSTRI/2.
2173 C DOPT=SQRT((4.*PAR(71)*VFR(2))/WSTRI)
2174 C PTOPT=(DOPT*WSTRI)/(2.*VFR(2))
2175 C IF(PTG.GT.PTOPT) IORDER=IORDER-1
2176 C IF(PTG.GT.PTOPT) GOTO 1
2177 C ENDIF
2178 C-----
2179  IF(sm.GE.hipr1(5)) goto 40
2180 
2181 2 k(1,1)=2
2182  k(1,3)=0
2183  k(1,4)=0
2184  k(1,5)=0
2185  k(n,1)=1
2186  k(n,3)=0
2187  k(n,4)=0
2188  k(n,5)=0
2189 
2190  RETURN
2191  END
2192 
2193 
2194  SUBROUTINE ar3jet(S,X1,X3,JL)
2195 C
2196  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2197  SAVE /hiparnt/
2198  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
2199  SAVE /lujets/
2200  common/ranseed/nseed
2201  SAVE /ranseed/
2202 C
2203  c=1./3.
2204  IF(k(jl,2).NE.21 .AND. k(jl+1,2).NE.21) c=8./27.
2205  exp1=3
2206  exp3=3
2207  IF(k(jl,2).NE.21) exp1=2
2208  IF(k(jl+1,2).NE.21) exp3=2
2209  a=0.24**2/s
2210  yma=alog(.5/sqrt(a)+sqrt(.25/a-1))
2211  d=4.*c*yma
2212  sm1=p(jl,5)**2/s
2213  sm3=p(jl+1,5)**2/s
2214  xt2m=(1.-2.*sqrt(sm1)+sm1-sm3)*(1.-2.*sqrt(sm3)-sm1+sm3)
2215  xt2m=min(.25,xt2m)
2216  ntry=0
2217 1 IF(ntry.EQ.5000) THEN
2218  x1=.5*(2.*sqrt(sm1)+1.+sm1-sm3)
2219  x3=.5*(2.*sqrt(sm3)+1.-sm1+sm3)
2220  RETURN
2221  ENDIF
2222  ntry=ntry+1
2223 
2224  xt2=a*(xt2m/a)**(rlu(0)**(1./d))
2225 
2226  ymax=alog(.5/sqrt(xt2)+sqrt(.25/xt2-1.))
2227  y=(2.*rlu(0)-1.)*ymax
2228  x1=1.-sqrt(xt2)*exp(y)
2229  x3=1.-sqrt(xt2)*exp(-y)
2230  x2=2.-x1-x3
2231  neg=0
2232  IF(k(jl,2).NE.21 .OR. k(jl+1,2).NE.21) THEN
2233  IF((1.-x1)*(1.-x2)*(1.-x3)-x2*sm1*(1.-x1)-x2*sm3*(1.-x3).
2234  & le.0..OR.x1.LE.2.*sqrt(sm1)-sm1+sm3.OR.x3.LE.2.*sqrt(sm3)
2235  & -sm3+sm1) neg=1
2236  x1=x1+sm1-sm3
2237  x3=x3-sm1+sm3
2238  ENDIF
2239  IF(neg.EQ.1) goto 1
2240 
2241  fg=2.*ymax*c*(x1**exp1+x3**exp3)/d
2242  xt2m=xt2
2243  IF(fg.LT.rlu(0)) goto 1
2244 
2245  RETURN
2246  END
2247 C*************************************************************
2248 
2249 
2250  SUBROUTINE arorie(S,X1,X3,JL)
2251 C
2252  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2253  SAVE /hiparnt/
2254  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
2255  SAVE /lujets/
2256  common/ranseed/nseed
2257  SAVE /ranseed/
2258 C
2259  w=sqrt(s)
2260  x2=2.-x1-x3
2261  e1=.5*x1*w
2262  e3=.5*x3*w
2263  p1=sqrt(e1**2-p(jl,5)**2)
2264  p3=sqrt(e3**2-p(jl+1,5)**2)
2265  cbet=1.
2266  IF(p1.GT.0..AND.p3.GT.0.) cbet=(p(jl,5)**2
2267  & +p(jl+1,5)**2+2.*e1*e3-s*(1.-x2))/(2.*p1*p3)
2268  IF(abs(cbet).GT.1.0) cbet=max(-1.,min(1.,cbet))
2269  bet=acos(cbet)
2270 
2271 C.....MINIMIZE PT1-SQUARED PLUS PT3-SQUARED.....
2272  IF(p1.GE.p3) THEN
2273  psi=.5*ulangl(p1**2+p3**2*cos(2.*bet),-p3**2*sin(2.*bet))
2274  pt1=p1*sin(psi)
2275  pz1=p1*cos(psi)
2276  pt3=p3*sin(psi+bet)
2277  pz3=p3*cos(psi+bet)
2278  ELSE IF(p3.GT.p1) THEN
2279  psi=.5*ulangl(p3**2+p1**2*cos(2.*bet),-p1**2*sin(2.*bet))
2280  pt1=p1*sin(bet+psi)
2281  pz1=-p1*cos(bet+psi)
2282  pt3=p3*sin(psi)
2283  pz3=-p3*cos(psi)
2284  ENDIF
2285 
2286  del=2.0*hipr1(40)*rlu(0)
2287  p(jl,4)=e1
2288  p(jl,1)=pt1*sin(del)
2289  p(jl,2)=-pt1*cos(del)
2290  p(jl,3)=pz1
2291  p(jl+2,4)=e3
2292  p(jl+2,1)=pt3*sin(del)
2293  p(jl+2,2)=-pt3*cos(del)
2294  p(jl+2,3)=pz3
2295  p(jl+1,4)=w-e1-e3
2296  p(jl+1,1)=-p(jl,1)-p(jl+2,1)
2297  p(jl+1,2)=-p(jl,2)-p(jl+2,2)
2298  p(jl+1,3)=-p(jl,3)-p(jl+2,3)
2299  RETURN
2300  END
2301 
2302 
2303 C
2304 C*******************************************************************
2305 C make boost and rotation to entries from IMIN to IMAX
2306 C*******************************************************************
2307  SUBROUTINE atrobo(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
2308  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
2309  SAVE /lujets/
2310  dimension rot(3,3),pv(3)
2311  DOUBLE PRECISION dp(4),dbex,dbey,dbez,dga,dga2,dbep,dgabep
2312  ierror=0
2313 
2314  IF(imin.LE.0 .OR. imax.GT.n .OR. imin.GT.imax) RETURN
2315 
2316  IF(the**2+phi**2.GT.1e-20) THEN
2317 C...ROTATE (TYPICALLY FROM Z AXIS TO DIRECTION THETA,PHI)
2318  rot(1,1)=cos(the)*cos(phi)
2319  rot(1,2)=-sin(phi)
2320  rot(1,3)=sin(the)*cos(phi)
2321  rot(2,1)=cos(the)*sin(phi)
2322  rot(2,2)=cos(phi)
2323  rot(2,3)=sin(the)*sin(phi)
2324  rot(3,1)=-sin(the)
2325  rot(3,2)=0.
2326  rot(3,3)=cos(the)
2327  DO 120 i=imin,imax
2328 C************** IF(MOD(K(I,1)/10000,10).GE.6) GOTO 120
2329  DO 100 j=1,3
2330  100 pv(j)=p(i,j)
2331  DO 110 j=1,3
2332  110 p(i,j)=rot(j,1)*pv(1)+rot(j,2)*pv(2)
2333  & +rot(j,3)*pv(3)
2334  120 CONTINUE
2335  ENDIF
2336 
2337  IF(bex**2+bey**2+bez**2.GT.1e-20) THEN
2338 C...LORENTZ BOOST (TYPICALLY FROM REST TO MOMENTUM/ENERGY=BETA)
2339  dbex=bex
2340  dbey=bey
2341  dbez=bez
2342  dga2=1d0-dbex**2-dbey**2-dbez**2
2343  IF(dga2.LE.0d0) THEN
2344  ierror=1
2345  RETURN
2346  ENDIF
2347  dga=1d0/dsqrt(dga2)
2348  DO 140 i=imin,imax
2349 C************* IF(MOD(K(I,1)/10000,10).GE.6) GOTO 140
2350  DO 130 j=1,4
2351 130 dp(j)=p(i,j)
2352  dbep=dbex*dp(1)+dbey*dp(2)+dbez*dp(3)
2353  dgabep=dga*(dga*dbep/(1d0+dga)+dp(4))
2354  p(i,1)=dp(1)+dgabep*dbex
2355  p(i,2)=dp(2)+dgabep*dbey
2356  p(i,3)=dp(3)+dgabep*dbez
2357  p(i,4)=dga*(dp(4)+dbep)
2358 140 CONTINUE
2359  ENDIF
2360 
2361  RETURN
2362  END
2363 C
2364 C
2365 C
2366  SUBROUTINE hijhrd(JP,JT,JOUT,JFLG,IOPJET0)
2367 C
2368 C IOPTJET=1, ALL JET WILL FORM SINGLE STRING SYSTEM
2369 C 0, ONLY Q-QBAR JET FORM SINGLE STRING SYSTEM
2370 C*******Perform jets production and fragmentation when JP JT *******
2371 C scatter. JOUT-> number of hard scatterings precede this one *
2372 C for the the same pair(JP,JT). JFLG->a flag to show whether *
2373 C jets can be produced (with valence quark=1,gluon=2, q-qbar=3)*
2374 C or not(0). Information of jets are in COMMON/ATTJET and *
2375 C /MINJET. ABS(NFP(JP,6)) is the total number jets produced by *
2376 C JP. If NFP(JP,6)<0 JP can not produce jet anymore. *
2377 C*******************************************************************
2378  dimension ip(100,2),ipq(50),ipb(50),it(100,2),itq(50),itb(50)
2379  common/hijcrdn/yp(3,300),yt(3,300)
2380  SAVE /hijcrdn/
2381  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
2382  SAVE /hiparnt/
2383  common/hijdat/hidat0(10,10),hidat(10)
2384  SAVE /hijdat/
2385  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
2386  SAVE /histrng/
2387  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
2388  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
2389  & pjpm(300,500),ntj(300),kftj(300,500),
2390  & pjtx(300,500),pjty(300,500),pjtz(300,500),
2391  & pjte(300,500),pjtm(300,500)
2392  SAVE /hijjet1/
2393  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
2394  & k2sg(900,100),pxsg(900,100),pysg(900,100),
2395  & pzsg(900,100),pesg(900,100),pmsg(900,100)
2396  SAVE /hijjet2/
2397  common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5)
2398  SAVE /hijjet4/
2399  common/ranseed/nseed
2400  SAVE /ranseed/
2401 C************************************ HIJING common block
2402  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
2403  SAVE /lujets/
2404  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
2405  SAVE /ludat1/
2406  common/pysubs/msel,msub(200),kfin(2,-40:40),ckin(200)
2407  SAVE /pysubs/
2408  common/pypars/mstp(200),parp(200),msti(200),pari(200)
2409  SAVE /pypars/
2410  common/pyint1/mint(400),vint(400)
2411  SAVE /pyint1/
2412  common/pyint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
2413  SAVE /pyint2/
2414  common/pyint5/ngen(0:200,3),xsec(0:200,3)
2415  SAVE /pyint5/
2416  common/hipyint/mint4,mint5,atco(200,20),atxs(0:200)
2417  SAVE /hipyint/
2418 C*********************************** LU common block
2419  mxjt=500
2420 C SIZE OF COMMON BLOCK FOR # OF PARTON PER STRING
2421  mxsg=900
2422 C SIZE OF COMMON BLOCK FOR # OF SINGLE STRINGS
2423  mxsj=100
2424 C SIZE OF COMMON BLOCK FOR # OF PARTON PER SINGLE
2425 C STRING
2426  jflg=0
2427  ihnt2(11)=jp
2428  ihnt2(12)=jt
2429 C
2430  iopjet=iopjet0
2431  IF(iopjet.EQ.1.AND.(nfp(jp,6).NE.0.OR.nft(jt,6).NE.0))
2432  & iopjet=0
2433  IF(jp.GT.ihnt2(1) .OR. jt.GT.ihnt2(3)) RETURN
2434  IF(nfp(jp,6).LT.0 .OR. nft(jt,6).LT.0) RETURN
2435 C ******** JP or JT can not produce jet anymore
2436 C
2437  IF(jout.EQ.0) THEN
2438  epp=pp(jp,4)+pp(jp,3)
2439  epm=pp(jp,4)-pp(jp,3)
2440  etp=pt(jt,4)+pt(jt,3)
2441  etm=pt(jt,4)-pt(jt,3)
2442  IF(epp.LT.0.0) go to 1000
2443  IF(epm.LT.0.0) go to 1000
2444  IF(etp.LT.0.0) go to 1000
2445  IF(etm.LT.0.0) go to 1000
2446  IF(epp/(epm+0.01).LE.etp/(etm+0.01)) RETURN
2447  ENDIF
2448 C ********for the first hard scattering of (JP,JT)
2449 C have collision only when Ycm(JP)>Ycm(JT)
2450 
2451  ecut1=hipr1(1)+hipr1(8)+pp(jp,14)+pp(jp,15)
2452  ecut2=hipr1(1)+hipr1(8)+pt(jt,14)+pt(jt,15)
2453  IF(pp(jp,4).LE.ecut1) THEN
2454  nfp(jp,6)=-abs(nfp(jp,6))
2455  RETURN
2456  ENDIF
2457  IF(pt(jt,4).LE.ecut2) THEN
2458  nft(jt,6)=-abs(nft(jt,6))
2459  RETURN
2460  ENDIF
2461 C *********must have enough energy to produce jets
2462 
2463  miss=0
2464  misp=0
2465  mist=0
2466 C
2467  IF(nfp(jp,10).EQ.0 .AND. nft(jt,10).EQ.0) THEN
2468  mint(44)=mint4
2469  mint(45)=mint5
2470  xsec(0,1)=atxs(0)
2471  xsec(11,1)=atxs(11)
2472  xsec(12,1)=atxs(12)
2473  xsec(28,1)=atxs(28)
2474  DO 120 i=1,20
2475  coef(11,i)=atco(11,i)
2476  coef(12,i)=atco(12,i)
2477  coef(28,i)=atco(28,i)
2478 120 CONTINUE
2479  ELSE
2480  isub11=0
2481  isub12=0
2482  isub28=0
2483  IF(xsec(11,1).NE.0) isub11=1
2484  IF(xsec(12,1).NE.0) isub12=1
2485  IF(xsec(28,1).NE.0) isub28=1
2486  mint(44)=mint4-isub11-isub12-isub28
2487  mint(45)=mint5-isub11-isub12-isub28
2488  xsec(0,1)=atxs(0)-atxs(11)-atxs(12)-atxs(28)
2489  xsec(11,1)=0.0
2490  xsec(12,1)=0.0
2491  xsec(28,1)=0.0
2492  DO 110 i=1,20
2493  coef(11,i)=0.0
2494  coef(12,i)=0.0
2495  coef(28,i)=0.0
2496 110 CONTINUE
2497  ENDIF
2498 C ********Scatter the valence quarks only once per NN
2499 C collision,
2500 C afterwards only gluon can have hard scattering.
2501  155 CALL pythia
2502  jj=mint(31)
2503  IF(jj.NE.1) go to 155
2504 C *********one hard collision at a time
2505  IF(k(7,2).EQ.-k(8,2)) THEN
2506  qmass2=(p(7,4)+p(8,4))**2-(p(7,1)+p(8,1))**2
2507  & -(p(7,2)+p(8,2))**2-(p(7,3)+p(8,3))**2
2508  qm=ulmass(k(7,2))
2509  IF(qmass2.LT.(2.0*qm+hipr1(1))**2) go to 155
2510  ENDIF
2511 C ********q-qbar jets must has minimum mass HIPR1(1)
2512  pxp=pp(jp,1)-p(3,1)
2513  pyp=pp(jp,2)-p(3,2)
2514  pzp=pp(jp,3)-p(3,3)
2515  pep=pp(jp,4)-p(3,4)
2516  pxt=pt(jt,1)-p(4,1)
2517  pyt=pt(jt,2)-p(4,2)
2518  pzt=pt(jt,3)-p(4,3)
2519  pet=pt(jt,4)-p(4,4)
2520 
2521  IF(pep.LE.ecut1) THEN
2522  misp=misp+1
2523  IF(misp.LT.50) go to 155
2524  nfp(jp,6)=-abs(nfp(jp,6))
2525  RETURN
2526  ENDIF
2527  IF(pet.LE.ecut2) THEN
2528  mist=mist+1
2529  IF(mist.LT.50) go to 155
2530  nft(jt,6)=-abs(nft(jt,6))
2531  RETURN
2532  ENDIF
2533 C ******** if the remain energy<ECUT the proj or targ
2534 C can not produce jet anymore
2535 
2536  wp=pep+pzp+pet+pzt
2537  wm=pep-pzp+pet-pzt
2538  IF(wp.LT.0.0 .OR. wm.LT.0.0) THEN
2539  miss=miss+1
2540  IF(miss.LT.50) go to 155
2541  RETURN
2542  ENDIF
2543 C ********the total W+, W- must be positive
2544  sw=wp*wm
2545  ampx=sqrt((ecut1-hipr1(8))**2+pxp**2+pyp**2+0.01)
2546  amtx=sqrt((ecut2-hipr1(8))**2+pxt**2+pyt**2+0.01)
2547  sxx=(ampx+amtx)**2
2548  IF(sw.LT.sxx.OR.vint(43).LT.hipr1(1)) THEN
2549  miss=miss+1
2550  IF(miss.LT.50) go to 155
2551  RETURN
2552  ENDIF
2553 C ********the proj and targ remnants must have at least
2554 C a CM energy that can produce two strings
2555 C with minimum mass HIPR1(1)(see HIJSFT HIJFRG)
2556 C
2557  hint1(41)=p(7,1)
2558  hint1(42)=p(7,2)
2559  hint1(43)=p(7,3)
2560  hint1(44)=p(7,4)
2561  hint1(45)=p(7,5)
2562  hint1(46)=sqrt(p(7,1)**2+p(7,2)**2)
2563  hint1(51)=p(8,1)
2564  hint1(52)=p(8,2)
2565  hint1(53)=p(8,3)
2566  hint1(54)=p(8,4)
2567  hint1(55)=p(8,5)
2568  hint1(56)=sqrt(p(8,1)**2+p(8,2)**2)
2569  ihnt2(14)=k(7,2)
2570  ihnt2(15)=k(8,2)
2571 C
2572  pinirad=(1.0-exp(-2.0*(vint(47)-hidat(1))))
2573  & /(1.0+exp(-2.0*(vint(47)-hidat(1))))
2574  i_inirad=0
2575  IF(rlu(0).LE.pinirad) i_inirad=1
2576  IF(k(7,2).EQ.-k(8,2)) go to 190
2577  IF(k(7,2).EQ.21.AND.k(8,2).EQ.21.AND.iopjet.EQ.1) go to 190
2578 C*******************************************************************
2579 C gluon jets are going to be connectd with
2580 C the final leading string of quark-aintquark
2581 C*******************************************************************
2582  jflg=2
2583  jpp=0
2584  lpq=0
2585  lpb=0
2586  jtt=0
2587  ltq=0
2588  ltb=0
2589  is7=0
2590  is8=0
2591  hint1(47)=0.0
2592  hint1(48)=0.0
2593  hint1(49)=0.0
2594  hint1(50)=0.0
2595  hint1(67)=0.0
2596  hint1(68)=0.0
2597  hint1(69)=0.0
2598  hint1(70)=0.0
2599  DO 180 i=9,n
2600  IF(k(i,3).EQ.1 .OR. k(i,3).EQ.2.OR.
2601  & abs(k(i,2)).GT.30) go to 180
2602 C************************************************************
2603  IF(k(i,3).EQ.7) THEN
2604  hint1(47)=hint1(47)+p(i,1)
2605  hint1(48)=hint1(48)+p(i,2)
2606  hint1(49)=hint1(49)+p(i,3)
2607  hint1(50)=hint1(50)+p(i,4)
2608  ENDIF
2609  IF(k(i,3).EQ.8) THEN
2610  hint1(67)=hint1(67)+p(i,1)
2611  hint1(68)=hint1(68)+p(i,2)
2612  hint1(69)=hint1(69)+p(i,3)
2613  hint1(70)=hint1(70)+p(i,4)
2614  ENDIF
2615 C************************modifcation made on Apr 10. 1996*****
2616  IF(k(i,2).GT.21.AND.k(i,2).LE.30) THEN
2617  ndr=ndr+1
2618  iadr(ndr,1)=jp
2619  iadr(ndr,2)=jt
2620  kfdr(ndr)=k(i,2)
2621  pdr(ndr,1)=p(i,1)
2622  pdr(ndr,2)=p(i,2)
2623  pdr(ndr,3)=p(i,3)
2624  pdr(ndr,4)=p(i,4)
2625  pdr(ndr,5)=p(i,5)
2626 C************************************************************
2627  go to 180
2628 C************************correction made on Oct. 14,1994*****
2629  ENDIF
2630  IF(k(i,3).EQ.7.OR.k(i,3).EQ.3) THEN
2631  IF(k(i,3).EQ.7.AND.k(i,2).NE.21.AND.k(i,2).EQ.k(7,2)
2632  & .AND.is7.EQ.0) THEN
2633  pp(jp,10)=p(i,1)
2634  pp(jp,11)=p(i,2)
2635  pp(jp,12)=p(i,3)
2636  pzp=pzp+p(i,3)
2637  pep=pep+p(i,4)
2638  nfp(jp,10)=1
2639  is7=1
2640  go to 180
2641  ENDIF
2642  IF(k(i,3).EQ.3.AND.(k(i,2).NE.21.OR.
2643  & i_inirad.EQ.0)) THEN
2644  pxp=pxp+p(i,1)
2645  pyp=pyp+p(i,2)
2646  pzp=pzp+p(i,3)
2647  pep=pep+p(i,4)
2648  go to 180
2649  ENDIF
2650  jpp=jpp+1
2651  ip(jpp,1)=i
2652  ip(jpp,2)=0
2653  IF(k(i,2).NE.21) THEN
2654  IF(k(i,2).GT.0) THEN
2655  lpq=lpq+1
2656  ipq(lpq)=jpp
2657  ip(jpp,2)=lpq
2658  ELSE IF(k(i,2).LT.0) THEN
2659  lpb=lpb+1
2660  ipb(lpb)=jpp
2661  ip(jpp,2)=-lpb
2662  ENDIF
2663  ENDIF
2664  ELSE IF(k(i,3).EQ.8.OR.k(i,3).EQ.4) THEN
2665  IF(k(i,3).EQ.8.AND.k(i,2).NE.21.AND.k(i,2).EQ.k(8,2)
2666  & .AND.is8.EQ.0) THEN
2667  pt(jt,10)=p(i,1)
2668  pt(jt,11)=p(i,2)
2669  pt(jt,12)=p(i,3)
2670  pzt=pzt+p(i,3)
2671  pet=pet+p(i,4)
2672  nft(jt,10)=1
2673  is8=1
2674  go to 180
2675  ENDIF
2676  IF(k(i,3).EQ.4.AND.(k(i,2).NE.21.OR.
2677  & i_inirad.EQ.0)) THEN
2678  pxt=pxt+p(i,1)
2679  pyt=pyt+p(i,2)
2680  pzt=pzt+p(i,3)
2681  pet=pet+p(i,4)
2682  go to 180
2683  ENDIF
2684  jtt=jtt+1
2685  it(jtt,1)=i
2686  it(jtt,2)=0
2687  IF(k(i,2).NE.21) THEN
2688  IF(k(i,2).GT.0) THEN
2689  ltq=ltq+1
2690  itq(ltq)=jtt
2691  it(jtt,2)=ltq
2692  ELSE IF(k(i,2).LT.0) THEN
2693  ltb=ltb+1
2694  itb(ltb)=jtt
2695  it(jtt,2)=-ltb
2696  ENDIF
2697  ENDIF
2698  ENDIF
2699  180 CONTINUE
2700 c
2701 c
2702  IF(lpq.NE.lpb .OR. ltq.NE.ltb) THEN
2703  miss=miss+1
2704  IF(miss.LE.50) go to 155
2705  WRITE(6,*) ' Q -QBAR NOT MATCHED IN HIJHRD'
2706  jflg=0
2707  RETURN
2708  ENDIF
2709 C****The following will rearrange the partons so that a quark is***
2710 C****allways followed by an anti-quark ****************************
2711 
2712  j=0
2713 181 j=j+1
2714  IF(j.GT.jpp) go to 182
2715  IF(ip(j,2).EQ.0) THEN
2716  go to 181
2717  ELSE IF(ip(j,2).NE.0) THEN
2718  lp=abs(ip(j,2))
2719  ip1=ip(j,1)
2720  ip2=ip(j,2)
2721  ip(j,1)=ip(ipq(lp),1)
2722  ip(j,2)=ip(ipq(lp),2)
2723  ip(ipq(lp),1)=ip1
2724  ip(ipq(lp),2)=ip2
2725  IF(ip2.GT.0) THEN
2726  ipq(ip2)=ipq(lp)
2727  ELSE IF(ip2.LT.0) THEN
2728  ipb(-ip2)=ipq(lp)
2729  ENDIF
2730 C ********replace J with a quark
2731  ip1=ip(j+1,1)
2732  ip2=ip(j+1,2)
2733  ip(j+1,1)=ip(ipb(lp),1)
2734  ip(j+1,2)=ip(ipb(lp),2)
2735  ip(ipb(lp),1)=ip1
2736  ip(ipb(lp),2)=ip2
2737  IF(ip2.GT.0) THEN
2738  ipq(ip2)=ipb(lp)
2739  ELSE IF(ip2.LT.0) THEN
2740  ipb(-ip2)=ipb(lp)
2741  ENDIF
2742 C ******** replace J+1 with anti-quark
2743  j=j+1
2744  go to 181
2745  ENDIF
2746 
2747 182 j=0
2748 183 j=j+1
2749  IF(j.GT.jtt) go to 184
2750  IF(it(j,2).EQ.0) THEN
2751  go to 183
2752  ELSE IF(it(j,2).NE.0) THEN
2753  lt=abs(it(j,2))
2754  it1=it(j,1)
2755  it2=it(j,2)
2756  it(j,1)=it(itq(lt),1)
2757  it(j,2)=it(itq(lt),2)
2758  it(itq(lt),1)=it1
2759  it(itq(lt),2)=it2
2760  IF(it2.GT.0) THEN
2761  itq(it2)=itq(lt)
2762  ELSE IF(it2.LT.0) THEN
2763  itb(-it2)=itq(lt)
2764  ENDIF
2765 C ********replace J with a quark
2766  it1=it(j+1,1)
2767  it2=it(j+1,2)
2768  it(j+1,1)=it(itb(lt),1)
2769  it(j+1,2)=it(itb(lt),2)
2770  it(itb(lt),1)=it1
2771  it(itb(lt),2)=it2
2772  IF(it2.GT.0) THEN
2773  itq(it2)=itb(lt)
2774  ELSE IF(it2.LT.0) THEN
2775  itb(-it2)=itb(lt)
2776  ENDIF
2777 C ******** replace J+1 with anti-quark
2778  j=j+1
2779  go to 183
2780 
2781  ENDIF
2782 
2783 184 CONTINUE
2784  IF(npj(jp)+jpp.GT.mxjt.OR.ntj(jt)+jtt.GT.mxjt) THEN
2785  jflg=0
2786  WRITE(6,*) 'number of partons per string exceeds'
2787  WRITE(6,*) 'the common block size'
2788  RETURN
2789  ENDIF
2790 C ********check the bounds of common blocks
2791  DO 186 j=1,jpp
2792  kfpj(jp,npj(jp)+j)=k(ip(j,1),2)
2793  pjpx(jp,npj(jp)+j)=p(ip(j,1),1)
2794  pjpy(jp,npj(jp)+j)=p(ip(j,1),2)
2795  pjpz(jp,npj(jp)+j)=p(ip(j,1),3)
2796  pjpe(jp,npj(jp)+j)=p(ip(j,1),4)
2797  pjpm(jp,npj(jp)+j)=p(ip(j,1),5)
2798 186 CONTINUE
2799  npj(jp)=npj(jp)+jpp
2800  DO 188 j=1,jtt
2801  kftj(jt,ntj(jt)+j)=k(it(j,1),2)
2802  pjtx(jt,ntj(jt)+j)=p(it(j,1),1)
2803  pjty(jt,ntj(jt)+j)=p(it(j,1),2)
2804  pjtz(jt,ntj(jt)+j)=p(it(j,1),3)
2805  pjte(jt,ntj(jt)+j)=p(it(j,1),4)
2806  pjtm(jt,ntj(jt)+j)=p(it(j,1),5)
2807 188 CONTINUE
2808  ntj(jt)=ntj(jt)+jtt
2809  go to 900
2810 C*****************************************************************
2811 CThis is the case of a quark-antiquark jet it will fragment alone
2812 C****************************************************************
2813 190 jflg=3
2814  IF(k(7,2).NE.21.AND.k(8,2).NE.21.AND.
2815  & k(7,2)*k(8,2).GT.0) go to 155
2816  jpp=0
2817  lpq=0
2818  lpb=0
2819  DO 200 i=9,n
2820  IF(k(i,3).EQ.1.OR.k(i,3).EQ.2.OR.
2821  & abs(k(i,2)).GT.30) go to 200
2822  IF(k(i,2).GT.21.AND.k(i,2).LE.30) THEN
2823  ndr=ndr+1
2824  iadr(ndr,1)=jp
2825  iadr(ndr,2)=jt
2826  kfdr(ndr)=k(i,2)
2827  pdr(ndr,1)=p(i,1)
2828  pdr(ndr,2)=p(i,2)
2829  pdr(ndr,3)=p(i,3)
2830  pdr(ndr,4)=p(i,4)
2831  pdr(ndr,5)=p(i,5)
2832 C************************************************************
2833  go to 200
2834 C************************correction made on Oct. 14,1994*****
2835  ENDIF
2836  IF(k(i,3).EQ.3.AND.(k(i,2).NE.21.OR.
2837  & i_inirad.EQ.0)) THEN
2838  pxp=pxp+p(i,1)
2839  pyp=pyp+p(i,2)
2840  pzp=pzp+p(i,3)
2841  pep=pep+p(i,4)
2842  go to 200
2843  ENDIF
2844  IF(k(i,3).EQ.4.AND.(k(i,2).NE.21.OR.
2845  & i_inirad.EQ.0)) THEN
2846  pxt=pxt+p(i,1)
2847  pyt=pyt+p(i,2)
2848  pzt=pzt+p(i,3)
2849  pet=pet+p(i,4)
2850  go to 200
2851  ENDIF
2852  jpp=jpp+1
2853  ip(jpp,1)=i
2854  ip(jpp,2)=0
2855  IF(k(i,2).NE.21) THEN
2856  IF(k(i,2).GT.0) THEN
2857  lpq=lpq+1
2858  ipq(lpq)=jpp
2859  ip(jpp,2)=lpq
2860  ELSE IF(k(i,2).LT.0) THEN
2861  lpb=lpb+1
2862  ipb(lpb)=jpp
2863  ip(jpp,2)=-lpb
2864  ENDIF
2865  ENDIF
2866 200 CONTINUE
2867  IF(lpq.NE.lpb) THEN
2868  miss=miss+1
2869  IF(miss.LE.50) go to 155
2870  WRITE(6,*) lpq,lpb, 'Q-QBAR NOT CONSERVED OR NOT MATCHED'
2871  jflg=0
2872  RETURN
2873  ENDIF
2874 
2875 C**** The following will rearrange the partons so that a quark is***
2876 C**** allways followed by an anti-quark ****************************
2877  j=0
2878 220 j=j+1
2879  IF(j.GT.jpp) go to 222
2880  IF(ip(j,2).EQ.0) go to 220
2881  lp=abs(ip(j,2))
2882  ip1=ip(j,1)
2883  ip2=ip(j,2)
2884  ip(j,1)=ip(ipq(lp),1)
2885  ip(j,2)=ip(ipq(lp),2)
2886  ip(ipq(lp),1)=ip1
2887  ip(ipq(lp),2)=ip2
2888  IF(ip2.GT.0) THEN
2889  ipq(ip2)=ipq(lp)
2890  ELSE IF(ip2.LT.0) THEN
2891  ipb(-ip2)=ipq(lp)
2892  ENDIF
2893  ipq(lp)=j
2894 C ********replace J with a quark
2895  ip1=ip(j+1,1)
2896  ip2=ip(j+1,2)
2897  ip(j+1,1)=ip(ipb(lp),1)
2898  ip(j+1,2)=ip(ipb(lp),2)
2899  ip(ipb(lp),1)=ip1
2900  ip(ipb(lp),2)=ip2
2901  IF(ip2.GT.0) THEN
2902  ipq(ip2)=ipb(lp)
2903  ELSE IF(ip2.LT.0) THEN
2904  ipb(-ip2)=ipb(lp)
2905  ENDIF
2906 C ******** replace J+1 with an anti-quark
2907  ipb(lp)=j+1
2908  j=j+1
2909  go to 220
2910 
2911 222 CONTINUE
2912  IF(lpq.GE.1) THEN
2913  DO 240 l0=2,lpq
2914  ip1=ip(2*l0-3,1)
2915  ip2=ip(2*l0-3,2)
2916  ip(2*l0-3,1)=ip(ipq(l0),1)
2917  ip(2*l0-3,2)=ip(ipq(l0),2)
2918  ip(ipq(l0),1)=ip1
2919  ip(ipq(l0),2)=ip2
2920  IF(ip2.GT.0) THEN
2921  ipq(ip2)=ipq(l0)
2922  ELSE IF(ip2.LT.0) THEN
2923  ipb(-ip2)=ipq(l0)
2924  ENDIF
2925  ipq(l0)=2*l0-3
2926 C
2927  ip1=ip(2*l0-2,1)
2928  ip2=ip(2*l0-2,2)
2929  ip(2*l0-2,1)=ip(ipb(l0),1)
2930  ip(2*l0-2,2)=ip(ipb(l0),2)
2931  ip(ipb(l0),1)=ip1
2932  ip(ipb(l0),2)=ip2
2933  IF(ip2.GT.0) THEN
2934  ipq(ip2)=ipb(l0)
2935  ELSE IF(ip2.LT.0) THEN
2936  ipb(-ip2)=ipb(l0)
2937  ENDIF
2938  ipb(l0)=2*l0-2
2939 240 CONTINUE
2940 C ********move all the qqbar pair to the front of
2941 C the list, except the first pair
2942  ip1=ip(2*lpq-1,1)
2943  ip2=ip(2*lpq-1,2)
2944  ip(2*lpq-1,1)=ip(ipq(1),1)
2945  ip(2*lpq-1,2)=ip(ipq(1),2)
2946  ip(ipq(1),1)=ip1
2947  ip(ipq(1),2)=ip2
2948  IF(ip2.GT.0) THEN
2949  ipq(ip2)=ipq(1)
2950  ELSE IF(ip2.LT.0) THEN
2951  ipb(-ip2)=ipq(1)
2952  ENDIF
2953  ipq(1)=2*lpq-1
2954 C ********move the first quark to the beginning of
2955 C the last string system
2956  ip1=ip(jpp,1)
2957  ip2=ip(jpp,2)
2958  ip(jpp,1)=ip(ipb(1),1)
2959  ip(jpp,2)=ip(ipb(1),2)
2960  ip(ipb(1),1)=ip1
2961  ip(ipb(1),2)=ip2
2962  IF(ip2.GT.0) THEN
2963  ipq(ip2)=ipb(1)
2964  ELSE IF(ip2.LT.0) THEN
2965  ipb(-ip2)=ipb(1)
2966  ENDIF
2967  ipb(1)=jpp
2968 C ********move the first anti-quark to the end of the
2969 C last string system
2970  ENDIF
2971  IF(nsg.GE.mxsg) THEN
2972  jflg=0
2973  WRITE(6,*) 'number of jets forming single strings exceeds'
2974  WRITE(6,*) 'the common block size'
2975  RETURN
2976  ENDIF
2977  IF(jpp.GT.mxsj) THEN
2978  jflg=0
2979  WRITE(6,*) 'number of partons per single jet system'
2980  WRITE(6,*) 'exceeds the common block size'
2981  RETURN
2982  ENDIF
2983 C ********check the bounds of common block size
2984  nsg=nsg+1
2985  njsg(nsg)=jpp
2986  iasg(nsg,1)=jp
2987  iasg(nsg,2)=jt
2988  iasg(nsg,3)=0
2989  DO 300 i=1,jpp
2990  k1sg(nsg,i)=2
2991  k2sg(nsg,i)=k(ip(i,1),2)
2992  IF(k2sg(nsg,i).LT.0) k1sg(nsg,i)=1
2993  pxsg(nsg,i)=p(ip(i,1),1)
2994  pysg(nsg,i)=p(ip(i,1),2)
2995  pzsg(nsg,i)=p(ip(i,1),3)
2996  pesg(nsg,i)=p(ip(i,1),4)
2997  pmsg(nsg,i)=p(ip(i,1),5)
2998 300 CONTINUE
2999  k1sg(nsg,1)=2
3000  k1sg(nsg,jpp)=1
3001 C******* reset the energy-momentum of incoming particles ********
3002 900 pp(jp,1)=pxp
3003  pp(jp,2)=pyp
3004  pp(jp,3)=pzp
3005  pp(jp,4)=pep
3006  pp(jp,5)=0.0
3007  pt(jt,1)=pxt
3008  pt(jt,2)=pyt
3009  pt(jt,3)=pzt
3010  pt(jt,4)=pet
3011  pt(jt,5)=0.0
3012 
3013  nfp(jp,6)=nfp(jp,6)+1
3014  nft(jt,6)=nft(jt,6)+1
3015  RETURN
3016 C
3017 1000 jflg=-1
3018  IF(ihpr2(10).EQ.0) RETURN
3019  WRITE(6,*) 'Fatal HIJHRD error'
3020  WRITE(6,*) jp, ' proj E+,E-',epp,epm,' status',nfp(jp,5)
3021  WRITE(6,*) jt, ' targ E+,E_',etp,etm,' status',nft(jt,5)
3022  RETURN
3023  END
3024 C
3025 C
3026 C
3027 C
3028 C
3029  SUBROUTINE jetini(JP,JT,I_TRIG)
3030 C*******Initialize PYTHIA for jet production**********************
3031 C I_TRIG=0: for normal processes
3032 C I_TRIG=1: for triggered processes
3033 C JP: sequence number of the projectile
3034 C JT: sequence number of the target
3035 C For A+A collisions, one has to initilize pythia
3036 C separately for each type of collisions, pp, pn,np and nn,
3037 C or hp and hn for hA collisions. In this subroutine we use the following
3038 C catalogue for different type of collisions:
3039 C h+h: h+h (I_TYPE=1)
3040 C h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
3041 C A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
3042 C A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
3043 C*****************************************************************
3044  CHARACTER beam*16,targ*16
3045  dimension xsec0(8,0:200),coef0(8,200,20),ini(8),
3046  & mint44(8),mint45(8)
3047  common/hijcrdn/yp(3,300),yt(3,300)
3048  SAVE /hijcrdn/
3049  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3050  SAVE /hiparnt/
3051  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3052  SAVE /histrng/
3053  common/hipyint/mint4,mint5,atco(200,20),atxs(0:200)
3054  SAVE /hipyint/
3055 C
3056  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
3057  SAVE /ludat1/
3058  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
3059  SAVE /ludat3/
3060  common/pysubs/msel,msub(200),kfin(2,-40:40),ckin(200)
3061  SAVE /pysubs/
3062  common/pypars/mstp(200),parp(200),msti(200),pari(200)
3063  SAVE /pypars/
3064  common/pyint1/mint(400),vint(400)
3065  SAVE /pyint1/
3066  common/pyint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
3067  SAVE /pyint2/
3068  common/pyint5/ngen(0:200,3),xsec(0:200,3)
3069  SAVE /pyint5/
3070  DATA ini/8*0/i_last/-1/
3071 
3072  save ini, i_last !uzhi
3073 
3074  ihnt2(11)=jp
3075  ihnt2(12)=jt
3076  IF(ihnt2(5).NE.0 .AND. ihnt2(6).NE.0) THEN
3077  i_type=1
3078  ELSE IF(ihnt2(5).NE.0 .AND. ihnt2(6).EQ.0) THEN
3079  i_type=1
3080  IF(nft(jt,4).EQ.2112) i_type=2
3081  ELSE IF(ihnt2(5).EQ.0 .AND. ihnt2(6).NE.0) THEN
3082  i_type=1
3083  IF(nfp(jp,4).EQ.2112) i_type=2
3084  ELSE
3085  IF(nfp(jp,4).EQ.2212 .AND. nft(jt,4).EQ.2212) THEN
3086  i_type=1
3087  ELSE IF(nfp(jp,4).EQ.2212 .AND. nft(jt,4).EQ.2112) THEN
3088  i_type=2
3089  ELSE IF(nfp(jp,4).EQ.2112 .AND. nft(jt,4).EQ.2212) THEN
3090  i_type=3
3091  ELSE
3092  i_type=4
3093  ENDIF
3094  ENDIF
3095 c
3096  IF(i_trig.NE.0) go to 160
3097  IF(i_trig.EQ.i_last) go to 150
3098  mstp(2)=2
3099 c ********second order running alpha_strong
3100  mstp(33)=1
3101  parp(31)=hipr1(17)
3102 C ********inclusion of K factor
3103  mstp(51)=3
3104 C ********Duke-Owens set 1 structure functions
3105  mstp(61)=1
3106 C ********INITIAL STATE RADIATION
3107  mstp(71)=1
3108 C ********FINAL STATE RADIATION
3109  IF(ihpr2(2).EQ.0.OR.ihpr2(2).EQ.2) mstp(61)=0
3110  IF(ihpr2(2).EQ.0.OR.ihpr2(2).EQ.1) mstp(71)=0
3111 c
3112  mstp(81)=0
3113 C ******** NO MULTIPLE INTERACTION
3114  mstp(82)=1
3115 C *******STRUCTURE OF MUTLIPLE INTERACTION
3116  mstp(111)=0
3117 C ********frag off(have to be done by local call)
3118  IF(ihpr2(10).EQ.0) mstp(122)=0
3119 C ********No printout of initialization information
3120  parp(81)=hipr1(8)
3121  ckin(5)=hipr1(8)
3122  ckin(3)=hipr1(8)
3123  ckin(4)=hipr1(9)
3124  IF(hipr1(9).LE.hipr1(8)) ckin(4)=-1.0
3125  ckin(9)=-10.0
3126  ckin(10)=10.0
3127  msel=0
3128  DO 100 isub=1,200
3129  msub(isub)=0
3130  100 CONTINUE
3131  msub(11)=1
3132  msub(12)=1
3133  msub(13)=1
3134  msub(28)=1
3135  msub(53)=1
3136  msub(68)=1
3137  msub(81)=1
3138  msub(82)=1
3139  DO 110 j=1,min(8,mdcy(21,3))
3140  110 mdme(mdcy(21,2)+j-1,1)=0
3141  isel=4
3142  IF(hint1(1).GE.20.0 .and. ihpr2(18).EQ.1) isel=5
3143  mdme(mdcy(21,2)+isel-1,1)=1
3144 C ********QCD subprocesses
3145  msub(14)=1
3146  msub(18)=1
3147  msub(29)=1
3148 C ******* direct photon production
3149  150 IF(ini(i_type).NE.0) go to 800
3150  go to 400
3151 C
3152 C *****triggered subprocesses, jet, photon, heavy quark and DY
3153 C
3154  160 i_type=4+i_type
3155  IF(i_trig.EQ.i_last) go to 260
3156  parp(81)=abs(hipr1(10))-0.25
3157  ckin(5)=abs(hipr1(10))-0.25
3158  ckin(3)=abs(hipr1(10))-0.25
3159  ckin(4)=abs(hipr1(10))+0.25
3160  IF(hipr1(10).LT.hipr1(8)) ckin(4)=-1.0
3161 c
3162  msel=0
3163  DO 101 isub=1,200
3164  msub(isub)=0
3165  101 CONTINUE
3166  IF(ihpr2(3).EQ.1) THEN
3167  msub(11)=1
3168  msub(12)=1
3169  msub(13)=1
3170  msub(28)=1
3171  msub(53)=1
3172  msub(68)=1
3173  msub(81)=1
3174  msub(82)=1
3175  msub(14)=1
3176  msub(18)=1
3177  msub(29)=1
3178  DO 102 j=1,min(8,mdcy(21,3))
3179  102 mdme(mdcy(21,2)+j-1,1)=0
3180  isel=4
3181  IF(hint1(1).GE.20.0 .and. ihpr2(18).EQ.1) isel=5
3182  mdme(mdcy(21,2)+isel-1,1)=1
3183 C ********QCD subprocesses
3184  ELSE IF(ihpr2(3).EQ.2) THEN
3185  msub(14)=1
3186  msub(18)=1
3187  msub(29)=1
3188 C ********Direct photon production
3189 c q+qbar->g+gamma,q+qbar->gamma+gamma, q+g->q+gamma
3190  ELSE IF(ihpr2(3).EQ.3) THEN
3191  ckin(3)=max(0.0,hipr1(10))
3192  ckin(5)=hipr1(8)
3193  parp(81)=hipr1(8)
3194  msub(81)=1
3195  msub(82)=1
3196  DO 105 j=1,min(8,mdcy(21,3))
3197  105 mdme(mdcy(21,2)+j-1,1)=0
3198  isel=4
3199  IF(hint1(1).GE.20.0 .and. ihpr2(18).EQ.1) isel=5
3200  mdme(mdcy(21,2)+isel-1,1)=1
3201 C **********Heavy quark production
3202  ENDIF
3203 260 IF(ini(i_type).NE.0) go to 800
3204 C
3205 C
3206 400 ini(i_type)=1
3207  IF(ihpr2(10).EQ.0) mstp(122)=0
3208  IF(nfp(jp,4).EQ.2212) THEN
3209  beam='P'
3210  ELSE IF(nfp(jp,4).EQ.-2212) THEN
3211  beam='P~'
3212  ELSE IF(nfp(jp,4).EQ.2112) THEN
3213  beam='N'
3214  ELSE IF(nfp(jp,4).EQ.-2112) THEN
3215  beam='N~'
3216  ELSE IF(nfp(jp,4).EQ.211) THEN
3217  beam='PI+'
3218  ELSE IF(nfp(jp,4).EQ.-211) THEN
3219  beam='PI-'
3220  ELSE IF(nfp(jp,4).EQ.321) THEN
3221  beam='PI+'
3222  ELSE IF(nfp(jp,4).EQ.-321) THEN
3223  beam='PI-'
3224  ELSE
3225  WRITE(6,*) 'unavailable beam type', nfp(jp,4)
3226  ENDIF
3227  IF(nft(jt,4).EQ.2212) THEN
3228  targ='P'
3229  ELSE IF(nft(jt,4).EQ.-2212) THEN
3230  targ='P~'
3231  ELSE IF(nft(jt,4).EQ.2112) THEN
3232  targ='N'
3233  ELSE IF(nft(jt,4).EQ.-2112) THEN
3234  targ='N~'
3235  ELSE IF(nft(jt,4).EQ.211) THEN
3236  targ='PI+'
3237  ELSE IF(nft(jt,4).EQ.-211) THEN
3238  targ='PI-'
3239  ELSE IF(nft(jt,4).EQ.321) THEN
3240  targ='PI+'
3241  ELSE IF(nft(jt,4).EQ.-321) THEN
3242  targ='PI-'
3243  ELSE
3244  WRITE(6,*) 'unavailable target type', nft(jt,4)
3245  ENDIF
3246 C
3247  ihnt2(16)=1
3248 C ******************indicate for initialization use when
3249 C structure functions are called in PYTHIA
3250 C
3251  CALL pyinit('CMS',beam,targ,hint1(1))
3252  mint4=mint(44)
3253  mint5=mint(45)
3254  mint44(i_type)=mint(44)
3255  mint45(i_type)=mint(45)
3256  atxs(0)=xsec(0,1)
3257  xsec0(i_type,0)=xsec(0,1)
3258  DO 500 i=1,200
3259  atxs(i)=xsec(i,1)
3260  xsec0(i_type,i)=xsec(i,1)
3261  DO 500 j=1,20
3262  atco(i,j)=coef(i,j)
3263  coef0(i_type,i,j)=coef(i,j)
3264 500 CONTINUE
3265 C
3266  ihnt2(16)=0
3267 C
3268  RETURN
3269 C ********Store the initialization information for
3270 C late use
3271 C
3272 C
3273 800 mint(44)=mint44(i_type)
3274  mint(45)=mint45(i_type)
3275  mint4=mint(44)
3276  mint5=mint(45)
3277  xsec(0,1)=xsec0(i_type,0)
3278  atxs(0)=xsec(0,1)
3279  DO 900 i=1,200
3280  xsec(i,1)=xsec0(i_type,i)
3281  atxs(i)=xsec(i,1)
3282  DO 900 j=1,20
3283  coef(i,j)=coef0(i_type,i,j)
3284  atco(i,j)=coef(i,j)
3285 900 CONTINUE
3286  i_last=i_trig
3287  mint(11)=nfp(jp,4)
3288  mint(12)=nft(jt,4)
3289  RETURN
3290  END
3291 C
3292 C
3293 C
3294  SUBROUTINE hijini
3295  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3296  SAVE /hiparnt/
3297  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3298  SAVE /histrng/
3299  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
3300  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
3301  & pjpm(300,500),ntj(300),kftj(300,500),
3302  & pjtx(300,500),pjty(300,500),pjtz(300,500),
3303  & pjte(300,500),pjtm(300,500)
3304  SAVE /hijjet1/
3305  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
3306  & k2sg(900,100),pxsg(900,100),pysg(900,100),
3307  & pzsg(900,100),pesg(900,100),pmsg(900,100)
3308  SAVE /hijjet2/
3309  common/hijjet4/ndr,iadr(900,2),kfdr(900),pdr(900,5)
3310  SAVE /hijjet4/
3311  common/ranseed/nseed
3312  SAVE /ranseed/
3313 C****************Reset the momentum of initial particles************
3314 C and assign flavors to the proj and targ string *
3315 C*******************************************************************
3316  nsg=0
3317  ndr=0
3318  ipp=2212
3319  ipt=2212
3320  IF(ihnt2(5).NE.0) ipp=ihnt2(5)
3321  IF(ihnt2(6).NE.0) ipt=ihnt2(6)
3322 C ********in case the proj or targ is a hadron.
3323 C
3324  DO 100 i=1,ihnt2(1)
3325  pp(i,1)=0.0
3326  pp(i,2)=0.0
3327  pp(i,3)=sqrt(hint1(1)**2/4.0-hint1(8)**2)
3328  pp(i,4)=hint1(1)/2
3329  pp(i,5)=hint1(8)
3330  pp(i,6)=0.0
3331  pp(i,7)=0.0
3332  pp(i,8)=0.0
3333  pp(i,9)=0.0
3334  pp(i,10)=0.0
3335  nfp(i,3)=ipp
3336  nfp(i,4)=ipp
3337  nfp(i,5)=0
3338  nfp(i,6)=0
3339  nfp(i,7)=0
3340  nfp(i,8)=0
3341  nfp(i,9)=0
3342  nfp(i,10)=0
3343  nfp(i,11)=0
3344  npj(i)=0
3345  IF(i.GT.abs(ihnt2(2))) nfp(i,3)=2112
3346  CALL attflv(nfp(i,3),idq,idqq)
3347  nfp(i,1)=idq
3348  nfp(i,2)=idqq
3349  nfp(i,15)=-1
3350  IF(abs(idq).GT.1000.OR.(abs(idq*idqq).LT.100.AND.
3351  & rlu(0).LT.0.5)) nfp(i,15)=1
3352  pp(i,14)=ulmass(idq)
3353  pp(i,15)=ulmass(idqq)
3354 100 CONTINUE
3355 C
3356  DO 200 i=1,ihnt2(3)
3357  pt(i,1)=0.0
3358  pt(i,2)=0.0
3359  pt(i,3)=-sqrt(hint1(1)**2/4.0-hint1(9)**2)
3360  pt(i,4)=hint1(1)/2.0
3361  pt(i,5)=hint1(9)
3362  pt(i,6)=0.0
3363  pt(i,7)=0.0
3364  pt(i,8)=0.0
3365  pt(i,9)=0.0
3366  pt(i,10)=0.0
3367  nft(i,3)=ipt
3368  nft(i,4)=ipt
3369  nft(i,5)=0
3370  nft(i,6)=0
3371  nft(i,7)=0
3372  nft(i,8)=0
3373  nft(i,9)=0
3374  nft(i,10)=0
3375  nft(i,11)=0
3376  ntj(i)=0
3377  IF(i.GT.abs(ihnt2(4))) nft(i,3)=2112
3378  CALL attflv(nft(i,3),idq,idqq)
3379  nft(i,1)=idq
3380  nft(i,2)=idqq
3381  nft(i,15)=1
3382  IF(abs(idq).GT.1000.OR.(abs(idq*idqq).LT.100.AND.
3383  & rlu(0).LT.0.5)) nft(i,15)=-1
3384  pt(i,14)=ulmass(idq)
3385  pt(i,15)=ulmass(idqq)
3386 200 CONTINUE
3387  RETURN
3388  END
3389 C
3390 C
3391 C
3392  SUBROUTINE attflv(ID,IDQ,IDQQ)
3393  common/ranseed/nseed
3394  SAVE /ranseed/
3395 C
3396  IF(abs(id).LT.100) THEN
3397  nsign=1
3398  idq=id/100
3399  idqq=-id/10+idq*10
3400  IF(abs(idq).EQ.3) nsign=-1
3401  idq=nsign*idq
3402  idqq=nsign*idqq
3403  IF(idq.LT.0) THEN
3404  id0=idq
3405  idq=idqq
3406  idqq=id0
3407  ENDIF
3408  RETURN
3409  ENDIF
3410 C ********return ID of quark(IDQ) and anti-quark(IDQQ)
3411 C for pions and kaons
3412 c
3413 C Return LU ID for quarks and diquarks for proton(ID=2212)
3414 C anti-proton(ID=-2212) and nuetron(ID=2112)
3415 C LU ID for d=1,u=2, (ud)0=2101, (ud)1=2103,
3416 C (dd)1=1103,(uu)1=2203.
3417 C Use SU(6) weight proton=1/3d(uu)1 + 1/6u(ud)1 + 1/2u(ud)0
3418 C nurtron=1/3u(dd)1 + 1/6d(ud)1 + 1/2d(ud)0
3419 C
3420  idq=2
3421  IF(abs(id).EQ.2112) idq=1
3422  idqq=2101
3423  x=rlu(0)
3424  IF(x.LE.0.5) go to 30
3425  IF(x.GT.0.666667) go to 10
3426  idqq=2103
3427  go to 30
3428 10 idq=1
3429  idqq=2203
3430  IF(abs(id).EQ.2112) THEN
3431  idq=2
3432  idqq=1103
3433  ENDIF
3434 30 IF(id.LT.0) THEN
3435  id00=idqq
3436  idqq=-idq
3437  idq=-id00
3438  ENDIF
3439  RETURN
3440  END
3441 C
3442 C*******************************************************************
3443 C This subroutine performs elastic scatterings and possible
3444 C elastic cascading within their own nuclei
3445 c*******************************************************************
3446  SUBROUTINE hijcsc(JP,JT)
3447  dimension psc1(5),psc2(5)
3448  common/hijcrdn/yp(3,300),yt(3,300)
3449  SAVE /hijcrdn/
3450  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3451  SAVE /hiparnt/
3452  common/ranseed/nseed
3453  SAVE /ranseed/
3454  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3455  SAVE /histrng/
3456  IF(jp.EQ.0 .OR. jt.EQ.0) go to 25
3457  DO 10 i=1,5
3458  psc1(i)=pp(jp,i)
3459  psc2(i)=pt(jt,i)
3460 10 CONTINUE
3461  CALL hijels(psc1,psc2)
3462  dpp1=psc1(1)-pp(jp,1)
3463  dpp2=psc1(2)-pp(jp,2)
3464  dpt1=psc2(1)-pt(jt,1)
3465  dpt2=psc2(2)-pt(jt,2)
3466  pp(jp,6)=pp(jp,6)+dpp1/2.0
3467  pp(jp,7)=pp(jp,7)+dpp2/2.0
3468  pp(jp,8)=pp(jp,8)+dpp1/2.0
3469  pp(jp,9)=pp(jp,9)+dpp2/2.0
3470  pt(jt,6)=pt(jt,6)+dpt1/2.0
3471  pt(jt,7)=pt(jt,7)+dpt2/2.0
3472  pt(jt,8)=pt(jt,8)+dpt1/2.0
3473  pt(jt,9)=pt(jt,9)+dpt2/2.0
3474  DO 20 i=1,4
3475  pp(jp,i)=psc1(i)
3476  pt(jt,i)=psc2(i)
3477 20 CONTINUE
3478  nfp(jp,5)=max(1,nfp(jp,5))
3479  nft(jt,5)=max(1,nft(jt,5))
3480 C ********Perform elastic scattering between JP and JT
3481  RETURN
3482 C ********The following is for possible elastic cascade
3483 c
3484 25 IF(jp.EQ.0) go to 45
3485  pabs=sqrt(pp(jp,1)**2+pp(jp,2)**2+pp(jp,3)**2)
3486  bx=pp(jp,1)/pabs
3487  by=pp(jp,2)/pabs
3488  bz=pp(jp,3)/pabs
3489  DO 40 i=1,ihnt2(1)
3490  IF(i.EQ.jp) go to 40
3491  dx=yp(1,i)-yp(1,jp)
3492  dy=yp(2,i)-yp(2,jp)
3493  dz=yp(3,i)-yp(3,jp)
3494  dis=dx*bx+dy*by+dz*bz
3495  IF(dis.LE.0) go to 40
3496  bb=dx**2+dy**2+dz**2-dis**2
3497  r2=bb*hipr1(40)/hipr1(31)/0.1
3498 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
3499  gs=1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3500  & *romg(r2))**2
3501  gs0=1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3502  & *romg(0.0))**2
3503  IF(rlu(0).GT.gs/gs0) go to 40
3504  DO 30 k=1,5
3505  psc1(k)=pp(jp,k)
3506  psc2(k)=pp(i,k)
3507 30 CONTINUE
3508  CALL hijels(psc1,psc2)
3509  dpp1=psc1(1)-pp(jp,1)
3510  dpp2=psc1(2)-pp(jp,2)
3511  dpt1=psc2(1)-pp(i,1)
3512  dpt2=psc2(2)-pp(i,2)
3513  pp(jp,6)=pp(jp,6)+dpp1/2.0
3514  pp(jp,7)=pp(jp,7)+dpp2/2.0
3515  pp(jp,8)=pp(jp,8)+dpp1/2.0
3516  pp(jp,9)=pp(jp,9)+dpp2/2.0
3517  pp(i,6)=pp(i,6)+dpt1/2.0
3518  pp(i,7)=pp(i,7)+dpt2/2.0
3519  pp(i,8)=pp(i,8)+dpt1/2.0
3520  pp(i,9)=pp(i,9)+dpt2/2.0
3521  DO 35 k=1,5
3522  pp(jp,k)=psc1(k)
3523  pp(i,k)=psc2(k)
3524 35 CONTINUE
3525  nfp(i,5)=max(1,nfp(i,5))
3526  go to 45
3527 40 CONTINUE
3528 45 IF(jt.EQ.0) go to 80
3529 50 pabs=sqrt(pt(jt,1)**2+pt(jt,2)**2+pt(jt,3)**2)
3530  bx=pt(jt,1)/pabs
3531  by=pt(jt,2)/pabs
3532  bz=pt(jt,3)/pabs
3533  DO 70 i=1,ihnt2(3)
3534  IF(i.EQ.jt) go to 70
3535  dx=yt(1,i)-yt(1,jt)
3536  dy=yt(2,i)-yt(2,jt)
3537  dz=yt(3,i)-yt(3,jt)
3538  dis=dx*bx+dy*by+dz*bz
3539  IF(dis.LE.0) go to 70
3540  bb=dx**2+dy**2+dz**2-dis**2
3541  r2=bb*hipr1(40)/hipr1(31)/0.1
3542 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
3543  gs=(1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3544  & *romg(r2)))**2
3545  gs0=(1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
3546  & *romg(0.0)))**2
3547  IF(rlu(0).GT.gs/gs0) go to 70
3548  DO 60 k=1,5
3549  psc1(k)=pt(jt,k)
3550  psc2(k)=pt(i,k)
3551 60 CONTINUE
3552  CALL hijels(psc1,psc2)
3553  dpp1=psc1(1)-pt(jt,1)
3554  dpp2=psc1(2)-pt(jt,2)
3555  dpt1=psc2(1)-pt(i,1)
3556  dpt2=psc2(2)-pt(i,2)
3557  pt(jt,6)=pt(jt,6)+dpp1/2.0
3558  pt(jt,7)=pt(jt,7)+dpp2/2.0
3559  pt(jt,8)=pt(jt,8)+dpp1/2.0
3560  pt(jt,9)=pt(jt,9)+dpp2/2.0
3561  pt(i,6)=pt(i,6)+dpt1/2.0
3562  pt(i,7)=pt(i,7)+dpt2/2.0
3563  pt(i,8)=pt(i,8)+dpt1/2.0
3564  pt(i,9)=pt(i,9)+dpt2/2.0
3565  DO 65 k=1,5
3566  pt(jt,k)=psc1(k)
3567  pt(i,k)=psc2(k)
3568 65 CONTINUE
3569  nft(i,5)=max(1,nft(i,5))
3570  go to 80
3571 70 CONTINUE
3572 80 RETURN
3573  END
3574 C
3575 C
3576 C*******************************************************************
3577 CThis subroutine performs elastic scattering between two nucleons
3578 C
3579 C*******************************************************************
3580  SUBROUTINE hijels(PSC1,PSC2)
3581  IMPLICIT DOUBLE PRECISION(d)
3582  dimension psc1(5),psc2(5)
3583  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3584  SAVE /hiparnt/
3585  common/ranseed/nseed
3586  SAVE /ranseed/
3587 C
3588  cc=1.0-hint1(12)/hint1(13)
3589  rr=(1.0-cc)*hint1(13)/hint1(12)/(1.0-hipr1(33))-1.0
3590  bb=0.5*(3.0+rr+sqrt(9.0+10.0*rr+rr**2))
3591  ep=sqrt((psc1(1)-psc2(1))**2+(psc1(2)-psc2(2))**2
3592  & +(psc1(3)-psc2(3))**2)
3593  IF(ep.LE.0.1) RETURN
3594  els0=98.0/ep+52.0*(1.0+rr)**2
3595  pcm1=psc1(1)+psc2(1)
3596  pcm2=psc1(2)+psc2(2)
3597  pcm3=psc1(3)+psc2(3)
3598  ecm=psc1(4)+psc2(4)
3599  am1=psc1(5)**2
3600  am2=psc2(5)**2
3601  amm=ecm**2-pcm1**2-pcm2**2-pcm3**2
3602  IF(amm.LE.psc1(5)+psc2(5)) RETURN
3603 C ********elastic scattering only when approaching
3604 C to each other
3605  pmax=(amm**2+am1**2+am2**2-2.0*amm*am1-2.0*amm*am2
3606  & -2.0*am1*am2)/4.0/amm
3607  pmax=abs(pmax)
3608 20 tt=rlu(0)*min(pmax,1.5)
3609  els=98.0*exp(-2.8*tt)/ep
3610  & +52.0*exp(-9.2*tt)*(1.0+rr*exp(-4.6*(bb-1.0)*tt))**2
3611  IF(rlu(0).GT.els/els0) go to 20
3612  phi=2.0*hipr1(40)*rlu(0)
3613 C
3614  dbx=pcm1/ecm
3615  dby=pcm2/ecm
3616  dbz=pcm3/ecm
3617  db=sqrt(dbx**2+dby**2+dbz**2)
3618  IF(db.GT.0.99999999d0) THEN
3619  dbx=dbx*(0.99999999d0/db)
3620  dby=dby*(0.99999999d0/db)
3621  dbz=dbz*(0.99999999d0/db)
3622  db=0.99999999d0
3623  WRITE(6,*) ' (HIJELS) boost vector too large'
3624 C ********Rescale boost vector if too close to unity.
3625  ENDIF
3626  dga=1d0/sqrt(1d0-db**2)
3627 C
3628  dp1=sqrt(tt)*sin(phi)
3629  dp2=sqrt(tt)*cos(phi)
3630  dp3=sqrt(pmax-tt)
3631  dp4=sqrt(pmax+am1)
3632  dbp=dbx*dp1+dby*dp2+dbz*dp3
3633  dgabp=dga*(dga*dbp/(1d0+dga)+dp4)
3634  psc1(1)=dp1+dgabp*dbx
3635  psc1(2)=dp2+dgabp*dby
3636  psc1(3)=dp3+dgabp*dbz
3637  psc1(4)=dga*(dp4+dbp)
3638 C
3639  dp1=-sqrt(tt)*sin(phi)
3640  dp2=-sqrt(tt)*cos(phi)
3641  dp3=-sqrt(pmax-tt)
3642  dp4=sqrt(pmax+am2)
3643  dbp=dbx*dp1+dby*dp2+dbz*dp3
3644  dgabp=dga*(dga*dbp/(1d0+dga)+dp4)
3645  psc2(1)=dp1+dgabp*dbx
3646  psc2(2)=dp2+dgabp*dby
3647  psc2(3)=dp3+dgabp*dbz
3648  psc2(4)=dga*(dp4+dbp)
3649  RETURN
3650  END
3651 C
3652 C
3653 C*******************************************************************
3654 C *
3655 C Subroutine HIJSFT *
3656 C *
3657 C Scatter two excited strings, JP from proj and JT from target *
3658 C*******************************************************************
3659  SUBROUTINE hijsft(JP,JT,JOUT,IERROR)
3660  common/hijcrdn/yp(3,300),yt(3,300)
3661  SAVE /hijcrdn/
3662  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
3663  SAVE /hiparnt/
3664  common/hijdat/hidat0(10,10),hidat(10)
3665  SAVE /hijdat/
3666  common/ranseed/nseed
3667  SAVE /ranseed/
3668  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
3669  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
3670  & pjpm(300,500),ntj(300),kftj(300,500),
3671  & pjtx(300,500),pjty(300,500),pjtz(300,500),
3672  & pjte(300,500),pjtm(300,500)
3673  SAVE /hijjet1/
3674  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
3675  & k2sg(900,100),pxsg(900,100),pysg(900,100),
3676  & pzsg(900,100),pesg(900,100),pmsg(900,100)
3677  SAVE /hijjet2/
3678  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
3679  SAVE /histrng/
3680  common/dpmcom1/jjp,jjt,amp,amt,apx0,atx0,ampn,amtn,amp0,amt0,
3681  & nfdp,nfdt,wp,wm,sw,xremp,xremt,dpkc1,dpkc2,pp11,pp12,
3682  & pt11,pt12,ptp2,ptt2
3683  SAVE /dpmcom1/
3684  common/dpmcom2/ndpm,kdpm(20,2),pdpm1(20,5),pdpm2(20,5)
3685  SAVE /dpmcom2/
3686 C*******************************************************************
3687 C JOUT-> the number
3688 C of hard scatterings preceding this soft collision.
3689 C IHNT2(13)-> 1=
3690 C double diffrac 2=single diffrac, 3=non-single diffrac.
3691 C*******************************************************************
3692  ierror=0
3693  jjp=jp
3694  jjt=jt
3695  ndpm=0
3696  iopmain=0
3697  IF(jp.GT.ihnt2(1) .OR. jt.GT.ihnt2(3)) RETURN
3698 
3699  epp=pp(jp,4)+pp(jp,3)
3700  epm=pp(jp,4)-pp(jp,3)
3701  etp=pt(jt,4)+pt(jt,3)
3702  etm=pt(jt,4)-pt(jt,3)
3703 
3704  wp=epp+etp
3705  wm=epm+etm
3706  sw=wp*wm
3707 C ********total W+,W- and center-of-mass energy
3708 
3709  IF(wp.LT.0.0 .OR. wm.LT.0.0) go to 1000
3710 
3711  IF(jout.EQ.0) THEN
3712  IF(epp.LT.0.0) go to 1000
3713  IF(epm.LT.0.0) go to 1000
3714  IF(etp.LT.0.0) go to 1000
3715  IF(etm.LT.0.0) go to 1000
3716  IF(epp/(epm+0.01).LE.etp/(etm+0.01)) RETURN
3717  ENDIF
3718 C ********For strings which does not follow a jet-prod,
3719 C scatter only if Ycm(JP)>Ycm(JT). When jets
3720 C are produced just before this collision
3721 C this requirement has already be enforced
3722 C (see SUBROUTINE HIJHRD)
3723  ihnt2(11)=jp
3724  ihnt2(12)=jt
3725 C
3726 C
3727 C
3728  miss=0
3729  pkc1=0.0
3730  pkc2=0.0
3731  pkc11=0.0
3732  pkc12=0.0
3733  pkc21=0.0
3734  pkc22=0.0
3735  dpkc11=0.0
3736  dpkc12=0.0
3737  dpkc21=0.0
3738  dpkc22=0.0
3739  IF(nfp(jp,10).EQ.1.OR.nft(jt,10).EQ.1) THEN
3740  IF(nfp(jp,10).EQ.1) THEN
3741  phi1=ulangl(pp(jp,10),pp(jp,11))
3742  ppjet=sqrt(pp(jp,10)**2+pp(jp,11)**2)
3743  pkc1=ppjet
3744  pkc11=pp(jp,10)
3745  pkc12=pp(jp,11)
3746  ENDIF
3747  IF(nft(jt,10).EQ.1) THEN
3748  phi2=ulangl(pt(jt,10),pt(jt,11))
3749  ptjet=sqrt(pt(jt,10)**2+pt(jt,11)**2)
3750  pkc2=ptjet
3751  pkc21=pt(jt,10)
3752  pkc22=pt(jt,11)
3753  ENDIF
3754  IF(ihpr2(4).GT.0.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1) THEN
3755  IF(nfp(jp,10).EQ.0) THEN
3756  phi=-phi2
3757  ELSE IF(nft(jt,10).EQ.0) THEN
3758  phi=phi1
3759  ELSE
3760  phi=(phi1+phi2-hipr1(40))/2.0
3761  ENDIF
3762  bx=hint1(19)*cos(hint1(20))
3763  by=hint1(19)*sin(hint1(20))
3764  xp0=yp(1,jp)
3765  yp0=yp(2,jp)
3766  xt0=yt(1,jt)+bx
3767  yt0=yt(2,jt)+by
3768  r1=max(1.2*ihnt2(1)**0.3333333,
3769  & sqrt(xp0**2+yp0**2))
3770  r2=max(1.2*ihnt2(3)**0.3333333,
3771  & sqrt((xt0-bx)**2+(yt0-by)**2))
3772  IF(abs(cos(phi)).LT.1.0e-5) THEN
3773  dd1=r1
3774  dd2=r1
3775  dd3=abs(by+sqrt(r2**2-(xp0-bx)**2)-yp0)
3776  dd4=abs(by-sqrt(r2**2-(xp0-bx)**2)-yp0)
3777  go to 5
3778  ENDIF
3779  bb=2.0*sin(phi)*(cos(phi)*yp0-sin(phi)*xp0)
3780  cc=(yp0**2-r1**2)*cos(phi)**2+xp0*sin(phi)*(
3781  & xp0*sin(phi)-2.0*yp0*cos(phi))
3782  dd=bb**2-4.0*cc
3783  IF(dd.LT.0.0) go to 10
3784  xx1=(-bb+sqrt(dd))/2.0
3785  xx2=(-bb-sqrt(dd))/2.0
3786  dd1=abs((xx1-xp0)/cos(phi))
3787  dd2=abs((xx2-xp0)/cos(phi))
3788 C
3789  bb=2.0*sin(phi)*(cos(phi)*(yt0-by)-sin(phi)*xt0)-2.0*bx
3790  cc=(bx**2+(yt0-by)**2-r2**2)*cos(phi)**2+xt0*sin(phi)
3791  & *(xt0*sin(phi)-2.0*cos(phi)*(yt0-by))
3792  & -2.0*bx*sin(phi)*(cos(phi)*(yt0-by)-sin(phi)*xt0)
3793  dd=bb**2-4.0*cc
3794  IF(dd.LT.0.0) go to 10
3795  xx1=(-bb+sqrt(dd))/2.0
3796  xx2=(-bb-sqrt(dd))/2.0
3797  dd3=abs((xx1-xt0)/cos(phi))
3798  dd4=abs((xx2-xt0)/cos(phi))
3799 C
3800  5 dd1=min(dd1,dd3)
3801  dd2=min(dd2,dd4)
3802  IF(dd1.LT.hipr1(13)) dd1=0.0
3803  IF(dd2.LT.hipr1(13)) dd2=0.0
3804  IF(nfp(jp,10).EQ.1.AND.ppjet.GT.hipr1(11)) THEN
3805  dp1=dd1*hipr1(14)/2.0
3806  dp1=min(dp1,ppjet-hipr1(11))
3807  pkc1=ppjet-dp1
3808  dpx1=cos(phi1)*dp1
3809  dpy1=sin(phi1)*dp1
3810  pkc11=pp(jp,10)-dpx1
3811  pkc12=pp(jp,11)-dpy1
3812  IF(dp1.GT.0.0) THEN
3813  cthep=pp(jp,12)/sqrt(pp(jp,12)**2+ppjet**2)
3814  dpz1=dp1*cthep/sqrt(1.0-cthep**2)
3815  dpe1=sqrt(dpx1**2+dpy1**2+dpz1**2)
3816  eppprm=pp(jp,4)+pp(jp,3)-dpe1-dpz1
3817  epmprm=pp(jp,4)-pp(jp,3)-dpe1+dpz1
3818  IF(eppprm.LE.0.0.OR.epmprm.LE.0.0) go to 15
3819  epp=eppprm
3820  epm=epmprm
3821  pp(jp,10)=pkc11
3822  pp(jp,11)=pkc12
3823  npj(jp)=npj(jp)+1
3824  kfpj(jp,npj(jp))=21
3825  pjpx(jp,npj(jp))=dpx1
3826  pjpy(jp,npj(jp))=dpy1
3827  pjpz(jp,npj(jp))=dpz1
3828  pjpe(jp,npj(jp))=dpe1
3829  pjpm(jp,npj(jp))=0.0
3830  pp(jp,3)=pp(jp,3)-dpz1
3831  pp(jp,4)=pp(jp,4)-dpe1
3832  ENDIF
3833  ENDIF
3834  15 IF(nft(jt,10).EQ.1.AND.ptjet.GT.hipr1(11)) THEN
3835  dp2=dd2*hipr1(14)/2.0
3836  dp2=min(dp2,ptjet-hipr1(11))
3837  pkc2=ptjet-dp2
3838  dpx2=cos(phi2)*dp2
3839  dpy2=sin(phi2)*dp2
3840  pkc21=pt(jt,10)-dpx2
3841  pkc22=pt(jt,11)-dpy2
3842  IF(dp2.GT.0.0) THEN
3843  cthet=pt(jt,12)/sqrt(pt(jt,12)**2+ptjet**2)
3844  dpz2=dp2*cthet/sqrt(1.0-cthet**2)
3845  dpe2=sqrt(dpx2**2+dpy2**2+dpz2**2)
3846  etpprm=pt(jt,4)+pt(jt,3)-dpe2-dpz2
3847  etmprm=pt(jt,4)-pt(jt,3)-dpe2+dpz2
3848  IF(etpprm.LE.0.0.OR.etmprm.LE.0.0) go to 16
3849  etp=etpprm
3850  etm=etmprm
3851  pt(jt,10)=pkc21
3852  pt(jt,11)=pkc22
3853  ntj(jt)=ntj(jt)+1
3854  kftj(jt,ntj(jt))=21
3855  pjtx(jt,ntj(jt))=dpx2
3856  pjty(jt,ntj(jt))=dpy2
3857  pjtz(jt,ntj(jt))=dpz2
3858  pjte(jt,ntj(jt))=dpe2
3859  pjtm(jt,ntj(jt))=0.0
3860  pt(jt,3)=pt(jt,3)-dpz2
3861  pt(jt,4)=pt(jt,4)-dpe2
3862  ENDIF
3863  ENDIF
3864  16 dpkc11=-(pp(jp,10)-pkc11)/2.0
3865  dpkc12=-(pp(jp,11)-pkc12)/2.0
3866  dpkc21=-(pt(jt,10)-pkc21)/2.0
3867  dpkc22=-(pt(jt,11)-pkc22)/2.0
3868  wp=epp+etp
3869  wm=epm+etm
3870  sw=wp*wm
3871  ENDIF
3872  ENDIF
3873 C ********If jet is quenched the pt from valence quark
3874 C hard scattering has to reduced by d*kapa
3875 C
3876 C
3877 10 ptp02=pp(jp,1)**2+pp(jp,2)**2
3878  ptt02=pt(jt,1)**2+pt(jt,2)**2
3879 C
3880  amq=max(pp(jp,14)+pp(jp,15),pt(jt,14)+pt(jt,15))
3881  amx=hipr1(1)+amq
3882 C ********consider mass cut-off for strings which
3883 C must also include quark's mass
3884  amp0=amx
3885  dpm0=amx
3886  nfdp=0
3887  IF(nfp(jp,5).LE.2.AND.nfp(jp,3).NE.0) THEN
3888  amp0=ulmass(nfp(jp,3))
3889  nfdp=nfp(jp,3)+2*nfp(jp,3)/abs(nfp(jp,3))
3890  dpm0=ulmass(nfdp)
3891  IF(dpm0.LE.0.0) THEN
3892  nfdp=nfdp-2*nfdp/abs(nfdp)
3893  dpm0=ulmass(nfdp)
3894  ENDIF
3895  ENDIF
3896  amt0=amx
3897  dtm0=amx
3898  nfdt=0
3899  IF(nft(jt,5).LE.2.AND.nft(jt,3).NE.0) THEN
3900  amt0=ulmass(nft(jt,3))
3901  nfdt=nft(jt,3)+2*nft(jt,3)/abs(nft(jt,3))
3902  dtm0=ulmass(nfdt)
3903  IF(dtm0.LE.0.0) THEN
3904  nfdt=nfdt-2*nfdt/abs(nfdt)
3905  dtm0=ulmass(nfdt)
3906  ENDIF
3907  ENDIF
3908 C
3909  ampn=sqrt(amp0**2+ptp02)
3910  amtn=sqrt(amt0**2+ptt02)
3911  snn=(ampn+amtn)**2+0.001
3912 C
3913  IF(sw.LT.snn+0.001) go to 4000
3914 C ********Scatter only if SW>SNN
3915 C*****give some PT kick to the two exited strings******************
3916 20 swptn=4.0*(max(amp0,amt0)**2+max(ptp02,ptt02))
3917  swptd=4.0*(max(dpm0,dtm0)**2+max(ptp02,ptt02))
3918  swptx=4.0*(amx**2+max(ptp02,ptt02))
3919  IF(sw.LE.swptn) THEN
3920  pkcmx=0.0
3921  ELSE IF(sw.GT.swptn .AND. sw.LE.swptd
3922  & .AND.npj(jp).EQ.0.AND.ntj(jt).EQ.0) THEN
3923  pkcmx=sqrt(sw/4.0-max(amp0,amt0)**2)
3924  & -sqrt(max(ptp02,ptt02))
3925  ELSE IF(sw.GT.swptd .AND. sw.LE.swptx
3926  & .AND.npj(jp).EQ.0.AND.ntj(jt).EQ.0) THEN
3927  pkcmx=sqrt(sw/4.0-max(dpm0,dtm0)**2)
3928  & -sqrt(max(ptp02,ptt02))
3929  ELSE IF(sw.GT.swptx) THEN
3930  pkcmx=sqrt(sw/4.0-amx**2)-sqrt(max(ptp02,ptt02))
3931  ENDIF
3932 C ********maximun PT kick
3933 C*********************************************************
3934 C
3935  IF(nfp(jp,10).EQ.1.OR.nft(jt,10).EQ.1) THEN
3936  IF(pkc1.GT.pkcmx) THEN
3937  pkc1=pkcmx
3938  pkc11=pkc1*cos(phi1)
3939  pkc12=pkc1*sin(phi1)
3940  dpkc11=-(pp(jp,10)-pkc11)/2.0
3941  dpkc12=-(pp(jp,11)-pkc12)/2.0
3942  ENDIF
3943  IF(pkc2.GT.pkcmx) THEN
3944  pkc2=pkcmx
3945  pkc21=pkc2*cos(phi2)
3946  pkc22=pkc2*sin(phi2)
3947  dpkc21=-(pt(jt,10)-pkc21)/2.0
3948  dpkc22=-(pt(jt,11)-pkc22)/2.0
3949  ENDIF
3950  dpkc1=dpkc11+dpkc21
3951  dpkc2=dpkc12+dpkc22
3952  nfp(jp,10)=-nfp(jp,10)
3953  nft(jt,10)=-nft(jt,10)
3954  go to 40
3955  ENDIF
3956 C ********If the valence quarks had a hard-collision
3957 C the pt kick is the pt from hard-collision.
3958  i_sng=0
3959  IF(ihpr2(13).NE.0 .AND. rlu(0).LE.hidat(4)) i_sng=1
3960  IF((nfp(jp,5).EQ.3 .OR.nft(jt,5).EQ.3).OR.
3961  & (npj(jp).NE.0.OR.nfp(jp,10).NE.0).OR.
3962  & (ntj(jt).NE.0.OR.nft(jt,10).NE.0)) i_sng=0
3963 C
3964 C ********decite whether to have single-diffractive
3965  IF(ihpr2(5).EQ.0) THEN
3966  pkc=hipr1(2)*sqrt(-alog(1.0-rlu(0)
3967  & *(1.0-exp(-pkcmx**2/hipr1(2)**2))))
3968  go to 30
3969  ENDIF
3970  pkc=hirnd2(3,0.0,pkcmx**2)
3971  pkc=sqrt(pkc)
3972  IF(pkc.GT.hipr1(20))
3973  & pkc=hipr1(2)*sqrt(-alog(exp(-hipr1(20)**2/hipr1(2)**2)
3974  & -rlu(0)*(exp(-hipr1(20)**2/hipr1(2)**2)-
3975  & exp(-pkcmx**2/hipr1(2)**2))))
3976 C
3977  IF(i_sng.EQ.1) pkc=0.65*sqrt(
3978  & -alog(1.0-rlu(0)*(1.0-exp(-pkcmx**2/0.65**2))))
3979 C ********select PT kick
3980 30 phi0=2.0*hipr1(40)*rlu(0)
3981  pkc11=pkc*sin(phi0)
3982  pkc12=pkc*cos(phi0)
3983  pkc21=-pkc11
3984  pkc22=-pkc12
3985  dpkc1=0.0
3986  dpkc2=0.0
3987 40 pp11=pp(jp,1)+pkc11-dpkc1
3988  pp12=pp(jp,2)+pkc12-dpkc2
3989  pt11=pt(jt,1)+pkc21-dpkc1
3990  pt12=pt(jt,2)+pkc22-dpkc2
3991  ptp2=pp11**2+pp12**2
3992  ptt2=pt11**2+pt12**2
3993 C
3994  ampn=sqrt(amp0**2+ptp2)
3995  amtn=sqrt(amt0**2+ptt2)
3996  snn=(ampn+amtn)**2+0.001
3997 C***************************************
3998  wp=epp+etp
3999  wm=epm+etm
4000  sw=wp*wm
4001 C****************************************
4002  IF(sw.LT.snn) THEN
4003  miss=miss+1
4004  IF(miss.LE.100) then
4005  pkc=0.0
4006  go to 30
4007  ENDIF
4008  IF(ihpr2(10).NE.0)
4009  & WRITE(6,*) 'Error occured in Pt kick section of HIJSFT'
4010  go to 4000
4011  ENDIF
4012 C******************************************************************
4013  ampd=sqrt(dpm0**2+ptp2)
4014  amtd=sqrt(dtm0**2+ptt2)
4015 
4016  ampx=sqrt(amx**2+ptp2)
4017  amtx=sqrt(amx**2+ptt2)
4018 
4019  dpn=ampn**2/sw
4020  dtn=amtn**2/sw
4021  dpd=ampd**2/sw
4022  dtd=amtd**2/sw
4023  dpx=ampx**2/sw
4024  dtx=amtx**2/sw
4025 C
4026  spntd=(ampn+amtd)**2
4027  spntx=(ampn+amtx)**2
4028 C ********CM energy if proj=N,targ=N*
4029  spdtn=(ampd+amtn)**2
4030  spxtn=(ampx+amtn)**2
4031 C ********CM energy if proj=N*,targ=N
4032  spdtx=(ampd+amtx)**2
4033  spxtd=(ampx+amtd)**2
4034  sdd=(ampd+amtd)**2
4035  sxx=(ampx+amtx)**2
4036 
4037 C
4038 C
4039 C ********CM energy if proj=delta, targ=delta
4040 C****************There are many different cases**********
4041 c IF(IHPR2(15).EQ.1) GO TO 500
4042 C
4043 C ********to have DPM type soft interactions
4044 C
4045  45 CONTINUE
4046  IF(sw.GT.sxx+0.001) THEN
4047  IF(i_sng.EQ.0) THEN
4048  d1=dpx
4049  d2=dtx
4050  nfp3=0
4051  nft3=0
4052  go to 400
4053  ELSE
4054 c**** 5/30/1998 this is identical to the above statement. Added to
4055 c**** avoid questional branching to block.
4056  IF((nfp(jp,5).EQ.3 .AND.nft(jt,5).EQ.3).OR.
4057  & (npj(jp).NE.0.OR.nfp(jp,10).NE.0).OR.
4058  & (ntj(jt).NE.0.OR.nft(jt,10).NE.0)) THEN
4059  d1=dpx
4060  d2=dtx
4061  nfp3=0
4062  nft3=0
4063  go to 400
4064  ENDIF
4065 C ********do not allow excited strings to have
4066 C single-diffr
4067  IF(rlu(0).GT.0.5.OR.(nft(jt,5).GT.2.OR.
4068  & ntj(jt).NE.0.OR.nft(jt,10).NE.0)) THEN
4069  d1=dpn
4070  d2=dtx
4071  nfp3=nfp(jp,3)
4072  nft3=0
4073  go to 220
4074  ELSE
4075  d1=dpx
4076  d2=dtn
4077  nfp3=0
4078  nft3=nft(jt,3)
4079  go to 240
4080  ENDIF
4081 C ********have single diffractive collision
4082  ENDIF
4083  ELSE IF(sw.GT.max(spdtx,spxtd)+0.001 .AND.
4084  & sw.LE.sxx+0.001) THEN
4085  IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0.AND.
4086  & rlu(0).GT.0.5).OR.(npj(jp).EQ.0
4087  & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2) THEN
4088  d1=dpd
4089  d2=dtx
4090  nfp3=nfdp
4091  nft3=0
4092  go to 220
4093  ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2) THEN
4094  d1=dpx
4095  d2=dtd
4096  nfp3=0
4097  nft3=nfdt
4098  go to 240
4099  ENDIF
4100  go to 4000
4101  ELSE IF(sw.GT.min(spdtx,spxtd)+0.001.AND.
4102  & sw.LE.max(spdtx,spxtd)+0.001) THEN
4103  IF(spdtx.LE.spxtd.AND.npj(jp).EQ.0
4104  & .AND.nfp(jp,5).LE.2) THEN
4105  d1=dpd
4106  d2=dtx
4107  nfp3=nfdp
4108  nft3=0
4109  go to 220
4110  ELSE IF(spdtx.GT.spxtd.AND.ntj(jt).EQ.0
4111  & .AND.nft(jt,5).LE.2) THEN
4112  d1=dpx
4113  d2=dtd
4114  nfp3=0
4115  nft3=nfdt
4116  go to 240
4117  ENDIF
4118 c*** 5/30/1998 added to avoid questional branching to another block
4119 c*** this is identical to the statement following the next ELSE IF
4120  IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0
4121  & .AND.rlu(0).GT.0.5).OR.(npj(jp).EQ.0
4122  & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2) THEN
4123  d1=dpn
4124  d2=dtx
4125  nfp3=nfp(jp,3)
4126  nft3=0
4127  go to 220
4128  ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2) THEN
4129  d1=dpx
4130  d2=dtn
4131  nfp3=0
4132  nft3=nft(jt,3)
4133  go to 240
4134  ENDIF
4135  go to 4000
4136  ELSE IF(sw.GT.max(spntx,spxtn)+0.001 .AND.
4137  & sw.LE.min(spdtx,spxtd)+0.001) THEN
4138  IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0
4139  & .AND.rlu(0).GT.0.5).OR.(npj(jp).EQ.0
4140  & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2) THEN
4141  d1=dpn
4142  d2=dtx
4143  nfp3=nfp(jp,3)
4144  nft3=0
4145  go to 220
4146  ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2) THEN
4147  d1=dpx
4148  d2=dtn
4149  nfp3=0
4150  nft3=nft(jt,3)
4151  go to 240
4152  ENDIF
4153  go to 4000
4154  ELSE IF(sw.GT.min(spntx,spxtn)+0.001 .AND.
4155  & sw.LE.max(spntx,spxtn)+0.001) THEN
4156  IF(spntx.LE.spxtn.AND.npj(jp).EQ.0
4157  & .AND.nfp(jp,5).LE.2) THEN
4158  d1=dpn
4159  d2=dtx
4160  nfp3=nfp(jp,3)
4161  nft3=0
4162  go to 220
4163  ELSEIF(spntx.GT.spxtn.AND.ntj(jt).EQ.0
4164  & .AND.nft(jt,5).LE.2) THEN
4165  d1=dpx
4166  d2=dtn
4167  nfp3=0
4168  nft3=nft(jt,3)
4169  go to 240
4170  ENDIF
4171  go to 4000
4172  ELSE IF(sw.LE.min(spntx,spxtn)+0.001 .AND.
4173  & (npj(jp).NE.0 .OR.ntj(jt).NE.0)) THEN
4174  go to 4000
4175  ELSE IF(sw.LE.min(spntx,spxtn)+0.001 .AND.
4176  & nfp(jp,5).GT.2.AND.nft(jt,5).GT.2) THEN
4177  go to 4000
4178  ELSE IF(sw.GT.sdd+0.001.AND.sw.LE.
4179  & min(spntx,spxtn)+0.001) THEN
4180  d1=dpd
4181  d2=dtd
4182  nfp3=nfdp
4183  nft3=nfdt
4184  go to 100
4185  ELSE IF(sw.GT.max(spntd,spdtn)+0.001
4186  & .AND. sw.LE.sdd+0.001) THEN
4187  IF(rlu(0).GT.0.5) THEN
4188  d1=dpd
4189  d2=dtn
4190  nfp3=nfdp
4191  nft3=nft(jt,3)
4192  go to 100
4193  ELSE
4194  d1=dpn
4195  d2=dtd
4196  nfp3=nfp(jp,3)
4197  nft3=nfdt
4198  go to 100
4199  ENDIF
4200  ELSE IF(sw.GT.min(spntd,spdtn)+0.001
4201  & .AND. sw.LE.max(spntd,spdtn)+0.001) THEN
4202  IF(spntd.GT.spdtn) THEN
4203  d1=dpd
4204  d2=dtn
4205  nfp3=nfdp
4206  nft3=nft(jt,3)
4207  go to 100
4208  ELSE
4209  d1=dpn
4210  d2=dtd
4211  nfp3=nfp(jp,3)
4212  nft3=nfdt
4213  go to 100
4214  ENDIF
4215  ELSE IF(sw.LE.min(spntd,spdtn)+0.001) THEN
4216  d1=dpn
4217  d2=dtn
4218  nfp3=nfp(jp,3)
4219  nft3=nft(jt,3)
4220  go to 100
4221  ENDIF
4222  WRITE(6,*) ' Error in HIJSFT: There is no path to here'
4223  RETURN
4224 C
4225 C*************** elastic scattering ***************
4226 C this is like elastic, both proj and targ mass
4227 C must be fixed
4228 C***************************************************
4229 100 nfp5=max(2,nfp(jp,5))
4230  nft5=max(2,nft(jt,5))
4231  bb1=1.0+d1-d2
4232  bb2=1.0+d2-d1
4233  IF(bb1**2.LT.4.0*d1 .OR. bb2**2.LT.4.0*d2) THEN
4234  miss=miss+1
4235  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4236  pkc=pkc*0.5
4237  go to 30
4238  ENDIF
4239  IF(rlu(0).LT.0.5) THEN
4240  x1=(bb1-sqrt(bb1**2-4.0*d1))/2.0
4241  x2=(bb2-sqrt(bb2**2-4.0*d2))/2.0
4242  ELSE
4243  x1=(bb1+sqrt(bb1**2-4.0*d1))/2.0
4244  x2=(bb2+sqrt(bb2**2-4.0*d2))/2.0
4245  ENDIF
4246  ihnt2(13)=2
4247  go to 600
4248 C
4249 C********** Single diffractive ***********************
4250 C either proj or targ's mass is fixed
4251 C*****************************************************
4252 220 nfp5=max(2,nfp(jp,5))
4253  nft5=3
4254  IF(nfp3.EQ.0) nfp5=3
4255  bb2=1.0+d2-d1
4256  IF(bb2**2.LT.4.0*d2) THEN
4257  miss=miss+1
4258  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4259  pkc=pkc*0.5
4260  go to 30
4261  ENDIF
4262  xmin=(bb2-sqrt(bb2**2-4.0*d2))/2.0
4263  xmax=(bb2+sqrt(bb2**2-4.0*d2))/2.0
4264  miss4=0
4265 222 x2=hirnd2(6,xmin,xmax)
4266  x1=d1/(1.0-x2)
4267  IF(x2*(1.0-x1).LT.(d2+1.e-4/sw)) THEN
4268  miss4=miss4+1
4269  IF(miss4.LE.1000) go to 222
4270  go to 5000
4271  ENDIF
4272  ihnt2(13)=2
4273  go to 600
4274 C ********Fix proj mass*********
4275 240 nfp5=3
4276  nft5=max(2,nft(jt,5))
4277  IF(nft3.EQ.0) nft5=3
4278  bb1=1.0+d1-d2
4279  IF(bb1**2.LT.4.0*d1) THEN
4280  miss=miss+1
4281  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4282  pkc=pkc*0.5
4283  go to 30
4284  ENDIF
4285  xmin=(bb1-sqrt(bb1**2-4.0*d1))/2.0
4286  xmax=(bb1+sqrt(bb1**2-4.0*d1))/2.0
4287  miss4=0
4288 242 x1=hirnd2(6,xmin,xmax)
4289  x2=d2/(1.0-x1)
4290  IF(x1*(1.0-x2).LT.(d1+1.e-4/sw)) THEN
4291  miss4=miss4+1
4292  IF(miss4.LE.1000) go to 242
4293  go to 5000
4294  ENDIF
4295  ihnt2(13)=2
4296  go to 600
4297 C ********Fix targ mass*********
4298 C
4299 C*************non-single diffractive**********************
4300 C both proj and targ may not be fixed in mass
4301 C*********************************************************
4302 C
4303 400 nfp5=3
4304  nft5=3
4305  bb1=1.0+d1-d2
4306  bb2=1.0+d2-d1
4307  IF(bb1**2.LT.4.0*d1 .OR. bb2**2.LT.4.0*d2) THEN
4308  miss=miss+1
4309  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
4310  pkc=pkc*0.5
4311  go to 30
4312  ENDIF
4313  xmin1=(bb1-sqrt(bb1**2-4.0*d1))/2.0
4314  xmax1=(bb1+sqrt(bb1**2-4.0*d1))/2.0
4315  xmin2=(bb2-sqrt(bb2**2-4.0*d2))/2.0
4316  xmax2=(bb2+sqrt(bb2**2-4.0*d2))/2.0
4317  miss4=0
4318 410 x1=hirnd2(4,xmin1,xmax1)
4319  x2=hirnd2(4,xmin2,xmax2)
4320  IF(nfp(jp,5).EQ.3.OR.nft(jt,5).EQ.3) THEN
4321  x1=hirnd2(6,xmin1,xmax1)
4322  x2=hirnd2(6,xmin2,xmax2)
4323  ENDIF
4324 C ********
4325  IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000.OR.
4326  & abs(nfp(jp,1)*nfp(jp,2)).LT.100) THEN
4327  x1=hirnd2(5,xmin1,xmax1)
4328  ENDIF
4329  IF(abs(nft(jt,1)*nft(jt,2)).GT.1000000.OR.
4330  & abs(nft(jt,1)*nft(jt,2)).LT.100) THEN
4331  x2=hirnd2(5,xmin2,xmax2)
4332  ENDIF
4333 c IF(IOPMAIN.EQ.3) X1=HIRND2(6,XMIN1,XMAX1)
4334 c IF(IOPMAIN.EQ.2) X2=HIRND2(6,XMIN2,XMAX2)
4335 C ********For q-qbar or (qq)-(qq)bar system use symetric
4336 C distribution, for q-(qq) or qbar-(qq)bar use
4337 C unsymetrical distribution
4338 C
4339  IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000) x1=1.0-x1
4340  xxp=x1*(1.0-x2)
4341  xxt=x2*(1.0-x1)
4342  IF(xxp.LT.(d1+1.e-4/sw) .OR. xxt.LT.(d2+1.e-4/sw)) THEN
4343  miss4=miss4+1
4344  IF(miss4.LE.1000) go to 410
4345  go to 5000
4346  ENDIF
4347  ihnt2(13)=3
4348 C***************************************************
4349 C***************************************************
4350 600 CONTINUE
4351  IF(x1*(1.0-x2).LT.(ampn**2-1.e-4)/sw.OR.
4352  & x2*(1.0-x1).LT.(amtn**2-1.e-4)/sw) THEN
4353  miss=miss+1
4354  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 2000
4355  pkc=0.0
4356  go to 30
4357  ENDIF
4358 C
4359  epp=(1.0-x2)*wp
4360  epm=x1*wm
4361  etp=x2*wp
4362  etm=(1.0-x1)*wm
4363  pp(jp,3)=(epp-epm)/2.0
4364  pp(jp,4)=(epp+epm)/2.0
4365  IF(epp*epm-ptp2.LT.0.0) go to 6000
4366  pp(jp,5)=sqrt(epp*epm-ptp2)
4367  nfp(jp,3)=nfp3
4368  nfp(jp,5)=nfp5
4369 
4370  pt(jt,3)=(etp-etm)/2.0
4371  pt(jt,4)=(etp+etm)/2.0
4372  IF(etp*etm-ptt2.LT.0.0) go to 6000
4373  pt(jt,5)=sqrt(etp*etm-ptt2)
4374  nft(jt,3)=nft3
4375  nft(jt,5)=nft5
4376 C*****recoil PT from hard-inter is shared by two end-partons
4377 C so that pt=p1+p2
4378  pp(jp,1)=pp11-pkc11
4379  pp(jp,2)=pp12-pkc12
4380 
4381  kickdip=1
4382  kickdit=1
4383  IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000.OR.
4384  & abs(nfp(jp,1)*nfp(jp,2)).LT.100) THEN
4385  kickdip=0
4386  ENDIF
4387  IF(abs(nft(jt,1)*nft(jt,2)).GT.1000000.OR.
4388  & abs(nft(jt,1)*nft(jt,2)).LT.100) THEN
4389  kickdit=0
4390  ENDIF
4391  IF((kickdip.EQ.0.AND.rlu(0).LT.0.5)
4392  & .OR.(kickdip.NE.0.AND.rlu(0)
4393  & .LT.0.5/(1.0+(pkc11**2+pkc12**2)/hipr1(22)**2))) THEN
4394  pp(jp,6)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0+pp(jp,6)
4395  pp(jp,7)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0+pp(jp,7)
4396  pp(jp,8)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0
4397  & +pp(jp,8)+pkc11
4398  pp(jp,9)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0
4399  & +pp(jp,9)+pkc12
4400  ELSE
4401  pp(jp,8)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0+pp(jp,8)
4402  pp(jp,9)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0+pp(jp,9)
4403  pp(jp,6)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0
4404  & +pp(jp,6)+pkc11
4405  pp(jp,7)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0
4406  & +pp(jp,7)+pkc12
4407  ENDIF
4408  pp(jp,1)=pp(jp,6)+pp(jp,8)
4409  pp(jp,2)=pp(jp,7)+pp(jp,9)
4410 C ********pt kick for proj
4411  pt(jt,1)=pt11-pkc21
4412  pt(jt,2)=pt12-pkc22
4413  IF((kickdit.EQ.0.AND.rlu(0).LT.0.5)
4414  & .OR.(kickdit.NE.0.AND.rlu(0)
4415  & .LT.0.5/(1.0+(pkc21**2+pkc22**2)/hipr1(22)**2))) THEN
4416  pt(jt,6)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0+pt(jt,6)
4417  pt(jt,7)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0+pt(jt,7)
4418  pt(jt,8)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0
4419  & +pt(jt,8)+pkc21
4420  pt(jt,9)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0
4421  & +pt(jt,9)+pkc22
4422  ELSE
4423  pt(jt,8)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0+pt(jt,8)
4424  pt(jt,9)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0+pt(jt,9)
4425  pt(jt,6)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0
4426  & +pt(jt,6)+pkc21
4427  pt(jt,7)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0
4428  & +pt(jt,7)+pkc22
4429  ENDIF
4430  pt(jt,1)=pt(jt,6)+pt(jt,8)
4431  pt(jt,2)=pt(jt,7)+pt(jt,9)
4432 C ********pt kick for targ
4433 
4434  IF(npj(jp).NE.0) nfp(jp,5)=3
4435  IF(ntj(jt).NE.0) nft(jt,5)=3
4436 C ********jets must be connected to string
4437  IF(epp/(epm+0.0001).LT.etp/(etm+0.0001).AND.
4438  & abs(nfp(jp,1)*nfp(jp,2)).LT.1000000)THEN
4439  DO 620 jsb=1,15
4440  psb=pp(jp,jsb)
4441  pp(jp,jsb)=pt(jt,jsb)
4442  pt(jt,jsb)=psb
4443  nsb=nfp(jp,jsb)
4444  nfp(jp,jsb)=nft(jt,jsb)
4445  nft(jt,jsb)=nsb
4446 620 CONTINUE
4447 C ********when Ycm(JP)<Ycm(JT) after the collision
4448 C exchange the positions of the two
4449  ENDIF
4450 C
4451  RETURN
4452 C**************************************************
4453 C**************************************************
4454 1000 ierror=1
4455  IF(ihpr2(10).EQ.0) RETURN
4456  WRITE(6,*) ' Fatal HIJSFT start error,abandon this event'
4457  WRITE(6,*) ' PROJ E+,E-,W+',epp,epm,wp
4458  WRITE(6,*) ' TARG E+,E-,W-',etp,etm,wm
4459  WRITE(6,*) ' W+*W-, (APN+ATN)^2',sw,snn
4460  RETURN
4461 2000 ierror=0
4462  IF(ihpr2(10).EQ.0) RETURN
4463  WRITE(6,*) ' (2)energy partition fail,'
4464  WRITE(6,*) ' HIJSFT not performed, but continue'
4465  WRITE(6,*) ' MP1,MPN',x1*(1.0-x2)*sw,ampn**2
4466  WRITE(6,*) ' MT2,MTN',x2*(1.0-x1)*sw,amtn**2
4467  RETURN
4468 3000 ierror=0
4469  IF(ihpr2(10).EQ.0) RETURN
4470  WRITE(6,*) ' (3)something is wrong with the pt kick, '
4471  WRITE(6,*) ' HIJSFT not performed, but continue'
4472  WRITE(6,*) ' D1=',d1,' D2=',d2,' SW=',sw
4473  WRITE(6,*) ' HISTORY NFP5=',nfp(jp,5),' NFT5=',nft(jt,5)
4474  WRITE(6,*) ' THIS COLLISON NFP5=',nfp5, ' NFT5=',nft5
4475  WRITE(6,*) ' # OF JET IN PROJ',npj(jp),' IN TARG',ntj(jt)
4476  RETURN
4477 4000 ierror=0
4478  IF(ihpr2(10).EQ.0) RETURN
4479  WRITE(6,*) ' (4)unable to choose process, but not harmful'
4480  WRITE(6,*) ' HIJSFT not performed, but continue'
4481  WRITE(6,*) ' PTP=',sqrt(ptp2),' PTT=',sqrt(ptt2),' SW=',sw
4482  WRITE(6,*) ' AMCUT=',amx,' JP=',jp,' JT=',jt
4483  WRITE(6,*) ' HISTORY NFP5=',nfp(jp,5),' NFT5=',nft(jt,5)
4484  RETURN
4485 5000 ierror=0
4486  IF(ihpr2(10).EQ.0) RETURN
4487  WRITE(6,*) ' energy partition failed(5),for limited try'
4488  WRITE(6,*) ' HIJSFT not performed, but continue'
4489  WRITE(6,*) ' NFP5=',nfp5,' NFT5=',nft5
4490  WRITE(6,*) ' D1',d1,' X1(1-X2)',x1*(1.0-x2)
4491  WRITE(6,*) ' D2',d2,' X2(1-X1)',x2*(1.0-x1)
4492  RETURN
4493 6000 pkc=0.0
4494  miss=miss+1
4495  IF(miss.LT.100) go to 30
4496  ierror=1
4497  IF(ihpr2(10).EQ.0) RETURN
4498  WRITE(6,*) ' ERROR OCCURED, HIJSFT NOT PERFORMED'
4499  WRITE(6,*) ' Abort this event'
4500  WRITE(6,*) 'MTP,PTP2',epp*epm,ptp2,' MTT,PTT2',etp*etm,ptt2
4501  RETURN
4502  END
4503 C
4504 C
4505 C***************************************
4506  SUBROUTINE hijflv(ID)
4507  common/ranseed/nseed
4508  SAVE /ranseed/
4509  id=1
4510  rnid=rlu(0)
4511  IF(rnid.GT.0.43478) THEN
4512  id=2
4513  IF(rnid.GT.0.86956) id=3
4514  ENDIF
4515  RETURN
4516  END
4517 C
4518 C
4519 C
4520  SUBROUTINE hiptdi(PT,PTMAX,IOPT)
4521  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4522  SAVE /hiparnt/
4523  common/ranseed/nseed
4524  SAVE /ranseed/
4525  IF(iopt.EQ.2) THEN
4526  pt=hirnd2(7,0.0,ptmax)
4527  IF(pt.GT.hipr1(8))
4528  & pt=hipr1(2)*sqrt(-alog(exp(-hipr1(8)**2/hipr1(2)**2)
4529  & -rlu(0)*(exp(-hipr1(8)**2/hipr1(2)**2)-
4530  & exp(-ptmax**2/hipr1(2)**2))))
4531 
4532  ELSE
4533  pt=hipr1(2)*sqrt(-alog(1.0-rlu(0)*
4534  & (1.0-exp(-ptmax**2/hipr1(2)**2))))
4535  ENDIF
4536  ptmax0=max(ptmax,0.01)
4537  pt=min(ptmax0-0.01,pt)
4538  RETURN
4539  END
4540 C*************************
4541 C
4542 C
4543 C
4544 C
4545 C ********************************************************
4546 C ************************ WOOD-SAX
4547  SUBROUTINE hijwds(IA,IDH,XHIGH)
4548 C SETS UP HISTOGRAM IDH WITH RADII FOR
4549 C NUCLEUS IA DISTRIBUTED ACCORDING TO THREE PARAM WOOD SAXON
4550  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4551  SAVE /hiparnt/
4552  common/wood/r,d,fnorm,w
4553  SAVE /wood/
4554  dimension iaa(20),rr(20),dd(20),ww(20),rms(20)
4555  EXTERNAL rwdsax,wdsax
4556 C
4557 C PARAMETERS OF SPECIAL NUCLEI FROM ATOMIC DATA AND NUC DATA TABLES
4558 C VOL 14, 5-6 1974
4559  DATA iaa/2,4,12,16,27,32,40,56,63,93,184,197,208,7*0./
4560  DATA rr/0.01,.964,2.355,2.608,2.84,3.458,3.766,3.971,4.214,
4561  1 4.87,6.51,6.38,6.624,7*0./
4562  DATA dd/0.5882,.322,.522,.513,.569,.61,.586,.5935,.586,.573,
4563  1 .535,.535,.549,7*0./
4564  DATA ww/0.0,.517,-0.149,-0.051,0.,-0.208,-0.161,13*0./
4565  DATA rms/2.11,1.71,2.46,2.73,3.05,3.247,3.482,3.737,3.925,4.31,
4566  1 5.42,5.33,5.521,7*0./
4567 
4568  SAVE iaa, rr, dd, ww, rms
4569 C
4570  a=ia
4571 C
4572 C ********SET WOOD-SAX PARAMS FIRST AS IN DATE ET AL
4573  d=0.54
4574 C ********D IS WOOD SAX DIFFUSE PARAM IN FM
4575  r=1.19*a**(1./3.) - 1.61*a**(-1./3.)
4576 C ********R IS RADIUS PARAM
4577  w=0.
4578 C ********W IS The third of three WOOD-SAX PARAM
4579 C
4580 C ********CHECK TABLE FOR SPECIAL CASES
4581  DO 10 i=1,13
4582  IF (ia.EQ.iaa(i)) THEN
4583  r=rr(i)
4584  d=dd(i)
4585  w=ww(i)
4586  rs=rms(i)
4587  END IF
4588 10 CONTINUE
4589 C ********FNORM is the normalize factor
4590  fnorm=1.0
4591  xlow=0.
4592  xhigh=r+ 12.*d
4593  IF (w.LT.-0.01) THEN
4594  IF (xhigh.GT.r/sqrt(abs(w))) xhigh=r/sqrt(abs(w))
4595  END IF
4596  fgaus=gauss1(rwdsax,xlow,xhigh,0.001)
4597  fnorm=1./fgaus
4598 C
4599  IF (idh.EQ.1) THEN
4600  hint1(72)=r
4601  hint1(73)=d
4602  hint1(74)=w
4603  hint1(75)=fnorm/4.0/hipr1(40)
4604  ELSE IF (idh.EQ.2) THEN
4605  hint1(76)=r
4606  hint1(77)=d
4607  hint1(78)=w
4608  hint1(79)=fnorm/4.0/hipr1(40)
4609  ENDIF
4610 C
4611 C NOW SET UP HBOOK FUNCTIONS IDH FOR R**2*RHO(R)
4612 C THESE HISTOGRAMS ARE USED TO GENERATE RANDOM RADII
4613  CALL hifun(idh,xlow,xhigh,rwdsax)
4614  RETURN
4615  END
4616 C
4617 C
4618  FUNCTION wdsax(X)
4619 C ********THREE PARAMETER WOOD SAXON
4620  common/wood/r,d,fnorm,w
4621  SAVE /wood/
4622  wdsax=fnorm*(1.+w*(x/r)**2)/(1+exp((x-r)/d))
4623  IF (w.LT.0.) THEN
4624  IF (x.GE.r/sqrt(abs(w))) wdsax=0.
4625  ENDIF
4626  RETURN
4627  END
4628 C
4629 C
4630  FUNCTION rwdsax(X)
4631  rwdsax=x*x*wdsax(x)
4632  RETURN
4633  END
4634 C
4635 C
4636 C
4637 C
4638  FUNCTION wdsax1(X)
4639 C ********THREE PARAMETER WOOD SAXON
4640 C FOR PROJECTILE
4641 C HINT1(72)=R, HINT1(73)=D, HINT1(74)=W, HINT1(75)=FNORM
4642 C
4643  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4644  SAVE /hiparnt/
4645  wdsax1=hint1(75)*(1.+hint1(74)*(x/hint1(72))**2)/
4646  & (1+exp((x-hint1(72))/hint1(73)))
4647  IF (hint1(74).LT.0.) THEN
4648  IF (x.GE.hint1(72)/sqrt(abs(hint1(74)))) wdsax1=0.
4649  ENDIF
4650  RETURN
4651  END
4652 C
4653 C
4654  FUNCTION wdsax2(X)
4655 C ********THREE PARAMETER WOOD SAXON
4656 C FOR TARGET
4657 C HINT1(76)=R,HINT1(77)=D, HINT1(78)=W, HINT1(79)=FNORM
4658 C
4659  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4660  SAVE /hiparnt/
4661  wdsax2=hint1(79)*(1.+hint1(78)*(x/hint1(76))**2)/
4662  & (1+exp((x-hint1(76))/hint1(77)))
4663  IF (hint1(78).LT.0.) THEN
4664  IF (x.GE.hint1(76)/sqrt(abs(hint1(78)))) wdsax2=0.
4665  ENDIF
4666  RETURN
4667  END
4668 C
4669 C
4670 C THIS FUNCTION IS TO CALCULATE THE NUCLEAR PROFILE FUNCTION
4671 C OF THE COLLIDERING SYSTEM (IN UNITS OF 1/mb)
4672 C
4673  FUNCTION profile(XB)
4674  common/pact/bb,b1,phi,z1
4675  SAVE /pact/
4676  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4677  SAVE /hiparnt/
4678  EXTERNAL flap, flap1, flap2
4679 C
4680  bb=xb
4681  profile=1.0
4682  IF(ihnt2(1).GT.1 .AND. ihnt2(3).GT.1) THEN
4683  profile=float(ihnt2(1))*float(ihnt2(3))*0.1*
4684  & gauss1(flap,0.0,hipr1(34),0.01)
4685  ELSE IF(ihnt2(1).EQ.1 .AND. ihnt2(3).GT.1) THEN
4686  profile=0.2*float(ihnt2(3))*
4687  & gauss1(flap2,0.0,hipr1(35),0.001)
4688  ELSE IF(ihnt2(1).GT.1 .AND. ihnt2(3).EQ.1) THEN
4689  profile=0.2*float(ihnt2(1))*
4690  & gauss1(flap1,0.0,hipr1(34),0.001)
4691  ENDIF
4692  RETURN
4693  END
4694 C
4695 C
4696  FUNCTION flap(X)
4697  common/pact/bb,b1,phi,z1
4698  SAVE /pact/
4699  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4700  SAVE /hiparnt/
4701  EXTERNAL fgp1
4702  b1=x
4703  flap=gauss2(fgp1,0.0,2.0*hipr1(40),0.01)
4704  RETURN
4705  END
4706 C
4707  FUNCTION fgp1(X)
4708  common/pact/bb,b1,phi,z1
4709  SAVE /pact/
4710  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4711  SAVE /hiparnt/
4712  EXTERNAL fgp2
4713  phi=x
4714  fgp1=2.0*gauss3(fgp2,0.0,hipr1(34),0.01)
4715  RETURN
4716  END
4717 C
4718  FUNCTION fgp2(X)
4719  common/pact/bb,b1,phi,z1
4720  SAVE /pact/
4721  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4722  SAVE /hiparnt/
4723  EXTERNAL fgp3
4724  z1=x
4725  fgp2=2.0*gauss4(fgp3,0.0,hipr1(35),0.01)
4726  RETURN
4727  END
4728 C
4729  FUNCTION fgp3(X)
4730  common/pact/bb,b1,phi,z1
4731  SAVE /pact/
4732  r1=sqrt(b1**2+z1**2)
4733  r2=sqrt(bb**2+b1**2-2.0*b1*bb*cos(phi)+x**2)
4734  fgp3=b1*wdsax1(r1)*wdsax2(r2)
4735  RETURN
4736  END
4737 C
4738 C
4739  FUNCTION flap1(X)
4740  common/pact/bb,b1,phi,z1
4741  SAVE /pact/
4742  r=sqrt(bb**2+x**2)
4743  flap1=wdsax1(r)
4744  RETURN
4745  END
4746 C
4747 C
4748  FUNCTION flap2(X)
4749  common/pact/bb,b1,phi,z1
4750  SAVE /pact/
4751  r=sqrt(bb**2+x**2)
4752  flap2=wdsax2(r)
4753  RETURN
4754  END
4755 C
4756 C
4757 C The next three subroutines are for Monte Carlo generation
4758 C according to a given function FHB. One calls first HIFUN
4759 C with assigned channel number I, low and up limits. Then to
4760 C generate the distribution one can call HIRND(I) which gives
4761 C you a random number generated according to the given function.
4762 C
4763  SUBROUTINE hifun(I,XMIN,XMAX,FHB)
4764  common/hijhb/rr(10,201),xx(10,201)
4765  SAVE /hijhb/
4766  EXTERNAL fhb
4767  fnorm=gauss1(fhb,xmin,xmax,0.001)
4768  DO 100 j=1,201
4769  xx(i,j)=xmin+(xmax-xmin)*(j-1)/200.0
4770  xdd=xx(i,j)
4771  rr(i,j)=gauss1(fhb,xmin,xdd,0.001)/fnorm
4772 100 CONTINUE
4773  RETURN
4774  END
4775 C
4776 C
4777 C
4778  FUNCTION hirnd(I)
4779  common/hijhb/rr(10,201),xx(10,201)
4780  SAVE /hijhb/
4781  common/ranseed/nseed
4782  SAVE /ranseed/
4783  rx=rlu(0)
4784  jl=0
4785  ju=202
4786 10 IF(ju-jl.GT.1) THEN
4787  jm=(ju+jl)/2
4788  IF((rr(i,201).GT.rr(i,1)).EQV.(rx.GT.rr(i,jm))) THEN
4789  jl=jm
4790  ELSE
4791  ju=jm
4792  ENDIF
4793  go to 10
4794  ENDIF
4795  j=jl
4796  IF(j.LT.1) j=1
4797  IF(j.GE.201) j=200
4798  hirnd=(xx(i,j)+xx(i,j+1))/2.0
4799  RETURN
4800  END
4801 C
4802 C
4803 C
4804 C
4805 C This generate random number between XMIN and XMAX
4806  FUNCTION hirnd2(I,XMIN,XMAX)
4807  common/hijhb/rr(10,201),xx(10,201)
4808  SAVE /hijhb/
4809  common/ranseed/nseed
4810  SAVE /ranseed/
4811  IF(xmin.LT.xx(i,1)) xmin=xx(i,1)
4812  IF(xmax.GT.xx(i,201)) xmax=xx(i,201)
4813  jmin=1+200*(xmin-xx(i,1))/(xx(i,201)-xx(i,1))
4814  jmax=1+200*(xmax-xx(i,1))/(xx(i,201)-xx(i,1))
4815  rx=rr(i,jmin)+(rr(i,jmax)-rr(i,jmin))*rlu(0)
4816  jl=0
4817  ju=202
4818 10 IF(ju-jl.GT.1) THEN
4819  jm=(ju+jl)/2
4820  IF((rr(i,201).GT.rr(i,1)).EQV.(rx.GT.rr(i,jm))) THEN
4821  jl=jm
4822  ELSE
4823  ju=jm
4824  ENDIF
4825  go to 10
4826  ENDIF
4827  j=jl
4828  IF(j.LT.1) j=1
4829  IF(j.GE.201) j=200
4830  hirnd2=(xx(i,j)+xx(i,j+1))/2.0
4831  RETURN
4832  END
4833 C
4834 C
4835 C
4836 C
4837  SUBROUTINE hijcrs
4838 C THIS IS TO CALCULATE THE CROSS SECTIONS OF JET PRODUCTION AND
4839 C THE TOTAL INELASTIC CROSS SECTIONS.
4840  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4841  SAVE /hiparnt/
4842  common/njet/n,ip_crs
4843  SAVE /njet/
4844  EXTERNAL fhin,ftot,fnjet,ftotjet,ftotrig
4845  IF(hint1(1).GE.10.0) CALL crsjet
4846 C ********calculate jet cross section(in mb)
4847 C
4848  aphx1=hipr1(6)*(ihnt2(1)**0.3333333-1.0)
4849  aphx2=hipr1(6)*(ihnt2(3)**0.3333333-1.0)
4850  hint1(11)=hint1(14)-aphx1*hint1(15)
4851  & -aphx2*hint1(16)+aphx1*aphx2*hint1(17)
4852  hint1(10)=gauss1(ftotjet,0.0,20.0,0.01)
4853  hint1(12)=gauss1(fhin,0.0,20.0,0.01)
4854  hint1(13)=gauss1(ftot,0.0,20.0,0.01)
4855  hint1(60)=hint1(61)-aphx1*hint1(62)
4856  & -aphx2*hint1(63)+aphx1*aphx2*hint1(64)
4857  hint1(59)=gauss1(ftotrig,0.0,20.0,0.01)
4858  IF(hint1(59).EQ.0.0) hint1(59)=hint1(60)
4859  IF(hint1(1).GE.10.0) Then
4860  DO 20 i=0,20
4861  n=i
4862  hint1(80+i)=gauss1(fnjet,0.0,20.0,0.01)/hint1(12)
4863  20 CONTINUE
4864  ENDIF
4865  hint1(10)=hint1(10)*hipr1(31)
4866  hint1(12)=hint1(12)*hipr1(31)
4867  hint1(13)=hint1(13)*hipr1(31)
4868  hint1(59)=hint1(59)*hipr1(31)
4869 C ********Total and Inel cross section are calculated
4870 C by Gaussian integration.
4871  IF(ihpr2(13).NE.0) THEN
4872  hipr1(33)=1.36*(1.0+36.0/hint1(1)**2)
4873  & *alog(0.6+0.1*hint1(1)**2)
4874  hipr1(33)=hipr1(33)/hint1(12)
4875  ENDIF
4876 C ********Parametrized cross section for single
4877 C diffractive reaction(Goulianos)
4878  RETURN
4879  END
4880 C
4881 C
4882 C
4883 C
4884  FUNCTION ftot(X)
4885  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4886  SAVE /hiparnt/
4887  omg=omg0(x)*(hipr1(30)+hint1(11))/hipr1(31)/2.0
4888  ftot=2.0*(1.0-exp(-omg))
4889  RETURN
4890  END
4891 C
4892 C
4893 C
4894  FUNCTION fhin(X)
4895  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4896  SAVE /hiparnt/
4897  omg=omg0(x)*(hipr1(30)+hint1(11))/hipr1(31)/2.0
4898  fhin=1.0-exp(-2.0*omg)
4899  RETURN
4900  END
4901 C
4902 C
4903 C
4904  FUNCTION ftotjet(X)
4905  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4906  SAVE /hiparnt/
4907  omg=omg0(x)*hint1(11)/hipr1(31)/2.0
4908  ftotjet=1.0-exp(-2.0*omg)
4909  RETURN
4910  END
4911 C
4912 C
4913 C
4914  FUNCTION ftotrig(X)
4915  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4916  SAVE /hiparnt/
4917  omg=omg0(x)*hint1(60)/hipr1(31)/2.0
4918  ftotrig=1.0-exp(-2.0*omg)
4919  RETURN
4920  END
4921 C
4922 C
4923 C
4924 C
4925  FUNCTION fnjet(X)
4926  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4927  SAVE /hiparnt/
4928  common/njet/n,ip_crs
4929  SAVE /njet/
4930  omg1=omg0(x)*hint1(11)/hipr1(31)
4931  c0=exp(n*alog(omg1)-sgmin(n+1))
4932  IF(n.EQ.0) c0=1.0-exp(-2.0*omg0(x)*hipr1(30)/hipr1(31)/2.0)
4933  fnjet=c0*exp(-omg1)
4934  RETURN
4935  END
4936 C
4937 C
4938 C
4939 C
4940 C
4941  FUNCTION sgmin(N)
4942  ga=0.
4943  IF(n.LE.2) go to 20
4944  DO 10 i=1,n-1
4945  z=i
4946  ga=ga+alog(z)
4947 10 CONTINUE
4948 20 sgmin=ga
4949  RETURN
4950  END
4951 C
4952 C
4953 C
4954  FUNCTION omg0(X)
4955  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
4956  SAVE /hiparnt/
4957  COMMON /besel/x4
4958  SAVE /besel/
4959  EXTERNAL bk
4960  x4=hipr1(32)*sqrt(x)
4961  omg0=hipr1(32)**2*gauss2(bk,x4,x4+20.0,0.01)/96.0
4962 c write(*,*)'OMG0 start',x,x4,OMG0
4963  RETURN
4964  END
4965 C
4966 C
4967 C
4968  FUNCTION romg(X)
4969 C ********This gives the eikonal function from a table
4970 C calculated in the first call
4971  dimension fr(0:1000)
4972 c --------------Khaled
4973 c we need to save FR because it is calculated only once
4974  common/eikonal/fr
4975 c --------------->end Khaled
4976 
4977  DATA i0/0/
4978  save i0 !Uzhi
4979  IF(i0.NE.0) go to 100
4980  DO 50 i=1,1001
4981  xr=(i-1)*0.01
4982  fr(i-1)=omg0(xr)
4983 
4984 50 CONTINUE
4985 100 i0=1
4986 
4987  IF(x.GE.10.0) THEN
4988  romg=0.0
4989  RETURN
4990  ENDIF
4991  ix=int(x*100)
4992 
4993  romg=(fr(ix)*((ix+1)*0.01-x)+fr(ix+1)*(x-ix*0.01))/0.01
4994 
4995  RETURN
4996  END
4997 C
4998 C
4999 C
5000 
5001  FUNCTION bk(X)
5002  COMMON /besel/x4
5003  SAVE /besel/
5004  bk=exp(-x)*(x**2-x4**2)**2.50/15.0
5005  RETURN
5006  END
5007 C
5008 C
5009 C THIS PROGRAM IS TO CALCULATE THE JET CROSS SECTION
5010 C THE INTEGRATION IS DONE BY USING VEGAS
5011 C
5012  SUBROUTINE crsjet
5013  IMPLICIT REAL*8(a-h,o-z)
5014  REAL hipr1(100),hint1(100)
5015  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5016  SAVE /hiparnt/
5017  common/njet/n,ip_crs
5018  SAVE /njet/
5019  common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn
5020  SAVE /bveg1/
5021  common/bveg2/xi(50,10),si,si2,swgt,schi,ndo,it
5022  SAVE /bveg2/
5023  common/bveg3/f,ti,tsi
5024  SAVE /bveg3/
5025  common/seedvax/num1
5026  SAVE /seedvax/
5027  EXTERNAL fjet,fjetrig
5028 C
5029 c************************
5030 c NCALL give the number of inner-iteration, ITMX
5031 C gives the limit of out-iteration. Nprn is an option
5032 C ( 1: print the integration process. 0: do not print)
5033 C
5034  ndim=3
5035  ip_crs=0
5036  CALL vegas(fjet,avgi,sd,chi2a)
5037  hint1(14)=avgi/2.5682
5038  IF(ihpr2(6).EQ.1 .AND. ihnt2(1).GT.1) THEN
5039  ip_crs=1
5040  CALL vegas(fjet,avgi,sd,chi2a)
5041  hint1(15)=avgi/2.5682
5042  ENDIF
5043  IF(ihpr2(6).EQ.1 .AND. ihnt2(3).GT.1) THEN
5044  ip_crs=2
5045  CALL vegas(fjet,avgi,sd,chi2a)
5046  hint1(16)=avgi/2.5682
5047  ENDIF
5048  IF(ihpr2(6).EQ.1.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1) THEN
5049  ip_crs=3
5050  CALL vegas(fjet,avgi,sd,chi2a)
5051  hint1(17)=avgi/2.5682
5052  ENDIF
5053 C ********Total inclusive jet cross section(Pt>P0)
5054 C
5055  IF(ihpr2(3).NE.0) THEN
5056  ip_crs=0
5057  CALL vegas(fjetrig,avgi,sd,chi2a)
5058  hint1(61)=avgi/2.5682
5059  IF(ihpr2(6).EQ.1 .AND. ihnt2(1).GT.1) THEN
5060  ip_crs=1
5061  CALL vegas(fjetrig,avgi,sd,chi2a)
5062  hint1(62)=avgi/2.5682
5063  ENDIF
5064  IF(ihpr2(6).EQ.1 .AND. ihnt2(3).GT.1) THEN
5065  ip_crs=2
5066  CALL vegas(fjetrig,avgi,sd,chi2a)
5067  hint1(63)=avgi/2.5682
5068  ENDIF
5069  IF(ihpr2(6).EQ.1.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1) THEN
5070  ip_crs=3
5071  CALL vegas(fjetrig,avgi,sd,chi2a)
5072  hint1(64)=avgi/2.5682
5073  ENDIF
5074  ENDIF
5075 C ********cross section of trigger jet
5076 C
5077  RETURN
5078  END
5079 C
5080 C
5081 C
5082  FUNCTION fjet(X,WGT)
5083  IMPLICIT REAL*8(a-h,o-z)
5084  REAL hipr1(100),hint1(100)
5085  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5086  SAVE /hiparnt/
5087  dimension x(10)
5088  pt2=(hint1(1)**2/4.0-hipr1(8)**2)*x(1)+hipr1(8)**2
5089  xt=2.0*dsqrt(pt2)/hint1(1)
5090  ymx1=dlog(1.0/xt+dsqrt(1.0/xt**2-1.0))
5091  y1=2.0*ymx1*x(2)-ymx1
5092  ymx2=dlog(2.0/xt-dexp(y1))
5093  ymn2=dlog(2.0/xt-dexp(-y1))
5094  y2=(ymx2+ymn2)*x(3)-ymn2
5095  fjet=2.0*ymx1*(ymx2+ymn2)*(hint1(1)**2/4.0-hipr1(8)**2)
5096  & *g(y1,y2,pt2)/2.0
5097  RETURN
5098  END
5099 C
5100 C
5101 C
5102  FUNCTION fjetrig(X,WGT)
5103  IMPLICIT REAL*8(a-h,o-z)
5104  REAL hipr1(100),hint1(100),ptmax,ptmin
5105  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5106  SAVE /hiparnt/
5107  dimension x(10)
5108  ptmin=abs(hipr1(10))-0.25
5109  ptmin=max(ptmin,hipr1(8))
5110  am2=0.d0
5111  IF(ihpr2(3).EQ.3) THEN
5112  am2=hipr1(7)**2
5113  ptmin=max(0.0,hipr1(10))
5114  ENDIF
5115  ptmax=abs(hipr1(10))+0.25
5116  IF(hipr1(10).LE.0.0) ptmax=hint1(1)/2.0-am2
5117  IF(ptmax.LE.ptmin) ptmax=ptmin+0.25
5118  pt2=(ptmax**2-ptmin**2)*x(1)+ptmin**2
5119  amt2=pt2+am2
5120  xt=2.0*dsqrt(amt2)/hint1(1)
5121  ymx1=dlog(1.0/xt+dsqrt(1.0/xt**2-1.0))
5122  y1=2.0*ymx1*x(2)-ymx1
5123  ymx2=dlog(2.0/xt-dexp(y1))
5124  ymn2=dlog(2.0/xt-dexp(-y1))
5125  y2=(ymx2+ymn2)*x(3)-ymn2
5126  IF(ihpr2(3).EQ.3) THEN
5127  gtrig=2.0*ghvq(y1,y2,amt2)
5128  ELSE IF(ihpr2(3).EQ.2) THEN
5129  gtrig=2.0*gphoton(y1,y2,pt2)
5130  ELSE
5131  gtrig=g(y1,y2,pt2)
5132  ENDIF
5133  fjetrig=2.0*ymx1*(ymx2+ymn2)*(ptmax**2-ptmin**2)
5134  & *gtrig/2.0
5135  RETURN
5136  END
5137 C
5138 C
5139 C
5140  FUNCTION ghvq(Y1,Y2,AMT2)
5141  IMPLICIT REAL*8 (a-h,o-z)
5142  REAL hipr1(100),hint1(100)
5143  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5144  SAVE /hiparnt/
5145  dimension f(2,7)
5146  xt=2.0*dsqrt(amt2)/hint1(1)
5147  x1=0.50*xt*(dexp(y1)+dexp(y2))
5148  x2=0.50*xt*(dexp(-y1)+dexp(-y2))
5149  ss=x1*x2*hint1(1)**2
5150  af=4.0
5151  IF(ihpr2(18).NE.0) af=5.0
5152  dlam=hipr1(15)
5153  aph=12.0*3.1415926/(33.0-2.0*af)/dlog(amt2/dlam**2)
5154 C
5155  CALL parton(f,x1,x2,amt2)
5156 C
5157  gqq=4.0*(cosh(y1-y2)+hipr1(7)**2/amt2)/(1.d0+cosh(y1-y2))/9.0
5158  & *(f(1,1)*f(2,2)+f(1,2)*f(2,1)+f(1,3)*f(2,4)
5159  & +f(1,4)*f(2,3)+f(1,5)*f(2,6)+f(1,6)*f(2,5))
5160  ggg=(8.d0*cosh(y1-y2)-1.d0)*(cosh(y1-y2)+2.0*hipr1(7)**2/amt2
5161  & -2.0*hipr1(7)**4/amt2**2)/(1.0+cosh(y1-y2))/24.d0
5162  & *f(1,7)*f(2,7)
5163 C
5164  ghvq=(gqq+ggg)*hipr1(23)*3.14159*aph**2/ss**2
5165  RETURN
5166  END
5167 C
5168 C
5169 C
5170  FUNCTION gphoton(Y1,Y2,PT2)
5171  IMPLICIT REAL*8 (a-h,o-z)
5172  REAL hipr1(100),hint1(100)
5173  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5174  SAVE /hiparnt/
5175  dimension f(2,7)
5176  xt=2.0*dsqrt(pt2)/hint1(1)
5177  x1=0.50*xt*(dexp(y1)+dexp(y2))
5178  x2=0.50*xt*(dexp(-y1)+dexp(-y2))
5179  z=dsqrt(1.d0-xt**2/x1/x2)
5180  ss=x1*x2*hint1(1)**2
5181  t=-(1.0-z)/2.0
5182  u=-(1.0+z)/2.0
5183  af=3.0
5184  dlam=hipr1(15)
5185  aph=12.0*3.1415926/(33.0-2.0*af)/dlog(pt2/dlam**2)
5186  aphem=1.0/137.0
5187 C
5188  CALL parton(f,x1,x2,pt2)
5189 C
5190  g11=-(u**2+1.0)/u/3.0*f(1,7)*(4.0*f(2,1)+4.0*f(2,2)
5191  & +f(2,3)+f(2,4)+f(2,5)+f(2,6))/9.0
5192  g12=-(t**2+1.0)/t/3.0*f(2,7)*(4.0*f(1,1)+4.0*f(1,2)
5193  & +f(1,3)+f(1,4)+f(1,5)+f(1,6))/9.0
5194  g2=8.0*(u**2+t**2)/u/t/9.0*(4.0*f(1,1)*f(2,2)
5195  & +4.0*f(1,2)*f(2,1)+f(1,3)*f(2,4)+f(1,4)*f(2,3)
5196  & +f(1,5)*f(2,6)+f(1,6)*f(2,5))/9.0
5197 C
5198  gphoton=(g11+g12+g2)*hipr1(23)*3.14159*aph*aphem/ss**2
5199  RETURN
5200  END
5201 C
5202 C
5203 C
5204 C
5205  FUNCTION g(Y1,Y2,PT2)
5206  IMPLICIT REAL*8 (a-h,o-z)
5207  REAL hipr1(100),hint1(100)
5208  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5209  SAVE /hiparnt/
5210  dimension f(2,7)
5211  xt=2.0*dsqrt(pt2)/hint1(1)
5212  x1=0.50*xt*(dexp(y1)+dexp(y2))
5213  x2=0.50*xt*(dexp(-y1)+dexp(-y2))
5214  z=dsqrt(1.d0-xt**2/x1/x2)
5215  ss=x1*x2*hint1(1)**2
5216  t=-(1.0-z)/2.0
5217  u=-(1.0+z)/2.0
5218  af=3.0
5219  dlam=hipr1(15)
5220  aph=12.0*3.1415926/(33.0-2.0*af)/dlog(pt2/dlam**2)
5221 C
5222  CALL parton(f,x1,x2,pt2)
5223 C
5224  g11=( (f(1,1)+f(1,2))*(f(2,3)+f(2,4)+f(2,5)+f(2,6))
5225  & +(f(1,3)+f(1,4))*(f(2,5)+f(2,6)) )*subcrs1(t,u)
5226 C
5227  g12=( (f(2,1)+f(2,2))*(f(1,3)+f(1,4)+f(1,5)+f(1,6))
5228  & +(f(2,3)+f(2,4))*(f(1,5)+f(1,6)) )*subcrs1(u,t)
5229 C
5230  g13=(f(1,1)*f(2,1)+f(1,2)*f(2,2)+f(1,3)*f(2,3)+f(1,4)*f(2,4)
5231  & +f(1,5)*f(2,5)+f(1,6)*f(2,6))*(subcrs1(u,t)
5232  & +subcrs1(t,u)-8.d0/t/u/27.d0)
5233 C
5234  g2=(af-1)*(f(1,1)*f(2,2)+f(2,1)*f(1,2)+f(1,3)*f(2,4)
5235  & +f(2,3)*f(1,4)+f(1,5)*f(2,6)+f(2,5)*f(1,6))*subcrs2(t,u)
5236 C
5237  g31=(f(1,1)*f(2,2)+f(1,3)*f(2,4)+f(1,5)*f(2,6))*subcrs3(t,u)
5238  g32=(f(2,1)*f(1,2)+f(2,3)*f(1,4)+f(2,5)*f(1,6))*subcrs3(u,t)
5239 C
5240  g4=(f(1,1)*f(2,2)+f(2,1)*f(1,2)+f(1,3)*f(2,4)+f(2,3)*f(1,4)+
5241  1 f(1,5)*f(2,6)+f(2,5)*f(1,6))*subcrs4(t,u)
5242 C
5243  g5=af*f(1,7)*f(2,7)*subcrs5(t,u)
5244 C
5245  g61=f(1,7)*(f(2,1)+f(2,2)+f(2,3)+f(2,4)+f(2,5)
5246  & +f(2,6))*subcrs6(t,u)
5247  g62=f(2,7)*(f(1,1)+f(1,2)+f(1,3)+f(1,4)+f(1,5)
5248  & +f(1,6))*subcrs6(u,t)
5249 C
5250  g7=f(1,7)*f(2,7)*subcrs7(t,u)
5251 C
5252  g=(g11+g12+g13+g2+g31+g32+g4+g5+g61+g62+g7)*hipr1(17)*
5253  1 3.14159d0*aph**2/ss**2
5254  RETURN
5255  END
5256 C
5257 C
5258 C
5259  FUNCTION subcrs1(T,U)
5260  IMPLICIT REAL*8(a-h,o-z)
5261  subcrs1=4.d0/9.d0*(1.d0+u**2)/t**2
5262  RETURN
5263  END
5264 C
5265 C
5266  FUNCTION subcrs2(T,U)
5267  IMPLICIT REAL*8(a-h,o-z)
5268  subcrs2=4.d0/9.d0*(t**2+u**2)
5269  RETURN
5270  END
5271 C
5272 C
5273  FUNCTION subcrs3(T,U)
5274  IMPLICIT REAL*8(a-h,o-z)
5275  subcrs3=4.d0/9.d0*(t**2+u**2+(1.d0+u**2)/t**2
5276  1 -2.d0*u**2/3.d0/t)
5277  RETURN
5278  END
5279 C
5280 C
5281  FUNCTION subcrs4(T,U)
5282  IMPLICIT REAL*8(a-h,o-z)
5283  subcrs4=8.d0/3.d0*(t**2+u**2)*(4.d0/9.d0/t/u-1.d0)
5284  RETURN
5285  END
5286 C
5287 C
5288 C
5289  FUNCTION subcrs5(T,U)
5290  IMPLICIT REAL*8(a-h,o-z)
5291  subcrs5=3.d0/8.d0*(t**2+u**2)*(4.d0/9.d0/t/u-1.d0)
5292  RETURN
5293  END
5294 C
5295 C
5296  FUNCTION subcrs6(T,U)
5297  IMPLICIT REAL*8(a-h,o-z)
5298  subcrs6=(1.d0+u**2)*(1.d0/t**2-4.d0/u/9.d0)
5299  RETURN
5300  END
5301 C
5302 C
5303  FUNCTION subcrs7(T,U)
5304  IMPLICIT REAL*8(a-h,o-z)
5305  subcrs7=9.d0/2.d0*(3.d0-t*u-u/t**2-t/u**2)
5306  RETURN
5307  END
5308 C
5309 C
5310 C
5311  SUBROUTINE parton(F,X1,X2,QQ)
5312  IMPLICIT REAL*8(a-h,o-z)
5313  REAL hipr1(100),hint1(100)
5314  common/hiparnt/hipr1,ihpr2(50),hint1,ihnt2(50)
5315  SAVE /hiparnt/
5316  common/njet/n,ip_crs
5317  SAVE /njet/
5318  dimension f(2,7)
5319  dlam=hipr1(15)
5320  q0=hipr1(16)
5321  s=dlog(dlog(qq/dlam**2)/dlog(q0**2/dlam**2))
5322  IF(ihpr2(7).EQ.2) go to 200
5323 C*******************************************************
5324  at1=0.419+0.004*s-0.007*s**2
5325  at2=3.460+0.724*s-0.066*s**2
5326  gmud=4.40-4.86*s+1.33*s**2
5327  at3=0.763-0.237*s+0.026*s**2
5328  at4=4.00+0.627*s-0.019*s**2
5329  gmd=-0.421*s+0.033*s**2
5330 C*******************************************************
5331  cas=1.265-1.132*s+0.293*s**2
5332  as=-0.372*s-0.029*s**2
5333  bs=8.05+1.59*s-0.153*s**2
5334  aphs=6.31*s-0.273*s**2
5335  btas=-10.5*s-3.17*s**2
5336  gms=14.7*s+9.80*s**2
5337 C********************************************************
5338 C CAC=0.135*S-0.075*S**2
5339 C AC=-0.036-0.222*S-0.058*S**2
5340 C BC=6.35+3.26*S-0.909*S**2
5341 C APHC=-3.03*S+1.50*S**2
5342 C BTAC=17.4*S-11.3*S**2
5343 C GMC=-17.9*S+15.6*S**2
5344 C***********************************************************
5345  cag=1.56-1.71*s+0.638*s**2
5346  ag=-0.949*s+0.325*s**2
5347  bg=6.0+1.44*s-1.05*s**2
5348  aphg=9.0-7.19*s+0.255*s**2
5349  btag=-16.5*s+10.9*s**2
5350  gmg=15.3*s-10.1*s**2
5351  go to 300
5352 C********************************************************
5353 200 at1=0.374+0.014*s
5354  at2=3.33+0.753*s-0.076*s**2
5355  gmud=6.03-6.22*s+1.56*s**2
5356  at3=0.761-0.232*s+0.023*s**2
5357  at4=3.83+0.627*s-0.019*s**2
5358  gmd=-0.418*s+0.036*s**2
5359 C************************************
5360  cas=1.67-1.92*s+0.582*s**2
5361  as=-0.273*s-0.164*s**2
5362  bs=9.15+0.530*s-0.763*s**2
5363  aphs=15.7*s-2.83*s**2
5364  btas=-101.0*s+44.7*s**2
5365  gms=223.0*s-117.0*s**2
5366 C*********************************
5367 C CAC=0.067*S-0.031*S**2
5368 C AC=-0.120-0.233*S-0.023*S**2
5369 C BC=3.51+3.66*S-0.453*S**2
5370 C APHC=-0.474*S+0.358*S**2
5371 C BTAC=9.50*S-5.43*S**2
5372 C GMC=-16.6*S+15.5*S**2
5373 C**********************************
5374  cag=0.879-0.971*s+0.434*s**2
5375  ag=-1.16*s+0.476*s**2
5376  bg=4.0+1.23*s-0.254*s**2
5377  aphg=9.0-5.64*s-0.817*s**2
5378  btag=-7.54*s+5.50*s**2
5379  gmg=-0.596*s+1.26*s**2
5380 C*********************************
5381 300 b12=dexp(gmre(at1)+gmre(at2+1.d0)-gmre(at1+at2+1.d0))
5382  b34=dexp(gmre(at3)+gmre(at4+1.d0)-gmre(at3+at4+1.d0))
5383  cnud=3.d0/b12/(1.d0+gmud*at1/(at1+at2+1.d0))
5384  cnd=1.d0/b34/(1.d0+gmd*at3/(at3+at4+1.d0))
5385 C********************************************************
5386 C FUD=X*(U+D)
5387 C FS=X*2(UBAR+DBAR+SBAR) AND UBAR=DBAR=SBAR
5388 C*******************************************************
5389  fud1=cnud*x1**at1*(1.d0-x1)**at2*(1.d0+gmud*x1)
5390  fs1=cas*x1**as*(1.d0-x1)**bs*(1.d0+aphs*x1
5391  & +btas*x1**2+gms*x1**3)
5392  f(1,3)=cnd*x1**at3*(1.d0-x1)**at4*(1.d0+gmd*x1)+fs1/6.d0
5393  f(1,1)=fud1-f(1,3)+fs1/3.d0
5394  f(1,2)=fs1/6.d0
5395  f(1,4)=fs1/6.d0
5396  f(1,5)=fs1/6.d0
5397  f(1,6)=fs1/6.d0
5398  f(1,7)=cag*x1**ag*(1.d0-x1)**bg*(1.d0+aphg*x1
5399  & +btag*x1**2+gmg*x1**3)
5400 C
5401  fud2=cnud*x2**at1*(1.d0-x2)**at2*(1.d0+gmud*x2)
5402  fs2=cas*x2**as*(1.d0-x2)**bs*(1.d0+aphs*x2
5403  & +btas*x2**2+gms*x2**3)
5404  f(2,3)=cnd*x2**at3*(1.d0-x2)**at4*(1.d0+gmd*x2)+fs2/6.d0
5405  f(2,1)=fud2-f(2,3)+fs2/3.d0
5406  f(2,2)=fs2/6.d0
5407  f(2,4)=fs2/6.d0
5408  f(2,5)=fs2/6.d0
5409  f(2,6)=fs2/6.d0
5410  f(2,7)=cag*x2**ag*(1.d0-x2)**bg*(1.d0+aphg*x2
5411  & +btag*x2**2+gmg*x2**3)
5412 C***********Nuclear effect on the structure function****************
5413 C
5414  IF(ihpr2(6).EQ.1 .AND. ihnt2(1).GT.1) THEN
5415  aax=1.193*alog(float(ihnt2(1)))**0.16666666
5416  rrx=aax*(x1**3-1.2*x1**2+0.21*x1)+1.0
5417  & +1.079*(float(ihnt2(1))**0.33333333-1.0)
5418  & /dlog(ihnt2(1)+1.0d0)*dsqrt(x1)*dexp(-x1**2/0.01)
5419  IF(ip_crs.EQ.1 .OR.ip_crs.EQ.3) rrx=dexp(-x1**2/0.01)
5420  DO 400 i=1,7
5421  f(1,i)=rrx*f(1,i)
5422  400 CONTINUE
5423  ENDIF
5424  IF(ihpr2(6).EQ.1 .AND. ihnt2(3).GT.1) THEN
5425  aax=1.193*alog(float(ihnt2(3)))**0.16666666
5426  rrx=aax*(x2**3-1.2*x2**2+0.21*x2)+1.0
5427  & +1.079*(float(ihnt2(3))**0.33333-1.0)
5428  & /dlog(ihnt2(3)+1.0d0)*dsqrt(x2)*dexp(-x2**2/0.01)
5429  IF(ip_crs.EQ.2 .OR. ip_crs.EQ.3) rrx=dexp(-x2**2/0.01)
5430  DO 500 i=1,7
5431  f(2,i)=rrx*f(2,i)
5432  500 CONTINUE
5433  ENDIF
5434 c
5435  RETURN
5436  END
5437 C
5438 C
5439 C
5440  FUNCTION gmre(X)
5441  IMPLICIT REAL*8(a-h,o-z)
5442  z=x
5443  IF(x.GT.3.0d0) go to 10
5444  z=x+3.d0
5445 10 gmre=0.5d0*dlog(2.d0*3.14159265d0/z)+z*dlog(z)-z+dlog(1.d0
5446  1 +1.d0/12.d0/z+1.d0/288.d0/z**2-139.d0/51840.d0/z**3
5447  1 -571.d0/2488320.d0/z**4)
5448  IF(z.EQ.x) go to 20
5449  gmre=gmre-dlog(z-1.d0)-dlog(z-2.d0)-dlog(z-3.d0)
5450 20 CONTINUE
5451  RETURN
5452  END
5453 c
5454 C
5455 C
5456  FUNCTION gmin(N)
5457  IMPLICIT REAL*8(a-h,o-z)
5458  ga=0.
5459  IF(n.LE.2) go to 20
5460  DO 10 i=1,n-1
5461  z=i
5462  ga=ga+dlog(z)
5463 10 CONTINUE
5464 20 gmin=ga
5465  RETURN
5466  END
5467 C
5468 C
5469 C***************************************************************
5470 
5471  BLOCK DATA hidata
5472  REAL*8 xl(10),xu(10),acc
5473  common/bveg1/xl,xu,acc,ndim,ncall,itmx,nprn
5474  SAVE /bveg1/
5475  common/seedvax/num1
5476  SAVE /seedvax/
5477  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
5478  SAVE /hiparnt/
5479  common/ranseed/nseed
5480  SAVE /ranseed/
5481  common/himain1/ natt,eatt,jatt,nt,np,n0,n01,n10,n11
5482  SAVE /himain1/
5483  common/himain2/katt(130000,4),patt(130000,4)
5484  SAVE /himain2/
5485  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
5486  SAVE /histrng/
5487  common/hijcrdn/yp(3,300),yt(3,300)
5488  SAVE /hijcrdn/
5489  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
5490  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
5491  & pjpm(300,500),ntj(300),kftj(300,500),
5492  & pjtx(300,500),pjty(300,500),pjtz(300,500),
5493  & pjte(300,500),pjtm(300,500)
5494  SAVE /hijjet1/
5495  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
5496  & k2sg(900,100),pxsg(900,100),pysg(900,100),
5497  & pzsg(900,100),pesg(900,100),pmsg(900,100)
5498  SAVE /hijjet2/
5499  common/hijdat/hidat0(10,10),hidat(10)
5500  SAVE /hijdat/
5501  common/hipyint/mint4,mint5,atco(200,20),atxs(0:200)
5502  SAVE /hipyint/
5503  DATA num1/30123984/,xl/10*0.d0/,xu/10*1.d0/
5504  DATA ncall/1000/itmx/100/acc/0.01/nprn/0/
5505 C...give all the switchs and parameters the default values
5506 
5507  DATA nseed/74769375/
5508  DATA hipr1/
5509  & 1.5, 0.35, 0.5, 0.9, 2.0, 0.1, 1.5, 2.0, -1.0, -2.25,
5510  & 2.0, 0.5, 1.0, 2.0, 0.2, 2.0, 2.5, 0.3, 0.1, 1.4,
5511  & 1.6, 1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 57.0,
5512  & 28.5, 3.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
5513  & 3.141592654,
5514  & 0.0, 0.4, 0.1, 1.5, 0.1, 0.25, 0.0, 0.5, 0.0, 0.0,
5515  & 50*0.0/
5516 
5517  DATA ihpr2/
5518  & 1, 3, 0, 1, 1, 1, 1, 10, 0, 0,
5519  & 1, 1, 1, 1, 0, 0, 1, 0, 0, 1,
5520  & 30*0/
5521 
5522  DATA hint1/100*0/
5523  DATA ihnt2/50*0/
5524 
5525 C...initialize all the data common blocks
5526  DATA natt/0/eatt/0.0/jatt/0/nt/0/np/0/n0/0/n01/0/n10/0/n11/0/
5527  DATA katt/520000*0/patt/520000*0.0/
5528 
5529  DATA nfp/4500*0/pp/4500*0.0/nft/4500*0/pt/4500*0.0/
5530 
5531  DATA yp/900*0.0/yt/900*0.0/
5532 
5533  DATA npj/300*0/kfpj/150000*0/pjpx/150000*0.0/pjpy/150000*0.0/
5534  & pjpz/150000*0.0/pjpe/150000*0.0/pjpm/150000*0.0/
5535  DATA ntj/300*0/kftj/150000*0/pjtx/150000*0.0/pjty/150000*0.0/
5536  & pjtz/150000*0.0/pjte/150000*0.0/pjtm/150000*0.0/
5537 
5538  DATA nsg/0/njsg/900*0/iasg/2700*0/k1sg/90000*0/k2sg/90000*0/
5539  & pxsg/90000*0.0/pysg/90000*0.0/pzsg/90000*0.0/pesg/90000*0.0/
5540  & pmsg/90000*0.0/
5541  DATA mint4/0/mint5/0/atco/4000*0.0/atxs/201*0.0/
5542  DATA (hidat0(1,i),i=1,10)/0.0,0.0,0.0,0.0,0.0,0.0,2.25,
5543  & 2.5,4.0,4.1/
5544  DATA (hidat0(2,i),i=1,10)/2.0,3.0,5.0,6.0,7.0,8.0,8.0,10.0,
5545  & 10.0,10.0/
5546  DATA (hidat0(3,i),i=1,10)/1.0,0.8,0.8,0.7,0.45,0.215,
5547  & 0.21,0.19,0.19,0.19/
5548  DATA (hidat0(4,i),i=1,10)/0.35,0.35,0.3,0.3,0.3,0.3,
5549  & 0.5,0.6,0.6,0.6/
5550  DATA (hidat0(5,i),i=1,10)/23.8,24.0,26.0,26.2,27.0,28.5,28.5,
5551  & 28.5,28.5,28.5/
5552  DATA ((hidat0(j,i),i=1,10),j=6,9)/40*0.0/
5553  DATA (hidat0(10,i),i=1,10)/5.0,20.0,53.0,62.0,100.0,200.0,
5554  & 546.0,900.0,1800.0,4000.0/
5555  DATA hidat/10*0.0/
5556  END
5557 C*******************************************************************
5558 C
5559 C
5560 C
5561 C
5562 C*******************************************************************
5563 C SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
5564 C - BY G.P. LEPAGE SEPT 1976/(REV)APR 1978
5565 C*******************************************************************
5566 C
5567  SUBROUTINE vegas(FXN,AVGI,SD,CHI2A)
5568  IMPLICIT REAL*8(a-h,o-z)
5569  common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn
5570  SAVE /bveg1/
5571  common/bveg2/xi(50,10),si,si2,swgt,schi,ndo,it
5572  SAVE /bveg2/
5573  common/bveg3/f,ti,tsi
5574  SAVE /bveg3/
5575  EXTERNAL fxn
5576  dimension d(50,10),di(50,10),xin(50),r(50),dx(10),dt(10),x(10)
5577  1 ,kg(10),ia(10)
5578  REAL*4 qran(10)
5579  DATA ndmx/50/,alph/1.5d0/,one/1.d0/,mds/-1/
5580 
5581  save ndmx, alph, one, mds !uzhi
5582 C
5583  ndo=1
5584  DO 1 j=1,ndim
5585 1 xi(1,j)=one
5586 C
5587  entry vegas1(fxn,avgi,sd,chi2a)
5588 C - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
5589  it=0
5590  si=0.
5591  si2=si
5592  swgt=si
5593  schi=si
5594 C
5595  entry vegas2(fxn,avgi,sd,chi2a)
5596 C - NO INITIALIZATION
5597  nd=ndmx
5598  ng=1
5599  IF(mds.EQ.0) go to 2
5600  ng=(ncall/2.)**(1./ndim)
5601  mds=1
5602  IF((2*ng-ndmx).LT.0) go to 2
5603  mds=-1
5604  npg=ng/ndmx+1
5605  nd=ng/npg
5606  ng=npg*nd
5607 2 k=ng**ndim
5608  npg=ncall/k
5609  IF(npg.LT.2) npg=2
5610  calls=npg*k
5611  dxg=one/ng
5612  dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-one)
5613  xnd=nd
5614  ndm=nd-1
5615  dxg=dxg*xnd
5616  xjac=one/calls
5617  DO 3 j=1,ndim
5618 c***this is the line 50
5619  dx(j)=xu(j)-xl(j)
5620 3 xjac=xjac*dx(j)
5621 C
5622 C REBIN PRESERVING BIN DENSITY
5623 C
5624  IF(nd.EQ.ndo) go to 8
5625  rc=ndo/xnd
5626  DO 7 j=1,ndim
5627  k=0
5628  xn=0.
5629  dr=xn
5630  i=k
5631 4 k=k+1
5632  dr=dr+one
5633  xo=xn
5634  xn=xi(k,j)
5635 5 IF(rc.GT.dr) go to 4
5636  i=i+1
5637  dr=dr-rc
5638  xin(i)=xn-(xn-xo)*dr
5639  IF(i.LT.ndm) go to 5
5640  DO 6 i=1,ndm
5641 6 xi(i,j)=xin(i)
5642 7 xi(nd,j)=one
5643  ndo=nd
5644 C
5645 8 CONTINUE
5646  IF(nprn.NE.0) WRITE(16,200) ndim,calls,it,itmx,acc,mds,nd
5647  1 ,(xl(j),xu(j),j=1,ndim)
5648 C
5649  entry vegas3(fxn,avgi,sd,chi2a)
5650 C - MAIN INTEGRATION LOOP
5651 9 it=it+1
5652  ti=0.
5653  tsi=ti
5654  DO 10 j=1,ndim
5655  kg(j)=1
5656  DO 10 i=1,nd
5657  d(i,j)=ti
5658 10 di(i,j)=ti
5659 C
5660 11 fb=0.
5661  f2b=fb
5662  k=0
5663 12 k=k+1
5664  CALL aran9(qran,ndim)
5665  wgt=xjac
5666  DO 15 j=1,ndim
5667  xn=(kg(j)-qran(j))*dxg+one
5668 c*****this is the line 100
5669  ia(j)=xn
5670  IF(ia(j).GT.1) go to 13
5671  xo=xi(ia(j),j)
5672  rc=(xn-ia(j))*xo
5673  go to 14
5674 13 xo=xi(ia(j),j)-xi(ia(j)-1,j)
5675  rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
5676 14 x(j)=xl(j)+rc*dx(j)
5677  wgt=wgt*xo*xnd
5678 15 CONTINUE
5679 C
5680  f=wgt
5681  f=f*fxn(x,wgt)
5682  f2=f*f
5683  fb=fb+f
5684  f2b=f2b+f2
5685  DO 16 j=1,ndim
5686  di(ia(j),j)=di(ia(j),j)+f
5687 16 IF(mds.GE.0) d(ia(j),j)=d(ia(j),j)+f2
5688  IF(k.LT.npg) go to 12
5689 C
5690  f2b=dsqrt(f2b*npg)
5691  f2b=(f2b-fb)*(f2b+fb)
5692  ti=ti+fb
5693  tsi=tsi+f2b
5694  IF(mds.GE.0) go to 18
5695  DO 17 j=1,ndim
5696 17 d(ia(j),j)=d(ia(j),j)+f2b
5697 18 k=ndim
5698 19 kg(k)=mod(kg(k),ng)+1
5699  IF(kg(k).NE.1) go to 11
5700  k=k-1
5701  IF(k.GT.0) go to 19
5702 C
5703 C FINAL RESULTS FOR THIS ITERATION
5704 C
5705  tsi=tsi*dv2g
5706  ti2=ti*ti
5707  wgt=ti2/(tsi+1.0d-37)
5708  si=si+ti*wgt
5709  si2=si2+ti2
5710  swgt=swgt+wgt
5711  swgt=swgt+1.0d-37
5712  si2=si2+1.0d-37
5713  schi=schi+ti2*wgt
5714  avgi=si/(swgt)
5715  sd=swgt*it/(si2)
5716  chi2a=sd*(schi/swgt-avgi*avgi)/(it-.999)
5717  sd=dsqrt(one/sd)
5718 C****this is the line 150
5719  IF(nprn.EQ.0) go to 21
5720  tsi=dsqrt(tsi)
5721  WRITE(16,201) it,ti,tsi,avgi,sd,chi2a
5722  IF(nprn.GE.0) go to 21
5723  DO 20 j=1,ndim
5724 20 WRITE(16,202) j,(xi(i,j),di(i,j),d(i,j),i=1,nd)
5725 C
5726 C REFINE GRID
5727 C
5728 21 DO 23 j=1,ndim
5729  xo=d(1,j)
5730  xn=d(2,j)
5731  d(1,j)=(xo+xn)/2.
5732  dt(j)=d(1,j)
5733  DO 22 i=2,ndm
5734  d(i,j)=xo+xn
5735  xo=xn
5736  xn=d(i+1,j)
5737  d(i,j)=(d(i,j)+xn)/3.
5738 22 dt(j)=dt(j)+d(i,j)
5739  d(nd,j)=(xn+xo)/2.
5740 23 dt(j)=dt(j)+d(nd,j)
5741 C
5742  DO 28 j=1,ndim
5743  rc=0.
5744  DO 24 i=1,nd
5745  r(i)=0.
5746  IF (dt(j).GE.1.0d18) THEN
5747  WRITE(6,*) '************** A SINGULARITY >1.0D18'
5748 C WRITE(5,1111)
5749 C1111 FORMAT(1X,'**************IMPORTANT NOTICE***************')
5750 C WRITE(5,1112)
5751 C1112 FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
5752 C WRITE(5,1113)
5753 C1113 FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
5754 C WRITE(5,1114)
5755 C1114 FORMAT(1X,'**************END NOTICE*************')
5756  END IF
5757  IF(d(i,j).LE.1.0d-18) go to 24
5758  xo=dt(j)/d(i,j)
5759  r(i)=((xo-one)/xo/dlog(xo))**alph
5760 24 rc=rc+r(i)
5761  rc=rc/xnd
5762  k=0
5763  xn=0.
5764  dr=xn
5765  i=k
5766 25 k=k+1
5767  dr=dr+r(k)
5768  xo=xn
5769 c****this is the line 200
5770  xn=xi(k,j)
5771 26 IF(rc.GT.dr) go to 25
5772  i=i+1
5773  dr=dr-rc
5774  xin(i)=xn-(xn-xo)*dr/(r(k)+1.0d-30)
5775  IF(i.LT.ndm) go to 26
5776  DO 27 i=1,ndm
5777 27 xi(i,j)=xin(i)
5778 28 xi(nd,j)=one
5779 C
5780  IF(it.LT.itmx.AND.acc*dabs(avgi).LT.sd) go to 9
5781 200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',i3,' NCALL=',f8.0
5782  1 /28x,' IT=',i5,' ITMX=',i5/28x,' ACC=',g9.3
5783  2 /28x,' MDS=',i3,' ND=',i4/28x,' (XL,XU)=',
5784  3 (t40,'( ',g12.6,' , ',g12.6,' )'))
5785 201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',i3,
5786  1 ': INTEGRAL =',g14.8/21x,'STD DEV =',g10.4 /
5787  2 ' ACCUMULATED RESULTS: INTEGRAL =',g14.8 /
5788  3 24x,'STD DEV =',g10.4 / 24x,'CHI**2 PER IT''N =',g10.4)
5789 202 FORMAT('0DATA FOR AXIS',i2 / ' ',6x,'X',7x,' DELT I ',
5790  1 2x,' CONV''CE ',11x,'X',7x,' DELT I ',2x,' CONV''CE '
5791  2 ,11x,'X',7x,' DELT I ',2x,' CONV''CE ' /
5792  2 (' ',3g12.4,5x,3g12.4,5x,3g12.4))
5793  RETURN
5794  END
5795 C
5796 C
5797  SUBROUTINE aran9(QRAN,NDIM)
5798  dimension qran(10)
5799  common/seedvax/num1
5800  DO 1 i=1,ndim
5801  1 qran(i)=rlu(0)
5802  RETURN
5803  END
5804 
5805 C
5806 C
5807 C*********GAUSSIAN ONE-DIMENSIONAL INTEGRATION PROGRAM*************
5808 C
5809  FUNCTION gauss1(F,A,B,EPS)
5810  EXTERNAL f
5811  dimension w(12),x(12)
5812  DATA const/1.0e-12/
5813  DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5814  & .0622535,0.0951585,.1246290,.1495960,.1691565,
5815  & .1826034,.1894506/
5816  DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5817  & .9445750,0.8656312,.7554044,.6178762,.4580168,
5818  & .2816036,.0950125/
5819  delta=const*abs(a-b)
5820  gauss1=0.0
5821  aa=a
5822 5 y=b-aa
5823  IF(abs(y).LE.delta) RETURN
5824 2 bb=aa+y
5825  c1=0.5*(aa+bb)
5826  c2=c1-aa
5827  s8=0.0
5828  s16=0.0
5829  DO 1 i=1,4
5830  u=x(i)*c2
5831 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
5832  DO 3 i=5,12
5833  u=x(i)*c2
5834 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
5835  s8=s8*c2
5836  s16=s16*c2
5837  IF(abs(s16-s8).GT.eps*(1.+abs(s16))) goto 4
5838  gauss1=gauss1+s16
5839  aa=bb
5840  goto 5
5841 4 y=0.5*y
5842  IF(abs(y).GT.delta) goto 2
5843  WRITE(6,7)
5844  gauss1=0.0
5845  RETURN
5846 7 FORMAT(1x,'GAUSS1....TOO HIGH ACURACY REQUIRED')
5847  END
5848 C
5849 C
5850 C
5851  FUNCTION gauss2(F,A,B,EPS)
5852  EXTERNAL f
5853  dimension w(12),x(12)
5854  DATA const/1.0e-12/
5855  DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5856  & .0622535,0.0951585,.1246290,.1495960,.1691565,
5857  & .1826034,.1894506/
5858  DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5859  & .9445750,0.8656312,.7554044,.6178762,.4580168,
5860  & .2816036,.0950125/
5861  delta=const*abs(a-b)
5862  gauss2=0.0
5863  aa=a
5864 5 y=b-aa
5865  IF(abs(y).LE.delta) RETURN
5866 2 bb=aa+y
5867  c1=0.5*(aa+bb)
5868  c2=c1-aa
5869  s8=0.0
5870  s16=0.0
5871  DO 1 i=1,4
5872  u=x(i)*c2
5873 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
5874  DO 3 i=5,12
5875  u=x(i)*c2
5876 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
5877  s8=s8*c2
5878  s16=s16*c2
5879  IF(abs(s16-s8).GT.eps*(1.+abs(s16))) goto 4
5880  gauss2=gauss2+s16
5881  aa=bb
5882  goto 5
5883 4 y=0.5*y
5884  IF(abs(y).GT.delta) goto 2
5885  WRITE(6,7)
5886  gauss2=0.0
5887  RETURN
5888 7 FORMAT(1x,'GAUSS2....TOO HIGH ACURACY REQUIRED')
5889  END
5890 C
5891 C
5892 C
5893  FUNCTION gauss3(F,A,B,EPS)
5894  EXTERNAL f
5895  dimension w(12),x(12)
5896  DATA const/1.0e-12/
5897  DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5898  & .0622535,0.0951585,.1246290,.1495960,.1691565,
5899  & .1826034,.1894506/
5900  DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5901  & .9445750,0.8656312,.7554044,.6178762,.4580168,
5902  & .2816036,.0950125/
5903  delta=const*abs(a-b)
5904  gauss3=0.0
5905  aa=a
5906 5 y=b-aa
5907  IF(abs(y).LE.delta) RETURN
5908 2 bb=aa+y
5909  c1=0.5*(aa+bb)
5910  c2=c1-aa
5911  s8=0.0
5912  s16=0.0
5913  DO 1 i=1,4
5914  u=x(i)*c2
5915 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
5916  DO 3 i=5,12
5917  u=x(i)*c2
5918 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
5919  s8=s8*c2
5920  s16=s16*c2
5921  IF(abs(s16-s8).GT.eps*(1.+abs(s16))) goto 4
5922  gauss3=gauss3+s16
5923  aa=bb
5924  goto 5
5925 4 y=0.5*y
5926  IF(abs(y).GT.delta) goto 2
5927  WRITE(6,7)
5928  gauss3=0.0
5929  RETURN
5930 7 FORMAT(1x,'GAUSS3....TOO HIGH ACURACY REQUIRED')
5931  END
5932 C
5933 C
5934 C
5935 C
5936  FUNCTION gauss4(F,A,B,EPS)
5937  EXTERNAL f
5938  dimension w(12),x(12)
5939  DATA const/1.0e-12/
5940  DATA w/0.1012285,.2223810,.3137067,.3623838,.0271525,
5941  & .0622535,0.0951585,.1246290,.1495960,.1691565,
5942  & .1826034,.1894506/
5943  DATA x/0.9602899,.7966665,.5255324,.1834346,.9894009,
5944  & .9445750,0.8656312,.7554044,.6178762,.4580168,
5945  & .2816036,.0950125/
5946  delta=const*abs(a-b)
5947  gauss4=0.0
5948  aa=a
5949 5 y=b-aa
5950  IF(abs(y).LE.delta) RETURN
5951 2 bb=aa+y
5952  c1=0.5*(aa+bb)
5953  c2=c1-aa
5954  s8=0.0
5955  s16=0.0
5956  DO 1 i=1,4
5957  u=x(i)*c2
5958 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
5959  DO 3 i=5,12
5960  u=x(i)*c2
5961 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
5962  s8=s8*c2
5963  s16=s16*c2
5964  IF(abs(s16-s8).GT.eps*(1.+abs(s16))) goto 4
5965  gauss4=gauss4+s16
5966  aa=bb
5967  goto 5
5968 4 y=0.5*y
5969  IF(abs(y).GT.delta) goto 2
5970  WRITE(6,7)
5971  gauss4=0.0
5972  RETURN
5973 7 FORMAT(1x,'GAUSS4....TOO HIGH ACURACY REQUIRED')
5974  END
5975 C
5976 C
5977 C
5978 C
5979 C
5980  SUBROUTINE title
5981  WRITE(6,200)
5982 200 FORMAT(//10x,
5983  & '**************************************************'/10x,
5984  & '* | \\ _______ / ------/ *'/10x,
5985  & '* ----- ------ |_____| /_/ / *'/10x,
5986  & '* ||\\ / |_____| / / \\ *'/10x,
5987  & '* /| \\ /_/ /_______ /_ / \\_ *'/10x,
5988  & '* / | / / / / / | ------- *'/10x,
5989  & '* | / /\\ / / | / | *'/10x,
5990  & '* | / / \\ / / \\_| / ------- *'/10x,
5991  & '* *'/10x,
5992  & '**************************************************'/10x,
5993  & ' HIJING '/10x,
5994  & ' Heavy Ion Jet INteraction Generator '/10x,
5995  & ' by '/10x,
5996  & ' X. N. Wang and M. Gyulassy '/10x,
5997  & ' Lawrence Berkeley Laboratory '//)
5998  RETURN
5999  END