program BUNCH9 ! ! Find number of bunched muons inside cooling momentum acceptance ! parameter ( version = 1.11 ) real*4 xp(3),pp(3),bfld(3),efld(3),pol(3) real*4 sarr(2,100,5) ! statistics array real nnmu,ngood,nplo,nphi,ntlo,nthi,nflg,nbox,nbeam real ntotal,ngoodmax,nb character(80) header character(20) fname ! do 20 i=1,2 ! initialize statistics array do 20 j=1,100 do k=1,3 sarr(i,j,k) = 0. end do sarr(i,j,4) = 1e9 sarr(i,j,5) = -1e9 20 continue ! open(unit=8,file='bunch9.inp',status='old') read(8,'(a20)') fname ! for009 input file read(8,*) freq ! [MHz] read(8,*) plow,phigh ! momentum band [GeV/c] read(8,*) nb1,nb2 ! bunch # range for new 009 & 003 files read(8,*) toff ! start time offset [s] read(8,*) boff ! bunchlength offset from period/2 [s] ! ! File 9 is the ICOOL postprocessor output open(unit=9,file=fname,status='old',iostat=ioc) open(unit=10,file='bunch9.dat',status='unknown',iostat=ioc) open(unit=13,file='bunch9.f03',status='unknown') open(unit=19,file='bunch9.f09',status='unknown') ! per = 1. / (freq*1e6) ! write(*,'(a,f5.2)') 'Welcome to program BUNCH9, version ', & version write(10,'(a,f5.2)') 'Welcome to program BUNCH9, version ', & version write(10,'(a20)') fname write(10,*) freq write(10,*) plow,phigh write(10,*) nb1,nb2 write(10,*) toff write(10,*) boff write(*,50) write(10,50) 50 format(t4,'#',t11,'t0',t25,'good',t37,'flg',t43,'nomu',t51, & 'plo',t57,'phi',t64,'tlo',t71,'thi') itmax = 0 ngoodmax = 0 tlow = 1e9 tref = 0. ! --------------------------------------------------------- do 500 it=1,18 ! Set beginning of adjustable bunching window inside rf wavelength tbeg = (it-1)*per/18. ! nnmu = 0 ngood = 0 nplo = 0 nphi = 0 ntlo = 0 nthi = 0 nflg = 0 ! read(9,*) ! skip header records read(9,*) read(9,*) nrec = 0 ! ......................................................... 100 continue nrec = nrec + 1 ! ......................................................... read(9,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol if( ioc .gt. 0 ) then write(*,'(a,i10)') 'Error reading (9) record ',nrec go to 900 else if( ioc .eq. -1 ) then ! EOF go to 300 end if ! ! --------------------------------------------------------- ! ! Make selections ! if( ievt .le. 0 ) then ! reference particle tref = tp go to 100 ! skip reference particle end if ! if( ipflg .ne. 0 ) then nflg = nflg + evtwt go to 100 end if ! if( ABS(iptyp) .ne. 2 ) then nnmu = nnmu + evtwt go to 100 ! muons only end if ! ......................................................... ! ! Local definitions px = pp(1) py = pp(2) pz = pp(3) pm = SQRT( px*px + py*py + pz*pz ) if( pm .lt. plow ) then nplo = nplo + evtwt go to 100 end if ! if( pm .gt. phigh ) then nphi = nphi + evtwt go to 100 end if ! ! Transform all particle times into a single period tp = tp - tref + toff count = tp / per if (count .lt. 0.) count = count - 1. u = tp - INT(count) * per ! ! The following cut throws out particles that are not in the selected bucket if( tbeg .le. per/2. ) then if( u .lt. tbeg ) then ntlo = ntlo + evtwt go to 100 end if ! if( u .gt. tbeg+per/2.+boff ) then nthi = nthi + evtwt go to 100 end if ! else ! tbeg > per/2 if( u.lt.tbeg .and. u.gt.tbeg-per/2.+boff) then ntlo = ntlo + evtwt go to 100 end if ! end if ! ! NGOOD counts all "good" particles in the full bunch train ngood = ngood + evtwt ! This tlow has the momentum cut applied if( tp .lt. tlow ) tlow = tp ! go to 100 ! get next record ! ......................................................... 300 continue write(*,490) it,tbeg,ngood,nflg,nnmu,nplo,nphi,ntlo,nthi write(10,490) it,tbeg,ngood,nflg,nnmu,nplo,nphi,ntlo,nthi 490 format(i5,2x,e12.6,1x,f8.0,4x,6f7.0) ntotal = ngood+nflg+nnmu+nplo+nphi+ntlo+nthi if( ngood .gt. ngoodmax ) then ngoodmax = ngood itmax = it end if rewind(9) 500 continue ! end of loop on IT ! frac = ngoodmax / ntotal write(*,505) 'NGOOD,NTOTAL,FRAC: ',ngoodmax,ntotal,frac 505 format(a,2f10.0,f12.4) write(10,505) 'NGOOD,NTOTAL,FRAC: ',ngoodmax,ntotal,frac ! --------------------------------------------------------- ! ! Extract information about optimally bunched beam ! read(9,'(a80)') header write(19,'(a80)') header read(9,'(a80)') header write(19,'(a80)') header read(9,'(a80)') header write(19,'(a80)') header nbox = 0 nrec = 0 tbeg = (itmax-1)* per/18. ! ......................................................... 510 continue nrec = nrec + 1 read(9,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol if( ioc .gt. 0 ) then write(*,'(a,i10)') 'Error reading (9) record ',nrec go to 900 else if( ioc .eq. -1 ) then ! EOF go to 600 end if ! ! Make selections ! if( ievt .le. 0 ) then ! write refpar to file, but don't use in analysis write(19,530)ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol go to 510 end if ! if( ipflg .ne. 0 ) go to 510 if( ABS(iptyp) .ne. 2 ) go to 510 ! px = pp(1) py = pp(2) pz = pp(3) sx = pol(1) sy = pol(2) sz = pol(3) pm = SQRT( px*px + py*py + pz*pz ) hel = (sx*px + sy*py + sz*pz) / pm ! if( pm .lt. plow ) go to 510 if( pm .gt. phigh ) go to 510 ! ! The following cut throws out particles that are not in the optimal bucket tp = tp - tref + toff count = tp / per if (count .lt. 0.) count = count - 1. u = tp - INT(count) * per if( tbeg .le. per/2. ) then if( u .lt. tbeg ) go to 510 if( u .gt. tbeg+per/2.+boff ) go to 510 else ! tbeg > per/2 if( u.lt.tbeg .and. u.gt.tbeg-per/2.+boff) go to 510 end if ! jb = (tp-tlow)/per + 1 ! bunch number if( jb.gt.0 .and. jb.le.100 ) then call STATS(1,jb,pm,evtwt,sarr) ! skip muons at production (no spin value) if( ABS(ipnum) .ne. 0 ) then call STATS(2,jb,hel,evtwt,sarr) end if ! if( jb.ge.nb1 .and. jb.le.nb2 ) then nbox = nbox + evtwt write(19,530)ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol 530 format(i7,2i4,i5,i6,1p,e18.10,6d23.15,7e13.5,e16.8,3e13.5) end if end if ! go to 510 ! get next record ! --------------------------------------------------------- ! ! Get statistics on the momentum and helicity distribution in each bunch ! 600 continue write(10,602) 602 format(t3,'IB',t9,'numb',t21,'Pmean',t34,'dP',t48,'Plow', & t62,'Phigh',t76,'Hmn',t90,'dH') ! do ib=1,100 call OUT_STAT(1,ib,sarr) nb = sarr(1,ib,1) pmmn = sarr(1,ib,2) pmsd = sarr(1,ib,3) pmlo = sarr(1,ib,4) pmhi = sarr(1,ib,5) call OUT_STAT(2,ib,sarr) hmn = sarr(2,ib,2) hsd = sarr(2,ib,3) write(10,'(i4,f8.0,6f14.6)') ib,nb,pmmn,pmsd,pmlo,pmhi,hmn,hsd end do ! --------------------------------------------------------- ! ! Write out beam file (for003) for truncated bunch train ! rewind(9) read(9,'(a80)') header write(13,'(a80)') header read(9,'(a80)') header read(9,'(a80)') header nbeam = 0 nrec = 0 ! ......................................................... 610 continue nrec = nrec + 1 read(9,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol if( ioc .gt. 0 ) then write(*,'(a,i10)') 'Error reading (9) record ',nrec go to 900 else if( ioc .eq. -1 ) then ! EOF go to 700 end if ! ! Make selections ! if( ievt .eq. 0 ) then write(13,'(7e14.6)') xp(3),pp(3),tp, 0.,0.,0., 2. go to 610 end if if( ievt .eq. -1 ) go to 610 if( ipflg .ne. 0 ) go to 610 if( ABS(iptyp) .lt. 2 ) go to 610 ! kill electrons ! ! n.b. tp can be < tlow because the momentum cut is not applied here ! This gives JB < 1 jb = (tp-tref+toff+tbeg-tlow)/per + 1 ! bunch number if( jb.ge.nb1 .and. jb.le.nb2 ) then nbeam = nbeam + evtwt write(13,640) ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp,pol 640 format(2i8,2i6,2e14.6,9e14.6) end if ! go to 610 ! get next record ! --------------------------------------------------------- 700 continue write(*,705) 'NBOX,NBEAM: ',nbox,nbeam 705 format(a,2f10.1) write(10,705) 'NBOX,NBEAM: ',nbox,nbeam ! 900 continue end ! ********************************************************* subroutine STATS(ih,jb,v,w,sarr) ! ! Accumulate statistics ! real*4 sarr(2,100,5) ! statistics array ! sarr(ih,jb,1) = sarr(ih,jb,1) + w sarr(ih,jb,2) = sarr(ih,jb,2) + w*v sarr(ih,jb,3) = sarr(ih,jb,3) + w*v*v if( v .lt. sarr(ih,jb,4)) sarr(ih,jb,4) = v if( v .gt. sarr(ih,jb,5)) sarr(ih,jb,5) = v ! end ! ********************************************************* subroutine OUT_STAT(ih,jb,sarr) ! ! Calculate mean and error on mean ! real*4 sarr(2,100,5) ! statistics array ! wtot = sarr(ih,jb,1) swv = sarr(ih,jb,2) swvv = sarr(ih,jb,3) ! if( wtot .gt. 0. ) then sarr(ih,jb,2) = swv / wtot term = swvv / wtot - (swv/wtot)**2 if( term .gt. 0. ) then sd = SQRT(term) sarr(ih,jb,3) = sd / SQRT(wtot) else sarr(ih,jb,3) = -998. end if else sarr(ih,jb,2) = -999. sarr(ih,jb,3) = -999. end if ! end