]>
Commit | Line | Data |
---|---|---|
1 | c.................... linana.f | |
2 | c======================================================================= | |
3 | c 10/26/01 update freezeout positions in case of interactions: | |
4 | clin-3/2009 Note: freezeout spacetime values cannot be trusted for K0S & K0L | |
5 | c as K0S/K0L are converted from K+/K- by hand at the end of hadron cascade. | |
6 | cms | |
7 | cms dlw & gsfs Commented out output file writing | |
8 | cms | |
9 | subroutine hbtout(nnew,nt,ntmax) | |
10 | c | |
11 | PARAMETER (MAXSTR=150001,MAXR=1) | |
12 | clin-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 | |
16 | cc SAVE /para7/ | |
17 | COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast | |
18 | cc SAVE /hbt/ | |
19 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
20 | cc SAVE /input1/ | |
21 | COMMON /AA/ R(3,MAXSTR) | |
22 | cc SAVE /AA/ | |
23 | COMMON /BB/ P(3,MAXSTR) | |
24 | cc SAVE /BB/ | |
25 | COMMON /CC/ E(MAXSTR) | |
26 | cc SAVE /CC/ | |
27 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
28 | cc SAVE /EE/ | |
29 | common /lastt/itimeh,bimp | |
30 | cc SAVE /lastt/ | |
31 | COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR) | |
32 | cc SAVE /tdecay/ | |
33 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
34 | cc SAVE /AREVT/ | |
35 | common/snn/efrm,npart1,npart2 | |
36 | cc SAVE /snn/ | |
37 | COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP | |
38 | cc 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) | |
43 | clin-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 | |
48 | c | |
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 | |
55 | c for each of the particles, search the freezeout record (common /hbt/) | |
56 | c 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 | |
65 | clin-5/2008 modified below to the above in case we have perturbative particles: | |
66 | c 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) | |
70 | c 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) | |
76 | c find particles with dp=0 and dr<0.01, considered to be those | |
77 | c without any interactions during this timestep, | |
78 | c thus keep their last positions and time: | |
79 | if(dr.le.0.01) then | |
80 | lastkp(iplast)=1 | |
81 | newkp(ip)=1 | |
82 | c if(lb(ip).eq.41) then | |
83 | c write(95,*) 'nt,ip,px,x=',nt,ip,p(1,ip),r(1,ip),ftsv(ip) | |
84 | c write(95,*) 'xnew=',xnew(1),xnew(2),xnew(3),xlast(4,ip) | |
85 | c endif | |
86 | clin-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 | |
94 | c for current particles with interactions, fill their current info in | |
95 | c 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 | |
100 | ctest off: write collision info | |
101 | c if(lb(ip).eq.41) then | |
102 | c write(95,*) 'nt,lb(ip)=',nt,lb(ip) | |
103 | c write(95,*) ' last p=',plast(1,iplast), | |
104 | c 1 plast(2,iplast),plast(3,iplast),plast(4,iplast) | |
105 | c write(95,*) ' after p=',p(1,ip),p(2,ip),p(3,ip),e(ip) | |
106 | c write(95,*) 'after x=',r(1,ip),r(2,ip),r(3,ip),ftsv(ip) | |
107 | c endif | |
108 | c | |
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 | |
113 | c | |
114 | if(nt.eq.ntmax) then | |
115 | c freezeout time for decay daughters at the last timestep | |
116 | c 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) | |
119 | c freezeout time for particles unformed at the next-to-last timestep | |
120 | c 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 | |
131 | clin-5/2008: | |
132 | dplast(iplast)=dpertp(ip) | |
133 | goto 150 | |
134 | endif | |
135 | 1005 continue | |
136 | endif | |
137 | 150 continue | |
138 | c if the current particle list is shorter than the freezeout record, | |
139 | c condense the last-collision record by filling new record from 1 to nnew, | |
140 | c 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 | |
156 | clin-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 | |
165 | ctest off look inside each NT timestep (for debugging purpose): | |
166 | c do ip=1,nlast | |
167 | c if(nt.eq.5000) then | |
168 | c write(95,*) ' p ',nt,ip,INVFLV(lblast(ip)),plast(1,ip), | |
169 | c 1 plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip) | |
170 | c write(95,*) ' x ',nt,ip,INVFLV(lblast(ip)),xlast(1,ip), | |
171 | c 1 xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip) | |
172 | c endif | |
173 | c enddo | |
174 | c | |
175 | if(nt.eq.ntmax) then | |
176 | clin-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 | |
184 | c | |
185 | c write(16,190) IAEVT,IARUN,nlast,bimp,npart1,npart2, | |
186 | c 1 NELP,NINP,NELT,NINTHJ | |
187 | c write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2, | |
188 | c 1 NELP,NINP,NELT,NINTHJ | |
189 | clin-5/2008 write out perturbatively-produced particles (deuterons only): | |
190 | c if(idpert.eq.1.or.idpert.eq.2) | |
191 | c 1 write(90,190) IAEVT,IARUN,ndpert,bimp,npart1,npart2, | |
192 | c 2 NELP,NINP,NELT,NINTHJ | |
193 | c do 1007 ip=1,nlast | |
194 | clin-12/14/03 No formation time for spectator projectile or target nucleons, | |
195 | c see ARINI1 in 'amptsub.f': | |
196 | clin-3/2009 To be consistent with new particles produced in hadron cascade | |
197 | c that are limited by the time-resolution (DT) of the hadron cascade, | |
198 | c freezeout time of spectator projectile or target nucleons is written as | |
199 | c DT as they are read at the 1st timestep and then propagated to time DT: | |
200 | c if(plast(1,ip).eq.0.and.plast(2,ip).eq.0 | |
201 | c 1 .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1)) | |
202 | c 2 .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then | |
203 | clin-5/2008 perturbatively-produced particles (currently only deuterons) | |
204 | c are written to ana/ampt_pert.dat (without the column for the mass); | |
205 | c ana/ampt.dat has regularly-produced particles (including deuterons); | |
206 | c these two sets of deuteron data are close to each other(but not the same | |
207 | c because of the bias from triggering the perturbative production); | |
208 | c ONLY use one data set for analysis to avoid double-counting: | |
209 | c if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then | |
210 | c write(16,200) INVFLV(lblast(ip)), plast(1,ip), | |
211 | c 1 plast(2,ip),plast(3,ip),plast(4,ip), | |
212 | c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip), | |
213 | c 3 xlast(4,ip) | |
214 | clin-12/14/03-end | |
215 | c else | |
216 | c if(idpert.eq.1.or.idpert.eq.2) then | |
217 | c write(90,250) INVFLV(lblast(ip)), plast(1,ip), | |
218 | c 1 plast(2,ip),plast(3,ip), | |
219 | c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip), | |
220 | c 3 xlast(4,ip) | |
221 | c else | |
222 | c write(99,*) 'Unexpected perturbative particles' | |
223 | c endif | |
224 | c endif | |
225 | c elseif(amax1(abs(xlast(1,ip)),abs(xlast(2,ip)), | |
226 | c 1 abs(xlast(3,ip)),abs(xlast(4,ip))).lt.9999) then | |
227 | c if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then | |
228 | c write(16,200) INVFLV(lblast(ip)), plast(1,ip), | |
229 | c 1 plast(2,ip),plast(3,ip),plast(4,ip), | |
230 | c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip) | |
231 | c else | |
232 | c if(idpert.eq.1.or.idpert.eq.2) then | |
233 | c write(90,250) INVFLV(lblast(ip)),plast(1,ip), | |
234 | c 1 plast(2,ip),plast(3,ip), | |
235 | c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip), | |
236 | c 3 dplast(ip) | |
237 | c else | |
238 | c write(99,*) 'Unexpected perturbative particles' | |
239 | c endif | |
240 | c endif | |
241 | c else | |
242 | c change format for large numbers: | |
243 | c if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then | |
244 | c write(16,201) INVFLV(lblast(ip)), plast(1,ip), | |
245 | c 1 plast(2,ip),plast(3,ip),plast(4,ip), | |
246 | c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip) | |
247 | c else | |
248 | c if(idpert.eq.1.or.idpert.eq.2) then | |
249 | c write(90,251) INVFLV(lblast(ip)), plast(1,ip), | |
250 | c 1 plast(2,ip),plast(3,ip), | |
251 | c 2 xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip), | |
252 | c 3 dplast(ip) | |
253 | c else | |
254 | c write(99,*) 'Unexpected perturbative particles' | |
255 | c endif | |
256 | c endif | |
257 | c endif | |
258 | c 1007 continue | |
259 | if(ioscar.eq.1) call hoscar | |
260 | endif | |
261 | 190 format(3(i7),f10.4,5x,6(i4)) | |
262 | clin-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) | |
267 | c | |
268 | return | |
269 | end | |
270 | ||
271 | c======================================================================= | |
272 | SUBROUTINE decomp(px0,py0,pz0,xm0,i,itq1) | |
273 | c | |
274 | IMPLICIT DOUBLE PRECISION(D) | |
275 | DOUBLE PRECISION enenew, pxnew, pynew, pznew | |
276 | DOUBLE PRECISION de0, beta2, gam | |
277 | common /lora/ enenew, pxnew, pynew, pznew | |
278 | cc SAVE /lora/ | |
279 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
280 | cc SAVE /HPARNT/ | |
281 | common /decom/ptwo(2,5) | |
282 | cc SAVE /decom/ | |
283 | COMMON/RNDF77/NSEED | |
284 | cc SAVE /RNDF77/ | |
285 | COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 | |
286 | common/embed/iembed,pxqembd,pyqembd,xembd,yembd | |
287 | SAVE | |
288 | c | |
289 | itq1=itq1 | |
290 | dcth=dble(RANART(NSEED))*2.d0-1.d0 | |
291 | dPHI=dble(RANART(NSEED)*HIPR1(40))*2.d0 | |
292 | clin-6/2009 Added if embedding a high-Pt quark pair after string melting: | |
293 | if(iembed.eq.1) then | |
294 | c Decompose the parent high-Pt pion to q and qbar with an internal momentum | |
295 | c parallel to the pion direction so that one parton has ~the same hight Pt | |
296 | c and the other parton has a very soft Pt: | |
297 | c 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 | |
304 | c | |
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) | |
313 | c | |
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 | |
318 | c 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 | |
324 | c | |
325 | call lorenza(de1,dpx,dpy,dpz,-dbex,-dbey,-dbez) | |
326 | ptwo(1,1)=sngl(pxnew) | |
327 | ptwo(1,2)=sngl(pynew) | |
328 | ptwo(1,3)=sngl(pznew) | |
329 | ptwo(1,4)=sngl(enenew) | |
330 | call lorenza(de2,-dpx,-dpy,-dpz,-dbex,-dbey,-dbez) | |
331 | ptwo(2,1)=sngl(pxnew) | |
332 | ptwo(2,2)=sngl(pynew) | |
333 | ptwo(2,3)=sngl(pznew) | |
334 | ptwo(2,4)=sngl(enenew) | |
335 | c | |
336 | RETURN | |
337 | END | |
338 | ||
339 | c======================================================================= | |
340 | SUBROUTINE HTOP | |
341 | c | |
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) | |
350 | cwei DOUBLE PRECISION PATT | |
351 | cc SAVE /HMAIN2/ | |
352 | COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 | |
353 | cc SAVE /HMAIN1/ | |
354 | COMMON /PARA1/ MUL | |
355 | cc 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) | |
359 | cc SAVE /prec1/ | |
360 | COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN) | |
361 | cc 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) | |
366 | cc SAVE /ARPRC/ | |
367 | common /decom/ptwo(2,5) | |
368 | cc SAVE /decom/ | |
369 | COMMON/RNDF77/NSEED | |
370 | cc 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) | |
375 | cc SAVE /NOPREC/ | |
376 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
377 | cc SAVE /HPARNT/ | |
378 | c 7/20/01: use double precision | |
379 | c otherwise sometimes beta>1 and gamma diverge in lorenza(): | |
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) | |
384 | cc SAVE /SOFT/ | |
385 | common/anim/nevent,isoft,isflag,izpc | |
386 | cc SAVE /anim/ | |
387 | DOUBLE PRECISION vxp0,vyp0,vzp0 | |
388 | common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN) | |
389 | cc SAVE /precpa/ | |
390 | common /para7/ ioscar,nsmbbbar,nsmmeson | |
391 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
392 | SAVE | |
393 | c | |
394 | npar=0 | |
395 | nnozpc=0 | |
396 | clin-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 | |
408 | c proj or targ nucleons without interactions, do not enter ZPC: | |
409 | elseif(idabs.gt.1000.and.i2.ne.0) then | |
410 | c 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 | |
414 | c mesons to be converted to q/qbar: | |
415 | nsmmeson=nsmmeson+1 | |
416 | endif | |
417 | enddo | |
418 | ||
419 | clin-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 | |
424 | c write(92,*) iaevt, 3*nsmbbbar+2*nsmmeson | |
425 | c write(92,*) ' event#, total # of initial partons after string | |
426 | c 1 melting' | |
427 | c write(92,*) 'String melting converts ',nsmbbbar, ' baryons &' | |
428 | c 1, nsmmeson, ' mesons' | |
429 | c write(92,*) 'Total # of initial particles= ',natt | |
430 | c write(92,*) 'Total # of initial particles (gamma,e,muon,...) | |
431 | c 1 not entering ZPC= ',natt-nsmbbbar-nsmmeson | |
432 | endif | |
433 | clin-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 | |
448 | c | |
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 | |
451 | c proj or targ nucleons without interactions, do not enter ZPC: | |
452 | inozpc=1 | |
453 | elseif(idabs.gt.1000.and.i2.ne.0) then | |
454 | c 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 | |
506 | c antibaryons: | |
507 | if(id.lt.0) then | |
508 | it(1)=-it(1) | |
509 | it(2)=-it(2) | |
510 | endif | |
511 | c 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 | |
519 | c 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 | |
543 | c save other particles (leptons and photons) outside of ZPC: | |
544 | inozpc=1 | |
545 | endif | |
546 | c | |
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 | |
586 | cyy 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2)) | |
587 | cyy 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2)) | |
588 | c | |
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)) | |
595 | c 8/19/02 avoid actual argument in common blocks of DECOMP: | |
596 | c 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)) | |
601 | c | |
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 | |
621 | c | |
622 | endif | |
623 | 100 continue | |
624 | MUL=NPAR | |
625 | c | |
626 | clin-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 | |
636 | c | |
637 | RETURN | |
638 | END | |
639 | ||
640 | c======================================================================= | |
641 | SUBROUTINE PTOH | |
642 | c | |
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) | |
654 | cc SAVE /loclco/ | |
655 | COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 | |
656 | cc SAVE /HMAIN1/ | |
657 | COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4) | |
658 | cc 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) | |
662 | cc SAVE /HJJET2/ | |
663 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
664 | cc 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) | |
669 | cc 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) | |
674 | cc SAVE /SOFT/ | |
675 | COMMON/RNDF77/NSEED | |
676 | cc SAVE /RNDF77/ | |
677 | common/anim/nevent,isoft,isflag,izpc | |
678 | cc SAVE /anim/ | |
679 | common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom | |
680 | cc SAVE /prtn23/ | |
681 | common /nzpc/nattzp | |
682 | cc SAVE /nzpc/ | |
683 | common /lora/ enenew, pxnew, pynew, pznew | |
684 | cc SAVE /lora/ | |
685 | COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) | |
686 | cc SAVE /LUDAT1A/ | |
687 | clin 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 | |
692 | c | |
693 | dimension xmdiag(MAXSTR),indx(MAXSTR),ndiag(MAXSTR) | |
694 | cwei DOUBLE PRECISION PATT | |
695 | SAVE | |
696 | c | |
697 | call coales | |
698 | c 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. | |
706 | c 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) | |
720 | c 5/02/01 try lowest spin states as first choices, | |
721 | c 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 | |
745 | c***** 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 | |
762 | c | |
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 | |
773 | c | |
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 | |
781 | c 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 | |
786 | c form a decuplet baryon if the q+diquark mass is closer to its mass | |
787 | c (and if the diquark has spin 1): | |
788 | cc for now only include Delta, which is present in ART: | |
789 | cc 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 | |
795 | clin-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 | |
814 | c | |
815 | if(ki.eq.kj.and.ki.eq.kk) then | |
816 | c 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 | |
820 | c form Lambda or Sigma according to 3-quark mass, | |
821 | c 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 | |
831 | cc 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 | |
838 | c***** mesons: | |
839 | else | |
840 | if(k1abs.eq.k2abs) then | |
841 | if(k1abs.le.2) then | |
842 | c treat diagonal mesons later in the subroutine: | |
843 | kf=0 | |
844 | elseif(k1abs.le.3) then | |
845 | c 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 | |
860 | c 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 | |
877 | c 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 | |
882 | c 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 | |
890 | c | |
891 | if(nuudd.gt.1) then | |
892 | call index1(MAXSTR,nuudd,xmdiag,indx) | |
893 | else | |
894 | indx(1)=1 | |
895 | end if | |
896 | c | |
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 | |
905 | c 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 | |
915 | c determine hadron formation time, position and momentum: | |
916 | inatt=0 | |
917 | clin-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) | |
923 | c | |
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) | |
939 | c | |
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) | |
980 | c | |
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 | |
999 | c | |
1000 | call lorenza(ftavg0,gxavg0,gyavg0,gzavg0,-bex,-bey,-bez) | |
1001 | GXAR(inatt)=sngl(pxnew) | |
1002 | GYAR(inatt)=sngl(pynew) | |
1003 | GZAR(inatt)=sngl(pznew) | |
1004 | FTAR(inatt)=sngl(enenew) | |
1005 | clin 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) | |
1015 | c | |
1016 | endif | |
1017 | 1006 CONTINUE | |
1018 | c number of hadrons formed from partons inside ZPC: | |
1019 | nattzp=natt | |
1020 | MSTJ(24)=mstj24 | |
1021 | c | |
1022 | RETURN | |
1023 | END | |
1024 | ||
1025 | c======================================================================= | |
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) | |
1039 | cc SAVE /SOFT/ | |
1040 | common /coal/dpcoal,drcoal,ecritl | |
1041 | cc 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) | |
1044 | cc 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) | |
1048 | cc SAVE /HJJET2/ | |
1049 | SAVE | |
1050 | c | |
1051 | do 1001 ISG=1, NSG | |
1052 | IOVER(ISG)=0 | |
1053 | 1001 continue | |
1054 | C1 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 | |
1057 | C 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 | |
1062 | c | |
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 | |
1076 | c 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))) | |
1079 | c | |
1080 | do 120 JSG=1,NSG | |
1081 | c 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))) | |
1098 | c 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) | |
1110 | c skip if outside of spatial radius: | |
1111 | if(drlocl.gt.drcoal) goto 120 | |
1112 | c 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 | |
1123 | c | |
1124 | C2 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 | |
1127 | C 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))) | |
1143 | c | |
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))) | |
1161 | c 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) | |
1173 | c skip if outside of spatial radius: | |
1174 | if(drlocl.gt.drcoal) goto 220 | |
1175 | c 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 | |
1186 | c | |
1187 | C3 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) | |
1191 | C 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))) | |
1207 | c | |
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))) | |
1221 | c | |
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))) | |
1244 | c 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) | |
1256 | c skip if outside of spatial radius: | |
1257 | if(drlocl.gt.drcoal) goto 320 | |
1258 | c 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 | |
1282 | c | |
1283 | RETURN | |
1284 | END | |
1285 | ||
1286 | c======================================================================= | |
1287 | SUBROUTINE exchge(isg,ipi,jsg,ipj) | |
1288 | c | |
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) | |
1295 | cc SAVE /SOFT/ | |
1296 | SAVE | |
1297 | c | |
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 | |
1331 | c | |
1332 | RETURN | |
1333 | END | |
1334 | ||
1335 | c======================================================================= | |
1336 | SUBROUTINE locldr(icall,drlocl) | |
1337 | c | |
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) | |
1342 | cc SAVE /loclco/ | |
1343 | common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom | |
1344 | cc SAVE /prtn23/ | |
1345 | common /lora/ enenew, pxnew, pynew, pznew | |
1346 | cc SAVE /lora/ | |
1347 | SAVE | |
1348 | c 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 | |
1354 | c 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 | |
1367 | c | |
1368 | call lorenza(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez) | |
1369 | gxp0(j)=pxnew | |
1370 | gyp0(j)=pynew | |
1371 | gzp0(j)=pznew | |
1372 | ftp0(j)=enenew | |
1373 | call lorenza(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez) | |
1374 | pxp0(j)=pxnew | |
1375 | pyp0(j)=pynew | |
1376 | pzp0(j)=pznew | |
1377 | pep0(j)=enenew | |
1378 | 1001 continue | |
1379 | c | |
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) | |
1388 | c | |
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) | |
1396 | c 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 | |
1407 | c boost the reference frame down by beta to get to the 3-parton rest frame: | |
1408 | do 1002 j=1,3 | |
1409 | call lorenza(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez) | |
1410 | gxp0(j)=pxnew | |
1411 | gyp0(j)=pynew | |
1412 | gzp0(j)=pznew | |
1413 | ftp0(j)=enenew | |
1414 | call lorenza(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez) | |
1415 | pxp0(j)=pxnew | |
1416 | pyp0(j)=pynew | |
1417 | pzp0(j)=pznew | |
1418 | pep0(j)=enenew | |
1419 | 1002 continue | |
1420 | c | |
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) | |
1429 | c | |
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 | |
1443 | c | |
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 | |
1451 | c | |
1452 | RETURN | |
1453 | END | |
1454 | ||
1455 | c======================================================================= | |
1456 | subroutine hoscar | |
1457 | c | |
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 | |
1462 | cc SAVE /snn/ | |
1463 | common /lastt/itimeh,bimp | |
1464 | cc SAVE /lastt/ | |
1465 | COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast | |
1466 | cc SAVE /hbt/ | |
1467 | common/oscar1/iap,izp,iat,izt | |
1468 | cc SAVE /oscar1/ | |
1469 | common/oscar2/FRAME,amptvn | |
1470 | cc SAVE /oscar2/ | |
1471 | SAVE | |
1472 | data nff/0/ | |
1473 | c | |
1474 | c 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 | |
1499 | c comment | |
1500 | c event header | |
1501 | write (19, 103) ievent, nlast, bimp, phi | |
1502 | c 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') | |
1519 | c | |
1520 | return | |
1521 | end | |
1522 | ||
1523 | c======================================================================= | |
1524 | subroutine getnp | |
1525 | ||
1526 | PARAMETER (MAXSTR=150001) | |
1527 | COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 | |
1528 | cc SAVE /HMAIN1/ | |
1529 | COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4) | |
1530 | cc SAVE /HMAIN2/ | |
1531 | COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50) | |
1532 | cc SAVE /HPARNT/ | |
1533 | common/snn/efrm,npart1,npart2 | |
1534 | cc SAVE /snn/ | |
1535 | cwei 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 | |
1545 | c | |
1546 | PZPROJ=SQRT(HINT1(6)**2-HINT1(8)**2) | |
1547 | PZTARG=SQRT(HINT1(7)**2-HINT1(9)**2) | |
1548 | epsiln=0.01 | |
1549 | c | |
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 | ||
1568 | c======================================================================= | |
1569 | c 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 | |
1577 | cc SAVE /INPUT2/ | |
1578 | COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5) | |
1579 | cc SAVE /LUJETSA/ | |
1580 | COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) | |
1581 | cc SAVE /LUDAT1A/ | |
1582 | COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) | |
1583 | cc SAVE /LUDAT2A/ | |
1584 | COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5) | |
1585 | cc SAVE /LUDAT3A/ | |
1586 | COMMON /CC/ E(MAXSTR) | |
1587 | cc SAVE /CC/ | |
1588 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
1589 | cc SAVE /EE/ | |
1590 | COMMON /PA/RPION(3,MAXSTR,MAXR) | |
1591 | cc SAVE /PA/ | |
1592 | COMMON /PB/PPION(3,MAXSTR,MAXR) | |
1593 | cc SAVE /PB/ | |
1594 | COMMON /PC/EPION(MAXSTR,MAXR) | |
1595 | cc SAVE /PC/ | |
1596 | COMMON /PD/LPION(MAXSTR,MAXR) | |
1597 | cc SAVE /PD/ | |
1598 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
1599 | cc SAVE /input1/ | |
1600 | common/resdcy/NSAV,iksdcy | |
1601 | cc SAVE /resdcy/ | |
1602 | common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl, | |
1603 | 1 px1n,py1n,pz1n,dp1n | |
1604 | cc SAVE /leadng/ | |
1605 | EXTERNAL IARFLV, INVFLV | |
1606 | COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR) | |
1607 | cc 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) | |
1612 | cc 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 | c & .or.lb1.eq.28.or.lb1.eq.29 | |
1618 | & .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9) | |
1619 | & .or.iabs(lb1).eq.16) then | |
1620 | kf=INVFLV(lb1) | |
1621 | else | |
1622 | return | |
1623 | endif | |
1624 | c | |
1625 | IP=1 | |
1626 | c label as undecayed and the only particle in the record: | |
1627 | N=1 | |
1628 | K(IP,1)=1 | |
1629 | K(IP,3)=0 | |
1630 | K(IP,4)=0 | |
1631 | K(IP,5)=0 | |
1632 | c | |
1633 | K(IP,2)=kf | |
1634 | P(IP,1)=px1 | |
1635 | P(IP,2)=py1 | |
1636 | P(IP,3)=pz1 | |
1637 | em1a=em1 | |
1638 | c eta or omega in ART may be below or too close to (pi+pi-pi0) mass, | |
1639 | c causing LUDECY error,thus increase their mass ADDM above this thresh, | |
1640 | c noting that rho (m=0.281) too close to 2pi thrshold fails to decay: | |
1641 | if((lb1.eq.0.or.lb1.eq.28).and.em1.lt.(2*APICH+API0+ADDM)) then | |
1642 | em1=2*APICH+API0+ADDM | |
1643 | c rho | |
1644 | elseif(lb1.ge.25.and.lb1.le.27.and.em1.lt.(2*APICH+ADDM)) then | |
1645 | em1=2*APICH+ADDM | |
1646 | c K* | |
1647 | elseif(iabs(lb1).eq.30.and.em1.lt.(APICH+AK0+ADDM)) then | |
1648 | em1=APICH+AK0+ADDM | |
1649 | c Delta created in ART may be below (n+pich) mass, causing LUDECY error: | |
1650 | elseif(iabs(lb1).ge.6.and.iabs(lb1).le.9 | |
1651 | 1 .and.em1.lt.(APICH+AN+ADDM)) then | |
1652 | em1=APICH+AN+ADDM | |
1653 | endif | |
1654 | c if(em1.ge.(em1a+0.01)) write (6,*) | |
1655 | c 1 'Mass increase in resdec():',nt,em1-em1a,lb1 | |
1656 | e1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) | |
1657 | P(IP,4)=e1 | |
1658 | P(IP,5)=em1 | |
1659 | clin-5/2008: | |
1660 | dpdecp=dpertp(i1) | |
1661 | call ludecy(IP) | |
1662 | c add decay time to daughter's formation time at the last timestep: | |
1663 | if(nt.eq.ntmax) then | |
1664 | tau0=hbarc/wid | |
1665 | taudcy=tau0*(-1.)*alog(1.-RANART(NSEED)) | |
1666 | ndaut=n-nsav | |
1667 | if(ndaut.le.1) then | |
1668 | write(10,*) 'note: ndaut(<1)=',ndaut | |
1669 | write(89,*) 'note: ndaut(<1)=',ndaut | |
1670 | call lulist(2) | |
1671 | stop | |
1672 | endif | |
1673 | c lorentz boost: | |
1674 | taudcy=taudcy*e1/em1 | |
1675 | tfnl=tfnl+taudcy | |
1676 | xfnl=xfnl+px1/e1*taudcy | |
1677 | yfnl=yfnl+py1/e1*taudcy | |
1678 | zfnl=zfnl+pz1/e1*taudcy | |
1679 | c at the last timestep, assign rho, K0S or eta (decay daughter) | |
1680 | c to lb(i1) only (not to lpion) in order to decay them again: | |
1681 | if(n.ge.(nsav+2)) then | |
1682 | do 1001 idau=nsav+2,n | |
1683 | kdaut=K(idau,2) | |
1684 | if(kdaut.eq.221.or.kdaut.eq.113 | |
1685 | 1 .or.kdaut.eq.213.or.kdaut.eq.-213 | |
1686 | 2 .or.kdaut.eq.310) then | |
1687 | c switch idau and i1(nsav+1): | |
1688 | ksave=kdaut | |
1689 | pxsave=p(idau,1) | |
1690 | pysave=p(idau,2) | |
1691 | pzsave=p(idau,3) | |
1692 | esave=p(idau,4) | |
1693 | xmsave=p(idau,5) | |
1694 | K(idau,2)=K(nsav+1,2) | |
1695 | p(idau,1)=p(nsav+1,1) | |
1696 | p(idau,2)=p(nsav+1,2) | |
1697 | p(idau,3)=p(nsav+1,3) | |
1698 | p(idau,4)=p(nsav+1,4) | |
1699 | p(idau,5)=p(nsav+1,5) | |
1700 | K(nsav+1,2)=ksave | |
1701 | p(nsav+1,1)=pxsave | |
1702 | p(nsav+1,2)=pysave | |
1703 | p(nsav+1,3)=pzsave | |
1704 | p(nsav+1,4)=esave | |
1705 | p(nsav+1,5)=xmsave | |
1706 | c note: phi decay may produce rho, K0s or eta, N*(1535) decay may produce | |
1707 | c eta, but only one daughter may be rho, K0s or eta: | |
1708 | goto 111 | |
1709 | endif | |
1710 | 1001 continue | |
1711 | endif | |
1712 | 111 continue | |
1713 | c | |
1714 | enet=0. | |
1715 | do 1002 idau=nsav+1,n | |
1716 | enet=enet+p(idau,4) | |
1717 | 1002 continue | |
1718 | c if(abs(enet-e1).gt.0.02) | |
1719 | c 1 write(93,*) 'resdec(): nt=',nt,enet-e1,lb1 | |
1720 | endif | |
1721 | ||
1722 | cyy 200 format(a20,3(1x,i6)) | |
1723 | cyy 210 format(i6,5(1x,f8.3)) | |
1724 | cyy 220 format(a2,i5,5(1x,f8.3)) | |
1725 | ||
1726 | do 1003 idau=nsav+1,n | |
1727 | kdaut=K(idau,2) | |
1728 | lbdaut=IARFLV(kdaut) | |
1729 | c K0S and K0L are named K+/K- during hadron cascade, and only | |
1730 | c at the last timestep they keep their real LB # before output; | |
1731 | c K0/K0bar (from K* decay) converted to K0S and K0L at the last timestep: | |
1732 | if(nt.eq.ntmax.and.(kdaut.eq.130.or.kdaut.eq.310 | |
1733 | 1 .or.iabs(kdaut).eq.311)) then | |
1734 | if(kdaut.eq.130) then | |
1735 | lbdaut=22 | |
1736 | elseif(kdaut.eq.310) then | |
1737 | lbdaut=24 | |
1738 | elseif(iabs(kdaut).eq.311) then | |
1739 | if(RANART(NSEED).lt.0.5) then | |
1740 | lbdaut=22 | |
1741 | else | |
1742 | lbdaut=24 | |
1743 | endif | |
1744 | endif | |
1745 | endif | |
1746 | c | |
1747 | if(idau.eq.(nsav+1)) then | |
1748 | LB(i1)=lbdaut | |
1749 | E(i1)=p(idau,5) | |
1750 | px1n=p(idau,1) | |
1751 | py1n=p(idau,2) | |
1752 | pz1n=p(idau,3) | |
1753 | clin-5/2008: | |
1754 | dp1n=dpdecp | |
1755 | else | |
1756 | nnn=nnn+1 | |
1757 | LPION(NNN,IRUN)=lbdaut | |
1758 | EPION(NNN,IRUN)=p(idau,5) | |
1759 | PPION(1,NNN,IRUN)=p(idau,1) | |
1760 | PPION(2,NNN,IRUN)=p(idau,2) | |
1761 | PPION(3,NNN,IRUN)=p(idau,3) | |
1762 | RPION(1,NNN,IRUN)=xfnl | |
1763 | RPION(2,NNN,IRUN)=yfnl | |
1764 | RPION(3,NNN,IRUN)=zfnl | |
1765 | tfdpi(NNN,IRUN)=tfnl | |
1766 | clin-5/2008: | |
1767 | dppion(NNN,IRUN)=dpdecp | |
1768 | endif | |
1769 | 1003 continue | |
1770 | cyy 230 format(a2,i5,5(1x,e8.2)) | |
1771 | return | |
1772 | end | |
1773 | ||
1774 | c======================================================================= | |
1775 | subroutine inidcy | |
1776 | ||
1777 | COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5) | |
1778 | cc SAVE /LUJETSA/ | |
1779 | common/resdcy/NSAV,iksdcy | |
1780 | cc SAVE /resdcy/ | |
1781 | SAVE | |
1782 | N=1 | |
1783 | NSAV=N | |
1784 | return | |
1785 | end | |
1786 | ||
1787 | c======================================================================= | |
1788 | clin-6/06/02 local parton freezeout motivated from critical density: | |
1789 | subroutine local(t) | |
1790 | c | |
1791 | implicit double precision (a-h, o-z) | |
1792 | PARAMETER (MAXPTN=400001) | |
1793 | PARAMETER (r0=1d0) | |
1794 | COMMON /para1/ mul | |
1795 | cc SAVE /para1/ | |
1796 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
1797 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
1798 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
1799 | cc SAVE /prec2/ | |
1800 | common /frzprc/ | |
1801 | & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN), | |
1802 | & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN), | |
1803 | & xmfrz(MAXPTN), | |
1804 | & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast | |
1805 | cc SAVE /frzprc/ | |
1806 | common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN) | |
1807 | cc SAVE /prec4/ | |
1808 | common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN) | |
1809 | cc SAVE /prec5/ | |
1810 | common /coal/dpcoal,drcoal,ecritl | |
1811 | cc SAVE /coal/ | |
1812 | SAVE | |
1813 | c | |
1814 | do 1001 it=1,301 | |
1815 | if(t.ge.tfrz(it).and.t.lt.tfrz(it+1)) then | |
1816 | if(it.eq.itlast) then | |
1817 | return | |
1818 | else | |
1819 | itlast=it | |
1820 | goto 50 | |
1821 | endif | |
1822 | endif | |
1823 | 1001 continue | |
1824 | write(1,*) 'local time out of range in LOCAL, stop',t,it | |
1825 | stop | |
1826 | 50 continue | |
1827 | c | |
1828 | do 200 ip=1,mul | |
1829 | c skip partons which have frozen out: | |
1830 | if(ifrz(ip).eq.1) goto 200 | |
1831 | if(it.eq.301) then | |
1832 | c freezeout all the left partons beyond the time of 3000 fm: | |
1833 | etcrit=1d6 | |
1834 | goto 150 | |
1835 | else | |
1836 | c freezeout when transverse energy density < etcrit: | |
1837 | etcrit=(ecritl*2d0/3d0) | |
1838 | endif | |
1839 | c skip partons which have not yet formed: | |
1840 | if(t.lt.FT5(ip)) goto 200 | |
1841 | rap0=rap(ip) | |
1842 | eta0=eta(ip) | |
1843 | x0=GX5(ip)+vx(ip)*(t-FT5(ip)) | |
1844 | y0=GY5(ip)+vy(ip)*(t-FT5(ip)) | |
1845 | detdy=0d0 | |
1846 | do 100 itest=1,mul | |
1847 | c skip self and partons which have not yet formed: | |
1848 | if(itest.eq.ip.or.t.lt.FT5(itest)) goto 100 | |
1849 | ettest=eta(itest) | |
1850 | xtest=GX5(itest)+vx(itest)*(t-FT5(itest)) | |
1851 | ytest=GY5(itest)+vy(itest)*(t-FT5(itest)) | |
1852 | drt=sqrt((xtest-x0)**2+(ytest-y0)**2) | |
1853 | c count partons within drt<1 and -1<(eta-eta0)<1: | |
1854 | if(dabs(ettest-eta0).le.1d0.and.drt.le.r0) | |
1855 | 1 detdy=detdy+dsqrt(PX5(itest)**2+PY5(itest)**2 | |
1856 | 2 +XMASS5(itest)**2)*0.5d0 | |
1857 | 100 continue | |
1858 | detdy=detdy*(dcosh(eta0)**2)/(t*3.1416d0*r0**2*dcosh(rap0)) | |
1859 | c when density is below critical density for phase transition, freeze out: | |
1860 | 150 if(detdy.le.etcrit) then | |
1861 | ifrz(ip)=1 | |
1862 | idfrz(ip)=ITYP5(ip) | |
1863 | pxfrz(ip)=PX5(ip) | |
1864 | pyfrz(ip)=PY5(ip) | |
1865 | pzfrz(ip)=PZ5(ip) | |
1866 | efrz(ip)=E5(ip) | |
1867 | xmfrz(ip)=XMASS5(ip) | |
1868 | if(t.gt.FT5(ip)) then | |
1869 | gxfrz(ip)=x0 | |
1870 | gyfrz(ip)=y0 | |
1871 | gzfrz(ip)=GZ5(ip)+vz(ip)*(t-FT5(ip)) | |
1872 | ftfrz(ip)=t | |
1873 | else | |
1874 | c if this freezeout time < formation time, use formation time & positions. | |
1875 | c This ensures the recovery of default hadron when e_crit=infty: | |
1876 | gxfrz(ip)=GX5(ip) | |
1877 | gyfrz(ip)=GY5(ip) | |
1878 | gzfrz(ip)=GZ5(ip) | |
1879 | ftfrz(ip)=FT5(ip) | |
1880 | endif | |
1881 | endif | |
1882 | 200 continue | |
1883 | c | |
1884 | return | |
1885 | end | |
1886 | ||
1887 | c======================================================================= | |
1888 | clin-6/06/02 initialization for local parton freezeout | |
1889 | subroutine inifrz | |
1890 | c | |
1891 | implicit double precision (a-h, o-z) | |
1892 | PARAMETER (MAXPTN=400001) | |
1893 | common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge | |
1894 | cc SAVE /ilist5/ | |
1895 | common /frzprc/ | |
1896 | & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN), | |
1897 | & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN), | |
1898 | & xmfrz(MAXPTN), | |
1899 | & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast | |
1900 | cc SAVE /frzprc/ | |
1901 | SAVE | |
1902 | c | |
1903 | c for freezeout time 0-10fm, use interval of 0.1fm; | |
1904 | c for 10-100fm, use interval of 1fm; | |
1905 | c for 100-1000fm, use interval of 10fm; | |
1906 | c for 1000-3000fm, use interval of 100fm: | |
1907 | step1=0.1d0 | |
1908 | step2=1d0 | |
1909 | step3=10d0 | |
1910 | step4=100d0 | |
1911 | c | |
1912 | do 1001 it=1,101 | |
1913 | tfrz(it)=0d0+dble(it-1)*step1 | |
1914 | 1001 continue | |
1915 | do 1002 it=102,191 | |
1916 | tfrz(it)=10d0+dble(it-101)*step2 | |
1917 | 1002 continue | |
1918 | do 1003 it=192,281 | |
1919 | tfrz(it)=100d0+dble(it-191)*step3 | |
1920 | 1003 continue | |
1921 | do 1004 it=282,301 | |
1922 | tfrz(it)=1000d0+dble(it-281)*step4 | |
1923 | 1004 continue | |
1924 | tfrz(302)=tlarge | |
1925 | c | |
1926 | return | |
1927 | end | |
1928 | ||
1929 | clin-5/2009 ctest on v2 analysis | |
1930 | c======================================================================= | |
1931 | c idd=0,1,2,3 specifies different subroutines for partonic flow analysis. | |
1932 | subroutine flowp(idd) | |
1933 | c | |
1934 | implicit double precision (a-h, o-z) | |
1935 | real dt | |
1936 | parameter (MAXPTN=400001) | |
1937 | csp | |
1938 | parameter (bmt=0.05d0) | |
1939 | dimension nlfile(3),nsfile(3),nmfile(3) | |
1940 | c | |
1941 | dimension v2pp(3),xnpp(3),v2psum(3),v2p2sm(3),nfile(3) | |
1942 | dimension tsp(31),v2pevt(3),v2pavg(3),varv2p(3) | |
1943 | common /ilist1/ | |
1944 | & iscat, jscat, next(MAXPTN), last(MAXPTN), | |
1945 | & ictype, icsta(MAXPTN), | |
1946 | & nic(MAXPTN), icels(MAXPTN) | |
1947 | cc SAVE /ilist1/ | |
1948 | COMMON /para1/ mul | |
1949 | cc SAVE /para1/ | |
1950 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
1951 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
1952 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
1953 | cc SAVE /prec2/ | |
1954 | COMMON /pflow/ v2p(30,3),xnpart(30,3),etp(30,3), | |
1955 | 1 s2p(30,3),v2p2(30,3),nevt(30) | |
1956 | cc SAVE /pflow/ | |
1957 | COMMON /pflowf/ v2pf(30,3),xnpf(30,3),etpf(30,3), | |
1958 | 1 xncoll(30),s2pf(30,3),v2pf2(30,3) | |
1959 | cc SAVE /pflowf/ | |
1960 | COMMON /pfrz/ v2pfrz(30,3),xnpfrz(30,3),etpfrz(30,3), | |
1961 | 1 s2pfrz(30,3),v2p2fz(30,3),tscatt(31), | |
1962 | 2 nevtfz(30),iscatt(30) | |
1963 | cc SAVE /pfrz/ | |
1964 | COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3), | |
1965 | 1 v2h2(30,3),s2h(30,3) | |
1966 | cc SAVE /hflow/ | |
1967 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
1968 | cc SAVE /AREVT/ | |
1969 | common/anim/nevent,isoft,isflag,izpc | |
1970 | cc SAVE /anim/ | |
1971 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
1972 | cc SAVE /input1/ | |
1973 | COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, | |
1974 | & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL | |
1975 | cc SAVE /INPUT2/ | |
1976 | cc SAVE itimep,iaevtp,v2pp,xnpp,v2psum,v2p2sm | |
1977 | cc SAVE nfile,itanim,nlfile,nsfile,nmfile | |
1978 | SAVE | |
1979 | csp | |
1980 | dimension etpl(30,3),etps(30,3),etplf(30,3),etpsf(30,3), | |
1981 | & etlfrz(30,3),etsfrz(30,3), | |
1982 | & xnpl(30,3),xnps(30,3),xnplf(30,3),xnpsf(30,3), | |
1983 | & xnlfrz(30,3),xnsfrz(30,3), | |
1984 | & v2pl(30,3),v2ps(30,3),v2plf(30,3),v2psf(30,3), | |
1985 | & s2pl(30,3),s2ps(30,3),s2plf(30,3),s2psf(30,3), | |
1986 | & DMYil(50,3),DMYfl(50,3), | |
1987 | & DMYis(50,3),DMYfs(50,3) | |
1988 | data tsp/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, | |
1989 | & 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, | |
1990 | & 2 , 3, 4, 5, 6, 7, 8, 9, 10, 20, 30/ | |
1991 | c idd=0: initialization for flow analysis, called by artdri.f: | |
1992 | if(idd.eq.0) then | |
1993 | nfile(1)=60 | |
1994 | nfile(2)=64 | |
1995 | nfile(3)=20 | |
1996 | cms OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN') | |
1997 | cms OPEN (nfile(1)+1, | |
1998 | cms 1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN') | |
1999 | cms OPEN (nfile(1)+2, | |
2000 | cms 1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN') | |
2001 | cms OPEN (nfile(1)+3, | |
2002 | cms 1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN') | |
2003 | cms OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN') | |
2004 | cms OPEN (nfile(2)+1, | |
2005 | cms 1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN') | |
2006 | cms OPEN (nfile(2)+2, | |
2007 | cms 1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN') | |
2008 | cms OPEN (nfile(2)+3, | |
2009 | cms 1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN') | |
2010 | cms OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN') | |
2011 | cms OPEN (nfile(3)+1, | |
2012 | cms 1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN') | |
2013 | cms OPEN (nfile(3)+2, | |
2014 | cms 1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN') | |
2015 | cms OPEN (nfile(3)+3, | |
2016 | cms 1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN') | |
2017 | cms OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN') | |
2018 | cms write(49, *) ' ievt, v2p, v2p_y2, v2p_y1' | |
2019 | c | |
2020 | cms OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN') | |
2021 | cms OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN') | |
2022 | cms OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN') | |
2023 | cms OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN') | |
2024 | cms write(88, *) ' ievt, v2h, v2h_y2, v2h_y1' | |
2025 | csp07/05 | |
2026 | nlfile(1)=70 | |
2027 | nlfile(2)=72 | |
2028 | nlfile(3)=74 | |
2029 | cms OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN') | |
2030 | cms OPEN (nlfile(1)+1, | |
2031 | cms 1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN') | |
2032 | cms OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN') | |
2033 | cms OPEN (nlfile(2)+1, | |
2034 | cms 1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN') | |
2035 | cms OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN') | |
2036 | cms OPEN (nlfile(3)+1, | |
2037 | cms 1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN') | |
2038 | nsfile(1)=76 | |
2039 | nsfile(2)=78 | |
2040 | nsfile(3)=80 | |
2041 | cms OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN') | |
2042 | cms OPEN (nsfile(1)+1, | |
2043 | cms 1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN') | |
2044 | cms OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN') | |
2045 | cms OPEN (nsfile(2)+1, | |
2046 | cms 1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN') | |
2047 | cms OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN') | |
2048 | cms OPEN (nsfile(3)+1, | |
2049 | cms 1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN') | |
2050 | nmfile(1)=82 | |
2051 | nmfile(2)=83 | |
2052 | nmfile(3)=84 | |
2053 | cms OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN') | |
2054 | cms OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN') | |
2055 | cms OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN') | |
2056 | clin-11/27/00 for animation: | |
2057 | if(nevent.eq.1) then | |
2058 | cms OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN') | |
2059 | cms write(91,*) ntmax, dt | |
2060 | cms OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN') | |
2061 | cms OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN') | |
2062 | endif | |
2063 | c | |
2064 | itimep=0 | |
2065 | itanim=0 | |
2066 | iaevtp=0 | |
2067 | csp | |
2068 | do 1002 ii=1,50 | |
2069 | do 1001 iy=1,3 | |
2070 | DMYil(ii,iy) = 0d0 | |
2071 | DMYfl(ii,iy) = 0d0 | |
2072 | DMYis(ii,iy) = 0d0 | |
2073 | DMYfs(ii,iy) = 0d0 | |
2074 | 1001 continue | |
2075 | 1002 continue | |
2076 | c | |
2077 | do 1003 ii=1,31 | |
2078 | tscatt(ii)=0d0 | |
2079 | 1003 continue | |
2080 | do 1005 ii=1,30 | |
2081 | nevt(ii)=0 | |
2082 | xncoll(ii)=0d0 | |
2083 | nevtfz(ii)=0 | |
2084 | iscatt(ii)=0 | |
2085 | do 1004 iy=1,3 | |
2086 | v2p(ii,iy)=0d0 | |
2087 | v2p2(ii,iy)=0d0 | |
2088 | s2p(ii,iy)=0d0 | |
2089 | etp(ii,iy)=0d0 | |
2090 | xnpart(ii,iy)=0d0 | |
2091 | v2pf(ii,iy)=0d0 | |
2092 | v2pf2(ii,iy)=0d0 | |
2093 | s2pf(ii,iy)=0d0 | |
2094 | etpf(ii,iy)=0d0 | |
2095 | xnpf(ii,iy)=0d0 | |
2096 | v2pfrz(ii,iy)=0d0 | |
2097 | v2p2fz(ii,iy)=0d0 | |
2098 | s2pfrz(ii,iy)=0d0 | |
2099 | etpfrz(ii,iy)=0d0 | |
2100 | xnpfrz(ii,iy)=0d0 | |
2101 | csp07/05 | |
2102 | etpl(ii,iy)=0d0 | |
2103 | etps(ii,iy)=0d0 | |
2104 | etplf(ii,iy)=0d0 | |
2105 | etpsf(ii,iy)=0d0 | |
2106 | etlfrz(ii,iy)=0d0 | |
2107 | etsfrz(ii,iy)=0d0 | |
2108 | xnpl(ii,iy)=0d0 | |
2109 | xnps(ii,iy)=0d0 | |
2110 | xnplf(ii,iy)=0d0 | |
2111 | xnpsf(ii,iy)=0d0 | |
2112 | xnlfrz(ii,iy)=0d0 | |
2113 | xnsfrz(ii,iy)=0d0 | |
2114 | v2pl(ii,iy)=0d0 | |
2115 | v2ps(ii,iy)=0d0 | |
2116 | v2plf(ii,iy)=0d0 | |
2117 | v2psf(ii,iy)=0d0 | |
2118 | s2pl(ii,iy)=0d0 | |
2119 | s2ps(ii,iy)=0d0 | |
2120 | s2plf(ii,iy)=0d0 | |
2121 | s2psf(ii,iy)=0d0 | |
2122 | 1004 continue | |
2123 | 1005 continue | |
2124 | do 1006 iy=1,3 | |
2125 | v2pevt(iy)=0d0 | |
2126 | v2pavg(iy)=0d0 | |
2127 | varv2p(iy)=0d0 | |
2128 | v2pp(iy)=0.d0 | |
2129 | xnpp(iy)=0d0 | |
2130 | v2psum(iy)=0.d0 | |
2131 | v2p2sm(iy)=0.d0 | |
2132 | 1006 continue | |
2133 | c idd=1: calculate parton elliptic flow, called by zpc.f: | |
2134 | else if(idd.eq.1) then | |
2135 | t2time = FT5(iscat) | |
2136 | do 1008 ianp = 1, 30 | |
2137 | if (t2time.lt.tsp(ianp+1).and.t2time.ge.tsp(ianp)) then | |
2138 | c write flow info only once at each fixed time: | |
2139 | xncoll(ianp)=xncoll(ianp)+1d0 | |
2140 | c to prevent an earlier t2time comes later in the same event | |
2141 | c and mess up nevt: | |
2142 | if(ianp.le.itimep.and.iaevt.eq.iaevtp) goto 101 | |
2143 | nevt(ianp)=nevt(ianp)+1 | |
2144 | tscatt(ianp+1)=t2time | |
2145 | iscatt(ianp)=1 | |
2146 | nevtfz(ianp)=nevtfz(ianp)+1 | |
2147 | do 100 I=1,mul | |
2148 | ypartn=0.5d0*dlog((E5(i)+PZ5(i))/(E5(i)-PZ5(i)+1.d-8)) | |
2149 | pt2=PX5(I)**2+PY5(I)**2 | |
2150 | ctest off: initial (pt,y) and (x,y) distribution: | |
2151 | c idtime=1 | |
2152 | c if(ianp.eq.idtime) then | |
2153 | c iityp=iabs(ITYP5(I)) | |
2154 | c if(iityp.eq.1.or.iityp.eq.2) then | |
2155 | c write(651,*) dsqrt(pt2),ypartn | |
2156 | c write(654,*) GX5(I),GY5(I) | |
2157 | c elseif(iityp.eq.1103.or.iityp.eq.2101 | |
2158 | c 1 .or.iityp.eq.2103.or.iityp.eq.2203. | |
2159 | c 2 .or.iityp.eq.3101.or.iityp.eq.3103. | |
2160 | c 3 .or.iityp.eq.3201.or.iityp.eq.3203.or.iityp.eq.3303) | |
2161 | c 4 then | |
2162 | c write(652,*) dsqrt(pt2),ypartn | |
2163 | c write(655,*) GX5(I),GY5(I) | |
2164 | c elseif(iityp.eq.21) then | |
2165 | c write(653,*) dsqrt(pt2),ypartn | |
2166 | c write(656,*) GX5(I),GY5(I) | |
2167 | c endif | |
2168 | c endif | |
2169 | ctest-end | |
2170 | ctest off density with 2fm radius and z:(-0.1*t,0.1*t): | |
2171 | c gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i) | |
2172 | c gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i) | |
2173 | c gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i) | |
2174 | c rt_now=dsqrt(gx_now**2+gy_now**2) | |
2175 | c zmax=0.1d0*t2time | |
2176 | c volume=3.1416d0*(2d0**2)*(2*zmax) | |
2177 | c if(rt_now.gt.2d0.or.dabs(gz_now).gt.zmax) | |
2178 | c 1 goto 100 | |
2179 | ctest-end | |
2180 | iloop=1 | |
2181 | if(dabs(ypartn).le.1d0) then | |
2182 | iloop=2 | |
2183 | if(dabs(ypartn).le.0.5d0) then | |
2184 | iloop=3 | |
2185 | endif | |
2186 | endif | |
2187 | do 50 iy=1,iloop | |
2188 | if(pt2.gt.0.) then | |
2189 | v2prtn=(PX5(I)**2-PY5(I)**2)/pt2 | |
2190 | if(abs(v2prtn).gt.1.) | |
2191 | 1 write(nfile(iy),*) 'v2prtn>1',v2prtn | |
2192 | v2p(ianp,iy)=v2p(ianp,iy)+v2prtn | |
2193 | v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2 | |
2194 | endif | |
2195 | xperp2=GX5(I)**2+GY5(I)**2 | |
2196 | if(xperp2.gt.0.) | |
2197 | 1 s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2 | |
2198 | xnpart(ianp,iy)=xnpart(ianp,iy)+1d0 | |
2199 | etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) | |
2200 | ctest off density: | |
2201 | c etp(ianp,iy)=etp(ianp,iy) | |
2202 | c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume | |
2203 | clin-2/22/00 to write out parton info only for formed ones: | |
2204 | if(FT5(I).le.t2time) then | |
2205 | v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn | |
2206 | v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2 | |
2207 | if(xperp2.gt.0.) | |
2208 | 1 s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2 | |
2209 | xnpf(ianp,iy)=xnpf(ianp,iy)+1d0 | |
2210 | etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) | |
2211 | ctest off density: | |
2212 | c etpf(ianp,iy)=etpf(ianp,iy) | |
2213 | c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume | |
2214 | endif | |
2215 | 50 continue | |
2216 | 100 continue | |
2217 | itimep=ianp | |
2218 | iaevtp=iaevt | |
2219 | clin-3/30/00 ebe v2 variables: | |
2220 | if(ianp.eq.30) then | |
2221 | do 1007 iy=1,3 | |
2222 | npartn=IDINT(xnpart(ianp,iy)-xnpp(iy)) | |
2223 | if(npartn.ne.0) then | |
2224 | v2pevt(iy)=(v2p(ianp,iy)-v2pp(iy))/npartn | |
2225 | v2psum(iy)=v2psum(iy)+v2pevt(iy) | |
2226 | v2p2sm(iy)=v2p2sm(iy)+v2pevt(iy)**2 | |
2227 | v2pp(iy)=v2p(ianp,iy) | |
2228 | xnpp(iy)=xnpart(ianp,iy) | |
2229 | endif | |
2230 | 1007 continue | |
2231 | write(49, 160) iaevt,v2pevt | |
2232 | endif | |
2233 | goto 101 | |
2234 | endif | |
2235 | 1008 continue | |
2236 | clin-11/28/00 for animation: | |
2237 | 101 if(nevent.eq.1) then | |
2238 | do 110 nt = 1, ntmax | |
2239 | time1=dble(nt*dt) | |
2240 | time2=dble((nt+1)*dt) | |
2241 | if (t2time.lt.time2.and.t2time.ge.time1) then | |
2242 | if(nt.le.itanim) return | |
2243 | c write(92,*) t2time | |
2244 | iform=0 | |
2245 | do 1009 I=1,mul | |
2246 | c write out parton info only for formed ones: | |
2247 | if(FT5(I).le.t2time) then | |
2248 | iform=iform+1 | |
2249 | endif | |
2250 | 1009 continue | |
2251 | c write(92,*) iform | |
2252 | do 120 I=1,mul | |
2253 | if(FT5(I).le.t2time) then | |
2254 | clin-11/29/00-ctest off calculate parton coordinates after propagation: | |
2255 | c gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i) | |
2256 | c gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i) | |
2257 | c gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i) | |
2258 | c write(92,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i) | |
2259 | c write(92,180) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i), | |
2260 | c 1 PX5(i),PY5(i),PZ5(i),E5(i) | |
2261 | ctest-end | |
2262 | endif | |
2263 | 120 continue | |
2264 | itanim=nt | |
2265 | endif | |
2266 | 110 continue | |
2267 | endif | |
2268 | c | |
2269 | 140 format(i10,4(2x,f7.2)) | |
2270 | 160 format(i10,3(2x,f9.5)) | |
2271 | c 180 format(i6,8(1x,f7.2)) | |
2272 | clin-5/17/01 calculate v2 for active partons (which have not frozen out): | |
2273 | c idd=3, called at end of zpc.f: | |
2274 | else if(idd.eq.3) then | |
2275 | do 1010 ianp=1,30 | |
2276 | if(iscatt(ianp).eq.0) tscatt(ianp+1)=tscatt(ianp) | |
2277 | 1010 continue | |
2278 | do 350 I=1,mul | |
2279 | ypartn=0.5d0*dlog((E5(i)+PZ5(i)+1.d-8) | |
2280 | 1 /(E5(i)-PZ5(i)+1.d-8)) | |
2281 | pt2=PX5(I)**2+PY5(I)**2 | |
2282 | iloop=1 | |
2283 | if(dabs(ypartn).le.1d0) then | |
2284 | iloop=2 | |
2285 | if(dabs(ypartn).le.0.5d0) then | |
2286 | iloop=3 | |
2287 | endif | |
2288 | endif | |
2289 | c | |
2290 | do 325 ianp=1,30 | |
2291 | if(iscatt(ianp).ne.0) then | |
2292 | if(FT5(I).lt.tscatt(ianp+1) | |
2293 | 1 .and.FT5(I).ge.tscatt(ianp)) then | |
2294 | do 1011 iy=1,iloop | |
2295 | if(pt2.gt.0.) then | |
2296 | v2prtn=(PX5(I)**2-PY5(I)**2)/pt2 | |
2297 | v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn | |
2298 | v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2 | |
2299 | endif | |
2300 | xperp2=GX5(I)**2+GY5(I)**2 | |
2301 | if(xperp2.gt.0.) s2pfrz(ianp,iy)= | |
2302 | 1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2 | |
2303 | etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) | |
2304 | xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0 | |
2305 | ctest off density: | |
2306 | c etpfrz(ianp,iy)=etpfrz(ianp,iy) | |
2307 | c 1 +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume | |
2308 | csp07/05 | |
2309 | if(ITYP5(I).eq.1.or.ITYP5(I).eq.2)then | |
2310 | etlfrz(ianp,iy)=etlfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) | |
2311 | xnlfrz(ianp,iy)=xnlfrz(ianp,iy)+1d0 | |
2312 | elseif(ITYP5(I).eq.3)then | |
2313 | etsfrz(ianp,iy)=etsfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2) | |
2314 | xnsfrz(ianp,iy)=xnsfrz(ianp,iy)+1d0 | |
2315 | endif | |
2316 | csp07/05 end | |
2317 | 1011 continue | |
2318 | c parton freezeout info taken, proceed to next parton: | |
2319 | goto 350 | |
2320 | endif | |
2321 | endif | |
2322 | 325 continue | |
2323 | 350 continue | |
2324 | c idd=2: calculate average partonic elliptic flow, called from artdri.f, | |
2325 | else if(idd.eq.2) then | |
2326 | do 1012 i=1,3 | |
2327 | cms write(nfile(i),*) ' tsp, v2p, v2p2, '// | |
2328 | cms 1 ' s2p, etp, xmult, nevt, xnctot' | |
2329 | cms write ((nfile(i)+1),*) ' tsp, v2pf, v2pf2, '// | |
2330 | cms 1 ' s2pf, etpf, xnform, xcrate' | |
2331 | cms write ((nfile(i)+2),*) ' tsp, v2pa, v2pa2, '// | |
2332 | cms 1 ' s2pa, etpa, xmult_ap, xnform, nevt' | |
2333 | cms write ((nfile(i)+3),*) ' tsph, v2ph, v2ph2, '// | |
2334 | cms 1 ' s2ph, etph, xmult_(ap/2+h),xmult_ap/2,nevt' | |
2335 | csp | |
2336 | cms write(nlfile(i),*) ' tsp, v2, s2, etp, xmul' | |
2337 | cms write(nsfile(i),*) ' tsp, v2, s2, etp, xmul' | |
2338 | cms write(nlfile(i)+1,*) ' tsp, v2, s2, etp, xmul' | |
2339 | cms write(nsfile(i)+1,*) ' tsp, v2, s2, etp, xmul' | |
2340 | c | |
2341 | 1012 continue | |
2342 | clin-3/30/00 ensemble average & variance of v2 (over particles in all events): | |
2343 | do 150 ii=1, 30 | |
2344 | if(nevt(ii).eq.0) goto 150 | |
2345 | do 1014 iy=1,3 | |
2346 | if(xnpart(ii,iy).gt.1) then | |
2347 | v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy) | |
2348 | v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy) | |
2349 | 1 -v2p(ii,iy)**2)/(xnpart(ii,iy)-1)) | |
2350 | s2p(ii,iy)=s2p(ii,iy)/xnpart(ii,iy) | |
2351 | c xmult and etp are multiplicity and et for an averaged event: | |
2352 | xmult=dble(xnpart(ii,iy)/dble(nevt(ii))) | |
2353 | etp(ii,iy)=etp(ii,iy)/dble(nevt(ii)) | |
2354 | csp | |
2355 | etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii)) | |
2356 | etps(ii,iy)=etps(ii,iy)/dble(nevt(ii)) | |
2357 | c | |
2358 | xnctot=0d0 | |
2359 | do 1013 inum=1,ii | |
2360 | xnctot=xnctot+xncoll(inum) | |
2361 | 1013 continue | |
2362 | if(nevt(1).ne.0) xnctot=xnctot/nevt(1) | |
2363 | cms write (nfile(iy),200) tsp(ii),v2p(ii,iy), | |
2364 | cms 1 v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot | |
2365 | endif | |
2366 | if(nevt(ii).ne.0) | |
2367 | 1 xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii) | |
2368 | c | |
2369 | if(xnpf(ii,iy).gt.1) then | |
2370 | v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy) | |
2371 | v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy) | |
2372 | 1 -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1)) | |
2373 | s2pf(ii,iy)=s2pf(ii,iy)/xnpf(ii,iy) | |
2374 | xnform=dble(xnpf(ii,iy)/dble(nevt(ii))) | |
2375 | etpf(ii,iy)=etpf(ii,iy)/dble(nevt(ii)) | |
2376 | csp | |
2377 | etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii)) | |
2378 | etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii)) | |
2379 | c | |
2380 | cms write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy), | |
2381 | cms 1 v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate | |
2382 | endif | |
2383 | csp | |
2384 | if(xnpl(ii,iy).gt.1) then | |
2385 | v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy) | |
2386 | s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy) | |
2387 | xmult=dble(xnpl(ii,iy)/dble(nevt(ii))) | |
2388 | etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii)) | |
2389 | cms write (nlfile(iy),201) tsp(ii),v2pl(ii,iy), | |
2390 | cms 1 s2pl(ii,iy),etpl(ii,iy),xmult | |
2391 | endif | |
2392 | if(xnps(ii,iy).gt.1) then | |
2393 | v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy) | |
2394 | s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy) | |
2395 | xmult=dble(xnps(ii,iy)/dble(nevt(ii))) | |
2396 | etps(ii,iy)=etps(ii,iy)/dble(nevt(ii)) | |
2397 | cms write (nsfile(iy),201) tsp(ii),v2ps(ii,iy), | |
2398 | cms 1 s2ps(ii,iy),etps(ii,iy),xmult | |
2399 | endif | |
2400 | if(xnplf(ii,iy).gt.1) then | |
2401 | v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy) | |
2402 | s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy) | |
2403 | xmult=dble(xnplf(ii,iy)/dble(nevt(ii))) | |
2404 | etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii)) | |
2405 | cms write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy), | |
2406 | cms 1 s2plf(ii,iy),etplf(ii,iy),xmult | |
2407 | endif | |
2408 | if(xnpsf(ii,iy).gt.1) then | |
2409 | v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy) | |
2410 | s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy) | |
2411 | xmult=dble(xnpsf(ii,iy)/dble(nevt(ii))) | |
2412 | etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii)) | |
2413 | cms write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy), | |
2414 | cms 1 s2psf(ii,iy),etpsf(ii,iy),xmult | |
2415 | endif | |
2416 | csp-end | |
2417 | 1014 continue | |
2418 | 150 continue | |
2419 | csp07/05 initial & final mt distrb | |
2420 | scalei=0d0 | |
2421 | scalef=0d0 | |
2422 | if(nevt(1).ne.0) SCALEi = 1d0 / dble(nevt(1)) / BMT | |
2423 | if(nevt(30).ne.0) SCALEf = 1d0 / dble(nevt(30)) / BMT | |
2424 | do 1016 iy=2,3 | |
2425 | yra = 1d0 | |
2426 | if(iy .eq. 2)yra = 2d0 | |
2427 | do 1015 i=1,50 | |
2428 | cms WRITE(nmfile(iy),251) BMT*dble(I - 0.5), | |
2429 | cms & SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra, | |
2430 | cms & SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra | |
2431 | 1015 continue | |
2432 | 1016 continue | |
2433 | csp07/05 end | |
2434 | clin-3/30/00 event-by-event average & variance of v2: | |
2435 | if(nevt(30).ge.1) then | |
2436 | do 1017 iy=1,3 | |
2437 | v2pavg(iy)=v2psum(iy)/nevt(30) | |
2438 | v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2 | |
2439 | if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0) | |
2440 | 1017 continue | |
2441 | write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg | |
2442 | write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p | |
2443 | endif | |
2444 | clin-11/28/00 for animation: | |
2445 | if(nevent.eq.1) then | |
2446 | do 1018 I=1,mul | |
2447 | if(FT5(I).le.t2time) then | |
2448 | write(93,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i) | |
2449 | endif | |
2450 | 1018 continue | |
2451 | clin-11/29/00 signal the end of animation file: | |
2452 | c write(91,*) -10. | |
2453 | c write(91,*) 0 | |
2454 | c write(92,*) -10. | |
2455 | c write(92,*) 0 | |
2456 | c close (91) | |
2457 | c close (92) | |
2458 | close (93) | |
2459 | endif | |
2460 | clin-5/18/01 calculate v2 for active partons: | |
2461 | do 450 ianp=1,30 | |
2462 | do 400 iy=1,3 | |
2463 | v2pact=0d0 | |
2464 | v2p2ac=0d0 | |
2465 | s2pact=0d0 | |
2466 | etpact=0d0 | |
2467 | xnacti=0d0 | |
2468 | if(xnpf(ianp,iy).gt.1) then | |
2469 | c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons: | |
2470 | v2pact=v2pf(ianp,iy)*xnpf(ianp,iy) | |
2471 | v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1) | |
2472 | 1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy) | |
2473 | s2pact=s2pf(ianp,iy)*xnpf(ianp,iy) | |
2474 | etpact=etpf(ianp,iy)*dble(nevt(ianp)) | |
2475 | xnpact=xnpf(ianp,iy) | |
2476 | c | |
2477 | do 1019 kanp=1,ianp | |
2478 | v2pact=v2pact-v2pfrz(kanp,iy) | |
2479 | v2p2ac=v2p2ac-v2p2fz(kanp,iy) | |
2480 | s2pact=s2pact-s2pfrz(kanp,iy) | |
2481 | etpact=etpact-etpfrz(kanp,iy) | |
2482 | xnpact=xnpact-xnpfrz(kanp,iy) | |
2483 | 1019 continue | |
2484 | c save the sum of v2p, v2p2, s2p, etp, and xnp for formed partons: | |
2485 | v2ph=v2pact | |
2486 | v2ph2=v2p2ac | |
2487 | s2ph=s2pact | |
2488 | etph=etpact | |
2489 | xnp2=xnpact/2d0 | |
2490 | c | |
2491 | if(xnpact.gt.1.and.nevt(ianp).ne.0) then | |
2492 | v2pact=v2pact/xnpact | |
2493 | v2p2ac=dsqrt((v2p2ac/xnpact | |
2494 | 1 -v2pact**2)/(xnpact-1)) | |
2495 | s2pact=s2pact/xnpact | |
2496 | xnacti=dble(xnpact/dble(nevt(ianp))) | |
2497 | etpact=etpact/dble(nevt(ianp)) | |
2498 | cms write (nfile(iy)+2, 250) tsp(ianp),v2pact, | |
2499 | cms 1 v2p2ac,s2pact,etpact,xnacti, | |
2500 | cms 2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp) | |
2501 | endif | |
2502 | endif | |
2503 | c To calculate combined v2 for active partons plus formed hadrons, | |
2504 | c add the sum of v2h, v2h2, s2h, eth, and xnh for formed hadrons: | |
2505 | c scale the hadron part in case nevt(ianp) != nevent: | |
2506 | shadr=dble(nevt(ianp))/dble(nevent) | |
2507 | ianh=ianp | |
2508 | v2ph=v2ph+v2h(ianh,iy)*xnhadr(ianh,iy)*shadr | |
2509 | v2ph2=v2ph2+(v2h2(ianh,iy)**2*(xnhadr(ianh,iy)-1) | |
2510 | 1 +v2h(ianh,iy)**2)*xnhadr(ianh,iy)*shadr | |
2511 | s2ph=s2ph+s2h(ianh,iy)*xnhadr(ianh,iy)*shadr | |
2512 | etph=etph+eth(ianh,iy)*dble(nevent)*shadr | |
2513 | xnph=xnpact+xnhadr(ianh,iy)*shadr | |
2514 | xnp2h=xnp2+xnhadr(ianh,iy)*shadr | |
2515 | if(xnph.gt.1) then | |
2516 | v2ph=v2ph/xnph | |
2517 | v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1)) | |
2518 | s2ph=s2ph/xnph | |
2519 | etph=etph/dble(nevt(ianp)) | |
2520 | xnp2=xnp2/dble(nevt(ianp)) | |
2521 | xnp2h=xnp2h/dble(nevent) | |
2522 | cms if(tsp(ianp).le.(ntmax*dt)) | |
2523 | cms 1 write (nfile(iy)+3, 250) tsp(ianp),v2ph, | |
2524 | cms 2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp) | |
2525 | endif | |
2526 | c | |
2527 | 400 continue | |
2528 | 450 continue | |
2529 | do 550 ianp=1,30 | |
2530 | do 500 iy=1,3 | |
2531 | v2pact=0d0 | |
2532 | v2p2ac=0d0 | |
2533 | s2pact=0d0 | |
2534 | etpact=0d0 | |
2535 | xnacti=0d0 | |
2536 | c reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons: | |
2537 | v2pact=v2pf(ianp,iy)*xnpf(ianp,iy) | |
2538 | v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1) | |
2539 | 1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy) | |
2540 | s2pact=s2pf(ianp,iy)*xnpf(ianp,iy) | |
2541 | etpact=etpf(ianp,iy)*dble(nevt(ianp)) | |
2542 | xnpact=xnpf(ianp,iy) | |
2543 | 500 continue | |
2544 | 550 continue | |
2545 | cms close (620) | |
2546 | cms close (630) | |
2547 | do 1021 nf=1,3 | |
2548 | do 1020 ifile=0,3 | |
2549 | cms close(nfile(nf)+ifile) | |
2550 | 1020 continue | |
2551 | 1021 continue | |
2552 | do 1022 nf=1,3 | |
2553 | cms close(740+nf) | |
2554 | 1022 continue | |
2555 | endif | |
2556 | cyy 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2) | |
2557 | cyy 210 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2)) | |
2558 | 240 format(a30,3(2x,f9.5)) | |
2559 | cyy 250 format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6) | |
2560 | csp | |
2561 | cyy 201 format(2x,f5.2,4(2x,f9.2)) | |
2562 | cyy 251 format(5e15.5) | |
2563 | c | |
2564 | return | |
2565 | end | |
2566 | ||
2567 | c======================================================================= | |
2568 | c Calculate flow from formed hadrons, called by art1e.f: | |
2569 | c Note: numbers in art not in double precision! | |
2570 | subroutine flowh(ct) | |
2571 | PARAMETER (MAXSTR=150001, MAXR=1) | |
2572 | dimension tsh(31) | |
2573 | DOUBLE PRECISION v2h,xnhadr,eth,v2h2,s2h | |
2574 | DOUBLE PRECISION v2hp,xnhadp,v2hsum,v2h2sm,v2hevt(3) | |
2575 | DOUBLE PRECISION pt2, v2hadr | |
2576 | COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3), | |
2577 | 1 v2h2(30,3),s2h(30,3) | |
2578 | cc SAVE /hflow/ | |
2579 | common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3) | |
2580 | cc SAVE /ebe/ | |
2581 | common /lastt/itimeh,bimp | |
2582 | cc SAVE /lastt/ | |
2583 | COMMON /RUN/ NUM | |
2584 | cc SAVE /RUN/ | |
2585 | COMMON /AA/ R(3,MAXSTR) | |
2586 | cc SAVE /AA/ | |
2587 | COMMON /BB/ P(3,MAXSTR) | |
2588 | cc SAVE /BB/ | |
2589 | COMMON /CC/ E(MAXSTR) | |
2590 | cc SAVE /CC/ | |
2591 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
2592 | cc SAVE /EE/ | |
2593 | COMMON /RR/ MASSR(0:MAXR) | |
2594 | cc SAVE /RR/ | |
2595 | common/anim/nevent,isoft,isflag,izpc | |
2596 | cc SAVE /anim/ | |
2597 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
2598 | cc SAVE /AREVT/ | |
2599 | SAVE | |
2600 | c | |
2601 | do 1001 ii = 1, 31 | |
2602 | tsh(ii)=float(ii-1) | |
2603 | 1001 continue | |
2604 | c | |
2605 | do 1004 ianh = 1, 30 | |
2606 | if ((ct+0.0001).lt.tsh(ianh+1) | |
2607 | 1 .and.(ct+0.0001).ge.tsh(ianh)) then | |
2608 | if(ianh.eq.itimeh) goto 101 | |
2609 | IA = 0 | |
2610 | DO 1002 J = 1, NUM | |
2611 | mult=MASSR(J) | |
2612 | IA = IA + MASSR(J - 1) | |
2613 | DO 100 IC = 1, mult | |
2614 | I = IA + IC | |
2615 | c 5/04/01 exclude leptons and photons: | |
2616 | if(iabs(LB(I)-10000).lt.100) goto 100 | |
2617 | px=p(1,i) | |
2618 | py=p(2,i) | |
2619 | pt2=dble(px)**2+dble(py)**2 | |
2620 | c 2/18/00 Note: e(i) gives the mass in ART: | |
2621 | ene=sqrt(e(i)**2+sngl(pt2)+p(3,i)**2) | |
2622 | RAP=0.5*alog((ene+p(3,i))/(ene-p(3,i))) | |
2623 | ctest off density with 2fm radius and z:(-0.1*t,0.1*t): | |
2624 | c rt_now=sqrt(r(1,i)**2+r(2,i)**2) | |
2625 | c gz_now=r(3,i) | |
2626 | c zmax=0.1*ct | |
2627 | c volume=3.1416*(2.**2)*2*zmax | |
2628 | c if(rt_now.gt.2.or.abs(gz_now).gt.zmax) | |
2629 | c 1 goto 100 | |
2630 | iloop=1 | |
2631 | if(abs(rap).le.1) then | |
2632 | iloop=2 | |
2633 | if(abs(rap).le.0.5) then | |
2634 | iloop=3 | |
2635 | endif | |
2636 | endif | |
2637 | do 50 iy=1,iloop | |
2638 | if(pt2.gt.0d0) then | |
2639 | v2hadr=(dble(px)**2-dble(py)**2)/pt2 | |
2640 | v2h(ianh,iy)=v2h(ianh,iy)+v2hadr | |
2641 | v2h2(ianh,iy)=v2h2(ianh,iy)+v2hadr**2 | |
2642 | if(dabs(v2hadr).gt.1d0) | |
2643 | 1 write(1,*) 'v2hadr>1',v2hadr,px,py | |
2644 | endif | |
2645 | xperp2=r(1,I)**2+r(2,I)**2 | |
2646 | if(xperp2.gt.0.) | |
2647 | 1 s2h(ianh,iy)=s2h(ianh,iy)+dble((r(1,I)**2-r(2,I)**2)/xperp2) | |
2648 | eth(ianh,iy)=eth(ianh,iy)+dble(SQRT(e(i)**2+sngl(pt2))) | |
2649 | ctest off density: | |
2650 | c eth(ianh,iy)=eth(ianh,iy) | |
2651 | c 1 +dble(SQRT(e(i)**2+sngl(pt2)+p(3,i)**2))/volume | |
2652 | xnhadr(ianh,iy)=xnhadr(ianh,iy)+1d0 | |
2653 | 50 continue | |
2654 | 100 continue | |
2655 | 1002 CONTINUE | |
2656 | itimeh=ianh | |
2657 | clin-5/04/01 ebe v2 variables: | |
2658 | if(ianh.eq.30) then | |
2659 | do 1003 iy=1,3 | |
2660 | nhadrn=IDINT(xnhadr(ianh,iy)-xnhadp(iy)) | |
2661 | if(nhadrn.ne.0) then | |
2662 | v2hevt(iy)=(v2h(ianh,iy)-v2hp(iy))/dble(nhadrn) | |
2663 | v2hsum(iy)=v2hsum(iy)+v2hevt(iy) | |
2664 | v2h2sm(iy)=v2h2sm(iy)+v2hevt(iy)**2 | |
2665 | v2hp(iy)=v2h(ianh,iy) | |
2666 | xnhadp(iy)=xnhadr(ianh,iy) | |
2667 | endif | |
2668 | 1003 continue | |
2669 | write(88, 160) iaevt,v2hevt | |
2670 | endif | |
2671 | goto 101 | |
2672 | endif | |
2673 | 1004 continue | |
2674 | 160 format(i10,3(2x,f9.5)) | |
2675 | clin-11/27/00 for animation: | |
2676 | 101 if(nevent.eq.1) then | |
2677 | IA = 0 | |
2678 | do 1005 J = 1, NUM | |
2679 | mult=MASSR(J) | |
2680 | IA = IA + MASSR(J - 1) | |
2681 | c write(91,*) ct | |
2682 | c write(91,*) mult | |
2683 | DO 150 IC = 1, mult | |
2684 | I = IA + IC | |
2685 | c write(91,210) LB(i),r(1,i),r(2,i),r(3,i), | |
2686 | c 1 p(1,i),p(2,i),p(3,i),e(i) | |
2687 | 150 continue | |
2688 | 1005 continue | |
2689 | return | |
2690 | endif | |
2691 | c 210 format(i6,7(1x,f8.3)) | |
2692 | return | |
2693 | end | |
2694 | ||
2695 | c======================================================================= | |
2696 | subroutine flowh0(NEVNT,idd) | |
2697 | c | |
2698 | dimension tsh(31) | |
2699 | DOUBLE PRECISION v2h,xnhadr,eth,v2h2,s2h | |
2700 | DOUBLE PRECISION v2hp,xnhadp,v2hsum,v2h2sm, | |
2701 | 1 v2havg(3),varv2h(3) | |
2702 | COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3), | |
2703 | 1 v2h2(30,3),s2h(30,3) | |
2704 | cc SAVE /hflow/ | |
2705 | common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3) | |
2706 | cc SAVE /ebe/ | |
2707 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
2708 | cc SAVE /input1/ | |
2709 | COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, | |
2710 | & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL | |
2711 | cc SAVE /INPUT2/ | |
2712 | common /lastt/itimeh,bimp | |
2713 | cc SAVE /lastt/ | |
2714 | SAVE | |
2715 | ||
2716 | c idd=0: initialization for flow analysis, called by artdri.f:: | |
2717 | if(idd.eq.0) then | |
2718 | itimeh=0 | |
2719 | c | |
2720 | do 1001 ii = 1, 31 | |
2721 | tsh(ii)=float(ii-1) | |
2722 | 1001 continue | |
2723 | c | |
2724 | do 1003 ii=1,30 | |
2725 | do 1002 iy=1,3 | |
2726 | v2h(ii,iy)=0d0 | |
2727 | xnhadr(ii,iy)=0d0 | |
2728 | eth(ii,iy)=0d0 | |
2729 | v2h2(ii,iy)=0d0 | |
2730 | s2h(ii,iy)=0d0 | |
2731 | 1002 continue | |
2732 | 1003 continue | |
2733 | do 1004 iy=1,3 | |
2734 | v2hp(iy)=0d0 | |
2735 | xnhadp(iy)=0d0 | |
2736 | v2hsum(iy)=0d0 | |
2737 | v2h2sm(iy)=0d0 | |
2738 | if(iy.eq.1) then | |
2739 | nunit=59 | |
2740 | elseif(iy.eq.2) then | |
2741 | nunit=68 | |
2742 | else | |
2743 | nunit=69 | |
2744 | endif | |
2745 | write(nunit,*) ' tsh, v2h, v2h2, s2h, '// | |
2746 | 1 ' eth, xmulth' | |
2747 | 1004 continue | |
2748 | c idd=2: calculate average hadronic elliptic flow, called by artdri.f: | |
2749 | else if(idd.eq.2) then | |
2750 | do 100 ii=1, 30 | |
2751 | do 1005 iy=1,3 | |
2752 | if(xnhadr(ii,iy).eq.0) then | |
2753 | xmulth=0. | |
2754 | elseif(xnhadr(ii,iy).gt.1) then | |
2755 | v2h(ii,iy)=v2h(ii,iy)/xnhadr(ii,iy) | |
2756 | eth(ii,iy)=eth(ii,iy)/dble(NEVNT) | |
2757 | v2h2(ii,iy)=dsqrt((v2h2(ii,iy)/xnhadr(ii,iy) | |
2758 | 1 -v2h(ii,iy)**2)/(xnhadr(ii,iy)-1)) | |
2759 | s2h(ii,iy)=s2h(ii,iy)/xnhadr(ii,iy) | |
2760 | xmulth=sngl(xnhadr(ii,iy)/NEVNT) | |
2761 | endif | |
2762 | if(iy.eq.1) then | |
2763 | nunit=59 | |
2764 | elseif(iy.eq.2) then | |
2765 | nunit=68 | |
2766 | else | |
2767 | nunit=69 | |
2768 | endif | |
2769 | if(tsh(ii).le.(ntmax*dt)) | |
2770 | 1 write (nunit,200) tsh(ii),v2h(ii,iy), | |
2771 | 2 v2h2(ii,iy),s2h(ii,iy),eth(ii,iy),xmulth | |
2772 | 1005 continue | |
2773 | 100 continue | |
2774 | c event-by-event average & variance of v2h: | |
2775 | do 1006 iy=1,3 | |
2776 | v2havg(iy)=v2hsum(iy)/dble(NEVNT) | |
2777 | varv2h(iy)=dsqrt(v2h2sm(iy)/dble(NEVNT)-v2havg(iy)**2) | |
2778 | 1006 continue | |
2779 | write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): avg=', v2havg | |
2780 | write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): var=', varv2h | |
2781 | endif | |
2782 | 200 format(2x,f5.2,3(2x,f7.4),2(2x,f9.2)) | |
2783 | 240 format(a30,3(2x,f9.5)) | |
2784 | return | |
2785 | end | |
2786 | ||
2787 | c======================================================================= | |
2788 | c 2/23/00 flow from all initial hadrons just before entering ARTMN: | |
2789 | subroutine iniflw(NEVNT,idd) | |
2790 | PARAMETER (MAXSTR=150001, MAXR=1) | |
2791 | DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult, | |
2792 | 1 v2bi,s2bi,xbmult | |
2793 | COMMON /RUN/ NUM | |
2794 | cc SAVE /RUN/ | |
2795 | COMMON /ARERC1/MULTI1(MAXR) | |
2796 | cc SAVE /ARERC1/ | |
2797 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
2798 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
2799 | & FT1(MAXSTR, MAXR), | |
2800 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
2801 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
2802 | cc SAVE /ARPRC1/ | |
2803 | COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult | |
2804 | cc SAVE /iflow/ | |
2805 | SAVE | |
2806 | c | |
2807 | if(idd.eq.0) then | |
2808 | v2i=0d0 | |
2809 | eti=0d0 | |
2810 | xmulti=0d0 | |
2811 | v2mi=0d0 | |
2812 | s2mi=0d0 | |
2813 | xmmult=0d0 | |
2814 | v2bi=0d0 | |
2815 | s2bi=0d0 | |
2816 | xbmult=0d0 | |
2817 | else if(idd.eq.1) then | |
2818 | do 1002 J = 1, NUM | |
2819 | do 1001 I = 1, MULTI1(J) | |
2820 | ITYP = ITYP1(I, J) | |
2821 | c all hadrons: | |
2822 | IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 100 | |
2823 | xmulti=xmulti+1.d0 | |
2824 | PX = PX1(I, J) | |
2825 | PY = PY1(I, J) | |
2826 | XM = XM1(I, J) | |
2827 | pt2=px**2+py**2 | |
2828 | xh=gx1(I,J) | |
2829 | yh=gy1(I,J) | |
2830 | xt2=xh**2+yh**2 | |
2831 | if(pt2.gt.0) v2i=v2i+dble((px**2-py**2)/pt2) | |
2832 | eti=eti+dble(SQRT(PX ** 2 + PY ** 2 + XM ** 2)) | |
2833 | c baryons only: | |
2834 | IF (ITYP .LT. -1000 .or. ITYP .GT. 1000) then | |
2835 | xbmult=xbmult+1.d0 | |
2836 | if(pt2.gt.0) v2bi=v2bi+dble((px**2-py**2)/pt2) | |
2837 | if(xt2.gt.0) s2bi=s2bi+dble((xh**2-yh**2)/xt2) | |
2838 | c mesons only: | |
2839 | else | |
2840 | xmmult=xmmult+1.d0 | |
2841 | if(pt2.gt.0) v2mi=v2mi+dble((px**2-py**2)/pt2) | |
2842 | if(xt2.gt.0) s2mi=s2mi+dble((xh**2-yh**2)/xt2) | |
2843 | endif | |
2844 | 100 continue | |
2845 | 1001 continue | |
2846 | 1002 continue | |
2847 | else if(idd.eq.2) then | |
2848 | if(xmulti.ne.0) v2i=v2i/xmulti | |
2849 | eti=eti/dble(NEVNT) | |
2850 | xmulti=xmulti/dble(NEVNT) | |
2851 | if(xmmult.ne.0) then | |
2852 | v2mi=v2mi/xmmult | |
2853 | s2mi=s2mi/xmmult | |
2854 | endif | |
2855 | xmmult=xmmult/dble(NEVNT) | |
2856 | if(xbmult.ne.0) then | |
2857 | v2bi=v2bi/xbmult | |
2858 | s2bi=s2bi/xbmult | |
2859 | endif | |
2860 | xbmult=xbmult/dble(NEVNT) | |
2861 | endif | |
2862 | c | |
2863 | return | |
2864 | end | |
2865 | ||
2866 | c======================================================================= | |
2867 | c 2/25/00 dN/dt analysis for production (before ZPCMN) | |
2868 | c and freezeout (right after ZPCMN) for all partons. | |
2869 | subroutine frztm(NEVNT,idd) | |
2870 | c | |
2871 | implicit double precision (a-h, o-z) | |
2872 | PARAMETER (MAXPTN=400001) | |
2873 | dimension tsf(31) | |
2874 | COMMON /PARA1/ MUL | |
2875 | cc SAVE /PARA1/ | |
2876 | COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN), | |
2877 | & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN), | |
2878 | & XMASS0(MAXPTN), ITYP0(MAXPTN) | |
2879 | cc SAVE /prec1/ | |
2880 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
2881 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
2882 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
2883 | cc SAVE /prec2/ | |
2884 | COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30), | |
2885 | & dnprod(30),detpro(30),dnfrz(30),detfrz(30) | |
2886 | cc SAVE /frzout/ | |
2887 | SAVE | |
2888 | data tsf/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, | |
2889 | & 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, | |
2890 | & 2 , 3, 4, 5, 6, 7, 8, 9, 10, 20, 30/ | |
2891 | c | |
2892 | if(idd.eq.0) then | |
2893 | do 1001 ii=1,30 | |
2894 | xnprod(ii)=0d0 | |
2895 | etprod(ii)=0d0 | |
2896 | xnfrz(ii)=0d0 | |
2897 | etfrz(ii)=0d0 | |
2898 | dnprod(ii)=0d0 | |
2899 | detpro(ii)=0d0 | |
2900 | dnfrz(ii)=0d0 | |
2901 | detfrz(ii)=0d0 | |
2902 | 1001 continue | |
2903 | cms OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN') | |
2904 | cms OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN') | |
2905 | else if(idd.eq.1) then | |
2906 | DO 100 ip = 1, MUL | |
2907 | do 1002 ii=1,30 | |
2908 | eth0=dSQRT(PX0(ip)**2+PY0(ip)**2+XMASS0(ip)**2) | |
2909 | eth2=dSQRT(PX5(ip)**2+PY5(ip)**2+XMASS5(ip)**2) | |
2910 | c total number and Et produced by time tsf(ii): | |
2911 | if (ft0(ip).lt.tsf(ii+1)) then | |
2912 | xnprod(ii)=xnprod(ii)+1d0 | |
2913 | etprod(ii)=etprod(ii)+eth0 | |
2914 | c number and Et produced from time tsf(ii) to tsf(ii+1): | |
2915 | if (ft0(ip).ge.tsf(ii)) then | |
2916 | dnprod(ii)=dnprod(ii)+1d0 | |
2917 | detpro(ii)=detpro(ii)+eth0 | |
2918 | endif | |
2919 | endif | |
2920 | c total number and Et freezed out by time tsf(ii): | |
2921 | if (FT5(ip).lt.tsf(ii+1)) then | |
2922 | xnfrz(ii)=xnfrz(ii)+1d0 | |
2923 | etfrz(ii)=etfrz(ii)+eth2 | |
2924 | c number and Et freezed out from time tsf(ii) to tsf(ii+1): | |
2925 | if (FT5(ip).ge.tsf(ii)) then | |
2926 | dnfrz(ii)=dnfrz(ii)+1d0 | |
2927 | detfrz(ii)=detfrz(ii)+eth2 | |
2928 | endif | |
2929 | endif | |
2930 | 1002 continue | |
2931 | 100 continue | |
2932 | else if(idd.eq.2) then | |
2933 | cms write (86,*) ' t, np, dnp/dt, etp '// | |
2934 | cms 1 ' detp/dt' | |
2935 | cms write (87,*) ' t, nf, dnf/dt, etf '// | |
2936 | cms 1 ' detf/dt' | |
2937 | do 1003 ii=1,30 | |
2938 | xnp=xnprod(ii)/dble(NEVNT) | |
2939 | xnf=xnfrz(ii)/dble(NEVNT) | |
2940 | etp=etprod(ii)/dble(NEVNT) | |
2941 | etf=etfrz(ii)/dble(NEVNT) | |
2942 | dxnp=dnprod(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) | |
2943 | dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) | |
2944 | detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) | |
2945 | detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii)) | |
2946 | cms write (86, 200) | |
2947 | cms 1 tsf(ii+1),xnp,dxnp,etp,detp | |
2948 | cms write (87, 200) | |
2949 | cms 1 tsf(ii+1),xnf,dxnf,etf,detf | |
2950 | 1003 continue | |
2951 | endif | |
2952 | cyy 200 format(2x,f9.2,4(2x,f10.2)) | |
2953 | c | |
2954 | return | |
2955 | end | |
2956 | ||
2957 | c======================================================================= | |
2958 | clin-6/2009 write out initial minijet information | |
2959 | c before propagating to its formation time: | |
2960 | subroutine minijet_out(BB) | |
2961 | PARAMETER (MAXSTR=150001) | |
2962 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
2963 | COMMON/hjcrdn/YP(3,300),YT(3,300) | |
2964 | COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500), | |
2965 | & PJPY(300,500),PJPZ(300,500),PJPE(300,500), | |
2966 | & PJPM(300,500),NTJ(300),KFTJ(300,500), | |
2967 | & PJTX(300,500),PJTY(300,500),PJTZ(300,500), | |
2968 | & PJTE(300,500),PJTM(300,500) | |
2969 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
2970 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
2971 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
2972 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
2973 | common /para7/ ioscar,nsmbbbar,nsmmeson | |
2974 | common/phidcy/iphidcy,pttrig,ntrig,maxmiss | |
2975 | SAVE | |
2976 | ntrig=0 | |
2977 | do I = 1, IHNT2(1) | |
2978 | do J = 1, NPJ(I) | |
2979 | pt=sqrt(PJPX(I,J)**2+PJPY(I,J)**2) | |
2980 | if(pt.ge.pttrig) ntrig=ntrig+1 | |
2981 | enddo | |
2982 | enddo | |
2983 | do I = 1, IHNT2(3) | |
2984 | do J = 1, NTJ(I) | |
2985 | pt=sqrt(PJTX(I,J)**2+PJTY(I,J)**2) | |
2986 | if(pt.ge.pttrig) ntrig=ntrig+1 | |
2987 | enddo | |
2988 | enddo | |
2989 | do I = 1, NSG | |
2990 | do J = 1, NJSG(I) | |
2991 | pt=sqrt(PXSG(I,J)**2+PYSG(I,J)**2) | |
2992 | if(pt.ge.pttrig) ntrig=ntrig+1 | |
2993 | enddo | |
2994 | enddo | |
2995 | c Require at least 1 initial minijet parton above the trigger Pt value: | |
2996 | if(ntrig.eq.0) return | |
2997 | ||
2998 | c.....transfer data from HIJING to ZPC | |
2999 | if(ioscar.eq.3) write(96,*) IAEVT,MISS,IHNT2(1),IHNT2(3) | |
3000 | DO 1008 I = 1, IHNT2(1) | |
3001 | DO 1007 J = 1, NPJ(I) | |
3002 | ityp=KFPJ(I,J) | |
3003 | c write out not only gluons: | |
3004 | c if(ityp.ne.21) goto 1007 | |
3005 | gx=YP(1,I)+0.5*BB | |
3006 | gy=YP(2,I) | |
3007 | gz=0d0 | |
3008 | ft=0d0 | |
3009 | px=PJPX(I,J) | |
3010 | py=PJPY(I,J) | |
3011 | pz=PJPZ(I,J) | |
3012 | xmass=PJPM(I,J) | |
3013 | if(ioscar.eq.3) then | |
3014 | if(amax1(abs(gx),abs(gy), | |
3015 | 1 abs(gz),abs(ft)).lt.9999) then | |
3016 | write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,1 | |
3017 | else | |
3018 | write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,1 | |
3019 | endif | |
3020 | endif | |
3021 | 1007 CONTINUE | |
3022 | 1008 CONTINUE | |
3023 | DO 1010 I = 1, IHNT2(3) | |
3024 | DO 1009 J = 1, NTJ(I) | |
3025 | ityp=KFTJ(I,J) | |
3026 | c if(ityp.ne.21) goto 1009 | |
3027 | gx=YT(1,I)-0.5*BB | |
3028 | gy=YT(2,I) | |
3029 | gz=0d0 | |
3030 | ft=0d0 | |
3031 | px=PJTX(I,J) | |
3032 | py=PJTY(I,J) | |
3033 | pz=PJTZ(I,J) | |
3034 | xmass=PJTM(I,J) | |
3035 | if(ioscar.eq.3) then | |
3036 | if(amax1(abs(gx),abs(gy), | |
3037 | 1 abs(gz),abs(ft)).lt.9999) then | |
3038 | write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,2 | |
3039 | else | |
3040 | write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,2 | |
3041 | endif | |
3042 | endif | |
3043 | 1009 CONTINUE | |
3044 | 1010 CONTINUE | |
3045 | DO 1012 I = 1, NSG | |
3046 | DO 1011 J = 1, NJSG(I) | |
3047 | ityp=K2SG(I,J) | |
3048 | c if(ityp.ne.21) goto 1011 | |
3049 | gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2))) | |
3050 | gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2))) | |
3051 | gz=0d0 | |
3052 | ft=0d0 | |
3053 | px=PXSG(I,J) | |
3054 | py=PYSG(I,J) | |
3055 | pz=PZSG(I,J) | |
3056 | xmass=PMSG(I,J) | |
3057 | if(ioscar.eq.3) then | |
3058 | if(amax1(abs(gx),abs(gy), | |
3059 | 1 abs(gz),abs(ft)).lt.9999) then | |
3060 | write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,3 | |
3061 | else | |
3062 | write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,3 | |
3063 | endif | |
3064 | endif | |
3065 | 1011 CONTINUE | |
3066 | 1012 CONTINUE | |
3067 | 200 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,f8.2),2(2x,f2.0),2x,I2) | |
3068 | 201 format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,e8.2),2(2x,f2.0),2x,I2) | |
3069 | c | |
3070 | return | |
3071 | end | |
3072 | ||
3073 | c======================================================================= | |
3074 | clin-6/2009 embed back-to-back high-Pt quark/antiquark pair | |
3075 | c via embedding back-to-back high-Pt pion pair then melting the pion pair | |
3076 | c by generating the internal quark and antiquark momentum parallel to | |
3077 | c the pion momentum (in order to produce a high-Pt and a low Pt parton): | |
3078 | subroutine embedHighPt | |
3079 | PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140) | |
3080 | common/embed/iembed,pxqembd,pyqembd,xembd,yembd | |
3081 | COMMON/RNDF77/NSEED | |
3082 | COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11 | |
3083 | COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4) | |
3084 | COMMON /ARPRC/ ITYPAR(MAXSTR), | |
3085 | & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR), | |
3086 | & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR), | |
3087 | & XMAR(MAXSTR) | |
3088 | cwei DOUBLE PRECISION PATT | |
3089 | SAVE | |
3090 | c | |
3091 | if(iembed.ne.1) return | |
3092 | ptq=sqrt(pxqembd**2+pyqembd**2) | |
3093 | if(ptq.lt.(pichmass/2.)) then | |
3094 | print *, 'Embedded quark transverse momentum is too small' | |
3095 | stop | |
3096 | endif | |
3097 | c Randomly embed u/ubar or d/dbar at high Pt: | |
3098 | idqembd=1+int(2*RANART(NSEED)) | |
3099 | c Flavor content for the charged pion that contains the leading quark: | |
3100 | if(idqembd.eq.1) then | |
3101 | idqsoft=-2 | |
3102 | idpi=-211 | |
3103 | elseif(idqembd.eq.2) then | |
3104 | idqsoft=-1 | |
3105 | idpi=211 | |
3106 | else | |
3107 | print *, 'Wrong quark flavor embedded' | |
3108 | stop | |
3109 | endif | |
3110 | c Caculate transverse momentum of the parent charged pion: | |
3111 | xmq=ulmass(idqembd) | |
3112 | xmqsoft=ulmass(idqsoft) | |
3113 | ptpi=((pichmass**2+xmq**2-xmqsoft**2)*ptq | |
3114 | 1 -sqrt((xmq**2+ptq**2)*(pichmass**4 | |
3115 | 2 -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2))) | |
3116 | 3 /(2.*xmq**2) | |
3117 | pxpi=ptpi*pxqembd/ptq | |
3118 | pypi=ptpi*pyqembd/ptq | |
3119 | c Embedded quark/antiquark are assumed to have pz=0: | |
3120 | pzpi=0. | |
3121 | c Insert the two parent charged pions, | |
3122 | c ipion=1 for the pion containing the leading quark, | |
3123 | c ipion=2 for the pion containing the leading antiquark of the same flavor: | |
3124 | do ipion=1,2 | |
3125 | if(ipion.eq.2) then | |
3126 | idpi=-idpi | |
3127 | pxpi=-pxpi | |
3128 | pypi=-pypi | |
3129 | pzpi=-pzpi | |
3130 | endif | |
3131 | NATT=NATT+1 | |
3132 | KATT(NATT,1)=idpi | |
3133 | KATT(NATT,2)=40 | |
3134 | KATT(NATT,3)=0 | |
3135 | PATT(NATT,1)=pxpi | |
3136 | PATT(NATT,2)=pypi | |
3137 | PATT(NATT,3)=pzpi | |
3138 | PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2) | |
3139 | EATT=EATT+PATT(NATT,4) | |
3140 | GXAR(NATT)=xembd | |
3141 | GYAR(NATT)=yembd | |
3142 | GZAR(NATT)=0. | |
3143 | FTAR(NATT)=0. | |
3144 | ITYPAR(NATT)=KATT(NATT,1) | |
3145 | PXAR(NATT)=PATT(NATT,1) | |
3146 | PYAR(NATT)=PATT(NATT,2) | |
3147 | PZAR(NATT)=PATT(NATT,3) | |
3148 | PEAR(NATT)=PATT(NATT,4) | |
3149 | XMAR(NATT)=pichmass | |
3150 | enddo | |
3151 | c | |
3152 | return | |
3153 | end |