[u/mrichter/AliRoot.git] / GEANT321 / ghits / giplan.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:10  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.20  by  S.Giani
*-- Author :
      SUBROUTINE GIPLAN(YC,X1,X2,S1,S2,IC,XINT,SINT,PZINT,IFLAG)
C.
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       Calculates intersection of track (X1,X2)                 *
C.    *       with plane parallel to (X-Z)                             *
C.    *        The track is approximated by a cubic in the             *
C.    *       track length.                                            *
C.    *       To improve stability, the coordinate system              *
C.    *       is shifted.                                              *
C.    *       input parameters                                         *
C.    *        YC    = Y COORDINATE OF PLANE                           *
C.    *        X1    = X,Y,Z,XP,YP,ZP OF 1ST POINT                     *
C.    *        X2    =                   2ND                           *
C.    *        S1(2) = S AT 1ST(2ND) POINT                             *
C.    *        IC    = 1 STRAIGHT LINE DEFINED BY X+XP                 *
C.    *        IC    = 2 STRAIGHT LINE DEFINED BY X1+X2                *
C.    *        IC    = 3 CUBIC MODEL                                   *
C.    *                                                                *
C.    *      output parameters                                         *
C.    *        XINT  = X,Y,Z,XP,YP,ZP AT INTERSECTION POINT            *
C.    *        SINT  = S AT INTERSECTION POINT                         *
C.    *        PZINT = PHI,Z,DPHI/DR,DZ/DR                             *
C.    *        IFLAG = 1 IF TRACK INTERSECTS PLANE                     *
C.    *              = 0 IF NOT                                        *
C.    *                                                                *
C.    *      Warning : the default accuracy is 10 microns. The value   *
C.    *      of EPSI must be changed for a better precision            *
C.    *                                                                *
C.    *    ==>Called by : <USER>, GUDIGI                               *
C.    *                                                                *
C.    *        Authors: R.BRUN/JJ.DUMONT from an original routine by   *
C.    *       H. BOERNER  KEK  OCTOBER 1982                            *
C.    *                                                                *
C.    *                                                                *
C.    ******************************************************************
C.
      DIMENSION X1(6),X2(6),XINT(6),PZINT(4),A(4),B(4),C(4)
C
      DATA MAXHIT/100/
      DATA EPSI/0.001/
C.
C.
C.    ------------------------------------------------------------------
C.
C.
      IFLAG  = 1
      DRCTN  = 1.
C
C             Track crossing the plane from above or below ?
C
      IF (X2(2).LT.X1(2))                        GO TO 5
      IF (YC.LT.X1(2))                           GO TO 90
      IF (YC.GT.X2(2))                           GO TO 90
      IF (IC.EQ.2) GO TO 30
      IF (IC.EQ.3) GO TO 7
C
      S=S1
      DXDS=X1(4)
      DYDS=X1(5)
      DZDS=X1(6)
      BX=X1(1)-DXDS*(X1(2)-YC)/DYDS
      BZ=X1(3)-DZDS*(X1(2)-YC)/DYDS
      TRL2=(BX-X1(1))**2+(X1(2)-YC)**2+(BZ-X1(3))**2
      GO TO 40
C
   5  IF (YC.LT.X2(2))                           GO TO 90
      IF (YC.GT.X1(2))                           GO TO 90
      IF(IC.EQ.2) GO TO 30
      DRCTN  = - 1.
C
      IF(IC.EQ.3) GOTO 7
      S=S2
      DXDS=X2(4)
      DYDS=X2(5)
      DZDS=X2(6)
      BX=X2(1)-DXDS*(X2(2)-YC)/DYDS
      BZ=X2(3)-DZDS*(X2(2)-YC)/DYDS
      TRL2=(BX-X2(1))**2+(X2(2)-YC)**2+(BZ-X2(3))**2
      GOTO 40
