]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/erpremc/trprfn.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
the MIXT geometry (IHEP+GPS2) has been introduced
[u/mrichter/AliRoot.git] / GEANT321 / erpremc / trprfn.F
CommitLineData
fe4da5cc 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#include "geant321/pilot.h"
10*CMZ : 3.21/02 29/03/94 15.41.49 by S.Giani
11*-- Author :
12 SUBROUTINE TRPRFN(X1,P1,H1,X2,P2,H2,CH,XL,R,MVAR,IFLAG,ITRAN,IERR)
13************************************************************************
14*
15*
16* SUBR. TRPRFN(X1,P1,H1,X2,P2,H2,CH,XL,*R*,MVAR,IFLAG,ITRAN,IERR*)
17*
18* Origin W.Wittek EMCSW/81/18
19*
20* Finite step length case coded by V.Innocente ( Feb. 88 )
21*
22* code improved: V.Innocente ( April. 90 )
23* inline code replaces external function
24* code improved: V.Innocente ( January 91 )
25* effect of energy loss added
26*
27*_______________________________________________________________________
28*
29* *** ERROR PROPAGATION ALONG A PARTICLE TRAJECTORY IN A MAGNETIC FIELD
30* ROUTINE ASSUMES THAT IN THE INTERVAL (X1,X2) THE QUANTITIES 1/P
31* AND (HX,HY,HZ) ARE CONSTANT.
32*
33* *** IFLAG = -1 INITIALIZATION, TRANSFORMATION OF ERROR MATRIX FROM
34* EXTERNAL TO SC VARIABLES
35* = 0 ERROR PROPAGATION FROM X1 TO X2
36* = 1 TRANSFORMATION OF ERROR MATRIX FROM SC TO
37* EXTERNAL VARIABLES
38*
39* ITRAN USED FOR IFLAG = 0 OR 1 ONLY
40* = 0 TRANSFORMATION MATRIX IS UPDATED ,BUT ERROR MATRIX IS NOT
41* TRANSFORMED
42* = 1 TRANSF. MATRIX IS UPDATED AND ERROR MATRIX IS TRANSFORMED
43*
44* MVAR SPECIFIES TYPE OF EXTERNAL VARIABLES
45* = 0 ( 1/P,LAMBDA,PHI,YT, ZT ; SC )
46* = 1 ( 1/P, Y', Z', Y, Z ; SPLINE )
47*
48* *** X1, P1, H1 X,Y,Z COMPONENTS OF POSITION, MOMENTUM AND MAGNETIC INPUT
49* FIELD VECTOR/GRADIENT AT STARTING POINT OF INTERVAL
50* X2, P2, H2 ...... AT END POINT OF INTERVAL INPUT
51* CH CHARGE OF PARTICLE INPUT
52* XL PATHLENGTH FROM X1 TO X2 ( NEGATIVE IF OPPOSITE
53* TO ACTUAL MOVEMENT OF PARTICLE ) INPUT
54* R ERROR MATRIX (TRIANGLE) INPUT/OUTPUT
55* B 5 * 5 TRANSFORMATION MATRIX FOR ERRORS IN
56* SC VARIABLES OUTPUT
57*
58* *** IERR = 1 ILLEGAL VALUE OF MVAR OUTPUT
59* 2 MOMENTUM IS ZERO
60* 3 H*ALFA/P AT X1 AND X2 DIFFER TOO MUCH
61* 4 PARTICLE MOVES IN Z - DIRECTION
62*
63************************************************************************
64*
65#if !defined(CERNLIB_SINGLE)
66 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
67 REAL X1,P1,H1,X2,P2,H2,R,CH,PS,PC,XL,SPX
68#endif
69#include "geant321/trcom3.inc"
70#include "geant321/gcunit.inc"
71 DIMENSION X1(3),P1(3),H1(9),X2(3),P2(3),H2(9)
72 DIMENSION R(15),PS(3),PC(3)
73*
74 DIMENSION T1(3),T2(3),U1(3),U2(3),V1(3),V2(3),HN(9)
75 DIMENSION AN1(3),AN2(3),DX(3)
76 DIMENSION HV1(3),HU1(3)
77*
78 SAVE INIT,DELHP6,CFACT8
79*
80 DATA INIT/0/
81#if !defined(CERNLIB_SINGLE)
82 DATA DELHP6/300.D0/
83*
84 DATA CFACT8 / 2.997925 D-4 /
85#endif
86#if defined(CERNLIB_SINGLE)
87 DATA DELHP6/300./
88*
89 DATA CFACT8 / 2.997925 E-4 /
90#endif
91*
92*____________________________________________________________________
93*
94 IERR=0
95 IF(IFLAG) 10, 20, 80
96*
97* *** TRANSFORM ERROR MATRIX FROM EXTERNAL TO INTERNAL VARIABLES;
98*
99 10 NEW=1
100 IF(MVAR.NE.1) GO TO 11
101 PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
102 IF(PA1.EQ.0.) GO TO 902
103 PS(1)=1./PA1
104 IF(P1(1).EQ.0.) GO TO 904
105 PS(2)=P1(2)/P1(1)
106 PS(3)=P1(3)/P1(1)
107 SPX=1.
108 IF(P1(1).LT.0.) SPX=-1.
109 CALL TRSPSC(PS,R,PC,R,H1,CH,IERR,SPX)
110 GO TO 19
111*
112 11 IF(MVAR.NE.0) GO TO 901
113 19 GO TO 900
114*
115* *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
116
117*
118 20 PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
119 PA2=SQRT(P2(1)**2+P2(2)**2+P2(3)**2)
120 IF(PA1*PA2.EQ.0.) GO TO 902
121C DPA = PA2 - PA1
122 PM1=1./PA1
123 PM2=1./PA2
124 DPM = PM2 - PM1
125*
126 DO 201 I=1,3
127 T1(I) = P1(I)*PM1
128 T2(I) = P2(I)*PM2
129201 CONTINUE
130*
131 SINL=T2(3)
132 SINL0=T1(3)
133*
134 COSL=SQRT(ABS(1.-SINL**2))
135 IF(COSL.EQ.0.) GO TO 904
136 COSL1=1./COSL
137 COSL0=SQRT(ABS(1.-SINL0**2))
138*
139* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
140* *** NEUTRAL PARTICLE OR FIELDFREE REGION
141*
142 DO 26 I=1,5
143 DO 15 K=1,5
144 A(I,K)=0.
145 15 CONTINUE
146 A(I,I)=1.
147 26 CONTINUE
148 A(4,3)=XL*COSL
149 A(5,2)=XL
150*
151 IF(CH.EQ.0.) GO TO 45
152 HA1=SQRT(H1(1)**2+H1(2)**2+H1(3)**2)
153 HA2=SQRT(H2(1)**2+H2(2)**2+H2(3)**2)
154 HAM1=HA1*PM1
155 HAM2=HA2*PM2
156 HAMX=MAX(HAM1,HAM2)
157 IF(HAMX.EQ.0.) GO TO 45
158*
159*
160*
161* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
162*
163*
164 IF(HA2.NE.0.) THEN
165 GAM=(H2(1)*T2(1)+H2(2)*T2(2)+H2(3)*T2(3))/HA2
166 ELSE
167 GAM=(H1(1)*T1(1)+H1(2)*T1(2)+H1(3)*T1(3))/HA1
168 ENDIF
169*
170 ALFA2=1.-GAM**2
171*
172 DH2=(H1(1)*PM1-H2(1)*PM2)**2+
173 1 (H1(2)*PM1-H2(2)*PM2)**2+
174 1 (H1(3)*PM1-H2(3)*PM2)**2
175 IF(DH2*ALFA2.GT.DELHP6**2) GO TO 903
176*
177* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
178*
179 PM12=(PM1+PM2)*0.5
180 P12=1./(2.*PM12)
181 HN(1)=(H1(1)*PM1+H2(1)*PM2)*P12*CH*CFACT8
182 HN(2)=(H1(2)*PM1+H2(2)*PM2)*P12*CH*CFACT8
183 HN(3)=(H1(3)*PM1+H2(3)*PM2)*P12*CH*CFACT8
184CC HN(4)=(H1(4)*PM1+H2(4)*PM2)*P12*CH*CFACT8
185CC HN(5)=(H1(5)*PM1+H2(5)*PM2)*P12*CH*CFACT8
186CC HN(6)=(H1(6)*PM1+H2(6)*PM2)*P12*CH*CFACT8
187CC HN(7)=(H1(7)*PM1+H2(7)*PM2)*P12*CH*CFACT8
188CC HN(8)=(H1(8)*PM1+H2(8)*PM2)*P12*CH*CFACT8
189CC HN(9)=(H1(9)*PM1+H2(9)*PM2)*P12*CH*CFACT8
190*
191 HM = SQRT(HN(1)**2+HN(2)**2+HN(3)**2)
192 OVER = 1./HM
193 HN(1) = OVER*HN(1)
194 HN(2) = OVER*HN(2)
195 HN(3) = OVER*HN(3)
196 PAV = .5*(PA1+PA2)
197 Q = - HM/PAV
198 THETA = Q*XL
199 SINT = SIN(THETA)
200 COST = COS(THETA)
201 GAMMA=HN(1)*T2(1)+HN(2)*T2(2)+HN(3)*T2(3)
202 AN2(1) = HN(2)*T2(3)-HN(3)*T2(2)
203 AN2(2) = HN(3)*T2(1)-HN(1)*T2(3)
204 AN2(3) = HN(1)*T2(2)-HN(2)*T2(1)
205*
206 AU = 1./SQRT(T1(1)**2+T1(2)**2)
207 U1(1) = -AU*T1(2)
208 U1(2) = AU*T1(1)
209 U1(3) = 0.D0
210 V1(1) = -T1(3)*U1(2)
211 V1(2) = T1(3)*U1(1)
212 V1(3) = T1(1)*U1(2)-T1(2)*U1(1)
213*
214 AU = 1./SQRT(T2(1)**2+T2(2)**2)
215 U2(1) = -AU*T2(2)
216 U2(2) = AU*T2(1)
217 U2(3) = 0.D0
218 V2(1) = -T2(3)*U2(2)
219 V2(2) = T2(3)*U2(1)
220 V2(3) = T2(1)*U2(2)-T2(2)*U2(1)
221*
222 DX(1) = X1(1) - X2(1)
223 DX(2) = X1(2) - X2(2)
224 DX(3) = X1(3) - X2(3)
225*
226*
227* *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
228* *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
229* *** TAKEN INTO ACCOUNT
230*
231 30 CONTINUE
232 QP = Q*PAV
233 ANV = -(HN(1)*U2(1)+HN(2)*U2(2) )
234 ANU = (HN(1)*V2(1)+HN(2)*V2(2)+HN(3)*V2(3))
235 OMCOST = 1.-COST
236 TMSINT = THETA-SINT
237*
238 HU1(1) = -HN(3)*U1(2)
239 HU1(2) = HN(3)*U1(1)
240 HU1(3) = HN(1)*U1(2)-HN(2)*U1(1)
241*
242 HV1(1) = HN(2)*V1(3)-HN(3)*V1(2)
243 HV1(2) = HN(3)*V1(1)-HN(1)*V1(3)
244 HV1(3) = HN(1)*V1(2)-HN(2)*V1(1)
245*
246*** 1/P
247*
248 A(1,1) = 1.-DPM*PAV*(1.+(T2(1)*DX(1)+T2(2)*DX(2)+T2(3)*DX(3))/XL)
249 + +2.*DPM*PAV
250*
251 A(1,2) = -DPM/THETA*
252 1 ( TMSINT*GAMMA*(HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3)) +
253 2 SINT*(V1(1)*T2(1)+V1(2)*T2(2)+V1(3)*T2(3)) +
254 3 OMCOST*(HV1(1)*T2(1)+HV1(2)*T2(2)+HV1(3)*T2(3)) )
255*
256 A(1,3) = -COSL0*DPM/THETA*
257 1 ( TMSINT*GAMMA*(HN(1)*U1(1)+HN(2)*U1(2) ) +
258 2 SINT*(U1(1)*T2(1)+U1(2)*T2(2) ) +
259 3 OMCOST*(HU1(1)*T2(1)+HU1(2)*T2(2)+HU1(3)*T2(3)) )
260*
261 A(1,4) = -DPM/XL*(U1(1)*T2(1)+U1(2)*T2(2) )
262*
263 A(1,5) = -DPM/XL*(V1(1)*T2(1)+V1(2)*T2(2)+V1(3)*T2(3))
264*
265*** Lambda
266*
267 A(2,1) = -QP*ANV*(T2(1)*DX(1)+T2(2)*DX(2)+T2(3)*DX(3))
268 + *(1.+DPM*PAV)
269*
270 A(2,2) = COST*(V1(1)*V2(1)+V1(2)*V2(2)+V1(3)*V2(3)) +
271 + SINT*(HV1(1)*V2(1)+HV1(2)*V2(2)+HV1(3)*V2(3)) +
272 1 OMCOST*(HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3))*
273 A (HN(1)*V2(1)+HN(2)*V2(2)+HN(3)*V2(3)) +
274 2 ANV*( -SINT*(V1(1)*T2(1)+V1(2)*T2(2)+V1(3)*T2(3)) +
275 3 OMCOST*(V1(1)*AN2(1)+V1(2)*AN2(2)+V1(3)*AN2(3)) -
276 4 TMSINT*GAMMA*(HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3)) )
277*
278 A(2,3) = COST*(U1(1)*V2(1)+U1(2)*V2(2) ) +
279 + SINT*(HU1(1)*V2(1)+HU1(2)*V2(2)+HU1(3)*V2(3)) +
280 1 OMCOST*(HN(1)*U1(1)+HN(2)*U1(2) )*
281 A (HN(1)*V2(1)+HN(2)*V2(2)+HN(3)*V2(3)) +
282 2 ANV*( -SINT*(U1(1)*T2(1)+U1(2)*T2(2) ) +
283 3 OMCOST*(U1(1)*AN2(1)+U1(2)*AN2(2) ) -
284 4 TMSINT*GAMMA*(HN(1)*U1(1)+HN(2)*U1(2) ) )
285 A(2,3) = COSL0*A(2,3)
286*
287 A(2,4) = -Q*ANV*(U1(1)*T2(1)+U1(2)*T2(2) )
288*
289 A(2,5) = -Q*ANV*(V1(1)*T2(1)+V1(2)*T2(2)+V1(3)*T2(3))
290*
291*** Phi
292*
293 A(3,1) = -QP*ANU*(T2(1)*DX(1)+T2(2)*DX(2)+T2(3)*DX(3))*COSL1
294 + *(1.+DPM*PAV)
295*
296 A(3,2) = COST*(V1(1)*U2(1)+V1(2)*U2(2) ) +
297 + SINT*(HV1(1)*U2(1)+HV1(2)*U2(2) ) +
298 1 OMCOST*(HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3))*
299 A (HN(1)*U2(1)+HN(2)*U2(2) ) +
300 2 ANU*( -SINT*(V1(1)*T2(1)+V1(2)*T2(2)+V1(3)*T2(3)) +
301 3 OMCOST*(V1(1)*AN2(1)+V1(2)*AN2(2)+V1(3)*AN2(3)) -
302 4 TMSINT*GAMMA*(HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3)) )
303 A(3,2) = COSL1*A(3,2)
304*
305 A(3,3) = COST*(U1(1)*U2(1)+U1(2)*U2(2) ) +
306 + SINT*(HU1(1)*U2(1)+HU1(2)*U2(2) ) +
307 1 OMCOST*(HN(1)*U1(1)+HN(2)*U1(2) )*
308 A (HN(1)*U2(1)+HN(2)*U2(2) ) +
309 2 ANU*( -SINT*(U1(1)*T2(1)+U1(2)*T2(2) ) +
310 3 OMCOST*(U1(1)*AN2(1)+U1(2)*AN2(2) ) -
311 4 TMSINT*GAMMA*(HN(1)*U1(1)+HN(2)*U1(2) ) )
312 A(3,3) = COSL1*COSL0*A(3,3)
313*
314 A(3,4) = -Q*ANU*(U1(1)*T2(1)+U1(2)*T2(2) )*COSL1
315*
316 A(3,5) = -Q*ANU*(V1(1)*T2(1)+V1(2)*T2(2)+V1(3)*T2(3))*COSL1
317*
318*** Yt
319*
320 A(4,1) = PAV*(U2(1)*DX(1)+U2(2)*DX(2) )
321 + *(1.+DPM*PAV)
322*
323 A(4,2) = ( SINT*(V1(1)*U2(1)+V1(2)*U2(2) ) +
324 1 OMCOST*(HV1(1)*U2(1)+HV1(2)*U2(2) ) +
325 2 TMSINT*(HN(1)*U2(1)+HN(2)*U2(2) )*
326 3 (HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3)) )/Q
327*
328 A(4,3) = ( SINT*(U1(1)*U2(1)+U1(2)*U2(2) ) +
329 1 OMCOST*(HU1(1)*U2(1)+HU1(2)*U2(2) ) +
330 2 TMSINT*(HN(1)*U2(1)+HN(2)*U2(2) )*
331 3 (HN(1)*U1(1)+HN(2)*U1(2) ) )*COSL0/Q
332*
333 A(4,4) = (U1(1)*U2(1)+U1(2)*U2(2) )
334*
335 A(4,5) = (V1(1)*U2(1)+V1(2)*U2(2) )
336*
337*** Zt
338*
339 A(5,1) = PAV*(V2(1)*DX(1)+V2(2)*DX(2)+V2(3)*DX(3))
340 + *(1.+DPM*PAV)
341*
342 A(5,2) = ( SINT*(V1(1)*V2(1)+V1(2)*V2(2)+V1(3)*V2(3)) +
343 1 OMCOST*(HV1(1)*V2(1)+HV1(2)*V2(2)+HV1(3)*V2(3)) +
344 2 TMSINT*(HN(1)*V2(1)+HN(2)*V2(2)+HN(3)*V2(3))*
345 3 (HN(1)*V1(1)+HN(2)*V1(2)+HN(3)*V1(3)) )/Q
346*
347 A(5,3) = ( SINT*(U1(1)*V2(1)+U1(2)*V2(2) ) +
348 1 OMCOST*(HU1(1)*V2(1)+HU1(2)*V2(2)+HU1(3)*V2(3)) +
349 2 TMSINT*(HN(1)*V2(1)+HN(2)*V2(2)+HN(3)*V2(3))*
350 3 (HN(1)*U1(1)+HN(2)*U1(2) ) )*COSL0/Q
351*
352 A(5,4) = (U1(1)*V2(1)+U1(2)*V2(2) )
353*
354 A(5,5) = (V1(1)*V2(1)+V1(2)*V2(2)+V1(3)*V2(3))
355 45 CONTINUE
356*
357* *** NEW = 0 TRANSFORMATION MATRIX IS UPDATED
358* 1 TRANSFORMATION MATRIX IS INITIALIZED
359*
360 IF(NEW.EQ.0) GO TO 23
361 NEW=0
362 DO 25 I=1,5
363 DO 24 K=1,5
364 B(I,K)=A(I,K)
365 24 CONTINUE
366 25 CONTINUE
367 GO TO 27
368 23 CONTINUE
369*
370 CALL XMM55(A,B,B)
371*
372 27 CONTINUE
373 80 IF(ITRAN.EQ.0) GO TO 90
374*
375*
376 J=0
377 DO 22 I=1,5
378 DO 21 K=I,5
379 J=J+1
380 S(J)=R(J)
381 21 CONTINUE
382 22 CONTINUE
383*
384*
385* *** TRANSFORM ERROR MATRIX
386*
387 CALL SSMT5T(B,S,S)
388*
389 NEW=1
390
391 J=0
392 DO 41 I=1,5
393 DO 40 K=I,5
394 J=J+1
395 R(J)=S(J)
396 40 CONTINUE
397 41 CONTINUE
398*
399 90 IF(IFLAG.LE.0) GO TO 900
400*
401*
402* *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES;
403*
404*
405 NEW=1
406 IF(MVAR.NE.1) GO TO 91
407 PC(1)=PM2
408 PC(2)=ASIN(P2(3)*PC(1))
409 IF (ABS (P2(1)) .LT. 1.E-30) P2(1) = 1.E-30
410 PC(3)=ATAN2(P2(2),P2(1))
411 CALL TRSCSP(PC,R,PS,R,H2,CH,IERR,SPX)
412 GO TO 900
413*
414 91 IF(MVAR.NE.0) GO TO 901
415 GO TO 900
416*
417* *** ERROR EXITS
418*
419 901 IERR=1
420 GO TO 999
421 902 IERR=2
422 GO TO 999
423 903 IERR=3
424C IF(INIT.NE.0) GO TO 30
425* WRITE (LOUT, 998) DH2,ALFA2,XL
426 998 FORMAT('0',' *** S/R TRPROP DELTA(H*ALFA/P)',5X
427 1,'EXCEEDS TOLERANCE '/'0',3E12.5//' ********** ',///)
428 INIT=1
429 GO TO 30
430 904 IERR=4
431 999 WRITE (LOUT, 1000) IERR
432 1000 FORMAT(1H ,' *** S/R ERPROP IERR =',I5)
433*
434 900 CONTINUE
435 END
436