--- /dev/null
+
+C*********************************************************************
+
+ SUBROUTINE LUTHRU(THR,OBL)
+
+C...Purpose: to perform thrust analysis to give thrust, oblateness
+C...and the related event axes.
+ COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
+ COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
+ COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
+ SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
+ DIMENSION TDI(3),TPR(3)
+
+C...Take copy of particles that are to be considered in thrust analysis.
+ NP=0
+ PS=0.
+ 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=LUCOMP(K(I,2))
+ IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
+ & KC.EQ.18) GOTO 100
+ IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
+ & GOTO 100
+ ENDIF
+ IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN
+ CALL LUERRM(11,'(LUTHRU:) no more memory left in LUJETS')
+ THR=-2.
+ OBL=-2.
+ 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)=1.
+ IF(ABS(PARU(42)-1.).GT.0.001) P(N+NP,5)=P(N+NP,4)**(PARU(42)-1.)
+ 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 LUERRM(8,'(LUTHRU:) too few particles for analysis')
+ THR=-1.
+ OBL=-1.
+ 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=ULANGL(P(N+NP+1,1),P(N+NP+1,2))
+ MSTU(33)=1
+ CALL LUDBRB(N+1,N+NP+1,0.,-PHI,0D0,0D0,0D0)
+ THE=ULANGL(P(N+NP+1,3),P(N+NP+1,1))
+ CALL LUDBRB(N+1,N+NP+1,-THE,0.,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)=0.
+ 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)=0.
+ 170 CONTINUE
+ NC=2**(MIN(MSTU(44),NP)-1)
+ DO 250 ILC=1,NC
+ DO 180 J=1,3
+ TDI(J)=0.
+ 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)=0.
+ ILG=0
+ 260 ILG=ILG+1
+ THP=0.
+ 270 THPS=THP
+ DO 280 J=1,3
+ IF(THP.LE.1E-10) TDI(J)=P(N+NP+MSTU(44)+4+ILG,J)
+ IF(THP.GT.1E-10) TDI(J)=TPR(J)
+ TPR(J)=0.
+ 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=(-1.)**INT(RLU(0)+0.5)
+ 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)=0.
+ 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=(-1.)**INT(RLU(0)+0.5)
+ 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)=0.
+ THP=0.
+ 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)=0.
+
+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)=0.
+ 340 CONTINUE
+ 350 CONTINUE
+ CALL LUDBRB(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