!+ ! Subroutine track1_custom (orbit, ele, param, track, err_flag, finished) ! ! For scattering in injection line "markers". And checking special apertures like inflector "D" ! And for recording electrons that intercept TRACKER planes (defined as markers) ! ! Note: This routine is not to be confused with track1_custom2. ! See the Bmad manual for more details. ! ! General rule: Your code may NOT modify any argument that is not listed as ! an output agument below." ! ! Modules Needed: ! use bmad ! ! Input: ! orbit -- Coord_struct: Starting position. ! ele -- Ele_struct: Element. ! param -- lat_param_struct: Lattice parameters. ! ! Output: ! orbit -- coord_struct: End position. ! track -- track_struct, optional: Structure holding the track information if the ! tracking method does tracking step-by-step. ! err_flg -- Logical: Set true if there is an error. False otherwise. !- subroutine track1_custom (orbit, ele, param, err_flag, finished, track) use bmad !use bmad_struct, only:ele_struct, coord_struct, lat_param_struct, track_struct !use bmad_interface use materials_mod use parameters_bmad !includes eloss, us_scatter, ds_scatter, tilt implicit none type (coord_struct) :: orbit, start_orb type (coord_struct) :: end_orb type (coord_struct) temp_orb, delta,temp_orb2 type (ele_struct) :: ele type (lat_param_struct) :: param type (track_struct), optional :: track logical hit/.false./,hit_y/.false./ logical err_flag logical energy_loss logical withinInflectorAperture logical first/.true./ !logical first1/.true./,first2/.true./,first3/.true./ !logical fbfirstx1/.true./,fbfirstx2/.true./,first_plane_detector/.true./ logical :: first_calo = .true., first_harp = .true., first_ibms=.true. logical finished !logical first_muonx1/.true./,first_muony1/.true./,first_muonx2/.true./,first_muony2/.true./ real(rp) limit_in/6.905/,limit_out/7.010/ !radial limits for the trackers real(rp), save :: t_first,t_first2,t_first3,t_first4 !first time a muon hit a fiber monitor real(rp) radius/7.112/ !radius of the ring real(rp) apRadius/0.045/ !radius of apeture real(rp) fbdia/0.0005/,spacing/0.013/ !diameter and length of intervals of the fibers real(rp),dimension(1:7):: xfiberPosition !positions of fibers measuring x real(rp)::tcycle = 2*Pi*rMagic/(0.9994*speedoflight) !time for muon to travel for one turn real(rp) thickness real(rp), save :: ibms_vec(1:3,1:3) real(rp), save:: dl real(rp) dx,dy character(32) :: r_name = 'track1_custom' character(26) :: filename character(3) :: turn integer lun,k,j integer nturn/128/ integer, save:: lun51,lun52 integer number_turns_dum/1/ integer caloNum ! calorimeter number (1 to 24) integer harpNum ! harp monitor number (1 to 4) integer muonNum integer fiber integer ibmsNum real(rp) :: x_hit, y_hit ! coords where particle hits calorimeter real(rp), parameter :: calo_edge_x1 = -0.31481 ! coords to define edges of calo front face real(rp), parameter :: calo_edge_x2 = -0.08981 real(rp), parameter :: calo_edge_y1 = -0.075 real(rp), parameter :: calo_edge_y2 = 0.075 ! !call out_io (s_fatal$, r_name, 'THIS DUMMY ROUTINE SHOULD NOT HAVE BEEN CALLED IN THE FIRST PLACE.') err_flag = .false. ! Remember to also set end_orb%t withinInflectorAperture = .true. energy_loss = eloss start_orb = orbit end_orb = start_orb temp_orb = start_orb if(first)then lun51=lunget() open(unit=lun51,file = trim(directory)//'/'//'UpstreamInflectorScattering.dat') lun52=lunget() open(unit=lun52,file = trim(directory)//'/'//'DownstreamInflectorScattering.dat') write(lun51,'(a)')'Change in coordinate vector due to upstream inflector scattering' write(lun51,'(6a12)')'x','xp','y','yp','z','zp' write(lun52,'(a)')'Change in coordinate vector due to downstream inflector scattering' write(lun52,'(6a12)')'x','xp','y','yp','z','zp' lun=lunget() open (unit=lun, file = trim(directory)//'/'//"log.dat",access='append') write(lun, '(a)')'TRACK1_CUSTOM: Scattering in inflector ends' write(lun, '(1x,a15,l)')' energy_loss = ', energy_loss write(lun,'(1x,a16)')ele%name first=.false. close(unit=lun) endif !print '(a,1x,2(a,2es12.4),2es12.4)',ele%name,'start_orb = ',start_orb%vec(1:2), ' end_orb = ', end_orb%vec(1:2), start_orb%s, ele%s if(ele%name == 'MARK_CRYO_US' .and. us_scatter)then call scatter( Al, 0.001_rp, temp_orb) ! end_orb = temp_orb end_orb%s = ele%s return endif if(ele%name == 'MARK_INFLECTOR_US')then withinInflectorAperture = .false. ! print '(a2,1x,6es12.4)','bs',temp_orb%vec(1:6) if(us_scatter)call inflector_scatter1(temp_orb, .true., .false., energy_loss) delta%vec = temp_orb%vec - start_orb%vec write(lun51,'(6es12.4)')delta%vec(1:6) ! print '(a2,1x,12es12.4,a,es12.4)','us',temp_orb%vec(1:6),delta%vec(1:6), ' inflector_width= ', inflector_width if(inflector_width <= 0.00901)then call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.) else call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.) endif if(.not. withinInflectorAperture)temp_orb%state=lost$ end_orb = temp_orb end_orb%s = ele%s orbit = end_orb return endif if(ele%name == 'MARK_INFLECTOR_DS')then withinInflectorAperture = .false. if(ds_scatter)call inflector_scatter1(temp_orb, .false., .true., energy_loss) delta%vec = temp_orb%vec - start_orb%vec write(lun52,'(6es12.4)')delta%vec(1:6) ! print '(a2,1x,6es12.4)','ds',delta%vec(1:6) if(inflector_width <= 0.00901)then call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.) else call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.) endif if(.not. withinInflectorAperture)temp_orb%state = lost$ temp_orb%vec(1:6) = temp_orb%vec(1:6) + delta%vec(1:6)*extra_inf_scatter_factor end_orb = temp_orb end_orb%s = ele%s orbit = end_orb return endif !========================================== ! fiber monitor section ! added by dlr, September 2016 ! revised for multiple scattering January 2017 !========================================== if (index(ele%name, 'FIBER')/=0) then ! check coords to see if particle hit harp front face x_hit = start_orb%vec(1) y_hit = start_orb%vec(3) if (first_harp) then ! create file and write header info lun = lunget() open(unit=lun, file=trim(directory)//'/'//"harp_plane_hits.dat", status="new", action="write") write(unit=lun, fmt='(a6,8a12)') " harp", & "t ","s ", & "x ","px ", & "y ","py ", & "z ","pz " close(unit=lun) open(unit=lun, file=trim(directory)//'/'//"fiber_hits.dat", status="new", action="write") write(unit=lun, fmt='(a6,1x,a6,9a12)') " harp","fiber", & "t ","s ", & "x ","px ", & "y ","py ", & "z ","pz ",'thickness' close(unit=lun) first_harp = .false. end if ! get harp num (integer) from element name read (ele%name(6:6),'(i1)') harpNum call which_fiber(harpNum, start_orb, fiber, thickness) ! write information about this hit to the file lun = lunget() open(unit=lun, file=trim(directory)//'/'//"harp_plane_hits.dat", status="old", access="append", action="write") write(unit=lun, fmt='(i5,8es18.8)') harpNum, & start_orb%t, start_orb%s, & start_orb%vec(1), start_orb%vec(2), & start_orb%vec(3), start_orb%vec(4), & start_orb%vec(5), start_orb%vec(6) close(unit=lun) if(fiber /= 0)then open(unit=lun, file=trim(directory)//'/'//"fiber_hits.dat", status="old", access="append", action="write") write(unit=lun, fmt='(i5,1x,i5,9es18.8)') harpNum, fiber, & start_orb%t, start_orb%s, & start_orb%vec(1), start_orb%vec(2), & start_orb%vec(3), start_orb%vec(4), & start_orb%vec(5), start_orb%vec(6), thickness close(unit=lun) end_orb = start_orb orbit = end_orb if(FiberScatter) call scatter(SciFiber, thickness, end_orb) ! call EnergyLoss(SciFiber,thickness,end_orb) if(FiberEnergyLoss) call fiber_energy_loss(thickness, end_orb) if(FiberEnergyLoss .or. FiberScatter) write(111,'(7es12.4)')thickness,start_orb%vec - end_orb%vec endif ! fiber hit return end if ! harp plane !======================== ! end fiber monitor section !======================== !========================================== ! calorimeter section ! added by Robin Bjorkquist, September 2015 !========================================== number_turns_dum=1 !ix_end_flag = 1 if(index(ele%name,'QUADM')/=0 .and. ix_end_flag /=0) call compute_moments_vs_time_calo(number_turns_dum, ix_end_flag, start_orb, ele%name) !if(index(ele%name,'QUADM')/=0 .and. ix_end_flag /=0) call compute_moments_vs_time_calo_sum(number_turns_dum, ix_end_flag, start_orb, ele%name) if (start_orb%species == positron$) then if (index(ele%name, 'CALO')/=0) then ! check coords to see if particle hit calorimeter front face x_hit = start_orb%vec(1) y_hit = start_orb%vec(3) if (calo_edge_x1 < x_hit .and. x_hit < calo_edge_x2 .and. calo_edge_y1 < y_hit .and. y_hit < calo_edge_y2) then if (first_calo) then ! create file and write header info lun = lunget() open(unit=lun, file=trim(directory)//'/'//"calo_hits.dat", status="new", action="write") write(unit=lun, fmt='(a10,a5,8a12)') " muon_num"," calo", & "t ","s ", & "x ","px ", & "y ","py ", & "z ","pz " close(unit=lun) first_calo = .false. end if ! get calo num (integer) from element name read (ele%name(5:6),'(i2)') caloNum ! get muon num, which is stored in the ix_user field of the positron's coord_struct muonNum = start_orb%ix_user ! write information about this hit to the file lun = lunget() open(unit=lun, file=trim(directory)//'/'//"calo_hits.dat", status="old", access="append", action="write") write(unit=lun, fmt='(i10,i5,8es12.4)') muonNum, & caloNum, & start_orb%t, start_orb%s, & start_orb%vec(1), start_orb%vec(2), & start_orb%vec(3), start_orb%vec(4), & start_orb%vec(5), start_orb%vec(6) close(unit=lun) ! kill positron track if(start_orb%species == positron$)temp_orb%state = lost$ end if ! particle hit front face end if ! CALO end_orb = temp_orb end_orb%s = ele%s orbit = end_orb return end if ! positron !======================== ! end calorimeter section !======================== !===================== !IBMS !===================== if(index(ele%name, 'IBMS')/= 0 .and. ibms)then if(first_ibms)then lun=lunget() open(unit = lun, file = trim(directory)//'/'//"IBMS_hits.dat",status="new",action="write") write(lun,'(a11,6a12)')' detector ',' x ','xp',' y ','yp',' s ',' t ' close(unit=lun) open(unit=lun, file=trim(directory)//'/'//"ibms_fiber_hits.dat", status="new", action="write") write(unit=lun, fmt='(a6,1x,a6,9a12)') "ibms","fiber", & "t ","s ", & "x ","px ", & "y ","py ", & "z ","pz ",'thickness' close(unit=lun) first_harp = .false. ! !MacCoy Semptember 27, 2022 ibms_vec(1,1)=-0.045 ibms_vec(2,1)= 0.0109 ! ibms_vec(3,1)= 0.0681 !at inflector exit !7180.1 - 7112 ibms_vec(3,1)= 0.0 !in ring !7180.1 - 7112 ibms_vec(1,2)=0. ibms_vec(2,2)= 0.0021 ! ibms_vec(3,2)=0.0021 !at inflector exit ibms_vec(3,2)=0.0 !in ring ibms_vec(1,3)=-4.397 !this is 0.097m upstream of start dl = ibms_vec(1,3) + 4.3 ibms_vec(2,3)= -1.9505 ibms_vec(3,3)= 1.016 first_ibms=.false. endif read (ele%name(5:5),'(i1)') ibmsNum dx=0. dy=0. if(ibmsNum == 1)then dx = dl * start_orb%vec(2) dy = dl*start_orb%vec(4) endif end_orb=start_orb orbit = end_orb !end_orb%vec(1) = start_orb%vec(1)-ibms_vec(ibmsNum,1)+dx !end_orb%vec(3) = start_orb%vec(3)-ibms_vec(ibmsNum,2)+dy call which_ibms_fiber(ibmsNum, end_orb, fiber, thickness) end_orb = start_orb orbit = end_orb lun=lunget() open(unit = lun, file = trim(directory)//'/'//"IBMS_hits.dat",position="append", access="sequential") write(lun,'(i11,1x,6es12.4)')ibmsNum, start_orb%vec(1)-ibms_vec(ibmsNum,1),start_orb%vec(2), start_orb%vec(3)-ibms_vec(ibmsNum,2),start_orb%vec(4), & ibms_vec(ibmsNum,3), start_orb%t ! write(lun,'(i11,1x,6es12.4)')ibmsNum, start_orb%vec(1),start_orb%vec(2), start_orb%vec(3),start_orb%vec(4), & ! ibms_vec(ibmsNum,3), start_orb%t close(unit=lun) if(fiber /= 0)then open(unit=lun, file=trim(directory)//'/'//"ibms_fiber_hits.dat", status="old", access="append", action="write") write(unit=lun, fmt='(i5,1x,i5,9es18.8)') ibmsNum, fiber, & start_orb%t, start_orb%s, & start_orb%vec(1), start_orb%vec(2), & start_orb%vec(3), start_orb%vec(4), & start_orb%vec(5), start_orb%vec(6), thickness close(unit=lun) endif endif !if(end_orb%state == lost$)print '(a12,1x,a16,7es12.4)','Lost in',ele%name,end_orb%s,end_orb%vec(1:6) end subroutine