PROGRAM BUU ****************************************************************************** * * * Authors: * * Wolfgang Bauer (BUU) * * and students: * * Bao-An Li (Pion Production) * * Dietrich Klakow (Momentum Dependent Mean Field) * * Catherine Mader (Correlations) * * * *----------------------------------------------------------------------------* * * * REFERENCES: * * W. Bauer et al., Phys. Rev. C 34, 2127 (1986) * * W. Bauer, Nucl. Phys. A471, 604 (1987) * * W. Bauer, Phys. Rev. Lett. 61, 2534 (1988) * * B.A. Li and W. Bauer, Phys. Lett. B254, 335 (1991) * * B.A. Li and W. Bauer, Phys. Rev. C 44, 450 (1991) * * B.A. Li, W. Bauer and G.F. Bertsch, Phys. Rev. C44, 2095 (1991) * * W. Bauer, C.K. Gelbke, S. Pratt, Annu. Rev. Nucl. Part. Sci. * * 42, 77 (1992) * ****************************************************************************** * PARAMETERS: * * MAXPAR - MAXIMUM NUMBER OF TESTPARTICLES PROGRAM CAN HANDLE * * MAXX - NUMBER OF MESHPOINTS IN X AND Y DIRECTION = 2 MAXX + 1 * * MAXZ - NUMBER OF MESHPOINTS IN Z DIRECTION = 2 MAXZ + 1 * * AMU - 1 ATOMIC MASS UNIT "GEV/C**2" * * MX,MY,MZ - MESH SIZES IN COORDINATE SPACE [FM] FOR PAULI LATTICE * * MPX,MPY,MPZ- MESH SIZES IN MOMENTUM SPACE [GEV/C] FOR PAULI LATTICE * *----------------------------------------------------------------------* IMPLICIT NONE INTEGER MAXPAR, MAXX, MAXZ, ISUM, MX, MY, MZ, MPX, MPY, MPZ INTEGER MPZP, MAXR REAL AMU PARAMETER (MAXPAR=100000, AMU=0.938919) PARAMETER (MAXX=40, MAXZ=81, MAXR=1000, ISUM=15) PARAMETER (MX=4, MY=4, MZ=15, MPX=4, MPY=4, MPZ=5, MPZP=10) *----------------------------------------------------------------------* C Read variables: CHARACTER*20 NAMEOUT, NAMERES CHARACTER*50 DIRECT INTEGER ZTA, ZPR, MASSTA, MASSPR, ISEED, NTMAX, ICOLL INTEGER NSTAR, NDIRECT, ICOU, SIGF, NUM, INSYS, IPOT, MODE INTEGER IAVOID, IMOMEN, ICOOR, INOUT, IPCOR, IPICOR INTEGER RESTART, NFREQ, NBIG1, NBIG2, NBIG3, NDMP, NTD(100) REAL ELAB, ZEROPT, B, DT, DBOX, R10, DENCUT, SIGRED REAL DIR, DX, DY, DZ, DPX, DPY, DPZ, CYCBOX c Common block variables: INTEGER MASS, ID, LB, MASSR REAL CALLS, ZET, R, P, E, RHO, GRADX, GRADY, GRADZ REAL GRADPX, GRADPY, GRADPZ CHARACTER*1 FF COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /FILE/ DIRECT, NAMEOUT, NAMERES COMMON /GG/ DX, DY, DZ, DPX, DPY, DPZ COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) COMMON /II/ GRADX(MAXPAR), GRADY(MAXPAR), GRADZ(MAXPAR), 1 GRADPX(MAXPAR), GRADPY(MAXPAR), GRADPZ(MAXPAR) COMMON /INPUT/ NSTAR, NDIRECT, DIR COMMON /MASS/ MASS COMMON /PAUL/ CALLS COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM COMMON /ZZ/ ZTA, ZPR, ZET(11) C 'Local variables' CHARACTER*4 CHAR4 INTEGER IN, LPR1, LNE1, LP1, LP2, LP3, LD1, LD2, LD3 INTEGER LD4, LN1, LN2, NCOLL, NBLOC, NCNNE, NCNND, NCNDN INTEGER NDIRT, NDECAY, NDECBL, NRES, LCOLL, LBLOC, LCNNE INTEGER LCNND, LCNDN, LDIRT, LDECAY, LDECBL, LRES, CHARGE INTEGER NLOST, OUTPAR, LF, IX, IY, IZ INTEGER NT, ISO, NRUN, J, I REAL RADTA, RADPR, ZDIST, ETA, PZTA, BETATA, GAMMTA, EPR REAL PZPR, BETAPR, GAMMPR, S, PSQARE, ETOTAL, GRADX2 REAL GRADY2, GRADZ2, FNORM, PXW, RD REAL PXPI, PYPI, PZPI, EPI C REAL T, V, VP, E0 C REAL GRADXN, GRADYN, GRADZN c Common block for pion correlation function output ! added 2-12 INTEGER nnntot, MAXPI, LBPI, LBD PARAMETER (MAXPI = 100*MAXR) INTEGER LPI(MAXPI), TPI(MAXPI), PICNT(0:MAXPI) REAL RPI(3,MAXPI), PPI(3,MAXPI) COMMON /PIHBT/ IPICOR, nnntot, TPI, PICNT, LPI, RPI, PPI *----------------------------------------------------------------------* * INPUT-SECTION: * * * * 1) TARGET-RELATED QUANTITIES * * MASSTA, zta - TARGET MASS IN AMU, target charge (INTEGER) * * * * 2) PROJECTILE-RELATED QUANTITIES * * MASSPR, zpr - PROJECTILE MASS IN AMU, proj. charge(INTEGER) * * ELAB - BEAM ENERGY IN [MEV/NUCLEON] (REAL) * * ZEROPT - DISPLACEMENT OF THE SYSTEM IN Z-DIREC. [FM](REAL) * * B - IMPACT PARAMETER [FM] (REAL) * * * * 3) PROGRAM-CONTROL PARAMETERS * * ISEED - SEED FOR RANDOM NUMBER GENERATOR (INTEGER) * * DT - TIME-STEP-SIZE [FM/C] (REAL) * * NTMAX - TOTAL NUMBER OF TIMESTEPS (INTEGER) * * dbox - boxsize in fermi (real) * * r10 - width of gausian (real) * * dencut - density cut for coorelations (real) * * iCOLL - (= 1 -> MEAN FIELD ONLY, * * - =-1 -> CACADE ONLY, ELSE FULL BUU) * * sigred - between 0 and 1 => sig_nn ->sigred*signn (REAL) * * Attention: sigred disabled in this version !!!!! * * nstar - 1= include N* resonance (integer) * * ndir - 1= include direct process (integer) * * dir - Percentage of direct pion production (real) * * icou - (= 1 -> coulomb, otherwise not) (integer) * * sigf - 0: cugnon 1: particle data group fit * * NUM - NUMBER OF TESTPARTICLES PER NUCLEON (INTEGER) * * INSYS - (=0 -> LAB-SYSTEM, ELSE C.M. SYSTEM) (INTEGER) * * IPOT - 1 -> SIGMA=2; 2 -> SIGMA=4/3; 3 -> SIGMA=7/6 * * IN MEAN FIELD POTENTIAL (INTEGER) * * MODE - (=1 -> interpolation for pauli-blocking, * * =2 -> local lookup, other -> unblocked)(integer) * * DX,DY,DZ - widths of cell for paulat in coor. sp. [fm](real) * * DPX,DPY,DPZ-widths of cell for paulat in mom. sp.[GeV/c](real) * * IAVOID - (=1 -> AVOID FIRST COLL. WITHIN SAME NUCL. * * =0 -> ALLOW THEM) (INTEGER) * * IMOMEN - FLAG FOR CHOICE OF MOMENTUM DISTRIBUTION * * (=1 -> WOODS-SAXON DENSITY AND LOCAL THOMAS-FERMI * * =2 -> NUCLEAR MATTER DEN. AND LOCAL THOMAS-FERMI * * =3 -> COHERENT BOOST IN Z-DIRECTION) (INTEGER) * * ICOOR - FLAG FOR CHOICE OF COORDINATE DISTRIBUTION * * 1: spheres * * 2: from experimental charge distributions * * CYCBOX - ne.0 => cyclic boundary conditions;boxsize CYCBOX * * INOUT - =1 -> WRITE OUT FINAL NUCLEON MOMENTA (INTEGER) * * =2 -> WRITE OUT FINAL NUCLEON COORDINATES * * =3 -> WRITE OUT MOMENTA AND COORDINATES * * IPCOR - =1 -> WRITE OUT EMISSION TIME, MOMENTA AND * * COORDINATES AT EMISSION TIME (INTEGER) * * IPICOR - =1 -> WRITE OUT EMISSION TIME, MOMENTA AND * * COORDINATES AT EMISSION TIME (INTEGER) * * NAMERES: name of file for reload * * * * 4) CONTROL-PRINTOUT OPTIONS * * NFREQ - NUMBER OF TIMSTEPS AFTER WHICH PRINTOUT * * IS REQUIRED (INTEGER) * * NBIG1 - (=1 -> XZ DENS AT Y=0, =2 -> XZ DENS FOR ALL Y, * * ELSE NO OUTPUT) (INTEGER) * * NBIG2 - (=1 -> PX-PZ DENS AT Y=0, ELSE NO OUTP.)(INTEGER) * * NBIG3 - (=1 -> Z-PZ PHASE-DENS, ELSE NO OUTP.) (INTEGER) * * NDMP - number of times for dump of coordinates and momenta * * NDMP > 0 => all, NDMP < 0 => only pions and reson. * * nt(i) - times for dump * * * *----------------------------------------------------------------------* in=5 read (in,*) direct read (in,*) nameout open (15,status='new',file=direct//nameout//'.d15') cdm for delta masses, remove first 3 characters from all lines cdm beginning with cdm... cdm OPEN (62, STATUS='NEW', FILE=DIRECT//NAMEOUT//'.dMASS') READ (in,*) MASSTA, zta write(isum,*) massta, zta, ' massta, zta' READ (in,*) MASSPR, zpr write(isum,*) masspr, zpr, ' masspr, zpr' READ (in,*) ELAB write(isum,*) elab, ' elab' READ (in,*) ZEROPT write(isum,*) zeropt, ' zeropt' READ (in,*) B write(isum,*) b, ' b' READ (in,*) ISEED write(isum,*) iseed, ' iseed' READ (in,*) DT write(isum,*) dt, ' dt' READ (in,*) NTMAX write(isum,*) ntmax, ' ntmax' READ (in,*) Dbox write(isum,*) dbox, ' dbox' READ (in,*) r10 write(isum,*) r10, ' r10' READ (in,*) dencut write(isum,*) dencut, ' dencut' READ (in,*) icoll write(isum,*) icoll, ' icoll' READ (in,*) sigred write(isum,*) sigred, ' sigred' read (in,*) nstar write (isum,*) nstar,' nstar' read (in,*) ndirect write (isum,*) ndirect,' direct' read (in,*) dir write (isum,*) dir,' dir' READ (in,*) icou write(isum,*) icou, ' icou' READ (in,*) sigf write(isum,*) sigf, ' sigf' READ (in,*) NUM write(isum,*) num, ' num' READ (in,*) INSYS write(isum,*) insys, ' insys' READ (in,*) IPOT write(isum,*) ipot, ' ipot' READ (in,*) mode write(isum,*) mode, ' mode' READ (in,*) dx,dy,dz write(isum,*) dx,dy,dz,' dx,dy,dz' READ (in,*) dpx,dpy,dpz write(isum,*) dpx,dpy,dpz,' dpx,dpy,dpz' READ (in,*) IAVOID write(isum,*) iavoid, ' iavoid' READ (in,*) IMOMEN write(isum,*) imomen, ' imomen' READ (in,*) ICOOR write(isum,*) icoor, ' icoor' READ (in,*) CYCBOX write(isum,*) cycbox, ' cycbox' READ (in,*) inout write(isum,*) inout, ' inout' READ (in,*) ipcor write(isum,*) ipcor, ' ipcor' READ (in,*) ipicor write(isum,*) ipicor, ' ipicor' READ (in,*) restart write(isum,*) restart, ' restart' READ (in,*) nameres write(isum,*) nameres, ' nameres' READ (in,*) NFREQ write(isum,*) nfreq, ' nfreq' READ (in,*) NBIG1 write(isum,*) nbig1, ' nbig1' READ (in,*) NBIG2 write(isum,*) nbig2, ' nbig2' READ (in,*) NBIG3 write(isum,*) nbig3, ' nbig3' READ (in,*) ndmp write(isum,*) ndmp, ' ndmp' do i=1, abs(ndmp) READ (in,*) ntd(i) write(isum,*) i, ntd(i), ' ntd' end do *----------------------------------------------------------------------* write(isum,'(/10x,''Phase space boundaries for Pauli-blocking'')') write(isum,'(10x,f7.3,'' < x < '',f7.3)') -dx*mx,dx*mx write(isum,'(10x,f7.3,'' < y < '',f7.3)') -dy*my,dy*my write(isum,'(10x,f7.3,'' < z < '',f7.3)') -dz*mz,dz*mz write(isum,'(10x,f7.3,'' < px < '',f7.3)') -dpx*mpx,dpx*mpx write(isum,'(10x,f7.3,'' < py < '',f7.3)') -dpy*mpy,dpy*mpy write(isum,'(10x,f7.3,'' < pz < '',f7.3/)') -dpz*mpz,dpz*mpzp CALL FRONT(ISUM,MASSTA,MASSPR,ELAB) *----------------------------------------------------------------------* RADTA = 1.124 * FLOAT(MASSTA)**(1./3.) RADPR = 1.124 * FLOAT(MASSPR)**(1./3.) ZDIST = RADTA + RADPR + 4.0 if ( cycbox.ne.0 ) zdist=0 MASS = MASSTA + MASSPR * IF (num*mass.GT.MAXPAR .or. num.gt.maxr) THEN WRITE(ISUM,'(//10X,''**** FATAL ERROR: TOO MANY TEST PART. **** & '')') STOP 'BUU: too many test particles' END IF * *----------------------------------------------------------------------* * RELATIVISTIC KINEMATICS * * 1) LABSYSTEM * ETA = FLOAT(MASSTA) * AMU PZTA = 0.0 BETATA = 0.0 GAMMTA = 1.0 * EPR = FLOAT(MASSPR) * (AMU + 0.001 * ELAB) PZPR = SQRT( EPR**2 - (AMU * FLOAT(MASSPR))**2 ) BETAPR = PZPR / EPR GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 ) * WRITE(ISUM,'(/27X,''**** KINEMATICAL PARAMETERS ****''/)') WRITE(ISUM,'(27X,''1) LAB-FRAME: TARGET PROJECTILE'')') WRITE(ISUM,'(27X,'' ETOTAL "GEV" '',2F11.4)') ETA, EPR WRITE(ISUM,'(27X,'' P "GEV/C" '',2F11.4)') PZTA, PZPR WRITE(ISUM,'(27X,'' BETA '',2F11.4)') BETATA, BETAPR WRITE(ISUM,'(27X,'' GAMMA '',2F11.4)') GAMMTA, GAMMPR IF (INSYS .NE. 0) THEN * * 2) C.M. SYSTEM * S = (EPR+ETA)**2 - PZPR**2 PSQARE =real( & (1d0*S)**2 + ((1d0*MASSTA)**4 +(1d0*MASSPR)**4)*(1d0*amu)**4 & - 2d0 * S * (1d0*AMU)**2 * (MASSTA**2 + MASSPR**2) & - 2d0 * (1d0*MASSTA**2 * MASSPR**2) * (1d0*AMU)**4) & / (4.0 * S) * ETA = SQRT ( PSQARE + (FLOAT(MASSTA) * AMU)**2 ) PZTA = - SQRT(ABS(PSQARE)) BETATA = PZTA / ETA GAMMTA = 1.0 / SQRT( 1.0 - BETATA**2 ) * EPR = SQRT ( PSQARE + (FLOAT(MASSPR) * AMU)**2 ) PZPR = SQRT(ABS(PSQARE)) BETAPR = PZPR/ EPR GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 ) * WRITE(ISUM,'(27X,''2) C.M.-FRAME: '')') WRITE(ISUM,'(27X,'' ETOTAL "GEV" '',2F11.4)') ETA, EPR WRITE(ISUM,'(27X,'' P "GEV/C" '',2F11.4)') PZTA, PZPR WRITE(ISUM,'(27X,'' BETA '',2F11.4)') BETATA, BETAPR WRITE(ISUM,'(27X,'' GAMMA '',2F11.4)') GAMMTA, GAMMPR WRITE(ISUM,'(27X,''S "GEV**2" '',F11.4)') S WRITE(ISUM,'(27X,''PSQARE "GEV/C"2 '',F11.4)') PSQARE WRITE(ISUM,'(/27X,''*** CALCULATION DONE IN CM-FRAME ***''/)') ELSE WRITE(ISUM,'(/27X,''*** CALCULATION DONE IN LAB-FRAME ***''/)') END IF * MOMENTUM PER PARTICLE PZTA = PZTA / FLOAT(MASSTA) PZPR = PZPR / FLOAT(MASSPR) *----------------------------------------------------------------------- * PROVIDE BETACM FOR QZZ-ROUTINE * BETACM = VELOCITY WITH WHICH CM-SYSTEM MOVES WITH RESPECT TO * SYSTEM IN WHICH HEAVY ION COLLISION IS CALCULATED * c IF(INSYS .NE. 0) THEN c BETACM = 0. c ELSE c SSS = (EPR+ETA)**2 - (PZPR*REAL(MASSPR))**2 c PSQ = real( c & (1d0*SSS)**2 + ((1d0*MASSTA)**4 + (1d0*MASSPR)**4)*amu**4 c & - 2d0 * SSS * (1d0*AMU)**2 * (MASSTA**2 + MASSPR**2) c & - 2d0 * (MASSTA**2 * MASSPR**2) * (1d0*AMU)**4) c & / (4.0 * SSS) c BETACM = SQRT(PSQ) / SQRT ( PSQ + (FLOAT(MASSTA) * AMU)**2 ) c END IF *----------------------------------------------------------------------- * WRITE(ISUM,'(24X,''IMPACT PARAMETER B FOR THIS RUN:'', & F14.6,'' FM'')/)') B * if (mode .eq. 1) then write(isum,'(/t25,''Pauli-principle on the lattice'' & '' - interpolation'')') else if (mode .eq. 2) then write(isum,'(/t25,''Pauli-principle on the lattice'' & '' - local evaluation'')') else write(isum,'(/t25,''No Pauli-principle in effect'')') end if *----------------------------------------------------------------------- * INITIALIZATION OF PHASE-SPACE * call densin(dbox, r10, gammta, gammpr) massr(0)=0 do i=1, num massr(i)=mass end do call coulin(masspr,massta,iseed) call gradin(ipot, icoll) c------------------------------------------------------------------------ CALL INIT(1 ,MASSTA ,RADTA ,B/2. ,ZEROPT+ZDIST/2. , 1 PZTA ,GAMMTA ,ISEED ,MASSTA ,ZTA ,IMOMEN ,ICOOR) * CALL INIT(1+MASSTA,MASS ,RADPR ,-B/2. ,ZEROPT-ZDIST/2. , 1 PZPR ,GAMMPR ,ISEED ,MASSPR ,ZPR ,IMOMEN ,ICOOR) c reload coordinates, momenta, energies and particle identity if ( restart.eq.1 ) call 1 reload(LPR1, LNE1, LP1, LP2, LP3, LD1, LD2, LD3, LD4, LN1, LN2) * OUTPAR = 0 CALL DENS(OUTPAR) * *----------------------------------------------------------------------- * CONTROL PRINTOUT OF INITIAL CONFIGURATION * WRITE(ISUM,'(''1'')') WRITE(ISUM,'(''********** INITIAL CONFIGURATION **********''/)') CALL SUMMRY(ISUM,MASSPR,MASSTA) CALL PRIBIG(ISUM, NBIG1, NBIG2, NBIG3, 0, PZPR, PZTA, dbox) * *----------------------------------------------------------------------- * CALCULATE MOMENTA FOR T = 0.5 * DT (TO OBTAIN 2ND DEGREE * ACCURACY!) * "COMPARE: J. AICHELIN ET AL., PHYS. REV. 31, 1730 (1985)" * calls=0. IF (ICOLL .NE. -1) THEN call platin(mode, num, dx, dy, dz, dpx, dpy, dpz, fnorm) c DO I = 1, num*mass c CALL GRADRU(r(1,i), r(2,i), r(3,i), p(1,i), p(2,i), p(3,i), c 1 GRADXn, GRADYn, GRADZn, lf, insys, dbox) c P(1,I) = P(1,I) - (0.5 * DT) * GRADXn c P(2,I) = P(2,I) - (0.5 * DT) * GRADYn c P(3,I) = P(3,I) - (0.5 * DT) * GRADZn c end do END IF *----------------------------------------------------------------------- * PAULI - INITALIZATION * call platin(mode, num, dx, dy, dz, dpx, dpy, dpz, fnorm) call gradu1(lf, insys) c call rapidd(nt) *5 INITIALIZATION OF TIME-LOOP VARIABLES *5.1 COLLISION COUNTERS NCNNE = 0 NCNND = 0 NCNDN = 0 NCOLL = 0 NBLOC = 0 NDIRT = 0 NDECAY = 0 NRES = 0 NT = 0 IF (IPCOR .EQ. 1) CALL PPINIT * initialize pion phase space for correlation functions ! added 2-12 IF (IPICOR .EQ. 1) THEN nnntot = 0 do j=1,maxpar lPI(j) = 0 tPI(j) = 0 rPI(1,j)= 0 rPI(2,j)= 0 rPI(3,j)= 0 pPI(1,j)= 0 pPI(2,j)= 0 pPI(3,j)= 0 picnt(j)= 0 end do END IF c call energy(t, v, vp, e0, fnorm, dbox) c write (isum,*) 'nt, t, v, vp, e0, etotal' c write (isum,130) nt, t, v, vp, e0,(t-e0+v+vp)/mass*1000 c if ( ipot.eq.4 ) stop '2' write (*,*) ' ==== BUU STARTED ====' write (*,*) write (*,*) ' t [fm] pions deltas N* Col./fm', & ' elast. inel. wPx R' * ============== LOOP OVER ALL TIME STEPS ================ * * * do nt=1,ntmax call countp(lpr1,lne1,lp1,lp2,lp3,ld1,ld2,ld3,ld4,ln1,ln2) IF (ICOLL .NE. 1) THEN * * *----------------------------------------------------------------------- c check time for dump do j=1, abs(ndmp) if ( ntd(j).eq.nt ) then c dump positions ... if ( nt .lt. 10 ) then write(char4,'(''000'',i1)') nt open (60,status='unknown',file=direct//nameout//'.t'//char4) else if ( nt .lt. 100 ) then write(char4,'(''00'',i2)') nt open (60,status='unknown',file=direct//nameout//'.t'//char4) else if ( nt .lt. 1000 ) then write(char4,'(''0'',i3)') nt open (60,status='unknown',file=direct//nameout//'.t'//char4) else if ( nt .lt. 10000 ) then write(char4,'(i4)') nt open (60,status='unknown',file=direct//nameout//'.t'//char4) else char4 = 'AAAA' open (60,status='unknown',file=direct//nameout//'.t'//char4) end if rewind (60) iso=0 if (NDMP .gt. 0) then * write out all particles: do nrun=1,num iso=iso+massr(nrun-1) do i=1+iso,massr(nrun)+iso write(60,'(i3,i5,f8.4,6f10.5)') lb(i),id(i),e(i) 1 ,p(1,i),p(2,i),p(3,i),r(1,i),r(2,i),r(3,i) end do end do else * only write out pions and resonances: do nrun=1,num iso=iso+massr(nrun-1) do i=1+iso,massr(nrun)+iso if (lb(i) .gt. 2) then write(60,'(i3,i5,f8.4,6f10.5,i4)') lb(i),id(i),e(i) 1 ,p(1,i),p(2,i),p(3,i),r(1,i),r(2,i),r(3,i),nrun end if end do end do end if close (60) end if end do *----------------------------------------------------------------------- CALL RELCOL(ISEED, IAVOID, DT, SIGRED, DBOX, LCOLL, LBLOC, 1 LCNNE, LCNND, LCNDN, LDIRT, LDECAY, LDECBL, LRES, NT, SIGF) NCOLL = NCOLL + LCOLL NBLOC = NBLOC + LBLOC NCNNE = NCNNE + LCNNE NCNND = NCNND + LCNND NCNDN = NCNDN + LCNDN NDIRT = NDIRT + LDIRT NDECAY= NDECAY+ LDECAY NDECBL= NDECBL+ LDECBL NRES = NRES + LRES END IF * * UPDATE POSITIONS * if (ipot .gt. 3) call gradu1(lf, insys) charge = 0. iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1 + iso, mass + iso charge = charge + zet(lb(i)) ETOTAL = SQRT( E(I)**2 +P(1,I)**2 +P(2,I)**2 +P(3,I)**2 ) R(1,I) = R(1,I) + DT * (P(1,I) / ETOTAL+gradpx(i)) R(2,I) = R(2,I) + DT * (P(2,I) / ETOTAL+gradpy(i)) R(3,I) = R(3,I) + DT * (P(3,I) / ETOTAL+gradpz(i)) c use cyclic boundary conitions C if ( cycbox.ne.0 ) then C if ( r(1,i).gt. cycbox/2 ) r(1,i)=r(1,i)-cycbox C if ( r(1,i).le.-cycbox/2 ) r(1,i)=r(1,i)+cycbox C if ( r(2,i).gt. cycbox/2 ) r(2,i)=r(2,i)-cycbox C if ( r(2,i).le.-cycbox/2 ) r(2,i)=r(2,i)+cycbox C if ( r(3,i).gt. cycbox/2 ) r(3,i)=r(3,i)-cycbox C if ( r(3,i).le.-cycbox/2 ) r(3,i)=r(3,i)+cycbox C end if end do do i = 1 + mass + iso, massr(nrun) + iso charge = charge + zet(lb(i)) ETOTAL = SQRT( E(I)**2 +P(1,I)**2 +P(2,I)**2 +P(3,I)**2 ) c pions don't see a mean field R(1,I) = R(1,I) + DT * P(1,I) / ETOTAL R(2,I) = R(2,I) + DT * P(2,I) / ETOTAL R(3,I) = R(3,I) + DT * P(3,I) / ETOTAL c use cyclic boundary conitions C if ( cycbox.ne.0 ) then C if ( r(1,i).gt. cycbox/2 ) r(1,i)=r(1,i)-cycbox C if ( r(1,i).le.-cycbox/2 ) r(1,i)=r(1,i)+cycbox C if ( r(2,i).gt. cycbox/2 ) r(2,i)=r(2,i)-cycbox C if ( r(2,i).le.-cycbox/2 ) r(2,i)=r(2,i)+cycbox C if ( r(3,i).gt. cycbox/2 ) r(3,i)=r(3,i)-cycbox C if ( r(3,i).le.-cycbox/2 ) r(3,i)=r(3,i)+cycbox C end if end do end do if ( num*(zta+zpr).ne.charge ) then write (*,*) charge, 'charge is not conserved' stop 'BUU:3' end if * UPDATE THE DELTA, N* AND PION COUNTERS call countp(lpr1,lne1,lp1,lp2,lp3,ld1,ld2,ld3,ld4,ln1,ln2) * * UPDATE DENSITY * CALL DENS(OUTPAR) * IF (IPCOR .EQ. 1) CALL PPCORR(NT, DENCUT, DBOX) * * UPDATE MOMENTA * if (icou .eq. 1) CALL COULPOT(DT) IF (ICOLL .NE. -1) THEN * Mean field interaction NLOST = 0 iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1 + iso, mass + iso * no mean field for pions IX = NINT( R(1,I)/dbox ) IY = NINT( R(2,I)/dbox ) IZ = NINT( R(3,I)/dbox ) IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND. & ABS(IZ) .LT. MAXZ ) THEN call gradu2(ix, iy, iz, gradx2, grady2, gradz2, dbox) P(1,I) = P(1,I) - DT * (GRADX(i) + gradx2) P(2,I) = P(2,I) - DT * (GRADY(i) + grady2) P(3,I) = P(3,I) - DT * (GRADZ(i) + gradz2) ELSE NLOST = NLOST + 1 END IF end do end do * IF (NLOST .NE. 0 .AND. (NT/NFREQ)*NFREQ .EQ. NT) THEN WRITE(ISUM,'(5X,''***'',I7,'' TESTPARTICLES LOST AFTER '', & ''TIME STEP NUMBER'',I4)') NLOST, NT END IF END IF * * update phase space density * call platin(mode, num, dx, dy, dz, dpx, dpy, dpz, fnorm) * * CONTROL-PRINTOUT OF CONFIGURATION (IF REQUIRED) * IF ( (NT/NFREQ)*NFREQ .EQ. NT ) THEN c write other output WRITE(ISUM,'(''1'')') WRITE(ISUM,'(//''********** TIME STEP ='',I5,'' 1 **********''/)') NT WRITE(ISUM,'('' COLLISIONS: ALL ='',I6, 1 ''; BLOCKED COLLISIONS ='',I6,''; ELAST. ='' ,I6,/ 1 '' N+N -> N+DELTA ='',I6,''; N+DELTA '' ,''-> N+N ='', 1 I6,'';N+N->N+N+PION ='',I6,/ '' DELTA OR N* DECAY ='', 1 I6,'';RESONANCE ='',I6)') NCOLL, NBLOC, NCNNE, NCNND, 1 NCNDN, NDIRT, NDECAY, NRES * * C CALL COLCOUNT c call energy(t, v, vp, e0, fnorm, dbox) C call rapidd(nt) CALL SUMMRY(ISUM, MASSPR, MASSTA) CALL PRIBIG(ISUM, NBIG1, NBIG2, NBIG3, NT, PZPR, PZTA, DBOX) END IF * * * Output after each timestep: c call wrstep(nt) CALL SUMRYN(PXW, RD) write(*,'(f8.2,8f9.3)') nt*dt, float(lp1+lp2+lp3)/float(num), 1 float(ld1+ld2+ld3+ld4)/float(num), float(ln1+ln2)/float(num), 1 FLOAT(lcoll)/dt/FLOAT(num), FLOAT(lcnne)/dt/FLOAT(num), 1 FLOAT(lcnnd+lcndn)/dt/FLOAT(num), -1000.*pxw, rd c if ( (rd.gt.2.*zdist).and.(cycbox.eq.0) ) c & stop 'distance of nuclei large enough: done' end do * * * ============== END OF TIME STEP LOOP ================ * *----------------------------------------------------------------------- *----------------------------------------------------------------------- * END OF JOB * IF (IPCOR .EQ. 1) CALL PPRINT(NTMAX, DT, B) *----------------------------------------------------------------------- * write out pion phase space for correlation functions * file .d97 has creation points of emitted pions * file .d98 has creation points of all pions (including reabsorbed * ones) * IF (IPICOR .EQ. 1) THEN c open (97,status='new',file=direct//nameout//'.d97') open (98,status='new',file=direct//nameout//'.d98') i = 0 do nrun = 1, num do iso = 1, massr(nrun) - mass i = i + 1 j = picnt(i) c write(97,'(i7,f6.2,i5,f7.2,6f9.2)') j, b, lPI(j), dt*tPI(j), c 1 (rPI(IX,j),iX=1,3), (pPI(IX,j)*1000.,IX=1,3) write (98,'(2i7,f6.2,i5,f7.2,6f9.2)') j, nrun, b, lPI(j), 1 dt*tPI(j), (rPI(IX,j),iX=1,3), (pPI(IX,j)*1000.,IX=1,3) C write (98,'(i7,f6.2,i5,f7.2,6f9.2)') i, b, lPI(i), dt*tPI(i), C 1 (rPI(IX,i),IX=1,3), (pPI(IX,i)*1000.,IX=1,3) end do end do close (97) close (98) END IF write(*,*) write(*,*) i, ' = number of pions.' *----------------------------------------------------------------------- * WRITE OUT FINAL STATE NUCLEON MOMENTA AND/OR COORDINATES if (inout .eq. 1) then open (94,status='unknown',file=direct//nameout//'.d94') rewind (94) write (94,*) 'timestep',nt iso=0 do nrun=1,num iso=iso+massr(nrun-1) do i=1+iso,massr(nrun)+iso write(94,'(i3,4f8.3)') lb(i),e(i),p(1,i),p(2,i),p(3,i) end do end do else if (inout .eq. 2) then open (94,status='unknown',file=direct//nameout//'.d94') rewind (94) write (94,*) 'timestep',nt iso=0 do nrun=1,num iso=iso+massr(nrun-1) do i=1+iso,massr(nrun)+iso write(94,'(i3,3f7.3)') lb(i),r(1,i),r(2,i),r(3,i) end do end do else if (inout .eq. 3) then open (94,status='unknown',file=direct//nameout//'.d94') rewind (94) write (94,*) 'timestep',nt iso=0 do nrun=1,num iso=iso+massr(nrun-1) do i=1+iso,massr(nrun)+iso write(94,'(i3,i5,f8.4,6f10.5)') lb(i),id(i),e(i) & ,p(1,i),p(2,i),p(3,i) & ,r(1,i),r(2,i),r(3,i) end do end do end if close (94) * UPDATE THE DELTA, N* AND PION COUNTERS call countp(lpr1,lne1,lp1,lp2,lp3,ld1,ld2,ld3,ld4,ln1,ln2) write(*,'(/'' Prot. Neutr. pi+ pi0 pi-'', & '' D++ D+ D0 D- N*+ N*0'')') write (*,'(11i7)') lpr1, lne1, lp1, lp2, lp3, ld1, ld2, ld3, & ld4, ln1, ln2 *----------------------------------------------------------------------- call platout *----------------------------------------------------------------------- cdm close (62) close (15) *----------------------------------------------------------------------- *Write out pion spectrum * open (63,status='new',file=direct//nameout//'.d63') iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1+iso, iso+massr(nrun) if (lb(i) .ge. 6) then lbd = lb(i) call decayt(p(1,i), p(2,i), p(3,i), e(i), lb(i), pxpi, pypi, 1 pzpi, epi, lbpi, iseed) write (63,'(2i7,f6.2,i5,4f9.2,i5)') i, nrun, b, lbpi, 1 1000.0*pxpi, 1000.0*pypi, 1000.0*pzpi, 1000.0*epi, lbd end if if ((lb(i).gt.2) .and. (lb(i).lt.6)) 1 write (63,'(2i7,f6.2,i5,4f9.2,i5)') i, nrun, b, lb(i), 1 (1000.0*p(ix,i), ix=1,3), 1000.0*e(i), lb(i) end do end do close (unit=63) STOP ' ==== BUU COMPLETED ====' END ************************************************************************ SUBROUTINE DIRECT(SRT,ISEED,PX3,PY3,PZ3,PX4,PY4,PZ4,PPX,PPY,PPZ) * PURPOSE : CALCULATE MOMENTUM OF PARTICLES IN THE FINAL STATE FROM * THE DIRECT PROCESS N+N--->N+N+PION * DATE : APRIL 20,1990 ************************************************************************ IMPLICIT NONE INTEGER ISEED REAL SRT, PX3, PY3, PZ3, PX4, PY4, PZ4, PPX, PPY, PPZ REAL AVMASS, AVPI PARAMETER (AVMASS=0.938919, AVPI=0.138037) REAL PM, EM, S, PPM2, PPM, PP3, E3, E41, E42, AL, A, B, C REAL ALFA, BETA, GARMA, PPZ1, PPZ2, PZ41, PZ42, P41, P42 PM = SQRT((SRT-AVPI)**2/4.-AVMASS**2) EM = SRT-2*AVMASS-AVPI S = SRT**2 PPM2 = (S+AVPI**2-4.*AVMASS**2)**2/(4.*S)-AVPI**2 PPM = SQRT(PPM2) C WE CONSIDER THE PROCESS N1+N2-->N3+N4+PION C GENERATE THE MOMENTUM OF N3 51 PX3 = 1.-2.*RAN(ISEED) PY3 = 1.-2.*RAN(ISEED) PZ3 = 1.-2.*RAN(ISEED) PP3 = PX3**2+PY3**2+PZ3**2 IF(PP3.GT.1.) GO TO 51 PX3 = PX3*PM PY3 = PY3*PM PZ3 = PZ3*PM E3 = SQRT(AVMASS**2+PX3**2+PY3**2+PZ3**2) c T3 = E3-AVMASS C GENERATE PX4 AND PY4 PX4 = PM*(1.-2.*RAN(ISEED)) PY4 = PM*(1.-2.*RAN(ISEED)) ALFA = AVMASS**2+PX4**2+PY4**2 C GET PPX AND PPY FOR PION FROM THE MOMENTUM CONSERVATION PPX = -(PX3+PX4) PPY = -(PY3+PY4) BETA = AVPI**2+PPX**2+PPY**2 GARMA = SRT-E3 AL = GARMA**2+BETA-ALFA-PZ3**2 C COEFFICIENTS FOR SOLVING PPZ A = (PZ3/GARMA)**2-1. B = -AL*PZ3/(GARMA**2) C = AL**2/(4.*GARMA**2)-BETA IF((B**2-4.*A*C).LT.0.)GO TO 51 c Find PPZ PPZ1 = (-B+SQRT(B**2-4.*A*C))/(2.*A) PPZ2 = (-B-SQRT(B**2-4.*A*C))/(2.*A) PZ41 = -(PPZ1+PZ3) PZ42 = -(PPZ2+PZ3) P41 = SQRT(PX4**2+PY4**2+PZ41**2) P42 = SQRT(PX4**2+PY4**2+PZ42**2) E41 = SQRT(PX4**2+PY4**2+PZ41**2+AVMASS**2)-AVMASS E42 = SQRT(PX4**2+PY4**2+PZ42**2+AVMASS**2)-AVMASS c EP1 = SQRT(AVPI**2+PPX**2+PPY**2+PPZ1**2) c EP2 = SQRT(AVPI**2+PPX**2+PPY**2+PPZ2**2) IF(ABS(PPZ1).GT.PPM.OR.ABS(PPZ2).GT.PPM)GO TO 10 IF(P41.GT.PM.OR.P42.GT.PM)GO TO 10 IF(E41.GT.EM.OR.E42.GT.EM)GO TO 10 * BOTH SOLUTIONS ARE ALL RIGHT. WE TAKE THE FIRST SOLUTION ONLY PPZ = PPZ1 PZ4 = PZ41 10 IF (ABS(PPZ1) .LE. PPM .AND. P41 .LE. PM .AND. ABS(PPZ2) .GT. PPM 1 .OR. P42 .GT. PM) THEN PPZ = PPZ1 PZ4 = PZ41 ENDIF IF (ABS(PPZ2) .LE. PPM .AND. P42 .LE. PM .AND. ABS(PPZ1) .GT. PPM 1 .OR. P41 .GT. PM) THEN PPZ = PPZ2 PZ4 = PZ42 ENDIF RETURN END ************************************************************************ subroutine dirpro(i1, i2, iseed, irun, iblock, srt, ntag) * * * purpose: direct proecess N+N->N+N+pion * * Warning: this part has never been properly checked !!! * ************************************************************************ IMPLICIT NONE INTEGER I1, I2, ISEED, IBLOCK, NTAG, IRUN REAL SRT INTEGER MAXPAR, MAXP, MAXR, MX, MY, MZ, MPX, MPY, MPZ, MPZP REAL AMN, AMP, AP1, AP2 PARAMETER (MAXPAR=100000, MAXP=100, MAXR=1000, AMN=0.939566, 1 AMP=0.938272, AP1=0.134974, AP2=0.139568) PARAMETER (MX=4, MY=4, MZ=15, MPX=4, MPY=4, MPZ=5, MPZP=10) INTEGER NNN, ID, LB, LPION, IPION REAL DX, DY, DZ, DPX, DPY, DPZ, BETAX, BETAY, BETAZ, GAMMA REAL RPION, PPION, EPION, R, P, E CHARACTER*1 FF COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) COMMON /GG/ DX, DY, DZ, DPX, DPY, DPZ COMMON /NN/ NNN COMMON /BG/ BETAX, BETAY, BETAZ, GAMMA COMMON /PA/ RPION(3,MAXP,MAXR) COMMON /PB/ PPION(3,MAXP,MAXR) COMMON /PC/ EPION(MAXP,MAXR) COMMON /PD/ IPION(MAXP,MAXR), LPION(MAXP,MAXR) c local variables INTEGER LB1, LB2, ID1, ID2 INTEGER IX1, IY1, IZ1, IPX1P, IPY1P, IPZ1P INTEGER IX2, IY2, IZ2, IPX2P, IPY2P, IPZ2P INTEGER IPX1, IPY1, IPZ1, IPX2, IPY2, IPZ2 REAL EM1, EM2, SIG1, SIG2, PX1, PY1, PZ1, PX2, PY2, PZ2 REAL PX3, PY3, PZ3, PX4, PY4, PZ4, PPX, PPY, PPZ, PPBETA REAL E1CM, E2CM, P1BETA, P2BETA, TRANSF, OCCUP, EPCM c functions REAL SIGMA PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) PX2 = P(1,I2) PY2 = P(2,I2) PZ2 = P(3,I2) LB1 = LB(I1) LB2 = LB(I2) *keep all coordinates for possible phase space change ix1 = nint(r(1,i1)/dx) iy1 = nint(r(2,i1)/dy) iz1 = nint(r(3,i1)/dz) ix2 = nint(r(1,i2)/dx) iy2 = nint(r(2,i2)/dy) iz2 = nint(r(3,i2)/dz) ipx1 = nint(px1/dpx) ipy1 = nint(py1/dpy) ipz1 = nint(pz1/dpz) ipx2 = nint(px2/dpx) ipy2 = nint(py2/dpy) ipz2 = nint(pz2/dpz) CALL DIRECT(SRT,ISEED,PX3,PY3,PZ3,PX4,PY4,PZ4,PPX,PPY,PPZ) * FIND THE MOMENTUM OF PARTICLES IN THE FINAL STATE IN THE NUCLEUS- * NUCLEUS CMS. FRAME AND CHECKING THE PAULI BLOCKING. * WAS COLLISION PAULI-FORBIDEN? IF YES, NTAG = -1 NTAG = 0 * LORENTZ-TRANSFORMATION INTO LAB FRAME E1CM = SQRT (EM1**2 + PX3**2 + PY3**2 + PZ3**2) P1BETA = PX3*BETAX + PY3*BETAY + PZ3*BETAZ TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) + E1CM ) P(1,I1) = BETAX * TRANSF + PX3 P(2,I1) = BETAY * TRANSF + PY3 P(3,I1) = BETAZ * TRANSF + PZ3 *CHECK IF PARTICLE #1 IS PAULI BLOCKED CALL PAULat(I1,occup) if (ran(iseed) .lt. occup) then ntag = -1 else ntag = 0 end if *IF PARTICLE #1 IS NOT PAULI BLOCKED IF (NTAG .NE. -1) THEN E2CM = SQRT (EM2**2 + PX4**2 + PY4**2 + PZ4**2) P2BETA = PX4*BETAX+PY4*BETAY+PZ4*BETAZ TRANSF = GAMMA * (GAMMA*P2BETA / (GAMMA + 1.) + E2CM) P(1,I2) = BETAX * TRANSF + PX4 P(2,I2) = BETAY * TRANSF + PY4 P(3,I2) = BETAZ * TRANSF + PZ4 *CHECK IF PARTICLE #2 IS PAULI BLOCKED CALL PAULat(I2,occup) if (ran(iseed) .lt. occup) then ntag = -1 else ntag = 0 end if * END IF * IF COLLISION IS BLOCKED,RESTORE THE MOMENTUM COORDINATES * AND LABELS OF I1 AND I2 IF (NTAG .EQ. -1) THEN P(1,I1) = PX1 P(2,I1) = PY1 P(3,I1) = PZ1 P(1,I2) = PX2 P(2,I2) = PY2 P(3,I2) = PZ2 E(I1) = EM1 E(I2) = EM2 LB(I1) = LB1 LB(I2) = LB2 ELSE PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) id1 = id(i1) id2 = id(i2) if ( id1.ne.0 ) ID(I1) = id1+id1/abs(id1) if ( id2.ne.0 ) ID(I2) = id2+id2/abs(id2) IBLOCK = 4 NNN = NNN + 1 * DETERMINE THE CHARGE STATES OF PARTICLES AFTER DIRECT PION PRODUCTION * (1) FOR P+P IF(LB(I1)*LB(I2).EQ.1)THEN IF (RAN(ISEED) .LT. SIGMA(2,1,SRT)/(2.*SIGMA(2,1,SRT) + 1 SIGMA(3,1,SRT)))THEN * (1.1)P+P-->P+P+PION(0) LPION(NNN,IRUN) = 4 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP1 ELSE * IN PION(0) PRODUCTION PROCESS NO CHARGE STATE CHANGE * (1.2)P+P -->P+N+PION(+) LPION(NNN,IRUN) = 5 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP2 IF(RAN(ISEED).LE.0.5)THEN LB(I1) = 1 LB(I2) = 2 E(I1) = AMP E(I2) = AMN ELSE LB(I1) = 2 LB(I2) = 1 E(I1) = AMN E(I2) = AMP ENDIF ENDIF ENDIF * (2)FOR N+N IF(LB(I1)*LB(I2).EQ.4)THEN IF (RAN(ISEED) .LT. SIGMA(2,1,SRT)/(2.*SIGMA(2,1,SRT) + 1 SIGMA(3,1,SRT)))THEN * (2.1)N+N-->N+N+PION(0) LPION(NNN,IRUN) = 4 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP1 ELSE * (2.2)N+N-->N+P+PION(-) LPION(NNN,IRUN) = 3 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP2 IF(RAN(ISEED).GT.0.5)THEN LB(I1) = 2 LB(I2) = 1 E(I1) = AMN E(I2) = AMP ELSE LB(I1) = 1 LB(I2) = 2 E(I1) = AMP E(I2) = AMN ENDIF ENDIF ENDIF * (3)FOR N+P IF(LB(I1)*LB(I2).EQ.2)THEN SIG1 = SIGMA(3,1,SRT)+SIGMA(4,2,SRT) SIG2 = SIGMA(2,1,SRT)+SIGMA(4,2,SRT) IF(RAN(ISEED).LT.SIG1/(SIG1+2.*SIG2))THEN * (3.1)N+P-->N+P+PION(0) LPION(NNN,IRUN) = 4 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP1 ELSE IF(RAN(ISEED).LT.0.5)THEN * (3.2)N+P-->N+N+PION(+) LPION(NNN,IRUN) = 5 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP2 LB(I1) = 2 LB(I2) = 2 E(I1) = AMN E(I2) = AMN ELSE * (3.3)N+P-->P+P+PION(-) LPION(NNN,IRUN) = 3 IPION(NNN,IRUN) = 0 EPION(NNN,IRUN) = AP2 LB(I1) = 1 LB(I2) = 1 E(I1) = AMP E(I2) = AMP ENDIF ENDIF ENDIF EM1 = E(I1) EM2 = E(I2) EPCM = SQRT(EPION(NNN,IRUN)**2+PPX**2+PPY**2+PPZ**2) PPBETA = PPX*BETAX+PPY*BETAY+PPZ*BETAZ TRANSF = GAMMA*(GAMMA*PPBETA/(GAMMA+1.)+EPCM) PPION(1,NNN,IRUN) = BETAX*TRANSF+PPX PPION(2,NNN,IRUN) = BETAY*TRANSF+PPY PPION(3,NNN,IRUN) = BETAZ*TRANSF+PPZ RPION(1,NNN,IRUN) = R(1,I1) RPION(2,NNN,IRUN) = R(2,I1) RPION(3,NNN,IRUN) = R(3,I1) if ((abs(ix1).le.mx).and.(abs(iy1).le.my).and.(abs(iz1).le.mz)) 1 then ipx1p = nint(p(1,i1)/dpx) ipy1p = nint(p(2,i1)/dpy) ipz1p = nint(p(3,i1)/dpz) if ((ipx1p.ne.ipx1) .or. (ipy1p.ne.ipy1) .or. (ipz1p.ne.ipz1)) 1 then if ((abs(ipx1).le.mpx) .and. (abs(ipy1).le.my) .and. 1 (ipz1.ge.-mpz) .and. (ipz1.le.mpzp)) 1 ff(ipz1,ipy1,ipx1,ix1,iy1,iz1) = 1 char(ichar(ff(ipz1,ipy1,ipx1,ix1,iy1,iz1))-1) if ((abs(ipx1p).le.mpx) .and. (abs(ipy1p).le.my) .and. 1 (ipz1p.ge.-mpz).and. (ipz1p.le.mpzp)) 1 ff(ipz1p,ipy1p,ipx1p,ix1,iy1,iz1) = 1 char(ichar(ff(ipz1p,ipy1p,ipx1p,ix1,iy1,iz1))+1) ipx1 = IPX1P ipy1 = IPY1P ipz1 = IPZ1P end if end if if ((abs(ix2).le.mx) .and. (abs(iy2).le.my) .and. 1 (abs(iz2).le.mz)) then ipx2p = nint(p(1,i2)/dpx) ipy2p = nint(p(2,i2)/dpy) ipz2p = nint(p(3,i2)/dpz) if ((ipx2p.ne.ipx2) .or. (ipy2p.ne.ipy2) .or. (ipz2p.ne.ipz2)) 1 then if ((abs(ipx2).le.mpx) .and. (abs(ipy2).le.my) .and. 1 (ipz2.ge.-mpz) .and. (ipz2.le.mpzp)) 1 ff(ipz2,ipy2,ipx2,ix2,iy2,iz2) = 1 char(ichar(ff(ipz2,ipy2,ipx2,ix2,iy2,iz2))-1) if ((abs(ipx2p).le.mpx) .and. (abs(ipy2p).le.my) .and. 1 (ipz2p.ge.-mpz) .and. (ipz2p.le.mpzp)) 1 ff(ipz2p,ipy2p,ipx2p,ix2,iy2,iz2) = 1 char(ichar(ff(ipz2p,ipy2p,ipx2p,ix2,iy2,iz2))+1) ipx2 = IPX2P ipy2 = IPY2P ipz2 = IPZ2P end if end if ENDIF RETURN end ************************************************************************ block data particles * set particle ID's, charges and masses implicit none integer proton, neutron, pionm, pion0, pionp, zta, zpr integer deltam, delta0, deltap, deltapp, ns0, nsp real masspa(11) real zet common /zz/ zta, zpr, zet(11) common /particle/ masspa, proton, neutron, pionm, pion0, pionp, 1 deltam, delta0, deltap, deltapp, ns0, nsp data zet /1.0 , 0.0 ,-1.0 , 0.0 , 1.0 , 1 -1.0 , 0.0 , 1.0 , 2.0 , 0.0 ,1.0/ data masspa/0.938272, 0.939566, 0.139568, 0.134974, 0.139568, 1 1.232 , 1.232 , 1.232 , 1.232 , 1.47 ,1.47/ data proton /1/ data neutron /2/ data pionm /3/ data pion0 /4/ data pionp /5/ data deltam /6/ data delta0 /7/ data deltap /8/ data deltapp /9/ data ns0 /10/ data nsp /11/ end ************************************************************************ SUBROUTINE BABACO(LB1, I1, LB2, I2, PXCM, PYCM, PZCM, SIGF, DONFL, 1 ISEED, IRUN, LDIRT, LBLOC, LCOLL, NTAG, IBLOCK, 1 SRT, DT, EM1, EM2, SIGRED, DBOX) * * * purpose: do baryon-baryon collision * * IBLOCK = 0 ; NOTHING HAS HAPPENED * IBLOCK = 1 ; ELASTIC N-N COLLISION * IBLOCK = 2 ; N + N -> N + DELTA * IBLOCK = 3 ; N + DELTA -> N + N * IBLOCK = 4 ; N + N _> N + N + PION,DIRECT PROCESS ************************************************************************ IMPLICIT NONE INTEGER LB1, I1, LB2, I2, SIGF, DONFL, ISEED, IRUN, LDIRT INTEGER LBLOC, LCOLL, NTAG, IBLOCK REAL SRT, DT, EM1, EM2, SIGRED, DBOX, PXCM, PYCM, PZCM INTEGER MAXPAR, MAXX, MAXZ REAL PI, AMP, AP1 PARAMETER (MAXPAR=100000, MAXX=40, MAXZ=81, PI=3.14159, 1 AMP=0.938272, AP1=0.134974) REAL A, B, C, SIGMA, LAMBDA, RHO0 REAL BETAX, BETAY, BETAZ, GAMMA, R, P, RHO COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /BG/ BETAX, BETAY, BETAZ, GAMMA COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 c local variables INTEGER IC REAL SIG, SIGNN, DS, DELTAR, TRANSF, ECM, PBETA REAL CUTOFF, NC, OCCUP c functions REAL PPELAST, NPTOTAL, NPINEL, NUNUEL LOGICAL NUCLEON, DELTA, RESONANCE * Prepare the elastic cross section for baryon+baryon collisions * com: parametrization of energy-dependence of n-n-xsect. * with low-energy-cutoff if (sigf.eq.1 .and. nucleon(lb1).and.nucleon(lb2) ) then c particle data group fit for nucleons CUTOFF = 1.89 sig = 100. if ( srt.lt.cutoff ) then signn = sig else if ( lb1.eq.lb2 ) then signn = ppelast(srt) else signn = nptotal(srt) - npinel(srt) end if end if if ( signn.ge.sig ) signn = sig else c cugnon- Cutoff energy determined by mass of lightest nucleon CUTOFF = (2*AMP + 0.02) SIG = 55. SIGNN = nunuel(srt,cutoff) END IF c reduce crosssection c write (*,*) (1.+sigred*(rho(nint(r(1,i1)/dbox), c & nint(r(2,i1)/dbox),nint(r(3,i1)/dbox))/rho0)) c signn = (1.+sigred*(rho(nint(r(1,i1)/dbox), c & nint(r(2,i1)/dbox),nint(r(3,i1)/dbox))/rho0))*signn donfl = 1 DS = SQRT(SIG/(10.*PI)) DELTAR = 1.8 IF ( (nucleon(LB1) .AND. resonance(LB2)) .OR. 1 (resonance(LB1) .AND. nucleon(LB2)) ) then * If nucleon+baryon resonance collision CUTOFF = (2*amp + ap1 + 0.02)**2 else IF ( (nucleon(LB1) .AND. nucleon(LB2)) .or. 1 (resonance(LB1) .AND. resonance(LB2))) then * if nucleon+nucleon or baryon resonance+baryon resonance collision * low energy cutoff is determined by lightest nucleon CUTOFF = (2*AMP + 0.02)**2 else stop 'collision not possible' end if * Check distance between particles and type of collision CALL DISTANCE (I1, I2, DELTAR, DS, DT, CUTOFF, IC, PXCM,PYCM,PZCM) IF (IC.EQ.-1) RETURN CALL CROSS(IRUN, PXCM, PYCM, PZCM, SRT, I1, I2, ISEED, IBLOCK, 1 NTAG, SIGNN, SIG) IF (IBLOCK.EQ.0) RETURN * A COLLISION HAS TAKEN PLACE !! LCOLL = LCOLL +1 *COM: FOR DIRECT PROCESS WE HAVE TREATED THE PAULI BLOCKING AND FIND * THE MOMENTUM OF PARTICLES IN THE ''LAB'' FRAME. SO RETURN TO * RELCOL AND GOTO LINE 400 IF (IBLOCK.EQ.4) THEN IF (NTAG.EQ.-1) THEN LBLOC = LBLOC+1 ELSE LDIRT = LDIRT+1 END IF return ENDIF * Otherwise, was collision pauli-forbiden? if yes, ntag = -1 NTAG = 0 * * Lorentz-transformation to lab frame & check if part 1 is pauli blocked ECM = SQRT (EM1**2 + PXCM**2 + PYCM**2 + PZCM**2) PBETA = PXCM*BETAX + PYCM*BETAY + PZCM*BETAZ TRANSF = GAMMA * ( GAMMA * PBETA / (GAMMA + 1) + ECM ) P(1,I1) = BETAX * TRANSF + PXCM P(2,I1) = BETAY * TRANSF + PYCM P(3,I1) = BETAZ * TRANSF + PZCM CALL PAULAT(I1, OCCUP) IF (RAN(ISEED) .LT. OCCUP) NTAG = -1 * if part 1 is not Pauli blocked, check if part 2 is Pauli blocked IF (NTAG .NE. -1) THEN ECM = SQRT (EM2**2 + PXCM**2 + PYCM**2 + PZCM**2) TRANSF = GAMMA * (-GAMMA*PBETA / (GAMMA + 1.) + ECM) P(1,I2) = BETAX * TRANSF - PXCM P(2,I2) = BETAY * TRANSF - PYCM P(3,I2) = BETAZ * TRANSF - PZCM CALL PAULAT(I2, OCCUP) IF (RAN(ISEED) .LT. OCCUP) NTAG = -1 END IF DONFL=0 return end ************************************************************************ subroutine changph(i, ix, iy, iz, ipx, ipy, ipz, am, em) * * * purpose: change wigner function f(r,p) after scattering * * of particle i * ************************************************************************ implicit none integer i, ix, iy, iz, ipx, ipy, ipz real am, em integer maxpar, mx, my, mz, mpx, mpy, mpz, mpzp PARAMETER (MAXPAR=100000) parameter (MX=4, MY=4, mz=15, MPX=4, MPY=4, mpz=5, mpzp=10) real p, dx, dy, dz, dpx, dpy, dpz character*1 ff COMMON /BB/ P(3,MAXPAR) common /hh/ ff(-mpz:mpzp,-mpy:mpy,-mpx:mpx,-mx:mx,-my:my,-mz:mz) common /gg/ dx, dy, dz, dpx, dpy, dpz integer ipxp, ipyp, ipzp if ((abs(ix).le.mx).and.(abs(iy).le.my).and.(abs(iz).le.mz)) then ipxp = nint(p(1,i)/dpx) ipyp = nint(p(2,i)/dpy) ipzp = nint(p(3,i)/dpz) if ((ipxp.ne.ipx) .or. (ipyp.ne.ipy) .or. (ipzp.ne.ipz)) then if ((abs(ipx).le.mpx) .and. (abs(ipy).le.mpy) .and. 1 (ipz.ge.-mpz) .and. (ipz.le.mpzp) .AND. (AM.LT.1.)) 1 ff(ipz,ipy,ipx,ix,iy,iz) = 1 char(ichar(ff(ipz,ipy,ipx,ix,iy,iz)) - 1) ipx = ipxp ipy = ipyp ipz = ipzp if ((abs(ipx).le.mpx) .and. (abs(ipy).le.mpy) .and. 1 (ipz.ge.-mpz) .and. (ipz.le.mpzp) .AND. (EM.LT.1.)) 1 ff(ipz,ipy,ipx,ix,iy,iz) = 1 char(ichar(ff(ipz,ipy,ipx,ix,iy,iz)) + 1) end if end if return end ************************************************************************ SUBROUTINE CMS(I1, I2, PX1CM, PY1CM, PZ1CM, SRT) * * calculate: - CM Momentum and Energy of pair of colliding particles * _ Lorentz transformation to reference system of pair ************************************************************************ IMPLICIT NONE INTEGER I1, I2 REAL PX1CM, PY1CM, PZ1CM, SRT INTEGER MAXPAR PARAMETER (MAXPAR=100000) REAL BETAX, BETAY, BETAZ, GAMMA, P, E COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /BG/ BETAX, BETAY, BETAZ, GAMMA REAL PX1, PY1, PZ1, PX2, PY2, PZ2, EM1, EM2, E1, E2 REAL S, ETOTAL, P1BETA, TRANSF PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) PX2 = P(1,I2) PY2 = P(2,I2) PZ2 = P(3,I2) EM1 = E(I1) EM2 = E(I2) E1 = SQRT(EM1**2+PX1**2+PY1**2+PZ1**2) E2 = SQRT(EM2**2 + PX2**2 + PY2**2 + PZ2**2 ) S = (E1+E2)**2-(PX1+PX2)**2-(PY1+PY2)**2-(PZ1+PZ2)**2 SRT = SQRT(S) *LORENTZ-TRANSFORMATION IN I1-I2-C.M. SYSTEM ETOTAL = E1 + E2 BETAX = (PX1+PX2) / ETOTAL BETAY = (PY1+PY2) / ETOTAL BETAZ = (PZ1+PZ2) / ETOTAL GAMMA = 1.0 / SQRT(1.0-BETAX**2-BETAY**2-BETAZ**2) *TRANSFORMATION OF MOMENTA (PX1CM = - PX2CM) P1BETA = PX1*BETAX + PY1*BETAY + PZ1 * BETAZ TRANSF = GAMMA * ( GAMMA * P1BETA / (GAMMA + 1) - E1 ) PX1CM = BETAX * TRANSF + PX1 PY1CM = BETAY * TRANSF + PY1 PZ1CM = BETAZ * TRANSF + PZ1 RETURN END ************************************************************************ SUBROUTINE COLCOUNT ************************************************************************ IMPLICIT NONE INTEGER ISUM, MAXPAR, MAXR PARAMETER (MAXPAR=100000, MAXR=1000, ISUM=15) INTEGER NUM, MASS, MASSR, ID, LB COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /MASS/ MASS COMMON /RUN/ NUM COMMON /RR/ MASSR(0:MAXR) INTEGER ISTOSS, JSTOSS, ISO, NRUN, STOSS(0:49) DO ISTOSS = 0,49 STOSS(ISTOSS) = 0 END DO ISO = 0 DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO ISTOSS = 1+ISO, MASS+ISO JSTOSS = ABS(ID(ISTOSS))-1 IF (JSTOSS .GT. 49) JSTOSS = 49 IF (JSTOSS .lT. 0) JSTOSS = 0 STOSS(JSTOSS) = STOSS(JSTOSS) + 1 END DO END DO WRITE(ISUM,'('' DISTRIBUTION OF COLLISION NUMBERS:'')') WRITE(ISUM,'(10(I7,'':'',I5))') (ISTOSS,STOSS(ISTOSS),ISTOSS=0,49) RETURN END ************************************************************************ SUBROUTINE COULIN(MASSPR,MASSTA,ISEED) * * * purpose: initialization of array lb() for all runs * * lb(i) = 1 => proton * * lb(i) = 2 => neutron * ************************************************************************ IMPLICIT NONE INTEGER MASSPR, MASSTA, ISEED INTEGER MAXPAR PARAMETER (MAXPAR=100000) INTEGER ZTA, ZPR, ID, LB, MASS, NUM REAL ZET COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /ZZ/ ZTA, ZPR, ZET(11) COMMON /MASS/ MASS COMMON /RUN/ NUM INTEGER I, IRUN, J1, J2, TEMPLB, IMIN, IMAX, MASSR DO IRUN=1,NUM MASSR = (IRUN-1)*MASS IMIN = MASSR + 1 IMAX = MASSR + ZTA DO I = IMIN, IMAX LB(i) = 1 END DO IMIN = IMAX+1 IMAX = MASSR + MASSTA DO I = IMIN, IMAX LB(i) = 2 END DO IMIN = IMAX+1 IMAX = MASSR + MASSTA + ZPR DO I = IMIN, IMAX LB(i) = 1 END DO IMIN = IMAX+1 IMAX = MASSR + MASSTA + MASSPR DO I = IMIN, IMAX LB(i) = 2 END DO c mix neutrons and protons do i=1, 100 IMIN = MASSR + 1 IMAX = MASSR + MASSTA do j1 = imin, IMAX j2=nint(ran(iseed)*(massta-1)+1.)+MASSR TEMPLB=lb(j1) lb(j1)=lb(j2) lb(j2)=TEMPLB end do IMIN = MASSR + MASSTA + 1 IMAX = MASSR + MASSTA + MASSPR do j1 = IMIN, IMAX j2=nint(ran(iseed)*(masspr- 1)+1.+massta)+MASSR TEMPLB=lb(j1) lb(j1)=lb(j2) lb(j2)=TEMPLB end do end do end do return end ************************************************************************ SUBROUTINE COULPOT(DT) * * Finds Coulomb potential ************************************************************************ IMPLICIT NONE INTEGER MAXR, MAXPAR PARAMETER (MAXR=1000, MAXPAR=100000) INTEGER NUM, ZTA, ZPR, MASSR, ID, LB REAL DT, ZET, R, P COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM COMMON /ZZ/ ZTA, ZPR, ZET(11) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) INTEGER ISO, IRUN, IL, I, J, JL, IDIR REAL GRP, RDIFF, DDX, DDY, DDZ, temp(3,maxr) iso = 0 do irun = 1, num iso = iso + massr(irun-1) do il = 1, massr(irun) temp(1,il) = 0. temp(2,il) = 0. temp(3,il) = 0. end do do il = 1, massr(irun) i = iso + il if (zet(lb(i)).ne.0) then do jl = 1, il-1 j = iso + jl if (zet(lb(j)).ne.0) then ddx=r(1,i)-r(1,j) ddy=r(2,i)-r(2,j) ddz=r(3,i)-r(3,j) rdiff = sqrt(ddx**2+ddy**2+ddz**2) if (rdiff .le. 1.) rdiff = 1. grp=zet(lb(i))*zet(lb(j))/rdiff**3 ddx=ddx*grp ddy=ddy*grp ddz=ddz*grp temp(1,il)=temp(1,il)+ddx temp(2,il)=temp(2,il)+ddy temp(3,il)=temp(3,il)+ddz temp(1,jl)=temp(1,jl)-ddx temp(2,jl)=temp(2,jl)-ddy temp(3,jl)=temp(3,jl)-ddz end if end do end if end do do il = 1, massr(irun) I = iso + il if (zet(lb(i)).ne.0) then do idir = 1,3 p(idir,i) = p(idir,i) + temp(idir,il)* dt * 0.00144 end do end if end do end do return end ************************************************************************ subroutine countp(lpr1,lne1,lp1,lp2,lp3,ld1,ld2,ld3,ld4,ln1,ln2) * * * purpose: count Nucleons, Pions, Deltas and N*'s * ************************************************************************ IMPLICIT NONE INTEGER LPR1, LNE1, LP1, LP2, LP3, LD1, LD2, LD3, LD4, LN1 INTEGER MAXR, MAXPAR, LN2 PARAMETER (MAXR=1000, MAXPAR=100000) INTEGER NUM, MASSR, ID, LB COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM INTEGER ISO, NRUN, I, LB1 *5.1 NUCLEON COUNTERS LPR1 = 0 LNE1 = 0 *5.2 PION COUNTERS : LP1,LP2 AND LP3 ARE THE NO. OF P+,P0 AND P- LP1 = 0 LP2 = 0 LP3 = 0 *5.3 DELTA COUNTERS : LD1,LD2,LD3 AND LD4 ARE THE NO. OF D++,D+,D0 AND D- LD1 = 0 LD2 = 0 LD3 = 0 LD4 = 0 *5.4 N*(1440) COUNTERS : LN1 AND LN2 ARE THE NO. OF N*+ AND N*0 LN1 = 0 LN2 = 0 iso = 0 do nrun = 1,num iso = iso + massr(nrun-1) do i = 1 + iso, massr(nrun) + iso LB1 = LB(I) IF ( LB1.EQ.1) then LPR1 = LPR1 + 1 else IF ( LB1.EQ.2) then LNE1 = LNE1 + 1 else IF ( LB1.EQ.9) then LD1 = LD1 + 1 else IF ( LB1.EQ.8) then LD2 = LD2 + 1 else IF ( LB1.EQ.7) then LD3 = LD3 + 1 else IF ( LB1.EQ.6) then LD4 = LD4 + 1 else IF ( LB1.EQ.11) then LN1 = LN1 + 1 else IF ( LB1.EQ.10) then LN2 = LN2 + 1 else IF ( LB1.EQ.5) then LP1 = LP1 + 1 else IF ( LB1.EQ.4) then LP2 = LP2 + 1 else IF ( LB1.EQ.3) then LP3 = LP3 + 1 end if end do end do return end ************************************************************************ SUBROUTINE CROSS(IRUN, PX, PY, PZ, SRT, I1, I2, ISEED, IBLOCK, 1 NTAG, SIGNN, SIG) * PURPOSE: * * DEALING WITH BARYON-BARYON COLLISIONS * * NOTE : * * VALID ONLY FOR BARYON-BARYON-DISTANCES LESS THAN 1.32 FM * * (1.32 = 2 * HARD-CORE-RADIUS [HRC] ) * * QUANTITIES: * * PX,PY,PZ - MOMENTUM COORDINATES OF ONE PARTICLE IN CM FRAME* * (output) * * SRT - SQRT OF S * * NSTAR =1 INCLUDING N* RESONANCE,ELSE NOT * * NDIRECT=1 INCLUDING DIRECT PION PRODUCTION PROCESS * * IBLOCK - THE INFORMATION BACK * * 0-> COLLISION CANNOT HAPPEN * * 1-> N-N ELASTIC COLLISION * * 2-> N+N->N+DELTA,OR N+N->N+N* REACTION * * 3-> N+DELTA->N+N OR N+N*->N+N REACTION * * 4-> N+N->N+N+PION,DIRTCT PROCESS * * N12 - IS USED TO SPECIFY BARYON-BARYON REACTION * * CHANNELS. M12 IS THE REVERSAL CHANNEL OF N12 * * N12, * * M12=1 FOR p+n-->delta(+)+ n * * 2 p+n-->delta(0)+ p * * A 3 p+p-->delta(++)+n * * 4 p+p-->delta(+)+p * * 5 n+n-->delta(0)+n * * 6 n+n-->delta(-)+p * * 7 n+p-->N*(0)+p * * 8 n+p-->N*(+)+n * * NOTE ABOUT N* RESONANCE: * * As it has been discussed in VerWest's paper,I= 1 (initial isospin) * channel can all be attributed to delta RESONANCE while I= 0 * * channel can all be attribured to N* RESONANCE.Only in n+p * * one can have I=0 channel so is the N* RESONANCE * * REFERENCES: J. CUGNON ET AL., NUCL. PHYS. A352, 505 (1981) * * Y. KITAZOE ET AL., PHYS. LETT. 166B, 35 (1986) * * B. VerWest el al., PHYS. PRV. C25 (1982)1979 * * Gy.Wolf ET AL., PREPRINT UNIV. GISSEN -90-3 * * CUTOFF = 2 * amp + 20 MEV * * * ************************************************************************ IMPLICIT NONE INTEGER IRUN, I1, I2, ISEED, IBLOCK, NTAG REAL PX, PY, PZ, SRT, SIGNN, SIG INTEGER MAXPAR REAL PI, AVMASS, amp, ap1 PARAMETER (MAXPAR=100000, PI=3.14159, AVMASS=0.938919, 1 amp=0.938272, ap1=0.134974) INTEGER NSTAR, NDIRECT, LB, ID REAL DIR, E COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /INPUT/ NSTAR, NDIRECT, DIR c local variables INTEGER N12, M12, DIRFL REAL EM1, EM2, PR, C2, X1, AS, A, TA, X, T1, T2, C1, DM REAL PR2, S1, S2, CT1, CT2, ST1, ST2, SS, PRF2, RENORM c functions LOGICAL NUCLEON, RESONANCE *----------------------------------------------------------------------- IBLOCK = 0 NTAG = 0 EM1 = E(I1) EM2 = E(I2) PR = SQRT( PX**2 + PY**2 + PZ**2 ) C2 = PZ / PR X1 = RAN(ISEED) *----------------------------------------------------------------------- *COM: Test for elastic scattering (N-N, DELTA-DELTA, N-DELTA OR N*-N*) c param. of Cugnon paper IF (X1 .LE. SIGNN/SIG) THEN AS = ( 3.65 * (SRT - 1.8766) )**6 A = 6.0 * AS / (1.0 + AS) TA = -2.0 * PR**2 X = RAN(ISEED) T1 = ALOG( EXP(A*TA) - X*(2.)*exp(A*ta/2.)*sinh(a*ta/2.) )/A C1 = 1.0 - T1/TA T1 = 2.0 * PI * RAN(ISEED) IBLOCK = 1 ELSE * Test for inelastic scattering * If the available energy is less than the pion-mass, nothing * can happen any more ==> RETURN (2*amp + ap1 + 0.02) * if there were 2 delta and they didn't scatt. elast., RETURN IF (SRT .LT. 2*amp + ap1 + 0.02) RETURN IF ( resonance(lb(i1)) .AND. resonance(lb(i2)) ) RETURN IF ( NUCLEON(LB(I1)).AND.NUCLEON(LB(I2))) THEN CALL NNDNCH(I1, I2, SRT, SIG, SIGNN, X1, NDIRECT, DIR, NSTAR, 1 DIRFL, N12, ISEED) IF ( N12.EQ.0 ) RETURN IF ( DIRFL.EQ.0 ) THEN c N+N->N+Resonance * Determine delta or N* mass via rejection method. CALL RESMASS(N12, ISEED, SRT, DM) cdm WRITE(62,*) DM, n12 CALL RELBAR (I1, I2, DM, N12) EM1 = E(I1) EM2 = E(I2) PR2 = (SRT**2 - EM1**2 - EM2**2)**2 - 4.0 * (EM1*EM2)**2 IF (PR2.LE.0.) PR2 = 0.00000001 PR = SQRT(PR2)/(2.*SRT) IBLOCK = 2 else C Direct process; dirpro changes momenta of baryons! call dirpro(i1, i2, iseed, irun, iblock, srt, ntag) return end if ELSE c N+Resonance-> N+N PRF2 = 0.25*SRT**2-AVMASS**2 IF (EM1.GT.1.) THEN RENORM = EM1*prf2/pr ELSE RENORM = EM2*prf2/pr END IF * DETERMINE THE REACTION CHANNELS IN THE FOLLOWING call ndnnch(i1, i2, srt, sig, signn, x1, nstar, m12, renorm) IF ( M12.EQ.0 ) RETURN CALL RELNUC (I1, I2, M12) PR = sqrt(PRF2) IBLOCK=3 END IF C1 = 1.0 - 2.0 * RAN(ISEED) T1 = 2.0 * PI * RAN(ISEED) END IF *COM: SET THE NEW MOMENTUM COORDINATES IF (PX .EQ. 0.0 .AND. PY .EQ. 0.0) THEN T2 = 0.0 ELSE T2=ATAN2(PY,PX) END IF S1 = SQRT( 1.0 - C1**2 ) S2 = SQRT( 1.0 - C2**2 ) CT1 = COS(T1) ST1 = SIN(T1) CT2 = COS(T2) ST2 = SIN(T2) SS = C2 * S1 * CT1 + S2 * C1 PX = PR * ( SS*CT2 - S1*ST1*ST2 ) PY = PR * ( SS*ST2 + S1*ST1*CT2 ) PZ = PR * ( C1*C2 - S1*S2* CT1 ) RETURN END ************************************************************************ SUBROUTINE DECAY(IRUN, I, NNN, ISEED) * PURPOSE:1. SORT DELTA OR N* DECAY PRODUCTS * 2. DETERMINE THE MOMENTUM AND COORDINATES OF NUCLEON AND PION * AFTER THE DELTA OR N* DECAYING * DATE : JAN. 24, 1990 ************************************************************************ IMPLICIT NONE INTEGER IRUN, I, NNN, ISEED INTEGER MAXPAR, MAXP, MAXR REAL AP1, AP2, AMP, AMN PARAMETER (MAXPAR=100000, MAXP=100, MAXR=1000, AMN=0.939566, 1 AMP=0.938272, AP1=0.134974, AP2=0.139568) INTEGER ID, LB, LPION, IPION REAL EPION COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /PC/ EPION(MAXP,MAXR) COMMON /PD/ IPION(MAXP,MAXR), LPION(MAXP,MAXR) REAL X, AM * DETERMINE THE DECAY PRODUCTS * 1. FOR DELTA++ DECAY IF(LB(I).EQ.9)THEN LB(I) = 1 AM = AMP LPION(NNN,IRUN) = 5 EPION(NNN,IRUN) = AP2 * 2. FOR DELTA+ DECAY ELSE IF(LB(I).EQ.8)THEN X = RAN(ISEED) IF(X.LT.(1./3.))THEN LB(I) = 2 AM = AMN LPION(NNN,IRUN) = 5 EPION(NNN,IRUN) = AP2 ELSE LB(I) = 1 AM = AMP LPION(NNN,IRUN) = 4 EPION(NNN,IRUN) = AP1 ENDIF * 3. FOR DELTA0 DECAY ELSE IF(LB(I).EQ.7) THEN X = RAN(ISEED) IF(X.LT.(1./3.))THEN LB(I) = 1 AM = AMP LPION(NNN,IRUN) = 3 EPION(NNN,IRUN) = AP2 ELSE LB(I) = 2 AM = AMN LPION(NNN,IRUN) = 4 EPION(NNN,IRUN) = AP1 ENDIF * 4. FOR DELTA- DECAY ELSE IF(LB(I).EQ.6) THEN LB(I) = 2 AM = AMN LPION(NNN,IRUN) = 3 EPION(NNN,IRUN) = AP2 *5. FOR N*+(1440) DECAY ELSE IF(LB(I).EQ.11)THEN X = RAN(ISEED) IF(X.GT.(1./3.))THEN LB(I) = 2 AM = AMN LPION(NNN,IRUN) = 5 EPION(NNN,IRUN) = AP2 ELSE LB(I) = 1 AM = AMP LPION(NNN,IRUN) = 4 EPION(NNN,IRUN) = AP1 ENDIF *6 FOR N*0(1440) DECAY ELSE IF(LB(I).EQ.10)THEN X = RAN(ISEED) IF(X.GT.(1./3.))THEN LB(I) = 1 AM = AMP LPION(NNN,IRUN) = 3 EPION(NNN,IRUN) = AP2 ELSE LB(I) = 2 AM = AMN LPION(NNN,IRUN) = 4 EPION(NNN,IRUN) = AP1 ENDIF ENDIF CALL DKINEMATICS(IRUN, I, NNN, AM, ISEED) IPION(NNN,IRUN) = ID(I) RETURN END ************************************************************************ SUBROUTINE DECAYT(p1x,p1y,p1z,e1,lb1,p2x,p2y,p2z,e2,lb2,ISEED) * * PURPOSE:1. SORT DELTA OR N* DECAY PRODUCTS * 2. DETERMINE THE MOMENTUM AND COORDINATES OF NUCLEON AND PION * AFTER THE DELTA OR N* DECAYING * This routine combines the routines decay and dkinematics. Since it * doesn't over-write the existing arrays unless the array variables are * sent in the call statement. This routine is used to determine pion * spectra at the end of the run. p1, e1 are for baryon variables and * p2, e2 are the pion coordinates (cmm 2-26-93) * DATE : JAN. 24,1990 ************************************************************************ IMPLICIT NONE INTEGER ISEED, LB2, LB1 real P1X, P1Y, P1Z, E1, P2X, P2Y, P2Z, E2 REAL AMN, AMP, AP1, AP2, PI PARAMETER (AMN=0.939566, AMP=0.938272, AP1=0.134974, 1 AP2=0.139568, PI=3.14159) INTEGER I, NUM, LPION, IPION, X REAL PX, PY, PZ, EDELTA, DM, P2, Q2, EP, EN REAL BDX, BDY, BDZ, BP PX = P1X PY = P1Y PZ = P1Z DM = E1 EDELTA = SQRT(DM**2 + PX**2 + PY**2 + PZ**2) BDX = PX/DM BDY = PY/DM BDZ = PZ/DM * DETERMINE THE DECAY PRODUCTS * 1. FOR DELTA++ DECAY IF(LB1.EQ.9)THEN LB1 = 1 E1 = AMP LB2 = 5 E2 = AP2 * 2. FOR DELTA+ DECAY ELSE IF(LB1.EQ.8)THEN X = RAN(ISEED) IF(X.LT.(1./3.))THEN LB1 = 2 E1 = AMN LB2 = 5 E2 = AP2 ELSE LB1 = 1 E1 = AMP LB2 = 4 E2 = AP1 ENDIF * 3. FOR DELTA0 DECAY ELSE IF(LB1.EQ.7) THEN X = RAN(ISEED) IF(X.LT.(1./3.))THEN LB1 = 1 E1 = AMP LB2 = 3 E2 = AP2 ELSE LB1 = 2 E1 = AMN LB2 = 4 E2 = AP1 ENDIF * 4. FOR DELTA- DECAY ELSE IF(LB1.EQ.6) THEN LB1 = 2 E1 = AMN LB2 = 3 E2 = AP2 *5. FOR N*+(1440) DECAY ELSE IF(LB1.EQ.11)THEN X = RAN(ISEED) IF(X.GT.(1./3.))THEN LB1 = 2 E1 = AMN LB2 = 5 E2 = AP2 ELSE LB1 = 1 E1 = AMP LB2 = 4 E2 = AP1 ENDIF *6 FOR N*0(1440) DECAY ELSE IF(LB1.EQ.10)THEN X = RAN(ISEED) IF(X.GT.(1./3.))THEN LB1 = 1 E1 = AMP LB2 = 3 E2 = AP2 ELSE LB1 = 2 E1 = AMN LB2 = 4 E2 = AP1 ENDIF ENDIF Q2 = ((DM**2 - E1**2 + E2**2)/(2.*DM))**2 - E2**2 IF(Q2.LE.0.)Q2 = 0.0000001 11 PX = 1. - 2.*RAN(ISEED) PY = 1. - 2.*RAN(ISEED) PZ = 1. - 2.*RAN(ISEED) P2 = PX**2 + PY**2 + PZ**2 IF(P2.GT.1.) GO TO 11 P2 = SQRT(Q2/P2) PX = P2*PX PY = P2*PY PZ = P2*PZ EP = SQRT(Q2 + E2**2) EN = SQRT(Q2 + E1**2) BP = (BDX*PX + BDY*PY + BDZ*PZ)/(1. + EDELTA/DM) P1X =-PX + BDX*(-BP + EN) P1Y =-PY + BDY*(-BP + EN) P1Z =-PZ + BDZ*(-BP + EN) P2X = PX + BDX*( BP + EP) P2Y = PY + BDY*( BP + EP) P2Z = PZ + BDZ*( BP + EP) RETURN END ************************************************************************ SUBROUTINE DENS(NESC) * * * PURPOSE: CALCULATION OF NUCLEAR DENSITY FROM SPATIAL * * DISTRIBUTION OF TESTPARTICLES * * * * VARIABLES (ALL INPUT, ALL INTEGER) * * NUM - NUMBER OF TESTPARTICLES PER NUCLEON * * MASS - MASS of system * * NESC - NUMBER OF ESCAPED PARTICLES (INTEGER,OUTPUT) * * dbox - size of box in coordinatespace * * * ************************************************************************ c in this Version: lorentz contraction of gaussians IMPLICIT NONE INTEGER NESC, MAXPAR, MAXX, MAXZ, DMAX, MAXR PARAMETER (MAXPAR=100000, MAXR=1000) PARAMETER (MAXX=40, MAXZ=81, DMAX=8) * INTEGER IX, IY, IZ, I, LX, LY, LZ, ISO, NRUN * INTEGER NUM, MASS, LMAX, MASSR, ID, LB REAL DBOXI, NORMTA, NORMPR, DRHOTA, DRHOPR, R, RHO COMMON /AA/ R(3,MAXPAR) COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /EE/ ID(MAXPAR), lb(maxpar) COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM COMMON /MASS/ MASS COMMON /SMEAR/ LMAX, DBOXI, NORMTA, NORMPR 1 ,DRHOTA(-DMAX:DMAX,-DMAX:DMAX,-DMAX:DMAX) 1 ,DRHOPR(-DMAX:DMAX,-DMAX:DMAX,-DMAX:DMAX) * DO IZ = -MAXZ,MAXZ DO IY = -MAXX,MAXX DO IX = -MAXX,MAXX RHO(IX,IY,IZ) = 0.0 end do end do end do * NESC = 0 * iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1+iso, mass+iso IX = NINT( R(1,I) *dboxi) IY = NINT( R(2,I) *dboxi) IZ = NINT( R(3,I) *dboxi) IF( IX .LE. -MAXX+lmax .OR. IX .GE. MAXX-lmax .OR. & IY .LE. -MAXX+lmax .OR. IY .GE. MAXX-lmax .OR. & IZ .LE. -MAXZ+lmax .OR. IZ .GE. MAXZ-lmax ) THEN NESC = NESC + 1 ELSE if ( id(i).ge.1 ) then c add to the contribution of the testparticle to the density; target do lz = iz-lmax, iz+lmax do ly = iy-lmax, iy+lmax do lx = ix-lmax, ix+lmax rho(lx,ly,lz) = rho(lx,ly,lz) + 1 drhota(lx-ix,ly-iy,lz-iz)*normta end do end do end do else c add to the contribution of the testparticle to the density; projectile do lz = iz-lmax, iz+lmax do ly = iy-lmax, iy+lmax do lx = ix-lmax, ix+lmax rho(lx,ly,lz) = rho(lx,ly,lz) + 1 drhopr(lx-ix,ly-iy,lz-iz)*normpr end do end do end do end if END IF end do end do * RETURN END ************************************************************************ SUBROUTINE DENSIN(DBOX, R10, GAMMTA, GAMMPR) * * This subroutine initializes the nucleon distributions used in * DENS. It could be made into an entry of the subroutine DENS and * the common block smear could be removed. * ************************************************************************ IMPLICIT NONE INTEGER DMAX REAL DBOX, R10, GAMMTA, GAMMPR PARAMETER (DMAX=8) INTEGER LMAX, NUM REAL DBOXI, NORMTA, NORMPR, DRHOTA, DRHOPR COMMON /RUN/ NUM COMMON /SMEAR/ LMAX, DBOXI, NORMTA, NORMPR 1 ,DRHOTA(-DMAX:DMAX,-DMAX:DMAX,-DMAX:DMAX) 1 ,DRHOPR(-DMAX:DMAX,-DMAX:DMAX,-DMAX:DMAX) INTEGER LZ, LY, LX REAL FACTOR, LZGAMTA2, LZGAMPR2, LY2, LX2 * IF (DBOX.LE.R10/(DMAX-1)) STOP 'SUBROUTINE dens: BOX TOO SMALL' DBOXI = 1. / DBOX LMAX = INT(R10*DBOXI + 0.5) NORMTA = 0.0 NORMPR = 0.0 FACTOR = -0.6631*(DBOX*1.87/R10)**2 DO LZ = -LMAX, LMAX LZGAMTA2 = 1.0*LZ*LZ*GAMMTA**2*FACTOR LZGAMPR2 = 1.0*LZ*LZ*GAMMPR**2*FACTOR DO LY = -LMAX, LMAX LY2 = 1.0*LY*LY*FACTOR DO LX = -LMAX, LMAX LX2 = 1.0*LX*LX*FACTOR c each testparticle is smeared using gaussian, that falls from .9 c to .1 in the same way as the surface of the nucleus DRHOTA(LX,LY,LZ) = EXP(LX2 + LY2 + LZGAMTA2) IF ( DRHOTA(LX,LY,LZ).LT.0.08 ) DRHOTA(LX,LY,LZ)=0. NORMTA = NORMTA + DRHOTA(LX,LY,LZ) DRHOPR(LX,LY,LZ) = EXP(LX2 + LY2 + LZGAMPR2) IF ( DRHOPR(LX,LY,LZ).LT.0.08 ) DRHOPR(LX,LY,LZ)=0. NORMPR = NORMPR + DRHOPR(LX,LY,LZ) END DO END DO END DO NORMTA = DBOXI**3/NUM/NORMTA NORMPR = DBOXI**3/NUM/NORMPR RETURN END ************************************************************************ SUBROUTINE DISTANCE(I1,I2,DELTAR,DS,DT,EC,IC,PX1CM,PY1CM,PZ1CM) * PURPOSE : CHECK IF THE COLLISION BETWEEN TWO PARTICLES CAN HAPPEN * BY CHECKING * (1) IF THE DISTANCE BETWEEN THEM IS SMALLER * THAN THE MAXIMUM DISTANCE DETERMINED FROM THE CROSS SECTION. * (2) IF PARTICLE WILL PASS EACH OTHER WITHIN * TWO HARD CORE RADIUS. * (3) IF PARTICLES WILL GET CLOSER. * VARIABLES : * IC=1 COLLISION HAPPENED * IC=-1 COLLISION CAN NOT HAPPEN ************************************************************************ IMPLICIT NONE INTEGER I1, I2, IC REAL DELTAR, DS, DT, EC, PX1CM, PY1CM, PZ1CM INTEGER MAXPAR, I PARAMETER (MAXPAR=100000) REAL R, P, E, BETAX, BETAY, BETAZ, GAMMA REAL GRADX, GRADY, GRADZ, GRADPX, GRADPY, GRADPZ COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /II/ GRADX(MAXPAR), GRADY(MAXPAR), GRADZ(MAXPAR), 1 GRADPX(MAXPAR), GRADPY(MAXPAR), GRADPZ(MAXPAR) COMMON /BG/ BETAX, BETAY, BETAZ, GAMMA REAL DX, DY, DZ, PX1, PX2, PY1, PY2, PZ1, PZ2 REAL EM1, EM2, E1, E2, S, PRCM, DRBETA, TRANSF REAL DXCM, DYCM, DZCM, DZZN, DZZ, BBB, RSQARE, DRCM c SAVE I C if ( i.eq.0 ) i=555 IC = -1 DZ = R(3,I1) - R(3,I2) IF (ABS(DZ) .GT. DELTAR) RETURN DY = R(2,I1) - R(2,I2) IF (ABS(DY) .GT. DELTAR) RETURN DX = R(1,I1) - R(1,I2) IF (ABS(DX) .GT. DELTAR) RETURN RSQARE = DX**2 + DY**2 + DZ**2 IF (RSQARE .GT. DELTAR**2) RETURN * NOW THERE IS ENOUGH ENERGY AVAILABLE ! PRCM = SQRT (PX1CM**2 + PY1CM**2 + PZ1CM**2) IF (PRCM .LE. 0.00001) RETURN PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) PX2 = P(1,I2) PY2 = P(2,I2) PZ2 = P(3,I2) EM1 = E(I1) EM2 = E(I2) *NOW PARTICLES ARE CLOSE ENOUGH TO EACH OTHER ! E1 = SQRT(EM1**2 + PX1**2 + PY1**2 + PZ1**2) E2 = SQRT(EM2**2 + PX2**2 + PY2**2 + PZ2**2) S = (e1+e2)**2-(px1+px2)**2-(py1+py2)**2-(pz1+pz2)**2 * LOW ENERGY CUTOFF ! IF (S .LT. ec) RETURN * TRANSFORMATION OF SPATIAL DISTANCE DRBETA = BETAX*DX + BETAY*DY + BETAZ*DZ TRANSF = GAMMA * GAMMA * DRBETA / (GAMMA + 1) DXCM = BETAX * TRANSF + DX DYCM = BETAY * TRANSF + DY DZCM = BETAZ * TRANSF + DZ * DETERMINING IF THIS IS THE POINT OF CLOSEST APPROACH DRCM = SQRT (DXCM**2 + DYCM**2 + DZCM**2 ) DZZ = (PX1CM*DXCM + PY1CM*DYCM + PZ1CM*DZCM) / PRCM BBB = DRCM**2 - DZZ**2 * WILL PARTICLE PASS EACH OTHER WITHIN 2 * HARD CORE RADIUS ? IF (BBB .GT. DS*DS) RETURN dx = dx + DT * (px1/e1 - px2/e2 + gradpx(i1) - gradpx(i2)) dy = dy + DT * (py1/e1 - py2/e2 + gradpy(i1) - gradpy(i2)) dz = dz + DT * (pz1/e1 - pz2/e2 + gradpz(i1) - gradpz(i2)) DRBETA = BETAX*DX + BETAY*DY + BETAZ*DZ TRANSf = GAMMA * GAMMA * DRBETa / (GAMMA + 1) DXCM = BETAX * TRANSf + DX DYCM = BETAY * TRANSf + DY DZCM = BETAZ * TRANSf + DZ DZZN = (PX1CM*DXCM + PY1CM*DYCM + PZ1CM*DZCM)/PRCM * WILL PARTICLES GET CLOSER ? c if ( ran(i).gt. .975 ) IC = 1 IF ( dzzn*dzz.gt.0. ) RETURN IC=1 RETURN END ************************************************************************ SUBROUTINE DKINEMATICS(IRUN, I, NNN, AM, ISEED) * PURPOSE: * CALCULATE THE MOMENTUM OF NUCLEON AND PION IN THE LAB. FRAME * AFTER DELTA OR N* DECAY * DATE : JAN. 24, 1990 ************************************************************************ IMPLICIT NONE INTEGER IRUN, I, NNN, ISEED REAL AM INTEGER MAXPAR, MAXP, MAXR PARAMETER (MAXPAR=100000, MAXP=100, MAXR=1000) REAL P, E, PPION, EPION COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /PB/ PPION(3,MAXP,MAXR) COMMON /PC/ EPION(MAXP,MAXR) REAL PX, PY, PZ, EDELTA, DM, Q, Q2, QX, QY, QZ REAL EP, EN, BDX, BDY, BDZ, BP, PM PX = P(1,I) PY = P(2,I) PZ = P(3,I) DM = E(I) EDELTA = SQRT(DM**2 + PX**2 + PY**2 + PZ**2) PM = EPION(NNN,IRUN) Q2 = ((DM**2-AM**2 + PM**2)/(2.*DM))**2-PM**2 IF(Q2.LE.0.)Q2 = 0.0000001 11 QX = 1.-2.*RAN(ISEED) QY = 1.-2.*RAN(ISEED) QZ = 1.-2.*RAN(ISEED) Q = QX**2 + QY**2 + QZ**2 IF(Q.GT.1.) GO TO 11 Q = SQRT(Q2/Q) QX = Q*QX QY = Q*QY QZ = Q*QZ EP = SQRT(Q2 + PM**2) EN = SQRT(Q2 + AM**2) BDX = PX/DM BDY = PY/DM BDZ = PZ/DM BP = (PX*QX + PY*QY + PZ*QZ)/(DM + EDELTA) E(I) = AM P(1,I) =-QX + BDX*(-BP + EN) P(2,I) =-QY + BDY*(-BP + EN) P(3,I) =-QZ + BDZ*(-BP + EN) PPION(1,NNN,IRUN) = QX + BDX*( BP + EP) PPION(2,NNN,IRUN) = QY + BDY*( BP + EP) PPION(3,NNN,IRUN) = QZ + BDZ*( BP + EP) RETURN END ************************************************************************ SUBROUTINE DRESONANCE(I1,I2) * PURPOSE : CALCULATE THE MASS AND MOMENTUM OF DELTA AFTER PION BEING * ABSORBED BY A NUCLEON * NOTE : THIS PROGRAM IS ONLY VALID FOR PION-NUCLEON DISTANCE SMALLER * THAN 2.11 FM * DATE : JAN.29,1990 ************************************************************************ IMPLICIT NONE INTEGER I1, I2, MAXPAR PARAMETER (MAXPAR=100000) INTEGER LB, ID REAL P, E COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) LOGICAL PION INTEGER I REAL E10, E20, DM * 1. DETERMINE THE MOMENTUM COMPONENT OF DELTA IN THE LAB. FRAME E10 = SQRT(E(I1)**2 + P(1,I1)**2 + P(2,I1)**2 + P(3,I1)**2) E20 = SQRT(E(I2)**2 + P(1,I2)**2 + P(2,I2)**2 + P(3,I2)**2) IF( pion(lb(i1)))THEN E(I1) = 0. I = I2 ELSE if ( pion(lb(i2)) ) then E(I2) = 0. I = I1 else stop 'there is no pion to be absorbed' ENDIF P(1,I) = P(1,I1) + P(1,I2) P(2,I) = P(2,I1) + P(2,I2) P(3,I) = P(3,I1) + P(3,I2) * 2. DETERMINE THE MASS OF DELTA BY USING OF THE REACTION KINEMATICS DM = SQRT((E10 + E20)**2 - P(1,I)**2 - P(2,I)**2 - P(3,I)**2) E(I) = DM RETURN END ************************************************************************ subroutine energy(t, v, vp, e0, fnorm, dbox) * * * calculate total energy * * * ************************************************************************ IMPLICIT NONE REAL T, V, VP, E0, FNORM, DBOX INTEGER MAXPAR, MAXX, MAXZ, MAXR INTEGER MX, MY, MZ, MPX, MPY, MPZ, MPZP PARAMETER (MAXPAR=100000, MAXX=40, MAXZ=81, MAXR=1000) PARAMETER (MX=4, MY=4, MZ=15, MPX=4, MPY=4, MPZ=5, MPZP=10) INTEGER NUM, MASS, ID, LB, MASSR REAL A, B, C, SIGMA, LAMBDA, RHO0 REAL R, P, E, RHO CHARACTER*1 FF COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /EE/ ID(MAXPAR), lb(maxpar) COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM COMMON /MASS/ MASS INTEGER ISO, IRUN, I, JSO, JRUN, J INTEGER ILX, ILY, ILZ, IX, IY, IZ, JX, JY, JZ REAL R1, R2, R3, L2, RM, DX t=0. e0=0. iso=0 do irun=1,num iso=iso+massr(irun-1) do i=1+iso,massr(irun)+iso t=t+sqrt(e(i)**2+p(1,i)**2+p(2,i)**2+p(3,i)**2) e0=e0+e(i) end do end do r1=0. r2=0. r3=0. do ix=-maxx,maxx do iy=-maxx,maxx do iz=-maxz,maxz r1=r1+rho(ix,iy,iz)**2. r2=r2+rho(ix,iy,iz)**(sigma+1.) end do end do end do r1=r1*dbox**3 r2=r2*dbox**3 c a and b are already divided by 2 ! v=a/rho0*r1+2.*b/(sigma+1.)/rho0**sigma*r2 dx=1.29 if ( c.ne.0 ) then l2=lambda**2 c loop over all boxes; Integral dr^3 do ilx=-9,9 do ily=-9,9 do ilz=-16,16 c loop i: all particles, is particle i in box? iso=0 do irun=1,num iso=iso+massr(irun-1) do i=1+iso,mass+iso ix=nint(r(1,i)/dx) if ( ix.eq.ilx ) then iy=nint(r(2,i)/dx) if ( iy.eq.ily ) then iz=nint(r(3,i)/dx) if ( iz.eq.ilz ) then c loop j: all particles, is particle j in box? jso=0 do jrun=1,num jso=jso+massr(jrun-1) do j=1+jso,mass+jso jx=nint(r(1,j)/dx) if ( ix.eq.jx ) then jy=nint(r(2,j)/dx) if ( iy.eq.jy ) then jz=nint(r(3,j)/dx) if ( iz.eq.jz ) then r3=r3+l2/(l2+(p(1,i)-p(1,j))**2 & +(p(2,i)-p(2,j))**2+(p(3,i)-p(3,j))**2) end if end if end if end do end do end if end if end if end do end do end do end do end do end if vp=c/rho0*r3/dx**3/num/num t=t/num e0=e0/num return end ************************************************************************ * * SUBROUTINE FRONT(IO, MASSTA, MASSPR, ELAB) * * * PURPOSE: WRITING FRONT PAGE OF OUTPUT * * VARIABLES: IO - OUTPUT-UNIT (INTEGER,INPUT) * * MASSTA - TARGET MASS (INTEGER,INPUT) * * MASSPR - PROJECTILE MASS (INTEGER,INPUT) * * ELAB - BEAM ENERGY "MEV/A" (REAL,INPUT) * * * ************************************************************************ IMPLICIT NONE INTEGER IO, MASSPR, MASSTA REAL ELAB WRITE(IO,'(''1'')') WRITE(IO,'(///27X,''BBBBBBBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBBBBBBBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBB BBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBB BBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBBBBBBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBBBBBBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBB BBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBB BBB UUU UUU UUU UUU'')') WRITE(IO,'(27X,''BBBBBBBBB UUUUUUUU UUUUUUUU'')') WRITE(IO,'(27X,''BBBBBBBB UUUUUU UUUUUU'')') WRITE(IO,'(27X,'' ''//)') WRITE(IO,'(10X,67(''*'')/10X,''*'',65X,''*''/ & 10X,''*'',10X,''REACTION:'',46X,''*''/ & 10X,''*'',65X,''*''/10X,''*'',22X,''MASS'',I4, & '' ==> MASS'',I4,22X,''*''/ & 10X,''*'',65X,''*''/10X,''*'',15X, & ''BEAM ENERGY = '',F9.2,'' MEV/NUCLEON'',15X,''*''/ & 10X,''*'',65X,''*''/10X,67(''*''))') MASSPR, & MASSTA, ELAB RETURN END ************************************************************************ SUBROUTINE GRADIN(IPOT, ICOLL) * * * provide parameters for potential * * * ************************************************************************ IMPLICIT NONE INTEGER IPOT, ICOLL, N PARAMETER (N=5) REAL A, B, C, SIGMA, LAMBDA, RHO0 COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 REAL AD(N), BD(N), CD(N), SIGD(N), LAMD(N), RHO0D(N) c value for A and B are already divided by 2 DATA AD /-0.062 , -0.109 , -0.178 , -0.05522, -0.002945/ DATA BD / 0.03525, 0.082 , 0.1515 , 0.07045, 0.01811/ DATA CD / 0.0 , 0.0 , 0.0 , -0.06495, -0.06491/ DATA SIGD / 2.0 , 1.333333, 1.166667, 1.24 , 2.45/ DATA LAMD / 0.0 , 0.0 , 0.0 , 0.41561, 0.41561/ DATA RHO0D / 0.168 , 0.168 , 0.168 , 0.16 , 0.16/ IF ( IPOT.LT.1.OR.IPOT.GT.N ) STOP 'POTENTIAL UNKNOWN' A = AD(IPOT) B = BD(IPOT) C = CD(IPOT) SIGMA = SIGD(IPOT) LAMBDA = LAMD(IPOT) RHO0 = RHO0D(IPOT) IF ( ICOLL.EQ.-1 ) then A = 0 B = 0 C = 0 end if RETURN END ************************************************************************ SUBROUTINE GRADU1(LF, INSYS) * * * PURPOSE: DETERMINE the momentum dependent part of * * GRADP/R(U(RHO(X,Y,Z))) only and store it * * in the arrays GRADX1 ... * * For momentum dependent potential U * * VARIABLES: * * RX, RY, RZ, PX, PY, PZ - COORDINATES OF POINT * * IN PHASESPACE (INTEGER,INPUT) * * GRADX, GRADY, GRADZ - GRADIENT OF U (REAL,OUTPUT) * * * ************************************************************************ IMPLICIT NONE INTEGER LF, INSYS INTEGER MAXX, MAXZ, MAXPAR, MAXR INTEGER MX, MY, MZ, MPX, MPY, MPZ, MPZP PARAMETER (MAXPAR=100000, MAXR=1000, MAXX=40, MAXZ=81) PARAMETER (MX=4, MY=4, MZ=15, MPX=4, MPY=4, MPZ=5, MPZP=10) CHARACTER*1 FF INTEGER NUM, MASS, MASSR REAL DX, DY, DZ, DPX, DPY, DPZ REAL A, B, C, SIGMA, LAMBDA, RHO0 REAL GRADX, GRADY, GRADZ, GRADPX, GRADPY, GRADPZ REAL R, P COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /GG/ DX, DY, DZ, DPX, DPY, DPZ COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) COMMON /II/ GRADX(MAXPAR), GRADY(MAXPAR), GRADZ(MAXPAR), 1 GRADPX(MAXPAR), GRADPY(MAXPAR), GRADPZ(MAXPAR) COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM COMMON /MASS/ MASS COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 integer ix, iy, iz, ipx, ipy, ipz, i, j, mpze, iso, nrun real px, py, pz, dpnx, dpny, dpn, cc, l2 real ddpx(-mpx:mpx), ddpy(-mpy:mpy), ddpz(-mpz:mpzp) real ddpx2(-mpx:mpx), ddpy2(-mpy:mpy), ddpz2(-mpz:mpzp) *----------------------------------------------------------------------- * * POTENTIAL USED: * U = 2*a * RHO/RHO0 + 2*b (RHO/RHO0)**(sigma) GEV * +momentum dependent part * Momentum independent part is calculated in gradu2 * ----------------------------------------------------------------------- * Momentum dependent part of the interaction if requested * see C.Gale et.al. Phys.Rev.C Vol.41 1545 (1990) * lf = 0 c loop over all baryonic test particles ISO = 0 L2 = LAMBDA*LAMBDA IF ( INSYS.NE.0 ) THEN MPZE = MPZ ELSE MPZE = MPZP END IF DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO J = MASS+ISO+1, massr(nrun) gradx(j) = 0 grady(j) = 0 gradz(j) = 0 gradpx(j) = 0 gradpy(j) = 0 gradpz(j) = 0 END DO IF (ABS(C).le.1E-5 ) THEN DO J = 1+ISO, MASS+ISO gradx(j) = 0. grady(j) = 0. gradz(j) = 0. gradpx(j) = 0. gradpy(j) = 0. gradpz(j) = 0. end do else DO J = 1+ISO, MASS+ISO ix = nint(r(1,j)/dx) iy = nint(r(2,j)/dy) iz = nint(r(3,j)/dz) px = p(1,j) py = p(2,j) pz = p(3,j) if ( abs(ix)+1.gt.mx .or. abs(iy)+1.gt.my .or. 1 abs(iz+1).gt.mz ) then c particle left the f-array lf = lf+1 else c take most of the multiplications out of the integral-loop ddpx(-mpx) = px+dpx*mpx ddpy(-mpy) = py+dpy*mpy ddpx2(-mpx) = ddpx(-mpx)**2 ddpy2(-mpy) = ddpy(-mpy)**2 do i = -mpx+1,mpx ddpx(i) = ddpx(i-1)-dpx ddpy(i) = ddpy(i-1)-dpy ddpx2(i) = ddpx(i)**2 ddpy2(i) = ddpy(i)**2 end do ddpz(-mpz) = pz+dpz*mpz ddpz2(-mpz) = ddpz(-mpz)**2 do i = -mpz+1,mpze ddpz(i) = ddpz(i-1)-dpz ddpz2(i) = ddpz(i)**2 end do c the dp^3 integral do ipx = -mpx,mpx dpnx = l2+ddpx2(ipx) do ipy = -mpy,mpy dpny = dpnx+ddpy2(ipy) do ipz = -mpz,mpze c grad_r(U) dpn = 1/(dpny+ddpz2(ipz)) gradx(j) = gradx(j) + 1 real(ichar(ff(ipz,ipy,ipx,ix+1,iy,iz)) - 1 ichar(ff(ipz,ipy,ipx,ix-1,iy,iz)))*dpn grady(j) = grady(j) + 1 real(ichar(ff(ipz,ipy,ipx,ix,iy+1,iz)) - 1 ichar(ff(ipz,ipy,ipx,ix,iy-1,iz)))*dpn gradz(j) = gradz(j) + 1 real(ichar(ff(ipz,ipy,ipx,ix,iy,iz+1)) - 1 ichar(ff(ipz,ipy,ipx,ix,iy,iz-1)))*dpn c grad_p(U) dpn = real(ichar(ff(ipz,ipy,ipx,ix,iy,iz)))*dpn**2 gradpx(j) = gradpx(j) + dpn*ddpx(ipx) gradpy(j) = gradpy(j) + dpn*ddpy(ipy) gradpz(j) = gradpz(j) + dpn*ddpz(ipz) end do end do end do end if c take care of the normalisationfactors cc = c/rho0*l2/dx/dy/dz/num gradx(j) = cc/dx*gradx(j) grady(j) = cc/dy*grady(j) gradz(j) = cc/dz*gradz(j) cc = -4.*cc gradpx(j) = cc*gradpx(j) gradpy(j) = cc*gradpy(j) gradpz(j) = cc*gradpz(j) end do end if end do RETURN END ************************************************************************ SUBROUTINE GRADU2(IX, IY, IZ, GRADX, GRADY, GRADZ, dbox) * * * PURPOSE: DETERMINE the density dependent, momentum * * independent part of GRAD(U(RHO(X,Y,Z))) * * for momentum dependent potential U * * VARIABLES: * * RX, RY, RZ, PX, PY, PZ - COORDINATES OF POINT * * IN PHASESPACE (INTEGER,INPUT) * * GRADX, GRADY, GRADZ - GRADIENT OF U (REAL,OUTPUT) * * * ************************************************************************ IMPLICIT NONE INTEGER IX, IY, IZ, MAXX, MAXZ REAL GRADX, GRADY, GRADZ, DBOX PARAMETER (MAXX=40, MAXZ=81) REAL A, B, C, SIGMA, LAMBDA, RHO0, RHO COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 REAL RXPLUS, RXMINS, RYPLUS, RYMINS, RZPLUS, RZMINS *----------------------------------------------------------------------- * * POTENTIAL USED: * U = 2*a * RHO/RHO0 + 2*b (RHO/RHO0)**(sigma) GEV * +momentum dependent part * Momentum independent part RXPLUS = RHO(IX+1,IY, IZ ) / RHO0 RXMINS = RHO(IX-1,IY, IZ ) / RHO0 RYPLUS = RHO(IX, IY+1,IZ ) / RHO0 RYMINS = RHO(IX, IY-1,IZ ) / RHO0 RZPLUS = RHO(IX, IY, IZ+1) / RHO0 RZMINS = RHO(IX, IY, IZ-1) / RHO0 GRADX = a * (RXPLUS - RXMINS) + b * (RXPLUS**sigma-RXMINS**sigma) GRADY = a * (RYPLUS - RYMINS) + b * (RYPLUS**sigma-RYMINS**sigma) GRADZ = a * (RZPLUS - RZMINS) + b * (RZPLUS**sigma-RZMINS**sigma) gradx = gradx/dbox grady = grady/dbox gradz = gradz/dbox *----------------------------------------------------------------------- * Momentum dependent part of the interaction if requested * see C.Gale et.al. Phys.Rev.C Vol.41 1545 (1990) * is already calculated in gradu1 and stored in an array RETURN END ************************************************************************ SUBROUTINE GRADPU(RX,RY,RZ,PX,PY,PZ,GRADPX,GRADPY,GRADPZ,lf,insys) * * * PURPOSE: DETERMINE GRADP(U(RHO(X,Y,Z))) * * For momentum dependent potential U * * to update the positions * * VARIABLES: * * RX, RY, RZ, PX, PY, PZ - COORDINATES OF POINT * * IN PHASESPACE (INTEGER,INPUT) * * GRADPX, GRADPY, GRADPZ - GRADIENT P OF U (REAL,OUTPUT) * * * ************************************************************************ IMPLICIT NONE INTEGER LF, INSYS REAL RX, RY, RZ, PX, PY, PZ, GRADPX, GRADPY, GRADPZ INTEGER MX, MY, MZ, MPX, MPY, MPZ, MPZP, MPZE PARAMETER (MX=4, MY=4, MZ=15, MPX=4, MPY=4, MPZ=5, MPZP=10) INTEGER NUM REAL A, B, C, SIGMA, LAMBDA, RHO0 REAL DX, DY, DZ, DPX, DPY, DPZ CHARACTER*1 FF COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 COMMON /RUN/ NUM COMMON /GG/ DX, DY, DZ, DPX, DPY, DPZ COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) INTEGER IX, IY, IZ, IPX, IPY, IPZ, I REAL DPN, CC, L2, DPNX, DPNY REAL DDPX(-MPX:MPX), DDPY(-MPY:MPY), DDPZ(-MPZ:MPZP) REAL DDPX2(-MPX:MPX), DDPY2(-MPY:MPY), DDPZ2(-MPZ:MPZP) lf=0.0 *----------------------------------------------------------------------- * * POTENTIAL USED: * U = 2*a * RHO/RHO0 + 2*b (RHO/RHO0)**(sigma) GEV * +momentum dependent part gradpx = 0. gradpy = 0. gradpz = 0. *----------------------------------------------------------------------- * Momentum dependent part of the interaction if requested * see C.Gale et.al. Phys.Rev.C Vol.41 1545 (1990) * if (abs(c).le.1e-5 ) return ix = nint(rx/dx) iy = nint(ry/dy) iz = nint(rz/dz) if ( abs(ix).gt.mx .or. abs(iy).gt.my .or. abs(iz).gt.mz ) then lf = 1. else if ( insys.ne.0 ) then mpze = mpz ELSE mpze = mpzp end if l2 = lambda*lambda ddpx(-mpx) = px+dpx*mpx ddpy(-mpy) = py+dpy*mpy ddpx2(-mpx) = ddpx(-mpx)**2 ddpy2(-mpy) = ddpy(-mpy)**2 do i = -mpx+1, mpx ddpx(i) = ddpx(i-1)-dpx ddpy(i) = ddpy(i-1)-dpy ddpx2(i) = ddpx(i)**2 ddpy2(i) = ddpy(i)**2 end do ddpz(-mpz) = pz+dpz*mpz ddpz2(-mpz) = ddpz(-mpz)**2 do i = -mpz+1,mpze ddpz(i) = ddpz(i-1)-dpz ddpz2(i) = ddpz(i)**2 end do do ipx = -mpx,mpx dpnx = l2+ddpx2(ipx) do ipy = -mpy,mpy dpny = dpnx+ddpy2(ipy) do ipz = -mpz,mpze dpn = real(ichar(ff(ipz,ipy,ipx,ix,iy,iz))) 1 /(dpny + ddpz2(ipz))**2 gradpx = gradpx+dpn*ddpx(ipx) gradpy = gradpy+dpn*ddpy(ipy) gradpz = gradpz+dpn*ddpz(ipz) end do end do end do cc = -4.*c/rho0*l2/dx/dy/dz/num gradpx = cc*gradpx gradpy = cc*gradpy gradpz = cc*gradpz end if RETURN END ************************************************************************ SUBROUTINE GRADRU(RX, RY, RZ, PX, PY, PZ, GRADX, GRADY, GRADZ, lf, 1 insys, dbox) * * * PURPOSE: DETERMINE GRADR(U(RHO(X,Y,Z))) * * For momentum dependent potential U * * VARIABLES: * * RX, RY, RZ, PX, PY, PZ - COORDINATES OF POINT * * IN PHASESPACE (INTEGER,INPUT) * * GRADX, GRADY, GRADZ - GRADIENT OF U (REAL,OUTPUT) * *!!!THIS ROUTINE IS CURRENTLY NOT USED!!!! * ************************************************************************ IMPLICIT NONE INTEGER INSYS, LF REAL RX, RY, RZ, PX, PY, PZ, GRADX, GRADY, GRADZ, DBOX INTEGER MAXX, MAXZ, MX, MY, MZ, MPX, MPY, MPZ, MPZP, MPZE PARAMETER (MAXX=40, MAXZ=81) PARAMETER (MX=4, MY=4, MZ=15, MPX=4, MPY=4, MPZ=5, MPZP=10) INTEGER NUM REAL A, B, C, SIGMA, LAMBDA, RHO0 REAL DX, DY, DZ, DPX, DPY, DPZ, RHO CHARACTER*1 FF COMMON /POTPARA/ A, B, C, SIGMA, LAMBDA, RHO0 COMMON /RUN/ NUM COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /GG/ DX, DY, DZ, DPX, DPY, DPZ COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) INTEGER IX, IY, IZ, IPX, IPY, IPZ, I REAL RXPLUS,RXMINS,RYPLUS,RYMINS,RZPLUS,RZMINS REAL GRADXI, GRADYI, GRADZI, DPN, CC,L2, DPNX, DPNY REAL DDPX(-MPX:MPX), DDPY(-MPY:MPY), DDPZ(-MPZ:MPZP) REAL DDPX2(-MPX:MPX), DDPY2(-MPY:MPY), DDPZ2(-MPZ:MPZP) lf=0 *----------------------------------------------------------------------- * * POTENTIAL USED: * U = 2*a * RHO/RHO0 + 2*b (RHO/RHO0)**(sigma) GEV * +momentum dependent part * Momentum independent part IX = NINT(RX/DBOX) IY = NINT(RY/DBOX) IZ = NINT(RZ/DBOX) RXPLUS = RHO(IX+1,IY, IZ ) / RHO0 RXMINS = RHO(IX-1,IY, IZ ) / RHO0 RYPLUS = RHO(IX, IY+1,IZ ) / RHO0 RYMINS = RHO(IX, IY-1,IZ ) / RHO0 RZPLUS = RHO(IX, IY, IZ+1) / RHO0 RZMINS = RHO(IX, IY, IZ-1) / RHO0 GRADX = a * (RXPLUS - RXMINS) + b * (RXPLUS**sigma-RXMINS**sigma) GRADY = a * (RYPLUS - RYMINS) + b * (RYPLUS**sigma-RYMINS**sigma) GRADZ = a * (RZPLUS - RZMINS) + b * (RZPLUS**sigma-RZMINS**sigma) GRADX = GRADX/DBOX GRADY = GRADY/DBOX GRADZ = GRADZ/DBOX *----------------------------------------------------------------------- * Momentum dependent part of the interaction if requested * see C.Gale et.al. Phys.Rev.C Vol.41 1545 (1990) * IF (ABS(C).LE.1E-5 ) RETURN IX = NINT(RX/DX) IY = NINT(RY/DY) IZ = NINT(RZ/DZ) IF(ABS(IX).GE.MX .OR. ABS(IY).GE.MY .OR. ABS(IZ).GE.MZ) THEN LF = 1 ELSE IF ( INSYS.NE.0 ) THEN MPZE = MPZ ELSE MPZE = MPZP END IF GRADXI = 0. GRADYI = 0. GRADZI = 0. L2 = LAMBDA*LAMBDA DDPX(-MPX) = PX+DPX*MPX DDPY(-MPY) = PY+DPY*MPY DDPX2(-MPX) = DDPX(-MPX)**2 DDPY2(-MPY) = DDPY(-MPY)**2 DO I = -MPX+1, MPX DDPX(I) = DDPX(I-1)-DPX DDPY(I) = DDPY(I-1)-DPY DDPX2(I) = DDPX(I)**2 DDPY2(I) = DDPY(I)**2 END DO ddpz(-mpz) = pz+dpz*mpz ddpz2(-mpz) = ddpz(-mpz)**2 do i = -mpz+1, mpze ddpz(i) = ddpz(i-1)-dpz ddpz2(i) = ddpz(i)**2 end do do ipx = -mpx, mpx dpnx = l2 + ddpx2(ipx) do ipy = -mpy, mpy dpny = dpnx + ddpy2(ipy) do ipz = -mpz, mpze dpn = dpny + ddpz2(ipz) gradxi = gradxi + real(ichar(ff(ipz,ipy,ipx,ix+1,iy,iz)) 1 -ichar(ff(ipz,ipy,ipx,ix-1,iy,iz)))/dpn gradyi = gradyi + real(ichar(ff(ipz,ipy,ipx,ix,iy+1,iz)) 1 -ichar(ff(ipz,ipy,ipx,ix,iy-1,iz)))/dpn gradzi = gradzi + real(ichar(ff(ipz,ipy,ipx,ix,iy,iz+1)) 1 -ichar(ff(ipz,ipy,ipx,ix,iy,iz-1)))/dpn end do end do end do cc = c/rho0*l2/dx/dy/dz/num c 1.90598=h**3 in GeVfm/c; 4=spin-isospin c cc=c/rho0*l2/dx/dy/dz/num*1.90598/4. gradx = gradx + cc/dx*gradxi grady = grady + cc/dy*gradyi gradz = gradz + cc/dz*gradzi end if RETURN END ************************************************************************ SUBROUTINE INIT(MINNUM, MAXNUM, RADIUS, X0, Z0, P0, GAMMA, ISEED, 1 massnu, znu, IOPT, ioptc) * * * PURPOSE: PROVIDING INITIAL CONDITIONS FOR PHASE-SPACE * * DISTRIBUTION OF TESTPARTICLES * * VARIABLES: (ALL INPUT) * * MINNUM - FIRST TESTPARTICLE TREATED IN ONE RUN (INTEGER) * * MAXNUM - LAST TESTPARTICLE TREATED IN ONE RUN (INTEGER) * * NUM - NUMBER OF TESTPARTICLES PER NUCLEON (INTEGER) * * RADIUS - RADIUS OF NUCLEUS "FM" (REAL) * * X0,Z0 - DISPLACEMENT OF CENTER OF NUCLEUS IN X,Z- * * DIRECTION "FM" (REAL) * * P0 - MOMENTUM-BOOST IN C.M. FRAME "GEV/C" (REAL) * * GAMMA - RELATIVISTIC GAMMA-FACTOR (REAL) * * ISEED - SEED FOR RANDOM-NUMBER GENERATOR (INTEGER) * * MASS - TOTAL MASS OF THE SYSTEM (INTEGER) * * MASSNU - MASS OF THE PARTICLE (INTEGER) * * ZNU - CHARGE OF THE PARTICLE (INTEGER) * * IOPT - OPTION FOR DIFFERENT OCCUPATION OF MOMENTUM * * SPACE (INTEGER) * * IOPTC - OPTION FOR DIFFERENT OCCUPATION OF COORDINATE * * SPACE (INTEGER) * * * ************************************************************************ IMPLICIT NONE INTEGER MINNUM, MAXNUM, ISEED, MASSNU, ZNU, IOPT, IOPTC REAL RADIUS, X0, Z0, P0, GAMMA INTEGER MAXPAR, MAXX, MAXZ REAL AMU, PI, HC PARAMETER (MAXPAR=100000, MAXX=40, MAXZ=81) PARAMETER (PI=3.14159, HC=0.197327, AMU=0.938919) INTEGER NUM, MASS, ID, LB REAL E, R, P COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /MASS/ MASS COMMON /RUN/ NUM * functions called REAL NURHO * local variables INTEGER I, IDIR, IDNUM, IRUN, NPART REAL BETA, EPART, PFERMI, PX, PY, PZ REAL RDIST, RHOW0, RHOWS, SIGN, X, Y, Z REAL RL, R0, MAXRHO, A(17), PTOT(3) *----------------------------------------------------------------------- * PREPARATION FOR LORENTZ-TRANSFORMATIONS * IF (P0 .NE. 0.) THEN SIGN = P0 / ABS(P0) ELSE SIGN = 0. END IF BETA = SIGN * SQRT(GAMMA**2-1.)/GAMMA *----------------------------------------------------------------------- * TARGET-ID = 1 AND PROJECTILE-ID = -1 * IF (MINNUM .EQ. 1) THEN IDNUM = 1 ELSE IDNUM = -1 END IF if ( ioptc.eq.2 ) then call initnu(massnu, znu, a, r0) maxrho=0. do rl=0.,r0,.01 if ( nurho(rl,a,r0).gt.maxrho ) then maxrho=nurho(rl,a,r0) end if end do end if *----------------------------------------------------------------------- * IDENTIFICATION OF TESTPARTICLES AND ASSIGMENT OF RESTMASS * * LOOP OVER ALL PARALLEL RUNS: DO 400 IRUN = 1,NUM DO I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS ID(I) = IDNUM E(I) = AMU end do *----------------------------------------------------------------------- * OCCUPATION OF COORDINATE-SPACE * DO 300 I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS if ( ioptc.eq.1 ) then 200 CONTINUE X = 1.0 - 2.0 * RAN(ISEED) Y = 1.0 - 2.0 * RAN(ISEED) Z = 1.0 - 2.0 * RAN(ISEED) IF (X*X+Y*Y+Z*Z .GT. 1.0) GOTO 200 R(1,I) = X * RADIUS R(2,I) = Y * RADIUS R(3,I) = Z * RADIUS else 250 X = r0*(1.0 - 2.0 * RAN(ISEED)) Y = r0*(1.0 - 2.0 * RAN(ISEED)) Z = r0*(1.0 - 2.0 * RAN(ISEED)) RL= sqrt(X*X+Y*Y+Z*Z) if (rl.gt.r0) goto 250 IF ( ran(iseed).gt.nurho(rl,a,r0)/maxrho) GOTO 250 R(1,I) = X R(2,I) = Y R(3,I) = Z end if 300 CONTINUE 400 CONTINUE *======================================================================= IF (IOPT .NE. 3) THEN *----- * OPTION 1: USE WOODS-SAXON PARAMETRIZATION FOR DENSITY AND *----- CALCULATE LOCAL FERMI-MOMENTUM * RHOW0 = 0.168 DO 1000 IRUN = 1,NUM DO 600 I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS 500 CONTINUE PX = 1.0 - 2.0 * RAN(ISEED) PY = 1.0 - 2.0 * RAN(ISEED) PZ = 1.0 - 2.0 * RAN(ISEED) IF (PX*PX+PY*PY+PZ*PZ .GT. 1.0) GOTO 500 RDIST = SQRT( R(1,I)**2 + R(2,I)**2 + R(3,I)**2 ) RHOWS = RHOW0 / ( 1.0 + EXP( (RDIST-RADIUS) / 0.55 ) ) PFERMI = hc * (1.5 * PI*PI * RHOWS)**(1./3.) *----- * OPTION 2: NUCLEAR MATTER CASE IF(IOPT.EQ.2) PFERMI=0.27 *----- c proton doesn't have fermi motion if ( (maxnum-minnum).eq.0 ) pfermi=0. P(1,I) = PFERMI * PX P(2,I) = PFERMI * PY P(3,I) = PFERMI * PZ 600 CONTINUE * * SET TOTAL MOMENTUM TO 0 IN REST FRAME AND BOOST * DO 700 IDIR = 1,3 PTOT(IDIR) = 0.0 700 CONTINUE NPART = 0 DO 900 I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS NPART = NPART + 1 DO 800 IDIR = 1,3 PTOT(IDIR) = PTOT(IDIR) + P(IDIR,I) 800 CONTINUE 900 CONTINUE DO 950 I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS DO 925 IDIR = 1,3 P(IDIR,I) = P(IDIR,I) - PTOT(IDIR) / FLOAT(NPART) 925 CONTINUE * BOOST IF (IOPT .EQ. 1) THEN EPART = SQRT(P(1,I)**2+P(2,I)**2+P(3,I)**2+AMU**2) P(3,I) = GAMMA*(P(3,I) + BETA*EPART) ELSE P(3,I) = P(3,I) + P0 END IF 950 CONTINUE 1000 CONTINUE *----- ELSE *----- * OPTION 3: GIVE ALL NUCLEONS JUST A Z-MOMENTUM ACCORDING TO * THE BOOST OF THE NUCLEI * DO 1200 IRUN = 1,NUM DO 1100 I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS P(1,I) = 0.0 P(2,I) = 0.0 P(3,I) = P0 1100 CONTINUE 1200 CONTINUE *----- END IF *======================================================================= * PUT PARTICLES IN THEIR POSITION IN COORDINATE-SPACE * (SHIFT AND RELATIVISTIC CONTRACTION) * DO 1400 IRUN = 1,NUM DO 1300 I = MINNUM+(IRUN-1)*MASS,MAXNUM+(IRUN-1)*MASS R(1,I) = R(1,I) + X0 R(3,I) = R(3,I) / GAMMA + Z0 1300 CONTINUE 1400 CONTINUE * RETURN END ************************************************************************ subroutine initnu(mass, z, a, r0) * * * input: massnumber and charge * * output: Fourier-Bessel C. and radius * * * ************************************************************************ IMPLICIT NONE INTEGER MASS, Z REAL A(17), R0 INTEGER N PARAMETER (N=1) integer i, j, masss(n), zs(n) real as(17,n), r0s(n) data masss /12/ data zs /6/ data r0s /8./ data as/ 0.15721e-1, 0.38732e-1, 0.36808e-1, 0.14671e-1, 1 -0.43277e-2, -0.97752e-2, -0.68908e-2, -0.27631e-2, 1 -0.63568e-3, 0.71808e-4, 0.18441e-3, 0.75066e-4, 1 0.51069e-4, 0.14308e-4, 0.23170e-5, 0.68465e-6, 0.0/ r0 = 0. do i = 1, n if ( mass.eq.masss(i).and.z.eq.zs(i) ) then r0 = r0s(i) do j = 1,17 a(j) = as(j,i) end do end if end do if ( r0.eq.0. ) then write (*,*) 'No data for nucleus',mass,z stop 'mz' end if return end ************************************************************************ SUBROUTINE MATRI2(IO, FIELD, DIM1, DIM2, TITLE) * * * PURPOSE: THIS SUBROUTINE GIVES A PLOT-LIKE MATRIX TO THE PRINTER * * IN WHICH THE FIELD-ELEMENTS ARE REPRESENTED BY 1-DIGIT * * INTEGER-NUMBERS * * INPUT: IO - UNIT-NUMBER FOR OUTPUT (INTEGER) * * FIELD - 2-DIMENSIONAL ARRAY TO BE PLOTTED (REAL) * * DIM1 - NUMBER OF POINTS IN X-DIRECTION (INTEGER) * * DIM2 - NUMBER OF POINTS IN Y-DIRECTION (INTEGER) * * TITLE - HEADER OF THE PLOT (CHARACTER * 60) * * * ************************************************************************ CHARACTER * 2 ZEILE(60) CHARACTER *60 TITLE, OUTFOR REAL FIELD(60,90) INTEGER DIM1, DIM2 * * /* 1.: FIND MAXIMUM OF FIELD HIGH = 0.0 DO 200 J = 1,DIM2 DO 100 I = 1,DIM1 IF(ABS(FIELD(I,J)).GT.HIGH) HIGH = ABS(FIELD(I,J)) 100 CONTINUE 200 CONTINUE * /* ERROR-CONDITION-EXIT IF(HIGH.EQ.0.0) THEN WRITE(IO,'(///'' == E ===> ----- SUBROUTINE MATRI2 -----'')') WRITE(IO,'('' == E ===> ATTEMPTED PLOT OF:'')') WRITE(IO,'('' == E ===> '',A60)') TITLE WRITE(IO,'('' == E ===> ROUTINE ENDED ABNORMALLY DUE TO '', & ''ERROR, PLOT-FIELD EMPTY''///)') RETURN END IF * /* 2.: WRITE HEADLINES WRITE(IO,'('' => ** SUBROUTINE MATRI2 ** (1-DIGIT)'',T90,''<='')') WRITE(IO,'('' => PLOT OF: '',A60,T90,''<='')') TITLE WRITE(IO,'('' => NUMBER OF GRID-POINTS: '',I2,'' * '',I2, 1 T90,''<='')') DIM1, DIM2 WRITE(IO,'(T10,''***** CONTOUR-HEIGHT-LEVELS: *****'')') WRITE(IO,'(T10,''9 : '',G13.6,'' - '',G13.6)') HIGH, HIGH*.9 WRITE(IO,'(T10,''7 : '',G13.6,'' - '',G13.6)') HIGH*.8, HIGH*.7 WRITE(IO,'(T10,''5 : '',G13.6,'' - '',G13.6)') HIGH*0.6, HIGH*0.5 WRITE(IO,'(T10,''3 : '',G13.6,'' - '',G13.6)') HIGH*0.4, HIGH*0.3 WRITE(IO,'(T10,''1 : '',G13.6,'' - '',G13.6)') HIGH*0.2, HIGH*0.1 WRITE(IO,'(T10,''- : > 0'')') * /* 3.: PLOT MATRIX WRITE(IO,'(''1'')') WRITE(OUTFOR,605) (DIM1+2)*2 605 FORMAT('(4X,',I3,'(''-''))') WRITE(IO,OUTFOR) * WRITE(OUTFOR,610) DIM1 610 FORMAT('(1X,I2,'' I '',',I2,'A2,'' I'')') DO 700 J = DIM2,1,-1 DO 650 I = 1,DIM1 HELP = FIELD(I,J) WRITE(ZEILE(I),'(I2.0)') INT(HELP*9.999/HIGH) IF (HELP .LT. HIGH*0.1 .AND. HELP .GT. 0.0001*HIGH) THEN ZEILE(I) = ' -' END IF 650 CONTINUE WRITE(IO,OUTFOR) J,(ZEILE(I),I=1,DIM1) 700 CONTINUE * WRITE(OUTFOR,710) (DIM1+2)*2, DIM1/2 710 FORMAT('(4X,',I3,'(''-'')/3X,'' '',',I2,'I4,'' '')') WRITE(IO,OUTFOR) (I,I=2,DIM1,2) * RETURN END ************************************************************************ subroutine ndnnch(i1, i2, srt, sig, signn, x1, nstar, m12, renorm) * * * purpose: determine channel for N+Resonance->N+N reactions * * channel is given in M12 (output) * * input: i1, i2: particle number * * srt : CM energy * * ndirect: Flag to allow direct procesess (percentage dir)* ************************************************************************ IMPLICIT NONE INTEGER I1, I2, NSTAR, M12 REAL SRT, SIG, SIGNN, X1, RENORM INTEGER MAXPAR PARAMETER (MAXPAR=100000) INTEGER PROTON, NEUTRON, PIONM, PION0, PIONP, DELTAM INTEGER DELTA0, DELTAP, DELTAPP, NS0, NSP, ID, LB REAL MASSPA(11) COMMON /PARTICLE/ MASSPA, PROTON, NEUTRON, PIONM, PION0, PIONP, 1 DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP COMMON /EE/ ID(MAXPAR), LB(MAXPAR) LOGICAL DELT REAL SIGND, SIGDN c functions REAL SIGMA, DENOM m12 = 0 * FOR n+delta(++)-->P+P: if ( (lb(i1).eq.neutron.and.lb(i2).eq.deltapp) .or. & (lb(i2).eq.neutron.and.lb(i1).eq.deltapp)) THEN SIGND = SIGMA(3,1,SRT) + 0.5*SIGMA(2,1,SRT) DELT = .TRUE. SIGDN = 0.25*SIGND*RENORM/DENOM(SRT,DELT) IF (X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 3 * FOR p+delta(-)-->n+n ELSE IF ((lb(i1).eq.proton .and.lb(i2).eq.deltam ) .or. 1 (lb(i2).eq.proton .and.lb(i1).eq.deltam ) ) THEN SIGND = SIGMA(3,1,SRT) + 0.5*SIGMA(2,1,SRT) DELT = .TRUE. SIGDN = 0.25*SIGND*RENORM/dENOM(SRT,DELT) IF (X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 6 * FOR p+delta(+)-->p+p else if ( (lb(i1).eq.proton .and.lb(i2).eq.deltap ) .or. & (lb(i2).eq.proton .and.lb(i1).eq.deltap ) ) then SIGND = 1.5*SIGMA(2,1,SRT) DELT = .TRUE. SIGDN = 0.25*SIGND*RENORM/DENOM(SRT,DELT) IF(X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 4 * FOR n+delta(0)-->n+n ELSE if ((lb(i1).eq.neutron.and.lb(i2).eq.delta0 ) .or. & (lb(i2).eq.neutron.and.lb(i1).eq.delta0 ) ) THEN SIGND = 1.5*SIGMA(2,1,SRT) DELT = .TRUE. SIGDN = 0.25*SIGND*RENORM/DENOM(SRT,DELT) IF(X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 5 * FOR n+delta(+)-->n+p else if ( (lb(i1).eq.neutron.and.lb(i2).eq.deltap) .or. & (lb(i2).eq.neutron.and.lb(i1).eq.deltap)) THEN SIGND = 0.5*SIGMA(2,1,SRT) + 0.25*SIGMA(3,1,SRT) DELT = .TRUE. SIGDN = 0.5*SIGND*RENORM/DENOM(SRT,DELT) IF (X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 1 * p+delta(0)-->p+n ELSE if ((lb(i1).eq.proton .and.lb(i2).eq.delta0 ) .or. & (lb(i2).eq.proton .and.lb(i1).eq.delta0 ) ) THEN SIGND = 0.5*SIGMA(2,1,SRT) + 0.25*SIGMA(3,1,SRT) DELT = .TRUE. SIGDN = 0.5*SIGND*RENORM/DENOM(SRT,DELT) IF (X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 2 * FOR p+N*(0)-->p+n else if ( (lb(i1).eq.proton .and.lb(i2).eq.ns0) .or. & (lb(i2).eq.proton .and.lb(i1).eq.ns0) ) then IF (NSTAR.EQ.1) THEN SIGND = (3./4.)*SIGMA(4,2,SRT) DELT = .FALSE. SIGDN = SIGND*RENORM/DENOM(SRT,DELT) IF (X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 7 else stop 'NDNNCH: N* detected, but NSTAR=0!' end if * n+N*(+)-->n+p ELSE if ( (lb(i1).eq.neutron.and.lb(i2).eq.nsp) .or. 1 (lb(i2).eq.neutron.and.lb(i1).eq.nsp) ) THEN IF (NSTAR.EQ.1) THEN SIGND = (3./4.)*SIGMA(4,2,SRT) DELT = .FALSE. SIGDN = SIGND*RENORM/DENOM(SRT,DELT) IF (X1.GT.(SIGNN + SIGDN)/SIG) RETURN M12 = 8 else stop 'NDNNCH: N* detected, but NSTAR=0!' end if END IF return end ************************************************************************ subroutine nndnch(i1, i2, srt, sig, signn, x1, ndirect, dir, 1 nstar, dirfl, n12, iseed) * * * purpose: determine channel for N+N->N+Resonance reactions * * channel is given in N12 (output) * * dirfl=1 (output) is direct process happens * * input: i1, i2: particle number * * srt : CM energy * * ndirect: Flag to allow direct procesess (percentage dir)* * N12=1 FOR p+n-->delta(+)+ n * * 2 p+n-->delta(0)+ p * * 3 p+p-->delta(++)+n * * 4 p+p-->delta(+)+p * * 5 n+n-->delta(0)+n * * 6 n+n-->delta(-)+p * * 7 n+p-->N*(0)+p * * 8 n+p-->N*(+)+n * ************************************************************************ IMPLICIT NONE INTEGER I1, I2, NDIRECT, NSTAR, DIRFL, N12, ISEED REAL SRT, DIR, SIGNN, SIG INTEGER MAXPAR PARAMETER (MAXPAR=100000) INTEGER DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP INTEGER PROTON, NEUTRON, PIONM, PION0, PIONP, ID, LB REAL MASSPA(11) COMMON /PARTICLE/ MASSPA, PROTON, NEUTRON, PIONM, PION0, PIONP, 1 DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP COMMON /EE/ ID(MAXPAR), LB(MAXPAR) REAL X1, SIG1, SIG2, SIGND, X2 c functions REAL SIGMA dirfl = 1 x2 = ran(iseed) n12 = 0 IF (LB(I1) .eq. proton .and. LB(I2) .eq. proton) THEN C p+p collision SIG1 = SIGMA(3,1,SRT) + 0.5*SIGMA(2,1,SRT) SIG2 = 1.5*SIGMA(2,1,SRT) SIGND = SIG1 + SIG2 IF (X1 .GT. (SIGNN + SIGND)/SIG) RETURN IF ((NDIRECT.EQ.1) .and. (RAN(ISEED).LE.DIR)) RETURN IF (X2.GT.SIG1/SIGND) THEN N12 = 4 ELSE N12 = 3 END IF ELSE IF(LB(I1) .eq. neutron .and. LB(I2) .EQ. neutron) THEN C n+n collision SIG1 = SIGMA(3,1,SRT) + 0.5*SIGMA(2,1,SRT) SIG2 = 1.5*SIGMA(2,1,SRT) SIGND = SIG1 + SIG2 IF (X1 .GT. (SIGNN + SIGND)/SIG) RETURN IF ((NDIRECT.EQ.1) .and. (RAN(ISEED).LE.DIR)) RETURN IF (X2.GT.SIG1/SIGND) THEN N12 = 5 ELSE N12 = 6 END IF ELSE IF( (LB(I1) .eq. neutron .and. LB(I2) .EQ. proton) .or. & (LB(I2) .eq. neutron .and. LB(I1) .EQ. proton) ) THEN C n+p collision SIG1 = 0.5*SIGMA(2,1,SRT) + 0.25*SIGMA(3,1,SRT) IF(NSTAR.EQ.1)THEN SIG2 = (3./4.)*SIGMA(4,2,SRT) ELSE SIG2 = 0. ENDIF SIGND = 2.0*(SIG1 + SIG2) IF (X1 .GT. (SIGNN + SIGND)/SIG) RETURN IF ((NDIRECT.EQ.1) .and. (RAN(ISEED).LE.DIR)) RETURN IF (X2.LT.SIG1/(SIG1 + SIG2) ) THEN N12 = 2 IF (RAN(ISEED).GE.0.5) N12 = 1 ELSE IF (NSTAR.EQ.1) THEN N12 = 8 IF (RAN(ISEED).GE.0.5) N12 = 7 ENDIF ENDIF ELSE C error!! TYPE*, I1, I2 TYPE*, LB(I1), LB(I2) STOP 'NNDNCH: channel not available' ENDIF DIRFL = 0 RETURN END ************************************************************************ * * SUBROUTINE PAULat(I,occup) * * * this is the new version which treats Pauli-blocking on a lattice * * * * VARIABLES: * * I - NUMBER OF PARTICLE (INTEGER, INPUT) * * occup - occupation of phase space element (real,output) * * [0,1] is the interval in which occup should be * * * ************************************************************************ implicit none integer ix,iy,iz,ipx,ipy,ipz integer ixm,iym,izm,ipxm,ipym,ipzm integer ixp,iyp,izp,ipxp,ipyp,ipzp integer mx,my,mz,mpx,mpy,mpz,mpzp integer i,j,kk,maxpar,nruns,mode integer nruns_in,mode_in integer num, mass integer iso, nrun, massr, maxr, lb, id real x,y,z,px,py,pz real dx,dy,dz,dpx,dpy,dpz real x_in,y_in,z_in,px_in,py_in,pz_in real dx_in,dy_in,dz_in,dpx_in,dpy_in,dpz_in,fnorm_out real rx,ry,rz,rpx,rpy,rpz real amu,e,r,p,occup, fnorm real occup_max,occup_min,histo(-2:16),calls,sumocc,sum0 character*1 ff, ffzero PARAMETER (AMU=0.938919, MAXPAR=100000) parameter (maxr=1000) parameter (mx=4, my=4, mz=15, mpx=4, mpy=4, mpz=5, mpzp=10) COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), lb(maxpar) common /gg/ dx, dy, dz, dpx, dpy, dpz common /hh/ ff(-mpz:mpzp,-mpy:mpy,-mpx:mpx,-mx:mx,-my:my,-mz:mz) COMMON /RR/ MASSR(0:MAXR) common /run/ num common /mass/ mass common /paul/ calls * calls = calls + 1. IF (abs(E(i)-AMU) .GT. 0.01) then * This particle is not a nucleon: no Pauli blocking histo(0) = histo(0) + 1. sum0 = sum0 + 1. occup = 0. return end if * * get coordinates of particle: X = R(1,I) Y = R(2,I) Z = R(3,I) PX = P(1,I) PY = P(2,I) PZ = P(3,I) go to 1 *----------------------------------------------------------------------- * with this entry point one can also use any set of phase space * coordinates and determine its pauli-blocking factor: entry PAULI1(px_in,py_in,pz_in,x_in,y_in,z_in,occup) calls = calls + 1. x = x_in y = y_in z = z_in px = px_in py = py_in pz = pz_in * 1 ix = nint(x/dx) iy = nint(y/dy) iz = nint(z/dz) ipx = nint(px/dpx) ipy = nint(py/dpy) ipz = nint(pz/dpz) * if ( (abs(ix).gt.mx) .or. (abs(iy).gt.my) .or. (abs(iz).gt.mz) & .or. (abs(ipx).gt.mpx) .or. (abs(ipy).gt.mpy) & .or. (ipz.lt.-mpz) .or. (ipz.gt.mpzp)) then occup = 0. histo(-2) = histo(-2) + 1. else if ( mode.eq.3) then occup=0 else if (mode .eq. 2) then * local evaluation of f at nearest point occup = real(ichar(ff(ipz,ipy,ipx,ix,iy,iz))) * fnorm else if (mode .eq. 1) then * interpolation using the surrounding 2**6=64 points * 1) get neighboring points in all coordinates ixm = nint(x/dx-0.5) iym = nint(y/dy-0.5) izm = nint(z/dz-0.5) ipxm = nint(px/dpx-0.5) ipym = nint(py/dpy-0.5) ipzm = nint(pz/dpz-0.5) ixp = ixm + 1 iyp = iym + 1 izp = izm + 1 ipxp = ipxm + 1 ipyp = ipym + 1 ipzp = ipzm + 1 rx = x/dx - real(ixm) ry = y/dy - real(iym) rz = z/dz - real(izm) rpx = px/dpx - real(ipxm) rpy = py/dpy - real(ipym) rpz = pz/dpz - real(ipzm) * 2) take care of boundaries if (ixm .lt. -mx) then ixm = -mx rx = 1. end if if (iym .lt. -my) then iym = -my ry = 1. end if if (izm .lt. -mz) then izm = -mz rz = 1. end if if (ipxm .lt. -mpx) then ipxm = -mpx rpx = 1. end if if (ipym .lt. -mpy) then ipym = -mpy rpy = 1. end if if (ipzm .lt. -mpz) then ipzm = -mpz rpz = 1. end if if (ixp .gt. mx) then ixp = mx rx = 0. end if if (iyp .gt. my) then iyp = my ry = 0. end if if (izp .gt. mz) then izp = mz rz = 0. end if if (ipxp .gt. mpx) then ipxp = mpx rpx = 0. end if if (ipyp .gt. mpy) then ipyp = mpy rpy = 0. end if if (ipzp .gt. mpzp) then ipzp = mpzp rpz = 0. end if * 3) interpolate occup = (1.-rx)*(1.-ry)*(1.-rz)*(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)*(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)*(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixm,iyp,izm))) & + rx * ry *(1.-rz)*(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz *(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixm,iym,izp))) & + rx *(1.-ry)* rz *(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixp,iym,izp))) & + (1.-rx)* ry * rz *(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixm,iyp,izp))) & + rx * ry * rz *(1.-rpx)*(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxm,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)* rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)* rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)* rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixm,iyp,izm))) & + rx * ry *(1.-rz)* rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz * rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixm,iym,izp))) & + rx *(1.-ry)* rz * rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixp,iym,izp))) & + (1.-rx)* ry * rz * rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixm,iyp,izp))) & + rx * ry * rz * rpx *(1.-rpy)*(1.-rpz)* & real(ichar(ff(ipzm,ipym,ipxp,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)*(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)*(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)*(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixm,iyp,izm))) & + rx * ry *(1.-rz)*(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz *(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixm,iym,izp))) & + rx *(1.-ry)* rz *(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixp,iym,izp))) & + (1.-rx)* ry * rz *(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixm,iyp,izp))) & + rx * ry * rz *(1.-rpx)* rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxm,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)* rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)* rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)* rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixm,iyp,izm))) & + rx * ry *(1.-rz)* rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz * rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixm,iym,izp))) & + rx *(1.-ry)* rz * rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixp,iym,izp))) & + (1.-rx)* ry * rz * rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixm,iyp,izp))) & + rx * ry * rz * rpx * rpy *(1.-rpz)* & real(ichar(ff(ipzm,ipyp,ipxp,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)*(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)*(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)*(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixm,iyp,izm))) & + rx * ry *(1.-rz)*(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz *(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixm,iym,izp))) & + rx *(1.-ry)* rz *(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixp,iym,izp))) & + (1.-rx)* ry * rz *(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixm,iyp,izp))) & + rx * ry * rz *(1.-rpx)*(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxm,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)* rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)* rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)* rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixm,iyp,izm))) & + rx * ry *(1.-rz)* rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz * rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixm,iym,izp))) & + rx *(1.-ry)* rz * rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixp,iym,izp))) & + (1.-rx)* ry * rz * rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixm,iyp,izp))) & + rx * ry * rz * rpx *(1.-rpy)* rpz * & real(ichar(ff(ipzp,ipym,ipxp,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)*(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)*(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)*(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixm,iyp,izm))) & + rx * ry *(1.-rz)*(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz *(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixm,iym,izp))) & + rx *(1.-ry)* rz *(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixp,iym,izp))) & + (1.-rx)* ry * rz *(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixm,iyp,izp))) & + rx * ry * rz *(1.-rpx)* rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxm,ixp,iyp,izp))) occup = occup & + (1.-rx)*(1.-ry)*(1.-rz)* rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixm,iym,izm))) & + rx *(1.-ry)*(1.-rz)* rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixp,iym,izm))) & + (1.-rx)* ry *(1.-rz)* rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixm,iyp,izm))) & + rx * ry *(1.-rz)* rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixp,iyp,izm))) occup = occup & + (1.-rx)*(1.-ry)* rz * rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixm,iym,izp))) & + rx *(1.-ry)* rz * rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixp,iym,izp))) & + (1.-rx)* ry * rz * rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixm,iyp,izp))) & + rx * ry * rz * rpx * rpy * rpz * & real(ichar(ff(ipzp,ipyp,ipxp,ixp,iyp,izp))) occup = occup * fnorm else * no pauli blocking occup = 0. end if end if sumocc = sumocc + occup if (occup .lt. occup_min) occup_min = occup if (occup .gt. occup_max) occup_max = occup if (occup .lt. 0.) then histo(-1) = histo(-1) + 1. else if (occup .gt. 1.55) then histo(16) = histo(16) + 1. else histo(nint(occup/0.1)) = histo(nint(occup/0.1)) + 1. end if * RETURN *----------------------------------------------------------------------- entry platout if (calls .gt. 0.) then Write(*,'(/'' Occupation statistics'')') Write(*,'( '' ---------------------'')') write(*,'(/'' occupation number'')') write(*,'( '' out '',f11.1)') histo(-2) write(*,'( '' < 0 '',f11.1)') histo(-1) do kk = 0,15 write(*,'(f13.2,f15.1)') real(kk)*0.1, histo(kk) end do write(*,'( '' >1.5 '',f11.1)') histo(16) write(*,'(/'' Maximum occupation ='',f13.4)') occup_max write(*,'( '' Minimum occupation ='',f13.4)') occup_min write(*,'( '' Average occupation ='',f13.4)') sumocc/calls write(*,'( '' No nucleon (=sum0) ='',f13.4)') sum0 write(*,'( '' fnorm ='',f13.4)') fnorm else write(*,'(/'' PAULAT: No calls yet !'')') end if return * *----------------------------------------------------------------------* ENTRY platin(mode_in,nruns_in, & dx_in,dy_in,dz_in,dpx_in,dpy_in,dpz_in,fnorm_out) * * * mode - =1 -> interpolation, =2 -> local * * other -> unblocked (integer,input) * * nruns - number of test particle per nucleon (integer,input) * * dx,dy,dz,dpx,dpy,dpz * * - spacing in 6 dimensions of phase space (real,input) * * units: dx (fm), dpx (GeV/c) * * fnorm - normalization factor for f, dimensionless, specifies * * how many nucleons should be in a fully occupied * * phase space cell (real,output) * *----------------------------------------------------------------------* * keep values of dummy variables: nruns = nruns_in mode = mode_in dx = dx_in dy = dy_in dz = dz_in dpx = dpx_in dpy = dpy_in dpz = dpz_in * ffzero = char(0) * clean ff-lattice: do iz = -mz,mz do iy = -my,my do ix = -mx,mx do ipx = -mpx,mpx do ipy = -mpy,mpy do ipz = -mpz,mpzp ff(ipz,ipy,ipx,ix,iy,iz) = ffzero end do end do end do end do end do end do * * occupy ff-lattice: iso=0 do nrun=1, num iso=iso+massr(nrun-1) do j=1+iso, mass+iso if ( lb(j).le.2 ) then ix = nint(r(1,j)/dx) if (abs(ix) .gt. mx) goto 1111 iy = nint(r(2,j)/dy) if (abs(iy) .gt. my) goto 1111 iz = nint(r(3,j)/dz) if (abs(iz) .gt. mz) goto 1111 ipx = nint(p(1,j)/dpx) if (abs(ipx) .gt. mpx) goto 1111 ipy = nint(p(2,j)/dpy) if (abs(ipy) .gt. mpy) goto 1111 ipz = nint(p(3,j)/dpz) if ((ipz.ge.-mpz) .and. (ipz.le.mpzp)) 1 ff(ipz,ipy,ipx,ix,iy,iz) = 1 char(ichar(ff(ipz,ipy,ipx,ix,iy,iz))+1) 1111 continue end if end do end do * * calculate f-norm (1.90598=h**3 in GeVfm/c, 4 for spin-isospin-deg.): fnorm = 1.90598 / (4.*real(nruns)*dx*dy*dz*dpx*dpy*dpz) fnorm_out = fnorm * RETURN END ************************************************************************ * * SUBROUTINE PDENS(P0, OUT, UNIT, FACTOR) * * * PURPOSE: DETERMINE MOMENTUMSPACE-DENSITY FROM MOMENTA OF * * PSEUDOPARTICLES (ONLY IN PX-PZ PLANE) * * VARIABLES: * * P0 - BOOST IN THE C.M.-SYSTEM "GEV/C" (REAL,INPUT) * * NUM - NUMBER OF PSEUDOPARTICLES/NUCLEON (INTEGER,INPUT) * * OUT - NUMBER OF PSEUDOPART. OUT OF RANGE (INTEGER,OUTPUT) * * UNIT - STRING CONTAINING THE SPACING OF MESHPOINTS * * IN PRHO (CHAR.*9,OUTPUT) * * FACTOR - CONVERSION FACTOR (REAL,OUTPUT) * * * ************************************************************************ IMPLICIT NONE INTEGER OUT REAL P0, FACTOR CHARACTER * 9 UNIT INTEGER MAXPAR, MAXR PARAMETER (MAXPAR=100000, MAXR=1000) INTEGER NUM, MASS, MASSR REAL P, PRHO COMMON /BB/ P(3,MAXPAR) COMMON /PP/ PRHO(-20:20,-24:24) COMMON /MASS/ MASS COMMON /RUN/ NUM COMMON /RR/ MASSR(0:MAXR) INTEGER IX, IZ, ISO, NRUN, I * IF (0.0 .LE. P0 .AND. P0 .LE. 0.85) THEN UNIT = ' 50 MEV/C' FACTOR = 20.0 ELSE IF (P0 .LE. 2.1) THEN UNIT = '100 MEV/C' FACTOR = 10.0 ELSE IF (P0 .LE. 4.6) THEN UNIT = '200 MEV/C' FACTOR = 5.0 ELSE IF (P0 .LE. 9.4) THEN UNIT = '400 MEV/C' FACTOR = 2.5 ELSE IF (P0 .LE. 19.4) THEN UNIT = '800 MEV/C' FACTOR = 1.25 ELSE UNIT = 'OVERFLOW ' FACTOR = 0.0 END IF * DO IZ = -24,24 DO IX = -20,20 PRHO(IX,IZ) = 0.0 end do end do * OUT = 0 * iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1 + iso, mass + iso IX = NINT( FACTOR * P(1,I) ) IZ = NINT( FACTOR * P(3,I) ) IF( IX .LE. -20 .OR. IX .GE. 20 .OR. & IZ .LE. -24 .OR. IZ .GE. 24 ) THEN OUT = OUT + 1 ELSE PRHO(IX,IZ) = PRHO(IX,IZ) + 1.0 END IF end do end do * DO IZ = -24,24 DO IX = -20,20 PRHO(IX,IZ) = PRHO(IX,IZ) / FLOAT(NUM) end do end do * RETURN END ************************************************************************ * * SUBROUTINE PHDENS(FACTOR, OUT, DBOX) * * * PURPOSE: DETERMINE PHASE-SPACE-DENSITY IN Z-PZ-PLANE * * VARIABLES: * * NUM - NUMBER OF PSEUDOPARTICLES/NUCLEON (INTEGER,INPUT) * * FACTOR - CONVERSION FACTOR (REAL,OUTPUT) * * OUT - NUMBER OF PSEUDOPART. OUT OF RANGE (INTEGER,OUTPUT) * * * ************************************************************************ IMPLICIT NONE INTEGER OUT REAL DBOX, FACTOR INTEGER MAXPAR, MAXR, MAXZ PARAMETER (MAXPAR=100000, MAXR=1000, MAXZ=81) INTEGER MASS, NUM, MASSR REAL R, P, PHRHO COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /QQ/ PHRHO(-MAXZ:MAXZ,-24:24) COMMON /MASS/ MASS COMMON /RUN/ NUM COMMON /RR/ MASSR(0:MAXR) INTEGER ISO, NRUN, I1, I2, I REAL FRACTN * DO I2 = -24,24 DO I1 = -MAXZ,MAXZ PHRHO(I1,I2) = 0.0 end do end do * write (*,*) 'num,dbox',num,dbox FRACTN = 1.0 / FLOAT(NUM) / dbox**3 * iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1 + iso, mass + iso I1 = NINT( R(3,I)/dbox) I2 = NINT( FACTOR * P(3,I) ) IF( I1 .LE. -MAXZ .OR. I1 .GE. MAXZ .OR. & I2 .LE. -24 .OR. I2 .GE. 24 ) THEN OUT = OUT + 1 ELSE PHRHO(I1,I2) = PHRHO(I1,I2) + FRACTN END IF end do end do * RETURN END ************************************************************************ subroutine pincol(i1, lb1, i2, lb2, nstar, donfl, lres, iseed, & ix1, iy1, iz1, ipx1, ipy1, ipz1, & ix2, iy2, iz2, ipx2, ipy2, ipz2, & pcx, pcy, pcz, dt) * * * purpose: perform pion-nucleon collision * * nstar: =1 N* Resonance included * * donfl: =1: no collision; do next particle * ************************************************************************ IMPLICIT NONE INTEGER I1, LB1, I2, LB2, NSTAR, DONFL, LRES, ISEED INTEGER IX1, IY1, IZ1, IPX1, IPY1, IPZ1 INTEGER IX2, IY2, IZ2, IPX2, IPY2, IPZ2 REAL PCX, PCY, PCZ, DT INTEGER MAXPAR, MX, MY, MZ, MPX, MPY, MPZ, MPZP REAL PI, amp, ap1 PARAMETER (MAXPAR=100000) PARAMETER (PI=3.14159, amp=0.938272, ap1=0.134974) PARAMETER (mx=4, my=4, MZ=15, mpx=4, mpy=4, MPZ=5, MPZP=10) INTEGER PROTON, NEUTRON, PIONM, PION0, PIONP INTEGER DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP, ID, LB REAL MASSPA(11) CHARACTER*1 FF COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) COMMON /PARTICLE/ MASSPA, PROTON, NEUTRON, PIONM, PION0, PIONP, 1 DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP c local variables INTEGER I, J, IC, NSRESO REAL XMAX, EC, XNPIOND, DELTAR, ds, x1 REAL XMAXN, XNPIONN, XX, LBP, CLEB2 c functions REAL XNPION LOGICAL RESONANCE, PION c donfl=0 if subroutine completed succsessfull donfl = 1 *IF PION-NUCLEON COLLISION THEN *1. PION(+)+PROTON-->DELTA++,PION(-)+NEUTRON-->DELTA(-) *2. PION(-)+PROTON-->DELTA0 ,PION(+)+NEUTRON-->DELTA+ *3. PION0 +PROTON-->DELTA+ ,PION0 +NEUTRON-->DELTA0 IF ( (lb1.eq.pionp .and. lb2.eq.proton ) .or. & (lb1.eq.pionm .and. lb2.eq.neutron) .or. & (lb2.eq.pionp .and. lb1.eq.proton ) .or. & (lb2.eq.pionm .and. lb1.eq.neutron) ) then XMAX = 198.3 cleb2 = 1 else IF ( (lb1.eq.pionm .and. lb2.eq.proton ) .or. & (lb1.eq.pionp .and. lb2.eq.neutron) .or. & (lb2.eq.pionm .and. lb1.eq.proton ) .or. & (lb2.eq.pionp .and. lb1.eq.neutron) ) then XMAX = 66.1 cleb2 = 0.33333333333 else IF ( (lb1.eq.pion0 .and. lb2.eq.proton ) .or. & (lb1.eq.pion0 .and. lb2.eq.neutron) .or. & (lb2.eq.pion0 .and. lb1.eq.proton ) .or. & (lb2.eq.pion0 .and. lb1.eq.neutron) ) then XMAX = 132.2 cleb2 = 0.6666666666 else stop 'PINCOL: this is not a pion-nuclon collision' end if * FOR PION+NUCLEON SYSTEM, THE MINIMUM S is determined using lightest * nucleon and pion masses (EC = PION MASS+NUCLEON MASS+20MEV)**2 EC = (ap1 + amp + .02)**2 XNPIOND = XNPION(I1,I2,1,cleb2) IF (XNPIOND.LE.0.000001) RETURN ic = -1 NSRESO = 0 IF ((NSTAR.NE.1) .or. ( cleb2.eq.1 )) THEN * ONLY DELTA RESONANCE IS POSSIBLE DELTAR = SQRT(XMAX /(10.*PI)) x1 = (XNPIOND*10.)/XMAX ELSE * NOW BOTH DELTA AND N* RESONANCE ARE POSSIBLE * DETERMINE THE RESONANT STATE BY USING OF THE MONTE CARLO METHOD * FOR N* RESONANCE IS 50MB XMAXN = 50. XNPIONN = XNPION(I1,I2,0,XMAXN) XX = XNPIONN/(XNPIONN+XNPIOND) IF (RAN(ISEED).LT.XX) THEN C N* RESONANCE IS SELECTED DELTAR = SQRT(XMAXN/(10.*PI)) x1 = (XNPIONN*10.)/XMAXN NSRESO = 1 ELSE c DELTA RESONANCE IS SELECTED DELTAR = SQRT(XMAX /(10.*PI)) x1 = (XNPIOND*10.)/XMAX END IF END IF ds = 1.5*deltar CALL DISTANCE(I1,I2,Ds,DELTAR,DT,EC,IC,PCX,PCY,PCZ) if (ic .eq. -1) return IF ( RAN(ISEED).GE.x1) return LRES = LRES + 1 if (pion(lb1)) then i = i2 j = i1 else I = I1 j = i2 end if lbp = lb(j) IF ( NSRESO.EQ.0 ) THEN *DELTA IS PRODUCED;IN THE FOLLOWING WE DETERMINE THE *CHARGE STATE OF THE PRODUCED DELTA * (1) p+pion(+)-->DELTA(++) IF ( (lb1.eq.pionp .and. lb2.eq.proton ) .or. & (lb2.eq.pionp .and. lb1.eq.proton ) ) then LB(I)=deltapp * (2) p+pion(0)-->delta(+) else IF ( (lb1.eq.pion0 .and. lb2.eq.proton ) .or. & (lb2.eq.pion0 .and. lb1.eq.proton ) ) then LB(I)=deltap * (3) n+pion(+)-->delta(+) else IF ( (lb1.eq.pionp .and. lb2.eq.neutron ) .or. & (lb2.eq.pionp .and. lb1.eq.neutron ) ) then LB(I)=deltap * (4) n+pion(0)-->delta(0) else IF ( (lb1.eq.pion0 .and. lb2.eq.neutron ) .or. & (lb2.eq.pion0 .and. lb1.eq.neutron ) ) then LB(I)=delta0 * (5) p+pion(-)-->delta(0) else IF ( (lb1.eq.pionm .and. lb2.eq.proton ) .or. & (lb2.eq.pionm .and. lb1.eq.proton ) ) then LB(I)=delta0 * (6) n+pion(-)-->delta(-) else IF ( (lb1.eq.pionm .and. lb2.eq.neutron ) .or. & (lb2.eq.pionm .and. lb1.eq.neutron ) ) then LB(I)=deltam END IF else *N* IS PRODUCED,WE DETERMINE THE CHARGE STATE OF THE PRODUCED N* * (1.1) n+pion(+)-->N*(+) IF ( (lb1.eq.pionp .and. lb2.eq.neutron ) .or. & (lb2.eq.pionp .and. lb1.eq.neutron ) ) then LB(I)=nsp * (1.2) p+pion(0)-->N*(+) else IF ( (lb1.eq.pion0 .and. lb2.eq.proton ) .or. & (lb2.eq.pion0 .and. lb1.eq.proton ) ) then LB(I)=nsp * (1.3) n+pion(0)-->N*(0) else IF ( (lb1.eq.pion0 .and. lb2.eq.neutron ) .or. & (lb2.eq.pion0 .and. lb1.eq.neutron ) ) then LB(I)=ns0 * (1.4) p+pion(-)-->N*(0) else IF ( (lb1.eq.pionm .and. lb2.eq.proton ) .or. & (lb2.eq.pionm .and. lb1.eq.proton ) ) then LB(I)=ns0 END IF END IF c calculate mass and momenta after pion is absorbed CALL DRESONANCE(I1,I2) * PAULI BLOCKING FOR DELTA OR N* IS NEGLECTED. CHANGE PHASE SPACE * DENSITY OF THE NUCLEON DUE TO RESONANCE . IF ( resonance(i1) ) THEN if ((abs(ix1).le.mx) .and. (abs(iy1).le.my) .and. 1 (abs(iz1).le.mz)) then if ((abs(ipx1).le.mpx) .and. (abs(ipy1).le.mpy) & .and. (ipz1.ge.-mpz) .and. (ipz1.le.mpzp)) & ff(ipz1,ipy1,ipx1,ix1,iy1,iz1) = & char(ichar(ff(ipz1,ipy1,ipx1,ix1,iy1,iz1)) - 1) end if end if IF ( resonance(i2) ) THEN if ((abs(ix2).le.mx) .and. (abs(iy2).le.my) .and. & (abs(iz2).le.mz)) then if ((abs(ipx2).le.mpx) .and. (abs(ipy2).le.mpy) & .and. (ipz2.ge.-mpz) .and. (ipz2.le.mpzp)) & ff(ipz2,ipy2,ipx2,ix2,iy2,iz2) = & char(ichar(ff(ipz2,ipy2,ipx2,ix2,iy2,iz2)) - 1) end if end if donfl=0 return end ************************************************************************ SUBROUTINE PPCORR(NT, DENCUT, DBOX) * * * This routine finds which time the individual particle last * * experienced a density greater than dencut and store the * * r and p coordinates, energy and lb of the particle * * * ************************************************************************ IMPLICIT NONE INTEGER NT, NTMAX, MAXPAR, MAXX, MAXZ, MAXR REAL DENCUT, DBOX, DT, B PARAMETER (MAXPAR=100000, MAXR=1000, MAXX=40, MAXZ=81) CHARACTER*20 NAMEOUT, NAMERES CHARACTER*50 DIRECT INTEGER NUM, MASS, ID, LB, MASSR REAL R, P, E, RHO COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /FILE/ DIRECT, NAMEOUT, NAMERES COMMON /MASS/ MASS COMMON /RR/ MASSR(0:MAXR) COMMON /RUN/ NUM INTEGER ISO, I, J, NRUN, IX, IY, IZ INTEGER TIMES(MAXPAR), LBOUT(MAXPAR), idout(maxpar) REAL DICHTE REAL ROUT(3,MAXPAR), POUT(3,MAXPAR), EOUT(MAXPAR) SAVE TIMES, LBOUT, idout, ROUT, POUT, EOUT ISO = 0 J = 0 DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO I = 1 + ISO, MASS + ISO J = J+1 IX = NINT(R(1,I)/DBOX) IY = NINT(R(2,I)/DBOX) IZ = NINT(R(3,I)/DBOX) IF (ABS(IX).LE.MAXX .AND. ABS(IY).LE.MAXX .AND. 1 ABS(IZ).LE.MAXZ) THEN DICHTE = RHO(IX, IY, IZ) IF ((DICHTE .GE. DENCUT) .or. (lb(i) .ne. lbout(j)) .or. 1 (id(i) .ne. idout(j))) THEN TIMES(J) = NT idOUT(J) = id(I) LBOUT(J) = LB(I) EOUT(J) = E(I) ROUT(1,J) = R(1,I) ROUT(2,J) = R(2,I) ROUT(3,J) = R(3,I) POUT(1,J) = P(1,I) POUT(2,J) = P(2,I) POUT(3,J) = P(3,I) END IF END IF END DO END DO RETURN ENTRY PPINIT ISO = 0 J = 0 DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO I = 1 + ISO, MASS + ISO J = J + 1 TIMES(J) = 0 idOUT(J) = id(I) LBOUT(J) = LB(I) EOUT(J) = E(I) ROUT(1,J) = R(1,I) ROUT(2,J) = R(2,I) ROUT(3,J) = R(3,I) POUT(1,J) = P(1,I) POUT(2,J) = P(2,I) POUT(3,J) = P(3,I) END DO END DO RETURN ENTRY PPRINT(NTMAX, DT, B) *----------------------------------------------------------------------- * Write out distribution of free nucleons at their last * interaction times. * x,y,z [fm], px,py,pz [MeV/c], times [fm/c] OPEN (95,STATUS='NEW',FILE=DIRECT//NAMEOUT//'.D95') ISO = 0 J = 0 iy = 0 DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO I = 1 + ISO, MASS + ISO J = J + 1 IF ((TIMES(J) .LT. NTMAX ) ) THEN iy = iy + 1 WRITE (95,'(2I7,F6.2,I5,F7.2,6F9.2,i10)') J, NRUN, B, 1 LBOUT(J), DT*TIMES(J), (ROUT(IX,J),IX=1,3), 1 (POUT(IX,J)*1000.,IX=1,3) END IF END DO END DO CLOSE (95) CLOSE (96) WRITE(*,*) WRITE(*,*) B, 100.*B/FLOAT(NUM), iy, ' = B[FM],WEIGHT,#PART.' RETURN END ************************************************************************ * * SUBROUTINE PRIBIG(ISUM, RFLAG, PFLAG, RPFLAG, NT, PZPR, PZTA,dbox) * * * PURPOSE: PROVIDE LARGE CONTROL-OUTPUT USING SUBROUTINE MATRI2 * * AND PLOT1D * * VARIABLES (ALL INTEGER, ALL INPUT) * * ISUM - OUTPUT UNIT * * RFLAG - (=1 -> PRINTOUT OF X-Z-DENSITY PLOT FOR Y=0) * * (=2 -> PRINTOUT OF X-Z-DENSITY INTEGRATED OVER Y) * * PFLAG - (=1 -> PRINTOUT OF PX-PZ-MOMENTUM-DENSITY) * * RPFLAG - (=1 -> PRINTOUT OF Z-PZ-PHASE-SPACE-PLOT) * * NT - TIME STEP * * NUM - NUMBER OF PARALLEL RUNS * * PZPR - PROJECTILE MOMENTUM (REAL) * * PZTA - TARGET MOMENTUM (REAL) * * * ************************************************************************ IMPLICIT NONE INTEGER ISUM, RFLAG, PFLAG, RPFLAG, NT REAL PZPR, PZTA, DBOX INTEGER MAXX, MAXZ, MAXXP, MAXZP PARAMETER (MAXX=40, MAXZ=81) PARAMETER (MAXXP=20, MAXZP=40) INTEGER IX, IY, IZ, NOUT, NX, NZ, OUTPAR REAL FACTOR, RELP, PLOT(60,90) CHARACTER* 9 UNITS CHARACTER*20 CTIME CHARACTER*60 TITLE INTEGER NUM REAL RHO, PRHO, PHRHO COMMON /DD/ RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ) COMMON /PP/ PRHO(-20:20,-24:24) COMMON /QQ/ PHRHO(-MAXZ:MAXZ,-24:24) COMMON /RUN/ NUM *----------------------------------------------------------------------- FACTOR = 0.0 NOUT = 0 OUTPAR = 0 WRITE(CTIME,'(''(TIME STEP = '',I6,'')'')') NT *----------------------------------------------------------------------- IF (RFLAG .EQ. 1) THEN DO IZ = -MAXZP,MAXZP NZ = IZ + MAXZP + 1 DO IX = -MAXXP,MAXXP NX = IX + MAXXP + 1 PLOT(NX,NZ) = RHO(IX,0,IZ) end do end do TITLE = 'X-Z-DENSITY-DISTRIBUTION AT Y=0 '//CTIME CALL MATRI2(ISUM,PLOT,2*MAXXP+1,2*MAXZP+1,TITLE) ELSE IF (RFLAG .EQ. 2) THEN DO IZ = 1,60 DO IX = 1,60 PLOT(IX,IZ) = 0.0 end do end do DO IZ = -MAXZP,MAXZP NZ = IZ + MAXZP + 1 DO IY = -MAXX,MAXX DO IX = -MAXXP,MAXXP NX = IX + MAXXP + 1 PLOT(NX,NZ) = PLOT(NX,NZ) + RHO(IX,IY,IZ) end do end do end do TITLE = 'X-Z-DENSITY-DISTRIBUTION, ALL Y '//CTIME CALL MATRI2(ISUM,PLOT,2*MAXXP+1,2*MAXZP+1,TITLE) END IF *----------------------------------------------------------------------- IF (PFLAG .EQ. 1) THEN RELP = ABS(PZTA) IF (ABS(PZPR) .GT. ABS(PZTA)) RELP = ABS(PZPR) CALL PDENS(RELP, OUTPAR, UNITS, FACTOR) IF (OUTPAR .NE. 0) THEN WRITE(ISUM,'('' WARNING: '',I7,'' PSEUDPARTICLES NOT IN THE '', & ''FOLLOWING MOMENTUM-DENSITY PLOT'')') OUTPAR END IF DO IZ = -24,24 NZ = IZ + 24 + 1 DO IX = -20,20 NX = IX + 20 + 1 PLOT(NX,NZ) = PRHO(IX,IZ) end do end do TITLE = 'PX-PZ-DENSITY FOR ALL PY "'//UNITS//'" '//CTIME CALL MATRI2(ISUM,PLOT,41,49,TITLE) END IF *----------------------------------------------------------------------- IF (RPFLAG .EQ. 1) THEN CALL PHDENS(FACTOR, NOUT, DBOX) IF (NOUT .NE. 0) THEN WRITE(ISUM,'('' WARNING: '',I7,'' PSEUDPARTICLES NOT IN THE '', & ''FOLLOWING PHASE-SPACE-DENSITY PLOT'')') NOUT END IF * DO IZ = -24,24 NZ = IZ + 24 + 1 DO IX = -MAXZP,MAXZP NX = IX + MAXZP + 1 PLOT(NX,NZ) = PHRHO(IX,IZ) end do end do * TITLE = 'PHASE-SPACE-DENSITY; X: Z-COORD.; Y: PZ "'//UNITS//'"' CALL MATRI2(ISUM,PLOT,2*MAXZP+1,49,TITLE) END IF * RETURN END ************************************************************************ SUBROUTINE RAPIDD(NT) * * * calculates the rapidity distribution for given testparticle * * momenta * * * * OUTPUT: ny(y) - Number of tespaticles with rapidity y/10 * * per(y) - Momentum distribution * * perr(y) - reduced Momentum distribution * * pex - <|Px|> * * pxw - * * pez - <|Pz|> * ************************************************************************ IMPLICIT NONE INTEGER NT INTEGER MAXPAR, NB, MAXR REAL AVMASS PARAMETER (MAXPAR=100000, nb=200, maxr=1000, AVMASS=0.938919) INTEGER MASS, NUM, MASSR REAL P COMMON /RR/ MASSR(0:MAXR) COMMON /MASS/ MASS COMMON /RUN/ NUM COMMON /BB/ P(3, MAXPAR) INTEGER I, J, NY(NB), MAX, MIN, ISO, NRUN REAL Y, E, PZ, PER(NB), PEX, PEZ, PXW, PT, PX, P2 REAL PEXER, PXWER, PEZER, PTER, PP, PPER, ERR(NB) do i = 1,nb ny(i) = 0 per(i) = 0. err(i) = 0. end do pp = 0. pex = 0. pez = 0. pxw = 0. pt = 0. pper = 0. pexer = 0. pezer = 0. pxwer = 0. pter = 0. ISO = 0 DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO i = 1 + ISO, MASS + ISO px = p(1,i) pz = p(3,i) p2 = px**2 + p(2,i)**2 pp = pp + sqrt(p2) pper = pper + p2 pex = pex + abs(px) pexer = pexer + px**2 pt = pt + px pter = pter + pt**2 pez = pez + abs(pz) pezer = pezer + pz**2 pxw = pxw + px*sign(1.,pz) pxwer = pxwer + (px*pz)**2 p2 = p2 + pz**2 e = sqrt(AVMASS**2 + p2) * Rapidity for testparticle i y = 0.5*log((e+pz)/(e-pz)) if ( abs(y).lt.2.5 ) then j = int(40.*(y+2.5)+1.) ny(j) = ny(j)+1 per(j) = per(j)+px err(j) = err(j)+px*px end if end do end do * Write the Rapidity and the meanvalues to file write (20,*) write (20,*) 'time = ',nt write (20,*) write (20,*) ' Rapidity st. error Num. of Part' do i = 1,nb if (ny(i).ne.0) then per(i) = per(i)/ny(i) err(i) = sqrt((err(i)/ny(i)-per(i)**2+1e-6)/ny(i)) max = i end if end do do i = nb,1,-1 if (ny(i).ne.0) then min = i end if end do do i = min,max write (20,100) .025*i-2.6,.025*i-2.5,per(i),err(i),' ',ny(i) 100 format (2f6.2,2f10.6,a3,i5) end do pp = pp/num/mass pex = pex/num/mass pez = pez/num/mass pxw = pxw/num/mass pt = pt/num/mass pper = sqrt((pper/num/mass-pp**2)/num/mass) pexer = sqrt((pexer/num/mass-pex**2)/num/mass) pezer = sqrt((pezer/num/mass-pez**2)/num/mass) pxwer = sqrt((pxwer/num/mass-pxw**2)/num/mass) pter = sqrt((pter /num/mass-pt**2 )/num/mass) write (20,*) write (20,*) ' <|Px|> <|Pz|> ' write (20,'(5f11.7)') pex,pxw,pez,pt,pp write (20,'(5f11.7)') pexer,pxwer,pezer,pter, pper write (20,*) return end ************************************************************************ SUBROUTINE RELBAR(I1, I2, DM, N12) * * * purpose: change energy & particle label according to channel N12 * * for nucleon+nucleon->nucleon+delta * * i1, i2: particle number * * n12: channel * ************************************************************************ IMPLICIT NONE INTEGER I1, I2, N12, MAXPAR PARAMETER (MAXPAR=100000) INTEGER PROTON, NEUTRON, PIONM, PION0, PIONP, ID, LB INTEGER DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP REAL MASSPA(11), E, DM COMMON /PARTICLE/ MASSPA, PROTON, NEUTRON, PIONM, PION0, PIONP, 1 DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) LOGICAL RESONANCE *1. p+n-->delta(+)+n IF (N12.EQ.1) THEN IF (LB(I1).EQ.proton) THEN LB(I1)=deltap ELSE LB(I2)=deltap END IF *2 p+n-->delta(0)+p ELSE IF (N12.EQ.2) THEN IF (LB(I1).EQ.neutron) THEN LB(I1)=delta0 ELSE LB(I2)=delta0 ENDIF *3 p+p-->delta(++)+n ELSE IF (N12.EQ.3) THEN LB(I1)=deltapp LB(I2)=neutron *4 p+p-->delta(+)+p ELSE IF (N12.EQ.4) THEN LB(I1)=deltap *5 n+n--> delta(0)+n ELSE IF (N12.EQ.5) THEN LB(I1)=delta0 *6 n+n--> delta(-)+p ELSE IF (N12.EQ.6) THEN LB(I1)=deltam LB(I2)=proton *7 n+p--> N*(0)+p ELSE IF (N12.EQ.7) THEN IF (LB(I1).EQ.proton) THEN LB(I2)=ns0 ELSE LB(I1)=ns0 END IF *8 n+p--> N*(+)+n ELSE IF (N12.EQ.8) THEN IF (LB(I1).EQ.proton) THEN LB(I1)=nsp ELSE LB(I2)=nsp END IF ELSE stop 'RELBAR: channel not available' END IF if ( resonance(lb(i1)) ) then e(i1)=dm e(i2)=masspa(lb(i2)) else if ( resonance(lb(i2)) ) then e(i1)=masspa(lb(i1)) e(i2)=dm else stop 'RELBAR: channel not available (2)' end if return end ************************************************************************ * * SUBROUTINE RELCOL(ISEED, IAVOID, DT, SIGRED, DBOX, LCOLL, LBLOC, 1 LCNNE, LCNND, LCNDN, LDIRT, LDECAY, LDECBL, LRES, NT, SIGF) * * * PURPOSE: CALCULATING THE KINEMATICS IN A COLLISION BETWEEN * * TWO PARTICLES - RELATIVISTIC FORMULA USED * * * * REFERENCES: HAGEDORN, RELATIVISTIC KINEMATICS (1963) * * * * VARIABLES: * * NUM - NUMBER OF TESTPARTICLES PER NUCLEON(INTEGER,INPUT) * * ISEED - SEED FOR RANDOM NUMBER GENERATOR (INTEGER,INPUT) * * IAVOID - (= 1 => AVOID FIRST CLLISIONS WITHIN THE SAME * * NUCLEUS, ELSE ALL COLLISIONS) (INTEGER,INPUT) * * DELTAR - MAXIMUM SPATIAL DISTANCE FOR WHICH A COLLISION * * STILL CAN OCCUR (REAL,INPUT) * * DT - TIME STEP SIZE (REAL,INPUT) * * LCOLL - NUMBER OF COLLISIONS (INTEGER,OUTPUT) * * LBLOC - NUMBER OF PULI-BLOCKED COLLISIONS (INTEGER,OUTPUT) * * LCNNE - NUMBER OF ELASTIC COLLISION (INTEGER,OUTPUT) * * LCNND - NUMBER OF N+N->N+DELTA REACTION (INTEGER,OUTPUT) * * LCNDN - NUMBER OF N+DELTA->N+N REACTION (INTEGER,OUTPUT) * * * * LB(I) IS USED TO LABEL PARTICLE'S CHARGE STATE * LB(I) = 1 PROTON * 2 NUETRON * 3 PION- * 4 PION0 * 5 PION+ * 6 DELTA-1 * 7 DELTA0 * 8 DELTA+1 * 9 DELTA+2 * 10 N*0 * 11 N*(+1) * NSTAR=1 INCLUDING N* RESONANCE * ELSE DELTA RESONANCE ONLY * NDIRECT=1 INCLUDING DIRECT PROCESS,ELSE NOT * DIR - PERCENTAGE OF DIRECT PION PRODUCTION PROCESS * sigf=0 -Cugnon elastic cross section =1 particle data group ************************************************************************ IMPLICIT NONE INTEGER ISEED, IAVOID, LCOLL, LBLOC, LCNNE, LCNND, LCNDN INTEGER LDIRT, LDECAY, LDECBL, LRES, NT, SIGF REAL DT, SIGRED, DBOX INTEGER MAXPAR, MAXP, MAXR, MX, MY, MZ, MPX, MPY, MPZ, MPZP REAL HC PARAMETER (MAXPAR=100000, MAXP=100, MAXR=1000, HC=0.197327) PARAMETER (MX=4, MY=4, mz=15, MPX=4, MPY=4, MPZ=5, MPZP=10) INTEGER NUM, MASS, ZTA, ZPR, NSTAR, NDIRECT, ID, LB, MASSR REAL DIR, DX, DY, DZ, DPX, DPY, DPZ, ZET, R, P, E CHARACTER*1 FF COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /HH/ FF(-MPZ:MPZP,-MPY:MPY,-MPX:MPX,-MX:MX,-MY:MY,-MZ:MZ) COMMON /GG/ DX, DY, DZ, DPX, DPY, DPZ COMMON /INPUT/ NSTAR, NDIRECT, DIR COMMON /RR/ MASSR(0:MAXR) COMMON /MASS/ MASS COMMON /RUN/ NUM COMMON /ZZ/ ZTA, ZPR, ZET(11) c common block connecting RELCOL, CROSS and other parts of the collision term INTEGER NNN, IPION, LPION REAL BETAX, BETAY, BETAZ, GAMMA, RPION, PPION, EPION COMMON /BG/ BETAX, BETAY, BETAZ, GAMMA COMMON /NN/ NNN COMMON /PA/ RPION(3,MAXP,MAXR) COMMON /PB/ PPION(3,MAXP,MAXR) COMMON /PC/ EPION(MAXP,MAXR) COMMON /PD/ IPION(MAXP,MAXR), LPION(MAXP,MAXR) c Common block for pion correlation function output ! added 2-12 INTEGER IPICOR, NNNTOT, MAXPI, NPICNT PARAMETER (MAXPI = MAXP*MAXR) INTEGER LPI(MAXPI), TPI(MAXPI), PICNT(0:MAXPI), 1 PICNT2(MAXP,MAXR) REAL RPI(3,MAXPI), PPI(3,MAXPI) COMMON /PIHBT/ IPICOR, nnntot, TPI, PICNT, LPI, RPI, PPI c local variables INTEGER LT(MAXPAR), IT(MAXPAR), ISWAP(MAXPAR), MASSRN(0:MAXR) REAL RT(3,MAXPAR), PT(3,MAXPAR), ET(MAXPAR) INTEGER IRUN, DONFL, I, I0, I1, I2, J1, J2, ISO INTEGER IX1, IY1, IZ1, IPX1, IPY1, IPZ1 INTEGER IX2, IY2, IZ2, IPX2, IPY2, IPZ2 INTEGER ID1, ID2, LB1, LB2, LB1T, LB2T, LBPI INTEGER IBLOCK, NTAG, NNNPI REAL EM1, EM2, AM1, AM2, OCCUP REAL PX1, PY1, PZ1, PX2, PY2, PZ2 REAL X1, Y1, Z1, X2, Y2, Z2, PCX, PCY, PCZ REAL SRT, T0, PDECAY, XDECAY c functions REAL WIDTH LOGICAL NUCLEON, PION, DELTA, RESONANCE *----------------------------------------------------------------------- * INITIALIZATION OF COUNTING VARIABLES LCOLL = 0 LBLOC = 0 LCNNE = 0 LCNND = 0 LCNDN = 0 LDIRT = 0 LDECAY = 0 LDECBL = 0 LRES = 0 ISO = 0 NNNPI = 0 MASSRN(0) = 0 *----------------------------------------------------------------------- * LOOP OVER ALL PARALLEL RUNS DO 1000 IRUN = 1, NUM c mix particle pointer iso = iso + massr(irun-1) do i = 1, massr(irun) iswap(i + iso) = i + iso end do do I1 = 1, 10 do i = 1, massr(irun) j1 = iswap(i+iso) j2 = nint(ran(iseed)*(massr(irun) - 1) + 1.) + iso iswap(i+iso) = iswap(j2) iswap(j2) = j1 end do end do NNN = 0 * LOOP OVER ALL PSEUDOPARTICLES 1 IN THE SAME RUN DO 800 J1 = 2, MASSR(IRUN) I1 = iswap(J1 + iso) IF (E(I1) .EQ. 0.0) GOTO 800 X1 = R(1,I1) Y1 = R(2,I1) Z1 = R(3,I1) PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) *keep all coordinates for possible phase space change ix1 = nint(x1/dx) iy1 = nint(y1/dy) iz1 = nint(z1/dz) ipx1 = nint(px1/dpx) ipy1 = nint(py1/dpy) ipz1 = nint(pz1/dpz) EM1 = E(I1) AM1 = EM1 ID1 = ID(I1) LB1 = LB(I1) IF ( resonance(lb1) ) then c Decay of Delta and N* resonance to nucleons and pions T0 = hc/0.2 IF ( delta (lb1) ) T0 = hc/WIDTH(EM1)/em1*sqrt(em1**2 1 + px1**2 + py1**2 + pz1**2) PDECAY = 1.-EXP(-DT/T0) XDECAY = RAN(ISEED) c decay actually happens IF (XDECAY.LT.PDECAY) THEN NNN = NNN + 1 CALL DECAY(IRUN, I1, NNN, ISEED) CALL PAULAT(I1, OCCUP) IF (RAN(ISEED).LT.OCCUP) THEN LDECBL = LDECBL+1 P(1,I1) = PX1 P(2,I1) = PY1 P(3,I1) = PZ1 E(I1) = EM1 LB(I1) = LB1 NNN = NNN-1 ELSE c increase decay counters LDECAY = LDECAY+1 * SAVE ALL THE COORDINATES FOR POSSIBLE CHANGE IN THE FOLLOWING COLLISION PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) EM1 = E(I1) AM1 = EM1 LB1 = LB(I1) IPX1 = NINT(PX1/DPX) IPY1 = NINT(PY1/DPY) IPZ1 = NINT(PZ1/DPZ) c take care of wigner function if ((abs(ix1).le.mx) .and. (abs(iy1).le.my) .and. 1 (abs(iz1).le.mz) .AND. (abs(ipx1).le.mpx) .and. 1 (abs(ipy1).le.mpy) .and. (ipz1.ge.-mpz).and. 1 (ipz1.le.mpzp)) ff(ipz1,ipy1,ipx1,ix1,iy1,iz1) = 1 char(ichar(ff(ipz1,ipy1,ipx1,ix1,iy1,iz1)) + 1) * New addition: Set pion coordinates RPION(1,NNN,IRUN) = X1 RPION(2,NNN,IRUN) = Y1 RPION(3,NNN,IRUN) = Z1 * take care of correlation function output if necessary nnntot = nnntot + 1 if (ipicor .eq. 1) then ! added 2-12 RPI(1,nnntot) = rpion(1,nnn,irun) RPI(2,nnntot) = rpion(2,nnn,irun) RPI(3,nnntot) = rpion(3,nnn,irun) PPI(1,nnntot) = ppion(1,nnn,irun) PPI(2,nnntot) = ppion(2,nnn,irun) PPI(3,nnntot) = ppion(3,nnn,irun) lpi(nnntot) = lpion(nnn,irun) TPI(nnntot) = nt picnt2(nnn,irun) = nnntot end if ENDIF ENDIF ENDIF * end of decay part * LOOP OVER ALL PSEUDOPARTICLES 2 IN THE SAME RUN DO 600 J2 = 1,J1-1 I2 = iswap(J2 + iso) c maximum cross section determines critical distance c ( pion+nucleon sigma=200mb) Z2 = R(3,I2) if ( abs(z1-z2).gt.3 ) GOTO 400 X2=R(1,I2) if ( abs(x1-x2).gt.3 ) GOTO 400 Y2=R(2,I2) if ( abs(y1-y2).gt.3 ) GOTO 400 if ( (x1-x2)**2+(y1-y2)**2+(z1-z2)**2.gt.9 ) GOTO 400 IF (E(I2).EQ.0.) GOTO 600 ID2 = ID(I2) IF (ID1*ID2.eq.1.and.IAVOID.eq.1) GOTO 400 LB2 = LB(I2) PX2 = P(1,I2) PY2 = P(2,I2) PZ2 = P(3,I2) EM2 = E(I2) AM2 = EM2 * KEEP ALL COORDINATES FOR POSSIBLE PHASE SPACE CHANGE ix2 = nint(x2/dx) iy2 = nint(y2/dy) iz2 = nint(z2/dz) ipx2 = nint(px2/dpx) ipy2 = nint(py2/dpy) ipz2 = nint(pz2/dpz) IF( pion(lb1) .AND. pion(lb2)) GOTO 400 IF( pion(lb1) .AND. resonance(lb2)) GOTO 400 IF( resonance(lb1).AND. pion(lb2)) GOTO 400 CALL CMS(I1,I2,PCX,PCY,PCZ,SRT) * SEPARATE NUCLEON+NUCLEON( BARYON RESONANCE+ BARYON RESONANCE ELASTIC * COLLISION), BARYON RESONANCE+NUCLEON AND BARYON-PION COLLISIONS if ( pion(lb1).or.pion(lb2) ) then * PION+NUCLEON COLLISION call pincol(i1, lb1, i2, lb2, nstar, donfl, lres, iseed, & ix1, iy1, iz1, ipx1, ipy1, ipz1, & ix2, iy2, iz2, ipx2, ipy2, ipz2, & pcx, pcy, pcz, dt) IF (E(I1).EQ.0.) GOTO 800 IF (donfl.eq.1 ) GOTO 400 PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) EM1 = E(I1) AM1 = EM1 LB1 = LB(I1) IPX1 = NINT(PX1/DPX) IPY1 = NINT(PY1/DPY) IPZ1 = NINT(PZ1/DPZ) IF ( ID1.NE.0 ) id(i1) = ID1+ID1/ABS(ID1) IF (E(I2).EQ.0.) GOTO 400 stop 'RELCOL: pion not absorbed' else if ( (nucleon(lb1) .or. resonance(lb1)) .and. 1 (nucleon(lb2) .or. resonance(lb2)) ) then c baryon-baryon collision call babaco(lb1, i1, lb2, i2, pcx, pcy, pcz, sigf, donfl, 1 iseed, irun, ldirt, lbloc, lcoll, ntag, 1 iblock, srt, dt, em1, em2, sigred, dbox) if ( donfl.eq.1 ) GOTO 400 * IF COLLISION BLOCKED, RESTORE THE MOMENTUM, MASSES & LABELS OF I1 & I2 IF (NTAG .EQ. -1) THEN LBLOC = LBLOC + 1 P(1,I1) = PX1 P(2,I1) = PY1 P(3,I1) = PZ1 P(1,I2) = PX2 P(2,I2) = PY2 P(3,I2) = PZ2 E(I1) = EM1 E(I2) = EM2 LB(I1) = LB1 LB(I2) = LB2 ELSE IF (IBLOCK.EQ.1) LCNNE = LCNNE + 1 IF (IBLOCK.EQ.2) LCNND = LCNND + 1 c increase NN->DN counters IF (IBLOCK.EQ.3) LCNDN = LCNDN + 1 c increase ND->NN counters PX1 = P(1,I1) PY1 = P(2,I1) PZ1 = P(3,I1) EM1 = E(I1) EM2 = E(I2) LB1 = LB(I1) LB2 = LB(I2) if ( id1.ne.0 ) ID(I1) = id1+id1/abs(id1) if ( id2.ne.0 ) ID(I2) = id2+id2/abs(id2) * change phase space density FOR NUCLEONS INVOLVED : * NOTE THAT f is the phase space distribution function for nucleons only call changph(i1,ix1,iy1,iz1,ipx1,ipy1,ipz1,am1,em1) call changph(i2,ix2,iy2,iz2,ipx2,ipy2,ipz2,am2,em2) AM1=EM1 AM2=EM2 END IF else write (*,*) lb1, lb2 stop 'relcol: collision channel not available' end if 400 CONTINUE 600 CONTINUE 800 CONTINUE * RELABEL PIONS LEFT IN THIS RUN EXCLUDING THOSE PIONS CREATED DURING * THIS TIME STEP AND COUNT THE TOTAL NO. OF PARTICLES IN THIS RUN I0 = ISO + MASS DO I = I0+1, MASSR(IRUN)+ISO IF ( E(I).GT.0.) then NNN = NNN+1 RPION(1,NNN,IRUN) = R(1,I) RPION(2,NNN,IRUN) = R(2,I) RPION(3,NNN,IRUN) = R(3,I) PPION(1,NNN,IRUN) = P(1,I) PPION(2,NNN,IRUN) = P(2,I) PPION(3,NNN,IRUN) = P(3,I) EPION(NNN,IRUN) = E(I) LPION(NNN,IRUN) = LB(I) if (ipicor .eq. 1) then npicnt = i - mass*irun picnt2(nnn, irun) = picnt(npicnt) end if end if END DO MASSRN(IRUN) = NNN + MASS NNNPI = NNNPI + MASSRN(IRUN) 1000 CONTINUE * RELABEL ALL THE PARTICLES EXISTING AFTER THIS TIME STEP IF (NNNTOT .EQ. 0) RETURN ! NO REARRANGING IF NO PIONS IF (NNNPI .GT. MAXPAR) STOP 'RELCOL: too many pions produced' NPICNT = 0 ISO = 0 I0 = 0 DO IRUN=1, NUM ISO = ISO + MASSR(IRUN-1) I0 = I0 + MASSRN(IRUN-1) I1 = I0 DO J1 = 1 + ISO, MASS + ISO I1 = I1 + 1 RT(1,I1) = R(1,J1) RT(2,I1) = R(2,J1) RT(3,I1) = R(3,J1) PT(1,I1) = P(1,J1) PT(2,I1) = P(2,J1) PT(3,I1) = P(3,J1) ET(I1) = E(J1) LT(I1) = LB(J1) IT(I1) = ID(J1) if ( et(I1).eq.0 ) stop 'ba:et(I1)=0' end do DO J1 = 1, MASSRN(IRUN) - MASS NPICNT = NPICNT + 1 I1 = I1 + 1 RT(1,I1) = RPION(1,J1,IRUN) RT(2,I1) = RPION(2,J1,IRUN) RT(3,I1) = RPION(3,J1,IRUN) PT(1,I1) = PPION(1,J1,IRUN) PT(2,I1) = PPION(2,J1,IRUN) PT(3,I1) = PPION(3,J1,IRUN) ET(I1) = EPION(J1,IRUN) LT(I1) = LPION(J1,IRUN) IT(I1) = 0 if (ipicor .eq. 1) then picnt(npicnt) = picnt2(j1, irun) end if if ( et(I1).eq.0 ) stop 'pi:et(I1)=0' end do end do ISO = 0 DO IRUN=1,NUM MASSR(IRUN) = MASSRN(IRUN) ISO = ISO + MASSR(IRUN-1) DO I1 = 1 + ISO, MASSR(IRUN) + ISO R(1,I1) = RT(1,I1) R(2,I1) = RT(2,I1) R(3,I1) = RT(3,I1) P(1,I1) = PT(1,I1) P(2,I1) = PT(2,I1) P(3,I1) = PT(3,I1) E(I1) = ET(I1) LB(I1) = LT(I1) ID(I1) = IT(I1) end do end do RETURN END ************************************************************************ SUBROUTINE RELNUC(I1, I2, M12) * * * purpose: change energy & particle label according to channel M12 * * for delta+nucleon->nucleon+nucleon * * i1, i2: particle number * * m12: channel * ************************************************************************ IMPLICIT NONE INTEGER I1, I2, M12, MAXPAR PARAMETER (MAXPAR=100000) INTEGER ID, LB INTEGER PROTON, NEUTRON, PIONM, PION0, PIONP INTEGER DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP REAL MASSPA(11), E COMMON /PARTICLE/ MASSPA, PROTON, NEUTRON, PIONM, PION0, PIONP, 1 DELTAM, DELTA0, DELTAP, DELTAPP, NS0, NSP COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) *(1) n+delta(+)-->n+p IF (M12.EQ.1) THEN IF (LB(I1).EQ.deltap) THEN LB(I1)=proton ELSE LB(I2)=proton END IF *(2) p+delta(0)-->p+n else IF (M12.EQ.2) THEN IF (LB(I1).EQ.delta0) THEN LB(I1)=neutron ELSE LB(I2)=neutron END IF *(3) n+delta(++)-->p+p OR *(4) p+delta(+)-->p+p else IF ( (M12.EQ.3) .OR. (M12.EQ.4) ) THEN LB(I1)=proton LB(I2)=proton *(5) n+delta(0)-->n+n OR *(6) p+delta(-)-->n+n else IF ( (M12.EQ.5) .OR. (M12.EQ.6) ) THEN LB(I1)=neutron LB(I2)=neutron *(7) p+N*(0)-->n+p OR *(8) n+N*(+)-->n+p else IF ( (M12.EQ.7) .OR. (M12.EQ.8) ) THEN LB(I1)=proton LB(I2)=neutron else stop 'RELNUC: channel not available' END IF e(i1)=masspa(lb(i1)) e(i2)=masspa(lb(i2)) return end ************************************************************************ SUBROUTINE RELOAD(LPR1,LNE1,LP1,LP2,LP3,LD1,LD2,LD3,LD4,LN1,LN2) C This subroutine will reload lb, id, e, p and r from data files C ************************************************************************ IMPLICIT NONE INTEGER LPR1, LNE1, LP1, LP2, LP3, LD1, LD2, LD3, LD4, LN1 INTEGER LN2, MAXPAR, MAXR PARAMETER (MAXPAR=100000, MAXR=1000) CHARACTER*20 NAMERES, NAMEOUT CHARACTER*50 DIRECT INTEGER MASS, ID, LB, MASSR REAL R, P, E COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /RR/ MASSR(0:MAXR) COMMON /MASS/ MASS COMMON /FILE/ DIRECT, NAMEOUT, NAMERES INTEGER I, NRUN, LB1, LB2 open (60,status='old',file=direct//nameres//'.d94') rewind (60) i = 0 nrun = 1 massr(0) = 0 massr(1) = 0 lb1 = 1 lb2 = 1 60 i = i+1 c read (60,'(i3,7f8.3)',err = 61) lb(i),e(i) read (60,'(i3,i5,f8.4,6f10.5)',err = 61) lb(i),id(i),e(i) & ,p(1,i),p(2,i),p(3,i) & ,r(1,i),r(2,i),r(3,i) c write (*,'(i3,7f8.3)') lb(i),e(i) c & ,p(1,i),p(2,i),p(3,i) c & ,r(1,i),r(2,i),r(3,i) massr(nrun) = massr(nrun)+1 lb1 = lb2 lb2 = lb(i) if ((lb1.le.5 .and. lb1.ge.3 .and. (lb2.le.2 .or. lb2.ge.6)) 1 .or. ((lb2.le.2 .or. lb2.ge.6) .and. massr(nrun).gt.mass)) then massr(nrun) = massr(nrun) - 1 nrun = nrun + 1 massr(nrun) = 1 end if goto 60 61 close (60) call countp(lpr1, lne1, lp1, lp2, lp3, ld1, ld2, ld3, ld4,ln1,ln2) RETURN END ************************************************************************ SUBROUTINE RESMASS(N12, ISEED, SRT, DM) * Determine delta mass via rejection method. * Parametrization of delta resonance shape according to Kitazoe's * or J.D. Jackson's mass formula and Breit Wigner formula for N*. ************************************************************************ IMPLICIT NONE INTEGER N12, ISEED, NTRY REAL SRT, avpi, avmass, AD0, AN0 PARAMETER (avmass=0.938919, avpi=0.138037, AD0=1.232, AN0=1.44) REAL DMAX, DMIN, AM0, WID0, FM, XD, WIDTH, FRES, DM DMAX = SRT - avmass - 0.005 DMIN = avmass + avpi if ( dmax.le.dmin ) RETURN IF (N12.LT.7) THEN AM0 = AD0 IF (DMAX.LT.AM0) THEN WID0 = WIDTH(DMAX) FM = FRES(DMAX, SRT, am0, DMIN, WID0) ELSE WID0 = WIDTH(AM0) FM = FRES(AM0 , SRT, am0, AM0 , WID0) END IF ELSE AM0 = AN0 WID0 = 0.2 IF (DMAX.LT.AM0) THEN FM = FRES(DMAX, SRT, am0, DMIN, WID0) ELSE FM = FRES(AM0 , SRT, am0, AM0 , WID0) END IF END IF IF (FM.EQ.0.) FM=0.00000001 NTRY=0 10 DM = RAN(ISEED) * (DMAX-DMIN) + DMIN NTRY = NTRY + 1 if (n12 .lt. 7) WID0 = WIDTH(DM) IF ((RAN(ISEED) .GT. FRES(DM, SRT, am0, DM, WID0)/FM) .AND. 1 (NTRY.LE.10)) GOTO 10 RETURN END ************************************************************************ * * SUBROUTINE SUMMRY(IO, MASSPR, MASSTA) * * * PURPOSE: PROVIDING STATISTICAL SUMMARY * * VARIABLES: * * IO - OUTPUT-UNIT (INTEGER,INPUT) * * MASSPR - PROJECTILE MASS (INTEGER,INPUT) * * MASSTA - TARGET MASS (INTEGER,INPUT) * * NUM - NUMBER OF PARALLEL RUNS (INTEGER,INPUT) * * DIMENSIONS: * * ALL COORDINATES R IN FM * * ALL MOMENTA P IN GEV/C PER NUCLEON * * ALL ENERGIES E IN GEV PER NUCLEON * * * ************************************************************************ IMPLICIT NONE INTEGER IO, MASSPR, MASSTA INTEGER MAXPAR, MAXR PARAMETER (MAXPAR=100000, maxr=1000) INTEGER MASS, NUM, ID, MASSR, LB REAL E, R, P COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) COMMON /EE/ ID(MAXPAR), lb(maxpar) COMMON /RR/ MASSR(0:MAXR) COMMON /MASS/ MASS COMMON /RUN/ NUM INTEGER I, IDMAXP, IDMAXT, IDP, IDT, IDTOTP, IDTOTT, IRUN INTEGER NUPR, NUTA, ISO, NRUN REAL E0, EKIN, EKINP, EKINT, ETOT, ETOTP, ETOTT, P2 REAL PX, PY, PZ, TOTPXP, TOTPXT, TOTPYP, TOTPYT, PIONS REAL TOTPZP, TOTPZT, TOTXP, TOTYP, TOTZP REAL TOTXT, TOTYT, TOTZT, X, Y, Z * * INITIALIZE COUNTING VARIABLES * nuta = 0 nupr = 0 TOT PX T = 0.0 TOT PY T = 0.0 TOT PZ T = 0.0 TOT PX P = 0.0 TOT PY P = 0.0 TOT PZ P = 0.0 TOT X T = 0.0 TOT Y T = 0.0 TOT Z T = 0.0 TOT X P = 0.0 TOT Y P = 0.0 TOT Z P = 0.0 E TOT P = 0.0 E TOT T = 0.0 E KIN P = 0.0 E KIN T = 0.0 IDMAXT = 0 IDTOTT = 0 IDMAXP = 0 IDTOTP = 0 * * FIND OUT TARGET AND PROJECTILE MEAN MOMENTUM AND KINETIC ENERGY, * RADIUS AND NUMBER OF COLLISIONS * iso = 0 do nrun = 1, num iso = iso + massr(nrun-1) do i = 1+iso, mass+iso if ( id(i).gt. 0 ) then IDT = ID(I) IF (IDT .GT. IDMAXT) IDMAXT = IDT PX = P(1,I) PY = P(2,I) PZ = P(3,I) X = R(1,I) Y = R(2,I) Z = R(3,I) P2 = PX**2 + PY**2 + PZ**2 E0 = E(I) ETOT = SQRT( E0**2 + P2 ) EKIN = ETOT - E0 TOT PX T = TOT PX T + PX TOT PY T = TOT PY T + PY TOT PZ T = TOT PZ T + PZ TOT X T = TOT X T + X TOT Y T = TOT Y T + Y TOT Z T = TOT Z T + Z E TOT T = E TOT T + ETOT E KIN T = E KIN T + EKIN nuta = nuta + 1 else IDP = - ID(I) IF (IDP .GT. IDMAXP) IDMAXP = IDP PX = P(1,I) PY = P(2,I) PZ = P(3,I) X = R(1,I) Y = R(2,I) Z = R(3,I) P2 = PX**2 + PY**2 + PZ**2 E0 = E(I) ETOT = SQRT( E0**2 + P2 ) EKIN = ETOT - E0 TOT PX P = TOT PX P + PX TOT PY P = TOT PY P + PY TOT PZ P = TOT PZ P + PZ TOT X P = TOT X P + X TOT Y P = TOT Y P + Y TOT Z P = TOT Z P + Z E TOT P = E TOT P + ETOT E KIN P = E KIN P + EKIN nupr = nupr + 1 end if end do end do * IF (nuta .NE. 0) THEN TOT PX T = TOT PX T / nuta TOT PY T = TOT PY T / nuta TOT PZ T = TOT PZ T / nuta TOT X T = TOT X T / nuta TOT Y T = TOT Y T / nuta TOT Z T = TOT Z T / nuta E TOT T = E TOT T / nuta E KIN T = E KIN T / nuta END IF * IF ( nupr .NE. 0) THEN TOT PX P = TOT PX P / nupr TOT PY P = TOT PY P / nupr TOT PZ P = TOT PZ P / nupr TOT X P = TOT X P / nupr TOT Y P = TOT Y P / nupr TOT Z P = TOT Z P / nupr E TOT P = E TOT P / nupr E KIN P = E KIN P / nupr END IF * * * OUTPUT * write (io,*) 'Number of baryons in traget/projectile',nuta,nupr IF (nuta .NE. 0) & WRITE(IO,'('' TA: P=('',F8.4,'';'',F8.4,'';'',F8.4,''); R=('', & F8.4,'';'',F8.4,'';'',F8.4,''); EKIN='',F8.4,'', ETOT'' & ,''='',F8.4,/t10,'' COLS(ALL;MAX)='',I6,'';'',I3)') & TOT PX T, TOT PY T, TOT PZ T, TOT X T, TOT Y T, & TOT Z T, E KIN T, E TOT T, IDTOTT, IDMAXT-1 IF (nupr .NE. 0) & WRITE(IO,'('' PR: P=('',F8.4,'';'',F8.4,'';'',F8.4,''); R=('', & F8.4,'';'',F8.4,'';'',F8.4,''); EKIN='',F8.4,'', ETOT'' & ,''='',F8.4,/t10,'' COLS(ALL;MAX)='',I6,'';'',I3)') & TOT PX P, TOT PY P, TOT PZ P, TOT X P, TOT Y P, & TOT Z P, E KIN P, E TOT P, IDTOTP, IDMAXP-1 * RETURN END ************************************************************************ subroutine sumryn(pxw, rd) * * * calculates transverse momentum and distance of the two nuclei * ************************************************************************ IMPLICIT NONE REAL PXW, RD INTEGER MAXPAR, MAXR PARAMETER (MAXPAR=100000, MAXR=1000) INTEGER NUM, MASS, MASSR, ID, LB REAL R, P COMMON /AA/ R(3,MAXPAR) COMMON /BB/ P(3,MAXPAR) COMMON /RUN/ NUM COMMON /MASS/ MASS COMMON /EE/ ID(MAXPAR), LB(MAXPAR) COMMON /RR/ MASSR(0:MAXR) INTEGER I, ISO, NRUN, NU1, NU2 REAL X1, X2, Y1, Y2, Z1, Z2, PX, PZ pxw = 0. x1 = 0. y1 = 0. z1 = 0. x2 = 0. y2 = 0. z2 = 0. nu1 = 0 nu2 = 0 ISO = 0 DO NRUN = 1, NUM ISO = ISO + MASSR(NRUN-1) DO i = 1 + ISO, MASS + ISO px = p(1,i) pz = p(3,i) pxw = pxw + px*sign(1.,pz) if ( id(i).gt.0 ) then x1 = x1 + r(1,i) y1 = y1 + r(2,i) z1 = z1 + r(3,i) nu1 = nu1 + 1 else if ( id(i).lt.0 ) then x2 = x2 + r(1,i) y2 = y2 + r(2,i) z2 = z2 + r(3,i) nu2 = nu2 + 1 end if end do end do pxw = pxw/num/mass rd = 0. if ( nu1.ne.0.and.nu2.ne.0 ) then x1 = x1/nu1 y1 = y1/nu1 z1 = z1/nu1 x2 = x2/nu2 y2 = y2/nu2 z2 = z2/nu2 rd = sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2) end if return end ************************************************************************ * * SUBROUTINE WRSTEP(I) * * * PURPOSE: THIS SUBROUTINE WRITES OUT TIME STEP TO LOG-FILE * * IN A SPACE SAVING WAY * * INPUT: I - TIME STEP NUMBER * * * ************************************************************************ IF(I/20*20.EQ.I) THEN WRITE(*,'(''+'',I3)') I ELSE IF((I-1)/20*20.EQ.(I-1)) THEN WRITE(*,'(I5,'',''$)') I ELSE WRITE(*,'(''+'',I3,'',''$)') I END IF RETURN END ************************************************************************ REAL FUNCTION DENOM(SRT, DELT) * DATE : NOV. 15, 1991 ************************************************************************ IMPLICIT NONE REAL SRT REAL PI, avmass, avpi, AN0, AD0 PARAMETER (Avpi=0.138037, PI=3.14159, avmass=0.938919, 1 AN0=1.44, AD0=1.232) INTEGER I REAL AM0, AMIN, AMAX, DMASS, AVMASS2, SPLUSM2 REAL TQAM0, A1, P0, S, DM, WEIGHT, Q2 LOGICAL DELT IF (DELT) THEN AM0 = AD0 ELSE AM0 = AN0 END IF AVMASS2 = AVMASS*AVMASS AMAX = SRT-AVMASS AMIN = AVMASS+AVPI DMASS = 0.001*(AMAX-AMIN) DENOM = 0.0 DM = AMIN S = SRT**2 SPLUSM2 = S + AVMASS2 TQAM0 = 0.2*AM0 WEIGHT = 0.5 DO I = 0, 999 IF (AMAX .LE. DM) GOTO 12 IF (DELT) THEN Q2 = ((DM**2-AVMASS**2+AVPI**2)/(2.*DM))**2-AVPI**2 TQAM0 = 0.47*AM0*(Q2**1.5)/(AVPI**2+0.6*Q2) END IF A1 = 4.0*TQAM0*AM0/(TQAM0**2 + (DM**2-AM0**2)**2) P0 = 0.25*(SPLUSM2 - DM*DM)**2/S - AVMASS2 IF (P0.LT.0.) stop 'denom: p0 neg' DENOM = DENOM + WEIGHT*DM*A1*SQRT(P0) DM = DM + DMASS WEIGHT = 1.0 END DO 12 DENOM = DENOM*DMASS/(2.*PI) RETURN END ************************************************************************ REAL FUNCTION FRES(DM, SRT, AM0, dm2, WIDTH) * Function FRES(DMASS) gives resonance mass distribution ************************************************************************ IMPLICIT NONE REAL dm2, DM, SRT, AM0, WIDTH, AVMASS PARAMETER (AVMASS=0.938919) REAL FD, P1, M0WIDTH M0WIDTH = AM0*WIDTH FD = 4.0*AM0*M0WIDTH/((DM**2-AM0**2)**2 + M0WIDTH**2) P1 = (SRT**2+DM2**2-AVMASS**2)**2/(4.*SRT**2)-DM2**2 FRES = FD*SQRT(P1)*dm2 RETURN END ************************************************************************ real function nurho(r, a, r0) * * * Density of a nucleus given by Fourier-Bessel Coefficients a * * H.De Vries Atomic Data and Nuclear Data Tables 36, p.495(1987) * * * ************************************************************************ IMPLICIT NONE INTEGER I REAL R, A(17), R0, PI, J0, S, X PARAMETER (PI=3.14159) if ( r.gt.r0 ) then nurho = 0. else s = 0. x = pi*r/r0 do i = 1,17 if ( x.le.1e-6) then j0 = 1. else j0 = sin(x*i)/x/i end if s = s + a(i)*j0 end do nurho = s end if return end ************************************************************************ REAL FUNCTION SIGMA(I, ID, SRT) * * *PURPOSE : THIS IS THE PROGRAM TO CALCULATE THE ISOSPIN DECOMPOSED * * CROSS SECTION BY USING OF B.J.VerWEST AND R.A.ARNDT'S * * PARAMETERIZATION * *REFERENCE: PHYS. REV. C25(1982)1979 * *QUANTITIES: i: 1=sigma_10^d,2=sigma_11, 3=sigma_10, 4=sigma_01 * * id: =1 FOR DELTA RESONANCE =2 FOR N* RESONANCE * * * * use same masses as defined in paper: avmass=0.9389, avpi=0.1380 * * cutoff at srt=2.0159 * ************************************************************************ IMPLICIT NONE INTEGER I, ID REAL SRT, AVMASS, AVPI, PI, HC PARAMETER (AVMASS=0.9389, AVPI=0.1380, PI=3.14159, 1 HC=0.197327) REAL ALPHA(4), BETA(4), M0(4), GAMMA(4) DATA ALPHA /6.03 , 3.772 , 15.28 , 146.3 / DATA BETA /1.7 , 1.262 , 0.0 , 0.0 / DATA M0 /1.203 , 1.188 , 1.245 , 1.472 / DATA GAMMA /0.1343, 0.09902, 0.1374, 0.02649/ real s0, p02, pr2, p2, spin, s, ss, gamma0, q2, q02 REAL amass0, amass, zplus, zminus, ss0 c define a function which finds momentum squared for a given s, m1 C and m2 REAL P2FUN,SF,M1F,M2F P2FUN(SF,M1F,M2F) = (SF - (M1F - M2F)**2) 1 *(SF - (M1F + M2F)**2)*0.25/SF IF ( SRT .LE.1E-6 ) STOP 'SIGMA: S less than threshold' S = SRT*SRT SIGMA = 0.0 S0 = (AVMASS + M0(I))**2 P02 = 0.25*S0 - AVMASS**2 IF (P02 .LE. 0.0) RETURN P2 = 0.25*S - AVMASS**2 IF (P2 .LT. 0.0) RETURN IF ( I.EQ.1 ) THEN PR2 = P2FUN(s,2.*AVMASS,AVPI) IF (PR2 .LT. 0.0) RETURN SPIN = (SRT - AVMASS)**2 SIGMA = 1.0/((SPIN-M0(I)**2)**2+M0(I)**2*GAMMA(I)**2) ELSE IF(ID.EQ.1)THEN ! delta AMASS0 = 1.22 GAMMA0 = 0.12 ELSE ! N* AMASS0 = 1.43 GAMMA0 = 0.2 END IF ZPLUS = (SRT - AVMASS - AMASS0)*2./GAMMA0 ZMINUS = (AVMASS + AVPI - AMASS0)*2./GAMMA0 AMASS = AMASS0+(gamma0*0.25)*LOG((1.+ZPLUS**2)/(1.+ZMINUS**2)) 1 /(ATAN(ZPLUS)-ATAN(ZMINUS)) PR2 = P2FUN(S,AVMASS,AMASS) IF (PR2 .LT. 0.0) RETURN SS = AMASS**2 Q2 = P2FUN(SS,AVMASS,AVPI) IF (Q2 .LT. 0.0) RETURN SS0 = M0(i)**2 Q02 = P2FUN(SS0,AVMASS,AVPI) IF (Q02 .LE. 0.0) RETURN SIGMA =(Q2/Q02)**1.5/((SS-M0(i)**2)**2+M0(i)**2*GAMMA(i)**2)*10. END IF SIGMA = SIGMA*0.5*HC**2*PI/P2*ALPHA(I)*(PR2/P02)**(0.5*BETA(I)) 1 *M0(I)**2*GAMMA(I)**2 RETURN END ************************************************************************ REAL FUNCTION WIDTH(DMASS) * Function WIDTH(DMASS) gives the delta decay width for * given delta mass, the formula given by Kitazoe is used ************************************************************************ IMPLICIT NONE REAL DMASS REAL AVMASS, AVPI, QAVAIL PARAMETER (AVMASS=0.938919, AVPI=0.138037) QAVAIL = 0.25*(DMASS**2-AVMASS**2-AVPI**2)**2-(AVMASS*AVPI)**2 WIDTH = 0.0 IF (QAVAIL .LE. 0.) RETURN QAVAIL = QAVAIL/DMASS/DMASS WIDTH = 0.47*QAVAIL**1.5 /(AVPI**2 + 0.6*QAVAIL) RETURN END ************************************************************************ REAL FUNCTION XNPION(I1, I2, LA, XMAX) * PURPOSE : CALCULATE THE PION+NUCLEON CROSS SECTION ACCORDING TO THE * * BREIT-WIGNER FORMULA/(p*)**2 * * VARIABLE : LA = 1 FOR DELTA RESONANCE * * LA = 0 FOR N* RESONANCE * * DATE : JAN.29,1990 * * Refernce: Danielewicz et al. Nuc. Phys A533, (1991) 712 (4.17) * ************************************************************************ IMPLICIT NONE INTEGER I1, I2, LA REAL XMAX INTEGER MAXPAR, MAXP, MAXR REAL AVMASS, AVPI, AD0, AN0, PI, HC PARAMETER (MAXPAR=100000, MAXP=100, MAXR=1000) PARAMETER (AVMASS=0.938919, AVPI=0.138037, AD0=1.232, AN0=1.44, 1 PI=3.14159, HC=0.197327) REAL P, E COMMON /BB/ P(3,MAXPAR) COMMON /CC/ E(MAXPAR) REAL E10, E20, P1, P2, P3, DM, GAM2, PDELTA2, F1, PSTAR2, 1 WIDTH, AM02, FACTOR XNPION = 0.0 * 1. DETERMINE THE MOMENTUM COMPONENT OF DELTA IN THE LAB. FRAME E10 = SQRT(E(I1)**2 + P(1,I1)**2 + P(2,I1)**2 + P(3,I1)**2) E20 = SQRT(E(I2)**2 + P(1,I2)**2 + P(2,I2)**2 + P(3,I2)**2) P1 = P(1,I1) + P(1,I2) P2 = P(2,I1) + P(2,I2) P3 = P(3,I1) + P(3,I2) * 2. DETERMINE THE MASS OF DELTA BY USING OF THE REACTION KINEMATICS DM = SQRT((E10 + E20)**2 - P1**2 - P2**2 - P3**2) IF (DM.LE.(avmass+avpi)) RETURN * 3. DETERMINE THE PION+NUCLEON CROSS SECTION ACCORDING TO THE * BREIT-WIGNER FORMULA IN UNIT OF FM**2 IF (LA.EQ.1) THEN GAM2 = WIDTH(DM)**2 AM02 = AD0*AD0 PDELTA2 = 0.051622 FACTOR = 4.*PI*HC*HC*4./2.0 ELSE GAM2 = 0.2**2 AM02 = AN0*AN0 PDELTA2 = 0.157897 FACTOR = PDELTA2/10. ENDIF F1 = AM02*GAM2/((DM**2 - AM02)**2 + AM02*GAM2) PSTAR2 = ((DM**2 - AVMASS**2 + AVPI**2)/(2.*DM))**2 - AVPI**2 IF (PSTAR2.LE.0.) RETURN XNPION = F1/PSTAR2*FACTOR*XMAX RETURN END ************************************************************************ * * * cross-sections np-total, pp-elastic and pp-inelastic are fits * * to Particle Data Group 1988 Phys. Lett. B april 1988 P.126 * * * * Input: E_cm in GeV Output: cross-section in mB * * * ************************************************************************ real function nptotal(ss) implicit none integer n real avmass parameter (n=10, avmass=0.9383) real x, ss, para(n) data para /143.7252 , 122.8609 , 94.89669, 39.97737 , 1 28.68538 , 18.14496, -13.12878, -1.713403, 1 -0.4569092, 9.672062e-2/ if ( ss.gt.2.*avmass ) then x=0.5*log(ss*ss*((ss/2.0/avmass)**2 - 1.0)) nptotal = exp((para(1) + para(2)*x + para(3)*x**2 + 1 para(7)*x**3)/(para(4) + para(5)*x + para(6)*x**2 + 1 para(8)*x**3) + para(9)*x + para(10)*x*x) else nptotal=1e33 end if return end ************************************************************************ real function npinel(ss) implicit none real ss, sigma, nunuin if ( ss.le.2.0159 ) then npinel=0. else if ( ss.le.2.74 ) then npinel = 0.5*(sigma(1,1,ss) + sigma(3,1,ss) + 2.0*sigma(2,1,ss) 1 + 3.0*sigma(4,2,ss)) else npinel=nunuin(ss) end if return end ************************************************************************ real function ppinel(ss) implicit none integer n parameter (n=10) real ss, para(n), sig data para /-44.72405 , -3.342227 , 13.00564 , 3.588756 , 1 -3.493346 , 0.8526935, 0.17149462, 1.5889808e-3, 1 -0.2923224, -9.5134848e-5/ if ( ss.le.2.015 ) then sig=0. else sig = log(ss/2.015)*((para(1) + para(2)*ss + para(3)*ss**2 + 1 para(9)*ss**3)/(para(4) + para(5)*ss + para(6)*ss**2 + 1 para(10)*ss**3) + para(7)*ss + para(8)*ss**2) end if if ( sig.lt.0 ) sig=0. ppinel=sig return end ************************************************************************ real function ppelast(ss) implicit none integer n real avmass parameter (n=10, avmass=0.9383) real x, ss, para(n) data para /-3346.801 , -136.029 , -6437.515 , 2.5812924, 1 -200.7209, -51.38086, -0.3318678, -329.8205 , 1 -658.6221, 3.3176016e-2/ if ( ss.gt.2.*avmass ) then x=.5*log(ss*ss*((ss/2.0/avmass)**2 - 1.0)) ppelast = exp((para(1) + para(5)*x + para(8)*x**2)/(para(9)*x**3 1 + para(2)*x**2 + para(6)*x + para(3)) + para(4) + 1 para(7)*x + para(10)*x**2) else ppelast=1e33 end if return end ************************************************************************ * Cugnon parametrization of cross sections ************************************************************************ real function nunuel(ss,cutoff) implicit none real ss, cutoff if ( ss.le.cutoff ) then nunuel=55. else nunuel=35./(1.+100.*(ss-cutoff))+20. end if return end ************************************************************************ real function nunuin(ss) implicit none real ss if ( ss.lt.2.0159 ) then nunuin=0. else nunuin=20./(.015/(ss-2.0159)**2+1.) end if return end ************************************************************************ * purpose: functions to identify particles * ************************************************************************ logical function nucleon(lb1) implicit none integer lb1 if ( lb1.le.2 ) then nucleon=.true. else nucleon=.false. end if return end c----------------------------------------------------------------------- logical function pion(lb1) implicit none integer lb1 if ( lb1.le.5.and.lb1.ge.3 ) then pion=.true. else pion=.false. end if return end c----------------------------------------------------------------------- logical function delta(lb1) implicit none integer lb1 if ( lb1.le.9.and.lb1.ge.6 ) then delta=.true. else delta=.false. end if return end c----------------------------------------------------------------------- logical function resonance(lb1) implicit none integer lb1 if ( lb1.le.11.and.lb1.ge.6 ) then resonance=.true. else resonance=.false. end if return end