mod_Kinematics.F90

Go to the documentation of this file.

MODULE modkinematics
implicit none
save
 
public :: getbwpropagator,reweightbwpropagator,setrunningscales,setpdfs!,EvalAlphaS
 
CONTAINS
 
 
SUBROUTINE shiftmass(p1,p2,m1,m2,p1hat,p2hat,MassWeight)
use modmisc
implicit none
real(8),intent(in) :: p1(1:4),p2(1:4)
real(8) :: m1,m2,p1hat(1:4),p2hat(1:4)
real(8),optional :: MassWeight
real(8) :: xi,eta,a,b,c,p1sq,p2sq,p1p2
real(8) :: p1hatsq, p2hatsq, p12hatsq
 
  p1sq = p1(1:4).dot.p1(1:4)
  p2sq = p2(1:4).dot.p2(1:4)
  p1p2 = p1(1:4).dot.p2(1:4)
 
  a = ( p1sq*p2(1:4) - p2sq*p1(1:4) + p1p2*(p2(1:4)-p1(1:4)) ).dot.( p1sq*p2(1:4) - p2sq*p1(1:4) + p1p2*(p2(1:4)-p1(1:4)) )
  b = ( p1sq+p2sq+2d0*p1p2+m2**2-m1**2 ) * ( p1p2**2 - p1sq*p2sq )
  c = 0.25d0*( p1sq+p2sq+2d0*p1p2+m2**2-m1**2 )**2*p1sq - (p1sq+p1p2)**2*m2**2
  eta = 1d0/2d0/a * ( -b - dsqrt( dabs(b**2 -4d0*a*c) ) )
  xi = ( p1sq+p2sq+2d0*p1p2 + m2**2 - m1**2 - 2d0*eta*(p2sq+p1p2) )/2d0/( p1sq + p1p2 )
 
  p2hat(1:4) = xi*p1(1:4) + eta*p2(1:4)
  p1hat(1:4) = (1d0-xi)*p1(1:4) + (1d0-eta)*p2(1:4)
 
  if( present(massweight) ) then
     p1hatsq = p1hat(1:4).dot.p1hat(1:4)
     p2hatsq = p2hat(1:4).dot.p2hat(1:4)
     p12hatsq = (p1hat(1:4)+p2hat(1:4)).dot.(p1hat(1:4)+p2hat(1:4)) ! Should be p1p2*2d0, but better to avoid un-anticipated uses
 
     ! Below is the same as 2d0/pi/g_d(p12hatsq,p1hatsq,p2hatsq)
     massweight = (p12hatsq**2+p1hatsq**2+p2hatsq**2-2d0*(p1hatsq*p2hatsq+p1hatsq*p12hatsq+p12hatsq*p2hatsq)) ! Writing this way instead of get_MInv should avoid the issue of - vs + invariant masses
     if(massweight.ge.0d0 .and. p12hatsq.ne.0d0) then
        massweight = sqrt(massweight/(p12hatsq**2))
     else
        massweight = 0d0
     endif
  endif
 
! if( dabs( (p1hat.dot.p1hat)-m1**2 )/m1**2.gt.1d-3 ) then
!     print *, "1",p1hat.dot.p1hat , m1**2
!     print *, p1
!     print *, p1hat
!     print *, a,b,c,eta,xi,p1sq + p1p2
!     pause
! endif
! if( dabs( (p2hat.dot.p2hat)-m2**2 )/m2**2.gt.1d-3 ) then
!     print *, "2",p2hat.dot.p2hat , m2**2
!     print *, p2
!     print *, p2hat
!     print *, a,b,c,eta,xi,p1sq + p1p2
!     pause
! endif
 
 
RETURN
END SUBROUTINE
 
 
 
 
FUNCTION zlepbranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zlepbranching
 
 
  if( xrnd .le. brlept_z_ee/(100d0*percent-brlept_z_tt) ) then
      zlepbranching = elm_
  elseif(xrnd .le. (brlept_z_ee+brlept_z_mm)/(100d0*percent-brlept_z_tt) ) then
      zlepbranching = mum_
  else
      print *, "error ",xrnd
      stop
  endif
 
!print *, "checker 2",(Brlept_Z_ee)/(100d0*percent-Brlept_Z_tt)
!print *, "checker 2",(Brlept_Z_ee+Brlept_Z_mm)/(100d0*percent-Brlept_Z_tt)
 
RETURN
END FUNCTION
 
FUNCTION zlepbranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zlepbranching_flat
 
 
  if( xrnd .le. 0.5d0 ) then
      zlepbranching_flat = elm_
  else
      zlepbranching_flat = mum_
  endif
 
RETURN
END FUNCTION
 
FUNCTION zlepplustaubranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zlepplustaubranching
 
  if( xrnd .le. brlept_z_ee ) then
      zlepplustaubranching = elm_
  elseif(xrnd .le. brlept_z_ee+brlept_z_mm ) then
      zlepplustaubranching = mum_
  elseif(xrnd .le. brlept_z_ee+brlept_z_mm+brlept_z_tt ) then
      zlepplustaubranching = tam_
  else
      print *, "error ",xrnd
      stop
  endif
 
! print *, "checker 3",Brlept_Z_ee
! print *, "checker 3",Brlept_Z_ee+Brlept_Z_mm
! print *, "checker 3",Brlept_Z_ee+Brlept_Z_mm+Brlept_Z_tt
 
RETURN
END FUNCTION
 
FUNCTION zlepplustaubranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zlepplustaubranching_flat
 
  if( xrnd .le. (1d0/3d0) ) then
      zlepplustaubranching_flat = elm_
  elseif(xrnd .le. (2d0/3d0) ) then
      zlepplustaubranching_flat = mum_
  else
      zlepplustaubranching_flat = tam_
  endif
 
RETURN
END FUNCTION
 
FUNCTION znubranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: znubranching
 
  if( xrnd .le. brlept_z_nn ) then
      znubranching = nue_
  elseif(xrnd .le. brlept_z_nn+brlept_z_nn ) then
      znubranching = num_
  elseif(xrnd .le. brlept_z_nn+brlept_z_nn+brlept_z_nn ) then
      znubranching = nut_
  else
      print *, "error ",xrnd
      stop
  endif
 
!print *, "checker 4",Brlept_Z_nn
!print *, "checker 4",Brlept_Z_nn+Brlept_Z_nn
!print *, "checker 4",Brlept_Z_nn+Brlept_Z_nn+Brlept_Z_nn
 
RETURN
END FUNCTION
 
FUNCTION znubranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: znubranching_flat
 
  if( xrnd .le. (1d0/3d0) ) then
      znubranching_flat = nue_
  elseif(xrnd .le. (2d0/3d0) ) then
      znubranching_flat = num_
  else
      znubranching_flat = nut_
  endif
 
RETURN
END FUNCTION
 
FUNCTION zquabranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zquabranching_flat
real(8),parameter :: ncol=3d0
real(8),parameter :: xxxx=1d0/15d0
real(8),parameter :: yyyy=ncol*xxxx
 
  if( xrnd .le. yyyy ) then
      zquabranching_flat = up_
  elseif(xrnd .le. yyyy+yyyy) then
      zquabranching_flat = chm_
  elseif(xrnd .le. yyyy+yyyy+yyyy) then
      zquabranching_flat = dn_
  elseif(xrnd .le. yyyy+yyyy+yyyy+yyyy) then
      zquabranching_flat = str_
  elseif(xrnd .le. yyyy+yyyy+yyyy+yyyy+yyyy) then
      zquabranching_flat = bot_
  else
      print *, "error ",xrnd
      stop
  endif
 
!print *, "checker 1",Brhadr_Z_uu,Brhadr_Z_dd
!print *, "checker 1",Brhadr_Z_uu+Brhadr_Z_cc+Brhadr_Z_dd+Brhadr_Z_ss+Brhadr_Z_bb
 
RETURN
END FUNCTION
 
FUNCTION zquabranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zquabranching
real(8),parameter :: ncol=3d0
real(8),parameter :: xxxx=1d0/15d0
real(8),parameter :: yyyy=ncol*xxxx
 
  if( xrnd .le. brhadr_z_uu ) then
      zquabranching = up_
  elseif(xrnd .le. brhadr_z_uu+brhadr_z_cc) then
      zquabranching = chm_
  elseif(xrnd .le. brhadr_z_uu+brhadr_z_cc+brhadr_z_dd) then
      zquabranching = dn_
  elseif(xrnd .le. brhadr_z_uu+brhadr_z_cc+brhadr_z_dd+brhadr_z_ss) then
      zquabranching = str_
  elseif(xrnd .le. brhadr_z_uu+brhadr_z_cc+brhadr_z_dd+brhadr_z_ss+brhadr_z_bb) then
      zquabranching = bot_
  else
      print *, "error ",xrnd
      stop
  endif
 
 
RETURN
END FUNCTION
 
FUNCTION zanybranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zanybranching_flat
real(8),parameter :: ncol=3d0
real(8),parameter :: xx=1d0/21d0
real(8),parameter :: yy=ncol*xx
 
