!+ ! Subroutine track1_custom (start_orb, ele, param, end_orb, 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: ! start_orb -- Coord_struct: Starting position. ! ele -- Ele_struct: Element. ! param -- lat_param_struct: Lattice parameters. ! ! Output: ! end_orb -- 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 (start_orb, ele, param, end_orb, err_flag, finished, track) use bmad_struct, only:ele_struct, coord_struct, lat_param_struct, track_struct use bmad_interface, dummy => track1_custom use materials_mod use parameters_bmad !includes eloss, us_scatter, ds_scatter, tilt implicit none type (coord_struct) :: 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. 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 character(32) :: r_name = 'track1_custom' character(26) :: filename character(3) :: turn integer lun,k,j integer nturn/128/ integer, save:: lun51,lun52, lun53 integer caloNum ! calorimeter number (1 to 24) integer harpNum ! harp monitor number (1 to 4) integer muonNum integer fiber 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 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') lun53=lunget() open(unit=lun53,file = trim(directory)//'/'//'CryoBeginScattering.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' write(lun53,'(a)')'Change in coordinate vector due to downstream inflector scattering' write(lun53,'(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 offset_particle(ele, set$, temp_orb) call scatter( Al, 0.001_rp, temp_orb) call offset_particle(ele, unset$, temp_orb) delta%vec = temp_orb%vec - start_orb%vec write(lun53,'(6es12.4)')delta%vec(1:6) 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)then call offset_particle(ele, set$, temp_orb) call inflector_scatter1(temp_orb, .true., .false., energy_loss) call offset_particle(ele, unset$, temp_orb) delta%vec = temp_orb%vec - start_orb%vec write(lun51,'(6es12.4)')delta%vec(1:6) endif ! print '(a2,1x,12es12.4,a,es12.4)','us',temp_orb%vec(1:6),delta%vec(1:6), ' inflector_width= ', inflector_width ! this next conditional determines apertures if custom_tracking through the inflector is used. ! Set inflector_width=0 in input.dat to skip if aperture is set in lattice file if(inflector_width <= 0.00901 .and. inflector_width > 0)then call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.) elseif(inflector_width > 0.00901)then call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.) endif if(inflector_width > 0 .and. .not. withinInflectorAperture)temp_orb%state=lost$ end_orb = temp_orb end_orb%s = ele%s return endif if(ele%name == 'MARK_INFLECTOR_DS')then withinInflectorAperture = .false. if(ds_scatter)then call offset_particle(ele, set$, temp_orb) call inflector_scatter1(temp_orb, .false., .true., energy_loss) call offset_particle(ele, unset$, temp_orb) delta%vec = temp_orb%vec - start_orb%vec write(lun52,'(6es12.4)')delta%vec(1:6) endif ! print '(a2,1x,6es12.4)','ds',delta%vec(1:6) ! this next conditional determines apertures if custom_tracking through the inflector is used. ! Set inflector_width=0 in input.dat to skip if aperture is set in lattice file if(inflector_width <= 0.00901 .and. inflector_width > 0)then call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.) elseif(inflector_width > 0.00901)then call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.) endif if(inflector_width > 0 .and. .not. withinInflectorAperture)temp_orb%state = lost$ end_orb = temp_orb end_orb%s = ele%s 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 calorimeter 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 if(FiberScatter)then call scatter(SciFiber, thickness, end_orb) call EnergyLoss(SciFiber,thickness,end_orb) endif 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 !========================================== 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 return end if ! positron !======================== ! end calorimeter section !======================== !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