pythrg.f

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

 
C*********************************************************************
 
C...PYTHRG
C...Calculates the mass eigenstates of the third generation sfermions.
C...Created:  5-31-96
 
      SUBROUTINE pythrg
 
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/pymssm/imss(0:99),rmss(0:99)
      common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
     &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
      SAVE /pydat1/,/pydat2/,/pymssm/,/pyssmt/
 
C...Local variables.
      DOUBLE PRECISION beta
      DOUBLE PRECISION am2(2,2),rt(2,2),di(2,2)
      DOUBLE PRECISION xmz2,xmw2,tanb,xmu,cos2b,xmql2,xmqr2
      DOUBLE PRECISION xmf,xmf2,diff,same,xmf12,xmf22,small
      DOUBLE PRECISION atr,amqr,amql
      INTEGER id1(3),id2(3),id3(3),id4(3)
      INTEGER if,i,j,ii,jj,it,l
      LOGICAL dterm
      DATA small/1d-3/
      DATA id1/10,10,13/
      DATA id2/5,6,15/
      DATA id3/15,16,17/
      DATA id4/11,12,14/
      DATA dterm/.true./
 
      xmz2=pmas(23,1)**2
      xmw2=pmas(24,1)**2
      tanb=rmss(5)
      xmu=-rmss(4)
      beta=atan(tanb)
      cos2b=cos(2d0*beta)
 
C...OPTION TO FIX T1, T2, B1 MASSES AND MIXINGS
 
      iopt=imss(5)
      IF(iopt.EQ.1) THEN
        ctt=dcos(rmss(27))
        ctt2=ctt**2
        stt=dsin(rmss(27))
        stt2=stt**2
        xm12=rmss(10)**2
        xm22=rmss(12)**2
        xmql2=ctt2*xm12+stt2*xm22
        xmqr2=stt2*xm12+ctt2*xm22
        xmf2=pymrun(6,pmas(6,1)**2)**2
        atop=-xmu/tanb+ctt*stt*(xm12-xm22)/sqrt(xmf2)
        rmss(16)=atop
C......SUBTRACT OUT D-TERM AND FERMION MASS
        xmql2=xmql2-xmf2-(4d0*xmw2-xmz2)*cos2b/6d0
        xmqr2=xmqr2-xmf2+(xmw2-xmz2)*cos2b*2d0/3d0
        IF(xmql2.GE.0d0) THEN
          rmss(10)=sqrt(xmql2)
        ELSE
          rmss(10)=-sqrt(-xmql2)
        ENDIF
        IF(xmqr2.GE.0d0) THEN
          rmss(12)=sqrt(xmqr2)
        ELSE
          rmss(12)=-sqrt(-xmqr2)
        ENDIF
 
C SAME FOR BOTTOM SQUARK
        ctt=dcos(rmss(26))
        ctt2=ctt**2
        stt=dsin(rmss(26))
        stt2=stt**2
        xm22=rmss(11)**2
        xmf2=pymrun(5,pmas(6,1)**2)**2
        xmql2=sign(rmss(10)**2,rmss(10))-(2d0*xmw2+xmz2)*cos2b/6d0+xmf2
        IF(abs(ctt).GE..9999d0) THEN
          abot=-xmu*tanb
          xmqr2=rmss(11)**2
        ELSEIF(abs(ctt).LE.1d-4) THEN
          abot=-xmu*tanb
          xmqr2=rmss(11)**2
        ELSE
          xm12=(xmql2-stt2*xm22)/ctt2
          xmqr2=stt2*xm12+ctt2*xm22
          abot=-xmu*tanb+ctt*stt*(xm12-xm22)/sqrt(xmf2)
        ENDIF
        rmss(15)=abot
C......SUBTRACT OUT D-TERM AND FERMION MASS
        xmqr2=xmqr2-(xmw2-xmz2)*cos2b/3d0-xmf2
        IF(xmqr2.GE.0d0) THEN
          rmss(11)=sqrt(xmqr2)
        ELSE
          rmss(11)=-sqrt(-xmqr2)
        ENDIF
C SAME FOR TAU SLEPTON
        ctt=dcos(rmss(28))
        ctt2=ctt**2
        stt=dsin(rmss(28))
        stt2=stt**2
        xm12=rmss(13)**2
        xm22=rmss(14)**2
        xmql2=ctt2*xm12+stt2*xm22
        xmqr2=stt2*xm12+ctt2*xm22
        xmfr=pmas(15,1)
        xmf2=xmfr**2
        atau=-xmu*tanb+ctt*stt*(xm12-xm22)/sqrt(xmf2)
        rmss(17)=atau
