hijing/src/luthru.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUTHRU(THR,OBL)
0005
0006 C...Purpose: to perform thrust analysis to give thrust, oblateness
0007 C...and the related event axes.
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       DIMENSION TDI(3),TPR(3)
0015
0016 C...Take copy of particles that are to be considered in thrust analysis.
0017       NP=0
0018       PS=0.
0019       DO 100 I=1,N
0020       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 100
0021       IF(MSTU(41).GE.2) THEN
0022         KC=LUCOMP(K(I,2))
0023         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0024      &  KC.EQ.18) GOTO 100
0025         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0026      &  GOTO 100
0027       ENDIF
0028       IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN
0029         CALL LUERRM(11,'(LUTHRU:) no more memory left in LUJETS')
0030         THR=-2.
0031         OBL=-2.
0032         RETURN
0033       ENDIF
0034       NP=NP+1
0035       K(N+NP,1)=23
0036       P(N+NP,1)=P(I,1)
0037       P(N+NP,2)=P(I,2)
0038       P(N+NP,3)=P(I,3)
0039       P(N+NP,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0040       P(N+NP,5)=1.
0041       IF(ABS(PARU(42)-1.).GT.0.001) P(N+NP,5)=P(N+NP,4)**(PARU(42)-1.)
0042       PS=PS+P(N+NP,4)*P(N+NP,5)
0043   100 CONTINUE
0044
0045 C...Very low multiplicities (0 or 1) not considered.
0046       IF(NP.LE.1) THEN
0047         CALL LUERRM(8,'(LUTHRU:) too few particles for analysis')
0048         THR=-1.
0049         OBL=-1.
0050         RETURN
0051       ENDIF
0052
0053 C...Loop over thrust and major. T axis along z direction in latter case.
0054       DO 280 ILD=1,2
0055       IF(ILD.EQ.2) THEN
0056         K(N+NP+1,1)=31
0057         PHI=ULANGL(P(N+NP+1,1),P(N+NP+1,2))
0058         CALL LUDBRB(N+1,N+NP+1,0.,-PHI,0D0,0D0,0D0)
0059         THE=ULANGL(P(N+NP+1,3),P(N+NP+1,1))
0060         CALL LUDBRB(N+1,N+NP+1,-THE,0.,0D0,0D0,0D0)
0061       ENDIF
0062
0063 C...Find and order particles with highest p (pT for major).
0064       DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4
0065   110 P(ILF,4)=0.
0066       DO 150 I=N+1,N+NP
0067       IF(ILD.EQ.2) P(I,4)=SQRT(P(I,1)**2+P(I,2)**2)
0068       DO 120 ILF=N+NP+MSTU(44)+3,N+NP+4,-1
0069       IF(P(I,4).LE.P(ILF,4)) GOTO 130
0070       DO 120 J=1,5
0071   120 P(ILF+1,J)=P(ILF,J)
0072       ILF=N+NP+3
0073   130 DO 140 J=1,5
0074   140 P(ILF+1,J)=P(I,J)
0075   150 CONTINUE
0076
0077 C...Find and order initial axes with highest thrust (major).
0078       DO 160 ILG=N+NP+MSTU(44)+5,N+NP+MSTU(44)+15
0079   160 P(ILG,4)=0.
0080       NC=2**(MIN(MSTU(44),NP)-1)
0081       DO 220 ILC=1,NC
0082       DO 170 J=1,3
0083   170 TDI(J)=0.
0084       DO 180 ILF=1,MIN(MSTU(44),NP)
0085       SGN=P(N+NP+ILF+3,5)
0086       IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
0087       DO 180 J=1,4-ILD
0088   180 TDI(J)=TDI(J)+SGN*P(N+NP+ILF+3,J)
0089       TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
0090       DO 190 ILG=N+NP+MSTU(44)+MIN(ILC,10)+4,N+NP+MSTU(44)+5,-1
0091       IF(TDS.LE.P(ILG,4)) GOTO 200
0092       DO 190 J=1,4
0093   190 P(ILG+1,J)=P(ILG,J)
0094       ILG=N+NP+MSTU(44)+4
0095   200 DO 210 J=1,3
0096   210 P(ILG+1,J)=TDI(J)
0097       P(ILG+1,4)=TDS
0098   220 CONTINUE
0099
0100 C...Iterate direction of axis until stable maximum.
0101       P(N+NP+ILD,4)=0.
0102       ILG=0
0103   230 ILG=ILG+1
0104       THP=0.
0105   240 THPS=THP
0106       DO 250 J=1,3
0107       IF(THP.LE.1E-10) TDI(J)=P(N+NP+MSTU(44)+4+ILG,J)
0108       IF(THP.GT.1E-10) TDI(J)=TPR(J)
0109   250 TPR(J)=0.
0110       DO 260 I=N+1,N+NP
0111       SGN=SIGN(P(I,5),TDI(1)*P(I,1)+TDI(2)*P(I,2)+TDI(3)*P(I,3))
0112       DO 260 J=1,4-ILD
0113   260 TPR(J)=TPR(J)+SGN*P(I,J)
0114       THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
0115       IF(THP.GE.THPS+PARU(48)) GOTO 240
0116
0117 C...Save good axis. Try new initial axis until a number of tries agree.
0118       IF(THP.LT.P(N+NP+ILD,4)-PARU(48).AND.ILG.LT.MIN(10,NC)) GOTO 230
0119       IF(THP.GT.P(N+NP+ILD,4)+PARU(48)) THEN
0120         IAGR=0
0121         SGN=(-1.)**INT(RLU(0)+0.5)
0122         DO 270 J=1,3
0123   270   P(N+NP+ILD,J)=SGN*TPR(J)/(PS*THP)
0124         P(N+NP+ILD,4)=THP
0125         P(N+NP+ILD,5)=0.
0126       ENDIF
0127       IAGR=IAGR+1
0128   280 IF(IAGR.LT.MSTU(45).AND.ILG.LT.MIN(10,NC)) GOTO 230
0129
0130 C...Find minor axis and value by orthogonality.
0131       SGN=(-1.)**INT(RLU(0)+0.5)
0132       P(N+NP+3,1)=-SGN*P(N+NP+2,2)
0133       P(N+NP+3,2)=SGN*P(N+NP+2,1)
0134       P(N+NP+3,3)=0.
0135       THP=0.
0136       DO 290 I=N+1,N+NP
0137   290 THP=THP+P(I,5)*ABS(P(N+NP+3,1)*P(I,1)+P(N+NP+3,2)*P(I,2))
0138       P(N+NP+3,4)=THP/PS
0139       P(N+NP+3,5)=0.
0140
0141 C...Fill axis information. Rotate back to original coordinate system.
0142       DO 300 ILD=1,3
0143       K(N+ILD,1)=31
0144       K(N+ILD,2)=96
0145       K(N+ILD,3)=ILD
0146       K(N+ILD,4)=0
0147       K(N+ILD,5)=0
0148       DO 300 J=1,5
0149       P(N+ILD,J)=P(N+NP+ILD,J)
0150   300 V(N+ILD,J)=0.
0151       CALL LUDBRB(N+1,N+3,THE,PHI,0D0,0D0,0D0)
0152
0153 C...Select storing option. Calculate thurst and oblateness.
0154       MSTU(61)=N+1
0155       MSTU(62)=NP
0156       IF(MSTU(43).LE.1) MSTU(3)=3
0157       IF(MSTU(43).GE.2) N=N+3
0158       THR=P(N+1,4)
0159       OBL=P(N+2,4)-P(N+3,4)
0160
0161       RETURN
0162       END