]>
Commit | Line | Data |
---|---|---|
cb220f83 | 1 | *---------------------------------------------------------------------- |
2 | * | |
3 | * Filename : HYDJET1_1.F | |
4 | * | |
5 | * Author : Igor Lokhtin | |
6 | * Version : HYDJET1_1.f | |
7 | * Last revision : 27-MAR-2006 | |
8 | * | |
9 | *====================================================================== | |
10 | * | |
11 | * Description : Fast event generator for AA collisons at the LHC | |
12 | * | |
13 | * Method I.P. Lokhtin and A.M. Snigirev, Eur. Phys. J C 45 (2006) 211 | |
14 | * | |
15 | *====================================================================== | |
16 | ||
17 | SUBROUTINE HYDRO(A,ifb,bmin,bmax,bfix,nh) | |
18 | real hsin,hgauss,hftaa | |
19 | real AW | |
20 | real A,bmin,bmax,bfix | |
21 | integer numjet,numpar | |
22 | integer ifb,nh,np | |
23 | external hsin,hgauss,hftaa,numjet,numpar,hyhard,hipsear | |
24 | external ludata | |
25 | common /lujets/ n,k(150000,5),p(150000,5),v(150000,5) | |
26 | common /hyjets/ nl,kl(150000,5),pl(150000,5),vl(150000,5) | |
27 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
28 | common /hyfpar/ bgen,nbcol,npart,npyt,nhyd | |
29 | common /hyflow/ ytfl,ylfl,fpart | |
30 | common /hyjpar/ nhsel,ptmin,njet | |
31 | save/lujets/,/hyjets/,/hyipar/,/hyfpar/,/hyflow/,/hyjpar/ | |
32 | ||
33 | integer nnhyd, khyd | |
34 | real phyd, vhyd | |
35 | common /hyd/ nnhyd, khyd(150000,5),phyd(150000,5),vhyd(150000,5) | |
36 | save /hyd/ | |
37 | ||
38 | * reset lujets and hyd arrays before event generation | |
39 | ||
40 | n=0 | |
41 | nnhyd=0 | |
42 | do ncl=1,150000 | |
43 | do j=1,5 | |
44 | p(ncl,j)=0. | |
45 | phyd(ncl,j)=0. | |
46 | v(ncl,j)=0. | |
47 | vhyd(ncl,j)=0. | |
48 | k(ncl,j)=0 | |
49 | khyd(ncl,j)=0 | |
50 | enddo | |
51 | end do | |
52 | ||
53 | * set initial beam paramters: atomic weigth and radius in fm | |
54 | AW=A | |
55 | RA=1.15*AW**0.333333 | |
56 | * | |
57 | pi=3.14159 | |
58 | ||
59 | * generate impact parameter of A-A collision | |
60 | ||
61 | if(ifb.eq.0) then | |
62 | if(bfix.lt.0.) then | |
63 | write(6,*) 'Impact parameter less than zero!' | |
64 | bfix=0. | |
65 | end if | |
66 | if (bfix.gt.2.) then | |
67 | write(6,*) 'Impact parameter larger than two nuclear radius!' | |
68 | bfix=2. | |
69 | end if | |
70 | b1=bfix*RA | |
71 | bgen=bfix | |
72 | else | |
73 | if(bmin.lt.0.) then | |
74 | write(6,*) 'Impact parameter less than zero!' | |
75 | bmin=0. | |
76 | end if | |
77 | if(bmax.gt.2.) then | |
78 | write(6,*) 'Impact parameter larger than two nuclear radius!' | |
79 | bmax=2. | |
80 | end if | |
81 | bminh=bmin | |
82 | bmaxh=bmax | |
83 | call hipsear(fmax1,xmin1) | |
84 | fmax=fmax1 | |
85 | xmin=xmin1 | |
86 | 3 bb1=xmin*rlu(0)+bminh | |
87 | ff1=fmax*rlu(0) | |
88 | fb=hsin(bb1) | |
89 | if(ff1.gt.fb) goto 3 | |
90 | b1=bb1*RA | |
91 | bgen=bb1 | |
92 | end if | |
93 | ||
94 | * set flow parameters | |
95 | Tf=0.14 ! freeze-out temperature | |
96 | if (ylfl.lt.0.01.or.ylfl.gt.7.) ylfl=5. | |
97 | etmax=ylfl ! longitudinal flow rapidity | |
98 | if (ytfl.lt.0.01.or.ytfl.gt.3.) ytfl=1. | |
99 | ytmax=ytfl ! transverse flow rapidity | |
100 | ||
101 | * set inelastic NN cross section, mb | |
102 | sigin=58. | |
103 | ||
104 | * calculate number of nucelons-participants | |
105 | bb=bgen*RA ! impact parameter | |
106 | npart=numpar(bb) ! Npart(b) | |
107 | npar0=411 ! Npart(Pb,b=0) | |
108 | ||
109 | * calculate number of binary NN sub-collisions | |
110 | br=max(1.e-10,0.25*bgen*bgen) | |
111 | factor=9.*sigin*AW*AW/(80.*pi*RA*RA) | |
112 | nbcol=int(factor*hftaa(br)) ! Nsub(b) | |
113 | nbco0=1923 ! Nsub(Pb,b=0) | |
114 | ||
115 | * generate total multiplicity in event, np, | |
116 | * fpart - fraction of soft multiplicity proportional to # of participants, | |
117 | * fbcol=1-fpart - fraction of multiplicity proprtional to # of NN subcollisions | |
118 | if(fpart.le.0.or.fpart.gt.1.) fpart=1. | |
119 | fbcol=1.-fpart | |
120 | rnp=nh*(fpart*npart+fbcol*nbcol)/(fpart*npar0+fbcol*nbco0) | |
121 | np=int(rnp) | |
122 | sign=sqrt(rnp) | |
123 | 1 if(nhsel.lt.4.and.np.gt.0) np=max(0,int(hgauss(rnp,sign))) | |
124 | if(np.gt.150000) then | |
125 | write(6,*) 'Warning, soft multiplicity too large!' | |
126 | goto 1 | |
127 | end if | |
128 | ||
129 | * generate hard parton-parton scattering (Q>ptmin) 'njet' times with PYTHIA | |
130 | if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4) | |
131 | > nhsel=0 | |
132 | njet=0 | |
133 | if(nhsel.ne.0) then | |
134 | if(ptmin.lt.5.or.ptmin.gt.500.) ptmin=10. | |
135 | q=ptmin | |
136 | njet=numjet(q) | |
137 | call hyhard | |
138 | end if | |
139 | ||
140 | npyt=n | |
141 | nhyd=np | |
142 | c if(nhsel.lt.3) then | |
143 | c nhyd=max(0,np-npyt) | |
144 | c else | |
145 | c nhyd=0 | |
146 | c np=n | |
147 | c end if | |
148 | if(nhyd.eq.0) goto 4 | |
149 | ||
150 | * generate sort of hadrons (pions, kaons, nucleons) in jetset7* format | |
151 | do ip=npyt+1,npyt+np !cycle on particles | |
152 | yy=49.*rlu(0) | |
153 | if(yy.lt.11.83333333) then | |
154 | kf=211 | |
155 | elseif(yy.lt.23.66666667) then | |
156 | kf=-211 | |
157 | elseif(yy.lt.35.5) then | |
158 | kf=111 | |
159 | elseif(yy.lt.38.375) then | |
160 | kf=321 | |
161 | elseif(yy.lt.41.25) then | |
162 | kf=-321 | |
163 | elseif(yy.lt.44.125) then | |
164 | kf=310 | |
165 | elseif(yy.lt.47.) then | |
166 | kf=130 | |
167 | elseif(yy.lt.47.5) then | |
168 | kf=2212 | |
169 | elseif(yy.lt.48.) then | |
170 | kf=-2212 | |
171 | elseif(yy.lt.48.5) then | |
172 | kf=2112 | |
173 | else | |
174 | kf=-2112 | |
175 | end if | |
176 | n=n+1 | |
177 | k(n,1)=1 | |
178 | k(n,2)=kf | |
179 | do j=3,5 | |
180 | k(n,j)=0 | |
181 | end do | |
182 | ||
183 | do j=1,5 | |
184 | v(n,j)=0. | |
185 | enddo | |
186 | kfa=iabs(kf) | |
187 | p(n,5)=ulmass(kfa) | |
188 | ||
189 | * generate uniform distribution in system of a fluid element | |
190 | 2 ep0=-1.*Tf*(log(max(1.e-10,rlu(0)))+log(max(1.e-10,rlu(0))) | |
191 | > +log(max(1.e-10,rlu(0)))) | |
192 | if(ep0.le.p(n,5)) go to 2 | |
193 | pp0=sqrt(abs(ep0**2-abs(p(n,5)**2))) | |
194 | probt=pp0/ep0 | |
195 | if(rlu(0).gt.probt) go to 2 | |
196 | ctp0=2.*rlu(0)-1. | |
197 | stp0=sqrt(abs(1.-ctp0**2)) | |
198 | php0=2.*pi*rlu(0) | |
199 | ||
200 | * generate coordinates of a fluid element | |
201 | c etaf=etmax*(2.*rlu(0)-1.) ! flat initial eta-spectrum | |
202 | etaf=hgauss(0.,etmax) ! gaussian initial eta-spectrum | |
203 | phif=2.*pi*rlu(0) | |
204 | rm1=sqrt(abs(RA*RA-b1*b1/4.*(sin(phif)**2)))+b1*cos(phif)/2. | |
205 | rm2=sqrt(abs(RA*RA-b1*b1/4.*(sin(phif)**2)))-b1*cos(phif)/2. | |
206 | RF=min(rm1,rm2) | |
207 | rrf=RF*RF*rlu(0) | |
208 | ||
209 | * generate four-velocity of a fluid element | |
210 | vradf=sinh(ytmax) | |
211 | sb=RA*RA*(pi-2.*asin(b1/RA/2.))-b1*sqrt(abs(RA*RA-b1*b1/4.)) | |
212 | reff=sqrt(sqrt(sb/pi)*RA) | |
213 | urf=vradf*sqrt(abs(rrf))/reff ! linear transverse profile | |
214 | c urf=vradf*rrf/reff**2 !square transverse profile | |
215 | utf=cosh(etaf)*sqrt(abs(1.+urf*urf)) | |
216 | uzf=sinh(etaf)*sqrt(abs(1.+urf*urf)) | |
217 | ||
218 | * boost in the laboratory system | |
219 | uipi=pp0*(urf*stp0*cos(phif-php0)+uzf*ctp0) | |
220 | bfac=uipi/(utf+1.)+ep0 | |
221 | p(n,1)=pp0*stp0*sin(php0)+urf*sin(phif)*bfac | |
222 | p(n,2)=pp0*stp0*cos(php0)+urf*cos(phif)*bfac | |
223 | p(n,3)=pp0*ctp0+uzf*bfac | |
224 | p(n,4)=sqrt(p(n,1)**2+p(n,2)**2+p(n,3)**2+p(n,5)**2) | |
225 | ||
226 | end do | |
227 | ||
228 | 4 continue | |
229 | ||
230 | * write(*,*) 'NHYD, NPYT, NTOT',nhyd,npyt,nhyd+npyt | |
231 | ||
232 | * fill array 'hyd' | |
233 | ||
234 | nnhyd = nhyd+npyt | |
235 | ||
236 | do ih=1,n | |
237 | do jh=1,5 | |
238 | phyd(ih,jh)=p(ih,jh) | |
239 | khyd(ih,jh)=k(ih,jh) | |
240 | vhyd(ih,jh)=v(ih,jh) | |
241 | end do | |
242 | end do | |
243 | ||
244 | return | |
245 | end | |
246 | ||
247 | ********************************* HYHARD *************************** | |
248 | SUBROUTINE HYHARD | |
249 | * generate 'njet' number of hard parton-parton scatterings with PYTHIA | |
250 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
251 | IMPLICIT INTEGER(I-N) | |
252 | REAL ptmin,pj,vj,pl,vl,bminh,bmaxh,AW,RA,sigin,bgen | |
253 | INTEGER PYK,PYCHGE,PYCOMP | |
254 | external pydata | |
255 | external pyp,pyr,pyk,pyquen | |
256 | common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5) | |
257 | common /lujets/ nj,kj(150000,5),pj(150000,5),vj(150000,5) | |
258 | common /hyjets/ nl,kl(150000,5),pl(150000,5),vl(150000,5) | |
259 | COMMON /PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) | |
260 | COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) | |
261 | COMMON /PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) | |
262 | COMMON /PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) | |
263 | common /pypars/ mstp(200),parp(200),msti(200),pari(200) | |
264 | common /hyjpar/ nhsel,ptmin,njet | |
265 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
266 | common /hyfpar/ bgen,nbcol,npart,npyt,nhyd | |
267 | save /pyjets/,/lujets/,/hyjets/,/pydat1/,/pydat2/,/pydat3/, | |
268 | + /pysubs/,/hyjpar/,/hyipar/,/hyfpar/ | |
269 | ||
270 | * reset array of partons in 'hyjets' | |
271 | nl=0 | |
272 | do i=1,150000 | |
273 | do j=1,5 | |
274 | pl(i,j)=0. | |
275 | end do | |
276 | do j=1,5 | |
277 | vl(i,j)=0. | |
278 | end do | |
279 | do j=1,5 | |
280 | kl(i,j)=0 | |
281 | end do | |
282 | end do | |
283 | ||
284 | * generate 'njet' PYTHIA events and fill arrays for partons and hadrons | |
285 | if(njet.ge.1) then | |
286 | mdcy(pycomp(111),1)=0 ! no pi0 decay | |
287 | mdcy(pycomp(310),1)=0 ! no K_S0 decay | |
288 | do ihard=1,njet | |
289 | mstp(111)=0 | |
290 | c mstj(41)=0 ! vacuum showering off | |
291 | call pyevnt ! generate hard scattering | |
292 | ||
293 | * generate quenched jets with PYQUEN if nhcel=2 | |
294 | if(nhsel.eq.2.or.nhsel.eq.4) then | |
295 | ifbp=0 | |
296 | Ap=dble(AW) | |
297 | bfixp=dble(RA*bgen) | |
298 | call pyquen(Ap,ifbp,bfixp) | |
299 | end if | |
300 | ||
301 | * fill array of partons | |
302 | nl=nl+n | |
303 | if(nl.gt.150000-np) goto 51 | |
304 | do i=nl-n+1,nl | |
305 | ip=i+n-nl | |
306 | do j=1,5 | |
307 | pl(i,j)=p(ip,j) | |
308 | end do | |
309 | do j=1,5 | |
310 | vl(i,j)=v(ip,j) | |
311 | end do | |
312 | do j=1,5 | |
313 | kl(i,j)=k(ip,j) | |
314 | end do | |
315 | do j=3,5 | |
316 | kk=kl(i,j) | |
317 | if(kk.gt.0) kl(i,j)=kk+nl-n | |
318 | end do | |
319 | end do | |
320 | 51 continue | |
321 | ||
322 | call pyexec ! hadronization done | |
323 | c call pyedit(2) ! remove partons & leave hadrons | |
324 | ||
325 | * fill array of final particles | |
326 | nu=nj+n | |
327 | if(nu.gt.150000-np) then | |
328 | write(6,*) 'Warning, multiplicity too large! Cut hard part.' | |
329 | goto 52 | |
330 | end if | |
331 | nj=nu | |
332 | do i=nj-n+1,nj | |
333 | ip=i+n-nj | |
334 | do j=1,5 | |
335 | pj(i,j)=p(ip,j) | |
336 | end do | |
337 | do j=1,5 | |
338 | vj(i,j)=v(ip,j) | |
339 | end do | |
340 | do j=1,5 | |
341 | kj(i,j)=k(ip,j) | |
342 | end do | |
343 | do j=3,5 | |
344 | kk=kj(i,j) | |
345 | if(kk.gt.0) then | |
346 | kj(i,j)=kk+nj-n | |
347 | end if | |
348 | end do | |
349 | end do | |
350 | ||
351 | end do | |
352 | 52 njet=ihard-1 | |
353 | end if | |
354 | ||
355 | ||
356 | return | |
357 | end | |
358 | ****************************** END HYHARD ************************** | |
359 | ||
360 | ********************************* HIPSEAR *************************** | |
361 | SUBROUTINE HIPSEAR (fmax,xmin) | |
362 | * find maximum and 'sufficient minimum' of differential inelasic AA cross | |
363 | * section as a function of impact paramater (xm, fm are outputs) | |
364 | real hsin | |
365 | external hsin | |
366 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
367 | save /hyipar/ | |
368 | ||
369 | xmin=bmaxh-bminh | |
370 | fmax=0. | |
371 | do j=1,1000 | |
372 | x=bminh+xmin*(j-1)/999. | |
373 | f=hsin(x) | |
374 | if(f.gt.fmax) then | |
375 | fmax=f | |
376 | endif | |
377 | end do | |
378 | ||
379 | return | |
380 | end | |
381 | ****************************** END HIPSEAR ************************** | |
382 | ||
383 | ************************* HARINV ********************************** | |
384 | SUBROUTINE HARINV(X,A,F,N,R) | |
385 | * gives interpolation of function F(X) with arrays A(N) and F(N) | |
386 | DIMENSION A(N),F(N) | |
387 | IF(X.LT.A(1))GO TO 11 | |
388 | IF(X.GT.A(N))GO TO 4 | |
389 | K1=1 | |
390 | K2=N | |
391 | 2 K3=K2-K1 | |
392 | IF(K3.LE.1)GO TO 6 | |
393 | K3=K1+K3/2 | |
394 | IF(A(K3)-X) 7,8,9 | |
395 | 7 K1=K3 | |
396 | GOTO2 | |
397 | 9 K2=K3 | |
398 | GOTO2 | |
399 | 8 P=F(K3) | |
400 | RETURN | |
401 | 3 B1=A(K1) | |
402 | B2=A(K1+1) | |
403 | B3=A(K1+2) | |
404 | B4=F(K1) | |
405 | B5=F(K1+1) | |
406 | B6=F(K1+2) | |
407 | R=B4*((X-B2)*(X-B3))/((B1-B2)*(B1-B3))+B5*((X-B1)*(X-B3))/ | |
408 | 1 ((B2-B1)*(B2-B3))+B6*((X-B1)*(X-B2))/((B3-B1)*(B3-B2)) | |
409 | RETURN | |
410 | 6 IF(K2.NE.N)GO TO 3 | |
411 | K1=N-2 | |
412 | GOTO3 | |
413 | 4 C=ABS(X-A(N)) | |
414 | IF(C.LT.0.1) GO TO 5 | |
415 | K1=N-2 | |
416 | 13 CONTINUE | |
417 | C13 PRINT 41,X | |
418 | C41 FORMAT(25H X IS OUT OF THE INTERVAL,3H X=,F15.9) | |
419 | GO TO 3 | |
420 | 5 R=F(N) | |
421 | RETURN | |
422 | 11 C=ABS(X-A(1)) | |
423 | IF(C.LT.0.1) GO TO 12 | |
424 | K1=1 | |
425 | GOTO 13 | |
426 | 12 R=F(1) | |
427 | RETURN | |
428 | END | |
429 | C************************** END HARINV ************************************* | |
430 | ||
431 | **************************** HSIMPA ********************************** | |
432 | SUBROUTINE HSIMPA (A1,B1,H1,REPS1,AEPS1,FUNCT,X, | |
433 | 1 AI,AIH,AIABS) | |
434 | * calculate intergal of function FUNCT(X) on the interval from A1 to B1 | |
435 | IMPLICIT REAL * 4 ( A-H,O-Z ) | |
436 | IMPLICIT INTEGER(I-N) | |
437 | DIMENSION F(7), P(5) | |
438 | H=SIGN ( H1, B1-A1 ) | |
439 | S=SIGN (1., H ) | |
440 | A=A1 | |
441 | B=B1 | |
442 | AI=0.0 | |
443 | AIH=0.0 | |
444 | AIABS=0.0 | |
445 | P(2)=4. | |
446 | P(4)=4. | |
447 | P(3)=2. | |
448 | P(5)=1. | |
449 | IF(B-A)1,2,1 | |
450 | 1 REPS=ABS(REPS1) | |
451 | AEPS=ABS(AEPS1) | |
452 | DO 3 K=1,7 | |
453 | 3 F(K)=10.E16 | |
454 | X=A | |
455 | C=0. | |
456 | F(1)=FUNCT(X)/3. | |
457 | 4 X0=X | |
458 | IF( (X0+4.*H-B)*S)5,5,6 | |
459 | 6 H=(B-X0)/4. | |
460 | IF ( H ) 7,2,7 | |
461 | 7 DO 8 K=2,7 | |
462 | 8 F(K)=10.E16 | |
463 | C=1. | |
464 | 5 DI2=F (1) | |
465 | DI3=ABS( F(1) ) | |
466 | DO 9 K=2,5 | |
467 | X=X+H | |
468 | IF((X-B)*S)23,24,24 | |
469 | 24 X=B | |
470 | 23 IF(F(K)-10.E16)10,11,10 | |
471 | 11 F(K)=FUNCT(X)/3. | |
472 | 10 DI2=DI2+P(K)*F(K) | |
473 | 9 DI3=DI3+P(K)*ABS(F(K)) | |
474 | DI1=(F(1)+4.*F(3)+F(5))*2.*H | |
475 | DI2=DI2*H | |
476 | DI3=DI3*H | |
477 | IF (REPS) 12,13,12 | |
478 | 13 IF (AEPS) 12,14,12 | |
479 | 12 EPS=ABS((AIABS+DI3)*REPS) | |
480 | IF(EPS-AEPS)15,16,16 | |
481 | 15 EPS=AEPS | |
482 | 16 DELTA=ABS(DI2-DI1) | |
483 | IF(DELTA-EPS)20,21,21 | |
484 | 20 IF(DELTA-EPS/8.)17,14,14 | |
485 | 17 H=2.*H | |
486 | F(1)=F(5) | |
487 | F(2)=F(6) | |
488 | F(3)=F(7) | |
489 | DO 19 K=4,7 | |
490 | 19 F(K)=10.E16 | |
491 | GO TO 18 | |
492 | 14 F(1)=F(5) | |
493 | F(3)=F(6) | |
494 | F(5)=F(7) | |
495 | F(2)=10.E16 | |
496 | F(4)=10.E16 | |
497 | F(6)=10.E16 | |
498 | F(7)=10.E16 | |
499 | 18 DI1=DI2+(DI2-DI1)/15. | |
500 | AI=AI+DI1 | |
501 | AIH=AIH+DI2 | |
502 | AIABS=AIABS+DI3 | |
503 | GO TO 22 | |
504 | 21 H=H/2. | |
505 | F(7)=F(5) | |
506 | F(6)=F(4) | |
507 | F(5)=F(3) | |
508 | F(3)=F(2) | |
509 | F(2)=10.E16 | |
510 | F(4)=10.E16 | |
511 | X=X0 | |
512 | C=0. | |
513 | GO TO 5 | |
514 | 22 IF(C)2,4,2 | |
515 | 2 RETURN | |
516 | END | |
517 | ************************* END HSIMPA ******************************* | |
518 | ||
519 | **************************** HSIMPB ********************************** | |
520 | SUBROUTINE HSIMPB (A1,B1,H1,REPS1,AEPS1,FUNCT,X, | |
521 | 1 AI,AIH,AIABS) | |
522 | * calculate intergal of function FUNCT(X) on the interval from A1 to B1 | |
523 | IMPLICIT REAL * 4 ( A-H,O-Z ) | |
524 | IMPLICIT INTEGER(I-N) | |
525 | DIMENSION F(7), P(5) | |
526 | H=SIGN ( H1, B1-A1 ) | |
527 | S=SIGN (1., H ) | |
528 | A=A1 | |
529 | B=B1 | |
530 | AI=0.0 | |
531 | AIH=0.0 | |
532 | AIABS=0.0 | |
533 | P(2)=4. | |
534 | P(4)=4. | |
535 | P(3)=2. | |
536 | P(5)=1. | |
537 | IF(B-A)1,2,1 | |
538 | 1 REPS=ABS(REPS1) | |
539 | AEPS=ABS(AEPS1) | |
540 | DO 3 K=1,7 | |
541 | 3 F(K)=10.E16 | |
542 | X=A | |
543 | C=0. | |
544 | F(1)=FUNCT(X)/3. | |
545 | 4 X0=X | |
546 | IF( (X0+4.*H-B)*S)5,5,6 | |
547 | 6 H=(B-X0)/4. | |
548 | IF ( H ) 7,2,7 | |
549 | 7 DO 8 K=2,7 | |
550 | 8 F(K)=10.E16 | |
551 | C=1. | |
552 | 5 DI2=F (1) | |
553 | DI3=ABS( F(1) ) | |
554 | DO 9 K=2,5 | |
555 | X=X+H | |
556 | IF((X-B)*S)23,24,24 | |
557 | 24 X=B | |
558 | 23 IF(F(K)-10.E16)10,11,10 | |
559 | 11 F(K)=FUNCT(X)/3. | |
560 | 10 DI2=DI2+P(K)*F(K) | |
561 | 9 DI3=DI3+P(K)*ABS(F(K)) | |
562 | DI1=(F(1)+4.*F(3)+F(5))*2.*H | |
563 | DI2=DI2*H | |
564 | DI3=DI3*H | |
565 | IF (REPS) 12,13,12 | |
566 | 13 IF (AEPS) 12,14,12 | |
567 | 12 EPS=ABS((AIABS+DI3)*REPS) | |
568 | IF(EPS-AEPS)15,16,16 | |
569 | 15 EPS=AEPS | |
570 | 16 DELTA=ABS(DI2-DI1) | |
571 | IF(DELTA-EPS)20,21,21 | |
572 | 20 IF(DELTA-EPS/8.)17,14,14 | |
573 | 17 H=2.*H | |
574 | F(1)=F(5) | |
575 | F(2)=F(6) | |
576 | F(3)=F(7) | |
577 | DO 19 K=4,7 | |
578 | 19 F(K)=10.E16 | |
579 | GO TO 18 | |
580 | 14 F(1)=F(5) | |
581 | F(3)=F(6) | |
582 | F(5)=F(7) | |
583 | F(2)=10.E16 | |
584 | F(4)=10.E16 | |
585 | F(6)=10.E16 | |
586 | F(7)=10.E16 | |
587 | 18 DI1=DI2+(DI2-DI1)/15. | |
588 | AI=AI+DI1 | |
589 | AIH=AIH+DI2 | |
590 | AIABS=AIABS+DI3 | |
591 | GO TO 22 | |
592 | 21 H=H/2. | |
593 | F(7)=F(5) | |
594 | F(6)=F(4) | |
595 | F(5)=F(3) | |
596 | F(3)=F(2) | |
597 | F(2)=10.E16 | |
598 | F(4)=10.E16 | |
599 | X=X0 | |
600 | C=0. | |
601 | GO TO 5 | |
602 | 22 IF(C)2,4,2 | |
603 | 2 RETURN | |
604 | END | |
605 | ************************* END HSIMPB ******************************* | |
606 | ||
607 | * function to calculate differential inelastic AA cross section | |
608 | real function hsin(x) | |
609 | external hftaa | |
610 | real hftaa | |
611 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
612 | save /hyipar/ | |
613 | sigp=sigin | |
614 | taapb0=33.16 | |
615 | br=max(1.e-10,0.25*x*x) | |
616 | hsin=x*(1.-exp(-1.*hftaa(br)*sigp*taapb0)) | |
617 | return | |
618 | end | |
619 | ||
620 | * set of functions to calculate number of nucleons-participants at im.par.b | |
621 | integer function numpar(c) | |
622 | IMPLICIT REAL * 4 (A-H,O-Z) | |
623 | real HFUNC1 | |
624 | external HFUNC1 | |
625 | common /hynup1/ bp,x | |
626 | EPS=0.01 | |
627 | A=0. | |
628 | B=6.28318 | |
629 | H=0.01*(B-A) | |
630 | bp=c | |
631 | CALL HSIMPA(A,B,H,EPS,1.E-8,HFUNC1,X,RES,AIH,AIABS) | |
632 | numpar=int(2.*RES) | |
633 | return | |
634 | end | |
635 | * | |
636 | real function HFUNC1(x) | |
637 | IMPLICIT REAL * 4 (A-H,O-Z) | |
638 | real HFUNC2 | |
639 | external HFUNC2 | |
640 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
641 | common /hynup1/ bp,xx | |
642 | save /hyipar/ | |
643 | xx=x | |
644 | b2=0.5*bp | |
645 | r1=abs(sqrt(abs(RA*RA-(b2*sin(xx))**2))+b2*cos(xx)) | |
646 | r2=abs(sqrt(abs(RA*RA-(b2*sin(xx))**2))-b2*cos(xx)) | |
647 | EPS=0.01 | |
648 | A=0. | |
649 | B=min(r1,r2) | |
650 | H=0.01*(B-A) | |
651 | CALL HSIMPB(A,B,H,EPS,1.E-8,HFUNC2,Y,RES,AIH,AIABS) | |
652 | HFUNC1=RES | |
653 | return | |
654 | end | |
655 | * | |
656 | real function HFUNC2(y) | |
657 | IMPLICIT REAL * 4 (A-H,O-Z) | |
658 | real hythik | |
659 | external hythik | |
660 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
661 | common /hynup1/ bp,x | |
662 | save /hyipar/ | |
663 | r1=sqrt(abs(y*y+bp*bp/4.+y*bp*cos(x))) | |
664 | r2=sqrt(abs(y*y+bp*bp/4.-y*bp*cos(x))) | |
665 | s=1.-exp(-0.1*sigin*hythik(r2)) | |
666 | HFUNC2=y*hythik(r2)*s | |
667 | return | |
668 | end | |
669 | ||
670 | * to calculate nuclear thikness function | |
671 | real function hythik(r) | |
672 | IMPLICIT REAL * 4 (A-H,O-Z) | |
673 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
674 | save /hyipar/ | |
675 | hythik=0.477465*AW*sqrt(abs(RA*RA-r*r))/(RA**3) | |
676 | return | |
677 | end | |
678 | ||
679 | * to calculate nuclear overlap function | |
680 | real function hftaa(br) | |
681 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
682 | save /hyipar/ | |
683 | sbr=sqrt(abs(1.-br)) | |
684 | taa=1.-br*(1.+(1.-0.25*br)*log(1./br)+2.*(1.-0.25*br)* | |
685 | + (log(1.+sbr)-sbr/(1.+sbr))-0.5*br*(1.-br)/((1.+sbr)**2)) | |
686 | hftaa=taa*((AW/207)**1.33333333) | |
687 | return | |
688 | end | |
689 | ||
690 | * function to calculate number of hard parton-parton scatterings with | |
691 | * Q>x at sqrt{s}=5.5 TeV | |
692 | integer function numjet(x) | |
693 | common /hyipar/ bminh,bmaxh,AW,RA,sigin,np | |
694 | common /hyfpar/ bgen,nbcol,npart,npyt,nhyd | |
695 | save /hyipar/,/hyfpar/ | |
696 | dimension ptj(30),sgj(30) | |
697 | data ptj/5.,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.,23.,26., | |
698 | + 30.,35.,40.,50.,60.,70.,80.,90.,100.,120.,150.,200., | |
699 | + 250.,300.,400.,450.,500./ | |
700 | data sgj/31.,16.8,9.5,6.1,4.1,2.7,1.4,0.8,0.46,0.3,0.2, | |
701 | + 0.11,6.8e-2,3.9e-2,2.06e-2,1.16e-2,4.46e-3,2.e-3,1.e-3, | |
702 | + 5.3e-4,3.1e-4,1.9e-4,7.7e-5,2.6e-5,5.9e-6,1.6e-6,5.9e-7, | |
703 | + 1.e-7,4.8e-8,2.3e-8/ | |
704 | pt=x | |
705 | i=0 | |
706 | 31 i=i+1 | |
707 | if(i.le.30) then | |
708 | dx=abs(ptj(i)-pt) | |
709 | if(dx.le.0.1) then | |
710 | sjet=sgj(i) | |
711 | goto 32 | |
712 | else | |
713 | goto 31 | |
714 | end if | |
715 | else | |
716 | call harinv(pt,ptj,sgj,30,sjet) | |
717 | end if | |
718 | 32 pjet=sjet/sigin | |
719 | ijet=0 | |
720 | do i=1,nbcol | |
721 | if(rlu(0).lt.pjet) ijet=ijet+1 | |
722 | end do | |
723 | numjet=ijet | |
724 | return | |
725 | end | |
726 | ||
727 | * function to generate gauss distribution | |
728 | real function hgauss(x0,sig) | |
729 | 41 u1=rlu(0) | |
730 | u2=rlu(0) | |
731 | v1=2.*u1-1. | |
732 | v2=2.*u2-1. | |
733 | s=v1**2+v2**2 | |
734 | if(s.gt.1) go to 41 | |
735 | hgauss=v1*sqrt(-2.*log(s)/s)*sig+x0 | |
736 | return | |
737 | end | |
738 | ||
739 | ||
740 | ************************************************************************** |