C......SUBTRACT OUT D-TERM AND FERMION MASS
        xmql2=xmql2-xmf2+(-.5d0*xmz2+xmw2)*cos2b
        xmqr2=xmqr2-xmf2+(xmz2-xmw2)*cos2b
        IF(xmql2.GE.0d0) THEN
          rmss(13)=sqrt(xmql2)
        ELSE
          rmss(13)=-sqrt(-xmql2)
        ENDIF
        IF(xmqr2.GE.0d0) THEN
          rmss(14)=sqrt(xmqr2)
        ELSE
          rmss(14)=-sqrt(-xmqr2)
        ENDIF
      ENDIF
      DO 170 l=1,3
        amql=rmss(id1(l))
        IF(amql.LT.0d0) THEN
          xmql2=-amql**2
        ELSE
          xmql2=amql**2
        ENDIF
        atr=rmss(id3(l))
        amqr=rmss(id4(l))
        IF(amqr.LT.0d0) THEN
          xmqr2=-amqr**2
        ELSE
          xmqr2=amqr**2
        ENDIF
        if=id2(l)
        xmf=pymrun(IF,pmas(6,1)**2)
        xmf2=xmf**2
        am2(1,1)=xmql2+xmf2
        am2(2,2)=xmqr2+xmf2
        IF(am2(1,1).EQ.am2(2,2)) am2(2,2)=am2(2,2)*1.00001d0
        IF(dterm) THEN
          IF(l.EQ.1) THEN
            am2(1,1)=am2(1,1)-(2d0*xmw2+xmz2)*cos2b/6d0
            am2(2,2)=am2(2,2)+(xmw2-xmz2)*cos2b/3d0
            am2(1,2)=xmf*(atr+xmu*tanb)
          ELSEIF(l.EQ.2) THEN
            am2(1,1)=am2(1,1)+(4d0*xmw2-xmz2)*cos2b/6d0
            am2(2,2)=am2(2,2)-(xmw2-xmz2)*cos2b*2d0/3d0
            am2(1,2)=xmf*(atr+xmu/tanb)
          ELSEIF(l.EQ.3) THEN
            IF(imss(8).EQ.1) THEN
              am2(1,1)=rmss(6)**2
              am2(2,2)=rmss(7)**2
              am2(1,2)=0d0
              rmss(13)=rmss(6)
              rmss(14)=rmss(7)
            ELSE
              am2(1,1)=am2(1,1)-(-.5d0*xmz2+xmw2)*cos2b
              am2(2,2)=am2(2,2)-(xmz2-xmw2)*cos2b
              am2(1,2)=xmf*(atr+xmu*tanb)
            ENDIF
          ENDIF
        ENDIF
        am2(2,1)=am2(1,2)
        detm=am2(1,1)*am2(2,2)-am2(2,1)**2
        IF(detm.LT.0d0) THEN
          WRITE(mstu(11),*) id2(l),detm,am2
          CALL pyerrm(30,' NEGATIVE**2 MASS FOR SFERMION IN PYTHRG ')
        ENDIF
        same=0.5d0*(am2(1,1)+am2(2,2))
        diff=0.5d0*sqrt((am2(1,1)-am2(2,2))**2+4d0*am2(1,2)*am2(2,1))
        xmf12=same-diff
        xmf22=same+diff
        it=0
        IF(xmf22-xmf12.GT.0d0) THEN
          rt(1,1) = sqrt(max(0d0,(xmf22-am2(1,1))/(xmf22-xmf12)))
          rt(2,2) = rt(1,1)
          rt(1,2) = -sign(sqrt(max(0d0,1d0-rt(1,1)**2)),
     &    am2(1,2)/(xmf22-xmf12))
          rt(2,1) = -rt(1,2)
        ELSE
          rt(1,1) = 1d0
          rt(2,2) = rt(1,1)
          rt(1,2) = 0d0
          rt(2,1) = -rt(1,2)
        ENDIF
  100   CONTINUE
        it=it+1
 
        DO 140 i=1,2
          DO 130 jj=1,2
            di(i,jj)=0d0
            DO 120 ii=1,2
              DO 110 j=1,2
                di(i,jj)=di(i,jj)+rt(i,j)*am2(j,ii)*rt(jj,ii)
  110         CONTINUE
  120       CONTINUE
  130     CONTINUE
  140   CONTINUE
 
        IF(di(1,1).GT.di(2,2)) THEN
          WRITE(mstu(11),*) ' ERROR IN DIAGONALIZATION '
          WRITE(mstu(11),*) l,sqrt(xmf12),sqrt(xmf22)
          WRITE(mstu(11),*) am2
          WRITE(mstu(11),*) di
          WRITE(mstu(11),*) rt
          di(1,1)=-rt(2,1)
          di(2,2)=rt(1,2)
          di(1,2)=-rt(2,2)
          di(2,1)=rt(1,1)
          DO 160 i=1,2
            DO 150 j=1,2
              rt(i,j)=di(i,j)
  150       CONTINUE
  160     CONTINUE
          goto 100
        ELSEIF(abs(di(1,2)*di(2,1)/di(1,1)/di(2,2)).GT.small) THEN
          WRITE(mstu(11),*) ' ERROR IN DIAGONALIZATION,'//
     &    ' OFF DIAGONAL ELEMENTS '
          WRITE(mstu(11),*) 'MASSES = ',l,sqrt(xmf12),sqrt(xmf22)
          WRITE(mstu(11),*) di
          WRITE(mstu(11),*) ' ROTATION = ',rt
C...STOP
        ELSEIF(di(1,1).LT.0d0.OR.di(2,2).LT.0d0) THEN
          WRITE(mstu(11),*) ' ERROR IN DIAGONALIZATION,'//
     &    ' NEGATIVE MASSES '
          CALL pystop(111)
        ENDIF
        pmas(pycomp(ksusy1+if),1)=sqrt(xmf12)
        pmas(pycomp(ksusy2+if),1)=sqrt(xmf22)
        sfmix(IF,1)=rt(1,1)
        sfmix(IF,2)=rt(1,2)
        sfmix(IF,3)=rt(2,1)
        sfmix(IF,4)=rt(2,2)
  170 CONTINUE
 
C.....TAU SNEUTRINO MASS...L=3
 
      xarg=am2(1,1)+xmw2*cos2b
      IF(xarg.LT.0d0) THEN
        WRITE(mstu(11),*) ' PYTHRG:: TAU SNEUTRINO MASS IS NEGATIVE'//
     &  ' FROM THE SUM RULE. '
        WRITE(mstu(11),*) '  TRY A SMALLER VALUE OF TAN(BETA). '
        RETURN
      ELSE
        pmas(pycomp(ksusy1+16),1)=sqrt(xarg)
      ENDIF
 
      RETURN
      END