! Find all phase space files ! -command line input: first two arguments are initial and final timestamped subdirectory. ! To get all dated subdirectories '00000000_000000' '99999999_999999' ! The next narg-2 argmuments are the phase_space files to be combined (anything with 'phase_space" in the name) ! -namelist input: if there are zero arguments on the command line, read data from 'energy_spin_vs_time_input.dat' (an example is appended to this file program energy_spin_vs_time use nr use energy_vs_time_mod use energy_vs_time_interface implicit none type muon_at_struct type(muon_alive_struct), allocatable :: muon_save(:) end type muon_at_struct type(muon_at_struct) muon_at(20) type(muon_alive_struct)muon character*300 string, new_string, phase_space_file(20) character*300 dir_file character*300 dir(300) character*20 dir_min, dir_max character*200 cwd character*24000 file_string/' '/, file_string_3/' '/, file_string_4/' '/, file_string_5/' '/, file_string_6/' '/, file_string_7/' '/ character*30 format_string character*100 heading(20) character*100 trimmed_file character*300 line real(rp) slope, offset real(rp) a,b,chi2, siga, sigb, q, phase(20), t(20), err(20) real(rp) phi_ave real(rp) vec5/0./,location/0./, jx/0./, jy/0./ real(rp) phase_sum(20), time_sum(20) real(rp) delta_ix_mu logical itexists, first/.true./, first3/.true./ logical write_phi_fixed/.true./, first_write/.true./ logical op integer reason, i, j, n_left(20)/20*0/ integer lun, lun2, lun3, lun4, lun5,lun6/0/ integer number_files integer n integer k integer nevents/1000000/ integer min_date,max_date, min_time, max_time, date, time integer nargs integer ix_mu(20),ix integer ix_numb integer, allocatable:: muon_number(:) integer ind integer total integer ios integer number(20) namelist/energy_spin_vs_time_input/nargs, dir_min, dir_max, phase_space_file dir(1)=' ' phase_space_file(1)='END000' phase_space_file(2)='END028' 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) if(nargs >= 4)then ix=nargs-2 do i=3,nargs call cesr_getarg(i,phase_space_file(i-2)) end do endif endif if(nargs == 0) then OPEN (UNIT=5, FILE='energy_spin_vs_time_input.dat', STATUS='old', IOSTAT=ios) READ (5, NML=energy_spin_vs_time_input, IOSTAT=ios) print *, 'ios=', ios rewind(unit=5) CLOSE(5) WRITE(6,NML=energy_spin_vs_time_input) ix=nargs-2 endif if(nargs == 2 .or. nargs >=4)then 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 endif phase_sum(1:ix)=0 time_sum(1:ix)=0 number(1:ix)=0 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_number(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(nargs >4)then lun4=lunget() open(lun4,file=trim(dir(n))//'/'//'Energy_spin_time.dat') lun5=lunget() open (lun5,file=trim(dir(n))//'/'//'omega_a.dat') endif if(first3)then if(nargs >=4)then write(lun3,'(20x,a,36x,80x,a)') trim(phase_space_file(1)), trim(phase_space_file(2)) write(lun3,'(a10,1x,20a16)')'muon index', 'x','px/p','y','py/p','dp/p' ,'time','spin_x','spin_y','spin_z','acos(s.p)', 'x','px/p','y','py/p','dp/p' ,'time','spin_x','spin_y','spin_z','acos(s.p)' write(lun4,'(a1,a)') ('#',trim(phase_space_file(i)), i=1,ix) if(ix<10)write(format_string,'("(a1,i9,",i1,"(6i12))")')ix if(ix>=10)write(format_string,'("(a1,i9,",i2,"(6i12))")')ix write(lun4, format_string)'#',(j, j=1,6*ix+1) heading(1:20)= 'dp/p time spin_x spin_y spin_z acos(s.p)' if(ix<10)write(format_string,'("(a10,",i1,"(6a12))")' )ix if(ix>=10)write(format_string,'("(a10,a12,",i2,"(6a12))")' )ix write(lun4, format_string)'muon index','omega_a (b)', (heading(j)(1:4), heading(j)(6:9), heading(j)(11:16), heading(j)(18:23), heading(j)(25:30), heading(j)(32:40), j=1,ix) write(lun5,'(a10,a12,a16,17a12)')'muon_numb','offset','slope','chi2','phi','spin_x','spin_y','spin_z','phi','spin_x','spin_y','spin_z','phi','spin_x','spin_y','spin_z','phi','spin_x','spin_y','spin_z' endif first3=.false. endif !first3 do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit do i = 1,ix if(index(phase_space_file(i),new_string)/=0)then dir_file = trim(dir(n))//'/'//trim(new_string) if(i == ix-1)then line='cat '//trim(dir_file)//' > all_'//trim(new_string) print '(a)', line call execute_command_line(line) endif call get_data_from_END_phase_space(dir_file, muon_at(i)%muon_save, n_left(i)) phi_ave = sum(muon_at(i)%muon_save(1:n_left(i))%phi,1,muon_at(i)%muon_save(1:n_left(i))%spin(3)/=0) total=0 do j=1,n_left(i) if(muon_at(i)%muon_save(j)%spin(3)/=0 .and. muon_at(i)%muon_save(j)%spin(2)/=0) total = total +1 end do phi_ave = phi_ave/total first_write = .true. do j=1,n_left(i) ! if(muon_at(i)%muon_save(j)%spin(3) /= 0)write(112,'(es12.4)')muon_at(i)%muon_save(j)%phi if(muon_at(i)%muon_save(j)%phi - phi_ave > twopi/2)muon_at(i)%muon_save(j)%phi = muon_at(i)%muon_save(j)%phi-twopi if(muon_at(i)%muon_save(j)%phi - phi_ave < -twopi/2)muon_at(i)%muon_save(j)%phi = muon_at(i)%muon_save(j)%phi+twopi if(write_phi_fixed)then if(first_write)then inquire(lun6,opened=op) if(op)close(lun6) lun6=lunget() trimmed_file = trim(trim(dir_file)//'_trimmed') open(lun6,file = trimmed_file) write(lun6,'(a10,15a12)')'muon','jx','jy','x','xp','y','yp','vec5','dp/p','time','location','spin_x','spin_y','spin_z','phi' first_write = .false. endif muon=muon_at(i)%muon_save(j) write(lun6,'(i10,15es12.4)')muon%muon_numb, jx, jy, muon%x, muon%xp, muon%y,muon%yp, vec5,muon%pz,muon%time,location,muon%spin, muon%phi endif end do if(write_phi_fixed)close(lun6) exit endif end do end do close(unit=lun) if(any(n_left(1:ix) <=0))then j=minloc(n_left(1:ix) ,1) print '(a,a,a1,a)',' no events in file ',dir(n),'/',trim(phase_space_file(j)) print *,trim(phase_space_file(j)) stop endif do j=1,ix-1 muon_at(j)%muon_save(:)%mask=.true. end do do i = 1, n_left(ix) ix_numb = muon_at(ix)%muon_save(i)%muon_numb ix_mu(ix) = i do j=1,ix-1 ix_mu(j)=minloc(abs(muon_at(j)%muon_save(:)%muon_numb - ix_numb),1,muon_at(j)%muon_save(:)%mask) muon_at(j)%muon_save(ix_mu(j))%mask = .false. delta_ix_mu = muon_at(j)%muon_save(ix_mu(j))%muon_numb - ix_numb if(delta_ix_mu /= 0)then print *, 'muon numbers out of synch' stop endif end do if(k+1>nevents) then print *,' more than ',nevents,' events' cycle endif k=k+1 muon_number(k) = muon_at(1)%muon_save(i)%muon_numb ind = ix_mu(1) write(lun2,'(i10, 1x, 4es16.8,1x,i16)')muon_at(1)%muon_save(ind)%muon_numb, muon_at(1)%muon_save(ind)%pz, & muon_at(1)%muon_save(ind)%time, muon_at(1)%muon_save(ind)%x, muon_at(1)%muon_save(ind)%xp, muon_at(1)%muon_save(ind)%ix write(lun3,'(i10, 1x, 20es16.8)')ix_numb, (muon_at(j)%muon_save(ix_mu(j))%x, & muon_at(j)%muon_save(ix_mu(j))%xp,muon_at(j)%muon_save(ix_mu(j))%y,muon_at(j)%muon_save(ix_mu(j))%yp, muon_at(j)%muon_save(ix_mu(j))%pz, & muon_at(j)%muon_save(ix_mu(j))%time, muon_at(j)%muon_save(ix_mu(j))%spin, muon_at(j)%muon_save(ix_mu(j))%phi, j=1,2) do j=1,ix t(j) = muon_at(j)%muon_save(ix_mu(j))%time phase(j) = muon_at(j)%muon_save(ix_mu(j))%phi end do call fit(t(3:ix),phase(3:ix),a,b,siga,sigb,chi2,q) do j=1,ix if(chi2 >0 .and. chi2 < 10)then phase_sum(j) = phase_sum(j)+ phase(j) time_sum(j) = time_sum(j) + t(j) number(j) = number(j)+1 endif end do write(format_string,'("(i10,es12.4,",i1,"(6es12.4))")') ix write(lun4, format_string) ix_numb, b, & (muon_at(j)%muon_save(ix_mu(j))%pz,muon_at(j)%muon_save(ix_mu(j))%time, muon_at(j)%muon_save(ix_mu(j))%spin, muon_at(j)%muon_save(ix_mu(j))%phi , j=1,ix) write(lun5,'(i10,es12.4,es16.8,17es12.4)')ix_numb, a,b,chi2,phase(1), muon_at(1)%muon_save(ix_mu(1))%spin, phase(2), muon_at(2)%muon_save(ix_mu(2))%spin, & phase(3), muon_at(3)%muon_save(ix_mu(3))%spin, phase(4), muon_at(4)%muon_save(ix_mu(4))%spin end do n_left = 0 close(unit=lun2) close(unit=lun3) close(unit=lun4) close(unit=lun5) if(ix == 2)then file_string = trim(file_string)//' '//trim(dir(n))//'/'//'Energy_vs_time_0.dat' endif if(ix >=2)then file_string_3 = trim(file_string_3)//' '//trim(dir(n))//'/'//'Energy_vs_time_spin_0.dat' file_string_4 = trim(file_string_4)//' '//trim(dir(n))//'/'//'Energy_spin_time.dat' file_string_5 = trim(file_string_5)//' '//trim(dir(n))//'/'//'omega_a.dat' file_string_6 = trim(file_string_6)//' '//trim(dir(n))//'/'//'MARK_INFLECTOR_DS_phase_space.dat_trimmed' file_string_7 = trim(file_string_7)//' '//trim(dir(n))//'/'//'DistStartInjLine_phase_space.dat_trimmed' endif 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_at(1)%muon_save, k, slope, offset) open(lun6,file= 'phase_vs_time.dat') write(lun6,'(3a12)') 'time','phase','number' do i=1,ix write(lun6,'(3es12.4)')time_sum(i)/number(i), phase_sum(i)/number(i) end do if(ix <=2)then 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') endif if(ix >= 2)then 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') print '(a)','cat '//trim(file_string_4)//' > all_Energy_spin_time.dat' call execute_command_line('cat '//trim(file_string_4)//' > all_Energy_spin_time.dat') print '(a)','cat '//trim(file_string_5)//' > all_omega_a.dat' call execute_command_line('cat '//trim(file_string_5)//' > all_omega_a.dat') print '(a)','cat '//trim(file_string_6)//' > all_MARK_INFLECTOR_DS_phase_space.dat_trimmed' call execute_command_line('cat '//trim(file_string_6)//' > all_MARK_INFLECTOR_DS_phase_space.dat_trimmed') print '(a)','cat '//trim(file_string_7)//' > all_DistStartInjLine_phase_space.dat_trimmed' call execute_command_line('cat '//trim(file_string_7)//' > all_DistStartInjLine_phase_space.dat_trimmed') endif end program energy_spin_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 real(rp) sxp(3), pvec(3) real(rp) sdotp real(rp) svec(3) 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 .or. index(long_string,'#')/=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 svec(1) = x elseif( i == 13)then svec(2) = x elseif( i == 14)then svec(3) = x elseif( i == 15)then acos_sp = x if(index(dir_file,'INFLECTOR') /=0 .or. index(dir_file,'DistStart')/=0 .or. index(dir_file,'BACKLEG')/=0)then pvec(1:2)=(/ang, ang_y/) pvec(3) = sqrt((1+pz)**2-pvec(1)**2-pvec(2)**2) spin= svec if(dot_product(spin,spin)>0)then call RightOrLeft(pvec,spin,sxp) sdotp= dot_product(svec,pvec) acos_sp = atan2(sxp(2),sdotp) if(acos_sp < 0)acos_sp = acos_sp + twopi endif endif 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 = svec 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 !example namelist input file ! &energy_spin_vs_time_input ! nargs = 22 ! dir_min = '00000000_000000' ! dir_max = '99999999_999999' ! phase_space_fille(1) = 'DistStartInjLine_phase_space.dat' ! phase_space_file(2) = 'MARK_INFLECTOR_DS_phase_space.dat ! phase_space_file(3) ='phase_space_at_01098_ns.dat' ! phase_space_file(4) ='phase_space_at_02216_ns.dat' ! phase_space_file(5) ='phase_space_at_03334_ns.dat' ! phase_space_file(6) ='phase_space_at_04452_ns.dat' ! phase_space_file(7) ='phase_space_at_05570_ns.dat' ! phase_space_file(8) ='phase_space_at_06688_ns.dat' ! phase_space_file(9) ='phase_space_at_07806_ns.dat' ! phase_space_file(10) ='phase_space_at_08924_ns.dat' ! phase_space_file(11) ='phase_space_at_10042_ns.dat' ! phase_space_file(12) ='phase_space_at_11160_ns.dat' ! phase_space_file(13) ='phase_space_at_12278_ns.dat' ! phase_space_file(14) ='phase_space_at_13396_ns.dat' ! phase_space_file(15) ='phase_space_at_14514_ns.dat' ! phase_space_file(16) ='phase_space_at_15632_ns.dat' ! phase_space_file(17) ='phase_space_at_16750_ns.dat' ! phase_space_file(18) ='phase_space_at_17868_ns.dat' ! phase_space_file(19) ='phase_space_at_18986_ns.dat' ! phase_space_file(20) ='phase_space_at_20104_ns.dat' ! phase_space_file(21) ='phase_space_at_21222_ns.dat' ! phase_space_file(22) ='phase_space_at_22340_ns.dat' !/