! real(8),parameter :: xx=1d0/11d0
! real(8),parameter :: yy=xx
 
 
  if( xrnd .le. yy ) then
      zanybranching_flat = up_
  elseif(xrnd .le. 2d0*yy ) then
      zanybranching_flat = chm_
  elseif(xrnd .le. 3d0*yy ) then
      zanybranching_flat = dn_
  elseif(xrnd .le. 4d0*yy ) then
      zanybranching_flat = str_
  elseif(xrnd .le. 5d0*yy ) then
      zanybranching_flat = bot_
  elseif(xrnd .le. 5d0*yy + xx ) then
      zanybranching_flat = elm_
  elseif(xrnd .le. 5d0*yy + 2d0*xx ) then
      zanybranching_flat = mum_
  elseif(xrnd .le. 5d0*yy + 3d0*xx ) then
      zanybranching_flat = tam_
  elseif(xrnd .le. 5d0*yy + 4d0*xx ) then
      zanybranching_flat = nue_
  elseif(xrnd .le. 5d0*yy + 5d0*xx ) then
      zanybranching_flat = num_
  elseif(xrnd .le. 5d0*yy + 6d0*xx ) then
      zanybranching_flat = nut_
  else
      print *, "error ",xrnd
      stop
  endif
 
 
 
 
 
!   if( xRnd .le. yy*scale_alpha_Z_uu ) then
!       ZAnyBranching = Up_
!   elseif(xRnd .le. 2*yy*scale_alpha_Z_uu ) then
!       ZAnyBranching = Chm_
!   elseif(xRnd .le. 2*yy*scale_alpha_Z_uu+yy*scale_alpha_Z_dd ) then
!       ZAnyBranching = Dn_
!   elseif(xRnd .le. 2*yy*scale_alpha_Z_uu+2*yy*scale_alpha_Z_dd ) then
!       ZAnyBranching = Str_
!   elseif(xRnd .le. 2*yy*scale_alpha_Z_uu+3*yy*scale_alpha_Z_dd ) then
!       ZAnyBranching = Bot_
!   elseif(xRnd .le. yy*(2*scale_alpha_Z_uu+3*scale_alpha_Z_dd) + xx*scale_alpha_Z_ll ) then
!       ZAnyBranching = ElM_
!   elseif(xRnd .le. yy*(2*scale_alpha_Z_uu+3*scale_alpha_Z_dd) + xx*2*scale_alpha_Z_ll ) then
!       ZAnyBranching = MuM_
!   elseif(xRnd .le. yy*(2*scale_alpha_Z_uu+3*scale_alpha_Z_dd) + xx*3*scale_alpha_Z_ll ) then
!       ZAnyBranching = TaM_
!   elseif(xRnd .le. yy*(2*scale_alpha_Z_uu+3*scale_alpha_Z_dd)+xx*3*scale_alpha_Z_ll + xx*scale_alpha_Z_nn ) then
!       ZAnyBranching = NuE_
!   elseif(xRnd .le. yy*(2*scale_alpha_Z_uu+3*scale_alpha_Z_dd)+xx*3*scale_alpha_Z_ll + xx*2*scale_alpha_Z_nn ) then
!       ZAnyBranching = NuM_
!   elseif(xRnd .le. yy*(2*scale_alpha_Z_uu+3*scale_alpha_Z_dd)+xx*3*scale_alpha_Z_ll + xx*3*scale_alpha_Z_nn ) then
!       ZAnyBranching = NuT_
!   else
!       print *, "error ",xRnd
!       stop
!   endif
 
 
 
RETURN
END FUNCTION
 
FUNCTION zanybranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: zanybranching
 
 
  if( xrnd .le. br_z_uu ) then
      zanybranching = up_
  elseif(xrnd .le. br_z_uu+br_z_cc) then
      zanybranching = chm_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd) then
      zanybranching = dn_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss) then
      zanybranching = str_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb) then
      zanybranching = bot_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb+br_z_ee) then
      zanybranching = elm_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb+br_z_ee+br_z_mm) then
      zanybranching = mum_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb+br_z_ee+br_z_mm+br_z_tt) then
      zanybranching = tam_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb+br_z_ee+br_z_mm+br_z_tt+br_z_nn) then
      zanybranching = nue_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb+br_z_ee+br_z_mm+br_z_tt+br_z_nn+br_z_nn) then
      zanybranching = num_
  elseif(xrnd .le. br_z_uu+br_z_cc+br_z_dd+br_z_ss+br_z_bb+br_z_ee+br_z_mm+br_z_tt+br_z_nn+br_z_nn+br_z_nn) then
      zanybranching = nut_
  else
      print *, "error ",xrnd
      stop
  endif
 
 
!   if( xRnd .le. scale_alpha_Z_uu*Br_Z_uu ) then
!       ZAnyBranching = Up_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc) then
!       ZAnyBranching = Chm_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd) then
!       ZAnyBranching = Dn_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss) then
!       ZAnyBranching = Str_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb) then
!       ZAnyBranching = Bot_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee) then
!       ZAnyBranching = ElM_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee+scale_alpha_Z_ll*Br_Z_mm) then
!       ZAnyBranching = MuM_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee+scale_alpha_Z_ll*Br_Z_mm+scale_alpha_Z_ll*Br_Z_tt) then
!       ZAnyBranching = TaM_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee+scale_alpha_Z_ll*Br_Z_mm+scale_alpha_Z_ll*Br_Z_tt+scale_alpha_Z_nn*Br_Z_nn) then
!       ZAnyBranching = NuE_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee+scale_alpha_Z_ll*Br_Z_mm+scale_alpha_Z_ll*Br_Z_tt+scale_alpha_Z_nn*Br_Z_nn+scale_alpha_Z_nn*Br_Z_nn) then
!       ZAnyBranching = NuM_
!   elseif(xRnd .le. scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee+scale_alpha_Z_ll*Br_Z_mm+scale_alpha_Z_ll*Br_Z_tt+scale_alpha_Z_nn*Br_Z_nn+scale_alpha_Z_nn*Br_Z_nn+scale_alpha_Z_nn*Br_Z_nn) then
!       ZAnyBranching = NuT_
!   else
!       print *, "error ",xRnd
!       stop
!   endif
 
 
! print *, "checker ",Br_Z_uu+Br_Z_cc+Br_Z_dd+Br_Z_ss+Br_Z_bb+Br_Z_ee+Br_Z_mm+Br_Z_tt+Br_Z_nn+Br_Z_nn+Br_Z_nn
! print *, "checker ",scale_alpha_Z_uu*Br_Z_uu+scale_alpha_Z_uu*Br_Z_cc+scale_alpha_Z_dd*Br_Z_dd+scale_alpha_Z_dd*Br_Z_ss+scale_alpha_Z_dd*Br_Z_bb+scale_alpha_Z_ll*Br_Z_ee+scale_alpha_Z_ll*Br_Z_mm+scale_alpha_Z_ll*Br_Z_tt+scale_alpha_Z_nn*Br_Z_nn+scale_alpha_Z_nn*Br_Z_nn+scale_alpha_Z_nn*Br_Z_nn
! pause
 
 
RETURN
END FUNCTION
 
FUNCTION wlepbranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wlepbranching
 
  if( xrnd .le. brlept_w_en /(100d0*percent-brlept_w_tn) ) then
      wlepbranching = elm_
  elseif(xrnd .le. (brlept_w_en+brlept_w_mn)/(100d0*percent-brlept_w_tn) ) then
      wlepbranching = mum_
  else
      print *, "error ",xrnd
      stop
  endif
 
!print *, "checker 6",Brlept_W_en /(100d0*percent-Brlept_W_tn)
!print *, "checker 6",(Brlept_W_en+Brlept_W_mn)/(100d0*percent-Brlept_W_tn)
 
RETURN
END FUNCTION
 
FUNCTION wlepbranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wlepbranching_flat
 
  if( xrnd .le. 0.5d0 ) then
      wlepbranching_flat = elm_
  else
      wlepbranching_flat = mum_
  endif
 
RETURN
END FUNCTION
 
FUNCTION wlepplustaubranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wlepplustaubranching
 
  if( xrnd .le. brlept_w_en ) then
      wlepplustaubranching = elm_
  elseif(xrnd .le. brlept_w_en+brlept_w_mn ) then
      wlepplustaubranching = mum_
  elseif(xrnd .le. brlept_w_en+brlept_w_mn+brlept_w_tn ) then
      wlepplustaubranching = tam_
  else
      print *, "error ",xrnd
      stop
  endif
 
! print *, "checker 7",Brlept_W_en
! print *, "checker 7",Brlept_W_en+Brlept_W_mn
! print *, "checker 7",Brlept_W_en+Brlept_W_mn+Brlept_W_tn
 
RETURN
END FUNCTION
 
FUNCTION wlepplustaubranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wlepplustaubranching_flat
 
  if( xrnd .le. (1d0/3d0) ) then
      wlepplustaubranching_flat = elm_
  elseif(xrnd .le. (2d0/3d0) ) then
      wlepplustaubranching_flat = mum_
  else
      wlepplustaubranching_flat = tam_
  endif
 
