pyjmas.f

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

 
C*********************************************************************
 
C...PYJMAS
C...Determines, approximately, the two jet masses that minimize
C...the sum m_H^2 + m_L^2, a la Clavelli and Wyler.
 
      SUBROUTINE pyjmas(PMH,PML)
 
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/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pyjets/,/pydat1/,/pydat2/
C...Local arrays.
      dimension sm(3,3),sax(3),ps(3,5)
 
C...Reset.
      np=0
      DO 120 j1=1,3
        DO 100 j2=j1,3
          sm(j1,j2)=0d0
  100   CONTINUE
        DO 110 j2=1,4
          ps(j1,j2)=0d0
  110   CONTINUE
  120 CONTINUE
      pss=0d0
      pimass=pmas(pycomp(211),1)
 
C...Take copy of particles that are to be considered in mass analysis.
      DO 170 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 170
        IF(mstu(41).GE.2) THEN
          kc=pycomp(k(i,2))
          IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
     &    kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
     &    k(i,2).EQ.ksusy1+39) goto 170
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2)).EQ.0)
     &    goto 170
        ENDIF
        IF(n+np+1.GE.mstu(4)-mstu(32)-5) THEN
          CALL pyerrm(11,'(PYJMAS:) no more memory left in PYJETS')
          pmh=-2d0
          pml=-2d0
          RETURN
        ENDIF
        np=np+1
        DO 130 j=1,5
          p(n+np,j)=p(i,j)
  130   CONTINUE
        IF(mstu(42).EQ.0) p(n+np,5)=0d0
        IF(mstu(42).EQ.1.AND.k(i,2).NE.22) p(n+np,5)=pimass
        p(n+np,4)=sqrt(p(n+np,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
 
C...Fill information in sphericity tensor and total momentum vector.
        DO 150 j1=1,3
          DO 140 j2=j1,3
            sm(j1,j2)=sm(j1,j2)+p(i,j1)*p(i,j2)
  140     CONTINUE
  150   CONTINUE
        pss=pss+(p(i,1)**2+p(i,2)**2+p(i,3)**2)
        DO 160 j=1,4
          ps(3,j)=ps(3,j)+p(n+np,j)
  160   CONTINUE
  170 CONTINUE
 
C...Very low multiplicities (0 or 1) not considered.
      IF(np.LE.1) THEN
        CALL pyerrm(8,'(PYJMAS:) too few particles for analysis')
        pmh=-1d0
        pml=-1d0
        RETURN
      ENDIF
      paru(61)=sqrt(max(0d0,ps(3,4)**2-ps(3,1)**2-ps(3,2)**2-
     &ps(3,3)**2))
 
C...Find largest eigenvalue to matrix (third degree equation).
      DO 190 j1=1,3
        DO 180 j2=j1,3
          sm(j1,j2)=sm(j1,j2)/pss
  180   CONTINUE
  190 CONTINUE
      sq=(sm(1,1)*sm(2,2)+sm(1,1)*sm(3,3)+sm(2,2)*sm(3,3)-
     &sm(1,2)**2-sm(1,3)**2-sm(2,3)**2)/3d0-1d0/9d0
      sr=-0.5d0*(sq+1d0/9d0+sm(1,1)*sm(2,3)**2+sm(2,2)*sm(1,3)**2+
     &sm(3,3)*sm(1,2)**2-sm(1,1)*sm(2,2)*sm(3,3))+
     &sm(1,2)*sm(1,3)*sm(2,3)+1d0/27d0
      sp=cos(acos(max(min(sr/sqrt(-sq**3),1d0),-1d0))/3d0)
      sma=1d0/3d0+sqrt(-sq)*max(2d0*sp,sqrt(3d0*(1d0-sp**2))-sp)
 
C...Find largest eigenvector by solving equation system.
      DO 210 j1=1,3
        sm(j1,j1)=sm(j1,j1)-sma
        DO 200 j2=j1+1,3
          sm(j2,j1)=sm(j1,j2)
  200   CONTINUE
  210 CONTINUE
      smax=0d0
      DO 230 j1=1,3
        DO 220 j2=1,3
          IF(abs(sm(j1,j2)).LE.smax) goto 220
          ja=j1
          jb=j2
          smax=abs(sm(j1,j2))
  220   CONTINUE
  230 CONTINUE
      smax=0d0
      DO 250 j3=ja+1,ja+2
        j1=j3-3*((j3-1)/3)
        rl=sm(j1,jb)/sm(ja,jb)
        DO 240 j2=1,3
          sm(j1,j2)=sm(j1,j2)-rl*sm(ja,j2)
          IF(abs(sm(j1,j2)).LE.smax) goto 240
          jc=j1
          smax=abs(sm(j1,j2))
  240   CONTINUE
  250 CONTINUE
      jb1=jb+1-3*(jb/3)
      jb2=jb+2-3*((jb+1)/3)
      sax(jb1)=-sm(jc,jb2)
      sax(jb2)=sm(jc,jb1)
      sax(jb)=-(sm(ja,jb1)*sax(jb1)+sm(ja,jb2)*sax(jb2))/sm(ja,jb)
 
C...Divide particles into two initial clusters by hemisphere.
      DO 270 i=n+1,n+np
        psax=p(i,1)*sax(1)+p(i,2)*sax(2)+p(i,3)*sax(3)
        is=1
        IF(psax.LT.0d0) is=2
        k(i,3)=is
        DO 260 j=1,4
          ps(is,j)=ps(is,j)+p(i,j)
  260   CONTINUE
  270 CONTINUE
      pms=max(1d-10,ps(1,4)**2-ps(1,1)**2-ps(1,2)**2-ps(1,3)**2)+
     &max(1d-10,ps(2,4)**2-ps(2,1)**2-ps(2,2)**2-ps(2,3)**2)
 
C...Reassign one particle at a time; find maximum decrease of m^2 sum.
  280 pmd=0d0
      im=0
      DO 290 j=1,4
        ps(3,j)=ps(1,j)-ps(2,j)
  290 CONTINUE
      DO 300 i=n+1,n+np
        pps=p(i,4)*ps(3,4)-p(i,1)*ps(3,1)-p(i,2)*ps(3,2)-p(i,3)*ps(3,3)
        IF(k(i,3).EQ.1) pmdi=2d0*(p(i,5)**2-pps)
        IF(k(i,3).EQ.2) pmdi=2d0*(p(i,5)**2+pps)
        IF(pmdi.LT.pmd) THEN
          pmd=pmdi
          im=i
        ENDIF
  300 CONTINUE
 
C...Loop back if significant reduction in sum of m^2.
      IF(pmd.LT.-paru(48)*pms) THEN
        pms=pms+pmd
        is=k(im,3)
        DO 310 j=1,4
          ps(is,j)=ps(is,j)-p(im,j)
          ps(3-is,j)=ps(3-is,j)+p(im,j)
  310   CONTINUE
        k(im,3)=3-is
        goto 280
      ENDIF
 
C...Final masses and output.
      mstu(61)=n+1
      mstu(62)=np
      ps(1,5)=sqrt(max(0d0,ps(1,4)**2-ps(1,1)**2-ps(1,2)**2-ps(1,3)**2))
      ps(2,5)=sqrt(max(0d0,ps(2,4)**2-ps(2,1)**2-ps(2,2)**2-ps(2,3)**2))
      pmh=max(ps(1,5),ps(2,5))
      pml=min(ps(1,5),ps(2,5))
 
      RETURN
      END