pysubh.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysubh.f

 
C*********************************************************************
 
C...PYSUBH
C...This routine computes the renormalization group improved
C...values of Higgs masses and couplings in the MSSM.
 
C...Program based on the work by M. Carena, J.R. Espinosa,
c...M. Quiros and C.E.M. Wagner, CERN-preprint CERN-TH/95-45
 
C...Input: MA,TANB = TAN(BETA),MQ,MUR,MTOP,AU,AD,MU
C...All masses in GeV units. MA is the CP-odd Higgs mass,
C...MTOP is the physical top mass, MQ and MUR are the soft
C...supersymmetry breaking mass parameters of left handed
C...and right handed stops respectively, AU and AD are the
C...stop and sbottom trilinear soft breaking terms,
C...respectively,  and MU is the supersymmetric
C...Higgs mass parameter. We use the  conventions from
C...the physics report of Haber and Kane: left right
C...stop mixing term proportional to (AU - MU/TANB)
C...We use as input TANB defined at the scale MTOP
 
C...Output: MH,HM,MHCH, SA = SIN(ALPHA), CA= COS(ALPHA), TANBA
C...where MH and HM are the lightest and heaviest CP-even
C...Higgs masses, MHCH is the charged Higgs mass and
C...ALPHA is the Higgs mixing angle
C...TANBA is the angle TANB at the CP-odd Higgs mass scale
 
C...Range of validity:
C...(STOP1**2 - STOP2**2)/(STOP2**2 + STOP1**2) < 0.5
C...(SBOT1**2 - SBOT2**2)/(SBOT2**2 + SBOT2**2) < 0.5
C...where STOP1, STOP2, SBOT1 and SBOT2 are the stop and
C...are the sbottom  mass eigenvalues, respectively. This
C...range automatically excludes the existence of tachyons.
C...For the charged Higgs mass computation, the method is
C...valid if
C...2 * |MB * AD* TANB|  < M_SUSY**2,  2 * |MTOP * AU| < M_SUSY**2
C...2 * |MB * MU * TANB| < M_SUSY**2,  2 * |MTOP * MU| < M_SUSY**2
C...where M_SUSY**2 is the average of the squared stop mass
C...eigenvalues, M_SUSY**2 = (STOP1**2 + STOP2**2)/2. The sbottom
C...masses have been assumed to be of order of the stop ones
C...M_SUSY**2 = (MQ**2 + MUR**2)*0.5 + MTOP**2
 
      SUBROUTINE pysubh (XMA,TANB,XMQ,XMUR,XMTOP,AU,AD,XMU,XMH,XHM,
     &xmhch,sa,ca,tanba)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement to help give large particle numbers.
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pyhtri/hhh(7)
      SAVE /pydat1/,/pydat2/
 
C...Local variables.
      DOUBLE PRECISION pyalem,pyalps
      DOUBLE PRECISION tanb,xmq,xmur,xmtop,au,ad,xmu,xmh,xhm
      DOUBLE PRECISION xmhch,sa,ca
      DOUBLE PRECISION xma,aem,alp1,alp2,alph3z,v,pi
      DOUBLE PRECISION q02
      DOUBLE PRECISION tanba,tanbt,xmb,alp3
      DOUBLE PRECISION rmtop,xms,t,sinb,cosb
      DOUBLE PRECISION xlam1,xlam2,xlam3,xlam4,xlam5,xlam6
      DOUBLE PRECISION xlam7,xau,xad,g1,g2,g3,hu,hd,hu2
      DOUBLE PRECISION hd2,hu4,hd4,sinbt,cosbt
      DOUBLE PRECISION trm2,detm2,xmh2,xhm2,xmhch2
      DOUBLE PRECISION sinalp,cosalp,aud,pi2,xms2,xms4,ad2
      DOUBLE PRECISION au2,xmu2,xmz,xms3
 
      xmz = pmas(23,1)
      q02=xmz**2
      aem=pyalem(q02)
      alp1=aem/(1d0-paru(102))
      alp2=aem/paru(102)
      alph3z=pyalps(q02)
 
      alp1 = 0.0101d0
      alp2 = 0.0337d0
      alph3z = 0.12d0
 
      v = 174.1d0
      pi = paru(1)
      tanba = tanb
      tanbt = tanb
 
C...MBOTTOM(MTOP) = 3. GEV
      xmb = pymrun(5,xmtop**2)
      alp3 = alph3z/(1d0 +(11d0 - 10d0/3d0)/4d0/pi*alph3z*
     &log(xmtop**2/xmz**2))
 