RETURN
END FUNCTION
 
FUNCTION wquaupbranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wquaupbranching
 
  if( xrnd .le. brhadr_w_ud ) then
      wquaupbranching = up_
  elseif(xrnd .le. brhadr_w_ud+brhadr_w_cs ) then
      wquaupbranching = chm_
  else
      print *, "error ",xrnd
      stop
  endif
 
 
!print *, "checker 8",Brhadr_W_ud
!print *, "checker 8",Brhadr_W_ud+Brhadr_W_cs
 
RETURN
END FUNCTION
 
FUNCTION wquaupbranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wquaupbranching_flat
 
  if( xrnd .le. 0.5d0 ) then
      wquaupbranching_flat = up_
  else
      wquaupbranching_flat = chm_
  endif
 
RETURN
END FUNCTION
 
FUNCTION wanybranching_flat(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wanybranching_flat
real(8),parameter :: ncol=3d0
real(8),parameter :: xx=1d0/9d0
real(8),parameter :: yy=ncol*xx
 
 
  if( xrnd .le.        yy ) then
      wanybranching_flat = up_
  elseif(xrnd .le. 2d0*yy ) then
      wanybranching_flat = chm_
  elseif(xrnd .le. 2d0*yy +     xx ) then
      wanybranching_flat = elm_
  elseif(xrnd .le. 2d0*yy + 2d0*xx ) then
      wanybranching_flat = mum_
  elseif(xrnd .le. 2d0*yy + 3d0*xx ) then
      wanybranching_flat = tam_
  else
      print *, "error ",xrnd
      stop
  endif
 
RETURN
END FUNCTION
 
FUNCTION wanybranching(xRnd)
use modparameters
implicit none
real(8) :: xrnd
integer :: wanybranching
 
 
  if( xrnd .le. br_w_ud ) then
      wanybranching = up_
  elseif(xrnd .le. br_w_ud+br_w_cs ) then
      wanybranching = chm_
  elseif(xrnd .le. br_w_ud+br_w_cs+br_w_en ) then
      wanybranching = elm_
  elseif(xrnd .le. br_w_ud+br_w_cs+br_w_en+br_w_mn ) then
      wanybranching = mum_
  elseif(xrnd .le. br_w_ud+br_w_cs+br_w_en+br_w_mn+br_w_tn ) then
      wanybranching = tam_
  else
      print *, "error ",xrnd
      stop
  endif
 
 
!   if( xRnd .le. scale_alpha_W_ud*Br_W_ud ) then
!       WAnyBranching = Up_
!   elseif(xRnd .le. scale_alpha_W_ud*Br_W_ud+scale_alpha_W_cs*Br_W_cs ) then
!       WAnyBranching = Chm_
!   elseif(xRnd .le. scale_alpha_W_ud*Br_W_ud+scale_alpha_W_cs*Br_W_cs+scale_alpha_W_ln*Br_W_en ) then
!       WAnyBranching = ElM_
!   elseif(xRnd .le. scale_alpha_W_ud*Br_W_ud+scale_alpha_W_cs*Br_W_cs+scale_alpha_W_ln*Br_W_en+scale_alpha_W_ln*Br_W_mn ) then
!       WAnyBranching = MuM_
!   elseif(xRnd .le. scale_alpha_W_ud*Br_W_ud+scale_alpha_W_cs*Br_W_cs+scale_alpha_W_ln*Br_W_en+scale_alpha_W_ln*Br_W_mn+scale_alpha_W_ln*Br_W_tn ) then
!       WAnyBranching = TaM_
!   else
!       print *, "error ",xRnd
!       stop
!   endif
 
! print *, "checker 2",Br_W_ud+Br_W_cs+Br_W_en+Br_W_mn+Br_W_tn
! print *, "checker 2",scale_alpha_W_ud*Br_W_ud+scale_alpha_W_cs*Br_W_cs+scale_alpha_W_ln*Br_W_en+scale_alpha_W_ln*Br_W_mn+scale_alpha_W_ln*Br_W_tn
! pause
 
RETURN
END FUNCTION
 
! VBranching and VVBranchings also take into account color multiplicity, so it should be used only when color is not taken into account inside the ME!
! This is why some final states count 3 times instead of just once.
! Rule of thumb for counting final states in this function is e,mu,ta,nus=1, d,u,s,c,b=3
! When a CKM partner is called, since sum(VCKM**@)=1 along a row, counting should not be unaffected.
SUBROUTINE vbranching(DecayMode,MY_IDUP,ICOLUP,CombWeight,ColorBase)
use modparameters
use modmisc
implicit none
integer, intent(in) :: DecayMode
integer :: MY_IDUP(1:3),ICOLUP(1:2,1:2),DKFlavor,ICOLUP_Base
integer, optional ::ColorBase
real(8), intent(out) :: CombWeight
real(8) :: DKRnd
 
!    particle associations:
!
!    IDUP(6)  -->  MomDK(:,2)  -->     v-spinor
!    IDUP(7)  -->  MomDK(:,1)  -->  ubar-spinor
!    IDUP(8)  -->  MomDK(:,4)  -->     v-spinor
!    IDUP(9)  -->  MomDK(:,3)  -->  ubar-spinor
 
   if (present(colorbase)) then
       icolup_base = colorbase
   else
       icolup_base = 800
   endif
 
   icolup(:,:) = 0
   combweight = 1d0
   if( decaymode.eq.0 ) then! Z1->2l
        call random_number(dkrnd)
        my_idup(1) = z0_
        dkflavor = zlepbranching_flat( dkrnd )!= ElM or MuM
        my_idup(2) =-dkflavor
        my_idup(3) =+dkflavor
        combweight = combweight*2d0
   elseif( decaymode.eq.1 ) then! Z1->2q
        call random_number(dkrnd)
        my_idup(1) = z0_
        dkflavor = zquabranching_flat( dkrnd )!= Up,Dn,Chm,Str,Bot
        my_idup(2) =-dkflavor
        my_idup(3) =+dkflavor
        icolup(1:2,1) = (/            0,icolup_base+3/)
        icolup(1:2,2) = (/icolup_base+3,            0/)
        combweight = combweight*5d0*3d0
   elseif( decaymode.eq.2 ) then! Z1->2tau
        my_idup(1) = z0_
        my_idup(2) = tap_
        my_idup(3) = tam_
   elseif( decaymode.eq.3 ) then! Z1->2nu
        call random_number(dkrnd)
        my_idup(1) = z0_
        dkflavor = znubranching_flat( dkrnd )!= NuE,NuM,NuT
        my_idup(2) =-dkflavor
        my_idup(3) =+dkflavor
        combweight = combweight*3d0
   elseif( decaymode.eq.4 ) then! W1(+)->lnu
        call random_number(dkrnd)
        my_idup(1) = wp_
        dkflavor = wlepbranching_flat( dkrnd )!= ElM or MuM
        my_idup(2) = +abs(dkflavor)     ! lepton(+)
        my_idup(3) = +abs(dkflavor)+7   ! neutrino
        combweight = combweight*2d0
   elseif( decaymode.eq.5 ) then! W1(+)->2q
        call random_number(dkrnd)
        my_idup(1) = wp_
        dkflavor = wquaupbranching_flat( dkrnd )!= Up,Chm
        my_idup(3) = +abs(dkflavor)           ! up flavor
        my_idup(2) = getckmpartner(my_idup(3))! anti-dn flavor
        icolup(1:2,1) = (/            0,icolup_base+3/)
        icolup(1:2,2) = (/icolup_base+3,            0/)
        combweight = combweight*2d0*3d0
   elseif( decaymode.eq.6 ) then! W1(+)->taunu
        my_idup(1) = wp_
        my_idup(2) = tap_
        my_idup(3) = nut_
   elseif( decaymode.eq.7 ) then! photon
        my_idup(1) = pho_
        my_idup(2) = not_a_particle_
        my_idup(3) = not_a_particle_
   elseif( decaymode.eq.8 ) then! Z1->2l+2tau
        call random_number(dkrnd)
        my_idup(1) = z0_
        dkflavor = zlepplustaubranching_flat( dkrnd )!= ElM or MuM or TaM
        my_idup(2) =-dkflavor
        my_idup(3) =+dkflavor
        combweight = combweight*3d0
   elseif( decaymode.eq.9 ) then! Z1-> anything
        call random_number(dkrnd)
        my_idup(1) = z0_
        dkflavor = zanybranching_flat( dkrnd )
        my_idup(2) =-dkflavor
        my_idup(3) =+dkflavor
        if(isaquark(dkflavor)) then
           icolup(1:2,1) = (/            0,icolup_base+3/)
           icolup(1:2,2) = (/icolup_base+3,            0/)
        endif
        combweight = combweight*21d0!*(6d0 + 5d0*3d0)
   elseif( decaymode.eq.10 ) then! W1(+)->l+tau  +nu
        call random_number(dkrnd)
        my_idup(1) = wp_
        dkflavor = wlepplustaubranching_flat( dkrnd )!= ElM or MuM or TaM
        my_idup(2) = +abs(dkflavor)     ! lepton(+)
        my_idup(3) = +abs(dkflavor)+7   ! neutrino
        combweight = combweight*3d0
   elseif( decaymode.eq.11 ) then! W1(+)-> anything
        call random_number(dkrnd)
        my_idup(1) = wp_
        dkflavor = wanybranching_flat( dkrnd )
        if(isaquark(dkflavor)) then
           my_idup(3) = +abs(dkflavor)           ! up flavor
           my_idup(2) = getckmpartner(my_idup(3))! anti-dn flavor
           icolup(1:2,1) = (/            0,icolup_base+3/)
           icolup(1:2,2) = (/icolup_base+3,            0/)
        else
           my_idup(2) = +abs(dkflavor)     ! lepton(+)
           my_idup(3) = +abs(dkflavor)+7   ! neutrino
        endif
        combweight = combweight*9d0!*(3d0 + 2d0*3d0)
   ! Exclusive Z1 decay modes
   elseif( &
      decaymode.eq.-2*2   .or. & ! Z->ee
      decaymode.eq.-3*3   .or. & ! Z->mumu
      decaymode.eq.-5*5   .or. & ! Z->dd
      decaymode.eq.-7*7   .or. & ! Z->uu
      decaymode.eq.-11*11 .or. & ! Z->ss
      decaymode.eq.-13*13 .or. & ! Z->cc
      decaymode.eq.-17*17      & ! Z->bb
      ) then
        my_idup(1) = z0_
        dkflavor = (                      &
            elm_*logicaltointeger(decaymode.eq.-2*2)   + & ! Z->ee
            mum_*logicaltointeger(decaymode.eq.-3*3)   + & ! Z->mumu
            dn_*logicaltointeger(decaymode.eq.-5*5)    + & ! Z->dd
            up_*logicaltointeger(decaymode.eq.-7*7)    + & ! Z->uu
            str_*logicaltointeger(decaymode.eq.-11*11) + & ! Z->ss
            chm_*logicaltointeger(decaymode.eq.-13*13) + & ! Z->cc
            bot_*logicaltointeger(decaymode.eq.-17*17)   & ! Z->bb
        )
        my_idup(2) =-dkflavor
        my_idup(3) =+dkflavor
        if(isaquark(dkflavor)) then
           icolup(1:2,1) = (/            0,icolup_base+3/)
           icolup(1:2,2) = (/icolup_base+3,            0/)
        endif
   ! Exclusive W+ decay modes
   elseif( &
      decaymode.eq.-2*1   .or. & ! W->enu
      decaymode.eq.-3*1   .or. & ! W->munu
      decaymode.eq.-5*7   .or. & ! W->du
      decaymode.eq.-5*13  .or. & ! W->dc
      decaymode.eq.-11*7  .or. & ! W->su
      decaymode.eq.-11*13 .or. & ! W->sc
      decaymode.eq.-17*7  .or. & ! W->bu
      decaymode.eq.-17*13      & ! W->bc
      ) then
        my_idup(1) = wp_
        my_idup(2) = (                    &
           elp_*logicaltointeger(decaymode.eq.-2*1)    + & ! W->enu
           mup_*logicaltointeger(decaymode.eq.-3*1)    + & ! W->munu
           adn_*logicaltointeger(decaymode.eq.-5*7)    + & ! W->du
           adn_*logicaltointeger(decaymode.eq.-5*13)   + & ! W->dc
           astr_*logicaltointeger(decaymode.eq.-11*7)  + & ! W->su
           astr_*logicaltointeger(decaymode.eq.-11*13) + & ! W->sc
           abot_*logicaltointeger(decaymode.eq.-17*7)  + & ! W->bu
           abot_*logicaltointeger(decaymode.eq.-17*13)   & ! W->bc
        )
        my_idup(3) = (                   &
           nue_*logicaltointeger(decaymode.eq.-2*1)   + & ! W->enu
           num_*logicaltointeger(decaymode.eq.-3*1)   + & ! W->munu
           up_*logicaltointeger(decaymode.eq.-5*7)    + & ! W->du
           chm_*logicaltointeger(decaymode.eq.-5*13)  + & ! W->dc
           up_*logicaltointeger(decaymode.eq.-11*7)   + & ! W->su
           chm_*logicaltointeger(decaymode.eq.-11*13) + & ! W->sc
           up_*logicaltointeger(decaymode.eq.-17*7)   + & ! W->bu
           chm_*logicaltointeger(decaymode.eq.-17*13)   & ! W->bc
        )
        if(isaquark(dkflavor)) then
           icolup(1:2,1) = (/            0,icolup_base+3/)
           icolup(1:2,2) = (/icolup_base+3,            0/)
        endif
   endif
 
 
RETURN
END SUBROUTINE
 
SUBROUTINE vvbranchings(MY_IDUP,ICOLUP,CombWeight,ColorBase)
use modparameters
use modmisc
implicit none
integer :: MY_IDUP(4:9),ICOLUP(1:2,6:9),ICOLUP_Base
integer, optional ::ColorBase
real(8), intent(out) :: CombWeight
integer :: tmp_idup(1:3),tmp_icolup(1:2,1:2)
real(8) :: tmp_CombWeight
real(8) :: DKRnd
 
!    particle associations:
!
!    IDUP(6)  -->  MomDK(:,2)  -->     v-spinor
!    IDUP(7)  -->  MomDK(:,1)  -->  ubar-spinor
!    IDUP(8)  -->  MomDK(:,4)  -->     v-spinor
!    IDUP(9)  -->  MomDK(:,3)  -->  ubar-spinor
 
   if (present(colorbase)) then
       icolup_base = colorbase
   else
       icolup_base = 800
   endif
 
   icolup(:,:) = 0
   combweight = 1d0
 
   call vbranching(decaymode1, tmp_idup, tmp_icolup, tmp_combweight, icolup_base)
   my_idup(4) = tmp_idup(1)
   my_idup(6) = tmp_idup(2)
   my_idup(7) = tmp_idup(3)
   icolup(1:2,6) = tmp_icolup(1:2,1)
   icolup(1:2,7) = tmp_icolup(1:2,2)
   combweight = combweight * tmp_combweight
 
 
   icolup_base = icolup_base+1
   call vbranching(decaymode2, tmp_idup, tmp_icolup, tmp_combweight, icolup_base)
   my_idup(5) = tmp_idup(1)
   my_idup(8) = tmp_idup(2)
   my_idup(9) = tmp_idup(3)
   icolup(1:2,8) = tmp_icolup(1:2,1)
   icolup(1:2,9) = tmp_icolup(1:2,2)
   if (isawdecay(decaymode2)) then
      my_idup(5) = -my_idup(5)
      call swap(my_idup(8),my_idup(9))
      my_idup(8) = -my_idup(8)
      my_idup(9) = -my_idup(9)
   endif
   if( process.lt.110 .or. process.gt.114) then ! for tHq processes only one V branching!
     combweight = combweight * tmp_combweight
   endif
 
RETURN
END SUBROUTINE
 
FUNCTION getckmpartner( Flavor )
use modmisc
use modparameters
implicit none
integer :: flavor,getckmpartner
real(8) :: flavorrnd,sumckm,vsq(1:3)
 
    call random_number(flavorrnd)
 
 
    if( abs(flavor).eq.abs(up_) ) then
        vsq(1) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(dn_)) ))**2
        vsq(2) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(str_)) ))**2
        vsq(3) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(bot_)) ))**2
        vsq(:) = vsq(:)/scale_alpha_w_ud
 
        sumckm = vsq(1)+vsq(2)+vsq(3)
        flavorrnd = flavorrnd*sumckm
 
        if( flavorrnd.le.vsq(1) ) then!  u-->d
           getckmpartner = -sign(1,flavor) * abs(dn_)
        elseif( flavorrnd.le.(vsq(2)+vsq(1)) ) then!  u-->s
           getckmpartner = -sign(1,flavor) * abs(str_)
        else!  u-->b
           getckmpartner = -sign(1,flavor) * abs(bot_)
        endif
 
    elseif( abs(flavor).eq.abs(chm_) ) then
        vsq(1) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(dn_)) ))**2
        vsq(2) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(str_)) ))**2
        vsq(3) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(bot_)) ))**2
        vsq(:) = vsq(:)/scale_alpha_w_cs
 
        sumckm = vsq(1)+vsq(2)+vsq(3)
        flavorrnd = flavorrnd*sumckm
 
        if( flavorrnd.le.vsq(2) ) then!  c-->s
           getckmpartner = -sign(1,flavor) * abs(str_)
        elseif( flavorrnd.le.(vsq(1)+vsq(2)) ) then!  c-->d
           getckmpartner = -sign(1,flavor) * abs(dn_)
        else!  c-->b
           getckmpartner = -sign(1,flavor) * abs(bot_)
        endif
 
    elseif( abs(flavor).eq.abs(top_) ) then
        vsq(1) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(dn_)) ))**2
        vsq(2) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(str_)) ))**2
        vsq(3) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(bot_)) ))**2
 
        sumckm = vsq(1)+vsq(2)+vsq(3)
        flavorrnd = flavorrnd*sumckm
 
        if( flavorrnd.le.vsq(3) ) then!  t-->b
           getckmpartner = -sign(1,flavor) * abs(bot_)
        elseif( flavorrnd.le.(vsq(2)+vsq(3)) ) then!  t-->s
           getckmpartner = -sign(1,flavor) * abs(str_)
        else!  t-->d
           getckmpartner = -sign(1,flavor) * abs(dn_)
        endif
 
 
    elseif( abs(flavor).eq.abs(dn_) ) then
        vsq(1) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(up_)) ))**2
        vsq(2) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(chm_)) ))**2
        vsq(3) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(top_)) ))**2
        vsq(:) = vsq(:)/scale_alpha_w_ud
 
        sumckm = vsq(1)+vsq(2)+vsq(3)
        flavorrnd = flavorrnd*sumckm
 
        if( flavorrnd.le.vsq(1) ) then!  d-->u
           getckmpartner = -sign(1,flavor) * abs(up_)
        elseif( flavorrnd.le.(vsq(2)+vsq(1)) ) then!  d-->c
           getckmpartner = -sign(1,flavor) * abs(chm_)
        else!  d-->t
           getckmpartner = -sign(1,flavor) * abs(top_)
        endif
 
    elseif( abs(flavor).eq.abs(str_) ) then
        vsq(1) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(up_)) ))**2
        vsq(2) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(chm_)) ))**2
        vsq(3) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(top_)) ))**2
        vsq(:) = vsq(:)/scale_alpha_w_cs
 
        sumckm = vsq(1)+vsq(2)+vsq(3)
        flavorrnd = flavorrnd*sumckm
 
        if( flavorrnd.le.vsq(2) ) then!  s-->c
           getckmpartner = -sign(1,flavor) * abs(chm_)
        elseif( flavorrnd.le.(vsq(1)+vsq(2)) ) then!  s-->u
           getckmpartner = -sign(1,flavor) * abs(up_)
        else!  s-->t
           getckmpartner = -sign(1,flavor) * abs(top_)
        endif
 
    elseif( abs(flavor).eq.abs(bot_) ) then
        vsq(1) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(up_)) ))**2
        vsq(2) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(chm_)) ))**2
        vsq(3) = (ckm( convertlhereverse(abs(flavor)), convertlhereverse(abs(top_)) ))**2
 
        sumckm = vsq(1)+vsq(2)+vsq(3)
        flavorrnd = flavorrnd*sumckm
 
        if( flavorrnd.le.vsq(3) ) then!  b-->t
           getckmpartner = -sign(1,flavor) * abs(top_)
        elseif( flavorrnd.le.(vsq(2)+vsq(3)) ) then!  b-->c
           getckmpartner = -sign(1,flavor) * abs(chm_)
        else!  b -->u
           getckmpartner = -sign(1,flavor) * abs(up_)
        endif
 
    else
        call error("Dn flavor conversion not yet implemented")
    endif
 
