pytecm.f

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

 
C*********************************************************************
 
C...PYTECM
C...Finds the s-hat dependent eigenvalues of the inverse propagator
C...matrix for gamma, Z, techni-rho, and techni-omega to optimize the
C...phase space generation.
 
      SUBROUTINE pytecm(S1,S2)
 
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/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pytcsm/itcm(0:99),rtcm(0:99)
      SAVE /pydat1/,/pydat2/,/pypars/,/pytcsm/
 
C...Local variables.
      DOUBLE PRECISION ar(4,4),wr(4),zr(4,4),zi(4,4),work(12,12),
     &at(4,4),wi(4),fv1(4),fv2(4),fv3(4),sh,aem,tanw,ct2w,qupd,alprht,
     &far,fao,fzr,fzo,shr,r1,r2,s1,s2,wdtp(0:400),wdte(0:400,0:5)
      INTEGER i,j,ierr
 
      sh=pmas(pycomp(ktechn+113),1)**2
      aem=pyalem(sh)
 
      tanw=sqrt(paru(102)/(1d0-paru(102)))
      ct2w=(1d0-2d0*paru(102))/(2d0*paru(102)/tanw)
      qupd=2d0*rtcm(2)-1d0
 
      alprht=2.91d0*(3d0/dble(itcm(1)))
      far=sqrt(aem/alprht)
      fao=far*qupd
      fzr=far*ct2w
      fzo=-fao*tanw
 
      ar(1,1) = sh
      ar(2,2) = sh-pmas(23,1)**2
      ar(3,3) = sh-pmas(pycomp(ktechn+113),1)**2
      ar(4,4) = sh-pmas(pycomp(ktechn+223),1)**2
      ar(1,2) = 0d0
      ar(2,1) = 0d0
      ar(1,3) = -sh*far
      ar(3,1) = ar(1,3)
      ar(1,4) = -sh*fao
      ar(4,1) = ar(1,4)
      ar(2,3) = -sh*fzr
      ar(3,2) = ar(2,3)
      ar(2,4) = -sh*fzo
      ar(4,2) = ar(2,4)
      ar(3,4) = 0d0
      ar(4,3) = 0d0
CCCCCCCC
      DO 110 i=1,4
        DO 100 j=1,4
          at(i,j)=0d0
  100   CONTINUE
  110 CONTINUE
      shr=sqrt(sh)
      CALL pywidt(23,sh,wdtp,wdte)
      at(2,2) = wdtp(0)*shr
      CALL pywidt(ktechn+113,sh,wdtp,wdte)
      at(3,3) = wdtp(0)*shr
      CALL pywidt(ktechn+223,sh,wdtp,wdte)
      at(4,4) = wdtp(0)*shr
CCCC
      CALL pyeicg(4,4,ar,at,wr,wi,0,zr,zi,fv1,fv2,fv3,ierr)
      DO 120 i=1,4
        wi(i)=sqrt(abs(sh-wr(i)))
        wr(i)=abs(wr(i))
  120 CONTINUE
      r1=min(wr(1),wr(2),wr(3),wr(4))
      r2=1d20
      s1=0d0
      s2=0d0
      DO 130 i=1,4
        IF(abs(wr(i)-r1).LT.1d-6) THEN
          s1=wi(i)
          goto 130
        ENDIF
        IF(wr(i).LE.r2) THEN
          r2=wr(i)
          s2=wi(i)
        ENDIF
  130 CONTINUE
      s1=s1**2
      s2=s2**2
      RETURN
      END