1 c....................amptsub.f
2 c.....this file contains 4 sections:
3 c.....1. ART subroutines;
4 c.....2. ART functions;
5 c.....3. ART block data;
6 c.....4. subprocesses borrowed from other codes.
7 c.....5. the previous artana.f
8 c.....6. the previous zpcsub.f
9 c.....7. subroutine getnp
10 c.....Note that Parts1-4 are the previous artsub.f
12 c=======================================================================
13 c.....subroutine to set up ART parameters and analysis files
14 c.....before looping different events
16 cms dlw & gsfs 8/2009 commented out lots of output files
20 PARAMETER (AMU= 0.9383)
21 double precision dpcoal,drcoal,ecritl
23 common /gg/ dx,dy,dz,dpx,dpy,dpz
25 c "SAVE " (without argument) is used for most subroutines and functions,
26 c this is important for the success when using "f77" to compile:
32 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
34 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
35 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
37 COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
39 common /imulst/ iperts
41 common /coal/dpcoal,drcoal,ecritl
42 common/anim/nevent,isoft,isflag,izpc
43 common /para7/ ioscar,nsmbbbar,nsmmeson
45 clin-10/03/03 ecritl: local energy density below which a parton
46 c will freeze out (in GeV/fm^3), for improvements on string melting,
47 c not used in this version of AMPT:
52 c combine ART initialization into ampt.ini:
53 c (Note that the following values are relics from the old ART structure)
54 c.....input parameter file
55 c OPEN(13, FILE = 'art1.ini', STATUS = 'UNKNOWN')
56 c READ (13, *) MASSTA, ZTA
59 c write(12,*) massta, zta, ' massta, zta'
60 c READ (13, *) MASSPR, ZPR
63 c write(12,*) masspr, zpr, ' masspr, zpr'
64 c READ (13, *) PLAB, IPLAB
67 c write(12,*) plab, iplab, ' plab, iplab'
69 elab=sqrt(plab**2+amu**2)-amu
76 c write(12,*) zeropt, ' zeropt'
77 clin-10/03/03 ISEED was used as a seed for random number inside ART,
80 c 0/1: (Normal or Perturbative) multistrange partice production.
81 c Perturbative option is disabled for now:
83 c READ (13, *) MANYB, B0, BI, BM
84 c 2/04/00 MANYB MUST BE SET TO 1 !
85 c in order to skip impact parameter setting by ART, then B0 has no effect.
90 c write(12,*) manyb, b0, bi, bm, ' manyb, b0, bi, bm'
92 c write(12,*) iseed, ' iseed'
94 c write(12,*) dt, ' dt'
96 c write(12,*) ntmax, ' ntmax'
99 c write(12,*) icoll, ' icoll'
101 c 2/11/03 run events without test particles for now:
103 c write(12,*) num, ' num'
106 c write(12,*) insys, ' insys'
109 c write(12,*) ipot, ' ipot'
112 IF(ICOLL.EQ.-1)IPOT=0
113 c write(12,*) mode, ' mode'
114 c READ (13, *) DX, DY, DZ
118 c write(12,*) dx,dy,dz,' dx,dy,dz'
119 c READ (13, *) DPX, DPY, DPZ
123 c write(12,*) dpx,dpy,dpz,' dpx,dpy,dpz'
124 c READ (13, *) IAVOID
126 c write(12,*) iavoid, ' iavoid'
127 c READ (13, *) IMOMEN
129 c write(12,*) imomen, ' imomen'
130 if(icoll.eq.-1)imomen=3
133 c write(12,*) nfreq, ' nfreq'
134 c READ (13, *) ICFLOW
136 c write(12,*) ICFLOW, ' ICFLOW'
139 c write(12,*) ICRHO, ' ICRHO'
142 c write(12,*)icou, ' icou'
143 * kaon potential control parameter
144 * KMUL IS A MULTIPLIER TO THE STANDARD K-N SCATTERING LENGTH
145 c READ (13, *) KPOTEN, KMUL
148 c write(12,*)kpoten,kmul, ' kpoten, kmul'
149 * mean field control parameter FOR BARYONS
150 * no mean filed is used for baryons if their
151 * local density is higher than dencut.
152 c READ (13, *) DENCUT
154 c write(12,*)dencut, ' dencut'
155 * test reactions in a box of side-length cycbox
157 c READ (13, *) CYCBOX
159 c write(12,*) cycbox, ' cycbox'
162 c if(ioscar.eq.2) then
163 if(ioscar.eq.2.or.ioscar.eq.3) then
164 cms OPEN (92,FILE='ana/parton-initial-afterPropagation.dat',
165 cms 1 STATUS = 'UNKNOWN')
168 clin-6/2009 write out full parton collision history:
169 cms OPEN (95,FILE='ana/parton-collisionsHistory.dat',
170 cms 1 STATUS='UNKNOWN')
171 clin-6/2009 write out initial minijet information:
172 cms OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat',
173 cms 1 STATUS='UNKNOWN')
174 clin-6/2009 write out parton info after coalescence:
175 if(isoft.eq.4.or.isoft.eq.5) then
176 cms OPEN (85,FILE='ana/parton-after-coalescence.dat',
177 cms 1 STATUS='UNKNOWN')
180 clin-6/2009 write out initial transverse positions of initial nucleons:
181 cms OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN')
186 c-----------------------------------------------------------------------
188 c.....subroutine to initialize cascade.
192 c.....before invoking ARINI:
193 c.....IAPAR2(1), IAINT2(1) must be set.
194 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
198 ctest off for resonance (phi, K*) studies:
199 c OPEN (89, FILE = 'ana/decay_rec.dat', STATUS = 'UNKNOWN')
202 GOTO (200, 200, 300) IFLG
204 c.....error choice of initialization
205 PRINT *, 'IAPAR2(1) must be 1, 2, or 3'
208 c.....to use default initial conditions generated by the cascade,
209 c.....or to read in initial conditions.
212 c.....to generate formation time and the position at formation time from
213 c.....read-in initial conditions with an averaged formation proper time.
215 c.....ordering the particle label according to increasing order of
216 c.....formation time.
222 c-----------------------------------------------------------------------
224 c.....subroutine to generate formation time and position at formation time
225 c.....from read-in initial conditions with an averaged formation proper
230 c.....before invoking ARINI1:
231 c.....ARPAR1(1), IAINT2(1) must be set:
232 PARAMETER (MAXSTR=150001)
233 double precision smearp,smearh
235 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
237 COMMON /ARPRC/ ITYPAR(MAXSTR),
238 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
239 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
242 COMMON /smearz/smearp,smearh
244 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
246 common/anim/nevent,isoft,isflag,izpc
250 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
254 common /para8/ idpert,npertd,idxsec
256 clin-5/2008 for perturbatively-produced hadrons (currently only deuterons):
257 cms OPEN (91, FILE = 'ana/deuteron_processes.dat',
258 cms 1 STATUS = 'UNKNOWN')
259 if(idpert.eq.1.or.idpert.eq.2) then
260 cms OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN')
262 c.....generate formation time and position at formation time.
265 clin-7/10/01 initial positions already given for hadrons
266 c formed from partons inside ZPC (from string melting):
267 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
268 if(NP.le.nattzp) return
269 do 1001 I = nattzp+1, NP
270 IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
271 PRINT *, ' IN ARINI1'
272 PRINT *, 'ABS(PZ) .GE. EE for particle ', I
273 PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I),
275 PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
276 PRINT *, ' XM = ', XMAR(I)
280 RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I)))
282 VX = PXAR(I) / PEAR(I)
283 VY = PYAR(I) / PEAR(I)
284 FTAR(I) = TAU0 * COSH(RAP)
285 GXAR(I) = GXAR(I) + VX * FTAR(I)
286 GYAR(I) = GYAR(I) + VY * FTAR(I)
287 GZAR(I) = TAU0 * SINH(RAP)
288 clin-5/2009 No formation time for spectator projectile or target nucleons:
289 if(PXAR(I).eq.0.and.PYAR(I).eq.0
290 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99
291 2 .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
293 FTAR(I)=TAUI*COSH(RAP)
294 GZAR(I)=TAUI*SINH(RAP)
298 clin-3/2009 cleanup of program flow:
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),
306 PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
307 PRINT *, ' XM = ', XMAR(I)
312 RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I)))
314 VX = PXAR(I) / PEAR(I)
315 VY = PYAR(I) / PEAR(I)
316 c.....give initial formation time shift
317 TAUI = FTAR(I) + TAU0
318 FTAR(I) = TAUI * COSH(RAP)
319 GXAR(I) = GXAR(I) + VX * TAU0 * COSH(RAP)
320 GYAR(I) = GYAR(I) + VY * TAU0 * COSH(RAP)
321 c 4/25/03: hadron z-position upon formation determined the same way as x,y:
322 GZAR(I) = TAUI * SINH(RAP)
323 c the old prescription:
324 c GZAR(I) = GZAR(I) + TAU0 * SINH(RAP)
325 zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
326 GZAR(I)=GZAR(I)+zsmear
328 c 10/05/01 no formation time for spectator projectile or target nucleons:
329 if(PXAR(I).eq.0.and.PYAR(I).eq.0
330 1 .and.(PEAR(I)*2/HINT1(1)).gt.0.99
331 2 .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
335 FTAR(I)=TAUI*COSH(RAP)
336 GZAR(I)=TAUI*SINH(RAP)+zsmear
339 clin-3/2009 cleanup of program flow:
342 clin-3/2009 Add initial hadrons before the hadron cascade starts:
348 c-----------------------------------------------------------------------
350 c.....subroutine to order particle labels according to increasing
355 c.....before invoking ARTORD:
356 c.....IAINT2(1) must be set:
357 PARAMETER (MAXSTR=150001,MAXR=1)
358 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
360 COMMON /ARPRC/ ITYPAR(MAXSTR),
361 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
362 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
365 clin-3/2009 Take care of particle weights when user inserts initial hadrons:
366 COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
367 1 dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
368 2 dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
369 DIMENSION dptemp(MAXSTR)
371 DIMENSION ITYP0(MAXSTR),
372 & GX0(MAXSTR), GY0(MAXSTR), GZ0(MAXSTR), FT0(MAXSTR),
373 & PX0(MAXSTR), PY0(MAXSTR), PZ0(MAXSTR), EE0(MAXSTR),
375 DIMENSION INDX(MAXSTR)
393 dptemp(I) = dpertp(I)
395 CALL ARINDX(MAXSTR, NP, FT0, INDX)
398 c IF (ITYP0(INDX(I)) .EQ. 211) THEN
399 c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 321) THEN
400 c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
401 c & ITYP0(INDX(I)) .EQ. 2112 .OR. ITYP0(INDX(I)) .EQ. -211 .OR.
402 c & ITYP0(INDX(I)) .EQ. 111) THEN
403 c IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
404 c & ITYP0(INDX(I)) .EQ. 2112) THEN
406 c ITYPAR(I) = ITYP0(INDX(I))
407 c GXAR(I) = GX0(INDX(I))
408 c GYAR(I) = GY0(INDX(I))
409 c GZAR(I) = GZ0(INDX(I))
410 c FTAR(I) = FT0(INDX(I))
411 c PXAR(I) = PX0(INDX(I))
412 c PYAR(I) = PY0(INDX(I))
413 c PZAR(I) = PZ0(INDX(I))
414 c PEAR(I) = EE0(INDX(I))
415 c XMAR(I) = XM0(INDX(I))
416 ITYPAR(NPAR) = ITYP0(INDX(I))
417 GXAR(NPAR) = GX0(INDX(I))
418 GYAR(NPAR) = GY0(INDX(I))
419 GZAR(NPAR) = GZ0(INDX(I))
420 FTAR(NPAR) = FT0(INDX(I))
421 PXAR(NPAR) = PX0(INDX(I))
422 PYAR(NPAR) = PY0(INDX(I))
423 PZAR(NPAR) = PZ0(INDX(I))
424 PEAR(NPAR) = EE0(INDX(I))
425 XMAR(NPAR) = XM0(INDX(I))
427 dpertp(NPAR)=dptemp(INDX(I))
436 c-----------------------------------------------------------------------
438 c.....subroutine to copy individually generated particle record into
439 c.....particle record for many test particle runs.
443 PARAMETER (MAXSTR=150001,MAXR=1)
444 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
446 COMMON /ARPRC/ ITYPAR(MAXSTR),
447 & GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
448 & PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
451 COMMON /ARERC1/MULTI1(MAXR)
453 COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
454 & GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
456 & PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
457 & EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
459 COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
461 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
463 COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
464 & IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
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)
473 MULTI1(K) = IAINT2(1)
474 DO 1001 I = 1, MULTI1(K)
475 ITYP1(I, K) = ITYPAR(I)
485 clin-3/2009 hadron weights are initialized in addhad():
486 clin-5/2008 all hadrons not perturbatively-produced have the weight of 1:
491 c initialize final time of each particle to ntmax*dt except for
492 c decay daughters, which have values given by tfdcy() and >(ntmax*dt):
500 tfdpi(ip,irun)=NTMAX*DT
507 c=======================================================================
509 c.....function to convert PDG flavor code into ART flavor code.
511 FUNCTION IARFLV(IPDG)
513 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
520 IF (IPDG .EQ. -1114) THEN
526 IF (IPDG .EQ. -2114) THEN
532 IF (IPDG .EQ. -2214) THEN
538 IF (IPDG .EQ. -2224) THEN
545 IF (IPDG .EQ. -2212) THEN
551 IF (IPDG .EQ. -2112) THEN
558 IF (IPDG .EQ. 221) THEN
564 IF (IPDG .EQ. 2212) THEN
570 IF (IPDG .EQ. 2112) THEN
576 IF (IPDG .EQ. -211) THEN
582 IF (IPDG .EQ. 111) THEN
588 IF (IPDG .EQ. 211) THEN
594 IF (IPDG .EQ. 1114) THEN
600 IF (IPDG .EQ. 2114) THEN
606 IF (IPDG .EQ. 2214) THEN
612 IF (IPDG .EQ. 2224) THEN
618 IF (IPDG .EQ. 3122) THEN
624 IF (IPDG .EQ. -3122) THEN
630 IF (IPDG .EQ. 3112) THEN
636 IF (IPDG .EQ. -3112) THEN
642 IF (IPDG .EQ. 3212) THEN
648 IF (IPDG .EQ. -3212) THEN
654 IF (IPDG .EQ. 3222) THEN
660 IF (IPDG .EQ. -3222) THEN
666 IF (IPDG .EQ. -321) THEN
672 IF (IPDG .EQ. 321) THEN
677 c.....temporary entry for K0
678 IF (IPDG .EQ. 311) THEN
683 c.....temporary entry for K0bar
684 IF (IPDG .EQ. -311) THEN
689 c.....temporary entry for K0S and K0L
690 IF (IPDG .EQ. 310 .OR. IPDG .EQ. 130) THEN
701 IF (IPDG .EQ. -213) THEN
707 IF (IPDG .EQ. 113) THEN
713 IF (IPDG .EQ. 213) THEN
719 IF (IPDG .EQ. 223) THEN
725 IF (IPDG .EQ. 333) THEN
731 IF (IPDG .EQ. 323) THEN
736 IF (IPDG .EQ. -323) THEN
740 c.....temporary entry for K*0
741 IF (IPDG .EQ. 313) THEN
745 c.....temporary entry for K*0bar
746 IF (IPDG .EQ. -313) THEN
752 IF (IPDG .EQ. 331) THEN
758 c IF (IPDG .EQ. 777) THEN
764 IF (IPDG .EQ. 3312) THEN
770 IF (IPDG .EQ. -3312) THEN
776 IF (IPDG .EQ. 3322) THEN
782 IF (IPDG .EQ. -3322) THEN
788 IF (IPDG .EQ. 3334) THEN
794 IF (IPDG .EQ. -3334) THEN
800 IF (IPDG .EQ. 6666) THEN
806 clin-3/2009 keep the same ID numbers in case there are initial deuterons:
807 IF (IPDG .EQ. 42 .or. IPDG .EQ. -42) THEN
813 IARFLV = IPDG + 10000
818 c-----------------------------------------------------------------------
820 c.....function to convert ART flavor code into PDG flavor code.
822 FUNCTION INVFLV(IART)
824 common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
831 IF (IART .EQ. -6) THEN
837 IF (IART .EQ. -7) THEN
843 IF (IART .EQ. -8) THEN
849 IF (IART .EQ. -9) THEN
856 IF (IART .EQ. -1) THEN
862 IF (IART .EQ. -2) THEN
869 IF (IART .EQ. 0) THEN
875 IF (IART .EQ. 1) THEN
881 IF (IART .EQ. 2) THEN
887 IF (IART .EQ. 3) THEN
893 IF (IART .EQ. 4) THEN
899 IF (IART .EQ. 5) THEN
905 IF (IART .EQ. 6) THEN
911 IF (IART .EQ. 7) THEN
917 IF (IART .EQ. 8) THEN
923 IF (IART .EQ. 9) THEN
928 cc.....N*(1440), N*(1535) temporary entry
929 c IF (IART .GE. 10 .AND. IART .LE.13) THEN
935 IF (IART .EQ. 14) THEN
940 IF (IART .EQ. -14) THEN
946 c.....temporary entry for Sigma's
947 c IF (IART .EQ. 15) THEN
949 c IF (R .GT. 2. / 3.) THEN
951 c ELSE IF (R .GT. 1./ 3. .AND. R .LE. 2. / 3.) THEN
960 IF (IART .EQ. 15) THEN
966 IF (IART .EQ. -15) THEN
972 IF (IART .EQ. 16) THEN
978 IF (IART .EQ. -16) THEN
984 IF (IART .EQ. 17) THEN
990 IF (IART .EQ. -17) THEN
995 clin-2/23/03 K0S and K0L are generated at the last timestep:
996 c.....temporary entry for K- and K0bar
997 IF (IART .EQ. 21) THEN
999 c IF (R .GT. 0.5) THEN
1004 c IF (R .GT. 0.5) THEN
1013 c.....temporary entry for K+ and K0
1014 IF (IART .EQ. 23) THEN
1016 c IF (R .GT. 0.5) THEN
1021 c IF (R .GT. 0.5) THEN
1031 IF (IART .EQ. 22) THEN
1036 IF (IART .EQ. 24) THEN
1042 IF (IART .EQ. 25) THEN
1048 IF (IART .EQ. 26) THEN
1054 IF (IART .EQ. 27) THEN
1060 IF (IART .EQ. 28) THEN
1066 IF (IART .EQ. 29) THEN
1071 c.....temporary entry for K*+ and K*0
1072 IF (IART .EQ. 30) THEN
1074 IF (RANART(NSEED).GT.0.5) INVFLV = 313
1078 c.....temporary entry for K*- and K*0bar
1079 IF (IART .EQ. -30) THEN
1081 IF (RANART(NSEED).GT.0.5) INVFLV = -313
1085 c... eta-prime (bar)
1086 IF (IART .EQ. 31) THEN
1092 IF (IART .EQ. 32) THEN
1098 IF (IART .EQ. 40) THEN
1104 IF (IART .EQ. -40) THEN
1110 IF (IART .EQ. 41) THEN
1116 IF (IART .EQ. -41) THEN
1122 IF (IART .EQ. 45) THEN
1128 IF (IART .EQ. -45) THEN
1134 IF (IART .EQ. 44) THEN
1140 clin-5/2008 deuteron ID numbers in ART and ampt.dat:
1141 IF (IART .EQ. 42) THEN
1142 INVFLV = 1000000000+10020
1144 ELSEIF (IART .EQ. -42) THEN
1145 INVFLV = -1000000000-10020
1150 INVFLV = IART - 10000
1155 c=======================================================================
1159 COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
1162 DATA ARPAR1/1.19, 99 * 0.0/
1163 DATA IAPAR2/3, 49 * 0/
1164 DATA ARINT1/100 * 0.0/
1169 c=======================================================================
1171 c.....Routine borrowed from ZPC.
1172 c.....double precision is modified to real*4.
1175 c subroutine index1(n, m, arrin, indx)
1176 subroutine arindx(n, m, arrin, indx)
1178 c indexes the first m elements of ARRIN of length n, i.e., outputs INDX
1179 c such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m
1181 c implicit real*4 (a-h, o-z)
1183 dimension arrin(n), indx(n)
1207 20 if (j .le. ir) then
1209 if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1
1211 if (q .lt. arrin(indx(j))) then
1225 c-----------------------------------------------------------------------
1227 c.....extracted from G. Song's ART expasion including K- interactions
1228 c.....file `NEWKAON.FOR'
1230 c 5/01/03 send iblock value into art1f.f, necessary for resonance studies:
1231 c subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1232 c & srt,pcx,pcy,pcz)
1233 subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1234 & srt,pcx,pcy,pcz,iblock)
1235 PARAMETER (MAXSTR=150001,MAXR=1)
1236 PARAMETER (AKA=0.498)
1237 COMMON /AA/ R(3,MAXSTR)
1239 COMMON /BB/ P(3,MAXSTR)
1241 COMMON /CC/ E(MAXSTR)
1243 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1245 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
1251 COMMON /PA/RPION(3,MAXSTR,MAXR)
1253 COMMON /PB/PPION(3,MAXSTR,MAXR)
1255 COMMON /PC/EPION(MAXSTR,MAXR)
1257 COMMON /PD/LPION(MAXSTR,MAXR)
1263 logical lb1bn, lb2bn,lb1mn,lb2mn
1265 c logical lb1bn1, lb2bayon1, lb1bn0, lb2bn0
1266 logical lb1bn1, lb2bn1, lb1bn0, lb2bn0
1267 cbz3/7/99 neutralk end
1268 logical lb1mn0, lb2mn0, lb1mn1, lb2mn1
1269 logical lb1mn2, lb2mn2
1271 c icase: flag for the type of reaction that is going to happen.
1272 c icase=-1, no desired reaction, return to main program.
1275 c 3, K(-) absorption.
1277 c nchrg: Net charges of the two incoming particles.
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
1300 c 1. consider N+N, N+Resonance, R + R reactions
1301 if(lb1bn.and.lb2bn) then
1304 c total cross section
1306 if(lb1.eq.9.and.lb2.eq.9) then
1309 if((lb1bn1.and.lb2.eq.9)
1310 & .or.(lb2bn1.and.lb1.eq.9))then
1313 if((lb1bn0.and.lb2.eq.9)
1314 & .or.(lb2bn0.and.lb1.eq.9)
1315 & .or.(lb1bn1.and.lb2bn1)) then
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
1323 if((lb1bn0.and.lb2bn0).or.(lb1bn1.and.lb2.eq.6)
1324 & .or.(lb2bn1.and.lb1.eq.6)) then
1327 if((lb1bn0.and.lb2.eq.6)
1328 & .or.(lb2bn0.and.lb1.eq.6))then
1331 if(lb1.eq.6.and.lb2.eq.6) then
1334 c brsig = x2kaon_no_isospin(srt)
1335 if(nchrg.ge.-1.and.nchrg.le.2) then
1336 c K,Kbar prduction x sect.
1340 c if(nchrg.eq.-2.or.nchrg.eq.3) then
1341 c brsig = x2kaon(srt+0.938-1.232)
1344 c brsig = x2kaon(srt+2.*(0.938-1.232))
1350 cbz3/7/99 neutralk end
1354 c 2. consider PI(meson:eta,omega,rho,phi) + N(N*,D)
1355 if((lb1bn.and.lb2mn).OR.(lb2bn.and.lb1mn)) then
1359 sigma0 = piNsg0(srt)
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
1367 c if(lb1bn1.or.lb2bn1) brsig=2.0*sigma0
1368 c if(lb1bn0.or.lb2bn0) brsig=0.5*sigma0
1369 if(lb1bn1.or.lb2bn1) brsig=0.5*sigma0
1370 if(lb1bn0.or.lb2bn0) brsig=2.0*sigma0
1372 c if(lb1.eq.9.or.lb2.eq.9) brsig=1.5*sigma0
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
1379 if(lb1bn1.or.lb2bn1) then
1389 c ratiok: the ratio of channels: ->nK+k- vs. -> pK0K-
1391 if(lb1bn0.or.lb2bn0) then
1398 c if(lb1.eq.6.or.lb2.eq.6) then
1404 if( (lb1bn0.and.lb2mn2)
1405 & .or.(lb2bn0.and.lb1mn2)
1406 & .or.(lb1.eq.6.and.lb2mn0).or.(lb2.eq.6.and.lb1mn0)) then
1408 if(lb1bn0.or.lb2bn0) brsig=sigma0
1409 c if(lb1.eq.6.or.lb2.eq.6) brsig=sigma0
1411 c if((lb1.eq.6.and.lb2mn2).or.(lb2.eq.6.and.lb1mn2))then
1414 c if((lb1bn1.and.lb2mn1).or.(lb2bn1.and.lb1mn1)
1415 c & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1420 if((lb1.eq.6.and.lb2mn2)
1421 & .or.(lb2.eq.6.and.lb1mn2))then
1426 if((lb1bn1.and.lb2mn1)
1427 & .or.(lb2bn1.and.lb1mn1)
1428 & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1431 cbz3/8/99 neutralk end
1434 IF (NCHRG .GE. -2 .AND. NCHRG .LE. 2) THEN
1435 BRSIG = 3.0 * SIGMA0
1437 cbz3/7/99 neutralk end
1441 c 3. consider K- + N(N*,D) absorption.
1442 c if((lb1bn.and.lb2.eq.21).OR.(lb2bn.and.lb1.eq.21)) then
1443 if( (lb1bn.and.(lb2.eq.21.or.lb2.eq.-30)).OR.
1444 & (lb2bn.and.(lb1.eq.21.or.lb1.eq.-30)) )then
1447 if(srt.le.(bmass+aka)) then
1449 c write(100,*)'--lb1,lb2,em1,em2,srt',lb1,lb2,em1,em2,srt
1453 pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
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
1461 sigsgm=3.*akPsgm(pkaon)
1462 sig=sigela+sigsgm+akPlam(pkaon)
1464 if(lb1.eq.2.or.lb2.eq.2.or.lb1.eq.7.or.lb2.eq.7.or.
1465 & lb1.eq.10.or.lb2.eq.10.or.lb1.eq.12.or.lb2.eq.12) then
1466 c K- + (D0, N*0)n ->
1469 sigsgm=2.*akNsgm(pkaon)
1470 sig=sigela+sigsgm+akNlam(pkaon)
1472 if(lb1.eq.6.or.lb2.eq.6) then
1476 sigsgm=akNsgm(pkaon)
1479 if(lb1.eq.9.or.lb2.eq.9) then
1483 sigsgm=2.*akPsgm(pkaon)
1484 sig=sigela+sigsgm+akPlam(pkaon)
1488 sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON))
1489 SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1490 SIG = sigela + SIGSGM + AKPLAM(PKAON)
1491 cbz3/8/99 neutralk end
1493 if(sig.gt.1.e-7) then
1494 c K(-) + N reactions
1498 c branch_lambda=akNlam(pkaon)/sig
1503 c 4. meson + hyperon -> K- + N
1504 c if(((lb1.ge.14.and.lb1.le.17).and.lb2mn).OR.
1505 c & ((lb2.ge.14.and.lb2.le.17).and.lb1mn)) then
1506 if(((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.3.and.lb2.le.5)).OR.
1507 & ((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.3.and.lb1.le.5)))then
1508 c first classify the reactions due to total charge.
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
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
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
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
1541 if(nchrg.ne.-100.and.srt.gt.(aka+bmass)) then
1542 c PI+sigma or PI + Lambda => Kbar + N reactions
1544 c pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
1545 pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2)
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)
1553 if(nchrg.ge.0) sigma0=akPsgm(pkaon)
1555 if(nchrg.lt.0) sigma0=akNsgm(pkaon)
1558 SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1559 cbz3/8/99 neutralk end
1562 sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/
1563 & (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0
1565 c if(nchrg.eq.-2.or.nchrg.eq.1) sig=2.*sig K-D++, K-D-
1567 if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig
1569 cbz3/8/99 neutralk end
1571 c the factor 2 comes from spin of delta, which is 3/2
1572 c detailed balance. copy from Page 423 of N.P. A614 1997
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
1580 SIG = 4.0 / 9.0 * SIG
1582 cbz3/8/99 neutralk end
1584 if(sig.lt.1.e-7) sig = 1.e-7
1587 * comment icase=4 statement below if only inelastic
1588 c PI+L/Si => Kbar + N OR ELASTIC SCATTERING
1591 c elastic xsecn of 10mb
1599 c if(em2.lt.0.2.and.em1.lt.0.2) then
1602 c assumed PI PI total x section.
1607 c if(srt.gt.s0) brsig = 2.7*(1.-s0**2/srt**2)**0.76
1608 c x section for PIPI->KKbar PRC43 (1991) 1881
1610 if(icase.eq.-1) then
1619 ec=(em1+em2+0.02)**2
1621 c if((e(i1).ge.1.).and.(e(i2).ge.1.)) ec = 4.75
1623 call distce(i1,i2,dsr,ds,dt,ec,srt,ic,px1cm,py1cm,pz1cm)
1625 c no anti-kaon production
1628 c write(60,*)'--------------distance-----',in
1632 clin-10/24/02 set to 0: ik,ik0-3,il,im,im3-4,in,inpion,ipipi,
1633 c sgsum,sgsum1,sgsum3:
1651 if(srt.gt.2.8639) then
1653 if(em1.lt.1.0.and.em2.lt.1.0) then
1656 c ratio_1=sgsum1/ik1/40.
1658 if(em1.gt.1.0.and.em2.gt.1.0) then
1661 c ratio_3=sgsum3/ik3/40.
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
1666 c ratio=sgsum/ik0/40.
1669 if(icase.eq.2) inpion=inpion+1
1670 if(icase.eq.5) ipipi=ipipi+1
1671 c write(62,*)'ik1,ik2,ik3',ik1,ik2,ik3,ratio_1,ratio_3,ratio
1672 c write(62,*)'inpion,ipipi',inpion,ipipi
1673 if(RANART(NSEED).gt.(brsig/sig)) then
1674 c no anti-kaon production
1679 c kaons could be created now.
1682 c write(60,*)'------in,s2kaon,sig=',in,brsig,sig,lb1,lb2
1683 call nnkaon(irun,iseed,
1684 & ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
1688 c call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1689 c & pcx,pcy,pcz,nchrg,ratiok)
1690 call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1691 & pcx,pcy,pcz,nchrg,ratiok,iblock)
1696 c write(63,*)'im3,lb1,lb2,pkaon',im3,lb1,lb2,pkaon
1697 c write(63,*)'sig,el,sigma',sig,brel,brsgm
1698 c write(63,*)'srt,pcx,pcy,pcz,em1,em2',srt,pcx,pcy,pcz,em1,em2
1699 call kaonN(brel,brsgm,irun,iseed,dt,nt,ictrl,
1700 & i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
1701 c this subroutine format is diff. since three final states are possible
1707 c write(64,*)'im4,sigma0,branch,sig=',im4,sigma0,brsig,sig
1708 c write(64,*)'lb1,lb2,em1,em2,pkaon=',lb1,lb2,em1,em2,pkaon
1711 if(RANART(NSEED).lt.brel) then
1716 c call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1717 c & pcx,pcy,pcz,nchrg)
1718 call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1719 & pcx,pcy,pcz,nchrg,iblock)
1723 c if(icase.eq.5) then
1725 c write(65,*)'---im5,s2kaon,sig=',im5,brsig,sig
1726 c call pipikaon(irun,iseed,dt,nt,ictrl,i1,i2,srt,pcx,pcy,pcz)
1729 c write(101,*)lb1,lb2,lb(i1),lb(i2)
1730 c write(101,*)em1,em2,e(i1),e(i2),srt
1736 ******************************************
1737 * for pp-->pp + kaon + anti-kaon
1738 c real*4 function X2kaon(srt)
1739 real function X2kaon(srt)
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 *
1745 ******************************************
1746 c minimum c.m.s. energy to create 2 kaon. = 2*(mp+mk)
1749 if(srt.lt.smin)return
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
1761 real function piNsg0(srt)
1763 * cross section in mb for PI- + P -> P + K0 + K-
1765 srt0 = 0.938 + 2.*0.498
1766 if(srt.lt.srt0) then
1770 ratio = srt0**2/srt**2
1771 piNsg0=1.121*(1.-ratio)**1.86*ratio**2
1775 real function akNel(pkaon)
1777 *cross section in mb for K- + N reactions.
1778 c the following data come from PRC 41 (1701)
1779 c sigma1: K(-) + neutron elastic
1780 if(pkaon.lt.0.5.or. pkaon.ge.4.0) sigma1=0.
1781 if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma1=20.*pkaon**2.74
1782 if(pkaon.ge.1.0.and.pkaon.lt.4.0) sigma1=20.*pkaon**(-1.8)
1787 real function akPel(pkaon)
1789 *cross section in mb for K- + N reactions.
1790 c the following data come from PRC 41 (1701)
1791 c sigma2: K(-) + proton elastic
1792 if(pkaon.lt.0.25.or. pkaon.ge.4.0) sigma2=0.
1793 if(pkaon.ge.0.25.and.pkaon.lt.4.0) sigma2=13.*pkaon**(-0.9)
1798 real function akNsgm(pkaon)
1800 *cross section in mb for K- + N reactions.
1801 c sigma2: x section for K- + n -> sigma0 + PI-
1802 if(pkaon.lt.0.5.or. pkaon.ge.6.0) sigma2=0.
1803 if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma2=1.2*pkaon**(-1.3)
1804 if(pkaon.ge.1.0.and.pkaon.lt.6.0) sigma2=1.2*pkaon**(-2.3)
1809 real function akPsgm(pkaon)
1811 *cross section in mb for K- + N reactions.
1812 c sigma1: x section for K- + p -> sigma0 + PI0
1813 if(pkaon.lt.0.2.or. pkaon.ge.1.5) sigma1=0.
1814 if(pkaon.ge.0.2.and.pkaon.lt.1.5) sigma1=0.6*pkaon**(-1.8)
1819 real function akPlam(pkaon)
1821 *cross section in mb for K- + N reactions.
1822 c sigma: x section for K- + p -> lambda + PI0
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)
1831 real function akNlam(pkaon)
1833 *cross section in mb for K- + N reactions.
1834 akNlam=akPlam(pkaon)
1838 * GQ Li parametrization (without resonance)
1839 real function akNPsg(pkaon)
1841 *cross section in mb for K- + N reactions.
1842 c sigma1: x section for K- + p/n -> sigma0 + PI0
1843 if(pkaon.le.0.345)then
1844 sigma1=0.624*pkaon**(-1.83)
1846 sigma1=0.7*pkaon**(-2.09)
1852 c-----------------------------------------------------------------------
1854 c.....extracted from G. Song's ART expasion including K- interactions
1855 c.....file `NEWNNK.FOR'
1857 subroutine nnkaon(irun,iseed,ictrl,i1,i2,iblock,
1858 & srt,pcx,pcy,pcz,nchrg)
1859 c <pt>=0.27+0.037*log(srt) was changed to 0.632 + ... on Aug. 14, 1997
1860 c CANCELED also alpha=1 changed to alpha=3 to decrease the leadng effect.
1861 PARAMETER (MAXSTR=150001,MAXR=1)
1862 PARAMETER (AKA=0.498)
1863 COMMON /AA/ R(3,MAXSTR)
1865 COMMON /BB/ P(3,MAXSTR)
1867 COMMON /CC/ E(MAXSTR)
1869 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1871 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
1877 COMMON /PA/RPION(3,MAXSTR,MAXR)
1879 COMMON /PB/PPION(3,MAXSTR,MAXR)
1881 COMMON /PC/EPION(MAXSTR,MAXR)
1883 COMMON /PD/LPION(MAXSTR,MAXR)
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)
1894 c 10/24/02 initialize n to 0:
1898 c if(nchrg.eq.-2.or.nchrg.ge.3) dm3=1.232
1899 c if(nchrg.eq.4) dm4=1.232
1900 if(nchrg.le.-1.or.nchrg.ge.3) dm3=1.232
1901 if(nchrg.eq.-2.or.nchrg.eq.4) dm4=1.232
1902 cbz3/11/99 neutralk end
1904 call fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
1906 c write(60,*)'------------final state fail-------',n
1907 c no anti-kaon production
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
1920 c 10/28/02 get rid of argument usage mismatch in rotate():
1924 c call rotate(pcx,pcy,pcz,px(1),py(1),pz(1))
1925 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1933 c call rotate(pcx,pcy,pcz,px(2),py(2),pz(2))
1934 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1942 c call rotate(pcx,pcy,pcz,px(3),py(3),pz(3))
1943 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1951 c call rotate(pcx,pcy,pcz,px(4),py(4),pz(4))
1952 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1960 if(nchrg.eq.-1.or.nchrg.eq.-2) then
1961 c To keep charge conservation. D-n->nnK0K-, D-D- -> nD-K0K-
1964 c lpion(nnn,irun)=24 ! K0
1965 cbz3/7/99 neutralk end
1968 c aka: rest mass of K
1971 lpion(nnn-1,irun)=21
1972 c aka: rest mass of K
1973 epion(nnn-1,irun)=aka
1974 * Find the momenta of particles in the final state in the nucleus_nucleus
1975 * cms frame. Lorentz transformation into lab frame.
1976 e1cm = sqrt(dm3**2 + px(1)**2 + py(1)**2 + pz(1)**2)
1977 p1beta = px(1)*betax + py(1)*betay + pz(1)*betaz
1978 transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
1979 pt1i1 = betax*transf + px(1)
1980 pt2i1 = betay*transf + py(1)
1981 pt3i1 = betaz*transf + pz(1)
1985 if(nchrg.ge.-2.and.nchrg.le.1) lb1=2
1988 if (nchrg .eq. -2 .or. nchrg .eq. -1) then
1991 cbz3/7/99 neutralk end
1994 c if(nchrg.eq.2.or.nchrg.eq.3) lb1=1
1995 c if(nchrg.eq.4) lb1=9
1996 if(nchrg.eq.1.or.nchrg.eq.2) lb1=1
1997 if(nchrg.eq.3.or.nchrg.eq.4) lb1=9
1998 cbz3/11/99 neutralk end
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)
2012 c if(nchrg.eq.-1.or.nchrg.eq.0) lb2=2
2013 c if(nchrg.eq. 2.or.nchrg.eq.1) lb2=1
2014 c if(nchrg.eq. 4.or.nchrg.eq.3) lb2=9
2015 c if(nchrg.eq.-2) lb2=6
2016 if(nchrg.ge.-1.or.nchrg.le.1) lb2=2
2017 if(nchrg.eq. 2.or.nchrg.eq.3) lb2=1
2018 if(nchrg.eq. 4) lb2=9
2019 if(nchrg.eq.-2) lb2=6
2020 cbz3/11/99 neutralk end
2022 c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2041 c iblock = 101 ! K(+)K(-) production
2042 * Get anti-kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2043 epcmk = sqrt(epion(nnn-1,irun)**2 + px(3)**2+py(3)**2+pz(3)**2)
2044 betak = px(3)*betax + py(3)*betay + pz(3)*betaz
2045 transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
2046 ppion(1,nnn-1,irun)=betax*transf + px(3)
2047 ppion(2,nnn-1,irun)=betay*transf + py(3)
2048 ppion(3,nnn-1,irun)=betaz*transf + pz(3)
2049 rpion(1,nnn-1,irun)=r(1,i1)
2050 rpion(2,nnn-1,irun)=r(2,i1)
2051 rpion(3,nnn-1,irun)=r(3,i1)
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)
2065 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2069 subroutine lorntz(ilo,b,pi,pj)
2070 c It uses to perform Lorentz (or inverse Lorentz) transformation
2071 dimension pi(4),pj(4),b(3)
2074 bb=b(1)*b(1)+b(2)*b(2)+b(3)*b(3)
2076 if(deno3.eq.0.)deno3=1.e-10
2079 if(ilo.eq.1) goto 100
2080 c Lorentz transformation
2081 pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
2082 pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
2083 c drb=drd(1)*b(1)+drd(2)*b(2)+drd(3)*b(3)
2084 c drdb=db(1)*b(1)+db(2)*b(2)+db(3)*b(3)
2086 pi(i)=pi(i)+b(i)*(ga*pib-gam*pi(4))
2087 pj(i)=pj(i)+b(i)*(ga*pjb-gam*pj(4))
2088 c drd(i)=drd(i)+b(i)*ga*drb
2089 c db(i)=db(i)+b(i)*ga*drdb
2091 pi(4)=gam*(pi(4)-pib)
2092 pj(4)=gam*(pj(4)-pjb)
2095 c inverse Lorentz transformation
2096 pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
2097 pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
2099 pi(i)=pi(i)+b(i)*(ga*pib+gam*pi(4))
2100 pj(i)=pj(i)+b(i)*(ga*pjb+gam*pj(4))
2102 pi(4)=gam*(pi(4)+pib)
2103 pj(4)=gam*(pj(4)+pjb)
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)
2116 c iflag=-1: fail to find momenta
2129 ekmax=(srt-dm3-dm4)/2.
2130 if(ekmax.le.aka) return
2131 pkmax=sqrt(ekmax**2-aka**2)
2133 if(dm3.le.0.0.or.dm4.le.0.0) then
2134 write(1,*)'error: minus mass!!!'
2138 c after we have the momenta for both nucleus, we sample the
2139 c transverse momentum for K-.
2140 c dsigma/dpt**2 = exp(-4.145*pt**2) obtained by fitting data on
2144 if(icount.gt.10) return
2145 ptkmi2=-1./4.145*alog(RANART(NSEED))
2150 if(rsq.ge.1.0.or.rsq.le.0.) goto 3
2151 fac=sqrt(-2.*alog(rsq)/rsq)
2153 if(guass.ge.5.) goto 3
2156 ekm=sqrt(aka**2+pzkm**2+ptkm**2)
2157 if(RANART(NSEED).gt.aka/ekm) goto 50
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
2164 150 ptkpl2=-1./3.68*alog(RANART(NSEED))
2169 if(rsq.ge.1.0.or.rsq.le.0.) goto 13
2170 fac=sqrt(-2.*alog(rsq)/rsq)
2172 if(guass.ge.3.25) goto 13
2175 ekp=sqrt(aka**2+pzkp**2+ptkp**2)
2176 if(RANART(NSEED).gt.aka/ekp) goto 150
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
2184 resten=srt-pe(3)-pe(4)
2187 if(resten.le.abs(restpz)) goto 50
2188 restms=sqrt(resten**2-restpz**2)
2190 if(restms.lt.(dm3+dm4)) goto 50
2191 ptp2=-1./2.76*alog(RANART(NSEED))
2194 px(2)=ptp*cos(2.*pio*bbb)
2195 py(2)=ptp*sin(2.*pio*bbb)
2196 px(1)=-1.*(px(4)+px(3)+px(2))
2197 py(1)=-1.*(py(4)+py(3)+py(2))
2198 c transverse mass for K-
2199 rmt3=sqrt(dm3**2+px(1)**2+py(1)**2)
2200 c transverse mass for K+
2201 rmt4=sqrt(dm4**2+px(2)**2+py(2)**2)
2202 if(restms.lt.(rmt3+rmt4)) goto 50
2203 c else: sampling success!
2204 pzcms=sqrt((restms**2-(rmt3+rmt4)**2)*
2205 & (restms**2-(rmt3-rmt4)**2))/2./restms
2206 if(RANART(NSEED).gt.0.5) then
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)
2224 c-----------------------------------------------------------------------
2226 c.....extracted from G. Song's ART expasion including K- interactions
2227 c.....file `NPIK.FOR'
2229 ****************************************
2230 c subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2231 c & pcx,pcy,pcz,nchrg,ratiok)
2232 subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2233 & pcx,pcy,pcz,nchrg,ratiok,iblock)
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)
2249 COMMON /BB/ P(3,MAXSTR)
2251 COMMON /CC/ E(MAXSTR)
2253 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2255 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
2261 COMMON /PA/RPION(3,MAXSTR,MAXR)
2263 COMMON /PB/PPION(3,MAXSTR,MAXR)
2265 COMMON /PC/EPION(MAXSTR,MAXR)
2267 COMMON /PD/LPION(MAXSTR,MAXR)
2269 dimension bb(3),p1(4),p2(4),p3(4),px(4),py(4),pz(4)
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)
2292 c k1 must be bayron. k2 be meson. If not, exchange.
2293 if(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13)) then
2299 c LB(I1) = 1 + 2 * RANART(NSEED)
2301 LB(k1) = 1 + int(2*RANART(NSEED))
2303 c pkmax=sqrt((srt**2-(aka+0.938+aka)**2)*(srt**2-(aka+0.938-aka)**2))
2305 pkmax=sqrt((srt**2-(aka+0.938+aka)**2)
2306 & *(srt**2-(aka+0.938-aka)**2))/2./srt
2307 pk = RANART(NSEED)*pkmax
2308 c-----------------------------------------------------
2309 css=1.-2.*RANART(NSEED)
2311 fai=2*3.1415926*RANART(NSEED)
2312 p3(1)=pk*sss*cos(fai)
2313 p3(2)=pk*sss*sin(fai)
2315 eip = srt - sqrt(aka**2 + pk**2)
2316 rmnp=sqrt(eip**2-pk**2)
2318 bb(i) = -1.*p3(i)/eip
2320 c bb: velocity of the other two particles as a whole.
2321 pznp=sqrt((rmnp**2-(aka+0.938)**2)
2322 c *(rmnp**2-(0.938-aka)**2))/2./rmnp
2323 c-----------------------------------------------------
2324 css=1.-2.*RANART(NSEED)
2326 fai=2*3.1415926*RANART(NSEED)
2327 p1(1)=pznp*sss*cos(fai)
2328 p1(2)=pznp*sss*sin(fai)
2330 p1(4)=sqrt(0.938**2+pznp**2)
2331 p2(4)=sqrt(aka**2+pznp**2)
2335 c p1,p2: the momenta of the two particles in their cms
2336 c p1: momentum of N or P
2337 c p2: momentum of anti_kaon
2338 c p3: momentum of K0 or K+
2340 c write(61,*)'--------p1,p2',p1,p2
2341 c write(61,*)'--------bb',bb
2342 call lorntz(ilo,bb,p1,p2)
2343 c******* Checking *************
2344 c pxsum = p1(1)+p2(1)+p3(1)
2345 c pysum = p1(2)+p2(2)+p3(2)
2346 c pzsum = p1(3)+p2(3)+p3(3)
2347 c pesum = p1(4)+p2(4)+sqrt(p3(1)**2+p3(2)**2+p3(3)**2+aka**2)-srt
2348 c write(61,*)'---p1,pxsum',p1,pxsum
2349 c write(61,*)'---p2,pysum',p2,pysum
2350 c write(61,*)'---p3,pzsum',p3,pzsum
2351 c write(61,*)'---pesum',pesum
2352 c***********************************
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
2359 c 10/28/02 get rid of argument usage mismatch in rotate():
2363 c call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2364 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2372 c call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2373 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2381 c call rotate(pcx,pcy,pcz,p3(1),p3(2),p3(3))
2382 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2390 c aka: rest mass of K
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)
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)
2413 c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2414 * k1 stand for nucleon, k2 stand for kaon. lpion stand for Kbar.
2433 c K(+)K(-) production
2435 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2436 c p2: momentum of anti-kaon.
2437 c epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2438 epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2+p2(2)**2+p2(3)**2)
2439 betak = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2440 transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
2441 ppion(1,nnn,irun)=betax*transf + p2(1)
2442 ppion(2,nnn,irun)=betay*transf + p2(2)
2443 ppion(3,nnn,irun)=betaz*transf + p2(3)
2445 dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2447 c write(400,*)'2 ', ppion(1,nnn,irun), ppion(2,nnn,irun),
2448 c & ppion(3,nnn,irun), dt*nt, srt
2450 c write(420,*)ppion(1,nnn,irun), ppion(2,nnn,irun),
2451 c & ppion(3,nnn,irun), dt*nt, srt
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)
2460 c-----------------------------------------------------------------------
2462 c.....extracted from G. Song's ART expasion including K- interactions
2463 c.....file `PIHYPN.FOR'
2465 ******************************************
2466 subroutine pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,
2467 & srt,pcx,pcy,pcz,nchrg,iblock)
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 ******************************************
2473 c NOTE: for PI + Hyperon: the produced kaons have mass 0.498
2474 PARAMETER (MAXSTR=150001,MAXR=1,PI=3.1415926)
2475 PARAMETER (AKA=0.498)
2476 COMMON /AA/ R(3,MAXSTR)
2478 COMMON /BB/ P(3,MAXSTR)
2480 COMMON /CC/ E(MAXSTR)
2482 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2484 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
2490 COMMON /PA/RPION(3,MAXSTR,MAXR)
2492 COMMON /PB/PPION(3,MAXSTR,MAXR)
2494 COMMON /PC/EPION(MAXSTR,MAXR)
2496 COMMON /PD/LPION(MAXSTR,MAXR)
2498 dimension p1(4),p2(4)
2511 if(ielstc .eq. 1) then
2512 * L/Si + meson -> L/Si + meson
2521 c PI + Sigma(or Lambda) -> Kbar + N
2524 c k1 must be bayron! So if I1 is PI, exchange k1 & k2.
2525 if(lb(i1).lt.14.or.lb(i1).gt.17) then
2530 LB(K1) = 1 + int(2*RANART(NSEED))
2531 if(nchrg.eq.-2) lb(k1)=6
2532 c if(nchrg.eq.-1) lb(k1)=2
2533 c if(nchrg.eq. 0) lb(k1)=1
2534 c if(nchrg.eq. 1) lb(k1)=9
2535 IF (NCHRG .EQ. 2) LB(K1) = 9
2536 cbz3/8/99 neutralk end
2541 if(nchrg.eq.-2.or.nchrg.eq.1) dm3=1.232
2543 c dm3,dm4: the mass of final state particles.
2546 ********Now, antikaon will be created.
2547 c call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2548 c pkmax: the maximum momentum of anti-kaon
2549 pkmax=sqrt((srt**2-(dm3+dm4)**2)*(srt**2-(dm3-dm4)**2))
2552 c-----------------------------------------------------
2553 css=1.-2.*RANART(NSEED)
2555 fai=2*3.1415926*RANART(NSEED)
2556 p1(1)=pk*sss*cos(fai)
2557 p1(2)=pk*sss*sin(fai)
2562 c p1,p2: the momenta of the two particles in their cms
2563 c p1: momentum of kaon
2564 c p2: momentum of Kbar
2566 * Rotate the momenta of particles in the cms of I1 & I2
2567 clin-10/28/02 get rid of argument usage mismatch in rotate():
2571 c call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2572 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2580 c call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2581 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
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)
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)
2606 c write(400,*)'3 ', pt1i2, pt2i2, pt3i2, dt*nt, srt
2608 c write(430,*)pt1i2, pt2i2, pt3i2, dt*nt, srt
2612 c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2615 * k1 stand for nucleon, k2 stand for kaon.
2627 cc iblock = 101 ! K(+)K(-) production
2628 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2632 c-----------------------------------------------------------------------
2634 c.....extracted from G. Song's ART expasion including K- interactions
2635 c.....file `KAONN.FOR'
2637 ****************************************
2638 subroutine kaonN(brel,brsgm,irun,iseed,dt,nt,
2639 & ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
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)
2648 COMMON /BB/ P(3,MAXSTR)
2650 COMMON /CC/ E(MAXSTR)
2652 COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
2654 COMMON /BG/BETAX,BETAY,BETAZ,GAMMA
2660 COMMON /PA/RPION(3,MAXSTR,MAXR)
2662 COMMON /PB/PPION(3,MAXSTR,MAXR)
2664 COMMON /PC/EPION(MAXSTR,MAXR)
2666 COMMON /PD/LPION(MAXSTR,MAXR)
2668 dimension p1(4),p2(4)
2681 c ratio: used for isospin decision.
2684 c k1 must be bayron! So if I1 is Kaon, exchange k1 & k2.
2685 if(e(i1).lt.0.5.and.e(i1).gt.0.01) then
2689 *** note: for print out only *******************************
2690 c record kaon's mass
2692 *** end **************
2694 if(rrr.lt.brel) then
2695 c Kbar + N -> Kbar + N
2703 if(rrr.lt.(brel+brsgm)) then
2704 c nchrg: Net charges of the two incoming particles.
2705 c Kbar + N -> Sigma + PI
2710 LB1 = 15 + int(3*RANART(NSEED))
2711 LB2 = 3 + int(3*RANART(NSEED))
2713 c Kbar + N -> Lambda + PI
2719 LB2 = 3 + int(3*RANART(NSEED))
2720 c if(nchrg.eq.1) lb2=5 ! K- + D++ -> Lambda + PI+
2721 c if(nchrg.eq.0) lb2=4 ! K- + p(D+,N*+) -> Lambda + PI0
2722 c if(nchrg.eq.-1) lb2=3 ! K- + n(D,N*) -> Lambda + PI-
2730 ********Now, antikaon will be created.
2731 c call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2732 c pkmax: the maximum momentum of anti-kaon
2733 c write(63,*)'srt,em1,em2',srt,em1,em2
2734 c write(63,*)'-srt,em1,em2',srt,em1,em2
2735 pkmax=sqrt((srt**2-(em1+em2)**2)*(srt**2-(em1-em2)**2))
2738 c-----------------------------------------------------
2739 css=1.-2.*RANART(NSEED)
2741 fai=2*3.1415926*RANART(NSEED)
2742 p1(1)=pk*sss*cos(fai)
2743 p1(2)=pk*sss*sin(fai)
2748 c p1,p2: the momenta of the two particles in their cms
2749 c p1: momentum of kaon
2750 c p2: momentum of Kbar
2752 * Rotate the momenta of particles in the cms of I1 & I2
2754 clin-10/28/02 get rid of argument usage mismatch in rotate():
2758 c call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2759 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2767 c call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2768 call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
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)
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)
2793 c if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2796 * k1 stand for bayron, k2 stand for meson.
2806 cc iblock = 101 ! K(+)K(-) production
2807 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2811 c=======================================================================
2813 clin Below is the previous artana.f:
2814 c=======================================================================
2816 c.....analysis subroutine before the hadronic space-time evolution
2819 PARAMETER (MAXSTR=150001, MAXR=1)
2820 c.....y cut for mt spectrum
2822 c PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
2823 PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
2825 c.....bin width for mt spectrum and y spectrum
2826 clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
2827 c PARAMETER (BMT = 0.05, BY = 0.2)
2828 PARAMETER (BMT = 0.05, BY = 0.4)
2831 COMMON /ARERC1/MULTI1(MAXR)
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)
2839 c & dm1k0s(50), DMT1LA(50), DMT1LB(50)
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)
2856 DO 1001 I = 1, MULTI1(J)
2863 c 2/24/03 leptons and photons:
2864 if(xm.lt.0.01) goto 200
2865 ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
2866 eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
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'
2873 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
2876 PRINT *, ' PZ = ', PZ, ' EE = ', EE
2878 PRINT *, ' XM = ', XM
2885 Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
2886 c.....rapidity cut for the rapidity distribution
2887 IF (ABS(Y) .GE. 10.0) GOTO 100
2888 clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
2889 c IY = 1 + int(ABS(Y) / BY)
2890 c Ieta = 1 + int(ABS(eta) / BY)
2891 IF (ABS(eta) .GE. 10.0) GOTO 100
2892 IY = 1 + int((Y+10.) / BY)
2893 Ieta = 1 + int((eta+10.) / BY)
2895 IF (ITYP .LT. -1000) THEN
2896 dy1ntb(IY) = dy1ntb(IY) - 1.0
2898 IF (ITYP .GT. 1000) THEN
2899 dy1ntb(IY) = dy1ntb(IY) + 1.0
2901 IF (ITYP .EQ. -2212) THEN
2902 dy1ntp(IY) = dy1ntp(IY) - 1.0
2904 IF (ITYP .EQ. 2212) THEN
2905 dy1ntp(IY) = dy1ntp(IY) + 1.0
2907 c IF (ITYP .EQ. -211 .OR. ITYP .EQ. -321 .OR.
2908 c & ITYP .EQ. -2212) THEN
2909 IF (ITYP .EQ. -2112) THEN
2910 DY1HM(IY) = DY1HM(IY) + 1.0
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
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
2927 IF (ITYP .EQ. 211) THEN
2928 DY1PIP(IY) = DY1PIP(IY) + 1.0
2930 IF (ITYP .EQ. -211) THEN
2931 DY1PIM(IY) = DY1PIM(IY) + 1.0
2933 IF (ITYP .EQ. 111) THEN
2934 DY1PI0(IY) = DY1PI0(IY) + 1.0
2936 IF (ITYP .EQ. 2212) THEN
2937 DY1PR(IY) = DY1PR(IY) + 1.0
2939 IF (ITYP .EQ. -2212) THEN
2940 DY1PB(IY) = DY1PB(IY) + 1.0
2943 IF (ITYP .EQ. 321) THEN
2944 DY1KP(IY) = DY1KP(IY) + 1.0
2946 IF (ITYP .EQ. -321) THEN
2947 DY1KM(IY) = DY1KM(IY) + 1.0
2949 clin-4/24/03 evaluate K0L instead of K0S, since sometimes we may decay K0S:
2950 c IF (ITYP .EQ. 310) THEN
2951 IF (ITYP .EQ. 130) THEN
2952 DY1K0S(IY) = DY1K0S(IY) + 1.0
2954 IF (ITYP .EQ. 3122) THEN
2955 DY1LA(IY) = DY1LA(IY) + 1.0
2957 IF (ITYP .EQ. -3122) THEN
2958 DY1LB(IY) = DY1LB(IY) + 1.0
2960 IF (ITYP .EQ. 333) THEN
2961 DY1PHI(IY) = DY1PHI(IY) + 1.0
2964 c.....insert rapidity cut for mt spectrum here
2965 100 IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
2966 IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
2967 IMT = 1 + int(DXMT / BMT)
2968 IF (ITYP .EQ. 211) THEN
2969 dm1pip(IMT) = dm1pip(IMT) + 1.0 / XMT
2971 IF (ITYP .EQ. -211) THEN
2972 dm1pim(IMT) = dm1pim(IMT) +
2975 IF (ITYP .EQ. 2212) THEN
2976 DMT1PR(IMT) = DMT1PR(IMT) + 1.0 / XMT
2978 IF (ITYP .EQ. -2212) THEN
2979 DMT1PB(IMT) = DMT1PB(IMT) + 1.0 / XMT
2981 IF (ITYP .EQ. 321) THEN
2982 DMT1KP(IMT) = DMT1KP(IMT) + 1.0 / XMT
2984 IF (ITYP .EQ. -321) THEN
2985 dm1km(IMT) = dm1km(IMT) + 1.0 / XMT
2988 c IF (ITYP .EQ. 310) THEN
2989 IF (ITYP .EQ. 130) THEN
2990 dm1k0s(IMT) = dm1k0s(IMT) + 1.0 / XMT
2992 IF (ITYP .EQ. 3122) THEN
2993 DMT1LA(IMT) = DMT1LA(IMT) + 1.0 / XMT
2995 IF (ITYP .EQ. -3122) THEN
2996 DMT1LB(IMT) = DMT1LB(IMT) + 1.0 / XMT
3006 c-----------------------------------------------------------------------
3008 c.....analysis subroutine after the hadronic space-time evolution
3012 PARAMETER (MAXSTR=150001, MAXR=1)
3013 c.....y cut for mt spectrum
3015 c PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3016 PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3018 c.....bin width for mt spectrum and y spectrum
3019 c PARAMETER (BMT = 0.05, BY = 0.2)
3020 PARAMETER (BMT = 0.05, BY = 0.4)
3023 COMMON /ARERC1/MULTI1(MAXR)
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)
3031 c & dm2k0s(50), DMT2LA(50), DMT2LB(50)
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)
3048 DO 1001 I = 1, MULTI1(J)
3055 XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
3056 c 2/24/03 leptons and photons:
3057 if(xm.lt.0.01) goto 200
3058 ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
3059 eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
3061 IF (ABS(PZ) .GE. EE) THEN
3062 PRINT *, 'IN ARTAN2'
3063 PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
3065 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3068 PRINT *, ' PZ = ', PZ, ' EE = ', EE
3070 PRINT *, ' XM = ', XM
3077 Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
3078 c.....rapidity cut for the rapidity distribution
3079 IF (ABS(Y) .GE. 10.0) GOTO 100
3080 c IY = 1 + int(ABS(Y) / BY)
3081 c Ieta = 1 + int(ABS(eta) / BY)
3082 IF (ABS(eta) .GE. 10.0) GOTO 100
3083 IY = 1 + int((Y+10.) / BY)
3084 Ieta = 1 + int((eta+10.) / BY)
3086 IF (ITYP .LT. -1000) THEN
3087 dy2ntb(IY) = dy2ntb(IY) - 1.0
3089 IF (ITYP .GT. 1000) THEN
3090 dy2ntb(IY) = dy2ntb(IY) + 1.0
3092 IF (ITYP .EQ. -2212) THEN
3093 dy2ntp(IY) = dy2ntp(IY) - 1.0
3095 IF (ITYP .EQ. 2212) THEN
3096 dy2ntp(IY) = dy2ntp(IY) + 1.0
3098 IF (ITYP .EQ. -2112) THEN
3099 DY2HM(IY) = DY2HM(IY) + 1.0
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
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
3116 IF (ITYP .EQ. 211) THEN
3117 DY2PIP(IY) = DY2PIP(IY) + 1.0
3119 IF (ITYP .EQ. -211) THEN
3120 DY2PIM(IY) = DY2PIM(IY) + 1.0
3122 IF (ITYP .EQ. 111) THEN
3123 DY2PI0(IY) = DY2PI0(IY) + 1.0
3125 IF (ITYP .EQ. 2212) THEN
3126 DY2PR(IY) = DY2PR(IY) + 1.0
3128 IF (ITYP .EQ. -2212) THEN
3129 DY2PB(IY) = DY2PB(IY) + 1.0
3132 IF (ITYP .EQ. 321) THEN
3133 DY2KP(IY) = DY2KP(IY) + 1.0
3135 IF (ITYP .EQ. -321) THEN
3136 DY2KM(IY) = DY2KM(IY) + 1.0
3139 c IF (ITYP .EQ. 310) THEN
3140 IF (ITYP .EQ. 130) THEN
3141 DY2K0S(IY) = DY2K0S(IY) + 1.0
3143 IF (ITYP .EQ. 3122) THEN
3144 DY2LA(IY) = DY2LA(IY) + 1.0
3146 IF (ITYP .EQ. -3122) THEN
3147 DY2LB(IY) = DY2LB(IY) + 1.0
3149 IF (ITYP .EQ. 333) THEN
3150 DY2PHI(IY) = DY2PHI(IY) + 1.0
3153 c.....insert rapidity cut for mt spectrum here
3154 100 IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
3155 IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
3156 IMT = 1 + int(DXMT / BMT)
3157 IF (ITYP .EQ. 211) THEN
3158 dm2pip(IMT) = dm2pip(IMT) + 1.0 / XMT
3160 IF (ITYP .EQ. -211) THEN
3161 dm2pim(IMT) = dm2pim(IMT) +
3164 IF (ITYP .EQ. 2212) THEN
3165 DMT2PR(IMT) = DMT2PR(IMT) + 1.0 / XMT
3167 IF (ITYP .EQ. -2212) THEN
3168 DMT2PB(IMT) = DMT2PB(IMT) + 1.0 / XMT
3170 IF (ITYP .EQ. 321) THEN
3171 DMT2KP(IMT) = DMT2KP(IMT) + 1.0 / XMT
3173 IF (ITYP .EQ. -321) THEN
3174 dm2km(IMT) = dm2km(IMT) + 1.0 / XMT
3177 c IF (ITYP .EQ. 310) THEN
3178 IF (ITYP .EQ. 130) THEN
3179 dm2k0s(IMT) = dm2k0s(IMT) + 1.0 / XMT
3181 IF (ITYP .EQ. 3122) THEN
3182 DMT2LA(IMT) = DMT2LA(IMT) + 1.0 / XMT
3184 IF (ITYP .EQ. -3122) THEN
3185 DMT2LB(IMT) = DMT2LB(IMT) + 1.0 / XMT
3195 c-----------------------------------------------------------------------
3197 c.....output analysis results at the end of the simulation
3199 SUBROUTINE ARTOUT(NEVNT)
3201 PARAMETER (MAXSTR=150001, MAXR=1)
3202 c.....y cut for mt spectrum
3204 c PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3205 PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3207 c.....bin width for mt spectrum and y spectrum
3208 c PARAMETER (BMT = 0.05, BY = 0.2)
3209 PARAMETER (BMT = 0.05, BY = 0.4)
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)
3218 c & dm1k0s(50), DMT1LA(50), DMT1LB(50)
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)
3233 c & dm2k0s(50), DMT2LA(50), DMT2LB(50)
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)
3247 cms OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN')
3248 cms OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN')
3249 cms OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN')
3250 cms OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN')
3251 cms OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN')
3252 cms OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN')
3253 cms OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN')
3254 cms OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN')
3255 cms OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN')
3257 c OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN')
3258 cms OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN')
3259 cms OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN')
3260 cms OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN')
3261 cms OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN')
3263 cms OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN')
3264 cms OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN')
3265 cms OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN')
3266 cms OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN')
3267 cms OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN')
3268 cms OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN')
3271 cms OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN')
3272 cms OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN')
3273 cms OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN')
3274 cms OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN')
3275 cms OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN')
3276 cms OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN')
3277 cms OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN')
3278 cms OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN')
3279 cms OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN')
3280 clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
3281 c SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0
3282 SCALE1 = 1. / REAL(NEVNT * NUM) / BY
3283 SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1)
3286 ymid=-10.+BY * (I - 0.5)
3287 cms WRITE (30, 333) ymid, SCALE1 * dy1ntb(I)
3288 cms WRITE (31, 333) ymid, SCALE1 * dy1ntp(I)
3289 cms WRITE (32, 333) ymid, SCALE1 * DY1HM(I)
3290 cms WRITE (37, 333) ymid, SCALE1 * DY1KP(I)
3291 cms WRITE (38, 333) ymid, SCALE1 * DY1KM(I)
3292 cms WRITE (39, 333) ymid, SCALE1 * DY1K0S(I)
3293 cms WRITE (40, 333) ymid, SCALE1 * DY1LA(I)
3294 cms WRITE (41, 333) ymid, SCALE1 * DY1LB(I)
3295 cms WRITE (42, 333) ymid, SCALE1 * DY1PHI(I)
3296 cms WRITE (33, 333) ymid, SCALE1 * DY1NEG(I)
3297 cms WRITE (34, 333) ymid, SCALE1 * DY1CH(I)
3298 cms WRITE (35, 333) ymid, SCALE1 * DE1NEG(I)
3299 cms WRITE (36, 333) ymid, SCALE1 * DE1CH(I)
3300 cms WRITE (43, 333) ymid, SCALE1 * dy1msn(I)
3301 cms WRITE (44, 333) ymid, SCALE1 * DY1PIP(I)
3302 cms WRITE (45, 333) ymid, SCALE1 * DY1PIM(I)
3303 cms WRITE (46, 333) ymid, SCALE1 * DY1PI0(I)
3304 cms WRITE (47, 333) ymid, SCALE1 * DY1PR(I)
3305 cms WRITE (48, 333) ymid, SCALE1 * DY1PB(I)
3307 IF (dm1pip(I) .NE. 0.0) THEN
3308 cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I)
3310 IF (dm1pim(I) .NE. 0.0) THEN
3311 cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 *
3314 IF (DMT1PR(I) .NE. 0.0) THEN
3315 cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I)
3317 IF (DMT1PB(I) .NE. 0.0) THEN
3318 cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I)
3320 IF (DMT1KP(I) .NE. 0.0) THEN
3321 cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I)
3323 IF (dm1km(I) .NE. 0.0) THEN
3324 cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3327 IF (dm1k0s(I) .NE. 0.0) THEN
3328 cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I)
3330 IF (DMT1LA(I) .NE. 0.0) THEN
3331 cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I)
3333 IF (DMT1LB(I) .NE. 0.0) THEN
3334 cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I)
3339 cms WRITE (I, *) 'after hadron evolution'
3342 cms WRITE (I, *) 'after hadron evolution'
3346 ymid=-10.+BY * (I - 0.5)
3347 cms WRITE (30, 333) ymid, SCALE1 * dy2ntb(I)
3348 cms WRITE (31, 333) ymid, SCALE1 * dy2ntp(I)
3349 cms WRITE (32, 333) ymid, SCALE1 * DY2HM(I)
3350 cms WRITE (37, 333) ymid, SCALE1 * DY2KP(I)
3351 cms WRITE (38, 333) ymid, SCALE1 * DY2KM(I)
3352 cms WRITE (39, 333) ymid, SCALE1 * DY2K0S(I)
3353 cms WRITE (40, 333) ymid, SCALE1 * DY2LA(I)
3354 cms WRITE (41, 333) ymid, SCALE1 * DY2LB(I)
3355 cms WRITE (42, 333) ymid, SCALE1 * DY2PHI(I)
3356 cms WRITE (33, 333) ymid, SCALE1 * DY2NEG(I)
3357 cms WRITE (34, 333) ymid, SCALE1 * DY2CH(I)
3358 cms WRITE (35, 333) ymid, SCALE1 * DE2NEG(I)
3359 cms WRITE (36, 333) ymid, SCALE1 * DE2CH(I)
3360 cms WRITE (43, 333) ymid, SCALE1 * dy2msn(I)
3361 cms WRITE (44, 333) ymid, SCALE1 * DY2PIP(I)
3362 cms WRITE (45, 333) ymid, SCALE1 * DY2PIM(I)
3363 cms WRITE (46, 333) ymid, SCALE1 * DY2PI0(I)
3364 cms WRITE (47, 333) ymid, SCALE1 * DY2PR(I)
3365 cms WRITE (48, 333) ymid, SCALE1 * DY2PB(I)
3367 IF (dm2pip(I) .NE. 0.0) THEN
3368 cms WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I)
3370 IF (dm2pim(I) .NE. 0.0) THEN
3371 cms WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 *
3374 IF (DMT2PR(I) .NE. 0.0) THEN
3375 cms WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I)
3377 IF (DMT2PB(I) .NE. 0.0) THEN
3378 cms WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I)
3380 IF (DMT2KP(I) .NE. 0.0) THEN
3381 cms WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I)
3383 IF (dm2km(I) .NE. 0.0) THEN
3384 cms WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3387 IF (dm2k0s(I) .NE. 0.0) THEN
3388 cms WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I)
3390 IF (DMT2LA(I) .NE. 0.0) THEN
3391 cms WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I)
3393 IF (DMT2LB(I) .NE. 0.0) THEN
3394 cms WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I)
3397 cms 333 format(2(f12.5,1x))
3402 c-----------------------------------------------------------------------
3404 c.....analysis subroutine in HIJING before parton cascade evolution
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),
3417 DIMENSION dyg1c(50), dmyg1c(50), deyg1c(50)
3418 DIMENSION snrpj1(50), snrtg1(50), snrin1(50),
3420 DIMENSION snyg1c(50), smyg1c(50), seyg1c(50)
3421 DOUBLE PRECISION GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
3425 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3427 COMMON/hjcrdn/YP(3,300),YT(3,300)
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)
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)
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)
3443 COMMON /AREVT/ IAEVT, IARUN, MISS
3449 IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
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)
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)
3498 DO 1006 I = 1, IHNT2(1)
3499 DO 1005 J = 1, NPJ(I)
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
3513 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
3514 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
3536 DO 1008 I = 1, IHNT2(3)
3537 DO 1007 J = 1, NTJ(I)
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
3551 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
3552 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
3575 DO 1009 J = 1, NJSG(I)
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
3589 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
3590 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
3612 DO 1011 I = 1, IHNT2(1)
3613 YR = SQRT(YP(1, I) ** 2 + YP(2, I) ** 2)
3614 IR = 1 + int(YR / DR)
3615 clin-4/2008 protect against out-of-bound errors:
3616 c IF (IR .GT. 50) GOTO 601
3617 IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601
3618 dnrpj1(IR) = dnrpj1(IR) + 1.0
3619 dnrtt1(IR) = dnrtt1(IR) + 1.0
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
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
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
3657 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
3658 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
3671 c.....count number of particles
3672 DO 1016 I = 1, IHNT2(1)
3673 DO 1015 J = 1, NPJ(I)
3675 IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1
3679 DO 1018 I = 1, IHNT2(3)
3680 DO 1017 J = 1, NTJ(I)
3682 IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1
3687 DO 1019 J = 1, NJSG(I)
3689 IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
3693 IF (IOUT .EQ. 1) THEN
3695 c PRINT *, ' in HJANA1 '
3696 c PRINT *, ' total number of partons = ', nsubp
3697 c PRINT *, ' total number of gluons = ', nsubg, MUL
3698 c PRINT *, ' number of projectile strings = ', IHNT2(1)
3699 c PRINT *, ' number of target strings = ', IHNT2(3)
3700 c PRINT *, ' number of independent strings = ', NSG
3701 PRINT *, ' in HJANA1 '
3702 PRINT *, ' total number of partons = ', nsubp / IW
3703 PRINT *, ' total number of gluons = ', nsubg / IW
3704 c PRINT *, ' number of projectile strings = ', IHNT2(1)
3705 c PRINT *, ' number of target strings = ', IHNT2(3)
3706 PRINT *, ' number of independent strings = ', NISG / IW
3713 c-----------------------------------------------------------------------
3715 c.....analysis subroutine in ZPC after generation of additional initial
3716 c.....phase space distributions.
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
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)
3730 COMMON /AREVT/ IAEVT, IARUN, MISS
3737 IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3739 dgxg1a(I) = sgxg1a(I)
3740 dgyg1a(I) = sgyg1a(I)
3745 sgxg1a(I) = dgxg1a(I)
3746 sgyg1a(I) = dgyg1a(I)
3755 IGX = 1 + int(sngl(ABS(GX5(I))) / DGX)
3756 clin-4/2008 protect against out-of-bound errors:
3757 c IF (IGX .GT. 50) GOTO 100
3758 IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100
3759 dgxg1a(IGX) = dgxg1a(IGX) + 1.0
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
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
3775 c-----------------------------------------------------------------------
3777 c.....analysis subroutine in HJAN1A
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
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)
3791 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
3793 COMMON /SREC1/ NSP, NST, NSI
3795 COMMON/hjcrdn/YP(3,300),YT(3,300)
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)
3801 COMMON /AREVT/ IAEVT, IARUN, MISS
3808 IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3810 DNRG1B(I) = SNRG1B(I)
3815 SNRG1B(I) = DNRG1B(I)
3826 IF (J .LE. NSP) THEN
3830 ELSE IF (J .LE. NSP + NST) THEN
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)))
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
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
3854 c-----------------------------------------------------------------------
3856 c.....analysis subroutine in HIJING after parton cascade evolution
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),
3872 DIMENSION dtpj2(50), dttg2(50), dtin2(50),
3874 DIMENSION dyg2c(50), dmyg2c(50), deyg2c(50)
3875 DIMENSION snrpj2(50), snrtg2(50), snrin2(50),
3877 DIMENSION stpj2(50), sttg2(50), stin2(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
3884 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3886 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
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)
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)
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)
3902 COMMON /AREVT/ IAEVT, IARUN, MISS
3906 common/anim/nevent,isoft,isflag,izpc
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)
3916 IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3927 dnrpj2(I) = snrpj2(I)
3928 dnrtg2(I) = snrtg2(I)
3929 dnrin2(I) = snrin2(I)
3930 dnrtt2(I) = snrtt2(I)
3934 dttot2(I) = sttot2(I)
3935 dyg2c(I) = snyg2c(I)
3936 dmyg2c(I) = smyg2c(I)
3937 deyg2c(I) = seyg2c(I)
3953 snrpj2(I) = dnrpj2(I)
3954 snrtg2(I) = dnrtg2(I)
3955 snrin2(I) = dnrin2(I)
3956 snrtt2(I) = dnrtt2(I)
3960 sttot2(I) = dttot2(I)
3961 snyg2c(I) = dyg2c(I)
3962 smyg2c(I) = dmyg2c(I)
3963 seyg2c(I) = deyg2c(I)
3974 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 510
3977 DO 1006 I = 1, IHNT2(1)
3978 DO 1005 J = 1, NPJ(I)
3986 c IF (ABS(PZ) .GE. PE) GOTO 200
3987 IF (ABS(PZ) .GE. PE) THEN
3988 PRINT *, ' IN HJANA2, PROJ STR ', I, ' PART ', J
3989 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3990 PRINT *, ' PZ = ', PZ, ' EE = ', PE
3991 PRINT *, ' XM = ', PM
3995 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
3996 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
4018 DO 1008 I = 1, IHNT2(3)
4019 DO 1007 J = 1, NTJ(I)
4027 c IF (ABS(PZ) .GE. PE) GOTO 400
4028 IF (ABS(PZ) .GE. PE) THEN
4029 PRINT *, ' IN HJANA2, TARG STR ', I, ' PART ', J
4030 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4031 PRINT *, ' PZ = ', PZ, ' EE = ', PE
4032 PRINT *, ' XM = ', PM
4036 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
4037 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
4066 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I)
4077 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
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))
4088 c IF (ABS(PZ) .GE. PE) GOTO 600
4089 IF (ABS(PZ) .GE. PE) THEN
4090 PRINT *, ' IN HJANA2, INDP STR ', I, ' PART ', J
4091 PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4092 PRINT *, ' PZ = ', PZ, ' EE = ', PE
4093 PRINT *, ' XM = ', PM
4097 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
4098 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
4121 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 520
4123 DO 1011 I = 1, IHNT2(1)
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
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
4138 DO 1012 I = 1, IHNT2(3)
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
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
4157 J = I + IHNT2(1) + IHNT2(3)
4159 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) J = I
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
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
4180 PM = sngl(XMASS5(I))
4182 c IF (ABS(PZ) .GE. PE) GOTO 800
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
4193 RAP = 0.5 * LOG((PE + PZ) / (PE - PZ))
4194 XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
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
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
4209 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 530
4211 c.....count number of particles
4212 DO 1016 I = 1, IHNT2(1)
4213 DO 1015 J = 1, NPJ(I)
4215 IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1
4219 DO 1018 I = 1, IHNT2(3)
4220 DO 1017 J = 1, NTJ(I)
4222 IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1
4233 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I)
4240 c IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
4241 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
4242 IF(K2SGS(I, J) .EQ. 21) nsubg = nsubg + 1
4244 IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
4252 IF (IOUT .EQ. 1) THEN
4255 c PRINT *, ' in HJANA2 '
4256 c PRINT *, ' total number of partons = ', nsubp
4257 c PRINT *, ' total number of gluons = ', nsubg, MUL
4258 c PRINT *, ' number of projectile strings = ', IHNT2(1)
4259 c PRINT *, ' number of target strings = ', IHNT2(3)
4260 c PRINT *, ' number of independent strings = ', NSG
4261 PRINT *, ' in HJANA2 '
4262 PRINT *, ' total number of partons = ', nsubp / IW
4263 PRINT *, ' total number of gluons = ', nsubg / IW
4264 c PRINT *, ' number of projectile strings = ', IHNT2(1)
4265 c PRINT *, ' number of target strings = ', IHNT2(3)
4266 PRINT *, ' number of independent strings = ', NISG / IW
4275 c-----------------------------------------------------------------------
4277 c.....subroutine called by HJANA2
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
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)
4293 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4295 COMMON/hjcrdn/YP(3,300),YT(3,300)
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)
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)
4307 COMMON /AREVT/ IAEVT, IARUN, MISS
4314 IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
4316 dgxp2a(I) = sgxp2a(I)
4317 dgyp2a(I) = sgyp2a(I)
4319 dgxg2a(I) = sgxg2a(I)
4320 dgyg2a(I) = sgyg2a(I)
4325 sgxp2a(I) = dgxp2a(I)
4326 sgyp2a(I) = dgyp2a(I)
4328 sgxg2a(I) = dgxg2a(I)
4329 sgyg2a(I) = dgyg2a(I)
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
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
4349 dtp2a(IT) = dtp2a(IT) + 1.0
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
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
4366 dtp2a(IT) = dtp2a(IT) + 1.0
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
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
4385 dtp2a(IT) = dtp2a(IT) + 1.0
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
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
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
4411 c-----------------------------------------------------------------------
4413 c.....analysis subroutine in HJANA2
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
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)
4430 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
4432 COMMON /SREC1/ NSP, NST, NSI
4434 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
4436 COMMON/hjcrdn/YP(3,300),YT(3,300)
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)
4442 COMMON /AREVT/ IAEVT, IARUN, MISS
4449 IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
4451 DNRG2B(I) = SNRG2B(I)
4452 dtg2b(I - 25) = stg2b(I - 25)
4456 SNRG2B(I) = DNRG2B(I)
4457 stg2b(I - 25) = dtg2b(I - 25)
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
4473 TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2))
4474 DTAU=TAU7 - sngl(ATAUI(J))
4475 IT = 1 + int(DTAU / DT)
4477 c IF (ABS(IT) .GT. 25) GOTO 200
4478 IF (IT .GT. 25 .OR. IT .LT. -24) GOTO 200
4480 dtg2b(IT) = dtg2b(IT) + 1.0
4487 c-----------------------------------------------------------------------
4489 c.....analysis subroutine before ARTMN
4492 PARAMETER (MAXSTR=150001, MAXR=1)
4493 c.....y cut for mt spectrum
4494 PARAMETER (YMIN = -1.0, YMAX = 1.0)
4496 c.....bin width for mt spectrum and y spectrum
4497 PARAMETER (DMT = 0.05, DY = 0.2)
4498 DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult,
4500 DIMENSION dndyh3(50), DMYH3(50), DEYH3(50)
4503 COMMON /ARERC1/MULTI1(MAXR)
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)
4513 COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
4520 DO 1001 I = 1, MULTI1(J)
4522 IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200
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
4538 Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
4539 c.....rapidity cut for the rapidity distribution
4540 c IY = 1 + int(ABS(Y) / DY)
4541 IY = 1 + int((Y+10.) / DY)
4542 IF (IY .GT. 50) GOTO 100
4543 dndyh3(IY) = dndyh3(IY) + 1.0
4544 DEYH3(IY) = DEYH3(IY) + XMT
4546 c.....insert rapidity cut for mt spectrum here
4547 IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200
4548 IMT = 1 + int(DXMT / DMT)
4549 IF (IMT .GT. 50) GOTO 200
4550 DMYH3(IMT) = DMYH3(IMT) + 1.0 / XMT
4558 c-----------------------------------------------------------------------
4560 c.....analysis subroutine after ARTMN
4562 PARAMETER (MAXSTR=150001, MAXR=1)
4563 c.....y cut for mt spectrum
4565 c PARAMETER (YMIN = -0.5, YMAX = 0.5)
4566 PARAMETER (YMIN = -1.0, YMAX = 1.0)
4568 c.....bin width for mt spectrum and y spectrum
4569 PARAMETER (DMT = 0.05, DY = 0.2)
4571 DIMENSION dndyh4(50), DMYH4(50), DEYH4(50)
4574 COMMON /ARERC1/MULTI1(MAXR)
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)
4584 COMMON /fflow/ v2f,etf,xmultf,v2fpi,xmulpi
4591 DO 1001 I = 1, MULTI1(J)
4593 IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200
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
4609 Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
4610 c.....rapidity cut for the rapidity distribution
4611 c IY = 1 + int(ABS(Y) / DY)
4612 IY = 1 + int((Y+10.) / DY)
4613 IF (IY .GT. 50) GOTO 100
4614 dndyh4(IY) = dndyh4(IY) + 1.0
4615 DEYH4(IY) = DEYH4(IY) + XMT
4617 c.....insert rapidity cut for mt spectrum here
4618 IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200
4619 IMT = 1 + int(DXMT / DMT)
4620 IF (IMT .GT. 50) GOTO 200
4621 DMYH4(IMT) = DMYH4(IMT) + 1.0 / XMT
4629 c=======================================================================
4631 c.....subroutine to get average values for different strings
4635 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4636 PARAMETER (MAXPTN=400001)
4637 PARAMETER (MAXSTR=150001)
4638 c REAL*4 YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB
4639 REAL YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB
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)
4647 COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
4649 COMMON /SREC1/ NSP, NST, NSI
4651 COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
4653 COMMON/hjcrdn/YP(3,300),YT(3,300)
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)
4660 COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4662 cbz6/28/99 flow1 end
4663 common/anim/nevent,isoft,isflag,izpc
4665 common/strg/np(maxstr)
4667 clin-6/06/02 test local freezeout:
4669 & gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
4670 & pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
4672 & tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
4676 clin-6/06/02 test local freezeout for string melting,
4677 c use space-time values at local freezeout saved in /frzprc/:
4694 DO 1002 I = 1, MAXSTR
4698 clin-4/25/03 add zt3(I) to track longitudinal positions of partons/strings:
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
4712 NSTR = NSP + NST + NSI
4714 clin-7/03/01 correct averaging on transverse coordinates, no shift needed:
4715 if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
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)
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)
4735 IF (I .LE. NSP) THEN
4737 ZT1(I) = dble(YP(1, J))
4738 ZT2(I) = dble(YP(2, J))
4740 ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN
4742 ZT1(I) = dble(YT(1, J))
4743 ZT2(I) = dble(YT(2, J))
4748 1 dble((YP(1, IASG(J, 1)) + YT(1, IASG(J, 2))))
4750 1 dble((YP(2, IASG(J, 1)) + YT(2, IASG(J, 2))))
4759 IF (NP(I).NE.0) THEN
4762 SHIFT=0.5d0*dble(BB)
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
4770 cbz6/28/99 flow1 end
4776 c Initialize hadron weights;
4777 c Can add initial hadrons before the hadron cascade starts (but after ZPC).
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),
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
4791 common /para8/ idpert,npertd,idxsec
4793 c All hadrons at the start of hadron cascade have the weight of 1
4794 c except those inserted by the user in this subroutine:
4799 c Specify number, species, weight, initial x,p,m for inserted hadrons here:
4802 DO 100 i=np0+1,np0+nadd
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))
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)))
4818 c.....give initial formation time shift and boost according to rapidity:
4820 FTAR(I)=TAUI*COSH(RAP)
4821 GXAR(I)=GXAR(I)+VX*TAU0*COSH(RAP)
4822 GYAR(I)=GYAR(I)+VY*TAU0*COSH(RAP)
4823 c Allow the intial z-position to be different from the Bjorken picture:
4824 GZAR(I)=TAUI*SINH(RAP)+GZAR(I)
4825 c GZAR(I)=TAUI*SINH(RAP)
4826 zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
4827 GZAR(I)=GZAR(I)+zsmear
4829 IAINT2(1)=IAINT2(1)+nadd
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'
4836 if(IAINT2(1).gt.MAXSTR) then
4837 write(16,*) 'Too many initial hadrons, array size is exceeded!'