RETURN
END FUNCTION
 
FUNCTION getckmpartner_flat( Flavor )
use modmisc
use modparameters
implicit none
integer :: flavor,getckmpartner_flat
real(8) :: flavorrnd,sumckm,vsq(1:3)
 
    call random_number(flavorrnd)
    vsq(:) = 1d0
    sumckm = vsq(1)+vsq(2)+vsq(3)
    flavorrnd = flavorrnd*sumckm
 
    if( abs(flavor).eq.abs(up_) ) then
        if( flavorrnd.le.vsq(1) ) then!  u-->d
           getckmpartner_flat = -sign(1,flavor) * abs(dn_)
        elseif( flavorrnd.le.(vsq(2)+vsq(1)) ) then!  u-->s
           getckmpartner_flat = -sign(1,flavor) * abs(str_)
        else!  u-->b
           getckmpartner_flat = -sign(1,flavor) * abs(bot_)
        endif
 
    elseif( abs(flavor).eq.abs(chm_) ) then
        if( flavorrnd.le.vsq(2) ) then!  c-->s
           getckmpartner_flat = -sign(1,flavor) * abs(str_)
        elseif( flavorrnd.le.(vsq(1)+vsq(2)) ) then!  c-->d
           getckmpartner_flat = -sign(1,flavor) * abs(dn_)
        else!  c-->b
           getckmpartner_flat = -sign(1,flavor) * abs(bot_)
        endif
 
    elseif( abs(flavor).eq.abs(top_) ) then
        if( flavorrnd.le.vsq(3) ) then!  t-->b
           getckmpartner_flat = -sign(1,flavor) * abs(bot_)
        elseif( flavorrnd.le.(vsq(2)+vsq(3)) ) then!  t-->s
           getckmpartner_flat = -sign(1,flavor) * abs(str_)
        else!  t-->d
           getckmpartner_flat = -sign(1,flavor) * abs(dn_)
        endif
 
 
    elseif( abs(flavor).eq.abs(dn_) ) then
        if( flavorrnd.le.vsq(1) ) then!  d-->u
           getckmpartner_flat = -sign(1,flavor) * abs(up_)
        elseif( flavorrnd.le.(vsq(2)+vsq(1)) ) then!  d-->c
           getckmpartner_flat = -sign(1,flavor) * abs(chm_)
        else!  d-->t
           getckmpartner_flat = -sign(1,flavor) * abs(top_)
        endif
 
    elseif( abs(flavor).eq.abs(str_) ) then
        if( flavorrnd.le.vsq(2) ) then!  s-->c
           getckmpartner_flat = -sign(1,flavor) * abs(chm_)
        elseif( flavorrnd.le.(vsq(1)+vsq(2)) ) then!  s-->u
           getckmpartner_flat = -sign(1,flavor) * abs(up_)
        else!  s-->t
           getckmpartner_flat = -sign(1,flavor) * abs(top_)
        endif
 
    elseif( abs(flavor).eq.abs(bot_) ) then
        if( flavorrnd.le.vsq(3) ) then!  b-->t
           getckmpartner_flat = -sign(1,flavor) * abs(top_)
        elseif( flavorrnd.le.(vsq(2)+vsq(3)) ) then!  b-->c
           getckmpartner_flat = -sign(1,flavor) * abs(chm_)
        else!  b -->u
           getckmpartner_flat = -sign(1,flavor) * abs(up_)
        endif
 
    else
        call error("Dn flavor conversion not yet implemented")
    endif
 
