3 C*********************************************************************
5 SUBROUTINE LUTHRU_HIJING(THR,OBL)
7 C...Purpose: to perform thrust analysis to give thrust, oblateness
8 C...and the related event axes.
9 #include "lujets_hijing.inc"
10 #include "ludat1_hijing.inc"
11 #include "ludat2_hijing.inc"
12 DIMENSION TDI(3),TPR(3)
14 C...Take copy of particles that are to be considered in thrust analysis.
18 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 100
19 IF(MSTU(41).GE.2) THEN
20 KC=LUCOMP_HIJING(K(I,2))
21 IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
23 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
26 IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN
28 $ ,'(LUTHRU_HIJING:) no more memory left in LUJETS_HIJING')
38 P(N+NP,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
40 IF(ABS(PARU(42)-1.).GT.0.001) P(N+NP,5)=P(N+NP,4)**(PARU(42)-1.)
41 PS=PS+P(N+NP,4)*P(N+NP,5)
44 C...Very low multiplicities (0 or 1) not considered.
47 $ ,'(LUTHRU_HIJING:) too few particles for analysis')
53 C...Loop over thrust and major. T axis along z direction in latter case.
57 PHI=ULANGL_HIJING(P(N+NP+1,1),P(N+NP+1,2))
58 CALL LUDBRB_HIJING(N+1,N+NP+1,0.,-PHI,0D0,0D0,0D0)
59 THE=ULANGL_HIJING(P(N+NP+1,3),P(N+NP+1,1))
60 CALL LUDBRB_HIJING(N+1,N+NP+1,-THE,0.,0D0,0D0,0D0)
63 C...Find and order particles with highest p (pT for major).
64 DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4
67 IF(ILD.EQ.2) P(I,4)=SQRT(P(I,1)**2+P(I,2)**2)
68 DO 120 ILF=N+NP+MSTU(44)+3,N+NP+4,-1
69 IF(P(I,4).LE.P(ILF,4)) GOTO 130
71 120 P(ILF+1,J)=P(ILF,J)
77 C...Find and order initial axes with highest thrust (major).
78 DO 160 ILG=N+NP+MSTU(44)+5,N+NP+MSTU(44)+15
80 NC=2**(MIN(MSTU(44),NP)-1)
84 DO 180 ILF=1,MIN(MSTU(44),NP)
86 IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
88 180 TDI(J)=TDI(J)+SGN*P(N+NP+ILF+3,J)
89 TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
90 DO 190 ILG=N+NP+MSTU(44)+MIN(ILC,10)+4,N+NP+MSTU(44)+5,-1
91 IF(TDS.LE.P(ILG,4)) GOTO 200
93 190 P(ILG+1,J)=P(ILG,J)
100 C...Iterate direction of axis until stable maximum.
107 IF(THP.LE.1E-10) TDI(J)=P(N+NP+MSTU(44)+4+ILG,J)
108 IF(THP.GT.1E-10) TDI(J)=TPR(J)
111 SGN=SIGN(P(I,5),TDI(1)*P(I,1)+TDI(2)*P(I,2)+TDI(3)*P(I,3))
113 260 TPR(J)=TPR(J)+SGN*P(I,J)
114 THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
115 IF(THP.GE.THPS+PARU(48)) GOTO 240
117 C...Save good axis. Try new initial axis until a number of tries agree.
118 IF(THP.LT.P(N+NP+ILD,4)-PARU(48).AND.ILG.LT.MIN(10,NC)) GOTO 230
119 IF(THP.GT.P(N+NP+ILD,4)+PARU(48)) THEN
121 SGN=(-1.)**INT(RLU_HIJING(0)+0.5)
123 270 P(N+NP+ILD,J)=SGN*TPR(J)/(PS*THP)
128 280 IF(IAGR.LT.MSTU(45).AND.ILG.LT.MIN(10,NC)) GOTO 230
130 C...Find minor axis and value by orthogonality.
131 SGN=(-1.)**INT(RLU_HIJING(0)+0.5)
132 P(N+NP+3,1)=-SGN*P(N+NP+2,2)
133 P(N+NP+3,2)=SGN*P(N+NP+2,1)
137 290 THP=THP+P(I,5)*ABS(P(N+NP+3,1)*P(I,1)+P(N+NP+3,2)*P(I,2))
141 C...Fill axis information. Rotate back to original coordinate system.
149 P(N+ILD,J)=P(N+NP+ILD,J)
151 CALL LUDBRB_HIJING(N+1,N+3,THE,PHI,0D0,0D0,0D0)
153 C...Select storing option. Calculate thurst and oblateness.
156 IF(MSTU(43).LE.1) MSTU(3)=3
157 IF(MSTU(43).GE.2) N=N+3
159 OBL=P(N+2,4)-P(N+3,4)