! ******************************************************** subroutine CANONICAL(pxc,pyc) ! ! Correct Px and Py for vector potential of solenoid field ! implicit double precision (a-h,o-z) include 'icommon.inc' ! pxc = pp(1) - 0.5d0 * pcharge * xp(2) * bfld(3) pyc = pp(2) + 0.5d0 * pcharge * xp(1) * bfld(3) ! return end ! ******************************************************** subroutine CLEAR c c Clear statistics arrays c implicit double precision (a-h,o-z) include 'icommon.inc' c do i=1,20 do j=1,5 st(j,i) = 0. end do st(6,i) = 1.e6 st(7,i) = -1.e6 end do c do i=1,20 do j=1,14 stsp(j,i) = 0. end do stsp(6,i) = 1.e6 stsp(7,i) = -1.e6 stsp(12,i) = 1.e6 stsp(13,i) = -1.e6 end do ! return end c ********************************************************* subroutine COVARIANCE ! ! Compute covariance matrices of particle distribution ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension q(6),r(6,6),s(6,6,6),cor2(6,6),cor3(6,6,6) dimension cov(6,6),sd(6),v(6) real*8 n character*2 lab(6) ! data lab/'x','y','ct','Px','Py','Pz'/ ! if( ncovar .eq. 0 ) go to 900 ! do 600 iz=1,ncovar ! Test is this is a desired covariance plane if( icvdes(iz) .ne. j7rec ) go to 600 ! do 50 i=1,6 q(i) = 0. do 50 j=1,6 r(i,j) = 0. do 50 k=1,6 s(i,j,k) = 0. 50 continue n = 0. ! do 100 ip=1,npart call READ_EVT(ip) ! if( ipflg .ne. 0 ) go to 100 ! v(1) = xp(1) v(2) = xp(2) v(3) = tp * cc v(4) = pp(1) v(5) = pp(2) v(6) = pp(3) ! if( refpar .ne. 0 ) then ! v(3) = (tp -tpref) * cc ! v(6) = pp(3) - ppref(3) ! end if ! ! Accumulate sums n = n + evtwt ! do i=1,6 q(i) = q(i) + v(i) * evtwt do j=1,6 r(i,j) = r(i,j) + v(i)*v(j) * evtwt do k=1,6 s(i,j,k) = s(i,j,k) + v(i)*v(j)*v(k) * evtwt end do end do end do ! 100 continue ! if( n .eq. 0. ) go to 900 ! ! Get expectation values do i=1,6 q(i) = q(i) / n ! 1st moment (mean) do j=1,6 r(i,j) = r(i,j) / n ! 2nd moment do k=1,6 s(i,j,k) = s(i,j,k) / n ! 3rd moment end do end do end do ! ! Get standard deviations do i=1,6 arg = r(i,i) - q(i)**2 if( arg .gt. 0. ) then sd(i) = SQRT( arg ) else sd(i) = -998. end if end do ! ! Convert moments to covariances and correlations do i=1,6 do j=1,6 cov(i,j) = r(i,j) - q(i)*q(j) fac = sd(i) * sd(j) if( fac .gt. 1e-10 ) then cor2(i,j) = cov(i,j) / fac else cor2(i,j) = -999. end if do k=1,6 fac = sd(i) * sd(j) * sd(k) if( fac .gt. 1e-10 ) then cor3(i,j,k) = s(i,j,k) - q(i)*r(j,k) - q(j)*r(i,k) & - q(k)*r(i,j) + 2.*q(i)*q(j)*q(k) cor3(i,j,k) = cor3(i,j,k) / fac else cor3(i,j,k) = -999. end if end do end do end do ! ! Compute dispersions if( sd(6) .gt. 0. ) then disx = cor2(1,6) * sd(1) * q(6) / sd(6) disy = cor2(2,6) * sd(2) * q(6) / sd(6) disxp = cor2(4,6) * sd(4) / sd(6) disyp = cor2(5,6) * sd(5) / sd(6) else disx = -999. disy = -999. disxp = -999. disyp = -999. end if ! write(2,*) write(2,110) j7rec 110 format(15('*'),t18,' PARTICLE COVARIANCE MATRIX, PLANE ',i5, & t63,15('*')) write(2,115) 115 format(t12,'Mean values') do i=1,6 if( i .eq. 3 ) then write(2,118) lab(i),q(i),'c*tpref',cc*tpref else if( i .eq. 6 ) then write(2,118) lab(i),q(i),'pzref',ppref(3) else write(2,118) lab(i),q(i) 118 format(t5,a2,t10,e13.6,t30,a7,t38,e13.6) end if end do ! write(2,*) write(2,120) 120 format(t10,'SD',t30,'Second order correlations') write(2,130) (lab(i),i=1,5) 130 format(t25,a2,t35,a2,t45,a2,t55,a2,t65,a2) ! do i=1,6 write(2,150) lab(i),sd(i),(cor2(i,j),j=1,i-1) 150 format(a3,e12.4,2x,5f10.3) end do ! write(2,*) write(2,'(a)')' Third order correlations' write(2,160) (lab(i),i=1,6) 160 format(t14,a2,t24,a2,t34,a2,t44,a2,t54,a2,t64,a2) ! do i=1,6 do j=i,6 if( i.eq. 1) then write(2,'(2a3,6f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j) else if( i .eq. 2 ) then write(2,'(2a3,10x,5f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j) else if( i .eq. 3 ) then write(2,'(2a3,20x,4f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j) else if( i .eq. 4 ) then write(2,'(2a3,30x,3f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j) else if( i .eq. 5 ) then write(2,'(2a3,40x,2f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j) else if( i .eq. 6 ) then write(2,'(2a3,50x,f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j) end if end do write(2,*) end do ! write(2,'(a)')' Dispersions [m]' write(2,165) 'D(X)= ',disx,'D(Y)= ',disy,"D'(X)= ", & disxp,"D'(Y)= ",disyp 165 format(a, e12.5, t20, a, e12.5, t40, a, e12.5, t60, a, e12.5) ! 600 continue ! end of covariance plane loop ! 900 continue end c ********************************************************* subroutine DISP_CORR(n,emx,emy,xarg,yarg) ! ! Compute 2D emittances in presence of dispersion ! Uses correction scheme of D. Neuffer ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension emx(5),emy(5),emd(6),evx(5),evy(5) real*8 n ! do i=1,6 emd(i) = 0. end do ! do 100 ip=1,npart ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 100 ! pzc = pp(3) x = xp(1) y = xp(2) if( pxycorr ) then call CANONICAL(pxc,pyc) else pxc = pp(1) pyc = pp(2) end if ! ! cpx = pxc / pzc ! x' cpx = pxc ! cpy = pyc / pzc ! y' cpy = pyc ! emd(1) = emd(1) + evtwt * pzc emd(2) = emd(2) + evtwt * pzc * pzc emd(3) = emd(3) + evtwt * x * pzc emd(4) = emd(4) + evtwt * cpx * pzc emd(5) = emd(5) + evtwt * y * pzc emd(6) = emd(6) + evtwt * cpy * pzc ! 100 continue ! ! Get expectation values do i=1,5 evx(i) = emx(i) / n evy(i) = emy(i) / n end do do i=1,6 emd(i) = emd(i) / n end do ! vp = emd(2) - emd(1)*emd(1) if( vp .gt. 0. ) then ! get dispersions fac = emd(1) / vp dx = (emd(3)-evx(1)*emd(1)) * fac dxp = (emd(4)-evx(3)*emd(1)) * fac dy = (emd(5)-evy(1)*emd(1)) * fac dyp = (emd(6)-evy(3)*emd(1)) * fac else xarg = -1. ! flag failure yarg = -1. go to 900 end if ! ! Correct the x and y distributions fac = vp / emd(1) / emd(1) evx(2) = evx(2) - dx*dx*fac evx(4) = evx(4) - dxp*dxp*fac evx(5) = evx(5) - dx*dxp*fac evy(2) = evy(2) - dy*dy*fac evy(4) = evy(4) - dyp*dyp*fac evy(5) = evy(5) - dy*dyp*fac ! xarg = evx(2)*evx(4) - evx(5)*evx(5) yarg = evy(2)*evy(4) - evy(5)*evy(5) ! 900 continue end c ********************************************************* subroutine EMITCOV ! ! Compute covariance emittances of beam at this plane ! implicit double precision (a-h,o-z) include 'icommon.inc' dimension q(6),r(6,6) dimension cov(6,6),v(6) real*8 n integer indx(6) ! do 600 iz=1,nemit ! Test if this is one of the desired emittance planes if( iemdes(iz) .ne. j7rec ) go to 600 ! do 50 i=1,6 q(i) = 0. do 50 j=1,6 r(i,j) = 0. 50 continue n = 0. ! do 100 ip=1,npart call READ_EVT(ip) ! if( ipflg .ne. 0 ) go to 100 ! v(1) = xp(1) v(2) = xp(2) v(3) = tp * cc v(4) = pp(1) v(5) = pp(2) v(6) = pp(3) ! ! Accumulate sums n = n + evtwt ! do i=1,6 q(i) = q(i) + v(i) * evtwt do j=1,6 r(i,j) = r(i,j) + v(i)*v(j) * evtwt end do end do ! 100 continue ! if( n .eq. 0. ) go to 900 ! ! Get expectation values do i=1,6 q(i) = q(i) / n ! 1st moment (mean) do j=1,6 r(i,j) = r(i,j) / n ! 2nd moment end do end do ! ! Convert moments to covariances do i=1,6 do j=1,6 cov(i,j) = r(i,j) - q(i)*q(j) end do end do ! ! Get LBNL value of 6-D emittance call LUDCMP(cov,6,6,indx,d) ! if( d .eq. 0. ) then ! flag singular matrix ! d = -998. ! go to 580 go to 900 end if ! do i=1,6 d = d * cov(i,i) end do if( d .gt. 0. ) then d = SQRT( d ) else ! d = -998. go to 900 end if ! !580 continue ! Normalization = 1 / (mc)^3 ! d = d / pmass**3 ! write(2,*) ! write(2,'(a,e13.6)') '6-D emittance [m^3] = ',d write(2,*) ! call LBNL_EMIT(iz) ! 600 continue ! 900 continue end c ********************************************************* subroutine EMIT2D ! ! Compute 2D emittances of beam at this plane ! including possible tail cuts and Bz corrections ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension emx(5),emy(5),emz(5),stem(2,6) dimension xmean(6),xstand(6),ystand(6) real*8 n integer ncut(6) ! data iteste/20/,sig_thre/1e-4/ ! if( nemit .eq. 0 ) go to 900 ! --------------------------------------------------------- do 600 iz=1,nemit ! Test if this is one of the desired emittance planes if( iemdes(iz) .ne. j7rec ) go to 600 ! iflg = 0 ! ......................................................... ! Tail cutting logic if( sigmacut ) then knt = 0 ! count iterations do i=1,6 ystand(i) = 1e6 ncut(i) = 10 end do ! 100 continue ! enter here for iterations on tail cut ! knt = knt + 1 if( knt .gt. iteste ) then iflg = 1 go to 300 end if ! do j=1,6 ! zero array for {x...Pz} statistics do i=1,2 stem(i,j) = 0. end do end do ! stev = 0. !.......................................................... ! Read back particle information for this plane do 200 ip=1,npart ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 200 ! if( knt .eq. 1 ) emevt(iz) = emevt(iz) + evtwt ! ! Correct Px and Py for vector potential if( pxycorr ) then call CANONICAL(pxc,pyc) else pxc = pp(1) pyc = pp(2) end if ! energy = SQRT( pxc**2 + pyc**2 + pp(3)**2 + pmass**2 ) z = pp(3) / energy * cc * tp ! stev = stev + evtwt stem(1,1) = stem(1,1) +evtwt* xp(1) stem(2,1) = stem(2,1) +evtwt* xp(1)**2 stem(1,2) = stem(1,2) +evtwt* xp(2) stem(2,2) = stem(2,2) +evtwt* xp(2)**2 stem(1,3) = stem(1,3) +evtwt* z stem(2,3) = stem(2,3) +evtwt* z*z stem(1,4) = stem(1,4) +evtwt* pxc stem(2,4) = stem(2,4) +evtwt* pxc**2 stem(1,5) = stem(1,5) +evtwt* pyc stem(2,5) = stem(2,5) +evtwt* pyc**2 stem(1,6) = stem(1,6) +evtwt* pp(3) stem(2,6) = stem(2,6) +evtwt* pp(3)**2 ! 200 continue ! end of particle loop ! ......................................................... if( stev .eq. 0. ) then iflg = 2 go to 300 end if ! call MEANSD(stev,stem,xmean,xstand) ! ! Test if tail cutting process has converged ! Stop if fractional change in sigma < sig_thre and ncuts=0 do 250 i=1,6 denom = xstand(i) if( denom .le. 0. ) go to 250 ! arg = ABS(ystand(i) - xstand(i)) / denom if( arg.gt.sig_thre .or. ncut(i).ne.0 ) then ! If any of the 6 variables requires a tail cut, load and cut all 6 do j=1,6 ystand(j) = xstand(j) end do ! call SIGMA_CUT(xmean,xstand,ncut) ! cut tails on {x...Pz} ! go to 100 end if ! 250 continue end if ! end of test on sigmacut ! --------------------------------------------------------- 300 continue ! ! Accumulate final statistics for emittance calculation ! ! We do x and y first because we may need them for corrections ! to the z emittance ! do i=1,5 emx(i) = 0.d0 emy(i) = 0.d0 end do ! embg = 0.d0 ! ......................................................... do 350 ip=1,npart ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 350 emevcut(iz) = emevcut(iz) + evtwt if( .not. sigmacut ) emevt(iz) = emevt(iz) + evtwt ! if( pxycorr ) then call CANONICAL(pxc,pyc) else pxc = pp(1) pyc = pp(2) end if ! pzc = pp(3) x = xp(1) y = xp(2) ! cpx = pxc / pzc cpx = pxc ! cpy = pyc / pzc cpy = pyc pm = SQRT( pxc**2 + pyc**2 + pzc**2 ) ! embg = embg + evtwt* pm / pmass ! emx(1) = emx(1) + evtwt* x emx(2) = emx(2) + evtwt* x * x emx(3) = emx(3) + evtwt* cpx emx(4) = emx(4) + evtwt* cpx * cpx emx(5) = emx(5) + evtwt* x * cpx ! emy(1) = emy(1) + evtwt* y emy(2) = emy(2) + evtwt* y * y emy(3) = emy(3) + evtwt* cpy emy(4) = emy(4) + evtwt* cpy * cpy emy(5) = emy(5) + evtwt* y * cpy ! 350 continue ! ......................................................... ! Calculate x and y normalized emittances ! For statistical definition, cf. J. Buon, CERN 91-04, p. 48 ! emitx(iz) = -999. emity(iz) = -999. n = emevcut(iz) if( n .lt. 2 ) go to 550 betgam = embg / n ! mean value of beta * gamma ! if( discorr ) then call DISP_CORR(n,emx,emy,xarg,yarg) if( xarg .gt. 0. )then ! emitx(iz) = SQRT( xarg ) * betgam emitx(iz) = 1.d0/pmass * SQRT( xarg ) end if if( yarg .gt. 0. )then ! emity(iz) = SQRT( yarg ) * betgam emity(iz) = 1.d0/pmass * SQRT( yarg ) end if ! else ! no dispersion correction qx = emx(1) / n vx = emx(2)/n - qx*qx if(vx .lt. 0.) go to 410 qpx = emx(3) / n vpx = emx(4)/n - qpx*qpx if(vpx .lt. 0.) go to 410 cx = emx(5) / n - qx*qpx arg = vx * vpx - cx * cx ! if( arg .gt. 0. )then ! emitx(iz) = SQRT( arg ) * betgam emitx(iz) = 1.d0/pmass * SQRT( arg ) end if c 410 continue qx = emy(1) / n vx = emy(2)/n - qx*qx if(vx .lt. 0.) go to 420 qpy = emy(3) / n vpx = emy(4)/n - qpy*qpy if(vpx .lt. 0.) go to 420 cx = emy(5) / n - qx*qpy arg = vx * vpx - cx * cx ! if( arg .gt. 0. )then ! emity(iz) = SQRT( arg ) * betgam emity(iz) = 1.d0/pmass * SQRT( arg ) end if end if ! end of test on discorr ! --------------------------------------------------------- 420 continue ! Now accumulate final statistics for z emittance calculation ! if( ipzcor .ne. 0 ) then ! get Pz-amp correlation value if( emitx(iz).eq.0. .or. emity(iz).eq.0. ) go to 550 betax = betgam * emx(2)/n / emitx(iz) betay = betgam * emy(2)/n / emity(iz) if( betax.gt.0. .and. betay.gt.0. ) then ! call PZ_AMP_COR(betax,betay,slp) ! end if end if ! emitz(iz) = -999. do i=1,5 emz(i) = 0.d0 end do ! ......................................................... do 450 ip=1,npart ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 450 ! if( pxycorr ) then call CANONICAL(pxc,pyc) else pxc = pp(1) pyc = pp(2) end if ! if( ipzcor .eq. 1 ) then ! use corrected Pz for emittance amp = SQRT( (pxc/pp(3))**2 + (pyc/pp(3))**2 + & (xp(1)/betax)**2 + (xp(2)/betay)**2 ) pzc = pp(3) - slp*amp ! else if( ipzcor .eq. 2 ) then r2 = xp(1)**2 + xp(2)**2 pt2 = pxc**2 + pyc**2 pm2 = pmass**2 amp = pt2/pm2 + ((pcharge*bfld(3))**2*r2)/(4.*pm2) amp = SQRT(amp) pzc = pp(3) - slp*amp ! else pzc = pp(3) end if ! energy = SQRT( pxc**2 + pyc**2 + pzc**2 + pmass**2 ) ! if( refpar.ne.0 .and. diagref) then ! with reference particle cpz = pzc - ppref(3) z = pp(3)/energy*cc * (tp - tpref) else ! no reference particle cpz = pzc z = pp(3)/energy*cc * tp end if ! emz(1) = emz(1) + evtwt* z emz(2) = emz(2) + evtwt* z * z emz(3) = emz(3) + evtwt* cpz emz(4) = emz(4) + evtwt* cpz * cpz emz(5) = emz(5) + evtwt* z * cpz ! 450 continue ! ......................................................... ! Calculate normalized z emittance ! emitz(iz) = -999. n = emevcut(iz) if( n .lt. 2 ) go to 550 ! qx = emz(1) / n ! qx = vx = emz(2)/n - qx*qx if(vx .lt. 0.) go to 550 qpz = emz(3) / n vpx = emz(4)/n - qpz*qpz if(vpx .lt. 0.) go to 550 cx = emz(5) / n - qx*qpz arg = vx * vpx - cx * cx ! if( arg .gt. 0.)then emitz(iz) = 1.d0/pmass * SQRT( arg ) end if c ! --------------------------------------------------------- 550 continue if( sigmacut ) then if( iflg .ne. 0 ) then ! flag problem cases if(emitx(iz) .gt. 0.) emitx(iz) = -emitx(iz) if(emity(iz) .gt. 0.) emity(iz) = -emity(iz) if(emitz(iz) .gt. 0.) emitz(iz) = -emitz(iz) end if ! ! Reset particle flag for particles in tails do ip=1,npart ! call READ_EVT(ip) ! if( ipflg .eq. 10 ) then ! tail particle ipflg = 0 call WRITE_EVT(ip) end if ! end do end if ! 600 continue ! end of emittance plane loop ! 900 continue return end c ******************************************************** subroutine FILL_HIST(izp) C c Fills histogram arrays C implicit double precision (a-h,o-z) include 'icommon.inc' ! Common variables: var c if( nhist .eq. 0 ) go to 900 call LOADVAR( izp ) c DO 110 I=1,NHISt if(ihdes(i) .ne. izp) go to 110 val = var( ihvar(i) ) call STATH(val,i) ! accumulate statistics vhi = hxmin(i) + nhbins(i)*hdx(i) ! if(val .lt. hxmin(i))then ihund(i) = ihund(i) + 1 go to 110 else if(val .gt. vhi)then ihovf(i) = ihovf(i) + 1 go to 110 else ihcon(i) = ihcon(i) + 1 ip = lhpnt(i) + int( (val-hxmin(i))/hdx(i) ) hval(ip) = hval(ip) + 1 end if ! 110 continue c 900 continue return end c ********************************************************* subroutine FILL_POL(izp) ! ! Fill polarization arrays at plane IZP ! implicit double precision (a-h,o-z) include 'icommon.inc' ! do 100 i=1,nemit if( iemdes(i) .ne. izp ) go to 100 ! if( izp .ge. 0 ) then ! normal plane if( jpar.gt.0 .and. ABS(iptyp).eq.2 .and. spintrk.gt.0 ) then px = pol(1) py = pol(2) pz = pol(3) hel = HELICITY( pol ) else px = 0. py = 0. pz = 0. hel = pol(3) ! DECAY sets this when spintrk=F end if ! else ! production if( jpar.gt.0 .and. ABS(iptyp).eq.2 ) then px = polinit(1) py = polinit(2) pz = polinit(3) hel = HELICITY( polinit ) end if end if ! polx(i) = polx(i) + px*evtwt poly(i) = poly(i) + py*evtwt polz(i) = polz(i) + pz*evtwt helic(i) = helic(i) + hel*evtwt 100 continue ! return end c ******************************************************** SUBROUTINE FILL_SCAT(izp) C c Fills scatterplot arrays C implicit double precision (a-h,o-z) include 'icommon.inc' ! Common variables: var c ! dimension QST(20) c IF(nscat .EQ. 0) GO TO 900 c call LOADVAR( izp ) c DO 200 I=1,Nscat if(isxdes(i).ne.izp .and. isydes(i).ne.izp) go to 200 ! idx = isxdes(i) ! idy = isydes(i) c c check if one of the plot quantities is downstream c if so, save in qst array ! if(idx .gt. izp)then ! x variable coming later ! qst(i) = var( isyvar(i) ) ! save y variable from this place ! go to 200 ! end if ! ! if(idy .gt. izp)then ! y variable coming later ! qst(i) = var( isxvar(i) ) ! save x variable from this place ! go to 200 ! end if c c check if one of the plot quantities was upstream ! if(idx .lt. izp)then ! previously saved the x variable ! vx = qst(i) ! vy = var( isyvar(i) ) ! else if(idy .lt. izp)then ! previously saved the y variable ! vx = var( isxvar(i) ) ! vy = qst(i) ! else ! both quantites are given at this plane vx = var( isxvar(i) ) vy = var( isyvar(i) ) ! end if c call STATSP(vx,vy,i) ! accumulate statistics ! c check limits vxhi = sxmin(i) + nsxbin(i)*sdx(i) if(vx .lt. sxmin(i))then isxund(i) = isxund(i) + 1 go to 200 else if(vx .gt. vxhi)then isxovf(i) = isxovf(i) + 1 go to 200 end if c vyhi = symin(i) + nsybin(i)*sdy(i) if(vy .lt. symin(i))then isyund(i) = isyund(i) + 1 go to 200 else if(vy .gt. vyhi)then isyovf(i) = isyovf(i) + 1 go to 200 end if c iscon(i) = iscon(i) + 1 ip = int( (vy-symin(i))/sdy(i) ) * nsxbin(i) ip = ip + int( (vx-sxmin(i))/sdx(i) ) + lspnt(i) sval(ip) = sval(ip) + 1 200 continue c 900 continue return end c ******************************************************** subroutine FILL_ZHIST ! ! Fill z-history arrays ! implicit double precision (a-h,o-z) include 'icommon.inc' ! Common variables: var ! Keep zhxval a function of particle # for print purposes ! call LOADVAR( 1 ) zp = xp(3) ! do i=1,nzhist if( jpar .le. ABS(nzhpar(i)) ) then if( zp.ge.zhxcur(jpar,i) .and. zp.lt.zhxmax(i) ) then ! izhknt(jpar,i) = izhknt(jpar,i) + 1 knt = izhknt(jpar,i) if( knt .gt. 70 ) knt = 70 zhxval( knt,jpar,i ) = zp zhyval( knt,jpar,i ) = var( izhyvar(i) ) ! update current pointer zhxcur(jpar,i) = zhxcur(jpar,i) + zhdx(i) end if ! end if end do ! return end ! ********************************************************* subroutine LBNL_EMIT(iz) ! ! Compute emittances, covariances, and related quantities ! following the wisdom of the LBNL group ! ! Adapted from G. Penn, 3 Nov 99 ! ! if reference particle used, ! only particles of same type used in calculations ! Sigma cuts implemented, but different from EMITTANCE ! "pzcorr" only affects longitudinal emittance calculation ! "pxycorr" has no effect ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension xtemp(3) dimension w(6),moments(6,6),mcovar(6,6),mtrans(4,4),avg(6) dimension xpc(3),ppc(3),mcovar2(6,6),sdev(6) dimension avglong(3),mlong(3,3),mclong(3,3),mc2(3,3) real*8 kappa,moments,mcovar,mtrans,mcovar2 real*8 Lmech,Lcanon,Ldim,mlong,mclong,mc2 character*8 axes(6) integer indperp(4),ishift(4),indx(6),indlong(3) integer usepar(mxpar),knt,iteste,ncut logical fset,longcorr,loopflag ! LOGICAL FIRST ! data axes/'x','y','tp-tpref','Px','Py','Ep-Epref'/ data iteste/20/ ! DATA FIRST/.TRUE./ ! ! IF( FIRST ) THEN ! OPEN(UNIT=18,FILE='MTRANS.DAT',STATUS='UNKNOWN') ! FIRST = .FALSE. ! END IF fset = .false. ! if( refpar .ne. 0 ) then call READ_EVT(0) refenergy = DSQRT(ppref(3)**2 + pmass**2) reftime = tpref zloc = xp(3) call FIELD(xp,tp) bzaxis = bfld(3) charge = pcharge fset = .true. else refenergy = 0. reftime = 0. bzaxis = 0. charge = 0. zloc = 0. end if ! ishift(1) = 1 ishift(2) = 2 ishift(3) = 4 ishift(4) = 5 ! dtrans = 0. dlong = 0. ! ! longcorr = pzcorr if( ipzcor .ne. 0 ) then longcorr = .true. else longcorr = .false. end if ! do ip=1,npart usepar(ip)=0 end do ! loopflag = .true. ! test for iterations knt = 0 ! count iterations ! ! --------------------------------------------------------- ! do while (loopflag) ! ! Particle cutting logic ncut = 0 do 200 ip=1,npart if (usepar(ip) .lt. 0) go to 200 ! call READ_EVT(ip) ! energy in GeV energy = DSQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + pmass**2) if(ipflg .ne. 0) usepar(ip)=-1 if((refpar .ne. 0) .and. & (ABS(iptyp) .ne. iptypref)) usepar(ip)=-2 ! RCF 19 APR 2010 if (usepar(ip) .lt. 0) go to 200 ! if (.not. fset) then ! find z and Bz zloc = xp(3) xtemp(1)=0. xtemp(2)=0. xtemp(3)=xp(3) call FIELD(xtemp,tp) bzaxis = bfld(3) charge = pcharge fset = .true. end if ! if (sigmacut) then ! check if amplitudes too large Atrans = 0. if (dtrans .gt. 0.) then kappa = 0.5 * charge * bzaxis / avgpz do i=1,3 xpc(i)=xp(i)-avg(i) ppc(i)=pp(i)-avg(i+3) end do Atrans = (1/pmass)*( (csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2) + & (xpc(1)**2 + xpc(2)**2)* avgpz * csgamma + & 2*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) + & 2*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) ) if (Atrans .gt. 2 * sig_cut**2 * dtrans) then ! factor of two to make up for combining 4D at once usepar(ip)=-3 ncut = ncut + 1 go to 200 end if end if if (dlong .gt. 0.) then ctime = tp - reftime - avglong(1) cenergy = energy - refenergy - avglong(2) if (longcorr) then ctime = ctime - corr_t * (Atrans-avglong(3)) cenergy = cenergy - corr_E * (Atrans-avglong(3)) end if Along = (cc/pmass) * ( ctime**2/delta + & delta * (cenergy - alphalong * ctime / delta)**2 ) if (Along .gt. sig_cut**2 * dlong) then usepar(ip)=-4 ncut = ncut + 1 go to 200 end if end if end if 200 continue ! end of particle loop ! ! calculate beam statistics ! initialize avgpz = 0. do 250 i=1,6 avg(i) = 0. sdev(i) = 0. do 250 j=1,6 moments(i,j) = 0. 250 continue stev = 0. ! first loop do 300 ip=1,npart if (usepar(ip) .lt. 0) go to 300 ! call READ_EVT(ip) w(1) = xp(1) w(2) = xp(2) w(3) = tp - reftime w(4) = pp(1) w(5) = pp(2) ! energy in GeV energy = DSQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + pmass**2) w(6) = energy - refenergy ! ! Accumulate sums stev = stev + evtwt avgpz = avgpz + pp(3) * evtwt do i=1,6 avg(i) = avg(i) + w(i) * evtwt do j=1,6 moments(i,j) = moments(i,j) + w(i)*w(j) * evtwt end do end do ! 300 continue ! end of first particle summation ! if( stev .eq. 0. ) go to 900 ! ! Get expectation values avgpz = avgpz / stev kappa = 0.5 * charge * bzaxis / avgpz do i=1,6 avg(i) = avg(i) / stev ! 1st moment (mean) do j=1,6 moments(i,j) = moments(i,j) / stev ! 2nd moment end do end do ! ! Convert moments to covariances and standard deviations do i=1,6 do j=1,6 mcovar(i,j) = moments(i,j) - avg(i)*avg(j) mcovar2(i,j) = mcovar(i,j) end do if( mcovar(i,i) .lt. 0. ) go to 900 ! RCF 22 Nov 99 sdev(i)=DSQRT(mcovar(i,i)) end do ! ! fill mtrans with x, y, Px, Py data (corr to index 1,2,4,5) do i=1,4 do j=1,4 mtrans(i,j) = mcovar(ishift(i),ishift(j)) end do end do ! ! get 6-D emittance call LUDCMP(mcovar,6,6,indx,dfull) ! if( dfull .eq. 0. ) then ! flag singular matrix dfull = -998. go to 680 end if ! do i=1,6 dfull = dfull * mcovar(i,i) end do if( dfull .gt. 0. ) then dfull = DSQRT( dfull ) else dfull = -998. go to 680 end if ! Normalization = c / ( (mc)^2 (mc^2)) dfull = dfull * cc / pmass**3 ! third pmass is actually mc^2 in MeV 680 continue ! ! get transverse emittance (square root of 4D emittance) ! IF( J7REC.EQ.3 .OR. J7REC.EQ.52 ) THEN ! WRITE(18,'(I5)')J7REC ! DO I=1,4 ! WRITE(18,'(4E14.6)') (MTRANS(I,J),J=1,4) ! END DO ! END IF call LUDCMP(mtrans,4,4,indperp,dtrans) ! if( dtrans .eq. 0. ) then ! flag singular matrix dtrans = -998. go to 780 end if ! ! IF( J7REC.EQ.3 .OR. J7REC.EQ.52 ) THEN ! WRITE(18,'(I5)')J7REC ! DO I=1,4 ! WRITE(18,'(4E14.6)') (MTRANS(I,J),J=1,4) ! END DO ! END IF do i=1,4 dtrans = dtrans * mtrans(i,i) end do if( dtrans .gt. 0. ) then ! start with 4D transverse emittance dtrans = DSQRT( dtrans ) ! take another square root to get pseudo 2D emittance dtrans = DSQRT( dtrans ) else dtrans = -998. go to 780 end if ! ! Normalization = 1 / (mc) dtrans = dtrans / pmass 780 continue ! ! Lmech = moments(1,5) - moments(2,4) Lcanon = Lmech + 0.5 * charge * bzaxis & * (moments(1,1) + moments(2,2)) if( dtrans.gt.0. ) then Ldim = (mcovar2(1,5) - mcovar2(2,4) + & 0.5 * charge * bzaxis * (mcovar2(1,1) + mcovar2(2,2))) & / ( 2 * pmass * dtrans) csalpha = - (mcovar2(1,4) + mcovar2(2,5)) & / (2 * pmass * dtrans) csbeta = (mcovar2(1,1) + mcovar2(2,2)) * avgpz & / (2 * pmass * dtrans) csgamma = (mcovar2(4,4) + mcovar2(5,5)) & / (2 * pmass * dtrans * avgpz) altgamma = (1 / csbeta) * (1 + csalpha**2 & + (csbeta * kappa - Ldim)**2) else longcorr = .false. Ldim = 0. csbeta = 0. csalpha = 0. csgamma = 0. altgamma = 0. end if ! ! longitudinal beam parameters ! initialize do 350 i=1,3 avglong(i) = 0. do 350 j=1,3 mlong(i,j) = 0. 350 continue stev = 0. ! ! second loop do 400 ip=1,npart if (usepar(ip) .lt. 0) go to 400 ! call READ_EVT(ip) w(1) = tp - reftime energy = DSQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + pmass**2) w(2) = energy - refenergy if (longcorr) then do i=1,3 xpc(i)=xp(i)-avg(i) ppc(i)=pp(i)-avg(i+3) end do Atrans = (1/pmass)*( (csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2) + & (xpc(1)**2 + xpc(2)**2)* avgpz * csgamma + & 2*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) + & 2*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) ) w(3) = Atrans else w(3) = 0. end if ! ! Accumulate sums stev = stev + evtwt do i=1,3 avglong(i) = avglong(i) + w(i) * evtwt do j=1,3 mlong(i,j) = mlong(i,j) + w(i)*w(j) * evtwt end do end do ! 400 continue ! if( stev .eq. 0. ) go to 900 ! ! Get expectation values do i=1,3 avglong(i) = avglong(i) / stev ! 1st moment (mean) do j=1,3 mlong(i,j) = mlong(i,j) / stev ! 2nd moment end do end do ! ! Convert moments to covariances and standard deviations do i=1,3 do j=1,3 mclong(i,j) = mlong(i,j) - avglong(i)*avglong(j) mc2(i,j) = mclong(i,j) end do if( mclong(i,i) .lt. 0. ) go to 900 end do ! ! check some amplitude spread; otherwise, no correlation sigmaA = DSQRT(mclong(3,3)) if (sigmaA .eq. 0.) longcorr = .false. if (longcorr) then ! get 3-D emittance call LUDCMP(mclong,3,3,indlong,dlong) ! if( dlong .eq. 0. ) then ! flag singular matrix dlong = -998. go to 880 end if ! do i=1,3 dlong = dlong * mclong(i,i) end do if( dlong .gt. 0. ) then dlong = DSQRT( dlong / mc2(3,3)) else dlong = -998. go to 880 end if else dlong = mclong(1,1)*mclong(2,2)-mclong(1,2)**2 if( dlong .gt. 0. ) then dlong = DSQRT(dlong) else dlong = -998. go to 880 end if end if ! normalize to c/(mc^2) dlong = dlong * cc / pmass ! pmass is actually mc^2 in MeV 880 continue ! if(dlong .gt. 0.) then if (longcorr) then alphalong = (mc2(1,2) - mc2(1,3)*mc2(2,3)/mc2(3,3)) & / (dlong * pmass / cc) delta = (mc2(1,1) - mc2(1,3)**2/mc2(3,3)) & / (dlong * pmass / cc) corr_t = mc2(1,3)/mc2(3,3) corr_E = mc2(2,3)/mc2(3,3) else alphalong = mc2(1,2) / (dlong * pmass / cc) delta = mc2(1,1) / (dlong * pmass / cc) corr_t = 0. corr_E = 0. end if else alphalong = 0. delta = 0. corr_t = 0. corr_E = 0. end if ! ! check whether should repeat loop knt = knt + 1 if ((.not. sigmacut) .or. (knt .gt. iteste)) then loopflag=.false. else ! Test if tail cutting process has converged ! Stop if ncut = 0 (need to test changes in parameters?) if ((ncut .eq. 0) .and. (knt .gt. 1)) loopflag = .false. end if ! end do ! ! ! --------------------------------------------------------- ! OUTPUTS if( edetail ) then ! give more detailed information write(2,*) write(2,210) j7rec 210 format(15('*'),t18,' PARTICLE MOMENTS, PLANE ',i5, & t63,15('*')) write(2,212) 212 format(t17,'Mean values',t35,'Standard Deviation') do i=1,6 write(2,214) axes(i),avg(i),sdev(i) 214 format(t5,a8,t16,e13.6,t34,e13.6) end do write(2,216) 216 format(t2,' Reference Particle:',t30,'time',t45,'Energy', & t60,'Pz') write(2,218) reftime,refenergy,ppref(3) 218 format(t24,e15.8,t39,e15.8,t54,e15.8) write(2,*) ! write(2,310) 310 format(t2,'Emittances: ',t22,'6D',t30,'Transverse', & t42,'Corrected Longitudinal') write(2,320) dfull,dtrans,dlong 320 format(t16,3e12.4) write(2,'(a,e12.4)')' eL*eT*eT = ',dtrans**2 * dlong write(2,330) 330 format(t2,'Angular Momentum: ',t21,'Mechanical',t33,'Canonical', & t45,'Dimensionless = Lcanon/2mceT') write(2,340) Lmech,Lcanon,Ldim 340 format(t20,3e12.4) write(2,*) write(2,420) 420 format(t10,'Modified Courant-Snyder Parameters: ', & 'ASSUMES A CYLINDRICAL BEAM') write(2,432) 432 format(t5,'beta_perp',t17,'alpha_perp',t29,'gamma_perp') write(2,434) csbeta,csalpha,csgamma 434 format(t3,3e12.4) write(2,436) altgamma 436 format(t5,'consistency check:',t27,e12.4) write(2,*) write(2,438) 438 format('Long:',t7,'alphaL',t22,'delta',t37,'corr_t',t52,'corr_E', & t67,'sigma Atrans') write(2,440) alphalong,delta,corr_t,corr_E,sigmaA 440 format(5e15.4) write(2,*) write(2,442) 442 format('MISC:',t7,'Z',t22,'N (weighted)',t37,'Bz on axis', & t52,'kappa',t67,'') write(2,444) zloc,stev,bzaxis,kappa,avgpz 444 format(5e15.4) write(2,*) end if ! ! Save some quantities for a summary table emit6(iz) = dfull emit4(iz) = dtrans emitl(iz) = dlong canang(iz) = Lcanon betper(iz) = csbeta ! 900 continue end ! ********************************************************* subroutine LOADVAR(izp) c c Load kinematic variables for diagnostics c implicit double precision (a-h,o-z) include 'icommon.inc' ! if( izp .le. -1 ) then do i=1,3 var(i) = xpint(i) var(i+3) = ppint(i) end do var(7) = tpint * cc else do i=1,3 var(i) = xp(i) var(i+3) = pp(i) end do var(7) = tp * cc if( refpar.ne.0 .and. diagref ) then ! redefine pz and t var(6) = var(6) - ppref(3) var(7) = var(7) - tpref*cc end if end if ! call VMAG(pp(1),var(8)) var(9) = SQRT(var(8)**2 + pmass**2) ! E gam = var(9) / pmass var(10) = pmass * (gam - 1.) ! KE var(11) = pp(1) / pp(3) ! x' var(12) = pp(2) / pp(3) ! y' var(13) = ACOS( pp(3) / var(8) ) ! space angle var(14) = SQRT( xp(1)**2 + xp(2)**2 ) ! r ! if( .not.(xp(1).eq.0. .and. xp(2).eq.0.) ) then phi = ATAN2( xp(2),xp(1) ) ! phi else phi = 0. end if ! var(15) = phi cph = COS( phi ) sph = SIN( phi ) var(16) = pp(1)*cph + pp(2)*sph var(17) = pp(2)*cph - pp(1)*sph ! x = xp(1) vlx = xp(2)*pp(3) - xp(3)*pp(2) vly = xp(3)*pp(1) - x*pp(3) vlz = x*pp(2) - xp(2)*pp(1) vl2 = vlx*vlx + vly*vly + vlz*vlz var(18) = vlz ! Lz var(19) = vl2 ! L2 var(20) = sarc ! var(21) = SQRT(pp(1)**2 + pp(2)**2) ! pT ! do i=1,3 ! var(20+i) = xpint(i) ! var(23+i) = ppint(i) ! end do ! var(27) = tpint * cc ! if( betaperp .ne. 0. ) then ! Palmer A2 var(28) = (var(14)/betaperp)**2 + var(11)**2 + var(12)**2 else var(28) = 0. end if c var(29) = var(14)**2 ! r2 ! hlab = 0. if( spin .and. jpar.gt.0 .and. ABS(iptyp).eq.2 .and. & spintrk.gt.0 ) then hlab = HELICITY( pol ) else hlab = pol(3) end if var(30) = hlab ! do i=1,3 var(30+i) = bfld(i) var(33+i) = efld(i) var(36+i) = pol(i) end do ! var(40) = phaserf ! return end ! ********************************************************* subroutine MEANSD(stev,stem,xmean,xstand) ! ! Compute mean and standard deviation of stored array ! implicit double precision (a-h,o-z) dimension stem(2,6),xmean(6),xstand(6) ! do i=1,6 qx = stem(1,i) / stev xmean(i) = qx vx = stem(2,i)/stev - qx*qx if( vx .ge. 0. ) then xstand(i) = SQRT( vx ) else xstand(i) = -999. end if end do ! return end c ********************************************************* subroutine OUT_EMIT c c Print out normalized emittances c implicit double precision (a-h,o-z) include 'icommon.inc' c real*8 n0 c if( nemit .eq. 0 ) go to 900 ! write(2,*) write(2,40) 40 format(18('*'),t24,'Weighted NORMALIZED 2-D EMITTANCES', & t61,18('*')) write(2,*) write(2,70) 70 format(t3,'des',t10,'No',t17,'N',t26,'emitX',t38,'emitY', + t50,'emitZ',t61,'Tr',t69,'fET',t77,'fEL') write(2,75) 75 format(t27,'[m]',t39,'[m]',t51,'[m]') ! do 200 i=1,nemit if( i .eq. 1 ) then n0 = emevt(i) ex0 = emitx(i) ey0 = emity(i) ez0 = emitz(i) end if ! if( n0*ex0*ey0*ez0 .eq. 0. ) go to 800 ! ftr = emevcut(i) / n0 fet = 0.5 * (emitx(i)/ex0 + emity(i)/ey0) fez = emitz(i) / ez0 write(2,150)iemdes(i),emevt(i),emevcut(i),emitx(i),emity(i), & emitz(i),ftr,fet,fez 150 format(i5,2f7.0,1x,3e12.4,3f8.3) 200 continue ! 800 continue write(2,*) write(2,805) 805 format(15('*'),t20,'LBNL NORMALIZED DET(COV) EMITTANCES', & t60,19('*')) write(2,*) write(2,810) 810 format(t2,'des',t12,'emitT',t25,'emitL',t38,'emit6D', & t51,'betaT',t64,'Lcanon') write(2,815) 815 format(t12,'[m]',t25,'[m]',t38,'[m^3]',t51,'[m]',t64, & '[m GeV/c]') write(2,*) ! do 850 i=1,nemit write(2,'(i5,2x,5e13.5)') iemdes(i),emit4(i),emitl(i),emit6(i), & betper(i),canang(i) 850 continue ! 900 continue return end c ********************************************************* subroutine OUT_HIST c c Generate histograms in character format c c maximum histogram size is 50 bins c implicit double precision (a-h,o-z) include 'icommon.inc' ! Common variables: ! /dvars/ labvar c integer lowb,highb,xrange character*80 scr(25) character*10 cstr*10,fname character*6 csti c if(nhist .le. 0) go to 900 c if( hcprn.ge.20 .and. hcprn.le.99 ) then ! write out histogram contents write(fname,20) hcprn 20 format('for0',i2,'.dat') open(unit=hcprn,file=fname,status='unknown') write(hcprn,'(a)')' HISTOGRAM CONTENTS' do ih=1,nhist ip = lhpnt(ih) - 1 nhb = nhbins(ih) write(hcprn,'(a,i5)')' Histogram # ',ih write(hcprn,'(2e12.4,2x,i4)') hxmin(ih),hdx(ih),nhbins(ih) write(hcprn,'(10i7)') (hval(ip+i),i=1,nhb) end do close(hcprn) end if c for each histogram, prepare virtual page in array SCR c do 700 ih=1,nhist if( ihcon(ih) .eq. 0 ) go to 700 ip = lhpnt(ih) - 1 nhb = nhbins(ih) c clear out SCR do i=1,25 do j=1,80 scr(i)(j:j) = ' ' end do end do c c find range of occupied x bins lowb = 0 do i = 1,nhb ival = hval(ip+i) if(ival .gt. 0)then lowb = i go to 30 end if end do c 30 continue highb = 51 do i=nhb,1,-1 ival = hval(ip+i) if(ival .gt. 0)then highb = i go to 40 end if end do c 40 continue xrange = highb - lowb + 1 c c left and right vertical lines do i=2,24 scr(i)(10:10) = '|' k = nhb + 11 scr(i)(k:k) = '|' end do c c top and bottom horizontal lines do i=1,nhb k = i + 10 scr(1)(k:k) = '-' scr(25)(k:k) = '-' end do c c histogram parameters scr(1)(63:69) = 'ICOOL ' ! write(cstr,'(f10.2)')version ! scr(1)(71:80) = cstr scr(1)(71:80) = version do i=62,80 scr(2)(i:i) = '-' end do scr(3)(63:73) = 'histogram =' write(csti,'(i6)') ih scr(3)(74:79) = csti scr(4)(63:73) = 'variable =' scr(4)(77:79) = labvar( ihvar(ih) ) scr(5)(63:73) = 'z plane =' ival = ihdes(ih) write(csti,'(i6)') ival scr(5)(74:79) = csti c scr(7)(63:70) = 'lo lim =' write(cstr,'(e10.3)') hxmin(ih) scr(7)(71:80) = cstr scr(8)(63:70) = 'hi lim =' xhi = hxmin(ih) + real(nhb)*hdx(ih) write(cstr,'(e10.3)') xhi scr(8)(71:80) = cstr scr(9)(63:70) = 'step =' write(cstr,'(e10.3)') hdx(ih) scr(9)(71:80) = cstr c scr(11)(63:73) = 'contents =' ival = ihcon(ih) write(csti,'(i6)') ival scr(11)(74:79) = csti scr(12)(63:73) = 'underflow =' ival = ihund(ih) write(csti,'(i6)') ival scr(12)(74:79) = csti scr(13)(63:73) = 'overflow =' ival = ihovf(ih) write(csti,'(i6)') ival scr(13)(74:79) = csti c scr(15)(63:73) = 'x range = ' write(csti,'(i6)') xrange scr(15)(74:79) = csti c do i=62,80 scr(16)(i:i) = '-' end do do i=1,16 scr(i)(61:61) = '|' end do c c y axis values c c find maximum value mxval = 0 do i=1,nhb ival = hval(ip+i) if( ival .gt. mxval ) mxval = ival end do if(hauto .and. mxval.gt.23)then c scale bin contents to maximum value c label y axis values do i=1,23 val = real(mxval)*real(i)/23. write(cstr,'(f9.1)')val scr(25-i)(1:9) = cstr end do else c label y axis from 1 to 23 do i=1,23 write(csti,'(i6)') i scr(25-i)(3:8) = csti end do mxval = 23 end if c step = real(mxval) / 23. c c fill in the histogram data c loop over columns do 200 i=1,nhb val = hval(ip+i) k = 10 + i c loop over rows do 150 j=1,23 height = (j-1)*step if(val .gt. height)then scr(25-j)(k:k) = 'X' else go to 170 end if 150 continue c 170 continue c mark contents that are off scale when not using autoscaling if(int(val) .gt. mxval)scr(2)(k:k) = '^' 200 continue c ! Fit two plots per printed page if( ffcr .and. mod(ih,2).eq.1 ) then call FORMFEED else if( ffcr .and. mod(ih,2).eq.0 ) then write(2,*) end if ! do i=1,25 write(2,'(a80)')scr(i) end do 700 continue c 900 continue return end c ********************************************************* subroutine OUT_POL ! ! Print summary of polarization values at "emittance" planes ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 n ! write(2,*) write(2,40) 40 format(20('*'),t28,'Weighted POLARIZATION',t58,20('*')) write(2,50) 50 format(t3,'des',t12,'N',t25,'polX',t37,'polY', & t49,'polZ',t60,'helicity') ! do 200 i=1,nemit n = emevt(i) if( n .lt. 1 ) go to 200 px = polx(i) / n py = poly(i) / n pz = polz(i) / n hel = helic(i) / n write(2,150) iemdes(i),n,px,py,pz,hel 150 format(i5,f8.0,5x,4f12.4) 200 continue ! return end ! ********************************************************* subroutine OUT_RHIST ! ! Generate R-history plots in character format ! c plots have 70 horizontal bins c plots have 23 vertical bins c implicit double precision (a-h,o-z) include 'icommon.inc' c integer*2 r,c character*80 scr(25) character*10 cstr character*6 csti c if(nrhist .le. 0) go to 900 c if( rhprin ) then ! print out R-history contents write(2,'(a)')' R-HISTORY CONTENTS' do ih=1,nrhist write(2,'(a,i5)')' R-History # ',ih do iz=1,70 izp = irhzmin(ih) + iz*irhdz(ih) - 1 write(2,'(i4,4e12.4)') izp,(rhyval(iz,j,ih),j=1,4) end do end do end if ! c for each R-history, prepare virtual page in array SCR c do 700 ih=1,nrhist ! c clear out SCR do i=1,25 do j=1,80 scr(i)(j:j) = ' ' end do end do c c Left vertical line do i=1,24 scr(i)(10:10) = '|' end do scr(12)(10:10) = '+' scr(12)(3:5) = labvar( irhyvar(ih) ) c ! Horizontal line do i=10,80 scr(24)(i:i) = '-' end do scr(24)(27:27) = '+' scr(24)(45:45) = '+' scr(24)(62:62) = '+' scr(25)(42:48) = 'Plane #' ! c Plot parameters scr(5)(3:7) = 'ICOOL' ! write(cstr,'(f10.2)')version ! scr(6)(1:7) = cstr(4:10) scr(6)(1:7) = version scr(8)(1:9) = 'R history' write(csti,'(i6)') ih scr(9)(1:5) = csti(2:6) do i=1,9 scr(4)(i:i) = '-' scr(10)(i:i) = '-' end do ! c if rhauto = true => force all points to fill plot c if(rhauto)then ymin = 1e9 ymax = -1e9 c Find maximum and minimum values ! do ip=1,4 do 80 iz=1,70 y = rhyval(iz,ip,ih) if( y .eq. 0. ) go to 75 if( y .gt. ymax) ymax = y if( y .lt. ymin) ymin = y 75 continue 80 end do end do ! else ! use input file limits ymin = rhymin(ih) ymax = rhymax(ih) end if ! dy = (ymax-ymin) / 23. ! ! Label minimum and maximum values on plot izp = irhzmin(ih) write(csti,'(i6)') izp scr(25)(5:10) = csti write(csti,'(i6)') izp+69*irhdz(ih) scr(25)(75:80) = csti write(cstr,'(e10.3)') ymin scr(24)(1:10) = cstr write(cstr,'(e10.3)') ymax scr(1)(1:10) = cstr ! c Fill in the R-history data ! if( dy .gt. 0. ) then do ip=1,4 ! do 180 iz=1,70 c = 10 + iz y = rhyval(iz,ip,ih) if( y .eq. 0. ) go to 175 r = 24. - (y-ymin) / dy if( r .lt. 1 ) go to 175 if( r .gt. 24 ) go to 175 ! if( ip .eq. 1 ) then scr(r)(c:c) = '-' else if( ip .eq. 2 ) then scr(r)(c:c) = '+' else if( ip .eq. 3 ) then scr(r)(c:c) = 'm' else scr(r)(c:c) = 's' end if ! 175 continue 180 end do end do ! end if ! ! Fit two plots per printed page if( ffcr .and. mod(ih,2).eq.1 ) then call FORMFEED else if( ffcr .and. mod(ih,2).eq.0 ) then write(2,*) end if ! ! Write out the virtual page ! do i=1,25 write(2,'(a80)')scr(i) end do ! 700 continue ! end of loop on R-histories c 900 continue return end c ********************************************************* subroutine OUT_SCAT c c Generate scatterplots in character format c c plots can have a maximum of 50 horizontal bins c plots can have a maximum of 23 vertical bins c implicit double precision (a-h,o-z) include 'icommon.inc' ! Common variables: ! /dvars/ labvar c integer*2 r,c character*80 scr(25) character*10 cstr character*6 csti c if(nscat .le. 0) go to 900 c c for each scatterplot, prepare virtual page in array SCR c do 700 ih=1,nscat if( iscon(ih) .eq. 0 ) go to 700 ip = lspnt(ih) - 1 nxb = nsxbin(ih) nyb = nsybin(ih) c clear out SCR do i=1,25 do j=1,80 scr(i)(j:j) = ' ' end do end do c c left and right vertical lines do i=1,nyb scr(25-i)(10:10) = '|' k = nxb + 11 scr(25-i)(k:k) = '|' end do c c top and bottom horizontal lines do i=1,nxb k = i + 10 scr(25)(k:k) = '-' scr(24-nyb)(k:k) = '-' end do c c scatterplot parameters scr(1)(63:69) = 'ICOOL ' ! write(cstr,'(f10.2)')version ! scr(1)(71:80) = cstr scr(1)(71:80) = version do i=62,80 scr(2)(i:i) = '-' end do scr(3)(63:73) = 'plot =' write(csti,'(i6)') ih scr(3)(74:79) = csti scr(4)(63:73) = 'var(x) =' scr(4)(77:79) = labvar( isxvar(ih) ) scr(5)(63:73) = 'dest(x) =' ival = isxdes(ih) write(csti,'(i6)') ival scr(5)(74:79) = csti scr(6)(63:73) = 'var(y) =' scr(6)(77:79) = labvar( isyvar(ih) ) scr(7)(63:73) = 'dest(y) =' ival = isydes(ih) write(csti,'(i6)') ival scr(7)(74:79) = csti c scr(9)(63:70) = 'lo x =' write(cstr,'(e10.3)') sxmin(ih) scr(9)(71:80) = cstr scr(10)(63:70) = 'hi x =' xhi = sxmin(ih) + real( nxb*sdx(ih) ) write(cstr,'(e10.3)') xhi scr(10)(71:80) = cstr scr(11)(63:70) = 'step x =' write(cstr,'(e10.3)') sdx(ih) scr(11)(71:80) = cstr scr(12)(63:70) = 'lo y =' write(cstr,'(e10.3)') symin(ih) scr(12)(71:80) = cstr scr(13)(63:70) = 'hi y =' yhi = symin(ih) + real( nyb*sdy(ih) ) write(cstr,'(e10.3)') yhi scr(13)(71:80) = cstr scr(14)(63:70) = 'step y =' write(cstr,'(e10.3)') sdy(ih) scr(14)(71:80) = cstr c scr(16)(63:73) = 'contents =' ival = iscon(ih) write(csti,'(i6)') ival scr(16)(74:79) = csti scr(17)(63:73) = 'under x =' ival = isxund(ih) write(csti,'(i6)') ival scr(17)(74:79) = csti scr(18)(63:73) = 'over x =' ival = isxovf(ih) write(csti,'(i6)') ival scr(18)(74:79) = csti scr(19)(63:73) = 'under y =' ival = isyund(ih) write(csti,'(i6)') ival scr(19)(74:79) = csti scr(20)(63:73) = 'over y =' ival = isyovf(ih) write(csti,'(i6)') ival scr(20)(74:79) = csti c do i=62,80 scr(21)(i:i) = '-' end do do i=1,21 scr(i)(61:61) = '|' end do c ! if sauto = true and maxval > 61, ! then scale all contents to a maximum of 61 c if(sauto)then mxval = 0 k = 0 c find maximum value do i=1,nxb do j=1,nyb k = k + 1 ival = sval(ip+k) if( ival .gt. mxval ) mxval = ival end do end do if(mxval .gt. 61)then k = 0 do i=1,nxb do j=1,nyb k = k + 1 sval(ip+k) = sval(ip+k)*61/mxval end do end do end if end if c c fill in the scatterplot data k = 0 do 200 j=1,nyb r = 25 - j do 150 i=1,nxb c = 10 + i k = k + 1 ival = sval(ip+k) ! if(ival .eq. 0)then ! use bin content to determine plot symbol scr(r)(c:c) = ' ' else if(ival .eq. 1)then scr(r)(c:c) = '+' else if(ival.gt.1 .and. ival.lt.10)then scr(r)(c:c) = CHAR(ival + 48) else if(ival.gt.9 .and. ival.lt.36)then scr(r)(c:c) = CHAR(ival + 87) else if(ival.gt.35 .and. ival.lt.62)then scr(r)(c:c) = CHAR(ival + 29) else if(ival .gt. 61)then scr(r)(c:c) = '*' end if c 150 continue 200 continue c ! Fit two plots per printed page if( ffcr .and. MOD(ih,2).eq.1 ) then call FORMFEED else if( ffcr .and. MOD(ih,2).eq.0 ) then write(2,*) end if ! do i=1,25 ! write out the array write(2,'(a80)')scr(i) end do 700 continue c 900 continue return end c ********************************************************* subroutine OUT_STATS c c Compute and list statistics for histogram and scatter variables c implicit double precision (a-h,o-z) include 'icommon.inc' c real*8 kurtx,kurty real*8 m1x(20),m2x(20),m3x(20),m4x(20) real*8 m1y(20),m2y(20),m3y(20),m4y(20) real*8 n,sd_x(20),sd_y(20) ! ! --------------------------------------------------------- if( nhist .eq. 0 ) go to 232 ! write(2,*) write(2,40) 40 format(20('*'),t28,'HISTOGRAM STATISTICS',t58,20('*')) write(2,100) 100 format(' ih',t7,'var',t13,'des',t22,'N',t31,'mean',t43,'sd', + t55,'lo',t69,'hi') c do 200 i=1,nhist n = st(1,i) if(n .ge. 1.)then m1x(i) = st(2,i) / n m2x(i) = st(3,i) / n m3x(i) = st(4,i) / n m4x(i) = st(5,i) / n vr = m2x(i) - m1x(i)*m1x(i) if(vr .ge. 0.)then sd_x(i) = sqrt(vr) else sd_x(i) = -998. end if else m1x(i) = -999. sd_x(i) = -999. end if write(2,150)i,ihvar(i),ihdes(i),n,m1x(i),sd_x(i),st(6,i),st(7,i) 150 format(2i4,i6,f11.1,4e12.4) 200 continue ! --------------------------------------------------------- write(2,*) write(2,205) 205 format(20('*'),t30,'HISTOGRAM MOMENTS',t58,20('*')) write(2,210) 210 format(' ih',t7,'var',t13,'des',t19,'Skewness',t31,'Kurtosis') c do 230 i=1,nhist n = st(1,i) if( n .eq. 0 ) go to 230 if( sd_x(i) .le. 0. ) go to 230 skewx = m3x(i) - 3.*m1x(i)*m2x(i) + 2.*m1x(i)**3 skewx = skewx / sd_x(i)**3 kurtx = m4x(i) - 4.*m1x(i)*m3x(i) + 6.*m1x(i)**2*m2x(i) & -3.*m1x(i)**4 kurtx = kurtx / sd_x(i)**4 - 3. write(2,220)i,ihvar(i),ihdes(i),skewx,kurtx 220 format(2i4,i6,2f12.4) 230 continue ! --------------------------------------------------------- 232 continue if( nscat .eq. 0 ) go to 900 ! write(2,*) write(2,235) 235 format(20('*'),t28,'SCATTERPLOT STATISTICS',t58,20('*')) write(2,240) 240 format(' is',t6,'var',t12,'des',t16,'ax',t27,'N',t35,'mean' & ,t47,'sd',t59,'lo',t71,'hi') c do 300 i=1,nscat n = stsp(1,i) if(n .ge. 1.)then m1x(i) = stsp(2,i) / n m2x(i) = stsp(3,i) / n m3x(i) = stsp(4,i) / n m4x(i) = stsp(5,i) / n vr = m2x(i) - m1x(i)*m1x(i) if(vr .ge. 0.)then sd_x(i) = sqrt(vr) else sd_x(i) = -998. end if else m1x(i) = -999. sd_x(i) = -999. end if write(2,250)i,isxvar(i),isxdes(i),n,m1x(i),sd_x(i), + stsp(6,i),stsp(7,i) 250 format(2i4,i6,2x,'x',2x,f12.1,4e12.4) ! ........................................................ if(n .ge. 1.)then m1y(i) = stsp(8,i) / n m2y(i) = stsp(9,i) / n m3y(i) = stsp(10,i) / n m4y(i) = stsp(11,i) / n vr = m2y(i) - m1y(i)*m1y(i) if(vr .ge. 0.)then sd_y(i) = sqrt(vr) else sd_y(i) = -998. end if else m1y(i) = -999. sd_y(i) = -999. end if write(2,270)i,isyvar(i),isydes(i),n,m1y(i),sd_y(i), + stsp(12,i),stsp(13,i) 270 format(2i4,i6,2x,'y',2x,f12.1,4e12.4) 300 continue ! --------------------------------------------------------- ! write(2,*) write(2,305) 305 format(20('*'),t29,'SCATTERPLOT MOMENTS',t58,20('*')) write(2,310) 310 format(' is',t7,'vx',t13,'dx',t17,'vy',t23,'dy', & t33,'r',t42,'skewx',t52,'kurtx',t62,'skewy',t72,'kurty') c do 330 i=1,nscat n = stsp(1,i) if( n .eq. 0 ) go to 330 if( sd_x(i).le.0. .or. sd_y(i).le.0. ) go to 330 ! skewx = m3x(i) - 3.*m1x(i)*m2x(i) + 2.*m1x(i)**3 skewx = skewx / sd_x(i)**3 kurtx = m4x(i) - 4.*m1x(i)*m3x(i) + 6.*m1x(i)**2*m2x(i) & -3.*m1x(i)**4 kurtx = kurtx / sd_x(i)**4 - 3. ! skewy = m3y(i) - 3.*m1y(i)*m2y(i) + 2.*m1y(i)**3 skewy = skewy / sd_y(i)**3 kurty = m4y(i) - 4.*m1y(i)*m3y(i) + 6.*m1y(i)**2*m2y(i) & -3.*m1y(i)**4 kurty = kurty / sd_y(i)**4 - 3. ! fac = sd_x(i) * sd_y(i) if( fac .le. 0. ) go to 330 r = (stsp(14,i) / n - m1x(i) * m1y(i)) / fac ! write(2,320)i,isxvar(i),isxdes(i),isyvar(i),isydes(i), & r,skewx,kurtx,skewy,kurty 320 format(2i4,i6,i4,i6,f12.3,4f10.4) 330 continue write(2,*) ! --------------------------------------------------------- 900 continue return end c ********************************************************* subroutine OUT_ZHIST c ! Generate Z-history plots in character format c c plots have 70 horizontal bins c plots have 23 vertical bins c implicit double precision (a-h,o-z) include 'icommon.inc' ! Common variables: ! /dvars/ labvar c integer*2 r,c character*80 scr(25) character*10 cstr character*6 csti c if(nzhist .le. 0) go to 900 c if( zhprin ) then ! print out z-history contents write(2,'(a)')' Z-HISTORY CONTENTS' do ih=1,nzhist write(2,'(a,i5)')' Z-History # ',ih do ip=0,nzhpar(ih) write(2,'(a,i5)')' Particle # ',ip do iz=1,izhknt(ip,ih) write(2,'(i4,2e12.4)') iz,zhxval(iz,ip,ih),zhyval(iz,ip,ih) end do end do end do end if ! c for each z-history, prepare virtual page in array SCR c do 700 ih=1,nzhist ! c clear out SCR do i=1,25 do j=1,80 scr(i)(j:j) = ' ' end do end do c c Left vertical line do i=1,24 scr(i)(10:10) = '|' end do scr(12)(10:10) = '+' scr(12)(3:5) = labvar( izhyvar(ih) ) c ! Horizontal line do i=10,80 scr(24)(i:i) = '-' end do scr(24)(27:27) = '+' scr(24)(45:45) = '+' scr(24)(62:62) = '+' scr(25)(45:45) = 'Z' ! c Scatterplot parameters scr(5)(3:7) = 'ICOOL' ! write(cstr,'(f10.2)')version ! scr(6)(1:7) = cstr(4:10) scr(6)(1:7) = version scr(8)(1:9) = 'Z history' write(csti,'(i6)') ih scr(9)(1:5) = csti(2:6) do i=1,9 scr(4)(i:i) = '-' scr(10)(i:i) = '-' end do ! xmin = zhxmin(ih) xmax = zhxmax(ih) ! c if zhauto = true => force all points to fill plot c if(zhauto)then dx = (xmax - xmin) / 70. xmaxp = xmax ! ymin = 1e9 ymax = -1e9 c Find maximum and minimum values ! do ip=0,nzhpar(ih) do 100 iz=1,izhknt(ip,ih) if( iz .gt. 70 ) go to 100 y = zhyval(iz,ip,ih) if( y .gt. ymax) ymax = y if( y .lt. ymin) ymin = y 100 continue end do ! else ! use input file limits dx = zhdx(ih) xmaxp = xmin + 70 * dx ymin = zhymin(ih) ymax = zhymax(ih) end if ! dy = (ymax-ymin) / 23. ! ! Label minimum and maximum values on plot write(cstr,'(e10.3)') xmin scr(25)(10:19) = cstr write(cstr,'(e10.3)') xmaxp scr(25)(71:80) = cstr write(cstr,'(e10.3)') ymin scr(24)(1:10) = cstr write(cstr,'(e10.3)') ymax scr(1)(1:10) = cstr ! c Fill in the z-history data ! if( dx .gt. 0. .and. dy .gt. 0. ) then do 250 ip=0,nzhpar(ih) ! do 200 iz=1,izhknt(ip,ih) if( iz .gt. 70 ) go to 200 x = zhxval(iz,ip,ih) c = (x-xmin) / dx + 10 if( c .lt. 11 ) c = 11 if( c .gt. 80 ) c = 80 y = zhyval(iz,ip,ih) r = 24. - (y-ymin) / dy if( r .lt. 1 ) go to 200 if( r .gt. 24 ) go to 200 ! if( ip .lt. 10 ) then scr(r)(c:c) = CHAR(ip + 48) else scr(r)(c:c) = CHAR(ip + 55) end if ! 200 continue ! 250 continue c end if ! ! Fit two plots per printed page if( ffcr .and. mod(ih,2).eq.1 ) then call FORMFEED else if( ffcr .and. mod(ih,2).eq.0 ) then write(2,*) end if ! ! Write out the virtual page ! do i=1,25 write(2,'(a80)')scr(i) end do ! 700 continue ! end of loop on z-histories c 900 continue return end ! ********************************************************* subroutine PZ_AMP_COR(betax,betay,slp) ! ! Determine Pz - amplitude correlation ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 v(5) ! do i=1,5 v(i) = 0. end do ! do 200 ip=1,npart ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 200 ! if( pxycorr ) then call CANONICAL(pxc,pyc) else pxc = pp(1) pyc = pp(2) end if ! if( ipzcor .eq. 1 ) then xpr = pxc / pp(3) ypr = pyc / pp(3) amp = SQRT( xpr**2 + ypr**2 + (xp(1)/betax)**2 & + (xp(2)/betay)**2 ) ! else if( ipzcor .eq. 2 ) then r2 = xp(1)**2 + xp(2)**2 pt2 = pxc**2 + pyc**2 pm2 = pmass**2 amp = pt2/pm2 + ((pcharge*bfld(3))**2*r2)/(4.*pm2) amp = SQRT(amp) ! else amp = 0. end if ! v(1) = v(1) + evtwt * pp(3) v(2) = v(2) + evtwt * amp v(3) = v(3) + evtwt * amp*amp v(4) = v(4) + evtwt * amp*pp(3) v(5) = v(5) + evtwt 200 continue ! ! Get parameters of least square fit to straight line (cf. NR p. 656) s = v(5) sx = v(2) sy = v(1) sxx = v(3) sxy = v(4) del = s*sxx - sx*sx if( del .ne. 0. ) then ! off = (sxx*sy - sx*sxy) / del slp = (s*sxy - sx*sy) / del end if ! return end ! ********************************************************* subroutine RHISTORY ! ! Compute and load R-history data into array ! implicit double precision (a-h,o-z) include 'icommon.inc' ! do 600 ih=1,nrhist if( j7rec.lt.irhzmin(ih) .or. j7rec.gt.irhzmax(ih) ) go to 600 izp = (j7rec - irhzmin(ih))/irhdz(ih) + 1 ! stev = 0. stem = 0. stem2 = 0. stemlo = 1e20 stemhi = -1e20 ! ! Accumulate statistics; get min and max values do 100 ip=1,npart call READ_EVT(ip) if( ipflg .ne. 0 ) go to 100 call LOADVAR(j7rec) v = var( irhyvar(ih) ) stev = stev + 1. stem = stem + v stem2 = stem2 + v*v if( v .lt. stemlo ) stemlo = v if( v .gt. stemhi ) stemhi = v 100 continue ! ! Calculate mean, standard deviation; load array if( stev .gt. 0. ) then qx = stem / stev vx = stem2/stev - qx*qx else qx = 0. vx = 0. end if if( vx .gt. 0. ) then sd = SQRT( vx ) else sd = 0. end if ! rhyval(izp,1,ih) = stemlo rhyval(izp,2,ih) = stemhi rhyval(izp,3,ih) = qx rhyval(izp,4,ih) = sd ! 600 continue ! return end ! ********************************************************* subroutine SIGMA_CUT(xmean,xstand,ncut) ! ! Flag particle with {x...Pz} in tails of distributions ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 xmean(6),xstand(6),v(6) integer ncut(6) ! do i=1,6 ncut(i) = 0 end do ! do 200 ip=1,npart ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 200 ! ! Correct Px and Py for vector potential if( pxycorr ) then call CANONICAL(pxc,pyc) else pxc = pp(1) pyc = pp(2) end if ! energy = SQRT( pxc**2 + pyc**2 + pp(3)**2 + pmass**2 ) z = pp(3) / energy * cc * tp ! v(1) = xp(1) v(2) = xp(2) v(3) = z v(4) = pxc v(5) = pyc v(6) = pp(3) ! do i=1,6 arg = v(i) - xmean(i) ! difference from mean if( ABS(arg).gt.sig_cut*xstand(i) & .and. xstand(i).gt.0. ) then ncut(i) = ncut(i) + 1 ipflg = 10 call WRITE_EVT(ip) go to 200 end if end do ! 200 continue ! return end c ********************************************************* subroutine STATH(v,ih) c c Accumulate statistics on the histogram variables c implicit double precision (a-h,o-z) include 'icommon.inc' c st(1,ih) = st(1,ih) + 1. st(2,ih) = st(2,ih) + v st(3,ih) = st(3,ih) + v*v st(4,ih) = st(4,ih) + v*v*v st(5,ih) = st(5,ih) + v*v*v*v if(v .lt. st(6,ih)) st(6,ih) = v if(v .gt. st(7,ih)) st(7,ih) = v c return end c ********************************************************* subroutine STATSP(vx,vy,is) c c Accumulate statistics on the scatterplot variables c implicit double precision (a-h,o-z) include 'icommon.inc' c stsp(1,is) = stsp(1,is) + 1. c stsp(2,is) = stsp(2,is) + vx stsp(3,is) = stsp(3,is) + vx*vx stsp(4,is) = stsp(4,is) + vx*vx*vx stsp(5,is) = stsp(5,is) + vx*vx*vx*vx if(vx .lt. stsp(6,is)) stsp(6,is) = vx if(vx .gt. stsp(7,is)) stsp(7,is) = vx c stsp(8,is) = stsp(8,is) + vy stsp(9,is) = stsp(9,is) + vy*vy stsp(10,is) = stsp(10,is) + vy*vy*vy stsp(11,is) = stsp(11,is) + vy*vy*vy*vy if(vy .lt. stsp(12,is)) stsp(12,is) = vy if(vy .gt. stsp(13,is)) stsp(13,is) = vy c stsp(14,is) = stsp(14,is) + vx * vy c return end