RETURN
END FUNCTION
 
 
 
FUNCTION getbwpropagator(sHat, scheme)
use modmisc
use modparameters
implicit none
real(8) :: getbwpropagator,shat
real(8) :: mhb, ghb, biggamma
integer :: scheme
double precision decaymass, qqq, qqq0
 
    if( scheme.eq.1 ) then! running width
        getbwpropagator =  1d0/( (shat-m_reso**2)**2 + (shat*ga_reso/m_reso)**2 )
    elseif( scheme.eq.2 ) then! fixed width
        getbwpropagator = 1d0/( (shat-m_reso**2)**2 + (m_reso*ga_reso)**2 )
    elseif( scheme.eq.3 ) then! Passarino's CPS
        if( mubarh.lt.0d0 .or. gabarh.lt.0d0 ) then
          call call_hto(m_reso/gev, m_top/gev, mhb, ghb)
          if( isnan(mubarh).or.isnan(gabarh) ) then
            print *, "Passarino's CALL_HTO returned a NaN"
            print *, "(gabarH,Ehat)",gabarh,dsqrt(dabs(shat))/gev
            stop 1
            RETURN
          endif
          mhb = mhb*gev
          ghb = ghb*gev
 
          mubarh = sqrt(mhb**2/(1d0+(ghb/mhb)**2))
          gabarh = mubarh/mhb*ghb
        endif
 
        getbwpropagator = 1d0/( (shat-mubarh**2)**2 + (mubarh*gabarh)**2 )
 
        !call HTO_gridHt(dsqrt(dabs(sHat))/GeV,BigGamma)
        !BigGamma = BigGamma*GeV
 
        !print *, dsqrt(dabs(sHat))/GeV, gabarH/GeV, BigGamma/GeV
    elseif( scheme.eq.4) then !alternate running width
        if ( m_zprime.ne.-1 ) then !if you're using ZPrime then it will consider the ZPrime mass
            decaymass = m_zprime
        else
            decaymass = m_z !if you're doing Z-ZPrime stuff this scheme won't work anyways so the point is moot
        end if
        ! PRINT *, decayMass, decayMass**2, sHat, 0.25d0*sHat, M_Reso, 0.25d0*M_Reso**2
        if(0.25d0*(shat) < decaymass**2) then !sHat is the mass squared!
            qqq = 0
        else
            qqq = sqrt(0.25d0*(shat) - decaymass**2)
        endif
        qqq0=sqrt(0.25d0*(m_reso**2) - decaymass**2)
        getbwpropagator =  1d0/( (shat-m_reso**2)**2 + (m_reso*ga_reso*qqq/qqq0)**2 ) !new style running width
    elseif( scheme.eq.0 ) then  !remove the propagator completely
        getbwpropagator = 1d0
    else
        print *, "Invalid scheme: ", scheme
        stop 1
    endif
 
 
