MODULE materials_mod USE precision_def ! use bmad IMPLICIT NONE TYPE material CHARACTER(LEN=20) :: name CHARACTER(LEN=20) :: symbol REAL(RP) :: Z ! nuclear charge [e] REAL(RP) :: A ! atomic mass [g/mol] REAL(RP) :: density ! density [g/cm3] REAL(RP) :: radl ! radiation length [cm] REAL(RP) :: nucintl ! nuclear interaction length [cm] REAL(RP) :: imean ! mean excitation energy [eV] END TYPE material ! Vacuum TYPE (material), TARGET :: Vacuum ! Elements (arranged by increasing nuclear charge) TYPE (material), TARGET :: Be ! Z=04 TYPE (material), TARGET :: Al ! Z=13 TYPE (material), TARGET :: Si ! z=14 TYPE (material), TARGET :: Ca ! z=6 TYPE (material), TARGET :: SciFiber ! z=6.5 TYPE (material), TARGET :: Ti ! z=22 TYPE (material), TARGET :: Cu ! z=29 TYPE (material), TARGET :: Nb ! z=41 ! Compounds (arranged by alphabetical order) TYPE (material) :: Kapton TYPE (material) :: InflectorCoilE821 TYPE (material) :: InflectorEndE821 TYPE (material) :: NbTi ! Convenient list of all available materials TYPE (material), TARGET, DIMENSION(13) :: materials INTEGER :: num_materials ! counter/index ! Pointers for materials in tracking TYPE (material), POINTER :: material_ptr ! general purpose TYPE (material), POINTER :: material_InflectorCryoWindowUS ! upstream TYPE (material), POINTER :: material_InflectorFlangeWindowUS ! | TYPE (material), POINTER :: material_InflectorCoilUS ! | TYPE (material), POINTER :: material_InflectorMandrelWindowUS ! | TYPE (material), POINTER :: material_InflectorMandrelWindowDS ! | TYPE (material), POINTER :: material_InflectorCoilDS ! | TYPE (material), POINTER :: material_InflectorFlangeWindowDS ! v TYPE (material), POINTER :: material_InflectorCryoWindowDS ! downstream TYPE (material), POINTER :: material_fiber logical quad_plate ! if true scatter when passing through quad plate - wall_hit_handler_custom CONTAINS SUBROUTINE InitializeMaterials() IMPLICIT NONE ! Counter num_materials = 0 ! Vacuum (all we really care about is the name) Vacuum%name = 'Vacuum' Vacuum%symbol = 'Vacuum' num_materials = num_materials+1 materials(num_materials) = Vacuum ! Elements (arranged by increasing nuclear charge) Be%name = 'Beryllium' Be%symbol = 'Be' Be%Z = 4. Be%A = 9.01 ! g/mol Be%density = 1.848 ! g/cm3 Be%radL = 35.276 ! cm Be%nucIntL = 39.449 ! cm Be%Imean = 63.7 ! eV num_materials = num_materials+1 materials(num_materials) = Be Al%name = 'Aluminum' Al%symbol = 'Al' Al%Z = 13. Al%A = 26.98 ! g/mol Al%density = 2.699 ! g/cm3 Al%radL = 8.896 ! cm Al%nucIntL = 38.877 ! cm Al%Imean = 166.0 ! eV num_materials = num_materials+1 materials(num_materials) = Al Si%name = 'Silicon' Si%symbol = 'Si' Si%Z = 14. Si%A = 28.09 ! g/mol Si%density = 2.330 ! g/cm3 Si%radL = 9.366 ! cm Si%nucIntL = 45.635 ! cm Si%Imean = 173.0 ! eV num_materials = num_materials+1 materials(num_materials) = Si Ca%name = 'Carbon' Ca%symbol = 'Ca' Ca%Z = 6 Ca%A = 12 ! g/mol Ca%density = 2.330 ! g/cm3 Ca%radL = 9.366 ! cm Ca%nucIntL = 45.635 ! cm Ca%Imean = 173.0 ! eV num_materials = num_materials+1 materials(num_materials) = Ca SciFiber%name = 'Kuraray' SciFiber%symbol = 'SciFiber' SciFiber%Z = 6 SciFiber%A = 12 ! g/mol SciFiber%density = 1.05 ! g/cm3 SciFiber%radL = 42.4 ! cm SciFiber%nucIntL = 55. ! cm SciFiber%Imean = 64.7 ! eV num_materials = num_materials+1 materials(num_materials) = SciFiber Ti%name = 'Titanium' Ti%symbol = 'Ti' Ti%Z = 22. Ti%A = 47.87 ! g/mol Ti%density = 4.540 ! g/cm3 Ti%radL = 3.560 ! cm Ti%nucIntL = 27.971 ! cm Ti%Imean = 233.0 ! eV num_materials = num_materials+1 materials(num_materials) = Ti Cu%name = 'Copper' Cu%symbol = 'Cu' Cu%Z = 29. Cu%A = 63.55 ! g/mol Cu%density = 8.96 ! g/cm3 Cu%radL = 1.436 ! cm Cu%nucIntL = 15.576 ! cm Cu%Imean = 322.0 ! eV num_materials = num_materials+1 materials(num_materials) = Cu Nb%name = 'Niobium' Nb%symbol = 'Nb' Nb%Z = 41. Nb%A = 92.91 ! g/mol Nb%density = 8.57 ! g/cm3 Nb%radL = 1.158 ! cm Nb%nucIntL = 18.485 ! cm Nb%Imean = 417.0 ! eV num_materials = num_materials+1 materials(num_materials) = Nb ! Compounds (arranged by alphabetical order) ! The superconducting inflector coil wires are insulated with Kapton Kapton%name = 'Kapton' Kapton%symbol = 'Kapton' Kapton%Z = 5.02779 ! calculated from number fractions (stoichiometry) and substituent nuclear charges Kapton%A = 12.7014 ! g/mol, calculated from mass fractions and substituent atomic weights Kapton%density = 1.420 ! g/cm3, GEANT/NIST database Kapton%radL = 28.575 ! cm, GEANT/NIST database Kapton%nucIntL = 55.841 ! cm, GEANT/NIST database Kapton%Imean = 79.6 ! eV, GEANT/NIST database num_materials = num_materials+1 materials(num_materials) = Kapton ! Kapton-wrapped, aluminum-stabilized superconducing NbTi/Cu inflector wires InflectorCoilE821%name = 'InflectorCoilE821' InflectorCoilE821%symbol = 'InflectorCoilE821' InflectorCoilE821%Z = 15.9748 ! calculated from number fractions (stoichiometry) and substituent nuclear charges InflectorCoilE821%A = 47.1471 ! g/mol, calculated from mass fractions and substituent atomic weights InflectorCoilE821%density = 3.33358 ! g/cm3, calculated from substituent densities and volumes InflectorCoilE821%radL = 4.950 ! cm, GEANT "G4Material" InflectorCoilE821%nucIntL = 35.388 ! cm, GEANT "G4Material" InflectorCoilE821%Imean = 218.16 ! eV, GEANT "G4Material" num_materials = num_materials+1 materials(num_materials) = InflectorCoilE821 ! Same as above, with the 1.5mm Al flange and 1.5mm Al mandrel window included InflectorEndE821%name = 'InflectorEndE821' InflectorEndE821%symbol = 'InflectorEndE821' InflectorEndE821%Z = 14.5278 ! calculated from number fractions (stoichiometry) and substituent nuclear charges InflectorEndE821%A = 38.5968 ! g/mol, calculated from mass fractions and substituent atomic weights InflectorEndE821%density = 3.03140 ! g/cm3, calculated from substituent densities and volumes InflectorEndE821%radL = 6.275 ! cm, GEANT "G4Material" InflectorEndE821%nucIntL = 36.968 ! cm, GEANT "G4Material" InflectorEndE821%Imean = 193.95 ! eV, GEANT "G4Material" num_materials = num_materials+1 materials(num_materials) = InflectorEndE821 ! NbTi (1:1, MASS fractions) NbTi%name = 'NbTi' NbTi%symbol = 'NbTi' NbTi%Z = 28.461 NbTi%A = 70.39 ! g/mol, calculated from mass fractions and substituent atomic weights NbTi%density = 5.9356 ! g/cm3, calculated from substituent densities and volumes NbTi%radL = 2.072 ! cm, GEANT "G4Material" NbTi%nucIntL = 23.750 ! cm, GEANT "G4Material" NbTi%Imean = 309.87 ! eV, GEANT "G4Material" num_materials = num_materials+1 materials(num_materials) = NbTi call PrintMaterials ! call AssignMaterialPointer(material_InflectorCryoWindowUS, 'Al') call AssignMaterialPointer(material_InflectorFlangeWindowUS, 'Al') call AssignMaterialPointer(material_InflectorCoilUS, 'InflectorCoilE821') call AssignMaterialPointer(material_InflectorMandrelWindowUS,'Al') call AssignMaterialPointer(material_InflectorMandrelWindowDS,'Al') call AssignMaterialPointer(material_InflectorCoilDS, 'InflectorCoilE821') call AssignMaterialPointer(material_InflectorFlangeWindowDS, 'Al') call AssignMaterialPointer(material_InflectorCryoWindowDS, 'Al') call AssignMaterialPointer(material_fiber, 'SciFiber') ! ! call Test_InflectorE821EndScatter ! Single E821 closed inflector end RETURN END SUBROUTINE InitializeMaterials SUBROUTINE PrintMaterials IMPLICIT NONE INTEGER :: i ! Print header WRITE(*,'(/,/,a)') repeat('=',150) WRITE(*,'(2x,a)') 'AVAILABLE MATERIALS (SCATTERING/ENERGY-LOSS)' WRITE(*,'(a)') repeat('=',150) WRITE(*,'(a3,2x,2a20,6a16)') '#', 'Name', 'Symbol', 'Z', 'A (g/mol)', 'dens.(g/cm3)', 'radL (cm)', 'nucIntL (cm)', 'Imean (eV)' WRITE(*,'(a)') repeat('-',150) ! Print material DO i=1, num_materials WRITE(*,'(i3,2x,2a20,6f16.3)') i, trim(materials(i)%name), trim(materials(i)%symbol), materials(i)%Z, materials(i)%A, materials(i)%density, materials(i)%radL, materials(i)%nucIntL, materials(i)%Imean ENDDO ! Print footer WRITE(*,'(a,/)') repeat('=',150) RETURN END SUBROUTINE PrintMaterials SUBROUTINE AssignMaterialPointer(ptr,symb) IMPLICIT NONE ! Subroutine arguments TYPE (material), POINTER :: ptr ! material pointer (see above) CHARACTER(LEN=*), INTENT(IN) :: symb ! material symbol (see above) ! Subroutine internal variables INTEGER i ! Looping index ! Loop over available materials DO i=1, num_materials ! Find the material (by symbol, e.g. "Cu") IF (symb==trim(materials(i)%symbol)) THEN ! Assign the pointer ptr => materials(i) RETURN ENDIF ENDDO ! If the material is not found, print an error and exit the program WRITE(*,'(/,/,a,a,a,/,/)') 'ERROR: MATERIAL SYMBOL ', symb, ' NOT RECOGNIZED' STOP END SUBROUTINE AssignMaterialPointer SUBROUTINE Scatter(mat,ds,co,surface_normal) USE bmad USE parameters_bmad IMPLICIT NONE ! Subroutine arguments TYPE(material), INTENT(IN) :: mat ! material to traverse REAL(rp), INTENT(IN) :: ds ! scattering distance TYPE(coord_struct), INTENT(INOUT) :: co ! coordinates to scatter REAL(rp), OPTIONAL, INTENT(IN) :: surface_normal(3) ! surface orientation ! Subroutine internal variables REAL(rp) :: ds_rescaled ! scattering length depends on angle of incidence, thickness REAL(rp) :: surface_dir(3) ! surface direction (default=(0,0,1)=longitudinal) REAL(rp) :: particle_dir(3) ! incident particle direction (default=(0,0,1)=longitudinal) REAL(rp) :: tpx(3), tpy(3) ! basis vectors of transverse scattering plane REAL(rp) :: pmag, betagamma, rapidity ! relativistic kinematics (scattering depends on momentum/energy) REAL(rp) :: theta0, rand1, rand2 ! scattering variables (@PDG) REAL(rp) :: dx, dy, deltaX(3) ! change in position due to scattering REAL(rp) :: dxp, dyp, deltaXP(3) ! change in momentum direction due to scattering REAL(rp) :: ex(3), ey(3), ez(3) ! basis vectors of reference trajectory (for convenience) ! No scattering in vacuum IF (index(mat%name,'Vacuum')/=0) RETURN ! IMPORTANT: PDG formulas have length units of centimeters, whereas default BMAD length units are meters. ds_rescaled = ds*100. ! BMAD[m] => PDG[cm] ! Surface orientation (NB: With respect to reference trajectory) surface_dir = [0.,0.,1.] ! longitudinal by default if (present(surface_normal)) surface_dir = surface_normal surface_dir = surface_dir/sqrt(dot_product(surface_dir,surface_dir)) ! renormalize just to be sure ! Calculate particle direction. (NB: With respect to reference trajectory.) particle_dir(1) = cos(co%vec(4)) * sin(co%vec(2)) ! radial particle_dir(2) = sin(co%vec(4)) ! vertical particle_dir(3) = cos(co%vec(4)) * cos(co%vec(2)) ! longitudinal particle_dir = particle_dir/sqrt(dot_product(particle_dir,particle_dir)) ! Scattering only affects the transverse coordinates of the particle. Define the basis vectors that ! span the transverse plane. (See an introductory text on differential geometry if this unclear.) tpx = [ cos(co%vec(2))*cos(co%vec(4)), 0.0_rp, -sin(co%vec(2))*cos(co%vec(4)) ] tpy = [ -sin(co%vec(2))*sin(co%vec(4)), cos(co%vec(4)), -cos(co%vec(2))*sin(co%vec(4)) ] ! When the particle enters the material at an angle, we need to rescale the scattering length. The rescale factor is chosen as follows: ! If we project the particle's unit normal onto the surface's unit normal, and multiply by some scaling factor (to be determined), we MUST ! recover the thickness of the material. That is dot_product(particle_dir,surface_dir)*scale_factor = ds = material thickness. ! print '(a,3es12.4,a,3es12.4)','particle dir =', particle_dir(1:3), 'surface_dir = ', surface_dir(1:3) ! print '(a,es12.4)','ds_rescaled = ', ds_rescaled ds_rescaled = ds_rescaled/dot_product(particle_dir,surface_dir) !this gives a divide by zero. Come back to repair later ! print '(a,es12.4)','ds_rescaled = ', ds_rescaled ! The RMS scattering angle depends on the energy/momentum pmag = ( 1+co%vec(6) )*pMagic ! momentum magnitude (eV) betagamma = pmag/mmu ! Lorentz beta*gamma rapidity = asinh(betagamma) ! boost rapidity ! "Multiple scattering through small angles," Particle Data Group (2013), Eq. (31.15). Our momentum ! units here are eV, not eV/c, so there is no need to multiply the momentum by "c" in the PDG formula. theta0 = 13.6e6/( tanh(rapidity)*pmag ) * sqrt(ds_rescaled/mat%radl) * ( 1+0.038*log(ds_rescaled/mat%radl) ) !print *,'pmag = ', pmag !print *,' betagamma = ', betagamma !print *,' rapidity = ',rapidity !print *,' mat%radl = ',mat%radl !print *,' theta0 = ',theta0 ! Apply transverse scattering perpendicular to the path of the particle (taking into account that changes in position/momentum are correlated) call ran_gauss(rand1) call ran_gauss(rand2) dx = rand1*ds_rescaled*theta0/sqrt(12.) + rand2*ds_rescaled*theta0/2. ! PDG Eq. (31.22) dxp = rand2*theta0 ! PDG Eq. (31.23) ! Apply transverse scattering perpendicular to the path of the particle (taking into account that changes in position/momentum are correlated) call ran_gauss(rand1) call ran_gauss(rand2) dy = rand1*ds_rescaled*theta0/sqrt(12.) + rand2*ds_rescaled*theta0/2. ! PDG Eq. (31.22) dyp = rand2*theta0 ! PDG Eq. (31.23) ! Convert from PDG[cm] => BMAD[m] dx = dx/100. dy = dy/100. ! The coordinates of the particle are defined with respect to the reference trajectory. Therefore, we need ! to project the changes in position/momenta due to scattering we calculated above onto the reference orbit, ! and add these changes to the particle's original coordinate vector to get the particle's final coordinate ! vector after scattering. ! Basis vectors of reference trajectory (for convenience) ex = [ 1., 0., 0. ] ! radial ey = [ 0., 1., 0. ] ! vertical ez = [ 0., 0., 1. ] ! downstream ! Change in position/momentum in transverse plane of particle deltaX = [ dx, dy, 0.0_rp ] deltaXP = [ dxp, dyp, 0.0_rp ] !print '(a,6es12.4)','scatter: deltaX, deltaXP', deltaX, deltaXP ! Project the change in position/momentum in the transverse plane of the particle back onto the reference orbit deltaX = dx *( dot_product(tpx,ex)*ex + dot_product(tpx,ey)*ey + dot_product(tpx,ez)*ez ) + & dy *( dot_product(tpy,ex)*ex + dot_product(tpy,ey)*ey + dot_product(tpy,ez)*ez ) deltaXP = dxp*( dot_product(tpx,ex)*ex + dot_product(tpx,ey)*ey + dot_product(tpx,ez)*ez ) + & dyp*( dot_product(tpy,ex)*ex + dot_product(tpy,ey)*ey + dot_product(tpy,ez)*ez ) !print '(a,6es12.4)','scatter: deltaX, deltaXP', deltaX, deltaXP ! Change the particle's position co%vec(1) = co%vec(1) + deltaX(1) ! x co%vec(3) = co%vec(3) + deltaX(2) ! y co%vec(5) = co%vec(5) + deltaX(3) ! z ! Change the particle's momentum co%vec(2) = co%vec(2) + deltaXP(1) ! xp co%vec(4) = co%vec(4) + deltaXP(2) ! yp co%vec(6) = co%vec(6) ! stays same -- no energy loss RETURN END SUBROUTINE Scatter SUBROUTINE EnergyLoss(mat,ds,co,surface_normal) USE bmad USE parameters_bmad IMPLICIT NONE ! Subroutine arguments TYPE(material), INTENT(IN) :: mat REAL(rp), INTENT(IN) :: ds TYPE(coord_struct), INTENT(INOUT) :: co REAL(rp), OPTIONAL, INTENT(IN) :: surface_normal(3) ! The following are used in the "Bethe Equation," PDG (2013) Eq. (31.5), ! and the Landau-Vavilov distribution PDG (2013) Figure 31.7 real(rp) :: K = 0.307075 ! Table 31.1 [MeV*cm2/mol] real(rp) :: Wmax ! Eq. (31.4) [eV] real(rp) :: plasmaEnergy ! Table 31.1 [eV] real(rp) :: delta_2 ! Eq. (31.6) [unitless] real(rp) :: MostProbableEnergyLoss ! Eq. (31.11) [MeV] real(rp) :: xi ! Paragraph just below Eq. (31.11) [MeV] real(rp) :: transitionEnergy ! Eq. (31.47) ! Relativistic kinematics, energy loss, etc. real(rp) :: p, betagamma, rapidity real(rp) :: E, dE, dE_avg, dE_dx ! NB: The actual length of material traversed depends on the angle of ! incidence of the particle relative to the material's surface normal real(rp) :: ds_rescaled, particle_dir(3), surface_dir(3) ! Variables for generating random energy-loss fluctuations from the Landau-Vavilov probability distribution integer :: i ! index in parameters_bmad::LandauVavilov_integral(-5:100) real(rp) :: rand, frac_low, frac_high, interpolated_index ! linear interpolation between indices of table (more precise) ! No energy loss in vacuum if (mat%name=='Vacuum') return ! IMPORTANT: PDG formulas have length units of centimeters, whereas default BMAD length units are meters. ds_rescaled = ds*100. ! BMAD[m] => PDG[cm] ! Surface orientation (NB: With respect to reference trajectory) surface_dir = [0.,0.,1.] ! longitudinal by default if (present(surface_normal)) surface_dir = surface_normal surface_dir = surface_dir/sqrt(dot_product(surface_dir,surface_dir)) ! Calculate particle direction. (NB: With respect to reference trajectory.) particle_dir(1) = cos(co%vec(4)) * sin(co%vec(2)) ! radial particle_dir(2) = sin(co%vec(4)) ! vertical particle_dir(3) = cos(co%vec(4)) * cos(co%vec(2)) ! longitudinal particle_dir = particle_dir/sqrt(dot_product(particle_dir,particle_dir)) ! When the particle enters the material at an angle, we need to rescale the scattering length. The rescale factor is chosen as follows: ! If we project the particle's unit normal onto the surface's unit normal, and multiply by some scaling factor (to be determined), we MUST ! recover the thickness of the material. That is dot_product(particle_dir,surface_dir)*scale_factor = ds = material thickness. ds_rescaled = ds_rescaled/dot_product(particle_dir,surface_dir) ! Energy loss depends on energy/momentum p = ( 1+co%vec(6) )*pMagic ! momentum magnitude (eV) betagamma = p/mmu ! Lorentz betagamma rapidity = asinh(betagamma) ! hyperbolic rotation angle E = cosh(rapidity)*mmu ! energy (eV) ! The Bethe Equation (31.5) gives the average energy loss, and the Landau-Vavilov distribution gives the energy-loss distribution. ! According to the PDG, "It is the very rare high-energy-transfer collisions ... that drive the mean of the distribution into the ! tail.... The most-probable energy loss should be used." The most-probable energy loss is given by PDG Eq. (31.11). Let's do ! a few preliminary helper calculations first, then assemble the pieces. (We'll calculate the Bethe Equation as a diagnostic.) ! The maximum energy transferred in a single collision is given be Eq. (31.4) Wmax = 2*me*betagamma**2/( 1 + 2*cosh(rapidity)*me/mmu + (me/mmu)**2 ) ! eV ! "Plasma energy," Table 31.1 plasmaEnergy = sqrt( mat%density * mat%Z / mat%A )*28.816 ! eV ! "Density effect," Eq. (31.6) delta_2 = log(plasmaEnergy/mat%Imean) + log(betagamma) - 0.5 ! dimensionless ! "Bethe Equation," Eq. (31.5). Note: The PDG detector thickness "x" has units of [g/cm2], so the units of "dE/dx" are [MeV/(g/cm2)]. ! Dividing "x" by the density will make it have units of [cm], and will make dE/dx have units of [MeV/cm] dE_dx = -K*mat%Z/mat%A/tanh(rapidity)**2 * ( 0.5*log(2*me*betagamma**2*Wmax/mat%Imean**2) - tanh(rapidity)**2 - delta_2 ) * mat%density ! MeV/cm ! Calculate the average energy loss from the Bethe Equation given the scattering length. dE_avg = dE_dx * ds_rescaled * (1.e6) ! eV ! So far, we have only calculated the /average/ energy shift from the Bethe Equation. How do we include ! fluctuations? The answer is shown in PDG (2013) Figure 31.7, i.e. the Landa-Vavilov distribution. ! "Fluctuations in energy loss," Eq. (31.11). Again, the PDG material thickness needs to have units of [g/cm2], so we multiply by the density. xi = ( K/2. )*( mat%Z/mat%A )*( ds_rescaled*mat%density/tanh(rapidity)**2 ) ![MeV] See text immediately beneath Eq. (31.11) ! Most probable energy loss for the Landau-Vavilov distribution, Eq. (31.11) MostProbableEnergyLoss = xi*( log(2*me*betagamma**2/mat%Imean) + log(xi/mat%Imean) + 0.200 - tanh(rapidity)**2 - 2.*delta_2 ) ! [MeV] dE=0 ! The method to generate a random energy loss based on the Landau-Vavilov distribution is neat: Look at Fig. (31.7) of the ! PDG (2013), which shows the Landau-Vavilov distribution. Now imagine taking the integral of the distribuiton. Since the ! distribution is always positive, the integral starts at zero and increases monotonically. In terms of the integral, the peak of ! the distribution will be where the integral is changing the fastest. So here's what we do: We throw a uniform random number, find ! where it is in the (tabulated) integral, and map it to an energy. The uniform random number represents the value of the integral ! itself. The corresponding index of the random number in the integral table therefore corresponds to the independent variable. ! Then we simply rescale the index to convert to an energy. (In order to be even more precise, we linearly interpolate between indices.) call ran_uniform(rand) ! This number corresponds to the value of the integral of the normalized Landau-Vavilov distribution DO i=-5, 100-1 ! Here are the bounds of the integral table. The index "i" corresponds to the independent variable IF ( LandauVavilov_integral(i)