pythiaeRHIC/pythia/pythru.f

0001
0002 C*********************************************************************
0003
0004 C...PYTHRU
0005 C...Performs thrust analysis to give thrust, oblateness
0006 C...and the related event axes.
0007
0008       SUBROUTINE PYTHRU(THR,OBL)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Parameter statement to help give large particle numbers.
0015       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0016      &KEXCIT=4000000,KDIMEN=5000000)
0017 C...Commonblocks.
0018       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0022 C...Local arrays.
0023       DIMENSION TDI(3),TPR(3)
0024
0025 C...Take copy of particles that are to be considered in thrust analysis.
0026       NP=0
0027       PS=0D0
0028       DO 100 I=1,N
0029         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 100
0030         IF(MSTU(41).GE.2) THEN
0031           KC=PYCOMP(K(I,2))
0032           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0033      &    KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0034      &    K(I,2).EQ.KSUSY1+39) GOTO 100
0035           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0)
0036      &    GOTO 100
0037         ENDIF
0038         IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN
0039           CALL PYERRM(11,'(PYTHRU:) no more memory left in PYJETS')
0040           THR=-2D0
0041           OBL=-2D0
0042           RETURN
0043         ENDIF
0044         NP=NP+1
0045         K(N+NP,1)=23
0046         P(N+NP,1)=P(I,1)
0047         P(N+NP,2)=P(I,2)
0048         P(N+NP,3)=P(I,3)
0049         P(N+NP,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0050         P(N+NP,5)=1D0
0051         IF(ABS(PARU(42)-1D0).GT.0.001D0) P(N+NP,5)=
0052      &  P(N+NP,4)**(PARU(42)-1D0)
0053         PS=PS+P(N+NP,4)*P(N+NP,5)
0054   100 CONTINUE
0055
0056 C...Very low multiplicities (0 or 1) not considered.
0057       IF(NP.LE.1) THEN
0058         CALL PYERRM(8,'(PYTHRU:) too few particles for analysis')
0059         THR=-1D0
0060         OBL=-1D0
0061         RETURN
0062       ENDIF
0063
0064 C...Loop over thrust and major. T axis along z direction in latter case.
0065       DO 320 ILD=1,2
0066         IF(ILD.EQ.2) THEN
0067           K(N+NP+1,1)=31
0068           PHI=PYANGL(P(N+NP+1,1),P(N+NP+1,2))
0069           MSTU(33)=1
0070           CALL PYROBO(N+1,N+NP+1,0D0,-PHI,0D0,0D0,0D0)
0071           THE=PYANGL(P(N+NP+1,3),P(N+NP+1,1))
0072           CALL PYROBO(N+1,N+NP+1,-THE,0D0,0D0,0D0,0D0)
0073         ENDIF
0074
0075 C...Find and order particles with highest p (pT for major).
0076         DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4
0077           P(ILF,4)=0D0
0078   110   CONTINUE
0079         DO 160 I=N+1,N+NP
0080           IF(ILD.EQ.2) P(I,4)=SQRT(P(I,1)**2+P(I,2)**2)
0081           DO 130 ILF=N+NP+MSTU(44)+3,N+NP+4,-1
0082             IF(P(I,4).LE.P(ILF,4)) GOTO 140
0083             DO 120 J=1,5
0084               P(ILF+1,J)=P(ILF,J)
0085   120       CONTINUE
0086   130     CONTINUE
0087           ILF=N+NP+3
0088   140     DO 150 J=1,5
0089             P(ILF+1,J)=P(I,J)
0090   150     CONTINUE
0091   160   CONTINUE
0092
0093 C...Find and order initial axes with highest thrust (major).
0094         DO 170 ILG=N+NP+MSTU(44)+5,N+NP+MSTU(44)+15
0095           P(ILG,4)=0D0
0096   170   CONTINUE
0097         NC=2**(MIN(MSTU(44),NP)-1)
0098         DO 250 ILC=1,NC
0099           DO 180 J=1,3
0100             TDI(J)=0D0
0101   180     CONTINUE
0102           DO 200 ILF=1,MIN(MSTU(44),NP)
0103             SGN=P(N+NP+ILF+3,5)
0104             IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
0105             DO 190 J=1,4-ILD
0106               TDI(J)=TDI(J)+SGN*P(N+NP+ILF+3,J)
0107   190       CONTINUE
0108   200     CONTINUE
0109           TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
0110           DO 220 ILG=N+NP+MSTU(44)+MIN(ILC,10)+4,N+NP+MSTU(44)+5,-1
0111             IF(TDS.LE.P(ILG,4)) GOTO 230
0112             DO 210 J=1,4
0113               P(ILG+1,J)=P(ILG,J)
0114   210       CONTINUE
0115   220     CONTINUE
0116           ILG=N+NP+MSTU(44)+4
0117   230     DO 240 J=1,3
0118             P(ILG+1,J)=TDI(J)
0119   240     CONTINUE
0120           P(ILG+1,4)=TDS
0121   250   CONTINUE
0122
0123 C...Iterate direction of axis until stable maximum.
0124         P(N+NP+ILD,4)=0D0
0125         ILG=0
0126   260   ILG=ILG+1
0127         THP=0D0
0128   270   THPS=THP
0129         DO 280 J=1,3
0130           IF(THP.LE.1D-10) TDI(J)=P(N+NP+MSTU(44)+4+ILG,J)
0131           IF(THP.GT.1D-10) TDI(J)=TPR(J)
0132           TPR(J)=0D0
0133   280   CONTINUE
0134         DO 300 I=N+1,N+NP
0135           SGN=SIGN(P(I,5),TDI(1)*P(I,1)+TDI(2)*P(I,2)+TDI(3)*P(I,3))
0136           DO 290 J=1,4-ILD
0137             TPR(J)=TPR(J)+SGN*P(I,J)
0138   290     CONTINUE
0139   300   CONTINUE
0140         THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
0141         IF(THP.GE.THPS+PARU(48)) GOTO 270
0142
0143 C...Save good axis. Try new initial axis until a number of tries agree.
0144         IF(THP.LT.P(N+NP+ILD,4)-PARU(48).AND.ILG.LT.MIN(10,NC)) GOTO 260
0145         IF(THP.GT.P(N+NP+ILD,4)+PARU(48)) THEN
0146           IAGR=0
0147           SGN=(-1D0)**INT(PYR(0)+0.5D0)
0148           DO 310 J=1,3
0149             P(N+NP+ILD,J)=SGN*TPR(J)/(PS*THP)
0150   310     CONTINUE
0151           P(N+NP+ILD,4)=THP
0152           P(N+NP+ILD,5)=0D0
0153         ENDIF
0154         IAGR=IAGR+1
0155         IF(IAGR.LT.MSTU(45).AND.ILG.LT.MIN(10,NC)) GOTO 260
0156   320 CONTINUE
0157
0158 C...Find minor axis and value by orthogonality.
0159       SGN=(-1D0)**INT(PYR(0)+0.5D0)
0160       P(N+NP+3,1)=-SGN*P(N+NP+2,2)
0161       P(N+NP+3,2)=SGN*P(N+NP+2,1)
0162       P(N+NP+3,3)=0D0
0163       THP=0D0
0164       DO 330 I=N+1,N+NP
0165         THP=THP+P(I,5)*ABS(P(N+NP+3,1)*P(I,1)+P(N+NP+3,2)*P(I,2))
0166   330 CONTINUE
0167       P(N+NP+3,4)=THP/PS
0168       P(N+NP+3,5)=0D0
0169
0170 C...Fill axis information. Rotate back to original coordinate system.
0171       DO 350 ILD=1,3
0172         K(N+ILD,1)=31
0173         K(N+ILD,2)=96
0174         K(N+ILD,3)=ILD
0175         K(N+ILD,4)=0
0176         K(N+ILD,5)=0
0177         DO 340 J=1,5
0178           P(N+ILD,J)=P(N+NP+ILD,J)
0179           V(N+ILD,J)=0D0
0180   340   CONTINUE
0181   350 CONTINUE
0182       CALL PYROBO(N+1,N+3,THE,PHI,0D0,0D0,0D0)
0183
0184 C...Calculate thrust and oblateness. Select storing option.
0185       THR=P(N+1,4)
0186       OBL=P(N+2,4)-P(N+3,4)
0187       MSTU(61)=N+1
0188       MSTU(62)=NP
0189       IF(MSTU(43).LE.1) MSTU(3)=3
0190       IF(MSTU(43).GE.2) N=N+3
0191
0192       RETURN
0193       END