| 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 |
| 192 | rnd1=max(RAN0(NSEED),1.0e-20) |
| 193 | rnd2=max(RAN0(NSEED),1.0e-20) |
| 194 | rnd3=max(RAN0(NSEED),1.0e-20) |
| 195 | R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0 |
| 196 | & +4.38*0.85*log(rnd3)/(4.38+0.85)) |
| 197 | endif |
| 198 | c |
| 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 |
| 249 | rnd1=max(RAN0(NSEED),1.0e-20) |
| 250 | rnd2=max(RAN0(NSEED),1.0e-20) |
| 251 | rnd3=max(RAN0(NSEED),1.0e-20) |
| 252 | R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0 |
| 253 | & +4.38*0.85*log(rnd3)/(4.38+0.85)) |
| 254 | endif |
| 255 | c |
| 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 |
| 324 | |
| 325 | C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS |
| 326 | C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET |
| 327 | C BY THE ANGLE PHI FOR EACH COLLISION.****************** |
| 328 | C |
| 329 | BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2)) |
| 330 | PHI=2.0*HIPR1(40)*RLU_HIJING(0) |
| 331 | BBX=BB*COS(PHI) |
| 332 | BBY=BB*SIN(PHI) |
| 333 | HINT1(19)=BB |
| 334 | HINT1(20)=PHI |
| 335 | C |
| 336 | DO 70 JP=1,IHNT2(1) |
| 337 | DO 70 JT=1,IHNT2(3) |
| 338 | SCIP(JP,JT)=-1.0 |
| 339 | B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2 |
| 340 | R2=B2*HIPR1(40)/HIPR1(31)/0.1 |
| 341 | C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb |
| 342 | RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2) |
| 343 | & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0) |
| 344 | RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2) |
| 345 | & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0) |
| 346 | APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0) |
| 347 | & *SQRT(1.0-RRB1) |
| 348 | APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0) |
| 349 | & *SQRT(1.0-RRB2) |
| 350 | HINT1(18)=HINT1(14)-APHX1*HINT1(15) |
| 351 | & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17) |
| 352 | IF(IHPR2(14).EQ.0.OR. |
| 353 | & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN |
| 354 | GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31)) |
| 355 | RANTOT=RLU_HIJING(0) |
| 356 | IF(RANTOT.GT.GS) GO TO 70 |
| 357 | GO TO 65 |
| 358 | ENDIF |
| 359 | GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18)) |
| 360 | & /HIPR1(31)/2.0*ROMG(0.0))) |
| 361 | R2=R2/GSTOT_0 |
| 362 | GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2)) |
| 363 | GSTOT=2.0*(1.0-SQRT(1.0-GS)) |
| 364 | RANTOT=RLU_HIJING(0)*GSTOT_0 |
| 365 | IF(RANTOT.GT.GSTOT) GO TO 70 |
| 366 | IF(RANTOT.GT.GS) THEN |
| 367 | CALL HIJCSC(JP,JT) |
| 368 | GO TO 70 |
| 369 | C ********perform elastic collisions |
| 370 | ENDIF |
| 371 | 65 SCIP(JP,JT)=R2 |
| 372 | RNIP(JP,JT)=RANTOT |
| 373 | SJIP(JP,JT)=HINT1(18) |
| 374 | NCOLT=NCOLT+1 |
| 375 | IPCOL(NCOLT)=JP |
| 376 | ITCOL(NCOLT)=JT |
| 377 | 70 CONTINUE |
| 378 | C ********total number interactions proj and targ has |
| 379 | C suffered |
| 380 | IF(NCOLT.EQ.0) THEN |
| 381 | NLOP=NLOP+1 |
| 382 | IF(NLOP.LE.20.OR. |
| 383 | & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60 |
| 46f33c91 |
384 | |
| 385 | |
| e74335a4 |
386 | RETURN |
| 387 | ENDIF |
| 388 | C ********At large impact parameter, there maybe no |
| 389 | C interaction at all. For NN collision |
| 390 | C repeat the event until interaction happens |
| 391 | C |
| 392 | IF(IHPR2(3).NE.0) THEN |
| 393 | NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5) |
| 394 | NHARD=MIN(NHARD,NCOLT) |
| 395 | JPHARD=IPCOL(NHARD) |
| 396 | JTHARD=ITCOL(NHARD) |
| 397 | ENDIF |
| 398 | C |
| 399 | IF(IHPR2(9).EQ.1) THEN |
| 400 | NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5) |
| 401 | NMINI=MIN(NMINI,NCOLT) |
| 402 | JPMINI=IPCOL(NMINI) |
| 403 | JTMINI=ITCOL(NMINI) |
| 404 | ENDIF |
| 405 | C ********Specifying the location of the hard and |
| 406 | C minijet if they are enforced by user |
| 407 | C |
| 408 | DO 200 JP=1,IHNT2(1) |
| 409 | DO 200 JT=1,IHNT2(3) |
| 410 | IF(SCIP(JP,JT).EQ.-1.0) GO TO 200 |
| 411 | NFP(JP,11)=NFP(JP,11)+1 |
| 412 | NFT(JT,11)=NFT(JT,11)+1 |
| 413 | IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN |
| 414 | NP=NP+1 |
| 415 | N01=N01+1 |
| 416 | ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN |
| 417 | NT=NT+1 |
| 418 | N10=N10+1 |
| 419 | ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN |
| 420 | NP=NP+1 |
| 421 | NT=NT+1 |
| 422 | N0=N0+1 |
| 423 | ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN |
| 424 | N11=N11+1 |
| 425 | ENDIF |
| 426 | JOUT=0 |
| 427 | NFP(JP,10)=0 |
| 428 | NFT(JT,10)=0 |
| 429 | C***************************************************************** |
| 430 | IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160 |
| 431 | C ********When IHPR2(8)=0 no jets are produced |
| 432 | IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160 |
| 433 | C ********jets can not be produced for (JP,JT) |
| 434 | C because not enough energy avaible for |
| 435 | C JP or JT |
| 436 | R2=SCIP(JP,JT) |
| 437 | HINT1(18)=SJIP(JP,JT) |
| 438 | TT=ROMG(R2)*HINT1(18)/HIPR1(31) |
| 439 | TTS=HIPR1(30)*ROMG(R2)/HIPR1(31) |
| 440 | NJET=0 |
| 441 | IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN |
| 442 | CALL JETINI(JP,JT,1) |
| 443 | CALL HIJHRD(JP,JT,0,JFLG,0) |
| 444 | HINT1(26)=HINT1(47) |
| 445 | HINT1(27)=HINT1(48) |
| 446 | HINT1(28)=HINT1(49) |
| 447 | HINT1(29)=HINT1(50) |
| 448 | HINT1(36)=HINT1(67) |
| 449 | HINT1(37)=HINT1(68) |
| 450 | HINT1(38)=HINT1(69) |
| 451 | HINT1(39)=HINT1(70) |
| 452 | C |
| 453 | IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1 |
| 454 | IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1 |
| 455 | IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND. |
| 456 | & JFLG.GE.3) IASG(NSG,3)=1 |
| 457 | IHNT2(9)=IHNT2(14) |
| 458 | IHNT2(10)=IHNT2(15) |
| 459 | DO 105 I05=1,5 |
| 460 | HINT1(20+I05)=HINT1(40+I05) |
| 461 | HINT1(30+I05)=HINT1(50+I05) |
| 462 | 105 CONTINUE |
| 463 | JOUT=1 |
| 464 | IF(IHPR2(8).EQ.0) GO TO 160 |
| 465 | RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2 |
| 466 | & /REAL(IHNT2(1))**0.6666667,1.0) |
| 467 | RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2 |
| 468 | & /REAL(IHNT2(3))**0.6666667,1.0) |
| 469 | APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0) |
| 470 | & *SQRT(1.0-RRB1) |
| 471 | APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0) |
| 472 | & *SQRT(1.0-RRB2) |
| 473 | HINT1(65)=HINT1(61)-APHX1*HINT1(62) |
| 474 | & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64) |
| 475 | TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31) |
| 476 | NJET=-1 |
| 477 | C ********subtract the trigger jet from total number |
| 478 | C of jet production to be done since it has |
| 479 | C already been produced here |
| 480 | XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG))) |
| 481 | 106 NJET=NJET+1 |
| 482 | XR1=XR1-ALOG(RLU_HIJING(0)) |
| 483 | IF(XR1.LT.TTRIG) GO TO 106 |
| 484 | XR=0.0 |
| 485 | 107 NJET=NJET+1 |
| 486 | XR=XR-ALOG(RLU_HIJING(0)) |
| 487 | IF(XR.LT.TT-TTRIG) GO TO 107 |
| 488 | NJET=NJET-1 |
| 489 | GO TO 112 |
| 490 | ENDIF |
| 491 | C ********create a hard interaction with specified P_T |
| 492 | c when IHPR2(3)>0 |
| 493 | IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110 |
| 494 | C ********create at least one pair of mini jets |
| 495 | C when IHPR2(9)=1 |
| 496 | C |
| 497 | IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)* |
| 498 | & (1.0-EXP(-TTS))) GO TO 160 |
| 499 | C ********this is the probability for no jet production |
| 500 | 110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT))) |
| 501 | 111 NJET=NJET+1 |
| 502 | XR=XR-ALOG(RLU_HIJING(0)) |
| 503 | IF(XR.LT.TT) GO TO 111 |
| 504 | 112 NJET=MIN(NJET,IHPR2(8)) |
| 505 | IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8)) |
| 506 | C ******** Determine number of mini jet production |
| 507 | C |
| 508 | DO 150 I_JET=1,NJET |
| 509 | CALL JETINI(JP,JT,0) |
| 510 | CALL HIJHRD(JP,JT,JOUT,JFLG,1) |
| 511 | C ********JFLG=1 jets valence quarks, JFLG=2 with |
| 512 | C gluon jet, JFLG=3 with q-qbar prod for |
| 513 | C (JP,JT). If JFLG=0 jets can not be produced |
| 514 | C this time. If JFLG=-1, error occured abandon |
| 515 | C this event. JOUT is the total hard scat for |
| 516 | C (JP,JT) up to now. |
| 517 | IF(JFLG.EQ.0) GO TO 160 |
| 518 | IF(JFLG.LT.0) THEN |
| 519 | IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD' |
| 520 | GO TO 50 |
| 521 | ENDIF |
| 522 | JOUT=JOUT+1 |
| 523 | IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1 |
| 524 | IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1 |
| 525 | IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND. |
| 526 | & JFLG.GE.3) IASG(NSG,3)=1 |
| 527 | C ******** jet with PT>HIPR1(11) will be quenched |
| 528 | 150 CONTINUE |
| 529 | 160 CONTINUE |
| 530 | CALL HIJSFT(JP,JT,JOUT,IERROR) |
| 531 | IF(IERROR.NE.0) THEN |
| 532 | IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT' |
| 533 | GO TO 50 |
| 534 | ENDIF |
| 535 | C |
| 536 | C ********conduct soft scattering between JP and JT |
| 537 | JATT=JATT+JOUT |
| 538 | |
| 539 | 200 CONTINUE |
| 540 | |
| 8e529af3 |
541 | c |
| 542 | c************************** |
| 543 | c |
| 544 | DO 201 JP=1,IHNT2(1) |
| 545 | IF(NFP(JP,5).GT.2) THEN |
| 546 | NINP=NINP+1 |
| 547 | ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN |
| 548 | NELP=NELP+1 |
| 549 | ENDIF |
| 550 | 201 continue |
| 551 | DO 202 JT=1,IHNT2(3) |
| 552 | IF(NFT(JT,5).GT.2) THEN |
| 553 | NINT=NINT+1 |
| 554 | ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN |
| 555 | NELT=NELT+1 |
| 556 | ENDIF |
| 557 | 202 continue |
| 558 | c |
| 559 | c******************************* |
| 560 | |
| e74335a4 |
561 | C********perform jet quenching for jets with PT>HIPR1(11)********** |
| 562 | |
| 563 | IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND. |
| 564 | & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN |
| 565 | DO 271 I=1,IHNT2(1) |
| 566 | IF(NFP(I,7).EQ.1) CALL QUENCH(I,1) |
| 567 | 271 CONTINUE |
| 568 | DO 272 I=1,IHNT2(3) |
| 569 | IF(NFT(I,7).EQ.1) CALL QUENCH(I,2) |
| 570 | 272 CONTINUE |
| 571 | DO 273 ISG=1,NSG |
| 572 | IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3) |
| 573 | 273 CONTINUE |
| 574 | ENDIF |
| 575 | C |
| 576 | C**************fragment all the string systems in the following***** |
| 577 | C |
| 578 | C********N_ST is where particle information starts |
| 579 | C********N_STR+1 is the number of strings in fragmentation |
| 580 | C********the number of strings before a line is stored in K(I,4) |
| 581 | C********IDSTR is id number of the string system (91,92 or 93) |
| 582 | C |
| cb3642a4 |
583 | GAMMA0=HINT1(1)/2.0/0.93827 |
| 584 | BETA0=sqrt(GAMMA0**2-1.0)/GAMMA0 |
| 585 | c |
| e74335a4 |
586 | IF(IHPR2(20).NE.0) THEN |
| 587 | DO 360 ISG=1,NSG |
| 588 | CALL HIJFRG(ISG,3,IERROR) |
| 589 | IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN |
| 590 | MSTU(24)=0 |
| 591 | MSTU(28)=0 |
| 592 | IF(IHPR2(10).NE.0) THEN |
| 593 | call LULIST_HIJING(1) |
| 594 | WRITE(6,*) 'error occured, repeat the event' |
| 595 | ENDIF |
| 596 | GO TO 50 |
| 597 | ENDIF |
| 598 | C ********Check errors |
| 599 | C |
| 600 | N_ST=1 |
| 601 | IDSTR=92 |
| 602 | IF(IHPR2(21).EQ.0) THEN |
| 603 | CALL LUEDIT_HIJING(2) |
| 604 | ELSE |
| 605 | 351 N_ST=N_ST+1 |
| 606 | IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351 |
| 607 | IDSTR=K(N_ST,2) |
| 608 | N_ST=N_ST+1 |
| 609 | ENDIF |
| 610 | C |
| 611 | IF(FRAME.EQ.'LAB') THEN |
| 612 | CALL HIBOOST |
| 613 | ENDIF |
| 614 | C ******** boost back to lab frame(if it was in) |
| 615 | C |
| 616 | N_STR=0 |
| 617 | DO 360 I=N_ST,N |
| 618 | IF(K(I,2).EQ.IDSTR) THEN |
| 619 | N_STR=N_STR+1 |
| 620 | GO TO 360 |
| 621 | ENDIF |
| 622 | K(I,4)=N_STR |
| 623 | NATT=NATT+1 |
| 624 | KATT(NATT,1)=K(I,2) |
| 625 | KATT(NATT,2)=20 |
| 626 | KATT(NATT,4)=K(I,1) |
| 627 | IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN |
| 628 | KATT(NATT,3)=0 |
| 629 | ELSE |
| 630 | KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4) |
| 631 | ENDIF |
| 632 | C ****** identify the mother particle |
| 633 | PATT(NATT,1)=P(I,1) |
| 634 | PATT(NATT,2)=P(I,2) |
| 635 | PATT(NATT,3)=P(I,3) |
| 636 | PATT(NATT,4)=P(I,4) |
| 637 | EATT=EATT+P(I,4) |
| cb3642a4 |
638 | VATT01=0.5*(YP(1,IASG(ISG,1))+YT(1,IASG(ISG,2))) |
| 639 | VATT02=0.5*(YP(2,IASG(ISG,1))+YT(2,IASG(ISG,2))) |
| 640 | VATT03=0.5*(YP(3,IASG(ISG,1))+YT(3,IASG(ISG,2))) |
| 641 | & /GAMMA0 |
| 642 | VATT04=-0.5*(YP(3,IASG(ISG,1)) |
| 643 | & -YT(3,IASG(ISG,2)))/GAMMA0/BETA0 |
| 644 | RARB=1.12*(IHNT2(1)**0.33333+IHNT2(3)**0.33333) |
| 645 | V3MIN1=RARB/GAMMA0 |
| 646 | V3MIN2=1.0/MAX(1.0,5.08*ABS(PATT(I,3))) |
| 647 | VATT_MIN=MAX(V3MIN1,V3MIN2) |
| 648 | VATT03=VATT03+(0.5-RLU(0))*VATT_MIN |
| 649 | amt2=P(I,1)**2+P(I,2)**2+P(I,5)**2 |
| 650 | IF(amt2.GT.0.0) THEN |
| 651 | tauf=0.2*2.0*P(I,3)/amt2 |
| 652 | VATT(NATT,1)=VATT01 |
| 653 | VATT(NATT,2)=VATT02 |
| 654 | VATT(NATT,3)=VATT03+tauf |
| 655 | VATT(NATT,4)=abs(VATT(NATT,3)) |
| 656 | ELSE |
| 657 | VATT(NATT,4)=abs(VATT03) |
| 658 | VATT(NATT,1)=VATT01 |
| 659 | VATT(NATT,2)=VATT02 |
| 660 | VATT(NATT,3)=VATT03 |
| 661 | ENDIF |
| e74335a4 |
662 | 360 CONTINUE |
| 663 | C ********Fragment the q-qbar jets systems ***** |
| 664 | C |
| 665 | JTP(1)=IHNT2(1) |
| 666 | JTP(2)=IHNT2(3) |
| 667 | DO 400 NTP=1,2 |
| 668 | DO 400 J_JTP=1,JTP(NTP) |
| 669 | CALL HIJFRG(J_JTP,NTP,IERROR) |
| 670 | IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN |
| 671 | MSTU(24)=0 |
| 672 | MSTU(28)=0 |
| 673 | IF(IHPR2(10).NE.0) THEN |
| 674 | call LULIST_HIJING(1) |
| 675 | WRITE(6,*) 'error occured, repeat the event' |
| 676 | ENDIF |
| 677 | GO TO 50 |
| 678 | ENDIF |
| 679 | C ********check errors |
| 680 | C |
| 681 | N_ST=1 |
| 682 | IDSTR=92 |
| 46f33c91 |
683 | |
| e74335a4 |
684 | IF(IHPR2(21).EQ.0) THEN |
| 685 | CALL LUEDIT_HIJING(2) |
| cb3642a4 |
686 | ELSE |
| e74335a4 |
687 | 381 N_ST=N_ST+1 |
| 688 | IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381 |
| 689 | IDSTR=K(N_ST,2) |
| 690 | N_ST=N_ST+1 |
| 46f33c91 |
691 | ENDIF |
| e74335a4 |
692 | IF(FRAME.EQ.'LAB') THEN |
| 693 | CALL HIBOOST |
| 694 | ENDIF |
| 695 | C ******** boost back to lab frame(if it was in) |
| 696 | C |
| cb3642a4 |
697 | NFTP=NFP(J_JTP,5) |
| 698 | IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5) |
| e74335a4 |
699 | N_STR=0 |
| 700 | DO 390 I=N_ST,N |
| 701 | IF(K(I,2).EQ.IDSTR) THEN |
| 702 | N_STR=N_STR+1 |
| 703 | GO TO 390 |
| 704 | ENDIF |
| 705 | K(I,4)=N_STR |
| 706 | NATT=NATT+1 |
| 707 | KATT(NATT,1)=K(I,2) |
| 708 | KATT(NATT,2)=NFTP |
| 709 | KATT(NATT,4)=K(I,1) |
| cb3642a4 |
710 | IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN |
| e74335a4 |
711 | KATT(NATT,3)=0 |
| cb3642a4 |
712 | ELSE |
| 713 | KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4) |
| e74335a4 |
714 | ENDIF |
| 715 | C ****** identify the mother particle |
| 716 | PATT(NATT,1)=P(I,1) |
| 717 | PATT(NATT,2)=P(I,2) |
| 718 | PATT(NATT,3)=P(I,3) |
| 719 | PATT(NATT,4)=P(I,4) |
| 720 | EATT=EATT+P(I,4) |
| cb3642a4 |
721 | C ******* the following calculate the formation time and position |
| 722 | C of the produced hadrons |
| 723 | IF(NTP.EQ.1) THEN |
| 724 | VATT01=YP(1,J_JTP) |
| 725 | VATT02=YP(2,J_JTP) |
| 726 | VATT03=YP(3,J_JTP)/GAMMA0 |
| 727 | ELSE IF(NTP.EQ.2) THEN |
| 728 | VATT01=YT(1,J_JTP) |
| 729 | VATT02=YT(2,J_JTP) |
| 730 | VATT03=YT(3,J_JTP)/GAMMA0 |
| 731 | ENDIF |
| 732 | RARB=1.12*(IHNT2(1)**0.33333+IHNT2(3)**0.33333) |
| 733 | V3MIN1=RARB/GAMMA0 |
| 734 | V3MIN2=1.0/MAX(1.0,5.08*ABS(PATT(I,3))) |
| 735 | VATT_MIN=MAX(V3MIN1,V3MIN2) |
| 736 | VATT03=VATT03+(0.5-RLU(0))*VATT_MIN |
| 737 | amt2=P(I,1)**2+P(I,2)**2+P(I,5)**2 |
| 738 | IF(amt2.GT.0.0) THEN |
| 739 | tauf=0.2*2.0*P(I,3)/amt2 |
| 740 | VATT(NATT,1)=VATT01 |
| 741 | VATT(NATT,2)=VATT02 |
| 742 | VATT(NATT,3)=VATT03+tauf |
| 743 | VATT(NATT,4)=abs(VATT(NATT,3)) |
| 744 | ELSE |
| 745 | VATT(NATT,4)=abs(VATT03) |
| 746 | VATT(NATT,1)=VATT01 |
| 747 | VATT(NATT,2)=VATT02 |
| 748 | VATT(NATT,3)=VATT03 |
| 749 | ENDIF |
| e74335a4 |
750 | 390 CONTINUE |
| 751 | 400 CONTINUE |
| 752 | C ********Fragment the q-qq related string systems |
| 753 | ENDIF |
| 754 | |
| 755 | DO 450 I=1,NDR |
| 756 | NATT=NATT+1 |
| 757 | KATT(NATT,1)=KFDR(I) |
| 758 | KATT(NATT,2)=40 |
| 759 | KATT(NATT,3)=0 |
| 760 | PATT(NATT,1)=PDR(I,1) |
| 761 | PATT(NATT,2)=PDR(I,2) |
| 762 | PATT(NATT,3)=PDR(I,3) |
| 763 | PATT(NATT,4)=PDR(I,4) |
| 764 | EATT=EATT+PDR(I,4) |
| cb3642a4 |
765 | amt2=P(I,1)**2+P(I,2)**2+P(I,5)**2 |
| 766 | IF(amt2.GT.0.0) THEN |
| 767 | tauf=0.2*2.0*P(I,3)/amt2 |
| 768 | VATT(NATT,1)=0.0 |
| 769 | VATT(NATT,2)=0.0 |
| 770 | VATT(NATT,3)=tauf |
| 771 | VATT(NATT,4)=abs(VATT(NATT,3)) |
| 772 | ELSE |
| 773 | VATT(NATT,4)=0.0 |
| 774 | VATT(NATT,1)=0.0 |
| 775 | VATT(NATT,2)=0.0 |
| 776 | VATT(NATT,3)=0.0 |
| 777 | ENDIF |
| e74335a4 |
778 | 450 CONTINUE |
| 779 | C ********store the direct-produced particles |
| 780 | C |
| 781 | DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0 |
| 782 | IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0 |
| 783 | & .AND.IHPR2(21).EQ.0) THEN |
| 784 | IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the |
| 785 | & event' |
| 786 | C call LULIST_HIJING(1) |
| 787 | GO TO 50 |
| 788 | ENDIF |
| 46f33c91 |
789 | |
| e74335a4 |
790 | RETURN |
| 791 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame