Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luthru_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE LUTHRU_HIJING(THR,OBL)    
6     
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)   
13     
14 C...Take copy of particles that are to be considered in thrust analysis.    
15       NP=0  
16       PS=0. 
17       DO 100 I=1,N  
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.  
22      &  KC.EQ.18) GOTO 100  
23         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
24      $       .EQ.0)GOTO 100    
25       ENDIF 
26       IF(N+NP+MSTU(44)+15.GE.MSTU(4)-MSTU(32)-5) THEN   
27          CALL LUERRM_HIJING(11
28      $        ,'(LUTHRU_HIJING:) no more memory left in LUJETS_HIJING')   
29         THR=-2. 
30         OBL=-2. 
31         RETURN  
32       ENDIF 
33       NP=NP+1   
34       K(N+NP,1)=23  
35       P(N+NP,1)=P(I,1)  
36       P(N+NP,2)=P(I,2)  
37       P(N+NP,3)=P(I,3)  
38       P(N+NP,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) 
39       P(N+NP,5)=1.  
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) 
42   100 CONTINUE  
43     
44 C...Very low multiplicities (0 or 1) not considered.    
45       IF(NP.LE.1) THEN  
46          CALL LUERRM_HIJING(8
47      $        ,'(LUTHRU_HIJING:) too few particles for analysis')   
48         THR=-1. 
49         OBL=-1. 
50         RETURN  
51       ENDIF 
52     
53 C...Loop over thrust and major. T axis along z direction in latter case.    
54       DO 280 ILD=1,2    
55       IF(ILD.EQ.2) THEN 
56         K(N+NP+1,1)=31  
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) 
61       ENDIF 
62     
63 C...Find and order particles with highest p (pT for major). 
64       DO 110 ILF=N+NP+4,N+NP+MSTU(44)+4 
65   110 P(ILF,4)=0.   
66       DO 150 I=N+1,N+NP 
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   
70       DO 120 J=1,5  
71   120 P(ILF+1,J)=P(ILF,J)   
72       ILF=N+NP+3    
73   130 DO 140 J=1,5  
74   140 P(ILF+1,J)=P(I,J) 
75   150 CONTINUE  
76     
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   
79   160 P(ILG,4)=0.   
80       NC=2**(MIN(MSTU(44),NP)-1)    
81       DO 220 ILC=1,NC   
82       DO 170 J=1,3  
83   170 TDI(J)=0. 
84       DO 180 ILF=1,MIN(MSTU(44),NP) 
85       SGN=P(N+NP+ILF+3,5)   
86       IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN    
87       DO 180 J=1,4-ILD  
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  
92       DO 190 J=1,4  
93   190 P(ILG+1,J)=P(ILG,J)   
94       ILG=N+NP+MSTU(44)+4   
95   200 DO 210 J=1,3  
96   210 P(ILG+1,J)=TDI(J) 
97       P(ILG+1,4)=TDS    
98   220 CONTINUE  
99     
100 C...Iterate direction of axis until stable maximum. 
101       P(N+NP+ILD,4)=0.  
102       ILG=0 
103   230 ILG=ILG+1 
104       THP=0.    
105   240 THPS=THP  
106       DO 250 J=1,3  
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)    
109   250 TPR(J)=0. 
110       DO 260 I=N+1,N+NP 
111       SGN=SIGN(P(I,5),TDI(1)*P(I,1)+TDI(2)*P(I,2)+TDI(3)*P(I,3))    
112       DO 260 J=1,4-ILD  
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 
116     
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    
120         IAGR=0  
121         SGN=(-1.)**INT(RLU_HIJING(0)+0.5)  
122         DO 270 J=1,3    
123   270   P(N+NP+ILD,J)=SGN*TPR(J)/(PS*THP)   
124         P(N+NP+ILD,4)=THP   
125         P(N+NP+ILD,5)=0.    
126       ENDIF 
127       IAGR=IAGR+1   
128   280 IF(IAGR.LT.MSTU(45).AND.ILG.LT.MIN(10,NC)) GOTO 230   
129     
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)   
134       P(N+NP+3,3)=0.    
135       THP=0.    
136       DO 290 I=N+1,N+NP 
137   290 THP=THP+P(I,5)*ABS(P(N+NP+3,1)*P(I,1)+P(N+NP+3,2)*P(I,2)) 
138       P(N+NP+3,4)=THP/PS    
139       P(N+NP+3,5)=0.    
140     
141 C...Fill axis information. Rotate back to original coordinate system.   
142       DO 300 ILD=1,3    
143       K(N+ILD,1)=31 
144       K(N+ILD,2)=96 
145       K(N+ILD,3)=ILD    
146       K(N+ILD,4)=0  
147       K(N+ILD,5)=0  
148       DO 300 J=1,5  
149       P(N+ILD,J)=P(N+NP+ILD,J)  
150   300 V(N+ILD,J)=0. 
151       CALL LUDBRB_HIJING(N+1,N+3,THE,PHI,0D0,0D0,0D0)  
152     
153 C...Select storing option. Calculate thurst and oblateness. 
154       MSTU(61)=N+1  
155       MSTU(62)=NP   
156       IF(MSTU(43).LE.1) MSTU(3)=3   
157       IF(MSTU(43).GE.2) N=N+3   
158       THR=P(N+1,4)  
159       OBL=P(N+2,4)-P(N+3,4) 
160     
161       RETURN    
162       END