| 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), |
| 7dd84b41 |
157 | & IPCOL(90000),ITCOL(90000),SCIP2(300,300) |
| 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 |
| 7dd84b41 |
345 | SCIP2(JP,JT)=R2 |
| e74335a4 |
346 | C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb |
| 347 | RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2) |
| 348 | & /1.2**2/REAL(IHNT2(1))**0.6666667,1.0) |
| 349 | RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2) |
| 350 | & /1.2**2/REAL(IHNT2(3))**0.6666667,1.0) |
| 351 | APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0) |
| 352 | & *SQRT(1.0-RRB1) |
| 353 | APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0) |
| 354 | & *SQRT(1.0-RRB2) |
| 355 | HINT1(18)=HINT1(14)-APHX1*HINT1(15) |
| 356 | & -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17) |
| 357 | IF(IHPR2(14).EQ.0.OR. |
| 358 | & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN |
| 359 | GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31)) |
| 360 | RANTOT=RLU_HIJING(0) |
| 361 | IF(RANTOT.GT.GS) GO TO 70 |
| 362 | GO TO 65 |
| 363 | ENDIF |
| 364 | GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18)) |
| 365 | & /HIPR1(31)/2.0*ROMG(0.0))) |
| 366 | R2=R2/GSTOT_0 |
| 367 | GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2)) |
| 368 | GSTOT=2.0*(1.0-SQRT(1.0-GS)) |
| 369 | RANTOT=RLU_HIJING(0)*GSTOT_0 |
| 370 | IF(RANTOT.GT.GSTOT) GO TO 70 |
| 371 | IF(RANTOT.GT.GS) THEN |
| 372 | CALL HIJCSC(JP,JT) |
| 373 | GO TO 70 |
| 374 | C ********perform elastic collisions |
| 375 | ENDIF |
| 376 | 65 SCIP(JP,JT)=R2 |
| 377 | RNIP(JP,JT)=RANTOT |
| 378 | SJIP(JP,JT)=HINT1(18) |
| 379 | NCOLT=NCOLT+1 |
| 7dd84b41 |
380 | if (R2.GT.2.D0) THEN |
| 381 | write (8,*) R2 |
| 382 | ENDIF |
| e74335a4 |
383 | IPCOL(NCOLT)=JP |
| 384 | ITCOL(NCOLT)=JT |
| 385 | 70 CONTINUE |
| db0ed1b7 |
386 | |
| db0ed1b7 |
387 | c *** cl glauber *** |
| db0ed1b7 |
388 | npart=0 |
| 389 | xmeana=0D0 |
| 390 | ymeana=0D0 |
| 391 | xmeanb=0D0 |
| 392 | ymeanb=0D0 |
| 393 | xmeanp=0D0 |
| 394 | ymeanp=0D0 |
| 395 | xm2=0D0 |
| 396 | ym2=0D0 |
| 397 | xym=0D0 |
| 398 | |
| fea3949c |
399 | IF(NCOLT>0) THEN |
| 400 | DO 1110 JP=1,IHNT2(1) |
| 401 | YP(1,JP)=YP(1,JP)+BBX |
| 402 | YP(2,JP)=YP(2,JP)+BBY |
| 403 | xmeana=xmeana+YP(1,JP) |
| 404 | ymeana=ymeana+YP(2,JP) |
| 405 | DO 1120 JT=1,IHNT2(3) |
| 7dd84b41 |
406 | IF(SCIP2(JP,JT).LT.2.0D0) THEN |
| fea3949c |
407 | npart=npart+1 |
| 408 | xmeanp=xmeanp+YP(1,JP) |
| 409 | ymeanp=ymeanp+YP(2,JP) |
| 410 | xm2=xm2+YP(1,JP)*YP(1,JP) |
| 411 | ym2=ym2+YP(2,JP)*YP(2,JP) |
| 412 | xym=xym+YP(1,JP)*YP(2,JP) |
| 413 | goto 1110 |
| 414 | end if |
| 415 | 1120 continue |
| 416 | 1110 continue |
| db0ed1b7 |
417 | |
| fea3949c |
418 | DO 1130 JT=1,IHNT2(3) |
| 419 | xmeanb=xmeanb+YT(1,JT) |
| 420 | ymeanb=ymeanb+YT(2,JT) |
| 421 | DO 1140 JP=1,IHNT2(1) |
| 7dd84b41 |
422 | IF(SCIP2(JP,JT).LT.2.0D0) THEN |
| fea3949c |
423 | npart=npart+1 |
| 424 | xmeanp=xmeanp+YT(1,JT) |
| 425 | ymeanp=ymeanp+YT(2,JT) |
| 426 | xm2=xm2+YT(1,JT)*YT(1,JT) |
| 427 | ym2=ym2+YT(2,JT)*YT(2,JT) |
| 428 | xym=xym+YT(1,JT)*YT(2,JT) |
| 429 | goto 1130 |
| 430 | end if |
| 431 | 1140 continue |
| 432 | 1130 continue |
| db0ed1b7 |
433 | |
| 8832289a |
434 | IF (npart.GT.0) THEN |
| 435 | xmeana=xmeana/IHNT2(1) |
| 436 | ymeana=ymeana/IHNT2(1) |
| 437 | xmeanb=xmeanb/IHNT2(3) |
| 438 | ymeanb=ymeanb/IHNT2(3) |
| 439 | xmeanp=xmeanp/npart |
| 440 | ymeanp=ymeanp/npart |
| 441 | xm2=xm2/npart |
| 442 | ym2=ym2/npart |
| 443 | xym=xym/npart |
| db0ed1b7 |
444 | |
| 8832289a |
445 | sx2=xm2-xmeanp*xmeanp |
| 446 | sy2=ym2-ymeanp*ymeanp |
| 447 | sxy=xym-xmeanp*ymeanp |
| fea3949c |
448 | |
| 8832289a |
449 | delx=xmeanb-xmeana |
| 450 | dely=ymeanb-ymeana |
| 451 | dtmp=delx**2+dely**2 |
| 452 | bbtrue=sqrt(dtmp) |
| 453 | dnumt=(sy2-sx2)*(delx**2-dely**2)-4D0*sxy*delx*dely |
| 454 | ddent=(sy2+sx2)*bbtrue**2 |
| 455 | eccrp=dnumt/ddent |
| 456 | dtmp=(sy2-sx2)*(sy2-sx2)+4D0*sxy*sxy |
| 457 | eccpart=sqrt(dtmp)/(sx2+sy2) |
| 458 | eccmc=(sy2-sx2)/(sy2+sx2) |
| 459 | write(*,*),'HOUT: ',bb,' ',bbtrue,' ',ncolt,' ',npart, |
| 460 | 1 ' ',eccrp,' ',eccpart |
| 461 | end if |
| fea3949c |
462 | end if |
| db0ed1b7 |
463 | |
| e74335a4 |
464 | C ********total number interactions proj and targ has |
| 465 | C suffered |
| 7dd84b41 |
466 | |
| e74335a4 |
467 | IF(NCOLT.EQ.0) THEN |
| 468 | NLOP=NLOP+1 |
| 469 | IF(NLOP.LE.20.OR. |
| 470 | & (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60 |
| 46f33c91 |
471 | |
| 472 | |
| e74335a4 |
473 | RETURN |
| 474 | ENDIF |
| 475 | C ********At large impact parameter, there maybe no |
| 476 | C interaction at all. For NN collision |
| 477 | C repeat the event until interaction happens |
| 478 | C |
| 479 | IF(IHPR2(3).NE.0) THEN |
| 480 | NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5) |
| 481 | NHARD=MIN(NHARD,NCOLT) |
| 482 | JPHARD=IPCOL(NHARD) |
| 483 | JTHARD=ITCOL(NHARD) |
| 484 | ENDIF |
| 485 | C |
| 486 | IF(IHPR2(9).EQ.1) THEN |
| 487 | NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5) |
| 488 | NMINI=MIN(NMINI,NCOLT) |
| 489 | JPMINI=IPCOL(NMINI) |
| 490 | JTMINI=ITCOL(NMINI) |
| 491 | ENDIF |
| 492 | C ********Specifying the location of the hard and |
| 493 | C minijet if they are enforced by user |
| 494 | C |
| 495 | DO 200 JP=1,IHNT2(1) |
| 496 | DO 200 JT=1,IHNT2(3) |
| 497 | IF(SCIP(JP,JT).EQ.-1.0) GO TO 200 |
| 498 | NFP(JP,11)=NFP(JP,11)+1 |
| 499 | NFT(JT,11)=NFT(JT,11)+1 |
| 500 | IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN |
| 501 | NP=NP+1 |
| 502 | N01=N01+1 |
| 503 | ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN |
| 504 | NT=NT+1 |
| 505 | N10=N10+1 |
| 506 | ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN |
| 507 | NP=NP+1 |
| 508 | NT=NT+1 |
| 509 | N0=N0+1 |
| 510 | ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN |
| 511 | N11=N11+1 |
| 512 | ENDIF |
| 513 | JOUT=0 |
| 514 | NFP(JP,10)=0 |
| 515 | NFT(JT,10)=0 |
| 516 | C***************************************************************** |
| 517 | IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160 |
| 518 | C ********When IHPR2(8)=0 no jets are produced |
| 519 | IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160 |
| 520 | C ********jets can not be produced for (JP,JT) |
| 521 | C because not enough energy avaible for |
| 522 | C JP or JT |
| 523 | R2=SCIP(JP,JT) |
| 524 | HINT1(18)=SJIP(JP,JT) |
| 525 | TT=ROMG(R2)*HINT1(18)/HIPR1(31) |
| 526 | TTS=HIPR1(30)*ROMG(R2)/HIPR1(31) |
| 527 | NJET=0 |
| 528 | IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN |
| 529 | CALL JETINI(JP,JT,1) |
| 530 | CALL HIJHRD(JP,JT,0,JFLG,0) |
| 531 | HINT1(26)=HINT1(47) |
| 532 | HINT1(27)=HINT1(48) |
| 533 | HINT1(28)=HINT1(49) |
| 534 | HINT1(29)=HINT1(50) |
| 535 | HINT1(36)=HINT1(67) |
| 536 | HINT1(37)=HINT1(68) |
| 537 | HINT1(38)=HINT1(69) |
| 538 | HINT1(39)=HINT1(70) |
| 539 | C |
| 540 | IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1 |
| 541 | IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1 |
| 542 | IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND. |
| 543 | & JFLG.GE.3) IASG(NSG,3)=1 |
| 544 | IHNT2(9)=IHNT2(14) |
| 545 | IHNT2(10)=IHNT2(15) |
| 546 | DO 105 I05=1,5 |
| 547 | HINT1(20+I05)=HINT1(40+I05) |
| 548 | HINT1(30+I05)=HINT1(50+I05) |
| 549 | 105 CONTINUE |
| 550 | JOUT=1 |
| 551 | IF(IHPR2(8).EQ.0) GO TO 160 |
| 552 | RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2 |
| 553 | & /REAL(IHNT2(1))**0.6666667,1.0) |
| 554 | RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2 |
| 555 | & /REAL(IHNT2(3))**0.6666667,1.0) |
| 556 | APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0) |
| 557 | & *SQRT(1.0-RRB1) |
| 558 | APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0) |
| 559 | & *SQRT(1.0-RRB2) |
| 560 | HINT1(65)=HINT1(61)-APHX1*HINT1(62) |
| 561 | & -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64) |
| 562 | TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31) |
| 563 | NJET=-1 |
| 564 | C ********subtract the trigger jet from total number |
| 565 | C of jet production to be done since it has |
| 566 | C already been produced here |
| 567 | XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG))) |
| 568 | 106 NJET=NJET+1 |
| 569 | XR1=XR1-ALOG(RLU_HIJING(0)) |
| 570 | IF(XR1.LT.TTRIG) GO TO 106 |
| 571 | XR=0.0 |
| 572 | 107 NJET=NJET+1 |
| 573 | XR=XR-ALOG(RLU_HIJING(0)) |
| 574 | IF(XR.LT.TT-TTRIG) GO TO 107 |
| 575 | NJET=NJET-1 |
| 576 | GO TO 112 |
| 577 | ENDIF |
| 578 | C ********create a hard interaction with specified P_T |
| 579 | c when IHPR2(3)>0 |
| 580 | IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110 |
| 581 | C ********create at least one pair of mini jets |
| 582 | C when IHPR2(9)=1 |
| 583 | C |
| 584 | IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)* |
| 585 | & (1.0-EXP(-TTS))) GO TO 160 |
| 586 | C ********this is the probability for no jet production |
| 587 | 110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT))) |
| 588 | 111 NJET=NJET+1 |
| 589 | XR=XR-ALOG(RLU_HIJING(0)) |
| 590 | IF(XR.LT.TT) GO TO 111 |
| 591 | 112 NJET=MIN(NJET,IHPR2(8)) |
| 592 | IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8)) |
| 593 | C ******** Determine number of mini jet production |
| 594 | C |
| 595 | DO 150 I_JET=1,NJET |
| 596 | CALL JETINI(JP,JT,0) |
| 597 | CALL HIJHRD(JP,JT,JOUT,JFLG,1) |
| 598 | C ********JFLG=1 jets valence quarks, JFLG=2 with |
| 599 | C gluon jet, JFLG=3 with q-qbar prod for |
| 600 | C (JP,JT). If JFLG=0 jets can not be produced |
| 601 | C this time. If JFLG=-1, error occured abandon |
| 602 | C this event. JOUT is the total hard scat for |
| 603 | C (JP,JT) up to now. |
| 604 | IF(JFLG.EQ.0) GO TO 160 |
| 605 | IF(JFLG.LT.0) THEN |
| 606 | IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD' |
| 607 | GO TO 50 |
| 608 | ENDIF |
| 609 | JOUT=JOUT+1 |
| 610 | IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1 |
| 611 | IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1 |
| 612 | IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND. |
| 613 | & JFLG.GE.3) IASG(NSG,3)=1 |
| 614 | C ******** jet with PT>HIPR1(11) will be quenched |
| 615 | 150 CONTINUE |
| 616 | 160 CONTINUE |
| 617 | CALL HIJSFT(JP,JT,JOUT,IERROR) |
| 618 | IF(IERROR.NE.0) THEN |
| 619 | IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT' |
| 620 | GO TO 50 |
| 621 | ENDIF |
| 622 | C |
| 623 | C ********conduct soft scattering between JP and JT |
| 624 | JATT=JATT+JOUT |
| 625 | |
| 626 | 200 CONTINUE |
| 627 | |
| 8e529af3 |
628 | c |
| 629 | c************************** |
| 630 | c |
| 631 | DO 201 JP=1,IHNT2(1) |
| 65a17253 |
632 | c write(6,*) JP, NFP(JP,3), NFP(JP,4), NFP(JP,5) |
| 8e529af3 |
633 | IF(NFP(JP,5).GT.2) THEN |
| 634 | NINP=NINP+1 |
| 635 | ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN |
| 636 | NELP=NELP+1 |
| 637 | ENDIF |
| 65a17253 |
638 | |
| 639 | IF(NFP(JP,5).LE.2) THEN |
| 640 | IF (NFP(JP,3) .EQ. 2212) THEN |
| 641 | NPSPECP = NPSPECP + 1 |
| 642 | ELSE IF (NFP(JP,3) .EQ. 2112) THEN |
| 643 | NNSPECP = NNSPECP + 1 |
| 644 | ENDIF |
| 645 | ENDIF |
| 8e529af3 |
646 | 201 continue |
| 647 | DO 202 JT=1,IHNT2(3) |
| 648 | IF(NFT(JT,5).GT.2) THEN |
| 649 | NINT=NINT+1 |
| 650 | ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN |
| 651 | NELT=NELT+1 |
| 652 | ENDIF |
| 65a17253 |
653 | |
| 654 | IF(NFT(JT,5).LE.2) THEN |
| 655 | IF (NFT(JT,3) .EQ. 2212) THEN |
| 656 | NPSPECT = NPSPECT + 1 |
| 657 | ELSE IF (NFT(JT,3) .EQ. 2112) THEN |
| 658 | NNSPECT = NNSPECT + 1 |
| 659 | ENDIF |
| 660 | ENDIF |
| 8e529af3 |
661 | 202 continue |
| 662 | c |
| 663 | c******************************* |
| 664 | |
| e74335a4 |
665 | C********perform jet quenching for jets with PT>HIPR1(11)********** |
| 666 | |
| 667 | IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND. |
| 668 | & IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN |
| 669 | DO 271 I=1,IHNT2(1) |
| 670 | IF(NFP(I,7).EQ.1) CALL QUENCH(I,1) |
| 671 | 271 CONTINUE |
| 672 | DO 272 I=1,IHNT2(3) |
| 673 | IF(NFT(I,7).EQ.1) CALL QUENCH(I,2) |
| 674 | 272 CONTINUE |
| 675 | DO 273 ISG=1,NSG |
| 676 | IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3) |
| 677 | 273 CONTINUE |
| 678 | ENDIF |
| 679 | C |
| 680 | C**************fragment all the string systems in the following***** |
| 681 | C |
| 682 | C********N_ST is where particle information starts |
| 683 | C********N_STR+1 is the number of strings in fragmentation |
| 684 | C********the number of strings before a line is stored in K(I,4) |
| 685 | C********IDSTR is id number of the string system (91,92 or 93) |
| 686 | C |
| cb3642a4 |
687 | c |
| e74335a4 |
688 | IF(IHPR2(20).NE.0) THEN |
| 689 | DO 360 ISG=1,NSG |
| 690 | CALL HIJFRG(ISG,3,IERROR) |
| 691 | IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN |
| 692 | MSTU(24)=0 |
| 693 | MSTU(28)=0 |
| 694 | IF(IHPR2(10).NE.0) THEN |
| 695 | call LULIST_HIJING(1) |
| 696 | WRITE(6,*) 'error occured, repeat the event' |
| 697 | ENDIF |
| 698 | GO TO 50 |
| 699 | ENDIF |
| 700 | C ********Check errors |
| 701 | C |
| 702 | N_ST=1 |
| 703 | IDSTR=92 |
| 704 | IF(IHPR2(21).EQ.0) THEN |
| 705 | CALL LUEDIT_HIJING(2) |
| dca9e44f |
706 | ELSE |
| db0ed1b7 |
707 | 351 N_ST=N_ST+1 |
| e74335a4 |
708 | IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351 |
| db0ed1b7 |
709 | IDSTR=K(N_ST,2) |
| 710 | N_ST=N_ST+1 |
| e74335a4 |
711 | ENDIF |
| 712 | C |
| 713 | IF(FRAME.EQ.'LAB') THEN |
| 714 | CALL HIBOOST |
| 715 | ENDIF |
| 716 | C ******** boost back to lab frame(if it was in) |
| 717 | C |
| 718 | N_STR=0 |
| 719 | DO 360 I=N_ST,N |
| 720 | IF(K(I,2).EQ.IDSTR) THEN |
| 721 | N_STR=N_STR+1 |
| 722 | GO TO 360 |
| 723 | ENDIF |
| 724 | K(I,4)=N_STR |
| 725 | NATT=NATT+1 |
| 726 | KATT(NATT,1)=K(I,2) |
| 727 | KATT(NATT,2)=20 |
| 728 | KATT(NATT,4)=K(I,1) |
| 65a17253 |
729 | IF(K(I,3).EQ.0 ) THEN |
| 730 | KATT(NATT,3)=0 |
| 731 | ELSE |
| 732 | IF(K(K(I,3),2).EQ.IDSTR) THEN |
| 733 | KATT(NATT,3)=0 |
| 734 | ELSE |
| 735 | KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4) |
| 736 | ENDIF |
| 737 | ENDIF |
| e74335a4 |
738 | C ****** identify the mother particle |
| 739 | PATT(NATT,1)=P(I,1) |
| 740 | PATT(NATT,2)=P(I,2) |
| 741 | PATT(NATT,3)=P(I,3) |
| 742 | PATT(NATT,4)=P(I,4) |
| 1a0dd691 |
743 | VATT(NATT,1)=V(I,1) |
| 744 | VATT(NATT,2)=V(I,2) |
| 745 | VATT(NATT,3)=V(I,3) |
| 746 | VATT(NATT,4)=V(I,4) |
| 747 | |
| e74335a4 |
748 | EATT=EATT+P(I,4) |
| 749 | 360 CONTINUE |
| 750 | C ********Fragment the q-qbar jets systems ***** |
| 751 | C |
| 752 | JTP(1)=IHNT2(1) |
| 753 | JTP(2)=IHNT2(3) |
| 754 | DO 400 NTP=1,2 |
| 755 | DO 400 J_JTP=1,JTP(NTP) |
| 756 | CALL HIJFRG(J_JTP,NTP,IERROR) |
| 757 | IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN |
| 758 | MSTU(24)=0 |
| 759 | MSTU(28)=0 |
| 760 | IF(IHPR2(10).NE.0) THEN |
| 761 | call LULIST_HIJING(1) |
| 762 | WRITE(6,*) 'error occured, repeat the event' |
| 763 | ENDIF |
| 764 | GO TO 50 |
| 765 | ENDIF |
| 766 | C ********check errors |
| 767 | C |
| 768 | N_ST=1 |
| 769 | IDSTR=92 |
| 46f33c91 |
770 | |
| 65a17253 |
771 | NFTP=NFP(J_JTP,5) |
| 772 | IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5) |
| 773 | |
| e74335a4 |
774 | IF(IHPR2(21).EQ.0) THEN |
| 775 | CALL LUEDIT_HIJING(2) |
| 65a17253 |
776 | ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN |
| e74335a4 |
777 | 381 N_ST=N_ST+1 |
| 778 | IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381 |
| 779 | IDSTR=K(N_ST,2) |
| 780 | N_ST=N_ST+1 |
| 46f33c91 |
781 | ENDIF |
| e74335a4 |
782 | IF(FRAME.EQ.'LAB') THEN |
| 783 | CALL HIBOOST |
| 784 | ENDIF |
| 785 | C ******** boost back to lab frame(if it was in) |
| 786 | C |
| 65a17253 |
787 | |
| e74335a4 |
788 | N_STR=0 |
| 789 | DO 390 I=N_ST,N |
| 790 | IF(K(I,2).EQ.IDSTR) THEN |
| 791 | N_STR=N_STR+1 |
| 792 | GO TO 390 |
| 793 | ENDIF |
| 794 | K(I,4)=N_STR |
| 795 | NATT=NATT+1 |
| 796 | KATT(NATT,1)=K(I,2) |
| 797 | KATT(NATT,2)=NFTP |
| 798 | KATT(NATT,4)=K(I,1) |
| cb3642a4 |
799 | IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN |
| e74335a4 |
800 | KATT(NATT,3)=0 |
| cb3642a4 |
801 | ELSE |
| 802 | KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4) |
| e74335a4 |
803 | ENDIF |
| 804 | C ****** identify the mother particle |
| 805 | PATT(NATT,1)=P(I,1) |
| 806 | PATT(NATT,2)=P(I,2) |
| 807 | PATT(NATT,3)=P(I,3) |
| 808 | PATT(NATT,4)=P(I,4) |
| 809 | EATT=EATT+P(I,4) |
| 4181e24f |
810 | VATT(NATT,1)=V(I,1) |
| 811 | VATT(NATT,2)=V(I,2) |
| 812 | VATT(NATT,3)=V(I,3) |
| 813 | VATT(NATT,4)=V(I,4) |
| e74335a4 |
814 | 390 CONTINUE |
| 815 | 400 CONTINUE |
| 816 | C ********Fragment the q-qq related string systems |
| 817 | ENDIF |
| 818 | |
| 819 | DO 450 I=1,NDR |
| 820 | NATT=NATT+1 |
| 821 | KATT(NATT,1)=KFDR(I) |
| 822 | KATT(NATT,2)=40 |
| 823 | KATT(NATT,3)=0 |
| 824 | PATT(NATT,1)=PDR(I,1) |
| 825 | PATT(NATT,2)=PDR(I,2) |
| 826 | PATT(NATT,3)=PDR(I,3) |
| 827 | PATT(NATT,4)=PDR(I,4) |
| 4181e24f |
828 | VATT(NATT,1)=V(I,1) |
| 829 | VATT(NATT,2)=V(I,2) |
| 830 | VATT(NATT,3)=V(I,3) |
| 831 | VATT(NATT,4)=V(I,4) |
| e74335a4 |
832 | EATT=EATT+PDR(I,4) |
| 833 | 450 CONTINUE |
| 834 | C ********store the direct-produced particles |
| 835 | C |
| 836 | DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0 |
| 837 | IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0 |
| 838 | & .AND.IHPR2(21).EQ.0) THEN |
| 839 | IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the |
| 840 | & event' |
| 841 | C call LULIST_HIJING(1) |
| 842 | GO TO 50 |
| 843 | ENDIF |
| 46f33c91 |
844 | |
| e74335a4 |
845 | RETURN |
| 846 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame