[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gnotub.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:53  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
*-- Author :
      SUBROUTINE GNOTUB(X,P,IACT,IFL,SNEXT,SNXT,SAFE)
C.    ******************************************************************
C.    *                                                                *
C.    *      Compute distance to intersection with boundary surface of *
C     *      volume TUBE or TUBS, from point X(1),X(2),X(3) outside    *
C     *      the volume along track with direction cosines X(4),X(5),  *
C     *      X(6)                                                      *
C.    *      P     (input)  : volume parameters                        *
C.    *      IACT  (input)  : action flag                              *
C.    *         = 0  Compute SAFE only                                 *
C.    *         = 1  Compute SAFE, compute SNXT only if SAFE.LT.SNEXT  *
C.    *         = 2  Compute both SAFE and SNXT                        *
C.    *         = 3  Compute SNXT only                                 *
C.    *      IFL   (input)  : 1 for TUBE, 2 for PHI segmented TUBE     *
C.    *      SNEXT (input)  : see IACT = 1                             *
C.    *      SNXT  (output) : distance to volume boundary along track  *
C.    *      SAFE  (output) : not larger than scalar distance to       *
C.    *                       volume boundaray                         *
C.    *      Called by : GNEXT, GNOPCO, GTNEXT                         *
C.    *                                                                *
C.    *      Authors   : Michel Maire and Rolf Nierhaus   21-JUN-1990  *
C.    *                                                                *
C.    ******************************************************************
C.    *                                                                *
C.    * 'TUBE'    is a tube. It has 3 parameters, the inside radius,   *
C.    *           the outside radius and the half length in z.         *
C.    * 'TUBS'    is a phi segment of a  tube.  It has 5 parameters,   *
C.    *           the same  3 as  'TUBE' plus  the phi  limits.  The   *
C.    *           segment  starts at  the first  limit and  includes   *
C.    *           increasing phi  value up  to the  second limit  or   *
C.    *           that plus 360 degrees.                               *
C.    *                                                                *
C.    ******************************************************************
#if !defined(CERNLIB_SINGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#endif
#include "geant321/gconsp.inc"
      REAL X(6),P(5),SNEXT,SNXT,SAFE
*
*     this part has to be moved outside the routine
      IF (IFL.EQ.2) THEN
         P4=P(4)*DEGRAD
         P5=P(5)*DEGRAD
         IF (P5.LT.P4) P5=P5+TWOPI
         C1=COS(P4)
         S1=SIN(P4)
         C2=COS(P5)
         S2=SIN(P5)
         FIO=0.5*(P5+P4)
         CFIO=COS(FIO)
         SFIO=SIN(FIO)
         DFI=0.5*(P5-P4)
         CDFI=COS(DFI)
*        SDFI=SIN(DFI)
      END IF
*
      SNXT=1.E10
      R=SQRT(X(1)**2+X(2)**2)
*
*     Compute SAFE radius
      IF (IACT.LT.3) THEN
         SAF1=P(1)-R
         SAF2=R-P(2)
         SAF3=ABS(X(3))-P(3)
         SAF4=0.
         IF (IFL.EQ.2.AND.R.GT.0.) THEN
            CPSI=(X(1)*CFIO+X(2)*SFIO)/R
            IF (CPSI.LT.CDFI) THEN
               IF ((X(2)*CFIO-X(1)*SFIO).LE.0.) THEN
                  SAF4=ABS(X(1)*S1-X(2)*C1)
               ELSE
                  SAF4=ABS(X(1)*S2-X(2)*C2)
               END IF
            END IF
         ENDIF
         SAFE=MAX(SAF1,SAF2,SAF3,SAF4)
         IF (IACT.EQ.0) GO TO 999
         IF (IACT.EQ.1.AND.SNEXT.LE.SAFE) GO TO 999
      END IF
*
*     Intersection with z-plane
*     only points outside the z range need to be considered
      IF (ABS(X(3)).GE.P(3)) THEN
         IF (X(3)*X(6).LT.0.) THEN
            S=(ABS(X(3))-P(3))/ABS(X(6))
            XI=X(1)+S*X(4)
            YI=X(2)+S*X(5)
            R2=XI**2+YI**2
            IF (P(1)**2.LE.R2.AND.R2.LE.P(2)**2) THEN
               IF (IFL.EQ.1) GO TO 101
               IF (R2.EQ.0.) GO TO 101
               CPSI=(XI*CFIO+YI*SFIO)/SQRT(R2)
               IF (CPSI.GE.CDFI) GO TO 101
            END IF
         END IF
      END IF
*
*     Intersection with cylinders
*     Intersection point (x,y,z)
*     (x,y,z) is on track :    x=X(1)+t*X(4)
*                              y=X(2)+t*X(5)
*                              z=X(3)+t*X(6)
*     (x,y,z) is on cylinder : x**2 + y**2 = R**2
*
*     (X(4)**2+X(5)**2)*t**2
*     +2.*(X(1)*X(4)+X(2)*X(5))*t
*     +X(1)**2+X(2)**2-R**2=0
*
      T1=X(4)**2+X(5)**2
      T2=(X(1)*X(4)+X(2)*X(5))
      T3=X(1)**2+X(2)**2
*     track parallel to the z axis ?
***** 21-JUN-1990
*     IF (T1.EQ.0.) GO TO 999
      IF (ABS(T1).LT.1.E-32) GO TO 999
      B =T2/T1
*
*     Intersection with outer cylinder
*     only points outside the outer cylinder need to be considered
      IF (R.GE.P(2)) THEN
         C=(T3-P(2)**2)/T1
         D=B**2-C
         IF (D.GE.0.) THEN
            S=-B-SQRT(D)
            IF (S.GE.0.) THEN
               ZI=X(3)+S*X(6)
               IF (ABS(ZI).LE.P(3)) THEN
                  IF (IFL.EQ.1) GO TO 101
                  XI=X(1)+S*X(4)
                  YI=X(2)+S*X(5)
                  CPSI=(XI*CFIO+YI*SFIO)/P(2)
                  IF (CPSI.GE.CDFI) GO TO 101
               END IF
            END IF
         END IF
      END IF
*     Intersection with inner cylinder
      IF (P(1).GT.0.) THEN
         C=(T3-P(1)**2)/T1
         D=B**2-C
         IF (D.GE.0.) THEN
            S=-B+SQRT(D)
            IF (S.GE.0.) THEN
               ZI=X(3)+S*X(6)
               IF (ABS(ZI).LE.P(3)) THEN
                  IF (IFL.EQ.1) GO TO 101
                  XI=X(1)+S*X(4)
                  YI=X(2)+S*X(5)
                  CPSI=(XI*CFIO+YI*SFIO)/P(1)
                  IF (CPSI.GE.CDFI) SNXT=S
               END IF
            END IF
         END IF
      END IF
*
*     Intersection with phi-planes
*     x=r*cos(phi)=X(1)+t*X(4)
*     y=r*sin(phi)=X(2)+t*X(5)
*     z           =X(3)+t*X(6)
*     t=(X(2)*cos(phi)-X(1)*sin(phi))/(X(4)*sin(phi)-X(5)*cos(phi))
      IF (IFL.EQ.2) THEN
*        track not parallel to the phi1 plane ?
         UN=X(4)*S1-X(5)*C1
         IF (UN.NE.0.) THEN
            S=(X(2)*C1-X(1)*S1)/UN
            IF (S.GE.0.) THEN
               ZI=X(3)+S*X(6)
               IF (ABS(ZI).LE.P(3)) THEN
                  XI=X(1)+S*X(4)
                  YI=X(2)+S*X(5)
                  R2=XI**2+YI**2
                  IF (P(1)**2.LE.R2.AND.R2.LE.P(2)**2) THEN
                     IF ((YI*CFIO-XI*SFIO).LE.0.) THEN
                        IF (S.LT.SNXT) SNXT=S
                     END IF
                  END IF
               END IF
            END IF
         END IF
*        track not parallel to the phi2 plane ?
         UN=X(4)*S2-X(5)*C2
         IF (UN.NE.0.) THEN
            S=(X(2)*C2-X(1)*S2)/UN
            IF (S.GE.0.) THEN
               ZI=X(3)+S*X(6)
               IF (ABS(ZI).LE.P(3)) THEN
                  XI=X(1)+S*X(4)
                  YI=X(2)+S*X(5)
                  R2=XI**2+YI**2
                  IF (P(1)**2.LE.R2.AND.R2.LE.P(2)**2) THEN
                     IF ((YI*CFIO-XI*SFIO).GE.0.) THEN
                        IF (S.LT.SNXT) SNXT=S
                     END IF
                  END IF
               END IF
            END IF
         END IF
      END IF
      GO TO 999
*
  101 SNXT=S
  999 END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:20:53 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani
	11	*-- Author :
	12	SUBROUTINE GNOTUB(X,P,IACT,IFL,SNEXT,SNXT,SAFE)
	13	C. ******************************************************************
	14	C. * *
	15	C. * Compute distance to intersection with boundary surface of *
	16	C * volume TUBE or TUBS, from point X(1),X(2),X(3) outside *
	17	C * the volume along track with direction cosines X(4),X(5), *
	18	C * X(6) *
	19	C. * P (input) : volume parameters *
	20	C. * IACT (input) : action flag *
	21	C. * = 0 Compute SAFE only *
	22	C. * = 1 Compute SAFE, compute SNXT only if SAFE.LT.SNEXT *
	23	C. * = 2 Compute both SAFE and SNXT *
	24	C. * = 3 Compute SNXT only *
	25	C. * IFL (input) : 1 for TUBE, 2 for PHI segmented TUBE *
	26	C. * SNEXT (input) : see IACT = 1 *
	27	C. * SNXT (output) : distance to volume boundary along track *
	28	C. * SAFE (output) : not larger than scalar distance to *
	29	C. * volume boundaray *
	30	C. * Called by : GNEXT, GNOPCO, GTNEXT *
	31	C. * *
	32	C. * Authors : Michel Maire and Rolf Nierhaus 21-JUN-1990 *
	33	C. * *
	34	C. ******************************************************************
	35	C. * *
	36	C. * 'TUBE' is a tube. It has 3 parameters, the inside radius, *
	37	C. * the outside radius and the half length in z. *
	38	C. * 'TUBS' is a phi segment of a tube. It has 5 parameters, *
	39	C. * the same 3 as 'TUBE' plus the phi limits. The *
	40	C. * segment starts at the first limit and includes *
	41	C. * increasing phi value up to the second limit or *
	42	C. * that plus 360 degrees. *
	43	C. * *
	44	C. ******************************************************************
	45	#if !defined(CERNLIB_SINGLE)
	46	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	47	#endif
	48	#include "geant321/gconsp.inc"
	49	REAL X(6),P(5),SNEXT,SNXT,SAFE
	50	*
	51	* this part has to be moved outside the routine
	52	IF (IFL.EQ.2) THEN
	53	P4=P(4)*DEGRAD
	54	P5=P(5)*DEGRAD
	55	IF (P5.LT.P4) P5=P5+TWOPI
	56	C1=COS(P4)
	57	S1=SIN(P4)
	58	C2=COS(P5)
	59	S2=SIN(P5)
	60	FIO=0.5*(P5+P4)
	61	CFIO=COS(FIO)
	62	SFIO=SIN(FIO)
	63	DFI=0.5*(P5-P4)
	64	CDFI=COS(DFI)
65	* SDFI=SIN(DFI)
66	END IF
67	*
68	SNXT=1.E10
69	R=SQRT(X(1)2+X(2)2)
70	*
71	* Compute SAFE radius
72	IF (IACT.LT.3) THEN
73	SAF1=P(1)-R
74	SAF2=R-P(2)
75	SAF3=ABS(X(3))-P(3)
76	SAF4=0.
77	IF (IFL.EQ.2.AND.R.GT.0.) THEN
78	CPSI=(X(1)CFIO+X(2)SFIO)/R
79	IF (CPSI.LT.CDFI) THEN
80	IF ((X(2)CFIO-X(1)SFIO).LE.0.) THEN
81	SAF4=ABS(X(1)S1-X(2)C1)
82	ELSE
83	SAF4=ABS(X(1)S2-X(2)C2)
84	END IF
85	END IF
86	ENDIF
87	SAFE=MAX(SAF1,SAF2,SAF3,SAF4)
88	IF (IACT.EQ.0) GO TO 999
89	IF (IACT.EQ.1.AND.SNEXT.LE.SAFE) GO TO 999
90	END IF
91	*
92	* Intersection with z-plane
93	* only points outside the z range need to be considered
94	IF (ABS(X(3)).GE.P(3)) THEN
95	IF (X(3)*X(6).LT.0.) THEN
96	S=(ABS(X(3))-P(3))/ABS(X(6))
97	XI=X(1)+S*X(4)
98	YI=X(2)+S*X(5)
99	R2=XI2+YI2
100	IF (P(1)2.LE.R2.AND.R2.LE.P(2)2) THEN
101	IF (IFL.EQ.1) GO TO 101
102	IF (R2.EQ.0.) GO TO 101
103	CPSI=(XICFIO+YISFIO)/SQRT(R2)
104	IF (CPSI.GE.CDFI) GO TO 101
105	END IF
106	END IF
107	END IF
108	*
109	* Intersection with cylinders
110	* Intersection point (x,y,z)
111	* (x,y,z) is on track : x=X(1)+t*X(4)
112	* y=X(2)+t*X(5)
113	* z=X(3)+t*X(6)
114	* (x,y,z) is on cylinder : x2 + y2 = R**2
115	*
116	* (X(4)2+X(5)2)t*2
117	* +2.(X(1)X(4)+X(2)X(5))t
118	* +X(1)2+X(2)2-R**2=0
119	*
120	T1=X(4)2+X(5)2
121	T2=(X(1)X(4)+X(2)X(5))
122	T3=X(1)2+X(2)2
123	* track parallel to the z axis ?
124	***** 21-JUN-1990
125	* IF (T1.EQ.0.) GO TO 999
126	IF (ABS(T1).LT.1.E-32) GO TO 999
127	B =T2/T1
128	*
129	* Intersection with outer cylinder
130	* only points outside the outer cylinder need to be considered
131	IF (R.GE.P(2)) THEN
132	C=(T3-P(2)**2)/T1
133	D=B**2-C
134	IF (D.GE.0.) THEN
135	S=-B-SQRT(D)
136	IF (S.GE.0.) THEN
137	ZI=X(3)+S*X(6)
138	IF (ABS(ZI).LE.P(3)) THEN
139	IF (IFL.EQ.1) GO TO 101
140	XI=X(1)+S*X(4)
141	YI=X(2)+S*X(5)
142	CPSI=(XICFIO+YISFIO)/P(2)
143	IF (CPSI.GE.CDFI) GO TO 101
144	END IF
145	END IF
146	END IF
147	END IF
148	* Intersection with inner cylinder
149	IF (P(1).GT.0.) THEN
150	C=(T3-P(1)**2)/T1
151	D=B**2-C
152	IF (D.GE.0.) THEN
153	S=-B+SQRT(D)
154	IF (S.GE.0.) THEN
155	ZI=X(3)+S*X(6)
156	IF (ABS(ZI).LE.P(3)) THEN
157	IF (IFL.EQ.1) GO TO 101
158	XI=X(1)+S*X(4)
159	YI=X(2)+S*X(5)
160	CPSI=(XICFIO+YISFIO)/P(1)
161	IF (CPSI.GE.CDFI) SNXT=S
162	END IF
163	END IF
164	END IF
165	END IF
166	*
167	* Intersection with phi-planes
168	* x=rcos(phi)=X(1)+tX(4)
169	* y=rsin(phi)=X(2)+tX(5)
170	* z =X(3)+t*X(6)
171	* t=(X(2)cos(phi)-X(1)sin(phi))/(X(4)sin(phi)-X(5)cos(phi))
172	IF (IFL.EQ.2) THEN
173	* track not parallel to the phi1 plane ?
174	UN=X(4)S1-X(5)C1
175	IF (UN.NE.0.) THEN
176	S=(X(2)C1-X(1)S1)/UN
177	IF (S.GE.0.) THEN
178	ZI=X(3)+S*X(6)
179	IF (ABS(ZI).LE.P(3)) THEN
180	XI=X(1)+S*X(4)
181	YI=X(2)+S*X(5)
182	R2=XI2+YI2
183	IF (P(1)2.LE.R2.AND.R2.LE.P(2)2) THEN
184	IF ((YICFIO-XISFIO).LE.0.) THEN
185	IF (S.LT.SNXT) SNXT=S
186	END IF
187	END IF
188	END IF
189	END IF
190	END IF
191	* track not parallel to the phi2 plane ?
192	UN=X(4)S2-X(5)C2
193	IF (UN.NE.0.) THEN
194	S=(X(2)C2-X(1)S2)/UN
195	IF (S.GE.0.) THEN
196	ZI=X(3)+S*X(6)
197	IF (ABS(ZI).LE.P(3)) THEN
198	XI=X(1)+S*X(4)
199	YI=X(2)+S*X(5)
200	R2=XI2+YI2
201	IF (P(1)2.LE.R2.AND.R2.LE.P(2)2) THEN
202	IF ((YICFIO-XISFIO).GE.0.) THEN
203	IF (S.LT.SNXT) SNXT=S
204	END IF
205	END IF
206	END IF
207	END IF
208	END IF
209	END IF
210	GO TO 999
211	*
212	101 SNXT=S
213	999 END
214