C
   30 DX=X2(1)-X1(1)
      DY=X2(2)-X1(2)
      DZ=X2(3)-X1(3)
      DS=SQRT(DX*DX+DY*DY+DZ*DZ)
      S=S1
      DXDS=DX/DS
      DYDS=DY/DS
      DZDS=DZ/DS
      BX=X1(1)-DX*(X1(2)-YC)/DY
      BZ=X1(3)-DZ*(X1(2)-YC)/DY
      TRL2=(BX-X1(1))**2+(X1(2)-YC)**2+(BZ-X1(3))**2
C
   40 TRLEN=SQRT(TRL2)*DRCTN+S
      XINT(1)=BX
      XINT(2)=YC
      XINT(3)=BZ
      SINT=TRLEN
      XINT(4)=DXDS
      XINT(5)=DYDS
      XINT(6)=DZDS
      GO TO 200
C
C               Shift coordinate system such that center of gravity=0
C
   7  IF(YC.LE.0.)                                GO TO 90
      SHIFTX = (X1(1) + X2(1)) * 0.5
      SHIFTY = (X1(2) + X2(2)) * 0.5
      SHIFTZ = (X1(3) + X2(3)) * 0.5
      SHIFTS = (S1 + S2) * 0.5
C
C             Only one value necessary since X1= -X2 etc...
C
      XSHFT  = X1(1) - SHIFTX
      YSHFT  = X1(2) - SHIFTY
      ZSHFT  = X1(3) - SHIFTZ
      SSHFT  = S1 - SHIFTS
C
      PABS1  = SQRT(X1(4)**2 + X1(5)**2 + X1(6)**2)
      PABS2  = SQRT(X2(4)**2 + X2(5)**2 + X2(6)**2)
      IF (PABS1.EQ.0..OR.PABS2.EQ.0.)            GO TO 90
C
C              Parametrize the track by a cubic through X1, X2
C
      CALL GCUBS(SSHFT,XSHFT,X1(4)/PABS1,X2(4)/PABS2,A)
      CALL GCUBS(SSHFT,YSHFT,X1(5)/PABS1,X2(5)/PABS2,B)
      CALL GCUBS(SSHFT,ZSHFT,X1(6)/PABS1,X2(6)/PABS2,C)
C
C              Iterate to find the track length corresponding to
C              the intersection of track and plane.
C              Start at S=0. middle of the shifted interval.
C
      DINTER = ABS(S2 - S1) * 0.5
      S      = 0.
C
      DO 10 I = 1,MAXHIT
         Y = SHIFTY + B(1) + S * (B(2) + S * (B(3) + S * B(4)))
         DR=(YC-Y)*DRCTN
         IF (ABS(DR).LT.EPSI)                     GO TO 20
         DINTER = DINTER * 0.5
         IF (DR.LT.0.)S = S - DINTER
         IF (DR.GE.0.)S = S + DINTER
  10  CONTINUE
C
C              Compute intersection in original coordinates
C
  20  CONTINUE
      XINT(1) = SHIFTX + A(1) + S * (A(2) + S * (A(3) + S * A(4)))
      XINT(2)=YC
      XINT(3) = SHIFTZ + C(1) + S * (C(2) + S * (C(3) + S * C(4)))
      XINT(4) = A(2) + S * (2. * A(3) + 3. * S * A(4))
      XINT(5) = B(2) + S * (2. * B(3) + 3. * S * B(4))
      XINT(6) = C(2) + S * (2. * C(3) + 3. * S * C(4))
C
C             Compute PHIHIT,ZHIT and corresponding derivatives
C
      SINT   = S + SHIFTS
  200 TERM   = 1. / (XINT(4) * XINT(1) + XINT(5) * XINT(2))
      PZINT(1) = ATAN2(XINT(2),XINT(1))
      PZINT(2) = XINT(3)
      PZINT(3) = (XINT(1) * XINT(5) - XINT(2) * XINT(4)) * TERM / YC
      PZINT(4) = TERM * XINT(6) * YC
      RETURN
