]>
Commit | Line | Data |
---|---|---|
1 | * | |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.2 2001/02/24 10:46:32 fca | |
6 | * Moving to double precision and introducing an additional security check | |
7 | * | |
8 | * Revision 1.1.1.1 1999/05/18 15:55:17 fca | |
9 | * AliRoot sources | |
10 | * | |
11 | * Revision 1.1.1.1 1995/10/24 10:20:54 cernlib | |
12 | * Geant | |
13 | * | |
14 | * | |
15 | #include "geant321/pilot.h" | |
16 | *CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani | |
17 | *-- Author : | |
18 | SUBROUTINE GNPGON (X, PAR, IACT, SNEXT, SNXT, SAFE) | |
19 | C. | |
20 | C. ****************************************************************** | |
21 | C. * * | |
22 | C. * COMPUTE DISTANCE UP TO INTERSECTION WITH 'PGON' VOLUME, * | |
23 | C. * FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6) * | |
24 | C. * * | |
25 | C. * PAR (input) : volume parameters * | |
26 | C. * IACT (input) : action flag * | |
27 | C. * = 0 Compute SAFE only * | |
28 | C. * = 1 Compute SAFE, and SNXT only if SNEXT .GT.new SAFE * | |
29 | C. * = 2 Compute both SAFE and SNXT * | |
30 | C. * = 3 Compute SNXT only * | |
31 | C. * SNEXT (input) : see IACT = 1 * | |
32 | C. * SNXT (output) : distance to volume boundary * | |
33 | C. * SAFE (output) : shortest distance to any boundary * | |
34 | C. * * | |
35 | C. * ==>Called by : GNEXT, GTNEXT * | |
36 | C. * Author A.McPherson, P.Weidhaas ********* * | |
37 | C. * * | |
38 | C. ****************************************************************** | |
39 | C. | |
40 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
41 | #include "geant321/gconsp.inc" | |
42 | ||
43 | REAL X(6), PAR(100), SNEXT, SNXT, SAFE | |
44 | DIMENSION TPAR(100) | |
45 | ||
46 | C----------------------------------------------------------------- | |
47 | ||
48 | SNXT = BIG | |
49 | SAFSEG = BIG | |
50 | R2 = X(1)*X(1) + X(2)*X(2) | |
51 | R = SQRT (R2) | |
52 | ||
53 | PDIV = PAR(2) / PAR(3) | |
54 | DELPHI = PDIV * DEGRAD | |
55 | DPHI2 = 0.5 * DELPHI | |
56 | CSDPH2 = COS (DPHI2) | |
57 | ||
58 | TPAR(1) = PAR(1) | |
59 | TPAR(2) = PAR(2) | |
60 | TPAR(3) = PAR(3) | |
61 | TPAR(4) = PAR(4) | |
62 | NZ = PAR(4) | |
63 | ||
64 | DO 5 I=1, NZ | |
65 | I3 = 3*(I-1) | |
66 | TPAR(5+I3) = PAR(5+I3) | |
67 | TPAR(6+I3) = PAR(6+I3) / CSDPH2 | |
68 | TPAR(7+I3) = PAR(7+I3) | |
69 | 5 CONTINUE | |
70 | ||
71 | C************************************************************************** | |
72 | C | |
73 | C...... Here we start the logic from "GNPCON" (which for reasons of | |
74 | C...... efficiency and clarity has been implemented inline). | |
75 | C | |
76 | C************************************************************************** | |
77 | ||
78 | ZMIN = TPAR(5) | |
79 | ZMAX = TPAR(3*NZ+2) | |
80 | SAFZ1 = X(3) - ZMIN | |
81 | SAFZ2 = ZMAX - X(3) | |
82 | SAFEZ = MIN (SAFZ1, SAFZ2) | |
83 | ||
84 | C | |
85 | C...... First determine in which segment the particle is located. | |
86 | C | |
87 | DO 10 JPH=8, 3*NZ-1, 3 | |
88 | IF (X(3) .LT. TPAR(JPH)) THEN | |
89 | IPH=JPH | |
90 | GO TO 20 | |
91 | ENDIF | |
92 | 10 CONTINUE | |
93 | IPH = 3*NZ+2 | |
94 | ||
95 | 20 CONTINUE | |
96 | C | |
97 | C...... The particle is in the segment bounded by z-planes at | |
98 | C...... Z1=PAR(IPL) and Z2=PAR(IPH), i.e., Z1 < X(3) < Z2. | |
99 | C | |
100 | C...... Set parameters for this segment and translate z-coordinate | |
101 | C...... of point relative to center of this segment. this is done in | |
102 | C...... preparation of invoking the algorithms used in "GNTUBE" and | |
103 | C...... "GNCONE" (which for reasons of efficiency and clarity are | |
104 | C...... implemented inline). | |
105 | C | |
106 | IPL = IPH - 3 | |
107 | DZ = 0.5 * (TPAR(IPH) - TPAR(IPL)) | |
108 | PT2 = TPAR(IPL+1) | |
109 | PT3 = TPAR(IPL+2) | |
110 | PT4 = TPAR(IPH+1) | |
111 | PT5 = TPAR(IPH+2) | |
112 | PT6 = TPAR(1) | |
113 | PT7 = TPAR(1) + TPAR(2) | |
114 | IF (PT7 .GT. 360.0) PT7 = PT7 - 360.0 | |
115 | ||
116 | XT3 = X(3) - 0.5 * (TPAR(IPL) + TPAR(IPH)) | |
117 | ||
118 | SAFZ2 = DZ + XT3 | |
119 | ZLENI = 0.5 / DZ | |
120 | ||
121 | PHI1 = PT6 * DEGRAD | |
122 | PHI2 = PT7 * DEGRAD | |
123 | IF (PHI2.LE.PHI1) PHI2=PHI2+TWOPI | |
124 | ||
125 | IF (IACT .LT. 3) THEN | |
126 | ||
127 | C ------------------------------------------------- | |
128 | C | Compute safety-distance 'SAFE' (P.Weidhaas) | | |
129 | C ------------------------------------------------- | |
130 | ||
131 | IFLAG = 2 | |
132 | IF (TPAR(2) .EQ. 360.0) IFLAG = 1 | |
133 | SAFZ1 = DZ - ABS(XT3) | |
134 | SAFEZ = MIN (SAFEZ,SAFZ1) | |
135 | C | |
136 | C...... Next determine whether the segment is a tube or a cone. | |
137 | C | |
138 | IF (PT2 .NE. PT4) GO TO 50 | |
139 | IF (PT3 .NE. PT5) GO TO 50 | |
140 | ||
141 | C********************************************************* | |
142 | C | |
143 | C...... The segment is a tube; invoke the algorithm | |
144 | C...... from routine "GNTUBE" inline to get "SAFER". | |
145 | C | |
146 | C********************************************************* | |
147 | ||
148 | SAFR1 = BIG | |
149 | IF(PT2.GT.0.) SAFR1 = R - PT2 | |
150 | SAFR2 = PT3 - R | |
151 | SAFER = MIN (SAFR1, SAFR2) | |
152 | ||
153 | IF (IFLAG .EQ. 2) GO TO 70 | |
154 | ||
155 | GO TO 100 | |
156 | ||
157 | ||
158 | 50 CONTINUE | |
159 | ||
160 | C********************************************************* | |
161 | C | |
162 | C...... The segment is a cone; invoke the algorithm | |
163 | C...... from routine "GNCONE" inline to get "SAFER". | |
164 | C | |
165 | C********************************************************* | |
166 | ||
167 | ||
168 | C...... Compute radial distance to inner wall. | |
169 | ||
170 | IF(PT2+PT4.GT.0.) THEN | |
171 | FACT = (PT4 - PT2) * ZLENI | |
172 | RAD1 = PT2 + FACT * SAFZ2 | |
173 | SAFR1 = (R - RAD1) / SQRT(1.0 + FACT*FACT) | |
174 | ELSE | |
175 | SAFR1 = BIG | |
176 | ENDIF | |
177 | ||
178 | C...... Compute radial distance to outer wall. | |
179 | ||
180 | FACT = (PT5 - PT3) * ZLENI | |
181 | RAD2 = PT3 + FACT * SAFZ2 | |
182 | SAFR2 = (RAD2 - R) / SQRT(1.0 + FACT*FACT) | |
183 | ||
184 | SAFER = MIN (SAFR1, SAFR2) | |
185 | ||
186 | IF (IFLAG .EQ. 1) GO TO 100 | |
187 | ||
188 | 70 CONTINUE | |
189 | ||
190 | C******************************************************************** | |
191 | C...... Here we handle the case of a phi-segment of a tube or cone. | |
192 | C...... in addition to the radial distances (SAFR1, SAFR2) and the | |
193 | C...... axial distances (SAFZ1, SAFZ2) we compute here the distance | |
194 | C...... to the phi-segment boundary that is closest to the point: | |
195 | C | |
196 | C...... For each phi-boundary we find the distance from the given | |
197 | C...... point to the outer (at R2) point of the segment boundary | |
198 | C...... (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define | |
199 | C...... "SAFSEG" to be the distance to segment PHI1, else we set | |
200 | C...... "SAFSEG" to be the distance to segment PHI2. | |
201 | C********************************************************************* | |
202 | ||
203 | ||
204 | COSPH1 = COS (PHI1) | |
205 | SINPH1 = SIN (PHI1) | |
206 | COSPH2 = COS (PHI2) | |
207 | SINPH2 = SIN (PHI2) | |
208 | ||
209 | C...... Get coordinates of outer endpoints (at R2) of both phi-segments. | |
210 | ||
211 | XS1 = R * COSPH1 | |
212 | YS1 = R * SINPH1 | |
213 | XS2 = R * COSPH2 | |
214 | YS2 = R * SINPH2 | |
215 | ||
216 | C...... Get distances (squared) from given point to each endpoint. | |
217 | ||
218 | DISTS1 = (X(1) - XS1)**2 + (X(2) - YS1)**2 | |
219 | DISTS2 = (X(1) - XS2)**2 + (X(2) - YS2)**2 | |
220 | ||
221 | C...... Get distance to that phi-segment whose endpoint | |
222 | C...... is closest to the given point. | |
223 | ||
224 | IF (DISTS1 .LE. DISTS2) THEN | |
225 | SAFSEG = ABS(SINPH1 * X(1) - COSPH1 * X(2)) | |
226 | ELSE | |
227 | SAFSEG = ABS(SINPH2 * X(1) - COSPH2 * X(2)) | |
228 | ENDIF | |
229 | ||
230 | 100 CONTINUE | |
231 | ||
232 | ||
233 | IF (SAFER .LE. 0.0) THEN | |
234 | ||
235 | C--------------------------------------------------------------------------- | |
236 | C | |
237 | C...... Here we handle the case in which SAFER < 0, i.e., the point is | |
238 | C...... inside the polygon but outside the inscribed polycone. We must | |
239 | C...... do an accurate calculation of "SAFER". | |
240 | C | |
241 | C--------------------------------------------------------------------------- | |
242 | ||
243 | FACT = SAFZ2 * ZLENI | |
244 | RAD1 = PT2 + FACT * (PT4 - PT2) | |
245 | RAD2 = PT3 + FACT * (PT5 - PT3) | |
246 | RR1 = RAD1 * RAD1 | |
247 | RR2 = RAD2 * RAD2 | |
248 | ||
249 | IF(X(1).EQ.0.)THEN | |
250 | PHI = ATAN2(X(2), 1.E-8) | |
251 | ELSE | |
252 | PHI = ATAN2(X(2), X(1)) | |
253 | ENDIF | |
254 | IF (PHI .LT. PHI1) PHI = PHI + TWOPI | |
255 | ||
256 | DIST=0. | |
257 | IF (PHI .GE. PHI1 .AND. PHI .LE. PHI2) THEN | |
258 | PHIREL = PHI - PHI1 | |
259 | NSECTR = INT(PHIREL / DELPHI) + 1 | |
260 | PHICTR = PHI1 + (2.0*NSECTR - 1.0) * DPHI2 | |
261 | COSPHC = COS (PHICTR) | |
262 | SINPHC = SIN (PHICTR) | |
263 | ||
264 | ||
265 | IF (R2 .GE. RR2) THEN | |
266 | SR2 = RAD2 | |
267 | DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SR2) | |
268 | FACTC = (PT5 - PT3) * ZLENI | |
269 | ELSEIF (R2 .LE. RR1) THEN | |
270 | RIN = RAD1 * CSDPH2 | |
271 | SRIN = RIN | |
272 | DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SRIN) | |
273 | FACTC = (PT4 - PT2) * ZLENI | |
274 | ELSE | |
275 | PRINT *, 'GNPGON: FACTC not calculated!!!' | |
276 | FACTC=0 | |
277 | ENDIF | |
278 | ||
279 | ENDIF | |
280 | ||
281 | SAFER = DIST / SQRT(1.0 + FACTC*FACTC) | |
282 | ||
283 | ENDIF | |
284 | ||
285 | SAFE = MIN (SAFEZ, SAFER, SAFSEG) | |
286 | ||
287 | ||
288 | IF (IACT .EQ. 0) GO TO 999 | |
289 | IF (IACT .EQ. 1) THEN | |
290 | IF (SNEXT .LT. SAFE) GO TO 999 | |
291 | ENDIF | |
292 | ENDIF | |
293 | ||
294 | ||
295 | C ------------------------------------------------ | |
296 | C | Compute vector-distance 'SNXT' (Nierhaus) | | |
297 | C ------------------------------------------------ | |
298 | ||
299 | ||
300 | CALL GNPGO1(X,PAR,SNXT) | |
301 | C | |
302 | 999 CONTINUE | |
303 | END | |
304 |