]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TAmpt/AMPT/linana.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Filling of empty histogams corrected, Disp2 is not part of Disp any more - corrected
[u/mrichter/AliRoot.git] / TAmpt / AMPT / linana.f
CommitLineData
0119ef9a 1c.................... linana.f
2c=======================================================================
3c 10/26/01 update freezeout positions in case of interactions:
4clin-3/2009 Note: freezeout spacetime values cannot be trusted for K0S & K0L
5c as K0S/K0L are converted from K+/K- by hand at the end of hadron cascade.
6cms
7cms dlw & gsfs Commented out output file writing
8cms
9 subroutine hbtout(nnew,nt,ntmax)
10c
11 PARAMETER (MAXSTR=150001,MAXR=1)
12clin-5/2008 give tolerance to regular particles (perturbative probability 1):
13 PARAMETER (oneminus=0.99999,oneplus=1.00001)
14 dimension lastkp(MAXSTR), newkp(MAXSTR),xnew(3)
15 common /para7/ ioscar,nsmbbbar,nsmmeson
16cc SAVE /para7/
17 COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
18cc SAVE /hbt/
19 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
20cc SAVE /input1/
21 COMMON /AA/ R(3,MAXSTR)
22cc SAVE /AA/
23 COMMON /BB/ P(3,MAXSTR)
24cc SAVE /BB/
25 COMMON /CC/ E(MAXSTR)
26cc SAVE /CC/
27 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
28cc SAVE /EE/
29 common /lastt/itimeh,bimp
30cc SAVE /lastt/
31 COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
32cc SAVE /tdecay/
33 COMMON /AREVT/ IAEVT, IARUN, MISS
34cc SAVE /AREVT/
35 common/snn/efrm,npart1,npart2
36cc SAVE /snn/
37 COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
38cc SAVE /HJGLBR/
39 COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
40 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
41 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
42 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
43clin-12/14/03:
44 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
45 EXTERNAL IARFLV, INVFLV
46 common /para8/ idpert,npertd,idxsec
47 SAVE
48c
49 do 1001 i=1,max0(nlast,nnew)
50 lastkp(i)=0
51 1001 continue
52 do 1002 i=1,nnew
53 newkp(i)=0
54 1002 continue
55c for each of the particles, search the freezeout record (common /hbt/)
56c to find & keep those which do not have interactions during this timestep:
57 do 100 ip=1,nnew
58 do 1004 iplast=1,nlast
59 if(p(1,ip).eq.plast(1,iplast).and.
60 1 p(2,ip).eq.plast(2,iplast).and.
61 2 p(3,ip).eq.plast(3,iplast).and.
62 3 e(ip).eq.plast(4,iplast).and.
63 4 lb(ip).eq.lblast(iplast).and.
64 5 dpertp(ip).eq.dplast(iplast).and.lastkp(iplast).eq.0) then
65clin-5/2008 modified below to the above in case we have perturbative particles:
66c 5 lastkp(iplast).eq.0) then
67 deltat=nt*dt-xlast(4,iplast)
68 ene=sqrt(plast(1,iplast)**2+plast(2,iplast)**2
69 1 +plast(3,iplast)**2+plast(4,iplast)**2)
70c xnew gives the coordinate if a particle free-streams to current time:
71 do 1003 ii=1,3
72 xnew(ii)=xlast(ii,iplast)+plast(ii,iplast)/ene*deltat
73 1003 continue
74 dr=sqrt((r(1,ip)-xnew(1))**2+(r(2,ip)-xnew(2))**2
75 1 +(r(3,ip)-xnew(3))**2)
76c find particles with dp=0 and dr<0.01, considered to be those
77c without any interactions during this timestep,
78c thus keep their last positions and time:
79 if(dr.le.0.01) then
80 lastkp(iplast)=1
81 newkp(ip)=1
82c if(lb(ip).eq.41) then
83c write(95,*) 'nt,ip,px,x=',nt,ip,p(1,ip),r(1,ip),ftsv(ip)
84c write(95,*) 'xnew=',xnew(1),xnew(2),xnew(3),xlast(4,ip)
85c endif
86clin-5/2009 Take care of formation time of particles read in at nt=ntmax-1:
87 if(nt.eq.ntmax.and.ftsv(ip).gt.((ntmax-1)*dt))
88 1 xlast(4,iplast)=ftsv(ip)
89 goto 100
90 endif
91 endif
92 1004 continue
93 100 continue
94c for current particles with interactions, fill their current info in
95c the freezeout record (if that record entry needs not to be kept):
96 do 150 ip=1,nnew
97 if(newkp(ip).eq.0) then
98 do 1005 iplast=1,nnew
99 if(lastkp(iplast).eq.0) then
100ctest off: write collision info
101c if(lb(ip).eq.41) then
102c write(95,*) 'nt,lb(ip)=',nt,lb(ip)
103c write(95,*) ' last p=',plast(1,iplast),
104c 1 plast(2,iplast),plast(3,iplast),plast(4,iplast)
105c write(95,*) ' after p=',p(1,ip),p(2,ip),p(3,ip),e(ip)
106c write(95,*) 'after x=',r(1,ip),r(2,ip),r(3,ip),ftsv(ip)
107c endif
108c
109 xlast(1,iplast)=r(1,ip)
110 xlast(2,iplast)=r(2,ip)
111 xlast(3,iplast)=r(3,ip)
112 xlast(4,iplast)=nt*dt
113c
114 if(nt.eq.ntmax) then
115c freezeout time for decay daughters at the last timestep
116c needs to include the decay time of the parent:
117 if(tfdcy(ip).gt.(ntmax*dt+0.001)) then
118 xlast(4,iplast)=tfdcy(ip)
119c freezeout time for particles unformed at the next-to-last timestep
120c needs to be their formation time instead of (ntmax*dt):
121 elseif(ftsv(ip).gt.((ntmax-1)*dt)) then
122 xlast(4,iplast)=ftsv(ip)
123 endif
124 endif
125 plast(1,iplast)=p(1,ip)
126 plast(2,iplast)=p(2,ip)
127 plast(3,iplast)=p(3,ip)
128 plast(4,iplast)=e(ip)
129 lblast(iplast)=lb(ip)
130 lastkp(iplast)=1
131clin-5/2008:
132 dplast(iplast)=dpertp(ip)
133 goto 150
134 endif
135 1005 continue
136 endif
137 150 continue
138c if the current particle list is shorter than the freezeout record,
139c condense the last-collision record by filling new record from 1 to nnew,
140c and label these entries as keep:
141 if(nnew.lt.nlast) then
142 do 170 iplast=1,nlast
143 if(lastkp(iplast).eq.0) then
144 do 1006 ip2=iplast+1,nlast
145 if(lastkp(ip2).eq.1) then
146 xlast(1,iplast)=xlast(1,ip2)
147 xlast(2,iplast)=xlast(2,ip2)
148 xlast(3,iplast)=xlast(3,ip2)
149 xlast(4,iplast)=xlast(4,ip2)
150 plast(1,iplast)=plast(1,ip2)
151 plast(2,iplast)=plast(2,ip2)
152 plast(3,iplast)=plast(3,ip2)
153 plast(4,iplast)=plast(4,ip2)
154 lblast(iplast)=lblast(ip2)
155 lastkp(iplast)=1
156clin-5/2008:
157 dplast(iplast)=dplast(ip2)
158 goto 170
159 endif
160 1006 continue
161 endif
162 170 continue
163 endif
164 nlast=nnew
165ctest off look inside each NT timestep (for debugging purpose):
166c do ip=1,nlast
167c if(nt.eq.5000) then
168c write(95,*) ' p ',nt,ip,INVFLV(lblast(ip)),plast(1,ip),
169c 1 plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip)
170c write(95,*) ' x ',nt,ip,INVFLV(lblast(ip)),xlast(1,ip),
171c 1 xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip)
172c endif
173c enddo
174c
175 if(nt.eq.ntmax) then
176clin-5/2008 find final number of perturbative particles (deuterons only):
177 ndpert=0
178 do ip=1,nlast
179 if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
180 else
181 ndpert=ndpert+1
182 endif
183 enddo
184c
185c write(16,190) IAEVT,IARUN,nlast,bimp,npart1,npart2,
186c 1 NELP,NINP,NELT,NINTHJ
187c write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2,
188c 1 NELP,NINP,NELT,NINTHJ
189clin-5/2008 write out perturbatively-produced particles (deuterons only):
190c if(idpert.eq.1.or.idpert.eq.2)
191c 1 write(90,190) IAEVT,IARUN,ndpert,bimp,npart1,npart2,
192c 2 NELP,NINP,NELT,NINTHJ
193c do 1007 ip=1,nlast
194clin-12/14/03 No formation time for spectator projectile or target nucleons,
195c see ARINI1 in 'amptsub.f':
196clin-3/2009 To be consistent with new particles produced in hadron cascade
197c that are limited by the time-resolution (DT) of the hadron cascade,
198c freezeout time of spectator projectile or target nucleons is written as
199c DT as they are read at the 1st timestep and then propagated to time DT:
200c if(plast(1,ip).eq.0.and.plast(2,ip).eq.0
201c 1 .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1))
202c 2 .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then
203clin-5/2008 perturbatively-produced particles (currently only deuterons)
204c are written to ana/ampt_pert.dat (without the column for the mass);
205c ana/ampt.dat has regularly-produced particles (including deuterons);
206c these two sets of deuteron data are close to each other(but not the same
207c because of the bias from triggering the perturbative production);
208c ONLY use one data set for analysis to avoid double-counting:
209c if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
210c write(16,200) INVFLV(lblast(ip)), plast(1,ip),
211c 1 plast(2,ip),plast(3,ip),plast(4,ip),
212c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),
213c 3 xlast(4,ip)
214clin-12/14/03-end
215c else
216c if(idpert.eq.1.or.idpert.eq.2) then
217c write(90,250) INVFLV(lblast(ip)), plast(1,ip),
218c 1 plast(2,ip),plast(3,ip),
219c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),
220c 3 xlast(4,ip)
221c else
222c write(99,*) 'Unexpected perturbative particles'
223c endif
224c endif
225c elseif(amax1(abs(xlast(1,ip)),abs(xlast(2,ip)),
226c 1 abs(xlast(3,ip)),abs(xlast(4,ip))).lt.9999) then
227c if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
228c write(16,200) INVFLV(lblast(ip)), plast(1,ip),
229c 1 plast(2,ip),plast(3,ip),plast(4,ip),
230c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
231c else
232c if(idpert.eq.1.or.idpert.eq.2) then
233c write(90,250) INVFLV(lblast(ip)),plast(1,ip),
234c 1 plast(2,ip),plast(3,ip),
235c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
236c 3 dplast(ip)
237c else
238c write(99,*) 'Unexpected perturbative particles'
239c endif
240c endif
241c else
242c change format for large numbers:
243c if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
244c write(16,201) INVFLV(lblast(ip)), plast(1,ip),
245c 1 plast(2,ip),plast(3,ip),plast(4,ip),
246c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
247c else
248c if(idpert.eq.1.or.idpert.eq.2) then
249c write(90,251) INVFLV(lblast(ip)), plast(1,ip),
250c 1 plast(2,ip),plast(3,ip),
251c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
252c 3 dplast(ip)
253c else
254c write(99,*) 'Unexpected perturbative particles'
255c endif
256c endif
257c endif
258c 1007 continue
259 if(ioscar.eq.1) call hoscar
260 endif
261 190 format(3(i7),f10.4,5x,6(i4))
262clin-3/2009 improve the output accuracy of Pz
263 200 format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,f8.2))
264 201 format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,e8.2))
265 250 format(I5,2(1x,f8.3),1x,f10.3,2(1x,f7.1),1x,f8.2,1x,f7.2,1x,e10.4)
266 251 format(I5,2(1x,f8.3),1x,f10.3,4(1x,e8.2),1x,e10.4)
267c
268 return
269 end
270
271c=======================================================================
272 SUBROUTINE decomp(px0,py0,pz0,xm0,i,itq1)
273c
274 IMPLICIT DOUBLE PRECISION(D)
275 DOUBLE PRECISION enenew, pxnew, pynew, pznew
276 DOUBLE PRECISION de0, beta2, gam
3006c44b 277 common /lora/ enenew, pxnew, pynew, pznew
278cc SAVE /lora/
0119ef9a 279 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
280cc SAVE /HPARNT/
281 common /decom/ptwo(2,5)
282cc SAVE /decom/
283 COMMON/RNDF77/NSEED
284cc SAVE /RNDF77/
285 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
286 common/embed/iembed,pxqembd,pyqembd,xembd,yembd
287 SAVE
288c
289 itq1=itq1
290 dcth=dble(RANART(NSEED))*2.d0-1.d0
291 dPHI=dble(RANART(NSEED)*HIPR1(40))*2.d0
292clin-6/2009 Added if embedding a high-Pt quark pair after string melting:
293 if(iembed.eq.1) then
294c Decompose the parent high-Pt pion to q and qbar with an internal momentum
295c parallel to the pion direction so that one parton has ~the same hight Pt
296c and the other parton has a very soft Pt:
297c Note: htop() decomposes a meson to q as it(1) followed by qbar as it(2):
298 if(i.eq.(natt-1).or.i.eq.natt) then
299 dcth=0.d0
300 dphi=dble(acos(pxqembd/sqrt(pxqembd**2+pyqembd**2)))
301 if(pyqembd.lt.0) dphi=HIPR1(40)*2.d0-dphi
302 endif
303 endif
304c
305 ds=dble(xm0)**2
306 dpcm=dsqrt((ds-dble(ptwo(1,5)+ptwo(2,5))**2)
307 1 *(ds-dble(ptwo(1,5)-ptwo(2,5))**2)/ds/4d0)
308 dpz=dpcm*dcth
309 dpx=dpcm*dsqrt(1.d0-dcth**2)*dcos(dphi)
310 dpy=dpcm*dsqrt(1.d0-dcth**2)*dsin(dphi)
311 de1=dsqrt(dble(ptwo(1,5))**2+dpcm**2)
312 de2=dsqrt(dble(ptwo(2,5))**2+dpcm**2)
313c
314 de0=dsqrt(dble(px0)**2+dble(py0)**2+dble(pz0)**2+dble(xm0)**2)
315 dbex=dble(px0)/de0
316 dbey=dble(py0)/de0
317 dbez=dble(pz0)/de0
318c boost the reference frame up by beta (pznew=gam(pz+beta e)):
319 beta2 = dbex ** 2 + dbey ** 2 + dbez ** 2
320 gam = 1.d0 / dsqrt(1.d0 - beta2)
321 if(beta2.ge.0.9999999999999d0) then
322 write(6,*) '1',dbex,dbey,dbez,beta2,gam
323 endif
324c
e2670fbe 325 call lorenza(de1,dpx,dpy,dpz,-dbex,-dbey,-dbez)
0119ef9a 326 ptwo(1,1)=sngl(pxnew)
327 ptwo(1,2)=sngl(pynew)
328 ptwo(1,3)=sngl(pznew)
329 ptwo(1,4)=sngl(enenew)
e2670fbe 330 call lorenza(de2,-dpx,-dpy,-dpz,-dbex,-dbey,-dbez)
0119ef9a 331 ptwo(2,1)=sngl(pxnew)
332 ptwo(2,2)=sngl(pynew)
333 ptwo(2,3)=sngl(pznew)
334 ptwo(2,4)=sngl(enenew)
335c
336 RETURN
337 END
338
339c=======================================================================
340 SUBROUTINE HTOP
341c
342 PARAMETER (MAXSTR=150001)
343 PARAMETER (MAXPTN=400001)
344 PARAMETER (MAXIDL=4001)
345 DOUBLE PRECISION GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
346 DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
347 1 GXSGS,GYSGS,GZSGS,FTSGS
348 dimension it(4)
349 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
350cwei DOUBLE PRECISION PATT
351cc SAVE /HMAIN2/
352 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
353cc SAVE /HMAIN1/
354 COMMON /PARA1/ MUL
355cc SAVE /PARA1/
356 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
357 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
358 & XMASS0(MAXPTN), ITYP0(MAXPTN)
359cc SAVE /prec1/
360 COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
361cc SAVE /ilist7/
362 COMMON /ARPRC/ ITYPAR(MAXSTR),
363 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
364 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
365 & XMAR(MAXSTR)
366cc SAVE /ARPRC/
367 common /decom/ptwo(2,5)
368cc SAVE /decom/
369 COMMON/RNDF77/NSEED
370cc SAVE /RNDF77/
371 COMMON /NOPREC/ NNOZPC, ITYPN(MAXIDL),
372 & GXN(MAXIDL), GYN(MAXIDL), GZN(MAXIDL), FTN(MAXIDL),
373 & PXN(MAXIDL), PYN(MAXIDL), PZN(MAXIDL), EEN(MAXIDL),
374 & XMN(MAXIDL)
375cc SAVE /NOPREC/
376 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
377cc SAVE /HPARNT/
378c 7/20/01: use double precision
e2670fbe 379c otherwise sometimes beta>1 and gamma diverge in lorenza():
0119ef9a 380 COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
381 & PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
382 & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
383 & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
384cc SAVE /SOFT/
385 common/anim/nevent,isoft,isflag,izpc
386cc SAVE /anim/
387 DOUBLE PRECISION vxp0,vyp0,vzp0
388 common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
389cc SAVE /precpa/
390 common /para7/ ioscar,nsmbbbar,nsmmeson
391 COMMON /AREVT/ IAEVT, IARUN, MISS
392 SAVE
393c
394 npar=0
395 nnozpc=0
396clin-5b/2008 calculate the number of hadrons to be converted to q/qbar:
397 if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
398 1 then
399 nsmbbbar=0
400 nsmmeson=0
401 do i=1,natt
402 id=ITYPAR(i)
403 idabs=iabs(id)
404 i2=MOD(idabs/10,10)
405 if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
406 1 .ge.(HINT1(1)/2*0.99).and.
407 2 ((id.eq.2112).or.(id.eq.2212))) then
408c proj or targ nucleons without interactions, do not enter ZPC:
409 elseif(idabs.gt.1000.and.i2.ne.0) then
410c baryons to be converted to q/qbar:
411 nsmbbbar=nsmbbbar+1
412 elseif((idabs.gt.100.and.idabs.lt.1000)
413 1 .or.idabs.gt.10000) then
414c mesons to be converted to q/qbar:
415 nsmmeson=nsmmeson+1
416 endif
417 enddo
418
419clin-6/2009:
420 if(ioscar.eq.2.or.ioscar.eq.3) then
421 write(92,*) iaevt,miss,3*nsmbbbar+2*nsmmeson,
422 1 nsmbbbar,nsmmeson,natt,natt-nsmbbbar-nsmmeson
423 endif
424c write(92,*) iaevt, 3*nsmbbbar+2*nsmmeson
425c write(92,*) ' event#, total # of initial partons after string
426c 1 melting'
427c write(92,*) 'String melting converts ',nsmbbbar, ' baryons &'
428c 1, nsmmeson, ' mesons'
429c write(92,*) 'Total # of initial particles= ',natt
430c write(92,*) 'Total # of initial particles (gamma,e,muon,...)
431c 1 not entering ZPC= ',natt-nsmbbbar-nsmmeson
432 endif
433clin-5b/2008-over
434 do 100 i=1,natt
435 id=ITYPAR(i)
436 idabs=iabs(id)
437 i4=MOD(idabs/1000,10)
438 i3=MOD(idabs/100,10)
439 i2=MOD(idabs/10,10)
440 i1=MOD(idabs,10)
441 rnum=RANART(NSEED)
442 ftime=0.197*PEAR(i)/(PXAR(i)**2+PYAR(i)**2+XMAR(i)**2)
443 inozpc=0
444 it(1)=0
445 it(2)=0
446 it(3)=0
447 it(4)=0
448c
449 if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
450 1 .ge.(HINT1(1)/2*0.99).and.((id.eq.2112).or.(id.eq.2212))) then
451c proj or targ nucleons without interactions, do not enter ZPC:
452 inozpc=1
453 elseif(idabs.gt.1000.and.i2.ne.0) then
454c baryons:
455 if(((i4.eq.1.or.i4.eq.2).and.i4.eq.i3)
456 1 .or.(i4.eq.3.and.i3.eq.3)) then
457 if(i1.eq.2) then
458 if(rnum.le.(1./2.)) then
459 it(1)=i4
460 it(2)=i3*1000+i2*100+1
461 elseif(rnum.le.(2./3.)) then
462 it(1)=i4
463 it(2)=i3*1000+i2*100+3
464 else
465 it(1)=i2
466 it(2)=i4*1000+i3*100+3
467 endif
468 elseif(i1.eq.4) then
469 if(rnum.le.(2./3.)) then
470 it(1)=i4
471 it(2)=i3*1000+i2*100+3
472 else
473 it(1)=i2
474 it(2)=i4*1000+i3*100+3
475 endif
476 endif
477 elseif(i4.eq.1.or.i4.eq.2) then
478 if(i1.eq.2) then
479 if(rnum.le.(1./2.)) then
480 it(1)=i2
481 it(2)=i4*1000+i3*100+1
482 elseif(rnum.le.(2./3.)) then
483 it(1)=i2
484 it(2)=i4*1000+i3*100+3
485 else
486 it(1)=i4
487 it(2)=i3*1000+i2*100+3
488 endif
489 elseif(i1.eq.4) then
490 if(rnum.le.(2./3.)) then
491 it(1)=i2
492 it(2)=i4*1000+i3*100+3
493 else
494 it(1)=i4
495 it(2)=i3*1000+i2*100+3
496 endif
497 endif
498 elseif(i4.ge.3) then
499 it(1)=i4
500 if(i3.lt.i2) then
501 it(2)=i2*1000+i3*100+1
502 else
503 it(2)=i3*1000+i2*100+3
504 endif
505 endif
506c antibaryons:
507 if(id.lt.0) then
508 it(1)=-it(1)
509 it(2)=-it(2)
510 endif
511c isoft=4or5 decompose diquark flavor it(2) to two quarks it(3)&(4):
512 if(isoft.eq.4.or.isoft.eq.5) then
513 it(3)=MOD(it(2)/1000,10)
514 it(4)=MOD(it(2)/100,10)
515 endif
516
517 elseif((idabs.gt.100.and.idabs.lt.1000)
518 1 .or.idabs.gt.10000) then
519c mesons:
520 if(i3.eq.i2) then
521 if(i3.eq.1.or.i3.eq.2) then
522 if(rnum.le.0.5) then
523 it(1)=1
524 it(2)=-1
525 else
526 it(1)=2
527 it(2)=-2
528 endif
529 else
530 it(1)=i3
531 it(2)=-i3
532 endif
533 else
534 if((isign(1,id)*(-1)**i3).eq.1) then
535 it(1)=i3
536 it(2)=-i2
537 else
538 it(1)=i2
539 it(2)=-i3
540 endif
541 endif
542 else
543c save other particles (leptons and photons) outside of ZPC:
544 inozpc=1
545 endif
546c
547 if(inozpc.eq.1) then
548 NJSGS(i)=0
549 nnozpc=nnozpc+1
550 itypn(nnozpc)=ITYPAR(i)
551 pxn(nnozpc)=PXAR(i)
552 pyn(nnozpc)=PYAR(i)
553 pzn(nnozpc)=PZAR(i)
554 een(nnozpc)=PEAR(i)
555 xmn(nnozpc)=XMAR(i)
556 gxn(nnozpc)=GXAR(i)
557 gyn(nnozpc)=GYAR(i)
558 gzn(nnozpc)=GZAR(i)
559 ftn(nnozpc)=FTAR(i)
560 else
561 NJSGS(i)=2
562 ptwo(1,5)=ulmass(it(1))
563 ptwo(2,5)=ulmass(it(2))
564 call decomp(patt(i,1),patt(i,2),patt(i,3),XMAR(i),i,it(1))
565 ipamax=2
566 if((isoft.eq.4.or.isoft.eq.5)
567 1 .and.iabs(it(2)).gt.1000) ipamax=1
568 do 1001 ipar=1,ipamax
569 npar=npar+1
570 ityp0(npar)=it(ipar)
571 px0(npar)=dble(ptwo(ipar,1))
572 py0(npar)=dble(ptwo(ipar,2))
573 pz0(npar)=dble(ptwo(ipar,3))
574 e0(npar)=dble(ptwo(ipar,4))
575 xmass0(npar)=dble(ptwo(ipar,5))
576 gx0(npar)=dble(GXAR(i))
577 gy0(npar)=dble(GYAR(i))
578 gz0(npar)=dble(GZAR(i))
579 ft0(npar)=dble(ftime)
580 lstrg0(npar)=i
581 lpart0(npar)=ipar
582 vxp0(npar)=dble(patt(i,1)/patt(i,4))
583 vyp0(npar)=dble(patt(i,2)/patt(i,4))
584 vzp0(npar)=dble(patt(i,3)/patt(i,4))
585 1001 continue
586cyy 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
587cyy 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
588c
589 if((isoft.eq.4.or.isoft.eq.5)
590 1 .and.iabs(it(2)).gt.1000) then
591 NJSGS(i)=3
592 xmdq=ptwo(2,5)
593 ptwo(1,5)=ulmass(it(3))
594 ptwo(2,5)=ulmass(it(4))
595c 8/19/02 avoid actual argument in common blocks of DECOMP:
596c call decomp(ptwo(2,1),ptwo(2,2),ptwo(2,3),xmdq)
597 ptwox=ptwo(2,1)
598 ptwoy=ptwo(2,2)
599 ptwoz=ptwo(2,3)
600 call decomp(ptwox,ptwoy,ptwoz,xmdq,i,it(1))
601c
602 do 1002 ipar=1,2
603 npar=npar+1
604 ityp0(npar)=it(ipar+2)
605 px0(npar)=dble(ptwo(ipar,1))
606 py0(npar)=dble(ptwo(ipar,2))
607 pz0(npar)=dble(ptwo(ipar,3))
608 e0(npar)=dble(ptwo(ipar,4))
609 xmass0(npar)=dble(ptwo(ipar,5))
610 gx0(npar)=dble(GXAR(i))
611 gy0(npar)=dble(GYAR(i))
612 gz0(npar)=dble(GZAR(i))
613 ft0(npar)=dble(ftime)
614 lstrg0(npar)=i
615 lpart0(npar)=ipar+1
616 vxp0(npar)=dble(patt(i,1)/patt(i,4))
617 vyp0(npar)=dble(patt(i,2)/patt(i,4))
618 vzp0(npar)=dble(patt(i,3)/patt(i,4))
619 1002 continue
620 endif
621c
622 endif
623 100 continue
624 MUL=NPAR
625c
626clin-5b/2008:
627 if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
628 1 then
629 if((natt-nsmbbbar-nsmmeson).ne.nnozpc)
630 1 write(92,*) 'Problem with the total # of initial particles
631 2 (gamma,e,muon,...) not entering ZPC'
632 if((3*nsmbbbar+2*nsmmeson).ne.npar)
633 1 write(92,*) 'Problem with the total # of initial partons
634 2 after string melting'
635 endif
636c
637 RETURN
638 END
639
640c=======================================================================
641 SUBROUTINE PTOH
642c
643 PARAMETER (MAXSTR=150001)
644 DOUBLE PRECISION gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
645 DOUBLE PRECISION gxp0,gyp0,gzp0,ft0fom,drlocl
646 DOUBLE PRECISION enenew, pxnew, pynew, pznew, beta2, gam
647 DOUBLE PRECISION ftavg0,gxavg0,gyavg0,gzavg0,bex,bey,bez
648 DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
649 1 GXSGS,GYSGS,GZSGS,FTSGS
650 DOUBLE PRECISION xmdiag,px1,py1,pz1,e1,px2,py2,pz2,e2,
651 1 px3,py3,pz3,e3,xmpair,etot
652 common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
653 1 pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
654cc SAVE /loclco/
655 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
656cc SAVE /HMAIN1/
657 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
658cc SAVE /HMAIN2/
659 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
660 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
661 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
662cc SAVE /HJJET2/
663 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
664cc SAVE /ARPRNT/
665 COMMON /ARPRC/ ITYPAR(MAXSTR),
666 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
667 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
668 & XMAR(MAXSTR)
669cc SAVE /ARPRC/
670 COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
671 & PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
672 & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
673 & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
674cc SAVE /SOFT/
675 COMMON/RNDF77/NSEED
676cc SAVE /RNDF77/
677 common/anim/nevent,isoft,isflag,izpc
678cc SAVE /anim/
679 common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
680cc SAVE /prtn23/
681 common /nzpc/nattzp
682cc SAVE /nzpc/
3006c44b 683 common /lora/ enenew, pxnew, pynew, pznew
684cc SAVE /lora/
0119ef9a 685 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
686cc SAVE /LUDAT1A/
687clin 4/19/2006
688 common /lastt/itimeh,bimp
689 COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
690 COMMON /AREVT/ IAEVT, IARUN, MISS
691 common /para7/ ioscar,nsmbbbar,nsmmeson
692c
693 dimension xmdiag(MAXSTR),indx(MAXSTR),ndiag(MAXSTR)
694cwei DOUBLE PRECISION PATT
695 SAVE
696c
697 call coales
698c obtain particle mass here without broadening by Breit-Wigner width:
699 mstj24=MSTJ(24)
700 MSTJ(24)=0
701 nuudd=0
702 npich=0
703 nrhoch=0
704 ppi0=1.
705 prho0=0.
706c determine hadron flavor except (pi0,rho0,eta,omega):
707 DO 1001 ISG = 1, NSG
708 if(NJSGS(ISG).ne.0) then
709 NATT=NATT+1
710 K1=K2SGS(ISG,1)
711 k1abs=iabs(k1)
712 PX1=PXSGS(ISG,1)
713 PY1=PYSGS(ISG,1)
714 PZ1=PZSGS(ISG,1)
715 K2=K2SGS(ISG,2)
716 k2abs=iabs(k2)
717 PX2=PXSGS(ISG,2)
718 PY2=PYSGS(ISG,2)
719 PZ2=PZSGS(ISG,2)
720c 5/02/01 try lowest spin states as first choices,
721c i.e. octet baryons and pseudoscalar mesons (ibs=2*baryonspin+1):
722 e1=PESGS(ISG,1)
723 e2=PESGS(ISG,2)
724 xmpair=dsqrt((e1+e2)**2-(px1+px2)**2-(py1+py2)**2
725 1 -(pz1+pz2)**2)
726 ibs=2
727 imspin=0
728 if(k1.eq.-k2.and.iabs(k1).le.2.
729 1 and.NJSGS(ISG).eq.2) then
730 nuudd=nuudd+1
731 xmdiag(nuudd)=xmpair
732 ndiag(nuudd)=natt
733 endif
734 K3=0
735 if((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3) then
736 K3=K2SGS(ISG,3)
737 k3abs=iabs(k3)
738 PX3=PXSGS(ISG,3)
739 PY3=PYSGS(ISG,3)
740 PZ3=PZSGS(ISG,3)
741 e3=PESGS(ISG,3)
742 xmpair=dsqrt((e1+e2+e3)**2-(px1+px2+px3)**2
743 1 -(py1+py2+py3)**2-(pz1+pz2+pz3)**2)
744 endif
745c***** isoft=3 baryon decomposition is different:
746 if(isoft.eq.3.and.
747 1 (k1abs.gt.1000.or.k2abs.gt.1000)) then
748 if(k1abs.gt.1000) then
749 kdq=k1abs
750 kk=k2abs
751 else
752 kdq=k2abs
753 kk=k1abs
754 endif
755 ki=MOD(kdq/1000,10)
756 kj=MOD(kdq/100,10)
757 if(MOD(kdq,10).eq.1) then
758 idqspn=0
759 else
760 idqspn=1
761 endif
762c
763 if(kk.gt.ki) then
764 ktemp=kk
765 kk=kj
766 kj=ki
767 ki=ktemp
768 elseif(kk.gt.kj) then
769 ktemp=kk
770 kk=kj
771 kj=ktemp
772 endif
773c
774 if(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
775 if(idqspn.eq.0) then
776 kf=1000*ki+100*kk+10*kj+ibs
777 else
778 kf=1000*ki+100*kj+10*kk+ibs
779 endif
780 elseif(ki.eq.kj.and.ki.eq.kk) then
781c can only be decuplet baryons:
782 kf=1000*ki+100*kj+10*kk+4
783 else
784 kf=1000*ki+100*kj+10*kk+ibs
785 endif
786c form a decuplet baryon if the q+diquark mass is closer to its mass
787c (and if the diquark has spin 1):
788cc for now only include Delta, which is present in ART:
789cc if(idqspn.eq.1.and.MOD(kf,10).eq.2) then
790 if(kf.eq.2112.or.kf.eq.2212) then
791 if(abs(sngl(xmpair)-ULMASS(kf)).gt.
792 1 abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
793 endif
794 if(k1.lt.0) kf=-kf
795clin-6/22/01 isoft=4or5 baryons:
796 elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
797 1 then
798 if(k1abs.gt.k2abs) then
799 ki=k1abs
800 kk=k2abs
801 else
802 ki=k2abs
803 kk=k1abs
804 endif
805 if(k3abs.gt.ki) then
806 kj=ki
807 ki=k3abs
808 elseif(k3abs.lt.kk) then
809 kj=kk
810 kk=k3abs
811 else
812 kj=k3abs
813 endif
814c
815 if(ki.eq.kj.and.ki.eq.kk) then
816c can only be decuplet baryons (Delta-,++, Omega):
817 ibs=4
818 kf=1000*ki+100*kj+10*kk+ibs
819 elseif(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
820c form Lambda or Sigma according to 3-quark mass,
821c for now neglect decuplet (Sigma*0 etc) which is absent in ART:
822 ibs=2
823 kf1=1000*ki+100*kj+10*kk+ibs
824 kf2=1000*ki+100*kk+10*kj+ibs
825 kf=kf1
826 if(abs(sngl(xmpair)-ULMASS(kf1)).gt.
827 1 abs(sngl(xmpair)-ULMASS(kf2))) kf=kf2
828 else
829 ibs=2
830 kf=1000*ki+100*kj+10*kk+ibs
831cc for now only include Delta0,+ as decuplets, which are present in ART:
832 if(kf.eq.2112.or.kf.eq.2212) then
833 if(abs(sngl(xmpair)-ULMASS(kf)).gt.
834 1 abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
835 endif
836 endif
837 if(k1.lt.0) kf=-kf
838c***** mesons:
839 else
840 if(k1abs.eq.k2abs) then
841 if(k1abs.le.2) then
842c treat diagonal mesons later in the subroutine:
843 kf=0
844 elseif(k1abs.le.3) then
845c do not form eta', only form phi from s-sbar, since no eta' in ART:
846 kf=333
847 else
848 kf=100*k1abs+10*k1abs+2*imspin+1
849 endif
850 else
851 if(k1abs.gt.k2abs) then
852 kmax=k1abs
853 kmin=k2abs
854 elseif(k1abs.lt.k2abs) then
855 kmax=k2abs
856 kmin=k1abs
857 endif
858 kf=(100*kmax+10*kmin+2*imspin+1)
859 1 *isign(1,k1+k2)*(-1)**kmax
860c form a vector meson if the q+qbar mass is closer to its mass:
861 if(MOD(iabs(kf),10).eq.1) then
862 if(abs(sngl(xmpair)-ULMASS(iabs(kf))).gt.
863 1 abs(sngl(xmpair)-ULMASS(iabs(kf)+2)))
864 2 kf=(iabs(kf)+2)*isign(1,kf)
865 endif
866 endif
867 endif
868 ITYPAR(NATT)=kf
869 KATT(NATT,1)=kf
870 if(iabs(kf).eq.211) then
871 npich=npich+1
872 elseif(iabs(kf).eq.213) then
873 nrhoch=nrhoch+1
874 endif
875 endif
876 1001 CONTINUE
877c assume Npi0=(Npi+ + Npi-)/2, Nrho0=(Nrho+ + Nrho-)/2 on the average:
878 if(nuudd.ne.0) then
879 ppi0=float(npich/2)/float(nuudd)
880 prho0=float(nrhoch/2)/float(nuudd)
881 endif
882c determine diagonal mesons (pi0,rho0,eta and omega) from uubar/ddbar:
883 npi0=0
884 DO 1002 ISG = 1, NSG
885 if(K2SGS(ISG,1).eq.-K2SGS(ISG,2)
886 1 .and.iabs(K2SGS(ISG,1)).le.2.and.NJSGS(ISG).eq.2) then
887 if(RANART(NSEED).le.ppi0) npi0=npi0+1
888 endif
889 1002 CONTINUE
890c
891 if(nuudd.gt.1) then
892 call index1(MAXSTR,nuudd,xmdiag,indx)
893 else
894 indx(1)=1
895 end if
896c
897 DO 1003 ix=1,nuudd
898 iuudd=indx(ix)
899 inatt=ndiag(iuudd)
900 if(ix.le.npi0) then
901 kf=111
902 elseif(RANART(NSEED).le.(prho0/(1-ppi0+0.00001))) then
903 kf=113
904 else
905c at T=150MeV, thermal weights for eta and omega(spin1) are about the same:
906 if(RANART(NSEED).le.0.5) then
907 kf=221
908 else
909 kf=223
910 endif
911 endif
912 ITYPAR(inatt)=kf
913 KATT(inatt,1)=kf
914 1003 CONTINUE
915c determine hadron formation time, position and momentum:
916 inatt=0
917clin-6/2009 write out parton info after coalescence:
918 if(ioscar.eq.3) then
919 WRITE (85, 395) IAEVT, 3*nsmbbbar+2*nsmmeson,nsmbbbar,nsmmeson,
920 1 bimp, NELP,NINP,NELT,NINTHJ,MISS
921 endif
922 395 format(4I8,f10.4,5I5)
923c
924 DO 1006 ISG = 1, NSG
925 if(NJSGS(ISG).ne.0) then
926 inatt=inatt+1
927 K1=K2SGS(ISG,1)
928 k1abs=iabs(k1)
929 PX1=PXSGS(ISG,1)
930 PY1=PYSGS(ISG,1)
931 PZ1=PZSGS(ISG,1)
932 K2=K2SGS(ISG,2)
933 k2abs=iabs(k2)
934 PX2=PXSGS(ISG,2)
935 PY2=PYSGS(ISG,2)
936 PZ2=PZSGS(ISG,2)
937 e1=PESGS(ISG,1)
938 e2=PESGS(ISG,2)
939c
940 if(NJSGS(ISG).eq.2) then
941 PXAR(inatt)=sngl(px1+px2)
942 PYAR(inatt)=sngl(py1+py2)
943 PZAR(inatt)=sngl(pz1+pz2)
944 PATT(inatt,1)=PXAR(inatt)
945 PATT(inatt,2)=PYAR(inatt)
946 PATT(inatt,3)=PZAR(inatt)
947 etot=e1+e2
948 elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
949 1 then
950 PX3=PXSGS(ISG,3)
951 PY3=PYSGS(ISG,3)
952 PZ3=PZSGS(ISG,3)
953 e3=PESGS(ISG,3)
954 PXAR(inatt)=sngl(px1+px2+px3)
955 PYAR(inatt)=sngl(py1+py2+py3)
956 PZAR(inatt)=sngl(pz1+pz2+pz3)
957 PATT(inatt,1)=PXAR(inatt)
958 PATT(inatt,2)=PYAR(inatt)
959 PATT(inatt,3)=PZAR(inatt)
960 etot=e1+e2+e3
961 endif
962 XMAR(inatt)=ULMASS(ITYPAR(inatt))
963 PEAR(inatt)=sqrt(PXAR(inatt)**2+PYAR(inatt)**2
964 1 +PZAR(inatt)**2+XMAR(inatt)**2)
965 PATT(inatt,4)=PEAR(inatt)
966 EATT=EATT+PEAR(inatt)
967 ipartn=NJSGS(ISG)
968 DO 1004 i=1,ipartn
969 ftp(i)=ftsgs(isg,i)
970 gxp(i)=gxsgs(isg,i)
971 gyp(i)=gysgs(isg,i)
972 gzp(i)=gzsgs(isg,i)
973 pxp(i)=pxsgs(isg,i)
974 pyp(i)=pysgs(isg,i)
975 pzp(i)=pzsgs(isg,i)
976 pmp(i)=pmsgs(isg,i)
977 pep(i)=pesgs(isg,i)
978 1004 CONTINUE
979 call locldr(ipartn,drlocl)
980c
981 tau0=ARPAR1(1)
982 ftavg0=ft0fom+dble(tau0)
983 gxavg0=0d0
984 gyavg0=0d0
985 gzavg0=0d0
986 DO 1005 i=1,ipartn
987 gxavg0=gxavg0+gxp0(i)/ipartn
988 gyavg0=gyavg0+gyp0(i)/ipartn
989 gzavg0=gzavg0+gzp0(i)/ipartn
990 1005 CONTINUE
991 bex=dble(PXAR(inatt))/etot
992 bey=dble(PYAR(inatt))/etot
993 bez=dble(PZAR(inatt))/etot
994 beta2 = bex ** 2 + bey ** 2 + bez ** 2
995 gam = 1.d0 / dsqrt(1.d0 - beta2)
996 if(beta2.ge.0.9999999999999d0) then
997 write(6,*) '2',bex,bey,bez,beta2,gam
998 endif
999c
e2670fbe 1000 call lorenza(ftavg0,gxavg0,gyavg0,gzavg0,-bex,-bey,-bez)
0119ef9a 1001 GXAR(inatt)=sngl(pxnew)
1002 GYAR(inatt)=sngl(pynew)
1003 GZAR(inatt)=sngl(pznew)
1004 FTAR(inatt)=sngl(enenew)
1005clin 4/19/2006 write out parton info after coalescence:
1006 if(ioscar.eq.3) then
1007 WRITE (85, 312) K2SGS(ISG,1),px1,py1,pz1,PMSGS(ISG,1),
1008 1 inatt,katt(inatt,1)
1009 WRITE (85, 312) K2SGS(ISG,2),px2,py2,pz2,PMSGS(ISG,2),
1010 1 inatt,katt(inatt,1)
1011 if(NJSGS(ISG).eq.3) WRITE (85, 312) K2SGS(ISG,3),
1012 1 px3,py3,pz3,PMSGS(ISG,3),inatt,katt(inatt,1)
1013 endif
1014 312 FORMAT(I6,4(1X,F10.3),1X,I6,1X,I6)
1015c
1016 endif
1017 1006 CONTINUE
1018c number of hadrons formed from partons inside ZPC:
1019 nattzp=natt
1020 MSTJ(24)=mstj24
1021c
1022 RETURN
1023 END
1024
1025c=======================================================================
1026 SUBROUTINE coales
1027
1028 PARAMETER (MAXSTR=150001)
1029 IMPLICIT DOUBLE PRECISION(D)
1030 DOUBLE PRECISION gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
1031 DIMENSION IOVER(MAXSTR),dp1(2:3),dr1(2:3)
1032 DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
1033 1 GXSGS,GYSGS,GZSGS,FTSGS
1034 double precision dpcoal,drcoal,ecritl
1035 COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
1036 & PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
1037 & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
1038 & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
1039cc SAVE /SOFT/
1040 common /coal/dpcoal,drcoal,ecritl
1041cc SAVE /coal/
1042 common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1043 1 pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
1044cc SAVE /loclco/
1045 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
1046 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
1047 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
1048cc SAVE /HJJET2/
1049 SAVE
1050c
1051 do 1001 ISG=1, NSG
1052 IOVER(ISG)=0
1053 1001 continue
1054C1 meson q coalesce with all available qbar:
1055 do 150 ISG=1,NSG
1056 if(NJSGS(ISG).ne.2.or.IOVER(ISG).eq.1) goto 150
1057C DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1058 if(K2SGS(ISG,1).lt.0) then
1059 write(6,*) 'Antiquark appears in quark loop; stop'
1060 stop
1061 endif
1062c
1063 do 1002 j=1,2
1064 ftp(j)=ftsgs(isg,j)
1065 gxp(j)=gxsgs(isg,j)
1066 gyp(j)=gysgs(isg,j)
1067 gzp(j)=gzsgs(isg,j)
1068 pxp(j)=pxsgs(isg,j)
1069 pyp(j)=pysgs(isg,j)
1070 pzp(j)=pzsgs(isg,j)
1071 pmp(j)=pmsgs(isg,j)
1072 pep(j)=pesgs(isg,j)
1073 1002 continue
1074 call locldr(2,drlocl)
1075 dr0=drlocl
1076c dp0^2 defined as (p1+p2)^2-(m1+m2)^2:
1077 dp0=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1078 & -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1079c
1080 do 120 JSG=1,NSG
1081c skip default or unavailable antiquarks:
1082 if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 120
1083 if(NJSGS(JSG).eq.2) then
1084 ipmin=2
1085 ipmax=2
1086 elseif(NJSGS(JSG).eq.3.and.K2SGS(JSG,1).lt.0) then
1087 ipmin=1
1088 ipmax=3
1089 else
1090 goto 120
1091 endif
1092 do 100 ip=ipmin,ipmax
1093 dplocl=dsqrt(2*(pep(1)*pesgs(jsg,ip)
1094 1 -pxp(1)*pxsgs(jsg,ip)
1095 2 -pyp(1)*pysgs(jsg,ip)
1096 3 -pzp(1)*pzsgs(jsg,ip)
1097 4 -pmp(1)*pmsgs(jsg,ip)))
1098c skip if outside of momentum radius:
1099 if(dplocl.gt.dpcoal) goto 120
1100 ftp(2)=ftsgs(jsg,ip)
1101 gxp(2)=gxsgs(jsg,ip)
1102 gyp(2)=gysgs(jsg,ip)
1103 gzp(2)=gzsgs(jsg,ip)
1104 pxp(2)=pxsgs(jsg,ip)
1105 pyp(2)=pysgs(jsg,ip)
1106 pzp(2)=pzsgs(jsg,ip)
1107 pmp(2)=pmsgs(jsg,ip)
1108 pep(2)=pesgs(jsg,ip)
1109 call locldr(2,drlocl)
1110c skip if outside of spatial radius:
1111 if(drlocl.gt.drcoal) goto 120
1112c q_isg coalesces with qbar_jsg:
1113 if((dp0.gt.dpcoal.or.dr0.gt.drcoal)
1114 1 .or.(drlocl.lt.dr0)) then
1115 dp0=dplocl
1116 dr0=drlocl
1117 call exchge(isg,2,jsg,ip)
1118 endif
1119 100 continue
1120 120 continue
1121 if(dp0.le.dpcoal.and.dr0.le.drcoal) IOVER(ISG)=1
1122 150 continue
1123c
1124C2 meson qbar coalesce with all available q:
1125 do 250 ISG=1,NSG
1126 if(NJSGS(ISG).ne.2.or.IOVER(ISG).eq.1) goto 250
1127C DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1128 do 1003 j=1,2
1129 ftp(j)=ftsgs(isg,j)
1130 gxp(j)=gxsgs(isg,j)
1131 gyp(j)=gysgs(isg,j)
1132 gzp(j)=gzsgs(isg,j)
1133 pxp(j)=pxsgs(isg,j)
1134 pyp(j)=pysgs(isg,j)
1135 pzp(j)=pzsgs(isg,j)
1136 pmp(j)=pmsgs(isg,j)
1137 pep(j)=pesgs(isg,j)
1138 1003 continue
1139 call locldr(2,drlocl)
1140 dr0=drlocl
1141 dp0=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1142 & -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1143c
1144 do 220 JSG=1,NSG
1145 if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 220
1146 if(NJSGS(JSG).eq.2) then
1147 ipmin=1
1148 ipmax=1
1149 elseif(NJSGS(JSG).eq.3.and.K2SGS(JSG,1).gt.0) then
1150 ipmin=1
1151 ipmax=3
1152 else
1153 goto 220
1154 endif
1155 do 200 ip=ipmin,ipmax
1156 dplocl=dsqrt(2*(pep(2)*pesgs(jsg,ip)
1157 1 -pxp(2)*pxsgs(jsg,ip)
1158 2 -pyp(2)*pysgs(jsg,ip)
1159 3 -pzp(2)*pzsgs(jsg,ip)
1160 4 -pmp(2)*pmsgs(jsg,ip)))
1161c skip if outside of momentum radius:
1162 if(dplocl.gt.dpcoal) goto 220
1163 ftp(1)=ftsgs(jsg,ip)
1164 gxp(1)=gxsgs(jsg,ip)
1165 gyp(1)=gysgs(jsg,ip)
1166 gzp(1)=gzsgs(jsg,ip)
1167 pxp(1)=pxsgs(jsg,ip)
1168 pyp(1)=pysgs(jsg,ip)
1169 pzp(1)=pzsgs(jsg,ip)
1170 pmp(1)=pmsgs(jsg,ip)
1171 pep(1)=pesgs(jsg,ip)
1172 call locldr(2,drlocl)
1173c skip if outside of spatial radius:
1174 if(drlocl.gt.drcoal) goto 220
1175c qbar_isg coalesces with q_jsg:
1176 if((dp0.gt.dpcoal.or.dr0.gt.drcoal)
1177 1 .or.(drlocl.lt.dr0)) then
1178 dp0=dplocl
1179 dr0=drlocl
1180 call exchge(isg,1,jsg,ip)
1181 endif
1182 200 continue
1183 220 continue
1184 if(dp0.le.dpcoal.and.dr0.le.drcoal) IOVER(ISG)=1
1185 250 continue
1186c
1187C3 baryon q (antibaryon qbar) coalesce with all available q (qbar):
1188 do 350 ISG=1,NSG
1189 if(NJSGS(ISG).ne.3.or.IOVER(ISG).eq.1) goto 350
1190 ibaryn=K2SGS(ISG,1)
1191C DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1192 do 1004 j=1,2
1193 ftp(j)=ftsgs(isg,j)
1194 gxp(j)=gxsgs(isg,j)
1195 gyp(j)=gysgs(isg,j)
1196 gzp(j)=gzsgs(isg,j)
1197 pxp(j)=pxsgs(isg,j)
1198 pyp(j)=pysgs(isg,j)
1199 pzp(j)=pzsgs(isg,j)
1200 pmp(j)=pmsgs(isg,j)
1201 pep(j)=pesgs(isg,j)
1202 1004 continue
1203 call locldr(2,drlocl)
1204 dr1(2)=drlocl
1205 dp1(2)=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1206 & -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1207c
1208 ftp(2)=ftsgs(isg,3)
1209 gxp(2)=gxsgs(isg,3)
1210 gyp(2)=gysgs(isg,3)
1211 gzp(2)=gzsgs(isg,3)
1212 pxp(2)=pxsgs(isg,3)
1213 pyp(2)=pysgs(isg,3)
1214 pzp(2)=pzsgs(isg,3)
1215 pmp(2)=pmsgs(isg,3)
1216 pep(2)=pesgs(isg,3)
1217 call locldr(2,drlocl)
1218 dr1(3)=drlocl
1219 dp1(3)=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1220 & -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1221c
1222 do 320 JSG=1,NSG
1223 if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 320
1224 if(NJSGS(JSG).eq.2) then
1225 if(ibaryn.gt.0) then
1226 ipmin=1
1227 else
1228 ipmin=2
1229 endif
1230 ipmax=ipmin
1231 elseif(NJSGS(JSG).eq.3.and.
1232 1 (ibaryn*K2SGS(JSG,1)).gt.0) then
1233 ipmin=1
1234 ipmax=3
1235 else
1236 goto 320
1237 endif
1238 do 300 ip=ipmin,ipmax
1239 dplocl=dsqrt(2*(pep(1)*pesgs(jsg,ip)
1240 1 -pxp(1)*pxsgs(jsg,ip)
1241 2 -pyp(1)*pysgs(jsg,ip)
1242 3 -pzp(1)*pzsgs(jsg,ip)
1243 4 -pmp(1)*pmsgs(jsg,ip)))
1244c skip if outside of momentum radius:
1245 if(dplocl.gt.dpcoal) goto 320
1246 ftp(2)=ftsgs(jsg,ip)
1247 gxp(2)=gxsgs(jsg,ip)
1248 gyp(2)=gysgs(jsg,ip)
1249 gzp(2)=gzsgs(jsg,ip)
1250 pxp(2)=pxsgs(jsg,ip)
1251 pyp(2)=pysgs(jsg,ip)
1252 pzp(2)=pzsgs(jsg,ip)
1253 pmp(2)=pmsgs(jsg,ip)
1254 pep(2)=pesgs(jsg,ip)
1255 call locldr(2,drlocl)
1256c skip if outside of spatial radius:
1257 if(drlocl.gt.drcoal) goto 320
1258c q_isg may coalesce with q_jsg for a baryon:
1259 ipi=0
1260 if(dp1(2).gt.dpcoal.or.dr1(2).gt.drcoal) then
1261 ipi=2
1262 if((dp1(3).gt.dpcoal.or.dr1(3).gt.drcoal)
1263 1 .and.dr1(3).gt.dr1(2)) ipi=3
1264 elseif(dp1(3).gt.dpcoal.or.dr1(3).gt.drcoal) then
1265 ipi=3
1266 elseif(dr1(2).lt.dr1(3)) then
1267 if(drlocl.lt.dr1(3)) ipi=3
1268 elseif(dr1(3).le.dr1(2)) then
1269 if(drlocl.lt.dr1(2)) ipi=2
1270 endif
1271 if(ipi.ne.0) then
1272 dp1(ipi)=dplocl
1273 dr1(ipi)=drlocl
1274 call exchge(isg,ipi,jsg,ip)
1275 endif
1276 300 continue
1277 320 continue
1278 if(dp1(2).le.dpcoal.and.dr1(2).le.drcoal
1279 1 .and.dp1(3).le.dpcoal.and.dr1(3).le.drcoal)
1280 2 IOVER(ISG)=1
1281 350 continue
1282c
1283 RETURN
1284 END
1285
1286c=======================================================================
1287 SUBROUTINE exchge(isg,ipi,jsg,ipj)
1288c
1289 implicit double precision (a-h, o-z)
1290 PARAMETER (MAXSTR=150001)
1291 COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
1292 & PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
1293 & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
1294 & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
1295cc SAVE /SOFT/
1296 SAVE
1297c
1298 k1=K1SGS(isg,ipi)
1299 k2=K2SGS(isg,ipi)
1300 px=PXSGS(isg,ipi)
1301 py=PYSGS(isg,ipi)
1302 pz=PZSGS(isg,ipi)
1303 pe=PESGS(isg,ipi)
1304 pm=PMSGS(isg,ipi)
1305 gx=GXSGS(isg,ipi)
1306 gy=GYSGS(isg,ipi)
1307 gz=GZSGS(isg,ipi)
1308 ft=FTSGS(isg,ipi)
1309 K1SGS(isg,ipi)=K1SGS(jsg,ipj)
1310 K2SGS(isg,ipi)=K2SGS(jsg,ipj)
1311 PXSGS(isg,ipi)=PXSGS(jsg,ipj)
1312 PYSGS(isg,ipi)=PYSGS(jsg,ipj)
1313 PZSGS(isg,ipi)=PZSGS(jsg,ipj)
1314 PESGS(isg,ipi)=PESGS(jsg,ipj)
1315 PMSGS(isg,ipi)=PMSGS(jsg,ipj)
1316 GXSGS(isg,ipi)=GXSGS(jsg,ipj)
1317 GYSGS(isg,ipi)=GYSGS(jsg,ipj)
1318 GZSGS(isg,ipi)=GZSGS(jsg,ipj)
1319 FTSGS(isg,ipi)=FTSGS(jsg,ipj)
1320 K1SGS(jsg,ipj)=k1
1321 K2SGS(jsg,ipj)=k2
1322 PXSGS(jsg,ipj)=px
1323 PYSGS(jsg,ipj)=py
1324 PZSGS(jsg,ipj)=pz
1325 PESGS(jsg,ipj)=pe
1326 PMSGS(jsg,ipj)=pm
1327 GXSGS(jsg,ipj)=gx
1328 GYSGS(jsg,ipj)=gy
1329 GZSGS(jsg,ipj)=gz
1330 FTSGS(jsg,ipj)=ft
1331c
1332 RETURN
1333 END
1334
1335c=======================================================================
1336 SUBROUTINE locldr(icall,drlocl)
1337c
1338 implicit double precision (a-h, o-z)
1339 dimension ftp0(3),pxp0(3),pyp0(3),pzp0(3),pep0(3)
1340 common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1341 1 pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
1342cc SAVE /loclco/
1343 common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
1344cc SAVE /prtn23/
3006c44b 1345 common /lora/ enenew, pxnew, pynew, pznew
1346cc SAVE /lora/
0119ef9a 1347 SAVE
1348c for 2-body kinematics:
1349 if(icall.eq.2) then
1350 etot=pep(1)+pep(2)
1351 bex=(pxp(1)+pxp(2))/etot
1352 bey=(pyp(1)+pyp(2))/etot
1353 bez=(pzp(1)+pzp(2))/etot
1354c boost the reference frame down by beta to get to the pair rest frame:
1355 do 1001 j=1,2
1356 beta2 = bex ** 2 + bey ** 2 + bez ** 2
1357 gam = 1.d0 / dsqrt(1.d0 - beta2)
1358 if(beta2.ge.0.9999999999999d0) then
1359 write(6,*) '4',pxp(1),pxp(2),pyp(1),pyp(2),
1360 1 pzp(1),pzp(2),pep(1),pep(2),pmp(1),pmp(2),
1361 2 dsqrt(pxp(1)**2+pyp(1)**2+pzp(1)**2+pmp(1)**2)/pep(1),
1362 3 dsqrt(pxp(1)**2+pyp(1)**2+pzp(1)**2)/pep(1)
1363 write(6,*) '4a',pxp(1)+pxp(2),pyp(1)+pyp(2),
1364 1 pzp(1)+pzp(2),etot
1365 write(6,*) '4b',bex,bey,bez,beta2,gam
1366 endif
1367c
e2670fbe 1368 call lorenza(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez)
0119ef9a 1369 gxp0(j)=pxnew
1370 gyp0(j)=pynew
1371 gzp0(j)=pznew
1372 ftp0(j)=enenew
e2670fbe 1373 call lorenza(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez)
0119ef9a 1374 pxp0(j)=pxnew
1375 pyp0(j)=pynew
1376 pzp0(j)=pznew
1377 pep0(j)=enenew
1378 1001 continue
1379c
1380 if(ftp0(1).ge.ftp0(2)) then
1381 ilate=1
1382 iearly=2
1383 else
1384 ilate=2
1385 iearly=1
1386 endif
1387 ft0fom=ftp0(ilate)
1388c
1389 dt0=ftp0(ilate)-ftp0(iearly)
1390 gxp0(iearly)=gxp0(iearly)+pxp0(iearly)/pep0(iearly)*dt0
1391 gyp0(iearly)=gyp0(iearly)+pyp0(iearly)/pep0(iearly)*dt0
1392 gzp0(iearly)=gzp0(iearly)+pzp0(iearly)/pep0(iearly)*dt0
1393 drlocl=dsqrt((gxp0(ilate)-gxp0(iearly))**2
1394 1 +(gyp0(ilate)-gyp0(iearly))**2
1395 2 +(gzp0(ilate)-gzp0(iearly))**2)
1396c for 3-body kinematics, used for baryons formation:
1397 elseif(icall.eq.3) then
1398 etot=pep(1)+pep(2)+pep(3)
1399 bex=(pxp(1)+pxp(2)+pxp(3))/etot
1400 bey=(pyp(1)+pyp(2)+pyp(3))/etot
1401 bez=(pzp(1)+pzp(2)+pzp(3))/etot
1402 beta2 = bex ** 2 + bey ** 2 + bez ** 2
1403 gam = 1.d0 / dsqrt(1.d0 - beta2)
1404 if(beta2.ge.0.9999999999999d0) then
1405 write(6,*) '5',bex,bey,bez,beta2,gam
1406 endif
1407c boost the reference frame down by beta to get to the 3-parton rest frame:
1408 do 1002 j=1,3
e2670fbe 1409 call lorenza(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez)
0119ef9a 1410 gxp0(j)=pxnew
1411 gyp0(j)=pynew
1412 gzp0(j)=pznew
1413 ftp0(j)=enenew
e2670fbe 1414 call lorenza(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez)
0119ef9a 1415 pxp0(j)=pxnew
1416 pyp0(j)=pynew
1417 pzp0(j)=pznew
1418 pep0(j)=enenew
1419 1002 continue
1420c
1421 if(ftp0(1).gt.ftp0(2)) then
1422 ilate=1
1423 if(ftp0(3).gt.ftp0(1)) ilate=3
1424 else
1425 ilate=2
1426 if(ftp0(3).ge.ftp0(2)) ilate=3
1427 endif
1428 ft0fom=ftp0(ilate)
1429c
1430 if(ilate.eq.1) then
1431 imin=2
1432 imax=3
1433 istep=1
1434 elseif(ilate.eq.2) then
1435 imin=1
1436 imax=3
1437 istep=2
1438 elseif(ilate.eq.3) then
1439 imin=1
1440 imax=2
1441 istep=1
1442 endif
1443c
1444 do 1003 iearly=imin,imax,istep
1445 dt0=ftp0(ilate)-ftp0(iearly)
1446 gxp0(iearly)=gxp0(iearly)+pxp0(iearly)/pep0(iearly)*dt0
1447 gyp0(iearly)=gyp0(iearly)+pyp0(iearly)/pep0(iearly)*dt0
1448 gzp0(iearly)=gzp0(iearly)+pzp0(iearly)/pep0(iearly)*dt0
1449 1003 continue
1450 endif
1451c
1452 RETURN
1453 END
1454
1455c=======================================================================
1456 subroutine hoscar
1457c
1458 parameter (MAXSTR=150001,AMN=0.939457,AMP=0.93828)
1459 character*8 code, reffra, FRAME
1460 character*25 amptvn
1461 common/snn/efrm,npart1,npart2
1462cc SAVE /snn/
1463 common /lastt/itimeh,bimp
1464cc SAVE /lastt/
1465 COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
1466cc SAVE /hbt/
1467 common/oscar1/iap,izp,iat,izt
1468cc SAVE /oscar1/
1469 common/oscar2/FRAME,amptvn
1470cc SAVE /oscar2/
1471 SAVE
1472 data nff/0/
1473c
1474c file header
1475 if(nff.eq.0) then
1476 write (19, 101) 'OSCAR1997A'
1477 write (19, 111) 'final_id_p_x'
1478 code = 'AMPT'
1479 if(FRAME.eq.'CMS') then
1480 reffra = 'nncm'
1481 xmp=(amp*izp+amn*(iap-izp))/iap
1482 xmt=(amp*izt+amn*(iat-izt))/iat
1483 ebeam=(efrm**2-xmp**2-xmt**2)/2./xmt
1484 elseif(FRAME.eq.'LAB') then
1485 reffra = 'lab'
1486 ebeam=efrm
1487 else
1488 reffra = 'unknown'
1489 ebeam=0.
1490 endif
1491 ntestp = 1
1492 write (19, 102) code, amptvn, iap, izp, iat, izt,
1493 & reffra, ebeam, ntestp
1494 nff = 1
1495 ievent = 1
1496 phi = 0.
1497 if(FRAME.eq.'CMS') write(19,112) efrm
1498 endif
1499c comment
1500c event header
1501 write (19, 103) ievent, nlast, bimp, phi
1502c particles
1503 do 99 i = 1, nlast
1504 ene=sqrt(plast(1,i)**2+plast(2,i)**2+plast(3,i)**2
1505 1 +plast(4,i)**2)
1506 write (19, 104) i, INVFLV(lblast(i)), plast(1,i),
1507 1 plast(2,i),plast(3,i),ene,plast(4,i),
1508 2 xlast(1,i),xlast(2,i),xlast(3,i),xlast(4,i)
1509 99 continue
1510 ievent = ievent + 1
1511 101 format (a10)
1512 111 format (a12)
1513 102 format (a4,1x,a20,1x,'(', i3, ',', i3, ')+(', i3, ',',
1514 & i3, ')', 2x, a4, 2x, e10.4, 2x, i8)
1515 103 format (i10, 2x, i10, 2x, f8.3, 2x, f8.3)
1516 104 format (i10, 2x, i10, 2x, 9(e12.6, 2x))
1517 112 format ('# Center-of-mass energy/nucleon-pair is',
1518 & f12.3,'GeV')
1519c
1520 return
1521 end
1522
1523c=======================================================================
1524 subroutine getnp
1525
1526 PARAMETER (MAXSTR=150001)
1527 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
1528cc SAVE /HMAIN1/
1529 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
1530cc SAVE /HMAIN2/
1531 COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50)
1532cc SAVE /HPARNT/
1533 common/snn/efrm,npart1,npart2
1534cc SAVE /snn/
1535cwei DOUBLE PRECISION PATT
1536 COMMON/EXTRACENT/BB
1537 SAVE
1538
1539 BB = HINT1(19)
1540 if(NATT.eq.0) then
1541 npart1=0
1542 npart2=0
1543 return
1544 endif
1545c
1546 PZPROJ=SQRT(HINT1(6)**2-HINT1(8)**2)
1547 PZTARG=SQRT(HINT1(7)**2-HINT1(9)**2)
1548 epsiln=0.01
1549c
1550 nspec1=0
1551 nspec2=0
1552 DO 1000 I = 1, NATT
1553 if((KATT(I,1).eq.2112.or.KATT(I,1).eq.2212)
1554 1 .and.PATT(I, 1).eq.0.and.PATT(I, 2).eq.0) then
1555 if(PATT(I, 3).gt.amax1(0.,PZPROJ-epsiln)) then
1556 nspec1=nspec1+1
1557 elseif(PATT(I, 3).lt.(-PZTARG+epsiln)) then
1558 nspec2=nspec2+1
1559 endif
1560 endif
1561 1000 CONTINUE
1562 npart1=IHNT2(1)-nspec1
1563 npart2=IHNT2(3)-nspec2
1564
1565 return
1566 end
1567
1568c=======================================================================
1569c 2/18/03 use PYTHIA to decay eta,rho,omega,k*,phi and Delta
1570 subroutine resdec(i1,nt,nnn,wid,idecay)
1571
1572 PARAMETER (hbarc=0.19733)
1573 PARAMETER (AK0=0.498,APICH=0.140,API0=0.135,AN=0.940,ADDM=0.02)
1574 PARAMETER (MAXSTR=150001, MAXR=1)
1575 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1576 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1577cc SAVE /INPUT2/
1578 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1579cc SAVE /LUJETSA/
1580 COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1581cc SAVE /LUDAT1A/
1582 COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
1583cc SAVE /LUDAT2A/
1584 COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
1585cc SAVE /LUDAT3A/
1586 COMMON /CC/ E(MAXSTR)
1587cc SAVE /CC/
1588 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1589cc SAVE /EE/
1590 COMMON /PA/RPION(3,MAXSTR,MAXR)
1591cc SAVE /PA/
1592 COMMON /PB/PPION(3,MAXSTR,MAXR)
1593cc SAVE /PB/
1594 COMMON /PC/EPION(MAXSTR,MAXR)
1595cc SAVE /PC/
1596 COMMON /PD/LPION(MAXSTR,MAXR)
1597cc SAVE /PD/
1598 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1599cc SAVE /input1/
1600 common/resdcy/NSAV,iksdcy
1601cc SAVE /resdcy/
1602 common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
1603 1 px1n,py1n,pz1n,dp1n
1604cc SAVE /leadng/
1605 EXTERNAL IARFLV, INVFLV
1606 COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
1607cc SAVE /tdecay/
1608 COMMON/RNDF77/NSEED
1609 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1610 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1611 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1612cc SAVE /RNDF77/
1613 SAVE
1614 irun=idecay
1615 if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1616 & .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
05cdcf94 1617c & .or.lb1.eq.28.or.lb1.eq.29
0119ef9a 1618 & .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1619 & .or.iabs(lb1).eq.16) then
1620 kf=INVFLV(lb1)
1621 else
1622 return
1623 endif
1624c
1625 IP=1
1626c label as undecayed and the only particle in the record:
1627 N=1
1628 K(IP,1)=1
1629 K(IP,3)=0
1630 K(IP,4)=0
1631 K(IP,5)=0
1632c
1633 K(IP,2)=kf
1634 P(IP,1)=px1
1635 P(IP,2)=py1
1636 P(IP,3)=pz1
1637 em1a=em1
1638c eta or omega in ART may be below or too close to (pi+pi-pi0) mass,
1639c causing LUDECY error,thus increase their mass ADDM above this thresh,
1640c noting that rho (m=0.281) too close to 2pi thrshold fails to decay:
1641 if((lb1.eq.0.or.lb1.eq.28).and.em1.lt.(2*APICH+API0+ADDM)) then
1642 em1=2*APICH+API0+ADDM
1643c rho
1644 elseif(lb1.ge.25.and.lb1.le.27.and.em1.lt.(2*APICH+ADDM)) then
1645 em1=2*APICH+ADDM
1646c K*
1647 elseif(iabs(lb1).eq.30.and.em1.lt.(APICH+AK0+ADDM)) then
1648 em1=APICH+AK0+ADDM
1649c Delta created in ART may be below (n+pich) mass, causing LUDECY error:
1650 elseif(iabs(lb1).ge.6.and.iabs(lb1).le.9
1651 1 .and.em1.lt.(APICH+AN+ADDM)) then
1652 em1=APICH+AN+ADDM
1653 endif
1654c if(em1.ge.(em1a+0.01)) write (6,*)
1655c 1 'Mass increase in resdec():',nt,em1-em1a,lb1
1656 e1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1657 P(IP,4)=e1
1658 P(IP,5)=em1
1659clin-5/2008:
1660 dpdecp=dpertp(i1)
1661 call ludecy(IP)
1662c add decay time to daughter's formation time at the last timestep:
1663 if(nt.eq.ntmax) then
1664 tau0=hbarc/wid
1665 taudcy=tau0*(-1.)*alog(1.-RANART(NSEED))
1666 ndaut=n-nsav
1667 if(ndaut.le.1) then
1668 write(10,*) 'note: ndaut(<1)=',ndaut
1669 write(89,*) 'note: ndaut(<1)=',ndaut
1670 call lulist(2)
1671 stop
1672 endif
1673c lorentz boost:
1674 taudcy=taudcy*e1/em1
1675 tfnl=tfnl+taudcy
1676 xfnl=xfnl+px1/e1*taudcy
1677 yfnl=yfnl+py1/e1*taudcy
1678 zfnl=zfnl+pz1/e1*taudcy
1679c at the last timestep, assign rho, K0S or eta (decay daughter)
1680c to lb(i1) only (not to lpion) in order to decay them again:
1681 if(n.ge.(nsav+2)) then
1682 do 1001 idau=nsav+2,n
1683 kdaut=K(idau,2)
1684 if(kdaut.eq.221.or.kdaut.eq.113
1685 1 .or.kdaut.eq.213.or.kdaut.eq.-213
1686 2 .or.kdaut.eq.310) then
1687c switch idau and i1(nsav+1):
1688 ksave=kdaut
1689 pxsave=p(idau,1)
1690 pysave=p(idau,2)
1691 pzsave=p(idau,3)
1692 esave=p(idau,4)
1693 xmsave=p(idau,5)
1694 K(idau,2)=K(nsav+1,2)
1695 p(idau,1)=p(nsav+1,1)
1696 p(idau,2)=p(nsav+1,2)
1697 p(idau,3)=p(nsav+1,3)
1698 p(idau,4)=p(nsav+1,4)
1699 p(idau,5)=p(nsav+1,5)
1700 K(nsav+1,2)=ksave
1701 p(nsav+1,1)=pxsave
1702 p(nsav+1,2)=pysave
1703 p(nsav+1,3)=pzsave
1704 p(nsav+1,4)=esave
1705 p(nsav+1,5)=xmsave
1706c note: phi decay may produce rho, K0s or eta, N*(1535) decay may produce
1707c eta, but only one daughter may be rho, K0s or eta:
1708 goto 111
1709 endif
1710 1001 continue
1711 endif
1712 111 continue
1713c
1714 enet=0.
1715 do 1002 idau=nsav+1,n
1716 enet=enet+p(idau,4)
1717 1002 continue
1718c if(abs(enet-e1).gt.0.02)
1719c 1 write(93,*) 'resdec(): nt=',nt,enet-e1,lb1
1720 endif
1721
1722cyy 200 format(a20,3(1x,i6))
1723cyy 210 format(i6,5(1x,f8.3))
1724cyy 220 format(a2,i5,5(1x,f8.3))
1725
1726 do 1003 idau=nsav+1,n
1727 kdaut=K(idau,2)
1728 lbdaut=IARFLV(kdaut)
1729c K0S and K0L are named K+/K- during hadron cascade, and only
1730c at the last timestep they keep their real LB # before output;
1731c K0/K0bar (from K* decay) converted to K0S and K0L at the last timestep:
1732 if(nt.eq.ntmax.and.(kdaut.eq.130.or.kdaut.eq.310
1733 1 .or.iabs(kdaut).eq.311)) then
1734 if(kdaut.eq.130) then
1735 lbdaut=22
1736 elseif(kdaut.eq.310) then
1737 lbdaut=24
1738 elseif(iabs(kdaut).eq.311) then
1739 if(RANART(NSEED).lt.0.5) then
1740 lbdaut=22
1741 else
1742 lbdaut=24
1743 endif
1744 endif
1745 endif
1746c
1747 if(idau.eq.(nsav+1)) then
1748 LB(i1)=lbdaut
1749 E(i1)=p(idau,5)
1750 px1n=p(idau,1)
1751 py1n=p(idau,2)
1752 pz1n=p(idau,3)
1753clin-5/2008:
1754 dp1n=dpdecp
1755 else
1756 nnn=nnn+1
1757 LPION(NNN,IRUN)=lbdaut
1758 EPION(NNN,IRUN)=p(idau,5)
1759 PPION(1,NNN,IRUN)=p(idau,1)
1760 PPION(2,NNN,IRUN)=p(idau,2)
1761 PPION(3,NNN,IRUN)=p(idau,3)
1762 RPION(1,NNN,IRUN)=xfnl
1763 RPION(2,NNN,IRUN)=yfnl
1764 RPION(3,NNN,IRUN)=zfnl
1765 tfdpi(NNN,IRUN)=tfnl
1766clin-5/2008:
1767 dppion(NNN,IRUN)=dpdecp
1768 endif
1769 1003 continue
1770cyy 230 format(a2,i5,5(1x,e8.2))
1771 return
1772 end
1773
1774c=======================================================================
1775 subroutine inidcy
1776
1777 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1778cc SAVE /LUJETSA/
1779 common/resdcy/NSAV,iksdcy
1780cc SAVE /resdcy/
1781 SAVE
1782 N=1
1783 NSAV=N
1784 return
1785 end
1786
1787c=======================================================================
1788clin-6/06/02 local parton freezeout motivated from critical density:
1789 subroutine local(t)
1790c
1791 implicit double precision (a-h, o-z)
1792 PARAMETER (MAXPTN=400001)
1793 PARAMETER (r0=1d0)
1794 COMMON /para1/ mul
1795cc SAVE /para1/
1796 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1797 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1798 & XMASS5(MAXPTN), ITYP5(MAXPTN)
1799cc SAVE /prec2/
1800 common /frzprc/
1801 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1802 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1803 & xmfrz(MAXPTN),
1804 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1805cc SAVE /frzprc/
1806 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1807cc SAVE /prec4/
1808 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1809cc SAVE /prec5/
1810 common /coal/dpcoal,drcoal,ecritl
1811cc SAVE /coal/
1812 SAVE
1813c
1814 do 1001 it=1,301
1815 if(t.ge.tfrz(it).and.t.lt.tfrz(it+1)) then
1816 if(it.eq.itlast) then
1817 return
1818 else
1819 itlast=it
1820 goto 50
1821 endif
1822 endif
1823 1001 continue
1824 write(1,*) 'local time out of range in LOCAL, stop',t,it
1825 stop
1826 50 continue
1827c
1828 do 200 ip=1,mul
1829c skip partons which have frozen out:
1830 if(ifrz(ip).eq.1) goto 200
1831 if(it.eq.301) then
1832c freezeout all the left partons beyond the time of 3000 fm:
1833 etcrit=1d6
1834 goto 150
1835 else
1836c freezeout when transverse energy density < etcrit:
1837 etcrit=(ecritl*2d0/3d0)
1838 endif
1839c skip partons which have not yet formed:
1840 if(t.lt.FT5(ip)) goto 200
1841 rap0=rap(ip)
1842 eta0=eta(ip)
1843 x0=GX5(ip)+vx(ip)*(t-FT5(ip))
1844 y0=GY5(ip)+vy(ip)*(t-FT5(ip))
1845 detdy=0d0
1846 do 100 itest=1,mul
1847c skip self and partons which have not yet formed:
1848 if(itest.eq.ip.or.t.lt.FT5(itest)) goto 100
1849 ettest=eta(itest)
1850 xtest=GX5(itest)+vx(itest)*(t-FT5(itest))
1851 ytest=GY5(itest)+vy(itest)*(t-FT5(itest))
1852 drt=sqrt((xtest-x0)**2+(ytest-y0)**2)
1853c count partons within drt<1 and -1<(eta-eta0)<1:
1854 if(dabs(ettest-eta0).le.1d0.and.drt.le.r0)
1855 1 detdy=detdy+dsqrt(PX5(itest)**2+PY5(itest)**2
1856 2 +XMASS5(itest)**2)*0.5d0
1857 100 continue
1858 detdy=detdy*(dcosh(eta0)**2)/(t*3.1416d0*r0**2*dcosh(rap0))
1859c when density is below critical density for phase transition, freeze out:
1860 150 if(detdy.le.etcrit) then
1861 ifrz(ip)=1
1862 idfrz(ip)=ITYP5(ip)
1863 pxfrz(ip)=PX5(ip)
1864 pyfrz(ip)=PY5(ip)
1865 pzfrz(ip)=PZ5(ip)
1866 efrz(ip)=E5(ip)
1867 xmfrz(ip)=XMASS5(ip)
1868 if(t.gt.FT5(ip)) then
1869 gxfrz(ip)=x0
1870 gyfrz(ip)=y0
1871 gzfrz(ip)=GZ5(ip)+vz(ip)*(t-FT5(ip))
1872 ftfrz(ip)=t
1873 else
1874c if this freezeout time < formation time, use formation time & positions.
1875c This ensures the recovery of default hadron when e_crit=infty:
1876 gxfrz(ip)=GX5(ip)
1877 gyfrz(ip)=GY5(ip)
1878 gzfrz(ip)=GZ5(ip)
1879 ftfrz(ip)=FT5(ip)
1880 endif
1881 endif
1882 200 continue
1883c
1884 return
1885 end
1886
1887c=======================================================================
1888clin-6/06/02 initialization for local parton freezeout
1889 subroutine inifrz
1890c
1891 implicit double precision (a-h, o-z)
1892 PARAMETER (MAXPTN=400001)
1893 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1894cc SAVE /ilist5/
1895 common /frzprc/
1896 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1897 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1898 & xmfrz(MAXPTN),
1899 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1900cc SAVE /frzprc/
1901 SAVE
1902c
1903c for freezeout time 0-10fm, use interval of 0.1fm;
1904c for 10-100fm, use interval of 1fm;
1905c for 100-1000fm, use interval of 10fm;
1906c for 1000-3000fm, use interval of 100fm:
1907 step1=0.1d0
1908 step2=1d0
1909 step3=10d0
1910 step4=100d0
1911c
1912 do 1001 it=1,101
1913 tfrz(it)=0d0+dble(it-1)*step1
1914 1001 continue
1915 do 1002 it=102,191
1916 tfrz(it)=10d0+dble(it-101)*step2
1917 1002 continue
1918 do 1003 it=192,281
1919 tfrz(it)=100d0+dble(it-191)*step3
1920 1003 continue
1921 do 1004 it=282,301
1922 tfrz(it)=1000d0+dble(it-281)*step4
1923 1004 continue
1924 tfrz(302)=tlarge
1925c
1926 return
1927 end
1928
1929clin-5/2009 ctest on v2 analysis
1930c=======================================================================
1931c idd=0,1,2,3 specifies different subroutines for partonic flow analysis.
1932 subroutine flowp(idd)
1933c
1934 implicit double precision (a-h, o-z)
1935 real dt
1936 parameter (MAXPTN=400001)
1937csp
1938 parameter (bmt=0.05d0)
1939 dimension nlfile(3),nsfile(3),nmfile(3)
1940c
1941 dimension v2pp(3),xnpp(3),v2psum(3),v2p2sm(3),nfile(3)
1942 dimension tsp(31),v2pevt(3),v2pavg(3),varv2p(3)
1943 common /ilist1/
1944 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1945 & ictype, icsta(MAXPTN),
1946 & nic(MAXPTN), icels(MAXPTN)
1947cc SAVE /ilist1/
1948 COMMON /para1/ mul
1949cc SAVE /para1/
1950 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1951 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1952 & XMASS5(MAXPTN), ITYP5(MAXPTN)
1953cc SAVE /prec2/
1954 COMMON /pflow/ v2p(30,3),xnpart(30,3),etp(30,3),
1955 1 s2p(30,3),v2p2(30,3),nevt(30)
1956cc SAVE /pflow/
1957 COMMON /pflowf/ v2pf(30,3),xnpf(30,3),etpf(30,3),
1958 1 xncoll(30),s2pf(30,3),v2pf2(30,3)
1959cc SAVE /pflowf/
1960 COMMON /pfrz/ v2pfrz(30,3),xnpfrz(30,3),etpfrz(30,3),
1961 1 s2pfrz(30,3),v2p2fz(30,3),tscatt(31),
1962 2 nevtfz(30),iscatt(30)
1963cc SAVE /pfrz/
1964 COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
1965 1 v2h2(30,3),s2h(30,3)
1966cc SAVE /hflow/
1967 COMMON /AREVT/ IAEVT, IARUN, MISS
1968cc SAVE /AREVT/
1969 common/anim/nevent,isoft,isflag,izpc
1970cc SAVE /anim/
1971 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1972cc SAVE /input1/
1973 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1974 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1975cc SAVE /INPUT2/
1976cc SAVE itimep,iaevtp,v2pp,xnpp,v2psum,v2p2sm
1977cc SAVE nfile,itanim,nlfile,nsfile,nmfile
1978 SAVE
1979csp
1980 dimension etpl(30,3),etps(30,3),etplf(30,3),etpsf(30,3),
1981 & etlfrz(30,3),etsfrz(30,3),
1982 & xnpl(30,3),xnps(30,3),xnplf(30,3),xnpsf(30,3),
1983 & xnlfrz(30,3),xnsfrz(30,3),
1984 & v2pl(30,3),v2ps(30,3),v2plf(30,3),v2psf(30,3),
1985 & s2pl(30,3),s2ps(30,3),s2plf(30,3),s2psf(30,3),
1986 & DMYil(50,3),DMYfl(50,3),
1987 & DMYis(50,3),DMYfs(50,3)
1988 data tsp/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
1989 & 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
1990 & 2 , 3, 4, 5, 6, 7, 8, 9, 10, 20, 30/
1991c idd=0: initialization for flow analysis, called by artdri.f:
1992 if(idd.eq.0) then
1993 nfile(1)=60
1994 nfile(2)=64
1995 nfile(3)=20
1996cms OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN')
1997cms OPEN (nfile(1)+1,
1998cms 1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN')
1999cms OPEN (nfile(1)+2,
2000cms 1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN')
2001cms OPEN (nfile(1)+3,
2002cms 1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN')
2003cms OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN')
2004cms OPEN (nfile(2)+1,
2005cms 1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN')
2006cms OPEN (nfile(2)+2,
2007cms 1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN')
2008cms OPEN (nfile(2)+3,
2009cms 1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN')
2010cms OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN')
2011cms OPEN (nfile(3)+1,
2012cms 1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN')
2013cms OPEN (nfile(3)+2,
2014cms 1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN')
2015cms OPEN (nfile(3)+3,
2016cms 1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN')
2017cms OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN')
2018cms write(49, *) ' ievt, v2p, v2p_y2, v2p_y1'
2019c
2020cms OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN')
2021cms OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN')
2022cms OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN')
2023cms OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN')
2024cms write(88, *) ' ievt, v2h, v2h_y2, v2h_y1'
2025csp07/05
2026 nlfile(1)=70
2027 nlfile(2)=72
2028 nlfile(3)=74
2029cms OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN')
2030cms OPEN (nlfile(1)+1,
2031cms 1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN')
2032cms OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN')
2033cms OPEN (nlfile(2)+1,
2034cms 1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN')
2035cms OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN')
2036cms OPEN (nlfile(3)+1,
2037cms 1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN')
2038 nsfile(1)=76
2039 nsfile(2)=78
2040 nsfile(3)=80
2041cms OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN')
2042cms OPEN (nsfile(1)+1,
2043cms 1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN')
2044cms OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN')
2045cms OPEN (nsfile(2)+1,
2046cms 1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN')
2047cms OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN')
2048cms OPEN (nsfile(3)+1,
2049cms 1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN')
2050 nmfile(1)=82
2051 nmfile(2)=83
2052 nmfile(3)=84
2053cms OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN')
2054cms OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN')
2055cms OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN')
2056clin-11/27/00 for animation:
2057 if(nevent.eq.1) then
2058cms OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN')
2059cms write(91,*) ntmax, dt
2060cms OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN')
2061cms OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN')
2062 endif
2063c
2064 itimep=0
2065 itanim=0
2066 iaevtp=0
2067csp
2068 do 1002 ii=1,50
2069 do 1001 iy=1,3
2070 DMYil(ii,iy) = 0d0
2071 DMYfl(ii,iy) = 0d0
2072 DMYis(ii,iy) = 0d0
2073 DMYfs(ii,iy) = 0d0
2074 1001 continue
2075 1002 continue
2076c
2077 do 1003 ii=1,31
2078 tscatt(ii)=0d0
2079 1003 continue
2080 do 1005 ii=1,30
2081 nevt(ii)=0
2082 xncoll(ii)=0d0
2083 nevtfz(ii)=0
2084 iscatt(ii)=0
2085 do 1004 iy=1,3
2086 v2p(ii,iy)=0d0
2087 v2p2(ii,iy)=0d0
2088 s2p(ii,iy)=0d0
2089 etp(ii,iy)=0d0
2090 xnpart(ii,iy)=0d0
2091 v2pf(ii,iy)=0d0
2092 v2pf2(ii,iy)=0d0
2093 s2pf(ii,iy)=0d0
2094 etpf(ii,iy)=0d0
2095 xnpf(ii,iy)=0d0
2096 v2pfrz(ii,iy)=0d0
2097 v2p2fz(ii,iy)=0d0
2098 s2pfrz(ii,iy)=0d0
2099 etpfrz(ii,iy)=0d0
2100 xnpfrz(ii,iy)=0d0
2101csp07/05
2102 etpl(ii,iy)=0d0
2103 etps(ii,iy)=0d0
2104 etplf(ii,iy)=0d0
2105 etpsf(ii,iy)=0d0
2106 etlfrz(ii,iy)=0d0
2107 etsfrz(ii,iy)=0d0
2108 xnpl(ii,iy)=0d0
2109 xnps(ii,iy)=0d0
2110 xnplf(ii,iy)=0d0
2111 xnpsf(ii,iy)=0d0
2112 xnlfrz(ii,iy)=0d0
2113 xnsfrz(ii,iy)=0d0
2114 v2pl(ii,iy)=0d0
2115 v2ps(ii,iy)=0d0
2116 v2plf(ii,iy)=0d0
2117 v2psf(ii,iy)=0d0
2118 s2pl(ii,iy)=0d0
2119 s2ps(ii,iy)=0d0
2120 s2plf(ii,iy)=0d0
2121 s2psf(ii,iy)=0d0
2122 1004 continue
2123 1005 continue
2124 do 1006 iy=1,3
2125 v2pevt(iy)=0d0
2126 v2pavg(iy)=0d0
2127 varv2p(iy)=0d0
2128 v2pp(iy)=0.d0
2129 xnpp(iy)=0d0
2130 v2psum(iy)=0.d0
2131 v2p2sm(iy)=0.d0
2132 1006 continue
2133c idd=1: calculate parton elliptic flow, called by zpc.f:
2134 else if(idd.eq.1) then
2135 t2time = FT5(iscat)
2136 do 1008 ianp = 1, 30
2137 if (t2time.lt.tsp(ianp+1).and.t2time.ge.tsp(ianp)) then
2138c write flow info only once at each fixed time:
2139 xncoll(ianp)=xncoll(ianp)+1d0
2140c to prevent an earlier t2time comes later in the same event
2141c and mess up nevt:
2142 if(ianp.le.itimep.and.iaevt.eq.iaevtp) goto 101
2143 nevt(ianp)=nevt(ianp)+1
2144 tscatt(ianp+1)=t2time
2145 iscatt(ianp)=1
2146 nevtfz(ianp)=nevtfz(ianp)+1
2147 do 100 I=1,mul
2148 ypartn=0.5d0*dlog((E5(i)+PZ5(i))/(E5(i)-PZ5(i)+1.d-8))
2149 pt2=PX5(I)**2+PY5(I)**2
2150ctest off: initial (pt,y) and (x,y) distribution:
2151c idtime=1
2152c if(ianp.eq.idtime) then
2153c iityp=iabs(ITYP5(I))
2154c if(iityp.eq.1.or.iityp.eq.2) then
2155c write(651,*) dsqrt(pt2),ypartn
2156c write(654,*) GX5(I),GY5(I)
2157c elseif(iityp.eq.1103.or.iityp.eq.2101
2158c 1 .or.iityp.eq.2103.or.iityp.eq.2203.
2159c 2 .or.iityp.eq.3101.or.iityp.eq.3103.
2160c 3 .or.iityp.eq.3201.or.iityp.eq.3203.or.iityp.eq.3303)
2161c 4 then
2162c write(652,*) dsqrt(pt2),ypartn
2163c write(655,*) GX5(I),GY5(I)
2164c elseif(iityp.eq.21) then
2165c write(653,*) dsqrt(pt2),ypartn
2166c write(656,*) GX5(I),GY5(I)
2167c endif
2168c endif
2169ctest-end
2170ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2171c gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2172c gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2173c gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2174c rt_now=dsqrt(gx_now**2+gy_now**2)
2175c zmax=0.1d0*t2time
2176c volume=3.1416d0*(2d0**2)*(2*zmax)
2177c if(rt_now.gt.2d0.or.dabs(gz_now).gt.zmax)
2178c 1 goto 100
2179ctest-end
2180 iloop=1
2181 if(dabs(ypartn).le.1d0) then
2182 iloop=2
2183 if(dabs(ypartn).le.0.5d0) then
2184 iloop=3
2185 endif
2186 endif
2187 do 50 iy=1,iloop
2188 if(pt2.gt.0.) then
2189 v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2190 if(abs(v2prtn).gt.1.)
2191 1 write(nfile(iy),*) 'v2prtn>1',v2prtn
2192 v2p(ianp,iy)=v2p(ianp,iy)+v2prtn
2193 v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2
2194 endif
2195 xperp2=GX5(I)**2+GY5(I)**2
2196 if(xperp2.gt.0.)
2197 1 s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2198 xnpart(ianp,iy)=xnpart(ianp,iy)+1d0
2199 etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2200ctest off density:
2201c etp(ianp,iy)=etp(ianp,iy)
2202c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2203clin-2/22/00 to write out parton info only for formed ones:
2204 if(FT5(I).le.t2time) then
2205 v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn
2206 v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2
2207 if(xperp2.gt.0.)
2208 1 s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2209 xnpf(ianp,iy)=xnpf(ianp,iy)+1d0
2210 etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2211ctest off density:
2212c etpf(ianp,iy)=etpf(ianp,iy)
2213c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2214 endif
2215 50 continue
2216 100 continue
2217 itimep=ianp
2218 iaevtp=iaevt
2219clin-3/30/00 ebe v2 variables:
2220 if(ianp.eq.30) then
2221 do 1007 iy=1,3
2222 npartn=IDINT(xnpart(ianp,iy)-xnpp(iy))
2223 if(npartn.ne.0) then
2224 v2pevt(iy)=(v2p(ianp,iy)-v2pp(iy))/npartn
2225 v2psum(iy)=v2psum(iy)+v2pevt(iy)
2226 v2p2sm(iy)=v2p2sm(iy)+v2pevt(iy)**2
2227 v2pp(iy)=v2p(ianp,iy)
2228 xnpp(iy)=xnpart(ianp,iy)
2229 endif
2230 1007 continue
2231 write(49, 160) iaevt,v2pevt
2232 endif
2233 goto 101
2234 endif
2235 1008 continue
2236clin-11/28/00 for animation:
2237 101 if(nevent.eq.1) then
2238 do 110 nt = 1, ntmax
2239 time1=dble(nt*dt)
2240 time2=dble((nt+1)*dt)
2241 if (t2time.lt.time2.and.t2time.ge.time1) then
2242 if(nt.le.itanim) return
2243c write(92,*) t2time
2244 iform=0
2245 do 1009 I=1,mul
2246c write out parton info only for formed ones:
2247 if(FT5(I).le.t2time) then
2248 iform=iform+1
2249 endif
2250 1009 continue
2251c write(92,*) iform
2252 do 120 I=1,mul
2253 if(FT5(I).le.t2time) then
2254clin-11/29/00-ctest off calculate parton coordinates after propagation:
2255c gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2256c gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2257c gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2258c write(92,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2259c write(92,180) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i),
2260c 1 PX5(i),PY5(i),PZ5(i),E5(i)
2261ctest-end
2262 endif
2263 120 continue
2264 itanim=nt
2265 endif
2266 110 continue
2267 endif
2268c
2269 140 format(i10,4(2x,f7.2))
2270 160 format(i10,3(2x,f9.5))
2271c 180 format(i6,8(1x,f7.2))
2272clin-5/17/01 calculate v2 for active partons (which have not frozen out):
2273c idd=3, called at end of zpc.f:
2274 else if(idd.eq.3) then
2275 do 1010 ianp=1,30
2276 if(iscatt(ianp).eq.0) tscatt(ianp+1)=tscatt(ianp)
2277 1010 continue
2278 do 350 I=1,mul
2279 ypartn=0.5d0*dlog((E5(i)+PZ5(i)+1.d-8)
2280 1 /(E5(i)-PZ5(i)+1.d-8))
2281 pt2=PX5(I)**2+PY5(I)**2
2282 iloop=1
2283 if(dabs(ypartn).le.1d0) then
2284 iloop=2
2285 if(dabs(ypartn).le.0.5d0) then
2286 iloop=3
2287 endif
2288 endif
2289c
2290 do 325 ianp=1,30
2291 if(iscatt(ianp).ne.0) then
2292 if(FT5(I).lt.tscatt(ianp+1)
2293 1 .and.FT5(I).ge.tscatt(ianp)) then
2294 do 1011 iy=1,iloop
2295 if(pt2.gt.0.) then
2296 v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2297 v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn
2298 v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2
2299 endif
2300 xperp2=GX5(I)**2+GY5(I)**2
2301 if(xperp2.gt.0.) s2pfrz(ianp,iy)=
2302 1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2303 etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2304 xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0
2305ctest off density:
2306c etpfrz(ianp,iy)=etpfrz(ianp,iy)
2307c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2308csp07/05
2309 if(ITYP5(I).eq.1.or.ITYP5(I).eq.2)then
2310 etlfrz(ianp,iy)=etlfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2311 xnlfrz(ianp,iy)=xnlfrz(ianp,iy)+1d0
2312 elseif(ITYP5(I).eq.3)then
2313 etsfrz(ianp,iy)=etsfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2314 xnsfrz(ianp,iy)=xnsfrz(ianp,iy)+1d0
2315 endif
2316csp07/05 end
2317 1011 continue
2318c parton freezeout info taken, proceed to next parton:
2319 goto 350
2320 endif
2321 endif
2322 325 continue
2323 350 continue
2324c idd=2: calculate average partonic elliptic flow, called from artdri.f,
2325 else if(idd.eq.2) then
2326 do 1012 i=1,3
2327cms write(nfile(i),*) ' tsp, v2p, v2p2, '//
2328cms 1 ' s2p, etp, xmult, nevt, xnctot'
2329cms write ((nfile(i)+1),*) ' tsp, v2pf, v2pf2, '//
2330cms 1 ' s2pf, etpf, xnform, xcrate'
2331cms write ((nfile(i)+2),*) ' tsp, v2pa, v2pa2, '//
2332cms 1 ' s2pa, etpa, xmult_ap, xnform, nevt'
2333cms write ((nfile(i)+3),*) ' tsph, v2ph, v2ph2, '//
2334cms 1 ' s2ph, etph, xmult_(ap/2+h),xmult_ap/2,nevt'
2335csp
2336cms write(nlfile(i),*) ' tsp, v2, s2, etp, xmul'
2337cms write(nsfile(i),*) ' tsp, v2, s2, etp, xmul'
2338cms write(nlfile(i)+1,*) ' tsp, v2, s2, etp, xmul'
2339cms write(nsfile(i)+1,*) ' tsp, v2, s2, etp, xmul'
2340c
2341 1012 continue
2342clin-3/30/00 ensemble average & variance of v2 (over particles in all events):
2343 do 150 ii=1, 30
2344 if(nevt(ii).eq.0) goto 150
2345 do 1014 iy=1,3
2346 if(xnpart(ii,iy).gt.1) then
2347 v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy)
2348 v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy)
2349 1 -v2p(ii,iy)**2)/(xnpart(ii,iy)-1))
2350 s2p(ii,iy)=s2p(ii,iy)/xnpart(ii,iy)
2351c xmult and etp are multiplicity and et for an averaged event:
2352 xmult=dble(xnpart(ii,iy)/dble(nevt(ii)))
2353 etp(ii,iy)=etp(ii,iy)/dble(nevt(ii))
2354csp
2355 etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2356 etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2357c
2358 xnctot=0d0
2359 do 1013 inum=1,ii
2360 xnctot=xnctot+xncoll(inum)
2361 1013 continue
2362 if(nevt(1).ne.0) xnctot=xnctot/nevt(1)
2363cms write (nfile(iy),200) tsp(ii),v2p(ii,iy),
2364cms 1 v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot
2365 endif
2366 if(nevt(ii).ne.0)
2367 1 xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii)
2368c
2369 if(xnpf(ii,iy).gt.1) then
2370 v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy)
2371 v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy)
2372 1 -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1))
2373 s2pf(ii,iy)=s2pf(ii,iy)/xnpf(ii,iy)
2374 xnform=dble(xnpf(ii,iy)/dble(nevt(ii)))
2375 etpf(ii,iy)=etpf(ii,iy)/dble(nevt(ii))
2376csp
2377 etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2378 etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2379c
2380cms write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy),
2381cms 1 v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate
2382 endif
2383csp
2384 if(xnpl(ii,iy).gt.1) then
2385 v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy)
2386 s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy)
2387 xmult=dble(xnpl(ii,iy)/dble(nevt(ii)))
2388 etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2389cms write (nlfile(iy),201) tsp(ii),v2pl(ii,iy),
2390cms 1 s2pl(ii,iy),etpl(ii,iy),xmult
2391 endif
2392 if(xnps(ii,iy).gt.1) then
2393 v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy)
2394 s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy)
2395 xmult=dble(xnps(ii,iy)/dble(nevt(ii)))
2396 etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2397cms write (nsfile(iy),201) tsp(ii),v2ps(ii,iy),
2398cms 1 s2ps(ii,iy),etps(ii,iy),xmult
2399 endif
2400 if(xnplf(ii,iy).gt.1) then
2401 v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy)
2402 s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy)
2403 xmult=dble(xnplf(ii,iy)/dble(nevt(ii)))
2404 etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2405cms write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy),
2406cms 1 s2plf(ii,iy),etplf(ii,iy),xmult
2407 endif
2408 if(xnpsf(ii,iy).gt.1) then
2409 v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy)
2410 s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy)
2411 xmult=dble(xnpsf(ii,iy)/dble(nevt(ii)))
2412 etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2413cms write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy),
2414cms 1 s2psf(ii,iy),etpsf(ii,iy),xmult
2415 endif
2416csp-end
2417 1014 continue
2418 150 continue
2419csp07/05 initial & final mt distrb
2420 scalei=0d0
2421 scalef=0d0
2422 if(nevt(1).ne.0) SCALEi = 1d0 / dble(nevt(1)) / BMT
2423 if(nevt(30).ne.0) SCALEf = 1d0 / dble(nevt(30)) / BMT
2424 do 1016 iy=2,3
2425 yra = 1d0
2426 if(iy .eq. 2)yra = 2d0
2427 do 1015 i=1,50
2428cms WRITE(nmfile(iy),251) BMT*dble(I - 0.5),
2429cms & SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra,
2430cms & SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra
2431 1015 continue
2432 1016 continue
2433csp07/05 end
2434clin-3/30/00 event-by-event average & variance of v2:
2435 if(nevt(30).ge.1) then
2436 do 1017 iy=1,3
2437 v2pavg(iy)=v2psum(iy)/nevt(30)
2438 v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2
2439 if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0)
2440 1017 continue
2441 write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg
2442 write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p
2443 endif
2444clin-11/28/00 for animation:
2445 if(nevent.eq.1) then
2446 do 1018 I=1,mul
2447 if(FT5(I).le.t2time) then
2448 write(93,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2449 endif
2450 1018 continue
2451clin-11/29/00 signal the end of animation file:
2452c write(91,*) -10.
2453c write(91,*) 0
2454c write(92,*) -10.
2455c write(92,*) 0
2456c close (91)
2457c close (92)
2458 close (93)
2459 endif
2460clin-5/18/01 calculate v2 for active partons:
2461 do 450 ianp=1,30
2462 do 400 iy=1,3
2463 v2pact=0d0
2464 v2p2ac=0d0
2465 s2pact=0d0
2466 etpact=0d0
2467 xnacti=0d0
2468 if(xnpf(ianp,iy).gt.1) then
2469c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2470 v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2471 v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2472 1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2473 s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2474 etpact=etpf(ianp,iy)*dble(nevt(ianp))
2475 xnpact=xnpf(ianp,iy)
2476c
2477 do 1019 kanp=1,ianp
2478 v2pact=v2pact-v2pfrz(kanp,iy)
2479 v2p2ac=v2p2ac-v2p2fz(kanp,iy)
2480 s2pact=s2pact-s2pfrz(kanp,iy)
2481 etpact=etpact-etpfrz(kanp,iy)
2482 xnpact=xnpact-xnpfrz(kanp,iy)
2483 1019 continue
2484c save the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2485 v2ph=v2pact
2486 v2ph2=v2p2ac
2487 s2ph=s2pact
2488 etph=etpact
2489 xnp2=xnpact/2d0
2490c
2491 if(xnpact.gt.1.and.nevt(ianp).ne.0) then
2492 v2pact=v2pact/xnpact
2493 v2p2ac=dsqrt((v2p2ac/xnpact
2494 1 -v2pact**2)/(xnpact-1))
2495 s2pact=s2pact/xnpact
2496 xnacti=dble(xnpact/dble(nevt(ianp)))
2497 etpact=etpact/dble(nevt(ianp))
2498cms write (nfile(iy)+2, 250) tsp(ianp),v2pact,
2499cms 1 v2p2ac,s2pact,etpact,xnacti,
2500cms 2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp)
2501 endif
2502 endif
2503c To calculate combined v2 for active partons plus formed hadrons,
2504c add the sum of v2h, v2h2, s2h, eth, and xnh for formed hadrons:
2505c scale the hadron part in case nevt(ianp) != nevent:
2506 shadr=dble(nevt(ianp))/dble(nevent)
2507 ianh=ianp
2508 v2ph=v2ph+v2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2509 v2ph2=v2ph2+(v2h2(ianh,iy)**2*(xnhadr(ianh,iy)-1)
2510 1 +v2h(ianh,iy)**2)*xnhadr(ianh,iy)*shadr
2511 s2ph=s2ph+s2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2512 etph=etph+eth(ianh,iy)*dble(nevent)*shadr
2513 xnph=xnpact+xnhadr(ianh,iy)*shadr
2514 xnp2h=xnp2+xnhadr(ianh,iy)*shadr
2515 if(xnph.gt.1) then
2516 v2ph=v2ph/xnph
2517 v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1))
2518 s2ph=s2ph/xnph
2519 etph=etph/dble(nevt(ianp))
2520 xnp2=xnp2/dble(nevt(ianp))
2521 xnp2h=xnp2h/dble(nevent)
2522cms if(tsp(ianp).le.(ntmax*dt))
2523cms 1 write (nfile(iy)+3, 250) tsp(ianp),v2ph,
2524cms 2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp)
2525 endif
2526c
2527 400 continue
2528 450 continue
2529 do 550 ianp=1,30
2530 do 500 iy=1,3
2531 v2pact=0d0
2532 v2p2ac=0d0
2533 s2pact=0d0
2534 etpact=0d0
2535 xnacti=0d0
2536c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2537 v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2538 v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2539 1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2540 s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2541 etpact=etpf(ianp,iy)*dble(nevt(ianp))
2542 xnpact=xnpf(ianp,iy)
2543 500 continue
2544 550 continue
2545cms close (620)
2546cms close (630)
2547 do 1021 nf=1,3
2548 do 1020 ifile=0,3
2549cms close(nfile(nf)+ifile)
2550 1020 continue
2551 1021 continue
2552 do 1022 nf=1,3
2553cms close(740+nf)
2554 1022 continue
2555 endif
2556cyy 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2)
2557cyy 210 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2))
2558 240 format(a30,3(2x,f9.5))
2559cyy 250 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6)
2560csp
2561cyy 201 format(2x,f5.2,4(2x,f9.2))
2562cyy 251 format(5e15.5)
2563c
2564 return
2565 end
2566
2567c=======================================================================
2568c Calculate flow from formed hadrons, called by art1e.f:
2569c Note: numbers in art not in double precision!
2570 subroutine flowh(ct)
2571 PARAMETER (MAXSTR=150001, MAXR=1)
2572 dimension tsh(31)
2573 DOUBLE PRECISION v2h,xnhadr,eth,v2h2,s2h
2574 DOUBLE PRECISION v2hp,xnhadp,v2hsum,v2h2sm,v2hevt(3)
2575 DOUBLE PRECISION pt2, v2hadr
2576 COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2577 1 v2h2(30,3),s2h(30,3)
2578cc SAVE /hflow/
2579 common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2580cc SAVE /ebe/
2581 common /lastt/itimeh,bimp
2582cc SAVE /lastt/
2583 COMMON /RUN/ NUM
2584cc SAVE /RUN/
2585 COMMON /AA/ R(3,MAXSTR)
2586cc SAVE /AA/
2587 COMMON /BB/ P(3,MAXSTR)
2588cc SAVE /BB/
2589 COMMON /CC/ E(MAXSTR)
2590cc SAVE /CC/
2591 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2592cc SAVE /EE/
2593 COMMON /RR/ MASSR(0:MAXR)
2594cc SAVE /RR/
2595 common/anim/nevent,isoft,isflag,izpc
2596cc SAVE /anim/
2597 COMMON /AREVT/ IAEVT, IARUN, MISS
2598cc SAVE /AREVT/
2599 SAVE
2600c
2601 do 1001 ii = 1, 31
2602 tsh(ii)=float(ii-1)
2603 1001 continue
2604c
2605 do 1004 ianh = 1, 30
2606 if ((ct+0.0001).lt.tsh(ianh+1)
2607 1 .and.(ct+0.0001).ge.tsh(ianh)) then
2608 if(ianh.eq.itimeh) goto 101
2609 IA = 0
2610 DO 1002 J = 1, NUM
2611 mult=MASSR(J)
2612 IA = IA + MASSR(J - 1)
2613 DO 100 IC = 1, mult
2614 I = IA + IC
2615c 5/04/01 exclude leptons and photons:
2616 if(iabs(LB(I)-10000).lt.100) goto 100
2617 px=p(1,i)
2618 py=p(2,i)
2619 pt2=dble(px)**2+dble(py)**2
2620c 2/18/00 Note: e(i) gives the mass in ART:
2621 ene=sqrt(e(i)**2+sngl(pt2)+p(3,i)**2)
2622 RAP=0.5*alog((ene+p(3,i))/(ene-p(3,i)))
2623ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2624c rt_now=sqrt(r(1,i)**2+r(2,i)**2)
2625c gz_now=r(3,i)
2626c zmax=0.1*ct
2627c volume=3.1416*(2.**2)*2*zmax
2628c if(rt_now.gt.2.or.abs(gz_now).gt.zmax)
2629c 1 goto 100
2630 iloop=1
2631 if(abs(rap).le.1) then
2632 iloop=2
2633 if(abs(rap).le.0.5) then
2634 iloop=3
2635 endif
2636 endif
2637 do 50 iy=1,iloop
2638 if(pt2.gt.0d0) then
2639 v2hadr=(dble(px)**2-dble(py)**2)/pt2
2640 v2h(ianh,iy)=v2h(ianh,iy)+v2hadr
2641 v2h2(ianh,iy)=v2h2(ianh,iy)+v2hadr**2
2642 if(dabs(v2hadr).gt.1d0)
2643 1 write(1,*) 'v2hadr>1',v2hadr,px,py
2644 endif
2645 xperp2=r(1,I)**2+r(2,I)**2
2646 if(xperp2.gt.0.)
2647 1 s2h(ianh,iy)=s2h(ianh,iy)+dble((r(1,I)**2-r(2,I)**2)/xperp2)
2648 eth(ianh,iy)=eth(ianh,iy)+dble(SQRT(e(i)**2+sngl(pt2)))
2649ctest off density:
2650c eth(ianh,iy)=eth(ianh,iy)
2651c 1 +dble(SQRT(e(i)**2+sngl(pt2)+p(3,i)**2))/volume
2652 xnhadr(ianh,iy)=xnhadr(ianh,iy)+1d0
2653 50 continue
2654 100 continue
2655 1002 CONTINUE
2656 itimeh=ianh
2657clin-5/04/01 ebe v2 variables:
2658 if(ianh.eq.30) then
2659 do 1003 iy=1,3
2660 nhadrn=IDINT(xnhadr(ianh,iy)-xnhadp(iy))
2661 if(nhadrn.ne.0) then
2662 v2hevt(iy)=(v2h(ianh,iy)-v2hp(iy))/dble(nhadrn)
2663 v2hsum(iy)=v2hsum(iy)+v2hevt(iy)
2664 v2h2sm(iy)=v2h2sm(iy)+v2hevt(iy)**2
2665 v2hp(iy)=v2h(ianh,iy)
2666 xnhadp(iy)=xnhadr(ianh,iy)
2667 endif
2668 1003 continue
2669 write(88, 160) iaevt,v2hevt
2670 endif
2671 goto 101
2672 endif
2673 1004 continue
2674 160 format(i10,3(2x,f9.5))
2675clin-11/27/00 for animation:
2676 101 if(nevent.eq.1) then
2677 IA = 0
2678 do 1005 J = 1, NUM
2679 mult=MASSR(J)
2680 IA = IA + MASSR(J - 1)
2681c write(91,*) ct
2682c write(91,*) mult
2683 DO 150 IC = 1, mult
2684 I = IA + IC
2685c write(91,210) LB(i),r(1,i),r(2,i),r(3,i),
2686c 1 p(1,i),p(2,i),p(3,i),e(i)
2687 150 continue
2688 1005 continue
2689 return
2690 endif
2691c 210 format(i6,7(1x,f8.3))
2692 return
2693 end
2694
2695c=======================================================================
2696 subroutine flowh0(NEVNT,idd)
2697c
2698 dimension tsh(31)
2699 DOUBLE PRECISION v2h,xnhadr,eth,v2h2,s2h
2700 DOUBLE PRECISION v2hp,xnhadp,v2hsum,v2h2sm,
2701 1 v2havg(3),varv2h(3)
2702 COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2703 1 v2h2(30,3),s2h(30,3)
2704cc SAVE /hflow/
2705 common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2706cc SAVE /ebe/
2707 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
2708cc SAVE /input1/
2709 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
2710 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
2711cc SAVE /INPUT2/
2712 common /lastt/itimeh,bimp
2713cc SAVE /lastt/
2714 SAVE
2715
2716c idd=0: initialization for flow analysis, called by artdri.f::
2717 if(idd.eq.0) then
2718 itimeh=0
2719c
2720 do 1001 ii = 1, 31
2721 tsh(ii)=float(ii-1)
2722 1001 continue
2723c
2724 do 1003 ii=1,30
2725 do 1002 iy=1,3
2726 v2h(ii,iy)=0d0
2727 xnhadr(ii,iy)=0d0
2728 eth(ii,iy)=0d0
2729 v2h2(ii,iy)=0d0
2730 s2h(ii,iy)=0d0
2731 1002 continue
2732 1003 continue
2733 do 1004 iy=1,3
2734 v2hp(iy)=0d0
2735 xnhadp(iy)=0d0
2736 v2hsum(iy)=0d0
2737 v2h2sm(iy)=0d0
2738 if(iy.eq.1) then
2739 nunit=59
2740 elseif(iy.eq.2) then
2741 nunit=68
2742 else
2743 nunit=69
2744 endif
2745 write(nunit,*) ' tsh, v2h, v2h2, s2h, '//
2746 1 ' eth, xmulth'
2747 1004 continue
2748c idd=2: calculate average hadronic elliptic flow, called by artdri.f:
2749 else if(idd.eq.2) then
2750 do 100 ii=1, 30
2751 do 1005 iy=1,3
2752 if(xnhadr(ii,iy).eq.0) then
2753 xmulth=0.
2754 elseif(xnhadr(ii,iy).gt.1) then
2755 v2h(ii,iy)=v2h(ii,iy)/xnhadr(ii,iy)
2756 eth(ii,iy)=eth(ii,iy)/dble(NEVNT)
2757 v2h2(ii,iy)=dsqrt((v2h2(ii,iy)/xnhadr(ii,iy)
2758 1 -v2h(ii,iy)**2)/(xnhadr(ii,iy)-1))
2759 s2h(ii,iy)=s2h(ii,iy)/xnhadr(ii,iy)
2760 xmulth=sngl(xnhadr(ii,iy)/NEVNT)
2761 endif
2762 if(iy.eq.1) then
2763 nunit=59
2764 elseif(iy.eq.2) then
2765 nunit=68
2766 else
2767 nunit=69
2768 endif
2769 if(tsh(ii).le.(ntmax*dt))
2770 1 write (nunit,200) tsh(ii),v2h(ii,iy),
2771 2 v2h2(ii,iy),s2h(ii,iy),eth(ii,iy),xmulth
2772 1005 continue
2773 100 continue
2774c event-by-event average & variance of v2h:
2775 do 1006 iy=1,3
2776 v2havg(iy)=v2hsum(iy)/dble(NEVNT)
2777 varv2h(iy)=dsqrt(v2h2sm(iy)/dble(NEVNT)-v2havg(iy)**2)
2778 1006 continue
2779 write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): avg=', v2havg
2780 write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): var=', varv2h
2781 endif
2782 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2))
2783 240 format(a30,3(2x,f9.5))
2784 return
2785 end
2786
2787c=======================================================================
2788c 2/23/00 flow from all initial hadrons just before entering ARTMN:
2789 subroutine iniflw(NEVNT,idd)
2790 PARAMETER (MAXSTR=150001, MAXR=1)
2791 DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult,
2792 1 v2bi,s2bi,xbmult
2793 COMMON /RUN/ NUM
2794cc SAVE /RUN/
2795 COMMON /ARERC1/MULTI1(MAXR)
2796cc SAVE /ARERC1/
2797 COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2798 & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
2799 & FT1(MAXSTR, MAXR),
2800 & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2801 & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2802cc SAVE /ARPRC1/
2803 COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
2804cc SAVE /iflow/
2805 SAVE
2806c
2807 if(idd.eq.0) then
2808 v2i=0d0
2809 eti=0d0
2810 xmulti=0d0
2811 v2mi=0d0
2812 s2mi=0d0
2813 xmmult=0d0
2814 v2bi=0d0
2815 s2bi=0d0
2816 xbmult=0d0
2817 else if(idd.eq.1) then
2818 do 1002 J = 1, NUM
2819 do 1001 I = 1, MULTI1(J)
2820 ITYP = ITYP1(I, J)
2821c all hadrons:
2822 IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 100
2823 xmulti=xmulti+1.d0
2824 PX = PX1(I, J)
2825 PY = PY1(I, J)
2826 XM = XM1(I, J)
2827 pt2=px**2+py**2
2828 xh=gx1(I,J)
2829 yh=gy1(I,J)
2830 xt2=xh**2+yh**2
2831 if(pt2.gt.0) v2i=v2i+dble((px**2-py**2)/pt2)
2832 eti=eti+dble(SQRT(PX ** 2 + PY ** 2 + XM ** 2))
2833c baryons only:
2834 IF (ITYP .LT. -1000 .or. ITYP .GT. 1000) then
2835 xbmult=xbmult+1.d0
2836 if(pt2.gt.0) v2bi=v2bi+dble((px**2-py**2)/pt2)
2837 if(xt2.gt.0) s2bi=s2bi+dble((xh**2-yh**2)/xt2)
2838c mesons only:
2839 else
2840 xmmult=xmmult+1.d0
2841 if(pt2.gt.0) v2mi=v2mi+dble((px**2-py**2)/pt2)
2842 if(xt2.gt.0) s2mi=s2mi+dble((xh**2-yh**2)/xt2)
2843 endif
2844 100 continue
2845 1001 continue
2846 1002 continue
2847 else if(idd.eq.2) then
2848 if(xmulti.ne.0) v2i=v2i/xmulti
2849 eti=eti/dble(NEVNT)
2850 xmulti=xmulti/dble(NEVNT)
2851 if(xmmult.ne.0) then
2852 v2mi=v2mi/xmmult
2853 s2mi=s2mi/xmmult
2854 endif
2855 xmmult=xmmult/dble(NEVNT)
2856 if(xbmult.ne.0) then
2857 v2bi=v2bi/xbmult
2858 s2bi=s2bi/xbmult
2859 endif
2860 xbmult=xbmult/dble(NEVNT)
2861 endif
2862c
2863 return
2864 end
2865
2866c=======================================================================
2867c 2/25/00 dN/dt analysis for production (before ZPCMN)
2868c and freezeout (right after ZPCMN) for all partons.
2869 subroutine frztm(NEVNT,idd)
2870c
2871 implicit double precision (a-h, o-z)
2872 PARAMETER (MAXPTN=400001)
2873 dimension tsf(31)
2874 COMMON /PARA1/ MUL
2875cc SAVE /PARA1/
2876 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
2877 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
2878 & XMASS0(MAXPTN), ITYP0(MAXPTN)
2879cc SAVE /prec1/
2880 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
2881 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
2882 & XMASS5(MAXPTN), ITYP5(MAXPTN)
2883cc SAVE /prec2/
2884 COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
2885 & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
2886cc SAVE /frzout/
2887 SAVE
2888 data tsf/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2889 & 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
2890 & 2 , 3, 4, 5, 6, 7, 8, 9, 10, 20, 30/
2891c
2892 if(idd.eq.0) then
2893 do 1001 ii=1,30
2894 xnprod(ii)=0d0
2895 etprod(ii)=0d0
2896 xnfrz(ii)=0d0
2897 etfrz(ii)=0d0
2898 dnprod(ii)=0d0
2899 detpro(ii)=0d0
2900 dnfrz(ii)=0d0
2901 detfrz(ii)=0d0
2902 1001 continue
2903cms OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN')
2904cms OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN')
2905 else if(idd.eq.1) then
2906 DO 100 ip = 1, MUL
2907 do 1002 ii=1,30
2908 eth0=dSQRT(PX0(ip)**2+PY0(ip)**2+XMASS0(ip)**2)
2909 eth2=dSQRT(PX5(ip)**2+PY5(ip)**2+XMASS5(ip)**2)
2910c total number and Et produced by time tsf(ii):
2911 if (ft0(ip).lt.tsf(ii+1)) then
2912 xnprod(ii)=xnprod(ii)+1d0
2913 etprod(ii)=etprod(ii)+eth0
2914c number and Et produced from time tsf(ii) to tsf(ii+1):
2915 if (ft0(ip).ge.tsf(ii)) then
2916 dnprod(ii)=dnprod(ii)+1d0
2917 detpro(ii)=detpro(ii)+eth0
2918 endif
2919 endif
2920c total number and Et freezed out by time tsf(ii):
2921 if (FT5(ip).lt.tsf(ii+1)) then
2922 xnfrz(ii)=xnfrz(ii)+1d0
2923 etfrz(ii)=etfrz(ii)+eth2
2924c number and Et freezed out from time tsf(ii) to tsf(ii+1):
2925 if (FT5(ip).ge.tsf(ii)) then
2926 dnfrz(ii)=dnfrz(ii)+1d0
2927 detfrz(ii)=detfrz(ii)+eth2
2928 endif
2929 endif
2930 1002 continue
2931 100 continue
2932 else if(idd.eq.2) then
2933cms write (86,*) ' t, np, dnp/dt, etp '//
2934cms 1 ' detp/dt'
2935cms write (87,*) ' t, nf, dnf/dt, etf '//
2936cms 1 ' detf/dt'
2937 do 1003 ii=1,30
2938 xnp=xnprod(ii)/dble(NEVNT)
2939 xnf=xnfrz(ii)/dble(NEVNT)
2940 etp=etprod(ii)/dble(NEVNT)
2941 etf=etfrz(ii)/dble(NEVNT)
2942 dxnp=dnprod(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2943 dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2944 detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2945 detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2946cms write (86, 200)
2947cms 1 tsf(ii+1),xnp,dxnp,etp,detp
2948cms write (87, 200)
2949cms 1 tsf(ii+1),xnf,dxnf,etf,detf
2950 1003 continue
2951 endif
2952cyy 200 format(2x,f9.2,4(2x,f10.2))
2953c
2954 return
2955 end
2956
2957c=======================================================================
2958clin-6/2009 write out initial minijet information
2959c before propagating to its formation time:
2960 subroutine minijet_out(BB)
2961 PARAMETER (MAXSTR=150001)
2962 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2963 COMMON/hjcrdn/YP(3,300),YT(3,300)
2964 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2965 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2966 & PJPM(300,500),NTJ(300),KFTJ(300,500),
2967 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2968 & PJTE(300,500),PJTM(300,500)
2969 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2970 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2971 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2972 COMMON /AREVT/ IAEVT, IARUN, MISS
2973 common /para7/ ioscar,nsmbbbar,nsmmeson
2974 common/phidcy/iphidcy,pttrig,ntrig,maxmiss
2975 SAVE
2976 ntrig=0
2977 do I = 1, IHNT2(1)
2978 do J = 1, NPJ(I)
2979 pt=sqrt(PJPX(I,J)**2+PJPY(I,J)**2)
2980 if(pt.ge.pttrig) ntrig=ntrig+1
2981 enddo
2982 enddo
2983 do I = 1, IHNT2(3)
2984 do J = 1, NTJ(I)
2985 pt=sqrt(PJTX(I,J)**2+PJTY(I,J)**2)
2986 if(pt.ge.pttrig) ntrig=ntrig+1
2987 enddo
2988 enddo
2989 do I = 1, NSG
2990 do J = 1, NJSG(I)
2991 pt=sqrt(PXSG(I,J)**2+PYSG(I,J)**2)
2992 if(pt.ge.pttrig) ntrig=ntrig+1
2993 enddo
2994 enddo
2995c Require at least 1 initial minijet parton above the trigger Pt value:
2996 if(ntrig.eq.0) return
2997
2998c.....transfer data from HIJING to ZPC
2999 if(ioscar.eq.3) write(96,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
3000 DO 1008 I = 1, IHNT2(1)
3001 DO 1007 J = 1, NPJ(I)
3002 ityp=KFPJ(I,J)
3003c write out not only gluons:
3004c if(ityp.ne.21) goto 1007
3005 gx=YP(1,I)+0.5*BB
3006 gy=YP(2,I)
3007 gz=0d0
3008 ft=0d0
3009 px=PJPX(I,J)
3010 py=PJPY(I,J)
3011 pz=PJPZ(I,J)
3012 xmass=PJPM(I,J)
3013 if(ioscar.eq.3) then
3014 if(amax1(abs(gx),abs(gy),
3015 1 abs(gz),abs(ft)).lt.9999) then
3016 write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3017 else
3018 write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3019 endif
3020 endif
3021 1007 CONTINUE
3022 1008 CONTINUE
3023 DO 1010 I = 1, IHNT2(3)
3024 DO 1009 J = 1, NTJ(I)
3025 ityp=KFTJ(I,J)
3026c if(ityp.ne.21) goto 1009
3027 gx=YT(1,I)-0.5*BB
3028 gy=YT(2,I)
3029 gz=0d0
3030 ft=0d0
3031 px=PJTX(I,J)
3032 py=PJTY(I,J)
3033 pz=PJTZ(I,J)
3034 xmass=PJTM(I,J)
3035 if(ioscar.eq.3) then
3036 if(amax1(abs(gx),abs(gy),
3037 1 abs(gz),abs(ft)).lt.9999) then
3038 write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3039 else
3040 write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3041 endif
3042 endif
3043 1009 CONTINUE
3044 1010 CONTINUE
3045 DO 1012 I = 1, NSG
3046 DO 1011 J = 1, NJSG(I)
3047 ityp=K2SG(I,J)
3048c if(ityp.ne.21) goto 1011
3049 gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2)))
3050 gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2)))
3051 gz=0d0
3052 ft=0d0
3053 px=PXSG(I,J)
3054 py=PYSG(I,J)
3055 pz=PZSG(I,J)
3056 xmass=PMSG(I,J)
3057 if(ioscar.eq.3) then
3058 if(amax1(abs(gx),abs(gy),
3059 1 abs(gz),abs(ft)).lt.9999) then
3060 write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3061 else
3062 write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3063 endif
3064 endif
3065 1011 CONTINUE
3066 1012 CONTINUE
3067 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,f8.2),2(2x,f2.0),2x,I2)
3068 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,e8.2),2(2x,f2.0),2x,I2)
3069c
3070 return
3071 end
3072
3073c=======================================================================
3074clin-6/2009 embed back-to-back high-Pt quark/antiquark pair
3075c via embedding back-to-back high-Pt pion pair then melting the pion pair
3076c by generating the internal quark and antiquark momentum parallel to
3077c the pion momentum (in order to produce a high-Pt and a low Pt parton):
3078 subroutine embedHighPt
3079 PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140)
3080 common/embed/iembed,pxqembd,pyqembd,xembd,yembd
3081 COMMON/RNDF77/NSEED
3082 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
3083 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
3084 COMMON /ARPRC/ ITYPAR(MAXSTR),
3085 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
3086 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
3087 & XMAR(MAXSTR)
3088cwei DOUBLE PRECISION PATT
3089 SAVE
3090c
3091 if(iembed.ne.1) return
3092 ptq=sqrt(pxqembd**2+pyqembd**2)
3093 if(ptq.lt.(pichmass/2.)) then
3094 print *, 'Embedded quark transverse momentum is too small'
3095 stop
3096 endif
3097c Randomly embed u/ubar or d/dbar at high Pt:
3098 idqembd=1+int(2*RANART(NSEED))
3099c Flavor content for the charged pion that contains the leading quark:
3100 if(idqembd.eq.1) then
3101 idqsoft=-2
3102 idpi=-211
3103 elseif(idqembd.eq.2) then
3104 idqsoft=-1
3105 idpi=211
3106 else
3107 print *, 'Wrong quark flavor embedded'
3108 stop
3109 endif
3110c Caculate transverse momentum of the parent charged pion:
3111 xmq=ulmass(idqembd)
3112 xmqsoft=ulmass(idqsoft)
3113 ptpi=((pichmass**2+xmq**2-xmqsoft**2)*ptq
3114 1 -sqrt((xmq**2+ptq**2)*(pichmass**4
3115 2 -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2)))
3116 3 /(2.*xmq**2)
3117 pxpi=ptpi*pxqembd/ptq
3118 pypi=ptpi*pyqembd/ptq
3119c Embedded quark/antiquark are assumed to have pz=0:
3120 pzpi=0.
3121c Insert the two parent charged pions,
3122c ipion=1 for the pion containing the leading quark,
3123c ipion=2 for the pion containing the leading antiquark of the same flavor:
3124 do ipion=1,2
3125 if(ipion.eq.2) then
3126 idpi=-idpi
3127 pxpi=-pxpi
3128 pypi=-pypi
3129 pzpi=-pzpi
3130 endif
3131 NATT=NATT+1
3132 KATT(NATT,1)=idpi
3133 KATT(NATT,2)=40
3134 KATT(NATT,3)=0
3135 PATT(NATT,1)=pxpi
3136 PATT(NATT,2)=pypi
3137 PATT(NATT,3)=pzpi
3138 PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2)
3139 EATT=EATT+PATT(NATT,4)
3140 GXAR(NATT)=xembd
3141 GYAR(NATT)=yembd
3142 GZAR(NATT)=0.
3143 FTAR(NATT)=0.
3144 ITYPAR(NATT)=KATT(NATT,1)
3145 PXAR(NATT)=PATT(NATT,1)
3146 PYAR(NATT)=PATT(NATT,2)
3147 PZAR(NATT)=PATT(NATT,3)
3148 PEAR(NATT)=PATT(NATT,4)
3149 XMAR(NATT)=pichmass
3150 enddo
3151c
3152 return
3153 end