!
!  GeGetRapidity
!  returns rapidity (not pseudorapidity) according to a distribution
!  given using a scaling (sigma_x = sigma_pi0 * ymax_x/ymax_pi0)
!
function GeGetRapidity(ipart,bpar,pt)
    use genesis_m, self=>GeGetRapidity
    use f90_kind
    implicit none
    integer, intent(in)          :: ipart
    real(kind=double)            :: GeGetRapidity,pt
    real, parameter 		 :: anuc = 0.938, mass_pi0=0.1349743
    real, parameter 		 :: hbymin = -2., hbymax = 8.
    integer, parameter 		 :: nhby = 1000
    real 			 :: ebeam,sqrts,y_0,gamma_proj,bpar
    real			 :: sigma_landau,emax_pi0,plmax_pi0,ymax_pi0
    real 			 :: emax,plmax,ymax,ymin,sig,ymean
    real 			 :: y, y_weight,ybin,aproj, atarg
    integer 			 :: iploop, yId,ibin
    logical, save 		 :: init = .true.

    if (init) then 
       yId = workingID
       ebeam = projectile%p%p%pz
       sqrts = sqrt( 2.*ebeam*anuc+2.*anuc**2 )
       y_0 = 0.5 * log( 2.*ebeam/anuc )
       gamma_proj = 0.5*sqrts/anuc
       sigma_landau = sqrt( log( gamma_proj ) )
       emax_pi0 = (sqrts**2.-4.*anuc**2+mass_pi0**2)/(2.*sqrts)
       plmax_pi0 = sqrt( emax_pi0**2 - mass_pi0**2 )
       ymax_pi0 = log ( (emax_pi0+plmax_pi0) / mass_pi0 )
       ymax_pi0 = log ( sqrts / mass_pi0 )

       do iploop = 1, size(particleList)
         if (particleList(iploop)%yHbookId == 0) then
           emax = (sqrts**2-4.*anuc**2+particleList(iploop)%mass**2)/2./sqrts
           plmax = sqrt( emax**2 - particleList(iploop)%mass**2 )
	   if (particleList(iploop)%mass.gt.0) then
              ymax = log ( (emax+plmax) / particleList(iploop)%mass )
              ymax = log ( sqrts / particleList(iploop)%mass )
	   else 
	      Print *,'GeGetRapidity: Using pseudorapidity instead'
	   endif
           ymin = - ymax
        
           sig = sigma_landau * ymax / ymax_pi0

!     Move to lab system
           ymax = ymax + y_0
           ymin = ymin + y_0
!     Shift of the rapidity distribution
           aproj = projectile%a
           atarg = target%a

           ymean = y_0 - y_shift (ebeam, aproj, atarg,bpar)
!     No y shift for J/psi and psi'
           if (particleList(iploop)%id.eq.441.or. &
               particleList(iploop)%id.eq.442) then
               ymean = y_0
           endif
           yId =  GeFreeHbookId(yId)
           particleList(iploop)%yHbookId = yId
           call hbook1(yId,'Rapidity '//particleList(iploop)%name,nhby,hbymin,hbymax,0.)
           ybin = (hbymax-hbymin)/float(nhby)
           y = hbymin - 0.5*ybin
           do ibin = 1,nhby
              y = y + ybin
              y_weight = dsigdy (y_0, ymean, ymax, ymin, sig, y)
              call hf1 (yId,y,y_weight)
           enddo
           call HCOPY(particleList(iploop)%yHbookId, GeFreeHbookId(originalID),' ')
         endif
       enddo
       init = .false.
     endif
      
   GeGetRapidity = hrndm1(particleList(ipart)%yHbookId)

   contains

!
!------------------------------------------------------------------------
!
     real function y_shift (ebeam, aproj, atarg,bpar)
!
     implicit none 
     real, parameter :: dens=0.17, anuc=0.938
     real	     :: ebeam, aproj,atarg,bpar
     real	     :: rmin, rmax,hmax,rcyl,vmax,hmic,vmin,rapb
     real 	     :: t_part, p_part,rapc,epsilon,sigma,a
     real, parameter :: pigrec = 3.141592653589793238

     y_shift = 0.

      if (aproj.lt.2.) then
         if (atarg.lt.12.) return
         y_shift = 1.
         return
      else if (abs(atarg/aproj-1.).lt..1) then
         return
      endif

      rmin=(3./4.*aproj/pigrec/dens)**(1./3.)
      rmax=(3./4.*atarg/pigrec/dens)**(1./3.)
      if (bpar.gt.(rmax+rmin)) return

      hmax = 2.*sqrt(rmax**2.-rmin**2.)

      if (bpar.lt.(rmax-rmin)) then
         rcyl=rmin
         vmax=2.*pigrec*rcyl**2*hmax
         hmic=rmax-sqrt(rmax**2-rmin**2)
         rapb=(rmax-hmic)/rmax
         vmin=2./3.*pigrec*rmax**3*(1.+0.5*rapb**3-1.5*rapb)
         t_part=dens*(2.*vmin+vmax)
         p_part=aproj
      else
         rcyl=0.5*(rmax+rmin-bpar)
         vmax=2.*pigrec*rcyl**2*hmax
         hmic=rmax-sqrt(rmax**2-rmin**2)
         rapb=(rmax-hmic)/rmax
         vmin=2./3.*pigrec*rmax**3*(1.+0.5*rapb**3-1.5*rapb)
         rapc=(rmin-2.*rcyl)/rmin
         t_part=dens*(2./3.*pigrec*rcyl**2*hmic+vmax)
         p_part=dens*2./3.*pigrec*rmin**3*(1.+0.5*rapc**3-1.5*rapc)
      endif
      y_shift = 0.0055 * (t_part-p_part)

      end function y_shift
!
!------------------------------------------------------------------------
!
      real function dsigdy (y_0, ymean, ymax, ymin, sig, y)
!
      real   :: y_0, ymean, ymax, ymin, sig, y
      real   :: epsilon,sigma,a

      dsigdy = 0.
      if (y.lt.ymin.or.y.gt.ymax) return

      epsilon = (y_0-ymean)**2/(ymax-ymin)
      if (y.lt.ymean) then
         sigma = sig * (1.-epsilon)
      else
         sigma = sig * (1.+epsilon)
      endif

      a = 0.75 / sigma
      dsigdy = 1. / (exp(a*(y-ymean)) + exp(-a*(y-ymean)))**2

      end function dsigdy

end function GeGetRapidity