pytbbn.f

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

 
 
C*********************************************************************
 
C...PYTBBN
C...Calculates the three-body decay of gluinos into
C...neutralinos and third generation fermions.
 
      SUBROUTINE pytbbn(I,NN,E,XMGLU,GAM)
 
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.
      EXTERNAL pysimp,pylamf
      DOUBLE PRECISION pysimp,pylamf
      INTEGER lin,nn
      DOUBLE PRECISION cosd,sind,cosd2,sind2,cos2d,sin2d
      DOUBLE PRECISION hl,hr,fl,fr,hl2,hr2,fl2,fr2
      DOUBLE PRECISION xms2(2),xm,xm2,xmg,xmg2,xmr,xmr2
      DOUBLE PRECISION sbar,smin,smax,xmqa,w,grs,g(0:6),summe(0:100)
      DOUBLE PRECISION ff,hh,hfl,hfr,hrfl,hlfr,xmq4,xm24
      DOUBLE PRECISION xln1,xln2,b1,b2
      DOUBLE PRECISION e,xmglu,gam
      DOUBLE PRECISION hrb(4),hlb(4),flb(4),frb(4)
      SAVE hrb,hlb,flb,frb
      DOUBLE PRECISION alphaw,alphas
      DOUBLE PRECISION hlt(4),hrt(4),flt(4),frt(4)
      SAVE hlt,hrt,flt,frt
      DOUBLE PRECISION amn(4),an(4,4),zn(3)
      SAVE amn,an,zn
      DOUBLE PRECISION ambot,sinc,cosc
      DOUBLE PRECISION amtop,sina,cosa
      DOUBLE PRECISION sinw,cosw,tanw
      DOUBLE PRECISION rot1(4,4)
      LOGICAL ifirst
      SAVE ifirst
      DATA ifirst/.true./
 
      tanb=rmss(5)
      sinb=tanb/sqrt(1d0+tanb**2)
      cosb=sinb/tanb
      xw=paru(102)
      sinw=sqrt(xw)
      cosw=sqrt(1d0-xw)
      tanw=sinw/cosw
      amw=pmas(24,1)
      cosc=sfmix(5,1)
      sinc=sfmix(5,3)
      cosa=sfmix(6,1)
      sina=sfmix(6,3)
      ambot=pymrun(5,xmglu**2)
      amtop=pymrun(6,xmglu**2)
      w2=sqrt(2d0)
      fakt1=ambot/w2/amw/cosb
      fakt2=amtop/w2/amw/sinb
      IF(ifirst) THEN
        DO 110 ii=1,4
          amn(ii)=smz(ii)
          DO 100 j=1,4
            rot1(ii,j)=0d0
            an(ii,j)=0d0
  100     CONTINUE
  110   CONTINUE
        rot1(1,1)=cosw
        rot1(1,2)=-sinw
        rot1(2,1)=-rot1(1,2)
        rot1(2,2)=rot1(1,1)
        rot1(3,3)=cosb
        rot1(3,4)=sinb
        rot1(4,3)=-rot1(3,4)
        rot1(4,4)=rot1(3,3)
        DO 140 ii=1,4
          DO 130 j=1,4
            DO 120 jj=1,4
              an(ii,j)=an(ii,j)+zmix(ii,jj)*rot1(jj,j)
  120       CONTINUE
  130     CONTINUE
  140   CONTINUE
        DO 150 j=1,4
          zn(1)=-fakt2*(-sinb*an(j,3)+cosb*an(j,4))
          zn(2)=-2d0*w2/3d0*sinw*(tanw*an(j,2)-an(j,1))
          zn(3)=-2*w2/3d0*sinw*an(j,1)-w2*(0.5d0-2d0/3d0*
     &    xw)*an(j,2)/cosw
          hrt(j)=zn(1)*cosa-zn(3)*sina
          hlt(j)=zn(1)*cosa+zn(2)*sina
          flt(j)=zn(3)*cosa+zn(1)*sina
          frt(j)=zn(2)*cosa-zn(1)*sina
