!  GeFormFactor.f90
!
!  Calculates the form factor for Dalitz decays A->B+l+l
!  according to either the Vector Dominance Model or data.
!  Returns: |F(q^2)|^2
!
!  Warning: this functions assumes that the given decay
!           is a Dalitz decay. No checks are performed.
!
!  References:
!  L.G. Landsberg, Physics Reports 128 No.6 (1985) 301-376. 
!
!  Thomas Ullrich, Uni Heidelberg, March 1995
!
real function GeFormFactor(q2, decaytype)
    use f90_kind
    use genesis_m
    use FormDecay_m, self=>GeFormFactor
    implicit none
    type(decay)                  :: decaytype
    real(kind=double)            :: q2
    integer                      :: i
    complex(kind=double)         :: F
    logical, save                :: initialize = .true.
    real(kind=double), parameter :: w_pi0(3)   = (/ 0.5,    0.5,     0. /)
    real(kind=double), parameter :: w_eta(3)   = (/ 1.2649, 0.1405, -0.4055 /)
    real(kind=double), parameter :: w_etap(3)  = (/ 0.6985, 0.0776,  0.2239 /)
    real(kind=double), parameter :: w_omega(3) = (/ 1.,     0.,      0. /)
    real(kind=double), save      :: mass(3)    = (/ 0.    , 0.,      0. /)
    real(kind=double), save      :: width(3)   = (/ 0.    , 0.,      0. /)
    
    ! Get masses and width of vector mesons
    if (initialize) then
       i = GeParticleIndex("rho")
       mass(1)  = particleList(i)%mass
       width(1) = particleList(i)%width
       i = GeParticleIndex("omega")
       mass(2)  = particleList(i)%mass
       width(2) = particleList(i)%width
       i = GeParticleIndex("phi")
       mass(3)  = particleList(i)%mass
       width(3) = particleList(i)%width
       initialize = .false.
       do i = 1, 3
          if (mass(i) <= 0 .or. width(i) <= 0 ) &
          call GeErrorMessage("GeFormFactor",                       &
	                      "cannot get proper mass/width for "// &
                              "'"//trim(particleList(i)%name)//"'", &
			      fatal=.true.)			     
       enddo
    endif
    
    
    select case (decaytype%parent%id)   ! identify decay by name of parent
    !
    ! pi0   -> e+ e- gamma
    !
    case (110)                          
       select case(decaytype%formFactor)
       case(USE_NONE)
          GeFormFactor = 1.
       case(USE_VDM)                    
	  F = (0.d0,0.d0)
	  do i = 1, 3
	     F = F + w_pi0(i)*mass(i)**2/ &
	         (mass(i)**2-q2-(0.,1.)*width(i)*mass(i))
	  enddo
          GeFormFactor = real(F*conjg(F))
       case(USE_LEPTON_G)
          GeFormFactor = (1./(1-5.5*q2))**2       ! pole approximation
       end select
    !
    ! eta   -> l+ l- gamma
    !
    case (220)                          
       select case(decaytype%formFactor)
       case(USE_NONE)
          GeFormFactor = 1.
       case(USE_VDM)                    
	  F = (0.d0,0.d0)
	  do i = 1, 3
	     F = F + w_eta(i)*mass(i)**2/ &
	         (mass(i)**2-q2-(0.,1.)*width(i)*mass(i))
	  enddo
          GeFormFactor = real(F*conjg(F))
       case(USE_LEPTON_G)
          GeFormFactor = (1./(1.-1.9*q2))**2      ! pole approximation
       end select
    !
    ! eta'  -> l+ l- gamma
    !
    case (330)                          
       select case(decaytype%formFactor)
       case(USE_NONE)
          GeFormFactor = 1.
       case(USE_VDM)
	  F = (0.d0,0.d0)
	  do i = 1, 3
	     F = F + w_etap(i)*mass(i)**2/ &
	         (mass(i)**2-q2-(0.,1.)*width(i)*mass(i))
	  enddo
          GeFormFactor = real(F*conjg(F))
       case(USE_LEPTON_G)
          GeFormFactor = 0.764**4/ &             ! Breit-Wigner fitted to data
	  ((0.764**2-q2)**2+0.1020**2*0.764**2)
       end select
    !
    ! omega -> l+ l- pi0
    !
    case (221)                          
       select case(decaytype%formFactor)
       case(USE_NONE)
          GeFormFactor = 1.
       case(USE_VDM)                    
	  F = (0.d0,0.d0)
	  do i = 1, 3
	     F = F + w_omega(i)*mass(i)**2/ &
	         (mass(i)**2-q2-(0.,1.)*width(i)*mass(i))
	  enddo
          GeFormFactor = real(F*conjg(F))
       case(USE_LEPTON_G)
          GeFormFactor = 0.6519**4/ &            ! Breit-Wigner fitted to data
	  ((0.6519**2-q2)**2+0.04198**2*.6519**2)
       end select
    !
    ! not implemented 
    !
    case default
       GeFormFactor = 1.
    end select
    
end function GeFormFactor