C...RMTOP= RUNNING TOP QUARK MASS
      rmtop = xmtop/(1d0+4d0*alp3/3d0/pi)
      xms = ((xmq**2 + xmur**2)/2d0 + xmtop**2)**0.5d0
      t = log(xms**2/xmtop**2)
      sinb = tanb/((1d0 + tanb**2)**0.5d0)
      cosb = sinb/tanb
C...IF(MA.LE.XMTOP) TANBA = TANBT
      IF(xma.GT.xmtop)
     &tanba = tanbt*(1d0-3d0/32d0/pi**2*
     &(rmtop**2/v**2/sinb**2-xmb**2/v**2/cosb**2)*
     &log(xma**2/xmtop**2))
 
      sinbt = tanbt/sqrt(1d0 + tanbt**2)
      cosbt = 1d0/sqrt(1d0 + tanbt**2)
C      COS2BT = (TANBT**2 - 1D0)/(TANBT**2 + 1D0)
      g1 = sqrt(alp1*4d0*pi)
      g2 = sqrt(alp2*4d0*pi)
      g3 = sqrt(alp3*4d0*pi)
      hu = rmtop/v/sinbt
      hd =  xmb/v/cosbt
      hu2=hu*hu
      hd2=hd*hd
      hu4=hu2*hu2
      hd4=hd2*hd2
      au2=au**2
      ad2=ad**2
      xms2=xms**2
      xms3=xms**3
      xms4=xms2*xms2
      xmu2=xmu*xmu
      pi2=pi*pi
 
      xau = (2d0*au2/xms2)*(1d0 - au2/12d0/xms2)
      xad = (2d0*ad2/xms2)*(1d0 - ad2/12d0/xms2)
      aud = (-6d0*xmu2/xms2 - ( xmu2- ad*au)**2/xms4
     &+ 3d0*(au + ad)**2/xms2)/6d0
      xlam1 = ((g1**2 + g2**2)/4d0)*(1d0-3d0*hd2*t/8d0/pi2)
     &+(3d0*hd4/8d0/pi2) * (t + xad/2d0 + (3d0*hd2/2d0 + hu2/2d0
     &- 8d0*g3**2) * (xad*t + t**2)/16d0/pi2)
     &-(3d0*hu4* xmu**4/96d0/pi2/xms4) * (1+ (9d0*hu2 -5d0* hd2
     &-  16d0*g3**2) *t/16d0/pi2)
      xlam2 = ((g1**2 + g2**2)/4d0)*(1d0-3d0*hu2*t/8d0/pi2)
     &+(3d0*hu4/8d0/pi2) * (t + xau/2d0 + (3d0*hu2/2d0 + hd2/2d0
     &- 8d0*g3**2) * (xau*t + t**2)/16d0/pi2)
     &-(3d0*hd4* xmu**4/96d0/pi2/xms4) * (1+ (9d0*hd2 -5d0* hu2
     &-  16d0*g3**2) *t/16d0/pi2)
      xlam3 = ((g2**2 - g1**2)/4d0)*(1d0-3d0*
     &(hu2 + hd2)*t/16d0/pi2)
     &+(6d0*hu2*hd2/16d0/pi2) * (t + aud/2d0 + (hu2 + hd2
     &- 8d0*g3**2) * (aud*t + t**2)/16d0/pi2)
     &+(3d0*hu4/96d0/pi2) * (3d0*xmu2/xms2 - xmu2*au2/
     &xms4)* (1d0+ (6d0*hu2 -2d0* hd2/2d0
     &-  16d0*g3**2) *t/16d0/pi2)
     &+(3d0*hd4/96d0/pi2) * (3d0*xmu2/xms2 - xmu2*ad2/
     &xms4)*(1d0+ (6d0*hd2 -2d0* hu2
     &-  16d0*g3**2) *t/16d0/pi2)
      xlam4 = (- g2**2/2d0)*(1d0-3d0*(hu2 + hd2)*t/16d0/pi2)
     &-(6d0*hu2*hd2/16d0/pi2) * (t + aud/2d0 + (hu2 + hd2
     &- 8d0*g3**2) * (aud*t + t**2)/16d0/pi2)
     &+(3d0*hu4/96d0/pi2) * (3d0*xmu2/xms2 - xmu2*au2/
     &xms4)*
     &(1+ (6d0*hu2 -2d0* hd2
     &-  16d0*g3**2) *t/16d0/pi2)
     &+(3d0*hd4/96d0/pi2) * (3d0*xmu2/xms2 - xmu2*ad2/
     &xms4)*
     &(1+ (6d0*hd2 -2d0* hu2/2d0
     &-  16d0*g3**2) *t/16d0/pi2)
      xlam5 = -(3d0*hu4* xmu2*au2/96d0/pi2/xms4) *
     &(1- (2d0*hd2 -6d0* hu2 + 16d0*g3**2) *t/16d0/pi2)
     &-(3d0*hd4* xmu2*ad2/96d0/pi2/xms4) *
     &(1- (2d0*hu2 -6d0* hd2 + 16d0*g3**2) *t/16d0/pi2)
      xlam6 = (3d0*hu4* xmu**3*au/96d0/pi2/xms4) *
     &(1- (7d0*hd2/2d0 -15d0* hu2/2d0 + 16d0*g3**2) *t/16d0/pi2)
     &+(3d0*hd4* xmu *(ad**3/xms3 - 6d0*ad/xms )/96d0/pi2/xms) *
     &(1- (hu2/2d0 -9d0* hd2/2d0 + 16d0*g3**2) *t/16d0/pi2)
      xlam7 = (3d0*hd4* xmu**3*ad/96d0/pi2/xms4) *
     &(1- (7d0*hu2/2d0 -15d0* hd2/2d0 + 16d0*g3**2) *t/16d0/pi2)
     &+(3d0*hu4* xmu *(au**3/xms3 - 6d0*au/xms )/96d0/pi2/xms) *
     &(1- (hd2/2d0 -9d0* hu2/2d0 + 16d0*g3**2) *t/16d0/pi2)
      hhh(1)=xlam1
      hhh(2)=xlam2
      hhh(3)=xlam3
      hhh(4)=xlam4
      hhh(5)=xlam5
      hhh(6)=xlam6
      hhh(7)=xlam7
      trm2 = xma**2 + 2d0*v**2* (xlam1* cosbt**2 +
     &2d0* xlam6*sinbt*cosbt
     &+ xlam5*sinbt**2 + xlam2* sinbt**2 + 2d0* xlam7*sinbt*cosbt
     &+ xlam5*cosbt**2)
      detm2 = 4d0*v**4*(-(sinbt*cosbt*(xlam3 + xlam4) +
     &xlam6*cosbt**2
     &+ xlam7* sinbt**2)**2 + (xlam1* cosbt**2 +
     &2d0* xlam6* cosbt*sinbt
     &+ xlam5*sinbt**2)*(xlam2* sinbt**2 +2d0* xlam7* cosbt*sinbt
     &+ xlam5*cosbt**2)) + xma**2*2d0*v**2 *
     &((xlam1* cosbt**2 +2d0*
     &xlam6* cosbt*sinbt + xlam5*sinbt**2)*cosbt**2 +
     &(xlam2* sinbt**2 +2d0* xlam7* cosbt*sinbt + xlam5*cosbt**2)
     &*sinbt**2
     &+2d0*sinbt*cosbt* (sinbt*cosbt*(xlam3
     &+ xlam4) + xlam6*cosbt**2
     &+ xlam7* sinbt**2))
 
      xmh2 = (trm2 - sqrt(trm2**2 - 4d0* detm2))/2d0
      xhm2 = (trm2 + sqrt(trm2**2 - 4d0* detm2))/2d0
      xhm = sqrt(xhm2)
      xmh = sqrt(xmh2)
      xmhch2 = xma**2 + (xlam5 - xlam4)* v**2
      xmhch = sqrt(xmhch2)
 
      sinalp = sqrt(((trm2**2 - 4d0* detm2)**0.5d0) -
     &((2d0*v**2*(xlam1* cosbt**2 + 2d0*
     &xlam6* cosbt*sinbt
     &+ xlam5*sinbt**2) + xma**2*sinbt**2)
     &- (2d0*v**2*(xlam2* sinbt**2 +2d0* xlam7* cosbt*sinbt
     &+ xlam5*cosbt**2) + xma**2*cosbt**2)))/
     &sqrt(((trm2**2 - 4d0* detm2)**0.5d0))/2d0**0.5d0
 
      cosalp = (2d0*(2d0*v**2*(sinbt*cosbt*(xlam3 + xlam4) +
     &xlam6*cosbt**2 + xlam7* sinbt**2) -
     &xma**2*sinbt*cosbt))/2d0**0.5d0/
     &sqrt(((trm2**2 - 4d0* detm2)**0.5d0)*
     &(((trm2**2 - 4d0* detm2)**0.5d0) -
     &((2d0*v**2*(xlam1* cosbt**2 + 2d0*
     &xlam6* cosbt*sinbt
     &+ xlam5*sinbt**2) + xma**2*sinbt**2)
     &- (2d0*v**2*(xlam2* sinbt**2 +2d0* xlam7* cosbt*sinbt
     &+ xlam5*cosbt**2) + xma**2*cosbt**2))))
 
      sa = -sinalp
      ca = -cosalp
 
  100 CONTINUE
 
      RETURN
      END