pythiaeRHIC/pythia/pyjmas.f

0001
0002 C*********************************************************************
0003
0004 C...PYJMAS
0005 C...Determines, approximately, the two jet masses that minimize
0006 C...the sum m_H^2 + m_L^2, a la Clavelli and Wyler.
0007
0008       SUBROUTINE PYJMAS(PMH,PML)
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 SM(3,3),SAX(3),PS(3,5)
0024
0025 C...Reset.
0026       NP=0
0027       DO 120 J1=1,3
0028         DO 100 J2=J1,3
0029           SM(J1,J2)=0D0
0030   100   CONTINUE
0031         DO 110 J2=1,4
0032           PS(J1,J2)=0D0
0033   110   CONTINUE
0034   120 CONTINUE
0035       PSS=0D0
0036       PIMASS=PMAS(PYCOMP(211),1)
0037
0038 C...Take copy of particles that are to be considered in mass analysis.
0039       DO 170 I=1,N
0040         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 170
0041         IF(MSTU(41).GE.2) THEN
0042           KC=PYCOMP(K(I,2))
0043           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0044      &    KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0045      &    K(I,2).EQ.KSUSY1+39) GOTO 170
0046           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0)
0047      &    GOTO 170
0048         ENDIF
0049         IF(N+NP+1.GE.MSTU(4)-MSTU(32)-5) THEN
0050           CALL PYERRM(11,'(PYJMAS:) no more memory left in PYJETS')
0051           PMH=-2D0
0052           PML=-2D0
0053           RETURN
0054         ENDIF
0055         NP=NP+1
0056         DO 130 J=1,5
0057           P(N+NP,J)=P(I,J)
0058   130   CONTINUE
0059         IF(MSTU(42).EQ.0) P(N+NP,5)=0D0
0060         IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PIMASS
0061         P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0062
0063 C...Fill information in sphericity tensor and total momentum vector.
0064         DO 150 J1=1,3
0065           DO 140 J2=J1,3
0066             SM(J1,J2)=SM(J1,J2)+P(I,J1)*P(I,J2)
0067   140     CONTINUE
0068   150   CONTINUE
0069         PSS=PSS+(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0070         DO 160 J=1,4
0071           PS(3,J)=PS(3,J)+P(N+NP,J)
0072   160   CONTINUE
0073   170 CONTINUE
0074
0075 C...Very low multiplicities (0 or 1) not considered.
0076       IF(NP.LE.1) THEN
0077         CALL PYERRM(8,'(PYJMAS:) too few particles for analysis')
0078         PMH=-1D0
0079         PML=-1D0
0080         RETURN
0081       ENDIF
0082       PARU(61)=SQRT(MAX(0D0,PS(3,4)**2-PS(3,1)**2-PS(3,2)**2-
0083      &PS(3,3)**2))
0084
0085 C...Find largest eigenvalue to matrix (third degree equation).
0086       DO 190 J1=1,3
0087         DO 180 J2=J1,3
0088           SM(J1,J2)=SM(J1,J2)/PSS
0089   180   CONTINUE
0090   190 CONTINUE
0091       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-
0092      &SM(1,2)**2-SM(1,3)**2-SM(2,3)**2)/3D0-1D0/9D0
0093       SR=-0.5D0*(SQ+1D0/9D0+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+
0094      &SM(3,3)*SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+
0095      &SM(1,2)*SM(1,3)*SM(2,3)+1D0/27D0
0096       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1D0),-1D0))/3D0)
0097       SMA=1D0/3D0+SQRT(-SQ)*MAX(2D0*SP,SQRT(3D0*(1D0-SP**2))-SP)
0098
0099 C...Find largest eigenvector by solving equation system.
0100       DO 210 J1=1,3
0101         SM(J1,J1)=SM(J1,J1)-SMA
0102         DO 200 J2=J1+1,3
0103           SM(J2,J1)=SM(J1,J2)
0104   200   CONTINUE
0105   210 CONTINUE
0106       SMAX=0D0
0107       DO 230 J1=1,3
0108         DO 220 J2=1,3
0109           IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 220
0110           JA=J1
0111           JB=J2
0112           SMAX=ABS(SM(J1,J2))
0113   220   CONTINUE
0114   230 CONTINUE
0115       SMAX=0D0
0116       DO 250 J3=JA+1,JA+2
0117         J1=J3-3*((J3-1)/3)
0118         RL=SM(J1,JB)/SM(JA,JB)
0119         DO 240 J2=1,3
0120           SM(J1,J2)=SM(J1,J2)-RL*SM(JA,J2)
0121           IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 240
0122           JC=J1
0123           SMAX=ABS(SM(J1,J2))
0124   240   CONTINUE
0125   250 CONTINUE
0126       JB1=JB+1-3*(JB/3)
0127       JB2=JB+2-3*((JB+1)/3)
0128       SAX(JB1)=-SM(JC,JB2)
0129       SAX(JB2)=SM(JC,JB1)
0130       SAX(JB)=-(SM(JA,JB1)*SAX(JB1)+SM(JA,JB2)*SAX(JB2))/SM(JA,JB)
0131
0132 C...Divide particles into two initial clusters by hemisphere.
0133       DO 270 I=N+1,N+NP
0134         PSAX=P(I,1)*SAX(1)+P(I,2)*SAX(2)+P(I,3)*SAX(3)
0135         IS=1
0136         IF(PSAX.LT.0D0) IS=2
0137         K(I,3)=IS
0138         DO 260 J=1,4
0139           PS(IS,J)=PS(IS,J)+P(I,J)
0140   260   CONTINUE
0141   270 CONTINUE
0142       PMS=MAX(1D-10,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2)+
0143      &MAX(1D-10,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2)
0144
0145 C...Reassign one particle at a time; find maximum decrease of m^2 sum.
0146   280 PMD=0D0
0147       IM=0
0148       DO 290 J=1,4
0149         PS(3,J)=PS(1,J)-PS(2,J)
0150   290 CONTINUE
0151       DO 300 I=N+1,N+NP
0152         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)
0153         IF(K(I,3).EQ.1) PMDI=2D0*(P(I,5)**2-PPS)
0154         IF(K(I,3).EQ.2) PMDI=2D0*(P(I,5)**2+PPS)
0155         IF(PMDI.LT.PMD) THEN
0156           PMD=PMDI
0157           IM=I
0158         ENDIF
0159   300 CONTINUE
0160
0161 C...Loop back if significant reduction in sum of m^2.
0162       IF(PMD.LT.-PARU(48)*PMS) THEN
0163         PMS=PMS+PMD
0164         IS=K(IM,3)
0165         DO 310 J=1,4
0166           PS(IS,J)=PS(IS,J)-P(IM,J)
0167           PS(3-IS,J)=PS(3-IS,J)+P(IM,J)
0168   310   CONTINUE
0169         K(IM,3)=3-IS
0170         GOTO 280
0171       ENDIF
0172
0173 C...Final masses and output.
0174       MSTU(61)=N+1
0175       MSTU(62)=NP
0176       PS(1,5)=SQRT(MAX(0D0,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2))
0177       PS(2,5)=SQRT(MAX(0D0,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2))
0178       PMH=MAX(PS(1,5),PS(2,5))
0179       PML=MIN(PS(1,5),PS(2,5))
0180
0181       RETURN
0182       END