! Find all phase space files with 3 digit number that indicates turn (using write_phase_space_every_n_turns) ! Read files for turns 0 and turn 50 and save muon number time, pz ! Write new phase space file with muon numbers from file with 50 turns and time, momentum from file 0 or a reference !if nothing on command line check all subdirectories that are named by timestamp !If there is a single directory on the command line read just that one !if there are two directories on the command line use all subdirectories that fall between those timestamps program energy_vs_time use energy_vs_time_mod use energy_vs_time_interface implicit none type(muon_alive_struct), allocatable :: muon_000(:), muon_050(:), muon_save(:) character*300 string, new_string, phase_space_file_1, phase_space_file_2 character*300 dir_file character*300 dir(300) character*20 dir_min, dir_max character*200 cwd character*24000 file_string/' '/, file_string_3/' '/ real(rp) slope, offset logical itexists, first/.true./, first3/.true./ integer reason, i, j, n_left/0/, n_left_0/0/ integer lun, lun2, lun3 integer number_files integer n integer k integer nevents/1000000/ integer min_date,max_date, min_time, max_time, date, time integer nargs integer indexx dir(1)=' ' phase_space_file_1='END000' phase_space_file_2='END140' nargs = command_argument_count() if (nargs == 1)then call cesr_getarg(1,dir(1)) ! the data is all in a single directory input directory name on command line ! print *, 'len(trim(dir(1)))= ',len(trim(dir(1))) number_files=1 print '(i10,1x,a19)', number_files,dir(1) else if( nargs == 2 .or. nargs==4)then call cesr_getarg(1,dir_min) call cesr_getarg(2,dir_max) read(dir_min(1:8),*)min_date read(dir_min(10:15),*)min_time read(dir_max(1:8),*)max_date read(dir_max(10:15),*)max_time if(nargs == 4)then call cesr_getarg(3,phase_space_file_1) call cesr_getarg(4,phase_space_file_2) endif endif print '(a,1x,a)',trim(phase_space_file_1), trim(phase_space_file_2) call getcwd(cwd) !if all directories are to be combined, input nothing on command line string='ls > out.dat' call execute_command_line(string) lun=lunget() open(unit=lun,file='out.dat', status='old') n=0 do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit if( (index(new_string,'2020')==0) .and. (index(new_string,'2019')==0) .and. (index(new_string,'2021')==0) .and. (index(new_string,'2022')== 0).and. (index(new_string,'2023')== 0).and. (index(new_string,'2024')== 0) )cycle n=n+1 dir(n)=trim(new_string) print '(i10,1x,a19)',n, dir(n) number_files = n end do close(lun) endif !_______________________________________________________________________ k=0 allocate(muon_save(1:nevents)) do n=1, number_files if(nargs == 2 .or. nargs == 4)then read(dir(n)(1:8),*)date read(dir(n)(10:15),*)time if(datemax_date)cycle if(time>max_time)cycle endif string='ls '//trim(dir(n))//' > '//trim(dir(n))//'/out.dat' call execute_command_line(string) lun=lunget() inquire (file = trim(dir(n))//'/'//'out.dat', exist = itexists) if(.not. itexists)then print *, trim(dir(n))//'/'//'out.dat', 'does not exist' cycle endif open(unit=lun,file=trim(dir(n))//'/'//'out.dat', status='old') lun2=lunget() open (unit=lun2,file = trim(dir(n))//'/'//'Energy_vs_time_0.dat') if(first)write(lun2,'(a10,1x,4a16,1x,a10)')'muon index', 'momentum','time',' initial disp',' inital angle','ix' first=.false. lun3=lunget() open (unit=lun3,file = trim(dir(n))//'/'//'Energy_vs_time_spin_0.dat') if(first3)then write(lun3,'(20x,a,36x,80x,a)') trim(phase_space_file_1), trim(phase_space_file_2) write(lun3,'(a10,1x,19a16)')'muon index', 'x','px/p','y','py/p','dp/p' ,'time','spin_x','spin_y','spin_z', 'x','px/p','y','py/p','dp/p' ,'time','spin_x','spin_y','spin_z','acos(s.p)' first3=.false. endif do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit indexx= index(new_string,phase_space_file_1) if(index(new_string,phase_space_file_1)/=0 .or. trim(new_string(1:6)) == trim(phase_space_file_1(1:6)))then dir_file = trim(dir(n))//'/'//trim(new_string) call get_data_from_END_phase_space(dir_file, muon_000, n_left_0) elseif(index(new_string,phase_space_file_2)/=0 .or. trim(new_string(1:6)) == trim(phase_space_file_2(1:6)))then dir_file = trim(dir(n))//'/'//trim(new_string) call get_data_from_END_phase_space(dir_file, muon_050,n_left) else cycle endif end do close(unit=lun) if(n_left <=0)cycle do i = 1, n_left do j=1, n_left_0 if(muon_050(i)%muon_numb == muon_000(j)%muon_numb)then if(k+1>nevents) then print *,' more than ',nevents,' events' cycle endif k=k+1 muon_save(k) = muon_000(j) write(lun2,'(i10, 1x, 4es16.8,1x,i16)')muon_000(j)%muon_numb, muon_000(j)%pz, muon_000(j)%time, muon_000(j)%x, muon_000(j)%xp, muon_000(j)%ix write(lun3,'(i10, 1x, 20es16.8)')muon_000(j)%muon_numb, muon_000(j)%x, muon_000(j)%xp,muon_000(j)%y,muon_000(j)%yp, muon_000(j)%pz, muon_000(j)%time, muon_000(j)%spin, & muon_050(i)%x, muon_050(i)%xp,muon_050(i)%y,muon_050(i)%yp, muon_050(i)%pz, muon_050(i)%time, muon_050(i)%spin, muon_050(i)%phi exit endif end do end do n_left = 0 close(unit=lun2) close(unit=lun3) file_string = trim(file_string)//' '//trim(dir(n))//'/'//'Energy_vs_time_0.dat' file_string_3 = trim(file_string_3)//' '//trim(dir(n))//'/'//'Energy_vs_time_spin_0.dat' if(any(dir(1:n-1) == dir(n)))then print *, ' directory ', n, ' same as ', ' earlier' endif print '(i10, 1x,a)', n, trim(dir(n)) end do call fit_line(muon_save, k, slope, offset) print '(a)','cat '//trim(file_string)//' > all_Energy_vs_time_0.dat' call execute_command_line('cat '//trim(file_string)//' > all_Energy_vs_time_0.dat') print '(a)','cat '//trim(file_string_3)//' > all_Energy_vs_time_spin_.dat' call execute_command_line('cat '//trim(file_string_3)//' > all_Energy_vs_time_spin.dat') end program energy_vs_time subroutine get_data_from_END_phase_space(dir_file,muon_alive, n_left) use energy_vs_time_mod use energy_vs_time_interface, dummy_except => get_data_from_END_phase_space implicit none type(muon_alive_struct), allocatable :: muon_alive(:) character*300 dir_file character*600 long_string logical itexists real(rp) x, pz, time,disp, ang real(rp) disp_y, ang_y real(rp) spin(3) real(rp) acos_sp integer ix,i,n,muon_number, n_left if(.not. allocated(muon_alive))allocate(muon_alive(1:1000000)) inquire (file=dir_file, exist = itexists) if(.not. itexists)return print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) n=0 do while(.true.) read(5,'(a)',end=99)long_string if(index(long_string,'muon')/=0)cycle i=0 ix=0 do while (i<15) call string_trim(long_string(ix+1:), long_string, ix) i=i+1 read(long_string(1:ix),*)x if(i == 1) then muon_number = int(x) n=n+1 elseif(i==4)then disp = x elseif(i==5)then ang = x elseif(i==6)then disp_y = x elseif(i==7)then ang_y = x elseif( i == 10)then time = x elseif( i == 9)then pz = x elseif( i == 12)then spin(1) = x elseif( i == 13)then spin(2) = x elseif( i == 14)then spin(3) = x elseif( i == 15)then acos_sp = x endif end do muon_alive(n)%muon_numb = muon_number muon_alive(n)%time = time muon_alive(n)%pz = pz muon_alive(n)%x= disp muon_alive(n)%yp = ang_y muon_alive(n)%y= disp_y muon_alive(n)%xp= ang muon_alive(n)%spin=spin muon_alive(n)%phi = acos_sp n_left = n end do 99 continue close(unit=5) return end subroutine get_data_from_END_phase_space subroutine fit_line(muon_save, k, slope, offset) use sim_utils use energy_vs_time_mod use energy_vs_time_interface, dummy_except => fit_line implicit none type(muon_alive_struct), allocatable :: muon_save(:) integer k,ix integer lun real(rp) slope, offset real(rp) sum_xp, sum_pp, sum_x, sum_p,sum_xsteps, det lun=lunget() open(unit=lun,file='line_fit.dat') ! then try x vs p sum_xp = sum( muon_save(:)%pz * muon_save(:)%x) sum_pp = sum(muon_save(:)%pz*muon_save(:)%pz) sum_x = sum(muon_save(:)%x) sum_p = sum(muon_save(:)%pz) det = k * sum_pp - sum_p**2 slope = (k*sum_xp - sum_x*sum_p)/det offset = (-sum_p * sum_xp + sum_pp * sum_x)/det print '(4(a,es12.4),a,i10)',' sum_xp =',sum_xp,' sum_pp =',sum_pp,' sum_x =',sum_x,' sum_p=',sum_p,' k=', k print '((a,es12.4))','det =', det print '(a,es12.4)', 'slope = ', slope, '; offset = ', offset write(lun,'(2(a,es12.4))') 'slope_x =',slope,'; offset_x= ', offset ! then try xp vs p sum_xp = sum( muon_save(:)%pz * muon_save(:)%xp) sum_pp = sum(muon_save(:)%pz*muon_save(:)%pz) sum_x = sum(muon_save(:)%xp) sum_p = sum(muon_save(:)%pz) det = k * sum_pp - sum_p**2 slope = (k*sum_xp - sum_x*sum_p)/det offset = (-sum_p * sum_xp + sum_pp * sum_x)/det print '(4(a,es12.4),a,i10)',' sum_xp =',sum_xp,' sum_pp =',sum_pp,' sum_x =',sum_x,' sum_p=',sum_p,' k=', k print '((a,es12.4))','det =', det print '(a,es12.4)', 'slope = ', slope, '; offset = ', offset write(lun,'(2(a,es12.4))') 'slope_xp =',slope,'; offset_xp= ', offset close(lun) return end subroutine fit_line