RETURN
END FUNCTION
 
FUNCTION reweightbwpropagator(sHat)! sHat is the resonance inv. mass squared
use modmisc
use modparameters
implicit none
real(8) :: reweightbwpropagator,shat
real(8) :: breitwigner,breitwigner_run
 
 
     reweightbwpropagator = 1d0
     breitwigner = getbwpropagator(shat, 2)
     breitwigner_run = getbwpropagator(shat, widthscheme)
 
    reweightbwpropagator = breitwigner_run/breitwigner
 
 
RETURN
END FUNCTION
 
 
 
 
SUBROUTINE boost2lab(x1,x2,NumPart,Mom)
implicit none
real(8) Mom(1:4,1:NumPart)
real(8) x1,x2
real(8) gamma,betagamma,MomTmp1,MomTmp4
integer :: i,NumPart
 
  gamma     = (x1+x2)/2d0/dsqrt(x1*x2)
  betagamma = (x2-x1)/2d0/dsqrt(x1*x2)
 
  do i=1,numpart
      momtmp1=mom(1,i)
      momtmp4=mom(4,i)
      mom(1,i)= gamma*momtmp1 - betagamma*momtmp4
      mom(4,i)= gamma*momtmp4 - betagamma*momtmp1
  enddo
 
RETURN
END SUBROUTINE
 
!LORENTZ.F
!VERSION 20130123
!
!A subroutine that performs a general boost to a four vector
!(vector) based on another four vector (boost). The primed and
!unprimed frames have their axes in parallel to one another.
!Rotation is not performed by this subroutine.
      subroutine lorentz(vector, boost)
 
      implicit none
 
      double precision vector(4), boost(4)
      double precision lambdaMtrx(4,4), vector_copy(4)
      double precision beta(2:4), beta_sq, gamma
      integer i,j
      double precision, parameter :: epsilon = 1d-13 !a small quantity slightly above machine precision
 
!      double precision KRONECKER_DELTA
!      external KRONECKER_DELTA
 
      do i=2,4
        beta(i) = boost(i)/boost(1)
      enddo
 
      beta_sq = beta(2)**2+beta(3)**2+beta(4)**2
 
  if(beta_sq.ge.epsilon)then
 
      gamma = 1d0/dsqrt(1d0-beta_sq)
 
      lambdamtrx(1,1) = gamma
 
      do i=2,4
        lambdamtrx(1,i) = gamma*beta(i)
        lambdamtrx(i,1) = lambdamtrx(1,i)
      enddo
 
      do i=2,4
      do j=2,4
        lambdamtrx(i,j) = (gamma-1d0)*beta(i)*beta(j)/beta_sq + kronecker_delta(i,j)
      enddo
      enddo
 
 
!apply boost to vector1
      vector_copy = vector
      vector = 0d0
      do i=1,4
      do j=1,4
        vector(i) = vector(i) + lambdamtrx(i,j)*vector_copy(j)
      enddo
      enddo
  endif
 
      return
      END subroutine lorentz
! KRONECKER_DELTA.F
!
! KRONECKER_DELTA(i,j)
! A function that returns 1 if i=j, and 0 otherwise.
      double precision function kronecker_delta(i,j)
      integer i,j
      if(i.eq.j)then
        kronecker_delta = 1d0
      else
        kronecker_delta = 0d0
      endif
 
      return
      end function kronecker_delta
 
 
 
 
 
subroutine setrunningscales(p,id) ! p in JHU-GeV, id in JHUGen conventions
use modparameters
use modmisc
implicit none
real(dp), intent(in) :: p(1:4,4:6) ! No need to run the second index from 4 to 7: pH, pJ1, pJ2
integer, intent(in) :: id(4:7) ! id_JJH/id_JJVV, id_J1, id_J2, id_JJ (if applicable)
real(8) :: polemass(3:7) ! mJJH, mH, mJ1, mJ2, mJJ (if applicable)
real(8) :: pjjhstar(4),phstar(4),pj(4,2),pjj(4),pjhstar(4),ptjet(5:6),maxptjet,minptjet
integer idx,ip
 
   phstar(:) = 0d0
   pjj(:) = 0d0
   ptjet(:) = 0d0
   polemass(3) = getmass(id(4)) ! Pole mass of the JJH system
   polemass(4) = m_reso
   do idx=4,6
      if(idx.eq.4) then
         do ip=1,4
            phstar(ip) = phstar(ip) + p(ip,idx)
         enddo
      else
         polemass(idx) = getmass(id(idx))
         do ip=1,4
            pjj(ip) = pjj(ip) + p(ip,idx)
         enddo
         ptjet(idx) = get_pt(p(1:4,idx))
      endif
   enddo
   maxptjet = maxval(ptjet)
   minptjet = minval(ptjet, mask=.not.all(p(1:4,5:6).eq.0d0, 1))
   polemass(7) = getmass(id(7)) ! Pole mass of the JJ system
 
   pjjhstar = pjj + phstar
   if(polemass(5).gt.polemass(6)) then
      pj(:,1)=p(:,5)
      pj(:,2)=p(:,6)
   else
      pj(:,1)=p(:,6)
      pj(:,2)=p(:,5)
      call swapr(polemass(5),polemass(6)) ! will use polemass(5) as the greater mass below
   endif
   pjhstar(1:4) = pj(1:4,1) + phstar(1:4)
 
   ! Determine the appropriate factorization scale for the chosen scheme from pole and invariant masses
   if(facscheme .eq. krenfacscheme_mhstar) then
      mu_fact = get_minv(phstar(1:4))
 
   elseif(facscheme .eq. -krenfacscheme_mhstar) then
      mu_fact = polemass(4)
 
   elseif(facscheme .eq. krenfacscheme_mjjhstar) then
      mu_fact = get_minv(pjjhstar(1:4))
   elseif(facscheme .eq. -krenfacscheme_mjjhstar) then
      mu_fact = polemass(3)
 
   elseif(facscheme .eq. krenfacscheme_mjj_mhstar) then
      mu_fact = get_minv(pjj(1:4))+get_minv(phstar(1:4))
   elseif(facscheme .eq. -krenfacscheme_mjj_mhstar) then
      mu_fact = polemass(4)+polemass(7)
   elseif(facscheme .eq. krenfacscheme_mj_mj_mhstar) then
      mu_fact = get_minv(pj(1:4,1))+get_minv(pj(1:4,2))+get_minv(phstar(1:4))
   elseif(facscheme .eq. -krenfacscheme_mj_mj_mhstar) then
      mu_fact = polemass(4)+polemass(5)+polemass(6)
 
   elseif(facscheme .eq. krenfacscheme_mjj) then
      mu_fact = get_minv(pjj(1:4))
   elseif(facscheme .eq. -krenfacscheme_mjj) then
      mu_fact = polemass(7)
   elseif(facscheme .eq. krenfacscheme_mj_mj) then
      mu_fact = get_minv(pjj(1:4))
   elseif(facscheme .eq. -krenfacscheme_mj_mj) then
      mu_fact = polemass(5)+polemass(6)
 
   elseif(facscheme .eq. krenfacscheme_mjhstar) then
      mu_fact = get_minv(pjhstar(1:4))
   elseif(facscheme .eq. krenfacscheme_mj_mhstar) then
      mu_fact = get_minv(pj(1:4,1))+get_minv(phstar(1:4))
   elseif((facscheme .eq. -krenfacscheme_mjhstar) .or. (facscheme .eq. -krenfacscheme_mj_mhstar)) then
      mu_fact = polemass(4)+polemass(5)
   elseif(facscheme .eq. krenfacscheme_mj) then
      mu_fact = get_minv(pj(1:4,1))
   elseif(facscheme .eq. -krenfacscheme_mj) then
      mu_fact = polemass(5)
   elseif(facscheme .eq. krenfacscheme_maxptj) then
      mu_fact = maxptjet
   elseif(facscheme .eq. krenfacscheme_minptj) then
      mu_fact = minptjet
   elseif(facscheme .eq. krenfacscheme_custom_scale) then
      call random_number(mu_fact)
      mu_fact = (customupperscalebound - customlowerscalebound)*mu_fact + customlowerscalebound
      mu_fact = mu_fact/100 !correction to make sure that you are operating in single GeVs, not hundreds of GeV
   else
      call error("This should never be able to happen.")
   endif
 
   ! Do the same for the renormalization scale
   if(renscheme .eq. krenfacscheme_mhstar) then
      mu_ren = get_minv(phstar(1:4))
 
   elseif(renscheme .eq. -krenfacscheme_mhstar) then
      mu_ren = polemass(4)
 
   elseif(renscheme .eq. krenfacscheme_mjjhstar) then
      mu_ren = get_minv(pjjhstar(1:4))
   elseif(renscheme .eq. -krenfacscheme_mjjhstar) then
      mu_ren = polemass(3)
 
   elseif(renscheme .eq. krenfacscheme_mjj_mhstar) then
      mu_ren = get_minv(pjj(1:4))+get_minv(phstar(1:4))
   elseif(renscheme .eq. -krenfacscheme_mjj_mhstar) then
      mu_ren = polemass(4)+polemass(7)
   elseif(renscheme .eq. krenfacscheme_mj_mj_mhstar) then
      mu_ren = get_minv(pj(1:4,1))+get_minv(pj(1:4,2))+get_minv(phstar(1:4))
   elseif(renscheme .eq. -krenfacscheme_mj_mj_mhstar) then
      mu_ren = polemass(4)+polemass(5)+polemass(6)
 
   elseif(renscheme .eq. krenfacscheme_mjj) then
      mu_ren = get_minv(pjj(1:4))
   elseif(renscheme .eq. -krenfacscheme_mjj) then
      mu_ren = polemass(7)
   elseif(renscheme .eq. krenfacscheme_mj_mj) then
      mu_ren = get_minv(pjj(1:4))
   elseif(renscheme .eq. -krenfacscheme_mj_mj) then
      mu_ren = polemass(5)+polemass(6)
 
   elseif(renscheme .eq. krenfacscheme_mjhstar) then
      mu_ren = get_minv(pjhstar(1:4))
   elseif(renscheme .eq. krenfacscheme_mj_mhstar) then
      mu_ren = get_minv(pj(1:4,1))+get_minv(phstar(1:4))
   elseif((renscheme .eq. -krenfacscheme_mjhstar) .or. (renscheme .eq. -krenfacscheme_mj_mhstar)) then
      mu_ren = polemass(4)+polemass(5)
   elseif(renscheme .eq. krenfacscheme_mj) then
      mu_ren = get_minv(pj(1:4,1))
   elseif(renscheme .eq. -krenfacscheme_mj) then
      mu_ren = polemass(5)
   elseif(renscheme .eq. krenfacscheme_maxptj) then
      mu_ren = maxptjet
   elseif(renscheme .eq. krenfacscheme_minptj) then
      mu_ren = minptjet
   else
      call error("This should never be able to happen.")
   endif
 
   ! Never ever allow the scales to go negative
   mu_fact = abs(mu_fact) * mufacmultiplier
   mu_ren = abs(mu_ren) * murenmultiplier
 
