This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / erpremc / trprop.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/03/06 15:37:35  mclareni
6 * Add geane321 source directories
7 *
8 *
9 *      MODIFIED FOR MAGFIELD GRADIENT 02-APR-82
10 #include "geant321/pilot.h"
11 *CMZ :  3.21/02 29/03/94  15.41.49  by  S.Giani
12 *-- Author :
13 C
14       SUBROUTINE TRPROP(X1,P1,H1,X2,P2,H2,CH,XL,R,MVAR,IFLAG,ITRAN,IERR)
15 C
16 C *** ERROR PROPAGATION ALONG A PARTICLE TRAJECTORY IN A MAGNETIC FIELD
17 C     ROUTINE ASSUMES THAT IN THE INTERVAL (X1,X2) THE QUANTITIES 1/P
18 C     AND (HX,HY,HZ) ARE RATHER CONSTANT. DELTA(PHI) MUST NOT BE TOO LARGE
19 C
20 C     Authors: A. Haas and W. Wittek
21 C
22 C *** IFLAG  =  -1   INITIALIZATION, TRANSFORMATION OF ERROR MATRIX FROM
23 C                    EXTERNAL TO SC VARIABLES
24 C            =   0   ERROR PROPAGATION FROM X1 TO X2
25 C            =   1   TRANSFORMATION OF ERROR MATRIX FROM SC TO
26 C                    EXTERNAL VARIABLES
27 C
28 C     ITRAN          USED FOR IFLAG = 0 OR 1 ONLY
29 C            =   0   TRANSFORMATION MATRIX IS UPDATED ,BUT ERROR MATRIX IS NOT
30 C                    TRANSFORMED
31 C           =    1   TRANSF. MATRIX IS UPDATED  AND ERROR MATRIX IS TRANSFORMED
32 C
33 C     MVAR           SPECIFIES TYPE OF EXTERNAL VARIABLES
34 C            =   0   ( 1/P,LAMBDA,PHI,YT, ZT ;   SC   )
35 C            =   1   ( 1/P,  Y',  Z',  Y,  Z ; SPLINE )
36 C
37 C *** X1, P1, H1     X,Y,Z COMPONENTS OF POSITION, MOMENTUM AND MAGNETIC   INPUT
38 C                    FIELD VECTOR/GRADIENT AT STARTING POINT OF INTERVAL
39 C     X2, P2, H2     ......  AT END POINT OF INTERVAL                      INPUT
40 C     CH             CHARGE OF PARTICLE                                    INPUT
41 C     XL             PATHLENGTH FROM X1 TO X2   ( NEGATIVE IF OPPOSITE
42 C                    TO ACTUAL MOVEMENT OF PARTICLE )                      INPUT
43 C     R              ERROR MATRIX  (TRIANGLE)                       INPUT/OUTPUT
44 C     B              5 * 5 TRANSFORMATION MATRIX FOR ERRORS IN
45 C                    SC VARIABLES                                         OUTPUT
46 C
47 C *** IERR   =  1    ILLEGAL VALUE OF MVAR                                OUTPUT
48 C               2    MOMENTUM IS ZERO
49 C               3    H*ALFA/P AT X1 AND X2 DIFFER TOO MUCH
50 C                    OR DELTA PHI IS TOO LARGE
51 C               4    PARTICLE MOVES IN Z - DIRECTION
52 C
53 #if !defined(CERNLIB_SINGLE)
54       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
55       REAL  X1,P1,H1,X2,P2,H2,R,CH,PS,PC,XL,SPX
56 #endif
57 #include "geant321/trcom3.inc"
58 #include "geant321/gcunit.inc"
59       DIMENSION X1(3),P1(3),H1(9),X2(3),P2(3),H2(9)
60       DIMENSION R(15),PS(3),PC(3),HN(9)
61 C
62       SAVE INIT,DELHP6,CFACT8
63 *
64       DATA INIT/0/
65 #if defined(CERNLIB_SINGLE)
66       DATA DELHP6/300./,DELFI6/0.1/
67       DATA CFACT8 / 2.997925 E-4 /
68 #endif
69 #if !defined(CERNLIB_SINGLE)
70       DATA DELHP6/300.D0/,DELFI6/0.1D0/
71       DATA CFACT8 / 2.997925 D-4 /
72 #endif
73 C
74       IERR=0
75       IF(IFLAG) 10, 20, 80
76 C
77 C *** TRANSFORM ERROR MATRIX FROM EXTERNAL TO INTERNAL VARIABLES;
78 C
79    10 NEW=1
80       IF(MVAR.NE.1) GO TO 11
81       PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
82       IF(PA1.EQ.0.) GO TO 902
83       PS(1)=1./PA1
84       IF(P1(1).EQ.0.) GO TO 904
85       PS(2)=P1(2)/P1(1)
86       PS(3)=P1(3)/P1(1)
87       SPX=1.
88       IF(P1(1).LT.0.) SPX=-1.
89       CALL TRSPSC(PS,R,PC,R,H1,CH,IERR,SPX)
90       GO TO 19
91 C
92    11 IF(MVAR.NE.0) GO TO 901
93    19 GO TO 900
94 C
95 C *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
96  
97 C
98    20 PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
99       PA2=SQRT(P2(1)**2+P2(2)**2+P2(3)**2)
100       IF(PA1*PA2.EQ.0.) GO TO 902
101       PM1=1./PA1
102       PM2=1./PA2
103 C
104       TN(1)=P1(1)+P2(1)
105       TN(2)=P1(2)+P2(2)
106       TN(3)=P1(3)+P2(3)
107       PM12=1./SQRT(TN(1)**2+TN(2)**2+TN(3)**2)
108       TN(1)=TN(1)*PM12
109       TN(2)=TN(2)*PM12
110       TN(3)=TN(3)*PM12
111 C
112       SINL=TN(3)
113       COSL=SQRT(ABS(1.-SINL**2))
114       IF(COSL.EQ.0.) GO TO 904
115       COSL1=1./COSL
116       SINP=TN(2)*COSL1
117       COSP=TN(1)*COSL1
118 C
119 C *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
120 C *** NEUTRAL PARTICLE OR FIELDFREE REGION
121 C
122       DO 26 I=1,5
123          DO 15 K=1,5
124             A(I,K)=0.
125    15    CONTINUE
126          A(I,I)=1.
127    26 CONTINUE
128       A(4,3)=XL*COSL
129       A(5,2)=XL
130 C
131       IF(CH.EQ.0.) GO TO 45
132       HA1=SQRT(H1(1)**2+H1(2)**2+H1(3)**2)
133       HA2=SQRT(H2(1)**2+H2(2)**2+H2(3)**2)
134       HAM1=HA1*PM1
135       HAM2=HA2*PM2
136       HAMX = MAX(HAM1,HAM2)
137       IF(HAMX.EQ.0.) GO TO 45
138 C
139 C *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
140 C
141       PM12=(PM1+PM2)*0.5
142       P12=1./(2.*PM12)
143       HN(1)=(H1(1)*PM1+H2(1)*PM2)*P12*CH*CFACT8
144       HN(2)=(H1(2)*PM1+H2(2)*PM2)*P12*CH*CFACT8
145       HN(3)=(H1(3)*PM1+H2(3)*PM2)*P12*CH*CFACT8
146       HN(4)=(H1(4)*PM1+H2(4)*PM2)*P12*CH*CFACT8
147       HN(5)=(H1(5)*PM1+H2(5)*PM2)*P12*CH*CFACT8
148       HN(6)=(H1(6)*PM1+H2(6)*PM2)*P12*CH*CFACT8
149       HN(7)=(H1(7)*PM1+H2(7)*PM2)*P12*CH*CFACT8
150       HN(8)=(H1(8)*PM1+H2(8)*PM2)*P12*CH*CFACT8
151       HN(9)=(H1(9)*PM1+H2(9)*PM2)*P12*CH*CFACT8
152 C
153       B0=HN(1)*COSP+HN(2)*SINP
154       B2=-HN(1)*SINP+HN(2)*COSP
155       B3=-B0*SINL+HN(3)*COSL
156       TGL=SINL*COSL1
157 C
158 C
159 C *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
160 C     AND WHETHER CHANGE OF TRACK DIRECTION DUE TO MAG.FIELD IS TOO LARGE
161 C
162 C
163       IF(HA2.EQ.0.) GO TO 29
164  
165       GAM=(H2(1)*TN(1)+H2(2)*TN(2)+H2(3)*TN(3))/HA2
166       GO TO 28
167    29 GAM=(H1(1)*TN(1)+H1(2)*TN(2)+H1(3)*TN(3))/HA1
168    28 CONTINUE
169       ALFA=SQRT(ABS(1.-GAM**2))
170 C
171       DH2=(H1(1)*PM1-H2(1)*PM2)**2+
172      1    (H1(2)*PM1-H2(2)*PM2)**2+
173      1    (H1(3)*PM1-H2(3)*PM2)**2
174       IF(DH2*ALFA**2.GT.DELHP6**2) GO TO 903
175       ALFAQ=-ALFA*CFACT8*(HAM1+HAM2)*0.5
176       DFI=ABS(XL*ALFAQ)
177       IF(DFI.GT.DELFI6) GO TO 903
178 C
179 C *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
180 C *** TAKING INTO ACCOUNT  FIELD GRADIENT PERPENDICULAR TO TRACK
181 C
182    30 COSP2=COSP*COSP
183       SINP2=SINP*SINP
184       COSIP=COSP*SINP
185 C
186       G22=SINP2*HN(9)+COSP2*HN(8)-2.0*COSIP*HN(7)
187       G33=SINL*SINL*(COSP2*HN(9)+SINP2*HN(8)+2.0*COSIP*HN(7))
188      ++COSL*(COSL*HN(6)-2.0*SINL*(COSP*HN(4)+SINP*HN(5)))
189       G23=SINL*(COSIP*(HN(9)-HN(8))+(SINP2-COSP2)*HN(7))
190      ++COSL*(COSP*HN(5)-SINP*HN(4))
191 C
192       A(2,1)=XL*B2
193       A(2,3)=-B0*XL*PM12
194       A(2,4)=(B2*B3*PM12+G22)*XL*PM12
195       A(2,5)=(-B2*B2*PM12+G23)*XL*PM12
196 C
197       A(3,1)=-XL*B3*COSL1
198       A(3,2)=B0*XL*PM12*COSL1**2
199       A(3,3)=1.+TGL*B2*XL*PM12
200       A(3,4)=(-B3*B3*PM12-G23)*XL*PM12*COSL1
201       A(3,5)=(B3*B2*PM12-G33)*XL*PM12*COSL1
202 C
203       A(4,5)=-B3*TGL*XL*PM12
204       A(5,4)=B3*TGL*XL*PM12
205 C
206    45 CONTINUE
207 C
208 C *** NEW = 0  TRANSFORMATION MATRIX IS UPDATED
209 C           1  TRANSFORMATION MATRIX IS INITIALIZED
210 C
211       IF(NEW.EQ.0) GO TO 23
212       NEW=0
213       DO 25 I=1,5
214          DO 24 K=1,5
215             B(I,K)=A(I,K)
216    24    CONTINUE
217    25 CONTINUE
218       GO TO 27
219    23 CONTINUE
220 C
221       CALL XMM55(A,B,B)
222 C
223    27 CONTINUE
224    80 IF(ITRAN.EQ.0) GO TO 90
225 C
226 C
227       J=0
228       DO 22 I=1,5
229          DO 21 K=I,5
230             J=J+1
231             S(J)=R(J)
232    21    CONTINUE
233    22 CONTINUE
234 C
235 C
236 C *** TRANSFORM ERROR MATRIX
237 C
238       CALL SSMT5T(B,S,S)
239 C
240       NEW=1
241  
242       J=0
243       DO 41 I=1,5
244          DO 40 K=I,5
245             J=J+1
246             R(J)=S(J)
247    40    CONTINUE
248    41 CONTINUE
249 C
250    90 IF(IFLAG.LE.0) GO TO 900
251 C
252 C
253 C *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES;
254 C
255 C
256       NEW=1
257       IF(MVAR.NE.1) GO TO 91
258       PC(1)=PM2
259       PC(2)=ASIN(P2(3)*PC(1))
260       IF (ABS (P2(1)) .LT. 1.E-30) P2(1) = 1.E-30
261       PC(3) = ATAN2 (P2(2),P2(1))
262       CALL TRSCSP(PC,R,PS,R,H2,CH,IERR,SPX)
263       GO TO 900
264 C
265    91 IF(MVAR.NE.0) GO TO 901
266       GO TO 900
267 C
268 C *** ERROR EXITS
269 C
270   901 IERR=1
271       GO TO 999
272   902 IERR=2
273       GO TO 999
274   903 IERR=3
275       IF(INIT.NE.0) GO TO 30
276 C     WRITE (LOUT, 998) DH2,DFI,ALFA,XL
277 C 998 FORMAT(1H0,48H *** S/R TRPROP   DELTA(H*ALFA/P)  OR DELTA(PHI),5X
278 C    1,22HEXCEEDS TOLERANCE     /1H0,4E12.5//16H **********    ,
279 C    251HATTENTION !   NO FURTHER WARNINGS WILL BE GIVEN    ///)
280       INIT=1
281       GO TO 30
282   904 IERR=4
283   999 WRITE (LOUT, 1000) IERR
284  1000 FORMAT(1H ,' *** S/R ERPROP   IERR =',I5)
285 C
286   900 RETURN
287       END