]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TAmpt/AMPT/amptsub.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Also rename PYTHIA call
[u/mrichter/AliRoot.git] / TAmpt / AMPT / amptsub.f
CommitLineData
0119ef9a 1c....................amptsub.f
2c.....this file contains 4 sections:
3c.....1. ART subroutines;
4c.....2. ART functions;
5c.....3. ART block data;
6c.....4. subprocesses borrowed from other codes.
7c.....5. the previous artana.f
8c.....6. the previous zpcsub.f
9c.....7. subroutine getnp
10c.....Note that Parts1-4 are the previous artsub.f
11c
12c=======================================================================
13c.....subroutine to set up ART parameters and analysis files
14c.....before looping different events
15cms
16cms dlw & gsfs 8/2009 commented out lots of output files
17cms
18 SUBROUTINE ARTSET
19c
20 PARAMETER (AMU= 0.9383)
21 double precision dpcoal,drcoal,ecritl
22 INTEGER ZTA, ZPR
23 common /gg/ dx,dy,dz,dpx,dpy,dpz
24clin-10/03/03
25c "SAVE " (without argument) is used for most subroutines and functions,
26c this is important for the success when using "f77" to compile:
27cc SAVE /gg/
28 common /zz/ zta,zpr
29cc SAVE /zz/
30 COMMON /RUN/ NUM
31cc SAVE /RUN/
32 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
33cc SAVE /input1/
34 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
35 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
36cc SAVE /INPUT2/
37 COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
38cc SAVE /INPUT3/
39 common /imulst/ iperts
40cc SAVE /imulst/
41 common /coal/dpcoal,drcoal,ecritl
42 common/anim/nevent,isoft,isflag,izpc
43 common /para7/ ioscar,nsmbbbar,nsmmeson
44 SAVE
45clin-10/03/03 ecritl: local energy density below which a parton
46c will freeze out (in GeV/fm^3), for improvements on string melting,
47c not used in this version of AMPT:
48clin-4/2008
49c data ecritl/1.d0/
50 ecritl=1.d0
51c
52c combine ART initialization into ampt.ini:
53c (Note that the following values are relics from the old ART structure)
54c.....input parameter file
55c OPEN(13, FILE = 'art1.ini', STATUS = 'UNKNOWN')
56c READ (13, *) MASSTA, ZTA
57 MASSTA=1
58 ZTA=1
59c write(12,*) massta, zta, ' massta, zta'
60c READ (13, *) MASSPR, ZPR
61 MASSPR=1
62 ZPR=1
63c write(12,*) masspr, zpr, ' masspr, zpr'
64c READ (13, *) PLAB, IPLAB
65 PLAB=14.6
66 IPLAB=2
67c 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.
74c READ (13, *) ZEROPT
75 ZEROPT=0.
76c write(12,*) zeropt, ' zeropt'
77clin-10/03/03 ISEED was used as a seed for random number inside ART,
78c not used in AMPT:
79 ISEED=700721
80c 0/1: (Normal or Perturbative) multistrange partice production.
81c Perturbative option is disabled for now:
82 iperts=0
83c READ (13, *) MANYB, B0, BI, BM
84c 2/04/00 MANYB MUST BE SET TO 1 !
85c 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
90c write(12,*) manyb, b0, bi, bm, ' manyb, b0, bi, bm'
91c READ (13, *) ISEED
92c write(12,*) iseed, ' iseed'
93c READ (13, *) DT
94c write(12,*) dt, ' dt'
95c READ (13, *) NTMAX
96c write(12,*) ntmax, ' ntmax'
97c READ (13, *) ICOLL
98 ICOLL=-1
99c write(12,*) icoll, ' icoll'
100c READ (13, *) NUM
101c 2/11/03 run events without test particles for now:
102 NUM=1
103c write(12,*) num, ' num'
104c READ (13, *) INSYS
105 INSYS=1
106c write(12,*) insys, ' insys'
107c READ (13, *) IPOT
108 IPOT=3
109c write(12,*) ipot, ' ipot'
110c READ (13, *) MODE
111 MODE=0
112 IF(ICOLL.EQ.-1)IPOT=0
113c write(12,*) mode, ' mode'
114c READ (13, *) DX, DY, DZ
115 DX=2.73
116 DY=2.73
117 DZ=2.73
118c write(12,*) dx,dy,dz,' dx,dy,dz'
119c READ (13, *) DPX, DPY, DPZ
120 DPX=0.6
121 DPY=0.6
122 DPZ=0.6
123c write(12,*) dpx,dpy,dpz,' dpx,dpy,dpz'
124c READ (13, *) IAVOID
125 IAVOID=1
126c write(12,*) iavoid, ' iavoid'
127c READ (13, *) IMOMEN
128 IMOMEN=1
129c write(12,*) imomen, ' imomen'
130 if(icoll.eq.-1)imomen=3
131c READ (13, *) NFREQ
132 NFREQ=10
133c write(12,*) nfreq, ' nfreq'
134c READ (13, *) ICFLOW
135 ICFLOW=0
136c write(12,*) ICFLOW, ' ICFLOW'
137c READ (13, *) ICRHO
138 ICRHO=0
139c write(12,*) ICRHO, ' ICRHO'
140c READ (13, *) ICOU
141 ICOU=0
142c write(12,*)icou, ' icou'
143* kaon potential control parameter
144* KMUL IS A MULTIPLIER TO THE STANDARD K-N SCATTERING LENGTH
145c READ (13, *) KPOTEN, KMUL
146 KPOTEN=0
147 KMUL=1
148c 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.
152c READ (13, *) DENCUT
153 DENCUT=15
154c write(12,*)dencut, ' dencut'
155* test reactions in a box of side-length cycbox
156* input cycbox
157c READ (13, *) CYCBOX
158 CYCBOX=0
159c write(12,*) cycbox, ' cycbox'
160c
161clin-5b/2008
162c if(ioscar.eq.2) then
163 if(ioscar.eq.2.or.ioscar.eq.3) then
164cms OPEN (92,FILE='ana/parton-initial-afterPropagation.dat',
165cms 1 STATUS = 'UNKNOWN')
166 endif
167 if(ioscar.eq.3) then
168clin-6/2009 write out full parton collision history:
169cms OPEN (95,FILE='ana/parton-collisionsHistory.dat',
170cms 1 STATUS='UNKNOWN')
171clin-6/2009 write out initial minijet information:
172cms OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat',
173cms 1 STATUS='UNKNOWN')
174clin-6/2009 write out parton info after coalescence:
175 if(isoft.eq.4.or.isoft.eq.5) then
176cms OPEN (85,FILE='ana/parton-after-coalescence.dat',
177cms 1 STATUS='UNKNOWN')
178 endif
179 endif
180clin-6/2009 write out initial transverse positions of initial nucleons:
181cms OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN')
182
183 RETURN
184 END
185
186c-----------------------------------------------------------------------
187
188c.....subroutine to initialize cascade.
189
190 SUBROUTINE ARINI
191
192c.....before invoking ARINI:
193c.....IAPAR2(1), IAINT2(1) must be set.
194 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
195cc SAVE /ARPRNT/
196 SAVE
197
198ctest off for resonance (phi, K*) studies:
199c OPEN (89, FILE = 'ana/decay_rec.dat', STATUS = 'UNKNOWN')
200
201 IFLG = IAPAR2(1)
202 GOTO (200, 200, 300) IFLG
203
204c.....error choice of initialization
205 PRINT *, 'IAPAR2(1) must be 1, 2, or 3'
206 STOP
207
208c.....to use default initial conditions generated by the cascade,
209c.....or to read in initial conditions.
210 200 RETURN
211
212c.....to generate formation time and the position at formation time from
213c.....read-in initial conditions with an averaged formation proper time.
214 300 CALL ARINI1
215c.....ordering the particle label according to increasing order of
216c.....formation time.
217 CALL ARTORD
218 RETURN
219
220 END
221
222c-----------------------------------------------------------------------
223
224c.....subroutine to generate formation time and position at formation time
225c.....from read-in initial conditions with an averaged formation proper
226c.....time.
227
228 SUBROUTINE ARINI1
229
230c.....before invoking ARINI1:
231c.....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)
236cc 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)
241cc SAVE /ARPRC/
242 COMMON /smearz/smearp,smearh
243cc SAVE /smearz/
244 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
245cc SAVE /input1/
246 common/anim/nevent,isoft,isflag,izpc
247cc SAVE /anim/
248 common /nzpc/nattzp
249cc SAVE /nzpc/
250 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
251cc SAVE /HPARNT/
252 COMMON/RNDF77/NSEED
253cc SAVE /RNDF77/
254 common /para8/ idpert,npertd,idxsec
255 SAVE
256clin-5/2008 for perturbatively-produced hadrons (currently only deuterons):
257cms OPEN (91, FILE = 'ana/deuteron_processes.dat',
258cms 1 STATUS = 'UNKNOWN')
259 if(idpert.eq.1.or.idpert.eq.2) then
260cms OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN')
261 endif
262c.....generate formation time and position at formation time.
263 TAU0 = ARPAR1(1)
264 NP = IAINT2(1)
265clin-7/10/01 initial positions already given for hadrons
266c 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)
288clin-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
297clin-7/10/01-end
298clin-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
310c 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)
316c.....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)
321c 4/25/03: hadron z-position upon formation determined the same way as x,y:
322 GZAR(I) = TAUI * SINH(RAP)
323c the old prescription:
324c GZAR(I) = GZAR(I) + TAU0 * SINH(RAP)
325 zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
326 GZAR(I)=GZAR(I)+zsmear
327cbz1/28/99end
328c 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
332clin-5/2008:
333c 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
339clin-3/2009 cleanup of program flow:
340 endif
341
342clin-3/2009 Add initial hadrons before the hadron cascade starts:
343 call addhad
344
345 RETURN
346 END
347
348c-----------------------------------------------------------------------
349
350c.....subroutine to order particle labels according to increasing
351c.....formation time
352
353 SUBROUTINE ARTORD
354
355c.....before invoking ARTORD:
356c.....IAINT2(1) must be set:
357 PARAMETER (MAXSTR=150001,MAXR=1)
358 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
359cc 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)
364cc SAVE /ARPRC/
365clin-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)
370c
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
378c
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)
392clin-3/2009:
393 dptemp(I) = dpertp(I)
394 1001 CONTINUE
395 CALL ARINDX(MAXSTR, NP, FT0, INDX)
396 DO 1002 I = 1, NP
397cbz12/3/98
398c IF (ITYP0(INDX(I)) .EQ. 211) THEN
399c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 321) THEN
400c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
401c & ITYP0(INDX(I)) .EQ. 2112 .OR. ITYP0(INDX(I)) .EQ. -211 .OR.
402c & ITYP0(INDX(I)) .EQ. 111) THEN
403c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
404c & ITYP0(INDX(I)) .EQ. 2112) THEN
405 NPAR = NPAR + 1
406c ITYPAR(I) = ITYP0(INDX(I))
407c GXAR(I) = GX0(INDX(I))
408c GYAR(I) = GY0(INDX(I))
409c GZAR(I) = GZ0(INDX(I))
410c FTAR(I) = FT0(INDX(I))
411c PXAR(I) = PX0(INDX(I))
412c PYAR(I) = PY0(INDX(I))
413c PZAR(I) = PZ0(INDX(I))
414c PEAR(I) = EE0(INDX(I))
415c 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))
426clin-3/2009:
427 dpertp(NPAR)=dptemp(INDX(I))
428c END IF
429cbz12/3/98end
430 1002 CONTINUE
431 IAINT2(1) = NPAR
432c
433 RETURN
434 END
435
436c-----------------------------------------------------------------------
437
438c.....subroutine to copy individually generated particle record into
439c.....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)
445cc 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)
450cc SAVE /ARPRC/
451 COMMON /ARERC1/MULTI1(MAXR)
452cc 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)
458cc SAVE /ARPRC1/
459 COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
460cc SAVE /tdecay/
461 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
462cc SAVE /input1/
463 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
464 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
465cc 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)
470cc 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)
485clin-3/2009 hadron weights are initialized in addhad():
486clin-5/2008 all hadrons not perturbatively-produced have the weight of 1:
487c dpp1(I,K)=1.
488 dpp1(I,K)=dpertp(I)
489 1001 CONTINUE
490
491c initialize final time of each particle to ntmax*dt except for
492c 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
497c
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
507c=======================================================================
508
509c.....function to convert PDG flavor code into ART flavor code.
510
511 FUNCTION IARFLV(IPDG)
512
513 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
514cc SAVE /input1/
515 COMMON/RNDF77/NSEED
516cc SAVE /RNDF77/
517 SAVE
518
519c.....anti-Delta-
520 IF (IPDG .EQ. -1114) THEN
521 IARFLV = -6
522 RETURN
523 END IF
524
525c.....anti-Delta0
526 IF (IPDG .EQ. -2114) THEN
527 IARFLV = -7
528 RETURN
529 END IF
530
531c.....anti-Delta+
532 IF (IPDG .EQ. -2214) THEN
533 IARFLV = -8
534 RETURN
535 END IF
536
537c.....anti-Delta++
538 IF (IPDG .EQ. -2224) THEN
539 IARFLV = -9
540 RETURN
541 END IF
542
543cbzdbg2/23/99
544c.....anti-proton
545 IF (IPDG .EQ. -2212) THEN
546 IARFLV = -1
547 RETURN
548 END IF
549
550c.....anti-neutron
551 IF (IPDG .EQ. -2112) THEN
552 IARFLV = -2
553 RETURN
554 END IF
555cbzdbg2/23/99end
556
557c.....eta
558 IF (IPDG .EQ. 221) THEN
559 IARFLV = 0
560 RETURN
561 END IF
562
563c.....proton
564 IF (IPDG .EQ. 2212) THEN
565 IARFLV = 1
566 RETURN
567 END IF
568
569c.....neutron
570 IF (IPDG .EQ. 2112) THEN
571 IARFLV = 2
572 RETURN
573 END IF
574
575c.....pi-
576 IF (IPDG .EQ. -211) THEN
577 IARFLV = 3
578 RETURN
579 END IF
580
581c.....pi0
582 IF (IPDG .EQ. 111) THEN
583 IARFLV = 4
584 RETURN
585 END IF
586
587c.....pi+
588 IF (IPDG .EQ. 211) THEN
589 IARFLV = 5
590 RETURN
591 END IF
592
593c.....Delta-
594 IF (IPDG .EQ. 1114) THEN
595 IARFLV = 6
596 RETURN
597 END IF
598
599c.....Delta0
600 IF (IPDG .EQ. 2114) THEN
601 IARFLV = 7
602 RETURN
603 END IF
604
605c.....Delta+
606 IF (IPDG .EQ. 2214) THEN
607 IARFLV = 8
608 RETURN
609 END IF
610
611c.....Delta++
612 IF (IPDG .EQ. 2224) THEN
613 IARFLV = 9
614 RETURN
615 END IF
616
617c.....Lambda
618 IF (IPDG .EQ. 3122) THEN
619 IARFLV = 14
620 RETURN
621 END IF
622
623c.....Lambda-bar
624 IF (IPDG .EQ. -3122) THEN
625 IARFLV = -14
626 RETURN
627 END IF
628
629c.....Sigma-
630 IF (IPDG .EQ. 3112) THEN
631 IARFLV = 15
632 RETURN
633 END IF
634
635c.....Sigma-bar
636 IF (IPDG .EQ. -3112) THEN
637 IARFLV = -15
638 RETURN
639 END IF
640
641c.....Sigma0
642 IF (IPDG .EQ. 3212) THEN
643 IARFLV = 16
644 RETURN
645 END IF
646
647c.....Sigma0-bar
648 IF (IPDG .EQ. -3212) THEN
649 IARFLV = -16
650 RETURN
651 END IF
652
653c.....Sigma+
654 IF (IPDG .EQ. 3222) THEN
655 IARFLV = 17
656 RETURN
657 END IF
658
659c.....Sigma+ -bar
660 IF (IPDG .EQ. -3222) THEN
661 IARFLV = -17
662 RETURN
663 END IF
664
665c.....K-
666 IF (IPDG .EQ. -321) THEN
667 IARFLV = 21
668 RETURN
669 END IF
670
671c.....K+
672 IF (IPDG .EQ. 321) THEN
673 IARFLV = 23
674 RETURN
675 END IF
676
677c.....temporary entry for K0
678 IF (IPDG .EQ. 311) THEN
679 IARFLV = 23
680 RETURN
681 END IF
682
683c.....temporary entry for K0bar
684 IF (IPDG .EQ. -311) THEN
685 IARFLV = 21
686 RETURN
687 END IF
688
689c.....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
700c.....rho-
701 IF (IPDG .EQ. -213) THEN
702 IARFLV = 25
703 RETURN
704 END IF
705
706c.....rho0
707 IF (IPDG .EQ. 113) THEN
708 IARFLV = 26
709 RETURN
710 END IF
711
712c.....rho+
713 IF (IPDG .EQ. 213) THEN
714 IARFLV = 27
715 RETURN
716 END IF
717
718c.....omega
719 IF (IPDG .EQ. 223) THEN
720 IARFLV = 28
721 RETURN
722 END IF
723
724c.....phi
725 IF (IPDG .EQ. 333) THEN
726 IARFLV = 29
727 RETURN
728 END IF
729
730c.....K*+
731 IF (IPDG .EQ. 323) THEN
732 IARFLV = 30
733 RETURN
734 END IF
735c.....K*-
736 IF (IPDG .EQ. -323) THEN
737 IARFLV = -30
738 RETURN
739 END IF
740c.....temporary entry for K*0
741 IF (IPDG .EQ. 313) THEN
742 IARFLV = 30
743 RETURN
744 END IF
745c.....temporary entry for K*0bar
746 IF (IPDG .EQ. -313) THEN
747 IARFLV = -30
748 RETURN
749 END IF
750
751c...... eta-prime
752 IF (IPDG .EQ. 331) THEN
753 IARFLV = 31
754 RETURN
755 END IF
756
757c...... a1
758c IF (IPDG .EQ. 777) THEN
759c IARFLV = 32
760c RETURN
761c END IF
762
763c... cascade-
764 IF (IPDG .EQ. 3312) THEN
765 IARFLV = 40
766 RETURN
767 END IF
768
769c... cascade+ (bar)
770 IF (IPDG .EQ. -3312) THEN
771 IARFLV = -40
772 RETURN
773 END IF
774
775c... cascade0
776 IF (IPDG .EQ. 3322) THEN
777 IARFLV = 41
778 RETURN
779 END IF
780
781c... cascade0 -bar
782 IF (IPDG .EQ. -3322) THEN
783 IARFLV = -41
784 RETURN
785 END IF
786
787c... Omega-
788 IF (IPDG .EQ. 3334) THEN
789 IARFLV = 45
790 RETURN
791 END IF
792
793c... Omega+ (bar)
794 IF (IPDG .EQ. -3334) THEN
795 IARFLV = -45
796 RETURN
797 END IF
798
799c... Di-Omega
800 IF (IPDG .EQ. 6666) THEN
801 IARFLV = 44
802 RETURN
803 END IF
804c sp06/05/01 end
805
806clin-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
812c.....other
813 IARFLV = IPDG + 10000
814
815 RETURN
816 END
817
818c-----------------------------------------------------------------------
819
820c.....function to convert ART flavor code into PDG flavor code.
821
822 FUNCTION INVFLV(IART)
823
824 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
825cc SAVE /input1/
826 COMMON/RNDF77/NSEED
827cc SAVE /RNDF77/
828 SAVE
829
830c.....anti-Delta-
831 IF (IART .EQ. -6) THEN
832 INVFLV = -1114
833 RETURN
834 END IF
835
836c.....anti-Delta0
837 IF (IART .EQ. -7) THEN
838 INVFLV = -2114
839 RETURN
840 END IF
841
842c.....anti-Delta+
843 IF (IART .EQ. -8) THEN
844 INVFLV = -2214
845 RETURN
846 END IF
847
848c.....anti-Delta++
849 IF (IART .EQ. -9) THEN
850 INVFLV = -2224
851 RETURN
852 END IF
853
854cbzdbg2/23/99
855c.....anti-proton
856 IF (IART .EQ. -1) THEN
857 INVFLV = -2212
858 RETURN
859 END IF
860
861c.....anti-neutron
862 IF (IART .EQ. -2) THEN
863 INVFLV = -2112
864 RETURN
865 END IF
866cbzdbg2/23/99end
867
868c.....eta
869 IF (IART .EQ. 0) THEN
870 INVFLV = 221
871 RETURN
872 END IF
873
874c.....proton
875 IF (IART .EQ. 1) THEN
876 INVFLV = 2212
877 RETURN
878 END IF
879
880c.....neutron
881 IF (IART .EQ. 2) THEN
882 INVFLV = 2112
883 RETURN
884 END IF
885
886c.....pi-
887 IF (IART .EQ. 3) THEN
888 INVFLV = -211
889 RETURN
890 END IF
891
892c.....pi0
893 IF (IART .EQ. 4) THEN
894 INVFLV = 111
895 RETURN
896 END IF
897
898c.....pi+
899 IF (IART .EQ. 5) THEN
900 INVFLV = 211
901 RETURN
902 END IF
903
904c.....Delta-
905 IF (IART .EQ. 6) THEN
906 INVFLV = 1114
907 RETURN
908 END IF
909
910c.....Delta0
911 IF (IART .EQ. 7) THEN
912 INVFLV = 2114
913 RETURN
914 END IF
915
916c.....Delta+
917 IF (IART .EQ. 8) THEN
918 INVFLV = 2214
919 RETURN
920 END IF
921
922c.....Delta++
923 IF (IART .EQ. 9) THEN
924 INVFLV = 2224
925 RETURN
926 END IF
927
928cc.....N*(1440), N*(1535) temporary entry
929c IF (IART .GE. 10 .AND. IART .LE.13) THEN
930c INVFLV = 0
931c RETURN
932c END IF
933
934c.....Lambda
935 IF (IART .EQ. 14) THEN
936 INVFLV = 3122
937 RETURN
938 END IF
939c.....Lambda-bar
940 IF (IART .EQ. -14) THEN
941 INVFLV = -3122
942 RETURN
943 END IF
944
945cbz3/12/99
946c.....temporary entry for Sigma's
947c IF (IART .EQ. 15) THEN
948c R = RANART(NSEED)
949c IF (R .GT. 2. / 3.) THEN
950c INVFLV = 3112
951c ELSE IF (R .GT. 1./ 3. .AND. R .LE. 2. / 3.) THEN
952c INVFLV = 3212
953c ELSE
954c INVFLV = 3222
955c END IF
956c RETURN
957c END IF
958
959c.....Sigma-
960 IF (IART .EQ. 15) THEN
961 INVFLV = 3112
962 RETURN
963 END IF
964
965c.....Sigma- bar
966 IF (IART .EQ. -15) THEN
967 INVFLV = -3112
968 RETURN
969 END IF
970
971c.....Sigma0
972 IF (IART .EQ. 16) THEN
973 INVFLV = 3212
974 RETURN
975 END IF
976
977c.....Sigma0 -bar
978 IF (IART .EQ. -16) THEN
979 INVFLV = -3212
980 RETURN
981 END IF
982
983c.....Sigma+
984 IF (IART .EQ. 17) THEN
985 INVFLV = 3222
986 RETURN
987 END IF
988
989c.....Sigma+ -bar
990 IF (IART .EQ. -17) THEN
991 INVFLV = -3222
992 RETURN
993 END IF
994
995clin-2/23/03 K0S and K0L are generated at the last timestep:
996c.....temporary entry for K- and K0bar
997 IF (IART .EQ. 21) THEN
998c R = RANART(NSEED)
999c IF (R .GT. 0.5) THEN
1000 INVFLV = -321
1001c ELSE
1002c INVFLV = -311
1003c R = RANART(NSEED)
1004c IF (R .GT. 0.5) THEN
1005c INVFLV = 310
1006c ELSE
1007c INVFLV = 130
1008c END IF
1009c END IF
1010 RETURN
1011 END IF
1012
1013c.....temporary entry for K+ and K0
1014 IF (IART .EQ. 23) THEN
1015c R = RANART(NSEED)
1016c IF (R .GT. 0.5) THEN
1017 INVFLV = 321
1018c ELSE
1019c INVFLV = 311
1020c R = RANART(NSEED)
1021c IF (R .GT. 0.5) THEN
1022c INVFLV = 310
1023c ELSE
1024c INVFLV = 130
1025c END IF
1026c END IF
1027 RETURN
1028 END IF
1029
1030c.....K0Long:
1031 IF (IART .EQ. 22) THEN
1032 INVFLV = 130
1033 RETURN
1034 ENDIF
1035c.....K0Short:
1036 IF (IART .EQ. 24) THEN
1037 INVFLV = 310
1038 RETURN
1039 ENDIF
1040
1041c.....rho-
1042 IF (IART .EQ. 25) THEN
1043 INVFLV = -213
1044 RETURN
1045 END IF
1046
1047c.....rho0
1048 IF (IART .EQ. 26) THEN
1049 INVFLV = 113
1050 RETURN
1051 END IF
1052
1053c.....rho+
1054 IF (IART .EQ. 27) THEN
1055 INVFLV = 213
1056 RETURN
1057 END IF
1058
1059c.....omega
1060 IF (IART .EQ. 28) THEN
1061 INVFLV = 223
1062 RETURN
1063 END IF
1064
1065c.....phi
1066 IF (IART .EQ. 29) THEN
1067 INVFLV = 333
1068 RETURN
1069 END IF
1070
1071c.....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
1078c.....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
1085c... eta-prime (bar)
1086 IF (IART .EQ. 31) THEN
1087 INVFLV = 331
1088 RETURN
1089 END IF
1090
1091c... a1
1092 IF (IART .EQ. 32) THEN
1093 INVFLV = 777
1094 RETURN
1095 END IF
1096
1097c... cascade-
1098 IF (IART .EQ. 40) THEN
1099 INVFLV = 3312
1100 RETURN
1101 END IF
1102
1103c... cascade+ (bar)
1104 IF (IART .EQ. -40) THEN
1105 INVFLV = -3312
1106 RETURN
1107 END IF
1108
1109c... cascade0
1110 IF (IART .EQ. 41) THEN
1111 INVFLV = 3322
1112 RETURN
1113 END IF
1114
1115c... cascade0 -bar
1116 IF (IART .EQ. -41) THEN
1117 INVFLV = -3322
1118 RETURN
1119 END IF
1120
1121c... Omega-
1122 IF (IART .EQ. 45) THEN
1123 INVFLV = 3334
1124 RETURN
1125 END IF
1126
1127c... Omega+ (bar)
1128 IF (IART .EQ. -45) THEN
1129 INVFLV = -3334
1130 RETURN
1131 END IF
1132
1133c... Di-Omega
1134 IF (IART .EQ. 44) THEN
1135 INVFLV = 6666
1136 RETURN
1137 END IF
1138c sp 12/19/00 end
1139
1140clin-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
1148c
1149c.....other
1150 INVFLV = IART - 10000
1151
1152 RETURN
1153 END
1154
1155c=======================================================================
1156
1157 BLOCK DATA ARDATA
1158
1159 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
1160cc 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
1169c=======================================================================
1170
1171c.....Routine borrowed from ZPC.
1172c.....double precision is modified to real*4.
1173
1174cbz1/29/99
1175c subroutine index1(n, m, arrin, indx)
1176 subroutine arindx(n, m, arrin, indx)
1177cbz1/29/99end
1178c indexes the first m elements of ARRIN of length n, i.e., outputs INDX
1179c such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m
1180
1181c 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
1225c-----------------------------------------------------------------------
1226
1227c.....extracted from G. Song's ART expasion including K- interactions
1228c.....file `NEWKAON.FOR'
1229
1230c 5/01/03 send iblock value into art1f.f, necessary for resonance studies:
1231c subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1232c & 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)
1238cc SAVE /AA/
1239 COMMON /BB/ P(3,MAXSTR)
1240cc SAVE /BB/
1241 COMMON /CC/ E(MAXSTR)
1242cc SAVE /CC/
1243 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1244cc SAVE /EE/
1245 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
1246cc SAVE /BG/
1247 COMMON /NN/NNN
1248cc SAVE /NN/
1249 COMMON /RUN/NUM
1250cc SAVE /RUN/
1251 COMMON /PA/RPION(3,MAXSTR,MAXR)
1252cc SAVE /PA/
1253 COMMON /PB/PPION(3,MAXSTR,MAXR)
1254cc SAVE /PB/
1255 COMMON /PC/EPION(MAXSTR,MAXR)
1256cc SAVE /PC/
1257 COMMON /PD/LPION(MAXSTR,MAXR)
1258cc SAVE /PD/
1259 COMMON/RNDF77/NSEED
1260cc SAVE /RNDF77/
1261 SAVE
1262c
1263 logical lb1bn, lb2bn,lb1mn,lb2mn
1264cbz3/7/99 neutralk
1265c logical lb1bn1, lb2bayon1, lb1bn0, lb2bn0
1266 logical lb1bn1, lb2bn1, lb1bn0, lb2bn0
1267cbz3/7/99 neutralk end
1268 logical lb1mn0, lb2mn0, lb1mn1, lb2mn1
1269 logical lb1mn2, lb2mn2
1270 icase=-1
1271c icase: flag for the type of reaction that is going to happen.
1272c icase=-1, no desired reaction, return to main program.
1273c 1, NN,ND,DD
1274c 2, PI+N, PI+D
1275c 3, K(-) absorption.
1276 nchrg=-100
1277c 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
1300c 1. consider N+N, N+Resonance, R + R reactions
1301 if(lb1bn.and.lb2bn) then
1302c NN,ND,DD:
1303 icase=1
1304c 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
1334c brsig = x2kaon_no_isospin(srt)
1335 if(nchrg.ge.-1.and.nchrg.le.2) then
1336c K,Kbar prduction x sect.
1337 brsig = x2kaon(srt)
1338 else
1339 brsig=0.0
1340c if(nchrg.eq.-2.or.nchrg.eq.3) then
1341c brsig = x2kaon(srt+0.938-1.232)
1342c else
1343c nchrg=4
1344c brsig = x2kaon(srt+2.*(0.938-1.232))
1345c endif
1346 endif
1347
1348cbz3/7/99 neutralk
1349 BRSIG = 2.0 * BRSIG
1350cbz3/7/99 neutralk end
1351
1352 endif
1353
1354c 2. consider PI(meson:eta,omega,rho,phi) + N(N*,D)
1355 if((lb1bn.and.lb2mn).OR.(lb2bn.and.lb1mn)) then
1356c 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
1366cbz3/2/99/song
1367c if(lb1bn1.or.lb2bn1) brsig=2.0*sigma0
1368c 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
1371cbz3/2/99/song end
1372c 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
1380cbz3/2/99/song
1381c brsig=1.5*sigma0
1382 brsig=3.0*sigma0
1383cbz3/2/99/song end
1384cbz3/11/99/song
1385c ratiok = 1./3.
1386 ratiok = 2./3.
1387cbz3/11/99/song end
1388
1389c ratiok: the ratio of channels: ->nK+k- vs. -> pK0K-
1390 endif
1391 if(lb1bn0.or.lb2bn0) then
1392 brsig=2.5*sigma0
1393cbz3/2/99/song
1394c ratiok = 0.8
1395 ratiok = 0.2
1396cbz3/2/99/song end
1397 endif
1398c if(lb1.eq.6.or.lb2.eq.6) then
1399c lb=6 : D-
1400c brsig=1.5*sigma0
1401c ratiok = 0.5
1402c 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
1409c if(lb1.eq.6.or.lb2.eq.6) brsig=sigma0
1410 endif
1411c if((lb1.eq.6.and.lb2mn2).or.(lb2.eq.6.and.lb1mn2))then
1412c nchrg=-2
1413c endif
1414c if((lb1bn1.and.lb2mn1).or.(lb2bn1.and.lb1mn1)
1415c & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1416c nchrg=2
1417c endif
1418
1419cbz3/11/99 neutralk
1420 if((lb1.eq.6.and.lb2mn2)
1421 & .or.(lb2.eq.6.and.lb1mn2))then
1422 nchrg=-2
1423 endif
1424cbz3/11/99 neutralk
1425cbz3/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
1431cbz3/8/99 neutralk end
1432
1433cbz3/7/99 neutralk
1434 IF (NCHRG .GE. -2 .AND. NCHRG .LE. 2) THEN
1435 BRSIG = 3.0 * SIGMA0
1436 END IF
1437cbz3/7/99 neutralk end
1438
1439 endif
1440
1441c 3. consider K- + N(N*,D) absorption.
1442c 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
1445c bmass=em1+em2-aka
1446 bmass=0.938
1447 if(srt.le.(bmass+aka)) then
1448cbz3/2/99
1449c write(100,*)'--lb1,lb2,em1,em2,srt',lb1,lb2,em1,em2,srt
1450cbz3/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
1458c 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
1466c 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
1473c 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
1480c K- + D++
1481 nchrg=1
1482 sigela=akPel(pkaon)
1483 sigsgm=2.*akPsgm(pkaon)
1484 sig=sigela+sigsgm+akPlam(pkaon)
1485 endif
1486
1487cbz3/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)
1491cbz3/8/99 neutralk end
1492
1493 if(sig.gt.1.e-7) then
1494c K(-) + N reactions
1495 icase=3
1496 brel=sigela/sig
1497 brsgm=sigsgm/sig
1498c branch_lambda=akNlam(pkaon)/sig
1499 brsig = sig
1500 endif
1501 endif
1502
1503c 4. meson + hyperon -> K- + N
1504c if(((lb1.ge.14.and.lb1.le.17).and.lb2mn).OR.
1505c & ((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
1508c 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
1513c 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
1520c 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
1530c 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
1537c D++
1538 bmass=1.232
1539 endif
1540 sig = 0.
1541 if(nchrg.ne.-100.and.srt.gt.(aka+bmass)) then
1542c PI+sigma or PI + Lambda => Kbar + N reactions
1543 icase=4
1544c 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)
1546c 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)
1550c sigma + pi
1551 else
1552c K-p or K-D++
1553 if(nchrg.ge.0) sigma0=akPsgm(pkaon)
1554c K-n or K-D-
1555 if(nchrg.lt.0) sigma0=akNsgm(pkaon)
1556
1557cbz3/8/99 neutralk
1558 SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1559cbz3/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
1564cbz3/8/99 neutralk
1565c if(nchrg.eq.-2.or.nchrg.eq.1) sig=2.*sig K-D++, K-D-
1566c K0barD++, K-D-
1567 if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig
1568
1569cbz3/8/99 neutralk end
1570
1571c the factor 2 comes from spin of delta, which is 3/2
1572c detailed balance. copy from Page 423 of N.P. A614 1997
1573
1574cbz3/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
1582cbz3/8/99 neutralk end
1583 brsig = sig
1584 if(sig.lt.1.e-7) sig = 1.e-7
1585 endif
1586csp05/07/01
1587* comment icase=4 statement below if only inelastic
1588c PI+L/Si => Kbar + N OR ELASTIC SCATTERING
1589 icase=4
1590 brsig = sig
1591c elastic xsecn of 10mb
1592 sigela = 10.
1593 sig = sig + sigela
1594 brel = sigela/sig
1595cc brsig = sig
1596csp05/07/01 end
1597 endif
1598c
1599c if(em2.lt.0.2.and.em1.lt.0.2) then
1600c PI + PI
1601c icase=5
1602c assumed PI PI total x section.
1603c sig=50.
1604c Mk + Mkbar
1605c s0=aka+aka
1606c brsig = 0.
1607c if(srt.gt.s0) brsig = 2.7*(1.-s0**2/srt**2)**0.76
1608c x section for PIPI->KKbar PRC43 (1991) 1881
1609c 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
1620c ec=3.59709
1621c 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
1625c no anti-kaon production
1626 ictrl = -1
1627c in=in+1
1628c write(60,*)'--------------distance-----',in
1629 return
1630 endif
1631
1632clin-10/24/02 set to 0: ik,ik0-3,il,im,im3-4,in,inpion,ipipi,
1633c 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
1656c 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
1661c 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
1666c ratio=sgsum/ik0/40.
1667 endif
1668 endif
1669 if(icase.eq.2) inpion=inpion+1
1670 if(icase.eq.5) ipipi=ipipi+1
1671c write(62,*)'ik1,ik2,ik3',ik1,ik2,ik3,ratio_1,ratio_3,ratio
1672c write(62,*)'inpion,ipipi',inpion,ipipi
1673 if(RANART(NSEED).gt.(brsig/sig)) then
1674c no anti-kaon production
1675 ictrl = -1
1676 return
1677 endif
1678 il=il+1
1679c kaons could be created now.
1680 if(icase.eq.1) then
1681 in=in+1
1682c 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
1688c call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1689c & 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
1693c
1694 if(icase.eq.3) then
1695 im3=im3+1
1696c write(63,*)'im3,lb1,lb2,pkaon',im3,lb1,lb2,pkaon
1697c write(63,*)'sig,el,sigma',sig,brel,brsgm
1698c 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)
1701c this subroutine format is diff. since three final states are possible
1702 endif
1703c
1704
1705 if(icase.eq.4) then
1706 im4=im4+1
1707c write(64,*)'im4,sigma0,branch,sig=',im4,sigma0,brsig,sig
1708c write(64,*)'lb1,lb2,em1,em2,pkaon=',lb1,lb2,em1,em2,pkaon
1709
1710csp06/07/01
1711 if(RANART(NSEED).lt.brel) then
1712 ielstc = 1
1713 else
1714 ielstc = 0
1715 endif
1716c call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1717c & pcx,pcy,pcz,nchrg)
1718 call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1719 & pcx,pcy,pcz,nchrg,iblock)
1720
1721csp06/07/01 end
1722 endif
1723c if(icase.eq.5) then
1724c im5=im5+1
1725c write(65,*)'---im5,s2kaon,sig=',im5,brsig,sig
1726c call pipikaon(irun,iseed,dt,nt,ictrl,i1,i2,srt,pcx,pcy,pcz)
1727c endif
1728cbz3/2/99
1729c write(101,*)lb1,lb2,lb(i1),lb(i2)
1730c write(101,*)em1,em2,e(i1),e(i2),srt
1731cbz3/2/99end
1732
1733 return
1734 end
1735
1736******************************************
1737* for pp-->pp + kaon + anti-kaon
1738c 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******************************************
1746c 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-
1764c 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.
1778c the following data come from PRC 41 (1701)
1779c 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.
1790c the following data come from PRC 41 (1701)
1791c 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.
1801c 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.
1812c 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.
1822c 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.
1842c 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
1852c-----------------------------------------------------------------------
1853
1854c.....extracted from G. Song's ART expasion including K- interactions
1855c.....file `NEWNNK.FOR'
1856
1857 subroutine nnkaon(irun,iseed,ictrl,i1,i2,iblock,
1858 & srt,pcx,pcy,pcz,nchrg)
1859c <pt>=0.27+0.037*log(srt) was changed to 0.632 + ... on Aug. 14, 1997
1860c 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)
1864cc SAVE /AA/
1865 COMMON /BB/ P(3,MAXSTR)
1866cc SAVE /BB/
1867 COMMON /CC/ E(MAXSTR)
1868cc SAVE /CC/
1869 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1870cc SAVE /EE/
1871 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
1872cc SAVE /BG/
1873 COMMON /NN/NNN
1874cc SAVE /NN/
1875 COMMON /RUN/NUM
1876cc SAVE /RUN/
1877 COMMON /PA/RPION(3,MAXSTR,MAXR)
1878cc SAVE /PA/
1879 COMMON /PB/PPION(3,MAXSTR,MAXR)
1880cc SAVE /PB/
1881 COMMON /PC/EPION(MAXSTR,MAXR)
1882cc SAVE /PC/
1883 COMMON /PD/LPION(MAXSTR,MAXR)
1884cc 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
1890c dm1=e(i1)
1891c dm2=e(i2)
1892 dm3=0.938
1893 dm4=0.938
1894c 10/24/02 initialize n to 0:
1895 n=0
1896
1897cbz3/11/99 neutralk
1898c if(nchrg.eq.-2.or.nchrg.ge.3) dm3=1.232
1899c 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
1902cbz3/11/99 neutralk end
1903 iblock = 0
1904 call fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
1905 if(iflag.lt.0) then
1906c write(60,*)'------------final state fail-------',n
1907c 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
1920c 10/28/02 get rid of argument usage mismatch in rotate():
1921 pxrota=px(1)
1922 pyrota=py(1)
1923 pzrota=pz(1)
1924c 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
1929c
1930 pxrota=px(2)
1931 pyrota=py(2)
1932 pzrota=pz(2)
1933c 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
1938c
1939 pxrota=px(3)
1940 pyrota=py(3)
1941 pzrota=pz(3)
1942c 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
1947c
1948 pxrota=px(4)
1949 pyrota=py(4)
1950 pzrota=pz(4)
1951c 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
1958c K+
1959 lpion(nnn,irun)=23
1960 if(nchrg.eq.-1.or.nchrg.eq.-2) then
1961c To keep charge conservation. D-n->nnK0K-, D-D- -> nD-K0K-
1962
1963cbz3/7/99 neutralk
1964c lpion(nnn,irun)=24 ! K0
1965cbz3/7/99 neutralk end
1966
1967 endif
1968c aka: rest mass of K
1969 epion(nnn,irun)=aka
1970c K-
1971 lpion(nnn-1,irun)=21
1972c 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
1983c lb1 = lb(i1)
1984 lb1 = 2
1985 if(nchrg.ge.-2.and.nchrg.le.1) lb1=2
1986
1987cbz3/7/99 neutralk
1988 if (nchrg .eq. -2 .or. nchrg .eq. -1) then
1989 lb1 = 6
1990 end if
1991cbz3/7/99 neutralk end
1992
1993cbz3/11/99 neutralk
1994c if(nchrg.eq.2.or.nchrg.eq.3) lb1=1
1995c 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
1998cbz3/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
2008c lb2 = lb(i2)
2009 lb2 = 2
2010
2011cbz3/11/99 neutralk
2012c if(nchrg.eq.-1.or.nchrg.eq.0) lb2=2
2013c if(nchrg.eq. 2.or.nchrg.eq.1) lb2=1
2014c if(nchrg.eq. 4.or.nchrg.eq.3) lb2=9
2015c 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
2020cbz3/11/99 neutralk end
2021
2022c 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
2034c px1 = p(1,i1)
2035c py1 = p(2,i1)
2036c pz1 = p(3,i1)
2037c em1 = e(i1)
2038c id(i1) = 2
2039c id(i2) = 2
2040c id1 = id(i1)
2041c 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)
2052clin-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)
2064clin-5/2008:
2065 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2066 return
2067 end
2068
2069 subroutine lorntz(ilo,b,pi,pj)
2070c It uses to perform Lorentz (or inverse Lorentz) transformation
2071 dimension pi(4),pj(4),b(3)
2072 SAVE
2073c 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
2080c 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)
2083c drb=drd(1)*b(1)+drd(2)*b(2)+drd(3)*b(3)
2084c 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))
2088c drd(i)=drd(i)+b(i)*ga*drb
2089c 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
2094100 continue
2095c 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
2111cc SAVE /RNDF77/
2112 SAVE
2113
2114 iseed=iseed
2115 iflag=-1
2116c iflag=-1: fail to find momenta
2117c = 1: success
2118 pio=3.1415926
2119 aka=0.498
2120c v=0.43
2121c w=-0.84
2122c b=3.78
2123c c=0.47
2124c d=3.60
2125c fmax=1.056
2126c 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
2138c after we have the momenta for both nucleus, we sample the
2139c transverse momentum for K-.
2140c dsigma/dpt**2 = exp(-4.145*pt**2) obtained by fitting data on
2141c page 72, fig 23i.
214250 continue
2143 icount=icount+1
2144 if(icount.gt.10) return
2145 ptkmi2=-1./4.145*alog(RANART(NSEED))
2146 ptkm=sqrt(ptkmi2)
21473 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
2164150 ptkpl2=-1./3.68*alog(RANART(NSEED))
2165 ptkp=sqrt(ptkpl2)
216613 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)
2186c resample
2187 if(resten.le.abs(restpz)) goto 50
2188 restms=sqrt(resten**2-restpz**2)
2189c 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))
2198c transverse mass for K-
2199 rmt3=sqrt(dm3**2+px(1)**2+py(1)**2)
2200c transverse mass for K+
2201 rmt4=sqrt(dm4**2+px(2)**2+py(2)**2)
2202 if(restms.lt.(rmt3+rmt4)) goto 50
2203c 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
2224c-----------------------------------------------------------------------
2225
2226c.....extracted from G. Song's ART expasion including K- interactions
2227c.....file `NPIK.FOR'
2228
2229****************************************
2230c subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2231c & 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)
2248cc SAVE /AA/
2249 COMMON /BB/ P(3,MAXSTR)
2250cc SAVE /BB/
2251 COMMON /CC/ E(MAXSTR)
2252cc SAVE /CC/
2253 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2254cc SAVE /EE/
2255 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
2256cc SAVE /BG/
2257 COMMON /NN/NNN
2258cc SAVE /NN/
2259 COMMON /RUN/NUM
2260cc SAVE /RUN/
2261 COMMON /PA/RPION(3,MAXSTR,MAXR)
2262cc SAVE /PA/
2263 COMMON /PB/PPION(3,MAXSTR,MAXR)
2264cc SAVE /PB/
2265 COMMON /PC/EPION(MAXSTR,MAXR)
2266cc SAVE /PC/
2267 COMMON /PD/LPION(MAXSTR,MAXR)
2268cc SAVE /PD/
2269 dimension bb(3),p1(4),p2(4),p3(4),px(4),py(4),pz(4)
2270 COMMON/RNDF77/NSEED
2271cc 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
2292c 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
2297cbz3/8/99 neutralk
2298cbz10/12/99
2299c LB(I1) = 1 + 2 * RANART(NSEED)
2300c LB(I2) = 23
2301 LB(k1) = 1 + int(2*RANART(NSEED))
2302 LB(k2) = 23
2303c pkmax=sqrt((srt**2-(aka+0.938+aka)**2)*(srt**2-(aka+0.938-aka)**2))
2304c & /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
2308c-----------------------------------------------------
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
2320c 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
2323c-----------------------------------------------------
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
2335c p1,p2: the momenta of the two particles in their cms
2336c p1: momentum of N or P
2337c p2: momentum of anti_kaon
2338c p3: momentum of K0 or K+
2339 ilo=1
2340c write(61,*)'--------p1,p2',p1,p2
2341c write(61,*)'--------bb',bb
2342 call lorntz(ilo,bb,p1,p2)
2343c******* Checking *************
2344c pxsum = p1(1)+p2(1)+p3(1)
2345c pysum = p1(2)+p2(2)+p3(2)
2346c pzsum = p1(3)+p2(3)+p3(3)
2347c pesum = p1(4)+p2(4)+sqrt(p3(1)**2+p3(2)**2+p3(3)**2+aka**2)-srt
2348c write(61,*)'---p1,pxsum',p1,pxsum
2349c write(61,*)'---p2,pysum',p2,pysum
2350c write(61,*)'---p3,pzsum',p3,pzsum
2351c write(61,*)'---pesum',pesum
2352c***********************************
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
2359c 10/28/02 get rid of argument usage mismatch in rotate():
2360 pxrota=p1(1)
2361 pyrota=p1(2)
2362 pzrota=p1(3)
2363c 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
2368c
2369 pxrota=p2(1)
2370 pyrota=p2(2)
2371 pzrota=p2(3)
2372c 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
2377c
2378 pxrota=p3(1)
2379 pyrota=p3(2)
2380 pzrota=p3(3)
2381c 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
2388c K(-)
2389 lpion(nnn,irun)=21
2390c 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
2413c 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
2426c px1 = p(1,i1)
2427c py1 = p(2,i1)
2428c pz1 = p(3,i1)
2429c em1 = e(i1)
2430c id(i1) = 2
2431c id(i2) = 2
2432c id1 = id(i1)
2433c K(+)K(-) production
2434 iblock = 101
2435* Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2436c p2: momentum of anti-kaon.
2437c 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)
2444clin-5/2008:
2445 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2446cbz3/2/99
2447c write(400,*)'2 ', ppion(1,nnn,irun), ppion(2,nnn,irun),
2448c & ppion(3,nnn,irun), dt*nt, srt
2449cbz3/2/99end
2450c write(420,*)ppion(1,nnn,irun), ppion(2,nnn,irun),
2451c & 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
2460c-----------------------------------------------------------------------
2461
2462c.....extracted from G. Song's ART expasion including K- interactions
2463c.....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
2473c 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)
2477cc SAVE /AA/
2478 COMMON /BB/ P(3,MAXSTR)
2479cc SAVE /BB/
2480 COMMON /CC/ E(MAXSTR)
2481cc SAVE /CC/
2482 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2483cc SAVE /EE/
2484 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
2485cc SAVE /BG/
2486 COMMON /NN/NNN
2487cc SAVE /NN/
2488 COMMON /RUN/NUM
2489cc SAVE /RUN/
2490 COMMON /PA/RPION(3,MAXSTR,MAXR)
2491cc SAVE /PA/
2492 COMMON /PB/PPION(3,MAXSTR,MAXR)
2493cc SAVE /PB/
2494 COMMON /PC/EPION(MAXSTR,MAXR)
2495cc SAVE /PC/
2496 COMMON /PD/LPION(MAXSTR,MAXR)
2497cc SAVE /PD/
2498 dimension p1(4),p2(4)
2499 COMMON/RNDF77/NSEED
2500cc 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
2510csp06/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
2520csp06/07/01 end
2521c PI + Sigma(or Lambda) -> Kbar + N
2522 k1=i1
2523 k2=i2
2524c 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
2529cbz3/8/99 neutralk
2530 LB(K1) = 1 + int(2*RANART(NSEED))
2531 if(nchrg.eq.-2) lb(k1)=6
2532c if(nchrg.eq.-1) lb(k1)=2
2533c if(nchrg.eq. 0) lb(k1)=1
2534c if(nchrg.eq. 1) lb(k1)=9
2535 IF (NCHRG .EQ. 2) LB(K1) = 9
2536cbz3/8/99 neutralk end
2537
2538c K-
2539 lb(k2)=21
2540 dm3=0.938
2541 if(nchrg.eq.-2.or.nchrg.eq.1) dm3=1.232
2542 dm4=aka
2543c dm3,dm4: the mass of final state particles.
2544 endif
2545
2546********Now, antikaon will be created.
2547c call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2548c 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
2552c-----------------------------------------------------
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
2562c p1,p2: the momenta of the two particles in their cms
2563c p1: momentum of kaon
2564c p2: momentum of Kbar
2565
2566* Rotate the momenta of particles in the cms of I1 & I2
2567clin-10/28/02 get rid of argument usage mismatch in rotate():
2568 pxrota=p1(1)
2569 pyrota=p1(2)
2570 pzrota=p1(3)
2571c 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
2576c
2577 pxrota=p2(1)
2578 pyrota=p2(2)
2579 pzrota=p2(3)
2580c 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
2585clin-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)
2605cbz3/2/99
2606c write(400,*)'3 ', pt1i2, pt2i2, pt3i2, dt*nt, srt
2607cbz3/2/99end
2608c write(430,*)pt1i2, pt2i2, pt3i2, dt*nt, srt
2609 eti2 = dm4
2610 lb2 = lb(k2)
2611
2612c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2613c k1=i1
2614c 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
2627cc iblock = 101 ! K(+)K(-) production
2628* Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2629 return
2630 end
2631
2632c-----------------------------------------------------------------------
2633
2634c.....extracted from G. Song's ART expasion including K- interactions
2635c.....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)
2647cc SAVE /AA/
2648 COMMON /BB/ P(3,MAXSTR)
2649cc SAVE /BB/
2650 COMMON /CC/ E(MAXSTR)
2651cc SAVE /CC/
2652 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2653cc SAVE /EE/
2654 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
2655cc SAVE /BG/
2656 COMMON /NN/NNN
2657cc SAVE /NN/
2658 COMMON /RUN/NUM
2659cc SAVE /RUN/
2660 COMMON /PA/RPION(3,MAXSTR,MAXR)
2661cc SAVE /PA/
2662 COMMON /PB/PPION(3,MAXSTR,MAXR)
2663cc SAVE /PB/
2664 COMMON /PC/EPION(MAXSTR,MAXR)
2665cc SAVE /PC/
2666 COMMON /PD/LPION(MAXSTR,MAXR)
2667cc SAVE /PD/
2668 dimension p1(4),p2(4)
2669 COMMON/RNDF77/NSEED
2670cc 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
2681c ratio: used for isospin decision.
2682 k1=i1
2683 k2=i2
2684c 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 *******************************
2690c record kaon's mass
2691 eee=e(k2)
2692*** end **************
2693 rrr=RANART(NSEED)
2694 if(rrr.lt.brel) then
2695c 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
2704c nchrg: Net charges of the two incoming particles.
2705c Kbar + N -> Sigma + PI
2706 em1=asa
2707 em2=0.138
2708
2709cbz3/8/99 neutralk
2710 LB1 = 15 + int(3*RANART(NSEED))
2711 LB2 = 3 + int(3*RANART(NSEED))
2712 else
2713c Kbar + N -> Lambda + PI
2714 em1=ala
2715 em2=0.138
2716c LAmbda
2717 lb1=14
2718cbz3/8/99 neutralk
2719 LB2 = 3 + int(3*RANART(NSEED))
2720c if(nchrg.eq.1) lb2=5 ! K- + D++ -> Lambda + PI+
2721c if(nchrg.eq.0) lb2=4 ! K- + p(D+,N*+) -> Lambda + PI0
2722c if(nchrg.eq.-1) lb2=3 ! K- + n(D,N*) -> Lambda + PI-
2723cbz3/8/99 neutralk
2724
2725 endif
2726 endif
2727 lb(k1)=lb1
2728 lb(k2)=lb2
2729
2730********Now, antikaon will be created.
2731c call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2732c pkmax: the maximum momentum of anti-kaon
2733c write(63,*)'srt,em1,em2',srt,em1,em2
2734c 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
2738c-----------------------------------------------------
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
2748c p1,p2: the momenta of the two particles in their cms
2749c p1: momentum of kaon
2750c p2: momentum of Kbar
2751
2752* Rotate the momenta of particles in the cms of I1 & I2
2753
2754clin-10/28/02 get rid of argument usage mismatch in rotate():
2755 pxrota=p1(1)
2756 pyrota=p1(2)
2757 pzrota=p1(3)
2758c 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
2763c
2764 pxrota=p2(1)
2765 pyrota=p2(2)
2766 pzrota=p2(3)
2767c 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
2772clin-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
2793c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2794c k1=i1
2795c 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
2806cc iblock = 101 ! K(+)K(-) production
2807* Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2808 return
2809 end
2810
2811c=======================================================================
2812
2813clin Below is the previous artana.f:
2814c=======================================================================
2815
2816c.....analysis subroutine before the hadronic space-time evolution
2817
2818 SUBROUTINE ARTAN1
2819 PARAMETER (MAXSTR=150001, MAXR=1)
2820c.....y cut for mt spectrum
2821cbz3/17/99
2822c PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
2823 PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
2824cbz3/17/99 end
2825c.....bin width for mt spectrum and y spectrum
2826clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
2827c PARAMETER (BMT = 0.05, BY = 0.2)
2828 PARAMETER (BMT = 0.05, BY = 0.4)
2829 COMMON /RUN/ NUM
2830cc SAVE /RUN/
2831 COMMON /ARERC1/MULTI1(MAXR)
2832cc 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)
2838cbz3/17/99
2839c & dm1k0s(50), DMT1LA(50), DMT1LB(50)
2840cc 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)
2851cc SAVE /ARANA1/
2852 SAVE
2853
2854cbz3/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)
2863c 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'
2872cbzdbg2/16/99
2873 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
2874cbzdbg2/16/99
2875cbzdbg2/15/99
2876 PRINT *, ' PZ = ', PZ, ' EE = ', EE
2877cbzdbg2/16/99
2878 PRINT *, ' XM = ', XM
2879cbzdbg2/16/99end
2880 GOTO 200
2881c STOP
2882cbzdbg2/15/99end
2883 END IF
2884 DXMT = XMT - XM
2885 Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
2886c.....rapidity cut for the rapidity distribution
2887 IF (ABS(Y) .GE. 10.0) GOTO 100
2888clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
2889c IY = 1 + int(ABS(Y) / BY)
2890c 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
2907c IF (ITYP .EQ. -211 .OR. ITYP .EQ. -321 .OR.
2908c & ITYP .EQ. -2212) THEN
2909 IF (ITYP .EQ. -2112) THEN
2910 DY1HM(IY) = DY1HM(IY) + 1.0
2911 END IF
2912c
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
2922cbz3/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
2942cbz3/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
2949clin-4/24/03 evaluate K0L instead of K0S, since sometimes we may decay K0S:
2950c 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
2964c.....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
2987clin-4/24/03:
2988c 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
3006c-----------------------------------------------------------------------
3007
3008c.....analysis subroutine after the hadronic space-time evolution
3009
3010 SUBROUTINE ARTAN2
3011
3012 PARAMETER (MAXSTR=150001, MAXR=1)
3013c.....y cut for mt spectrum
3014cbz3/17/99
3015c PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3016 PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3017cbz3/17/99 end
3018c.....bin width for mt spectrum and y spectrum
3019c PARAMETER (BMT = 0.05, BY = 0.2)
3020 PARAMETER (BMT = 0.05, BY = 0.4)
3021 COMMON /RUN/ NUM
3022cc SAVE /RUN/
3023 COMMON /ARERC1/MULTI1(MAXR)
3024cc 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)
3030cbz3/17/99
3031c & dm2k0s(50), DMT2LA(50), DMT2LB(50)
3032cc 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)
3043cbz3/17/99 end
3044cc 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)
3056c 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'
3064cbzdbg2/16/99
3065 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3066cbzdbg2/16/99
3067cbzdbg2/15/99
3068 PRINT *, ' PZ = ', PZ, ' EE = ', EE
3069cbzdbg2/16/99
3070 PRINT *, ' XM = ', XM
3071cbzdbg2/16/99end
3072 GOTO 200
3073c STOP
3074cbzdbg2/15/99end
3075 END IF
3076 DXMT = XMT - XM
3077 Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
3078c.....rapidity cut for the rapidity distribution
3079 IF (ABS(Y) .GE. 10.0) GOTO 100
3080c IY = 1 + int(ABS(Y) / BY)
3081c 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
3111cbz3/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
3131cbz3/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
3138clin-4/24/03:
3139c 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
3153c.....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
3176clin-4/24/03:
3177c 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
3195c-----------------------------------------------------------------------
3196
3197c.....output analysis results at the end of the simulation
3198
3199 SUBROUTINE ARTOUT(NEVNT)
3200
3201 PARAMETER (MAXSTR=150001, MAXR=1)
3202c.....y cut for mt spectrum
3203cbz3/17/99
3204c PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3205 PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3206cbz3/17/99 end
3207c.....bin width for mt spectrum and y spectrum
3208c PARAMETER (BMT = 0.05, BY = 0.2)
3209 PARAMETER (BMT = 0.05, BY = 0.4)
3210 COMMON /RUN/ NUM
3211cc 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)
3217cbz3/17/99
3218c & dm1k0s(50), DMT1LA(50), DMT1LB(50)
3219cc 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)
3230cbz3/17/99 end
3231cc SAVE /ARANA1/
3232cbz3/17/99
3233c & 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)
3244cc SAVE /ARANA2/
3245 SAVE
3246cbz3/17/99 end
3247cms OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN')
3248cms OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN')
3249cms OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN')
3250cms OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN')
3251cms OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN')
3252cms OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN')
3253cms OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN')
3254cms OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN')
3255cms OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN')
3256clin-4/24/03
3257c OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN')
3258cms OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN')
3259cms OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN')
3260cms OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN')
3261cms OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN')
3262cbz3/17/99
3263cms OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN')
3264cms OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN')
3265cms OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN')
3266cms OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN')
3267cms OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN')
3268cms OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN')
3269cbz3/17/99 end
3270
3271cms OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN')
3272cms OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN')
3273cms OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN')
3274cms OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN')
3275cms OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN')
3276cms OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN')
3277cms OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN')
3278cms OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN')
3279cms OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN')
3280clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
3281c SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0
3282 SCALE1 = 1. / REAL(NEVNT * NUM) / BY
3283 SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1)
3284c
3285 DO 1001 I = 1, 50
3286 ymid=-10.+BY * (I - 0.5)
3287cms WRITE (30, 333) ymid, SCALE1 * dy1ntb(I)
3288cms WRITE (31, 333) ymid, SCALE1 * dy1ntp(I)
3289cms WRITE (32, 333) ymid, SCALE1 * DY1HM(I)
3290cms WRITE (37, 333) ymid, SCALE1 * DY1KP(I)
3291cms WRITE (38, 333) ymid, SCALE1 * DY1KM(I)
3292cms WRITE (39, 333) ymid, SCALE1 * DY1K0S(I)
3293cms WRITE (40, 333) ymid, SCALE1 * DY1LA(I)
3294cms WRITE (41, 333) ymid, SCALE1 * DY1LB(I)
3295cms WRITE (42, 333) ymid, SCALE1 * DY1PHI(I)
3296cms WRITE (33, 333) ymid, SCALE1 * DY1NEG(I)
3297cms WRITE (34, 333) ymid, SCALE1 * DY1CH(I)
3298cms WRITE (35, 333) ymid, SCALE1 * DE1NEG(I)
3299cms WRITE (36, 333) ymid, SCALE1 * DE1CH(I)
3300cms WRITE (43, 333) ymid, SCALE1 * dy1msn(I)
3301cms WRITE (44, 333) ymid, SCALE1 * DY1PIP(I)
3302cms WRITE (45, 333) ymid, SCALE1 * DY1PIM(I)
3303cms WRITE (46, 333) ymid, SCALE1 * DY1PI0(I)
3304cms WRITE (47, 333) ymid, SCALE1 * DY1PR(I)
3305cms WRITE (48, 333) ymid, SCALE1 * DY1PB(I)
3306
3307 IF (dm1pip(I) .NE. 0.0) THEN
3308cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I)
3309 END IF
3310 IF (dm1pim(I) .NE. 0.0) THEN
3311cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 *
3312cms & dm1pim(I)
3313 END IF
3314 IF (DMT1PR(I) .NE. 0.0) THEN
3315cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I)
3316 END IF
3317 IF (DMT1PB(I) .NE. 0.0) THEN
3318cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I)
3319 END IF
3320 IF (DMT1KP(I) .NE. 0.0) THEN
3321cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I)
3322 END IF
3323 IF (dm1km(I) .NE. 0.0) THEN
3324cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3325cms & dm1km(I)
3326 END IF
3327 IF (dm1k0s(I) .NE. 0.0) THEN
3328cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I)
3329 END IF
3330 IF (DMT1LA(I) .NE. 0.0) THEN
3331cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I)
3332 END IF
3333 IF (DMT1LB(I) .NE. 0.0) THEN
3334cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I)
3335 END IF
3336 1001 CONTINUE
3337c
3338 DO 1002 I = 30, 48
3339cms WRITE (I, *) 'after hadron evolution'
3340 1002 CONTINUE
3341 DO 1003 I = 50, 58
3342cms WRITE (I, *) 'after hadron evolution'
3343 1003 CONTINUE
3344
3345 DO 1004 I = 1, 50
3346 ymid=-10.+BY * (I - 0.5)
3347cms WRITE (30, 333) ymid, SCALE1 * dy2ntb(I)
3348cms WRITE (31, 333) ymid, SCALE1 * dy2ntp(I)
3349cms WRITE (32, 333) ymid, SCALE1 * DY2HM(I)
3350cms WRITE (37, 333) ymid, SCALE1 * DY2KP(I)
3351cms WRITE (38, 333) ymid, SCALE1 * DY2KM(I)
3352cms WRITE (39, 333) ymid, SCALE1 * DY2K0S(I)
3353cms WRITE (40, 333) ymid, SCALE1 * DY2LA(I)
3354cms WRITE (41, 333) ymid, SCALE1 * DY2LB(I)
3355cms WRITE (42, 333) ymid, SCALE1 * DY2PHI(I)
3356cms WRITE (33, 333) ymid, SCALE1 * DY2NEG(I)
3357cms WRITE (34, 333) ymid, SCALE1 * DY2CH(I)
3358cms WRITE (35, 333) ymid, SCALE1 * DE2NEG(I)
3359cms WRITE (36, 333) ymid, SCALE1 * DE2CH(I)
3360cms WRITE (43, 333) ymid, SCALE1 * dy2msn(I)
3361cms WRITE (44, 333) ymid, SCALE1 * DY2PIP(I)
3362cms WRITE (45, 333) ymid, SCALE1 * DY2PIM(I)
3363cms WRITE (46, 333) ymid, SCALE1 * DY2PI0(I)
3364cms WRITE (47, 333) ymid, SCALE1 * DY2PR(I)
3365cms WRITE (48, 333) ymid, SCALE1 * DY2PB(I)
3366c
3367 IF (dm2pip(I) .NE. 0.0) THEN
3368cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I)
3369 END IF
3370 IF (dm2pim(I) .NE. 0.0) THEN
3371cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 *
3372cms & dm2pim(I)
3373 END IF
3374 IF (DMT2PR(I) .NE. 0.0) THEN
3375cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I)
3376 END IF
3377 IF (DMT2PB(I) .NE. 0.0) THEN
3378cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I)
3379 END IF
3380 IF (DMT2KP(I) .NE. 0.0) THEN
3381cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I)
3382 END IF
3383 IF (dm2km(I) .NE. 0.0) THEN
3384cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3385cms & dm2km(I)
3386 END IF
3387 IF (dm2k0s(I) .NE. 0.0) THEN
3388cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I)
3389 END IF
3390 IF (DMT2LA(I) .NE. 0.0) THEN
3391cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I)
3392 END IF
3393 IF (DMT2LB(I) .NE. 0.0) THEN
3394cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I)
3395 END IF
3396 1004 CONTINUE
3397cms 333 format(2(f12.5,1x))
3398
3399 RETURN
3400 END
3401
3402c-----------------------------------------------------------------------
3403
3404c.....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
3424cc SAVE /PARA1/
3425 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3426cc SAVE /HPARNT/
3427 COMMON/hjcrdn/YP(3,300),YT(3,300)
3428cc 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)
3434cc 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)
3438cc 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)
3442cc SAVE /prec1/
3443 COMMON /AREVT/ IAEVT, IARUN, MISS
3444cc SAVE /AREVT/
3445 COMMON /AROUT/ IOUT
3446cc 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
3497c.....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)
3615clin-4/2008 protect against out-of-bound errors:
3616c 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
3671c.....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
3694cbzdbg2/16/99
3695c PRINT *, ' in HJANA1 '
3696c PRINT *, ' total number of partons = ', nsubp
3697c PRINT *, ' total number of gluons = ', nsubg, MUL
3698c PRINT *, ' number of projectile strings = ', IHNT2(1)
3699c PRINT *, ' number of target strings = ', IHNT2(3)
3700c 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
3704c PRINT *, ' number of projectile strings = ', IHNT2(1)
3705c PRINT *, ' number of target strings = ', IHNT2(3)
3706 PRINT *, ' number of independent strings = ', NISG / IW
3707cbzdbg2/16/99end
3708 END IF
3709c
3710 RETURN
3711 END
3712
3713c-----------------------------------------------------------------------
3714
3715c.....analysis subroutine in ZPC after generation of additional initial
3716c.....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
3725cc 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)
3729cc SAVE /prec2/
3730 COMMON /AREVT/ IAEVT, IARUN, MISS
3731cc SAVE /AREVT/
3732 COMMON /AROUT/ IOUT
3733cc 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
3753c.....analysis
3754 DO 1003 I = 1, MUL
3755 IGX = 1 + int(sngl(ABS(GX5(I))) / DGX)
3756clin-4/2008 protect against out-of-bound errors:
3757c 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
3771c
3772 RETURN
3773 END
3774
3775c-----------------------------------------------------------------------
3776
3777c.....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
3786cc 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)
3790cc SAVE /prec2/
3791 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
3792cc SAVE /ilist8/
3793 COMMON /SREC1/ NSP, NST, NSI
3794cc SAVE /SREC1/
3795 COMMON/hjcrdn/YP(3,300),YT(3,300)
3796cc 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)
3800cc SAVE /HJJET2/
3801 COMMON /AREVT/ IAEVT, IARUN, MISS
3802cc SAVE /AREVT/
3803 COMMON /AROUT/ IOUT
3804cc 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
3822c.....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
3850c
3851 RETURN
3852 END
3853
3854c-----------------------------------------------------------------------
3855
3856c.....analysis subroutine in HIJING after parton cascade evolution
3857 SUBROUTINE HJANA2
3858c
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
3883cc SAVE /PARA1/
3884 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3885cc SAVE /HPARNT/
3886 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
3887cc 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)
3893cc 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)
3897cc 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)
3901cc SAVE /prec2/
3902 COMMON /AREVT/ IAEVT, IARUN, MISS
3903cc SAVE /AREVT/
3904 COMMON /AROUT/ IOUT
3905cc SAVE /AROUT/
3906 common/anim/nevent,isoft,isflag,izpc
3907cc 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)
3912cc 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
3973clin-4/28/01:
3974 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 510
3975
3976c.....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)
3985cbzdbg2/16/99
3986c 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
3994cbzdbg2/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)
4026cbzdbg2/16/99
4027c 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
4035cbzdbg2/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
4059clin-4/28/01:
4060 510 continue
4061
4062 DO 1010 I = 1, NSG
4063clin-4/25/01 soft3:
4064c 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
4068clin-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)
4076clin-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
4085clin-4/25/01-end
4086
4087cbzdbg2/16/99
4088c 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
4096cbzdbg2/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
4120clin-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
4153clin-4/28/01:
4154 520 continue
4155
4156 DO 1013 I = 1, NSG
4157 J = I + IHNT2(1) + IHNT2(3)
4158clin-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))
4181cbzdbg2/16/99
4182c 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
4192cbzdbg2/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
4208clin-4/25/01 soft3:
4209 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 530
4210
4211c.....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
4226clin-4/25/01 soft3:
4227 530 continue
4228
4229 DO 1020 I = 1, NSG
4230clin-4/25/01 soft3:
4231c 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
4235clin-4/25/01-end
4236
4237 nsubp = nsubp + 1
4238
4239clin-4/25/01
4240c 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
4246clin-4/25/01-end
4247 1019 CONTINUE
4248 1020 CONTINUE
4249cbzdbg2/16/99
4250 NISG = NISG + NSG
4251
4252 IF (IOUT .EQ. 1) THEN
4253cbzdbg2/16/99end
4254cbzdbg2/16/99
4255c PRINT *, ' in HJANA2 '
4256c PRINT *, ' total number of partons = ', nsubp
4257c PRINT *, ' total number of gluons = ', nsubg, MUL
4258c PRINT *, ' number of projectile strings = ', IHNT2(1)
4259c PRINT *, ' number of target strings = ', IHNT2(3)
4260c 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
4264c PRINT *, ' number of projectile strings = ', IHNT2(1)
4265c 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
4275c-----------------------------------------------------------------------
4276
4277c.....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
4288cc 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)
4292cc SAVE /prec2/
4293 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4294cc SAVE /HPARNT/
4295 COMMON/hjcrdn/YP(3,300),YT(3,300)
4296cc 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)
4302cc 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)
4306cc SAVE /HJJET2/
4307 COMMON /AREVT/ IAEVT, IARUN, MISS
4308cc SAVE /AREVT/
4309 COMMON /AROUT/ IOUT
4310cc 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
4336c.....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
4407c
4408 RETURN
4409 END
4410
4411c-----------------------------------------------------------------------
4412
4413c.....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
4425cc 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)
4429cc SAVE /prec2/
4430 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
4431cc SAVE /ilist8/
4432 COMMON /SREC1/ NSP, NST, NSI
4433cc SAVE /SREC1/
4434 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
4435cc SAVE /SREC2/
4436 COMMON/hjcrdn/YP(3,300),YT(3,300)
4437cc 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)
4441cc SAVE /HJJET2/
4442 COMMON /AREVT/ IAEVT, IARUN, MISS
4443cc SAVE /AREVT/
4444 COMMON /AROUT/ IOUT
4445cc 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
4463c.....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)
4476cbzdbg2/21/99
4477c IF (ABS(IT) .GT. 25) GOTO 200
4478 IF (IT .GT. 25 .OR. IT .LT. -24) GOTO 200
4479cbzdbg2/21/99end
4480 dtg2b(IT) = dtg2b(IT) + 1.0
4481 200 CONTINUE
4482 1003 CONTINUE
4483c
4484 RETURN
4485 END
4486
4487c-----------------------------------------------------------------------
4488
4489c.....analysis subroutine before ARTMN
4490 SUBROUTINE HJANA3
4491c
4492 PARAMETER (MAXSTR=150001, MAXR=1)
4493c.....y cut for mt spectrum
4494 PARAMETER (YMIN = -1.0, YMAX = 1.0)
4495cbz11/7/99 end
4496c.....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
4502cc SAVE /RUN/
4503 COMMON /ARERC1/MULTI1(MAXR)
4504cc 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)
4510cc SAVE /ARPRC1/
4511 COMMON /AROUT/ IOUT
4512cc SAVE /AROUT/
4513 COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
4514cc 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))
4539c.....rapidity cut for the rapidity distribution
4540c 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
4546c.....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
4554c
4555 RETURN
4556 END
4557
4558c-----------------------------------------------------------------------
4559
4560c.....analysis subroutine after ARTMN
4561 SUBROUTINE HJANA4
4562 PARAMETER (MAXSTR=150001, MAXR=1)
4563c.....y cut for mt spectrum
4564cbz11/7/99
4565c PARAMETER (YMIN = -0.5, YMAX = 0.5)
4566 PARAMETER (YMIN = -1.0, YMAX = 1.0)
4567cbz11/7/99 end
4568c.....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
4573cc SAVE /RUN/
4574 COMMON /ARERC1/MULTI1(MAXR)
4575cc 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)
4581cc SAVE /ARPRC1/
4582 COMMON /AROUT/ IOUT
4583cc SAVE /AROUT/
4584 COMMON /fflow/ v2f,etf,xmultf,v2fpi,xmulpi
4585cc 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))
4610c.....rapidity cut for the rapidity distribution
4611c 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
4617c.....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
4625c
4626 RETURN
4627 END
4628
4629c=======================================================================
4630
4631c.....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)
4638c 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
4642cc 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)
4646cc SAVE /prec2/
4647 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
4648cc SAVE /ilist8/
4649 COMMON /SREC1/ NSP, NST, NSI
4650cc SAVE /SREC1/
4651 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
4652cc SAVE /SREC2/
4653 COMMON/hjcrdn/YP(3,300),YT(3,300)
4654cc 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)
4658cc SAVE /HJJET2/
4659cbz6/28/99 flow1
4660 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4661cc SAVE /HPARNT/
4662cbz6/28/99 flow1 end
4663 common/anim/nevent,isoft,isflag,izpc
4664cc SAVE /anim/
4665 common/strg/np(maxstr)
4666cc SAVE /strg/
4667clin-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
4673cc SAVE /frzprc/
4674 SAVE
4675
4676clin-6/06/02 test local freezeout for string melting,
4677c 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
4692clin-6/06/02-end
4693
4694 DO 1002 I = 1, MAXSTR
4695 ATAUI(I) = 0d0
4696 ZT1(I) = 0d0
4697 ZT2(I) = 0d0
4698clin-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
4714clin-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
4726clin-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
4756cbz6/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
4770cbz6/28/99 flow1 end
4771c
4772 RETURN
4773 END
4774
4775clin-3/2009
4776c Initialize hadron weights;
4777c 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
4793c All hadrons at the start of hadron cascade have the weight of 1
4794c except those inserted by the user in this subroutine:
4795 np0=IAINT2(1)
4796 DO i=1,np0
4797 dpertp(I)=1.
4798 ENDDO
4799c 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
4813c
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)
4818c.....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)
4823c Allow the intial z-position to be different from the Bjorken picture:
4824 GZAR(I)=TAUI*SINH(RAP)+GZAR(I)
4825c 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
4830c
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
4840c
4841 return
4842 end