return
end subroutine setrunningscales
 
SUBROUTINE setpdfs(x1,x2,pdf)
use modparameters
implicit none
real(8) :: x1,x2,pdfscale
real(8) :: upv(1:2),dnv(1:2),usea(1:2),dsea(1:2),str(1:2),chm(1:2),bot(1:2),glu(1:2),phot(1:2),sbar(1:2),cbar(1:2),bbar(1:2)
integer,parameter :: swpdf_u=1, swpdf_d=1, swpdf_c=1, swpdf_s=1, swpdf_b=1, swpdf_g=1
real(8) :: pdf(-6:6,1:2),nnpdf(1:2,-6:7)
 
        pdfscale=mu_fact/gev
        pdf(:,:) = 0d0
 
#if useLHAPDF==1
        call evolvepdf(x1,pdfscale,nnpdf(1,-6:7))
        call evolvepdf(x2,pdfscale,nnpdf(2,-6:7))
            nnpdf(1,-6:7) = nnpdf(1,-6:7)/x1
            nnpdf(2,-6:7) = nnpdf(2,-6:7)/x2
 
            pdf(up_,1)   = nnpdf(1,+2)         * swpdf_u
            pdf(aup_,1)  = nnpdf(1,-2)         * swpdf_u
            pdf(dn_,1)   = nnpdf(1,+1)         * swpdf_d
            pdf(adn_,1)  = nnpdf(1,-1)         * swpdf_d
            pdf(chm_,1)  = nnpdf(1,+4)         * swpdf_c
            pdf(achm_,1) = nnpdf(1,-4)         * swpdf_c
            pdf(str_,1)  = nnpdf(1,+3)         * swpdf_s
            pdf(astr_,1) = nnpdf(1,-3)         * swpdf_s
            pdf(bot_,1)  = nnpdf(1,+5)         * swpdf_b
            pdf(abot_,1) = nnpdf(1,-5)         * swpdf_b
            pdf(0,1)     = nnpdf(1,+0)         * swpdf_g
 
            pdf(up_,2)   = nnpdf(2,+2)         * swpdf_u
            pdf(aup_,2)  = nnpdf(2,-2)         * swpdf_u
            pdf(dn_,2)   = nnpdf(2,+1)         * swpdf_d
            pdf(adn_,2)  = nnpdf(2,-1)         * swpdf_d
            pdf(chm_,2)  = nnpdf(2,+4)         * swpdf_c
            pdf(achm_,2) = nnpdf(2,-4)         * swpdf_c
            pdf(str_,2)  = nnpdf(2,+3)         * swpdf_s
            pdf(astr_,2) = nnpdf(2,-3)         * swpdf_s
            pdf(bot_,2)  = nnpdf(2,+5)         * swpdf_b
            pdf(abot_,2) = nnpdf(2,-5)         * swpdf_b
            pdf(0,2)     = nnpdf(2,+0)         * swpdf_g
 
            pdf(:,:) = dabs(pdf(:,:))
 
            RETURN
 
#else
        if( pdfset.eq.1 ) then
            call cteq6(x1,pdfscale,upv(1),dnv(1),usea(1),dsea(1),str(1),chm(1),bot(1),glu(1))
            call cteq6(x2,pdfscale,upv(2),dnv(2),usea(2),dsea(2),str(2),chm(2),bot(2),glu(2))
        elseif( pdfset.eq.2 ) then
            call getallpdfs("pdfs/mstw2008lo",0,x1,pdfscale,upv(1),dnv(1),usea(1),dsea(1),str(1),sbar(1),chm(1),cbar(1),bot(1),bbar(1),glu(1),phot(1))
            str(1)= (str(1)+sbar(1))/2d0
            chm(1)= (chm(1)+cbar(1))/2d0
            bot(1)= (bot(1)+bbar(1))/2d0
            upv(1)=upv(1)/x1
            dnv(1)=dnv(1)/x1
            usea(1)=usea(1)/x1
            dsea(1)=dsea(1)/x1
            str(1)=str(1)/x1
            chm(1)=chm(1)/x1
            bot(1)=bot(1)/x1
            glu(1)=glu(1)/x1
            phot(1)=phot(1)/x1
 
            call getallpdfs("pdfs/mstw2008lo",0,x2,pdfscale,upv(2),dnv(2),usea(2),dsea(2),str(2),sbar(2),chm(2),cbar(2),bot(2),bbar(2),glu(2),phot(2))
            str(2)= (str(2)+sbar(2))/2d0
            chm(2)= (chm(2)+cbar(2))/2d0
            bot(2)= (bot(2)+bbar(2))/2d0
            upv(2)=upv(2)/x2
            dnv(2)=dnv(2)/x2
            usea(2)=usea(2)/x2
            dsea(2)=dsea(2)/x2
            str(2)=str(2)/x2
            chm(2)=chm(2)/x2
            bot(2)=bot(2)/x2
            glu(2)=glu(2)/x2
            phot(2)=phot(2)/x2
 
        elseif( pdfset.ge.201 .and. pdfset.le.240) then
            call getallpdfs("pdfs/mstw2008lo.90cl",pdfset-200,x1,pdfscale,upv(1),dnv(1),usea(1),dsea(1),str(1),sbar(1),chm(1),cbar(1),bot(1),bbar(1),glu(1),phot(1))
            str(1)= (str(1)+sbar(1))/2d0
            chm(1)= (chm(1)+cbar(1))/2d0
            bot(1)= (bot(1)+bbar(1))/2d0
            upv(1)=upv(1)/x1
            dnv(1)=dnv(1)/x1
            usea(1)=usea(1)/x1
            dsea(1)=dsea(1)/x1
            str(1)=str(1)/x1
            chm(1)=chm(1)/x1
            bot(1)=bot(1)/x1
            glu(1)=glu(1)/x1
            phot(1)=phot(1)/x1
 
            call getallpdfs("pdfs/mstw2008lo.90cl",pdfset-200,x2,pdfscale,upv(2),dnv(2),usea(2),dsea(2),str(2),sbar(2),chm(2),cbar(2),bot(2),bbar(2),glu(2),phot(2))
            str(2)= (str(2)+sbar(2))/2d0
            chm(2)= (chm(2)+cbar(2))/2d0
            bot(2)= (bot(2)+bbar(2))/2d0
            upv(2)=upv(2)/x2
            dnv(2)=dnv(2)/x2
            usea(2)=usea(2)/x2
            dsea(2)=dsea(2)/x2
            str(2)=str(2)/x2
            chm(2)=chm(2)/x2
            bot(2)=bot(2)/x2
            glu(2)=glu(2)/x2
            phot(2)=phot(2)/x2
        elseif( pdfset.eq.3 ) then
 
            call nnevolvepdf(x1,pdfscale,nnpdf(1,-6:7))
            call nnevolvepdf(x2,pdfscale,nnpdf(2,-6:7))
            nnpdf(1,-6:7) = nnpdf(1,-6:7)/x1
            nnpdf(2,-6:7) = nnpdf(2,-6:7)/x2
 
    !       PROTON CONTENT
            pdf(up_,1)   = nnpdf(1,+2)         * swpdf_u
            pdf(aup_,1)  = nnpdf(1,-2)         * swpdf_u
            pdf(dn_,1)   = nnpdf(1,+1)         * swpdf_d
            pdf(adn_,1)  = nnpdf(1,-1)         * swpdf_d
            pdf(chm_,1)  = nnpdf(1,+4)         * swpdf_c
            pdf(achm_,1) = nnpdf(1,-4)         * swpdf_c
            pdf(str_,1)  = nnpdf(1,+3)         * swpdf_s
            pdf(astr_,1) = nnpdf(1,-3)         * swpdf_s
            pdf(bot_,1)  = nnpdf(1,+5)         * swpdf_b
            pdf(abot_,1) = nnpdf(1,-5)         * swpdf_b
            pdf(0,1)     = nnpdf(1,+0)         * swpdf_g
 
            pdf(up_,2)   = nnpdf(2,+2)         * swpdf_u
            pdf(aup_,2)  = nnpdf(2,-2)         * swpdf_u
            pdf(dn_,2)   = nnpdf(2,+1)         * swpdf_d
            pdf(adn_,2)  = nnpdf(2,-1)         * swpdf_d
            pdf(chm_,2)  = nnpdf(2,+4)         * swpdf_c
            pdf(achm_,2) = nnpdf(2,-4)         * swpdf_c
            pdf(str_,2)  = nnpdf(2,+3)         * swpdf_s
            pdf(astr_,2) = nnpdf(2,-3)         * swpdf_s
            pdf(bot_,2)  = nnpdf(2,+5)         * swpdf_b
            pdf(abot_,2) = nnpdf(2,-5)         * swpdf_b
            pdf(0,2)     = nnpdf(2,+0)         * swpdf_g
 
            pdf(:,:) = dabs(pdf(:,:))
            RETURN
        else
            print *, "PDFSet",pdfset,"not available!"
            stop
        endif