C
  90  IFLAG  = 0
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:21:10 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.20 by S.Giani
	11	*-- Author :
	12	SUBROUTINE GIPLAN(YC,X1,X2,S1,S2,IC,XINT,SINT,PZINT,IFLAG)
	13	C.
	14	C.
	15	C. ******************************************************************
	16	C. * *
	17	C. * Calculates intersection of track (X1,X2) *
	18	C. * with plane parallel to (X-Z) *
	19	C. * The track is approximated by a cubic in the *
	20	C. * track length. *
	21	C. * To improve stability, the coordinate system *
	22	C. * is shifted. *
	23	C. * input parameters *
	24	C. * YC = Y COORDINATE OF PLANE *
	25	C. * X1 = X,Y,Z,XP,YP,ZP OF 1ST POINT *
	26	C. * X2 = 2ND *
	27	C. * S1(2) = S AT 1ST(2ND) POINT *
	28	C. * IC = 1 STRAIGHT LINE DEFINED BY X+XP *
	29	C. * IC = 2 STRAIGHT LINE DEFINED BY X1+X2 *
	30	C. * IC = 3 CUBIC MODEL *
	31	C. * *
	32	C. * output parameters *
	33	C. * XINT = X,Y,Z,XP,YP,ZP AT INTERSECTION POINT *
	34	C. * SINT = S AT INTERSECTION POINT *
	35	C. * PZINT = PHI,Z,DPHI/DR,DZ/DR *
	36	C. * IFLAG = 1 IF TRACK INTERSECTS PLANE *
	37	C. * = 0 IF NOT *
	38	C. * *
	39	C. * Warning : the default accuracy is 10 microns. The value *
	40	C. * of EPSI must be changed for a better precision *
	41	C. * *
	42	C. * ==>Called by : <USER>, GUDIGI *
	43	C. * *
	44	C. * Authors: R.BRUN/JJ.DUMONT from an original routine by *
	45	C. * H. BOERNER KEK OCTOBER 1982 *
	46	C. * *
	47	C. * *
	48	C. ******************************************************************
	49	C.
	50	DIMENSION X1(6),X2(6),XINT(6),PZINT(4),A(4),B(4),C(4)
	51	C
	52	DATA MAXHIT/100/
	53	DATA EPSI/0.001/
	54	C.
	55	C.
	56	C. ------------------------------------------------------------------
	57	C.
	58	C.
	59	IFLAG = 1
	60	DRCTN = 1.
	61	C
	62	C Track crossing the plane from above or below ?
	63	C
	64	IF (X2(2).LT.X1(2)) GO TO 5
