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 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) Q = coord%spin !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 = coord 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 co_electron%ix_user = 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 real(rp):: energymaxr= 5.2e7 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 !print '(a,3es12.4,a,3es12.4)',' P=',P,' Q=',Q 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) !write(901,'(3es12.4)')energy,theta,phi !obtain electron momentum with respect to rest frame polarizaton pem = sqrt(energy**2-me**2) !dlr 20180713 Pe(1) = pem*cos(theta) !dlr 20180713 Pe(2) = pem*sin(theta)*cos(phi) !dlr 20180713 Pe(3) = pem*sin(theta)*sin(phi) Pe(3) = pem*cos(theta) Pe(1) = pem*sin(theta)*cos(phi) Pe(2) = pem*sin(theta)*sin(phi) y=energy/energymaxr asym = (2*y-1)/(3-2*y) write(112,'(8es12.4)')energy,theta,phi,pe(1:3),y,asym !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)) !dlr 20180713 Ry = reshape((/cos(thetar-pi/2),0.0_rp,-sin(thetar-pi/2),0.0_rp,1.0_rp,0.0_rp,sin(thetar-pi/2),0.0_rp,cos(thetar-pi/2)/),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(111,'(2es12.4,es12.4,4es12.4,2es12.4)') y,asym,thetal,pe4,energy,theta !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 precision_def use random_mod use parameters_bmad 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 ! enforce the minimum energy (can't have energy < rest mass) if (energy < me) energy = me 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 real(rp) a real(rp) energy_local, y_local energy_local = energy y_local = y ! a=0.5*(2*y-1)/(3-2*y) a=(2*y-1)/(3-2*y) !dlr 20180724 ! n = y**2*(3-2*y) !dlr 20180724 ! MuonDecayAngDist = 0.5*(1-x)*(1-a-a*x) ! MuonDecayAngDist = (1+x + 0.5*a*(1-x**2))/2 !dlr 20180716 MuonDecayAngDist = (1+x + 0.5*a*(-1+x**2))/2 !dlr 20181114 end function end subroutine electronDist