This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / glisur.F
1 *PATCH,GEXAM1.
2 *CMZ :  3.21/02 29/03/94  15.41.35  by  S.Giani
3 *DECK,GLISUR
4 *CMZ :          12/09/95  11.22.49  by  S.Ravndal
5 *CMZ :  3.21/02 03/07/94  19.08.44  by  S.Giani
6 *-- Author :
7       SUBROUTINE GLISUR(X0,X1,MEDI0,MEDI1,U,PDOTU,IERR)
8 C.
9 C.    ******************************************************************
10 C.    *                                                                *
11 C.    *  This routine simulates the surface profile of a boundary      *
12 C.    *  between two media as seen by an approaching particle with     *
13 C.    *  coordinates and direction given by X0. The surface is         *
14 C.    *  identified by the arguments MEDI0 and MEDI1 which are the     *
15 C.    *  media indices of the region in which the track is presently   *
16 C.    *  and the one which it approaches, respectively. The input      *
17 C.    *  vector X1 contains the coordinates of a point on the other    *
18 C.    *  side of the boundary from X0, and lying within medium MEDI1.  *
19 C.    *  The result is a unit vector U normal to the surface of        *
20 C.    *  reflection at X0 and pointing into the medium from which the  *
21 C.    *  track is approaching. The quality of the surface finish is    *
22 C.    *  described by the parameter POLISH, which varies from 0 to 1.  *
23 C.    *  POLISH=0 means maximum roughness, with effective plane of     *
24 C.    *  reflection distributed as cos(alpha), where alpha is the      *
25 C.    *  angle between the unit normal to the effective plane of       *
26 C.    *  reflection and the normal to the nominal medium boundary at   *
27 C.    *  X0. POLISH=1 means perfect smoothness. In between the surface *
28 C.    *  is modeled by a bell-shaped distribution in alpha, with       *
29 C.    *  limits defined by                                             *
30 C.    *              sin(alpha) = +/- (1-POLISH)                       *
31 C.    *  At the interface between two media, the surface POLISH        *
32 C.    *  parameter is taken from a user routine GUPLSH(MEDI0,MEDI1).   *
33 C.    *  The indices MEDI0 and MEDI1 refer to the media declared by    *
34 C.    *  the user. If GGPERP returns an error, there was a precision   *
35 C.    *  problem with the tracking, and point X0 is on the same side   *
36 C.    *  of the surface as X1. In this case (or any other error        *
37 C.    *  condition from GGPERP) a return is made with IERR ^=0. If     *
38 C.    *  IERR=0 on return then U contains a valid unit vector.         *
39 C.    *                                                                *
40 C.    *    ==>Called by : <USER>, UGINIT                               *
41 C.    *       Author    F.Carminati, R.Jones ***********               *
42 C.    *                                                                *
43 C.    ******************************************************************
44 C.
45 #include "geant321/gcvolu.inc"
46 #include "geant321/gcvol1.inc"
47 C
48       REAL X0(6),X1(3),U(3)
49       REAL D(3),RNDM(3)
50
51 *
52 * *** Decide which volume defines the surface
53       IF (NLEVEL.GE.NLEVL1) THEN
54          LVOL = 1
55       ELSE
56          LVOL = -1
57          DO 10 I=2,NLEVEL
58             IF ((NAMES1(I).NE.NAMES(I)).OR. (NUMBR1(I).NE.NUMBER(I)))
59      +      LVOL = 1
60    10    CONTINUE
61       END IF
62 *
63 * *** Find the nominal unit normal to the surface
64       IF (LVOL.EQ.1) THEN
65 C
66 C Instead of X0, now X1 is used by GGPERP, in order to calculate
67 C the right perpendicular even near to a corner,
68 C see the modifactions in GTCKOV concerning VBOU
69 C
70          CALL GGPERP(X1,U,IERR)
71       ELSE
72          LVOL = NLEVEL
73          NLEVEL = NLEVL1
74          CALL GGPERP(X1,U,IERR)
75          U(1) = -U(1)
76          U(2) = -U(2)
77          U(3) = -U(3)
78          NLEVEL = LVOL
79       END IF
80       IF (IERR.NE.0) GO TO 999
81       PDOTU = X0(4)*U(1) + X0(5)*U(2) + X0(6)*U(3)
82       POLISH=GUPLSH(MEDI0,MEDI1)
83       IF (POLISH.LT.1.) THEN
84 *
85 * *** Generate distortion vector D within a uniform sphere
86          DIAM = 2.*(1.-POLISH)
87          DIA2 = DIAM**2
88          IROUND = 0
89    20    IF(IROUND.GT.5) GO TO 999
90          IROUND = IROUND+1
91          CALL GRNDM(RNDM,3)
92          D(1) = RNDM(1)-0.5
93          D(2) = RNDM(2)-0.5
94          D(3) = RNDM(3)-0.5
95          D2 = D(1)**2+D(2)**2+D(3)**2
96          IF (D2.GT.0.25) GO TO 20
97          D(1) = D(1)*DIAM
98          D(2) = D(2)*DIAM
99          D(3) = D(3)*DIAM
100          D2 = D2*DIA2
101 *
102 * *** Insure that V=U+D will cause reflection away from surface
103          PDOTD = X0(4)*D(1) + X0(5)*D(2) + X0(6)*D(3)
104          PDOTV = PDOTU+PDOTD
105          UDOTD = U(1)*D(1) + U(2)*D(2) + U(3)*D(3)
106          V2 = 1+D2+2*UDOTD
107          IF (PDOTD*V2+PDOTV*(1.-D2).GT.0.) GO TO 20
108 *
109 * *** Normalize V and call it U
110          VABS = 1./SQRT(V2)
111          U(1) = (U(1)+D(1))*VABS
112          U(2) = (U(2)+D(2))*VABS
113          U(3) = (U(3)+D(3))*VABS
114          PDOTU = PDOTV*VABS
115       END IF
116   999 END