modkinematics Module Reference

Functions/Subroutines
subroutine	shiftmass (p1, p2, m1, m2, p1hat, p2hat, MassWeight)
integer function	zlepbranching (xRnd)
integer function	zlepbranching_flat (xRnd)
integer function	zlepplustaubranching (xRnd)
integer function	zlepplustaubranching_flat (xRnd)
integer function	znubranching (xRnd)
integer function	znubranching_flat (xRnd)
integer function	zquabranching_flat (xRnd)
integer function	zquabranching (xRnd)
integer function	zanybranching_flat (xRnd)
integer function	zanybranching (xRnd)
integer function	wlepbranching (xRnd)
integer function	wlepbranching_flat (xRnd)
integer function	wlepplustaubranching (xRnd)
integer function	wlepplustaubranching_flat (xRnd)
integer function	wquaupbranching (xRnd)
integer function	wquaupbranching_flat (xRnd)
integer function	wanybranching_flat (xRnd)
integer function	wanybranching (xRnd)
subroutine	vbranching (DecayMode, MY_IDUP, ICOLUP, CombWeight, ColorBase)
subroutine	vvbranchings (MY_IDUP, ICOLUP, CombWeight, ColorBase)
integer function	getckmpartner (Flavor)
integer function	getckmpartner_flat (Flavor)
real(8) function, public	getbwpropagator (sHat, scheme)
real(8) function, public	reweightbwpropagator (sHat)
subroutine	boost2lab (x1, x2, NumPart, Mom)
subroutine	lorentz (vector, boost)
double precision function	kronecker_delta (i, j)
subroutine, public	setrunningscales (p, id)
subroutine, public	setpdfs (x1, x2, pdf)
subroutine	cteq6 (X, SCALE, UPV, DNV, USEA, DSEA, STR, CHM, BOT, GLU)

Function/Subroutine Documentation

boost2lab()

subroutine modkinematics::boost2lab	(	real(8)	x1,
		real(8)	x2,
		integer	NumPart,
		real(8), dimension(1:4,1:numpart)	Mom
	)

Definition at line 1099 of file mod_Kinematics.F90.

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

cteq6()

subroutine modkinematics::cteq6	(	double precision	X,
		double precision	SCALE,
		double precision	UPV,
		double precision	DNV,
		double precision	USEA,
		double precision	DSEA,
		double precision	STR,
		double precision	CHM,
		double precision	BOT,
		double precision	GLU
	)

Definition at line 1552 of file mod_Kinematics.F90.

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

getbwpropagator()

real(8) function, public modkinematics::getbwpropagator	(	real(8)	sHat,
		integer	scheme
	)

Definition at line 1018 of file mod_Kinematics.F90.

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

getckmpartner()

integer function modkinematics::getckmpartner

(

integer

Flavor

)

Definition at line 824 of file mod_Kinematics.F90.

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

getckmpartner_flat()

integer function modkinematics::getckmpartner_flat

(

integer

Flavor

)

Definition at line 942 of file mod_Kinematics.F90.

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

kronecker_delta()

double precision function modkinematics::kronecker_delta	(	integer	i,
		integer	j
	)

Definition at line 1179 of file mod_Kinematics.F90.

      integer i,j
      if(i.eq.j)then
        kronecker_delta = 1d0
      else
        kronecker_delta = 0d0
      endif
 
      return

lorentz()

subroutine modkinematics::lorentz	(	double precision, dimension(4)	vector,
		double precision, dimension(4)	boost
	)

Definition at line 1126 of file mod_Kinematics.F90.

 
      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

reweightbwpropagator()

real(8) function, public modkinematics::reweightbwpropagator

(

real(8)

sHat

)

Definition at line 1078 of file mod_Kinematics.F90.

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

setpdfs()

subroutine, public modkinematics::setpdfs	(	real(8)	x1,
		real(8)	x2,
		real(8), dimension(-6:6,1:2)	pdf
	)

Definition at line 1344 of file mod_Kinematics.F90.

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

setrunningscales()

subroutine, public modkinematics::setrunningscales	(	real(dp), dimension(1:4,4:6), intent(in)	p,
		integer, dimension(4:7), intent(in)	id
	)

