1 c.................... hijing1.383_ampt.f
3 c The variables isng in HIJSFT and JL in ATTRAD were not initialized.
4 c The version initialize them. (as found by Fernando Marroquim)
9 c Nuclear distribution for deuteron is taken as the Hulthen wave
10 c function as provided by Brian Cole (Columbia)
11 clin used my own implementation of impact parameter
12 clin & proton-neutron distance within a deuteron.
17 c The parameters for Wood-Saxon distribution for deuteron are
18 c constrained to give the right rms ratius 2.116 fm
24 c The following common block is added to record the number of elastic
25 c (NELT, NELP) and inelastic (NINT, NINP) participants
27 c COMMON/HJGLBR/NELT,NINT,NELP,NINP
32 c A bug in the quenching subroutine is corrected. When calculating the
33 c distance between two wounded nucleons, the displacement of the
34 c impact parameter was not inculded. This bug was discovered by
35 c Dr. V.Uzhinskii JINR, Dubna, Russia
40 c Modification Oct. 8, 1998. In hijing, log(ran(nseed)) occasionally
41 c causes overfloat. It is modified to log(max(ran(nseed),1.0e-20)).
44 C Nothing important has been changed here. A few 'garbage' has been
45 C cleaned up here, like common block HJJET3 for the sea quark strings
46 C which were originally created to implement the DPM scheme which
47 C later was abadoned in the final version. The lines which operate
48 C on these data are also deleted in the program.
52 C There are some changes in the program: subroutine HARDJET is now
53 C consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
54 C for the triggered hard processes. Now that is done altogether
55 C with other normal hard processes in modified JETINI. In the new
56 C version one calls JETINI every time one calls HIJHRD. In the new
57 C version the effect of the isospin of the nucleon on hard processes,
58 C especially direct photons is correctly considered.
59 C For A+A collisions, one has to initilize pythia
60 C separately for each type of collisions, pp, pn,np and nn,
61 C or hp and hn for hA collisions. In JETINI we use the following
62 C catalogue for different types of collisions:
64 C h+A: h+p (itype=1), h+n (itype=2)
65 C A+h: p+h (itype=1), n+h (itype=2)
66 C A+A: p+p (itype=1), p+n (itype=2), n+p (itype=3), n+n (itype=4)
67 C*****************************************************************
71 C Last modification on January 5, 1998. Two mistakes are corrected in
72 C function G. A Mistake in the subroutine Parton is also corrected.
73 C (These are pointed out by Ysushi Nara).
76 C Last modifcation on April 10, 1996. To conduct final
77 C state radiation, PYTHIA reorganize the two scattered
78 C partons and their final momenta will be a little
79 C different. The summed total momenta of the partons
80 C from the final state radiation are stored in HINT1(26-29)
81 C and HINT1(36-39) which are little different from
82 C HINT1(21-24) and HINT1(41-44).
86 C Last modfication on September 11, 1995. When HIJING and
87 C PYTHIA are initialized, the shadowing is evaluated at
88 C b=0 which is the maximum. This will cause overestimate
89 C of shadowing for peripheral interactions. To correct this
90 C problem, shadowing is set to zero when initializing. Then
91 C use these maximum cross section without shadowing as a
92 C normalization of the Monte Carlo. This however increase
93 C the computing time. IHNT2(16) is used to indicate whether
94 C the sturcture function is called for (IHNT2(16)=1) initialization
95 C or for (IHNT2(16)=0)normal collisions simulation
97 C Last modification on Aagust 28, 1994. Two bugs associate
98 C with the impact parameter dependence of the shadowing is
102 c Last modification on October 14, 1994. One bug is corrected
103 c in the direct photon production option in subroutine
104 C HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
105 C Another bug associated with keeping the decay history
106 C in the particle information is also corrected.(this problem
107 C was reported by Matt Bloomer)
110 C Last modification on July 15, 1994. The option to trig on
111 C heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1)
112 C is added. To do this, set IHPR2(3)=3. For inclusive production,
113 C one should reset HIPR1(10)=0.0. One can also trig larger pt
114 C QQbar production by giving HIPR1(10) a nonvanishing value.
115 C The mass of the heavy quark in the calculation of the cross
116 C section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
117 C default is the charm mass D=1.5). We also include a separate
118 C K-factor for heavy quark and direct photon production by
121 C Last modification on May 24, 1994. The option to
122 C retain the information of all particles including those
123 C who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is
124 C added to contain the line number of the parent particle
125 C of the current line which is produced via a decay.
126 C KATT(I,4) is the status number of the particle: 11=particle
127 C which has decayed; 1=finally produced particle.
130 C Last modification on May 24, 1994( in HIJSFT when valence quark
131 C is quenched, the following error is corrected. 1.2*IHNT2(1) -->
132 C 1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
135 C Last modification on March 16, 1994 (heavy flavor production
136 C processes MSUB(81)=1 MSUB(82)=1 have been switched on,
137 C charm production is the default, B-quark option is
138 C IHPR2(18), when it is switched on, charm quark is
142 C Last modification on March 23, 1994 (an error is corrected
143 C in the impact parameter dependence of the jet cross section)
145 C Last modification Oct. 1993 to comply with non-vax
148 C*********************************************
149 C LAST MODIFICATION April 5, 1991
150 CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B
151 C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
152 C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
153 C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18),
154 C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
156 C gluon jets can form a single string system
158 C initial state radiation is included
160 C all QCD subprocesses are included
162 c direct particles production is included(currently only direct
165 C Effect of high P_T trigger bias on multiple jets distribution
167 C******************************************************************
169 C Heavy Ion Jet INteraction Generator *
171 C X. N. Wang and M. Gyulassy *
172 C Lawrence Berkeley Laboratory *
174 C******************************************************************
176 C******************************************************************
177 C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of *
178 C proj, NFP(K,4) original ID of proj. NFP(K,5)=colli status(0=no,*
179 C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
180 C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it *
181 C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
182 C (0=has not been,1=is going to be, -1=has already been scattered *
183 C NFP(k,11) total number of interactions this proj has suffered *
184 C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant *
185 C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and *
186 C diquark,PP(K,10)=PT of the hard scattering between the valence *
187 C quarks; PP(K,14,15)=the mass of quark,diquark. *
188 C******************************************************************
190 C****************************************************************
194 C****************************************************************
195 SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
197 cgsfs Added following for consistency with AMPT call
198 double precision BMIN0, BMAX0
201 PARAMETER (MAXPTN=400001)
202 clin-4/20/01 PARAMETER (MAXSTR = 1600)
203 PARAMETER (MAXSTR=150001)
206 PARAMETER (MAXIDL=4001)
209 DOUBLE PRECISION GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
210 DOUBLE PRECISION GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
211 DOUBLE PRECISION ATAUI, ZT1, ZT2, ZT3
212 DOUBLE PRECISION xnprod,etprod,xnfrz,etfrz,
213 & dnprod,detpro,dnfrz,detfrz
218 DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
219 & IPCOL(90000),ITCOL(90000)
220 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
223 COMMON/hjcrdn/YP(3,300),YT(3,300)
225 clin-7/16/03 NINT is a intrinsic fortran function, rename it to NINTHJ
226 c COMMON/HJGLBR/NELT,NINT,NELP,NINP
227 COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
229 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
232 c COMMON/HMAIN2/KATT(130000,4),PATT(130000,4)
233 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
235 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
237 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
238 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
239 & PJPM(300,500),NTJ(300),KFTJ(300,500),
240 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
241 & PJTE(300,500),PJTM(300,500)
244 c COMMON/HJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
245 c & K2SG(900,100),PXSG(900,100),PYSG(900,100),
246 c & PZSG(900,100),PESG(900,100),PMSG(900,100)
247 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
248 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
249 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
251 COMMON/HJJET4/NDR,IADR(MAXSTR,2),KFDR(MAXSTR),PDR(MAXSTR,5)
253 c common/xydr/rtdr(900,2)
254 common/xydr/rtdr(MAXSTR,2)
259 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
261 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
264 clin-9/29/03 changed name in order to distinguish from /prec2/
265 COMMON /ARPRC/ ITYPAR(MAXSTR),
266 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
267 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
270 c COMMON /ARPRC/ ITYP(MAXSTR),
271 c & GX(MAXSTR), GY(MAXSTR), GZ(MAXSTR), FT(MAXSTR),
272 c & PX(MAXSTR), PY(MAXSTR), PZ(MAXSTR), EE(MAXSTR),
280 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
281 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
282 & XMASS0(MAXPTN), ITYP0(MAXPTN)
284 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
285 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
286 & XMASS5(MAXPTN), ITYP5(MAXPTN)
288 COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
290 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
292 COMMON /SREC1/ NSP, NST, NSI
294 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
299 COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
300 & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
303 common/anim/nevent,isoft,isflag,izpc
306 DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
307 1 GXSGS,GYSGS,GZSGS,FTSGS
308 COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
309 & PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
310 & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
311 & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
313 clin-4/26/01 lepton and photon info:
314 COMMON /NOPREC/ NNOZPC, ITYPN(MAXIDL),
315 & GXN(MAXIDL), GYN(MAXIDL), GZN(MAXIDL), FTN(MAXIDL),
316 & PXN(MAXIDL), PYN(MAXIDL), PZN(MAXIDL), EEN(MAXIDL),
320 common /lastt/itimeh,bimp
322 COMMON /AREVT/ IAEVT, IARUN, MISS
323 common/phidcy/iphidcy,pttrig,ntrig,maxmiss
324 cwei DOUBLE PRECISION PATT
327 cgsfs WRITE(*,*) "IN Hijing, FRAME=",FRAME
328 cgsfs WRITE(*,*) "IN Hijing, BMIN=",BMIN0
329 cgsfs WRITE(*,*) "IN Hijing, BMAX=",BMAX0
331 BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
333 IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
335 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
337 C ********HIPR1(31) is in mb =0.1fm**2
338 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
339 C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
341 cgsfs WRITE(*,*) "IN Hijing, Modified BMIN=",BMIN
342 cgsfs WRITE(*,*) "IN Hijing, Modified BMAX=",BMAX
346 IF(IHNT2(1).LE.1) GO TO 14
352 C ********choose theta from uniform cos(theta) distr
353 PHI=RANART(NSEED)*2.0*HIPR1(40)
354 C ********choose phi form uniform phi distr 0 to 2*pi
355 YP(1,KP)=R*SX*COS(PHI)
356 YP(2,KP)=R*SX*SIN(PHI)
358 IF(HIPR1(29).EQ.0.0) GO TO 10
360 DNBP1=(YP(1,KP)-YP(1,KP2))**2
361 DNBP2=(YP(2,KP)-YP(2,KP2))**2
362 DNBP3=(YP(3,KP)-YP(3,KP2))**2
363 DNBP=DNBP1+DNBP2+DNBP3
364 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
365 C ********two neighbors cannot be closer than
370 clin-1/27/03 Hulthen wavefn for deuteron borrowed from hijing1.382.f,
371 c but modified [divide by 2, & x(p)=-x(n)]:
372 c (Note: hijing1.383.f has corrected this bug in hijing1.382.f)
373 if(IHNT2(1).EQ.2) then
374 rnd1=max(RANART(NSEED),1.0e-20)
375 rnd2=max(RANART(NSEED),1.0e-20)
376 rnd3=max(RANART(NSEED),1.0e-20)
377 R=-(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
378 & +4.38*0.85*log(rnd3)/(4.38+0.85))
382 PHI=RANART(NSEED)*2.0*HIPR1(40)
383 c R above is the relative distance between p & n in a deuteron:
385 YP(1,1)=R*SX*COS(PHI)
386 YP(2,1)=R*SX*SIN(PHI)
388 c p & n has opposite coordinates in the deuteron frame:
396 IF(YP(3,I).GT.YP(3,J)) GO TO 12
408 C******************************
412 IF(IHNT2(3).LE.1) GO TO 24
418 C ********choose theta from uniform cos(theta) distr
419 PHI=RANART(NSEED)*2.0*HIPR1(40)
420 C ********chose phi form uniform phi distr 0 to 2*pi
421 YT(1,KT)=R*SX*COS(PHI)
422 YT(2,KT)=R*SX*SIN(PHI)
424 IF(HIPR1(29).EQ.0.0) GO TO 20
426 DNBT1=(YT(1,KT)-YT(1,KT2))**2
427 DNBT2=(YT(2,KT)-YT(2,KT2))**2
428 DNBT3=(YT(3,KT)-YT(3,KT2))**2
429 DNBT=DNBT1+DNBT2+DNBT3
430 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
431 C ********two neighbors cannot be closer than
436 clin-1/27/03 Hulthen wavefn for deuteron borrowed from hijing1.382.f,
437 c but modified [divide by 2, & x(p)=-x(n)]:
438 if(IHNT2(3).EQ.2) then
439 rnd1=max(RANART(NSEED),1.0e-20)
440 rnd2=max(RANART(NSEED),1.0e-20)
441 rnd3=max(RANART(NSEED),1.0e-20)
442 R=-(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
443 & +4.38*0.85*log(rnd3)/(4.38+0.85))
447 PHI=RANART(NSEED)*2.0*HIPR1(40)
449 YT(1,1)=R*SX*COS(PHI)
450 YT(2,1)=R*SX*SIN(PHI)
459 IF(YT(3,I).LT.YT(3,J)) GO TO 22
471 C********************
476 c IF(MISS.GT.50) THEN
477 IF(MISS.GT.maxmiss) THEN
478 WRITE(6,*) 'infinite loop happened in HIJING'
490 C ********Initialize for a new event
504 C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
505 C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
506 C BY THE ANGLE PHI FOR EACH COLLISION.******************
508 BB=SQRT(BMIN**2+RANART(NSEED)*(BMAX**2-BMIN**2))
510 PHI=2.0*HIPR1(40)*RANART(NSEED)
521 B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
522 R2=B2*HIPR1(40)/HIPR1(31)/0.1
523 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
524 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
525 & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
526 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
527 & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
528 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
530 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
532 HINT1(18)=HINT1(14)-APHX1*HINT1(15)
533 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
534 IF(IHPR2(14).EQ.0.OR.
535 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
536 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
538 IF(RANTOT.GT.GS) GO TO 70
541 GSTOT0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
542 & /HIPR1(31)/2.0*ROMG(0.0)))
544 GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
545 GSTOT=2.0*(1.0-SQRT(1.0-GS))
546 RANTOT=RANART(NSEED)*GSTOT0
547 IF(RANTOT.GT.GSTOT) GO TO 70
548 IF(RANTOT.GT.GS) THEN
551 C ********perform elastic collisions
555 SJIP(JP,JT)=HINT1(18)
560 C ********total number interactions proj and targ has
563 clin-5/22/01 write impact parameter:
565 c write(6,*) '#impact parameter,nlop,ncolt=',bimp,nlop,ncolt
570 & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
573 C ********At large impact parameter, there maybe no
574 C interaction at all. For NN collision
575 C repeat the event until interaction happens
577 IF(IHPR2(3).NE.0) THEN
578 NHARD=1+INT(RANART(NSEED)*(NCOLT-1)+0.5)
579 NHARD=MIN(NHARD,NCOLT)
582 clin-6/2009 ctest off:
583 c write(99,*) IAEVT,NHARD,NCOLT,JPHARD,JTHARD
586 IF(IHPR2(9).EQ.1) THEN
587 NMINI=1+INT(RANART(NSEED)*(NCOLT-1)+0.5)
588 NMINI=MIN(NMINI,NCOLT)
592 C ********Specifying the location of the hard and
593 C minijet if they are enforced by user
597 IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
598 NFP(JP,11)=NFP(JP,11)+1
599 NFT(JT,11)=NFT(JT,11)+1
600 IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
603 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
606 ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
610 ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
616 C*****************************************************************
617 IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
618 C ********When IHPR2(8)=0 no jets are produced
619 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
620 C ********jets can not be produced for (JP,JT)
621 C because not enough energy avaible for
624 HINT1(18)=SJIP(JP,JT)
625 TT=ROMG(R2)*HINT1(18)/HIPR1(31)
626 TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
629 IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
631 CALL HIJHRD(JP,JT,0,JFLG,0)
641 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
642 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
643 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
644 & JFLG.GE.3) IASG(NSG,3)=1
648 HINT1(20+I05)=HINT1(40+I05)
649 HINT1(30+I05)=HINT1(50+I05)
651 clin-6/2009 ctest off:
652 c write(99,*) jp,jt,IHPR2(3),HIPR1(10),njet,
653 c 1 ihnt2(9),hint1(21),hint1(22),hint1(23),
654 c 2 ihnt2(10),hint1(31),hint1(32),hint1(33)
657 IF(IHPR2(8).EQ.0) GO TO 160
658 RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
659 & /REAL(IHNT2(1))**0.6666667,1.0)
660 RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
661 & /REAL(IHNT2(3))**0.6666667,1.0)
662 APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
664 APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
666 HINT1(65)=HINT1(61)-APHX1*HINT1(62)
667 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
668 TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
670 C ********subtract the trigger jet from total number
671 C of jet production to be done since it has
672 C already been produced here
673 XR1=-ALOG(EXP(-TTRIG)+RANART(NSEED)*(1.0-EXP(-TTRIG)))
675 XR1=XR1-ALOG(max(RANART(NSEED),1.0e-20))
676 IF(XR1.LT.TTRIG) GO TO 106
679 XR=XR-ALOG(max(RANART(NSEED),1.0e-20))
680 IF(XR.LT.TT-TTRIG) GO TO 107
684 C ********create a hard interaction with specified P_T
686 IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
687 C ********create at least one pair of mini jets
690 IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
691 & (1.0-EXP(-TTS))) GO TO 160
692 C ********this is the probability for no jet production
693 110 XR=-ALOG(EXP(-TT)+RANART(NSEED)*(1.0-EXP(-TT)))
695 XR=XR-ALOG(max(RANART(NSEED),1.0e-20))
696 IF(XR.LT.TT) GO TO 111
697 112 NJET=MIN(NJET,IHPR2(8))
698 IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
699 C ******** Determine number of mini jet production
703 CALL HIJHRD(JP,JT,JOUT,JFLG,1)
704 C ********JFLG=1 jets valence quarks, JFLG=2 with
705 C gluon jet, JFLG=3 with q-qbar prod for
706 C (JP,JT). If JFLG=0 jets can not be produced
707 C this time. If JFLG=-1, error occured abandon
708 C this event. JOUT is the total hard scat for
710 IF(JFLG.EQ.0) GO TO 160
712 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
716 IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
717 IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
718 IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
719 & JFLG.GE.3) IASG(NSG,3)=1
720 C ******** jet with PT>HIPR1(11) will be quenched
724 CALL HIJSFT(JP,JT,JOUT,IERROR)
726 IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
730 C ********conduct soft scattering between JP and JT
734 c**************************
736 clin-6/2009 write out initial minijet information:
738 if(pttrig.gt.0.and.ntrig.eq.0) goto 50
739 clin-6/2009 write out initial transverse positions of initial nucleons:
740 c write(94,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
743 c write(94,203) YP(1,JP)+0.5*BB, YP(2,JP), JP, NFP(JP,5)
744 IF(NFP(JP,5).GT.2) THEN
746 ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
751 clin-6/2009 target nucleon # has a minus sign for distinction from projectile:
752 c write(94,203) YT(1,JT)-0.5*BB, YT(2,JT), -JT, NFT(JT,5)
753 IF(NFT(JT,5).GT.2) THEN
755 ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
759 203 format(f10.3,1x,f10.3,2(1x,I5))
761 c*******************************
764 C********perform jet quenching for jets with PT>HIPR1(11)**********
766 IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
767 & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
769 IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
772 IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
775 IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
779 clin*****4/09/01-soft1, default way of treating strings:
781 clin-4/16/01 allow fragmentation:
785 c.....transfer data from HIJING to ZPC
791 DO 1008 I = 1, IHNT2(1)
793 DO 1007 J = 1, NPJ(I)
795 c.....for now only consider gluon cascade
796 IF (KFPJ(I, J) .EQ. 21) THEN
802 ITYP0(NPAR) = KFPJ(I, J)
804 clin-7/20/01 add dble or sngl to make precisions consistent
805 c GX0(NPAR) = YP(1, I)
806 GX0(NPAR) = dble(YP(1, I) + 0.5 * BB)
808 GY0(NPAR) = dble(YP(2, I))
811 PX0(NPAR) = dble(PJPX(I, J))
812 PY0(NPAR) = dble(PJPY(I, J))
813 PZ0(NPAR) = dble(PJPZ(I, J))
814 XMASS0(NPAR) = dble(PJPM(I, J))
815 c E0(NPAR) = dble(PJPE(I, J))
816 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
817 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
821 c.....end gluon selection
826 DO 1010 I = 1, IHNT2(3)
828 DO 1009 J = 1, NTJ(I)
830 c.....for now only consider gluon cascade
831 IF (KFTJ(I, J) .EQ. 21) THEN
836 ITYP0(NPAR) = KFTJ(I, J)
838 clin-7/20/01 add dble or sngl to make precisions consistent
839 c GX0(NPAR) = YT(1, I)
840 GX0(NPAR) = dble(YT(1, I) - 0.5 * BB)
842 GY0(NPAR) = dble(YT(2, I))
845 PX0(NPAR) = dble(PJTX(I, J))
846 PY0(NPAR) = dble(PJTY(I, J))
847 PZ0(NPAR) = dble(PJTZ(I, J))
848 XMASS0(NPAR) = dble(PJTM(I, J))
849 c E0(NPAR) = dble(PJTE(I, J))
850 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
851 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
854 c.....end gluon selection
861 DO 1011 J = 1, NJSG(I)
863 c.....for now only consider gluon cascade
864 IF (K2SG(I, J) .EQ. 21) THEN
869 ITYP0(NPAR) = K2SG(I, J)
870 clin-7/20/01 add dble or sngl to make precisions consistent:
872 1 dble(YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))
874 2 dble(YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))
877 PX0(NPAR) = dble(PXSG(I, J))
878 PY0(NPAR) = dble(PYSG(I, J))
879 PZ0(NPAR) = dble(PZSG(I, J))
880 XMASS0(NPAR) = dble(PMSG(I, J))
881 c E0(NPAR) = dble(PESG(I, J))
882 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
883 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
885 c.....end gluon selection
897 if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
898 c.....call ZPC for parton cascade
901 c write out parton and wounded nucleon information to ana/zpc1.mom:
903 c WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
904 c WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
906 cc WRITE (14, 411) PX5(I), PY5(I), PZ5(I), ITYP5(I),
908 if(dmax1(abs(GX5(I)),abs(GY5(I)),abs(GZ5(I)),abs(FT5(I)))
910 c write(14,210) ITYP5(I), PX5(I), PY5(I), PZ5(I), XMASS5(I),
911 c 1 GX5(I), GY5(I), GZ5(I), FT5(I)
913 c change format for large numbers:
914 c write(14,211) ITYP5(I), PX5(I), PY5(I), PZ5(I), XMASS5(I),
915 c 1 GX5(I), GY5(I), GZ5(I), FT5(I)
919 210 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
920 211 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
921 395 format(3I8,f10.4,4I5)
925 c 411 FORMAT(1X, 3F10.3, I6, 2F10.3)
928 clin-5/2009 ctest off:
931 c.....transfer data back from ZPC to HIJING
933 IF (LSTRG1(I) .LE. NSP) THEN
936 KFPJ(NSTRG, NPART) = ITYP5(I)
937 clin-7/20/01 add dble or sngl to make precisions consistent
938 PJPX(NSTRG, NPART) = sngl(PX5(I))
939 PJPY(NSTRG, NPART) = sngl(PY5(I))
940 PJPZ(NSTRG, NPART) = sngl(PZ5(I))
941 PJPE(NSTRG, NPART) = sngl(E5(I))
942 PJPM(NSTRG, NPART) = sngl(XMASS5(I))
943 ELSE IF (LSTRG1(I) .LE. NSP + NST) THEN
944 NSTRG = LSTRG1(I) - NSP
946 KFTJ(NSTRG, NPART) = ITYP5(I)
947 PJTX(NSTRG, NPART) = sngl(PX5(I))
948 PJTY(NSTRG, NPART) = sngl(PY5(I))
949 PJTZ(NSTRG, NPART) = sngl(PZ5(I))
950 PJTE(NSTRG, NPART) = sngl(E5(I))
951 PJTM(NSTRG, NPART) = sngl(XMASS5(I))
953 NSTRG = LSTRG1(I) - NSP - NST
955 K2SG(NSTRG, NPART) = ITYP5(I)
956 PXSG(NSTRG, NPART) = sngl(PX5(I))
957 PYSG(NSTRG, NPART) = sngl(PY5(I))
958 PZSG(NSTRG, NPART) = sngl(PZ5(I))
959 PESG(NSTRG, NPART) = sngl(E5(I))
960 PMSG(NSTRG, NPART) = sngl(XMASS5(I))
969 clin*****4/09/01-soft2, put q+dq+X in strings into ZPC:
970 elseif(isoft.eq.2) then
978 clin No fragmentation to hadrons, only on parton level,
979 c and transfer minijet and string data from HIJING to ZPC:
981 clin-4/12/01 forbid soft radiation before ZPC to avoid small-mass strings,
982 c and forbid jet order reversal before ZPC to avoid unphysical flavors:
986 IF(IHPR2(20).NE.0) THEN
988 DO 310 jjtp=1,IHNT2(2*NTP-1)
990 c change: do gluon kink only once: either here or in fragmentation.
991 CALL HIJFRG(jjtp,NTP,IERROR)
995 NPJ(jjtp)=MAX0(N-2,0)
997 clin-4/12/01: NPJ(jjtp)=MAX0(ipartn-2,0)
1000 NTJ(jjtp)=MAX0(N-2,0)
1001 clin-4/12/01: NTJ(jjtp)=MAX0(ipartn-2,0)
1008 ITYP0(NPAR) = K(II,2)
1011 clin-7/20/01 add dble or sngl to make precisions consistent
1012 PX0(NPAR) = dble(P(II,1))
1013 PY0(NPAR) = dble(P(II,2))
1014 PZ0(NPAR) = dble(P(II,3))
1015 XMASS0(NPAR) = dble(P(II,5))
1016 c E0(NPAR) = dble(P(II,4))
1017 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
1018 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
1019 IF (NTP .EQ. 1) THEN
1020 clin-7/20/01 add dble or sngl to make precisions consistent
1021 GX0(NPAR) = dble(YP(1, jjtp)+0.5 * BB)
1022 GY0(NPAR) = dble(YP(2, jjtp))
1025 if(IITYP.eq.2112.or.IITYP.eq.2212) then
1026 elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1027 1 (II.eq.1.or.II.eq.N)) then
1028 PP(nstrg,6)=sngl(PX0(NPAR))
1029 PP(nstrg,7)=sngl(PY0(NPAR))
1030 PP(nstrg,14)=sngl(XMASS0(NPAR))
1031 elseif((IITYP.eq.1103.or.IITYP.eq.2101
1032 1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1033 2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1034 3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1035 4 .and.(II.eq.1.or.II.eq.N)) then
1036 PP(nstrg,8)=sngl(PX0(NPAR))
1037 PP(nstrg,9)=sngl(PY0(NPAR))
1038 PP(nstrg,15)=sngl(XMASS0(NPAR))
1040 NPART = LPART0(NPAR)-1
1041 KFPJ(NSTRG, NPART) = ITYP0(NPAR)
1042 PJPX(NSTRG, NPART) = sngl(PX0(NPAR))
1043 PJPY(NSTRG, NPART) = sngl(PY0(NPAR))
1044 PJPZ(NSTRG, NPART) = sngl(PZ0(NPAR))
1045 PJPE(NSTRG, NPART) = sngl(E0(NPAR))
1046 PJPM(NSTRG, NPART) = sngl(XMASS0(NPAR))
1049 GX0(NPAR) = dble(YT(1, jjtp)-0.5 * BB)
1050 GY0(NPAR) = dble(YT(2, jjtp))
1052 nstrg=LSTRG0(NPAR)-NSP
1053 if(IITYP.eq.2112.or.IITYP.eq.2212) then
1054 elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1055 1 (II.eq.1.or.II.eq.N)) then
1056 PT(nstrg,6)=sngl(PX0(NPAR))
1057 PT(nstrg,7)=sngl(PY0(NPAR))
1058 PT(nstrg,14)=sngl(XMASS0(NPAR))
1059 elseif((IITYP.eq.1103.or.IITYP.eq.2101
1060 1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1061 2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1062 3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1063 4 .and.(II.eq.1.or.II.eq.N)) then
1064 PT(nstrg,8)=sngl(PX0(NPAR))
1065 PT(nstrg,9)=sngl(PY0(NPAR))
1066 PT(nstrg,15)=sngl(XMASS0(NPAR))
1068 NPART = LPART0(NPAR)-1
1069 KFTJ(NSTRG, NPART) = ITYP0(NPAR)
1070 PJTX(NSTRG, NPART) = sngl(PX0(NPAR))
1071 PJTY(NSTRG, NPART) = sngl(PY0(NPAR))
1072 PJTZ(NSTRG, NPART) = sngl(PZ0(NPAR))
1073 PJTE(NSTRG, NPART) = sngl(E0(NPAR))
1074 PJTM(NSTRG, NPART) = sngl(XMASS0(NPAR))
1082 CALL HIJFRG(ISG,3,IERROR)
1091 ITYP0(NPAR) = K(II,2)
1093 1 dble(YP(1,IASG(ISG,1))+YT(1,IASG(ISG,2)))
1095 2 dble(YP(2,IASG(ISG,1))+YT(2,IASG(ISG,2)))
1098 PX0(NPAR) = dble(P(II,1))
1099 PY0(NPAR) = dble(P(II,2))
1100 PZ0(NPAR) = dble(P(II,3))
1101 XMASS0(NPAR) = dble(P(II,5))
1102 c E0(NPAR) = dble(P(II,4))
1103 E0(NPAR) = dsqrt(PX0(NPAR)**2+PY0(NPAR)**2
1104 1 +PZ0(NPAR)**2+XMASS0(NPAR)**2)
1114 if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
1115 c.....call ZPC for parton cascade
1119 c WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
1120 WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
1124 c WRITE (14, 311) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1125 c & XMASS5(I), E5(I)
1126 WRITE (14, 312) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1127 & XMASS5(I), E5(I),LSTRG1(I), LPART1(I)
1129 c 311 FORMAT(1X, 3F10.4, I6, 2F10.4)
1130 312 FORMAT(1X, 3F10.3, I6, 2F10.3,1X,I6,1X,I3)
1133 clin-5/2009 ctest off:
1136 clin-4/13/01 initialize four momenta and invariant mass of strings after ZPC:
1149 IF (LSTRG1(I) .LE. NSP) THEN
1151 c nucleons without interactions:
1152 if(IITYP.eq.2112.or.IITYP.eq.2212) then
1153 clin-7/20/01 add dble or sngl to make precisions consistent
1154 PP(nstrg,1)=sngl(PX5(I))
1155 PP(nstrg,2)=sngl(PY5(I))
1156 PP(nstrg,3)=sngl(PZ5(I))
1157 PP(nstrg,4)=sngl(E5(I))
1158 PP(nstrg,5)=sngl(XMASS5(I))
1160 elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1161 1 (LPART1(I).eq.1.or.LPART1(I).eq.(NPJ(NSTRG)+2))) then
1162 PP(nstrg,6)=sngl(PX5(I))
1163 PP(nstrg,7)=sngl(PY5(I))
1164 PP(nstrg,14)=sngl(XMASS5(I))
1165 PP(nstrg,1)=PP(nstrg,1)+sngl(PX5(I))
1166 PP(nstrg,2)=PP(nstrg,2)+sngl(PY5(I))
1167 PP(nstrg,3)=PP(nstrg,3)+sngl(PZ5(I))
1168 PP(nstrg,4)=PP(nstrg,4)+sngl(E5(I))
1169 PP(nstrg,5)=sqrt(PP(nstrg,4)**2-PP(nstrg,1)**2
1170 1 -PP(nstrg,2)**2-PP(nstrg,3)**2)
1172 elseif((IITYP.eq.1103.or.IITYP.eq.2101
1173 1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1174 2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1175 3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1176 4 .and.(LPART1(I).eq.1.or.LPART1(I).eq.(NPJ(NSTRG)+2))) then
1177 PP(nstrg,8)=sngl(PX5(I))
1178 PP(nstrg,9)=sngl(PY5(I))
1179 PP(nstrg,15)=sngl(XMASS5(I))
1180 PP(nstrg,1)=PP(nstrg,1)+sngl(PX5(I))
1181 PP(nstrg,2)=PP(nstrg,2)+sngl(PY5(I))
1182 PP(nstrg,3)=PP(nstrg,3)+sngl(PZ5(I))
1183 PP(nstrg,4)=PP(nstrg,4)+sngl(E5(I))
1184 PP(nstrg,5)=sqrt(PP(nstrg,4)**2-PP(nstrg,1)**2
1185 1 -PP(nstrg,2)**2-PP(nstrg,3)**2)
1186 c partons in projectile or target strings:
1189 KFPJ(NSTRG, NPART) = ITYP5(I)
1190 PJPX(NSTRG, NPART) = sngl(PX5(I))
1191 PJPY(NSTRG, NPART) = sngl(PY5(I))
1192 PJPZ(NSTRG, NPART) = sngl(PZ5(I))
1193 PJPE(NSTRG, NPART) = sngl(E5(I))
1194 PJPM(NSTRG, NPART) = sngl(XMASS5(I))
1196 ELSE IF (LSTRG1(I) .LE. NSP + NST) THEN
1197 NSTRG = LSTRG1(I) - NSP
1198 if(IITYP.eq.2112.or.IITYP.eq.2212) then
1199 PT(nstrg,1)=sngl(PX5(I))
1200 PT(nstrg,2)=sngl(PY5(I))
1201 PT(nstrg,3)=sngl(PZ5(I))
1202 PT(nstrg,4)=sngl(E5(I))
1203 PT(nstrg,5)=sngl(XMASS5(I))
1204 elseif((IITYP.eq.1.or.IITYP.eq.2).and.
1205 1 (LPART1(I).eq.1.or.LPART1(I).eq.(NTJ(NSTRG)+2))) then
1206 PT(nstrg,6)=sngl(PX5(I))
1207 PT(nstrg,7)=sngl(PY5(I))
1208 PT(nstrg,14)=sngl(XMASS5(I))
1209 PT(nstrg,1)=PT(nstrg,1)+sngl(PX5(I))
1210 PT(nstrg,2)=PT(nstrg,2)+sngl(PY5(I))
1211 PT(nstrg,3)=PT(nstrg,3)+sngl(PZ5(I))
1212 PT(nstrg,4)=PT(nstrg,4)+sngl(E5(I))
1213 PT(nstrg,5)=sqrt(PT(nstrg,4)**2-PT(nstrg,1)**2
1214 1 -PT(nstrg,2)**2-PT(nstrg,3)**2)
1215 elseif((IITYP.eq.1103.or.IITYP.eq.2101
1216 1 .or.IITYP.eq.2103.or.IITYP.eq.2203.
1217 2 .or.IITYP.eq.3101.or.IITYP.eq.3103.
1218 3 .or.IITYP.eq.3201.or.IITYP.eq.3203.or.IITYP.eq.3303)
1219 4 .and.(LPART1(I).eq.1.or.LPART1(I).eq.(NTJ(NSTRG)+2))) then
1220 PT(nstrg,8)=sngl(PX5(I))
1221 PT(nstrg,9)=sngl(PY5(I))
1222 PT(nstrg,15)=sngl(XMASS5(I))
1223 PT(nstrg,1)=PT(nstrg,1)+sngl(PX5(I))
1224 PT(nstrg,2)=PT(nstrg,2)+sngl(PY5(I))
1225 PT(nstrg,3)=PT(nstrg,3)+sngl(PZ5(I))
1226 PT(nstrg,4)=PT(nstrg,4)+sngl(E5(I))
1227 PT(nstrg,5)=sqrt(PT(nstrg,4)**2-PT(nstrg,1)**2
1228 1 -PT(nstrg,2)**2-PT(nstrg,3)**2)
1231 KFTJ(NSTRG, NPART) = ITYP5(I)
1232 PJTX(NSTRG, NPART) = sngl(PX5(I))
1233 PJTY(NSTRG, NPART) = sngl(PY5(I))
1234 PJTZ(NSTRG, NPART) = sngl(PZ5(I))
1235 PJTE(NSTRG, NPART) = sngl(E5(I))
1236 PJTM(NSTRG, NPART) = sngl(XMASS5(I))
1239 NSTRG = LSTRG1(I) - NSP - NST
1241 K2SG(NSTRG, NPART) = ITYP5(I)
1242 PXSG(NSTRG, NPART) = sngl(PX5(I))
1243 PYSG(NSTRG, NPART) = sngl(PY5(I))
1244 PZSG(NSTRG, NPART) = sngl(PZ5(I))
1245 PESG(NSTRG, NPART) = sngl(E5(I))
1246 PMSG(NSTRG, NPART) = sngl(XMASS5(I))
1251 clin-4/09/01 turn on fragmentation with soft radiation
1252 c and jet order reversal to form hadrons after ZPC:
1256 clin-4/13/01 allow small mass strings (D=1.5GeV):
1263 clin-4/19/01-soft3, fragment strings, then convert hadrons to partons
1265 elseif(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
1266 clin-4/24/01 normal fragmentation first:
1269 IF(IHPR2(20).NE.0) THEN
1271 CALL HIJFRG(ISG,3,IERROR)
1275 IF(IHPR2(21).EQ.0) THEN
1279 IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO 551
1284 IF(FRAME.EQ.'LAB') THEN
1287 C ******** boost back to lab frame(if it was in)
1291 IF(K(I,2).EQ.IDSTR) THEN
1300 c from Yasushi, to avoid violation of array limits:
1301 c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1302 clin-4/2008 to avoid out-of-bound error in K():
1303 c IF(K(I,3).EQ.0 .OR.
1304 c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1306 IF(K(I,3).EQ.0) THEN
1308 ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1312 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1315 C ****** identify the mother particle
1321 GXAR(NATT) = 0.5 * (YP(1, IASG(ISG, 1)) +
1322 & YT(1, IASG(ISG, 2)))
1323 GYAR(NATT) = 0.5 * (YP(2, IASG(ISG, 1)) +
1324 & YT(2, IASG(ISG, 2)))
1327 ITYPAR(NATT) = K(I, 2)
1328 PXAR(NATT) = P(I, 1)
1329 PYAR(NATT) = P(I, 2)
1330 PZAR(NATT) = P(I, 3)
1331 PEAR(NATT) = P(I, 4)
1332 XMAR(NATT) = P(I, 5)
1336 C ********Fragment the q-qbar jets systems *****
1341 DO 600 jjtp=1,JTP(NTP)
1342 CALL HIJFRG(jjtp,NTP,IERROR)
1346 IF(IHPR2(21).EQ.0) THEN
1350 IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO 581
1354 IF(FRAME.EQ.'LAB') THEN
1357 C ******** boost back to lab frame(if it was in)
1360 IF(NTP.EQ.2) NFTP=10+NFT(jjtp,5)
1363 IF(K(I,2).EQ.IDSTR) THEN
1372 c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1374 c IF(K(I,3).EQ.0 .OR.
1375 c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1377 IF(K(I,3).EQ.0) THEN
1379 ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1383 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1386 C ****** identify the mother particle
1392 IF (NTP .EQ. 1) THEN
1393 GXAR(NATT) = YP(1, jjtp)+0.5 * BB
1394 GYAR(NATT) = YP(2, jjtp)
1396 GXAR(NATT) = YT(1, jjtp)-0.5 * BB
1397 GYAR(NATT) = YT(2, jjtp)
1401 ITYPAR(NATT) = K(I, 2)
1402 PXAR(NATT) = P(I, 1)
1403 PYAR(NATT) = P(I, 2)
1404 PZAR(NATT) = P(I, 3)
1405 PEAR(NATT) = P(I, 4)
1406 XMAR(NATT) = P(I, 5)
1411 C ********Fragment the q-qq related string systems
1413 clin-4/2008 check for zero NDR value:
1418 KATT(NATT,1)=KFDR(I)
1421 PATT(NATT,1)=PDR(I,1)
1422 PATT(NATT,2)=PDR(I,2)
1423 PATT(NATT,3)=PDR(I,3)
1424 PATT(NATT,4)=PDR(I,4)
1426 clin-11/11/03 set direct photons positions and time at formation:
1427 GXAR(NATT) = rtdr(I,1)
1428 GYAR(NATT) = rtdr(I,2)
1431 ITYPAR(NATT) =KATT(NATT,1)
1432 PXAR(NATT) = PATT(NATT,1)
1433 PYAR(NATT) = PATT(NATT,2)
1434 PZAR(NATT) = PATT(NATT,3)
1435 PEAR(NATT) = PATT(NATT,4)
1436 XMAR(NATT) = PDR(I,5)
1440 clin-6/2009 ctest on:
1445 clin-4/19/01 convert hadrons to partons for ZPC (with GX0 given):
1448 clin-7/03/01 move up, used in zpstrg (otherwise not set and incorrect):
1456 if(ioscar.eq.3) WRITE (95, *) IAEVT, mul
1457 c.....call ZPC for parton cascade
1460 c WRITE (14, 395) ITEST, MUL, bimp, NELP,NINP,NELT,NINTHJ
1461 WRITE (14, 395) IAEVT, MISS, MUL, bimp, NELP,NINP,NELT,NINTHJ
1465 c WRITE (14, 511) PX5(I), PY5(I), PZ5(I), ITYP5(I),
1466 c & XMASS5(I), E5(I)
1467 WRITE (14, 512) ITYP5(I), PX5(I), PY5(I), PZ5(I),
1468 & XMASS5(I), LSTRG1(I), LPART1(I), FT5(I)
1471 c 511 FORMAT(1X, 3F10.4, I6, 2F10.4)
1472 512 FORMAT(I6,4(1X,F10.3),1X,I6,1X,I3,1X,F10.3)
1473 c 513 FORMAT(1X, 4F10.4)
1475 clin-5/2009 ctest off:
1478 clin save data after ZPC for fragmentation purpose:
1479 c.....transfer data back from ZPC to HIJING
1480 DO 1018 I = 1, MAXSTR
1499 K2SGS(NSTRG, NPART) = ITYP5(I)
1500 PXSGS(NSTRG, NPART) = PX5(I)
1501 PYSGS(NSTRG, NPART) = PY5(I)
1502 PZSGS(NSTRG, NPART) = PZ5(I)
1503 PMSGS(NSTRG, NPART) = XMASS5(I)
1504 clin-7/20/01 E5(I) does no include the finite parton mass XMASS5(I),
1505 c so define it anew:
1506 c PESGS(NSTRG, NPART) = E5(I)
1507 c if(abs(PZ5(i)/E5(i)).gt.0.9999999d0)
1508 c 1 write(91,*) 'a',PX5(i),PY5(i),XMASS5(i),PZ5(i),E5(i)
1509 E5(I)=dsqrt(PX5(I)**2+PY5(I)**2+PZ5(I)**2+XMASS5(I)**2)
1510 PESGS(NSTRG, NPART) = E5(I)
1511 c if(abs(PZ5(i)/E5(i)).gt.0.9999999d0)
1512 c 1 write(91,*) 'b: new E5(I)=',E5(i)
1514 GXSGS(NSTRG, NPART) = GX5(I)
1515 GYSGS(NSTRG, NPART) = GY5(I)
1516 GZSGS(NSTRG, NPART) = GZ5(I)
1517 FTSGS(NSTRG, NPART) = FT5(I)
1527 C**************fragment all the string systems in the following*****
1529 C********nsbst is where particle information starts
1530 C********nsbstR+1 is the number of strings in fragmentation
1531 C********the number of strings before a line is stored in K(I,4)
1532 C********IDSTR is id number of the string system (91,92 or 93)
1534 clin-4/30/01 convert partons to hadrons after ZPC:
1535 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
1541 KATT(NATT,1)=ITYPN(I)
1551 ITYPAR(NATT)=ITYPN(I)
1561 IF(IHPR2(20).NE.0) THEN
1563 CALL HIJFRG(ISG,3,IERROR)
1564 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
1567 IF(IHPR2(10).NE.0) THEN
1569 WRITE(6,*) 'error occured ISG, repeat the event'
1575 C ********Check errors
1579 IF(IHPR2(21).EQ.0) THEN
1583 IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO 351
1588 IF(FRAME.EQ.'LAB') THEN
1591 C ******** boost back to lab frame(if it was in)
1595 IF(K(I,2).EQ.IDSTR) THEN
1604 c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1606 c IF(K(I,3).EQ.0 .OR.
1607 c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1609 IF(K(I,3).EQ.0) THEN
1611 ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1615 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1618 C ****** identify the mother particle
1627 c GXAR(NATT) = 0.5 * (YP(1, IASG(ISG, 1)) +
1628 c & YT(1, IASG(ISG, 2)))
1629 c GYAR(NATT) = 0.5 * (YP(2, IASG(ISG, 1)) +
1630 c & YT(2, IASG(ISG, 2)))
1631 LSG = NSP + NST + ISG
1632 GXAR(NATT) = sngl(ZT1(LSG))
1633 GYAR(NATT) = sngl(ZT2(LSG))
1634 GZAR(NATT) = sngl(ZT3(LSG))
1635 FTAR(NATT) = sngl(ATAUI(LSG))
1637 ITYPAR(NATT) = K(I, 2)
1638 PXAR(NATT) = P(I, 1)
1639 PYAR(NATT) = P(I, 2)
1640 PZAR(NATT) = P(I, 3)
1641 PEAR(NATT) = P(I, 4)
1642 XMAR(NATT) = P(I, 5)
1646 C ********Fragment the q-qbar jets systems *****
1651 DO 400 jjtp=1,JTP(NTP)
1652 CALL HIJFRG(jjtp,NTP,IERROR)
1653 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
1656 IF(IHPR2(10).NE.0) THEN
1658 WRITE(6,*) 'error occured P&T, repeat the event'
1660 clin-6/2009 when this happens, the event will be repeated,
1661 c and another record for the same event number will be written into
1662 c zpc.dat, zpc.res, minijet-initial-beforePropagation.dat,
1663 c parton-initial-afterPropagation.dat, parton-after-coalescence.dat,
1664 c and parton-collisionsHistory.dat.
1668 C ********check errors
1672 IF(IHPR2(21).EQ.0) THEN
1676 IF(K(nsbst,2).LT.91.OR.K(nsbst,2).GT.93) GO TO 381
1680 IF(FRAME.EQ.'LAB') THEN
1683 C ******** boost back to lab frame(if it was in)
1686 IF(NTP.EQ.2) NFTP=10+NFT(jjtp,5)
1689 IF(K(I,2).EQ.IDSTR) THEN
1698 c IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
1700 c IF(K(I,3).EQ.0 .OR.
1701 c 1 (K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR)) THEN
1703 IF(K(I,3).EQ.0) THEN
1705 ELSEIF(K(I,3).ne.0.and.K(K(I,3),2).EQ.IDSTR) THEN
1709 KATT(NATT,3)=NATT-I+K(I,3)+nsbstR-K(K(I,3),4)
1711 C ****** identify the mother particle
1719 c IF (NTP .EQ. 1) THEN
1720 c GXAR(NATT) = YP(1, jjtp)
1722 c GXAR(NATT) = YT(1, jjtp)
1724 c IF (NTP .EQ. 1) THEN
1725 c GYAR(NATT) = YP(2, jjtp)
1727 c GYAR(NATT) = YT(2, jjtp)
1729 IF (NTP .EQ. 1) THEN
1734 GXAR(NATT) = sngl(ZT1(LSG))
1735 GYAR(NATT) = sngl(ZT2(LSG))
1736 GZAR(NATT) = sngl(ZT3(LSG))
1737 FTAR(NATT) = sngl(ATAUI(LSG))
1739 ITYPAR(NATT) = K(I, 2)
1740 PXAR(NATT) = P(I, 1)
1741 PYAR(NATT) = P(I, 2)
1742 PZAR(NATT) = P(I, 3)
1743 PEAR(NATT) = P(I, 4)
1744 XMAR(NATT) = P(I, 5)
1749 C ********Fragment the q-qq related string systems
1754 KATT(NATT,1)=KFDR(I)
1757 PATT(NATT,1)=PDR(I,1)
1758 PATT(NATT,2)=PDR(I,2)
1759 PATT(NATT,3)=PDR(I,3)
1760 PATT(NATT,4)=PDR(I,4)
1762 clin-11/11/03 set direct photons positions and time at formation:
1763 GXAR(NATT) = rtdr(I,1)
1764 GYAR(NATT) = rtdr(I,2)
1767 ITYPAR(NATT) =KATT(NATT,1)
1768 PXAR(NATT) = PATT(NATT,1)
1769 PYAR(NATT) = PATT(NATT,2)
1770 PZAR(NATT) = PATT(NATT,3)
1771 PEAR(NATT) = PATT(NATT,4)
1772 XMAR(NATT) = PDR(I,5)
1775 C ********store the direct-produced particles
1781 DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
1782 IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
1783 & .AND.IHPR2(21).EQ.0) THEN
1785 & WRITE(6,*) 'Energy not conserved, repeat the event'
1787 write(6,*) 'violated:EATT,NATT,B=',EATT,NATT,bimp
1790 c write(6,*) 'satisfied:EATT,NATT,B=',EATT,NATT,bimp
1798 SUBROUTINE HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)
1799 CHARACTER FRAME*4,PROJ*4,TARG*4,EFRAME*4
1800 DOUBLE PRECISION DD1,DD2,DD3,DD4
1801 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
1803 COMMON/hjcrdn/YP(3,300),YT(3,300)
1805 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
1807 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
1809 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1811 EXTERNAL FNKICK,FNKC2,FNSTRU,FNSTRM,FNSTRS
1822 HINT1(8)=MAX(ULMASS(2112),ULMASS(2212))
1825 IF(PROJ.NE.'A') THEN
1826 IF(PROJ.EQ.'P') THEN
1828 ELSE IF(PROJ.EQ.'PBAR') THEN
1830 ELSE IF(PROJ.EQ.'PI+') THEN
1832 ELSE IF(PROJ.EQ.'PI-') THEN
1834 ELSE IF(PROJ.EQ.'K+') THEN
1836 ELSE IF(PROJ.EQ.'K-') THEN
1838 ELSE IF(PROJ.EQ.'N') THEN
1840 ELSE IF(PROJ.EQ.'NBAR') THEN
1843 WRITE(6,*) PROJ, 'wrong or unavailable proj name'
1846 HINT1(8)=ULMASS(IHNT2(5))
1848 IF(TARG.NE.'A') THEN
1849 IF(TARG.EQ.'P') THEN
1851 ELSE IF(TARG.EQ.'PBAR') THEN
1853 ELSE IF(TARG.EQ.'PI+') THEN
1855 ELSE IF(TARG.EQ.'PI-') THEN
1857 ELSE IF(TARG.EQ.'K+') THEN
1859 ELSE IF(TARG.EQ.'K-') THEN
1861 ELSE IF(TARG.EQ.'N') THEN
1863 ELSE IF(TARG.EQ.'NBAR') THEN
1866 WRITE(6,*) TARG,'wrong or unavailable targ name'
1869 HINT1(9)=ULMASS(IHNT2(6))
1872 C...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
1873 IF(IHPR2(12).GT.0) THEN
1874 CALL LUGIVE('MDCY(C221,1)=0')
1875 clin-11/07/00 no K* decays:
1876 CALL LUGIVE('MDCY(C313,1)=0')
1877 CALL LUGIVE('MDCY(C-313,1)=0')
1878 CALL LUGIVE('MDCY(C323,1)=0')
1879 CALL LUGIVE('MDCY(C-323,1)=0')
1880 clin-1/04/01 no K0 and K0bar decays so K0L and K0S do not appear,
1881 c this way the K/Kbar difference is accounted for exactly:
1882 CALL LUGIVE('MDCY(C311,1)=0')
1883 CALL LUGIVE('MDCY(C-311,1)=0')
1884 clin-11/08/00 no Delta decays:
1885 CALL LUGIVE('MDCY(C1114,1)=0')
1886 CALL LUGIVE('MDCY(C2114,1)=0')
1887 CALL LUGIVE('MDCY(C2214,1)=0')
1888 CALL LUGIVE('MDCY(C2224,1)=0')
1889 CALL LUGIVE('MDCY(C-1114,1)=0')
1890 CALL LUGIVE('MDCY(C-2114,1)=0')
1891 CALL LUGIVE('MDCY(C-2214,1)=0')
1892 CALL LUGIVE('MDCY(C-2224,1)=0')
1895 CALL LUGIVE('MDCY(C213,1)=0')
1896 CALL LUGIVE('MDCY(C-213,1)=0')
1897 CALL LUGIVE('MDCY(C113,1)=0')
1898 CALL LUGIVE('MDCY(C223,1)=0')
1899 CALL LUGIVE('MDCY(C333,1)=0')
1901 CALL LUGIVE('MDCY(C111,1)=0')
1902 CALL LUGIVE('MDCY(C310,1)=0')
1903 CALL LUGIVE('MDCY(C411,1)=0;MDCY(C-411,1)=0')
1904 CALL LUGIVE('MDCY(C421,1)=0;MDCY(C-421,1)=0')
1905 CALL LUGIVE('MDCY(C431,1)=0;MDCY(C-431,1)=0')
1906 CALL LUGIVE('MDCY(C511,1)=0;MDCY(C-511,1)=0')
1907 CALL LUGIVE('MDCY(C521,1)=0;MDCY(C-521,1)=0')
1908 CALL LUGIVE('MDCY(C531,1)=0;MDCY(C-531,1)=0')
1909 CALL LUGIVE('MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
1910 CALL LUGIVE('MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
1911 CALL LUGIVE('MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
1912 CALL LUGIVE('MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
1913 CALL LUGIVE('MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
1914 CALL LUGIVE('MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
1915 CALL LUGIVE('MDCY(C3334,1)=0;MDCY(C-3334,1)=0')
1919 IF(IHPR2(10).EQ.0) THEN
1925 clin parj(41) and (42) are a, b parameters in Lund, read from input.ampt:
1931 clin 2 popcorn parameters read from input.ampt:
1936 clin parj(21) gives the mean gaussian width for hadron Pt:
1938 clin parj(2) is gamma_s=P(s)/P(u), kappa propto 1/b/(2+a) assumed.
1939 rkp=HIPR1(4)*(2+HIPR1(3))/PARJ(42)/(2+PARJ(41))
1940 PARJ(2)=PARJ(2)**(1./rkp)
1941 PARJ(21)=PARJ(21)*sqrt(rkp)
1942 clin-10/31/00 update when string tension is changed:
1945 C ******** set up for jetset
1946 IF(FRAME.EQ.'LAB') THEN
1950 HINT1(1)=SQRT(HINT1(8)**2+2.0*HINT1(9)*EFRM+HINT1(9)**2)
1951 DD4=DSQRT(DD1**2-DD2**2)/(DD1+DD3)
1953 HINT1(3)=0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1954 DD4=DSQRT(DD1**2-DD2**2)/DD1
1955 HINT1(4)=0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1959 ELSE IF(FRAME.EQ.'CMS') THEN
1966 DD4=DSQRT(1.D0-4.D0*DD2**2/DD1**2)
1967 HINT1(4)=0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1968 DD4=DSQRT(1.D0-4.D0*DD3**2/DD1**2)
1969 HINT1(5)=-0.5*sngl(DLOG((1.D0+DD4)/(1.D0-DD4)))
1970 HINT1(6)=HINT1(1)/2.0
1971 HINT1(7)=HINT1(1)/2.0
1973 C ********define Lorentz transform to lab frame
1975 C ********calculate the cross sections involved with
1976 C nucleon collisions.
1977 IF(IHNT2(1).GT.1) THEN
1978 CALL HIJWDS(IHNT2(1),1,RMAX)
1980 C ********set up Wood-Sax distr for proj.
1982 IF(IHNT2(3).GT.1) THEN
1983 CALL HIJWDS(IHNT2(3),2,RMAX)
1985 C ********set up Wood-Sax distr for targ.
1991 IF(I.EQ.10) GO TO 30
1992 IF(HIDAT0(10,I).LE.HINT1(1)) GO TO 20
1995 HIDAT(J)=HIDAT0(J,I-1)+(HIDAT0(J,I)-HIDAT0(J,I-1))
1996 & *(HINT1(1)-HIDAT0(10,I-1))/(HIDAT0(10,I)-HIDAT0(10,I-1))
1999 HIPR1(30)=2.0*HIDAT(5)
2004 IF(IHPR2(5).NE.0) THEN
2005 CALL HIFUN(3,0.0,36.0,FNKICK)
2006 C ********booking for generating pt**2 for pt kick
2008 CALL HIFUN(7,0.0,6.0,FNKC2)
2009 CALL HIFUN(4,0.0,1.0,FNSTRU)
2010 CALL HIFUN(5,0.0,1.0,FNSTRM)
2011 CALL HIFUN(6,0.0,1.0,FNSTRS)
2012 C ********booking for x distribution of valence quarks
2014 IF(FRAME.EQ.'LAB') EFRAME='Elab'
2015 WRITE(6,100) EFRAME,EFRM,PROJ,IHNT2(1),IHNT2(2),
2016 & TARG,IHNT2(3),IHNT2(4)
2018 & 10X,'**************************************************'/
2020 & 10X,'* HIJING has been initialized at *'/
2021 & 10X,'*',13X,A4,'= ',F10.2,' GeV/n',13X,'*'/
2023 & 10X,'*',8X,'for ',
2024 & A4,'(',I3,',',I3,')',' + ',A4,'(',I3,',',I3,')',7X,'*'/
2025 & 10X,'**************************************************')
2032 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2035 FNKICK=1.0/(X+HIPR1(19)**2)/(X+HIPR1(20)**2)
2036 & /(1+EXP((SQRT(X)-HIPR1(20))/0.4))
2042 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2045 FNKC2=X*EXP(-2.0*X/HIPR1(42))
2052 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2055 FNSTRU=(1.0-X)**HIPR1(44)/
2056 & (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
2063 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2066 FNSTRM=1.0/((1.0-X)**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
2067 & /(X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
2073 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2076 FNSTRS=(1.0-X)**HIPR1(47)/
2077 & (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(48)
2085 IMPLICIT DOUBLE PRECISION(D)
2086 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
2088 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2090 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2094 DBETA=dble(P(I,3)/P(I,4))
2095 IF(ABS(DBETA).GE.1.D0) THEN
2097 IF(DB.GT.0.99999999D0) THEN
2098 C ********Rescale boost vector if too close to unity.
2099 WRITE(6,*) '(HIBOOT:) boost vector too large'
2102 DGA=1D0/SQRT(1D0-DB**2)
2105 P(I,3)=sngl((DP3+DB*DP4)*DGA)
2106 P(I,4)=sngl((DP4+DB*DP3)*DGA)
2109 Y=0.5*sngl(DLOG((1.D0+DBETA)/(1.D0-DBETA)))
2110 AMT=SQRT(P(I,1)**2+P(I,2)**2+P(I,5)**2)
2111 P(I,3)=AMT*SINH(Y+HINT1(3))
2112 P(I,4)=AMT*COSH(Y+HINT1(3))
2120 SUBROUTINE QUENCH(JPJT,NTP)
2121 PARAMETER (MAXSTR=150001)
2122 DIMENSION RDP(300),LQP(300),RDT(300),LQT(300)
2123 COMMON/hjcrdn/YP(3,300),YT(3,300)
2125 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2128 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2129 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2130 & PJPM(300,500),NTJ(300),KFTJ(300,500),
2131 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2132 & PJTE(300,500),PJTM(300,500)
2134 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2135 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2136 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2138 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2150 IF(NTP.EQ.2) GO TO 400
2151 IF(NTP.EQ.3) GO TO 2000
2152 C*******************************************************
2153 C Jet interaction for proj jet in the direction PHIP
2154 C******************************************************
2156 IF(NFP(JPJT,7).NE.1) RETURN
2160 PTJET0=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2)
2161 IF(PTJET0.LE.HIPR1(11)) GO TO 290
2162 PTOT=SQRT(PTJET0*PTJET0+PJPZ(JP,I)**2)
2163 IF(PTOT.LT.HIPR1(8)) GO TO 290
2164 PHIP=ULANGL(PJPX(JP,I),PJPY(JP,I))
2165 C******* find the wounded proj which can interact with jet***
2167 DO 100 I2=1,IHNT2(1)
2168 IF(NFP(I2,5).NE.3 .OR. I2.EQ.JP) GO TO 100
2169 DX=YP(1,I2)-YP(1,JP)
2170 DY=YP(2,I2)-YP(2,JP)
2174 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2175 IF(DPHI.GE.HIPR1(40)/2.0) GO TO 100
2176 RD0=SQRT(DX*DX+DY*DY)
2177 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 100
2180 RDP(KP)=COS(DPHI)*RD0
2182 C******* rearrange according decending rd************
2185 IF(RDP(I2).LT.RDP(J2)) GO TO 110
2193 C****** find wounded targ which can interact with jet********
2195 DO 120 I2=1,IHNT2(3)
2196 IF(NFT(I2,5).NE.3) GO TO 120
2197 DX=YT(1,I2)-YP(1,JP)-BBX
2198 DY=YT(2,I2)-YP(2,JP)-BBY
2202 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2203 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 120
2204 RD0=SQRT(DX*DX+DY*DY)
2205 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 120
2208 RDT(KT)=COS(DPHI)*RD0
2210 C******* rearrange according decending rd************
2213 IF(RDT(I2).LT.RDT(J2)) GO TO 130
2227 PTOT=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2+PJPZ(JP,I)**2)
2232 200 RN=RANART(NSEED)
2233 210 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 290
2234 IF(MT.GE.KT) GO TO 220
2235 IF(MP.GE.KP) GO TO 240
2236 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 240
2239 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
2241 IF(KFPJ(JP,I).NE.21) DP=0.5*DP
2242 C ********string tension of quark jet is 0.5 of gluon's
2243 IF(DP.LE.0.2) GO TO 210
2244 IF(PTOT.LE.0.4) GO TO 290
2245 IF(PTOT.LE.DP) DP=PTOT-0.2
2248 IF(KFPJ(JP,I).NE.21) THEN
2249 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
2251 DE=SQRT(PJPM(JP,I)**2+PTOT**2)
2252 & -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
2253 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
2255 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
2256 PP(LQP(MP),4)=SQRT(ERSHU)
2257 PP(LQP(MP),5)=SQRT(AMSHU)
2259 C ********reshuffle the energy when jet has mass
2264 C ********momentum and energy transfer from jet
2266 NPJ(LQP(MP))=NPJ(LQP(MP))+1
2267 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
2268 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
2269 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
2270 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
2271 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
2272 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
2277 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
2279 IF(DP.LE.0.2) GO TO 210
2280 IF(PTOT.LE.0.4) GO TO 290
2281 IF(PTOT.LE.DP) DP=PTOT-0.2
2284 IF(KFPJ(JP,I).NE.21) THEN
2285 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
2287 DE=SQRT(PJPM(JP,I)**2+PTOT**2)
2288 & -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
2289 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
2291 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
2292 PT(LQT(MT),4)=SQRT(ERSHU)
2293 PT(LQT(MT),5)=SQRT(AMSHU)
2295 C ********reshuffle the energy when jet has mass
2301 C ********momentum and energy transfer from jet
2302 NTJ(LQT(MT))=NTJ(LQT(MT))+1
2303 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
2304 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
2305 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
2306 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
2307 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
2308 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
2310 260 PJPX(JP,I)=(PTOT-DP)*V1
2311 PJPY(JP,I)=(PTOT-DP)*V2
2312 PJPZ(JP,I)=(PTOT-DP)*V3
2313 PJPE(JP,I)=PJPE(JP,I)-DE
2322 C*******************************************************
2323 C Jet interaction for target jet in the direction PHIT
2324 C******************************************************
2326 C******* find the wounded proj which can interact with jet***
2328 400 IF(NFT(JPJT,7).NE.1) RETURN
2331 PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
2332 IF(PTJET0.LE.HIPR1(11)) GO TO 690
2333 PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
2334 IF(PTOT.LT.HIPR1(8)) GO TO 690
2335 PHIT=ULANGL(PJTX(JT,I),PJTY(JT,I))
2337 DO 500 I2=1,IHNT2(1)
2338 IF(NFP(I2,5).NE.3) GO TO 500
2339 DX=YP(1,I2)+BBX-YT(1,JT)
2340 DY=YP(2,I2)+BBY-YT(2,JT)
2344 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2345 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 500
2346 RD0=SQRT(DX*DX+DY*DY)
2347 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 500
2350 RDP(KP)=COS(DPHI)*RD0
2352 C******* rearrange according to decending rd************
2355 IF(RDP(I2).LT.RDP(J2)) GO TO 510
2363 C****** find wounded targ which can interact with jet********
2365 DO 520 I2=1,IHNT2(3)
2366 IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 520
2367 DX=YT(1,I2)-YT(1,JT)
2368 DY=YT(2,I2)-YT(2,JT)
2372 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2373 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 520
2374 RD0=SQRT(DX*DX+DY*DY)
2375 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 520
2378 RDT(KT)=COS(DPHI)*RD0
2380 C******* rearrange according to decending rd************
2383 IF(RDT(I2).LT.RDT(J2)) GO TO 530
2397 PTOT=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2+PJTZ(JT,I)**2)
2402 600 RN=RANART(NSEED)
2403 610 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 690
2404 IF(MT.GE.KT) GO TO 620
2405 IF(MP.GE.KP) GO TO 640
2406 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 640
2409 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
2411 IF(KFTJ(JT,I).NE.21) DP=0.5*DP
2412 C ********string tension of quark jet is 0.5 of gluon's
2413 IF(DP.LE.0.2) GO TO 610
2414 IF(PTOT.LE.0.4) GO TO 690
2415 IF(PTOT.LE.DP) DP=PTOT-0.2
2418 IF(KFTJ(JT,I).NE.21) THEN
2419 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
2421 DE=SQRT(PJTM(JT,I)**2+PTOT**2)
2422 & -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
2423 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
2425 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
2426 PP(LQP(MP),4)=SQRT(ERSHU)
2427 PP(LQP(MP),5)=SQRT(AMSHU)
2429 C ********reshuffle the energy when jet has mass
2435 C ********momentum and energy transfer from jet
2436 NPJ(LQP(MP))=NPJ(LQP(MP))+1
2437 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
2438 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
2439 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
2440 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
2441 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
2442 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
2448 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
2450 IF(DP.LE.0.2) GO TO 610
2451 IF(PTOT.LE.0.4) GO TO 690
2452 IF(PTOT.LE.DP) DP=PTOT-0.2
2455 IF(KFTJ(JT,I).NE.21) THEN
2456 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
2458 DE=SQRT(PJTM(JT,I)**2+PTOT**2)
2459 & -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
2460 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
2462 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
2463 PT(LQT(MT),4)=SQRT(ERSHU)
2464 PT(LQT(MT),5)=SQRT(AMSHU)
2466 C ********reshuffle the energy when jet has mass
2472 C ********momentum and energy transfer from jet
2473 NTJ(LQT(MT))=NTJ(LQT(MT))+1
2474 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
2475 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
2476 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
2477 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
2478 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
2479 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
2481 660 PJTX(JT,I)=(PTOT-DP)*V1
2482 PJTY(JT,I)=(PTOT-DP)*V2
2483 PJTZ(JT,I)=(PTOT-DP)*V3
2484 PJTE(JT,I)=PJTE(JT,I)-DE
2491 C********************************************************
2492 C Q-QBAR jet interaction
2493 C********************************************************
2495 IF(IASG(ISG,3).NE.1) RETURN
2499 XJ=(YP(1,JP)+BBX+YT(1,JT))/2.0
2500 YJ=(YP(2,JP)+BBY+YT(2,JT))/2.0
2501 DO 2690 I=1,NJSG(ISG)
2502 PTJET0=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2)
2503 IF(PTJET0.LE.HIPR1(11).OR.PESG(ISG,I).LT.HIPR1(1))
2505 PTOT=SQRT(PTJET0*PTJET0+PZSG(ISG,I)**2)
2506 IF(PTOT.LT.MAX(HIPR1(1),HIPR1(8))) GO TO 2690
2507 PHIQ=ULANGL(PXSG(ISG,I),PYSG(ISG,I))
2509 DO 2500 I2=1,IHNT2(1)
2510 IF(NFP(I2,5).NE.3.OR.I2.EQ.JP) GO TO 2500
2516 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2517 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2500
2518 RD0=SQRT(DX*DX+DY*DY)
2519 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2500
2522 RDP(KP)=COS(DPHI)*RD0
2524 C******* rearrange according to decending rd************
2527 IF(RDP(I2).LT.RDP(J2)) GO TO 2510
2535 C****** find wounded targ which can interact with jet********
2537 DO 2520 I2=1,IHNT2(3)
2538 IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 2520
2544 IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI
2545 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2520
2546 RD0=SQRT(DX*DX+DY*DY)
2547 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2520
2550 RDT(KT)=COS(DPHI)*RD0
2552 C******* rearrange according to decending rd************
2555 IF(RDT(I2).LT.RDT(J2)) GO TO 2530
2569 PTOT=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2
2575 2600 RN=RANART(NSEED)
2576 2610 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 2690
2577 IF(MT.GE.KT) GO TO 2620
2578 IF(MP.GE.KP) GO TO 2640
2579 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 2640
2582 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
2583 DP=DRR*HIPR1(14)/2.0
2584 IF(DP.LE.0.2) GO TO 2610
2585 IF(PTOT.LE.0.4) GO TO 2690
2586 IF(PTOT.LE.DP) DP=PTOT-0.2
2589 IF(K2SG(ISG,I).NE.21) THEN
2590 IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
2591 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
2593 DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
2594 & -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
2595 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
2597 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
2598 PP(LQP(MP),4)=SQRT(ERSHU)
2599 PP(LQP(MP),5)=SQRT(AMSHU)
2601 C ********reshuffle the energy when jet has mass
2607 C ********momentum and energy transfer from jet
2608 NPJ(LQP(MP))=NPJ(LQP(MP))+1
2609 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
2610 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
2611 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
2612 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
2613 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
2614 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
2620 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
2622 IF(DP.LE.0.2) GO TO 2610
2623 IF(PTOT.LE.0.4) GO TO 2690
2624 IF(PTOT.LE.DP) DP=PTOT-0.2
2627 IF(K2SG(ISG,I).NE.21) THEN
2628 IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
2629 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
2631 DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
2632 & -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
2633 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
2635 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
2636 PT(LQT(MT),4)=SQRT(ERSHU)
2637 PT(LQT(MT),5)=SQRT(AMSHU)
2639 C ********reshuffle the energy when jet has mass
2645 C ********momentum and energy transfer from jet
2646 NTJ(LQT(MT))=NTJ(LQT(MT))+1
2647 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
2648 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
2649 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
2650 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
2651 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
2652 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
2654 2660 PXSG(ISG,I)=(PTOT-DP)*V1
2655 PYSG(ISG,I)=(PTOT-DP)*V2
2656 PZSG(ISG,I)=(PTOT-DP)*V3
2657 PESG(ISG,I)=PESG(ISG,I)-DE
2670 SUBROUTINE HIJFRG(JTP,NTP,IERROR)
2671 C NTP=1, fragment proj string, NTP=2, targ string,
2672 C NTP=3, independent
2673 C strings from jets. JTP is the line number of the string
2674 C*******Fragment all leadng strings of proj and targ**************
2675 C IHNT2(1)=atomic #, IHNT2(2)=proton #(=-1 if anti-proton) *
2676 C******************************************************************
2677 PARAMETER (MAXSTR=150001)
2678 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2680 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
2682 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
2684 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2685 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2686 & PJPM(300,500),NTJ(300),KFTJ(300,500),
2687 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2688 & PJTE(300,500),PJTM(300,500)
2690 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2691 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2692 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2695 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
2697 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2702 common/anim/nevent,isoft,isflag,izpc
2707 c.....set up fragmentation function according to the number of collisions
2708 c.....a wounded nucleon has suffered
2709 c IF (NTP .EQ. 1) THEN
2710 c NCOLL = NFP(JTP, 11)
2711 c ELSE IF (NTP .EQ. 2) THEN
2712 c NCOLL = NFT(JTP, 11)
2713 c ELSE IF (NTP .EQ. 3) THEN
2714 c NCOLL = (NFP(IASG(JTP,1), 11) + NFT(IASG(JTP,2), 11)) / 2
2716 c IF (NCOLL .LE. 1) THEN
2718 c ELSE IF (NCOLL .EQ. 2) THEN
2720 c ELSE IF (NCOLL .EQ. 3) THEN
2722 c ELSE IF (NCOLL .EQ. 4) THEN
2724 c ELSE IF (NCOLL .EQ. 5) THEN
2726 c ELSE IF (NCOLL .GE. 6) THEN
2735 C ********initialize the document lines
2739 DO 100 I=1,NJSG(ISG)
2749 C IF(IHPR2(1).GT.0) CALL ATTRAD(IERROR)
2750 c IF(IERROR.NE.0) RETURN
2752 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2756 IF(NTP.EQ.2) GO TO 200
2757 IF(JTP.GT.IHNT2(1)) RETURN
2758 IF(NFP(JTP,5).NE.3.AND.NFP(JTP,3).NE.0
2759 & .AND.NPJ(JTP).EQ.0.AND.NFP(JTP,10).EQ.0) GO TO 1000
2760 IF(NFP(JTP,15).EQ.-1) THEN
2780 C ********for NFP(JTP,15)=-1 NFP(JTP,1) IS IN -Z DIRECTION
2786 IF((ABS(PB1-PP(JTP,1)).GT.0.01.OR.
2787 & ABS(PB2-PP(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
2788 & WRITE(6,*) ' Pt of Q and QQ do not sum to the total',jtp
2789 & ,ntp,pq11,pq21,pb1,'*',pq12,pq22,pb2,'*',pp(JTP,1),pp(JTP,2)
2792 200 IF(JTP.GT.IHNT2(3)) RETURN
2793 IF(NFT(JTP,5).NE.3.AND.NFT(JTP,3).NE.0
2794 & .AND.NTJ(JTP).EQ.0.AND.NFT(JTP,10).EQ.0) GO TO 1200
2795 IF(NFT(JTP,15).EQ.1) THEN
2814 C ********for NFT(JTP,15)=1 NFT(JTP,1) IS IN +Z DIRECTION
2821 IF((ABS(PB1-PT(JTP,1)).GT.0.01.OR.
2822 & ABS(PB2-PT(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
2823 & WRITE(6,*) ' Pt of Q and QQ do not sum to the total',jtp
2824 & ,ntp,pq11,pq21,pb1,'*',pq12,pq22,pb2,'*',pt(JTP,1),pt(JTP,2)
2825 300 IF(PECM.LT.HIPR1(1)) THEN
2827 IF(IHPR2(10).EQ.0) RETURN
2828 WRITE(6,*) ' ECM=',PECM,' energy of the string is too small'
2830 write (6,*) 'JTP,NTP,pq=',JTP,NTP,pq11,pq12,pq21,pq22
2833 AMT=PECM**2+PB1**2+PB2**2
2834 AMT1=AM1**2+PQ11**2+PQ12**2
2835 AMT2=AM2**2+PQ21**2+PQ22**2
2836 PZCM=SQRT(ABS(AMT**2+AMT1**2+AMT2**2-2.0*AMT*AMT1
2837 & -2.0*AMT*AMT2-2.0*AMT1*AMT2))/2.0/SQRT(AMT)
2838 C *******PZ of end-partons in c.m. frame of the string
2844 P(1,4)=SQRT(AMT1+PZCM**2)
2851 P(2,4)=SQRT(AMT2+PZCM**2)
2855 CALL HIROBO(0.0,0.0,0.0,0.0,BTZ)
2857 IF((PQ21**2+PQ22**2).GT.(PQ11**2+PQ12**2)) THEN
2870 ELSE IF(NTP.EQ.2) THEN
2875 C*******************attach produced jets to the leadng partons****
2876 IF(NTP.EQ.1.AND.NPJ(JTP).NE.0) THEN
2878 C IF(NPJ(JTP).GE.2) CALL HIJSRT(JTP,1)
2879 C ********sort jets in order of y
2881 IF((ABS(KF1).GT.1000.AND.KF1.LT.0)
2882 & .OR.(ABS(KF1).LT.1000.AND.KF1.GT.0)) IEX=1
2897 clin-4/12/01: IF(IEX.EQ.1) I0=NPJ(JTP)-I+1
2898 IF(IEX.EQ.1.and.(isoft.ne.2.or.isflag.ne.0))
2900 C ********reverse the order of jets
2904 IF(KK1.NE.21 .AND. KK1.NE.0) K(I+1,1)=
2905 & 1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
2906 P(I+1,1)=PJPX(JTP,I0)
2907 P(I+1,2)=PJPY(JTP,I0)
2908 P(I+1,3)=PJPZ(JTP,I0)
2909 P(I+1,4)=PJPE(JTP,I0)
2910 P(I+1,5)=PJPM(JTP,I0)
2913 ELSE IF(NTP.EQ.2.AND.NTJ(JTP).NE.0) THEN
2915 c IF(NTJ(JTP).GE.2) CALL HIJSRT(JTP,2)
2916 C ********sort jets in order of y
2918 IF((ABS(KF2).GT.1000.AND.KF2.LT.0)
2919 & .OR.(ABS(KF2).LT.1000.AND.KF2.GT.0)) IEX=0
2933 clin-4/12/01: IF(IEX.EQ.1) I0=NTJ(JTP)-I+1
2934 IF(IEX.EQ.1.and.(isoft.ne.2.or.isflag.ne.0))
2936 C ********reverse the order of jets
2940 IF(KK1.NE.21 .AND. KK1.NE.0) K(I+1,1)=
2941 & 1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
2942 P(I+1,1)=PJTX(JTP,I0)
2943 P(I+1,2)=PJTY(JTP,I0)
2944 P(I+1,3)=PJTZ(JTP,I0)
2945 P(I+1,4)=PJTE(JTP,I0)
2946 P(I+1,5)=PJTM(JTP,I0)
2950 IF(IHPR2(1).GT.0.AND.RANART(NSEED).LE.HIDAT(3)) THEN
2953 IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
2955 IF(HINT1(1).GE.1000.0.AND.JETOT.EQ.0)THEN
2962 ELSE IF(JETOT.EQ.0.AND.IHPR2(1).GT.0.AND.
2963 & HINT1(1).GE.1000.0.AND.
2964 & RANART(NSEED).LE.0.8) THEN
2969 IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
2975 IF(IERROR.NE.0) RETURN
2976 C ******** conduct soft radiations
2977 C****************************
2982 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
2992 C ********proj remain as a nucleon or delta
2995 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
3006 C ********targ remain as a nucleon or delta
3009 if(isoft.ne.2.or.isflag.ne.0) CALL LUEXEC
3018 C****************************************************************
3019 C conduct soft radiation according to dipole approxiamtion
3020 C****************************************************************
3021 SUBROUTINE ATTRAD(IERROR)
3023 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3025 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3027 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3034 C.....S INVARIANT MASS-SQUARED BETWEEN PARTONS I AND I+1......
3035 C.....SM IS THE LARGEST MASS-SQUARED....
3040 S=2.*(P(I,4)*P(I+1,4)-P(I,1)*P(I+1,1)-P(I,2)*P(I+1,2)
3041 & -P(I,3)*P(I+1,3))+P(I,5)**2+P(I+1,5)**2
3043 WP=SQRT(S)-1.5*(P(I,5)+P(I+1,5))
3045 PBT1=P(I,1)+P(I+1,1)
3046 PBT2=P(I,2)+P(I+1,2)
3047 PBT3=P(I,3)+P(I+1,3)
3048 PBT4=P(I,4)+P(I+1,4)
3049 BTT=(PBT1**2+PBT2**2+PBT3**2)/PBT4**2
3050 IF(BTT.GE.1.0-1.0E-10) GO TO 30
3051 IF((I.NE.1.OR.I.NE.N-1).AND.
3052 & (K(I,2).NE.21.AND.K(I+1,2).NE.21)) GO TO 30
3057 S=(SM+1.5*(P(JL,5)+P(JL+1,5)))**2
3058 IF(SM.LT.HIPR1(5)) GOTO 2
3060 C.....MAKE PLACE FOR ONE GLUON.....
3061 IF(JL+1.EQ.N) GOTO 190
3070 C.....BOOST TO REST SYSTEM FOR PARTICLES JL AND JL+1.....
3071 P1=P(JL,1)+P(JL+1,1)
3072 P2=P(JL,2)+P(JL+1,2)
3073 P3=P(JL,3)+P(JL+1,3)
3074 P4=P(JL,4)+P(JL+1,4)
3080 CALL ATROBO(0.,0.,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
3081 IF(IERROR.NE.0) RETURN
3082 C.....ROTATE TO Z-AXIS....
3083 CTH=P(JL,3)/SQRT(P(JL,4)**2-P(JL,5)**2)
3084 IF(ABS(CTH).GT.1.0) CTH=MAX(-1.,MIN(1.,CTH))
3086 PHI=ULANGL(P(JL,1),P(JL,2))
3087 CALL ATROBO(0.,-PHI,0.,0.,0.,IMIN,IMAX,IERROR)
3088 CALL ATROBO(-THETA,0.,0.,0.,0.,IMIN,IMAX,IERROR)
3090 C.....CREATE ONE GLUON AND ORIENTATE.....
3092 1 CALL AR3JET(S,X1,X3,JL)
3093 CALL ARORIE(S,X1,X3,JL)
3094 IF(HIDAT(2).GT.0.0) THEN
3095 PTG1=SQRT(P(JL,1)**2+P(JL,2)**2)
3096 PTG2=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
3097 PTG3=SQRT(P(JL+2,1)**2+P(JL+2,2)**2)
3098 PTG=MAX(PTG1,PTG2,PTG3)
3099 IF(PTG.GT.HIDAT(2)) THEN
3100 FMFACT=EXP(-(PTG**2-HIDAT(2)**2)/HIPR1(2)**2)
3101 IF(RANART(NSEED).GT.FMFACT) GO TO 1
3104 C.....ROTATE AND BOOST BACK.....
3107 CALL ATROBO(THETA,PHI,-BEX,-BEY,-BEZ,IMIN,IMAX,IERROR)
3108 IF(IERROR.NE.0) RETURN
3109 C.....ENUMERATE THE GLUONS.....
3122 C----THETA FUNCTION DAMPING OF THE EMITTED GLUONS. FOR HADRON-HADRON.
3124 C IF(VFR(2).GT.0.) THEN
3125 C PTG=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
3127 C DOPT=SQRT((4.*PAR(71)*VFR(2))/WSTRI)
3128 C PTOPT=(DOPT*WSTRI)/(2.*VFR(2))
3129 C IF(PTG.GT.PTOPT) IORDER=IORDER-1
3130 C IF(PTG.GT.PTOPT) GOTO 1
3133 IF(SM.GE.HIPR1(5)) GOTO 40
3148 SUBROUTINE AR3JET(S,X1,X3,JL)
3150 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3152 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3159 IF(K(JL,2).NE.21 .AND. K(JL+1,2).NE.21) C=8./27.
3162 IF(K(JL,2).NE.21) EXP1=2
3163 IF(K(JL+1,2).NE.21) EXP3=2
3165 YMA=ALOG(.5/SQRT(A)+SQRT(.25/A-1))
3169 XT2M=(1.-2.*SQRT(SM1)+SM1-SM3)*(1.-2.*SQRT(SM3)-SM1+SM3)
3172 1 IF(NTRY.EQ.5000) THEN
3173 X1=.5*(2.*SQRT(SM1)+1.+SM1-SM3)
3174 X3=.5*(2.*SQRT(SM3)+1.-SM1+SM3)
3179 XT2=A*(XT2M/A)**(RANART(NSEED)**(1./D))
3181 YMAX=ALOG(.5/SQRT(XT2)+SQRT(.25/XT2-1.))
3182 Y=(2.*RANART(NSEED)-1.)*YMAX
3183 X1=1.-SQRT(XT2)*EXP(Y)
3184 X3=1.-SQRT(XT2)*EXP(-Y)
3187 IF(K(JL,2).NE.21 .OR. K(JL+1,2).NE.21) THEN
3188 IF((1.-X1)*(1.-X2)*(1.-X3)-X2*SM1*(1.-X1)-X2*SM3*(1.-X3).
3189 & LE.0..OR.X1.LE.2.*SQRT(SM1)-SM1+SM3.OR.X3.LE.2.*SQRT(SM3)
3196 FG=2.*YMAX*C*(X1**EXP1+X3**EXP3)/D
3198 IF(FG.LT.RANART(NSEED)) GOTO 1
3202 C*************************************************************
3205 SUBROUTINE ARORIE(S,X1,X3,JL)
3207 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3209 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3219 P1=SQRT(E1**2-P(JL,5)**2)
3220 P3=SQRT(E3**2-P(JL+1,5)**2)
3222 IF(P1.GT.0..AND.P3.GT.0.) CBET=(P(JL,5)**2
3223 & +P(JL+1,5)**2+2.*E1*E3-S*(1.-X2))/(2.*P1*P3)
3224 IF(ABS(CBET).GT.1.0) CBET=MAX(-1.,MIN(1.,CBET))
3227 C.....MINIMIZE PT1-SQUARED PLUS PT3-SQUARED.....
3229 PSI=.5*ULANGL(P1**2+P3**2*COS(2.*BET),-P3**2*SIN(2.*BET))
3234 ELSE IF(P3.GT.P1) THEN
3235 PSI=.5*ULANGL(P3**2+P1**2*COS(2.*BET),-P1**2*SIN(2.*BET))
3237 PZ1=-P1*COS(BET+PSI)
3242 DEL=2.0*HIPR1(40)*RANART(NSEED)
3244 P(JL,1)=PT1*SIN(DEL)
3245 P(JL,2)=-PT1*COS(DEL)
3248 P(JL+2,1)=PT3*SIN(DEL)
3249 P(JL+2,2)=-PT3*COS(DEL)
3252 P(JL+1,1)=-P(JL,1)-P(JL+2,1)
3253 P(JL+1,2)=-P(JL,2)-P(JL+2,2)
3254 P(JL+1,3)=-P(JL,3)-P(JL+2,3)
3260 C*******************************************************************
3261 C make boost and rotation to entries from IMIN to IMAX
3262 C*******************************************************************
3263 SUBROUTINE ATROBO(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
3264 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3266 DIMENSION ROT(3,3),PV(3)
3267 DOUBLE PRECISION DP(4),DBEX,DBEY,DBEZ,DGA,DGA2,DBEP,DGABEP
3271 IF(IMIN.LE.0 .OR. IMAX.GT.N .OR. IMIN.GT.IMAX) RETURN
3273 IF(THE**2+PHI**2.GT.1E-20) THEN
3274 C...ROTATE (TYPICALLY FROM Z AXIS TO DIRECTION THETA,PHI)
3275 ROT(1,1)=COS(THE)*COS(PHI)
3277 ROT(1,3)=SIN(THE)*COS(PHI)
3278 ROT(2,1)=COS(THE)*SIN(PHI)
3280 ROT(2,3)=SIN(THE)*SIN(PHI)
3285 C************** IF(MOD(K(I,1)/10000,10).GE.6) GOTO 120
3289 110 P(I,J)=ROT(J,1)*PV(1)+ROT(J,2)*PV(2)
3294 IF(BEX**2+BEY**2+BEZ**2.GT.1E-20) THEN
3295 C...LORENTZ BOOST (TYPICALLY FROM REST TO MOMENTUM/ENERGY=BETA)
3299 DGA2=1D0-DBEX**2-DBEY**2-DBEZ**2
3300 IF(DGA2.LE.0D0) THEN
3306 C************* IF(MOD(K(I,1)/10000,10).GE.6) GOTO 140
3308 130 DP(J)=dble(P(I,J))
3309 DBEP=DBEX*DP(1)+DBEY*DP(2)+DBEZ*DP(3)
3310 DGABEP=DGA*(DGA*DBEP/(1D0+DGA)+DP(4))
3311 P(I,1)=sngl(DP(1)+DGABEP*DBEX)
3312 P(I,2)=sngl(DP(2)+DGABEP*DBEY)
3313 P(I,3)=sngl(DP(3)+DGABEP*DBEZ)
3314 P(I,4)=sngl(DGA*(DP(4)+DBEP))
3323 SUBROUTINE HIJHRD(JP,JT,JOUT,JFLG,IOPJT0)
3325 C IOPTJET=1, ALL JET WILL FORM SINGLE STRING SYSTEM
3326 C 0, ONLY Q-QBAR JET FORM SINGLE STRING SYSTEM
3327 C*******Perform jets production and fragmentation when JP JT *******
3328 C scatter. JOUT-> number of hard scatterings precede this one *
3329 C for the the same pair(JP,JT). JFLG->a flag to show whether *
3330 C jets can be produced (with valence quark=1,gluon=2, q-qbar=3)*
3331 C or not(0). Information of jets are in COMMON/ATTJET and *
3332 C /MINJET. ABS(NFP(JP,6)) is the total number jets produced by *
3333 C JP. If NFP(JP,6)<0 JP can not produce jet anymore. *
3334 C*******************************************************************
3335 PARAMETER (MAXSTR=150001)
3336 DIMENSION IP(100,2),IPQ(50),IPB(50),IT(100,2),ITQ(50),ITB(50)
3337 COMMON/hjcrdn/YP(3,300),YT(3,300)
3339 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3341 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
3343 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
3345 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3346 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3347 & PJPM(300,500),NTJ(300),KFTJ(300,500),
3348 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3349 & PJTE(300,500),PJTM(300,500)
3351 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
3352 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
3353 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
3355 c COMMON/HJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
3356 COMMON/HJJET4/NDR,IADR(MAXSTR,2),KFDR(MAXSTR),PDR(MAXSTR,5)
3357 common/xydr/rtdr(MAXSTR,2)
3361 C************************************ HIJING common block
3362 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3364 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3366 COMMON/PYSUBSA/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
3368 COMMON/PYPARSA/MSTP(200),PARP(200),MSTI(200),PARI(200)
3370 COMMON/PYINT1A/MINT(400),VINT(400)
3372 COMMON/PYINT2A/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
3374 COMMON/PYINT5A/NGEN(0:200,3),XSEC(0:200,3)
3376 COMMON/HPINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
3379 C*********************************** LU common block
3381 C SIZE OF COMMON BLOCK FOR # OF PARTON PER STRING
3383 C SIZE OF COMMON BLOCK FOR # OF SINGLE STRINGS
3385 C SIZE OF COMMON BLOCK FOR # OF PARTON PER SINGLE
3392 IF(IOPJET.EQ.1.AND.(NFP(JP,6).NE.0.OR.NFT(JT,6).NE.0))
3394 IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
3395 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) RETURN
3396 C ******** JP or JT can not produce jet anymore
3399 EPP=PP(JP,4)+PP(JP,3)
3400 EPM=PP(JP,4)-PP(JP,3)
3401 ETP=PT(JT,4)+PT(JT,3)
3402 ETM=PT(JT,4)-PT(JT,3)
3403 IF(EPP.LT.0.0) GO TO 1000
3404 IF(EPM.LT.0.0) GO TO 1000
3405 IF(ETP.LT.0.0) GO TO 1000
3406 IF(ETM.LT.0.0) GO TO 1000
3407 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
3409 C ********for the first hard scattering of (JP,JT)
3410 C have collision only when Ycm(JP)>Ycm(JT)
3412 ECUT1=HIPR1(1)+HIPR1(8)+PP(JP,14)+PP(JP,15)
3413 ECUT2=HIPR1(1)+HIPR1(8)+PT(JT,14)+PT(JT,15)
3414 IF(PP(JP,4).LE.ECUT1) THEN
3415 NFP(JP,6)=-ABS(NFP(JP,6))
3418 IF(PT(JT,4).LE.ECUT2) THEN
3419 NFT(JT,6)=-ABS(NFT(JT,6))
3422 C *********must have enough energy to produce jets
3428 IF(NFP(JP,10).EQ.0 .AND. NFT(JT,10).EQ.0) THEN
3436 COEF(11,I)=ATCO(11,I)
3437 COEF(12,I)=ATCO(12,I)
3438 COEF(28,I)=ATCO(28,I)
3444 IF(XSEC(11,1).NE.0) ISUB11=1
3445 IF(XSEC(12,1).NE.0) ISUB12=1
3446 IF(XSEC(28,1).NE.0) ISUB28=1
3447 MINT(44)=MINT4-ISUB11-ISUB12-ISUB28
3448 MINT(45)=MINT5-ISUB11-ISUB12-ISUB28
3449 XSEC(0,1)=ATXS(0)-ATXS(11)-ATXS(12)-ATXS(28)
3459 C ********Scatter the valence quarks only once per NN
3461 C afterwards only gluon can have hard scattering.
3464 IF(JJ.NE.1) GO TO 155
3465 C *********one hard collision at a time
3466 IF(K(7,2).EQ.-K(8,2)) THEN
3467 QMASS2=(P(7,4)+P(8,4))**2-(P(7,1)+P(8,1))**2
3468 & -(P(7,2)+P(8,2))**2-(P(7,3)+P(8,3))**2
3470 IF(QMASS2.LT.(2.0*QM+HIPR1(1))**2) GO TO 155
3472 C ********q-qbar jets must has minimum mass HIPR1(1)
3482 IF(PEP.LE.ECUT1) THEN
3484 IF(MISP.LT.50) GO TO 155
3485 NFP(JP,6)=-ABS(NFP(JP,6))
3488 IF(PET.LE.ECUT2) THEN
3490 IF(MIST.LT.50) GO TO 155
3491 NFT(JT,6)=-ABS(NFT(JT,6))
3494 C ******** if the remain energy<ECUT the proj or targ
3495 C can not produce jet anymore
3499 IF(WP.LT.0.0 .OR. WM.LT.0.0) THEN
3501 clin-6/2009 Let user set the limit when selecting high-Pt events
3502 c because more attempts may be needed:
3503 c IF(MISS.LT.50) GO TO 155
3504 if(pttrig.gt.0) then
3505 if(MISS.LT.maxmiss) then
3506 write(6,*) 'Failed to generate minijet Pt>',pttrig,'GeV'
3510 IF(MISS.LT.50) GO TO 155
3515 C ********the total W+, W- must be positive
3517 AMPX=SQRT((ECUT1-HIPR1(8))**2+PXP**2+PYP**2+0.01)
3518 AMTX=SQRT((ECUT2-HIPR1(8))**2+PXT**2+PYT**2+0.01)
3520 IF(SW.LT.SXX.OR.VINT(43).LT.HIPR1(1)) THEN
3522 clin-6/2009 ctest on
3523 c IF(MISS.LT.50) GO TO 155
3524 IF(MISS.GT.maxmiss) GO TO 155
3527 C ********the proj and targ remnants must have at least
3528 C a CM energy that can produce two strings
3529 C with minimum mass HIPR1(1)(see HIJSFT HIJFRG)
3536 HINT1(46)=SQRT(P(7,1)**2+P(7,2)**2)
3542 HINT1(56)=SQRT(P(8,1)**2+P(8,2)**2)
3546 PINIRD=(1.0-EXP(-2.0*(VINT(47)-HIDAT(1))))
3547 & /(1.0+EXP(-2.0*(VINT(47)-HIDAT(1))))
3549 IF(RANART(NSEED).LE.PINIRD) IINIRD=1
3550 IF(K(7,2).EQ.-K(8,2)) GO TO 190
3551 IF(K(7,2).EQ.21.AND.K(8,2).EQ.21.AND.IOPJET.EQ.1) GO TO 190
3552 C*******************************************************************
3553 C gluon jets are going to be connectd with
3554 C the final leadng string of quark-aintquark
3555 C*******************************************************************
3574 IF(K(I,3).EQ.1 .OR. K(I,3).EQ.2.OR.
3575 & ABS(K(I,2)).GT.30) GO TO 180
3576 C************************************************************
3577 IF(K(I,3).EQ.7) THEN
3578 HINT1(47)=HINT1(47)+P(I,1)
3579 HINT1(48)=HINT1(48)+P(I,2)
3580 HINT1(49)=HINT1(49)+P(I,3)
3581 HINT1(50)=HINT1(50)+P(I,4)
3583 IF(K(I,3).EQ.8) THEN
3584 HINT1(67)=HINT1(67)+P(I,1)
3585 HINT1(68)=HINT1(68)+P(I,2)
3586 HINT1(69)=HINT1(69)+P(I,3)
3587 HINT1(70)=HINT1(70)+P(I,4)
3589 C************************modifcation made on Apr 10. 1996*****
3590 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
3600 rtdr(NDR,1)=0.5*(YP(1,JP)+YT(1,JT))
3601 rtdr(NDR,2)=0.5*(YP(2,JP)+YT(2,JT))
3602 C************************************************************
3604 C************************correction made on Oct. 14,1994*****
3606 IF(K(I,3).EQ.7.OR.K(I,3).EQ.3) THEN
3607 IF(K(I,3).EQ.7.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(7,2)
3608 & .AND.IS7.EQ.0) THEN
3618 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
3619 & IINIRD.EQ.0)) THEN
3629 IF(K(I,2).NE.21) THEN
3630 IF(K(I,2).GT.0) THEN
3634 ELSE IF(K(I,2).LT.0) THEN
3640 ELSE IF(K(I,3).EQ.8.OR.K(I,3).EQ.4) THEN
3641 IF(K(I,3).EQ.8.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(8,2)
3642 & .AND.IS8.EQ.0) THEN
3652 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
3653 & IINIRD.EQ.0)) THEN
3663 IF(K(I,2).NE.21) THEN
3664 IF(K(I,2).GT.0) THEN
3668 ELSE IF(K(I,2).LT.0) THEN
3678 IF(LPQ.NE.LPB .OR. LTQ.NE.LTB) THEN
3680 clin-6/2009 ctest on
3681 c IF(MISS.LE.50) GO TO 155
3682 IF(MISS.LE.maxmiss) GO TO 155
3683 WRITE(6,*) ' Q -QBAR NOT MATCHED IN HIJHRD'
3687 C****The following will rearrange the partons so that a quark is***
3688 C****allways followed by an anti-quark ****************************
3692 IF(J.GT.JPP) GO TO 182
3693 IF(IP(J,2).EQ.0) THEN
3695 ELSE IF(IP(J,2).NE.0) THEN
3699 IP(J,1)=IP(IPQ(LP),1)
3700 IP(J,2)=IP(IPQ(LP),2)
3705 ELSE IF(IP2.LT.0) THEN
3708 C ********replace J with a quark
3711 IP(J+1,1)=IP(IPB(LP),1)
3712 IP(J+1,2)=IP(IPB(LP),2)
3717 ELSE IF(IP2.LT.0) THEN
3720 C ******** replace J+1 with anti-quark
3727 IF(J.GT.JTT) GO TO 184
3728 IF(IT(J,2).EQ.0) THEN
3730 ELSE IF(IT(J,2).NE.0) THEN
3734 IT(J,1)=IT(ITQ(LT),1)
3735 IT(J,2)=IT(ITQ(LT),2)
3740 ELSE IF(IT2.LT.0) THEN
3743 C ********replace J with a quark
3746 IT(J+1,1)=IT(ITB(LT),1)
3747 IT(J+1,2)=IT(ITB(LT),2)
3752 ELSE IF(IT2.LT.0) THEN
3755 C ******** replace J+1 with anti-quark
3762 IF(NPJ(JP)+JPP.GT.MXJT.OR.NTJ(JT)+JTT.GT.MXJT) THEN
3764 WRITE(6,*) 'number of partons per string exceeds'
3765 WRITE(6,*) 'the common block size'
3768 C ********check the bounds of common blocks
3770 KFPJ(JP,NPJ(JP)+J)=K(IP(J,1),2)
3771 PJPX(JP,NPJ(JP)+J)=P(IP(J,1),1)
3772 PJPY(JP,NPJ(JP)+J)=P(IP(J,1),2)
3773 PJPZ(JP,NPJ(JP)+J)=P(IP(J,1),3)
3774 PJPE(JP,NPJ(JP)+J)=P(IP(J,1),4)
3775 PJPM(JP,NPJ(JP)+J)=P(IP(J,1),5)
3779 KFTJ(JT,NTJ(JT)+J)=K(IT(J,1),2)
3780 PJTX(JT,NTJ(JT)+J)=P(IT(J,1),1)
3781 PJTY(JT,NTJ(JT)+J)=P(IT(J,1),2)
3782 PJTZ(JT,NTJ(JT)+J)=P(IT(J,1),3)
3783 PJTE(JT,NTJ(JT)+J)=P(IT(J,1),4)
3784 PJTM(JT,NTJ(JT)+J)=P(IT(J,1),5)
3788 C*****************************************************************
3789 CThis is the case of a quark-antiquark jet it will fragment alone
3790 C****************************************************************
3792 IF(K(7,2).NE.21.AND.K(8,2).NE.21.AND.
3793 & K(7,2)*K(8,2).GT.0) GO TO 155
3798 IF(K(I,3).EQ.1.OR.K(I,3).EQ.2.OR.
3799 & ABS(K(I,2)).GT.30) GO TO 200
3800 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
3810 rtdr(NDR,1)=0.5*(YP(1,JP)+YT(1,JT))
3811 rtdr(NDR,2)=0.5*(YP(2,JP)+YT(2,JT))
3812 C************************************************************
3814 C************************correction made on Oct. 14,1994*****
3816 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
3817 & IINIRD.EQ.0)) THEN
3824 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
3825 & IINIRD.EQ.0)) THEN
3835 IF(K(I,2).NE.21) THEN
3836 IF(K(I,2).GT.0) THEN
3840 ELSE IF(K(I,2).LT.0) THEN
3849 clin-6/2009 ctest on
3850 c IF(MISS.LE.50) GO TO 155
3851 IF(MISS.LE.maxmiss) GO TO 155
3852 WRITE(6,*) LPQ,LPB, 'Q-QBAR NOT CONSERVED OR NOT MATCHED'
3857 C**** The following will rearrange the partons so that a quark is***
3858 C**** allways followed by an anti-quark ****************************
3861 IF(J.GT.JPP) GO TO 222
3862 IF(IP(J,2).EQ.0) GO TO 220
3866 IP(J,1)=IP(IPQ(LP),1)
3867 IP(J,2)=IP(IPQ(LP),2)
3872 ELSE IF(IP2.LT.0) THEN
3876 C ********replace J with a quark
3879 IP(J+1,1)=IP(IPB(LP),1)
3880 IP(J+1,2)=IP(IPB(LP),2)
3885 ELSE IF(IP2.LT.0) THEN
3888 C ******** replace J+1 with an anti-quark
3898 IP(2*L0-3,1)=IP(IPQ(L0),1)
3899 IP(2*L0-3,2)=IP(IPQ(L0),2)
3904 ELSE IF(IP2.LT.0) THEN
3911 IP(2*L0-2,1)=IP(IPB(L0),1)
3912 IP(2*L0-2,2)=IP(IPB(L0),2)
3917 ELSE IF(IP2.LT.0) THEN
3922 C ********move all the qqbar pair to the front of
3923 C the list, except the first pair
3926 IP(2*LPQ-1,1)=IP(IPQ(1),1)
3927 IP(2*LPQ-1,2)=IP(IPQ(1),2)
3932 ELSE IF(IP2.LT.0) THEN
3936 C ********move the first quark to the beginning of
3937 C the last string system
3940 IP(JPP,1)=IP(IPB(1),1)
3941 IP(JPP,2)=IP(IPB(1),2)
3946 ELSE IF(IP2.LT.0) THEN
3950 C ********move the first anti-quark to the end of the
3951 C last string system
3953 IF(NSG.GE.MXSG) THEN
3955 WRITE(6,*) 'number of jets forming single strings exceeds'
3956 WRITE(6,*) 'the common block size'
3959 IF(JPP.GT.MXSJ) THEN
3961 WRITE(6,*) 'number of partons per single jet system'
3962 WRITE(6,*) 'exceeds the common block size'
3965 C ********check the bounds of common block size
3973 K2SG(NSG,I)=K(IP(I,1),2)
3974 IF(K2SG(NSG,I).LT.0) K1SG(NSG,I)=1
3975 PXSG(NSG,I)=P(IP(I,1),1)
3976 PYSG(NSG,I)=P(IP(I,1),2)
3977 PZSG(NSG,I)=P(IP(I,1),3)
3978 PESG(NSG,I)=P(IP(I,1),4)
3979 PMSG(NSG,I)=P(IP(I,1),5)
3983 C******* reset the energy-momentum of incoming particles ********
3995 NFP(JP,6)=NFP(JP,6)+1
3996 NFT(JT,6)=NFT(JT,6)+1
4000 IF(IHPR2(10).EQ.0) RETURN
4001 WRITE(6,*) 'Fatal HIJHRD error'
4002 WRITE(6,*) JP, ' proj E+,E-',EPP,EPM,' status',NFP(JP,5)
4003 WRITE(6,*) JT, ' targ E+,E_',ETP,ETM,' status',NFT(JT,5)
4011 SUBROUTINE JETINI(JP,JT,itrig)
4012 C*******Initialize PYTHIA for jet production**********************
4013 C itrig=0: for normal processes
4014 C itrig=1: for triggered processes
4015 C JP: sequence number of the projectile
4016 C JT: sequence number of the target
4017 C For A+A collisions, one has to initilize pythia
4018 C separately for each type of collisions, pp, pn,np and nn,
4019 C or hp and hn for hA collisions. In this subroutine we use the following
4020 C catalogue for different type of collisions:
4021 C h+h: h+h (itype=1)
4022 C h+A: h+p (itype=1), h+n (itype=2)
4023 C A+h: p+h (itype=1), n+h (itype=2)
4024 C A+A: p+p (itype=1), p+n (itype=2), n+p (itype=3), n+n (itype=4)
4025 C*****************************************************************
4026 CHARACTER BEAM*16,TARG*16
4027 DIMENSION XSEC0(8,0:200),COEF0(8,200,20),INI(8),
4028 & MINT44(8),MINT45(8)
4029 COMMON/hjcrdn/YP(3,300),YT(3,300)
4031 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4033 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
4035 COMMON/HPINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
4038 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
4040 COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
4042 COMMON/PYSUBSA/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
4044 COMMON/PYPARSA/MSTP(200),PARP(200),MSTI(200),PARI(200)
4046 COMMON/PYINT1A/MINT(400),VINT(400)
4048 COMMON/PYINT2A/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
4050 COMMON/PYINT5A/NGEN(0:200,3),XSEC(0:200,3)
4053 clin DATA INI/8*0/ilast/-1/
4054 DATA INI/8*0/,ilast/-1/
4058 IF(IHNT2(5).NE.0 .AND. IHNT2(6).NE.0) THEN
4060 ELSE IF(IHNT2(5).NE.0 .AND. IHNT2(6).EQ.0) THEN
4062 IF(NFT(JT,4).EQ.2112) itype=2
4063 ELSE IF(IHNT2(5).EQ.0 .AND. IHNT2(6).NE.0) THEN
4065 IF(NFP(JP,4).EQ.2112) itype=2
4067 IF(NFP(JP,4).EQ.2212 .AND. NFT(JT,4).EQ.2212) THEN
4069 ELSE IF(NFP(JP,4).EQ.2212 .AND. NFT(JT,4).EQ.2112) THEN
4071 ELSE IF(NFP(JP,4).EQ.2112 .AND. NFT(JT,4).EQ.2212) THEN
4078 IF(itrig.NE.0) GO TO 160
4079 IF(itrig.EQ.ilast) GO TO 150
4081 c ********second order running alpha_strong
4084 C ********inclusion of K factor
4086 C ********Duke-Owens set 1 structure functions
4088 C ********INITIAL STATE RADIATION
4090 C ********FINAL STATE RADIATION
4091 IF(IHPR2(2).EQ.0.OR.IHPR2(2).EQ.2) MSTP(61)=0
4092 IF(IHPR2(2).EQ.0.OR.IHPR2(2).EQ.1) MSTP(71)=0
4095 C ******** NO MULTIPLE INTERACTION
4097 C *******STRUCTURE OF MUTLIPLE INTERACTION
4099 C ********frag off(have to be done by local call)
4100 IF(IHPR2(10).EQ.0) MSTP(122)=0
4101 C ********No printout of initialization information
4106 IF(HIPR1(9).LE.HIPR1(8)) CKIN(4)=-1.0
4121 DO 110 J=1,MIN(8,MDCY(21,3))
4122 110 MDME(MDCY(21,2)+J-1,1)=0
4124 IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
4125 MDME(MDCY(21,2)+ISEL-1,1)=1
4126 C ********QCD subprocesses
4130 C ******* direct photon production
4131 150 IF(INI(itype).NE.0) GO TO 800
4134 C *****triggered subprocesses, jet, photon, heavy quark and DY
4137 IF(itrig.EQ.ilast) GO TO 260
4138 PARP(81)=ABS(HIPR1(10))-0.25
4139 CKIN(5)=ABS(HIPR1(10))-0.25
4140 CKIN(3)=ABS(HIPR1(10))-0.25
4141 CKIN(4)=ABS(HIPR1(10))+0.25
4142 IF(HIPR1(10).LT.HIPR1(8)) CKIN(4)=-1.0
4148 IF(IHPR2(3).EQ.1) THEN
4160 DO 102 J=1,MIN(8,MDCY(21,3))
4161 102 MDME(MDCY(21,2)+J-1,1)=0
4163 IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
4164 MDME(MDCY(21,2)+ISEL-1,1)=1
4165 C ********QCD subprocesses
4166 ELSE IF(IHPR2(3).EQ.2) THEN
4170 C ********Direct photon production
4171 c q+qbar->g+gamma,q+qbar->gamma+gamma, q+g->q+gamma
4172 ELSE IF(IHPR2(3).EQ.3) THEN
4173 CKIN(3)=MAX(0.0,HIPR1(10))
4178 DO 105 J=1,MIN(8,MDCY(21,3))
4179 105 MDME(MDCY(21,2)+J-1,1)=0
4181 IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
4182 MDME(MDCY(21,2)+ISEL-1,1)=1
4183 C **********Heavy quark production
4185 260 IF(INI(itype).NE.0) GO TO 800
4189 IF(IHPR2(10).EQ.0) MSTP(122)=0
4190 IF(NFP(JP,4).EQ.2212) THEN
4192 ELSE IF(NFP(JP,4).EQ.-2212) THEN
4194 ELSE IF(NFP(JP,4).EQ.2112) THEN
4196 ELSE IF(NFP(JP,4).EQ.-2112) THEN
4198 ELSE IF(NFP(JP,4).EQ.211) THEN
4200 ELSE IF(NFP(JP,4).EQ.-211) THEN
4202 ELSE IF(NFP(JP,4).EQ.321) THEN
4204 ELSE IF(NFP(JP,4).EQ.-321) THEN
4207 WRITE(6,*) 'unavailable beam type', NFP(JP,4)
4209 IF(NFT(JT,4).EQ.2212) THEN
4211 ELSE IF(NFT(JT,4).EQ.-2212) THEN
4213 ELSE IF(NFT(JT,4).EQ.2112) THEN
4215 ELSE IF(NFT(JT,4).EQ.-2112) THEN
4217 ELSE IF(NFT(JT,4).EQ.211) THEN
4219 ELSE IF(NFT(JT,4).EQ.-211) THEN
4221 ELSE IF(NFT(JT,4).EQ.321) THEN
4223 ELSE IF(NFT(JT,4).EQ.-321) THEN
4226 WRITE(6,*) 'unavailable target type', NFT(JT,4)
4230 C ******************indicate for initialization use when
4231 C structure functions are called in PYTHIA
4233 CALL PYINITA('CMS',BEAM,TARG,HINT1(1))
4236 MINT44(itype)=MINT(44)
4237 MINT45(itype)=MINT(45)
4239 XSEC0(itype,0)=XSEC(0,1)
4242 XSEC0(itype,I)=XSEC(I,1)
4245 COEF0(itype,I,J)=COEF(I,J)
4251 C ********Store the initialization information for
4255 800 MINT(44)=MINT44(itype)
4256 MINT(45)=MINT45(itype)
4259 XSEC(0,1)=XSEC0(itype,0)
4262 XSEC(I,1)=XSEC0(itype,I)
4265 COEF(I,J)=COEF0(itype,I,J)
4277 PARAMETER (MAXSTR=150001)
4278 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4280 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
4282 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
4283 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
4284 & PJPM(300,500),NTJ(300),KFTJ(300,500),
4285 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
4286 & PJTE(300,500),PJTM(300,500)
4288 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
4289 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
4290 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
4292 c COMMON/HJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
4293 COMMON/HJJET4/NDR,IADR(MAXSTR,2),KFDR(MAXSTR),PDR(MAXSTR,5)
4298 C****************Reset the momentum of initial particles************
4299 C and assign flavors to the proj and targ string *
4300 C*******************************************************************
4305 IF(IHNT2(5).NE.0) IPP=IHNT2(5)
4306 IF(IHNT2(6).NE.0) IPT=IHNT2(6)
4307 C ********in case the proj or targ is a hadron.
4312 PP(I,3)=SQRT(HINT1(1)**2/4.0-HINT1(8)**2)
4335 IF(I.GT.ABS(IHNT2(2))) NFP(I,3)=2112
4336 CALL ATTFLV(NFP(I,3),IDQ,IDQQ)
4340 IF(ABS(IDQ).GT.1000.OR.(ABS(IDQ*IDQQ).LT.100.AND.
4341 & RANART(NSEED).LT.0.5)) NFP(I,15)=1
4342 PP(I,14)=ULMASS(IDQ)
4343 PP(I,15)=ULMASS(IDQQ)
4349 PT(I,3)=-SQRT(HINT1(1)**2/4.0-HINT1(9)**2)
4350 PT(I,4)=HINT1(1)/2.0
4372 IF(I.GT.ABS(IHNT2(4))) NFT(I,3)=2112
4373 CALL ATTFLV(NFT(I,3),IDQ,IDQQ)
4377 IF(ABS(IDQ).GT.1000.OR.(ABS(IDQ*IDQQ).LT.100.AND.
4378 & RANART(NSEED).LT.0.5)) NFT(I,15)=-1
4379 PT(I,14)=ULMASS(IDQ)
4380 PT(I,15)=ULMASS(IDQQ)
4387 SUBROUTINE ATTFLV(ID,IDQ,IDQQ)
4392 IF(ABS(ID).LT.100) THEN
4396 IF(ABS(IDQ).EQ.3) NSIGN=-1
4406 C ********return ID of quark(IDQ) and anti-quark(IDQQ)
4407 C for pions and kaons
4409 C Return LU ID for quarks and diquarks for proton(ID=2212)
4410 C anti-proton(ID=-2212) and nuetron(ID=2112)
4411 C LU ID for d=1,u=2, (ud)0=2101, (ud)1=2103,
4412 C (dd)1=1103,(uu)1=2203.
4413 C Use SU(6) weight proton=1/3d(uu)1 + 1/6u(ud)1 + 1/2u(ud)0
4414 C nurtron=1/3u(dd)1 + 1/6d(ud)1 + 1/2d(ud)0
4417 IF(ABS(ID).EQ.2112) IDQ=1
4420 IF(X.LE.0.5) GO TO 30
4421 IF(X.GT.0.666667) GO TO 10
4426 IF(ABS(ID).EQ.2112) THEN
4438 C*******************************************************************
4439 C This subroutine performs elastic scatterings and possible
4440 C elastic cascading within their own nuclei
4441 c*******************************************************************
4442 SUBROUTINE HIJCSC(JP,JT)
4443 DIMENSION PSC1(5),PSC2(5)
4444 COMMON/hjcrdn/YP(3,300),YT(3,300)
4446 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4450 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
4453 IF(JP.EQ.0 .OR. JT.EQ.0) GO TO 25
4458 CALL HIJELS(PSC1,PSC2)
4459 DPP1=PSC1(1)-PP(JP,1)
4460 DPP2=PSC1(2)-PP(JP,2)
4461 DPT1=PSC2(1)-PT(JT,1)
4462 DPT2=PSC2(2)-PT(JT,2)
4463 PP(JP,6)=PP(JP,6)+DPP1/2.0
4464 PP(JP,7)=PP(JP,7)+DPP2/2.0
4465 PP(JP,8)=PP(JP,8)+DPP1/2.0
4466 PP(JP,9)=PP(JP,9)+DPP2/2.0
4467 PT(JT,6)=PT(JT,6)+DPT1/2.0
4468 PT(JT,7)=PT(JT,7)+DPT2/2.0
4469 PT(JT,8)=PT(JT,8)+DPT1/2.0
4470 PT(JT,9)=PT(JT,9)+DPT2/2.0
4475 NFP(JP,5)=MAX(1,NFP(JP,5))
4476 NFT(JT,5)=MAX(1,NFT(JT,5))
4477 C ********Perform elastic scattering between JP and JT
4479 C ********The following is for possible elastic cascade
4481 25 IF(JP.EQ.0) GO TO 45
4482 PABS=SQRT(PP(JP,1)**2+PP(JP,2)**2+PP(JP,3)**2)
4487 IF(I.EQ.JP) GO TO 40
4491 DIS=DX*BX+DY*BY+DZ*BZ
4492 IF(DIS.LE.0) GO TO 40
4493 BB=DX**2+DY**2+DZ**2-DIS**2
4494 R2=BB*HIPR1(40)/HIPR1(31)/0.1
4495 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
4496 GS=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
4498 GS0=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
4500 IF(RANART(NSEED).GT.GS/GS0) GO TO 40
4505 CALL HIJELS(PSC1,PSC2)
4506 DPP1=PSC1(1)-PP(JP,1)
4507 DPP2=PSC1(2)-PP(JP,2)
4508 DPT1=PSC2(1)-PP(I,1)
4509 DPT2=PSC2(2)-PP(I,2)
4510 PP(JP,6)=PP(JP,6)+DPP1/2.0
4511 PP(JP,7)=PP(JP,7)+DPP2/2.0
4512 PP(JP,8)=PP(JP,8)+DPP1/2.0
4513 PP(JP,9)=PP(JP,9)+DPP2/2.0
4514 PP(I,6)=PP(I,6)+DPT1/2.0
4515 PP(I,7)=PP(I,7)+DPT2/2.0
4516 PP(I,8)=PP(I,8)+DPT1/2.0
4517 PP(I,9)=PP(I,9)+DPT2/2.0
4522 NFP(I,5)=MAX(1,NFP(I,5))
4525 45 IF(JT.EQ.0) GO TO 80
4526 clin 50 PABS=SQRT(PT(JT,1)**2+PT(JT,2)**2+PT(JT,3)**2)
4527 PABS=SQRT(PT(JT,1)**2+PT(JT,2)**2+PT(JT,3)**2)
4532 IF(I.EQ.JT) GO TO 70
4536 DIS=DX*BX+DY*BY+DZ*BZ
4537 IF(DIS.LE.0) GO TO 70
4538 BB=DX**2+DY**2+DZ**2-DIS**2
4539 R2=BB*HIPR1(40)/HIPR1(31)/0.1
4540 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
4541 GS=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
4543 GS0=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
4545 IF(RANART(NSEED).GT.GS/GS0) GO TO 70
4550 CALL HIJELS(PSC1,PSC2)
4551 DPP1=PSC1(1)-PT(JT,1)
4552 DPP2=PSC1(2)-PT(JT,2)
4553 DPT1=PSC2(1)-PT(I,1)
4554 DPT2=PSC2(2)-PT(I,2)
4555 PT(JT,6)=PT(JT,6)+DPP1/2.0
4556 PT(JT,7)=PT(JT,7)+DPP2/2.0
4557 PT(JT,8)=PT(JT,8)+DPP1/2.0
4558 PT(JT,9)=PT(JT,9)+DPP2/2.0
4559 PT(I,6)=PT(I,6)+DPT1/2.0
4560 PT(I,7)=PT(I,7)+DPT2/2.0
4561 PT(I,8)=PT(I,8)+DPT1/2.0
4562 PT(I,9)=PT(I,9)+DPT2/2.0
4567 NFT(I,5)=MAX(1,NFT(I,5))
4574 C*******************************************************************
4575 CThis subroutine performs elastic scattering between two nucleons
4577 C*******************************************************************
4578 SUBROUTINE HIJELS(PSC1,PSC2)
4579 IMPLICIT DOUBLE PRECISION(D)
4580 DIMENSION PSC1(5),PSC2(5)
4581 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4587 CC=1.0-HINT1(12)/HINT1(13)
4588 RR=(1.0-CC)*HINT1(13)/HINT1(12)/(1.0-HIPR1(33))-1.0
4589 BB=0.5*(3.0+RR+SQRT(9.0+10.0*RR+RR**2))
4590 EP=SQRT((PSC1(1)-PSC2(1))**2+(PSC1(2)-PSC2(2))**2
4591 & +(PSC1(3)-PSC2(3))**2)
4592 IF(EP.LE.0.1) RETURN
4593 ELS0=98.0/EP+52.0*(1.0+RR)**2
4594 PCM1=PSC1(1)+PSC2(1)
4595 PCM2=PSC1(2)+PSC2(2)
4596 PCM3=PSC1(3)+PSC2(3)
4600 AMM=ECM**2-PCM1**2-PCM2**2-PCM3**2
4601 IF(AMM.LE.PSC1(5)+PSC2(5)) RETURN
4602 C ********elastic scattering only when approaching
4604 PMAX=(AMM**2+AM1**2+AM2**2-2.0*AMM*AM1-2.0*AMM*AM2
4605 & -2.0*AM1*AM2)/4.0/AMM
4607 20 TT=RANART(NSEED)*MIN(PMAX,1.5)
4608 ELS=98.0*EXP(-2.8*TT)/EP
4609 & +52.0*EXP(-9.2*TT)*(1.0+RR*EXP(-4.6*(BB-1.0)*TT))**2
4610 IF(RANART(NSEED).GT.ELS/ELS0) GO TO 20
4611 PHI=2.0*HIPR1(40)*RANART(NSEED)
4616 DB=dSQRT(DBX**2+DBY**2+DBZ**2)
4617 IF(DB.GT.0.99999999D0) THEN
4618 DBX=DBX*(0.99999999D0/DB)
4619 DBY=DBY*(0.99999999D0/DB)
4620 DBZ=DBZ*(0.99999999D0/DB)
4622 WRITE(6,*) ' (HIJELS) boost vector too large'
4623 C ********Rescale boost vector if too close to unity.
4625 DGA=1D0/SQRT(1D0-DB**2)
4627 DP1=dble(SQRT(TT)*SIN(PHI))
4628 DP2=dble(SQRT(TT)*COS(PHI))
4629 DP3=dble(SQRT(PMAX-TT))
4630 DP4=dble(SQRT(PMAX+AM1))
4631 DBP=DBX*DP1+DBY*DP2+DBZ*DP3
4632 DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4)
4633 PSC1(1)=sngl(DP1+DGABP*DBX)
4634 PSC1(2)=sngl(DP2+DGABP*DBY)
4635 PSC1(3)=sngl(DP3+DGABP*DBZ)
4636 PSC1(4)=sngl(DGA*(DP4+DBP))
4638 DP1=-dble(SQRT(TT)*SIN(PHI))
4639 DP2=-dble(SQRT(TT)*COS(PHI))
4640 DP3=-dble(SQRT(PMAX-TT))
4641 DP4=dble(SQRT(PMAX+AM2))
4642 DBP=DBX*DP1+DBY*DP2+DBZ*DP3
4643 DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4)
4644 PSC2(1)=sngl(DP1+DGABP*DBX)
4645 PSC2(2)=sngl(DP2+DGABP*DBY)
4646 PSC2(3)=sngl(DP3+DGABP*DBZ)
4647 PSC2(4)=sngl(DGA*(DP4+DBP))
4652 C*******************************************************************
4654 C Subroutine HIJSFT *
4656 C Scatter two excited strings, JP from proj and JT from target *
4657 C*******************************************************************
4658 SUBROUTINE HIJSFT(JP,JT,JOUT,IERROR)
4659 PARAMETER (MAXSTR=150001)
4660 COMMON/hjcrdn/YP(3,300),YT(3,300)
4662 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4664 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
4668 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
4669 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
4670 & PJPM(300,500),NTJ(300),KFTJ(300,500),
4671 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
4672 & PJTE(300,500),PJTM(300,500)
4675 c COMMON/HJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
4676 c & K2SG(900,100),PXSG(900,100),PYSG(900,100),
4677 c & PZSG(900,100),PESG(900,100),PMSG(900,100)
4679 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
4681 COMMON/DPMCM1/JJP,JJT,AMP,AMT,APX0,ATX0,AMPN,AMTN,AMP0,AMT0,
4682 & NFDP,NFDT,WP,WM,SW,XREMP,XREMT,DPKC1,DPKC2,PP11,PP12,
4683 & PT11,PT12,PTP2,PTT2
4685 COMMON/DPMCM2/NDPM,KDPM(20,2),PDPM1(20,5),PDPM2(20,5)
4688 C*******************************************************************
4690 C of hard scatterings preceding this soft collision.
4692 C double diffrac 2=single diffrac, 3=non-single diffrac.
4693 C*******************************************************************
4699 IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
4701 EPP=PP(JP,4)+PP(JP,3)
4702 EPM=PP(JP,4)-PP(JP,3)
4703 ETP=PT(JT,4)+PT(JT,3)
4704 ETM=PT(JT,4)-PT(JT,3)
4709 C ********total W+,W- and center-of-mass energy
4711 IF(WP.LT.0.0 .OR. WM.LT.0.0) GO TO 1000
4714 IF(EPP.LT.0.0) GO TO 1000
4715 IF(EPM.LT.0.0) GO TO 1000
4716 IF(ETP.LT.0.0) GO TO 1000
4717 IF(ETM.LT.0.0) GO TO 1000
4718 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
4720 C ********For strings which does not follow a jet-prod,
4721 C scatter only if Ycm(JP)>Ycm(JT). When jets
4722 C are produced just before this collision
4723 C this requirement has already be enforced
4724 C (see SUBROUTINE HIJHRD)
4741 IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
4742 IF(NFP(JP,10).EQ.1) THEN
4743 PHI1=ULANGL(PP(JP,10),PP(JP,11))
4744 PPJET=SQRT(PP(JP,10)**2+PP(JP,11)**2)
4749 IF(NFT(JT,10).EQ.1) THEN
4750 PHI2=ULANGL(PT(JT,10),PT(JT,11))
4751 PTJET=SQRT(PT(JT,10)**2+PT(JT,11)**2)
4756 IF(IHPR2(4).GT.0.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
4757 IF(NFP(JP,10).EQ.0) THEN
4759 ELSE IF(NFT(JT,10).EQ.0) THEN
4762 PHI=(PHI1+PHI2-HIPR1(40))/2.0
4764 BX=HINT1(19)*COS(HINT1(20))
4765 BY=HINT1(19)*SIN(HINT1(20))
4770 R1=MAX(1.2*IHNT2(1)**0.3333333,
4771 & SQRT(XP0**2+YP0**2))
4772 R2=MAX(1.2*IHNT2(3)**0.3333333,
4773 & SQRT((XT0-BX)**2+(YT0-BY)**2))
4774 IF(ABS(COS(PHI)).LT.1.0E-5) THEN
4777 DD3=ABS(BY+SQRT(R2**2-(XP0-BX)**2)-YP0)
4778 DD4=ABS(BY-SQRT(R2**2-(XP0-BX)**2)-YP0)
4781 BB=2.0*SIN(PHI)*(COS(PHI)*YP0-SIN(PHI)*XP0)
4782 CC=(YP0**2-R1**2)*COS(PHI)**2+XP0*SIN(PHI)*(
4783 & XP0*SIN(PHI)-2.0*YP0*COS(PHI))
4785 IF(DD.LT.0.0) GO TO 10
4786 XX1=(-BB+SQRT(DD))/2.0
4787 XX2=(-BB-SQRT(DD))/2.0
4788 DD1=ABS((XX1-XP0)/COS(PHI))
4789 DD2=ABS((XX2-XP0)/COS(PHI))
4791 BB=2.0*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)-2.0*BX
4792 CC=(BX**2+(YT0-BY)**2-R2**2)*COS(PHI)**2+XT0*SIN(PHI)
4793 & *(XT0*SIN(PHI)-2.0*COS(PHI)*(YT0-BY))
4794 & -2.0*BX*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)
4796 IF(DD.LT.0.0) GO TO 10
4797 XX1=(-BB+SQRT(DD))/2.0
4798 XX2=(-BB-SQRT(DD))/2.0
4799 DD3=ABS((XX1-XT0)/COS(PHI))
4800 DD4=ABS((XX2-XT0)/COS(PHI))
4804 IF(DD1.LT.HIPR1(13)) DD1=0.0
4805 IF(DD2.LT.HIPR1(13)) DD2=0.0
4806 IF(NFP(JP,10).EQ.1.AND.PPJET.GT.HIPR1(11)) THEN
4807 DP1=DD1*HIPR1(14)/2.0
4808 DP1=MIN(DP1,PPJET-HIPR1(11))
4812 PKC11=PP(JP,10)-DPX1
4813 PKC12=PP(JP,11)-DPY1
4815 CTHEP=PP(JP,12)/SQRT(PP(JP,12)**2+PPJET**2)
4816 DPZ1=DP1*CTHEP/SQRT(1.0-CTHEP**2)
4817 DPE1=SQRT(DPX1**2+DPY1**2+DPZ1**2)
4818 EPPPRM=PP(JP,4)+PP(JP,3)-DPE1-DPZ1
4819 EPMPRM=PP(JP,4)-PP(JP,3)-DPE1+DPZ1
4820 IF(EPPPRM.LE.0.0.OR.EPMPRM.LE.0.0) GO TO 15
4827 PJPX(JP,NPJ(JP))=DPX1
4828 PJPY(JP,NPJ(JP))=DPY1
4829 PJPZ(JP,NPJ(JP))=DPZ1
4830 PJPE(JP,NPJ(JP))=DPE1
4831 PJPM(JP,NPJ(JP))=0.0
4832 PP(JP,3)=PP(JP,3)-DPZ1
4833 PP(JP,4)=PP(JP,4)-DPE1
4836 15 IF(NFT(JT,10).EQ.1.AND.PTJET.GT.HIPR1(11)) THEN
4837 DP2=DD2*HIPR1(14)/2.0
4838 DP2=MIN(DP2,PTJET-HIPR1(11))
4842 PKC21=PT(JT,10)-DPX2
4843 PKC22=PT(JT,11)-DPY2
4845 CTHET=PT(JT,12)/SQRT(PT(JT,12)**2+PTJET**2)
4846 DPZ2=DP2*CTHET/SQRT(1.0-CTHET**2)
4847 DPE2=SQRT(DPX2**2+DPY2**2+DPZ2**2)
4848 ETPPRM=PT(JT,4)+PT(JT,3)-DPE2-DPZ2
4849 ETMPRM=PT(JT,4)-PT(JT,3)-DPE2+DPZ2
4850 IF(ETPPRM.LE.0.0.OR.ETMPRM.LE.0.0) GO TO 16
4857 PJTX(JT,NTJ(JT))=DPX2
4858 PJTY(JT,NTJ(JT))=DPY2
4859 PJTZ(JT,NTJ(JT))=DPZ2
4860 PJTE(JT,NTJ(JT))=DPE2
4861 PJTM(JT,NTJ(JT))=0.0
4862 PT(JT,3)=PT(JT,3)-DPZ2
4863 PT(JT,4)=PT(JT,4)-DPE2
4866 16 DPKC11=-(PP(JP,10)-PKC11)/2.0
4867 DPKC12=-(PP(JP,11)-PKC12)/2.0
4868 DPKC21=-(PT(JT,10)-PKC21)/2.0
4869 DPKC22=-(PT(JT,11)-PKC22)/2.0
4875 C ********If jet is quenched the pt from valence quark
4876 C hard scattering has to reduced by d*kapa
4879 10 PTP02=PP(JP,1)**2+PP(JP,2)**2
4880 PTT02=PT(JT,1)**2+PT(JT,2)**2
4882 AMQ=MAX(PP(JP,14)+PP(JP,15),PT(JT,14)+PT(JT,15))
4884 C ********consider mass cut-off for strings which
4885 C must also include quark's mass
4889 IF(NFP(JP,5).LE.2.AND.NFP(JP,3).NE.0) THEN
4890 AMP0=ULMASS(NFP(JP,3))
4891 NFDP=NFP(JP,3)+2*NFP(JP,3)/ABS(NFP(JP,3))
4893 IF(DPM0.LE.0.0) THEN
4894 NFDP=NFDP-2*NFDP/ABS(NFDP)
4901 IF(NFT(JT,5).LE.2.AND.NFT(JT,3).NE.0) THEN
4902 AMT0=ULMASS(NFT(JT,3))
4903 NFDT=NFT(JT,3)+2*NFT(JT,3)/ABS(NFT(JT,3))
4905 IF(DTM0.LE.0.0) THEN
4906 NFDT=NFDT-2*NFDT/ABS(NFDT)
4911 AMPN=SQRT(AMP0**2+PTP02)
4912 AMTN=SQRT(AMT0**2+PTT02)
4913 SNN=(AMPN+AMTN)**2+0.001
4915 IF(SW.LT.SNN+0.001) GO TO 4000
4916 C ********Scatter only if SW>SNN
4917 C*****give some PT kick to the two exited strings******************
4918 clin 20 SWPTN=4.0*(MAX(AMP0,AMT0)**2+MAX(PTP02,PTT02))
4919 SWPTN=4.0*(MAX(AMP0,AMT0)**2+MAX(PTP02,PTT02))
4920 SWPTD=4.0*(MAX(DPM0,DTM0)**2+MAX(PTP02,PTT02))
4921 SWPTX=4.0*(AMX**2+MAX(PTP02,PTT02))
4922 IF(SW.LE.SWPTN) THEN
4924 ELSE IF(SW.GT.SWPTN .AND. SW.LE.SWPTD
4925 & .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
4926 PKCMX=SQRT(SW/4.0-MAX(AMP0,AMT0)**2)
4927 & -SQRT(MAX(PTP02,PTT02))
4928 ELSE IF(SW.GT.SWPTD .AND. SW.LE.SWPTX
4929 & .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
4930 PKCMX=SQRT(SW/4.0-MAX(DPM0,DTM0)**2)
4931 & -SQRT(MAX(PTP02,PTT02))
4932 ELSE IF(SW.GT.SWPTX) THEN
4933 PKCMX=SQRT(SW/4.0-AMX**2)-SQRT(MAX(PTP02,PTT02))
4935 C ********maximun PT kick
4936 C*********************************************************
4938 IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
4939 IF(PKC1.GT.PKCMX) THEN
4941 PKC11=PKC1*COS(PHI1)
4942 PKC12=PKC1*SIN(PHI1)
4943 DPKC11=-(PP(JP,10)-PKC11)/2.0
4944 DPKC12=-(PP(JP,11)-PKC12)/2.0
4946 IF(PKC2.GT.PKCMX) THEN
4948 PKC21=PKC2*COS(PHI2)
4949 PKC22=PKC2*SIN(PHI2)
4950 DPKC21=-(PT(JT,10)-PKC21)/2.0
4951 DPKC22=-(PT(JT,11)-PKC22)/2.0
4955 NFP(JP,10)=-NFP(JP,10)
4956 NFT(JT,10)=-NFT(JT,10)
4959 C ********If the valence quarks had a hard-collision
4960 C the pt kick is the pt from hard-collision.
4962 IF(IHPR2(13).NE.0 .AND. RANART(NSEED).LE.HIDAT(4)) isng=1
4963 IF((NFP(JP,5).EQ.3 .OR.NFT(JT,5).EQ.3).OR.
4964 & (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
4965 & (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) isng=0
4967 C ********decite whether to have single-diffractive
4968 IF(IHPR2(5).EQ.0) THEN
4969 PKC=HIPR1(2)*SQRT(-ALOG(1.0-RANART(NSEED)
4970 & *(1.0-EXP(-PKCMX**2/HIPR1(2)**2))))
4974 clin-10/28/02 get rid of argument usage mismatch in HIRND2():
4975 c PKC=HIRND2(3,0.0,PKCMX**2)
4978 PKC=HIRND2(3,xminhi,xmaxhi)
4981 IF(PKC.GT.HIPR1(20))
4982 & PKC=HIPR1(2)*SQRT(-ALOG(EXP(-HIPR1(20)**2/HIPR1(2)**2)
4983 & -RANART(NSEED)*(EXP(-HIPR1(20)**2/HIPR1(2)**2)-
4984 & EXP(-PKCMX**2/HIPR1(2)**2))))
4986 IF(isng.EQ.1) PKC=0.65*SQRT(
4987 & -ALOG(1.0-RANART(NSEED)*(1.0-EXP(-PKCMX**2/0.65**2))))
4988 C ********select PT kick
4989 30 PHI0=2.0*HIPR1(40)*RANART(NSEED)
4996 40 PP11=PP(JP,1)+PKC11-DPKC1
4997 PP12=PP(JP,2)+PKC12-DPKC2
4998 PT11=PT(JT,1)+PKC21-DPKC1
4999 PT12=PT(JT,2)+PKC22-DPKC2
5000 PTP2=PP11**2+PP12**2
5001 PTT2=PT11**2+PT12**2
5003 AMPN=SQRT(AMP0**2+PTP2)
5004 AMTN=SQRT(AMT0**2+PTT2)
5005 SNN=(AMPN+AMTN)**2+0.001
5006 C***************************************
5010 C****************************************
5013 IF(MISS.LE.100) then
5018 & WRITE(6,*) 'Error occured in Pt kick section of HIJSFT'
5021 C******************************************************************
5022 AMPD=SQRT(DPM0**2+PTP2)
5023 AMTD=SQRT(DTM0**2+PTT2)
5025 AMPX=SQRT(AMX**2+PTP2)
5026 AMTX=SQRT(AMX**2+PTT2)
5035 SPNTD=(AMPN+AMTD)**2
5036 SPNTX=(AMPN+AMTX)**2
5037 C ********CM energy if proj=N,targ=N*
5038 SPDTN=(AMPD+AMTN)**2
5039 SPXTN=(AMPX+AMTN)**2
5040 C ********CM energy if proj=N*,targ=N
5041 SPDTX=(AMPD+AMTX)**2
5042 SPXTD=(AMPX+AMTD)**2
5048 C ********CM energy if proj=delta, targ=delta
5049 C****************There are many different cases**********
5050 c IF(IHPR2(15).EQ.1) GO TO 500
5052 C ********to have DPM type soft interactions
5055 IF(SW.GT.SXX+0.001) THEN
5063 c**** 5/30/1998 this is identical to the above statement. Added to
5064 c**** avoid questional branching to block.
5065 IF((NFP(JP,5).EQ.3 .AND.NFT(JT,5).EQ.3).OR.
5066 & (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
5067 & (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
5074 C ********do not allow excited strings to have
5076 IF(RANART(NSEED).GT.0.5.OR.(NFT(JT,5).GT.2.OR.
5077 & NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
5090 C ********have single diffractive collision
5092 ELSE IF(SW.GT.MAX(SPDTX,SPXTD)+0.001 .AND.
5093 & SW.LE.SXX+0.001) THEN
5094 IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0.AND.
5095 & RANART(NSEED).GT.0.5).OR.(NPJ(JP).EQ.0
5096 & .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
5102 ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
5110 ELSE IF(SW.GT.MIN(SPDTX,SPXTD)+0.001.AND.
5111 & SW.LE.MAX(SPDTX,SPXTD)+0.001) THEN
5112 IF(SPDTX.LE.SPXTD.AND.NPJ(JP).EQ.0
5113 & .AND.NFP(JP,5).LE.2) THEN
5119 ELSE IF(SPDTX.GT.SPXTD.AND.NTJ(JT).EQ.0
5120 & .AND.NFT(JT,5).LE.2) THEN
5127 c*** 5/30/1998 added to avoid questional branching to another block
5128 c*** this is identical to the statement following the next ELSE IF
5129 IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
5130 & .AND.RANART(NSEED).GT.0.5).OR.(NPJ(JP).EQ.0
5131 & .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
5137 ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
5145 ELSE IF(SW.GT.MAX(SPNTX,SPXTN)+0.001 .AND.
5146 & SW.LE.MIN(SPDTX,SPXTD)+0.001) THEN
5147 IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
5148 & .AND.RANART(NSEED).GT.0.5).OR.(NPJ(JP).EQ.0
5149 & .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
5155 ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
5163 ELSE IF(SW.GT.MIN(SPNTX,SPXTN)+0.001 .AND.
5164 & SW.LE.MAX(SPNTX,SPXTN)+0.001) THEN
5165 IF(SPNTX.LE.SPXTN.AND.NPJ(JP).EQ.0
5166 & .AND.NFP(JP,5).LE.2) THEN
5172 ELSEIF(SPNTX.GT.SPXTN.AND.NTJ(JT).EQ.0
5173 & .AND.NFT(JT,5).LE.2) THEN
5181 ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
5182 & (NPJ(JP).NE.0 .OR.NTJ(JT).NE.0)) THEN
5184 ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
5185 & NFP(JP,5).GT.2.AND.NFT(JT,5).GT.2) THEN
5187 ELSE IF(SW.GT.SDD+0.001.AND.SW.LE.
5188 & MIN(SPNTX,SPXTN)+0.001) THEN
5194 ELSE IF(SW.GT.MAX(SPNTD,SPDTN)+0.001
5195 & .AND. SW.LE.SDD+0.001) THEN
5196 IF(RANART(NSEED).GT.0.5) THEN
5209 ELSE IF(SW.GT.MIN(SPNTD,SPDTN)+0.001
5210 & .AND. SW.LE.MAX(SPNTD,SPDTN)+0.001) THEN
5211 IF(SPNTD.GT.SPDTN) THEN
5224 ELSE IF(SW.LE.MIN(SPNTD,SPDTN)+0.001) THEN
5231 WRITE(6,*) ' Error in HIJSFT: There is no path to here'
5234 C*************** elastic scattering ***************
5235 C this is like elastic, both proj and targ mass
5237 C***************************************************
5238 100 NFP5=MAX(2,NFP(JP,5))
5239 NFT5=MAX(2,NFT(JT,5))
5242 IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
5244 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
5248 IF(RANART(NSEED).LT.0.5) THEN
5249 X1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
5250 X2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
5252 X1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
5253 X2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
5258 C********** Single diffractive ***********************
5259 C either proj or targ's mass is fixed
5260 C*****************************************************
5261 220 NFP5=MAX(2,NFP(JP,5))
5263 IF(NFP3.EQ.0) NFP5=3
5265 IF(BB2**2.LT.4.0*D2) THEN
5267 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
5271 XMIN=(BB2-SQRT(BB2**2-4.0*D2))/2.0
5272 XMAX=(BB2+SQRT(BB2**2-4.0*D2))/2.0
5274 222 X2=HIRND2(6,XMIN,XMAX)
5276 IF(X2*(1.0-X1).LT.(D2+1.E-4/SW)) THEN
5278 IF(MISS4.LE.1000) GO TO 222
5283 C ********Fix proj mass*********
5285 NFT5=MAX(2,NFT(JT,5))
5286 IF(NFT3.EQ.0) NFT5=3
5288 IF(BB1**2.LT.4.0*D1) THEN
5290 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
5294 XMIN=(BB1-SQRT(BB1**2-4.0*D1))/2.0
5295 XMAX=(BB1+SQRT(BB1**2-4.0*D1))/2.0
5297 242 X1=HIRND2(6,XMIN,XMAX)
5299 IF(X1*(1.0-X2).LT.(D1+1.E-4/SW)) THEN
5301 IF(MISS4.LE.1000) GO TO 242
5306 C ********Fix targ mass*********
5308 C*************non-single diffractive**********************
5309 C both proj and targ may not be fixed in mass
5310 C*********************************************************
5316 IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
5318 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
5322 XMIN1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
5323 XMAX1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
5324 XMIN2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
5325 XMAX2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
5327 410 X1=HIRND2(4,XMIN1,XMAX1)
5328 X2=HIRND2(4,XMIN2,XMAX2)
5329 IF(NFP(JP,5).EQ.3.OR.NFT(JT,5).EQ.3) THEN
5330 X1=HIRND2(6,XMIN1,XMAX1)
5331 X2=HIRND2(6,XMIN2,XMAX2)
5334 IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
5335 & ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
5336 X1=HIRND2(5,XMIN1,XMAX1)
5338 IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
5339 & ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
5340 X2=HIRND2(5,XMIN2,XMAX2)
5342 c IF(IOPMAIN.EQ.3) X1=HIRND2(6,XMIN1,XMAX1)
5343 c IF(IOPMAIN.EQ.2) X2=HIRND2(6,XMIN2,XMAX2)
5344 C ********For q-qbar or (qq)-(qq)bar system use symetric
5345 C distribution, for q-(qq) or qbar-(qq)bar use
5346 C unsymetrical distribution
5348 IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000) X1=1.0-X1
5351 IF(XXP.LT.(D1+1.E-4/SW) .OR. XXT.LT.(D2+1.E-4/SW)) THEN
5353 IF(MISS4.LE.1000) GO TO 410
5357 C***************************************************
5358 C***************************************************
5360 IF(X1*(1.0-X2).LT.(AMPN**2-1.E-4)/SW.OR.
5361 & X2*(1.0-X1).LT.(AMTN**2-1.E-4)/SW) THEN
5363 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 2000
5372 PP(JP,3)=(EPP-EPM)/2.0
5373 PP(JP,4)=(EPP+EPM)/2.0
5374 IF(EPP*EPM-PTP2.LT.0.0) GO TO 6000
5375 PP(JP,5)=SQRT(EPP*EPM-PTP2)
5379 PT(JT,3)=(ETP-ETM)/2.0
5380 PT(JT,4)=(ETP+ETM)/2.0
5381 IF(ETP*ETM-PTT2.LT.0.0) GO TO 6000
5382 PT(JT,5)=SQRT(ETP*ETM-PTT2)
5385 C*****recoil PT from hard-inter is shared by two end-partons
5392 IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
5393 & ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
5396 IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
5397 & ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
5400 IF((KCDIP.EQ.0.AND.RANART(NSEED).LT.0.5)
5401 & .OR.(KCDIP.NE.0.AND.RANART(NSEED)
5402 & .LT.0.5/(1.0+(PKC11**2+PKC12**2)/HIPR1(22)**2))) THEN
5403 PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,6)
5404 PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,7)
5405 PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
5407 PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
5410 PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,8)
5411 PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,9)
5412 PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
5414 PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
5417 PP(JP,1)=PP(JP,6)+PP(JP,8)
5418 PP(JP,2)=PP(JP,7)+PP(JP,9)
5419 C ********pt kick for proj
5422 IF((KCDIT.EQ.0.AND.RANART(NSEED).LT.0.5)
5423 & .OR.(KCDIT.NE.0.AND.RANART(NSEED)
5424 & .LT.0.5/(1.0+(PKC21**2+PKC22**2)/HIPR1(22)**2))) THEN
5425 PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,6)
5426 PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,7)
5427 PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
5429 PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
5432 PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,8)
5433 PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,9)
5434 PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
5436 PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
5439 PT(JT,1)=PT(JT,6)+PT(JT,8)
5440 PT(JT,2)=PT(JT,7)+PT(JT,9)
5441 C ********pt kick for targ
5443 IF(NPJ(JP).NE.0) NFP(JP,5)=3
5444 IF(NTJ(JT).NE.0) NFT(JT,5)=3
5445 C ********jets must be connected to string
5446 IF(EPP/(EPM+0.0001).LT.ETP/(ETM+0.0001).AND.
5447 & ABS(NFP(JP,1)*NFP(JP,2)).LT.1000000)THEN
5450 PP(JP,JSB)=PT(JT,JSB)
5453 NFP(JP,JSB)=NFT(JT,JSB)
5456 C ********when Ycm(JP)<Ycm(JT) after the collision
5457 C exchange the positions of the two
5461 C**************************************************
5462 C**************************************************
5464 IF(IHPR2(10).EQ.0) RETURN
5465 WRITE(6,*) ' Fatal HIJSFT start error,abandon this event'
5466 WRITE(6,*) ' PROJ E+,E-,W+',EPP,EPM,WP
5467 WRITE(6,*) ' TARG E+,E-,W-',ETP,ETM,WM
5468 WRITE(6,*) ' W+*W-, (APN+ATN)^2',SW,SNN
5471 IF(IHPR2(10).EQ.0) RETURN
5472 WRITE(6,*) ' (2)energy partition fail,'
5473 WRITE(6,*) ' HIJSFT not performed, but continue'
5474 WRITE(6,*) ' MP1,MPN',X1*(1.0-X2)*SW,AMPN**2
5475 WRITE(6,*) ' MT2,MTN',X2*(1.0-X1)*SW,AMTN**2
5478 IF(IHPR2(10).EQ.0) RETURN
5479 WRITE(6,*) ' (3)something is wrong with the pt kick, '
5480 WRITE(6,*) ' HIJSFT not performed, but continue'
5481 WRITE(6,*) ' D1=',D1,' D2=',D2,' SW=',SW
5482 WRITE(6,*) ' HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
5483 WRITE(6,*) ' THIS COLLISON NFP5=',NFP5, ' NFT5=',NFT5
5484 WRITE(6,*) ' # OF JET IN PROJ',NPJ(JP),' IN TARG',NTJ(JT)
5487 IF(IHPR2(10).EQ.0) RETURN
5488 WRITE(6,*) ' (4)unable to choose process, but not harmful'
5489 WRITE(6,*) ' HIJSFT not performed, but continue'
5490 WRITE(6,*) ' PTP=',SQRT(PTP2),' PTT=',SQRT(PTT2),' SW=',SW
5491 WRITE(6,*) ' AMCUT=',AMX,' JP=',JP,' JT=',JT
5492 WRITE(6,*) ' HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
5495 IF(IHPR2(10).EQ.0) RETURN
5496 WRITE(6,*) ' energy partition failed(5),for limited try'
5497 WRITE(6,*) ' HIJSFT not performed, but continue'
5498 WRITE(6,*) ' NFP5=',NFP5,' NFT5=',NFT5
5499 WRITE(6,*) ' D1',D1,' X1(1-X2)',X1*(1.0-X2)
5500 WRITE(6,*) ' D2',D2,' X2(1-X1)',X2*(1.0-X1)
5504 IF(MISS.LT.100) GO TO 30
5506 IF(IHPR2(10).EQ.0) RETURN
5507 WRITE(6,*) ' ERROR OCCURED, HIJSFT NOT PERFORMED'
5508 WRITE(6,*) ' Abort this event'
5509 WRITE(6,*) 'MTP,PTP2',EPP*EPM,PTP2,' MTT,PTT2',ETP*ETM,PTT2
5515 C ********************************************************
5516 C ************************ WOOD-SAX
5517 SUBROUTINE HIJWDS(IA,IDH,XHIGH)
5518 C SETS UP HISTOGRAM IDH WITH RADII FOR
5519 C NUCLEUS IA DISTRIBUTED ACCORDING TO THREE PARAM WOOD SAXON
5520 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5522 COMMON/WOOD/R,D,FNORM,W
5524 c DIMENSION IAA(20),RR(20),DD(20),WW(20),RMS(20)
5525 DIMENSION IAA(20),RR(20),DD(20),WW(20)
5526 EXTERNAL RWDSAX,WDSAX
5529 C PARAMETERS OF SPECIAL NUCLEI FROM ATOMIC DATA AND NUC DATA TABLES
5531 DATA IAA/2,4,12,16,27,32,40,56,63,93,184,197,208,7*0./
5532 DATA RR/0.01,.964,2.355,2.608,2.84,3.458,3.766,3.971,4.214,
5533 1 4.87,6.51,6.38,6.624,7*0./
5534 DATA DD/0.5882,.322,.522,.513,.569,.61,.586,.5935,.586,.573,
5535 1 .535,.535,.549,7*0./
5536 DATA WW/0.0,.517,-0.149,-0.051,0.,-0.208,-0.161,13*0./
5537 c DATA RMS/2.11,1.71,2.46,2.73,3.05,3.247,3.482,3.737,3.925,4.31,
5538 c 1 5.42,5.33,5.521,7*0./
5542 C ********SET WOOD-SAX PARAMS FIRST AS IN DATE ET AL
5544 C ********D IS WOOD SAX DIFFUSE PARAM IN FM
5545 R=1.19*A**(1./3.) - 1.61*A**(-1./3.)
5546 C ********R IS RADIUS PARAM
5548 C ********W IS The third of three WOOD-SAX PARAM
5550 C ********CHECK TABLE FOR SPECIAL CASES
5552 IF (IA.EQ.IAA(I)) THEN
5556 clin RS not used RS=RMS(I)
5559 C ********FNORM is the normalize factor
5563 IF (W.LT.-0.01) THEN
5564 IF (XHIGH.GT.R/SQRT(ABS(W))) XHIGH=R/SQRT(ABS(W))
5566 FGAUS=GAUSS1(RWDSAX,XLOW,XHIGH,0.001)
5573 HINT1(75)=FNORM/4.0/HIPR1(40)
5574 ELSE IF (IDH.EQ.2) THEN
5578 HINT1(79)=FNORM/4.0/HIPR1(40)
5581 C NOW SET UP HBOOK FUNCTIONS IDH FOR R**2*RHO(R)
5582 C THESE HISTOGRAMS ARE USED TO GENERATE RANDOM RADII
5583 CALL HIFUN(IDH,XLOW,XHIGH,RWDSAX)
5589 C ********THREE PARAMETER WOOD SAXON
5590 COMMON/WOOD/R,D,FNORM,W
5593 WDSAX=FNORM*(1.+W*(X/R)**2)/(1+EXP((X-R)/D))
5595 IF (X.GE.R/SQRT(ABS(W))) WDSAX=0.
5610 C The next three subroutines are for Monte Carlo generation
5611 C according to a given function FHB. One calls first HIFUN
5612 C with assigned channel number I, low and up limits. Then to
5613 C generate the distribution one can call HIRND(I) which gives
5614 C you a random number generated according to the given function.
5616 SUBROUTINE HIFUN(I,XMIN,XMAX,FHB)
5617 COMMON/HIJHB/RR(10,201),XX(10,201)
5621 FNORM=GAUSS1(FHB,XMIN,XMAX,0.001)
5623 XX(I,J)=XMIN+(XMAX-XMIN)*(J-1)/200.0
5625 RR(I,J)=GAUSS1(FHB,XMIN,XDD,0.001)/FNORM
5633 COMMON/HIJHB/RR(10,201),XX(10,201)
5641 10 IF(JU-JL.GT.1) THEN
5643 IF((RR(I,201).GT.RR(I,1)).EQV.(RX.GT.RR(I,JM))) THEN
5653 HIRND=(XX(I,J)+XX(I,J+1))/2.0
5660 C This generate random number between XMIN and XMAX
5661 FUNCTION HIRND2(I,XMIN,XMAX)
5662 COMMON/HIJHB/RR(10,201),XX(10,201)
5667 IF(XMIN.LT.XX(I,1)) XMIN=XX(I,1)
5668 IF(XMAX.GT.XX(I,201)) XMAX=XX(I,201)
5669 JMIN=1+int(200*(XMIN-XX(I,1))/(XX(I,201)-XX(I,1)))
5670 JMAX=1+int(200*(XMAX-XX(I,1))/(XX(I,201)-XX(I,1)))
5671 RX=RR(I,JMIN)+(RR(I,JMAX)-RR(I,JMIN))*RANART(NSEED)
5674 10 IF(JU-JL.GT.1) THEN
5676 IF((RR(I,201).GT.RR(I,1)).EQV.(RX.GT.RR(I,JM))) THEN
5686 HIRND2=(XX(I,J)+XX(I,J+1))/2.0
5694 C THIS IS TO CALCULATE THE CROSS SECTIONS OF JET PRODUCTION AND
5695 C THE TOTAL INELASTIC CROSS SECTIONS.
5696 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5700 EXTERNAL FHIN,FTOT,FNJET,FTOTJT,FTOTRG
5702 IF(HINT1(1).GE.10.0) CALL CRSJET
5703 C ********calculate jet cross section(in mb)
5705 APHX1=HIPR1(6)*(IHNT2(1)**0.3333333-1.0)
5706 APHX2=HIPR1(6)*(IHNT2(3)**0.3333333-1.0)
5707 HINT1(11)=HINT1(14)-APHX1*HINT1(15)
5708 & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
5709 HINT1(10)=GAUSS1(FTOTJT,0.0,20.0,0.01)
5710 HINT1(12)=GAUSS1(FHIN,0.0,20.0,0.01)
5711 HINT1(13)=GAUSS1(FTOT,0.0,20.0,0.01)
5712 HINT1(60)=HINT1(61)-APHX1*HINT1(62)
5713 & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
5714 HINT1(59)=GAUSS1(FTOTRG,0.0,20.0,0.01)
5715 IF(HINT1(59).EQ.0.0) HINT1(59)=HINT1(60)
5716 IF(HINT1(1).GE.10.0) Then
5719 HINT1(80+I)=GAUSS1(FNJET,0.0,20.0,0.01)/HINT1(12)
5722 HINT1(10)=HINT1(10)*HIPR1(31)
5723 HINT1(12)=HINT1(12)*HIPR1(31)
5724 HINT1(13)=HINT1(13)*HIPR1(31)
5725 HINT1(59)=HINT1(59)*HIPR1(31)
5726 C ********Total and Inel cross section are calculated
5727 C by Gaussian integration.
5728 IF(IHPR2(13).NE.0) THEN
5729 HIPR1(33)=1.36*(1.0+36.0/HINT1(1)**2)
5730 & *ALOG(0.6+0.1*HINT1(1)**2)
5731 HIPR1(33)=HIPR1(33)/HINT1(12)
5733 C ********Parametrized cross section for single
5734 C diffractive reaction(Goulianos)
5742 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5745 OMG=OMG0(X)*(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
5746 FTOT=2.0*(1.0-EXP(-OMG))
5753 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5756 OMG=OMG0(X)*(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
5757 FHIN=1.0-EXP(-2.0*OMG)
5764 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5767 OMG=OMG0(X)*HINT1(11)/HIPR1(31)/2.0
5768 FTOTJT=1.0-EXP(-2.0*OMG)
5775 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5778 OMG=OMG0(X)*HINT1(60)/HIPR1(31)/2.0
5779 FTOTRG=1.0-EXP(-2.0*OMG)
5787 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5792 OMG1=OMG0(X)*HINT1(11)/HIPR1(31)
5793 C0=EXP(N*ALOG(OMG1)-SGMIN(N+1))
5794 IF(N.EQ.0) C0=1.0-EXP(-2.0*OMG0(X)*HIPR1(30)/HIPR1(31)/2.0)
5818 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
5824 X4=HIPR1(32)*SQRT(X)
5825 OMG0=HIPR1(32)**2*GAUSS2(BK,X4,X4+20.0,0.01)/96.0
5832 C ********This gives the eikonal function from a table
5833 C calculated in the first call
5834 DIMENSION FR(0:1000)
5835 clin-10/29/02 unsaved FR causes wrong values for ROMG with f77 compiler:
5840 IF(I0.NE.0) GO TO 100
5851 ROMG=(FR(IX)*((IX+1)*0.01-X)+FR(IX+1)*(X-IX*0.01))/0.01
5861 BK=EXP(-X)*(X**2-X4**2)**2.50/15.0
5866 C THIS PROGRAM IS TO CALCULATE THE JET CROSS SECTION
5867 C THE INTEGRATION IS DONE BY USING VEGAS
5870 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5871 REAL HIPR1(100),HINT1(100)
5872 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5876 COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
5878 COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
5880 COMMON/BVEG3/F,TI,TSI
5884 EXTERNAL FJET,FJETRG
5887 c************************
5888 c NCALL give the number of inner-iteration, ITMX
5889 C gives the limit of out-iteration. Nprn is an option
5890 C ( 1: print the integration process. 0: do not print)
5894 CALL VEGAS(FJET,AVGI,SD,CHI2A)
5895 HINT1(14)=sngl(AVGI)/2.5682
5896 IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
5898 CALL VEGAS(FJET,AVGI,SD,CHI2A)
5899 HINT1(15)=sngl(AVGI)/2.5682
5901 IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
5903 CALL VEGAS(FJET,AVGI,SD,CHI2A)
5904 HINT1(16)=sngl(AVGI)/2.5682
5906 IF(IHPR2(6).EQ.1.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
5908 CALL VEGAS(FJET,AVGI,SD,CHI2A)
5909 HINT1(17)=sngl(AVGI)/2.5682
5911 C ********Total inclusive jet cross section(Pt>P0)
5913 IF(IHPR2(3).NE.0) THEN
5915 CALL VEGAS(FJETRG,AVGI,SD,CHI2A)
5916 HINT1(61)=sngl(AVGI)/2.5682
5917 IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
5919 CALL VEGAS(FJETRG,AVGI,SD,CHI2A)
5920 HINT1(62)=sngl(AVGI)/2.5682
5922 IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
5924 CALL VEGAS(FJETRG,AVGI,SD,CHI2A)
5925 HINT1(63)=sngl(AVGI)/2.5682
5927 IF(IHPR2(6).EQ.1.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
5929 CALL VEGAS(FJETRG,AVGI,SD,CHI2A)
5930 HINT1(64)=sngl(AVGI)/2.5682
5933 C ********cross section of trigger jet
5940 FUNCTION FJET(X,WGT)
5941 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5942 REAL HIPR1(100),HINT1(100)
5943 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5948 PT2=dble(HINT1(1)**2/4.0-HIPR1(8)**2)*X(1)+dble(HIPR1(8))**2
5949 XT=2.0d0*DSQRT(PT2)/dble(HINT1(1))
5950 YMX1=DLOG(1.0d0/XT+DSQRT(1.0d0/XT**2-1.0d0))
5951 Y1=2.0d0*YMX1*X(2)-YMX1
5952 YMX2=DLOG(2.0d0/XT-DEXP(Y1))
5953 YMN2=DLOG(2.0d0/XT-DEXP(-Y1))
5954 Y2=(YMX2+YMN2)*X(3)-YMN2
5955 FJET=2.0d0*YMX1*(YMX2+YMN2)*dble(HINT1(1)**2/4.0-HIPR1(8)**2)
5956 & *G(Y1,Y2,PT2)/2.0d0
5962 FUNCTION FJETRG(X,WGT)
5963 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5964 REAL HIPR1(100),HINT1(100),PTMAX,PTMIN
5965 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
5970 PTMIN=ABS(HIPR1(10))-0.25
5971 PTMIN=MAX(PTMIN,HIPR1(8))
5973 IF(IHPR2(3).EQ.3) THEN
5974 AM2=dble(HIPR1(7)**2)
5975 PTMIN=MAX(0.0,HIPR1(10))
5977 PTMAX=ABS(HIPR1(10))+0.25
5978 IF(HIPR1(10).LE.0.0) PTMAX=HINT1(1)/2.0-sngl(AM2)
5979 IF(PTMAX.LE.PTMIN) PTMAX=PTMIN+0.25
5980 PT2=dble(PTMAX**2-PTMIN**2)*X(1)+dble(PTMIN)**2
5982 XT=2.0d0*DSQRT(AMT2)/dble(HINT1(1))
5983 YMX1=DLOG(1.0d0/XT+DSQRT(1.0d0/XT**2-1.0d0))
5984 Y1=2.0d0*YMX1*X(2)-YMX1
5985 YMX2=DLOG(2.0d0/XT-DEXP(Y1))
5986 YMN2=DLOG(2.0d0/XT-DEXP(-Y1))
5987 Y2=(YMX2+YMN2)*X(3)-YMN2
5988 IF(IHPR2(3).EQ.3) THEN
5989 GTRIG=2.0d0*GHVQ(Y1,Y2,AMT2)
5990 ELSE IF(IHPR2(3).EQ.2) THEN
5991 GTRIG=2.0d0*GPHOTN(Y1,Y2,PT2)
5995 FJETRG=2.0d0*YMX1*(YMX2+YMN2)*dble(PTMAX**2-PTMIN**2)
6002 FUNCTION GHVQ(Y1,Y2,AMT2)
6003 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6004 REAL HIPR1(100),HINT1(100)
6005 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
6009 XT=2.0d0*DSQRT(AMT2)/dble(HINT1(1))
6010 X1=0.5d0*XT*(DEXP(Y1)+DEXP(Y2))
6011 X2=0.5d0*XT*(DEXP(-Y1)+DEXP(-Y2))
6012 SS=X1*X2*dble(HINT1(1))**2
6014 IF(IHPR2(18).NE.0) AF=5.0d0
6015 DLAM=dble(HIPR1(15))
6016 APH=12.0d0*3.1415926d0/(33.d0-2.d0*AF)/DLOG(AMT2/DLAM**2)
6018 CALL PARTON(F,X1,X2,AMT2)
6020 Gqq=4.d0*(DCOSH(Y1-Y2)+dble(HIPR1(7))**2/AMT2)
6021 & /(1.D0+DCOSH(Y1-Y2))
6022 & /9.d0*(F(1,1)*F(2,2)+F(1,2)*F(2,1)+F(1,3)*F(2,4)
6023 & +F(1,4)*F(2,3)+F(1,5)*F(2,6)+F(1,6)*F(2,5))
6024 Ggg=(8.D0*DCOSH(Y1-Y2)-1.D0)
6025 & *(DCOSH(Y1-Y2)+2.d0*dble(HIPR1(7))**2
6026 & /AMT2-2.d0*dble(HIPR1(7))**4/AMT2**2)/(1.d0+DCOSH(Y1-Y2))
6027 & /24.d0*F(1,7)*F(2,7)
6029 GHVQ=(Gqq+Ggg)*dble(HIPR1(23))*3.14159d0*APH**2/SS**2
6035 FUNCTION GPHOTN(Y1,Y2,PT2)
6036 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6037 REAL HIPR1(100),HINT1(100)
6038 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
6042 XT=2.d0*DSQRT(PT2)/dble(HINT1(1))
6043 X1=0.5d0*XT*(DEXP(Y1)+DEXP(Y2))
6044 X2=0.5d0*XT*(DEXP(-Y1)+DEXP(-Y2))
6045 Z=DSQRT(1.D0-XT**2/X1/X2)
6046 SS=X1*X2*dble(HINT1(1))**2
6050 DLAM=dble(HIPR1(15))
6051 APH=12.d0*3.1415926d0/(33.d0-2.d0*AF)/DLOG(PT2/DLAM**2)
6054 CALL PARTON(F,X1,X2,PT2)
6056 G11=-(U**2+1.d0)/U/3.d0*F(1,7)*(4.d0*F(2,1)+4.d0*F(2,2)
6057 & +F(2,3)+F(2,4)+F(2,5)+F(2,6))/9.d0
6058 G12=-(T**2+1.d0)/T/3.d0*F(2,7)*(4.d0*F(1,1)+4.d0*F(1,2)
6059 & +F(1,3)+F(1,4)+F(1,5)+F(1,6))/9.d0
6060 G2=8.d0*(U**2+T**2)/U/T/9.d0*(4.d0*F(1,1)*F(2,2)
6061 & +4.d0*F(1,2)*F(2,1)+F(1,3)*F(2,4)+F(1,4)*F(2,3)
6062 & +F(1,5)*F(2,6)+F(1,6)*F(2,5))/9.d0
6064 GPHOTN=(G11+G12+G2)*dble(HIPR1(23))*3.14159d0*APH*APHEM/SS**2
6071 FUNCTION G(Y1,Y2,PT2)
6072 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6073 REAL HIPR1(100),HINT1(100)
6074 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
6078 XT=2.d0*DSQRT(PT2)/dble(HINT1(1))
6079 X1=0.5d0*XT*(DEXP(Y1)+DEXP(Y2))
6080 X2=0.5d0*XT*(DEXP(-Y1)+DEXP(-Y2))
6081 Z=DSQRT(1.D0-XT**2/X1/X2)
6082 SS=X1*X2*dble(HINT1(1))**2
6086 DLAM=dble(HIPR1(15))
6087 APH=12.d0*3.1415926d0/(33.d0-2.d0*AF)/DLOG(PT2/DLAM**2)
6089 CALL PARTON(F,X1,X2,PT2)
6091 G11=( (F(1,1)+F(1,2))*(F(2,3)+F(2,4)+F(2,5)+F(2,6))
6092 & +(F(1,3)+F(1,4))*(F(2,5)+F(2,6)) )*SUBCR1(T,U)
6094 G12=( (F(2,1)+F(2,2))*(F(1,3)+F(1,4)+F(1,5)+F(1,6))
6095 & +(F(2,3)+F(2,4))*(F(1,5)+F(1,6)) )*SUBCR1(U,T)
6097 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)
6098 & +F(1,5)*F(2,5)+F(1,6)*F(2,6))*(SUBCR1(U,T)
6099 & +SUBCR1(T,U)-8.D0/T/U/27.D0)
6101 G2=(AF-1)*(F(1,1)*F(2,2)+F(2,1)*F(1,2)+F(1,3)*F(2,4)
6102 & +F(2,3)*F(1,4)+F(1,5)*F(2,6)+F(2,5)*F(1,6))*SUBCR2(T,U)
6104 G31=(F(1,1)*F(2,2)+F(1,3)*F(2,4)+F(1,5)*F(2,6))*SUBCR3(T,U)
6105 G32=(F(2,1)*F(1,2)+F(2,3)*F(1,4)+F(2,5)*F(1,6))*SUBCR3(U,T)
6107 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)+
6108 1 F(1,5)*F(2,6)+F(2,5)*F(1,6))*SUBCR4(T,U)
6110 G5=AF*F(1,7)*F(2,7)*SUBCR5(T,U)
6112 G61=F(1,7)*(F(2,1)+F(2,2)+F(2,3)+F(2,4)+F(2,5)
6113 & +F(2,6))*SUBCR6(T,U)
6114 G62=F(2,7)*(F(1,1)+F(1,2)+F(1,3)+F(1,4)+F(1,5)
6115 & +F(1,6))*SUBCR6(U,T)
6117 G7=F(1,7)*F(2,7)*SUBCR7(T,U)
6119 G=(G11+G12+G13+G2+G31+G32+G4+G5+G61+G62+G7)*dble(HIPR1(17))*
6120 1 3.14159D0*APH**2/SS**2
6126 FUNCTION SUBCR1(T,U)
6127 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6128 SUBCR1=4.D0/9.D0*(1.D0+U**2)/T**2
6133 FUNCTION SUBCR2(T,U)
6134 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6135 SUBCR2=4.D0/9.D0*(T**2+U**2)
6140 FUNCTION SUBCR3(T,U)
6141 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6142 SUBCR3=4.D0/9.D0*(T**2+U**2+(1.D0+U**2)/T**2
6143 1 -2.D0*U**2/3.D0/T)
6148 FUNCTION SUBCR4(T,U)
6149 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6150 SUBCR4=8.D0/3.D0*(T**2+U**2)*(4.D0/9.D0/T/U-1.D0)
6156 FUNCTION SUBCR5(T,U)
6157 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6158 SUBCR5=3.D0/8.D0*(T**2+U**2)*(4.D0/9.D0/T/U-1.D0)
6163 FUNCTION SUBCR6(T,U)
6164 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6165 SUBCR6=(1.D0+U**2)*(1.D0/T**2-4.D0/U/9.D0)
6170 FUNCTION SUBCR7(T,U)
6171 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6172 SUBCR7=9.D0/2.D0*(3.D0-T*U-U/T**2-T/U**2)
6178 SUBROUTINE PARTON(F,X1,X2,QQ)
6179 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6180 REAL HIPR1(100),HINT1(100)
6181 COMMON/HPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
6187 DLAM=dble(HIPR1(15))
6189 S=DLOG(DLOG(QQ/DLAM**2)/DLOG(Q0**2/DLAM**2))
6190 IF(IHPR2(7).EQ.2) GO TO 200
6191 C*******************************************************
6192 AT1=0.419d0+0.004d0*S-0.007d0*S**2
6193 AT2=3.460d0+0.724d0*S-0.066d0*S**2
6194 GMUD=4.40d0-4.86d0*S+1.33d0*S**2
6195 AT3=0.763d0-0.237d0*S+0.026d0*S**2
6196 AT4=4.00d0+0.627d0*S-0.019d0*S**2
6197 GMD=-0.421d0*S+0.033d0*S**2
6198 C*******************************************************
6199 CAS=1.265d0-1.132d0*S+0.293d0*S**2
6200 AS=-0.372d0*S-0.029d0*S**2
6201 BS=8.05d0+1.59d0*S-0.153d0*S**2
6202 APHS=6.31d0*S-0.273d0*S**2
6203 BTAS=-10.5d0*S-3.17d0*S**2
6204 GMS=14.7d0*S+9.80d0*S**2
6205 C********************************************************
6206 C CAC=0.135*S-0.075*S**2
6207 C AC=-0.036-0.222*S-0.058*S**2
6208 C BC=6.35+3.26*S-0.909*S**2
6209 C APHC=-3.03*S+1.50*S**2
6210 C BTAC=17.4*S-11.3*S**2
6211 C GMC=-17.9*S+15.6*S**2
6212 C***********************************************************
6213 CAG=1.56d0-1.71d0*S+0.638d0*S**2
6214 AG=-0.949d0*S+0.325d0*S**2
6215 BG=6.0d0+1.44d0*S-1.05d0*S**2
6216 APHG=9.0d0-7.19d0*S+0.255d0*S**2
6217 BTAG=-16.5d0*S+10.9d0*S**2
6218 GMG=15.3d0*S-10.1d0*S**2
6220 C********************************************************
6221 200 AT1=0.374d0+0.014d0*S
6222 AT2=3.33d0+0.753d0*S-0.076d0*S**2
6223 GMUD=6.03d0-6.22d0*S+1.56d0*S**2
6224 AT3=0.761d0-0.232d0*S+0.023d0*S**2
6225 AT4=3.83d0+0.627d0*S-0.019d0*S**2
6226 GMD=-0.418d0*S+0.036d0*S**2
6227 C************************************
6228 CAS=1.67d0-1.92d0*S+0.582d0*S**2
6229 AS=-0.273d0*S-0.164d0*S**2
6230 BS=9.15d0+0.530d0*S-0.763d0*S**2
6231 APHS=15.7d0*S-2.83d0*S**2
6232 BTAS=-101.0d0*S+44.7d0*S**2
6233 GMS=223.0d0*S-117.0d0*S**2
6234 C*********************************
6235 C CAC=0.067*S-0.031*S**2
6236 C AC=-0.120-0.233*S-0.023*S**2
6237 C BC=3.51+3.66*S-0.453*S**2
6238 C APHC=-0.474*S+0.358*S**2
6239 C BTAC=9.50*S-5.43*S**2
6240 C GMC=-16.6*S+15.5*S**2
6241 C**********************************
6242 CAG=0.879d0-0.971d0*S+0.434d0*S**2
6243 AG=-1.16d0*S+0.476d0*S**2
6244 BG=4.0d0+1.23d0*S-0.254d0*S**2
6245 APHG=9.0d0-5.64d0*S-0.817d0*S**2
6246 BTAG=-7.54d0*S+5.50d0*S**2
6247 GMG=-0.596d0*S+1.26d0*S**2
6248 C*********************************
6249 300 B12=DEXP(GMRE(AT1)+GMRE(AT2+1.D0)-GMRE(AT1+AT2+1.D0))
6250 B34=DEXP(GMRE(AT3)+GMRE(AT4+1.D0)-GMRE(AT3+AT4+1.D0))
6251 CNUD=3.D0/B12/(1.D0+GMUD*AT1/(AT1+AT2+1.D0))
6252 CND=1.D0/B34/(1.D0+GMD*AT3/(AT3+AT4+1.D0))
6253 C********************************************************
6255 C FS=X*2(UBAR+DBAR+SBAR) AND UBAR=DBAR=SBAR
6256 C*******************************************************
6257 FUD1=CNUD*X1**AT1*(1.D0-X1)**AT2*(1.D0+GMUD*X1)
6258 FS1=CAS*X1**AS*(1.D0-X1)**BS*(1.D0+APHS*X1
6259 & +BTAS*X1**2+GMS*X1**3)
6260 F(1,3)=CND*X1**AT3*(1.D0-X1)**AT4*(1.D0+GMD*X1)+FS1/6.D0
6261 F(1,1)=FUD1-F(1,3)+FS1/3.D0
6266 F(1,7)=CAG*X1**AG*(1.D0-X1)**BG*(1.D0+APHG*X1
6267 & +BTAG*X1**2+GMG*X1**3)
6269 FUD2=CNUD*X2**AT1*(1.D0-X2)**AT2*(1.D0+GMUD*X2)
6270 FS2=CAS*X2**AS*(1.D0-X2)**BS*(1.D0+APHS*X2
6271 & +BTAS*X2**2+GMS*X2**3)
6272 F(2,3)=CND*X2**AT3*(1.D0-X2)**AT4*(1.D0+GMD*X2)+FS2/6.D0
6273 F(2,1)=FUD2-F(2,3)+FS2/3.D0
6278 F(2,7)=CAG*X2**AG*(1.D0-X2)**BG*(1.D0+APHG*X2
6279 & +BTAG*X2**2+GMG*X2**3)
6280 C***********Nuclear effect on the structure function****************
6282 IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
6283 AAX=1.193d0*dble(ALOG(FLOAT(IHNT2(1)))**0.16666666)
6284 RRX=AAX*(X1**3-1.2d0*X1**2+0.21d0*X1)+1.d0
6285 & +dble(1.079*(FLOAT(IHNT2(1))**0.33333333-1.0))
6286 & /dble(ALOG(float(IHNT2(1))+1.0))*DSQRT(X1)
6287 & *DEXP(-X1**2/0.01d0)
6288 c & /DLOG(IHNT2(1)+1.0D0)*(DSQRT(X1)*DEXP(-X1**2/0.01)
6289 IF(ipcrs.EQ.1 .OR.ipcrs.EQ.3) RRX=DEXP(-X1**2/0.01d0)
6294 IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
6295 AAX=1.193d0*dble(ALOG(FLOAT(IHNT2(3)))**0.16666666)
6296 RRX=AAX*(X2**3-1.2d0*X2**2+0.21d0*X2)+1.d0
6297 & +dble(1.079*(FLOAT(IHNT2(3))**0.33333-1.0))
6298 & /dble(ALOG(float(IHNT2(3))+1.0))*DSQRT(X2)
6299 & *DEXP(-X2**2/0.01d0)
6300 c & /DLOG(IHNT2(3)+1.0D0)*DSQRT(X2)*DEXP(-X2**2/0.01)
6301 IF(ipcrs.EQ.2 .OR. ipcrs.EQ.3) RRX=DEXP(-X2**2/0.01d0)
6313 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6316 IF(X.GT.3.0D0) GO TO 10
6318 10 GMRE=0.5D0*DLOG(2.D0*3.14159265D0/Z)+Z*DLOG(Z)-Z+DLOG(1.D0
6319 1 +1.D0/12.D0/Z+1.D0/288.D0/Z**2-139.D0/51840.D0/Z**3
6320 1 -571.D0/2488320.D0/Z**4)
6322 GMRE=GMRE-DLOG(Z-1.D0)-DLOG(Z-2.D0)-DLOG(Z-3.D0)
6329 C***************************************************************
6332 PARAMETER (MAXSTR=150001)
6333 DOUBLE PRECISION XL(10),XU(10),ACC
6334 COMMON/BVEG1/XL,XU,ACC,NDIM,NCALL,ITMX,NPRN
6338 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
6340 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
6342 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
6344 COMMON/HSTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
6346 COMMON/hjcrdn/YP(3,300),YT(3,300)
6348 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
6349 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
6350 & PJPM(300,500),NTJ(300),KFTJ(300,500),
6351 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
6352 & PJTE(300,500),PJTM(300,500)
6354 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
6355 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
6356 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
6358 COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
6360 COMMON/HPINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
6362 cwei DOUBLE PRECISION PATT
6364 DATA NUM1/30123984/,XL/10*0.D0/,XU/10*1.D0/
6365 DATA NCALL/1000/,ITMX/100/,ACC/0.01/,NPRN/0/
6366 C...give all the switchs and parameters the default values
6367 clin-4/2008 input.ampt provides NSEED for AMPT:
6368 c DATA NSEED/74769375/
6370 & 1.5, 0.35, 0.5, 0.9, 2.0, 0.1, 1.5, 2.0, -1.0, -2.25,
6371 & 2.0, 0.5, 1.0, 2.0, 0.2, 2.0, 2.5, 0.3, 0.1, 1.4,
6372 & 1.6, 1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 57.0,
6373 & 28.5, 3.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
6375 & 0.0, 0.4, 0.1, 1.5, 0.1, 0.25, 0.0, 0.5, 0.0, 0.0,
6379 & 1, 3, 0, 1, 1, 1, 1, 10, 0, 0,
6380 & 1, 1, 1, 1, 0, 0, 1, 0, 0, 1,
6386 C...initialize all the data common blocks
6387 DATA NATT/0/,EATT/0.0/,JATT/0/,NT/0/,NP/0/,
6388 1 N0/0/,N01/0/,N10/0/,N11/0/
6390 c DATA KATT/520000*0/PATT/520000*0.0/
6391 DATA KATT/600004*0/,PATT/600004*0.0/
6393 DATA NFP/4500*0/,PP/4500*0.0/,NFT/4500*0/,PT/4500*0.0/
6395 DATA YP/900*0.0/,YT/900*0.0/
6397 DATA NPJ/300*0/,KFPJ/150000*0/,PJPX/150000*0.0/,PJPY/150000*0.0/
6398 & ,PJPZ/150000*0.0/,PJPE/150000*0.0/,PJPM/150000*0.0/
6399 DATA NTJ/300*0/,KFTJ/150000*0/,PJTX/150000*0.0/,PJTY/150000*0.0/
6400 & ,PJTZ/150000*0.0/,PJTE/150000*0.0/,PJTM/150000*0.0/
6403 c DATA NSG/0/,NJSG/900*0/,IASG/2700*0/,K1SG/90000*0/,K2SG/90000*0/
6404 c & ,PXSG/90000*0.0/,PYSG/90000*0.0/,PZSG/90000*0.0/
6405 c & ,PESG/90000*0.0/,PMSG/90000*0.0/
6406 DATA NSG/0/,NJSG/150001*0/,IASG/450003*0/,
6407 & K1SG/15000100*0/,K2SG/15000100*0/,
6408 & PXSG/15000100*0.0/,PYSG/15000100*0.0/,PZSG/15000100*0.0/,
6409 & PESG/15000100*0.0/,PMSG/15000100*0.0/
6410 DATA MINT4/0/,MINT5/0/,ATCO/4000*0.0/,ATXS/201*0.0/
6411 DATA (HIDAT0(1,I),I=1,10)/0.0,0.0,0.0,0.0,0.0,0.0,2.25,
6413 DATA (HIDAT0(2,I),I=1,10)/2.0,3.0,5.0,6.0,7.0,8.0,8.0,10.0,
6415 DATA (HIDAT0(3,I),I=1,10)/1.0,0.8,0.8,0.7,0.45,0.215,
6416 & 0.21,0.19,0.19,0.19/
6417 DATA (HIDAT0(4,I),I=1,10)/0.35,0.35,0.3,0.3,0.3,0.3,
6419 DATA (HIDAT0(5,I),I=1,10)/23.8,24.0,26.0,26.2,27.0,28.5,28.5,
6421 DATA ((HIDAT0(J,I),I=1,10),J=6,9)/40*0.0/
6422 DATA (HIDAT0(10,I),I=1,10)/5.0,20.0,53.0,62.0,100.0,200.0,
6423 & 546.0,900.0,1800.0,4000.0/
6426 C*******************************************************************
6431 C*******************************************************************
6432 C SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
6433 C - BY G.P. LEPAGE SEPT 1976/(REV)APR 1978
6434 C*******************************************************************
6436 SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A)
6437 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6438 COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
6440 COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
6442 COMMON/BVEG3/F,TI,TSI
6445 DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10)
6450 DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/
6456 ENTRY VEGAS1(FXN,AVGI,SD,CHI2A)
6457 C - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
6464 ENTRY VEGAS2(FXN,AVGI,SD,CHI2A)
6465 C - NO INITIALIZATION
6468 IF(MDS.EQ.0) GO TO 2
6469 NG=int((real(NCALL)/2.)**(1./real(NDIM)))
6471 IF((2*NG-NDMX).LT.0) GO TO 2
6481 DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
6487 c***this is the line 50
6491 C REBIN PRESERVING BIN DENSITY
6493 IF(ND.EQ.NDO) GO TO 8
6504 5 IF(RC.GT.DR) GO TO 4
6507 XIN(I)=XN-(XN-XO)*DR
6508 IF(I.LT.NDM) GO TO 5
6515 c IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND
6516 c 1 ,(XL(J),XU(J),J=1,NDIM)
6518 ENTRY VEGAS3(FXN,AVGI,SD,CHI2A)
6519 C - MAIN INTEGRATION LOOP
6533 CALL ARAN9(QRAN,NDIM)
6536 XN=dble(float(KG(J))-QRAN(J))*DXG+ONE
6537 c*****this is the line 100
6539 IF(IA(J).GT.1) GO TO 13
6543 13 XO=XI(IA(J),J)-XI(IA(J)-1,J)
6544 RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
6545 14 X(J)=XL(J)+RC*DX(J)
6555 DI(IA(J),J)=DI(IA(J),J)+F
6556 16 IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
6557 IF(K.LT.NPG) GO TO 12
6560 F2B=(F2B-FB)*(F2B+FB)
6563 IF(MDS.GE.0) GO TO 18
6565 17 D(IA(J),J)=D(IA(J),J)+F2B
6567 19 KG(K)=MOD(KG(K),NG)+1
6568 IF(KG(K).NE.1) GO TO 11
6572 C FINAL RESULTS FOR THIS ITERATION
6576 WGT=TI2/(TSI+1.0d-37)
6585 CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/dble(float(IT)-.999)
6587 C****this is the line 150
6588 IF(NPRN.EQ.0) GO TO 21
6590 c WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A
6591 c IF(NPRN.GE.0) GO TO 21
6593 c20 WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND)
6606 D(I,J)=(D(I,J)+XN)/3.d0
6607 22 DT(J)=DT(J)+D(I,J)
6608 D(ND,J)=(XN+XO)/2.d0
6609 23 DT(J)=DT(J)+D(ND,J)
6615 IF (DT(J).GE.1.0D18) THEN
6616 WRITE(6,*) '************** A SINGULARITY >1.0D18'
6618 C1111 FORMAT(1X,'**************IMPORTANT NOTICE***************')
6620 C1112 FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
6622 C1113 FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
6624 C1114 FORMAT(1X,'**************END NOTICE*************')
6626 IF(D(I,J).LE.1.0D-18) GO TO 24
6628 R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH
6638 c****this is the line 200
6640 26 IF(RC.GT.DR) GO TO 25
6643 XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30)
6644 IF(I.LT.NDM) GO TO 26
6649 IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9
6650 c200 FORMAT('0INPUT PARAMETERS FOR VEGAS: NDIM=',I3,' NCALL=',F8.0
6651 c 1 /28X,' IT=',I5,' ITMX=',I5/28X,' ACC=',G9.3
6652 c 2 /28X,' MDS=',I3,' ND=',I4/28X,' (XL,XU)=',
6653 c 3 (T40,'( ',G12.6,' , ',G12.6,' )'))
6654 c201 FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3,
6655 c 1 ': INTEGRAL =',G14.8/21X,'STD DEV =',G10.4 /
6656 c 2 ' ACCUMULATED RESULTS: INTEGRAL =',G14.8 /
6657 c 3 24X,'STD DEV =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4)
6658 c202 FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,' DELT I ',
6659 c 1 2X,' CONV''CE ',11X,'X',7X,' DELT I ',2X,' CONV''CE '
6660 c 2 ,11X,'X',7X,' DELT I ',2X,' CONV''CE ' /
6661 c 2 (' ',3G12.4,5X,3G12.4,5X,3G12.4))
6666 SUBROUTINE ARAN9(QRAN,NDIM)
6671 1 QRAN(I)=RANART(NUM1)
6677 C*********GAUSSIAN ONE-DIMENSIONAL INTEGRATION PROGRAM*************
6679 FUNCTION GAUSS1(F,A,B,EPS)
6681 DIMENSION W(12),X(12)
6684 DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
6685 & .0622535,0.0951585,.1246290,.1495960,.1691565,
6686 & .1826034,.1894506/
6687 DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
6688 & .9445750,0.8656312,.7554044,.6178762,.4580168,
6689 & .2816036,.0950125/
6691 DELTA=CONST*ABS(A-B)
6695 IF(ABS(Y).LE.DELTA) RETURN
6703 1 S8=S8+W(I)*(F(C1+U)+F(C1-U))
6706 3 S16=S16+W(I)*(F(C1+U)+F(C1-U))
6709 IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
6714 IF(ABS(Y).GT.DELTA) GOTO 2
6718 7 FORMAT(1X,'GAUSS1....TOO HIGH ACURACY REQUIRED')
6723 FUNCTION GAUSS2(F,A,B,EPS)
6725 DIMENSION W(12),X(12)
6728 DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
6729 & .0622535,0.0951585,.1246290,.1495960,.1691565,
6730 & .1826034,.1894506/
6731 DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
6732 & .9445750,0.8656312,.7554044,.6178762,.4580168,
6733 & .2816036,.0950125/
6735 DELTA=CONST*ABS(A-B)
6739 IF(ABS(Y).LE.DELTA) RETURN
6747 1 S8=S8+W(I)*(F(C1+U)+F(C1-U))
6750 3 S16=S16+W(I)*(F(C1+U)+F(C1-U))
6753 IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
6758 IF(ABS(Y).GT.DELTA) GOTO 2
6762 7 FORMAT(1X,'GAUSS2....TOO HIGH ACURACY REQUIRED')
6778 c & '**************************************************'/10X,
6779 c & '* | \ _______ / ------/ *'/10X,
6780 c & '* ----- ------ |_____| /_/ / *'/10X,
6781 c & '* ||\ / |_____| / / \ *'/10X,
6782 c & '* /| \ /_/ /_______ /_ / \_ *'/10X,
6783 c & '* / | / / / / / | ------- *'/10X,
6784 c & '* | / /\ / / | / | *'/10X,
6785 c & '* | / / \ / / \_| / ------- *'/10X,
6787 & '**************************************************'/10X,
6788 & '* | | _______ / ------/ *'/10X,
6789 & '* ----- ------ |_____| /_/ / *'/10X,
6790 & '* ||| / |_____| / / | *'/10X,
6791 & '* /| | /_/ /_______ /_ / | *'/10X,
6792 & '* / | / / / / / | ------- *'/10X,
6793 & '* | / /| / / | / | *'/10X,
6794 & '* | / / | / / _| / ------- *'/10X,
6796 & '**************************************************'/10X,
6798 & ' Heavy Ion Jet INteraction Generator '/10X,
6800 & ' X. N. Wang and M. Gyulassy '/10X,
6801 & ' Lawrence Berkeley Laboratory '//)