subroutine DecayElectron(coord,co_electron,n) use bmad use parameters_bmad use spin_mod !main subroutine implicit none type(coord_struct) coord,co_electron !input and output coord_struct real(rp),dimension(6):: evector !phase space coordinates of the generated electron real(rp),dimension(3)::P,Q !muon momentum and polarization vector in lab frame real(rp),dimension(3)::Pe !generated electron momentum in lab frame type(pauli_struct),dimension(3):: sig !pauli matrices real(rp):: p0c = 3.094353005E9 !reference/magical momentum real(rp):: pm !muon momentum magnitude in lab frame real(rp):: gammae,gammam !gamma for the electron and muon real(rp):: betae,betam !beta for the electron and muon integer j,n !n is the index of muon ! extract momentum and polarization from coord_struct P(1) = coord%vec(2)*p0c P(2) = coord%vec(4)*p0c pm = coord%vec(6)*p0c+p0c P(3) = sqrt(pm**2-P(1)**2-P(2)**2) sig(1)%sigma = reshape((/(0.0_rp,0.0_rp),(1._rp,0.0_rp),(1._rp,0.0_rp),(0.0_rp,0.0_rp)/),(/2,2/)) sig(2)%sigma = reshape((/(0.0_rp,0.0_rp),(0.0_rp,1._rp),(0.0_rp,-1._rp),(0.0_rp,0.0_rp)/),(/2,2/)) sig(3)%sigma = reshape((/(1._rp,0.0_rp),(0.0_rp,0.0_rp),(0.0_rp,0.0_rp),(-1._rp,0.0_rp)/),(/2,2/)) do j=1,3 Q(j) = dot_product(coord%spin,matmul(coord%spin,transpose(sig(j)%sigma))) end do !generate electron momentum call GenerateDist(P,Q,Pe) !convert to phase space coordinates evector(1) = coord%vec(1) evector(3) = coord%vec(3) gammam = sqrt(dot_product(P,P)+mmu**2)/mmu gammae = sqrt(dot_product(Pe,Pe)+me**2)/me betae = sqrt(dot_product(Pe,Pe))/me/gammae betam = pm/mmu/gammam evector(5) = coord%vec(5)*betae/betam evector(2) = Pe(1)/p0c evector(4) = Pe(2)/p0c evector(6) = (sqrt(dot_product(Pe,Pe))-p0c)/p0c co_electron%vec = evector co_electron%t = coord%t co_electron%s = coord%s co_electron%species = positron$ !electron$ co_electron%p0c = coord%p0c co_electron%direction = coord%direction co_electron%location = coord%location co_electron%beta = betae co_electron%ix_ele = coord%ix_ele co_electron%state = alive$ co_electron%path_len = n !store the index of muon in path_len !print *, 'charge',coord%charge,'p0c',coord%p0c,'beta',coord%beta,'direction',coord%direction,'species',coord%species,'location',coord%location, & ! 'ix_ele',coord%ix_ele,'state',coord%state return end subroutine DecayElectron !+ !subroutine that generates the decay electron's momentum according to !the given momentum and polarization of the muon subroutine GenerateDist(P,Q,Pe) use bmad use parameters_bmad implicit none real(rp),dimension(3)::P,Q !muon momentum and polarization vector in lab frame real(rp),dimension(3)::Pe !generated electron momentum real(rp),dimension(4)::Pe4 !electron four-momentum real(rp),dimension(3,3):: Rz,Ry !rotation matrices real(rp),dimension(4,4):: Lz !matrix for Lorentz boost real(rp) a,b,c !rotation angles real(rp) thetal,thetar !polarization angle in lab frame and rest frame real(rp) cosThetal real(rp) gamma,beta real(rp) energy,theta,phi !parameters of the generated electron real(rp) pem !magnitude of electron momentum in rest frame real(rp) y,asym !fractional energy and asymmetry real(rp):: energyMaxl = 3.094353005E9 !maximum energy in lab frame integer i !rotate P onto xz plane if (P(1) == 0) then ! take care of the special case if (P(2) == 0) then a = 0 else a = -P(2)/abs(P(2))*pi/2 end if else a = -atan(P(2)/P(1)) end if Rz = reshape((/cos(a),sin(a),0.0_rp,-sin(a),cos(a),0.0_rp,0.0_rp,0.0_rp,1.0_rp/),shape(Rz)) P = matmul(P,transpose(Rz)) Q = matmul(Q,transpose(Rz)) !rotate P to align it with the z axis b = -atan(P(1)/P(3)) Ry = reshape((/cos(b),0.0_rp,-sin(b),0.0_rp,1.0_rp,0.0_rp,sin(b),0.0_rp,cos(b)/),shape(Ry)) P = matmul(P,transpose(Ry)) Q = matmul(Q,transpose(Ry)) !for convenience, rotate Q around z axis to ensure it is on xz plane if (Q(1) == 0) then ! take care of the special case if (Q(2) == 0) then c = 0 else c = -Q(2)/abs(Q(2))*pi/2 end if else c = -atan(Q(2)/Q(1)) end if Rz = reshape((/cos(c),sin(c),0.0_rp,-sin(c),cos(c),0.0_rp,0.0_rp,0.0_rp,1.0_rp/),shape(Rz)) Q = matmul(Q,transpose(Rz)) !determine the polarization angle thetal = acos(dot_product(P,Q)/sqrt(dot_product(P,P)*dot_product(Q,Q))) if (Q(1) < 0) thetal = -thetal !attach the right sign gamma = sqrt(dot_product(P,P)+mmu**2)/mmu !boost to the rest frame thetar = atan(gamma*tan(thetal)) if (Q(3) < 0) thetar = thetar+pi !get the proper angle !generate energy and angular distribution of the decay electron call electronDist(energy,theta,phi) !obtain electron momentum with respect to rest frame polarizaton pem = sqrt(energy**2-me**2) Pe(1) = pem*cos(theta) Pe(2) = pem*sin(theta)*cos(phi) Pe(3) = pem*sin(theta)*sin(phi) !obtain elctron momentum with respect to lab frame muon momentum Ry = reshape((/cos(thetar),0.0_rp,-sin(thetar),0.0_rp,1.0_rp,0.0_rp,sin(thetar),0.0_rp,cos(thetar)/),shape(Ry)) Pe = matmul(Pe,transpose(Ry)) Pe4 = (/energy,Pe/) !boost back to lab frame beta = sqrt(dot_product(P,P))/mmu/gamma Lz = reshape((/gamma,0.0_rp,0.0_rp,beta*gamma,0.0_rp,1.0_rp,0.0_rp,0.0_rp,0.0_rp,0.0_rp,1.0_rp,0.0_rp,beta*gamma,0.0_rp,0.0_rp,gamma/), shape(Lz)) Pe4 = matmul(Pe4,transpose(Lz)) Pe = Pe4(2:4) !Calculate and write the parameters for testing !y = Pe4(1)/energyMaxl !asym = (2*y-1)/(3-2*y) !write(11,'(2es12.4)') y,asym !rotate back Rz = reshape((/cos(-c),sin(-c),0.0_rp,-sin(-c),cos(-c),0.0_rp,0.0_rp,0.0_rp,1.0_rp/),shape(Rz)) Pe = matmul(Pe,transpose(Rz)) Ry = reshape((/cos(-b),0.0_rp,-sin(-b),0.0_rp,1.0_rp,0.0_rp,sin(-b),0.0_rp,cos(-b)/),shape(Ry)) Pe = matmul(Pe,transpose(Ry)) Rz = reshape((/cos(-a),sin(-a),0.0_rp,-sin(-a),cos(-a),0.0_rp,0.0_rp,0.0_rp,1.0_rp/),shape(Rz)) Pe = matmul(Pe,transpose(Rz)) return end subroutine GenerateDist !+ !subroutine that generates the energy and the angles thea and phi of the decay electron !in the muon rest frame with respect to muon polarization subroutine electronDist(energy,theta,phi) use bmad use random_mod integer :: seed = 1. integer i real(rp) harvest real(rp) x,y,costheta,theta,phi real(rp):: energy,energymax = 52.8e6 real(rp) x1, x2, tol ! real(rp) MuonToElectronEnergy, MuonDecayAngDist !generate distribution of electrons assuming muon at rest and polarization P in the +z direction !initialize randomizer's seed (harvest) ! ran_seed = 0 <--> seed randomizer w/ current CPU time !call ran_seed_put(seed) call ran_uniform(harvest) x=harvest x1=0. x2=1. tol = 0.0001 y = inverse(MuonToElectronEnergy,x, x1,x2,tol) energy = energymax*y call ran_uniform(harvest) x=harvest costheta = inverse(MuonDecayAngDist,x,-x2,x2,tol) theta = acos(costheta) call ran_uniform(harvest) phi = harvest*2*pi return contains !+ ! Function MuonToElectronEnergy(x) ! ! Routine to compute the probability that the muon decay electron will have fraction energy x=E/Emax ! ! Defining relation: ! MuonToElectronEnergy(x) = x**3*(2-x) ! function MuonToElectronEnergy(x) use precision_def implicit none real(rp) MuonToElectronEnergy, x MuonToElectronEnergy = x**3 * (2-x) energy = x end function !+ ! Function MuonDecayAngDist(x,y) ! ! Routine to compute the probability that the decay electron momenum is at and angle acos(x) with respect to muon spin ! ! Defining relation ! MuonDecayAngDist(x) = 0.5*(1-x)*(1-a-a*x) ! where a = 0.5*(2*y-1)/(3-2*y) ! function MuonDecayAngDist(x) use precision_def implicit none real(rp) MuonDecayAngDist, x,y real(rp) a a=0.5*(2*energy-1)/(3-2*energy) MuonDecayAngDist = 0.5*(1-x)*(1-a-a*x) end function end subroutine electronDist