!  GeGetPt.f90
!
!  Returns a pt according to a pt distribution.
!
!  Detlef Irmscher, Thomas Ullrich, Uni Heidelberg, September 1993
!  Last update: 08/05/96 13:10:39 tu
!
function GeGetPt (ipart, weight)
    use genesis_m, self=>GeGetPt
    use f90_kind
    implicit none
    real(kind=double)              :: GeGetPt
    integer, intent(in)            :: ipart
    logical, save                  :: init = .true.
    real                           :: pt, mass, mt
    real                           :: wt
    integer                        :: iploop, ibin
    integer                        :: ptid
    integer, save                  :: pionPtOriginalId, pionIndex
    integer, save                  :: etaIndex,rhoIndex,omegaIndex
    integer                        :: pionWorkingId
    real(kind=double), intent(out) :: weight

    if (init) then
        !
	!  First fill the pion0 pt distribution.
	!  
	pionPtOriginalId = GeFreeHbookId(originalID)      ! id not to be integrated
	pionWorkingId    = GeFreeHbookId(workingID)       ! id which will be integrated
	pionIndex = GeParticleIndex("pion0")
	etaIndex = GeParticleIndex("eta")
        rhoIndex = GeParticleIndex("rho")
	omegaIndex = GeParticleIndex("omega")

	if (particleList(pionIndex)%ptHbookId == 0) then
	   !
	   !  No pion pt distribution provided by the user
	   !
           call HBOOK1 (pionWorkingId, 'pt '//particleList(pionIndex)%name,    &
                        1000, 0.005, 4.005, 0.)
           if (projectile%a.eq.208) then
              print *,'GeGetPt> use default pt-spectra for Pb-beam'
	      do ibin = 1, 1000
	         pt   = 4./1000.*(float(ibin)-.5)
                 mass = particleList(pionIndex)%mass
                 mt   = sqrt(pt**2+mass**2)
                 wt   = pt*GetMtLead(1,mt)
		 call HFILL (pionWorkingId, pt, 0.,wt)
              enddo	   
           else
              print *,'GeGetPt> use default pt-spectra for p-beam'
	      do ibin = 1, 1000
	         pt   = 4./1000.*(float(ibin)-.5)
                 mass = particleList(pionIndex)%mass
                 mt   = sqrt(pt**2+mass**2)
                 wt   = pt*GetMtProton(1,mt)
                 call HFILL (pionWorkingId, pt, 0.,wt)
	      enddo	   
           endif
	else
	   !
	   !  Pion pt distribution exist
	   !

	   pionWorkingId = particleList(pionIndex)%ptHbookId
	endif
        particleList(pionIndex)%ptHbookId = pionWorkingId
	call HCOPY(pionWorkingId, pionPtOriginalId,' ')
	
	!
	!  Now all the other pt distributions
	!
        ptid = workingID
        do iploop = 1, size(particleList)
            ! pt distribution given for this particle?
            if (particleList(iploop)%ptHbookId == 0) then
                ! get free hbook id
                ptid = GeFreeHbookId(ptid)
                ! fill pt input distribution
                call HBOOK1 (ptid, 'pt '//particleList(iploop)%name,    &
                             1000, 0.005, 4.005, 0.)
               
                mass = particleList(iploop)%mass

                if (projectile%a.eq.208) then
                   do ibin = 1, 1000
                      pt = 4./1000.*(float(ibin)-.5)
                      mt = sqrt(pt**2+mass**2)
                      wt = pt*GetMtLead(2,mt) 
                      call HFILL (ptid, pt, 0.,wt)                           
                   enddo
                else
                   do ibin = 1, 1000
                      pt = 4./1000.*(float(ibin)-.5)
                      mt   = sqrt(pt**2+mass**2)
                      if (iploop == etaIndex) then
                          wt = pt*GetMtProton(2,mt)
                      else if  (iploop == rhoIndex .or. iploop == omegaIndex) then
                          wt = pt*GetMtProton(3,mt) 
                      else
                          wt = pt*GetMtProton(4,mt) 
                      endif
                      call HFILL (ptid, pt, 0.,wt)                           
                   enddo
                endif
                ! keep hbook id for particle
                particleList(iploop)%ptHbookId = ptid
		! copy histogram to keep original distribution before integration 
		call HCOPY(particleList(iploop)%ptHbookId, GeFreeHbookId(originalID),' ')
            endif
        enddo
        init = .false.
    endif

    ! pt distribution from histogram
    if (is_pt_weight /= 1) then             ! cumulative distribution
        pt     = HRNDM1(particleList(ipart)%ptHbookId)
        weight = 1.
    else                                    ! white distribution weighted
        call random_number(pt)
        pt     = pt*4.0
        weight = HX(particleList(ipart)%ptHbookId, pt)/                 &
                 HSUM(particleList(ipart)%ptHbookId)
    endif
    GeGetPt = pt

