| 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" |
| 162 | #include "himain1.inc" |
| 163 | #include "himain2.inc" |
| 164 | #include "histrng.inc" |
| 165 | #include "hijjet1.inc" |
| 166 | #include "hijjet2.inc" |
| 167 | #include "hijjet4.inc" |
| 168 | C |
| 169 | #include "lujets_hijing.inc" |
| 170 | #include "ludat1_hijing.inc" |
| 171 | SAVE |
| 172 | |
| 173 | BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35)) |
| 174 | BMIN=MIN(BMIN0,BMAX) |
| 175 | IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN |
| 176 | BMIN=0.0 |
| 177 | BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40)) |
| 178 | ENDIF |
| 179 | C ********HIPR1(31) is in mb =0.1fm**2 |
| 180 | C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS |
| 181 | C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm) |
| 182 | C |
| 183 | YP(1,1)=0.0 |
| 184 | YP(2,1)=0.0 |
| 185 | YP(3,1)=0.0 |
| 186 | IF(IHNT2(1).LE.1) GO TO 14 |
| 187 | DO 10 KP=1,IHNT2(1) |
| 188 | 5 R=HIRND(1) |
| 189 | X=RLU_HIJING(0) |
| 190 | CX=2.0*X-1.0 |
| 191 | SX=SQRT(1.0-CX*CX) |
| 192 | C ********choose theta from uniform cos(theta) distr |
| 193 | PHI=RLU_HIJING(0)*2.0*HIPR1(40) |
| 194 | C ********choose phi form uniform phi distr 0 to 2*pi |
| 195 | YP(1,KP)=R*SX*COS(PHI) |
| 196 | YP(2,KP)=R*SX*SIN(PHI) |
| 197 | YP(3,KP)=R*CX |
| 198 | IF(HIPR1(29).EQ.0.0) GO TO 10 |
| 199 | DO 8 KP2=1,KP-1 |
| 200 | DNBP1=(YP(1,KP)-YP(1,KP2))**2 |
| 201 | DNBP2=(YP(2,KP)-YP(2,KP2))**2 |
| 202 | DNBP3=(YP(3,KP)-YP(3,KP2))**2 |
| 203 | DNBP=DNBP1+DNBP2+DNBP3 |
| 204 | IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5 |
| 205 | C ********two neighbors cannot be closer than |
| 206 | C HIPR1(29) |
| 207 | 8 CONTINUE |
| 208 | 10 CONTINUE |
| 209 | DO 12 I=1,IHNT2(1)-1 |
| 210 | DO 12 J=I+1,IHNT2(1) |
| 211 | IF(YP(3,I).GT.YP(3,J)) GO TO 12 |
| 212 | Y1=YP(1,I) |
| 213 | Y2=YP(2,I) |
| 214 | Y3=YP(3,I) |
| 215 | YP(1,I)=YP(1,J) |
| 216 | YP(2,I)=YP(2,J) |
| 217 | YP(3,I)=YP(3,J) |
| 218 | YP(1,J)=Y1 |
| 219 | YP(2,J)=Y2 |
| 220 | YP(3,J)=Y3 |
| 221 | 12 CONTINUE |
| 222 | C |
| 223 | C****************************** |
| 224 | 14 YT(1,1)=0.0 |
| 225 | YT(2,1)=0.0 |
| 226 | YT(3,1)=0.0 |
| 227 | IF(IHNT2(3).LE.1) GO TO 24 |
| 228 | DO 20 KT=1,IHNT2(3) |
| 229 | 15 R=HIRND(2) |
| 230 | X=RLU_HIJING(0) |
| 231 | CX=2.0*X-1.0 |
| 232 | SX=SQRT(1.0-CX*CX) |
| 233 | C ********choose theta from uniform cos(theta) distr |
| 234 | PHI=RLU_HIJING(0)*2.0*HIPR1(40) |
| 235 | C ********chose phi form uniform phi distr 0 to 2*pi |
| 236 | YT(1,KT)=R*SX*COS(PHI) |
| 237 | YT(2,KT)=R*SX*SIN(PHI) |
| 238 | YT(3,KT)=R*CX |
| 239 | IF(HIPR1(29).EQ.0.0) GO TO 20 |
| 240 | DO 18 KT2=1,KT-1 |
| 241 | DNBT1=(YT(1,KT)-YT(1,KT2))**2 |
| 242 | DNBT2=(YT(2,KT)-YT(2,KT2))**2 |
| 243 | DNBT3=(YT(3,KT)-YT(3,KT2))**2 |
| 244 | DNBT=DNBT1+DNBT2+DNBT3 |
| 245 | IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15 |
| 246 | C ********two neighbors cannot be closer than |
| 247 | C HIPR1(29) |
| 248 | 18 CONTINUE |
| 249 | 20 CONTINUE |
| 250 | DO 22 I=1,IHNT2(3)-1 |
| 251 | DO 22 J=I+1,IHNT2(3) |
| 252 | IF(YT(3,I).LT.YT(3,J)) GO TO 22 |
| 253 | Y1=YT(1,I) |
| 254 | Y2=YT(2,I) |
| 255 | Y3=YT(3,I) |
| 256 | YT(1,I)=YT(1,J) |
| 257 | YT(2,I)=YT(2,J) |
| 258 | YT(3,I)=YT(3,J) |
| 259 | YT(1,J)=Y1 |
| 260 | YT(2,J)=Y2 |
| 261 | YT(3,J)=Y3 |
| 262 | 22 CONTINUE |
| 263 | |
| 264 | C******************** |
| 265 | 24 MISS=-1 |
| 266 | |
| 267 | 50 MISS=MISS+1 |
| 268 | IF(MISS.GT.50) THEN |
| 269 | WRITE(6,*) 'infinite loop happened in HIJING' |
| 270 | STOP |
| 271 | ENDIF |
| 272 | |
| 273 | NATT=0 |
| 274 | JATT=0 |
| 275 | EATT=0.0 |
| 276 | CALL HIJINI |
| 277 | NLOP=0 |
| 278 | C ********Initialize for a new event |
| 279 | 60 NT=0 |
| 280 | NP=0 |
| 281 | N0=0 |
| 282 | N01=0 |
| 283 | N10=0 |
| 284 | N11=0 |
| 285 | NSG=0 |
| 286 | NCOLT=0 |
| 287 | |
| 288 | C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS |
| 289 | C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET |
| 290 | C BY THE ANGLE PHI FOR EACH COLLISION.****************** |
| 291 | C |
| 292 | BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2)) |
| 293 | PHI=2.0*HIPR1(40)*RLU_HIJING(0) |
| 294 | BBX=BB*COS(PHI) |
| 295 | BBY=BB*SIN(PHI) |
| 296 | HINT1(19)=BB |
| 297 | HINT1(20)=PHI |
| 298 | C |
| 299 | DO 70 JP=1,IHNT2(1) |
| 300 | DO 70 JT=1,IHNT2(3) |
| 301 | SCIP(JP,JT)=-1.0 |
| 302 | B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2 |
| 303 | R2=B2*HIPR1(40)/HIPR1(31)/0.1 |
| 304 | C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb |
| 305 | RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2) |
| 306 | & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0) |
| 307 | RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2) |
| 308 | & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0) |
| 309 | APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0) |
| 310 | & *SQRT(1.0-RRB1) |
| 311 | APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0) |
| 312 | & *SQRT(1.0-RRB2) |
| 313 | HINT1(18)=HINT1(14)-APHX1*HINT1(15) |
| 314 | & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17) |
| 315 | IF(IHPR2(14).EQ.0.OR. |
| 316 | & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN |
| 317 | GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31)) |
| 318 | RANTOT=RLU_HIJING(0) |
| 319 | IF(RANTOT.GT.GS) GO TO 70 |
| 320 | GO TO 65 |
| 321 | ENDIF |
| 322 | GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18)) |
| 323 | & /HIPR1(31)/2.0*ROMG(0.0))) |
| 324 | R2=R2/GSTOT_0 |
| 325 | GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2)) |
| 326 | GSTOT=2.0*(1.0-SQRT(1.0-GS)) |
| 327 | RANTOT=RLU_HIJING(0)*GSTOT_0 |
| 328 | IF(RANTOT.GT.GSTOT) GO TO 70 |
| 329 | IF(RANTOT.GT.GS) THEN |
| 330 | CALL HIJCSC(JP,JT) |
| 331 | GO TO 70 |
| 332 | C ********perform elastic collisions |
| 333 | ENDIF |
| 334 | 65 SCIP(JP,JT)=R2 |
| 335 | RNIP(JP,JT)=RANTOT |
| 336 | SJIP(JP,JT)=HINT1(18) |
| 337 | NCOLT=NCOLT+1 |
| 338 | IPCOL(NCOLT)=JP |
| 339 | ITCOL(NCOLT)=JT |
| 340 | 70 CONTINUE |
| 341 | C ********total number interactions proj and targ has |
| 342 | C suffered |
| 343 | IF(NCOLT.EQ.0) THEN |
| 344 | NLOP=NLOP+1 |
| 345 | IF(NLOP.LE.20.OR. |
| 346 | & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60 |
| 46f33c91 |
347 | |
| 348 | |
| e74335a4 |
349 | RETURN |
| 350 | ENDIF |
| 351 | C ********At large impact parameter, there maybe no |
| 352 | C interaction at all. For NN collision |
| 353 | C repeat the event until interaction happens |
| 354 | C |
| 355 | IF(IHPR2(3).NE.0) THEN |
| 356 | NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5) |
| 357 | NHARD=MIN(NHARD,NCOLT) |
| 358 | JPHARD=IPCOL(NHARD) |
| 359 | JTHARD=ITCOL(NHARD) |
| 360 | ENDIF |
| 361 | C |
| 362 | IF(IHPR2(9).EQ.1) THEN |
| 363 | NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5) |
| 364 | NMINI=MIN(NMINI,NCOLT) |
| 365 | JPMINI=IPCOL(NMINI) |
| 366 | JTMINI=ITCOL(NMINI) |
| 367 | ENDIF |
| 368 | C ********Specifying the location of the hard and |
| 369 | C minijet if they are enforced by user |
| 370 | C |
| 371 | DO 200 JP=1,IHNT2(1) |
| 372 | DO 200 JT=1,IHNT2(3) |
| 373 | IF(SCIP(JP,JT).EQ.-1.0) GO TO 200 |
| 374 | NFP(JP,11)=NFP(JP,11)+1 |
| 375 | NFT(JT,11)=NFT(JT,11)+1 |
| 376 | IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN |
| 377 | NP=NP+1 |
| 378 | N01=N01+1 |
| 379 | ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN |
| 380 | NT=NT+1 |
| 381 | N10=N10+1 |
| 382 | ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN |
| 383 | NP=NP+1 |
| 384 | NT=NT+1 |
| 385 | N0=N0+1 |
| 386 | ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN |
| 387 | N11=N11+1 |
| 388 | ENDIF |
| 389 | JOUT=0 |
| 390 | NFP(JP,10)=0 |
| 391 | NFT(JT,10)=0 |
| 392 | C***************************************************************** |
| 393 | IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160 |
| 394 | C ********When IHPR2(8)=0 no jets are produced |
| 395 | IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160 |
| 396 | C ********jets can not be produced for (JP,JT) |
| 397 | C because not enough energy avaible for |
| 398 | C JP or JT |
| 399 | R2=SCIP(JP,JT) |
| 400 | HINT1(18)=SJIP(JP,JT) |
| 401 | TT=ROMG(R2)*HINT1(18)/HIPR1(31) |
| 402 | TTS=HIPR1(30)*ROMG(R2)/HIPR1(31) |
| 403 | NJET=0 |
| 404 | IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN |
| 405 | CALL JETINI(JP,JT,1) |
| 406 | CALL HIJHRD(JP,JT,0,JFLG,0) |
| 407 | HINT1(26)=HINT1(47) |
| 408 | HINT1(27)=HINT1(48) |
| 409 | HINT1(28)=HINT1(49) |
| 410 | HINT1(29)=HINT1(50) |
| 411 | HINT1(36)=HINT1(67) |
| 412 | HINT1(37)=HINT1(68) |
| 413 | HINT1(38)=HINT1(69) |
| 414 | HINT1(39)=HINT1(70) |
| 415 | C |
| 416 | IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1 |
| 417 | IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1 |
| 418 | IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND. |
| 419 | & JFLG.GE.3) IASG(NSG,3)=1 |
| 420 | IHNT2(9)=IHNT2(14) |
| 421 | IHNT2(10)=IHNT2(15) |
| 422 | DO 105 I05=1,5 |
| 423 | HINT1(20+I05)=HINT1(40+I05) |
| 424 | HINT1(30+I05)=HINT1(50+I05) |
| 425 | 105 CONTINUE |
| 426 | JOUT=1 |
| 427 | IF(IHPR2(8).EQ.0) GO TO 160 |
| 428 | RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2 |
| 429 | & /REAL(IHNT2(1))**0.6666667,1.0) |
| 430 | RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2 |
| 431 | & /REAL(IHNT2(3))**0.6666667,1.0) |
| 432 | APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0) |
| 433 | & *SQRT(1.0-RRB1) |
| 434 | APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0) |
| 435 | & *SQRT(1.0-RRB2) |
| 436 | HINT1(65)=HINT1(61)-APHX1*HINT1(62) |
| 437 | & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64) |
| 438 | TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31) |
| 439 | NJET=-1 |
| 440 | C ********subtract the trigger jet from total number |
| 441 | C of jet production to be done since it has |
| 442 | C already been produced here |
| 443 | XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG))) |
| 444 | 106 NJET=NJET+1 |
| 445 | XR1=XR1-ALOG(RLU_HIJING(0)) |
| 446 | IF(XR1.LT.TTRIG) GO TO 106 |
| 447 | XR=0.0 |
| 448 | 107 NJET=NJET+1 |
| 449 | XR=XR-ALOG(RLU_HIJING(0)) |
| 450 | IF(XR.LT.TT-TTRIG) GO TO 107 |
| 451 | NJET=NJET-1 |
| 452 | GO TO 112 |
| 453 | ENDIF |
| 454 | C ********create a hard interaction with specified P_T |
| 455 | c when IHPR2(3)>0 |
| 456 | IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110 |
| 457 | C ********create at least one pair of mini jets |
| 458 | C when IHPR2(9)=1 |
| 459 | C |
| 460 | IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)* |
| 461 | & (1.0-EXP(-TTS))) GO TO 160 |
| 462 | C ********this is the probability for no jet production |
| 463 | 110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT))) |
| 464 | 111 NJET=NJET+1 |
| 465 | XR=XR-ALOG(RLU_HIJING(0)) |
| 466 | IF(XR.LT.TT) GO TO 111 |
| 467 | 112 NJET=MIN(NJET,IHPR2(8)) |
| 468 | IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8)) |
| 469 | C ******** Determine number of mini jet production |
| 470 | C |
| 471 | DO 150 I_JET=1,NJET |
| 472 | CALL JETINI(JP,JT,0) |
| 473 | CALL HIJHRD(JP,JT,JOUT,JFLG,1) |
| 474 | C ********JFLG=1 jets valence quarks, JFLG=2 with |
| 475 | C gluon jet, JFLG=3 with q-qbar prod for |
| 476 | C (JP,JT). If JFLG=0 jets can not be produced |
| 477 | C this time. If JFLG=-1, error occured abandon |
| 478 | C this event. JOUT is the total hard scat for |
| 479 | C (JP,JT) up to now. |
| 480 | IF(JFLG.EQ.0) GO TO 160 |
| 481 | IF(JFLG.LT.0) THEN |
| 482 | IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD' |
| 483 | GO TO 50 |
| 484 | ENDIF |
| 485 | JOUT=JOUT+1 |
| 486 | IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1 |
| 487 | IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1 |
| 488 | IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND. |
| 489 | & JFLG.GE.3) IASG(NSG,3)=1 |
| 490 | C ******** jet with PT>HIPR1(11) will be quenched |
| 491 | 150 CONTINUE |
| 492 | 160 CONTINUE |
| 493 | CALL HIJSFT(JP,JT,JOUT,IERROR) |
| 494 | IF(IERROR.NE.0) THEN |
| 495 | IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT' |
| 496 | GO TO 50 |
| 497 | ENDIF |
| 498 | C |
| 499 | C ********conduct soft scattering between JP and JT |
| 500 | JATT=JATT+JOUT |
| 501 | |
| 502 | 200 CONTINUE |
| 503 | |
| 504 | C********perform jet quenching for jets with PT>HIPR1(11)********** |
| 505 | |
| 506 | IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND. |
| 507 | & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN |
| 508 | DO 271 I=1,IHNT2(1) |
| 509 | IF(NFP(I,7).EQ.1) CALL QUENCH(I,1) |
| 510 | 271 CONTINUE |
| 511 | DO 272 I=1,IHNT2(3) |
| 512 | IF(NFT(I,7).EQ.1) CALL QUENCH(I,2) |
| 513 | 272 CONTINUE |
| 514 | DO 273 ISG=1,NSG |
| 515 | IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3) |
| 516 | 273 CONTINUE |
| 517 | ENDIF |
| 518 | C |
| 519 | C**************fragment all the string systems in the following***** |
| 520 | C |
| 521 | C********N_ST is where particle information starts |
| 522 | C********N_STR+1 is the number of strings in fragmentation |
| 523 | C********the number of strings before a line is stored in K(I,4) |
| 524 | C********IDSTR is id number of the string system (91,92 or 93) |
| 525 | C |
| 526 | IF(IHPR2(20).NE.0) THEN |
| 527 | DO 360 ISG=1,NSG |
| 528 | CALL HIJFRG(ISG,3,IERROR) |
| 529 | IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN |
| 530 | MSTU(24)=0 |
| 531 | MSTU(28)=0 |
| 532 | IF(IHPR2(10).NE.0) THEN |
| 533 | call LULIST_HIJING(1) |
| 534 | WRITE(6,*) 'error occured, repeat the event' |
| 535 | ENDIF |
| 536 | GO TO 50 |
| 537 | ENDIF |
| 538 | C ********Check errors |
| 539 | C |
| 540 | N_ST=1 |
| 541 | IDSTR=92 |
| 542 | IF(IHPR2(21).EQ.0) THEN |
| 543 | CALL LUEDIT_HIJING(2) |
| 544 | ELSE |
| 545 | 351 N_ST=N_ST+1 |
| 546 | IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351 |
| 547 | IDSTR=K(N_ST,2) |
| 548 | N_ST=N_ST+1 |
| 549 | ENDIF |
| 550 | C |
| 551 | IF(FRAME.EQ.'LAB') THEN |
| 552 | CALL HIBOOST |
| 553 | ENDIF |
| 554 | C ******** boost back to lab frame(if it was in) |
| 555 | C |
| 556 | N_STR=0 |
| 557 | DO 360 I=N_ST,N |
| 558 | IF(K(I,2).EQ.IDSTR) THEN |
| 559 | N_STR=N_STR+1 |
| 560 | GO TO 360 |
| 561 | ENDIF |
| 562 | K(I,4)=N_STR |
| 563 | NATT=NATT+1 |
| 564 | KATT(NATT,1)=K(I,2) |
| 565 | KATT(NATT,2)=20 |
| 566 | KATT(NATT,4)=K(I,1) |
| 567 | IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN |
| 568 | KATT(NATT,3)=0 |
| 569 | ELSE |
| 570 | KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4) |
| 571 | ENDIF |
| 572 | C ****** identify the mother particle |
| 573 | PATT(NATT,1)=P(I,1) |
| 574 | PATT(NATT,2)=P(I,2) |
| 575 | PATT(NATT,3)=P(I,3) |
| 576 | PATT(NATT,4)=P(I,4) |
| 2818cab6 |
577 | VATT(NATT,1)=V(I,1) |
| 578 | VATT(NATT,2)=V(I,2) |
| 579 | VATT(NATT,3)=V(I,3) |
| 580 | VATT(NATT,4)=V(I,4) |
| 581 | |
| e74335a4 |
582 | EATT=EATT+P(I,4) |
| 583 | 360 CONTINUE |
| 584 | C ********Fragment the q-qbar jets systems ***** |
| 585 | C |
| 586 | JTP(1)=IHNT2(1) |
| 587 | JTP(2)=IHNT2(3) |
| 588 | DO 400 NTP=1,2 |
| 589 | DO 400 J_JTP=1,JTP(NTP) |
| 590 | CALL HIJFRG(J_JTP,NTP,IERROR) |
| 591 | IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN |
| 592 | MSTU(24)=0 |
| 593 | MSTU(28)=0 |
| 594 | IF(IHPR2(10).NE.0) THEN |
| 595 | call LULIST_HIJING(1) |
| 596 | WRITE(6,*) 'error occured, repeat the event' |
| 597 | ENDIF |
| 598 | GO TO 50 |
| 599 | ENDIF |
| 600 | C ********check errors |
| 601 | C |
| 602 | N_ST=1 |
| 603 | IDSTR=92 |
| 46f33c91 |
604 | |
| 605 | NFTP=NFP(J_JTP,5) |
| 606 | IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5) |
| 607 | |
| e74335a4 |
608 | IF(IHPR2(21).EQ.0) THEN |
| 609 | CALL LUEDIT_HIJING(2) |
| 46f33c91 |
610 | ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN |
| e74335a4 |
611 | 381 N_ST=N_ST+1 |
| 612 | IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381 |
| 613 | IDSTR=K(N_ST,2) |
| 614 | N_ST=N_ST+1 |
| 46f33c91 |
615 | ENDIF |
| e74335a4 |
616 | IF(FRAME.EQ.'LAB') THEN |
| 617 | CALL HIBOOST |
| 618 | ENDIF |
| 619 | C ******** boost back to lab frame(if it was in) |
| 620 | C |
| e74335a4 |
621 | N_STR=0 |
| 622 | DO 390 I=N_ST,N |
| 623 | IF(K(I,2).EQ.IDSTR) THEN |
| 624 | N_STR=N_STR+1 |
| 625 | GO TO 390 |
| 626 | ENDIF |
| 627 | K(I,4)=N_STR |
| 628 | NATT=NATT+1 |
| 629 | KATT(NATT,1)=K(I,2) |
| 630 | KATT(NATT,2)=NFTP |
| 631 | KATT(NATT,4)=K(I,1) |
| 632 | IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN |
| 633 | KATT(NATT,3)=0 |
| 634 | ELSE |
| 635 | KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4) |
| 636 | ENDIF |
| 637 | C ****** identify the mother particle |
| 638 | PATT(NATT,1)=P(I,1) |
| 639 | PATT(NATT,2)=P(I,2) |
| 640 | PATT(NATT,3)=P(I,3) |
| 641 | PATT(NATT,4)=P(I,4) |
| 2818cab6 |
642 | VATT(NATT,1)=V(I,1) |
| 643 | VATT(NATT,2)=V(I,2) |
| 644 | VATT(NATT,3)=V(I,3) |
| 645 | VATT(NATT,4)=V(I,4) |
| 646 | |
| e74335a4 |
647 | EATT=EATT+P(I,4) |
| 648 | 390 CONTINUE |
| 649 | 400 CONTINUE |
| 650 | C ********Fragment the q-qq related string systems |
| 651 | ENDIF |
| 652 | |
| 653 | DO 450 I=1,NDR |
| 654 | NATT=NATT+1 |
| 655 | KATT(NATT,1)=KFDR(I) |
| 656 | KATT(NATT,2)=40 |
| 657 | KATT(NATT,3)=0 |
| 658 | PATT(NATT,1)=PDR(I,1) |
| 659 | PATT(NATT,2)=PDR(I,2) |
| 660 | PATT(NATT,3)=PDR(I,3) |
| 661 | PATT(NATT,4)=PDR(I,4) |
| 2818cab6 |
662 | VATT(NATT,1)=V(I,1) |
| 663 | VATT(NATT,2)=V(I,2) |
| 664 | VATT(NATT,3)=V(I,3) |
| 665 | VATT(NATT,4)=V(I,4) |
| 666 | |
| e74335a4 |
667 | EATT=EATT+PDR(I,4) |
| 668 | 450 CONTINUE |
| 669 | C ********store the direct-produced particles |
| 670 | C |
| 671 | DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0 |
| 672 | IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0 |
| 673 | & .AND.IHPR2(21).EQ.0) THEN |
| 674 | IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the |
| 675 | & event' |
| 676 | C call LULIST_HIJING(1) |
| 677 | GO TO 50 |
| 678 | ENDIF |
| 46f33c91 |
679 | |
| e74335a4 |
680 | RETURN |
| 681 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame