5 * Revision 1.1.1.1 1995/10/24 10:21:44 cernlib
9 #include "geant321/pilot.h"
10 #if !defined(CERNLIB_OLD)
11 *CMZ : 3.21/02 29/03/94 15.22.08 by S.Giani
13 SUBROUTINE GTMEDI (X, NUMED)
15 C. ******************************************************************
17 C. * Finds in which volume/medium the point X is, and updates the *
18 C. * common /GCVOLU/ and the structure JGPAR accordingly. *
20 C. * NUMED returns the tracking medium number, or 0 if point is *
21 C. * outside the experimental setup. *
23 C. * Note : For INWVOL = 2, INFROM set to a positive number is *
24 C. * interpreted by GTMEDI as the number IN of the content *
25 C. * just left by the current track within the mother volume *
26 C. * where the point X is assumed to be. *
28 C. * Note : INFROM is set correctly by this routine but it is *
29 C. * used on entrance only in the case GSNEXT has been called *
30 C. * by the user. In other words the value of INFROM received *
31 C. * on entrance is not considered necessarily valid. This *
32 C. * assumption has been made for safety. A wrong value of *
33 C. * INFROM can cause wrong tracking. *
35 C. * Called by : GTRACK *
36 C. * Authors : S.Banerjee, R.Brun, F.Bruyant, A.McPherson *
39 C. * Modified by S.Giani (1993) to perform the search according *
40 C. * to the new 'virtual divisions' algorithm and to build the *
41 C. * stack of the 'possible overlapping volumes' in the case of *
42 C. * MANY volumes. Any kind of boolean operation is now possible.*
43 C. * Divisions along arbitrary axis are now possible. *
45 C. ******************************************************************
47 #include "geant321/gcflag.inc"
48 #include "geant321/gckine.inc"
49 #include "geant321/gcbank.inc"
50 #include "geant321/gcvolu.inc"
51 #include "geant321/gctrak.inc"
52 #if defined(CERNLIB_USRJMP)
53 #include "geant321/gcjump.inc"
55 #include "geant321/gchvir.inc"
56 #include "geant321/gcvdma.inc"
65 C. ------------------------------------------------------------------
72 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
73 if(ingoto.eq.-1.and.q(jvo+3).lt.0)then
75 elseif(ingoto.eq.0)then
76 call ggperp(x,veccos,ierr)
88 * SECTION I: The /GCVOLU/ table contains the initial guess for a path
89 * in the geometry tree on which X may be found. Look along this
90 * path until X is found inside. This is the starting position.
91 * If this is an ONLY volume with no daughters, we are done;
92 * otherwise reset search record variables, proceed to section II.
94 * *** Check if point is in current volume
100 * *** LSAMVL is a logical variable that indicates whether we are still
101 * *** in the current volume or not. It is used in GTRACK to detect
102 * *** precision problems.
104 C***** Code Expanded From Routine: GTRNSF
106 100 IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
107 XC(1) = X(1) - GTRAN(1,NLEVEL)
108 XC(2) = X(2) - GTRAN(2,NLEVEL)
109 XC(3) = X(3) - GTRAN(3,NLEVEL)
112 XL1 = X(1) - GTRAN(1,NLEVEL)
113 XL2 = X(2) - GTRAN(2,NLEVEL)
114 XL3 = X(3) - GTRAN(3,NLEVEL)
115 XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*
117 XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*
119 XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*
126 C***** End of Code Expanded From Routine: GTRNSF
128 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
130 * Note: At entry the variable INGOTO may contain the index of a volume
131 * contained within the current one at NLEVEL. If so, begin by checking
132 * if X lies inside. This improves the search speed over that of GMEDIA.
135 if(raytra.eq.1..and.imyse.eq.1)then
136 CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
137 CALL GFIND(NAME,'SEEN',ISSEEN)
138 if(isseen.eq.-2.or.isseen.eq.-1)goto 189
140 IF ((INGOTO.LE.0).OR.(INGOTO.GT.NIN)) THEN
144 * *** Entrance in content INGOTO predicted by GTNEXT
148 JVOT = LQ(JVOLUM-IVOT)
149 JPAR = LQ(JGPAR-NLEVEL-1)
152 C***** Code Expanded From Routine: GITRAN
154 C. ------------------------------------------------------------------
156 IF (IROTT .EQ. 0) THEN
157 XT(1) = XC(1) - Q(5+JIN)
158 XT(2) = XC(2) - Q(6+JIN)
159 XT(3) = XC(3) - Q(7+JIN)
162 XL1 = XC(1) - Q(5+JIN)
163 XL2 = XC(2) - Q(6+JIN)
164 XL3 = XC(3) - Q(7+JIN)
166 XT(1) = XL1*Q(JR+1) + XL2*Q(JR+2) + XL3*Q(JR+3)
167 XT(2) = XL1*Q(JR+4) + XL2*Q(JR+5) + XL3*Q(JR+6)
168 XT(3) = XL1*Q(JR+7) + XL2*Q(JR+8) + XL3*Q(JR+9)
171 C***** End of Code Expanded From Routine: GITRAN
173 * * Check if point is in content
175 CALL GINME (XT, Q(JVOT+2), Q(JPAR+1), IYES)
178 * If so, prepare information for volume retrieval, and return
183 NAMES(NL1) = IQ(JVOLUM+IVOT)
184 NUMBER(NL1) = Q(JIN+3)
186 LINMX(NL1) = Q(JVO+3)
187 GONLY(NL1) = Q(JIN+8)
188 IF (LQ(LQ(JVOLUM-IVOT)).EQ.0) THEN
189 NLDEV(NL1) = NLDEV(NLEVEL)
193 CALL GTRMUL (GTRAN(1,NLEVEL), GRMAT(1,NLEVEL), Q(JIN+5),
194 + IROTT, GTRAN(1,NL1), GRMAT(1,NL1))
202 call ggperp(x,veccos,ierr)
213 * End of INGOTO processing
215 189 JPAR = LQ(JGPAR-NLEVEL)
216 CALL GINME (XC, Q(JVO+2), Q(JPAR+1), IYES)
219 * ** Point not in current volume, go up the tree
223 IF (NLEVEL.GT.1) THEN
225 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
228 INFROM=LINDEX(NLEVEL+1)
236 * * Point is outside setup
243 * * Point in current volume but not in INGOTO. We block the
244 * * corresponding volume
246 IF (INGOTO.GT.0) THEN
249 IQ(JIN) = IBSET(IQ(JIN),4)
253 * * Found a volume up the tree which contains our point. We block
254 * * the branch we came up from.
258 IQ(JIN) = IBSET(IQ(JIN),4)
262 * ** Point is in current volume
266 IF (INWVOL.NE.2) INFROM = 0
269 * SECTION II: X is found inside current node at NLEVEL in /GCVOLU/.
270 * Search all contents recursively for any containing X.
271 * Take the first one found, if any, and continue at that
272 * level, incrementing NLEVEL and extending /GCVOLU/ tables.
273 * This is continued until a level is reached where X is not
274 * found in any of the contents, or there are no contents.
275 * Note: Since Section II is re-entered from Section III, a blocking word
276 * is used to mark those contents already checked. Upon exit from Section
277 * II, these blocking words are cleared at NLEVEL, but may remain set in
278 * levels between NLEVEL-1 and NLMIN, if any. They must be cleared at exit.
280 * ** Check contents, if any
282 200 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
284 if(raytra.eq.1..and.imyse.eq.1)then
285 CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
286 CALL GFIND(NAME,'SEEN',ISSEEN)
287 if(isseen.eq.-2.or.isseen.eq.-1)goto 300
290 * * Case with no contents
295 * * Case with contents defined by division
297 ELSEIF (NIN.LT.0) THEN
298 CALL GMEDIV (JVO, IN, XC, 1)
300 if(raytra.eq.1..and.neufla.eq.1)then
302 call ggperp(x,veccos,ierr)
316 * * Case with contents positioned
320 clmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+3)
321 chmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+4)
322 ndivto=q(jvirt+4*(LVOLUM(NLEVEL)-1)+2)
323 iaxis =q(jvirt+4*(LVOLUM(NLEVEL)-1)+1)
325 ivdiv=((xc(iaxis)-clmoth)*ndivto/(chmoth-clmoth))+1
328 elseif(ivdiv.gt.ndivto)then
332 call gfcoor(xc,iaxis,cx)
334 if((cx-clmoth).lt.-1.)then
336 elseif((cx-chmoth).gt.1.)then
341 elseif(cx.lt.clmoth)then
345 ivdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
348 elseif(ivdiv.gt.ndivto)then
352 jvdiv=lq(jvirt-LVOLUM(NLEVEL))
353 iofset=iq(jvdiv+ivdiv)
354 ncont=iq(jvdiv+iofset+1)
356 if(ncont.eq.0)goto 260
358 JCONT = LQ(JVO-NIN-1)+1
362 * For each selected content in turn, check if point is inside
372 * If the value IQ(JCONT+ICONT)=0 then we are back in the mother.
373 * So jump to 260, the search is finished. Clean-up should be done
374 * only up to ICONT-1, so we set:
380 IF (.NOT.BTEST(IQ(JIN),4)) THEN
381 CALL GMEPOS (JVO, IN, XC, 1)
384 IF (GONLY(NLEVEL).NE.0.) THEN
395 IQ(JIN) = IBSET(IQ(JIN),4)
401 260 IF(NCONT.EQ.NIN) THEN
404 IQ(JIN) = IBCLR(IQ(JIN),4)
414 IQ(JIN) = IBCLR(IQ(JIN),4)
418 IQ(JIN) = IBCLR(IQ(JIN),4)
423 IQ(JIN) = IBCLR(IQ(JIN),4)
430 * SECTION III: X is found at current node (NLEVEL in /GCVOLU) but not in
431 * any of its contents, if any. If this is a MANY volume,
432 * save it as a candidate best-choice, and continue the search
433 * by backing up the tree one node and proceed to Section II.
434 * If this is an ONLY volume, proceed to Section IV.
436 * *** Point is in current volume/medium, and not in any content
438 300 IF (GONLY(NLEVEL).EQ.0.) THEN
440 * ** Lowest level is 'NOT ONLY'
442 IF (NLEVEL.GT.NLMANY) THEN
449 manyle(nvmany)=nlevel
451 manyna(nvmany,i)=names(i)
452 manynu(nvmany,i)=number(i)
456 * * Go up the tree up to a volume with positioned contents
459 310 INFROM = LINDEX(NLEVEL)
461 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
463 IF (NIN.LT.0) GO TO 310
465 C***** Code Expanded From Routine: GTRNSF
467 IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
468 XC(1) = X(1) - GTRAN(1,NLEVEL)
469 XC(2) = X(2) - GTRAN(2,NLEVEL)
470 XC(3) = X(3) - GTRAN(3,NLEVEL)
473 XL1 = X(1) - GTRAN(1,NLEVEL)
474 XL2 = X(2) - GTRAN(2,NLEVEL)
475 XL3 = X(3) - GTRAN(3,NLEVEL)
476 XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) +
477 + XL3* GRMAT(3,NLEVEL)
478 XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) +
479 + XL3* GRMAT(6,NLEVEL)
480 XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) +
481 + XL3* GRMAT(9,NLEVEL)
484 C***** End of Code Expanded From Routine: GTRNSF
488 IQ(JIN) = IBSET(IQ(JIN),4)
489 NLMIN = MIN(NLEVEL,NLMIN)
493 * SECTION IV: This is the end of the search. The current node (NLEVEL
494 * in /GCVOLU/) is the lowest ONLY volume in which X is found.
495 * If X was also found in any of its contents, they are MANY
496 * volumes: the best-choice is the one among them at the greatest
497 * level in the tree, and it is stored. Otherwise the current
498 * volume is the solution. Before exit, all of the blocking
499 * words leftover in the tree must be reset to zero.
500 * Note: A valid structure is assumed, in which no ONLY volumes overlap.
501 * If this rule is violated, or if a daughter is not entirely contained
502 * within the mother volume, the results are unpredictable.
504 DO 419 NL=NLMIN,NLEVEL-1
505 JVO = LQ(JVOLUM-LVOLUM(NL))
509 IQ(JIN) = IBCLR(IQ(JIN),4)
517 IF (NLMANY.GT.0) CALL GFCVOL
518 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
519 IF(JVIN.NE.0) IQ(JVIN) = IBCLR(IQ(JVIN),4)
522 999 IF(JGSTAT.NE.0) CALL GFSTAT(4)