5 SUBROUTINE HIJHRD(JP,JT,JOUT,JFLG,IOPJET0)
7 C IOPTJET=1, ALL JET WILL FORM SINGLE STRING SYSTEM
8 C 0, ONLY Q-QBAR JET FORM SINGLE STRING SYSTEM
9 C*******Perform jets production and fragmentation when JP JT *******
10 C scatter. JOUT-> number of hard scatterings precede this one *
11 C for the the same pair(JP,JT). JFLG->a flag to show whether *
12 C jets can be produced (with valence quark=1,gluon=2, q-qbar=3)*
13 C or not(0). Information of jets are in COMMON/ATTJET and *
14 C /MINJET. ABS(NFP(JP,6)) is the total number jets produced by *
15 C JP. If NFP(JP,6)<0 JP can not produce jet anymore. *
16 C*******************************************************************
17 DIMENSION IP(100,2),IPQ(50),IPB(50),IT(100,2),ITQ(50),ITB(50)
18 #include "hijcrdn.inc"
19 #include "hiparnt.inc"
21 #include "histrng.inc"
22 #include "hijjet1.inc"
23 #include "hijjet2.inc"
24 #include "hijjet4.inc"
25 C************************************ HIJING common block
26 #include "lujets_hijing.inc"
27 #include "ludat1_hijing.inc"
28 #include "pysubs_hijing.inc"
29 #include "pypars_hijing.inc"
30 #include "pyint1_hijing.inc"
31 #include "pyint2_hijing.inc"
32 #include "pyint5_hijing.inc"
33 #include "hipyint.inc"
35 C*********************************** LU common block
37 C SIZE OF COMMON BLOCK FOR # OF PARTON PER STRING
39 C SIZE OF COMMON BLOCK FOR # OF SINGLE STRINGS
41 C SIZE OF COMMON BLOCK FOR # OF PARTON PER SINGLE
48 IF(IOPJET.EQ.1.AND.(NFP(JP,6).NE.0.OR.NFT(JT,6).NE.0))
50 IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
51 IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) RETURN
52 C ******** JP or JT can not produce jet anymore
59 IF(EPP.LT.0.0) GO TO 1000
60 IF(EPM.LT.0.0) GO TO 1000
61 IF(ETP.LT.0.0) GO TO 1000
62 IF(ETM.LT.0.0) GO TO 1000
63 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
65 C ********for the first hard scattering of (JP,JT)
66 C have collision only when Ycm(JP)>Ycm(JT)
68 ECUT1=HIPR1(1)+HIPR1(8)+PP(JP,14)+PP(JP,15)
69 ECUT2=HIPR1(1)+HIPR1(8)+PT(JT,14)+PT(JT,15)
70 IF(PP(JP,4).LE.ECUT1) THEN
71 NFP(JP,6)=-ABS(NFP(JP,6))
74 IF(PT(JT,4).LE.ECUT2) THEN
75 NFT(JT,6)=-ABS(NFT(JT,6))
78 C *********must have enough energy to produce jets
84 IF(NFP(JP,10).EQ.0 .AND. NFT(JT,10).EQ.0) THEN
100 IF(XSEC(11,1).NE.0) ISUB11=1
101 IF(XSEC(12,1).NE.0) ISUB12=1
102 IF(XSEC(28,1).NE.0) ISUB28=1
103 MINT(44)=MINT4-ISUB11-ISUB12-ISUB28
104 MINT(45)=MINT5-ISUB11-ISUB12-ISUB28
105 XSEC(0,1)=ATXS(0)-ATXS(11)-ATXS(12)-ATXS(28)
115 C ********Scatter the valence quarks only once per NN
117 C afterwards only gluon can have hard scattering.
118 155 CALL PYTHIA_HIJING
120 IF(JJ.NE.1) GO TO 155
121 C *********one hard collision at a time
122 IF(K(7,2).EQ.-K(8,2)) THEN
123 QMASS2=(P(7,4)+P(8,4))**2-(P(7,1)+P(8,1))**2
124 & -(P(7,2)+P(8,2))**2-(P(7,3)+P(8,3))**2
125 QM=ULMASS_HIJING(K(7,2))
126 IF(QMASS2.LT.(2.0*QM+HIPR1(1))**2) GO TO 155
128 C ********q-qbar jets must has minimum mass HIPR1(1)
138 IF(PEP.LE.ECUT1) THEN
140 IF(MISP.LT.50) GO TO 155
141 NFP(JP,6)=-ABS(NFP(JP,6))
144 IF(PET.LE.ECUT2) THEN
146 IF(MIST.LT.50) GO TO 155
147 NFT(JT,6)=-ABS(NFT(JT,6))
150 C ******** if the remain energy<ECUT the proj or targ
151 C can not produce jet anymore
155 IF(WP.LT.0.0 .OR. WM.LT.0.0) THEN
157 IF(MISS.LT.50) GO TO 155
160 C ********the total W+, W- must be positive
162 AMPX=SQRT((ECUT1-HIPR1(8))**2+PXP**2+PYP**2+0.01)
163 AMTX=SQRT((ECUT2-HIPR1(8))**2+PXT**2+PYT**2+0.01)
165 IF(SW.LT.SXX.OR.VINT(43).LT.HIPR1(1)) THEN
167 IF(MISS.LT.50) GO TO 155
170 C ********the proj and targ remnants must have at least
171 C a CM energy that can produce two strings
172 C with minimum mass HIPR1(1)(see HIJSFT HIJFRG)
179 HINT1(46)=SQRT(P(7,1)**2+P(7,2)**2)
185 HINT1(56)=SQRT(P(8,1)**2+P(8,2)**2)
189 PINIRAD=(1.0-EXP(-2.0*(VINT(47)-HIDAT(1))))
190 & /(1.0+EXP(-2.0*(VINT(47)-HIDAT(1))))
192 IF(RLU_HIJING(0).LE.PINIRAD) I_INIRAD=1
193 IF(K(7,2).EQ.-K(8,2)) GO TO 190
194 IF(K(7,2).EQ.21.AND.K(8,2).EQ.21.AND.IOPJET.EQ.1) GO TO 190
195 C*******************************************************************
196 C gluon jets are going to be connectd with
197 C the final leading string of quark-aintquark
198 C*******************************************************************
217 IF(K(I,3).EQ.1 .OR. K(I,3).EQ.2.OR.
218 & ABS(K(I,2)).GT.30) GO TO 180
219 C************************************************************
221 HINT1(47)=HINT1(47)+P(I,1)
222 HINT1(48)=HINT1(48)+P(I,2)
223 HINT1(49)=HINT1(49)+P(I,3)
224 HINT1(50)=HINT1(50)+P(I,4)
227 HINT1(67)=HINT1(67)+P(I,1)
228 HINT1(68)=HINT1(68)+P(I,2)
229 HINT1(69)=HINT1(69)+P(I,3)
230 HINT1(70)=HINT1(70)+P(I,4)
232 C************************modifcation made on Apr 10. 1996*****
233 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
243 C************************************************************
245 C************************correction made on Oct. 14,1994*****
247 IF(K(I,3).EQ.7.OR.K(I,3).EQ.3) THEN
248 IF(K(I,3).EQ.7.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(7,2)
249 & .AND.IS7.EQ.0) THEN
259 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
260 & I_INIRAD.EQ.0)) THEN
270 IF(K(I,2).NE.21) THEN
275 ELSE IF(K(I,2).LT.0) THEN
281 ELSE IF(K(I,3).EQ.8.OR.K(I,3).EQ.4) THEN
282 IF(K(I,3).EQ.8.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(8,2)
283 & .AND.IS8.EQ.0) THEN
293 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
294 & I_INIRAD.EQ.0)) THEN
304 IF(K(I,2).NE.21) THEN
309 ELSE IF(K(I,2).LT.0) THEN
319 IF(LPQ.NE.LPB .OR. LTQ.NE.LTB) THEN
321 IF(MISS.LE.50) GO TO 155
322 WRITE(6,*) ' Q -QBAR NOT MATCHED IN HIJHRD'
326 C****The following will rearrange the partons so that a quark is***
327 C****allways followed by an anti-quark ****************************
331 IF(J.GT.JPP) GO TO 182
332 IF(IP(J,2).EQ.0) THEN
334 ELSE IF(IP(J,2).NE.0) THEN
338 IP(J,1)=IP(IPQ(LP),1)
339 IP(J,2)=IP(IPQ(LP),2)
344 ELSE IF(IP2.LT.0) THEN
347 C ********replace J with a quark
350 IP(J+1,1)=IP(IPB(LP),1)
351 IP(J+1,2)=IP(IPB(LP),2)
356 ELSE IF(IP2.LT.0) THEN
359 C ******** replace J+1 with anti-quark
366 IF(J.GT.JTT) GO TO 184
367 IF(IT(J,2).EQ.0) THEN
369 ELSE IF(IT(J,2).NE.0) THEN
373 IT(J,1)=IT(ITQ(LT),1)
374 IT(J,2)=IT(ITQ(LT),2)
379 ELSE IF(IT2.LT.0) THEN
382 C ********replace J with a quark
385 IT(J+1,1)=IT(ITB(LT),1)
386 IT(J+1,2)=IT(ITB(LT),2)
391 ELSE IF(IT2.LT.0) THEN
394 C ******** replace J+1 with anti-quark
401 IF(NPJ(JP)+JPP.GT.MXJT.OR.NTJ(JT)+JTT.GT.MXJT) THEN
403 WRITE(6,*) 'number of partons per string exceeds'
404 WRITE(6,*) 'the common block size'
407 C ********check the bounds of common blocks
409 KFPJ(JP,NPJ(JP)+J)=K(IP(J,1),2)
410 PJPX(JP,NPJ(JP)+J)=P(IP(J,1),1)
411 PJPY(JP,NPJ(JP)+J)=P(IP(J,1),2)
412 PJPZ(JP,NPJ(JP)+J)=P(IP(J,1),3)
413 PJPE(JP,NPJ(JP)+J)=P(IP(J,1),4)
414 PJPM(JP,NPJ(JP)+J)=P(IP(J,1),5)
418 KFTJ(JT,NTJ(JT)+J)=K(IT(J,1),2)
419 PJTX(JT,NTJ(JT)+J)=P(IT(J,1),1)
420 PJTY(JT,NTJ(JT)+J)=P(IT(J,1),2)
421 PJTZ(JT,NTJ(JT)+J)=P(IT(J,1),3)
422 PJTE(JT,NTJ(JT)+J)=P(IT(J,1),4)
423 PJTM(JT,NTJ(JT)+J)=P(IT(J,1),5)
427 C*****************************************************************
428 CThis is the case of a quark-antiquark jet it will fragment alone
429 C****************************************************************
431 IF(K(7,2).NE.21.AND.K(8,2).NE.21.AND.
432 & K(7,2)*K(8,2).GT.0) GO TO 155
437 IF(K(I,3).EQ.1.OR.K(I,3).EQ.2.OR.
438 & ABS(K(I,2)).GT.30) GO TO 200
439 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
449 C************************************************************
451 C************************correction made on Oct. 14,1994*****
453 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
454 & I_INIRAD.EQ.0)) THEN
461 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
462 & I_INIRAD.EQ.0)) THEN
472 IF(K(I,2).NE.21) THEN
477 ELSE IF(K(I,2).LT.0) THEN
486 IF(MISS.LE.50) GO TO 155
487 WRITE(6,*) LPQ,LPB, 'Q-QBAR NOT CONSERVED OR NOT MATCHED'
492 C**** The following will rearrange the partons so that a quark is***
493 C**** allways followed by an anti-quark ****************************
496 IF(J.GT.JPP) GO TO 222
497 IF(IP(J,2).EQ.0) GO TO 220
501 IP(J,1)=IP(IPQ(LP),1)
502 IP(J,2)=IP(IPQ(LP),2)
507 ELSE IF(IP2.LT.0) THEN
511 C ********replace J with a quark
514 IP(J+1,1)=IP(IPB(LP),1)
515 IP(J+1,2)=IP(IPB(LP),2)
520 ELSE IF(IP2.LT.0) THEN
523 C ******** replace J+1 with an anti-quark
533 IP(2*L0-3,1)=IP(IPQ(L0),1)
534 IP(2*L0-3,2)=IP(IPQ(L0),2)
539 ELSE IF(IP2.LT.0) THEN
546 IP(2*L0-2,1)=IP(IPB(L0),1)
547 IP(2*L0-2,2)=IP(IPB(L0),2)
552 ELSE IF(IP2.LT.0) THEN
557 C ********move all the qqbar pair to the front of
558 C the list, except the first pair
561 IP(2*LPQ-1,1)=IP(IPQ(1),1)
562 IP(2*LPQ-1,2)=IP(IPQ(1),2)
567 ELSE IF(IP2.LT.0) THEN
571 C ********move the first quark to the beginning of
572 C the last string system
575 IP(JPP,1)=IP(IPB(1),1)
576 IP(JPP,2)=IP(IPB(1),2)
581 ELSE IF(IP2.LT.0) THEN
585 C ********move the first anti-quark to the end of the
590 WRITE(6,*) 'number of jets forming single strings exceeds'
591 WRITE(6,*) 'the common block size'
596 WRITE(6,*) 'number of partons per single jet system'
597 WRITE(6,*) 'exceeds the common block size'
600 C ********check the bounds of common block size
608 K2SG(NSG,I)=K(IP(I,1),2)
609 IF(K2SG(NSG,I).LT.0) K1SG(NSG,I)=1
610 PXSG(NSG,I)=P(IP(I,1),1)
611 PYSG(NSG,I)=P(IP(I,1),2)
612 PZSG(NSG,I)=P(IP(I,1),3)
613 PESG(NSG,I)=P(IP(I,1),4)
614 PMSG(NSG,I)=P(IP(I,1),5)
618 C******* reset the energy-momentum of incoming particles ********
630 NFP(JP,6)=NFP(JP,6)+1
631 NFT(JT,6)=NFT(JT,6)+1
635 IF(IHPR2(10).EQ.0) RETURN
636 WRITE(6,*) 'Fatal HIJHRD error'
637 WRITE(6,*) JP, ' proj E+,E-',EPP,EPM,' status',NFP(JP,5)
638 WRITE(6,*) JT, ' targ E+,E_',ETP,ETM,' status',NFT(JT,5)