module mia_input use mia_types use haw_toolbox use cesr_utils use tbt_mod use file_number_mod use tbt_struct type offset real(rp) :: nextS !S position where turn offset changes next integer :: offsetBy !# turns to offset by end type offset contains subroutine read_data(data,iset,nset) !Subroutine opens a file of bpm data and reads it into the module. !Routine reads in how many processors are active and assigns that !number to the variable bpmproc implicit none type(data_set) data !Data set integer :: iset, nset !Number of current set integer :: num_turns, num_bpms,num_bunches character(30) :: input integer :: i, j, k real(rp) :: sumh,sumv,aveh,avev character(30) :: Plane(2) !See next line: data Plane/' data file 1.',' data file 2.'/ type (tbt_all_data_struct) :: tbt logical err_flag if (fileread(iset)==.false.) then write (*, '(1x, 3a)', advance = "no") "Enter number for ", Plane(iset) read *, input if (is_number(trim(adjustl(input)))) then Read(input(1:len_trim(input)),*) data%filenum else Print *, trim(adjustl(input)), " is not a valid number!" call exit(1) endif write (*,'(1x, a48)', advance = "no") "Is this a paired file? (Y/N) " read *, input if(input /= 'Y' .and. input /= 'y') then call open_tbt_file(data%filenum, .true., tbt) else call open_tbt_file(data%filenum-1, .true., tbt) call open_tbt_file(data%filenum, .false., tbt) end if else if (isPaired == .true.) then call open_tbt_file(data%filenum-1, .true., tbt) call open_tbt_file(data%filenum, .false., tbt) else call open_tbt_file(data%filenum, .true., tbt) end if num_bpms = 0 do i = 0, 150 if (tbt%bunch(bunchOfInterest)%bpm(i)%data_ok == .true.) num_bpms = num_bpms + 1 end do print *, "Num bpms = ", num_bpms num_bunches = size(tbt%bunch) print *,"Num bunches = ", num_bunches num_turns = size(tbt%shaker) print *,"Num turns = ", num_turns if (.not. ALLOCATED(data%proc)) call allocate_data(data,num_turns,num_bpms) data%numturns = num_turns data%bpmproc = num_bpms k = 0 do i = 0, 150 if (tbt%bunch(bunchOfInterest)%bpm(i)%data_ok == .false.) cycle k = k + 1 data%proc(k)%label = tbt%bpm(i)%name !Calculate and subtract average from data; copy into data structure sumh = 0.; sumv = 0. do j = 1, num_turns sumh = sumh + tbt%bunch(bunchOfInterest)%bpm(i)%turn(j)%x sumv = sumv + tbt%bunch(bunchOfInterest)%bpm(i)%turn(j)%y end do aveh = sumh / num_turns avev = sumv / num_turns do j = 1, num_turns !Factor of 1000 to convert m to mm data%poshis(j,2*k-1) = 1000. * (tbt%bunch(bunchOfInterest)%bpm(i)%turn(j)%x - aveh) data%poshis(j,2*k) = 1000. * (tbt%bunch(bunchOfInterest)%bpm(i)%turn(j)%y - avev) end do end do end subroutine read_data subroutine open_tbt_file(number, initiate, tbt) integer number type (tbt_all_data_struct) :: tbt logical initiate, err_flag Print *, "Attempting to open file #", number call tbt_read_data_file(number, gain_and_pedestal_correct$, initiate, tbt, err_flag) if (err_flag) call exit(1) end subroutine open_tbt_file subroutine read_orbit_data(data1,data2) use bmad use cesr_read_data_mod type(data_set) data1, data2 !Data set integer :: num_orbits = 0 !Number of orbits to read integer, allocatable :: orbitNumbers(:) !Array of orbit numbers integer :: i, j, k, file_num character(140) filein logical err_flag type (butns_struct), allocatable :: butns(:) real(rp) sumh, sumv, aveh, avev type (lat_struct), target :: lat character(200) :: lat_file = '/nfs/cesr/online/machine_data/lattice/cesr/bmad/chess-u_6000mev_20190904.lat' type (db_struct) :: cesr, db call bmad_parser (lat_file, lat) call bmad_to_cesr (lat, cesr, err_flag) call read_orbit_list(num_orbits, orbitNumbers) allocate(butns(num_orbits)) do i = 1, num_orbits call number_to_file_name (orbitNumbers(i), 'orbit', filein, file_num, err_flag) if (err_flag) return ! Error call read_butns_file (filein, butns(i), db, err_flag) if (err_flag) return end do !Assume every orbit has same number of BPMs num_bpms = 0 do i = 0, 150 if (butns(1)%det(i)%ok == .true.) num_bpms = num_bpms + 1 end do print *, "Num bpms = ", num_bpms num_turns = num_orbits print *,"Num orbits = ", num_turns if (.not. ALLOCATED(data1%proc)) call allocate_data(data1,num_turns,num_bpms) data1%numturns = num_turns data1%bpmproc = num_bpms if (.not. ALLOCATED(data2%proc)) call allocate_data(data2,num_turns,num_bpms) data2%numturns = num_turns data2%bpmproc = num_bpms k = 0 do i = 0, 150 if (butns(1)%det(i)%ok == .false.) cycle k = k+1 data1%proc(k)%label = cesr%detector(i)%bmad_name data2%proc(k)%label = cesr%detector(i)%bmad_name data1%proc(k)%eleNum = i data2%proc(k)%eleNum = i !Calculate and subtract average from data; copy into data structure sumh = 0.; sumv = 0. do j = 1, num_turns sumh = sumh + butns(j)%det(i)%x_orb sumv = sumv + butns(j)%det(i)%y_orb end do aveh = sumh / num_turns avev = sumv / num_turns do j = 1, num_turns !Factor of 1000 to convert m to mm data1%poshis(j,2*k-1) = 1000. * (butns(j)%det(i)%x_orb - aveh) data1%poshis(j,2*k) = 1000. * (butns(j)%det(i)%y_orb - avev) data2%poshis(j,2*k-1) = 1000. * (butns(j)%det(i)%x_orb - aveh) data2%poshis(j,2*k) = 1000. * (butns(j)%det(i)%y_orb - avev) end do end do deallocate(orbitNumbers, butns) end subroutine read_orbit_data subroutine read_orbit_list(num_orbits, orbitNumbers) integer :: num_orbits !Number of orbits to read integer, allocatable :: orbitNumbers(:) !Array of orbit numbers character(40) :: longstring, first, second integer ::loc, iostatus, number1, number2, counter, j if(orbitList == 'none') then write (*, '(a)', advance = "no") "Enter filename containing list of orbit numbers: " read *, orbitList endif open (10,file=trim(adjustl(orbitList))) num_orbits = 0 do read(10,'(a)',iostat=iostatus) longstring !Read the whole line into a single character variable if (iostatus/=0) exit read(longstring,*) first !Use an internal read to get the first variable loc=index(longstring,',') !Get the location of the comma after the first parameter if(loc == 0) then num_orbits = num_orbits+1 else second=longstring(loc+1:len_trim(longstring)) !Make a new string without the first parameter read(first(1:len_trim(first)),*) number1 read(second(1:len_trim(second)),*) number2 num_orbits = num_orbits+(number2-number1)+1 endif end do close(10) open (11,file=trim(adjustl(orbitList))) allocate(orbitNumbers(num_orbits)) counter = 0 do read(11,'(a)',iostat=iostatus) longstring !Read the whole line into a single character variable if (iostatus/=0) exit read(longstring,*) first !Use an internal read to get the first variable read(first(1:len_trim(first)),*) number1 loc=index(longstring,',') !Get the location of the comma after the first parameter if(loc == 0) then counter = counter+1 orbitNumbers(counter) = number1 else second=longstring(loc+1:len_trim(longstring)) !Make a new string without the first parameter read(second(1:len_trim(second)),*) number2 do j = number1, number2 counter = counter+1 orbitNumbers(counter) = j end do endif end do close(11) end subroutine read_orbit_list subroutine getTurnOffset(data) type(data_set) :: data character(120) :: line integer :: inputstatus, currentBunch,& numOffsets, & commaComa logical :: endOfTheLine type(offset) :: offsets(5) integer :: i numOffsets = 1 do i=1, 5 offsets(i)%nextS = 8314.99 offsets(i)%offsetBy = 0 end do open (unit = 14, file = './data/turnoffset.dat', status = "old", & iostat = inputstatus) Read (14, *, iostat = inputstatus) line if (inputstatus > 0) stop "*** INPUT ERROR ***" if (inputstatus < 0) then Print *, "Could not read turn offset file (turnoffset.dat)" Print *, "If this is machine data you may see a phase error" return end if !Do loop parses a line that has been read in, and !reads in the next line at the end. do while (inputstatus == 0) write (currentBunch,*), line(1:scan(line,',')) !Remove portion just read in: line = adjustl(line(scan(line,','):len(line))) !Find next comma (or set to end of line in case there are no more commas): commaComa = scan(line,',') if (commaComa == 0) then commaComa = len_trim(line) endOfTheLine = .true. end if ! if (currentBunch == data(1)%bunch) then ! write (offsets(numOffsets),*), ! end if ! if (currentBunch == data(2)%bunch) then ! end if Read (14, *, iostat = inputstatus) line end do offsets(1)%nextS = 502.5 offsets(1)%offsetBy = -1 offsets(2)%nextS = 90000.4 offsets(2)%offsetBy = 1 ! call correctTurnOffset(data, offsets) end subroutine getTurnOffset subroutine correctTurnOffset(data, lastTime) type(data_set) :: data logical :: lastTime !True if file 2 has been read in ! type(offset) :: offsets(*) integer :: offsetAmount, newNumTurns real :: offsetStart real(rp), allocatable :: tempPosHis(:,:) integer :: i, j, k, count, & offsetCount, & startTurn, endTurn offsetStart = 502.5 offsetAmount = 1 newNumTurns = NUM_TURNS - offsetAmount !Copy posHis into a temporary matrix before resizing allocate(tempPosHis(NUM_TURNS,2*NUM_BPMS)) do i=1, NUM_TURNS do j=1, 2*NUM_BPMS tempPosHis(i,j) = data%poshis(i,j) end do end do deallocate(data%poshis) allocate(data%poshis(newNumTurns, 2*NUM_BPMS)) !Here comes a shift.... do i=1, NUM_BPMS if (data_struc%proc(i)%sPos >= offsetStart) then !Cut turns from beginning startTurn = 1+offsetAmount endTurn = newNumTurns+offsetAmount else !Cut turns from end startTurn = 1 endTurn = newNumTurns end if count = 1 do j = startTurn, endTurn data%poshis(count,2*i-1) = tempPosHis(j,2*i-1) data%poshis(count,2*i) = tempPosHis(j,2*i) count = count+1 end do end do data%numturns = newNumTurns !Change NUM_TURNS only on last call to this subroutine if (lastTime) then NUM_TURNS = newNumTurns Print *, "new num turns: ", NUM_TURNS end if end subroutine correctTurnOffset subroutine getLatticeAndTune ! !Gets lattice name and integer tunes for the current set of files. !Not required; these are only for the CESRV output file !and to check for phase errors. The integer tune is useful !for MIA to check for large gaps between detectors which may !cause phase errors (if the phase difference is greater than 2*pi !MIA will not have the correct total phase.) ! character(120) :: line integer :: inputstatus open (unit = 81, file = './data/currentMiaLat.dat', status="old", & iostat = inputstatus) if (inputstatus > 0) stop "*** INPUT ERROR ***" if (inputstatus < 0) then Print *, "Could not find lattice name. If you are importing this", & "data into CESRV, you will need to add the lattice name manually." latticeName = " " integerTune(1) = 0 integerTune(2) = 0 return end if Read (81, *, iostat = inputstatus) line latticeName = trim(adjustl(line)) Read (81,'(2i)'), integerTune(1) Read (81,'(2i)'), integerTune(2) end subroutine getLatticeAndTune subroutine bpm_ops(data,iset, nset, first_run) ! !Calls functions that locate and match BPMs !Remove if bpm numbers, pairs, etc. are not obtained !from input files (ex. from CESRV or another program) ! type(data_set) data(*) integer :: iset, nset logical :: first_run Print *, "Calling locate_bpm" call locate_bpm (data(iset)) !****This is called for both data files--not necessary if (iset > 1) then if (first_run) then call findPair(nset,data) end if call match_processors (data) call get_ele_sPos(data) endif end subroutine bpm_ops subroutine locate_bpm (data) ! !Creates a list of the bpms that are active. Keeps track of bpms in West !or East !Assigns BPMs a number in addition to their label. This is positive !in the west and negative in the east ! implicit none type(data_set)data !Data set integer :: n_bpm !Number of BPMs Print *, "You have entered locate_bpm" do n_bpm = 1, NUM_BPMS call parseDet(data%proc(n_bpm)%label, data%proc(n_bpm)%number,& data%proc(n_bpm)%is_west) if (n_bpm > 1) then !Check for extra bpms (ex BPM 101 which comes after 12W) ! Switch to cesrv bpm index: !data%proc(n_bpm)%eleNum = cesrv_det_idx(data%proc(n_bpm)%eleNum) if (data%proc(n_bpm)%number>100 .and. data%proc(n_bpm-1)%number<50) then Print *, "locate_bpm: possibly a 100+ BPM. Label: ", data%proc(n_bpm)%label data%proc(n_bpm)%number = data%proc(n_bpm-1)%number + 0.2 data%proc(n_bpm)%is_west = data%proc(n_bpm-1)%is_west end if end if end do end subroutine locate_bpm subroutine parseDet(label, number, isWest) ! !Parses det label into a number for easier comparison !Should get rid of isWest someday.... ! integer :: nchar, & !Place were a certain string was found. n_bpm,& !Number of BPMs i !Counters character(5) :: temp !Holds temporary file number character (13) :: label real(rp) :: number logical :: isWest temp = trim(adjustl(label)) if ((scan(temp,"W")+scan(temp,"w")) > 0) then isWest = .true. else if ((scan(temp,"E")+scan(temp,"e")) > 0) then isWest = .false. else isWest = .true. !Extra detectors without E/W label are on the west side, I believe Print *, "BONUS DETECTOR DETECTED! Label: ", temp end if number = extract_int(label) !Check for BPM0E and change its number to be different than 0W: if (number == 0 .and. & .not. isWest) then number = -0.001 end if end subroutine parseDet subroutine findPair(nset,data) ! !Uses list of bpms from prior subroutine to find l's that are included in !the data file being studied ! implicit none integer :: inputstatus, openstatus,& !Status indicators num, i, j, i_file, nset, k, li, &!Counters and nset--number of sets nfile, rnum, q !Counter and actual number of BPM pairs logical :: found_one(2) !If a BPM pair has been found type(data_set) data(*) !All data type(known_spacings), allocatable :: bpm_old(:) !BPM data read in character (30) :: temp !knownl.inp contains a list of BPMs with known spacing. !Change to accept different locations for knownl.inp... open (unit = 2, file = "./knownl.inp", status = "old", iostat = openstatus) if (openstatus .le. 0) print *, "Opened local copy of knownl.inp" if (openstatus > 0) then open (unit = 2, file = "/nfs/cesr/offline/instr/anal/cbpmII/tune/knownl.inp", status = "old", iostat = openstatus) if (openstatus .le. 0) print *, "Opened central copy of knownl.inp" if (openstatus > 0) Stop "*** Cannot open file knownl.inp ***" end if 60 format (1x, i5,/) 61 format (1x, a13, 3x, a13, 2x, f7.3) read (2,60, iostat = inputstatus) num if (inputstatus > 0) stop "*** INPUT ERROR (knownl.inp)***" if (inputstatus < 0) stop "*** Not enough data (knownl.inp)***" allocate (bpm_old(num)) call allocate_bpm_pairs(num) bpm_old(1)%number = num !Change bpm_pairs%number not to be in every bpm of the array? bpm_old(:)%in_use = .false. rnum = 0 !Actual number of BPMs in use do i = 1, num read (2,61, iostat = inputstatus) (bpm_old(i)%bpm_name(j),j=1,2), & bpm_old(i)%length if (inputstatus > 0) stop "*** INPUT ERROR (knownl.inp)***" if (inputstatus < 0) stop "*** Not enough data (knownl.inp)***" bpm_old(i)%in_use = .true. do q=1,2 temp = bpm_old(i)%bpm_name(q) temp = trim(adjustl(temp(9:len_trim(temp)))) bpm_old(i)%bpm_name(q) = temp end do end do !This could be made more efficient by finding a BPM number for !the pairs and scanning for that. !Also, this checks all BPMs vs all pairs--change to stop when a !match is found. bpm_old(1)%has_one = .true. do i_file = 1, nset !File do i = 1, num !BPM pairs bpm_old(i)%bpm_pntr(1)=0 bpm_old(i)%bpm_pntr(2)=0 ! do j = 1, NUM_BPMS !All BPMs do j = 1, data(i_file)%bpmproc !All BPMs in data set do k = 1, 2 !First and second of the BPM pair if (data(i_file)%proc(j)%label == bpm_old(i)%bpm_name(k)) then bpm_old(i)%bpm_pntr(k) = j end if end do end do if (bpm_old(i)%bpm_pntr(1) == 0 .or. & bpm_old(i)%bpm_pntr(2) == 0) then bpm_old(i)%in_use = .false. else bpm_old(i)%in_use = .true. found_one(i_file) = .true. end if end do if (.not. found_one(i_file)) then Print *, "**** Warning: No detector pairs were found ****" Print *, "**** Analysis will not be complete ****" end if !Unnecessary? bpm_old(1)%has_one = bpm_old(1)%has_one .and. found_one(i_file) end do !Copies BPMs pairs that were found into bpm_pairs. !Not all BPMs in knownl.inp may be in use. ! do nfile = 1, nset !File do li = 1, num !BPM pairs from knownl.inp if (bpm_old(li)%in_use == .true.) then rnum = rnum+1 !Files 1 and 2 should be the same !Remove file from bpm_pairs? !Can check consistancy in previous do loops. bpm_pairs(rnum) = bpm_old(li) bpm_pairs(rnum) = bpm_old(li) endif enddo if (rnum > 0) then bpm_pairs(1)%has_one = .true. else bpm_pairs(1)%has_one = .false. endif ! enddo bpm_pairs(1)%number = rnum deallocate (bpm_old) close(2) end subroutine findPair subroutine get_ele_sPos(data) ! !Reads in a file containing bpms with their s positions and !element numbers and matches them with bpms in the input file. ! type(data_set) :: data(*) character(10) :: detName integer :: eleNum, openstatus, i, inStat, bpmInt real(rp) :: lastNum real(rp) :: sPos, bpmNum logical :: continue, theEnd noSPos = .false. open (unit = 273, file = "./one_ring.det", status = "old", iostat = openstatus) if (openstatus .le. 0) print *, "Opened local copy of one_ring.det" if (openstatus > 0) then open (unit = 273, file = "/nfs/cesr/offline/instr/anal/cbpmII/tune/one_ring.det", status = "old", iostat = openstatus) if (openstatus .le. 0) print *, "Opened central copy of one_ring.det" if (openstatus > 0) then Print *, "*** Cannot open file one_ring.det ***" Print *, "*** No output file for CESRV will be generated. ***" noSPos = .true. return end if end if 188 format (1x,a7,5x,f10.6,5x,i3) 56 format (2x,i2) theEnd=.false. !It's just the beginning continue = .true. do while (continue) read (273,188,iostat=inStat) detName, sPos, eleNum if (inStat>0) stop "*** Error reading one_ring.det ***" if (inStat < 0) then continue = .false. cycle endif detName = trim(adjustL(detName)) print *, detName select case(detName(1:2)) case ('EN') !Location of the end of the ring: the length of the machine theEnd=.true. endLoc=sPos case ('DE') read (detName,56) bpmInt bpmNum = bpmInt if (detName(8:8) == "A") then bpmNum = bpmNum+0.5 endif if ( detName(7:7) == "E" .or. detName(5:5) == "E") then bpmNum = -bpmNum if (bpmNum == 0) then bpmNum = -0.001 !Since there is no -0...something close. endif endif case default Print *, "Error in reading BPM S Positions" exit end select if (.not. theEnd) then do i=1,num_bpms if (data_struc%proc(i)%number == bpmNum) then data_struc%proc(i)%sPos = sPos data_struc%proc(i)%eleNum = eleNum if (i <= data(1)%bpmproc) then data(1)%proc(i)%sPos = sPos data(1)%proc(i)%eleNum = eleNum end if if (i <= data(2)%bpmproc) then data(2)%proc(i)%sPos = sPos data(2)%proc(i)%eleNum = eleNum end if end if end do end if print *, data(1)%proc(i)%eleNum end do call check_sPos() end subroutine get_ele_sPos subroutine check_sPos() ! !Looks for detectors whose position and element # were not found by !the get_ele_sPos subroutine. Position is only used for plotting, !so if a BPM is not found its position is approximated as being !halfway between the two adjacent detectors. The element number !passed to CESRV will be 0. ! integer :: i do i=1, NUM_BPMS if (data_struc%proc(i)%sPos == -999) then Print *, "Position and element number not found for detector ", & data_struc%proc(i)%label if (i>1) then data_struc%proc(i)%eleNum = 0 data_struc%proc(i)%sPos = (data_struc%proc(i-1)%sPos - & data_struc%proc(i+1)%sPos) / 2.0 end if end if end do end subroutine check_sPos subroutine match_processors(data) ! !This routine makes sure the data being compared from !each file has matching processors (if there are multiple files) ! integer :: nbpm,& !BPM number n_file_1, n_file_2, & !Counters for # of BPMs in each file i !counter type(data_set) :: data(*) !Data set character(120) :: vetoString character(10) :: detName integer :: length n_file_1 = 1 n_file_2 = 1 nbpm = 1 do while (n_file_1 <= data(1)%bpmproc) if (data(1)%proc(n_file_1)%label == & data(2)%proc(n_file_2)%label) then data_struc%proc(nbpm) = data(1)%proc(n_file_1) n_file_1 = n_file_1 + 1 n_file_2 = n_file_2 + 1 nbpm = nbpm + 1 else do i = n_file_2 +1, data(2)%bpmproc if (data(1)%proc(n_file_1)%label == & data(2)%proc(i)%label) then data_struc%proc(nbpm) = data(1)%proc(n_file_1) n_file_1 = n_file_1 + 1 n_file_2 = i + 1 nbpm = nbpm + 1 cycle else !Test for detectors not matched and !put them on the veto list if (i == data(2)%bpmproc) then Print *, "Detector ", data(1)%proc(n_file_1)%label,& " was not found in file 2. Vetoing this detector." end if end if end do n_file_1 = n_file_1 + 1 end if end do Print *, data(1)%bpmproc, data(2)%bpmproc, "nbpm: ", nbpm if (data(1)%bpmproc /= nbpm - 1) then print *, "Number of Processors don't match; Mode not supported" stop else do i=1, data(1)%bpmproc data_struc%proc(i)%label = data(1)%proc(i)%label data_struc%proc(i)%number = data(1)%proc(i)%number end do end if end subroutine match_processors end module mia_input