]> 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
- Adding sparse histograms to analysis
[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
1617 & .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1618 & .or.iabs(lb1).eq.16) then
1619 kf=INVFLV(lb1)
1620 else
1621 return
1622 endif
1623c
1624 IP=1
1625c label as undecayed and the only particle in the record:
1626 N=1
1627 K(IP,1)=1
1628 K(IP,3)=0
1629 K(IP,4)=0
1630 K(IP,5)=0
1631c
1632 K(IP,2)=kf
1633 P(IP,1)=px1
1634 P(IP,2)=py1
1635 P(IP,3)=pz1
1636 em1a=em1
1637c eta or omega in ART may be below or too close to (pi+pi-pi0) mass,
1638c causing LUDECY error,thus increase their mass ADDM above this thresh,
1639c noting that rho (m=0.281) too close to 2pi thrshold fails to decay:
1640 if((lb1.eq.0.or.lb1.eq.28).and.em1.lt.(2*APICH+API0+ADDM)) then
1641 em1=2*APICH+API0+ADDM
1642c rho
1643 elseif(lb1.ge.25.and.lb1.le.27.and.em1.lt.(2*APICH+ADDM)) then
1644 em1=2*APICH+ADDM
1645c K*
1646 elseif(iabs(lb1).eq.30.and.em1.lt.(APICH+AK0+ADDM)) then
1647 em1=APICH+AK0+ADDM
1648c Delta created in ART may be below (n+pich) mass, causing LUDECY error:
1649 elseif(iabs(lb1).ge.6.and.iabs(lb1).le.9
1650 1 .and.em1.lt.(APICH+AN+ADDM)) then
1651 em1=APICH+AN+ADDM
1652 endif
1653c if(em1.ge.(em1a+0.01)) write (6,*)
1654c 1 'Mass increase in resdec():',nt,em1-em1a,lb1
1655 e1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1656 P(IP,4)=e1
1657 P(IP,5)=em1
1658clin-5/2008:
1659 dpdecp=dpertp(i1)
1660 call ludecy(IP)
1661c add decay time to daughter's formation time at the last timestep:
1662 if(nt.eq.ntmax) then
1663 tau0=hbarc/wid
1664 taudcy=tau0*(-1.)*alog(1.-RANART(NSEED))
1665 ndaut=n-nsav
1666 if(ndaut.le.1) then
1667 write(10,*) 'note: ndaut(<1)=',ndaut
1668 write(89,*) 'note: ndaut(<1)=',ndaut
1669 call lulist(2)
1670 stop
1671 endif
1672c lorentz boost:
1673 taudcy=taudcy*e1/em1
1674 tfnl=tfnl+taudcy
1675 xfnl=xfnl+px1/e1*taudcy
1676 yfnl=yfnl+py1/e1*taudcy
1677 zfnl=zfnl+pz1/e1*taudcy
1678c at the last timestep, assign rho, K0S or eta (decay daughter)
1679c to lb(i1) only (not to lpion) in order to decay them again:
1680 if(n.ge.(nsav+2)) then
1681 do 1001 idau=nsav+2,n
1682 kdaut=K(idau,2)
1683 if(kdaut.eq.221.or.kdaut.eq.113
1684 1 .or.kdaut.eq.213.or.kdaut.eq.-213
1685 2 .or.kdaut.eq.310) then
1686c switch idau and i1(nsav+1):
1687 ksave=kdaut
1688 pxsave=p(idau,1)
1689 pysave=p(idau,2)
1690 pzsave=p(idau,3)
1691 esave=p(idau,4)
1692 xmsave=p(idau,5)
1693 K(idau,2)=K(nsav+1,2)
1694 p(idau,1)=p(nsav+1,1)
1695 p(idau,2)=p(nsav+1,2)
1696 p(idau,3)=p(nsav+1,3)
1697 p(idau,4)=p(nsav+1,4)
1698 p(idau,5)=p(nsav+1,5)
1699 K(nsav+1,2)=ksave
1700 p(nsav+1,1)=pxsave
1701 p(nsav+1,2)=pysave
1702 p(nsav+1,3)=pzsave
1703 p(nsav+1,4)=esave
1704 p(nsav+1,5)=xmsave
1705c note: phi decay may produce rho, K0s or eta, N*(1535) decay may produce
1706c eta, but only one daughter may be rho, K0s or eta:
1707 goto 111
1708 endif
1709 1001 continue
1710 endif
1711 111 continue
1712c
1713 enet=0.
1714 do 1002 idau=nsav+1,n
1715 enet=enet+p(idau,4)
1716 1002 continue
1717c if(abs(enet-e1).gt.0.02)
1718c 1 write(93,*) 'resdec(): nt=',nt,enet-e1,lb1
1719 endif
1720
1721cyy 200 format(a20,3(1x,i6))
1722cyy 210 format(i6,5(1x,f8.3))
1723cyy 220 format(a2,i5,5(1x,f8.3))
1724
1725 do 1003 idau=nsav+1,n
1726 kdaut=K(idau,2)
1727 lbdaut=IARFLV(kdaut)
1728c K0S and K0L are named K+/K- during hadron cascade, and only
1729c at the last timestep they keep their real LB # before output;
1730c K0/K0bar (from K* decay) converted to K0S and K0L at the last timestep:
1731 if(nt.eq.ntmax.and.(kdaut.eq.130.or.kdaut.eq.310
1732 1 .or.iabs(kdaut).eq.311)) then
1733 if(kdaut.eq.130) then
1734 lbdaut=22
1735 elseif(kdaut.eq.310) then
1736 lbdaut=24
1737 elseif(iabs(kdaut).eq.311) then
1738 if(RANART(NSEED).lt.0.5) then
1739 lbdaut=22
1740 else
1741 lbdaut=24
1742 endif
1743 endif
1744 endif
1745c
1746 if(idau.eq.(nsav+1)) then
1747 LB(i1)=lbdaut
1748 E(i1)=p(idau,5)
1749 px1n=p(idau,1)
1750 py1n=p(idau,2)
1751 pz1n=p(idau,3)
1752clin-5/2008:
1753 dp1n=dpdecp
1754 else
1755 nnn=nnn+1
1756 LPION(NNN,IRUN)=lbdaut
1757 EPION(NNN,IRUN)=p(idau,5)
1758 PPION(1,NNN,IRUN)=p(idau,1)
1759 PPION(2,NNN,IRUN)=p(idau,2)
1760 PPION(3,NNN,IRUN)=p(idau,3)
1761 RPION(1,NNN,IRUN)=xfnl
1762 RPION(2,NNN,IRUN)=yfnl
1763 RPION(3,NNN,IRUN)=zfnl
1764 tfdpi(NNN,IRUN)=tfnl
1765clin-5/2008:
1766 dppion(NNN,IRUN)=dpdecp
1767 endif
1768 1003 continue
1769cyy 230 format(a2,i5,5(1x,e8.2))
1770 return
1771 end
1772
1773c=======================================================================
1774 subroutine inidcy
1775
1776 COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1777cc SAVE /LUJETSA/
1778 common/resdcy/NSAV,iksdcy
1779cc SAVE /resdcy/
1780 SAVE
1781 N=1
1782 NSAV=N
1783 return
1784 end
1785
1786c=======================================================================
1787clin-6/06/02 local parton freezeout motivated from critical density:
1788 subroutine local(t)
1789c
1790 implicit double precision (a-h, o-z)
1791 PARAMETER (MAXPTN=400001)
1792 PARAMETER (r0=1d0)
1793 COMMON /para1/ mul
1794cc SAVE /para1/
1795 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1796 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1797 & XMASS5(MAXPTN), ITYP5(MAXPTN)
1798cc SAVE /prec2/
1799 common /frzprc/
1800 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1801 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1802 & xmfrz(MAXPTN),
1803 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1804cc SAVE /frzprc/
1805 common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1806cc SAVE /prec4/
1807 common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1808cc SAVE /prec5/
1809 common /coal/dpcoal,drcoal,ecritl
1810cc SAVE /coal/
1811 SAVE
1812c
1813 do 1001 it=1,301
1814 if(t.ge.tfrz(it).and.t.lt.tfrz(it+1)) then
1815 if(it.eq.itlast) then
1816 return
1817 else
1818 itlast=it
1819 goto 50
1820 endif
1821 endif
1822 1001 continue
1823 write(1,*) 'local time out of range in LOCAL, stop',t,it
1824 stop
1825 50 continue
1826c
1827 do 200 ip=1,mul
1828c skip partons which have frozen out:
1829 if(ifrz(ip).eq.1) goto 200
1830 if(it.eq.301) then
1831c freezeout all the left partons beyond the time of 3000 fm:
1832 etcrit=1d6
1833 goto 150
1834 else
1835c freezeout when transverse energy density < etcrit:
1836 etcrit=(ecritl*2d0/3d0)
1837 endif
1838c skip partons which have not yet formed:
1839 if(t.lt.FT5(ip)) goto 200
1840 rap0=rap(ip)
1841 eta0=eta(ip)
1842 x0=GX5(ip)+vx(ip)*(t-FT5(ip))
1843 y0=GY5(ip)+vy(ip)*(t-FT5(ip))
1844 detdy=0d0
1845 do 100 itest=1,mul
1846c skip self and partons which have not yet formed:
1847 if(itest.eq.ip.or.t.lt.FT5(itest)) goto 100
1848 ettest=eta(itest)
1849 xtest=GX5(itest)+vx(itest)*(t-FT5(itest))
1850 ytest=GY5(itest)+vy(itest)*(t-FT5(itest))
1851 drt=sqrt((xtest-x0)**2+(ytest-y0)**2)
1852c count partons within drt<1 and -1<(eta-eta0)<1:
1853 if(dabs(ettest-eta0).le.1d0.and.drt.le.r0)
1854 1 detdy=detdy+dsqrt(PX5(itest)**2+PY5(itest)**2
1855 2 +XMASS5(itest)**2)*0.5d0
1856 100 continue
1857 detdy=detdy*(dcosh(eta0)**2)/(t*3.1416d0*r0**2*dcosh(rap0))
1858c when density is below critical density for phase transition, freeze out:
1859 150 if(detdy.le.etcrit) then
1860 ifrz(ip)=1
1861 idfrz(ip)=ITYP5(ip)
1862 pxfrz(ip)=PX5(ip)
1863 pyfrz(ip)=PY5(ip)
1864 pzfrz(ip)=PZ5(ip)
1865 efrz(ip)=E5(ip)
1866 xmfrz(ip)=XMASS5(ip)
1867 if(t.gt.FT5(ip)) then
1868 gxfrz(ip)=x0
1869 gyfrz(ip)=y0
1870 gzfrz(ip)=GZ5(ip)+vz(ip)*(t-FT5(ip))
1871 ftfrz(ip)=t
1872 else
1873c if this freezeout time < formation time, use formation time & positions.
1874c This ensures the recovery of default hadron when e_crit=infty:
1875 gxfrz(ip)=GX5(ip)
1876 gyfrz(ip)=GY5(ip)
1877 gzfrz(ip)=GZ5(ip)
1878 ftfrz(ip)=FT5(ip)
1879 endif
1880 endif
1881 200 continue
1882c
1883 return
1884 end
1885
1886c=======================================================================
1887clin-6/06/02 initialization for local parton freezeout
1888 subroutine inifrz
1889c
1890 implicit double precision (a-h, o-z)
1891 PARAMETER (MAXPTN=400001)
1892 common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1893cc SAVE /ilist5/
1894 common /frzprc/
1895 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1896 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1897 & xmfrz(MAXPTN),
1898 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1899cc SAVE /frzprc/
1900 SAVE
1901c
1902c for freezeout time 0-10fm, use interval of 0.1fm;
1903c for 10-100fm, use interval of 1fm;
1904c for 100-1000fm, use interval of 10fm;
1905c for 1000-3000fm, use interval of 100fm:
1906 step1=0.1d0
1907 step2=1d0
1908 step3=10d0
1909 step4=100d0
1910c
1911 do 1001 it=1,101
1912 tfrz(it)=0d0+dble(it-1)*step1
1913 1001 continue
1914 do 1002 it=102,191
1915 tfrz(it)=10d0+dble(it-101)*step2
1916 1002 continue
1917 do 1003 it=192,281
1918 tfrz(it)=100d0+dble(it-191)*step3
1919 1003 continue
1920 do 1004 it=282,301
1921 tfrz(it)=1000d0+dble(it-281)*step4
1922 1004 continue
1923 tfrz(302)=tlarge
1924c
1925 return
1926 end
1927
1928clin-5/2009 ctest on v2 analysis
1929c=======================================================================
1930c idd=0,1,2,3 specifies different subroutines for partonic flow analysis.
1931 subroutine flowp(idd)
1932c
1933 implicit double precision (a-h, o-z)
1934 real dt
1935 parameter (MAXPTN=400001)
1936csp
1937 parameter (bmt=0.05d0)
1938 dimension nlfile(3),nsfile(3),nmfile(3)
1939c
1940 dimension v2pp(3),xnpp(3),v2psum(3),v2p2sm(3),nfile(3)
1941 dimension tsp(31),v2pevt(3),v2pavg(3),varv2p(3)
1942 common /ilist1/
1943 & iscat, jscat, next(MAXPTN), last(MAXPTN),
1944 & ictype, icsta(MAXPTN),
1945 & nic(MAXPTN), icels(MAXPTN)
1946cc SAVE /ilist1/
1947 COMMON /para1/ mul
1948cc SAVE /para1/
1949 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1950 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1951 & XMASS5(MAXPTN), ITYP5(MAXPTN)
1952cc SAVE /prec2/
1953 COMMON /pflow/ v2p(30,3),xnpart(30,3),etp(30,3),
1954 1 s2p(30,3),v2p2(30,3),nevt(30)
1955cc SAVE /pflow/
1956 COMMON /pflowf/ v2pf(30,3),xnpf(30,3),etpf(30,3),
1957 1 xncoll(30),s2pf(30,3),v2pf2(30,3)
1958cc SAVE /pflowf/
1959 COMMON /pfrz/ v2pfrz(30,3),xnpfrz(30,3),etpfrz(30,3),
1960 1 s2pfrz(30,3),v2p2fz(30,3),tscatt(31),
1961 2 nevtfz(30),iscatt(30)
1962cc SAVE /pfrz/
1963 COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
1964 1 v2h2(30,3),s2h(30,3)
1965cc SAVE /hflow/
1966 COMMON /AREVT/ IAEVT, IARUN, MISS
1967cc SAVE /AREVT/
1968 common/anim/nevent,isoft,isflag,izpc
1969cc SAVE /anim/
1970 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1971cc SAVE /input1/
1972 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1973 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1974cc SAVE /INPUT2/
1975cc SAVE itimep,iaevtp,v2pp,xnpp,v2psum,v2p2sm
1976cc SAVE nfile,itanim,nlfile,nsfile,nmfile
1977 SAVE
1978csp
1979 dimension etpl(30,3),etps(30,3),etplf(30,3),etpsf(30,3),
1980 & etlfrz(30,3),etsfrz(30,3),
1981 & xnpl(30,3),xnps(30,3),xnplf(30,3),xnpsf(30,3),
1982 & xnlfrz(30,3),xnsfrz(30,3),
1983 & v2pl(30,3),v2ps(30,3),v2plf(30,3),v2psf(30,3),
1984 & s2pl(30,3),s2ps(30,3),s2plf(30,3),s2psf(30,3),
1985 & DMYil(50,3),DMYfl(50,3),
1986 & DMYis(50,3),DMYfs(50,3)
1987 data tsp/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
1988 & 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
1989 & 2 , 3, 4, 5, 6, 7, 8, 9, 10, 20, 30/
1990c idd=0: initialization for flow analysis, called by artdri.f:
1991 if(idd.eq.0) then
1992 nfile(1)=60
1993 nfile(2)=64
1994 nfile(3)=20
1995cms OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN')
1996cms OPEN (nfile(1)+1,
1997cms 1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN')
1998cms OPEN (nfile(1)+2,
1999cms 1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN')
2000cms OPEN (nfile(1)+3,
2001cms 1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN')
2002cms OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN')
2003cms OPEN (nfile(2)+1,
2004cms 1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN')
2005cms OPEN (nfile(2)+2,
2006cms 1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN')
2007cms OPEN (nfile(2)+3,
2008cms 1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN')
2009cms OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN')
2010cms OPEN (nfile(3)+1,
2011cms 1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN')
2012cms OPEN (nfile(3)+2,
2013cms 1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN')
2014cms OPEN (nfile(3)+3,
2015cms 1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN')
2016cms OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN')
2017cms write(49, *) ' ievt, v2p, v2p_y2, v2p_y1'
2018c
2019cms OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN')
2020cms OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN')
2021cms OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN')
2022cms OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN')
2023cms write(88, *) ' ievt, v2h, v2h_y2, v2h_y1'
2024csp07/05
2025 nlfile(1)=70
2026 nlfile(2)=72
2027 nlfile(3)=74
2028cms OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN')
2029cms OPEN (nlfile(1)+1,
2030cms 1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN')
2031cms OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN')
2032cms OPEN (nlfile(2)+1,
2033cms 1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN')
2034cms OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN')
2035cms OPEN (nlfile(3)+1,
2036cms 1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN')
2037 nsfile(1)=76
2038 nsfile(2)=78
2039 nsfile(3)=80
2040cms OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN')
2041cms OPEN (nsfile(1)+1,
2042cms 1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN')
2043cms OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN')
2044cms OPEN (nsfile(2)+1,
2045cms 1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN')
2046cms OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN')
2047cms OPEN (nsfile(3)+1,
2048cms 1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN')
2049 nmfile(1)=82
2050 nmfile(2)=83
2051 nmfile(3)=84
2052cms OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN')
2053cms OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN')
2054cms OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN')
2055clin-11/27/00 for animation:
2056 if(nevent.eq.1) then
2057cms OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN')
2058cms write(91,*) ntmax, dt
2059cms OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN')
2060cms OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN')
2061 endif
2062c
2063 itimep=0
2064 itanim=0
2065 iaevtp=0
2066csp
2067 do 1002 ii=1,50
2068 do 1001 iy=1,3
2069 DMYil(ii,iy) = 0d0
2070 DMYfl(ii,iy) = 0d0
2071 DMYis(ii,iy) = 0d0
2072 DMYfs(ii,iy) = 0d0
2073 1001 continue
2074 1002 continue
2075c
2076 do 1003 ii=1,31
2077 tscatt(ii)=0d0
2078 1003 continue
2079 do 1005 ii=1,30
2080 nevt(ii)=0
2081 xncoll(ii)=0d0
2082 nevtfz(ii)=0
2083 iscatt(ii)=0
2084 do 1004 iy=1,3
2085 v2p(ii,iy)=0d0
2086 v2p2(ii,iy)=0d0
2087 s2p(ii,iy)=0d0
2088 etp(ii,iy)=0d0
2089 xnpart(ii,iy)=0d0
2090 v2pf(ii,iy)=0d0
2091 v2pf2(ii,iy)=0d0
2092 s2pf(ii,iy)=0d0
2093 etpf(ii,iy)=0d0
2094 xnpf(ii,iy)=0d0
2095 v2pfrz(ii,iy)=0d0
2096 v2p2fz(ii,iy)=0d0
2097 s2pfrz(ii,iy)=0d0
2098 etpfrz(ii,iy)=0d0
2099 xnpfrz(ii,iy)=0d0
2100csp07/05
2101 etpl(ii,iy)=0d0
2102 etps(ii,iy)=0d0
2103 etplf(ii,iy)=0d0
2104 etpsf(ii,iy)=0d0
2105 etlfrz(ii,iy)=0d0
2106 etsfrz(ii,iy)=0d0
2107 xnpl(ii,iy)=0d0
2108 xnps(ii,iy)=0d0
2109 xnplf(ii,iy)=0d0
2110 xnpsf(ii,iy)=0d0
2111 xnlfrz(ii,iy)=0d0
2112 xnsfrz(ii,iy)=0d0
2113 v2pl(ii,iy)=0d0
2114 v2ps(ii,iy)=0d0
2115 v2plf(ii,iy)=0d0
2116 v2psf(ii,iy)=0d0
2117 s2pl(ii,iy)=0d0
2118 s2ps(ii,iy)=0d0
2119 s2plf(ii,iy)=0d0
2120 s2psf(ii,iy)=0d0
2121 1004 continue
2122 1005 continue
2123 do 1006 iy=1,3
2124 v2pevt(iy)=0d0
2125 v2pavg(iy)=0d0
2126 varv2p(iy)=0d0
2127 v2pp(iy)=0.d0
2128 xnpp(iy)=0d0
2129 v2psum(iy)=0.d0
2130 v2p2sm(iy)=0.d0
2131 1006 continue
2132c idd=1: calculate parton elliptic flow, called by zpc.f:
2133 else if(idd.eq.1) then
2134 t2time = FT5(iscat)
2135 do 1008 ianp = 1, 30
2136 if (t2time.lt.tsp(ianp+1).and.t2time.ge.tsp(ianp)) then
2137c write flow info only once at each fixed time:
2138 xncoll(ianp)=xncoll(ianp)+1d0
2139c to prevent an earlier t2time comes later in the same event
2140c and mess up nevt:
2141 if(ianp.le.itimep.and.iaevt.eq.iaevtp) goto 101
2142 nevt(ianp)=nevt(ianp)+1
2143 tscatt(ianp+1)=t2time
2144 iscatt(ianp)=1
2145 nevtfz(ianp)=nevtfz(ianp)+1
2146 do 100 I=1,mul
2147 ypartn=0.5d0*dlog((E5(i)+PZ5(i))/(E5(i)-PZ5(i)+1.d-8))
2148 pt2=PX5(I)**2+PY5(I)**2
2149ctest off: initial (pt,y) and (x,y) distribution:
2150c idtime=1
2151c if(ianp.eq.idtime) then
2152c iityp=iabs(ITYP5(I))
2153c if(iityp.eq.1.or.iityp.eq.2) then
2154c write(651,*) dsqrt(pt2),ypartn
2155c write(654,*) GX5(I),GY5(I)
2156c elseif(iityp.eq.1103.or.iityp.eq.2101
2157c 1 .or.iityp.eq.2103.or.iityp.eq.2203.
2158c 2 .or.iityp.eq.3101.or.iityp.eq.3103.
2159c 3 .or.iityp.eq.3201.or.iityp.eq.3203.or.iityp.eq.3303)
2160c 4 then
2161c write(652,*) dsqrt(pt2),ypartn
2162c write(655,*) GX5(I),GY5(I)
2163c elseif(iityp.eq.21) then
2164c write(653,*) dsqrt(pt2),ypartn
2165c write(656,*) GX5(I),GY5(I)
2166c endif
2167c endif
2168ctest-end
2169ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2170c gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2171c gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2172c gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2173c rt_now=dsqrt(gx_now**2+gy_now**2)
2174c zmax=0.1d0*t2time
2175c volume=3.1416d0*(2d0**2)*(2*zmax)
2176c if(rt_now.gt.2d0.or.dabs(gz_now).gt.zmax)
2177c 1 goto 100
2178ctest-end
2179 iloop=1
2180 if(dabs(ypartn).le.1d0) then
2181 iloop=2
2182 if(dabs(ypartn).le.0.5d0) then
2183 iloop=3
2184 endif
2185 endif
2186 do 50 iy=1,iloop
2187 if(pt2.gt.0.) then
2188 v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2189 if(abs(v2prtn).gt.1.)
2190 1 write(nfile(iy),*) 'v2prtn>1',v2prtn
2191 v2p(ianp,iy)=v2p(ianp,iy)+v2prtn
2192 v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2
2193 endif
2194 xperp2=GX5(I)**2+GY5(I)**2
2195 if(xperp2.gt.0.)
2196 1 s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2197 xnpart(ianp,iy)=xnpart(ianp,iy)+1d0
2198 etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2199ctest off density:
2200c etp(ianp,iy)=etp(ianp,iy)
2201c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2202clin-2/22/00 to write out parton info only for formed ones:
2203 if(FT5(I).le.t2time) then
2204 v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn
2205 v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2
2206 if(xperp2.gt.0.)
2207 1 s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2208 xnpf(ianp,iy)=xnpf(ianp,iy)+1d0
2209 etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2210ctest off density:
2211c etpf(ianp,iy)=etpf(ianp,iy)
2212c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2213 endif
2214 50 continue
2215 100 continue
2216 itimep=ianp
2217 iaevtp=iaevt
2218clin-3/30/00 ebe v2 variables:
2219 if(ianp.eq.30) then
2220 do 1007 iy=1,3
2221 npartn=IDINT(xnpart(ianp,iy)-xnpp(iy))
2222 if(npartn.ne.0) then
2223 v2pevt(iy)=(v2p(ianp,iy)-v2pp(iy))/npartn
2224 v2psum(iy)=v2psum(iy)+v2pevt(iy)
2225 v2p2sm(iy)=v2p2sm(iy)+v2pevt(iy)**2
2226 v2pp(iy)=v2p(ianp,iy)
2227 xnpp(iy)=xnpart(ianp,iy)
2228 endif
2229 1007 continue
2230 write(49, 160) iaevt,v2pevt
2231 endif
2232 goto 101
2233 endif
2234 1008 continue
2235clin-11/28/00 for animation:
2236 101 if(nevent.eq.1) then
2237 do 110 nt = 1, ntmax
2238 time1=dble(nt*dt)
2239 time2=dble((nt+1)*dt)
2240 if (t2time.lt.time2.and.t2time.ge.time1) then
2241 if(nt.le.itanim) return
2242c write(92,*) t2time
2243 iform=0
2244 do 1009 I=1,mul
2245c write out parton info only for formed ones:
2246 if(FT5(I).le.t2time) then
2247 iform=iform+1
2248 endif
2249 1009 continue
2250c write(92,*) iform
2251 do 120 I=1,mul
2252 if(FT5(I).le.t2time) then
2253clin-11/29/00-ctest off calculate parton coordinates after propagation:
2254c gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2255c gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2256c gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2257c write(92,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2258c write(92,180) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i),
2259c 1 PX5(i),PY5(i),PZ5(i),E5(i)
2260ctest-end
2261 endif
2262 120 continue
2263 itanim=nt
2264 endif
2265 110 continue
2266 endif
2267c
2268 140 format(i10,4(2x,f7.2))
2269 160 format(i10,3(2x,f9.5))
2270c 180 format(i6,8(1x,f7.2))
2271clin-5/17/01 calculate v2 for active partons (which have not frozen out):
2272c idd=3, called at end of zpc.f:
2273 else if(idd.eq.3) then
2274 do 1010 ianp=1,30
2275 if(iscatt(ianp).eq.0) tscatt(ianp+1)=tscatt(ianp)
2276 1010 continue
2277 do 350 I=1,mul
2278 ypartn=0.5d0*dlog((E5(i)+PZ5(i)+1.d-8)
2279 1 /(E5(i)-PZ5(i)+1.d-8))
2280 pt2=PX5(I)**2+PY5(I)**2
2281 iloop=1
2282 if(dabs(ypartn).le.1d0) then
2283 iloop=2
2284 if(dabs(ypartn).le.0.5d0) then
2285 iloop=3
2286 endif
2287 endif
2288c
2289 do 325 ianp=1,30
2290 if(iscatt(ianp).ne.0) then
2291 if(FT5(I).lt.tscatt(ianp+1)
2292 1 .and.FT5(I).ge.tscatt(ianp)) then
2293 do 1011 iy=1,iloop
2294 if(pt2.gt.0.) then
2295 v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2296 v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn
2297 v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2
2298 endif
2299 xperp2=GX5(I)**2+GY5(I)**2
2300 if(xperp2.gt.0.) s2pfrz(ianp,iy)=
2301 1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2302 etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2303 xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0
2304ctest off density:
2305c etpfrz(ianp,iy)=etpfrz(ianp,iy)
2306c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2307csp07/05
2308 if(ITYP5(I).eq.1.or.ITYP5(I).eq.2)then
2309 etlfrz(ianp,iy)=etlfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2310 xnlfrz(ianp,iy)=xnlfrz(ianp,iy)+1d0
2311 elseif(ITYP5(I).eq.3)then
2312 etsfrz(ianp,iy)=etsfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2313 xnsfrz(ianp,iy)=xnsfrz(ianp,iy)+1d0
2314 endif
2315csp07/05 end
2316 1011 continue
2317c parton freezeout info taken, proceed to next parton:
2318 goto 350
2319 endif
2320 endif
2321 325 continue
2322 350 continue
2323c idd=2: calculate average partonic elliptic flow, called from artdri.f,
2324 else if(idd.eq.2) then
2325 do 1012 i=1,3
2326cms write(nfile(i),*) ' tsp, v2p, v2p2, '//
2327cms 1 ' s2p, etp, xmult, nevt, xnctot'
2328cms write ((nfile(i)+1),*) ' tsp, v2pf, v2pf2, '//
2329cms 1 ' s2pf, etpf, xnform, xcrate'
2330cms write ((nfile(i)+2),*) ' tsp, v2pa, v2pa2, '//
2331cms 1 ' s2pa, etpa, xmult_ap, xnform, nevt'
2332cms write ((nfile(i)+3),*) ' tsph, v2ph, v2ph2, '//
2333cms 1 ' s2ph, etph, xmult_(ap/2+h),xmult_ap/2,nevt'
2334csp
2335cms write(nlfile(i),*) ' tsp, v2, s2, etp, xmul'
2336cms write(nsfile(i),*) ' tsp, v2, s2, etp, xmul'
2337cms write(nlfile(i)+1,*) ' tsp, v2, s2, etp, xmul'
2338cms write(nsfile(i)+1,*) ' tsp, v2, s2, etp, xmul'
2339c
2340 1012 continue
2341clin-3/30/00 ensemble average & variance of v2 (over particles in all events):
2342 do 150 ii=1, 30
2343 if(nevt(ii).eq.0) goto 150
2344 do 1014 iy=1,3
2345 if(xnpart(ii,iy).gt.1) then
2346 v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy)
2347 v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy)
2348 1 -v2p(ii,iy)**2)/(xnpart(ii,iy)-1))
2349 s2p(ii,iy)=s2p(ii,iy)/xnpart(ii,iy)
2350c xmult and etp are multiplicity and et for an averaged event:
2351 xmult=dble(xnpart(ii,iy)/dble(nevt(ii)))
2352 etp(ii,iy)=etp(ii,iy)/dble(nevt(ii))
2353csp
2354 etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2355 etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2356c
2357 xnctot=0d0
2358 do 1013 inum=1,ii
2359 xnctot=xnctot+xncoll(inum)
2360 1013 continue
2361 if(nevt(1).ne.0) xnctot=xnctot/nevt(1)
2362cms write (nfile(iy),200) tsp(ii),v2p(ii,iy),
2363cms 1 v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot
2364 endif
2365 if(nevt(ii).ne.0)
2366 1 xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii)
2367c
2368 if(xnpf(ii,iy).gt.1) then
2369 v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy)
2370 v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy)
2371 1 -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1))
2372 s2pf(ii,iy)=s2pf(ii,iy)/xnpf(ii,iy)
2373 xnform=dble(xnpf(ii,iy)/dble(nevt(ii)))
2374 etpf(ii,iy)=etpf(ii,iy)/dble(nevt(ii))
2375csp
2376 etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2377 etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2378c
2379cms write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy),
2380cms 1 v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate
2381 endif
2382csp
2383 if(xnpl(ii,iy).gt.1) then
2384 v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy)
2385 s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy)
2386 xmult=dble(xnpl(ii,iy)/dble(nevt(ii)))
2387 etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2388cms write (nlfile(iy),201) tsp(ii),v2pl(ii,iy),
2389cms 1 s2pl(ii,iy),etpl(ii,iy),xmult
2390 endif
2391 if(xnps(ii,iy).gt.1) then
2392 v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy)
2393 s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy)
2394 xmult=dble(xnps(ii,iy)/dble(nevt(ii)))
2395 etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2396cms write (nsfile(iy),201) tsp(ii),v2ps(ii,iy),
2397cms 1 s2ps(ii,iy),etps(ii,iy),xmult
2398 endif
2399 if(xnplf(ii,iy).gt.1) then
2400 v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy)
2401 s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy)
2402 xmult=dble(xnplf(ii,iy)/dble(nevt(ii)))
2403 etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2404cms write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy),
2405cms 1 s2plf(ii,iy),etplf(ii,iy),xmult
2406 endif
2407 if(xnpsf(ii,iy).gt.1) then
2408 v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy)
2409 s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy)
2410 xmult=dble(xnpsf(ii,iy)/dble(nevt(ii)))
2411 etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2412cms write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy),
2413cms 1 s2psf(ii,iy),etpsf(ii,iy),xmult
2414 endif
2415csp-end
2416 1014 continue
2417 150 continue
2418csp07/05 initial & final mt distrb
2419 scalei=0d0
2420 scalef=0d0
2421 if(nevt(1).ne.0) SCALEi = 1d0 / dble(nevt(1)) / BMT
2422 if(nevt(30).ne.0) SCALEf = 1d0 / dble(nevt(30)) / BMT
2423 do 1016 iy=2,3
2424 yra = 1d0
2425 if(iy .eq. 2)yra = 2d0
2426 do 1015 i=1,50
2427cms WRITE(nmfile(iy),251) BMT*dble(I - 0.5),
2428cms & SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra,
2429cms & SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra
2430 1015 continue
2431 1016 continue
2432csp07/05 end
2433clin-3/30/00 event-by-event average & variance of v2:
2434 if(nevt(30).ge.1) then
2435 do 1017 iy=1,3
2436 v2pavg(iy)=v2psum(iy)/nevt(30)
2437 v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2
2438 if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0)
2439 1017 continue
2440 write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg
2441 write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p
2442 endif
2443clin-11/28/00 for animation:
2444 if(nevent.eq.1) then
2445 do 1018 I=1,mul
2446 if(FT5(I).le.t2time) then
2447 write(93,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2448 endif
2449 1018 continue
2450clin-11/29/00 signal the end of animation file:
2451c write(91,*) -10.
2452c write(91,*) 0
2453c write(92,*) -10.
2454c write(92,*) 0
2455c close (91)
2456c close (92)
2457 close (93)
2458 endif
2459clin-5/18/01 calculate v2 for active partons:
2460 do 450 ianp=1,30
2461 do 400 iy=1,3
2462 v2pact=0d0
2463 v2p2ac=0d0
2464 s2pact=0d0
2465 etpact=0d0
2466 xnacti=0d0
2467 if(xnpf(ianp,iy).gt.1) then
2468c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2469 v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2470 v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2471 1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2472 s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2473 etpact=etpf(ianp,iy)*dble(nevt(ianp))
2474 xnpact=xnpf(ianp,iy)
2475c
2476 do 1019 kanp=1,ianp
2477 v2pact=v2pact-v2pfrz(kanp,iy)
2478 v2p2ac=v2p2ac-v2p2fz(kanp,iy)
2479 s2pact=s2pact-s2pfrz(kanp,iy)
2480 etpact=etpact-etpfrz(kanp,iy)
2481 xnpact=xnpact-xnpfrz(kanp,iy)
2482 1019 continue
2483c save the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2484 v2ph=v2pact
2485 v2ph2=v2p2ac
2486 s2ph=s2pact
2487 etph=etpact
2488 xnp2=xnpact/2d0
2489c
2490 if(xnpact.gt.1.and.nevt(ianp).ne.0) then
2491 v2pact=v2pact/xnpact
2492 v2p2ac=dsqrt((v2p2ac/xnpact
2493 1 -v2pact**2)/(xnpact-1))
2494 s2pact=s2pact/xnpact
2495 xnacti=dble(xnpact/dble(nevt(ianp)))
2496 etpact=etpact/dble(nevt(ianp))
2497cms write (nfile(iy)+2, 250) tsp(ianp),v2pact,
2498cms 1 v2p2ac,s2pact,etpact,xnacti,
2499cms 2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp)
2500 endif
2501 endif
2502c To calculate combined v2 for active partons plus formed hadrons,
2503c add the sum of v2h, v2h2, s2h, eth, and xnh for formed hadrons:
2504c scale the hadron part in case nevt(ianp) != nevent:
2505 shadr=dble(nevt(ianp))/dble(nevent)
2506 ianh=ianp
2507 v2ph=v2ph+v2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2508 v2ph2=v2ph2+(v2h2(ianh,iy)**2*(xnhadr(ianh,iy)-1)
2509 1 +v2h(ianh,iy)**2)*xnhadr(ianh,iy)*shadr
2510 s2ph=s2ph+s2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2511 etph=etph+eth(ianh,iy)*dble(nevent)*shadr
2512 xnph=xnpact+xnhadr(ianh,iy)*shadr
2513 xnp2h=xnp2+xnhadr(ianh,iy)*shadr
2514 if(xnph.gt.1) then
2515 v2ph=v2ph/xnph
2516 v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1))
2517 s2ph=s2ph/xnph
2518 etph=etph/dble(nevt(ianp))
2519 xnp2=xnp2/dble(nevt(ianp))
2520 xnp2h=xnp2h/dble(nevent)
2521cms if(tsp(ianp).le.(ntmax*dt))
2522cms 1 write (nfile(iy)+3, 250) tsp(ianp),v2ph,
2523cms 2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp)
2524 endif
2525c
2526 400 continue
2527 450 continue
2528 do 550 ianp=1,30
2529 do 500 iy=1,3
2530 v2pact=0d0
2531 v2p2ac=0d0
2532 s2pact=0d0
2533 etpact=0d0
2534 xnacti=0d0
2535c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2536 v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2537 v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2538 1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2539 s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2540 etpact=etpf(ianp,iy)*dble(nevt(ianp))
2541 xnpact=xnpf(ianp,iy)
2542 500 continue
2543 550 continue
2544cms close (620)
2545cms close (630)
2546 do 1021 nf=1,3
2547 do 1020 ifile=0,3
2548cms close(nfile(nf)+ifile)
2549 1020 continue
2550 1021 continue
2551 do 1022 nf=1,3
2552cms close(740+nf)
2553 1022 continue
2554 endif
2555cyy 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2)
2556cyy 210 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2))
2557 240 format(a30,3(2x,f9.5))
2558cyy 250 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6)
2559csp
2560cyy 201 format(2x,f5.2,4(2x,f9.2))
2561cyy 251 format(5e15.5)
2562c
2563 return
2564 end
2565
2566c=======================================================================
2567c Calculate flow from formed hadrons, called by art1e.f:
2568c Note: numbers in art not in double precision!
2569 subroutine flowh(ct)
2570 PARAMETER (MAXSTR=150001, MAXR=1)
2571 dimension tsh(31)
2572 DOUBLE PRECISION v2h,xnhadr,eth,v2h2,s2h
2573 DOUBLE PRECISION v2hp,xnhadp,v2hsum,v2h2sm,v2hevt(3)
2574 DOUBLE PRECISION pt2, v2hadr
2575 COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2576 1 v2h2(30,3),s2h(30,3)
2577cc SAVE /hflow/
2578 common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2579cc SAVE /ebe/
2580 common /lastt/itimeh,bimp
2581cc SAVE /lastt/
2582 COMMON /RUN/ NUM
2583cc SAVE /RUN/
2584 COMMON /AA/ R(3,MAXSTR)
2585cc SAVE /AA/
2586 COMMON /BB/ P(3,MAXSTR)
2587cc SAVE /BB/
2588 COMMON /CC/ E(MAXSTR)
2589cc SAVE /CC/
2590 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2591cc SAVE /EE/
2592 COMMON /RR/ MASSR(0:MAXR)
2593cc SAVE /RR/
2594 common/anim/nevent,isoft,isflag,izpc
2595cc SAVE /anim/
2596 COMMON /AREVT/ IAEVT, IARUN, MISS
2597cc SAVE /AREVT/
2598 SAVE
2599c
2600 do 1001 ii = 1, 31
2601 tsh(ii)=float(ii-1)
2602 1001 continue
2603c
2604 do 1004 ianh = 1, 30
2605 if ((ct+0.0001).lt.tsh(ianh+1)
2606 1 .and.(ct+0.0001).ge.tsh(ianh)) then
2607 if(ianh.eq.itimeh) goto 101
2608 IA = 0
2609 DO 1002 J = 1, NUM
2610 mult=MASSR(J)
2611 IA = IA + MASSR(J - 1)
2612 DO 100 IC = 1, mult
2613 I = IA + IC
2614c 5/04/01 exclude leptons and photons:
2615 if(iabs(LB(I)-10000).lt.100) goto 100
2616 px=p(1,i)
2617 py=p(2,i)
2618 pt2=dble(px)**2+dble(py)**2
2619c 2/18/00 Note: e(i) gives the mass in ART:
2620 ene=sqrt(e(i)**2+sngl(pt2)+p(3,i)**2)
2621 RAP=0.5*alog((ene+p(3,i))/(ene-p(3,i)))
2622ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2623c rt_now=sqrt(r(1,i)**2+r(2,i)**2)
2624c gz_now=r(3,i)
2625c zmax=0.1*ct
2626c volume=3.1416*(2.**2)*2*zmax
2627c if(rt_now.gt.2.or.abs(gz_now).gt.zmax)
2628c 1 goto 100
2629 iloop=1
2630 if(abs(rap).le.1) then
2631 iloop=2
2632 if(abs(rap).le.0.5) then
2633 iloop=3
2634 endif
2635 endif
2636 do 50 iy=1,iloop
2637 if(pt2.gt.0d0) then
2638 v2hadr=(dble(px)**2-dble(py)**2)/pt2
2639 v2h(ianh,iy)=v2h(ianh,iy)+v2hadr
2640 v2h2(ianh,iy)=v2h2(ianh,iy)+v2hadr**2
2641 if(dabs(v2hadr).gt.1d0)
2642 1 write(1,*) 'v2hadr>1',v2hadr,px,py
2643 endif
2644 xperp2=r(1,I)**2+r(2,I)**2
2645 if(xperp2.gt.0.)
2646 1 s2h(ianh,iy)=s2h(ianh,iy)+dble((r(1,I)**2-r(2,I)**2)/xperp2)
2647 eth(ianh,iy)=eth(ianh,iy)+dble(SQRT(e(i)**2+sngl(pt2)))
2648ctest off density:
2649c eth(ianh,iy)=eth(ianh,iy)
2650c 1 +dble(SQRT(e(i)**2+sngl(pt2)+p(3,i)**2))/volume
2651 xnhadr(ianh,iy)=xnhadr(ianh,iy)+1d0
2652 50 continue
2653 100 continue
2654 1002 CONTINUE
2655 itimeh=ianh
2656clin-5/04/01 ebe v2 variables:
2657 if(ianh.eq.30) then
2658 do 1003 iy=1,3
2659 nhadrn=IDINT(xnhadr(ianh,iy)-xnhadp(iy))
2660 if(nhadrn.ne.0) then
2661 v2hevt(iy)=(v2h(ianh,iy)-v2hp(iy))/dble(nhadrn)
2662 v2hsum(iy)=v2hsum(iy)+v2hevt(iy)
2663 v2h2sm(iy)=v2h2sm(iy)+v2hevt(iy)**2
2664 v2hp(iy)=v2h(ianh,iy)
2665 xnhadp(iy)=xnhadr(ianh,iy)
2666 endif
2667 1003 continue
2668 write(88, 160) iaevt,v2hevt
2669 endif
2670 goto 101
2671 endif
2672 1004 continue
2673 160 format(i10,3(2x,f9.5))
2674clin-11/27/00 for animation:
2675 101 if(nevent.eq.1) then
2676 IA = 0
2677 do 1005 J = 1, NUM
2678 mult=MASSR(J)
2679 IA = IA + MASSR(J - 1)
2680c write(91,*) ct
2681c write(91,*) mult
2682 DO 150 IC = 1, mult
2683 I = IA + IC
2684c write(91,210) LB(i),r(1,i),r(2,i),r(3,i),
2685c 1 p(1,i),p(2,i),p(3,i),e(i)
2686 150 continue
2687 1005 continue
2688 return
2689 endif
2690c 210 format(i6,7(1x,f8.3))
2691 return
2692 end
2693
2694c=======================================================================
2695 subroutine flowh0(NEVNT,idd)
2696c
2697 dimension tsh(31)
2698 DOUBLE PRECISION v2h,xnhadr,eth,v2h2,s2h
2699 DOUBLE PRECISION v2hp,xnhadp,v2hsum,v2h2sm,
2700 1 v2havg(3),varv2h(3)
2701 COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2702 1 v2h2(30,3),s2h(30,3)
2703cc SAVE /hflow/
2704 common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2705cc SAVE /ebe/
2706 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
2707cc SAVE /input1/
2708 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
2709 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
2710cc SAVE /INPUT2/
2711 common /lastt/itimeh,bimp
2712cc SAVE /lastt/
2713 SAVE
2714
2715c idd=0: initialization for flow analysis, called by artdri.f::
2716 if(idd.eq.0) then
2717 itimeh=0
2718c
2719 do 1001 ii = 1, 31
2720 tsh(ii)=float(ii-1)
2721 1001 continue
2722c
2723 do 1003 ii=1,30
2724 do 1002 iy=1,3
2725 v2h(ii,iy)=0d0
2726 xnhadr(ii,iy)=0d0
2727 eth(ii,iy)=0d0
2728 v2h2(ii,iy)=0d0
2729 s2h(ii,iy)=0d0
2730 1002 continue
2731 1003 continue
2732 do 1004 iy=1,3
2733 v2hp(iy)=0d0
2734 xnhadp(iy)=0d0
2735 v2hsum(iy)=0d0
2736 v2h2sm(iy)=0d0
2737 if(iy.eq.1) then
2738 nunit=59
2739 elseif(iy.eq.2) then
2740 nunit=68
2741 else
2742 nunit=69
2743 endif
2744 write(nunit,*) ' tsh, v2h, v2h2, s2h, '//
2745 1 ' eth, xmulth'
2746 1004 continue
2747c idd=2: calculate average hadronic elliptic flow, called by artdri.f:
2748 else if(idd.eq.2) then
2749 do 100 ii=1, 30
2750 do 1005 iy=1,3
2751 if(xnhadr(ii,iy).eq.0) then
2752 xmulth=0.
2753 elseif(xnhadr(ii,iy).gt.1) then
2754 v2h(ii,iy)=v2h(ii,iy)/xnhadr(ii,iy)
2755 eth(ii,iy)=eth(ii,iy)/dble(NEVNT)
2756 v2h2(ii,iy)=dsqrt((v2h2(ii,iy)/xnhadr(ii,iy)
2757 1 -v2h(ii,iy)**2)/(xnhadr(ii,iy)-1))
2758 s2h(ii,iy)=s2h(ii,iy)/xnhadr(ii,iy)
2759 xmulth=sngl(xnhadr(ii,iy)/NEVNT)
2760 endif
2761 if(iy.eq.1) then
2762 nunit=59
2763 elseif(iy.eq.2) then
2764 nunit=68
2765 else
2766 nunit=69
2767 endif
2768 if(tsh(ii).le.(ntmax*dt))
2769 1 write (nunit,200) tsh(ii),v2h(ii,iy),
2770 2 v2h2(ii,iy),s2h(ii,iy),eth(ii,iy),xmulth
2771 1005 continue
2772 100 continue
2773c event-by-event average & variance of v2h:
2774 do 1006 iy=1,3
2775 v2havg(iy)=v2hsum(iy)/dble(NEVNT)
2776 varv2h(iy)=dsqrt(v2h2sm(iy)/dble(NEVNT)-v2havg(iy)**2)
2777 1006 continue
2778 write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): avg=', v2havg
2779 write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): var=', varv2h
2780 endif
2781 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2))
2782 240 format(a30,3(2x,f9.5))
2783 return
2784 end
2785
2786c=======================================================================
2787c 2/23/00 flow from all initial hadrons just before entering ARTMN:
2788 subroutine iniflw(NEVNT,idd)
2789 PARAMETER (MAXSTR=150001, MAXR=1)
2790 DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult,
2791 1 v2bi,s2bi,xbmult
2792 COMMON /RUN/ NUM
2793cc SAVE /RUN/
2794 COMMON /ARERC1/MULTI1(MAXR)
2795cc SAVE /ARERC1/
2796 COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2797 & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
2798 & FT1(MAXSTR, MAXR),
2799 & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2800 & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2801cc SAVE /ARPRC1/
2802 COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
2803cc SAVE /iflow/
2804 SAVE
2805c
2806 if(idd.eq.0) then
2807 v2i=0d0
2808 eti=0d0
2809 xmulti=0d0
2810 v2mi=0d0
2811 s2mi=0d0
2812 xmmult=0d0
2813 v2bi=0d0
2814 s2bi=0d0
2815 xbmult=0d0
2816 else if(idd.eq.1) then
2817 do 1002 J = 1, NUM
2818 do 1001 I = 1, MULTI1(J)
2819 ITYP = ITYP1(I, J)
2820c all hadrons:
2821 IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 100
2822 xmulti=xmulti+1.d0
2823 PX = PX1(I, J)
2824 PY = PY1(I, J)
2825 XM = XM1(I, J)
2826 pt2=px**2+py**2
2827 xh=gx1(I,J)
2828 yh=gy1(I,J)
2829 xt2=xh**2+yh**2
2830 if(pt2.gt.0) v2i=v2i+dble((px**2-py**2)/pt2)
2831 eti=eti+dble(SQRT(PX ** 2 + PY ** 2 + XM ** 2))
2832c baryons only:
2833 IF (ITYP .LT. -1000 .or. ITYP .GT. 1000) then
2834 xbmult=xbmult+1.d0
2835 if(pt2.gt.0) v2bi=v2bi+dble((px**2-py**2)/pt2)
2836 if(xt2.gt.0) s2bi=s2bi+dble((xh**2-yh**2)/xt2)
2837c mesons only:
2838 else
2839 xmmult=xmmult+1.d0
2840 if(pt2.gt.0) v2mi=v2mi+dble((px**2-py**2)/pt2)
2841 if(xt2.gt.0) s2mi=s2mi+dble((xh**2-yh**2)/xt2)
2842 endif
2843 100 continue
2844 1001 continue
2845 1002 continue
2846 else if(idd.eq.2) then
2847 if(xmulti.ne.0) v2i=v2i/xmulti
2848 eti=eti/dble(NEVNT)
2849 xmulti=xmulti/dble(NEVNT)
2850 if(xmmult.ne.0) then
2851 v2mi=v2mi/xmmult
2852 s2mi=s2mi/xmmult
2853 endif
2854 xmmult=xmmult/dble(NEVNT)
2855 if(xbmult.ne.0) then
2856 v2bi=v2bi/xbmult
2857 s2bi=s2bi/xbmult
2858 endif
2859 xbmult=xbmult/dble(NEVNT)
2860 endif
2861c
2862 return
2863 end
2864
2865c=======================================================================
2866c 2/25/00 dN/dt analysis for production (before ZPCMN)
2867c and freezeout (right after ZPCMN) for all partons.
2868 subroutine frztm(NEVNT,idd)
2869c
2870 implicit double precision (a-h, o-z)
2871 PARAMETER (MAXPTN=400001)
2872 dimension tsf(31)
2873 COMMON /PARA1/ MUL
2874cc SAVE /PARA1/
2875 COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
2876 & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
2877 & XMASS0(MAXPTN), ITYP0(MAXPTN)
2878cc SAVE /prec1/
2879 COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
2880 & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
2881 & XMASS5(MAXPTN), ITYP5(MAXPTN)
2882cc SAVE /prec2/
2883 COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
2884 & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
2885cc SAVE /frzout/
2886 SAVE
2887 data tsf/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2888 & 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
2889 & 2 , 3, 4, 5, 6, 7, 8, 9, 10, 20, 30/
2890c
2891 if(idd.eq.0) then
2892 do 1001 ii=1,30
2893 xnprod(ii)=0d0
2894 etprod(ii)=0d0
2895 xnfrz(ii)=0d0
2896 etfrz(ii)=0d0
2897 dnprod(ii)=0d0
2898 detpro(ii)=0d0
2899 dnfrz(ii)=0d0
2900 detfrz(ii)=0d0
2901 1001 continue
2902cms OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN')
2903cms OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN')
2904 else if(idd.eq.1) then
2905 DO 100 ip = 1, MUL
2906 do 1002 ii=1,30
2907 eth0=dSQRT(PX0(ip)**2+PY0(ip)**2+XMASS0(ip)**2)
2908 eth2=dSQRT(PX5(ip)**2+PY5(ip)**2+XMASS5(ip)**2)
2909c total number and Et produced by time tsf(ii):
2910 if (ft0(ip).lt.tsf(ii+1)) then
2911 xnprod(ii)=xnprod(ii)+1d0
2912 etprod(ii)=etprod(ii)+eth0
2913c number and Et produced from time tsf(ii) to tsf(ii+1):
2914 if (ft0(ip).ge.tsf(ii)) then
2915 dnprod(ii)=dnprod(ii)+1d0
2916 detpro(ii)=detpro(ii)+eth0
2917 endif
2918 endif
2919c total number and Et freezed out by time tsf(ii):
2920 if (FT5(ip).lt.tsf(ii+1)) then
2921 xnfrz(ii)=xnfrz(ii)+1d0
2922 etfrz(ii)=etfrz(ii)+eth2
2923c number and Et freezed out from time tsf(ii) to tsf(ii+1):
2924 if (FT5(ip).ge.tsf(ii)) then
2925 dnfrz(ii)=dnfrz(ii)+1d0
2926 detfrz(ii)=detfrz(ii)+eth2
2927 endif
2928 endif
2929 1002 continue
2930 100 continue
2931 else if(idd.eq.2) then
2932cms write (86,*) ' t, np, dnp/dt, etp '//
2933cms 1 ' detp/dt'
2934cms write (87,*) ' t, nf, dnf/dt, etf '//
2935cms 1 ' detf/dt'
2936 do 1003 ii=1,30
2937 xnp=xnprod(ii)/dble(NEVNT)
2938 xnf=xnfrz(ii)/dble(NEVNT)
2939 etp=etprod(ii)/dble(NEVNT)
2940 etf=etfrz(ii)/dble(NEVNT)
2941 dxnp=dnprod(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2942 dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2943 detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2944 detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2945cms write (86, 200)
2946cms 1 tsf(ii+1),xnp,dxnp,etp,detp
2947cms write (87, 200)
2948cms 1 tsf(ii+1),xnf,dxnf,etf,detf
2949 1003 continue
2950 endif
2951cyy 200 format(2x,f9.2,4(2x,f10.2))
2952c
2953 return
2954 end
2955
2956c=======================================================================
2957clin-6/2009 write out initial minijet information
2958c before propagating to its formation time:
2959 subroutine minijet_out(BB)
2960 PARAMETER (MAXSTR=150001)
2961 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2962 COMMON/hjcrdn/YP(3,300),YT(3,300)
2963 COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2964 & PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2965 & PJPM(300,500),NTJ(300),KFTJ(300,500),
2966 & PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2967 & PJTE(300,500),PJTM(300,500)
2968 COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2969 & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2970 & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2971 COMMON /AREVT/ IAEVT, IARUN, MISS
2972 common /para7/ ioscar,nsmbbbar,nsmmeson
2973 common/phidcy/iphidcy,pttrig,ntrig,maxmiss
2974 SAVE
2975 ntrig=0
2976 do I = 1, IHNT2(1)
2977 do J = 1, NPJ(I)
2978 pt=sqrt(PJPX(I,J)**2+PJPY(I,J)**2)
2979 if(pt.ge.pttrig) ntrig=ntrig+1
2980 enddo
2981 enddo
2982 do I = 1, IHNT2(3)
2983 do J = 1, NTJ(I)
2984 pt=sqrt(PJTX(I,J)**2+PJTY(I,J)**2)
2985 if(pt.ge.pttrig) ntrig=ntrig+1
2986 enddo
2987 enddo
2988 do I = 1, NSG
2989 do J = 1, NJSG(I)
2990 pt=sqrt(PXSG(I,J)**2+PYSG(I,J)**2)
2991 if(pt.ge.pttrig) ntrig=ntrig+1
2992 enddo
2993 enddo
2994c Require at least 1 initial minijet parton above the trigger Pt value:
2995 if(ntrig.eq.0) return
2996
2997c.....transfer data from HIJING to ZPC
2998 if(ioscar.eq.3) write(96,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
2999 DO 1008 I = 1, IHNT2(1)
3000 DO 1007 J = 1, NPJ(I)
3001 ityp=KFPJ(I,J)
3002c write out not only gluons:
3003c if(ityp.ne.21) goto 1007
3004 gx=YP(1,I)+0.5*BB
3005 gy=YP(2,I)
3006 gz=0d0
3007 ft=0d0
3008 px=PJPX(I,J)
3009 py=PJPY(I,J)
3010 pz=PJPZ(I,J)
3011 xmass=PJPM(I,J)
3012 if(ioscar.eq.3) then
3013 if(amax1(abs(gx),abs(gy),
3014 1 abs(gz),abs(ft)).lt.9999) then
3015 write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3016 else
3017 write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3018 endif
3019 endif
3020 1007 CONTINUE
3021 1008 CONTINUE
3022 DO 1010 I = 1, IHNT2(3)
3023 DO 1009 J = 1, NTJ(I)
3024 ityp=KFTJ(I,J)
3025c if(ityp.ne.21) goto 1009
3026 gx=YT(1,I)-0.5*BB
3027 gy=YT(2,I)
3028 gz=0d0
3029 ft=0d0
3030 px=PJTX(I,J)
3031 py=PJTY(I,J)
3032 pz=PJTZ(I,J)
3033 xmass=PJTM(I,J)
3034 if(ioscar.eq.3) then
3035 if(amax1(abs(gx),abs(gy),
3036 1 abs(gz),abs(ft)).lt.9999) then
3037 write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3038 else
3039 write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3040 endif
3041 endif
3042 1009 CONTINUE
3043 1010 CONTINUE
3044 DO 1012 I = 1, NSG
3045 DO 1011 J = 1, NJSG(I)
3046 ityp=K2SG(I,J)
3047c if(ityp.ne.21) goto 1011
3048 gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2)))
3049 gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2)))
3050 gz=0d0
3051 ft=0d0
3052 px=PXSG(I,J)
3053 py=PYSG(I,J)
3054 pz=PZSG(I,J)
3055 xmass=PMSG(I,J)
3056 if(ioscar.eq.3) then
3057 if(amax1(abs(gx),abs(gy),
3058 1 abs(gz),abs(ft)).lt.9999) then
3059 write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3060 else
3061 write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3062 endif
3063 endif
3064 1011 CONTINUE
3065 1012 CONTINUE
3066 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,f8.2),2(2x,f2.0),2x,I2)
3067 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,e8.2),2(2x,f2.0),2x,I2)
3068c
3069 return
3070 end
3071
3072c=======================================================================
3073clin-6/2009 embed back-to-back high-Pt quark/antiquark pair
3074c via embedding back-to-back high-Pt pion pair then melting the pion pair
3075c by generating the internal quark and antiquark momentum parallel to
3076c the pion momentum (in order to produce a high-Pt and a low Pt parton):
3077 subroutine embedHighPt
3078 PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140)
3079 common/embed/iembed,pxqembd,pyqembd,xembd,yembd
3080 COMMON/RNDF77/NSEED
3081 COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
3082 COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
3083 COMMON /ARPRC/ ITYPAR(MAXSTR),
3084 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
3085 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
3086 & XMAR(MAXSTR)
3087cwei DOUBLE PRECISION PATT
3088 SAVE
3089c
3090 if(iembed.ne.1) return
3091 ptq=sqrt(pxqembd**2+pyqembd**2)
3092 if(ptq.lt.(pichmass/2.)) then
3093 print *, 'Embedded quark transverse momentum is too small'
3094 stop
3095 endif
3096c Randomly embed u/ubar or d/dbar at high Pt:
3097 idqembd=1+int(2*RANART(NSEED))
3098c Flavor content for the charged pion that contains the leading quark:
3099 if(idqembd.eq.1) then
3100 idqsoft=-2
3101 idpi=-211
3102 elseif(idqembd.eq.2) then
3103 idqsoft=-1
3104 idpi=211
3105 else
3106 print *, 'Wrong quark flavor embedded'
3107 stop
3108 endif
3109c Caculate transverse momentum of the parent charged pion:
3110 xmq=ulmass(idqembd)
3111 xmqsoft=ulmass(idqsoft)
3112 ptpi=((pichmass**2+xmq**2-xmqsoft**2)*ptq
3113 1 -sqrt((xmq**2+ptq**2)*(pichmass**4
3114 2 -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2)))
3115 3 /(2.*xmq**2)
3116 pxpi=ptpi*pxqembd/ptq
3117 pypi=ptpi*pyqembd/ptq
3118c Embedded quark/antiquark are assumed to have pz=0:
3119 pzpi=0.
3120c Insert the two parent charged pions,
3121c ipion=1 for the pion containing the leading quark,
3122c ipion=2 for the pion containing the leading antiquark of the same flavor:
3123 do ipion=1,2
3124 if(ipion.eq.2) then
3125 idpi=-idpi
3126 pxpi=-pxpi
3127 pypi=-pypi
3128 pzpi=-pzpi
3129 endif
3130 NATT=NATT+1
3131 KATT(NATT,1)=idpi
3132 KATT(NATT,2)=40
3133 KATT(NATT,3)=0
3134 PATT(NATT,1)=pxpi
3135 PATT(NATT,2)=pypi
3136 PATT(NATT,3)=pzpi
3137 PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2)
3138 EATT=EATT+PATT(NATT,4)
3139 GXAR(NATT)=xembd
3140 GYAR(NATT)=yembd
3141 GZAR(NATT)=0.
3142 FTAR(NATT)=0.
3143 ITYPAR(NATT)=KATT(NATT,1)
3144 PXAR(NATT)=PATT(NATT,1)
3145 PYAR(NATT)=PATT(NATT,2)
3146 PZAR(NATT)=PATT(NATT,3)
3147 PEAR(NATT)=PATT(NATT,4)
3148 XMAR(NATT)=pichmass
3149 enddo
3150c
3151 return
3152 end