program ENDOF9 ! ! Removes failed tracks from ICOOL file FOR009.DAT ! Writes file containing all events that make it to the end ! (or writes tracks specified by ECALC9F or EIGEMIT) ! parameter ( mxgood = 100000 ) real*8 tp,xp(3),pp(3),bfld(3),evtwt,efld(3),sarc,pol(3) integer igde(mxgood),igdp(mxgood),f9dp character*80 line,fmt data igde/mxgood*0/,igdp/mxgood*0/ data idec/0/,iexp/0/ ! version = 1.12 write(*,'(a,f6.2)') 'Welcome to program ENDOF9, version ',version ! open(unit=9,file='for009.dat',status='old',iostat=ioc) if(ioc .ne. 0) then write(*,*)'Couldnt open FOR009.DAT' go to 900 end if open(unit=10,file='endof9.dat',status='unknown',recl=512) ! ! --------------------------------------------------------- ! ECALC9F cut option open(unit=2,file='ecalc9f.cut',status='old',iostat=ioc) if( ioc .eq. 0 ) then ! file exists, use it write(*,*) 'Selecting particles using ECALC9F.CUT' kk = 0 ! Read whole file to find last region number 20 continue read(2,*,iostat=ios) ive,ivp,ivr if( ios .eq. -1 ) then ! EOF backspace(unit=2) read(2,*,iostat=ios) ive,ivp,ivr lreg = ivr rewind(unit=2) go to 50 else go to 20 end if ! Go thru file again, select events in last region 50 continue read(2,*,iostat=ios) ive,ivp,ivr if( ios .eq. 0 ) then ! good read if( ivr .ne. lreg ) go to 50 kk = kk + 1 if( kk .gt. mxgood ) then write(*,*) 'Too many good tracks' go to 900 end if igde(kk) = ive igdp(kk) = ivp go to 50 else go to 230 end if end if ! ! --------------------------------------------------------- ! EIGEMIT cut option open(unit=2,file='eigemit.acc',status='old',iostat=ioc) if( ioc .eq. 0 ) then ! file exists, use it write(*,*) 'Selecting particles using EIGEMIT.ACC' 70 continue read(2,*,iostat=ios) ive,ivp,ivr if( ios .eq. 0 ) then ! good read kk = kk + 1 if( kk .gt. mxgood ) then write(*,*) 'Too many good tracks' go to 900 end if igde(kk) = ive igdp(kk) = ivp go to 70 else ! assume we have EOF go to 230 end if end if ! ! --------------------------------------------------------- ! Normal option: no external acceptance file, just work with for009.dat knt = 0 ! count lines in final region ksrg = -1 ! flag for region change ! read(9,*) ! skip header records read(9,*) read(9,*) ! 100 continue read(9,*,iostat=ioc) ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol ! if( ioc .gt. 0 ) then write(*,*)'Error reading FOR009.DAT' go to 900 else if( ioc .eq. -1 ) then ! EOF ! Position for009 at first particle in the last region do i=1,knt+1 backspace(unit=9) end do go to 200 end if ! if( jsrg .ne. ksrg ) then ! new region knt = 1 ksrg = jsrg else ! same region knt = knt + 1 end if ! go to 100 ! --------------------------------------------------------- 200 continue ! Make list of good tracks in the last region on the file kk = 0 jevt = -999 ! 220 continue read(9,*,iostat=ioc) ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol ! if( ioc .gt. 0 ) then write(*,*)'Error reading FOR009.DAT' go to 900 else if( ioc .eq. -1 ) then ! EOF go to 230 end if ! if( ipflg .ne. 0 ) go to 220 if( ievt .eq. jevt ) go to 220 ! another entry for same track ! kk = kk + 1 if( kk .gt. mxgood ) then write(*,*) 'Too many good tracks' go to 900 end if igde(kk) = ievt igdp(kk) = ipnum jevt = ievt go to 220 ! ! --------------------------------------------------------- ! All options converge here 230 continue write(*,235) kk 235 format(i6,' good tracks at end of FOR009.DAT') ! Write new file with bad tracks eliminated ! rewind(9) read(9,'(a80)') line ! title write(10,'(a80)') line read(9,'(a80)') line ! units write(10,'(a80)') line ! read(9,*) ! skip line of column labels read(9,'(a80)') line ! read beginning of 1st data row as characters do i=1,80 if( line(i:i) .eq. '.' ) idec = i if( line(i:i) .eq. 'E' ) iexp = i if( line(i:i).eq.' ' .and. iexp.gt.0 ) then if( i-iexp .lt. 5 ) then f9dp = iexp - idec - 1 else f9dp = 17 end if go to 250 end if end do 250 continue if( f9dp .eq. 4 ) then write(10,254) 254 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t43,'x',t55,'y',t67,'z',t79,'Px',t91,'Py', & t103,'Pz',t115,'Bx',t127,'By',t139,'Bz',t151,'wt', & t163,'Ex',t175,'Ey',t187,'Ez',t199,'arclength',t211,'polX', & t223,'polY',t235,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e12.4)' else if( f9dp .eq. 6 ) then write(10,256) 256 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t45,'x',t59,'y',t73,'z',t87,'Px',t101,'Py', & t115,'Pz',t129,'Bx',t145,'By',t157,'Bz',t171,'wt', & t185,'Ex',t199,'Ey',t213,'Ez',t227,'arclength',t241,'polX', & t255,'polY',t269,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e14.6)' else if( f9dp .eq. 8 ) then write(10,258) 258 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t47,'x',t63,'y',t79,'z',t95,'Px',t111,'Py', & t127,'Pz',t143,'Bx',t159,'By',t175,'Bz',t191,'wt', & t207,'Ex',t223,'Ey',t239,'Ez',t255,'arclength',t271,'polX', & t287,'polY',t303,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e16.8)' else if( f9dp .eq. 10 ) then write(10,260) 260 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t49,'x',t67,'y',t85,'z',t103,'Px',t121,'Py', & t139,'Pz',t157,'Bx',t175,'By',t193,'Bz',t211,'wt', & t229,'Ex',t247,'Ey',t265,'Ez',t283,'arclength',t301,'polX', & t319,'polY',t337,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e18.10)' else if( f9dp .eq. 12 ) then write(10,262) 262 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t51,'x',t71,'y',t91,'z',t111,'Px',t131,'Py', & t151,'Pz',t171,'Bx',t191,'By',t211,'Bz',t231,'wt', & t251,'Ex',t271,'Ey',t291,'Ez',t311,'arclength',t331,'polX', & t351,'polY',t371,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e20.12)' else if( f9dp .eq. 14 ) then write(10,264) 264 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t53,'x',t75,'y',t97,'z',t119,'Px',t141,'Py', & t163,'Pz',t185,'Bx',t207,'By',t229,'Bz',t251,'wt', & t273,'Ex',t295,'Ey',t317,'Ez',t339,'arclength',t361,'polX', & t383,'polY',t405,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e22.14)' else if( f9dp .eq. 16 ) then write(10,266) 266 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t55,'x',t79,'y',t103,'z',t127,'Px',t151,'Py', & t175,'Pz',t199,'Bx',t223,'By',t247,'Bz',t271,'wt', & t295,'Ex',t319,'Ey',t343,'Ez',t367,'arclength',t391,'polX', & t415,'polY',t439,'polZ') fmt = '(i7,2i4,i5,i6,1x,1p,18e24.16)' else if( f9dp .eq. 17 ) then write(10,267) 267 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t56,'x',t81,'y',t106,'z',t131,'Px',t156,'Py', & t181,'Pz',t206,'Bx',t231,'By',t256,'Bz',t281,'wt', & t306,'Ex',t331,'Ey',t356,'Ez',t381,'arclength',t406,'polX', & t431,'polY',t456,'polZ') fmt = '(i7, 2i4, i5, i6, 1x, 1p, 18e25.16E3)' else write(10,266) fmt = '(i7,2i4,i5,i6,1x,1p,18e24.16)' end if ! backspace(unit=9) ! rewind 1st data line ! 300 continue read(9,*,end=350) ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol ! do i=1,kk if( ievt.eq.igde(i) .and. ipnum.eq.igdp(i) ) then write(10,fmt) ievt,ipnum,iptyp,ipflg,jsrg, & tp,xp,pp,bfld,evtwt,efld,sarc,pol go to 300 ! leave loop if you find good one end if end do ! go to 300 ! 350 continue ! if( kk .gt. 0 ) then write(*,*) 'Set of example good tracks' mxt = MIN(kk,100) do i=1,mxt,10 write(*,'(10i6)') (igde(j),j=i,i+9) end do end if ! close(9) close(10) ! 900 continue end