!........................................................................ !+ ! program : injtrack ! ! Description: track injected particles and calculate the inject efficiency ! ! ! Arguments : param.in ! ! Mod/Commons: ! ! Calls : ! ! Author : Suntao Wang ! ! Modified : 2015.11.09. !- !........................................................................ ! program injtrack use bmad use precision_def use quick_plot use random_mod use bmadz_interface use bsim_cesr_interface use sim_utils use cesr_lrbbi_mod use setup_group_mod use cesr_group_mod use cesr_db_mod use bookkeeper_mod ! use reverse_mod implicit none character*150 file_name, out_file, screen_file, lossname, und_file character*150 lat_file,line, text_stat, position_file, check_ele_name real(rp) init_orbit(6), init_sigma(6), stored_beam_sigma(6),temp(6) ! input initial orbits and its Gaussian sigma real(rp), allocatable :: inj_init_orbit(:,:) ! the initial injection orbits of all the particles real(rp), allocatable :: lost_info(:,:) ! the orbits at the lost position real(rp), allocatable :: lost_init(:,:) ! the initial injection orbits for all the lost particles real(rp), allocatable :: und_cords(:,:,:,:) ! the cordinates which enter and exit two KYMA undulators real(rp) harvest(6) ! random generator temp variable real(rp) frev ! revolution frequency real(rp) injeff, scv_eff(2), rad_mon(2) real(rp) emit_x, emit_y real(rp) betay_col43w, betay_col43e, betax_34w,betay_34w type (lat_struct), target :: ring, ring_rev ! This structure holds the lattice info type (ele_struct), pointer :: ele type (ele_pointer_struct), allocatable :: eles(:) type (coord_struct), allocatable :: orbit(:), orb(:) type (normal_modes_struct) modes integer n_arg, n_turns, n_particles, track_specie integer i, j, k, ix, ix_start, ix_end, ix_branch, track_state, ix_u ! indices integer, allocatable :: lost_ix(:,:) integer col43w_idx, col43e_idx, ccu_idx(2), injs34w_idx, injs34e_idx, check_ele_idx integer i_repeat, repeat_times,random_seed integer i_col43w, i_col43e, i_ccuw1, i_ccuw2 logical err, rand_init, out_screen, save_und_cord logical ok, rad_damp, use_stored_beam_size, lrbbi_off real(rp) target_tunes(3), tunes_hz(3), temp_tunes(3), delta_var, emit_ratio integer int_Qx, int_Qy, Qtune_steps, ix_var, var_steps integer ix_Qx, ix_Qy, n_x, n_y real(rp) dQx, dQy, dQz real(rp) Q_x0, Q_x1, dQ_x, Q_y0, Q_y1, dQ_y, Q_z character(8) :: date character(10) :: time character(5) :: zone integer,dimension(8) :: values n_particles=1 frev = 390.136 Qtune_steps = 10 rad_damp = .false. use_stored_beam_size = .false. lrbbi_off = .false. stored_beam_sigma(:)= 0 var_steps = 1 delta_var = 0 ix_var = 1 emit_ratio = 1.0 out_screen = .true. save_und_cord = .false. check_ele_name = 'IP_L0' ! input list parameters namelist / input / lat_file, n_particles, n_turns, init_orbit, init_sigma, track_specie, & random_seed, repeat_times,rand_init, Qtune_steps, rad_damp, use_stored_beam_size, & var_steps, ix_var, delta_var, emit_ratio, lrbbi_off, Q_x0, Q_x1, dQ_x, Q_y0, Q_y1, dQ_y, Q_z, & out_screen, save_und_cord, check_ele_name ! read in the paramter file n_arg = cesr_iargc() if (n_arg > 1) then print *, 'Usage: injtrack ' print *, 'Default: = param.in' stop endif if (n_arg == 1) call cesr_getarg(1, file_name) if (n_arg == 0) then print '(a, $)', ' Input command file : ' read (*, '(a)') file_name endif if (file_name == '') file_name = 'param.in' print *, 'Opening: ', trim(file_name) open (unit= 1, file = file_name, status = 'old') read(1, nml = input) close (unit = 1) ! out put file name ix = index(file_name,'.') out_file = file_name(1:ix)//'out' screen_file = file_name(1:ix)//'txt' call date_and_time(date,time,zone,values) allocate(inj_init_orbit(n_particles,6)) allocate(lost_init(n_particles,6)) allocate(lost_info(n_particles,6)) allocate(lost_ix(n_particles,4)) ! random_seed = 0 call ran_seed_put (random_seed) !use the same random pattern ! Read in a lattice and calculate the twiss parameters. bmad_com%auto_bookkeeper = .false. ! bmad_com%rel_tol_adaptive_tracking = 1e-9 ! set Runga tracking tolerance lower ! bmad_com%abs_tol_adaptive_tracking = 1e-11 call bmad_parser (lat_file, ring) ! Read in a lattice. ! CESRV output a lattice file in which the quads are used for qtuneing ! Their slave_status contributes is contro_slave$ (Qtuning), have to remove Qtuning groups in order to be used in choose_quads() ! Otherwise those quads won't be chosen in the routine. do i = 1, ring%n_ele_max if (ring%ele(i)%key == group$ .and. ring%ele(i)%name(1:12) == 'CSR_QTUNEING') then ring%ele(i)%ix_ele = -1 endif enddo call remove_eles_from_lat(ring) if (rad_damp) then call set_on_off(rfcavity$, ring, on$) bmad_com%radiation_damping_on = .true. !turn on radiation damping bmad_com%radiation_fluctuations_on = .true. endif if (lrbbi_off) then call set_on_off(beambeam$, ring, off$) endif if ( track_specie .eq. 1 ) then ring%param%particle = positron$ !set the ref charge to be positron, default tracking charge should be the same. else ring%param%particle = positron$ !set the ref charge to be positron !call lat_reverse (ring, ring_rev) !lat_reverse keep positron$ B field but set the tracking specie to be electron call reverse_lat (ring, ring_rev, track_antiparticle = .true.) ring = ring_rev endif call closed_orbit_calc (ring, orb, 6) call lat_make_mat6 (ring, -1, orb) call twiss_at_start (ring) call twiss_propagate_all (ring) call calc_z_tune(ring) ! find the indices of the scrapers and the undulators do i=1,ring%n_ele_track if (ring%ele(i)%name .eq. 'COL43W') then col43w_idx = i elseif (ring%ele(i)%name .eq. 'COL43E') then col43e_idx = i ! elseif (ring%ele(i)%name .eq. 'CCU_W1') then elseif (ring%ele(i)%name .eq. 'CCU_FLINE_STRAIGHT') then ccu_idx(1)=i ! elseif (ring%ele(i)%name .eq. 'CCU_W2') then elseif (ring%ele(i)%name .eq. 'CCU_GLINE_STRAIGHT') then ccu_idx(2)=i elseif (ring%ele(i)%name .eq. 'INJS34W') then injs34w_idx = i elseif (ring%ele(i)%name .eq. 'INJS34E') then injs34e_idx = i elseif (ring%ele(i)%name .eq. check_ele_name) then check_ele_idx = i endif enddo if (save_und_cord) then allocate(und_cords(n_particles, n_turns, 6, ring%n_ele_track)) ! save the closed orbit data open(unit=100, file='closed_orb.dat', status='replace') do i=1,ring%n_ele_track write(100, '(i8, 7es18.5e5)') i,ring%ele(i)%s, orb(i)%vec(:) enddo close(100) endif open(unit=101, file=out_file, status='replace') write(101, '(18a10)') '! var_step', 'init_orb', 'i_repeat', 'inj_eff', 'Qx','Qy','Qz','ccu_w1','ccu_w2','col43w','col43e','Qx_real','Qy_real','Qz_real','by_43w','by_43e','bx_34w','by_34w' open(unit=200, file=screen_file, status='replace') write(200, '(a15, a,2x,a,2x,a)') 'Starting time: ', date, time, zone write (200, '(a)') '==========================================================================================' write (200, '(a6, 6a12,a7,a12)') ' ', 'CCU_W1', 'CCU_W2', 'Septum34W', 'col43W', 'col43E', 'Septum34E', ' ', check_ele_name write (200, '(a6, 7I12)') 'Index:', ccu_idx(1), ccu_idx(2), injs34w_idx, col43w_idx, col43e_idx, injs34e_idx, check_ele_idx write (200, '(a)') '==========================================================================================' write (200,'(a16,i5,a21,i5)') 'No. of elements=',ring%n_ele_max,'No.of tracking lat=',ring%n_ele_track if ( track_specie .eq. 1) then ! positron write (200,'(a8,i5,a15,i5,a7,a7,i5,a6)') 'Track of',n_particles,' positrons for',n_turns,'turns,','repeat',repeat_times,'times' else ! electron write (200,'(a8,i5,a15,i5,a7,a7,i5,a6)') 'Track of',n_particles,' electrons for',n_turns,'turns,','repeat',repeat_times,'times' endif ! Scan tunes for injection effiency calculation if(dQ_x .gt. 0) then n_x = anint(abs((Q_x1 - Q_x0) / dQ_x)) else n_x = 0 endif if(dQ_y .gt. 0) then n_y = anint(abs((Q_y1 - Q_y0) / dQ_y)) else n_y = 0 endif int_Qx = int(ring%ele(ring%n_ele_track)%a%phi / twopi) int_Qy = int(ring%ele(ring%n_ele_track)%b%phi / twopi) write (200, '(a28,f12.6,a9,f12.6,a9,f12.6)') 'Before changing Qtune: Qx = ',ring%a%tune/twopi*frev, & ' Qy = ',ring%b%tune/twopi*frev,' Qz = ',ring%z%tune/twopi*frev do ix_Qy = 0, n_y tunes_hz(2) = Q_y0 + ix_Qy*dQ_y temp_tunes(2)=tunes_hz(2) do ix_Qx = 0, n_x tunes_hz(1) = Q_x0 + ix_Qx*dQ_x tunes_hz(3)=Q_z ! set the output und_cords file name for each tune if (save_und_cord) then write(und_file,'(A,I3,A,I3,A)') 'und_',floor(tunes_hz(1)),'_',floor(tunes_hz(2)),'.dat' open(unit=300, file=und_file, ACTION="write", STATUS="replace", FORM="unformatted") do k =1, n_particles do i = 1, n_turns do ix_u = 1, ring%n_ele_track und_cords(k, i, :, ix_u) = -1 enddo enddo enddo endif ! Set tunes write (200,'(a28,f12.6,a9,f12.6,a9,f12.6)') 'Target tunes: Qx = ',tunes_hz(1), & ' Qy = ',tunes_hz(2),' Qz = ', tunes_hz(3) temp_tunes(1)=tunes_hz(1) temp_tunes(3)=tunes_hz(3) target_tunes(1) = tunes_hz(1)/frev target_tunes(2) = tunes_hz(2)/frev target_tunes(3) = tunes_hz(3)/frev dQx = ( target_tunes(1) - ring%a%tune/twopi ) / Qtune_steps dQy = ( target_tunes(2) - ring%b%tune/twopi ) / Qtune_steps dQz = ( target_tunes(3) - ring%z%tune/twopi ) / Qtune_steps target_tunes(1) = ring%a%tune/twopi + int_Qx target_tunes(2) = ring%b%tune/twopi + int_Qy target_tunes(3) = ring%z%tune/twopi global_com%exit_on_error=.false. do i=1,Qtune_steps target_tunes(1) = target_tunes(1) + dQx target_tunes(2) = target_tunes(2) + dQy target_tunes(3) = target_tunes(3) + dQz call set_tune3(ring, target_tunes, orb, ok, .false.) enddo global_com%exit_on_error=.true. call lattice_bookkeeper(ring) call closed_orbit_calc (ring, orb, 6) call lat_make_mat6 (ring, -1, orb) call twiss_at_start (ring) call twiss_propagate_all (ring) call radiation_integrals (ring, orb, modes) tunes_hz(1)=ring%a%tune/twopi*frev tunes_hz(2)=ring%b%tune/twopi*frev tunes_hz(3)=ring%z%tune/twopi*frev write (200,'(a28,f12.6,a9,f12.6,a9,f12.6)')'After Qtune: Qx = ',tunes_hz(1),' Qy = ',tunes_hz(2),' Qz = ',tunes_hz(3) betay_col43w=ring%ele(col43w_idx)%b%beta betay_col43e=ring%ele(col43e_idx)%b%beta betax_34w=ring%ele(injs34w_idx)%a%beta betay_34w=ring%ele(injs34w_idx)%b%beta if ( track_specie .eq. 1) then ! positron ix_start = injs34w_idx write (200,'(a39)') 'The stored beam orbit at Septum 34W is ' write (200,'(6a14)'), "x","x'","y","y'","z","z'" write (200,'(6es14.5)') orb(injs34w_idx)%vec(:) else ! electron ix_start = injs34e_idx write (200,'(a39)') 'The stored beam orbit at Septum 34E is ' write (200,'(6a14)') "x","x'","y","y'","z","z'" write (200,'(6es14.5)') orb(injs34e_idx)%vec(:) endif ix_end = ix_start ! calcuate the stored beam dimension at the injection point emit_x = modes%a%emittance emit_y = emit_x * emit_ratio ! write (200,'(a14,f10.2,a5, a14,f10.2,a2)') 'Horizon emit:', emit_x*1E9,'nm ','Vertical emit:', emit_y*1E9,'nm' if (use_stored_beam_size) then stored_beam_sigma(1) = sqrt(emit_x * ring%ele(ix_start)%a%beta + (ring%ele(ix_start)%a%eta * modes%sigE_E)**2) stored_beam_sigma(3) = sqrt(emit_y * ring%ele(ix_start)%b%beta + (ring%ele(ix_start)%b%eta * modes%sigE_E)**2) stored_beam_sigma(5) = modes%sig_z init_sigma(:)= stored_beam_sigma(:) endif write (200,'(a26, 6es10.2)') 'The inital injection orbit: ',init_orbit(:) write (200,'(a26, 6es10.2)') 'The inital injection sigma: ',init_sigma(:) write (200,'(7a10,a12)') 'Steps', 'Repeat', 'Particle','Lost At', 'Turns','Ele idx','Plane', 'Ele name' if(out_screen) then write (*,'(7a10,a12)') 'Steps', 'Variable', 'Particle','Lost At', 'Turns','Ele idx','Plane', 'Ele name' endif temp(:)=init_orbit(:) do j = 1, var_steps init_orbit(ix_var)=temp(ix_var) + delta_var*(j- floor(var_steps/2.0_rp)-1) do i_repeat= 1,repeat_times do k = 1, n_particles lost_ix(k,:)=-1 ! initialize the orbit data enddo call reallocate_coord (orbit, ring%n_ele_max) orbit(ix_start)%vec = init_orbit ! init orbit at the injection point do k = 1, n_particles call ran_gauss(harvest) ! generate gaussian distributed random numbers with unit sigma if(rand_init) then inj_init_orbit(k,1)=init_orbit(1)+init_sigma(1)*harvest(1) ! record all the initial coordinates inj_init_orbit(k,2)=init_orbit(2)+init_sigma(2)*harvest(2) inj_init_orbit(k,3)=init_orbit(3)+init_sigma(3)*harvest(3) inj_init_orbit(k,4)=init_orbit(4)+init_sigma(4)*harvest(4) inj_init_orbit(k,5)=init_orbit(5)+init_sigma(5)*harvest(5) inj_init_orbit(k,6)=init_orbit(6)+init_sigma(6)*harvest(6) endif orbit(ix_start)%vec = inj_init_orbit(k,:) ! transfer the initial orbits to the starting element do i =1, n_turns if (track_specie .eq. 1) then call track_many (ring, orbit, ix_start, ix_end, 1, 0, track_state) ! positron else ! ring%param%rel_tracking_charge = -1 ! call track_reverse (ring, orbit, ix_start, ix_end, 0, track_state) ! electron call track_many (ring, orbit, ix_start, ix_end, 1, 0, track_state) ! electron endif if (save_und_cord) then do ix_u = 1, ring%n_ele_track und_cords(k, i, :, ix_u) = orbit(ix_u)%vec(:) enddo endif if (track_state /= moving_forward$) then ! The particle is lost lost_ix(k,1)=track_state ! record at which element the particle is lost lost_ix(k,2)=i ! record at which turn the particle is lost if (track_state .eq. ccu_idx(1)) then lost_ix(k,3) = 1 elseif (track_state .eq. ccu_idx(2)) then lost_ix(k,3) = 2 elseif (track_state .eq. col43w_idx) then lost_ix(k,3) = 3 elseif (track_state .eq. col43e_idx) then lost_ix(k,3) = 4 endif lost_ix(k,4) = orbit(lost_ix(k,1))%state ! record the particle lost at which plane: (3:-x; 4:+x; 5:-y; 6:+y; 7:z) lost_info(k,:)=orbit(lost_ix(k,1))%vec(:) ! record the lost particle info at the lost position exit endif orbit(ix_start) = orbit(ix_end) lost_ix(k,1) = -1 ! The particle is survived. lost_ix(k,2) = n_turns+1 enddo ! i= n_turns if (lost_ix(k,1) .eq. -1) then lossname='noloss' else lossname=ring%ele(lost_ix(k,1))%name endif if (out_screen) then write (*,'(I10, es10.2, 5I10, 4x, 1a10)') j, init_orbit(ix_var), k, lost_ix(k,1), lost_ix(k,2), lost_ix(k,3), lost_ix(k,4), lossname endif write (200,'(7I10, 4x, 1a10)') j, i_repeat, k, lost_ix(k,1), lost_ix(k,2), lost_ix(k,3), lost_ix(k,4), lossname enddo ! k=n_particles ! count the particles and do the statistics analysis ix=0 ! counts for survived particles i_col43w=0 i_col43e=0 i_ccuw1=0 i_ccuw2=0 do k=1, n_particles if ( lost_ix(k,1) .eq. -1) then ix=ix+1 elseif (lost_ix(k,1) .eq. col43w_idx) then i_col43w=i_col43w+1 elseif (lost_ix(k,1) .eq. col43e_idx) then i_col43e=i_col43e+1 elseif ((ring%ele(lost_ix(k,1))%name(1:4) .eq. 'CCU_') .or. (lost_ix(k,1) .eq. ccu_idx(1)) .or. (ring%ele(lost_ix(k,1))%name(1:4) .eq. 'und_') ) then i_ccuw1=i_ccuw1+1 !count for all particles lost at CCUs elseif ((ring%ele(lost_ix(k,1))%name(1:5) .eq. 'CCU_W') .or. (lost_ix(k,1) .eq. ccu_idx(2)) .or. (ring%ele(lost_ix(k,1))%name(1:8) .eq. 'und_08wn')) then i_ccuw2=i_ccuw2+1 ! count for all particles lost endif enddo injeff= ix*100.0_rp/n_particles scv_eff(1)=i_col43w*100.0_rp/n_particles scv_eff(2)=i_col43e*100.0_rp/n_particles rad_mon(1)=i_ccuw1*100.0_rp/n_particles rad_mon(2)=i_ccuw2*100.0_rp/n_particles write (200,'(a15, f6.2, a1,a15, f6.2, a1,a15, f6.2, a1)') 'Injection eff: ', injeff, '%', & 'Scrap eff',scv_eff(1)+scv_eff(2),'%', 'Rad_mon',rad_mon(1),'%' if(.not. rand_init) exit write(101, '(i10, es10.2, i10, 15f10.2)') j, init_orbit(ix_var), i_repeat, injeff, temp_tunes(1), temp_tunes(2), temp_tunes(3), & rad_mon(1), rad_mon(2), scv_eff(1), scv_eff(2), tunes_hz(1), tunes_hz(2), tunes_hz(3), betay_col43w, betay_col43e,betax_34w,betay_34w if (save_und_cord) then write(300) und_cords close(unit=300) endif enddo !i_repeat = repeat_times enddo ! var_step enddo !ix_Qx enddo !ix_Qy if (save_und_cord) then deallocate(und_cords) endif close(101) call date_and_time(date,time,zone,values) write(200, '(a16, a,2x,a,2x,a)') 'Finishing time: ', date, time, zone close(200) end program