! GeDalitzDecay.f90 ! ! Dalitz decay A->B+l+l ! ! References: ! N.Kroll, W.Wada, Phys. Rev. 98 No.2 1955 (1355) ! D.W. Joseph, Nuovo Cimento XVI, No.6 1960 (2021) ! C.H. Lai, C. Quigg, Preprint FNAL FN-296, Batavia 1976 ! ! Detlef Irmscher, Thomas Ullrich, Uni Heidelberg, September 1993 ! Last update: 24.3.95 tu - major update of form factor handling ! - enabled for muons ! subroutine GeDalitzDecay (treeIndex) use genesis_m, self=>GeDalitzDecay use f90_kind use mom_m use FormDecay_m implicit none integer, intent(in) :: treeIndex type(tree), pointer :: Dalitz, other, lepton1, lepton2 logical, save :: init = .true. integer, save :: angleid integer :: ibody, ilist integer :: ptid, ibinx integer, parameter :: nxbins = 10000 real :: boffx, rand(2) real(kind=double) :: mass_a, mass_b, mass_l, mvirt_2 real(kind=double) :: delta, epsilon real(kind=double) :: xd, dalval, bracket real(kind=double) :: costheta, sintheta, phi, cosphi, sinphi real(kind=double) :: e1, p1, e3, p3 type(mom4) :: p4boost real :: xlow,xhigh real, external :: GeCos2 ! ! Initialize ! Fill histograms with the proper pair mass distributions. ! if (init) then ptid = workingID do ilist = 1, size(decayList) ! identify dalitz decays in decayList if (.not.GeIsDalitzDecay(decayList(ilist))) cycle if (.not.decayList(ilist)%enabled) cycle ! get free hbook id ptid = GeFreeHbookId(ptid) ! book 1-dim histogram for x-dalitz boffx = 0.5/float(nxbins) call hbook1 (ptid, 'dalitz decay '//decayList(ilist)%name, & nxbins, 0.+boffx, 1.+boffx, 0.) ! get masses of particles mass_a = decayList(ilist)%parent%mass mass_l = 0. mass_b = 0. do ibody = 1, 3 if (abs(decayList(ilist)%child(ibody)%id) == 12 .or. & ! electron abs(decayList(ilist)%child(ibody)%id) == 14 ) then ! muon mass_l = decayList(ilist)%child(ibody)%mass else mass_b = decayList(ilist)%child(ibody)%mass endif enddo ! parameters for dalitz decay delta = mass_b**2/mass_a**2 epsilon = mass_l**2/mass_a**2 ! fill histogram do ibinx = 1, nxbins xd = 1./float(nxbins)*ibinx if (xd <= 4.*epsilon) cycle ! kinem. limit bracket = (1. + xd/(1. - delta))**2 & - 4.*xd/(1. - delta)**2 if (bracket < 0.) cycle ! kinem. limit dalval = bracket**1.5 /xd * & sqrt(1-4*epsilon/xd)*(1 + 2*epsilon/xd)* & GeFormFactor(xd*mass_a**2,decayList(ilist)) call HFILL (ptid, sngl(xd), 0., sngl(dalval)) enddo ! keep hbook id for decay decayList(ilist)%HbookId = ptid ! copy histogram to keep original distribution before integration call HCOPY(decayList(ilist)%HbookId, GeFreeHbookId(originalID),' ') enddo ! ! generate 1+cos**2 angular distribution ! angleid = GeFreeHbookId(ptid) xlow =-1 xhigh = 1 call hbfun1(angleid,'1+cos(t)**2',1000,xlow,xhigh,GeCos2) init = .false. endif ! Set up tree, identify decay childs Dalitz => treeList(treeIndex) ! Dalitz = pointer to parent particle do ibody = 1, 3 select case (Dalitz%decaytype%child(ibody)%id) case (-12) lepton1 => treeList(Dalitz%child(ibody)) ! pointer to 1st lepton case (-14) lepton1 => treeList(Dalitz%child(ibody)) case (12) lepton2 => treeList(Dalitz%child(ibody)) ! pointer to 2nd lepton case (14) lepton2 => treeList(Dalitz%child(ibody)) case default other => treeList(Dalitz%child(ibody)) ! pointer to 3rd part. end select enddo ! get masses of decay particles mass_a = Dalitz%parttype%mass mass_b = other%parttype%mass mass_l = lepton1%parttype%mass ! get Dalitz x from histogram do xd = HRNDM1 (Dalitz%decaytype%HbookId) ! calculate mass of virtual photon mvirt_2 = xd*mass_a**2 ! check kinematics if (mass_a - mass_b - 2.*mass_l .gt. sqrt(mvirt_2) .and. & sqrt(mvirt_2)/2. .gt. mass_l) exit enddo ! decay virtual photon in two leptons in their CMS ! calculate CM energies and momenta of children e1 = sqrt(mvirt_2)/2. p1 = sqrt((e1 + mass_l)*(e1 - mass_l)) ! the angle of the leptons to the z axis is generated uniform if (mass_b < 0.01) then ! true Dalitz decay gamma e+ e- costheta = hrndm1 (angleid) else call random_number(rand) costheta = 2.*rand(1) - 1. endif sintheta = sqrt((1. + costheta)*(1. - costheta)) phi = TWOPI * rand(2) sinphi = sin(phi) cosphi = cos(phi) ! fill 4-vectors of leptons lepton1%p4%E = e1 lepton1%p4%p%px = p1*sintheta*cosphi lepton1%p4%p%py = p1*sintheta*sinphi lepton1%p4%p%pz = p1*costheta lepton2%p4%p = -lepton1%p4%p lepton2%p4%E = e1 ! calculate components of non-lepton in CMS of parent (pion, eta...) e3 = (mass_a**2 + mass_b**2 - mvirt_2)/(2.*mass_a) p3 = sqrt((e3 + mass_b)*(e3 - mass_b)) call random_number(rand) costheta = 2.*rand(1) - 1. sintheta = sqrt((1. + costheta)*(1. - costheta)) phi = TWOPI * rand(2) sinphi = sin(phi) cosphi = cos(phi) ! fill 4-vector of non-lepton other%p4%E = e3 other%p4%p%px = p3*sintheta*cosphi other%p4%p%py = p3*sintheta*sinphi other%p4%p%pz = p3*costheta ! rotate lepton vectors into coordinate system of virtual photon sintheta = -sintheta cosphi = -cosphi sinphi = -sinphi lepton1%p4%p = GeP3rot(lepton1%p4%p, costheta, sintheta, cosphi, sinphi) lepton2%p4%p = GeP3rot(lepton2%p4%p, costheta, sintheta, cosphi, sinphi) ! boost leptons with virtual photon p4boost%E = sqrt(p3**2 + mvirt_2) p4boost%p = -other%p4%p lepton1%p4 = GeBoost(lepton1%p4, p4boost) lepton2%p4 = GeBoost(lepton2%p4, p4boost) ! boost into lab system lepton1%p4 = GeBoost(lepton1%p4, Dalitz%p4) lepton2%p4 = GeBoost(lepton2%p4, Dalitz%p4) other%p4 = GeBoost(other%p4, Dalitz%p4) !contains ! ! ! ! ! Rotate vector into new coordinate system ! ! ! function GeP3rot (p3old, costh, sinth, cosph, sinph) ! implicit none ! type(mom3) :: GeP3rot ! type (mom3), intent(in) :: p3old ! real (kind=double), intent(in) :: costh, sinth, cosph, sinph ! ! ! make new vector by rotating old one with anti-clockwise Euler angles ! GeP3rot%px = p3old%px*costh*cosph - p3old%py*sinph + p3old%pz*sinth*cosph ! GeP3rot%py = p3old%px*costh*sinph + p3old%py*cosph + p3old%pz*sinth*sinph ! GeP3rot%pz = -p3old%px*sinth + p3old%pz*costh ! end function GeP3rot end subroutine GeDalitzDecay