contains

   !
   ! mt-scaling if not pion.
   ! Use Bourquin-Gaillard for want of better things. NP B114 (1976) 334
   ! fx(pt)/fpi(pt) = ((mt(pi)+2) / (mt(x)+2) )**12.3
   !
   function GeBourquinGaillard (pt, ipart)
       implicit none
       real(kind=double)   :: GeBourquinGaillard
       real, intent(in)    :: pt
       integer, intent(in) :: ipart
       real(kind=double)   :: pimass2, partmass2

       ! get mass**2 of pion and particle
       pimass2   = particleList(GeParticleIndex("pion0"))%mass**2
       partmass2 = particleList(ipart)%mass**2
       GeBourquinGaillard = ((sqrt(pt**2+pimass2)   + 2.) /                             &
                            (sqrt(pt**2+partmass2) + 2.))** 12.3
   end function GeBourquinGaillard
   
   !
   !   Parametrization of WA80  pi0 spectrum for want of better things.
   !   central S+Au events (7.7%)
   !   Reference: Georg Hoelker, PhD thesis, University Mainz (1993)
   !
   real function GePtPi0 (pt)
       implicit none
       real, intent(in)    :: pt
       GePtPi0 = pt*(9.3/(pt+9.3))**47 
   end function GePtPi0

   !
   !  Parameterization of CERES/TAPS for 450 GeV p-beam  
   !  G.Agakichiev et al., submitted to Z.Phys C (8-10-97 )
   !
   ! implemented only for sqrt(s) = 29.1
   !
   real function GetMtProton (id,mt)
      
      integer :: id                ! particle identifier
      real    :: mt                ! transvers mass
      real    :: p1,p2,p3          ! fit parameters depending on id
     
      if (id == 1) then            ! for pions
         p1 = 1.
         p2 = 7.5
         p3 = 0.02
      else if (id == 2) then       ! for eta
         p1 = 0.15
         p2 = 6.5
         p3 = 0.0111
      else if (id == 3) then       ! for rho,omega
         p1 = 0.15
         p2 = 6.5
         p3 = 0.02
      else                         ! for phi,eta'
         p1 = 0.15
         p2 = 6.5
         p3 = 0.0111
      endif     
   
      GetMtProton = p1*exp(-mt*p2) + p3*(1-2*mt/29.1)**7.9/(1+mt**2)**4

   end function GetMtProton    

   !
   ! Parameterization for lead-beam based on NA44,NA49 and WA98 data
   ! 
   !  assume mt-scaling for all particles but treat pions different
   !
   real function GetMtLead(id,mt)
     
      implicit none
    
      real            :: mt  ! transvers mass
      integer         :: id  ! particle identifier
      real, parameter :: t1 = 0.19   ! universal slope for all particles (for NA49 Kaons)
      real, parameter :: t2 = 0.073  ! low mt component for pions 
      real, parameter :: a  = .3707  ! relative weight of both components for pions 

      if (id == 1) then
         GetMtLead =  a*exp(-mt/t1)+exp(-mt/t2)
      else
         GetMtLead =  exp(-mt/t1)
      endif

  end function GetMtLead

end function GeGetPt