#endif
 
IF( collider.EQ.1 ) THEN
!       PROTON CONTENT
        pdf(up_,1)   = (upv(1) + usea(1))  * swpdf_u
        pdf(aup_,1)  = usea(1)             * swpdf_u
        pdf(dn_,1)   = (dnv(1) + dsea(1))  * swpdf_d
        pdf(adn_,1)  = dsea(1)             * swpdf_d
        pdf(chm_,1)  = chm(1)              * swpdf_c
        pdf(achm_,1) = chm(1)              * swpdf_c
        pdf(str_,1)  = str(1)              * swpdf_s
        pdf(astr_,1) = str(1)              * swpdf_s
        pdf(bot_,1)  = bot(1)              * swpdf_b
        pdf(abot_,1) = bot(1)              * swpdf_b
        pdf(0,1)     = glu(1)              * swpdf_g
 
!       PROTON CONTENT
        pdf(up_,2)   = (upv(2) + usea(2))  * swpdf_u
        pdf(aup_,2)  = usea(2)             * swpdf_u
        pdf(dn_,2)   = (dnv(2) + dsea(2))  * swpdf_d
        pdf(adn_,2)  = dsea(2)             * swpdf_d
        pdf(chm_,2)  = chm(2)              * swpdf_c
        pdf(achm_,2) = chm(2)              * swpdf_c
        pdf(str_,2)  = str(2)              * swpdf_s
        pdf(astr_,2) = str(2)              * swpdf_s
        pdf(bot_,2)  = bot(2)              * swpdf_b
        pdf(abot_,2) = bot(2)              * swpdf_b
        pdf(0,2)     = glu(2)              * swpdf_g
 
ELSEIF( collider.EQ.2 ) THEN
!       PROTON CONTENT
        pdf(up_,1)   = (upv(1) + usea(1))  * swpdf_u
        pdf(aup_,1)  = usea(1)             * swpdf_u
        pdf(dn_,1)   = (dnv(1) + dsea(1))  * swpdf_d
        pdf(adn_,1)  = dsea(1)             * swpdf_d
        pdf(chm_,1)  = chm(1)              * swpdf_c
        pdf(achm_,1) = chm(1)              * swpdf_c
        pdf(str_,1)  = str(1)              * swpdf_s
        pdf(astr_,1) = str(1)              * swpdf_s
        pdf(bot_,1)  = bot(1)              * swpdf_b
        pdf(abot_,1) = bot(1)              * swpdf_b
        pdf(0,1)     = glu(1)              * swpdf_g
 
!       ANTI-PROTON CONTENT
        pdf(up_,2)   = usea(2)             * swpdf_u
        pdf(aup_,2)  = (upv(2)+usea(2))    * swpdf_u
        pdf(dn_,2)   = dsea(2)             * swpdf_d
        pdf(adn_,2)  = (dnv(2) + dsea(2))  * swpdf_d
        pdf(chm_,2)  = chm(2)              * swpdf_c
        pdf(achm_,2) = chm(2)              * swpdf_c
        pdf(str_,2)  = str(2)              * swpdf_s
        pdf(astr_,2) = str(2)              * swpdf_s
        pdf(bot_,2)  = bot(2)              * swpdf_b
        pdf(abot_,2) = bot(2)              * swpdf_b
        pdf(0,2)     = glu(2)              * swpdf_g
 
ENDIF
 
pdf(:,:) = dabs(pdf(:,:))
 
 
RETURN
END SUBROUTINE
 
SUBROUTINE cteq6(X,SCALE,UPV,DNV,USEA,DSEA,STR,CHM,BOT,GLU)
implicit none
double precision X,SCALE,UPV,DNV,USEA,DSEA,STR,CHM,BOT,GLU
double precision Q,xsave,qsave,Ctq6Pdf,D,U
 
         q=scale
         xsave=x
         qsave=q
         u =         ctq6pdf(1,x,q)
         d =         ctq6pdf(2,x,q)
         usea =      ctq6pdf(-1,x,q)
         dsea =      ctq6pdf(-2,x,q)
         str =       ctq6pdf(3,x,q)
         chm =       ctq6pdf(4,x,q)
         bot =       ctq6pdf(5,x,q)
         glu  =      ctq6pdf(0,x,q)
         upv=u-usea
         dnv=d-dsea
         x=xsave
         q=qsave
RETURN
END SUBROUTINE
 
!! QCD scale from MCFM
!! Implementation into JHUGen by Ulascan Sarica, Dec. 2015
!subroutine EvalAlphaS()
!   use ModParameters
!   IMPLICIT NONE
!#if useLHAPDF==1
!!--- This is simply a wrapper to the LHAPDF implementation of the running coupling alphas, in the style of the native MCFM routine
!   DOUBLE PRECISION alphasPDF
!   REAL(DP) :: Q
!      Q = Mu_Ren/GeV
!      alphas=alphasPDF(Q)
!#else
!!     Evaluation of strong coupling constant alphas
!!     Original Author: R.K. Ellis
!!     q -- Scale at which alpha_s is to be evaluated
!!     alphas_mz -- ModParameters value of alpha_s at the mass of the Z-boson
!!     nloops_pdf -- ModParameters value of the number of loops (1,2, or 3) at which the beta function is evaluated to determine running.
!!     If you somehow need a more complete implementation, check everything at or before commit 28472c5bfee128dde458fd4929b4d3ece9519ab8
!   INTEGER, PARAMETER :: NF6=6
!   INTEGER, PARAMETER :: NF5=5
!   INTEGER, PARAMETER :: NF4=4
!   INTEGER, PARAMETER :: NF3=3
!   INTEGER, PARAMETER :: NF2=2
!   INTEGER, PARAMETER :: NF1=1
!
!      IF (Mu_Ren .LE. 0d0) THEN
!         WRITE(6,*) 'ModKinematics::EvalAlphaS: Mu_Ren .le. 0, Mu_Ren (GeV) = ',(Mu_Ren*GeV)
!         stop
!      ENDIF
!      IF (nQflavors_pdf .NE. NF5) THEN
!         WRITE(6,*) 'ModKinematics::EvalAlphaS: nQflavors_pdf invalid, nQflavors_pdf = ',nQflavors_pdf
!         WRITE(6,*) 'ModKinematics::EvalAlphaS: Check 28472c5bfee128dde458fd4929b4d3ece9519ab8'
!         stop
!      ENDIF
!      IF (nloops_pdf .NE. 1) THEN
!         WRITE(6,*) 'ModKinematics::EvalAlphaS: nloops_pdf invalid, nloops_pdf = ',nloops_pdf
!         WRITE(6,*) 'ModKinematics::EvalAlphaS: Check 28472c5bfee128dde458fd4929b4d3ece9519ab8'
!         stop
!      ENDIF
!
!      alphas=alphas_mz/(1.0_dp+alphas_mz*B0_PDF(NF5)*2.0_dp*dlog((Mu_Ren/zmass_pdf)))
!#endif
!      ! Calculate the derived couplings
!      call ComputeQCDVariables()
!   RETURN
!end subroutine EvalAlphaS
 
 
 
 
END MODULE