Definition at line 1194 of file mod_Kinematics.F90.

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

shiftmass()

subroutine modkinematics::shiftmass	(	real(8), dimension(1:4), intent(in)	p1,
		real(8), dimension(1:4), intent(in)	p2,
		real(8)	m1,
		real(8)	m2,
		real(8), dimension(1:4)	p1hat,
		real(8), dimension(1:4)	p2hat,
		real(8), optional	MassWeight
	)

Definition at line 11 of file mod_Kinematics.F90.

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

vbranching()

subroutine modkinematics::vbranching	(	integer, intent(in)	DecayMode,
		integer, dimension(1:3)	MY_IDUP,
		integer, dimension(1:2,1:2)	ICOLUP,
		real(8), intent(out)	CombWeight,
		integer, optional	ColorBase
	)

Definition at line 585 of file mod_Kinematics.F90.

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

vvbranchings()

subroutine modkinematics::vvbranchings	(	integer, dimension(4:9)	MY_IDUP,
		integer, dimension(1:2,6:9)	ICOLUP,
		real(8), intent(out)	CombWeight,
		integer, optional	ColorBase
	)

Definition at line 768 of file mod_Kinematics.F90.

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

wanybranching()

integer function modkinematics::wanybranching

(

real(8)

xRnd

)

Definition at line 536 of file mod_Kinematics.F90.

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

wanybranching_flat()

integer function modkinematics::wanybranching_flat

(

real(8)

xRnd

)

Definition at line 508 of file mod_Kinematics.F90.

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

wlepbranching()

integer function modkinematics::wlepbranching

(

real(8)

xRnd

)

Definition at line 394 of file mod_Kinematics.F90.

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

wlepbranching_flat()

integer function modkinematics::wlepbranching_flat

(

real(8)

xRnd

)

Definition at line 415 of file mod_Kinematics.F90.

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

wlepplustaubranching()

integer function modkinematics::wlepplustaubranching

(

real(8)

xRnd

)

Definition at line 430 of file mod_Kinematics.F90.

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

wlepplustaubranching_flat()

integer function modkinematics::wlepplustaubranching_flat

(

real(8)

xRnd

)

Definition at line 454 of file mod_Kinematics.F90.

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

wquaupbranching()

integer function modkinematics::wquaupbranching

(

real(8)

xRnd

)

Definition at line 471 of file mod_Kinematics.F90.

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

wquaupbranching_flat()

integer function modkinematics::wquaupbranching_flat

(

real(8)

xRnd

)

Definition at line 493 of file mod_Kinematics.F90.

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

zanybranching()

integer function modkinematics::zanybranching

(

real(8)

xRnd

)

Definition at line 323 of file mod_Kinematics.F90.

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

zanybranching_flat()

integer function modkinematics::zanybranching_flat

(

real(8)

xRnd

)

Definition at line 247 of file mod_Kinematics.F90.

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

zlepbranching()

integer function modkinematics::zlepbranching

(

real(8)

xRnd

)

Definition at line 69 of file mod_Kinematics.F90.

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

zlepbranching_flat()

integer function modkinematics::zlepbranching_flat

(

real(8)

xRnd

)

Definition at line 91 of file mod_Kinematics.F90.

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

zlepplustaubranching()

integer function modkinematics::zlepplustaubranching

(

real(8)

xRnd

)

Definition at line 107 of file mod_Kinematics.F90.

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

zlepplustaubranching_flat()

integer function modkinematics::zlepplustaubranching_flat

(

real(8)

xRnd

)

Definition at line 131 of file mod_Kinematics.F90.

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

znubranching()

integer function modkinematics::znubranching

(

real(8)

xRnd

)

Definition at line 148 of file mod_Kinematics.F90.

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

znubranching_flat()

integer function modkinematics::znubranching_flat

(

real(8)

xRnd

)

Definition at line 172 of file mod_Kinematics.F90.

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

zquabranching()

integer function modkinematics::zquabranching

(

real(8)

xRnd

)

Definition at line 219 of file mod_Kinematics.F90.

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

zquabranching_flat()

integer function modkinematics::zquabranching_flat

(

real(8)

xRnd

)

Definition at line 189 of file mod_Kinematics.F90.

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