Fix for Npart counting. The bug was causing FPE (Constantin)
[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / hijing.F
1 * $Id$
2 C     Version 1.36
3 C     Nothing important has been changed here. A few 'garbage' has been
4 C     cleaned up here, like common block HIJJET3 for the sea quark strings
5 C     which were originally created to implement the DPM scheme which
6 C     later was abadoned in the final version. The lines which operate
7 C     on these data are also deleted in the program.
8 C
9 C
10 C     Version 1.35
11 C     There are some changes in the program: subroutine HARDJET is now
12 C     consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
13 C     for the triggered hard processes. Now that is done  altogether
14 C     with other normal hard processes in modified JETINI. In the new
15 C     version one calls JETINI every time one calls HIJHRD. In the new
16 C     version the effect of the isospin of the nucleon on hard processes,
17 C     especially direct photons is correctly considered.
18 C     For A+A collisions, one has to initilize pythia
19 C     separately for each type of collisions, pp, pn,np and nn,
20 C     or hp and hn for hA collisions. In JETINI we use the following
21 C     catalogue for different types of collisions:
22 C     h+h: h+h (I_TYPE=1)
23 C     h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
24 C     A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
25 C     A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
26 C*****************************************************************
27 c
28 C
29 C     Version 1.34
30 C     Last modification on January 5, 1998. Two misstakes are corrected in
31 C     function G. A Misstake in the subroutine Parton is also corrected.
32 C     (These are pointed out by Ysushi Nara).
33 C
34 C
35 C       Last modifcation on April 10, 1996. To conduct final
36 C       state radiation, PYTHIA reorganize the two scattered
37 C       partons and their final momenta will be a little
38 C       different. The summed total momenta of the partons
39 C       from the final state radiation are stored in HINT1(26-29)
40 C       and HINT1(36-39) which are little different from 
41 C       HINT1(21-24) and HINT1(41-44).
42 C
43 C       Version 1.33
44 C
45 C       Last modfication  on September 11, 1995. When HIJING and
46 C       PYTHIA are initialized, the shadowing is evaluated at
47 C       b=0 which is the maximum. This will cause overestimate
48 C       of shadowing for peripheral interactions. To correct this
49 C       problem, shadowing is set to zero when initializing. Then
50 C       use these maximum  cross section without shadowing as a
51 C       normalization of the Monte Carlo. This however increase
52 C       the computing time. IHNT2(16) is used to indicate whether
53 C       the sturcture function is called for (IHNT2(16)=1) initialization
54 C       or for (IHNT2(16)=0)normal collisions simulation
55 C
56 C       Last modification on Aagust 28, 1994. Two bugs associate
57 C       with the impact parameter dependence of the shadowing is
58 C       corrected.
59 C
60 C
61 c       Last modification on October 14, 1994. One bug is corrected
62 c       in the direct photon production option in subroutine
63 C       HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
64 C       Another bug associated with keeping the decay history
65 C       in the particle information is also corrected.(this problem
66 C       was reported by Matt Bloomer)
67 C
68 C
69 C       Last modification on July 15, 1994. The option to trig on
70 C       heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1) 
71 C       is added. To do this, set IHPR2(3)=3. For inclusive production,
72 C       one should reset HIPR1(10)=0.0. One can also trig larger pt
73 C       QQbar production by giving HIPR1(10) a nonvanishing value.
74 C       The mass of the heavy quark in the calculation of the cross
75 C       section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
76 C       default is the charm mass D=1.5). We also include a separate
77 C       K-factor for heavy quark and direct photon production by
78 C       HIPR1(23)(D=2.0).
79 C
80 C       Last modification on May 24, 1994.  The option to
81 C       retain the information of all particles including those
82 C       who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is 
83 C       added to contain the line number of the parent particle 
84 C       of the current line which is produced via a decay. 
85 C       KATT(I,4) is the status number of the particle: 11=particle
86 C       which has decayed; 1=finally produced particle.
87 C
88 C
89 C       Last modification on May 24, 1994( in HIJSFT when valence quark
90 C       is quenched, the following error is corrected. 1.2*IHNT2(1) --> 
91 C       1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
92 C
93 C
94 C       Last modification on March 16, 1994 (heavy flavor production
95 C       processes MSUB(81)=1 MSUB(82)=1 have been switched on,
96 C       charm production is the default, B-quark option is
97 C       IHPR2(18), when it is switched on, charm quark is 
98 C       automatically off)
99 C
100 C
101 C       Last modification on March 23, 1994 (an error is corrected
102 C       in the impact parameter dependence of the jet cross section)
103 C
104 C       Last modification Oct. 1993 to comply with non-vax
105 C       machines' compiler 
106 C
107 C*********************************************
108 C       LAST MODIFICATION April 5, 1991
109 CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B 
110 C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
111 C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
112 C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18), 
113 C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
114 C
115 C       gluon jets can form a single string system
116 C
117 C       initial state radiation is included
118 C       
119 C       all QCD subprocesses are included
120 c
121 c       direct particles production is included(currently only direct
122 C               photon)
123 c
124 C       Effect of high P_T trigger bias on multiple jets distribution
125 c
126 C******************************************************************
127 C                               HIJING.10                         *
128 C                 Heavy Ion Jet INteraction Generator             *
129 C                                  by                             *
130 C                  X. N. Wang      and   M. Gyulassy              *
131 C                     Lawrence Berkeley Laboratory                *
132 C                                                                 *
133 C******************************************************************
134 C
135 C******************************************************************
136 C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of  *
137 C proj, NFP(K,4) original ID of proj.  NFP(K,5)=colli status(0=no,*
138 C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
139 C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it   *
140 C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
141 C (0=has not been,1=is going to be, -1=has already been scattered *
142 C NFP(k,11) total number of interactions this proj has suffered   *
143 C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant  *
144 C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and     *
145 C diquark,PP(K,10)=PT of the hard scattering between the valence  *
146 C quarks; PP(K,14,15)=the mass of quark,diquark.                  * 
147 C******************************************************************
148 C
149 C****************************************************************
150 C
151 C       SUBROUTINE HIJING
152 C
153 C****************************************************************
154         SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
155         CHARACTER FRAME*8
156         DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
157      &                  IPCOL(90000),ITCOL(90000),SCIP2(300,300)
158 #define BLANKET_SAVE
159 #include "hiparnt.inc"
160 C
161 #include "hijcrdn.inc"
162 #include "hijglbr.inc"
163 #include "himain1.inc"
164 #include "himain2.inc"
165 #include "histrng.inc"
166 #include "hijjet1.inc"
167 #include "hijjet2.inc"
168 #include "hijjet4.inc"
169 C
170 #include "lujets_hijing.inc"
171 #include "ludat1_hijing.inc"
172         SAVE
173
174         BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
175         BMIN=MIN(BMIN0,BMAX)
176         IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
177                 BMIN=0.0
178                 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
179         ENDIF
180 C                       ********HIPR1(31) is in mb =0.1fm**2
181 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS 
182 C       BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
183 C
184         YP(1,1)=0.0
185         YP(2,1)=0.0
186         YP(3,1)=0.0
187         IF(IHNT2(1).LE.1) GO TO 14
188         DO 10 KP=1,IHNT2(1)
189 5       R=HIRND(1)
190 c
191         if(IHNT2(1).EQ.2) then
192            rnd1=max(RLU_HIJING(NSEED),1.0e-20)
193            rnd2=max(RLU_HIJING(NSEED),1.0e-20)
194            rnd3=max(RLU_HIJING(NSEED),1.0e-20)
195            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
196      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
197         endif
198 c
199         X=RLU_HIJING(0)
200         CX=2.0*X-1.0
201         SX=SQRT(1.0-CX*CX)
202 C               ********choose theta from uniform cos(theta) distr
203         PHI=RLU_HIJING(0)*2.0*HIPR1(40)
204 C               ********choose phi form uniform phi distr 0 to 2*pi
205         YP(1,KP)=R*SX*COS(PHI)
206         YP(2,KP)=R*SX*SIN(PHI)
207         YP(3,KP)=R*CX
208         IF(HIPR1(29).EQ.0.0) GO TO 10
209         DO 8  KP2=1,KP-1
210                 DNBP1=(YP(1,KP)-YP(1,KP2))**2
211                 DNBP2=(YP(2,KP)-YP(2,KP2))**2
212                 DNBP3=(YP(3,KP)-YP(3,KP2))**2
213                 DNBP=DNBP1+DNBP2+DNBP3
214                 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
215 C                       ********two neighbors cannot be closer than 
216 C                               HIPR1(29)
217 8       CONTINUE
218 10      CONTINUE
219 c*******************************
220         if(IHNT2(1).EQ.2) then
221            YP(1,2)=-YP(1,1)
222            YP(2,2)=-YP(2,1)
223            YP(3,2)=-YP(3,1)
224         endif
225 c********************************
226         DO 12 I=1,IHNT2(1)-1
227         DO 12 J=I+1,IHNT2(1)
228         IF(YP(3,I).GT.YP(3,J)) GO TO 12
229         Y1=YP(1,I)
230         Y2=YP(2,I)
231         Y3=YP(3,I)
232         YP(1,I)=YP(1,J)
233         YP(2,I)=YP(2,J)
234         YP(3,I)=YP(3,J)
235         YP(1,J)=Y1
236         YP(2,J)=Y2
237         YP(3,J)=Y3
238 12      CONTINUE
239 C
240 C******************************
241 14      YT(1,1)=0.0
242         YT(2,1)=0.0
243         YT(3,1)=0.0
244         IF(IHNT2(3).LE.1) GO TO 24
245         DO 20 KT=1,IHNT2(3)
246 15      R=HIRND(2)
247 c
248          if(IHNT2(3).EQ.2) then
249             rnd1=max(RLU_HIJING(NSEED),1.0e-20)
250             rnd2=max(RLU_HIJING(NSEED),1.0e-20)
251             rnd3=max(RLU_HIJING(NSEED),1.0e-20)
252             R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
253      &           +4.38*0.85*log(rnd3)/(4.38+0.85))
254          endif
255 c
256         X=RLU_HIJING(0)
257         CX=2.0*X-1.0
258         SX=SQRT(1.0-CX*CX)
259 C               ********choose theta from uniform cos(theta) distr
260         PHI=RLU_HIJING(0)*2.0*HIPR1(40)
261 C               ********chose phi form uniform phi distr 0 to 2*pi
262         YT(1,KT)=R*SX*COS(PHI)
263         YT(2,KT)=R*SX*SIN(PHI)
264         YT(3,KT)=R*CX
265         IF(HIPR1(29).EQ.0.0) GO TO 20
266         DO 18  KT2=1,KT-1
267                 DNBT1=(YT(1,KT)-YT(1,KT2))**2
268                 DNBT2=(YT(2,KT)-YT(2,KT2))**2
269                 DNBT3=(YT(3,KT)-YT(3,KT2))**2
270                 DNBT=DNBT1+DNBT2+DNBT3
271                 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
272 C                       ********two neighbors cannot be closer than 
273 C                               HIPR1(29)
274 18      CONTINUE
275 20      CONTINUE
276 c**********************************
277          if(IHNT2(3).EQ.2) then
278             YT(1,2)=-YT(1,1)
279             YT(2,2)=-YT(2,1)
280             YT(3,2)=-YT(3,1)
281          endif
282 c*********************************
283         DO 22 I=1,IHNT2(3)-1
284         DO 22 J=I+1,IHNT2(3)
285         IF(YT(3,I).LT.YT(3,J)) GO TO 22
286         Y1=YT(1,I)
287         Y2=YT(2,I)
288         Y3=YT(3,I)
289         YT(1,I)=YT(1,J)
290         YT(2,I)=YT(2,J)
291         YT(3,I)=YT(3,J)
292         YT(1,J)=Y1
293         YT(2,J)=Y2
294         YT(3,J)=Y3
295 22      CONTINUE
296
297 C********************
298 24      MISS=-1
299
300 50      MISS=MISS+1
301         IF(MISS.GT.50) THEN
302            WRITE(6,*) 'infinite loop happened in  HIJING'
303            STOP
304         ENDIF
305
306         NATT=0
307         JATT=0
308         EATT=0.0
309         CALL HIJINI
310         NLOP=0
311 C                       ********Initialize for a new event
312 60      NT=0
313         NP=0
314         N0=0
315         N01=0
316         N10=0
317         N11=0
318         NELT=0
319         NINT=0
320         NELP=0
321         NINP=0
322         NSG=0
323         NCOLT=0
324         NPSPECP=0
325         NNSPECP=0
326         NPSPECT=0
327         NNSPECT=0 
328
329 C****   BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS 
330 C       RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET 
331 C       BY THE ANGLE PHI  FOR EACH COLLISION.******************
332 C
333         BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
334         PHI=2.0*HIPR1(40)*RLU_HIJING(0)
335         BBX=BB*COS(PHI)
336         BBY=BB*SIN(PHI)
337         HINT1(19)=BB
338         HINT1(20)=PHI
339 C
340         DO 70 JP=1,IHNT2(1)
341         DO 70 JT=1,IHNT2(3)
342            SCIP(JP,JT)=-1.0
343            B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
344            R2=B2*HIPR1(40)/HIPR1(31)/0.1
345            SCIP2(JP,JT)=R2
346 C               ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
347            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
348      &          /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
349            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
350      &          /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
351            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
352      &           *SQRT(1.0-RRB1)
353            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
354      &           *SQRT(1.0-RRB2)
355            HINT1(18)=HINT1(14)-APHX1*HINT1(15)
356      &                  -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
357            IF(IHPR2(14).EQ.0.OR.
358      &          (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
359               GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
360               RANTOT=RLU_HIJING(0)
361               IF(RANTOT.GT.GS) GO TO 70
362               GO TO 65
363            ENDIF
364            GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
365      &             /HIPR1(31)/2.0*ROMG(0.0)))
366            R2=R2/GSTOT_0
367            GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
368            GSTOT=2.0*(1.0-SQRT(1.0-GS))
369            RANTOT=RLU_HIJING(0)*GSTOT_0
370            IF(RANTOT.GT.GSTOT) GO TO 70
371            IF(RANTOT.GT.GS) THEN
372               CALL HIJCSC(JP,JT)
373               GO TO 70
374 C                       ********perform elastic collisions
375            ENDIF
376  65        SCIP(JP,JT)=R2
377            RNIP(JP,JT)=RANTOT
378            SJIP(JP,JT)=HINT1(18)
379            NCOLT=NCOLT+1
380            if (R2.GT.2.D0) THEN
381               write (8,*) R2
382            ENDIF
383            IPCOL(NCOLT)=JP
384            ITCOL(NCOLT)=JT
385 70      CONTINUE
386
387 c *** cl glauber ***
388         npart=0
389         xmeana=0D0
390         ymeana=0D0
391         xmeanb=0D0
392         ymeanb=0D0
393         xmeanp=0D0
394         ymeanp=0D0
395         xm2=0D0
396         ym2=0D0
397         xym=0D0
398
399         IF(NCOLT>0) THEN 
400            DO 1110 JP=1,IHNT2(1)
401               YP(1,JP)=YP(1,JP)+BBX
402               YP(2,JP)=YP(2,JP)+BBY
403               xmeana=xmeana+YP(1,JP)
404               ymeana=ymeana+YP(2,JP)
405               DO 1120 JT=1,IHNT2(3)
406                  IF(SCIP2(JP,JT).LT.2.0D0) THEN
407                     npart=npart+1
408                     xmeanp=xmeanp+YP(1,JP)
409                     ymeanp=ymeanp+YP(2,JP)
410                     xm2=xm2+YP(1,JP)*YP(1,JP)
411                     ym2=ym2+YP(2,JP)*YP(2,JP)
412                     xym=xym+YP(1,JP)*YP(2,JP)
413                     goto 1110
414                  end if
415  1120         continue
416  1110      continue
417
418            DO 1130 JT=1,IHNT2(3)
419               xmeanb=xmeanb+YT(1,JT)
420               ymeanb=ymeanb+YT(2,JT)
421               DO 1140 JP=1,IHNT2(1)
422                  IF(SCIP2(JP,JT).LT.2.0D0) THEN
423                     npart=npart+1
424                     xmeanp=xmeanp+YT(1,JT)
425                     ymeanp=ymeanp+YT(2,JT)
426                     xm2=xm2+YT(1,JT)*YT(1,JT)
427                     ym2=ym2+YT(2,JT)*YT(2,JT)
428                     xym=xym+YT(1,JT)*YT(2,JT)
429                     goto 1130
430                  end if
431  1140         continue
432  1130      continue
433
434            IF (npart.GT.0) THEN
435               xmeana=xmeana/IHNT2(1)
436               ymeana=ymeana/IHNT2(1)
437               xmeanb=xmeanb/IHNT2(3)
438               ymeanb=ymeanb/IHNT2(3)
439               xmeanp=xmeanp/npart
440               ymeanp=ymeanp/npart
441               xm2=xm2/npart
442               ym2=ym2/npart
443               xym=xym/npart
444
445               sx2=xm2-xmeanp*xmeanp
446               sy2=ym2-ymeanp*ymeanp
447               sxy=xym-xmeanp*ymeanp
448            
449               delx=xmeanb-xmeana
450               dely=ymeanb-ymeana
451               dtmp=delx**2+dely**2
452               bbtrue=sqrt(dtmp)
453               dnumt=(sy2-sx2)*(delx**2-dely**2)-4D0*sxy*delx*dely
454               ddent=(sy2+sx2)*bbtrue**2
455               eccrp=dnumt/ddent
456               dtmp=(sy2-sx2)*(sy2-sx2)+4D0*sxy*sxy
457               eccpart=sqrt(dtmp)/(sx2+sy2)
458               eccmc=(sy2-sx2)/(sy2+sx2)
459               write(*,*),'HOUT: ',bb,' ',bbtrue,' ',ncolt,' ',npart,
460      1             ' ',eccrp,' ',eccpart
461            end if
462         end if
463
464 C               ********total number interactions proj and targ has
465 C                               suffered
466
467         IF(NCOLT.EQ.0) THEN
468            NLOP=NLOP+1
469            IF(NLOP.LE.20.OR.
470      &           (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
471
472
473            RETURN
474         ENDIF
475 C               ********At large impact parameter, there maybe no
476 C                       interaction at all. For NN collision
477 C                       repeat the event until interaction happens
478 C
479         IF(IHPR2(3).NE.0) THEN
480            NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
481            NHARD=MIN(NHARD,NCOLT)
482            JPHARD=IPCOL(NHARD)
483            JTHARD=ITCOL(NHARD)
484         ENDIF
485 C
486         IF(IHPR2(9).EQ.1) THEN
487                 NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
488                 NMINI=MIN(NMINI,NCOLT)
489                 JPMINI=IPCOL(NMINI)
490                 JTMINI=ITCOL(NMINI)
491         ENDIF
492 C               ********Specifying the location of the hard and
493 C                       minijet if they are enforced by user
494 C
495         DO 200 JP=1,IHNT2(1)
496         DO 200 JT=1,IHNT2(3)
497         IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
498                 NFP(JP,11)=NFP(JP,11)+1
499                 NFT(JT,11)=NFT(JT,11)+1
500         IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
501                 NP=NP+1
502                 N01=N01+1
503         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
504                 NT=NT+1
505                 N10=N10+1
506         ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
507                 NP=NP+1
508                 NT=NT+1
509                 N0=N0+1
510         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
511                 N11=N11+1
512         ENDIF
513         JOUT=0
514         NFP(JP,10)=0
515         NFT(JT,10)=0
516 C*****************************************************************
517         IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
518 C               ********When IHPR2(8)=0 no jets are produced
519         IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
520 C               ********jets can not be produced for (JP,JT)
521 C                       because not enough energy avaible for 
522 C                               JP or JT 
523         R2=SCIP(JP,JT)
524         HINT1(18)=SJIP(JP,JT)
525         TT=ROMG(R2)*HINT1(18)/HIPR1(31)
526         TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
527         NJET=0
528         IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
529            CALL JETINI(JP,JT,1)
530            CALL HIJHRD(JP,JT,0,JFLG,0)
531            HINT1(26)=HINT1(47)
532            HINT1(27)=HINT1(48)
533            HINT1(28)=HINT1(49)
534            HINT1(29)=HINT1(50)
535            HINT1(36)=HINT1(67)
536            HINT1(37)=HINT1(68)
537            HINT1(38)=HINT1(69)
538            HINT1(39)=HINT1(70)
539 C
540            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
541            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
542            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
543      &                          JFLG.GE.3) IASG(NSG,3)=1
544            IHNT2(9)=IHNT2(14)
545            IHNT2(10)=IHNT2(15)
546            DO 105 I05=1,5
547               HINT1(20+I05)=HINT1(40+I05)
548               HINT1(30+I05)=HINT1(50+I05)
549  105       CONTINUE
550            JOUT=1
551            IF(IHPR2(8).EQ.0) GO TO 160
552            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
553      &          /REAL(IHNT2(1))**0.6666667,1.0)
554            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
555      &          /REAL(IHNT2(3))**0.6666667,1.0)
556            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
557      &           *SQRT(1.0-RRB1)
558            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
559      &           *SQRT(1.0-RRB2)
560            HINT1(65)=HINT1(61)-APHX1*HINT1(62)
561      &                  -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
562            TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
563            NJET=-1
564 C               ********subtract the trigger jet from total number
565 C                       of jet production  to be done since it has
566 C                               already been produced here
567            XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
568  106       NJET=NJET+1
569            XR1=XR1-ALOG(RLU_HIJING(0))
570            IF(XR1.LT.TTRIG) GO TO 106
571            XR=0.0
572  107       NJET=NJET+1
573            XR=XR-ALOG(RLU_HIJING(0))
574            IF(XR.LT.TT-TTRIG) GO TO 107
575            NJET=NJET-1
576            GO TO 112
577         ENDIF
578 C               ********create a hard interaction with specified P_T
579 c                                when IHPR2(3)>0
580         IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
581 C               ********create at least one pair of mini jets 
582 C                       when IHPR2(9)=1
583 C
584         IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
585      &          (1.0-EXP(-TTS))) GO TO 160
586 C               ********this is the probability for no jet production
587 110     XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
588 111     NJET=NJET+1
589         XR=XR-ALOG(RLU_HIJING(0))
590         IF(XR.LT.TT) GO TO 111
591 112     NJET=MIN(NJET,IHPR2(8))
592         IF(IHPR2(8).LT.0)  NJET=ABS(IHPR2(8))
593 C               ******** Determine number of mini jet production
594 C
595         DO 150 I_JET=1,NJET
596            CALL JETINI(JP,JT,0)
597            CALL HIJHRD(JP,JT,JOUT,JFLG,1)
598 C               ********JFLG=1 jets valence quarks, JFLG=2 with 
599 C                       gluon jet, JFLG=3 with q-qbar prod for
600 C                       (JP,JT). If JFLG=0 jets can not be produced 
601 C                       this time. If JFLG=-1, error occured abandon
602 C                       this event. JOUT is the total hard scat for
603 C                       (JP,JT) up to now.
604            IF(JFLG.EQ.0) GO TO 160
605            IF(JFLG.LT.0) THEN
606               IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
607               GO TO 50
608            ENDIF
609            JOUT=JOUT+1
610            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
611            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
612            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
613      &                  JFLG.GE.3) IASG(NSG,3)=1
614 C               ******** jet with PT>HIPR1(11) will be quenched
615  150    CONTINUE
616  160    CONTINUE
617         CALL HIJSFT(JP,JT,JOUT,IERROR)
618         IF(IERROR.NE.0) THEN
619            IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
620            GO TO 50
621         ENDIF
622 C
623 C               ********conduct soft scattering between JP and JT
624         JATT=JATT+JOUT
625
626 200     CONTINUE
627
628 c
629 c**************************
630 c
631         DO 201 JP=1,IHNT2(1)
632 c           write(6,*) JP, NFP(JP,3), NFP(JP,4), NFP(JP,5)
633            IF(NFP(JP,5).GT.2) THEN
634               NINP=NINP+1
635            ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
636               NELP=NELP+1
637            ENDIF
638
639            IF(NFP(JP,5).LE.2) THEN
640               IF (NFP(JP,3) .EQ. 2212) THEN
641                  NPSPECP = NPSPECP + 1
642               ELSE IF (NFP(JP,3) .EQ. 2112) THEN
643                  NNSPECP = NNSPECP + 1
644               ENDIF
645            ENDIF
646  201    continue
647         DO 202 JT=1,IHNT2(3)
648            IF(NFT(JT,5).GT.2) THEN
649               NINT=NINT+1
650            ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
651               NELT=NELT+1
652            ENDIF
653
654            IF(NFT(JT,5).LE.2) THEN
655               IF (NFT(JT,3) .EQ. 2212) THEN
656                  NPSPECT = NPSPECT + 1
657               ELSE IF (NFT(JT,3) .EQ. 2112) THEN
658                  NNSPECT = NNSPECT + 1
659               ENDIF
660            ENDIF
661  202    continue
662 c     
663 c*******************************
664
665 C********perform jet quenching for jets with PT>HIPR1(11)**********
666
667         IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
668      &                  IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
669                 DO 271 I=1,IHNT2(1)
670                         IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
671 271             CONTINUE
672                 DO 272 I=1,IHNT2(3)
673                         IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
674 272             CONTINUE
675                 DO 273 ISG=1,NSG
676                         IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
677 273             CONTINUE
678         ENDIF
679 C
680 C**************fragment all the string systems in the following*****
681 C
682 C********N_ST is where particle information starts
683 C********N_STR+1 is the number of strings in fragmentation
684 C********the number of strings before a line is stored in K(I,4)
685 C********IDSTR is id number of the string system (91,92 or 93)
686 C
687 c
688         IF(IHPR2(20).NE.0) THEN
689            DO 360 ISG=1,NSG
690                 CALL HIJFRG(ISG,3,IERROR)
691                 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
692                    MSTU(24)=0
693                    MSTU(28)=0
694                    IF(IHPR2(10).NE.0) THEN
695                       call LULIST_HIJING(1)
696                       WRITE(6,*) 'error occured, repeat the event'
697                    ENDIF
698                    GO TO 50
699                 ENDIF
700 C                       ********Check errors
701 C
702                 N_ST=1
703                 IDSTR=92
704                 IF(IHPR2(21).EQ.0) THEN
705                    CALL LUEDIT_HIJING(2)
706                 ELSE
707 351                N_ST=N_ST+1
708                    IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  351
709                    IDSTR=K(N_ST,2)
710                    N_ST=N_ST+1
711                 ENDIF
712 C
713                 IF(FRAME.EQ.'LAB') THEN
714                         CALL HIBOOST
715                 ENDIF
716 C               ******** boost back to lab frame(if it was in)
717 C
718                 N_STR=0
719                 DO 360 I=N_ST,N
720                    IF(K(I,2).EQ.IDSTR) THEN
721                       N_STR=N_STR+1
722                       GO TO 360
723                    ENDIF
724                    K(I,4)=N_STR
725                    NATT=NATT+1
726                    KATT(NATT,1)=K(I,2)
727                    KATT(NATT,2)=20
728                    KATT(NATT,4)=K(I,1)
729                    IF(K(I,3).EQ.0 ) THEN
730                       KATT(NATT,3)=0
731                    ELSE
732                       IF(K(K(I,3),2).EQ.IDSTR) THEN
733                          KATT(NATT,3)=0
734                       ELSE
735                          KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
736                       ENDIF
737                    ENDIF
738 C       ****** identify the mother particle
739                    PATT(NATT,1)=P(I,1)
740                    PATT(NATT,2)=P(I,2)
741                    PATT(NATT,3)=P(I,3)
742                    PATT(NATT,4)=P(I,4)
743                    VATT(NATT,1)=V(I,1)
744                    VATT(NATT,2)=V(I,2)
745                    VATT(NATT,3)=V(I,3)
746                    VATT(NATT,4)=V(I,4)
747
748                    EATT=EATT+P(I,4)
749 360        CONTINUE
750 C               ********Fragment the q-qbar jets systems *****
751 C
752            JTP(1)=IHNT2(1)
753            JTP(2)=IHNT2(3)
754            DO 400 NTP=1,2
755            DO 400 J_JTP=1,JTP(NTP)
756                 CALL HIJFRG(J_JTP,NTP,IERROR)
757                 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
758                    MSTU(24)=0
759                    MSTU(28)=0
760                    IF(IHPR2(10).NE.0) THEN
761                       call LULIST_HIJING(1)
762                       WRITE(6,*) 'error occured, repeat the event'
763                    ENDIF
764                    GO TO 50
765                 ENDIF
766 C                       ********check errors
767 C
768                 N_ST=1
769                 IDSTR=92
770
771                 NFTP=NFP(J_JTP,5)
772                 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
773
774                 IF(IHPR2(21).EQ.0) THEN
775                    CALL LUEDIT_HIJING(2)
776                 ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN
777 381                N_ST=N_ST+1
778                    IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  381
779                    IDSTR=K(N_ST,2)
780                    N_ST=N_ST+1
781                 ENDIF
782                 IF(FRAME.EQ.'LAB') THEN
783                         CALL HIBOOST
784                 ENDIF
785 C               ******** boost back to lab frame(if it was in)
786 C
787
788                 N_STR=0
789                 DO 390 I=N_ST,N
790                    IF(K(I,2).EQ.IDSTR) THEN
791                       N_STR=N_STR+1
792                       GO TO 390
793                    ENDIF
794                    K(I,4)=N_STR
795                    NATT=NATT+1
796                    KATT(NATT,1)=K(I,2)
797                    KATT(NATT,2)=NFTP
798                    KATT(NATT,4)=K(I,1)
799                    IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
800                       KATT(NATT,3)=0
801                    ELSE
802                       KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
803                    ENDIF
804 C       ****** identify the mother particle
805                    PATT(NATT,1)=P(I,1)
806                    PATT(NATT,2)=P(I,2)
807                    PATT(NATT,3)=P(I,3)
808                    PATT(NATT,4)=P(I,4)
809                    EATT=EATT+P(I,4)
810                    VATT(NATT,1)=V(I,1)  
811                    VATT(NATT,2)=V(I,2)  
812                    VATT(NATT,3)=V(I,3)  
813                    VATT(NATT,4)=V(I,4)
814 390             CONTINUE 
815 400        CONTINUE
816 C               ********Fragment the q-qq related string systems
817         ENDIF
818
819         DO 450 I=1,NDR
820                 NATT=NATT+1
821                 KATT(NATT,1)=KFDR(I)
822                 KATT(NATT,2)=40
823                 KATT(NATT,3)=0
824                 PATT(NATT,1)=PDR(I,1)
825                 PATT(NATT,2)=PDR(I,2)
826                 PATT(NATT,3)=PDR(I,3)
827                 PATT(NATT,4)=PDR(I,4)
828                 VATT(NATT,1)=V(I,1)  
829                 VATT(NATT,2)=V(I,2)  
830                 VATT(NATT,3)=V(I,3)  
831                 VATT(NATT,4)=V(I,4)
832                 EATT=EATT+PDR(I,4)
833 450     CONTINUE
834 C                       ********store the direct-produced particles
835 C
836         DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
837         IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
838      &     .AND.IHPR2(21).EQ.0) THEN
839         IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
840      &     event'
841 C               call LULIST_HIJING(1)
842                 GO TO 50
843         ENDIF
844
845         RETURN
846         END