pythru.f

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

 
C*********************************************************************
 
C...PYTHRU
C...Performs thrust analysis to give thrust, oblateness
C...and the related event axes.
 
      SUBROUTINE pythru(THR,OBL)
 
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 tdi(3),tpr(3)
 
C...Take copy of particles that are to be considered in thrust analysis.
      np=0
      ps=0d0
      DO 100 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 100
        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 100
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2)).EQ.0)
     &    goto 100
        ENDIF
        IF(n+np+mstu(44)+15.GE.mstu(4)-mstu(32)-5) THEN
          CALL pyerrm(11,'(PYTHRU:) no more memory left in PYJETS')
          thr=-2d0
          obl=-2d0
          RETURN
        ENDIF
        np=np+1
        k(n+np,1)=23
        p(n+np,1)=p(i,1)
        p(n+np,2)=p(i,2)
        p(n+np,3)=p(i,3)
        p(n+np,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
        p(n+np,5)=1d0
        IF(abs(paru(42)-1d0).GT.0.001d0) p(n+np,5)=
     &  p(n+np,4)**(paru(42)-1d0)
        ps=ps+p(n+np,4)*p(n+np,5)
  100 CONTINUE
 
C...Very low multiplicities (0 or 1) not considered.
      IF(np.LE.1) THEN
        CALL pyerrm(8,'(PYTHRU:) too few particles for analysis')
        thr=-1d0
        obl=-1d0
        RETURN
      ENDIF
 
C...Loop over thrust and major. T axis along z direction in latter case.
      DO 320 ild=1,2
        IF(ild.EQ.2) THEN
          k(n+np+1,1)=31
          phi=pyangl(p(n+np+1,1),p(n+np+1,2))
          mstu(33)=1
          CALL pyrobo(n+1,n+np+1,0d0,-phi,0d0,0d0,0d0)
          the=pyangl(p(n+np+1,3),p(n+np+1,1))
          CALL pyrobo(n+1,n+np+1,-the,0d0,0d0,0d0,0d0)
        ENDIF
 
C...Find and order particles with highest p (pT for major).
        DO 110 ilf=n+np+4,n+np+mstu(44)+4
          p(ilf,4)=0d0
  110   CONTINUE
        DO 160 i=n+1,n+np
          IF(ild.EQ.2) p(i,4)=sqrt(p(i,1)**2+p(i,2)**2)
          DO 130 ilf=n+np+mstu(44)+3,n+np+4,-1
            IF(p(i,4).LE.p(ilf,4)) goto 140
            DO 120 j=1,5
              p(ilf+1,j)=p(ilf,j)
  120       CONTINUE
  130     CONTINUE
          ilf=n+np+3
  140     DO 150 j=1,5
            p(ilf+1,j)=p(i,j)
  150     CONTINUE
  160   CONTINUE
 
C...Find and order initial axes with highest thrust (major).
        DO 170 ilg=n+np+mstu(44)+5,n+np+mstu(44)+15
          p(ilg,4)=0d0
  170   CONTINUE
        nc=2**(min(mstu(44),np)-1)
        DO 250 ilc=1,nc
          DO 180 j=1,3
            tdi(j)=0d0
  180     CONTINUE
          DO 200 ilf=1,min(mstu(44),np)
            sgn=p(n+np+ilf+3,5)
            IF(2**ilf*((ilc+2**(ilf-1)-1)/2**ilf).GE.ilc) sgn=-sgn
            DO 190 j=1,4-ild
              tdi(j)=tdi(j)+sgn*p(n+np+ilf+3,j)
  190       CONTINUE
  200     CONTINUE
          tds=tdi(1)**2+tdi(2)**2+tdi(3)**2
          DO 220 ilg=n+np+mstu(44)+min(ilc,10)+4,n+np+mstu(44)+5,-1
            IF(tds.LE.p(ilg,4)) goto 230
            DO 210 j=1,4
              p(ilg+1,j)=p(ilg,j)
  210       CONTINUE
  220     CONTINUE
          ilg=n+np+mstu(44)+4
  230     DO 240 j=1,3
            p(ilg+1,j)=tdi(j)
  240     CONTINUE
          p(ilg+1,4)=tds
  250   CONTINUE
 
C...Iterate direction of axis until stable maximum.
        p(n+np+ild,4)=0d0
        ilg=0
  260   ilg=ilg+1
        thp=0d0
  270   thps=thp
        DO 280 j=1,3
          IF(thp.LE.1d-10) tdi(j)=p(n+np+mstu(44)+4+ilg,j)
          IF(thp.GT.1d-10) tdi(j)=tpr(j)
          tpr(j)=0d0
  280   CONTINUE
        DO 300 i=n+1,n+np
          sgn=sign(p(i,5),tdi(1)*p(i,1)+tdi(2)*p(i,2)+tdi(3)*p(i,3))
          DO 290 j=1,4-ild
            tpr(j)=tpr(j)+sgn*p(i,j)
  290     CONTINUE
  300   CONTINUE
        thp=sqrt(tpr(1)**2+tpr(2)**2+tpr(3)**2)/ps
        IF(thp.GE.thps+paru(48)) goto 270
 
C...Save good axis. Try new initial axis until a number of tries agree.
        IF(thp.LT.p(n+np+ild,4)-paru(48).AND.ilg.LT.min(10,nc)) goto 260
        IF(thp.GT.p(n+np+ild,4)+paru(48)) THEN
          iagr=0
          sgn=(-1d0)**int(pyr(0)+0.5d0)
          DO 310 j=1,3
            p(n+np+ild,j)=sgn*tpr(j)/(ps*thp)
  310     CONTINUE
          p(n+np+ild,4)=thp
          p(n+np+ild,5)=0d0
        ENDIF
        iagr=iagr+1
        IF(iagr.LT.mstu(45).AND.ilg.LT.min(10,nc)) goto 260
  320 CONTINUE
 
C...Find minor axis and value by orthogonality.
      sgn=(-1d0)**int(pyr(0)+0.5d0)
      p(n+np+3,1)=-sgn*p(n+np+2,2)
      p(n+np+3,2)=sgn*p(n+np+2,1)
      p(n+np+3,3)=0d0
      thp=0d0
      DO 330 i=n+1,n+np
        thp=thp+p(i,5)*abs(p(n+np+3,1)*p(i,1)+p(n+np+3,2)*p(i,2))
  330 CONTINUE
      p(n+np+3,4)=thp/ps
      p(n+np+3,5)=0d0
 
C...Fill axis information. Rotate back to original coordinate system.
      DO 350 ild=1,3
        k(n+ild,1)=31
        k(n+ild,2)=96
        k(n+ild,3)=ild
        k(n+ild,4)=0
        k(n+ild,5)=0
        DO 340 j=1,5
          p(n+ild,j)=p(n+np+ild,j)
          v(n+ild,j)=0d0
  340   CONTINUE
  350 CONTINUE
      CALL pyrobo(n+1,n+3,the,phi,0d0,0d0,0d0)
 
C...Calculate thrust and oblateness. Select storing option.
      thr=p(n+1,4)
      obl=p(n+2,4)-p(n+3,4)
      mstu(61)=n+1
      mstu(62)=np
      IF(mstu(43).LE.1) mstu(3)=3
      IF(mstu(43).GE.2) n=n+3
 
      RETURN
      END