65	IF (YC.LT.X1(2)) GO TO 90
66	IF (YC.GT.X2(2)) GO TO 90
67	IF (IC.EQ.2) GO TO 30
68	IF (IC.EQ.3) GO TO 7
69	C
70	S=S1
71	DXDS=X1(4)
72	DYDS=X1(5)
73	DZDS=X1(6)
74	BX=X1(1)-DXDS*(X1(2)-YC)/DYDS
75	BZ=X1(3)-DZDS*(X1(2)-YC)/DYDS
76	TRL2=(BX-X1(1))2+(X1(2)-YC)2+(BZ-X1(3))**2
77	GO TO 40
78	C
79	5 IF (YC.LT.X2(2)) GO TO 90
80	IF (YC.GT.X1(2)) GO TO 90
81	IF(IC.EQ.2) GO TO 30
82	DRCTN = - 1.
83	C
84	IF(IC.EQ.3) GOTO 7
85	S=S2
86	DXDS=X2(4)
87	DYDS=X2(5)
88	DZDS=X2(6)
89	BX=X2(1)-DXDS*(X2(2)-YC)/DYDS
90	BZ=X2(3)-DZDS*(X2(2)-YC)/DYDS
91	TRL2=(BX-X2(1))2+(X2(2)-YC)2+(BZ-X2(3))**2
92	GOTO 40
93	C
94	30 DX=X2(1)-X1(1)
95	DY=X2(2)-X1(2)
96	DZ=X2(3)-X1(3)
97	DS=SQRT(DXDX+DYDY+DZ*DZ)
98	S=S1
99	DXDS=DX/DS
100	DYDS=DY/DS
101	DZDS=DZ/DS
102	BX=X1(1)-DX*(X1(2)-YC)/DY
103	BZ=X1(3)-DZ*(X1(2)-YC)/DY
104	TRL2=(BX-X1(1))2+(X1(2)-YC)2+(BZ-X1(3))**2
105	C
106	40 TRLEN=SQRT(TRL2)*DRCTN+S
107	XINT(1)=BX
108	XINT(2)=YC
109	XINT(3)=BZ
110	SINT=TRLEN
111	XINT(4)=DXDS
112	XINT(5)=DYDS
113	XINT(6)=DZDS
114	GO TO 200
115	C
116	C Shift coordinate system such that center of gravity=0
117	C
118	7 IF(YC.LE.0.) GO TO 90
119	SHIFTX = (X1(1) + X2(1)) * 0.5
120	SHIFTY = (X1(2) + X2(2)) * 0.5
121	SHIFTZ = (X1(3) + X2(3)) * 0.5
122	SHIFTS = (S1 + S2) * 0.5
123	C
124	C Only one value necessary since X1= -X2 etc...
125	C
126	XSHFT = X1(1) - SHIFTX
127	YSHFT = X1(2) - SHIFTY
128	ZSHFT = X1(3) - SHIFTZ
129	SSHFT = S1 - SHIFTS
130	C
131	PABS1 = SQRT(X1(4)2 + X1(5)2 + X1(6)**2)
132	PABS2 = SQRT(X2(4)2 + X2(5)2 + X2(6)**2)
133	IF (PABS1.EQ.0..OR.PABS2.EQ.0.) GO TO 90
134	C
135	C Parametrize the track by a cubic through X1, X2
136	C
137	CALL GCUBS(SSHFT,XSHFT,X1(4)/PABS1,X2(4)/PABS2,A)
138	CALL GCUBS(SSHFT,YSHFT,X1(5)/PABS1,X2(5)/PABS2,B)
139	CALL GCUBS(SSHFT,ZSHFT,X1(6)/PABS1,X2(6)/PABS2,C)
140	C
141	C Iterate to find the track length corresponding to
142	C the intersection of track and plane.
143	C Start at S=0. middle of the shifted interval.
144	C
145	DINTER = ABS(S2 - S1) * 0.5
146	S = 0.
147	C
148	DO 10 I = 1,MAXHIT
149	Y = SHIFTY + B(1) + S * (B(2) + S * (B(3) + S * B(4)))
150	DR=(YC-Y)*DRCTN
151	IF (ABS(DR).LT.EPSI) GO TO 20
152	DINTER = DINTER * 0.5
153	IF (DR.LT.0.)S = S - DINTER
154	IF (DR.GE.0.)S = S + DINTER
155	10 CONTINUE
156	C
157	C Compute intersection in original coordinates
158	C
159	20 CONTINUE
160	XINT(1) = SHIFTX + A(1) + S * (A(2) + S * (A(3) + S * A(4)))
161	XINT(2)=YC
162	XINT(3) = SHIFTZ + C(1) + S * (C(2) + S * (C(3) + S * C(4)))
163	XINT(4) = A(2) + S * (2. * A(3) + 3. * S * A(4))
164	XINT(5) = B(2) + S * (2. * B(3) + 3. * S * B(4))
165	XINT(6) = C(2) + S * (2. * C(3) + 3. * S * C(4))
166	C
167	C Compute PHIHIT,ZHIT and corresponding derivatives
168	C
169	SINT = S + SHIFTS
170	200 TERM = 1. / (XINT(4) * XINT(1) + XINT(5) * XINT(2))
171	PZINT(1) = ATAN2(XINT(2),XINT(1))
172	PZINT(2) = XINT(3)
173	PZINT(3) = (XINT(1) * XINT(5) - XINT(2) * XINT(4)) * TERM / YC
174	PZINT(4) = TERM * XINT(6) * YC
175	RETURN
176	C
177	90 IFLAG = 0
178	END