]>
Commit | Line | Data |
---|---|---|
0119ef9a | 1 | c....................amptsub.f |
2 | c.....this file contains 4 sections: | |
3 | c.....1. ART subroutines; | |
4 | c.....2. ART functions; | |
5 | c.....3. ART block data; | |
6 | c.....4. subprocesses borrowed from other codes. | |
7 | c.....5. the previous artana.f | |
8 | c.....6. the previous zpcsub.f | |
9 | c.....7. subroutine getnp | |
10 | c.....Note that Parts1-4 are the previous artsub.f | |
11 | c | |
12 | c======================================================================= | |
13 | c.....subroutine to set up ART parameters and analysis files | |
14 | c.....before looping different events | |
15 | cms | |
16 | cms dlw & gsfs 8/2009 commented out lots of output files | |
17 | cms | |
18 | SUBROUTINE ARTSET | |
19 | c | |
20 | PARAMETER (AMU= 0.9383) | |
21 | double precision dpcoal,drcoal,ecritl | |
22 | INTEGER ZTA, ZPR | |
23 | common /gg/ dx,dy,dz,dpx,dpy,dpz | |
24 | clin-10/03/03 | |
25 | c "SAVE " (without argument) is used for most subroutines and functions, | |
26 | c this is important for the success when using "f77" to compile: | |
27 | cc SAVE /gg/ | |
28 | common /zz/ zta,zpr | |
29 | cc SAVE /zz/ | |
30 | COMMON /RUN/ NUM | |
31 | cc SAVE /RUN/ | |
32 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
33 | cc SAVE /input1/ | |
34 | COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, | |
35 | & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL | |
36 | cc SAVE /INPUT2/ | |
37 | COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX | |
38 | cc SAVE /INPUT3/ | |
39 | common /imulst/ iperts | |
40 | cc SAVE /imulst/ | |
41 | common /coal/dpcoal,drcoal,ecritl | |
42 | common/anim/nevent,isoft,isflag,izpc | |
43 | common /para7/ ioscar,nsmbbbar,nsmmeson | |
44 | SAVE | |
45 | clin-10/03/03 ecritl: local energy density below which a parton | |
46 | c will freeze out (in GeV/fm^3), for improvements on string melting, | |
47 | c not used in this version of AMPT: | |
48 | clin-4/2008 | |
49 | c data ecritl/1.d0/ | |
50 | ecritl=1.d0 | |
51 | c | |
52 | c combine ART initialization into ampt.ini: | |
53 | c (Note that the following values are relics from the old ART structure) | |
54 | c.....input parameter file | |
55 | c OPEN(13, FILE = 'art1.ini', STATUS = 'UNKNOWN') | |
56 | c READ (13, *) MASSTA, ZTA | |
57 | MASSTA=1 | |
58 | ZTA=1 | |
59 | c write(12,*) massta, zta, ' massta, zta' | |
60 | c READ (13, *) MASSPR, ZPR | |
61 | MASSPR=1 | |
62 | ZPR=1 | |
63 | c write(12,*) masspr, zpr, ' masspr, zpr' | |
64 | c READ (13, *) PLAB, IPLAB | |
65 | PLAB=14.6 | |
66 | IPLAB=2 | |
67 | c write(12,*) plab, iplab, ' plab, iplab' | |
68 | if(iplab.eq.2)then | |
69 | elab=sqrt(plab**2+amu**2)-amu | |
70 | else | |
71 | elab=plab | |
72 | endif | |
73 | elab=elab*1000. | |
74 | c READ (13, *) ZEROPT | |
75 | ZEROPT=0. | |
76 | c write(12,*) zeropt, ' zeropt' | |
77 | clin-10/03/03 ISEED was used as a seed for random number inside ART, | |
78 | c not used in AMPT: | |
79 | ISEED=700721 | |
80 | c 0/1: (Normal or Perturbative) multistrange partice production. | |
81 | c Perturbative option is disabled for now: | |
82 | iperts=0 | |
83 | c READ (13, *) MANYB, B0, BI, BM | |
84 | c 2/04/00 MANYB MUST BE SET TO 1 ! | |
85 | c in order to skip impact parameter setting by ART, then B0 has no effect. | |
86 | MANYB=1 | |
87 | B0=1 | |
88 | BI=0 | |
89 | BM=0 | |
90 | c write(12,*) manyb, b0, bi, bm, ' manyb, b0, bi, bm' | |
91 | c READ (13, *) ISEED | |
92 | c write(12,*) iseed, ' iseed' | |
93 | c READ (13, *) DT | |
94 | c write(12,*) dt, ' dt' | |
95 | c READ (13, *) NTMAX | |
96 | c write(12,*) ntmax, ' ntmax' | |
97 | c READ (13, *) ICOLL | |
98 | ICOLL=-1 | |
99 | c write(12,*) icoll, ' icoll' | |
100 | c READ (13, *) NUM | |
101 | c 2/11/03 run events without test particles for now: | |
102 | NUM=1 | |
103 | c write(12,*) num, ' num' | |
104 | c READ (13, *) INSYS | |
105 | INSYS=1 | |
106 | c write(12,*) insys, ' insys' | |
107 | c READ (13, *) IPOT | |
108 | IPOT=3 | |
109 | c write(12,*) ipot, ' ipot' | |
110 | c READ (13, *) MODE | |
111 | MODE=0 | |
112 | IF(ICOLL.EQ.-1)IPOT=0 | |
113 | c write(12,*) mode, ' mode' | |
114 | c READ (13, *) DX, DY, DZ | |
115 | DX=2.73 | |
116 | DY=2.73 | |
117 | DZ=2.73 | |
118 | c write(12,*) dx,dy,dz,' dx,dy,dz' | |
119 | c READ (13, *) DPX, DPY, DPZ | |
120 | DPX=0.6 | |
121 | DPY=0.6 | |
122 | DPZ=0.6 | |
123 | c write(12,*) dpx,dpy,dpz,' dpx,dpy,dpz' | |
124 | c READ (13, *) IAVOID | |
125 | IAVOID=1 | |
126 | c write(12,*) iavoid, ' iavoid' | |
127 | c READ (13, *) IMOMEN | |
128 | IMOMEN=1 | |
129 | c write(12,*) imomen, ' imomen' | |
130 | if(icoll.eq.-1)imomen=3 | |
131 | c READ (13, *) NFREQ | |
132 | NFREQ=10 | |
133 | c write(12,*) nfreq, ' nfreq' | |
134 | c READ (13, *) ICFLOW | |
135 | ICFLOW=0 | |
136 | c write(12,*) ICFLOW, ' ICFLOW' | |
137 | c READ (13, *) ICRHO | |
138 | ICRHO=0 | |
139 | c write(12,*) ICRHO, ' ICRHO' | |
140 | c READ (13, *) ICOU | |
141 | ICOU=0 | |
142 | c write(12,*)icou, ' icou' | |
143 | * kaon potential control parameter | |
144 | * KMUL IS A MULTIPLIER TO THE STANDARD K-N SCATTERING LENGTH | |
145 | c READ (13, *) KPOTEN, KMUL | |
146 | KPOTEN=0 | |
147 | KMUL=1 | |
148 | c write(12,*)kpoten,kmul, ' kpoten, kmul' | |
149 | * mean field control parameter FOR BARYONS | |
150 | * no mean filed is used for baryons if their | |
151 | * local density is higher than dencut. | |
152 | c READ (13, *) DENCUT | |
153 | DENCUT=15 | |
154 | c write(12,*)dencut, ' dencut' | |
155 | * test reactions in a box of side-length cycbox | |
156 | * input cycbox | |
157 | c READ (13, *) CYCBOX | |
158 | CYCBOX=0 | |
159 | c write(12,*) cycbox, ' cycbox' | |
160 | c | |
161 | clin-5b/2008 | |
162 | c if(ioscar.eq.2) then | |
163 | if(ioscar.eq.2.or.ioscar.eq.3) then | |
164 | cms OPEN (92,FILE='ana/parton-initial-afterPropagation.dat', | |
165 | cms 1 STATUS = 'UNKNOWN') | |
166 | endif | |
167 | if(ioscar.eq.3) then | |
168 | clin-6/2009 write out full parton collision history: | |
169 | cms OPEN (95,FILE='ana/parton-collisionsHistory.dat', | |
170 | cms 1 STATUS='UNKNOWN') | |
171 | clin-6/2009 write out initial minijet information: | |
172 | cms OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat', | |
173 | cms 1 STATUS='UNKNOWN') | |
174 | clin-6/2009 write out parton info after coalescence: | |
175 | if(isoft.eq.4.or.isoft.eq.5) then | |
176 | cms OPEN (85,FILE='ana/parton-after-coalescence.dat', | |
177 | cms 1 STATUS='UNKNOWN') | |
178 | endif | |
179 | endif | |
180 | clin-6/2009 write out initial transverse positions of initial nucleons: | |
181 | cms OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN') | |
182 | ||
183 | RETURN | |
184 | END | |
185 | ||
186 | c----------------------------------------------------------------------- | |
187 | ||
188 | c.....subroutine to initialize cascade. | |
189 | ||
190 | SUBROUTINE ARINI | |
191 | ||
192 | c.....before invoking ARINI: | |
193 | c.....IAPAR2(1), IAINT2(1) must be set. | |
194 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
195 | cc SAVE /ARPRNT/ | |
196 | SAVE | |
197 | ||
198 | ctest off for resonance (phi, K*) studies: | |
199 | c OPEN (89, FILE = 'ana/decay_rec.dat', STATUS = 'UNKNOWN') | |
200 | ||
201 | IFLG = IAPAR2(1) | |
202 | GOTO (200, 200, 300) IFLG | |
203 | ||
204 | c.....error choice of initialization | |
205 | PRINT *, 'IAPAR2(1) must be 1, 2, or 3' | |
206 | STOP | |
207 | ||
208 | c.....to use default initial conditions generated by the cascade, | |
209 | c.....or to read in initial conditions. | |
210 | 200 RETURN | |
211 | ||
212 | c.....to generate formation time and the position at formation time from | |
213 | c.....read-in initial conditions with an averaged formation proper time. | |
214 | 300 CALL ARINI1 | |
215 | c.....ordering the particle label according to increasing order of | |
216 | c.....formation time. | |
217 | CALL ARTORD | |
218 | RETURN | |
219 | ||
220 | END | |
221 | ||
222 | c----------------------------------------------------------------------- | |
223 | ||
224 | c.....subroutine to generate formation time and position at formation time | |
225 | c.....from read-in initial conditions with an averaged formation proper | |
226 | c.....time. | |
227 | ||
228 | SUBROUTINE ARINI1 | |
229 | ||
230 | c.....before invoking ARINI1: | |
231 | c.....ARPAR1(1), IAINT2(1) must be set: | |
232 | PARAMETER (MAXSTR=150001) | |
233 | double precision smearp,smearh | |
234 | ||
235 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
236 | cc SAVE /ARPRNT/ | |
237 | COMMON /ARPRC/ ITYPAR(MAXSTR), | |
238 | & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR), | |
239 | & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR), | |
240 | & XMAR(MAXSTR) | |
241 | cc SAVE /ARPRC/ | |
242 | COMMON /smearz/smearp,smearh | |
243 | cc SAVE /smearz/ | |
244 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
245 | cc SAVE /input1/ | |
246 | common/anim/nevent,isoft,isflag,izpc | |
247 | cc SAVE /anim/ | |
248 | common /nzpc/nattzp | |
249 | cc SAVE /nzpc/ | |
250 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
251 | cc SAVE /HPARNT/ | |
252 | COMMON/RNDF77/NSEED | |
253 | cc SAVE /RNDF77/ | |
254 | common /para8/ idpert,npertd,idxsec | |
255 | SAVE | |
256 | clin-5/2008 for perturbatively-produced hadrons (currently only deuterons): | |
257 | cms OPEN (91, FILE = 'ana/deuteron_processes.dat', | |
258 | cms 1 STATUS = 'UNKNOWN') | |
259 | if(idpert.eq.1.or.idpert.eq.2) then | |
260 | cms OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN') | |
261 | endif | |
262 | c.....generate formation time and position at formation time. | |
263 | TAU0 = ARPAR1(1) | |
264 | NP = IAINT2(1) | |
265 | clin-7/10/01 initial positions already given for hadrons | |
266 | c formed from partons inside ZPC (from string melting): | |
267 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then | |
268 | if(NP.le.nattzp) return | |
269 | do 1001 I = nattzp+1, NP | |
270 | IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN | |
271 | PRINT *, ' IN ARINI1' | |
272 | PRINT *, 'ABS(PZ) .GE. EE for particle ', I | |
273 | PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), | |
274 | & ' PY = ', PYAR(I) | |
275 | PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I) | |
276 | PRINT *, ' XM = ', XMAR(I) | |
277 | RAP = 1000000.0 | |
278 | GOTO 50 | |
279 | END IF | |
280 | RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I))) | |
281 | 50 CONTINUE | |
282 | VX = PXAR(I) / PEAR(I) | |
283 | VY = PYAR(I) / PEAR(I) | |
284 | FTAR(I) = TAU0 * COSH(RAP) | |
285 | GXAR(I) = GXAR(I) + VX * FTAR(I) | |
286 | GYAR(I) = GYAR(I) + VY * FTAR(I) | |
287 | GZAR(I) = TAU0 * SINH(RAP) | |
288 | clin-5/2009 No formation time for spectator projectile or target nucleons: | |
289 | if(PXAR(I).eq.0.and.PYAR(I).eq.0 | |
290 | 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99 | |
291 | 2 .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then | |
292 | TAUI=1.E-20 | |
293 | FTAR(I)=TAUI*COSH(RAP) | |
294 | GZAR(I)=TAUI*SINH(RAP) | |
295 | endif | |
296 | 1001 continue | |
297 | clin-7/10/01-end | |
298 | clin-3/2009 cleanup of program flow: | |
299 | else | |
300 | DO 1002 I = 1, NP | |
301 | IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN | |
302 | PRINT *, ' IN ARINI1' | |
303 | PRINT *, 'ABS(PZ) .GE. EE for particle ', I | |
304 | PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), | |
305 | & ' PY = ', PYAR(I) | |
306 | PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I) | |
307 | PRINT *, ' XM = ', XMAR(I) | |
308 | RAP = 1000000.0 | |
309 | GOTO 100 | |
310 | c STOP | |
311 | END IF | |
312 | RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I))) | |
313 | 100 CONTINUE | |
314 | VX = PXAR(I) / PEAR(I) | |
315 | VY = PYAR(I) / PEAR(I) | |
316 | c.....give initial formation time shift | |
317 | TAUI = FTAR(I) + TAU0 | |
318 | FTAR(I) = TAUI * COSH(RAP) | |
319 | GXAR(I) = GXAR(I) + VX * TAU0 * COSH(RAP) | |
320 | GYAR(I) = GYAR(I) + VY * TAU0 * COSH(RAP) | |
321 | c 4/25/03: hadron z-position upon formation determined the same way as x,y: | |
322 | GZAR(I) = TAUI * SINH(RAP) | |
323 | c the old prescription: | |
324 | c GZAR(I) = GZAR(I) + TAU0 * SINH(RAP) | |
325 | zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.) | |
326 | GZAR(I)=GZAR(I)+zsmear | |
327 | cbz1/28/99end | |
328 | c 10/05/01 no formation time for spectator projectile or target nucleons: | |
329 | if(PXAR(I).eq.0.and.PYAR(I).eq.0 | |
330 | 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99 | |
331 | 2 .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then | |
332 | clin-5/2008: | |
333 | c TAUI=0.00001 | |
334 | TAUI=1.E-20 | |
335 | FTAR(I)=TAUI*COSH(RAP) | |
336 | GZAR(I)=TAUI*SINH(RAP)+zsmear | |
337 | endif | |
338 | 1002 CONTINUE | |
339 | clin-3/2009 cleanup of program flow: | |
340 | endif | |
341 | ||
342 | clin-3/2009 Add initial hadrons before the hadron cascade starts: | |
343 | call addhad | |
344 | ||
345 | RETURN | |
346 | END | |
347 | ||
348 | c----------------------------------------------------------------------- | |
349 | ||
350 | c.....subroutine to order particle labels according to increasing | |
351 | c.....formation time | |
352 | ||
353 | SUBROUTINE ARTORD | |
354 | ||
355 | c.....before invoking ARTORD: | |
356 | c.....IAINT2(1) must be set: | |
357 | PARAMETER (MAXSTR=150001,MAXR=1) | |
358 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
359 | cc SAVE /ARPRNT/ | |
360 | COMMON /ARPRC/ ITYPAR(MAXSTR), | |
361 | & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR), | |
362 | & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR), | |
363 | & XMAR(MAXSTR) | |
364 | cc SAVE /ARPRC/ | |
365 | clin-3/2009 Take care of particle weights when user inserts initial hadrons: | |
366 | COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR), | |
367 | 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), | |
368 | 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) | |
369 | DIMENSION dptemp(MAXSTR) | |
370 | c | |
371 | DIMENSION ITYP0(MAXSTR), | |
372 | & GX0(MAXSTR), GY0(MAXSTR), GZ0(MAXSTR), FT0(MAXSTR), | |
373 | & PX0(MAXSTR), PY0(MAXSTR), PZ0(MAXSTR), EE0(MAXSTR), | |
374 | & XM0(MAXSTR) | |
375 | DIMENSION INDX(MAXSTR) | |
376 | EXTERNAL ARINDX | |
377 | SAVE | |
378 | c | |
379 | NPAR = 0 | |
380 | NP = IAINT2(1) | |
381 | DO 1001 I = 1, NP | |
382 | ITYP0(I) = ITYPAR(I) | |
383 | GX0(I) = GXAR(I) | |
384 | GY0(I) = GYAR(I) | |
385 | GZ0(I) = GZAR(I) | |
386 | FT0(I) = FTAR(I) | |
387 | PX0(I) = PXAR(I) | |
388 | PY0(I) = PYAR(I) | |
389 | PZ0(I) = PZAR(I) | |
390 | EE0(I) = PEAR(I) | |
391 | XM0(I) = XMAR(I) | |
392 | clin-3/2009: | |
393 | dptemp(I) = dpertp(I) | |
394 | 1001 CONTINUE | |
395 | CALL ARINDX(MAXSTR, NP, FT0, INDX) | |
396 | DO 1002 I = 1, NP | |
397 | cbz12/3/98 | |
398 | c IF (ITYP0(INDX(I)) .EQ. 211) THEN | |
399 | c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 321) THEN | |
400 | c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR. | |
401 | c & ITYP0(INDX(I)) .EQ. 2112 .OR. ITYP0(INDX(I)) .EQ. -211 .OR. | |
402 | c & ITYP0(INDX(I)) .EQ. 111) THEN | |
403 | c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR. | |
404 | c & ITYP0(INDX(I)) .EQ. 2112) THEN | |
405 | NPAR = NPAR + 1 | |
406 | c ITYPAR(I) = ITYP0(INDX(I)) | |
407 | c GXAR(I) = GX0(INDX(I)) | |
408 | c GYAR(I) = GY0(INDX(I)) | |
409 | c GZAR(I) = GZ0(INDX(I)) | |
410 | c FTAR(I) = FT0(INDX(I)) | |
411 | c PXAR(I) = PX0(INDX(I)) | |
412 | c PYAR(I) = PY0(INDX(I)) | |
413 | c PZAR(I) = PZ0(INDX(I)) | |
414 | c PEAR(I) = EE0(INDX(I)) | |
415 | c XMAR(I) = XM0(INDX(I)) | |
416 | ITYPAR(NPAR) = ITYP0(INDX(I)) | |
417 | GXAR(NPAR) = GX0(INDX(I)) | |
418 | GYAR(NPAR) = GY0(INDX(I)) | |
419 | GZAR(NPAR) = GZ0(INDX(I)) | |
420 | FTAR(NPAR) = FT0(INDX(I)) | |
421 | PXAR(NPAR) = PX0(INDX(I)) | |
422 | PYAR(NPAR) = PY0(INDX(I)) | |
423 | PZAR(NPAR) = PZ0(INDX(I)) | |
424 | PEAR(NPAR) = EE0(INDX(I)) | |
425 | XMAR(NPAR) = XM0(INDX(I)) | |
426 | clin-3/2009: | |
427 | dpertp(NPAR)=dptemp(INDX(I)) | |
428 | c END IF | |
429 | cbz12/3/98end | |
430 | 1002 CONTINUE | |
431 | IAINT2(1) = NPAR | |
432 | c | |
433 | RETURN | |
434 | END | |
435 | ||
436 | c----------------------------------------------------------------------- | |
437 | ||
438 | c.....subroutine to copy individually generated particle record into | |
439 | c.....particle record for many test particle runs. | |
440 | ||
441 | SUBROUTINE ARINI2(K) | |
442 | ||
443 | PARAMETER (MAXSTR=150001,MAXR=1) | |
444 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
445 | cc SAVE /ARPRNT/ | |
446 | COMMON /ARPRC/ ITYPAR(MAXSTR), | |
447 | & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR), | |
448 | & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR), | |
449 | & XMAR(MAXSTR) | |
450 | cc SAVE /ARPRC/ | |
451 | COMMON /ARERC1/MULTI1(MAXR) | |
452 | cc SAVE /ARERC1/ | |
453 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
454 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
455 | & FT1(MAXSTR, MAXR), | |
456 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
457 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
458 | cc SAVE /ARPRC1/ | |
459 | COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR) | |
460 | cc SAVE /tdecay/ | |
461 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
462 | cc SAVE /input1/ | |
463 | COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, | |
464 | & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL | |
465 | cc SAVE /INPUT2/ | |
466 | COMMON/RNDF77/NSEED | |
467 | COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR), | |
468 | 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), | |
469 | 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) | |
470 | cc SAVE /RNDF77/ | |
471 | SAVE | |
472 | ||
473 | MULTI1(K) = IAINT2(1) | |
474 | DO 1001 I = 1, MULTI1(K) | |
475 | ITYP1(I, K) = ITYPAR(I) | |
476 | GX1(I, K) = GXAR(I) | |
477 | GY1(I, K) = GYAR(I) | |
478 | GZ1(I, K) = GZAR(I) | |
479 | FT1(I, K) = FTAR(I) | |
480 | PX1(I, K) = PXAR(I) | |
481 | PY1(I, K) = PYAR(I) | |
482 | PZ1(I, K) = PZAR(I) | |
483 | EE1(I, K) = PEAR(I) | |
484 | XM1(I, K) = XMAR(I) | |
485 | clin-3/2009 hadron weights are initialized in addhad(): | |
486 | clin-5/2008 all hadrons not perturbatively-produced have the weight of 1: | |
487 | c dpp1(I,K)=1. | |
488 | dpp1(I,K)=dpertp(I) | |
489 | 1001 CONTINUE | |
490 | ||
491 | c initialize final time of each particle to ntmax*dt except for | |
492 | c decay daughters, which have values given by tfdcy() and >(ntmax*dt): | |
493 | do 1002 ip=1,MAXSTR | |
494 | tfdcy(ip)=NTMAX*DT | |
495 | tft(ip)=NTMAX*DT | |
496 | 1002 continue | |
497 | c | |
498 | do 1004 irun=1,MAXR | |
499 | do 1003 ip=1,MAXSTR | |
500 | tfdpi(ip,irun)=NTMAX*DT | |
501 | 1003 continue | |
502 | 1004 continue | |
503 | ||
504 | RETURN | |
505 | END | |
506 | ||
507 | c======================================================================= | |
508 | ||
509 | c.....function to convert PDG flavor code into ART flavor code. | |
510 | ||
511 | FUNCTION IARFLV(IPDG) | |
512 | ||
513 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
514 | cc SAVE /input1/ | |
515 | COMMON/RNDF77/NSEED | |
516 | cc SAVE /RNDF77/ | |
517 | SAVE | |
518 | ||
519 | c.....anti-Delta- | |
520 | IF (IPDG .EQ. -1114) THEN | |
521 | IARFLV = -6 | |
522 | RETURN | |
523 | END IF | |
524 | ||
525 | c.....anti-Delta0 | |
526 | IF (IPDG .EQ. -2114) THEN | |
527 | IARFLV = -7 | |
528 | RETURN | |
529 | END IF | |
530 | ||
531 | c.....anti-Delta+ | |
532 | IF (IPDG .EQ. -2214) THEN | |
533 | IARFLV = -8 | |
534 | RETURN | |
535 | END IF | |
536 | ||
537 | c.....anti-Delta++ | |
538 | IF (IPDG .EQ. -2224) THEN | |
539 | IARFLV = -9 | |
540 | RETURN | |
541 | END IF | |
542 | ||
543 | cbzdbg2/23/99 | |
544 | c.....anti-proton | |
545 | IF (IPDG .EQ. -2212) THEN | |
546 | IARFLV = -1 | |
547 | RETURN | |
548 | END IF | |
549 | ||
550 | c.....anti-neutron | |
551 | IF (IPDG .EQ. -2112) THEN | |
552 | IARFLV = -2 | |
553 | RETURN | |
554 | END IF | |
555 | cbzdbg2/23/99end | |
556 | ||
557 | c.....eta | |
558 | IF (IPDG .EQ. 221) THEN | |
559 | IARFLV = 0 | |
560 | RETURN | |
561 | END IF | |
562 | ||
563 | c.....proton | |
564 | IF (IPDG .EQ. 2212) THEN | |
565 | IARFLV = 1 | |
566 | RETURN | |
567 | END IF | |
568 | ||
569 | c.....neutron | |
570 | IF (IPDG .EQ. 2112) THEN | |
571 | IARFLV = 2 | |
572 | RETURN | |
573 | END IF | |
574 | ||
575 | c.....pi- | |
576 | IF (IPDG .EQ. -211) THEN | |
577 | IARFLV = 3 | |
578 | RETURN | |
579 | END IF | |
580 | ||
581 | c.....pi0 | |
582 | IF (IPDG .EQ. 111) THEN | |
583 | IARFLV = 4 | |
584 | RETURN | |
585 | END IF | |
586 | ||
587 | c.....pi+ | |
588 | IF (IPDG .EQ. 211) THEN | |
589 | IARFLV = 5 | |
590 | RETURN | |
591 | END IF | |
592 | ||
593 | c.....Delta- | |
594 | IF (IPDG .EQ. 1114) THEN | |
595 | IARFLV = 6 | |
596 | RETURN | |
597 | END IF | |
598 | ||
599 | c.....Delta0 | |
600 | IF (IPDG .EQ. 2114) THEN | |
601 | IARFLV = 7 | |
602 | RETURN | |
603 | END IF | |
604 | ||
605 | c.....Delta+ | |
606 | IF (IPDG .EQ. 2214) THEN | |
607 | IARFLV = 8 | |
608 | RETURN | |
609 | END IF | |
610 | ||
611 | c.....Delta++ | |
612 | IF (IPDG .EQ. 2224) THEN | |
613 | IARFLV = 9 | |
614 | RETURN | |
615 | END IF | |
616 | ||
617 | c.....Lambda | |
618 | IF (IPDG .EQ. 3122) THEN | |
619 | IARFLV = 14 | |
620 | RETURN | |
621 | END IF | |
622 | ||
623 | c.....Lambda-bar | |
624 | IF (IPDG .EQ. -3122) THEN | |
625 | IARFLV = -14 | |
626 | RETURN | |
627 | END IF | |
628 | ||
629 | c.....Sigma- | |
630 | IF (IPDG .EQ. 3112) THEN | |
631 | IARFLV = 15 | |
632 | RETURN | |
633 | END IF | |
634 | ||
635 | c.....Sigma-bar | |
636 | IF (IPDG .EQ. -3112) THEN | |
637 | IARFLV = -15 | |
638 | RETURN | |
639 | END IF | |
640 | ||
641 | c.....Sigma0 | |
642 | IF (IPDG .EQ. 3212) THEN | |
643 | IARFLV = 16 | |
644 | RETURN | |
645 | END IF | |
646 | ||
647 | c.....Sigma0-bar | |
648 | IF (IPDG .EQ. -3212) THEN | |
649 | IARFLV = -16 | |
650 | RETURN | |
651 | END IF | |
652 | ||
653 | c.....Sigma+ | |
654 | IF (IPDG .EQ. 3222) THEN | |
655 | IARFLV = 17 | |
656 | RETURN | |
657 | END IF | |
658 | ||
659 | c.....Sigma+ -bar | |
660 | IF (IPDG .EQ. -3222) THEN | |
661 | IARFLV = -17 | |
662 | RETURN | |
663 | END IF | |
664 | ||
665 | c.....K- | |
666 | IF (IPDG .EQ. -321) THEN | |
667 | IARFLV = 21 | |
668 | RETURN | |
669 | END IF | |
670 | ||
671 | c.....K+ | |
672 | IF (IPDG .EQ. 321) THEN | |
673 | IARFLV = 23 | |
674 | RETURN | |
675 | END IF | |
676 | ||
677 | c.....temporary entry for K0 | |
678 | IF (IPDG .EQ. 311) THEN | |
679 | IARFLV = 23 | |
680 | RETURN | |
681 | END IF | |
682 | ||
683 | c.....temporary entry for K0bar | |
684 | IF (IPDG .EQ. -311) THEN | |
685 | IARFLV = 21 | |
686 | RETURN | |
687 | END IF | |
688 | ||
689 | c.....temporary entry for K0S and K0L | |
690 | IF (IPDG .EQ. 310 .OR. IPDG .EQ. 130) THEN | |
691 | R = RANART(NSEED) | |
692 | IF (R .GT. 0.5) THEN | |
693 | IARFLV = 23 | |
694 | ELSE | |
695 | IARFLV = 21 | |
696 | END IF | |
697 | RETURN | |
698 | END IF | |
699 | ||
700 | c.....rho- | |
701 | IF (IPDG .EQ. -213) THEN | |
702 | IARFLV = 25 | |
703 | RETURN | |
704 | END IF | |
705 | ||
706 | c.....rho0 | |
707 | IF (IPDG .EQ. 113) THEN | |
708 | IARFLV = 26 | |
709 | RETURN | |
710 | END IF | |
711 | ||
712 | c.....rho+ | |
713 | IF (IPDG .EQ. 213) THEN | |
714 | IARFLV = 27 | |
715 | RETURN | |
716 | END IF | |
717 | ||
718 | c.....omega | |
719 | IF (IPDG .EQ. 223) THEN | |
720 | IARFLV = 28 | |
721 | RETURN | |
722 | END IF | |
723 | ||
724 | c.....phi | |
725 | IF (IPDG .EQ. 333) THEN | |
726 | IARFLV = 29 | |
727 | RETURN | |
728 | END IF | |
729 | ||
730 | c.....K*+ | |
731 | IF (IPDG .EQ. 323) THEN | |
732 | IARFLV = 30 | |
733 | RETURN | |
734 | END IF | |
735 | c.....K*- | |
736 | IF (IPDG .EQ. -323) THEN | |
737 | IARFLV = -30 | |
738 | RETURN | |
739 | END IF | |
740 | c.....temporary entry for K*0 | |
741 | IF (IPDG .EQ. 313) THEN | |
742 | IARFLV = 30 | |
743 | RETURN | |
744 | END IF | |
745 | c.....temporary entry for K*0bar | |
746 | IF (IPDG .EQ. -313) THEN | |
747 | IARFLV = -30 | |
748 | RETURN | |
749 | END IF | |
750 | ||
751 | c...... eta-prime | |
752 | IF (IPDG .EQ. 331) THEN | |
753 | IARFLV = 31 | |
754 | RETURN | |
755 | END IF | |
756 | ||
757 | c...... a1 | |
758 | c IF (IPDG .EQ. 777) THEN | |
759 | c IARFLV = 32 | |
760 | c RETURN | |
761 | c END IF | |
762 | ||
763 | c... cascade- | |
764 | IF (IPDG .EQ. 3312) THEN | |
765 | IARFLV = 40 | |
766 | RETURN | |
767 | END IF | |
768 | ||
769 | c... cascade+ (bar) | |
770 | IF (IPDG .EQ. -3312) THEN | |
771 | IARFLV = -40 | |
772 | RETURN | |
773 | END IF | |
774 | ||
775 | c... cascade0 | |
776 | IF (IPDG .EQ. 3322) THEN | |
777 | IARFLV = 41 | |
778 | RETURN | |
779 | END IF | |
780 | ||
781 | c... cascade0 -bar | |
782 | IF (IPDG .EQ. -3322) THEN | |
783 | IARFLV = -41 | |
784 | RETURN | |
785 | END IF | |
786 | ||
787 | c... Omega- | |
788 | IF (IPDG .EQ. 3334) THEN | |
789 | IARFLV = 45 | |
790 | RETURN | |
791 | END IF | |
792 | ||
793 | c... Omega+ (bar) | |
794 | IF (IPDG .EQ. -3334) THEN | |
795 | IARFLV = -45 | |
796 | RETURN | |
797 | END IF | |
798 | ||
799 | c... Di-Omega | |
800 | IF (IPDG .EQ. 6666) THEN | |
801 | IARFLV = 44 | |
802 | RETURN | |
803 | END IF | |
804 | c sp06/05/01 end | |
805 | ||
806 | clin-3/2009 keep the same ID numbers in case there are initial deuterons: | |
807 | IF (IPDG .EQ. 42 .or. IPDG .EQ. -42) THEN | |
808 | IARFLV = IPDG | |
809 | RETURN | |
810 | END IF | |
811 | ||
812 | c.....other | |
813 | IARFLV = IPDG + 10000 | |
814 | ||
815 | RETURN | |
816 | END | |
817 | ||
818 | c----------------------------------------------------------------------- | |
819 | ||
820 | c.....function to convert ART flavor code into PDG flavor code. | |
821 | ||
822 | FUNCTION INVFLV(IART) | |
823 | ||
824 | common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT | |
825 | cc SAVE /input1/ | |
826 | COMMON/RNDF77/NSEED | |
827 | cc SAVE /RNDF77/ | |
828 | SAVE | |
829 | ||
830 | c.....anti-Delta- | |
831 | IF (IART .EQ. -6) THEN | |
832 | INVFLV = -1114 | |
833 | RETURN | |
834 | END IF | |
835 | ||
836 | c.....anti-Delta0 | |
837 | IF (IART .EQ. -7) THEN | |
838 | INVFLV = -2114 | |
839 | RETURN | |
840 | END IF | |
841 | ||
842 | c.....anti-Delta+ | |
843 | IF (IART .EQ. -8) THEN | |
844 | INVFLV = -2214 | |
845 | RETURN | |
846 | END IF | |
847 | ||
848 | c.....anti-Delta++ | |
849 | IF (IART .EQ. -9) THEN | |
850 | INVFLV = -2224 | |
851 | RETURN | |
852 | END IF | |
853 | ||
854 | cbzdbg2/23/99 | |
855 | c.....anti-proton | |
856 | IF (IART .EQ. -1) THEN | |
857 | INVFLV = -2212 | |
858 | RETURN | |
859 | END IF | |
860 | ||
861 | c.....anti-neutron | |
862 | IF (IART .EQ. -2) THEN | |
863 | INVFLV = -2112 | |
864 | RETURN | |
865 | END IF | |
866 | cbzdbg2/23/99end | |
867 | ||
868 | c.....eta | |
869 | IF (IART .EQ. 0) THEN | |
870 | INVFLV = 221 | |
871 | RETURN | |
872 | END IF | |
873 | ||
874 | c.....proton | |
875 | IF (IART .EQ. 1) THEN | |
876 | INVFLV = 2212 | |
877 | RETURN | |
878 | END IF | |
879 | ||
880 | c.....neutron | |
881 | IF (IART .EQ. 2) THEN | |
882 | INVFLV = 2112 | |
883 | RETURN | |
884 | END IF | |
885 | ||
886 | c.....pi- | |
887 | IF (IART .EQ. 3) THEN | |
888 | INVFLV = -211 | |
889 | RETURN | |
890 | END IF | |
891 | ||
892 | c.....pi0 | |
893 | IF (IART .EQ. 4) THEN | |
894 | INVFLV = 111 | |
895 | RETURN | |
896 | END IF | |
897 | ||
898 | c.....pi+ | |
899 | IF (IART .EQ. 5) THEN | |
900 | INVFLV = 211 | |
901 | RETURN | |
902 | END IF | |
903 | ||
904 | c.....Delta- | |
905 | IF (IART .EQ. 6) THEN | |
906 | INVFLV = 1114 | |
907 | RETURN | |
908 | END IF | |
909 | ||
910 | c.....Delta0 | |
911 | IF (IART .EQ. 7) THEN | |
912 | INVFLV = 2114 | |
913 | RETURN | |
914 | END IF | |
915 | ||
916 | c.....Delta+ | |
917 | IF (IART .EQ. 8) THEN | |
918 | INVFLV = 2214 | |
919 | RETURN | |
920 | END IF | |
921 | ||
922 | c.....Delta++ | |
923 | IF (IART .EQ. 9) THEN | |
924 | INVFLV = 2224 | |
925 | RETURN | |
926 | END IF | |
927 | ||
928 | cc.....N*(1440), N*(1535) temporary entry | |
929 | c IF (IART .GE. 10 .AND. IART .LE.13) THEN | |
930 | c INVFLV = 0 | |
931 | c RETURN | |
932 | c END IF | |
933 | ||
934 | c.....Lambda | |
935 | IF (IART .EQ. 14) THEN | |
936 | INVFLV = 3122 | |
937 | RETURN | |
938 | END IF | |
939 | c.....Lambda-bar | |
940 | IF (IART .EQ. -14) THEN | |
941 | INVFLV = -3122 | |
942 | RETURN | |
943 | END IF | |
944 | ||
945 | cbz3/12/99 | |
946 | c.....temporary entry for Sigma's | |
947 | c IF (IART .EQ. 15) THEN | |
948 | c R = RANART(NSEED) | |
949 | c IF (R .GT. 2. / 3.) THEN | |
950 | c INVFLV = 3112 | |
951 | c ELSE IF (R .GT. 1./ 3. .AND. R .LE. 2. / 3.) THEN | |
952 | c INVFLV = 3212 | |
953 | c ELSE | |
954 | c INVFLV = 3222 | |
955 | c END IF | |
956 | c RETURN | |
957 | c END IF | |
958 | ||
959 | c.....Sigma- | |
960 | IF (IART .EQ. 15) THEN | |
961 | INVFLV = 3112 | |
962 | RETURN | |
963 | END IF | |
964 | ||
965 | c.....Sigma- bar | |
966 | IF (IART .EQ. -15) THEN | |
967 | INVFLV = -3112 | |
968 | RETURN | |
969 | END IF | |
970 | ||
971 | c.....Sigma0 | |
972 | IF (IART .EQ. 16) THEN | |
973 | INVFLV = 3212 | |
974 | RETURN | |
975 | END IF | |
976 | ||
977 | c.....Sigma0 -bar | |
978 | IF (IART .EQ. -16) THEN | |
979 | INVFLV = -3212 | |
980 | RETURN | |
981 | END IF | |
982 | ||
983 | c.....Sigma+ | |
984 | IF (IART .EQ. 17) THEN | |
985 | INVFLV = 3222 | |
986 | RETURN | |
987 | END IF | |
988 | ||
989 | c.....Sigma+ -bar | |
990 | IF (IART .EQ. -17) THEN | |
991 | INVFLV = -3222 | |
992 | RETURN | |
993 | END IF | |
994 | ||
995 | clin-2/23/03 K0S and K0L are generated at the last timestep: | |
996 | c.....temporary entry for K- and K0bar | |
997 | IF (IART .EQ. 21) THEN | |
998 | c R = RANART(NSEED) | |
999 | c IF (R .GT. 0.5) THEN | |
1000 | INVFLV = -321 | |
1001 | c ELSE | |
1002 | c INVFLV = -311 | |
1003 | c R = RANART(NSEED) | |
1004 | c IF (R .GT. 0.5) THEN | |
1005 | c INVFLV = 310 | |
1006 | c ELSE | |
1007 | c INVFLV = 130 | |
1008 | c END IF | |
1009 | c END IF | |
1010 | RETURN | |
1011 | END IF | |
1012 | ||
1013 | c.....temporary entry for K+ and K0 | |
1014 | IF (IART .EQ. 23) THEN | |
1015 | c R = RANART(NSEED) | |
1016 | c IF (R .GT. 0.5) THEN | |
1017 | INVFLV = 321 | |
1018 | c ELSE | |
1019 | c INVFLV = 311 | |
1020 | c R = RANART(NSEED) | |
1021 | c IF (R .GT. 0.5) THEN | |
1022 | c INVFLV = 310 | |
1023 | c ELSE | |
1024 | c INVFLV = 130 | |
1025 | c END IF | |
1026 | c END IF | |
1027 | RETURN | |
1028 | END IF | |
1029 | ||
1030 | c.....K0Long: | |
1031 | IF (IART .EQ. 22) THEN | |
1032 | INVFLV = 130 | |
1033 | RETURN | |
1034 | ENDIF | |
1035 | c.....K0Short: | |
1036 | IF (IART .EQ. 24) THEN | |
1037 | INVFLV = 310 | |
1038 | RETURN | |
1039 | ENDIF | |
1040 | ||
1041 | c.....rho- | |
1042 | IF (IART .EQ. 25) THEN | |
1043 | INVFLV = -213 | |
1044 | RETURN | |
1045 | END IF | |
1046 | ||
1047 | c.....rho0 | |
1048 | IF (IART .EQ. 26) THEN | |
1049 | INVFLV = 113 | |
1050 | RETURN | |
1051 | END IF | |
1052 | ||
1053 | c.....rho+ | |
1054 | IF (IART .EQ. 27) THEN | |
1055 | INVFLV = 213 | |
1056 | RETURN | |
1057 | END IF | |
1058 | ||
1059 | c.....omega | |
1060 | IF (IART .EQ. 28) THEN | |
1061 | INVFLV = 223 | |
1062 | RETURN | |
1063 | END IF | |
1064 | ||
1065 | c.....phi | |
1066 | IF (IART .EQ. 29) THEN | |
1067 | INVFLV = 333 | |
1068 | RETURN | |
1069 | END IF | |
1070 | ||
1071 | c.....temporary entry for K*+ and K*0 | |
1072 | IF (IART .EQ. 30) THEN | |
1073 | INVFLV = 323 | |
1074 | IF (RANART(NSEED).GT.0.5) INVFLV = 313 | |
1075 | RETURN | |
1076 | END IF | |
1077 | ||
1078 | c.....temporary entry for K*- and K*0bar | |
1079 | IF (IART .EQ. -30) THEN | |
1080 | INVFLV = -323 | |
1081 | IF (RANART(NSEED).GT.0.5) INVFLV = -313 | |
1082 | RETURN | |
1083 | END IF | |
1084 | ||
1085 | c... eta-prime (bar) | |
1086 | IF (IART .EQ. 31) THEN | |
1087 | INVFLV = 331 | |
1088 | RETURN | |
1089 | END IF | |
1090 | ||
1091 | c... a1 | |
1092 | IF (IART .EQ. 32) THEN | |
1093 | INVFLV = 777 | |
1094 | RETURN | |
1095 | END IF | |
1096 | ||
1097 | c... cascade- | |
1098 | IF (IART .EQ. 40) THEN | |
1099 | INVFLV = 3312 | |
1100 | RETURN | |
1101 | END IF | |
1102 | ||
1103 | c... cascade+ (bar) | |
1104 | IF (IART .EQ. -40) THEN | |
1105 | INVFLV = -3312 | |
1106 | RETURN | |
1107 | END IF | |
1108 | ||
1109 | c... cascade0 | |
1110 | IF (IART .EQ. 41) THEN | |
1111 | INVFLV = 3322 | |
1112 | RETURN | |
1113 | END IF | |
1114 | ||
1115 | c... cascade0 -bar | |
1116 | IF (IART .EQ. -41) THEN | |
1117 | INVFLV = -3322 | |
1118 | RETURN | |
1119 | END IF | |
1120 | ||
1121 | c... Omega- | |
1122 | IF (IART .EQ. 45) THEN | |
1123 | INVFLV = 3334 | |
1124 | RETURN | |
1125 | END IF | |
1126 | ||
1127 | c... Omega+ (bar) | |
1128 | IF (IART .EQ. -45) THEN | |
1129 | INVFLV = -3334 | |
1130 | RETURN | |
1131 | END IF | |
1132 | ||
1133 | c... Di-Omega | |
1134 | IF (IART .EQ. 44) THEN | |
1135 | INVFLV = 6666 | |
1136 | RETURN | |
1137 | END IF | |
1138 | c sp 12/19/00 end | |
1139 | ||
1140 | clin-5/2008 deuteron ID numbers in ART and ampt.dat: | |
1141 | IF (IART .EQ. 42) THEN | |
1142 | INVFLV = 1000000000+10020 | |
1143 | RETURN | |
1144 | ELSEIF (IART .EQ. -42) THEN | |
1145 | INVFLV = -1000000000-10020 | |
1146 | RETURN | |
1147 | END IF | |
1148 | c | |
1149 | c.....other | |
1150 | INVFLV = IART - 10000 | |
1151 | ||
1152 | RETURN | |
1153 | END | |
1154 | ||
1155 | c======================================================================= | |
1156 | ||
1157 | BLOCK DATA ARDATA | |
1158 | ||
1159 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
1160 | cc SAVE /ARPRNT/ | |
1161 | SAVE | |
1162 | DATA ARPAR1/1.19, 99 * 0.0/ | |
1163 | DATA IAPAR2/3, 49 * 0/ | |
1164 | DATA ARINT1/100 * 0.0/ | |
1165 | DATA IAINT2/50 * 0/ | |
1166 | ||
1167 | END | |
1168 | ||
1169 | c======================================================================= | |
1170 | ||
1171 | c.....Routine borrowed from ZPC. | |
1172 | c.....double precision is modified to real*4. | |
1173 | ||
1174 | cbz1/29/99 | |
1175 | c subroutine index1(n, m, arrin, indx) | |
1176 | subroutine arindx(n, m, arrin, indx) | |
1177 | cbz1/29/99end | |
1178 | c indexes the first m elements of ARRIN of length n, i.e., outputs INDX | |
1179 | c such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m | |
1180 | ||
1181 | c implicit real*4 (a-h, o-z) | |
1182 | ||
1183 | dimension arrin(n), indx(n) | |
1184 | SAVE | |
1185 | do 1001 j = 1, m | |
1186 | indx(j) = j | |
1187 | 1001 continue | |
1188 | l = m / 2 + 1 | |
1189 | ir = m | |
1190 | 10 continue | |
1191 | if (l .gt. 1) then | |
1192 | l = l - 1 | |
1193 | indxt = indx(l) | |
1194 | q = arrin(indxt) | |
1195 | else | |
1196 | indxt = indx(ir) | |
1197 | q = arrin(indxt) | |
1198 | indx(ir) = indx(1) | |
1199 | ir = ir - 1 | |
1200 | if (ir .eq. 1) then | |
1201 | indx(1) = indxt | |
1202 | return | |
1203 | end if | |
1204 | end if | |
1205 | i = l | |
1206 | j = l + l | |
1207 | 20 if (j .le. ir) then | |
1208 | if (j .lt. ir) then | |
1209 | if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1 | |
1210 | end if | |
1211 | if (q .lt. arrin(indx(j))) then | |
1212 | indx(i) = indx(j) | |
1213 | i = j | |
1214 | j = j + j | |
1215 | else | |
1216 | j = ir + 1 | |
1217 | end if | |
1218 | goto 20 | |
1219 | end if | |
1220 | indx(i) = indxt | |
1221 | goto 10 | |
1222 | ||
1223 | end | |
1224 | ||
1225 | c----------------------------------------------------------------------- | |
1226 | ||
1227 | c.....extracted from G. Song's ART expasion including K- interactions | |
1228 | c.....file `NEWKAON.FOR' | |
1229 | ||
1230 | c 5/01/03 send iblock value into art1f.f, necessary for resonance studies: | |
1231 | c subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2, | |
1232 | c & srt,pcx,pcy,pcz) | |
1233 | subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2, | |
1234 | & srt,pcx,pcy,pcz,iblock) | |
1235 | PARAMETER (MAXSTR=150001,MAXR=1) | |
1236 | PARAMETER (AKA=0.498) | |
1237 | COMMON /AA/ R(3,MAXSTR) | |
1238 | cc SAVE /AA/ | |
1239 | COMMON /BB/ P(3,MAXSTR) | |
1240 | cc SAVE /BB/ | |
1241 | COMMON /CC/ E(MAXSTR) | |
1242 | cc SAVE /CC/ | |
1243 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
1244 | cc SAVE /EE/ | |
1245 | COMMON /BG/BETAX,BETAY,BETAZ,GAMMA | |
1246 | cc SAVE /BG/ | |
1247 | COMMON /NN/NNN | |
1248 | cc SAVE /NN/ | |
1249 | COMMON /RUN/NUM | |
1250 | cc SAVE /RUN/ | |
1251 | COMMON /PA/RPION(3,MAXSTR,MAXR) | |
1252 | cc SAVE /PA/ | |
1253 | COMMON /PB/PPION(3,MAXSTR,MAXR) | |
1254 | cc SAVE /PB/ | |
1255 | COMMON /PC/EPION(MAXSTR,MAXR) | |
1256 | cc SAVE /PC/ | |
1257 | COMMON /PD/LPION(MAXSTR,MAXR) | |
1258 | cc SAVE /PD/ | |
1259 | COMMON/RNDF77/NSEED | |
1260 | cc SAVE /RNDF77/ | |
1261 | SAVE | |
1262 | c | |
1263 | logical lb1bn, lb2bn,lb1mn,lb2mn | |
1264 | cbz3/7/99 neutralk | |
1265 | c logical lb1bn1, lb2bayon1, lb1bn0, lb2bn0 | |
1266 | logical lb1bn1, lb2bn1, lb1bn0, lb2bn0 | |
1267 | cbz3/7/99 neutralk end | |
1268 | logical lb1mn0, lb2mn0, lb1mn1, lb2mn1 | |
1269 | logical lb1mn2, lb2mn2 | |
1270 | icase=-1 | |
1271 | c icase: flag for the type of reaction that is going to happen. | |
1272 | c icase=-1, no desired reaction, return to main program. | |
1273 | c 1, NN,ND,DD | |
1274 | c 2, PI+N, PI+D | |
1275 | c 3, K(-) absorption. | |
1276 | nchrg=-100 | |
1277 | c nchrg: Net charges of the two incoming particles. | |
1278 | ictrl = 1 | |
1279 | lb1=lb(i1) | |
1280 | lb2=lb(i2) | |
1281 | em1=e(i1) | |
1282 | em2=e(i2) | |
1283 | lb1bn=lb1.eq.1.or.lb1.eq.2.or.(lb1.gt.5.and.lb1.le.13) | |
1284 | lb2bn=lb2.eq.1.or.lb2.eq.2.or.(lb2.gt.5.and.lb2.le.13) | |
1285 | lb1bn0=lb1.eq.2.or.lb1.eq.7.or.lb1.eq.10.or.lb1.eq.12 | |
1286 | lb2bn0=lb2.eq.2.or.lb2.eq.7.or.lb2.eq.10.or.lb2.eq.12 | |
1287 | lb1bn1=lb1.eq.1.or.lb1.eq.8.or.lb1.eq.11.or.lb1.eq.13 | |
1288 | lb2bn1=lb2.eq.1.or.lb2.eq.8.or.lb2.eq.11.or.lb2.eq.13 | |
1289 | lb1mn=em1.lt.0.2.or.lb1.eq.0.or.(lb1.ge.25.and.lb1.le.29) | |
1290 | lb2mn=em2.lt.0.2.or.lb2.eq.0.or.(lb2.ge.25.and.lb2.le.29) | |
1291 | lb1mn0=lb1.eq.0.or.lb1.eq.4.or.lb1.eq.26.or. | |
1292 | & lb1.eq.28.or.lb1.eq.29 | |
1293 | lb2mn0=lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or. | |
1294 | & lb2.eq.28.or.lb2.eq.29 | |
1295 | lb1mn1= lb1.eq.5.or.lb1.eq.27 | |
1296 | lb2mn1= lb2.eq.5.or.lb2.eq.27 | |
1297 | lb1mn2=lb1.eq.3.or.lb1.eq.25 | |
1298 | lb2mn2=lb2.eq.3.or.lb2.eq.25 | |
1299 | ||
1300 | c 1. consider N+N, N+Resonance, R + R reactions | |
1301 | if(lb1bn.and.lb2bn) then | |
1302 | c NN,ND,DD: | |
1303 | icase=1 | |
1304 | c total cross section | |
1305 | sig=40. | |
1306 | if(lb1.eq.9.and.lb2.eq.9) then | |
1307 | nchrg=4 | |
1308 | endif | |
1309 | if((lb1bn1.and.lb2.eq.9) | |
1310 | & .or.(lb2bn1.and.lb1.eq.9))then | |
1311 | nchrg=3 | |
1312 | endif | |
1313 | if((lb1bn0.and.lb2.eq.9) | |
1314 | & .or.(lb2bn0.and.lb1.eq.9) | |
1315 | & .or.(lb1bn1.and.lb2bn1)) then | |
1316 | nchrg=2 | |
1317 | endif | |
1318 | if((lb1bn1.and.lb2bn0).or.(lb1.eq.6.and.lb2.eq.9) | |
1319 | & .or.(lb2bn1.and.lb1bn0) | |
1320 | & .or.(lb2.eq.6.and.lb1.eq.9))then | |
1321 | nchrg=1 | |
1322 | endif | |
1323 | if((lb1bn0.and.lb2bn0).or.(lb1bn1.and.lb2.eq.6) | |
1324 | & .or.(lb2bn1.and.lb1.eq.6)) then | |
1325 | nchrg=0 | |
1326 | endif | |
1327 | if((lb1bn0.and.lb2.eq.6) | |
1328 | & .or.(lb2bn0.and.lb1.eq.6))then | |
1329 | nchrg=-1 | |
1330 | endif | |
1331 | if(lb1.eq.6.and.lb2.eq.6) then | |
1332 | nchrg=-2 | |
1333 | endif | |
1334 | c brsig = x2kaon_no_isospin(srt) | |
1335 | if(nchrg.ge.-1.and.nchrg.le.2) then | |
1336 | c K,Kbar prduction x sect. | |
1337 | brsig = x2kaon(srt) | |
1338 | else | |
1339 | brsig=0.0 | |
1340 | c if(nchrg.eq.-2.or.nchrg.eq.3) then | |
1341 | c brsig = x2kaon(srt+0.938-1.232) | |
1342 | c else | |
1343 | c nchrg=4 | |
1344 | c brsig = x2kaon(srt+2.*(0.938-1.232)) | |
1345 | c endif | |
1346 | endif | |
1347 | ||
1348 | cbz3/7/99 neutralk | |
1349 | BRSIG = 2.0 * BRSIG | |
1350 | cbz3/7/99 neutralk end | |
1351 | ||
1352 | endif | |
1353 | ||
1354 | c 2. consider PI(meson:eta,omega,rho,phi) + N(N*,D) | |
1355 | if((lb1bn.and.lb2mn).OR.(lb2bn.and.lb1mn)) then | |
1356 | c PN,PD | |
1357 | icase=2 | |
1358 | sig=20. | |
1359 | sigma0 = piNsg0(srt) | |
1360 | brsig=0.0 | |
1361 | if((lb1bn1.and.lb2mn0) | |
1362 | & .or.(lb2bn1.and.lb1mn0). | |
1363 | & or.(lb1bn0.and.lb2mn1).or.(lb2bn0.and.lb1mn1). | |
1364 | & or.(lb1.eq.9.and.lb2mn2).or.(lb2.eq.9.and.lb1mn2))then | |
1365 | nchrg=1 | |
1366 | cbz3/2/99/song | |
1367 | c if(lb1bn1.or.lb2bn1) brsig=2.0*sigma0 | |
1368 | c if(lb1bn0.or.lb2bn0) brsig=0.5*sigma0 | |
1369 | if(lb1bn1.or.lb2bn1) brsig=0.5*sigma0 | |
1370 | if(lb1bn0.or.lb2bn0) brsig=2.0*sigma0 | |
1371 | cbz3/2/99/song end | |
1372 | c if(lb1.eq.9.or.lb2.eq.9) brsig=1.5*sigma0 | |
1373 | endif | |
1374 | if( (lb1bn0.and.lb2mn0 ) | |
1375 | & .or.(lb2bn0.and.lb1mn0) | |
1376 | & .or.(lb1bn1.and.lb2mn2).or.(lb2bn1.and.lb1mn2) | |
1377 | & .or.(lb1.eq.6.and.lb2mn1).or.(lb2.eq.6.and.lb1mn1)) then | |
1378 | nchrg=0 | |
1379 | if(lb1bn1.or.lb2bn1) then | |
1380 | cbz3/2/99/song | |
1381 | c brsig=1.5*sigma0 | |
1382 | brsig=3.0*sigma0 | |
1383 | cbz3/2/99/song end | |
1384 | cbz3/11/99/song | |
1385 | c ratiok = 1./3. | |
1386 | ratiok = 2./3. | |
1387 | cbz3/11/99/song end | |
1388 | ||
1389 | c ratiok: the ratio of channels: ->nK+k- vs. -> pK0K- | |
1390 | endif | |
1391 | if(lb1bn0.or.lb2bn0) then | |
1392 | brsig=2.5*sigma0 | |
1393 | cbz3/2/99/song | |
1394 | c ratiok = 0.8 | |
1395 | ratiok = 0.2 | |
1396 | cbz3/2/99/song end | |
1397 | endif | |
1398 | c if(lb1.eq.6.or.lb2.eq.6) then | |
1399 | c lb=6 : D- | |
1400 | c brsig=1.5*sigma0 | |
1401 | c ratiok = 0.5 | |
1402 | c endif | |
1403 | endif | |
1404 | if( (lb1bn0.and.lb2mn2) | |
1405 | & .or.(lb2bn0.and.lb1mn2) | |
1406 | & .or.(lb1.eq.6.and.lb2mn0).or.(lb2.eq.6.and.lb1mn0)) then | |
1407 | nchrg=-1 | |
1408 | if(lb1bn0.or.lb2bn0) brsig=sigma0 | |
1409 | c if(lb1.eq.6.or.lb2.eq.6) brsig=sigma0 | |
1410 | endif | |
1411 | c if((lb1.eq.6.and.lb2mn2).or.(lb2.eq.6.and.lb1mn2))then | |
1412 | c nchrg=-2 | |
1413 | c endif | |
1414 | c if((lb1bn1.and.lb2mn1).or.(lb2bn1.and.lb1mn1) | |
1415 | c & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then | |
1416 | c nchrg=2 | |
1417 | c endif | |
1418 | ||
1419 | cbz3/11/99 neutralk | |
1420 | if((lb1.eq.6.and.lb2mn2) | |
1421 | & .or.(lb2.eq.6.and.lb1mn2))then | |
1422 | nchrg=-2 | |
1423 | endif | |
1424 | cbz3/11/99 neutralk | |
1425 | cbz3/8/99 neutralk | |
1426 | if((lb1bn1.and.lb2mn1) | |
1427 | & .or.(lb2bn1.and.lb1mn1) | |
1428 | & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then | |
1429 | nchrg=2 | |
1430 | endif | |
1431 | cbz3/8/99 neutralk end | |
1432 | ||
1433 | cbz3/7/99 neutralk | |
1434 | IF (NCHRG .GE. -2 .AND. NCHRG .LE. 2) THEN | |
1435 | BRSIG = 3.0 * SIGMA0 | |
1436 | END IF | |
1437 | cbz3/7/99 neutralk end | |
1438 | ||
1439 | endif | |
1440 | ||
1441 | c 3. consider K- + N(N*,D) absorption. | |
1442 | c if((lb1bn.and.lb2.eq.21).OR.(lb2bn.and.lb1.eq.21)) then | |
1443 | if( (lb1bn.and.(lb2.eq.21.or.lb2.eq.-30)).OR. | |
1444 | & (lb2bn.and.(lb1.eq.21.or.lb1.eq.-30)) )then | |
1445 | c bmass=em1+em2-aka | |
1446 | bmass=0.938 | |
1447 | if(srt.le.(bmass+aka)) then | |
1448 | cbz3/2/99 | |
1449 | c write(100,*)'--lb1,lb2,em1,em2,srt',lb1,lb2,em1,em2,srt | |
1450 | cbz3/2/99end | |
1451 | pkaon=0. | |
1452 | else | |
1453 | pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2) | |
1454 | endif | |
1455 | sig=0. | |
1456 | if(lb1.eq.1.or.lb2.eq.1.or.lb1.eq.8.or.lb2.eq.8.or. | |
1457 | & lb1.eq.11.or.lb2.eq.11.or.lb1.eq.13.or.lb2.eq.13) then | |
1458 | c K- + (D+,N*+)p -> | |
1459 | nchrg=0 | |
1460 | sigela=akPel(pkaon) | |
1461 | sigsgm=3.*akPsgm(pkaon) | |
1462 | sig=sigela+sigsgm+akPlam(pkaon) | |
1463 | endif | |
1464 | if(lb1.eq.2.or.lb2.eq.2.or.lb1.eq.7.or.lb2.eq.7.or. | |
1465 | & lb1.eq.10.or.lb2.eq.10.or.lb1.eq.12.or.lb2.eq.12) then | |
1466 | c K- + (D0, N*0)n -> | |
1467 | nchrg=-1 | |
1468 | sigela=akNel(pkaon) | |
1469 | sigsgm=2.*akNsgm(pkaon) | |
1470 | sig=sigela+sigsgm+akNlam(pkaon) | |
1471 | endif | |
1472 | if(lb1.eq.6.or.lb2.eq.6) then | |
1473 | c K- + D- | |
1474 | nchrg=-2 | |
1475 | sigela=akNel(pkaon) | |
1476 | sigsgm=akNsgm(pkaon) | |
1477 | sig=sigela+sigsgm | |
1478 | endif | |
1479 | if(lb1.eq.9.or.lb2.eq.9) then | |
1480 | c K- + D++ | |
1481 | nchrg=1 | |
1482 | sigela=akPel(pkaon) | |
1483 | sigsgm=2.*akPsgm(pkaon) | |
1484 | sig=sigela+sigsgm+akPlam(pkaon) | |
1485 | endif | |
1486 | ||
1487 | cbz3/8/99 neutralk | |
1488 | sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON)) | |
1489 | SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON) | |
1490 | SIG = sigela + SIGSGM + AKPLAM(PKAON) | |
1491 | cbz3/8/99 neutralk end | |
1492 | ||
1493 | if(sig.gt.1.e-7) then | |
1494 | c K(-) + N reactions | |
1495 | icase=3 | |
1496 | brel=sigela/sig | |
1497 | brsgm=sigsgm/sig | |
1498 | c branch_lambda=akNlam(pkaon)/sig | |
1499 | brsig = sig | |
1500 | endif | |
1501 | endif | |
1502 | ||
1503 | c 4. meson + hyperon -> K- + N | |
1504 | c if(((lb1.ge.14.and.lb1.le.17).and.lb2mn).OR. | |
1505 | c & ((lb2.ge.14.and.lb2.le.17).and.lb1mn)) then | |
1506 | if(((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.3.and.lb2.le.5)).OR. | |
1507 | & ((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.3.and.lb1.le.5)))then | |
1508 | c first classify the reactions due to total charge. | |
1509 | nchrg=-100 | |
1510 | if((lb1.eq.15.and.(lb2.eq.3.or.lb2.eq.25)).OR. | |
1511 | & (lb2.eq.15.and.(lb1.eq.3.or.lb1.eq.25))) then | |
1512 | nchrg=-2 | |
1513 | c D- | |
1514 | bmass=1.232 | |
1515 | endif | |
1516 | if((lb1.eq.15.and.lb2mn0).or.(lb2.eq.15.and.lb1mn0).OR. | |
1517 | & ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.3.or.lb2.eq.25)).OR. | |
1518 | & ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.3.or.lb1.eq.25)))then | |
1519 | nchrg=-1 | |
1520 | c n | |
1521 | bmass=0.938 | |
1522 | endif | |
1523 | if((lb1.eq.15.and.(lb2.eq.5.or.lb2.eq.27)).OR. | |
1524 | & (lb2.eq.15.and.(lb1.eq.5.or.lb1.eq.27)).or. | |
1525 | & (lb1.eq.17.and.(lb2.eq.3.or.lb2.eq.25)).OR. | |
1526 | & (lb2.eq.17.and.(lb1.eq.3.or.lb1.eq.25)).or. | |
1527 | & ((lb1.eq.14.or.lb1.eq.16).and.lb2mn0).OR. | |
1528 | & ((lb2.eq.14.or.lb2.eq.16).and.lb1mn0)) then | |
1529 | nchrg=0 | |
1530 | c p | |
1531 | bmass=0.938 | |
1532 | endif | |
1533 | if((lb1.eq.17.and.lb2mn0).or.(lb2.eq.17.and.lb1mn0).OR. | |
1534 | & ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.5.or.lb2.eq.27)).OR. | |
1535 | & ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.5.or.lb1.eq.27)))then | |
1536 | nchrg=1 | |
1537 | c D++ | |
1538 | bmass=1.232 | |
1539 | endif | |
1540 | sig = 0. | |
1541 | if(nchrg.ne.-100.and.srt.gt.(aka+bmass)) then | |
1542 | c PI+sigma or PI + Lambda => Kbar + N reactions | |
1543 | icase=4 | |
1544 | c pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2) | |
1545 | pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2) | |
1546 | c lambda + Pi | |
1547 | if(lb1.eq.14.or.lb2.eq.14) then | |
1548 | if(nchrg.ge.0) sigma0=akPlam(pkaon) | |
1549 | if(nchrg.lt.0) sigma0=akNlam(pkaon) | |
1550 | c sigma + pi | |
1551 | else | |
1552 | c K-p or K-D++ | |
1553 | if(nchrg.ge.0) sigma0=akPsgm(pkaon) | |
1554 | c K-n or K-D- | |
1555 | if(nchrg.lt.0) sigma0=akNsgm(pkaon) | |
1556 | ||
1557 | cbz3/8/99 neutralk | |
1558 | SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON) | |
1559 | cbz3/8/99 neutralk end | |
1560 | ||
1561 | endif | |
1562 | sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/ | |
1563 | & (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0 | |
1564 | cbz3/8/99 neutralk | |
1565 | c if(nchrg.eq.-2.or.nchrg.eq.1) sig=2.*sig K-D++, K-D- | |
1566 | c K0barD++, K-D- | |
1567 | if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig | |
1568 | ||
1569 | cbz3/8/99 neutralk end | |
1570 | ||
1571 | c the factor 2 comes from spin of delta, which is 3/2 | |
1572 | c detailed balance. copy from Page 423 of N.P. A614 1997 | |
1573 | ||
1574 | cbz3/8/99 neutralk | |
1575 | IF (LB1 .EQ. 14 .OR. LB2 .EQ. 14) THEN | |
1576 | SIG = 4.0 / 3.0 * SIG | |
1577 | ELSE IF (NCHRG .EQ. -2 .OR. NCHRG .EQ. 2) THEN | |
1578 | SIG = 8.0 / 9.0 * SIG | |
1579 | ELSE | |
1580 | SIG = 4.0 / 9.0 * SIG | |
1581 | END IF | |
1582 | cbz3/8/99 neutralk end | |
1583 | brsig = sig | |
1584 | if(sig.lt.1.e-7) sig = 1.e-7 | |
1585 | endif | |
1586 | csp05/07/01 | |
1587 | * comment icase=4 statement below if only inelastic | |
1588 | c PI+L/Si => Kbar + N OR ELASTIC SCATTERING | |
1589 | icase=4 | |
1590 | brsig = sig | |
1591 | c elastic xsecn of 10mb | |
1592 | sigela = 10. | |
1593 | sig = sig + sigela | |
1594 | brel = sigela/sig | |
1595 | cc brsig = sig | |
1596 | csp05/07/01 end | |
1597 | endif | |
1598 | c | |
1599 | c if(em2.lt.0.2.and.em1.lt.0.2) then | |
1600 | c PI + PI | |
1601 | c icase=5 | |
1602 | c assumed PI PI total x section. | |
1603 | c sig=50. | |
1604 | c Mk + Mkbar | |
1605 | c s0=aka+aka | |
1606 | c brsig = 0. | |
1607 | c if(srt.gt.s0) brsig = 2.7*(1.-s0**2/srt**2)**0.76 | |
1608 | c x section for PIPI->KKbar PRC43 (1991) 1881 | |
1609 | c endif | |
1610 | if(icase.eq.-1) then | |
1611 | ictrl = -1 | |
1612 | return | |
1613 | endif | |
1614 | px1cm=pcx | |
1615 | py1cm=pcy | |
1616 | pz1cm=pcz | |
1617 | ds=sqrt(sig/31.4) | |
1618 | dsr=ds+0.1 | |
1619 | ec=(em1+em2+0.02)**2 | |
1620 | c ec=3.59709 | |
1621 | c if((e(i1).ge.1.).and.(e(i2).ge.1.)) ec = 4.75 | |
1622 | ||
1623 | call distce(i1,i2,dsr,ds,dt,ec,srt,ic,px1cm,py1cm,pz1cm) | |
1624 | if(ic.eq.-1) then | |
1625 | c no anti-kaon production | |
1626 | ictrl = -1 | |
1627 | c in=in+1 | |
1628 | c write(60,*)'--------------distance-----',in | |
1629 | return | |
1630 | endif | |
1631 | ||
1632 | clin-10/24/02 set to 0: ik,ik0-3,il,im,im3-4,in,inpion,ipipi, | |
1633 | c sgsum,sgsum1,sgsum3: | |
1634 | ik=0 | |
1635 | ik0=0 | |
1636 | ik1=0 | |
1637 | ik2=0 | |
1638 | ik3=0 | |
1639 | il=0 | |
1640 | im=0 | |
1641 | im3=0 | |
1642 | im4=0 | |
1643 | in=0 | |
1644 | inpion=0 | |
1645 | ipipi=0 | |
1646 | sgsum=0. | |
1647 | sgsum1=0. | |
1648 | sgsum3=0. | |
1649 | if(icase.eq.1) then | |
1650 | ik=ik+1 | |
1651 | if(srt.gt.2.8639) then | |
1652 | ik0=ik0+1 | |
1653 | if(em1.lt.1.0.and.em2.lt.1.0) then | |
1654 | ik1=ik1+1 | |
1655 | sgsum1=sgsum1+brsig | |
1656 | c ratio_1=sgsum1/ik1/40. | |
1657 | endif | |
1658 | if(em1.gt.1.0.and.em2.gt.1.0) then | |
1659 | ik3=ik3+1 | |
1660 | sgsum3=sgsum3+brsig | |
1661 | c ratio_3=sgsum3/ik3/40. | |
1662 | endif | |
1663 | if(em1.gt.1.0.and.em2.lt.1.0) ik2=ik2+1 | |
1664 | if(em1.lt.1.0.and.em2.gt.1.0) ik2=ik2+1 | |
1665 | sgsum=sgsum+brsig | |
1666 | c ratio=sgsum/ik0/40. | |
1667 | endif | |
1668 | endif | |
1669 | if(icase.eq.2) inpion=inpion+1 | |
1670 | if(icase.eq.5) ipipi=ipipi+1 | |
1671 | c write(62,*)'ik1,ik2,ik3',ik1,ik2,ik3,ratio_1,ratio_3,ratio | |
1672 | c write(62,*)'inpion,ipipi',inpion,ipipi | |
1673 | if(RANART(NSEED).gt.(brsig/sig)) then | |
1674 | c no anti-kaon production | |
1675 | ictrl = -1 | |
1676 | return | |
1677 | endif | |
1678 | il=il+1 | |
1679 | c kaons could be created now. | |
1680 | if(icase.eq.1) then | |
1681 | in=in+1 | |
1682 | c write(60,*)'------in,s2kaon,sig=',in,brsig,sig,lb1,lb2 | |
1683 | call nnkaon(irun,iseed, | |
1684 | & ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg) | |
1685 | endif | |
1686 | if(icase.eq.2) then | |
1687 | im=im+1 | |
1688 | c call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt, | |
1689 | c & pcx,pcy,pcz,nchrg,ratiok) | |
1690 | call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt, | |
1691 | & pcx,pcy,pcz,nchrg,ratiok,iblock) | |
1692 | endif | |
1693 | c | |
1694 | if(icase.eq.3) then | |
1695 | im3=im3+1 | |
1696 | c write(63,*)'im3,lb1,lb2,pkaon',im3,lb1,lb2,pkaon | |
1697 | c write(63,*)'sig,el,sigma',sig,brel,brsgm | |
1698 | c write(63,*)'srt,pcx,pcy,pcz,em1,em2',srt,pcx,pcy,pcz,em1,em2 | |
1699 | call kaonN(brel,brsgm,irun,iseed,dt,nt,ictrl, | |
1700 | & i1,i2,iblock,srt,pcx,pcy,pcz,nchrg) | |
1701 | c this subroutine format is diff. since three final states are possible | |
1702 | endif | |
1703 | c | |
1704 | ||
1705 | if(icase.eq.4) then | |
1706 | im4=im4+1 | |
1707 | c write(64,*)'im4,sigma0,branch,sig=',im4,sigma0,brsig,sig | |
1708 | c write(64,*)'lb1,lb2,em1,em2,pkaon=',lb1,lb2,em1,em2,pkaon | |
1709 | ||
1710 | csp06/07/01 | |
1711 | if(RANART(NSEED).lt.brel) then | |
1712 | ielstc = 1 | |
1713 | else | |
1714 | ielstc = 0 | |
1715 | endif | |
1716 | c call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt, | |
1717 | c & pcx,pcy,pcz,nchrg) | |
1718 | call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt, | |
1719 | & pcx,pcy,pcz,nchrg,iblock) | |
1720 | ||
1721 | csp06/07/01 end | |
1722 | endif | |
1723 | c if(icase.eq.5) then | |
1724 | c im5=im5+1 | |
1725 | c write(65,*)'---im5,s2kaon,sig=',im5,brsig,sig | |
1726 | c call pipikaon(irun,iseed,dt,nt,ictrl,i1,i2,srt,pcx,pcy,pcz) | |
1727 | c endif | |
1728 | cbz3/2/99 | |
1729 | c write(101,*)lb1,lb2,lb(i1),lb(i2) | |
1730 | c write(101,*)em1,em2,e(i1),e(i2),srt | |
1731 | cbz3/2/99end | |
1732 | ||
1733 | return | |
1734 | end | |
1735 | ||
1736 | ****************************************** | |
1737 | * for pp-->pp + kaon + anti-kaon | |
1738 | c real*4 function X2kaon(srt) | |
1739 | real function X2kaon(srt) | |
1740 | SAVE | |
1741 | * This function contains the experimental total pp->pp+K(+)K(-) Xsections * | |
1742 | * srt = DSQRT(s) in GeV * | |
1743 | * xsec = production cross section in mb * | |
1744 | * * | |
1745 | ****************************************** | |
1746 | c minimum c.m.s. energy to create 2 kaon. = 2*(mp+mk) | |
1747 | smin = 2.8639 | |
1748 | x2kaon=0.0000001 | |
1749 | if(srt.lt.smin)return | |
1750 | sigma1 = 2.8 | |
1751 | sigma2 = 7.7 | |
1752 | sigma3 = 3.9 | |
1753 | x = srt**2/smin**2 + 0.0000001 | |
1754 | f1 = (1.+1./sqrt(x))*alog(x) - 4.*(1.-1./sqrt(x)) | |
1755 | f2 = 1. - (1./sqrt(x))*(1.+alog(sqrt(x))) | |
1756 | f3 = ((x-1.)/x**2)**3.5 | |
1757 | x2kaon = (1.-1./x)**3*(sigma1*f1 + sigma2*f2) + sigma3*f3 | |
1758 | return | |
1759 | END | |
1760 | ||
1761 | real function piNsg0(srt) | |
1762 | SAVE | |
1763 | * cross section in mb for PI- + P -> P + K0 + K- | |
1764 | c Mn + 2* Mk | |
1765 | srt0 = 0.938 + 2.*0.498 | |
1766 | if(srt.lt.srt0) then | |
1767 | piNsg0 = 0.0 | |
1768 | return | |
1769 | endif | |
1770 | ratio = srt0**2/srt**2 | |
1771 | piNsg0=1.121*(1.-ratio)**1.86*ratio**2 | |
1772 | return | |
1773 | end | |
1774 | ||
1775 | real function akNel(pkaon) | |
1776 | SAVE | |
1777 | *cross section in mb for K- + N reactions. | |
1778 | c the following data come from PRC 41 (1701) | |
1779 | c sigma1: K(-) + neutron elastic | |
1780 | if(pkaon.lt.0.5.or. pkaon.ge.4.0) sigma1=0. | |
1781 | if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma1=20.*pkaon**2.74 | |
1782 | if(pkaon.ge.1.0.and.pkaon.lt.4.0) sigma1=20.*pkaon**(-1.8) | |
1783 | akNel=sigma1 | |
1784 | return | |
1785 | end | |
1786 | ||
1787 | real function akPel(pkaon) | |
1788 | SAVE | |
1789 | *cross section in mb for K- + N reactions. | |
1790 | c the following data come from PRC 41 (1701) | |
1791 | c sigma2: K(-) + proton elastic | |
1792 | if(pkaon.lt.0.25.or. pkaon.ge.4.0) sigma2=0. | |
1793 | if(pkaon.ge.0.25.and.pkaon.lt.4.0) sigma2=13.*pkaon**(-0.9) | |
1794 | akPel=sigma2 | |
1795 | return | |
1796 | end | |
1797 | ||
1798 | real function akNsgm(pkaon) | |
1799 | SAVE | |
1800 | *cross section in mb for K- + N reactions. | |
1801 | c sigma2: x section for K- + n -> sigma0 + PI- | |
1802 | if(pkaon.lt.0.5.or. pkaon.ge.6.0) sigma2=0. | |
1803 | if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma2=1.2*pkaon**(-1.3) | |
1804 | if(pkaon.ge.1.0.and.pkaon.lt.6.0) sigma2=1.2*pkaon**(-2.3) | |
1805 | akNsgm=sigma2 | |
1806 | return | |
1807 | end | |
1808 | ||
1809 | real function akPsgm(pkaon) | |
1810 | SAVE | |
1811 | *cross section in mb for K- + N reactions. | |
1812 | c sigma1: x section for K- + p -> sigma0 + PI0 | |
1813 | if(pkaon.lt.0.2.or. pkaon.ge.1.5) sigma1=0. | |
1814 | if(pkaon.ge.0.2.and.pkaon.lt.1.5) sigma1=0.6*pkaon**(-1.8) | |
1815 | akPsgm=sigma1 | |
1816 | return | |
1817 | end | |
1818 | ||
1819 | real function akPlam(pkaon) | |
1820 | SAVE | |
1821 | *cross section in mb for K- + N reactions. | |
1822 | c sigma: x section for K- + p -> lambda + PI0 | |
1823 | p=pkaon | |
1824 | if(pkaon.lt.0.2.or. pkaon.ge.10.0) sigma=0. | |
1825 | if(pkaon.ge.0.2.and.pkaon.lt.0.9) sigma=50.*p**2-67.*p+24. | |
1826 | if(pkaon.ge.0.9.and.pkaon.lt.10.0) sigma=3.0*pkaon**(-2.6) | |
1827 | akPlam=sigma | |
1828 | return | |
1829 | end | |
1830 | ||
1831 | real function akNlam(pkaon) | |
1832 | SAVE | |
1833 | *cross section in mb for K- + N reactions. | |
1834 | akNlam=akPlam(pkaon) | |
1835 | return | |
1836 | end | |
1837 | ||
1838 | * GQ Li parametrization (without resonance) | |
1839 | real function akNPsg(pkaon) | |
1840 | SAVE | |
1841 | *cross section in mb for K- + N reactions. | |
1842 | c sigma1: x section for K- + p/n -> sigma0 + PI0 | |
1843 | if(pkaon.le.0.345)then | |
1844 | sigma1=0.624*pkaon**(-1.83) | |
1845 | else | |
1846 | sigma1=0.7*pkaon**(-2.09) | |
1847 | endif | |
1848 | akNPsg=sigma1 | |
1849 | return | |
1850 | end | |
1851 | ||
1852 | c----------------------------------------------------------------------- | |
1853 | ||
1854 | c.....extracted from G. Song's ART expasion including K- interactions | |
1855 | c.....file `NEWNNK.FOR' | |
1856 | ||
1857 | subroutine nnkaon(irun,iseed,ictrl,i1,i2,iblock, | |
1858 | & srt,pcx,pcy,pcz,nchrg) | |
1859 | c <pt>=0.27+0.037*log(srt) was changed to 0.632 + ... on Aug. 14, 1997 | |
1860 | c CANCELED also alpha=1 changed to alpha=3 to decrease the leadng effect. | |
1861 | PARAMETER (MAXSTR=150001,MAXR=1) | |
1862 | PARAMETER (AKA=0.498) | |
1863 | COMMON /AA/ R(3,MAXSTR) | |
1864 | cc SAVE /AA/ | |
1865 | COMMON /BB/ P(3,MAXSTR) | |
1866 | cc SAVE /BB/ | |
1867 | COMMON /CC/ E(MAXSTR) | |
1868 | cc SAVE /CC/ | |
1869 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
1870 | cc SAVE /EE/ | |
1871 | COMMON /BG/BETAX,BETAY,BETAZ,GAMMA | |
1872 | cc SAVE /BG/ | |
1873 | COMMON /NN/NNN | |
1874 | cc SAVE /NN/ | |
1875 | COMMON /RUN/NUM | |
1876 | cc SAVE /RUN/ | |
1877 | COMMON /PA/RPION(3,MAXSTR,MAXR) | |
1878 | cc SAVE /PA/ | |
1879 | COMMON /PB/PPION(3,MAXSTR,MAXR) | |
1880 | cc SAVE /PB/ | |
1881 | COMMON /PC/EPION(MAXSTR,MAXR) | |
1882 | cc SAVE /PC/ | |
1883 | COMMON /PD/LPION(MAXSTR,MAXR) | |
1884 | cc SAVE /PD/ | |
1885 | dimension px(4),py(4),pz(4) | |
1886 | COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR), | |
1887 | 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), | |
1888 | 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) | |
1889 | SAVE | |
1890 | c dm1=e(i1) | |
1891 | c dm2=e(i2) | |
1892 | dm3=0.938 | |
1893 | dm4=0.938 | |
1894 | c 10/24/02 initialize n to 0: | |
1895 | n=0 | |
1896 | ||
1897 | cbz3/11/99 neutralk | |
1898 | c if(nchrg.eq.-2.or.nchrg.ge.3) dm3=1.232 | |
1899 | c if(nchrg.eq.4) dm4=1.232 | |
1900 | if(nchrg.le.-1.or.nchrg.ge.3) dm3=1.232 | |
1901 | if(nchrg.eq.-2.or.nchrg.eq.4) dm4=1.232 | |
1902 | cbz3/11/99 neutralk end | |
1903 | iblock = 0 | |
1904 | call fstate(iseed,srt,dm3,dm4,px,py,pz,iflag) | |
1905 | if(iflag.lt.0) then | |
1906 | c write(60,*)'------------final state fail-------',n | |
1907 | c no anti-kaon production | |
1908 | ictrl = -1 | |
1909 | n=n+1 | |
1910 | return | |
1911 | endif | |
1912 | iblock = 12 | |
1913 | * Rotate the momenta of particles in the cms of I1 & I2 | |
1914 | * px(1), py(1), pz(1): momentum of I1 | |
1915 | * px(2), py(2), pz(2): momentum of I2 | |
1916 | * px(3), py(3), pz(3): momentum of anti-kaon | |
1917 | * px(4), py(4), pz(4): momentum of kaon | |
1918 | ||
1919 | ||
1920 | c 10/28/02 get rid of argument usage mismatch in rotate(): | |
1921 | pxrota=px(1) | |
1922 | pyrota=py(1) | |
1923 | pzrota=pz(1) | |
1924 | c call rotate(pcx,pcy,pcz,px(1),py(1),pz(1)) | |
1925 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
1926 | px(1)=pxrota | |
1927 | py(1)=pyrota | |
1928 | pz(1)=pzrota | |
1929 | c | |
1930 | pxrota=px(2) | |
1931 | pyrota=py(2) | |
1932 | pzrota=pz(2) | |
1933 | c call rotate(pcx,pcy,pcz,px(2),py(2),pz(2)) | |
1934 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
1935 | px(2)=pxrota | |
1936 | py(2)=pyrota | |
1937 | pz(2)=pzrota | |
1938 | c | |
1939 | pxrota=px(3) | |
1940 | pyrota=py(3) | |
1941 | pzrota=pz(3) | |
1942 | c call rotate(pcx,pcy,pcz,px(3),py(3),pz(3)) | |
1943 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
1944 | px(3)=pxrota | |
1945 | py(3)=pyrota | |
1946 | pz(3)=pzrota | |
1947 | c | |
1948 | pxrota=px(4) | |
1949 | pyrota=py(4) | |
1950 | pzrota=pz(4) | |
1951 | c call rotate(pcx,pcy,pcz,px(4),py(4),pz(4)) | |
1952 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
1953 | px(4)=pxrota | |
1954 | py(4)=pyrota | |
1955 | pz(4)=pzrota | |
1956 | ||
1957 | nnn=nnn+2 | |
1958 | c K+ | |
1959 | lpion(nnn,irun)=23 | |
1960 | if(nchrg.eq.-1.or.nchrg.eq.-2) then | |
1961 | c To keep charge conservation. D-n->nnK0K-, D-D- -> nD-K0K- | |
1962 | ||
1963 | cbz3/7/99 neutralk | |
1964 | c lpion(nnn,irun)=24 ! K0 | |
1965 | cbz3/7/99 neutralk end | |
1966 | ||
1967 | endif | |
1968 | c aka: rest mass of K | |
1969 | epion(nnn,irun)=aka | |
1970 | c K- | |
1971 | lpion(nnn-1,irun)=21 | |
1972 | c aka: rest mass of K | |
1973 | epion(nnn-1,irun)=aka | |
1974 | * Find the momenta of particles in the final state in the nucleus_nucleus | |
1975 | * cms frame. Lorentz transformation into lab frame. | |
1976 | e1cm = sqrt(dm3**2 + px(1)**2 + py(1)**2 + pz(1)**2) | |
1977 | p1beta = px(1)*betax + py(1)*betay + pz(1)*betaz | |
1978 | transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm) | |
1979 | pt1i1 = betax*transf + px(1) | |
1980 | pt2i1 = betay*transf + py(1) | |
1981 | pt3i1 = betaz*transf + pz(1) | |
1982 | eti1 = dm3 | |
1983 | c lb1 = lb(i1) | |
1984 | lb1 = 2 | |
1985 | if(nchrg.ge.-2.and.nchrg.le.1) lb1=2 | |
1986 | ||
1987 | cbz3/7/99 neutralk | |
1988 | if (nchrg .eq. -2 .or. nchrg .eq. -1) then | |
1989 | lb1 = 6 | |
1990 | end if | |
1991 | cbz3/7/99 neutralk end | |
1992 | ||
1993 | cbz3/11/99 neutralk | |
1994 | c if(nchrg.eq.2.or.nchrg.eq.3) lb1=1 | |
1995 | c if(nchrg.eq.4) lb1=9 | |
1996 | if(nchrg.eq.1.or.nchrg.eq.2) lb1=1 | |
1997 | if(nchrg.eq.3.or.nchrg.eq.4) lb1=9 | |
1998 | cbz3/11/99 neutralk end | |
1999 | ||
2000 | * For second nulceon, same | |
2001 | e2cm = sqrt(dm4**2 + px(2)**2 + py(2)**2 + pz(2)**2) | |
2002 | p2beta = px(2)*betax + py(2)*betay + pz(2)*betaz | |
2003 | transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm) | |
2004 | pt1i2 = betax*transf + px(2) | |
2005 | pt2i2 = betay*transf + py(2) | |
2006 | pt3i2 = betaz*transf + pz(2) | |
2007 | eti2 = dm4 | |
2008 | c lb2 = lb(i2) | |
2009 | lb2 = 2 | |
2010 | ||
2011 | cbz3/11/99 neutralk | |
2012 | c if(nchrg.eq.-1.or.nchrg.eq.0) lb2=2 | |
2013 | c if(nchrg.eq. 2.or.nchrg.eq.1) lb2=1 | |
2014 | c if(nchrg.eq. 4.or.nchrg.eq.3) lb2=9 | |
2015 | c if(nchrg.eq.-2) lb2=6 | |
2016 | if(nchrg.ge.-1.or.nchrg.le.1) lb2=2 | |
2017 | if(nchrg.eq. 2.or.nchrg.eq.3) lb2=1 | |
2018 | if(nchrg.eq. 4) lb2=9 | |
2019 | if(nchrg.eq.-2) lb2=6 | |
2020 | cbz3/11/99 neutralk end | |
2021 | ||
2022 | c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then | |
2023 | p(1,i1)=pt1i1 | |
2024 | p(2,i1)=pt2i1 | |
2025 | p(3,i1)=pt3i1 | |
2026 | e(i1)=eti1 | |
2027 | lb(i1)=lb1 | |
2028 | p(1,i2)=pt1i2 | |
2029 | p(2,i2)=pt2i2 | |
2030 | p(3,i2)=pt3i2 | |
2031 | e(i2)=eti2 | |
2032 | lb(i2)=lb2 | |
2033 | ||
2034 | c px1 = p(1,i1) | |
2035 | c py1 = p(2,i1) | |
2036 | c pz1 = p(3,i1) | |
2037 | c em1 = e(i1) | |
2038 | c id(i1) = 2 | |
2039 | c id(i2) = 2 | |
2040 | c id1 = id(i1) | |
2041 | c iblock = 101 ! K(+)K(-) production | |
2042 | * Get anti-kaons' momenta and coordinates in nucleus-nucleus cms. frame. | |
2043 | epcmk = sqrt(epion(nnn-1,irun)**2 + px(3)**2+py(3)**2+pz(3)**2) | |
2044 | betak = px(3)*betax + py(3)*betay + pz(3)*betaz | |
2045 | transf= gamma*(gamma*betak/(gamma+1.) + epcmk) | |
2046 | ppion(1,nnn-1,irun)=betax*transf + px(3) | |
2047 | ppion(2,nnn-1,irun)=betay*transf + py(3) | |
2048 | ppion(3,nnn-1,irun)=betaz*transf + pz(3) | |
2049 | rpion(1,nnn-1,irun)=r(1,i1) | |
2050 | rpion(2,nnn-1,irun)=r(2,i1) | |
2051 | rpion(3,nnn-1,irun)=r(3,i1) | |
2052 | clin-5/2008: | |
2053 | dppion(nnn-1,irun)=dpertp(i1)*dpertp(i2) | |
2054 | * Same thing for kaon ************************************** | |
2055 | epcmak = sqrt(epion(nnn,irun)**2 + px(4)**2 +py(4)**2+pz(4)**2) | |
2056 | betaak = px(4)*betax + py(4)*betay + pz(4)*betaz | |
2057 | transf= gamma*(gamma*betaak/(gamma+1.) + epcmak) | |
2058 | ppion(1,nnn,irun)=betax*transf + px(4) | |
2059 | ppion(2,nnn,irun)=betay*transf + py(4) | |
2060 | ppion(3,nnn,irun)=betaz*transf + pz(4) | |
2061 | rpion(1,nnn,irun)=r(1,i2) | |
2062 | rpion(2,nnn,irun)=r(2,i2) | |
2063 | rpion(3,nnn,irun)=r(3,i2) | |
2064 | clin-5/2008: | |
2065 | dppion(nnn,irun)=dpertp(i1)*dpertp(i2) | |
2066 | return | |
2067 | end | |
2068 | ||
2069 | subroutine lorntz(ilo,b,pi,pj) | |
2070 | c It uses to perform Lorentz (or inverse Lorentz) transformation | |
2071 | dimension pi(4),pj(4),b(3) | |
2072 | SAVE | |
2073 | c dimension db(3) | |
2074 | bb=b(1)*b(1)+b(2)*b(2)+b(3)*b(3) | |
2075 | deno3=sqrt(1.-bb) | |
2076 | if(deno3.eq.0.)deno3=1.e-10 | |
2077 | gam=1./deno3 | |
2078 | ga=gam*gam/(gam+1.) | |
2079 | if(ilo.eq.1) goto 100 | |
2080 | c Lorentz transformation | |
2081 | pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3) | |
2082 | pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3) | |
2083 | c drb=drd(1)*b(1)+drd(2)*b(2)+drd(3)*b(3) | |
2084 | c drdb=db(1)*b(1)+db(2)*b(2)+db(3)*b(3) | |
2085 | do 1001 i=1,3 | |
2086 | pi(i)=pi(i)+b(i)*(ga*pib-gam*pi(4)) | |
2087 | pj(i)=pj(i)+b(i)*(ga*pjb-gam*pj(4)) | |
2088 | c drd(i)=drd(i)+b(i)*ga*drb | |
2089 | c db(i)=db(i)+b(i)*ga*drdb | |
2090 | 1001 continue | |
2091 | pi(4)=gam*(pi(4)-pib) | |
2092 | pj(4)=gam*(pj(4)-pjb) | |
2093 | return | |
2094 | 100 continue | |
2095 | c inverse Lorentz transformation | |
2096 | pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3) | |
2097 | pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3) | |
2098 | do 1002 i=1,3 | |
2099 | pi(i)=pi(i)+b(i)*(ga*pib+gam*pi(4)) | |
2100 | pj(i)=pj(i)+b(i)*(ga*pjb+gam*pj(4)) | |
2101 | 1002 continue | |
2102 | pi(4)=gam*(pi(4)+pib) | |
2103 | pj(4)=gam*(pj(4)+pjb) | |
2104 | return | |
2105 | end | |
2106 | ||
2107 | subroutine fstate(iseed,srt,dm3,dm4,px,py,pz,iflag) | |
2108 | * function: decide final momentum for N,N,K(+),and K(-) | |
2109 | dimension px(4), py(4), pz(4), pe(4) | |
2110 | COMMON/RNDF77/NSEED | |
2111 | cc SAVE /RNDF77/ | |
2112 | SAVE | |
2113 | ||
2114 | iseed=iseed | |
2115 | iflag=-1 | |
2116 | c iflag=-1: fail to find momenta | |
2117 | c = 1: success | |
2118 | pio=3.1415926 | |
2119 | aka=0.498 | |
2120 | c v=0.43 | |
2121 | c w=-0.84 | |
2122 | c b=3.78 | |
2123 | c c=0.47 | |
2124 | c d=3.60 | |
2125 | c fmax=1.056 | |
2126 | c gmax=1.+c | |
2127 | ||
2128 | icount=0 | |
2129 | ekmax=(srt-dm3-dm4)/2. | |
2130 | if(ekmax.le.aka) return | |
2131 | pkmax=sqrt(ekmax**2-aka**2) | |
2132 | ||
2133 | if(dm3.le.0.0.or.dm4.le.0.0) then | |
2134 | write(1,*)'error: minus mass!!!' | |
2135 | return | |
2136 | endif | |
2137 | ||
2138 | c after we have the momenta for both nucleus, we sample the | |
2139 | c transverse momentum for K-. | |
2140 | c dsigma/dpt**2 = exp(-4.145*pt**2) obtained by fitting data on | |
2141 | c page 72, fig 23i. | |
2142 | 50 continue | |
2143 | icount=icount+1 | |
2144 | if(icount.gt.10) return | |
2145 | ptkmi2=-1./4.145*alog(RANART(NSEED)) | |
2146 | ptkm=sqrt(ptkmi2) | |
2147 | 3 v1=RANART(NSEED) | |
2148 | v2=RANART(NSEED) | |
2149 | rsq=v1**2+v2**2 | |
2150 | if(rsq.ge.1.0.or.rsq.le.0.) goto 3 | |
2151 | fac=sqrt(-2.*alog(rsq)/rsq) | |
2152 | guass=v1*fac | |
2153 | if(guass.ge.5.) goto 3 | |
2154 | xstar=guass/5. | |
2155 | pzkm=pkmax*xstar | |
2156 | ekm=sqrt(aka**2+pzkm**2+ptkm**2) | |
2157 | if(RANART(NSEED).gt.aka/ekm) goto 50 | |
2158 | bbb=RANART(NSEED) | |
2159 | px(3)=ptkm*cos(2.*pio*bbb) | |
2160 | py(3)=ptkm*sin(2.*pio*bbb) | |
2161 | if(RANART(NSEED).gt.0.5) pzkm=-1.*pzkm | |
2162 | pz(3)=pzkm | |
2163 | pe(3)=ekm | |
2164 | 150 ptkpl2=-1./3.68*alog(RANART(NSEED)) | |
2165 | ptkp=sqrt(ptkpl2) | |
2166 | 13 v1=RANART(NSEED) | |
2167 | v2=RANART(NSEED) | |
2168 | rsq=v1**2+v2**2 | |
2169 | if(rsq.ge.1.0.or.rsq.le.0.) goto 13 | |
2170 | fac=sqrt(-2.*alog(rsq)/rsq) | |
2171 | guass=v1*fac | |
2172 | if(guass.ge.3.25) goto 13 | |
2173 | xstar=guass/3.25 | |
2174 | pzkp=pkmax*xstar | |
2175 | ekp=sqrt(aka**2+pzkp**2+ptkp**2) | |
2176 | if(RANART(NSEED).gt.aka/ekp) goto 150 | |
2177 | bbb=RANART(NSEED) | |
2178 | px(4)=ptkp*cos(2.*pio*bbb) | |
2179 | py(4)=ptkp*sin(2.*pio*bbb) | |
2180 | if(RANART(NSEED).gt.0.5) pzkp=-1.*pzkp | |
2181 | pz(4)=pzkp | |
2182 | pe(4)=ekp | |
2183 | ||
2184 | resten=srt-pe(3)-pe(4) | |
2185 | restpz=-pz(3)-pz(4) | |
2186 | c resample | |
2187 | if(resten.le.abs(restpz)) goto 50 | |
2188 | restms=sqrt(resten**2-restpz**2) | |
2189 | c resample | |
2190 | if(restms.lt.(dm3+dm4)) goto 50 | |
2191 | ptp2=-1./2.76*alog(RANART(NSEED)) | |
2192 | ptp=sqrt(ptp2) | |
2193 | bbb=RANART(NSEED) | |
2194 | px(2)=ptp*cos(2.*pio*bbb) | |
2195 | py(2)=ptp*sin(2.*pio*bbb) | |
2196 | px(1)=-1.*(px(4)+px(3)+px(2)) | |
2197 | py(1)=-1.*(py(4)+py(3)+py(2)) | |
2198 | c transverse mass for K- | |
2199 | rmt3=sqrt(dm3**2+px(1)**2+py(1)**2) | |
2200 | c transverse mass for K+ | |
2201 | rmt4=sqrt(dm4**2+px(2)**2+py(2)**2) | |
2202 | if(restms.lt.(rmt3+rmt4)) goto 50 | |
2203 | c else: sampling success! | |
2204 | pzcms=sqrt((restms**2-(rmt3+rmt4)**2)* | |
2205 | & (restms**2-(rmt3-rmt4)**2))/2./restms | |
2206 | if(RANART(NSEED).gt.0.5) then | |
2207 | pz(1)=pzcms | |
2208 | pz(2)=-pzcms | |
2209 | else | |
2210 | pz(1)=-pzcms | |
2211 | pz(2)=pzcms | |
2212 | endif | |
2213 | beta=restpz/resten | |
2214 | gama=1./sqrt(1.-beta**2) | |
2215 | pz(1)=pz(1)*gama + beta*gama*sqrt(rmt3**2+pz(1)**2) | |
2216 | pz(2)=pz(2)*gama + beta*gama*sqrt(rmt4**2+pz(2)**2) | |
2217 | pe(1)=sqrt(rmt3**2+pz(1)**2) | |
2218 | pe(2)=sqrt(rmt4**2+pz(2)**2) | |
2219 | ||
2220 | iflag=1 | |
2221 | return | |
2222 | end | |
2223 | ||
2224 | c----------------------------------------------------------------------- | |
2225 | ||
2226 | c.....extracted from G. Song's ART expasion including K- interactions | |
2227 | c.....file `NPIK.FOR' | |
2228 | ||
2229 | **************************************** | |
2230 | c subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt, | |
2231 | c & pcx,pcy,pcz,nchrg,ratiok) | |
2232 | subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt, | |
2233 | & pcx,pcy,pcz,nchrg,ratiok,iblock) | |
2234 | * | |
2235 | * Process: PI + N -> K(-) + ANYTHING | |
2236 | * 1. PI- + P -> P + K0 + K- | |
2237 | * 2. PI+ + N -> P + K+ + K- | |
2238 | * 3. PI0 + P -> P + K+ + K- | |
2239 | * 4. PI0 + N -> P + K0 + K- | |
2240 | * 5. PI0 + N -> N + K+ + K- | |
2241 | * 6. PI- + P -> N + K+ + K- | |
2242 | * 7. PI- + N -> N + K0 + K- | |
2243 | * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494 | |
2244 | **************************************** | |
2245 | PARAMETER (MAXSTR=150001,MAXR=1,PI=3.1415926) | |
2246 | PARAMETER (AKA=0.498) | |
2247 | COMMON /AA/ R(3,MAXSTR) | |
2248 | cc SAVE /AA/ | |
2249 | COMMON /BB/ P(3,MAXSTR) | |
2250 | cc SAVE /BB/ | |
2251 | COMMON /CC/ E(MAXSTR) | |
2252 | cc SAVE /CC/ | |
2253 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
2254 | cc SAVE /EE/ | |
2255 | COMMON /BG/BETAX,BETAY,BETAZ,GAMMA | |
2256 | cc SAVE /BG/ | |
2257 | COMMON /NN/NNN | |
2258 | cc SAVE /NN/ | |
2259 | COMMON /RUN/NUM | |
2260 | cc SAVE /RUN/ | |
2261 | COMMON /PA/RPION(3,MAXSTR,MAXR) | |
2262 | cc SAVE /PA/ | |
2263 | COMMON /PB/PPION(3,MAXSTR,MAXR) | |
2264 | cc SAVE /PB/ | |
2265 | COMMON /PC/EPION(MAXSTR,MAXR) | |
2266 | cc SAVE /PC/ | |
2267 | COMMON /PD/LPION(MAXSTR,MAXR) | |
2268 | cc SAVE /PD/ | |
2269 | dimension bb(3),p1(4),p2(4),p3(4),px(4),py(4),pz(4) | |
2270 | COMMON/RNDF77/NSEED | |
2271 | cc SAVE /RNDF77/ | |
2272 | COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR), | |
2273 | 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), | |
2274 | 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) | |
2275 | SAVE | |
2276 | iseed=iseed | |
2277 | dt=dt | |
2278 | nchrg=nchrg | |
2279 | nt=nt | |
2280 | ratiok=ratiok | |
2281 | px(1)=px(1) | |
2282 | py(1)=py(1) | |
2283 | pz(1)=pz(1) | |
2284 | px1cm=pcx | |
2285 | py1cm=pcy | |
2286 | pz1cm=pcz | |
2287 | ictrl = 1 | |
2288 | lb1=lb(i1) | |
2289 | lb2=lb(i2) | |
2290 | k1=i1 | |
2291 | k2=i2 | |
2292 | c k1 must be bayron. k2 be meson. If not, exchange. | |
2293 | if(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13)) then | |
2294 | k1=i2 | |
2295 | k2=i1 | |
2296 | endif | |
2297 | cbz3/8/99 neutralk | |
2298 | cbz10/12/99 | |
2299 | c LB(I1) = 1 + 2 * RANART(NSEED) | |
2300 | c LB(I2) = 23 | |
2301 | LB(k1) = 1 + int(2*RANART(NSEED)) | |
2302 | LB(k2) = 23 | |
2303 | c pkmax=sqrt((srt**2-(aka+0.938+aka)**2)*(srt**2-(aka+0.938-aka)**2)) | |
2304 | c & /2./srt | |
2305 | pkmax=sqrt((srt**2-(aka+0.938+aka)**2) | |
2306 | & *(srt**2-(aka+0.938-aka)**2))/2./srt | |
2307 | pk = RANART(NSEED)*pkmax | |
2308 | c----------------------------------------------------- | |
2309 | css=1.-2.*RANART(NSEED) | |
2310 | sss=sqrt(1.-css**2) | |
2311 | fai=2*3.1415926*RANART(NSEED) | |
2312 | p3(1)=pk*sss*cos(fai) | |
2313 | p3(2)=pk*sss*sin(fai) | |
2314 | p3(3)=pk*css | |
2315 | eip = srt - sqrt(aka**2 + pk**2) | |
2316 | rmnp=sqrt(eip**2-pk**2) | |
2317 | do 1001 i= 1, 3 | |
2318 | bb(i) = -1.*p3(i)/eip | |
2319 | 1001 continue | |
2320 | c bb: velocity of the other two particles as a whole. | |
2321 | pznp=sqrt((rmnp**2-(aka+0.938)**2) | |
2322 | c *(rmnp**2-(0.938-aka)**2))/2./rmnp | |
2323 | c----------------------------------------------------- | |
2324 | css=1.-2.*RANART(NSEED) | |
2325 | sss=sqrt(1.-css**2) | |
2326 | fai=2*3.1415926*RANART(NSEED) | |
2327 | p1(1)=pznp*sss*cos(fai) | |
2328 | p1(2)=pznp*sss*sin(fai) | |
2329 | p1(3)=pznp*css | |
2330 | p1(4)=sqrt(0.938**2+pznp**2) | |
2331 | p2(4)=sqrt(aka**2+pznp**2) | |
2332 | do 1002 i=1,3 | |
2333 | p2(i)=-1.*p1(i) | |
2334 | 1002 continue | |
2335 | c p1,p2: the momenta of the two particles in their cms | |
2336 | c p1: momentum of N or P | |
2337 | c p2: momentum of anti_kaon | |
2338 | c p3: momentum of K0 or K+ | |
2339 | ilo=1 | |
2340 | c write(61,*)'--------p1,p2',p1,p2 | |
2341 | c write(61,*)'--------bb',bb | |
2342 | call lorntz(ilo,bb,p1,p2) | |
2343 | c******* Checking ************* | |
2344 | c pxsum = p1(1)+p2(1)+p3(1) | |
2345 | c pysum = p1(2)+p2(2)+p3(2) | |
2346 | c pzsum = p1(3)+p2(3)+p3(3) | |
2347 | c pesum = p1(4)+p2(4)+sqrt(p3(1)**2+p3(2)**2+p3(3)**2+aka**2)-srt | |
2348 | c write(61,*)'---p1,pxsum',p1,pxsum | |
2349 | c write(61,*)'---p2,pysum',p2,pysum | |
2350 | c write(61,*)'---p3,pzsum',p3,pzsum | |
2351 | c write(61,*)'---pesum',pesum | |
2352 | c*********************************** | |
2353 | ||
2354 | * Rotate the momenta of particles in the cms of I1 & I2 | |
2355 | * px(1), py(1), pz(1): momentum of I1 | |
2356 | * px(2), py(2), pz(2): momentum of I2 | |
2357 | * px(3), py(3), pz(3): momentum of anti-kaon | |
2358 | ||
2359 | c 10/28/02 get rid of argument usage mismatch in rotate(): | |
2360 | pxrota=p1(1) | |
2361 | pyrota=p1(2) | |
2362 | pzrota=p1(3) | |
2363 | c call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3)) | |
2364 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2365 | p1(1)=pxrota | |
2366 | p1(2)=pyrota | |
2367 | p1(3)=pzrota | |
2368 | c | |
2369 | pxrota=p2(1) | |
2370 | pyrota=p2(2) | |
2371 | pzrota=p2(3) | |
2372 | c call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3)) | |
2373 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2374 | p2(1)=pxrota | |
2375 | p2(2)=pyrota | |
2376 | p2(3)=pzrota | |
2377 | c | |
2378 | pxrota=p3(1) | |
2379 | pyrota=p3(2) | |
2380 | pzrota=p3(3) | |
2381 | c call rotate(pcx,pcy,pcz,p3(1),p3(2),p3(3)) | |
2382 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2383 | p3(1)=pxrota | |
2384 | p3(2)=pyrota | |
2385 | p3(3)=pzrota | |
2386 | ||
2387 | nnn=nnn+1 | |
2388 | c K(-) | |
2389 | lpion(nnn,irun)=21 | |
2390 | c aka: rest mass of K | |
2391 | epion(nnn,irun)=aka | |
2392 | * Find the momenta of particles in the final state in the nucleus_nucleus | |
2393 | * cms frame. Lorentz transformation into lab frame. | |
2394 | e1cm = sqrt(0.938**2 + p1(1)**2 + p1(2)**2 + p1(3)**2) | |
2395 | p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz | |
2396 | transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm) | |
2397 | pt1i1 = betax*transf + p1(1) | |
2398 | pt2i1 = betay*transf + p1(2) | |
2399 | pt3i1 = betaz*transf + p1(3) | |
2400 | eti1 = 0.938 | |
2401 | lb1 = lb(k1) | |
2402 | ||
2403 | * For second nulceon, same | |
2404 | e2cm = sqrt(aka**2 + p3(1)**2 + p3(2)**2 + p3(3)**2) | |
2405 | p2beta = p3(1)*betax + p3(2)*betay + p3(3)*betaz | |
2406 | transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm) | |
2407 | pt1i2 = betax*transf + p3(1) | |
2408 | pt2i2 = betay*transf + p3(2) | |
2409 | pt3i2 = betaz*transf + p3(3) | |
2410 | eti2 = aka | |
2411 | lb2 = lb(k2) | |
2412 | ||
2413 | c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then | |
2414 | * k1 stand for nucleon, k2 stand for kaon. lpion stand for Kbar. | |
2415 | p(1,k1)=pt1i1 | |
2416 | p(2,k1)=pt2i1 | |
2417 | p(3,k1)=pt3i1 | |
2418 | e(k1)=eti1 | |
2419 | lb(k1)=lb1 | |
2420 | p(1,k2)=pt1i2 | |
2421 | p(2,k2)=pt2i2 | |
2422 | p(3,k2)=pt3i2 | |
2423 | e(k2)=eti2 | |
2424 | lb(k2)=lb2 | |
2425 | ||
2426 | c px1 = p(1,i1) | |
2427 | c py1 = p(2,i1) | |
2428 | c pz1 = p(3,i1) | |
2429 | c em1 = e(i1) | |
2430 | c id(i1) = 2 | |
2431 | c id(i2) = 2 | |
2432 | c id1 = id(i1) | |
2433 | c K(+)K(-) production | |
2434 | iblock = 101 | |
2435 | * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame. | |
2436 | c p2: momentum of anti-kaon. | |
2437 | c epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2 + p2(2)**2 + p2(3)**2) | |
2438 | epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2+p2(2)**2+p2(3)**2) | |
2439 | betak = p2(1)*betax + p2(2)*betay + p2(3)*betaz | |
2440 | transf= gamma*(gamma*betak/(gamma+1.) + epcmk) | |
2441 | ppion(1,nnn,irun)=betax*transf + p2(1) | |
2442 | ppion(2,nnn,irun)=betay*transf + p2(2) | |
2443 | ppion(3,nnn,irun)=betaz*transf + p2(3) | |
2444 | clin-5/2008: | |
2445 | dppion(nnn,irun)=dpertp(i1)*dpertp(i2) | |
2446 | cbz3/2/99 | |
2447 | c write(400,*)'2 ', ppion(1,nnn,irun), ppion(2,nnn,irun), | |
2448 | c & ppion(3,nnn,irun), dt*nt, srt | |
2449 | cbz3/2/99end | |
2450 | c write(420,*)ppion(1,nnn,irun), ppion(2,nnn,irun), | |
2451 | c & ppion(3,nnn,irun), dt*nt, srt | |
2452 | k=i2 | |
2453 | if(lb(i1).eq.1.or.lb(i1).eq.2) k=i1 | |
2454 | rpion(1,nnn,irun)=r(1,k) | |
2455 | rpion(2,nnn,irun)=r(2,k) | |
2456 | rpion(3,nnn,irun)=r(3,k) | |
2457 | return | |
2458 | end | |
2459 | ||
2460 | c----------------------------------------------------------------------- | |
2461 | ||
2462 | c.....extracted from G. Song's ART expasion including K- interactions | |
2463 | c.....file `PIHYPN.FOR' | |
2464 | ||
2465 | ****************************************** | |
2466 | subroutine pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2, | |
2467 | & srt,pcx,pcy,pcz,nchrg,iblock) | |
2468 | * | |
2469 | * Process: PI + sigma(or Lambda) -> Kbar + N | |
2470 | * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494 | |
2471 | ****************************************** | |
2472 | ||
2473 | c NOTE: for PI + Hyperon: the produced kaons have mass 0.498 | |
2474 | PARAMETER (MAXSTR=150001,MAXR=1,PI=3.1415926) | |
2475 | PARAMETER (AKA=0.498) | |
2476 | COMMON /AA/ R(3,MAXSTR) | |
2477 | cc SAVE /AA/ | |
2478 | COMMON /BB/ P(3,MAXSTR) | |
2479 | cc SAVE /BB/ | |
2480 | COMMON /CC/ E(MAXSTR) | |
2481 | cc SAVE /CC/ | |
2482 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
2483 | cc SAVE /EE/ | |
2484 | COMMON /BG/BETAX,BETAY,BETAZ,GAMMA | |
2485 | cc SAVE /BG/ | |
2486 | COMMON /NN/NNN | |
2487 | cc SAVE /NN/ | |
2488 | COMMON /RUN/NUM | |
2489 | cc SAVE /RUN/ | |
2490 | COMMON /PA/RPION(3,MAXSTR,MAXR) | |
2491 | cc SAVE /PA/ | |
2492 | COMMON /PB/PPION(3,MAXSTR,MAXR) | |
2493 | cc SAVE /PB/ | |
2494 | COMMON /PC/EPION(MAXSTR,MAXR) | |
2495 | cc SAVE /PC/ | |
2496 | COMMON /PD/LPION(MAXSTR,MAXR) | |
2497 | cc SAVE /PD/ | |
2498 | dimension p1(4),p2(4) | |
2499 | COMMON/RNDF77/NSEED | |
2500 | cc SAVE /RNDF77/ | |
2501 | SAVE | |
2502 | irun=irun | |
2503 | iseed=iseed | |
2504 | nt=nt | |
2505 | dt=dt | |
2506 | px1cm=pcx | |
2507 | py1cm=pcy | |
2508 | pz1cm=pcz | |
2509 | ictrl = 1 | |
2510 | csp06/07/01 | |
2511 | if(ielstc .eq. 1) then | |
2512 | * L/Si + meson -> L/Si + meson | |
2513 | k1=i1 | |
2514 | k2=i2 | |
2515 | dm3=e(k1) | |
2516 | dm4=e(k2) | |
2517 | iblock = 10 | |
2518 | else | |
2519 | iblock = 12 | |
2520 | csp06/07/01 end | |
2521 | c PI + Sigma(or Lambda) -> Kbar + N | |
2522 | k1=i1 | |
2523 | k2=i2 | |
2524 | c k1 must be bayron! So if I1 is PI, exchange k1 & k2. | |
2525 | if(lb(i1).lt.14.or.lb(i1).gt.17) then | |
2526 | k1=i2 | |
2527 | k2=i1 | |
2528 | endif | |
2529 | cbz3/8/99 neutralk | |
2530 | LB(K1) = 1 + int(2*RANART(NSEED)) | |
2531 | if(nchrg.eq.-2) lb(k1)=6 | |
2532 | c if(nchrg.eq.-1) lb(k1)=2 | |
2533 | c if(nchrg.eq. 0) lb(k1)=1 | |
2534 | c if(nchrg.eq. 1) lb(k1)=9 | |
2535 | IF (NCHRG .EQ. 2) LB(K1) = 9 | |
2536 | cbz3/8/99 neutralk end | |
2537 | ||
2538 | c K- | |
2539 | lb(k2)=21 | |
2540 | dm3=0.938 | |
2541 | if(nchrg.eq.-2.or.nchrg.eq.1) dm3=1.232 | |
2542 | dm4=aka | |
2543 | c dm3,dm4: the mass of final state particles. | |
2544 | endif | |
2545 | ||
2546 | ********Now, antikaon will be created. | |
2547 | c call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1) | |
2548 | c pkmax: the maximum momentum of anti-kaon | |
2549 | pkmax=sqrt((srt**2-(dm3+dm4)**2)*(srt**2-(dm3-dm4)**2)) | |
2550 | & /2./srt | |
2551 | pk=pkmax | |
2552 | c----------------------------------------------------- | |
2553 | css=1.-2.*RANART(NSEED) | |
2554 | sss=sqrt(1.-css**2) | |
2555 | fai=2*3.1415926*RANART(NSEED) | |
2556 | p1(1)=pk*sss*cos(fai) | |
2557 | p1(2)=pk*sss*sin(fai) | |
2558 | p1(3)=pk*css | |
2559 | do 1001 i=1,3 | |
2560 | p2(i)=-1.*p1(i) | |
2561 | 1001 continue | |
2562 | c p1,p2: the momenta of the two particles in their cms | |
2563 | c p1: momentum of kaon | |
2564 | c p2: momentum of Kbar | |
2565 | ||
2566 | * Rotate the momenta of particles in the cms of I1 & I2 | |
2567 | clin-10/28/02 get rid of argument usage mismatch in rotate(): | |
2568 | pxrota=p1(1) | |
2569 | pyrota=p1(2) | |
2570 | pzrota=p1(3) | |
2571 | c call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3)) | |
2572 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2573 | p1(1)=pxrota | |
2574 | p1(2)=pyrota | |
2575 | p1(3)=pzrota | |
2576 | c | |
2577 | pxrota=p2(1) | |
2578 | pyrota=p2(2) | |
2579 | pzrota=p2(3) | |
2580 | c call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3)) | |
2581 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2582 | p2(1)=pxrota | |
2583 | p2(2)=pyrota | |
2584 | p2(3)=pzrota | |
2585 | clin-10/28/02-end | |
2586 | ||
2587 | * Find the momenta of particles in the final state in the nucleus_nucleus | |
2588 | * cms frame. Lorentz transformation into lab frame. | |
2589 | e1cm = sqrt(dm3**2 + p1(1)**2 + p1(2)**2 + p1(3)**2) | |
2590 | p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz | |
2591 | transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm) | |
2592 | pt1i1 = betax*transf + p1(1) | |
2593 | pt2i1 = betay*transf + p1(2) | |
2594 | pt3i1 = betaz*transf + p1(3) | |
2595 | eti1 = dm3 | |
2596 | lb1 = lb(k1) | |
2597 | ||
2598 | * For second kaon, same | |
2599 | e2cm = sqrt(dm4**2 + p2(1)**2 + p2(2)**2 + p2(3)**2) | |
2600 | p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz | |
2601 | transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm) | |
2602 | pt1i2 = betax*transf + p2(1) | |
2603 | pt2i2 = betay*transf + p2(2) | |
2604 | pt3i2 = betaz*transf + p2(3) | |
2605 | cbz3/2/99 | |
2606 | c write(400,*)'3 ', pt1i2, pt2i2, pt3i2, dt*nt, srt | |
2607 | cbz3/2/99end | |
2608 | c write(430,*)pt1i2, pt2i2, pt3i2, dt*nt, srt | |
2609 | eti2 = dm4 | |
2610 | lb2 = lb(k2) | |
2611 | ||
2612 | c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then | |
2613 | c k1=i1 | |
2614 | c k2=i2 | |
2615 | * k1 stand for nucleon, k2 stand for kaon. | |
2616 | p(1,k1)=pt1i1 | |
2617 | p(2,k1)=pt2i1 | |
2618 | p(3,k1)=pt3i1 | |
2619 | e(k1)=eti1 | |
2620 | lb(k1)=lb1 | |
2621 | p(1,k2)=pt1i2 | |
2622 | p(2,k2)=pt2i2 | |
2623 | p(3,k2)=pt3i2 | |
2624 | e(k2)=eti2 | |
2625 | lb(k2)=lb2 | |
2626 | ||
2627 | cc iblock = 101 ! K(+)K(-) production | |
2628 | * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame. | |
2629 | return | |
2630 | end | |
2631 | ||
2632 | c----------------------------------------------------------------------- | |
2633 | ||
2634 | c.....extracted from G. Song's ART expasion including K- interactions | |
2635 | c.....file `KAONN.FOR' | |
2636 | ||
2637 | **************************************** | |
2638 | subroutine kaonN(brel,brsgm,irun,iseed,dt,nt, | |
2639 | & ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg) | |
2640 | * | |
2641 | * Process: PI + sigma(or Lambda) <- Kbar + N | |
2642 | * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494 | |
2643 | **************************************** | |
2644 | PARAMETER (MAXSTR=150001,MAXR=1,PI=3.1415926) | |
2645 | PARAMETER (AKA=0.498,ALA=1.1157,ASA=1.1974) | |
2646 | COMMON /AA/ R(3,MAXSTR) | |
2647 | cc SAVE /AA/ | |
2648 | COMMON /BB/ P(3,MAXSTR) | |
2649 | cc SAVE /BB/ | |
2650 | COMMON /CC/ E(MAXSTR) | |
2651 | cc SAVE /CC/ | |
2652 | COMMON /EE/ ID(MAXSTR),LB(MAXSTR) | |
2653 | cc SAVE /EE/ | |
2654 | COMMON /BG/BETAX,BETAY,BETAZ,GAMMA | |
2655 | cc SAVE /BG/ | |
2656 | COMMON /NN/NNN | |
2657 | cc SAVE /NN/ | |
2658 | COMMON /RUN/NUM | |
2659 | cc SAVE /RUN/ | |
2660 | COMMON /PA/RPION(3,MAXSTR,MAXR) | |
2661 | cc SAVE /PA/ | |
2662 | COMMON /PB/PPION(3,MAXSTR,MAXR) | |
2663 | cc SAVE /PB/ | |
2664 | COMMON /PC/EPION(MAXSTR,MAXR) | |
2665 | cc SAVE /PC/ | |
2666 | COMMON /PD/LPION(MAXSTR,MAXR) | |
2667 | cc SAVE /PD/ | |
2668 | dimension p1(4),p2(4) | |
2669 | COMMON/RNDF77/NSEED | |
2670 | cc SAVE /RNDF77/ | |
2671 | SAVE | |
2672 | dt=dt | |
2673 | irun=irun | |
2674 | iseed=iseed | |
2675 | nchrg=nchrg | |
2676 | nt=nt | |
2677 | px1cm=pcx | |
2678 | py1cm=pcy | |
2679 | pz1cm=pcz | |
2680 | ictrl = 1 | |
2681 | c ratio: used for isospin decision. | |
2682 | k1=i1 | |
2683 | k2=i2 | |
2684 | c k1 must be bayron! So if I1 is Kaon, exchange k1 & k2. | |
2685 | if(e(i1).lt.0.5.and.e(i1).gt.0.01) then | |
2686 | k1=i2 | |
2687 | k2=i1 | |
2688 | endif | |
2689 | *** note: for print out only ******************************* | |
2690 | c record kaon's mass | |
2691 | eee=e(k2) | |
2692 | *** end ************** | |
2693 | rrr=RANART(NSEED) | |
2694 | if(rrr.lt.brel) then | |
2695 | c Kbar + N -> Kbar + N | |
2696 | lb1=lb(k1) | |
2697 | lb2=lb(k2) | |
2698 | em1=e(k1) | |
2699 | em2=e(k2) | |
2700 | iblock = 10 | |
2701 | else | |
2702 | iblock = 12 | |
2703 | if(rrr.lt.(brel+brsgm)) then | |
2704 | c nchrg: Net charges of the two incoming particles. | |
2705 | c Kbar + N -> Sigma + PI | |
2706 | em1=asa | |
2707 | em2=0.138 | |
2708 | ||
2709 | cbz3/8/99 neutralk | |
2710 | LB1 = 15 + int(3*RANART(NSEED)) | |
2711 | LB2 = 3 + int(3*RANART(NSEED)) | |
2712 | else | |
2713 | c Kbar + N -> Lambda + PI | |
2714 | em1=ala | |
2715 | em2=0.138 | |
2716 | c LAmbda | |
2717 | lb1=14 | |
2718 | cbz3/8/99 neutralk | |
2719 | LB2 = 3 + int(3*RANART(NSEED)) | |
2720 | c if(nchrg.eq.1) lb2=5 ! K- + D++ -> Lambda + PI+ | |
2721 | c if(nchrg.eq.0) lb2=4 ! K- + p(D+,N*+) -> Lambda + PI0 | |
2722 | c if(nchrg.eq.-1) lb2=3 ! K- + n(D,N*) -> Lambda + PI- | |
2723 | cbz3/8/99 neutralk | |
2724 | ||
2725 | endif | |
2726 | endif | |
2727 | lb(k1)=lb1 | |
2728 | lb(k2)=lb2 | |
2729 | ||
2730 | ********Now, antikaon will be created. | |
2731 | c call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1) | |
2732 | c pkmax: the maximum momentum of anti-kaon | |
2733 | c write(63,*)'srt,em1,em2',srt,em1,em2 | |
2734 | c write(63,*)'-srt,em1,em2',srt,em1,em2 | |
2735 | pkmax=sqrt((srt**2-(em1+em2)**2)*(srt**2-(em1-em2)**2)) | |
2736 | & /2./srt | |
2737 | pk=pkmax | |
2738 | c----------------------------------------------------- | |
2739 | css=1.-2.*RANART(NSEED) | |
2740 | sss=sqrt(1.-css**2) | |
2741 | fai=2*3.1415926*RANART(NSEED) | |
2742 | p1(1)=pk*sss*cos(fai) | |
2743 | p1(2)=pk*sss*sin(fai) | |
2744 | p1(3)=pk*css | |
2745 | do 1001 i=1,3 | |
2746 | p2(i)=-1.*p1(i) | |
2747 | 1001 continue | |
2748 | c p1,p2: the momenta of the two particles in their cms | |
2749 | c p1: momentum of kaon | |
2750 | c p2: momentum of Kbar | |
2751 | ||
2752 | * Rotate the momenta of particles in the cms of I1 & I2 | |
2753 | ||
2754 | clin-10/28/02 get rid of argument usage mismatch in rotate(): | |
2755 | pxrota=p1(1) | |
2756 | pyrota=p1(2) | |
2757 | pzrota=p1(3) | |
2758 | c call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3)) | |
2759 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2760 | p1(1)=pxrota | |
2761 | p1(2)=pyrota | |
2762 | p1(3)=pzrota | |
2763 | c | |
2764 | pxrota=p2(1) | |
2765 | pyrota=p2(2) | |
2766 | pzrota=p2(3) | |
2767 | c call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3)) | |
2768 | call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota) | |
2769 | p2(1)=pxrota | |
2770 | p2(2)=pyrota | |
2771 | p2(3)=pzrota | |
2772 | clin-10/28/02-end | |
2773 | ||
2774 | * Find the momenta of particles in the final state in the nucleus_nucleus | |
2775 | * cms frame. Lorentz transformation into lab frame. | |
2776 | e1cm = sqrt(em1**2 + p1(1)**2 + p1(2)**2 + p1(3)**2) | |
2777 | p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz | |
2778 | transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm) | |
2779 | pt1i1 = betax*transf + p1(1) | |
2780 | pt2i1 = betay*transf + p1(2) | |
2781 | pt3i1 = betaz*transf + p1(3) | |
2782 | eti1 = em1 | |
2783 | ||
2784 | * For second kaon, same | |
2785 | e2cm = sqrt(em2**2 + p2(1)**2 + p2(2)**2 + p2(3)**2) | |
2786 | p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz | |
2787 | transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm) | |
2788 | pt1i2 = betax*transf + p2(1) | |
2789 | pt2i2 = betay*transf + p2(2) | |
2790 | pt3i2 = betaz*transf + p2(3) | |
2791 | eti2 = em2 | |
2792 | ||
2793 | c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then | |
2794 | c k1=i1 | |
2795 | c k2=i2 | |
2796 | * k1 stand for bayron, k2 stand for meson. | |
2797 | p(1,k1)=pt1i1 | |
2798 | p(2,k1)=pt2i1 | |
2799 | p(3,k1)=pt3i1 | |
2800 | e(k1)=eti1 | |
2801 | p(1,k2)=pt1i2 | |
2802 | p(2,k2)=pt2i2 | |
2803 | p(3,k2)=pt3i2 | |
2804 | e(k2)=eti2 | |
2805 | ||
2806 | cc iblock = 101 ! K(+)K(-) production | |
2807 | * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame. | |
2808 | return | |
2809 | end | |
2810 | ||
2811 | c======================================================================= | |
2812 | ||
2813 | clin Below is the previous artana.f: | |
2814 | c======================================================================= | |
2815 | ||
2816 | c.....analysis subroutine before the hadronic space-time evolution | |
2817 | ||
2818 | SUBROUTINE ARTAN1 | |
2819 | PARAMETER (MAXSTR=150001, MAXR=1) | |
2820 | c.....y cut for mt spectrum | |
2821 | cbz3/17/99 | |
2822 | c PARAMETER (YMT1 = -0.4, YMT2 = 0.4) | |
2823 | PARAMETER (YMT1 = -1.0, YMT2 = 1.0) | |
2824 | cbz3/17/99 end | |
2825 | c.....bin width for mt spectrum and y spectrum | |
2826 | clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat: | |
2827 | c PARAMETER (BMT = 0.05, BY = 0.2) | |
2828 | PARAMETER (BMT = 0.05, BY = 0.4) | |
2829 | COMMON /RUN/ NUM | |
2830 | cc SAVE /RUN/ | |
2831 | COMMON /ARERC1/MULTI1(MAXR) | |
2832 | cc SAVE /ARERC1/ | |
2833 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
2834 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
2835 | & FT1(MAXSTR, MAXR), | |
2836 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
2837 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
2838 | cbz3/17/99 | |
2839 | c & dm1k0s(50), DMT1LA(50), DMT1LB(50) | |
2840 | cc SAVE /ARPRC1/ | |
2841 | COMMON /ARANA1/ | |
2842 | & dy1ntb(50), dy1ntp(50), DY1HM(50), | |
2843 | & DY1KP(50), DY1KM(50), DY1K0S(50), | |
2844 | & DY1LA(50), DY1LB(50), DY1PHI(50), | |
2845 | & dm1pip(50), dm1pim(50), DMT1PR(50), | |
2846 | & DMT1PB(50), DMT1KP(50), dm1km(50), | |
2847 | & dm1k0s(50), DMT1LA(50), DMT1LB(50), | |
2848 | & dy1msn(50), DY1PIP(50), DY1PIM(50), | |
2849 | & DY1PI0(50), DY1PR(50), DY1PB(50) | |
2850 | & ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50) | |
2851 | cc SAVE /ARANA1/ | |
2852 | SAVE | |
2853 | ||
2854 | cbz3/17/99 end | |
2855 | DO 1002 J = 1, NUM | |
2856 | DO 1001 I = 1, MULTI1(J) | |
2857 | ITYP = ITYP1(I, J) | |
2858 | PX = PX1(I, J) | |
2859 | PY = PY1(I, J) | |
2860 | PZ = PZ1(I, J) | |
2861 | EE = EE1(I, J) | |
2862 | XM = XM1(I, J) | |
2863 | c 2/24/03 leptons and photons: | |
2864 | if(xm.lt.0.01) goto 200 | |
2865 | ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2) | |
2866 | eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5)) | |
2867 | ||
2868 | XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) | |
2869 | IF (ABS(PZ) .GE. EE) THEN | |
2870 | PRINT *, 'IN ARTAN1' | |
2871 | PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR' | |
2872 | cbzdbg2/16/99 | |
2873 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
2874 | cbzdbg2/16/99 | |
2875 | cbzdbg2/15/99 | |
2876 | PRINT *, ' PZ = ', PZ, ' EE = ', EE | |
2877 | cbzdbg2/16/99 | |
2878 | PRINT *, ' XM = ', XM | |
2879 | cbzdbg2/16/99end | |
2880 | GOTO 200 | |
2881 | c STOP | |
2882 | cbzdbg2/15/99end | |
2883 | END IF | |
2884 | DXMT = XMT - XM | |
2885 | Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) | |
2886 | c.....rapidity cut for the rapidity distribution | |
2887 | IF (ABS(Y) .GE. 10.0) GOTO 100 | |
2888 | clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat: | |
2889 | c IY = 1 + int(ABS(Y) / BY) | |
2890 | c Ieta = 1 + int(ABS(eta) / BY) | |
2891 | IF (ABS(eta) .GE. 10.0) GOTO 100 | |
2892 | IY = 1 + int((Y+10.) / BY) | |
2893 | Ieta = 1 + int((eta+10.) / BY) | |
2894 | ||
2895 | IF (ITYP .LT. -1000) THEN | |
2896 | dy1ntb(IY) = dy1ntb(IY) - 1.0 | |
2897 | END IF | |
2898 | IF (ITYP .GT. 1000) THEN | |
2899 | dy1ntb(IY) = dy1ntb(IY) + 1.0 | |
2900 | END IF | |
2901 | IF (ITYP .EQ. -2212) THEN | |
2902 | dy1ntp(IY) = dy1ntp(IY) - 1.0 | |
2903 | END IF | |
2904 | IF (ITYP .EQ. 2212) THEN | |
2905 | dy1ntp(IY) = dy1ntp(IY) + 1.0 | |
2906 | END IF | |
2907 | c IF (ITYP .EQ. -211 .OR. ITYP .EQ. -321 .OR. | |
2908 | c & ITYP .EQ. -2212) THEN | |
2909 | IF (ITYP .EQ. -2112) THEN | |
2910 | DY1HM(IY) = DY1HM(IY) + 1.0 | |
2911 | END IF | |
2912 | c | |
2913 | IF (LUCHGE(ITYP).ne.0) THEN | |
2914 | DY1CH(IY) = DY1CH(IY) + 1.0 | |
2915 | DE1CH(Ieta) = DE1CH(Ieta) + 1.0 | |
2916 | IF (LUCHGE(ITYP).lt.0) THEN | |
2917 | DY1NEG(IY) = DY1NEG(IY) + 1.0 | |
2918 | DE1NEG(Ieta) = DE1NEG(Ieta) + 1.0 | |
2919 | endif | |
2920 | END IF | |
2921 | ||
2922 | cbz3/17/99 | |
2923 | IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR. | |
2924 | & (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN | |
2925 | dy1msn(IY) = dy1msn(IY) + 1.0 | |
2926 | END IF | |
2927 | IF (ITYP .EQ. 211) THEN | |
2928 | DY1PIP(IY) = DY1PIP(IY) + 1.0 | |
2929 | END IF | |
2930 | IF (ITYP .EQ. -211) THEN | |
2931 | DY1PIM(IY) = DY1PIM(IY) + 1.0 | |
2932 | END IF | |
2933 | IF (ITYP .EQ. 111) THEN | |
2934 | DY1PI0(IY) = DY1PI0(IY) + 1.0 | |
2935 | END IF | |
2936 | IF (ITYP .EQ. 2212) THEN | |
2937 | DY1PR(IY) = DY1PR(IY) + 1.0 | |
2938 | END IF | |
2939 | IF (ITYP .EQ. -2212) THEN | |
2940 | DY1PB(IY) = DY1PB(IY) + 1.0 | |
2941 | END IF | |
2942 | cbz3/17/99 end | |
2943 | IF (ITYP .EQ. 321) THEN | |
2944 | DY1KP(IY) = DY1KP(IY) + 1.0 | |
2945 | END IF | |
2946 | IF (ITYP .EQ. -321) THEN | |
2947 | DY1KM(IY) = DY1KM(IY) + 1.0 | |
2948 | END IF | |
2949 | clin-4/24/03 evaluate K0L instead of K0S, since sometimes we may decay K0S: | |
2950 | c IF (ITYP .EQ. 310) THEN | |
2951 | IF (ITYP .EQ. 130) THEN | |
2952 | DY1K0S(IY) = DY1K0S(IY) + 1.0 | |
2953 | END IF | |
2954 | IF (ITYP .EQ. 3122) THEN | |
2955 | DY1LA(IY) = DY1LA(IY) + 1.0 | |
2956 | END IF | |
2957 | IF (ITYP .EQ. -3122) THEN | |
2958 | DY1LB(IY) = DY1LB(IY) + 1.0 | |
2959 | END IF | |
2960 | IF (ITYP .EQ. 333) THEN | |
2961 | DY1PHI(IY) = DY1PHI(IY) + 1.0 | |
2962 | END IF | |
2963 | ||
2964 | c.....insert rapidity cut for mt spectrum here | |
2965 | 100 IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200 | |
2966 | IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200 | |
2967 | IMT = 1 + int(DXMT / BMT) | |
2968 | IF (ITYP .EQ. 211) THEN | |
2969 | dm1pip(IMT) = dm1pip(IMT) + 1.0 / XMT | |
2970 | END IF | |
2971 | IF (ITYP .EQ. -211) THEN | |
2972 | dm1pim(IMT) = dm1pim(IMT) + | |
2973 | & 1.0 / XMT | |
2974 | END IF | |
2975 | IF (ITYP .EQ. 2212) THEN | |
2976 | DMT1PR(IMT) = DMT1PR(IMT) + 1.0 / XMT | |
2977 | END IF | |
2978 | IF (ITYP .EQ. -2212) THEN | |
2979 | DMT1PB(IMT) = DMT1PB(IMT) + 1.0 / XMT | |
2980 | END IF | |
2981 | IF (ITYP .EQ. 321) THEN | |
2982 | DMT1KP(IMT) = DMT1KP(IMT) + 1.0 / XMT | |
2983 | END IF | |
2984 | IF (ITYP .EQ. -321) THEN | |
2985 | dm1km(IMT) = dm1km(IMT) + 1.0 / XMT | |
2986 | END IF | |
2987 | clin-4/24/03: | |
2988 | c IF (ITYP .EQ. 310) THEN | |
2989 | IF (ITYP .EQ. 130) THEN | |
2990 | dm1k0s(IMT) = dm1k0s(IMT) + 1.0 / XMT | |
2991 | END IF | |
2992 | IF (ITYP .EQ. 3122) THEN | |
2993 | DMT1LA(IMT) = DMT1LA(IMT) + 1.0 / XMT | |
2994 | END IF | |
2995 | IF (ITYP .EQ. -3122) THEN | |
2996 | DMT1LB(IMT) = DMT1LB(IMT) + 1.0 / XMT | |
2997 | END IF | |
2998 | ||
2999 | 200 CONTINUE | |
3000 | 1001 CONTINUE | |
3001 | 1002 CONTINUE | |
3002 | ||
3003 | RETURN | |
3004 | END | |
3005 | ||
3006 | c----------------------------------------------------------------------- | |
3007 | ||
3008 | c.....analysis subroutine after the hadronic space-time evolution | |
3009 | ||
3010 | SUBROUTINE ARTAN2 | |
3011 | ||
3012 | PARAMETER (MAXSTR=150001, MAXR=1) | |
3013 | c.....y cut for mt spectrum | |
3014 | cbz3/17/99 | |
3015 | c PARAMETER (YMT1 = -0.4, YMT2 = 0.4) | |
3016 | PARAMETER (YMT1 = -1.0, YMT2 = 1.0) | |
3017 | cbz3/17/99 end | |
3018 | c.....bin width for mt spectrum and y spectrum | |
3019 | c PARAMETER (BMT = 0.05, BY = 0.2) | |
3020 | PARAMETER (BMT = 0.05, BY = 0.4) | |
3021 | COMMON /RUN/ NUM | |
3022 | cc SAVE /RUN/ | |
3023 | COMMON /ARERC1/MULTI1(MAXR) | |
3024 | cc SAVE /ARERC1/ | |
3025 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
3026 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
3027 | & FT1(MAXSTR, MAXR), | |
3028 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
3029 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
3030 | cbz3/17/99 | |
3031 | c & dm2k0s(50), DMT2LA(50), DMT2LB(50) | |
3032 | cc SAVE /ARPRC1/ | |
3033 | COMMON /ARANA2/ | |
3034 | & dy2ntb(50), dy2ntp(50), DY2HM(50), | |
3035 | & DY2KP(50), DY2KM(50), DY2K0S(50), | |
3036 | & DY2LA(50), DY2LB(50), DY2PHI(50), | |
3037 | & dm2pip(50), dm2pim(50), DMT2PR(50), | |
3038 | & DMT2PB(50), DMT2KP(50), dm2km(50), | |
3039 | & dm2k0s(50), DMT2LA(50), DMT2LB(50), | |
3040 | & dy2msn(50), DY2PIP(50), DY2PIM(50), | |
3041 | & DY2PI0(50), DY2PR(50), DY2PB(50) | |
3042 | & ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50) | |
3043 | cbz3/17/99 end | |
3044 | cc SAVE /ARANA2/ | |
3045 | SAVE | |
3046 | ||
3047 | DO 1002 J = 1, NUM | |
3048 | DO 1001 I = 1, MULTI1(J) | |
3049 | ITYP = ITYP1(I, J) | |
3050 | PX = PX1(I, J) | |
3051 | PY = PY1(I, J) | |
3052 | PZ = PZ1(I, J) | |
3053 | EE = EE1(I, J) | |
3054 | XM = XM1(I, J) | |
3055 | XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) | |
3056 | c 2/24/03 leptons and photons: | |
3057 | if(xm.lt.0.01) goto 200 | |
3058 | ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2) | |
3059 | eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5)) | |
3060 | ||
3061 | IF (ABS(PZ) .GE. EE) THEN | |
3062 | PRINT *, 'IN ARTAN2' | |
3063 | PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR' | |
3064 | cbzdbg2/16/99 | |
3065 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
3066 | cbzdbg2/16/99 | |
3067 | cbzdbg2/15/99 | |
3068 | PRINT *, ' PZ = ', PZ, ' EE = ', EE | |
3069 | cbzdbg2/16/99 | |
3070 | PRINT *, ' XM = ', XM | |
3071 | cbzdbg2/16/99end | |
3072 | GOTO 200 | |
3073 | c STOP | |
3074 | cbzdbg2/15/99end | |
3075 | END IF | |
3076 | DXMT = XMT - XM | |
3077 | Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) | |
3078 | c.....rapidity cut for the rapidity distribution | |
3079 | IF (ABS(Y) .GE. 10.0) GOTO 100 | |
3080 | c IY = 1 + int(ABS(Y) / BY) | |
3081 | c Ieta = 1 + int(ABS(eta) / BY) | |
3082 | IF (ABS(eta) .GE. 10.0) GOTO 100 | |
3083 | IY = 1 + int((Y+10.) / BY) | |
3084 | Ieta = 1 + int((eta+10.) / BY) | |
3085 | ||
3086 | IF (ITYP .LT. -1000) THEN | |
3087 | dy2ntb(IY) = dy2ntb(IY) - 1.0 | |
3088 | END IF | |
3089 | IF (ITYP .GT. 1000) THEN | |
3090 | dy2ntb(IY) = dy2ntb(IY) + 1.0 | |
3091 | END IF | |
3092 | IF (ITYP .EQ. -2212) THEN | |
3093 | dy2ntp(IY) = dy2ntp(IY) - 1.0 | |
3094 | END IF | |
3095 | IF (ITYP .EQ. 2212) THEN | |
3096 | dy2ntp(IY) = dy2ntp(IY) + 1.0 | |
3097 | END IF | |
3098 | IF (ITYP .EQ. -2112) THEN | |
3099 | DY2HM(IY) = DY2HM(IY) + 1.0 | |
3100 | END IF | |
3101 | ||
3102 | IF (LUCHGE(ITYP).ne.0) THEN | |
3103 | DY2CH(IY) = DY2CH(IY) + 1.0 | |
3104 | DE2CH(Ieta) = DE2CH(Ieta) + 1.0 | |
3105 | IF (LUCHGE(ITYP).lt.0) THEN | |
3106 | DY2NEG(IY) = DY2NEG(IY) + 1.0 | |
3107 | DE2NEG(Ieta) = DE2NEG(Ieta) + 1.0 | |
3108 | endif | |
3109 | END IF | |
3110 | ||
3111 | cbz3/17/99 | |
3112 | IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR. | |
3113 | & (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN | |
3114 | dy2msn(IY) = dy2msn(IY) + 1.0 | |
3115 | END IF | |
3116 | IF (ITYP .EQ. 211) THEN | |
3117 | DY2PIP(IY) = DY2PIP(IY) + 1.0 | |
3118 | END IF | |
3119 | IF (ITYP .EQ. -211) THEN | |
3120 | DY2PIM(IY) = DY2PIM(IY) + 1.0 | |
3121 | END IF | |
3122 | IF (ITYP .EQ. 111) THEN | |
3123 | DY2PI0(IY) = DY2PI0(IY) + 1.0 | |
3124 | END IF | |
3125 | IF (ITYP .EQ. 2212) THEN | |
3126 | DY2PR(IY) = DY2PR(IY) + 1.0 | |
3127 | END IF | |
3128 | IF (ITYP .EQ. -2212) THEN | |
3129 | DY2PB(IY) = DY2PB(IY) + 1.0 | |
3130 | END IF | |
3131 | cbz3/17/99 end | |
3132 | IF (ITYP .EQ. 321) THEN | |
3133 | DY2KP(IY) = DY2KP(IY) + 1.0 | |
3134 | END IF | |
3135 | IF (ITYP .EQ. -321) THEN | |
3136 | DY2KM(IY) = DY2KM(IY) + 1.0 | |
3137 | END IF | |
3138 | clin-4/24/03: | |
3139 | c IF (ITYP .EQ. 310) THEN | |
3140 | IF (ITYP .EQ. 130) THEN | |
3141 | DY2K0S(IY) = DY2K0S(IY) + 1.0 | |
3142 | END IF | |
3143 | IF (ITYP .EQ. 3122) THEN | |
3144 | DY2LA(IY) = DY2LA(IY) + 1.0 | |
3145 | END IF | |
3146 | IF (ITYP .EQ. -3122) THEN | |
3147 | DY2LB(IY) = DY2LB(IY) + 1.0 | |
3148 | END IF | |
3149 | IF (ITYP .EQ. 333) THEN | |
3150 | DY2PHI(IY) = DY2PHI(IY) + 1.0 | |
3151 | END IF | |
3152 | ||
3153 | c.....insert rapidity cut for mt spectrum here | |
3154 | 100 IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200 | |
3155 | IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200 | |
3156 | IMT = 1 + int(DXMT / BMT) | |
3157 | IF (ITYP .EQ. 211) THEN | |
3158 | dm2pip(IMT) = dm2pip(IMT) + 1.0 / XMT | |
3159 | END IF | |
3160 | IF (ITYP .EQ. -211) THEN | |
3161 | dm2pim(IMT) = dm2pim(IMT) + | |
3162 | & 1.0 / XMT | |
3163 | END IF | |
3164 | IF (ITYP .EQ. 2212) THEN | |
3165 | DMT2PR(IMT) = DMT2PR(IMT) + 1.0 / XMT | |
3166 | END IF | |
3167 | IF (ITYP .EQ. -2212) THEN | |
3168 | DMT2PB(IMT) = DMT2PB(IMT) + 1.0 / XMT | |
3169 | END IF | |
3170 | IF (ITYP .EQ. 321) THEN | |
3171 | DMT2KP(IMT) = DMT2KP(IMT) + 1.0 / XMT | |
3172 | END IF | |
3173 | IF (ITYP .EQ. -321) THEN | |
3174 | dm2km(IMT) = dm2km(IMT) + 1.0 / XMT | |
3175 | END IF | |
3176 | clin-4/24/03: | |
3177 | c IF (ITYP .EQ. 310) THEN | |
3178 | IF (ITYP .EQ. 130) THEN | |
3179 | dm2k0s(IMT) = dm2k0s(IMT) + 1.0 / XMT | |
3180 | END IF | |
3181 | IF (ITYP .EQ. 3122) THEN | |
3182 | DMT2LA(IMT) = DMT2LA(IMT) + 1.0 / XMT | |
3183 | END IF | |
3184 | IF (ITYP .EQ. -3122) THEN | |
3185 | DMT2LB(IMT) = DMT2LB(IMT) + 1.0 / XMT | |
3186 | END IF | |
3187 | ||
3188 | 200 CONTINUE | |
3189 | 1001 CONTINUE | |
3190 | 1002 CONTINUE | |
3191 | ||
3192 | RETURN | |
3193 | END | |
3194 | ||
3195 | c----------------------------------------------------------------------- | |
3196 | ||
3197 | c.....output analysis results at the end of the simulation | |
3198 | ||
3199 | SUBROUTINE ARTOUT(NEVNT) | |
3200 | ||
3201 | PARAMETER (MAXSTR=150001, MAXR=1) | |
3202 | c.....y cut for mt spectrum | |
3203 | cbz3/17/99 | |
3204 | c PARAMETER (YMT1 = -0.4, YMT2 = 0.4) | |
3205 | PARAMETER (YMT1 = -1.0, YMT2 = 1.0) | |
3206 | cbz3/17/99 end | |
3207 | c.....bin width for mt spectrum and y spectrum | |
3208 | c PARAMETER (BMT = 0.05, BY = 0.2) | |
3209 | PARAMETER (BMT = 0.05, BY = 0.4) | |
3210 | COMMON /RUN/ NUM | |
3211 | cc SAVE /RUN/ | |
3212 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
3213 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
3214 | & FT1(MAXSTR, MAXR), | |
3215 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
3216 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
3217 | cbz3/17/99 | |
3218 | c & dm1k0s(50), DMT1LA(50), DMT1LB(50) | |
3219 | cc SAVE /ARPRC1/ | |
3220 | COMMON /ARANA1/ | |
3221 | & dy1ntb(50), dy1ntp(50), DY1HM(50), | |
3222 | & DY1KP(50), DY1KM(50), DY1K0S(50), | |
3223 | & DY1LA(50), DY1LB(50), DY1PHI(50), | |
3224 | & dm1pip(50), dm1pim(50), DMT1PR(50), | |
3225 | & DMT1PB(50), DMT1KP(50), dm1km(50), | |
3226 | & dm1k0s(50), DMT1LA(50), DMT1LB(50), | |
3227 | & dy1msn(50), DY1PIP(50), DY1PIM(50), | |
3228 | & DY1PI0(50), DY1PR(50), DY1PB(50) | |
3229 | & ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50) | |
3230 | cbz3/17/99 end | |
3231 | cc SAVE /ARANA1/ | |
3232 | cbz3/17/99 | |
3233 | c & dm2k0s(50), DMT2LA(50), DMT2LB(50) | |
3234 | COMMON /ARANA2/ | |
3235 | & dy2ntb(50), dy2ntp(50), DY2HM(50), | |
3236 | & DY2KP(50), DY2KM(50), DY2K0S(50), | |
3237 | & DY2LA(50), DY2LB(50), DY2PHI(50), | |
3238 | & dm2pip(50), dm2pim(50), DMT2PR(50), | |
3239 | & DMT2PB(50), DMT2KP(50), dm2km(50), | |
3240 | & dm2k0s(50), DMT2LA(50), DMT2LB(50), | |
3241 | & dy2msn(50), DY2PIP(50), DY2PIM(50), | |
3242 | & DY2PI0(50), DY2PR(50), DY2PB(50) | |
3243 | & ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50) | |
3244 | cc SAVE /ARANA2/ | |
3245 | SAVE | |
3246 | cbz3/17/99 end | |
3247 | cms OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN') | |
3248 | cms OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN') | |
3249 | cms OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN') | |
3250 | cms OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN') | |
3251 | cms OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN') | |
3252 | cms OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN') | |
3253 | cms OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN') | |
3254 | cms OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN') | |
3255 | cms OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN') | |
3256 | clin-4/24/03 | |
3257 | c OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN') | |
3258 | cms OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN') | |
3259 | cms OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN') | |
3260 | cms OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN') | |
3261 | cms OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN') | |
3262 | cbz3/17/99 | |
3263 | cms OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN') | |
3264 | cms OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN') | |
3265 | cms OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN') | |
3266 | cms OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN') | |
3267 | cms OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN') | |
3268 | cms OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN') | |
3269 | cbz3/17/99 end | |
3270 | ||
3271 | cms OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN') | |
3272 | cms OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN') | |
3273 | cms OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN') | |
3274 | cms OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN') | |
3275 | cms OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN') | |
3276 | cms OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN') | |
3277 | cms OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN') | |
3278 | cms OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN') | |
3279 | cms OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN') | |
3280 | clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat: | |
3281 | c SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0 | |
3282 | SCALE1 = 1. / REAL(NEVNT * NUM) / BY | |
3283 | SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1) | |
3284 | c | |
3285 | DO 1001 I = 1, 50 | |
3286 | ymid=-10.+BY * (I - 0.5) | |
3287 | cms WRITE (30, 333) ymid, SCALE1 * dy1ntb(I) | |
3288 | cms WRITE (31, 333) ymid, SCALE1 * dy1ntp(I) | |
3289 | cms WRITE (32, 333) ymid, SCALE1 * DY1HM(I) | |
3290 | cms WRITE (37, 333) ymid, SCALE1 * DY1KP(I) | |
3291 | cms WRITE (38, 333) ymid, SCALE1 * DY1KM(I) | |
3292 | cms WRITE (39, 333) ymid, SCALE1 * DY1K0S(I) | |
3293 | cms WRITE (40, 333) ymid, SCALE1 * DY1LA(I) | |
3294 | cms WRITE (41, 333) ymid, SCALE1 * DY1LB(I) | |
3295 | cms WRITE (42, 333) ymid, SCALE1 * DY1PHI(I) | |
3296 | cms WRITE (33, 333) ymid, SCALE1 * DY1NEG(I) | |
3297 | cms WRITE (34, 333) ymid, SCALE1 * DY1CH(I) | |
3298 | cms WRITE (35, 333) ymid, SCALE1 * DE1NEG(I) | |
3299 | cms WRITE (36, 333) ymid, SCALE1 * DE1CH(I) | |
3300 | cms WRITE (43, 333) ymid, SCALE1 * dy1msn(I) | |
3301 | cms WRITE (44, 333) ymid, SCALE1 * DY1PIP(I) | |
3302 | cms WRITE (45, 333) ymid, SCALE1 * DY1PIM(I) | |
3303 | cms WRITE (46, 333) ymid, SCALE1 * DY1PI0(I) | |
3304 | cms WRITE (47, 333) ymid, SCALE1 * DY1PR(I) | |
3305 | cms WRITE (48, 333) ymid, SCALE1 * DY1PB(I) | |
3306 | ||
3307 | IF (dm1pip(I) .NE. 0.0) THEN | |
3308 | cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I) | |
3309 | END IF | |
3310 | IF (dm1pim(I) .NE. 0.0) THEN | |
3311 | cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * | |
3312 | cms & dm1pim(I) | |
3313 | END IF | |
3314 | IF (DMT1PR(I) .NE. 0.0) THEN | |
3315 | cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I) | |
3316 | END IF | |
3317 | IF (DMT1PB(I) .NE. 0.0) THEN | |
3318 | cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I) | |
3319 | END IF | |
3320 | IF (DMT1KP(I) .NE. 0.0) THEN | |
3321 | cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I) | |
3322 | END IF | |
3323 | IF (dm1km(I) .NE. 0.0) THEN | |
3324 | cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * | |
3325 | cms & dm1km(I) | |
3326 | END IF | |
3327 | IF (dm1k0s(I) .NE. 0.0) THEN | |
3328 | cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I) | |
3329 | END IF | |
3330 | IF (DMT1LA(I) .NE. 0.0) THEN | |
3331 | cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I) | |
3332 | END IF | |
3333 | IF (DMT1LB(I) .NE. 0.0) THEN | |
3334 | cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I) | |
3335 | END IF | |
3336 | 1001 CONTINUE | |
3337 | c | |
3338 | DO 1002 I = 30, 48 | |
3339 | cms WRITE (I, *) 'after hadron evolution' | |
3340 | 1002 CONTINUE | |
3341 | DO 1003 I = 50, 58 | |
3342 | cms WRITE (I, *) 'after hadron evolution' | |
3343 | 1003 CONTINUE | |
3344 | ||
3345 | DO 1004 I = 1, 50 | |
3346 | ymid=-10.+BY * (I - 0.5) | |
3347 | cms WRITE (30, 333) ymid, SCALE1 * dy2ntb(I) | |
3348 | cms WRITE (31, 333) ymid, SCALE1 * dy2ntp(I) | |
3349 | cms WRITE (32, 333) ymid, SCALE1 * DY2HM(I) | |
3350 | cms WRITE (37, 333) ymid, SCALE1 * DY2KP(I) | |
3351 | cms WRITE (38, 333) ymid, SCALE1 * DY2KM(I) | |
3352 | cms WRITE (39, 333) ymid, SCALE1 * DY2K0S(I) | |
3353 | cms WRITE (40, 333) ymid, SCALE1 * DY2LA(I) | |
3354 | cms WRITE (41, 333) ymid, SCALE1 * DY2LB(I) | |
3355 | cms WRITE (42, 333) ymid, SCALE1 * DY2PHI(I) | |
3356 | cms WRITE (33, 333) ymid, SCALE1 * DY2NEG(I) | |
3357 | cms WRITE (34, 333) ymid, SCALE1 * DY2CH(I) | |
3358 | cms WRITE (35, 333) ymid, SCALE1 * DE2NEG(I) | |
3359 | cms WRITE (36, 333) ymid, SCALE1 * DE2CH(I) | |
3360 | cms WRITE (43, 333) ymid, SCALE1 * dy2msn(I) | |
3361 | cms WRITE (44, 333) ymid, SCALE1 * DY2PIP(I) | |
3362 | cms WRITE (45, 333) ymid, SCALE1 * DY2PIM(I) | |
3363 | cms WRITE (46, 333) ymid, SCALE1 * DY2PI0(I) | |
3364 | cms WRITE (47, 333) ymid, SCALE1 * DY2PR(I) | |
3365 | cms WRITE (48, 333) ymid, SCALE1 * DY2PB(I) | |
3366 | c | |
3367 | IF (dm2pip(I) .NE. 0.0) THEN | |
3368 | cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I) | |
3369 | END IF | |
3370 | IF (dm2pim(I) .NE. 0.0) THEN | |
3371 | cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * | |
3372 | cms & dm2pim(I) | |
3373 | END IF | |
3374 | IF (DMT2PR(I) .NE. 0.0) THEN | |
3375 | cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I) | |
3376 | END IF | |
3377 | IF (DMT2PB(I) .NE. 0.0) THEN | |
3378 | cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I) | |
3379 | END IF | |
3380 | IF (DMT2KP(I) .NE. 0.0) THEN | |
3381 | cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I) | |
3382 | END IF | |
3383 | IF (dm2km(I) .NE. 0.0) THEN | |
3384 | cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * | |
3385 | cms & dm2km(I) | |
3386 | END IF | |
3387 | IF (dm2k0s(I) .NE. 0.0) THEN | |
3388 | cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I) | |
3389 | END IF | |
3390 | IF (DMT2LA(I) .NE. 0.0) THEN | |
3391 | cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I) | |
3392 | END IF | |
3393 | IF (DMT2LB(I) .NE. 0.0) THEN | |
3394 | cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I) | |
3395 | END IF | |
3396 | 1004 CONTINUE | |
3397 | cms 333 format(2(f12.5,1x)) | |
3398 | ||
3399 | RETURN | |
3400 | END | |
3401 | ||
3402 | c----------------------------------------------------------------------- | |
3403 | ||
3404 | c.....analysis subroutine in HIJING before parton cascade evolution | |
3405 | SUBROUTINE HJANA1 | |
3406 | ||
3407 | PARAMETER (YMAX = 1.0, YMIN = -1.0) | |
3408 | PARAMETER (DMT = 0.05, DY = 0.2) | |
3409 | PARAMETER (DR = 0.2) | |
3410 | PARAMETER (MAXPTN=400001,MAXSTR=150001) | |
3411 | DIMENSION dyp1(50), DMYP1(200), DEYP1(50) | |
3412 | DIMENSION dyg1(50), DMYG1(200), DEYG1(50) | |
3413 | DIMENSION SNYP1(50), SMYP1(200), SEYP1(50) | |
3414 | DIMENSION SNYG1(50), SMYG1(200), SEYG1(50) | |
3415 | DIMENSION dnrpj1(50), dnrtg1(50), dnrin1(50), | |
3416 | & dnrtt1(50) | |
3417 | DIMENSION dyg1c(50), dmyg1c(50), deyg1c(50) | |
3418 | DIMENSION snrpj1(50), snrtg1(50), snrin1(50), | |
3419 | & snrtt1(50) | |
3420 | DIMENSION snyg1c(50), smyg1c(50), seyg1c(50) | |
3421 | DOUBLE PRECISION GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0 | |
3422 | ||
3423 | COMMON /PARA1/ MUL | |
3424 | cc SAVE /PARA1/ | |
3425 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
3426 | cc SAVE /HPARNT/ | |
3427 | COMMON/hjcrdn/YP(3,300),YT(3,300) | |
3428 | cc SAVE /hjcrdn/ | |
3429 | COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500), | |
3430 | & PJPY(300,500),PJPZ(300,500),PJPE(300,500), | |
3431 | & PJPM(300,500),NTJ(300),KFTJ(300,500), | |
3432 | & PJTX(300,500),PJTY(300,500),PJTZ(300,500), | |
3433 | & PJTE(300,500),PJTM(300,500) | |
3434 | cc SAVE /HJJET1/ | |
3435 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
3436 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
3437 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
3438 | cc SAVE /HJJET2/ | |
3439 | COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN), | |
3440 | & PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN), | |
3441 | & XMASS0(MAXPTN), ITYP0(MAXPTN) | |
3442 | cc SAVE /prec1/ | |
3443 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
3444 | cc SAVE /AREVT/ | |
3445 | COMMON /AROUT/ IOUT | |
3446 | cc SAVE /AROUT/ | |
3447 | SAVE | |
3448 | DATA IW/0/ | |
3449 | IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN | |
3450 | DO 1001 I = 1, 200 | |
3451 | DMYP1(I) = SMYP1(I) | |
3452 | DMYG1(I) = SMYG1(I) | |
3453 | 1001 CONTINUE | |
3454 | ||
3455 | DO 1002 I = 1, 50 | |
3456 | dyp1(I) = SNYP1(I) | |
3457 | DEYP1(I) = SEYP1(I) | |
3458 | dyg1(I) = SNYG1(I) | |
3459 | DEYG1(I) = SEYG1(I) | |
3460 | dnrpj1(I) = snrpj1(I) | |
3461 | dnrtg1(I) = snrtg1(I) | |
3462 | dnrin1(I) = snrin1(I) | |
3463 | dnrtt1(I) = snrtt1(I) | |
3464 | dyg1c(I) = snyg1c(I) | |
3465 | dmyg1c(I) = smyg1c(I) | |
3466 | deyg1c(I) = seyg1c(I) | |
3467 | 1002 CONTINUE | |
3468 | nsubp = nsubpS | |
3469 | nsubg = nsubgS | |
3470 | NISG = NISGS | |
3471 | ELSE | |
3472 | DO 1003 I = 1, 200 | |
3473 | SMYP1(I) = DMYP1(I) | |
3474 | SMYG1(I) = DMYG1(I) | |
3475 | 1003 CONTINUE | |
3476 | ||
3477 | DO 1004 I = 1, 50 | |
3478 | SNYP1(I) = dyp1(I) | |
3479 | SEYP1(I) = DEYP1(I) | |
3480 | SNYG1(I) = dyg1(I) | |
3481 | SEYG1(I) = DEYG1(I) | |
3482 | snrpj1(I) = dnrpj1(I) | |
3483 | snrtg1(I) = dnrtg1(I) | |
3484 | snrin1(I) = dnrin1(I) | |
3485 | snrtt1(I) = dnrtt1(I) | |
3486 | snyg1c(I) = dyg1c(I) | |
3487 | smyg1c(I) = dmyg1c(I) | |
3488 | seyg1c(I) = deyg1c(I) | |
3489 | 1004 CONTINUE | |
3490 | nsubpS = nsubp | |
3491 | nsubgS = nsubg | |
3492 | NISGS = NISG | |
3493 | isevt = IAEVT | |
3494 | isrun = IARUN | |
3495 | IW = IW + 1 | |
3496 | END IF | |
3497 | c.....analysis | |
3498 | DO 1006 I = 1, IHNT2(1) | |
3499 | DO 1005 J = 1, NPJ(I) | |
3500 | ITYP = KFPJ(I, J) | |
3501 | PX = PJPX(I, J) | |
3502 | PY = PJPY(I, J) | |
3503 | PZ = PJPZ(I, J) | |
3504 | PE = PJPE(I, J) | |
3505 | PM = PJPM(I, J) | |
3506 | IF (ABS(PZ) .GE. PE) THEN | |
3507 | PRINT *, ' IN HJANA1, PROJ STR ', I, ' PART ', J | |
3508 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
3509 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
3510 | PRINT *, ' XM = ', PM | |
3511 | GOTO 200 | |
3512 | END IF | |
3513 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
3514 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
3515 | DXMT = XMT - PM | |
3516 | IY = 1 + int(ABS(RAP) / DY) | |
3517 | IF (IY .GT. 50) GOTO 100 | |
3518 | dyp1(IY) = dyp1(IY) + 1.0 | |
3519 | DEYP1(IY) = DEYP1(IY) + XMT | |
3520 | IF (ITYP .EQ. 21) THEN | |
3521 | dyg1(IY) = dyg1(IY) + 1.0 | |
3522 | DEYG1(IY) = DEYG1(IY) + XMT | |
3523 | END IF | |
3524 | 100 CONTINUE | |
3525 | IMT = 1 + int(DXMT / DMT) | |
3526 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 200 | |
3527 | IF (IMT .GT. 200) GOTO 200 | |
3528 | DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT | |
3529 | IF (ITYP .EQ. 21) THEN | |
3530 | DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT | |
3531 | END IF | |
3532 | 200 CONTINUE | |
3533 | 1005 CONTINUE | |
3534 | 1006 CONTINUE | |
3535 | ||
3536 | DO 1008 I = 1, IHNT2(3) | |
3537 | DO 1007 J = 1, NTJ(I) | |
3538 | ITYP = KFTJ(I, J) | |
3539 | PX = PJTX(I, J) | |
3540 | PY = PJTY(I, J) | |
3541 | PZ = PJTZ(I, J) | |
3542 | PE = PJTE(I, J) | |
3543 | PM = PJTM(I, J) | |
3544 | IF (ABS(PZ) .GE. PE) THEN | |
3545 | PRINT *, ' IN HJANA1, TARG STR ', I, ' PART ', J | |
3546 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
3547 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
3548 | PRINT *, ' XM = ', PM | |
3549 | GOTO 400 | |
3550 | END IF | |
3551 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
3552 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
3553 | DXMT = XMT - PM | |
3554 | IY = 1 + int(ABS(RAP) / DY) | |
3555 | IF (IY .GT. 50) GOTO 300 | |
3556 | dyp1(IY) = dyp1(IY) + 1.0 | |
3557 | DEYP1(IY) = DEYP1(IY) + XMT | |
3558 | IF (ITYP .EQ. 21) THEN | |
3559 | dyg1(IY) = dyg1(IY) + 1.0 | |
3560 | DEYG1(IY) = DEYG1(IY) + XMT | |
3561 | END IF | |
3562 | 300 CONTINUE | |
3563 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 400 | |
3564 | IMT = 1 + int(DXMT / DMT) | |
3565 | IF (IMT .GT. 200) GOTO 400 | |
3566 | DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT | |
3567 | IF (ITYP .EQ. 21) THEN | |
3568 | DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT | |
3569 | END IF | |
3570 | 400 CONTINUE | |
3571 | 1007 CONTINUE | |
3572 | 1008 CONTINUE | |
3573 | ||
3574 | DO 1010 I = 1, NSG | |
3575 | DO 1009 J = 1, NJSG(I) | |
3576 | ITYP = K2SG(I, J) | |
3577 | PX = PXSG(I, J) | |
3578 | PY = PYSG(I, J) | |
3579 | PZ = PZSG(I, J) | |
3580 | PE = PESG(I, J) | |
3581 | PM = PMSG(I, J) | |
3582 | IF (ABS(PZ) .GE. PE) THEN | |
3583 | PRINT *, ' IN HJANA1, INDP STR ', I, ' PART ', J | |
3584 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
3585 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
3586 | PRINT *, ' XM = ', PM | |
3587 | GOTO 600 | |
3588 | END IF | |
3589 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
3590 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
3591 | DXMT = XMT - PM | |
3592 | IY = 1 + int(ABS(RAP) / DY) | |
3593 | IF (IY .GT. 50) GOTO 500 | |
3594 | dyp1(IY) = dyp1(IY) + 1.0 | |
3595 | DEYP1(IY) = DEYP1(IY) + XMT | |
3596 | IF (ITYP .EQ. 21) THEN | |
3597 | dyg1(IY) = dyg1(IY) + 1.0 | |
3598 | DEYG1(IY) = DEYG1(IY) + XMT | |
3599 | END IF | |
3600 | 500 CONTINUE | |
3601 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 600 | |
3602 | IMT = 1 + int(DXMT / DMT) | |
3603 | IF (IMT .GT. 200) GOTO 600 | |
3604 | DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT | |
3605 | IF (ITYP .EQ. 21) THEN | |
3606 | DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT | |
3607 | END IF | |
3608 | 600 CONTINUE | |
3609 | 1009 CONTINUE | |
3610 | 1010 CONTINUE | |
3611 | ||
3612 | DO 1011 I = 1, IHNT2(1) | |
3613 | YR = SQRT(YP(1, I) ** 2 + YP(2, I) ** 2) | |
3614 | IR = 1 + int(YR / DR) | |
3615 | clin-4/2008 protect against out-of-bound errors: | |
3616 | c IF (IR .GT. 50) GOTO 601 | |
3617 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601 | |
3618 | dnrpj1(IR) = dnrpj1(IR) + 1.0 | |
3619 | dnrtt1(IR) = dnrtt1(IR) + 1.0 | |
3620 | 601 CONTINUE | |
3621 | 1011 CONTINUE | |
3622 | ||
3623 | DO 1012 I = 1, IHNT2(3) | |
3624 | YR = SQRT(YT(1, I) ** 2 + YT(2, I) ** 2) | |
3625 | IR = 1 + int(YR / DR) | |
3626 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 602 | |
3627 | dnrtg1(IR) = dnrtg1(IR) + 1.0 | |
3628 | dnrtt1(IR) = dnrtt1(IR) + 1.0 | |
3629 | 602 CONTINUE | |
3630 | 1012 CONTINUE | |
3631 | ||
3632 | DO 1013 I = 1, NSG | |
3633 | Y1 = 0.5 * (YP(1, IASG(I, 1)) + YT(1, IASG(I, 2))) | |
3634 | Y2 = 0.5 * (YP(2, IASG(I, 1)) + YT(2, IASG(I, 2))) | |
3635 | YR = SQRT(Y1 ** 2 + Y2 ** 2) | |
3636 | IR = 1 + int(YR / DR) | |
3637 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 603 | |
3638 | dnrin1(IR) = dnrin1(IR) + 1.0 | |
3639 | dnrtt1(IR) = dnrtt1(IR) + 1.0 | |
3640 | 603 CONTINUE | |
3641 | 1013 CONTINUE | |
3642 | ||
3643 | DO 1014 I = 1, MUL | |
3644 | ITYP = ITYP0(I) | |
3645 | PX = sngl(PX0(I)) | |
3646 | PY = sngl(PY0(I)) | |
3647 | PZ = sngl(PZ0(I)) | |
3648 | PE = sngl(E0(I)) | |
3649 | PM = sngl(XMASS0(I)) | |
3650 | IF (ABS(PZ) .GE. PE) THEN | |
3651 | PRINT *, ' IN HJANA1, GLUON ', I | |
3652 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
3653 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
3654 | PRINT *, ' XM = ', PM | |
3655 | GOTO 800 | |
3656 | END IF | |
3657 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
3658 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
3659 | DXMT = XMT - PM | |
3660 | IY = 1 + int(ABS(RAP) / DY) | |
3661 | IF (IY .GT. 50) GOTO 700 | |
3662 | dyg1c(IY) = dyg1c(IY) + 1.0 | |
3663 | deyg1c(IY) = deyg1c(IY) + XMT | |
3664 | 700 CONTINUE | |
3665 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 800 | |
3666 | IMT = 1 + int(DXMT / DMT) | |
3667 | IF (IMT .GT. 50) GOTO 800 | |
3668 | dmyg1c(IMT) = dmyg1c(IMT) + 1.0 / XMT | |
3669 | 800 CONTINUE | |
3670 | 1014 CONTINUE | |
3671 | c.....count number of particles | |
3672 | DO 1016 I = 1, IHNT2(1) | |
3673 | DO 1015 J = 1, NPJ(I) | |
3674 | nsubp = nsubp + 1 | |
3675 | IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1 | |
3676 | 1015 CONTINUE | |
3677 | 1016 CONTINUE | |
3678 | ||
3679 | DO 1018 I = 1, IHNT2(3) | |
3680 | DO 1017 J = 1, NTJ(I) | |
3681 | nsubp = nsubp + 1 | |
3682 | IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1 | |
3683 | 1017 CONTINUE | |
3684 | 1018 CONTINUE | |
3685 | ||
3686 | DO 1020 I = 1, NSG | |
3687 | DO 1019 J = 1, NJSG(I) | |
3688 | nsubp = nsubp + 1 | |
3689 | IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1 | |
3690 | 1019 CONTINUE | |
3691 | 1020 CONTINUE | |
3692 | NISG = NISG + NSG | |
3693 | IF (IOUT .EQ. 1) THEN | |
3694 | cbzdbg2/16/99 | |
3695 | c PRINT *, ' in HJANA1 ' | |
3696 | c PRINT *, ' total number of partons = ', nsubp | |
3697 | c PRINT *, ' total number of gluons = ', nsubg, MUL | |
3698 | c PRINT *, ' number of projectile strings = ', IHNT2(1) | |
3699 | c PRINT *, ' number of target strings = ', IHNT2(3) | |
3700 | c PRINT *, ' number of independent strings = ', NSG | |
3701 | PRINT *, ' in HJANA1 ' | |
3702 | PRINT *, ' total number of partons = ', nsubp / IW | |
3703 | PRINT *, ' total number of gluons = ', nsubg / IW | |
3704 | c PRINT *, ' number of projectile strings = ', IHNT2(1) | |
3705 | c PRINT *, ' number of target strings = ', IHNT2(3) | |
3706 | PRINT *, ' number of independent strings = ', NISG / IW | |
3707 | cbzdbg2/16/99end | |
3708 | END IF | |
3709 | c | |
3710 | RETURN | |
3711 | END | |
3712 | ||
3713 | c----------------------------------------------------------------------- | |
3714 | ||
3715 | c.....analysis subroutine in ZPC after generation of additional initial | |
3716 | c.....phase space distributions. | |
3717 | ||
3718 | SUBROUTINE HJAN1A | |
3719 | PARAMETER (MAXPTN=400001) | |
3720 | PARAMETER (DGX = 0.2, DGY = 0.2, DT = 0.2) | |
3721 | DIMENSION dgxg1a(50), dgyg1a(50), dtg1a(50) | |
3722 | DIMENSION sgxg1a(50), sgyg1a(50), stg1a(50) | |
3723 | DOUBLE PRECISION GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5 | |
3724 | COMMON /PARA1/ MUL | |
3725 | cc SAVE /PARA1/ | |
3726 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
3727 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
3728 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
3729 | cc SAVE /prec2/ | |
3730 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
3731 | cc SAVE /AREVT/ | |
3732 | COMMON /AROUT/ IOUT | |
3733 | cc SAVE /AROUT/ | |
3734 | SAVE | |
3735 | DATA IW/0/ | |
3736 | ||
3737 | IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN | |
3738 | DO 1001 I = 1, 50 | |
3739 | dgxg1a(I) = sgxg1a(I) | |
3740 | dgyg1a(I) = sgyg1a(I) | |
3741 | dtg1a(I) = stg1a(I) | |
3742 | 1001 CONTINUE | |
3743 | ELSE | |
3744 | DO 1002 I = 1, 50 | |
3745 | sgxg1a(I) = dgxg1a(I) | |
3746 | sgyg1a(I) = dgyg1a(I) | |
3747 | stg1a(I) = dtg1a(I) | |
3748 | 1002 CONTINUE | |
3749 | isevt = IAEVT | |
3750 | isrun = IARUN | |
3751 | IW = IW + 1 | |
3752 | END IF | |
3753 | c.....analysis | |
3754 | DO 1003 I = 1, MUL | |
3755 | IGX = 1 + int(sngl(ABS(GX5(I))) / DGX) | |
3756 | clin-4/2008 protect against out-of-bound errors: | |
3757 | c IF (IGX .GT. 50) GOTO 100 | |
3758 | IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100 | |
3759 | dgxg1a(IGX) = dgxg1a(IGX) + 1.0 | |
3760 | 100 CONTINUE | |
3761 | IGY = 1 + int(sngl(ABS(GY5(I))) / DGY) | |
3762 | IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 200 | |
3763 | dgyg1a(IGY) = dgyg1a(IGY) + 1.0 | |
3764 | 200 CONTINUE | |
3765 | IT = 1 + int(sngl(SQRT(FT5(I) ** 2 - GZ5(I) ** 2)) / DT) | |
3766 | IF (IT .GT. 50 .or. IT .LT. 1) GOTO 300 | |
3767 | dtg1a(IT) = dtg1a(IT) + 1.0 | |
3768 | 300 CONTINUE | |
3769 | 1003 CONTINUE | |
3770 | CALL HJAN1B | |
3771 | c | |
3772 | RETURN | |
3773 | END | |
3774 | ||
3775 | c----------------------------------------------------------------------- | |
3776 | ||
3777 | c.....analysis subroutine in HJAN1A | |
3778 | ||
3779 | SUBROUTINE HJAN1B | |
3780 | PARAMETER (MAXPTN=400001,MAXSTR=150001) | |
3781 | PARAMETER (DR = 0.2, DT = 0.2) | |
3782 | DIMENSION DNRG1B(50), dtg1b(50) | |
3783 | DIMENSION SNRG1B(50), stg1b(50) | |
3784 | DOUBLE PRECISION GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5 | |
3785 | COMMON /PARA1/ MUL | |
3786 | cc SAVE /PARA1/ | |
3787 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
3788 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
3789 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
3790 | cc SAVE /prec2/ | |
3791 | COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN) | |
3792 | cc SAVE /ilist8/ | |
3793 | COMMON /SREC1/ NSP, NST, NSI | |
3794 | cc SAVE /SREC1/ | |
3795 | COMMON/hjcrdn/YP(3,300),YT(3,300) | |
3796 | cc SAVE /hjcrdn/ | |
3797 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
3798 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
3799 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
3800 | cc SAVE /HJJET2/ | |
3801 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
3802 | cc SAVE /AREVT/ | |
3803 | COMMON /AROUT/ IOUT | |
3804 | cc SAVE /AROUT/ | |
3805 | SAVE | |
3806 | DATA IW/0/ | |
3807 | ||
3808 | IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN | |
3809 | DO 1001 I = 1, 50 | |
3810 | DNRG1B(I) = SNRG1B(I) | |
3811 | dtg1b(I) = stg1b(I) | |
3812 | 1001 CONTINUE | |
3813 | ELSE | |
3814 | DO 1002 I = 1, 50 | |
3815 | SNRG1B(I) = DNRG1B(I) | |
3816 | stg1b(I) = dtg1b(I) | |
3817 | 1002 CONTINUE | |
3818 | isevt = IAEVT | |
3819 | isrun = IARUN | |
3820 | IW = IW + 1 | |
3821 | END IF | |
3822 | c.....analysis | |
3823 | DO 1003 I = 1, MUL | |
3824 | J = LSTRG1(I) | |
3825 | ||
3826 | IF (J .LE. NSP) THEN | |
3827 | K = J | |
3828 | GX0 = YP(1, J) | |
3829 | GY0 = YP(2, J) | |
3830 | ELSE IF (J .LE. NSP + NST) THEN | |
3831 | K = J - NSP | |
3832 | GX0 = YT(1, K) | |
3833 | GY0 = YT(2, K) | |
3834 | ELSE | |
3835 | K = J - NSP - NST | |
3836 | GX0 = 0.5 * (YP(1, IASG(K, 1)) + YT(1, IASG(K, 2))) | |
3837 | GY0 = 0.5 * (YP(2, IASG(K, 1)) + YT(2, IASG(K, 2))) | |
3838 | END IF | |
3839 | R0 = SQRT((sngl(GX5(I)) - GX0)**2 + (sngl(GY5(I)) - GY0)**2) | |
3840 | IR = 1 + int(R0 / DR) | |
3841 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 100 | |
3842 | DNRG1B(IR) = DNRG1B(IR) + 1.0 | |
3843 | 100 CONTINUE | |
3844 | TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2)) | |
3845 | IT = 1 + int(TAU7 / DT) | |
3846 | IF (IT .GT. 50 .or. IT .LT. 1) GOTO 200 | |
3847 | dtg1b(IT) = dtg1b(IT) + 1.0 | |
3848 | 200 CONTINUE | |
3849 | 1003 CONTINUE | |
3850 | c | |
3851 | RETURN | |
3852 | END | |
3853 | ||
3854 | c----------------------------------------------------------------------- | |
3855 | ||
3856 | c.....analysis subroutine in HIJING after parton cascade evolution | |
3857 | SUBROUTINE HJANA2 | |
3858 | c | |
3859 | PARAMETER (YMAX = 1.0, YMIN = -1.0) | |
3860 | PARAMETER (DMT = 0.05, DY = 0.2) | |
3861 | PARAMETER (DR = 0.2, DT = 0.2) | |
3862 | PARAMETER (MAXPTN=400001) | |
3863 | PARAMETER (MAXSTR=150001) | |
3864 | DOUBLE PRECISION PXSGS,PYSGS,PZSGS,PESGS,PMSGS, | |
3865 | 1 GXSGS,GYSGS,GZSGS,FTSGS | |
3866 | DIMENSION dyp2(50), DMYP2(200), DEYP2(50) | |
3867 | DIMENSION dyg2(50), DMYG2(200), DEYG2(50) | |
3868 | DIMENSION SNYP2(50), SMYP2(200), SEYP2(50) | |
3869 | DIMENSION SNYG2(50), SMYG2(200), SEYG2(50) | |
3870 | DIMENSION dnrpj2(50), dnrtg2(50), dnrin2(50), | |
3871 | & dnrtt2(50) | |
3872 | DIMENSION dtpj2(50), dttg2(50), dtin2(50), | |
3873 | & dttot2(50) | |
3874 | DIMENSION dyg2c(50), dmyg2c(50), deyg2c(50) | |
3875 | DIMENSION snrpj2(50), snrtg2(50), snrin2(50), | |
3876 | & snrtt2(50) | |
3877 | DIMENSION stpj2(50), sttg2(50), stin2(50), | |
3878 | & sttot2(50) | |
3879 | DIMENSION snyg2c(50), smyg2c(50), seyg2c(50) | |
3880 | DOUBLE PRECISION ATAUI, ZT1, ZT2, ZT3 | |
3881 | DOUBLE PRECISION GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5 | |
3882 | COMMON /PARA1/ MUL | |
3883 | cc SAVE /PARA1/ | |
3884 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
3885 | cc SAVE /HPARNT/ | |
3886 | COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR) | |
3887 | cc SAVE /SREC2/ | |
3888 | COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500), | |
3889 | & PJPY(300,500),PJPZ(300,500),PJPE(300,500), | |
3890 | & PJPM(300,500),NTJ(300),KFTJ(300,500), | |
3891 | & PJTX(300,500),PJTY(300,500),PJTZ(300,500), | |
3892 | & PJTE(300,500),PJTM(300,500) | |
3893 | cc SAVE /HJJET1/ | |
3894 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
3895 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
3896 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
3897 | cc SAVE /HJJET2/ | |
3898 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
3899 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
3900 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
3901 | cc SAVE /prec2/ | |
3902 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
3903 | cc SAVE /AREVT/ | |
3904 | COMMON /AROUT/ IOUT | |
3905 | cc SAVE /AROUT/ | |
3906 | common/anim/nevent,isoft,isflag,izpc | |
3907 | cc SAVE /anim/ | |
3908 | COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3), | |
3909 | & PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3), | |
3910 | & GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3), | |
3911 | & K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR) | |
3912 | cc SAVE /SOFT/ | |
3913 | SAVE | |
3914 | DATA IW/0/ | |
3915 | ||
3916 | IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN | |
3917 | DO 1001 I = 1, 200 | |
3918 | DMYP2(I) = SMYP2(I) | |
3919 | DMYG2(I) = SMYG2(I) | |
3920 | 1001 CONTINUE | |
3921 | ||
3922 | DO 1002 I = 1, 50 | |
3923 | dyp2(I) = SNYP2(I) | |
3924 | DEYP2(I) = SEYP2(I) | |
3925 | dyg2(I) = SNYG2(I) | |
3926 | DEYG2(I) = SEYG2(I) | |
3927 | dnrpj2(I) = snrpj2(I) | |
3928 | dnrtg2(I) = snrtg2(I) | |
3929 | dnrin2(I) = snrin2(I) | |
3930 | dnrtt2(I) = snrtt2(I) | |
3931 | dtpj2(I) = stpj2(I) | |
3932 | dttg2(I) = sttg2(I) | |
3933 | dtin2(I) = stin2(I) | |
3934 | dttot2(I) = sttot2(I) | |
3935 | dyg2c(I) = snyg2c(I) | |
3936 | dmyg2c(I) = smyg2c(I) | |
3937 | deyg2c(I) = seyg2c(I) | |
3938 | 1002 CONTINUE | |
3939 | nsubp = nsubpS | |
3940 | nsubg = nsubgS | |
3941 | NISG = NISGS | |
3942 | ELSE | |
3943 | DO 1003 I = 1, 200 | |
3944 | SMYP2(I) = DMYP2(I) | |
3945 | SMYG2(I) = DMYG2(I) | |
3946 | 1003 CONTINUE | |
3947 | ||
3948 | DO 1004 I = 1, 50 | |
3949 | SNYP2(I) = dyp2(I) | |
3950 | SEYP2(I) = DEYP2(I) | |
3951 | SNYG2(I) = dyg2(I) | |
3952 | SEYG2(I) = DEYG2(I) | |
3953 | snrpj2(I) = dnrpj2(I) | |
3954 | snrtg2(I) = dnrtg2(I) | |
3955 | snrin2(I) = dnrin2(I) | |
3956 | snrtt2(I) = dnrtt2(I) | |
3957 | stpj2(I) = dtpj2(I) | |
3958 | sttg2(I) = dttg2(I) | |
3959 | stin2(I) = dtin2(I) | |
3960 | sttot2(I) = dttot2(I) | |
3961 | snyg2c(I) = dyg2c(I) | |
3962 | smyg2c(I) = dmyg2c(I) | |
3963 | seyg2c(I) = deyg2c(I) | |
3964 | 1004 CONTINUE | |
3965 | nsubpS = nsubp | |
3966 | nsubgS = nsubg | |
3967 | NISGS = NISG | |
3968 | isevt = IAEVT | |
3969 | isrun = IARUN | |
3970 | IW = IW + 1 | |
3971 | END IF | |
3972 | ||
3973 | clin-4/28/01: | |
3974 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 510 | |
3975 | ||
3976 | c.....analysis | |
3977 | DO 1006 I = 1, IHNT2(1) | |
3978 | DO 1005 J = 1, NPJ(I) | |
3979 | ITYP = KFPJ(I, J) | |
3980 | PX = PJPX(I, J) | |
3981 | PY = PJPY(I, J) | |
3982 | PZ = PJPZ(I, J) | |
3983 | PE = PJPE(I, J) | |
3984 | PM = PJPM(I, J) | |
3985 | cbzdbg2/16/99 | |
3986 | c IF (ABS(PZ) .GE. PE) GOTO 200 | |
3987 | IF (ABS(PZ) .GE. PE) THEN | |
3988 | PRINT *, ' IN HJANA2, PROJ STR ', I, ' PART ', J | |
3989 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
3990 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
3991 | PRINT *, ' XM = ', PM | |
3992 | GOTO 200 | |
3993 | END IF | |
3994 | cbzdbg2/16/99end | |
3995 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
3996 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
3997 | DXMT = XMT - PM | |
3998 | IY = 1 + int(ABS(RAP) / DY) | |
3999 | IF (IY .GT. 50) GOTO 100 | |
4000 | dyp2(IY) = dyp2(IY) + 1.0 | |
4001 | DEYP2(IY) = DEYP2(IY) + XMT | |
4002 | IF (ITYP .EQ. 21) THEN | |
4003 | dyg2(IY) = dyg2(IY) + 1.0 | |
4004 | DEYG2(IY) = DEYG2(IY) + XMT | |
4005 | END IF | |
4006 | 100 CONTINUE | |
4007 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 200 | |
4008 | IMT = 1 + int(DXMT / DMT) | |
4009 | IF (IMT .GT. 200) GOTO 200 | |
4010 | DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT | |
4011 | IF (ITYP .EQ. 21) THEN | |
4012 | DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT | |
4013 | END IF | |
4014 | 200 CONTINUE | |
4015 | 1005 CONTINUE | |
4016 | 1006 CONTINUE | |
4017 | ||
4018 | DO 1008 I = 1, IHNT2(3) | |
4019 | DO 1007 J = 1, NTJ(I) | |
4020 | ITYP = KFTJ(I, J) | |
4021 | PX = PJTX(I, J) | |
4022 | PY = PJTY(I, J) | |
4023 | PZ = PJTZ(I, J) | |
4024 | PE = PJTE(I, J) | |
4025 | PM = PJTM(I, J) | |
4026 | cbzdbg2/16/99 | |
4027 | c IF (ABS(PZ) .GE. PE) GOTO 400 | |
4028 | IF (ABS(PZ) .GE. PE) THEN | |
4029 | PRINT *, ' IN HJANA2, TARG STR ', I, ' PART ', J | |
4030 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
4031 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
4032 | PRINT *, ' XM = ', PM | |
4033 | GOTO 400 | |
4034 | END IF | |
4035 | cbzdbg2/16/99end | |
4036 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
4037 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
4038 | DXMT = XMT - PM | |
4039 | IY = 1 + int(ABS(RAP) / DY) | |
4040 | IF (IY .GT. 50) GOTO 300 | |
4041 | dyp2(IY) = dyp2(IY) + 1.0 | |
4042 | DEYP2(IY) = DEYP2(IY) + XMT | |
4043 | IF (ITYP .EQ. 21) THEN | |
4044 | dyg2(IY) = dyg2(IY) + 1.0 | |
4045 | DEYG2(IY) = DEYG2(IY) + XMT | |
4046 | END IF | |
4047 | 300 CONTINUE | |
4048 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 400 | |
4049 | IMT = 1 + int(DXMT / DMT) | |
4050 | IF (IMT .GT. 200) GOTO 400 | |
4051 | DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT | |
4052 | IF (ITYP .EQ. 21) THEN | |
4053 | DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT | |
4054 | END IF | |
4055 | 400 CONTINUE | |
4056 | 1007 CONTINUE | |
4057 | 1008 CONTINUE | |
4058 | ||
4059 | clin-4/28/01: | |
4060 | 510 continue | |
4061 | ||
4062 | DO 1010 I = 1, NSG | |
4063 | clin-4/25/01 soft3: | |
4064 | c DO J = 1, NJSG(I) | |
4065 | NJ=NJSG(I) | |
4066 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I) | |
4067 | DO 1009 J = 1, NJ | |
4068 | clin-4/25/01-end | |
4069 | ||
4070 | ITYP = K2SG(I, J) | |
4071 | PX = PXSG(I, J) | |
4072 | PY = PYSG(I, J) | |
4073 | PZ = PZSG(I, J) | |
4074 | PE = PESG(I, J) | |
4075 | PM = PMSG(I, J) | |
4076 | clin-4/25/01 soft3: | |
4077 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then | |
4078 | ITYP = K2SGS(I, J) | |
4079 | PX = sngl(PXSGS(I, J)) | |
4080 | PY = sngl(PYSGS(I, J)) | |
4081 | PZ = sngl(PZSGS(I, J)) | |
4082 | PE = sngl(PESGS(I, J)) | |
4083 | PM = sngl(PMSGS(I, J)) | |
4084 | endif | |
4085 | clin-4/25/01-end | |
4086 | ||
4087 | cbzdbg2/16/99 | |
4088 | c IF (ABS(PZ) .GE. PE) GOTO 600 | |
4089 | IF (ABS(PZ) .GE. PE) THEN | |
4090 | PRINT *, ' IN HJANA2, INDP STR ', I, ' PART ', J | |
4091 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
4092 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
4093 | PRINT *, ' XM = ', PM | |
4094 | GOTO 600 | |
4095 | END IF | |
4096 | cbzdbg2/16/99end | |
4097 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
4098 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
4099 | DXMT = XMT - PM | |
4100 | IY = 1 + int(ABS(RAP) / DY) | |
4101 | IF (IY .GT. 50) GOTO 500 | |
4102 | dyp2(IY) = dyp2(IY) + 1.0 | |
4103 | DEYP2(IY) = DEYP2(IY) + XMT | |
4104 | IF (ITYP .EQ. 21) THEN | |
4105 | dyg2(IY) = dyg2(IY) + 1.0 | |
4106 | DEYG2(IY) = DEYG2(IY) + XMT | |
4107 | END IF | |
4108 | 500 CONTINUE | |
4109 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 600 | |
4110 | IMT = 1 + int(DXMT / DMT) | |
4111 | IF (IMT .GT. 200) GOTO 600 | |
4112 | DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT | |
4113 | IF (ITYP .EQ. 21) THEN | |
4114 | DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT | |
4115 | END IF | |
4116 | 600 CONTINUE | |
4117 | 1009 CONTINUE | |
4118 | 1010 CONTINUE | |
4119 | ||
4120 | clin-4/28/01: | |
4121 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 520 | |
4122 | ||
4123 | DO 1011 I = 1, IHNT2(1) | |
4124 | J = I | |
4125 | YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2)) | |
4126 | IR = 1 + int(YR / DR) | |
4127 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601 | |
4128 | dnrpj2(IR) = dnrpj2(IR) + 1.0 | |
4129 | dnrtt2(IR) = dnrtt2(IR) + 1.0 | |
4130 | 601 CONTINUE | |
4131 | IT = 1 + int(sngl(ATAUI(J)) / DT) | |
4132 | IF (IT .GT. 50 .or. IT .LT. 1) GOTO 602 | |
4133 | dtpj2(IT) = dtpj2(IT) + 1.0 | |
4134 | dttot2(IT) = dttot2(IT) + 1.0 | |
4135 | 602 CONTINUE | |
4136 | 1011 CONTINUE | |
4137 | ||
4138 | DO 1012 I = 1, IHNT2(3) | |
4139 | J = I + IHNT2(1) | |
4140 | YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2)) | |
4141 | IR = 1 + int(YR / DR) | |
4142 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 603 | |
4143 | dnrtg2(IR) = dnrtg2(IR) + 1.0 | |
4144 | dnrtt2(IR) = dnrtt2(IR) + 1.0 | |
4145 | 603 CONTINUE | |
4146 | IT = 1 + int(sngl(ATAUI(J)) / DT) | |
4147 | IF (IT .GT. 50 .or. IT .LT. 1) GOTO 604 | |
4148 | dttg2(IT) = dttg2(IT) + 1.0 | |
4149 | dttot2(IT) = dttot2(IT) + 1.0 | |
4150 | 604 CONTINUE | |
4151 | 1012 CONTINUE | |
4152 | ||
4153 | clin-4/28/01: | |
4154 | 520 continue | |
4155 | ||
4156 | DO 1013 I = 1, NSG | |
4157 | J = I + IHNT2(1) + IHNT2(3) | |
4158 | clin-4/28/01: | |
4159 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) J = I | |
4160 | ||
4161 | YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2)) | |
4162 | IR = 1 + int(YR / DR) | |
4163 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 605 | |
4164 | dnrin2(IR) = dnrin2(IR) + 1.0 | |
4165 | dnrtt2(IR) = dnrtt2(IR) + 1.0 | |
4166 | 605 CONTINUE | |
4167 | IT = 1 + int(sngl(ATAUI(J)) / DT) | |
4168 | IF (IT .GT. 50 .or. IT .LT. 1) GOTO 606 | |
4169 | dtin2(IT) = dtin2(IT) + 1.0 | |
4170 | dttot2(IT) = dttot2(IT) + 1.0 | |
4171 | 606 CONTINUE | |
4172 | 1013 CONTINUE | |
4173 | ||
4174 | DO 1014 I = 1, MUL | |
4175 | ITYP = ITYP5(I) | |
4176 | PX = sngl(PX5(I)) | |
4177 | PY = sngl(PY5(I)) | |
4178 | PZ = sngl(PZ5(I)) | |
4179 | PE = sngl(E5(I)) | |
4180 | PM = sngl(XMASS5(I)) | |
4181 | cbzdbg2/16/99 | |
4182 | c IF (ABS(PZ) .GE. PE) GOTO 800 | |
4183 | ||
4184 | IF (ABS(PZ) .GE. PE) THEN | |
4185 | PRINT *, ' IN HJANA2, GLUON ', I | |
4186 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
4187 | PRINT *, ' PZ = ', PZ, ' EE = ', PE | |
4188 | PRINT *, ' XM = ', PM | |
4189 | GOTO 800 | |
4190 | END IF | |
4191 | ||
4192 | cbzdbg2/16/99end | |
4193 | RAP = 0.5 * LOG((PE + PZ) / (PE - PZ)) | |
4194 | XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2) | |
4195 | DXMT = XMT - PM | |
4196 | IY = 1 + int(ABS(RAP) / DY) | |
4197 | IF (IY .GT. 50) GOTO 700 | |
4198 | dyg2c(IY) = dyg2c(IY) + 1.0 | |
4199 | deyg2c(IY) = deyg2c(IY) + XMT | |
4200 | 700 CONTINUE | |
4201 | IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 800 | |
4202 | IMT = 1 + int(DXMT / DMT) | |
4203 | IF (IMT .GT. 50) GOTO 800 | |
4204 | dmyg2c(IMT) = dmyg2c(IMT) + 1.0 / XMT | |
4205 | 800 CONTINUE | |
4206 | 1014 CONTINUE | |
4207 | ||
4208 | clin-4/25/01 soft3: | |
4209 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 530 | |
4210 | ||
4211 | c.....count number of particles | |
4212 | DO 1016 I = 1, IHNT2(1) | |
4213 | DO 1015 J = 1, NPJ(I) | |
4214 | nsubp = nsubp + 1 | |
4215 | IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1 | |
4216 | 1015 CONTINUE | |
4217 | 1016 CONTINUE | |
4218 | ||
4219 | DO 1018 I = 1, IHNT2(3) | |
4220 | DO 1017 J = 1, NTJ(I) | |
4221 | nsubp = nsubp + 1 | |
4222 | IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1 | |
4223 | 1017 CONTINUE | |
4224 | 1018 CONTINUE | |
4225 | ||
4226 | clin-4/25/01 soft3: | |
4227 | 530 continue | |
4228 | ||
4229 | DO 1020 I = 1, NSG | |
4230 | clin-4/25/01 soft3: | |
4231 | c DO J = 1, NJSG(I) | |
4232 | NJ=NJSG(I) | |
4233 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I) | |
4234 | DO 1019 J = 1, NJ | |
4235 | clin-4/25/01-end | |
4236 | ||
4237 | nsubp = nsubp + 1 | |
4238 | ||
4239 | clin-4/25/01 | |
4240 | c IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1 | |
4241 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then | |
4242 | IF(K2SGS(I, J) .EQ. 21) nsubg = nsubg + 1 | |
4243 | else | |
4244 | IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1 | |
4245 | endif | |
4246 | clin-4/25/01-end | |
4247 | 1019 CONTINUE | |
4248 | 1020 CONTINUE | |
4249 | cbzdbg2/16/99 | |
4250 | NISG = NISG + NSG | |
4251 | ||
4252 | IF (IOUT .EQ. 1) THEN | |
4253 | cbzdbg2/16/99end | |
4254 | cbzdbg2/16/99 | |
4255 | c PRINT *, ' in HJANA2 ' | |
4256 | c PRINT *, ' total number of partons = ', nsubp | |
4257 | c PRINT *, ' total number of gluons = ', nsubg, MUL | |
4258 | c PRINT *, ' number of projectile strings = ', IHNT2(1) | |
4259 | c PRINT *, ' number of target strings = ', IHNT2(3) | |
4260 | c PRINT *, ' number of independent strings = ', NSG | |
4261 | PRINT *, ' in HJANA2 ' | |
4262 | PRINT *, ' total number of partons = ', nsubp / IW | |
4263 | PRINT *, ' total number of gluons = ', nsubg / IW | |
4264 | c PRINT *, ' number of projectile strings = ', IHNT2(1) | |
4265 | c PRINT *, ' number of target strings = ', IHNT2(3) | |
4266 | PRINT *, ' number of independent strings = ', NISG / IW | |
4267 | END IF | |
4268 | ||
4269 | CALL HJAN2A | |
4270 | CALL HJAN2B | |
4271 | ||
4272 | RETURN | |
4273 | END | |
4274 | ||
4275 | c----------------------------------------------------------------------- | |
4276 | ||
4277 | c.....subroutine called by HJANA2 | |
4278 | SUBROUTINE HJAN2A | |
4279 | ||
4280 | PARAMETER (DGX = 0.2, DGY = 0.2, DT = 0.2) | |
4281 | PARAMETER (MAXPTN=400001,MAXSTR=150001) | |
4282 | DIMENSION dgxp2a(50), dgyp2a(50), dtp2a(50) | |
4283 | DIMENSION dgxg2a(50), dgyg2a(50), dtg2a(50) | |
4284 | DIMENSION sgxp2a(50), sgyp2a(50), stp2a(50) | |
4285 | DIMENSION sgxg2a(50), sgyg2a(50), stg2a(50) | |
4286 | DOUBLE PRECISION GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5 | |
4287 | COMMON /PARA1/ MUL | |
4288 | cc SAVE /PARA1/ | |
4289 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
4290 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
4291 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
4292 | cc SAVE /prec2/ | |
4293 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
4294 | cc SAVE /HPARNT/ | |
4295 | COMMON/hjcrdn/YP(3,300),YT(3,300) | |
4296 | cc SAVE /hjcrdn/ | |
4297 | COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500), | |
4298 | & PJPY(300,500),PJPZ(300,500),PJPE(300,500), | |
4299 | & PJPM(300,500),NTJ(300),KFTJ(300,500), | |
4300 | & PJTX(300,500),PJTY(300,500),PJTZ(300,500), | |
4301 | & PJTE(300,500),PJTM(300,500) | |
4302 | cc SAVE /HJJET1/ | |
4303 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
4304 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
4305 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
4306 | cc SAVE /HJJET2/ | |
4307 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
4308 | cc SAVE /AREVT/ | |
4309 | COMMON /AROUT/ IOUT | |
4310 | cc SAVE /AROUT/ | |
4311 | SAVE | |
4312 | DATA IW/0/ | |
4313 | ||
4314 | IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN | |
4315 | DO 1001 I = 1, 50 | |
4316 | dgxp2a(I) = sgxp2a(I) | |
4317 | dgyp2a(I) = sgyp2a(I) | |
4318 | dtp2a(I) = stp2a(I) | |
4319 | dgxg2a(I) = sgxg2a(I) | |
4320 | dgyg2a(I) = sgyg2a(I) | |
4321 | dtg2a(I) = stg2a(I) | |
4322 | 1001 CONTINUE | |
4323 | ELSE | |
4324 | DO 1002 I = 1, 50 | |
4325 | sgxp2a(I) = dgxp2a(I) | |
4326 | sgyp2a(I) = dgyp2a(I) | |
4327 | stp2a(I) = dtp2a(I) | |
4328 | sgxg2a(I) = dgxg2a(I) | |
4329 | sgyg2a(I) = dgyg2a(I) | |
4330 | stg2a(I) = dtg2a(I) | |
4331 | 1002 CONTINUE | |
4332 | isevt = IAEVT | |
4333 | isrun = IARUN | |
4334 | IW = IW + 1 | |
4335 | END IF | |
4336 | c.....analysis | |
4337 | DO 1004 I = 1, IHNT2(1) | |
4338 | DO 1003 J = 1, NPJ(I) | |
4339 | IF (KFPJ(I, J) .NE. 21) THEN | |
4340 | IGX = 1 + int(ABS(YP(1, I)) / DGX) | |
4341 | IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100 | |
4342 | dgxp2a(IGX) = dgxp2a(IGX) + 1.0 | |
4343 | 100 CONTINUE | |
4344 | IGY = 1 + int(ABS(YP(2, I)) / DGY) | |
4345 | IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 200 | |
4346 | dgyp2a(IGY) = dgyp2a(IGY) + 1.0 | |
4347 | 200 CONTINUE | |
4348 | IT = 1 | |
4349 | dtp2a(IT) = dtp2a(IT) + 1.0 | |
4350 | END IF | |
4351 | 1003 CONTINUE | |
4352 | 1004 CONTINUE | |
4353 | ||
4354 | DO 1006 I = 1, IHNT2(3) | |
4355 | DO 1005 J = 1, NTJ(I) | |
4356 | IF (KFTJ(I, J) .NE. 21) THEN | |
4357 | IGX = 1 + int(ABS(YT(1, I)) / DGX) | |
4358 | IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 300 | |
4359 | dgxp2a(IGX) = dgxp2a(IGX) + 1.0 | |
4360 | 300 CONTINUE | |
4361 | IGY = 1 + int(ABS(YT(2, I)) / DGY) | |
4362 | IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 400 | |
4363 | dgyp2a(IGY) = dgyp2a(IGY) + 1.0 | |
4364 | 400 CONTINUE | |
4365 | IT = 1 | |
4366 | dtp2a(IT) = dtp2a(IT) + 1.0 | |
4367 | END IF | |
4368 | 1005 CONTINUE | |
4369 | 1006 CONTINUE | |
4370 | ||
4371 | DO 1008 I = 1, NSG | |
4372 | DO 1007 J = 1, NJSG(I) | |
4373 | IF (K2SG(I, J) .NE. 21) THEN | |
4374 | IGX = 1 + int(ABS(0.5 * | |
4375 | & (YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))) / DGX) | |
4376 | IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 500 | |
4377 | dgxp2a(IGX) = dgxp2a(IGX) + 1.0 | |
4378 | 500 CONTINUE | |
4379 | IGY = 1 + int(ABS(0.5 * | |
4380 | & (YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))) / DGY) | |
4381 | IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 600 | |
4382 | dgyp2a(IGY) = dgyp2a(IGY) + 1.0 | |
4383 | 600 CONTINUE | |
4384 | IT = 1 | |
4385 | dtp2a(IT) = dtp2a(IT) + 1.0 | |
4386 | END IF | |
4387 | 1007 CONTINUE | |
4388 | 1008 CONTINUE | |
4389 | ||
4390 | DO 1009 I = 1, MUL | |
4391 | IGX = 1 + int(ABS(sngl(GX5(I))) / DGX) | |
4392 | IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 700 | |
4393 | dgxg2a(IGX) = dgxg2a(IGX) + 1.0 | |
4394 | dgxp2a(IGX) = dgxp2a(IGX) + 1.0 | |
4395 | 700 CONTINUE | |
4396 | IGY = 1 + int(ABS(sngl(GY5(I))) / DGY) | |
4397 | IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 800 | |
4398 | dgyg2a(IGY) = dgyg2a(IGY) + 1.0 | |
4399 | dgyp2a(IGY) = dgyp2a(IGY) + 1.0 | |
4400 | 800 CONTINUE | |
4401 | IT = 1 + int(SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2)) / DT) | |
4402 | IF (IT .GT. 50 .or. IT .LT. 1) GOTO 900 | |
4403 | dtg2a(IT) = dtg2a(IT) + 1.0 | |
4404 | dtp2a(IT) = dtp2a(IT) + 1.0 | |
4405 | 900 CONTINUE | |
4406 | 1009 CONTINUE | |
4407 | c | |
4408 | RETURN | |
4409 | END | |
4410 | ||
4411 | c----------------------------------------------------------------------- | |
4412 | ||
4413 | c.....analysis subroutine in HJANA2 | |
4414 | ||
4415 | SUBROUTINE HJAN2B | |
4416 | ||
4417 | PARAMETER (MAXPTN=400001) | |
4418 | PARAMETER (MAXSTR=150001) | |
4419 | PARAMETER (DR = 0.2, DT = 0.2) | |
4420 | DIMENSION DNRG2B(50), dtg2b(-24:25) | |
4421 | DIMENSION SNRG2B(50), stg2b(-24:25) | |
4422 | DOUBLE PRECISION GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5 | |
4423 | DOUBLE PRECISION ATAUI, ZT1, ZT2, ZT3 | |
4424 | COMMON /PARA1/ MUL | |
4425 | cc SAVE /PARA1/ | |
4426 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
4427 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
4428 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
4429 | cc SAVE /prec2/ | |
4430 | COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN) | |
4431 | cc SAVE /ilist8/ | |
4432 | COMMON /SREC1/ NSP, NST, NSI | |
4433 | cc SAVE /SREC1/ | |
4434 | COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR) | |
4435 | cc SAVE /SREC2/ | |
4436 | COMMON/hjcrdn/YP(3,300),YT(3,300) | |
4437 | cc SAVE /hjcrdn/ | |
4438 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
4439 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
4440 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
4441 | cc SAVE /HJJET2/ | |
4442 | COMMON /AREVT/ IAEVT, IARUN, MISS | |
4443 | cc SAVE /AREVT/ | |
4444 | COMMON /AROUT/ IOUT | |
4445 | cc SAVE /AROUT/ | |
4446 | SAVE | |
4447 | DATA IW/0/ | |
4448 | ||
4449 | IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN | |
4450 | DO 1001 I = 1, 50 | |
4451 | DNRG2B(I) = SNRG2B(I) | |
4452 | dtg2b(I - 25) = stg2b(I - 25) | |
4453 | 1001 CONTINUE | |
4454 | ELSE | |
4455 | DO 1002 I = 1, 50 | |
4456 | SNRG2B(I) = DNRG2B(I) | |
4457 | stg2b(I - 25) = dtg2b(I - 25) | |
4458 | 1002 CONTINUE | |
4459 | isevt = IAEVT | |
4460 | isrun = IARUN | |
4461 | IW = IW + 1 | |
4462 | END IF | |
4463 | c.....analysis | |
4464 | DO 1003 I = 1, MUL | |
4465 | J = LSTRG1(I) | |
4466 | GX0 = sngl(ZT1(J)) | |
4467 | GY0 = sngl(ZT2(J)) | |
4468 | R0 = SQRT((sngl(GX5(I)) - GX0)**2 + (sngl(GY5(I)) - GY0)**2) | |
4469 | IR = 1 + int(R0 / DR) | |
4470 | IF (IR .GT. 50 .or. IR .LT. 1) GOTO 100 | |
4471 | DNRG2B(IR) = DNRG2B(IR) + 1.0 | |
4472 | 100 CONTINUE | |
4473 | TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2)) | |
4474 | DTAU=TAU7 - sngl(ATAUI(J)) | |
4475 | IT = 1 + int(DTAU / DT) | |
4476 | cbzdbg2/21/99 | |
4477 | c IF (ABS(IT) .GT. 25) GOTO 200 | |
4478 | IF (IT .GT. 25 .OR. IT .LT. -24) GOTO 200 | |
4479 | cbzdbg2/21/99end | |
4480 | dtg2b(IT) = dtg2b(IT) + 1.0 | |
4481 | 200 CONTINUE | |
4482 | 1003 CONTINUE | |
4483 | c | |
4484 | RETURN | |
4485 | END | |
4486 | ||
4487 | c----------------------------------------------------------------------- | |
4488 | ||
4489 | c.....analysis subroutine before ARTMN | |
4490 | SUBROUTINE HJANA3 | |
4491 | c | |
4492 | PARAMETER (MAXSTR=150001, MAXR=1) | |
4493 | c.....y cut for mt spectrum | |
4494 | PARAMETER (YMIN = -1.0, YMAX = 1.0) | |
4495 | cbz11/7/99 end | |
4496 | c.....bin width for mt spectrum and y spectrum | |
4497 | PARAMETER (DMT = 0.05, DY = 0.2) | |
4498 | DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult, | |
4499 | 1 v2bi,s2bi,xbmult | |
4500 | DIMENSION dndyh3(50), DMYH3(50), DEYH3(50) | |
4501 | COMMON /RUN/ NUM | |
4502 | cc SAVE /RUN/ | |
4503 | COMMON /ARERC1/MULTI1(MAXR) | |
4504 | cc SAVE /ARERC1/ | |
4505 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
4506 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
4507 | & FT1(MAXSTR, MAXR), | |
4508 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
4509 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
4510 | cc SAVE /ARPRC1/ | |
4511 | COMMON /AROUT/ IOUT | |
4512 | cc SAVE /AROUT/ | |
4513 | COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult | |
4514 | cc SAVE /iflow/ | |
4515 | SAVE | |
4516 | DATA IW/0/ | |
4517 | ||
4518 | IW = IW + 1 | |
4519 | DO 1002 J = 1, NUM | |
4520 | DO 1001 I = 1, MULTI1(J) | |
4521 | ITYP = ITYP1(I, J) | |
4522 | IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200 | |
4523 | PX = PX1(I, J) | |
4524 | PY = PY1(I, J) | |
4525 | PZ = PZ1(I, J) | |
4526 | EE = EE1(I, J) | |
4527 | XM = XM1(I, J) | |
4528 | XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) | |
4529 | IF (ABS(PZ) .GE. EE) THEN | |
4530 | PRINT *, 'IN HJANA3' | |
4531 | PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR' | |
4532 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
4533 | PRINT *, ' PZ = ', PZ, ' EE = ', EE | |
4534 | PRINT *, ' XM = ', XM | |
4535 | GOTO 200 | |
4536 | END IF | |
4537 | DXMT = XMT - XM | |
4538 | Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) | |
4539 | c.....rapidity cut for the rapidity distribution | |
4540 | c IY = 1 + int(ABS(Y) / DY) | |
4541 | IY = 1 + int((Y+10.) / DY) | |
4542 | IF (IY .GT. 50) GOTO 100 | |
4543 | dndyh3(IY) = dndyh3(IY) + 1.0 | |
4544 | DEYH3(IY) = DEYH3(IY) + XMT | |
4545 | 100 CONTINUE | |
4546 | c.....insert rapidity cut for mt spectrum here | |
4547 | IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200 | |
4548 | IMT = 1 + int(DXMT / DMT) | |
4549 | IF (IMT .GT. 50) GOTO 200 | |
4550 | DMYH3(IMT) = DMYH3(IMT) + 1.0 / XMT | |
4551 | 200 CONTINUE | |
4552 | 1001 CONTINUE | |
4553 | 1002 CONTINUE | |
4554 | c | |
4555 | RETURN | |
4556 | END | |
4557 | ||
4558 | c----------------------------------------------------------------------- | |
4559 | ||
4560 | c.....analysis subroutine after ARTMN | |
4561 | SUBROUTINE HJANA4 | |
4562 | PARAMETER (MAXSTR=150001, MAXR=1) | |
4563 | c.....y cut for mt spectrum | |
4564 | cbz11/7/99 | |
4565 | c PARAMETER (YMIN = -0.5, YMAX = 0.5) | |
4566 | PARAMETER (YMIN = -1.0, YMAX = 1.0) | |
4567 | cbz11/7/99 end | |
4568 | c.....bin width for mt spectrum and y spectrum | |
4569 | PARAMETER (DMT = 0.05, DY = 0.2) | |
4570 | ||
4571 | DIMENSION dndyh4(50), DMYH4(50), DEYH4(50) | |
4572 | COMMON /RUN/ NUM | |
4573 | cc SAVE /RUN/ | |
4574 | COMMON /ARERC1/MULTI1(MAXR) | |
4575 | cc SAVE /ARERC1/ | |
4576 | COMMON /ARPRC1/ITYP1(MAXSTR, MAXR), | |
4577 | & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), | |
4578 | & FT1(MAXSTR, MAXR), | |
4579 | & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR), | |
4580 | & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR) | |
4581 | cc SAVE /ARPRC1/ | |
4582 | COMMON /AROUT/ IOUT | |
4583 | cc SAVE /AROUT/ | |
4584 | COMMON /fflow/ v2f,etf,xmultf,v2fpi,xmulpi | |
4585 | cc SAVE /fflow/ | |
4586 | SAVE | |
4587 | DATA IW/0/ | |
4588 | ||
4589 | IW = IW + 1 | |
4590 | DO 1002 J = 1, NUM | |
4591 | DO 1001 I = 1, MULTI1(J) | |
4592 | ITYP = ITYP1(I, J) | |
4593 | IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200 | |
4594 | PX = PX1(I, J) | |
4595 | PY = PY1(I, J) | |
4596 | PZ = PZ1(I, J) | |
4597 | EE = EE1(I, J) | |
4598 | XM = XM1(I, J) | |
4599 | XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2) | |
4600 | IF (ABS(PZ) .GE. EE) THEN | |
4601 | PRINT *, 'IN HJANA4' | |
4602 | PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR' | |
4603 | PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY | |
4604 | PRINT *, ' PZ = ', PZ, ' EE = ', EE | |
4605 | PRINT *, ' XM = ', XM | |
4606 | GOTO 200 | |
4607 | END IF | |
4608 | DXMT = XMT - XM | |
4609 | Y = 0.5 * LOG((EE + PZ) / (EE - PZ)) | |
4610 | c.....rapidity cut for the rapidity distribution | |
4611 | c IY = 1 + int(ABS(Y) / DY) | |
4612 | IY = 1 + int((Y+10.) / DY) | |
4613 | IF (IY .GT. 50) GOTO 100 | |
4614 | dndyh4(IY) = dndyh4(IY) + 1.0 | |
4615 | DEYH4(IY) = DEYH4(IY) + XMT | |
4616 | 100 CONTINUE | |
4617 | c.....insert rapidity cut for mt spectrum here | |
4618 | IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200 | |
4619 | IMT = 1 + int(DXMT / DMT) | |
4620 | IF (IMT .GT. 50) GOTO 200 | |
4621 | DMYH4(IMT) = DMYH4(IMT) + 1.0 / XMT | |
4622 | 200 CONTINUE | |
4623 | 1001 CONTINUE | |
4624 | 1002 CONTINUE | |
4625 | c | |
4626 | RETURN | |
4627 | END | |
4628 | ||
4629 | c======================================================================= | |
4630 | ||
4631 | c.....subroutine to get average values for different strings | |
4632 | ||
4633 | SUBROUTINE zpstrg | |
4634 | ||
4635 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
4636 | PARAMETER (MAXPTN=400001) | |
4637 | PARAMETER (MAXSTR=150001) | |
4638 | c REAL*4 YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB | |
4639 | REAL YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB | |
4640 | ||
4641 | COMMON /PARA1/ MUL | |
4642 | cc SAVE /PARA1/ | |
4643 | COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN), | |
4644 | & PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN), | |
4645 | & XMASS5(MAXPTN), ITYP5(MAXPTN) | |
4646 | cc SAVE /prec2/ | |
4647 | COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN) | |
4648 | cc SAVE /ilist8/ | |
4649 | COMMON /SREC1/ NSP, NST, NSI | |
4650 | cc SAVE /SREC1/ | |
4651 | COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR) | |
4652 | cc SAVE /SREC2/ | |
4653 | COMMON/hjcrdn/YP(3,300),YT(3,300) | |
4654 | cc SAVE /hjcrdn/ | |
4655 | COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100), | |
4656 | & K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100), | |
4657 | & PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100) | |
4658 | cc SAVE /HJJET2/ | |
4659 | cbz6/28/99 flow1 | |
4660 | COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50) | |
4661 | cc SAVE /HPARNT/ | |
4662 | cbz6/28/99 flow1 end | |
4663 | common/anim/nevent,isoft,isflag,izpc | |
4664 | cc SAVE /anim/ | |
4665 | common/strg/np(maxstr) | |
4666 | cc SAVE /strg/ | |
4667 | clin-6/06/02 test local freezeout: | |
4668 | common /frzprc/ | |
4669 | & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN), | |
4670 | & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN), | |
4671 | & xmfrz(MAXPTN), | |
4672 | & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast | |
4673 | cc SAVE /frzprc/ | |
4674 | SAVE | |
4675 | ||
4676 | clin-6/06/02 test local freezeout for string melting, | |
4677 | c use space-time values at local freezeout saved in /frzprc/: | |
4678 | if(isoft.eq.5) then | |
4679 | do 1001 I = 1, MUL | |
4680 | ITYP5(i)=idfrz(i) | |
4681 | GX5(i)=gxfrz(i) | |
4682 | GY5(i)=gyfrz(i) | |
4683 | GZ5(i)=gzfrz(i) | |
4684 | FT5(i)=ftfrz(i) | |
4685 | PX5(i)=pxfrz(i) | |
4686 | PY5(i)=pyfrz(i) | |
4687 | PZ5(i)=pzfrz(i) | |
4688 | E5(i)=efrz(i) | |
4689 | XMASS5(i)=xmfrz(i) | |
4690 | 1001 continue | |
4691 | endif | |
4692 | clin-6/06/02-end | |
4693 | ||
4694 | DO 1002 I = 1, MAXSTR | |
4695 | ATAUI(I) = 0d0 | |
4696 | ZT1(I) = 0d0 | |
4697 | ZT2(I) = 0d0 | |
4698 | clin-4/25/03 add zt3(I) to track longitudinal positions of partons/strings: | |
4699 | ZT3(I) = 0d0 | |
4700 | NP(I) = 0 | |
4701 | 1002 CONTINUE | |
4702 | DO 1003 I = 1, MUL | |
4703 | ISTRG = LSTRG1(I) | |
4704 | TAU7 = SQRT(FT5(I) ** 2 - GZ5(I) ** 2) | |
4705 | ATAUI(ISTRG) = ATAUI(ISTRG) + TAU7 | |
4706 | ZT1(ISTRG) = ZT1(ISTRG) + GX5(I) | |
4707 | ZT2(ISTRG) = ZT2(ISTRG) + GY5(I) | |
4708 | ZT3(ISTRG) = ZT3(ISTRG) + GZ5(I) | |
4709 | NP(ISTRG) = NP(ISTRG) + 1 | |
4710 | 1003 CONTINUE | |
4711 | ||
4712 | NSTR = NSP + NST + NSI | |
4713 | ||
4714 | clin-7/03/01 correct averaging on transverse coordinates, no shift needed: | |
4715 | if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then | |
4716 | DO 1004 I = 1, NSTR | |
4717 | IF (NP(I) .NE. 0) THEN | |
4718 | ATAUI(I) = ATAUI(I) / NP(I) | |
4719 | ZT1(I) = ZT1(I) / NP(I) | |
4720 | ZT2(I) = ZT2(I) / NP(I) | |
4721 | ZT3(I) = ZT3(I) / NP(I) | |
4722 | ENDIF | |
4723 | 1004 CONTINUE | |
4724 | return | |
4725 | endif | |
4726 | clin-7/03/01-end | |
4727 | ||
4728 | DO 1005 I = 1, NSTR | |
4729 | IF (NP(I) .NE. 0) THEN | |
4730 | ATAUI(I) = ATAUI(I) / NP(I) | |
4731 | ZT1(I) = ZT1(I) / NP(I) | |
4732 | ZT2(I) = ZT2(I) / NP(I) | |
4733 | ZT3(I) = ZT3(I) / NP(I) | |
4734 | ELSE | |
4735 | IF (I .LE. NSP) THEN | |
4736 | J = I | |
4737 | ZT1(I) = dble(YP(1, J)) | |
4738 | ZT2(I) = dble(YP(2, J)) | |
4739 | ZT3(I) = 0d0 | |
4740 | ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN | |
4741 | J = I - NSP | |
4742 | ZT1(I) = dble(YT(1, J)) | |
4743 | ZT2(I) = dble(YT(2, J)) | |
4744 | ZT3(I) = 0d0 | |
4745 | ELSE | |
4746 | J = I - NSP - NST | |
4747 | ZT1(I) = 0.5d0* | |
4748 | 1 dble((YP(1, IASG(J, 1)) + YT(1, IASG(J, 2)))) | |
4749 | ZT2(I) = 0.5d0 * | |
4750 | 1 dble((YP(2, IASG(J, 1)) + YT(2, IASG(J, 2)))) | |
4751 | ZT3(I) = 0d0 | |
4752 | END IF | |
4753 | END IF | |
4754 | 1005 CONTINUE | |
4755 | ||
4756 | cbz6/28/99 flow1 | |
4757 | BB = HINT1(19) | |
4758 | DO 1006 I = 1, NSTR | |
4759 | IF (NP(I).NE.0) THEN | |
4760 | SHIFT=0d0 | |
4761 | ELSE | |
4762 | SHIFT=0.5d0*dble(BB) | |
4763 | END IF | |
4764 | IF (I .LE. NSP) THEN | |
4765 | ZT1(I) = ZT1(I) + SHIFT | |
4766 | ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN | |
4767 | ZT1(I) = ZT1(I) - SHIFT | |
4768 | END IF | |
4769 | 1006 CONTINUE | |
4770 | cbz6/28/99 flow1 end | |
4771 | c | |
4772 | RETURN | |
4773 | END | |
4774 | ||
4775 | clin-3/2009 | |
4776 | c Initialize hadron weights; | |
4777 | c Can add initial hadrons before the hadron cascade starts (but after ZPC). | |
4778 | subroutine addhad | |
4779 | PARAMETER (MAXSTR=150001,MAXR=1,xmd=1.8756) | |
4780 | double precision smearp,smearh | |
4781 | COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50) | |
4782 | COMMON /ARPRC/ ITYPAR(MAXSTR), | |
4783 | & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR), | |
4784 | & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR), | |
4785 | & XMAR(MAXSTR) | |
4786 | COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR), | |
4787 | 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR), | |
4788 | 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR) | |
4789 | COMMON /smearz/smearp,smearh | |
4790 | COMMON/RNDF77/NSEED | |
4791 | common /para8/ idpert,npertd,idxsec | |
4792 | SAVE | |
4793 | c All hadrons at the start of hadron cascade have the weight of 1 | |
4794 | c except those inserted by the user in this subroutine: | |
4795 | np0=IAINT2(1) | |
4796 | DO i=1,np0 | |
4797 | dpertp(I)=1. | |
4798 | ENDDO | |
4799 | c Specify number, species, weight, initial x,p,m for inserted hadrons here: | |
4800 | nadd=0 | |
4801 | tau0=ARPAR1(1) | |
4802 | DO 100 i=np0+1,np0+nadd | |
4803 | ITYPAR(I)=42 | |
4804 | dpertp(I)=1d0/dble(nadd) | |
4805 | GXAR(I)=5.*(1.-2.*RANART(NSEED)) | |
4806 | GYAR(I)=5.*(1.-2.*RANART(NSEED)) | |
4807 | GZAR(I)=2.*(1.-2.*RANART(NSEED)) | |
4808 | FTAR(I)=0. | |
4809 | PXAR(I)=1. | |
4810 | PYAR(I)=0. | |
4811 | PZAR(I)=1. | |
4812 | XMAR(I)=xmd | |
4813 | c | |
4814 | PEAR(I)=sqrt(PXAR(I)**2+PYAR(I)**2+PZAR(I)**2+XMAR(I)**2) | |
4815 | RAP=0.5*alog((PEAR(I)+PZAR(I))/(PEAR(I)-PZAR(I))) | |
4816 | VX=PXAR(I)/PEAR(I) | |
4817 | VY=PYAR(I)/PEAR(I) | |
4818 | c.....give initial formation time shift and boost according to rapidity: | |
4819 | TAUI=FTAR(I)+TAU0 | |
4820 | FTAR(I)=TAUI*COSH(RAP) | |
4821 | GXAR(I)=GXAR(I)+VX*TAU0*COSH(RAP) | |
4822 | GYAR(I)=GYAR(I)+VY*TAU0*COSH(RAP) | |
4823 | c Allow the intial z-position to be different from the Bjorken picture: | |
4824 | GZAR(I)=TAUI*SINH(RAP)+GZAR(I) | |
4825 | c GZAR(I)=TAUI*SINH(RAP) | |
4826 | zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.) | |
4827 | GZAR(I)=GZAR(I)+zsmear | |
4828 | 100 CONTINUE | |
4829 | IAINT2(1)=IAINT2(1)+nadd | |
4830 | c | |
4831 | if(nadd.ge.1.and.idpert.ne.1.and.idpert.ne.2) then | |
4832 | write(16,*) 'IDPERT must be 1 or 2 to add initial hadrons, | |
4833 | 1 set NPERTD to 0 if you do not need perturbative deuterons' | |
4834 | stop | |
4835 | endif | |
4836 | if(IAINT2(1).gt.MAXSTR) then | |
4837 | write(16,*) 'Too many initial hadrons, array size is exceeded!' | |
4838 | stop | |
4839 | endif | |
4840 | c | |
4841 | return | |
4842 | end |