C          FLU(J)=ZN(3)
C          FRU(J)=ZN(2)
          zn(1)=-fakt1*(cosb*an(j,3)+sinb*an(j,4))
          zn(2)=w2/3d0*sinw*(tanw*an(j,2)-an(j,1))
          zn(3)=w2/3d0*sinw*an(j,1)+w2*(0.5d0-xw/3d0)*an(j,2)/cosw
          hrb(j)=zn(1)*cosc-zn(3)*sinc
          hlb(j)=zn(1)*cosc+zn(2)*sinc
          flb(j)=zn(3)*cosc+zn(1)*sinc
          frb(j)=zn(2)*cosc-zn(1)*sinc
C          FLD(J)=ZN(3)
C          FRD(J)=ZN(2)
  150   CONTINUE
C        AMST(1)=PMAS(PYCOMP(KSUSY1+6),1)
C        AMST(2)=PMAS(PYCOMP(KSUSY2+6),1)
C        AMSB(1)=PMAS(PYCOMP(KSUSY1+5),1)
C        AMSB(2)=PMAS(PYCOMP(KSUSY2+5),1)
        ifirst=.false.
      ENDIF
 
      IF(nint(3d0*e).EQ.2) THEN
        hl=hlt(i)
        hr=hrt(i)
        fl=flt(i)
        fr=frt(i)
        cosd=sfmix(6,1)
        sind=sfmix(6,3)
        xms2(1)=pmas(pycomp(ksusy1+6),1)**2
        xms2(2)=pmas(pycomp(ksusy2+6),1)**2
        xm=pmas(6,1)
      ELSE
        hl=hlb(i)
        hr=hrb(i)
        fl=flb(i)
        fr=frb(i)
        cosd=sfmix(5,1)
        sind=sfmix(5,3)
        xms2(1)=pmas(pycomp(ksusy1+5),1)**2
        xms2(2)=pmas(pycomp(ksusy2+5),1)**2
        xm=pmas(5,1)
      ENDIF
      cosd2=cosd*cosd
      sind2=sind*sind
      cos2d=cosd2-sind2
      sin2d=sind*cosd*2d0
      hl2=hl*hl
      hr2=hr*hr
      fl2=fl*fl
      fr2=fr*fr
      ff=fl*fr
      hh=hl*hr
      hfl=hl*fl
      hfr=hr*fr
      hrfl=hr*fl
      hlfr=hl*fr
      xm2=xm*xm
      xmg=xmglu
      xmg2=xmg*xmg
      alphaw=pyalem(xmg2)
      alphas=pyalps(xmg2)
      xmr=amn(i)
      xmr2=xmr*xmr
      xmq4=xmg*xm2*xmr
      xm24=(xmg2+xm2)*(xm2+xmr2)
      smin=4d0*xm2
      smax=(xmg-abs(xmr))**2
      xmqa=xmg2+2d0*xm2+xmr2
      DO 170 lin=1,nn-1
        sbar=smin+dble(lin)*(smax-smin)/dble(nn)
        grs=sbar-xmqa
        w=pylamf(xmg2,xmr2,sbar)*(0.25d0-xm2/sbar)
        w=dsqrt(w)
        xln1=log(abs((grs/2d0+xms2(1)-w)/(grs/2d0+xms2(1)+w)))
        xln2=log(abs((grs/2d0+xms2(2)-w)/(grs/2d0+xms2(2)+w)))
        b1=1d0/(grs/2d0+xms2(1)-w)-1d0/(grs/2d0+xms2(1)+w)
        b2=1d0/(grs/2d0+xms2(2)-w)-1d0/(grs/2d0+xms2(2)+w)
        g(0)=-2d0*(hl2+fl2+hr2+fr2+(hfr-hfl)*sin2d
     &  +2d0*(ff*sind2-hh*cosd2))*w
        g(1)=((hl2+fl2)*(xmqa-2d0*xms2(1)-2d0*xm*xmg*sin2d)
     &  +4d0*hfl*xm*xmr)*xln1
     &  +((hl2+fl2)*((xmqa-xms2(1))*xms2(1)-xm24
     &  +2d0*xm*xmg*(xm2+xmr2-xms2(1))*sin2d)
     &  -4d0*hfl*xmr*xm*(xmg2+xm2-xms2(1))
     &  +8d0*hfl*xmq4*sin2d)*b1
        g(2)=((hr2+fr2)*(xmqa-2d0*xms2(2)+2d0*xm*xmg*sin2d)
     &  +4d0*hfr*xmr*xm)*xln2
     &  +((hr2+fr2)*((xmqa-xms2(2))*xms2(2)-xm24
     &  +2d0*xmg*xm*sin2d*(xms2(2)-xm2-xmr2))
     &  +4d0*hfr*xm*xmr*(xms2(2)-xmg2-xm2)
     &  -8d0*hfr*xmq4*sin2d)*b2
        g(3)=(2d0*hfl*sin2d*(xms2(1)*(grs+xms2(1))+xm2*(sbar-xmg2-xmr2)
     &  +xmg2*xmr2+xm2*xm2)-2d0*xmr*xmg*(hl2*sind2+fl2*cosd2)*sbar
     &  -2d0*xmg*xm*hfl*(sbar+xmr2-xmg2)
     &  +xmr*xm*(hl2+fl2)*sin2d*(sbar+xmg2-xmr2)
     &  -4d0*xmq4*(hl2-fl2)*cos2d)/(grs+2d0*xms2(1))*xln1
        g(4)=4d0*cos2d*xm*xmg/(xms2(1)-xms2(2))*
     &  (((hlfr+hrfl)*(xm2+xmr2)+2d0*xm*xmr*(hh+ff))*(xln1-xln2)
     &  +(hlfr+hrfl)*(xms2(2)*xln2-xms2(1)*xln1))
        g(5)=(2d0*(hh*cosd2-ff*sind2)
     &  *((xms2(2)*(xms2(2)+grs)+xm2*xm2+xmg2*xmr2)*xln2
     &  +(xms2(1)*(xms2(1)+grs)+xm2*xm2+xmg2*xmr2)*xln1)
     &  +xm*((hh-ff)*sin2d*xmg-(hrfl-hlfr)*xmr)
     &  *((grs+xms2(1)*2d0)*xln1-(grs+xms2(2)*2d0)*xln2)
     &  +((hrfl-hlfr)*xmr*(sin2d*xmg*(sbar-4d0*xm2)
     &  +cos2d*xm*(sbar+xmg2-xmr2))
     &  +2d0*(ff*cosd2-hh*sind2)*xm2*(sbar-xmg2-xmr2))
     &  *(xln1+xln2))/(grs+xms2(1)+xms2(2))
        g(6)=(-2d0*hfr*sin2d*(xms2(2)*(grs+xms2(2))+xm2*(sbar-xmg2-xmr2)
     &  +xmg2*xmr2+xm2*xm2)-2d0*xmr*xmg*(hr2*sind2+fr2*cosd2)*sbar
     &  -2d0*xmg*xm*hfr*(sbar+xmr2-xmg2)
     &  -xmr*xm*(hr2+fr2)*sin2d*(sbar+xmg2-xmr2)
     &  -4d0*xmq4*(hr2-fr2)*cos2d)/(grs+2d0*xms2(2))*xln2
        summe(lin)=0d0
        DO 160 j=0,6
          summe(lin)=summe(lin)+g(j)
  160   CONTINUE
  170 CONTINUE
      summe(0)=0d0
      summe(nn)=0d0
      gam = alphaw * alphas * pysimp(summe,smin,smax,nn)
     &/ (16d0 * paru(1) * paru(102) * xmglu**3)
 
      RETURN
      END