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