c eeloss.f c get electron eloss implicit double precision (a-z) c character*8 keyin, keymat c real me, mpi, ipot, ne, na, min, k, mmu, mp c=2.99792458e10 ! cm/s pi=3.1415926 mp=938.27231 me=0.51099906 mpi=139.5678 mmu=105.658389 na=6.022136e23 de=5.0989e-25 ! MeV-cm**2 = 4pi r_e**2 mc**2 k=.307075 !MeV/gm/cm**2 min=me zin=1. tin=845. tin=4045. tin=400. tin=3245. c write(6,*) 'Enter electron kinetic energy in MeV' c read(5,*) tin ein=tin+min pin=sqrt(ein**2-min**2) tau=tin/me beta=pin/ein gamma=1./sqrt(1-beta**2) tmax=2.*me*beta**2*gamma**2 tmax=tmax/(1.+2.*gamma*me/min+(me/min)**2) c beta=((tau*(tau+2.))**2)/(tau+1.) rho=19.35 ! Tungsten z=74. a=183.85 t=10./rho rho=8.96 ! Copper z=29. a=63.55 rho=0.0708 ! LH2 g/cm**3 z=1. a=1. thick=rho*1. ! 1cm thickness in g/cm*2 c if(min.eq.me) keyin='electron' c if(min.eq.mpi) keyin='pion' c if(min.eq.mp) keyin='proton' adjI=16.0e-6*z**0.9 ! MeV i=z*(9.76+58.8*z**-1.19) ! eV if(abs(a-1.).lt.0.02 .and. abs(z-1.).lt.0.02) then adjI=21.8e-6 ! H2 is special i=21.8 endif ne=z*na*rho/a x0=716.4*a/(z*(z+1.)*log(287./sqrt(z))) ! rad. length in g/cm**2 c now use Sternheimer PRB 3, 3681 (71) hnu=28.8*sqrt(rho*z/a) ! ev cbar=2.*log(i/hnu)+1. xa=cbar/4.606 m=3.0 c the following is for solids & liquids...gases have to be treated separately if(i.lt.100.) then x1bar=2. x0bar=0.2 if(cbar.ge. 3.681) x0bar=0.326*cbar-1.0 else x1bar=3. x0bar=0.2 if(cbar.ge.5.215) x0bar=0.326*cbar-1.5 endif littlea=4.606*(xa-x0bar)/(x1bar-x0bar)**3 x=log10(pin/min) delta=4.606*(x-xa) + 4.606*(xa-x0bar)*((x1bar-x)/(x1-x0bar))**m biga=0.1536*z/a bigb=log(me*1.e12/i**2) wmax=pin**2/(min*(min/2./me + me/2./min + ein/min)) sdedx=biga/beta**2*(bigb+0.693+2.*log(pin/min) * +log(wmax)-2.*beta**2-delta) sdedx_mprob=biga*thick/beta**2*(bigb+0.891+2.*log(pin/min) * +log(biga*thick/beta**2)-beta**2-delta) c now Bethe-Bloch bb=de*ne*(zin/beta)**2 bb=bb*(log(2.*me*beta**2*gamma**2/adjI) - beta**2) c delta=log(28.816*sqrt(rho*z/a)/adjI) + log(beta*gamma) - 0.5 bbpdg=k*z/a*(zin/beta)**2 bbpdg=bbpdg*(0.5*log(2.*me*beta**2*gamma**2*tmax/adjI**2) * - beta**2 - delta/2.) brems_pdg=ein/x0 ! ~bremss eloss for electrons in MeV/(g/cm^2) b=beta c beta=(tau*(tau+2.))**2/(tau+1.)/c f=1.-beta**2+(tau**2/8.-(2.*tau+1.)*log(2.))/(tau+1.)**2 dedx=log(((tau+2.)*tau**2)/(2.*(adjI/0.511)**2)) +f dedx=dedx*0.1535*rho*z/a/beta**2 beta=b fernow=de/2.*ne*(2.*log(2*me/adjI) + 3.*log(gamma) -1.95) ! MeV/cm krane=log( (tin*beta*beta*(tin+me)**2) / (2.*adjI**2*me) ) krane=krane+1.-beta**2 + (1.-sqrt(1.-beta**2))/8. krane=krane-(2.*sqrt(1.-beta**2) -1. +beta**2)*log(2.) krane=krane*(1.439976e-13)**2*2.*pi*na*z*rho/me/beta**2/a power=10.*sdedx*rho*1. ! I(uA)*dedx*rho*t=beam power into target write(6,*) ' i, hnu, cbar, xa:', i, hnu, cbar, xa write(6,*) ' x1bar, x0bar, a, x:', x1bar, x0bar, littlea, x write(6,*) ' A,B, delta: ', biga, bigb, delta write(6,*) ' Stern. average dedx (MeV/g/cm^2) = ', sdedx, * ' or in MeV/cm =', sdedx*rho write(6,*) ' Stern most prob dedx (MeV/g/cm^2) = ', * sdedx_mprob/rho, ' or in MeV/cm =', sdedx_mprob write(6,*) write(6,*) 'beta, tau, te, ne: ', beta, tau, tin, ne write(6,*) 'rho, z, a, adjI(eV): ', rho, z, a, adjI*1.e6 write(6,*) 'min, ein, tin, pin:', min, ein, tin, pin write(6,*) ' gam, gamma', ein/min, 1./sqrt(1-(pin/ein)**2) write(6,*) write(6,*) 'dedx, f (TKH)= ', dedx, f, dedx/rho write(6,*) ' krane dedx = ', krane, krane/rho write(6,*) 'fernow dedx (MeV/cm)= ', fernow, fernow/rho * , ' (MeV/g/cm^2)' c write(6,*) 'Bethe-Bloch de/dx (MeV/cm) = ', bb write(6,*) 'Bethe-Bloch de/dx (MeV/(g/cm**2)) ala pdg = ', bbpdg write(6,*) 'Bethe-Bloch de/dx (MeV/cm) ala pdg = ', bbpdg*rho write(6,*) write(6,*) 'Radiation Length (g/cm**2) ala pdg = ', x0 write(6,*) 'Radiation Length (cm) ala pdg = ', x0/rho write(6,*) 'Brems. de/dx (MeV/(g/cm^2)) ala pdg = ', brems_pdg write(6,*) 'Brems. de/dx (MeV/cm) ala pdg = ', brems_pdg*rho write(6,*) 'Brems. cone angle (deg) = ', 1./gamma*180./3.1415926 c 6.25e12 electrons/s per uA lum=6.25e12*na/a*thick/4./pi write(6,*) 'Luminosity is ', lum, " cm^-2 s-1 per cm per uA" c write(6,10) keyin, min, tin, pin, beta, gamma write(6,10) min, tin, pin, beta, gamma 10 format(/,"Incident Particle: ", !a10, * /, " Mass=", f8.3, " MeV, K.E.=", * f6.1, " MeV, momentum=", f6.1, " MeV/c", /, * " beta=", f9.7, ", gamma=", f8.3,/) write(6,11) z, a, rho, i, x0, x0/rho 11 format("Material:",/, " Z=", f4.0, ", A=", f9.4, ", rho=", f9.4, * " g/cm^3",/, " ion. pot.=", f7.2, " eV, r.l.=", f7.2, * " g/cm^2, or ", f7.2, " cm",/) write(6,12) sdedx, sdedx*rho, brems_pdg, brems_pdg*rho, * 1./gamma*180./3.1415926, sdedx+brems_pdg, * sdedx*rho+brems_pdg*rho 12 format("Ionization Eloss: ", f6.2, " MeV/g/cm^2, or ", f6.2, * " MeV/cm",/, * "Bremsstrahlung Eloss: ", f6.2, " MeV/g/cm^2, or ", * f6.2, " MeV/cm, ", /, " Cone Angle: ", f7.3, " deg.",/, * "Total Eloss: ", f6.2, " MeV/g/cm^2, or ", * f9.4, " MeV/cm",/ ) write(6,13) power, power*4, power*15 13 format( * " Target Power per cm per 10 uA: ", f6.1, " Watts", /, * " Target Power per 4 cm per 10 uA: ", f6.1, " Watts", /, * " Target Power per 15 cm per 10 uA: ", f6.1, " Watts",/) end